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