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