TPC/AliTPCtrackerParam.cxx

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