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