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