Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibBase.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
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  Base class for the calibration components using 
20 //  as input TPCseeds and ESDs
21 //  Event loop outside of the component
22 //
23 //
24 // Base functionality to be implemeneted by component 
25 /* 
26    //In some cases only one of this function to be implemented
27    virtual void     Process(AliESDEvent *event)
28    virtual void     Process(AliTPCseed *track)
29    //
30    virtual Long64_t Merge(TCollection *li);
31    virtual void     Analyze()
32    void             Terminate();
33 */
34 // Functionality provided by base class for Algorith debuging:
35 //  TTreeSRedirector * cstream =  GetDebugStreamer() - get debug streamer which can be use for numerical debugging
36 //                      
37
38
39
40 //  marian.ivanov@cern.ch
41 // 
42 #include "AliTPCcalibBase.h"
43 #include "TSystem.h"
44 #include "TFile.h"
45 #include "TTreeStream.h"
46 #include "TTimeStamp.h"
47 #include "TGraph.h"
48 #include "TGraphErrors.h"
49 #include "TF1.h"
50 #include "TH1.h"
51 #include "THnSparse.h"
52 #include "TH1D.h"
53 #include "TH2D.h"
54 #include "TAxis.h"
55 #include "AliRecoParam.h"
56
57
58 #include "AliLog.h"
59 #include "AliESDEvent.h"
60
61
62 ClassImp(AliTPCcalibBase)
63
64 AliTPCcalibBase::AliTPCcalibBase():
65     TNamed(),
66     fDebugStreamer(0),
67     fStreamLevel(0),   
68     fRun(0),                  //!  current Run number
69     fEvent(0),                //!  current Event number
70     fTime(0),                 //!  current Time
71     fTrigger(0),              //! current trigger type
72     fMagF(0),                 //! current magnetic field
73     fTriggerMaskReject(-1),   //trigger mask - reject trigger
74     fTriggerMaskAccept(-1),   //trigger mask - accept trigger
75     fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
76     fRejectLaser(kTRUE),                 //flag- reject laser
77     fTriggerClass(),
78     fCurrentEvent(0),           //! current event
79     fCurrentTrack(0),           //! current esd track
80     fCurrentFriendTrack(0),           //! current esd track
81     fCurrentSeed(0),            //! current seed
82     fDebugLevel(0)
83 {
84   //
85   // Constructor
86   //
87 }
88
89 AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
90   TNamed(name,title),
91   fDebugStreamer(0),
92   fStreamLevel(0),   
93   fRun(0),                  //!  current Run number
94   fEvent(0),                //!  current Event number
95   fTime(0),                 //!  current Time
96   fTrigger(0),              //! current trigger type
97   fMagF(0),                 //! current magnetic field
98   fTriggerMaskReject(-1),   //trigger mask - reject trigger
99   fTriggerMaskAccept(-1),   //trigger mask - accept trigger
100   fHasLaser(kFALSE),                    //flag the laser is overlayed with given event 
101   fRejectLaser(kTRUE),                 //flag- reject laser
102   fTriggerClass(),
103   fCurrentEvent(0),           //! current event
104   fCurrentTrack(0),           //! current esd track
105   fCurrentFriendTrack(0),           //! current esd track
106   fCurrentSeed(0),            //! current seed
107   fDebugLevel(0)
108 {
109   //
110   // Constructor
111   //
112 }
113
114 AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
115   TNamed(calib),
116   fDebugStreamer(0),
117   fStreamLevel(calib.fStreamLevel),
118   fRun(0),                  //!  current Run number
119   fEvent(0),                //!  current Event number
120   fTime(0),                 //!  current Time
121   fTrigger(0),              //! current trigger type
122   fMagF(0),                 //! current magnetic field
123   fTriggerMaskReject(calib.fTriggerMaskReject),   //trigger mask - reject trigger
124   fTriggerMaskAccept(calib.fTriggerMaskAccept),   //trigger mask - accept trigger
125   fHasLaser(calib.fHasLaser),                    //flag the laser is overlayed with given event
126   fRejectLaser(calib.fRejectLaser),                 //flag- reject laser
127   fTriggerClass(calib.fTriggerClass),
128   fCurrentEvent(0),           //! current event
129   fCurrentTrack(0),           //! current esd track
130   fCurrentFriendTrack(0),           //! current esd track
131   fCurrentSeed(0),            //! current seed
132   fDebugLevel(calib.fDebugLevel)
133 {
134   //
135   // copy constructor
136   //
137 }
138
139 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
140   //
141   // operator=
142   //
143   ((TNamed *)this)->operator=(calib);
144   fDebugStreamer=0;
145   fStreamLevel=calib.fStreamLevel;
146   fDebugLevel=calib.fDebugLevel;
147   return *this;
148 }
149
150
151 AliTPCcalibBase::~AliTPCcalibBase() {
152   //
153   // destructor
154   //
155   if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
156   if (fDebugStreamer) delete fDebugStreamer;
157   fDebugStreamer=0;
158 }
159
160 void  AliTPCcalibBase::Terminate(){
161   //
162   //
163   //
164   if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
165   if (fDebugStreamer) delete fDebugStreamer;
166   fDebugStreamer = 0;
167   return;
168 }
169
170 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
171   //
172   // Get Debug streamer
173   // In case debug streamer not yet initialized and StreamLevel>0 create new one
174   //
175   if (fStreamLevel==0) return 0;
176   if (fDebugStreamer) return fDebugStreamer;
177   TString dsName;
178   dsName=GetName();
179   dsName+="Debug.root";
180   dsName.ReplaceAll(" ",""); 
181   fDebugStreamer = new TTreeSRedirector(dsName.Data());
182   return fDebugStreamer;
183 }
184
185
186 void    AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
187   //
188   //
189   //
190   fRun     = event->GetRunNumber();
191   fEvent   = event->GetEventNumberInFile();
192   fTime    = event->GetTimeStamp();
193   fTrigger = event->GetTriggerMask();
194   fMagF    = event->GetMagneticField();
195   fTriggerClass = event->GetFiredTriggerClasses().Data();
196   fHasLaser = HasLaser(event); 
197   
198 }
199
200
201 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
202   //
203   //
204   //
205   // Use event specie
206   Bool_t isLaser=kFALSE;
207   UInt_t specie = event->GetEventSpecie();  // select only cosmic events
208   if (specie==AliRecoParam::kCalib) {
209     isLaser = kTRUE;
210   }
211   return isLaser;
212 }
213
214
215
216 Bool_t AliTPCcalibBase::AcceptTrigger(){
217   //
218   // Apply trigger mask - Don't do calibration for non proper triggers
219   // 
220   if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
221   if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
222   if (fHasLaser && fRejectLaser) return kFALSE;
223   return kTRUE;
224 }
225
226
227 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
228   //
229   // store  - copy debug output to the destination position
230   // currently ONLY for local copy
231   if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
232   if (fStreamLevel==0) return;
233   TString dsName;
234   dsName=GetName();
235   dsName+="Debug.root";
236   dsName.ReplaceAll(" ",""); 
237   TString dsName2=path;
238   gSystem->MakeDirectory(dsName2.Data());
239   dsName2+=gSystem->HostName();
240   gSystem->MakeDirectory(dsName2.Data());
241   dsName2+="/";
242   TTimeStamp s;
243   dsName2+=Int_t(s.GetNanoSec());
244   dsName2+="/";
245   gSystem->MakeDirectory(dsName2.Data());
246   dsName2+=dsName;
247   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
248   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
249   TFile::Cp(dsName.Data(),dsName2.Data());
250 }
251
252
253
254 TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp, Bool_t useMedian, TTreeSRedirector *cstream, Int_t ival){
255   //
256   // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
257   // 
258   TF1 funcGaus("funcGaus","gaus");
259   TH2D * hist = h->Projection(axisDim1, axisDim2);
260   Double_t *xvec = new Double_t[hist->GetNbinsX()];
261   Double_t *yvec = new Double_t[hist->GetNbinsX()];
262   Double_t *xerr = new Double_t[hist->GetNbinsX()];
263   Double_t *yerr = new Double_t[hist->GetNbinsX()];
264   Int_t counter  = 0;
265   TH1D * projectionHist =0;
266   //
267
268   for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
269     Int_t nsum=0;
270     Int_t imin   =  i;
271     Int_t imax   =  i;
272
273     for (Int_t idelta=0; idelta<nmaxBin; idelta++){
274       //
275       imin   =  TMath::Max(i-idelta,1);
276       imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
277       nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
278       if (nsum==0) break;
279       if (nsum>minEntries) break;
280     }
281     if (nsum<minEntries) continue;
282     //
283     hist->GetXaxis()->SetRange(imin,imax);
284     projectionHist = hist->ProjectionY("gain",imin,imax);
285     // Determine Median:
286     Float_t xMin = projectionHist->GetXaxis()->GetXmin();
287     Float_t xMax = projectionHist->GetXaxis()->GetXmax();
288     Float_t xMedian = (xMin+xMax)*0.5;
289     Float_t integral = 0;
290     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
291       integral+=projectionHist->GetBinContent(jbin);
292     }
293     //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
294     //
295     //
296     Float_t currentSum=0;
297     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
298       currentSum += projectionHist->GetBinContent(jbin);
299       if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
300       if (currentSum<fracUp*integral)  xMax = projectionHist->GetBinCenter(jbin+1);      
301       if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
302         xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
303                    projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
304           (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
305       }
306     }
307     //
308     Float_t rms  = projectionHist->GetRMS();
309     //    i += interval;
310     //
311     Double_t xcenter =  hist->GetMean(); 
312     Double_t xrms    =  hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.); 
313     Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
314     if (rms>0){
315       // cut on +- 4 RMS
316       projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
317       Double_t chi2 = funcGaus.GetChisquare();
318       //  
319       xvec[counter] = xcenter;
320       yvec[counter] = funcGaus.GetParameter(ival);
321       xerr[counter] = xrms;
322       yerr[counter] = funcGaus.GetParError(ival); 
323       if (useMedian) yvec[counter] = xMedian;
324       if (cstream){
325         (*cstream)<<"fitDebug"<<
326           "xcenter="<<xcenter<<
327           "xMin="<<xMin<<
328           "xMax="<<xMax<<
329           "xMedian="<<xMedian<<
330           "xFitM"<<yvec[counter]<<
331           "xFitS"<<yerr[counter]<<
332           "chi2="<<chi2<<         
333           "\n";
334       }
335       counter++;
336     }else{
337       xvec[counter] = xcenter;
338       yvec[counter] = xMedian;
339       xerr[counter] = xrms;
340       yerr[counter] = binWidth/TMath::Sqrt(12.); 
341       counter++;
342     }
343     delete projectionHist;
344   }
345   
346   TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
347   delete [] xvec;
348   delete [] yvec;
349   delete [] xerr;
350   delete [] yerr;
351   delete hist;
352   return graphErrors;
353 }
354
355
356 void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
357
358   // Method for the correct logarithmic binning of histograms
359
360   TAxis *axis = h->GetAxis(axisDim);
361   int bins = axis->GetNbins();
362
363   Double_t from = axis->GetXmin();
364   Double_t to = axis->GetXmax();
365   Double_t *new_bins = new Double_t[bins + 1];
366
367   new_bins[0] = from;
368   Double_t factor = pow(to/from, 1./bins);
369
370   for (int i = 1; i <= bins; i++) {
371    new_bins[i] = factor * new_bins[i-1];
372   }
373   axis->Set(bins, new_bins);
374   delete [] new_bins;
375
376 }
377 void AliTPCcalibBase::BinLogX(TH1 *h) {
378
379   // Method for the correct logarithmic binning of histograms
380
381   TAxis *axis = h->GetXaxis();
382   int bins = axis->GetNbins();
383
384   Double_t from = axis->GetXmin();
385   Double_t to = axis->GetXmax();
386   Double_t *new_bins = new Double_t[bins + 1];
387
388   new_bins[0] = from;
389   Double_t factor = pow(to/from, 1./bins);
390
391   for (int i = 1; i <= bins; i++) {
392    new_bins[i] = factor * new_bins[i-1];
393   }
394   axis->Set(bins, new_bins);
395   delete [] new_bins;
396
397 }