ReMIdDedx_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////
2 // \file ReMIdDedx_module.cc
3 // \brief A module for creating trees which can be used for kNN training and generating sample histograms for ReMId
4 // \version $Id: ReMIDTrain_module.cc,v 2013/02/21 16:30:29 raddatz Exp $
5 // \author Nicholas Raddatz - raddatz@physics.umn.edu
6 ////////////////////////////////////////////////////////////////////
7 #include <string>
8 
9 // ROOT includes
10 #include "TVector3.h"
11 #include "TTree.h"
12 #include "TH1D.h"
13 #include "TH2D.h"
14 #include "TMath.h"
15 
16 // Framework includes
17 #include "fhiclcpp/ParameterSet.h"
19 #include "Utilities/AssociationUtil.h"
24 
25 // NOvA includes
26 #include "Geometry/Geometry.h"
27 #include "Geometry/LiveGeometry.h"
29 #include "ReMId/ReMId.h"
30 #include "NumuEnergy/NumuE.h"
31 #include "QEEventFinder/QePId.h"
32 #include "GeometryObjects/Geo.h"
35 #include "MCCheater/BackTracker.h"
36 #include "RecoBase/Track.h"
37 #include "Simulation/Simulation.h"
39 #include "SummaryData/POTSum.h"
40 #include "SummaryData/SpillData.h"
41 #include "MCCheater/BackTracker.h"
42 
43 
44 
45 namespace remid {
46 
47  // A module to analyze ReMId objects for pid training
48  class ReMIdDedx : public art::EDAnalyzer {
49  public:
50  explicit ReMIdDedx(fhicl::ParameterSet const& pset);
51  ~ReMIdDedx();
52 
53  void analyze(art::Event const& evt);
54 
55  void beginJob();
56  void beginRun(const art::Run& run);
57  void endJob();
58  //virtual void beginSubRun(art::SubRun& sr);
59 
60  private:
62  std::string fSliceLabel; ///<Where to find reconstructed slices
63  std::string fTrackLabel; ///<Where to find recondtructed tracks
64  std::string fPidLabel; ///<Where to find the pids
70 
71  TH1D* fPOT; ///<Histogram of POT passing spill cuts
78 
79  double fTotalPOT;
80  bool fND; ///< Is detector ND?
81 
82  // variables used in filling the kNN training trees
83  float fTrackLength;
85  //float fDedxSep;
86  //float fScatSep;
87  //float fnuE;
88  //float fvisE;
89  //float frecoE;
90  //float fMeasFrac;
91 
92  TTree* DedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
93  TTree* ProtonDedxSample; ///< Tree for holding dE/dx variables needed for making signal sample histograms
94 
95  // variables used in filling the sample histogram trees
96  float fDedx;
97  int fnDedx;
99  float fDistDedx;
100  float fProtonDedx;
104  //float fCosScat;
105  //float fDistLastScat;
106  //float fScatVar;
107  //int fnScat;
108  //float fDistScat;
109  //int fnPlane;
110  //int fpdg;
111  //int fintType;
112  bool fVert;
113  bool fUsed;
114  float fzStart;
115  float fzEnd;
119  float fProtonzEnd;
132 
133 
138 
139  //novadaq::cnv::DetId fDetID; ///< Detector ID in nova daq convention namespace typedef
140 
141 
142  //unsigned int fdibfirst = 0;
143  //unsigned int fdiblast = 0;
144 
146 
147 
148  };
149 
150 }
151 
152 
153 namespace remid{
154 
155 
156  //.......................................................................
157 
159  EDAnalyzer(pset)
160  {
161  fHitModuleLabel = (pset.get< std::string >("HitModuleLabel"));
162  fSliceLabel = (pset.get< std::string >("SliceLabel"));
163  fTrackLabel = (pset.get< std::string >("TrackLabel"));
164  fPidLabel = (pset.get< std::string >("PidLabel" ));
165  fGeneratorLabel = pset.get< std::string >("GeneratorInput");
166  fNumiBeamLabel = pset.get< std::string >("NuMIBeamInput");
167  fSpillQualLabel = pset.get< std::string >("SpillQualInput");
168  fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
169  fQepidLabel = pset.get< std::string >("QepidInput");
170  }
171 
172  //.......................................................................
173 
175  {
176  }
177 
178  //.......................................................................
179 
181  {
183 
184  DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
185  DedxSample->Branch("dedx", &fDedx, "dedx/F");
186  DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
187  DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
188  DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
189  DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
190  DedxSample->Branch("vert", &fVert, "vert/O");
191  DedxSample->Branch("used", &fUsed, "used/O");
192  DedxSample->Branch("zStart",&fzStart,"zStart/F");
193  DedxSample->Branch("zEnd",&fzEnd,"zEnd/F");
194  DedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
195  DedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
196  DedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
197  DedxSample->Branch("passVertX",&fPassVertX,"passVertX/O");
198  DedxSample->Branch("passVertY",&fPassVertY,"passVertY/O");
199  DedxSample->Branch("passVertZ",&fPassVertZ,"passVertZ/O");
200  DedxSample->Branch("passProtonRemid",&fPassProtonRemid,"passProtonRemid/O");
201  DedxSample->Branch("passDedxRatio",&fPassDedxRatio,"passDedxRatio/O");
202  DedxSample->Branch("passEDiffZTest",&fPassEDiffZTest,"passEDiffZTestX/O");
203  DedxSample->Branch("trueNumuCCQE",&fTrueNumuCCQE,"trueNumuCCQE/O");
204 
205  ProtonDedxSample = tfs->make<TTree>("fProtonDedxSample","Variables for Creating Proton Dedx Sample Histograms");
206  ProtonDedxSample->Branch("Protondedx", &fProtonDedx, "Protondedx/F");
207  ProtonDedxSample->Branch("ProtonnDedx", &fProtonnDedx, "ProtonnDedx/I");
208  ProtonDedxSample->Branch("ProtonnDedxUsed", &fProtonnDedxUsed, "ProtonnDedxUsed/I");
209  ProtonDedxSample->Branch("ProtondistFromEnd", &fProtonDistDedx, "ProtondistFromEnd/F");
210  ProtonDedxSample->Branch("ProtontrackLength", &fProtonTrackLength, "ProtontrackLength/F");
211  ProtonDedxSample->Branch("Protonvert", &fProtonVert, "Protonvert/O");
212  ProtonDedxSample->Branch("Protonused", &fProtonUsed, "Protonused/O");
213  ProtonDedxSample->Branch("ProtonzStart",&fProtonzStart,"ProtonzStart/F");
214  ProtonDedxSample->Branch("ProtonzEnd",&fProtonzEnd,"ProtonzEnd/F");
215  ProtonDedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
216  ProtonDedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
217  ProtonDedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");
218  ProtonDedxSample->Branch("passVertX",&fPassVertX,"passVertX/O");
219  ProtonDedxSample->Branch("passVertY",&fPassVertY,"passVertY/O");
220  ProtonDedxSample->Branch("passVertZ",&fPassVertZ,"passVertZ/O");
221  ProtonDedxSample->Branch("passProtonRemid",&fPassProtonRemid,"passProtonRemid/O");
222  ProtonDedxSample->Branch("passDedxRatio",&fPassDedxRatio,"passDedxRatio/O");
223  ProtonDedxSample->Branch("passEDiffZTest",&fPassEDiffZTest,"passEDiffZTestX/O");
224  ProtonDedxSample->Branch("trueNumuCCQE",&fTrueNumuCCQE,"trueNumuCCQE/O");
225 
226  fPOT = tfs->make<TH1D>("POT", ";;Total POT",1, 0.0, 1.0);
227  fEventsPassingSpillCuts = tfs->make<TH1D>("EventsPassingSpillCuts", ";;Events Passing Spill Cuts",1, 0.0, 1.0);
228  fSlicesPassingSelectionCuts = tfs->make<TH1D>("SlicesPassingSelectionCuts", ";;Slices Passing Selection Cuts",1, 0.0, 1.0);
229  fTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("TrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
230  fSuperTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("SuperTrueSlicesPassingSelectionCuts", ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
231  fTrueProtonEFromTrackLength = tfs->make<TH2D>("TrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
232  fSuperTrueProtonEFromTrackLength = tfs->make<TH2D>("SuperTrueProtonEFromTrackLength", ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
233 
234  fTotalPOT = 0.0;
235 
236 
237 
238  }
239 
240  //......................................................................
241 
243 
244  {
245  fND = false;
246 
247  novadaq::cnv::DetId detID = fgeom->DetId();
248  if (detID == novadaq::cnv::kNEARDET) fND = true;
249 
250  return;
251  }
252 
253 
254  //.......................................................................
255 
257  {
258  //Only written for ND
259 
260  if (!fND) return;
261 
262  //if (evt.event() == 443183){
263  // std::cout<<"Right event!"<<std::endl;
264  //}
265 
266  //Does this event pass spill cuts?
268  if(bt->HaveTruthInfo()){
269  evt.getByLabel(fGeneratorLabel, spillPot);
270  }
271  else{
272  evt.getByLabel(fNumiBeamLabel, spillPot);
273  }
274 
275  if (spillPot.failedToGet()){
276  mf::LogError("Susie") << "Spill Data not found for real data event";
277  return;
278  }
279 
280  //SL -- only cut on beam quality for data
281  if (!bt->HaveTruthInfo()){
282  if (std::abs(spillPot->deltaspilltimensec) > 0.5e9) return;
283  if ((spillPot->spillpot)*1e12 < 2e12) return;
284  if ((spillPot->hornI < -202)||(spillPot->hornI > -198)) return;
285  if ((spillPot->posx < -2.00)||(spillPot->posx > 2.00)) return;
286  if ((spillPot->posy < -2.00)||(spillPot->posy > 2.00)) return;
287  if ((spillPot->widthx < 0.57)||(spillPot->widthx > 1.58)) return;
288  if ((spillPot->widthy < 0.57)||(spillPot->widthy > 1.58)) return;
289  }
290 
291  //Data quality spill cuts
292 
293  // Get the EventQuality information
295  evt.getByLabel(fSpillQualLabel, spillQual);
296 
297  if(spillQual.failedToGet()){
298  mf::LogError("CAFMaker") << "Spill EventQuality not found for real data event";
299  return;
300  }
301 
302  if (spillQual->fracdcm3hits > 0.45) return;
303  if (spillQual->nmissingdcms > 0) return;
304 
305 
306  //Filling spill pot after it passes the cuts
307  float spillpot = spillPot->spillpot;
308  if (!bt->HaveTruthInfo()) spillpot *= 1e12;
309  fTotalPOT += spillpot;
310 
311  fEventsPassingSpillCuts->Fill(0.5);
312 
313  //Getting event hits
315  evt.getByLabel(fHitModuleLabel, hithdl);
316 
318  for(size_t h = 0; h < hithdl->size(); ++h){
319  art::Ptr<rb::CellHit> hit(hithdl, h);
320  allhits.push_back(hit);
321  }
322 
323 
324 
325  // get the slices out of the event
327  evt.getByLabel(fSliceLabel,slicecol);
328 
329 
330  art::PtrVector<rb::Cluster> slicelist;
331 
332  for(unsigned int i = 0; i<slicecol->size();++i){
333  art::Ptr<rb::Cluster>slice(slicecol,i);
334  slicelist.push_back(slice);
335  }
336 
337  // get assosciations between slices and tracks
338  art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);
339 
340  // get assosciations between slices and numu energy
341  art::FindManyP<numue::NumuE> fndmnyNumuEnergy(slicecol, evt, fNumuEnergyLabel);
342  // get assosciations between slices and qepid
343  art::FindManyP<qeef::QePId> fndmnyQepid(slicecol, evt, fQepidLabel);
344 
345  //if (evt.event() == 443183){
346  // std::cout<<"Right event 2!"<<std::endl;
347  //}
348 
349  // loop over the slices
350  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
351 
352  // If slice is noise, skip it
353  if(slicelist[iSlice]->IsNoise()){ continue; }
354 
355  // get all of the tracks associated to this slice
356  std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
357  // if no tracks move on
358  if(tracks.size() == 0){ continue; }
359  //Require precisely two Kalman tracks
360  if(tracks.size() != 2){ continue; }
361 
362  //if ((evt.event() == 443183)&&(iSlice == 2)){
363  // std::cout<<"Right event and slice!"<<std::endl;
364  //}
365 
366  //Require BOTH are 3D
367  bool all3D = true;
368 
369  // get the pid that belongs to this recotrack
370  art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);
371 
372  int bestidx = -1;
373  double bestpid = -1;
374  // loop over all of the tracks and find the most muon like one;
375  for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
376  if(!tracks[itrk]->Is3D()){
377  all3D = false;
378  continue;
379  }
380  // get the pid that belongs to this track
381  art::Ptr<remid::ReMId> pid = foids.at(itrk);
382  if(pid->Value() > bestpid){
383  bestpid = pid->Value();
384  bestidx = itrk;
385  }
386  }//End of loop over tracks
387 
388  if (!all3D){ continue; }
389 
390  //Default all pass cuts to false
391  fPassOffTrackHits = false;
392  fPassOffTrackEnergy = false;
393  fPassHadEOnMuonTrack = false;
394  fPassVertX = false;
395  fPassVertY = false;
396  fPassVertZ = false;
397  fPassProtonRemid = false;
398  fPassDedxRatio = false;
399  fPassEDiffZTest = false;
400 
401  //Slightly harsh ReMId cut on muon track
402  if(bestidx == -1 || bestpid < 0.8){ continue; }
403 
404  // if ((evt.event() == 443183)&&(iSlice == 2)){
405  //std::cout<<"Right event and slice2!"<<std::endl;
406  //}
407 
408  if((bestidx != 0) && (bestidx !=1)){
409  std::cout<<"BUMMER! A weird track id for muon that shouldn't happen!!"<<bestidx<<std::endl;
410  continue;
411  }
412 
413  int protonidx = -1;
414  if (bestidx == 0){
415  protonidx = 1;
416  }
417  if (bestidx == 1){
418  protonidx = 0;
419  }
420 
421  if(protonidx == -1){
422  std::cout<<"BUMMER! A weird track id for proton that shouldn't happen!!"<<std::endl;
423  continue;
424  }
425 
426  //Filling muon track length properties
427  fzStart = tracks[bestidx]->Start().Z();
428  fzEnd = tracks[bestidx]->Stop().Z();
429  fTrackLength = tracks[bestidx]->TotalLength();
430 
431  //Filling proton track length properties
432  fProtonzStart = tracks[protonidx]->Start().Z();
433  fProtonzEnd = tracks[protonidx]->Stop().Z();
434  fProtonTrackLength = tracks[protonidx]->TotalLength();
435 
436  art::Ptr<remid::ReMId> pid = foids.at(bestidx);
437  art::Ptr<remid::ReMId> protonpid = foids.at(protonidx);
438 
439  //Looking at numu energy things
440  numue::NumuE numuE;
441  if (!fndmnyNumuEnergy.isValid()){ continue; }
442  std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
443  if (energyVec.empty()){ continue; }
444  numuE = *energyVec[0];
445 
446  //Require decent energy object; not too much hadronic energy in muon catcher
447  if (numuE.E() == -5){ continue; }
448  if (!(numuE.NDHadCalCat() + numuE.NDHadCalTran() < 0.03)){ continue; }
449  //Not much hadronic energy contamination on muon track
450  //if (!(numuE.HadTrkE() < 0.025)){ continue; }
451  if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }
452 
453  //if ((evt.event() == 443183)&&(iSlice == 2)){
454  // std::cout<<"Right event and slice3!"<<std::endl;
455  //}
456 
457  //Looking at qepid things
458  qeef::QePId qepid;
459  if (!fndmnyQepid.isValid()){ continue; }
460  std::vector<art::Ptr<qeef::QePId> > qVec = fndmnyQepid.at(iSlice);
461  if (qVec.empty()){ continue; }
462  qepid = *qVec[0];
463 
464  //Require ratio of average dE/dx of two tracks is good
465  //if (!(qepid.Dedx() > 1.75)){ continue; }
466  if (qepid.Dedx() > 1.75){ fPassDedxRatio = true; }
467  //Require Ztest of two methods of calculating QE neutrino energy is good
468  //if (!(qepid.EdiffZ() > -1.3)){ continue; }
469  if (qepid.EdiffZ() > -1.3){ fPassEDiffZTest = true; }
470 
471 
472  //Slice energy mostly on two tracks
473  double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy()+tracks[protonidx]->CalorimetricEnergy());
474  //if (!(calEDiff < 0.1)){ continue; }
475  if (calEDiff < 0.1){ fPassOffTrackEnergy = true; }
476 
477  //Slice hits mostly on two tracks
478  int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell()+tracks[protonidx]->NCell());
479  //if (!(hitDiff < 4)){ continue; }
480  if (hitDiff < 4){ fPassOffTrackHits = true; }
481 
482  //Two tracks start near each other
483  double diffX = fabs(tracks[bestidx]->Start().X()-tracks[protonidx]->Start().X());
484  double diffY = fabs(tracks[bestidx]->Start().Y()-tracks[protonidx]->Start().Y());
485  double diffZ = fabs(tracks[bestidx]->Start().Z()-tracks[protonidx]->Start().Z());
486  //if (!(diffX < 6.0)){ continue; }
487  //if (!(diffY < 6.0)){ continue; }
488  //if (!(diffZ < 10.0)){ continue; }
489  if (diffX < 6.0){ fPassVertX = true; }
490  if (diffY < 6.0){ fPassVertY = true; }
491  if (diffZ < 10.0){ fPassVertZ = true; }
492 
493  //Require proton to have a low ReMId value
494  //if (!(protonpid->Value() < 0.15)){ continue; }
495  if (protonpid->Value() < 0.15){ fPassProtonRemid = true; }
496 
497  //if ((evt.event() == 443183)&&(iSlice == 2)){
498  // std::cout<<"Right event and slice4!"<<std::endl;
499  //}
500 
501  double trueProtonE = -1.0;
502 
503  //Check if true signal for sim events
504  fTrueNumuCCQE = false;
505  fTrueNumuCCWithProton = false;
507  if (bt->HaveTruthInfo()){
508  std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),allhits);
509  if (!funcReturn.empty()){
510  const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();
511  if((test_neutrino.CCNC()==simb::kCC) && (std::abs(test_neutrino.Nu().PdgCode()) == 14) && (test_neutrino.InteractionType()==simb::kCCQE)){
512  fTrueNumuCCQE = true;
513  //Find if matching proton, true E of proton
514  const std::vector<const sim::Particle*> test_particle = bt->HitsToParticle(tracks[protonidx]->AllCells());
515  if (!test_particle.empty()){
516  int pdg = test_particle[0]->PdgCode();
517  //std::cout<<" with pdg: "<<pdg<<std::endl;
518  if ((pdg == 2212)||(pdg == -2212)){
519  fTrueNumuCCWithProton = true;
520  trueProtonE = test_particle[0]->E();
521  //Test to see if matches to 80% efficiency and purity
522  std::set<int> protonTrackID;
523  protonTrackID.insert(test_particle[0]->TrackId());
524  double trackProtonEff = bt->HitCollectionEfficiency(protonTrackID, tracks[protonidx]->AllCells(), allhits, geo::kXorY, 0, false);
525  double trackProtonPur = bt->HitCollectionPurity(protonTrackID, tracks[protonidx]->AllCells(), 0, 0, false);
526  if ((trackProtonEff > 0.8)&&(trackProtonPur > 0.8)){
528  }
529  }
530 
531  }
532  }
533  }
534  }
535 
536  if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){
537  //if ((evt.event() == 443183)&&(iSlice == 2)){
538  // std::cout<<"Right event and slice5!"<<std::endl;
539  //}
541  fSlicesPassingSelectionCuts->Fill(0.5);
545  }
549  }
550  }
551  // get all the dedx variables filled
552  fnDedx = pid->NDedx()-1;
554  // loop over all of the dE/dx measurements
555  for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
556  fDedx = pid->Dedx(idedx);
557  fVert = pid->DedxVertex(idedx);
558  fUsed = pid->DedxUsed(idedx);
559  // get the distance from the end of the track
560  fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[bestidx]->TotalLength();
561  // fill the dedx sample tree
562  DedxSample->Fill();
563  }//Loop over DedxSample
564 
565  // get all the proton dedx variables filled
566  fProtonnDedx = protonpid->NDedx()-1;
567  fProtonnDedxUsed = protonpid->NMeasurementPlanes();
568  // loop over all of the dE/dx measurements
569  for(unsigned int sidedx = 0; sidedx<protonpid->NDedx(); ++sidedx){
570  fProtonDedx = protonpid->Dedx(sidedx);
571  fProtonVert = protonpid->DedxVertex(sidedx);
572  fProtonUsed = protonpid->DedxUsed(sidedx);
573  // get the distance from the end of the track
574  fProtonDistDedx = (1.0 - protonpid->DedxLocation(sidedx))*tracks[protonidx]->TotalLength();
575  // fill the dedx sample tree
576  //std::cout<<"now we fill!"<<std::endl;
577  ProtonDedxSample->Fill();
578  }//Loop over DedxSample
579 
580  }//If passes contaiment
581 
582 
583  } // end loop over slices
584 
585 
586  }
587 
588 
589  //.......................................................................
590 
592  {
593  fPOT->Fill(.5, fTotalPOT);
594  }
595 
596 
597  //.......................................................................
599  {
600 
601  bool contained = false;
602 
603  // check that this is a contained event
604 
605  // get some useful information here
606  // 1) how close is the slice to the edge of the detector
607  unsigned int slcncellsfromedge = 999999999;
608  for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
609  const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
610  const int plane = chit->Plane();
611  unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
612  if(cellsfromedge < slcncellsfromedge){
613  slcncellsfromedge = cellsfromedge;
614  }
615  }
616  // 2) how close is the track to the edge of the detector
617  // in planes
618  //SL - NOPE, need to ask slice, not track
619  unsigned int firstpl = slice->MinPlane();
620  unsigned int lastpl = slice->MaxPlane();
621  // in projected cells, just define directions here actually get projected cell below
622  TVector3 dirbwk = -trk->Dir();
623  TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
624 
625  // 3) hits in slice
626  int sliceNHit = slice->NCell();
627 
628  // 4) slice mean time
629  double slcTime = slice->MeanTNS();
630 
631  // 5) hits on track
632  int trackNHit = trk->NCell();
633 
634  //std::cout<<"In contained event loop!"<<std::endl;
635  // cut the slice number of cells from edge
636  if((sliceNHit > 20)&&(slcncellsfromedge > 1)&&(slcTime>50000)&&(slcTime<450000)&&(trackNHit > 5)){
637  if(firstpl > 1 && lastpl < 212){
638  //std::cout<<"In contained event loop! 1"<<std::endl;
639  int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
640  int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
641  // also make sure the track doesn't start too close to the end of the active region
642  if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
643  // if the track ends in the muon catcher, check that it doesn't clip the corner
644  //std::cout<<"In contained event loop! 2"<<std::endl;
645  if(trk->Stop().Z() < 1275){
646  //std::cout<<"In contained event loop! 2a"<<std::endl;
647  contained = true;
648  }
649  else{
650  //std::cout<<"In contained event loop! 2b"<<std::endl;
651  /*
652  // find the trajectory point closest in z to the transition plane
653  double mindist = 99999999;
654  unsigned int mintrajid = 0;
655  for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
656  double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
657  if(dist > 0 && dist < mindist){
658  mindist = dist;
659  mintrajid = itraj;
660  }
661  }
662  // get the direction at this trajectory point
663  TVector3 trandir = trk->Dir();
664  if(mintrajid < trk->NTrajectoryPoints() -2){
665  TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
666  }
667  */
668  //std::cout<<"In contained event loop! 3"<<std::endl;
669  //double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
670  double ypostranKirk = flivegeom->YPositionAtTransition(trk->Start(),trk->Dir());
671  if( ypostranKirk < 55){
672  //std::cout<<"In contained event loop! 4"<<std::endl;
673  contained = true;
674  }
675 
676  } // end else in muon catcher
677  } // end if track ends projections contained
678  } // end if track end planes contained
679  // } // end if expected number of diblocks
680  } // end if slice contained
681 
682 
683  return contained;
684 
685  }
686 
687 
688 
689 
690 }
691 
692 
693 namespace remid{
694 
696 
697 }
std::string fQepidLabel
size_t NTrajectoryPoints() const
Definition: Track.h:83
double EdiffZ() const
Definition: QePId.cxx:93
back track the reconstruction to the simulation
double Dedx(unsigned int i) const
Definition: ReMId.cxx:126
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
art::ServiceHandle< geo::LiveGeometry > flivegeom
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
bool fND
Is detector ND?
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
float E() const
Definition: Energy.cxx:27
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double Dedx() const
Definition: QePId.cxx:99
std::string fNumiBeamLabel
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
std::string fSpillQualLabel
TH1D * fPOT
Histogram of POT passing spill cuts.
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
TH1D * fTrueSlicesPassingSelectionCuts
std::string fTrackLabel
Where to find recondtructed tracks.
art::ServiceHandle< cheat::BackTracker > bt
TH2D * fSuperTrueProtonEFromTrackLength
charged current quasi-elastic
Definition: MCNeutrino.h:96
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
int nmissingdcms
of missing DCMs
Definition: EventQuality.h:23
TH2D * fTrueProtonEFromTrackLength
float abs(float number)
Definition: d0nt_math.hpp:39
double Value() const
Definition: PID.h:22
std::string fPidLabel
Where to find the pids.
Definition: Run.h:31
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
Float_t Y
Definition: plot.C:38
bool DedxVertex(unsigned int i) const
Definition: ReMId.cxx:138
Float_t Z
Definition: plot.C:38
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
int InteractionType() const
Definition: MCNeutrino.h:150
TTree * ProtonDedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
std::string fNumuEnergyLabel
signed long long int deltaspilltimensec
Definition: SpillData.h:24
double fracdcm3hits
fraction of DCM3 hits in horizontal modules
Definition: EventQuality.h:24
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
int NMeasurementPlanes() const
Definition: ReMId.cxx:223
float HadTrkE() const
Definition: NumuE.cxx:284
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double widthx
mm
Definition: SpillData.h:36
int evt
std::string fSliceLabel
Where to find reconstructed slices.
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
bool DedxUsed(unsigned int i) const
Definition: ReMId.cxx:132
double hornI
kA
Definition: SpillData.h:27
art::ServiceHandle< nova::dbi::RunHistoryService > frh
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double posy
mm
Definition: SpillData.h:35
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
float NDHadCalTran() const
Definition: NumuE.cxx:334
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
OStream cout
Definition: OStream.cxx:6
unsigned int NDedx() const
Definition: ReMId.cxx:102
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
TH1D * fSuperTrueSlicesPassingSelectionCuts
std::string fHitModuleLabel
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
TH1D * fSlicesPassingSelectionCuts
std::string fGeneratorLabel
ReMIdDedx(fhicl::ParameterSet const &pset)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
T * make(ARGS...args) const
double DedxLocation(unsigned int i) const
Definition: ReMId.cxx:120
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
TRandom3 r(0)
TH1D * fEventsPassingSpillCuts
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
void analyze(art::Event const &evt)
A PID for muons.
Definition: FillPIDs.h:10
T min(const caf::Proxy< T > &a, T b)
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void beginRun(const art::Run &run)
Event generator information.
Definition: MCNeutrino.h:18
Float_t X
Definition: plot.C:38
art::ServiceHandle< geo::Geometry > fgeom
double posx
mm
Definition: SpillData.h:34
float NDHadCalCat() const
Definition: NumuE.cxx:339
double widthy
mm
Definition: SpillData.h:37
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
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.