TPC/AliTPCtrackerParam.cxx

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