MVAAlg.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////
2 // \file MVAAlg.cxx
3 // tylerdalion@gmail.com
4 ////////////////////////////////////////////////////////////////////
5 
6 #include "dune/FDSensOpt/MVAAlg/MVAAlg.h"
7 #include "larcore/Geometry/Geometry.h"
8 #include "larcorealg/Geometry/TPCGeo.h"
9 #include "lardataalg/DetectorInfo/DetectorProperties.h"
10 #include "lardata/ArtDataHelper/TrackUtils.h" // lar::utils::TrackPitchInView()
11 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
12 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
13 #include "art/Framework/Services/Registry/ServiceHandle.h"
14 #include "TPrincipal.h"
15 #include "TVectorD.h"
16 #include "TF1.h"
17 #include "canvas/Persistency/Common/FindManyP.h"
18 #include "canvas/Persistency/Common/FindOneP.h"
19 #include "messagefacility/MessageLogger/MessageLogger.h"
20 #include "lardata/Utilities/AssociationUtil.h"
21 
22 //--------------------------------------------------------------------------------
23 dunemva::MVAAlg::MVAAlg( fhicl::ParameterSet const& p )
24  : fReader("")
25  , fCalorimetryAlg (p.get<fhicl::ParameterSet>("CalorimetryAlg"))
26  , fMakeAnaTree (p.get<bool>("MakeAnaTree"))
27  , fMakeWeightTree (p.get<bool>("MakeWeightTree"))
28 {
29 
30  this->reconfigure(p);
31  if(fMakeAnaTree) this->MakeTree();
32 
33  fFile.open("events.txt"); // for displaying selected events
34 
35  if(!fMakeWeightTree){
36 
37  // branch, name, unit, type
38  fReader.AddVariable("evtcharge", &evtcharge);
39  fReader.AddVariable("ntrack", &ntrack);
40  fReader.AddVariable("maxtrklength", &maxtrklength);
41  fReader.AddVariable("avgtrklength", &avgtrklength);
42  fReader.AddVariable("trkdedx", &trkdedx);
43  fReader.AddVariable("trkrch", &trkrch);
44  fReader.AddVariable("trkrt", &trkrt);
45  fReader.AddVariable("trkfr", &trkfr);
46  fReader.AddVariable("trkpida", &trkpida_save);
47  fReader.AddVariable("fract_5_wires", &fract_5_wires);
48  fReader.AddVariable("fract_10_wires", &fract_10_wires);
49  fReader.AddVariable("fract_50_wires", &fract_50_wires);
50  fReader.AddVariable("fract_100_wires", &fract_100_wires);
51  fReader.AddVariable("trkcosx", &trkcosx);
52  fReader.AddVariable("trkcosy", &trkcosy);
53  fReader.AddVariable("trkcosz", &trkcosz);
54  fReader.AddVariable("ET", &ET);
55 
56  // Nice to plot and verify sig/back sample composition
57  //fReader.AddVariable("ccnc", &ccnc);
58  //fReader.AddVariable("NuPdg", &NuPdg); // must be floats
59 
60  if(fSelect=="nue"){
61  fReader.AddVariable("nshower", &nshower);
62  fReader.AddVariable("showerdedx", &showerdedx);
63  fReader.AddVariable("eshower", &eshower);
64  fReader.AddVariable("frshower", &frshower);
65  fReader.AddVariable("nhitspershw", &nhitspershw);
66  fReader.AddVariable("shwlength", &shwlength);
67  fReader.AddVariable("shwmax", &shwmax);
68  fReader.AddVariable("shwdisx", &shwdisx);
69  fReader.AddVariable("shwdisy", &shwdisy);
70  fReader.AddVariable("shwdisz", &shwdisz);
71  fReader.AddVariable("shwcosx", &shwcosx);
72  fReader.AddVariable("shwcosy", &shwcosy);
73  fReader.AddVariable("shwcosz", &shwcosz);
74  }
75 
76 
77  /*
78  fReader.AddVariable("evtcharge", &evtcharge);
79  fReader.AddVariable("ntrack", &ntrack);
80  fReader.AddVariable("maxtrklength", &maxtrklength);
81  fReader.AddVariable("trkdedx", &trkdedx);
82  fReader.AddVariable("trkrch", &trkrch);
83  fReader.AddVariable("trkrt", &trkrt);
84  fReader.AddVariable("trkfr", &trkfr);
85  fReader.AddVariable("trkpida", &trkpida_save);
86  fReader.AddVariable("fract_5_wires", &fract_5_wires);
87  fReader.AddVariable("fract_10_wires", &fract_10_wires);
88  fReader.AddVariable("fract_50_wires", &fract_50_wires);
89  fReader.AddVariable("fract_100_wires", &fract_100_wires);
90  fReader.AddVariable("trkcosx", &trkcosx);
91  fReader.AddVariable("trkcosy", &trkcosy);
92  fReader.AddVariable("trkcosz", &trkcosz);
93  fReader.AddVariable("ET", &ET);
94  fReader.AddVariable("nshower", &nshower);
95  fReader.AddVariable("showerdedx", &showerdedx);
96  fReader.AddVariable("eshower", &eshower);
97  fReader.AddVariable("frshower", &frshower);
98  fReader.AddVariable("nhitspershw", &nhitspershw);
99  fReader.AddVariable("shwlength", &shwlength);
100  fReader.AddVariable("shwmax", &shwmax);
101  fReader.AddVariable("shwdisx", &shwdisx);
102  fReader.AddVariable("shwdisy", &shwdisy);
103  fReader.AddVariable("shwdisz", &shwdisz);
104  fReader.AddVariable("shwcosx", &shwcosx);
105  fReader.AddVariable("shwcosy", &shwcosy);
106  fReader.AddVariable("shwcosz", &shwcosz);
107  */
108 
109 
110  fMVAMethod =p.get< std::string >("MVAMethod");
111 
112  std::string weightFileFull;
113  fWeightFile=p.get< std::string >("WeightFile");
114  cet::search_path sp("FW_SEARCH_PATH");
115  sp.find_file(fWeightFile, weightFileFull);
116  fReader.BookMVA(fMVAMethod, weightFileFull);
117 
118  }
119 
120  fWeightTree[0] = tfs->make<TTree>("all","a Tree for MVA");
121  fWeightTree[1] = tfs->make<TTree>("sig","a Tree for MVA");
122  fWeightTree[2] = tfs->make<TTree>("nc","a Tree for MVA");
123  fWeightTree[3] = tfs->make<TTree>("cc","a Tree for MVA");
124  for (int i = 0; i<4; ++i){
125  fWeightTree[i]->Branch("run",&run,"run/I");
126  fWeightTree[i]->Branch("subrun",&subrun,"subrun/I");
127  fWeightTree[i]->Branch("event",&event,"event/I");
128  fWeightTree[i]->Branch("itype",&itype,"itype/I");
129 
130  // all tmva reader variables must be floats
131  fWeightTree[i]->Branch("oscpro",&fOscPro,"oscpro/F");
132  fWeightTree[i]->Branch("weight",&weight,"weight/F");
133  fWeightTree[i]->Branch("evtcharge",&evtcharge,"evtcharge/F");
134  fWeightTree[i]->Branch("ntrack",&ntrack,"ntrack/F");
135  fWeightTree[i]->Branch("maxtrklength",&maxtrklength,"maxtrklength/F");
136  fWeightTree[i]->Branch("avgtrklength",&avgtrklength,"avgtrklength/F");
137  fWeightTree[i]->Branch("trkdedx",&trkdedx,"trkdedx/F");
138  fWeightTree[i]->Branch("trkrch",&trkrch,"trkrch/F");
139  fWeightTree[i]->Branch("trkrt",&trkrt,"trkrt/F");
140  fWeightTree[i]->Branch("trkfr",&trkfr,"trkfr/F");
141  fWeightTree[i]->Branch("trkpida",&trkpida_save,"trkpida/F");
142  fWeightTree[i]->Branch("nshower",&nshower,"nshower/F");
143  fWeightTree[i]->Branch("showerdedx",&showerdedx,"showerdedx/F");
144  fWeightTree[i]->Branch("eshower",&eshower,"eshower/F");
145  fWeightTree[i]->Branch("frshower",&frshower,"frshower/F");
146  fWeightTree[i]->Branch("nhitspershw",&nhitspershw,"nhitspershw/F");
147  fWeightTree[i]->Branch("shwlength",&shwlength,"shwlength/F");
148  fWeightTree[i]->Branch("shwmax",&shwmax,"shwmax/F");
149  fWeightTree[i]->Branch("fract_5_wires",&fract_5_wires,"fract_5_wires/F");
150  fWeightTree[i]->Branch("fract_10_wires",&fract_10_wires,"fract_10_wires/F");
151  fWeightTree[i]->Branch("fract_50_wires",&fract_50_wires,"fract_50_wires/F");
152  fWeightTree[i]->Branch("fract_100_wires",&fract_100_wires,"fract_100_wires/F");
153  fWeightTree[i]->Branch("shwdis",&shwdis,"shwdis/F");
154  fWeightTree[i]->Branch("shwdisx",&shwdisx,"shwdisx/F");
155  fWeightTree[i]->Branch("shwdisy",&shwdisy,"shwdisy/F");
156  fWeightTree[i]->Branch("shwdisz",&shwdisz,"shwdisz/F");
157  fWeightTree[i]->Branch("shwcosx",&shwcosx,"shwcosx/F");
158  fWeightTree[i]->Branch("shwcosy",&shwcosy,"shwcosy/F");
159  fWeightTree[i]->Branch("shwcosz",&shwcosz,"shwcosz/F");
160  fWeightTree[i]->Branch("trkcosx",&trkcosx,"trkcosx/F");
161  fWeightTree[i]->Branch("trkcosy",&trkcosy,"trkcosy/F");
162  fWeightTree[i]->Branch("trkcosz",&trkcosz,"trkcosz/F");
163  fWeightTree[i]->Branch("ET",&ET,"ET/F");
164 
165  // Nice to double check samples through to the end
166  //fWeightTree[i]->Branch("ccnc",&ccnc,"ccnc/F");
167  //fWeightTree[i]->Branch("NuPdg",&NuPdg,"NuPdg/F");
168  }
169 
170  int mva_bins = 80;
171  mva_nue_osc = tfs->make<TH1D>("mva_nue_osc","mva_nue_osc",mva_bins,-1,1);
172  mva_nc = tfs->make<TH1D>("mva_nc","mva_nc",mva_bins,-1,1);
173  mva_numu = tfs->make<TH1D>("mva_numu","mva_numu",mva_bins,-1,1);
174  mva_nue_beam = tfs->make<TH1D>("mva_nue_beam","mva_nue_beam",mva_bins,-1,1);
175  mva_nutau = tfs->make<TH1D>("mva_nutau","mva_nutau",mva_bins,-1,1);
176 
177  mva_numu_nue = tfs->make<TH1D>("mva_numu_nue","mva_numu_nue",mva_bins,-1,1);
178  mva_nue_nue = tfs->make<TH1D>("mva_nue_nue","mva_nue_nue",mva_bins,-1,1);
179  mva_numu_numu = tfs->make<TH1D>("mva_numu_numu","mva_numu_numu",mva_bins,-1,1);
180  mva_nue_numu = tfs->make<TH1D>("mva_nue_numu","mva_nue_numu",mva_bins,-1,1);
181  mva_numu_nutau = tfs->make<TH1D>("mva_numu_nutau","mva_numu_nutau",mva_bins,-1,1);
182  mva_nue_nutau = tfs->make<TH1D>("mva_nue_nutau","mva_nue_nutau",mva_bins,-1,1);
183 
184  enu_nc = tfs->make<TH1D>("enu_nc","enu_nc",100,0,100);
185  enu_numu_nue = tfs->make<TH1D>("enu_numu_nue","enu_numu_nue",100,0,100);
186  enu_nue_nue = tfs->make<TH1D>("enu_nue_nue","enu_nue_nue",100,0,100);
187  enu_numu_numu = tfs->make<TH1D>("enu_numu_numu","enu_numu_numu",100,0,100);
188  enu_nue_numu = tfs->make<TH1D>("enu_nue_numu","enu_nue_numu",100,0,100);
189  enu_numu_nutau = tfs->make<TH1D>("enu_numu_nutau","enu_numu_nutau",100,0,100);
190  enu_nue_nutau = tfs->make<TH1D>("enu_nue_nutau","enu_nue_nutau",100,0,100);
191 
192  mva_nue_osc->Sumw2();
193  mva_nc->Sumw2();
194  mva_numu->Sumw2();
195  mva_nue_beam->Sumw2();
196  mva_nutau->Sumw2();
197 
198  mva_numu_nue->Sumw2();
199  mva_nue_nue->Sumw2();
200  mva_numu_numu->Sumw2();
201  mva_nue_numu->Sumw2();
202  mva_numu_nutau->Sumw2();
203  mva_nue_nutau->Sumw2();
204 
205  enu_nc->Sumw2();
206  enu_numu_nue->Sumw2();
207  enu_nue_nue->Sumw2();
208  enu_numu_numu->Sumw2();
209  enu_nue_numu->Sumw2();
210  enu_numu_nutau->Sumw2();
211  enu_nue_nutau->Sumw2();
212 
213 
214  char name[5][100] = {"#nu_{e}^{osc}","NC","#nu_{#mu} CC","#nu_{e}^{beam}","#nu_{#tau} CC"};
215 
216  for (int i = 0; i<nSamples; ++i){
217  events_truth[i] = 0;
218  events_reco[i] = 0;
219  enu[i] = tfs->make<TH1D>(Form("enu_%d",i),name[i],80,0,20);
220  enu[i]->Sumw2();
221  enu_osc[i] = tfs->make<TH1D>(Form("enu_osc_%d",i),name[i],160,0,20);
222  enu_osc[i]->Sumw2();
223  hvtxx[i] = tfs->make<TH1D>(Form("hvtxx_%d",i), name[i],100,-10,10);
224  hvtxy[i] = tfs->make<TH1D>(Form("hvtxy_%d",i), name[i],100,-10,10);
225  hvtxz[i] = tfs->make<TH1D>(Form("hvtxz_%d",i), name[i],100,-10,10);
226  hvtxx[i]->Sumw2();
227  hvtxy[i]->Sumw2();
228  hvtxz[i]->Sumw2();
229  htrklen[i] = tfs->make<TH1D>(Form("htrklen_%d",i), name[i],100,0,600);
230  htrklen[i]->Sumw2();
231  htrkdedx[i] = tfs->make<TH1D>(Form("htrkdedx_%d",i), name[i], 100,0,20);
232  htrkdedx[i]->Sumw2();
233  hrch[i] = tfs->make<TH1D>(Form("hrch_%d",i),name[i],100,0,1);
234  hrch[i]->Sumw2();
235  hrt[i] = tfs->make<TH1D>(Form("hrt_%d",i),name[i],100,0,1);
236  hrt[i]->Sumw2();
237  hpida[i] = tfs->make<TH1D>(Form("hpida_%d",i),name[i],100,0,50);
238  hpida[i]->Sumw2();
239  hnshw[i] = tfs->make<TH1D>(Form("hnshw_%d",i),name[i],20,0,20);
240  hnshw[i]->Sumw2();
241  heshw[i] = tfs->make<TH1D>(Form("heshw_%d",i),name[i],100,0,10);
242  heshw[i]->Sumw2();
243  hshwdedx[i] = tfs->make<TH1D>(Form("hshwdedx_%d",i),name[i],100,0,20);
244  hshwdedx[i]->Sumw2();
245  hdisx[i] = tfs->make<TH1D>(Form("hdisx_%d",i),name[i],100,-20,20);
246  hdisx[i]->Sumw2();
247  hdisy[i] = tfs->make<TH1D>(Form("hdisy_%d",i),name[i],100,-20,20);
248  hdisy[i]->Sumw2();
249  hdisz[i] = tfs->make<TH1D>(Form("hdisz_%d",i),name[i],100,-20,20);
250  hdisz[i]->Sumw2();
251  hfrshower[i] = tfs->make<TH1D>(Form("hfrshower_%d",i),name[i],110,0,1.1);
252  hfrshower[i]->Sumw2();
253  hnhitspershw[i] = tfs->make<TH1D>(Form("hnhitspershw_%d",i),name[i],100,0,20);
254  hnhitspershw[i]->Sumw2();
255  hevtcharge[i] = tfs->make<TH1D>(Form("hevtcharge_%d",i),name[i],100,0,3e6);
256  hevtcharge[i]->Sumw2();
257  hfr100w[i] = tfs->make<TH1D>(Form("hfr100w_%d",i),name[i],110,0,1.1);
258  hfr100w[i]->Sumw2();
259  }
260 
261  fPOT = tfs->make<TTree>("pottree","pot tree");
262  fPOT->Branch("pot",&pot,"pot/D");
263  fPOT->Branch("run",&run,"run/I");
264  fPOT->Branch("subrun",&subrun,"subrun/I");
265 
266 }
267 
268 dunemva::MVAAlg::~MVAAlg(){
269  fFile.close();
270 }
271 
272 //--------------------------------------------------------------------------------
273 void dunemva::MVAAlg::reconfigure(fhicl::ParameterSet const& p){
274 
275  fRawDigitModuleLabel = p.get< std::string >("RawDigitModuleLabel");
276  fWireModuleLabel = p.get< std::string >("WireModuleLabel");
277  fHitsModuleLabel = p.get< std::string >("HitsModuleLabel");
278  fTrackModuleLabel = p.get< std::string >("TrackModuleLabel");
279  fShowerModuleLabel = p.get< std::string >("ShowerModuleLabel");
280  fClusterModuleLabel = p.get< std::string >("ClusterModuleLabel");
281  fVertexModuleLabel = p.get< std::string >("VertexModuleLabel");
282  fGenieGenModuleLabel = p.get< std::string >("GenieGenModuleLabel");
283  fPOTModuleLabel = p.get< std::string >("POTModuleLabel");
284  fFlashModuleLabel = p.get< std::string >("FlashModuleLabel");
285  fCalorimetryModuleLabel = p.get< std::string >("CalorimetryModuleLabel");
286  fSelect = p.get< std::string >("Select");
287  fBeamMode = p.get< std::string >("BeamMode","FHC");
288  fFidVolCut = p.get< double >("FidVolCut");
289 }
290 
291 //--------------------------------------------------------------------------------
292 void dunemva::MVAAlg::Run( const art::Event & evt, double& result, double& wgt ){
293 
294  this->ResetVars();
295  this->PrepareEvent(evt); // does not reset vars, make sure to reset after evaluating
296  this->CalculateInputs();
297  if(!fMakeWeightTree){
298 
299  if (isinfidvol){
300  result = fReader.EvaluateMVA(fMVAMethod);
301  mf::LogVerbatim("MVASelect") << fMVAMethod
302  << " returned " << result;
303 
304  // Fill histogram of MVA value for each type
305 
306  // itype is... (see instantiation of "name")
307  // 0 for oscillated NuE CC
308  // 1 for NC background
309  // 2 for NuMu CC (oscillated component negligible)
310  // 3 for beam NuE
311  // 4 for NuTau CC
312 
313  if (itype==0) mva_nue_osc->Fill(result,1.);
314  if (itype==1) mva_nc->Fill(result,1.);
315  if (itype==2) mva_numu->Fill(result,1.);
316  if (itype==3) mva_nue_beam->Fill(result,1.);
317  if (itype==4) mva_nutau->Fill(result,1.);
318 
319  if (ccnc_truth==1){
320  }
321  else if (std::abs(pntype_flux)==14&&std::abs(nuPDG_truth)==12){
322  mva_numu_nue->Fill(result, oscpro);
323  }
324  else if (std::abs(pntype_flux)==12&&std::abs(nuPDG_truth)==12){
325  mva_nue_nue->Fill(result, oscpro);
326  }
327  else if (std::abs(pntype_flux)==14&&std::abs(nuPDG_truth)==14){
328  mva_numu_numu->Fill(result, oscpro);
329  }
330  else if (std::abs(pntype_flux)==12&&std::abs(nuPDG_truth)==14){
331  mva_nue_numu->Fill(result, oscpro);
332  }
333  else if (std::abs(pntype_flux)==14&&std::abs(nuPDG_truth)==16){
334  mva_numu_nutau->Fill(result, oscpro);
335  }
336  else if (std::abs(pntype_flux)==12&&std::abs(nuPDG_truth)==16){
337  mva_nue_nutau->Fill(result, oscpro);
338  }
339  }
340  else{//not in fiducial volume
341  result = -2.;
342  }
343  if (isinfidvoltruth){
344  if (ccnc_truth==1){
345  enu_nc->Fill(enu_truth, oscpro);
346  }
347  else if (std::abs(pntype_flux)==14&&std::abs(nuPDG_truth)==12){
348  enu_numu_nue->Fill(enu_truth, oscpro);
349  }
350  else if (std::abs(pntype_flux)==12&&std::abs(nuPDG_truth)==12){
351  enu_nue_nue->Fill(enu_truth, oscpro);
352  }
353  else if (std::abs(pntype_flux)==14&&std::abs(nuPDG_truth)==14){
354  enu_numu_numu->Fill(enu_truth, oscpro);
355  }
356  else if (std::abs(pntype_flux)==12&&std::abs(nuPDG_truth)==14){
357  enu_nue_numu->Fill(enu_truth, oscpro);
358  }
359  else if (std::abs(pntype_flux)==14&&std::abs(nuPDG_truth)==16){
360  enu_numu_nutau->Fill(enu_truth, oscpro);
361  }
362  else if (std::abs(pntype_flux)==12&&std::abs(nuPDG_truth)==16){
363  enu_nue_nutau->Fill(enu_truth, oscpro);
364  }
365  }
366  }
367 
368  wgt = weight;
369 }
370 
371 
372 //--------------------------------------------------------------------------------
373 void dunemva::MVAAlg::endSubRun(const art::SubRun& sr){
374  auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
375  if (potListHandle)
376  pot = potListHandle->totpot;
377  else
378  pot = 0.;
379  if (fPOT) fPOT->Fill();
380 }
381 
382 
383 //--------------------------------------------------------------------------------
384 void dunemva::MVAAlg::CalculateInputs( ){
385 
386  itype = -99999;
387  fOscPro = -99999;
388  weight = -99999;
389  evtcharge = -99999;
390  ntrack = -99999;
391  maxtrklength = -99999;
392  avgtrklength = -99999;
393  trkdedx = -99999;
394  trkrch = -99999;
395  trkrt = -99999;
396  trkfr = -99999;
397  trkpida_save = -99999;
398  nshower = -99999;
399  showerdedx = -99999;
400  eshower = -99999;
401  frshower = -99999;
402  nhitspershw = -99999;
403  shwlength = -99999;
404  shwmax = -99999;
405  fract_5_wires = -99999;
406  fract_10_wires = -99999;
407  fract_50_wires = -99999;
408  fract_100_wires = -99999;
409  shwdis = -99999;
410  ET = 0;
411 
412  // itype is... (see instantiation of "name")
413  // 0 for NuE signal
414  // 1 for NC background
415  // 2 for NuMu CC
416  // 3 for beam NuE
417  // 4 for NuTau CC
418 
419  // ccnc = ccnc_truth;
420  // pntype_flux is the original type
421  // nuPDG_truth is the type at the FD
422  //NuPdg = std::abs(nuPDG_truth);
423 
424  itype = -1;
425  if (ccnc_truth==1){
426  itype = 1; // NC background
427  }
428  else if (std::abs(pntype_flux)!=12&&std::abs(nuPDG_truth)==12){
429  // oscillated NuE (from NuMu)
430  itype = 0;
431  }
432  else if (std::abs(nuPDG_truth)==12){
433  // intrinsic NuE
434  itype = 3;
435  }
436  else if (std::abs(nuPDG_truth)==14){
437  // intrinsic NuMu (oscillated from NuE is negligible)
438  itype = 2;
439  }
440  else if (std::abs(nuPDG_truth)==16){
441  // NuTau background
442  itype = 4;
443  }
444  if (itype == -1){
445  //std::cout << "Unknown type: "
446  // << pntype_flux << " to " << nuPDG_truth
447  // << " (ccnc=" << ccnc_truth << ")"
448  // << std::endl;
449  }
450 
451  oscpro = this->OscPro(ccnc_truth,pntype_flux,nuPDG_truth,enu_truth);
452  float norm = this->Norm(ccnc_truth,pntype_flux,nuPDG_truth, subrun);
453  weight = oscpro*norm;
454  fOscPro = oscpro;
455 
456 
457  // Choose the most upstream vertex, even if it is just
458  // the beginning of a single track
459  bool findvtx = false;
460  float vtxx = 0;
461  float vtxy = 0;
462  float vtxz = 100000;
463  if (nvtx>0){
464  for (int i = 0; i<nvtx&&i<kMaxVertices; ++i){
465  if (vtx[i][2]<vtxz){
466  vtxx = vtx[i][0];
467  vtxy = vtx[i][1];
468  vtxz = vtx[i][2];
469  findvtx = true;
470  }
471  }
472  }
473  else if (ntracks_reco>0){
474  for (int i = 0; i<ntracks_reco; ++i){
475  if (trkstartz[i]<vtxz){
476  vtxx = trkstartx[i];
477  vtxy = trkstarty[i];
478  vtxz = trkstartz[i];
479  findvtx = true;
480  }
481  }
482  }
483 
484  if (std::abs(nuvtxx_truth)<360-50&&
485  std::abs(nuvtxy_truth)<600-50&&
486  nuvtxz_truth>50&&nuvtxz_truth<1394-150){
487  isinfidvoltruth = 1;
488  //if (enu_truth>20) continue;
489  //if (enu_truth<0.5) continue;
490  //if (ccnc_truth==1&&Y_truth*enu_truth<0.5) continue;
491  events_truth[itype]+=oscpro*norm;
492  enu[itype]->Fill(enu_truth,norm);
493  enu_osc[itype]->Fill(enu_truth,oscpro*norm);
494  if (findvtx){
495  hvtxx[itype]->Fill(vtxx-nuvtxx_truth, oscpro*norm);
496  hvtxy[itype]->Fill(vtxy-nuvtxy_truth, oscpro*norm);
497  hvtxz[itype]->Fill(vtxz-nuvtxz_truth, oscpro*norm);
498  }
499  }
500 
501  // get avg trk length
502  avgtrklength = 0.;
503  for (int i = 0; i<ntracks_reco; ++i){
504  avgtrklength += trklen[i];
505  }
506  if(ntracks_reco)
507  avgtrklength /= ntracks_reco;
508  ntrack = 1.*ntracks_reco;
509 
510  // Try to approximate transverse energy
511  unsigned int ntracks_notvtx(0);
512  if(findvtx && ntracks_reco>0){
513 
514  double Edir[3] = {0.,0.,0.};
515  double Etot(0.);
516  for (int i = 0; i<nshws; ++i){
517 
518  // don't count particles not out of event vertex
519  // VERY rough approximation here..
520  if( std::sqrt( std::pow(vtxx-shwstartx[i],2)
521  + std::pow(vtxy-shwstarty[i],2)
522  + std::pow(vtxz-shwstartz[i],2) > 2. ) ) continue;
523 
524  Edir[0] += std::abs(shwenergy[i][shwbestplane[i]])*shwdcosx[i];
525  Edir[1] += std::abs(shwenergy[i][shwbestplane[i]])*shwdcosy[i];
526  Edir[2] += std::abs(shwenergy[i][shwbestplane[i]])*shwdcosz[i];
527  Etot += std::abs(shwenergy[i][shwbestplane[i]]);
528  }
529  for (int i = 0; i<ntracks_reco; ++i){
530 
531  // don't count particles not out of event vertex
532  // VERY rough approximation here..
533  if( std::sqrt( std::pow(vtxx-trkstartx[i],2)
534  + std::pow(vtxy-trkstarty[i],2)
535  + std::pow(vtxz-trkstartz[i],2) > 2. ) ){
536  //if(trklen[i]<140) ntracks_notvtx++; this made it worse
537  ntracks_notvtx++;
538  continue;
539  }
540 
541  Edir[0] += std::abs(trkke[i][trkbestplane[i]])*trkstartdcosx[i];
542  Edir[1] += std::abs(trkke[i][trkbestplane[i]])*trkstartdcosy[i];
543  Edir[2] += std::abs(trkke[i][trkbestplane[i]])*trkstartdcosz[i];
544  Etot += std::abs(trkke[i][trkbestplane[i]]);
545  }
546 
547  if(Etot!=0.){
548  Edir[0] /= Etot;
549  Edir[1] /= Etot;
550  Edir[2] /= Etot;
551  }
552 
553  ET = std::sqrt( Edir[0]*Edir[0] + Edir[1]*Edir[1] );
554  //ET /= Etot;
555  if(ET>1||ET<0){
556  std::cout << ET << ", " << Etot << std::endl;
557  }
558 
559  if( ET > 0.9 ){
560  fFile << run << " " << subrun << " " << event << "\n";
561  }
562 
563  } else {
564  //overflow
565  ET = 0;
566  }
567 
568  // making loose assumptions about origin and dimensions...
569  if( findvtx&&
570  std::abs(vtxx) < 360-50 &&
571  std::abs(vtxy) < 600-50 &&
572  vtxz > 50 &&
573  vtxz < 1394-150 ){//in the fiducial volume
574 
575  isinfidvol = 1;
576  events_reco[itype] += oscpro*norm;
577  //ntrack = ntracks_reco;
578  //ntrack = ntracks_notvtx;
579  maxtrklength = 0;
580  int itrk = -1;
581  for (int i = 0; i<ntracks_reco; ++i){
582  if (trklen[i]>maxtrklength){
583  maxtrklength = trklen[i];
584  itrk = i;
585  }
586  }
587  //if (itype==2&&maxtrklength<5) std::cout<<run<<" "<<subrun<<" "<<event<<std::endl;
588  htrklen[itype]->Fill(maxtrklength,oscpro*norm);
589  //if (maxtrklength>300) continue;
590 
591  //std::vector<std::map<int,float> > vtrkdedx(3); don't use until Reco saves dEds per hit
592  std::vector<std::vector<std::map<int,float> > > vtrktick(fGeom->NTPC());
593  for (size_t i = 0; i<fGeom->NTPC(); ++i){
594  vtrktick[i].resize(3);
595  }
596 
597  trkdedx = 0;
598  trkrch = 0;
599  trkrt = 0;
600  trkfr = 0;
601  trkpida_save = 0;
602  evtcharge = 0;
603 
604  nshower = 1.*nshws;
605  hnshw[itype]->Fill(nshws,oscpro*norm);
606 
607  eshower = 0;
608  showerdedx = 0;
609  int ishw = -1;
610  for (int i = 0; i<nshws; ++i){
611  if (shwbestplane[i]>=0&&shwbestplane[i]<3){
612  if (shwenergy[i][shwbestplane[i]]>eshower){
613  eshower = shwenergy[i][shwbestplane[i]];
614  showerdedx = shwdedx[i][shwbestplane[i]];
615  ishw = i;
616  }
617  }
618  }
619  shwdis = 1000;
620  shwdisx = 1000;
621  shwdisy = 1000;
622  shwdisz = 1000;
623  shwcosx = -2;
624  shwcosy = -2;
625  shwcosz = -2;
626  if (ishw!=-1){
627  shwdis = sqrt(pow(shwstartx[ishw]-vtxx,2)+
628  pow(shwstarty[ishw]-vtxy,2)+
629  pow(shwstartz[ishw]-vtxz,2));
630  shwdisx = shwstartx[ishw]-vtxx;
631  shwdisy = shwstarty[ishw]-vtxy;
632  shwdisz = shwstartz[ishw]-vtxz;
633  shwcosx = shwdcosx[ishw];
634  shwcosy = shwdcosy[ishw];
635  shwcosz = shwdcosz[ishw];
636  }
637  trkcosx = -2;
638  trkcosy = -2;
639  trkcosz = -2;
640  if (itrk!=-1){
641  trkcosx = trkstartdcosx[itrk];
642  trkcosy = trkstartdcosy[itrk];
643  trkcosz = trkstartdcosz[itrk];
644  }
645  int shwwire0 = 100000;
646  int shwwire1 = -1;
647  int offset = 0;
648  //int lastwire = -1;
649  for (int i = 0; i<nhits && i<40000; ++i){
650  if (itrk!=-1){
651  if (hit_trkkey[i]==itrk){
652  // dont use hit_dEds until is is saved in Reco
653  // for now it is alwasy -9999
654  //vtrkdedx[hit_plane[i]][hit_wire[i]] = hit_dEds[i];
655  vtrktick[hit_tpc[i]][hit_plane[i]][hit_wire[i]] = hit_peakT[i];
656  }
657  }
658  if (hit_plane[i]==2){
659  evtcharge += hit_charge[i]*exp(hit_peakT[i]*0.5/taulife);
660  //evtcharge += hit_charge[i];
661  // if (run==20000000&&event==4137){
662  // std::cout<<hit_wire[i]<<" "<<shwwire0<<" "<<shwwire1<<std::endl;
663  // }
664  if (hit_shwkey[i]==ishw){
665  //if (lastwire == -1) lastwire = hit_wire[i];
666  //if (hit_wire[i]<lastwire) offset = 479;
667  offset = (hit_tpc[i]/4)*fGeom->Nwires(2);
668  if (hit_wire[i]+offset<shwwire0){
669  shwwire0 = hit_wire[i]+offset;
670  }
671  if (hit_wire[i]+offset>shwwire1){
672  shwwire1 = hit_wire[i]+offset;
673  }
674  //lastwire = hit_wire[i];
675  // if (hit_trkkey[i]>=0){
676  // if (sqrt(pow(trkstartx[hit_trkkey[i]]-vtxx,2)+
677  // pow(trkstarty[hit_trkkey[i]]-vtxy,2)+
678  // pow(trkstartz[hit_trkkey[i]]-vtxz,2))<shwdis)
679  // shwdis = sqrt(pow(trkstartx[hit_trkkey[i]]-vtxx,2)+
680  // pow(trkstarty[hit_trkkey[i]]-vtxy,2)+
681  // pow(trkstartz[hit_trkkey[i]]-vtxz,2));
682  // if (sqrt(pow(trkendx[hit_trkkey[i]]-vtxx,2)+
683  // pow(trkendy[hit_trkkey[i]]-vtxy,2)+
684  // pow(trkendz[hit_trkkey[i]]-vtxz,2))<shwdis)
685  // shwdis = sqrt(pow(trkendx[hit_trkkey[i]]-vtxx,2)+
686  // pow(trkendy[hit_trkkey[i]]-vtxy,2)+
687  // pow(trkendz[hit_trkkey[i]]-vtxz,2));
688  // }
689  }
690  }
691  }
692  hevtcharge[itype]->Fill(evtcharge,oscpro*norm);
693  std::vector<float> shwph;
694  if (shwwire1>=shwwire0){
695  shwph.resize(shwwire1-shwwire0+1,0);
696  shwlength = shwwire1-shwwire0+1;
697  }
698  else{
699  shwlength = 0;
700  }
701  // if (run==20000001&&event==4394){
702  // std::cout<<ishw<<" "<<shwwire0<<" "<<shwwire1<<std::endl;
703  // }
704 
705  if (itrk!=-1){ // itrk is track index of longest track
706  int ipl[2] = {-1,-1}; // top 2 with most hits
707  int maxhits[2] = {0,0};
708 
709  // find the two planes with the most hits
710  for (int i = 0; i<3; ++i){
711  // set maxhits[0]
712  int tothits = 0;
713  for (size_t j = 0; j<fGeom->NTPC(); ++j){
714  tothits += vtrktick[j][i].size();
715  }
716  if (tothits>maxhits[0]){ // vtrktick[i] is a map from wire number to charge
717  if (ipl[1]!=-1){ // dont do first time (i=0)
718  maxhits[1] = maxhits[0];
719  ipl[1] = ipl[0];
720  }
721  maxhits[0] = tothits;
722  ipl[0] = i;
723  }
724  // set maxhits[1], isn't called first time (i=0)
725  else if (tothits>maxhits[1]){
726  maxhits[1] = tothits;
727  ipl[1] = i;
728  }
729  }
730 
731  //int totwires = 0;
732  // int startw[2] = {-1,-1};
733  // for (int i = 0; i<2; ++i){
734  // if (ipl[i]==-1) continue;
735  // int startw[i] = int(vtrkdedx[ipl[i]].begin()->first+0.3*(vtrkdedx[ipl[i]].rbegin()->first-vtrkdedx[ipl[i]].begin()->first));
736  // }
737 
738  /* can't calculate using vtrkdedx yet since PANDORA didn't save hit_dEds for each hit
739  for (int i = 0; i<2; ++i){
740  if (ipl[i]==-1) continue;
741  for (auto const &j : vtrkdedx[ipl[i]]){
742  trkdedx += j.second;
743  ++totwires;
744  }
745  }
746  if (totwires) trkdedx/=(totwires);
747  */
748 
749  if(itrk==-1) trkdedx = 0;
750  else trkdedx = trkke[itrk][trkbestplane[itrk]]/trklen[itrk];
751 
752  float qall = -1;
753  float qtrk = -1;
754  float pida = -1;
755  float ipl_evtcharge = 0;
756  int npida = 0;
757  std::vector<float> vhitq;
758  for (int i = 0; i<nhits && i<40000; ++i){ // only 40000 hits saved after Reco
759  for (int j = 0; j<2; ++j){ // planes with top 2 # of hits (ipl[j])
760  if (hit_plane[i] == ipl[j]){
761  if (hit_trkkey[i] == itrk){
762  if (hit_resrange[i]<30&&hit_resrange[i]>0.6){
763  pida += hit_dEds[i]*pow(hit_resrange[i],0.42);
764  ++npida;
765  }
766  }
767  if (vtrktick[hit_tpc[i]][ipl[j]].find(hit_wire[i])!=vtrktick[hit_tpc[i]][ipl[j]].end()){
768  if (std::abs(vtrktick[hit_tpc[i]][ipl[j]][hit_wire[i]]-hit_peakT[i])<200){
769  vhitq.push_back(hit_charge[i]);
770  qall+=hit_charge[i];
771  if (hit_trkkey[i] == itrk){
772  qtrk+=hit_charge[i];
773  } // if hit is actually in the longest track
774  } // if event-hit is within 200 ticks of track-hit on this wire
775  } // if event-hit wire is in track-hit map (wire-->peak time), track has a hit on this wire
776 
777  // sum the denominator to get trk charge fraction of event charge
778  // (normal evtcharge sums only on collection plane)
779  ipl_evtcharge += hit_charge[i];
780 
781  } // if hit is in top 2 planes
782  } // planes with top 2 # of hits
783  } // all nhits in event
784 
785  // sort all of the hit charges for hits within a 200tick width of the track
786  std::sort(vhitq.begin(), vhitq.end());
787  float q1 = 0;
788  float q2 = 0;
789  for (size_t i = 0; i<vhitq.size(); ++i){
790  if (i<vhitq.size()/2){
791  q1+=vhitq[i];
792  }
793  else{
794  q2+=vhitq[i];
795  }
796  }
797  hrch[itype]->Fill(q1/q2,oscpro*norm);
798  hrt[itype]->Fill(TMath::Min(float(0.999999),qtrk/qall),oscpro*norm);
799  hpida[itype]->Fill(pida/npida,oscpro*norm);
800  if(q2!=0) trkrch = q1/q2;
801  trkrt = qtrk/qall;
802  if(ipl_evtcharge!=0) trkfr = qtrk/ipl_evtcharge;
803 
804  if(trkfr>1.){
805  std::cout << "Erroneous track charge fraction: " << trkfr
806  << " = " << qtrk << " / " << ipl_evtcharge << std::endl;
807  }
808 
809  //trkpida = pida/npida;
810  trkpida_save = trkpida[itrk][trkbestplane[itrk]];
811  if (trkpida_save>100.) trkpida_save = 100.;
812  if (trkpida_save<0.) trkpida_save = 0.;
813  }//itrk!=-1
814  else{
815  hrch[itype]->Fill(0.,oscpro*norm);
816  hrt[itype]->Fill(0.,oscpro*norm);
817  hpida[itype]->Fill(0.,oscpro*norm);
818  }
819  if (trkdedx<0.) trkdedx = 0.;
820  if (trkdedx>100.) trkdedx = 100.;
821  htrkdedx[itype]->Fill(trkdedx, oscpro*norm);
822 
823  frshower = 0;
824  int totalshwhits = 0;
825  std::map<int,int> shwwires;
826  offset = 0;
827  //lastwire = -1;
828 
829  for (int i = 0; i<nhits && i<40000; ++i){
830  if (hit_plane[i]==2&&hit_shwkey[i]==ishw){
831  frshower+=hit_charge[i]*exp(hit_peakT[i]*0.5/taulife);
832  shwwires[hit_wire[i]] = 1;
833  ++totalshwhits;
834  offset = (hit_tpc[i]/4)*fGeom->Nwires(2);
835  //if (lastwire ==-1) lastwire = hit_wire[i];
836  //if (hit_wire[i]<lastwire) offset = 479;
837  shwph[hit_wire[i]+offset-shwwire0] += hit_charge[i]*exp(hit_peakT[i]*0.5/taulife);
838  //lastwire = hit_wire[i];
839  }
840  }
841  frshower/=evtcharge;
842  if (shwwires.size()){
843  nhitspershw = float(totalshwhits)/shwwires.size();
844  }
845  else
846  nhitspershw = 0;
847 
848  shwmax = 0;
849  float maxshw = 0;
850  fract_5_wires = 0;
851  fract_10_wires = 0;
852  fract_50_wires = 0;
853  fract_100_wires = 0;
854  float totshwph = 0;
855  for (size_t i = 0; i<shwph.size(); ++i){
856  if (shwph[i]>maxshw){
857  shwmax = i;
858  maxshw = shwph[i];
859  }
860  totshwph += shwph[i];
861  float ph_5_wires = 0;
862  float ph_10_wires = 0;
863  float ph_50_wires = 0;
864  float ph_100_wires = 0;
865  for (size_t j = i; j<i+100&&j<shwph.size(); ++j){
866  if (j<i+5){
867  ph_5_wires += shwph[j];
868  }
869  if (j<i+10){
870  ph_10_wires += shwph[j];
871  }
872  if (j<i+50){
873  ph_50_wires += shwph[j];
874  }
875  if (j<i+100){
876  ph_100_wires += shwph[j];
877  }
878  }
879  if (ph_5_wires>fract_5_wires){
880  fract_5_wires = ph_5_wires;
881  }
882  if (ph_10_wires>fract_10_wires){
883  fract_10_wires = ph_10_wires;
884  }
885  if (ph_50_wires>fract_50_wires){
886  fract_50_wires = ph_50_wires;
887  }
888  if (ph_100_wires>fract_100_wires){
889  fract_100_wires = ph_100_wires;
890  }
891  }
892  if (shwlength) shwmax/=shwlength;
893  if (totshwph) {
894  fract_5_wires/=totshwph;
895  fract_10_wires/=totshwph;
896  fract_50_wires/=totshwph;
897  fract_100_wires/=totshwph;
898  }
899  //if (frshower<0.1&&itype==0) std::cout<<run<<" "<<subrun<<" "<<event<<std::endl;
900  if (ishw!=-1){
901  heshw[itype]->Fill(eshower,oscpro*norm);
902  if (showerdedx<0.) showerdedx = 0.;
903  if (showerdedx>100.) showerdedx = 100.;
904  if (eshower>0.5){
905  hshwdedx[itype]->Fill(showerdedx,oscpro*norm);
906  hfrshower[itype]->Fill(frshower,oscpro*norm);
907  hnhitspershw[itype]->Fill(nhitspershw,oscpro*norm);
908  hfr100w[itype]->Fill(fract_100_wires,oscpro*norm);
909  hdisx[itype]->Fill(shwstartx[ishw]-trkstartx[itrk],oscpro*norm);
910  hdisy[itype]->Fill(shwstarty[ishw]-trkstarty[itrk],oscpro*norm);
911  hdisz[itype]->Fill(shwstartz[ishw]-trkstartz[itrk],oscpro*norm);
912  }
913  }
914 
915 
916  // if (maxtrklength<300&&
917  // evtcharge>1e5&&
918  // evtcharge<5e6&&
919  // nshower>0&&
920  // eshower>0.5&&
921  // subrun<666){
922  // if (itype==0) fWeightTree[1]->Fill();
923  // if (itype==1) fWeightTree[2]->Fill();
924  // if (itype==2) fWeightTree[3]->Fill();
925  // fWeightTree[0]->Fill();
926  // }
927 
928  if(fMakeWeightTree){
929 
930  // Fill Weight Trees
931  //
932  // fWeightTree[0] is all
933  // fWeightTree[1] is signal
934  // fWeightTree[2] is NC background
935  // fWeightTree[3] is CC background
936 
937  // itype is... (see instantiation of "name")
938  // 0 for oscillated NuE CC
939  // 1 for NC background
940  // 2 for NuMu CC (oscillated component negligible)
941  // 3 for beam NuE
942  // 4 for NuTau CC
943 
944  if(fSelect=="numu"){
945  if (itype==2) fWeightTree[1]->Fill(); // signal (numucc)
946  if (itype==1) fWeightTree[2]->Fill(); // NC background
947  if (itype==0) fWeightTree[3]->Fill(); // CC background (nuecc)
948  }
949  else if(fSelect=="nue"){
950  if (itype==0) fWeightTree[1]->Fill(); // signal (nuecc)
951  if (itype==1) fWeightTree[2]->Fill(); // NC background
952  if (itype==2) fWeightTree[3]->Fill(); // CC background (numucc)
953  }
954  fWeightTree[0]->Fill();
955 
956  } // if making weight tree
957 
958  } //in fiducial volume
959  else {
960  mf::LogVerbatim("MVASelect") << " Not found in fiducial volume. nvtx=" << nvtx << " True vtx = ("
961  << nuvtxx_truth <<", "<< nuvtxy_truth <<", "<< nuvtxz_truth << "), Reco vtx = ("<<vtxx<<", "<<vtxy<<", "<<vtxz<<")";
962  }
963 
964  if(itype == -99999) mf::LogVerbatim("MVASelect") << " itype not set";
965  if(weight == -9999) mf::LogVerbatim("MVASelect") << " weight not set";
966  if(evtcharge == -9999) mf::LogVerbatim("MVASelect") << " evtcharge not set";
967  if(ntrack == -9999) mf::LogVerbatim("MVASelect") << " ntrack not set";
968  if(maxtrklength == -9999) mf::LogVerbatim("MVASelect") << " maxtrklength not set";
969  if(avgtrklength == -9999) mf::LogVerbatim("MVASelect") << " avgtrklength not set";
970  if(trkdedx == -9999) mf::LogVerbatim("MVASelect") << " trkdedx not set";
971  if(trkrch == -9999) mf::LogVerbatim("MVASelect") << " trkrch not set";
972  if(trkrt == -9999) mf::LogVerbatim("MVASelect") << " trkrt not set";
973  if(trkfr == -9999) mf::LogVerbatim("MVASelect") << " trkfr not set";
974  if(trkpida_save == -9999) mf::LogVerbatim("MVASelect") << " trkpida_save not set";
975  if(nshower == -9999) mf::LogVerbatim("MVASelect") << " nshower not set";
976  if(showerdedx == -9999) mf::LogVerbatim("MVASelect") << " showerdedx not set";
977  if(eshower == -9999) mf::LogVerbatim("MVASelect") << " eshower not set";
978  if(frshower == -9999) mf::LogVerbatim("MVASelect") << " frshower not set";
979  if(nhitspershw == -9999) mf::LogVerbatim("MVASelect") << " nhitspershw not set";
980  if(shwlength == -9999) mf::LogVerbatim("MVASelect") << " shwlength not set";
981  if(shwmax == -9999) mf::LogVerbatim("MVASelect") << " shwmax not set";
982  if(fract_5_wires == -9999) mf::LogVerbatim("MVASelect") << " fract_5_wiresnot set";
983  if(fract_10_wires == -9999) mf::LogVerbatim("MVASelect") << " fract_10_wiresnot set";
984  if(fract_50_wires == -9999) mf::LogVerbatim("MVASelect") << " fract_50_wiresnot set";
985  if(fract_100_wires == -9999) mf::LogVerbatim("MVASelect") << " fract_100_wiresnot set";
986  if(shwdis == -9999) mf::LogVerbatim("MVASelect") << " shwdisnot set";
987  //if(ET == 0) mf::LogVerbatim("MVASelect") << " ET not set";
988 
989 
990 } // CalculateInputs()
991 
992 
993 //--------------------------------------------------------------------------------
994 float dunemva::MVAAlg::Norm(int ccnc, int nu0, int nu1, int subrun){
995 
996  float norm = 1;
997  float mass = 814.308*1.4e3/1e6; //kt
998 
999  float pot_nu = 1.09125e+23;
1000  float pot_nue = 1.05555e+23;
1001  float pot_nutau = 2.83966e+23;
1002  if(fBeamMode == "RHC"){
1003  pot_nu = 1.89753e+23;
1004  pot_nue = 1.80693e+23;
1005  pot_nutau = 4.19518e+23;
1006  }
1007 
1008  float potpermwyr = 1.47e21/1.07;
1009  float pot = 150*potpermwyr/mass;//kt*mw*yr
1010 
1011  if (ccnc==1){
1012  norm = pot/(pot_nu+pot_nue+pot_nutau);
1013  }
1014  else if (std::abs(nu0)==14&&std::abs(nu1)==12){
1015  norm = pot/pot_nue;
1016  }
1017  else if (std::abs(nu0)==12&&std::abs(nu1)==12){
1018  norm = pot/pot_nu;
1019  }
1020  else if (std::abs(nu0)==14&&std::abs(nu1)==14){
1021  norm = pot/pot_nu;
1022  }
1023  else if (std::abs(nu0)==12&&std::abs(nu1)==14){
1024  norm = pot/pot_nutau;
1025  }
1026  else if (std::abs(nu0)==14&&std::abs(nu1)==16){
1027  norm = pot/pot_nutau;
1028  }
1029  else if (std::abs(nu0)==12&&std::abs(nu1)==16){
1030  norm = pot/pot_nue;
1031  }
1032  else{
1033  std::cout << "Unknown oscillation: "
1034  << nu0 << " to " << nu1 << std::endl;
1035  }
1036  return norm;
1037 } // Norm()
1038 
1039 
1040 //--------------------------------------------------------------------------------
1041 float dunemva::MVAAlg::OscPro(int ccnc, int nu0, int nu1, float NuE){
1042 
1043  float osc_dm2=0.002457;
1044  float osc_L=1300.;
1045  //float sinth23=pow(sin(0.67),2);//sin(3.1415926535/4.0);
1046  //float sin2th13=0.094;
1047 
1048  //FIXME osceq is created on the heap which was memory leak bound. For now, a delete has been called at the end of the function. This needs a more long term solution
1049  TF1 *osceq = 0;
1050  if(!osceq){
1051  osceq = new TF1("f2","sin(1.267*[0]*[1]/x)*sin(1.267*[0]*[1]/x)",0.,120.);
1052  osceq->SetParameters(osc_dm2,osc_L);
1053  }
1054 
1055  float OscProb = 0 ;
1056  float NumuToNutau ;
1057  float NumuToNue ;
1058  float NueSurvival ;
1059  float NumuSurvival ;
1060  float NueToNutau ;
1061  float NueToNumu;
1062 
1063  //float th13 = 0.156;
1064  //float th23 = 0.670;
1065  float th13 = 0.148;
1066  float th23 = 0.738;
1067 
1068  //float th13 = 0.;
1069  // float th23 = 0.;
1070 
1071  //NumuToNutau = 4.*sinth23*(1.-sinth23)*pow(1-sin2th13/4,2) ;
1072  //NumuToNutau *= osceq->Eval(TMath::Abs(NuE)) ;
1073 
1074  //NumuToNue = sinth23*sin2th13*osceq->Eval(TMath::Abs(NuE)) ;
1075  NumuToNue = pow(sin(th23),2)*pow(sin(2*th13),2)*osceq->Eval(TMath::Abs(NuE)) ;
1076  //NueSurvival = 1.- sin2th13*osceq->Eval(TMath::Abs(NuE)) ;
1077  NueSurvival = 1-pow(sin(2*th13),2)*osceq->Eval(TMath::Abs(NuE)) ;
1078  //NueSurvival = 1.;
1079 
1080  //NumuSurvival = 1. - NumuToNutau - NumuToNue ;
1081  //NumuSurvival = 1.;
1082  NumuSurvival = 1-(pow(cos(th13),2)*pow(sin(2*th23),2)+pow(sin(th23),4)*pow(sin(2*th13),2))*osceq->Eval(TMath::Abs(NuE)) ;
1083 
1084  //NumuSurvival = 1.;
1085  //NueToNutau = (1.-sinth23)*sin2th13*osceq->Eval(TMath::Abs(NuE)) ;
1086 
1087  NumuToNutau = 1 - NumuSurvival - NumuToNue;
1088  NueToNutau = 1 - NueSurvival - NumuToNue;
1089 
1090  NueToNumu = NumuToNue;
1091 
1092  if (ccnc==1){
1093  OscProb = 1;
1094  }
1095  else if (std::abs(nu0)==14&&std::abs(nu1)==12){
1096  OscProb = NumuToNue;
1097  }
1098  else if (std::abs(nu0)==12&&std::abs(nu1)==12){
1099  OscProb = NueSurvival;
1100  }
1101  else if (std::abs(nu0)==14&&std::abs(nu1)==14){
1102  OscProb = NumuSurvival;
1103  }
1104  else if (std::abs(nu0)==12&&std::abs(nu1)==14){
1105  OscProb = NueToNumu;
1106  }
1107  else if (std::abs(nu0)==14&&std::abs(nu1)==16){
1108  OscProb = NumuToNutau;
1109  }
1110  else if (std::abs(nu0)==12&&std::abs(nu1)==16){
1111  OscProb = NueToNutau;
1112  }
1113  else{
1114  std::cout << "Unknown oscillation: "
1115  << nu0 << " to " << nu1 << std::endl;
1116  }
1117  delete osceq;
1118 
1119  return OscProb;
1120 } // OscPro()
1121 
1122 
1123 
1124 //--------------------------------------------------------------------------------
1125 void dunemva::MVAAlg::PrepareEvent(const art::Event& evt){
1126 
1127  //std::cout << " ~~~~~~~~~~~~~~~ MVA: Getting Event Reco ~~~~~~~~~~~~~~ " << std::endl;
1128 
1129  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
1130  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
1131  const sim::ParticleList& plist = pi_serv->ParticleList();
1132 
1133  run = evt.run();
1134  subrun = evt.subRun();
1135  event = evt.id().event();
1136  art::Timestamp ts = evt.time();
1137  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
1138  evttime = tts.AsDouble();
1139  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
1140  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
1141  taulife = detProp.ElectronLifetime();
1142  isdata = evt.isRealData();
1143 
1144  // * Raw Digits
1145  std::vector<art::Ptr<raw::RawDigit> > rawlist;
1146  auto rawListHandle = evt.getHandle<std::vector<raw::RawDigit> >(fRawDigitModuleLabel);
1147  if (rawListHandle)
1148  art::fill_ptr_vector(rawlist, rawListHandle);
1149 
1150  // * wires
1151  std::vector<art::Ptr<recob::Wire>> wirelist;
1152  auto wireListHandle = evt.getHandle< std::vector<recob::Wire>>(fWireModuleLabel);
1153  if (wireListHandle)
1154  art::fill_ptr_vector(wirelist, wireListHandle);
1155 
1156  // * hits
1157  std::vector<art::Ptr<recob::Hit> > hitlist;
1158  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
1159  if (hitListHandle)
1160  art::fill_ptr_vector(hitlist, hitListHandle);
1161 
1162  // * tracks
1163  std::vector<art::Ptr<recob::Track> > tracklist;
1164  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
1165  if (trackListHandle)
1166  art::fill_ptr_vector(tracklist, trackListHandle);
1167 
1168  // * vertices
1169  std::vector<art::Ptr<recob::Vertex> > vtxlist;
1170  auto vtxListHandle = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel);
1171  if (vtxListHandle)
1172  art::fill_ptr_vector(vtxlist, vtxListHandle);
1173 
1174  // * showers
1175  std::vector<art::Ptr<recob::Shower>> shwlist;
1176  auto shwListHandle = evt.getHandle<std::vector<recob::Shower>>(fShowerModuleLabel);
1177  if (shwListHandle)
1178  art::fill_ptr_vector(shwlist, shwListHandle);
1179 
1180  // * flashes
1181  std::vector<art::Ptr<recob::OpFlash> > flashlist;
1182  auto flashListHandle = evt.getHandle< std::vector<recob::OpFlash> >(fFlashModuleLabel);
1183  if (flashListHandle)
1184  art::fill_ptr_vector(flashlist, flashListHandle);
1185 
1186  // * associations
1187  art::FindManyP<recob::Hit> fmth(trackListHandle, evt, fTrackModuleLabel);
1188  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt, fTrackModuleLabel);
1189  art::FindManyP<recob::SpacePoint> fmhs(hitListHandle, evt, fTrackModuleLabel);
1190  art::FindMany<anab::Calorimetry> fmcal(trackListHandle, evt, fCalorimetryModuleLabel);
1191 
1192  // charge from raw digits
1193  rawcharge = 0;
1194  /* Comment for now as it is too slow
1195  for (size_t i = 0; i<rawlist.size(); ++i){
1196  if (fGeom->SignalType(rawlist[i]->Channel()) == geo::kCollection){
1197  double pedestal = rawlist[i]->GetPedestal();
1198  for (size_t j = 0; j<rawlist[i]->NADC(); ++j){
1199  rawcharge += rawlist[i]->ADC(j)-pedestal;
1200  }
1201  }
1202  }
1203  */
1204  //charge from wires
1205  wirecharge = 0;
1206  for (size_t i = 0; i<wirelist.size(); ++i){
1207  if (fGeom->SignalType(wirelist[i]->Channel()) == geo::kCollection){
1208  const recob::Wire::RegionsOfInterest_t& signalROI = wirelist[i]->SignalROI();
1209  for(const auto& range : signalROI.get_ranges()){
1210  const std::vector<float>& signal = range.data();
1211  raw::TDCtick_t roiFirstBinTick = range.begin_index();
1212  for (size_t j = 0; j<signal.size(); ++j){
1213  wirecharge += signal[j]*exp((j+roiFirstBinTick)*0.5/taulife);
1214  }
1215  }
1216  }
1217  }
1218 
1219  //hit information
1220  nhits = hitlist.size();
1221  nhits_stored = std::min(nhits, kMaxHits);
1222  for (int i = 0; i < nhits && i < kMaxHits ; ++i){//loop over hits
1223  hit_channel[i] = hitlist[i]->Channel();
1224  hit_plane[i] = hitlist[i]->WireID().Plane;
1225  hit_wire[i] = hitlist[i]->WireID().Wire;
1226  hit_tpc[i] = hitlist[i]->WireID().TPC;
1227  hit_peakT[i] = hitlist[i]->PeakTime();
1228  hit_charge[i] = hitlist[i]->Integral();
1229  hit_summedADC[i] = hitlist[i]->SummedADC();
1230  hit_startT[i] = hitlist[i]->PeakTimeMinusRMS();
1231  hit_endT[i] = hitlist[i]->PeakTimePlusRMS();
1232 
1233  }
1234 
1235  //track information
1236  ntracks_reco=tracklist.size();
1237 
1238  recob::Track::Vector_t larStart;
1239  recob::Track::Vector_t larEnd;
1240  for(int i=0; i<std::min(int(tracklist.size()),kMaxTrack);++i){
1241  recob::Track::Point_t trackStart, trackEnd;
1242  std::tie(trackStart, trackEnd) = tracklist[i]->Extent();
1243  larStart = tracklist[i]->VertexDirection();
1244  larEnd = tracklist[i]->EndDirection();
1245 
1246  trkid[i] = tracklist[i]->ID();
1247  trkstartx[i] = trackStart.X();
1248  trkstarty[i] = trackStart.Y();
1249  trkstartz[i] = trackStart.Z();
1250  trkendx[i] = trackEnd.X();
1251  trkendy[i] = trackEnd.Y();
1252  trkendz[i] = trackEnd.Z();
1253  trkstartdcosx[i] = larStart.X();
1254  trkstartdcosy[i] = larStart.Y();
1255  trkstartdcosz[i] = larStart.Z();
1256  trkenddcosx[i] = larEnd.X();
1257  trkenddcosy[i] = larEnd.Y();
1258  trkenddcosz[i] = larEnd.Z();
1259  trklen[i] = tracklist[i]->Length();
1260  if (fmthm.isValid()){
1261  auto vhit = fmthm.at(i);
1262  auto vmeta = fmthm.data(i);
1263  for (size_t h = 0; h < vhit.size(); ++h){
1264  if (vhit[h].key()<kMaxHits){
1265  hit_trkkey[vhit[h].key()] = tracklist[i].key();
1266  if (vmeta[h]->Dx()){
1267  hit_dQds[vhit[h].key()] = vhit[h]->Integral()*fCalorimetryAlg.LifetimeCorrection(clockData, detProp, vhit[h]->PeakTime())/vmeta[h]->Dx();
1268  hit_dEds[vhit[h].key()] = fCalorimetryAlg.dEdx_AREA(clockData, detProp, *vhit[h], vmeta[h]->Dx());
1269  }
1270  hit_resrange[vhit[h].key()] = tracklist[i]->Length(vmeta[h]->Index());
1271  }
1272  }//loop over all hits
1273  }//fmthm is valid
1274  else if (fmth.isValid()){
1275  std::vector< art::Ptr<recob::Hit> > vhit = fmth.at(i);
1276  for (size_t h = 0; h < vhit.size(); ++h){
1277  if (vhit[h].key()<kMaxHits){
1278  hit_trkkey[vhit[h].key()] = tracklist[i].key();
1279  }
1280  }
1281  }
1282  if (fmcal.isValid()){
1283  unsigned maxnumhits = 0;
1284  std::vector<const anab::Calorimetry*> calos = fmcal.at(i);
1285  for (auto const& calo : calos){
1286  if (calo->PlaneID().isValid){
1287  trkke[i][calo->PlaneID().Plane] = calo->KineticEnergy();
1288  if (calo->dEdx().size()>maxnumhits){
1289  maxnumhits = calo->dEdx().size();
1290  trkbestplane[i] = calo->PlaneID().Plane;
1291  }
1292  double pida = 0;
1293  int used_trkres = 0;
1294  for (size_t ip = 0; ip<calo->dEdx().size(); ++ip){
1295  if (calo->ResidualRange()[ip]<30){
1296  pida += calo->dEdx()[ip]*pow(calo->ResidualRange()[ip],0.42);
1297  ++used_trkres;
1298  }
1299  }
1300  if (used_trkres) pida/=used_trkres;
1301  trkpida[i][calo->PlaneID().Plane] = pida;
1302  }
1303  }
1304  }
1305  if (!isdata&&fmth.isValid()){
1306  // Find true track for each reconstructed track
1307  int TrackID = 0;
1308  std::vector< art::Ptr<recob::Hit> > allHits = fmth.at(i);
1309 
1310  std::map<int,double> trkide;
1311  for(size_t h = 0; h < allHits.size(); ++h){
1312  art::Ptr<recob::Hit> hit = allHits[h];
1313  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
1314  for(size_t e = 0; e < TrackIDs.size(); ++e){
1315  trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
1316  }
1317  }
1318  // Work out which IDE despoited the most charge in the hit if there was more than one.
1319  double maxe = -1;
1320  double tote = 0;
1321  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
1322  tote += ii->second;
1323  if ((ii->second)>maxe){
1324  maxe = ii->second;
1325  TrackID = ii->first;
1326  }
1327  }
1328  // Now have trackID, so get PdG code and T0 etc.
1329  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
1330  if (particle){
1331  trkg4id[i] = TrackID;
1332  trkg4pdg[i] = particle->PdgCode();
1333  trkg4startx[i] = particle->Vx();
1334  trkg4starty[i] = particle->Vy();
1335  trkg4startz[i] = particle->Vz();
1336  float sum_energy = 0;
1337  int numhits = 0;
1338  //std::map<float,float> hite;
1339  double x = 0;
1340  double y = 0;
1341  double z = 0;
1342  double mindis = 1e10;
1343  //find the closest point to the neutrino vertex
1344  for(size_t h = 0; h < allHits.size(); ++h){
1345  art::Ptr<recob::Hit> hit = allHits[h];
1346  if (hit->WireID().Plane==2){
1347  std::vector<art::Ptr<recob::SpacePoint> > spts = fmhs.at(hit.key());
1348  if (spts.size()){
1349  double dis = sqrt(pow(spts[0]->XYZ()[0]-trkg4startx[i],2)+
1350  pow(spts[0]->XYZ()[1]-trkg4starty[i],2)+
1351  pow(spts[0]->XYZ()[2]-trkg4startz[i],2));
1352  if (dis<mindis){
1353  mindis = dis;
1354  x = spts[0]->XYZ()[0];
1355  y = spts[0]->XYZ()[1];
1356  z = spts[0]->XYZ()[2];
1357  }
1358  }
1359  }
1360  }
1361  for(size_t h = 0; h < allHits.size(); ++h){
1362  art::Ptr<recob::Hit> hit = allHits[h];
1363  if (hit->WireID().Plane==2){
1364  std::vector<art::Ptr<recob::SpacePoint> > spts = fmhs.at(hit.key());
1365  if (spts.size()){
1366  if (sqrt(pow(spts[0]->XYZ()[0]-x,2)+
1367  pow(spts[0]->XYZ()[1]-y,2)+
1368  pow(spts[0]->XYZ()[2]-z,2))<3){
1369  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
1370  float toten = 0;
1371  for(size_t e = 0; e < TrackIDs.size(); ++e){
1372  //sum_energy += TrackIDs[e].energy;
1373  toten+=TrackIDs[e].energy;
1374  }
1375  if (toten){
1376  sum_energy += toten;
1377  ++numhits;
1378  }
1379  }
1380  }
1381  }
1382  }
1383 
1384  float pitch = 0;
1385  float dis1 = sqrt(pow(trkstartx[i]-trkg4startx[i],2)+
1386  pow(trkstarty[i]-trkg4starty[i],2)+
1387  pow(trkstartz[i]-trkg4startz[i],2));
1388  float dis2 = sqrt(pow(trkendx[i]-trkg4startx[i],2)+
1389  pow(trkendy[i]-trkg4starty[i],2)+
1390  pow(trkendz[i]-trkg4startz[i],2));
1391  if (dis1<dis2){
1392  try{
1393  pitch = lar::util::TrackPitchInView(*(tracklist[i]),geo::kZ,0);
1394  }
1395  catch(...){
1396  pitch = 0;
1397  }
1398  }
1399  else{
1400  try{
1401  pitch = lar::util::TrackPitchInView(*(tracklist[i]),geo::kZ,tracklist[i]->NumberTrajectoryPoints()-1);
1402  }
1403  catch(...){
1404  pitch = 0;
1405  }
1406  }
1407  if ( pitch && numhits ) {
1408  trkg4initdedx[i] = sum_energy/(numhits*pitch);
1409  }
1410  else{
1411  trkg4initdedx[i] = 0;
1412  }
1413  }//if (particle)
1414  }//MC
1415  }
1416 
1417  //vertex information
1418  nvtx = vtxlist.size();
1419  for (int i = 0; i < nvtx && i < kMaxVertices ; ++i){//loop over hits
1420  Double_t xyz[3] = {};
1421  vtxlist[i]->XYZ(xyz);
1422  for (size_t j = 0; j<3; ++j) vtx[i][j] = xyz[j];
1423  }
1424 
1425  //shower information
1426  if (shwListHandle.isValid()){
1427  art::FindManyP<recob::Hit> fmsh(shwListHandle, evt, fShowerModuleLabel);
1428 
1429  nshws = shwlist.size();
1430 
1431  for (int i = 0; i<std::min(int(shwlist.size()),kMaxShower); ++i){
1432  shwid[i] = shwlist[i]->ID();
1433  shwdcosx[i] = shwlist[i]->Direction().X();
1434  shwdcosy[i] = shwlist[i]->Direction().Y();
1435  shwdcosz[i] = shwlist[i]->Direction().Z();
1436  shwstartx[i] = shwlist[i]->ShowerStart().X();
1437  shwstarty[i] = shwlist[i]->ShowerStart().Y();
1438  shwstartz[i] = shwlist[i]->ShowerStart().Z();
1439  for (size_t j = 0; j<(shwlist[i]->Energy()).size(); ++j){
1440  shwenergy[i][j] = shwlist[i]->Energy()[j];
1441  }
1442  for (size_t j = 0; j<(shwlist[i]->dEdx()).size(); ++j){
1443  shwdedx[i][j] = shwlist[i]->dEdx()[j];
1444  }
1445  shwbestplane[i] = shwlist[i]->best_plane();
1446  if (fmsh.isValid()){
1447  auto vhit = fmsh.at(i);
1448  for (size_t h = 0; h < vhit.size(); ++h){
1449  if (vhit[h].key()<kMaxHits){
1450  hit_shwkey[vhit[h].key()] = shwlist[i].key();
1451  }
1452  }
1453  }
1454  if (!isdata&&fmsh.isValid()){
1455  // Find true track for each reconstructed track
1456  int TrackID = 0;
1457  std::vector< art::Ptr<recob::Hit> > allHits = fmsh.at(i);
1458  std::map<int,double> trkide;
1459  for(size_t h = 0; h < allHits.size(); ++h){
1460  art::Ptr<recob::Hit> hit = allHits[h];
1461  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
1462  for(size_t e = 0; e < TrackIDs.size(); ++e){
1463  trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
1464  }
1465  }
1466  // Work out which IDE despoited the most charge in the hit if there was more than one.
1467  double maxe = -1;
1468  double tote = 0;
1469  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
1470  tote += ii->second;
1471  if ((ii->second)>maxe){
1472  maxe = ii->second;
1473  TrackID = ii->first;
1474  }
1475  }
1476  // Now have trackID, so get PdG code and T0 etc.
1477  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
1478  if (particle){
1479  shwg4id[i] = TrackID;
1480  }
1481  }
1482  }
1483  }
1484 
1485  // flash information
1486  flash_total = flashlist.size();
1487  for ( int f = 0; f < std::min(flash_total,kMaxHits); ++f ) {
1488  flash_time[f] = flashlist[f]->Time();
1489  flash_width[f] = flashlist[f]->TimeWidth();
1490  flash_abstime[f] = flashlist[f]->AbsTime();
1491  flash_YCenter[f] = flashlist[f]->YCenter();
1492  flash_YWidth[f] = flashlist[f]->YWidth();
1493  flash_ZCenter[f] = flashlist[f]->ZCenter();
1494  flash_ZWidth[f] = flashlist[f]->ZWidth();
1495  flash_TotalPE[f] = flashlist[f]->TotalPE();
1496  }
1497 
1498  if (!isdata){
1499 
1500  // * MC truth information
1501  std::vector<art::Ptr<simb::MCTruth> > mclist;
1502  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
1503  if (mctruthListHandle)
1504  art::fill_ptr_vector(mclist, mctruthListHandle);
1505 
1506  std::vector<art::Ptr<simb::MCFlux> > fluxlist;
1507  auto mcfluxListHandle = evt.getHandle< std::vector<simb::MCFlux> >(fGenieGenModuleLabel);
1508  if (mcfluxListHandle)
1509  art::fill_ptr_vector(fluxlist, mcfluxListHandle);
1510 
1511 
1512  mcevts_truth=mclist.size();
1513  if (mcevts_truth){
1514  art::Ptr<simb::MCTruth> mctruth = mclist[0];
1515  if (mctruth->Origin() == simb::kBeamNeutrino){
1516  nuPDG_truth = mctruth->GetNeutrino().Nu().PdgCode();
1517  ccnc_truth = mctruth->GetNeutrino().CCNC();
1518  mode_truth = mctruth->GetNeutrino().Mode();
1519  Q2_truth = mctruth->GetNeutrino().QSqr();
1520  W_truth = mctruth->GetNeutrino().W();
1521  X_truth = mctruth->GetNeutrino().X();
1522  Y_truth = mctruth->GetNeutrino().Y();
1523  hitnuc_truth = mctruth->GetNeutrino().HitNuc();
1524  target_truth = mctruth->GetNeutrino().Target();
1525  enu_truth = mctruth->GetNeutrino().Nu().E();
1526  nuvtxx_truth = mctruth->GetNeutrino().Nu().Vx();
1527  nuvtxy_truth = mctruth->GetNeutrino().Nu().Vy();
1528  nuvtxz_truth = mctruth->GetNeutrino().Nu().Vz();
1529  if (mctruth->GetNeutrino().Nu().P()){
1530  nu_dcosx_truth = mctruth->GetNeutrino().Nu().Px()/mctruth->GetNeutrino().Nu().P();
1531  nu_dcosy_truth = mctruth->GetNeutrino().Nu().Py()/mctruth->GetNeutrino().Nu().P();
1532  nu_dcosz_truth = mctruth->GetNeutrino().Nu().Pz()/mctruth->GetNeutrino().Nu().P();
1533  }
1534  lep_mom_truth = mctruth->GetNeutrino().Lepton().P();
1535  if (mctruth->GetNeutrino().Lepton().P()){
1536  lep_dcosx_truth = mctruth->GetNeutrino().Lepton().Px()/mctruth->GetNeutrino().Lepton().P();
1537  lep_dcosy_truth = mctruth->GetNeutrino().Lepton().Py()/mctruth->GetNeutrino().Lepton().P();
1538  lep_dcosz_truth = mctruth->GetNeutrino().Lepton().Pz()/mctruth->GetNeutrino().Lepton().P();
1539  }
1540 
1541  if (mctruth->NParticles()){
1542  simb::MCParticle particle = mctruth->GetParticle(0);
1543  t0_truth = particle.T();
1544  }
1545 
1546 
1547  float mindist2 = 9999; // cm;
1548  TVector3 nuvtx(nuvtxx_truth, nuvtxy_truth, nuvtxz_truth);
1549  infidvol = insideFidVol(nuvtxx_truth, nuvtxy_truth, nuvtxz_truth);
1550  //find the closest reco vertex to the neutrino mc truth
1551  if (infidvol)
1552  {
1553  // vertex is when at least two tracks meet
1554  for(size_t i = 0; i < vtxlist.size(); ++i){ // loop over vertices
1555  Double_t xyz[3] = {};
1556  vtxlist[i]->XYZ(xyz);
1557  TVector3 vtxreco(xyz);
1558  float dist2 = pma::Dist2(vtxreco, nuvtx);
1559  if (dist2 < mindist2)
1560  {
1561  mindist2 = dist2;
1562  vtxrecomc = std::sqrt(dist2);
1563  vtxrecomcx = vtxreco.X() - nuvtxx_truth;
1564  vtxrecomcy = vtxreco.Y() - nuvtxy_truth;
1565  vtxrecomcz = vtxreco.Z() - nuvtxz_truth;
1566  }
1567  }
1568 
1569  // two endpoints of tracks are somehow also vertices...
1570  for (size_t i = 0; i < tracklist.size(); ++i){ // loop over tracks
1571  float dist2 = pma::Dist2(tracklist[i]->Vertex(), nuvtx);
1572  if (dist2 < mindist2)
1573  {
1574  mindist2 = dist2;
1575  vtxrecomc = std::sqrt(dist2);
1576  vtxrecomcx = tracklist[i]->Vertex().X() - nuvtxx_truth;
1577  vtxrecomcy = tracklist[i]->Vertex().Y() - nuvtxy_truth;
1578  vtxrecomcz = tracklist[i]->Vertex().Z() - nuvtxz_truth;
1579 
1580  }
1581  dist2 = pma::Dist2(tracklist[i]->End(), nuvtx);
1582  if (dist2 < mindist2)
1583  {
1584  mindist2 = dist2;
1585  vtxrecomc = std::sqrt(dist2);
1586  vtxrecomcx = tracklist[i]->End().X() - nuvtxx_truth;
1587  vtxrecomcy = tracklist[i]->End().Y() - nuvtxy_truth;
1588  vtxrecomcz = tracklist[i]->End().Z() - nuvtxz_truth;
1589 
1590  }
1591  }
1592  }
1593  }//is neutrino
1594  } // if numver of mctruth si not zero
1595 
1596  if (fluxlist.size()){
1597  ptype_flux = fluxlist[0]->fptype;
1598  pdpx_flux = fluxlist[0]->fpdpx;
1599  pdpy_flux = fluxlist[0]->fpdpy;
1600  pdpz_flux = fluxlist[0]->fpdpz;
1601  pntype_flux = fluxlist[0]->fntype;
1602  vx_flux = fluxlist[0]->fvx;
1603  vy_flux = fluxlist[0]->fvy;
1604  vz_flux = fluxlist[0]->fvz;
1605  }
1606 
1607  //save g4 particle information
1608  std::vector<const simb::MCParticle* > geant_part;
1609 
1610  // ### Looping over all the Geant4 particles from the BackTrackerService ###
1611  for(size_t p = 0; p < plist.size(); ++p)
1612  {
1613  // ### Filling the vector with MC Particles ###
1614  geant_part.push_back(plist.Particle(p));
1615  }
1616 
1617  //std::cout<<"No of geant part= "<<geant_part.size()<<std::endl;
1618 
1619  // ### Setting a string for primary ###
1620  std::string pri("primary");
1621 
1622  int primary=0;
1623  int geant_particle=0;
1624 
1625  // ############################################################
1626  // ### Determine the number of primary particles from geant ###
1627  // ############################################################
1628  for( unsigned int i = 0; i < geant_part.size(); ++i ){
1629  geant_particle++;
1630  // ### Counting the number of primary particles ###
1631  if(geant_part[i]->Process()==pri)
1632  { primary++;}
1633  }//<---End i loop
1634 
1635 
1636  // ### Saving the number of primary particles ###
1637  no_primaries=primary;
1638  // ### Saving the number of Geant4 particles ###
1639  geant_list_size=geant_particle;
1640 
1641  // ### Looping over all the Geant4 particles ###
1642  for( unsigned int i = 0; i < geant_part.size(); ++i ){
1643 
1644  // ### If this particle is primary, set = 1 ###
1645  if(geant_part[i]->Process()==pri)
1646  {process_primary[i]=1;}
1647  // ### If this particle is not-primary, set = 0 ###
1648  else
1649  {process_primary[i]=0;}
1650 
1651  // ### Saving the particles mother TrackID ###
1652  Mother[i]=geant_part[i]->Mother();
1653  // ### Saving the particles TrackID ###
1654  TrackId[i]=geant_part[i]->TrackId();
1655  // ### Saving the PDG Code ###
1656  pdg[i]=geant_part[i]->PdgCode();
1657  // ### Saving the particles Energy ###
1658  Eng[i]=geant_part[i]->E();
1659 
1660  // ### Saving the Px, Py, Pz info ###
1661  Px[i]=geant_part[i]->Px();
1662  Py[i]=geant_part[i]->Py();
1663  Pz[i]=geant_part[i]->Pz();
1664 
1665  // ### Saving the Start and End Point for this particle ###
1666  StartPointx[i]=geant_part[i]->Vx();
1667  StartPointy[i]=geant_part[i]->Vy();
1668  StartPointz[i]=geant_part[i]->Vz();
1669  EndPointx[i]=geant_part[i]->EndPosition()[0];
1670  EndPointy[i]=geant_part[i]->EndPosition()[1];
1671  EndPointz[i]=geant_part[i]->EndPosition()[2];
1672 
1673  // ### Saving the processes for this particle ###
1674  //std::cout<<"finding proc"<<std::endl;
1675  G4Process.push_back( geant_part[i]->Process() );
1676  G4FinalProcess.push_back( geant_part[i]->EndProcess() );
1677  //std::cout<<"found proc"<<std::endl;
1678  // std::cout << "ID " << TrackId[i] << ", pdg " << pdg[i] << ", Start X,Y,Z " << StartPointx[i] << ", " << StartPointy[i] << ", " << StartPointz[i]
1679  // << ", End XYZ " << EndPointx[i] << ", " << EndPointy[i] << ", " << EndPointz[i] << ", Start Proc " << G4Process[i] << ", End Proc " << G4FinalProcess[i]
1680  // << std::endl;
1681 
1682  // ### Saving the Start direction cosines for this particle ###
1683  Startdcosx[i] = geant_part[i]->Momentum(0).Px() / geant_part[i]->Momentum(0).P();
1684  Startdcosy[i] = geant_part[i]->Momentum(0).Py() / geant_part[i]->Momentum(0).P();
1685  Startdcosz[i] = geant_part[i]->Momentum(0).Pz() / geant_part[i]->Momentum(0).P();
1686  // ### Saving the number of Daughters for this particle ###
1687  NumberDaughters[i]=geant_part[i]->NumberDaughters();
1688 
1689  } //geant particles
1690 
1691 
1692  }//is neutrino
1693 
1694  if(fMakeAnaTree)
1695  fTree->Fill();
1696 
1697 }
1698 
1699 
1700 
1701 
1702 bool dunemva::MVAAlg::insideFidVol(const double posX, const double posY, const double posZ)
1703 {
1704 
1705  double vtx[3] = {posX, posY, posZ};
1706  bool inside = false;
1707 
1708  geo::TPCID idtpc = fGeom->FindTPCAtPosition(vtx);
1709 
1710  if (fGeom->HasTPC(idtpc))
1711  {
1712  const geo::TPCGeo& tpcgeo = fGeom->GetElement(idtpc);
1713  double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
1714  double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
1715  double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();
1716 
1717  for (size_t c = 0; c < fGeom->Ncryostats(); c++)
1718  {
1719  const geo::CryostatGeo& cryostat = fGeom->Cryostat(c);
1720  for (size_t t = 0; t < cryostat.NTPC(); t++)
1721  {
1722  const geo::TPCGeo& tpcg = cryostat.TPC(t);
1723  if (tpcg.MinX() < minx) minx = tpcg.MinX();
1724  if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX();
1725  if (tpcg.MinY() < miny) miny = tpcg.MinY();
1726  if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
1727  if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
1728  if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
1729  }
1730  }
1731 
1732 
1733  //x
1734  double dista = fabs(minx - posX);
1735  double distb = fabs(posX - maxx);
1736  if ((posX > minx) && (posX < maxx) &&
1737  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1738  //y
1739  dista = fabs(maxy - posY);
1740  distb = fabs(posY - miny);
1741  if (inside && (posY > miny) && (posY < maxy) &&
1742  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1743  else inside = false;
1744 
1745  //z
1746  dista = fabs(maxz - posZ);
1747  distb = fabs(posZ - minz);
1748  if (inside && (posZ > minz) && (posZ < maxz) &&
1749  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1750  else inside = false;
1751  }
1752 
1753  return inside;
1754 }
1755 
1756 void dunemva::MVAAlg::ResetVars(){
1757 
1758  G4Process.clear();
1759  G4FinalProcess.clear();
1760 
1761  run = -9999;
1762  subrun = -9999;
1763  event = -9999;
1764  evttime = -9999;
1765  taulife = 0;
1766  isdata = -9999;
1767 
1768  ntracks_reco = 0;
1769  for (int i = 0; i < kMaxTrack; ++i){
1770  trkid[i] = -9999;
1771  trkstartx[i] = -9999;
1772  trkstarty[i] = -9999;
1773  trkstartz[i] = -9999;
1774  trkendx[i] = -9999;
1775  trkendy[i] = -9999;
1776  trkendz[i] = -9999;
1777  trkstartdcosx[i] = -9999;
1778  trkstartdcosy[i] = -9999;
1779  trkstartdcosz[i] = -9999;
1780  trkenddcosx[i] = -9999;
1781  trkenddcosy[i] = -9999;
1782  trkenddcosz[i] = -9999;
1783  trklen[i] = -9999;
1784  trkbestplane[i] = -9999;
1785  for (int j = 0; j<3; ++j){
1786  trkke[i][j] = -9999;
1787  trkpida[i][j] = -9999;
1788  }
1789  trkg4id[i] = -9999;
1790  trkg4pdg[i] = -9999;
1791  trkg4startx[i] = -9999;
1792  trkg4starty[i] = -9999;
1793  trkg4startz[i] = -9999;
1794  trkg4initdedx[i] = -9999;
1795  }
1796 
1797  nshws = 0;
1798  for (int i = 0; i<kMaxShower; ++i){
1799  shwid[i] = -9999;
1800  shwdcosx[i] = -9999;
1801  shwdcosy[i] = -9999;
1802  shwdcosz[i] = -9999;
1803  shwstartx[i] = -9999;
1804  shwstarty[i] = -9999;
1805  shwstartz[i] = -9999;
1806  for (int j = 0; j<3; ++j){
1807  shwenergy[i][j] = -9999;
1808  shwdedx[i][j] = -9999;
1809  }
1810  shwbestplane[i] = -9999;
1811  trkg4id[i] = -9999;
1812  }
1813 
1814  flash_total = 0;
1815  for (int f = 0; f < kMaxFlash; ++f) {
1816  flash_time[f] = -9999;
1817  flash_width[f] = -9999;
1818  flash_abstime[f] = -9999;
1819  flash_YCenter[f] = -9999;
1820  flash_YWidth[f] = -9999;
1821  flash_ZCenter[f] = -9999;
1822  flash_ZWidth[f] = -9999;
1823  flash_TotalPE[f] = -9999;
1824  }
1825 
1826  nhits = 0;
1827  nhits_stored = 0;
1828  for (int i = 0; i<kMaxHits; ++i){
1829  hit_plane[i] = -9999;
1830  hit_wire[i] = -9999;
1831  hit_tpc[i] = -9999;
1832  hit_channel[i] = -9999;
1833  hit_peakT[i] = -9999;
1834  hit_charge[i] = -9999;
1835  hit_summedADC[i] = -9999;
1836  hit_startT[i] = -9999;
1837  hit_endT[i] = -9999;
1838  hit_trkkey[i] = -9999;
1839  hit_dQds[i] = -9999;
1840  hit_dEds[i] = -9999;
1841  hit_resrange[i] = -9999;
1842  hit_shwkey[i] = -9999;
1843  }
1844 
1845  infidvol = 0;
1846  nvtx = 0;
1847  for (int i = 0; i<kMaxVertices; ++i){
1848  vtx[i][0] = -9999;
1849  vtx[i][1] = -9999;
1850  vtx[i][2] = -9999;
1851  }
1852  vtxrecomc = 9999;
1853  vtxrecomcx = 9999;
1854  vtxrecomcy = 9999;
1855  vtxrecomcz = 9999;
1856 
1857  mcevts_truth = -9999;
1858  nuPDG_truth = -9999;
1859  ccnc_truth = -9999;
1860  mode_truth = -9999;
1861  enu_truth = -9999;
1862  Q2_truth = -9999;
1863  W_truth = -9999;
1864  X_truth = -9999;
1865  Y_truth = -9999;
1866  hitnuc_truth = -9999;
1867  target_truth = -9999;
1868  nuvtxx_truth = -9999;
1869  nuvtxy_truth = -9999;
1870  nuvtxz_truth = -9999;
1871  nu_dcosx_truth = -9999;
1872  nu_dcosy_truth = -9999;
1873  nu_dcosz_truth = -9999;
1874  lep_mom_truth = -9999;
1875  lep_dcosx_truth = -9999;
1876  lep_dcosy_truth = -9999;
1877  lep_dcosz_truth = -9999;
1878  t0_truth = -9999;
1879 
1880  no_primaries = -99999;
1881  geant_list_size=-9999;
1882  for (int i = 0; i<kMaxPrimaries; ++i){
1883  pdg[i] = -99999;
1884  Eng[i] = -99999;
1885  Px[i] = -99999;
1886  Py[i] = -99999;
1887  Pz[i] = -99999;
1888  StartPointx[i] = -99999;
1889  StartPointy[i] = -99999;
1890  StartPointz[i] = -99999;
1891  EndPointx[i] = -99999;
1892  EndPointy[i] = -99999;
1893  EndPointz[i] = -99999;
1894  Startdcosx[i]= -99999;
1895  Startdcosy[i]= -99999;
1896  Startdcosz[i]= -99999;
1897  NumberDaughters[i] = -99999;
1898  Mother[i] = -99999;
1899  TrackId[i] = -99999;
1900  process_primary[i] = -99999;
1901  }
1902 
1903  ptype_flux = -99999;
1904  pdpx_flux = -99999;
1905  pdpy_flux = -99999;
1906  pdpz_flux = -99999;
1907  pntype_flux = -99999;
1908  vx_flux = -99999;
1909  vy_flux = -99999;
1910  vz_flux = -99999;
1911 
1912  isinfidvol = 0;
1913  isinfidvoltruth = 0;
1914  oscpro = 0;
1915 }
1916 
1917 
1918 void dunemva::MVAAlg::MakeTree(){
1919 
1920  fTree = tfs->make<TTree>("nueana","analysis tree");
1921  fTree->Branch("run",&run,"run/I");
1922  fTree->Branch("subrun",&subrun,"subrun/I");
1923  fTree->Branch("event",&event,"event/I");
1924  fTree->Branch("evttime",&evttime,"evttime/F");
1925  fTree->Branch("taulife",&taulife,"taulife/F");
1926  fTree->Branch("isdata",&isdata,"isdata/S");
1927  fTree->Branch("ntracks_reco",&ntracks_reco,"ntracks_reco/I");
1928  fTree->Branch("trkid",trkid,"trkid[ntracks_reco]/I");
1929  fTree->Branch("trkstartx",trkstartx,"trkstartx[ntracks_reco]/F");
1930  fTree->Branch("trkstarty",trkstarty,"trkstarty[ntracks_reco]/F");
1931  fTree->Branch("trkstartz",trkstartz,"trkstartz[ntracks_reco]/F");
1932  fTree->Branch("trkendx",trkendx,"trkendx[ntracks_reco]/F");
1933  fTree->Branch("trkendy",trkendy,"trkendy[ntracks_reco]/F");
1934  fTree->Branch("trkendz",trkendz,"trkendz[ntracks_reco]/F");
1935  fTree->Branch("trkstartdcosx",trkstartdcosx,"trkstartdcosx[ntracks_reco]/F");
1936  fTree->Branch("trkstartdcosy",trkstartdcosy,"trkstartdcosy[ntracks_reco]/F");
1937  fTree->Branch("trkstartdcosz",trkstartdcosz,"trkstartdcosz[ntracks_reco]/F");
1938  fTree->Branch("trkenddcosx",trkenddcosx,"trkenddcosx[ntracks_reco]/F");
1939  fTree->Branch("trkenddcosy",trkenddcosy,"trkenddcosy[ntracks_reco]/F");
1940  fTree->Branch("trkenddcosz",trkenddcosz,"trkenddcosz[ntracks_reco]/F");
1941  fTree->Branch("trklen",trklen,"trklen[ntracks_reco]/F");
1942  fTree->Branch("trkbestplane",trkbestplane,"trkbestplane[ntracks_reco]/I");
1943  fTree->Branch("trkke",trkke,"trkke[ntracks_reco][3]/F");
1944  fTree->Branch("trkpida",trkpida,"trkpida[ntracks_reco][3]/F");
1945  fTree->Branch("trkg4id",trkg4id,"trkg4id[ntracks_reco]/I");
1946  fTree->Branch("trkg4pdg",trkg4pdg,"trkg4pdg[ntracks_reco]/I");
1947  fTree->Branch("trkg4startx",trkg4startx,"trkg4startx[ntracks_reco]/F");
1948  fTree->Branch("trkg4starty",trkg4starty,"trkg4starty[ntracks_reco]/F");
1949  fTree->Branch("trkg4startz",trkg4startz,"trkg4startz[ntracks_reco]/F");
1950  fTree->Branch("trkg4initdedx",trkg4initdedx,"trkg4initdedx[ntracks_reco]/F");
1951  fTree->Branch("nshws",&nshws,"nshws/I");
1952  fTree->Branch("shwid",shwid,"shwid[nshws]/I");
1953  fTree->Branch("shwdcosx",shwdcosx,"shwdcosx[nshws]/F");
1954  fTree->Branch("shwdcosy",shwdcosy,"shwdcosy[nshws]/F");
1955  fTree->Branch("shwdcosz",shwdcosz,"shwdcosz[nshws]/F");
1956  fTree->Branch("shwstartx",shwstartx,"shwstartx[nshws]/F");
1957  fTree->Branch("shwstarty",shwstarty,"shwstarty[nshws]/F");
1958  fTree->Branch("shwstartz",shwstartz,"shwstartz[nshws]/F");
1959  fTree->Branch("shwenergy",shwenergy,"shwenergy[nshws][3]/F");
1960  fTree->Branch("shwdedx",shwdedx,"shwdedx[nshws][3]/F");
1961  fTree->Branch("shwbestplane",shwbestplane,"shwbestplane[nshws]/I");
1962  fTree->Branch("shwg4id",shwg4id,"shwg4id[nshws]/I");
1963  fTree->Branch("flash_total" ,&flash_total ,"flash_total/I");
1964  fTree->Branch("flash_time" ,flash_time ,"flash_time[flash_total]/F");
1965  fTree->Branch("flash_width" ,flash_width ,"flash_width[flash_total]/F");
1966  fTree->Branch("flash_abstime",flash_abstime,"flash_abstime[flash_total]/F");
1967  fTree->Branch("flash_YCenter",flash_YCenter,"flash_YCenter[flash_total]/F");
1968  fTree->Branch("flash_YWidth" ,flash_YWidth ,"flash_YWidth[flash_total]/F");
1969  fTree->Branch("flash_ZCenter",flash_ZCenter,"flash_ZCenter[flash_total]/F");
1970  fTree->Branch("flash_ZWidth" ,flash_ZWidth ,"flash_ZWidth[flash_total]/F");
1971  fTree->Branch("flash_TotalPE",flash_TotalPE,"flash_TotalPE[flash_total]/F");
1972  fTree->Branch("nhits",&nhits,"nhits/I");
1973  fTree->Branch("nhits_stored",&nhits_stored,"nhits_stored/I");
1974  fTree->Branch("hit_plane",hit_plane,"hit_plane[nhits_stored]/S");
1975  fTree->Branch("hit_tpc",hit_tpc,"hit_tpc[nhits_stored]/S");
1976  fTree->Branch("hit_wire",hit_wire,"hit_wire[nhits_stored]/S");
1977  fTree->Branch("hit_channel",hit_channel,"hit_channel[nhits_stored]/I");
1978  fTree->Branch("hit_peakT",hit_peakT,"hit_peakT[nhits_stored]/F");
1979  fTree->Branch("hit_charge",hit_charge,"hit_charge[nhits_stored]/F");
1980  fTree->Branch("hit_summedADC",hit_summedADC,"hit_summedADC[nhits_stored]/F");
1981  fTree->Branch("hit_startT",hit_startT,"hit_startT[nhits_stored]/F");
1982  fTree->Branch("hit_endT",hit_endT,"hit_endT[nhits_stored]/F");
1983  fTree->Branch("hit_trkkey",hit_trkkey,"hit_trkkey[nhits_stored]/I");
1984  fTree->Branch("hit_dQds",hit_dQds,"hit_dQds[nhits_stored]/F");
1985  fTree->Branch("hit_dEds",hit_dEds,"hit_dEds[nhits_stored]/F");
1986  fTree->Branch("hit_resrange",hit_resrange,"hit_resrange[nhits_stored]/F");
1987  fTree->Branch("hit_shwkey",hit_shwkey,"hit_shwkey[nhits_stored]/I");
1988  fTree->Branch("infidvol",&infidvol,"infidvol/I");
1989  fTree->Branch("nvtx",&nvtx,"nvtx/S");
1990  fTree->Branch("vtx",vtx,"vtx[nvtx][3]/F");
1991  fTree->Branch("vtxrecomc",&vtxrecomc,"vtxrecomc/F");
1992  fTree->Branch("vtxrecomcx",&vtxrecomcx,"vtxrecomcx/F");
1993  fTree->Branch("vtxrecomcy",&vtxrecomcy,"vtxrecomcy/F");
1994  fTree->Branch("vtxrecomcz",&vtxrecomcz,"vtxrecomcz/F");
1995  fTree->Branch("mcevts_truth",&mcevts_truth,"mcevts_truth/I");
1996  fTree->Branch("nuPDG_truth",&nuPDG_truth,"nuPDG_truth/I");
1997  fTree->Branch("ccnc_truth",&ccnc_truth,"ccnc_truth/I");
1998  fTree->Branch("mode_truth",&mode_truth,"mode_truth/I");
1999  fTree->Branch("enu_truth",&enu_truth,"enu_truth/F");
2000  fTree->Branch("Q2_truth",&Q2_truth,"Q2_truth/F");
2001  fTree->Branch("W_truth",&W_truth,"W_truth/F");
2002  fTree->Branch("X_truth",&X_truth,"X_truth/F");
2003  fTree->Branch("Y_truth",&Y_truth,"Y_truth/F");
2004  fTree->Branch("hitnuc_truth",&hitnuc_truth,"hitnuc_truth/I");
2005  fTree->Branch("target_truth",&target_truth,"target_truth/I");
2006  fTree->Branch("nuvtxx_truth",&nuvtxx_truth,"nuvtxx_truth/F");
2007  fTree->Branch("nuvtxy_truth",&nuvtxy_truth,"nuvtxy_truth/F");
2008  fTree->Branch("nuvtxz_truth",&nuvtxz_truth,"nuvtxz_truth/F");
2009  fTree->Branch("nu_dcosx_truth",&nu_dcosx_truth,"nu_dcosx_truth/F");
2010  fTree->Branch("nu_dcosy_truth",&nu_dcosy_truth,"nu_dcosy_truth/F");
2011  fTree->Branch("nu_dcosz_truth",&nu_dcosz_truth,"nu_dcosz_truth/F");
2012  fTree->Branch("lep_mom_truth",&lep_mom_truth,"lep_mom_truth/F");
2013  fTree->Branch("lep_dcosx_truth",&lep_dcosx_truth,"lep_dcosx_truth/F");
2014  fTree->Branch("lep_dcosy_truth",&lep_dcosy_truth,"lep_dcosy_truth/F");
2015  fTree->Branch("lep_dcosz_truth",&lep_dcosz_truth,"lep_dcosz_truth/F");
2016  fTree->Branch("t0_truth",&t0_truth,"t0_truth/F");
2017  fTree->Branch("no_primaries",&no_primaries,"no_primaries/I");
2018  fTree->Branch("geant_list_size",&geant_list_size,"geant_list_size/I");
2019  fTree->Branch("pdg",pdg,"pdg[geant_list_size]/I");
2020  fTree->Branch("Eng",Eng,"Eng[geant_list_size]/F");
2021  fTree->Branch("Px",Px,"Px[geant_list_size]/F");
2022  fTree->Branch("Py",Py,"Py[geant_list_size]/F");
2023  fTree->Branch("Pz",Pz,"Pz[geant_list_size]/F");
2024  fTree->Branch("StartPointx",StartPointx,"StartPointx[geant_list_size]/F");
2025  fTree->Branch("StartPointy",StartPointy,"StartPointy[geant_list_size]/F");
2026  fTree->Branch("StartPointz",StartPointz,"StartPointz[geant_list_size]/F");
2027  fTree->Branch("EndPointx",EndPointx,"EndPointx[geant_list_size]/F");
2028  fTree->Branch("EndPointy",EndPointy,"EndPointy[geant_list_size]/F");
2029  fTree->Branch("EndPointz",EndPointz,"EndPointz[geant_list_size]/F");
2030  fTree->Branch("Startdcosx",Startdcosx,"Startdcosx[geant_list_size]/F");
2031  fTree->Branch("Startdcosy",Startdcosy,"Startdcosy[geant_list_size]/F");
2032  fTree->Branch("Startdcosz",Startdcosz,"Startdcosz[geant_list_size]/F");
2033  fTree->Branch("N(((((((((umberDaughters",NumberDaughters,"NumberDaughters[geant_list_size]/I");
2034  fTree->Branch("Mother",Mother,"Mother[geant_list_size]/I");
2035  fTree->Branch("TrackId",TrackId,"TrackId[geant_list_size]/I");
2036  fTree->Branch("process_primary",process_primary,"process_primary[geant_list_size]/I");
2037  fTree->Branch("G4Process",&G4Process);//,"G4Process[geant_list_size]");
2038  fTree->Branch("G4FinalProcess",&G4FinalProcess);//,"G4FinalProcess[geant_list_size]");
2039  fTree->Branch("ptype_flux",&ptype_flux,"ptype_flux/I");
2040  fTree->Branch("pdpx_flux",&pdpx_flux,"pdpx_flux/F");
2041  fTree->Branch("pdpy_flux",&pdpy_flux,"pdpy_flux/F");
2042  fTree->Branch("pdpz_flux",&pdpz_flux,"pdpz_flux/F");
2043  fTree->Branch("pntype_flux",&pntype_flux,"pntype_flux/I");
2044  fTree->Branch("vx_flux",&vx_flux,"vx_flux/F");
2045  fTree->Branch("vy_flux",&vy_flux,"vy_flux/F");
2046  fTree->Branch("vz_flux",&vz_flux,"vz_flux/F");
2047 
2048 }
name
static QCString name
Definition: declinfo.cpp:673
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
dunemva::MVAAlg::fClusterModuleLabel
std::string fClusterModuleLabel
Definition: MVAAlg.h:314
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
dunemva::MVAAlg::shwdcosy
Float_t shwdcosy[kMaxShower]
Definition: MVAAlg.h:229
dunemva::MVAAlg::maxtrklength
float maxtrklength
Definition: MVAAlg.h:99
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
dunemva::MVAAlg::oscpro
float oscpro
Definition: MVAAlg.h:327
dunemva::MVAAlg::shwg4id
int shwg4id[kMaxTrack]
Definition: MVAAlg.h:237
dunemva::MVAAlg::lep_mom_truth
Float_t lep_mom_truth
Definition: MVAAlg.h:268
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:148
dunemva::MVAAlg::htrkdedx
TH1D * htrkdedx[5]
Definition: MVAAlg.h:366
dunemva::MVAAlg::hvtxx
TH1D * hvtxx[5]
Definition: MVAAlg.h:362
art::Timestamp::timeLow
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
dunemva::MVAAlg::fFile
std::ofstream fFile
Definition: MVAAlg.h:138
dunemva::MVAAlg::flash_YWidth
Float_t flash_YWidth[kMaxFlash]
Definition: MVAAlg.h:245
dunemva::MVAAlg::wirecharge
double wirecharge
Definition: MVAAlg.h:130
simb::MCNeutrino::QSqr
double QSqr() const
Definition: MCNeutrino.h:157
dunemva::MVAAlg::trkg4startx
float trkg4startx[kMaxTrack]
Definition: MVAAlg.h:193
dunemva::MVAAlg::evttime
float evttime
Definition: MVAAlg.h:165
dunemva::MVAAlg::hit_startT
Float_t hit_startT[kMaxHits]
Definition: MVAAlg.h:207
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
dunemva::MVAAlg::fTree
TTree * fTree
Definition: MVAAlg.h:156
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
dunemva::MVAAlg::fWireModuleLabel
std::string fWireModuleLabel
Definition: MVAAlg.h:312
dunemva::MVAAlg::fPOTModuleLabel
std::string fPOTModuleLabel
Definition: MVAAlg.h:319
geo::GeometryCore::GetElement
CryostatGeo const & GetElement(geo::CryostatID const &cryoid) const
Definition: GeometryCore.h:2033
dunemva::MVAAlg::hit_channel
Int_t hit_channel[kMaxHits]
Definition: MVAAlg.h:202
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
dunemva::MVAAlg::trkg4id
int trkg4id[kMaxTrack]
Definition: MVAAlg.h:191
dunemva::MVAAlg::W_truth
Float_t W_truth
Definition: MVAAlg.h:257
dunemva::MVAAlg::trkstartz
float trkstartz[kMaxTrack]
Definition: MVAAlg.h:175
dunemva::MVAAlg::PrepareEvent
void PrepareEvent(const art::Event &event)
Definition: MVAAlg.cxx:1125
dunemva::MVAAlg::mva_numu_nutau
TH1D * mva_numu_nutau
Definition: MVAAlg.h:344
dunemva::MVAAlg::mva_nue_beam
TH1D * mva_nue_beam
Definition: MVAAlg.h:337
result
static QCString result
Definition: fortranscanner.cpp:56614
dunemva::MVAAlg::pdpz_flux
Float_t pdpz_flux
Definition: MVAAlg.h:302
dunemva::MVAAlg::hdisy
TH1D * hdisy[5]
Definition: MVAAlg.h:374
dunemva::MVAAlg::Q2_truth
Float_t Q2_truth
Definition: MVAAlg.h:256
dunemva::MVAAlg::fract_5_wires
float fract_5_wires
Definition: MVAAlg.h:112
dunemva::MVAAlg::EndPointz
Float_t EndPointz[kMaxPrimaries]
Definition: MVAAlg.h:287
cheat::BackTrackerService::HitToTrackIDEs
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Definition: BackTrackerService_service.cc:173
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
dunemva::MVAAlg::enu_numu_numu
TH1D * enu_numu_numu
Definition: MVAAlg.h:350
dunemva::MVAAlg::mva_nc
TH1D * mva_nc
Definition: MVAAlg.h:335
dunemva::MVAAlg::hnhitspershw
TH1D * hnhitspershw[5]
Definition: MVAAlg.h:377
dunemva::MVAAlg::EndPointx
Float_t EndPointx[kMaxPrimaries]
Definition: MVAAlg.h:285
dunemva::MVAAlg::mcevts_truth
Int_t mcevts_truth
Definition: MVAAlg.h:251
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
string
std::string string
Definition: nybbler.cc:12
TrackIDs
std::vector< TrackID > TrackIDs
Definition: ClusteringValidation_module.cc:59
dunemva::MVAAlg::hpida
TH1D * hpida[5]
Definition: MVAAlg.h:369
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService_service.cc:128
dunemva::MVAAlg::trkendy
float trkendy[kMaxTrack]
Definition: MVAAlg.h:177
dunemva::MVAAlg::nuvtxx_truth
Float_t nuvtxx_truth
Definition: MVAAlg.h:262
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
dunemva::MVAAlg::showerdedx
float showerdedx
Definition: MVAAlg.h:106
dunemva::MVAAlg::trkstartx
float trkstartx[kMaxTrack]
Definition: MVAAlg.h:173
dunemva::MVAAlg::flash_abstime
Float_t flash_abstime[kMaxFlash]
Definition: MVAAlg.h:243
dunemva::MVAAlg::hit_dQds
Float_t hit_dQds[kMaxHits]
Definition: MVAAlg.h:210
dunemva::MVAAlg::shwstartx
Float_t shwstartx[kMaxShower]
Definition: MVAAlg.h:231
dunemva::MVAAlg::enu_nc
TH1D * enu_nc
Definition: MVAAlg.h:347
dunemva::MVAAlg::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: MVAAlg.cxx:273
dunemva::MVAAlg::nhits_stored
int nhits_stored
Definition: MVAAlg.h:199
dunemva::MVAAlg::shwbestplane
int shwbestplane[kMaxShower]
Definition: MVAAlg.h:236
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
dunemva::MVAAlg::nvtx
Short_t nvtx
Definition: MVAAlg.h:217
dunemva::MVAAlg::Startdcosz
Float_t Startdcosz[kMaxPrimaries]
Definition: MVAAlg.h:290
dunemva::MVAAlg::shwmax
float shwmax
Definition: MVAAlg.h:111
dunemva::MVAAlg::fFidVolCut
double fFidVolCut
Definition: MVAAlg.h:323
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
simb::MCTruth::Origin
simb::Origin_t Origin() const
Definition: MCTruth.h:74
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
geo::BoxBoundedGeo::MinX
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
dunemva::MVAAlg::fOscPro
float fOscPro
Definition: MVAAlg.h:94
dunemva::MVAAlg::fBeamMode
std::string fBeamMode
Definition: MVAAlg.h:148
dunemva::MVAAlg::nuvtxy_truth
Float_t nuvtxy_truth
Definition: MVAAlg.h:263
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
dunemva::MVAAlg::Py
Float_t Py[kMaxPrimaries]
Definition: MVAAlg.h:280
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
dunemva::MVAAlg::rawcharge
double rawcharge
Definition: MVAAlg.h:129
simb::MCNeutrino::HitNuc
int HitNuc() const
Definition: MCNeutrino.h:152
kMaxPrimaries
constexpr int kMaxPrimaries
Definition: MVAAlg.h:64
kMaxTrack
constexpr int kMaxTrack
Definition: MVAAlg.h:60
dunemva::MVAAlg::StartPointx
Float_t StartPointx[kMaxPrimaries]
Definition: MVAAlg.h:282
dunemva::MVAAlg::TrackId
int TrackId[kMaxPrimaries]
Definition: MVAAlg.h:292
dunemva::MVAAlg::~MVAAlg
virtual ~MVAAlg()
Definition: MVAAlg.cxx:268
dunemva::MVAAlg::vtxrecomc
Float_t vtxrecomc
Definition: MVAAlg.h:220
dunemva::MVAAlg::MakeTree
void MakeTree()
Definition: MVAAlg.cxx:1918
dunemva::MVAAlg::Norm
float Norm(int ccnc, int nu0, int nu1, int subrun)
Definition: MVAAlg.cxx:994
lar::sparse_vector::get_ranges
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
Definition: sparse_vector.h:717
dunemva::MVAAlg::trkpida_save
float trkpida_save
Definition: MVAAlg.h:104
dunemva::MVAAlg::G4Process
std::vector< std::string > G4Process
Definition: MVAAlg.h:295
geo::kZ
Planes which measure Z direction.
Definition: geo_types.h:132
dunemva::MVAAlg::enu
TH1D * enu[5]
Definition: MVAAlg.h:360
dunemva::MVAAlg::lep_dcosx_truth
Float_t lep_dcosx_truth
Definition: MVAAlg.h:269
dunemva::MVAAlg::hit_summedADC
Float_t hit_summedADC[kMaxHits]
Definition: MVAAlg.h:206
dunemva::MVAAlg::trkg4starty
float trkg4starty[kMaxTrack]
Definition: MVAAlg.h:194
geo::BoxBoundedGeo::MaxX
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
dunemva::MVAAlg::flash_YCenter
Float_t flash_YCenter[kMaxFlash]
Definition: MVAAlg.h:244
dunemva::MVAAlg::ntracks_reco
int ntracks_reco
Definition: MVAAlg.h:171
dunemva::MVAAlg::shwdcosz
Float_t shwdcosz[kMaxShower]
Definition: MVAAlg.h:230
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
dunemva::MVAAlg::trkg4pdg
int trkg4pdg[kMaxTrack]
Definition: MVAAlg.h:192
dunemva::MVAAlg::hshwdedx
TH1D * hshwdedx[5]
Definition: MVAAlg.h:372
dunemva::MVAAlg::trkenddcosz
float trkenddcosz[kMaxTrack]
Definition: MVAAlg.h:184
geo::GeometryCore::SignalType
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
Definition: GeometryCore.cxx:578
Index
unsigned int Index
Definition: sspmapmaker_v1.c:14
dunemva::MVAAlg::flash_TotalPE
Float_t flash_TotalPE[kMaxFlash]
Definition: MVAAlg.h:248
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
dunemva::MVAAlg::trkke
float trkke[kMaxTrack][3]
Definition: MVAAlg.h:186
dunemva::MVAAlg::mva_nue_osc
TH1D * mva_nue_osc
Definition: MVAAlg.h:334
dunemva::MVAAlg::htrklen
TH1D * htrklen[5]
Definition: MVAAlg.h:365
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:1994
dunemva::MVAAlg::hdisx
TH1D * hdisx[5]
Definition: MVAAlg.h:373
dunemva::MVAAlg::fMakeAnaTree
bool fMakeAnaTree
Definition: MVAAlg.h:331
dunemva::MVAAlg::fWeightFile
std::string fWeightFile
Definition: MVAAlg.h:142
dunemva::MVAAlg::trkdedx
float trkdedx
Definition: MVAAlg.h:100
dunemva::MVAAlg::hit_peakT
Float_t hit_peakT[kMaxHits]
Definition: MVAAlg.h:204
dunemva::MVAAlg::tfs
art::ServiceHandle< art::TFileService > tfs
Definition: MVAAlg.h:145
dunemva::MVAAlg::trkg4initdedx
float trkg4initdedx[kMaxTrack]
Definition: MVAAlg.h:196
dunemva::MVAAlg::fract_10_wires
float fract_10_wires
Definition: MVAAlg.h:113
dunemva::MVAAlg::Run
void Run(const art::Event &evt, std::vector< double > &result, double &wgt)
geo::GeometryCore::Nwires
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
Definition: GeometryCore.h:3342
dunemva::MVAAlg::nuPDG_truth
Int_t nuPDG_truth
Definition: MVAAlg.h:252
Geometry.h
art framework interface to geometry description
dunemva::MVAAlg::trkcosy
float trkcosy
Definition: MVAAlg.h:124
dunemva::MVAAlg::MVAAlg
MVAAlg(fhicl::ParameterSet const &p)
Definition: MVAAlg.cxx:23
dunemva::MVAAlg::hfr100w
TH1D * hfr100w[5]
Definition: MVAAlg.h:379
dunemva::MVAAlg::fVertexModuleLabel
std::string fVertexModuleLabel
Definition: MVAAlg.h:317
dunemva::MVAAlg::lep_dcosz_truth
Float_t lep_dcosz_truth
Definition: MVAAlg.h:271
kMaxFlash
constexpr int kMaxFlash
Definition: MVAAlg.h:65
raw::TDCtick_t
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
dunemva::MVAAlg::fReader
TMVA::Reader fReader
Definition: MVAAlg.h:140
dunemva::MVAAlg::Pz
Float_t Pz[kMaxPrimaries]
Definition: MVAAlg.h:281
recob::Track::Vector_t
tracking::Vector_t Vector_t
Definition: Track.h:54
geo::GeometryCore::FindTPCAtPosition
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.h:2454
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
dunemva::MVAAlg::shwdisz
float shwdisz
Definition: MVAAlg.h:119
dunemva::MVAAlg::pntype_flux
Int_t pntype_flux
Definition: MVAAlg.h:303
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
dunemva::MVAAlg::shwcosx
float shwcosx
Definition: MVAAlg.h:120
dunemva::MVAAlg::trkstartdcosy
float trkstartdcosy[kMaxTrack]
Definition: MVAAlg.h:180
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
dunemva::MVAAlg::trkenddcosx
float trkenddcosx[kMaxTrack]
Definition: MVAAlg.h:182
dunemva::MVAAlg::hrt
TH1D * hrt[5]
Definition: MVAAlg.h:368
dunemva::MVAAlg::shwstartz
Float_t shwstartz[kMaxShower]
Definition: MVAAlg.h:233
dunemva::MVAAlg::fRawDigitModuleLabel
std::string fRawDigitModuleLabel
Definition: MVAAlg.h:311
dunemva::MVAAlg::flash_ZCenter
Float_t flash_ZCenter[kMaxFlash]
Definition: MVAAlg.h:246
e
const double e
Definition: gUpMuFluxGen.cxx:165
dunemva::MVAAlg::fCalorimetryModuleLabel
std::string fCalorimetryModuleLabel
Definition: MVAAlg.h:321
simb::MCNeutrino::Lepton
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
simb::MCNeutrino::W
double W() const
Definition: MCNeutrino.h:154
dunemva::MVAAlg::flash_time
Float_t flash_time[kMaxFlash]
Definition: MVAAlg.h:241
dunemva::MVAAlg::mva_nue_nutau
TH1D * mva_nue_nutau
Definition: MVAAlg.h:345
dunemva::MVAAlg::trkstartdcosx
float trkstartdcosx[kMaxTrack]
Definition: MVAAlg.h:179
art::DataViewImpl::time
Timestamp time() const
Definition: DataViewImpl.cc:106
dunemva::MVAAlg::fMakeWeightTree
bool fMakeWeightTree
Definition: MVAAlg.h:332
dunemva::MVAAlg::hit_wire
Short_t hit_wire[kMaxHits]
Definition: MVAAlg.h:201
dunemva::MVAAlg::fCalorimetryAlg
calo::CalorimetryAlg fCalorimetryAlg
Definition: MVAAlg.h:329
dunemva::MVAAlg::shwstarty
Float_t shwstarty[kMaxShower]
Definition: MVAAlg.h:232
simb::MCNeutrino::Y
double Y() const
Definition: MCNeutrino.h:156
dunemva::MVAAlg::fHitsModuleLabel
std::string fHitsModuleLabel
Definition: MVAAlg.h:313
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
dunemva::MVAAlg::fract_100_wires
float fract_100_wires
Definition: MVAAlg.h:115
dunemva::MVAAlg::pot
double pot
Definition: MVAAlg.h:168
dunemva::MVAAlg::fWeightTree
TTree * fWeightTree[6]
Definition: MVAAlg.h:355
dunemva::MVAAlg::hrch
TH1D * hrch[5]
Definition: MVAAlg.h:367
dunemva::MVAAlg::hit_dEds
Float_t hit_dEds[kMaxHits]
Definition: MVAAlg.h:211
dunemva::MVAAlg::pdpx_flux
Float_t pdpx_flux
Definition: MVAAlg.h:300
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
dunemva::MVAAlg::trkcosz
float trkcosz
Definition: MVAAlg.h:125
dunemva::MVAAlg::X_truth
Float_t X_truth
Definition: MVAAlg.h:258
dunemva::MVAAlg::EndPointy
Float_t EndPointy[kMaxPrimaries]
Definition: MVAAlg.h:286
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
art::Ptr::key
key_type key() const noexcept
Definition: Ptr.h:216
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
dunemva::MVAAlg::G4FinalProcess
std::vector< std::string > G4FinalProcess
Definition: MVAAlg.h:296
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
geo::BoxBoundedGeo::MinZ
double MinZ() const
Returns the world z coordinate of the start of the box.
Definition: BoxBoundedGeo.h:118
dunemva::MVAAlg::hitnuc_truth
Int_t hitnuc_truth
Definition: MVAAlg.h:260
dunemva::MVAAlg::nuvtxz_truth
Float_t nuvtxz_truth
Definition: MVAAlg.h:264
dunemva::MVAAlg::frshower
float frshower
Definition: MVAAlg.h:108
dunemva::MVAAlg::isdata
short isdata
Definition: MVAAlg.h:167
geo::CryostatGeo::NTPC
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
dunemva::MVAAlg::trkendz
float trkendz[kMaxTrack]
Definition: MVAAlg.h:178
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:321
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
dunemva::MVAAlg::hit_endT
Float_t hit_endT[kMaxHits]
Definition: MVAAlg.h:208
dunemva::MVAAlg::StartPointy
Float_t StartPointy[kMaxPrimaries]
Definition: MVAAlg.h:283
dunemva::MVAAlg::mva_nue_nue
TH1D * mva_nue_nue
Definition: MVAAlg.h:341
simb::MCNeutrino::X
double X() const
Definition: MCNeutrino.h:155
dunemva::MVAAlg::vz_flux
Float_t vz_flux
Definition: MVAAlg.h:306
dunemva::MVAAlg::mva_nue_numu
TH1D * mva_nue_numu
Definition: MVAAlg.h:343
dunemva::MVAAlg::shwdisy
float shwdisy
Definition: MVAAlg.h:118
dunemva::MVAAlg::ResetVars
void ResetVars()
Definition: MVAAlg.cxx:1756
dunemva::MVAAlg::hdisz
TH1D * hdisz[5]
Definition: MVAAlg.h:375
dunemva::MVAAlg::isinfidvoltruth
int isinfidvoltruth
Definition: MVAAlg.h:326
TrackID
TrackID
Definition: ClusteringValidation_module.cc:57
dunemva::MVAAlg::hit_shwkey
Int_t hit_shwkey[kMaxHits]
Definition: MVAAlg.h:213
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
dunemva::MVAAlg::mode_truth
Int_t mode_truth
Definition: MVAAlg.h:254
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
dunemva::MVAAlg::hit_plane
Short_t hit_plane[kMaxHits]
Definition: MVAAlg.h:200
dunemva::MVAAlg::process_primary
int process_primary[kMaxPrimaries]
Definition: MVAAlg.h:294
dunemva::MVAAlg::lep_dcosy_truth
Float_t lep_dcosy_truth
Definition: MVAAlg.h:270
dunemva::MVAAlg::flash_width
Float_t flash_width[kMaxFlash]
Definition: MVAAlg.h:242
dunemva::MVAAlg::events_truth
float events_truth[5]
Definition: MVAAlg.h:358
dunemva::MVAAlg::hfrshower
TH1D * hfrshower[5]
Definition: MVAAlg.h:376
dunemva::MVAAlg::hit_charge
Float_t hit_charge[kMaxHits]
Definition: MVAAlg.h:205
dunemva::MVAAlg::mva_nutau
TH1D * mva_nutau
Definition: MVAAlg.h:338
dunemva::MVAAlg::fGeom
art::ServiceHandle< geo::Geometry > fGeom
Definition: MVAAlg.h:144
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:1202
dunemva::MVAAlg::trkid
int trkid[kMaxTrack]
Definition: MVAAlg.h:172
dunemva::MVAAlg::trklen
float trklen[kMaxTrack]
Definition: MVAAlg.h:185
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
dunemva::MVAAlg::enu_nue_nue
TH1D * enu_nue_nue
Definition: MVAAlg.h:349
dunemva::MVAAlg::trkstarty
float trkstarty[kMaxTrack]
Definition: MVAAlg.h:174
dunemva::MVAAlg::ccnc_truth
Int_t ccnc_truth
Definition: MVAAlg.h:253
geo::BoxBoundedGeo::MaxY
double MaxY() const
Returns the world y coordinate of the end of the box.
Definition: BoxBoundedGeo.h:106
cheat::ParticleInventoryService::ParticleList
const sim::ParticleList & ParticleList() const
Definition: ParticleInventoryService_service.cc:114
dunemva::MVAAlg::shwcosy
float shwcosy
Definition: MVAAlg.h:121
geo::GeometryCore::NTPC
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:2303
simb::MCTruth::GetParticle
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
dunemva::MVAAlg::fShowerModuleLabel
std::string fShowerModuleLabel
Definition: MVAAlg.h:316
dunemva::MVAAlg::OscPro
float OscPro(int ccnc, int nu0, int nu1, float NuE)
Definition: MVAAlg.cxx:1041
test.p
p
Definition: test.py:223
geo::GeometryCore::HasTPC
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
Definition: GeometryCore.h:2397
dunemva::MVAAlg::enu_nue_numu
TH1D * enu_nue_numu
Definition: MVAAlg.h:351
dunemva::MVAAlg::heshw
TH1D * heshw[5]
Definition: MVAAlg.h:371
kMaxVertices
constexpr int kMaxVertices
Definition: AnalysisTree_module.cc:107
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
dunemva::MVAAlg::trkbestplane
int trkbestplane[kMaxTrack]
Definition: MVAAlg.h:188
kMaxHits
constexpr int kMaxHits
Definition: AnalysisTree_module.cc:104
lar::sparse_vector< float >
dunemva::MVAAlg::avgtrklength
float avgtrklength
Definition: MVAAlg.h:98
simb::MCNeutrino::Target
int Target() const
Definition: MCNeutrino.h:151
dunemva::MVAAlg::geant_list_size
int geant_list_size
Definition: MVAAlg.h:276
AssociationUtil.h
Utility object to perform functions of association.
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
dunemva::MVAAlg::hvtxz
TH1D * hvtxz[5]
Definition: MVAAlg.h:364
dunemva::MVAAlg::vtxrecomcz
Float_t vtxrecomcz
Definition: MVAAlg.h:223
dunemva::MVAAlg::shwdcosx
Float_t shwdcosx[kMaxShower]
Definition: MVAAlg.h:228
dunemva::MVAAlg::insideFidVol
bool insideFidVol(const double posX, const double posY, const double posZ)
Definition: MVAAlg.cxx:1702
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
End
void End(void)
Definition: gXSecComp.cxx:210
dunemva::MVAAlg::fSelect
std::string fSelect
Definition: MVAAlg.h:147
dunemva::MVAAlg::fFlashModuleLabel
std::string fFlashModuleLabel
Definition: MVAAlg.h:320
dunemva::MVAAlg::weight
float weight
Definition: MVAAlg.h:95
dunemva::MVAAlg::target_truth
Int_t target_truth
Definition: MVAAlg.h:261
dunemva::MVAAlg::shwdisx
float shwdisx
Definition: MVAAlg.h:117
geo::BoxBoundedGeo::MaxZ
double MaxZ() const
Returns the world z coordinate of the end of the box.
Definition: BoxBoundedGeo.h:121
dunemva::MVAAlg::StartPointz
Float_t StartPointz[kMaxPrimaries]
Definition: MVAAlg.h:284
dunemva::MVAAlg::events_reco
float events_reco[5]
Definition: MVAAlg.h:359
dunemva::MVAAlg::evtcharge
float evtcharge
Definition: MVAAlg.h:96
dunemva::MVAAlg::fract_50_wires
float fract_50_wires
Definition: MVAAlg.h:114
dunemva::MVAAlg::no_primaries
int no_primaries
Definition: MVAAlg.h:275
dunemva::MVAAlg::eshower
float eshower
Definition: MVAAlg.h:107
recob::Track::Point_t
tracking::Point_t Point_t
Definition: Track.h:53
dunemva::MVAAlg::mva_numu
TH1D * mva_numu
Definition: MVAAlg.h:336
lar::util::TrackPitchInView
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
Definition: TrackUtils.cxx:78
dunemva::MVAAlg::vtx
Float_t vtx[kMaxVertices][3]
Definition: MVAAlg.h:218
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
calo::CalorimetryAlg::dEdx_AREA
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
Definition: CalorimetryAlg.cxx:131
dunemva::MVAAlg::trkrch
float trkrch
Definition: MVAAlg.h:101
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
dunemva::MVAAlg::taulife
float taulife
Definition: MVAAlg.h:166
dunemva::MVAAlg::trkfr
float trkfr
Definition: MVAAlg.h:103
dunemva::MVAAlg::Mother
int Mother[kMaxPrimaries]
Definition: MVAAlg.h:293
dunemva::MVAAlg::vtxrecomcy
Float_t vtxrecomcy
Definition: MVAAlg.h:222
dunemva::MVAAlg::trkcosx
float trkcosx
Definition: MVAAlg.h:123
dunemva::MVAAlg::fMVAMethod
std::string fMVAMethod
Definition: MVAAlg.h:141
dunemva::MVAAlg::hevtcharge
TH1D * hevtcharge[5]
Definition: MVAAlg.h:378
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
dunemva::MVAAlg::vtxrecomcx
Float_t vtxrecomcx
Definition: MVAAlg.h:221
dunemva::MVAAlg::hit_tpc
Short_t hit_tpc[kMaxHits]
Definition: MVAAlg.h:203
dunemva::MVAAlg::enu_numu_nue
TH1D * enu_numu_nue
Definition: MVAAlg.h:348
train.x
list x
Definition: train.py:276
dunemva::MVAAlg::fPOT
TTree * fPOT
Definition: MVAAlg.h:157
art::detail::Vertex
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
dunemva::MVAAlg::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: MVAAlg.h:315
dunemva::MVAAlg::trkendx
float trkendx[kMaxTrack]
Definition: MVAAlg.h:176
dunemva::MVAAlg::t0_truth
Float_t t0_truth
Definition: MVAAlg.h:272
dunemva::MVAAlg::hnshw
TH1D * hnshw[5]
Definition: MVAAlg.h:370
dunemva::MVAAlg::trkg4startz
float trkg4startz[kMaxTrack]
Definition: MVAAlg.h:195
dunemva::MVAAlg::shwdedx
Float_t shwdedx[kMaxShower][3]
Definition: MVAAlg.h:235
dunemva::MVAAlg::Startdcosx
Float_t Startdcosx[kMaxPrimaries]
Definition: MVAAlg.h:288
dunemva::MVAAlg::hit_trkkey
Int_t hit_trkkey[kMaxHits]
Definition: MVAAlg.h:209
dunemva::MVAAlg::fGenieGenModuleLabel
std::string fGenieGenModuleLabel
Definition: MVAAlg.h:318
dunemva::MVAAlg::shwenergy
Float_t shwenergy[kMaxShower][3]
Definition: MVAAlg.h:234
dunemva::MVAAlg::hvtxy
TH1D * hvtxy[5]
Definition: MVAAlg.h:363
dunemva::MVAAlg::vx_flux
Float_t vx_flux
Definition: MVAAlg.h:304
dunemva::MVAAlg::hit_resrange
Float_t hit_resrange[kMaxHits]
Definition: MVAAlg.h:212
dunemva::MVAAlg::nhitspershw
float nhitspershw
Definition: MVAAlg.h:109
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
dunemva::MVAAlg::Startdcosy
Float_t Startdcosy[kMaxPrimaries]
Definition: MVAAlg.h:289
dunemva::MVAAlg::shwdis
float shwdis
Definition: MVAAlg.h:116
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
dunemva::MVAAlg::enu_truth
Float_t enu_truth
Definition: MVAAlg.h:255
calo::CalorimetryAlg::LifetimeCorrection
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Definition: CalorimetryAlg.cxx:246
dunemva::MVAAlg::pdg
int pdg[kMaxPrimaries]
Definition: MVAAlg.h:277
dunemva::MVAAlg::Y_truth
Float_t Y_truth
Definition: MVAAlg.h:259
dunemva::MVAAlg::nu_dcosx_truth
Float_t nu_dcosx_truth
Definition: MVAAlg.h:265
dunemva::MVAAlg::trkenddcosy
float trkenddcosy[kMaxTrack]
Definition: MVAAlg.h:183
dunemva::MVAAlg::vy_flux
Float_t vy_flux
Definition: MVAAlg.h:305
dunemva::MVAAlg::mva_numu_numu
TH1D * mva_numu_numu
Definition: MVAAlg.h:342
dunemva::MVAAlg::ntrack
float ntrack
Definition: MVAAlg.h:97
genie::units::sr
static constexpr double sr
Definition: Units.h:166
dunemva::MVAAlg::enu_nue_nutau
TH1D * enu_nue_nutau
Definition: MVAAlg.h:353
dunemva::MVAAlg::shwid
int shwid[kMaxShower]
Definition: MVAAlg.h:227
geo::BoxBoundedGeo::MinY
double MinY() const
Returns the world y coordinate of the start of the box.
Definition: BoxBoundedGeo.h:103
dunemva::MVAAlg::NumberDaughters
int NumberDaughters[kMaxPrimaries]
Definition: MVAAlg.h:291
simb::MCNeutrino::Mode
int Mode() const
Definition: MCNeutrino.h:149
bool
int bool
Definition: qglobal.h:345
dunemva::MVAAlg::trkpida
float trkpida[kMaxTrack][3]
Definition: MVAAlg.h:187
TrackUtils.h
Utility functions to extract information from recob::Track
dunemva::MVAAlg::endSubRun
void endSubRun(const art::SubRun &sr)
Definition: MVAAlg.cxx:373
dunemva::MVAAlg::Px
Float_t Px[kMaxPrimaries]
Definition: MVAAlg.h:279
art::Event::id
EventID id() const
Definition: Event.cc:34
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:222
dunemva::MVAAlg::enu_numu_nutau
TH1D * enu_numu_nutau
Definition: MVAAlg.h:352
cet::search_path
Definition: search_path.h:36
dunemva::MVAAlg::shwlength
float shwlength
Definition: MVAAlg.h:110
calo
calorimetry
Definition: CaloChecker_module.cc:43
dunemva::MVAAlg::nshower
float nshower
Definition: MVAAlg.h:105
dunemva::MVAAlg::flash_total
int flash_total
Definition: MVAAlg.h:240
dunemva::MVAAlg::shwcosz
float shwcosz
Definition: MVAAlg.h:122
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
event
Event finding and building.
Definition: EventCheater_module.cc:32
dunemva::MVAAlg::trkrt
float trkrt
Definition: MVAAlg.h:102
dunemva::MVAAlg::flash_ZWidth
Float_t flash_ZWidth[kMaxFlash]
Definition: MVAAlg.h:247
TPCGeo.h
Encapsulate the construction of a single detector plane.
kMaxShower
constexpr int kMaxShower
Definition: MVAAlg.h:61
dunemva::MVAAlg::trkstartdcosz
float trkstartdcosz[kMaxTrack]
Definition: MVAAlg.h:181
dunemva::MVAAlg::enu_osc
TH1D * enu_osc[5]
Definition: MVAAlg.h:361
dunemva::MVAAlg::nu_dcosy_truth
Float_t nu_dcosy_truth
Definition: MVAAlg.h:266
geo::kCollection
Signal from collection planes.
Definition: geo_types.h:146
dunemva::MVAAlg::Eng
Float_t Eng[kMaxPrimaries]
Definition: MVAAlg.h:278
simb::kBeamNeutrino
Beam neutrinos.
Definition: MCTruth.h:23
dunemva::MVAAlg::nu_dcosz_truth
Float_t nu_dcosz_truth
Definition: MVAAlg.h:267
dunemva::MVAAlg::pdpy_flux
Float_t pdpy_flux
Definition: MVAAlg.h:301
dunemva::MVAAlg::CalculateInputs
void CalculateInputs()
Definition: MVAAlg.cxx:384
dunemva::MVAAlg::ptype_flux
Int_t ptype_flux
Definition: MVAAlg.h:299
dunemva::MVAAlg::isinfidvol
int isinfidvol
Definition: MVAAlg.h:325
dunemva::MVAAlg::mva_numu_nue
TH1D * mva_numu_nue
Definition: MVAAlg.h:340
Generated by   doxygen 1.8.11