1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 Revision 1.5.4.1 2002/06/10 15:26:12 hristov
21 Revision 1.9 2002/05/13 09:53:08 hristov
22 Some frequent problems corrected: arrays with variable size have to be defined via the operator new, default values for the arguments have to be used only in the header files, etc.
24 Revision 1.8 2002/05/08 18:21:40 kowal2
25 Now the code is blind to the binning used for pulls, efficiencies etc.
27 Revision 1.7 2002/04/10 16:30:07 kowal2
33 /**************************************************************************
35 * This class builds AliTPCtrack objects from generated tracks to feed *
36 * ITS tracking (V2). The AliTPCtrack is built from its first hit in *
37 * the TPC. The track is assigned a Kalman-like covariance matrix *
38 * depending on its pT and pseudorapidity and track parameters are *
39 * smeared according to this covariance matrix. *
40 * Output file contains sorted tracks, ready for matching with ITS. *
43 * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm *
45 * Test macro is: AliBarrelRec_TPCparam.C *
47 * 2002/04/28: Major upgrade of the class *
48 * - Covariance matrices and pulls are now separeted for pions, kaons *
50 * - A parameterization of dE/dx in the TPC has been included and it is *
51 * used to assign a mass to each track according to a rough dE/dx PID. *
52 * - All the "numbers" have been moved to the file with the DataBase, *
53 * they are read as objects of the class AliTPCkineGrid, and assigned *
54 * to data memebers of the class AliTPCtrackerParam. *
55 * - All the code necessary to create a BataBase has been included in *
56 * class (see the macro AliTPCtrackingParamDB.C for the details). *
58 * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it *
60 **************************************************************************/
64 #include <TGraphErrors.h>
71 #include "AliGausCorr.h"
72 #include "AliKalmanTrack.h"
74 #include "AliMagFCM.h"
75 #include "AliTPCkineGrid.h"
76 #include "AliTPCtrack.h"
77 #include "AliTPCtrackerParam.h"
79 Double_t RegFunc(Double_t *x,Double_t *par) {
80 // This is the function used to regularize the covariance matrix
81 Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
82 TMath::Power(x[1],par[3]);
85 Double_t FitRegFunc(Double_t *x,Double_t *par) {
86 // This is the function used for the fit of the covariance
87 // matrix as a function of the momentum.
88 // For cosl the average value 0.87 is used.
89 Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
90 TMath::Power(0.87,par[3]);
94 // structure for DB building
104 Double_t dP0,dP1,dP2,dP3,dP4;
105 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
108 // cov matrix structure
110 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
113 ClassImp(AliTPCtrackerParam)
115 //-----------------------------------------------------------------------------
116 AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
118 //-----------------------------------------------------------------------------
119 // This is the class conctructor
120 //-----------------------------------------------------------------------------
122 fBz = Bz; // value of the z component of L3 field (Tesla)
123 fColl = coll; // collision code (0: PbPb6000)
124 fSelAndSmear = kTRUE; // by default selection and smearing are done
127 cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n";
128 cerr<<" Available: 0.4"<<endl;
131 cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n";
132 cerr<<" Available: 0 -> PbPb6000"<<endl;
135 fDBfileName = gSystem->Getenv("ALICE_ROOT");
136 fDBfileName.Append("/TPC/CovMatrixDB_");
137 if(fColl==0) fDBfileName.Append("PbPb6000");
138 if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
140 //-----------------------------------------------------------------------------
141 AliTPCtrackerParam::~AliTPCtrackerParam()
143 //-----------------------------------------------------------------------------
145 Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
147 //-----------------------------------------------------------------------------
148 // This function creates the TPC parameterized tracks
149 //-----------------------------------------------------------------------------
152 TTree *covTreePi[50];
153 TTree *covTreeKa[50];
154 TTree *covTreeEl[50];
157 cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
158 fDBfileName.Data()<<"\n+++\n";
159 // Read paramters from DB file
160 if(!ReadAllData(fDBfileName.Data())) {
161 cerr<<"AliTPCtrackerParam::BuildTPCtracks: \
162 Could not read data from DB\n\n"; return 1;
164 // Read the trees with regularized cov. matrices from DB
166 fileDB = TFile::Open(fDBfileName.Data());
167 Int_t nBinsPi = fDBgridPi.GetTotBins();
168 for(Int_t l=0; l<nBinsPi; l++) {
169 str = "/CovMatrices/Pions/CovTreePi_bin";
171 covTreePi[l] = (TTree*)fileDB->Get(str.Data());
173 Int_t nBinsKa = fDBgridKa.GetTotBins();
174 for(Int_t l=0; l<nBinsKa; l++) {
175 str = "/CovMatrices/Kaons/CovTreeKa_bin";
177 covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
179 Int_t nBinsEl = fDBgridEl.GetTotBins();
180 for(Int_t l=0; l<nBinsEl; l++) {
181 str = "/CovMatrices/Electrons/CovTreeEl_bin";
183 covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
186 } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
188 TFile *infile=(TFile*)inp;
191 // Get gAlice object from file
192 if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
193 cerr<<"gAlice has not been found on galice.root !\n";
197 // Check if value in the galice file is equal to selected one (fBz)
198 AliMagFCM *fiel = (AliMagFCM*)gAlice->Field();
199 Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
200 printf("Magnetic field is %6.2f Tesla\n",fieval);
202 cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!"<<endl;
203 cerr<<"Field selected is: "<<fBz<<" T\n";
204 cerr<<"Field found on file is: "<<fieval<<" T\n";
209 AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
210 Int_t ver = TPC->IsVersion();
211 cerr<<"+++ TPC version "<<ver<<" has been found !\n";
212 AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
215 digp = new AliTPCParamSR();
217 else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
219 if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
222 // Set the conversion constant between curvature and Pt
223 AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
227 AliTPCtrack *tpctrack=0;
228 Double_t hPx,hPy,hPz,hPt,hEta,xg,yg,zg,xl,yl,zl;
230 Float_t cosAlpha,sinAlpha;
231 Int_t bin,label,pdg,charge;
233 Int_t nParticles,nTracks,arrentr;
235 //Int_t nSel=0,nAcc=0;
237 // loop over first n events in file
238 for(Int_t evt=0; evt<n; evt++){
239 cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
241 // tree for TPC tracks
242 sprintf(tname,"TreeT_TPC_%d",evt);
243 TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
244 tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
246 // array for TPC tracks
247 TObjArray tarray(20000);
249 // get the particles stack
250 nParticles = gAlice->GetEvent(evt);
251 Bool_t * done = new Bool_t[nParticles];
252 Int_t * pdgCodes = new Int_t[nParticles];
254 // loop on particles and store pdg codes
255 for(Int_t l=0; l<nParticles; l++) {
256 Part = gAlice->Particle(l);
257 pdgCodes[l] = Part->GetPdgCode();
261 // Get TreeH with hits
262 TTree *TH=gAlice->TreeH();
263 nTracks=(Int_t)TH->GetEntries();
264 cerr<<"+++\n+++ Number of particles in event "<<evt<<": "<<nParticles<<
265 "\n+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
268 // loop over entries in TreeH
269 for(Int_t l=0; l<nTracks; l++) {
270 if(l%1000==0) cerr<<" --- Processing primary track "
271 <<l<<" of "<<nTracks<<" ---\r";
275 tpcHit=(AliTPChit*)TPC->FirstHit(-1);
276 for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
277 if(tpcHit->fQ !=0.) continue;
278 // Get particle momentum at hit
279 hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
280 hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
281 // reject hits with Pt<mag*0.45 GeV/c
282 if(hPt<(fBz*0.45)) continue;
285 label=tpcHit->Track();
286 // check if this track has already been processed
287 if(done[label]) continue;
288 // PDG code & electric charge
289 pdg = pdgCodes[label];
290 if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
291 else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
293 pdg = TMath::Abs(pdg);
294 if(pdg>3000) pdg=211;
295 if(fSelAndSmear) SetParticle(pdg);
297 if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
298 if(tpcHit->fQ != 0.) continue;
299 // Get global coordinates of hit
300 xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
301 if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
303 // Get TPC sector, Alpha angle and local coordinates
304 // cerr<<"Sector "<<tpcHit->fSector<<endl;
305 digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
306 alpha = TMath::ATan2(sinAlpha,cosAlpha);
307 xl = xg*cosAlpha + yg*sinAlpha;
308 yl =-xg*sinAlpha + yg*cosAlpha;
310 //printf("Alpha %f xl %f yl %f zl %f\n",Alpha,xl,yl,zl);
312 // reject tracks which are not in the TPC acceptance
313 if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
315 hEta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(hPz/hPt)));
317 // Apply selection according to TPC efficiency
318 //if(TMath::Abs(pdg)==211) nAcc++;
319 if(fSelAndSmear && !SelectedTrack(hPt,hEta)) continue;
320 //if(TMath::Abs(pdg)==211) nSel++;
322 // create AliTPCtrack object
323 BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge);
326 bin = fDBgrid->GetBin(hPt,hEta);
329 fCovTree = covTreePi[bin];
332 fCovTree = covTreeKa[bin];
335 fCovTree = covTreePi[bin];
338 fCovTree = covTreeEl[bin];
341 fCovTree = covTreePi[bin];
344 // deal with covariance matrix and smearing of parameters
347 // assign the track a dE/dx and make a rough PID
351 // put track in array
352 AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
353 iotrack->SetLabel(label);
354 tarray.AddLast(iotrack);
355 // Mark track as "done" and register the pdg code
360 } // loop over entries in TreeH
362 // sort array with TPC tracks (decreasing pT)
365 arrentr = tarray.GetEntriesFast();
366 for(Int_t l=0; l<arrentr; l++) {
367 tpctrack=(AliTPCtrack*)tarray.UncheckedAt(l);
371 // write the tree with tracks in the output file
379 printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
380 //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
384 if(fileDB) fileDB->Close();
388 //-----------------------------------------------------------------------------
389 void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
390 //-----------------------------------------------------------------------------
391 // This function computes the dE/dx for pions, kaons, protons and electrons,
392 // in the [pT,eta] bins.
393 // Input file is CovMatrix_AllEvts.root.
394 //-----------------------------------------------------------------------------
396 gStyle->SetOptStat(0);
397 gStyle->SetOptFit(10001);
399 Char_t *part="PIONS";
402 // create a chain with compared tracks
403 TChain cmptrkchain("cmptrktree");
404 cmptrkchain.Add("CovMatrix_AllEvts_1.root");
405 cmptrkchain.Add("CovMatrix_AllEvts_2.root");
406 cmptrkchain.Add("CovMatrix_AllEvts_3.root");
409 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
410 Int_t entries = (Int_t)cmptrkchain.GetEntries();
411 cerr<<" Number of entries: "<<entries<<endl;
413 InitializeKineGrid("DB","points");
414 InitializeKineGrid("dEdx","bins");
437 const Int_t nTotBins = fDBgrid->GetTotBins();
439 cerr<<" Fit bins: "<<nTotBins<<endl;
442 Int_t * n = new Int_t[nTotBins];
443 Double_t * p = new Double_t[nTotBins];
444 Double_t * ep = new Double_t[nTotBins];
445 Double_t * mean = new Double_t[nTotBins];
446 Double_t * sigma = new Double_t[nTotBins];
448 for(Int_t l=0; l<nTotBins; l++) {
449 n[l] = 1; // set to 1 to avoid divisions by 0
450 p[l] = mean[l] = sigma[l] = ep[l] = 0.;
453 // loop on chain entries for the mean
454 for(Int_t l=0; l<entries; l++) {
455 cmptrkchain.GetEvent(l);
456 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
457 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
459 mean[bin] += cmptrk.dEdx;
461 } // loop on chain entries
463 for(Int_t l=0; l<nTotBins; l++) {
466 n[l] = 1; // set to 1 to avoid divisions by 0
469 // loop on chain entries for the sigma
470 for(Int_t l=0; l<entries; l++) {
471 cmptrkchain.GetEvent(l);
472 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
473 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
474 if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
476 sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
477 } // loop on chain entries
479 for(Int_t l=0; l<nTotBins; l++) {
480 sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
484 TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700);
486 // create the graph for dEdx vs p
487 TGraphErrors *gr = new TGraphErrors(nTotBins,p,mean,ep,sigma);
488 TString title(" : dE/dx vs momentum"); title.Prepend(part);
489 TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
490 frame->SetXTitle("p [GeV/c]");
491 frame->SetYTitle("dE/dx [a.u.]");
498 for(Int_t i=0; i<nTotBins; i++) {
499 fdEdxMeanPi.SetParam(i,mean[i]);
500 fdEdxRMSPi.SetParam(i,sigma[i]);
504 for(Int_t i=0; i<nTotBins; i++) {
505 fdEdxMeanKa.SetParam(i,mean[i]);
506 fdEdxRMSKa.SetParam(i,sigma[i]);
510 for(Int_t i=0; i<nTotBins; i++) {
511 fdEdxMeanPr.SetParam(i,mean[i]);
512 fdEdxRMSPr.SetParam(i,sigma[i]);
516 for(Int_t i=0; i<nTotBins; i++) {
517 fdEdxMeanEl.SetParam(i,mean[i]);
518 fdEdxRMSEl.SetParam(i,sigma[i]);
523 // write results to file
524 WritedEdx(outName,pdg);
534 //-----------------------------------------------------------------------------
535 void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
536 //-----------------------------------------------------------------------------
537 // This function computes the pulls for pions, kaons and electrons,
538 // in the [pT,eta] bins.
539 // Input file is CovMatrix_AllEvts.root.
540 // Output file is pulls.root.
541 //-----------------------------------------------------------------------------
543 // create a chain with compared tracks
544 TChain cmptrkchain("cmptrktree");
545 cmptrkchain.Add("CovMatrix_AllEvts_1.root");
546 cmptrkchain.Add("CovMatrix_AllEvts_2.root");
547 cmptrkchain.Add("CovMatrix_AllEvts_3.root");
550 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
551 Int_t entries = (Int_t)cmptrkchain.GetEntries();
552 cerr<<" Number of entries: "<<entries<<endl;
555 Char_t hname[100], htitle[100];
558 AliTPCkineGrid pulls[5];
559 TH1F *hDum = new TH1F("name","title",100,-7.,7.);
560 TF1 *g = new TF1("g","gaus");
562 InitializeKineGrid("pulls","bins");
563 InitializeKineGrid("DB","points");
567 // loop on the particles Pi,Ka,El
568 for(Int_t part=0; part<3; part++) {
573 cerr<<" Processing pions ...\n";
577 cerr<<" Processing kaons ...\n";
582 cerr<<" Processing electrons ...\n";
586 SetParticle(thisPdg);
588 for(Int_t i=0;i<5;i++) {
589 pulls[i].~AliTPCkineGrid();
590 new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
592 nTotBins = fDBgrid->GetTotBins();
594 // create histograms for the all the bins
601 hPulls0_ = new TH1F[nTotBins];
602 hPulls1_ = new TH1F[nTotBins];
603 hPulls2_ = new TH1F[nTotBins];
604 hPulls3_ = new TH1F[nTotBins];
605 hPulls4_ = new TH1F[nTotBins];
608 for(Int_t i=0; i<nTotBins; i++) {
609 sprintf(hname,"hPulls0_%d",i);
610 sprintf(htitle,"P0 pulls for bin %d",i);
611 hDum->SetName(hname); hDum->SetTitle(htitle);
613 sprintf(hname,"hPulls1_%d",i);
614 sprintf(htitle,"P1 pulls for bin %d",i);
615 hDum->SetName(hname); hDum->SetTitle(htitle);
617 sprintf(hname,"hPulls2_%d",i);
618 sprintf(htitle,"P2 pulls for bin %d",i);
619 hDum->SetName(hname); hDum->SetTitle(htitle);
621 sprintf(hname,"hPulls3_%d",i);
622 sprintf(htitle,"P3 pulls for bin %d",i);
623 hDum->SetName(hname); hDum->SetTitle(htitle);
625 sprintf(hname,"hPulls4_%d",i);
626 sprintf(htitle,"P4 pulls for bin %d",i);
627 hDum->SetName(hname); hDum->SetTitle(htitle);
631 // loop on chain entries
632 for(Int_t i=0; i<entries; i++) {
633 cmptrkchain.GetEvent(i);
634 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
635 // fill histograms with the pulls
636 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
637 hPulls0_[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
638 hPulls1_[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
639 hPulls2_[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
640 hPulls3_[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
641 hPulls4_[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
642 } // loop on chain entries
644 // compute the sigma of the distributions
645 for(Int_t i=0; i<nTotBins; i++) {
646 if(hPulls0_[i].GetEntries()>1000) {
647 g->SetRange(-3.*hPulls0_[i].GetRMS(),3.*hPulls0_[i].GetRMS());
648 hPulls0_[i].Fit("g","R,Q,N");
649 pulls[0].SetParam(i,g->GetParameter(2));
650 } else pulls[0].SetParam(i,-1.);
651 if(hPulls1_[i].GetEntries()>1000) {
652 g->SetRange(-3.*hPulls1_[i].GetRMS(),3.*hPulls1_[i].GetRMS());
653 hPulls1_[i].Fit("g","R,Q,N");
654 pulls[1].SetParam(i,g->GetParameter(2));
655 } else pulls[1].SetParam(i,-1.);
656 if(hPulls2_[i].GetEntries()>1000) {
657 g->SetRange(-3.*hPulls2_[i].GetRMS(),3.*hPulls2_[i].GetRMS());
658 hPulls2_[i].Fit("g","R,Q,N");
659 pulls[2].SetParam(i,g->GetParameter(2));
660 } else pulls[2].SetParam(i,-1.);
661 if(hPulls3_[i].GetEntries()>1000) {
662 g->SetRange(-3.*hPulls3_[i].GetRMS(),3.*hPulls3_[i].GetRMS());
663 hPulls3_[i].Fit("g","R,Q,N");
664 pulls[3].SetParam(i,g->GetParameter(2));
665 } else pulls[3].SetParam(i,-1.);
666 if(hPulls4_[i].GetEntries()>1000) {
667 g->SetRange(-3.*hPulls4_[i].GetRMS(),3.*hPulls4_[i].GetRMS());
668 hPulls4_[i].Fit("g","R,Q,N");
669 pulls[4].SetParam(i,g->GetParameter(2));
670 } else pulls[4].SetParam(i,-1.);
676 for(Int_t i=0;i<5;i++) {
677 fPullsPi[i].~AliTPCkineGrid();
678 new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
682 for(Int_t i=0;i<5;i++) {
683 fPullsKa[i].~AliTPCkineGrid();
684 new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
688 for(Int_t i=0;i<5;i++) {
689 fPullsEl[i].~AliTPCkineGrid();
690 new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
701 } // loop on particle species
703 // write pulls to file
709 //-----------------------------------------------------------------------------
710 void AliTPCtrackerParam::BuildTrack(Double_t alpha,
711 Double_t x,Double_t y,Double_t z,
712 Double_t px,Double_t py,
713 Double_t pz,Double_t pt,
715 //-----------------------------------------------------------------------------
716 // This function uses GEANT info to set true track parameters
717 //-----------------------------------------------------------------------------
719 Double_t xx[5],cc[15];
720 cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
721 cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
724 TVector3 bfield(0.,0.,fBz);
727 // radius [cm] of track projection in (x,y)
728 Double_t rho = pt*100./0.299792458/bfield.Z();
729 // center of track projection in local reference frame
733 // position (local) and momentum (local) at the hit
734 // in the bending plane (z=0)
736 hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.);
737 TVector3 vrho = hmom.Cross(bfield);
741 TVector3 vcenter = hpos+vrho;
743 Double_t x0 = vcenter.X();
745 // fX = xref X-coordinate of this track (reference plane)
746 // fAlpha = Alpha Rotation angle the local (TPC sector)
747 // fP0 = YL Y-coordinate of a track
748 // fP1 = ZG Z-coordinate of a track
749 // fP2 = C*x0 x0 is center x in rotated frame
750 // fP3 = Tgl tangent of the track momentum dip angle
751 // fP4 = C track curvature
758 // create the object AliTPCtrack
759 AliTPCtrack track(0,xx,cc,xref,alpha);
760 new(&fTrack) AliTPCtrack(track);
764 //-----------------------------------------------------------------------------
765 void AliTPCtrackerParam::CompareTPCtracks(
766 const Char_t* galiceName,
767 const Char_t* trkGeaName,
768 const Char_t* trkKalName,
769 const Char_t* covmatName,
770 const Char_t* tpceffName) const {
771 //-----------------------------------------------------------------------------
772 // This function compares tracks from TPC Kalman Filter V2 with
773 // true tracks at TPC 1st hit. It gives:
774 // - a tree with Kalman cov. matrix elements, resolutions, dEdx
775 // - the efficiencies as a function of particle type, pT, eta
776 //-----------------------------------------------------------------------------
778 TFile *kalFile = TFile::Open(trkKalName);
779 TFile *geaFile = TFile::Open(trkGeaName);
780 TFile *galiceFile = TFile::Open(galiceName);
782 // particles from TreeK
783 AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
784 const Int_t nparticles = gAlice->GetEvent(0);
786 Int_t * kalLab = new Int_t[nparticles];
787 for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1;
790 // tracks from Kalman
791 TTree *kaltree=(TTree*)kalFile->Get("TreeT_TPC_0");
792 AliTPCtrack *kaltrack=new AliTPCtrack;
793 kaltree->SetBranchAddress("tracks",&kaltrack);
794 Int_t kalEntries = (Int_t)kaltree->GetEntries();
796 // tracks from 1st hit
797 TTree *geatree=(TTree*)geaFile->Get("TreeT_TPC_0");
798 AliTPCtrack *geatrack=new AliTPCtrack;
799 geatree->SetBranchAddress("tracks",&geatrack);
800 Int_t geaEntries = (Int_t)geatree->GetEntries();
802 cerr<<"+++\n+++ Number of tracks: TPC Kalman = "<<kalEntries<<endl<<"+++ TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
804 // set pointers for TPC tracks:
805 // the entry number of the track labelled l is stored in kalLab[l]
807 for (Int_t j=0; j<kalEntries; j++) {
808 kaltree->GetEvent(j);
809 if(kaltrack->GetLabel()>=0) {
810 if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
811 kalLab[kaltrack->GetLabel()] = j;
816 cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
817 cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
821 Double_t cc[15],dAlpha;
822 // --- [Pt,eta] binning ---
824 // |eta|<0.3 0.3<=|eta|<0.6 0.6<=|eta|
828 // 1.0<=Pt<1.5 9 10 11
829 // 1.5<=Pt<3.0 12 13 14
830 // 3.0<=Pt<5.0 15 16 17
831 // 5.0<=Pt<7.0 18 19 20
832 // 7.0<=Pt<15. 21 22 23
835 Int_t pi=0,ka=0,mu=0,el=0,pr=0;
836 Int_t geaPi[27],geaKa[27],geaPr[27],geaEl[27],geaMu[27];
837 Int_t kalPi[27],kalKa[27],kalPr[27],kalEl[27],kalMu[27];
838 Float_t effPi[27],effKa[27],effPr[27],effEl[27],effMu[27];
840 for(Int_t j=0; j<27; j++) {
841 geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
842 kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
843 effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
847 // create the tree for comparison results
848 TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
849 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");
851 cerr<<"Doing track comparison...\n";
852 // loop on tracks at 1st hit
853 for(Int_t j=0; j<geaEntries; j++) {
854 geatree->GetEvent(j);
856 label = geatrack->GetLabel();
857 Part = (TParticle*)gAlice->Particle(label);
858 cmptrk.pdg = Part->GetPdgCode();
859 cmptrk.eta = Part->Eta();
860 cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
862 cmptrk.pt = 1/TMath::Abs(geatrack->Get1Pt());
863 cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
864 cmptrk.p = cmptrk.pt/cmptrk.cosl;
867 bin = GetBin(cmptrk.pt,cmptrk.eta);
869 if(abs(cmptrk.pdg)==211) geaPi[bin]++;
870 if(abs(cmptrk.pdg)==321) geaKa[bin]++;
871 if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
872 if(abs(cmptrk.pdg)==11) geaEl[bin]++;
873 if(abs(cmptrk.pdg)==13) geaMu[bin]++;
876 // check if there is the corresponding track in the TPC kalman and get it
877 if(kalLab[label]==-1) continue;
878 kaltree->GetEvent(kalLab[label]);
880 // and go on only if it has xref = 84.57 cm (inner pad row)
881 if(kaltrack->GetX()>90.) continue;
883 if(abs(cmptrk.pdg)==211) { kalPi[bin]++; pi++; }
884 if(abs(cmptrk.pdg)==321) { kalKa[bin]++; ka++; }
885 if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
886 if(abs(cmptrk.pdg)==11) { kalEl[bin]++; el++; }
887 if(abs(cmptrk.pdg)==13) { kalMu[bin]++; mu++; }
889 kaltrack->PropagateTo(geatrack->GetX());
891 cmptrk.dEdx = kaltrack->GetdEdx();
893 // compute errors on parameters
894 dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
895 if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
897 cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
898 cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
899 cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
900 cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
901 cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
902 cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
904 // get covariance matrix
905 // beware: lines 3 and 4 are inverted!
906 kaltrack->GetCovariance(cc);
927 } // loop on tracks at TPC 1st hit
929 cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
930 cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
933 // Write tree to file
934 TFile *outfile = new TFile(covmatName,"recreate");
939 // Write efficiencies to file
940 FILE *effFile = fopen(tpceffName,"w");
941 //fprintf(effFile,"%d\n",kalEntries);
942 for(Int_t j=0; j<27; j++) {
943 if(geaPi[j]>=5) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
944 if(geaKa[j]>=5) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
945 if(geaPr[j]>=5) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
946 if(geaEl[j]>=5) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
947 if(geaMu[j]>=5) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
948 fprintf(effFile,"%f %f %f %f %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
951 for(Int_t j=0; j<27; j++) {
952 fprintf(effFile,"%d %d %d %d %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
954 for(Int_t j=0; j<27; j++) {
955 fprintf(effFile,"%d %d %d %d %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
959 // delete AliRun object
960 delete gAlice; gAlice=0;
962 // close all input files
971 //-----------------------------------------------------------------------------
972 void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
973 //-----------------------------------------------------------------------------
974 // This function assigns the track a dE/dx and makes a rough PID
975 //-----------------------------------------------------------------------------
977 Double_t mean = fdEdxMean->GetValueAt(pt,eta);
978 Double_t rms = fdEdxRMS->GetValueAt(pt,eta);
980 Double_t dEdx = gRandom->Gaus(mean,rms);
982 fTrack.SetdEdx(dEdx);
984 AliTPCtrackParam t(fTrack);
987 Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
990 if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
991 t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
993 if (dEdx < 39.+ 12./p/p) {
994 t.AssignMass(0.49368); new(&fTrack) AliTPCtrack(t); return;
996 t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
1000 if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
1001 t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1003 t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
1006 t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1008 //-----------------------------------------------------------------------------
1009 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
1010 //-----------------------------------------------------------------------------
1011 // This function deals with covariance matrix and smearing
1012 //-----------------------------------------------------------------------------
1016 Double_t trkKine[2],trkRegPar[4];
1017 Double_t xref,alpha,xx[5],xxsm[5],cc[15];
1020 fCovTree->SetBranchAddress("matrix",&covmat);
1022 // get random entry from the tree
1023 treeEntries = (Int_t)fCovTree->GetEntries();
1024 fCovTree->GetEvent(gRandom->Integer(treeEntries));
1026 // get P and Cosl from track
1027 cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1028 p = 1./TMath::Abs(fTrack.Get1Pt())/cosl;
1033 // get covariance matrix from regularized matrix
1034 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(0,l);
1035 cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1037 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(1,l);
1038 cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1039 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(2,l);
1040 cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1042 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(3,l);
1043 cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1045 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(4,l);
1046 cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1048 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(5,l);
1049 cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1050 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(6,l);
1051 cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1053 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(7,l);
1054 cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1056 for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(8,l);
1057 cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1059 TMatrixD covMatSmear(5,5);
1061 covMatSmear = GetSmearingMatrix(cc,pt,eta);
1063 // get track original parameters
1064 xref =fTrack.GetX();
1065 alpha=fTrack.GetAlpha();
1066 xx[0]=fTrack.GetY();
1067 xx[1]=fTrack.GetZ();
1068 xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp();
1069 xx[3]=fTrack.GetTgl();
1070 xx[4]=fTrack.GetC();
1072 // use smearing matrix to smear the original parameters
1073 SmearTrack(xx,xxsm,covMatSmear);
1075 AliTPCtrack track(0,xxsm,cc,xref,alpha);
1076 new(&fTrack) AliTPCtrack(track);
1080 //-----------------------------------------------------------------------------
1081 void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg) {
1082 //-----------------------------------------------------------------------------
1083 // This function draws the TPC efficiencies in the [pT,eta] bins
1084 //-----------------------------------------------------------------------------
1089 const Int_t n = fEff->GetPointsPt();
1090 Double_t * effsA = new Double_t[n];
1091 Double_t * effsB = new Double_t[n];
1092 Double_t * effsC = new Double_t[n];
1093 Double_t * pt = new Double_t[n];
1095 fEff->GetArrayPt(pt);
1096 for(Int_t i=0;i<n;i++) {
1097 effsA[i] = fEff->GetParam(i,0);
1098 effsB[i] = fEff->GetParam(i,1);
1099 effsC[i] = fEff->GetParam(i,2);
1102 TGraph *grA = new TGraph(n,pt,effsA);
1103 TGraph *grB = new TGraph(n,pt,effsB);
1104 TGraph *grC = new TGraph(n,pt,effsC);
1106 grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1107 grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1108 grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1110 TString title("Distribution of the TPC efficiencies");
1113 title.Prepend("PIONS - ");
1116 title.Prepend("KAONS - ");
1119 title.Prepend("PROTONS - ");
1122 title.Prepend("ELECTRONS - ");
1125 title.Prepend("MUONS - ");
1130 gStyle->SetOptStat(0);
1131 TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1136 TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1137 frame->SetXTitle("p_{T} [GeV/c]");
1143 TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1144 leg->AddEntry(grA,"|#eta|<0.3","p");
1145 leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1146 leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1147 leg->SetFillColor(0);
1157 //-----------------------------------------------------------------------------
1158 void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg,
1160 //-----------------------------------------------------------------------------
1161 // This function draws the pulls in the [pT,eta] bins
1162 //-----------------------------------------------------------------------------
1167 const Int_t n = (fPulls+par)->GetPointsPt();
1168 Double_t * pullsA = new Double_t[n];
1169 Double_t * pullsB = new Double_t[n];
1170 Double_t * pullsC = new Double_t[n];
1171 Double_t * pt = new Double_t[n];
1172 (fPulls+par)->GetArrayPt(pt);
1173 for(Int_t i=0;i<n;i++) {
1174 pullsA[i] = (fPulls+par)->GetParam(i,0);
1175 pullsB[i] = (fPulls+par)->GetParam(i,1);
1176 pullsC[i] = (fPulls+par)->GetParam(i,2);
1179 TGraph *grA = new TGraph(n,pt,pullsA);
1180 TGraph *grB = new TGraph(n,pt,pullsB);
1181 TGraph *grC = new TGraph(n,pt,pullsC);
1183 grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1184 grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1185 grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1187 TString title("Distribution of the pulls: ");
1190 title.Prepend("PIONS - ");
1193 title.Prepend("KAONS - ");
1196 title.Prepend("ELECTRONS - ");
1207 title.Append(" #eta");
1210 title.Append("tg #lambda");
1217 gStyle->SetOptStat(0);
1218 TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1223 TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1224 frame->SetXTitle("p_{T} [GeV/c]");
1230 TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1231 leg->AddEntry(grA,"|#eta|<0.3","p");
1232 leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1233 leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1234 leg->SetFillColor(0);
1244 //-----------------------------------------------------------------------------
1245 Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
1246 //-----------------------------------------------------------------------------
1247 // This function tells bin number in a grid (pT,eta)
1248 //-----------------------------------------------------------------------------
1249 if(TMath::Abs(eta)<0.3) {
1250 if(pt<0.3) return 0;
1251 if(pt>=0.3 && pt<0.5) return 3;
1252 if(pt>=0.5 && pt<1.) return 6;
1253 if(pt>=1. && pt<1.5) return 9;
1254 if(pt>=1.5 && pt<3.) return 12;
1255 if(pt>=3. && pt<5.) return 15;
1256 if(pt>=5. && pt<7.) return 18;
1257 if(pt>=7. && pt<15.) return 21;
1258 if(pt>=15.) return 24;
1260 if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
1261 if(pt<0.3) return 1;
1262 if(pt>=0.3 && pt<0.5) return 4;
1263 if(pt>=0.5 && pt<1.) return 7;
1264 if(pt>=1. && pt<1.5) return 10;
1265 if(pt>=1.5 && pt<3.) return 13;
1266 if(pt>=3. && pt<5.) return 16;
1267 if(pt>=5. && pt<7.) return 19;
1268 if(pt>=7. && pt<15.) return 22;
1269 if(pt>=15.) return 25;
1271 if(TMath::Abs(eta)>=0.6) {
1272 if(pt<0.3) return 2;
1273 if(pt>=0.3 && pt<0.5) return 5;
1274 if(pt>=0.5 && pt<1.) return 8;
1275 if(pt>=1. && pt<1.5) return 11;
1276 if(pt>=1.5 && pt<3.) return 14;
1277 if(pt>=3. && pt<5.) return 17;
1278 if(pt>=5. && pt<7.) return 20;
1279 if(pt>=7. && pt<15.) return 23;
1280 if(pt>=15.) return 26;
1286 //-----------------------------------------------------------------------------
1287 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1288 Double_t eta) const {
1289 //-----------------------------------------------------------------------------
1290 // This function stretches the covariance matrix according to the pulls
1291 //-----------------------------------------------------------------------------
1293 TMatrixD covMat(5,5);
1296 covMat(1,0)=cc[1]; covMat(0,1)=covMat(1,0);
1298 covMat(2,0)=cc[3]; covMat(0,2)=covMat(2,0);
1299 covMat(2,1)=cc[4]; covMat(1,2)=covMat(2,1);
1301 covMat(3,0)=cc[6]; covMat(0,3)=covMat(3,0);
1302 covMat(3,1)=cc[7]; covMat(1,3)=covMat(3,1);
1303 covMat(3,2)=cc[8]; covMat(2,3)=covMat(3,2);
1305 covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1306 covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1307 covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1308 covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1311 TMatrixD stretchMat(5,5);
1312 for(Int_t k=0;k<5;k++) {
1313 for(Int_t l=0;l<5;l++) {
1318 for(Int_t i=0;i<5;i++) stretchMat(i,i) =
1319 TMath::Sqrt((fPulls+i)->GetValueAt(pt,eta));
1321 TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1322 TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1326 //-----------------------------------------------------------------------------
1327 void AliTPCtrackerParam::InitializeKineGrid(Option_t* which,Option_t* how) {
1328 //-----------------------------------------------------------------------------
1329 // This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1330 // which = "DB" -> initialize fDBgrid... members
1331 // "eff" -> initialize fEff... members
1332 // "pulls" -> initialize fPulls... members
1333 // "dEdx" -> initialize fdEdx... members
1334 // how = "points" -> initialize with the points of the grid
1335 // "bins" -> initialize with the bins of the grid
1336 //-----------------------------------------------------------------------------
1338 const char *points = strstr(how,"points");
1339 const char *bins = strstr(how,"bins");
1340 const char *DB = strstr(which,"DB");
1341 const char *eff = strstr(which,"eff");
1342 const char *pulls = strstr(which,"pulls");
1343 const char *dEdx = strstr(which,"dEdx");
1346 Int_t nEta=0,nPtPi=0,nPtKa=0,nPtPr=0,nPtEl=0,nPtMu=0;
1348 Double_t etaPoints[2] = {0.3,0.6};
1349 Double_t etaBins[3] = {0.15,0.45,0.75};
1350 Double_t ptPoints8[8] = {0.3,0.5,1.,1.5,3.,5.,7.,15.};
1351 Double_t ptPoints5[5] = {0.3,0.5,1.,1.5,5.};
1352 Double_t ptBins9[9] = {0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
1353 Double_t ptBins6[6] = {0.244,0.390,0.676,1.190,2.36,10.};
1355 Double_t *eta=0,*ptPi=0,*ptKa=0,*ptPr=0,*ptEl=0,*ptMu=0;
1386 AliTPCkineGrid *dummy=0;
1389 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1390 new(&fDBgridPi) AliTPCkineGrid(*dummy);
1392 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1393 new(&fDBgridKa) AliTPCkineGrid(*dummy);
1395 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1396 new(&fDBgridEl) AliTPCkineGrid(*dummy);
1400 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1401 new(&fEffPi) AliTPCkineGrid(*dummy);
1403 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1404 new(&fEffKa) AliTPCkineGrid(*dummy);
1406 dummy = new AliTPCkineGrid(nPtPr,nEta,ptPr,eta);
1407 new(&fEffPr) AliTPCkineGrid(*dummy);
1409 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1410 new(&fEffEl) AliTPCkineGrid(*dummy);
1412 dummy = new AliTPCkineGrid(nPtMu,nEta,ptMu,eta);
1413 new(&fEffMu) AliTPCkineGrid(*dummy);
1417 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1418 for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1420 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1421 for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1423 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1424 for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1428 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1429 new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1430 new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1432 dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1433 new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1434 new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1436 dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1437 new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1438 new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1440 dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1441 new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1442 new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1448 //-----------------------------------------------------------------------------
1449 void AliTPCtrackerParam::MakeDataBase() {
1450 //-----------------------------------------------------------------------------
1451 // This function creates the DB file and store in it:
1452 // - TPC Efficiencies for pi,ka,pr,el,mu (AliTPCkineGrid class)
1453 // - Pulls for pi,ka,el (AliTPCkineGrid class)
1454 // - Regularization parameters for pi,ka,el (TMatrixD class)
1455 // - dE/dx parameterization for pi,ka,pr,el (AliTPCkineGrid class)
1456 // - Regularized cov. matrices for pi,ka,el (COVMATRIX structure)
1457 //-----------------------------------------------------------------------------
1459 // define some file names
1460 Char_t *effFile ="CovMatrixDB_partial.root";
1461 Char_t *pullsFile="CovMatrixDB_partial.root";
1462 Char_t *regPiFile="CovMatrixDB_partial.root";
1463 Char_t *regKaFile="CovMatrixDB_partial.root";
1464 Char_t *regElFile="CovMatrixDB_partial.root";
1465 Char_t *dEdxPiFile="dEdxPi.root";
1466 Char_t *dEdxKaFile="dEdxKa.root";
1467 Char_t *dEdxPrFile="dEdxPr.root";
1468 Char_t *dEdxElFile="dEdxEl.root";
1469 Char_t *cmFile1 ="CovMatrix_AllEvts_1.root";
1470 Char_t *cmFile2 ="CovMatrix_AllEvts_2.root";
1471 Char_t *cmFile3 ="CovMatrix_AllEvts_3.root";
1473 // store the effieciencies
1475 WriteEffs(fDBfileName.Data());
1478 ReadPulls(pullsFile);
1479 WritePulls(fDBfileName.Data());
1481 //store the regularization parameters
1482 ReadRegParams(regPiFile,211);
1483 WriteRegParams(fDBfileName.Data(),211);
1484 ReadRegParams(regKaFile,321);
1485 WriteRegParams(fDBfileName.Data(),321);
1486 ReadRegParams(regElFile,11);
1487 WriteRegParams(fDBfileName.Data(),11);
1489 // store the dEdx parameters
1490 ReaddEdx(dEdxPiFile,211);
1491 WritedEdx(fDBfileName.Data(),211);
1492 ReaddEdx(dEdxKaFile,321);
1493 WritedEdx(fDBfileName.Data(),321);
1494 ReaddEdx(dEdxPrFile,2212);
1495 WritedEdx(fDBfileName.Data(),2212);
1496 ReaddEdx(dEdxElFile,11);
1497 WritedEdx(fDBfileName.Data(),11);
1501 // store the regularized covariance matrices
1503 InitializeKineGrid("DB","points");
1505 const Int_t nBinsPi = fDBgridPi.GetTotBins();
1506 const Int_t nBinsKa = fDBgridKa.GetTotBins();
1507 const Int_t nBinsEl = fDBgridEl.GetTotBins();
1510 // create the trees for cov. matrices
1512 TTree *CovTreePi_ = NULL;
1513 CovTreePi_ = new TTree[nBinsPi];
1515 TTree *CovTreeKa_ = NULL;
1516 CovTreeKa_ = new TTree[nBinsKa];
1517 // trees for electrons
1518 TTree *CovTreeEl_ = NULL;
1519 CovTreeEl_ = new TTree[nBinsEl];
1521 Char_t hname[100], htitle[100];
1525 for(Int_t i=0; i<nBinsPi; i++) {
1526 sprintf(hname,"CovTreePi_bin%d",i);
1527 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1528 CovTreePi_[i].SetName(hname); CovTreePi_[i].SetTitle(htitle);
1529 CovTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1531 for(Int_t i=0; i<nBinsKa; i++) {
1532 sprintf(hname,"CovTreeKa_bin%d",i);
1533 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1534 CovTreeKa_[i].SetName(hname); CovTreeKa_[i].SetTitle(htitle);
1535 CovTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1537 for(Int_t i=0; i<nBinsEl; i++) {
1538 sprintf(hname,"CovTreeEl_bin%d",i);
1539 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1540 CovTreeEl_[i].SetName(hname); CovTreeEl_[i].SetTitle(htitle);
1541 CovTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1544 // create the chain with the compared tracks
1545 TChain cmptrkchain("cmptrktree");
1546 cmptrkchain.Add(cmFile1);
1547 cmptrkchain.Add(cmFile2);
1548 cmptrkchain.Add(cmFile3);
1551 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
1552 Int_t entries = (Int_t)cmptrkchain.GetEntries();
1553 cerr<<" Number of entries: "<<entries<<endl;
1555 Int_t trkPdg,trkBin;
1556 Double_t trkKine[2],trkRegPar[4];
1557 Int_t * nPerBinPi = new Int_t[nBinsPi];
1558 for(Int_t k=0;k<nBinsPi;k++) nPerBinPi[k]=0;
1559 Int_t * nPerBinKa = new Int_t[nBinsKa];
1560 for(Int_t k=0;k<nBinsKa;k++) nPerBinKa[k]=0;
1561 Int_t * nPerBinEl = new Int_t[nBinsEl];
1562 for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
1564 // loop on chain entries
1565 for(Int_t l=0; l<entries; l++) {
1566 if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1568 cmptrkchain.GetEvent(l);
1570 trkPdg = TMath::Abs(cmptrk.pdg);
1571 // use only pions, kaons, electrons
1572 if(trkPdg!=211 && trkPdg!=321 && trkPdg!=11) continue;
1573 SetParticle(trkPdg);
1574 trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1575 //cerr<<cmptrk.pt<<" "<<cmptrk.eta<<" "<<trkBin<<endl;
1577 if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1578 if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1579 if(trkPdg==11 && nPerBinEl[trkBin]>=5000) continue;
1581 trkKine[0] = cmptrk.p;
1582 trkKine[1] = cmptrk.cosl;
1584 // get regularized covariance matrix
1585 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(0,k);
1586 covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1587 covmat.c10 = cmptrk.c10;
1588 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(1,k);
1589 covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1590 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(2,k);
1591 covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1592 covmat.c21 = cmptrk.c21;
1593 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(3,k);
1594 covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
1595 covmat.c30 = cmptrk.c30;
1596 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(4,k);
1597 covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
1598 covmat.c32 = cmptrk.c32;
1599 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(5,k);
1600 covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
1601 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(6,k);
1602 covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
1603 covmat.c41 = cmptrk.c41;
1604 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(7,k);
1605 covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
1606 covmat.c43 = cmptrk.c43;
1607 for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(8,k);
1608 covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
1613 CovTreePi_[trkBin].Fill();
1614 nPerBinPi[trkBin]++;
1617 CovTreeKa_[trkBin].Fill();
1618 nPerBinKa[trkBin]++;
1620 case 11: // electrons
1621 CovTreeEl_[trkBin].Fill();
1622 nPerBinEl[trkBin]++;
1625 } // loop on chain entries
1627 // store all trees the DB file
1628 TFile *DBfile = new TFile(fDBfileName.Data(),"update");
1629 DBfile->mkdir("CovMatrices");
1630 gDirectory->cd("/CovMatrices");
1631 gDirectory->mkdir("Pions");
1632 gDirectory->mkdir("Kaons");
1633 gDirectory->mkdir("Electrons");
1635 gDirectory->cd("/CovMatrices/Pions");
1636 fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
1637 for(Int_t i=0;i<nBinsPi;i++) CovTreePi_[i].Write();
1639 gDirectory->cd("/CovMatrices/Kaons");
1640 fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
1641 for(Int_t i=0;i<nBinsKa;i++) CovTreeKa_[i].Write();
1643 gDirectory->cd("/CovMatrices/Electrons");
1644 fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
1645 for(Int_t i=0;i<nBinsEl;i++) CovTreeEl_[i].Write();
1648 delete [] nPerBinPi;
1649 delete [] nPerBinKa;
1650 delete [] nPerBinEl;
1654 //-----------------------------------------------------------------------------
1655 void AliTPCtrackerParam::MergeEvents(Int_t evFirst,Int_t evLast) {
1656 //-----------------------------------------------------------------------------
1657 // This function: 1) merges the files from track comparison
1658 // (beware: better no more than 100 events per file)
1659 // 2) computes the average TPC efficiencies
1660 //-----------------------------------------------------------------------------
1662 Char_t *outName="TPCeff.root";
1665 Float_t effPi,effKa,effPr,effEl,effMu;
1666 Float_t avEffPi[27],avEffKa[27],avEffPr[27],avEffEl[27],avEffMu[27];
1667 Int_t evtsPi[27],evtsKa[27],evtsPr[27],evtsEl[27],evtsMu[27];
1668 for(Int_t j=0; j<27; j++) {
1669 avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
1670 evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
1673 // create the chain for the tree of compared tracks
1674 TChain ch1("cmptrktree");
1675 TChain ch2("cmptrktree");
1676 TChain ch3("cmptrktree");
1679 for(Int_t j=evFirst; j<=evLast; j++) {
1680 cerr<<"Processing event "<<j<<endl;
1682 TString covName("CovMatrix.");
1683 TString effName("TPCeff.");
1686 covName.Append(".root");
1687 effName.Append(".dat");
1689 if(gSystem->AccessPathName(covName.Data(),kFileExists) ||
1690 gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
1692 if(j<=100) ch1.Add(covName.Data());
1693 if(j>100 && j<=200) ch2.Add(covName.Data());
1694 if(j>200) ch3.Add(covName.Data());
1697 FILE *effIn = fopen(effName.Data(),"r");
1698 for(Int_t k=0; k<27; k++) {
1699 nCol = fscanf(effIn,"%f %f %f %f %f",&effPi,&effKa,&effPr,&effEl,&effMu);
1700 if(effPi>=0.) {avEffPi[k]+=effPi; evtsPi[k]++;}
1701 if(effKa>=0.) {avEffKa[k]+=effKa; evtsKa[k]++;}
1702 if(effPr>=0.) {avEffPr[k]+=effPr; evtsPr[k]++;}
1703 if(effEl>=0.) {avEffEl[k]+=effEl; evtsEl[k]++;}
1704 if(effMu>=0.) {avEffMu[k]+=effMu; evtsMu[k]++;}
1710 // merge chain in one file
1712 covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
1713 ch1.Merge(covOut,1000000000);
1716 covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
1717 ch2.Merge(covOut,1000000000);
1720 covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
1721 ch3.Merge(covOut,1000000000);
1728 InitializeKineGrid("eff","bins");
1731 for(Int_t j=0; j<27; j++) {
1732 if(evtsPi[j]==0) evtsPi[j]++;
1733 fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
1738 for(Int_t k=0; k<27; k++) {
1739 if(k>14 && k!=21 && k!=22 && k!=23) continue;
1740 if(evtsKa[k]==0) evtsKa[k]++;
1741 fEffKa.SetParam(j,(Double_t)avEffKa[k]/evtsKa[k]);
1746 for(Int_t j=0; j<15; j++) {
1747 if(evtsPr[j]==0) evtsPr[j]++;
1748 fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
1752 for(Int_t j=0; j<27; j++) {
1753 if(evtsEl[j]==0) evtsEl[j]++;
1754 fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
1758 for(Int_t j=0; j<9; j++) {
1759 if(evtsMu[j]==0) evtsMu[j]++;
1760 fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
1763 // write efficiencies to a file
1768 //-----------------------------------------------------------------------------
1769 Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
1770 //-----------------------------------------------------------------------------
1771 // This function reads all parameters from the DB
1772 //-----------------------------------------------------------------------------
1774 if(!ReadEffs(inName)) return 0;
1775 if(!ReadPulls(inName)) return 0;
1776 if(!ReadRegParams(inName,211)) return 0;
1777 if(!ReadRegParams(inName,321)) return 0;
1778 if(!ReadRegParams(inName,11)) return 0;
1779 if(!ReaddEdx(inName,211)) return 0;
1780 if(!ReaddEdx(inName,321)) return 0;
1781 if(!ReaddEdx(inName,2212)) return 0;
1782 if(!ReaddEdx(inName,11)) return 0;
1783 if(!ReadDBgrid(inName)) return 0;
1787 //-----------------------------------------------------------------------------
1788 Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
1789 //-----------------------------------------------------------------------------
1790 // This function reads the dEdx parameters from the DB
1791 //-----------------------------------------------------------------------------
1793 if(gSystem->AccessPathName(inName,kFileExists)) {
1794 cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
1797 TFile *inFile = TFile::Open(inName);
1800 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
1801 fdEdxMeanPi.~AliTPCkineGrid();
1802 new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
1803 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
1804 fdEdxRMSPi.~AliTPCkineGrid();
1805 new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
1808 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
1809 fdEdxMeanKa.~AliTPCkineGrid();
1810 new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
1811 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
1812 fdEdxRMSKa.~AliTPCkineGrid();
1813 new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
1816 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
1817 fdEdxMeanPr.~AliTPCkineGrid();
1818 new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
1819 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
1820 fdEdxRMSPr.~AliTPCkineGrid();
1821 new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
1824 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
1825 fdEdxMeanEl.~AliTPCkineGrid();
1826 new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
1827 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
1828 fdEdxRMSEl.~AliTPCkineGrid();
1829 new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
1836 //-----------------------------------------------------------------------------
1837 Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
1838 //-----------------------------------------------------------------------------
1839 // This function reads the kine grid from the DB
1840 //-----------------------------------------------------------------------------
1842 if(gSystem->AccessPathName(inName,kFileExists)) {
1843 cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n";
1846 TFile *inFile = TFile::Open(inName);
1848 // first read the DB grid for the different particles
1849 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
1850 fDBgridPi.~AliTPCkineGrid();
1851 new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
1852 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
1853 fDBgridKa.~AliTPCkineGrid();
1854 new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
1855 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
1856 fDBgridEl.~AliTPCkineGrid();
1857 new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
1863 //-----------------------------------------------------------------------------
1864 Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
1865 //-----------------------------------------------------------------------------
1866 // This function reads the TPC efficiencies from the DB
1867 //-----------------------------------------------------------------------------
1869 if(gSystem->AccessPathName(inName,kFileExists)) {
1870 cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n";
1873 TFile *inFile = TFile::Open(inName);
1875 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
1876 fEffPi.~AliTPCkineGrid();
1877 new(&fEffPi) AliTPCkineGrid(*fEff);
1878 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
1879 fEffKa.~AliTPCkineGrid();
1880 new(&fEffKa) AliTPCkineGrid(*fEff);
1881 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
1882 fEffPr.~AliTPCkineGrid();
1883 new(&fEffPr) AliTPCkineGrid(*fEff);
1884 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
1885 fEffEl.~AliTPCkineGrid();
1886 new(&fEffEl) AliTPCkineGrid(*fEff);
1887 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
1888 fEffMu.~AliTPCkineGrid();
1889 new(&fEffMu) AliTPCkineGrid(*fEff);
1895 //-----------------------------------------------------------------------------
1896 Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
1897 //-----------------------------------------------------------------------------
1898 // This function reads the pulls from the DB
1899 //-----------------------------------------------------------------------------
1901 if(gSystem->AccessPathName(inName,kFileExists)) {
1902 cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n";
1905 TFile *inFile = TFile::Open(inName);
1907 for(Int_t i=0; i<5; i++) {
1908 TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
1909 TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
1910 TString el("/Pulls/Electrons/PullsEl_"); el+=i;
1911 fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
1912 fPullsPi[i].~AliTPCkineGrid();
1913 new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
1914 fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
1915 fPullsKa[i].~AliTPCkineGrid();
1916 new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
1917 fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
1918 fPullsEl[i].~AliTPCkineGrid();
1919 new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
1926 //-----------------------------------------------------------------------------
1927 Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
1928 //-----------------------------------------------------------------------------
1929 // This function reads the regularization parameters from the DB
1930 //-----------------------------------------------------------------------------
1932 if(gSystem->AccessPathName(inName,kFileExists)) {
1933 cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
1936 TFile *inFile = TFile::Open(inName);
1939 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
1940 new(&fRegParPi) TMatrixD(*fRegPar);
1943 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
1944 new(&fRegParKa) TMatrixD(*fRegPar);
1947 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
1948 new(&fRegParEl) TMatrixD(*fRegPar);
1955 //-----------------------------------------------------------------------------
1956 void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
1957 //-----------------------------------------------------------------------------
1958 // This function regularizes the elements of the covariance matrix
1959 // that show a momentum depence:
1960 // c00, c11, c22, c33, c44, c20, c24, c40, c31
1961 // The regularization is done separately for pions, kaons, electrons:
1962 // give "Pion","Kaon" or "Electron" as first argument.
1963 //-----------------------------------------------------------------------------
1965 gStyle->SetOptStat(0);
1966 gStyle->SetOptFit(10001);
1969 Char_t *part="Pions - ";
1971 InitializeKineGrid("DB","points");
1973 const Int_t fitbins = fDBgrid->GetBinsPt();
1974 cerr<<" Fit bins: "<<fitbins<<endl;
1980 cerr<<" Processing pions ...\n";
1985 cerr<<" Processing kaons ...\n";
1987 case 11: // electrons
1989 part="Electrons - ";
1990 cerr<<" Processing electrons ...\n";
1994 // create the chain with all compared tracks
1995 TChain cmptrkchain("cmptrktree");
1996 cmptrkchain.Add("CovMatrix_AllEvts_1.root");
1997 cmptrkchain.Add("CovMatrix_AllEvts_2.root");
1998 cmptrkchain.Add("CovMatrix_AllEvts_3.root");
2001 cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
2002 Int_t entries = (Int_t)cmptrkchain.GetEntries();
2006 Int_t * n = new Int_t[fitbins];
2007 Int_t * n00 = new Int_t[fitbins];
2008 Int_t * n11 = new Int_t[fitbins];
2009 Int_t * n20 = new Int_t[fitbins];
2010 Int_t * n22 = new Int_t[fitbins];
2011 Int_t * n31 = new Int_t[fitbins];
2012 Int_t * n33 = new Int_t[fitbins];
2013 Int_t * n40 = new Int_t[fitbins];
2014 Int_t * n42 = new Int_t[fitbins];
2015 Int_t * n44 = new Int_t[fitbins];
2016 Double_t * p = new Double_t[fitbins];
2017 Double_t * ep = new Double_t[fitbins];
2018 Double_t * mean00 = new Double_t[fitbins];
2019 Double_t * mean11 = new Double_t[fitbins];
2020 Double_t * mean20 = new Double_t[fitbins];
2021 Double_t * mean22 = new Double_t[fitbins];
2022 Double_t * mean31 = new Double_t[fitbins];
2023 Double_t * mean33 = new Double_t[fitbins];
2024 Double_t * mean40 = new Double_t[fitbins];
2025 Double_t * mean42 = new Double_t[fitbins];
2026 Double_t * mean44 = new Double_t[fitbins];
2027 Double_t * sigma00 = new Double_t[fitbins];
2028 Double_t * sigma11 = new Double_t[fitbins];
2029 Double_t * sigma20 = new Double_t[fitbins];
2030 Double_t * sigma22 = new Double_t[fitbins];
2031 Double_t * sigma31 = new Double_t[fitbins];
2032 Double_t * sigma33 = new Double_t[fitbins];
2033 Double_t * sigma40 = new Double_t[fitbins];
2034 Double_t * sigma42 = new Double_t[fitbins];
2035 Double_t * sigma44 = new Double_t[fitbins];
2036 Double_t * rmean = new Double_t[fitbins];
2037 Double_t * rsigma = new Double_t[fitbins];
2040 for(Int_t l=0; l<fitbins; l++) {
2042 n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
2044 mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
2045 sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
2048 // loop on chain entries for mean
2049 for(Int_t l=0; l<entries; l++) {
2050 cmptrkchain.GetEvent(l);
2051 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2052 pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2055 mean00[pbin]+=cmptrk.c00;
2056 mean11[pbin]+=cmptrk.c11;
2057 mean20[pbin]+=cmptrk.c20;
2058 mean22[pbin]+=cmptrk.c22;
2059 mean31[pbin]+=cmptrk.c31;
2060 mean33[pbin]+=cmptrk.c33;
2061 mean40[pbin]+=cmptrk.c40;
2062 mean42[pbin]+=cmptrk.c42;
2063 mean44[pbin]+=cmptrk.c44;
2064 } // loop on chain entries
2066 for(Int_t l=0; l<fitbins; l++) {
2079 // loop on chain entries for sigma
2080 for(Int_t l=0; l<entries; l++) {
2081 cmptrkchain.GetEvent(l);
2082 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2083 pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2084 if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2085 sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2086 if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2087 sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2088 if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2089 sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2090 if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2091 sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2092 if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2093 sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2094 if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2095 sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2096 if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2097 sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2098 if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2099 sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2100 if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2101 sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2102 } // loop on chain entries
2104 for(Int_t l=0; l<fitbins; l++) {
2105 sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2106 sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2107 sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2108 sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2109 sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2110 sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2111 sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2112 sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2113 sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2118 TF1 *func = new TF1("FitRegFunc",FitRegFunc,0.23,50.,4);
2119 func->SetParNames("A_meas","A_scatt","B1","B2");
2121 // line to draw on the plots
2122 TLine *lin = new TLine(-1,1,1.69,1);
2123 lin->SetLineStyle(2);
2124 lin->SetLineWidth(2);
2126 // matrix used to store fit results
2127 TMatrixD fitRes(9,4);
2131 // create the canvas
2132 TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
2133 canv00->Divide(1,2);
2134 // create the graph for cov matrix
2135 TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
2136 TString title00("C(y,y)"); title00.Prepend(part);
2137 TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
2138 frame00->SetXTitle("p [GeV/c]");
2139 canv00->cd(1); gPad->SetLogx();
2142 // Sets initial values for parameters
2143 func->SetParameters(1.6e-3,1.9e-4,1.5,0.);
2144 // Fit points in range defined by function
2145 gr00->Fit("FitRegFunc","R,Q");
2146 func->GetParameters(fitpar);
2147 for(Int_t i=0; i<4; i++) fitRes(0,i)=fitpar[i];
2148 for(Int_t l=0; l<fitbins; l++) {
2149 rmean[l] = mean00[l]/FitRegFunc(&p[l],fitpar);
2150 rsigma[l] = sigma00[l]/FitRegFunc(&p[l],fitpar);
2152 // create the graph the regularized cov. matrix
2153 TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2154 TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2155 regtitle00.Prepend(part);
2156 TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
2157 frame00reg->SetXTitle("p [GeV/c]");
2158 canv00->cd(2); gPad->SetLogx();
2166 // create the canvas
2167 TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
2168 canv11->Divide(1,2);
2169 // create the graph for cov matrix
2170 TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
2171 TString title11("C(z,z)"); title11.Prepend(part);
2172 TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
2173 frame11->SetXTitle("p [GeV/c]");
2174 canv11->cd(1); gPad->SetLogx();
2177 // Sets initial values for parameters
2178 func->SetParameters(1.2e-3,8.1e-4,1.,0.);
2179 // Fit points in range defined by function
2180 gr11->Fit("FitRegFunc","R,Q");
2181 func->GetParameters(fitpar);
2182 for(Int_t i=0; i<4; i++) fitRes(1,i)=fitpar[i];
2183 for(Int_t l=0; l<fitbins; l++) {
2184 rmean[l] = mean11[l]/FitRegFunc(&p[l],fitpar);
2185 rsigma[l] = sigma11[l]/FitRegFunc(&p[l],fitpar);
2187 // create the graph the regularized cov. matrix
2188 TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2189 TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2190 regtitle11.Prepend(part);
2191 TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
2192 frame11reg->SetXTitle("p [GeV/c]");
2193 canv11->cd(2); gPad->SetLogx();
2201 // create the canvas
2202 TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
2203 canv20->Divide(1,2);
2204 // create the graph for cov matrix
2205 TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
2206 TString title20("C(#eta, y)"); title20.Prepend(part);
2207 TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
2208 frame20->SetXTitle("p [GeV/c]");
2209 canv20->cd(1); gPad->SetLogx();
2212 // Sets initial values for parameters
2213 func->SetParameters(7.3e-5,1.2e-5,1.5,0.);
2214 // Fit points in range defined by function
2215 gr20->Fit("FitRegFunc","R,Q");
2216 func->GetParameters(fitpar);
2217 for(Int_t i=0; i<4; i++) fitRes(2,i)=fitpar[i];
2218 for(Int_t l=0; l<fitbins; l++) {
2219 rmean[l] = mean20[l]/FitRegFunc(&p[l],fitpar);
2220 rsigma[l] = sigma20[l]/FitRegFunc(&p[l],fitpar);
2222 // create the graph the regularized cov. matrix
2223 TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2224 TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2225 regtitle20.Prepend(part);
2226 TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
2227 frame20reg->SetXTitle("p [GeV/c]");
2228 canv20->cd(2); gPad->SetLogx();
2236 // create the canvas
2237 TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
2238 canv22->Divide(1,2);
2239 // create the graph for cov matrix
2240 TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
2241 TString title22("C(#eta, #eta)"); title22.Prepend(part);
2242 TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
2243 frame22->SetXTitle("p [GeV/c]");
2244 canv22->cd(1); gPad->SetLogx();
2247 // Sets initial values for parameters
2248 func->SetParameters(5.2e-6,1.1e-6,2.,1.);
2249 // Fit points in range defined by function
2250 gr22->Fit("FitRegFunc","R,Q");
2251 func->GetParameters(fitpar);
2252 for(Int_t i=0; i<4; i++) fitRes(3,i)=fitpar[i];
2253 for(Int_t l=0; l<fitbins; l++) {
2254 rmean[l] = mean22[l]/FitRegFunc(&p[l],fitpar);
2255 rsigma[l] = sigma22[l]/FitRegFunc(&p[l],fitpar);
2257 // create the graph the regularized cov. matrix
2258 TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2259 TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2260 regtitle22.Prepend(part);
2261 TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
2262 frame22reg->SetXTitle("p [GeV/c]");
2263 canv22->cd(2); gPad->SetLogx();
2271 // create the canvas
2272 TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
2273 canv31->Divide(1,2);
2274 // create the graph for cov matrix
2275 TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
2276 TString title31("C(tg #lambda,z)"); title31.Prepend(part);
2277 TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
2278 frame31->SetXTitle("p [GeV/c]");
2279 canv31->cd(1); gPad->SetLogx();
2282 // Sets initial values for parameters
2283 func->SetParameters(-1.2e-5,-1.2e-5,1.5,3.);
2284 // Fit points in range defined by function
2285 gr31->Fit("FitRegFunc","R,Q");
2286 func->GetParameters(fitpar);
2287 for(Int_t i=0; i<4; i++) fitRes(4,i)=fitpar[i];
2288 for(Int_t l=0; l<fitbins; l++) {
2289 rmean[l] = mean31[l]/FitRegFunc(&p[l],fitpar);
2290 rsigma[l] = -sigma31[l]/FitRegFunc(&p[l],fitpar);
2292 // create the graph the regularized cov. matrix
2293 TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2294 TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2295 regtitle31.Prepend(part);
2296 TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
2297 frame31reg->SetXTitle("p [GeV/c]");
2298 canv31->cd(2); gPad->SetLogx();
2306 // create the canvas
2307 TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
2308 canv33->Divide(1,2);
2309 // create the graph for cov matrix
2310 TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
2311 TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
2312 TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
2313 frame33->SetXTitle("p [GeV/c]");
2314 canv33->cd(1); gPad->SetLogx();
2317 // Sets initial values for parameters
2318 func->SetParameters(1.3e-7,4.6e-7,1.7,4.);
2319 // Fit points in range defined by function
2320 gr33->Fit("FitRegFunc","R,Q");
2321 func->GetParameters(fitpar);
2322 for(Int_t i=0; i<4; i++) fitRes(5,i)=fitpar[i];
2323 for(Int_t l=0; l<fitbins; l++) {
2324 rmean[l] = mean33[l]/FitRegFunc(&p[l],fitpar);
2325 rsigma[l] = sigma33[l]/FitRegFunc(&p[l],fitpar);
2327 // create the graph the regularized cov. matrix
2328 TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2329 TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2330 regtitle33.Prepend(part);
2331 TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
2332 frame33reg->SetXTitle("p [GeV/c]");
2333 canv33->cd(2); gPad->SetLogx();
2341 // create the canvas
2342 TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
2343 canv40->Divide(1,2);
2344 // create the graph for cov matrix
2345 TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
2346 TString title40("C(C,y)"); title40.Prepend(part);
2347 TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
2348 frame40->SetXTitle("p [GeV/c]");
2349 canv40->cd(1); gPad->SetLogx();
2352 // Sets initial values for parameters
2353 func->SetParameters(4.e-7,4.4e-8,1.5,0.);
2354 // Fit points in range defined by function
2355 gr40->Fit("FitRegFunc","R,Q");
2356 func->GetParameters(fitpar);
2357 for(Int_t i=0; i<4; i++) fitRes(6,i)=fitpar[i];
2358 for(Int_t l=0; l<fitbins; l++) {
2359 rmean[l] = mean40[l]/FitRegFunc(&p[l],fitpar);
2360 rsigma[l] = sigma40[l]/FitRegFunc(&p[l],fitpar);
2362 // create the graph the regularized cov. matrix
2363 TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2364 TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2365 regtitle40.Prepend(part);
2366 TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
2367 frame40reg->SetXTitle("p [GeV/c]");
2368 canv40->cd(2); gPad->SetLogx();
2376 // create the canvas
2377 TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
2378 canv42->Divide(1,2);
2379 // create the graph for cov matrix
2380 TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
2381 TString title42("C(C, #eta)"); title42.Prepend(part);
2382 TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
2383 frame42->SetXTitle("p [GeV/c]");
2384 canv42->cd(1); gPad->SetLogx();
2387 // Sets initial values for parameters
2388 func->SetParameters(3.e-8,8.2e-9,2.,1.);
2389 // Fit points in range defined by function
2390 gr42->Fit("FitRegFunc","R,Q");
2391 func->GetParameters(fitpar);
2392 for(Int_t i=0; i<4; i++) fitRes(7,i)=fitpar[i];
2393 for(Int_t l=0; l<fitbins; l++) {
2394 rmean[l] = mean42[l]/FitRegFunc(&p[l],fitpar);
2395 rsigma[l] = sigma42[l]/FitRegFunc(&p[l],fitpar);
2397 // create the graph the regularized cov. matrix
2398 TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2399 TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2400 regtitle42.Prepend(part);
2401 TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
2402 frame42reg->SetXTitle("p [GeV/c]");
2403 canv42->cd(2); gPad->SetLogx();
2411 // create the canvas
2412 TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
2413 canv44->Divide(1,2);
2414 // create the graph for cov matrix
2415 TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
2416 TString title44("C(C,C)"); title44.Prepend(part);
2417 TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
2418 frame44->SetXTitle("p [GeV/c]");
2419 canv44->cd(1); gPad->SetLogx();
2422 // Sets initial values for parameters
2423 func->SetParameters(1.8e-10,5.8e-11,2.,3.);
2424 // Fit points in range defined by function
2425 gr44->Fit("FitRegFunc","R,Q");
2426 func->GetParameters(fitpar);
2427 for(Int_t i=0; i<4; i++) fitRes(8,i)=fitpar[i];
2428 for(Int_t l=0; l<fitbins; l++) {
2429 rmean[l] = mean44[l]/FitRegFunc(&p[l],fitpar);
2430 rsigma[l] = sigma44[l]/FitRegFunc(&p[l],fitpar);
2432 // create the graph the regularized cov. matrix
2433 TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2434 TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2435 regtitle44.Prepend(part);
2436 TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
2437 frame44reg->SetXTitle("p [GeV/c]");
2438 canv44->cd(2); gPad->SetLogx();
2448 new(&fRegParPi) TMatrixD(fitRes);
2451 new(&fRegParKa) TMatrixD(fitRes);
2454 new(&fRegParEl) TMatrixD(fitRes);
2458 // write fit parameters to file
2459 WriteRegParams(outName,pdg);
2496 //-----------------------------------------------------------------------------
2497 Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
2498 //-----------------------------------------------------------------------------
2499 // This function makes a selection according to TPC tracking efficiency
2500 //-----------------------------------------------------------------------------
2504 eff = fEff->GetValueAt(pt,eta);
2506 if(gRandom->Rndm() < eff) return kTRUE;
2510 //-----------------------------------------------------------------------------
2511 void AliTPCtrackerParam::SetParticle(Int_t pdg) {
2512 //-----------------------------------------------------------------------------
2513 // This function set efficiencies, pulls, etc... for the particle
2514 // specie of the current particle
2515 //-----------------------------------------------------------------------------
2519 fDBgrid = &fDBgridPi;
2522 fRegPar = &fRegParPi;
2523 fdEdxMean = &fdEdxMeanPi;
2524 fdEdxRMS = &fdEdxRMSPi;
2527 fDBgrid = &fDBgridKa;
2530 fRegPar = &fRegParKa;
2531 fdEdxMean = &fdEdxMeanKa;
2532 fdEdxRMS = &fdEdxRMSKa;
2535 fDBgrid = &fDBgridPi;
2538 fRegPar = &fRegParPi;
2539 fdEdxMean = &fdEdxMeanPr;
2540 fdEdxRMS = &fdEdxRMSPr;
2543 fDBgrid = &fDBgridEl;
2546 fRegPar = &fRegParEl;
2547 fdEdxMean = &fdEdxMeanEl;
2548 fdEdxRMS = &fdEdxRMSEl;
2551 fDBgrid = &fDBgridPi;
2554 fRegPar = &fRegParPi;
2555 fdEdxMean = &fdEdxMeanPi;
2556 fdEdxRMS = &fdEdxRMSPi;
2562 //-----------------------------------------------------------------------------
2563 void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
2565 //-----------------------------------------------------------------------------
2566 // This function smears track parameters according to streched cov. matrix
2567 //-----------------------------------------------------------------------------
2568 AliGausCorr *corgen = new AliGausCorr(cov,5);
2570 corgen->GetGaussN(corr);
2574 for(Int_t l=0;l<5;l++) {
2575 xxsm[l] = xx[l]+corr[l];
2580 //-----------------------------------------------------------------------------
2581 Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
2582 //-----------------------------------------------------------------------------
2583 // This function writes the dEdx parameters to the DB
2584 //-----------------------------------------------------------------------------
2587 Char_t *dirName="Pions";
2588 Char_t *meanName="dEdxMeanPi";
2589 Char_t *RMSName="dEdxRMSPi";
2593 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2594 } else { opt="update"; }
2599 meanName="dEdxMeanPi";
2600 RMSName="dEdxRMSPi";
2604 meanName="dEdxMeanKa";
2605 RMSName="dEdxRMSKa";
2609 meanName="dEdxMeanPr";
2610 RMSName="dEdxRMSPr";
2613 dirName="Electrons";
2614 meanName="dEdxMeanEl";
2615 RMSName="dEdxRMSEl";
2619 TFile *outFile = new TFile(outName,opt);
2620 if(!gDirectory->cd("/dEdx")) {
2621 outFile->mkdir("dEdx");
2622 gDirectory->cd("/dEdx");
2624 TDirectory *dir2 = gDirectory->mkdir(dirName);
2626 fdEdxMean->SetName(meanName); fdEdxMean->Write();
2627 fdEdxRMS->SetName(RMSName); fdEdxRMS->Write();
2635 //-----------------------------------------------------------------------------
2636 Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
2637 //-----------------------------------------------------------------------------
2638 // This function writes the TPC efficiencies to the DB
2639 //-----------------------------------------------------------------------------
2644 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2645 } else { opt="update"; }
2647 TFile *outFile = new TFile(outName,opt);
2649 outFile->mkdir("Efficiencies");
2650 gDirectory->cd("/Efficiencies");
2651 gDirectory->mkdir("Pions");
2652 gDirectory->mkdir("Kaons");
2653 gDirectory->mkdir("Protons");
2654 gDirectory->mkdir("Electrons");
2655 gDirectory->mkdir("Muons");
2657 gDirectory->cd("/Efficiencies/Pions");
2658 fEffPi.SetName("EffPi");
2660 gDirectory->cd("/Efficiencies/Kaons");
2661 fEffKa.SetName("EffKa");
2663 gDirectory->cd("/Efficiencies/Protons");
2664 fEffPr.SetName("EffPr");
2666 gDirectory->cd("/Efficiencies/Electrons");
2667 fEffEl.SetName("EffEl");
2669 gDirectory->cd("/Efficiencies/Muons");
2670 fEffMu.SetName("EffMu");
2679 //-----------------------------------------------------------------------------
2680 Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
2681 //-----------------------------------------------------------------------------
2682 // This function writes the pulls to the DB
2683 //-----------------------------------------------------------------------------
2687 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2688 } else { opt="update"; }
2690 TFile *outFile = new TFile(outName,opt);
2692 outFile->mkdir("Pulls");
2693 gDirectory->cd("/Pulls");
2694 gDirectory->mkdir("Pions");
2695 gDirectory->mkdir("Kaons");
2696 gDirectory->mkdir("Electrons");
2698 for(Int_t i=0;i<5;i++) {
2699 TString pi("PullsPi_"); pi+=i;
2700 TString ka("PullsKa_"); ka+=i;
2701 TString el("PullsEl_"); el+=i;
2702 fPullsPi[i].SetName(pi.Data());
2703 fPullsKa[i].SetName(ka.Data());
2704 fPullsEl[i].SetName(el.Data());
2705 gDirectory->cd("/Pulls/Pions");
2706 fPullsPi[i].Write();
2707 gDirectory->cd("/Pulls/Kaons");
2708 fPullsKa[i].Write();
2709 gDirectory->cd("/Pulls/Electrons");
2710 fPullsEl[i].Write();
2717 //-----------------------------------------------------------------------------
2718 Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
2719 //-----------------------------------------------------------------------------
2720 // This function writes the regularization parameters to the DB
2721 //-----------------------------------------------------------------------------
2724 Char_t *dirName="Pions";
2725 Char_t *keyName="RegPions";
2729 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
2730 } else { opt="update"; }
2742 dirName="Electrons";
2743 keyName="RegElectrons";
2747 TFile *outFile = new TFile(outName,opt);
2748 if(!gDirectory->cd("/RegParams")) {
2749 outFile->mkdir("RegParams");
2750 gDirectory->cd("/RegParams");
2752 TDirectory *dir2 = gDirectory->mkdir(dirName);
2754 fRegPar->Write(keyName);