BreakPointProtonAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: BreakPointProtonAna
3 // Module Type: analyzer
4 // File: BreakPointProtonAna_module.cc
5 //
6 // Generated at Wed Sep 24 07:35:37 2014 by Hayes Merritt using artmod
7 // from cetpkgsupport v1_07_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <cmath>
11 #include <iostream>
12 #include <algorithm>
13 
14 // Histogram includes
15 #include "TH1.h"
16 #include "TH2.h"
17 #include "TTree.h"
18 #include "TVector.h"
19 #include "TProfile.h"
20 #include "TMath.h"
21 
22 // ART includes
23 #include "Calibrator/Calibrator.h"
37 #include "fhiclcpp/ParameterSet.h"
39 
40 #include "RecoBase/Cluster.h"
41 #include "RecoBase/Track.h"
42 #include "RecoBase/Vertex.h"
43 #include "RecoBase/FitSum.h"
44 #include "RecoBase/RecoHit.h"
45 #include "Utilities/AssociationUtil.h"
46 #include "Utilities/func/MathUtil.h"
47 
48 #include "MCCheater/BackTracker.h"
52 
54 #include "SummaryData/POTSum.h"
55 #include "SummaryData/SpillData.h"
56 
57 namespace bpfit {
58  class BreakPointProtonAna;
59 }
60 
62 public:
63  explicit BreakPointProtonAna(fhicl::ParameterSet const & p);
64  // The destructor generated by the compiler is fine for classes
65  // without bare pointers or other resource use.
66 
67  // Plugins should not be copied or assigned.
68  BreakPointProtonAna(BreakPointProtonAna const &) = delete;
72 
73  // Required functions.
74  void analyze(art::Event const & e) override;
75  void beginJob() override;
76  void endJob() override;
77  virtual void endSubRun(const art::SubRun &sr) override;
78 
79 private:
80 
81  // get BackTracker handle
83  // Grab calibrator
85 
86  // Declare member data here.
88  std::string fVertexLabel; ///< Where to find vertices
89  std::string fProngLabel; ///< Where to find prongs
90  std::string fProngInstance; ///< Which type of prong
91  std::string fTrackLabel; ///< Where to find tracks
92  std::string fFitSumLabel; ///< Where to find FitSum objects
93  std::string fMCTruthLabel; ///< Where to find MCTruth information
94  std::string fGeneratorLabel; ///< Where to find spill information in MC
95  std::string fSpillLabel; ///< Where to find spill information in data
96  bool fIsMC; ///< Is MC or data
97  double fTotalPOT; ///< Total POT for one subrun
98 
99  //TTree *fOutTree;
100 
101  ////////////////////
102  //// Histograms ////
103  ////////////////////
104  TH1D *fSRPOT;
105  TH1F *fTrueNuE;
106  TH1F *fRecoNuE;
107  TH1F *fRecoNuP;
108  TH1F *fPrDotXPr;
110  TH1F *fPrjxnDiff;
112  TH1F *fMuonPurity;
115  TH1F *fNuMomentum;
132  //TH1F *fShowerWidth;
171  //TH1F *fMuDelMomFracDiff;
172  TH1F *fPrCosTheta;
173  //TH1F *fMuCosTheta;
175  //unsigned int Run = -1;
176  //unsigned int SubRun = -1;
177  //unsigned int Event = -1;
178  //unsigned int Slice = -1;
179 };
180 
181 
183  :
184  EDAnalyzer(p),
185  fSlicerToken(consumes<std::vector<rb::Cluster>>(p.get< std::string >("SlicerLabel"))),
186  fVertexLabel (p.get< std::string >("VertexLabel" )),
187  fProngLabel (p.get< std::string >("ProngLabel" )),
188  fProngInstance (p.get< std::string >("ProngInstance" )),
189  fTrackLabel (p.get< std::string >("TrackLabel" )),
190  fFitSumLabel (p.get< std::string >("FitSumLabel" )),
191  fMCTruthLabel (p.get< std::string >("MCTruthLabel" )),
192  fGeneratorLabel(p.get< std::string >("GeneratorLabel")),
193  fSpillLabel (p.get< std::string >("SpillLabel" )),
194  fIsMC (p.get< bool >("IsMC" ))
195  // More initializers here.
196 {}
197 
199 {
201 
202  fTotalPOT = 0.0;
203 
204  fSRPOT = tfs->make<TH1D>("fSRPOT","Total POT For One SubRun",5,0,5);
205  fTrueNuE = tfs->make<TH1F>("fTrueNuE","True Neutrino Energy",120,0,6);
206  fRecoNuE = tfs->make<TH1F>("fRecoNuE","Reconstructed Neutrino Energy",120,0,6);
207  fRecoNuP = tfs->make<TH1F>("fRecoNuP","Reconstructed Neutrino Momentum",120,0,6);
208  fPrDotXPr = tfs->make<TH1F>("fPrDotXPr","Expected Proton Momentum Dot Product With Reco Proton Momentum",100,-1,1);
209  fPrDotXPrSlcE = tfs->make<TH1F>("fPrDotXPrSlcE","Expected Proton Momentum Dot Product With Reco Proton Momentum Using Slice Energy",100,-1,1);
210  fPrMomFracBias = tfs->make<TH1F>("fMomDiff","(Expected Proton Momentum - Reco Proton Momentum)/Expected",100,-1,1);
211  fMuonMomentum = tfs->make<TH1F>("fMuonMomentum","Reconstructed Muon Momentum",120,0,6);
212  fProtonMomentum = tfs->make<TH1F>("fProtonMomentum","Reconstructed Proton Momentum",120,0,6);
213  fMuonPurity = tfs->make<TH1F>("fMuonPurity","Selected Muon Purity By Truth ID",6000,-3000,3000);
214  fProtonPurity = tfs->make<TH1F>("fProtonPurity","Selected Proton Purity By Truth ID",6000,-3000,3000);
215  fNuMomentum = tfs->make<TH1F>("fNuMomentum","True Neutrino Momentum",120,0,6);
216  fRecoNuETrue = tfs->make<TH1F>("fRecoNuETrue","Reconstructed Neutrino Energy",120,0,6);
217  fRecoNuPTrue = tfs->make<TH1F>("fRecoNuPTrue","Reconstructed Neutrino Momentum",120,0,6);
218  fPrDotXPrTrue = tfs->make<TH1F>("fPrDotXPrTrue","Expected Proton Momentum Dot Product With Reco Proton Momentum",100,-1,1);
219  fPrMomFracBiasTrue = tfs->make<TH1F>("fMomDiffTrue","(Expected Proton Momentum - Reco Proton Momentum)/Expected",100,-1,1);
220  fMuonMomentumTrue = tfs->make<TH1F>("fMuonMomentumTrue","Reconstructed Muon Momentum",120,0,6);
221  fProtonMomentumTrue = tfs->make<TH1F>("fProtonMomentumTrue","Reconstructed Proton Momentum",120,0,6);
222  fNuMomentumTrue = tfs->make<TH1F>("fNuMomentumTrue","True Neutrino Momentum",120,0,6);
223  fRecoNuETrueQE = tfs->make<TH1F>("fRecoNuETrueQE","Reconstructed Neutrino Energy",120,0,6);
224  fRecoNuPTrueQE = tfs->make<TH1F>("fRecoNuPTrueQE","Reconstructed Neutrino Momentum",120,0,6);
225  fPrDotXPrTrueQE = tfs->make<TH1F>("fPrDotXPrTrueQE","Expected Proton Momentum Dot Product With Reco Proton Momentum",100,-1,1);
226  fPrDotXPrTrueQESlcE = tfs->make<TH1F>("fPrDotXPrTrueQESlcE","Expected Proton Momentum Dot Product With Reco Proton Momentum Using Slice Energy",100,-1,1);
227  fPrMomFracBiasTrueQE = tfs->make<TH1F>("fMomDiffTrueQE","(Expected Proton Momentum - Reco Proton Momentum)/Expected",100,-1,1);
228  fMuonMomentumTrueQE = tfs->make<TH1F>("fMuonMomentumTrueQE","Reconstructed Muon Momentum",120,0,6);
229  fProtonMomentumTrueQE = tfs->make<TH1F>("fProtonMomentumTrueQE","Reconstructed Proton Momentum",120,0,6);
230  fNuMomentumTrueQE = tfs->make<TH1F>("fNuMomentumTrueQE","True Neutrino Momentum",120,0,6);
231  fPrjxnVsRecoNuE = tfs->make<TH2D>("fPrjxnVsRecoNuE","E^{reco}_{#nu} Vs Projection",120,0,6,100,-1,1);
232  fPmuFracDiff = tfs->make<TH1F>("fPmuFracDiff","Fractional Momentum Difference",400,-2,2);
233  fPprFracDiff = tfs->make<TH1F>("fPprFracDiff","Fractional Momentum Difference",400,-2,2);
234  fPmuFracDiffTrue = tfs->make<TH1F>("fPmuFracDiffTrue","Fractional Momentum Difference",400,-2,2);
235  fPprFracDiffTrue = tfs->make<TH1F>("fPprFracDiffTrue","Fractional Momentum Difference",400,-2,2);
236  fPmuFracDiffTrueQE = tfs->make<TH1F>("fPmuFracDiffTrueQE","Fractional Momentum Difference",400,-2,2);
237  fPprFracDiffTrueQE = tfs->make<TH1F>("fPprFracDiffTrueQE","Fractional Momentum Difference",400,-2,2);
238  fPhotonAsMuonMotherType = tfs->make<TH1F>("fPhotonAsMuonMotherType","Photon Mother Type",6000,-3000,3000);
239  fPhotonAsProtonMotherType = tfs->make<TH1F>("fPhotonAsProtonMotherType","Photon Mother Type",6000,-3000,3000);
240  fMuonTrkLngth = tfs->make<TH1F>("fMuonTrkLngth","Muon Track Length",1000,0,6000);
241  fProtonTrkLngth = tfs->make<TH1F>("fProtonTrkLngth","Proton Track Length",1000,0,6000);
242  fMuonTrkLngthGood = tfs->make<TH1F>("fMuonTrkLngthGood","Muon Track Length",1000,0,6000);
243  fProtonTrkLngthGood = tfs->make<TH1F>("fProtonTrkLngthGood","Proton Track Length",1000,0,6000);
244  fEvtQEPurity1 = tfs->make<TH1F>("fEvtQEPurity1","Event Interaction Type",10000,0,10000);
245  fEvtQEPurity2 = tfs->make<TH1F>("fEvtQEPurity2","Event Interaction Type",10000,0,10000);
246  fEvtQEPurity3 = tfs->make<TH1F>("fEvtQEPurity3","Event Interaction Type",10000,0,10000);
247  fEvtQEPurity4 = tfs->make<TH1F>("fEvtQEPurity4","Event Interaction Type",10000,0,10000);
248  fEvtQEPurity5 = tfs->make<TH1F>("fEvtQEPurity5","Event Interaction Type",10000,0,10000);
249  fEvtQEPurity6 = tfs->make<TH1F>("fEvtQEPurity6","Event Interaction Type",10000,0,10000);
250  fEvtQEPurity7 = tfs->make<TH1F>("fEvtQEPurity7","Event Interaction Type",10000,0,10000);
251  fEvtQEPurity8 = tfs->make<TH1F>("fEvtQEPurity8","Event Interaction Type",10000,0,10000);
252  fEvtQEPurity9 = tfs->make<TH1F>("fEvtQEPurity9","Event Interaction Type",10000,0,10000);
253  fEvtQEPurity10 = tfs->make<TH1F>("fEvtQEPurity10","Event Interaction Type",10000,0,10000);
254  fEvtQEPurity1slcE = tfs->make<TH1F>("fEvtQEPurity1slcE","Event Interaction Type",10000,0,10000);
255  fEvtQEPurity2slcE = tfs->make<TH1F>("fEvtQEPurity2slcE","Event Interaction Type",10000,0,10000);
256  fEvtQEPurity3slcE = tfs->make<TH1F>("fEvtQEPurity3slcE","Event Interaction Type",10000,0,10000);
257  fEvtQEPurity4slcE = tfs->make<TH1F>("fEvtQEPurity4slcE","Event Interaction Type",10000,0,10000);
258  fEvtQEPurity5slcE = tfs->make<TH1F>("fEvtQEPurity5slcE","Event Interaction Type",10000,0,10000);
259  fEvtQEPurity6slcE = tfs->make<TH1F>("fEvtQEPurity6slcE","Event Interaction Type",10000,0,10000);
260  fEvtQEPurity7slcE = tfs->make<TH1F>("fEvtQEPurity7slcE","Event Interaction Type",10000,0,10000);
261  fEvtQEPurity8slcE = tfs->make<TH1F>("fEvtQEPurity8slcE","Event Interaction Type",10000,0,10000);
262  fEvtQEPurity9slcE = tfs->make<TH1F>("fEvtQEPurity9slcE","Event Interaction Type",10000,0,10000);
263  fEvtQEPurity10slcE = tfs->make<TH1F>("fEvtQEPurity10slcE","Event Interaction Type",10000,0,10000);
264  fMuonPurityGood = tfs->make<TH1F>("fMuonPurityGood","Selected Muon Purity By Truth ID",6000,-3000,3000);
265  fProtonPurityGood = tfs->make<TH1F>("fProtonPurityGood","Selected Proton Purity By Truth ID",6000,-3000,3000);
266  fMuonPurityGoodSlcE = tfs->make<TH1F>("fMuonPurityGoodSlcE","Selected Muon Purity By Truth ID",6000,-3000,3000);
267  fProtonPurityGoodSlcE = tfs->make<TH1F>("fProtonPurityGoodSlcE","Selected Proton Purity By Truth ID",6000,-3000,3000);
268  fPrDelMomFracDiff = tfs->make<TH1F>("fPrDelDiff","(True Proton Momentum - Reco Proton Momentum)/True",100,-1,1);
269  // fMuDelMomFracDiff = tfs->make<TH1F>("fMomDiffTrue","(Expected Proton Momentum - Reco Proton Momentum)/Expected",100,-1,1);
270  fPrCosTheta = tfs->make<TH1F>("fPrCosTheta","Proton Direction Cos(#theta)",200,-1,1);
271  fMuVPrCosTheta = tfs->make<TH2F>("fMuVPrCosTheta","Muon Cos(#theta) Vs. Proton Cos(#theta)",200,-1,1,200,-1,1);
272  fPrjxnDiff = tfs->make<TH1F>("fPrjxnDiff","Difference In Proton Projection (slice energy vs true neutrino energy)",200,-1,1);
273 }
274 
276 {
278 
279  if(fIsMC) sr.getByLabel("generator", pot);
280  else sr.getByLabel("ifdbspillinfo",pot);
281 
282  if(!pot.failedToGet()){
283  if(!fIsMC) fTotalPOT += (pot->totgoodpot*1.0e12);
284  else fTotalPOT += pot->totgoodpot;
285  }
286 }
287 
289 {
290  std::cout<<" \n\n\n Total POT: "<<fTotalPOT<<std::endl;
291  // Fill histogram containing POT for file/subrun
292  fSRPOT->Fill(3,fTotalPOT);
293 }
294 
296 {
297  ///////////////////////////////////////////////////////////////////////
298  /// For reference purposes InteractionType 1001(1002) is CCQE(NCQE).///
299  ///////////////////////////////////////////////////////////////////////
300 
301  // get slices
303  e.getByToken(fSlicerToken, hslice);
304 
305  // get associations between slices and tracks
306  art::FindManyP< rb::Track > fSlTrk(hslice, e, fTrackLabel);
307 
308  // get associations between slices and vertices.
309  art::FindManyP<rb::Vertex> fSlVrtx(hslice, e, fVertexLabel);
310 
311  // get associations between slices and 3D prongs
313  fSlPrng(hslice, e, art::InputTag(fProngLabel,fProngInstance));
314 
315  // Get slices
317  for(unsigned int i=0;i<hslice->size();i++){
318  slices.push_back(art::Ptr<rb::Cluster>(hslice,i));
319  }
320 
321  // Get truth neutrino stuff
322  std::vector< std::vector< cheat::NeutrinoEffPur > > NuEffPur;
323  if(fIsMC) NuEffPur = bt->SlicesToMCTruthsTable(slices);
324 
325  // Loop over number of slices
326  for(unsigned int islice=0; islice<slices.size();islice++){
327  // get vertices
328  std::vector < art::Ptr< rb::Vertex > > vtx;
329  if(fSlVrtx.isValid()) vtx = fSlVrtx.at(islice);
330  bool isCntdVtx = false;
331  art::Ptr<rb::Vertex> vtx_0;
332  // by design have a max of 1 vertex per slice
333  if(vtx.size()==1){
334  vtx_0 = vtx.at(0);
335  isCntdVtx = fabs(vtx_0->GetX())<150
336  && fabs(vtx_0->GetY())<150
337  && vtx_0->GetZ()>50;
338  }
339 
340  // get prongs
341  std::vector < art::Ptr< rb::Prong > > prongs;
342  if(fSlPrng.isValid()) prongs = fSlPrng.at(islice);
343  // get associations between fuzzyk 3d prongs and tracks
344  art::FindManyP< rb::Track > fPrngTrk(prongs, e, fTrackLabel);
345 
346  TLorentzVector Muon, Proton, TruMuon;
347  bool haveMuon = false;
348  bool haveProton = false;
349  bool isTruMuon = false;
350  bool isTruProton = false;
351  int proton_pdg = 0;
352  int muon_pdg = 0;
353  float proton_track_length = -10.;
354  float muon_track_length = -10.;
355  float true_proton_mom = -10.;
356 
357  // Make sure the vertex is contained
358  if(!isCntdVtx) continue;
359 
360  // Calculate calibrated slice energy
361  art::PtrVector<rb::CellHit> slchits = slices.at(islice)->AllCells();
362  double slice_GeV = 1.0E-3;
363  // Make reco hit for every cell hit in order to get calibrated cell energy.
364  for(unsigned int j=0;j<slchits.size();j++){
365  const rb::CellHit *cellhit = slices.at(islice)->Cell(j).get();
366  double Wtemp = slices.at(islice)->W(cellhit);
367  rb::RecoHit rhit(hCal->MakeRecoHit(*cellhit,Wtemp));
368  if(rhit.IsCalibrated()) slice_GeV += rhit.GeV();
369  }
370 
371  for(unsigned int iprong=0; iprong<prongs.size();iprong++){
372  // Grab the current prong to get truth information
373  // art::Ptr< rb::Prong > prongi = prongs.at(iprong);
374  // get tracks associated with prongs
375  std::vector < art::Ptr< rb::Track > > pr_bp_trks;
376  if(fPrngTrk.isValid()){pr_bp_trks = fPrngTrk.at(iprong);}
377  else std::cout<<" Prong Track isn't valid!!!!"<<std::endl;
378  art::FindManyP< rb::FitSum > fKin(pr_bp_trks, e, fFitSumLabel);
379  unsigned int index = -1;
380  unsigned int lngst_trk_index = -1;
381 
382  for(unsigned int itrk=0; itrk<pr_bp_trks.size();itrk++){
383  art::Ptr< rb::Track > tracki = pr_bp_trks.at(itrk);
384  // Get the kinematic information for each track
385  std::vector<art::Ptr<rb::FitSum>> ikinsum;
386  if(fKin.isValid()) ikinsum = fKin.at(itrk);
387 
388  art::Ptr< rb::FitSum > ifit = ikinsum.at(0);
389 
390  for(unsigned int jtrk=itrk+1; jtrk<pr_bp_trks.size();jtrk++){
391  art::Ptr< rb::Track > trackj = pr_bp_trks.at(jtrk);
392  // Get the kinematic information for each track
393  std::vector<art::Ptr<rb::FitSum>> jkinsum;
394  if(fKin.isValid()) jkinsum = fKin.at(jtrk);
395  // Grab kinematics of compared track
396  art::Ptr< rb::FitSum > jfit = jkinsum.at(0);
397 
398  if(itrk<1){
399  // Determine which particle assumption is best fit
400  float long_trk = tracki->TotalLength();
401  float cur_trk_len = trackj->TotalLength();
402  float min_X2 = (ifit->Chi2(2)+ifit->Chi2(3))/(ifit->Ndof(2)+ifit->Ndof(3)-4);
403  //+(ifit->Chi2(0)+ifit->Chi2(1))/(ifit->Ndof(0)+ifit->Ndof(1)-4);
404  float cur_X2 = (jfit->Chi2(2)+jfit->Chi2(3))/(jfit->Ndof(2)+jfit->Ndof(3)-4);
405  //+(jfit->Chi2(0)+jfit->Chi2(1))/(jfit->Ndof(0)+jfit->Ndof(1)-4);
406  if(cur_X2 < min_X2){
407  min_X2 = cur_X2;
408  index = jtrk;
409  }
410  else index = itrk;
411  // Store index of longest track.
412  if(cur_trk_len>long_trk){
413  long_trk = cur_trk_len;
414  lngst_trk_index = jtrk;
415  }
416  else lngst_trk_index = itrk;
417  }
418  } // end loop over jtrk
419 
420  if(fIsMC){
421  if(itrk==index){
422  const std::vector< const sim::Particle *> trupart =
423  bt->HitsToParticle((*tracki).AllCells());
424  if(ifit->PDG()==13){
425  fMuonTrkLngth->Fill(tracki->TotalLength());
426  if(fabs(tracki->Stop().X())<150 && fabs(tracki->Stop().Y())<150
427  && tracki->Stop().Z()<1200
428  && ((ifit->Chi2(2)+ifit->Chi2(3))/(ifit->Ndof(2)+ifit->Ndof(3)-4)) < 0.1
429  && itrk==lngst_trk_index){
430  Muon.SetPxPyPzE(ifit->FourMom().Px(), ifit->FourMom().Py(), ifit->FourMom().Pz(), ifit->FourMom().E());
431  fMuonPurity->Fill(trupart.at(0)->PdgCode());
432  fPmuFracDiff->Fill((Muon.P()-trupart.at(0)->P())/trupart.at(0)->P());
433  haveMuon = true;
434  muon_pdg = trupart.at(0)->PdgCode();
435  muon_track_length = tracki->TotalLength();
436 
437  if(fabs(trupart.at(0)->PdgCode())==13){
438  fPmuFracDiffTrue->Fill((Muon.P()-trupart.at(0)->P())/trupart.at(0)->P());
439  isTruMuon = true;
440  }
441  if(fabs(trupart.at(0)->PdgCode())==22){
442  if(trupart.at(0)->Mother()!=0){
443  const sim::Particle *motherpart = bt->TrackIDToParticle(trupart.at(0)->Mother());
444  fPhotonAsMuonMotherType->Fill(motherpart->PdgCode());
445  }
446  }
447 
448  break; // now that we have a muon we want to skip to the next prong and associated BP tracks.
449  }
450  } // end if
451  else if(ifit->PDG()==2212){
452  fProtonTrkLngth->Fill(tracki->TotalLength());
453  if(tracki->TotalLength() < 125.0
454  && ((ifit->Chi2(2)+ifit->Chi2(3))/(ifit->Ndof(2)+ifit->Ndof(3)-4)) < 0.02){
455  fProtonPurity->Fill(trupart.at(0)->PdgCode());
456  fPprFracDiff->Fill((Muon.P()-trupart.at(0)->P())/trupart.at(0)->P());
457  Proton.SetPxPyPzE(ifit->FourMom().Px(), ifit->FourMom().Py(), ifit->FourMom().Pz(), ifit->FourMom().E());
458  haveProton = true;
459  proton_pdg = trupart.at(0)->PdgCode();
460  proton_track_length = tracki->TotalLength();
461  true_proton_mom = trupart.at(0)->P();
462  if(fabs(trupart.at(0)->PdgCode())==2212){
463  fPprFracDiffTrue->Fill((Proton.P()-trupart.at(0)->P())/trupart.at(0)->P());
464  isTruProton = true;
465  }
466  if(fabs(trupart.at(0)->PdgCode())==22){
467  if(trupart.at(0)->Mother()!=0){
468  const sim::Particle *motherpart = bt->TrackIDToParticle(trupart.at(0)->Mother());
469  fPhotonAsProtonMotherType->Fill(motherpart->PdgCode());
470  }
471  }
472 
473  break; // now that we have a proton we want to skip to next prong and associated BP tracks.
474  }
475  }
476  } // end if best fit track index matches current track index statement
477  } //end IsMC
478  else{ // Data
479  if(itrk==index){
480  if(ifit->PDG()==13){
481  fMuonTrkLngth->Fill(tracki->TotalLength());
482  if(fabs(tracki->Stop().X())<150 && fabs(tracki->Stop().Y())<150
483  && tracki->Stop().Z()<1200
484  && ((ifit->Chi2(2)+ifit->Chi2(3))/(ifit->Ndof(2)+ifit->Ndof(3)-4)) < 0.1
485  && itrk==lngst_trk_index){
486  Muon.SetPxPyPzE(ifit->FourMom().Px(), ifit->FourMom().Py(), ifit->FourMom().Pz(), ifit->FourMom().E());
487  haveMuon = true;
488  muon_track_length = tracki->TotalLength();
489  break; // now that we have a muon we want to skip to the next prong and associated BP tracks.
490  }
491  }
492  if(ifit->PDG()==2212){
493  fProtonTrkLngth->Fill(tracki->TotalLength());
494  if(tracki->TotalLength() < 125.0
495  && ((ifit->Chi2(2)+ifit->Chi2(3))/(ifit->Ndof(2)+ifit->Ndof(3)-4)) < 0.02){
496  haveProton = true;
497  Proton.SetPxPyPzE(ifit->FourMom().Px(), ifit->FourMom().Py(), ifit->FourMom().Pz(), ifit->FourMom().E());
498  proton_track_length = tracki->TotalLength();
499 
500  break;
501  }
502  }
503  }
504  }
505 
506  } // end loop over itrk
507  } // end loop over 3d prongs
508 
509  if(prongs.size()==2 && haveMuon && haveProton && slices[islice]->NCell()>=20){
510  if(fIsMC){
511  // Sort neutrinos associated with each slice and sort them according
512  // to purity.
513  std::vector< cheat::NeutrinoEffPur > iNuEffPur = NuEffPur[islice];
514  //int purest_index = 0;
515  for(int ipur=0; ipur<int(iNuEffPur.size())-1;ipur++){
516  if(iNuEffPur[ipur].purity<iNuEffPur[ipur+1].purity){
517  cheat::NeutrinoEffPur temp = iNuEffPur[ipur];
518  iNuEffPur[ipur] = iNuEffPur[ipur+1];
519  iNuEffPur[ipur+1] = temp;
520  }
521  }
522  // Get most pure neutrino
523  simb::MCNeutrino nu = iNuEffPur[0].neutrinoInt->GetNeutrino();
524  simb::MCParticle nupart = nu.Nu();
525 
526  TVector3 neu, mu, expected_proton, proton, reconu;
527  // Instead of using truth neutrino energy use QE reconstructed neutrino energy.
528  // This allows the same values to be used in data and MC.
529  // Also, initially we set the neutrino to have unit magnitude as we need the angle
530  // between the neutrino and the reco muon to calculate the neutrino energy.
531  neu.SetXYZ(0,0,1);
532  mu.SetXYZ(Muon.Px(), Muon.Py(), Muon.Pz());
533  double mp = 0.938;
534  double mn = 0.9396;
535  double cosmunu = mu.Dot(neu)*1./(mu.Mag()*neu.Mag());
536  double reconuE = 0.5*(mp*mp - mn*mn + 2.*mn*Muon.E() - 0.1057*0.1057)*1./(mn - Muon.E() + mu.Mag()*cosmunu);
537  // Now set the neutrino 3 momentum magnitude equal to the reconstructed neutrino energy.
538  neu.SetXYZ(0,0,reconuE);
539  // Determine expected proton 3 momentum.
540  expected_proton = neu - mu;
541  //slcE_proton = slcE_neu - mu/;
542  proton.SetXYZ(Proton.Px(), Proton.Py(), Proton.Pz());
543  reconu = proton + mu;
544 
545  float projection = proton.Dot(expected_proton)*1./proton.Mag()*1./expected_proton.Mag();
546  // float slcEprojection = proton.Dot(slcE_proton)*1./proton.Mag()*1./slcE_proton.Mag();
547 
548  // fPrjxnDiff->Fill(projection - slcEprojection);
549  fPrDotXPr->Fill(projection);
550  // fPrDotXPrSlcE->Fill(slcEprojection);
551  fPrMomFracBias->Fill((expected_proton.Mag()-proton.Mag())/expected_proton.Mag());
552  if(projection > 0.99){
553  fProtonPurityGoodSlcE->Fill(proton_pdg);
554  fMuonPurityGoodSlcE->Fill(muon_pdg);
555  fRecoNuE->Fill(reconuE);
556  fProtonTrkLngthGood->Fill(proton_track_length);
557  fMuonTrkLngthGood->Fill(muon_track_length);
558  fMuonMomentum->Fill(mu.Mag());
559  fProtonMomentum->Fill(proton.Mag());
560  fPrDelMomFracDiff->Fill((true_proton_mom - proton.Mag())/true_proton_mom);
561  fPrCosTheta->Fill(proton.Z()/proton.Mag());
562  fMuVPrCosTheta->Fill(mu.Z()/mu.Mag(), proton.Z()/proton.Mag());
563  // Print to screen (out file) to create a list of events with protons
564  // std::cout<<"bubbles : "<< e.run()
565  //<<" "<< e.subRun()
566  // <<" "<< e.id().event()
567  // <<" "<< islice
568  // <<std::endl;
569  }
570  fNuMomentum->Fill(reconu.Mag());
571 
572  //===================================================================
573  // The projection of the two proton vectors could be used to help ID
574  // CCQE events.
575  //===================================================================
576  if(projection>0.90)fEvtQEPurity1->Fill(nu.InteractionType());
577  if(projection>0.91)fEvtQEPurity2->Fill(nu.InteractionType());
578  if(projection>0.92)fEvtQEPurity3->Fill(nu.InteractionType());
579  if(projection>0.93)fEvtQEPurity4->Fill(nu.InteractionType());
580  if(projection>0.94)fEvtQEPurity5->Fill(nu.InteractionType());
581  if(projection>0.95)fEvtQEPurity6->Fill(nu.InteractionType());
582  if(projection>0.96)fEvtQEPurity7->Fill(nu.InteractionType());
583  if(projection>0.97)fEvtQEPurity8->Fill(nu.InteractionType());
584  if(projection>0.98)fEvtQEPurity9->Fill(nu.InteractionType());
585  if(projection>0.99)fEvtQEPurity10->Fill(nu.InteractionType());
586  /*if(slcEprojection>0.90)fEvtQEPurity1slcE->Fill(nu.InteractionType());
587  if(slcEprojection>0.91)fEvtQEPurity2slcE->Fill(nu.InteractionType());
588  if(slcEprojection>0.92)fEvtQEPurity3slcE->Fill(nu.InteractionType());
589  if(slcEprojection>0.93)fEvtQEPurity4slcE->Fill(nu.InteractionType());
590  if(slcEprojection>0.94)fEvtQEPurity5slcE->Fill(nu.InteractionType());
591  if(slcEprojection>0.95)fEvtQEPurity6slcE->Fill(nu.InteractionType());
592  if(slcEprojection>0.96)fEvtQEPurity7slcE->Fill(nu.InteractionType());
593  if(slcEprojection>0.97)fEvtQEPurity8slcE->Fill(nu.InteractionType());
594  if(slcEprojection>0.98)fEvtQEPurity9slcE->Fill(nu.InteractionType());
595  if(slcEprojection>0.99)fEvtQEPurity10slcE->Fill(nu.InteractionType());*/
596  //===================================================================
597 
598  if(nu.InteractionType()==1001){
599  fPrjxnVsRecoNuE->Fill(reconuE,projection);
600  fTrueNuE->Fill(nupart.E());
601  fRecoNuETrueQE->Fill(reconuE);
602  fPrDotXPrTrueQE->Fill(projection);
603  fPrMomFracBiasTrueQE->Fill((expected_proton.Mag()-proton.Mag())/expected_proton.Mag());
604  fNuMomentumTrueQE->Fill(reconu.Mag());
605  fMuonMomentumTrueQE->Fill(mu.Mag());
606  fProtonMomentumTrueQE->Fill(proton.Mag());
607 
608  // Fill particle purity plots allowing for purity measurements when
609  // selecting only true CCQE events. (fPurityGood histograms)
610  fMuonPurityGood->Fill(muon_pdg);
611  fProtonPurityGood->Fill(proton_pdg);
612  if(isTruMuon && isTruProton){
613  fPrDotXPrTrue->Fill(projection);
614  fPrMomFracBiasTrue->Fill((expected_proton.Mag()-proton.Mag())/expected_proton.Mag());
615  fRecoNuETrue->Fill(reconuE);
616  fMuonMomentumTrue->Fill(mu.Mag());
617  fProtonMomentumTrue->Fill(proton.Mag());
618  fNuMomentumTrue->Fill(reconu.Mag());
619  }
620  }
621 
622  } // end IsMC
623  else{ // Data
624  // Instead of using truth neutrino energy use QE reconstructed neutrino energy.
625  // This allows the same values to be used in data and MC.
626  // Also, initially we set the neutrino to have unit magnitude as we need the angle
627  // between the neutrino and the reco muon to calculate the neutrino energy.
628  TVector3 neu, mu, expected_proton, proton, reconu;
629  neu.SetXYZ(0,0,1);
630  mu.SetXYZ(Muon.Px(), Muon.Py(), Muon.Pz());
631  double mp = 0.938;
632  double mn = 0.9396;
633  double cosmunu = mu.Dot(neu)*1./(mu.Mag()*neu.Mag());
634  double reconuE = 0.5*(mp*mp - mn*mn + 2.*mn*Muon.E() - 0.1057*0.1057)*1./(mn - Muon.E() + mu.Mag()*cosmunu);
635  // Now set the neutrino 3 momentum magnitude equal to the reconstructed neutrino energy.
636  neu.SetXYZ(0,0,reconuE);
637  // Determine expected proton 3 momentum.
638  expected_proton = neu - mu;
639  //slcE_proton = slcE_neu - mu/;
640  proton.SetXYZ(Proton.Px(), Proton.Py(), Proton.Pz());
641  reconu = proton + mu;
642  float projection = proton.Dot(expected_proton)*1./proton.Mag()*1./expected_proton.Mag();
643  fPrDotXPr->Fill(projection);
644 
645  if(projection > 0.99){
646  fRecoNuE->Fill(reconuE);
647  fProtonTrkLngthGood->Fill(proton_track_length);
648  fMuonTrkLngthGood->Fill(muon_track_length);
649  fMuonMomentum->Fill(mu.Mag());
650  fProtonMomentum->Fill(proton.Mag());
651  fPrCosTheta->Fill(proton.Z()/proton.Mag());
652  fMuVPrCosTheta->Fill(mu.Z()/mu.Mag(), proton.Z()/proton.Mag());
653  }
654  }
655 
656  } // end event selection
657 
658 
659  } // end loop over slices
660 
661 } // end analyze module
662 
664 
665 
666  /*
667  art::PtrVector<rb::CellHit> hits = track->AllCells();
668  double proton_GeV = 1.0E-3;
669  int n_uncalib_hits = 0;
670  // Make reco hit for every cell hit in order to get calibrated cell energy.
671  for(unsigned int j=0;j<hits.size();j++){
672  const rb::CellHit *cellhit = track->Cell(j).get();
673  double Wtemp = track->W(cellhit);
674  rb::RecoHit rhit(hCal->MakeRecoHit(*cellhit,Wtemp));
675  if(rhit.IsCalibrated()) proton_GeV += rhit.GeV();
676  else n_uncalib_hits++;
677  }
678  */
double E(const int i=0) const
Definition: MCParticle.h:232
back track the reconstruction to the simulation
art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
int PdgCode() const
Definition: MCParticle.h:211
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
art::ServiceHandle< cheat::BackTracker > bt
double GetX() const
Definition: Vertex.h:23
double fTotalPOT
Total POT for one subrun.
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::string fGeneratorLabel
Where to find spill information in MC.
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
const int Ndof(int index) const
Definition: FitSum.h:68
double GetY() const
Definition: Vertex.h:24
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
std::string fVertexLabel
Where to find vertices.
art::ServiceHandle< calib::Calibrator > hCal
std::string fSpillLabel
Where to find spill information in data.
double GetZ() const
Definition: Vertex.h:25
const int PDG() const
Definition: FitSum.h:73
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void analyze(art::Event const &e) override
std::string fProngLabel
Where to find prongs.
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
int InteractionType() const
Definition: MCNeutrino.h:150
std::string fProngInstance
Which type of prong.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
const TLorentzVector FourMom() const
Definition: FitSum.h:72
BreakPointProtonAna(fhicl::ParameterSet const &p)
#define pot
std::string fTrackLabel
Where to find tracks.
TString mp
Definition: loadincs.C:4
const double j
Definition: BetheBloch.cxx:29
reference at(size_type n)
Definition: PtrVector.h:365
const double Chi2(int index) const
Definition: FitSum.h:61
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
size_type size() const
Definition: PtrVector.h:308
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
float GeV() const
Definition: RecoHit.cxx:69
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fFitSumLabel
Where to find FitSum objects.
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
BreakPointProtonAna & operator=(BreakPointProtonAna const &)=delete
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
std::string fMCTruthLabel
Where to find MCTruth information.
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
ProductToken< T > consumes(InputTag const &)
bool failedToGet() const
Definition: Handle.h:196
static constexpr Double_t sr
Definition: Munits.h:164
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
virtual void endSubRun(const art::SubRun &sr) override