The script to add the FMD analysis to the train
[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");
592a0c8f 37TGraphErrors * gr = gain->GetGraphGainVsTime(0,500)
489dce1b 38
592a0c8f 39TH2D * GainTime = gain->GetHistGainTime()->Projection(0,1)
40GainTime->GetXaxis()->SetTimeDisplay(kTRUE)
41GainTime->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}")
42GainTime->Draw("colz")
489dce1b 43
44//
45// MakeSlineFit example
46//
592a0c8f 47AliSplineFit *fit = AliTPCcalibTimeGain::MakeSplineFit(gr)
48
489dce1b 49TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
50
51gr->SetMarkerStyle(25);
592a0c8f 52gr->Draw("lp");
489dce1b 53grfit->SetLineColor(2);
54grfit->Draw("lu");
55
56*/
57
58
59#include "Riostream.h"
60#include "TChain.h"
61#include "TTree.h"
62#include "TH1F.h"
63#include "TH2F.h"
64#include "TH3F.h"
65#include "THnSparse.h"
d0b31f57 66#include "TGraphErrors.h"
489dce1b 67#include "TList.h"
68#include "TMath.h"
69#include "TCanvas.h"
70#include "TFile.h"
71#include "TF1.h"
72#include "TVectorD.h"
73#include "TProfile.h"
74#include "TGraphErrors.h"
75#include "TCanvas.h"
76
77#include "AliTPCclusterMI.h"
78#include "AliTPCseed.h"
79#include "AliESDVertex.h"
80#include "AliESDEvent.h"
81#include "AliESDfriend.h"
82#include "AliESDInputHandler.h"
83#include "AliAnalysisManager.h"
84
85#include "AliTracker.h"
86#include "AliMagF.h"
87#include "AliTPCCalROC.h"
88
89#include "AliLog.h"
90
91#include "AliTPCcalibTimeGain.h"
92
93#include "TTreeStream.h"
94#include "AliTPCTracklet.h"
95#include "TTimeStamp.h"
96#include "AliTPCcalibDB.h"
97#include "AliTPCcalibLaser.h"
98#include "AliDCSSensorArray.h"
99#include "AliDCSSensor.h"
100
101ClassImp(AliTPCcalibTimeGain)
102
103
104AliTPCcalibTimeGain::AliTPCcalibTimeGain()
105 :AliTPCcalibBase(),
106 fHistGainTime(0),
107 fGainVsTime(0),
108 fHistDeDxTotal(0),
109 fIntegrationTimeDeDx(0),
110 fMIP(0),
592a0c8f 111 fUseMax(0),
489dce1b 112 fLowerTrunc(0),
113 fUpperTrunc(0),
114 fUseShapeNorm(0),
115 fUsePosNorm(0),
116 fUsePadNorm(0),
592a0c8f 117 fUseCookAnalytical(0),
d0b31f57 118 fIsCosmic(0),
119 fLowMemoryConsumption(0)
489dce1b 120{
121 AliInfo("Default Constructor");
122}
123
124
125AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
126 :AliTPCcalibBase(),
127 fHistGainTime(0),
128 fGainVsTime(0),
129 fHistDeDxTotal(0),
130 fIntegrationTimeDeDx(0),
131 fMIP(0),
592a0c8f 132 fUseMax(0),
489dce1b 133 fLowerTrunc(0),
134 fUpperTrunc(0),
135 fUseShapeNorm(0),
136 fUsePosNorm(0),
137 fUsePadNorm(0),
592a0c8f 138 fUseCookAnalytical(0),
d0b31f57 139 fIsCosmic(0),
140 fLowMemoryConsumption(0)
489dce1b 141 {
142
143 SetName(name);
144 SetTitle(title);
145
146 AliInfo("Non Default Constructor");
147
148 fIntegrationTimeDeDx = deltaIntegrationTimeGain;
149
150 Double_t deltaTime = EndTime - StartTime;
151
152
592a0c8f 153 // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available), run number
d0b31f57 154 Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
592a0c8f 155 Int_t binsGainTime[6] = {100, timeBins, 2, 25, 200, 10000000};
156 Double_t xminGainTime[6] = {0.5, StartTime, 0.5, 0, 0.1, -0.5};
157 Double_t xmaxGainTime[6] = { 4, EndTime, 2.5, 250, 50, 9999999.5};
158 fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx,time,type,driftlength,momenta,run number;dEdx",6,binsGainTime,xminGainTime,xmaxGainTime);
489dce1b 159 BinLogX(fHistGainTime, 4);
160 //
161 fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000);
162 BinLogX(fHistDeDxTotal);
163
164 // default values for dE/dx
165 fMIP = 50.;
592a0c8f 166 fUseMax = kTRUE;
489dce1b 167 fLowerTrunc = 0.0;
168 fUpperTrunc = 0.7;
169 fUseShapeNorm = kTRUE;
170 fUsePosNorm = kFALSE;
171 fUsePadNorm = kFALSE;
592a0c8f 172 fUseCookAnalytical = kFALSE;
489dce1b 173 //
174 fIsCosmic = kTRUE;
592a0c8f 175 fLowMemoryConsumption = kFALSE;
489dce1b 176 //
177
178 }
179
180
181
182AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
183 //
184 //
185 //
186}
187
188
189void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
190 //
191 // main track loop
192 //
193 if (!event) {
194 Printf("ERROR: ESD not available");
195 return;
d0b31f57 196 }
197
198 if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
199 ProcessCosmicEvent(event);
200 } else {
201 ProcessBeamEvent(event);
202 }
203
489dce1b 204
489dce1b 205
d0b31f57 206
207}
208
209
210void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
211
489dce1b 212 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
213 if (!ESDfriend) {
214 Printf("ERROR: ESDfriend not available");
215 return;
216 }
217 //
218 UInt_t time = event->GetTimeStamp();
d0b31f57 219 Int_t ntracks = event->GetNumberOfTracks();
592a0c8f 220 Int_t runNumber = event->GetRunNumber();
489dce1b 221 //
222 // track loop
223 //
224 for (Int_t i=0;i<ntracks;++i) {
225
226 AliESDtrack *track = event->GetTrack(i);
227 if (!track) continue;
228
229 const AliExternalTrackParam * trackIn = track->GetInnerParam();
230 const AliExternalTrackParam * trackOut = track->GetOuterParam();
231 if (!trackIn) continue;
232 if (!trackOut) continue;
233
234 // calculate necessary track parameters
489dce1b 235 Double_t meanP = trackIn->GetP();
236 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
489dce1b 237 Int_t NclsDeDx = track->GetTPCNcls();
238
489dce1b 239 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
d0b31f57 240 if (NclsDeDx < 60) continue;
489dce1b 241 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
242 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
243
244 // Get seeds
245 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
246 if (!friendTrack) continue;
247 TObject *calibObject;
248 AliTPCseed *seed = 0;
249 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
250 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
251 }
252
253 if (seed) {
592a0c8f 254 Double_t TPCsignal = GetTPCdEdx(seed);
255 fHistDeDxTotal->Fill(meanP, TPCsignal);
489dce1b 256 //
d0b31f57 257 if (fLowMemoryConsumption) {
258 if (meanP < 20) continue;
259 meanP = 30; // set all momenta to one in order to save memory
489dce1b 260 }
d0b31f57 261 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
592a0c8f 262 Double_t vec[6] = {TPCsignal,time,1,meanDrift,meanP,runNumber};
d0b31f57 263 fHistGainTime->Fill(vec);
264 }
489dce1b 265
d0b31f57 266 }
489dce1b 267
d0b31f57 268}
269
270
271
272void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
273
274 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
275 if (!ESDfriend) {
276 Printf("ERROR: ESDfriend not available");
277 return;
278 }
279 //
280 UInt_t time = event->GetTimeStamp();
281 Int_t ntracks = event->GetNumberOfTracks();
592a0c8f 282 Int_t runNumber = event->GetRunNumber();
d0b31f57 283 //
284 // track loop
285 //
286 for (Int_t i=0;i<ntracks;++i) {
287
288 AliESDtrack *track = event->GetTrack(i);
289 if (!track) continue;
290
291 const AliExternalTrackParam * trackIn = track->GetInnerParam();
292 const AliExternalTrackParam * trackOut = track->GetOuterParam();
293 if (!trackIn) continue;
294 if (!trackOut) continue;
295
296 // calculate necessary track parameters
297 Double_t meanP = trackIn->GetP();
298 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
299 Int_t NclsDeDx = track->GetTPCNcls();
300
301 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
302 if (NclsDeDx < 60) continue;
303 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
304 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
305
306 // Get seeds
307 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
308 if (!friendTrack) continue;
309 TObject *calibObject;
310 AliTPCseed *seed = 0;
311 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
312 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
313 }
314
315 if (seed) {
592a0c8f 316 Double_t TPCsignal = GetTPCdEdx(seed);
317 fHistDeDxTotal->Fill(meanP, TPCsignal);
d0b31f57 318 //
319 if (fLowMemoryConsumption) {
320 if (meanP > 0.5 || meanP < 0.4) continue;
321 meanP = 0.45; // set all momenta to one in order to save memory
322 }
323 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
592a0c8f 324 Double_t vec[6] = {TPCsignal,time,2,meanDrift,meanP,runNumber};
d0b31f57 325 fHistGainTime->Fill(vec);
489dce1b 326 }
327
328 }
d0b31f57 329
489dce1b 330}
331
332
592a0c8f 333Float_t AliTPCcalibTimeGain::GetTPCdEdx(AliTPCseed * seed) {
334
335 Double_t signal = 0;
336 //
337 if (!fUseCookAnalytical) {
338 signal = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,fUseMax,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
339 } else {
340 signal = (1/fMIP)*seed->CookdEdxAnalytical(fLowerTrunc,fUpperTrunc,fUseMax);
341 }
342 //
343 return signal;
344}
345
346
347void AliTPCcalibTimeGain::AnalyzeRun(Int_t minEntries) {
489dce1b 348 //
349 //
350 //
489dce1b 351 if (fIsCosmic) {
d0b31f57 352 fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
353 fHistGainTime->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons
489dce1b 354 } else {
d0b31f57 355 fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
356 fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
489dce1b 357 }
489dce1b 358 //
592a0c8f 359 fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,minEntries,10);
489dce1b 360 //
361 return;
362}
363
364
592a0c8f 365TGraphErrors * AliTPCcalibTimeGain::GetGraphGainVsTime(Int_t runNumber, Int_t minEntries) {
366 //
367 //
368 //
369 if (runNumber == 0) {
370 if (!fGainVsTime) {
371 AnalyzeRun(minEntries);
372 }
373 } else {
374 // 1st check if the current run was cosmic or beam event
375 fHistGainTime->GetAxis(5)->SetRangeUser(runNumber,runNumber);
376 AnalyzeRun(minEntries);
377 }
378 if (fGainVsTime->GetN() == 0) return 0;
379 return fGainVsTime;
380}
381
489dce1b 382Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
383
384 TIterator* iter = li->MakeIterator();
385 AliTPCcalibTimeGain* cal = 0;
386
387 while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
388 if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
389 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
390 return -1;
391 }
392
393 // add histograms here...
394 if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
395 if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
396
397 }
398
399 return 0;
400
401}
402
403
592a0c8f 404AliSplineFit * AliTPCcalibTimeGain::MakeSplineFit(TGraphErrors * graph) {
405 //
406 //
407 //
408 AliSplineFit *fit = new AliSplineFit();
409 fit->SetGraph(graph);
410 fit->SetMinPoints(graph->GetN()+1);
411 fit->InitKnots(graph,2,0,0.001);
412 fit->SplineFit(0);
413 return fit;
414
415}
416
489dce1b 417
418void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
419
420 // Method for the correct logarithmic binning of histograms
421
422 TAxis *axis = h->GetAxis(axisDim);
423 int bins = axis->GetNbins();
424
425 Double_t from = axis->GetXmin();
426 Double_t to = axis->GetXmax();
427 Double_t *new_bins = new Double_t[bins + 1];
428
429 new_bins[0] = from;
430 Double_t factor = pow(to/from, 1./bins);
431
432 for (int i = 1; i <= bins; i++) {
433 new_bins[i] = factor * new_bins[i-1];
434 }
435 axis->Set(bins, new_bins);
436 delete new_bins;
437
438}
439
440
441void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
442
443 // Method for the correct logarithmic binning of histograms
444
445 TAxis *axis = h->GetXaxis();
446 int bins = axis->GetNbins();
447
448 Double_t from = axis->GetXmin();
449 Double_t to = axis->GetXmax();
450 Double_t *new_bins = new Double_t[bins + 1];
451
452 new_bins[0] = from;
453 Double_t factor = pow(to/from, 1./bins);
454
455 for (int i = 1; i <= bins; i++) {
456 new_bins[i] = factor * new_bins[i-1];
457 }
458 axis->Set(bins, new_bins);
459 delete new_bins;
460
461}
462
463
464
465void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
466 //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
467 const Double_t sigma = 0.06;
468
d0b31f57 469 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 470 BinLogX(histBG);
471
d0b31f57 472 TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
489dce1b 473 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
474 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
475 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
476 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
477 foElectron->SetParameters(ini);
478 foMuon->SetParameters(ini);
479 foPion->SetParameters(ini);
480 foKaon->SetParameters(ini);
481 foProton->SetParameters(ini);
482
483 TCanvas *CanvCheck1 = new TCanvas();
484 hist->Draw("colz");
485 foElectron->Draw("same");
486 foMuon->Draw("same");
487 foPion->Draw("same");
488 foKaon->Draw("same");
489 foProton->Draw("same");
490
491 // Loop over all points of the input histogram
492
493 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
494 Double_t x = hist->GetXaxis()->GetBinCenter(i);
495 for(Int_t j=1; j < hist->GetNbinsY(); j++) {
496 Long64_t n = hist->GetBin(i, j);
497 Double_t y = hist->GetYaxis()->GetBinCenter(j);
498 Double_t entries = hist->GetBinContent(n);
499 Double_t mass = 0;
500
501 // 1. identify protons
502 mass = 0.938;
503 if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
504 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
505 }
506
507 // 2. identify electrons
508 mass = 0.000511;
509 if (fIsCosmic) {
510 if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
511 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
512 }
513 } else {
d0b31f57 514 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 515 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
516 }
517 }
518
519 // 3. identify either muons or pions depending on cosmic or not
520 if (fIsCosmic) {
521 mass = 0.1056;
522 if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
523 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
524 }
525 } else {
526 mass = 0.1396;
527 if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
528 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
529 }
530 }
531
532 // 4. for pp also Kaons must be included
533 if (!fIsCosmic) {
534 mass = 0.4936;
535 if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
536 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
537 }
538 }
539 }
540 }
541
542 // Fit Aleph-Parameters to the obtained profile
543 TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
544 funcAlephD->SetParameters(ini);
545
546 TCanvas *CanvCheck2 = new TCanvas();
547 histBG->Draw();
548
549 //FitSlices
550 TObjArray * arr = new TObjArray();
551 histBG->FitSlicesY(0,0,-1,0,"QN",arr);
552 TH1D * fitMean = (TH1D*) arr->At(1);
553 fitMean->Draw("same");
554
555 funcAlephD->SetParLimits(2,1e-3,1e-7);
556 funcAlephD->SetParLimits(3,0.5,3.5);
557 funcAlephD->SetParLimits(4,0.5,3.5);
558 fitMean->Fit(funcAlephD, "QNR");
559 funcAlephD->Draw("same");
560
561 for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
562
563 foElectron->SetParameters(ini);
564 foMuon->SetParameters(ini);
565 foPion->SetParameters(ini);
566 foKaon->SetParameters(ini);
567 foProton->SetParameters(ini);
568
569 TCanvas *CanvCheck3 = new TCanvas();
570 hist->Draw("colz");
571 foElectron->Draw("same");
572 foMuon->Draw("same");
573 foPion->Draw("same");
574 foKaon->Draw("same");
575 foProton->Draw("same");
576
577 CanvCheck1->Print();
578 CanvCheck2->Print();
579 CanvCheck3->Print();
580
581 return;
582
583
584}