TPC/AliTPCtrackerParam.cxx

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