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