New directory for macros
[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(),
ae28e92e 75 fDebugLevel(0)
3ab35fc5 76{
77 //
78 // Constructor
79 //
80}
81
8cf12c96 82AliTPCcalibBase::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
15e48021 95 fTriggerClass(),
8cf12c96 96 fDebugLevel(0)
97{
98 //
99 // Constructor
100 //
101}
102
ae28e92e 103AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
104 TNamed(calib),
105 fDebugStreamer(0),
106 fStreamLevel(calib.fStreamLevel),
e5cbd73f 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
c51653e8 112 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
113 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
2d22fa36 114 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
115 fRejectLaser(calib.fRejectLaser), //flag- reject laser
fe8454d9 116 fTriggerClass(calib.fTriggerClass),
ae28e92e 117 fDebugLevel(calib.fDebugLevel)
118{
119 //
120 // copy constructor
121 //
122}
123
124AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
125 //
126 //
127 //
128 ((TNamed *)this)->operator=(calib);
129 fDebugStreamer=0;
130 fStreamLevel=calib.fStreamLevel;
131 fDebugLevel=calib.fDebugLevel;
525062ee 132 return *this;
ae28e92e 133}
134
135
3ab35fc5 136AliTPCcalibBase::~AliTPCcalibBase() {
137 //
138 // destructor
139 //
f7f33dec 140 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
ae28e92e 141 if (fDebugStreamer) delete fDebugStreamer;
142 fDebugStreamer=0;
143}
144
145void AliTPCcalibBase::Terminate(){
146 //
147 //
148 //
f7f33dec 149 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
ae28e92e 150 if (fDebugStreamer) delete fDebugStreamer;
151 fDebugStreamer = 0;
152 return;
153}
154
155TTreeSRedirector *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;
3ab35fc5 168}
f7f33dec 169
170
108953e9 171void 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();
15e48021 180 fTriggerClass = event->GetFiredTriggerClasses().Data();
2d22fa36 181 fHasLaser = HasLaser(event);
15e48021 182
108953e9 183}
184
2d22fa36 185
186Bool_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
c51653e8 214Bool_t AliTPCcalibBase::AcceptTrigger(){
215 //
216 // Apply trigger mask - Don't do calibration for non proper triggers
217 //
2d22fa36 218 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
219 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
220 if (fHasLaser && fRejectLaser) return kFALSE;
c51653e8 221 return kTRUE;
222}
223
224
f7f33dec 225void 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+="/";
47328338 240 TTimeStamp s;
241 dsName2+=Int_t(s.GetNanoSec());
242 dsName2+="/";
243 gSystem->MakeDirectory(dsName2.Data());
f7f33dec 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}
52a5850f 249
250
251
252TGraphErrors * 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;
d0b31f57 286 Float_t integral = 0;
287 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
288 integral+=projectionHist->GetBinContent(jbin);
289 }
290 printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
52a5850f 291 //
292 //
293 Float_t currentSum=0;
294 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
295 currentSum += projectionHist->GetBinContent(jbin);
d0b31f57 296 if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
297 if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
298 if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
52a5850f 299 xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
300 projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
301 (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
302 }
303 }
304 //
305 Float_t rms = projectionHist->GetRMS();
306 // i += interval;
307 //
308 Double_t xcenter = hist->GetMean();
309 Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
310 Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
311 if (rms>0){
312 TF1 funcGaus("funcGaus","gaus");
313 // cut on +- 4 RMS
314 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
315 //
316 xvec[counter] = xcenter;
317 yvec[counter] = funcGaus.GetParameter(1);
318 xerr[counter] = xrms;
319 yerr[counter] = funcGaus.GetParError(1);
320 counter++;
321 }else{
322 xvec[counter] = xcenter;
323 yvec[counter] = xMedian;
324 xerr[counter] = xrms;
325 yerr[counter] = binWidth/TMath::Sqrt(12.);
326 counter++;
327 }
328 delete projectionHist;
329 }
330
331 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
332 delete [] xvec;
333 delete [] yvec;
334 delete [] xerr;
335 delete [] yerr;
336 delete hist;
337 return graphErrors;
338}