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