Line data Source code
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 : //
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
41 : //
42 : #include "AliTPCcalibBase.h"
43 : #include "TSystem.h"
44 : #include "TFile.h"
45 : #include "TTreeStream.h"
46 : #include "TTimeStamp.h"
47 : #include "TGraph.h"
48 : #include "TGraphErrors.h"
49 : #include "TF1.h"
50 : #include "TH1.h"
51 : #include "THnSparse.h"
52 : #include "TH1D.h"
53 : #include "TH2D.h"
54 : #include "TAxis.h"
55 : #include "AliRecoParam.h"
56 :
57 :
58 : #include "AliLog.h"
59 : #include "AliESDEvent.h"
60 :
61 :
62 6 : ClassImp(AliTPCcalibBase)
63 :
64 : AliTPCcalibBase::AliTPCcalibBase():
65 0 : TNamed(),
66 0 : fDebugStreamer(0),
67 0 : fStreamLevel(0),
68 0 : fRun(0), //! current Run number
69 0 : fEvent(0), //! current Event number
70 0 : fTime(0), //! current Time
71 0 : fTrigger(0), //! current trigger type
72 0 : fMagF(0), //! current magnetic field
73 0 : fTriggerMaskReject(-1), //trigger mask - reject trigger
74 0 : fTriggerMaskAccept(-1), //trigger mask - accept trigger
75 0 : fHasLaser(kFALSE), //flag the laser is overlayed with given event
76 0 : fRejectLaser(kTRUE), //flag- reject laser
77 0 : fTriggerClass(),
78 0 : fCurrentEvent(0), //! current event
79 0 : fCurrentTrack(0), //! current esd track
80 0 : fCurrentFriendTrack(0), //! current esd track
81 0 : fCurrentSeed(0), //! current seed
82 0 : fDebugLevel(0)
83 0 : {
84 : //
85 : // Constructor
86 : //
87 0 : }
88 :
89 : AliTPCcalibBase::AliTPCcalibBase(const char * name, const char * title):
90 0 : TNamed(name,title),
91 0 : fDebugStreamer(0),
92 0 : fStreamLevel(0),
93 0 : fRun(0), //! current Run number
94 0 : fEvent(0), //! current Event number
95 0 : fTime(0), //! current Time
96 0 : fTrigger(0), //! current trigger type
97 0 : fMagF(0), //! current magnetic field
98 0 : fTriggerMaskReject(-1), //trigger mask - reject trigger
99 0 : fTriggerMaskAccept(-1), //trigger mask - accept trigger
100 0 : fHasLaser(kFALSE), //flag the laser is overlayed with given event
101 0 : fRejectLaser(kTRUE), //flag- reject laser
102 0 : fTriggerClass(),
103 0 : fCurrentEvent(0), //! current event
104 0 : fCurrentTrack(0), //! current esd track
105 0 : fCurrentFriendTrack(0), //! current esd track
106 0 : fCurrentSeed(0), //! current seed
107 0 : fDebugLevel(0)
108 0 : {
109 : //
110 : // Constructor
111 : //
112 0 : }
113 :
114 : AliTPCcalibBase::AliTPCcalibBase(const AliTPCcalibBase&calib):
115 0 : TNamed(calib),
116 0 : fDebugStreamer(0),
117 0 : fStreamLevel(calib.fStreamLevel),
118 0 : fRun(0), //! current Run number
119 0 : fEvent(0), //! current Event number
120 0 : fTime(0), //! current Time
121 0 : fTrigger(0), //! current trigger type
122 0 : fMagF(0), //! current magnetic field
123 0 : fTriggerMaskReject(calib.fTriggerMaskReject), //trigger mask - reject trigger
124 0 : fTriggerMaskAccept(calib.fTriggerMaskAccept), //trigger mask - accept trigger
125 0 : fHasLaser(calib.fHasLaser), //flag the laser is overlayed with given event
126 0 : fRejectLaser(calib.fRejectLaser), //flag- reject laser
127 0 : fTriggerClass(calib.fTriggerClass),
128 0 : fCurrentEvent(0), //! current event
129 0 : fCurrentTrack(0), //! current esd track
130 0 : fCurrentFriendTrack(0), //! current esd track
131 0 : fCurrentSeed(0), //! current seed
132 0 : fDebugLevel(calib.fDebugLevel)
133 0 : {
134 : //
135 : // copy constructor
136 : //
137 0 : }
138 :
139 : AliTPCcalibBase &AliTPCcalibBase::operator=(const AliTPCcalibBase&calib){
140 : //
141 : // operator=
142 : //
143 0 : if (this == &calib) return (*this);
144 :
145 0 : ((TNamed *)this)->operator=(calib);
146 0 : fDebugStreamer=0;
147 0 : fStreamLevel=calib.fStreamLevel;
148 0 : fDebugLevel=calib.fDebugLevel;
149 0 : return *this;
150 0 : }
151 :
152 :
153 0 : AliTPCcalibBase::~AliTPCcalibBase() {
154 : //
155 : // destructor
156 : //
157 0 : if (fDebugLevel>0) printf("AliTPCcalibBase::~AliTPCcalibBase\n");
158 0 : if (fDebugStreamer) delete fDebugStreamer;
159 0 : fDebugStreamer=0;
160 0 : }
161 :
162 : void AliTPCcalibBase::Terminate(){
163 : //
164 : //
165 : //
166 0 : if (fDebugLevel>0) printf("AliTPCcalibBase::Terminate\n");
167 0 : if (fDebugStreamer) delete fDebugStreamer;
168 0 : fDebugStreamer = 0;
169 0 : return;
170 : }
171 :
172 : TTreeSRedirector *AliTPCcalibBase::GetDebugStreamer(){
173 : //
174 : // Get Debug streamer
175 : // In case debug streamer not yet initialized and StreamLevel>0 create new one
176 : //
177 0 : if (fStreamLevel==0) return 0;
178 0 : if (fDebugStreamer) return fDebugStreamer;
179 0 : TString dsName;
180 0 : dsName=GetName();
181 0 : dsName+="Debug.root";
182 0 : dsName.ReplaceAll(" ","");
183 0 : fDebugStreamer = new TTreeSRedirector(dsName.Data());
184 : return fDebugStreamer;
185 0 : }
186 :
187 :
188 : void AliTPCcalibBase::UpdateEventInfo(AliESDEvent * event){
189 : //
190 : //
191 : //
192 0 : fRun = event->GetRunNumber();
193 0 : fEvent = event->GetEventNumberInFile();
194 0 : fTime = event->GetTimeStamp();
195 0 : fTrigger = event->GetTriggerMask();
196 0 : fMagF = event->GetMagneticField();
197 0 : fTriggerClass = event->GetFiredTriggerClasses().Data();
198 0 : fHasLaser = HasLaser(event);
199 :
200 0 : }
201 :
202 :
203 : Bool_t AliTPCcalibBase::HasLaser(AliESDEvent *event){
204 : //
205 : //
206 : //
207 : // Use event specie
208 : Bool_t isLaser=kFALSE;
209 0 : UInt_t specie = event->GetEventSpecie(); // select only cosmic events
210 0 : if (specie==AliRecoParam::kCalib) {
211 : isLaser = kTRUE;
212 0 : }
213 0 : return isLaser;
214 : }
215 :
216 :
217 :
218 : Bool_t AliTPCcalibBase::AcceptTrigger(){
219 : //
220 : // Apply trigger mask - Don't do calibration for non proper triggers
221 : //
222 0 : if (fTriggerMaskReject==(Int_t)fTrigger) return kFALSE;
223 0 : if (fTriggerMaskAccept>0 && fTriggerMaskAccept!=(Int_t)fTrigger) return kFALSE;
224 0 : if (fHasLaser && fRejectLaser) return kFALSE;
225 0 : return kTRUE;
226 0 : }
227 :
228 :
229 : void AliTPCcalibBase::RegisterDebugOutput(const char *path){
230 : //
231 : // store - copy debug output to the destination position
232 : // currently ONLY for local copy
233 0 : if (fDebugLevel>0) printf("AliTPCcalibBase::RegisterDebugOutput(%s)\n",path);
234 0 : if (fStreamLevel==0) return;
235 0 : TString dsName;
236 0 : dsName=GetName();
237 0 : dsName+="Debug.root";
238 0 : dsName.ReplaceAll(" ","");
239 0 : TString dsName2=path;
240 0 : gSystem->MakeDirectory(dsName2.Data());
241 0 : dsName2+=gSystem->HostName();
242 0 : gSystem->MakeDirectory(dsName2.Data());
243 0 : dsName2+="/";
244 0 : TTimeStamp s;
245 0 : dsName2+=Int_t(s.GetNanoSec());
246 0 : dsName2+="/";
247 0 : gSystem->MakeDirectory(dsName2.Data());
248 0 : dsName2+=dsName;
249 0 : AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
250 0 : printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
251 0 : TFile::Cp(dsName.Data(),dsName2.Data());
252 0 : }
253 :
254 : 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){
255 : //
256 : // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
257 : //
258 0 : TH2D * hist = h->Projection(axisDim1, axisDim2);
259 0 : TGraphErrors *gr=FitSlices(hist, minEntries, nmaxBin, fracLow, fracUp, useMedian, cstream, ival);
260 0 : delete hist;
261 0 : return gr;
262 : }
263 :
264 : TGraphErrors * AliTPCcalibBase::FitSlices(TH2* hist, Int_t minEntries, Int_t nmaxBin, Float_t fracLow, Float_t fracUp, Bool_t useMedian, TTreeSRedirector *cstream, Int_t ival){
265 : //
266 : // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram
267 : //
268 : const Int_t kMinStat=100;
269 0 : TF1 funcGaus("funcGaus","gaus");
270 0 : Double_t *xvec = new Double_t[hist->GetNbinsX()];
271 0 : Double_t *yvec = new Double_t[hist->GetNbinsX()];
272 0 : Double_t *xerr = new Double_t[hist->GetNbinsX()];
273 0 : Double_t *yerr = new Double_t[hist->GetNbinsX()];
274 : Int_t counter = 0;
275 : TH1D * projectionHist =0;
276 : //
277 :
278 0 : for(Int_t i=1; i <= hist->GetNbinsX(); i++) {
279 : Int_t nsum=0;
280 : Int_t imin = i;
281 : Int_t imax = i;
282 :
283 0 : for (Int_t idelta=0; idelta<nmaxBin; idelta++){
284 : //
285 0 : imin = TMath::Max(i-idelta,1);
286 0 : imax = TMath::Min(i+idelta,hist->GetNbinsX());
287 0 : nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()));
288 0 : if (nsum==0) break;
289 0 : if (nsum>minEntries) break;
290 : }
291 0 : if (nsum<kMinStat) continue;
292 : //
293 0 : hist->GetXaxis()->SetRange(imin,imax);
294 0 : projectionHist = hist->ProjectionY("gain",imin,imax);
295 : // Determine Median:
296 0 : Float_t xMin = projectionHist->GetXaxis()->GetXmin();
297 0 : Float_t xMax = projectionHist->GetXaxis()->GetXmax();
298 0 : Float_t xMedian = (xMin+xMax)*0.5;
299 : Float_t integral = 0;
300 0 : for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
301 0 : integral+=projectionHist->GetBinContent(jbin);
302 : }
303 : //printf("Integral %f\t%f\n",integral, projectionHist->GetSum());
304 : //
305 : //
306 : Float_t currentSum=0;
307 0 : for(Int_t jbin=1; jbin<projectionHist->GetNbinsX()-1; jbin++) {
308 0 : currentSum += projectionHist->GetBinContent(jbin);
309 0 : if (currentSum<fracLow*integral) xMin = projectionHist->GetBinCenter(jbin);
310 0 : if (currentSum<fracUp*integral) xMax = projectionHist->GetBinCenter(jbin+1);
311 0 : if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){
312 0 : xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+
313 0 : projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/
314 0 : (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1));
315 0 : }
316 : }
317 : //
318 0 : Float_t rms = projectionHist->GetRMS();
319 : // i += interval;
320 : //
321 0 : Double_t xcenter = hist->GetMean();
322 0 : Double_t xrms = hist->GetRMS()+hist->GetXaxis()->GetBinWidth(1)/TMath::Sqrt(12.);
323 0 : Double_t binWidth = projectionHist->GetXaxis()->GetBinWidth(1);
324 0 : if (rms>0){
325 : // cut on +- 4 RMS
326 0 : projectionHist->Fit(&funcGaus,"QN","",xMin, xMax);
327 0 : Double_t chi2 = funcGaus.GetChisquare();
328 : //
329 0 : xvec[counter] = xcenter;
330 0 : yvec[counter] = funcGaus.GetParameter(ival);
331 0 : xerr[counter] = xrms;
332 0 : yerr[counter] = funcGaus.GetParError(ival);
333 0 : if (useMedian) yvec[counter] = xMedian;
334 0 : if (cstream){
335 0 : (*cstream)<<"fitDebug"<<
336 0 : "xcenter="<<xcenter<<
337 0 : "xMin="<<xMin<<
338 0 : "xMax="<<xMax<<
339 0 : "xMedian="<<xMedian<<
340 0 : "xFitM"<<yvec[counter]<<
341 0 : "xFitS"<<yerr[counter]<<
342 0 : "chi2="<<chi2<<
343 : "\n";
344 : }
345 0 : counter++;
346 0 : }else{
347 0 : xvec[counter] = xcenter;
348 0 : yvec[counter] = xMedian;
349 0 : xerr[counter] = xrms;
350 0 : yerr[counter] = binWidth/TMath::Sqrt(12.);
351 0 : counter++;
352 : }
353 0 : delete projectionHist;
354 0 : }
355 :
356 0 : TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
357 0 : delete [] xvec;
358 0 : delete [] yvec;
359 0 : delete [] xerr;
360 0 : delete [] yerr;
361 : return graphErrors;
362 0 : }
363 :
364 : TH2* AliTPCcalibBase::NormalizedProjection(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t normDim, Float_t minStatFrac)
365 : {
366 : // Create projection in axisDim1, axisDim2 in slices of normDim
367 : // Normalize the the entries in the dimension normDim and the sum up
368 : // Don't use bin for sum if statistics is less than minStatFrac of average number of entries
369 :
370 0 : if (!h) return 0x0;
371 : // ---| axis ranges of dimension in which to normalize |----------------------
372 0 : TAxis *axisNorm=h->GetAxis(normDim);
373 0 : const Int_t first=axisNorm->GetFirst();
374 0 : const Int_t last =axisNorm->GetLast();
375 0 : const Int_t nbins=last-first+1;
376 :
377 : // ---| final merged histogram |----------------------------------------------
378 0 : TH2 * hist = h->Projection(axisDim1, axisDim2);
379 0 : const Double_t statPerBin = hist->Integral()/Double_t(nbins);
380 0 : hist->Reset();
381 0 : hist->SetName(Form("%s_norm", hist->GetName()));
382 : // printf("statPerBin: %.2f\n", statPerBin);
383 : // printf("first, last, bins: %i, %i, %i\n", first, last, nbins);
384 :
385 : // ---| array with 2D projections per bin and statistics
386 0 : TObjArray arrProjections(nbins);
387 0 : arrProjections.SetOwner();
388 0 : Double_t *values=new Double_t[nbins];
389 :
390 : // ---| loop over normDim bins and do projections |---------------------------
391 : Int_t nused=0;
392 0 : for (Int_t ibin=first; ibin<=last; ++ibin) {
393 0 : axisNorm->SetRange(ibin,ibin);
394 0 : TH2 *proj=h->Projection(axisDim1, axisDim2);
395 0 : if (!proj) continue;
396 0 : proj->SetName(Form("%s_%d", h->GetName(), ibin));
397 : // proj->SetDirectory(0x0);
398 :
399 0 : proj->SetUniqueID(0); //steer nomalisation
400 :
401 0 : const Double_t stat=proj->Integral();
402 : // printf("bin: %i (%.2f)\n", ibin, stat);
403 0 : if (stat>minStatFrac*statPerBin) {
404 0 : proj->SetUniqueID(1);
405 0 : proj->Scale(1./stat);
406 0 : values[nused]=stat;
407 0 : ++nused;
408 : // printf(" used\n");
409 0 : }
410 0 : arrProjections.Add(proj);
411 0 : }
412 :
413 : // ---| Get Median of used bins, scale and sum up |---------------------------
414 0 : Double_t median=TMath::Median(nused, values);
415 : // printf("median: %.2f\n", median);
416 0 : for (Int_t ihist=0; ihist<nbins; ++ihist) {
417 0 : TH2 *proj=static_cast<TH2*>(arrProjections.At(ihist));
418 0 : if (!proj) continue;
419 0 : if (proj->GetUniqueID()) proj->Scale(median);
420 0 : hist->Add(proj);
421 0 : }
422 :
423 0 : axisNorm->SetRange(first,last);
424 0 : delete [] values;
425 :
426 : return hist;
427 0 : }
428 :
429 :
430 : void AliTPCcalibBase::BinLogX(THnSparse *h, Int_t axisDim) {
431 :
432 : // Method for the correct logarithmic binning of histograms
433 :
434 0 : TAxis *axis = h->GetAxis(axisDim);
435 0 : int bins = axis->GetNbins();
436 :
437 0 : Double_t from = axis->GetXmin();
438 0 : Double_t to = axis->GetXmax();
439 0 : Double_t *new_bins = new Double_t[bins + 1];
440 :
441 0 : new_bins[0] = from;
442 0 : Double_t factor = pow(to/from, 1./bins);
443 :
444 0 : for (int i = 1; i <= bins; i++) {
445 0 : new_bins[i] = factor * new_bins[i-1];
446 : }
447 0 : axis->Set(bins, new_bins);
448 0 : delete [] new_bins;
449 :
450 0 : }
451 : void AliTPCcalibBase::BinLogX(TH1 *h) {
452 :
453 : // Method for the correct logarithmic binning of histograms
454 :
455 0 : TAxis *axis = h->GetXaxis();
456 0 : int bins = axis->GetNbins();
457 :
458 0 : Double_t from = axis->GetXmin();
459 0 : Double_t to = axis->GetXmax();
460 0 : Double_t *new_bins = new Double_t[bins + 1];
461 :
462 0 : new_bins[0] = from;
463 0 : Double_t factor = pow(to/from, 1./bins);
464 :
465 0 : for (int i = 1; i <= bins; i++) {
466 0 : new_bins[i] = factor * new_bins[i-1];
467 : }
468 0 : axis->Set(bins, new_bins);
469 0 : delete [] new_bins;
470 :
471 0 : }
472 :
473 : void AliTPCcalibBase::BinLogX(TAxis *axis) {
474 :
475 : // Method for the correct logarithmic binning of histograms
476 :
477 0 : Int_t bins = axis->GetNbins();
478 :
479 0 : Double_t from = axis->GetXmin();
480 0 : Double_t to = axis->GetXmax();
481 0 : if (from<0) return;
482 0 : Double_t *new_bins = new Double_t[bins + 1];
483 :
484 0 : new_bins[0] = from;
485 0 : Double_t factor = pow(to/from, 1./bins);
486 :
487 0 : for (int i = 1; i <= bins; i++) {
488 0 : new_bins[i] = factor * new_bins[i-1];
489 : }
490 0 : axis->Set(bins, new_bins);
491 0 : delete [] new_bins;
492 0 : }
|
