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