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