plot_timing_resolution.C
Go to the documentation of this file.
1 #include <vector>
2 #include <algorithm>
3 #include <utility>
4 #include <map>
5 #include <iterator>
6 #include <iostream>
7 #include <stdio.h>
8 #include <cmath>
9 #include "TTree.h"
10 #include "TChain.h"
11 #include "TH1F.h"
12 #include "TH2F.h"
13 #include "TF1.h"
14 #include "TFile.h"
15 #include "TMath.h"
16 #include "TStyle.h"
17 #include "TColor.h"
18 #include "TCanvas.h"
19 #include "TProfile.h"
20 #include "TProfile2D.h"
21 #include "TGraph.h"
22 
23 #include <stdint.h>
24 #include <sstream>
25 #include <fstream>
26 #include <string>
27 
28 //Template macro for calculating timing resolution as a function of PE. This works by calculating time differences between hits on a
29 //thoughgoing cosmic track within a single DCM. This macro expects ntuple input from the TimingCalibration module. Track quality cuts are applied there
30 //and in the further upstream step FiberCalibration. Corrections for time-of-flight and distance to the readout are done to the hits within TimingCalibration.
31 //This macro can be adapted to look at resolution in different views
32 
33 using namespace std;
34 
35 //fit function for the timing resolution. Additional fit formats can be explored
36 Double_t twfit(Double_t *v, Double_t *par)
37 {
38  Double_t x = v[0];
39  Double_t fitval = par[0]/(par[1]+pow(x,par[2])) + par[3];
40 
41  return fitval;
42 }
43 
44 
46 
47  TH1::SetDefaultSumw2();
48  TH2::SetDefaultSumw2();
49  TProfile::SetDefaultSumw2();
50  TProfile2D::SetDefaultSumw2();
51 
52  TChain *timewalk = new TChain("timingcal/TimingCalib");
53  timewalk->Add("PATH_TO_FILES");
54 
55 
56  vector<float> *dcmhits = 0;
57  vector<float> *dcmpe = 0;
58  vector<float> *dcmstime = 0;
59 
60  //turn off all branches, only use the ones we care about
61  //timewalk->SetBranchStatus("*", 0);
62  //turn back on branches of interest
63  timewalk->SetBranchAddress("dcmhits",&dcmhits);
64  timewalk->SetBranchAddress("dcmpe", &dcmpe);
65  timewalk->SetBranchAddress("dcmstime",&dcmstime);
66 
67  TH2F* fHistoDiffSPE = new TH2F("fHistoDiffSPESingle","Corrected hit time difference vs PE;pe;diff (ns)",84,0.0,2100.0,2000,-3000.0,3000.0);
68 
69 
70  unsigned int entries = timewalk->GetEntries();
71  std::cout<<entries<<std::endl;
72  //if (entries > 900000) entries = 900000;
73 
74 
75  for(unsigned int i=0; i<entries; ++i){
76  timewalk->GetEntry(i);
77  if(i%1000 == 0) cout<<"through "<<i<<" of "<<entries<<" events."<<endl;
78  unsigned int count=0;
79  for (unsigned int m=0; m<dcmhits->size(); ++m){
80  unsigned int dhit = (*dcmhits)[m];
81 
82  for (unsigned int j=count; j<count+dhit; ++j){
83  for (unsigned int k=count; k<count+dhit; ++k){
84  if (j==k) continue;
85  float timea = (*dcmstime)[j];
86  float timeb = (*dcmstime)[k];
87  float pea = (*dcmpe)[j];
88  float peb = (*dcmpe)[k];
89  int bina = ((int)pea / 25) + 1;
90  int binb = ((int)peb / 25) + 1;
91  if (bina != binb) continue;
92  fHistoDiffSPE->Fill(pea, timea - timeb);
93  }
94  }
95  count += dhit;
96  }
97  }
98 
99  //optional second part, loop over bins to calculate rms of each bin
100  TH1F* fHistoTemp[64];
101  TH1F* fHistoFit = new TH1F("fHistoFit","Expected Time - TNS from fit vs pe;pe;mean diff (ns)",64,0.0,1600.0);
102  TH1F* fHistoRMS = new TH1F("fHistoRMS","Timing resolution vs pe;pe;resolution (ns)",64,0.0,1600.0);
103 
104  double yvec[64];
105  double xvec[64];
106  for (unsigned int i=1; i<65; ++i){
107  TString myname;
108  myname.Form("%d",i);
109  fHistoTemp[i-1] = new TH1F("fHistoTemp"+myname,"temp",2000,-3000.0,3000.0);
110  double max = 0.0;
111  double maxbin = 150.0;
112  double t = fHistoDiffSPE->GetBinContent(i,0);
113  fHistoTemp[i-1]->SetBinContent(0,t);
114  t = fHistoDiffSPE->GetBinContent(i,2001);
115  fHistoTemp[i-1]->SetBinContent(2001,t);
116  for(unsigned int j=1; j<2001; ++j){
117  double tot = fHistoDiffSPE->GetBinContent(i,j);
118  fHistoTemp[i-1]->SetBinContent(j,tot);
119  if (tot > max){
120  max = tot;
121  maxbin = j;
122  }
123  }
124  //truncate top and bottom 1% from fit
125  double areaM = fHistoTemp[i-1]->Integral(0,2001);
126  fHistoTemp[i-1]->Scale(100.0/areaM);
127  int low = 0;
128  int high = 0;
129  for(unsigned int j=1; j<2001; ++j){
130  Double_t area = fHistoTemp[i-1]->Integral(0,j);
131  if (area >=1.0){
132  low = j;
133  break;
134  }
135  }
136  for(unsigned int j=2000; j>1; --j){
137  Double_t area = fHistoTemp[i-1]->Integral(j,2001);
138  if (area >=1.0){
139  high = j;
140  break;
141  }
142  }
143  fHistoTemp[i-1]->Scale(areaM/100.0);
144  double tot = 0;
145  double avg = 0;
146  double avgsqr = 0;
147  double rms = 0;
148  std::cout<<" i: "<<i<<" low: "<<low<<" high: "<<high<<std::endl;
149  //calculate rms
150  for(int ix = 1; ix < fHistoTemp[i-1]->GetNbinsX()+1; ++ix){
151 
152  if ((ix >= low) && (ix <= high)){
153  const double y = fHistoTemp[i-1]->GetXaxis()->GetBinCenter(ix);
154  const double w = fHistoTemp[i-1]->GetBinContent(ix);
155  avg += w*y;
156  avgsqr += w*y*y;
157  tot += w;
158  }
159  }
160  if (tot > 0){
161  avg /= tot;
162  avgsqr /= tot;
163  //factor of sqrt(2) since each time difference is counted twice (+ and -)
164  rms = sqrt(avgsqr-avg*avg)/sqrt(2.);
165  }
166  delete fHistoTemp[i-1];
167  fHistoFit->SetBinContent(i,avg);
168  fHistoFit->SetBinError(i,rms);
169  fHistoRMS->SetBinContent(i,rms);
170  }
171 
172  TCanvas* c1 = new TCanvas("resolutionfit","resolutionfit",1);
173  TF1* fitfunc = new TF1("tw_first",twfit,5.0,1600.0,4);
174  fHistoRMS->Fit("tw_first","RW");
175  fHistoRMS->SetMarkerStyle(20);
176  fHistoRMS->SetTitle("");
177  fHistoRMS->GetXaxis()->SetTitle("photoelectrons");
178  fHistoRMS->SetMinimum(0.0);
179  fHistoRMS->Draw("HIST P");
180  fitfunc->Draw("same");
181 
182 
183  TFile f("resolution_nd_data_fixspeed_sp.root","recreate");
184  fHistoDiffSPE->Write();
185  c1->Write();
186  f.Close();
187 }
T max(const caf::Proxy< T > &a, T b)
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Int_t par
Definition: SimpleIterate.C:24
void plot_timing_resolution()
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
c1
Definition: demo5.py:24
Float_t w
Definition: plot.C:20
Double_t twfit(Double_t *v, Double_t *par)