TestParticleCorrections_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file TestParticleCorrections_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:
18 #include "Calibrator/Calibrator.h"
21 #include "Geometry/Geometry.h"
23 #include "MCCheater/BackTracker.h"
24 #include "Simulation/FLSHitList.h"
25 #include "Simulation/Particle.h"
26 #include "RecoBase/RecoHit.h"
27 
28 // ROOT
29 #include "TFile.h"
30 #include "TGraph.h"
31 #include "TH2.h"
32 
33 namespace calib
34 {
35  /// Make plots of threshold effect from PCHits files
37  {
38  public:
39  explicit TestParticleCorrections(const fhicl::ParameterSet& pset);
40  virtual ~TestParticleCorrections();
41 
42  virtual void beginJob() override;
43 
44  virtual void analyze(const art::Event& evt) override;
45 
46  virtual void reconfigure(const fhicl::ParameterSet& pset);
47 
48  protected:
49  TH1 *fE, *fMu, *fProt, *fPi;
51 
52  TH1* fCell;
53 
54  TGraph *fGE, *fGMu, *fGProt, *fGPi;
55  };
56 
57  //---------------------------------------------------------------------
59  : EDAnalyzer(pset)
60  {
61  reconfigure(pset);
62  }
63 
64  //---------------------------------------------------------------------
66  {
67  }
68 
69  //---------------------------------------------------------------------
71  {
72  }
73 
74  //---------------------------------------------------------------------
76  {
78 
79  fE = tfs->make<TH1F>("e", "True e^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
80  fMu = tfs->make<TH1F>("mu", "True #mu^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
81  fProt = tfs->make<TH1F>("prot", "True p;Energy estimate (GeV);Slices", 400, 0, 4);
82  fPi = tfs->make<TH1F>("pi", "True #pi^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
83 
84 
85  fEPre = tfs->make<TH1F>("e_pre", "True e^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
86  fMuPre = tfs->make<TH1F>("mu_pre", "True #mu^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
87  fProtPre = tfs->make<TH1F>("prot_pre", "True p;Energy estimate (GeV);Slices", 400, 0, 4);
88  fPiPre = tfs->make<TH1F>("pi_pre", "True #pi^{#pm};Energy estimate (GeV);Slices", 400, 0, 4);
89 
90  fE->SetLineColor(kRed);
91  fMu->SetLineColor(kBlue);
92  fProt->SetLineColor(kBlack);
93  fPi->SetLineColor(kMagenta);
94 
95  fEPre->SetLineColor(kRed);
96  fMuPre->SetLineColor(kBlue);
97  fProtPre->SetLineColor(kBlack);
98  fPiPre->SetLineColor(kMagenta);
99 
100  fCell = tfs->make<TH1F>("cell", ";Cell reco/true", 100, 0, 2);
101 
102  TFile* f = new TFile("pcorr_curves.root");
103  fGE = (TGraph*)f->Get("e");
104  fGMu = (TGraph*)f->Get("mu");
105  fGProt = (TGraph*)f->Get("p");
106  fGPi = (TGraph*)f->Get("pi");
107  }
108 
109  //---------------------------------------------------------------------
110  double corr(TGraph* g, double thresh)
111  {
112  for(int i = 0; i < g->GetN(); ++i){
113  double x, y;
114  g->GetPoint(i, x, y);
115 
116  if(x > thresh) return 100/y;
117  }
118 
119  return 1; // Uh oh
120  }
121 
122  //---------------------------------------------------------------------
124  {
125  const double hack = 4095/2170.; // why?
126 
129 
131  evt.getByLabel("slicer", slices);
132 
133  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
134  const art::Ptr<rb::Cluster> slice(slices, sliceIdx);
135  if(slice->IsNoise()) continue;
136 
137  for(unsigned int iCell = 0; iCell < slice->NCell(); ++iCell){
138  const double reco = slice->RecoHit( iCell).GeV()*hack;
139  const std::vector<sim::FLSHit> flss = bt->HitToFLSHit(slice->Cell(iCell));
140  double truth = 0;
141  for(auto f: flss) truth += f.GetEdep();
142 
143  fCell->Fill(reco/truth);
144  }
145 
146  const double est = slice->TotalGeV();
147 
148  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(slice->AllCells());
149 
150  if(parts.empty()) continue;
151  const sim::Particle* part = parts[0];
152 
154  dummy.SetPlane(slice->MinPlane());
155  dummy.SetCell(slice->MinCell(geo::kX));
156  dummy.SetADC(0, 25);
157  dummy.SetPE(cal->GetPE(&dummy));
158 
159  dummy.SetMC();
160 
161  const rb::RecoHit rhit = cal->MakeRecoHit(dummy, slice->MeanX()); // symmetric
162 
163  // Crash when not calibrated. Always should be in MC
164  const double thresh = rhit.GeV()*1000; // MeV
165  // std::cout << thresh << std::endl;
166 
167  // From gdml, correct to total density of soup
168  const double dead = (.85*.859+.15*1.49)/(.85*.859);
169 
170  switch(part->PdgCode()){
171  case -11:
172  case +11:
173  fE->Fill(est*corr(fGE, thresh)*dead*hack);
174  fEPre->Fill(est*hack);
175  break;
176  case -13:
177  case +13:
178  fMu->Fill(est*corr(fGMu, thresh)*dead*hack);
179  fMuPre->Fill(est*hack);
180  break;
181  case 2212:
182  fProt->Fill(est*corr(fGProt, thresh)*dead*hack);
183  fProtPre->Fill(est*hack);
184  break;
185  case -211:
186  case +211:
187  fPi->Fill(est*corr(fGPi, thresh)*dead*hack);
188  fPiPre->Fill(est*hack);
189  break;
190  }
191  }
192  }
193 
195 
196 } //end namespace
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
double corr(TGraph *g, double thresh)
DEFINE_ART_MODULE(TestTMapFile)
virtual void reconfigure(const fhicl::ParameterSet &pset)
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void SetPlane(unsigned short plane)
Definition: CellHit.h:53
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:231
void SetCell(unsigned short cell)
Definition: CellHit.h:52
void SetADC(uint32_t i, int16_t iADC)
Definition: RawDigit.cxx:108
TString part[npart]
Definition: Style.C:32
CDPStorage service.
void SetPE(float pe)
Definition: CellHit.h:55
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
int evt
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TestParticleCorrections(const fhicl::ParameterSet &pset)
void SetMC(bool isMC=true)
Definition: RawDigit.h:106
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
virtual void analyze(const art::Event &evt) override
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
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
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
Make plots of threshold effect from PCHits files.
float GetPE(rb::CellHit const &)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Encapsulate the geometry of one entire detector (near, far, ndos)
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.