1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // class T0 walk correction Alla Maevskaya alla@inr.ru 20.11.2007 //
22 ///////////////////////////////////////////////////////////////////////////////
24 #include "AliT0CalibWalk.h"
27 #include <TObjArray.h>
33 #include <Riostream.h>
34 #include <TSpectrum.h>
41 ClassImp(AliT0CalibWalk)
43 //________________________________________________________________
44 AliT0CalibWalk::AliT0CalibWalk(): TNamed(),
54 //________________________________________________________________
55 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
57 fAmpLEDRec(0), fQTC(0),
62 TString namst = "Calib_";
64 SetName(namst.Data());
65 SetTitle(namst.Data());
69 //________________________________________________________________
70 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
79 SetName(calibda.GetName());
80 SetTitle(calibda.GetName());
85 //________________________________________________________________
86 AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
88 // assignment operator
89 SetName(calibda.GetName());
90 SetTitle(calibda.GetName());
95 //________________________________________________________________
96 AliT0CalibWalk::~AliT0CalibWalk()
100 //________________________________________________________________
102 Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
104 //make walk corerction for preprocessor
105 Float_t sigma,cfdmean, qtmean, ledmean;
108 cout<<" @@@@ fCalibByData "<<fCalibByData<<endl;
110 gFile = TFile::Open(laserFile);
112 AliError("No input laser data found ");
117 TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
119 for (Int_t ibin=0; ibin<2000; ibin++) {
120 Float_t bincont = hAmp->GetBinContent(ibin);
122 mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
123 cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
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];
139 for (Int_t ii=0; ii<nmips; ii++)
140 x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0;
142 for (Int_t i=0; i<24; i++)
145 for (Int_t im=0; im<nmips; im++)
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);
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);
158 AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
160 AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
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 ) {
165 printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
168 if(hCFD && hCFD->GetEntries()>500 ) {
169 if( hCFD->GetRMS() >= 1.5)
170 GetMeanAndSigma(hCFD, cfdmean, sigma);
172 cfdmean = hCFD->GetMean();
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;
180 if(hQTC && hQTC->GetEntries()>500) {
181 GetMeanAndSigma(hQTC, qtmean, sigma);
186 printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
191 if( hLED && hLED->GetEntries()>500) {
192 GetMeanAndSigma(hLED, ledmean, sigma);
197 printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
202 xx2[im] = x2[nmips-im-1];
205 if (hQTC) delete hQTC;
206 if (hCFD) delete hCFD;
207 if (hLED) delete hLED;
211 for (Int_t imi=0; imi<nmips; imi++)
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;
218 if(i==0) cout<<"Making graphs..."<<endl;
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));
225 fWalk.AddAtAndExpand(grwalkqtc,i);
226 fAmpLEDRec.AddAtAndExpand(grwalkled,i);
227 // cout<<" add walk "<<i<<endl;
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;
237 cout<<"Graphs created..."<<endl;
240 Float_t xpoint, ypoint, xdata[250], ydata[250];
243 cout<<" read ingraph "<<endl;
244 ifstream ingraph ("calibfit.txt");
245 for (Int_t i=0; i<24; i++)
247 for (Int_t ip=0; ip<200; ip++)
249 ingraph>>ipmt>>xpoint>>ypoint;
250 // cout<<i<<" "<<ipmt<<" "<<ip<<" "<< xpoint<<" "<<ypoint<<endl;
254 for (Int_t ip=200; ip<250; ip++) {
255 xdata[ip] =xdata[ip-1]+10;
256 ydata[ip]=ydata[199];
259 TGraph *grwalkqtc = new TGraph (250,xdata,ydata);
261 grwalkqtc->SetTitle(Form("PMT%i",i));
262 fWalk.AddAtAndExpand(grwalkqtc,i);
263 for (Int_t ip=0; ip<250; ip++)
275 void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma)
278 const double window = 2.; //fit window
280 double meanEstimate, sigmaEstimate;
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");
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);