Increase version number
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibBase.cxx
CommitLineData
3ab35fc5 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//
ae28e92e 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
3ab35fc5 41//
42#include "AliTPCcalibBase.h"
f7f33dec 43#include "TSystem.h"
44#include "TFile.h"
ae28e92e 45#include "TTreeStream.h"
47328338 46#include "TTimeStamp.h"
52a5850f 47#include "TGraph.h"
48#include "TGraphErrors.h"
49#include "TF1.h"
3226e48f 50#include "TH1.h"
51#include "THnSparse.h"
52a5850f 52#include "TH1D.h"
53#include "TH2D.h"
3226e48f 54#include "TAxis.h"
52a5850f 55
56
57#include "AliLog.h"
108953e9 58#include "AliESDEvent.h"
f7f33dec 59
3ab35fc5 60
61ClassImp(AliTPCcalibBase)
62
ae28e92e 63AliTPCcalibBase::AliTPCcalibBase():
64 TNamed(),
65 fDebugStreamer(0),
108953e9 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
c51653e8 72 fTriggerMaskReject(-1), //trigger mask - reject trigger
73 fTriggerMaskAccept(-1), //trigger mask - accept trigger
2d22fa36 74 fHasLaser(kFALSE), //flag the laser is overlayed with given event
75 fRejectLaser(kTRUE), //flag- reject laser
15e48021 76 fTriggerClass(),
b842d904 77 fCurrentEvent(0), //! current event
78 fCurrentTrack(0), //! current esd track
4486a91f 79 fCurrentFriendTrack(0), //! current esd track
b842d904 80 fCurrentSeed(0), //! current seed
ae28e92e 81 fDebugLevel(0)
3ab35fc5 82{
83 //
84 // Constructor
85 //
86}
87
8cf12c96 88AliTPCcalibBase::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
15e48021 101 fTriggerClass(),
b842d904 102 fCurrentEvent(0), //! current event
103 fCurrentTrack(0), //! current esd track
4486a91f 104 fCurrentFriendTrack(0), //! current esd track
b842d904 105 fCurrentSeed(0), //! current seed
8cf12c96 106 fDebugLevel(0)
107{
108 //
109 // Constructor
110 //
111}
112
ae28e92e 113AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
114 TNamed(calib),
115 fDebugStreamer(0),
116 fStreamLevel(calib.fStreamLevel),
e5cbd73f 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
c51653e8 122 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
123 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
2d22fa36 124 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
125 fRejectLaser(calib.fRejectLaser), //flag- reject laser
fe8454d9 126 fTriggerClass(calib.fTriggerClass),
b842d904 127 fCurrentEvent(0), //! current event
128 fCurrentTrack(0), //! current esd track
4486a91f 129 fCurrentFriendTrack(0), //! current esd track
b842d904 130 fCurrentSeed(0), //! current seed
ae28e92e 131 fDebugLevel(calib.fDebugLevel)
132{
133 //
134 // copy constructor
135 //
136}
137
138AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
139 //
4486a91f 140 // operator=
ae28e92e 141 //
142 ((TNamed *)this)->operator=(calib);
143 fDebugStreamer=0;
144 fStreamLevel=calib.fStreamLevel;
145 fDebugLevel=calib.fDebugLevel;
525062ee 146 return *this;
ae28e92e 147}
148
149
3ab35fc5 150AliTPCcalibBase::~AliTPCcalibBase() {
151 //
152 // destructor
153 //
f7f33dec 154 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
ae28e92e 155 if (fDebugStreamer) delete fDebugStreamer;
156 fDebugStreamer=0;
157}
158
159void AliTPCcalibBase::Terminate(){
160 //
161 //
162 //
f7f33dec 163 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
ae28e92e 164 if (fDebugStreamer) delete fDebugStreamer;
165 fDebugStreamer = 0;
166 return;
167}
168
169TTreeSRedirector *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;
3ab35fc5 182}
f7f33dec 183
184
108953e9 185void 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();
15e48021 194 fTriggerClass = event->GetFiredTriggerClasses().Data();
2d22fa36 195 fHasLaser = HasLaser(event);
15e48021 196
108953e9 197}
198
2d22fa36 199
200Bool_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
c51653e8 228Bool_t AliTPCcalibBase::AcceptTrigger(){
229 //
230 // Apply trigger mask - Don't do calibration for non proper triggers
231 //
2d22fa36 232 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
233 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
234 if (fHasLaser && fRejectLaser) return kFALSE;
c51653e8 235 return kTRUE;
236}
237
238
f7f33dec 239void 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+="/";
47328338 254 TTimeStamp s;
255 dsName2+=Int_t(s.GetNanoSec());
256 dsName2+="/";
257 gSystem->MakeDirectory(dsName2.Data());
f7f33dec 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}
52a5850f 263
264
265
70b0f0bf 266TGraphErrors * 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){
52a5850f 267 //
268 // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
269 //
3d549531 270 TF1 funcGaus("funcGaus","gaus");
52a5850f 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;
d0b31f57 301 Float_t integral = 0;
302 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
303 integral+=projectionHist->GetBinContent(jbin);
304 }
7f766801 305 //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
52a5850f 306 //
307 //
308 Float_t currentSum=0;
309 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
310 currentSum += projectionHist->GetBinContent(jbin);
d0b31f57 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){
52a5850f 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){
52a5850f 327 // cut on +- 4 RMS
328 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
3d549531 329 Double_t chi2 = funcGaus.GetChisquare();
52a5850f 330 //
331 xvec[counter] = xcenter;
70b0f0bf 332 yvec[counter] = funcGaus.GetParameter(ival);
52a5850f 333 xerr[counter] = xrms;
70b0f0bf 334 yerr[counter] = funcGaus.GetParError(ival);
3d549531 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 }
52a5850f 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}
3226e48f 366
367
368void 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}
389void 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}