]> 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
set defaults for PbPb at LHC
[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 //
265}
266
267
268void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
269 //
270 // main track loop
271 //
272 if (!event) {
273 Printf("ERROR: ESD not available");
274 return;
d0b31f57 275 }
276
277 if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
278 ProcessCosmicEvent(event);
279 } else {
280 ProcessBeamEvent(event);
281 }
282
489dce1b 283
489dce1b 284
d0b31f57 285
286}
287
288
289void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
290
489dce1b 291 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
292 if (!ESDfriend) {
293 Printf("ERROR: ESDfriend not available");
294 return;
295 }
296 //
297 UInt_t time = event->GetTimeStamp();
d0b31f57 298 Int_t ntracks = event->GetNumberOfTracks();
592a0c8f 299 Int_t runNumber = event->GetRunNumber();
489dce1b 300 //
301 // track loop
302 //
303 for (Int_t i=0;i<ntracks;++i) {
304
305 AliESDtrack *track = event->GetTrack(i);
306 if (!track) continue;
307
308 const AliExternalTrackParam * trackIn = track->GetInnerParam();
309 const AliExternalTrackParam * trackOut = track->GetOuterParam();
310 if (!trackIn) continue;
311 if (!trackOut) continue;
312
313 // calculate necessary track parameters
489dce1b 314 Double_t meanP = trackIn->GetP();
315 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
489dce1b 316 Int_t NclsDeDx = track->GetTPCNcls();
317
489dce1b 318 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
d0b31f57 319 if (NclsDeDx < 60) continue;
489dce1b 320 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
321 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
322
323 // Get seeds
324 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
325 if (!friendTrack) continue;
326 TObject *calibObject;
327 AliTPCseed *seed = 0;
328 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
329 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
330 }
331
332 if (seed) {
592a0c8f 333 Double_t TPCsignal = GetTPCdEdx(seed);
3af3fbc4 334 if (NclsDeDx > 100) fHistDeDxTotal->Fill(meanP, TPCsignal);
489dce1b 335 //
d0b31f57 336 if (fLowMemoryConsumption) {
337 if (meanP < 20) continue;
338 meanP = 30; // set all momenta to one in order to save memory
489dce1b 339 }
d0b31f57 340 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
592a0c8f 341 Double_t vec[6] = {TPCsignal,time,1,meanDrift,meanP,runNumber};
d0b31f57 342 fHistGainTime->Fill(vec);
343 }
489dce1b 344
d0b31f57 345 }
489dce1b 346
d0b31f57 347}
348
349
350
351void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
352
353 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
354 if (!ESDfriend) {
355 Printf("ERROR: ESDfriend not available");
356 return;
357 }
358 //
359 UInt_t time = event->GetTimeStamp();
360 Int_t ntracks = event->GetNumberOfTracks();
592a0c8f 361 Int_t runNumber = event->GetRunNumber();
d0b31f57 362 //
363 // track loop
364 //
365 for (Int_t i=0;i<ntracks;++i) {
366
367 AliESDtrack *track = event->GetTrack(i);
368 if (!track) continue;
369
370 const AliExternalTrackParam * trackIn = track->GetInnerParam();
371 const AliExternalTrackParam * trackOut = track->GetOuterParam();
372 if (!trackIn) continue;
373 if (!trackOut) continue;
374
375 // calculate necessary track parameters
376 Double_t meanP = trackIn->GetP();
377 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
378 Int_t NclsDeDx = track->GetTPCNcls();
379
380 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
381 if (NclsDeDx < 60) continue;
382 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
383 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
384
385 // Get seeds
386 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
387 if (!friendTrack) continue;
388 TObject *calibObject;
389 AliTPCseed *seed = 0;
390 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
391 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
392 }
393
394 if (seed) {
592a0c8f 395 Double_t TPCsignal = GetTPCdEdx(seed);
396 fHistDeDxTotal->Fill(meanP, TPCsignal);
d0b31f57 397 //
398 if (fLowMemoryConsumption) {
399 if (meanP > 0.5 || meanP < 0.4) continue;
400 meanP = 0.45; // set all momenta to one in order to save memory
401 }
402 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
592a0c8f 403 Double_t vec[6] = {TPCsignal,time,2,meanDrift,meanP,runNumber};
d0b31f57 404 fHistGainTime->Fill(vec);
489dce1b 405 }
406
407 }
d0b31f57 408
489dce1b 409}
410
411
592a0c8f 412Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
413
414 Double_t signal = 0;
415 //
416 if (!fUseCookAnalytical) {
417 signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
418 } else {
419 signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
420 }
421 //
422 return signal;
423}
424
425
426void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
489dce1b 427 //
428 //
429 //
489dce1b 430 if (fIsCosmic) {
d0b31f57 431 fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
432 fHistGainTime->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons
489dce1b 433 } else {
d0b31f57 434 fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
435 fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
489dce1b 436 }
489dce1b 437 //
592a0c8f 438 fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,minEntries,10);
489dce1b 439 //
440 return;
441}
442
443
592a0c8f 444TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
445 //
446 //
447 //
448 if (runNumber == 0) {
449 if (!fGainVsTime) {
450 AnalyzeRun(minEntries);
451 }
452 } else {
453 // 1st check if the current run was cosmic or beam event
454 fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
455 AnalyzeRun(minEntries);
456 }
457 if (fGainVsTime->GetN() == 0) return 0;
458 return fGainVsTime;
459}
460
489dce1b 461Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
462
463 TIterator* iter = li->MakeIterator();
464 AliTPCcalibTimeGain* cal = 0;
465
466 while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
467 if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
468 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
469 return -1;
470 }
471
472 // add histograms here...
473 if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
474 if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
475
476 }
477
478 return 0;
479
480}
481
482
592a0c8f 483AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
484 //
485 //
486 //
487 AliSplineFit *fit = new AliSplineFit();
488 fit->SetGraph(graph);
489 fit->SetMinPoints(graph->GetN()+1);
490 fit->InitKnots(graph,2,0,0.001);
491 fit->SplineFit(0);
492 return fit;
493
494}
495
489dce1b 496
3af3fbc4 497
e55e5512 498TGraphErrors * AliTPCcalibTimeGain::GetGraphAttachment(Int_t minEntries, Int_t nmaxBin, Float_t /*fracLow*/, Float_t /*fracUp*/) {
3af3fbc4 499 //
500 // For each time bin the driftlength-dependence of the signal is fitted.
501 //
502 TH3D * hist = fHistGainTime->Projection(1, 0, 3);
503 Double_t *xvec = new Double_t[hist->GetNbinsX()];
504 Double_t *yvec = new Double_t[hist->GetNbinsX()];
505 Double_t *xerr = new Double_t[hist->GetNbinsX()];
506 Double_t *yerr = new Double_t[hist->GetNbinsX()];
507 Int_t counter = 0;
508 TH2D * projectionHist = 0x0;
509 //
510 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
511 Int_t nsum=0;
512 Int_t imin = i;
513 Int_t imax = i;
514 for (Int_t idelta=0; idelta<nmaxBin; idelta++){
515 //
516 imin = TMath::Max(i-idelta,1);
517 imax = TMath::Min(i+idelta,hist->GetNbinsX());
518 nsum = TMath::Nint(hist->Integral(imin,imax,1,hist->GetNbinsY()-1,1,hist->GetNbinsZ()-1));
519 if (nsum==0) break;
520 if (nsum>minEntries) break;
521 }
522 if (nsum<minEntries) continue;
523 //
524 hist->GetXaxis()->SetRange(imin,imax);
525 projectionHist = (TH2D*)hist->Project3D("yzNUFNOF");
526 //
527 TObjArray arr;
528 projectionHist->FitSlicesY(0,2, projectionHist->GetNbinsX()-2,0,"QNR",&arr);
529 TH1D * histAttach = (TH1D*)arr.At(1);
530 TF1 pol1("polynom1","pol1",10,240);
531 histAttach->Fit(&pol1,"QNR");
532 xvec[counter] = 0.5*(hist->GetXaxis()->GetBinCenter(imin)+hist->GetXaxis()->GetBinCenter(imax));
533 yvec[counter] = pol1.GetParameter(1)/pol1.GetParameter(0);
534 xerr[counter] = 0;
535 yerr[counter] = pol1.GetParError(1)/pol1.GetParameter(0);
536 counter++;
537 //
538 delete projectionHist;
539 }
540
541 TGraphErrors * graphErrors = new TGraphErrors(counter, xvec, yvec, xerr, yerr);
542 delete [] xvec;
543 delete [] yvec;
544 delete [] xerr;
545 delete [] yerr;
546 delete hist;
547 return graphErrors;
548
549}
550
551
552
489dce1b 553void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
554
555 // Method for the correct logarithmic binning of histograms
556
557 TAxis *axis = h->GetAxis(axisDim);
558 int bins = axis->GetNbins();
559
560 Double_t from = axis->GetXmin();
561 Double_t to = axis->GetXmax();
562 Double_t *new_bins = new Double_t[bins + 1];
563
564 new_bins[0] = from;
565 Double_t factor = pow(to/from, 1./bins);
566
567 for (int i = 1; i <= bins; i++) {
568 new_bins[i] = factor * new_bins[i-1];
569 }
570 axis->Set(bins, new_bins);
571 delete new_bins;
572
573}
574
575
576void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
577
578 // Method for the correct logarithmic binning of histograms
579
580 TAxis *axis = h->GetXaxis();
581 int bins = axis->GetNbins();
582
583 Double_t from = axis->GetXmin();
584 Double_t to = axis->GetXmax();
585 Double_t *new_bins = new Double_t[bins + 1];
586
587 new_bins[0] = from;
588 Double_t factor = pow(to/from, 1./bins);
589
590 for (int i = 1; i <= bins; i++) {
591 new_bins[i] = factor * new_bins[i-1];
592 }
593 axis->Set(bins, new_bins);
594 delete new_bins;
595
596}
597
598
599
600void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
601 //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
602 const Double_t sigma = 0.06;
603
d0b31f57 604 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 605 BinLogX(histBG);
606
d0b31f57 607 TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
489dce1b 608 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
609 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
610 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
611 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
612 foElectron->SetParameters(ini);
613 foMuon->SetParameters(ini);
614 foPion->SetParameters(ini);
615 foKaon->SetParameters(ini);
616 foProton->SetParameters(ini);
617
618 TCanvas *CanvCheck1 = new TCanvas();
619 hist->Draw("colz");
620 foElectron->Draw("same");
621 foMuon->Draw("same");
622 foPion->Draw("same");
623 foKaon->Draw("same");
624 foProton->Draw("same");
625
626 // Loop over all points of the input histogram
627
628 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
629 Double_t x = hist->GetXaxis()->GetBinCenter(i);
630 for(Int_t j=1; j < hist->GetNbinsY(); j++) {
631 Long64_t n = hist->GetBin(i, j);
632 Double_t y = hist->GetYaxis()->GetBinCenter(j);
633 Double_t entries = hist->GetBinContent(n);
634 Double_t mass = 0;
635
636 // 1. identify protons
637 mass = 0.938;
638 if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
639 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
640 }
641
642 // 2. identify electrons
643 mass = 0.000511;
644 if (fIsCosmic) {
645 if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
646 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
647 }
648 } else {
d0b31f57 649 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 650 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
651 }
652 }
653
654 // 3. identify either muons or pions depending on cosmic or not
655 if (fIsCosmic) {
656 mass = 0.1056;
657 if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
658 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
659 }
660 } else {
661 mass = 0.1396;
662 if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
663 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
664 }
665 }
666
667 // 4. for pp also Kaons must be included
668 if (!fIsCosmic) {
669 mass = 0.4936;
670 if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
671 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
672 }
673 }
674 }
675 }
676
677 // Fit Aleph-Parameters to the obtained profile
678 TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
679 funcAlephD->SetParameters(ini);
680
681 TCanvas *CanvCheck2 = new TCanvas();
682 histBG->Draw();
683
684 //FitSlices
685 TObjArray * arr = new TObjArray();
686 histBG->FitSlicesY(0,0,-1,0,"QN",arr);
687 TH1D * fitMean = (TH1D*) arr->At(1);
688 fitMean->Draw("same");
689
690 funcAlephD->SetParLimits(2,1e-3,1e-7);
691 funcAlephD->SetParLimits(3,0.5,3.5);
692 funcAlephD->SetParLimits(4,0.5,3.5);
693 fitMean->Fit(funcAlephD, "QNR");
694 funcAlephD->Draw("same");
695
696 for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
697
698 foElectron->SetParameters(ini);
699 foMuon->SetParameters(ini);
700 foPion->SetParameters(ini);
701 foKaon->SetParameters(ini);
702 foProton->SetParameters(ini);
703
704 TCanvas *CanvCheck3 = new TCanvas();
705 hist->Draw("colz");
706 foElectron->Draw("same");
707 foMuon->Draw("same");
708 foPion->Draw("same");
709 foKaon->Draw("same");
710 foProton->Draw("same");
711
712 CanvCheck1->Print();
713 CanvCheck2->Print();
714 CanvCheck3->Print();
715
716 return;
717
718
719}