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"
55 #include "AliRecoParam.h"
59 #include "AliESDEvent.h"
62 ClassImp(AliTPCcalibBase)
64 AliTPCcalibBase::AliTPCcalibBase():
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
73 fTriggerMaskReject(-1), //trigger mask - reject trigger
74 fTriggerMaskAccept(-1), //trigger mask - accept trigger
75 fHasLaser(kFALSE), //flag the laser is overlayed with given event
76 fRejectLaser(kTRUE), //flag- reject laser
78 fCurrentEvent(0), //! current event
79 fCurrentTrack(0), //! current esd track
80 fCurrentFriendTrack(0), //! current esd track
81 fCurrentSeed(0), //! current seed
89 AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
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
103 fCurrentEvent(0), //! current event
104 fCurrentTrack(0), //! current esd track
105 fCurrentFriendTrack(0), //! current esd track
106 fCurrentSeed(0), //! current seed
114 AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
117 fStreamLevel(calib.fStreamLevel),
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
123 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
124 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
125 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
126 fRejectLaser(calib.fRejectLaser), //flag- reject laser
127 fTriggerClass(calib.fTriggerClass),
128 fCurrentEvent(0), //! current event
129 fCurrentTrack(0), //! current esd track
130 fCurrentFriendTrack(0), //! current esd track
131 fCurrentSeed(0), //! current seed
132 fDebugLevel(calib.fDebugLevel)
139 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
143 if (this == &calib) return (*this);
145 ((TNamed *)this)->operator=(calib);
147 fStreamLevel=calib.fStreamLevel;
148 fDebugLevel=calib.fDebugLevel;
153 AliTPCcalibBase::~AliTPCcalibBase() {
157 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
158 if (fDebugStreamer) delete fDebugStreamer;
162 void AliTPCcalibBase::Terminate(){
166 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
167 if (fDebugStreamer) delete fDebugStreamer;
172 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
174 // Get Debug streamer
175 // In case debug streamer not yet initialized and StreamLevel>0 create new one
177 if (fStreamLevel==0) return 0;
178 if (fDebugStreamer) return fDebugStreamer;
181 dsName+="Debug.root";
182 dsName.ReplaceAll(" ","");
183 fDebugStreamer = new TTreeSRedirector(dsName.Data());
184 return fDebugStreamer;
188 void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
192 fRun = event->GetRunNumber();
193 fEvent = event->GetEventNumberInFile();
194 fTime = event->GetTimeStamp();
195 fTrigger = event->GetTriggerMask();
196 fMagF = event->GetMagneticField();
197 fTriggerClass = event->GetFiredTriggerClasses().Data();
198 fHasLaser = HasLaser(event);
203 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
208 Bool_t isLaser=kFALSE;
209 UInt_t specie = event->GetEventSpecie(); // select only cosmic events
210 if (specie==AliRecoParam::kCalib) {
218 Bool_t AliTPCcalibBase::AcceptTrigger(){
220 // Apply trigger mask - Don't do calibration for non proper triggers
222 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
223 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
224 if (fHasLaser && fRejectLaser) return kFALSE;
229 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
231 // store - copy debug output to the destination position
232 // currently ONLY for local copy
233 if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
234 if (fStreamLevel==0) return;
237 dsName+="Debug.root";
238 dsName.ReplaceAll(" ","");
239 TString dsName2=path;
240 gSystem->MakeDirectory(dsName2.Data());
241 dsName2+=gSystem->HostName();
242 gSystem->MakeDirectory(dsName2.Data());
245 dsName2+=Int_t(s.GetNanoSec());
247 gSystem->MakeDirectory(dsName2.Data());
249 AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
250 printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
251 TFile::Cp(dsName.Data(),dsName2.Data());
256 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){
258 // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
260 TF1 funcGaus("funcGaus","gaus");
261 TH2D * hist = h->Projection(axisDim1, axisDim2);
262 Double_t *xvec = new Double_t[hist->GetNbinsX()];
263 Double_t *yvec = new Double_t[hist->GetNbinsX()];
264 Double_t *xerr = new Double_t[hist->GetNbinsX()];
265 Double_t *yerr = new Double_t[hist->GetNbinsX()];
267 TH1D * projectionHist =0;
270 for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
275 for (Int_t idelta=0; idelta<nmaxBin; idelta++){
277 imin = TMath::Max(i-idelta,1);
278 imax = TMath::Min(i+idelta,hist->GetNbinsX());
279 nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
281 if (nsum>minEntries) break;
283 if (nsum<minEntries) continue;
285 hist->GetXaxis()->SetRange(imin,imax);
286 projectionHist = hist->ProjectionY("gain",imin,imax);
288 Float_t xMin = projectionHist->GetXaxis()->GetXmin();
289 Float_t xMax = projectionHist->GetXaxis()->GetXmax();
290 Float_t xMedian = (xMin+xMax)*0.5;
291 Float_t integral = 0;
292 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
293 integral+=projectionHist->GetBinContent(jbin);
295 //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
298 Float_t currentSum=0;
299 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
300 currentSum += projectionHist->GetBinContent(jbin);
301 if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
302 if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
303 if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
304 xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
305 projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
306 (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
310 Float_t rms = projectionHist->GetRMS();
313 Double_t xcenter = hist->GetMean();
314 Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
315 Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
318 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
319 Double_t chi2 = funcGaus.GetChisquare();
321 xvec[counter] = xcenter;
322 yvec[counter] = funcGaus.GetParameter(ival);
323 xerr[counter] = xrms;
324 yerr[counter] = funcGaus.GetParError(ival);
325 if (useMedian) yvec[counter] = xMedian;
327 (*cstream)<<"fitDebug"<<
328 "xcenter="<<xcenter<<
331 "xMedian="<<xMedian<<
332 "xFitM"<<yvec[counter]<<
333 "xFitS"<<yerr[counter]<<
339 xvec[counter] = xcenter;
340 yvec[counter] = xMedian;
341 xerr[counter] = xrms;
342 yerr[counter] = binWidth/TMath::Sqrt(12.);
345 delete projectionHist;
348 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
358 void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
360 // Method for the correct logarithmic binning of histograms
362 TAxis *axis = h->GetAxis(axisDim);
363 int bins = axis->GetNbins();
365 Double_t from = axis->GetXmin();
366 Double_t to = axis->GetXmax();
367 Double_t *new_bins = new Double_t[bins + 1];
370 Double_t factor = pow(to/from, 1./bins);
372 for (int i = 1; i <= bins; i++) {
373 new_bins[i] = factor * new_bins[i-1];
375 axis->Set(bins, new_bins);
379 void AliTPCcalibBase::BinLogX(TH1 *h) {
381 // Method for the correct logarithmic binning of histograms
383 TAxis *axis = h->GetXaxis();
384 int bins = axis->GetNbins();
386 Double_t from = axis->GetXmin();
387 Double_t to = axis->GetXmax();
388 Double_t *new_bins = new Double_t[bins + 1];
391 Double_t factor = pow(to/from, 1./bins);
393 for (int i = 1; i <= bins; i++) {
394 new_bins[i] = factor * new_bins[i-1];
396 axis->Set(bins, new_bins);
401 void AliTPCcalibBase::BinLogX(TAxis *axis) {
403 // Method for the correct logarithmic binning of histograms
405 Int_t bins = axis->GetNbins();
407 Double_t from = axis->GetXmin();
408 Double_t to = axis->GetXmax();
410 Double_t *new_bins = new Double_t[bins + 1];
413 Double_t factor = pow(to/from, 1./bins);
415 for (int i = 1; i <= bins; i++) {
416 new_bins[i] = factor * new_bins[i-1];
418 axis->Set(bins, new_bins);