MuonRemoveAna_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \file MuonRemoveAna_module.cc
3 /// \brief A module to assess the quality of muon removal
4 ///
5 /// \version $ Id: MuonRemoveAna_module.cc, v 1.11 2012-10-26 22:19:03 ksachdev Exp $
6 /// \author Kanika Sachdev (ksachdev@physics.umn.edu)
7 ////////////////////////////////////////////////////////////////////////
8 
9 
13 #include "MCCheater/BackTracker.h"
14 #include "RecoBase/CellHit.h"
15 #include "RecoBase/RecoHit.h"
16 #include "RecoBase/Track.h"
17 #include "RecoBase/Cluster.h"
18 
19 
28 #include <fhiclcpp/ParameterSet.h>
30 
31 
32 #include <map>
33 #include <vector>
34 #include <cmath>
35 #include "TNtuple.h"
36 #include "TH1F.h"
37 
38 namespace murem
39 {
40 
42  {
43 
44  public:
45  explicit MuonRemoveAna(fhicl::ParameterSet const& pset);
46  virtual ~MuonRemoveAna();
47 
48  void analyze(art::Event const& evt);
49  void beginJob();
50  void endJob();
51  void reconfigure(const fhicl::ParameterSet& p);
52 
53  private:
54 
57 
58  int fEvtId;
59  int fIsCC;
60  float fLength;
61  float fNTracks;
62  float fCosZ;
63  float fNHits;
64  float fPePerHit;
65  float fRmsX;
66  float fRmsY;
67  float fRmsZ;
74  float fTrackEff;
75  float fTrackPur;
76 
77  std::vector<float> fMuonHitMip;
78  std::vector<float> fMuonHitGeV;
79 
80  TTree* fSliceTree;
82  TH1F* fdedxmu;
83  TH1F* fdedxsh;
84  TH1F* fpurl;
85  TH1F* fpurh;
86  TH1F* feffmuh;
87  TH1F* feffmul;
88  TH1F* feffhadl;
89  TH1F* feffhadh;
90 
97 
99 
100  };
101 
102 }
103 namespace murem
104 {
105 
106  /*********************************************************************************************/
107  ////////////////////////////////////////////////////////////////////////
108 
109 
111  : EDAnalyzer(p),
112  fTrkCleanUpAlg(0),
113  fTrkCleanUpPSet(p.get< fhicl::ParameterSet >("TrkCleanUpPSet") )
114  {
115  this->reconfigure(p);
116  }
117 
118 
119  /*********************************************************************************************/
120  ////////////////////////////////////////////////////////////////////////
121 
123  {
124  if(fTrkCleanUpAlg) delete fTrkCleanUpAlg;
125  }
126 
127 
128  /*********************************************************************************************/
129  ////////////////////////////////////////////////////////////////////////
130 
131 
133  {
134  fTrackLabel = p.get< std::string>("TrackLabel");
135  fMRTrackLabel = p.get< std::string>("MRTrackLabel");
136  fMCTruthLabel = p.get< std::string>("fMCTruthLabel");
137  fSliceLabel = p.get< std::string>("SliceLabel");
138  fMRSliceLabel = p.get< std::string>("MRSliceLabel");
139  fMrccLabel = p.get< std::string>("MrccLabel");
140  }
141 
142 
143  /*********************************************************************************************/
144  ////////////////////////////////////////////////////////////////////////
145 
147  {
148 
150 
151  fSliceTree = tfs->make<TTree>("fSliceTree","fSliceTree");
152 
153  fSliceTree->Branch("event", &fEvtId);
154  fSliceTree->Branch("isCC", &fIsCC);
155  fSliceTree->Branch("trackLength", &fLength);
156  fSliceTree->Branch("cosZ", &fCosZ);
157  fSliceTree->Branch("nTracks", &fNTracks);
158  fSliceTree->Branch("trackEff", &fTrackEff);
159  fSliceTree->Branch("trackPur", &fTrackPur);
160  fSliceTree->Branch("nHits", &fNHits);
161  fSliceTree->Branch("pePerHit", &fPePerHit);
162  fSliceTree->Branch("rmsX", &fRmsX);
163  fSliceTree->Branch("rmsY", &fRmsY);
164  fSliceTree->Branch("rmsZ", &fRmsZ);
165  fSliceTree->Branch("hadEnergy", &fTotHadEnergy);
166  fSliceTree->Branch("hadFracRemoved", &fHadFracRemoved);
167  fSliceTree->Branch("muEnergy", &fTotMuonEnergy);
168  fSliceTree->Branch("muFracRemaining", &fMuFracRemaining);
169  fSliceTree->Branch("trueSharedEnergy",&fTrueSharedEnergy);
170  fSliceTree->Branch("tecoSharedEnergy",&fRecoSharedEnergy);
171  fSliceTree->Branch("muonHitMip", &fMuonHitMip);
172  fSliceTree->Branch("muonHitGev", &fMuonHitGeV);
173 
174 
175  fMrccSliceTree = tfs->make<TTree>("fMrccSliceTree","fMrccSliceTree");
176 
177  fMrccSliceTree->Branch("event", &fEvtId);
178  fMrccSliceTree->Branch("isCC", &fIsCC);
179  fMrccSliceTree->Branch("trackLength", &fLength);
180  fMrccSliceTree->Branch("cosZ", &fCosZ);
181  fMrccSliceTree->Branch("nTracks", &fNTracks);
182  fMrccSliceTree->Branch("nHits", &fNHits);
183  fMrccSliceTree->Branch("pePerHit", &fPePerHit);
184  fMrccSliceTree->Branch("rmsX", &fRmsX);
185  fMrccSliceTree->Branch("rmsY", &fRmsY);
186  fMrccSliceTree->Branch("rmsZ", &fRmsZ);
187 
188  fdedxmu = tfs->make<TH1F>("dedxmu", "dE/dx in a Plane; dE/dx per Plane [GeV/cm]; ", 200,0,0.005);
189  fdedxsh = tfs->make<TH1F>("dedxsh", "dE/dx in a Plane; dE/dx per Plane [GeV/cm]; ", 200,0,0.005);
190 
191  fpurl = tfs->make<TH1F>("purl", "Purity Vs Lower Bound Value; dE/dx per Plane [GeV/cm]; Purity", 200,0,0.005);
192  fpurh = tfs->make<TH1F>("purh", "Purity Vs Upper Bound Value; dE/dx per Plane [GeV/cm]; Purity", 200,0,0.005);
193 
194  feffmuh = tfs->make<TH1F>("effmuh", "Signal Efficiency Vs Upper Bound Value; dE/dx per Plane [GeV/cm]; Efficiency", 200,0,0.005);
195  feffmul = tfs->make<TH1F>("effmul", "Signal Efficiency Vs Lower Bound Value; dE/dx per Plane [GeV/cm]; Efficiency", 200,0,0.005);
196 
197  feffhadl = tfs->make<TH1F>("effhadl", "Background Efficiency Vs Lower Bound Value; dE/dx per Plane [GeV/cm]; Efficiency", 200,0,0.005);
198  feffhadh = tfs->make<TH1F>("effhadh", "Background Efficiency Vs Upper Bound Value; dE/dx per Plane [GeV/cm]; Efficiency", 200,0,0.005);
199  }
200 
201  /*********************************************************************************************/
202  ////////////////////////////////////////////////////////////////////////
203 
205  {
206 
208 
210  fEvtId = evt.id().event();
211 
212  fIsCC = -5;
213  fLength = -5;
214  fNTracks = -5;
215  fNHits = -5;
216  fPePerHit = -5;
217  fRmsX = -5;
218  fRmsY = -5;
219  fRmsZ = -5;
220  fMuFracRemaining = -5;
221  fHadFracRemoved = -5;
222  fTotMuonEnergy = -5;
223  fTotHadEnergy = -5;
224  fTrueSharedEnergy= -5;
225  fRecoSharedEnergy= -5;
226  fTrackEff = -5;
227  fTrackPur = -5;
228  // Grab all slices
229 
231  evt.getByLabel( fSliceLabel, slicecol);
233 
234  for(int i = 0; i < (int)slicecol->size(); i++){
235  art::Ptr<rb::Cluster> sli(slicecol, i);
236  slices.push_back(sli);
237  }
238  art::FindManyP<rb::Track> fmt(slicecol, evt, fTrackLabel);
239 
240 
241  // Grab all mrcc slices
242 
244  evt.getByLabel( fMRSliceLabel, mrccslicecol);
245  art::PtrVector<rb::Cluster > mrccslices;
246 
247  for(int i = 0; i < (int)mrccslicecol->size(); i++){
248  art::Ptr<rb::Cluster> sli(mrccslicecol, i);
249  mrccslices.push_back(sli);
250  }
251  // art::FindManyP<rb::Track> fmtmrcc(mrccslicecol, evt, fMRTrackLabel);
252 
253 
254  // Make a vector of all cellhits. It is needed for backtracker functions.
255 
257  for(int s = 0; s<(int)slices.size(); s++){
258  for(int c = 0; c<(int)slices[s]->NCell(); c++)
259  allHits.push_back( slices[s]->Cell(c) );
260  }
261 
262 
263  // Make a vector of all mrcc cellhits. It is needed for backtracker functions.
264 
265  art::PtrVector< rb::CellHit> mrccallHits;
266  for(int s = 0; s<(int)mrccslices.size(); s++){
267  for(int c = 0; c<(int)mrccslices[s]->NCell(); c++)
268  mrccallHits.push_back( mrccslices[s]->Cell(c) );
269  }
270 
271  // Get a list of removed muons.
272 
274  evt.getByLabel( fMrccLabel, removedMuonCol);
276 
277  for(int i = 0; i < (int)removedMuonCol->size(); i++){
278  art::Ptr<simb::MCParticle> mu( removedMuonCol, i );
279  removedMuons.push_back( mu );
280  }
281 
282  // Fill slice level variables for mrcc slices.
283  int nMrccSlices = mrccslices.size();
284  for( int iSli = 0; iSli < nMrccSlices; iSli++){
285 
286  if(!mrccslices[iSli]->IsNoise()) {
287 
288  std::vector< cheat::NeutrinoEffPur > nuEffPur =
289  bt->SliceToNeutrinoInteractions(mrccslices[iSli], mrccallHits);
290 
291  if (nuEffPur.size() > 0 && nuEffPur[0].purity>0.6){
292 
293  const simb::MCNeutrino& mcNu = nuEffPur[0].neutrinoInt->GetNeutrino();
294  if ( abs(mcNu.Nu().PdgCode()) == 14 ){
295 
296  fIsCC = mcNu.CCNC()==simb::kCC;
297  fNHits = mrccslices[iSli]->NCell();
298 
299  // std::vector< art::Ptr<rb::Track> > mrcctracks = fmtmrcc.at(iSli);
300 
301  // fNTracks = mrcctracks.size();
302  // rb::Track longestTrack;
303  // double maxLength = -1;
304  // double length;
305 
306  // int nMrccTracks = mrcctracks.size();
307  // for( int iTrk = 0; iTrk < nMrccTracks; iTrk++){
308  // length = mrcctracks[iTrk]->TotalLength();
309  // if( length > maxLength ){
310  // maxLength = length;
311  // longestTrack = *mrcctracks[iTrk];
312  // }
313  // }
314 
315  // fCosZ = longestTrack.Dir().Z();
316  // fLength = longestTrack.TotalLength();
317 
318  }// neutrino is muon type
319 
320  }// can match to a neutrino interaction
321  }// non-noise slice
322  }// end loop over mrc slices
323 
324 
325  fIsCC = -5;
326  fLength = -5;
327  fNTracks = -5;
328  fNHits = -5;
329  fPePerHit = -5;
330  fRmsX = -5;
331  fRmsY = -5;
332  fRmsZ = -5;
333 
334 
335  // Fill tree for original slices
336 
337  int nSlices = slices.size();
338  for( int iSli = 0; iSli < nSlices; iSli++){
339 
340  if(!slices[iSli]->IsNoise()) {
341  std::vector< cheat::NeutrinoEffPur > nuEffPur = bt->SliceToNeutrinoInteractions(slices[iSli], allHits);
342 
343  if (nuEffPur.size() > 0 && nuEffPur[0].purity > 0.6){
344 
345  const simb::MCNeutrino& mcNu = nuEffPur[0].neutrinoInt->GetNeutrino();
346  if ( mcNu.Nu().PdgCode() == 14 ){
347 
348  fIsCC = mcNu.CCNC()==simb::kCC;
349  fNHits = slices[iSli]->NCell();
350  std::vector< art::Ptr<rb::Track> > tracks = fmt.at(iSli);
351 
352  fNTracks = tracks.size();
353  rb::Track longestTrack;
354  double maxLength = -1;
355  double length;
356 
357  int nTracks = tracks.size();
358  for( int iTrk = 0; iTrk < nTracks; iTrk++){
359  length = tracks[iTrk]->TotalLength();
360  if( length > maxLength ){
361  maxLength = length;
362  longestTrack = *tracks[iTrk];
363  }
364  }
365 
366  fCosZ = longestTrack.Dir().Z();
367  fLength = longestTrack.TotalLength();
368 
369  // Code to calculate the muon fraction remaining and hadronic fraction removed
370  // Also, to compare the true shared energy with what the TrackCleanUpAlg calls
371  // shared energy, we need to find the track that was removed.
372 
373  if (fIsCC){
374 
375  std::vector< const sim::Particle*> particles = bt->MCTruthToParticles( nuEffPur[0].neutrinoInt );
376 
377  int muMotherId;
378  int michelElID = -1;
379  std::set<int> primMuTrackIDsSet;
380 
381  for( int p = 0; p < (int)particles.size(); p++)
382  if( abs(particles[p]->PdgCode()) == 13 &&
383  (particles[p]->Process() == "Primary" ||
384  particles[p]->Process() == "primary")){
385  muMotherId = particles[p]->TrackId();
386  primMuTrackIDsSet.insert(muMotherId);
387  }
388 
389  for( int p = 0; p < (int)particles.size(); p++){
390  if( particles[p]->Mother() == muMotherId ){
391 
392  if( abs(particles[p]->PdgCode()) == 11 &&
393  (particles[p]->Process() == "muMinusCaptureAtRest" ||
394  particles[p]->Process() == "Decay")){
395  michelElID = particles[p]->TrackId();
396 
397  }
398  else{
399  primMuTrackIDsSet.insert( particles[p]->TrackId() );
400  }
401  }
402  }// end loop over particles associated with this MCTruth
403 
404 
405  // Look for the track that was removed by MuonRemove
406 
407  rb::Track muonTrack;
408  for( int t = 0; t< (int)tracks.size(); t++){
409  if( !tracks[t]->Is3D() || tracks[t]->NCell()==0 )
410  continue;
411  TVector3 trkStart= tracks[t]->Start();
412  TVector3 trkStop = tracks[t]->Stop();
413 
414  for(int mu = 0; mu < (int)removedMuons.size(); mu++){
415 
416  if( removedMuons[mu]->Position(0).Vect() == trkStart &&
417  removedMuons[mu]->EndPosition().Vect() == trkStop ){
418  muonTrack = *(tracks[t]);
419 
420  }
421  }
422  }// end loop over tracks
423 
424 
425  float totEnergy = 0;
426  float muonEnergy = 0;
427 
428  int nMuCells = muonTrack.NCell();
429 
430  // vectors to store the indices for planes that have
431  // only muon energy in them, And those that have ha-
432  // -dronic contamination
433  std::vector<int> pureMuPlanes;
434  double totEInCell=0, muEInCell=0;
435 
436  for( int hit = 0; hit < nMuCells; hit++ ){
437 
438  std::vector< cheat::TrackIDE > flshitsID =
439  bt->HitToTrackIDE( muonTrack.Cell(hit) );
440 
441  for( int f = 0; f < (int)flshitsID.size(); f++){
442  totEnergy += flshitsID[f].energy;
443  totEInCell = flshitsID[f].energy;
444 
445  if( primMuTrackIDsSet.find( flshitsID[f].trackID) != primMuTrackIDsSet.end() ){
446  muonEnergy += flshitsID[f].energy;
447  muEInCell = flshitsID[f].energy;
448  }
449 
450  }//end loop over flh hits in cell
451 
452  if( muEInCell == totEInCell ){
453  int plane = muonTrack.Cell(hit)->Plane();
454  if( std::find( pureMuPlanes.begin(), pureMuPlanes.end(), plane ) == pureMuPlanes.end())
455  pureMuPlanes.push_back( plane);
456  }
457  }// end loop over muon hits
458 
459  fTrueSharedEnergy = totEnergy - muonEnergy;
460  if(muonTrack.NCell()>0){
461  fRecoSharedEnergy = fTrkCleanUpAlg->ExtraEOnTrack( muonTrack, *slices[iSli] );
462  std::map< unsigned int, double > dedxinplane = fTrkCleanUpAlg->DeDxInPlane( muonTrack );
463  std::map< unsigned int, double >::iterator iter;
464  for( iter = dedxinplane.begin(); iter != dedxinplane.end(); ++iter){
465  if( std::find( pureMuPlanes.begin(), pureMuPlanes.end(), iter->first) != pureMuPlanes.end() )
466  fdedxmu->Fill( iter->second );
467  else
468  fdedxsh->Fill( iter->second );
469  }
470 
471  }
472  float totMichelEnergy = 0;
473  float totSliceEnergy = 0;
474 
475 
476  fTotMuonEnergy = 0;
477  fTotHadEnergy = 0;
478  fMuFracRemaining = 0;
479  fHadFracRemoved = 0;
480 
481  int nCells = slices[iSli]->NCell();
482  for( int c = 0; c < nCells; c++){
483 
484  if(!bt->IsNoise( slices[iSli]->Cell(c) ) ){
485 
486  muonEnergy = 0;
487  totEnergy = 0;
488  float hadEnergy = 0;
489  float michelE = 0;
490  bool hasMuon = false;
491 
492  std::vector< cheat::TrackIDE> flshitsID =
493  bt->HitToTrackIDE( slices[iSli]->Cell(c) );
494 
495  // Check if this is a muon hit
496 
497  for( int f = 0; f < (int)flshitsID.size(); f++){
498  totEnergy += flshitsID[f].energy;
499 
500  const int testId = flshitsID[f].trackID;
501  if( primMuTrackIDsSet.find( testId ) != primMuTrackIDsSet.end() ){
502  muonEnergy += flshitsID[f].energy;
503  hasMuon = true;
504  }
505 
506  if( flshitsID[f].energy == michelElID)
507  michelE += flshitsID[f].energy;
508 
509  }// end loop over flshits in this cellhit
510 
511 
512  totEnergy = totEnergy - michelE;
513  totMichelEnergy += michelE;
514  totSliceEnergy += totEnergy;
515 
516  // if(hasMuon)
517  // fTrueSharedEnergy += totEnergy - muonEnergy;
518 
519  if( totEnergy == muonEnergy ){
520  rb::RecoHit rHit = slices[iSli]->RecoHit(c);
521  if( rHit.IsCalibrated() ){
522  fMuonHitMip.push_back( rHit.MIP() );
523  fMuonHitGeV.push_back( rHit.GeV() );
524  }
525  }
526  // See if this cell has a hit in the mrcc-cellhits.
527 
528  bool foundHit = false;
529  rb::CellHit matchedHit;
530 
531  int nMrccHits = mrccallHits.size();
532  for( int k = 0; k < nMrccHits; k++){
533  if( slices[iSli]->Cell(c)->Plane() == mrccallHits[k]->Plane() &&
534  slices[iSli]->Cell(c)->Cell() == mrccallHits[k]->Cell() &&
535  slices[iSli]->Cell(c)->TDC() == mrccallHits[k]->TDC() ){
536  foundHit = true;
537  matchedHit = *(mrccallHits[k]);
538  break;
539  }
540  }// end loop over mrcc cellhits
541 
542  float muadcRatio = 0;
543 
544  if(hasMuon){
545 
546  if(foundHit)
547  muadcRatio = (float)( slices[iSli]->Cell(c)->ADC(0)
548  - matchedHit.ADC(0) )/slices[iSli]->Cell(c)->ADC(0);
549  else
550  muadcRatio = 1;
551 
552  fMuFracRemaining += (muonEnergy - (totEnergy * muadcRatio));
553  fTotMuonEnergy += muonEnergy;
554  }
555 
556  hadEnergy = totEnergy - muonEnergy;
557 
558  if(hadEnergy != 0){
559  float adcRatio = 0;
560  if(foundHit)
561  adcRatio = (float)matchedHit.ADC(0)/slices[iSli]->Cell(c)->ADC(0);
562  else
563  adcRatio = 0;
564  fHadFracRemoved += hadEnergy - (totEnergy * adcRatio);
565  fTotHadEnergy += hadEnergy;
566  }
567 
568  }// cellhit is not noise
569  }// end loop over slice cellhits
570 
571  std::cout<<"truesharedE "<<fTrueSharedEnergy<<" recosharedE "<<fRecoSharedEnergy<<"\n";
572  if( totMichelEnergy/totSliceEnergy >= 0.9 || fTotMuonEnergy==0)
573  continue;
574 
575  if(fTotMuonEnergy != 0 )
577 
578  if(fTotHadEnergy != 0 )
580 
581  std::cout<<" tot mu energy "<<fTotMuonEnergy<<" tot had energy "<<fTotHadEnergy<<std::endl;
582  std::cout<<" mu rem frac "<<fMuFracRemaining<<" had frac rem "<<fHadFracRemoved<<std::endl;
583 
584 
585  fTrackEff = bt->HitCollectionEfficiency( primMuTrackIDsSet, muonTrack.AllCells(), allHits, geo::kXorY);
586  fTrackPur = bt->HitCollectionPurity( primMuTrackIDsSet, muonTrack.AllCells(), NULL, NULL, true );
587 
588  }// interaction is CC
589 
590  fSliceTree->Fill();
591 
592  }// neutrino is muon type
593  }// can match to a neutrino interaction
594  }// non-noise slices
595  }// end loop over slices
596 
597 
598 
599  }// End of analyser
600 
601 
602  /*********************************************************************************************/
603  ////////////////////////////////////////////////////////////////////////
604 
606  {
607  int nbins = fdedxmu->GetNbinsX();
608  float integralMu = fdedxmu->Integral(1, nbins+1);
609  float integralHad = fdedxsh->Integral(1, nbins+1);
610 
611  for( int ibin = 1; ibin <= nbins; ++ibin){
612  float leftIntegralMu = fdedxmu->Integral(1, ibin);
613  float rightIntegralMu = fdedxmu->Integral(ibin, nbins+1);
614  float leftIntegralHad = fdedxsh->Integral(1, ibin);
615  float rightIntegralHad = fdedxsh->Integral(ibin, nbins+1);
616 
617 
618  float purHighCut = 0;
619  if( (leftIntegralMu + leftIntegralHad) != 0)
620  purHighCut = 1 - ( leftIntegralHad/ (leftIntegralMu+leftIntegralHad) );
621 
622  float purLowCut = 0;
623  if( (rightIntegralMu + rightIntegralHad) != 0)
624  purLowCut = 1 - ( rightIntegralHad/ rightIntegralMu+rightIntegralHad);
625 
626  float effHighCutMu = leftIntegralMu / integralMu;
627  float effLowCutMu = rightIntegralMu / integralMu;
628  float effHighCutHad = rightIntegralHad / integralHad;
629  float effLowCutHad = leftIntegralHad / integralHad;
630 
631  float x = fdedxmu->GetBinCenter(ibin);
632  fpurh->Fill( x, purHighCut);
633  fpurl->Fill( x, purLowCut);
634  feffmuh->Fill( x, effHighCutMu);
635  feffmul->Fill( x, effLowCutMu);
636  feffhadh->Fill( x, effHighCutHad);
637  feffhadl->Fill( x, effLowCutHad);
638  }// end loop over bins
639  }
640 
641  /*********************************************************************************************/
642  ////////////////////////////////////////////////////////////////////////
643 
645 
646 }// End of namespace
647 
648 
649 /*********************************************************************************************/
650 ////////////////////////////////////////////////////////////////////////
art::ServiceHandle< geo::Geometry > fgeom
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
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
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::vector< float > fMuonHitGeV
void abs(TH1 *hist)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
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...
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
const XML_Char * s
Definition: expat.h:262
const int nbins
Definition: cellShifts.C:15
float ExtraEOnTrack(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Calculates the total hadronic energy on the provided track based on the assumption of MIP...
length
Definition: demo0.py:21
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TrackCleanUpAlg * fTrkCleanUpAlg
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
std::map< unsigned int, double > DeDxInPlane(const rb::Track &muonTrack)
A function to calculate the energy, dE/dx and tracklength of a track through a plane.
double energy
Definition: plottest35.C:25
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
void reconfigure(const fhicl::ParameterSet &p)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
size_type size() const
Definition: PtrVector.h:308
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
void analyze(art::Event const &evt)
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
std::vector< float > fMuonHitMip
float MIP() const
Definition: RecoHit.cxx:58
::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
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
fhicl::ParameterSet fTrkCleanUpPSet
Definition: event.h:1
Event generator information.
Definition: MCNeutrino.h:18
EventID id() const
Definition: Event.h:56
MuonRemoveAna(fhicl::ParameterSet const &pset)
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const