AbsCalib_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: AbsCalib
3 // Module Type: analyzer
4 // File: AbsCalib_module.cc
5 //
6 // Generated at Fri May 10 14:12:00 2013 by Abigail Waldron using artmod
7 // from cetpkgsupport v1_02_00.
8 ////////////////////////////////////////////////////////////////////////
9 
17 #include "fhiclcpp/ParameterSet.h"
18 #include <TH2D.h>
19 #include <TH1D.h>
20 #include <TF1.h>
23 #include <iostream>
24 
25 namespace calib{
26  class AbsCalib;
27 }
28 
29 namespace calib{
30  class AbsCalib : public art::EDAnalyzer {
31  public:
32  explicit AbsCalib(fhicl::ParameterSet const & p);
33  virtual ~AbsCalib();
34 
35  void analyze(art::Event const & e) override;
36 
37  void beginJob() override;
38  void endJob() override;
39  void beginRun(art::Run const & r) override;
40  void endRun(art::Run const & r) override;
41 
42 
43 
44  private:
45 
46  double fitMEU(TH1D*& proj, double& sigma, bool isTruth=false);
47 
48  // Declare member data here.
49  double fWindowStart;
50  double fWindowLength;
51 
54 
55  TH2D* fdEdx;
56  TH2D* fdEdx_minus;
57  TH2D* fdEdx_plus;
58  TH2D* fdEdx_true;
59  TH1D* fWindow;
61  TH1D* fWindow_plus;
62  TH1D* fWindow_true;
63  };
64 
65 
67  : EDAnalyzer(p)
68  {
69  fWindowStart = 100.;
70  fWindowLength = 100.;
71  }
72 
74  {
75  }
76 
78  {
79  }
80 
82  {
84  fdEdx = tfs->make<TH2D>("dEdx","; ; ;",
85  50,0.,500.,200,0.,200.);
86  fdEdx_minus = tfs->make<TH2D>("dEdx_minus","; ; ;",
87  50,0.,500.,200,0.,200.);
88  fdEdx_plus = tfs->make<TH2D>("dEdx_plus","; ; ;",
89  50,0.,500.,200,0.,200.);
90  fdEdx_true = tfs->make<TH2D>("dEdx_true","; ; ;",
91  50,0.,500.,200,0.,5.);
92  fWindow = tfs->make<TH1D>("window","window",
93  200,0.,200.);
94  fWindow_minus = tfs->make<TH1D>("window_minus","window_minus",
95  200,0.,200.);
96  fWindow_plus = tfs->make<TH1D>("window_plus","window_plus",
97  200,0.,200.);
98  fWindow_true = tfs->make<TH1D>("window_true","window_true",
99  200,0.,5.);
100  }
101 
103  {
104  // For now calculate MEU per job, maybe later move to run
105  fMEU = this->fitMEU(fWindow,fMEU_error);
109 
110 
111  double sys_err = fabs(fMEU_plus - fMEU);
112  sys_err = (fabs(fMEU - fMEU_minus) > sys_err ) ? fabs(fMEU - fMEU_minus) : sys_err;
113 
114  std::cout << "MEU = " << fMEU << " +/-" << fMEU_error << "(stat) PECorr/cm"
115  << " +/-" << sys_err << "(sys) PECorr/cm"
116  << std::endl;
117 
118 
119  std::cout << "True MEU = " << fMEU_true << " +/-" << fMEU_error_true << "(stat) MeV/cm"
120  << std::endl;
121 
122 
123  }
124 
125  double AbsCalib::fitMEU(TH1D*& proj, double& sigma, bool isTruth){
126 
127  Double_t parameters[3];
128  parameters[0] = proj->GetBinContent(proj->GetMaximumBin());
129  parameters[1] = proj->GetBinCenter(proj->GetMaximumBin());
130  parameters[2] = proj->GetRMS();
131 
132 
133  // only care about finding peak position, use gaussian truncated about mode
134  Double_t bound_size = (isTruth) ? 0.2 : 0.5; // narrower for truth, *not* MC reco... don't mess this up.
135 
136  Double_t lowerBound = std::max(0.,parameters[1] - bound_size*proj->GetRMS());
137  Double_t upperBound = parameters[1] + bound_size*proj->GetRMS();
138 
139 
140  TF1* fitfunction = new TF1("fitfunction","[0]*TMath::Gaus(x,[1],[2])",lowerBound,upperBound);
141  fitfunction->SetParameters(parameters);
142 
143 
144  // fitting
145  proj->Fit(fitfunction,"RQ");
146 
147 
148  // retrieve values
149  sigma = fitfunction->GetParError(1);
150  return fitfunction->GetParameter(1);
151  }
152 
153 
155  {
156  }
157 
159  {
160  art::Handle<TH2D> tmp, tmpMinus, tmpPlus, tmpTrue;
161  r.getByLabel("muondedx","dEdx",tmp);
162  r.getByLabel("muondedx","dEdxMinus",tmpMinus);
163  r.getByLabel("muondedx","dEdxPlus",tmpPlus);
164  r.getByLabel("muondedx","dEdxTrue",tmpTrue);
165 
166  for( int i=0; i<fdEdx->GetNbinsX(); ++i ){
167  for( int j=0; j<fdEdx->GetNbinsY(); ++j ){
168  double thing = fdEdx->GetBinContent(i+1,j+1);
169  double thingMinus = fdEdx_minus->GetBinContent(i+1,j+1);
170  double thingPlus = fdEdx_plus->GetBinContent(i+1,j+1);
171  double thingTrue = fdEdx_true->GetBinContent(i+1,j+1);
172  thing += tmp->GetBinContent(i+1,j+1);
173  thingMinus += tmpMinus->GetBinContent(i+1,j+1);
174  thingPlus += tmpPlus->GetBinContent(i+1,j+1);
175  thingTrue += tmpTrue->GetBinContent(i+1,j+1);
176  fdEdx->SetBinContent(i+1,j+1,thing);
177  fdEdx_minus->SetBinContent(i+1,j+1,thingMinus);
178  fdEdx_plus->SetBinContent(i+1,j+1,thingPlus);
179  fdEdx_true->SetBinContent(i+1,j+1,thingTrue);
180  }
181  }
182 
183  // project into track window
184  int startBin = fdEdx->GetXaxis()->FindBin(fWindowStart);
185  int endBin = 1 + fdEdx->GetXaxis()->FindBin(fWindowStart+fWindowLength);
186  fWindow->Add((TH1D*)tmp->ProjectionY("window",startBin,endBin));
187  fWindow_minus->Add((TH1D*)tmpMinus->ProjectionY("window_minus",startBin,endBin));
188  fWindow_plus->Add((TH1D*)tmpPlus->ProjectionY("window_plus",startBin,endBin));
189  fWindow_true->Add((TH1D*)tmpTrue->ProjectionY("window_true",startBin,endBin));
190  }
191 }// end namspace
192 
T max(const caf::Proxy< T > &a, T b)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double fitMEU(TH1D *&proj, double &sigma, bool isTruth=false)
const char * p
Definition: xmltok.h:285
Float_t tmp
Definition: plot.C:36
AbsCalib(fhicl::ParameterSet const &p)
DEFINE_ART_MODULE(TestTMapFile)
Definition: Run.h:31
void endJob() override
void beginJob() override
void beginRun(art::Run const &r) override
CDPStorage service.
void endRun(art::Run const &r) override
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
void analyze(art::Event const &e) override
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
TRandom3 r(0)
Float_t e
Definition: plot.C:35