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