Fixing warning.
[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"
50#include "TH1D.h"
51#include "TH2D.h"
52#include "THnSparse.h"
53
54
55#include "AliLog.h"
108953e9 56#include "AliESDEvent.h"
f7f33dec 57
3ab35fc5 58
59ClassImp(AliTPCcalibBase)
60
ae28e92e 61AliTPCcalibBase::AliTPCcalibBase():
62 TNamed(),
63 fDebugStreamer(0),
108953e9 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
c51653e8 70 fTriggerMaskReject(-1), //trigger mask - reject trigger
71 fTriggerMaskAccept(-1), //trigger mask - accept trigger
2d22fa36 72 fHasLaser(kFALSE), //flag the laser is overlayed with given event
73 fRejectLaser(kTRUE), //flag- reject laser
15e48021 74 fTriggerClass(),
b842d904 75 fCurrentEvent(0), //! current event
76 fCurrentTrack(0), //! current esd track
77 fCurrentSeed(0), //! current seed
ae28e92e 78 fDebugLevel(0)
3ab35fc5 79{
80 //
81 // Constructor
82 //
83}
84
8cf12c96 85AliTPCcalibBase::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
15e48021 98 fTriggerClass(),
b842d904 99 fCurrentEvent(0), //! current event
100 fCurrentTrack(0), //! current esd track
101 fCurrentSeed(0), //! current seed
8cf12c96 102 fDebugLevel(0)
103{
104 //
105 // Constructor
106 //
107}
108
ae28e92e 109AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
110 TNamed(calib),
111 fDebugStreamer(0),
112 fStreamLevel(calib.fStreamLevel),
e5cbd73f 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
c51653e8 118 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
119 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
2d22fa36 120 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
121 fRejectLaser(calib.fRejectLaser), //flag- reject laser
fe8454d9 122 fTriggerClass(calib.fTriggerClass),
b842d904 123 fCurrentEvent(0), //! current event
124 fCurrentTrack(0), //! current esd track
125 fCurrentSeed(0), //! current seed
ae28e92e 126 fDebugLevel(calib.fDebugLevel)
127{
128 //
129 // copy constructor
130 //
131}
132
133AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
134 //
135 //
136 //
137 ((TNamed *)this)->operator=(calib);
138 fDebugStreamer=0;
139 fStreamLevel=calib.fStreamLevel;
140 fDebugLevel=calib.fDebugLevel;
525062ee 141 return *this;
ae28e92e 142}
143
144
3ab35fc5 145AliTPCcalibBase::~AliTPCcalibBase() {
146 //
147 // destructor
148 //
f7f33dec 149 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
ae28e92e 150 if (fDebugStreamer) delete fDebugStreamer;
151 fDebugStreamer=0;
152}
153
154void AliTPCcalibBase::Terminate(){
155 //
156 //
157 //
f7f33dec 158 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
ae28e92e 159 if (fDebugStreamer) delete fDebugStreamer;
160 fDebugStreamer = 0;
161 return;
162}
163
164TTreeSRedirector *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;
3ab35fc5 177}
f7f33dec 178
179
108953e9 180void 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();
15e48021 189 fTriggerClass = event->GetFiredTriggerClasses().Data();
2d22fa36 190 fHasLaser = HasLaser(event);
15e48021 191
108953e9 192}
193
2d22fa36 194
195Bool_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
c51653e8 223Bool_t AliTPCcalibBase::AcceptTrigger(){
224 //
225 // Apply trigger mask - Don't do calibration for non proper triggers
226 //
2d22fa36 227 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
228 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
229 if (fHasLaser && fRejectLaser) return kFALSE;
c51653e8 230 return kTRUE;
231}
232
233
f7f33dec 234void 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+="/";
47328338 249 TTimeStamp s;
250 dsName2+=Int_t(s.GetNanoSec());
251 dsName2+="/";
252 gSystem->MakeDirectory(dsName2.Data());
f7f33dec 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}
52a5850f 258
259
260
70b0f0bf 261TGraphErrors * 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 262 //
263 // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
264 //
3d549531 265 TF1 funcGaus("funcGaus","gaus");
52a5850f 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;
d0b31f57 296 Float_t integral = 0;
297 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
298 integral+=projectionHist->GetBinContent(jbin);
299 }
7f766801 300 //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
52a5850f 301 //
302 //
303 Float_t currentSum=0;
304 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
305 currentSum += projectionHist->GetBinContent(jbin);
d0b31f57 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){
52a5850f 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){
52a5850f 322 // cut on +- 4 RMS
323 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
3d549531 324 Double_t chi2 = funcGaus.GetChisquare();
52a5850f 325 //
326 xvec[counter] = xcenter;
70b0f0bf 327 yvec[counter] = funcGaus.GetParameter(ival);
52a5850f 328 xerr[counter] = xrms;
70b0f0bf 329 yerr[counter] = funcGaus.GetParError(ival);
3d549531 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 }
52a5850f 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}