TPC/AliTPCtrackerParam.cxx

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