During simulation: fill STU region w/ non null time sums
[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   if (this == &calib) return (*this);
144
145   ((TNamed *)this)->operator=(calib);
146   fDebugStreamer=0;
147   fStreamLevel=calib.fStreamLevel;
148   fDebugLevel=calib.fDebugLevel;
149   return *this;
150 }
151
152
153 AliTPCcalibBase::~AliTPCcalibBase() {
154   //
155   // destructor
156   //
157   if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
158   if (fDebugStreamer) delete fDebugStreamer;
159   fDebugStreamer=0;
160 }
161
162 void  AliTPCcalibBase::Terminate(){
163   //
164   //
165   //
166   if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
167   if (fDebugStreamer) delete fDebugStreamer;
168   fDebugStreamer = 0;
169   return;
170 }
171
172 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
173   //
174   // Get Debug streamer
175   // In case debug streamer not yet initialized and StreamLevel>0 create new one
176   //
177   if (fStreamLevel==0) return 0;
178   if (fDebugStreamer) return fDebugStreamer;
179   TString dsName;
180   dsName=GetName();
181   dsName+="Debug.root";
182   dsName.ReplaceAll(" ",""); 
183   fDebugStreamer = new TTreeSRedirector(dsName.Data());
184   return fDebugStreamer;
185 }
186
187
188 void    AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
189   //
190   //
191   //
192   fRun     = event->GetRunNumber();
193   fEvent   = event->GetEventNumberInFile();
194   fTime    = event->GetTimeStamp();
195   fTrigger = event->GetTriggerMask();
196   fMagF    = event->GetMagneticField();
197   fTriggerClass = event->GetFiredTriggerClasses().Data();
198   fHasLaser = HasLaser(event); 
199   
200 }
201
202
203 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
204   //
205   //
206   //
207   // Use event specie
208   Bool_t isLaser=kFALSE;
209   UInt_t specie = event->GetEventSpecie();  // select only cosmic events
210   if (specie==AliRecoParam::kCalib) {
211     isLaser = kTRUE;
212   }
213   return isLaser;
214 }
215
216
217
218 Bool_t AliTPCcalibBase::AcceptTrigger(){
219   //
220   // Apply trigger mask - Don't do calibration for non proper triggers
221   // 
222   if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
223   if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
224   if (fHasLaser && fRejectLaser) return kFALSE;
225   return kTRUE;
226 }
227
228
229 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
230   //
231   // store  - copy debug output to the destination position
232   // currently ONLY for local copy
233   if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
234   if (fStreamLevel==0) return;
235   TString dsName;
236   dsName=GetName();
237   dsName+="Debug.root";
238   dsName.ReplaceAll(" ",""); 
239   TString dsName2=path;
240   gSystem->MakeDirectory(dsName2.Data());
241   dsName2+=gSystem->HostName();
242   gSystem->MakeDirectory(dsName2.Data());
243   dsName2+="/";
244   TTimeStamp s;
245   dsName2+=Int_t(s.GetNanoSec());
246   dsName2+="/";
247   gSystem->MakeDirectory(dsName2.Data());
248   dsName2+=dsName;
249   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
250   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
251   TFile::Cp(dsName.Data(),dsName2.Data());
252 }
253
254
255
256 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){
257   //
258   // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
259   // 
260   TF1 funcGaus("funcGaus","gaus");
261   TH2D * hist = h->Projection(axisDim1, axisDim2);
262   Double_t *xvec = new Double_t[hist->GetNbinsX()];
263   Double_t *yvec = new Double_t[hist->GetNbinsX()];
264   Double_t *xerr = new Double_t[hist->GetNbinsX()];
265   Double_t *yerr = new Double_t[hist->GetNbinsX()];
266   Int_t counter  = 0;
267   TH1D * projectionHist =0;
268   //
269
270   for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
271     Int_t nsum=0;
272     Int_t imin   =  i;
273     Int_t imax   =  i;
274
275     for (Int_t idelta=0; idelta<nmaxBin; idelta++){
276       //
277       imin   =  TMath::Max(i-idelta,1);
278       imax   =  TMath::Min(i+idelta,hist->GetNbinsX());
279       nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
280       if (nsum==0) break;
281       if (nsum>minEntries) break;
282     }
283     if (nsum<minEntries) continue;
284     //
285     hist->GetXaxis()->SetRange(imin,imax);
286     projectionHist = hist->ProjectionY("gain",imin,imax);
287     // Determine Median:
288     Float_t xMin = projectionHist->GetXaxis()->GetXmin();
289     Float_t xMax = projectionHist->GetXaxis()->GetXmax();
290     Float_t xMedian = (xMin+xMax)*0.5;
291     Float_t integral = 0;
292     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
293       integral+=projectionHist->GetBinContent(jbin);
294     }
295     //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
296     //
297     //
298     Float_t currentSum=0;
299     for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
300       currentSum += projectionHist->GetBinContent(jbin);
301       if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
302       if (currentSum<fracUp*integral)  xMax = projectionHist->GetBinCenter(jbin+1);      
303       if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
304         xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
305                    projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
306           (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
307       }
308     }
309     //
310     Float_t rms  = projectionHist->GetRMS();
311     //    i += interval;
312     //
313     Double_t xcenter =  hist->GetMean(); 
314     Double_t xrms    =  hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.); 
315     Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
316     if (rms>0){
317       // cut on +- 4 RMS
318       projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
319       Double_t chi2 = funcGaus.GetChisquare();
320       //  
321       xvec[counter] = xcenter;
322       yvec[counter] = funcGaus.GetParameter(ival);
323       xerr[counter] = xrms;
324       yerr[counter] = funcGaus.GetParError(ival); 
325       if (useMedian) yvec[counter] = xMedian;
326       if (cstream){
327         (*cstream)<<"fitDebug"<<
328           "xcenter="<<xcenter<<
329           "xMin="<<xMin<<
330           "xMax="<<xMax<<
331           "xMedian="<<xMedian<<
332           "xFitM"<<yvec[counter]<<
333           "xFitS"<<yerr[counter]<<
334           "chi2="<<chi2<<         
335           "\n";
336       }
337       counter++;
338     }else{
339       xvec[counter] = xcenter;
340       yvec[counter] = xMedian;
341       xerr[counter] = xrms;
342       yerr[counter] = binWidth/TMath::Sqrt(12.); 
343       counter++;
344     }
345     delete projectionHist;
346   }
347   
348   TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
349   delete [] xvec;
350   delete [] yvec;
351   delete [] xerr;
352   delete [] yerr;
353   delete hist;
354   return graphErrors;
355 }
356
357
358 void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
359
360   // Method for the correct logarithmic binning of histograms
361
362   TAxis *axis = h->GetAxis(axisDim);
363   int bins = axis->GetNbins();
364
365   Double_t from = axis->GetXmin();
366   Double_t to = axis->GetXmax();
367   Double_t *new_bins = new Double_t[bins + 1];
368
369   new_bins[0] = from;
370   Double_t factor = pow(to/from, 1./bins);
371
372   for (int i = 1; i <= bins; i++) {
373    new_bins[i] = factor * new_bins[i-1];
374   }
375   axis->Set(bins, new_bins);
376   delete [] new_bins;
377
378 }
379 void AliTPCcalibBase::BinLogX(TH1 *h) {
380
381   // Method for the correct logarithmic binning of histograms
382
383   TAxis *axis = h->GetXaxis();
384   int bins = axis->GetNbins();
385
386   Double_t from = axis->GetXmin();
387   Double_t to = axis->GetXmax();
388   Double_t *new_bins = new Double_t[bins + 1];
389
390   new_bins[0] = from;
391   Double_t factor = pow(to/from, 1./bins);
392
393   for (int i = 1; i <= bins; i++) {
394    new_bins[i] = factor * new_bins[i-1];
395   }
396   axis->Set(bins, new_bins);
397   delete [] new_bins;
398
399 }