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