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