TPC/AliTPCtrackerParam.cxx

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