Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /* $Id$ */
17 :
18 : ///////////////////////////////////////////////////////////////////////////////
19 : // //
20 : // class T0 walk correction Alla Maevskaya alla@inr.ru 20.11.2007 //
21 : // //
22 : ///////////////////////////////////////////////////////////////////////////////
23 :
24 : #include "AliT0CalibWalk.h"
25 : #include "AliLog.h"
26 :
27 : #include <TObjArray.h>
28 : #include <TGraph.h>
29 : #include <TFile.h>
30 : #include <TH2F.h>
31 : #include <TMath.h>
32 : #include <TSystem.h>
33 : #include <Riostream.h>
34 : #include <TSpectrum.h>
35 : #include <TProfile.h>
36 : #include <TF1.h>
37 :
38 : using std::cout;
39 : using std::endl;
40 : using std::ifstream;
41 20 : ClassImp(AliT0CalibWalk)
42 :
43 : //________________________________________________________________
44 3 : AliT0CalibWalk::AliT0CalibWalk(): TNamed(),
45 3 : fWalk(0),
46 3 : fAmpLEDRec(0),
47 3 : fQTC(0),
48 3 : fAmpLED(0),
49 3 : fCalibByData(kFALSE)
50 15 : {
51 : //
52 6 : }
53 :
54 : //________________________________________________________________
55 0 : AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
56 0 : fWalk(0),
57 0 : fAmpLEDRec(0), fQTC(0),
58 0 : fAmpLED(0),
59 0 : fCalibByData(kFALSE)
60 :
61 0 : {
62 0 : TString namst = "Calib_";
63 0 : namst += name;
64 0 : SetName(namst.Data());
65 0 : SetTitle(namst.Data());
66 :
67 0 : }
68 :
69 : //________________________________________________________________
70 : AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
71 0 : TNamed(calibda),
72 0 : fWalk(0),
73 0 : fAmpLEDRec(0),
74 0 : fQTC(0),
75 0 : fAmpLED(0) ,
76 0 : fCalibByData(kFALSE)
77 0 : {
78 : // copy constructor
79 0 : SetName(calibda.GetName());
80 0 : SetTitle(calibda.GetName());
81 :
82 :
83 0 : }
84 :
85 : //________________________________________________________________
86 : AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
87 : {
88 : // assignment operator
89 0 : SetName(calibda.GetName());
90 0 : SetTitle(calibda.GetName());
91 :
92 0 : return *this;
93 : }
94 :
95 : //________________________________________________________________
96 : AliT0CalibWalk::~AliT0CalibWalk()
97 0 : {
98 : //
99 0 : }
100 : //________________________________________________________________
101 :
102 : Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
103 : {
104 : //make walk corerction for preprocessor
105 0 : Float_t sigma,cfdmean, qtmean, ledmean;
106 : Bool_t ok=true;
107 0 : Float_t mips[50];
108 0 : cout<<" @@@@ fCalibByData "<<fCalibByData<<endl;
109 :
110 0 : gFile = TFile::Open(laserFile);
111 0 : if(!gFile) {
112 0 : AliError("No input laser data found ");
113 0 : }
114 : else
115 : {
116 : // gFile->ls();
117 0 : TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
118 : Int_t nmips=0;
119 0 : for (Int_t ibin=0; ibin<2000; ibin++) {
120 0 : Float_t bincont = hAmp->GetBinContent(ibin);
121 0 : if(bincont>0){
122 0 : mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
123 0 : cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
124 0 : nmips++;
125 0 : }
126 : }
127 : /*
128 : if (nmips<17) {
129 : ok=false;
130 : return ok;
131 : }
132 : */
133 0 : Float_t x1[50], y1[50];
134 0 : Float_t x2[50], xx2[50],y2[50];
135 0 : Float_t xx1[50],yy1[50], xx[50];
136 :
137 : Float_t cfd0 = 0;
138 :
139 0 : for (Int_t ii=0; ii<nmips; ii++)
140 0 : x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0;
141 :
142 0 : for (Int_t i=0; i<24; i++)
143 : {
144 : cfd0 = 0;
145 0 : for (Int_t im=0; im<nmips; im++)
146 : {
147 0 : TString cfd = Form("hCFD%i_%i",i+1,im+1);
148 0 : TString qtc = Form("hQTC%i_%i",i+1,im+1);
149 0 : TString led = Form("hLED%i_%i",i+1,im+1);
150 :
151 0 : TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
152 0 : TH1F *hLED = (TH1F*) gFile->Get(led.Data());
153 0 : TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
154 : // hCFD->SetDirectory(0);
155 : // hQTC->SetDirectory(0);
156 : // hLED->SetDirectory(0);
157 0 : if(!hCFD )
158 0 : AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
159 0 : if(!hQTC )
160 0 : AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
161 0 : if(!hLED)
162 0 : AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
163 0 : if( hCFD && hCFD->GetEntries()<500 ) {
164 : ok=false;
165 0 : printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
166 0 : return ok;
167 : }
168 0 : if(hCFD && hCFD->GetEntries()>500 ) {
169 0 : if( hCFD->GetRMS() >= 1.5)
170 0 : GetMeanAndSigma(hCFD, cfdmean, sigma);
171 : else
172 0 : cfdmean = hCFD->GetMean();
173 :
174 0 : Int_t maxBin = hCFD->GetMaximumBin();
175 0 : Double_t meanEstimate = hCFD->GetBinCenter( maxBin);
176 0 : if(TMath::Abs(meanEstimate - cfdmean) > 20 ) cfdmean = meanEstimate;
177 0 : if (im == 0) cfd0 = cfdmean;
178 0 : y1[im] = cfdmean - cfd0;
179 0 : }
180 0 : if(hQTC && hQTC->GetEntries()>500) {
181 0 : GetMeanAndSigma(hQTC, qtmean, sigma);
182 :
183 0 : x1[im] = qtmean;
184 0 : if( x1[im] == 0) {
185 : ok=false;
186 0 : printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
187 0 : return ok;
188 : }
189 : }
190 :
191 0 : if( hLED && hLED->GetEntries()>500) {
192 0 : GetMeanAndSigma(hLED, ledmean, sigma);
193 : }
194 : else
195 : {
196 : // ok=false;
197 0 : printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
198 0 : ledmean=0;
199 : // return ok;
200 : }
201 0 : x2[im] = ledmean;
202 0 : xx2[im] = x2[nmips-im-1];
203 :
204 0 : xx[im]=mips[im];
205 0 : if (hQTC) delete hQTC;
206 0 : if (hCFD) delete hCFD;
207 0 : if (hLED) delete hLED;
208 :
209 0 : }
210 :
211 0 : for (Int_t imi=0; imi<nmips; imi++)
212 : {
213 0 : yy1[imi] = Float_t (mips[nmips-imi-1]);
214 0 : xx1[imi] = x2[nmips-imi-1];
215 : // cout<<" LED rev "<<i<<" "<<imi<<" "<<xx1[imi]<<" "<< yy1[imi]<<"nmips-imi-1 = "<< nmips-imi-1<<endl;
216 : }
217 :
218 0 : if(i==0) cout<<"Making graphs..."<<endl;
219 :
220 0 : TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
221 0 : grwalkqtc->SetTitle(Form("PMT%i",i));
222 0 : TGraph *grwalkled = new TGraph (nmips,x2,y1);
223 0 : grwalkled->SetTitle(Form("PMT%i",i));
224 0 : if(!fCalibByData)
225 0 : fWalk.AddAtAndExpand(grwalkqtc,i);
226 0 : fAmpLEDRec.AddAtAndExpand(grwalkled,i);
227 : // cout<<" add walk "<<i<<endl;
228 :
229 : //fit amplitude graphs to make comparison wth new one
230 0 : TGraph *grampled = new TGraph (nmips,xx1,yy1);
231 0 : TGraph *grqtc = new TGraph (nmips,x1,xx);
232 0 : fQTC.AddAtAndExpand(grqtc,i);
233 0 : fAmpLED.AddAtAndExpand(grampled,i);
234 : // cout<<" add amp "<<i<<endl;
235 :
236 0 : if(i==23)
237 0 : cout<<"Graphs created..."<<endl;
238 : }
239 0 : }
240 0 : Float_t xpoint, ypoint, xdata[250], ydata[250];
241 0 : Int_t ipmt;
242 0 : if(fCalibByData) {
243 0 : cout<<" read ingraph "<<endl;
244 0 : ifstream ingraph ("calibfit.txt");
245 0 : for (Int_t i=0; i<24; i++)
246 : {
247 0 : for (Int_t ip=0; ip<200; ip++)
248 : {
249 0 : ingraph>>ipmt>>xpoint>>ypoint;
250 : // cout<<i<<" "<<ipmt<<" "<<ip<<" "<< xpoint<<" "<<ypoint<<endl;
251 0 : xdata[ip]=xpoint;
252 0 : ydata[ip]=ypoint;
253 : }
254 0 : for (Int_t ip=200; ip<250; ip++) {
255 0 : xdata[ip] =xdata[ip-1]+10;
256 0 : ydata[ip]=ydata[199];
257 : }
258 :
259 0 : TGraph *grwalkqtc = new TGraph (250,xdata,ydata);
260 0 : grwalkqtc->Print();
261 0 : grwalkqtc->SetTitle(Form("PMT%i",i) );
262 0 : fWalk.AddAtAndExpand(grwalkqtc,i);
263 0 : for (Int_t ip=0; ip<250; ip++)
264 : {
265 0 : xdata[ip]=0;
266 0 : ydata[ip]=0;
267 : }
268 : }
269 0 : }
270 :
271 :
272 0 : return ok;
273 0 : }
274 : void AliT0CalibWalk::SetWalk2015(TString filename)
275 : {
276 0 : printf("AliT0CalibWalk::SetWalk2015\n");
277 : //set zero LED correction
278 0 : Float_t ampled[200], walkled[200], ampqtc[200];;
279 0 : for (int i=0; i<200; i++) {
280 0 : ampled[i] = Float_t(i*10);
281 0 : walkled[i]=0;
282 0 : ampqtc[i] = Float_t (i*100);
283 : }
284 : Int_t nbins=0;
285 0 : Double_t xgr[350], ygr[350];
286 0 : Float_t ampcut10[24] = {1563, 1544, 1555, 1645, 1745, 1555,
287 : 1526, 1455, 1739, 1413, 1614, 1194,
288 : 1693, 1454, 1566, 1670, 1516, 1810,
289 : 1607, 1516, 1675, 1566, 1443, 1545 };
290 :
291 0 : Float_t meanmax[24]= {3077.13, 3069.43, 3070.27, 3123.97, 3150.55, 3127.22,
292 : 3110.7, 3088.24, 3164.09, 3154.94, 3142.93, 3103.78,
293 : 3091.14, 3055.07, 3043.82, 3051.3, 3062.58, 3087.11,
294 : 3055.03, 3072.61, 3064.25, 3120.34, 3115.8, 3095.93};
295 0 : TProfile * prQTC_CFD[24];
296 0 : TFile *f = new TFile(filename.Data());
297 0 : if (!f) {
298 0 : printf (" no file \n");
299 0 : return;
300 : }
301 0 : for (int i=0; i<24; i++) {
302 0 : prQTC_CFD[i]= (TProfile*) f->Get(Form("Slew%i",i+1) );
303 : }
304 : // collect graph
305 0 : for (Int_t i=0; i<24; i++)
306 : {
307 : nbins=0;
308 0 : cout<<" nachalo "<<i<<endl;
309 0 : for (Int_t j=10; j<350; j++) {
310 0 : Int_t nentr = prQTC_CFD[i]->GetBinEntries(j);
311 0 : Float_t prqt = prQTC_CFD[i]->GetBinContent(j);
312 0 : xgr[nbins]=prQTC_CFD[i]->GetBinCenter(j);
313 0 : cout<<" bin "<<j<<" cont "<<prqt<<endl;
314 0 : if (prQTC_CFD[i]->GetBinCenter(j)>ampcut10[i]) {
315 0 : if(nentr>1000) {
316 0 : ygr[nbins]=prqt-meanmax[i];
317 0 : }
318 : else
319 0 : ygr[nbins]= ygr[nbins-1];
320 0 : cout<<nbins<<" "<< xgr[nbins]<<" "<< ygr[nbins]<<endl;
321 0 : nbins++;
322 0 : }
323 : }
324 0 : TGraph *grwalkqtc = new TGraph (nbins,xgr,ygr);
325 0 : grwalkqtc->SetTitle(Form("PMT%i",i+1));
326 0 : fWalk.AddAtAndExpand(grwalkqtc,i);
327 0 : TGraph *grwalkled = new TGraph (100,ampled,walkled);
328 0 : fAmpLEDRec.AddAtAndExpand(grwalkled,i);
329 0 : grwalkled->SetTitle(Form("PMT%i",i+1));
330 :
331 : //fit amplitude graphs to make comparison wth new one
332 0 : TGraph *grampled = new TGraph (200,ampled,ampled);
333 0 : TGraph *grqtc = new TGraph (200,ampqtc,ampqtc);
334 0 : fQTC.AddAtAndExpand(grqtc,i);
335 0 : fAmpLED.AddAtAndExpand(grampled,i);
336 0 : cout<<" add amp "<<i<<endl;
337 :
338 :
339 : }
340 0 : printf(" AliT0CalibWalk return \n");
341 :
342 0 : }
343 : //__________________________________________________________________________________________
344 : void AliT0CalibWalk::SetWalkDima(TString filename)
345 : {
346 0 : Float_t ampled[200], walkled[200], ampqtc[200];;
347 0 : for (int i=0; i<200; i++) {
348 0 : ampled[i] = Float_t(i*10);
349 0 : walkled[i]=0;
350 0 : ampqtc[i] = Float_t (i*100);
351 : }
352 0 : TFile *file = new TFile(filename.Data());
353 0 : file->ls();
354 0 : if (!file) {
355 0 : printf (" no file \n");
356 0 : return;
357 : }
358 0 : TDirectoryFile *dr = (TDirectoryFile*)file->Get("resultGraphs");
359 : TGraph *currGraph ;
360 0 : TString aPMTname;
361 0 : for(Int_t iPMT = 0; iPMT < 24; iPMT++)
362 : {
363 0 : if(iPMT<12) aPMTname = Form("C_%02d_QTCCFDgraph", iPMT+1);
364 0 : else aPMTname = Form("A_%02d_QTCCFDgraph", iPMT-12+1);
365 0 : printf(" %s \n",aPMTname.Data());
366 0 : currGraph = (TGraph*)dr->Get(aPMTname.Data());
367 : // currGraph = (TGraph*)file->FindObjectAny(aPMTname.Data());
368 0 : currGraph->SetTitle(Form("PMT%i",iPMT+1));
369 0 : fWalk.AddAtAndExpand(currGraph,iPMT);
370 : // fWalk.At(iPMT)->Print();
371 0 : TGraph *grwalkled = new TGraph (100,ampled,walkled);
372 0 : fAmpLEDRec.AddAtAndExpand(grwalkled,iPMT);
373 :
374 0 : TGraph *grampled = new TGraph (200,ampled,ampled);
375 0 : TGraph *grqtc = new TGraph (200,ampqtc,ampqtc);
376 0 : fQTC.AddAtAndExpand(grqtc,iPMT);
377 0 : fAmpLED.AddAtAndExpand(grampled,iPMT);
378 0 : cout<<" add amp "<<iPMT<<endl;
379 : }
380 :
381 :
382 0 : }
383 : //__________________________________________________________________________________________
384 : void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma)
385 : {
386 :
387 : const double window = 2.; //fit window
388 :
389 : double meanEstimate, sigmaEstimate;
390 : int maxBin;
391 0 : maxBin = hist->GetMaximumBin(); //position of maximum
392 0 : meanEstimate = hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
393 0 : sigmaEstimate = hist->GetRMS();
394 0 : TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
395 0 : fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
396 0 : hist->Fit("fit","RQ","Q");
397 :
398 0 : mean = (Float_t) fit->GetParameter(1);
399 0 : sigma = (Float_t) fit->GetParameter(2);
400 : // printf(" mean %f max %f sigma %f \n",mean, meanEstimate , sigma);
401 :
402 0 : delete fit;
403 0 : }
404 :
405 : //______________________________________________________________________________
406 : void AliT0CalibWalk::SetWalkZero()
407 : {
408 0 : Float_t amp[100], walk[100];
409 0 : for (int i=0; i<100; i++) {
410 0 : amp[i] = Float_t(i*100);
411 0 : walk[i]=0;
412 : }
413 : TGraph *gramp, *grwalk;
414 0 : for (int igr=0; igr<24; igr++) {
415 0 : gramp = new TGraph(100, amp,amp);
416 0 : grwalk = new TGraph(100, amp,walk);
417 0 : grwalk->SetTitle(Form("PMT%i",igr+1));
418 0 : gramp->SetTitle(Form("PMT%i",igr+1));
419 0 : fWalk.AddAtAndExpand(grwalk,igr);
420 0 : fAmpLEDRec.AddAtAndExpand(grwalk,igr);
421 0 : fQTC.AddAtAndExpand(gramp,igr);
422 0 : fAmpLED.AddAtAndExpand(gramp,igr);
423 :
424 : }
425 0 : }
|
