LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
NeutrinoTrackingEff_module.cc
Go to the documentation of this file.
1 //
2 //**Tracking Efficiency module***
3 //The basic idea is to loop over the hits from a given track and call BackTracker
4 //then look at std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackID(hit);
5 //then associete the hits to a G4 track ID (particle) that generate those hits(electrons)
6 //It was developed for CC neutrio interactions, it also can handle proton decay events p->k+nu_bar
7 //And protons, pion and muons from particle cannon by using the option isNeutrinoInt = false;
8 //All histograms are as a function of truth info (momentum, length)
9 //
10 // A. Higuera
11 // ahiguera@central.uh.edu
12 
13 #ifndef NeutrinoTrackingEff_Module
14 #define NeutrinoTrackingEff_Module
15 
16 // LArSoft includes
17 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
18 #include "lardata/Utilities/GeometryUtilities.h"
19 #include "larsim/Simulation/LArG4Parameters.h"
20 #include "lardataobj/RecoBase/Track.h"
21 #include "nusimdata/SimulationBase/MCParticle.h"
22 #include "nusimdata/SimulationBase/MCTruth.h"
23 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
24 #include "larsim/MCCheater/BackTrackerService.h"
25 #include "larsim/MCCheater/ParticleInventoryService.h"
26 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
27 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
28 #include "lardataobj/RawData/ExternalTrigger.h"
29 
30 // Framework includes
31 #include "art/Framework/Core/EDAnalyzer.h"
32 #include "art/Framework/Principal/Event.h"
33 #include "art/Framework/Principal/Handle.h"
34 #include "art/Framework/Services/Registry/ServiceHandle.h"
35 #include "art/Framework/Services/Optional/TFileService.h"
36 #include "art/Framework/Core/ModuleMacros.h"
37 #include "canvas/Persistency/Common/FindManyP.h"
38 #include "canvas/Persistency/Common/PtrVector.h"
39 #include "messagefacility/MessageLogger/MessageLogger.h"
40 #include "fhiclcpp/ParameterSet.h"
41 
42 // ROOT includes
43 #include "TFile.h"
44 #include "TTree.h"
45 #include "TDirectory.h"
46 #include "TH1.h"
47 #include "TEfficiency.h"
48 #include "TGraphAsymmErrors.h"
49 
50 // C/C++ standard libraries
51 #include <vector>
52 
53 #define MAX_TRACKS 1000
54 using namespace std;
55 
56 //========================================================================
57 
58 namespace DUNE{
59 
60 class NeutrinoTrackingEff : public art::EDAnalyzer {
61 public:
62 
63  explicit NeutrinoTrackingEff(fhicl::ParameterSet const& pset);
64  virtual ~NeutrinoTrackingEff();
65 
66  void beginJob();
67  void endJob();
68  void beginRun(const art::Run& run);
69  void analyze(const art::Event& evt);
70 
71  void reconfigure(fhicl::ParameterSet const& pset);
72 
73  void processEff(const art::Event& evt, bool &isFiducial);
74  void truthMatcher( std::vector<art::Ptr<recob::Hit>> all_hits, std::vector<art::Ptr<recob::Hit>> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet);
75  double truthLength( const simb::MCParticle *MCparticle );
76  bool insideFV(double vertex[4]);
77  void doEfficiencies();
78 
79 private:
80 
81  // the parameters we'll read from the .fcl
82  std::string fMCTruthModuleLabel;
83  std::string fTrackModuleLabel;
84  int fNeutrinoPDGcode;
85  int fLeptonPDGcode;
86  double fMaxNeutrinoE;
87  double fMaxLeptonP;
88  bool fisNeutrinoInt;
89 
90  int MC_isCC;
91  int MC_incoming_PDG;
92  double MC_incoming_P[4];
93  double MC_vertex[4];
94  double MC_lepton_startMomentum[4];
95 
96  int MC_leading_protonID;
97  int MC_leading_PionPlusID;
98  int MC_leading_PionMinusID;
99  int MC_leptonID;
100  int MC_kaonID;
101  int MC_michelID;
102 
103  double MC_leptonP;
104  double MC_leading_PionPlusP;
105  double MC_leading_ProtonP;
106  double MC_leading_PionMinusP;
107  double MC_kaonP;
108  double MC_michelP;
109 
110  TH1D *h_Ev_den;
111  TH1D *h_Ev_num;
112  TH1D *h_Pmu_den;
113  TH1D *h_Pmu_num;
114  TH1D *h_theta_den;
115  TH1D *h_theta_num;
116  TH1D *h_Pproton_den;
117  TH1D *h_Pproton_num;
118  TH1D *h_Ppion_plus_den;
119  TH1D *h_Ppion_plus_num;
120  TH1D *h_Ppion_minus_den;
121  TH1D *h_Ppion_minus_num;
122 
123  TH1D *h_Efrac_lepton;
124  TH1D *h_Ecomplet_lepton;
125  TH1D *h_Efrac_proton;
126  TH1D *h_Ecomplet_proton;
127  TH1D *h_Efrac_pion_plus;
128  TH1D *h_Ecomplet_pion_plus;
129  TH1D *h_Efrac_pion_minus;
130  TH1D *h_Ecomplet_pion_minus;
131  TH1D *h_trackRes_lepton;
132  TH1D *h_trackRes_proton;
133  TH1D *h_trackRes_pion_plus;
134  TH1D *h_trackRes_pion_minus;
135 
136  TH1D *h_muon_length;
137  TH1D *h_proton_length;
138  TH1D *h_pionp_length;
139  TH1D *h_pionm_length;
140  TH1D *h_muonwtrk_length;
141  TH1D *h_protonwtrk_length;
142  TH1D *h_pionpwtrk_length;
143  TH1D *h_pionmwtrk_length;
144 
145  TEfficiency* h_Eff_Ev = 0;
146  TEfficiency* h_Eff_Pmu = 0;
147  TEfficiency* h_Eff_theta = 0;
148  TEfficiency* h_Eff_Pproton = 0;
149  TEfficiency* h_Eff_Ppion_plus = 0;
150  TEfficiency* h_Eff_Ppion_minus = 0;
151 
152  TEfficiency* h_Eff_Lmuon = 0;
153  TEfficiency* h_Eff_Lproton = 0;
154  TEfficiency* h_Eff_Lpion_plus = 0;
155  TEfficiency* h_Eff_Lpion_minus = 0;
156 
157 
158  //nucleon decay histograms
159  TH1D *h_Pkaon_den;
160  TH1D *h_Pkaon_num;
161  TH1D *h_Pmichel_e_den;
162  TH1D *h_Pmichel_e_num;
163  TH1D *h_Efrac_kaon;
164  TH1D *h_Ecomplet_kaon;
165  TH1D *h_trackRes_kaon;
166  TH1D *h_Efrac_michel;
167  TH1D *h_Ecomplet_michel;
168  TH1D *h_trackRes_michel;
169  TH1D *h_kaon_length;
170  TH1D *h_michel_length;
171  TH1D *h_kaonwtrk_length;
172  TH1D *h_michelwtrk_length;
173  TEfficiency* h_Eff_Pkaon =0;
174  TEfficiency* h_Eff_Pmichel =0;
175  TEfficiency* h_Eff_Lkaon = 0;
176  TEfficiency* h_Eff_Lmichel =0;
177 
178 
179  float fFidVolCutX;
180  float fFidVolCutY;
181  float fFidVolCutZ;
182 
183  float fFidVolXmin;
184  float fFidVolXmax;
185  float fFidVolYmin;
186  float fFidVolYmax;
187  float fFidVolZmin;
188  float fFidVolZmax;
189 
190  detinfo::DetectorProperties const *detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
191  detinfo::DetectorClocks const *ts = lar::providerFrom<detinfo::DetectorClocksService>();
192  double XDriftVelocity = detprop->DriftVelocity()*1e-3; //cm/ns
193  double WindowSize = detprop->NumberTimeSamples() * ts->TPCClock().TickPeriod() * 1e3;
194  art::ServiceHandle<geo::Geometry> geom;
195 
196  }; // class NeutrinoTrackingEff
197 
198 
199 //========================================================================
200 NeutrinoTrackingEff::NeutrinoTrackingEff(fhicl::ParameterSet const& parameterSet)
201  : EDAnalyzer(parameterSet)
202 {
203  reconfigure(parameterSet);
204 }
205 //========================================================================
206 NeutrinoTrackingEff::~NeutrinoTrackingEff(){
207  //destructor
208 }
209 //========================================================================
210 void NeutrinoTrackingEff::reconfigure(fhicl::ParameterSet const& p){
211 
212  fMCTruthModuleLabel = p.get<std::string>("MCTruthModuleLabel");
213  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel");
214  fisNeutrinoInt = p.get<bool>("isNeutrinoInt");
215  fLeptonPDGcode = p.get<int>("LeptonPDGcode");
216  fNeutrinoPDGcode = p.get<int>("NeutrinoPDGcode");
217  fMaxNeutrinoE = p.get<double>("MaxNeutrinoE");
218  fMaxLeptonP = p.get<double>("MaxLeptonP");
219  fFidVolCutX = p.get<float>("FidVolCutX");
220  fFidVolCutY = p.get<float>("FidVolCutY");
221  fFidVolCutZ = p.get<float>("FidVolCutZ");
222 }
223 //========================================================================
224 void NeutrinoTrackingEff::beginJob(){
225  std::cout<<"job begin..."<<std::endl;
226  // Get geometry.
227  auto const* geo = lar::providerFrom<geo::Geometry>();
228  // Define histogram boundaries (cm).
229  // For now only draw cryostat=0.
230  double minx = 1e9;
231  double maxx = -1e9;
232  double miny = 1e9;
233  double maxy = -1e9;
234  double minz = 1e9;
235  double maxz = -1e9;
236  for (size_t i = 0; i<geo->NTPC(); ++i){
237  double local[3] = {0.,0.,0.};
238  double world[3] = {0.,0.,0.};
239  const geo::TPCGeo &tpc = geo->TPC(i);
240  tpc.LocalToWorld(local,world);
241  if (minx>world[0]-geo->DetHalfWidth(i))
242  minx = world[0]-geo->DetHalfWidth(i);
243  if (maxx<world[0]+geo->DetHalfWidth(i))
244  maxx = world[0]+geo->DetHalfWidth(i);
245  if (miny>world[1]-geo->DetHalfHeight(i))
246  miny = world[1]-geo->DetHalfHeight(i);
247  if (maxy<world[1]+geo->DetHalfHeight(i))
248  maxy = world[1]+geo->DetHalfHeight(i);
249  if (minz>world[2]-geo->DetLength(i)/2.)
250  minz = world[2]-geo->DetLength(i)/2.;
251  if (maxz<world[2]+geo->DetLength(i)/2.)
252  maxz = world[2]+geo->DetLength(i)/2.;
253  }
254 
255  fFidVolXmin = minx + fFidVolCutX;
256  fFidVolXmax = maxx - fFidVolCutX;
257  fFidVolYmin = miny + fFidVolCutY;
258  fFidVolYmax = maxy - fFidVolCutY;
259  fFidVolZmin = minz + fFidVolCutZ;
260  fFidVolZmax = maxz - fFidVolCutZ;
261 
262  std::cout<<"Fiducial volume:"<<"\n"
263  <<fFidVolXmin<<"\t< x <\t"<<fFidVolXmax<<"\n"
264  <<fFidVolYmin<<"\t< y <\t"<<fFidVolYmax<<"\n"
265  <<fFidVolZmin<<"\t< z <\t"<<fFidVolZmax<<"\n";
266 
267  art::ServiceHandle<art::TFileService> tfs;
268 
269  double E_bins[21] ={0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4,4.5,5.0,5.5,6.0,7.0,8.0,10.0,12.0,14.0,17.0,20.0,25.0};
270  double theta_bin[44]= { 0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,22.,24.,26.,28.,30.,32.,34.,36.,38.,40.,42.,44.,46.,48.,50.,55.,60.,65.,70.,75.,80.,85.,90.};
271  double Pbins[18] ={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.5,3.0};
272 
273  for (int i = 0; i<21; ++i) E_bins[i] *= fMaxNeutrinoE/25.;
274  for (int i = 0; i<18; ++i) Pbins[i] *= fMaxLeptonP/3.0;
275 
276  h_Ev_den = tfs->make<TH1D>("h_Ev_den","Neutrino Energy; Neutrino Energy (GeV); Tracking Efficiency",20,E_bins);
277  h_Ev_num = tfs->make<TH1D>("h_Ev_num","Neutrino Energy; Neutrino Energy (GeV); Tracking Efficiency",20,E_bins);
278  h_Pmu_den = tfs->make<TH1D>("h_Pmu_den","Muon Momentum; Muon Momentum (GeV); Tracking Efficiency",20,E_bins);
279  h_Pmu_num = tfs->make<TH1D>("h_Pmu_num","Muon Momentum; Muon Momentum (GeV); Tracking Efficiency",20,E_bins);
280  h_theta_den = tfs->make<TH1D>("h_theta_den","Theta; Theta w.r.t beam direction (Degrees); Tracking Efficiency",43,theta_bin);
281  h_theta_num = tfs->make<TH1D>("h_theta_num","Theta; Theta w.r.t beam direction (Degrees); Tracking Efficiency",43,theta_bin);
282  h_Pproton_den = tfs->make<TH1D>("h_Pproton_den","Protons; Proton Momentum (GeV); Tracking Efficiency", 17, Pbins);
283  h_Pproton_num = tfs->make<TH1D>("h_Pproton_num","Protons; Proton Momentum (GeV); Tracking Efficiency", 17, Pbins);
284  h_Ppion_plus_den = tfs->make<TH1D>("h_Ppion_plus_den", "Pions Plus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
285  h_Ppion_plus_num = tfs->make<TH1D>("h_Ppion_plus_num", "Pions Plus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
286  h_Ppion_minus_den = tfs->make<TH1D>("h_Ppion_minus_den", "Pions Minus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
287  h_Ppion_minus_num = tfs->make<TH1D>("h_Ppion_minus_num", "Pions Minus; Pion Momentum (GeV); Tracking Efficiency", 17, Pbins);
288  h_Ev_den->Sumw2();
289  h_Ev_num->Sumw2();
290  h_Pmu_den->Sumw2();
291  h_Pmu_num->Sumw2();
292  h_theta_den->Sumw2();
293  h_theta_num->Sumw2();
294  h_Pproton_den->Sumw2();
295  h_Pproton_num->Sumw2();
296  h_Ppion_plus_den->Sumw2();
297  h_Ppion_plus_num->Sumw2();
298  h_Ppion_minus_den->Sumw2();
299  h_Ppion_minus_num->Sumw2();
300 
301  h_Efrac_lepton = tfs->make<TH1D>("h_Efrac_lepton","Efrac Lepton; Track Purity;",60,0,1.2);
302  h_Ecomplet_lepton = tfs->make<TH1D>("h_Ecomplet_lepton","Ecomplet Lepton; Track Completeness;",60,0,1.2);
303  h_Efrac_proton = tfs->make<TH1D>("h_Efrac_proton","Efrac Proton; Track Purity;",60,0,1.2);
304  h_Ecomplet_proton = tfs->make<TH1D>("h_Ecomplet_proton","Ecomplet Proton; Track Completeness;",60,0,1.2);
305  h_Efrac_pion_plus = tfs->make<TH1D>("h_Efrac_pion_plus","Efrac Pion +; Track Purity;",60,0,1.2);
306  h_Ecomplet_pion_plus = tfs->make<TH1D>("h_Ecomplet_pion_plus","Ecomplet Pion +; Track Completeness;",60,0,1.2);
307  h_Efrac_pion_minus = tfs->make<TH1D>("h_Efrac_pion_minus","Efrac Pion -; Track Purity;",60,0,1.2);
308  h_Ecomplet_pion_minus = tfs->make<TH1D>("h_Ecomplet_pion_minus","Ecomplet Pion -; Track Completeness;",60,0,1.2);
309  h_trackRes_lepton = tfs->make<TH1D>("h_trackRes_lepton", "Muon Residual; Truth length - Reco length (cm);",200,-100,100);
310  h_trackRes_proton = tfs->make<TH1D>("h_trackRes_proton", "Proton Residual; Truth length - Reco length (cm);",200,-100,100);
311  h_trackRes_pion_plus = tfs->make<TH1D>("h_trackRes_pion_plus", "Pion + Residual; Truth length - Reco length (cm);",200,-100,100);
312  h_trackRes_pion_minus = tfs->make<TH1D>("h_trackRes_pion_minus", "Pion - Residual; Truth length - Reco length (cm);",200,-100,100);
313  h_Efrac_lepton->Sumw2();
314  h_Ecomplet_lepton->Sumw2();
315  h_Efrac_proton->Sumw2();
316  h_Ecomplet_proton->Sumw2();
317  h_Efrac_pion_plus->Sumw2();
318  h_Ecomplet_pion_plus->Sumw2();
319  h_Efrac_pion_minus->Sumw2();
320  h_Ecomplet_pion_minus->Sumw2();
321  h_trackRes_lepton->Sumw2();
322  h_trackRes_proton->Sumw2();
323  h_trackRes_pion_plus->Sumw2();
324  h_trackRes_pion_minus->Sumw2();
325 
326  h_muon_length = tfs->make<TH1D>("h_muon_length","Muon Length; Muon Truth Length (cm)",40,0,100);
327  h_proton_length = tfs->make<TH1D>("h_proton_length","Proton Length; Proton Truth Length (cm)",40,0,100);
328  h_pionp_length = tfs->make<TH1D>("h_pionp_length","Pion + Length; Pion^{+} Truth Length (cm)",40,0,100);
329  h_pionm_length = tfs->make<TH1D>("h_pionm_length","Pion - Length; Pion^{-} Truth Length (cm)",40,0,100);
330 
331  h_muonwtrk_length = tfs->make<TH1D>("h_muonwtrk_length","Muon Length; Muon Truth Length (cm)",40,0,100);
332  h_protonwtrk_length = tfs->make<TH1D>("h_protonwtrk_length","Proton Length; Proton Truth Length (cm)",40,0,100);
333  h_pionpwtrk_length = tfs->make<TH1D>("h_pionpwtrk_length","Pion + Length; Pion^{+} Truth Length (cm)",40,0,100);
334  h_pionmwtrk_length = tfs->make<TH1D>("h_pionmwtrk_length","Pion - Length; Pion^{-} Truth Length (cm)",40,0,100);
335 
336  h_muon_length->Sumw2();
337  h_muonwtrk_length->Sumw2();
338  h_proton_length->Sumw2();
339  h_protonwtrk_length->Sumw2();
340  h_pionp_length->Sumw2();
341  h_pionpwtrk_length->Sumw2();
342  h_pionm_length->Sumw2();
343  h_pionmwtrk_length->Sumw2();
344 
345  h_Pkaon_den = tfs->make<TH1D>("h_Pkaon_den","Kaon; Kaon Momentum (GeV); Tracking Efficiency", 17, Pbins);
346  h_Pkaon_num = tfs->make<TH1D>("h_Pkaon_num","Kaon; Kaon Momentum (GeV); Tracking Efficiency", 17, Pbins);
347  h_Pmichel_e_den = tfs->make<TH1D>("h_Pmichel_e_den","Michel Electron; Michele e Momentum (GeV); Tracking Efficiency", 17, Pbins);
348  h_Pmichel_e_num = tfs->make<TH1D>("h_Pmichel_e_num","Michel Electron; Michele e Momentum (GeV); Tracking Efficiency", 17, Pbins);
349  h_Pkaon_den->Sumw2();
350  h_Pkaon_num->Sumw2();
351  h_Pmichel_e_den->Sumw2();
352  h_Pmichel_e_num->Sumw2();
353  h_Efrac_kaon = tfs->make<TH1D>("h_Efrac_kaon","Efrac Kaon; Track Purity;",60,0,1.2);
354  h_Ecomplet_kaon = tfs->make<TH1D>("h_Ecomplet_kaon","Ecomplet Kaon; Track Completeness;",60,0,1.2);
355  h_trackRes_kaon = tfs->make<TH1D>("h_trackRes_kaon","Kaon Residual; Truth length - Reco length (cm);",200,-100,100);
356  h_Efrac_michel = tfs->make<TH1D>("h_Efrac_michel","Efrac Michel; Track Energy fraction;",60,0,1.2);
357  h_Ecomplet_michel = tfs->make<TH1D>("h_Ecomplet_michel","Ecomplet Michel; Track Completeness;",60,0,1.2);
358  h_trackRes_michel = tfs->make<TH1D>("h_trackRes_michel","Michel Residual; Truth length - Reco length (cm);",200,-100,100);
359  h_kaon_length = tfs->make<TH1D>("h_kaon_length","Kaon Length; Kaon Truth Length (cm)",40,0,100);
360  h_kaonwtrk_length = tfs->make<TH1D>("h_kaonwtrk_length","Kaon Length; Kaon Truth Length (cm)",40,0,100);
361  h_michel_length = tfs->make<TH1D>("h_michel_length","Michel Length; Michel e Truth Length (cm)",40,0,100);
362  h_michelwtrk_length = tfs->make<TH1D>("h_michelwtrk_length","Michel Length; Michel e Truth Length (cm)",40,0,100);
363 
364  h_Efrac_kaon->Sumw2();
365  h_Ecomplet_kaon->Sumw2();
366  h_trackRes_kaon->Sumw2();
367  h_Efrac_michel->Sumw2();
368  h_Ecomplet_michel->Sumw2();
369  h_trackRes_michel->Sumw2();
370  h_kaon_length->Sumw2();
371  h_kaonwtrk_length->Sumw2();
372  h_michel_length->Sumw2();
373  h_michelwtrk_length->Sumw2();
374 
375 }
376 //========================================================================
377 void NeutrinoTrackingEff::endJob(){
378 
379  doEfficiencies();
380 
381 }
382 //========================================================================
383 void NeutrinoTrackingEff::beginRun(const art::Run& /*run*/){
384  mf::LogInfo("NeutrinoTrackingEff")<<"begin run..."<<std::endl;
385 }
386 //========================================================================
387 void NeutrinoTrackingEff::analyze( const art::Event& event ){
388  if (event.isRealData()) return;
389 
390  bool isFiducial = false;
391  processEff(event, isFiducial);
392 }
393 //========================================================================
394 void NeutrinoTrackingEff::processEff( const art::Event& event, bool &isFiducial){
395 
397  art::Handle<std::vector<simb::MCTruth>> MCtruthHandle;
398  event.getByLabel(fMCTruthModuleLabel, MCtruthHandle);
399  std::vector<art::Ptr<simb::MCTruth>> MCtruthlist;
400  art::fill_ptr_vector(MCtruthlist, MCtruthHandle);
401  art::Ptr<simb::MCTruth> MCtruth;
402 
403  //For now assume that there is only one neutrino interaction...
404  MCtruth = MCtruthlist[0];
405  if( MCtruth->NeutrinoSet() ){
406  simb::MCNeutrino nu = MCtruth->GetNeutrino();
407  if( nu.CCNC() == 0 ) MC_isCC = 1;
408  else if ( nu.CCNC() == 1 ) MC_isCC = 0;
409  simb::MCParticle neutrino = nu.Nu();
410  MC_incoming_PDG = nu.Nu().PdgCode();
411  const TLorentzVector& nu_momentum = nu.Nu().Momentum(0);
412  nu_momentum.GetXYZT(MC_incoming_P);
413  const TLorentzVector& vertex =neutrino.Position(0);
414  vertex.GetXYZT(MC_vertex);
415  }
416 
418 
419  double tmp_leadingPionPlusE = 0.0;
420  double tmp_leadingPionMinusE = 0.0;
421  double tmp_leadingProtonE = 0.0;
422 
423  simb::MCParticle *MClepton = NULL;
424  simb::MCParticle *MCproton = NULL;
425  simb::MCParticle *MCpion_plus = NULL;
426  simb::MCParticle *MCpion_minus = NULL;
427  simb::MCParticle *MCkaon = NULL;
428  simb::MCParticle *MCmichel = NULL;
429 
430  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
431  const sim::ParticleList& plist = pi_serv->ParticleList();
432  simb::MCParticle *particle=0;
433  int i=0; // particle index
434 
435  for( sim::ParticleList::const_iterator ipar = plist.begin(); ipar!=plist.end(); ++ipar){
436  particle = ipar->second;
437  if( particle->PdgCode() == fLeptonPDGcode && particle->Mother() == 0 ){ //primary lepton
438  const TLorentzVector& lepton_momentum =particle->Momentum(0);
439  lepton_momentum.GetXYZT(MC_lepton_startMomentum);
440  MC_leptonID = particle->TrackId();
441  MC_leptonP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
442  MClepton = particle;
443  }
444  if( particle->Mother() == 0 ){ //save primary particle i.e. from the neutrino interaction
445  //save leading pion and proton
446  if( particle->PdgCode() == 2212 ){
447  if(particle->Momentum().E() > tmp_leadingProtonE){
448  tmp_leadingProtonE = particle->Momentum().E();
449  MC_leading_protonID = particle->TrackId();
450  MC_leading_ProtonP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
451  MCproton = particle;
452  }
453  }
454  else if( particle->PdgCode() == 211 ){
455  if(particle->Momentum().E() > tmp_leadingPionPlusE){
456  tmp_leadingPionPlusE = particle->Momentum().E();
457  MC_leading_PionPlusID = particle->TrackId();
458  MC_leading_PionPlusP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
459  MCpion_plus = particle;
460  }
461  }
462  else if( particle->PdgCode() == -211 ){
463  if(particle->Momentum().E() > tmp_leadingPionMinusE){
464  tmp_leadingPionMinusE = particle->Momentum().E();
465  MC_leading_PionMinusID = particle->TrackId();
466  MC_leading_PionMinusP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
467  MCpion_minus = particle;
468  }
469  }
470  i++; //paticle index
471  }
472 
473  //=======================================================================================
474  //add Nucleon decay stuff and particle cannon
475  //=======================================================================================
476  if(!fisNeutrinoInt ){
477  if( particle->Mother() == 0 ){
478  const TLorentzVector& positionStart = particle->Position(0);
479  positionStart.GetXYZT(MC_vertex);
480  }
481  if( particle->PdgCode() == 321 ){ //save primary Kaon
482  MC_kaonID = particle->TrackId();
483  MC_kaonP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
484  MCkaon = particle;
485  }
486  else if( particle->PdgCode() == fLeptonPDGcode ){ // Particle cannon muon
487  MC_leptonID = particle->TrackId();
488  MC_leptonP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
489  MClepton = particle;
490  }
491  else if( particle->PdgCode() == 2212 ){
492  if(particle->Momentum().E() > tmp_leadingProtonE){
493  tmp_leadingProtonE = particle->Momentum().E();
494  MC_leading_protonID = particle->TrackId();
495  MC_leading_ProtonP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
496  MCproton = particle;
497  }
498  }
499  else if( particle->PdgCode() == 211 ){
500  if(particle->Momentum().E() > tmp_leadingPionPlusE){
501  tmp_leadingPionPlusE = particle->Momentum().E();
502  MC_leading_PionPlusID = particle->TrackId();
503  MC_leading_PionPlusP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
504  MCpion_plus = particle;
505  }
506  }
507  else if( particle->PdgCode() == -211 ){
508  if(particle->Momentum().E() > tmp_leadingPionMinusE){
509  tmp_leadingPionMinusE = particle->Momentum().E();
510  MC_leading_PionMinusID = particle->TrackId();
511  MC_leading_PionMinusP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
512  MCpion_minus = particle;
513  }
514  }
515  else if( particle->Process() =="Decay" && particle->PdgCode() == -11){ // michel electron from muon decay
516  MC_michelID = particle->TrackId();
517  MC_michelP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
518  MCmichel = particle;
519  }
520  else if( TMath::Abs(particle->PdgCode() == 321) ){ //save primary Kaon
521  MC_kaonID = particle->TrackId();
522  MC_kaonP = sqrt(pow(particle->Momentum().Px(),2)+pow(particle->Momentum().Py(),2)+pow(particle->Momentum().Pz(),2));
523  MCkaon = particle;
524  }
525  }
526  }
527  //===================================================================
528  //Saving denominator histograms
529  //===================================================================
530  isFiducial =insideFV( MC_vertex );
531  if( !isFiducial ) return;
532  double Pv = sqrt(pow(MC_incoming_P[0],2)+pow(MC_incoming_P[1],2)+pow(MC_incoming_P[2],2));
533  double theta_mu = acos((MC_incoming_P[0]*MC_lepton_startMomentum[0] + MC_incoming_P[1]*MC_lepton_startMomentum[1] +MC_incoming_P[2]*MC_lepton_startMomentum[2])/(Pv*MC_leptonP) );
534  theta_mu *= (180.0/3.14159);
535  double truth_lengthLepton = truthLength(MClepton);
536  double proton_length = truthLength(MCproton);
537  double pion_plus_length = truthLength(MCpion_plus);
538  double pion_minus_length = truthLength(MCpion_minus);
539  double kaonLength = truthLength(MCkaon);
540  double michelLength = truthLength(MCmichel);
541 
542  //save CC events within the fiducial volume with the favorite neutrino flavor
543  if( MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE) ){
544  if( MClepton ){
545  h_Ev_den->Fill(MC_incoming_P[3]);
546  h_Pmu_den->Fill(MC_leptonP);
547  h_theta_den->Fill(theta_mu);
548  h_muon_length->Fill(truth_lengthLepton);
549  }
550  if( MCproton ){
551  h_Pproton_den->Fill(MC_leading_ProtonP);
552  h_proton_length->Fill(proton_length);
553  }
554  if( MCpion_plus ){
555  h_Ppion_plus_den->Fill( MC_leading_PionPlusP);
556  h_pionp_length->Fill(pion_plus_length);
557  }
558  if( MCpion_minus ){
559  h_Ppion_minus_den->Fill( MC_leading_PionMinusP);
560  h_pionm_length->Fill(pion_minus_length);
561  }
562  if( MCkaon ){
563  h_Pkaon_den->Fill(MC_kaonP);
564  h_kaon_length->Fill(kaonLength);
565  }
566  }
567 
568  //save events for Nucleon decay and particle cannon
569  if(!fisNeutrinoInt ){
570  if( MClepton ){
571  h_Pmu_den->Fill(MC_leptonP);
572  h_muon_length->Fill(truth_lengthLepton);
573  }
574  if( MCkaon ){
575  h_Pkaon_den->Fill(MC_kaonP);
576  h_kaon_length->Fill(kaonLength);
577  }
578  if( MCproton ){
579  h_Pproton_den->Fill(MC_leading_ProtonP);
580  h_proton_length->Fill(proton_length);
581  }
582  if( MCpion_plus ){
583  h_Ppion_plus_den->Fill( MC_leading_PionPlusP);
584  h_pionp_length->Fill(pion_plus_length);
585  }
586  if( MCpion_minus ){
587  h_Ppion_minus_den->Fill( MC_leading_PionMinusP);
588  h_pionm_length->Fill(pion_minus_length);
589  }
590  if( MCmichel ){
591  h_Pmichel_e_den->Fill(MC_michelP);
592  h_michel_length->Fill(michelLength);
593  }
594  }
595 
596  //========================================================================
597  // Reco stuff, once we have selected a MC Particle let's find out if there is a track associated
598  //========================================================================
599  art::Handle<std::vector<recob::Track>> trackListHandle;
600  if(! event.getByLabel(fTrackModuleLabel, trackListHandle)) return;
601  std::vector<art::Ptr<recob::Track> > tracklist;
602  art::fill_ptr_vector(tracklist, trackListHandle);
603  int n_recoTrack = tracklist.size();
604 
605  art::FindManyP<recob::Hit> track_hits(trackListHandle, event, fTrackModuleLabel);
606  if( n_recoTrack == 0 ){
607  LOG_DEBUG("NeutrinoTrackingEff")<<"There are no reco tracks... bye";
608  return;
609  }
610  LOG_DEBUG("NeutrinoTrackingEff")<<"Found this many reco tracks "<<n_recoTrack;
611 
612 
613  double Efrac_lepton =0.0;
614  double Ecomplet_lepton =0.0;
615  double Efrac_proton =0.0;
616  double Ecomplet_proton =0.0;
617  double Efrac_pionplus =0.0;
618  double Ecomplet_pionplus =0.0;
619  double Efrac_pionminus =0.0;
620  double Ecomplet_pionminus =0.0;
621  double Efrac_kaon =0.0;
622  double Ecomplet_kaon =0.0;
623  double Efrac_michel =0.0;
624  double Ecomplet_michel =0.0;
625  double trackLength_lepton =0.0;
626  double trackLength_proton =0.0;
627  double trackLength_pion_plus =0.0;
628  double trackLength_pion_minus =0.0;
629  double trackLength_kaon =0.0;
630  double trackLength_michel =0.0;
631  const simb::MCParticle *MClepton_reco = NULL;
632  const simb::MCParticle *MCproton_reco = NULL;
633  const simb::MCParticle *MCpion_plus_reco = NULL;
634  const simb::MCParticle *MCpion_minus_reco = NULL;
635  const simb::MCParticle *MCkaon_reco = NULL;
636  const simb::MCParticle *MCmichel_reco = NULL;
637 
638  std::vector<art::Ptr<recob::Hit>> tmp_all_trackHits = track_hits.at(0);
639  std::vector<art::Ptr<recob::Hit>> all_hits;
640  art::Handle<std::vector<recob::Hit>> hithandle;
641  if(event.get(tmp_all_trackHits[0].id(), hithandle)) art::fill_ptr_vector(all_hits, hithandle);
642 
643  for(int i=0; i<n_recoTrack; i++) {
644  art::Ptr<recob::Track> track = tracklist[i];
645  std::vector<art::Ptr<recob::Hit>> all_trackHits = track_hits.at(i);
646  double tmpEfrac = 0;
647  double tmpEcomplet =0;
648  const simb::MCParticle *particle;
649  truthMatcher( all_hits, all_trackHits, particle, tmpEfrac, tmpEcomplet );
650  if (!particle) continue;
651  //std::cout<<particle->PdgCode()<<" "<<particle->TrackId()<<" Efrac "<<tmpEfrac<<std::endl;
652  if( (particle->PdgCode() == fLeptonPDGcode) && (particle->TrackId() == MC_leptonID) ){
653  //save the best track ... based on completeness if there is more than one track
654  //if( tmpEfrac > Efrac_lepton ){ ///this was base on purity
655  if( tmpEcomplet > Ecomplet_lepton ){
656  Ecomplet_lepton = tmpEcomplet;
657  Efrac_lepton = tmpEfrac;
658  trackLength_lepton = track->Length();
659  MClepton_reco = particle;
660  }
661  }
662  else if( (particle->PdgCode() == 2212) && (particle->TrackId() == MC_leading_protonID) ){
663  //save the best track ... based on completeness if there is more than one track
664  if( tmpEcomplet > Ecomplet_proton ){
665  Ecomplet_proton = tmpEcomplet;
666  Efrac_proton = tmpEfrac;
667  trackLength_proton = track->Length();
668  MCproton_reco = particle;
669  }
670  }
671  else if( (particle->PdgCode() == 211) && (particle->TrackId() == MC_leading_PionPlusID) ){
672  //save the best track ... based on completeness if there is more than one track
673  if( tmpEcomplet > Ecomplet_pionplus ){
674  Ecomplet_pionplus = tmpEcomplet;
675  Efrac_pionplus = tmpEfrac;
676  trackLength_pion_plus = track->Length();
677  MCpion_plus_reco = particle;
678  }
679  }
680  else if( (particle->PdgCode() == -211) && (particle->TrackId() == MC_leading_PionMinusID) ){
681  //save the best track ... based on completeness if there is more than one track
682  if( tmpEcomplet > Ecomplet_pionminus ){
683  Ecomplet_pionminus = tmpEcomplet;
684  Efrac_pionminus = tmpEfrac;
685  trackLength_pion_minus = track->Length();
686  MCpion_minus_reco = particle;
687  }
688  }
689  //kaon from nucleon decay
690  else if( (TMath::Abs(particle->PdgCode()) == 321) && (particle->TrackId() == MC_kaonID) ){
691  //save the best track ... based on completeness if there is more than one track
692  if( tmpEcomplet > Ecomplet_kaon ){
693  Ecomplet_kaon = tmpEcomplet;
694  Efrac_kaon = tmpEfrac;
695  trackLength_kaon = track->Length();
696  MCkaon_reco = particle;
697  }
698  }
699  //michel from nucleon decay
700  else if( (particle->PdgCode() == -11) && (particle->TrackId() == MC_michelID) ){
701  //save the best track ... based on completeness if there is more than one track
702  if( tmpEcomplet > Ecomplet_michel ){
703  Ecomplet_michel = tmpEcomplet;
704  Efrac_michel = tmpEfrac;
705  trackLength_michel = track->Length();
706  MCmichel_reco = particle;
707  }
708  }
709 
710  }
711 
712  double Reco_LengthRes = truth_lengthLepton-trackLength_lepton;
713  double Reco_LengthResProton = proton_length-trackLength_proton;
714  double Reco_LengthResPionPlus = pion_plus_length-trackLength_pion_plus;
715  double Reco_LengthResPionMinus = pion_minus_length-trackLength_pion_minus;
716 
717  if( MClepton_reco && MClepton ){
718  if( MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE) ){
719  h_Pmu_num->Fill(MC_leptonP);
720  h_Ev_num->Fill(MC_incoming_P[3]);
721  h_theta_num->Fill(theta_mu);
722  h_Efrac_lepton->Fill(Efrac_lepton);
723  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
724  h_trackRes_lepton->Fill(Reco_LengthRes);
725  h_muonwtrk_length->Fill(truth_lengthLepton);
726  }
727  }
728  if( MCproton_reco && MCproton ){
729  if( MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE) ){
730  h_Pproton_num->Fill(MC_leading_ProtonP);
731  h_Efrac_proton->Fill(Efrac_proton);
732  h_Ecomplet_proton->Fill(Ecomplet_proton);
733  h_trackRes_proton->Fill(Reco_LengthResProton);
734  h_protonwtrk_length->Fill(proton_length);
735  }
736  }
737  if( MCpion_plus_reco && MCpion_plus ){
738  if( MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE) ){
739  h_Ppion_plus_num->Fill(MC_leading_PionPlusP);
740  h_Efrac_pion_plus->Fill(Efrac_pionplus);
741  h_Ecomplet_pion_plus->Fill(Ecomplet_pionplus);
742  h_trackRes_pion_plus->Fill(Reco_LengthResPionPlus);
743  h_pionpwtrk_length->Fill(pion_plus_length);
744  }
745  }
746  if( MCpion_minus_reco && MCpion_minus ){
747  if( MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE) ) {
748  h_Ppion_minus_num->Fill(MC_leading_PionMinusP);
749  h_Efrac_pion_minus->Fill(Efrac_pionminus);
750  h_Ecomplet_pion_minus->Fill(Ecomplet_pionminus);
751  h_trackRes_pion_minus->Fill(Reco_LengthResPionMinus);
752  h_pionmwtrk_length->Fill(pion_minus_length);
753  }
754  }
755  if( MCkaon_reco && MCkaon ){
756  if( MC_isCC && (fNeutrinoPDGcode == MC_incoming_PDG) && (MC_incoming_P[3] <= fMaxNeutrinoE) ) {
757  h_Pkaon_num->Fill(MC_kaonP);
758  h_Efrac_kaon->Fill(Efrac_kaon);
759  h_Ecomplet_kaon->Fill(Ecomplet_kaon);
760  h_trackRes_kaon->Fill(kaonLength-trackLength_kaon);
761  h_kaonwtrk_length->Fill(kaonLength);
762  }
763  }
764  //Non neutrino events
765  //=========================================================
766  if(!fisNeutrinoInt ){
767  if( MClepton_reco && MClepton ){
768  h_Pmu_num->Fill(MC_leptonP);
769  h_Efrac_lepton->Fill(Efrac_lepton);
770  h_Ecomplet_lepton->Fill(Ecomplet_lepton);
771  h_trackRes_lepton->Fill(Reco_LengthRes);
772  h_muonwtrk_length->Fill(truth_lengthLepton);
773  }
774  if( MCkaon_reco && MCkaon ){
775  h_Pkaon_num->Fill(MC_kaonP);
776  h_Efrac_kaon->Fill(Efrac_kaon);
777  h_Ecomplet_kaon->Fill(Ecomplet_kaon);
778  h_trackRes_kaon->Fill(kaonLength-trackLength_kaon);
779  h_kaonwtrk_length->Fill(kaonLength);
780  }
781  if( MCproton_reco && MCproton ){
782  h_Pproton_num->Fill(MC_leading_ProtonP);
783  h_Efrac_proton->Fill(Efrac_proton);
784  h_Ecomplet_proton->Fill(Ecomplet_proton);
785  h_trackRes_proton->Fill(Reco_LengthResProton);
786  h_protonwtrk_length->Fill(proton_length);
787  }
788  if( MCpion_plus_reco && MCpion_plus ){
789  h_Ppion_plus_num->Fill(MC_leading_PionPlusP);
790  h_Efrac_pion_plus->Fill(Efrac_pionplus);
791  h_Ecomplet_pion_plus->Fill(Ecomplet_pionplus);
792  h_trackRes_pion_plus->Fill(Reco_LengthResPionPlus);
793  h_pionpwtrk_length->Fill(pion_plus_length);
794  }
795  if( MCpion_minus_reco && MCpion_minus ){
796  h_Ppion_minus_num->Fill(MC_leading_PionMinusP);
797  h_Efrac_pion_minus->Fill(Efrac_pionminus);
798  h_Ecomplet_pion_minus->Fill(Ecomplet_pionminus);
799  h_trackRes_pion_minus->Fill(Reco_LengthResPionMinus);
800  h_pionmwtrk_length->Fill(pion_minus_length);
801  }
802  if( MCmichel_reco && MCmichel ){
803  h_Pmichel_e_num->Fill(MC_michelP);
804  h_Efrac_michel->Fill(Efrac_michel);
805  h_Ecomplet_michel->Fill(Ecomplet_michel);
806  h_trackRes_michel->Fill(michelLength-trackLength_michel);
807  h_michelwtrk_length->Fill(michelLength);
808  }
809 
810  }
811 
812 }
813 //========================================================================
814 void NeutrinoTrackingEff::truthMatcher( std::vector<art::Ptr<recob::Hit>> all_hits, std::vector<art::Ptr<recob::Hit>> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet){
815 
816  //std::cout<<"truthMatcher..."<<std::endl;
817  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
818  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
819  std::map<int,double> trkID_E;
820  for(size_t j = 0; j < track_hits.size(); ++j){
821  art::Ptr<recob::Hit> hit = track_hits[j];
822  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(hit);
823  for(size_t k = 0; k < TrackIDs.size(); k++){
824  trkID_E[TrackIDs[k].trackID] += TrackIDs[k].energy;
825  }
826  }
827  double E_em =0.0;
828  double max_E = -999.0;
829  double total_E = 0.0;
830  int TrackID = -999;
831  double partial_E =0.0; // amount of energy deposited by the particle that deposited more energy... tomato potato... blabla
834  if( !trkID_E.size() ) {
835  MCparticle = 0;
836  return; //Ghost track???
837  }
838  for(std::map<int,double>::iterator ii = trkID_E.begin(); ii!=trkID_E.end(); ++ii){
839  total_E += ii->second;
840  if((ii->second)>max_E){
841  partial_E = ii->second;
842  max_E = ii->second;
843  TrackID = ii->first;
844  if( TrackID < 0 ) E_em += ii->second;
845  }
846  }
847  MCparticle = pi_serv->TrackIdToParticle_P(TrackID);
848 
849  //In the current simulation, we do not save EM Shower daughters in GEANT. But we do save the energy deposition in TrackIDEs. If the energy deposition is from a particle that is the daughter of
850  //an EM particle, the negative of the parent track ID is saved in TrackIDE for the daughter particle
851  //we don't want to track gammas or any other EM activity
852  if( TrackID < 0 ) return;
853 
854  //Efrac = (partial_E+E_em)/total_E;
855  Efrac = (partial_E)/total_E;
856 
857  //completeness
858  double totenergy =0;
859  for(size_t k = 0; k < all_hits.size(); ++k){
860  art::Ptr<recob::Hit> hit = all_hits[k];
861  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(hit);
862  for(size_t l = 0; l < TrackIDs.size(); ++l){
863  if(TrackIDs[l].trackID==TrackID) totenergy += TrackIDs[l].energy;
864  }
865  }
866  Ecomplet = partial_E/totenergy;
867 }
868 //========================================================================
869 double NeutrinoTrackingEff::truthLength( const simb::MCParticle *MCparticle ){
870  //calculate the truth length considering only the part that is inside the TPC
871  //Base on a peace of code from dune/TrackingAna/TrackingEfficiency_module.cc
872 
873  if( !MCparticle ) return -999.0;
874  int numberTrajectoryPoints = MCparticle->NumberTrajectoryPoints();
875  std::vector<double> TPCLengthHits(numberTrajectoryPoints, 0.0);
876  int FirstHit=0, LastHit=0;
877  double TPCLength = 0.0;
878  bool BeenInVolume = false;
879 
880  for(unsigned int MCHit=0; MCHit < TPCLengthHits.size(); ++MCHit) {
881  const TLorentzVector& tmpPosition= MCparticle->Position(MCHit);
882  double const tmpPosArray[]={tmpPosition[0],tmpPosition[1],tmpPosition[2]};
883  if (MCHit!=0) TPCLengthHits[MCHit] = sqrt( pow( (MCparticle->Vx(MCHit-1)-MCparticle->Vx(MCHit)),2)+ pow( (MCparticle->Vy(MCHit-1)-MCparticle->Vy(MCHit)),2)+ pow( (MCparticle->Vz(MCHit-1)-MCparticle->Vz(MCHit)),2));
884  geo::TPCID tpcid = geom->FindTPCAtPosition(tmpPosArray);
885  if(tpcid.isValid) {
886  // -- Check if hit is within drift window...
887  geo::CryostatGeo const& cryo = geom->Cryostat(tpcid.Cryostat);
888  geo::TPCGeo const& tpc = cryo.TPC(tpcid.TPC);
889  double XPlanePosition = tpc.PlaneLocation(0)[0];
890  double DriftTimeCorrection = fabs( tmpPosition[0] - XPlanePosition ) / XDriftVelocity;
891  double TimeAtPlane = MCparticle->T() + DriftTimeCorrection;
892  if( TimeAtPlane < detprop->TriggerOffset() || TimeAtPlane > detprop->TriggerOffset() + WindowSize ) continue;
893  LastHit = MCHit;
894  if( !BeenInVolume ) {
895  BeenInVolume = true;
896  FirstHit = MCHit;
897  }
898  }
899  }
900  for (int Hit = FirstHit+1; Hit <= LastHit; ++Hit ) TPCLength += TPCLengthHits[Hit];
901  return TPCLength;
902 }
903 //========================================================================
904 bool NeutrinoTrackingEff::insideFV( double vertex[4]){
905 
906  double x = vertex[0];
907  double y = vertex[1];
908  double z = vertex[2];
909 
910  if (x>fFidVolXmin && x<fFidVolXmax&&
911  y>fFidVolYmin && y<fFidVolYmax&&
912  z>fFidVolZmin && z<fFidVolZmax)
913  return true;
914  else
915  return false;
916 }
917 //========================================================================
918 void NeutrinoTrackingEff::doEfficiencies(){
919 
920  art::ServiceHandle<art::TFileService> tfs;
921 
922  if(TEfficiency::CheckConsistency(*h_Ev_num,*h_Ev_den)){
923  h_Eff_Ev = tfs->make<TEfficiency>(*h_Ev_num,*h_Ev_den);
924  TGraphAsymmErrors *grEff_Ev = h_Eff_Ev->CreateGraph();
925  grEff_Ev->Write("grEff_Ev");
926  h_Eff_Ev->Write("h_Eff_Ev");
927  }
928  if(TEfficiency::CheckConsistency(*h_Pmu_num,*h_Pmu_den)){
929  h_Eff_Pmu = tfs->make<TEfficiency>(*h_Pmu_num,*h_Pmu_den);
930  TGraphAsymmErrors *grEff_Pmu = h_Eff_Pmu->CreateGraph();
931  grEff_Pmu->Write("grEff_Pmu");
932  h_Eff_Pmu->Write("h_Eff_Pmu");
933  }
934  if(TEfficiency::CheckConsistency(*h_theta_num,*h_theta_den)){
935  h_Eff_theta = tfs->make<TEfficiency>(*h_theta_num,*h_theta_den);
936  TGraphAsymmErrors *grEff_theta = h_Eff_theta->CreateGraph();
937  grEff_theta->Write("grEff_theta");
938  h_Eff_theta->Write("h_Eff_theta");
939  }
940  if(TEfficiency::CheckConsistency(*h_Pproton_num,*h_Pproton_den)){
941  h_Eff_Pproton = tfs->make<TEfficiency>(*h_Pproton_num,*h_Pproton_den);
942  TGraphAsymmErrors *grEff_Pproton = h_Eff_Pproton->CreateGraph();
943  grEff_Pproton->Write("grEff_Pproton");
944  h_Eff_Pproton->Write("h_Eff_Pproton");
945  }
946  if(TEfficiency::CheckConsistency(*h_Ppion_plus_num,*h_Ppion_plus_den)){
947  h_Eff_Ppion_plus = tfs->make<TEfficiency>(*h_Ppion_plus_num,*h_Ppion_plus_den);
948  TGraphAsymmErrors *grEff_Ppion_plus = h_Eff_Ppion_plus->CreateGraph();
949  grEff_Ppion_plus->Write("grEff_Ppion_plus");
950  h_Eff_Ppion_plus->Write("h_Eff_Ppion_plus");
951  }
952  if(TEfficiency::CheckConsistency(*h_Ppion_minus_num,*h_Ppion_minus_den)){
953  h_Eff_Ppion_minus = tfs->make<TEfficiency>(*h_Ppion_minus_num,*h_Ppion_minus_den);
954  TGraphAsymmErrors *grEff_Ppion_minus = h_Eff_Ppion_minus->CreateGraph();
955  grEff_Ppion_minus->Write("grEff_Ppion_minus");
956  h_Eff_Ppion_minus->Write("h_Eff_Ppion_minus");
957  }
958  if(TEfficiency::CheckConsistency(*h_muonwtrk_length,*h_muon_length)){
959  h_Eff_Lmuon = tfs->make<TEfficiency>(*h_muonwtrk_length,*h_muon_length);
960  TGraphAsymmErrors *grEff_Lmuon = h_Eff_Lmuon->CreateGraph();
961  grEff_Lmuon->Write("grEff_Lmuon");
962  h_Eff_Lmuon->Write("h_Eff_Lmuon");
963  }
964  if(TEfficiency::CheckConsistency(*h_protonwtrk_length,*h_proton_length)){
965  h_Eff_Lproton = tfs->make<TEfficiency>(*h_protonwtrk_length,*h_proton_length);
966  TGraphAsymmErrors *grEff_Lproton = h_Eff_Lproton->CreateGraph();
967  grEff_Lproton->Write("grEff_Lproton");
968  h_Eff_Lproton->Write("h_Eff_Lproton");
969  }
970  if(TEfficiency::CheckConsistency(*h_pionpwtrk_length,*h_pionp_length)){
971  h_Eff_Lpion_plus = tfs->make<TEfficiency>(*h_pionpwtrk_length,*h_pionp_length);
972  TGraphAsymmErrors *grEff_Lpion_plus = h_Eff_Lpion_plus->CreateGraph();
973  grEff_Lpion_plus->Write("grEff_Lpion_plus");
974  h_Eff_Lpion_plus->Write("h_Eff_Lpion_plus");
975  }
976  if(TEfficiency::CheckConsistency(*h_pionpwtrk_length,*h_pionp_length)){
977  h_Eff_Lpion_minus = tfs->make<TEfficiency>(*h_pionmwtrk_length,*h_pionm_length);
978  TGraphAsymmErrors *grEff_Lpion_minus = h_Eff_Lpion_minus->CreateGraph();
979  grEff_Lpion_minus->Write("grEff_Lpion_minus");
980  h_Eff_Lpion_minus->Write("h_Eff_Lpion_minus");
981  }
982  if(TEfficiency::CheckConsistency(*h_Pkaon_num,*h_Pkaon_den)){
983  h_Eff_Pkaon = tfs->make<TEfficiency>(*h_Pkaon_num,*h_Pkaon_den);
984  TGraphAsymmErrors *grEff_Pkaon = h_Eff_Pkaon->CreateGraph();
985  grEff_Pkaon->Write("grEff_Pkaon");
986  h_Eff_Pkaon->Write("h_Eff_Pkaon");
987  }
988  if(TEfficiency::CheckConsistency(*h_kaonwtrk_length,*h_kaon_length)){
989  h_Eff_Lkaon = tfs->make<TEfficiency>(*h_kaonwtrk_length,*h_kaon_length);
990  TGraphAsymmErrors *grEff_Lkaon = h_Eff_Lkaon->CreateGraph();
991  grEff_Lkaon->Write("grEff_Lkaon");
992  h_Eff_Lkaon->Write("h_Eff_Lkaon");
993  }
994  if(TEfficiency::CheckConsistency(*h_Pmichel_e_num,*h_Pmichel_e_den)){
995  h_Eff_Pmichel = tfs->make<TEfficiency>(*h_Pmichel_e_num,*h_Pmichel_e_den);
996  TGraphAsymmErrors *grEff_Pmichel = h_Eff_Pmichel->CreateGraph();
997  grEff_Pmichel->Write("grEff_Pmichel");
998  h_Eff_Pmichel->Write("h_Eff_Pmichel");
999  }
1000  if(TEfficiency::CheckConsistency(*h_michelwtrk_length,*h_michel_length)){
1001  h_Eff_Lmichel = tfs->make<TEfficiency>(*h_michelwtrk_length,*h_michel_length);
1002  TGraphAsymmErrors *grEff_Lmichel = h_Eff_Lmichel->CreateGraph();
1003  grEff_Lmichel->Write("grEff_Lmichel");
1004  h_Eff_Lmichel->Write("h_Eff_Lmichel");
1005  }
1006 
1007 }
1008 //========================================================================
1009 DEFINE_ART_MODULE(NeutrinoTrackingEff)
1010 
1011 }
1012 
1013 #endif // NeutrinoTrackingEff_Module
DUNE::NeutrinoTrackingEff::truthLength
double truthLength(const simb::MCParticle *MCparticle)
Definition: NeutrinoTrackingEff_module.cc:869
x
Float_t x
Definition: compare.C:6
LArG4Parameters.h
Store parameters for running LArG4.
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int const &id)
Definition: ParticleInventoryService_service.cc:134
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:222
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:223
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:216
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:152
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:74
detinfo::DetectorClocks::TPCClock
virtual const ::detinfo::ElecClock & TPCClock() const =0
Lends a constant TPC clock with time set to trigger time.
detinfo::DetectorProperties::TriggerOffset
virtual int TriggerOffset() const =0
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
cheat::BackTrackerService::HitToTrackIDEs
const std::vector< sim::TrackIDE > HitToTrackIDEs(recob::Hit const &hit)
Definition: BackTrackerService_service.cc:145
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:212
TrackIDs
std::vector< TrackID > TrackIDs
Definition: ClusteringValidation_module.cc:75
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:217
y
Float_t y
Definition: compare.C:6
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:150
sim::ParticleList::const_iterator
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
z
Double_t z
Definition: plot.C:279
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:37
art::DataViewImpl::get
bool get(SelectorBase const &, Handle< PROD > &result) const
Definition: DataViewImpl.h:307
std
STL namespace.
art::Event::isRealData
bool isRealData() const
Definition: Event.h:83
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:36
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:219
MCParticle.h
Particle class.
art::Run
Definition: Run.h:30
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:214
Event.h
DUNE::NeutrinoTrackingEff::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Definition: NeutrinoTrackingEff_module.cc:210
geo::GeometryCore::FindTPCAtPosition
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.h:2417
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
DUNE::NeutrinoTrackingEff::truthMatcher
void truthMatcher(std::vector< art::Ptr< recob::Hit >> all_hits, std::vector< art::Ptr< recob::Hit >> track_hits, const simb::MCParticle *&MCparticle, double &Efrac, double &Ecomplet)
Definition: NeutrinoTrackingEff_module.cc:814
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:174
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
breakpoints::beginJob
void beginJob()
Definition: Breakpoints.cc:14
DUNE::NeutrinoTrackingEff::analyze
void analyze(const art::Event &evt)
Definition: NeutrinoTrackingEff_module.cc:387
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:228
sim::ParticleList::begin
iterator begin()
Definition: ParticleList.h:305
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:280
DUNE::NeutrinoTrackingEff::processEff
void processEff(const art::Event &evt, bool &isFiducial)
Definition: NeutrinoTrackingEff_module.cc:394
TrackID
TrackID
Definition: ClusteringValidation_module.cc:73
detinfo::DetectorProperties::NumberTimeSamples
virtual unsigned int NumberTimeSamples() const =0
DUNE::NeutrinoTrackingEff::beginRun
void beginRun(const art::Run &run)
Definition: NeutrinoTrackingEff_module.cc:383
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:195
Track.h
Provides recob::Track data product.
geo_types.h
Definition of data types for geometry description.
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
detinfo::DetectorClocks
Conversion of times between different formats and references.
Definition: DetectorClocks.h:298
DUNE::NeutrinoTrackingEff::geom
art::ServiceHandle< geo::Geometry > geom
Definition: NeutrinoTrackingEff_module.cc:194
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:225
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
detinfo::ElecClock::TickPeriod
double TickPeriod() const
A single tick period in microseconds.
Definition: ElecClock.h:342
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:174
detinfo::DetectorProperties::DriftVelocity
virtual double DriftVelocity(double efield=0., double temperature=0.) const =0
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:224
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:227
LOG_DEBUG
#define LOG_DEBUG(id)
Definition: MessageLogger.h:263
simb::MCTruth::NeutrinoSet
bool NeutrinoSet() const
Definition: MCTruth.h:75
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
e
Float_t e
Definition: plot.C:34
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
geo
Namespace collecting geometry-related classes utilities.
Definition: TrackContainmentAlg.hh:20
track
Float_t track
Definition: plot.C:34
geo::TPCGeo::LocalToWorld
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:490
geo::TPCGeo::PlaneLocation
const double * PlaneLocation(unsigned int p) const
Returns the coordinates of the center of the specified plane [cm].
Definition: TPCGeo.cxx:412
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:226
event
Event finding and building.
Definition: AggregateEvent_module.cc:55
vertex
vertex reconstruction
Definition: AggregateVertex_module.cc:43