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
82 AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
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
103 AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
106 fStreamLevel(calib.fStreamLevel),
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
112 fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
113 fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
114 fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
115 fRejectLaser(calib.fRejectLaser), //flag- reject laser
116 fTriggerClass(calib.fTriggerClass),
117 fDebugLevel(calib.fDebugLevel)
124 AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
128 ((TNamed *)this)->operator=(calib);
130 fStreamLevel=calib.fStreamLevel;
131 fDebugLevel=calib.fDebugLevel;
136 AliTPCcalibBase::~AliTPCcalibBase() {
140 if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
141 if (fDebugStreamer) delete fDebugStreamer;
145 void AliTPCcalibBase::Terminate(){
149 if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
150 if (fDebugStreamer) delete fDebugStreamer;
155 TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
157 // Get Debug streamer
158 // In case debug streamer not yet initialized and StreamLevel>0 create new one
160 if (fStreamLevel==0) return 0;
161 if (fDebugStreamer) return fDebugStreamer;
164 dsName+="Debug.root";
165 dsName.ReplaceAll(" ","");
166 fDebugStreamer = new TTreeSRedirector(dsName.Data());
167 return fDebugStreamer;
171 void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
175 fRun = event->GetRunNumber();
176 fEvent = event->GetEventNumberInFile();
177 fTime = event->GetTimeStamp();
178 fTrigger = event->GetTriggerMask();
179 fMagF = event->GetMagneticField();
180 fTriggerClass = event->GetFiredTriggerClasses().Data();
181 fHasLaser = HasLaser(event);
186 Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
190 // Thresholds more than 8 tracks with small dip angle
192 const Int_t kMinLaserTracks = 8;
193 const Float_t kThrLaser = 0.3;
194 const Float_t kLaserTgl = 0.01;
196 Int_t ntracks = event->GetNumberOfTracks();
197 if (ntracks<kMinLaserTracks) return kFALSE;
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;
205 if (TMath::Abs(track->GetTgl())<kLaserTgl) nlaser++;
207 if (nlaser>kMinLaserTracks) return kTRUE;
208 if (nall>0 && nlaser/nall>kThrLaser) return kTRUE;
214 Bool_t AliTPCcalibBase::AcceptTrigger(){
216 // Apply trigger mask - Don't do calibration for non proper triggers
218 if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
219 if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
220 if (fHasLaser && fRejectLaser) return kFALSE;
225 void AliTPCcalibBase::RegisterDebugOutput(const char *path){
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;
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());
241 dsName2+=Int_t(s.GetNanoSec());
243 gSystem->MakeDirectory(dsName2.Data());
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());
252 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){
254 // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
256 TF1 funcGaus("funcGaus","gaus");
257 TH2D * hist = h->Projection(axisDim1, axisDim2);
258 Double_t *xvec = new Double_t[hist->GetNbinsX()];
259 Double_t *yvec = new Double_t[hist->GetNbinsX()];
260 Double_t *xerr = new Double_t[hist->GetNbinsX()];
261 Double_t *yerr = new Double_t[hist->GetNbinsX()];
263 TH1D * projectionHist =0;
266 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
271 for (Int_t idelta=0; idelta<nmaxBin; idelta++){
273 imin = TMath::Max(i-idelta,1);
274 imax = TMath::Min(i+idelta,hist->GetNbinsX());
275 nsum = TMath::Nint(hist->Integral(imin,imax,0,hist->GetNbinsY()));
277 if (nsum>minEntries) break;
279 if (nsum<minEntries) continue;
281 hist->GetXaxis()->SetRange(imin,imax);
282 projectionHist = hist->ProjectionY("gain",imin,imax);
284 Float_t xMin = projectionHist->GetXaxis()->GetXmin();
285 Float_t xMax = projectionHist->GetXaxis()->GetXmax();
286 Float_t xMedian = (xMin+xMax)*0.5;
287 Float_t integral = 0;
288 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
289 integral+=projectionHist->GetBinContent(jbin);
291 //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
294 Float_t currentSum=0;
295 for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
296 currentSum += projectionHist->GetBinContent(jbin);
297 if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
298 if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
299 if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
300 xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
301 projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
302 (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
306 Float_t rms = projectionHist->GetRMS();
309 Double_t xcenter = hist->GetMean();
310 Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
311 Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
314 projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
315 Double_t chi2 = funcGaus.GetChisquare();
317 xvec[counter] = xcenter;
318 yvec[counter] = funcGaus.GetParameter(1);
319 xerr[counter] = xrms;
320 yerr[counter] = funcGaus.GetParError(1);
321 if (useMedian) yvec[counter] = xMedian;
323 (*cstream)<<"fitDebug"<<
324 "xcenter="<<xcenter<<
327 "xMedian="<<xMedian<<
328 "xFitM"<<yvec[counter]<<
329 "xFitS"<<yerr[counter]<<
335 xvec[counter] = xcenter;
336 yvec[counter] = xMedian;
337 xerr[counter] = xrms;
338 yerr[counter] = binWidth/TMath::Sqrt(12.);
341 delete projectionHist;
344 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);