Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
drawWaveforms.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TTree.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "TH3.h"
6 #include "TCanvas.h"
7 #include "TAxis.h"
8 
9 double par0 = 55800;
10 double par1 = -5.73e-6;
11 double par2 = -1.03e-10;
13 
14 Double_t myfunction(Double_t *x, Double_t *par)
15 {
16  Float_t xx =x[0];
17  Double_t f = 1;
18  if(xx > par[0])
19  f += (par[1]*(xx-par[0]) + par[2]*pow((xx - par[0]),2));
20 
21  return f;
22 }
23 
24 void drawWaveform(TString bcalEnd = "DS", int layer = 1, int run = 30894)
25 {
26  double lineWidth = 2.0;
27  double textSize = 0.07;
28 
29  bool fall2016 = false;
30 
31  TString figsDir = "./figs/"; figsDir+=run; figsDir+="/";
32 
33  TH2D* hIntegralRatio_UnsaturatedIntegral = new TH2D("IntegralRatio_UnsaturatedIntegral", "; Unsaturated Integral (fADC); Saturated/Unsaturated Integral", 750, 0, 120000, 100, 0.5, 1.05);
34  TH2D* hIntegralRatio_SaturatedIntegral = new TH2D("IntegralRatio_SaturatedIntegral", "; Saturated Integral (fADC); Saturated/Unsaturated Integral", 750, 0, 120000, 100, 0.5, 1.05);
35 
36  TH2D* hIntegralResidual_SaturatedIntegral = new TH2D("IntegralResidual_SaturatedIntegral", "; Saturated Integral (fADC); Integral Residual", 750, 0, 120000, 100, -0.4, 0.4);
37 
38  TF1 *fknown = new TF1("fKnown",myfunction,0,160000,3);
39  fknown->SetParameters(par0, par1, par2);
40 
41  TH2D* hIntegralRatio_UnsaturatedIntegralFull = new TH2D("IntegralRatio_UnsaturatedIntegralFull", "; Unsaturated Full Integral (fADC); Saturated/Unsaturated Full Integral", 500, 0, 200000, 100, 0.5, 1.05);
42 
43  TString fileName = "hd_root_";
44  fileName+=run; fileName+="_saturate.root";
45  TFile *f = TFile::Open(fileName);
46 
47  TH2I* hbcal_waveform_saturate_2D = (TH2I*)f->Get(Form("bcal_saturation/bcal%s_waveform_saturate_layer%d", bcalEnd.Data(), layer));
48  TH2I* hbcal_waveform_2D = (TH2I*)f->Get(Form("bcal_saturation/bcal%s_waveform_noSaturate_layer%d", bcalEnd.Data(), layer));
49  //h->GetYaxis()->SetNdivisions(4);
50  //h->GetXaxis()->SetNdivisions(6);
51  //h->SetLineWidth(lineWidth);
52 
53  TCanvas* aa = new TCanvas("aa", "aa", 600, 400);
54  hbcal_waveform_saturate_2D->Draw("colz");
55  hbcal_waveform_saturate_2D->SetMaximum(4096);
56  hbcal_waveform_saturate_2D->GetXaxis()->SetRangeUser(0,400);
57  aa->Print(Form("%sSaturate_%s_layer%d.pdf",figsDir.Data(),bcalEnd.Data(),layer));
58 
59  TCanvas* bb = new TCanvas("bb", "bb", 600, 400);
60  hbcal_waveform_2D->Draw("colz");
61 
62  TCanvas* cc = new TCanvas("cc", "cc", 600, 400);
63  TH1F* hBase = new TH1F("hBase", ";Samples; fADC", 100, 0, 100);
64  hBase->SetMaximum(6000.);
65  hBase->Draw();
66 
67  int nbins = hbcal_waveform_2D->GetXaxis()->GetNbins();
68  for(int i=0; i<nbins; i++) {
69 
70  TH1I* hbcal_waveform = (TH1I*)hbcal_waveform_2D->ProjectionY(Form("Event_%d",i),i+1,i+1);
71  if(hbcal_waveform->Integral() < 10) continue;
72 
73  hbcal_waveform->SetLineColor(kBlack);
74  hbcal_waveform->SetLineWidth(2);
75  int thresh = 120;
76  int threshCrossBin = 0;
77  for(int j=0; j<hbcal_waveform->GetXaxis()->GetNbins(); j++) {
78  if(hbcal_waveform->GetBinContent(j) > thresh){
79  threshCrossBin = j;
80  break;
81  }
82  }
83 
84  TH1I* hbcal_waveform_temp = (TH1I*)hbcal_waveform->Clone();
85  hbcal_waveform_temp->Reset();
86 
87  // put all waveforms on the same x-axis
88  int nSaturate = 0;
89  int sampleOffset = 20;
90  int binMax = 0;
91  double fADCmax = 0;
92  int nSamples = 0;
93  for(int j=0; j<hbcal_waveform->GetXaxis()->GetNbins(); j++) {
94  double fADC = hbcal_waveform->GetBinContent(j+1);
95  if(fADC > 4096) {
96  fADC = 4096;
97  nSaturate++;
98  }
99  hbcal_waveform_temp->SetBinContent(j+1 - threshCrossBin+sampleOffset, fADC);
100  hbcal_waveform_temp->SetBinError(j+1 - threshCrossBin+sampleOffset, 0.); //sqrt(fADC));
101 
102  // find maximum bin for outlier waveform shapes
103  if(fADC > fADCmax) {
104  binMax = j+1 - threshCrossBin;
105  fADCmax = fADC;
106  }
107 
108  // find how many samples in the readout window
109  if(fADC > thresh)
110  nSamples++;
111  }
112 
113  // remove waveforms with peaks too far after threshold crossing
114  if(binMax > 5)
115  continue;
116 
117  // remove waveforms with too few samples in readout window
118  if(nSamples < 26)
119  continue;
120 
121  cc->cd();
122  hbcal_waveform_temp->Draw("h same");
123 
124  int pedestal = 100;
125  int integralNSA = 26; // Spring 2016 = 26 and Fall 2016 = 13
126  int integralNSB = 1; // Spring 2016 = 1 and Fall 2016 = -2
127  if(fall2016) {
128  integralNSA = 13;
129  integralNSB = -2;
130  }
131 
132  int integralBinMin = sampleOffset - integralNSB;
133  int integralBinMax = sampleOffset + integralNSA;
134 
135  // demonstration with single waveform
136  if(i == 211) {
137  hbcal_waveform_temp->SetTitle("");
138  TH1I *hdemo_scale = (TH1I*)hbcal_waveform_temp->Clone();
139  double scale = 1.5;
140  hdemo_scale->Scale(scale);
141  TH1I *hdemo_saturate = (TH1I*)hdemo_scale->Clone();
142  for(int j=0; j<hdemo_saturate->GetXaxis()->GetNbins(); j++) {
143  if(j+1 < integralBinMin || j+1 >= integralBinMax)
144  hdemo_saturate->SetBinContent(j+1, 0);
145  if(hdemo_saturate->GetBinContent(j+1) > 4095)
146  hdemo_saturate->SetBinContent(j+1, 4096);
147  }
148  hdemo_saturate->SetFillColor(kRed);
149 
150  TCanvas *demo = new TCanvas("demo", "demo", 800, 400);
151  demo->Divide(2,1);
152  demo->cd(1);
153  hbcal_waveform_temp->SetMaximum(6000);
154  hbcal_waveform_temp->Draw();
155  demo->cd(2);
156  hdemo_scale->SetMaximum(6000);
157  hdemo_scale->Draw();
158  hdemo_saturate->Draw("same");
159 
160  demo->Print(Form("%sexampleSaturate_%0.1f.pdf", figsDir.Data(), scale));
161  }
162 
163  double integral_unscaled = 0;
164  double integral_unscaled_full = 0;
165  double integral_unsaturated = 0;
166  double integral_unsaturated_full = 0;
167 
168  // try different scale factors for integral
169  for(int j=0; j<100; j++) {
170  double scale = 1. + j*0.02;
171 
172  // calculate integral
173  double integral = 0;
174  double integral_full = 0;
175  for(int k=integralBinMin; k<hbcal_waveform_temp->GetXaxis()->GetNbins(); k++) {
176  double fADC = scale * (hbcal_waveform_temp->GetBinContent(k) - pedestal);
177  double fADC_truncated = fADC > (4095 - pedestal) ? (4095 - pedestal) : fADC;
178 
179  // only calculate integral over used samples
180  if(k < integralBinMax)
181  integral += fADC_truncated;
182 
183  // test using full integral (maybe relevant for mcsmear?)
184  integral_full += fADC_truncated;
185  }
186  if(j==0) {
187  integral_unscaled = integral;
188  integral_unscaled_full = integral_full;
189  }
190  else {
191  integral_unsaturated = integral_unscaled*scale;
192  integral_unsaturated_full = integral_unscaled_full*scale;
193  }
194 
195  // check for outlier fADC pulses
196  if(integral_unsaturated < 60000 && integral/integral_unsaturated < 0.9)
197  cout<<i<<endl;
198 
199  hIntegralRatio_UnsaturatedIntegral->Fill(integral_unsaturated, integral/integral_unsaturated);
200  hIntegralRatio_SaturatedIntegral->Fill(integral, integral/integral_unsaturated);
201  double integral_residual = (integral/fknown->Eval(integral) - integral_unsaturated)/integral_unsaturated;
202  hIntegralResidual_SaturatedIntegral->Fill(integral, integral_residual);
203 
204  hIntegralRatio_UnsaturatedIntegralFull->Fill(integral_unsaturated_full, integral_full/integral_unsaturated_full);
205  }
206 
207  //cc->Print(Form("figs/bcal_waveform_%03d.ps", i));
208  }
209 
210  double maxXaxis = 1.2e5;
211  if(fall2016)
212  maxXaxis = 10e4;
213 
214  TCanvas* dd = new TCanvas("dd", "dd", 1000, 400);
215  dd->Divide(3,1);
216  dd->cd(1);
217  hIntegralRatio_UnsaturatedIntegral->Draw("colz");
218 
219  TF1 *f1 = new TF1("myfunc",myfunction,0,160000,3);
220  f1->SetParameters(55000,-4e-6,0);
221  f1->SetLineColor(kRed);
222  f1->SetParNames("offset","linear","quad");
223  //hIntegralRatio_UnsaturatedIntegral->Fit(f1, "", "", 30000, 150000);
224  hIntegralRatio_UnsaturatedIntegral->Draw("colz");
225  hIntegralRatio_UnsaturatedIntegral->GetXaxis()->SetRangeUser(0., maxXaxis);
226 
227  dd->cd(2);
228  hIntegralRatio_SaturatedIntegral->Draw("colz");
229  hIntegralRatio_SaturatedIntegral->GetXaxis()->SetRangeUser(0., maxXaxis);
230 
231  TF1 *f2 = new TF1("myfunc",myfunction,0,160000,3);
232  f2->SetParameters(55000,-4e-6,0);
233  if(fall2016) {
234  f2->SetParameters(35000,0,-4e-6);
235  //f2->FixParameter(0, 34000); // Fall 2016
236  //f2->FixParameter(1, 0); // Fall 2016
237  }
238  f2->SetLineColor(kRed);
239  f2->SetParNames("offset","linear","quad");
240  if(fall2016)
241  hIntegralRatio_SaturatedIntegral->Fit(f2, "", "", 30000, 43000); // Fall 2016
242  else
243  hIntegralRatio_SaturatedIntegral->Fit(f2, "", "", 30000, 150000); // Spring 2016
244 
245  // set parameters for plotting residual
246  par0 = f2->GetParameter(0); par0err = f2->GetParError(0);
247  par1 = f2->GetParameter(1); par1err = f2->GetParError(1);
248  par2 = f2->GetParameter(2); par2err = f2->GetParError(2);
249 
250  hIntegralRatio_SaturatedIntegral->Draw("colz");
251 
252  dd->cd(3);
253  hIntegralResidual_SaturatedIntegral->Draw("colz");
254  hIntegralResidual_SaturatedIntegral->GetXaxis()->SetRangeUser(0., maxXaxis);
255 
256  dd->Print(Form("%s%s_layer%d.pdf",figsDir.Data(),bcalEnd.Data(),layer));
257 
258  return;
259 
260 }
261 
262 // main function to loop over BCAL layers and ends
263 void drawWaveforms(int run = 30894) {
264 
265  gSystem->Exec(Form("mkdir -p figs/%d", run));
266 
267  ofstream parameters;
268  string parFileName = Form("./figs/%d/parameters.txt", run);
269  //cout<<parFileName.data()<<endl;
270  parameters.open(parFileName);
271 
272  TH1F *h0 = new TH1F("par0", "Integral Offset Parameter", 8, 0, 8);
273  TH1F *h1 = new TH1F("par1", "Linear Parameter", 8, 0, 8);
274  TH1F *h2 = new TH1F("par2", "Quadratic Parameter", 8, 0, 8);
275 
276  TString bcalEnd[2] = {"US", "DS"};
277  for(int end=0; end<2; end++) {
278  for(int layer=1; layer<=4; layer++) {
279 
280  drawWaveform(bcalEnd[end], layer, run); // first fit to determine parameters
281  drawWaveform(bcalEnd[end], layer, run); // then fill residual histogram with fitted parameters
282 
283  parameters<<end<<" "<<layer<<" "<<par0<<" "<<par1<<" "<<par2<<endl;
284  h0->SetBinContent(end*4 + layer, par0);
285  h0->SetBinError(end*4 + layer, par0err);
286  h1->SetBinContent(end*4 + layer, par1);
287  h1->SetBinError(end*4 + layer, par1err);
288  h2->SetBinContent(end*4 + layer, par2);
289  h2->SetBinError(end*4 + layer, par2err);
290  }
291  }
292 
293  TCanvas *pars = new TCanvas("pars", "pars", 1200, 500);
294  pars->Divide(3,1);
295  pars->cd(1);
296  h0->Draw();
297  h0->Fit("pol0");
298  pars->cd(2);
299  h1->Draw();
300  h1->SetMinimum(-8e-6); h1->SetMaximum(-3e-6);
301  h1->Fit("pol0");
302  pars->cd(3);
303  h2->Draw();
304  h2->SetMinimum(-1.6e-10); h2->SetMaximum(-1.0e-11);
305  h2->Fit("pol0");
306  TString parFigName = "./figs/";
307  parFigName+=run; parFigName+="/pars.pdf";
308  pars->Print(parFigName);
309 
310  parameters.close();
311 
312  return;
313 }
314 
double par1
Definition: drawWaveforms.C:10
void drawWaveforms(int run=30894)
Int_t layer
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
double par2err
Definition: drawWaveforms.C:12
Double_t myfunction(Double_t *x, Double_t *par)
Definition: drawWaveforms.C:14
TF1 * f
Definition: FitGains.C:21
TFile f1(fnam)
TEllipse * e
void drawWaveform(TString bcalEnd="DS", int layer=1, int run=30894)
Definition: drawWaveforms.C:24
double par0err
Definition: drawWaveforms.C:12
double par2
Definition: drawWaveforms.C:11
TF1 * f2
Definition: FitGains.C:28
static TH1I * pedestal[nChan]
double par0
Definition: drawWaveforms.C:9
double par1err
Definition: drawWaveforms.C:12