Bug fix in the order of the Ds cuts (Sadhana, Francesco)
[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 #include <TF1.h>
37
38 ClassImp(AliT0CalibWalk)
39
40 //________________________________________________________________
41   AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
42                                       fWalk(0),
43                                       fAmpLEDRec(0),
44                                       fQTC(0),
45                                       fAmpLED(0)
46 {
47   //
48 }
49
50 //________________________________________________________________
51 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
52                                       fWalk(0),
53                                       fAmpLEDRec(0),                                                                  fQTC(0),
54                                       fAmpLED(0)
55     
56 {
57   TString namst = "Calib_";
58   namst += name;
59   SetName(namst.Data());
60   SetTitle(namst.Data());
61
62 }
63
64 //________________________________________________________________
65 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
66   TNamed(calibda),              
67   fWalk(0),
68   fAmpLEDRec(0),
69   fQTC(0),
70   fAmpLED(0)
71  
72 {
73 // copy constructor
74   SetName(calibda.GetName());
75   SetTitle(calibda.GetName());
76
77
78 }
79
80 //________________________________________________________________
81 AliT0CalibWalk &AliT0CalibWalk::operator =(const AliT0CalibWalk& calibda)
82 {
83 // assignment operator
84   SetName(calibda.GetName());
85   SetTitle(calibda.GetName());
86  
87   return *this;
88 }
89
90 //________________________________________________________________
91 AliT0CalibWalk::~AliT0CalibWalk()
92 {
93   //
94 }
95 //________________________________________________________________  
96
97 Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
98 {
99   //make walk corerction for preprocessor
100   Int_t npeaks = 20;
101   Int_t sigma=3.;
102   Bool_t down=false;
103   Int_t index[20];
104   Bool_t ok=true;
105   Float_t   mips[20];
106   Int_t nfound=0;
107   
108   gFile = TFile::Open(laserFile);
109   if(!gFile) {
110     AliError("No input laser data found ");
111   }
112   else
113     {
114       //      gFile->ls();
115       TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
116       Int_t nmips=0;
117       for (Int_t ibin=0; ibin<1000; ibin++) {
118         Float_t bincont = hAmp->GetBinContent(ibin);
119         if(bincont>0){ 
120           mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
121           cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
122           nmips++;
123         }       
124       }    
125       
126       if (nmips<15) ok=false; 
127        
128       Float_t x1[50], y1[50]; 
129       Float_t x2[50], xx2[50],y2[50];
130       Float_t xx1[50],yy1[50], xx[50];
131       
132       Float_t cfd0[24];
133       Int_t startim = 0;
134       
135       for (Int_t ii=0; ii<nmips; ii++)
136         x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0; 
137       
138       for (Int_t i=0; i<24; i++)
139         {
140           for (Int_t im=startim; im<nmips; im++)
141             {         
142               TString cfd = Form("hCFD%i_%i",i+1,im+1);
143               TString qtc = Form("hQTC%i_%i",i+1,im+1);
144               TString led = Form("hLED%i_%i",i+1,im+1);
145               
146               TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
147               TH1F *hLED = (TH1F*) gFile->Get(led.Data());
148               TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
149               hCFD->SetDirectory(0);
150               hQTC->SetDirectory(0);
151               hLED->SetDirectory(0);
152               if(!hCFD )
153                 AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
154               if(!hQTC )
155                 AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
156               if(!hLED)       
157                 AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
158               
159               if(hCFD ) {
160                 hCFD->GetXaxis()->SetRangeUser(32200,32450);
161                 TSpectrum *s = new TSpectrum(2*npeaks,1);
162                 nfound = s->Search(hCFD,sigma," ",0.1);
163                 if(nfound!=0){
164                   Float_t *xpeak = s->GetPositionX();
165                   TMath::Sort(nfound, xpeak, index,down);
166                   Float_t xp = xpeak[index[0]];
167                   Double_t hmax = xp+3*sigma;
168                   Double_t hmin = xp-3*sigma;
169                   hCFD->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
170                 }
171                 else
172                   {
173                     hCFD->Rebin(2);
174                     TSpectrum *s1 = new TSpectrum(2*npeaks,1);
175                     nfound = s1->Search(hCFD,sigma," ",0.1);
176                     if(nfound!=0){
177                       Float_t *xpeak = s1->GetPositionX();
178                       TMath::Sort(nfound, xpeak, index,down);
179                       Float_t xp = xpeak[index[0]];
180                       Double_t hmax = xp+3*sigma;
181                       Double_t hmin = xp-3*sigma;
182                       hCFD->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
183                     }
184                     else 
185                       ok=false;
186                   }
187
188                 if (im == 0) cfd0[i] = hCFD->GetMean();
189                 y1[im] =  hCFD->GetMean() - cfd0[i];
190               }
191               if( hQTC) x1[im] = hQTC->GetMean();
192               
193               if( hLED){
194                 TSpectrum *s = new TSpectrum(2*npeaks,1);
195                 nfound = s->Search(hLED,sigma," ",0.1);
196                 if(nfound!=0){
197                   Float_t *xpeak = s->GetPositionX();
198                   TMath::Sort(nfound, xpeak, index,down);
199                   Float_t xp = xpeak[index[0]];
200                   Double_t hmax = xp+10*sigma;
201                   Double_t hmin = xp-10*sigma;
202                   hLED->GetXaxis()->SetRangeUser(hmin-10,hmax+10);
203                 }
204                 else 
205                   ok=false;
206                 x2[im] = hLED->GetMean();
207                 xx2[im] = x2[nmips-im-1]; 
208               }
209               xx[im]=mips[im];
210               
211               if (hQTC) delete  hQTC;
212               if (hCFD) delete  hCFD;
213               if (hLED) delete  hLED;
214             }
215           
216           for (Int_t imi=0; imi<nmips; imi++)
217             {
218               yy1[imi] = Float_t (mips[nmips-imi-1]);
219               xx1[imi] = x2[nmips-imi-1]; 
220               //              cout<<" LED rev "<<i<<" "<<imi<<" "<<xx1[imi]<<" "<< yy1[imi]<<"nmips-imi-1 = "<< nmips-imi-1<<endl;
221             }
222           
223           if(i==0) cout<<"Making graphs..."<<endl;
224                   
225           TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
226           grwalkqtc->SetTitle(Form("PMT%i",i));
227           TGraph *grwalkled = new TGraph (nmips,x2,y1);
228           grwalkled->SetTitle(Form("PMT%i",i));
229           fWalk.AddAtAndExpand(grwalkqtc,i);
230           fAmpLEDRec.AddAtAndExpand(grwalkled,i);
231           //      cout<<" add walk "<<i<<endl;
232          
233           //fit amplitude graphs to make comparison wth new one   
234           TGraph *grampled = new TGraph (nmips,xx1,yy1);
235           TGraph *grqtc = new TGraph (nmips,x1,xx);
236           fQTC.AddAtAndExpand(grqtc,i);  
237           fAmpLED.AddAtAndExpand(grampled,i);
238           //      cout<<" add amp "<<i<<endl;
239
240           if(i==23)
241             cout<<"Graphs created..."<<endl;
242         }
243     } //if gFile exits
244   return ok;
245 }
246
247
248
249