1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Base class for the calibration components using
20 // as input TPCseeds and ESDs
21 // Event loop outside of the component
24 // Base functionality to be implemeneted by component
26 //In some cases only one of this function to be implemented
27 virtual void Process(AliESDEvent *event)
28 virtual void Process(AliTPCseed *track)
30 virtual Long64_t Merge(TCollection *li);
31 virtual void Analyze()
34 // Functionality provided by base class for Algorith debuging:
35 // TTreeSRedirector * cstream = GetDebugStreamer() - get debug streamer which can be use for numerical debugging
40 // marian.ivanov@cern.ch
42 #include "AliTPCcalibBase.h"
45 #include "TTreeStream.h"
46 #include "TTimeStamp.h"
48 #include "TGraphErrors.h"
51 #include "THnSparse.h"
58 #include "AliESDEvent.h"
61 ClassImp(AliTPCcalibBase)
63 AliTPCcalibBase::AliTPCcalibBase():
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
72 fTriggerMaskReject(-1), //trigger mask - reject trigger
73 fTriggerMaskAccept(-1), //trigger mask - accept trigger
74 fHasLaser(kFALSE), //flag the laser is overlayed with given event
75 fRejectLaser(kTRUE), //flag- reject laser
77 fCurrentEvent(0), //! current event
78 fCurrentTrack(0), //! current esd track
79 fCurrentSeed(0), //! current seed
87 AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
91 fRun(0), //! current Run number
92 fEvent(0), //! current Event number
93 fTime(0), //! current Time
94 fTrigger(0), //! current trigger type
95 fMagF(0), //! current magnetic field
96 fTriggerMaskReject(-1), //trigger mask - reject trigger
97 fTriggerMaskAccept(-1), //trigger mask - accept trigger
98 fHasLaser(kFALSE), //flag the laser is overlayed with given event
99 fRejectLaser(kTRUE), //flag- reject laser
101 fCurrentEvent(0), //! current event
102 fCurrentTrack(0), //! current esd track
103 fCurrentSeed(0), //! current seed
111 AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
114 fStreamLevel(calib.fStreamLevel),
115 fRun(0), //! current Run number
116 fEvent(0), //! current Event number
117 fTime(0), //! current Time
118 fTrigger(0), //! current trigger type
119 fMagF(0), //! current magnetic field
120 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
121 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
122 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
123 fRejectLaser(calib.fRejectLaser), //flag- reject laser
124 fTriggerClass(calib.fTriggerClass),
125 fCurrentEvent(0), //! current event
126 fCurrentTrack(0), //! current esd track
127 fCurrentSeed(0), //! current seed
128 fDebugLevel(calib.fDebugLevel)
135 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
139 ((TNamed *)this)->operator=(calib);
141 fStreamLevel=calib.fStreamLevel;
142 fDebugLevel=calib.fDebugLevel;
147 AliTPCcalibBase::~AliTPCcalibBase() {
151 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
152 if (fDebugStreamer) delete fDebugStreamer;
156 void AliTPCcalibBase::Terminate(){
160 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
161 if (fDebugStreamer) delete fDebugStreamer;
166 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
168 // Get Debug streamer
169 // In case debug streamer not yet initialized and StreamLevel>0 create new one
171 if (fStreamLevel==0) return 0;
172 if (fDebugStreamer) return fDebugStreamer;
175 dsName+="Debug.root";
176 dsName.ReplaceAll(" ","");
177 fDebugStreamer = new TTreeSRedirector(dsName.Data());
178 return fDebugStreamer;
182 void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
186 fRun = event->GetRunNumber();
187 fEvent = event->GetEventNumberInFile();
188 fTime = event->GetTimeStamp();
189 fTrigger = event->GetTriggerMask();
190 fMagF = event->GetMagneticField();
191 fTriggerClass = event->GetFiredTriggerClasses().Data();
192 fHasLaser = HasLaser(event);
197 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
201 // Thresholds more than 8 tracks with small dip angle
203 const Int_t kMinLaserTracks = 8;
204 const Float_t kThrLaser = 0.3;
205 const Float_t kLaserTgl = 0.01;
207 Int_t ntracks = event->GetNumberOfTracks();
208 if (ntracks<kMinLaserTracks) return kFALSE;
211 for (Int_t i=0;i<ntracks;++i) {
212 AliESDtrack *track=event->GetTrack(i);
213 if (!track) continue;
214 if (track->GetTPCNcls()<=0) continue;
216 if (TMath::Abs(track->GetTgl())<kLaserTgl) nlaser++;
218 if (nlaser>kMinLaserTracks) return kTRUE;
219 if (nall>0 && nlaser/nall>kThrLaser) return kTRUE;
225 Bool_t AliTPCcalibBase::AcceptTrigger(){
227 // Apply trigger mask - Don't do calibration for non proper triggers
229 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
230 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
231 if (fHasLaser && fRejectLaser) return kFALSE;
236 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
238 // store - copy debug output to the destination position
239 // currently ONLY for local copy
240 if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
241 if (fStreamLevel==0) return;
244 dsName+="Debug.root";
245 dsName.ReplaceAll(" ","");
246 TString dsName2=path;
247 gSystem->MakeDirectory(dsName2.Data());
248 dsName2+=gSystem->HostName();
249 gSystem->MakeDirectory(dsName2.Data());
252 dsName2+=Int_t(s.GetNanoSec());
254 gSystem->MakeDirectory(dsName2.Data());
256 AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
257 printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
258 TFile::Cp(dsName.Data(),dsName2.Data());
263 TGraphErrors * 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){
265 // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
267 TF1 funcGaus("funcGaus","gaus");
268 TH2D * hist = h->Projection(axisDim1, axisDim2);
269 Double_t *xvec = new Double_t[hist->GetNbinsX()];
270 Double_t *yvec = new Double_t[hist->GetNbinsX()];
271 Double_t *xerr = new Double_t[hist->GetNbinsX()];
272 Double_t *yerr = new Double_t[hist->GetNbinsX()];
274 TH1D * projectionHist =0;
277 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
282 for (Int_t idelta=0; idelta<nmaxBin; idelta++){
284 imin = TMath::Max(i-idelta,1);
285 imax = TMath::Min(i+idelta,hist->GetNbinsX());
286 nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
288 if (nsum>minEntries) break;
290 if (nsum<minEntries) continue;
292 hist->GetXaxis()->SetRange(imin,imax);
293 projectionHist = hist->ProjectionY("gain",imin,imax);
295 Float_t xMin = projectionHist->GetXaxis()->GetXmin();
296 Float_t xMax = projectionHist->GetXaxis()->GetXmax();
297 Float_t xMedian = (xMin+xMax)*0.5;
298 Float_t integral = 0;
299 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
300 integral+=projectionHist->GetBinContent(jbin);
302 //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
305 Float_t currentSum=0;
306 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
307 currentSum += projectionHist->GetBinContent(jbin);
308 if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
309 if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
310 if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
311 xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
312 projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
313 (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
317 Float_t rms = projectionHist->GetRMS();
320 Double_t xcenter = hist->GetMean();
321 Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
322 Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
325 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
326 Double_t chi2 = funcGaus.GetChisquare();
328 xvec[counter] = xcenter;
329 yvec[counter] = funcGaus.GetParameter(ival);
330 xerr[counter] = xrms;
331 yerr[counter] = funcGaus.GetParError(ival);
332 if (useMedian) yvec[counter] = xMedian;
334 (*cstream)<<"fitDebug"<<
335 "xcenter="<<xcenter<<
338 "xMedian="<<xMedian<<
339 "xFitM"<<yvec[counter]<<
340 "xFitS"<<yerr[counter]<<
346 xvec[counter] = xcenter;
347 yvec[counter] = xMedian;
348 xerr[counter] = xrms;
349 yerr[counter] = binWidth/TMath::Sqrt(12.);
352 delete projectionHist;
355 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
365 void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
367 // Method for the correct logarithmic binning of histograms
369 TAxis *axis = h->GetAxis(axisDim);
370 int bins = axis->GetNbins();
372 Double_t from = axis->GetXmin();
373 Double_t to = axis->GetXmax();
374 Double_t *new_bins = new Double_t[bins + 1];
377 Double_t factor = pow(to/from, 1./bins);
379 for (int i = 1; i <= bins; i++) {
380 new_bins[i] = factor * new_bins[i-1];
382 axis->Set(bins, new_bins);
386 void AliTPCcalibBase::BinLogX(TH1 *h) {
388 // Method for the correct logarithmic binning of histograms
390 TAxis *axis = h->GetXaxis();
391 int bins = axis->GetNbins();
393 Double_t from = axis->GetXmin();
394 Double_t to = axis->GetXmax();
395 Double_t *new_bins = new Double_t[bins + 1];
398 Double_t factor = pow(to/from, 1./bins);
400 for (int i = 1; i <= bins; i++) {
401 new_bins[i] = factor * new_bins[i-1];
403 axis->Set(bins, new_bins);