TrackingEfficiency_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: TrackingEfficiency
3 // Module Type: analyzer
4 // File: TrackingEfficiency_module.cc
5 //
6 // Generated at Sun Mar 24 09:05:02 2013 by Tingjun Yang using artmod
7 // from art v1_02_06.
8 //
9 // Using AnaTree as base, have written a module to compare reconstruction
10 // algorithm efficiencies.
11 //
12 // Thomas Karl Warburton
13 // k.warburton@sheffield.ac.uk
14 //
15 ////////////////////////////////////////////////////////////////////////
16 
17 // Framework includes
21 #include "fhiclcpp/ParameterSet.h"
26 #include "art_root_io/TFileService.h"
27 #include "art_root_io/TFileDirectory.h"
29 
30 // LArSoft includes
49 
52 
53 // ROOT includes
54 #include "TTree.h"
55 #include "TTimeStamp.h"
56 #include "TLorentzVector.h"
57 #include "TH2F.h"
58 #include "TEfficiency.h"
59 #include "TNtuple.h"
60 #include "TFile.h"
61 
62 //standard library includes
63 #include <map>
64 #include <iostream>
65 #include <iomanip>
66 #include <sstream>
67 #include <cmath>
68 #include <memory>
69 #include <limits> // std::numeric_limits<>
70 
71 namespace TrackingEfficiency {
72  class TrackingEfficiency;
73 }
74 
76 public:
77  explicit TrackingEfficiency(fhicl::ParameterSet const & p);
78  virtual ~TrackingEfficiency();
79 
80  void analyze(art::Event const & e) override;
81 
82  void beginRun(art::Run const& run) override;
83  void beginJob() override;
84  void endJob() override;
85  void endRun(art::Run const&) override;
86  //void reconfigure(fhicl::ParameterSet const & p) ;
87 
88  // ----------- Declare my structs ----------
89  struct EfficHists {
90  //EfficHists();
91  // EfficHists(const std::string& subdir);
92  TH1D *Length;
93  TH1D *Energy;
94  TH1D *Energy_Depos;
95  TH1D *Theta;
96  TH1D *Theta_XZ;
97  TH1D *Theta_YZ;
98  TH1D *Phi;
99  TH2D *Phi_v_Theta;
100  TH1D *Eta_XY;
101  TH1D *Eta_ZY;
103  struct TEfficiencies {
104  //TEfficiencies();
105  //TEfficiencies(const std::string& subdir);
106  TEfficiency *Effic_Length;
107  TEfficiency *Effic_Energy;
108  TEfficiency *Effic_Energy_Depos;
109  TEfficiency *Effic_Theta;
110  TEfficiency *Effic_Theta_XZ;
111  TEfficiency *Effic_Theta_YZ;
112  TEfficiency *Effic_Phi;
113  TEfficiency *Effic_Phi_v_Theta;
114  TEfficiency *Effic_Eta_XY;
115  TEfficiency *Effic_Eta_ZY;
117  // ----------- Declare my structs ----------
118 
119 private:
120 
121  void ResetVars();
122  void MCTruthInformation ( detinfo::DetectorClocksData const& clockData,
123  const simb::MCParticle *particle, double &Energy, double &EnergyDeposited, double &TPCLength,
124  double &Theta_XZ, double &Theta_YZ, double &Eta_XY, double &Eta_ZY, double &Theta, double &Phi,
125  int &MCPdgCode, int &MCTrackId );
127  double Energy, double EnergyDeposited, double TPCLength,
128  double Theta_XZ, double Theta_YZ, double Eta_XY, double Eta_ZY, double Theta, double Phi );
129  void Make_Efficiencies ( struct TrackingEfficiency::TrackingEfficiency::TEfficiencies *AllHistPtr, bool matched,
130  double Energy, double EnergyDeposited, double TPCLength,
131  double Theta_XZ, double Theta_YZ, double Eta_XY, double Eta_ZY, double Theta, double Phi );
132  void EfficPlot_1D ( TH1D *AllHist, TH1D *MatchHist, TEfficiency *EfficHist );
133  void EfficPlot_2D ( TH2D *AllHist, TH2D *MatchHist, TEfficiency *EfficHist );
134 
135  TTree* fTree;
136 
140  double TrackLength;
141 
146 
148  // Handles
151 
153  double WindowSize;
154  // Parameter List
159 
160  };
161 // ********************************** Begin Run *******************************************************
163 
164 }
165 // *********************************** Begin Job ********************************************************
167 {
168  // Implementation of optional member function here.
170  fTree = tfs->make<TTree>("TrackingTree","analysis tree");
171  fTree->Branch("MCTPCLength" ,&MCTPCLength ,"MCTPCLength/D" );
172  fTree->Branch("MCEnergy" ,&MCEnergy ,"MCEnergy/D" );
173  fTree->Branch("MCEnergyDeposited",&MCEnergyDeposited,"MCEnergyDeposited/D");
174  fTree->Branch("MCTheta" ,&MCTheta ,"MCTheta/D" );
175  fTree->Branch("MCTheta_XZ" ,&MCTheta_XZ ,"MCTheta_XZ/D" );
176  fTree->Branch("MCTheta_YZ" ,&MCTheta_YZ ,"MCTheta/D" );
177  fTree->Branch("MCPhi" ,&MCPhi ,"MCPhi/D" );
178  fTree->Branch("MCEta_XY" ,&MCEta_XY ,"MCEta_XY/D" );
179  fTree->Branch("MCEta_ZY" ,&MCEta_ZY ,"MCEta_ZY/D" );
180  fTree->Branch("MCMatched" ,&MCMatched ,"MCMatched/B" );
181  fTree->Branch("MCPdgCode" ,&MCPdgCode ,"MCPdgCode/I" );
182  fTree->Branch("MCTrackId" ,&MCTrackId ,"MCTrackId/I" );
183  fTree->Branch("MatchedTrackID" ,&MatchedTrackID ,"MatchedTrackID/I" );
184 
185 
186 
187  MCAllEffic.Length = tfs->make<TH1D>("Length_MCAll" ,"All charged Monte Carlo particle track lengths;Length (cm);Efficiency", 50, 0, 500 );
188  MCAllEffic.Energy = tfs->make<TH1D>("Energy_MCAll" ,"All charged Monte Carlo particle energies;Energy (GeV);Efficiency" , 50, 0, 20 );
189  MCAllEffic.Energy_Depos = tfs->make<TH1D>("Energy_Depos_MCAll","All charged Monte Carlo particle energy deposited;Energy Deposited (GeV);Efficiency",50, 0, 3.5 );
190  MCAllEffic.Theta = tfs->make<TH1D>("Theta_MCAll" ,"All charged Monte Carlo particle theta;Theta (rad);Efficiency", 50, 0, TMath::Pi() );
191  MCAllEffic.Theta_XZ = tfs->make<TH1D>("Theta_XZ_MCAll" ,"All charged Monte Carlo particle thetaXZ;ThetaxZ (rad);Efficiency",50, -TMath::Pi(), TMath::Pi() );
192  MCAllEffic.Theta_YZ = tfs->make<TH1D>("Theta_YZ_MCAll" ,"All charged Monte Carlo particle thetaYZ;ThetaYZ (rad);Efficiency",50, -TMath::Pi(), 0 );
193  MCAllEffic.Phi = tfs->make<TH1D>("Phi_MCAll" ,"All charged Monte Carlo particle phi;Phi (rad);Efficiency", 50, -TMath::Pi(), 0 );
194  MCAllEffic.Eta_XY = tfs->make<TH1D>("Eta_XY_MCAll" ,"All charged Monte Carlo particle etaXZ;EtaxY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
195  MCAllEffic.Eta_ZY = tfs->make<TH1D>("Eta_ZY_MCAll" ,"All charged Monte Carlo particle etaYZ;EtaZY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
196  MCAllEffic.Phi_v_Theta = tfs->make<TH2D>("Phi_v_Theta_MCAll" ,"All charged Monte Carlo phi versus theta;Theta (rad);Phi (rad)", 50, 0, TMath::Pi(), 50, -TMath::Pi(), 0);
197 
198  MatchedAllEffic.Length = tfs->make<TH1D>("Length_MatchAll" ,"Matched all charged Monte Carlo particle track lengths;Length (cm);Efficiency",50, 0, 500 );
199  MatchedAllEffic.Energy = tfs->make<TH1D>("Energy_MatchAll" ,"Matched all charged Monte Carlo particle energies;Energy (GeV);Efficiency" ,50, 0, 20 );
200  MatchedAllEffic.Energy_Depos = tfs->make<TH1D>("Energy_Depos_MatchAll","Matched all charged Monte Carlo particle energy deposited;Energy Deposited (GeV);Efficiency",50, 0, 3.5 );
201  MatchedAllEffic.Theta = tfs->make<TH1D>("Theta_MatchAll" ,"Matched all charged Monte Carlo particle theta;Theta (rad);Efficiency", 50, 0, TMath::Pi() );
202  MatchedAllEffic.Theta_XZ = tfs->make<TH1D>("Theta_XZ_MatchAll" ,"Matched all charged Monte Carlo particle thetaXZ;ThetaxZ (rad);Efficiency",50, -TMath::Pi(), TMath::Pi() );
203  MatchedAllEffic.Theta_YZ = tfs->make<TH1D>("Theta_YZ_MatchAll" ,"Matched all charged Monte Carlo particle thetaYZ;ThetaYZ (rad);Efficiency",50, -TMath::Pi(), 0 );
204  MatchedAllEffic.Phi = tfs->make<TH1D>("Phi_MatchAll" ,"Matched all charged Monte Carlo particle phi;Phi (rad);Efficiency", 50, -TMath::Pi(), 0 );
205  MatchedAllEffic.Eta_XY = tfs->make<TH1D>("Eta_XY_MatchAll" ,"Matched all charged Monte Carlo particle etaXZ;EtaxY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
206  MatchedAllEffic.Eta_ZY = tfs->make<TH1D>("Eta_ZY_MatchAll" ,"Matched all charged Monte Carlo particle etaYZ;EtaZY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
207  MatchedAllEffic.Phi_v_Theta = tfs->make<TH2D>("Phi_v_Theta_MatchAll" ,"Matched all charged Monte Carlo phi versus theta;Theta (rad);Phi (rad)", 50, 0, TMath::Pi(), 50, -TMath::Pi(), 0);
208 
209  AllEfficiencies.Effic_Length = tfs->make<TEfficiency>("Effic_Length_All" ,"Efficiency vs Particle Length for all charged particles;Length (cm);Efficiency", 50, 0, 500 );
210  AllEfficiencies.Effic_Energy = tfs->make<TEfficiency>("Effic_Energy_All" ,"Efficiency vs Particle Energy for all charged particles;Energy (GeV);Efficiency", 50, 0, 20 );
211  AllEfficiencies.Effic_Energy_Depos = tfs->make<TEfficiency>("Effic_Energy_Depos_All","Efficiency vs Particle Energy Deposited for all charged particles;Energy Deposited (GeV);Efficiency",50, 0, 3.5 );
212  AllEfficiencies.Effic_Theta = tfs->make<TEfficiency>("Effic_Theta_All" ,"Efficiency vs Particle Theta for all charged particles;Theta (rad);Efficiency", 50, 0, TMath::Pi() );
213  AllEfficiencies.Effic_Theta_XZ = tfs->make<TEfficiency>("Effic_Theta_XZ_All" ,"Efficiency vs Particle ThetaXZ for all charged particles;ThetaxZ (rad);Efficiency",50, -TMath::Pi(), TMath::Pi() );
214  AllEfficiencies.Effic_Theta_YZ = tfs->make<TEfficiency>("Effic_Theta_YZ_All" ,"Efficiency vs Particle ThetaYZ for all charged particles;ThetaYZ (rad);Efficiency",50, -TMath::Pi(), 0 );
215  AllEfficiencies.Effic_Phi = tfs->make<TEfficiency>("Effic_Phi_All" ,"Efficiency vs Particle Phi for all charged particles;Phi (rad);Efficiency", 50, -TMath::Pi(), 0 );
216  AllEfficiencies.Effic_Eta_XY = tfs->make<TEfficiency>("Effic_Eta_XY_All" ,"Efficiency vs Particle EtaXZ for all charged particles;EtaxY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
217  AllEfficiencies.Effic_Eta_ZY = tfs->make<TEfficiency>("Effic_Eta_ZY_All" ,"Efficiency vs Particle EtaYZ for all charged particles;EtaZY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
218  AllEfficiencies.Effic_Phi_v_Theta = tfs->make<TEfficiency>("Effic_Phi_v_Theta_All" ,"Efficiency of Phi versus theta for all charged particles;Theta (rad);Phi (rad)", 50, 0, TMath::Pi(), 50, -TMath::Pi(), 0);
219 
220  // --------------------- Now for Protons ---------------------------
221  MCProtonEffic.Length = tfs->make<TH1D>("Length_MCProton" ,"Proton charged Monte Carlo particle track lengths;Length (cm);Efficiency", 50, 0, 500 );
222  MCProtonEffic.Energy = tfs->make<TH1D>("Energy_MCProton" ,"Proton charged Monte Carlo particle energies;Energy (GeV);Efficiency" , 50, 0, 20 );
223  MCProtonEffic.Energy_Depos = tfs->make<TH1D>("Energy_Depos_MCProton","Proton charged Monte Carlo particle energy deposited;Energy Deposited (GeV);Efficiency",50, 0, 3.5 );
224  MCProtonEffic.Theta = tfs->make<TH1D>("Theta_MCProton" ,"Proton charged Monte Carlo particle theta;Theta (rad);Efficiency", 50, 0, TMath::Pi() );
225  MCProtonEffic.Theta_XZ = tfs->make<TH1D>("Theta_XZ_MCProton" ,"Proton charged Monte Carlo particle thetaXZ;ThetaxZ (rad);Efficiency",50, -TMath::Pi(), TMath::Pi() );
226  MCProtonEffic.Theta_YZ = tfs->make<TH1D>("Theta_YZ_MCProton" ,"Proton charged Monte Carlo particle thetaYZ;ThetaYZ (rad);Efficiency",50, -TMath::Pi(), 0 );
227  MCProtonEffic.Phi = tfs->make<TH1D>("Phi_MCProton" ,"Proton charged Monte Carlo particle phi;Phi (rad);Efficiency", 50, -TMath::Pi(), 0 );
228  MCProtonEffic.Eta_XY = tfs->make<TH1D>("Eta_XY_MCProton" ,"Proton charged Monte Carlo particle etaXZ;EtaxY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
229  MCProtonEffic.Eta_ZY = tfs->make<TH1D>("Eta_ZY_MCProton" ,"Proton charged Monte Carlo particle etaYZ;EtaZY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
230  MCProtonEffic.Phi_v_Theta = tfs->make<TH2D>("Phi_v_Theta_MCProton" ,"Proton charged Monte Carlo phi versus theta;Theta (rad);Phi (rad)", 50, 0, TMath::Pi(), 50, -TMath::Pi(), 0);
231 
232  MatchedProtonEffic.Length = tfs->make<TH1D>("Length_MatchProton" ,"Matched all charged Monte Carlo particle track lengths;Length (cm);Efficiency",50, 0, 500 );
233  MatchedProtonEffic.Energy = tfs->make<TH1D>("Energy_MatchProton" ,"Matched all charged Monte Carlo particle energies;Energy (GeV);Efficiency" ,50, 0, 20 );
234  MatchedProtonEffic.Energy_Depos = tfs->make<TH1D>("Energy_Depos_MatchProton","Matched all charged Monte Carlo particle energy deposited;Energy Deposited (GeV);Efficiency",50, 0, 3.5 );
235  MatchedProtonEffic.Theta = tfs->make<TH1D>("Theta_MatchProton" ,"Matched all charged Monte Carlo particle theta;Theta (rad);Efficiency", 50, 0, TMath::Pi() );
236  MatchedProtonEffic.Theta_XZ = tfs->make<TH1D>("Theta_XZ_MatchProton" ,"Matched all charged Monte Carlo particle thetaXZ;ThetaxZ (rad);Efficiency",50, -TMath::Pi(), TMath::Pi() );
237  MatchedProtonEffic.Theta_YZ = tfs->make<TH1D>("Theta_YZ_MatchProton" ,"Matched all charged Monte Carlo particle thetaYZ;ThetaYZ (rad);Efficiency",50, -TMath::Pi(), 0 );
238  MatchedProtonEffic.Phi = tfs->make<TH1D>("Phi_MatchProton" ,"Matched all charged Monte Carlo particle phi;Phi (rad);Efficiency", 50, -TMath::Pi(), 0 );
239  MatchedProtonEffic.Eta_XY = tfs->make<TH1D>("Eta_XY_MatchProton" ,"Matched all charged Monte Carlo particle etaXZ;EtaxY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
240  MatchedProtonEffic.Eta_ZY = tfs->make<TH1D>("Eta_ZY_MatchProton" ,"Matched all charged Monte Carlo particle etaYZ;EtaZY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
241  MatchedProtonEffic.Phi_v_Theta = tfs->make<TH2D>("Phi_v_Theta_MatchProton" ,"Matched all charged Monte Carlo phi versus theta;Theta (rad);Phi (rad)", 50, 0, TMath::Pi(), 50, -TMath::Pi(), 0);
242 
243  ProtonEfficiencies.Effic_Length = tfs->make<TEfficiency>("Effic_Length_Proton" ,"Efficiency vs Particle Length for all charged particles;Length (cm);Efficiency", 50, 0, 500 );
244  ProtonEfficiencies.Effic_Energy = tfs->make<TEfficiency>("Effic_Energy_Proton" ,"Efficiency vs Particle Energy for all charged particles;Energy (GeV);Efficiency", 50, 0, 20 );
245  ProtonEfficiencies.Effic_Energy_Depos = tfs->make<TEfficiency>("Effic_Energy_Depos_Proton","Efficiency vs Particle Energy Deposited for all charged particles;Energy Deposited (GeV);Efficiency",50, 0, 3.5 );
246  ProtonEfficiencies.Effic_Theta = tfs->make<TEfficiency>("Effic_Theta_Proton" ,"Efficiency vs Particle Theta for all charged particles;Theta (rad);Efficiency", 50, 0, TMath::Pi() );
247  ProtonEfficiencies.Effic_Theta_XZ = tfs->make<TEfficiency>("Effic_Theta_XZ_Proton" ,"Efficiency vs Particle ThetaXZ for all charged particles;ThetaxZ (rad);Efficiency",50, -TMath::Pi(), TMath::Pi() );
248  ProtonEfficiencies.Effic_Theta_YZ = tfs->make<TEfficiency>("Effic_Theta_YZ_Proton" ,"Efficiency vs Particle ThetaYZ for all charged particles;ThetaYZ (rad);Efficiency",50, -TMath::Pi(), 0 );
249  ProtonEfficiencies.Effic_Phi = tfs->make<TEfficiency>("Effic_Phi_Proton" ,"Efficiency vs Particle Phi for all charged particles;Phi (rad);Efficiency", 50, -TMath::Pi(), 0 );
250  ProtonEfficiencies.Effic_Eta_XY = tfs->make<TEfficiency>("Effic_Eta_XY_Proton" ,"Efficiency vs Particle EtaXZ for all charged particles;EtaxY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
251  ProtonEfficiencies.Effic_Eta_ZY = tfs->make<TEfficiency>("Effic_Eta_ZY_Proton" ,"Efficiency vs Particle EtaYZ for all charged particles;EtaZY (rad);Efficiency", 50, -TMath::Pi(), TMath::Pi() );
252  ProtonEfficiencies.Effic_Phi_v_Theta = tfs->make<TEfficiency>("Effic_Phi_v_Theta_Proton" ,"Efficiency of Phi versus theta for all charged particles;Theta (rad);Phi (rad)", 50, 0, TMath::Pi(), 50, -TMath::Pi(), 0);
253 
254 }
255 // ************************************ End Job *********************************************************
257 }
258 // ************************************ End Run *********************************************************
260 }
261 
262 // ********************************** pset param *******************************************************
264  : EDAnalyzer(pset)
265  , fHitsModuleLabel ( pset.get< std::string >("HitsModuleLabel"))
266  , fTrackModuleLabel ( pset.get< std::string >("TrackModuleLabel"))
267  , fMCTruthT0ModuleLabel ( pset.get< std::string >("MCTruthT0ModuleLabel"))
268  , fPhotonT0ModuleLabel ( pset.get< std::string >("PhotonT0ModuleLabel"))
269 {
270  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
271  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
272  XDriftVelocity = detProp.DriftVelocity()*1e-3; //cm/ns
273  WindowSize = detProp.NumberTimeSamples() * clockData.TPCClock().TickPeriod() * 1e3;
274 }
275 // ******************************************************************************************************
277 {
278  // Clean up dynamic memory and other resources here.
279 
280 }
281 // ************************************ Analyse *********************************************************
283 {
284  //std::cout << "\n\n************* New Event / Module running *************\n\n" << std::endl;
285  // Implementation of required member function here.
286  std::vector<art::Ptr<recob::Track> > tracklist;
287  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
288  if (trackListHandle)
289  art::fill_ptr_vector(tracklist, trackListHandle);
290 
291  const sim::ParticleList& plist = pi_serv->ParticleList();
292 
293  auto trackh = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
294 
295  auto trackvh = evt.getHandle< std::vector< art::PtrVector < recob::Track > > >(fTrackModuleLabel);
296 
297  int NPart = 0;
298 
299  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
300  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
301  if ( trackListHandle.isValid() ) { // Check that trackListHandle is valid.....
302  art::FindManyP<recob::Hit> fmht (trackListHandle, evt, fTrackModuleLabel);
303  art::FindMany<anab::T0> fmt0 (trackListHandle, evt, fMCTruthT0ModuleLabel);
304  art::FindMany<anab::T0> fmphot(trackListHandle, evt, fPhotonT0ModuleLabel);
305  int ntracks_reco=tracklist.size();
306 
307  // **************************************************************
308  // Want to loop through Monte Carlo Truth Particles.
309  // **************************************************************
310  for ( sim::ParticleList::const_iterator ipar = plist.begin(); ipar!=plist.end(); ++ipar){
311  simb::MCParticle *particle = ipar->second;
312 
313  MCMatched = false; // MCParticle not matched by defualt
314  AllChargedTrack = false; // MCParticle not charged by default
315  MatchedTrackID = -1; // No reco track has this ID
316 
317  // ----- Check that have a charged particle -----
318  if ( fabs(particle->PdgCode() ) == 13 ) AllChargedTrack = true;
319  else if ( particle->PdgCode() == 2212 ) AllChargedTrack = true;
320  /*
321  else if ( fabs(particle->PdgCode() == 11) ) AllChargedTrack = true; // Electron
322  else if ( fabs(particle->PdgCode() == 211) ) AllChargedTrack = true; // Pions
323  else if ( fabs(particle->PdgCode() == 321) ) AllChargedTrack = true; // Kaons
324  //*/
325  if (!AllChargedTrack) continue;
326 
327  // ---- Get MCTruth Information and check that MCParticle goes in TPC ---
328  MCTruthInformation ( clockData, particle, MCEnergy, MCEnergyDeposited, MCTPCLength,
329  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi,
330  MCPdgCode, MCTrackId );
331  if ( MCTPCLength < 1. ) continue;
332  ++NPart;
333  //std::cout << "Looking at MCParticle " << MCTrackId << ", with PdgCode " << MCPdgCode << ", MCTPCLength " << MCTPCLength << std::endl;
334 
335  // ----- Fill the histograms for all the Monte Carlo Particles -----
336  HistoFiller ( &MCAllEffic, MCEnergy, MCEnergyDeposited, MCTPCLength,
337  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi );
338  if ( MCPdgCode == 2212 ) {
339  HistoFiller ( &MCProtonEffic, MCEnergy, MCEnergyDeposited, MCTPCLength,
340  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi );
341  } // Proton
342 
343  for(int Track=0; Track < ntracks_reco; ++Track){
344  if ( MCMatched ) continue;
345  // *****************************************************************************
346  // T0 stuff - So can correct X Start/End positions and identify MCParticle!!
347  // *****************************************************************************
348  double TickT0 = 0;
349  if ( fmt0.isValid() ) {
350  std::vector<const anab::T0*> T0s = fmt0.at(Track);
351  for (size_t t0size =0; t0size < T0s.size(); t0size++) {
352  MCTruthT0 = T0s[t0size]->Time();
353  MCTruthTickT0 = MCTruthT0 / sampling_rate(clockData);
354  MCTruthTrackID = T0s[t0size]->TriggerBits();
355  } // T0 size
356  TickT0 = MCTruthTickT0;
357  } // T0 valid
358  if ( fmphot.isValid() ) {
359  std::vector<const anab::T0*> PhotT0 = fmphot.at(Track);
360  for (size_t T0it=0; T0it<PhotT0.size(); ++T0it) {
361  PhotonCounterT0 = PhotT0[T0it]->Time();
362  PhotonCounterTickT0 = PhotonCounterT0 / sampling_rate(clockData);
363  PhotonCounterID = PhotT0[T0it]->TriggerBits();
364  }
365  //TickT0 = PhotonCounterTickT0;
366  }
367  // ************** END T0 stuff ***************
368 
369  // ***** Check whether MCParticle made this track *****
370  if ( MCTrackId != MCTruthTrackID ) continue;
371 
372  // ******************************************************************************************
373  // Correct X and get track length etc now that we have matched a Track with an MCParticle!!
374  // ******************************************************************************************
375  // ----- Hit Stuff ------
376  std::vector< art::Ptr<recob::Hit> > allHits = fmht.at(Track);
377  double Hit_Size = allHits.size();
378  // ---- Correct X positions and get track positions.
379  recob::Track::Point_t trackStart, trackEnd;
380  std::tie(trackStart, trackEnd) = tracklist[Track]->Extent();
381  trackStart.SetX( trackStart.X() - detProp.ConvertTicksToX( TickT0, allHits[Hit_Size-1]->WireID().Plane, allHits[Hit_Size-1]->WireID().TPC, allHits[Hit_Size-1]->WireID().Cryostat )); // Correct X, last entry is first 'hit'
382  trackEnd.SetX( trackEnd.X() - detProp.ConvertTicksToX( TickT0, allHits[0]->WireID().Plane, allHits[0]->WireID().TPC, allHits[0]->WireID().Cryostat)); // Correct X, first entry is last 'hit'
383  // ---- Get lengths and angles.
384  art::Ptr<recob::Track> ptrack(trackh, Track);
385  const recob::Track& track = *ptrack;
386  TrackLength = track.Length();
387  int ntraj = track.NumberTrajectoryPoints();
388  if(ntraj > 0) {
389  auto dir = track.VertexDirection();
390  TrackTheta_XZ = std::atan2(dir.X(), dir.Z());
391  TrackTheta_YZ = std::atan2(dir.Y(), dir.Z());
392  TrackEta_XY = std::atan2(dir.X(), dir.Y());
393  TrackEta_ZY = std::atan2(dir.Z(), dir.Y());
394  TrackTheta = dir.Theta();
395  TrackPhi = dir.Phi();
396  }
397  // *********** Corrected X and track information **************
398 
399  // ***************************** Check how well reconstrcuted ******************************
400  //std::cout << "Track " << tracklist[Track]->ID() << " had reco Length " << TrackLength << ", MCLength " << MCTPCLength << " ---> " << TrackLength / MCTPCLength << std::endl;
401  if ( TrackLength / MCTPCLength > 0.5 &&
402  TrackLength / MCTPCLength < 1.5 &&
403  MCMatched != true && // If first time this MCParticle passed criteria
404  1 ) {
405 
406  MCMatched = true; // Stop MCParticle passing criteria with multiple tracks
407  MatchedTrackID = tracklist[Track]->ID(); // Which trackid matched this MCParticle?
408  //std::cout << "Particle matched so fill numerator and denominator\n" << std::endl;
409  // --- Track is well reconstructed, so fill the all charged histograms and TEfficiencies!!! ----
410  HistoFiller ( &MatchedAllEffic, MCEnergy, MCEnergyDeposited, MCTPCLength,
411  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi );
412  Make_Efficiencies ( &AllEfficiencies, true,
413  MCEnergy, MCEnergyDeposited, MCTPCLength,
414  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi);
415 
416  // --- If is a proton fill that efficiency plot too ----
417  if ( MCPdgCode == 2212 ) {
418  HistoFiller ( &MatchedProtonEffic, MCEnergy, MCEnergyDeposited, MCTPCLength,
419  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi );
420  Make_Efficiencies ( &ProtonEfficiencies, true,
421  MCEnergy, MCEnergyDeposited, MCTPCLength,
422  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi);
423  } // Proton
424 
425  } // If well matched track
426  // ***************************** Check how well reconstrcuted ******************************
427 
428  } // Track Loop
429  if (!MCMatched) { // If charged MCParticle doesn't get matched want to fill just denominator.
430  //std::cout << "MCParticle not matched so only fill denominator...\n" << std::endl;
431  Make_Efficiencies ( &AllEfficiencies, false,
432  MCEnergy, MCEnergyDeposited, MCTPCLength,
433  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi);
434  // --- If a proton, also want to fill that one too! ---
435  if ( MCPdgCode == 2212 ) {
436  Make_Efficiencies ( &ProtonEfficiencies, false,
437  MCEnergy, MCEnergyDeposited, MCTPCLength,
438  MCTheta_XZ, MCTheta_YZ, MCEta_XY, MCEta_ZY, MCTheta, MCPhi);
439  } // Proton
440  } // If not matched
441  // ******** Fill Tree for each MCParticle **********
442  fTree->Fill();
443  } // Loop over MCParticles
444  } // if trackListHandle.isValid()
445  //std::cout << "\nThis event had " << NPart << " particles." << std::endl;
446 }
447 // *********************************** Monte Carlo Truth Extraction ********************************************************
449  const simb::MCParticle *particle, double &Energy, double &EnergyDeposited, double &TPCLength,
450  double &Theta_XZ, double &Theta_YZ, double &Eta_XY, double &Eta_ZY, double &Theta, double &Phi,
451  int &MCPdgCode, int &MCTrackId ) {
452  int numberTrajectoryPoints = particle->NumberTrajectoryPoints(); // Looking at each MC hit
453  //double TPCLengthHits[numberTrajectoryPoints];
454  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0);
455  bool BeenInVolume = false;
456  int FirstHit=0, LastHit=0;
457  TPCLength = 0;
458 
459  MCPdgCode = particle->PdgCode();
460  MCTrackId = particle->TrackId();
461 
462  for(unsigned int MCHit=0; MCHit < TPCLengthHits.size(); ++MCHit) {
463  const TLorentzVector& tmpPosition=particle->Position(MCHit);
464  double const tmpPosArray[]={tmpPosition[0],tmpPosition[1],tmpPosition[2]};
465 
466  if (MCHit!=0) TPCLengthHits[MCHit] = pow ( pow( (particle->Vx(MCHit-1)-particle->Vx(MCHit)),2)
467  + pow( (particle->Vy(MCHit-1)-particle->Vy(MCHit)),2)
468  + pow( (particle->Vz(MCHit-1)-particle->Vz(MCHit)),2)
469  , 0.5 );
470  // --- Check if hit is in TPC...
471  geo::TPCID tpcid = geom->FindTPCAtPosition(tmpPosArray);
472  if (tpcid.isValid) {
473  // -- Check if hit is within drift window...
474  geo::CryostatGeo const& cryo = geom->Cryostat(tpcid.Cryostat);
475  geo::TPCGeo const& tpc = cryo.TPC(tpcid.TPC);
476  double XPlanePosition = tpc.PlaneLocation(0)[0];
477  double DriftTimeCorrection = fabs( tmpPosition[0] - XPlanePosition ) / XDriftVelocity;
478  double TimeAtPlane = particle->T() + DriftTimeCorrection;
479  if ( TimeAtPlane < trigger_offset(clockData)
480  || TimeAtPlane > trigger_offset(clockData) + WindowSize
481  ) continue;
482  // -- Good hit in TPC
483  LastHit = MCHit;
484  if ( !BeenInVolume ) {
485  BeenInVolume = true;
486  FirstHit = MCHit;
487  }
488  } // TPC.valid
489  /*
490  if (MCTrackId == 54 ) {
491  std::cout << MCHit << " of " << numberTrajectoryPoints << " x,y,z " << tmpPosition[0] << " " << tmpPosition[1] << " " << tmpPosition[2] << " " << TPCLengthHits[MCHit] << ", t " << particle->T() << ". In volume? " << BeenInVolume << " first/last hit " << FirstHit << " " << LastHit << std::endl;
492  }
493  //*/
494  } // MCTrajPoints
495 
496  Energy = particle->E(FirstHit);
497  EnergyDeposited = particle->E(FirstHit) - particle->E(LastHit);
498  for (int Hit = FirstHit+1; Hit <= LastHit; ++Hit ) {
499  TPCLength += TPCLengthHits[Hit];
500  //if (MCTrackId == 54 ) std::cout << FirstHit << " " << LastHit << ", " << Hit << ". dx " << TPCLengthHits[Hit] << " Len " << TPCLength << std::endl;
501  }
502  TLorentzVector& momentumStart = (TLorentzVector&)particle->Momentum(FirstHit);
503  TVector3 mcstartmom = particle->Momentum(FirstHit).Vect();
504  Theta_XZ = std::atan2(momentumStart.Px(), momentumStart.Pz());
505  Theta_YZ = std::atan2(momentumStart.Py(), momentumStart.Pz());
506  Eta_XY = std::atan2(momentumStart.Px(), momentumStart.Py());
507  Eta_ZY = std::atan2(momentumStart.Pz(), momentumStart.Py());
508  Theta = mcstartmom.Theta();
509  Phi = mcstartmom.Phi();
510 
511 } // MCTruthInformation
512 
513 // ******************************** Efficiency Filler ****************************************************
515  double Energy, double EnergyDeposited, double TPCLength,
516  double Theta_XZ, double Theta_YZ, double Eta_XY, double Eta_ZY, double Theta, double Phi ) {
517  HistPtr->Length -> Fill( TPCLength );
518  HistPtr->Energy -> Fill( Energy );
519  HistPtr->Energy_Depos -> Fill( EnergyDeposited );
520  HistPtr->Theta -> Fill( Theta );
521  HistPtr->Theta_XZ -> Fill( Theta_XZ );
522  HistPtr->Theta_YZ -> Fill( Theta_YZ );
523  HistPtr->Phi -> Fill( Phi );
524  HistPtr->Phi_v_Theta -> Fill( Theta, Phi );
525  HistPtr->Eta_XY -> Fill( Eta_XY );
526  HistPtr->Eta_ZY -> Fill( Eta_ZY );
527 } // HistoFiller
528 
529 // ******************************** Make Efficiencies ****************************************************
531  double Energy, double EnergyDeposited, double TPCLength,
532  double Theta_XZ, double Theta_YZ, double Eta_XY, double Eta_ZY, double Theta, double Phi ) {
533  AllHistPtr->Effic_Length -> Fill ( matched, TPCLength );
534  AllHistPtr->Effic_Energy -> Fill ( matched, Energy );
535  AllHistPtr->Effic_Energy_Depos-> Fill ( matched, EnergyDeposited );
536  AllHistPtr->Effic_Theta -> Fill ( matched, Theta );
537  AllHistPtr->Effic_Theta_XZ -> Fill ( matched, Theta_XZ );
538  AllHistPtr->Effic_Theta_YZ -> Fill ( matched, Theta_YZ );
539  AllHistPtr->Effic_Phi -> Fill ( matched, Phi );
540  AllHistPtr->Effic_Eta_XY -> Fill ( matched, Eta_XY );
541  AllHistPtr->Effic_Eta_ZY -> Fill ( matched, Eta_ZY );
542 
543  AllHistPtr->Effic_Phi_v_Theta -> Fill ( matched, Theta, Phi );
544 }
545 // ******************************** Reset Variables *****************************************************
547 }
548 // ******************************** Define Module *****************************************************
double E(const int i=0) const
Definition: MCParticle.h:233
void analyze(art::Event const &e) override
struct TrackingEfficiency::TrackingEfficiency::EfficHists MCProtonEffic
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
art::ServiceHandle< cheat::ParticleInventoryService > pi_serv
int PdgCode() const
Definition: MCParticle.h:212
struct TrackingEfficiency::TrackingEfficiency::EfficHists MatchedProtonEffic
Encapsulate the construction of a single cyostat.
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
struct TrackingEfficiency::TrackingEfficiency::EfficHists MCAllEffic
constexpr T pow(T x)
Definition: pow.h:72
struct TrackingEfficiency::TrackingEfficiency::EfficHists MatchedAllEffic
Geometry information for a single TPC.
Definition: TPCGeo.h:38
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
Vector_t VertexDirection() const
Definition: Track.h:132
void MCTruthInformation(detinfo::DetectorClocksData const &clockData, const simb::MCParticle *particle, double &Energy, double &EnergyDeposited, double &TPCLength, double &Theta_XZ, double &Theta_YZ, double &Eta_XY, double &Eta_ZY, double &Theta, double &Phi, int &MCPdgCode, int &MCTrackId)
STL namespace.
intermediate_table::const_iterator const_iterator
void HistoFiller(struct TrackingEfficiency::TrackingEfficiency::EfficHists *HistPtr, double Energy, double EnergyDeposited, double TPCLength, double Theta_XZ, double Theta_YZ, double Eta_XY, double Eta_ZY, double Theta, double Phi)
string dir
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
Particle class.
struct TrackingEfficiency::TrackingEfficiency::TEfficiencies AllEfficiencies
art framework interface to geometry description
Definition: Run.h:17
int TrackId() const
Definition: MCParticle.h:210
void EfficPlot_2D(TH2D *AllHist, TH2D *MatchHist, TEfficiency *EfficHist)
TrackingEfficiency(fhicl::ParameterSet const &p)
const double e
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
void EfficPlot_1D(TH1D *AllHist, TH1D *MatchHist, TEfficiency *EfficHist)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
double T(const int i=0) const
Definition: MCParticle.h:224
p
Definition: test.py:223
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
struct TrackingEfficiency::TrackingEfficiency::TEfficiencies ProtonEfficiencies
Encapsulate the geometry of a wire.
void Make_Efficiencies(struct TrackingEfficiency::TrackingEfficiency::TEfficiencies *AllHistPtr, bool matched, double Energy, double EnergyDeposited, double TPCLength, double Theta_XZ, double Theta_YZ, double Eta_XY, double Eta_ZY, double Theta, double Phi)
double Vx(const int i=0) const
Definition: MCParticle.h:221
Declaration of signal hit object.
art::ServiceHandle< geo::Geometry > geom
Encapsulate the construction of a single detector plane.
Contains all timing reference information for the detector.
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
tracking::Point_t Point_t
Definition: Track.h:53
Provides recob::Track data product.
double Vz(const int i=0) const
Definition: MCParticle.h:223
int trigger_offset(DetectorClocksData const &data)
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:1036
TCEvent evt
Definition: DataStructs.cxx:7
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Definition: TrackState.h:17
void beginRun(art::Run const &run) override
const double * PlaneLocation(unsigned int p) const
Definition: TPCGeo.cxx:382
double Vy(const int i=0) const
Definition: MCParticle.h:222
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:49
Encapsulate the construction of a single detector plane.