TPC/AliTPCtrackerParam.cxx

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