added weight option to onTheFly
[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
18//0. Libraries to lod
19gSystem->Load("libANALYSIS");
20gSystem->Load("libSTAT");
21gSystem->Load("libTPCcalib");
22
23
24//1. Do calibration ...
25//
26// compare reference
27
28//
d0b31f57 29//2. Visualize results
489dce1b 30//
31TFile fcalib("CalibObjects.root");
32TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
33AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain");
d0b31f57 34TGraphErrors * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10,0.2,0.8);
35gain->GetHistGainTime()->Projection(0,1)->Draw("colz")
36gr->SetMarkerStyle(25);
37gr->Draw("lp")
489dce1b 38
39
40//
41// MakeSlineFit example
42//
43AliSplineFit fit;
44fit.SetGraph(gr)
45fit->SetMinPoints(gr->GetN()+1);
46fit->InitKnots(gr,2,0,0.001)
47fit.SplineFit(0)
48fit.MakeDiffHisto(gr)->Draw();
49TGraph * grfit = fit.MakeGraph(gr->GetX()[0],gr->GetX()[gr->GetN()-1],50000,0);
50
51gr->SetMarkerStyle(25);
52gr->Draw("alp");
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),
111 fLowerTrunc(0),
112 fUpperTrunc(0),
113 fUseShapeNorm(0),
114 fUsePosNorm(0),
115 fUsePadNorm(0),
d0b31f57 116 fIsCosmic(0),
117 fLowMemoryConsumption(0)
489dce1b 118{
119 AliInfo("Default Constructor");
120}
121
122
123AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeGain)
124 :AliTPCcalibBase(),
125 fHistGainTime(0),
126 fGainVsTime(0),
127 fHistDeDxTotal(0),
128 fIntegrationTimeDeDx(0),
129 fMIP(0),
130 fLowerTrunc(0),
131 fUpperTrunc(0),
132 fUseShapeNorm(0),
133 fUsePosNorm(0),
134 fUsePadNorm(0),
d0b31f57 135 fIsCosmic(0),
136 fLowMemoryConsumption(0)
489dce1b 137 {
138
139 SetName(name);
140 SetTitle(title);
141
142 AliInfo("Non Default Constructor");
143
144 fIntegrationTimeDeDx = deltaIntegrationTimeGain;
145
146 Double_t deltaTime = EndTime - StartTime;
147
148
d0b31f57 149 // 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)
150 Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain);
151 Int_t binsGainTime[5] = {100, timeBins, 2, 25, 200};
152 Double_t xminGainTime[5] = {0.5, StartTime, 0.5, 0, 0.1};
489dce1b 153 Double_t xmaxGainTime[5] = { 4, EndTime, 2.5, 250, 50};
d0b31f57 154 fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx;dEdx,time,type,driftlength,momenta",5,binsGainTime,xminGainTime,xmaxGainTime);
489dce1b 155 BinLogX(fHistGainTime, 4);
156 //
157 fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000);
158 BinLogX(fHistDeDxTotal);
159
160 // default values for dE/dx
161 fMIP = 50.;
162 fLowerTrunc = 0.0;
163 fUpperTrunc = 0.7;
164 fUseShapeNorm = kTRUE;
165 fUsePosNorm = kFALSE;
166 fUsePadNorm = kFALSE;
167 //
168 fIsCosmic = kTRUE;
d0b31f57 169 fLowMemoryConsumption = kTRUE;
489dce1b 170 //
171
172 }
173
174
175
176AliTPCcalibTimeGain::~AliTPCcalibTimeGain(){
177 //
178 //
179 //
180}
181
182
183void AliTPCcalibTimeGain::Process(AliESDEvent *event) {
184 //
185 // main track loop
186 //
187 if (!event) {
188 Printf("ERROR: ESD not available");
189 return;
d0b31f57 190 }
191
192 if (fIsCosmic) { // this should be removed at some point based on trigger mask !?
193 ProcessCosmicEvent(event);
194 } else {
195 ProcessBeamEvent(event);
196 }
197
489dce1b 198
489dce1b 199
d0b31f57 200
201}
202
203
204void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) {
205
489dce1b 206 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
207 if (!ESDfriend) {
208 Printf("ERROR: ESDfriend not available");
209 return;
210 }
211 //
212 UInt_t time = event->GetTimeStamp();
d0b31f57 213 Int_t ntracks = event->GetNumberOfTracks();
489dce1b 214 //
215 // track loop
216 //
217 for (Int_t i=0;i<ntracks;++i) {
218
219 AliESDtrack *track = event->GetTrack(i);
220 if (!track) continue;
221
222 const AliExternalTrackParam * trackIn = track->GetInnerParam();
223 const AliExternalTrackParam * trackOut = track->GetOuterParam();
224 if (!trackIn) continue;
225 if (!trackOut) continue;
226
227 // calculate necessary track parameters
489dce1b 228 Double_t meanP = trackIn->GetP();
229 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
489dce1b 230 Int_t NclsDeDx = track->GetTPCNcls();
231
489dce1b 232 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
d0b31f57 233 if (NclsDeDx < 60) continue;
489dce1b 234 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
235 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
236
237 // Get seeds
238 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
239 if (!friendTrack) continue;
240 TObject *calibObject;
241 AliTPCseed *seed = 0;
242 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
243 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
244 }
245
246 if (seed) {
247 Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
489dce1b 248 fHistDeDxTotal->Fill(meanP, TPCsignalMax);
489dce1b 249 //
d0b31f57 250 if (fLowMemoryConsumption) {
251 if (meanP < 20) continue;
252 meanP = 30; // set all momenta to one in order to save memory
489dce1b 253 }
d0b31f57 254 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
255 Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP};
256 fHistGainTime->Fill(vec);
257 }
489dce1b 258
d0b31f57 259 }
489dce1b 260
d0b31f57 261}
262
263
264
265void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) {
266
267 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
268 if (!ESDfriend) {
269 Printf("ERROR: ESDfriend not available");
270 return;
271 }
272 //
273 UInt_t time = event->GetTimeStamp();
274 Int_t ntracks = event->GetNumberOfTracks();
275 //
276 // track loop
277 //
278 for (Int_t i=0;i<ntracks;++i) {
279
280 AliESDtrack *track = event->GetTrack(i);
281 if (!track) continue;
282
283 const AliExternalTrackParam * trackIn = track->GetInnerParam();
284 const AliExternalTrackParam * trackOut = track->GetOuterParam();
285 if (!trackIn) continue;
286 if (!trackOut) continue;
287
288 // calculate necessary track parameters
289 Double_t meanP = trackIn->GetP();
290 Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ());
291 Int_t NclsDeDx = track->GetTPCNcls();
292
293 // exclude tracks which do not look like primaries or are simply too short or on wrong sectors
294 if (NclsDeDx < 60) continue;
295 if (TMath::Abs(trackIn->GetTgl()) > 1) continue;
296 if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue;
297
298 // Get seeds
299 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
300 if (!friendTrack) continue;
301 TObject *calibObject;
302 AliTPCseed *seed = 0;
303 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
304 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
305 }
306
307 if (seed) {
308 Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0);
309 fHistDeDxTotal->Fill(meanP, TPCsignalMax);
310 //
311 if (fLowMemoryConsumption) {
312 if (meanP > 0.5 || meanP < 0.4) continue;
313 meanP = 0.45; // set all momenta to one in order to save memory
314 }
315 //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta
316 Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP};
317 fHistGainTime->Fill(vec);
489dce1b 318 }
319
320 }
d0b31f57 321
489dce1b 322}
323
324
325void AliTPCcalibTimeGain::Analyze() {
326 //
327 //
328 //
489dce1b 329 if (fIsCosmic) {
d0b31f57 330 fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics
331 fHistGainTime->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons
489dce1b 332 } else {
d0b31f57 333 fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data
334 fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions
489dce1b 335 }
489dce1b 336 //
d0b31f57 337 fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,2000,10);
489dce1b 338 //
339 return;
340}
341
342
343Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) {
344
345 TIterator* iter = li->MakeIterator();
346 AliTPCcalibTimeGain* cal = 0;
347
348 while ((cal = (AliTPCcalibTimeGain*)iter->Next())) {
349 if (!cal->InheritsFrom(AliTPCcalibTimeGain::Class())) {
350 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
351 return -1;
352 }
353
354 // add histograms here...
355 if (cal->GetHistGainTime()) fHistGainTime->Add(cal->GetHistGainTime());
356 if (cal->GetHistDeDxTotal()) fHistDeDxTotal->Add(cal->GetHistDeDxTotal());
357
358 }
359
360 return 0;
361
362}
363
364
489dce1b 365
366void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) {
367
368 // Method for the correct logarithmic binning of histograms
369
370 TAxis *axis = h->GetAxis(axisDim);
371 int bins = axis->GetNbins();
372
373 Double_t from = axis->GetXmin();
374 Double_t to = axis->GetXmax();
375 Double_t *new_bins = new Double_t[bins + 1];
376
377 new_bins[0] = from;
378 Double_t factor = pow(to/from, 1./bins);
379
380 for (int i = 1; i <= bins; i++) {
381 new_bins[i] = factor * new_bins[i-1];
382 }
383 axis->Set(bins, new_bins);
384 delete new_bins;
385
386}
387
388
389void AliTPCcalibTimeGain::BinLogX(TH1 *h) {
390
391 // Method for the correct logarithmic binning of histograms
392
393 TAxis *axis = h->GetXaxis();
394 int bins = axis->GetNbins();
395
396 Double_t from = axis->GetXmin();
397 Double_t to = axis->GetXmax();
398 Double_t *new_bins = new Double_t[bins + 1];
399
400 new_bins[0] = from;
401 Double_t factor = pow(to/from, 1./bins);
402
403 for (int i = 1; i <= bins; i++) {
404 new_bins[i] = factor * new_bins[i-1];
405 }
406 axis->Set(bins, new_bins);
407 delete new_bins;
408
409}
410
411
412
413void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) {
414 //{0.0762*MIP,10.632,1.34e-05,1.863,1.948}
415 const Double_t sigma = 0.06;
416
d0b31f57 417 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 418 BinLogX(histBG);
419
d0b31f57 420 TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100);
489dce1b 421 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100);
422 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100);
423 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100);
424 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.938,[0],[1],[2],[3],[4])",0.01,100);
425 foElectron->SetParameters(ini);
426 foMuon->SetParameters(ini);
427 foPion->SetParameters(ini);
428 foKaon->SetParameters(ini);
429 foProton->SetParameters(ini);
430
431 TCanvas *CanvCheck1 = new TCanvas();
432 hist->Draw("colz");
433 foElectron->Draw("same");
434 foMuon->Draw("same");
435 foPion->Draw("same");
436 foKaon->Draw("same");
437 foProton->Draw("same");
438
439 // Loop over all points of the input histogram
440
441 for(Int_t i=1; i < hist->GetNbinsX(); i++) {
442 Double_t x = hist->GetXaxis()->GetBinCenter(i);
443 for(Int_t j=1; j < hist->GetNbinsY(); j++) {
444 Long64_t n = hist->GetBin(i, j);
445 Double_t y = hist->GetYaxis()->GetBinCenter(j);
446 Double_t entries = hist->GetBinContent(n);
447 Double_t mass = 0;
448
449 // 1. identify protons
450 mass = 0.938;
451 if (TMath::Abs(y - foProton->Eval(x))/foProton->Eval(x) < 4*sigma && x < 0.65 ) {
452 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
453 }
454
455 // 2. identify electrons
456 mass = 0.000511;
457 if (fIsCosmic) {
458 if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 4*sigma && (x>0.25&&x<0.7) && fIsCosmic) {
459 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
460 }
461 } else {
d0b31f57 462 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 463 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
464 }
465 }
466
467 // 3. identify either muons or pions depending on cosmic or not
468 if (fIsCosmic) {
469 mass = 0.1056;
470 if (TMath::Abs(y - foMuon->Eval(x))/foMuon->Eval(x) < 4*sigma && x > 0.25 && x < 100 ) {
471 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
472 }
473 } else {
474 mass = 0.1396;
475 if (TMath::Abs(y - foPion->Eval(x))/foPion->Eval(x) < 4*sigma && x > 0.15 && x < 2) {
476 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
477 }
478 }
479
480 // 4. for pp also Kaons must be included
481 if (!fIsCosmic) {
482 mass = 0.4936;
483 if (TMath::Abs(y - foKaon->Eval(x))/foKaon->Eval(x) < 4*sigma && x > 0.2 && x < 0.3 ) {
484 for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y);
485 }
486 }
487 }
488 }
489
490 // Fit Aleph-Parameters to the obtained profile
491 TF1 * funcAlephD = new TF1("AlephParametrizationD", "AliExternalTrackParam::BetheBlochAleph(x,[0],[1],[2],[3],[4])",0.3,10000);
492 funcAlephD->SetParameters(ini);
493
494 TCanvas *CanvCheck2 = new TCanvas();
495 histBG->Draw();
496
497 //FitSlices
498 TObjArray * arr = new TObjArray();
499 histBG->FitSlicesY(0,0,-1,0,"QN",arr);
500 TH1D * fitMean = (TH1D*) arr->At(1);
501 fitMean->Draw("same");
502
503 funcAlephD->SetParLimits(2,1e-3,1e-7);
504 funcAlephD->SetParLimits(3,0.5,3.5);
505 funcAlephD->SetParLimits(4,0.5,3.5);
506 fitMean->Fit(funcAlephD, "QNR");
507 funcAlephD->Draw("same");
508
509 for(Int_t i=0;i<5;i++) ini[i] = funcAlephD->GetParameter(i);
510
511 foElectron->SetParameters(ini);
512 foMuon->SetParameters(ini);
513 foPion->SetParameters(ini);
514 foKaon->SetParameters(ini);
515 foProton->SetParameters(ini);
516
517 TCanvas *CanvCheck3 = new TCanvas();
518 hist->Draw("colz");
519 foElectron->Draw("same");
520 foMuon->Draw("same");
521 foPion->Draw("same");
522 foKaon->Draw("same");
523 foProton->Draw("same");
524
525 CanvCheck1->Print();
526 CanvCheck2->Print();
527 CanvCheck3->Print();
528
529 return;
530
531
532}