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