BetheBlochFit_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file BetheBlochFit_module.cc
3 // \brief Bethe-Bloch fit to muondEdx output
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
14 #include "fhiclcpp/ParameterSet.h"
15 
16 #include <iostream>
17 
18 #include "TF1.h"
19 #include "TFile.h"
20 #include "TGeoManager.h"
21 #include "TGraph.h"
22 #include "TH2.h"
23 #include "TProfile.h"
24 
28 
30 #include "Geometry/Geometry.h"
31 
32 
33 namespace calib
34 {
36  {
37  public:
38  explicit BetheBlochFit(const fhicl::ParameterSet& pset);
39  virtual ~BetheBlochFit();
40 
41  void analyze(const art::Event& evt) override;
42 
43  void beginJob() override;
44  void endJob() override;
45  void beginRun(const art::Run& run) override;
46  void endRun(const art::Run& run) override;
47 
48  protected:
49  TH1* ProfileMPV(TH2* h) const;
50 
51  double fWindowStart;
52  double fWindowLength;
53 
54  TH2* fdEdx;
55  TH2* fdEdx_true;
56  };
57 
58  //......................................................................
60  : EDAnalyzer(pset)
61  {
62  // TODO: get these from pset
63  fWindowStart = 100.;
64  fWindowLength = 100.;
65  }
66 
67  //......................................................................
69  {
70  }
71 
72  //......................................................................
74  {
75  }
76 
77  //......................................................................
79  {
81  fdEdx = tfs->make<TH2D>("dEdx", ";Distance to track end (cm);PECorr / cm",
82  50, 0, 500, 100, 0, 200);
83  fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm",
84  50, 0, 500, 100, 0, 5);
85  }
86 
87  //......................................................................
89  {
90  }
91 
92  //......................................................................
94  {
95  art::Handle<TH2D> dE, dEtrue;
96  run.getByLabel("muondedx", "dEdx", dE);
97  run.getByLabel("muondedx", "dEdxTrue", dEtrue);
98 
99  fdEdx->Add(&*dE);
100  fdEdx_true->Add(&*dEtrue);
101  }
102 
103  //......................................................................
104  TH1* BetheBlochFit::ProfileMPV(TH2* h) const
105  {
106  TH1* ret = new TH1F("", "", h->GetNbinsX(),
107  h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
108  for(int x = 0; x < h->GetNbinsX()+2; ++x){
109  double maxVal = -1e10;
110  double maxy = 0;
111  for(int y = 1; y < h->GetNbinsY()+1; ++y){
112  const double val = h->GetBinContent(x, y);
113  if(val > maxVal){
114  maxVal = val;
115  maxy = h->GetYaxis()->GetBinCenter(y);
116  }
117  }
118  TH1* proj = h->ProjectionY("", x, x);
119  proj->Fit("gaus", "0Q", "", maxy-0.2, maxy+0.2);
120  TF1* func = proj->GetFunction("gaus");
121  if(func){
122  ret->SetBinContent(x, func->GetParameter(1));
123  ret->SetBinError(x, func->GetParError(1));
124  }
125  }
126 
127  return ret;
128  }
129 
130  //......................................................................
132  {
134 
135  // This differs from soup by about 0.2%
136  // Representative of the whole detector
137  TGeoVolume* block = geom->ROOTGeoManager()->FindVolumeFast("vBlock");
138  const double avgDensity = 1e3*block->Weight()/block->Capacity();
139 
140  // const double avgDensity =
141  // geom->ROOTGeoManager()->GetMaterial("Soup")->GetDensity();
142 
143  const double scintDensity =
144  geom->ROOTGeoManager()->GetMaterial("Scintillator")->GetDensity();
145 
146  std::cout << "Detector average (soup) density is " << avgDensity << std::endl;
147  std::cout << "Scintillator density is " << scintDensity << std::endl;
148 
149  TProfile* prof = fdEdx->ProfileX();
150 
153 
154  EnergyLossVsDistance el_ana(bb_ana, scintDensity, avgDensity);
155  EnergyLossVsDistance el_tab(bb_tab, scintDensity, avgDensity);
156 
157  TF1* f_ana = el_ana.GetTF1();
158  prof->Fit(f_ana, "", "", 50, 500);
159  TF1* f_tab = el_tab.GetTF1();
160  prof->Fit(f_tab, "", "", 50, 500);
161 
162  // Average seen by muon travelling along z-axis
163  const double Lz = (geom->Plane(0)->Cell(4)->HalfD()+geom->Plane(1)->Cell(4)->HalfD())/2;
164 
165  const double pePerMeV = f_ana->GetParameter(0);
166 
167  std::cout << "1 MeV = " << pePerMeV << " PE -> 1 PE = " << 1/pePerMeV << " MeV" << std::endl;
168 
169  const double MeVPerCmPerMIP = bb_ana->MIP()*scintDensity;
170 
171  std::cout << "1 MIP = " << MeVPerCmPerMIP*Lz << " MeV [MeVScale]" << std::endl;
172 
173  const double pePerMIPPerCm = pePerMeV*MeVPerCmPerMIP;
174 
175  std::cout << "1 MIP = " << pePerMIPPerCm*Lz << " PE -> 1 PE = " << 1/pePerMIPPerCm/Lz << " MIP [MIPScale]" << std::endl;
176 
177  // Subvert the file service
179  tfs->file().cd("bb");
180  f_ana->Draw("c");
181  f_ana->GetHistogram()->Write("fit_ana");
182  f_tab->Draw("c");
183  f_tab->GetHistogram()->Write("fit_tab");
184 
185  el_ana.GetEnergyLossGraph()->Write("el_ana");
186  el_tab.GetEnergyLossGraph()->Write("el_tab");
187 
188  bb_tab->GetdEdxGraph(scintDensity)->Write("dEdx_curve");
189  bb_ana->GetdEdxGraph(scintDensity)->Write("dEdx_curve_ana");
190  bb_ana->GetMPVGraph(scintDensity, Lz)->Write("mpv_curve_ana");
191  bb_ana->GetMPVGraph(scintDensity, Lz*2)->Write("mpv_curve_ana2");
192  bb_ana->GetMPVGraph(scintDensity, Lz/2)->Write("mpv_curve_ana3");
193 
194  ProfileMPV(fdEdx)->Write("prof_mpv");
195  }
196 
198 
199 } // namespace
200 
void endRun(const art::Run &run) override
double HalfD() const
Definition: CellGeo.cxx:205
virtual double MIP() const =0
double maxy
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
TGeoManager * ROOTGeoManager() const
TFile & file() const
Definition: TFileService.h:63
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
double dE
Definition: Run.h:31
void analyze(const art::Event &evt) override
TH1 * ProfileMPV(TH2 *h) const
CDPStorage service.
block
print "ROW IS " print row
Definition: elec2geo.py:31
TGraph * GetMPVGraph(double localDensity, double dist) const
Definition: BetheBloch.cxx:45
TGraph * GetdEdxGraph(double localDensity) const
Definition: BetheBloch.cxx:33
int evt
double func(double x, double y)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void beginRun(const art::Run &run) override
Definition: run.py:1
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
Float_t proj
Definition: plot.C:35
void geom(int which=0)
Definition: geom.C:163
Integrate energy losses to calculate deposit near a track end.
Encapsulate the geometry of one entire detector (near, far, ndos)
BetheBlochFit(const fhicl::ParameterSet &pset)