]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - T0/AliT0CalibWalk.cxx
ALIROOT-5507 Adding root files to the OADB.par
[u/mrichter/AliRoot.git] / T0 / AliT0CalibWalk.cxx
... / ...
CommitLineData
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
38using std::cout;
39using std::endl;
40using std::ifstream;
41ClassImp(AliT0CalibWalk)
42
43//________________________________________________________________
44 AliT0CalibWalk::AliT0CalibWalk(): TNamed(),
45 fWalk(0),
46 fAmpLEDRec(0),
47 fQTC(0),
48 fAmpLED(0),
49 fCalibByData(kFALSE)
50{
51 //
52}
53
54//________________________________________________________________
55AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
56 fWalk(0),
57 fAmpLEDRec(0), fQTC(0),
58 fAmpLED(0),
59 fCalibByData(kFALSE)
60
61{
62 TString namst = "Calib_";
63 namst += name;
64 SetName(namst.Data());
65 SetTitle(namst.Data());
66
67}
68
69//________________________________________________________________
70AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
71 TNamed(calibda),
72 fWalk(0),
73 fAmpLEDRec(0),
74 fQTC(0),
75 fAmpLED(0) ,
76 fCalibByData(kFALSE)
77{
78// copy constructor
79 SetName(calibda.GetName());
80 SetTitle(calibda.GetName());
81
82
83}
84
85//________________________________________________________________
86AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
87{
88// assignment operator
89 SetName(calibda.GetName());
90 SetTitle(calibda.GetName());
91
92 return *this;
93}
94
95//________________________________________________________________
96AliT0CalibWalk::~AliT0CalibWalk()
97{
98 //
99}
100//________________________________________________________________
101
102Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
103{
104 //make walk corerction for preprocessor
105 Float_t sigma,cfdmean, qtmean, ledmean;
106 Bool_t ok=true;
107 Float_t mips[50];
108 cout<<" @@@@ fCalibByData "<<fCalibByData<<endl;
109
110 gFile = TFile::Open(laserFile);
111 if(!gFile) {
112 AliError("No input laser data found ");
113 }
114 else
115 {
116 // gFile->ls();
117 TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
118 Int_t nmips=0;
119 for (Int_t ibin=0; ibin<2000; ibin++) {
120 Float_t bincont = hAmp->GetBinContent(ibin);
121 if(bincont>0){
122 mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
123 cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
124 nmips++;
125 }
126 }
127 /*
128 if (nmips<17) {
129 ok=false;
130 return ok;
131 }
132 */
133 Float_t x1[50], y1[50];
134 Float_t x2[50], xx2[50],y2[50];
135 Float_t xx1[50],yy1[50], xx[50];
136
137 Float_t cfd0 = 0;
138
139 for (Int_t ii=0; ii<nmips; ii++)
140 x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0;
141
142 for (Int_t i=0; i<24; i++)
143 {
144 cfd0 = 0;
145 for (Int_t im=0; im<nmips; im++)
146 {
147 TString cfd = Form("hCFD%i_%i",i+1,im+1);
148 TString qtc = Form("hQTC%i_%i",i+1,im+1);
149 TString led = Form("hLED%i_%i",i+1,im+1);
150
151 TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
152 TH1F *hLED = (TH1F*) gFile->Get(led.Data());
153 TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
154 // hCFD->SetDirectory(0);
155 // hQTC->SetDirectory(0);
156 // hLED->SetDirectory(0);
157 if(!hCFD )
158 AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
159 if(!hQTC )
160 AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
161 if(!hLED)
162 AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
163 if( hCFD && hCFD->GetEntries()<500 ) {
164 ok=false;
165 printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
166 return ok;
167 }
168 if(hCFD && hCFD->GetEntries()>500 ) {
169 if( hCFD->GetRMS() >= 1.5)
170 GetMeanAndSigma(hCFD, cfdmean, sigma);
171 else
172 cfdmean = hCFD->GetMean();
173
174 Int_t maxBin = hCFD->GetMaximumBin();
175 Double_t meanEstimate = hCFD->GetBinCenter( maxBin);
176 if(TMath::Abs(meanEstimate - cfdmean) > 20 ) cfdmean = meanEstimate;
177 if (im == 0) cfd0 = cfdmean;
178 y1[im] = cfdmean - cfd0;
179 }
180 if(hQTC && hQTC->GetEntries()>500) {
181 GetMeanAndSigma(hQTC, qtmean, sigma);
182
183 x1[im] = qtmean;
184 if( x1[im] == 0) {
185 ok=false;
186 printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
187 return ok;
188 }
189 }
190
191 if( hLED && hLED->GetEntries()>500) {
192 GetMeanAndSigma(hLED, ledmean, sigma);
193 }
194 else
195 {
196 // ok=false;
197 printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
198 ledmean=0;
199 // return ok;
200 }
201 x2[im] = ledmean;
202 xx2[im] = x2[nmips-im-1];
203
204 xx[im]=mips[im];
205 if (hQTC) delete hQTC;
206 if (hCFD) delete hCFD;
207 if (hLED) delete hLED;
208
209 }
210
211 for (Int_t imi=0; imi<nmips; imi++)
212 {
213 yy1[imi] = Float_t (mips[nmips-imi-1]);
214 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 if(i==0) cout<<"Making graphs..."<<endl;
219
220 TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
221 grwalkqtc->SetTitle(Form("PMT%i",i));
222 TGraph *grwalkled = new TGraph (nmips,x2,y1);
223 grwalkled->SetTitle(Form("PMT%i",i));
224 if(!fCalibByData)
225 fWalk.AddAtAndExpand(grwalkqtc,i);
226 fAmpLEDRec.AddAtAndExpand(grwalkled,i);
227 // cout<<" add walk "<<i<<endl;
228
229 //fit amplitude graphs to make comparison wth new one
230 TGraph *grampled = new TGraph (nmips,xx1,yy1);
231 TGraph *grqtc = new TGraph (nmips,x1,xx);
232 fQTC.AddAtAndExpand(grqtc,i);
233 fAmpLED.AddAtAndExpand(grampled,i);
234 // cout<<" add amp "<<i<<endl;
235
236 if(i==23)
237 cout<<"Graphs created..."<<endl;
238 }
239 }
240 Float_t xpoint, ypoint, xdata[250], ydata[250];
241 Int_t ipmt;
242 if(fCalibByData) {
243 cout<<" read ingraph "<<endl;
244 ifstream ingraph ("calibfit.txt");
245 for (Int_t i=0; i<24; i++)
246 {
247 for (Int_t ip=0; ip<200; ip++)
248 {
249 ingraph>>ipmt>>xpoint>>ypoint;
250 // cout<<i<<" "<<ipmt<<" "<<ip<<" "<< xpoint<<" "<<ypoint<<endl;
251 xdata[ip]=xpoint;
252 ydata[ip]=ypoint;
253 }
254 for (Int_t ip=200; ip<250; ip++) {
255 xdata[ip] =xdata[ip-1]+10;
256 ydata[ip]=ydata[199];
257 }
258
259 TGraph *grwalkqtc = new TGraph (250,xdata,ydata);
260 grwalkqtc->Print();
261 grwalkqtc->SetTitle(Form("PMT%i",i));
262 fWalk.AddAtAndExpand(grwalkqtc,i);
263 for (Int_t ip=0; ip<250; ip++)
264 {
265 xdata[ip]=0;
266 ydata[ip]=0;
267 }
268 }
269 }
270
271
272return ok;
273}
274
275void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma)
276{
277
278 const double window = 2.; //fit window
279
280 double meanEstimate, sigmaEstimate;
281 int maxBin;
282 maxBin = hist->GetMaximumBin(); //position of maximum
283 meanEstimate = hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
284 sigmaEstimate = hist->GetRMS();
285 TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
286 fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
287 hist->Fit("fit","RQ","Q");
288
289 mean = (Float_t) fit->GetParameter(1);
290 sigma = (Float_t) fit->GetParameter(2);
291 // printf(" mean %f max %f sigma %f \n",mean, meanEstimate , sigma);
292
293 delete fit;
294}
295
296
297