]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibTimeGain.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
new module for hadron correlations
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTimeGain.cxx
CommitLineData
489dce1b 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
592a0c8f 18This class provides the calibration of the time dependence of the TPC gain due to pressure and temperature changes etc.
19
20
21//0. Libraries to load
489dce1b 22gSystem->Load("libANALYSIS");
23gSystem->Load("libSTAT");
24gSystem->Load("libTPCcalib");
25
26
27//1. Do calibration ...
28//
29// compare reference
30
31//
d0b31f57 32//2. Visualize results
489dce1b 33//
34TFile fcalib("CalibObjects.root");
35TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
36AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
3af3fbc4 37TGraphErrors * gr = gain->GetGraphGainVsTime(0,1000)
489dce1b 38
3af3fbc4 39// gain->GetHistGainTime()->GetAxis(1)->SetRangeUser(1213.8e6,1214.3e6)
592a0c8f 40TH2D * GainTime = gain->GetHistGainTime()->Projection(0,1)
41GainTime->GetXaxis()->SetTimeDisplay(kTRUE)
42GainTime->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}")
43GainTime->Draw("colz")
489dce1b 44
45//
46// MakeSlineFit example
47//
592a0c8f 48AliSplineFit *fit = AliTPCcalibTimeGain::MakeSplineFit(gr)
49
489dce1b 50TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
51
52gr->SetMarkerStyle(25);
592a0c8f 53gr->Draw("lp");
489dce1b 54grfit->SetLineColor(2);
55grfit->Draw("lu");
56
3af3fbc4 57//
58// QA - dE/dx resoultion as a function of time
59//TCa
60
61TGraph * grSigma = AliTPCcalibBase::FitSlicesSigma(gain->GetHistGainTime(),0,1,1800,5)
62
63TCanvas *c1 = new TCanvas("c1","transparent pad",200,10,700,500);
64 TPad *pad1 = new TPad("pad1","",0,0,1,1);
65 TPad *pad2 = new TPad("pad2","",0,0,1,1);
66 pad2->SetFillStyle(4000); //will be transparent
67 pad1->Draw();
68 pad1->cd();
69
70GainTime->Draw("colz")
71gr->Draw("lp")
72
73
74
75 c1->cd();
76 Double_t ymin = -0.04;
77 Double_t ymax = 0.12;
78Double_t dy = (ymax-ymin)/0.8;
79Double_t xmin = GainTime->GetXaxis()->GetXmin()
80Double_t xmax = GainTime->GetXaxis()->GetXmax()
81Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
82pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);
83 pad2->Draw();
84 pad2->cd();
85grSigma->SetLineColor(2);
86grSigma->SetLineWidth(2);
87grSigma->Draw("lp")
88TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
89 axis->SetLabelColor(kRed);
90 axis->SetTitle("dE/dx resolution #sigma_{dE/dx}");
91 axis->Draw();
92
93////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
94
95 ----> make Attachment study
96
97TFile fcalib("CalibObjects40366b.root");
98TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
99AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
100TGraphErrors * grAttach = gain->GetGraphAttachment(2000,4)
101
102TCanvas *c1 = new TCanvas("c1","transparent pad",200,10,700,500);
103 TPad *pad1 = new TPad("pad1","",0,0,1,1);
104 TPad *pad2 = new TPad("pad2","",0,0,1,1);
105 pad2->SetFillStyle(4000); //will be transparent
106 pad1->Draw();
107 pad1->cd();
108
109gain->GetHistGainTime()->GetAxis(1)->SetRangeUser(1213.8e6,1214.3e6)
110TH2D * GainTime = gain->GetHistGainTime()->Projection(0,1)
111GainTime->GetXaxis()->SetTimeDisplay(kTRUE)
112GainTime->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}")
113GainTime->Draw("colz")
114//gr->Draw("lp")
115
116 c1->cd();
117 Double_t ymin = -0.001;
118 Double_t ymax = 0.001;
119Double_t dy = (ymax-ymin)/0.8;
120Double_t xmin = GainTime->GetXaxis()->GetXmin()
121Double_t xmax = GainTime->GetXaxis()->GetXmax()
122Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
123pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);
124 pad2->Draw();
125 pad2->cd();
126grAttach->SetLineColor(2);
127grAttach->SetLineWidth(2);
128grAttach->Draw("lp")
129TGaxis *axis = new TGaxis(xmax,ymin,xmax,ymax,ymin,ymax,50510,"+L");
130 axis->SetLabelColor(kRed);
131 axis->SetTitle("attachment coefficient b");
132 axis->Draw();
133
134
489dce1b 135*/
136
137
138#include "Riostream.h"
139#include "TChain.h"
140#include "TTree.h"
141#include "TH1F.h"
142#include "TH2F.h"
143#include "TH3F.h"
144#include "THnSparse.h"
d0b31f57 145#include "TGraphErrors.h"
489dce1b 146#include "TList.h"
147#include "TMath.h"
148#include "TCanvas.h"
149#include "TFile.h"
150#include "TF1.h"
151#include "TVectorD.h"
152#include "TProfile.h"
153#include "TGraphErrors.h"
154#include "TCanvas.h"
155
156#include "AliTPCclusterMI.h"
157#include "AliTPCseed.h"
158#include "AliESDVertex.h"
159#include "AliESDEvent.h"
160#include "AliESDfriend.h"
161#include "AliESDInputHandler.h"
162#include "AliAnalysisManager.h"
163
164#include "AliTracker.h"
165#include "AliMagF.h"
166#include "AliTPCCalROC.h"
167
168#include "AliLog.h"
169
170#include "AliTPCcalibTimeGain.h"
171
172#include "TTreeStream.h"
173#include "AliTPCTracklet.h"
174#include "TTimeStamp.h"
175#include "AliTPCcalibDB.h"
176#include "AliTPCcalibLaser.h"
177#include "AliDCSSensorArray.h"
178#include "AliDCSSensor.h"
179
180ClassImp(AliTPCcalibTimeGain)
181
182
183AliTPCcalibTimeGain::AliTPCcalibTimeGain()
184 :AliTPCcalibBase(),
185 fHistGainTime(0),
186 fGainVsTime(0),
187 fHistDeDxTotal(0),
188 fIntegrationTimeDeDx(0),
189 fMIP(0),
592a0c8f 190 fUseMax(0),
489dce1b 191 fLowerTrunc(0),
192 fUpperTrunc(0),
193 fUseShapeNorm(0),
194 fUsePosNorm(0),
195 fUsePadNorm(0),
592a0c8f 196 fUseCookAnalytical(0),
d0b31f57 197 fIsCosmic(0),
198 fLowMemoryConsumption(0)
489dce1b 199{
200 AliInfo("Default Constructor");
201}
202
203
204AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
205 :AliTPCcalibBase(),
206 fHistGainTime(0),
207 fGainVsTime(0),
208 fHistDeDxTotal(0),
209 fIntegrationTimeDeDx(0),
210 fMIP(0),
592a0c8f 211 fUseMax(0),
489dce1b 212 fLowerTrunc(0),
213 fUpperTrunc(0),
214 fUseShapeNorm(0),
215 fUsePosNorm(0),
216 fUsePadNorm(0),
592a0c8f 217 fUseCookAnalytical(0),
d0b31f57 218 fIsCosmic(0),
219 fLowMemoryConsumption(0)
3af3fbc4 220{
489dce1b 221
222 SetName(name);
223 SetTitle(title);
3af3fbc4 224
489dce1b 225 AliInfo("Non Default Constructor");
3af3fbc4 226
489dce1b 227 fIntegrationTimeDeDx = deltaIntegrationTimeGain;
3af3fbc4 228
489dce1b 229 Double_t deltaTime = EndTime - StartTime;
230
3af3fbc4 231
592a0c8f 232 // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available), run number
d0b31f57 233 Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
3af3fbc4 234 Int_t binsGainTime[6] = {150, timeBins, 2, 25, 200, 10000000};
592a0c8f 235 Double_t xminGainTime[6] = {0.5, StartTime, 0.5, 0, 0.1, -0.5};
3af3fbc4 236 Double_t xmaxGainTime[6] = { 8, EndTime, 2.5, 250, 50, 9999999.5};
592a0c8f 237 fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number;dEdx",6,binsGainTime,xminGainTime,xmaxGainTime);
489dce1b 238 BinLogX(fHistGainTime, 4);
239 //
3af3fbc4 240 fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,8);
489dce1b 241 BinLogX(fHistDeDxTotal);
242
243 // default values for dE/dx
244 fMIP = 50.;
592a0c8f 245 fUseMax = kTRUE;
489dce1b 246 fLowerTrunc = 0.0;
247 fUpperTrunc = 0.7;
248 fUseShapeNorm = kTRUE;
249 fUsePosNorm = kFALSE;
250 fUsePadNorm = kFALSE;
592a0c8f 251 fUseCookAnalytical = kFALSE;
489dce1b 252 //
253 fIsCosmic = kTRUE;
3af3fbc4 254 fLowMemoryConsumption = kTRUE;
489dce1b 255 //
256
3af3fbc4 257}
489dce1b 258
259
260
261AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
262 //
263 //
264 //
2cd9f718 265 delete fHistGainTime;
266 delete fGainVsTime;
267 delete fHistDeDxTotal;
489dce1b 268}
269
270
271void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
272 //
273 // main track loop
274 //
275 if (!event) {
276 Printf("ERROR: ESD not available");
277 return;
d0b31f57 278 }
279
280 if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
281 ProcessCosmicEvent(event);
282 } else {
283 ProcessBeamEvent(event);
284 }
285
489dce1b 286
489dce1b 287
d0b31f57 288
289}
290
291
292void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
293
489dce1b 294 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
295 if (!ESDfriend) {
296 Printf("ERROR: ESDfriend not available");
297 return;
298 }
299 //
300 UInt_t time = event->GetTimeStamp();
d0b31f57 301 Int_t ntracks = event->GetNumberOfTracks();
592a0c8f 302 Int_t runNumber = event->GetRunNumber();
489dce1b 303 //
304 // track loop
305 //
306 for (Int_t i=0;i<ntracks;++i) {
307
308 AliESDtrack *track = event->GetTrack(i);
309 if (!track) continue;
310
311 const AliExternalTrackParam * trackIn = track->GetInnerParam();
312 const AliExternalTrackParam * trackOut = track->GetOuterParam();
313 if (!trackIn) continue;
314 if (!trackOut) continue;
315
316 // calculate necessary track parameters
489dce1b 317 Double_t meanP = trackIn->GetP();
318 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
489dce1b 319 Int_t NclsDeDx = track->GetTPCNcls();
320
489dce1b 321 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
d0b31f57 322 if (NclsDeDx < 60) continue;
489dce1b 323 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
324 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
325
326 // Get seeds
327 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
328 if (!friendTrack) continue;
329 TObject *calibObject;
330 AliTPCseed *seed = 0;
331 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
332 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
333 }
334
335 if (seed) {
592a0c8f 336 Double_t TPCsignal = GetTPCdEdx(seed);
3af3fbc4 337 if (NclsDeDx > 100) fHistDeDxTotal->Fill(meanP, TPCsignal);
489dce1b 338 //
d0b31f57 339 if (fLowMemoryConsumption) {
340 if (meanP < 20) continue;
341 meanP = 30; // set all momenta to one in order to save memory
489dce1b 342 }
d0b31f57 343 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
592a0c8f 344 Double_t vec[6] = {TPCsignal,time,1,meanDrift,meanP,runNumber};
d0b31f57 345 fHistGainTime->Fill(vec);
346 }
489dce1b 347
d0b31f57 348 }
489dce1b 349
d0b31f57 350}
351
352
353
354void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
355
356 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
357 if (!ESDfriend) {
358 Printf("ERROR: ESDfriend not available");
359 return;
360 }
361 //
362 UInt_t time = event->GetTimeStamp();
363 Int_t ntracks = event->GetNumberOfTracks();
592a0c8f 364 Int_t runNumber = event->GetRunNumber();
d0b31f57 365 //
366 // track loop
367 //
368 for (Int_t i=0;i<ntracks;++i) {
369
370 AliESDtrack *track = event->GetTrack(i);
371 if (!track) continue;
372
373 const AliExternalTrackParam * trackIn = track->GetInnerParam();
374 const AliExternalTrackParam * trackOut = track->GetOuterParam();
375 if (!trackIn) continue;
376 if (!trackOut) continue;
377
378 // calculate necessary track parameters
379 Double_t meanP = trackIn->GetP();
380 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
381 Int_t NclsDeDx = track->GetTPCNcls();
382
383 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
384 if (NclsDeDx < 60) continue;
385 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
386 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
387
388 // Get seeds
389 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
390 if (!friendTrack) continue;
391 TObject *calibObject;
392 AliTPCseed *seed = 0;
393 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
394 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
395 }
396
397 if (seed) {
592a0c8f 398 Double_t TPCsignal = GetTPCdEdx(seed);
399 fHistDeDxTotal->Fill(meanP, TPCsignal);
d0b31f57 400 //
401 if (fLowMemoryConsumption) {
402 if (meanP > 0.5 || meanP < 0.4) continue;
403 meanP = 0.45; // set all momenta to one in order to save memory
404 }
405 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
592a0c8f 406 Double_t vec[6] = {TPCsignal,time,2,meanDrift,meanP,runNumber};
d0b31f57 407 fHistGainTime->Fill(vec);
489dce1b 408 }
409
410 }
d0b31f57 411
489dce1b 412}
413
414
592a0c8f 415Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
416
417 Double_t signal = 0;
418 //
419 if (!fUseCookAnalytical) {
420 signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
421 } else {
422 signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
423 }
424 //
425 return signal;
426}
427
428
429void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
489dce1b 430 //
431 //
432 //
489dce1b 433 if (fIsCosmic) {
d0b31f57 434 fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
435 fHistGainTime->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons
489dce1b 436 } else {
d0b31f57 437 fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
438 fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
489dce1b 439 }
489dce1b 440 //
592a0c8f 441 fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,minEntries,10);
489dce1b 442 //
443 return;
444}
445
446
592a0c8f 447TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
448 //
449 //
450 //
451 if (runNumber == 0) {
452 if (!fGainVsTime) {
453 AnalyzeRun(minEntries);
454 }
455 } else {
456 // 1st check if the current run was cosmic or beam event
457 fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
458 AnalyzeRun(minEntries);
459 }
460 if (fGainVsTime->GetN() == 0) return 0;
461 return fGainVsTime;
462}
463
489dce1b 464Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
465
466 TIterator* iter = li->MakeIterator();
467 AliTPCcalibTimeGain* cal = 0;
468
469 while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
470 if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
471 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
472 return -1;
473 }
474
475 // add histograms here...
476 if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
477 if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
478
479 }
480
481 return 0;
482
483}
484
485
592a0c8f 486AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
487 //
488 //
489 //
490 AliSplineFit *fit = new AliSplineFit();
491 fit->SetGraph(graph);
492 fit->SetMinPoints(graph->GetN()+1);
493 fit->InitKnots(graph,2,0,0.001);
494 fit->SplineFit(0);
495 return fit;
496
497}
498
489dce1b 499
3af3fbc4 500
e55e5512 501TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t /*fracLow*/, Float_t /*fracUp*/) {
3af3fbc4 502 //
503 // For each time bin the driftlength-dependence of the signal is fitted.
504 //
505 TH3D * hist = fHistGainTime->Projection(1, 0, 3);
506 Double_t *xvec = new Double_t[hist->GetNbinsX()];
507 Double_t *yvec = new Double_t[hist->GetNbinsX()];
508 Double_t *xerr = new Double_t[hist->GetNbinsX()];
509 Double_t *yerr = new Double_t[hist->GetNbinsX()];
510 Int_t counter = 0;
511 TH2D * projectionHist = 0x0;
512 //
513 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
514 Int_t nsum=0;
515 Int_t imin = i;
516 Int_t imax = i;
517 for (Int_t idelta=0; idelta<nmaxBin; idelta++){
518 //
519 imin = TMath::Max(i-idelta,1);
520 imax = TMath::Min(i+idelta,hist->GetNbinsX());
521 nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()-1,1,hist->GetNbinsZ()-1));
522 if (nsum==0) break;
523 if (nsum>minEntries) break;
524 }
525 if (nsum<minEntries) continue;
526 //
527 hist->GetXaxis()->SetRange(imin,imax);
528 projectionHist = (TH2D*)hist->Project3D("yzNUFNOF");
529 //
530 TObjArray arr;
531 projectionHist->FitSlicesY(0,2, projectionHist->GetNbinsX()-2,0,"QNR",&arr);
532 TH1D * histAttach = (TH1D*)arr.At(1);
533 TF1 pol1("polynom1","pol1",10,240);
534 histAttach->Fit(&pol1,"QNR");
535 xvec[counter] = 0.5*(hist->GetXaxis()->GetBinCenter(imin)+hist->GetXaxis()->GetBinCenter(imax));
536 yvec[counter] = pol1.GetParameter(1)/pol1.GetParameter(0);
537 xerr[counter] = 0;
538 yerr[counter] = pol1.GetParError(1)/pol1.GetParameter(0);
539 counter++;
540 //
541 delete projectionHist;
542 }
543
544 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
545 delete [] xvec;
546 delete [] yvec;
547 delete [] xerr;
548 delete [] yerr;
549 delete hist;
550 return graphErrors;
551
552}
553
554
555
489dce1b 556void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
557
558 // Method for the correct logarithmic binning of histograms
559
560 TAxis *axis = h->GetAxis(axisDim);
561 int bins = axis->GetNbins();
562
563 Double_t from = axis->GetXmin();
564 Double_t to = axis->GetXmax();
565 Double_t *new_bins = new Double_t[bins + 1];
566
567 new_bins[0] = from;
568 Double_t factor = pow(to/from, 1./bins);
569
570 for (int i = 1; i <= bins; i++) {
571 new_bins[i] = factor * new_bins[i-1];
572 }
573 axis->Set(bins, new_bins);
574 delete new_bins;
575
576}
577
578
579void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
580
581 // Method for the correct logarithmic binning of histograms
582
583 TAxis *axis = h->GetXaxis();
584 int bins = axis->GetNbins();
585
586 Double_t from = axis->GetXmin();
587 Double_t to = axis->GetXmax();
588 Double_t *new_bins = new Double_t[bins + 1];
589
590 new_bins[0] = from;
591 Double_t factor = pow(to/from, 1./bins);
592
593 for (int i = 1; i <= bins; i++) {
594 new_bins[i] = factor * new_bins[i-1];
595 }
596 axis->Set(bins, new_bins);
597 delete new_bins;
598
599}
600
601
602
603void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
604 //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
605 const Double_t sigma = 0.06;
606
d0b31f57 607 TH2F *histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",TMath::Nint(0.5*hist->GetNbinsX()),0.1,5000.,hist->GetNbinsY(),hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
489dce1b 608 BinLogX(histBG);
609
d0b31f57 610 TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
489dce1b 611 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
612 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
613 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
614 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
615 foElectron->SetParameters(ini);
616 foMuon->SetParameters(ini);
617 foPion->SetParameters(ini);
618 foKaon->SetParameters(ini);
619 foProton->SetParameters(ini);
620
621 TCanvas *CanvCheck1 = new TCanvas();
622 hist->Draw("colz");
623 foElectron->Draw("same");
624 foMuon->Draw("same");
625 foPion->Draw("same");
626 foKaon->Draw("same");
627 foProton->Draw("same");
628
629 // Loop over all points of the input histogram
630
631 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
632 Double_t x = hist->GetXaxis()->GetBinCenter(i);
633 for(Int_t j=1; j < hist->GetNbinsY(); j++) {
634 Long64_t n = hist->GetBin(i, j);
635 Double_t y = hist->GetYaxis()->GetBinCenter(j);
636 Double_t entries = hist->GetBinContent(n);
637 Double_t mass = 0;
638
639 // 1. identify protons
640 mass = 0.938;
641 if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
642 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
643 }
644
645 // 2. identify electrons
646 mass = 0.000511;
647 if (fIsCosmic) {
648 if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
649 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
650 }
651 } else {
d0b31f57 652 if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 2*sigma && ((x>0.25&&x<0.35) || (x>1.5&&x<1.8) || (x>0.65&&x<0.7))) {
489dce1b 653 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
654 }
655 }
656
657 // 3. identify either muons or pions depending on cosmic or not
658 if (fIsCosmic) {
659 mass = 0.1056;
660 if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
661 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
662 }
663 } else {
664 mass = 0.1396;
665 if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
666 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
667 }
668 }
669
670 // 4. for pp also Kaons must be included
671 if (!fIsCosmic) {
672 mass = 0.4936;
673 if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
674 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
675 }
676 }
677 }
678 }
679
680 // Fit Aleph-Parameters to the obtained profile
681 TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
682 funcAlephD->SetParameters(ini);
683
684 TCanvas *CanvCheck2 = new TCanvas();
685 histBG->Draw();
686
687 //FitSlices
688 TObjArray * arr = new TObjArray();
689 histBG->FitSlicesY(0,0,-1,0,"QN",arr);
690 TH1D * fitMean = (TH1D*) arr->At(1);
691 fitMean->Draw("same");
692
693 funcAlephD->SetParLimits(2,1e-3,1e-7);
694 funcAlephD->SetParLimits(3,0.5,3.5);
695 funcAlephD->SetParLimits(4,0.5,3.5);
696 fitMean->Fit(funcAlephD, "QNR");
697 funcAlephD->Draw("same");
698
699 for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
700
701 foElectron->SetParameters(ini);
702 foMuon->SetParameters(ini);
703 foPion->SetParameters(ini);
704 foKaon->SetParameters(ini);
705 foProton->SetParameters(ini);
706
707 TCanvas *CanvCheck3 = new TCanvas();
708 hist->Draw("colz");
709 foElectron->Draw("same");
710 foMuon->Draw("same");
711 foPion->Draw("same");
712 foKaon->Draw("same");
713 foProton->Draw("same");
714
715 CanvCheck1->Print();
716 CanvCheck2->Print();
717 CanvCheck3->Print();
718
719 return;
720
721
722}