git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerParam.cxx
CommitLineData
6eb67451 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. *
c68184b6 14 **************************************************************************/
15
16/*
17$Log$
6eb69b9a 18Revision 1.7 2002/04/10 16:30:07 kowal2
19logs added
20
c68184b6 21*/
22
23
6eb69b9a 24/**************************************************************************
6eb67451 25 * *
26 * This class builds AliTPCtrack objects from generated tracks to feed *
27 * ITS tracking (V2). The AliTPCtrack is built from its first hit in *
28 * the TPC. The track is assigned a Kalman-like covariance matrix *
29 * depending on its pT and pseudorapidity and track parameters are *
70521312 30 * smeared according to this covariance matrix. *
6eb69b9a 31 * Output file contains sorted tracks, ready for matching with ITS. *
6eb67451 32 * *
e130146c 33 * For details: *
34 * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm *
35 * *
36 * Test macro is: AliBarrelRec_TPCparam.C *
6eb67451 37 * *
6eb69b9a 38 * 2002/04/28: Major upgrade of the class *
39 * - Covariance matrices and pulls are now separeted for pions, kaons *
40 * and electrons. *
41 * - A parameterization of dE/dx in the TPC has been included and it is *
42 * used to assign a mass to each track according to a rough dE/dx PID. *
43 * - All the "numbers" have been moved to the file with the DataBase, *
44 * they are read as objects of the class AliTPCkineGrid, and assigned *
45 * to data memebers of the class AliTPCtrackerParam. *
46 * - All the code necessary to create a BataBase has been included in *
47 * class (see the macro AliTPCtrackingParamDB.C for the details). *
48 * *
49 * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it *
6eb67451 50 * *
51 **************************************************************************/
6eb69b9a 52#include <TChain.h>
53#include <TF1.h>
54#include <TGraph.h>
55#include <TGraphErrors.h>
56#include <TLegend.h>
57#include <TLine.h>
58#include <TMatrixD.h>
59#include <TNtuple.h>
60#include <TSystem.h>
6eb67451 61#include "alles.h"
6eb69b9a 62#include "AliGausCorr.h"
6eb67451 63#include "AliKalmanTrack.h"
6eb69b9a 64#include "AliMagF.h"
70521312 65#include "AliMagFCM.h"
6eb69b9a 66#include "AliTPCkineGrid.h"
67#include "AliTPCtrack.h"
68#include "AliTPCtrackerParam.h"
69
70Double_t RegFunc(Double_t *x,Double_t *par) {
71// This is the function used to regularize the covariance matrix
72 Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
73 TMath::Power(x[1],par[3]);
74 return value;
75}
76Double_t FitRegFunc(Double_t *x,Double_t *par) {
77// This is the function used for the fit of the covariance
78// matrix as a function of the momentum.
79// For cosl the average value 0.87 is used.
80 Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
81 TMath::Power(0.87,par[3]);
82 return value;
83}
70521312 84
6eb69b9a 85// structure for DB building
86typedef struct {
87 Int_t pdg;
88 Int_t bin;
89 Double_t r;
90 Double_t p;
91 Double_t pt;
92 Double_t cosl;
93 Double_t eta;
94 Double_t dpt;
95 Double_t dP0,dP1,dP2,dP3,dP4;
96 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
97 Double_t dEdx;
98} COMPTRACK;
99// cov matrix structure
100typedef struct {
101 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
102} COVMATRIX;
70521312 103
6eb67451 104ClassImp(AliTPCtrackerParam)
105
6eb69b9a 106//-----------------------------------------------------------------------------
107AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
6eb67451 108{
6eb69b9a 109//-----------------------------------------------------------------------------
e130146c 110// This is the class conctructor
6eb69b9a 111//-----------------------------------------------------------------------------
e130146c 112
6eb69b9a 113 fBz = Bz; // value of the z component of L3 field (Tesla)
114 fColl = coll; // collision code (0: PbPb6000)
115 fSelAndSmear = kTRUE; // by default selection and smearing are done
116
117 if(fBz!=0.4) {
118 cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n";
119 cerr<<" Available: 0.4"<<endl;
120 }
121 if(fColl!=0) {
122 cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n";
123 cerr<<" Available: 0 -> PbPb6000"<<endl;
124 }
125
126 fDBfileName = gSystem->Getenv("ALICE_ROOT");
127 fDBfileName.Append("/TPC/CovMatrixDB_");
128 if(fColl==0) fDBfileName.Append("PbPb6000");
129 if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
6eb67451 130}
6eb69b9a 131//-----------------------------------------------------------------------------
6eb67451 132AliTPCtrackerParam::~AliTPCtrackerParam()
133{}
6eb69b9a 134//-----------------------------------------------------------------------------
6eb67451 135
136Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
137{
6eb69b9a 138//-----------------------------------------------------------------------------
e130146c 139// This function creates the TPC parameterized tracks
6eb69b9a 140//-----------------------------------------------------------------------------
141
142 TFile *fileDB=0;
143 TTree *covTreePi[50];
144 TTree *covTreeKa[50];
145 TTree *covTreeEl[50];
146
147 if(fSelAndSmear) {
148 cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
149 fDBfileName.Data()<<"\n+++\n";
150 // Read paramters from DB file
151 if(!ReadAllData(fDBfileName.Data())) {
152 cerr<<"AliTPCtrackerParam::BuildTPCtracks:
153 Could not read data from DB\n\n"; return 1;
154 }
155 // Read the trees with regularized cov. matrices from DB
156 TString str;
157 fileDB = TFile::Open(fDBfileName.Data());
158 Int_t nBinsPi = fDBgridPi.GetTotBins();
159 for(Int_t l=0; l<nBinsPi; l++) {
160 str = "/CovMatrices/Pions/CovTreePi_bin";
161 str+=l;
162 covTreePi[l] = (TTree*)fileDB->Get(str.Data());
163 }
164 Int_t nBinsKa = fDBgridKa.GetTotBins();
165 for(Int_t l=0; l<nBinsKa; l++) {
166 str = "/CovMatrices/Kaons/CovTreeKa_bin";
167 str+=l;
168 covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
169 }
170 Int_t nBinsEl = fDBgridEl.GetTotBins();
171 for(Int_t l=0; l<nBinsEl; l++) {
172 str = "/CovMatrices/Electrons/CovTreeEl_bin";
173 str+=l;
174 covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
175 }
e130146c 176
6eb69b9a 177 } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
6eb67451 178
179 TFile *infile=(TFile*)inp;
6eb69b9a 180 infile->cd();
6eb67451 181
6eb67451 182 // Get gAlice object from file
183 if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
184 cerr<<"gAlice has not been found on galice.root !\n";
185 return 1;
186 }
187
6eb69b9a 188 // Check if value in the galice file is equal to selected one (fBz)
70521312 189 AliMagFCM *fiel = (AliMagFCM*)gAlice->Field();
190 Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
191 printf("Magnetic field is %6.2f Tesla\n",fieval);
192 if(fBz!=fieval) {
193 cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!"<<endl;
194 cerr<<"Field selected is: "<<fBz<<" T\n";
195 cerr<<"Field found on file is: "<<fieval<<" T\n";
196 return 0;
197 }
198
6eb69b9a 199 // Get TPC detector
200 AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
201 Int_t ver = TPC->IsVersion();
202 cerr<<"+++ TPC version "<<ver<<" has been found !\n";
203 AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
204 if(digp){
205 delete digp;
206 digp = new AliTPCParamSR();
207 }
208 else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
209
210 if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
211 TPC->SetParam(digp);
212
213 // Set the conversion constant between curvature and Pt
70521312 214 AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
6eb67451 215
6eb69b9a 216 TParticle *Part=0;
217 AliTPChit *tpcHit=0;
218 AliTPCtrack *tpctrack=0;
219 Double_t hPx,hPy,hPz,hPt,hEta,xg,yg,zg,xl,yl,zl;
220 Double_t alpha;
221 Float_t cosAlpha,sinAlpha;
222 Int_t bin,label,pdg,charge;
223 Int_t tracks=0;
224 Int_t nParticles,nTracks,arrentr;
225 Char_t tname[100];
226 //Int_t nSel=0,nAcc=0;
6eb67451 227
228 // loop over first n events in file
229 for(Int_t evt=0; evt<n; evt++){
230 cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
231
6eb67451 232 // tree for TPC tracks
6eb67451 233 sprintf(tname,"TreeT_TPC_%d",evt);
234 TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
235 tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
236
237 // array for TPC tracks
238 TObjArray tarray(20000);
239
6eb69b9a 240 // get the particles stack
241 nParticles = gAlice->GetEvent(evt);
242 Bool_t done[nParticles];
243 Int_t pdgCodes[nParticles];
6eb67451 244
6eb69b9a 245 // loop on particles and store pdg codes
246 for(Int_t l=0; l<nParticles; l++) {
247 Part = gAlice->Particle(l);
248 pdgCodes[l] = Part->GetPdgCode();
249 done[l] = kFALSE;
7a09f434 250 }
6eb67451 251
252 // Get TreeH with hits
253 TTree *TH=gAlice->TreeH();
6eb69b9a 254 nTracks=(Int_t)TH->GetEntries();
255 cerr<<"+++\n+++ Number of particles in event "<<evt<<": "<<nParticles<<
256 "\n+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
257 "\n+++\n\n";
6eb67451 258
259 // loop over entries in TreeH
6eb69b9a 260 for(Int_t l=0; l<nTracks; l++) {
261 if(l%1000==0) cerr<<" --- Processing primary track "
262 <<l<<" of "<<nTracks<<" ---\r";
6eb67451 263 TPC->ResetHits();
6eb69b9a 264 TH->GetEvent(l);
6eb67451 265 // Get FirstHit
266 tpcHit=(AliTPChit*)TPC->FirstHit(-1);
267 for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
268 if(tpcHit->fQ !=0.) continue;
269 // Get particle momentum at hit
270 hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
271 hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
272 // reject hits with Pt<mag*0.45 GeV/c
273 if(hPt<(fBz*0.45)) continue;
274
275 // Get track label
276 label=tpcHit->Track();
277 // check if this track has already been processed
278 if(done[label]) continue;
6eb69b9a 279 // PDG code & electric charge
280 pdg = pdgCodes[label];
e130146c 281 if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
282 else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
6eb67451 283 else continue;
6eb69b9a 284 pdg = TMath::Abs(pdg);
285 if(pdg>3000) pdg=211;
286 if(fSelAndSmear) SetParticle(pdg);
6eb67451 287
288 if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
289 if(tpcHit->fQ != 0.) continue;
290 // Get global coordinates of hit
291 xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
292 if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
293
294 // Get TPC sector, Alpha angle and local coordinates
6eb69b9a 295 // cerr<<"Sector "<<tpcHit->fSector<<endl;
6eb67451 296 digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
e130146c 297 alpha = TMath::ATan2(sinAlpha,cosAlpha);
6eb67451 298 xl = xg*cosAlpha + yg*sinAlpha;
299 yl =-xg*sinAlpha + yg*cosAlpha;
300 zl = zg;
301 //printf("Alpha %f xl %f yl %f zl %f\n",Alpha,xl,yl,zl);
302
303 // reject tracks which are not in the TPC acceptance
304 if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
305
6eb69b9a 306 hEta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(hPz/hPt)));
6eb67451 307
308 // Apply selection according to TPC efficiency
e130146c 309 //if(TMath::Abs(pdg)==211) nAcc++;
6eb69b9a 310 if(fSelAndSmear && !SelectedTrack(hPt,hEta)) continue;
e130146c 311 //if(TMath::Abs(pdg)==211) nSel++;
6eb67451 312
6eb67451 313 // create AliTPCtrack object
6eb69b9a 314 BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge);
315
316 if(fSelAndSmear) {
317 bin = fDBgrid->GetBin(hPt,hEta);
318 switch (pdg) {
319 case 211:
320 fCovTree = covTreePi[bin];
321 break;
322 case 321:
323 fCovTree = covTreeKa[bin];
324 break;
325 case 2212:
326 fCovTree = covTreePi[bin];
327 break;
328 case 11:
329 fCovTree = covTreeEl[bin];
330 break;
331 case 13:
332 fCovTree = covTreePi[bin];
333 break;
334 }
335 // deal with covariance matrix and smearing of parameters
336 CookTrack(hPt,hEta);
337
338 // assign the track a dE/dx and make a rough PID
339 CookdEdx(hPt,hEta);
340 }
6eb67451 341
342 // put track in array
6eb69b9a 343 AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
344 iotrack->SetLabel(label);
345 tarray.AddLast(iotrack);
346 // Mark track as "done" and register the pdg code
347 done[label]=kTRUE;
348 tracks++;
6eb67451 349 }
350
351 } // loop over entries in TreeH
352
6eb67451 353 // sort array with TPC tracks (decreasing pT)
6eb69b9a 354 tarray.Sort();
6eb67451 355
6eb69b9a 356 arrentr = tarray.GetEntriesFast();
6eb67451 357 for(Int_t l=0; l<arrentr; l++) {
6eb69b9a 358 tpctrack=(AliTPCtrack*)tarray.UncheckedAt(l);
6eb67451 359 tracktree->Fill();
360 }
361
362 // write the tree with tracks in the output file
363 out->cd();
364 tracktree->Write();
365
366 delete tracktree;
367
368 printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
6eb69b9a 369 //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
6eb67451 370
371 } // loop on events
372
6eb69b9a 373 if(fileDB) fileDB->Close();
374
6eb67451 375 return 0;
376}
6eb69b9a 377//-----------------------------------------------------------------------------
378void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
379//-----------------------------------------------------------------------------
380// This function computes the dE/dx for pions, kaons, protons and electrons,
381// in the [pT,eta] bins.
382// Input file is CovMatrix_AllEvts.root.
383//-----------------------------------------------------------------------------
384
385 gStyle->SetOptStat(0);
386 gStyle->SetOptFit(10001);
387
388 Char_t *part="PIONS";
389 Double_t ymax=500.;
390
391 // create a chain with compared tracks
392 TChain cmptrkchain("cmptrktree");
393 cmptrkchain.Add("CovMatrix_AllEvts_1.root");
394 cmptrkchain.Add("CovMatrix_AllEvts_2.root");
395 cmptrkchain.Add("CovMatrix_AllEvts_3.root");
396
397 COMPTRACK cmptrk;
398 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
399 Int_t entries = (Int_t)cmptrkchain.GetEntries();
400 cerr<<" Number of entries: "<<entries<<endl;
401
402 InitializeKineGrid("DB","points");
403 InitializeKineGrid("dEdx","bins");
404
405 switch(pdg) {
406 case 211:
407 part = "PIONS";
408 ymax = 100.;
409 break;
410 case 321:
411 part = "KAONS";
412 ymax = 300.;
413 break;
414 case 2212:
415 part = "PROTONS";
416 ymax = 500.;
417 break;
418 case 11:
419 part = "ELECTRONS";
420 ymax = 100.;
421 break;
422 }
423
424 SetParticle(pdg);
425
426 const Int_t nTotBins = fDBgrid->GetTotBins();
427
428 cerr<<" Fit bins: "<<nTotBins<<endl;
429
430 Int_t bin=0;
431 Int_t n[nTotBins];
432 Double_t p[nTotBins],ep[nTotBins],mean[nTotBins],sigma[nTotBins];
433
434 for(Int_t l=0; l<nTotBins; l++) {
435 n[l] = 1; // set to 1 to avoid divisions by 0
436 p[l] = mean[l] = sigma[l] = ep[l] = 0.;
437 }
438
439 // loop on chain entries for the mean
440 for(Int_t l=0; l<entries; l++) {
441 cmptrkchain.GetEvent(l);
442 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
443 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
444 p[bin] += cmptrk.p;
445 mean[bin] += cmptrk.dEdx;
446 n[bin]++;
447 } // loop on chain entries
448
449 for(Int_t l=0; l<nTotBins; l++) {
450 p[l] /= n[l];
451 mean[l] /= n[l];
452 n[l] = 1; // set to 1 to avoid divisions by 0
453 }
454
455 // loop on chain entries for the sigma
456 for(Int_t l=0; l<entries; l++) {
457 cmptrkchain.GetEvent(l);
458 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
459 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
460 if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
461 n[bin]++;
462 sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
463 } // loop on chain entries
464
465 for(Int_t l=0; l<nTotBins; l++) {
466 sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
467 }
468
469 // create the canvas
470 TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700);
471
472 // create the graph for dEdx vs p
473 TGraphErrors *gr = new TGraphErrors(nTotBins,p,mean,ep,sigma);
474 TString title(" : dE/dx vs momentum"); title.Prepend(part);
475 TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
476 frame->SetXTitle("p [GeV/c]");
477 frame->SetYTitle("dE/dx [a.u.]");
478 canv->SetLogx();
479 frame->Draw();
480 gr->Draw("P");
481
482 switch(pdg) {
483 case 211:
484 for(Int_t i=0; i<nTotBins; i++) {
485 fdEdxMeanPi.SetParam(i,mean[i]);
486 fdEdxRMSPi.SetParam(i,sigma[i]);
487 }
488 break;
489 case 321:
490 for(Int_t i=0; i<nTotBins; i++) {
491 fdEdxMeanKa.SetParam(i,mean[i]);
492 fdEdxRMSKa.SetParam(i,sigma[i]);
493 }
494 break;
495 case 2212:
496 for(Int_t i=0; i<nTotBins; i++) {
497 fdEdxMeanPr.SetParam(i,mean[i]);
498 fdEdxRMSPr.SetParam(i,sigma[i]);
499 }
500 break;
501 case 11:
502 for(Int_t i=0; i<nTotBins; i++) {
503 fdEdxMeanEl.SetParam(i,mean[i]);
504 fdEdxRMSEl.SetParam(i,sigma[i]);
505 }
506 break;
507 }
508
509 // write results to file
510 WritedEdx(outName,pdg);
511
512
513 return;
514}
515//-----------------------------------------------------------------------------
516void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
517//-----------------------------------------------------------------------------
518// This function computes the pulls for pions, kaons and electrons,
519// in the [pT,eta] bins.
520// Input file is CovMatrix_AllEvts.root.
521// Output file is pulls.root.
522//-----------------------------------------------------------------------------
523
524 // create a chain with compared tracks
525 TChain cmptrkchain("cmptrktree");
526 cmptrkchain.Add("CovMatrix_AllEvts_1.root");
527 cmptrkchain.Add("CovMatrix_AllEvts_2.root");
528 cmptrkchain.Add("CovMatrix_AllEvts_3.root");
529
530 COMPTRACK cmptrk;
531 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
532 Int_t entries = (Int_t)cmptrkchain.GetEntries();
533 cerr<<" Number of entries: "<<entries<<endl;
534
535 Int_t thisPdg=0;
536 Char_t hname[100], htitle[100];
537 Int_t nTotBins,bin;
538
539 AliTPCkineGrid pulls[5];
540 TH1F *hDum = new TH1F("name","title",100,-7.,7.);
541 TF1 *g = new TF1("g","gaus");
542
543 InitializeKineGrid("pulls","bins");
544 InitializeKineGrid("DB","points");
545
546
547
548 // loop on the particles Pi,Ka,El
549 for(Int_t part=0; part<3; part++) {
550
551 switch (part) {
552 case 0: // pions
553 thisPdg=211;
554 cerr<<" Processing pions ...\n";
555 break;
556 case 1: // kaons
557 thisPdg=321;
558 cerr<<" Processing kaons ...\n";
559 break;
560
561 case 2: // electrons
562 thisPdg=11;
563 cerr<<" Processing electrons ...\n";
564 break;
565 }
566
567 SetParticle(thisPdg);
6eb67451 568
6eb69b9a 569 for(Int_t i=0;i<5;i++) {
570 pulls[i].~AliTPCkineGrid();
571 new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
572 }
573 nTotBins = fDBgrid->GetTotBins();
574
575 // create histograms for the all the bins
576 TH1F *hPulls0_=NULL;
577 TH1F *hPulls1_=NULL;
578 TH1F *hPulls2_=NULL;
579 TH1F *hPulls3_=NULL;
580 TH1F *hPulls4_=NULL;
581
582 hPulls0_ = new TH1F[nTotBins];
583 hPulls1_ = new TH1F[nTotBins];
584 hPulls2_ = new TH1F[nTotBins];
585 hPulls3_ = new TH1F[nTotBins];
586 hPulls4_ = new TH1F[nTotBins];
587
588
589 for(Int_t i=0; i<nTotBins; i++) {
590 sprintf(hname,"hPulls0_%d",i);
591 sprintf(htitle,"P0 pulls for bin %d",i);
592 hDum->SetName(hname); hDum->SetTitle(htitle);
593 hPulls0_[i] = *hDum;
594 sprintf(hname,"hPulls1_%d",i);
595 sprintf(htitle,"P1 pulls for bin %d",i);
596 hDum->SetName(hname); hDum->SetTitle(htitle);
597 hPulls1_[i] = *hDum;
598 sprintf(hname,"hPulls2_%d",i);
599 sprintf(htitle,"P2 pulls for bin %d",i);
600 hDum->SetName(hname); hDum->SetTitle(htitle);
601 hPulls2_[i] = *hDum;
602 sprintf(hname,"hPulls3_%d",i);
603 sprintf(htitle,"P3 pulls for bin %d",i);
604 hDum->SetName(hname); hDum->SetTitle(htitle);
605 hPulls3_[i] = *hDum;
606 sprintf(hname,"hPulls4_%d",i);
607 sprintf(htitle,"P4 pulls for bin %d",i);
608 hDum->SetName(hname); hDum->SetTitle(htitle);
609 hPulls4_[i] = *hDum;
610 }
611
612 // loop on chain entries
613 for(Int_t i=0; i<entries; i++) {
614 cmptrkchain.GetEvent(i);
615 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
616 // fill histograms with the pulls
617 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
618 hPulls0_[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
619 hPulls1_[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
620 hPulls2_[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
621 hPulls3_[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
622 hPulls4_[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
623 } // loop on chain entries
624
625 // compute the sigma of the distributions
626 for(Int_t i=0; i<nTotBins; i++) {
627 if(hPulls0_[i].GetEntries()>1000) {
628 g->SetRange(-3.*hPulls0_[i].GetRMS(),3.*hPulls0_[i].GetRMS());
629 hPulls0_[i].Fit("g","R,Q,N");
630 pulls[0].SetParam(i,g->GetParameter(2));
631 } else pulls[0].SetParam(i,-1.);
632 if(hPulls1_[i].GetEntries()>1000) {
633 g->SetRange(-3.*hPulls1_[i].GetRMS(),3.*hPulls1_[i].GetRMS());
634 hPulls1_[i].Fit("g","R,Q,N");
635 pulls[1].SetParam(i,g->GetParameter(2));
636 } else pulls[1].SetParam(i,-1.);
637 if(hPulls2_[i].GetEntries()>1000) {
638 g->SetRange(-3.*hPulls2_[i].GetRMS(),3.*hPulls2_[i].GetRMS());
639 hPulls2_[i].Fit("g","R,Q,N");
640 pulls[2].SetParam(i,g->GetParameter(2));
641 } else pulls[2].SetParam(i,-1.);
642 if(hPulls3_[i].GetEntries()>1000) {
643 g->SetRange(-3.*hPulls3_[i].GetRMS(),3.*hPulls3_[i].GetRMS());
644 hPulls3_[i].Fit("g","R,Q,N");
645 pulls[3].SetParam(i,g->GetParameter(2));
646 } else pulls[3].SetParam(i,-1.);
647 if(hPulls4_[i].GetEntries()>1000) {
648 g->SetRange(-3.*hPulls4_[i].GetRMS(),3.*hPulls4_[i].GetRMS());
649 hPulls4_[i].Fit("g","R,Q,N");
650 pulls[4].SetParam(i,g->GetParameter(2));
651 } else pulls[4].SetParam(i,-1.);
652 } // loop on bins
653
654
655 switch (part) {
656 case 0: // pions
657 for(Int_t i=0;i<5;i++) {
658 fPullsPi[i].~AliTPCkineGrid();
659 new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
660 }
661 break;
662 case 1: // kaons
663 for(Int_t i=0;i<5;i++) {
664 fPullsKa[i].~AliTPCkineGrid();
665 new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
666 }
667 break;
668 case 2: // electrons
669 for(Int_t i=0;i<5;i++) {
670 fPullsEl[i].~AliTPCkineGrid();
671 new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
672 }
673 break;
674 }
675
676 delete [] hPulls0_;
677 delete [] hPulls1_;
678 delete [] hPulls2_;
679 delete [] hPulls3_;
680 delete [] hPulls4_;
681
682 } // loop on particle species
683
684 // write pulls to file
685 WritePulls(outName);
686
687
688 return;
689}
690//-----------------------------------------------------------------------------
691void AliTPCtrackerParam::BuildTrack(Double_t alpha,
692 Double_t x,Double_t y,Double_t z,
693 Double_t px,Double_t py,
694 Double_t pz,Double_t pt,
695 Int_t ch) {
696//-----------------------------------------------------------------------------
e130146c 697// This function uses GEANT info to set true track parameters
6eb69b9a 698//-----------------------------------------------------------------------------
6eb67451 699 Double_t xref = x;
700 Double_t xx[5],cc[15];
701 cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
702 cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
703
704 // Magnetic field
e130146c 705 TVector3 bfield(0.,0.,fBz);
6eb67451 706
707
708 // radius [cm] of track projection in (x,y)
e130146c 709 Double_t rho = pt*100./0.299792458/bfield.Z();
6eb67451 710 // center of track projection in local reference frame
e130146c 711 TVector3 hmom,hpos;
6eb67451 712
713
714 // position (local) and momentum (local) at the hit
715 // in the bending plane (z=0)
e130146c 716 hpos.SetXYZ(x,y,0.);
717 hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.);
718 TVector3 vrho = hmom.Cross(bfield);
719 vrho *= ch;
720 vrho.SetMag(rho);
6eb67451 721
e130146c 722 TVector3 vcenter = hpos+vrho;
6eb67451 723
e130146c 724 Double_t x0 = vcenter.X();
6eb67451 725
726 // fX = xref X-coordinate of this track (reference plane)
727 // fAlpha = Alpha Rotation angle the local (TPC sector)
728 // fP0 = YL Y-coordinate of a track
729 // fP1 = ZG Z-coordinate of a track
730 // fP2 = C*x0 x0 is center x in rotated frame
731 // fP3 = Tgl tangent of the track momentum dip angle
732 // fP4 = C track curvature
733 xx[0] = y;
734 xx[1] = z;
e130146c 735 xx[3] = pz/pt;
6eb67451 736 xx[4] = -ch/rho;
737 xx[2] = xx[4]*x0;
738
739 // create the object AliTPCtrack
6eb69b9a 740 AliTPCtrack track(0,xx,cc,xref,alpha);
741 new(&fTrack) AliTPCtrack(track);
6eb67451 742
6eb69b9a 743 return;
6eb67451 744}
6eb69b9a 745//-----------------------------------------------------------------------------
746void AliTPCtrackerParam::CompareTPCtracks(
747 const Char_t* galiceName="galice.root",
748 const Char_t* trkGeaName="AliTPCtracksGeant.root",
749 const Char_t* trkKalName="AliTPCtracksSorted.root",
750 const Char_t* covmatName="CovMatrix.root",
751 const Char_t* tpceffName="TPCeff.dat") const {
752//-----------------------------------------------------------------------------
753// This function compares tracks from TPC Kalman Filter V2 with
754// true tracks at TPC 1st hit. It gives:
755// - a tree with Kalman cov. matrix elements, resolutions, dEdx
756// - the efficiencies as a function of particle type, pT, eta
757//-----------------------------------------------------------------------------
758
759 TFile *kalFile = TFile::Open(trkKalName);
760 TFile *geaFile = TFile::Open(trkGeaName);
761 TFile *galiceFile = TFile::Open(galiceName);
762
763 // particles from TreeK
764 AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
765 const Int_t nparticles = gAlice->GetEvent(0);
766
767 Int_t kalLab[nparticles];
768 for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1;
769
6eb67451 770
6eb69b9a 771 // tracks from Kalman
772 TTree *kaltree=(TTree*)kalFile->Get("TreeT_TPC_0");
773 AliTPCtrack *kaltrack=new AliTPCtrack;
774 kaltree->SetBranchAddress("tracks",&kaltrack);
775 Int_t kalEntries = (Int_t)kaltree->GetEntries();
776
777 // tracks from 1st hit
778 TTree *geatree=(TTree*)geaFile->Get("TreeT_TPC_0");
779 AliTPCtrack *geatrack=new AliTPCtrack;
780 geatree->SetBranchAddress("tracks",&geatrack);
781 Int_t geaEntries = (Int_t)geatree->GetEntries();
782
783 cerr<<"+++\n+++ Number of tracks: TPC Kalman = "<<kalEntries<<endl<<"+++ TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
784
785 // set pointers for TPC tracks:
786 // the entry number of the track labelled l is stored in kalLab[l]
787 Int_t fake=0,mult=0;
788 for (Int_t j=0; j<kalEntries; j++) {
789 kaltree->GetEvent(j);
790 if(kaltrack->GetLabel()>=0) {
791 if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
792 kalLab[kaltrack->GetLabel()] = j;
793 }
794 else fake++;
795 }
796
797 cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
798 cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
799
800 TParticle *Part;
801 Int_t label;
802 Double_t cc[15],dAlpha;
803 // --- [Pt,eta] binning ---
804 //
805 // |eta|<0.3 0.3<=|eta|<0.6 0.6<=|eta|
806 // Pt<0.3 0 1 2
807 // 0.3<=Pt<0.5 3 4 5
808 // 0.5<=Pt<1.0 6 7 8
809 // 1.0<=Pt<1.5 9 10 11
810 // 1.5<=Pt<3.0 12 13 14
811 // 3.0<=Pt<5.0 15 16 17
812 // 5.0<=Pt<7.0 18 19 20
813 // 7.0<=Pt<15. 21 22 23
814 // 15.<=Pt 24 25 26
815 //
816 Int_t pi=0,ka=0,mu=0,el=0,pr=0;
817 Int_t geaPi[27],geaKa[27],geaPr[27],geaEl[27],geaMu[27];
818 Int_t kalPi[27],kalKa[27],kalPr[27],kalEl[27],kalMu[27];
819 Float_t effPi[27],effKa[27],effPr[27],effEl[27],effMu[27];
820 Int_t bin;
821 for(Int_t j=0; j<27; j++) {
822 geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
823 kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
824 effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
825 }
e130146c 826
6eb69b9a 827 COMPTRACK cmptrk;
828 // create the tree for comparison results
829 TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
830 cmptrktree->Branch("comptracks",&cmptrk,"pdg/I:bin:r/D:p:pt:cosl:eta:dpt:dP0:dP1:dP2:dP3:dP4:c00:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44:dEdx");
6eb67451 831
6eb69b9a 832 cerr<<"Doing track comparison...\n";
833 // loop on tracks at 1st hit
834 for(Int_t j=0; j<geaEntries; j++) {
835 geatree->GetEvent(j);
836
837 label = geatrack->GetLabel();
838 Part = (TParticle*)gAlice->Particle(label);
839 cmptrk.pdg = Part->GetPdgCode();
840 cmptrk.eta = Part->Eta();
841 cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
842
843 cmptrk.pt = 1/TMath::Abs(geatrack->Get1Pt());
844 cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
845 cmptrk.p = cmptrk.pt/cmptrk.cosl;
846
847
848 bin = GetBin(cmptrk.pt,cmptrk.eta);
849 cmptrk.bin = bin;
850 if(abs(cmptrk.pdg)==211) geaPi[bin]++;
851 if(abs(cmptrk.pdg)==321) geaKa[bin]++;
852 if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
853 if(abs(cmptrk.pdg)==11) geaEl[bin]++;
854 if(abs(cmptrk.pdg)==13) geaMu[bin]++;
855
856
857 // check if there is the corresponding track in the TPC kalman and get it
858 if(kalLab[label]==-1) continue;
859 kaltree->GetEvent(kalLab[label]);
860
861 // and go on only if it has xref = 84.57 cm (inner pad row)
862 if(kaltrack->GetX()>90.) continue;
863
864 if(abs(cmptrk.pdg)==211) { kalPi[bin]++; pi++; }
865 if(abs(cmptrk.pdg)==321) { kalKa[bin]++; ka++; }
866 if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
867 if(abs(cmptrk.pdg)==11) { kalEl[bin]++; el++; }
868 if(abs(cmptrk.pdg)==13) { kalMu[bin]++; mu++; }
869
870 kaltrack->PropagateTo(geatrack->GetX());
871
872 cmptrk.dEdx = kaltrack->GetdEdx();
873
874 // compute errors on parameters
875 dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
876 if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
877
878 cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
879 cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
880 cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
881 cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
882 cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
883 cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
884
885 // get covariance matrix
886 // beware: lines 3 and 4 are inverted!
887 kaltrack->GetCovariance(cc);
888
889 cmptrk.c00 = cc[0];
890 cmptrk.c10 = cc[1];
891 cmptrk.c11 = cc[2];
892 cmptrk.c20 = cc[3];
893 cmptrk.c21 = cc[4];
894 cmptrk.c22 = cc[5];
895 cmptrk.c30 = cc[10];
896 cmptrk.c31 = cc[11];
897 cmptrk.c32 = cc[12];
898 cmptrk.c33 = cc[14];
899 cmptrk.c40 = cc[6];
900 cmptrk.c41 = cc[7];
901 cmptrk.c42 = cc[8];
902 cmptrk.c43 = cc[13];
903 cmptrk.c44 = cc[9];
904
905 // fill tree
906 cmptrktree->Fill();
907
908 } // loop on tracks at TPC 1st hit
909
910 cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
911 cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
912 cerr<<"+++"<<endl;
913
914 // Write tree to file
915 TFile *outfile = new TFile(covmatName,"recreate");
916 cmptrktree->Write();
917 outfile->Close();
918
919
920 // Write efficiencies to file
921 FILE *effFile = fopen(tpceffName,"w");
922 //fprintf(effFile,"%d\n",kalEntries);
923 for(Int_t j=0; j<27; j++) {
924 if(geaPi[j]>=5) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
925 if(geaKa[j]>=5) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
926 if(geaPr[j]>=5) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
927 if(geaEl[j]>=5) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
928 if(geaMu[j]>=5) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
929 fprintf(effFile,"%f %f %f %f %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
930 }
6eb67451 931
6eb69b9a 932 for(Int_t j=0; j<27; j++) {
933 fprintf(effFile,"%d %d %d %d %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
934 }
935 for(Int_t j=0; j<27; j++) {
936 fprintf(effFile,"%d %d %d %d %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
937 }
938 fclose(effFile);
6eb67451 939
6eb69b9a 940 // delete AliRun object
941 delete gAlice; gAlice=0;
942
943 // close all input files
944 kalFile->Close();
945 geaFile->Close();
946 galiceFile->Close();
6eb67451 947
6eb69b9a 948 return;
6eb67451 949}
6eb69b9a 950//-----------------------------------------------------------------------------
951void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
952//-----------------------------------------------------------------------------
953// This function assigns the track a dE/dx and makes a rough PID
954//-----------------------------------------------------------------------------
6eb67451 955
6eb69b9a 956 Double_t mean = fdEdxMean->GetValueAt(pt,eta);
957 Double_t rms = fdEdxRMS->GetValueAt(pt,eta);
958
959 Double_t dEdx = gRandom->Gaus(mean,rms);
960
961 fTrack.SetdEdx(dEdx);
962
963 AliTPCtrackParam t(fTrack);
964
965 //Very rough PID
966 Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
967
968 if (p<0.6) {
969 if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
970 t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
6eb67451 971 }
6eb69b9a 972 if (dEdx < 39.+ 12./p/p) {
973 t.AssignMass(0.49368); new(&fTrack) AliTPCtrack(t); return;
974 }
975 t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
976 }
6eb67451 977
6eb69b9a 978 if (p<1.2) {
979 if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
980 t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
6eb67451 981 }
6eb69b9a 982 t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
6eb67451 983 }
984
6eb69b9a 985 t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
6eb67451 986}
6eb69b9a 987//-----------------------------------------------------------------------------
988void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
989//-----------------------------------------------------------------------------
990// This function deals with covariance matrix and smearing
991//-----------------------------------------------------------------------------
992
993 COVMATRIX covmat;
994 Double_t p,cosl;
995 Double_t trkKine[2],trkRegPar[4];
996 Double_t xref,alpha,xx[5],xxsm[5],cc[15];
997 Int_t treeEntries;
998
999 fCovTree->SetBranchAddress("matrix",&covmat);
1000
1001 // get random entry from the tree
1002 treeEntries = (Int_t)fCovTree->GetEntries();
1003 fCovTree->GetEvent(gRandom->Integer(treeEntries));
1004
1005 // get P and Cosl from track
1006 cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1007 p = 1./TMath::Abs(fTrack.Get1Pt())/cosl;
1008
1009 trkKine[0] = p;
1010 trkKine[1] = cosl;
1011
1012 // get covariance matrix from regularized matrix
1013 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(0,l);
1014 cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1015 cc[1] = covmat.c10;
1016 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(1,l);
1017 cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1018 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(2,l);
1019 cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1020 cc[4] = covmat.c21;
1021 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(3,l);
1022 cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1023 cc[6] = covmat.c30;
1024 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(4,l);
1025 cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1026 cc[8] = covmat.c32;
1027 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(5,l);
1028 cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1029 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(6,l);
1030 cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1031 cc[11]= covmat.c41;
1032 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(7,l);
1033 cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1034 cc[13]= covmat.c43;
1035 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(8,l);
1036 cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1037
1038 TMatrixD covMatSmear(5,5);
1039
1040 covMatSmear = GetSmearingMatrix(cc,pt,eta);
1041
1042 // get track original parameters
1043 xref =fTrack.GetX();
1044 alpha=fTrack.GetAlpha();
1045 xx[0]=fTrack.GetY();
1046 xx[1]=fTrack.GetZ();
1047 xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp();
1048 xx[3]=fTrack.GetTgl();
1049 xx[4]=fTrack.GetC();
1050
1051 // use smearing matrix to smear the original parameters
1052 SmearTrack(xx,xxsm,covMatSmear);
1053
1054 AliTPCtrack track(0,xxsm,cc,xref,alpha);
1055 new(&fTrack) AliTPCtrack(track);
1056
1057 return;
1058}
1059//-----------------------------------------------------------------------------
1060void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg=211) {
1061//-----------------------------------------------------------------------------
1062// This function draws the TPC efficiencies in the [pT,eta] bins
1063//-----------------------------------------------------------------------------
1064
1065 ReadEffs(inName);
1066 SetParticle(pdg);
1067
1068 const Int_t n = fEff->GetPointsPt();
1069 Double_t effsA[n],effsB[n],effsC[n],pt[n];
1070 fEff->GetArrayPt(pt);
1071 for(Int_t i=0;i<n;i++) {
1072 effsA[i] = fEff->GetParam(i,0);
1073 effsB[i] = fEff->GetParam(i,1);
1074 effsC[i] = fEff->GetParam(i,2);
1075 }
1076
1077 TGraph *grA = new TGraph(n,pt,effsA);
1078 TGraph *grB = new TGraph(n,pt,effsB);
1079 TGraph *grC = new TGraph(n,pt,effsC);
1080
1081 grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1082 grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1083 grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1084
1085 TString title("Distribution of the TPC efficiencies");
1086 switch (pdg) {
1087 case 211:
1088 title.Prepend("PIONS - ");
1089 break;
1090 case 321:
1091 title.Prepend("KAONS - ");
1092 break;
1093 case 2212:
1094 title.Prepend("PROTONS - ");
1095 break;
1096 case 11:
1097 title.Prepend("ELECTRONS - ");
1098 break;
1099 case 13:
1100 title.Prepend("MUONS - ");
1101 break;
1102 }
1103
1104
1105 gStyle->SetOptStat(0);
1106 TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1107 c->SetLogx();
1108 c->SetGridx();
1109 c->SetGridy();
1110
1111 TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1112 frame->SetXTitle("p_{T} [GeV/c]");
1113 frame->Draw();
1114 grA->Draw("P");
1115 grB->Draw("P");
1116 grC->Draw("P");
1117
1118 TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1119 leg->AddEntry(grA,"|#eta|<0.3","p");
1120 leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1121 leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1122 leg->SetFillColor(0);
1123 leg->Draw();
1124
1125
1126 return;
1127}
1128//-----------------------------------------------------------------------------
1129void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg=211,
1130 Int_t par=0) {
1131//-----------------------------------------------------------------------------
1132// This function draws the pulls in the [pT,eta] bins
1133//-----------------------------------------------------------------------------
1134
1135 ReadPulls(inName);
1136 SetParticle(pdg);
1137
1138 const Int_t n = (fPulls+par)->GetPointsPt();
1139 Double_t pullsA[n],pullsB[n],pullsC[n],pt[n];
1140 (fPulls+par)->GetArrayPt(pt);
1141 for(Int_t i=0;i<n;i++) {
1142 pullsA[i] = (fPulls+par)->GetParam(i,0);
1143 pullsB[i] = (fPulls+par)->GetParam(i,1);
1144 pullsC[i] = (fPulls+par)->GetParam(i,2);
1145 }
1146
1147 TGraph *grA = new TGraph(n,pt,pullsA);
1148 TGraph *grB = new TGraph(n,pt,pullsB);
1149 TGraph *grC = new TGraph(n,pt,pullsC);
1150
1151 grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1152 grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1153 grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1154
1155 TString title("Distribution of the pulls: ");
1156 switch (pdg) {
1157 case 211:
1158 title.Prepend("PIONS - ");
1159 break;
1160 case 321:
1161 title.Prepend("KAONS - ");
1162 break;
1163 case 11:
1164 title.Prepend("ELECTRONS - ");
1165 break;
1166 }
1167 switch (par) {
1168 case 0:
1169 title.Append("y");
1170 break;
1171 case 1:
1172 title.Append("z");
1173 break;
1174 case 2:
1175 title.Append(" #eta");
1176 break;
1177 case 3:
1178 title.Append("tg #lambda");
1179 break;
1180 case 4:
1181 title.Append("C");
1182 break;
1183 }
1184
1185 gStyle->SetOptStat(0);
1186 TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1187 c->SetLogx();
1188 c->SetGridx();
1189 c->SetGridy();
1190
1191 TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1192 frame->SetXTitle("p_{T} [GeV/c]");
1193 frame->Draw();
1194 grA->Draw("P");
1195 grB->Draw("P");
1196 grC->Draw("P");
1197
1198 TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1199 leg->AddEntry(grA,"|#eta|<0.3","p");
1200 leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1201 leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1202 leg->SetFillColor(0);
1203 leg->Draw();
1204
6eb67451 1205
6eb69b9a 1206 return;
1207}
1208//-----------------------------------------------------------------------------
e130146c 1209Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
6eb69b9a 1210//-----------------------------------------------------------------------------
e130146c 1211// This function tells bin number in a grid (pT,eta)
6eb69b9a 1212//-----------------------------------------------------------------------------
6eb67451 1213 if(TMath::Abs(eta)<0.3) {
1214 if(pt<0.3) return 0;
1215 if(pt>=0.3 && pt<0.5) return 3;
1216 if(pt>=0.5 && pt<1.) return 6;
1217 if(pt>=1. && pt<1.5) return 9;
1218 if(pt>=1.5 && pt<3.) return 12;
1219 if(pt>=3. && pt<5.) return 15;
1220 if(pt>=5. && pt<7.) return 18;
1221 if(pt>=7. && pt<15.) return 21;
1222 if(pt>=15.) return 24;
1223 }
1224 if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
1225 if(pt<0.3) return 1;
1226 if(pt>=0.3 && pt<0.5) return 4;
1227 if(pt>=0.5 && pt<1.) return 7;
1228 if(pt>=1. && pt<1.5) return 10;
1229 if(pt>=1.5 && pt<3.) return 13;
1230 if(pt>=3. && pt<5.) return 16;
1231 if(pt>=5. && pt<7.) return 19;
1232 if(pt>=7. && pt<15.) return 22;
1233 if(pt>=15.) return 25;
1234 }
1235 if(TMath::Abs(eta)>=0.6) {
1236 if(pt<0.3) return 2;
1237 if(pt>=0.3 && pt<0.5) return 5;
1238 if(pt>=0.5 && pt<1.) return 8;
1239 if(pt>=1. && pt<1.5) return 11;
1240 if(pt>=1.5 && pt<3.) return 14;
1241 if(pt>=3. && pt<5.) return 17;
1242 if(pt>=5. && pt<7.) return 20;
1243 if(pt>=7. && pt<15.) return 23;
1244 if(pt>=15.) return 26;
1245 }
1246
1247 return -1;
1248
1249}
6eb69b9a 1250//-----------------------------------------------------------------------------
1251TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1252 Double_t eta) const {
1253//-----------------------------------------------------------------------------
e130146c 1254// This function stretches the covariance matrix according to the pulls
6eb69b9a 1255//-----------------------------------------------------------------------------
1256
e130146c 1257 TMatrixD covMat(5,5);
1258
1259 covMat(0,0)=cc[0];
1260 covMat(1,0)=cc[1]; covMat(0,1)=covMat(1,0);
1261 covMat(1,1)=cc[2];
1262 covMat(2,0)=cc[3]; covMat(0,2)=covMat(2,0);
1263 covMat(2,1)=cc[4]; covMat(1,2)=covMat(2,1);
1264 covMat(2,2)=cc[5];
1265 covMat(3,0)=cc[6]; covMat(0,3)=covMat(3,0);
1266 covMat(3,1)=cc[7]; covMat(1,3)=covMat(3,1);
1267 covMat(3,2)=cc[8]; covMat(2,3)=covMat(3,2);
1268 covMat(3,3)=cc[9];
1269 covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1270 covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1271 covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1272 covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1273 covMat(4,4)=cc[14];
1274
1275 TMatrixD stretchMat(5,5);
6eb67451 1276 for(Int_t k=0;k<5;k++) {
1277 for(Int_t l=0;l<5;l++) {
e130146c 1278 stretchMat(k,l)=0.;
6eb67451 1279 }
1280 }
1281
6eb69b9a 1282 for(Int_t i=0;i<5;i++) stretchMat(i,i) =
1283 TMath::Sqrt((fPulls+i)->GetValueAt(pt,eta));
1284
1285 TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1286 TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1287
1288 return covMatSmear;
1289}
1290//-----------------------------------------------------------------------------
1291void AliTPCtrackerParam::InitializeKineGrid(Option_t* which,Option_t* how) {
1292//-----------------------------------------------------------------------------
1293// This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1294// which = "DB" -> initialize fDBgrid... members
1295// "eff" -> initialize fEff... members
1296// "pulls" -> initialize fPulls... members
1297// "dEdx" -> initialize fdEdx... members
1298// how = "points" -> initialize with the points of the grid
1299// "bins" -> initialize with the bins of the grid
1300//-----------------------------------------------------------------------------
1301
1302 const char *points = strstr(how,"points");
1303 const char *bins = strstr(how,"bins");
1304 const char *DB = strstr(which,"DB");
1305 const char *eff = strstr(which,"eff");
1306 const char *pulls = strstr(which,"pulls");
1307 const char *dEdx = strstr(which,"dEdx");
1308
1309
1310 Int_t nEta=0,nPtPi=0,nPtKa=0,nPtPr=0,nPtEl=0,nPtMu=0;
1311
1312 Double_t etaPoints[2] = {0.3,0.6};
1313 Double_t etaBins[3] = {0.15,0.45,0.75};
1314 Double_t ptPoints8[8] = {0.3,0.5,1.,1.5,3.,5.,7.,15.};
1315 Double_t ptPoints5[5] = {0.3,0.5,1.,1.5,5.};
1316 Double_t ptBins9[9] = {0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
1317 Double_t ptBins6[6] = {0.244,0.390,0.676,1.190,2.36,10.};
1318
1319 Double_t *eta=0,*ptPi=0,*ptKa=0,*ptPr=0,*ptEl=0,*ptMu=0;
1320
1321 if(points) {
1322 nEta = 2;
1323 nPtPi = 8;
1324 nPtKa = 5;
1325 nPtPr = 4;
1326 nPtEl = 8;
1327 nPtMu = 2;
1328 eta = etaPoints;
1329 ptPi = ptPoints8;
1330 ptKa = ptPoints5;
1331 ptPr = ptPoints8;
1332 ptEl = ptPoints8;
1333 ptMu = ptPoints8;
1334 }
1335 if(bins) {
1336 nEta = 3;
1337 nPtPi = 9;
1338 nPtKa = 6;
1339 nPtPr = 5;
1340 nPtEl = 9;
1341 nPtMu = 3;
1342 eta = etaBins;
1343 ptPi = ptBins9;
1344 ptKa = ptBins6;
1345 ptPr = ptBins9;
1346 ptEl = ptBins9;
1347 ptMu = ptBins9;
1348 }
1349
1350 AliTPCkineGrid *dummy=0;
1351
1352 if(DB) {
1353 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1354 new(&fDBgridPi) AliTPCkineGrid(*dummy);
1355 delete dummy;
1356 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1357 new(&fDBgridKa) AliTPCkineGrid(*dummy);
1358 delete dummy;
1359 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1360 new(&fDBgridEl) AliTPCkineGrid(*dummy);
1361 delete dummy;
1362 }
1363 if(eff) {
1364 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1365 new(&fEffPi) AliTPCkineGrid(*dummy);
1366 delete dummy;
1367 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1368 new(&fEffKa) AliTPCkineGrid(*dummy);
1369 delete dummy;
1370 dummy = new AliTPCkineGrid(nPtPr,nEta,ptPr,eta);
1371 new(&fEffPr) AliTPCkineGrid(*dummy);
1372 delete dummy;
1373 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1374 new(&fEffEl) AliTPCkineGrid(*dummy);
1375 delete dummy;
1376 dummy = new AliTPCkineGrid(nPtMu,nEta,ptMu,eta);
1377 new(&fEffMu) AliTPCkineGrid(*dummy);
1378 delete dummy;
1379 }
1380 if(pulls) {
1381 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1382 for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1383 delete dummy;
1384 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1385 for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1386 delete dummy;
1387 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1388 for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1389 delete dummy;
1390 }
1391 if(dEdx) {
1392 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1393 new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1394 new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1395 delete dummy;
1396 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1397 new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1398 new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1399 delete dummy;
1400 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1401 new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1402 new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1403 delete dummy;
1404 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1405 new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1406 new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1407 delete dummy;
1408 }
1409
1410 return;
1411}
1412//-----------------------------------------------------------------------------
1413void AliTPCtrackerParam::MakeDataBase() {
1414//-----------------------------------------------------------------------------
1415// This function creates the DB file and store in it:
1416// - TPC Efficiencies for pi,ka,pr,el,mu (AliTPCkineGrid class)
1417// - Pulls for pi,ka,el (AliTPCkineGrid class)
1418// - Regularization parameters for pi,ka,el (TMatrixD class)
1419// - dE/dx parameterization for pi,ka,pr,el (AliTPCkineGrid class)
1420// - Regularized cov. matrices for pi,ka,el (COVMATRIX structure)
1421//-----------------------------------------------------------------------------
1422
1423 // define some file names
1424 Char_t *effFile ="CovMatrixDB_partial.root";
1425 Char_t *pullsFile="CovMatrixDB_partial.root";
1426 Char_t *regPiFile="CovMatrixDB_partial.root";
1427 Char_t *regKaFile="CovMatrixDB_partial.root";
1428 Char_t *regElFile="CovMatrixDB_partial.root";
1429 Char_t *dEdxPiFile="dEdxPi.root";
1430 Char_t *dEdxKaFile="dEdxKa.root";
1431 Char_t *dEdxPrFile="dEdxPr.root";
1432 Char_t *dEdxElFile="dEdxEl.root";
1433 Char_t *cmFile1 ="CovMatrix_AllEvts_1.root";
1434 Char_t *cmFile2 ="CovMatrix_AllEvts_2.root";
1435 Char_t *cmFile3 ="CovMatrix_AllEvts_3.root";
1436
1437 // store the effieciencies
1438 ReadEffs(effFile);
1439 WriteEffs(fDBfileName.Data());
1440
1441 // store the pulls
1442 ReadPulls(pullsFile);
1443 WritePulls(fDBfileName.Data());
1444
1445 //store the regularization parameters
1446 ReadRegParams(regPiFile,211);
1447 WriteRegParams(fDBfileName.Data(),211);
1448 ReadRegParams(regKaFile,321);
1449 WriteRegParams(fDBfileName.Data(),321);
1450 ReadRegParams(regElFile,11);
1451 WriteRegParams(fDBfileName.Data(),11);
1452
1453 // store the dEdx parameters
1454 ReaddEdx(dEdxPiFile,211);
1455 WritedEdx(fDBfileName.Data(),211);
1456 ReaddEdx(dEdxKaFile,321);
1457 WritedEdx(fDBfileName.Data(),321);
1458 ReaddEdx(dEdxPrFile,2212);
1459 WritedEdx(fDBfileName.Data(),2212);
1460 ReaddEdx(dEdxElFile,11);
1461 WritedEdx(fDBfileName.Data(),11);
1462
1463
1464 //
1465 // store the regularized covariance matrices
1466 //
1467 InitializeKineGrid("DB","points");
1468
1469 const Int_t nBinsPi = fDBgridPi.GetTotBins();
1470 const Int_t nBinsKa = fDBgridKa.GetTotBins();
1471 const Int_t nBinsEl = fDBgridEl.GetTotBins();
1472
1473
1474 // create the trees for cov. matrices
1475 // trees for pions
1476 TTree *CovTreePi_ = NULL;
1477 CovTreePi_ = new TTree[nBinsPi];
1478 // trees for kaons
1479 TTree *CovTreeKa_ = NULL;
1480 CovTreeKa_ = new TTree[nBinsKa];
1481 // trees for electrons
1482 TTree *CovTreeEl_ = NULL;
1483 CovTreeEl_ = new TTree[nBinsEl];
1484
1485 Char_t hname[100], htitle[100];
1486 COVMATRIX covmat;
1487
1488
1489 for(Int_t i=0; i<nBinsPi; i++) {
1490 sprintf(hname,"CovTreePi_bin%d",i);
1491 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1492 CovTreePi_[i].SetName(hname); CovTreePi_[i].SetTitle(htitle);
1493 CovTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1494 }
1495 for(Int_t i=0; i<nBinsKa; i++) {
1496 sprintf(hname,"CovTreeKa_bin%d",i);
1497 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1498 CovTreeKa_[i].SetName(hname); CovTreeKa_[i].SetTitle(htitle);
1499 CovTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1500 }
1501 for(Int_t i=0; i<nBinsEl; i++) {
1502 sprintf(hname,"CovTreeEl_bin%d",i);
1503 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1504 CovTreeEl_[i].SetName(hname); CovTreeEl_[i].SetTitle(htitle);
1505 CovTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1506 }
6eb67451 1507
6eb69b9a 1508 // create the chain with the compared tracks
1509 TChain cmptrkchain("cmptrktree");
1510 cmptrkchain.Add(cmFile1);
1511 cmptrkchain.Add(cmFile2);
1512 cmptrkchain.Add(cmFile3);
1513
1514 COMPTRACK cmptrk;
1515 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
1516 Int_t entries = (Int_t)cmptrkchain.GetEntries();
1517 cerr<<" Number of entries: "<<entries<<endl;
1518
1519 Int_t trkPdg,trkBin;
1520 Double_t trkKine[2],trkRegPar[4];
1521 Int_t nPerBinPi[nBinsPi]; for(Int_t k=0;k<nBinsPi;k++) nPerBinPi[k]=0;
1522 Int_t nPerBinKa[nBinsKa]; for(Int_t k=0;k<nBinsKa;k++) nPerBinKa[k]=0;
1523 Int_t nPerBinEl[nBinsEl]; for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
1524
1525 // loop on chain entries
1526 for(Int_t l=0; l<entries; l++) {
1527 if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1528 // get the event
1529 cmptrkchain.GetEvent(l);
1530 // get the pdg code
1531 trkPdg = TMath::Abs(cmptrk.pdg);
1532 // use only pions, kaons, electrons
1533 if(trkPdg!=211 && trkPdg!=321 && trkPdg!=11) continue;
1534 SetParticle(trkPdg);
1535 trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1536 //cerr<<cmptrk.pt<<" "<<cmptrk.eta<<" "<<trkBin<<endl;
1537
1538 if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1539 if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1540 if(trkPdg==11 && nPerBinEl[trkBin]>=5000) continue;
1541
1542 trkKine[0] = cmptrk.p;
1543 trkKine[1] = cmptrk.cosl;
1544
1545 // get regularized covariance matrix
1546 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(0,k);
1547 covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1548 covmat.c10 = cmptrk.c10;
1549 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(1,k);
1550 covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1551 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(2,k);
1552 covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1553 covmat.c21 = cmptrk.c21;
1554 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(3,k);
1555 covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
1556 covmat.c30 = cmptrk.c30;
1557 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(4,k);
1558 covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
1559 covmat.c32 = cmptrk.c32;
1560 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(5,k);
1561 covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
1562 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(6,k);
1563 covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
1564 covmat.c41 = cmptrk.c41;
1565 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(7,k);
1566 covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
1567 covmat.c43 = cmptrk.c43;
1568 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(8,k);
1569 covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
1570
1571 // fill the tree
1572 switch (trkPdg) {
1573 case 211: // pions
1574 CovTreePi_[trkBin].Fill();
1575 nPerBinPi[trkBin]++;
1576 break;
1577 case 321: // kaons
1578 CovTreeKa_[trkBin].Fill();
1579 nPerBinKa[trkBin]++;
1580 break;
1581 case 11: // electrons
1582 CovTreeEl_[trkBin].Fill();
1583 nPerBinEl[trkBin]++;
1584 break;
1585 }
1586 } // loop on chain entries
1587
1588 // store all trees the DB file
1589 TFile *DBfile = new TFile(fDBfileName.Data(),"update");
1590 DBfile->mkdir("CovMatrices");
1591 gDirectory->cd("/CovMatrices");
1592 gDirectory->mkdir("Pions");
1593 gDirectory->mkdir("Kaons");
1594 gDirectory->mkdir("Electrons");
1595 // store pions
1596 gDirectory->cd("/CovMatrices/Pions");
1597 fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
1598 for(Int_t i=0;i<nBinsPi;i++) CovTreePi_[i].Write();
1599 // store kaons
1600 gDirectory->cd("/CovMatrices/Kaons");
1601 fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
1602 for(Int_t i=0;i<nBinsKa;i++) CovTreeKa_[i].Write();
1603 // store electrons
1604 gDirectory->cd("/CovMatrices/Electrons");
1605 fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
1606 for(Int_t i=0;i<nBinsEl;i++) CovTreeEl_[i].Write();
1607
1608 DBfile->Close();
1609
1610 return;
1611}
1612//-----------------------------------------------------------------------------
1613void AliTPCtrackerParam::MergeEvents(Int_t evFirst=1,Int_t evLast=1) {
1614//-----------------------------------------------------------------------------
1615// This function: 1) merges the files from track comparison
1616// (beware: better no more than 100 events per file)
1617// 2) computes the average TPC efficiencies
1618//-----------------------------------------------------------------------------
1619
1620 Char_t *outName="TPCeff.root";
1621
1622 Int_t nCol;
1623 Float_t effPi,effKa,effPr,effEl,effMu;
1624 Float_t avEffPi[27],avEffKa[27],avEffPr[27],avEffEl[27],avEffMu[27];
1625 Int_t evtsPi[27],evtsKa[27],evtsPr[27],evtsEl[27],evtsMu[27];
1626 for(Int_t j=0; j<27; j++) {
1627 avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
1628 evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
1629 }
6eb67451 1630
6eb69b9a 1631 // create the chain for the tree of compared tracks
1632 TChain ch1("cmptrktree");
1633 TChain ch2("cmptrktree");
1634 TChain ch3("cmptrktree");
6eb67451 1635
6eb67451 1636
6eb69b9a 1637 for(Int_t j=evFirst; j<=evLast; j++) {
1638 cerr<<"Processing event "<<j<<endl;
6eb67451 1639
6eb69b9a 1640 TString covName("CovMatrix.");
1641 TString effName("TPCeff.");
1642 covName+=j;
1643 effName+=j;
1644 covName.Append(".root");
1645 effName.Append(".dat");
6eb67451 1646
6eb69b9a 1647 if(gSystem->AccessPathName(covName.Data(),kFileExists) ||
1648 gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
1649
1650 if(j<=100) ch1.Add(covName.Data());
1651 if(j>100 && j<=200) ch2.Add(covName.Data());
1652 if(j>200) ch3.Add(covName.Data());
1653
1654
1655 FILE *effIn = fopen(effName.Data(),"r");
1656 for(Int_t k=0; k<27; k++) {
1657 nCol = fscanf(effIn,"%f %f %f %f %f",&effPi,&effKa,&effPr,&effEl,&effMu);
1658 if(effPi>=0.) {avEffPi[k]+=effPi; evtsPi[k]++;}
1659 if(effKa>=0.) {avEffKa[k]+=effKa; evtsKa[k]++;}
1660 if(effPr>=0.) {avEffPr[k]+=effPr; evtsPr[k]++;}
1661 if(effEl>=0.) {avEffEl[k]+=effEl; evtsEl[k]++;}
1662 if(effMu>=0.) {avEffMu[k]+=effMu; evtsMu[k]++;}
1663 }
1664 fclose(effIn);
1665
1666 } // loop on events
1667
1668 // merge chain in one file
1669 TFile *covOut=0;
1670 covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
1671 ch1.Merge(covOut,1000000000);
1672 covOut->Close();
1673 delete covOut;
1674 covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
1675 ch2.Merge(covOut,1000000000);
1676 covOut->Close();
1677 delete covOut;
1678 covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
1679 ch3.Merge(covOut,1000000000);
1680 covOut->Close();
1681 delete covOut;
1682
1683
1684 // efficiencies
1685
1686 InitializeKineGrid("eff","bins");
1687
1688 // pions
1689 for(Int_t j=0; j<27; j++) {
1690 if(evtsPi[j]==0) evtsPi[j]++;
1691 fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
1692 }
1693
1694 // kaons
1695 Int_t j=0;
1696 for(Int_t k=0; k<27; k++) {
1697 if(k>14 && k!=21 && k!=22 && k!=23) continue;
1698 if(evtsKa[k]==0) evtsKa[k]++;
1699 fEffKa.SetParam(j,(Double_t)avEffKa[k]/evtsKa[k]);
1700 j++;
1701 }
1702
1703 // protons
1704 for(Int_t j=0; j<15; j++) {
1705 if(evtsPr[j]==0) evtsPr[j]++;
1706 fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
1707 }
1708
1709 // electrons
1710 for(Int_t j=0; j<27; j++) {
1711 if(evtsEl[j]==0) evtsEl[j]++;
1712 fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
1713 }
1714
1715 // muons
1716 for(Int_t j=0; j<9; j++) {
1717 if(evtsMu[j]==0) evtsMu[j]++;
1718 fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
1719 }
1720
1721 // write efficiencies to a file
1722 WriteEffs(outName);
1723
1724 return;
1725}
1726//-----------------------------------------------------------------------------
1727Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
1728//-----------------------------------------------------------------------------
1729// This function reads all parameters from the DB
1730//-----------------------------------------------------------------------------
1731
1732 if(!ReadEffs(inName)) return 0;
1733 if(!ReadPulls(inName)) return 0;
1734 if(!ReadRegParams(inName,211)) return 0;
1735 if(!ReadRegParams(inName,321)) return 0;
1736 if(!ReadRegParams(inName,11)) return 0;
1737 if(!ReaddEdx(inName,211)) return 0;
1738 if(!ReaddEdx(inName,321)) return 0;
1739 if(!ReaddEdx(inName,2212)) return 0;
1740 if(!ReaddEdx(inName,11)) return 0;
1741 if(!ReadDBgrid(inName)) return 0;
1742
1743 return 1;
6eb67451 1744}
6eb69b9a 1745//-----------------------------------------------------------------------------
1746Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
1747//-----------------------------------------------------------------------------
1748// This function reads the dEdx parameters from the DB
1749//-----------------------------------------------------------------------------
1750
1751 if(gSystem->AccessPathName(inName,kFileExists)) {
1752 cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
1753 return 0; }
1754
1755 TFile *inFile = TFile::Open(inName);
1756 switch (pdg) {
1757 case 211:
1758 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
1759 fdEdxMeanPi.~AliTPCkineGrid();
1760 new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
1761 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
1762 fdEdxRMSPi.~AliTPCkineGrid();
1763 new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
1764 break;
1765 case 321:
1766 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
1767 fdEdxMeanKa.~AliTPCkineGrid();
1768 new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
1769 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
1770 fdEdxRMSKa.~AliTPCkineGrid();
1771 new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
1772 break;
1773 case 2212:
1774 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
1775 fdEdxMeanPr.~AliTPCkineGrid();
1776 new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
1777 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
1778 fdEdxRMSPr.~AliTPCkineGrid();
1779 new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
1780 break;
1781 case 11:
1782 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
1783 fdEdxMeanEl.~AliTPCkineGrid();
1784 new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
1785 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
1786 fdEdxRMSEl.~AliTPCkineGrid();
1787 new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
1788 break;
1789 }
1790 inFile->Close();
6eb67451 1791
6eb69b9a 1792 return 1;
1793}
1794//-----------------------------------------------------------------------------
1795Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
1796//-----------------------------------------------------------------------------
1797// This function reads the kine grid from the DB
1798//-----------------------------------------------------------------------------
1799
1800 if(gSystem->AccessPathName(inName,kFileExists)) {
1801 cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n";
1802 return 0; }
1803
1804 TFile *inFile = TFile::Open(inName);
1805
1806 // first read the DB grid for the different particles
1807 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
1808 fDBgridPi.~AliTPCkineGrid();
1809 new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
1810 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
1811 fDBgridKa.~AliTPCkineGrid();
1812 new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
1813 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
1814 fDBgridEl.~AliTPCkineGrid();
1815 new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
1816
1817 inFile->Close();
1818
1819 return 1;
1820}
1821//-----------------------------------------------------------------------------
1822Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
1823//-----------------------------------------------------------------------------
1824// This function reads the TPC efficiencies from the DB
1825//-----------------------------------------------------------------------------
1826
1827 if(gSystem->AccessPathName(inName,kFileExists)) {
1828 cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n";
1829 return 0; }
1830
1831 TFile *inFile = TFile::Open(inName);
1832
1833 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
1834 fEffPi.~AliTPCkineGrid();
1835 new(&fEffPi) AliTPCkineGrid(*fEff);
1836 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
1837 fEffKa.~AliTPCkineGrid();
1838 new(&fEffKa) AliTPCkineGrid(*fEff);
1839 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
1840 fEffPr.~AliTPCkineGrid();
1841 new(&fEffPr) AliTPCkineGrid(*fEff);
1842 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
1843 fEffEl.~AliTPCkineGrid();
1844 new(&fEffEl) AliTPCkineGrid(*fEff);
1845 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
1846 fEffMu.~AliTPCkineGrid();
1847 new(&fEffMu) AliTPCkineGrid(*fEff);
1848
1849 inFile->Close();
1850
1851 return 1;
1852}
1853//-----------------------------------------------------------------------------
1854Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
1855//-----------------------------------------------------------------------------
1856// This function reads the pulls from the DB
1857//-----------------------------------------------------------------------------
1858
1859 if(gSystem->AccessPathName(inName,kFileExists)) {
1860 cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n";
1861 return 0; }
1862
1863 TFile *inFile = TFile::Open(inName);
1864
1865 for(Int_t i=0; i<5; i++) {
1866 TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
1867 TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
1868 TString el("/Pulls/Electrons/PullsEl_"); el+=i;
1869 fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
1870 fPullsPi[i].~AliTPCkineGrid();
1871 new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
1872 fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
1873 fPullsKa[i].~AliTPCkineGrid();
1874 new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
1875 fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
1876 fPullsEl[i].~AliTPCkineGrid();
1877 new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
1878 }
1879
1880 inFile->Close();
1881
1882 return 1;
1883}
1884//-----------------------------------------------------------------------------
1885Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
1886//-----------------------------------------------------------------------------
1887// This function reads the regularization parameters from the DB
1888//-----------------------------------------------------------------------------
1889
1890 if(gSystem->AccessPathName(inName,kFileExists)) {
1891 cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
1892 return 0; }
1893
1894 TFile *inFile = TFile::Open(inName);
1895 switch (pdg) {
1896 case 211:
1897 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
1898 new(&fRegParPi) TMatrixD(*fRegPar);
1899 break;
1900 case 321:
1901 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
1902 new(&fRegParKa) TMatrixD(*fRegPar);
1903 break;
1904 case 11:
1905 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
1906 new(&fRegParEl) TMatrixD(*fRegPar);
1907 break;
1908 }
1909 inFile->Close();
1910
1911 return 1;
1912}
1913//-----------------------------------------------------------------------------
1914void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
1915//-----------------------------------------------------------------------------
1916// This function regularizes the elements of the covariance matrix
1917// that show a momentum depence:
1918// c00, c11, c22, c33, c44, c20, c24, c40, c31
1919// The regularization is done separately for pions, kaons, electrons:
1920// give "Pion","Kaon" or "Electron" as first argument.
1921//-----------------------------------------------------------------------------
1922
1923 gStyle->SetOptStat(0);
1924 gStyle->SetOptFit(10001);
1925
1926 Int_t thisPdg=211;
1927 Char_t *part="Pions - ";
1928
1929 InitializeKineGrid("DB","points");
1930 SetParticle(pdg);
1931 const Int_t fitbins = fDBgrid->GetBinsPt();
1932 cerr<<" Fit bins: "<<fitbins<<endl;
1933
1934 switch (pdg) {
1935 case 211: // pions
1936 thisPdg=211;
1937 part="Pions - ";
1938 cerr<<" Processing pions ...\n";
1939 break;
1940 case 321: //kaons
1941 thisPdg=321;
1942 part="Kaons - ";
1943 cerr<<" Processing kaons ...\n";
1944 break;
1945 case 11: // electrons
1946 thisPdg= 11;
1947 part="Electrons - ";
1948 cerr<<" Processing electrons ...\n";
1949 break;
1950 }
1951
1952 // create the chain with all compared tracks
1953 TChain cmptrkchain("cmptrktree");
1954 cmptrkchain.Add("CovMatrix_AllEvts_1.root");
1955 cmptrkchain.Add("CovMatrix_AllEvts_2.root");
1956 cmptrkchain.Add("CovMatrix_AllEvts_3.root");
1957
1958 COMPTRACK cmptrk;
1959 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
1960 Int_t entries = (Int_t)cmptrkchain.GetEntries();
1961
1962
1963 Int_t pbin;
1964 Int_t n[fitbins];
1965 Int_t n00[fitbins],n11[fitbins],n20[fitbins],n22[fitbins],n31[fitbins],n33[fitbins],n40[fitbins],n42[fitbins],n44[fitbins];
1966 Double_t p[fitbins],ep[fitbins];
1967 Double_t mean00[fitbins],mean11[fitbins],mean20[fitbins],mean22[fitbins],mean31[fitbins],mean33[fitbins],mean40[fitbins],mean42[fitbins],mean44[fitbins];
1968 Double_t sigma00[fitbins],sigma11[fitbins],sigma20[fitbins],sigma22[fitbins],sigma31[fitbins],sigma33[fitbins],sigma40[fitbins],sigma42[fitbins],sigma44[fitbins];
1969 Double_t rmean[fitbins],rsigma[fitbins];
1970 Double_t fitpar[4];
1971
1972 for(Int_t l=0; l<fitbins; l++) {
1973 n[l]=1;
1974 n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
1975 p[l ]=ep[l]=0.;
1976 mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
1977 sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
1978 }
1979
1980 // loop on chain entries for mean
1981 for(Int_t l=0; l<entries; l++) {
1982 cmptrkchain.GetEvent(l);
1983 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
1984 pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
1985 n[pbin]++;
1986 p[pbin]+=cmptrk.p;
1987 mean00[pbin]+=cmptrk.c00;
1988 mean11[pbin]+=cmptrk.c11;
1989 mean20[pbin]+=cmptrk.c20;
1990 mean22[pbin]+=cmptrk.c22;
1991 mean31[pbin]+=cmptrk.c31;
1992 mean33[pbin]+=cmptrk.c33;
1993 mean40[pbin]+=cmptrk.c40;
1994 mean42[pbin]+=cmptrk.c42;
1995 mean44[pbin]+=cmptrk.c44;
1996 } // loop on chain entries
1997
1998 for(Int_t l=0; l<fitbins; l++) {
1999 p[l]/=n[l];
2000 mean00[l]/=n[l];
2001 mean11[l]/=n[l];
2002 mean20[l]/=n[l];
2003 mean22[l]/=n[l];
2004 mean31[l]/=n[l];
2005 mean33[l]/=n[l];
2006 mean40[l]/=n[l];
2007 mean42[l]/=n[l];
2008 mean44[l]/=n[l];
2009 }
2010
2011 // loop on chain entries for sigma
2012 for(Int_t l=0; l<entries; l++) {
2013 cmptrkchain.GetEvent(l);
2014 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2015 pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2016 if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2017 sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2018 if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2019 sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2020 if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2021 sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2022 if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2023 sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2024 if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2025 sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2026 if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2027 sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2028 if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2029 sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2030 if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2031 sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2032 if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2033 sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2034 } // loop on chain entries
2035
2036 for(Int_t l=0; l<fitbins; l++) {
2037 sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2038 sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2039 sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2040 sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2041 sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2042 sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2043 sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2044 sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2045 sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2046 }
2047
2048
2049 // Fit function
2050 TF1 *func = new TF1("FitRegFunc",FitRegFunc,0.23,50.,4);
2051 func->SetParNames("A_meas","A_scatt","B1","B2");
2052
2053 // line to draw on the plots
2054 TLine *lin = new TLine(-1,1,1.69,1);
2055 lin->SetLineStyle(2);
2056 lin->SetLineWidth(2);
2057
2058 // matrix used to store fit results
2059 TMatrixD fitRes(9,4);
2060
2061 // --- c00 ---
2062
2063 // create the canvas
2064 TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
2065 canv00->Divide(1,2);
2066 // create the graph for cov matrix
2067 TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
2068 TString title00("C(y,y)"); title00.Prepend(part);
2069 TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
2070 frame00->SetXTitle("p [GeV/c]");
2071 canv00->cd(1); gPad->SetLogx();
2072 frame00->Draw();
2073 gr00->Draw("P");
2074 // Sets initial values for parameters
2075 func->SetParameters(1.6e-3,1.9e-4,1.5,0.);
2076 // Fit points in range defined by function
2077 gr00->Fit("FitRegFunc","R,Q");
2078 func->GetParameters(fitpar);
2079 for(Int_t i=0; i<4; i++) fitRes(0,i)=fitpar[i];
2080 for(Int_t l=0; l<fitbins; l++) {
2081 rmean[l] = mean00[l]/FitRegFunc(&p[l],fitpar);
2082 rsigma[l] = sigma00[l]/FitRegFunc(&p[l],fitpar);
2083 }
2084 // create the graph the regularized cov. matrix
2085 TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2086 TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2087 regtitle00.Prepend(part);
2088 TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
2089 frame00reg->SetXTitle("p [GeV/c]");
2090 canv00->cd(2); gPad->SetLogx();
2091 frame00reg->Draw();
2092 gr00reg->Draw("P");
2093 lin->Draw("same");
2094
2095
2096 // --- c11 ---
2097
2098 // create the canvas
2099 TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
2100 canv11->Divide(1,2);
2101 // create the graph for cov matrix
2102 TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
2103 TString title11("C(z,z)"); title11.Prepend(part);
2104 TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
2105 frame11->SetXTitle("p [GeV/c]");
2106 canv11->cd(1); gPad->SetLogx();
2107 frame11->Draw();
2108 gr11->Draw("P");
2109 // Sets initial values for parameters
2110 func->SetParameters(1.2e-3,8.1e-4,1.,0.);
2111 // Fit points in range defined by function
2112 gr11->Fit("FitRegFunc","R,Q");
2113 func->GetParameters(fitpar);
2114 for(Int_t i=0; i<4; i++) fitRes(1,i)=fitpar[i];
2115 for(Int_t l=0; l<fitbins; l++) {
2116 rmean[l] = mean11[l]/FitRegFunc(&p[l],fitpar);
2117 rsigma[l] = sigma11[l]/FitRegFunc(&p[l],fitpar);
2118 }
2119 // create the graph the regularized cov. matrix
2120 TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2121 TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2122 regtitle11.Prepend(part);
2123 TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
2124 frame11reg->SetXTitle("p [GeV/c]");
2125 canv11->cd(2); gPad->SetLogx();
2126 frame11reg->Draw();
2127 gr11reg->Draw("P");
2128 lin->Draw("same");
2129
2130
2131 // --- c20 ---
2132
2133 // create the canvas
2134 TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
2135 canv20->Divide(1,2);
2136 // create the graph for cov matrix
2137 TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
2138 TString title20("C(#eta, y)"); title20.Prepend(part);
2139 TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
2140 frame20->SetXTitle("p [GeV/c]");
2141 canv20->cd(1); gPad->SetLogx();
2142 frame20->Draw();
2143 gr20->Draw("P");
2144 // Sets initial values for parameters
2145 func->SetParameters(7.3e-5,1.2e-5,1.5,0.);
2146 // Fit points in range defined by function
2147 gr20->Fit("FitRegFunc","R,Q");
2148 func->GetParameters(fitpar);
2149 for(Int_t i=0; i<4; i++) fitRes(2,i)=fitpar[i];
2150 for(Int_t l=0; l<fitbins; l++) {
2151 rmean[l] = mean20[l]/FitRegFunc(&p[l],fitpar);
2152 rsigma[l] = sigma20[l]/FitRegFunc(&p[l],fitpar);
2153 }
2154 // create the graph the regularized cov. matrix
2155 TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2156 TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2157 regtitle20.Prepend(part);
2158 TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
2159 frame20reg->SetXTitle("p [GeV/c]");
2160 canv20->cd(2); gPad->SetLogx();
2161 frame20reg->Draw();
2162 gr20reg->Draw("P");
2163 lin->Draw("same");
2164
2165
2166 // --- c22 ---
2167
2168 // create the canvas
2169 TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
2170 canv22->Divide(1,2);
2171 // create the graph for cov matrix
2172 TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
2173 TString title22("C(#eta, #eta)"); title22.Prepend(part);
2174 TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
2175 frame22->SetXTitle("p [GeV/c]");
2176 canv22->cd(1); gPad->SetLogx();
2177 frame22->Draw();
2178 gr22->Draw("P");
2179 // Sets initial values for parameters
2180 func->SetParameters(5.2e-6,1.1e-6,2.,1.);
2181 // Fit points in range defined by function
2182 gr22->Fit("FitRegFunc","R,Q");
2183 func->GetParameters(fitpar);
2184 for(Int_t i=0; i<4; i++) fitRes(3,i)=fitpar[i];
2185 for(Int_t l=0; l<fitbins; l++) {
2186 rmean[l] = mean22[l]/FitRegFunc(&p[l],fitpar);
2187 rsigma[l] = sigma22[l]/FitRegFunc(&p[l],fitpar);
2188 }
2189 // create the graph the regularized cov. matrix
2190 TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2191 TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2192 regtitle22.Prepend(part);
2193 TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
2194 frame22reg->SetXTitle("p [GeV/c]");
2195 canv22->cd(2); gPad->SetLogx();
2196 frame22reg->Draw();
2197 gr22reg->Draw("P");
2198 lin->Draw("same");
2199
2200
2201 // --- c31 ---
2202
2203 // create the canvas
2204 TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
2205 canv31->Divide(1,2);
2206 // create the graph for cov matrix
2207 TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
2208 TString title31("C(tg #lambda,z)"); title31.Prepend(part);
2209 TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
2210 frame31->SetXTitle("p [GeV/c]");
2211 canv31->cd(1); gPad->SetLogx();
2212 frame31->Draw();
2213 gr31->Draw("P");
2214 // Sets initial values for parameters
2215 func->SetParameters(-1.2e-5,-1.2e-5,1.5,3.);
2216 // Fit points in range defined by function
2217 gr31->Fit("FitRegFunc","R,Q");
2218 func->GetParameters(fitpar);
2219 for(Int_t i=0; i<4; i++) fitRes(4,i)=fitpar[i];
2220 for(Int_t l=0; l<fitbins; l++) {
2221 rmean[l] = mean31[l]/FitRegFunc(&p[l],fitpar);
2222 rsigma[l] = -sigma31[l]/FitRegFunc(&p[l],fitpar);
2223 }
2224 // create the graph the regularized cov. matrix
2225 TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2226 TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2227 regtitle31.Prepend(part);
2228 TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
2229 frame31reg->SetXTitle("p [GeV/c]");
2230 canv31->cd(2); gPad->SetLogx();
2231 frame31reg->Draw();
2232 gr31reg->Draw("P");
2233 lin->Draw("same");
2234
2235
2236 // --- c33 ---
2237
2238 // create the canvas
2239 TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
2240 canv33->Divide(1,2);
2241 // create the graph for cov matrix
2242 TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
2243 TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
2244 TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
2245 frame33->SetXTitle("p [GeV/c]");
2246 canv33->cd(1); gPad->SetLogx();
2247 frame33->Draw();
2248 gr33->Draw("P");
2249 // Sets initial values for parameters
2250 func->SetParameters(1.3e-7,4.6e-7,1.7,4.);
2251 // Fit points in range defined by function
2252 gr33->Fit("FitRegFunc","R,Q");
2253 func->GetParameters(fitpar);
2254 for(Int_t i=0; i<4; i++) fitRes(5,i)=fitpar[i];
2255 for(Int_t l=0; l<fitbins; l++) {
2256 rmean[l] = mean33[l]/FitRegFunc(&p[l],fitpar);
2257 rsigma[l] = sigma33[l]/FitRegFunc(&p[l],fitpar);
2258 }
2259 // create the graph the regularized cov. matrix
2260 TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2261 TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2262 regtitle33.Prepend(part);
2263 TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
2264 frame33reg->SetXTitle("p [GeV/c]");
2265 canv33->cd(2); gPad->SetLogx();
2266 frame33reg->Draw();
2267 gr33reg->Draw("P");
2268 lin->Draw("same");
2269
2270
2271 // --- c40 ---
2272
2273 // create the canvas
2274 TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
2275 canv40->Divide(1,2);
2276 // create the graph for cov matrix
2277 TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
2278 TString title40("C(C,y)"); title40.Prepend(part);
2279 TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
2280 frame40->SetXTitle("p [GeV/c]");
2281 canv40->cd(1); gPad->SetLogx();
2282 frame40->Draw();
2283 gr40->Draw("P");
2284 // Sets initial values for parameters
2285 func->SetParameters(4.e-7,4.4e-8,1.5,0.);
2286 // Fit points in range defined by function
2287 gr40->Fit("FitRegFunc","R,Q");
2288 func->GetParameters(fitpar);
2289 for(Int_t i=0; i<4; i++) fitRes(6,i)=fitpar[i];
2290 for(Int_t l=0; l<fitbins; l++) {
2291 rmean[l] = mean40[l]/FitRegFunc(&p[l],fitpar);
2292 rsigma[l] = sigma40[l]/FitRegFunc(&p[l],fitpar);
2293 }
2294 // create the graph the regularized cov. matrix
2295 TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2296 TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2297 regtitle40.Prepend(part);
2298 TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
2299 frame40reg->SetXTitle("p [GeV/c]");
2300 canv40->cd(2); gPad->SetLogx();
2301 frame40reg->Draw();
2302 gr40reg->Draw("P");
2303 lin->Draw("same");
2304
2305
2306 // --- c42 ---
2307
2308 // create the canvas
2309 TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
2310 canv42->Divide(1,2);
2311 // create the graph for cov matrix
2312 TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
2313 TString title42("C(C, #eta)"); title42.Prepend(part);
2314 TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
2315 frame42->SetXTitle("p [GeV/c]");
2316 canv42->cd(1); gPad->SetLogx();
2317 frame42->Draw();
2318 gr42->Draw("P");
2319 // Sets initial values for parameters
2320 func->SetParameters(3.e-8,8.2e-9,2.,1.);
2321 // Fit points in range defined by function
2322 gr42->Fit("FitRegFunc","R,Q");
2323 func->GetParameters(fitpar);
2324 for(Int_t i=0; i<4; i++) fitRes(7,i)=fitpar[i];
2325 for(Int_t l=0; l<fitbins; l++) {
2326 rmean[l] = mean42[l]/FitRegFunc(&p[l],fitpar);
2327 rsigma[l] = sigma42[l]/FitRegFunc(&p[l],fitpar);
2328 }
2329 // create the graph the regularized cov. matrix
2330 TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2331 TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2332 regtitle42.Prepend(part);
2333 TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
2334 frame42reg->SetXTitle("p [GeV/c]");
2335 canv42->cd(2); gPad->SetLogx();
2336 frame42reg->Draw();
2337 gr42reg->Draw("P");
2338 lin->Draw("same");
2339
2340
2341 // --- c44 ---
2342
2343 // create the canvas
2344 TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
2345 canv44->Divide(1,2);
2346 // create the graph for cov matrix
2347 TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
2348 TString title44("C(C,C)"); title44.Prepend(part);
2349 TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
2350 frame44->SetXTitle("p [GeV/c]");
2351 canv44->cd(1); gPad->SetLogx();
2352 frame44->Draw();
2353 gr44->Draw("P");
2354 // Sets initial values for parameters
2355 func->SetParameters(1.8e-10,5.8e-11,2.,3.);
2356 // Fit points in range defined by function
2357 gr44->Fit("FitRegFunc","R,Q");
2358 func->GetParameters(fitpar);
2359 for(Int_t i=0; i<4; i++) fitRes(8,i)=fitpar[i];
2360 for(Int_t l=0; l<fitbins; l++) {
2361 rmean[l] = mean44[l]/FitRegFunc(&p[l],fitpar);
2362 rsigma[l] = sigma44[l]/FitRegFunc(&p[l],fitpar);
2363 }
2364 // create the graph the regularized cov. matrix
2365 TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2366 TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2367 regtitle44.Prepend(part);
2368 TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
2369 frame44reg->SetXTitle("p [GeV/c]");
2370 canv44->cd(2); gPad->SetLogx();
2371 frame44reg->Draw();
2372 gr44reg->Draw("P");
2373 lin->Draw("same");
2374
2375
2376
2377
2378 switch (pdg) {
2379 case 211:
2380 new(&fRegParPi) TMatrixD(fitRes);
2381 break;
2382 case 321:
2383 new(&fRegParKa) TMatrixD(fitRes);
2384 break;
2385 case 11:
2386 new(&fRegParEl) TMatrixD(fitRes);
2387 break;
2388 }
2389
2390 // write fit parameters to file
2391 WriteRegParams(outName,pdg);
2392
2393 return;
2394}
2395//-----------------------------------------------------------------------------
2396Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
2397//-----------------------------------------------------------------------------
2398// This function makes a selection according to TPC tracking efficiency
2399//-----------------------------------------------------------------------------
2400
2401 Double_t eff=0.;
2402
2403 eff = fEff->GetValueAt(pt,eta);
2404
2405 if(gRandom->Rndm() < eff) return kTRUE;
2406
2407 return kFALSE;
2408}
2409//-----------------------------------------------------------------------------
2410void AliTPCtrackerParam::SetParticle(Int_t pdg) {
2411//-----------------------------------------------------------------------------
2412// This function set efficiencies, pulls, etc... for the particle
2413// specie of the current particle
2414//-----------------------------------------------------------------------------
2415
2416 switch (pdg) {
2417 case 211:
2418 fDBgrid = &fDBgridPi;
2419 fEff = &fEffPi;
2420 fPulls = fPullsPi;
2421 fRegPar = &fRegParPi;
2422 fdEdxMean = &fdEdxMeanPi;
2423 fdEdxRMS = &fdEdxRMSPi;
2424 break;
2425 case 321:
2426 fDBgrid = &fDBgridKa;
2427 fEff = &fEffKa;
2428 fPulls = fPullsKa;
2429 fRegPar = &fRegParKa;
2430 fdEdxMean = &fdEdxMeanKa;
2431 fdEdxRMS = &fdEdxRMSKa;
2432 break;
2433 case 2212:
2434 fDBgrid = &fDBgridPi;
2435 fEff = &fEffPr;
2436 fPulls = fPullsPi;
2437 fRegPar = &fRegParPi;
2438 fdEdxMean = &fdEdxMeanPr;
2439 fdEdxRMS = &fdEdxRMSPr;
2440 break;
2441 case 11:
2442 fDBgrid = &fDBgridEl;
2443 fEff = &fEffEl;
2444 fPulls = fPullsEl;
2445 fRegPar = &fRegParEl;
2446 fdEdxMean = &fdEdxMeanEl;
2447 fdEdxRMS = &fdEdxRMSEl;
2448 break;
2449 case 13:
2450 fDBgrid = &fDBgridPi;
2451 fEff = &fEffMu;
2452 fPulls = fPullsPi;
2453 fRegPar = &fRegParPi;
2454 fdEdxMean = &fdEdxMeanPi;
2455 fdEdxRMS = &fdEdxRMSPi;
2456 break;
2457 }
2458
2459 return;
2460}
2461//-----------------------------------------------------------------------------
2462void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
2463 const {
2464//-----------------------------------------------------------------------------
e130146c 2465// This function smears track parameters according to streched cov. matrix
6eb69b9a 2466//-----------------------------------------------------------------------------
70521312 2467 AliGausCorr *corgen = new AliGausCorr(cov,5);
2468 TArrayD corr(5);
2469 corgen->GetGaussN(corr);
2470 delete corgen;
2471 corgen = 0;
6eb67451 2472
2473 for(Int_t l=0;l<5;l++) {
2474 xxsm[l] = xx[l]+corr[l];
2475 }
2476
2477 return;
2478}
6eb69b9a 2479//-----------------------------------------------------------------------------
2480Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
2481//-----------------------------------------------------------------------------
2482// This function writes the dEdx parameters to the DB
2483//-----------------------------------------------------------------------------
2484
2485 Option_t *opt;
2486 Char_t *dirName="Pions";
2487 Char_t *meanName="dEdxMeanPi";
2488 Char_t *RMSName="dEdxRMSPi";
2489
2490 SetParticle(pdg);
2491
2492 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2493 } else { opt="update"; }
2494
2495 switch (pdg) {
2496 case 211:
2497 dirName="Pions";
2498 meanName="dEdxMeanPi";
2499 RMSName="dEdxRMSPi";
2500 break;
2501 case 321:
2502 dirName="Kaons";
2503 meanName="dEdxMeanKa";
2504 RMSName="dEdxRMSKa";
2505 break;
2506 case 2212:
2507 dirName="Protons";
2508 meanName="dEdxMeanPr";
2509 RMSName="dEdxRMSPr";
2510 break;
2511 case 11:
2512 dirName="Electrons";
2513 meanName="dEdxMeanEl";
2514 RMSName="dEdxRMSEl";
2515 break;
2516 }
6eb67451 2517
6eb69b9a 2518 TFile *outFile = new TFile(outName,opt);
2519 if(!gDirectory->cd("/dEdx")) {
2520 outFile->mkdir("dEdx");
2521 gDirectory->cd("/dEdx");
2522 }
2523 TDirectory *dir2 = gDirectory->mkdir(dirName);
2524 dir2->cd();
2525 fdEdxMean->SetName(meanName); fdEdxMean->Write();
2526 fdEdxRMS->SetName(RMSName); fdEdxRMS->Write();
6eb67451 2527
6eb69b9a 2528 outFile->Close();
2529 delete outFile;
6eb67451 2530
6eb67451 2531
6eb69b9a 2532 return 1;
2533}
2534//-----------------------------------------------------------------------------
2535Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
2536//-----------------------------------------------------------------------------
2537// This function writes the TPC efficiencies to the DB
2538//-----------------------------------------------------------------------------
2539
2540
2541 Option_t *opt;
2542
2543 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2544 } else { opt="update"; }
2545
2546 TFile *outFile = new TFile(outName,opt);
2547
2548 outFile->mkdir("Efficiencies");
2549 gDirectory->cd("/Efficiencies");
2550 gDirectory->mkdir("Pions");
2551 gDirectory->mkdir("Kaons");
2552 gDirectory->mkdir("Protons");
2553 gDirectory->mkdir("Electrons");
2554 gDirectory->mkdir("Muons");
2555
2556 gDirectory->cd("/Efficiencies/Pions");
2557 fEffPi.SetName("EffPi");
2558 fEffPi.Write();
2559 gDirectory->cd("/Efficiencies/Kaons");
2560 fEffKa.SetName("EffKa");
2561 fEffKa.Write();
2562 gDirectory->cd("/Efficiencies/Protons");
2563 fEffPr.SetName("EffPr");
2564 fEffPr.Write();
2565 gDirectory->cd("/Efficiencies/Electrons");
2566 fEffEl.SetName("EffEl");
2567 fEffEl.Write();
2568 gDirectory->cd("/Efficiencies/Muons");
2569 fEffMu.SetName("EffMu");
2570 fEffMu.Write();
2571
2572 outFile->Close();
2573
2574 delete outFile;
2575
2576 return 1;
2577}
2578//-----------------------------------------------------------------------------
2579Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
2580//-----------------------------------------------------------------------------
2581// This function writes the pulls to the DB
2582//-----------------------------------------------------------------------------
2583
2584 Option_t *opt;
2585
2586 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2587 } else { opt="update"; }
2588
2589 TFile *outFile = new TFile(outName,opt);
2590
2591 outFile->mkdir("Pulls");
2592 gDirectory->cd("/Pulls");
2593 gDirectory->mkdir("Pions");
2594 gDirectory->mkdir("Kaons");
2595 gDirectory->mkdir("Electrons");
2596
2597 for(Int_t i=0;i<5;i++) {
2598 TString pi("PullsPi_"); pi+=i;
2599 TString ka("PullsKa_"); ka+=i;
2600 TString el("PullsEl_"); el+=i;
2601 fPullsPi[i].SetName(pi.Data());
2602 fPullsKa[i].SetName(ka.Data());
2603 fPullsEl[i].SetName(el.Data());
2604 gDirectory->cd("/Pulls/Pions");
2605 fPullsPi[i].Write();
2606 gDirectory->cd("/Pulls/Kaons");
2607 fPullsKa[i].Write();
2608 gDirectory->cd("/Pulls/Electrons");
2609 fPullsEl[i].Write();
2610 }
2611 outFile->Close();
2612 delete outFile;
6eb67451 2613
6eb69b9a 2614 return 1;
2615}
2616//-----------------------------------------------------------------------------
2617Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
2618//-----------------------------------------------------------------------------
2619// This function writes the regularization parameters to the DB
2620//-----------------------------------------------------------------------------
2621
2622 Option_t *opt;
2623 Char_t *dirName="Pions";
2624 Char_t *keyName="RegPions";
2625
2626 SetParticle(pdg);
2627
2628 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2629 } else { opt="update"; }
2630
2631 switch (pdg) {
2632 case 211:
2633 dirName="Pions";
2634 keyName="RegPions";
2635 break;
2636 case 321:
2637 dirName="Kaons";
2638 keyName="RegKaons";
2639 break;
2640 case 11:
2641 dirName="Electrons";
2642 keyName="RegElectrons";
2643 break;
2644 }
70521312 2645
6eb69b9a 2646 TFile *outFile = new TFile(outName,opt);
2647 if(!gDirectory->cd("/RegParams")) {
2648 outFile->mkdir("RegParams");
2649 gDirectory->cd("/RegParams");
6eb67451 2650 }
6eb69b9a 2651 TDirectory *dir2 = gDirectory->mkdir(dirName);
2652 dir2->cd();
2653 fRegPar->Write(keyName);
6eb67451 2654
6eb69b9a 2655 outFile->Close();
2656 delete outFile;
70521312 2657
2658
6eb69b9a 2659 return 1;
6eb67451 2660}
6eb67451 2661
2662
2663
2664
2665
2666
2667
2668
2669