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