PWG1/AliComparisonEff.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliComparisonEff class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms) are stored in the output picture_eff.root file.
   8 //
   9 // Author: J.Otwinowski 04/02/2008
  10 //------------------------------------------------------------------------------
  11
  12 /*
  13   // after running analysis, read the file, and get component
  14   gSystem->Load("libPWG1.so");
  15   TFile f("Output.root");
  16   AliComparisonEff * comp = (AliComparisonEff*)f.Get("AliComparisonEff");
  17
  18
  19   // analyse comparison data (output stored in pictures_eff.root)
  20   comp->Analyse();
  21
  22   // paramtetrisation of the TPC track length (for information only)
  23   TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2);  // length function
  24   TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
  25   fl2.SetParameter(1,1);
  26   fl2.SetParameter(0,1);
  27 */
  28
  29
  30 #include <iostream>
  31
  32 #include "TFile.h"
  33 #include "TCint.h"
  34 #include "TH3F.h"
  35 #include "TH2F.h"
  36 #include "TF1.h"
  37 #include "TProfile.h"
  38 #include "TProfile2D.h"
  39 #include "TGraph2D.h"
  40 #include "TCanvas.h"
  41 #include "TGraph.h"
  42 //
  43 #include "AliESDEvent.h"
  44 #include "AliESD.h"
  45 #include "AliESDfriend.h"
  46 #include "AliESDfriendTrack.h"
  47 #include "AliRecInfoCuts.h"
  48 #include "AliMCInfoCuts.h"
  49 #include "AliLog.h"
  50 //
  51 #include "AliMathBase.h"
  52 #include "AliTreeDraw.h"
  53 #include "AliMagFMaps.h"
  54 #include "AliESDVertex.h"
  55 #include "AliExternalTrackParam.h"
  56 #include "AliTracker.h"
  57
  58 #include "AliMCInfo.h"
  59 #include "AliESDRecInfo.h"
  60 #include "AliComparisonEff.h"
  61
  62 using namespace std;
  63
  64
  65 ClassImp(AliComparisonEff)
  66
  67 //_____________________________________________________________________________
  68 AliComparisonEff::AliComparisonEff():
  69   TNamed("AliComparisonEff","AliComparisonEff"),
  70
  71   // histograms
  72
  73   fMCPt(0),
  74   fMCRecPt(0),
  75   fMCRecPrimPt(0),
  76   fMCRecSecPt(0),
  77
  78   fEffTPCPt(0),
  79   fEffTPCPtMC(0),
  80   fEffTPCPtF(0),
  81   //
  82   fEffTPCPt_P(0),
  83   fEffTPCPtMC_P(0),
  84   fEffTPCPtF_P(0),
  85   //
  86   fEffTPCPt_Pi(0),
  87   fEffTPCPtMC_Pi(0),
  88   fEffTPCPtF_Pi(0),
  89   //
  90   fEffTPCPt_K(0),
  91   fEffTPCPtMC_K(0),
  92   fEffTPCPtF_K(0),
  93
  94   fEffTPCTan(0),
  95   fEffTPCTanMC(0),
  96   fEffTPCTanF(0),
  97   //
  98   fEffTPCPtTan(0),
  99   fEffTPCPtTanMC(0),
 100   fEffTPCPtTanF(0),
 101
 102   // Cuts
 103   fCutsRC(0),
 104   fCutsMC(0),
 105
 106   fVertex(0)
 107 {
 108   // init vertex
 109   fVertex = new AliESDVertex();
 110   fVertex->SetXv(0.0); fVertex->SetYv(0.0); fVertex->SetZv(0.0);
 111
 112   for(Int_t i=0; i<4; ++i)
 113   {
 114     fTPCPtDCASigmaIdeal[i]=0;
 115     fTPCPtDCASigmaFull[i]=0;
 116     fTPCPtDCASigmaDay0[i]=0;
 117
 118     fTPCPtDCAXY[i]=0;
 119     fTPCPtDCAZ[i]=0;
 120
 121         fTPCPtDCASigmaIdealPid[i]=0;
 122         fTPCPtDCASigmaFullPid[i]=0;
 123         fTPCPtDCASigmaDay0Pid[i]=0;
 124
 125         fTPCPtDCAXYPid[i]=0;
 126         fTPCPtDCAZPid[i]=0;
 127   }
 128   InitHisto();
 129   InitCuts();
 130 }
 131
 132 //_____________________________________________________________________________
 133 AliComparisonEff::~AliComparisonEff(){
 134
 135   //
 136   if(fMCPt)  delete  fEffTPCPt; fEffTPCPt=0;
 137   if(fMCRecPt)  delete  fMCRecPt; fMCRecPt=0;
 138   if(fMCRecPrimPt)  delete  fMCRecPrimPt; fMCRecPrimPt=0;
 139   if(fMCRecSecPt)  delete  fMCRecSecPt; fMCRecSecPt=0;
 140
 141   //
 142   if(fEffTPCPt)   delete  fEffTPCPt;   fEffTPCPt=0;
 143   if(fEffTPCPtMC) delete  fEffTPCPtMC; fEffTPCPtMC=0;
 144   if(fEffTPCPtF)  delete  fEffTPCPtF;  fEffTPCPtF=0;
 145
 146   //
 147   if(fEffTPCPt_P)   delete  fEffTPCPt_P;   fEffTPCPt_P=0;
 148   if(fEffTPCPtMC_P) delete  fEffTPCPtMC_P; fEffTPCPtMC_P=0;
 149   if(fEffTPCPtF_P)  delete  fEffTPCPtF_P;  fEffTPCPtF_P=0;
 150
 151   //
 152   if(fEffTPCPt_Pi)   delete  fEffTPCPt_Pi;   fEffTPCPt_Pi=0;
 153   if(fEffTPCPtMC_Pi) delete  fEffTPCPtMC_Pi; fEffTPCPtMC_Pi=0;
 154   if(fEffTPCPtF_Pi)  delete  fEffTPCPtF_Pi;  fEffTPCPtF_Pi=0;
 155
 156   //
 157   if(fEffTPCPt_K)   delete  fEffTPCPt_K;   fEffTPCPt_K=0;
 158   if(fEffTPCPtMC_K) delete  fEffTPCPtMC_K; fEffTPCPtMC_K=0;
 159   if(fEffTPCPtF_K)  delete  fEffTPCPtF_K;  fEffTPCPtF_K=0;
 160
 161   //
 162   if(fEffTPCTan)   delete  fEffTPCTan;   fEffTPCTan=0;
 163   if(fEffTPCTanMC) delete  fEffTPCTanMC; fEffTPCTanMC=0;
 164   if(fEffTPCTanF)  delete  fEffTPCTanF;  fEffTPCTanF=0;
 165
 166   //
 167   if(fEffTPCPtTan)   delete  fEffTPCPtTan;   fEffTPCPtTan=0;
 168   if(fEffTPCPtTanMC) delete  fEffTPCPtTanMC; fEffTPCPtTanMC=0;
 169   if(fEffTPCPtTanF)  delete  fEffTPCPtTanF;  fEffTPCPtTanF=0;
 170
 171   for(Int_t i=0; i<4; ++i)
 172   {
 173     if(fTPCPtDCASigmaIdeal[i]) delete  fTPCPtDCASigmaIdeal[i];  fTPCPtDCASigmaIdeal[i]=0;
 174     if(fTPCPtDCASigmaFull[i]) delete  fTPCPtDCASigmaFull[i];  fTPCPtDCASigmaFull[i]=0;
 175     if(fTPCPtDCASigmaDay0[i]) delete  fTPCPtDCASigmaDay0[i];  fTPCPtDCASigmaDay0[i]=0;
 176
 177     if(fTPCPtDCAXY[i]) delete  fTPCPtDCAXY[i];  fTPCPtDCAXY[i]=0;
 178     if(fTPCPtDCAZ[i]) delete  fTPCPtDCAZ[i];  fTPCPtDCAZ[i]=0;
 179
 180         if(fTPCPtDCASigmaIdealPid[i]) delete  fTPCPtDCASigmaIdealPid[i];  fTPCPtDCASigmaIdealPid[i]=0;
 181         if(fTPCPtDCASigmaFullPid[i]) delete  fTPCPtDCASigmaFullPid[i];  fTPCPtDCASigmaFullPid[i]=0;
 182         if(fTPCPtDCASigmaDay0Pid[i]) delete  fTPCPtDCASigmaDay0Pid[i];  fTPCPtDCASigmaDay0Pid[i]=0;
 183
 184         if(fTPCPtDCAXYPid[i]) delete  fTPCPtDCAXYPid[i]; fTPCPtDCAXYPid[i]=0;
 185         if(fTPCPtDCAZPid[i]) delete   fTPCPtDCAZPid[i];  fTPCPtDCAZPid[i]=0;
 186   }
 187 }
 188
 189 //_____________________________________________________________________________
 190 void AliComparisonEff::InitHisto(){
 191
 192   // Init histograms
 193   //
 194   fMCPt = new TH1F("fMCPt","fMCPt",50,0.1,3);
 195   fMCPt->SetXTitle("p_{t}");
 196   fMCPt->SetYTitle("yield");
 197
 198   fMCRecPt = new TH1F("fMCRecPt","fMCRecPt",50,0.1,3);
 199   fMCRecPt->SetXTitle("p_{t}");
 200   fMCRecPt->SetYTitle("yield");
 201
 202   fMCRecPrimPt = new TH1F("fMCRecPrimPt","fMCRecPrimPt",50,0.1,3);
 203   fMCRecPrimPt->SetXTitle("p_{t}");
 204   fMCRecPrimPt->SetYTitle("yield");
 205
 206   fMCRecSecPt = new TH1F("fMCRecSecPt","fMCRecSecPt",50,0.1,3);
 207   fMCRecSecPt->SetXTitle("p_{t}");
 208   fMCRecSecPt->SetYTitle("yield");
 209
 210   // Efficiency as function of pt
 211   fEffTPCPt = new TProfile("Eff_pt","Eff_Pt",50,0.1,3);
 212   fEffTPCPt->SetXTitle("p_{t}");
 213   fEffTPCPt->SetYTitle("TPC Efficiency");
 214
 215   fEffTPCPtMC = new TProfile("MC_Eff_pt","MC_Eff_Pt",50,0.1,3);
 216   fEffTPCPtMC->SetXTitle("p_{t}");
 217   fEffTPCPtMC->SetYTitle("MC TPC Efficiency");
 218
 219   fEffTPCPtF = new TProfile("F_Eff_pt","F_Eff_Pt",50,0.1,3);
 220   fEffTPCPtF->SetXTitle("p_{t}");
 221   fEffTPCPtF->SetYTitle("TPC Findable Efficiency");
 222
 223   // Efficiency as function of pt protons
 224   fEffTPCPt_P = new TProfile("Eff_pt_P","Eff_Pt_P",50,0.1,3);
 225   fEffTPCPt_P->SetXTitle("p_{t}");
 226   fEffTPCPt_P->SetYTitle("TPC Efficiency");
 227
 228   fEffTPCPtMC_P = new TProfile("MC_Eff_pt_P","MC_Eff_Pt_P",50,0.1,3);
 229   fEffTPCPtMC_P->SetXTitle("p_{t}");
 230   fEffTPCPtMC_P->SetYTitle("MC TPC Efficiency");
 231
 232   fEffTPCPtF_P = new TProfile("F_Eff_pt_P","F_Eff_Pt_P",50,0.1,3);
 233   fEffTPCPtF_P->SetXTitle("p_{t}");
 234   fEffTPCPtF_P->SetYTitle("TPC Findable Efficiency");
 235
 236   // Efficiency as function of pt pions
 237   fEffTPCPt_Pi = new TProfile("Eff_pt_Pi","Eff_Pit_Pi",50,0.1,3);
 238   fEffTPCPt_Pi->SetXTitle("p_{t}");
 239   fEffTPCPt_Pi->SetYTitle("TPC Efficiency");
 240
 241   fEffTPCPtMC_Pi = new TProfile("MC_Eff_pt_Pi","MC_Eff_Pit_Pi",50,0.1,3);
 242   fEffTPCPtMC_Pi->SetXTitle("p_{t}");
 243   fEffTPCPtMC_Pi->SetYTitle("MC TPC Efficiency");
 244
 245   fEffTPCPtF_Pi = new TProfile("F_Eff_pt_Pi","F_Eff_Pit_Pi",50,0.1,3);
 246   fEffTPCPtF_Pi->SetXTitle("p_{t}");
 247   fEffTPCPtF_Pi->SetYTitle("TPC Findable Efficiency");
 248
 249   // Efficiency as function of pt kaons
 250   fEffTPCPt_K = new TProfile("Eff_pt_K","Eff_Kt_K",50,0.1,3);
 251   fEffTPCPt_K->SetXTitle("p_{t}");
 252   fEffTPCPt_K->SetYTitle("TPC Efficiency");
 253
 254   fEffTPCPtMC_K = new TProfile("MC_Eff_pt_K","MC_Eff_Kt_K",50,0.1,3);
 255   fEffTPCPtMC_K->SetXTitle("p_{t}");
 256   fEffTPCPtMC_K->SetYTitle("MC TPC Efficiency");
 257
 258   fEffTPCPtF_K = new TProfile("F_Eff_pt_K","F_Eff_Kt_K",50,0.1,3);
 259   fEffTPCPtF_K->SetXTitle("p_{t}");
 260   fEffTPCPtF_K->SetYTitle("TPC Findable Efficiency");
 261
 262   // Efficiency as function of tan(theta)
 263   fEffTPCTan = new TProfile("Eff_tan","Eff_tan",50,-2.5,2.5);
 264   fEffTPCTan->SetXTitle("tan(#theta)");
 265   fEffTPCTan->SetYTitle("TPC Efficiency");
 266
 267   fEffTPCTanMC = new TProfile("MC_Eff_tan","MC_Eff_tan",50,-2.5,2.5);
 268   fEffTPCTanMC->SetXTitle("tan(#theta)");
 269   fEffTPCTanMC->SetYTitle("MC TPC Efficiency");
 270
 271   fEffTPCTanF = new TProfile("F_Eff_tan","F_Eff_tan",50,-2.5,2.5);
 272   fEffTPCPtF->SetXTitle("tan(#theta)");
 273   fEffTPCPtF->SetYTitle("TPC Findable Efficiency");
 274
 275   // Efficiency as function of pt and tan(theta)
 276   fEffTPCPtTan = new TProfile2D("Eff_pt_tan","Eff_Pt_tan",10,0.1,3,20,-2.,2.);
 277   fEffTPCPtTan->SetXTitle("tan(#theta)");
 278   fEffTPCPtTan->SetYTitle("p_{t}");
 279
 280   fEffTPCPtTanMC = new TProfile2D("MC_Eff_pt_tan_MC","MC Eff Pt",10,0.1,3,20,-2.,2.);
 281   fEffTPCPtTanMC->SetXTitle("tan(#theta)");
 282   fEffTPCPtTanMC->SetYTitle("p_{t}");
 283
 284   fEffTPCPtTanF = new TProfile2D("MC_Eff_pt_tan_F","MC Eff Pt",10,0.1,3,20,-2.,2.);
 285   fEffTPCPtTanF->SetXTitle("tan(#theta)");
 286   fEffTPCPtTanF->SetYTitle("p_{t}");
 287
 288   char name[256];
 289   for(Int_t i=0; i<4; ++i)
 290   {
 291     sprintf(name, "fTPCPtDCASigmaIdeal_%d",i);
 292     fTPCPtDCASigmaIdeal[i] = new TH2F(name,name,50,0.1,3,100,0,100);
 293
 294     sprintf(name, "fTPCPtDCASigmaFull_%d",i);
 295     fTPCPtDCASigmaFull[i] = new TH2F(name,name,50,0.1,3,100,0,100);
 296
 297     sprintf(name, "fTPCPtDCASigmaDay0_%d",i);
 298     fTPCPtDCASigmaDay0[i] = new TH2F(name,name,50,0.1,3,100,0,100);
 299
 300     sprintf(name, "fTPCPtDCAXY_%d",i);
 301     fTPCPtDCAXY[i]= new TH2F(name,name,50,0.1,3,100,0,100);
 302     sprintf(name, "fTPCPtDCAZ_%d",i);
 303     fTPCPtDCAZ[i]= new TH2F(name,name,50,0.1,3,100,0,100);
 304
 305     sprintf(name, "fTPCPtDCASigmaIdealPid_%d",i);
 306         fTPCPtDCASigmaIdealPid[i] = new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
 307
 308     sprintf(name, "fTPCPtDCASigmaFullPid_%d",i);
 309         fTPCPtDCASigmaFullPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
 310
 311     sprintf(name, "fTPCPtDCASigmaDay0Pid_%d",i);
 312         fTPCPtDCASigmaDay0Pid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
 313
 314     sprintf(name, "fTPCPtDCAXYPid_%d",i);
 315         fTPCPtDCAXYPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
 316
 317     sprintf(name, "fTPCPtDCAZPid_%d",i);
 318         fTPCPtDCAZPid[i]= new TH3F(name,name,50,0.1,3,100,0,100,5,0,5);
 319   }
 320 }
 321
 322 //_____________________________________________________________________________
 323 void AliComparisonEff::InitCuts()
 324 {
 325
 326   if(!fCutsMC)
 327     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 328   if(!fCutsRC)
 329     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 330 }
 331
 332 //_____________________________________________________________________________
 333 void AliComparisonEff::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 334 {
 335   // Fill efficiency comparison information
 336
 337   AliExternalTrackParam *track = 0;
 338   Double_t kRadius    = 3.0;      // beam pipe radius
 339   Double_t kMaxStep   = 5.0;      // max step
 340   Double_t field      = 0.4;      // mag. field
 341   Double_t kMaxD      = 123456.0; // max distance
 342
 343   // systematics
 344   const Double_t kSigma2Full_xy  = 0.25; // ExB full systematics  [cm]
 345   const Double_t kSigma2Full_z  =  5.0;  // [cm]
 346
 347   const Double_t kSigma2Day0_xy  = 0.02; //  ExB  [cm]
 348   const Double_t kSigma2Day0_z  =  0.1; //  [cm] goofie
 349
 350   //
 351   Double_t      DCASigmaIdeal=0;
 352   Double_t  DCASigmaFull=0;
 353   Double_t  DCASigmaDay0=0;
 354
 355   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 356   const Double_t* dv;
 357
 358   // distance to Prim. vertex
 359   dv = infoMC->GetVDist();
 360
 361   Float_t mcpt = infoMC->GetParticle().Pt();
 362   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
 363   //Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz();
 364   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
 365
 366   // calculate and set prim. vertex
 367   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
 368   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
 369   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
 370
 371   // Check selection cuts
 372   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return;
 373
 374   // transform Pdg to Pid
 375   Double_t pid = -1;
 376   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetEP() ) pid = 0;
 377   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetMuP() ) pid = 1;
 378   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
 379   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
 380   if( TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
 381
 382   //cout << "dv[0] " << dv[0] << " dv[1] " << dv[1]  <<  " dv[2] " << dv[2] << endl;
 383   //cout << "v[0] " << fVertex->GetXv()  << " v[1] " << fVertex->GetYv()  <<  " v[2] " << fVertex->GetZv()<< endl;
 384   if (TMath::Abs(tantheta)<fCutsRC->GetMaxAbsTanTheta())
 385   {
 386
 387     // MC pt
 388     fMCPt->Fill(mcpt);
 389
 390
 391     if (infoRC->GetESDtrack()->GetTPCInnerParam())
 392     {
 393       if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
 394       {
 395         AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
 396         track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
 397
 398                 //
 399                 cov[2] = track->GetCovariance()[2];
 400
 401             // Eff = infoRC->GetStatus(1)==3 && isPrim / isPrim;
 402         // Pt vs ( dca[0]^2/cov[0]^2 + dca[1]^2/cov[2]^2 )
 403         // Pt vs ( dca[0]^2/(cov[0]^2 + kSigma2Full_xy)  + dca[1]^2/( cov[2]^2 + kSigma2Full_z )
 404         // Pt vs ( dca[0]^2/(cov[0]^2 + kSigma2_xy)  + dca[1]^2/( cov[2]^2 + kSigma2_z )
 405
 406          if(cov[0]>0.0 && cov[2]>0.0)
 407                  {
 408                          DCASigmaIdeal = TMath::Power(dca[0],2)/cov[0]
 409                                                                           + TMath::Power(dca[1],2)/cov[2];
 410
 411                          DCASigmaFull = TMath::Power(dca[0],2)/(cov[0]+kSigma2Full_xy)
 412                                                                          + TMath::Power(dca[1],2)/(cov[2]+kSigma2Full_z);
 413
 414                          DCASigmaDay0 = TMath::Power(dca[0],2)/(cov[0]+kSigma2Day0_xy)
 415                                                                          + TMath::Power(dca[1],2)/(cov[2]+kSigma2Day0_z);
 416
 417               //cout << "dca[0] " << dca[0]  << " dca[1] " << dca[1]  << endl;
 418               //cout << "cov[0] " << cov[0]  << " cov[2] " << cov[2]  << endl;
 419               //cout << "DCASigmaIdeal " << DCASigmaIdeal  << " DCASigmaFull " << DCASigmaFull  << " DCASigmaDay0 "  <<DCASigmaDay0 <<  endl;
 420             //cout << " -------------------------------------------------------- "<<  endl;
 421                  }
 422
 423          if(infoRC->GetStatus(1)==3) fMCRecPt->Fill(mcpt);
 424          if(infoRC->GetStatus(1)==3 && isPrim) fMCRecPrimPt->Fill(mcpt);
 425          if(infoRC->GetStatus(1)==3 && !isPrim) fMCRecSecPt->Fill(mcpt);
 426
 427          if(isPrim)
 428                  {
 429            fTPCPtDCASigmaIdeal[0]->Fill(mcpt,DCASigmaIdeal);
 430            fTPCPtDCASigmaFull[0]->Fill(mcpt,DCASigmaFull);
 431            fTPCPtDCASigmaDay0[0]->Fill(mcpt,DCASigmaDay0);
 432
 433            fTPCPtDCAXY[0]->Fill(mcpt,dca[0]);
 434            fTPCPtDCAZ[0]->Fill(mcpt,dca[1]);
 435
 436            fTPCPtDCASigmaIdealPid[0]->Fill(mcpt,DCASigmaIdeal,pid);
 437            fTPCPtDCASigmaFullPid[0]->Fill(mcpt,DCASigmaFull,pid);
 438            fTPCPtDCASigmaDay0Pid[0]->Fill(mcpt,DCASigmaDay0,pid);
 439
 440            fTPCPtDCAXYPid[0]->Fill(mcpt,dca[0],pid);
 441            fTPCPtDCAZPid[0]->Fill(mcpt,dca[1],pid);
 442
 443                    if(infoRC->GetStatus(1)==3)
 444                    {
 445              fTPCPtDCASigmaIdeal[1]->Fill(mcpt,DCASigmaIdeal);
 446              fTPCPtDCASigmaFull[1]->Fill(mcpt,DCASigmaFull);
 447              fTPCPtDCASigmaDay0[1]->Fill(mcpt,DCASigmaDay0);
 448
 449              fTPCPtDCAXY[1]->Fill(mcpt,dca[0]);
 450              fTPCPtDCAZ[1]->Fill(mcpt,dca[1]);
 451
 452              fTPCPtDCASigmaIdealPid[1]->Fill(mcpt,DCASigmaIdeal,pid);
 453              fTPCPtDCASigmaFullPid[1]->Fill(mcpt,DCASigmaFull,pid);
 454              fTPCPtDCASigmaDay0Pid[1]->Fill(mcpt,DCASigmaDay0,pid);
 455
 456              fTPCPtDCAXYPid[1]->Fill(mcpt,dca[0],pid);
 457              fTPCPtDCAZPid[1]->Fill(mcpt,dca[1],pid);
 458            }
 459                  }
 460
 461         // Cont = infoRC->GetStatus(1)==3 && !isPrim / infoRC->GetStatus(1)==3
 462         // Pt vs ( dz[0]^2/cov[0]^2 + dz[1]^2/cov[2]^2 )
 463         // Pt vs ( dz[0]^2/(cov[0]^2 + kSigma2Full_xy)  + dz[1]^2/( cov[2]^2 + kSigma2Full_z )
 464         // Pt vs ( dz[0]^2/(cov[0]^2 + kSigma2_xy)  + dz[1]^2/( cov[2]^2 + kSigma2_z )
 465
 466                  if(infoRC->GetStatus(1)==3)
 467                  {
 468                    fTPCPtDCASigmaIdeal[2]->Fill(mcpt,DCASigmaIdeal);
 469            fTPCPtDCASigmaFull[2]->Fill(mcpt,DCASigmaFull);
 470            fTPCPtDCASigmaDay0[2]->Fill(mcpt,DCASigmaDay0);
 471
 472            fTPCPtDCAXY[2]->Fill(mcpt,dca[0]);
 473            fTPCPtDCAZ[2]->Fill(mcpt,dca[1]);
 474
 475            fTPCPtDCASigmaIdealPid[2]->Fill(mcpt,DCASigmaIdeal,pid);
 476            fTPCPtDCASigmaFullPid[2]->Fill(mcpt,DCASigmaFull,pid);
 477            fTPCPtDCASigmaDay0Pid[2]->Fill(mcpt,DCASigmaDay0,pid);
 478
 479            fTPCPtDCAXYPid[2]->Fill(mcpt,dca[0],pid);
 480            fTPCPtDCAZPid[2]->Fill(mcpt,dca[1],pid);
 481
 482                    if(isPrim==0)
 483                    {
 484              fTPCPtDCASigmaIdeal[3]->Fill(mcpt,DCASigmaIdeal);
 485              fTPCPtDCASigmaFull[3]->Fill(mcpt,DCASigmaFull);
 486              fTPCPtDCASigmaDay0[3]->Fill(mcpt,DCASigmaDay0);
 487
 488              fTPCPtDCAXY[3]->Fill(mcpt,dca[0]);
 489              fTPCPtDCAZ[3]->Fill(mcpt,dca[1]);
 490
 491              fTPCPtDCASigmaIdealPid[3]->Fill(mcpt,DCASigmaIdeal,pid);
 492              fTPCPtDCASigmaFullPid[3]->Fill(mcpt,DCASigmaFull,pid);
 493              fTPCPtDCASigmaDay0Pid[3]->Fill(mcpt,DCASigmaDay0,pid);
 494
 495              fTPCPtDCAXYPid[3]->Fill(mcpt,dca[0],pid);
 496              fTPCPtDCAZPid[3]->Fill(mcpt,dca[1],pid);
 497            }
 498              }
 499            delete track;
 500        }
 501      }
 502          else
 503          {
 504        if(isPrim)
 505            {
 506              fTPCPtDCASigmaIdeal[0]->Fill(mcpt,0.0);
 507          fTPCPtDCASigmaFull[0]->Fill(mcpt,0.0);
 508          fTPCPtDCASigmaDay0[0]->Fill(mcpt,0.0);
 509
 510          fTPCPtDCAXY[0]->Fill(mcpt,0.0);
 511          fTPCPtDCAZ[0]->Fill(mcpt,0.0);
 512
 513          fTPCPtDCASigmaIdealPid[0]->Fill(mcpt,0.0,pid);
 514          fTPCPtDCASigmaFullPid[0]->Fill(mcpt,0.0,pid);
 515          fTPCPtDCASigmaDay0Pid[0]->Fill(mcpt,0.0,pid);
 516
 517          fTPCPtDCAXYPid[0]->Fill(mcpt,0.0,pid);
 518          fTPCPtDCAZPid[0]->Fill(mcpt,0.0,pid);
 519            }
 520      }
 521   }
 522
 523   // only primary particles
 524   if (!isPrim) return;
 525
 526   // pt
 527   if (TMath::Abs(tantheta)<fCutsRC->GetMaxAbsTanTheta()){
 528
 529     fEffTPCPt->Fill(mcpt, infoRC->GetStatus(1)==3);
 530     fEffTPCPtMC->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
 531     if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
 532       fEffTPCPtF->Fill(mcpt, infoRC->GetStatus(1)==3);
 533     }
 534
 535     // protons
 536     if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetProt()) {
 537            fEffTPCPt_P->Fill(mcpt, infoRC->GetStatus(1)==3);
 538        fEffTPCPtMC_P->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
 539
 540        if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
 541          fEffTPCPtF_P->Fill(mcpt, infoRC->GetStatus(1)==3);
 542        }
 543         }
 544
 545     // pions
 546     if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetPiP()) {
 547           fEffTPCPt_Pi->Fill(mcpt, infoRC->GetStatus(1)==3);
 548       fEffTPCPtMC_Pi->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
 549
 550        if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
 551          fEffTPCPtF_Pi->Fill(mcpt, infoRC->GetStatus(1)==3);
 552        }
 553         }
 554
 555         // kaons
 556     if(TMath::Abs(infoMC->GetParticle().GetPdgCode())==fCutsMC->GetKP()) {
 557           fEffTPCPt_K->Fill(mcpt, infoRC->GetStatus(1)==3);
 558       fEffTPCPtMC_K->Fill(mcpt, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
 559
 560        if(infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()) {
 561          fEffTPCPtF_K->Fill(mcpt, infoRC->GetStatus(1)==3);
 562        }
 563         }
 564   }
 565
 566   // theta
 567   if (TMath::Abs(mcpt)>fCutsRC->GetPtMin()){
 568     fEffTPCTan->Fill(tantheta, infoRC->GetStatus(1)==3);
 569     fEffTPCTanMC->Fill(tantheta, infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
 570     if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
 571       fEffTPCTanF->Fill(tantheta, infoRC->GetStatus(1)==3);
 572     }
 573   }
 574
 575   // pt-theta
 576   fEffTPCPtTan->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
 577   fEffTPCPtTanMC->Fill(mcpt,tantheta,infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits());
 578   if (infoMC->GetRowsWithDigits()>fCutsMC->GetMinRowsWithDigits()){
 579     fEffTPCPtTanF->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
 580   }
 581 }
 582
 583 //_____________________________________________________________________________
 584 void AliComparisonEff::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
 585 {
 586   // Process comparison information
 587   Process(infoMC,infoRC);
 588 }
 589
 590 //_____________________________________________________________________________
 591 Long64_t AliComparisonEff::Merge(TCollection* list)
 592 {
 593   // Merge list of objects (needed by PROOF)
 594
 595   if (!list)
 596   return 0;
 597
 598   if (list->IsEmpty())
 599   return 1;
 600
 601   TIterator* iter = list->MakeIterator();
 602   TObject* obj = 0;
 603
 604   // collection of generated histograms
 605
 606   Int_t count=0;
 607   while((obj = iter->Next()) != 0)
 608   {
 609     AliComparisonEff* entry = dynamic_cast<AliComparisonEff*>(obj);
 610     if (entry == 0)
 611       continue;
 612
 613     fMCPt->Add(entry->fMCPt);
 614     fMCRecPt->Add(entry->fMCRecPt);
 615     fMCRecPrimPt->Add(entry->fMCRecPrimPt);
 616     fMCRecSecPt->Add(entry->fMCRecSecPt);
 617
 618     fEffTPCPt->Add(entry->fEffTPCPt);
 619         fEffTPCPtMC->Add(entry->fEffTPCPtMC);
 620         fEffTPCPtF->Add(entry->fEffTPCPtF);
 621
 622     fEffTPCPt_P->Add(entry->fEffTPCPt_P);
 623         fEffTPCPtMC_P->Add(entry->fEffTPCPtMC_P);
 624         fEffTPCPtF_P->Add(entry->fEffTPCPtF_P);
 625
 626     fEffTPCPt_Pi->Add(entry->fEffTPCPt_Pi);
 627         fEffTPCPtMC_Pi->Add(entry->fEffTPCPtMC_Pi);
 628         fEffTPCPtF_Pi->Add(entry->fEffTPCPtF_Pi);
 629
 630     fEffTPCPt_K->Add(entry->fEffTPCPt_K);
 631         fEffTPCPtMC_K->Add(entry->fEffTPCPtMC_K);
 632         fEffTPCPtF_K->Add(entry->fEffTPCPtF_K);
 633
 634         fEffTPCTan->Add(entry->fEffTPCTan);
 635         fEffTPCTanMC->Add(entry->fEffTPCTanMC);
 636         fEffTPCTanF->Add(entry->fEffTPCTanF);
 637
 638         fEffTPCPtTan->Add(entry->fEffTPCPtTan);
 639         fEffTPCPtTanMC->Add(entry->fEffTPCPtTanMC);
 640         fEffTPCPtTanF->Add(entry->fEffTPCPtTanF);
 641
 642     for(Int_t i=0; i<4; ++i)
 643     {
 644       fTPCPtDCASigmaIdeal[i]->Add(entry->fTPCPtDCASigmaIdeal[i]);
 645       fTPCPtDCASigmaFull[i]->Add(entry->fTPCPtDCASigmaFull[i]);
 646       fTPCPtDCASigmaDay0[i]->Add(entry->fTPCPtDCASigmaDay0[i]);
 647
 648       fTPCPtDCAXY[i]->Add(entry->fTPCPtDCAXY[i]);
 649       fTPCPtDCAZ[i]->Add(entry->fTPCPtDCAXY[i]);
 650
 651       fTPCPtDCASigmaIdealPid[i]->Add(entry->fTPCPtDCASigmaIdealPid[i]);
 652       fTPCPtDCASigmaFullPid[i]->Add(entry->fTPCPtDCASigmaFullPid[i]);
 653       fTPCPtDCASigmaDay0Pid[i]->Add(entry->fTPCPtDCASigmaDay0Pid[i]);
 654
 655       fTPCPtDCAXYPid[i]->Add(entry->fTPCPtDCAXYPid[i]);
 656       fTPCPtDCAZPid[i]->Add(entry->fTPCPtDCAXYPid[i]);
 657     }
 658
 659   count++;
 660   }
 661
 662 return count;
 663 }
 664
 665 //_____________________________________________________________________________
 666 void AliComparisonEff::Analyse()
 667 {
 668   // Analyse output histograms
 669
 670   AliComparisonEff * comp=this;
 671
 672   // calculate efficiency and contamination (4 sigma)
 673
 674   TH1 *h_sigmaidealpid[20];
 675   TH1 *h_sigmafullpid[20];
 676   TH1 *h_sigmaday0pid[20];
 677
 678   //TH1 *h_sigmaday0pidclone[20];
 679
 680   TH1 *h_sigmaidealpidtot[4];
 681   TH1 *h_sigmafullpidtot[4];
 682   TH1 *h_sigmaday0pidtot[4];
 683
 684   //TH1 *h_sigmaday0pidtotclone[4];
 685
 686   char name[256];
 687   char name1[256];
 688   Int_t idx;
 689
 690   for(Int_t i=0; i<4; ++i)
 691   {
 692      //total
 693      comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange(1,4);
 694      sprintf(name,"h_sigmaidealpidtot_%d",i);
 695      h_sigmaidealpidtot[i] = comp->fTPCPtDCASigmaIdealPid[i]->Project3D();
 696      h_sigmaidealpidtot[i]->SetName(name);
 697
 698      comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange(1,4);
 699      sprintf(name,"h_sigmafullpidtot_%d",i);
 700      h_sigmafullpidtot[i] = comp->fTPCPtDCASigmaFullPid[i]->Project3D();
 701      h_sigmafullpidtot[i]->SetName(name);
 702
 703      comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange(1,4);
 704      sprintf(name,"h_sigmaday0pidtot_%d",i);
 705      h_sigmaday0pidtot[i] = comp->fTPCPtDCASigmaDay0Pid[i]->Project3D();
 706      h_sigmaday0pidtot[i]->SetName(name);
 707
 708      // pid wise
 709      for(Int_t j=0; j<5; ++j)
 710          {
 711        idx = i*5 + j;
 712
 713        comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange(1,4);
 714        comp->fTPCPtDCASigmaIdealPid[i]->GetZaxis()->SetRange(j+1,j+1);
 715
 716        sprintf(name,"h_sigmaidealpid_%d",idx);
 717        h_sigmaidealpid[idx] = comp->fTPCPtDCASigmaIdealPid[i]->Project3D();
 718        h_sigmaidealpid[idx]->SetName(name);
 719
 720
 721        comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange(1,4);
 722        comp->fTPCPtDCASigmaFullPid[i]->GetZaxis()->SetRange(j+1,j+1);
 723
 724        sprintf(name,"h_sigmafullpid_%d",idx);
 725        h_sigmafullpid[idx] = comp->fTPCPtDCASigmaFullPid[i]->Project3D();
 726        h_sigmafullpid[idx]->SetName(name);
 727
 728        comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange(1,4);
 729        comp->fTPCPtDCASigmaDay0Pid[i]->GetZaxis()->SetRange(j+1,j+1);
 730
 731        sprintf(name,"h_sigmaday0pid_%d",idx);
 732        h_sigmaday0pid[idx] = comp->fTPCPtDCASigmaDay0Pid[i]->Project3D();
 733        h_sigmaday0pid[idx]->SetName(name);
 734
 735         }
 736   }
 737
 738   // calculate efficiency and contamination (all pids)
 739   h_sigmaidealpidtot[0]->Sumw2();
 740   h_sigmaidealpidtot[1]->Divide(h_sigmaidealpidtot[0]);
 741   h_sigmaidealpidtot[2]->Sumw2();
 742   h_sigmaidealpidtot[3]->Divide(h_sigmaidealpidtot[2]);
 743
 744   h_sigmafullpidtot[0]->Sumw2();
 745   h_sigmafullpidtot[1]->Divide(h_sigmafullpidtot[0]);
 746   h_sigmafullpidtot[2]->Sumw2();
 747   h_sigmafullpidtot[3]->Divide(h_sigmafullpidtot[2]);
 748
 749   h_sigmaday0pidtot[0]->Sumw2();
 750   h_sigmaday0pidtot[1]->Divide(h_sigmaday0pidtot[0]);
 751   h_sigmaday0pidtot[2]->Sumw2();
 752   h_sigmaday0pidtot[3]->Divide(h_sigmaday0pidtot[2]);
 753
 754   // calculate efficiency pid wise
 755   for(Int_t idx = 0; idx<5; idx++)
 756   {
 757     h_sigmaidealpid[idx]->Sumw2();
 758     h_sigmaidealpid[idx+5]->Divide(h_sigmaidealpid[idx]);
 759
 760     h_sigmafullpid[idx]->Sumw2();
 761     h_sigmafullpid[idx+5]->Divide(h_sigmafullpid[idx]);
 762
 763     h_sigmaday0pid[idx]->Sumw2();
 764     h_sigmaday0pid[idx+5]->Divide(h_sigmaday0pid[idx]);
 765   }
 766
 767   // calculate cont. pid wise
 768   for(Int_t idx = 0; idx<5; idx++)
 769   {
 770     h_sigmaidealpid[idx+15]->Divide(h_sigmaidealpidtot[2]);
 771     h_sigmafullpid[idx+15]->Divide(h_sigmafullpidtot[2]);
 772     h_sigmaday0pid[idx+15]->Divide(h_sigmaday0pidtot[2]);
 773   }
 774
 775   // write results
 776   TFile *fp = new TFile("pictures_eff.root","recreate");
 777   fp->cd();
 778
 779   TCanvas * c = new TCanvas("Efficiency","Track efficiency");
 780   c->cd();
 781
 782   fMCPt->Write();
 783   fMCRecPt->Write();
 784   fMCRecPrimPt->Write();
 785   fMCRecSecPt->Write();
 786
 787   for(int i = 0; i<4;i++)
 788   {
 789    comp->fTPCPtDCASigmaIdealPid[i]->GetYaxis()->SetRange();
 790    comp->fTPCPtDCASigmaIdealPid[i]->GetZaxis()->SetRange();
 791    comp->fTPCPtDCASigmaFullPid[i]->GetYaxis()->SetRange();
 792    comp->fTPCPtDCASigmaFullPid[i]->GetZaxis()->SetRange();
 793    comp->fTPCPtDCASigmaDay0Pid[i]->GetYaxis()->SetRange();
 794    comp->fTPCPtDCASigmaDay0Pid[i]->GetZaxis()->SetRange();
 795
 796     comp->fTPCPtDCASigmaIdealPid[i]->Write();
 797     comp->fTPCPtDCASigmaFullPid[i]->Write();
 798     comp->fTPCPtDCASigmaDay0Pid[i]->Write();
 799   }
 800   //
 801   comp->fEffTPCTanF->SetXTitle("Tan(#theta)");
 802   comp->fEffTPCTanF->SetYTitle("eff_{findable}");
 803   comp->fEffTPCTanF->Write("EffTanFindable");
 804   //
 805   comp->fEffTPCTan->SetXTitle("Tan(#theta)");
 806   comp->fEffTPCTan->SetYTitle("eff_{all}");
 807   comp->fEffTPCTan->Write("EffTanAll");
 808
 809   h_sigmaidealpidtot[1]->Write("Eff_SigmaIdeal");
 810   h_sigmaidealpidtot[3]->Write("Cont_SigmaIdeal");
 811
 812   h_sigmafullpidtot[1]->Write("Eff_SigmaFull");
 813   h_sigmafullpidtot[3]->Write("Cont_SigmaFull");
 814
 815   h_sigmaday0pidtot[1]->Write("Eff_SigmaDay0");
 816   h_sigmaday0pidtot[3]->Write("Cont_SigmaDay0");
 817
 818   for(Int_t idx = 0; idx<5; idx++)
 819   {
 820     sprintf(name,"Eff_SigmaIdeal_%d",idx);
 821     sprintf(name1,"Cont_SigmaIdeal_%d",idx);
 822
 823     h_sigmaidealpid[idx+5]->Write(name);
 824     h_sigmaidealpid[idx+15]->Write(name1);
 825
 826     sprintf(name,"Eff_SigmaFull_%d",idx);
 827     sprintf(name1,"Cont_SigmaFull_%d",idx);
 828
 829     h_sigmafullpid[idx+5]->Write(name);
 830     h_sigmafullpid[idx+15]->Write(name1);
 831
 832     sprintf(name,"Eff_SigmaDay0_%d",idx);
 833     sprintf(name1,"Cont_SigmaDay0_%d",idx);
 834
 835     h_sigmaday0pid[idx+5]->Write(name);
 836     h_sigmaday0pid[idx+15]->Write(name1);
 837   }
 838
 839   //
 840   fp->Close();
 841 }