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