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