adding some cuts for systematic study
[u/mrichter/AliRoot.git] / T0 / AliT0CalibTimeEq.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18
19 ///////////////////////////////////////////////////////////////////////////////
20 //                                                                           //
21 // class for T0 calibration                       TM-AC-AM_6-02-2006  
22 // equalize time shift for each time CFD channel
23 //                                                                           //
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include "AliT0CalibTimeEq.h"
27 #include "AliLog.h"
28 #include <TFile.h>
29 #include <TMath.h>
30 #include <TF1.h>
31 #include <TSpectrum.h>
32 #include <TProfile.h>
33 #include <iostream>
34
35 ClassImp(AliT0CalibTimeEq)
36
37 //________________________________________________________________
38   AliT0CalibTimeEq::AliT0CalibTimeEq():TNamed(),
39                                        fMeanVertex(0),        
40                                        fRmsVertex(0)      
41 {
42   //
43   for(Int_t i=0; i<24; i++) {
44    fTimeEq[i] = 0;            // Time Equalized for OCDB         
45    fTimeEqRms[i] = -1;        // RMS of Time Equalized for OCDB  
46    for (Int_t ih=0; ih<5; ih++)   fCFDvalue[i][ih] = 0;
47   }
48 }
49
50 //________________________________________________________________
51 AliT0CalibTimeEq::AliT0CalibTimeEq(const char* name):TNamed(),
52                                        fMeanVertex(0),        
53                                        fRmsVertex(0)      
54 {
55   //constructor
56
57   TString namst = "Calib_";
58   namst += name;
59   SetName(namst.Data());
60   SetTitle(namst.Data());
61   for(Int_t i=0; i<24; i++) {
62    fTimeEq[i] = 0;            // Time Equalized for OCDB         
63    fTimeEqRms[i] = -1;        // RMS of Time Equalized for OCDB  
64    for (Int_t ih=0; ih<5; ih++)   fCFDvalue[i][ih] = 0;
65   }
66 }
67
68 //________________________________________________________________
69 AliT0CalibTimeEq::AliT0CalibTimeEq(const AliT0CalibTimeEq& calibda):TNamed(calibda),            
70                                        fMeanVertex(0),        
71                                        fRmsVertex(0)      
72 {
73 // copy constructor
74   SetName(calibda.GetName());
75   SetTitle(calibda.GetName());
76   ((AliT0CalibTimeEq &) calibda).Copy(*this);
77
78 }
79
80 //________________________________________________________________
81 AliT0CalibTimeEq &AliT0CalibTimeEq::operator =(const AliT0CalibTimeEq& calibda)
82 {
83 // assignment operator
84   SetName(calibda.GetName());
85   SetTitle(calibda.GetName());
86  
87   if (this != &calibda) (( AliT0CalibTimeEq &) calibda).Copy(*this);
88   return *this;
89 }
90
91 //________________________________________________________________
92 AliT0CalibTimeEq::~AliT0CalibTimeEq()
93 {
94   //
95   // destrictor
96 }
97 //________________________________________________________________
98 void AliT0CalibTimeEq::Reset()
99 {
100   //reset values
101
102   memset(fCFDvalue,0,120*sizeof(Float_t));
103   memset(fTimeEq,1,24*sizeof(Float_t));
104 }
105
106
107 //________________________________________________________________
108 void  AliT0CalibTimeEq::Print(Option_t*) const
109 {
110   // print time values
111
112   printf("\n    ----    PM Arrays       ----\n\n");
113   printf(" Time delay CFD \n");
114   for (Int_t i=0; i<24; i++) 
115     printf(" CFD  %f diff %f  \n",fCFDvalue[i][0],fTimeEq[i]);
116
117
118
119 //________________________________________________________________
120 Bool_t AliT0CalibTimeEq::ComputeOnlineParams(const char* filePhys)
121 {
122   // compute online equalized time
123   Float_t meandiff, sigmadiff, meanver, meancfdtime, sigmacfdtime;
124   meandiff = sigmadiff =  meanver = meancfdtime = sigmacfdtime =0;
125   Double_t rmsver=0;
126   Int_t nent=0;
127   Int_t okdiff=0;
128   Int_t oktime=0;
129   Bool_t ok=false;
130   //
131   gFile = TFile::Open(filePhys);
132   if(!gFile) {
133     AliError("No input PHYS data found ");
134   }
135   else
136     {
137       ok=true;
138       for (Int_t i=0; i<24; i++)
139         {
140           meandiff = sigmadiff =  meanver = meancfdtime = sigmacfdtime =0;
141           TH1F *cfd     = (TH1F*) gFile->Get(Form("CFD1minCFD%d",i+1));
142           TH1F *cfdtime = (TH1F*) gFile->Get(Form("CFD%d",i+1));
143           if(!cfd) {
144             AliWarning(Form("no Diff histograms collected by PHYS DA for channel %i", i));
145             okdiff++;
146             if(okdiff<4) {
147               meandiff = 0;
148               sigmadiff = 0;
149             }
150             else
151               ok = false; 
152           }
153           if(!cfdtime) {
154             AliWarning(Form("no CFD histograms collected by PHYS DA for channel %i", i));
155             oktime++;
156             if(oktime<4) {
157               meancfdtime = 0;
158             }
159             else
160               ok = false; 
161           }
162           if(cfd) {
163             nent = Int_t(cfd->GetEntries());
164             if( nent<=50) {
165               okdiff++;
166               //              printf(" pmt %i nent %i cfd->GetRMS() %f cfd->GetMean() %f \n",
167               //     i, nent, cfd->GetRMS(), cfd->GetMean() );
168               if(okdiff<4) {
169                 meandiff = 0;
170                 sigmadiff = 0;
171               }
172               else
173                 {
174                   // printf(" OK fsle:: pmt %i nent %i cfd->GetRMS() %f cfd->GetMean() %f \n",
175                   // i, nent, cfd->GetRMS(), cfd->GetMean() );
176                   AliWarning(Form(" Not  enouph data in PMT %i- PMT1:  %i ", i, nent));
177                   ok = false; 
178                 }
179             }
180             if(nent>50)  { 
181               if(cfd->GetRMS()>1.5 )
182                 GetMeanAndSigma(cfd, meandiff, sigmadiff);
183               if(cfd->GetRMS()<=1.5) 
184                 {
185                   meandiff = cfd->GetMean();
186                   sigmadiff=cfd->GetRMS();
187                 }
188               Int_t   maxBin = cfd->GetMaximumBin(); 
189               Double_t  meanEstimate = cfd->GetBinCenter( maxBin); 
190               if(TMath::Abs(meanEstimate - meandiff) > 20 ) meandiff = meanEstimate; 
191             }
192           }
193           
194           if(cfdtime) {
195             nent = Int_t(cfdtime->GetEntries());
196             if( nent<=50  ) {
197               oktime++;
198               if(oktime<4) {
199                 meancfdtime = 0;
200               }
201               else
202                 {
203                   AliWarning(Form(" Not  enouph data in PMT %i:  %i ", i, nent));
204                   ok = false; 
205                 }
206             }
207             if(nent > 50  )  { //!!!!!!!!!!
208               if(cfdtime->GetRMS()>1.5 )
209                 GetMeanAndSigma(cfdtime,meancfdtime, sigmacfdtime);
210               if(cfdtime->GetRMS()<=1.5) 
211                 {
212                   if(cfdtime->GetRMS()==0 ||cfdtime->GetMean()==0 ) {
213                     ok = false;
214                   }
215                   meancfdtime = cfdtime->GetMean();
216                   sigmacfdtime = cfdtime->GetRMS();
217                 }
218               Int_t   maxBin = cfdtime->GetMaximumBin(); 
219               Double_t  meanEstimate = cfdtime->GetBinCenter( maxBin); 
220               if(TMath::Abs(meanEstimate - meancfdtime) > 20 ) meancfdtime = meanEstimate; 
221             }
222           }
223           SetTimeEq(i,meandiff);
224           SetTimeEqRms(i,sigmadiff);
225           SetCFDvalue(i,0,meancfdtime);
226           if (cfd) delete cfd;
227           if (cfdtime) delete cfdtime;
228         }
229       TH1F *ver = (TH1F*) gFile->Get("hVertex") ;
230       if(ver) {
231         meanver = ver->GetMean();
232         rmsver = ver->GetRMS();
233       }
234       SetMeanVertex(meanver);
235       SetRmsVertex(rmsver);
236       gFile->Close();
237       delete gFile;
238     }
239     return ok; 
240 }
241
242 //________________________________________________________________
243 Int_t AliT0CalibTimeEq::ComputeOfflineParams(const char* filePhys, Float_t *timecdb, Float_t *cfdvalue, Int_t badpmt)
244 {
245   // fStatus implementation:
246   //ok means
247   // 1000  - can not open file 
248   // for timediff
249   // 20 >3 histos are empty
250   // -11  low statistics oe empty histos in few channels; for these channels OCDB value will be written back WARNING
251   // for cfd
252   // 20 >2 histos are empty or with low statistic 
253   // -11 if less 3 channels are empty or low statistics in few channels ;for these channels OCDB value will be written back WARNING
254    //
255   // compute offline equalized time
256   Float_t meandiff, sigmadiff, meanver, meancfdtime, sigmacfdtime;
257   meandiff = sigmadiff =  meanver = meancfdtime = sigmacfdtime =0;
258   Int_t nent=0;
259   Int_t ok = 0;
260   Int_t okcfd=0;
261   TH1F *cfddiff = NULL; 
262   TH1F *cfdtime = NULL;
263   TObjArray * tzeroObj = NULL;
264
265   gFile = TFile::Open(filePhys);
266   if(!gFile) {
267     AliError("No input PHYS data found ");
268     ok = 1000;
269     return ok;
270   }
271   else
272     {
273       meandiff = sigmadiff =  meanver = meancfdtime = sigmacfdtime =0;
274       //      TDirectory *dr = (TDirectory*) gFile->Get("T0Calib");
275       tzeroObj = dynamic_cast<TObjArray*>(gFile->Get("T0Calib"));
276       for (Int_t i=0; i<24; i++)
277         {
278           if (i != badpmt) {        
279             if(tzeroObj) {
280               cfddiff = (TH1F*) tzeroObj->FindObject(Form("CFD1minCFD%d",i+1));
281               cfdtime = (TH1F*)tzeroObj->FindObject(Form("CFD%d",i+1));
282             }
283             else
284               {
285                 cfddiff = (TH1F*)gFile->Get(Form("CFD1minCFD%d",i+1));
286                 cfdtime = (TH1F*)gFile->Get(Form("CFD%d",i+1));
287                 
288               }
289             if(!cfddiff ) {
290               AliWarning(Form("no  histograms collected by pass0 for diff channel %i\n", i));    
291                 meandiff = timecdb[i];
292                 sigmadiff = 0;
293               }
294             if(cfddiff) {
295               nent = Int_t(cfddiff->GetEntries());
296               if ( nent == 0  ) {
297                 AliWarning(Form("no  entries in histogram for diff channel %i\n", i));
298                 meandiff = timecdb[i];
299                 sigmadiff = 0;
300               }
301               if(nent<=100 && nent>0) { //!!!!!
302                 AliWarning(Form(" Not  enouph data in PMT %i- PMT1:  %i \n", i, nent));
303                 meandiff = timecdb[i];
304                 sigmadiff = 0; 
305               }
306               if(nent>=100 )  { //!!!!!
307                 {
308                   if(cfddiff->GetRMS()>1.5 )
309                     GetMeanAndSigma(cfddiff, meandiff, sigmadiff);
310                   if(cfddiff->GetRMS()<=1.5) 
311                     {
312                       meandiff = cfddiff->GetMean();
313                       sigmadiff = cfddiff->GetRMS();
314                     }
315                   Int_t   maxBin = cfddiff->GetMaximumBin(); 
316                   Double_t  meanEstimate = cfddiff->GetBinCenter( maxBin); 
317                   if(TMath::Abs(meanEstimate - meandiff) > 20 ) meandiff = meanEstimate; 
318                 }
319               }     
320             }
321             if(!cfdtime ) {
322               AliWarning(Form("no  histograms collected by pass0 for time channel %i", i));
323               meancfdtime = cfdvalue[i];
324               okcfd++;
325               if(okcfd<2) {
326                 meancfdtime = cfdvalue[i];
327                 //      ok = -11;
328               }
329               else {
330                 AliError(Form("no  histograms collected by pass0 for time %i channels ", okcfd));
331                 return 20; 
332               }
333             }
334             if(cfdtime) {
335               nent = Int_t(cfdtime->GetEntries());
336               if (nent == 0 || 
337                   (cfdtime->GetRMS() == 0 || cfdtime->GetMean() ==0 ) ) 
338                 {
339                   okcfd++;
340                   if(okcfd<2) {
341                     meancfdtime = cfdvalue[i];
342                     //   ok = -11;
343                     printf("!!!!bad data:: pmt %i nent%i RMS %f mean %f cdbtime %f \n",
344                            i, nent, cfdtime->GetRMS(), cfdtime->GetMean(), cfdvalue[i] );
345                   }
346                   else
347                     {
348                     printf("!!!!fatal data:: pmt %i nent%i RMS %f mean %f cdbtime %f \n",
349                            i, nent, cfdtime->GetRMS(), cfdtime->GetMean(), cfdvalue[i]);
350                       AliError(Form(" histograms collected by pass0 for time %i channels are empty", okcfd));
351                       return 20; 
352                     }
353                 }
354               
355                   
356               if(nent<=100 && nent>0 ) 
357                 {
358                   okcfd++;
359                 AliWarning(Form(" Not  enouph data in PMT in CFD peak %i - %i ", i, nent));
360                   meancfdtime = cfdvalue[i];
361                   ok = -11;
362                   printf("!!!!low statstics:: pmt %i nent%i RMS %f mean %f cdbtime %f \n",
363                          i, nent, cfdtime->GetRMS(), cfdtime->GetMean(),  cfdvalue[i]);
364                   if (okcfd>3) return ok;
365                 }
366                       
367               if( nent>100 )    { //!!!!!
368                 if(cfdtime->GetRMS()>1.5 )
369                   GetMeanAndSigma(cfdtime,meancfdtime, sigmacfdtime);
370                 if(cfdtime->GetRMS()<=1.5) 
371                   {
372                     meancfdtime = cfdtime->GetMean();
373                     sigmacfdtime=cfdtime->GetRMS();
374                   }
375                 Int_t   maxBin = cfdtime->GetMaximumBin(); 
376                 Double_t  meanEstimate = cfdtime->GetBinCenter( maxBin); 
377                 if(TMath::Abs(meanEstimate - meancfdtime) > 20 ) meancfdtime = meanEstimate; 
378             }
379           }
380           
381           SetTimeEq(i,meandiff);
382           SetTimeEqRms(i,sigmadiff);
383           SetCFDvalue(i,0, meancfdtime );
384           if (cfddiff) cfddiff->Reset();
385           if (cfdtime) cfdtime->Reset();
386           } //bad pmt
387         }      
388       gFile->Close();
389       delete gFile;
390     }
391     return ok; 
392    }
393
394 //________________________________________________________________________
395 void AliT0CalibTimeEq::GetMeanAndSigma(TH1F* hist,  Float_t &mean, Float_t &sigma) {
396   
397   const double window = 5.;  //fit window 
398   double meanEstimate, sigmaEstimate; 
399   int maxBin;
400   maxBin        =  hist->GetMaximumBin(); //position of maximum
401   meanEstimate  =  hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
402   sigmaEstimate = hist->GetRMS();
403   // sigmaEstimate = 10;
404   TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
405   fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
406   hist->Fit("fit","QR","");
407
408   mean  = (Float_t) fit->GetParameter(1);
409   sigma = (Float_t) fit->GetParameter(2);
410
411   delete fit;
412 }
413
414