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 /* $Id$ */
  17
  18 /**************************************************************************
  19  *                                                                        *
  20  * This class builds AliTPCtrack objects from generated tracks to feed    *
  21  * ITS tracking (V2). The AliTPCtrack is built from its first hit in      *
  22  * the TPC. The track is assigned a Kalman-like covariance matrix         *
  23  * depending on its pT and pseudorapidity and track parameters are        *
  24  * smeared according to this covariance matrix.                           *
  25  * Output file contains sorted tracks, ready for matching with ITS.       *
  26  *                                                                        *
  27  * For details:                                                           *
  28  * Alice Internal Note (submitted - get it from andrea.dainese@pd.infn.it)*
  29  * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm    *
  30  *                                                                        *
  31  * Test macro is: AliBarrelRec_TPCparam.C                                 *
  32  *                                                                        *
  33  * 2002/10/01: Introduction of the database for pp collisions (B=0.4 T)   *
  34  * - Everything done separately for pions, kaons, protons, electrons and  *
  35  *   muons.                                                               *
  36  * - Now (only for pp) the tracks are built from the AliTrackReferences   *
  37  *   which contain the position and momentum of all tracks at R = 83 cm;  *
  38  *   This eliminates the loss of tracks due the dead zone of the TPC      *
  39  *   where the 1st hit is not created.                                    *
  40  * - In AliBarrelRec_TPCparam.C there many possible ways of determining   *
  41  *   the z coordinate of the primary vertex in pp events (from pixels,    *
  42  *   from ITS tracks, smearing according to resolution given by tracks.   *
  43  *                                                                        *
  44  * 2002/04/28: Major upgrade of the class                                 *
  45  * - Covariance matrices and pulls are now separeted for pions, kaons     *
  46  *   and electrons.                                                       *
  47  * - A parameterization of dE/dx in the TPC has been included and it is   *
  48  *   used to assign a mass to each track according to a rough dE/dx PID.  *
  49  * - All the "numbers" have been moved to the file with the DataBase,     *
  50  *   they are read as objects of the class AliTPCkineGrid, and assigned   *
  51  *   to data memebers of the class AliTPCtrackerParam.                    *
  52  * - All the code necessary to create a BataBase has been included in     *
  53  *   class (see the macro AliTPCtrackingParamDB.C for the details).       *
  54  *                                                                        *
  55  *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
  56  *                                                                        *
  57  **************************************************************************/
  58 //------- Root headers --------
  59 #include <TChain.h>
  60 #include <TF1.h>
  61 #include <TGraph.h>
  62 #include <TGraphErrors.h>
  63 #include <TLegend.h>
  64 #include <TLine.h>
  65 #include <TMatrixD.h>
  66 #include <TNtuple.h>
  67 #include <TSystem.h>
  68 #include <TH2F.h>
  69 //------ AliRoot headers ------
  70 #include "alles.h"
  71 #include "AliGausCorr.h"
  72 #include "AliKalmanTrack.h"
  73 #include "AliMagF.h"
  74 #include "AliMagFCM.h"
  75 #include "AliRunLoader.h"
  76 #include "AliTPCLoader.h"
  77 #include "AliTPCkineGrid.h"
  78 #include "AliTPCtrack.h"
  79 #include "AliTrackReference.h"
  80 #include "AliTPCtrackerParam.h"
  81 //-----------------------------
  82
  83 Double_t RegFunc(Double_t *x,Double_t *par) {
  84 // This is the function used to regularize the covariance matrix
  85   Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2]);
  86   return value;
  87 }
  88
  89 // structure for DB building
  90 typedef struct {
  91   Int_t    pdg;
  92   Int_t    bin;
  93   Double_t r;
  94   Double_t p;
  95   Double_t pt;
  96   Double_t cosl;
  97   Double_t eta;
  98   Double_t dpt;
  99   Double_t dP0,dP1,dP2,dP3,dP4;
 100   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
 101   Double_t dEdx;
 102 } COMPTRACK;
 103 // cov matrix structure
 104 typedef struct {
 105   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
 106 } COVMATRIX;
 107
 108 ClassImp(AliTPCtrackerParam)
 109
 110 //-----------------------------------------------------------------------------
 111 AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz,
 112                                        const Int_t n, const char* evfoldname):
 113   fEvFolderName(evfoldname) {
 114 //-----------------------------------------------------------------------------
 115 // This is the class conctructor
 116 //-----------------------------------------------------------------------------
 117
 118   fNevents = n;         // events to be processed
 119   fBz = Bz;             // value of the z component of L3 field (Tesla)
 120   fColl = coll;         // collision code (0: PbPb6000; 1: pp)
 121   fSelAndSmear = kTRUE; // by default selection and smearing are done
 122
 123   if(fBz!=0.4) {
 124     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid field!\n";
 125     cerr<<"      Available:  0.4"<<endl;
 126   }
 127   if(fColl!=0 && fColl!=1) {
 128     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid collision!\n";
 129     cerr<<"      Available:  0   ->   PbPb6000"<<endl;
 130     cerr<<"                  1   ->   pp"<<endl;
 131   }
 132
 133   fDBfileName = gSystem->Getenv("ALICE_ROOT");
 134   fDBfileName.Append("/TPC/CovMatrixDB_");
 135   //fDBfileName = "CovMatrixDB_";
 136   if(fColl==0) fDBfileName.Append("PbPb6000");
 137   if(fColl==1) fDBfileName.Append("pp");
 138   if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
 139 }
 140 //-----------------------------------------------------------------------------
 141 AliTPCtrackerParam::~AliTPCtrackerParam() {}
 142 //-----------------------------------------------------------------------------
 143 Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out) {
 144 //-----------------------------------------------------------------------------
 145 // This function creates the TPC parameterized tracks
 146 //-----------------------------------------------------------------------------
 147
 148   Error("BuildTPCtracks","in and out parameters ignored. new io");
 149
 150 /********************************************/
 151   AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
 152   if (rl == 0x0)
 153    {
 154      Error("BuildTPCtracks","Can not get Run Loader from event folder named %s.",
 155            fEvFolderName.Data());
 156      return 2;
 157    }
 158   AliLoader* tpcloader = rl->GetLoader("TPCLoader");
 159   if (tpcloader == 0x0)
 160    {
 161      Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
 162      return 3;
 163    }
 164
 165 /********************************************/
 166
 167   TFile *fileDB=0;
 168   TTree *covTreePi[50];
 169   TTree *covTreeKa[50];
 170   TTree *covTreePr[50];
 171   TTree *covTreeEl[50];
 172   TTree *covTreeMu[50];
 173
 174   if(fSelAndSmear) {
 175     cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
 176       fDBfileName.Data()<<"\n+++\n";
 177     // Read paramters from DB file
 178     if(!ReadAllData(fDBfileName.Data())) {
 179       cerr<<"AliTPCtrackerParam::BuildTPCtracks: \
 180              Could not read data from DB\n\n"; return 1;
 181     }
 182
 183     // Read the trees with regularized cov. matrices from DB
 184     TString str;
 185     fileDB = TFile::Open(fDBfileName.Data());
 186     Int_t nBinsPi = fDBgridPi.GetTotBins();
 187     for(Int_t l=0; l<nBinsPi; l++) {
 188       str = "/CovMatrices/Pions/CovTreePi_bin";
 189       str+=l;
 190       covTreePi[l] = (TTree*)fileDB->Get(str.Data());
 191     }
 192     Int_t nBinsKa = fDBgridKa.GetTotBins();
 193     for(Int_t l=0; l<nBinsKa; l++) {
 194       str = "/CovMatrices/Kaons/CovTreeKa_bin";
 195       str+=l;
 196       covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
 197     }
 198     Int_t nBinsPr = fDBgridPr.GetTotBins();
 199     for(Int_t l=0; l<nBinsPr; l++) {
 200       if(fColl==0) {
 201         str = "/CovMatrices/Pions/CovTreePi_bin";
 202       } else {
 203         str = "/CovMatrices/Protons/CovTreePr_bin";
 204       }
 205       str+=l;
 206       covTreePr[l] = (TTree*)fileDB->Get(str.Data());
 207     }
 208     Int_t nBinsEl = fDBgridEl.GetTotBins();
 209     for(Int_t l=0; l<nBinsEl; l++) {
 210       str = "/CovMatrices/Electrons/CovTreeEl_bin";
 211       str+=l;
 212       covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
 213     }
 214     Int_t nBinsMu = fDBgridMu.GetTotBins();
 215     for(Int_t l=0; l<nBinsMu; l++) {
 216       if(fColl==0) {
 217         str = "/CovMatrices/Pions/CovTreePi_bin";
 218       } else {
 219         str = "/CovMatrices/Muons/CovTreeMu_bin";
 220       }
 221       str+=l;
 222       covTreeMu[l] = (TTree*)fileDB->Get(str.Data());
 223     }
 224
 225   } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
 226
 227   TFile *infile=(TFile*)inp;
 228   infile->cd();
 229
 230   // Get gAlice object from file
 231
 232   rl->LoadgAlice();
 233   rl->LoadHeader();
 234   tpcloader->LoadHits("read");
 235
 236   if(!(gAlice=rl->GetAliRun())) {
 237     cerr<<"Can not get gAlice from Run Loader !\n";
 238     return 1;
 239   }
 240
 241   // Check if value in the galice file is equal to selected one (fBz)
 242   AliMagF *fiel = (AliMagF*)gAlice->Field();
 243   Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
 244   printf("Magnetic field is %6.2f Tesla\n",fieval);
 245   if(fBz!=fieval) {
 246     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!"<<endl;
 247     cerr<<"Field selected is: "<<fBz<<" T\n";
 248     cerr<<"Field found on file is: "<<fieval<<" T\n";
 249     return 0;
 250   }
 251
 252   // Get TPC detector
 253   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
 254   Int_t ver = TPC->IsVersion();
 255   cerr<<"+++ TPC version "<<ver<<" has been found !\n";
 256
 257   rl->CdGAFile();
 258   AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
 259   if(digp){
 260     delete digp;
 261     digp = new AliTPCParamSR();
 262   }
 263   else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
 264
 265   if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
 266   TPC->SetParam(digp);
 267
 268   // Set the conversion constant between curvature and Pt
 269   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
 270
 271   TParticle       *Part=0;
 272   AliTPCseedGeant *seed=0;
 273   AliTPCtrack     *tpctrack=0;
 274   Double_t     sPt,sEta;
 275   Int_t        cl=0,bin,label,pdg,charge;
 276   Int_t        tracks;
 277   Int_t        nParticles,nSeeds,arrentr;
 278   Char_t       tname[100];
 279   //Int_t nSel=0,nAcc=0;
 280
 281
 282   // loop over first n events in file
 283   for(Int_t evt=0; evt<fNevents; evt++){
 284     cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
 285     tracks=0;
 286
 287     // tree for TPC tracks
 288     sprintf(tname,"TreeT_TPC_%d",evt);
 289     TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
 290     tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
 291
 292     // array for TPC tracks
 293     TObjArray tArray(20000);
 294
 295     // array for TPC seeds with geant info
 296     TObjArray sArray(20000);
 297
 298     // get the particles stack
 299     nParticles = (Int_t)gAlice->GetEvent(evt);
 300
 301     Bool_t   *done     = new Bool_t[nParticles];
 302     Int_t    *pdgCodes = new Int_t[nParticles];
 303     Double_t *ptkine   = new Double_t[nParticles];
 304     Double_t *pzkine   = new Double_t[nParticles];
 305
 306     // loop on particles and store pdg codes
 307     for(Int_t l=0; l<nParticles; l++) {
 308       Part        = (TParticle*)gAlice->Particle(l);
 309       pdgCodes[l] = Part->GetPdgCode();
 310       ptkine[l]   = Part->Pt();
 311       pzkine[l]   = Part->Pz();
 312       done[l]     = kFALSE;
 313     }
 314     cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nParticles<<
 315       "\n+++\n";
 316
 317     cerr<<"\n ********** MAKING SEEDS *****************"<<endl<<endl;
 318
 319
 320     // Create the seeds for the TPC tracks at the inner radius of TPC
 321     if(fColl==0) {
 322       // Get TreeH with hits
 323       TTree *TH = tpcloader->TreeH();
 324       MakeSeedsFromHits(TPC,TH,sArray);
 325     } else {
 326       // Get TreeTR with track references
 327       TTree *TTR = rl->TreeTR();
 328       MakeSeedsFromRefs(TTR,sArray);
 329     }
 330
 331
 332     nSeeds = sArray.GetEntries();
 333     cerr<<"\n\n+++\n+++ Number of seeds: "<<nSeeds<<"\n+++\n";
 334
 335
 336     cerr<<"\n ********** BUILDING TRACKS **************"<<endl<<endl;
 337
 338     // loop over entries in sArray
 339     for(Int_t l=0; l<nSeeds; l++) {
 340       //if(l%1000==0) cerr<<"  --- Processing seed "
 341       //                <<l<<" of "<<nSeeds<<" ---\r";
 342
 343       seed = (AliTPCseedGeant*)sArray.At(l);
 344
 345       // this is TEMPORARY: only for reconstruction of pp production for charm
 346       if(fColl==1) cl = CheckLabel(seed,nParticles,ptkine,pzkine);
 347       if(cl) continue;
 348
 349       // Get track label
 350       label = seed->GetLabel();
 351
 352       // check if this track has already been processed
 353       if(done[label]) continue;
 354       // PDG code & electric charge
 355       pdg = pdgCodes[label];
 356       if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
 357       else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
 358       else continue;
 359       pdg = TMath::Abs(pdg);
 360       if(pdg>3000) pdg=211;
 361
 362       if(fSelAndSmear) SetParticle(pdg);
 363
 364
 365       sPt  = seed->GetPt();
 366       sEta = seed->GetEta();
 367
 368       // Apply selection according to TPC efficiency
 369       //if(TMath::Abs(pdg)==211) nAcc++;
 370       if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue;
 371       //if(TMath::Abs(pdg)==211) nSel++;
 372
 373       // create AliTPCtrack object
 374       BuildTrack(seed,charge);
 375
 376       if(fSelAndSmear) {
 377         bin = fDBgrid->GetBin(sPt,sEta);
 378         switch (pdg) {
 379         case 211:
 380           fCovTree = covTreePi[bin];
 381           break;
 382         case 321:
 383           fCovTree = covTreeKa[bin];
 384           break;
 385         case 2212:
 386           fCovTree = covTreePr[bin];
 387           break;
 388         case 11:
 389           fCovTree = covTreeEl[bin];
 390           break;
 391         case 13:
 392           fCovTree = covTreeMu[bin];
 393           break;
 394         }
 395         // deal with covariance matrix and smearing of parameters
 396         CookTrack(sPt,sEta);
 397
 398         // assign the track a dE/dx and make a rough PID
 399         CookdEdx(sPt,sEta);
 400       }
 401
 402       // put track in array
 403       AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
 404       iotrack->SetLabel(label);
 405       tArray.AddLast(iotrack);
 406       // Mark track as "done" and register the pdg code
 407       done[label] = kTRUE;
 408       tracks++;
 409
 410     } // loop over entries in sArray
 411
 412
 413     // sort array with TPC tracks (decreasing pT)
 414     tArray.Sort();
 415
 416     arrentr = tArray.GetEntriesFast();
 417     for(Int_t l=0; l<arrentr; l++) {
 418       tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
 419       tracktree->Fill();
 420     }
 421
 422
 423     // write the tree with tracks in the output file
 424     out->cd();
 425     tracktree->Write();
 426
 427     delete tracktree;
 428     delete [] done;
 429     delete [] pdgCodes;
 430     delete [] ptkine;
 431     delete [] pzkine;
 432
 433     printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
 434     //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
 435
 436     sArray.Delete();
 437     tArray.Delete();
 438
 439     infile->cd();
 440   } // loop on events
 441
 442   if(fileDB) fileDB->Close();
 443
 444   return 0;
 445 }
 446 //-----------------------------------------------------------------------------
 447 void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
 448 //-----------------------------------------------------------------------------
 449 // This function computes the dE/dx for pions, kaons, protons and electrons,
 450 // in the [pT,eta] bins.
 451 // Input file is CovMatrix_AllEvts.root.
 452 //-----------------------------------------------------------------------------
 453
 454   gStyle->SetOptStat(0);
 455   gStyle->SetOptFit(10001);
 456
 457   Char_t *part="PIONS";
 458   Double_t ymax=500.;
 459
 460   /*
 461   // create a chain with compared tracks
 462   TChain *cmptrkchain = new ("cmptrktree");
 463   cmptrkchain.Add("CovMatrix_AllEvts.root");
 464   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
 465   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
 466   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
 467   */
 468
 469
 470   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
 471   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
 472
 473   COMPTRACK cmptrk;
 474   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
 475   Int_t entries = (Int_t)cmptrktree->GetEntries();
 476   cerr<<" Number of entries: "<<entries<<endl;
 477
 478   InitializeKineGrid("DB");
 479   InitializeKineGrid("dEdx");
 480
 481   switch(pdg) {
 482   case 211:
 483     part = "PIONS";
 484     ymax = 100.;
 485     break;
 486   case 321:
 487     part = "KAONS";
 488     ymax = 300.;
 489     break;
 490   case 2212:
 491     part = "PROTONS";
 492     ymax = 500.;
 493     break;
 494   case 11:
 495     part = "ELECTRONS";
 496     ymax = 100.;
 497     break;
 498   case 13:
 499     part = "MUONS";
 500     ymax = 100.;
 501     break;
 502   }
 503
 504   SetParticle(pdg);
 505
 506   const Int_t nTotBins = fDBgrid->GetTotBins();
 507
 508   cerr<<" Fit bins: "<<nTotBins<<endl;
 509
 510   Int_t bin=0;
 511   Int_t        *n = new Int_t[nTotBins];
 512   Double_t     *p = new Double_t[nTotBins];
 513   Double_t    *ep = new Double_t[nTotBins];
 514   Double_t  *mean = new Double_t[nTotBins];
 515   Double_t *sigma = new Double_t[nTotBins];
 516
 517   for(Int_t l=0; l<nTotBins; l++) {
 518     n[l] = 1; // set to 1 to avoid divisions by 0
 519     p[l] = mean[l] = sigma[l] = ep[l] = 0.;
 520   }
 521
 522   // loop on chain entries for the mean
 523   for(Int_t l=0; l<entries; l++) {
 524     cmptrktree->GetEvent(l);
 525     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
 526     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
 527     p[bin] += cmptrk.p;
 528     mean[bin] += cmptrk.dEdx;
 529     n[bin]++;
 530   } // loop on chain entries
 531
 532   for(Int_t l=0; l<nTotBins; l++) {
 533     p[l] /= n[l];
 534     mean[l] /= n[l];
 535     n[l] = 1; // set to 1 to avoid divisions by 0
 536   }
 537
 538   // loop on chain entries for the sigma
 539   for(Int_t l=0; l<entries; l++) {
 540     cmptrktree->GetEvent(l);
 541     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
 542     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
 543     if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
 544     n[bin]++;
 545     sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
 546   } // loop on chain entries
 547
 548   for(Int_t l=0; l<nTotBins; l++) {
 549     sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
 550   }
 551
 552   // create the canvas
 553   TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700);
 554
 555   // create the graph for dEdx vs p
 556   TGraphErrors *gr = new TGraphErrors(nTotBins,p,mean,ep,sigma);
 557   TString title("  : dE/dx vs momentum"); title.Prepend(part);
 558   TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
 559   frame->SetXTitle("p [GeV/c]");
 560   frame->SetYTitle("dE/dx [a.u.]");
 561   canv->SetLogx();
 562   frame->Draw();
 563   gr->Draw("P");
 564
 565   switch(pdg) {
 566   case 211:
 567     for(Int_t i=0; i<nTotBins; i++) {
 568       fdEdxMeanPi.SetParam(i,mean[i]);
 569       fdEdxRMSPi.SetParam(i,sigma[i]);
 570     }
 571     break;
 572   case 321:
 573     for(Int_t i=0; i<nTotBins; i++) {
 574       fdEdxMeanKa.SetParam(i,mean[i]);
 575       fdEdxRMSKa.SetParam(i,sigma[i]);
 576     }
 577     break;
 578   case 2212:
 579     for(Int_t i=0; i<nTotBins; i++) {
 580       fdEdxMeanPr.SetParam(i,mean[i]);
 581       fdEdxRMSPr.SetParam(i,sigma[i]);
 582     }
 583     break;
 584   case 11:
 585     for(Int_t i=0; i<nTotBins; i++) {
 586       fdEdxMeanEl.SetParam(i,mean[i]);
 587       fdEdxRMSEl.SetParam(i,sigma[i]);
 588     }
 589     break;
 590   case 13:
 591     for(Int_t i=0; i<nTotBins; i++) {
 592       fdEdxMeanMu.SetParam(i,mean[i]);
 593       fdEdxRMSMu.SetParam(i,sigma[i]);
 594     }
 595     break;
 596   }
 597
 598   // write results to file
 599   WritedEdx(outName,pdg);
 600
 601   delete [] n;
 602   delete [] p;
 603   delete [] ep;
 604   delete [] mean;
 605   delete [] sigma;
 606
 607   return;
 608 }
 609 //-----------------------------------------------------------------------------
 610 void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
 611 //-----------------------------------------------------------------------------
 612 // This function computes the pulls for pions, kaons and electrons,
 613 // in the [pT,eta] bins.
 614 // Input file is CovMatrix_AllEvts.root.
 615 // Output file is pulls.root.
 616 //-----------------------------------------------------------------------------
 617
 618   /*
 619   // create a chain with compared tracks
 620   TChain *cmptrkchain = new ("cmptrktree");
 621   cmptrkchain.Add("CovMatrix_AllEvts.root");
 622   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
 623   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
 624   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
 625   */
 626
 627
 628   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
 629   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
 630
 631   COMPTRACK cmptrk;
 632   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
 633   Int_t entries = (Int_t)cmptrktree->GetEntries();
 634   cerr<<" Number of entries: "<<entries<<endl;
 635
 636   Int_t thisPdg=0;
 637   Char_t hname[100], htitle[100];
 638   Int_t nTotBins,bin;
 639
 640   AliTPCkineGrid  pulls[5];
 641   TH1F *hDum = new TH1F("name","title",100,-7.,7.);
 642   TF1 *g = new TF1("g","gaus");
 643
 644   InitializeKineGrid("pulls");
 645   InitializeKineGrid("DB");
 646
 647
 648
 649   // loop on the particles Pi,Ka,Pr,El,Mu
 650   for(Int_t part=0; part<5; part++) {
 651
 652     switch (part) {
 653     case 0: // pions
 654       thisPdg=211;
 655       cerr<<" Processing pions ...\n";
 656       break;
 657     case 1: // kaons
 658       thisPdg=321;
 659       cerr<<" Processing kaons ...\n";
 660       break;
 661     case 2: // protons
 662       thisPdg=2212;
 663       cerr<<" Processing protons ...\n";
 664       break;
 665     case 3: // electrons
 666       thisPdg=11;
 667       cerr<<" Processing electrons ...\n";
 668       break;
 669     case 4: // muons
 670       thisPdg=13;
 671       cerr<<" Processing muons ...\n";
 672       break;
 673     }
 674
 675     SetParticle(thisPdg);
 676
 677     for(Int_t i=0;i<5;i++) {
 678       pulls[i].~AliTPCkineGrid();
 679       new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
 680     }
 681     nTotBins = fDBgrid->GetTotBins();
 682     cerr<<"nTotBins = "<<nTotBins<<endl;
 683
 684     // create histograms for the all the bins
 685     TH1F *hPulls0_=NULL;
 686     TH1F *hPulls1_=NULL;
 687     TH1F *hPulls2_=NULL;
 688     TH1F *hPulls3_=NULL;
 689     TH1F *hPulls4_=NULL;
 690
 691     hPulls0_ = new TH1F[nTotBins];
 692     hPulls1_ = new TH1F[nTotBins];
 693     hPulls2_ = new TH1F[nTotBins];
 694     hPulls3_ = new TH1F[nTotBins];
 695     hPulls4_ = new TH1F[nTotBins];
 696
 697
 698     for(Int_t i=0; i<nTotBins; i++) {
 699       sprintf(hname,"hPulls0_%d",i);
 700       sprintf(htitle,"P0 pulls for bin %d",i);
 701       hDum->SetName(hname); hDum->SetTitle(htitle);
 702       hPulls0_[i] = *hDum;
 703       sprintf(hname,"hPulls1_%d",i);
 704       sprintf(htitle,"P1 pulls for bin %d",i);
 705       hDum->SetName(hname); hDum->SetTitle(htitle);
 706       hPulls1_[i] = *hDum;
 707       sprintf(hname,"hPulls2_%d",i);
 708       sprintf(htitle,"P2 pulls for bin %d",i);
 709       hDum->SetName(hname); hDum->SetTitle(htitle);
 710       hPulls2_[i] = *hDum;
 711       sprintf(hname,"hPulls3_%d",i);
 712       sprintf(htitle,"P3 pulls for bin %d",i);
 713       hDum->SetName(hname); hDum->SetTitle(htitle);
 714       hPulls3_[i] = *hDum;
 715       sprintf(hname,"hPulls4_%d",i);
 716       sprintf(htitle,"P4 pulls for bin %d",i);
 717       hDum->SetName(hname); hDum->SetTitle(htitle);
 718       hPulls4_[i] = *hDum;
 719     }
 720
 721     // loop on chain entries
 722     for(Int_t i=0; i<entries; i++) {
 723       cmptrktree->GetEvent(i);
 724       if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
 725       // fill histograms with the pulls
 726       bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
 727       //cerr<<" pt "<<cmptrk.pt<<"   eta "<<cmptrk.eta<<"   bin "<<bin<<endl;
 728       hPulls0_[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
 729       hPulls1_[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
 730       hPulls2_[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
 731       hPulls3_[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
 732       hPulls4_[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
 733     } // loop on chain entries
 734
 735     // compute the sigma of the distributions
 736     for(Int_t i=0; i<nTotBins; i++) {
 737       if(hPulls0_[i].GetEntries()>10) {
 738         g->SetRange(-3.*hPulls0_[i].GetRMS(),3.*hPulls0_[i].GetRMS());
 739         hPulls0_[i].Fit("g","R,Q,N");
 740         pulls[0].SetParam(i,g->GetParameter(2));
 741       } else pulls[0].SetParam(i,-1.);
 742       if(hPulls1_[i].GetEntries()>10) {
 743         g->SetRange(-3.*hPulls1_[i].GetRMS(),3.*hPulls1_[i].GetRMS());
 744         hPulls1_[i].Fit("g","R,Q,N");
 745         pulls[1].SetParam(i,g->GetParameter(2));
 746       } else pulls[1].SetParam(i,-1.);
 747       if(hPulls2_[i].GetEntries()>10) {
 748         g->SetRange(-3.*hPulls2_[i].GetRMS(),3.*hPulls2_[i].GetRMS());
 749         hPulls2_[i].Fit("g","R,Q,N");
 750         pulls[2].SetParam(i,g->GetParameter(2));
 751       } else pulls[2].SetParam(i,-1.);
 752       if(hPulls3_[i].GetEntries()>10) {
 753         g->SetRange(-3.*hPulls3_[i].GetRMS(),3.*hPulls3_[i].GetRMS());
 754         hPulls3_[i].Fit("g","R,Q,N");
 755         pulls[3].SetParam(i,g->GetParameter(2));
 756       } else pulls[3].SetParam(i,-1.);
 757       if(hPulls4_[i].GetEntries()>10) {
 758         g->SetRange(-3.*hPulls4_[i].GetRMS(),3.*hPulls4_[i].GetRMS());
 759         hPulls4_[i].Fit("g","R,Q,N");
 760         pulls[4].SetParam(i,g->GetParameter(2));
 761       } else pulls[4].SetParam(i,-1.);
 762     } // loop on bins
 763
 764
 765     switch (part) {
 766     case 0: // pions
 767       for(Int_t i=0;i<5;i++) {
 768         fPullsPi[i].~AliTPCkineGrid();
 769         new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
 770       }
 771       break;
 772     case 1: // kaons
 773       for(Int_t i=0;i<5;i++) {
 774         fPullsKa[i].~AliTPCkineGrid();
 775         new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
 776       }
 777       break;
 778     case 2: // protons
 779       for(Int_t i=0;i<5;i++) {
 780         fPullsPr[i].~AliTPCkineGrid();
 781         new(&fPullsPr[i]) AliTPCkineGrid(pulls[i]);
 782       }
 783       break;
 784     case 3: // electrons
 785       for(Int_t i=0;i<5;i++) {
 786         fPullsEl[i].~AliTPCkineGrid();
 787         new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
 788       }
 789       break;
 790     case 4: // muons
 791       for(Int_t i=0;i<5;i++) {
 792         fPullsMu[i].~AliTPCkineGrid();
 793         new(&fPullsMu[i]) AliTPCkineGrid(pulls[i]);
 794         //cerr<<" mu  pulls "<<i<<"  "<<fPullsMu[i].GetParam(0)<<endl;
 795       }
 796       break;
 797     }
 798
 799     delete [] hPulls0_;
 800     delete [] hPulls1_;
 801     delete [] hPulls2_;
 802     delete [] hPulls3_;
 803     delete [] hPulls4_;
 804
 805   } // loop on particle species
 806
 807   // write pulls to file
 808   WritePulls(outName);
 809
 810
 811   return;
 812 }
 813 //-----------------------------------------------------------------------------
 814 void AliTPCtrackerParam::AnalyzeResolutions(Int_t pdg) {
 815 //-----------------------------------------------------------------------------
 816 // This function computes the resolutions:
 817 // - dy
 818 // - dC
 819 // - dPt/Pt
 820 // as a function of Pt
 821 // Input file is CovMatrix_AllEvts.root.
 822 //-----------------------------------------------------------------------------
 823
 824   /*
 825   // create a chain with compared tracks
 826   TChain *cmptrkchain = new ("cmptrktree");
 827   cmptrkchain.Add("CovMatrix_AllEvts.root");
 828   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
 829   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
 830   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
 831   */
 832
 833
 834   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
 835   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
 836
 837   COMPTRACK cmptrk;
 838   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
 839   Int_t entries = (Int_t)cmptrktree->GetEntries();
 840   cerr<<" Number of entries: "<<entries<<endl;
 841
 842
 843   Int_t bin = 0;
 844
 845   InitializeKineGrid("DB");
 846   InitializeKineGrid("eff");
 847
 848   SetParticle(pdg);
 849
 850   const Int_t nPtBins = fEff->GetPointsPt();
 851   cerr<<"nPtBins = "<<nPtBins<<endl;
 852   Double_t *dP0     = new Double_t[nPtBins];
 853   Double_t *dP4     = new Double_t[nPtBins];
 854   Double_t *dPtToPt = new Double_t[nPtBins];
 855   Double_t *pt      = new Double_t[nPtBins];
 856   fEff->GetArrayPt(pt);
 857
 858
 859   TH1F *hDumP0 = new TH1F("nameP0","dy",100,-0.3,0.3);
 860   TH1F *hDumP4 = new TH1F("nameP4","dC",100,-0.0005,0.0005);
 861   TH1F *hDumPt = new TH1F("namePt","dp_{T}/p_{T}",100,-0.5,0.5);
 862
 863   TF1 *g = new TF1("g","gaus");
 864
 865   // create histograms for the all the bins
 866   TH1F *hP0_=NULL;
 867   TH1F *hP4_=NULL;
 868   TH1F *hPt_=NULL;
 869
 870   hP0_ = new TH1F[nPtBins];
 871   hP4_ = new TH1F[nPtBins];
 872   hPt_ = new TH1F[nPtBins];
 873
 874   for(Int_t i=0; i<nPtBins; i++) {
 875     hP0_[i] = *hDumP0;
 876     hP4_[i] = *hDumP4;
 877     hPt_[i] = *hDumPt;
 878   }
 879
 880   // loop on chain entries
 881   for(Int_t i=0; i<entries; i++) {
 882     cmptrktree->GetEvent(i);
 883     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
 884     // fill histograms with the residuals
 885     bin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
 886     //cerr<<" pt "<<cmptrk.pt<<"   eta "<<cmptrk.eta<<"   bin "<<bin<<endl;
 887     hP0_[bin].Fill(cmptrk.dP0);
 888     hP4_[bin].Fill(cmptrk.dP4);
 889     hPt_[bin].Fill(cmptrk.dpt/cmptrk.pt);
 890   } // loop on chain entries
 891
 892
 893   TCanvas *cP0res = new TCanvas("cP0res","cP0res",0,0,1200,700);
 894   cP0res->Divide(5,2);
 895   TCanvas *cP4res = new TCanvas("cP4res","cP4res",0,0,1200,700);
 896   cP4res->Divide(5,2);
 897   TCanvas *cPtres = new TCanvas("cPtres","cPtres",0,0,1200,700);
 898   cPtres->Divide(5,2);
 899
 900   // Draw histograms
 901   for(Int_t i=0; i<nPtBins; i++) {
 902     cP0res->cd(i+1); hP0_[i].Draw();
 903     cP4res->cd(i+1); hP4_[i].Draw();
 904     cPtres->cd(i+1); hPt_[i].Draw();
 905   }
 906
 907
 908   // compute the sigma of the distributions
 909   for(Int_t i=0; i<nPtBins; i++) {
 910     if(hP0_[i].GetEntries()>10) {
 911       g->SetRange(-3.*hP0_[i].GetRMS(),3.*hP0_[i].GetRMS());
 912       hP0_[i].Fit("g","R,Q,N");
 913       dP0[i] = g->GetParameter(2);
 914     } else dP0[i] = 0.;
 915     if(hP4_[i].GetEntries()>10) {
 916       g->SetRange(-3.*hP4_[i].GetRMS(),3.*hP4_[i].GetRMS());
 917       hP4_[i].Fit("g","R,Q,N");
 918       dP4[i] = g->GetParameter(2);
 919     } else dP4[i] = 0.;
 920     if(hPt_[i].GetEntries()>10) {
 921       g->SetRange(-3.*hPt_[i].GetRMS(),3.*hPt_[i].GetRMS());
 922       hPt_[i].Fit("g","R,Q,N");
 923       dPtToPt[i] = 100.*g->GetParameter(2);
 924     } else dPtToPt[i] = 0.;
 925   } // loop on bins
 926
 927
 928   TGraph *grdP0 = new TGraph(nPtBins,pt,dP0);
 929   TGraph *grdP4 = new TGraph(nPtBins,pt,dP4);
 930   TGraph *grdPtToPt = new TGraph(nPtBins,pt,dPtToPt);
 931
 932   grdP0->SetMarkerStyle(20); grdP0->SetMarkerColor(2); grdP0->SetMarkerSize(1.5);
 933   grdP4->SetMarkerStyle(21); grdP4->SetMarkerColor(3); grdP4->SetMarkerSize(1.5);
 934   grdPtToPt->SetMarkerStyle(22); grdPtToPt->SetMarkerColor(4); grdPtToPt->SetMarkerSize(1.5);
 935
 936   // Plot Results
 937   gStyle->SetOptStat(0);
 938   TCanvas *c1 = new TCanvas("c1","dP0",0,0,900,700);
 939   c1->SetLogx();
 940   c1->SetGridx();
 941   c1->SetGridy();
 942
 943   TH2F *frame1 = new TH2F("frame1","y resolution VS p_{T} in TPC",2,0.1,30,2,0,0.1);
 944   frame1->SetXTitle("p_{T} [GeV/c]");
 945   frame1->SetYTitle("#sigma(y) [cm]");
 946   frame1->Draw();
 947   grdP0->Draw("P");
 948
 949
 950   TCanvas *c2 = new TCanvas("c2","dP4",0,0,900,700);
 951   c2->SetLogx();
 952   c2->SetGridx();
 953   c2->SetGridy();
 954
 955   TH2F *frame2 = new TH2F("frame2","C resolution VS p_{T} in TPC",2,0.1,30,2,0,0.0001);
 956   frame2->SetXTitle("p_{T} [GeV/c]");
 957   frame2->SetYTitle("#sigma(C) [1/cm]");
 958   frame2->Draw();
 959   grdP4->Draw("P");
 960
 961   TCanvas *c3 = new TCanvas("c3","dPtToPt",0,0,900,700);
 962   c3->SetLogx();
 963   c3->SetLogy();
 964   c3->SetGridx();
 965   c3->SetGridy();
 966
 967   TH2F *frame3 = new TH2F("frame3","Relative p_{T} resolution VS p_{T} in TPC",2,0.1,30,2,0.1,30.);
 968   frame3->SetXTitle("p_{T} [GeV/c]");
 969   frame3->SetYTitle("dp_{T}/p_{T} (%)");
 970   frame3->Draw();
 971   grdPtToPt->Draw("P");
 972
 973
 974   delete [] dP0;
 975   delete [] dP4;
 976   delete [] dPtToPt;
 977   delete [] pt;
 978
 979
 980   delete [] hP0_;
 981   delete [] hP4_;
 982   delete [] hPt_;
 983
 984   return;
 985 }
 986 //-----------------------------------------------------------------------------
 987 void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {
 988 //-----------------------------------------------------------------------------
 989 // This function uses GEANT info to set true track parameters
 990 //-----------------------------------------------------------------------------
 991   Double_t xref = s->GetXL();
 992   Double_t xx[5],cc[15];
 993   cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
 994   cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
 995
 996   // Magnetic field
 997   TVector3 bfield(0.,0.,fBz);
 998
 999
1000   // radius [cm] of track projection in (x,y)
1001   Double_t rho = s->GetPt()*100./0.299792458/bfield.Z();
1002   // center of track projection in local reference frame
1003   TVector3 sMom,sPos;
1004
1005
1006   // position (local) and momentum (local) at the seed
1007   // in the bending plane (z=0)
1008   sPos.SetXYZ(s->GetXL(),s->GetYL(),0.);
1009   sMom.SetXYZ(s->GetPx()*TMath::Cos(s->GetAlpha())+s->GetPy()*TMath::Sin(s->GetAlpha()),-s->GetPx()*TMath::Sin(s->GetAlpha())+s->GetPy()*TMath::Cos(s->GetAlpha()),0.);
1010   TVector3 vrho = sMom.Cross(bfield);
1011   vrho *= ch;
1012   vrho.SetMag(rho);
1013
1014   TVector3 vcenter = sPos+vrho;
1015
1016   Double_t x0 = vcenter.X();
1017
1018   // fX     = xref         X-coordinate of this track (reference plane)
1019   // fAlpha = Alpha        Rotation angle the local (TPC sector)
1020   // fP0    = YL           Y-coordinate of a track
1021   // fP1    = ZG           Z-coordinate of a track
1022   // fP2    = C*x0         x0 is center x in rotated frame
1023   // fP3    = Tgl          tangent of the track momentum dip angle
1024   // fP4    = C            track curvature
1025   xx[0] = s->GetYL();
1026   xx[1] = s->GetZL();
1027   xx[3] = s->GetPz()/s->GetPt();
1028   xx[4] = -ch/rho;
1029   xx[2] = xx[4]*x0;
1030
1031   // create the object AliTPCtrack
1032   AliTPCtrack track(0,xx,cc,xref,s->GetAlpha());
1033   new(&fTrack) AliTPCtrack(track);
1034
1035   return;
1036 }
1037 //-----------------------------------------------------------------------------
1038 Int_t AliTPCtrackerParam::CheckLabel(AliTPCseedGeant *s,Int_t nPart,
1039                                      Double_t *ptkine,Double_t *pzkine) const {
1040 //-----------------------------------------------------------------------------
1041 // This function checks if the label of the seed has been correctly
1042 // assigned (to do only for pp charm production with AliRoot v3-08-02)
1043 //-----------------------------------------------------------------------------
1044
1045   Int_t sLabel = s->GetLabel();
1046   Double_t sPt = s->GetPt();
1047   Double_t sPz = s->GetPz();
1048
1049   // check if the label is correct (comparing momentum)
1050   if(sLabel<nPart &&
1051      TMath::Abs(sPt-ptkine[sLabel])*
1052      TMath::Abs(sPz-pzkine[sLabel])<0.001) return 0;
1053
1054   if((sLabel-30)>=nPart) return 1;
1055
1056   Double_t diff=0,mindiff=1000.;
1057   Int_t bestLabel=0;
1058
1059   for(Int_t i=sLabel-30; i<sLabel; i++) {
1060     if(i<0 || i>=nPart) continue;
1061     diff = TMath::Abs(sPt-ptkine[i])*TMath::Abs(sPz-pzkine[i]);
1062     if(diff<mindiff) { mindiff = diff; bestLabel = i; }
1063   }
1064
1065   if(mindiff>0.001) return 1;
1066   s->SetLabel(bestLabel);
1067
1068   return 0;
1069 }
1070 //-----------------------------------------------------------------------------
1071 void AliTPCtrackerParam::CompareTPCtracks(
1072                            const Char_t* galiceName,
1073                            const Char_t* trkGeaName,
1074                            const Char_t* trkKalName,
1075                            const Char_t* covmatName,
1076                            const Char_t* tpceffasciiName,
1077                            const Char_t* tpceffrootName) {
1078 //-----------------------------------------------------------------------------
1079 // This function compares tracks from TPC Kalman Filter V2 with
1080 // geant tracks at TPC 1st hit. It gives:
1081 //   - a tree with Kalman cov. matrix elements, resolutions, dEdx
1082 //   - the efficiencies as a function of particle type, pT, eta
1083 //-----------------------------------------------------------------------------
1084
1085   TFile *kalFile    = TFile::Open(trkKalName);
1086   TFile *geaFile    = TFile::Open(trkGeaName);
1087   TFile *galiceFile = TFile::Open(galiceName);
1088
1089   // get the AliRun object
1090   AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
1091
1092
1093   // create the tree for comparison results
1094   COMPTRACK cmptrk;
1095   TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
1096   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");
1097
1098   InitializeKineGrid("eff");
1099   InitializeKineGrid("DB");
1100   Int_t   effBins = fEffPi.GetTotPoints();
1101   Int_t effBinsPt = fEffPi.GetPointsPt();
1102   Double_t    *pt = new Double_t[effBinsPt];
1103   fEffPi.GetArrayPt(pt);
1104
1105   TParticle *Part;
1106   Double_t ptgener;
1107   Bool_t   usethis;
1108   Int_t    label;
1109   Double_t cc[15],dAlpha;
1110   Int_t    pi=0,ka=0,mu=0,el=0,pr=0;
1111   Int_t   *geaPi = new Int_t[effBins];
1112   Int_t   *geaKa = new Int_t[effBins];
1113   Int_t   *geaPr = new Int_t[effBins];
1114   Int_t   *geaEl = new Int_t[effBins];
1115   Int_t   *geaMu = new Int_t[effBins];
1116   Int_t   *kalPi = new Int_t[effBins];
1117   Int_t   *kalKa = new Int_t[effBins];
1118   Int_t   *kalPr = new Int_t[effBins];
1119   Int_t   *kalEl = new Int_t[effBins];
1120   Int_t   *kalMu = new Int_t[effBins];
1121   Float_t *effPi = new Float_t[effBins];
1122   Float_t *effKa = new Float_t[effBins];
1123   Float_t *effPr = new Float_t[effBins];
1124   Float_t *effEl = new Float_t[effBins];
1125   Float_t *effMu = new Float_t[effBins];
1126   Int_t bin;
1127   for(Int_t j=0; j<effBins; j++) {
1128     geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
1129     kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
1130     effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
1131   }
1132
1133   Char_t tname[100];
1134
1135   // loop on events in file
1136   for(Int_t evt=0; evt<fNevents; evt++) {
1137     cerr<<"\n  --- Reading tracks for event "<<evt<<" ---\n\n";
1138     sprintf(tname,"TreeT_TPC_%d",evt);
1139
1140     // particles from TreeK
1141     const Int_t nparticles = gAlice->GetEvent(evt);
1142
1143     Int_t *kalLab = new Int_t[nparticles];
1144     for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1;
1145
1146
1147     // tracks from Kalman
1148     TTree *kaltree=(TTree*)kalFile->Get(tname);
1149     if(!kaltree) continue;
1150     AliTPCtrack *kaltrack=new AliTPCtrack;
1151     kaltree->SetBranchAddress("tracks",&kaltrack);
1152     Int_t kalEntries = (Int_t)kaltree->GetEntries();
1153
1154     // tracks from 1st hit
1155     TTree *geatree=(TTree*)geaFile->Get(tname);
1156     if(!geatree) continue;
1157     AliTPCtrack *geatrack=new AliTPCtrack;
1158     geatree->SetBranchAddress("tracks",&geatrack);
1159     Int_t geaEntries = (Int_t)geatree->GetEntries();
1160
1161     cerr<<"+++\n+++ Number of tracks:  TPC Kalman  = "<<kalEntries<<endl<<"+++                    TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
1162
1163     // set pointers for TPC tracks:
1164     // the entry number of the track labelled l is stored in kalLab[l]
1165     Int_t fake=0,mult=0;
1166     for (Int_t j=0; j<kalEntries; j++) {
1167       kaltree->GetEvent(j);
1168       if(kaltrack->GetLabel()>=0) {
1169         if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
1170         kalLab[kaltrack->GetLabel()] = j;
1171       }
1172       else fake++;
1173     }
1174
1175     cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
1176     cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
1177
1178     /*
1179     // Read the labels of the seeds
1180     char sname[100];
1181     Int_t sLabel,ncol;
1182     Bool_t *hasSeed = new Bool_t[nparticles];
1183     for(Int_t i=0; i<nparticles; i++) hasSeed[i] = kFALSE;
1184     sprintf(sname,"seedLabels.%d.dat",evt);
1185     FILE *seedFile = fopen(sname,"r");
1186     while(1) {
1187       ncol = fscanf(seedFile,"%d",&sLabel);
1188       if(ncol<1) break;
1189       if(sLabel>0) hasSeed[sLabel]=kTRUE;
1190     }
1191     fclose(seedFile);
1192     */
1193
1194     cerr<<"Doing track comparison...\n";
1195     // loop on tracks at 1st hit
1196     for(Int_t j=0; j<geaEntries; j++) {
1197       geatree->GetEvent(j);
1198
1199       label = geatrack->GetLabel();
1200       Part = (TParticle*)gAlice->Particle(label);
1201
1202       // use only injected tracks with fixed values of pT
1203       ptgener = Part->Pt();
1204       usethis = kFALSE;
1205       for(Int_t l=0; l<fEffPi.GetPointsPt(); l++) {
1206         if(TMath::Abs(ptgener-pt[l])<0.01) usethis = kTRUE;
1207       }
1208       if(!usethis) continue;
1209
1210       // check if it has the seed
1211       //if(!hasSeed[label]) continue;
1212
1213       /*
1214       // check if track is entirely contained in a TPC sector
1215       Bool_t out = kFALSE;
1216       for(Int_t l=0; l<10; l++) {
1217         Double_t x = 85. + (250.-85.)*(Double_t)l/9.;
1218         //cerr<<"x "<<x<<"    X "<<geatrack->GetX()<<endl;
1219         Double_t y = geatrack->GetY() + (
1220         TMath::Sqrt(1-(geatrack->GetC()*geatrack->GetX()-geatrack->GetEta())*
1221                       (geatrack->GetC()*geatrack->GetX()-geatrack->GetEta()))
1222        -TMath::Sqrt(1-(geatrack->GetC()*x-geatrack->GetEta())*
1223                       (geatrack->GetC()*x-geatrack->GetEta()))
1224         )/geatrack->GetC();
1225         y = TMath::Abs(y);
1226         //cerr<<"y "<<y<<"    Y "<<geatrack->GetY()<<endl;
1227         if(y > 0.8*x*TMath::Tan(TMath::Pi()/18)) { out = kTRUE; break; }
1228       }
1229       if(out) continue;
1230       */
1231
1232       cmptrk.pdg = Part->GetPdgCode();
1233       cmptrk.eta = Part->Eta();
1234       cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
1235
1236       cmptrk.pt   = 1/TMath::Abs(geatrack->Get1Pt());
1237       cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
1238       cmptrk.p    = cmptrk.pt/cmptrk.cosl;
1239
1240
1241       bin = fDBgridPi.GetBin(cmptrk.pt,cmptrk.eta);
1242       cmptrk.bin = bin;
1243       if(abs(cmptrk.pdg)==211)  geaPi[bin]++;
1244       if(abs(cmptrk.pdg)==321)  geaKa[bin]++;
1245       if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
1246       if(abs(cmptrk.pdg)==11)   geaEl[bin]++;
1247       if(abs(cmptrk.pdg)==13)   geaMu[bin]++;
1248
1249
1250       // check if there is the corresponding track in the TPC kalman and get it
1251       if(kalLab[label]==-1) continue;
1252       kaltree->GetEvent(kalLab[label]);
1253
1254       // and go on only if it has xref = 84.57 cm (inner pad row)
1255       if(kaltrack->GetX()>90.) continue;
1256
1257       if(abs(cmptrk.pdg)==211)  { kalPi[bin]++; pi++; }
1258       if(abs(cmptrk.pdg)==321)  { kalKa[bin]++; ka++; }
1259       if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
1260       if(abs(cmptrk.pdg)==11)   { kalEl[bin]++; el++; }
1261       if(abs(cmptrk.pdg)==13)   { kalMu[bin]++; mu++; }
1262
1263       kaltrack->PropagateTo(geatrack->GetX(),1.e-9,0.);
1264
1265       cmptrk.dEdx = kaltrack->GetdEdx();
1266
1267       // compute errors on parameters
1268       dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
1269       if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
1270
1271       cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
1272       cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
1273       cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
1274       cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
1275       cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
1276       cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
1277
1278       // get covariance matrix
1279       // beware: lines 3 and 4 in the matrix are inverted!
1280       kaltrack->GetCovariance(cc);
1281
1282       cmptrk.c00 = cc[0];
1283       cmptrk.c10 = cc[1];
1284       cmptrk.c11 = cc[2];
1285       cmptrk.c20 = cc[3];
1286       cmptrk.c21 = cc[4];
1287       cmptrk.c22 = cc[5];
1288       cmptrk.c30 = cc[10];
1289       cmptrk.c31 = cc[11];
1290       cmptrk.c32 = cc[12];
1291       cmptrk.c33 = cc[14];
1292       cmptrk.c40 = cc[6];
1293       cmptrk.c41 = cc[7];
1294       cmptrk.c42 = cc[8];
1295       cmptrk.c43 = cc[13];
1296       cmptrk.c44 = cc[9];
1297
1298       // fill tree
1299       cmptrktree->Fill();
1300
1301     } // loop on tracks at TPC 1st hit
1302
1303     delete [] kalLab;
1304     //delete [] hasSeed;
1305
1306   } // end loop on events in file
1307
1308
1309   cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
1310   cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
1311   cerr<<"+++"<<endl;
1312
1313   // Write tree to file
1314   TFile *outfile = new TFile(covmatName,"recreate");
1315   cmptrktree->Write();
1316   outfile->Close();
1317
1318
1319   // Write efficiencies to ascii file
1320   FILE *effFile = fopen(tpceffasciiName,"w");
1321   //fprintf(effFile,"%d\n",kalEntries);
1322   for(Int_t j=0; j<effBins; j++) {
1323     if(geaPi[j]>=100) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
1324     if(geaKa[j]>=100) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
1325     if(geaPr[j]>=100) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
1326     if(geaEl[j]>=500) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
1327     if(geaMu[j]>=100) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
1328     fprintf(effFile,"%f  %f  %f  %f  %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
1329   }
1330
1331   for(Int_t j=0; j<effBins; j++) {
1332     fprintf(effFile,"%d  %d  %d  %d  %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
1333   }
1334   for(Int_t j=0; j<effBins; j++) {
1335     fprintf(effFile,"%d  %d  %d  %d  %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
1336   }
1337   fclose(effFile);
1338
1339   // Write efficiencies to root file
1340   for(Int_t j=0; j<effBins; j++) {
1341     fEffPi.SetParam(j,(Double_t)effPi[j]);
1342     fEffKa.SetParam(j,(Double_t)effKa[j]);
1343     fEffPr.SetParam(j,(Double_t)effPr[j]);
1344     fEffEl.SetParam(j,(Double_t)effEl[j]);
1345     fEffMu.SetParam(j,(Double_t)effMu[j]);
1346   }
1347   WriteEffs(tpceffrootName);
1348
1349   // delete AliRun object
1350   delete gAlice; gAlice=0;
1351
1352   // close all input files
1353   kalFile->Close();
1354   geaFile->Close();
1355   galiceFile->Close();
1356
1357   delete [] pt;
1358   delete [] effPi;
1359   delete [] effKa;
1360   delete [] effPr;
1361   delete [] effEl;
1362   delete [] effMu;
1363   delete [] geaPi;
1364   delete [] geaKa;
1365   delete [] geaPr;
1366   delete [] geaEl;
1367   delete [] geaMu;
1368   delete [] kalPi;
1369   delete [] kalKa;
1370   delete [] kalPr;
1371   delete [] kalEl;
1372   delete [] kalMu;
1373
1374   return;
1375 }
1376 //-----------------------------------------------------------------------------
1377 void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
1378 //-----------------------------------------------------------------------------
1379 // This function assigns the track a dE/dx and makes a rough PID
1380 //-----------------------------------------------------------------------------
1381
1382   Double_t mean = fdEdxMean->GetValueAt(pt,eta);
1383   Double_t rms  = fdEdxRMS->GetValueAt(pt,eta);
1384
1385   Double_t dEdx = gRandom->Gaus(mean,rms);
1386
1387   fTrack.SetdEdx(dEdx);
1388
1389   AliTPCtrackParam t(fTrack);
1390
1391   //Very rough PID
1392   Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
1393
1394   if (p<0.6) {
1395     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
1396       t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1397     }
1398     if (dEdx < 39.+ 12./p/p) {
1399       t.AssignMass(0.49368); new(&fTrack) AliTPCtrack(t); return;
1400     }
1401     t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
1402   }
1403
1404   if (p<1.2) {
1405     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
1406       t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1407     }
1408     t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
1409   }
1410
1411   t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1412 }
1413 //-----------------------------------------------------------------------------
1414 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
1415 //-----------------------------------------------------------------------------
1416 // This function deals with covariance matrix and smearing
1417 //-----------------------------------------------------------------------------
1418
1419   COVMATRIX covmat;
1420   Double_t  p,cosl;
1421   Double_t  trkKine[1],trkRegPar[3];
1422   Double_t  xref,alpha,xx[5],xxsm[5],cc[15];
1423   Int_t     treeEntries;
1424
1425   fCovTree->SetBranchAddress("matrix",&covmat);
1426
1427   // get random entry from the tree
1428   treeEntries = (Int_t)fCovTree->GetEntries();
1429   fCovTree->GetEvent(gRandom->Integer(treeEntries));
1430
1431   // get P and Cosl from track
1432   cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1433   p    = 1./TMath::Abs(fTrack.Get1Pt())/cosl;
1434
1435   trkKine[0] = p;
1436
1437   // get covariance matrix from regularized matrix
1438   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(0,l);
1439   cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1440   cc[1] = covmat.c10;
1441   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(1,l);
1442   cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1443   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(2,l);
1444   cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1445   cc[4] = covmat.c21;
1446   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(3,l);
1447   cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1448   cc[6] = covmat.c30;
1449   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(4,l);
1450   cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1451   cc[8] = covmat.c32;
1452   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(5,l);
1453   cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1454   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(6,l);
1455   cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1456   cc[11]= covmat.c41;
1457   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(7,l);
1458   cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1459   cc[13]= covmat.c43;
1460   for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l);
1461   cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1462
1463   TMatrixD covMatSmear(5,5);
1464
1465   covMatSmear = GetSmearingMatrix(cc,pt,eta);
1466
1467   // get track original parameters
1468   xref =fTrack.GetX();
1469   alpha=fTrack.GetAlpha();
1470   xx[0]=fTrack.GetY();
1471   xx[1]=fTrack.GetZ();
1472   xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp();
1473   xx[3]=fTrack.GetTgl();
1474   xx[4]=fTrack.GetC();
1475
1476   // use smearing matrix to smear the original parameters
1477   SmearTrack(xx,xxsm,covMatSmear);
1478
1479   AliTPCtrack track(0,xxsm,cc,xref,alpha);
1480   new(&fTrack) AliTPCtrack(track);
1481
1482   return;
1483 }
1484 //-----------------------------------------------------------------------------
1485 void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg) {
1486 //-----------------------------------------------------------------------------
1487 // This function draws the TPC efficiencies in the [pT,eta] bins
1488 //-----------------------------------------------------------------------------
1489
1490   ReadEffs(inName);
1491   SetParticle(pdg);
1492
1493   const Int_t n = fEff->GetPointsPt();
1494   Double_t *effsA = new Double_t[n];
1495   Double_t *effsB = new Double_t[n];
1496   Double_t *effsC = new Double_t[n];
1497   Double_t *pt    = new Double_t[n];
1498
1499   fEff->GetArrayPt(pt);
1500   for(Int_t i=0;i<n;i++) {
1501     effsA[i] = fEff->GetParam(i,0);
1502     effsB[i] = fEff->GetParam(i,1);
1503     effsC[i] = fEff->GetParam(i,2);
1504   }
1505
1506   TGraph *grA = new TGraph(n,pt,effsA);
1507   TGraph *grB = new TGraph(n,pt,effsB);
1508   TGraph *grC = new TGraph(n,pt,effsC);
1509
1510   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1511   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1512   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1513
1514   TString title("Distribution of the TPC efficiencies");
1515   switch (pdg) {
1516   case 211:
1517     title.Prepend("PIONS - ");
1518     break;
1519   case 321:
1520     title.Prepend("KAONS - ");
1521     break;
1522   case 2212:
1523     title.Prepend("PROTONS - ");
1524     break;
1525   case 11:
1526     title.Prepend("ELECTRONS - ");
1527     break;
1528   case 13:
1529     title.Prepend("MUONS - ");
1530     break;
1531   }
1532
1533
1534   gStyle->SetOptStat(0);
1535   TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1536   c->SetLogx();
1537   c->SetGridx();
1538   c->SetGridy();
1539
1540   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1541   frame->SetXTitle("p_{T} [GeV/c]");
1542   frame->Draw();
1543   grA->Draw("P");
1544   grB->Draw("P");
1545   grC->Draw("P");
1546
1547   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1548   leg->AddEntry(grA,"|#eta|<0.3","p");
1549   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1550   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1551   leg->SetFillColor(0);
1552   leg->Draw();
1553
1554   delete [] effsA;
1555   delete [] effsB;
1556   delete [] effsC;
1557   delete [] pt;
1558
1559   return;
1560 }
1561 //-----------------------------------------------------------------------------
1562 void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg,
1563                                    Int_t par) {
1564 //-----------------------------------------------------------------------------
1565 // This function draws the pulls in the [pT,eta] bins
1566 //-----------------------------------------------------------------------------
1567
1568   ReadPulls(inName);
1569   SetParticle(pdg);
1570
1571   const Int_t n = (fPulls+par)->GetPointsPt();
1572   Double_t *pullsA = new Double_t[n];
1573   Double_t *pullsB = new Double_t[n];
1574   Double_t *pullsC = new Double_t[n];
1575   Double_t *pt     = new Double_t[n];
1576   (fPulls+par)->GetArrayPt(pt);
1577   for(Int_t i=0;i<n;i++) {
1578     pullsA[i] = (fPulls+par)->GetParam(i,0);
1579     pullsB[i] = (fPulls+par)->GetParam(i,1);
1580     pullsC[i] = (fPulls+par)->GetParam(i,2);
1581   }
1582
1583   TGraph *grA = new TGraph(n,pt,pullsA);
1584   TGraph *grB = new TGraph(n,pt,pullsB);
1585   TGraph *grC = new TGraph(n,pt,pullsC);
1586
1587   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1588   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1589   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1590
1591   TString title("Distribution of the pulls: ");
1592   switch (pdg) {
1593   case 211:
1594     title.Prepend("PIONS - ");
1595     break;
1596   case 321:
1597     title.Prepend("KAONS - ");
1598     break;
1599   case 2212:
1600     title.Prepend("PROTONS - ");
1601     break;
1602   case 11:
1603     title.Prepend("ELECTRONS - ");
1604     break;
1605   case 13:
1606     title.Prepend("MUONS - ");
1607     break;
1608   }
1609   switch (par) {
1610   case 0:
1611     title.Append("y");
1612     break;
1613   case 1:
1614     title.Append("z");
1615     break;
1616   case 2:
1617     title.Append("    #eta");
1618     break;
1619   case 3:
1620     title.Append("tg    #lambda");
1621     break;
1622   case 4:
1623     title.Append("C");
1624     break;
1625   }
1626
1627   gStyle->SetOptStat(0);
1628   TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1629   c->SetLogx();
1630   c->SetGridx();
1631   c->SetGridy();
1632
1633   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1634   frame->SetXTitle("p_{T} [GeV/c]");
1635   frame->Draw();
1636   grA->Draw("P");
1637   grB->Draw("P");
1638   grC->Draw("P");
1639
1640   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1641   leg->AddEntry(grA,"|#eta|<0.3","p");
1642   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1643   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1644   leg->SetFillColor(0);
1645   leg->Draw();
1646
1647   delete [] pullsA;
1648   delete [] pullsB;
1649   delete [] pullsC;
1650   delete [] pt;
1651
1652   return;
1653 }
1654 //-----------------------------------------------------------------------------
1655 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1656                                                Double_t eta) const {
1657 //-----------------------------------------------------------------------------
1658 // This function stretches the covariance matrix according to the pulls
1659 //-----------------------------------------------------------------------------
1660
1661   TMatrixD covMat(5,5);
1662
1663   covMat(0,0)=cc[0];
1664   covMat(1,0)=cc[1];  covMat(0,1)=covMat(1,0);
1665   covMat(1,1)=cc[2];
1666   covMat(2,0)=cc[3];  covMat(0,2)=covMat(2,0);
1667   covMat(2,1)=cc[4];  covMat(1,2)=covMat(2,1);
1668   covMat(2,2)=cc[5];
1669   covMat(3,0)=cc[6];  covMat(0,3)=covMat(3,0);
1670   covMat(3,1)=cc[7];  covMat(1,3)=covMat(3,1);
1671   covMat(3,2)=cc[8];  covMat(2,3)=covMat(3,2);
1672   covMat(3,3)=cc[9];
1673   covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1674   covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1675   covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1676   covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1677   covMat(4,4)=cc[14];
1678
1679
1680   TMatrixD stretchMat(5,5);
1681   for(Int_t k=0;k<5;k++) {
1682     for(Int_t l=0;l<5;l++) {
1683       stretchMat(k,l)=0.;
1684     }
1685   }
1686
1687   for(Int_t i=0;i<5;i++) {
1688     stretchMat(i,i) = (fPulls+i)->GetValueAt(pt,eta);
1689     if(stretchMat(i,i)==0.) stretchMat(i,i) = 1.;
1690   }
1691
1692   TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1693   TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1694
1695   return covMatSmear;
1696 }
1697 //-----------------------------------------------------------------------------
1698 void AliTPCtrackerParam::InitializeKineGrid(Option_t* which) {
1699 //-----------------------------------------------------------------------------
1700 // This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1701 // which = "DB"     -> initialize fDBgrid... members
1702 //         "eff"    -> initialize fEff... members
1703 //         "pulls"  -> initialize fPulls... members
1704 //         "dEdx"   -> initialize fdEdx... members
1705 //-----------------------------------------------------------------------------
1706
1707   const char *DB     = strstr(which,"DB");
1708   const char *eff    = strstr(which,"eff");
1709   const char *pulls  = strstr(which,"pulls");
1710   const char *dEdx   = strstr(which,"dEdx");
1711
1712
1713   Int_t nEta=0, nPt=0;
1714
1715   Double_t etaPoints[2] = {0.3,0.6};
1716   Double_t etaBins[3]   = {0.15,0.45,0.75};
1717
1718   Double_t ptPoints[9]  = {0.4,0.6,0.8,1.2,1.7,3.,5.,8.,15.};
1719   Double_t ptBins[10]  = {0.3,0.5,0.7,1.,1.5,2.,4.,6.,10.,20.};
1720
1721
1722   Double_t *eta=0,*pt=0;
1723
1724   if(DB) {
1725     nEta = 2;
1726     nPt  = 9;
1727     eta  = etaPoints;
1728     pt   = ptPoints;
1729   } else {
1730     nEta = 3;
1731     nPt  = 10;
1732     eta  = etaBins;
1733     pt   = ptBins;
1734   }
1735
1736   AliTPCkineGrid *dummy=0;
1737
1738   if(DB) {
1739     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1740     new(&fDBgridPi) AliTPCkineGrid(*dummy);
1741     new(&fDBgridKa) AliTPCkineGrid(*dummy);
1742     new(&fDBgridPr) AliTPCkineGrid(*dummy);
1743     new(&fDBgridEl) AliTPCkineGrid(*dummy);
1744     new(&fDBgridMu) AliTPCkineGrid(*dummy);
1745     delete dummy;
1746   }
1747   if(eff) {
1748     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1749     new(&fEffPi) AliTPCkineGrid(*dummy);
1750     new(&fEffKa) AliTPCkineGrid(*dummy);
1751     new(&fEffPr) AliTPCkineGrid(*dummy);
1752     new(&fEffEl) AliTPCkineGrid(*dummy);
1753     new(&fEffMu) AliTPCkineGrid(*dummy);
1754     delete dummy;
1755   }
1756   if(pulls) {
1757     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1758     for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1759     for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1760     for(Int_t i=0;i<5;i++) new(&fPullsPr[i]) AliTPCkineGrid(*dummy);
1761     for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1762     for(Int_t i=0;i<5;i++) new(&fPullsMu[i]) AliTPCkineGrid(*dummy);
1763     delete dummy;
1764   }
1765   if(dEdx) {
1766     dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
1767     new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1768     new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1769     new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1770     new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1771     new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1772     new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1773     new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1774     new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1775     new(&fdEdxMeanMu) AliTPCkineGrid(*dummy);
1776     new(&fdEdxRMSMu) AliTPCkineGrid(*dummy);
1777     delete dummy;
1778   }
1779
1780   return;
1781 }
1782 //-----------------------------------------------------------------------------
1783 void AliTPCtrackerParam::MakeDataBase() {
1784 //-----------------------------------------------------------------------------
1785 // This function creates the DB file and store in it:
1786 //  - TPC Efficiencies for pi,ka,pr,el,mu     (AliTPCkineGrid class)
1787 //  - Pulls for pi,ka,el                      (AliTPCkineGrid class)
1788 //  - Regularization parameters for pi,ka,el  (TMatrixD class)
1789 //  - dE/dx parameterization for pi,ka,pr,el  (AliTPCkineGrid class)
1790 //  - Regularized cov. matrices for pi,ka,el  (COVMATRIX structure)
1791 //-----------------------------------------------------------------------------
1792
1793   // define some file names
1794   Char_t *effFile   ="TPCeff.root";
1795   Char_t *pullsFile ="pulls.root";
1796   Char_t *regPiFile ="regPi.root";
1797   Char_t *regKaFile ="regKa.root";
1798   Char_t *regPrFile ="regPr.root";
1799   Char_t *regElFile ="regEl.root";
1800   Char_t *regMuFile ="regMu.root";
1801   Char_t *dEdxPiFile="dEdxPi.root";
1802   Char_t *dEdxKaFile="dEdxKa.root";
1803   Char_t *dEdxPrFile="dEdxPr.root";
1804   Char_t *dEdxElFile="dEdxEl.root";
1805   Char_t *dEdxMuFile="dEdxMu.root";
1806   Char_t *cmFile    ="CovMatrix_AllEvts.root";
1807   /*
1808   Char_t *cmFile1  ="CovMatrix_AllEvts_1.root";
1809   Char_t *cmFile2  ="CovMatrix_AllEvts_2.root";
1810   Char_t *cmFile3  ="CovMatrix_AllEvts_3.root";
1811   */
1812
1813   // store the effieciencies
1814   ReadEffs(effFile);
1815   WriteEffs(fDBfileName.Data());
1816
1817   // store the pulls
1818   ReadPulls(pullsFile);
1819   WritePulls(fDBfileName.Data());
1820
1821   //store the regularization parameters
1822   ReadRegParams(regPiFile,211);
1823   WriteRegParams(fDBfileName.Data(),211);
1824   ReadRegParams(regKaFile,321);
1825   WriteRegParams(fDBfileName.Data(),321);
1826   ReadRegParams(regPrFile,2212);
1827   WriteRegParams(fDBfileName.Data(),2212);
1828   ReadRegParams(regElFile,11);
1829   WriteRegParams(fDBfileName.Data(),11);
1830   ReadRegParams(regMuFile,13);
1831   WriteRegParams(fDBfileName.Data(),13);
1832
1833   // store the dEdx parameters
1834   ReaddEdx(dEdxPiFile,211);
1835   WritedEdx(fDBfileName.Data(),211);
1836   ReaddEdx(dEdxKaFile,321);
1837   WritedEdx(fDBfileName.Data(),321);
1838   ReaddEdx(dEdxPrFile,2212);
1839   WritedEdx(fDBfileName.Data(),2212);
1840   ReaddEdx(dEdxElFile,11);
1841   WritedEdx(fDBfileName.Data(),11);
1842   ReaddEdx(dEdxMuFile,13);
1843   WritedEdx(fDBfileName.Data(),13);
1844
1845
1846   //
1847   // store the regularized covariance matrices
1848   //
1849   InitializeKineGrid("DB");
1850
1851   const Int_t nBinsPi = fDBgridPi.GetTotBins();
1852   const Int_t nBinsKa = fDBgridKa.GetTotBins();
1853   const Int_t nBinsPr = fDBgridPr.GetTotBins();
1854   const Int_t nBinsEl = fDBgridEl.GetTotBins();
1855   const Int_t nBinsMu = fDBgridMu.GetTotBins();
1856
1857
1858   // create the trees for cov. matrices
1859   // trees for pions
1860   TTree *CovTreePi_ = NULL;
1861   CovTreePi_ = new TTree[nBinsPi];
1862   // trees for kaons
1863   TTree *CovTreeKa_ = NULL;
1864   CovTreeKa_ = new TTree[nBinsKa];
1865   // trees for protons
1866   TTree *CovTreePr_ = NULL;
1867   CovTreePr_ = new TTree[nBinsPr];
1868   // trees for electrons
1869   TTree *CovTreeEl_ = NULL;
1870   CovTreeEl_ = new TTree[nBinsEl];
1871   // trees for muons
1872   TTree *CovTreeMu_ = NULL;
1873   CovTreeMu_ = new TTree[nBinsMu];
1874
1875   Char_t hname[100], htitle[100];
1876   COVMATRIX covmat;
1877
1878
1879   for(Int_t i=0; i<nBinsPi; i++) {
1880     sprintf(hname,"CovTreePi_bin%d",i);
1881     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1882     CovTreePi_[i].SetName(hname); CovTreePi_[i].SetTitle(htitle);
1883     CovTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1884   }
1885   for(Int_t i=0; i<nBinsKa; i++) {
1886     sprintf(hname,"CovTreeKa_bin%d",i);
1887     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1888     CovTreeKa_[i].SetName(hname); CovTreeKa_[i].SetTitle(htitle);
1889     CovTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1890   }
1891   for(Int_t i=0; i<nBinsPr; i++) {
1892     sprintf(hname,"CovTreePr_bin%d",i);
1893     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1894     CovTreePr_[i].SetName(hname); CovTreePr_[i].SetTitle(htitle);
1895     CovTreePr_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1896   }
1897   for(Int_t i=0; i<nBinsEl; i++) {
1898     sprintf(hname,"CovTreeEl_bin%d",i);
1899     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1900     CovTreeEl_[i].SetName(hname); CovTreeEl_[i].SetTitle(htitle);
1901     CovTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1902   }
1903   for(Int_t i=0; i<nBinsMu; i++) {
1904     sprintf(hname,"CovTreeMu_bin%d",i);
1905     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1906     CovTreeMu_[i].SetName(hname); CovTreeMu_[i].SetTitle(htitle);
1907     CovTreeMu_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1908   }
1909
1910   /*
1911   // create the chain with the compared tracks
1912   TChain *cmptrktree = new TChain("cmptrktree");
1913   cmptrkchain.Add(cmFile1);
1914   cmptrkchain.Add(cmFile2);
1915   cmptrkchain.Add(cmFile3);
1916   */
1917
1918   TFile *infile = TFile::Open(cmFile);
1919   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
1920
1921   COMPTRACK cmptrk;
1922   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
1923   Int_t entries = (Int_t)cmptrktree->GetEntries();
1924   cerr<<" Number of entries: "<<entries<<endl;
1925
1926   Int_t trkPdg,trkBin;
1927   Double_t trkKine[1],trkRegPar[3];
1928   Int_t *nPerBinPi = new Int_t[nBinsPi];
1929   for(Int_t k=0;k<nBinsPi;k++) nPerBinPi[k]=0;
1930   Int_t *nPerBinKa = new Int_t[nBinsKa];
1931   for(Int_t k=0;k<nBinsKa;k++) nPerBinKa[k]=0;
1932   Int_t *nPerBinMu = new Int_t[nBinsMu];
1933   for(Int_t k=0;k<nBinsMu;k++) nPerBinMu[k]=0;
1934   Int_t *nPerBinEl = new Int_t[nBinsEl];
1935   for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
1936   Int_t *nPerBinPr = new Int_t[nBinsPr];
1937   for(Int_t k=0;k<nBinsPr;k++) nPerBinPr[k]=0;
1938
1939   // loop on chain entries
1940   for(Int_t l=0; l<entries; l++) {
1941     if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1942     // get the event
1943     cmptrktree->GetEvent(l);
1944     // get the pdg code
1945     trkPdg = TMath::Abs(cmptrk.pdg);
1946     // use only pions, kaons, protons, electrons, muons
1947     if(trkPdg!=211 && trkPdg!=321 && trkPdg!=2212 && trkPdg!=11 && trkPdg!=13) continue;
1948     SetParticle(trkPdg);
1949     trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1950     //cerr<<cmptrk.pt<<"  "<<cmptrk.eta<<"  "<<trkBin<<endl;
1951
1952     if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1953     if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1954     if(trkPdg==2212 && nPerBinPr[trkBin]>=5000) continue;
1955     if(trkPdg==11  && nPerBinEl[trkBin]>=5000) continue;
1956     if(trkPdg==13 && nPerBinMu[trkBin]>=5000) continue;
1957
1958     trkKine[0] = cmptrk.p;
1959
1960     // get regularized covariance matrix
1961     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(0,k);
1962     covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1963     covmat.c10 = cmptrk.c10;
1964     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(1,k);
1965     covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1966     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(2,k);
1967     covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1968     covmat.c21 = cmptrk.c21;
1969     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(3,k);
1970     covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
1971     covmat.c30 = cmptrk.c30;
1972     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(4,k);
1973     covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
1974     covmat.c32 = cmptrk.c32;
1975     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(5,k);
1976     covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
1977     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(6,k);
1978     covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
1979     covmat.c41 = cmptrk.c41;
1980     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(7,k);
1981     covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
1982     covmat.c43 = cmptrk.c43;
1983     for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(8,k);
1984     covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
1985
1986     // fill the tree
1987     switch (trkPdg) {
1988     case 211: // pions
1989       CovTreePi_[trkBin].Fill();
1990       nPerBinPi[trkBin]++;
1991       break;
1992     case 321: // kaons
1993       CovTreeKa_[trkBin].Fill();
1994       nPerBinKa[trkBin]++;
1995       break;
1996     case 2212: // protons
1997       CovTreePr_[trkBin].Fill();
1998       nPerBinPr[trkBin]++;
1999       break;
2000     case 11: // electrons
2001       CovTreeEl_[trkBin].Fill();
2002       nPerBinEl[trkBin]++;
2003       break;
2004     case 13: // muons
2005       CovTreeMu_[trkBin].Fill();
2006       nPerBinMu[trkBin]++;
2007       break;
2008     }
2009   } // loop on chain entries
2010
2011   // store all trees the DB file
2012   TFile *DBfile = new TFile(fDBfileName.Data(),"update");
2013   DBfile->mkdir("CovMatrices");
2014   gDirectory->cd("/CovMatrices");
2015   gDirectory->mkdir("Pions");
2016   gDirectory->mkdir("Kaons");
2017   gDirectory->mkdir("Protons");
2018   gDirectory->mkdir("Electrons");
2019   gDirectory->mkdir("Muons");
2020   // store pions
2021   gDirectory->cd("/CovMatrices/Pions");
2022   fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
2023   for(Int_t i=0;i<nBinsPi;i++) CovTreePi_[i].Write();
2024   // store kaons
2025   gDirectory->cd("/CovMatrices/Kaons");
2026   fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
2027   for(Int_t i=0;i<nBinsKa;i++) CovTreeKa_[i].Write();
2028   // store kaons
2029   gDirectory->cd("/CovMatrices/Protons");
2030   fDBgridPr.SetName("DBgridPr"); fDBgridPr.Write();
2031   for(Int_t i=0;i<nBinsPr;i++) CovTreePr_[i].Write();
2032   // store electrons
2033   gDirectory->cd("/CovMatrices/Electrons");
2034   fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
2035   for(Int_t i=0;i<nBinsEl;i++) CovTreeEl_[i].Write();
2036   // store kaons
2037   gDirectory->cd("/CovMatrices/Muons");
2038   fDBgridMu.SetName("DBgridMu"); fDBgridMu.Write();
2039   for(Int_t i=0;i<nBinsMu;i++) CovTreeMu_[i].Write();
2040
2041   DBfile->Close();
2042   delete [] nPerBinPi;
2043   delete [] nPerBinKa;
2044   delete [] nPerBinPr;
2045   delete [] nPerBinEl;
2046   delete [] nPerBinMu;
2047
2048   return;
2049 }
2050 //-----------------------------------------------------------------------------
2051 void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *TPC,TTree *TH,
2052                                            TObjArray &seedArray) const {
2053 //-----------------------------------------------------------------------------
2054 // This function makes the seeds for tracks from the 1st hits in the TPC
2055 //-----------------------------------------------------------------------------
2056
2057   Double_t xg,yg,zg,px,py,pz,pt;
2058   Int_t label;
2059   Int_t nTracks=(Int_t)TH->GetEntries();
2060
2061   cerr<<"+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
2062          "\n+++\n\n";
2063
2064   AliTPChit *tpcHit=0;
2065
2066   // loop over entries in TreeH
2067   for(Int_t l=0; l<nTracks; l++) {
2068     if(l%1000==0) cerr<<"  --- Processing primary track "
2069                       <<l<<" of "<<nTracks<<" ---\r";
2070     TPC->ResetHits();
2071     TH->GetEvent(l);
2072     // Get FirstHit
2073     tpcHit=(AliTPChit*)TPC->FirstHit(-1);
2074     for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
2075       if(tpcHit->fQ !=0.) continue;
2076       // Get particle momentum at hit
2077       px=tpcHit->X(); py=tpcHit->Y(); pz=tpcHit->Z();
2078
2079       pt=TMath::Sqrt(px*px+py*py);
2080       // reject hits with Pt<mag*0.45 GeV/c
2081       if(pt<(fBz*0.45)) continue;
2082
2083       // Get track label
2084       label=tpcHit->Track();
2085
2086       if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
2087       if(tpcHit->fQ != 0.) continue;
2088       // Get global coordinates of hit
2089       xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
2090       if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
2091
2092       // create the seed
2093       AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
2094
2095       // reject tracks which are not in the TPC acceptance
2096       if(!ioseed->InTPCAcceptance()) { delete ioseed; continue; }
2097
2098       // put seed in array
2099       seedArray.AddLast(ioseed);
2100     }
2101
2102   } // loop over entries in TreeH
2103
2104   return;
2105 }
2106 //-----------------------------------------------------------------------------
2107 void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *TTR,
2108                                            TObjArray &seedArray) const {
2109 //-----------------------------------------------------------------------------
2110 // This function makes the seeds for tracks from the track references
2111 //-----------------------------------------------------------------------------
2112
2113   Double_t xg,yg,zg,px,py,pz,pt;
2114   Int_t label;
2115   Int_t nnn,k;
2116
2117   TClonesArray *tkRefArray = new TClonesArray("AliTrackReference");
2118
2119   TBranch *b =(TBranch*)TTR->GetBranch("TPC");
2120   if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
2121   b->SetAddress(&tkRefArray);
2122   Int_t nTkRef = (Int_t)b->GetEntries();
2123   cerr<<"+++\n+++ Number of entries in TreeTR(TPC): "<<nTkRef<<
2124          "\n+++\n\n";
2125
2126   // loop on track references
2127   for(Int_t l=0; l<nTkRef; l++){
2128     if(l%1000==0) cerr<<"  --- Processing primary track "
2129                       <<l<<" of "<<nTkRef<<" ---\r";
2130     if(!b->GetEvent(l)) continue;
2131     nnn = tkRefArray->GetEntriesFast();
2132
2133     if(nnn <= 0) continue;
2134     for(k=0; k<nnn; k++) {
2135       AliTrackReference *tkref = (AliTrackReference*)tkRefArray->UncheckedAt(k);
2136
2137       xg = tkref->X();
2138       yg = tkref->Y();
2139       zg = tkref->Z();
2140       px = tkref->Px();
2141       py = tkref->Py();
2142       pz = tkref->Pz();
2143       label = tkref->GetTrack();
2144
2145       pt=TMath::Sqrt(px*px+py*py);
2146       // reject hits with Pt<mag*0.45 GeV/c
2147       if(pt<(fBz*0.45)) continue;
2148
2149       // create the seed
2150       AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
2151
2152       // reject if not at the inner part of TPC
2153       if(TMath::Abs(ioseed->GetXL()-82.9701) > 0.1) {
2154         delete ioseed; continue;
2155       }
2156
2157       // reject tracks which are not in the TPC acceptance
2158       if(!ioseed->InTPCAcceptance()) {
2159         delete ioseed; continue;
2160       }
2161
2162       // put seed in array
2163       seedArray.AddLast(ioseed);
2164
2165     }
2166
2167
2168   } // loop on track references
2169
2170   delete tkRefArray;
2171
2172
2173   return;
2174 }
2175 //-----------------------------------------------------------------------------
2176 void AliTPCtrackerParam::MergeEvents(Int_t evFirst,Int_t evLast) {
2177 //-----------------------------------------------------------------------------
2178 // This function: 1) merges the files from track comparison
2179 //                   (beware: better no more than 100 events per file)
2180 //                2) computes the average TPC efficiencies
2181 //-----------------------------------------------------------------------------
2182
2183   Char_t *outName="TPCeff.root";
2184
2185   // merge files with tracks
2186   cerr<<" ******** MERGING FILES **********\n\n";
2187
2188   // create the chain for the tree of compared tracks
2189   TChain ch1("cmptrktree");
2190   TChain ch2("cmptrktree");
2191   TChain ch3("cmptrktree");
2192
2193   for(Int_t j=evFirst; j<=evLast; j++) {
2194     cerr<<"Processing event "<<j<<endl;
2195
2196     TString covName("CovMatrix.");
2197     covName+=j;
2198     covName.Append(".root");
2199
2200     if(gSystem->AccessPathName(covName.Data(),kFileExists)) continue;
2201
2202
2203     if(j<=100) ch1.Add(covName.Data());
2204     if(j>100 && j<=200) ch2.Add(covName.Data());
2205     if(j>200) ch3.Add(covName.Data());
2206
2207   } // loop on events
2208
2209   // merge chain in one file
2210   TFile *covOut=0;
2211   covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
2212   ch1.Merge(covOut,1000000000);
2213   covOut->Close();
2214   delete covOut;
2215   covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
2216   ch2.Merge(covOut,1000000000);
2217   covOut->Close();
2218   delete covOut;
2219   covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
2220   ch3.Merge(covOut,1000000000);
2221   covOut->Close();
2222   delete covOut;
2223
2224
2225   // efficiencies
2226   cerr<<" ***** EFFICIENCIES ******\n\n";
2227
2228   ReadEffs("TPCeff.1.root");
2229
2230   Int_t n = fEffPi.GetTotPoints();
2231   Double_t *avEffPi = new Double_t[n];
2232   Double_t *avEffKa = new Double_t[n];
2233   Double_t *avEffPr = new Double_t[n];
2234   Double_t *avEffEl = new Double_t[n];
2235   Double_t *avEffMu = new Double_t[n];
2236   Int_t    *evtsPi  = new Int_t[n];
2237   Int_t    *evtsKa  = new Int_t[n];
2238   Int_t    *evtsPr  = new Int_t[n];
2239   Int_t    *evtsEl  = new Int_t[n];
2240   Int_t    *evtsMu  = new Int_t[n];
2241
2242   for(Int_t j=0; j<n; j++) {
2243     avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
2244     evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
2245   }
2246
2247   for(Int_t j=evFirst; j<=evLast; j++) {
2248     cerr<<"Processing event "<<j<<endl;
2249
2250     TString effName("TPCeff.");
2251     effName+=j;
2252     effName.Append(".root");
2253
2254     if(gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
2255
2256     ReadEffs(effName.Data());
2257
2258     for(Int_t k=0; k<n; k++) {
2259       if(fEffPi.GetParam(k)>=0.) {avEffPi[k]+=fEffPi.GetParam(k); evtsPi[k]++;}
2260       if(fEffKa.GetParam(k)>=0.) {avEffKa[k]+=fEffKa.GetParam(k); evtsKa[k]++;}
2261       if(fEffPr.GetParam(k)>=0.) {avEffPr[k]+=fEffPr.GetParam(k); evtsPr[k]++;}
2262       if(fEffEl.GetParam(k)>=0.) {avEffEl[k]+=fEffEl.GetParam(k); evtsEl[k]++;}
2263       if(fEffMu.GetParam(k)>=0.) {avEffMu[k]+=fEffMu.GetParam(k); evtsMu[k]++;}
2264     }
2265
2266   } // loop on events
2267
2268   // compute average efficiencies
2269   for(Int_t j=0; j<n; j++) {
2270     if(evtsPi[j]==0) evtsPi[j]++;
2271     fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
2272     if(evtsKa[j]==0) evtsKa[j]++;
2273     fEffKa.SetParam(j,(Double_t)avEffKa[j]/evtsKa[j]);
2274     if(evtsPr[j]==0) evtsPr[j]++;
2275     fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
2276     if(evtsEl[j]==0) evtsEl[j]++;
2277     fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
2278     if(evtsMu[j]==0) evtsMu[j]++;
2279     fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
2280   }
2281
2282   // write efficiencies to a file
2283   WriteEffs(outName);
2284
2285   delete [] avEffPi;
2286   delete [] avEffKa;
2287   delete [] avEffPr;
2288   delete [] avEffEl;
2289   delete [] avEffMu;
2290   delete [] evtsPi;
2291   delete [] evtsKa;
2292   delete [] evtsPr;
2293   delete [] evtsEl;
2294   delete [] evtsMu;
2295
2296   return;
2297 }
2298 //-----------------------------------------------------------------------------
2299 Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
2300 //-----------------------------------------------------------------------------
2301 // This function reads all parameters from the DB
2302 //-----------------------------------------------------------------------------
2303
2304   if(!ReadEffs(inName)) return 0;
2305   if(!ReadPulls(inName)) return 0;
2306   if(!ReadRegParams(inName,211)) return 0;
2307   if(!ReadRegParams(inName,321)) return 0;
2308   if(!ReadRegParams(inName,2212)) return 0;
2309   if(!ReadRegParams(inName,11)) return 0;
2310   if(!ReadRegParams(inName,13)) return 0;
2311   if(!ReaddEdx(inName,211)) return 0;
2312   if(!ReaddEdx(inName,321)) return 0;
2313   if(!ReaddEdx(inName,2212)) return 0;
2314   if(!ReaddEdx(inName,11)) return 0;
2315   if(!ReaddEdx(inName,13)) return 0;
2316   if(!ReadDBgrid(inName)) return 0;
2317
2318   return 1;
2319 }
2320 //-----------------------------------------------------------------------------
2321 Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
2322 //-----------------------------------------------------------------------------
2323 // This function reads the dEdx parameters from the DB
2324 //-----------------------------------------------------------------------------
2325
2326   if(gSystem->AccessPathName(inName,kFileExists)) {
2327     cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
2328     return 0; }
2329
2330   TFile *inFile = TFile::Open(inName);
2331   switch (pdg) {
2332   case 211:
2333     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
2334     fdEdxMeanPi.~AliTPCkineGrid();
2335     new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
2336     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
2337     fdEdxRMSPi.~AliTPCkineGrid();
2338     new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
2339     break;
2340   case 321:
2341     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
2342     fdEdxMeanKa.~AliTPCkineGrid();
2343     new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
2344     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
2345     fdEdxRMSKa.~AliTPCkineGrid();
2346     new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
2347     break;
2348   case 2212:
2349     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
2350     fdEdxMeanPr.~AliTPCkineGrid();
2351     new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
2352     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
2353     fdEdxRMSPr.~AliTPCkineGrid();
2354     new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
2355     break;
2356   case 11:
2357     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
2358     fdEdxMeanEl.~AliTPCkineGrid();
2359     new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
2360     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
2361     fdEdxRMSEl.~AliTPCkineGrid();
2362     new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
2363     break;
2364   case 13:
2365     if(fColl==0) {
2366     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
2367     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
2368     } else {
2369       fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxMeanMu");
2370       fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxRMSMu");
2371     }
2372     fdEdxMeanMu.~AliTPCkineGrid();
2373     new(&fdEdxMeanMu) AliTPCkineGrid(*fdEdxMean);
2374     fdEdxRMSMu.~AliTPCkineGrid();
2375     new(&fdEdxRMSMu) AliTPCkineGrid(*fdEdxRMS);
2376     break;
2377   }
2378   inFile->Close();
2379
2380   return 1;
2381 }
2382 //-----------------------------------------------------------------------------
2383 Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
2384 //-----------------------------------------------------------------------------
2385 // This function reads the kine grid from the DB
2386 //-----------------------------------------------------------------------------
2387
2388   if(gSystem->AccessPathName(inName,kFileExists)) {
2389     cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n";
2390     return 0; }
2391
2392   TFile *inFile = TFile::Open(inName);
2393
2394   // first read the DB grid for the different particles
2395   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2396   fDBgridPi.~AliTPCkineGrid();
2397   new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
2398   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
2399   fDBgridKa.~AliTPCkineGrid();
2400   new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
2401   if(fColl==0) {
2402     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2403   } else {
2404     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Protons/DBgridPr");
2405   }
2406   fDBgridPr.~AliTPCkineGrid();
2407   new(&fDBgridPr) AliTPCkineGrid(*fDBgrid);
2408   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
2409   fDBgridEl.~AliTPCkineGrid();
2410   new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
2411   if(fColl==0) {
2412     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
2413   } else {
2414     fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Muons/DBgridMu");
2415   }
2416   fDBgridMu.~AliTPCkineGrid();
2417   new(&fDBgridMu) AliTPCkineGrid(*fDBgrid);
2418
2419   inFile->Close();
2420
2421   return 1;
2422 }
2423 //-----------------------------------------------------------------------------
2424 Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
2425 //-----------------------------------------------------------------------------
2426 // This function reads the TPC efficiencies from the DB
2427 //-----------------------------------------------------------------------------
2428
2429   if(gSystem->AccessPathName(inName,kFileExists)) {
2430     cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n";
2431     return 0; }
2432
2433   TFile *inFile = TFile::Open(inName);
2434
2435   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
2436   fEffPi.~AliTPCkineGrid();
2437   new(&fEffPi) AliTPCkineGrid(*fEff);
2438   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
2439   fEffKa.~AliTPCkineGrid();
2440   new(&fEffKa) AliTPCkineGrid(*fEff);
2441   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
2442   fEffPr.~AliTPCkineGrid();
2443   new(&fEffPr) AliTPCkineGrid(*fEff);
2444   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
2445   fEffEl.~AliTPCkineGrid();
2446   new(&fEffEl) AliTPCkineGrid(*fEff);
2447   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
2448   fEffMu.~AliTPCkineGrid();
2449   new(&fEffMu) AliTPCkineGrid(*fEff);
2450
2451   inFile->Close();
2452
2453   return 1;
2454 }
2455 //-----------------------------------------------------------------------------
2456 Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
2457 //-----------------------------------------------------------------------------
2458 // This function reads the pulls from the DB
2459 //-----------------------------------------------------------------------------
2460
2461   if(gSystem->AccessPathName(inName,kFileExists)) {
2462     cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n";
2463     return 0; }
2464
2465   TFile *inFile = TFile::Open(inName);
2466
2467   for(Int_t i=0; i<5; i++) {
2468     TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
2469     TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
2470     TString pr("/Pulls/Protons/PullsPr_"); pr+=i;
2471     TString el("/Pulls/Electrons/PullsEl_"); el+=i;
2472     TString mu("/Pulls/Muons/PullsMu_"); mu+=i;
2473
2474     fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2475     fPullsPi[i].~AliTPCkineGrid();
2476     new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
2477
2478     fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
2479     fPullsKa[i].~AliTPCkineGrid();
2480     new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
2481
2482     if(fColl==0) {
2483       fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2484     } else {
2485       fPulls = (AliTPCkineGrid*)inFile->Get(pr.Data());
2486     }
2487     fPullsPr[i].~AliTPCkineGrid();
2488     new(&fPullsPr[i]) AliTPCkineGrid(*fPulls);
2489
2490     fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
2491     fPullsEl[i].~AliTPCkineGrid();
2492     new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
2493
2494     if(fColl==0) {
2495       fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
2496     } else {
2497       fPulls = (AliTPCkineGrid*)inFile->Get(mu.Data());
2498     }
2499     fPullsMu[i].~AliTPCkineGrid();
2500     new(&fPullsMu[i]) AliTPCkineGrid(*fPulls);
2501   }
2502
2503   inFile->Close();
2504
2505   return 1;
2506 }
2507 //-----------------------------------------------------------------------------
2508 Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
2509 //-----------------------------------------------------------------------------
2510 // This function reads the regularization parameters from the DB
2511 //-----------------------------------------------------------------------------
2512
2513   if(gSystem->AccessPathName(inName,kFileExists)) {
2514     cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
2515     return 0; }
2516
2517   TFile *inFile = TFile::Open(inName);
2518   switch (pdg) {
2519   case 211:
2520     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2521     new(&fRegParPi) TMatrixD(*fRegPar);
2522     break;
2523   case 321:
2524     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
2525     new(&fRegParKa) TMatrixD(*fRegPar);
2526     break;
2527   case 2212:
2528     if(fColl==0) {
2529       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2530     } else {
2531       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
2532     }
2533     new(&fRegParPr) TMatrixD(*fRegPar);
2534     break;
2535   case 11:
2536     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
2537     new(&fRegParEl) TMatrixD(*fRegPar);
2538     break;
2539   case 13:
2540     if(fColl==0) {
2541       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
2542     } else {
2543       fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
2544     }
2545     new(&fRegParMu) TMatrixD(*fRegPar);
2546     break;
2547   }
2548   inFile->Close();
2549
2550   return 1;
2551 }
2552 //-----------------------------------------------------------------------------
2553 void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
2554 //-----------------------------------------------------------------------------
2555 // This function regularizes the elements of the covariance matrix
2556 // that show a momentum depence:
2557 // c00, c11, c22, c33, c44, c20, c24, c40, c31
2558 // The regularization is done separately for pions, kaons, electrons:
2559 // give "Pion","Kaon" or "Electron" as first argument.
2560 //-----------------------------------------------------------------------------
2561
2562   gStyle->SetOptStat(0);
2563   gStyle->SetOptFit(10001);
2564
2565   Int_t thisPdg=211;
2566   Char_t *part="Pions - ";
2567
2568   InitializeKineGrid("DB");
2569   SetParticle(pdg);
2570   const Int_t fitbins = fDBgrid->GetBinsPt();
2571   cerr<<" Fit bins:  "<<fitbins<<endl;
2572
2573   switch (pdg) {
2574   case 211: // pions
2575     thisPdg=211;
2576     part="Pions - ";
2577     cerr<<" Processing pions ...\n";
2578     break;
2579   case 321: //kaons
2580     thisPdg=321;
2581     part="Kaons - ";
2582     cerr<<" Processing kaons ...\n";
2583     break;
2584   case 2212: //protons
2585     thisPdg=2212;
2586     part="Protons - ";
2587     cerr<<" Processing protons ...\n";
2588     break;
2589   case 11: // electrons
2590     thisPdg= 11;
2591     part="Electrons - ";
2592     cerr<<" Processing electrons ...\n";
2593     break;
2594   case 13: // muons
2595     thisPdg= 13;
2596     part="Muons - ";
2597     cerr<<" Processing muons ...\n";
2598     break;
2599   }
2600
2601
2602   /*
2603   // create a chain with compared tracks
2604   TChain *cmptrkchain = new ("cmptrktree");
2605   cmptrkchain.Add("CovMatrix_AllEvts.root");
2606   //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
2607   //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
2608   //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
2609   */
2610
2611
2612   TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
2613   TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
2614
2615   COMPTRACK cmptrk;
2616   cmptrktree->SetBranchAddress("comptracks",&cmptrk);
2617   Int_t entries = (Int_t)cmptrktree->GetEntries();
2618
2619
2620   Int_t pbin;
2621   Int_t    *n       = new Int_t[fitbins];
2622   Int_t    *n00     = new Int_t[fitbins];
2623   Int_t    *n11     = new Int_t[fitbins];
2624   Int_t    *n20     = new Int_t[fitbins];
2625   Int_t    *n22     = new Int_t[fitbins];
2626   Int_t    *n31     = new Int_t[fitbins];
2627   Int_t    *n33     = new Int_t[fitbins];
2628   Int_t    *n40     = new Int_t[fitbins];
2629   Int_t    *n42     = new Int_t[fitbins];
2630   Int_t    *n44     = new Int_t[fitbins];
2631   Double_t *p       = new Double_t[fitbins];
2632   Double_t *ep      = new Double_t[fitbins];
2633   Double_t *mean00  = new Double_t[fitbins];
2634   Double_t *mean11  = new Double_t[fitbins];
2635   Double_t *mean20  = new Double_t[fitbins];
2636   Double_t *mean22  = new Double_t[fitbins];
2637   Double_t *mean31  = new Double_t[fitbins];
2638   Double_t *mean33  = new Double_t[fitbins];
2639   Double_t *mean40  = new Double_t[fitbins];
2640   Double_t *mean42  = new Double_t[fitbins];
2641   Double_t *mean44  = new Double_t[fitbins];
2642   Double_t *sigma00 = new Double_t[fitbins];
2643   Double_t *sigma11 = new Double_t[fitbins];
2644   Double_t *sigma20 = new Double_t[fitbins];
2645   Double_t *sigma22 = new Double_t[fitbins];
2646   Double_t *sigma31 = new Double_t[fitbins];
2647   Double_t *sigma33 = new Double_t[fitbins];
2648   Double_t *sigma40 = new Double_t[fitbins];
2649   Double_t *sigma42 = new Double_t[fitbins];
2650   Double_t *sigma44 = new Double_t[fitbins];
2651   Double_t *rmean   = new Double_t[fitbins];
2652   Double_t *rsigma  = new Double_t[fitbins];
2653   Double_t fitpar[3];
2654
2655   for(Int_t l=0; l<fitbins; l++) {
2656     n[l]=1;
2657     n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
2658     p[l ]=ep[l]=0.;
2659     mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
2660     sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
2661   }
2662
2663   // loop on chain entries for mean
2664   for(Int_t l=0; l<entries; l++) {
2665     cmptrktree->GetEvent(l);
2666     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2667     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2668     n[pbin]++;
2669     p[pbin]+=cmptrk.p;
2670     mean00[pbin]+=cmptrk.c00;
2671     mean11[pbin]+=cmptrk.c11;
2672     mean20[pbin]+=cmptrk.c20;
2673     mean22[pbin]+=cmptrk.c22;
2674     mean31[pbin]+=cmptrk.c31;
2675     mean33[pbin]+=cmptrk.c33;
2676     mean40[pbin]+=cmptrk.c40;
2677     mean42[pbin]+=cmptrk.c42;
2678     mean44[pbin]+=cmptrk.c44;
2679   } // loop on chain entries
2680
2681   for(Int_t l=0; l<fitbins; l++) {
2682     p[l]/=n[l];
2683     mean00[l]/=n[l];
2684     mean11[l]/=n[l];
2685     mean20[l]/=n[l];
2686     mean22[l]/=n[l];
2687     mean31[l]/=n[l];
2688     mean33[l]/=n[l];
2689     mean40[l]/=n[l];
2690     mean42[l]/=n[l];
2691     mean44[l]/=n[l];
2692   }
2693
2694   // loop on chain entries for sigma
2695   for(Int_t l=0; l<entries; l++) {
2696     cmptrktree->GetEvent(l);
2697     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2698     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2699     if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2700       sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2701     if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2702       sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2703     if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2704       sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2705     if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2706       sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2707     if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2708       sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2709     if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2710       sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2711     if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2712       sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2713     if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2714       sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2715     if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2716       sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2717   } // loop on chain entries
2718
2719   for(Int_t l=0; l<fitbins; l++) {
2720     sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2721     sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2722     sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2723     sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2724     sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2725     sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2726     sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2727     sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2728     sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2729   }
2730
2731
2732   // Fit function
2733   TF1 *func = new TF1("RegFunc",RegFunc,0.23,50.,3);
2734   func->SetParNames("A_meas","A_scatt","B");
2735
2736   // line to draw on the plots
2737   TLine *lin = new TLine(-1,1,1.69,1);
2738   lin->SetLineStyle(2);
2739   lin->SetLineWidth(2);
2740
2741   // matrix used to store fit results
2742   TMatrixD fitRes(9,3);
2743
2744   //    --- c00 ---
2745
2746   // create the canvas
2747   TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
2748   canv00->Divide(1,2);
2749   // create the graph for cov matrix
2750   TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
2751   TString title00("C(y,y)"); title00.Prepend(part);
2752   TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
2753   frame00->SetXTitle("p [GeV/c]");
2754   canv00->cd(1);  gPad->SetLogx();
2755   frame00->Draw();
2756   gr00->Draw("P");
2757   // Sets initial values for parameters
2758   func->SetParameters(1.6e-3,1.9e-4,1.5);
2759   // Fit points in range defined by function
2760   gr00->Fit("RegFunc","R,Q");
2761   func->GetParameters(fitpar);
2762   for(Int_t i=0; i<3; i++) fitRes(0,i)=fitpar[i];
2763   for(Int_t l=0; l<fitbins; l++) {
2764     rmean[l]  = mean00[l]/RegFunc(&p[l],fitpar);
2765     rsigma[l] = sigma00[l]/RegFunc(&p[l],fitpar);
2766   }
2767   // create the graph the regularized cov. matrix
2768   TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2769   TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B})");
2770   regtitle00.Prepend(part);
2771   TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
2772   frame00reg->SetXTitle("p [GeV/c]");
2773   canv00->cd(2); gPad->SetLogx();
2774   frame00reg->Draw();
2775   gr00reg->Draw("P");
2776   lin->Draw("same");
2777
2778
2779   //    --- c11 ---
2780
2781   // create the canvas
2782   TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
2783   canv11->Divide(1,2);
2784   // create the graph for cov matrix
2785   TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
2786   TString title11("C(z,z)"); title11.Prepend(part);
2787   TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
2788   frame11->SetXTitle("p [GeV/c]");
2789   canv11->cd(1);  gPad->SetLogx();
2790   frame11->Draw();
2791   gr11->Draw("P");
2792   // Sets initial values for parameters
2793   func->SetParameters(1.2e-3,8.1e-4,1.);
2794   // Fit points in range defined by function
2795   gr11->Fit("RegFunc","R,Q");
2796   func->GetParameters(fitpar);
2797   for(Int_t i=0; i<3; i++) fitRes(1,i)=fitpar[i];
2798   for(Int_t l=0; l<fitbins; l++) {
2799     rmean[l]  = mean11[l]/RegFunc(&p[l],fitpar);
2800     rsigma[l] = sigma11[l]/RegFunc(&p[l],fitpar);
2801   }
2802   // create the graph the regularized cov. matrix
2803   TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2804   TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B})");
2805   regtitle11.Prepend(part);
2806   TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
2807   frame11reg->SetXTitle("p [GeV/c]");
2808   canv11->cd(2); gPad->SetLogx();
2809   frame11reg->Draw();
2810   gr11reg->Draw("P");
2811   lin->Draw("same");
2812
2813
2814   //    --- c20 ---
2815
2816   // create the canvas
2817   TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
2818   canv20->Divide(1,2);
2819   // create the graph for cov matrix
2820   TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
2821   TString title20("C(#eta, y)"); title20.Prepend(part);
2822   TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
2823   frame20->SetXTitle("p [GeV/c]");
2824   canv20->cd(1);  gPad->SetLogx();
2825   frame20->Draw();
2826   gr20->Draw("P");
2827   // Sets initial values for parameters
2828   func->SetParameters(7.3e-5,1.2e-5,1.5);
2829   // Fit points in range defined by function
2830   gr20->Fit("RegFunc","R,Q");
2831   func->GetParameters(fitpar);
2832   for(Int_t i=0; i<3; i++) fitRes(2,i)=fitpar[i];
2833   for(Int_t l=0; l<fitbins; l++) {
2834     rmean[l]  = mean20[l]/RegFunc(&p[l],fitpar);
2835     rsigma[l] = sigma20[l]/RegFunc(&p[l],fitpar);
2836   }
2837   // create the graph the regularized cov. matrix
2838   TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2839   TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B})");
2840   regtitle20.Prepend(part);
2841   TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
2842   frame20reg->SetXTitle("p [GeV/c]");
2843   canv20->cd(2); gPad->SetLogx();
2844   frame20reg->Draw();
2845   gr20reg->Draw("P");
2846   lin->Draw("same");
2847
2848
2849   //    --- c22 ---
2850
2851   // create the canvas
2852   TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
2853   canv22->Divide(1,2);
2854   // create the graph for cov matrix
2855   TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
2856   TString title22("C(#eta, #eta)"); title22.Prepend(part);
2857   TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
2858   frame22->SetXTitle("p [GeV/c]");
2859   canv22->cd(1);  gPad->SetLogx();
2860   frame22->Draw();
2861   gr22->Draw("P");
2862   // Sets initial values for parameters
2863   func->SetParameters(5.2e-6,1.1e-6,2.);
2864   // Fit points in range defined by function
2865   gr22->Fit("RegFunc","R,Q");
2866   func->GetParameters(fitpar);
2867   for(Int_t i=0; i<3; i++) fitRes(3,i)=fitpar[i];
2868   for(Int_t l=0; l<fitbins; l++) {
2869     rmean[l]  = mean22[l]/RegFunc(&p[l],fitpar);
2870     rsigma[l] = sigma22[l]/RegFunc(&p[l],fitpar);
2871   }
2872   // create the graph the regularized cov. matrix
2873   TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2874   TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B})");
2875   regtitle22.Prepend(part);
2876   TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
2877   frame22reg->SetXTitle("p [GeV/c]");
2878   canv22->cd(2); gPad->SetLogx();
2879   frame22reg->Draw();
2880   gr22reg->Draw("P");
2881   lin->Draw("same");
2882
2883
2884   //    --- c31 ---
2885
2886   // create the canvas
2887   TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
2888   canv31->Divide(1,2);
2889   // create the graph for cov matrix
2890   TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
2891   TString title31("C(tg #lambda,z)"); title31.Prepend(part);
2892   TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
2893   frame31->SetXTitle("p [GeV/c]");
2894   canv31->cd(1);  gPad->SetLogx();
2895   frame31->Draw();
2896   gr31->Draw("P");
2897   // Sets initial values for parameters
2898   func->SetParameters(-1.2e-5,-1.2e-5,1.5);
2899   // Fit points in range defined by function
2900   gr31->Fit("RegFunc","R,Q");
2901   func->GetParameters(fitpar);
2902   for(Int_t i=0; i<3; i++) fitRes(4,i)=fitpar[i];
2903   for(Int_t l=0; l<fitbins; l++) {
2904     rmean[l]  = mean31[l]/RegFunc(&p[l],fitpar);
2905     rsigma[l] = -sigma31[l]/RegFunc(&p[l],fitpar);
2906   }
2907   // create the graph the regularized cov. matrix
2908   TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2909   TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B})");
2910   regtitle31.Prepend(part);
2911   TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
2912   frame31reg->SetXTitle("p [GeV/c]");
2913   canv31->cd(2); gPad->SetLogx();
2914   frame31reg->Draw();
2915   gr31reg->Draw("P");
2916   lin->Draw("same");
2917
2918
2919   //    --- c33 ---
2920
2921   // create the canvas
2922   TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
2923   canv33->Divide(1,2);
2924   // create the graph for cov matrix
2925   TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
2926   TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
2927   TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
2928   frame33->SetXTitle("p [GeV/c]");
2929   canv33->cd(1);  gPad->SetLogx();
2930   frame33->Draw();
2931   gr33->Draw("P");
2932   // Sets initial values for parameters
2933   func->SetParameters(1.3e-7,4.6e-7,1.7);
2934   // Fit points in range defined by function
2935   gr33->Fit("RegFunc","R,Q");
2936   func->GetParameters(fitpar);
2937   for(Int_t i=0; i<3; i++) fitRes(5,i)=fitpar[i];
2938   for(Int_t l=0; l<fitbins; l++) {
2939     rmean[l]  = mean33[l]/RegFunc(&p[l],fitpar);
2940     rsigma[l] = sigma33[l]/RegFunc(&p[l],fitpar);
2941   }
2942   // create the graph the regularized cov. matrix
2943   TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2944   TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B})");
2945   regtitle33.Prepend(part);
2946   TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
2947   frame33reg->SetXTitle("p [GeV/c]");
2948   canv33->cd(2); gPad->SetLogx();
2949   frame33reg->Draw();
2950   gr33reg->Draw("P");
2951   lin->Draw("same");
2952
2953
2954   //    --- c40 ---
2955
2956   // create the canvas
2957   TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
2958   canv40->Divide(1,2);
2959   // create the graph for cov matrix
2960   TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
2961   TString title40("C(C,y)"); title40.Prepend(part);
2962   TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
2963   frame40->SetXTitle("p [GeV/c]");
2964   canv40->cd(1);  gPad->SetLogx();
2965   frame40->Draw();
2966   gr40->Draw("P");
2967   // Sets initial values for parameters
2968   func->SetParameters(4.e-7,4.4e-8,1.5);
2969   // Fit points in range defined by function
2970   gr40->Fit("RegFunc","R,Q");
2971   func->GetParameters(fitpar);
2972   for(Int_t i=0; i<3; i++) fitRes(6,i)=fitpar[i];
2973   for(Int_t l=0; l<fitbins; l++) {
2974     rmean[l]  = mean40[l]/RegFunc(&p[l],fitpar);
2975     rsigma[l] = sigma40[l]/RegFunc(&p[l],fitpar);
2976   }
2977   // create the graph the regularized cov. matrix
2978   TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2979   TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B})");
2980   regtitle40.Prepend(part);
2981   TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
2982   frame40reg->SetXTitle("p [GeV/c]");
2983   canv40->cd(2); gPad->SetLogx();
2984   frame40reg->Draw();
2985   gr40reg->Draw("P");
2986   lin->Draw("same");
2987
2988
2989   //    --- c42 ---
2990
2991   // create the canvas
2992   TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
2993   canv42->Divide(1,2);
2994   // create the graph for cov matrix
2995   TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
2996   TString title42("C(C, #eta)"); title42.Prepend(part);
2997   TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
2998   frame42->SetXTitle("p [GeV/c]");
2999   canv42->cd(1);  gPad->SetLogx();
3000   frame42->Draw();
3001   gr42->Draw("P");
3002   // Sets initial values for parameters
3003   func->SetParameters(3.e-8,8.2e-9,2.);
3004   // Fit points in range defined by function
3005   gr42->Fit("RegFunc","R,Q");
3006   func->GetParameters(fitpar);
3007   for(Int_t i=0; i<3; i++) fitRes(7,i)=fitpar[i];
3008   for(Int_t l=0; l<fitbins; l++) {
3009     rmean[l]  = mean42[l]/RegFunc(&p[l],fitpar);
3010     rsigma[l] = sigma42[l]/RegFunc(&p[l],fitpar);
3011   }
3012   // create the graph the regularized cov. matrix
3013   TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
3014   TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B})");
3015   regtitle42.Prepend(part);
3016   TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
3017   frame42reg->SetXTitle("p [GeV/c]");
3018   canv42->cd(2); gPad->SetLogx();
3019   frame42reg->Draw();
3020   gr42reg->Draw("P");
3021   lin->Draw("same");
3022
3023
3024   //    --- c44 ---
3025
3026   // create the canvas
3027   TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
3028   canv44->Divide(1,2);
3029   // create the graph for cov matrix
3030   TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
3031   TString title44("C(C,C)"); title44.Prepend(part);
3032   TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
3033   frame44->SetXTitle("p [GeV/c]");
3034   canv44->cd(1);  gPad->SetLogx();
3035   frame44->Draw();
3036   gr44->Draw("P");
3037   // Sets initial values for parameters
3038   func->SetParameters(1.8e-10,5.8e-11,2.);
3039   // Fit points in range defined by function
3040   gr44->Fit("RegFunc","R,Q");
3041   func->GetParameters(fitpar);
3042   for(Int_t i=0; i<3; i++) fitRes(8,i)=fitpar[i];
3043   for(Int_t l=0; l<fitbins; l++) {
3044     rmean[l]  = mean44[l]/RegFunc(&p[l],fitpar);
3045     rsigma[l] = sigma44[l]/RegFunc(&p[l],fitpar);
3046   }
3047   // create the graph the regularized cov. matrix
3048   TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
3049   TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B})");
3050   regtitle44.Prepend(part);
3051   TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
3052   frame44reg->SetXTitle("p [GeV/c]");
3053   canv44->cd(2); gPad->SetLogx();
3054   frame44reg->Draw();
3055   gr44reg->Draw("P");
3056   lin->Draw("same");
3057
3058
3059
3060
3061   switch (pdg) {
3062   case 211:
3063     new(&fRegParPi) TMatrixD(fitRes);
3064     break;
3065   case 321:
3066     new(&fRegParKa) TMatrixD(fitRes);
3067     break;
3068   case 2212:
3069     new(&fRegParPr) TMatrixD(fitRes);
3070     break;
3071   case 11:
3072     new(&fRegParEl) TMatrixD(fitRes);
3073     break;
3074   case 13:
3075     new(&fRegParMu) TMatrixD(fitRes);
3076     break;
3077   }
3078
3079   // write fit parameters to file
3080   WriteRegParams(outName,pdg);
3081
3082   delete [] n;
3083   delete [] n00;
3084   delete [] n11;
3085   delete [] n20;
3086   delete [] n22;
3087   delete [] n31;
3088   delete [] n33;
3089   delete [] n40;
3090   delete [] n42;
3091   delete [] n44;
3092   delete [] p;
3093   delete [] ep;
3094   delete [] mean00;
3095   delete [] mean11;
3096   delete [] mean20;
3097   delete [] mean22;
3098   delete [] mean31;
3099   delete [] mean33;
3100   delete [] mean40;
3101   delete [] mean42;
3102   delete [] mean44;
3103   delete [] sigma00;
3104   delete [] sigma11;
3105   delete [] sigma20;
3106   delete [] sigma22;
3107   delete [] sigma31;
3108   delete [] sigma33;
3109   delete [] sigma40;
3110   delete [] sigma42;
3111   delete [] sigma44;
3112   delete [] rmean;
3113   delete [] rsigma;
3114
3115   return;
3116 }
3117 //-----------------------------------------------------------------------------
3118 Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
3119 //-----------------------------------------------------------------------------
3120 // This function makes a selection according to TPC tracking efficiency
3121 //-----------------------------------------------------------------------------
3122
3123   Double_t eff=0.;
3124
3125   eff = fEff->GetValueAt(pt,eta);
3126
3127   if(gRandom->Rndm() < eff) return kTRUE;
3128
3129   return kFALSE;
3130 }
3131 //-----------------------------------------------------------------------------
3132 void AliTPCtrackerParam::SetParticle(Int_t pdg) {
3133 //-----------------------------------------------------------------------------
3134 // This function set efficiencies, pulls, etc... for the particle
3135 // specie of the current particle
3136 //-----------------------------------------------------------------------------
3137
3138   switch (pdg) {
3139   case 211:
3140     fDBgrid = &fDBgridPi;
3141     fEff    = &fEffPi;
3142     fPulls  =  fPullsPi;
3143     fRegPar = &fRegParPi;
3144     fdEdxMean = &fdEdxMeanPi;
3145     fdEdxRMS  = &fdEdxRMSPi;
3146     break;
3147   case 321:
3148     fDBgrid = &fDBgridKa;
3149     fEff    = &fEffKa;
3150     fPulls  =  fPullsKa;
3151     fRegPar = &fRegParKa;
3152     fdEdxMean = &fdEdxMeanKa;
3153     fdEdxRMS  = &fdEdxRMSKa;
3154     break;
3155   case 2212:
3156     fDBgrid = &fDBgridPr;
3157     fEff    = &fEffPr;
3158     fPulls  =  fPullsPr;
3159     fRegPar = &fRegParPr;
3160     fdEdxMean = &fdEdxMeanPr;
3161     fdEdxRMS  = &fdEdxRMSPr;
3162     break;
3163   case 11:
3164     fDBgrid = &fDBgridEl;
3165     fEff    = &fEffEl;
3166     fPulls  =  fPullsEl;
3167     fRegPar = &fRegParEl;
3168     fdEdxMean = &fdEdxMeanEl;
3169     fdEdxRMS  = &fdEdxRMSEl;
3170     break;
3171   case 13:
3172     fDBgrid = &fDBgridMu;
3173     fEff    = &fEffMu;
3174     fPulls  =  fPullsMu;
3175     fRegPar = &fRegParMu;
3176     fdEdxMean = &fdEdxMeanMu;
3177     fdEdxRMS  = &fdEdxRMSMu;
3178     break;
3179   }
3180
3181   return;
3182 }
3183 //-----------------------------------------------------------------------------
3184 void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
3185   const {
3186 //-----------------------------------------------------------------------------
3187 // This function smears track parameters according to streched cov. matrix
3188 //-----------------------------------------------------------------------------
3189   AliGausCorr *corgen = new AliGausCorr(cov,5);
3190   TArrayD corr(5);
3191   corgen->GetGaussN(corr);
3192   delete corgen;
3193   corgen = 0;
3194
3195   for(Int_t l=0;l<5;l++) {
3196     xxsm[l] = xx[l]+corr[l];
3197   }
3198
3199   return;
3200 }
3201 //-----------------------------------------------------------------------------
3202 Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
3203 //-----------------------------------------------------------------------------
3204 // This function writes the dEdx parameters to the DB
3205 //-----------------------------------------------------------------------------
3206
3207   Option_t *opt;
3208   Char_t *dirName="Pions";
3209   Char_t *meanName="dEdxMeanPi";
3210   Char_t *RMSName="dEdxRMSPi";
3211
3212   SetParticle(pdg);
3213
3214   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
3215   } else { opt="update"; }
3216
3217   switch (pdg) {
3218   case 211:
3219     dirName="Pions";
3220     meanName="dEdxMeanPi";
3221     RMSName="dEdxRMSPi";
3222     break;
3223   case 321:
3224     dirName="Kaons";
3225     meanName="dEdxMeanKa";
3226     RMSName="dEdxRMSKa";
3227     break;
3228   case 2212:
3229     dirName="Protons";
3230     meanName="dEdxMeanPr";
3231     RMSName="dEdxRMSPr";
3232     break;
3233   case 11:
3234     dirName="Electrons";
3235     meanName="dEdxMeanEl";
3236     RMSName="dEdxRMSEl";
3237     break;
3238   case 13:
3239     dirName="Muons";
3240     meanName="dEdxMeanMu";
3241     RMSName="dEdxRMSMu";
3242     break;
3243   }
3244
3245   TFile *outFile = new TFile(outName,opt);
3246   if(!gDirectory->cd("/dEdx")) {
3247     outFile->mkdir("dEdx");
3248     gDirectory->cd("/dEdx");
3249   }
3250   TDirectory *dir2 = gDirectory->mkdir(dirName);
3251   dir2->cd();
3252   fdEdxMean->SetName(meanName); fdEdxMean->Write();
3253   fdEdxRMS->SetName(RMSName);  fdEdxRMS->Write();
3254
3255   outFile->Close();
3256   delete outFile;
3257
3258
3259   return 1;
3260 }
3261 //-----------------------------------------------------------------------------
3262 Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
3263 //-----------------------------------------------------------------------------
3264 // This function writes the TPC efficiencies to the DB
3265 //-----------------------------------------------------------------------------
3266
3267
3268   Option_t *opt;
3269
3270   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
3271   } else { opt="update"; }
3272
3273   TFile *outFile = new TFile(outName,opt);
3274
3275   outFile->mkdir("Efficiencies");
3276   gDirectory->cd("/Efficiencies");
3277   gDirectory->mkdir("Pions");
3278   gDirectory->mkdir("Kaons");
3279   gDirectory->mkdir("Protons");
3280   gDirectory->mkdir("Electrons");
3281   gDirectory->mkdir("Muons");
3282
3283   gDirectory->cd("/Efficiencies/Pions");
3284   fEffPi.SetName("EffPi");
3285   fEffPi.Write();
3286   gDirectory->cd("/Efficiencies/Kaons");
3287   fEffKa.SetName("EffKa");
3288   fEffKa.Write();
3289   gDirectory->cd("/Efficiencies/Protons");
3290   fEffPr.SetName("EffPr");
3291   fEffPr.Write();
3292   gDirectory->cd("/Efficiencies/Electrons");
3293   fEffEl.SetName("EffEl");
3294   fEffEl.Write();
3295   gDirectory->cd("/Efficiencies/Muons");
3296   fEffMu.SetName("EffMu");
3297   fEffMu.Write();
3298
3299   outFile->Close();
3300
3301   delete outFile;
3302
3303   return 1;
3304 }
3305 //-----------------------------------------------------------------------------
3306 Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
3307 //-----------------------------------------------------------------------------
3308 // This function writes the pulls to the DB
3309 //-----------------------------------------------------------------------------
3310
3311   Option_t *opt;
3312
3313   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
3314   } else { opt="update"; }
3315
3316   TFile *outFile = new TFile(outName,opt);
3317
3318   outFile->mkdir("Pulls");
3319   gDirectory->cd("/Pulls");
3320   gDirectory->mkdir("Pions");
3321   gDirectory->mkdir("Kaons");
3322   gDirectory->mkdir("Protons");
3323   gDirectory->mkdir("Electrons");
3324   gDirectory->mkdir("Muons");
3325
3326   for(Int_t i=0;i<5;i++) {
3327     TString pi("PullsPi_"); pi+=i;
3328     TString ka("PullsKa_"); ka+=i;
3329     TString pr("PullsPr_"); pr+=i;
3330     TString el("PullsEl_"); el+=i;
3331     TString mu("PullsMu_"); mu+=i;
3332     fPullsPi[i].SetName(pi.Data());
3333     fPullsKa[i].SetName(ka.Data());
3334     fPullsPr[i].SetName(pr.Data());
3335     fPullsEl[i].SetName(el.Data());
3336     fPullsMu[i].SetName(mu.Data());
3337     gDirectory->cd("/Pulls/Pions");
3338     fPullsPi[i].Write();
3339     gDirectory->cd("/Pulls/Kaons");
3340     fPullsKa[i].Write();
3341     gDirectory->cd("/Pulls/Protons");
3342     fPullsPr[i].Write();
3343     gDirectory->cd("/Pulls/Electrons");
3344     fPullsEl[i].Write();
3345     gDirectory->cd("/Pulls/Muons");
3346     fPullsMu[i].Write();
3347   }
3348   outFile->Close();
3349   delete outFile;
3350
3351   return 1;
3352 }
3353 //-----------------------------------------------------------------------------
3354 Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
3355 //-----------------------------------------------------------------------------
3356 // This function writes the regularization parameters to the DB
3357 //-----------------------------------------------------------------------------
3358
3359   Option_t *opt;
3360   Char_t *dirName="Pions";
3361   Char_t *keyName="RegPions";
3362
3363   SetParticle(pdg);
3364
3365   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
3366   } else { opt="update"; }
3367
3368   switch (pdg) {
3369   case 211:
3370     dirName="Pions";
3371     keyName="RegPions";
3372     break;
3373   case 321:
3374     dirName="Kaons";
3375     keyName="RegKaons";
3376     break;
3377   case 2212:
3378     dirName="Protons";
3379     keyName="RegProtons";
3380     break;
3381   case 11:
3382     dirName="Electrons";
3383     keyName="RegElectrons";
3384     break;
3385   case 13:
3386     dirName="Muons";
3387     keyName="RegMuons";
3388     break;
3389   }
3390
3391   TFile *outFile = new TFile(outName,opt);
3392   if(!gDirectory->cd("/RegParams")) {
3393     outFile->mkdir("RegParams");
3394     gDirectory->cd("/RegParams");
3395   }
3396   TDirectory *dir2 = gDirectory->mkdir(dirName);
3397   dir2->cd();
3398   fRegPar->Write(keyName);
3399
3400   outFile->Close();
3401   delete outFile;
3402
3403
3404   return 1;
3405 }
3406 //-----------------------------------------------------------------------------
3407 //*****************************************************************************
3408 //-----------------------------------------------------------------------------
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419
3420
3421
3422
3423
3424
3425
3426
3427
3428
3429
3430
3431