LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
OpticalRecoAna_module.cc
Go to the documentation of this file.
1 //Analyzer for flash<--->track matching (testing against MC)
2 #include "OpticalRecoAna.h"
3 
4 // LArSoft includes
6 
7 // Framework includes
9 
10 namespace opreco {
11  DEFINE_ART_MODULE(OpticalRecoAna)
12 }
13 
14 
15 // ART includes
19 
20 // ROOT includes
21 #include "TTree.h"
22 
23 // C++ Includes
24 #include <iostream>
25 #include <sstream>
26 
27  // Constructor
29 {
30 
31  // Indicate that the Input Module comes from .fcl
32  fFlashModuleLabel = pset.get<std::string>("FlashModule");
33  fTrackModuleLabel = pset.get<std::string>("TrackModule");
34  fKineticEnergyMin = pset.get<float>("KineticEnergyMin");
35  fPEMin = pset.get<float>("PEMin");
36  fTimeMatchMax = pset.get<float>("TimeMatchMax");
37 
38  fFlash_match_vector.clear();
39  fTrack_match_vector.clear();
40  fParticle_match_vector.clear();
41 }
42 
43 
44 
45 // Destructor
47 {}
48 
49 
50 
51 // Do something here to setup the file (like make a TTree)
53 {
55  fTimeDiff = tfs->make<TH1F>("htdiff","Time difference between particles and flashes; t_diff (ns); flash/particle pairs",1e3,-5e6,5e6);
56  fTimeDiff_fine = tfs->make<TH1F>("htdiff_fine","Time difference between particles and flashes; t_diff (ns); flash/particle pairs",100,-1000,1000);
57  fNBeamFlashes = tfs->make<TH1I>("hNBeamFlashes","Number of flashes OnBeamTime per event; N_{Flashes}; Events",5,0,5);
58 
59  fMatchTree_PF = tfs->make<TTree>("MatchTree_PF","MatchTree_PF");
60  fMatchTree_PF->Branch("Run",&fRun,"Run/I");
61  fMatchTree_PF->Branch("Event",&fEvent,"Event/I");
62  fMatchTree_PF->Branch("p_id",&fParticleID,"p_id/I");
63  fMatchTree_PF->Branch("p_time",&fParticleTime,"p_time/F");
64  fMatchTree_PF->Branch("p_vx",&fParticleVx,"p_vx/F");
65  fMatchTree_PF->Branch("p_vy",&fParticleVy,"p_vy/F");
66  fMatchTree_PF->Branch("p_vz",&fParticleVz,"p_vz/F");
67  fMatchTree_PF->Branch("p_mid",&fParticleMother,"p_mid/I");
68  fMatchTree_PF->Branch("p_trackid",&fParticleTrackID,"p_trackid/I");
69  fMatchTree_PF->Branch("FlashID",&fFlashID,"flash_id/I");
70  fMatchTree_PF->Branch("f_time",&fFlashTime,"f_time/F");
71  fMatchTree_PF->Branch("f_timewidth",&fFlashTimeWidth,"f_timewidth/F");
72  fMatchTree_PF->Branch("f_y",&fFlashY,"f_y/F");
73  fMatchTree_PF->Branch("f_z",&fFlashZ,"f_z/F");
74  fMatchTree_PF->Branch("f_ywidth",&fFlashYWidth,"f_ywidth/F");
75  fMatchTree_PF->Branch("f_zwidth",&fFlashYWidth,"f_zwidth/F");
76  fMatchTree_PF->Branch("f_onbeamtime",&fFlashOnBeamTime,"f_onbeamtime/I");
77  fMatchTree_PF->Branch("f_pe",&fFlashPE,"f_pe/F");
78  fMatchTree_PF->Branch("matchIndex",&fMatchIndex,"matchIndex/I");
79 
80  fMatchTree_PF_NotNu = tfs->make<TTree>("MatchTree_PF_NotNu","MatchTree_PF_NotNu");
81  fMatchTree_PF_NotNu->Branch("Run",&fRun,"Run/I");
82  fMatchTree_PF_NotNu->Branch("Event",&fEvent,"Event/I");
83  fMatchTree_PF_NotNu->Branch("p_id",&fParticleID,"p_id/I");
84  fMatchTree_PF_NotNu->Branch("p_time",&fParticleTime,"p_time/F");
85  fMatchTree_PF_NotNu->Branch("p_vx",&fParticleVx,"p_vx/F");
86  fMatchTree_PF_NotNu->Branch("p_vy",&fParticleVy,"p_vy/F");
87  fMatchTree_PF_NotNu->Branch("p_vz",&fParticleVz,"p_vz/F");
88  fMatchTree_PF_NotNu->Branch("p_mid",&fParticleMother,"p_mid/I");
89  fMatchTree_PF_NotNu->Branch("p_trackid",&fParticleTrackID,"p_trackid/I");
90  fMatchTree_PF_NotNu->Branch("FlashID",&fFlashID,"flash_id/I");
91  fMatchTree_PF_NotNu->Branch("f_time",&fFlashTime,"f_time/F");
92  fMatchTree_PF_NotNu->Branch("f_timewidth",&fFlashTimeWidth,"f_timewidth/F");
93  fMatchTree_PF_NotNu->Branch("f_y",&fFlashY,"f_y/F");
94  fMatchTree_PF_NotNu->Branch("f_z",&fFlashZ,"f_z/F");
95  fMatchTree_PF_NotNu->Branch("f_ywidth",&fFlashYWidth,"f_ywidth/F");
96  fMatchTree_PF_NotNu->Branch("f_zwidth",&fFlashYWidth,"f_zwidth/F");
97  fMatchTree_PF_NotNu->Branch("f_onbeamtime",&fFlashOnBeamTime,"f_onbeamtime/I");
98  fMatchTree_PF_NotNu->Branch("f_pe",&fFlashPE,"f_pe/F");
99  fMatchTree_PF_NotNu->Branch("matchIndex",&fMatchIndex_NotNu,"matchIndex/I");
100 
101 }
102 
103 
104 
105 // The analyzer itself
107 {
108 
109  //get flag for MC
110  const bool is_MC= !(evt.isRealData());
111 
112  fRun = evt.id().run();
113  fEvent = evt.id().event();
114 
115  //first get out track, flash, and particle list handles
117  evt.getByLabel(fFlashModuleLabel, flash_handle);
118  std::vector<recob::OpFlash> const& flash_vector(*flash_handle);
119  fFlash_match_vector.resize(flash_vector.size());
120 
121  LOG_INFO ("OpticalRecoAna")
122  << "Number of flashes is " << flash_vector.size() << std::flush;
123 
124  //art::Handle< std::vector<recob::Track> > track_handle;
125  //evt.getByLabel(fTrackModuleLabel, track_handle);
126  //std::vector<recob::Track> const& track_vector(*track_handle);
127  //fTrack_match_vector.resize(track_vector.size());
128 
129  //LOG_INFO ("OpticalRecoAna")
130  //<< "Number of tracks is " << track_vector.size() << std::flush;
131 
132  //match_flashes_to_tracks(flash_vector, track_vector);
133 
134  //all this for the MC matching
135  if(is_MC){
136 
138  std::vector<simb::MCParticle> particle_list;
139  get_MC_particle_list(pi_serv->ParticleList(),particle_list);
140  fParticle_match_vector.resize(particle_list.size());
141 
142  mf::LogInfo("OpticalRecoAna")
143  << "Number of MC particles is " << particle_list.size();
144 
145 
147  const float ns_per_PMT_tick = 1e3;// already converted to microseconds//( 1e3 / odp->SampleFreq()) ; //SampleFreq is in MHz
148  art::ServiceHandle<geo::Geometry> geometry_handle;
149 
150  match_flashes_to_particles(flash_vector,
151  particle_list,
152  ns_per_PMT_tick,
153  *geometry_handle);
154  /*
155  FillMatchTree_PF(flash_vector,
156  particle_list,
157  ns_per_PMT_tick,
158  *geometry_handle);
159  */
160  int n_flashes = 0;
161  for(auto const& my_flash : flash_vector){
162  if(my_flash.OnBeamTime()==1) n_flashes++;
163  }
164  fNBeamFlashes->Fill(n_flashes);
165 
166  //check_flash_matches();
167 
168  //match_tracks_to_particles(track_handle,particle_list);
169 
170  }
171 
172 }
173 
176  return (pi_serv->TrackIdToMCTruth_P(particle.TrackId()))->Origin();
177 }
178 
179 void opreco::OpticalRecoAna::get_MC_particle_list(sim::ParticleList const& plist,std::vector<simb::MCParticle> & particle_vector) {
180 
181  for(sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
182 
183  const simb::MCParticle* particle = (*ipart).second;
184  /*
185  std::cout << "Particle " << particle_vector.size() << " (" << (*ipart).first << ") is " << particle->PdgCode() << " (mother=" << particle->Mother()
186  << ") with K_energy " << particle->E()-0.001*particle->Mass() << std::endl;
187  std::cout << "\t(X,Y,Z,T)=(" << particle->Vx() << "," << particle->Vy() << "," << particle->Vz() << "," << particle->T() << ")" << std::endl;
188  */
189  //do not store if it's below our energy cut
190  if( particle->E() < 0.001*particle->Mass() + fKineticEnergyMin ) continue;
191 
192  //check to see if it's a charged particle we expect to leave ionization
193  const int pdg = std::abs(particle->PdgCode());
194  if( pdg==11 // electron
195  || pdg==13 //muon
196  || pdg==15 //tau
197  || pdg==211 //charged pion
198  || pdg==321 //charged kaon
199  || pdg==2212 //proton
200  || pdg==2214 //delta
201  || pdg==2224 //delta
202  || pdg==3222 //sigma
203  || pdg==3112 //sigma
204  || pdg==3312 //xi
205  || pdg==3334 //omega
206  ) {
207  particle_vector.emplace_back(*particle);
208  //std::cout << "Particle " << particle_vector.size() << " (" << (*ipart).first << ") is " << pdg << " (status=" << particle->StatusCode()
209  // << ") with K_energy " << particle->E()-0.001*particle->Mass() << std::endl;
210  //std::cout << "\t(X,Y,Z,T)=(" << particle->Vx() << "," << particle->Vy() << "," << particle->Vz() << "," << particle->T() << ")" << std::endl;
211  }
212  }
213 
214 }
215 
216 float opreco::OpticalRecoAna::update_MC_particle_time(simb::MCParticle const& particle, bool & pass_check, geo::Geometry const& geometry){
217 
218  pass_check = false;
219 
220  unsigned int tpc = 0;
221  unsigned int cstat = 0;
222 
223  const size_t numtraj = particle.NumberTrajectoryPoints();
224  size_t t_iter = 0;
225  while(t_iter < numtraj){
226  try{
227  // check if the particle is inside a TPC
228  double pos[3] = {particle.Vx(t_iter), particle.Vy(t_iter), particle.Vz(t_iter)};
229  geometry.PositionToTPC(pos, tpc, cstat);
230  }
231  catch(cet::exception &e){
232  t_iter++;
233  continue;
234  }
235  break;
236  }
237 
238  if(t_iter == numtraj)
239  return 0;
240 
241  pass_check = true;
242  return particle.T(t_iter);
243 }
244 
245 void opreco::OpticalRecoAna::match_flashes_to_tracks(std::vector<recob::OpFlash> const& flash_vector,
246  std::vector<recob::Track> const& track_vector){
247  bool matching=false;
248 
249  for(size_t i_flash=0; i_flash < flash_vector.size(); i_flash++){
250 
251  recob::OpFlash const& my_flash( flash_vector.at(i_flash) );
252  if(my_flash.TotalPE() < fPEMin) continue;
253 
254  for(size_t i_track=0; i_track < track_vector.size(); i_track++){
255 
256  recob::Track const& my_track( track_vector.at(i_track) );
257 
258  matching=false;
259  compare_track_and_flash(my_track,my_flash,matching);
260  if(matching){
261  fFlash_match_vector.at(i_flash).track_indices.push_back(i_track);
262  fTrack_match_vector.at(i_track).flash_indices.push_back(i_flash);
263  }
264 
265  }//end inner loop over tracks
266 
267  }//end loop over flashes
268 
269 }//end match_flashes_to_tracks
270 
272  recob::OpFlash const& flash,
273  bool & matching){
274 }
275 
276 void opreco::OpticalRecoAna::match_flashes_to_particles(std::vector<recob::OpFlash> const& flash_vector,
277  std::vector<simb::MCParticle> const& particle_vector,
278  float const& ns_per_PMT_tick,
279  geo::Geometry const& geometry){
280 
281  fMatchIndex=0;
283  int lastFlashID = -1;
284 
285  for(size_t i_particle=0; i_particle < particle_vector.size(); i_particle++){
286 
287  simb::MCParticle const& my_particle( particle_vector.at(i_particle) );
288  bool pass_check = false;
289  const float particle_time = update_MC_particle_time(my_particle,pass_check,geometry);
290  if(!pass_check) continue;
291 
292  //std::cout << "Particle " << i_particle << " (" << my_particle.PdgCode() << ")" << particle_time << std::endl;
293 
294  for(size_t i_flash=0; i_flash < flash_vector.size(); i_flash++){
295 
296  recob::OpFlash const& my_flash( flash_vector.at(i_flash) );
297  if(my_flash.TotalPE() < fPEMin) continue;
298 
299  const float flash_time = my_flash.Time()*ns_per_PMT_tick;
300 
301  //if(my_flash.OnBeamTime())
302  //std::cout << "\tFlash " << i_flash << " time is " << flash_time << std::endl;
303 
304  fTimeDiff->Fill(particle_time-flash_time);
305  fTimeDiff_fine->Fill(particle_time-flash_time);
306 
307  fParticleID = i_particle;
308  fParticleTime = particle_time;
309  fParticleVx = my_particle.Vx(0);
310  fParticleVy = my_particle.Vy(0);
311  fParticleVz = my_particle.Vz(0);
312  fParticleMother = my_particle.Mother();
313  fParticleTrackID = my_particle.TrackId();
314  fFlashID = i_flash;
315  fFlashTime = flash_time;
316  fFlashTimeWidth = my_flash.TimeWidth()*ns_per_PMT_tick;
317  fFlashY = my_flash.YCenter();
318  fFlashZ = my_flash.ZCenter();
319  fFlashYWidth = my_flash.YWidth();
320  fFlashZWidth = my_flash.ZWidth();
321  fFlashPE = my_flash.TotalPE();
322  fFlashOnBeamTime = my_flash.OnBeamTime();
323 
324  bool beam_match = ((std::abs(particle_time - flash_time )/fFlashTimeWidth)<fTimeMatchMax &&
326  bool nonbeam_match = ((std::abs(particle_time - flash_time )/fFlashTimeWidth)<fTimeMatchMax &&
328 
329  if( beam_match ){
330  fMatchTree_PF->Fill();
331  fMatchIndex++;
332 
333  //std::cout << "\t\tParticle (x,y,z)=(" << my_particle.Vx() << "," << my_particle.Vy() << "," << my_particle.Vz() << std::endl;
334  //std::cout << "\t\tParticle PDG = " << my_particle.PdgCode() << " mother=" << my_particle.Mother() << std::endl;
335  //std::cout << "\t\tFlash (y,z) = ("
336  // << my_flash.YCenter() << " +/- " << my_flash.YWidth() << ","
337  // << my_flash.ZCenter() << " +/- " << my_flash.ZWidth() << std::endl;
338 
339  fFlash_match_vector.at(i_flash).particle_indices.push_back(i_particle);
340  fParticle_match_vector.at(i_particle).flash_indices.push_back(i_flash);
341  }
342  else if( nonbeam_match ){
343  if(lastFlashID!=fFlashID){
344  fMatchTree_PF_NotNu->Fill();
346  lastFlashID = fFlashID;
347  }
348  }
349 
350 
351  }//end inner loop over particles
352 
353  }//end loop over flashes
354 
355 }//end match_flashes_to_particles
356 
357 void opreco::OpticalRecoAna::FillMatchTree_PF(std::vector<recob::OpFlash> const& flash_vector,
358  std::vector<simb::MCParticle> const& particle_vector,
359  float const& ns_per_PMT_tick,
360  geo::Geometry const& geometry){
361 
362  for(size_t i_flash=0; i_flash < flash_vector.size(); i_flash++){
363 
364  recob::OpFlash const& my_flash( flash_vector.at(i_flash) );
365  const float flash_time = my_flash.Time()*ns_per_PMT_tick;
366 
367  for(auto i_particle : fFlash_match_vector.at(i_flash).particle_indices){
368 
369  simb::MCParticle const& my_particle( particle_vector.at(i_particle) );
370  bool pass_check = false;
371  const float particle_time = update_MC_particle_time(my_particle,pass_check,geometry);
372 
373  if(my_particle.Mother()==0 && my_particle.TrackId()>1e4){
374 
375  fParticleID = i_particle;
376  fParticleTime = particle_time;
377  fParticleVx = my_particle.Vx(0);
378  fParticleVy = my_particle.Vy(0);
379  fParticleVz = my_particle.Vz(0);
380  fFlashID = i_flash;
381  fFlashTime = flash_time;
382  fFlashY = my_flash.YCenter();
383  fFlashZ = my_flash.ZCenter();
384  fFlashYWidth = my_flash.YWidth();
385  fFlashZWidth = my_flash.ZWidth();
386 
387  fMatchTree_PF->Fill();
388 
389  return;
390  }
391 
392  }
393 
394  }
395 
396 }
397 
398 //currently unused...
400  recob::OpFlash const& flash,
401  bool & matching,
402  float const& ns_per_PMT_tick,
403  geo::Geometry const& geometry){
404 }
405 
406 void opreco::OpticalRecoAna::match_flashes_to_particles(std::vector<recob::Track> const& track_vector,
407  std::vector<simb::MCParticle> const& particle_vector,
408  geo::Geometry const& geometry){
409  bool matching=false;
410 
411  for(size_t i_track=0; i_track < track_vector.size(); i_track++){
412 
413  recob::Track const& my_track( track_vector.at(i_track) );
414 
415  for(size_t i_particle=0; i_particle < particle_vector.size(); i_particle++){
416 
417  simb::MCParticle const& my_particle( particle_vector.at(i_particle) );
418 
419  matching=false;
420  compare_particle_and_track(my_particle,my_track,matching,geometry);
421  if(matching){
422  fTrack_match_vector.at(i_track).particle_indices.push_back(i_track);
423  fParticle_match_vector.at(i_particle).track_indices.push_back(i_track);
424  }
425 
426  }//end inner loop over particles
427 
428  }//end loop over tracks
429 
430 }//end match_tracks_to_particles
431 
433  recob::Track const& track,
434  bool & matching,
435  geo::Geometry const& geometry){
436 }
double E(const int i=0) const
Definition: MCParticle.h:237
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:222
int PdgCode() const
Definition: MCParticle.h:216
void match_flashes_to_particles(std::vector< recob::OpFlash > const &, std::vector< simb::MCParticle > const &, float const &, geo::Geometry const &)
#define LOG_INFO(category)
std::string fFlashModuleLabel
OpticalRecoAna(const fhicl::ParameterSet &)
void compare_track_and_flash(recob::Track const &, recob::OpFlash const &, bool &)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
double Mass() const
Definition: MCParticle.h:243
enum simb::_ev_origin Origin_t
event origin types
void analyze(const art::Event &)
std::vector< flash_match > fFlash_match_vector
bool isRealData() const
Definition: Event.h:83
RunNumber_t run() const
Definition: EventID.h:99
int TrackId() const
Definition: MCParticle.h:214
geo::TPCGeo const & PositionToTPC(geo::Point_t const &point) const
Returns the TPC at specified location.
simb::Origin_t get_MC_particle_origin(simb::MCParticle const &)
void FillMatchTree_PF(std::vector< recob::OpFlash > const &, std::vector< simb::MCParticle > const &, float const &, geo::Geometry const &)
std::string fTrackModuleLabel
double Time() const
Definition: OpFlash.h:83
void compare_particle_and_flash(simb::MCParticle const &, recob::OpFlash const &, bool &, float const &, geo::Geometry const &)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
T get(std::string const &key) const
Definition: ParameterSet.h:231
double T(const int i=0) const
Definition: MCParticle.h:228
iterator begin()
Definition: ParticleList.h:305
The geometry of one entire detector, as served by art.
Definition: Geometry.h:110
void match_flashes_to_tracks(std::vector< recob::OpFlash > const &, std::vector< recob::Track > const &)
double Vx(const int i=0) const
Definition: MCParticle.h:225
T * make(ARGS...args) const
std::vector< particle_match > fParticle_match_vector
void get_MC_particle_list(sim::ParticleList const &, std::vector< simb::MCParticle > &)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void compare_particle_and_track(simb::MCParticle const &, recob::Track const &, bool &, geo::Geometry const &)
Int_t ipart
Definition: Style.C:10
double Vz(const int i=0) const
Definition: MCParticle.h:227
EventNumber_t event() const
Definition: EventID.h:117
float update_MC_particle_time(simb::MCParticle const &, bool &, geo::Geometry const &)
TCEvent evt
Definition: DataStructs.cxx:5
Float_t e
Definition: plot.C:34
Float_t track
Definition: plot.C:34
std::vector< track_match > fTrack_match_vector
EventID id() const
Definition: Event.h:56
double Vy(const int i=0) const
Definition: MCParticle.h:226
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:52
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int const &id)
Beam neutrinos.
Definition: MCTruth.h:21