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