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