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 /**************************************************************************
20 * This class builds AliTPCtrack objects from generated tracks to feed *
21 * ITS tracking (V2). The AliTPCtrack is built from its first hit in *
22 * the TPC. The track is assigned a Kalman-like covariance matrix *
23 * depending on its pT and pseudorapidity and track parameters are *
24 * smeared according to this covariance matrix. *
25 * Output file contains sorted tracks, ready for matching with ITS. *
28 * Alice Internal Note 2003-011 *
30 * Test macro is: AliBarrelRec_TPCparam.C *
32 * 2002/10/01: Introduction of the database for pp collisions (B=0.4 T) *
33 * - Everything done separately for pions, kaons, protons, electrons and *
35 * - Now (only for pp) the tracks are built from the AliTrackReferences *
36 * which contain the position and momentum of all tracks at R = 83 cm; *
37 * This eliminates the loss of tracks due the dead zone of the TPC *
38 * where the 1st hit is not created. *
39 * - In AliBarrelRec_TPCparam.C there many possible ways of determining *
40 * the z coordinate of the primary vertex in pp events (from pixels, *
41 * from ITS tracks, smearing according to resolution given by tracks. *
43 * 2002/04/28: Major upgrade of the class *
44 * - Covariance matrices and pulls are now separeted for pions, kaons *
46 * - A parameterization of dE/dx in the TPC has been included and it is *
47 * used to assign a mass to each track according to a rough dE/dx PID. *
48 * - All the "numbers" have been moved to the file with the DataBase, *
49 * they are read as objects of the class AliTPCkineGrid, and assigned *
50 * to data memebers of the class AliTPCtrackerParam. *
51 * - All the code necessary to create a BataBase has been included in *
52 * class (see the macro AliTPCtrackingParamDB.C for the details). *
54 * 2006/03/16: Adapted to ESD input/output *
56 * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it *
57 * adapted to ESD output by Marcello Lunardon, Padova *
58 **************************************************************************/
60 // This is a dummy comment
64 //------- Root headers --------
65 #include <Riostream.h>
70 #include <TGraphErrors.h>
75 #include <TParticle.h>
80 //------ AliRoot headers ------
81 #include "AliGausCorr.h"
82 #include "AliTracker.h"
86 #include "AliRunLoader.h"
88 #include "AliTPCParamSR.h"
89 #include "AliTPCkineGrid.h"
90 #include "AliTPCtrack.h"
91 #include "AliTPCtrackerParam.h"
92 #include "AliTrackReference.h"
93 #include "AliESDtrack.h"
94 //-----------------------------
96 Double_t RegFunc(Double_t *x,Double_t *par) {
97 // This is the function used to regularize the covariance matrix
98 Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2]);
102 // structure for DB building
112 Double_t dP0,dP1,dP2,dP3,dP4;
113 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
116 // cov matrix structure
118 Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
121 ClassImp(AliTPCtrackerParam)
123 //-----------------------------------------------------------------------------
124 AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz,
125 const char* evfoldname):TObject(),
126 fEvFolderName(evfoldname),
165 //-----------------------------------------------------------------------------
166 // This is the class conctructor
167 //-----------------------------------------------------------------------------
169 // fBz = kBz; // value of the z component of L3 field (Tesla)
170 // fColl = kcoll; // collision code (0: PbPb6000; 1: pp)
171 // fSelAndSmear = kTRUE; // by default selection and smearing are done
173 if(fBz!=0.4 && fBz!=0.5) {
174 Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n Available: 0.4 or 0.5");
176 if(fColl!=0 && fColl!=1) {
177 Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n Available: 0 -> PbPb6000\n 1 -> pp");
180 fDBfileName = gSystem->Getenv("ALICE_ROOT");
181 fDBfileName.Append("/TPC/CovMatrixDB_");
182 //fDBfileName = "CovMatrixDB_";
183 if(fColl==0) fDBfileName.Append("PbPb6000");
184 if(fColl==1) fDBfileName.Append("pp");
185 if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
186 // use same DB for 0.4 and 0.5 T; for 0.5 T, correction done in CookTrack()
187 if(fBz==0.5) fDBfileName.Append("_B0.4T.root");
189 //-----------------------------------------------------------------------------
190 AliTPCtrackerParam::~AliTPCtrackerParam() {}
191 //____________________________________________________________________________
192 AliTPCtrackerParam::AliTPCtrackerParam( const AliTPCtrackerParam& p)
233 // dummy copy constructor
235 //----------------------------------------------------------------------------
236 AliTPCtrackerParam::AliTPCseedGeant::AliTPCseedGeant(
237 Double_t x,Double_t y,Double_t z,
238 Double_t px,Double_t py,Double_t pz,
252 //----------------------------------------------------------------------------
253 // Constructor of the geant seeds
254 //----------------------------------------------------------------------------
256 Double_t a = TMath::ATan2(y,x)*180./TMath::Pi();
258 fSector = (Int_t)(a/20.);
259 fAlpha = 10.+20.*fSector;
261 fAlpha *= TMath::Pi();
263 //-----------------------------------------------------------------------------
264 Int_t AliTPCtrackerParam::Init() {
265 //-----------------------------------------------------------------------------
266 // This function reads the DB from file
267 //-----------------------------------------------------------------------------
270 printf("+++\n+++ Reading DataBase from:%s\n+++\n+++\n",fDBfileName.Data());
271 // Read paramters from DB file
272 if(!ReadAllData(fDBfileName.Data())) {
273 printf("AliTPCtrackerParam::BuildTPCtracks: \
274 Could not read data from DB\n\n"); return 1;
277 } else printf("\n ! Creating ALL TRUE tracks at TPC inner radius !\n\n");
280 // Check if value in the galice file is equal to selected one (fBz)
281 AliMagF *fiel = (AliMagF*)gAlice->Field();
282 Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
283 printf("Magnetic field is %6.2f Tesla\n",fieval);
285 printf("AliTPCtrackerParam::BuildTPCtracks: Invalid field!");
286 printf("Field selected is: %f T\n",fBz);
287 printf("Field found on file is: %f T\n",fieval);
291 // Set the conversion constant between curvature and Pt
292 AliTracker::SetFieldMap(fiel,kTRUE);
296 //-----------------------------------------------------------------------------
297 Int_t AliTPCtrackerParam::BuildTPCtracks(AliESDEvent *event) {
298 //-----------------------------------------------------------------------------
299 // This function creates the TPC parameterized tracks and writes them
300 // as AliESDtrack objects in the ESD event
301 //-----------------------------------------------------------------------------
304 if(!event) return -1;
306 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
308 Error("BuildTPCtracks","Can not get Run Loader from event folder named %s",
309 fEvFolderName.Data());
312 AliLoader* tpcloader = rl->GetLoader("TPCLoader");
313 if (tpcloader == 0x0) {
314 Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
318 // Get gAlice object from file
319 if(!gAlice) rl->LoadgAlice();
321 rl->LoadKinematics();
322 tpcloader->LoadHits("read");
324 if(!(gAlice=rl->GetAliRun())) {
325 printf("Cannot get gAlice from Run Loader !\n");
330 AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC");
333 AliTPCParam *digp=(AliTPCParam*)gDirectory->Get("75x40_100x60");
336 digp = new AliTPCParamSR();
338 else digp=(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60");
340 if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
344 AliTPCseedGeant *seed=0;
345 AliTPCtrack *tpctrack=0;
347 Int_t bin,label,pdg,charge;
349 Int_t nParticles,nSeeds,arrentr;
350 //Int_t nSel=0,nAcc=0;
352 Int_t evt=event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
356 // array for TPC tracks
357 TObjArray tArray(20000);
359 // array for TPC seeds with geant info
360 TObjArray sArray(20000);
362 // get the particles stack
363 nParticles = (Int_t)gAlice->GetEvent(evt);
365 Bool_t *done = new Bool_t[nParticles];
366 Int_t *pdgCodes = new Int_t[nParticles];
368 // loop on particles and store pdg codes
369 for(Int_t l=0; l<nParticles; l++) {
370 part = (TParticle*)gAlice->GetMCApp()->Particle(l);
371 pdgCodes[l] = part->GetPdgCode();
375 printf("+++ Number of particles: %d\n",nParticles);
377 // Create the seeds for the TPC tracks at the inner radius of TPC
379 // Get TreeH with hits
380 TTree *th = tpcloader->TreeH();
381 MakeSeedsFromHits(tpc,th,sArray);
383 // Get TreeTR with track references
385 TTree *ttr = rl->TreeTR();
387 MakeSeedsFromRefs(ttr,sArray);
390 nSeeds = sArray.GetEntries();
391 printf("+++ Number of seeds: %d\n",nSeeds);
393 // loop over entries in sArray
394 for(Int_t l=0; l<nSeeds; l++) {
395 //if(l%1==0) printf(" --- Processing seed %d of %d ---\n",l,nSeeds);
397 seed = (AliTPCseedGeant*)sArray.At(l);
400 label = seed->GetLabel();
402 // check if this track has already been processed
403 if(done[label]) continue;
405 // PDG code & electric charge
406 pdg = pdgCodes[label];
407 if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
408 else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
410 pdg = TMath::Abs(pdg);
411 if(pdg>3000) pdg=211;
413 if(fSelAndSmear) SetParticle(pdg);
416 sEta = seed->GetEta();
418 // Apply selection according to TPC efficiency
419 //if(TMath::Abs(pdg)==211) nAcc++;
420 if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue;
421 //if(TMath::Abs(pdg)==211) nSel++;
423 // create AliTPCtrack object
424 BuildTrack(seed,charge);
427 bin = fDBgrid->GetBin(sPt,sEta);
430 //fCovTree = &(fCovTreePi[bin]);
431 fCovTree = fCovTreePi[bin];
434 //fCovTree = &(fCovTreeKa[bin]);
435 fCovTree = fCovTreeKa[bin];
438 //fCovTree = &(fCovTreePr[bin]);
439 fCovTree = fCovTreePr[bin];
442 //fCovTree = &(fCovTreeEl[bin]);
443 fCovTree = fCovTreeEl[bin];
446 //fCovTree = &(fCovTreeMu[bin]);
447 fCovTree = fCovTreeMu[bin];
450 // deal with covariance matrix and smearing of parameters
453 // assign the track a dE/dx and make a rough PID
457 // put track in array
458 AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
459 iotrack->SetLabel(label);
460 tArray.AddLast(iotrack);
461 // Mark track as "done" and register the pdg code
465 } // loop over entries in sArray
467 // sort array with TPC tracks (decreasing pT)
470 // convert to AliESDtrack and write to AliESD
471 arrentr = tArray.GetEntriesFast();
474 for(Int_t l=0; l<arrentr; l++) {
475 tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
476 AliESDtrack ioESDtrack;
477 ioESDtrack.UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
479 wgts[0]=0.; wgts[1]=0.; wgts[2]=0.; wgts[3]=0.; wgts[4]=0.;
480 if(TMath::Abs(tpctrack->GetMass()-0.9)<0.1) {
482 } else if(TMath::Abs(tpctrack->GetMass()-0.5)<0.1) {
488 ioESDtrack.SetESDpid(wgts);
489 event->AddTrack(&ioESDtrack);
496 printf("+++ Number of TPC tracks: %d\n",tracks);
497 //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
504 //-----------------------------------------------------------------------------
505 void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
506 //-----------------------------------------------------------------------------
507 // This function computes the dE/dx for pions, kaons, protons and electrons,
508 // in the [pT,eta] bins.
509 // Input file is CovMatrix_AllEvts.root.
510 //-----------------------------------------------------------------------------
512 gStyle->SetOptStat(0);
513 gStyle->SetOptFit(10001);
515 const char *part="PIONS";
519 // create a chain with compared tracks
520 TChain *cmptrkchain = new ("cmptrktree");
521 cmptrkchain.Add("CovMatrix_AllEvts.root");
522 //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
523 //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
524 //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
528 TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
529 TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
532 cmptrktree->SetBranchAddress("comptracks",&cmptrk);
533 Int_t entries = (Int_t)cmptrktree->GetEntries();
534 cerr<<" Number of entries: "<<entries<<endl;
536 InitializeKineGrid("DB");
537 InitializeKineGrid("dEdx");
564 const Int_t knTotBins = fDBgrid->GetTotBins();
566 cerr<<" Fit bins: "<<knTotBins<<endl;
569 Int_t *n = new Int_t[knTotBins];
570 Double_t *p = new Double_t[knTotBins];
571 Double_t *ep = new Double_t[knTotBins];
572 Double_t *mean = new Double_t[knTotBins];
573 Double_t *sigma = new Double_t[knTotBins];
575 for(Int_t l=0; l<knTotBins; l++) {
576 n[l] = 1; // set to 1 to avoid divisions by 0
577 p[l] = mean[l] = sigma[l] = ep[l] = 0.;
580 // loop on chain entries for the mean
581 for(Int_t l=0; l<entries; l++) {
582 cmptrktree->GetEvent(l);
583 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
584 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
586 mean[bin] += cmptrk.dEdx;
588 } // loop on chain entries
590 for(Int_t l=0; l<knTotBins; l++) {
593 n[l] = 1; // set to 1 to avoid divisions by 0
596 // loop on chain entries for the sigma
597 for(Int_t l=0; l<entries; l++) {
598 cmptrktree->GetEvent(l);
599 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
600 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
601 if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
603 sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
604 } // loop on chain entries
606 for(Int_t l=0; l<knTotBins; l++) {
607 sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
611 TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700);
613 // create the graph for dEdx vs p
614 TGraphErrors *gr = new TGraphErrors(knTotBins,p,mean,ep,sigma);
615 TString title(" : dE/dx vs momentum"); title.Prepend(part);
616 TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
617 frame->SetXTitle("p [GeV/c]");
618 frame->SetYTitle("dE/dx [a.u.]");
625 for(Int_t i=0; i<knTotBins; i++) {
626 fdEdxMeanPi.SetParam(i,mean[i]);
627 fdEdxRMSPi.SetParam(i,sigma[i]);
631 for(Int_t i=0; i<knTotBins; i++) {
632 fdEdxMeanKa.SetParam(i,mean[i]);
633 fdEdxRMSKa.SetParam(i,sigma[i]);
637 for(Int_t i=0; i<knTotBins; i++) {
638 fdEdxMeanPr.SetParam(i,mean[i]);
639 fdEdxRMSPr.SetParam(i,sigma[i]);
643 for(Int_t i=0; i<knTotBins; i++) {
644 fdEdxMeanEl.SetParam(i,mean[i]);
645 fdEdxRMSEl.SetParam(i,sigma[i]);
649 for(Int_t i=0; i<knTotBins; i++) {
650 fdEdxMeanMu.SetParam(i,mean[i]);
651 fdEdxRMSMu.SetParam(i,sigma[i]);
656 // write results to file
657 WritedEdx(outName,pdg);
667 //-----------------------------------------------------------------------------
668 void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
669 //-----------------------------------------------------------------------------
670 // This function computes the pulls for pions, kaons and electrons,
671 // in the [pT,eta] bins.
672 // Input file is CovMatrix_AllEvts.root.
673 // Output file is pulls.root.
674 //-----------------------------------------------------------------------------
677 // create a chain with compared tracks
678 TChain *cmptrkchain = new ("cmptrktree");
679 cmptrkchain.Add("CovMatrix_AllEvts.root");
680 //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
681 //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
682 //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
686 TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
687 TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
690 cmptrktree->SetBranchAddress("comptracks",&cmptrk);
691 Int_t entries = (Int_t)cmptrktree->GetEntries();
692 cerr<<" Number of entries: "<<entries<<endl;
695 Char_t hname[100], htitle[100];
698 AliTPCkineGrid pulls[5];
699 TH1F *hDum = new TH1F("name","title",100,-7.,7.);
700 TF1 *g = new TF1("g","gaus");
702 InitializeKineGrid("pulls");
703 InitializeKineGrid("DB");
707 // loop on the particles Pi,Ka,Pr,El,Mu
708 for(Int_t part=0; part<5; part++) {
713 cerr<<" Processing pions ...\n";
717 cerr<<" Processing kaons ...\n";
721 cerr<<" Processing protons ...\n";
725 cerr<<" Processing electrons ...\n";
729 cerr<<" Processing muons ...\n";
733 SetParticle(thisPdg);
735 for(Int_t i=0;i<5;i++) {
736 pulls[i].~AliTPCkineGrid();
737 new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
739 nTotBins = fDBgrid->GetTotBins();
740 cerr<<"nTotBins = "<<nTotBins<<endl;
742 // create histograms for the all the bins
749 hPulls0 = new TH1F[nTotBins];
750 hPulls1 = new TH1F[nTotBins];
751 hPulls2 = new TH1F[nTotBins];
752 hPulls3 = new TH1F[nTotBins];
753 hPulls4 = new TH1F[nTotBins];
756 for(Int_t i=0; i<nTotBins; i++) {
757 sprintf(hname,"hPulls0%d",i);
758 sprintf(htitle,"P0 pulls for bin %d",i);
759 hDum->SetName(hname); hDum->SetTitle(htitle);
761 sprintf(hname,"hPulls1%d",i);
762 sprintf(htitle,"P1 pulls for bin %d",i);
763 hDum->SetName(hname); hDum->SetTitle(htitle);
765 sprintf(hname,"hPulls2%d",i);
766 sprintf(htitle,"P2 pulls for bin %d",i);
767 hDum->SetName(hname); hDum->SetTitle(htitle);
769 sprintf(hname,"hPulls3%d",i);
770 sprintf(htitle,"P3 pulls for bin %d",i);
771 hDum->SetName(hname); hDum->SetTitle(htitle);
773 sprintf(hname,"hPulls4%d",i);
774 sprintf(htitle,"P4 pulls for bin %d",i);
775 hDum->SetName(hname); hDum->SetTitle(htitle);
779 // loop on chain entries
780 for(Int_t i=0; i<entries; i++) {
781 cmptrktree->GetEvent(i);
782 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
783 // fill histograms with the pulls
784 bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
785 //cerr<<" pt "<<cmptrk.pt<<" eta "<<cmptrk.eta<<" bin "<<bin<<endl;
786 hPulls0[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
787 hPulls1[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
788 hPulls2[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
789 hPulls3[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
790 hPulls4[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
791 } // loop on chain entries
793 // compute the sigma of the distributions
794 for(Int_t i=0; i<nTotBins; i++) {
795 if(hPulls0[i].GetEntries()>10) {
796 g->SetRange(-3.*hPulls0[i].GetRMS(),3.*hPulls0[i].GetRMS());
797 hPulls0[i].Fit("g","R,Q,N");
798 pulls[0].SetParam(i,g->GetParameter(2));
799 } else pulls[0].SetParam(i,-1.);
800 if(hPulls1[i].GetEntries()>10) {
801 g->SetRange(-3.*hPulls1[i].GetRMS(),3.*hPulls1[i].GetRMS());
802 hPulls1[i].Fit("g","R,Q,N");
803 pulls[1].SetParam(i,g->GetParameter(2));
804 } else pulls[1].SetParam(i,-1.);
805 if(hPulls2[i].GetEntries()>10) {
806 g->SetRange(-3.*hPulls2[i].GetRMS(),3.*hPulls2[i].GetRMS());
807 hPulls2[i].Fit("g","R,Q,N");
808 pulls[2].SetParam(i,g->GetParameter(2));
809 } else pulls[2].SetParam(i,-1.);
810 if(hPulls3[i].GetEntries()>10) {
811 g->SetRange(-3.*hPulls3[i].GetRMS(),3.*hPulls3[i].GetRMS());
812 hPulls3[i].Fit("g","R,Q,N");
813 pulls[3].SetParam(i,g->GetParameter(2));
814 } else pulls[3].SetParam(i,-1.);
815 if(hPulls4[i].GetEntries()>10) {
816 g->SetRange(-3.*hPulls4[i].GetRMS(),3.*hPulls4[i].GetRMS());
817 hPulls4[i].Fit("g","R,Q,N");
818 pulls[4].SetParam(i,g->GetParameter(2));
819 } else pulls[4].SetParam(i,-1.);
825 for(Int_t i=0;i<5;i++) {
826 fPullsPi[i].~AliTPCkineGrid();
827 new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
831 for(Int_t i=0;i<5;i++) {
832 fPullsKa[i].~AliTPCkineGrid();
833 new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
837 for(Int_t i=0;i<5;i++) {
838 fPullsPr[i].~AliTPCkineGrid();
839 new(&fPullsPr[i]) AliTPCkineGrid(pulls[i]);
843 for(Int_t i=0;i<5;i++) {
844 fPullsEl[i].~AliTPCkineGrid();
845 new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
849 for(Int_t i=0;i<5;i++) {
850 fPullsMu[i].~AliTPCkineGrid();
851 new(&fPullsMu[i]) AliTPCkineGrid(pulls[i]);
852 //cerr<<" mu pulls "<<i<<" "<<fPullsMu[i].GetParam(0)<<endl;
863 } // loop on particle species
865 // write pulls to file
871 //-----------------------------------------------------------------------------
872 void AliTPCtrackerParam::AnalyzeResolutions(Int_t pdg) {
873 //-----------------------------------------------------------------------------
874 // This function computes the resolutions:
878 // as a function of Pt
879 // Input file is CovMatrix_AllEvts.root.
880 //-----------------------------------------------------------------------------
883 // create a chain with compared tracks
884 TChain *cmptrkchain = new ("cmptrktree");
885 cmptrkchain.Add("CovMatrix_AllEvts.root");
886 //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
887 //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
888 //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
892 TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
893 TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
896 cmptrktree->SetBranchAddress("comptracks",&cmptrk);
897 Int_t entries = (Int_t)cmptrktree->GetEntries();
898 cerr<<" Number of entries: "<<entries<<endl;
903 InitializeKineGrid("DB");
904 InitializeKineGrid("eff");
908 const Int_t knPtBins = fEff->GetPointsPt();
909 cerr<<"knPtBins = "<<knPtBins<<endl;
910 Double_t *dP0 = new Double_t[knPtBins];
911 Double_t *dP4 = new Double_t[knPtBins];
912 Double_t *dPtToPt = new Double_t[knPtBins];
913 Double_t *pt = new Double_t[knPtBins];
914 fEff->GetArrayPt(pt);
917 TH1F *hDumP0 = new TH1F("nameP0","dy",100,-0.3,0.3);
918 TH1F *hDumP4 = new TH1F("nameP4","dC",100,-0.0005,0.0005);
919 TH1F *hDumPt = new TH1F("namePt","dp_{T}/p_{T}",100,-0.5,0.5);
921 TF1 *g = new TF1("g","gaus");
923 // create histograms for the all the bins
928 hP0 = new TH1F[knPtBins];
929 hP4 = new TH1F[knPtBins];
930 hPt = new TH1F[knPtBins];
932 for(Int_t i=0; i<knPtBins; i++) {
938 // loop on chain entries
939 for(Int_t i=0; i<entries; i++) {
940 cmptrktree->GetEvent(i);
941 if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
942 // fill histograms with the residuals
943 bin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
944 //cerr<<" pt "<<cmptrk.pt<<" eta "<<cmptrk.eta<<" bin "<<bin<<endl;
945 hP0[bin].Fill(cmptrk.dP0);
946 hP4[bin].Fill(cmptrk.dP4);
947 hPt[bin].Fill(cmptrk.dpt/cmptrk.pt);
948 } // loop on chain entries
951 TCanvas *cP0res = new TCanvas("cP0res","cP0res",0,0,1200,700);
953 TCanvas *cP4res = new TCanvas("cP4res","cP4res",0,0,1200,700);
955 TCanvas *cPtres = new TCanvas("cPtres","cPtres",0,0,1200,700);
959 for(Int_t i=0; i<knPtBins; i++) {
960 cP0res->cd(i+1); hP0[i].Draw();
961 cP4res->cd(i+1); hP4[i].Draw();
962 cPtres->cd(i+1); hPt[i].Draw();
966 // compute the sigma of the distributions
967 for(Int_t i=0; i<knPtBins; i++) {
968 if(hP0[i].GetEntries()>10) {
969 g->SetRange(-3.*hP0[i].GetRMS(),3.*hP0[i].GetRMS());
970 hP0[i].Fit("g","R,Q,N");
971 dP0[i] = g->GetParameter(2);
973 if(hP4[i].GetEntries()>10) {
974 g->SetRange(-3.*hP4[i].GetRMS(),3.*hP4[i].GetRMS());
975 hP4[i].Fit("g","R,Q,N");
976 dP4[i] = g->GetParameter(2);
978 if(hPt[i].GetEntries()>10) {
979 g->SetRange(-3.*hPt[i].GetRMS(),3.*hPt[i].GetRMS());
980 hPt[i].Fit("g","R,Q,N");
981 dPtToPt[i] = 100.*g->GetParameter(2);
982 } else dPtToPt[i] = 0.;
986 TGraph *grdP0 = new TGraph(knPtBins,pt,dP0);
987 TGraph *grdP4 = new TGraph(knPtBins,pt,dP4);
988 TGraph *grdPtToPt = new TGraph(knPtBins,pt,dPtToPt);
990 grdP0->SetMarkerStyle(20); grdP0->SetMarkerColor(2); grdP0->SetMarkerSize(1.5);
991 grdP4->SetMarkerStyle(21); grdP4->SetMarkerColor(3); grdP4->SetMarkerSize(1.5);
992 grdPtToPt->SetMarkerStyle(22); grdPtToPt->SetMarkerColor(4); grdPtToPt->SetMarkerSize(1.5);
995 gStyle->SetOptStat(0);
996 TCanvas *c1 = new TCanvas("c1","dP0",0,0,900,700);
1001 TH2F *frame1 = new TH2F("frame1","y resolution VS p_{T} in TPC",2,0.1,30,2,0,0.1);
1002 frame1->SetXTitle("p_{T} [GeV/c]");
1003 frame1->SetYTitle("#sigma(y) [cm]");
1008 TCanvas *c2 = new TCanvas("c2","dP4",0,0,900,700);
1013 TH2F *frame2 = new TH2F("frame2","C resolution VS p_{T} in TPC",2,0.1,30,2,0,0.0001);
1014 frame2->SetXTitle("p_{T} [GeV/c]");
1015 frame2->SetYTitle("#sigma(C) [1/cm]");
1019 TCanvas *c3 = new TCanvas("c3","dPtToPt",0,0,900,700);
1025 TH2F *frame3 = new TH2F("frame3","Relative p_{T} resolution VS p_{T} in TPC",2,0.1,30,2,0.1,30.);
1026 frame3->SetXTitle("p_{T} [GeV/c]");
1027 frame3->SetYTitle("dp_{T}/p_{T} (%)");
1029 grdPtToPt->Draw("P");
1044 //-----------------------------------------------------------------------------
1045 void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {
1046 //-----------------------------------------------------------------------------
1047 // This function uses GEANT info to set true track parameters
1048 //-----------------------------------------------------------------------------
1049 Double_t xref = s->GetXL();
1050 Double_t xx[5],cc[15];
1051 cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=0.;
1052 cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
1055 TVector3 bfield(0.,0.,-fBz);
1058 // radius [cm] of track projection in (x,y)
1059 Double_t rho = s->GetPt()*100./0.299792458/TMath::Abs(bfield.Z());
1060 // center of track projection in local reference frame
1064 // position (local) and momentum (local) at the seed
1065 // in the bending plane (z=0)
1066 sPos.SetXYZ(s->GetXL(),s->GetYL(),0.);
1067 sMom.SetXYZ(s->GetPx()*TMath::Cos(s->GetAlpha())+s->GetPy()*TMath::Sin(s->GetAlpha()),-s->GetPx()*TMath::Sin(s->GetAlpha())+s->GetPy()*TMath::Cos(s->GetAlpha()),0.);
1068 TVector3 vrho = sMom.Cross(bfield);
1072 TVector3 vcenter = sPos+vrho;
1074 Double_t x0 = vcenter.X();
1076 // fX = xref X-coordinate of this track (reference plane)
1077 // fAlpha = Alpha Rotation angle the local (TPC sector)
1078 // fP0 = YL Y-coordinate of a track
1079 // fP1 = ZG Z-coordinate of a track
1080 // fP2 = sin(phi) sine of the (local) azimuthal angle
1081 // fP3 = Tgl tangent of the track momentum dip angle
1082 // fP4 = C track curvature
1085 xx[2] = ch/rho*(xref-x0);
1086 xx[3] = s->GetPz()/s->GetPt();
1089 // create the object AliTPCtrack
1090 AliTPCtrack track(xref,s->GetAlpha(),xx,cc,0);
1091 new(&fTrack) AliTPCtrack(track);
1095 //-----------------------------------------------------------------------------
1096 void AliTPCtrackerParam::CompareTPCtracks(
1098 const Char_t* galiceName,
1099 const Char_t* trkGeaName,
1100 const Char_t* trkKalName,
1101 const Char_t* covmatName,
1102 const Char_t* tpceffasciiName,
1103 const Char_t* tpceffrootName) {
1104 //-----------------------------------------------------------------------------
1105 // This function compares tracks from TPC Kalman Filter V2 with
1106 // geant tracks at TPC 1st hit. It gives:
1107 // - a tree with Kalman cov. matrix elements, resolutions, dEdx
1108 // - the efficiencies as a function of particle type, pT, eta
1109 //-----------------------------------------------------------------------------
1111 TFile *kalFile = TFile::Open(trkKalName);
1112 TFile *geaFile = TFile::Open(trkGeaName);
1113 TFile *galiceFile = TFile::Open(galiceName);
1115 // get the AliRun object
1116 AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
1119 // create the tree for comparison results
1121 TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
1122 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");
1124 InitializeKineGrid("eff");
1125 InitializeKineGrid("DB");
1126 Int_t effBins = fEffPi.GetTotPoints();
1127 Int_t effBinsPt = fEffPi.GetPointsPt();
1128 Double_t *pt = new Double_t[effBinsPt];
1129 fEffPi.GetArrayPt(pt);
1136 Int_t pi=0,ka=0,mu=0,el=0,pr=0;
1137 Int_t *geaPi = new Int_t[effBins];
1138 Int_t *geaKa = new Int_t[effBins];
1139 Int_t *geaPr = new Int_t[effBins];
1140 Int_t *geaEl = new Int_t[effBins];
1141 Int_t *geaMu = new Int_t[effBins];
1142 Int_t *kalPi = new Int_t[effBins];
1143 Int_t *kalKa = new Int_t[effBins];
1144 Int_t *kalPr = new Int_t[effBins];
1145 Int_t *kalEl = new Int_t[effBins];
1146 Int_t *kalMu = new Int_t[effBins];
1147 Float_t *effPi = new Float_t[effBins];
1148 Float_t *effKa = new Float_t[effBins];
1149 Float_t *effPr = new Float_t[effBins];
1150 Float_t *effEl = new Float_t[effBins];
1151 Float_t *effMu = new Float_t[effBins];
1153 for(Int_t j=0; j<effBins; j++) {
1154 geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
1155 kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
1156 effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
1161 // loop on events in file
1162 for(Int_t evt=0; evt<nEvents; evt++) {
1163 cerr<<"\n --- Reading tracks for event "<<evt<<" ---\n\n";
1164 sprintf(tname,"TreeT_TPC_%d",evt);
1166 // particles from TreeK
1167 const Int_t knparticles = gAlice->GetEvent(evt);
1169 Int_t *kalLab = new Int_t[knparticles];
1170 for(Int_t i=0; i<knparticles; i++) kalLab[i] = -1;
1173 // tracks from Kalman
1174 TTree *kaltree=(TTree*)kalFile->Get(tname);
1175 if(!kaltree) continue;
1176 AliTPCtrack *kaltrack=new AliTPCtrack;
1177 kaltree->SetBranchAddress("tracks",&kaltrack);
1178 Int_t kalEntries = (Int_t)kaltree->GetEntries();
1180 // tracks from 1st hit
1181 TTree *geatree=(TTree*)geaFile->Get(tname);
1182 if(!geatree) continue;
1183 AliTPCtrack *geatrack=new AliTPCtrack;
1184 geatree->SetBranchAddress("tracks",&geatrack);
1185 Int_t geaEntries = (Int_t)geatree->GetEntries();
1187 cerr<<"+++\n+++ Number of tracks: TPC Kalman = "<<kalEntries<<endl<<"+++ TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
1189 // set pointers for TPC tracks:
1190 // the entry number of the track labelled l is stored in kalLab[l]
1191 Int_t fake=0,mult=0;
1192 for (Int_t j=0; j<kalEntries; j++) {
1193 kaltree->GetEvent(j);
1194 if(kaltrack->GetLabel()>=0) {
1195 if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
1196 kalLab[kaltrack->GetLabel()] = j;
1201 cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
1202 cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
1205 // Read the labels of the seeds
1208 Bool_t *hasSeed = new Bool_t[knparticles];
1209 for(Int_t i=0; i<knparticles; i++) hasSeed[i] = kFALSE;
1210 sprintf(sname,"seedLabels.%d.dat",evt);
1211 FILE *seedFile = fopen(sname,"r");
1213 ncol = fscanf(seedFile,"%d",&sLabel);
1215 if(sLabel>0) hasSeed[sLabel]=kTRUE;
1220 cerr<<"Doing track comparison...\n";
1221 // loop on tracks at 1st hit
1222 for(Int_t j=0; j<geaEntries; j++) {
1223 geatree->GetEvent(j);
1225 label = geatrack->GetLabel();
1226 part = (TParticle*)gAlice->GetMCApp()->Particle(label);
1228 // use only injected tracks with fixed values of pT
1229 ptgener = part->Pt();
1231 for(Int_t l=0; l<fEffPi.GetPointsPt(); l++) {
1232 if(TMath::Abs(ptgener-pt[l])<0.01) usethis = kTRUE;
1234 if(!usethis) continue;
1236 // check if it has the seed
1237 //if(!hasSeed[label]) continue;
1240 // check if track is entirely contained in a TPC sector
1241 Bool_t out = kFALSE;
1242 for(Int_t l=0; l<10; l++) {
1243 Double_t x = 85. + (250.-85.)*(Double_t)l/9.;
1244 //cerr<<"x "<<x<<" X "<<geatrack->GetX()<<endl;
1245 Double_t y = geatrack->GetY() + (
1246 TMath::Sqrt(1-(geatrack->GetC()*geatrack->GetX()-geatrack->GetEta())*
1247 (geatrack->GetC()*geatrack->GetX()-geatrack->GetEta()))
1248 -TMath::Sqrt(1-(geatrack->GetC()*x-geatrack->GetEta())*
1249 (geatrack->GetC()*x-geatrack->GetEta()))
1252 //cerr<<"y "<<y<<" Y "<<geatrack->GetY()<<endl;
1253 if(y > 0.8*x*TMath::Tan(TMath::Pi()/18)) { out = kTRUE; break; }
1258 cmptrk.pdg = part->GetPdgCode();
1259 cmptrk.eta = part->Eta();
1260 cmptrk.r = TMath::Sqrt(part->Vx()*part->Vx()+part->Vy()*part->Vy());
1262 cmptrk.pt = geatrack->Pt();
1263 cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
1264 cmptrk.p = cmptrk.pt/cmptrk.cosl;
1267 bin = fDBgridPi.GetBin(cmptrk.pt,cmptrk.eta);
1269 if(abs(cmptrk.pdg)==211) geaPi[bin]++;
1270 if(abs(cmptrk.pdg)==321) geaKa[bin]++;
1271 if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
1272 if(abs(cmptrk.pdg)==11) geaEl[bin]++;
1273 if(abs(cmptrk.pdg)==13) geaMu[bin]++;
1276 // check if there is the corresponding track in the TPC kalman and get it
1277 if(kalLab[label]==-1) continue;
1278 kaltree->GetEvent(kalLab[label]);
1280 // and go on only if it has xref = 84.57 cm (inner pad row)
1281 if(kaltrack->GetX()>90.) continue;
1283 if(abs(cmptrk.pdg)==211) { kalPi[bin]++; pi++; }
1284 if(abs(cmptrk.pdg)==321) { kalKa[bin]++; ka++; }
1285 if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
1286 if(abs(cmptrk.pdg)==11) { kalEl[bin]++; el++; }
1287 if(abs(cmptrk.pdg)==13) { kalMu[bin]++; mu++; }
1289 kaltrack->PropagateTo(geatrack->GetX(),1.e-9,0.);
1291 cmptrk.dEdx = kaltrack->GetdEdx();
1293 // compute errors on parameters
1294 dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
1295 if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
1297 cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
1298 cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
1299 cmptrk.dP2 = kaltrack->GetSnp()-geatrack->GetSnp();
1300 cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
1301 cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
1302 cmptrk.dpt = 1/kaltrack->GetSigned1Pt()-1/geatrack->GetSigned1Pt();
1304 // get covariance matrix
1305 // beware: lines 3 and 4 in the matrix are inverted!
1306 //kaltrack->GetCovariance(cc);
1308 cmptrk.c00 = kaltrack->GetSigmaY2();
1309 cmptrk.c10 = kaltrack->GetSigmaZY();
1310 cmptrk.c11 = kaltrack->GetSigmaZ2();
1311 cmptrk.c20 = kaltrack->GetSigmaSnpY();
1312 cmptrk.c21 = kaltrack->GetSigmaSnpY();
1313 cmptrk.c22 = kaltrack->GetSigmaSnp2();
1314 cmptrk.c30 = kaltrack->GetSigmaTglY();
1315 cmptrk.c31 = kaltrack->GetSigmaTglZ();
1316 cmptrk.c32 = kaltrack->GetSigmaTglSnp();
1317 cmptrk.c33 = kaltrack->GetSigmaTgl2();
1318 cmptrk.c40 = kaltrack->GetSigma1PtY();
1319 cmptrk.c41 = kaltrack->GetSigma1PtZ();
1320 cmptrk.c42 = kaltrack->GetSigma1PtSnp();
1321 cmptrk.c43 = kaltrack->GetSigma1PtTgl();
1322 cmptrk.c44 = kaltrack->GetSigma1Pt2();
1327 } // loop on tracks at TPC 1st hit
1330 //delete [] hasSeed;
1332 } // end loop on events in file
1335 cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
1336 cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
1339 // Write tree to file
1340 TFile *outfile = new TFile(covmatName,"recreate");
1341 cmptrktree->Write();
1345 // Write efficiencies to ascii file
1346 FILE *effFile = fopen(tpceffasciiName,"w");
1347 //fprintf(effFile,"%d\n",kalEntries);
1348 for(Int_t j=0; j<effBins; j++) {
1349 if(geaPi[j]>=100) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
1350 if(geaKa[j]>=100) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
1351 if(geaPr[j]>=100) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
1352 if(geaEl[j]>=500) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
1353 if(geaMu[j]>=100) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
1354 fprintf(effFile,"%f %f %f %f %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
1357 for(Int_t j=0; j<effBins; j++) {
1358 fprintf(effFile,"%d %d %d %d %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
1360 for(Int_t j=0; j<effBins; j++) {
1361 fprintf(effFile,"%d %d %d %d %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
1365 // Write efficiencies to root file
1366 for(Int_t j=0; j<effBins; j++) {
1367 fEffPi.SetParam(j,(Double_t)effPi[j]);
1368 fEffKa.SetParam(j,(Double_t)effKa[j]);
1369 fEffPr.SetParam(j,(Double_t)effPr[j]);
1370 fEffEl.SetParam(j,(Double_t)effEl[j]);
1371 fEffMu.SetParam(j,(Double_t)effMu[j]);
1373 WriteEffs(tpceffrootName);
1375 // delete AliRun object
1376 delete gAlice; gAlice=0;
1378 // close all input files
1381 galiceFile->Close();
1402 //-----------------------------------------------------------------------------
1403 void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
1404 //-----------------------------------------------------------------------------
1405 // This function assigns the track a dE/dx and makes a rough PID
1406 //-----------------------------------------------------------------------------
1408 Double_t mean = fdEdxMean->GetValueAt(pt,eta);
1409 Double_t rms = fdEdxRMS->GetValueAt(pt,eta);
1411 Double_t dEdx = gRandom->Gaus(mean,rms);
1413 fTrack.SetdEdx(dEdx);
1415 AliTPCtrackParam t(fTrack);
1418 Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
1420 Double_t massPi = (Double_t)TDatabasePDG::Instance()->GetParticle(211)->Mass();
1421 Double_t massKa = (Double_t)TDatabasePDG::Instance()->GetParticle(321)->Mass();
1422 Double_t massPr = (Double_t)TDatabasePDG::Instance()->GetParticle(2212)->Mass();
1425 if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
1426 t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
1428 if (dEdx < 39.+ 12./p/p) {
1429 t.AssignMass(massKa); new(&fTrack) AliTPCtrack(t); return;
1431 t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return;
1435 if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
1436 t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
1438 t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return;
1441 t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
1443 //-----------------------------------------------------------------------------
1444 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
1445 //-----------------------------------------------------------------------------
1446 // This function deals with covariance matrix and smearing
1447 //-----------------------------------------------------------------------------
1451 Double_t trkKine[1],trkRegPar[3];
1452 Double_t xref,alpha,xx[5],xxsm[5],cc[15];
1455 fCovTree->SetBranchAddress("matrix",&covmat);
1457 // get random entry from the tree
1458 treeEntries = (Int_t)fCovTree->GetEntries();
1459 fCovTree->GetEvent(gRandom->Integer(treeEntries));
1461 // get P and Cosl from track
1462 cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1463 p = fTrack.Pt()/cosl;
1467 // get covariance matrix from regularized matrix
1468 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(0,l);
1469 cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1471 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(1,l);
1472 cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1473 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(2,l);
1474 cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1476 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(3,l);
1477 cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1479 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(4,l);
1480 cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1482 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(5,l);
1483 cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1484 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(6,l);
1485 cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1487 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(7,l);
1488 cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1490 for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l);
1491 cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1494 TMatrixD covMatSmear(5,5);
1495 covMatSmear = GetSmearingMatrix(cc,pt,eta);
1497 // get track original parameters
1498 xref =fTrack.GetX();
1499 alpha=fTrack.GetAlpha();
1500 xx[0]=fTrack.GetY();
1501 xx[1]=fTrack.GetZ();
1502 xx[2]=fTrack.GetSnp();
1503 xx[3]=fTrack.GetTgl();
1504 xx[4]=fTrack.GetC();
1506 // use smearing matrix to smear the original parameters
1508 SmearTrack(xx,xxsm,covMatSmear);
1510 AliTPCtrack track(xref,alpha,xxsm,cc,0);
1511 new(&fTrack) AliTPCtrack(track);
1515 //-----------------------------------------------------------------------------
1516 void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg) {
1517 //-----------------------------------------------------------------------------
1518 // This function draws the TPC efficiencies in the [pT,eta] bins
1519 //-----------------------------------------------------------------------------
1524 const Int_t kn = fEff->GetPointsPt();
1525 Double_t *effsA = new Double_t[kn];
1526 Double_t *effsB = new Double_t[kn];
1527 Double_t *effsC = new Double_t[kn];
1528 Double_t *pt = new Double_t[kn];
1530 fEff->GetArrayPt(pt);
1531 for(Int_t i=0;i<kn;i++) {
1532 effsA[i] = fEff->GetParam(i,0);
1533 effsB[i] = fEff->GetParam(i,1);
1534 effsC[i] = fEff->GetParam(i,2);
1537 TGraph *grA = new TGraph(kn,pt,effsA);
1538 TGraph *grB = new TGraph(kn,pt,effsB);
1539 TGraph *grC = new TGraph(kn,pt,effsC);
1541 grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1542 grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1543 grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1545 TString title("Distribution of the TPC efficiencies");
1548 title.Prepend("PIONS - ");
1551 title.Prepend("KAONS - ");
1554 title.Prepend("PROTONS - ");
1557 title.Prepend("ELECTRONS - ");
1560 title.Prepend("MUONS - ");
1565 gStyle->SetOptStat(0);
1566 TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1571 TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1572 frame->SetXTitle("p_{T} [GeV/c]");
1578 TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1579 leg->AddEntry(grA,"|#eta|<0.3","p");
1580 leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1581 leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1582 leg->SetFillColor(0);
1592 //-----------------------------------------------------------------------------
1593 void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg,
1595 //-----------------------------------------------------------------------------
1596 // This function draws the pulls in the [pT,eta] bins
1597 //-----------------------------------------------------------------------------
1602 const Int_t kn = (fPulls+par)->GetPointsPt();
1603 Double_t *pullsA = new Double_t[kn];
1604 Double_t *pullsB = new Double_t[kn];
1605 Double_t *pullsC = new Double_t[kn];
1606 Double_t *pt = new Double_t[kn];
1607 (fPulls+par)->GetArrayPt(pt);
1608 for(Int_t i=0;i<kn;i++) {
1609 pullsA[i] = (fPulls+par)->GetParam(i,0);
1610 pullsB[i] = (fPulls+par)->GetParam(i,1);
1611 pullsC[i] = (fPulls+par)->GetParam(i,2);
1614 TGraph *grA = new TGraph(kn,pt,pullsA);
1615 TGraph *grB = new TGraph(kn,pt,pullsB);
1616 TGraph *grC = new TGraph(kn,pt,pullsC);
1618 grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1619 grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1620 grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1622 TString title("Distribution of the pulls: ");
1625 title.Prepend("PIONS - ");
1628 title.Prepend("KAONS - ");
1631 title.Prepend("PROTONS - ");
1634 title.Prepend("ELECTRONS - ");
1637 title.Prepend("MUONS - ");
1648 title.Append(" #eta");
1651 title.Append("tg #lambda");
1658 gStyle->SetOptStat(0);
1659 TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1664 TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1665 frame->SetXTitle("p_{T} [GeV/c]");
1671 TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1672 leg->AddEntry(grA,"|#eta|<0.3","p");
1673 leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1674 leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1675 leg->SetFillColor(0);
1685 //-----------------------------------------------------------------------------
1686 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1687 Double_t eta) const {
1688 //-----------------------------------------------------------------------------
1689 // This function stretches the covariance matrix according to the pulls
1690 //-----------------------------------------------------------------------------
1692 TMatrixD covMat(5,5);
1695 covMat(1,0)=cc[1]; covMat(0,1)=covMat(1,0);
1697 covMat(2,0)=cc[3]; covMat(0,2)=covMat(2,0);
1698 covMat(2,1)=cc[4]; covMat(1,2)=covMat(2,1);
1700 covMat(3,0)=cc[6]; covMat(0,3)=covMat(3,0);
1701 covMat(3,1)=cc[7]; covMat(1,3)=covMat(3,1);
1702 covMat(3,2)=cc[8]; covMat(2,3)=covMat(3,2);
1704 covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1705 covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1706 covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1707 covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1711 TMatrixD stretchMat(5,5);
1712 for(Int_t k=0;k<5;k++) {
1713 for(Int_t l=0;l<5;l++) {
1718 for(Int_t i=0;i<5;i++) {
1719 stretchMat(i,i) = (fPulls+i)->GetValueAt(pt,eta);
1720 if(stretchMat(i,i)==0.) stretchMat(i,i) = 1.;
1723 TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1724 TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1728 //-----------------------------------------------------------------------------
1729 void AliTPCtrackerParam::InitializeKineGrid(Option_t* which) {
1730 //-----------------------------------------------------------------------------
1731 // This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1732 // which = "DB" -> initialize fDBgrid... members
1733 // "eff" -> initialize fEff... members
1734 // "pulls" -> initialize fPulls... members
1735 // "dEdx" -> initialize fdEdx... members
1736 //-----------------------------------------------------------------------------
1738 const char *db = strstr(which,"DB");
1739 const char *eff = strstr(which,"eff");
1740 const char *pulls = strstr(which,"pulls");
1741 const char *dEdx = strstr(which,"dEdx");
1744 Int_t nEta=0, nPt=0;
1746 Double_t etaPoints[2] = {0.3,0.6};
1747 Double_t etaBins[3] = {0.15,0.45,0.75};
1749 Double_t ptPoints[9] = {0.4,0.6,0.8,1.2,1.7,3.,5.,8.,15.};
1750 Double_t ptBins[10] = {0.3,0.5,0.7,1.,1.5,2.,4.,6.,10.,20.};
1753 Double_t *eta=0,*pt=0;
1767 AliTPCkineGrid *dummy=0;
1770 dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1771 new(&fDBgridPi) AliTPCkineGrid(*dummy);
1772 new(&fDBgridKa) AliTPCkineGrid(*dummy);
1773 new(&fDBgridPr) AliTPCkineGrid(*dummy);
1774 new(&fDBgridEl) AliTPCkineGrid(*dummy);
1775 new(&fDBgridMu) AliTPCkineGrid(*dummy);
1779 dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1780 new(&fEffPi) AliTPCkineGrid(*dummy);
1781 new(&fEffKa) AliTPCkineGrid(*dummy);
1782 new(&fEffPr) AliTPCkineGrid(*dummy);
1783 new(&fEffEl) AliTPCkineGrid(*dummy);
1784 new(&fEffMu) AliTPCkineGrid(*dummy);
1788 dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1789 for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1790 for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1791 for(Int_t i=0;i<5;i++) new(&fPullsPr[i]) AliTPCkineGrid(*dummy);
1792 for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1793 for(Int_t i=0;i<5;i++) new(&fPullsMu[i]) AliTPCkineGrid(*dummy);
1797 dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1798 new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1799 new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1800 new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1801 new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1802 new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1803 new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1804 new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1805 new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1806 new(&fdEdxMeanMu) AliTPCkineGrid(*dummy);
1807 new(&fdEdxRMSMu) AliTPCkineGrid(*dummy);
1813 //-----------------------------------------------------------------------------
1814 void AliTPCtrackerParam::MakeDataBase() {
1815 //-----------------------------------------------------------------------------
1816 // This function creates the DB file and store in it:
1817 // - TPC Efficiencies for pi,ka,pr,el,mu (AliTPCkineGrid class)
1818 // - Pulls for pi,ka,el (AliTPCkineGrid class)
1819 // - Regularization parameters for pi,ka,el (TMatrixD class)
1820 // - dE/dx parameterization for pi,ka,pr,el (AliTPCkineGrid class)
1821 // - Regularized cov. matrices for pi,ka,el (COVMATRIX structure)
1822 //-----------------------------------------------------------------------------
1824 // define some file names
1825 const char *effFile ="TPCeff.root";
1826 const char *pullsFile ="pulls.root";
1827 const char *regPiFile ="regPi.root";
1828 const char *regKaFile ="regKa.root";
1829 const char *regPrFile ="regPr.root";
1830 const char *regElFile ="regEl.root";
1831 const char *regMuFile ="regMu.root";
1832 const char *dEdxPiFile="dEdxPi.root";
1833 const char *dEdxKaFile="dEdxKa.root";
1834 const char *dEdxPrFile="dEdxPr.root";
1835 const char *dEdxElFile="dEdxEl.root";
1836 const char *dEdxMuFile="dEdxMu.root";
1837 const char *cmFile ="CovMatrix_AllEvts.root";
1839 Char_t *cmFile1 ="CovMatrix_AllEvts_1.root";
1840 Char_t *cmFile2 ="CovMatrix_AllEvts_2.root";
1841 Char_t *cmFile3 ="CovMatrix_AllEvts_3.root";
1844 // store the effieciencies
1846 WriteEffs(fDBfileName.Data());
1849 ReadPulls(pullsFile);
1850 WritePulls(fDBfileName.Data());
1852 //store the regularization parameters
1853 ReadRegParams(regPiFile,211);
1854 WriteRegParams(fDBfileName.Data(),211);
1855 ReadRegParams(regKaFile,321);
1856 WriteRegParams(fDBfileName.Data(),321);
1857 ReadRegParams(regPrFile,2212);
1858 WriteRegParams(fDBfileName.Data(),2212);
1859 ReadRegParams(regElFile,11);
1860 WriteRegParams(fDBfileName.Data(),11);
1861 ReadRegParams(regMuFile,13);
1862 WriteRegParams(fDBfileName.Data(),13);
1864 // store the dEdx parameters
1865 ReaddEdx(dEdxPiFile,211);
1866 WritedEdx(fDBfileName.Data(),211);
1867 ReaddEdx(dEdxKaFile,321);
1868 WritedEdx(fDBfileName.Data(),321);
1869 ReaddEdx(dEdxPrFile,2212);
1870 WritedEdx(fDBfileName.Data(),2212);
1871 ReaddEdx(dEdxElFile,11);
1872 WritedEdx(fDBfileName.Data(),11);
1873 ReaddEdx(dEdxMuFile,13);
1874 WritedEdx(fDBfileName.Data(),13);
1878 // store the regularized covariance matrices
1880 InitializeKineGrid("DB");
1882 const Int_t knBinsPi = fDBgridPi.GetTotBins();
1883 const Int_t knBinsKa = fDBgridKa.GetTotBins();
1884 const Int_t knBinsPr = fDBgridPr.GetTotBins();
1885 const Int_t knBinsEl = fDBgridEl.GetTotBins();
1886 const Int_t knBinsMu = fDBgridMu.GetTotBins();
1889 // create the trees for cov. matrices
1891 TTree *covTreePi1 = NULL;
1892 covTreePi1 = new TTree[knBinsPi];
1894 TTree *covTreeKa1 = NULL;
1895 covTreeKa1 = new TTree[knBinsKa];
1896 // trees for protons
1897 TTree *covTreePr1 = NULL;
1898 covTreePr1 = new TTree[knBinsPr];
1899 // trees for electrons
1900 TTree *covTreeEl1 = NULL;
1901 covTreeEl1 = new TTree[knBinsEl];
1903 TTree *covTreeMu1 = NULL;
1904 covTreeMu1 = new TTree[knBinsMu];
1906 Char_t hname[100], htitle[100];
1910 for(Int_t i=0; i<knBinsPi; i++) {
1911 sprintf(hname,"CovTreePi_bin%d",i);
1912 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1913 covTreePi1[i].SetName(hname); covTreePi1[i].SetTitle(htitle);
1914 covTreePi1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1916 for(Int_t i=0; i<knBinsKa; i++) {
1917 sprintf(hname,"CovTreeKa_bin%d",i);
1918 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1919 covTreeKa1[i].SetName(hname); covTreeKa1[i].SetTitle(htitle);
1920 covTreeKa1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1922 for(Int_t i=0; i<knBinsPr; i++) {
1923 sprintf(hname,"CovTreePr_bin%d",i);
1924 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1925 covTreePr1[i].SetName(hname); covTreePr1[i].SetTitle(htitle);
1926 covTreePr1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1928 for(Int_t i=0; i<knBinsEl; i++) {
1929 sprintf(hname,"CovTreeEl_bin%d",i);
1930 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1931 covTreeEl1[i].SetName(hname); covTreeEl1[i].SetTitle(htitle);
1932 covTreeEl1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1934 for(Int_t i=0; i<knBinsMu; i++) {
1935 sprintf(hname,"CovTreeMu_bin%d",i);
1936 sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1937 covTreeMu1[i].SetName(hname); covTreeMu1[i].SetTitle(htitle);
1938 covTreeMu1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1942 // create the chain with the compared tracks
1943 TChain *cmptrktree = new TChain("cmptrktree");
1944 cmptrkchain.Add(cmFile1);
1945 cmptrkchain.Add(cmFile2);
1946 cmptrkchain.Add(cmFile3);
1949 TFile *infile = TFile::Open(cmFile);
1950 TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
1953 cmptrktree->SetBranchAddress("comptracks",&cmptrk);
1954 Int_t entries = (Int_t)cmptrktree->GetEntries();
1955 cerr<<" Number of entries: "<<entries<<endl;
1957 Int_t trkPdg,trkBin;
1958 Double_t trkKine[1],trkRegPar[3];
1959 Int_t *nPerBinPi = new Int_t[knBinsPi];
1960 for(Int_t k=0;k<knBinsPi;k++) nPerBinPi[k]=0;
1961 Int_t *nPerBinKa = new Int_t[knBinsKa];
1962 for(Int_t k=0;k<knBinsKa;k++) nPerBinKa[k]=0;
1963 Int_t *nPerBinMu = new Int_t[knBinsMu];
1964 for(Int_t k=0;k<knBinsMu;k++) nPerBinMu[k]=0;
1965 Int_t *nPerBinEl = new Int_t[knBinsEl];
1966 for(Int_t k=0;k<knBinsEl;k++) nPerBinEl[k]=0;
1967 Int_t *nPerBinPr = new Int_t[knBinsPr];
1968 for(Int_t k=0;k<knBinsPr;k++) nPerBinPr[k]=0;
1970 // loop on chain entries
1971 for(Int_t l=0; l<entries; l++) {
1972 if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1974 cmptrktree->GetEvent(l);
1976 trkPdg = TMath::Abs(cmptrk.pdg);
1977 // use only pions, kaons, protons, electrons, muons
1978 if(trkPdg!=211 && trkPdg!=321 && trkPdg!=2212 && trkPdg!=11 && trkPdg!=13) continue;
1979 SetParticle(trkPdg);
1980 trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1981 //cerr<<cmptrk.pt<<" "<<cmptrk.eta<<" "<<trkBin<<endl;
1983 if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1984 if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1985 if(trkPdg==2212 && nPerBinPr[trkBin]>=5000) continue;
1986 if(trkPdg==11 && nPerBinEl[trkBin]>=5000) continue;
1987 if(trkPdg==13 && nPerBinMu[trkBin]>=5000) continue;
1989 trkKine[0] = cmptrk.p;
1991 // get regularized covariance matrix
1992 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(0,k);
1993 covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1994 covmat.c10 = cmptrk.c10;
1995 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(1,k);
1996 covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1997 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(2,k);
1998 covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1999 covmat.c21 = cmptrk.c21;
2000 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(3,k);
2001 covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
2002 covmat.c30 = cmptrk.c30;
2003 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(4,k);
2004 covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
2005 covmat.c32 = cmptrk.c32;
2006 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(5,k);
2007 covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
2008 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(6,k);
2009 covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
2010 covmat.c41 = cmptrk.c41;
2011 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(7,k);
2012 covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
2013 covmat.c43 = cmptrk.c43;
2014 for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(8,k);
2015 covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
2020 covTreePi1[trkBin].Fill();
2021 nPerBinPi[trkBin]++;
2024 covTreeKa1[trkBin].Fill();
2025 nPerBinKa[trkBin]++;
2027 case 2212: // protons
2028 covTreePr1[trkBin].Fill();
2029 nPerBinPr[trkBin]++;
2031 case 11: // electrons
2032 covTreeEl1[trkBin].Fill();
2033 nPerBinEl[trkBin]++;
2036 covTreeMu1[trkBin].Fill();
2037 nPerBinMu[trkBin]++;
2040 } // loop on chain entries
2042 // store all trees the DB file
2043 TFile *dbfile = new TFile(fDBfileName.Data(),"update");
2044 dbfile->mkdir("CovMatrices");
2045 gDirectory->cd("/CovMatrices");
2046 gDirectory->mkdir("Pions");
2047 gDirectory->mkdir("Kaons");
2048 gDirectory->mkdir("Protons");
2049 gDirectory->mkdir("Electrons");
2050 gDirectory->mkdir("Muons");
2052 gDirectory->cd("/CovMatrices/Pions");
2053 fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
2054 for(Int_t i=0;i<knBinsPi;i++) covTreePi1[i].Write();
2056 gDirectory->cd("/CovMatrices/Kaons");
2057 fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
2058 for(Int_t i=0;i<knBinsKa;i++) covTreeKa1[i].Write();
2060 gDirectory->cd("/CovMatrices/Protons");
2061 fDBgridPr.SetName("DBgridPr"); fDBgridPr.Write();
2062 for(Int_t i=0;i<knBinsPr;i++) covTreePr1[i].Write();
2064 gDirectory->cd("/CovMatrices/Electrons");
2065 fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
2066 for(Int_t i=0;i<knBinsEl;i++) covTreeEl1[i].Write();
2068 gDirectory->cd("/CovMatrices/Muons");
2069 fDBgridMu.SetName("DBgridMu"); fDBgridMu.Write();
2070 for(Int_t i=0;i<knBinsMu;i++) covTreeMu1[i].Write();
2073 delete [] nPerBinPi;
2074 delete [] nPerBinKa;
2075 delete [] nPerBinPr;
2076 delete [] nPerBinEl;
2077 delete [] nPerBinMu;
2081 //-----------------------------------------------------------------------------
2082 void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *tpc,TTree *th,
2083 TObjArray &seedArray) const {
2084 //-----------------------------------------------------------------------------
2085 // This function makes the seeds for tracks from the 1st hits in the TPC
2086 //-----------------------------------------------------------------------------
2088 Double_t xg,yg,zg,px,py,pz,pt;
2090 Int_t nTracks=(Int_t)th->GetEntries();
2092 cout<<"+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
2095 AliTPChit *tpcHit=0;
2097 // loop over entries in TreeH
2098 for(Int_t l=0; l<nTracks; l++) {
2099 if(l%1000==0) cerr<<" --- Processing primary track "
2100 <<l<<" of "<<nTracks<<" ---\r";
2104 tpcHit=(AliTPChit*)tpc->FirstHit(-1);
2105 for( ; tpcHit; tpcHit=(AliTPChit*)tpc->NextHit() ) {
2106 if(tpcHit->fQ !=0.) continue;
2107 // Get particle momentum at hit
2108 px=tpcHit->X(); py=tpcHit->Y(); pz=tpcHit->Z();
2110 pt=TMath::Sqrt(px*px+py*py);
2111 // reject hits with Pt<mag*0.45 GeV/c
2112 if(pt<(fBz*0.45)) continue;
2115 label=tpcHit->Track();
2117 if((tpcHit=(AliTPChit*)tpc->NextHit())==0) break;
2118 if(tpcHit->fQ != 0.) continue;
2119 // Get global coordinates of hit
2120 xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
2121 if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
2124 AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
2126 // reject tracks which are not in the TPC acceptance
2127 if(!ioseed->InTPCAcceptance()) { delete ioseed; continue; }
2129 // put seed in array
2130 seedArray.AddLast(ioseed);
2133 } // loop over entries in TreeH
2137 //-----------------------------------------------------------------------------
2138 void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr,
2139 TObjArray &seedArray) const {
2140 //-----------------------------------------------------------------------------
2141 // This function makes the seeds for tracks from the track references
2142 //-----------------------------------------------------------------------------
2144 Double_t xg,yg,zg,px,py,pz,pt;
2148 TClonesArray *tkRefArray = new TClonesArray("AliTrackReference");
2150 TBranch *b =(TBranch*)ttr->GetBranch("TPC");
2151 if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
2152 b->SetAddress(&tkRefArray);
2153 Int_t nTkRef = (Int_t)b->GetEntries();
2154 cerr<<"+++ Number of entries in TreeTR(TPC): "<<nTkRef<<"\n";
2156 // loop on track references
2157 for(Int_t l=0; l<nTkRef; l++){
2158 //if(l%1000==0) cerr<<" --- Processing primary track "
2159 // <<l<<" of "<<nTkRef<<" ---\r";
2160 if(!b->GetEvent(l)) continue;
2161 nnn = tkRefArray->GetEntriesFast();
2163 if(nnn <= 0) continue;
2164 for(k=0; k<nnn; k++) {
2165 AliTrackReference *tkref = (AliTrackReference*)tkRefArray->UncheckedAt(k);
2173 label = tkref->GetTrack();
2175 pt=TMath::Sqrt(px*px+py*py);
2176 // reject hits with Pt<mag*0.45 GeV/c
2177 if(pt<(fBz*0.45)) continue;
2180 AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
2182 // reject if not at the inner part of TPC
2183 if(TMath::Abs(ioseed->GetXL()-83.8) > 0.2) {
2184 //cerr<<"not at TPC inner part: XL = "<<ioseed->GetXL()<<endl;
2185 delete ioseed; continue;
2188 // reject tracks which are not in the TPC acceptance
2189 if(!ioseed->InTPCAcceptance()) {
2190 delete ioseed; continue;
2193 // put seed in array
2194 seedArray.AddLast(ioseed);
2199 } // loop on track references
2206 //-----------------------------------------------------------------------------
2207 void AliTPCtrackerParam::MergeEvents(Int_t evFirst,Int_t evLast) {
2208 //-----------------------------------------------------------------------------
2209 // This function: 1) merges the files from track comparison
2210 // (beware: better no more than 100 events per file)
2211 // 2) computes the average TPC efficiencies
2212 //-----------------------------------------------------------------------------
2214 const char *outName="TPCeff.root";
2216 // merge files with tracks
2217 cerr<<" ******** MERGING FILES **********\n\n";
2219 // create the chain for the tree of compared tracks
2220 TChain ch1("cmptrktree");
2221 TChain ch2("cmptrktree");
2222 TChain ch3("cmptrktree");
2224 for(Int_t j=evFirst; j<=evLast; j++) {
2225 cerr<<"Processing event "<<j<<endl;
2227 TString covName("CovMatrix.");
2229 covName.Append(".root");
2231 if(gSystem->AccessPathName(covName.Data(),kFileExists)) continue;
2234 if(j<=100) ch1.Add(covName.Data());
2235 if(j>100 && j<=200) ch2.Add(covName.Data());
2236 if(j>200) ch3.Add(covName.Data());
2240 // merge chain in one file
2242 covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
2243 ch1.Merge(covOut,1000000000);
2246 covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
2247 ch2.Merge(covOut,1000000000);
2250 covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
2251 ch3.Merge(covOut,1000000000);
2257 cerr<<" ***** EFFICIENCIES ******\n\n";
2259 ReadEffs("TPCeff.1.root");
2261 Int_t n = fEffPi.GetTotPoints();
2262 Double_t *avEffPi = new Double_t[n];
2263 Double_t *avEffKa = new Double_t[n];
2264 Double_t *avEffPr = new Double_t[n];
2265 Double_t *avEffEl = new Double_t[n];
2266 Double_t *avEffMu = new Double_t[n];
2267 Int_t *evtsPi = new Int_t[n];
2268 Int_t *evtsKa = new Int_t[n];
2269 Int_t *evtsPr = new Int_t[n];
2270 Int_t *evtsEl = new Int_t[n];
2271 Int_t *evtsMu = new Int_t[n];
2273 for(Int_t j=0; j<n; j++) {
2274 avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
2275 evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
2278 for(Int_t j=evFirst; j<=evLast; j++) {
2279 cerr<<"Processing event "<<j<<endl;
2281 TString effName("TPCeff.");
2283 effName.Append(".root");
2285 if(gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
2287 ReadEffs(effName.Data());
2289 for(Int_t k=0; k<n; k++) {
2290 if(fEffPi.GetParam(k)>=0.) {avEffPi[k]+=fEffPi.GetParam(k); evtsPi[k]++;}
2291 if(fEffKa.GetParam(k)>=0.) {avEffKa[k]+=fEffKa.GetParam(k); evtsKa[k]++;}
2292 if(fEffPr.GetParam(k)>=0.) {avEffPr[k]+=fEffPr.GetParam(k); evtsPr[k]++;}
2293 if(fEffEl.GetParam(k)>=0.) {avEffEl[k]+=fEffEl.GetParam(k); evtsEl[k]++;}
2294 if(fEffMu.GetParam(k)>=0.) {avEffMu[k]+=fEffMu.GetParam(k); evtsMu[k]++;}
2299 // compute average efficiencies
2300 for(Int_t j=0; j<n; j++) {
2301 if(evtsPi[j]==0) evtsPi[j]++;
2302 fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
2303 if(evtsKa[j]==0) evtsKa[j]++;
2304 fEffKa.SetParam(j,(Double_t)avEffKa[j]/evtsKa[j]);
2305 if(evtsPr[j]==0) evtsPr[j]++;
2306 fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
2307 if(evtsEl[j]==0) evtsEl[j]++;
2308 fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
2309 if(evtsMu[j]==0) evtsMu[j]++;
2310 fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
2313 // write efficiencies to a file
2329 //-----------------------------------------------------------------------------
2330 Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
2331 //-----------------------------------------------------------------------------
2332 // This function reads all parameters from the DB
2333 //-----------------------------------------------------------------------------
2335 if(!ReadEffs(inName)) return 0;
2336 if(!ReadPulls(inName)) return 0;
2337 if(!ReadRegParams(inName,211)) return 0;
2338 if(!ReadRegParams(inName,321)) return 0;
2339 if(!ReadRegParams(inName,2212)) return 0;
2340 if(!ReadRegParams(inName,11)) return 0;
2341 if(!ReadRegParams(inName,13)) return 0;
2342 if(!ReaddEdx(inName,211)) return 0;
2343 if(!ReaddEdx(inName,321)) return 0;
2344 if(!ReaddEdx(inName,2212)) return 0;
2345 if(!ReaddEdx(inName,11)) return 0;
2346 if(!ReaddEdx(inName,13)) return 0;
2347 if(!ReadDBgrid(inName)) return 0;
2348 if(!ReadCovTrees(inName)) return 0;
2352 //-----------------------------------------------------------------------------
2353 Int_t AliTPCtrackerParam::ReadCovTrees(const Char_t* inName) {
2354 //-----------------------------------------------------------------------------
2355 // This function reads the covariance matrix trees from the DB
2356 //-----------------------------------------------------------------------------
2358 if(gSystem->AccessPathName(inName,kFileExists)) {
2359 cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
2364 TFile *inFile = TFile::Open(inName);
2367 Int_t nBinsPi = fDBgridPi.GetTotBins();
2368 for(Int_t l=0; l<nBinsPi; l++) {
2369 str = "/CovMatrices/Pions/CovTreePi_bin";
2371 fCovTreePi[l] = (TTree*)inFile->Get(str.Data());
2372 //fCovTree = (TTree*)inFile->Get(str.Data());
2373 //fCovTreePi[l] = new TTree(*fCovTree);
2375 Int_t nBinsKa = fDBgridKa.GetTotBins();
2376 for(Int_t l=0; l<nBinsKa; l++) {
2377 str = "/CovMatrices/Kaons/CovTreeKa_bin";
2379 fCovTreeKa[l] = (TTree*)inFile->Get(str.Data());
2380 //fCovTree = (TTree*)inFile->Get(str.Data());
2381 //fCovTreeKa[l] = new TTree(*fCovTree);
2383 Int_t nBinsPr = fDBgridPr.GetTotBins();
2384 for(Int_t l=0; l<nBinsPr; l++) {
2386 str = "/CovMatrices/Pions/CovTreePi_bin";
2388 str = "/CovMatrices/Protons/CovTreePr_bin";
2391 fCovTreePr[l] = (TTree*)inFile->Get(str.Data());
2392 //fCovTree = (TTree*)inFile->Get(str.Data());
2393 //fCovTreePr[l] = new TTree(*fCovTree);
2395 Int_t nBinsEl = fDBgridEl.GetTotBins();
2396 for(Int_t l=0; l<nBinsEl; l++) {
2397 str = "/CovMatrices/Electrons/CovTreeEl_bin";
2399 fCovTreeEl[l] = (TTree*)inFile->Get(str.Data());
2400 //fCovTree = (TTree*)inFile->Get(str.Data());
2401 //fCovTreeEl[l] = new TTree(*fCovTree);
2403 Int_t nBinsMu = fDBgridMu.GetTotBins();
2404 for(Int_t l=0; l<nBinsMu; l++) {
2406 str = "/CovMatrices/Pions/CovTreePi_bin";
2408 str = "/CovMatrices/Muons/CovTreeMu_bin";
2411 fCovTreeMu[l] = (TTree*)inFile->Get(str.Data());
2412 //fCovTree = (TTree*)inFile->Get(str.Data());
2413 //fCovTreeMu[l] = new TTree(*fCovTree);
2420 //-----------------------------------------------------------------------------
2421 Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
2422 //-----------------------------------------------------------------------------
2423 // This function reads the dEdx parameters from the DB
2424 //-----------------------------------------------------------------------------
2426 if(gSystem->AccessPathName(inName,kFileExists)) {
2427 cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
2430 TFile *inFile = TFile::Open(inName);
2433 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
2434 fdEdxMeanPi.~AliTPCkineGrid();
2435 new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
2436 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
2437 fdEdxRMSPi.~AliTPCkineGrid();
2438 new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
2441 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
2442 fdEdxMeanKa.~AliTPCkineGrid();
2443 new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
2444 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
2445 fdEdxRMSKa.~AliTPCkineGrid();
2446 new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
2449 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
2450 fdEdxMeanPr.~AliTPCkineGrid();
2451 new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
2452 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
2453 fdEdxRMSPr.~AliTPCkineGrid();
2454 new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
2457 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
2458 fdEdxMeanEl.~AliTPCkineGrid();
2459 new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
2460 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
2461 fdEdxRMSEl.~AliTPCkineGrid();
2462 new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
2466 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
2467 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
2469 fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxMeanMu");
2470 fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxRMSMu");
2472 fdEdxMeanMu.~AliTPCkineGrid();
2473 new(&fdEdxMeanMu) AliTPCkineGrid(*fdEdxMean);
2474 fdEdxRMSMu.~AliTPCkineGrid();
2475 new(&fdEdxRMSMu) AliTPCkineGrid(*fdEdxRMS);
2482 //-----------------------------------------------------------------------------
2483 Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
2484 //-----------------------------------------------------------------------------
2485 // This function reads the kine grid from the DB
2486 //-----------------------------------------------------------------------------
2488 if(gSystem->AccessPathName(inName,kFileExists)) {
2489 cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n";
2492 TFile *inFile = TFile::Open(inName);
2494 // first read the DB grid for the different particles
2495 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2496 fDBgridPi.~AliTPCkineGrid();
2497 new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
2498 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
2499 fDBgridKa.~AliTPCkineGrid();
2500 new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
2502 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2504 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Protons/DBgridPr");
2506 fDBgridPr.~AliTPCkineGrid();
2507 new(&fDBgridPr) AliTPCkineGrid(*fDBgrid);
2508 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
2509 fDBgridEl.~AliTPCkineGrid();
2510 new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
2512 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2514 fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Muons/DBgridMu");
2516 fDBgridMu.~AliTPCkineGrid();
2517 new(&fDBgridMu) AliTPCkineGrid(*fDBgrid);
2523 //-----------------------------------------------------------------------------
2524 Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
2525 //-----------------------------------------------------------------------------
2526 // This function reads the TPC efficiencies from the DB
2527 //-----------------------------------------------------------------------------
2529 if(gSystem->AccessPathName(inName,kFileExists)) {
2530 cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n";
2533 TFile *inFile = TFile::Open(inName);
2535 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
2536 fEffPi.~AliTPCkineGrid();
2537 new(&fEffPi) AliTPCkineGrid(*fEff);
2538 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
2539 fEffKa.~AliTPCkineGrid();
2540 new(&fEffKa) AliTPCkineGrid(*fEff);
2541 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
2542 fEffPr.~AliTPCkineGrid();
2543 new(&fEffPr) AliTPCkineGrid(*fEff);
2544 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
2545 fEffEl.~AliTPCkineGrid();
2546 new(&fEffEl) AliTPCkineGrid(*fEff);
2547 fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
2548 fEffMu.~AliTPCkineGrid();
2549 new(&fEffMu) AliTPCkineGrid(*fEff);
2555 //-----------------------------------------------------------------------------
2556 Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
2557 //-----------------------------------------------------------------------------
2558 // This function reads the pulls from the DB
2559 //-----------------------------------------------------------------------------
2561 if(gSystem->AccessPathName(inName,kFileExists)) {
2562 cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n";
2565 TFile *inFile = TFile::Open(inName);
2567 for(Int_t i=0; i<5; i++) {
2568 TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
2569 TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
2570 TString pr("/Pulls/Protons/PullsPr_"); pr+=i;
2571 TString el("/Pulls/Electrons/PullsEl_"); el+=i;
2572 TString mu("/Pulls/Muons/PullsMu_"); mu+=i;
2574 fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2575 fPullsPi[i].~AliTPCkineGrid();
2576 new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
2578 fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
2579 fPullsKa[i].~AliTPCkineGrid();
2580 new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
2583 fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2585 fPulls = (AliTPCkineGrid*)inFile->Get(pr.Data());
2587 fPullsPr[i].~AliTPCkineGrid();
2588 new(&fPullsPr[i]) AliTPCkineGrid(*fPulls);
2590 fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
2591 fPullsEl[i].~AliTPCkineGrid();
2592 new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
2595 fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2597 fPulls = (AliTPCkineGrid*)inFile->Get(mu.Data());
2599 fPullsMu[i].~AliTPCkineGrid();
2600 new(&fPullsMu[i]) AliTPCkineGrid(*fPulls);
2607 //-----------------------------------------------------------------------------
2608 Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
2609 //-----------------------------------------------------------------------------
2610 // This function reads the regularization parameters from the DB
2611 //-----------------------------------------------------------------------------
2613 if(gSystem->AccessPathName(inName,kFileExists)) {
2614 cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
2616 // Introduced change to "reverse" the TMatrixD read from file.
2617 // Needed because storage mode of TMatrixD changed from column-wise
2618 // to rwo-wise in ROOT.
2620 // A.Dainese 03/06/05
2622 TMatrixD dummyMatrix(9,3);
2624 TFile *inFile = TFile::Open(inName);
2627 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2628 new(&fRegParPi) TMatrixD(*fRegPar);
2629 //printf("before\n");
2630 //for(Int_t jj=0;jj<9;jj++) printf("%13.10f %13.10f %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2));
2631 dummyMatrix = fRegParPi;
2632 fRegParPi(0,0) = dummyMatrix(0,0);
2633 fRegParPi(0,1) = dummyMatrix(0,1);
2634 fRegParPi(0,2) = dummyMatrix(0,2);
2635 fRegParPi(1,0) = dummyMatrix(3,0);
2636 fRegParPi(1,1) = dummyMatrix(1,1);
2637 fRegParPi(1,2) = dummyMatrix(1,2);
2638 fRegParPi(2,0) = dummyMatrix(6,0);
2639 fRegParPi(2,1) = dummyMatrix(3,2);
2640 fRegParPi(2,2) = dummyMatrix(2,2);
2641 fRegParPi(3,0) = dummyMatrix(1,0);
2642 fRegParPi(3,1) = dummyMatrix(4,0);
2643 fRegParPi(3,2) = dummyMatrix(7,0);
2644 fRegParPi(4,0) = dummyMatrix(3,1);
2645 fRegParPi(4,1) = dummyMatrix(4,1);
2646 fRegParPi(4,2) = dummyMatrix(7,1);
2647 fRegParPi(5,0) = dummyMatrix(6,1);
2648 fRegParPi(5,1) = dummyMatrix(4,2);
2649 fRegParPi(5,2) = dummyMatrix(7,2);
2650 fRegParPi(6,0) = dummyMatrix(2,0);
2651 fRegParPi(6,1) = dummyMatrix(5,0);
2652 fRegParPi(6,2) = dummyMatrix(8,0);
2653 fRegParPi(7,0) = dummyMatrix(2,1);
2654 fRegParPi(7,1) = dummyMatrix(5,1);
2655 fRegParPi(7,2) = dummyMatrix(8,1);
2656 fRegParPi(8,0) = dummyMatrix(6,2);
2657 fRegParPi(8,1) = dummyMatrix(5,2);
2658 fRegParPi(8,2) = dummyMatrix(8,2);
2659 //printf("after\n");
2660 //for(Int_t jj=0;jj<9;jj++) printf("%13.10f %13.10f %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2));
2663 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
2664 new(&fRegParKa) TMatrixD(*fRegPar);
2665 dummyMatrix = fRegParKa;
2666 fRegParKa(0,0) = dummyMatrix(0,0);
2667 fRegParKa(0,1) = dummyMatrix(0,1);
2668 fRegParKa(0,2) = dummyMatrix(0,2);
2669 fRegParKa(1,0) = dummyMatrix(3,0);
2670 fRegParKa(1,1) = dummyMatrix(1,1);
2671 fRegParKa(1,2) = dummyMatrix(1,2);
2672 fRegParKa(2,0) = dummyMatrix(6,0);
2673 fRegParKa(2,1) = dummyMatrix(3,2);
2674 fRegParKa(2,2) = dummyMatrix(2,2);
2675 fRegParKa(3,0) = dummyMatrix(1,0);
2676 fRegParKa(3,1) = dummyMatrix(4,0);
2677 fRegParKa(3,2) = dummyMatrix(7,0);
2678 fRegParKa(4,0) = dummyMatrix(3,1);
2679 fRegParKa(4,1) = dummyMatrix(4,1);
2680 fRegParKa(4,2) = dummyMatrix(7,1);
2681 fRegParKa(5,0) = dummyMatrix(6,1);
2682 fRegParKa(5,1) = dummyMatrix(4,2);
2683 fRegParKa(5,2) = dummyMatrix(7,2);
2684 fRegParKa(6,0) = dummyMatrix(2,0);
2685 fRegParKa(6,1) = dummyMatrix(5,0);
2686 fRegParKa(6,2) = dummyMatrix(8,0);
2687 fRegParKa(7,0) = dummyMatrix(2,1);
2688 fRegParKa(7,1) = dummyMatrix(5,1);
2689 fRegParKa(7,2) = dummyMatrix(8,1);
2690 fRegParKa(8,0) = dummyMatrix(6,2);
2691 fRegParKa(8,1) = dummyMatrix(5,2);
2692 fRegParKa(8,2) = dummyMatrix(8,2);
2696 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2698 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
2700 new(&fRegParPr) TMatrixD(*fRegPar);
2701 dummyMatrix = fRegParPr;
2702 fRegParPr(0,0) = dummyMatrix(0,0);
2703 fRegParPr(0,1) = dummyMatrix(0,1);
2704 fRegParPr(0,2) = dummyMatrix(0,2);
2705 fRegParPr(1,0) = dummyMatrix(3,0);
2706 fRegParPr(1,1) = dummyMatrix(1,1);
2707 fRegParPr(1,2) = dummyMatrix(1,2);
2708 fRegParPr(2,0) = dummyMatrix(6,0);
2709 fRegParPr(2,1) = dummyMatrix(3,2);
2710 fRegParPr(2,2) = dummyMatrix(2,2);
2711 fRegParPr(3,0) = dummyMatrix(1,0);
2712 fRegParPr(3,1) = dummyMatrix(4,0);
2713 fRegParPr(3,2) = dummyMatrix(7,0);
2714 fRegParPr(4,0) = dummyMatrix(3,1);
2715 fRegParPr(4,1) = dummyMatrix(4,1);
2716 fRegParPr(4,2) = dummyMatrix(7,1);
2717 fRegParPr(5,0) = dummyMatrix(6,1);
2718 fRegParPr(5,1) = dummyMatrix(4,2);
2719 fRegParPr(5,2) = dummyMatrix(7,2);
2720 fRegParPr(6,0) = dummyMatrix(2,0);
2721 fRegParPr(6,1) = dummyMatrix(5,0);
2722 fRegParPr(6,2) = dummyMatrix(8,0);
2723 fRegParPr(7,0) = dummyMatrix(2,1);
2724 fRegParPr(7,1) = dummyMatrix(5,1);
2725 fRegParPr(7,2) = dummyMatrix(8,1);
2726 fRegParPr(8,0) = dummyMatrix(6,2);
2727 fRegParPr(8,1) = dummyMatrix(5,2);
2728 fRegParPr(8,2) = dummyMatrix(8,2);
2731 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
2732 new(&fRegParEl) TMatrixD(*fRegPar);
2733 dummyMatrix = fRegParEl;
2734 fRegParEl(0,0) = dummyMatrix(0,0);
2735 fRegParEl(0,1) = dummyMatrix(0,1);
2736 fRegParEl(0,2) = dummyMatrix(0,2);
2737 fRegParEl(1,0) = dummyMatrix(3,0);
2738 fRegParEl(1,1) = dummyMatrix(1,1);
2739 fRegParEl(1,2) = dummyMatrix(1,2);
2740 fRegParEl(2,0) = dummyMatrix(6,0);
2741 fRegParEl(2,1) = dummyMatrix(3,2);
2742 fRegParEl(2,2) = dummyMatrix(2,2);
2743 fRegParEl(3,0) = dummyMatrix(1,0);
2744 fRegParEl(3,1) = dummyMatrix(4,0);
2745 fRegParEl(3,2) = dummyMatrix(7,0);
2746 fRegParEl(4,0) = dummyMatrix(3,1);
2747 fRegParEl(4,1) = dummyMatrix(4,1);
2748 fRegParEl(4,2) = dummyMatrix(7,1);
2749 fRegParEl(5,0) = dummyMatrix(6,1);
2750 fRegParEl(5,1) = dummyMatrix(4,2);
2751 fRegParEl(5,2) = dummyMatrix(7,2);
2752 fRegParEl(6,0) = dummyMatrix(2,0);
2753 fRegParEl(6,1) = dummyMatrix(5,0);
2754 fRegParEl(6,2) = dummyMatrix(8,0);
2755 fRegParEl(7,0) = dummyMatrix(2,1);
2756 fRegParEl(7,1) = dummyMatrix(5,1);
2757 fRegParEl(7,2) = dummyMatrix(8,1);
2758 fRegParEl(8,0) = dummyMatrix(6,2);
2759 fRegParEl(8,1) = dummyMatrix(5,2);
2760 fRegParEl(8,2) = dummyMatrix(8,2);
2764 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2766 fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
2768 new(&fRegParMu) TMatrixD(*fRegPar);
2769 dummyMatrix = fRegParMu;
2770 fRegParMu(0,0) = dummyMatrix(0,0);
2771 fRegParMu(0,1) = dummyMatrix(0,1);
2772 fRegParMu(0,2) = dummyMatrix(0,2);
2773 fRegParMu(1,0) = dummyMatrix(3,0);
2774 fRegParMu(1,1) = dummyMatrix(1,1);
2775 fRegParMu(1,2) = dummyMatrix(1,2);
2776 fRegParMu(2,0) = dummyMatrix(6,0);
2777 fRegParMu(2,1) = dummyMatrix(3,2);
2778 fRegParMu(2,2) = dummyMatrix(2,2);
2779 fRegParMu(3,0) = dummyMatrix(1,0);
2780 fRegParMu(3,1) = dummyMatrix(4,0);
2781 fRegParMu(3,2) = dummyMatrix(7,0);
2782 fRegParMu(4,0) = dummyMatrix(3,1);
2783 fRegParMu(4,1) = dummyMatrix(4,1);
2784 fRegParMu(4,2) = dummyMatrix(7,1);
2785 fRegParMu(5,0) = dummyMatrix(6,1);
2786 fRegParMu(5,1) = dummyMatrix(4,2);
2787 fRegParMu(5,2) = dummyMatrix(7,2);
2788 fRegParMu(6,0) = dummyMatrix(2,0);
2789 fRegParMu(6,1) = dummyMatrix(5,0);
2790 fRegParMu(6,2) = dummyMatrix(8,0);
2791 fRegParMu(7,0) = dummyMatrix(2,1);
2792 fRegParMu(7,1) = dummyMatrix(5,1);
2793 fRegParMu(7,2) = dummyMatrix(8,1);
2794 fRegParMu(8,0) = dummyMatrix(6,2);
2795 fRegParMu(8,1) = dummyMatrix(5,2);
2796 fRegParMu(8,2) = dummyMatrix(8,2);
2803 //-----------------------------------------------------------------------------
2804 void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
2805 //-----------------------------------------------------------------------------
2806 // This function regularizes the elements of the covariance matrix
2807 // that show a momentum depence:
2808 // c00, c11, c22, c33, c44, c20, c24, c40, c31
2809 // The regularization is done separately for pions, kaons, electrons:
2810 // give "Pion","Kaon" or "Electron" as first argument.
2811 //-----------------------------------------------------------------------------
2813 gStyle->SetOptStat(0);
2814 gStyle->SetOptFit(10001);
2817 const char *part="Pions - ";
2819 InitializeKineGrid("DB");
2821 const Int_t kfitbins = fDBgrid->GetBinsPt();
2822 cerr<<" Fit bins: "<<kfitbins<<endl;
2828 cerr<<" Processing pions ...\n";
2833 cerr<<" Processing kaons ...\n";
2835 case 2212: //protons
2838 cerr<<" Processing protons ...\n";
2840 case 11: // electrons
2842 part="Electrons - ";
2843 cerr<<" Processing electrons ...\n";
2848 cerr<<" Processing muons ...\n";
2854 // create a chain with compared tracks
2855 TChain *cmptrkchain = new ("cmptrktree");
2856 cmptrkchain.Add("CovMatrix_AllEvts.root");
2857 //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
2858 //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
2859 //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
2863 TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
2864 TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
2867 cmptrktree->SetBranchAddress("comptracks",&cmptrk);
2868 Int_t entries = (Int_t)cmptrktree->GetEntries();
2872 Int_t *n = new Int_t[kfitbins];
2873 Int_t *n00 = new Int_t[kfitbins];
2874 Int_t *n11 = new Int_t[kfitbins];
2875 Int_t *n20 = new Int_t[kfitbins];
2876 Int_t *n22 = new Int_t[kfitbins];
2877 Int_t *n31 = new Int_t[kfitbins];
2878 Int_t *n33 = new Int_t[kfitbins];
2879 Int_t *n40 = new Int_t[kfitbins];
2880 Int_t *n42 = new Int_t[kfitbins];
2881 Int_t *n44 = new Int_t[kfitbins];
2882 Double_t *p = new Double_t[kfitbins];
2883 Double_t *ep = new Double_t[kfitbins];
2884 Double_t *mean00 = new Double_t[kfitbins];
2885 Double_t *mean11 = new Double_t[kfitbins];
2886 Double_t *mean20 = new Double_t[kfitbins];
2887 Double_t *mean22 = new Double_t[kfitbins];
2888 Double_t *mean31 = new Double_t[kfitbins];
2889 Double_t *mean33 = new Double_t[kfitbins];
2890 Double_t *mean40 = new Double_t[kfitbins];
2891 Double_t *mean42 = new Double_t[kfitbins];
2892 Double_t *mean44 = new Double_t[kfitbins];
2893 Double_t *sigma00 = new Double_t[kfitbins];
2894 Double_t *sigma11 = new Double_t[kfitbins];
2895 Double_t *sigma20 = new Double_t[kfitbins];
2896 Double_t *sigma22 = new Double_t[kfitbins];
2897 Double_t *sigma31 = new Double_t[kfitbins];
2898 Double_t *sigma33 = new Double_t[kfitbins];
2899 Double_t *sigma40 = new Double_t[kfitbins];
2900 Double_t *sigma42 = new Double_t[kfitbins];
2901 Double_t *sigma44 = new Double_t[kfitbins];
2902 Double_t *rmean = new Double_t[kfitbins];
2903 Double_t *rsigma = new Double_t[kfitbins];
2906 for(Int_t l=0; l<kfitbins; l++) {
2908 n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
2910 mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
2911 sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
2914 // loop on chain entries for mean
2915 for(Int_t l=0; l<entries; l++) {
2916 cmptrktree->GetEvent(l);
2917 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2918 pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2921 mean00[pbin]+=cmptrk.c00;
2922 mean11[pbin]+=cmptrk.c11;
2923 mean20[pbin]+=cmptrk.c20;
2924 mean22[pbin]+=cmptrk.c22;
2925 mean31[pbin]+=cmptrk.c31;
2926 mean33[pbin]+=cmptrk.c33;
2927 mean40[pbin]+=cmptrk.c40;
2928 mean42[pbin]+=cmptrk.c42;
2929 mean44[pbin]+=cmptrk.c44;
2930 } // loop on chain entries
2932 for(Int_t l=0; l<kfitbins; l++) {
2945 // loop on chain entries for sigma
2946 for(Int_t l=0; l<entries; l++) {
2947 cmptrktree->GetEvent(l);
2948 if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2949 pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2950 if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2951 sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2952 if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2953 sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2954 if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2955 sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2956 if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2957 sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2958 if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2959 sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2960 if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2961 sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2962 if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2963 sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2964 if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2965 sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2966 if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2967 sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2968 } // loop on chain entries
2970 for(Int_t l=0; l<kfitbins; l++) {
2971 sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2972 sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2973 sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2974 sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2975 sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2976 sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2977 sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2978 sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2979 sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2984 TF1 *func = new TF1("RegFunc",RegFunc,0.23,50.,3);
2985 func->SetParNames("A_meas","A_scatt","B");
2987 // line to draw on the plots
2988 TLine *lin = new TLine(-1,1,1.69,1);
2989 lin->SetLineStyle(2);
2990 lin->SetLineWidth(2);
2992 // matrix used to store fit results
2993 TMatrixD fitRes(9,3);
2997 // create the canvas
2998 TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
2999 canv00->Divide(1,2);
3000 // create the graph for cov matrix
3001 TGraphErrors *gr00 = new TGraphErrors(kfitbins,p,mean00,ep,sigma00);
3002 TString title00("C(y,y)"); title00.Prepend(part);
3003 TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
3004 frame00->SetXTitle("p [GeV/c]");
3005 canv00->cd(1); gPad->SetLogx();
3008 // Sets initial values for parameters
3009 func->SetParameters(1.6e-3,1.9e-4,1.5);
3010 // Fit points in range defined by function
3011 gr00->Fit("RegFunc","R,Q");
3012 func->GetParameters(fitpar);
3013 for(Int_t i=0; i<3; i++) fitRes(0,i)=fitpar[i];
3014 for(Int_t l=0; l<kfitbins; l++) {
3015 rmean[l] = mean00[l]/RegFunc(&p[l],fitpar);
3016 rsigma[l] = sigma00[l]/RegFunc(&p[l],fitpar);
3018 // create the graph the regularized cov. matrix
3019 TGraphErrors *gr00reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3020 TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B})");
3021 regtitle00.Prepend(part);
3022 TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
3023 frame00reg->SetXTitle("p [GeV/c]");
3024 canv00->cd(2); gPad->SetLogx();
3032 // create the canvas
3033 TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
3034 canv11->Divide(1,2);
3035 // create the graph for cov matrix
3036 TGraphErrors *gr11 = new TGraphErrors(kfitbins,p,mean11,ep,sigma11);
3037 TString title11("C(z,z)"); title11.Prepend(part);
3038 TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
3039 frame11->SetXTitle("p [GeV/c]");
3040 canv11->cd(1); gPad->SetLogx();
3043 // Sets initial values for parameters
3044 func->SetParameters(1.2e-3,8.1e-4,1.);
3045 // Fit points in range defined by function
3046 gr11->Fit("RegFunc","R,Q");
3047 func->GetParameters(fitpar);
3048 for(Int_t i=0; i<3; i++) fitRes(1,i)=fitpar[i];
3049 for(Int_t l=0; l<kfitbins; l++) {
3050 rmean[l] = mean11[l]/RegFunc(&p[l],fitpar);
3051 rsigma[l] = sigma11[l]/RegFunc(&p[l],fitpar);
3053 // create the graph the regularized cov. matrix
3054 TGraphErrors *gr11reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3055 TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B})");
3056 regtitle11.Prepend(part);
3057 TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
3058 frame11reg->SetXTitle("p [GeV/c]");
3059 canv11->cd(2); gPad->SetLogx();
3067 // create the canvas
3068 TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
3069 canv20->Divide(1,2);
3070 // create the graph for cov matrix
3071 TGraphErrors *gr20 = new TGraphErrors(kfitbins,p,mean20,ep,sigma20);
3072 TString title20("C(#eta, y)"); title20.Prepend(part);
3073 TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
3074 frame20->SetXTitle("p [GeV/c]");
3075 canv20->cd(1); gPad->SetLogx();
3078 // Sets initial values for parameters
3079 func->SetParameters(7.3e-5,1.2e-5,1.5);
3080 // Fit points in range defined by function
3081 gr20->Fit("RegFunc","R,Q");
3082 func->GetParameters(fitpar);
3083 for(Int_t i=0; i<3; i++) fitRes(2,i)=fitpar[i];
3084 for(Int_t l=0; l<kfitbins; l++) {
3085 rmean[l] = mean20[l]/RegFunc(&p[l],fitpar);
3086 rsigma[l] = sigma20[l]/RegFunc(&p[l],fitpar);
3088 // create the graph the regularized cov. matrix
3089 TGraphErrors *gr20reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3090 TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B})");
3091 regtitle20.Prepend(part);
3092 TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
3093 frame20reg->SetXTitle("p [GeV/c]");
3094 canv20->cd(2); gPad->SetLogx();
3102 // create the canvas
3103 TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
3104 canv22->Divide(1,2);
3105 // create the graph for cov matrix
3106 TGraphErrors *gr22 = new TGraphErrors(kfitbins,p,mean22,ep,sigma22);
3107 TString title22("C(#eta, #eta)"); title22.Prepend(part);
3108 TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
3109 frame22->SetXTitle("p [GeV/c]");
3110 canv22->cd(1); gPad->SetLogx();
3113 // Sets initial values for parameters
3114 func->SetParameters(5.2e-6,1.1e-6,2.);
3115 // Fit points in range defined by function
3116 gr22->Fit("RegFunc","R,Q");
3117 func->GetParameters(fitpar);
3118 for(Int_t i=0; i<3; i++) fitRes(3,i)=fitpar[i];
3119 for(Int_t l=0; l<kfitbins; l++) {
3120 rmean[l] = mean22[l]/RegFunc(&p[l],fitpar);
3121 rsigma[l] = sigma22[l]/RegFunc(&p[l],fitpar);
3123 // create the graph the regularized cov. matrix
3124 TGraphErrors *gr22reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3125 TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B})");
3126 regtitle22.Prepend(part);
3127 TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
3128 frame22reg->SetXTitle("p [GeV/c]");
3129 canv22->cd(2); gPad->SetLogx();
3137 // create the canvas
3138 TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
3139 canv31->Divide(1,2);
3140 // create the graph for cov matrix
3141 TGraphErrors *gr31 = new TGraphErrors(kfitbins,p,mean31,ep,sigma31);
3142 TString title31("C(tg #lambda,z)"); title31.Prepend(part);
3143 TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
3144 frame31->SetXTitle("p [GeV/c]");
3145 canv31->cd(1); gPad->SetLogx();
3148 // Sets initial values for parameters
3149 func->SetParameters(-1.2e-5,-1.2e-5,1.5);
3150 // Fit points in range defined by function
3151 gr31->Fit("RegFunc","R,Q");
3152 func->GetParameters(fitpar);
3153 for(Int_t i=0; i<3; i++) fitRes(4,i)=fitpar[i];
3154 for(Int_t l=0; l<kfitbins; l++) {
3155 rmean[l] = mean31[l]/RegFunc(&p[l],fitpar);
3156 rsigma[l] = -sigma31[l]/RegFunc(&p[l],fitpar);
3158 // create the graph the regularized cov. matrix
3159 TGraphErrors *gr31reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3160 TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B})");
3161 regtitle31.Prepend(part);
3162 TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
3163 frame31reg->SetXTitle("p [GeV/c]");
3164 canv31->cd(2); gPad->SetLogx();
3172 // create the canvas
3173 TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
3174 canv33->Divide(1,2);
3175 // create the graph for cov matrix
3176 TGraphErrors *gr33 = new TGraphErrors(kfitbins,p,mean33,ep,sigma33);
3177 TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
3178 TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
3179 frame33->SetXTitle("p [GeV/c]");
3180 canv33->cd(1); gPad->SetLogx();
3183 // Sets initial values for parameters
3184 func->SetParameters(1.3e-7,4.6e-7,1.7);
3185 // Fit points in range defined by function
3186 gr33->Fit("RegFunc","R,Q");
3187 func->GetParameters(fitpar);
3188 for(Int_t i=0; i<3; i++) fitRes(5,i)=fitpar[i];
3189 for(Int_t l=0; l<kfitbins; l++) {
3190 rmean[l] = mean33[l]/RegFunc(&p[l],fitpar);
3191 rsigma[l] = sigma33[l]/RegFunc(&p[l],fitpar);
3193 // create the graph the regularized cov. matrix
3194 TGraphErrors *gr33reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3195 TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B})");
3196 regtitle33.Prepend(part);
3197 TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
3198 frame33reg->SetXTitle("p [GeV/c]");
3199 canv33->cd(2); gPad->SetLogx();
3207 // create the canvas
3208 TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
3209 canv40->Divide(1,2);
3210 // create the graph for cov matrix
3211 TGraphErrors *gr40 = new TGraphErrors(kfitbins,p,mean40,ep,sigma40);
3212 TString title40("C(C,y)"); title40.Prepend(part);
3213 TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
3214 frame40->SetXTitle("p [GeV/c]");
3215 canv40->cd(1); gPad->SetLogx();
3218 // Sets initial values for parameters
3219 func->SetParameters(4.e-7,4.4e-8,1.5);
3220 // Fit points in range defined by function
3221 gr40->Fit("RegFunc","R,Q");
3222 func->GetParameters(fitpar);
3223 for(Int_t i=0; i<3; i++) fitRes(6,i)=fitpar[i];
3224 for(Int_t l=0; l<kfitbins; l++) {
3225 rmean[l] = mean40[l]/RegFunc(&p[l],fitpar);
3226 rsigma[l] = sigma40[l]/RegFunc(&p[l],fitpar);
3228 // create the graph the regularized cov. matrix
3229 TGraphErrors *gr40reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3230 TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B})");
3231 regtitle40.Prepend(part);
3232 TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
3233 frame40reg->SetXTitle("p [GeV/c]");
3234 canv40->cd(2); gPad->SetLogx();
3242 // create the canvas
3243 TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
3244 canv42->Divide(1,2);
3245 // create the graph for cov matrix
3246 TGraphErrors *gr42 = new TGraphErrors(kfitbins,p,mean42,ep,sigma42);
3247 TString title42("C(C, #eta)"); title42.Prepend(part);
3248 TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
3249 frame42->SetXTitle("p [GeV/c]");
3250 canv42->cd(1); gPad->SetLogx();
3253 // Sets initial values for parameters
3254 func->SetParameters(3.e-8,8.2e-9,2.);
3255 // Fit points in range defined by function
3256 gr42->Fit("RegFunc","R,Q");
3257 func->GetParameters(fitpar);
3258 for(Int_t i=0; i<3; i++) fitRes(7,i)=fitpar[i];
3259 for(Int_t l=0; l<kfitbins; l++) {
3260 rmean[l] = mean42[l]/RegFunc(&p[l],fitpar);
3261 rsigma[l] = sigma42[l]/RegFunc(&p[l],fitpar);
3263 // create the graph the regularized cov. matrix
3264 TGraphErrors *gr42reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3265 TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B})");
3266 regtitle42.Prepend(part);
3267 TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
3268 frame42reg->SetXTitle("p [GeV/c]");
3269 canv42->cd(2); gPad->SetLogx();
3277 // create the canvas
3278 TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
3279 canv44->Divide(1,2);
3280 // create the graph for cov matrix
3281 TGraphErrors *gr44 = new TGraphErrors(kfitbins,p,mean44,ep,sigma44);
3282 TString title44("C(C,C)"); title44.Prepend(part);
3283 TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
3284 frame44->SetXTitle("p [GeV/c]");
3285 canv44->cd(1); gPad->SetLogx();
3288 // Sets initial values for parameters
3289 func->SetParameters(1.8e-10,5.8e-11,2.);
3290 // Fit points in range defined by function
3291 gr44->Fit("RegFunc","R,Q");
3292 func->GetParameters(fitpar);
3293 for(Int_t i=0; i<3; i++) fitRes(8,i)=fitpar[i];
3294 for(Int_t l=0; l<kfitbins; l++) {
3295 rmean[l] = mean44[l]/RegFunc(&p[l],fitpar);
3296 rsigma[l] = sigma44[l]/RegFunc(&p[l],fitpar);
3298 // create the graph the regularized cov. matrix
3299 TGraphErrors *gr44reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
3300 TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B})");
3301 regtitle44.Prepend(part);
3302 TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
3303 frame44reg->SetXTitle("p [GeV/c]");
3304 canv44->cd(2); gPad->SetLogx();
3314 new(&fRegParPi) TMatrixD(fitRes);
3317 new(&fRegParKa) TMatrixD(fitRes);
3320 new(&fRegParPr) TMatrixD(fitRes);
3323 new(&fRegParEl) TMatrixD(fitRes);
3326 new(&fRegParMu) TMatrixD(fitRes);
3330 // write fit parameters to file
3331 WriteRegParams(outName,pdg);
3368 //-----------------------------------------------------------------------------
3369 Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
3370 //-----------------------------------------------------------------------------
3371 // This function makes a selection according to TPC tracking efficiency
3372 //-----------------------------------------------------------------------------
3376 eff = fEff->GetValueAt(pt,eta);
3378 if(gRandom->Rndm() < eff) return kTRUE;
3382 //-----------------------------------------------------------------------------
3383 void AliTPCtrackerParam::SetParticle(Int_t pdg) {
3384 //-----------------------------------------------------------------------------
3385 // This function set efficiencies, pulls, etc... for the particle
3386 // specie of the current particle
3387 //-----------------------------------------------------------------------------
3391 fDBgrid = &fDBgridPi;
3394 fRegPar = &fRegParPi;
3395 fdEdxMean = &fdEdxMeanPi;
3396 fdEdxRMS = &fdEdxRMSPi;
3399 fDBgrid = &fDBgridKa;
3402 fRegPar = &fRegParKa;
3403 fdEdxMean = &fdEdxMeanKa;
3404 fdEdxRMS = &fdEdxRMSKa;
3407 fDBgrid = &fDBgridPr;
3410 fRegPar = &fRegParPr;
3411 fdEdxMean = &fdEdxMeanPr;
3412 fdEdxRMS = &fdEdxRMSPr;
3415 fDBgrid = &fDBgridEl;
3418 fRegPar = &fRegParEl;
3419 fdEdxMean = &fdEdxMeanEl;
3420 fdEdxRMS = &fdEdxRMSEl;
3423 fDBgrid = &fDBgridMu;
3426 fRegPar = &fRegParMu;
3427 fdEdxMean = &fdEdxMeanMu;
3428 fdEdxRMS = &fdEdxRMSMu;
3434 //-----------------------------------------------------------------------------
3435 void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
3437 //-----------------------------------------------------------------------------
3438 // This function smears track parameters according to streched cov. matrix
3439 //-----------------------------------------------------------------------------
3440 Double_t xref=xxsm[0]; xxsm[0]=0;
3442 AliGausCorr *corgen = new AliGausCorr(cov,5);
3444 corgen->GetGaussN(corr);
3445 // check on fP2(ESD) = fX*fP4-fP2(TPC)
3446 // with fX=xref (not smeared), fP4=xx[4]+corr[4] e fP2=xx[2]+corr[2];
3447 // if fP2(ESD)^2 > 1 -> resmear...
3448 Double_t fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]);
3449 while ( (fp2esd*fp2esd) > 1.0 ) {
3450 Warning("SmearTrack()","Badly smeared track, retrying...");
3451 corgen->GetGaussN(corr);
3452 fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]);
3458 for(Int_t l=0;l<5;l++) xxsm[l] = xx[l]+corr[l];
3462 //-----------------------------------------------------------------------------
3463 Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
3464 //-----------------------------------------------------------------------------
3465 // This function writes the dEdx parameters to the DB
3466 //-----------------------------------------------------------------------------
3469 const char *dirName="Pions";
3470 const char *meanName="dEdxMeanPi";
3471 const char *rmsName="dEdxRMSPi";
3475 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
3476 } else { opt="update"; }
3481 meanName="dEdxMeanPi";
3482 rmsName="dEdxRMSPi";
3486 meanName="dEdxMeanKa";
3487 rmsName="dEdxRMSKa";
3491 meanName="dEdxMeanPr";
3492 rmsName="dEdxRMSPr";
3495 dirName="Electrons";
3496 meanName="dEdxMeanEl";
3497 rmsName="dEdxRMSEl";
3501 meanName="dEdxMeanMu";
3502 rmsName="dEdxRMSMu";
3506 TFile *outFile = new TFile(outName,opt);
3507 if(!gDirectory->cd("/dEdx")) {
3508 outFile->mkdir("dEdx");
3509 gDirectory->cd("/dEdx");
3511 TDirectory *dir2 = gDirectory->mkdir(dirName);
3513 fdEdxMean->SetName(meanName); fdEdxMean->Write();
3514 fdEdxRMS->SetName(rmsName); fdEdxRMS->Write();
3522 //-----------------------------------------------------------------------------
3523 Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
3524 //-----------------------------------------------------------------------------
3525 // This function writes the TPC efficiencies to the DB
3526 //-----------------------------------------------------------------------------
3531 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
3532 } else { opt="update"; }
3534 TFile *outFile = new TFile(outName,opt);
3536 outFile->mkdir("Efficiencies");
3537 gDirectory->cd("/Efficiencies");
3538 gDirectory->mkdir("Pions");
3539 gDirectory->mkdir("Kaons");
3540 gDirectory->mkdir("Protons");
3541 gDirectory->mkdir("Electrons");
3542 gDirectory->mkdir("Muons");
3544 gDirectory->cd("/Efficiencies/Pions");
3545 fEffPi.SetName("EffPi");
3547 gDirectory->cd("/Efficiencies/Kaons");
3548 fEffKa.SetName("EffKa");
3550 gDirectory->cd("/Efficiencies/Protons");
3551 fEffPr.SetName("EffPr");
3553 gDirectory->cd("/Efficiencies/Electrons");
3554 fEffEl.SetName("EffEl");
3556 gDirectory->cd("/Efficiencies/Muons");
3557 fEffMu.SetName("EffMu");
3566 //-----------------------------------------------------------------------------
3567 Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
3568 //-----------------------------------------------------------------------------
3569 // This function writes the pulls to the DB
3570 //-----------------------------------------------------------------------------
3574 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
3575 } else { opt="update"; }
3577 TFile *outFile = new TFile(outName,opt);
3579 outFile->mkdir("Pulls");
3580 gDirectory->cd("/Pulls");
3581 gDirectory->mkdir("Pions");
3582 gDirectory->mkdir("Kaons");
3583 gDirectory->mkdir("Protons");
3584 gDirectory->mkdir("Electrons");
3585 gDirectory->mkdir("Muons");
3587 for(Int_t i=0;i<5;i++) {
3588 TString pi("PullsPi_"); pi+=i;
3589 TString ka("PullsKa_"); ka+=i;
3590 TString pr("PullsPr_"); pr+=i;
3591 TString el("PullsEl_"); el+=i;
3592 TString mu("PullsMu_"); mu+=i;
3593 fPullsPi[i].SetName(pi.Data());
3594 fPullsKa[i].SetName(ka.Data());
3595 fPullsPr[i].SetName(pr.Data());
3596 fPullsEl[i].SetName(el.Data());
3597 fPullsMu[i].SetName(mu.Data());
3598 gDirectory->cd("/Pulls/Pions");
3599 fPullsPi[i].Write();
3600 gDirectory->cd("/Pulls/Kaons");
3601 fPullsKa[i].Write();
3602 gDirectory->cd("/Pulls/Protons");
3603 fPullsPr[i].Write();
3604 gDirectory->cd("/Pulls/Electrons");
3605 fPullsEl[i].Write();
3606 gDirectory->cd("/Pulls/Muons");
3607 fPullsMu[i].Write();
3614 //-----------------------------------------------------------------------------
3615 Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
3616 //-----------------------------------------------------------------------------
3617 // This function writes the regularization parameters to the DB
3618 //-----------------------------------------------------------------------------
3621 const char *dirName="Pions";
3622 const char *keyName="RegPions";
3626 if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate";
3627 } else { opt="update"; }
3640 keyName="RegProtons";
3643 dirName="Electrons";
3644 keyName="RegElectrons";
3652 TFile *outFile = new TFile(outName,opt);
3653 if(!gDirectory->cd("/RegParams")) {
3654 outFile->mkdir("RegParams");
3655 gDirectory->cd("/RegParams");
3657 TDirectory *dir2 = gDirectory->mkdir(dirName);
3659 fRegPar->Write(keyName);
3667 //-----------------------------------------------------------------------------
3668 //*****************************************************************************
3669 //-----------------------------------------------------------------------------