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