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