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