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