ParticleCorrections_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ParticleCorrections_module.cc
3 // \brief Make plots of threshold effect from PCHits files
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 //C++ and C:
8 #include <cmath>
9 #include <vector>
10 
11 //Framework:
16 
17 //NovaSoft:
20 #include "Geometry/Geometry.h"
22 #include "MCCheater/BackTracker.h"
23 #include "Simulation/FLSHitList.h"
24 #include "Simulation/Particle.h"
25 
26 // ROOT
27 #include "TH2.h"
28 
29 namespace calib
30 {
31  /// Make plots of threshold effect from PCHits files
33  {
34  public:
35  explicit ParticleCorrections(const fhicl::ParameterSet& pset);
36  virtual ~ParticleCorrections();
37 
38  virtual void beginJob() override;
39 
40  virtual void analyze(const art::Event& evt) override;
41 
42  virtual void reconfigure(const fhicl::ParameterSet& pset);
43 
44  protected:
47 
48  TH1* fT;
49  };
50 
51  //---------------------------------------------------------------------
53  : EDAnalyzer(pset)
54  {
55  reconfigure(pset);
56  }
57 
58  //---------------------------------------------------------------------
60  {
61  }
62 
63  //---------------------------------------------------------------------
65  {
66  }
67 
68  //---------------------------------------------------------------------
70  {
72 
73  fDepE = tfs->make<TH1F>("depe", "True e^{-};True energy deposit (GeV);Cells", 1000, 0, .1);
74  fDepPos = tfs->make<TH1F>("deppos", "True e^{+};True energy deposit (GeV);Cells", 1000, 0, .1);
75  fDepMum = tfs->make<TH1F>("depmum", "True #mu^{-};True energy deposit (GeV);Cells", 1000, 0, .1);
76  fDepMup = tfs->make<TH1F>("depmup", "True #mu^{+};True energy deposit (GeV);Cells", 1000, 0, .1);
77  fDepProt = tfs->make<TH1F>("depprot", "True p;True energy deposit (GeV);Cells", 1000, 0, .1);
78  fDepPip = tfs->make<TH1F>("deppip", "True #pi^{+};True energy deposit (GeV);Cells", 1000, 0, .1);
79  fDepPim = tfs->make<TH1F>("deppim", "True #pi^{-};True energy deposit (GeV);Cells", 1000, 0, .1);
80 
81  fDepE2 = tfs->make<TH1F>("depe2", "True e^{-};True energy deposit (GeV);Cells", 1000, 0, .1);
82  fDepPos2 = tfs->make<TH1F>("deppos2", "True e^{+};True energy deposit (GeV);Cells", 1000, 0, .1);
83  fDepMum2 = tfs->make<TH1F>("depmum2", "True #mu^{-};True energy deposit (GeV);Cells", 1000, 0, .1);
84  fDepMup2 = tfs->make<TH1F>("depmup2", "True #mu^{+};True energy deposit (GeV);Cells", 1000, 0, .1);
85  fDepProt2 = tfs->make<TH1F>("depprot2", "True p;True energy deposit (GeV);Cells", 1000, 0, .1);
86  fDepPip2 = tfs->make<TH1F>("deppip2", "True #pi^{+};True energy deposit (GeV);Cells", 1000, 0, .1);
87  fDepPim2 = tfs->make<TH1F>("deppim2", "True #pi^{-};True energy deposit (GeV);Cells", 1000, 0, .1);
88 
89  fT = tfs->make<TH1F>("t", "", 500, 0, 500e3);
90  }
91 
92  //---------------------------------------------------------------------
94  {
96  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
97 
99  evt.getByLabel("geantgen", flsls);
100 
101  assert(flsls->size() == 1);
102 
103  const std::vector<sim::FLSHit>& flss = (*flsls)[0].fHits;
104 
105  std::map<int, std::map<geo::OfflineChan, double>> deps;
106 
107  for(const sim::FLSHit& fls: flss){
108  // Worried about Michels and pion decay
109  /*
110  if(!pnav.IsPrimary(fls.GetTrackID())){
111  std::cout << "Rejecting non-primary " << fls.GetTrackID() << " " << fls.GetPDG() << std::endl;
112  continue;
113  }
114  */
115 
116  const sim::Particle* part = bt->TrackIDToParticle(fls.GetTrackID());
117  // Only want the main particle
118  if(part->Process() == "Decay") continue;
119 
120  if(!part){
121  std::cout << "No particle" << std::endl;
122  continue;
123  }
124  // if(part->PdgCode() != fls.GetPDG()) continue;
125 
126  const int evePdg = bt->TrackIDToParticle(pnav.EveId(fls.GetTrackID()))->PdgCode();
127 
128  // Don't worry about Birks for now
129  deps[evePdg/*fls.GetPDG()*/][geo::OfflineChan(fls.GetPlaneID(),
130  fls.GetCellID())] += fls.GetEdep();
131 
132  // Make sure we don't have Michels etc
133  if(evePdg/*fls.GetPDG()*/ == 11) fT->Fill(fls.GetEntryT());
134  }
135 
136  for(auto it: deps[ +11]) fDepE ->Fill(it.second);
137  for(auto it: deps[ -11]) fDepPos ->Fill(it.second);
138  for(auto it: deps[ +13]) fDepMum ->Fill(it.second);
139  for(auto it: deps[ -13]) fDepMup ->Fill(it.second);
140  for(auto it: deps[2212]) fDepProt->Fill(it.second);
141  for(auto it: deps[+211]) fDepPip ->Fill(it.second);
142  for(auto it: deps[-211]) fDepPim ->Fill(it.second);
143 
144  for(auto it: deps[ +11]) fDepE2 ->Fill(it.second, it.second);
145  for(auto it: deps[ -11]) fDepPos2 ->Fill(it.second, it.second);
146  for(auto it: deps[ +13]) fDepMum2 ->Fill(it.second, it.second);
147  for(auto it: deps[ -13]) fDepMup2 ->Fill(it.second, it.second);
148  for(auto it: deps[2212]) fDepProt2->Fill(it.second, it.second);
149  for(auto it: deps[+211]) fDepPip2 ->Fill(it.second, it.second);
150  for(auto it: deps[-211]) fDepPim2 ->Fill(it.second, it.second);
151  }
152 
154 
155 } //end namespace
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
set< int >::iterator it
ParticleCorrections(const fhicl::ParameterSet &pset)
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
virtual void reconfigure(const fhicl::ParameterSet &pset)
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
Definition: MCParticle.h:214
TString part[npart]
Definition: Style.C:32
CDPStorage service.
int evt
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Make plots of threshold effect from PCHits files.
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
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
A (plane, cell) pair.
Definition: OfflineChan.h:17
assert(nhit_max >=nhit_nbins)
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual void analyze(const art::Event &evt) override
int EveId(const int trackID) const