SuperDDTEva_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file SuperDDTEva.cxx
3 //
4 // \brief First official DDT analyzer in offline!
5 //
6 // \author ricken@fnal.gov
7 ////////////////////////////////////////////////////////////////////////
15 #include "MCCheater/BackTracker.h"
17 
26 
27 #include "DAQChannelMap/DAQChannelMap.h"
28 #include "RecoBase/Track.h"
29 #include "RecoBase/Cluster.h"
30 #include "Simulation/FLSHit.h"
31 #include "Simulation/FLSHitList.h"
32 
33 
34 #include "TVector3.h"
35 #include "TLorentzVector.h"
36 #include "TNtuple.h"
37 #include "TH1I.h"
38 
39 #include <string>
40 #include <iostream>
41 #include <cmath>
42 #include <algorithm>
43 #include <vector>
44 #include <fstream>
45 
46 class TTree;
47 
48 namespace ddt {
49  class SuperDDTEva;
50  struct Chit;
51 }
52 
53 struct ddt::Chit
54 {
55  Chit(unsigned short p, unsigned short c, uint32_t t, int16_t a, int ssid, int tsid, int rnsid, int trid) :
56  plane(p), cell(c), tdc(t), adc(a), ssliceid(ssid), tsliceid(tsid), rnsliceid(rnsid), trackid(trid) {}
57 
58  unsigned short plane, cell;
59  uint32_t tdc;
60  int16_t adc;
62 };
63 
64 struct trueHit
65 {
66  trueHit(unsigned short p, unsigned short c, double t, bool f, double e, int _pdg) :
67  plane(p), cell(c), texit(t), found(f), energy(e), P_D_G(_pdg) {}
68  unsigned short plane, cell;
69  double texit;
70  bool found;
71  double energy;
72  int P_D_G;
73 };
74 
76 {
77 public:
78  explicit SuperDDTEva(fhicl::ParameterSet const &pset);
79  virtual ~SuperDDTEva();
80 
81  void analyze(const art::Event &evt);
82  void beginJob();
83  void reconfigure(fhicl::ParameterSet const &pset);
84 
85 private:
86  TTree* MCTruth_Tree;
87  TTree* HitList_Tree;
88  TTree* Track_Tree;
89  TTree* TrueHits_Tree;
90  /*
91  TH1I* nhits_hist;
92  TH1I* ntracks_hist;
93  */
94  TH1I* n_RNSlice;
95  TH1I* n_SpaceSlice;
96 
97  //Truth info
98  int pdg;
99  int nFLShits;
100  double Px_truth;
101  double Py_truth;
102  double Pz_truth;
105  double x0_truth;
106  double y0_truth;
107  double z0_truth;
108  double t0_truth;
109  double x1_truth;
110  double y1_truth;
111  double z1_truth;
112  double t1_truth;
113 
114  //Track info
115  int view;
116  double v0_ddt;
117  double z0_ddt;
118  double t0_ddt;
119  double v1_ddt;
120  double z1_ddt;
121  double t1_ddt;
122 
123  //HitList Info:
124  int Plane, Cell,ADC;
125  uint32_t TDC;
126  int trackID, TsliceID, RNsliceID, SsliceID,eventID;
127  bool isNoise, isMM;
128 
131 
134 
137 
141 
142  //True Hits Info
143  int Plane_true, Cell_true;
144  double Time_true, Energy_true;
145  bool Found_true;
146 
147 };
148 
149 
150 
152  : EDAnalyzer(pset)
153 {
154  reconfigure(pset);
155 }
156 
157 
158 //......................................................................
160 {
161  TimeSlice_ModuleLabel = pset.get< std::string >("TimeSlice_ModuleLabel" );
162  TimeSlice_InstanceLabel = pset.get< std::string >("TimeSlice_InstanceLabel" );
163  RNSlice_ModuleLabel = pset.get< std::string >("RNSlice_ModuleLabel" );
164  RNSlice_InstanceLabel = pset.get< std::string >("RNSlice_InstanceLabel" );
165  SpaceSlice_ModuleLabel = pset.get< std::string >("SpaceSlice_ModuleLabel" );
166  SpaceSlice_InstanceLabel = pset.get< std::string >("SpaceSlice_InstanceLabel" );
167  OnlineTracker_ModuleLabel = pset.get< std::string >("OnlineTracker_ModuleLabel" );
168  OnlineTracker_InstanceLabel = pset.get< std::string >("OnlineTracker_InstanceLabel" );
169  fFLSHitListLabel = pset.get< std::string >("FLSHitListLabel") ;
170 }
171 
172 //......................................................................
173 
175 
176 
177 //......................................................................
178 
180 {
182  n_RNSlice = tfs->make<TH1I>("n_RNSlice","number of remove noise slice",1000,-0.5,999.5);
183  n_SpaceSlice = tfs->make<TH1I>("n_SpaceSlice","number of space slice",100,-0.5,99.5);
184 
185  MCTruth_Tree = tfs->make<TTree>("MCTruth","MC_info");
186  // MCTruth_Tree -> Branch("PDG",&pdg,"PDG/I");
187  MCTruth_Tree -> Branch("nhits",&nFLShits,"nhits/I");
188  MCTruth_Tree -> Branch("Px",&Px_truth,"Px/D");
189  MCTruth_Tree -> Branch("Py",&Py_truth,"Py/D");
190  MCTruth_Tree -> Branch("Pz",&Pz_truth,"Pz/D");
191  MCTruth_Tree -> Branch("EnergyDump",&energydump_truth,"energydump/D");
192  MCTruth_Tree -> Branch("PathLength",&pathlength_truth,"pathlength/D");
193  MCTruth_Tree -> Branch("x0",&x0_truth,"x0/D");
194  MCTruth_Tree -> Branch("y0",&y0_truth,"y0/D");
195  MCTruth_Tree -> Branch("z0",&z0_truth,"z0/D");
196  MCTruth_Tree -> Branch("t0",&t0_truth,"t0/D");
197  MCTruth_Tree -> Branch("x1",&x1_truth,"x1/D");
198  MCTruth_Tree -> Branch("y1",&y1_truth,"y1/D");
199  MCTruth_Tree -> Branch("z1",&z1_truth,"z1/D");
200  MCTruth_Tree -> Branch("t1",&t1_truth,"t1/D");
201  MCTruth_Tree -> Branch("evtID",&eventID,"evtID/I");
202 
203  HitList_Tree = tfs->make<TTree>("HitList","HitList_info");
204  HitList_Tree -> Branch("Plane",&Plane,"Plane/I");
205  HitList_Tree -> Branch("Cell",&Cell,"Cell/I");
206  HitList_Tree -> Branch("TDC",&TDC,"TDC/i");
207  HitList_Tree -> Branch("ADC",&ADC,"ADC/I");
208  HitList_Tree -> Branch("trackID",&trackID,"trackID/I");
209  HitList_Tree -> Branch("TimesliceID",&TsliceID,"TimesliceID/I");
210  HitList_Tree -> Branch("RNsliceID",&RNsliceID,"RNsliceID/I");
211  HitList_Tree -> Branch("SsliceID",&SsliceID,"SsliceID/I");
212  HitList_Tree -> Branch("eventID",&eventID,"eventID/I");
213  HitList_Tree -> Branch("IsNoise",&isNoise,"IsNoise/O");
214  HitList_Tree -> Branch("IsMM",&isMM,"IsMM/O");
215 
216  Track_Tree = tfs->make<TTree>("Track_Tree","Track_info");
217  Track_Tree -> Branch("View",&view,"View/I");
218  Track_Tree -> Branch("v0",&v0_ddt,"v0/D");
219  Track_Tree -> Branch("z0",&z0_ddt,"z0/D");
220  Track_Tree -> Branch("t0",&t0_ddt,"t0/D");
221  Track_Tree -> Branch("v1",&v1_ddt,"v1/D");
222  Track_Tree -> Branch("z1",&z1_ddt,"z1/D");
223  Track_Tree -> Branch("t1",&t1_ddt,"t1/D");
224  Track_Tree -> Branch("eventID",&eventID,"eventID/I");
225 
226  TrueHits_Tree = tfs->make<TTree>("TrueHits_Tree","TrueHits_info");
227  TrueHits_Tree -> Branch("Plane",&Plane_true,"Plane/I");
228  TrueHits_Tree -> Branch("Cell",&Cell_true,"Cell/I");
229  TrueHits_Tree -> Branch("Energy",&Energy_true,"Energy/D");
230  TrueHits_Tree -> Branch("Time",&Time_true,"Time/D");
231  TrueHits_Tree -> Branch("Found",&Found_true,"Found/O");
232  TrueHits_Tree -> Branch("eventID",&eventID,"eventID/I");
233 
234 
235  /*
236  nhits_hist = tfs->make<TH1I>("nhitsHist","Number of Hits per Slice",1000,0,1000);
237  ntracks_hist = tfs->make<TH1I>("ntracks","Number of Tracks per Slice",50,0,50);
238  */
239 }
240 
242 {
243  eventID = evt.id().event();
244  std::vector<Chit> chits;
245  std::vector<trueHit> thits;
246 
247 
248  //Get Truth Info
249  //Make a backtracker
251  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
252 
253  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
254  const sim::Particle *p = (*i).second;
255  pdg = p->PdgCode();
256  Px_truth = p->Px();
257  Py_truth = p->Py();
258  Pz_truth = p->Pz();
259 
260  if (bt->ParticleToFLSHit(p->TrackId()).size()==0) continue;
261 
262  //Later info needs to be obtained from the FLSHit:
263  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
264 
265  pathlength_truth = flshits[0].GetTotalPathLength();
266  //pathlength_truth = 0;
267  energydump_truth = flshits[0].GetTotalEnergyLoss();
268  x0_truth = flshits[0].GetEntryX();
269  y0_truth = flshits[0].GetEntryY();
270  z0_truth = flshits[0].GetEntryZ();
271  t0_truth = flshits[0].GetEntryT();
272  x1_truth = flshits[0].GetExitX();
273  y1_truth = flshits[0].GetExitY();
274  z1_truth = flshits[0].GetExitZ();
275  t1_truth = flshits[0].GetExitT();
276 
277  thits.emplace_back(flshits[0].GetPlaneID(),flshits[0].GetCellID(),t1_truth,false,energydump_truth, pdg);
278 
279  nFLShits = flshits.size();
280  for (int h =1; h!= nFLShits; ++h) {
281 
282  // std::cout<<"plane: "<<flshits[h].GetPlaneID()<<" cell: "<<flshits[h].GetCellID()<<" time: "<<t1_truth<<" pdg: "<<pdg<<std::endl;
283  thits.emplace_back(flshits[h].GetPlaneID(),flshits[h].GetCellID(),t1_truth,false,flshits[h].GetTotalEnergyLoss(), pdg);
284  // std::cout<<"t_true:"<<t1_truth<<std::endl;
285  // std::cout<<"energy: "<<flshits[h].GetTotalEnergyLoss()<<std::endl;
286 
287  pathlength_truth+=flshits[h].GetTotalPathLength();
288  energydump_truth+=flshits[h].GetTotalEnergyLoss();
289  if (t0_truth > flshits[h].GetEntryT() ) {
290  x0_truth = flshits[h].GetEntryX();
291  y0_truth = flshits[h].GetEntryY();
292  z0_truth = flshits[h].GetEntryZ();
293  t0_truth = flshits[h].GetEntryT();
294  }
295  else if (t1_truth < flshits[h].GetExitT() ) {
296  x1_truth = flshits[h].GetExitX();
297  y1_truth = flshits[h].GetExitY();
298  z1_truth = flshits[h].GetExitZ();
299  t1_truth = flshits[h].GetExitT();
300  }
301  }
302  MCTruth_Tree->Fill();
303  }
304 
305 
306  //Get TimeSlice info:
310  if (ghl.n_groups()==0) return;
311  // std::cout<<"number of original time slices: "<< ghl.n_groups()<<std::endl;
312  for (unsigned i =0; i!= ghl.n_groups();++i) {
313  novaddt::HitList hits = ghl.get(i);
314 
315  for (auto hit_ptr = hits.begin(); hit_ptr!=hits.end();++hit_ptr){
316  chits.emplace_back(hit_ptr->Plane().val, hit_ptr->Cell().val,
317  hit_ptr->TDC().val, hit_ptr->ADC().val,
318  0,i+1, 0, 0 //spacesliceID, timesliceID, rnsliceID, trackID
319  );
320  }
321 
322  }
323 
324  //Flip the RNSliceID:
328  n_RNSlice -> Fill(RNghl.n_groups());
329  // std::cout<<"number of removenoise slice: "<< RNghl.n_groups()<<std::endl;
330  for (unsigned i =0; i!= RNghl.n_groups(); ++i) {
331  novaddt::HitList hits = RNghl.get(i);
332  for (auto hit_ptr = hits.begin(); hit_ptr!=hits.end();++hit_ptr){
333  for (auto chit_ptr= chits.begin(); chit_ptr!=chits.end();++chit_ptr ) {
334  if (chit_ptr->rnsliceid) continue;
335  if (chit_ptr->plane == hit_ptr->Plane().val &&
336  chit_ptr->cell == hit_ptr->Cell().val &&
337  chit_ptr->tdc == hit_ptr->TDC().val) {
338  chit_ptr->rnsliceid =i+1;
339  break;
340  }
341  }
342  }
343  }
344 
345 
346  //Flip the SpaceSliceID:
347  art::Handle<std::vector<novaddt::HitList> > hit_lists; // from HitList
349  n_SpaceSlice->Fill(hit_lists->size());
350  for (unsigned i=0; i!=hit_lists->size(); ++i) {
351  novaddt::HitList const &hits = hit_lists->at(i);
352  for (auto hit_ptr = hits.begin(); hit_ptr!=hits.end();++hit_ptr){
353  for (auto chit_ptr= chits.begin(); chit_ptr!=chits.end();++chit_ptr ) {
354  if (chit_ptr->ssliceid > 0) continue;
355  if (chit_ptr->plane == hit_ptr->Plane().val &&
356  chit_ptr->cell == hit_ptr->Cell().val &&
357  chit_ptr->tdc == hit_ptr->TDC().val) {
358  chit_ptr->ssliceid =i+1;
359  break;
360  }
361  }
362  }
363  }
364  /*
365  //Get Track Info
366  art::Handle<std::vector<novaddt::Track> > tracks;
367  evt.getByLabel(OnlineTracker_ModuleLabel,OnlineTracker_InstanceLabel, tracks);
368 
369  art::FindOneP<novaddt::HitList> fohl(tracks, evt, OnlineTracker_ModuleLabel);
370  for (unsigned i =0; i != tracks->size(); ++i) {
371  view = tracks->at(i).View();
372  v0_ddt = tracks->at(i).StartV();
373  z0_ddt = tracks->at(i).StartZ();
374  v1_ddt = tracks->at(i).EndV();
375  z1_ddt = tracks->at(i).EndZ();
376  t0_ddt = tracks->at(i).StartT();
377  t1_ddt = tracks->at(i).EndT();
378 
379  Track_Tree->Fill();
380  //Get associated hit list:
381  art::Ptr<novaddt::HitList> hl_ptr = fohl.at(i);
382  for (auto const& hit : (*hl_ptr) ) {
383  //loop over the vector of chits:
384  for (auto chit_ptr= chits.begin(); chit_ptr!=chits.end();++chit_ptr ) {
385  if (chit_ptr->trackid) continue;
386  if (chit_ptr->plane == hit.Plane().val &&
387  chit_ptr->cell == hit.Cell().val &&
388  chit_ptr->tdc == hit.TDC().val) {
389  chit_ptr->trackid =i+1;
390  break;
391  }
392  }
393  }
394 
395  }
396  */
397  // std::cout<<"Here is the list of the found true hits:"<<std::endl;
398 
399  //Fill the hitlist tree:
400  for (auto const& chit : chits) {
401  Plane = chit.plane;
402  Cell = chit.cell;
403  TDC = chit.tdc;
404  ADC = chit.adc;
405  TsliceID = chit.tsliceid;
406  RNsliceID = chit.rnsliceid;
407  SsliceID = chit.ssliceid;
408  trackID = chit.trackid;
409  isNoise = true;
410  isMM = false;
411  for (auto true_hit_ptr= thits.begin(); true_hit_ptr != thits.end(); ++true_hit_ptr) {
412  if (true_hit_ptr->plane == Plane && true_hit_ptr->cell== Cell) {
413 
414  if ( std::max(TDC,uint32_t(true_hit_ptr->texit / 15.625))
415  -std::min(TDC,uint32_t(true_hit_ptr->texit / 15.625)) <70 ) {
416  isNoise = false;
417 
418  // std::cout<<"plane: "<<Plane<<" cell:"<<Cell<<std::endl;
419 
420  // if (true_hit_ptr->found) std::cout<<"This hit has already been assigned to a DAQHit!"<<std::endl;
421  true_hit_ptr->found = true;
422  if (true_hit_ptr->P_D_G == 42) {
423  isMM = true;
424  break;
425  }
426  }
427  }
428  }
429 
430  HitList_Tree -> Fill();
431  }
432 
433  for (auto const& true_hit: thits) {
434  Plane_true = true_hit.plane;
435  Cell_true = true_hit.cell;
436  Energy_true = true_hit.energy;
437  Found_true = true_hit.found;
438  Time_true = true_hit.texit;
439  // std::cout<<"energy dump: "<< true_hit.energy<<std::endl;
440  TrueHits_Tree-> Fill();
441  }
442 
443  return;
444 }
445 
446 
T max(const caf::Proxy< T > &a, T b)
HitList get(unsigned group) const
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
void analyze(const art::Event &evt)
std::string RNSlice_ModuleLabel
double Py(const int i=0) const
Definition: MCParticle.h:230
Chit(unsigned short p, unsigned short c, uint32_t t, int16_t a, int ssid, int tsid, int rnsid, int trid)
trueHit(unsigned short p, unsigned short c, double t, bool f, double e, int _pdg)
void reconfigure(fhicl::ParameterSet const &pset)
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
std::string SpaceSlice_ModuleLabel
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
list_type::const_iterator const_iterator
double Px(const int i=0) const
Definition: MCParticle.h:229
Int_t eventID
Definition: plot.C:81
DEFINE_ART_MODULE(TestTMapFile)
int TrackId() const
Definition: MCParticle.h:209
std::string fFLSHitListLabel
void hits()
Definition: readHits.C:15
void beginJob()
unsigned short plane
std::string OnlineTracker_InstanceLabel
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
std::string RNSlice_InstanceLabel
double energy
Definition: plottest35.C:25
std::string TimeSlice_InstanceLabel
SuperDDTEva(fhicl::ParameterSet const &pset)
unsigned short cell
unsigned n_groups() const
EventNumber_t event() const
Definition: EventID.h:116
unsigned short plane
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double Pz(const int i=0) const
Definition: MCParticle.h:231
std::string SpaceSlice_InstanceLabel
T min(const caf::Proxy< T > &a, T b)
Int_t trackID
Definition: plot.C:84
Float_t e
Definition: plot.C:35
std::string OnlineTracker_ModuleLabel
EventID id() const
Definition: Event.h:56
std::string TimeSlice_ModuleLabel
enum BeamMode string