ClusterEva_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file ClusterEva.cxx
3 /// \brief TODO
4 /// \version $Id: ClusterEva.cxx,v 1.7 2012-09-25 06:30:12 ricken Exp $
5 /// \author ricken@fnal.gov
6 ////////////////////////////////////////////////////////////////////////
7 #include <cmath>
8 #include <vector>
9 #include <string>
10 
11 // ROOT includes
12 #include "TH1I.h"
13 #include "TNtuple.h"
14 #include "TBenchmark.h"
15 
16 // NOvA includes
17 #include "Calibrator/Calibrator.h"
19 #include "Geometry/Geometry.h"
21 #include "RecoBase/CellHit.h"
22 #include "RecoBase/Cluster.h"
23 #include "RecoBase/RecoHit.h"
24 #include "MCCheater/BackTracker.h"
26 
27 // Framework includes
34 
35 class TTree;
36 class TVector3;
37 
38 namespace geo { class Geometry; }
39 namespace rb { class CellHit; class Cluster; }
40 
41 ///tracking algorithms
42 namespace pa {
43  class ClusterEva : public art::EDAnalyzer {
44  public:
45  explicit ClusterEva(fhicl::ParameterSet const& pset);
46  virtual ~ClusterEva();
47 
48  void beginJob();
49  void analyze(const art::Event& evt);
50  void reconfigure(const fhicl::ParameterSet& p);
51 
52  private:
54  std::string fClusterInput; ///< Input folder from cluster reco
55  unsigned fNpart;
56  TTree* ClusterTruth;
57  bool _view;
58  double _tmin, _tmax;
59  int _nhits, _eventID;
60  int _ADC;
61  std::vector<float> _nTH;
62  std::vector<float> _nTHtot;
63  std::vector<float> _TADC;
64  std::vector<float> _TADCtot;
65  std::vector<float> _Ttmin;
66  std::vector<float> _Ttmax;
67  std::vector<float> _Tenergy;
68  };
69 }
70 
71 namespace pa{
72  //......................................................................
73  ClusterEva::ClusterEva(fhicl::ParameterSet const& pset)
74  : EDAnalyzer(pset)
75  {
76  reconfigure(pset);
77  }
78 
79  //......................................................................
81  {
82  }
83 
84  //......................................................................
86  {
87  fCellHitInput = pset.get< std::string >("CellHitInput");
88  fClusterInput = pset.get< std::string >("ClusterInput") ;
89 
90  }
91 
92  //......................................................................
94  {
96  ClusterTruth = tfs->make<TTree>("ClusterTruth","Hough Phase Space");
97 
98  ClusterTruth -> Branch("evtID", &_eventID, "evtID/I");
99  ClusterTruth -> Branch("view", &_view, "VIEW/O"); //0: X; 1: Y
100  ClusterTruth -> Branch("Tmin", &_tmin, "Tmin/D");
101  ClusterTruth -> Branch("Tmax", &_tmax, "Tmax/D");
102  ClusterTruth -> Branch("Nhits", &_nhits, "Nhits/I");
103  ClusterTruth -> Branch("totADC",&_ADC, "totADC/I");
104 
105  ClusterTruth -> Branch("TrueHits", &_nTH);
106  ClusterTruth -> Branch("totTrueHits", &_nTHtot);
107  ClusterTruth -> Branch("TrueADC", &_TADC);
108  ClusterTruth -> Branch("totTrueADC", &_TADCtot);
109  ClusterTruth -> Branch("TrueTmin", &_Ttmin);
110  ClusterTruth -> Branch("TrueTmax", &_Ttmax);
111  ClusterTruth -> Branch("TrueEnergy", &_Tenergy);
112  }
113 
115  {
116  _eventID = evt.id().event();
117 
118  //make a backtracker
120  // define vector holding the clusters
122  evt.getByLabel(fClusterInput, slices);
124  evt.getByLabel(fCellHitInput, hitcol);
125 
127 
128  for(unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx){
129  //Now you are inside a new cluster, you have to clear the vector:
130  _nTH.clear();
131  _nTHtot.clear();
132  _TADC.clear();
133  _TADCtot.clear();
134  _Ttmin.clear();
135  _Ttmax.clear();
136  _Tenergy.clear();
137 
138  const rb::Cluster& slice = (*slices)[sliceIdx];
139  _ADC = 0;
140  _tmin = (slice.Cell(0))->TNS();
141  _tmax =_tmin;
142  _nhits = slice.NCell();
143  if ((slice.Cell(0))-> View()==geo::kX) _view = 0;
144  else _view = 1;
145 
146  std::vector<int> track_ides;
147 
148  for(int hitIdx = 0; hitIdx != _nhits; ++hitIdx){
149 
150  const art::Ptr<rb::CellHit> chit = slice.Cell(hitIdx);
151  _ADC += chit->ADC();
152 
153  if (_tmin > chit->TNS() ) _tmin=chit->TNS();
154  if (_tmax < chit->TNS() ) _tmax=chit->TNS();
155 
156  //Now work on MC truth:
157  try {
158  if (!bt->IsNoise(chit)) {
159  const sim::Particle *part = bt->HitToParticle(chit);
160 
161  if ( find(track_ides.begin(),track_ides.end(),part->TrackId()) == track_ides.end())
162  {
163 
164  // std::cout<<"new particle track id:"<<part->TrackId()<<std::endl;
165  track_ides.push_back(part->TrackId() );
166  //Initialize items:
167  int cheater_nhits = 0;
168  int cheater_ADC = 0;
169  float cheater_tmin = 9.e19;
170  float cheater_tmax = -9.e10;
171 
172  //Now you have to loop over all the hits to check the true hits mainly caused by this particle:
173  for (unsigned i=0; i != hitcol->size(); ++i) {
174 
175  art::Ptr<rb::CellHit> hit_ptr(hitcol, i);
176  /*
177  if (!_view && hit_ptr->View()==geo::kY) continue;
178  else if (_view && hit_ptr->View()==geo::kX) continue;
179  */
180  // std::cout<<"pass view check"<<std::endl;
181  if (bt->HitToParticle(hit_ptr) == part) {
182  // std::cout<<"Get the particle ID:"<<(bt->HitToParticle(hit_ptr) )->TrackId()<<std::endl;
183  ++cheater_nhits;
184  cheater_ADC+=hit_ptr->ADC();
185 
186  // std::cout<<cheater_nhits<<"th hit with ADC "<<hit_ptr->ADC()<<std::endl;
187  if (cheater_tmin > hit_ptr->TNS() ) cheater_tmin = hit_ptr->TNS();
188  if (cheater_tmax < hit_ptr->TNS() ) cheater_tmax = hit_ptr->TNS();
189  }
190  }
191 
192  _nTH.push_back(1);
193  _TADC.push_back(chit->ADC());
194 
195  _Tenergy.push_back(part->E() );
196  _nTHtot.push_back(cheater_nhits);
197  _TADCtot.push_back(cheater_ADC);
198  _Ttmin.push_back(cheater_tmin);
199  _Ttmax.push_back(cheater_tmax);
200 
201  // std::cout<<"real energy of this particle: "<<part->E()<<" GeV"<<std::endl;
202 
203 
204  }
205  //Or this particle has been left trace in previous hits:
206  else {
207  int particle_idx = find(track_ides.begin(),track_ides.end(),part->TrackId()) - track_ides.begin();
208 
209  ++_nTH.at(particle_idx);
210  _TADC.at(particle_idx) += chit->ADC();
211 
212  }
213  }
214  }
215  catch (...) {
216  continue;
217  }
218  }//End looping over this slice
219 
220  //Fill the tree:
221  /*
222  std::cout<<"clustered hits Particle size: "<<_nTH.size()<<std::endl;
223  std::cout<<"clustered ADC Particle size: " <<_TADC.size()<<std::endl;
224  std::cout<<"total hits Particle size: " <<_nTHtot.size()<<std::endl;
225  std::cout<<"total ADC Particle size: " <<_TADC.size()<<std::endl;
226  std::cout<<"size of truth energy: " <<_Tenergy.size()<<std::endl;
227  */
228  ClusterTruth->Fill();
229 
230 
231  }// end looping over all slices
232 
233  return;
234 
235  } //end analyze
236 
238 
239 }
double E(const int i=0) const
Definition: MCParticle.h:232
float TNS() const
Definition: CellHit.h:46
std::vector< float > _TADC
std::string fClusterInput
Input folder from cluster reco.
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
tracking algorithms
std::vector< float > _Ttmin
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
DEFINE_ART_MODULE(TestTMapFile)
const sim::Particle * HitToParticle(art::Ptr< rb::CellHit > const &hit, bool quiet=false) const
Returns the sim::Particle object contributing the most light to the specified rb::CellHit.
std::string fCellHitInput
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int TrackId() const
Definition: MCParticle.h:209
TString part[npart]
Definition: Style.C:32
std::vector< float > _Tenergy
void beginJob()
std::vector< float > _TADCtot
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
Perform a "2 point" Hough transform on a collection of hits.
Definition: View.py:1
void reconfigure(const fhicl::ParameterSet &p)
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
std::vector< float > _nTHtot
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
std::vector< float > _nTH
std::vector< float > _Ttmax
Helper for AttenCurve.
Definition: Path.h:10
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
void analyze(const art::Event &evt)