130b3ca35f2a5d3271049db7d8b7871086240276
[u/mrichter/AliRoot.git] / T0 / AliT0CalibWalk.cxx
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
37
38 ClassImp(AliT0CalibWalk)
39
40 //________________________________________________________________
41   AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
42                                       fWalk(0),
43                                       fAmpLEDRec(0)
44 {
45   //
46 }
47
48 //________________________________________________________________
49 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
50                                       fWalk(0),
51                                       fAmpLEDRec(0)                                   
52 {
53   TString namst = "Calib_";
54   namst += name;
55   SetName(namst.Data());
56   SetTitle(namst.Data());
57
58 }
59
60 //________________________________________________________________
61 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
62   TNamed(calibda),              
63   fWalk(0),
64   fAmpLEDRec(0)
65
66 {
67 // copy constructor
68   SetName(calibda.GetName());
69   SetTitle(calibda.GetName());
70
71
72 }
73
74 //________________________________________________________________
75 AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
76 {
77 // assignment operator
78   SetName(calibda.GetName());
79   SetTitle(calibda.GetName());
80  
81   return *this;
82 }
83
84 //________________________________________________________________
85 AliT0CalibWalk::~AliT0CalibWalk()
86 {
87   //
88 }
89 //________________________________________________________________
90 void AliT0CalibWalk::SetWalk(Int_t ipmt)
91 {
92   //read QTC walk graph from external file
93
94   Int_t mv, ps; 
95   Int_t x[70000], y[70000], index[70000];
96   Float_t time[10000],amplitude[10000];
97   string buffer;
98   Bool_t down=false;
99   
100   const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-Amp.txt");
101   ifstream inFile(filename);
102   if(!inFile) AliError(Form("Cannot open file %s !",filename));
103   
104   Int_t i=0;
105   while(getline(inFile,buffer)){
106     inFile >> ps >> mv;
107
108     x[i]=ps; y[i]=mv;
109     i++;
110   }
111   inFile.close();
112   cout<<" number of data "<<i<<endl;
113  
114   TMath::Sort(i, y, index,down);
115   Int_t amp=0, iin=0, isum=0, sum=0;
116   Int_t ind=0;
117   for (Int_t ii=0; ii<i; ii++)
118     {
119       ind=index[ii];
120       if(y[ind] == amp)
121         {
122           sum +=x[ind];
123           iin++;
124         }
125       else
126         {
127           if(iin>0)
128             time[isum] = Float_t (sum/(iin));
129           else
130             time[isum] =Float_t (x[ind]);
131           amplitude[isum] = Float_t (amp);
132           amp=y[ind];
133           iin=0;
134           isum++;
135           sum=0;
136         }
137     }
138
139   inFile.close();
140
141   TGraph* gr = new TGraph(isum, amplitude, time);
142   fWalk.AddAtAndExpand(gr,ipmt);
143 }
144
145 //________________________________________________________________
146
147 void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
148 {
149   // read LED walk from external file
150   Float_t mv, ps; 
151   Float_t x[100], y[100];
152   string buffer;
153   
154   const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-LED.txt");
155    ifstream inFile(filename);
156   if(!inFile) {AliError(Form("Cannot open file %s !",filename));}
157   
158   inFile >> mv>>ps;
159   Int_t i=0;
160   
161   while(getline(inFile,buffer)){
162     x[i]=mv; y[i]=ps;   
163     inFile >> mv >> ps;
164     i++;
165   }
166   inFile.close();
167   Float_t y1[100], x1[100];
168   for (Int_t ir=0; ir<i; ir++){
169     y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
170   TGraph* gr = new TGraph(i,y1,x1);
171   fAmpLEDRec.AddAtAndExpand(gr,ipmt);
172   
173 }
174
175 //________________________________________________________________
176
177 void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
178 {
179   //make walk corerction for preprocessor
180
181   TFile *gFile = TFile::Open(laserFile);
182   //gSystem->Load("libSpectrum");
183  
184   Int_t npeaks = 20;
185   Int_t sigma=3.;
186   Bool_t down = false;
187
188   Int_t index[20];
189   Char_t buf1[10], buf2[10],buf3[20];
190   
191   
192   for (Int_t i=0; i<24; i++)
193     {
194       if(i>11)  {  sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
195
196       if(i<12)  { sprintf(buf1,"T0_C_%i_CFD",i+1); }
197
198       sprintf(buf2,"CFD_QTC%i",i+1);
199       sprintf(buf3,"CFD_LED%i",i+1);
200       cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
201       TH2F *qtc_cfd = (TH2F*) gFile->Get(buf2);
202       TH2F *led_cfd = (TH2F*) gFile->Get(buf3);
203       TH1F *cfd = (TH1F*) gFile->Get(buf1);
204       //       cfd->Draw();
205       TSpectrum *s = new TSpectrum(2*npeaks,1);
206       Int_t nfound = s->Search(cfd,sigma," ",0.05);
207       cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
208       if(nfound!=0){
209         Float_t *xpeak = s->GetPositionX();
210         TMath::Sort(nfound, xpeak, index,down);
211         Float_t xp = xpeak[index[0]];
212         Int_t xbin = cfd->GetXaxis()->FindBin(xp);
213         Float_t yp = cfd->GetBinContent(xbin);
214         cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
215         Float_t hmax = xp+10*sigma;
216         Float_t hmin = xp-10*sigma;
217         cout<<hmin<< " "<<hmax<<endl;
218         //          cfd->GetXaxis()->SetRange(hmin,hmax);
219         //          TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
220         // cfd->Fit("g1","R");
221         //      Int_t hminbin=  qtc_cfd->GetXaxis()->GetFirst();
222         //      Int_t hmaxbin=  qtc_cfd->GetXaxis()->GetLast();
223         Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
224         cout<<" nbins "<<nbins<<endl;
225         //      cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
226         //  qtc_cfd->Draw();
227         TProfile *pr_y = qtc_cfd->ProfileX();
228         //      Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
229         
230         pr_y->SetMaximum(hmax);
231         pr_y->SetMinimum(hmin);
232         Int_t np=nbins/20;
233         Double_t *xx = new Double_t[np];
234         Double_t *yy = new Double_t[np];
235         Int_t ng=0;
236         Double_t yg=0;
237         for (Int_t ip=1; ip<nbins; ip++)
238           {
239             if(ip%20 != 0 ) {
240               if (pr_y->GetBinContent(ip) !=0)
241                 yg +=pr_y->GetBinContent(ip);
242               //        cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
243               ng++;}
244             else {
245               xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
246               yy[ip/20] = yg/ng;
247               yg=0;
248               ng=0;
249               //                      cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
250             }
251           }
252         TGraph *gr = new TGraph(np,xx,yy);
253         gr->SetMinimum(hmin);
254         gr->SetMaximum(hmax);
255         fWalk.AddAtAndExpand(gr,i);       
256         //      fWalk.AddAtAndExpand(gr,i+12);
257
258         Int_t nbinsled= led_cfd->GetXaxis()->GetNbins();
259         cout<<" nbins led "<<nbinsled<<endl;
260         //      cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
261         //  qtc_cfd->Draw();
262         TProfile *prled_y = led_cfd->ProfileX();
263         //      Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
264         
265         prled_y->SetMaximum(hmax);
266         prled_y->SetMinimum(hmin);
267         Int_t npled=nbinsled/20;
268         Double_t *xxled = new Double_t[np];
269         Double_t *yyled = new Double_t[np];
270         Int_t ngled=0;
271         Double_t ygled=0;
272         for (Int_t ip=1; ip<nbinsled; ip++)
273           {
274             if(ip%20 != 0 ) {
275               if (prled_y->GetBinContent(ip) !=0)
276                 ygled +=prled_y->GetBinContent(ip);
277               //        cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
278               ngled++;}
279             else {
280               xxled[ip/20] = Float_t (prled_y->GetBinCenter(ip));
281               yyled[ip/20] = ygled/ngled;
282               ygled=0;
283               ngled=0;
284             }
285           }
286         TGraph *grled = new TGraph(npled,xxled,yyled);
287         grled->SetMinimum(hmin);
288         grled->SetMaximum(hmax);
289         fAmpLEDRec.AddAtAndExpand(grled,i);       
290         //      delete [] xx;
291         //      delete [] yy;
292         
293       }
294       
295     }
296   
297
298
299 }
300