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