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 /*
  17 $Log$
  18 Revision 1.5.4.1  2002/06/10 15:26:12  hristov
  19 Merged with v3-08-02
  20
  21 Revision 1.9  2002/05/13 09:53:08  hristov
  22 Some frequent problems corrected: arrays with variable size have to be defined via the operator new, default values for the arguments have to be  used only in the header files, etc.
  23
  24 Revision 1.8  2002/05/08 18:21:40  kowal2
  25 Now the code is blind to the binning used for pulls, efficiencies etc.
  26
  27 Revision 1.7  2002/04/10 16:30:07  kowal2
  28 logs added
  29
  30 */
  31
  32
  33 /**************************************************************************
  34  *                                                                        *
  35  * This class builds AliTPCtrack objects from generated tracks to feed    *
  36  * ITS tracking (V2). The AliTPCtrack is built from its first hit in      *
  37  * the TPC. The track is assigned a Kalman-like covariance matrix         *
  38  * depending on its pT and pseudorapidity and track parameters are        *
  39  * smeared according to this covariance matrix.                           *
  40  * Output file contains sorted tracks, ready for matching with ITS.       *
  41  *                                                                        *
  42  * For details:                                                           *
  43  * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm    *
  44  *                                                                        *
  45  * Test macro is: AliBarrelRec_TPCparam.C                                 *
  46  *                                                                        *
  47  * 2002/04/28: Major upgrade of the class                                 *
  48  * - Covariance matrices and pulls are now separeted for pions, kaons     *
  49  *   and electrons.                                                       *
  50  * - A parameterization of dE/dx in the TPC has been included and it is   *
  51  *   used to assign a mass to each track according to a rough dE/dx PID.  *
  52  * - All the "numbers" have been moved to the file with the DataBase,     *
  53  *   they are read as objects of the class AliTPCkineGrid, and assigned   *
  54  *   to data memebers of the class AliTPCtrackerParam.                    *
  55  * - All the code necessary to create a BataBase has been included in     *
  56  *   class (see the macro AliTPCtrackingParamDB.C for the details).       *
  57  *                                                                        *
  58  *  Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it    *
  59  *                                                                        *
  60  **************************************************************************/
  61 #include <TChain.h>
  62 #include <TF1.h>
  63 #include <TGraph.h>
  64 #include <TGraphErrors.h>
  65 #include <TLegend.h>
  66 #include <TLine.h>
  67 #include <TMatrixD.h>
  68 #include <TNtuple.h>
  69 #include <TSystem.h>
  70 #include "alles.h"
  71 #include "AliGausCorr.h"
  72 #include "AliKalmanTrack.h"
  73 #include "AliMagF.h"
  74 #include "AliMagFCM.h"
  75 #include "AliTPCkineGrid.h"
  76 #include "AliTPCtrack.h"
  77 #include "AliTPCtrackerParam.h"
  78
  79 Double_t RegFunc(Double_t *x,Double_t *par) {
  80 // This is the function used to regularize the covariance matrix
  81   Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
  82                    TMath::Power(x[1],par[3]);
  83   return value;
  84 }
  85 Double_t FitRegFunc(Double_t *x,Double_t *par) {
  86 // This is the function used for the fit of the covariance
  87 // matrix as a function of the momentum.
  88 // For cosl the average value 0.87 is used.
  89   Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
  90                    TMath::Power(0.87,par[3]);
  91   return value;
  92 }
  93
  94 // structure for DB building
  95 typedef struct {
  96   Int_t    pdg;
  97   Int_t    bin;
  98   Double_t r;
  99   Double_t p;
 100   Double_t pt;
 101   Double_t cosl;
 102   Double_t eta;
 103   Double_t dpt;
 104   Double_t dP0,dP1,dP2,dP3,dP4;
 105   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
 106   Double_t dEdx;
 107 } COMPTRACK;
 108 // cov matrix structure
 109 typedef struct {
 110   Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
 111 } COVMATRIX;
 112
 113 ClassImp(AliTPCtrackerParam)
 114
 115 //-----------------------------------------------------------------------------
 116 AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
 117 {
 118 //-----------------------------------------------------------------------------
 119 // This is the class conctructor
 120 //-----------------------------------------------------------------------------
 121
 122   fBz = Bz;             // value of the z component of L3 field (Tesla)
 123   fColl = coll;         // collision code (0: PbPb6000)
 124   fSelAndSmear = kTRUE; // by default selection and smearing are done
 125
 126   if(fBz!=0.4) {
 127     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid field!\n";
 128     cerr<<"      Available:  0.4"<<endl;
 129   }
 130   if(fColl!=0) {
 131     cerr<<"AliTPCtrackerParam::AliTPCtrackerParam:  Invalid collision!\n";
 132     cerr<<"      Available:  0   ->   PbPb6000"<<endl;
 133   }
 134
 135   fDBfileName = gSystem->Getenv("ALICE_ROOT");
 136   fDBfileName.Append("/TPC/CovMatrixDB_");
 137   if(fColl==0) fDBfileName.Append("PbPb6000");
 138   if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
 139 }
 140 //-----------------------------------------------------------------------------
 141 AliTPCtrackerParam::~AliTPCtrackerParam()
 142 {}
 143 //-----------------------------------------------------------------------------
 144
 145 Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
 146 {
 147 //-----------------------------------------------------------------------------
 148 // This function creates the TPC parameterized tracks
 149 //-----------------------------------------------------------------------------
 150
 151   TFile *fileDB=0;
 152   TTree *covTreePi[50];
 153   TTree *covTreeKa[50];
 154   TTree *covTreeEl[50];
 155
 156   if(fSelAndSmear) {
 157     cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
 158       fDBfileName.Data()<<"\n+++\n";
 159     // Read paramters from DB file
 160     if(!ReadAllData(fDBfileName.Data())) {
 161       cerr<<"AliTPCtrackerParam::BuildTPCtracks: \
 162              Could not read data from DB\n\n"; return 1;
 163     }
 164     // Read the trees with regularized cov. matrices from DB
 165     TString str;
 166     fileDB = TFile::Open(fDBfileName.Data());
 167     Int_t nBinsPi = fDBgridPi.GetTotBins();
 168     for(Int_t l=0; l<nBinsPi; l++) {
 169       str = "/CovMatrices/Pions/CovTreePi_bin";
 170       str+=l;
 171       covTreePi[l] = (TTree*)fileDB->Get(str.Data());
 172     }
 173     Int_t nBinsKa = fDBgridKa.GetTotBins();
 174     for(Int_t l=0; l<nBinsKa; l++) {
 175       str = "/CovMatrices/Kaons/CovTreeKa_bin";
 176       str+=l;
 177       covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
 178     }
 179     Int_t nBinsEl = fDBgridEl.GetTotBins();
 180     for(Int_t l=0; l<nBinsEl; l++) {
 181       str = "/CovMatrices/Electrons/CovTreeEl_bin";
 182       str+=l;
 183       covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
 184     }
 185
 186   } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
 187
 188   TFile *infile=(TFile*)inp;
 189   infile->cd();
 190
 191   // Get gAlice object from file
 192   if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
 193     cerr<<"gAlice has not been found on galice.root !\n";
 194     return 1;
 195   }
 196
 197   // Check if value in the galice file is equal to selected one (fBz)
 198   AliMagFCM *fiel = (AliMagFCM*)gAlice->Field();
 199   Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
 200   printf("Magnetic field is %6.2f Tesla\n",fieval);
 201   if(fBz!=fieval) {
 202     cerr<<"AliTPCtrackerParam::BuildTPCtracks:  Invalid field!"<<endl;
 203     cerr<<"Field selected is: "<<fBz<<" T\n";
 204     cerr<<"Field found on file is: "<<fieval<<" T\n";
 205     return 0;
 206   }
 207
 208   // Get TPC detector
 209   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
 210   Int_t ver = TPC->IsVersion();
 211   cerr<<"+++ TPC version "<<ver<<" has been found !\n";
 212   AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
 213   if(digp){
 214     delete digp;
 215     digp = new AliTPCParamSR();
 216   }
 217   else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
 218
 219   if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
 220   TPC->SetParam(digp);
 221
 222   // Set the conversion constant between curvature and Pt
 223   AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
 224
 225   TParticle   *Part=0;
 226   AliTPChit   *tpcHit=0;
 227   AliTPCtrack *tpctrack=0;
 228   Double_t     hPx,hPy,hPz,hPt,hEta,xg,yg,zg,xl,yl,zl;
 229   Double_t     alpha;
 230   Float_t      cosAlpha,sinAlpha;
 231   Int_t        bin,label,pdg,charge;
 232   Int_t        tracks=0;
 233   Int_t        nParticles,nTracks,arrentr;
 234   Char_t       tname[100];
 235   //Int_t nSel=0,nAcc=0;
 236
 237   // loop over first n events in file
 238   for(Int_t evt=0; evt<n; evt++){
 239     cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
 240
 241     // tree for TPC tracks
 242     sprintf(tname,"TreeT_TPC_%d",evt);
 243     TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
 244     tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
 245
 246     // array for TPC tracks
 247     TObjArray tarray(20000);
 248
 249     // get the particles stack
 250     nParticles = gAlice->GetEvent(evt);
 251     Bool_t * done = new Bool_t[nParticles];
 252     Int_t  * pdgCodes = new Int_t[nParticles];
 253
 254     // loop on particles and store pdg codes
 255     for(Int_t l=0; l<nParticles; l++) {
 256       Part        = gAlice->Particle(l);
 257       pdgCodes[l] = Part->GetPdgCode();
 258       done[l]     = kFALSE;
 259     }
 260
 261     // Get TreeH with hits
 262     TTree *TH=gAlice->TreeH();
 263     nTracks=(Int_t)TH->GetEntries();
 264     cerr<<"+++\n+++ Number of particles in event "<<evt<<":  "<<nParticles<<
 265       "\n+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
 266       "\n+++\n\n";
 267
 268     // loop over entries in TreeH
 269     for(Int_t l=0; l<nTracks; l++) {
 270       if(l%1000==0) cerr<<"  --- Processing primary track "
 271                         <<l<<" of "<<nTracks<<" ---\r";
 272       TPC->ResetHits();
 273       TH->GetEvent(l);
 274       // Get FirstHit
 275       tpcHit=(AliTPChit*)TPC->FirstHit(-1);
 276       for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
 277         if(tpcHit->fQ !=0.) continue;
 278         // Get particle momentum at hit
 279         hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
 280         hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
 281         // reject hits with Pt<mag*0.45 GeV/c
 282         if(hPt<(fBz*0.45)) continue;
 283
 284         // Get track label
 285         label=tpcHit->Track();
 286         // check if this track has already been processed
 287         if(done[label]) continue;
 288         // PDG code & electric charge
 289         pdg = pdgCodes[label];
 290         if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
 291         else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
 292         else continue;
 293         pdg = TMath::Abs(pdg);
 294         if(pdg>3000) pdg=211;
 295         if(fSelAndSmear) SetParticle(pdg);
 296
 297         if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
 298         if(tpcHit->fQ != 0.) continue;
 299         // Get global coordinates of hit
 300         xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
 301         if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
 302
 303         // Get TPC sector, Alpha angle and local coordinates
 304         // cerr<<"Sector "<<tpcHit->fSector<<endl;
 305         digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
 306         alpha = TMath::ATan2(sinAlpha,cosAlpha);
 307         xl = xg*cosAlpha + yg*sinAlpha;
 308         yl =-xg*sinAlpha + yg*cosAlpha;
 309         zl = zg;
 310         //printf("Alpha %f   xl %f  yl %f  zl %f\n",Alpha,xl,yl,zl);
 311
 312         // reject tracks which are not in the TPC acceptance
 313         if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
 314
 315         hEta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(hPz/hPt)));
 316
 317         // Apply selection according to TPC efficiency
 318         //if(TMath::Abs(pdg)==211) nAcc++;
 319         if(fSelAndSmear && !SelectedTrack(hPt,hEta)) continue;
 320         //if(TMath::Abs(pdg)==211) nSel++;
 321
 322         // create AliTPCtrack object
 323         BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge);
 324
 325         if(fSelAndSmear) {
 326           bin = fDBgrid->GetBin(hPt,hEta);
 327           switch (pdg) {
 328           case 211:
 329             fCovTree = covTreePi[bin];
 330             break;
 331           case 321:
 332             fCovTree = covTreeKa[bin];
 333             break;
 334           case 2212:
 335             fCovTree = covTreePi[bin];
 336             break;
 337           case 11:
 338             fCovTree = covTreeEl[bin];
 339             break;
 340           case 13:
 341             fCovTree = covTreePi[bin];
 342             break;
 343           }
 344           // deal with covariance matrix and smearing of parameters
 345           CookTrack(hPt,hEta);
 346
 347           // assign the track a dE/dx and make a rough PID
 348           CookdEdx(hPt,hEta);
 349         }
 350
 351         // put track in array
 352         AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
 353         iotrack->SetLabel(label);
 354         tarray.AddLast(iotrack);
 355         // Mark track as "done" and register the pdg code
 356         done[label]=kTRUE;
 357         tracks++;
 358       }
 359
 360     } // loop over entries in TreeH
 361
 362     // sort array with TPC tracks (decreasing pT)
 363     tarray.Sort();
 364
 365     arrentr = tarray.GetEntriesFast();
 366     for(Int_t l=0; l<arrentr; l++) {
 367       tpctrack=(AliTPCtrack*)tarray.UncheckedAt(l);
 368       tracktree->Fill();
 369     }
 370
 371     // write the tree with tracks in the output file
 372     out->cd();
 373     tracktree->Write();
 374
 375     delete tracktree;
 376     delete [] done;
 377     delete [] pdgCodes;
 378
 379     printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
 380     //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
 381
 382   } // loop on events
 383
 384   if(fileDB) fileDB->Close();
 385
 386   return 0;
 387 }
 388 //-----------------------------------------------------------------------------
 389 void AliTPCtrackerParam::AnalyzedEdx(const Char_t *outName,Int_t pdg) {
 390 //-----------------------------------------------------------------------------
 391 // This function computes the dE/dx for pions, kaons, protons and electrons,
 392 // in the [pT,eta] bins.
 393 // Input file is CovMatrix_AllEvts.root.
 394 //-----------------------------------------------------------------------------
 395
 396   gStyle->SetOptStat(0);
 397   gStyle->SetOptFit(10001);
 398
 399   Char_t *part="PIONS";
 400   Double_t ymax=500.;
 401
 402   // create a chain with compared tracks
 403   TChain cmptrkchain("cmptrktree");
 404   cmptrkchain.Add("CovMatrix_AllEvts_1.root");
 405   cmptrkchain.Add("CovMatrix_AllEvts_2.root");
 406   cmptrkchain.Add("CovMatrix_AllEvts_3.root");
 407
 408   COMPTRACK cmptrk;
 409   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
 410   Int_t entries = (Int_t)cmptrkchain.GetEntries();
 411   cerr<<" Number of entries: "<<entries<<endl;
 412
 413   InitializeKineGrid("DB","points");
 414   InitializeKineGrid("dEdx","bins");
 415
 416   switch(pdg) {
 417   case 211:
 418     part = "PIONS";
 419     ymax = 100.;
 420     break;
 421   case 321:
 422     part = "KAONS";
 423     ymax = 300.;
 424     break;
 425   case 2212:
 426     part = "PROTONS";
 427     ymax = 500.;
 428     break;
 429   case 11:
 430     part = "ELECTRONS";
 431     ymax = 100.;
 432     break;
 433   }
 434
 435   SetParticle(pdg);
 436
 437   const Int_t nTotBins = fDBgrid->GetTotBins();
 438
 439   cerr<<" Fit bins: "<<nTotBins<<endl;
 440
 441   Int_t bin=0;
 442   Int_t * n = new Int_t[nTotBins];
 443   Double_t * p = new Double_t[nTotBins];
 444   Double_t * ep = new Double_t[nTotBins];
 445   Double_t * mean = new Double_t[nTotBins];
 446   Double_t * sigma = new Double_t[nTotBins];
 447
 448   for(Int_t l=0; l<nTotBins; l++) {
 449     n[l] = 1; // set to 1 to avoid divisions by 0
 450     p[l] = mean[l] = sigma[l] = ep[l] = 0.;
 451   }
 452
 453   // loop on chain entries for the mean
 454   for(Int_t l=0; l<entries; l++) {
 455     cmptrkchain.GetEvent(l);
 456     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
 457     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
 458     p[bin] += cmptrk.p;
 459     mean[bin] += cmptrk.dEdx;
 460     n[bin]++;
 461   } // loop on chain entries
 462
 463   for(Int_t l=0; l<nTotBins; l++) {
 464     p[l] /= n[l];
 465     mean[l] /= n[l];
 466     n[l] = 1; // set to 1 to avoid divisions by 0
 467   }
 468
 469   // loop on chain entries for the sigma
 470   for(Int_t l=0; l<entries; l++) {
 471     cmptrkchain.GetEvent(l);
 472     if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
 473     bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
 474     if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
 475     n[bin]++;
 476     sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
 477   } // loop on chain entries
 478
 479   for(Int_t l=0; l<nTotBins; l++) {
 480     sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
 481   }
 482
 483   // create the canvas
 484   TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700);
 485
 486   // create the graph for dEdx vs p
 487   TGraphErrors *gr = new TGraphErrors(nTotBins,p,mean,ep,sigma);
 488   TString title("  : dE/dx vs momentum"); title.Prepend(part);
 489   TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
 490   frame->SetXTitle("p [GeV/c]");
 491   frame->SetYTitle("dE/dx [a.u.]");
 492   canv->SetLogx();
 493   frame->Draw();
 494   gr->Draw("P");
 495
 496   switch(pdg) {
 497   case 211:
 498     for(Int_t i=0; i<nTotBins; i++) {
 499       fdEdxMeanPi.SetParam(i,mean[i]);
 500       fdEdxRMSPi.SetParam(i,sigma[i]);
 501     }
 502     break;
 503   case 321:
 504     for(Int_t i=0; i<nTotBins; i++) {
 505       fdEdxMeanKa.SetParam(i,mean[i]);
 506       fdEdxRMSKa.SetParam(i,sigma[i]);
 507     }
 508     break;
 509   case 2212:
 510     for(Int_t i=0; i<nTotBins; i++) {
 511       fdEdxMeanPr.SetParam(i,mean[i]);
 512       fdEdxRMSPr.SetParam(i,sigma[i]);
 513     }
 514     break;
 515   case 11:
 516     for(Int_t i=0; i<nTotBins; i++) {
 517       fdEdxMeanEl.SetParam(i,mean[i]);
 518       fdEdxRMSEl.SetParam(i,sigma[i]);
 519     }
 520     break;
 521   }
 522
 523   // write results to file
 524   WritedEdx(outName,pdg);
 525
 526   delete [] n;
 527   delete [] p;
 528   delete [] ep;
 529   delete [] mean;
 530   delete [] sigma;
 531
 532   return;
 533 }
 534 //-----------------------------------------------------------------------------
 535 void AliTPCtrackerParam::AnalyzePulls(const Char_t *outName) {
 536 //-----------------------------------------------------------------------------
 537 // This function computes the pulls for pions, kaons and electrons,
 538 // in the [pT,eta] bins.
 539 // Input file is CovMatrix_AllEvts.root.
 540 // Output file is pulls.root.
 541 //-----------------------------------------------------------------------------
 542
 543   // create a chain with compared tracks
 544   TChain cmptrkchain("cmptrktree");
 545   cmptrkchain.Add("CovMatrix_AllEvts_1.root");
 546   cmptrkchain.Add("CovMatrix_AllEvts_2.root");
 547   cmptrkchain.Add("CovMatrix_AllEvts_3.root");
 548
 549   COMPTRACK cmptrk;
 550   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
 551   Int_t entries = (Int_t)cmptrkchain.GetEntries();
 552   cerr<<" Number of entries: "<<entries<<endl;
 553
 554   Int_t thisPdg=0;
 555   Char_t hname[100], htitle[100];
 556   Int_t nTotBins,bin;
 557
 558   AliTPCkineGrid  pulls[5];
 559   TH1F *hDum = new TH1F("name","title",100,-7.,7.);
 560   TF1 *g = new TF1("g","gaus");
 561
 562   InitializeKineGrid("pulls","bins");
 563   InitializeKineGrid("DB","points");
 564
 565
 566
 567   // loop on the particles Pi,Ka,El
 568   for(Int_t part=0; part<3; part++) {
 569
 570     switch (part) {
 571     case 0: // pions
 572       thisPdg=211;
 573       cerr<<" Processing pions ...\n";
 574       break;
 575     case 1: // kaons
 576       thisPdg=321;
 577       cerr<<" Processing kaons ...\n";
 578       break;
 579
 580     case 2: // electrons
 581       thisPdg=11;
 582       cerr<<" Processing electrons ...\n";
 583       break;
 584     }
 585
 586     SetParticle(thisPdg);
 587
 588     for(Int_t i=0;i<5;i++) {
 589       pulls[i].~AliTPCkineGrid();
 590       new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
 591     }
 592     nTotBins = fDBgrid->GetTotBins();
 593
 594     // create histograms for the all the bins
 595     TH1F *hPulls0_=NULL;
 596     TH1F *hPulls1_=NULL;
 597     TH1F *hPulls2_=NULL;
 598     TH1F *hPulls3_=NULL;
 599     TH1F *hPulls4_=NULL;
 600
 601     hPulls0_ = new TH1F[nTotBins];
 602     hPulls1_ = new TH1F[nTotBins];
 603     hPulls2_ = new TH1F[nTotBins];
 604     hPulls3_ = new TH1F[nTotBins];
 605     hPulls4_ = new TH1F[nTotBins];
 606
 607
 608     for(Int_t i=0; i<nTotBins; i++) {
 609       sprintf(hname,"hPulls0_%d",i);
 610       sprintf(htitle,"P0 pulls for bin %d",i);
 611       hDum->SetName(hname); hDum->SetTitle(htitle);
 612       hPulls0_[i] = *hDum;
 613       sprintf(hname,"hPulls1_%d",i);
 614       sprintf(htitle,"P1 pulls for bin %d",i);
 615       hDum->SetName(hname); hDum->SetTitle(htitle);
 616       hPulls1_[i] = *hDum;
 617       sprintf(hname,"hPulls2_%d",i);
 618       sprintf(htitle,"P2 pulls for bin %d",i);
 619       hDum->SetName(hname); hDum->SetTitle(htitle);
 620       hPulls2_[i] = *hDum;
 621       sprintf(hname,"hPulls3_%d",i);
 622       sprintf(htitle,"P3 pulls for bin %d",i);
 623       hDum->SetName(hname); hDum->SetTitle(htitle);
 624       hPulls3_[i] = *hDum;
 625       sprintf(hname,"hPulls4_%d",i);
 626       sprintf(htitle,"P4 pulls for bin %d",i);
 627       hDum->SetName(hname); hDum->SetTitle(htitle);
 628       hPulls4_[i] = *hDum;
 629     }
 630
 631     // loop on chain entries
 632     for(Int_t i=0; i<entries; i++) {
 633       cmptrkchain.GetEvent(i);
 634       if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
 635       // fill histograms with the pulls
 636       bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
 637       hPulls0_[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
 638       hPulls1_[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
 639       hPulls2_[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
 640       hPulls3_[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
 641       hPulls4_[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
 642     } // loop on chain entries
 643
 644     // compute the sigma of the distributions
 645     for(Int_t i=0; i<nTotBins; i++) {
 646       if(hPulls0_[i].GetEntries()>1000) {
 647         g->SetRange(-3.*hPulls0_[i].GetRMS(),3.*hPulls0_[i].GetRMS());
 648         hPulls0_[i].Fit("g","R,Q,N");
 649         pulls[0].SetParam(i,g->GetParameter(2));
 650       } else pulls[0].SetParam(i,-1.);
 651       if(hPulls1_[i].GetEntries()>1000) {
 652         g->SetRange(-3.*hPulls1_[i].GetRMS(),3.*hPulls1_[i].GetRMS());
 653         hPulls1_[i].Fit("g","R,Q,N");
 654         pulls[1].SetParam(i,g->GetParameter(2));
 655       } else pulls[1].SetParam(i,-1.);
 656       if(hPulls2_[i].GetEntries()>1000) {
 657         g->SetRange(-3.*hPulls2_[i].GetRMS(),3.*hPulls2_[i].GetRMS());
 658         hPulls2_[i].Fit("g","R,Q,N");
 659         pulls[2].SetParam(i,g->GetParameter(2));
 660       } else pulls[2].SetParam(i,-1.);
 661       if(hPulls3_[i].GetEntries()>1000) {
 662         g->SetRange(-3.*hPulls3_[i].GetRMS(),3.*hPulls3_[i].GetRMS());
 663         hPulls3_[i].Fit("g","R,Q,N");
 664         pulls[3].SetParam(i,g->GetParameter(2));
 665       } else pulls[3].SetParam(i,-1.);
 666       if(hPulls4_[i].GetEntries()>1000) {
 667         g->SetRange(-3.*hPulls4_[i].GetRMS(),3.*hPulls4_[i].GetRMS());
 668         hPulls4_[i].Fit("g","R,Q,N");
 669         pulls[4].SetParam(i,g->GetParameter(2));
 670       } else pulls[4].SetParam(i,-1.);
 671     } // loop on bins
 672
 673
 674     switch (part) {
 675     case 0: // pions
 676       for(Int_t i=0;i<5;i++) {
 677         fPullsPi[i].~AliTPCkineGrid();
 678         new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]);
 679       }
 680       break;
 681     case 1: // kaons
 682       for(Int_t i=0;i<5;i++) {
 683         fPullsKa[i].~AliTPCkineGrid();
 684         new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
 685       }
 686       break;
 687     case 2: // electrons
 688       for(Int_t i=0;i<5;i++) {
 689         fPullsEl[i].~AliTPCkineGrid();
 690         new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
 691       }
 692       break;
 693     }
 694
 695     delete [] hPulls0_;
 696     delete [] hPulls1_;
 697     delete [] hPulls2_;
 698     delete [] hPulls3_;
 699     delete [] hPulls4_;
 700
 701   } // loop on particle species
 702
 703   // write pulls to file
 704   WritePulls(outName);
 705
 706
 707   return;
 708 }
 709 //-----------------------------------------------------------------------------
 710 void AliTPCtrackerParam::BuildTrack(Double_t alpha,
 711                                     Double_t x,Double_t y,Double_t z,
 712                                     Double_t px,Double_t py,
 713                                     Double_t pz,Double_t pt,
 714                                     Int_t ch) {
 715 //-----------------------------------------------------------------------------
 716 // This function uses GEANT info to set true track parameters
 717 //-----------------------------------------------------------------------------
 718   Double_t xref = x;
 719   Double_t xx[5],cc[15];
 720   cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
 721   cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
 722
 723   // Magnetic field
 724   TVector3 bfield(0.,0.,fBz);
 725
 726
 727   // radius [cm] of track projection in (x,y)
 728   Double_t rho = pt*100./0.299792458/bfield.Z();
 729   // center of track projection in local reference frame
 730   TVector3 hmom,hpos;
 731
 732
 733   // position (local) and momentum (local) at the hit
 734   // in the bending plane (z=0)
 735   hpos.SetXYZ(x,y,0.);
 736   hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.);
 737   TVector3 vrho = hmom.Cross(bfield);
 738   vrho *= ch;
 739   vrho.SetMag(rho);
 740
 741   TVector3 vcenter = hpos+vrho;
 742
 743   Double_t x0 = vcenter.X();
 744
 745   // fX     = xref         X-coordinate of this track (reference plane)
 746   // fAlpha = Alpha        Rotation angle the local (TPC sector)
 747   // fP0    = YL           Y-coordinate of a track
 748   // fP1    = ZG           Z-coordinate of a track
 749   // fP2    = C*x0         x0 is center x in rotated frame
 750   // fP3    = Tgl          tangent of the track momentum dip angle
 751   // fP4    = C            track curvature
 752   xx[0] = y;
 753   xx[1] = z;
 754   xx[3] = pz/pt;
 755   xx[4] = -ch/rho;
 756   xx[2] = xx[4]*x0;
 757
 758   // create the object AliTPCtrack
 759   AliTPCtrack track(0,xx,cc,xref,alpha);
 760   new(&fTrack) AliTPCtrack(track);
 761
 762   return;
 763 }
 764 //-----------------------------------------------------------------------------
 765 void AliTPCtrackerParam::CompareTPCtracks(
 766                            const Char_t* galiceName,
 767                            const Char_t* trkGeaName,
 768                            const Char_t* trkKalName,
 769                            const Char_t* covmatName,
 770                            const Char_t* tpceffName) const {
 771 //-----------------------------------------------------------------------------
 772 // This function compares tracks from TPC Kalman Filter V2 with
 773 // true tracks at TPC 1st hit. It gives:
 774 //   - a tree with Kalman cov. matrix elements, resolutions, dEdx
 775 //   - the efficiencies as a function of particle type, pT, eta
 776 //-----------------------------------------------------------------------------
 777
 778   TFile *kalFile    = TFile::Open(trkKalName);
 779   TFile *geaFile    = TFile::Open(trkGeaName);
 780   TFile *galiceFile = TFile::Open(galiceName);
 781
 782   // particles from TreeK
 783   AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
 784   const Int_t nparticles = gAlice->GetEvent(0);
 785
 786   Int_t * kalLab = new Int_t[nparticles];
 787   for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1;
 788
 789
 790   // tracks from Kalman
 791   TTree *kaltree=(TTree*)kalFile->Get("TreeT_TPC_0");
 792   AliTPCtrack *kaltrack=new AliTPCtrack;
 793   kaltree->SetBranchAddress("tracks",&kaltrack);
 794   Int_t kalEntries = (Int_t)kaltree->GetEntries();
 795
 796   // tracks from 1st hit
 797   TTree *geatree=(TTree*)geaFile->Get("TreeT_TPC_0");
 798   AliTPCtrack *geatrack=new AliTPCtrack;
 799   geatree->SetBranchAddress("tracks",&geatrack);
 800   Int_t geaEntries = (Int_t)geatree->GetEntries();
 801
 802   cerr<<"+++\n+++ Number of tracks:  TPC Kalman  = "<<kalEntries<<endl<<"+++                    TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
 803
 804   // set pointers for TPC tracks:
 805   // the entry number of the track labelled l is stored in kalLab[l]
 806   Int_t fake=0,mult=0;
 807   for (Int_t j=0; j<kalEntries; j++) {
 808     kaltree->GetEvent(j);
 809     if(kaltrack->GetLabel()>=0) {
 810       if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
 811       kalLab[kaltrack->GetLabel()] = j;
 812     }
 813     else fake++;
 814   }
 815
 816   cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
 817   cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
 818
 819   TParticle *Part;
 820   Int_t    label;
 821   Double_t cc[15],dAlpha;
 822   //                 ---  [Pt,eta] binning ---
 823   //
 824   //               |eta|<0.3    0.3<=|eta|<0.6    0.6<=|eta|
 825   //      Pt<0.3       0              1                2
 826   // 0.3<=Pt<0.5       3              4                5
 827   // 0.5<=Pt<1.0       6              7                8
 828   // 1.0<=Pt<1.5       9             10               11
 829   // 1.5<=Pt<3.0      12             13               14
 830   // 3.0<=Pt<5.0      15             16               17
 831   // 5.0<=Pt<7.0      18             19               20
 832   // 7.0<=Pt<15.      21             22               23
 833   // 15.<=Pt          24             25               26
 834   //
 835   Int_t    pi=0,ka=0,mu=0,el=0,pr=0;
 836   Int_t    geaPi[27],geaKa[27],geaPr[27],geaEl[27],geaMu[27];
 837   Int_t    kalPi[27],kalKa[27],kalPr[27],kalEl[27],kalMu[27];
 838   Float_t  effPi[27],effKa[27],effPr[27],effEl[27],effMu[27];
 839   Int_t bin;
 840   for(Int_t j=0; j<27; j++) {
 841     geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
 842     kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
 843     effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
 844   }
 845
 846   COMPTRACK cmptrk;
 847   // create the tree for comparison results
 848   TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
 849   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");
 850
 851   cerr<<"Doing track comparison...\n";
 852   // loop on tracks at 1st hit
 853   for(Int_t j=0; j<geaEntries; j++) {
 854     geatree->GetEvent(j);
 855
 856     label = geatrack->GetLabel();
 857     Part = (TParticle*)gAlice->Particle(label);
 858     cmptrk.pdg = Part->GetPdgCode();
 859     cmptrk.eta = Part->Eta();
 860     cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
 861
 862     cmptrk.pt   = 1/TMath::Abs(geatrack->Get1Pt());
 863     cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
 864     cmptrk.p    = cmptrk.pt/cmptrk.cosl;
 865
 866
 867     bin = GetBin(cmptrk.pt,cmptrk.eta);
 868     cmptrk.bin = bin;
 869     if(abs(cmptrk.pdg)==211)  geaPi[bin]++;
 870     if(abs(cmptrk.pdg)==321)  geaKa[bin]++;
 871     if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
 872     if(abs(cmptrk.pdg)==11)   geaEl[bin]++;
 873     if(abs(cmptrk.pdg)==13)   geaMu[bin]++;
 874
 875
 876     // check if there is the corresponding track in the TPC kalman and get it
 877     if(kalLab[label]==-1) continue;
 878     kaltree->GetEvent(kalLab[label]);
 879
 880     // and go on only if it has xref = 84.57 cm (inner pad row)
 881     if(kaltrack->GetX()>90.) continue;
 882
 883     if(abs(cmptrk.pdg)==211)  { kalPi[bin]++; pi++; }
 884     if(abs(cmptrk.pdg)==321)  { kalKa[bin]++; ka++; }
 885     if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
 886     if(abs(cmptrk.pdg)==11)   { kalEl[bin]++; el++; }
 887     if(abs(cmptrk.pdg)==13)   { kalMu[bin]++; mu++; }
 888
 889     kaltrack->PropagateTo(geatrack->GetX());
 890
 891     cmptrk.dEdx = kaltrack->GetdEdx();
 892
 893     // compute errors on parameters
 894     dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
 895     if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
 896
 897     cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
 898     cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
 899     cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
 900     cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
 901     cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
 902     cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
 903
 904     // get covariance matrix
 905     // beware: lines 3 and 4 are inverted!
 906     kaltrack->GetCovariance(cc);
 907
 908     cmptrk.c00 = cc[0];
 909     cmptrk.c10 = cc[1];
 910     cmptrk.c11 = cc[2];
 911     cmptrk.c20 = cc[3];
 912     cmptrk.c21 = cc[4];
 913     cmptrk.c22 = cc[5];
 914     cmptrk.c30 = cc[10];
 915     cmptrk.c31 = cc[11];
 916     cmptrk.c32 = cc[12];
 917     cmptrk.c33 = cc[14];
 918     cmptrk.c40 = cc[6];
 919     cmptrk.c41 = cc[7];
 920     cmptrk.c42 = cc[8];
 921     cmptrk.c43 = cc[13];
 922     cmptrk.c44 = cc[9];
 923
 924     // fill tree
 925     cmptrktree->Fill();
 926
 927   } // loop on tracks at TPC 1st hit
 928
 929   cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
 930   cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
 931   cerr<<"+++"<<endl;
 932
 933   // Write tree to file
 934   TFile *outfile = new TFile(covmatName,"recreate");
 935   cmptrktree->Write();
 936   outfile->Close();
 937
 938
 939   // Write efficiencies to file
 940   FILE *effFile = fopen(tpceffName,"w");
 941   //fprintf(effFile,"%d\n",kalEntries);
 942   for(Int_t j=0; j<27; j++) {
 943     if(geaPi[j]>=5) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
 944     if(geaKa[j]>=5) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
 945     if(geaPr[j]>=5) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
 946     if(geaEl[j]>=5) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
 947     if(geaMu[j]>=5) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
 948     fprintf(effFile,"%f  %f  %f  %f  %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
 949   }
 950
 951   for(Int_t j=0; j<27; j++) {
 952     fprintf(effFile,"%d  %d  %d  %d  %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
 953   }
 954   for(Int_t j=0; j<27; j++) {
 955     fprintf(effFile,"%d  %d  %d  %d  %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
 956   }
 957   fclose(effFile);
 958
 959   // delete AliRun object
 960   delete gAlice; gAlice=0;
 961
 962   // close all input files
 963   kalFile->Close();
 964   geaFile->Close();
 965   galiceFile->Close();
 966
 967   delete [] kalLab;
 968
 969   return;
 970 }
 971 //-----------------------------------------------------------------------------
 972 void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) {
 973 //-----------------------------------------------------------------------------
 974 // This function assigns the track a dE/dx and makes a rough PID
 975 //-----------------------------------------------------------------------------
 976
 977   Double_t mean = fdEdxMean->GetValueAt(pt,eta);
 978   Double_t rms  = fdEdxRMS->GetValueAt(pt,eta);
 979
 980   Double_t dEdx = gRandom->Gaus(mean,rms);
 981
 982   fTrack.SetdEdx(dEdx);
 983
 984   AliTPCtrackParam t(fTrack);
 985
 986   //Very rough PID
 987   Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
 988
 989   if (p<0.6) {
 990     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
 991       t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
 992     }
 993     if (dEdx < 39.+ 12./p/p) {
 994       t.AssignMass(0.49368); new(&fTrack) AliTPCtrack(t); return;
 995     }
 996     t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
 997   }
 998
 999   if (p<1.2) {
1000     if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
1001       t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1002     }
1003     t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return;
1004   }
1005
1006   t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return;
1007 }
1008 //-----------------------------------------------------------------------------
1009 void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
1010 //-----------------------------------------------------------------------------
1011 // This function deals with covariance matrix and smearing
1012 //-----------------------------------------------------------------------------
1013
1014   COVMATRIX covmat;
1015   Double_t  p,cosl;
1016   Double_t  trkKine[2],trkRegPar[4];
1017   Double_t  xref,alpha,xx[5],xxsm[5],cc[15];
1018   Int_t     treeEntries;
1019
1020   fCovTree->SetBranchAddress("matrix",&covmat);
1021
1022   // get random entry from the tree
1023   treeEntries = (Int_t)fCovTree->GetEntries();
1024   fCovTree->GetEvent(gRandom->Integer(treeEntries));
1025
1026   // get P and Cosl from track
1027   cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl()));
1028   p    = 1./TMath::Abs(fTrack.Get1Pt())/cosl;
1029
1030   trkKine[0] = p;
1031   trkKine[1] = cosl;
1032
1033   // get covariance matrix from regularized matrix
1034   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(0,l);
1035   cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
1036   cc[1] = covmat.c10;
1037   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(1,l);
1038   cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
1039   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(2,l);
1040   cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
1041   cc[4] = covmat.c21;
1042   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(3,l);
1043   cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
1044   cc[6] = covmat.c30;
1045   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(4,l);
1046   cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
1047   cc[8] = covmat.c32;
1048   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(5,l);
1049   cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
1050   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(6,l);
1051   cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
1052   cc[11]= covmat.c41;
1053   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(7,l);
1054   cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
1055   cc[13]= covmat.c43;
1056   for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(8,l);
1057   cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
1058
1059   TMatrixD covMatSmear(5,5);
1060
1061   covMatSmear = GetSmearingMatrix(cc,pt,eta);
1062
1063   // get track original parameters
1064   xref =fTrack.GetX();
1065   alpha=fTrack.GetAlpha();
1066   xx[0]=fTrack.GetY();
1067   xx[1]=fTrack.GetZ();
1068   xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp();
1069   xx[3]=fTrack.GetTgl();
1070   xx[4]=fTrack.GetC();
1071
1072   // use smearing matrix to smear the original parameters
1073   SmearTrack(xx,xxsm,covMatSmear);
1074
1075   AliTPCtrack track(0,xxsm,cc,xref,alpha);
1076   new(&fTrack) AliTPCtrack(track);
1077
1078   return;
1079 }
1080 //-----------------------------------------------------------------------------
1081 void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg) {
1082 //-----------------------------------------------------------------------------
1083 // This function draws the TPC efficiencies in the [pT,eta] bins
1084 //-----------------------------------------------------------------------------
1085
1086   ReadEffs(inName);
1087   SetParticle(pdg);
1088
1089   const Int_t n = fEff->GetPointsPt();
1090   Double_t * effsA = new Double_t[n];
1091   Double_t * effsB = new Double_t[n];
1092   Double_t * effsC = new Double_t[n];
1093   Double_t * pt = new Double_t[n];
1094
1095   fEff->GetArrayPt(pt);
1096   for(Int_t i=0;i<n;i++) {
1097     effsA[i] = fEff->GetParam(i,0);
1098     effsB[i] = fEff->GetParam(i,1);
1099     effsC[i] = fEff->GetParam(i,2);
1100   }
1101
1102   TGraph *grA = new TGraph(n,pt,effsA);
1103   TGraph *grB = new TGraph(n,pt,effsB);
1104   TGraph *grC = new TGraph(n,pt,effsC);
1105
1106   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1107   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1108   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1109
1110   TString title("Distribution of the TPC efficiencies");
1111   switch (pdg) {
1112   case 211:
1113     title.Prepend("PIONS - ");
1114     break;
1115   case 321:
1116     title.Prepend("KAONS - ");
1117     break;
1118   case 2212:
1119     title.Prepend("PROTONS - ");
1120     break;
1121   case 11:
1122     title.Prepend("ELECTRONS - ");
1123     break;
1124   case 13:
1125     title.Prepend("MUONS - ");
1126     break;
1127   }
1128
1129
1130   gStyle->SetOptStat(0);
1131   TCanvas *c = new TCanvas("c","effs",0,0,900,700);
1132   c->SetLogx();
1133   c->SetGridx();
1134   c->SetGridy();
1135
1136   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1);
1137   frame->SetXTitle("p_{T} [GeV/c]");
1138   frame->Draw();
1139   grA->Draw("P");
1140   grB->Draw("P");
1141   grC->Draw("P");
1142
1143   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1144   leg->AddEntry(grA,"|#eta|<0.3","p");
1145   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1146   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1147   leg->SetFillColor(0);
1148   leg->Draw();
1149
1150   delete [] effsA;
1151   delete [] effsB;
1152   delete [] effsC;
1153   delete [] pt;
1154
1155   return;
1156 }
1157 //-----------------------------------------------------------------------------
1158 void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg,
1159                                    Int_t par) {
1160 //-----------------------------------------------------------------------------
1161 // This function draws the pulls in the [pT,eta] bins
1162 //-----------------------------------------------------------------------------
1163
1164   ReadPulls(inName);
1165   SetParticle(pdg);
1166
1167   const Int_t n = (fPulls+par)->GetPointsPt();
1168   Double_t * pullsA = new Double_t[n];
1169   Double_t * pullsB = new Double_t[n];
1170   Double_t * pullsC = new Double_t[n];
1171   Double_t * pt = new Double_t[n];
1172   (fPulls+par)->GetArrayPt(pt);
1173   for(Int_t i=0;i<n;i++) {
1174     pullsA[i] = (fPulls+par)->GetParam(i,0);
1175     pullsB[i] = (fPulls+par)->GetParam(i,1);
1176     pullsC[i] = (fPulls+par)->GetParam(i,2);
1177   }
1178
1179   TGraph *grA = new TGraph(n,pt,pullsA);
1180   TGraph *grB = new TGraph(n,pt,pullsB);
1181   TGraph *grC = new TGraph(n,pt,pullsC);
1182
1183   grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5);
1184   grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5);
1185   grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5);
1186
1187   TString title("Distribution of the pulls: ");
1188   switch (pdg) {
1189   case 211:
1190     title.Prepend("PIONS - ");
1191     break;
1192   case 321:
1193     title.Prepend("KAONS - ");
1194     break;
1195   case 11:
1196     title.Prepend("ELECTRONS - ");
1197     break;
1198   }
1199   switch (par) {
1200   case 0:
1201     title.Append("y");
1202     break;
1203   case 1:
1204     title.Append("z");
1205     break;
1206   case 2:
1207     title.Append("    #eta");
1208     break;
1209   case 3:
1210     title.Append("tg    #lambda");
1211     break;
1212   case 4:
1213     title.Append("C");
1214     break;
1215   }
1216
1217   gStyle->SetOptStat(0);
1218   TCanvas *c = new TCanvas("c","pulls",0,0,900,700);
1219   c->SetLogx();
1220   c->SetGridx();
1221   c->SetGridy();
1222
1223   TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2);
1224   frame->SetXTitle("p_{T} [GeV/c]");
1225   frame->Draw();
1226   grA->Draw("P");
1227   grB->Draw("P");
1228   grC->Draw("P");
1229
1230   TLegend *leg = new TLegend(0.2,0.2,0.4,0.4);
1231   leg->AddEntry(grA,"|#eta|<0.3","p");
1232   leg->AddEntry(grB,"0.3<|#eta|<0.6","p");
1233   leg->AddEntry(grC,"0.6<|#eta|<0.9","p");
1234   leg->SetFillColor(0);
1235   leg->Draw();
1236
1237   delete [] pullsA;
1238   delete [] pullsB;
1239   delete [] pullsC;
1240   delete [] pt;
1241
1242   return;
1243 }
1244 //-----------------------------------------------------------------------------
1245 Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
1246 //-----------------------------------------------------------------------------
1247 // This function tells bin number in a grid (pT,eta)
1248 //-----------------------------------------------------------------------------
1249   if(TMath::Abs(eta)<0.3) {
1250     if(pt<0.3)            return 0;
1251     if(pt>=0.3 && pt<0.5) return 3;
1252     if(pt>=0.5 && pt<1.)  return 6;
1253     if(pt>=1. && pt<1.5)  return 9;
1254     if(pt>=1.5 && pt<3.)  return 12;
1255     if(pt>=3. && pt<5.)   return 15;
1256     if(pt>=5. && pt<7.)   return 18;
1257     if(pt>=7. && pt<15.)  return 21;
1258     if(pt>=15.)           return 24;
1259   }
1260   if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
1261     if(pt<0.3)            return 1;
1262     if(pt>=0.3 && pt<0.5) return 4;
1263     if(pt>=0.5 && pt<1.)  return 7;
1264     if(pt>=1. && pt<1.5)  return 10;
1265     if(pt>=1.5 && pt<3.)  return 13;
1266     if(pt>=3. && pt<5.)   return 16;
1267     if(pt>=5. && pt<7.)   return 19;
1268     if(pt>=7. && pt<15.)  return 22;
1269     if(pt>=15.)           return 25;
1270   }
1271   if(TMath::Abs(eta)>=0.6) {
1272     if(pt<0.3)            return 2;
1273     if(pt>=0.3 && pt<0.5) return 5;
1274     if(pt>=0.5 && pt<1.)  return 8;
1275     if(pt>=1. && pt<1.5)  return 11;
1276     if(pt>=1.5 && pt<3.)  return 14;
1277     if(pt>=3. && pt<5.)   return 17;
1278     if(pt>=5. && pt<7.)   return 20;
1279     if(pt>=7. && pt<15.)  return 23;
1280     if(pt>=15.)           return 26;
1281   }
1282
1283   return -1;
1284
1285 }
1286 //-----------------------------------------------------------------------------
1287 TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
1288                                                Double_t eta) const {
1289 //-----------------------------------------------------------------------------
1290 // This function stretches the covariance matrix according to the pulls
1291 //-----------------------------------------------------------------------------
1292
1293   TMatrixD covMat(5,5);
1294
1295   covMat(0,0)=cc[0];
1296   covMat(1,0)=cc[1];  covMat(0,1)=covMat(1,0);
1297   covMat(1,1)=cc[2];
1298   covMat(2,0)=cc[3];  covMat(0,2)=covMat(2,0);
1299   covMat(2,1)=cc[4];  covMat(1,2)=covMat(2,1);
1300   covMat(2,2)=cc[5];
1301   covMat(3,0)=cc[6];  covMat(0,3)=covMat(3,0);
1302   covMat(3,1)=cc[7];  covMat(1,3)=covMat(3,1);
1303   covMat(3,2)=cc[8];  covMat(2,3)=covMat(3,2);
1304   covMat(3,3)=cc[9];
1305   covMat(4,0)=cc[10]; covMat(0,4)=covMat(4,0);
1306   covMat(4,1)=cc[11]; covMat(1,4)=covMat(4,1);
1307   covMat(4,2)=cc[12]; covMat(2,4)=covMat(4,2);
1308   covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
1309   covMat(4,4)=cc[14];
1310
1311   TMatrixD stretchMat(5,5);
1312   for(Int_t k=0;k<5;k++) {
1313     for(Int_t l=0;l<5;l++) {
1314       stretchMat(k,l)=0.;
1315     }
1316   }
1317
1318   for(Int_t i=0;i<5;i++) stretchMat(i,i) =
1319                            TMath::Sqrt((fPulls+i)->GetValueAt(pt,eta));
1320
1321   TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
1322   TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
1323
1324   return covMatSmear;
1325 }
1326 //-----------------------------------------------------------------------------
1327 void AliTPCtrackerParam::InitializeKineGrid(Option_t* which,Option_t* how) {
1328 //-----------------------------------------------------------------------------
1329 // This function initializes ([pt,eta] points) the data members AliTPCkineGrid
1330 // which = "DB"     -> initialize fDBgrid... members
1331 //         "eff"    -> initialize fEff... members
1332 //         "pulls"  -> initialize fPulls... members
1333 //         "dEdx"   -> initialize fdEdx... members
1334 // how =   "points" -> initialize with the points of the grid
1335 //         "bins"   -> initialize with the bins of the grid
1336 //-----------------------------------------------------------------------------
1337
1338   const char *points = strstr(how,"points");
1339   const char *bins   = strstr(how,"bins");
1340   const char *DB     = strstr(which,"DB");
1341   const char *eff    = strstr(which,"eff");
1342   const char *pulls  = strstr(which,"pulls");
1343   const char *dEdx   = strstr(which,"dEdx");
1344
1345
1346   Int_t nEta=0,nPtPi=0,nPtKa=0,nPtPr=0,nPtEl=0,nPtMu=0;
1347
1348   Double_t etaPoints[2] = {0.3,0.6};
1349   Double_t etaBins[3]   = {0.15,0.45,0.75};
1350   Double_t ptPoints8[8] = {0.3,0.5,1.,1.5,3.,5.,7.,15.};
1351   Double_t ptPoints5[5] = {0.3,0.5,1.,1.5,5.};
1352   Double_t ptBins9[9]   = {0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
1353   Double_t ptBins6[6]   = {0.244,0.390,0.676,1.190,2.36,10.};
1354
1355   Double_t *eta=0,*ptPi=0,*ptKa=0,*ptPr=0,*ptEl=0,*ptMu=0;
1356
1357   if(points) {
1358     nEta  = 2;
1359     nPtPi = 8;
1360     nPtKa = 5;
1361     nPtPr = 4;
1362     nPtEl = 8;
1363     nPtMu = 2;
1364     eta  = etaPoints;
1365     ptPi = ptPoints8;
1366     ptKa = ptPoints5;
1367     ptPr = ptPoints8;
1368     ptEl = ptPoints8;
1369     ptMu = ptPoints8;
1370   }
1371   if(bins) {
1372     nEta  = 3;
1373     nPtPi = 9;
1374     nPtKa = 6;
1375     nPtPr = 5;
1376     nPtEl = 9;
1377     nPtMu = 3;
1378     eta  = etaBins;
1379     ptPi = ptBins9;
1380     ptKa = ptBins6;
1381     ptPr = ptBins9;
1382     ptEl = ptBins9;
1383     ptMu = ptBins9;
1384   }
1385
1386   AliTPCkineGrid *dummy=0;
1387
1388   if(DB) {
1389     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1390     new(&fDBgridPi) AliTPCkineGrid(*dummy);
1391     delete dummy;
1392     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1393     new(&fDBgridKa) AliTPCkineGrid(*dummy);
1394     delete dummy;
1395     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1396     new(&fDBgridEl) AliTPCkineGrid(*dummy);
1397     delete dummy;
1398   }
1399   if(eff) {
1400     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1401     new(&fEffPi) AliTPCkineGrid(*dummy);
1402     delete dummy;
1403     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1404     new(&fEffKa) AliTPCkineGrid(*dummy);
1405     delete dummy;
1406     dummy = new AliTPCkineGrid(nPtPr,nEta,ptPr,eta);
1407     new(&fEffPr) AliTPCkineGrid(*dummy);
1408     delete dummy;
1409     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1410     new(&fEffEl) AliTPCkineGrid(*dummy);
1411     delete dummy;
1412     dummy = new AliTPCkineGrid(nPtMu,nEta,ptMu,eta);
1413     new(&fEffMu) AliTPCkineGrid(*dummy);
1414     delete dummy;
1415   }
1416   if(pulls) {
1417     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1418     for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
1419     delete dummy;
1420     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1421     for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
1422     delete dummy;
1423     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1424     for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
1425     delete dummy;
1426   }
1427   if(dEdx) {
1428     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1429     new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
1430     new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
1431     delete dummy;
1432     dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
1433     new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
1434     new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
1435     delete dummy;
1436     dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
1437     new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
1438     new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
1439     delete dummy;
1440     dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
1441     new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
1442     new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
1443     delete dummy;
1444   }
1445
1446   return;
1447 }
1448 //-----------------------------------------------------------------------------
1449 void AliTPCtrackerParam::MakeDataBase() {
1450 //-----------------------------------------------------------------------------
1451 // This function creates the DB file and store in it:
1452 //  - TPC Efficiencies for pi,ka,pr,el,mu     (AliTPCkineGrid class)
1453 //  - Pulls for pi,ka,el                      (AliTPCkineGrid class)
1454 //  - Regularization parameters for pi,ka,el  (TMatrixD class)
1455 //  - dE/dx parameterization for pi,ka,pr,el  (AliTPCkineGrid class)
1456 //  - Regularized cov. matrices for pi,ka,el  (COVMATRIX structure)
1457 //-----------------------------------------------------------------------------
1458
1459   // define some file names
1460   Char_t *effFile  ="CovMatrixDB_partial.root";
1461   Char_t *pullsFile="CovMatrixDB_partial.root";
1462   Char_t *regPiFile="CovMatrixDB_partial.root";
1463   Char_t *regKaFile="CovMatrixDB_partial.root";
1464   Char_t *regElFile="CovMatrixDB_partial.root";
1465   Char_t *dEdxPiFile="dEdxPi.root";
1466   Char_t *dEdxKaFile="dEdxKa.root";
1467   Char_t *dEdxPrFile="dEdxPr.root";
1468   Char_t *dEdxElFile="dEdxEl.root";
1469   Char_t *cmFile1  ="CovMatrix_AllEvts_1.root";
1470   Char_t *cmFile2  ="CovMatrix_AllEvts_2.root";
1471   Char_t *cmFile3  ="CovMatrix_AllEvts_3.root";
1472
1473   // store the effieciencies
1474   ReadEffs(effFile);
1475   WriteEffs(fDBfileName.Data());
1476
1477   // store the pulls
1478   ReadPulls(pullsFile);
1479   WritePulls(fDBfileName.Data());
1480
1481   //store the regularization parameters
1482   ReadRegParams(regPiFile,211);
1483   WriteRegParams(fDBfileName.Data(),211);
1484   ReadRegParams(regKaFile,321);
1485   WriteRegParams(fDBfileName.Data(),321);
1486   ReadRegParams(regElFile,11);
1487   WriteRegParams(fDBfileName.Data(),11);
1488
1489   // store the dEdx parameters
1490   ReaddEdx(dEdxPiFile,211);
1491   WritedEdx(fDBfileName.Data(),211);
1492   ReaddEdx(dEdxKaFile,321);
1493   WritedEdx(fDBfileName.Data(),321);
1494   ReaddEdx(dEdxPrFile,2212);
1495   WritedEdx(fDBfileName.Data(),2212);
1496   ReaddEdx(dEdxElFile,11);
1497   WritedEdx(fDBfileName.Data(),11);
1498
1499
1500   //
1501   // store the regularized covariance matrices
1502   //
1503   InitializeKineGrid("DB","points");
1504
1505   const Int_t nBinsPi = fDBgridPi.GetTotBins();
1506   const Int_t nBinsKa = fDBgridKa.GetTotBins();
1507   const Int_t nBinsEl = fDBgridEl.GetTotBins();
1508
1509
1510   // create the trees for cov. matrices
1511   // trees for pions
1512   TTree *CovTreePi_ = NULL;
1513   CovTreePi_ = new TTree[nBinsPi];
1514   // trees for kaons
1515   TTree *CovTreeKa_ = NULL;
1516   CovTreeKa_ = new TTree[nBinsKa];
1517   // trees for electrons
1518   TTree *CovTreeEl_ = NULL;
1519   CovTreeEl_ = new TTree[nBinsEl];
1520
1521   Char_t hname[100], htitle[100];
1522   COVMATRIX covmat;
1523
1524
1525   for(Int_t i=0; i<nBinsPi; i++) {
1526     sprintf(hname,"CovTreePi_bin%d",i);
1527     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1528     CovTreePi_[i].SetName(hname); CovTreePi_[i].SetTitle(htitle);
1529     CovTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
1530   }
1531   for(Int_t i=0; i<nBinsKa; i++) {
1532     sprintf(hname,"CovTreeKa_bin%d",i);
1533     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1534     CovTreeKa_[i].SetName(hname); CovTreeKa_[i].SetTitle(htitle);
1535     CovTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1536   }
1537   for(Int_t i=0; i<nBinsEl; i++) {
1538     sprintf(hname,"CovTreeEl_bin%d",i);
1539     sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
1540     CovTreeEl_[i].SetName(hname); CovTreeEl_[i].SetTitle(htitle);
1541     CovTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
1542   }
1543
1544   // create the chain with the compared tracks
1545   TChain cmptrkchain("cmptrktree");
1546   cmptrkchain.Add(cmFile1);
1547   cmptrkchain.Add(cmFile2);
1548   cmptrkchain.Add(cmFile3);
1549
1550   COMPTRACK cmptrk;
1551   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
1552   Int_t entries = (Int_t)cmptrkchain.GetEntries();
1553   cerr<<" Number of entries: "<<entries<<endl;
1554
1555   Int_t trkPdg,trkBin;
1556   Double_t trkKine[2],trkRegPar[4];
1557   Int_t * nPerBinPi = new Int_t[nBinsPi];
1558   for(Int_t k=0;k<nBinsPi;k++) nPerBinPi[k]=0;
1559   Int_t * nPerBinKa = new Int_t[nBinsKa];
1560   for(Int_t k=0;k<nBinsKa;k++) nPerBinKa[k]=0;
1561   Int_t * nPerBinEl = new Int_t[nBinsEl];
1562   for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
1563
1564   // loop on chain entries
1565   for(Int_t l=0; l<entries; l++) {
1566     if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
1567     // get the event
1568     cmptrkchain.GetEvent(l);
1569     // get the pdg code
1570     trkPdg = TMath::Abs(cmptrk.pdg);
1571     // use only pions, kaons, electrons
1572     if(trkPdg!=211 && trkPdg!=321 && trkPdg!=11) continue;
1573     SetParticle(trkPdg);
1574     trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
1575     //cerr<<cmptrk.pt<<"  "<<cmptrk.eta<<"  "<<trkBin<<endl;
1576
1577     if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
1578     if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
1579     if(trkPdg==11  && nPerBinEl[trkBin]>=5000) continue;
1580
1581     trkKine[0] = cmptrk.p;
1582     trkKine[1] = cmptrk.cosl;
1583
1584     // get regularized covariance matrix
1585     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(0,k);
1586     covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
1587     covmat.c10 = cmptrk.c10;
1588     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(1,k);
1589     covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
1590     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(2,k);
1591     covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
1592     covmat.c21 = cmptrk.c21;
1593     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(3,k);
1594     covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
1595     covmat.c30 = cmptrk.c30;
1596     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(4,k);
1597     covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
1598     covmat.c32 = cmptrk.c32;
1599     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(5,k);
1600     covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
1601     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(6,k);
1602     covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
1603     covmat.c41 = cmptrk.c41;
1604     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(7,k);
1605     covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
1606     covmat.c43 = cmptrk.c43;
1607     for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(8,k);
1608     covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
1609
1610     // fill the tree
1611     switch (trkPdg) {
1612     case 211: // pions
1613       CovTreePi_[trkBin].Fill();
1614       nPerBinPi[trkBin]++;
1615       break;
1616     case 321: // kaons
1617       CovTreeKa_[trkBin].Fill();
1618       nPerBinKa[trkBin]++;
1619       break;
1620     case 11: // electrons
1621       CovTreeEl_[trkBin].Fill();
1622       nPerBinEl[trkBin]++;
1623       break;
1624     }
1625   } // loop on chain entries
1626
1627   // store all trees the DB file
1628   TFile *DBfile = new TFile(fDBfileName.Data(),"update");
1629   DBfile->mkdir("CovMatrices");
1630   gDirectory->cd("/CovMatrices");
1631   gDirectory->mkdir("Pions");
1632   gDirectory->mkdir("Kaons");
1633   gDirectory->mkdir("Electrons");
1634   // store pions
1635   gDirectory->cd("/CovMatrices/Pions");
1636   fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
1637   for(Int_t i=0;i<nBinsPi;i++) CovTreePi_[i].Write();
1638   // store kaons
1639   gDirectory->cd("/CovMatrices/Kaons");
1640   fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
1641   for(Int_t i=0;i<nBinsKa;i++) CovTreeKa_[i].Write();
1642   // store electrons
1643   gDirectory->cd("/CovMatrices/Electrons");
1644   fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
1645   for(Int_t i=0;i<nBinsEl;i++) CovTreeEl_[i].Write();
1646
1647   DBfile->Close();
1648   delete [] nPerBinPi;
1649   delete [] nPerBinKa;
1650   delete [] nPerBinEl;
1651
1652   return;
1653 }
1654 //-----------------------------------------------------------------------------
1655 void AliTPCtrackerParam::MergeEvents(Int_t evFirst,Int_t evLast) {
1656 //-----------------------------------------------------------------------------
1657 // This function: 1) merges the files from track comparison
1658 //                   (beware: better no more than 100 events per file)
1659 //                2) computes the average TPC efficiencies
1660 //-----------------------------------------------------------------------------
1661
1662   Char_t *outName="TPCeff.root";
1663
1664   Int_t    nCol;
1665   Float_t effPi,effKa,effPr,effEl,effMu;
1666   Float_t avEffPi[27],avEffKa[27],avEffPr[27],avEffEl[27],avEffMu[27];
1667   Int_t    evtsPi[27],evtsKa[27],evtsPr[27],evtsEl[27],evtsMu[27];
1668   for(Int_t j=0; j<27; j++) {
1669     avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
1670     evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
1671   }
1672
1673   // create the chain for the tree of compared tracks
1674   TChain ch1("cmptrktree");
1675   TChain ch2("cmptrktree");
1676   TChain ch3("cmptrktree");
1677
1678
1679   for(Int_t j=evFirst; j<=evLast; j++) {
1680     cerr<<"Processing event "<<j<<endl;
1681
1682     TString covName("CovMatrix.");
1683     TString effName("TPCeff.");
1684     covName+=j;
1685     effName+=j;
1686     covName.Append(".root");
1687     effName.Append(".dat");
1688
1689     if(gSystem->AccessPathName(covName.Data(),kFileExists) ||
1690        gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
1691
1692     if(j<=100) ch1.Add(covName.Data());
1693     if(j>100 && j<=200) ch2.Add(covName.Data());
1694     if(j>200) ch3.Add(covName.Data());
1695
1696
1697     FILE *effIn = fopen(effName.Data(),"r");
1698     for(Int_t k=0; k<27; k++) {
1699       nCol = fscanf(effIn,"%f  %f  %f  %f  %f",&effPi,&effKa,&effPr,&effEl,&effMu);
1700       if(effPi>=0.) {avEffPi[k]+=effPi; evtsPi[k]++;}
1701       if(effKa>=0.) {avEffKa[k]+=effKa; evtsKa[k]++;}
1702       if(effPr>=0.) {avEffPr[k]+=effPr; evtsPr[k]++;}
1703       if(effEl>=0.) {avEffEl[k]+=effEl; evtsEl[k]++;}
1704       if(effMu>=0.) {avEffMu[k]+=effMu; evtsMu[k]++;}
1705     }
1706     fclose(effIn);
1707
1708   } // loop on events
1709
1710   // merge chain in one file
1711   TFile *covOut=0;
1712   covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
1713   ch1.Merge(covOut,1000000000);
1714   covOut->Close();
1715   delete covOut;
1716   covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
1717   ch2.Merge(covOut,1000000000);
1718   covOut->Close();
1719   delete covOut;
1720   covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
1721   ch3.Merge(covOut,1000000000);
1722   covOut->Close();
1723   delete covOut;
1724
1725
1726   // efficiencies
1727
1728   InitializeKineGrid("eff","bins");
1729
1730   // pions
1731   for(Int_t j=0; j<27; j++) {
1732     if(evtsPi[j]==0) evtsPi[j]++;
1733     fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
1734   }
1735
1736   // kaons
1737   Int_t j=0;
1738   for(Int_t k=0; k<27; k++) {
1739     if(k>14 && k!=21 && k!=22 && k!=23) continue;
1740     if(evtsKa[k]==0) evtsKa[k]++;
1741     fEffKa.SetParam(j,(Double_t)avEffKa[k]/evtsKa[k]);
1742     j++;
1743   }
1744
1745   // protons
1746   for(Int_t j=0; j<15; j++) {
1747     if(evtsPr[j]==0) evtsPr[j]++;
1748     fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
1749   }
1750
1751   // electrons
1752   for(Int_t j=0; j<27; j++) {
1753     if(evtsEl[j]==0) evtsEl[j]++;
1754     fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
1755   }
1756
1757   // muons
1758   for(Int_t j=0; j<9; j++) {
1759     if(evtsMu[j]==0) evtsMu[j]++;
1760     fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
1761   }
1762
1763   // write efficiencies to a file
1764   WriteEffs(outName);
1765
1766   return;
1767 }
1768 //-----------------------------------------------------------------------------
1769 Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
1770 //-----------------------------------------------------------------------------
1771 // This function reads all parameters from the DB
1772 //-----------------------------------------------------------------------------
1773
1774   if(!ReadEffs(inName)) return 0;
1775   if(!ReadPulls(inName)) return 0;
1776   if(!ReadRegParams(inName,211)) return 0;
1777   if(!ReadRegParams(inName,321)) return 0;
1778   if(!ReadRegParams(inName,11)) return 0;
1779   if(!ReaddEdx(inName,211)) return 0;
1780   if(!ReaddEdx(inName,321)) return 0;
1781   if(!ReaddEdx(inName,2212)) return 0;
1782   if(!ReaddEdx(inName,11)) return 0;
1783   if(!ReadDBgrid(inName)) return 0;
1784
1785   return 1;
1786 }
1787 //-----------------------------------------------------------------------------
1788 Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
1789 //-----------------------------------------------------------------------------
1790 // This function reads the dEdx parameters from the DB
1791 //-----------------------------------------------------------------------------
1792
1793   if(gSystem->AccessPathName(inName,kFileExists)) {
1794     cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
1795     return 0; }
1796
1797   TFile *inFile = TFile::Open(inName);
1798   switch (pdg) {
1799   case 211:
1800     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
1801     fdEdxMeanPi.~AliTPCkineGrid();
1802     new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
1803     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
1804     fdEdxRMSPi.~AliTPCkineGrid();
1805     new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
1806     break;
1807   case 321:
1808     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
1809     fdEdxMeanKa.~AliTPCkineGrid();
1810     new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
1811     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
1812     fdEdxRMSKa.~AliTPCkineGrid();
1813     new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
1814     break;
1815   case 2212:
1816     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
1817     fdEdxMeanPr.~AliTPCkineGrid();
1818     new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
1819     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
1820     fdEdxRMSPr.~AliTPCkineGrid();
1821     new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
1822     break;
1823   case 11:
1824     fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
1825     fdEdxMeanEl.~AliTPCkineGrid();
1826     new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
1827     fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
1828     fdEdxRMSEl.~AliTPCkineGrid();
1829     new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
1830     break;
1831   }
1832   inFile->Close();
1833
1834   return 1;
1835 }
1836 //-----------------------------------------------------------------------------
1837 Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
1838 //-----------------------------------------------------------------------------
1839 // This function reads the kine grid from the DB
1840 //-----------------------------------------------------------------------------
1841
1842   if(gSystem->AccessPathName(inName,kFileExists)) {
1843     cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n";
1844     return 0; }
1845
1846   TFile *inFile = TFile::Open(inName);
1847
1848   // first read the DB grid for the different particles
1849   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
1850   fDBgridPi.~AliTPCkineGrid();
1851   new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
1852   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
1853   fDBgridKa.~AliTPCkineGrid();
1854   new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
1855   fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
1856   fDBgridEl.~AliTPCkineGrid();
1857   new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
1858
1859   inFile->Close();
1860
1861   return 1;
1862 }
1863 //-----------------------------------------------------------------------------
1864 Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
1865 //-----------------------------------------------------------------------------
1866 // This function reads the TPC efficiencies from the DB
1867 //-----------------------------------------------------------------------------
1868
1869   if(gSystem->AccessPathName(inName,kFileExists)) {
1870     cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n";
1871     return 0; }
1872
1873   TFile *inFile = TFile::Open(inName);
1874
1875   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
1876   fEffPi.~AliTPCkineGrid();
1877   new(&fEffPi) AliTPCkineGrid(*fEff);
1878   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
1879   fEffKa.~AliTPCkineGrid();
1880   new(&fEffKa) AliTPCkineGrid(*fEff);
1881   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
1882   fEffPr.~AliTPCkineGrid();
1883   new(&fEffPr) AliTPCkineGrid(*fEff);
1884   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
1885   fEffEl.~AliTPCkineGrid();
1886   new(&fEffEl) AliTPCkineGrid(*fEff);
1887   fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
1888   fEffMu.~AliTPCkineGrid();
1889   new(&fEffMu) AliTPCkineGrid(*fEff);
1890
1891   inFile->Close();
1892
1893   return 1;
1894 }
1895 //-----------------------------------------------------------------------------
1896 Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
1897 //-----------------------------------------------------------------------------
1898 // This function reads the pulls from the DB
1899 //-----------------------------------------------------------------------------
1900
1901   if(gSystem->AccessPathName(inName,kFileExists)) {
1902     cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n";
1903     return 0; }
1904
1905   TFile *inFile = TFile::Open(inName);
1906
1907   for(Int_t i=0; i<5; i++) {
1908     TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
1909     TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
1910     TString el("/Pulls/Electrons/PullsEl_"); el+=i;
1911     fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
1912     fPullsPi[i].~AliTPCkineGrid();
1913     new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
1914     fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
1915     fPullsKa[i].~AliTPCkineGrid();
1916     new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
1917     fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
1918     fPullsEl[i].~AliTPCkineGrid();
1919     new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
1920   }
1921
1922   inFile->Close();
1923
1924   return 1;
1925 }
1926 //-----------------------------------------------------------------------------
1927 Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
1928 //-----------------------------------------------------------------------------
1929 // This function reads the regularization parameters from the DB
1930 //-----------------------------------------------------------------------------
1931
1932   if(gSystem->AccessPathName(inName,kFileExists)) {
1933     cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
1934     return 0; }
1935
1936   TFile *inFile = TFile::Open(inName);
1937   switch (pdg) {
1938   case 211:
1939     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
1940     new(&fRegParPi) TMatrixD(*fRegPar);
1941     break;
1942   case 321:
1943     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
1944     new(&fRegParKa) TMatrixD(*fRegPar);
1945     break;
1946   case 11:
1947     fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
1948     new(&fRegParEl) TMatrixD(*fRegPar);
1949     break;
1950   }
1951   inFile->Close();
1952
1953   return 1;
1954 }
1955 //-----------------------------------------------------------------------------
1956 void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
1957 //-----------------------------------------------------------------------------
1958 // This function regularizes the elements of the covariance matrix
1959 // that show a momentum depence:
1960 // c00, c11, c22, c33, c44, c20, c24, c40, c31
1961 // The regularization is done separately for pions, kaons, electrons:
1962 // give "Pion","Kaon" or "Electron" as first argument.
1963 //-----------------------------------------------------------------------------
1964
1965   gStyle->SetOptStat(0);
1966   gStyle->SetOptFit(10001);
1967
1968   Int_t thisPdg=211;
1969   Char_t *part="Pions - ";
1970
1971   InitializeKineGrid("DB","points");
1972   SetParticle(pdg);
1973   const Int_t fitbins = fDBgrid->GetBinsPt();
1974   cerr<<" Fit bins:  "<<fitbins<<endl;
1975
1976   switch (pdg) {
1977   case 211: // pions
1978     thisPdg=211;
1979     part="Pions - ";
1980     cerr<<" Processing pions ...\n";
1981     break;
1982   case 321: //kaons
1983     thisPdg=321;
1984     part="Kaons - ";
1985     cerr<<" Processing kaons ...\n";
1986     break;
1987   case 11: // electrons
1988     thisPdg= 11;
1989     part="Electrons - ";
1990     cerr<<" Processing electrons ...\n";
1991     break;
1992   }
1993
1994   // create the chain with all compared tracks
1995   TChain cmptrkchain("cmptrktree");
1996   cmptrkchain.Add("CovMatrix_AllEvts_1.root");
1997   cmptrkchain.Add("CovMatrix_AllEvts_2.root");
1998   cmptrkchain.Add("CovMatrix_AllEvts_3.root");
1999
2000   COMPTRACK cmptrk;
2001   cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
2002   Int_t entries = (Int_t)cmptrkchain.GetEntries();
2003
2004
2005   Int_t pbin;
2006   Int_t * n = new Int_t[fitbins];
2007   Int_t * n00 = new Int_t[fitbins];
2008   Int_t * n11 = new Int_t[fitbins];
2009   Int_t * n20 = new Int_t[fitbins];
2010   Int_t * n22 = new Int_t[fitbins];
2011   Int_t * n31 = new Int_t[fitbins];
2012   Int_t * n33 = new Int_t[fitbins];
2013   Int_t * n40 = new Int_t[fitbins];
2014   Int_t * n42 = new Int_t[fitbins];
2015   Int_t * n44 = new Int_t[fitbins];
2016   Double_t * p = new Double_t[fitbins];
2017   Double_t * ep = new Double_t[fitbins];
2018   Double_t * mean00 = new Double_t[fitbins];
2019   Double_t * mean11 = new Double_t[fitbins];
2020   Double_t * mean20 = new Double_t[fitbins];
2021   Double_t * mean22 = new Double_t[fitbins];
2022   Double_t * mean31 = new Double_t[fitbins];
2023   Double_t * mean33 = new Double_t[fitbins];
2024   Double_t * mean40 = new Double_t[fitbins];
2025   Double_t * mean42 = new Double_t[fitbins];
2026   Double_t * mean44 = new Double_t[fitbins];
2027   Double_t * sigma00 = new Double_t[fitbins];
2028   Double_t * sigma11 = new Double_t[fitbins];
2029   Double_t * sigma20 = new Double_t[fitbins];
2030   Double_t * sigma22 = new Double_t[fitbins];
2031   Double_t * sigma31 = new Double_t[fitbins];
2032   Double_t * sigma33 = new Double_t[fitbins];
2033   Double_t * sigma40 = new Double_t[fitbins];
2034   Double_t * sigma42 = new Double_t[fitbins];
2035   Double_t * sigma44 = new Double_t[fitbins];
2036   Double_t * rmean = new Double_t[fitbins];
2037   Double_t * rsigma = new Double_t[fitbins];
2038   Double_t fitpar[4];
2039
2040   for(Int_t l=0; l<fitbins; l++) {
2041     n[l]=1;
2042     n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
2043     p[l ]=ep[l]=0.;
2044     mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
2045     sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
2046   }
2047
2048   // loop on chain entries for mean
2049   for(Int_t l=0; l<entries; l++) {
2050     cmptrkchain.GetEvent(l);
2051     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2052     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2053     n[pbin]++;
2054     p[pbin]+=cmptrk.p;
2055     mean00[pbin]+=cmptrk.c00;
2056     mean11[pbin]+=cmptrk.c11;
2057     mean20[pbin]+=cmptrk.c20;
2058     mean22[pbin]+=cmptrk.c22;
2059     mean31[pbin]+=cmptrk.c31;
2060     mean33[pbin]+=cmptrk.c33;
2061     mean40[pbin]+=cmptrk.c40;
2062     mean42[pbin]+=cmptrk.c42;
2063     mean44[pbin]+=cmptrk.c44;
2064   } // loop on chain entries
2065
2066   for(Int_t l=0; l<fitbins; l++) {
2067     p[l]/=n[l];
2068     mean00[l]/=n[l];
2069     mean11[l]/=n[l];
2070     mean20[l]/=n[l];
2071     mean22[l]/=n[l];
2072     mean31[l]/=n[l];
2073     mean33[l]/=n[l];
2074     mean40[l]/=n[l];
2075     mean42[l]/=n[l];
2076     mean44[l]/=n[l];
2077   }
2078
2079   // loop on chain entries for sigma
2080   for(Int_t l=0; l<entries; l++) {
2081     cmptrkchain.GetEvent(l);
2082     if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
2083     pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
2084     if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
2085       sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
2086     if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
2087       sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
2088     if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
2089       sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
2090     if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
2091       sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
2092     if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
2093       sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
2094     if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
2095       sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
2096     if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
2097       sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
2098     if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
2099       sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
2100     if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
2101       sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
2102   } // loop on chain entries
2103
2104   for(Int_t l=0; l<fitbins; l++) {
2105     sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
2106     sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
2107     sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
2108     sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
2109     sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
2110     sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
2111     sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
2112     sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
2113     sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
2114   }
2115
2116
2117   // Fit function
2118   TF1 *func = new TF1("FitRegFunc",FitRegFunc,0.23,50.,4);
2119   func->SetParNames("A_meas","A_scatt","B1","B2");
2120
2121   // line to draw on the plots
2122   TLine *lin = new TLine(-1,1,1.69,1);
2123   lin->SetLineStyle(2);
2124   lin->SetLineWidth(2);
2125
2126   // matrix used to store fit results
2127   TMatrixD fitRes(9,4);
2128
2129   //    --- c00 ---
2130
2131   // create the canvas
2132   TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
2133   canv00->Divide(1,2);
2134   // create the graph for cov matrix
2135   TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
2136   TString title00("C(y,y)"); title00.Prepend(part);
2137   TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
2138   frame00->SetXTitle("p [GeV/c]");
2139   canv00->cd(1);  gPad->SetLogx();
2140   frame00->Draw();
2141   gr00->Draw("P");
2142   // Sets initial values for parameters
2143   func->SetParameters(1.6e-3,1.9e-4,1.5,0.);
2144   // Fit points in range defined by function
2145   gr00->Fit("FitRegFunc","R,Q");
2146   func->GetParameters(fitpar);
2147   for(Int_t i=0; i<4; i++) fitRes(0,i)=fitpar[i];
2148   for(Int_t l=0; l<fitbins; l++) {
2149     rmean[l]  = mean00[l]/FitRegFunc(&p[l],fitpar);
2150     rsigma[l] = sigma00[l]/FitRegFunc(&p[l],fitpar);
2151   }
2152   // create the graph the regularized cov. matrix
2153   TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2154   TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2155   regtitle00.Prepend(part);
2156   TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
2157   frame00reg->SetXTitle("p [GeV/c]");
2158   canv00->cd(2); gPad->SetLogx();
2159   frame00reg->Draw();
2160   gr00reg->Draw("P");
2161   lin->Draw("same");
2162
2163
2164   //    --- c11 ---
2165
2166   // create the canvas
2167   TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
2168   canv11->Divide(1,2);
2169   // create the graph for cov matrix
2170   TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
2171   TString title11("C(z,z)"); title11.Prepend(part);
2172   TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
2173   frame11->SetXTitle("p [GeV/c]");
2174   canv11->cd(1);  gPad->SetLogx();
2175   frame11->Draw();
2176   gr11->Draw("P");
2177   // Sets initial values for parameters
2178   func->SetParameters(1.2e-3,8.1e-4,1.,0.);
2179   // Fit points in range defined by function
2180   gr11->Fit("FitRegFunc","R,Q");
2181   func->GetParameters(fitpar);
2182   for(Int_t i=0; i<4; i++) fitRes(1,i)=fitpar[i];
2183   for(Int_t l=0; l<fitbins; l++) {
2184     rmean[l]  = mean11[l]/FitRegFunc(&p[l],fitpar);
2185     rsigma[l] = sigma11[l]/FitRegFunc(&p[l],fitpar);
2186   }
2187   // create the graph the regularized cov. matrix
2188   TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2189   TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2190   regtitle11.Prepend(part);
2191   TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
2192   frame11reg->SetXTitle("p [GeV/c]");
2193   canv11->cd(2); gPad->SetLogx();
2194   frame11reg->Draw();
2195   gr11reg->Draw("P");
2196   lin->Draw("same");
2197
2198
2199   //    --- c20 ---
2200
2201   // create the canvas
2202   TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
2203   canv20->Divide(1,2);
2204   // create the graph for cov matrix
2205   TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
2206   TString title20("C(#eta, y)"); title20.Prepend(part);
2207   TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
2208   frame20->SetXTitle("p [GeV/c]");
2209   canv20->cd(1);  gPad->SetLogx();
2210   frame20->Draw();
2211   gr20->Draw("P");
2212   // Sets initial values for parameters
2213   func->SetParameters(7.3e-5,1.2e-5,1.5,0.);
2214   // Fit points in range defined by function
2215   gr20->Fit("FitRegFunc","R,Q");
2216   func->GetParameters(fitpar);
2217   for(Int_t i=0; i<4; i++) fitRes(2,i)=fitpar[i];
2218   for(Int_t l=0; l<fitbins; l++) {
2219     rmean[l]  = mean20[l]/FitRegFunc(&p[l],fitpar);
2220     rsigma[l] = sigma20[l]/FitRegFunc(&p[l],fitpar);
2221   }
2222   // create the graph the regularized cov. matrix
2223   TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2224   TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2225   regtitle20.Prepend(part);
2226   TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
2227   frame20reg->SetXTitle("p [GeV/c]");
2228   canv20->cd(2); gPad->SetLogx();
2229   frame20reg->Draw();
2230   gr20reg->Draw("P");
2231   lin->Draw("same");
2232
2233
2234   //    --- c22 ---
2235
2236   // create the canvas
2237   TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
2238   canv22->Divide(1,2);
2239   // create the graph for cov matrix
2240   TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
2241   TString title22("C(#eta, #eta)"); title22.Prepend(part);
2242   TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
2243   frame22->SetXTitle("p [GeV/c]");
2244   canv22->cd(1);  gPad->SetLogx();
2245   frame22->Draw();
2246   gr22->Draw("P");
2247   // Sets initial values for parameters
2248   func->SetParameters(5.2e-6,1.1e-6,2.,1.);
2249   // Fit points in range defined by function
2250   gr22->Fit("FitRegFunc","R,Q");
2251   func->GetParameters(fitpar);
2252   for(Int_t i=0; i<4; i++) fitRes(3,i)=fitpar[i];
2253   for(Int_t l=0; l<fitbins; l++) {
2254     rmean[l]  = mean22[l]/FitRegFunc(&p[l],fitpar);
2255     rsigma[l] = sigma22[l]/FitRegFunc(&p[l],fitpar);
2256   }
2257   // create the graph the regularized cov. matrix
2258   TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2259   TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2260   regtitle22.Prepend(part);
2261   TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
2262   frame22reg->SetXTitle("p [GeV/c]");
2263   canv22->cd(2); gPad->SetLogx();
2264   frame22reg->Draw();
2265   gr22reg->Draw("P");
2266   lin->Draw("same");
2267
2268
2269   //    --- c31 ---
2270
2271   // create the canvas
2272   TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
2273   canv31->Divide(1,2);
2274   // create the graph for cov matrix
2275   TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
2276   TString title31("C(tg #lambda,z)"); title31.Prepend(part);
2277   TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
2278   frame31->SetXTitle("p [GeV/c]");
2279   canv31->cd(1);  gPad->SetLogx();
2280   frame31->Draw();
2281   gr31->Draw("P");
2282   // Sets initial values for parameters
2283   func->SetParameters(-1.2e-5,-1.2e-5,1.5,3.);
2284   // Fit points in range defined by function
2285   gr31->Fit("FitRegFunc","R,Q");
2286   func->GetParameters(fitpar);
2287   for(Int_t i=0; i<4; i++) fitRes(4,i)=fitpar[i];
2288   for(Int_t l=0; l<fitbins; l++) {
2289     rmean[l]  = mean31[l]/FitRegFunc(&p[l],fitpar);
2290     rsigma[l] = -sigma31[l]/FitRegFunc(&p[l],fitpar);
2291   }
2292   // create the graph the regularized cov. matrix
2293   TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2294   TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2295   regtitle31.Prepend(part);
2296   TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
2297   frame31reg->SetXTitle("p [GeV/c]");
2298   canv31->cd(2); gPad->SetLogx();
2299   frame31reg->Draw();
2300   gr31reg->Draw("P");
2301   lin->Draw("same");
2302
2303
2304   //    --- c33 ---
2305
2306   // create the canvas
2307   TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
2308   canv33->Divide(1,2);
2309   // create the graph for cov matrix
2310   TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
2311   TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
2312   TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
2313   frame33->SetXTitle("p [GeV/c]");
2314   canv33->cd(1);  gPad->SetLogx();
2315   frame33->Draw();
2316   gr33->Draw("P");
2317   // Sets initial values for parameters
2318   func->SetParameters(1.3e-7,4.6e-7,1.7,4.);
2319   // Fit points in range defined by function
2320   gr33->Fit("FitRegFunc","R,Q");
2321   func->GetParameters(fitpar);
2322   for(Int_t i=0; i<4; i++) fitRes(5,i)=fitpar[i];
2323   for(Int_t l=0; l<fitbins; l++) {
2324     rmean[l]  = mean33[l]/FitRegFunc(&p[l],fitpar);
2325     rsigma[l] = sigma33[l]/FitRegFunc(&p[l],fitpar);
2326   }
2327   // create the graph the regularized cov. matrix
2328   TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2329   TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2330   regtitle33.Prepend(part);
2331   TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
2332   frame33reg->SetXTitle("p [GeV/c]");
2333   canv33->cd(2); gPad->SetLogx();
2334   frame33reg->Draw();
2335   gr33reg->Draw("P");
2336   lin->Draw("same");
2337
2338
2339   //    --- c40 ---
2340
2341   // create the canvas
2342   TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
2343   canv40->Divide(1,2);
2344   // create the graph for cov matrix
2345   TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
2346   TString title40("C(C,y)"); title40.Prepend(part);
2347   TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
2348   frame40->SetXTitle("p [GeV/c]");
2349   canv40->cd(1);  gPad->SetLogx();
2350   frame40->Draw();
2351   gr40->Draw("P");
2352   // Sets initial values for parameters
2353   func->SetParameters(4.e-7,4.4e-8,1.5,0.);
2354   // Fit points in range defined by function
2355   gr40->Fit("FitRegFunc","R,Q");
2356   func->GetParameters(fitpar);
2357   for(Int_t i=0; i<4; i++) fitRes(6,i)=fitpar[i];
2358   for(Int_t l=0; l<fitbins; l++) {
2359     rmean[l]  = mean40[l]/FitRegFunc(&p[l],fitpar);
2360     rsigma[l] = sigma40[l]/FitRegFunc(&p[l],fitpar);
2361   }
2362   // create the graph the regularized cov. matrix
2363   TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2364   TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2365   regtitle40.Prepend(part);
2366   TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
2367   frame40reg->SetXTitle("p [GeV/c]");
2368   canv40->cd(2); gPad->SetLogx();
2369   frame40reg->Draw();
2370   gr40reg->Draw("P");
2371   lin->Draw("same");
2372
2373
2374   //    --- c42 ---
2375
2376   // create the canvas
2377   TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
2378   canv42->Divide(1,2);
2379   // create the graph for cov matrix
2380   TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
2381   TString title42("C(C, #eta)"); title42.Prepend(part);
2382   TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
2383   frame42->SetXTitle("p [GeV/c]");
2384   canv42->cd(1);  gPad->SetLogx();
2385   frame42->Draw();
2386   gr42->Draw("P");
2387   // Sets initial values for parameters
2388   func->SetParameters(3.e-8,8.2e-9,2.,1.);
2389   // Fit points in range defined by function
2390   gr42->Fit("FitRegFunc","R,Q");
2391   func->GetParameters(fitpar);
2392   for(Int_t i=0; i<4; i++) fitRes(7,i)=fitpar[i];
2393   for(Int_t l=0; l<fitbins; l++) {
2394     rmean[l]  = mean42[l]/FitRegFunc(&p[l],fitpar);
2395     rsigma[l] = sigma42[l]/FitRegFunc(&p[l],fitpar);
2396   }
2397   // create the graph the regularized cov. matrix
2398   TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2399   TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2400   regtitle42.Prepend(part);
2401   TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
2402   frame42reg->SetXTitle("p [GeV/c]");
2403   canv42->cd(2); gPad->SetLogx();
2404   frame42reg->Draw();
2405   gr42reg->Draw("P");
2406   lin->Draw("same");
2407
2408
2409   //    --- c44 ---
2410
2411   // create the canvas
2412   TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
2413   canv44->Divide(1,2);
2414   // create the graph for cov matrix
2415   TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
2416   TString title44("C(C,C)"); title44.Prepend(part);
2417   TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
2418   frame44->SetXTitle("p [GeV/c]");
2419   canv44->cd(1);  gPad->SetLogx();
2420   frame44->Draw();
2421   gr44->Draw("P");
2422   // Sets initial values for parameters
2423   func->SetParameters(1.8e-10,5.8e-11,2.,3.);
2424   // Fit points in range defined by function
2425   gr44->Fit("FitRegFunc","R,Q");
2426   func->GetParameters(fitpar);
2427   for(Int_t i=0; i<4; i++) fitRes(8,i)=fitpar[i];
2428   for(Int_t l=0; l<fitbins; l++) {
2429     rmean[l]  = mean44[l]/FitRegFunc(&p[l],fitpar);
2430     rsigma[l] = sigma44[l]/FitRegFunc(&p[l],fitpar);
2431   }
2432   // create the graph the regularized cov. matrix
2433   TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
2434   TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
2435   regtitle44.Prepend(part);
2436   TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
2437   frame44reg->SetXTitle("p [GeV/c]");
2438   canv44->cd(2); gPad->SetLogx();
2439   frame44reg->Draw();
2440   gr44reg->Draw("P");
2441   lin->Draw("same");
2442
2443
2444
2445
2446   switch (pdg) {
2447   case 211:
2448     new(&fRegParPi) TMatrixD(fitRes);
2449     break;
2450   case 321:
2451     new(&fRegParKa) TMatrixD(fitRes);
2452     break;
2453   case 11:
2454     new(&fRegParEl) TMatrixD(fitRes);
2455     break;
2456   }
2457
2458   // write fit parameters to file
2459   WriteRegParams(outName,pdg);
2460
2461   delete [] n;
2462   delete [] n00;
2463   delete [] n11;
2464   delete [] n20;
2465   delete [] n22;
2466   delete [] n31;
2467   delete [] n33;
2468   delete [] n40;
2469   delete [] n42;
2470   delete [] n44;
2471   delete [] p;
2472   delete [] ep;
2473   delete [] mean00;
2474   delete [] mean11;
2475   delete [] mean20;
2476   delete [] mean22;
2477   delete [] mean31;
2478   delete [] mean33;
2479   delete [] mean40;
2480   delete [] mean42;
2481   delete [] mean44;
2482   delete [] sigma00;
2483   delete [] sigma11;
2484   delete [] sigma20;
2485   delete [] sigma22;
2486   delete [] sigma31;
2487   delete [] sigma33;
2488   delete [] sigma40;
2489   delete [] sigma42;
2490   delete [] sigma44;
2491   delete [] rmean;
2492   delete [] rsigma;
2493
2494   return;
2495 }
2496 //-----------------------------------------------------------------------------
2497 Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
2498 //-----------------------------------------------------------------------------
2499 // This function makes a selection according to TPC tracking efficiency
2500 //-----------------------------------------------------------------------------
2501
2502   Double_t eff=0.;
2503
2504   eff = fEff->GetValueAt(pt,eta);
2505
2506   if(gRandom->Rndm() < eff) return kTRUE;
2507
2508   return kFALSE;
2509 }
2510 //-----------------------------------------------------------------------------
2511 void AliTPCtrackerParam::SetParticle(Int_t pdg) {
2512 //-----------------------------------------------------------------------------
2513 // This function set efficiencies, pulls, etc... for the particle
2514 // specie of the current particle
2515 //-----------------------------------------------------------------------------
2516
2517   switch (pdg) {
2518   case 211:
2519     fDBgrid = &fDBgridPi;
2520     fEff    = &fEffPi;
2521     fPulls  =  fPullsPi;
2522     fRegPar = &fRegParPi;
2523     fdEdxMean = &fdEdxMeanPi;
2524     fdEdxRMS  = &fdEdxRMSPi;
2525     break;
2526   case 321:
2527     fDBgrid = &fDBgridKa;
2528     fEff    = &fEffKa;
2529     fPulls  =  fPullsKa;
2530     fRegPar = &fRegParKa;
2531     fdEdxMean = &fdEdxMeanKa;
2532     fdEdxRMS  = &fdEdxRMSKa;
2533     break;
2534   case 2212:
2535     fDBgrid = &fDBgridPi;
2536     fEff    = &fEffPr;
2537     fPulls  =  fPullsPi;
2538     fRegPar = &fRegParPi;
2539     fdEdxMean = &fdEdxMeanPr;
2540     fdEdxRMS  = &fdEdxRMSPr;
2541     break;
2542   case 11:
2543     fDBgrid = &fDBgridEl;
2544     fEff    = &fEffEl;
2545     fPulls  =  fPullsEl;
2546     fRegPar = &fRegParEl;
2547     fdEdxMean = &fdEdxMeanEl;
2548     fdEdxRMS  = &fdEdxRMSEl;
2549     break;
2550   case 13:
2551     fDBgrid = &fDBgridPi;
2552     fEff    = &fEffMu;
2553     fPulls  =  fPullsPi;
2554     fRegPar = &fRegParPi;
2555     fdEdxMean = &fdEdxMeanPi;
2556     fdEdxRMS  = &fdEdxRMSPi;
2557     break;
2558   }
2559
2560   return;
2561 }
2562 //-----------------------------------------------------------------------------
2563 void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
2564   const {
2565 //-----------------------------------------------------------------------------
2566 // This function smears track parameters according to streched cov. matrix
2567 //-----------------------------------------------------------------------------
2568   AliGausCorr *corgen = new AliGausCorr(cov,5);
2569   TArrayD corr(5);
2570   corgen->GetGaussN(corr);
2571   delete corgen;
2572   corgen = 0;
2573
2574   for(Int_t l=0;l<5;l++) {
2575     xxsm[l] = xx[l]+corr[l];
2576   }
2577
2578   return;
2579 }
2580 //-----------------------------------------------------------------------------
2581 Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) {
2582 //-----------------------------------------------------------------------------
2583 // This function writes the dEdx parameters to the DB
2584 //-----------------------------------------------------------------------------
2585
2586   Option_t *opt;
2587   Char_t *dirName="Pions";
2588   Char_t *meanName="dEdxMeanPi";
2589   Char_t *RMSName="dEdxRMSPi";
2590
2591   SetParticle(pdg);
2592
2593   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
2594   } else { opt="update"; }
2595
2596   switch (pdg) {
2597   case 211:
2598     dirName="Pions";
2599     meanName="dEdxMeanPi";
2600     RMSName="dEdxRMSPi";
2601     break;
2602   case 321:
2603     dirName="Kaons";
2604     meanName="dEdxMeanKa";
2605     RMSName="dEdxRMSKa";
2606     break;
2607   case 2212:
2608     dirName="Protons";
2609     meanName="dEdxMeanPr";
2610     RMSName="dEdxRMSPr";
2611     break;
2612   case 11:
2613     dirName="Electrons";
2614     meanName="dEdxMeanEl";
2615     RMSName="dEdxRMSEl";
2616     break;
2617   }
2618
2619   TFile *outFile = new TFile(outName,opt);
2620   if(!gDirectory->cd("/dEdx")) {
2621     outFile->mkdir("dEdx");
2622     gDirectory->cd("/dEdx");
2623   }
2624   TDirectory *dir2 = gDirectory->mkdir(dirName);
2625   dir2->cd();
2626   fdEdxMean->SetName(meanName); fdEdxMean->Write();
2627   fdEdxRMS->SetName(RMSName);  fdEdxRMS->Write();
2628
2629   outFile->Close();
2630   delete outFile;
2631
2632
2633   return 1;
2634 }
2635 //-----------------------------------------------------------------------------
2636 Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) {
2637 //-----------------------------------------------------------------------------
2638 // This function writes the TPC efficiencies to the DB
2639 //-----------------------------------------------------------------------------
2640
2641
2642   Option_t *opt;
2643
2644   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
2645   } else { opt="update"; }
2646
2647   TFile *outFile = new TFile(outName,opt);
2648
2649   outFile->mkdir("Efficiencies");
2650   gDirectory->cd("/Efficiencies");
2651   gDirectory->mkdir("Pions");
2652   gDirectory->mkdir("Kaons");
2653   gDirectory->mkdir("Protons");
2654   gDirectory->mkdir("Electrons");
2655   gDirectory->mkdir("Muons");
2656
2657   gDirectory->cd("/Efficiencies/Pions");
2658   fEffPi.SetName("EffPi");
2659   fEffPi.Write();
2660   gDirectory->cd("/Efficiencies/Kaons");
2661   fEffKa.SetName("EffKa");
2662   fEffKa.Write();
2663   gDirectory->cd("/Efficiencies/Protons");
2664   fEffPr.SetName("EffPr");
2665   fEffPr.Write();
2666   gDirectory->cd("/Efficiencies/Electrons");
2667   fEffEl.SetName("EffEl");
2668   fEffEl.Write();
2669   gDirectory->cd("/Efficiencies/Muons");
2670   fEffMu.SetName("EffMu");
2671   fEffMu.Write();
2672
2673   outFile->Close();
2674
2675   delete outFile;
2676
2677   return 1;
2678 }
2679 //-----------------------------------------------------------------------------
2680 Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) {
2681 //-----------------------------------------------------------------------------
2682 // This function writes the pulls to the DB
2683 //-----------------------------------------------------------------------------
2684
2685   Option_t *opt;
2686
2687   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
2688   } else { opt="update"; }
2689
2690   TFile *outFile = new TFile(outName,opt);
2691
2692   outFile->mkdir("Pulls");
2693   gDirectory->cd("/Pulls");
2694   gDirectory->mkdir("Pions");
2695   gDirectory->mkdir("Kaons");
2696   gDirectory->mkdir("Electrons");
2697
2698   for(Int_t i=0;i<5;i++) {
2699     TString pi("PullsPi_"); pi+=i;
2700     TString ka("PullsKa_"); ka+=i;
2701     TString el("PullsEl_"); el+=i;
2702     fPullsPi[i].SetName(pi.Data());
2703     fPullsKa[i].SetName(ka.Data());
2704     fPullsEl[i].SetName(el.Data());
2705     gDirectory->cd("/Pulls/Pions");
2706     fPullsPi[i].Write();
2707     gDirectory->cd("/Pulls/Kaons");
2708     fPullsKa[i].Write();
2709     gDirectory->cd("/Pulls/Electrons");
2710     fPullsEl[i].Write();
2711   }
2712   outFile->Close();
2713   delete outFile;
2714
2715   return 1;
2716 }
2717 //-----------------------------------------------------------------------------
2718 Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) {
2719 //-----------------------------------------------------------------------------
2720 // This function writes the regularization parameters to the DB
2721 //-----------------------------------------------------------------------------
2722
2723   Option_t *opt;
2724   Char_t *dirName="Pions";
2725   Char_t *keyName="RegPions";
2726
2727   SetParticle(pdg);
2728
2729   if(gSystem->AccessPathName(outName,kFileExists))  { opt="recreate";
2730   } else { opt="update"; }
2731
2732   switch (pdg) {
2733   case 211:
2734     dirName="Pions";
2735     keyName="RegPions";
2736     break;
2737   case 321:
2738     dirName="Kaons";
2739     keyName="RegKaons";
2740     break;
2741   case 11:
2742     dirName="Electrons";
2743     keyName="RegElectrons";
2744     break;
2745   }
2746
2747   TFile *outFile = new TFile(outName,opt);
2748   if(!gDirectory->cd("/RegParams")) {
2749     outFile->mkdir("RegParams");
2750     gDirectory->cd("/RegParams");
2751   }
2752   TDirectory *dir2 = gDirectory->mkdir(dirName);
2753   dir2->cd();
2754   fRegPar->Write(keyName);
2755
2756   outFile->Close();
2757   delete outFile;
2758
2759
2760   return 1;
2761 }
2762
2763
2764
2765
2766
2767
2768
2769
2770