working version of preprocessor
[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  
183   Int_t npeaks = 20;
184   Double_t sigma=3.;
185   Bool_t down = false;
186
187   Int_t index[20];
188   Char_t buf1[10], buf2[10],buf3[20];
189   
190   
191   for (Int_t i=0; i<24; i++)
192     {
193       if(i>11)  {  sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
194
195       if(i<12)  { sprintf(buf1,"T0_C_%i_CFD",i+1); }
196
197       sprintf(buf2,"CFD_QTC%i",i+1);
198       sprintf(buf3,"CFD_LED%i",i+1);
199       cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
200       TH2F *qtcVScfd = (TH2F*) gFile->Get(buf2);
201       TH2F *ledVScfd = (TH2F*) gFile->Get(buf3);
202       TH1F *cfd = (TH1F*) gFile->Get(buf1);
203       TSpectrum *s = new TSpectrum(2*npeaks,1);
204       Int_t nfound = s->Search(cfd,sigma,"goff",0.05);
205       cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
206       if(nfound!=0){
207         Float_t *xpeak = s->GetPositionX();
208         TMath::Sort(nfound, xpeak, index,down);
209         Float_t xp = xpeak[index[0]];
210         Int_t xbin = cfd->GetXaxis()->FindBin(xp);
211         Float_t yp = cfd->GetBinContent(xbin);
212         cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
213         Float_t hmax = xp+10*sigma;
214         Float_t hmin = xp-10*sigma;
215         cout<<hmin<< " "<<hmax<<endl;
216
217         //QTC
218         Int_t nbins= qtcVScfd->GetXaxis()->GetNbins();
219         cout<<" nbins "<<nbins<<endl;
220         TProfile *prY = qtcVScfd->ProfileX();
221         
222         prY->SetMaximum(hmax);
223         prY->SetMinimum(hmin);
224         Int_t np=nbins/20;
225         Double_t *xx = new Double_t[np];
226         Double_t *yy = new Double_t[np];
227         Int_t ng=0;
228         Double_t yg=0;
229         for (Int_t ip=1; ip<nbins; ip++)
230           {
231             if(ip%20 != 0 ) {
232               if (prY->GetBinContent(ip) !=0)
233                 yg +=prY->GetBinContent(ip);
234               ng++;}
235             else {
236               xx[ip/20] = Float_t (prY->GetBinCenter(ip));
237               yy[ip/20] = yg/ng;
238               yg=0;
239               ng=0;
240             }
241           }
242         TGraph *gr = new TGraph(np,xx,yy);
243         gr->SetMinimum(hmin);
244         gr->SetMaximum(hmax);
245         fWalk.AddAtAndExpand(gr,i);       
246
247         //LED
248         Int_t nbinsled= ledVScfd->GetXaxis()->GetNbins();
249         cout<<" nbins led "<<nbinsled<<endl;
250         TProfile *prledY = ledVScfd->ProfileX();
251         
252         prledY->SetMaximum(hmax);
253         prledY->SetMinimum(hmin);
254         Int_t npled=nbinsled/20;
255         Double_t *xxled = new Double_t[np];
256         Double_t *yyled = new Double_t[np];
257         Int_t ngled=0;
258         Double_t ygled=0;
259         for (Int_t ip=1; ip<nbinsled; ip++)
260           {
261             if(ip%20 != 0 ) {
262               if (prledY->GetBinContent(ip) !=0)
263                 ygled +=prledY->GetBinContent(ip);
264               ngled++;}
265             else {
266               xxled[ip/20] = Float_t (prledY->GetBinCenter(ip));
267               yyled[ip/20] = ygled/ngled;
268               ygled=0;
269               ngled=0;
270             }
271           }
272         TGraph *grled = new TGraph(npled,xxled,yyled);
273         grled->SetMinimum(hmin);
274         grled->SetMaximum(hmax);
275         fAmpLEDRec.AddAtAndExpand(grled,i);       
276
277         delete [] xx;
278         delete [] yy;
279         delete [] xxled;
280         delete [] yyled;
281         delete prY;
282         delete prledY;
283         
284       }
285       delete cfd;
286       delete qtcVScfd;
287       delete ledVScfd;
288     }
289   
290 }
291