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"
52 #include "THnSparse.h"
56 #include "AliESDEvent.h"
59 ClassImp(AliTPCcalibBase)
61 AliTPCcalibBase::AliTPCcalibBase():
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
70 fTriggerMaskReject(-1), //trigger mask - reject trigger
71 fTriggerMaskAccept(-1), //trigger mask - accept trigger
72 fHasLaser(kFALSE), //flag the laser is overlayed with given event
73 fRejectLaser(kTRUE), //flag- reject laser
75 fCurrentEvent(0), //! current event
76 fCurrentTrack(0), //! current esd track
77 fCurrentSeed(0), //! current seed
85 AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
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
99 fCurrentEvent(0), //! current event
100 fCurrentTrack(0), //! current esd track
101 fCurrentSeed(0), //! current seed
109 AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
112 fStreamLevel(calib.fStreamLevel),
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
118 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
119 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
120 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
121 fRejectLaser(calib.fRejectLaser), //flag- reject laser
122 fTriggerClass(calib.fTriggerClass),
123 fCurrentEvent(0), //! current event
124 fCurrentTrack(0), //! current esd track
125 fCurrentSeed(0), //! current seed
126 fDebugLevel(calib.fDebugLevel)
133 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
137 ((TNamed *)this)->operator=(calib);
139 fStreamLevel=calib.fStreamLevel;
140 fDebugLevel=calib.fDebugLevel;
145 AliTPCcalibBase::~AliTPCcalibBase() {
149 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
150 if (fDebugStreamer) delete fDebugStreamer;
154 void AliTPCcalibBase::Terminate(){
158 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
159 if (fDebugStreamer) delete fDebugStreamer;
164 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
166 // Get Debug streamer
167 // In case debug streamer not yet initialized and StreamLevel>0 create new one
169 if (fStreamLevel==0) return 0;
170 if (fDebugStreamer) return fDebugStreamer;
173 dsName+="Debug.root";
174 dsName.ReplaceAll(" ","");
175 fDebugStreamer = new TTreeSRedirector(dsName.Data());
176 return fDebugStreamer;
180 void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
184 fRun = event->GetRunNumber();
185 fEvent = event->GetEventNumberInFile();
186 fTime = event->GetTimeStamp();
187 fTrigger = event->GetTriggerMask();
188 fMagF = event->GetMagneticField();
189 fTriggerClass = event->GetFiredTriggerClasses().Data();
190 fHasLaser = HasLaser(event);
195 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
199 // Thresholds more than 8 tracks with small dip angle
201 const Int_t kMinLaserTracks = 8;
202 const Float_t kThrLaser = 0.3;
203 const Float_t kLaserTgl = 0.01;
205 Int_t ntracks = event->GetNumberOfTracks();
206 if (ntracks<kMinLaserTracks) return kFALSE;
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;
214 if (TMath::Abs(track->GetTgl())<kLaserTgl) nlaser++;
216 if (nlaser>kMinLaserTracks) return kTRUE;
217 if (nall>0 && nlaser/nall>kThrLaser) return kTRUE;
223 Bool_t AliTPCcalibBase::AcceptTrigger(){
225 // Apply trigger mask - Don't do calibration for non proper triggers
227 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
228 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
229 if (fHasLaser && fRejectLaser) return kFALSE;
234 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
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;
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());
250 dsName2+=Int_t(s.GetNanoSec());
252 gSystem->MakeDirectory(dsName2.Data());
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());
261 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){
263 // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
265 TF1 funcGaus("funcGaus","gaus");
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()];
272 TH1D * projectionHist =0;
275 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
280 for (Int_t idelta=0; idelta<nmaxBin; idelta++){
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()));
286 if (nsum>minEntries) break;
288 if (nsum<minEntries) continue;
290 hist->GetXaxis()->SetRange(imin,imax);
291 projectionHist = hist->ProjectionY("gain",imin,imax);
293 Float_t xMin = projectionHist->GetXaxis()->GetXmin();
294 Float_t xMax = projectionHist->GetXaxis()->GetXmax();
295 Float_t xMedian = (xMin+xMax)*0.5;
296 Float_t integral = 0;
297 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
298 integral+=projectionHist->GetBinContent(jbin);
300 //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
303 Float_t currentSum=0;
304 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
305 currentSum += projectionHist->GetBinContent(jbin);
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){
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));
315 Float_t rms = projectionHist->GetRMS();
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);
323 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
324 Double_t chi2 = funcGaus.GetChisquare();
326 xvec[counter] = xcenter;
327 yvec[counter] = funcGaus.GetParameter(ival);
328 xerr[counter] = xrms;
329 yerr[counter] = funcGaus.GetParError(ival);
330 if (useMedian) yvec[counter] = xMedian;
332 (*cstream)<<"fitDebug"<<
333 "xcenter="<<xcenter<<
336 "xMedian="<<xMedian<<
337 "xFitM"<<yvec[counter]<<
338 "xFitS"<<yerr[counter]<<
344 xvec[counter] = xcenter;
345 yvec[counter] = xMedian;
346 xerr[counter] = xrms;
347 yerr[counter] = binWidth/TMath::Sqrt(12.);
350 delete projectionHist;
353 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);