Merge pull request #99 from emilelavaut/develop

space-point  and waveform add

Author: aolivier23

duneana/CalibAna/CalibAnaTree.h

@@ -206,192 +206,6 @@ class dune::CalibAnaTree : public art::EDAnalyzer {
 
   void DoTailFit();
 
-  //single hit fct
-  /*
-  bool AllSame(std::vector<int> v);
-
-  std::vector<std::string> GetGeneratorTag(art::Event const& e,
-    art::InputTag fG4producer, 
-    art::ServiceHandle<cheat::BackTrackerService> bt_serv);
-
-  bool Inside( int k , std::list<int> liste);
-
-  float GetDist( float x0 , float y0 , float z0 , 
-		 float x1 , float y1 , float z1 );
-
-  int NearOrFar( bool IsPDVD , bool IsPDHD , const recob::Hit & hit);
-
-  void GetTimeIsolation(art::Event const & ev, 
-    art::InputTag HitLabel, 
-    float const PeakTimeWdInt, 
-    float const PeakTimeWdExt, 
-    std::list<int> & index_list_single, 
-    std::list<int> & index_listIsolated);
-
-  void GetListOfTimeCoincidenceHit( bool IsPDVD , bool IsPDHD ,
-    art::Event const & ev,  
-    art::InputTag HitLabel, 
-    const float CoincidenceWd1_l, 
-    const float CoincidenceWd1_r,
-    const float CoincidenceWd2_l, 
-    const float CoincidenceWd2_r, 
-    const recob::Hit & HitCol,
-    std::list<geo::WireID> & WireInd1,
-    std::list<geo::WireID> & WireInd2,
-    std::list<int>   & ChannelInd1,
-    std::list<int>   & ChannelInd2,
-    std::list<float> & EInd1,
-    std::list<float> & EInd2,
-    std::list<float> & PTInd1,
-    std::list<float> & PTInd2,
-    std::list<float> & PAInd1,
-    std::list<float> & PAInd2);
-
-  void GetListOfCrossingChannel(bool IsPDVD , bool IsPDHD, 
-    float Ymin , float Ymax , float Zmin , float Zmax,
-    geo::WireID & WireCol, 
-    std::list<geo::WireID> & WireInd1,
-    std::list<geo::WireID> & WireInd2,
-    std::list<int>  & ChInd1, 
-    std::list<float> & EInd1,
-    std::list<float> & YInd1,
-    std::list<float> & ZInd1,
-    std::list<int>  & ChIntersectInd1,
-    std::list<float> & EIntersectInd1,
-    std::list<int>  & ChInd2, 
-    std::list<float> & EInd2, 
-    std::list<float> & YInd2, 
-    std::list<float> & ZInd2,
-    std::list<int>  & ChIntersectInd2,
-    std::list<float> & EIntersectInd2);
-
-  void GetListOf3ViewsPoint( float pitch , float alpha ,
-    std::list<int> & ChIntersectInd1,
-    std::list<float> YInd1,
-    std::list<float> ZInd1,
-    std::list<float> EIntersectInd1,
-    std::list<int> & ChIntersectInd2,
-    std::list<float> YInd2,
-    std::list<float> ZInd2,
-    std::list<float> EIntersectInd2,
-    std::list<float> & listYSP,
-    std::list<float> & listZSP,
-    std::list<float> & listEind1SP, 
-    std::list<float> & listEind2SP,
-    std::list<int> & listCh1SP,
-    std::list<int> & listCh2SP);
-
-  std::vector<int> GetXYZIsolatedPoint( 
-    std::vector<float> vYPoint,
-    std::vector<float> vZPoint,
-    std::vector<float> vPeakTimeCol,
-    std::vector<int>   vNOF,
-    float fElectronVelocity,
-    float fTickToMus,
-    float radiusInt,
-    float radiusExt);
-
-  bool IntersectOutsideOfTPC(
-    float Ymin , float Ymax , float Zmin , float Zmax, 
-    double ChInd_start_y,
-    double ChInd_start_z,
-    double ChInd_end_y,
-    double ChInd_end_z,
-    double ChCol_start_y,
-    double ChCol_start_z,
-    double ChCol_end_y,
-    double ChCol_end_z,
-    double& y , double& z );
-
-  point gen_yz(int size,
-    std::vector<int> vIndex,
-    std::vector<float> vY,
-    std::vector<float> vZ,
-    std::vector<int> vNOF);
-
-  float dist2(point a, point b);
-
-  float randf(float m);
-
-  int nearest(point pt, point cent, int n_cluster, float *d2);
-
-  int reallocate(point pt, std::vector<std::vector<float>> ClusterPosition , float threshold);
-
-  float GetDist2D(float y0,float z0,float y1,float z1);
-
-  float mean(float y,float z);
-
-  void kpp(point pts, int len, point cent, int n_cent);
-
-  std::vector<std::vector<float> > lloyd(point pts, int len, int n_cluster);
-
-  std::vector<std::vector<float> > GetData(int len,point data);
-
-  std::vector<int> CheckCompletude(
-    std::vector<std::vector<float> > &data,
-    std::vector<std::vector<float> > &cluster,
-    float RMS,
-    float mult);
-
-  std::vector<int> CheckClusters(
-    std::vector<std::vector<float> > &data,
-    std::vector<std::vector<float> > &cluster,
-    float RMS, 
-    float mult, 
-    float tmp);
-
-  std::vector<ClusterInfo*> GetCluster( 
-    bool fAsConverged,
-    float CRP_T0,
-    int n_point,
-    int n_cluster,
-    point p,
-    std::vector<float> vEInd1PointByEvent,
-    std::vector<float> vEInd2PointByEvent,           
-    std::vector<int> vChInd1PointByEvent,
-    std::vector<int> vChInd2PointByEvent,
-    std::vector<float> vEnergyColByEvent,
-    std::vector<float> vPeakTimeColByEvent,
-    std::vector<int> vChannelColByEvent,
-    bool truth,
-    std::vector<int> vMCPDGByEvent,
-    std::vector<int> vMCMOMpdgByEvent,
-    std::vector<float> vMCWeightByEvent,  
-    std::vector<std::string> vGeneratorTagByEvent,
-    std::vector<float> vMCXByEvent,
-    std::vector<float> vMCYByEvent,
-    std::vector<float> vMCZByEvent,
-    std::vector<int> vNoFByEvent,
-    std::vector<float> vMCEByEvent,
-   std::vector<int> vMCNeByEvent);
-
-  std::vector<ClusterInfo*> SingleHitAnalysis(
-    art::Event const& e,
-    art::InputTag fRDTLabel,
-    art::InputTag fG4producer,
-    art::InputTag fHITproducer, 
-    bool  bIsPDVD,
-    bool  bIsPDHD, 
-    float fCoincidenceWd1_left, 
-    float fCoincidenceWd1_right, 
-    float fCoincidenceWd2_left, 
-    float fCoincidenceWd2_right, 
-    bool  bIs3ViewsCoincidence, 
-    float fPitch, 
-    float fPitchMultiplier,
-    bool  fVerbose,
-    float fMinSizeCluster, 
-    float fMaxSizeCluster,
-    float fNumberInitClusters, 
-    float fRadiusInt,
-    float fRadiusExt, 
-    float fgeoYmin, 
-    float fgeoYmax, 
-    float fgeoZmin,
-    float fgeoZmax,
-    float fElectronVelocity,
-    float fTickTimeInMus);
-*/
   // declare truth utils, ported from CAFana in SBNCode
   std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>>
   PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels,
@@ -421,7 +235,7 @@ class dune::CalibAnaTree : public art::EDAnalyzer {
 
   float fRadiusInt;
   float fRadiusExt;
-  //float fElectronVelocity;
+  //float fElectronVelocity; //You have a local variable of the same name, and it's only used in analyze().  So commenting this out to make Clang happy.
   float fCoincidenceWd1_left;
   float fCoincidenceWd1_right;
   float fCoincidenceWd2_left;

duneana/CalibAna/CalibAnaTreeObj.h

@@ -74,6 +74,16 @@ namespace dune {
       end(-1) {}
   };
 
+  struct ClusterHitInfo {
+    HitInfo h; //!< Hit information by itself
+    std::vector<Vector3D> vSpacePoint; //!< Associated LE-SpacePoint created with SingleHit method
+                                       //!< Several LE-SpacePoint can be associated to one hit 
+				       //!< x coordinate is not reliable x = PeakTime*V_drift
+    ClusterHitInfo():
+      vSpacePoint( {Vector3D()} ) {}
+
+  };
+
   struct TrackHitInfo {
     HitInfo h; //!< Hit information by itself
     float pitch; //!< Pitch of track across wire the hit is on [cm]
@@ -330,6 +340,17 @@ namespace dune {
      int NumberOfInduction1; //!< number of induction1 wires/channels in cluster
      int NumberOfInduction2; //!< number of induction2 wires/channels in cluster
 
+
+     std::vector<ClusterHitInfo> vHitsInfo; //!< vector of HitInfo for all collection hits in the LE-Cluster
+                                            //!< verif vHisInfo.size() = NumberOfCollection
+
+     //Wire information/Wave-Forms
+     //<! time window for wave-form = peaktime +/- Rext/V_drift 
+     //<! wire windows spreading of cluster +/- fHitRawDigitsWireCollectWidth (same as for track)
+     std::vector<WireInfo> wires0; //!< induction 1
+     std::vector<WireInfo> wires1; //!< induction 2
+     std::vector<WireInfo> wires2; //!< collection
+
      //truth info
      bool               truth;      //!< yes if truth info no if not
      std::vector<float> vMC_energy; //!< vector of energy for collection hits in the cluster

duneana/CalibAna/CalibAnaTree_module.cc

@@ -69,8 +69,8 @@ dune::CalibAnaTree::CalibAnaTree(fhicl::ParameterSet const& p)
   fLowEnergyClusterAnalysis     = p.get<bool>("LowEnergyClusterAnalysis", true);
   fRDTLabel                     = p.get<std::string>("RDTLabel","tpcrawdecoder:daq");
 
-  fRadiusInt                    = p.get<float>("RadiusInt");
-  fRadiusExt                    = p.get<float>("RadiusExt");
+  fRadiusInt                    = p.get<float>("RadiusInt",4);
+  fRadiusExt                    = p.get<float>("RadiusExt",10);
 
   fCoincidenceWd1_left          = p.get<float>("CoincidenceWindow1_left"  ,3); //in mus,
   fCoincidenceWd1_right         = p.get<float>("CoincidenceWindow1_right",3); //in mus,
@@ -80,7 +80,7 @@ dune::CalibAnaTree::CalibAnaTree(fhicl::ParameterSet const& p)
   fPitch                        = p.get<float>("Pitch",0.5); //cm
   fPitchMultiplier              = p.get<float>("PitchMultiplier",1.2); // 20% error    
 
-  bIs3ViewsCoincidence          = p.get<bool>("Is3ViewsCoincidence");
+  bIs3ViewsCoincidence          = p.get<bool>("Is3ViewsCoincidence",false);
   bIsPDVD                       = p.get<bool>("IsPDVD",false);
   bIsPDHD                       = p.get<bool>("IsPDHD",false);
   bIsFDVD                       = p.get<bool>("IsFDVD",false);
@@ -272,6 +272,50 @@ void dune::CalibAnaTree::analyze(art::Event const& e)
       std::cout<<"Received "<<simchannels.size()<<" SimChannels for this event."<<std::endl;
   }
 
+  // Collect raw digits for saving hits
+  std::vector<art::Ptr<raw::RawDigit>> rawdigitlist;
+  for (const art::InputTag &t: fRawDigitproducers) {
+    try {
+      art::ValidHandle<std::vector<raw::RawDigit>> thisdigits = e.getValidHandle<std::vector<raw::RawDigit>>(t);
+      art::fill_ptr_vector(rawdigitlist, thisdigits);
+    }
+    catch(...) {
+      if (!fSilenceMissingDataProducts) throw;
+      else { if (fVerbose) {std::cout<<" no raw digits " <<std::endl;}} // Allow Raw Digits to not be present
+    }
+  }
+
+  // The raw digit list is not sorted, so make it into a map on the WireID
+  std::map<geo::WireID, art::Ptr<raw::RawDigit>> rawdigits;
+  for (const art::Ptr<raw::RawDigit> &d: rawdigitlist) {
+
+    std::vector<geo::WireID> wids;
+    // Handle bad channel ID
+    try {
+      wids = wireReadout->ChannelToWire(d->Channel());
+    }
+    catch(...) {
+      //Fail if something...fails...
+      throw cet::exception("CalibAnaTree_module") << "Raw digit " <<
+            d->Channel() << " has no wires";
+    }
+
+    //We shouldn't ever have a raw digit NOT mapped to a wire
+    if (wids.size() == 0) {
+      throw cet::exception("CalibAnaTree_module") << "Raw digit " <<
+            d->Channel() << " has no wires";
+    }
+
+    //A raw digit can come from multiple wire segments because some wires
+    //are wrapped
+    for (const auto & w : wids) {
+      // Ignore wires that are already mapped
+      if (rawdigits.count(w)) continue;
+      rawdigits[w] = d;
+    }
+  }
+
+
   if (fLowEnergyClusterAnalysis)
   {
     if (fVerbose)
@@ -297,6 +341,9 @@ void dune::CalibAnaTree::analyze(art::Event const& e)
       {
 	fCluster = vCluster[i];
 	fCluster->meta = fMeta;
+	//if (fVerbose) std::cout<< " tick window " << fRadiusExt/fElectronVelocity << " wire window " << fHitRawDigitsWireCollectWidth <<std::endl;
+	FillWireInfoForLECluster(fCluster , fWireReadout , rawdigits , fRadiusExt/fElectronVelocity , fHitRawDigitsWireCollectWidth); 
+	// fRadiusExt/fElectronVelocity is the size (in tt) of the full (cluster + veto)
         fTree_LE->Fill();
       }
     }
@@ -335,49 +382,6 @@ void dune::CalibAnaTree::analyze(art::Event const& e)
   art::FindManyP<recob::Hit, recob::TrackHitMeta> fmtrkHits(tracks, e, thm_label);
   art::FindManyP<anab::Calorimetry> fmCalo(tracks, e, fCALOproducer);
 
-  // Collect raw digits for saving hits
-  std::vector<art::Ptr<raw::RawDigit>> rawdigitlist;
-  for (const art::InputTag &t: fRawDigitproducers) {
-    try {
-      art::ValidHandle<std::vector<raw::RawDigit>> thisdigits = e.getValidHandle<std::vector<raw::RawDigit>>(t);
-      art::fill_ptr_vector(rawdigitlist, thisdigits);
-    }
-    catch(...) {
-      if (!fSilenceMissingDataProducts) throw;
-      else {} // Allow Raw Digits to not be present
-    }
-  }
-
-  // The raw digit list is not sorted, so make it into a map on the WireID
-  std::map<geo::WireID, art::Ptr<raw::RawDigit>> rawdigits;
-  for (const art::Ptr<raw::RawDigit> &d: rawdigitlist) {
-
-    std::vector<geo::WireID> wids;
-    // Handle bad channel ID
-    try {
-      wids = wireReadout->ChannelToWire(d->Channel());
-    }
-    catch(...) {
-      //Fail if something...fails...
-      throw cet::exception("CalibAnaTree_module") << "Raw digit " <<
-            d->Channel() << " has no wires";
-    }
-
-    //We shouldn't ever have a raw digit NOT mapped to a wire
-    if (wids.size() == 0) {
-      throw cet::exception("CalibAnaTree_module") << "Raw digit " <<
-            d->Channel() << " has no wires";
-    }
-
-    //A raw digit can come from multiple wire segments because some wires
-    //are wrapped
-    for (const auto & w : wids) {
-      // Ignore wires that are already mapped
-      if (rawdigits.count(w)) continue;
-      rawdigits[w] = d;
-    }
-  }
-
   // Collect all hits
   art::ValidHandle<std::vector<recob::Hit>> allhit_handle = e.getValidHandle<std::vector<recob::Hit>>(fHITproducer);
   std::vector<art::Ptr<recob::Hit>> allHits;