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