test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TProfile.h"
   6
   7 #include "TCanvas.h"
   8 #include "TList.h"
   9 #include "TParticle.h"
  10 #include "TParticlePDG.h"
  11 #include "TProfile.h"
  12
  13 #include "AliAnalysisTask.h"
  14 #include "AliAnalysisManager.h"
  15
  16 #include "AliESDEvent.h"
  17 #include "AliStack.h"
  18 #include "AliESDVertex.h"
  19 #include "AliESDInputHandler.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliMultiplicity.h"
  22
  23 #include "AliMCParticle.h"
  24 #include "AliMCEvent.h"
  25 #include "AliAnalysisTaskEfficiency.h"
  26 #include "AliExternalTrackParam.h"
  27 #include "AliTrackReference.h"
  28 #include "AliStack.h"
  29 #include "AliHeader.h"
  30 #include "AliGenEventHeader.h"
  31 #include "AliGenCocktailEventHeader.h"
  32 #include "AliGenPythiaEventHeader.h"
  33 #include "AliGenDPMjetEventHeader.h"
  34
  35 // Analysis Task for basic QA on the ESD
  36 // Authors: Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing, Veronica Canoa, AM, Eva Sicking
  37
  38 ClassImp(AliAnalysisTaskEfficiency)
  39
  40 //________________________________________________________________________
  41   AliAnalysisTaskEfficiency::AliAnalysisTaskEfficiency(const char *name)
  42     : AliAnalysisTaskSE(name)
  43     ,fFieldOn(kTRUE)
  44     ,fHists(0)
  45     ,fHistRECpt(0)
  46     ,fHistMCpt(0)
  47     ,fHistMCNRpt(0)
  48     ,fHistFAKEpt(0)
  49     ,fHistMULTpt(0)
  50     ,fHistRECeta(0)
  51     ,fHistMCeta(0)
  52     ,fHistMCNReta(0)
  53     ,fHistFAKEeta(0)
  54     ,fHistMULTeta(0)
  55     ,fHistRECphi(0)
  56     ,fHistMCphi(0)
  57     ,fHistMCNRphi(0)
  58     ,fHistFAKEphi(0)
  59     ,fHistMULTphi(0)
  60     ,fHistRECHPTeta(0)
  61     ,fHistMCHPTeta(0)
  62     ,fHistMCNRHPTeta(0)
  63     ,fHistFAKEHPTeta(0)
  64     ,fHistRECHPTphi(0)
  65     ,fHistMCHPTphi(0)
  66     ,fHistMCNRHPTphi(0)
  67     ,fHistFAKEHPTphi(0)
  68     ,fHistRecMult(0)
  69     ,fHistNCluster(0)
  70     ,fh1VtxEff(0)
  71     ,fh2VtxTpcSpd(0)
  72     ,fh2EtaCorrelation(0)
  73     ,fh2EtaCorrelationShift(0)
  74     ,fh2PhiCorrelation(0)
  75     ,fh2PhiCorrelationShift(0)
  76     ,fh2PtCorrelation(0)
  77     ,fh2PtCorrelationShift(0)
  78     ,fHistDeltaphiprimaries(0)
  79     ,fHistDeltaphisecondaries(0)
  80     ,fHistDeltaphireject(0)
  81     ,fHistTPCITSdeltax(0)
  82     ,fHistTPCITSdeltay(0)
  83     ,fHistTPCITSdeltaz(0)
  84     ,fHistTPCITSdeltar(0)
  85     ,fHistTPCTRDdeltax(0)
  86     ,fHistTPCTRDdeltay(0)
  87     ,fHistTPCTRDdeltaz(0)
  88     ,fHistTPCTRDdeltar(0)
  89     ,fHistTPCTRDdeltaphi(0)
  90     ,fPtPullsPos(0)
  91     ,fPtPullsNeg(0)
  92     ,fPtShiftPos(0)
  93     ,fPtShiftNeg(0)
  94     ,fTrackType(0)
  95     ,fCuts(0)
  96     ,fVertexRvsZ(0)
  97     ,fVertexRvsZC(0)
  98     ,fEtaMultiFluc(0)
  99     ,fEtaMulti(0)
 100     ,fEtaMultiH(0)
 101     ,fPhiMultiFluc(0)
 102     ,fPhiMulti(0)
 103     ,fPhiMultiH(0)
 104 {
 105   //
 106   // Constructor
 107   //
 108
 109   for(int i = 0;i< 2;++i){
 110     fh2MultSpdChips[i] = 0;
 111   }
 112
 113   for(int i = 0;i < kTotalAC;++i){
 114     fh2PhiPadRow[i] = fh2PhiLayer[i] = 0;
 115   }
 116
 117   for(int i = 0; i < kTotalVtx; ++i){
 118     fh2VertexCorrelation[i] =
 119       fh2VertexCorrelationShift[i] = 0;
 120     fh1VertexShift[i] =
 121       fh1VertexShiftNorm[i] = 0;
 122   }
 123
 124   for(int i = 0;i< 8;++i){
 125     fHistRECptCharge[i]=0;
 126     fHistMCptCharge[i]=0;
 127     fHistMCNRptCharge[i]=0;
 128     fHistFAKEptCharge[i]=0;
 129
 130     fHistRECetaCharge[i]=0;
 131     fHistMCetaCharge[i]=0;
 132     fHistMCNRetaCharge[i]=0;
 133     fHistFAKEetaCharge[i]=0;
 134
 135     fHistRECphiCharge[i]=0;
 136     fHistMCphiCharge[i]=0;
 137     fHistMCNRphiCharge[i]=0;
 138     fHistFAKEphiCharge[i]=0;
 139
 140   }
 141
 142   DefineOutput(1,  TList::Class());
 143 }
 144
 145
 146 //________________________________________________________________________
 147 void AliAnalysisTaskEfficiency::UserCreateOutputObjects()
 148 {
 149   // Create histograms
 150
 151
 152   fHists = new TList();
 153   fHistRECpt   = new TH1F("fHistRECpt",  " p_{T} distribution: all reconstructed",        100, 0., 20.);
 154   fHistMCpt    = new TH1F("fHistMCpt",   " p_{T} distribution: all MC",                   100, 0., 20.);
 155   fHistMCNRpt  = new TH1F("fHistMCNRpt", " p_{T} distribution: all not-reconstructed MC", 100, 0., 20.);
 156   fHistFAKEpt  = new TH1F("fHistFAKEpt", " p_{T} distribution: all Fake",                 100, 0., 20.);
 157   fHistMULTpt  = new TH1F("fHistMULTpt", " p_{T} distribution: multiply rec.",            100, 0., 20.);
 158
 159   fHistRECeta   = new TH1F("fHistRECeta",  " #eta distribution: all reconstructed",        100,-1.0,1.0);
 160   fHistMCeta    = new TH1F("fHistMCeta",   " #eta distribution: all MC",                   100,-1.0,1.0);
 161   fHistMCNReta  = new TH1F("fHistMCNReta", " #eta distribution: all not-reconstructed MC", 100,-1.0,1.0);
 162   fHistFAKEeta  = new TH1F("fHistFAKEeta", " #eta distribution: all Fake",                 100,-1.0,1.0);
 163   fHistMULTeta  = new TH1F("fHistMULTeta", " #eta distribution: multiply rec.",            100,-1.0,1.0);
 164
 165   fHistRECphi   = new TH1F("fHistRECphi",  " #phi distribution: all reconstructed",        314, 0., 6.28);
 166   fHistMCphi    = new TH1F("fHistMCphi",   " #phi distribution: all MC",                   314, 0., 6.28);
 167   fHistMCNRphi  = new TH1F("fHistMCNRphi", " #phi distribution: all not-reconstructed MC", 314, 0., 6.28);
 168   fHistFAKEphi  = new TH1F("fHistFAKEphi", " #phi distribution: all Fake",                 314, 0., 6.28);
 169   fHistMULTphi  = new TH1F("fHistMULTphi", " #phi distribution: multipli rec.",            314, 0., 6.28);
 170
 171   fHistRECHPTeta   = new TH1F("fHistRECHPTeta",  " #eta distribution: all reconstructed",        100,-1.0,1.0);
 172   fHistMCHPTeta    = new TH1F("fHistMCHPTeta",   " #eta distribution: all MC",                   100,-1.0,1.0);
 173   fHistMCNRHPTeta  = new TH1F("fHistMCNRHPTeta", " #eta distribution: all not-reconstructed MC", 100,-1.0,1.0);
 174   fHistFAKEHPTeta  = new TH1F("fHistFAKEHPTeta", " #eta distribution: all Fake",                 100,-1.0,1.0);
 175
 176   fHistRECHPTphi   = new TH1F("fHistRECHPTphi",  " #phi distribution: all reconstructed",        314, 0., 6.28);
 177   fHistMCHPTphi    = new TH1F("fHistMCHPTphi",   " #phi distribution: all MC",                   314, 0., 6.28);
 178   fHistMCNRHPTphi  = new TH1F("fHistMCNRHPTphi", " #phi distribution: all not-reconstructed MC", 314, 0., 6.28);
 179   fHistFAKEHPTphi  = new TH1F("fHistFAKEHPTphi", " #phi distribution: all Fake",                 314, 0., 6.28);
 180
 181   fHistRecMult     = new TH1F("fHistRecMult",  "Multiple reconstructed tracks", 50, 0., 50.);
 182   fHistNCluster    = new TH1F("fHistNCluster", "Number of clusters for suspicious tracks", 300, 0., 300.);
 183
 184   // CKB
 185
 186   fh1VtxEff = new TH1F("fh1VtxEff","VtxEff TPC",4,-0.5,3.5);
 187   fh1VtxEff->GetXaxis()->SetBinLabel(1,"NO TPC VTX");
 188   fh1VtxEff->GetXaxis()->SetBinLabel(2,"TPC VTX");
 189   fh1VtxEff->GetXaxis()->SetBinLabel(3,"NO SPD VTX");
 190   fh1VtxEff->GetXaxis()->SetBinLabel(4,"SPD VTX");
 191
 192   TString labels[kTotalAC];
 193   labels[kNegA]="NegA";
 194   labels[kPosA]="PosA";
 195   labels[kNegC]="NegC";
 196   labels[kPosC]="PosC";
 197
 198   for(int i = 0;i< kTotalAC;++i){
 199     fh2PhiPadRow[i] = new TH2F(Form("fh2PhiPadRow_%d",i),Form("Padrow vs phi AC %s;phi;padrow",labels[i].Data()),360,0.,360.,159,-0.5,158.5);
 200     fh2PhiLayer[i] = new TH2F(Form("fh2PhiLayer_%d",i),Form("layer vs phi AC %s;phi;layer",labels[i].Data()),360,0.,360.,6,-0.5,5.5);
 201   }
 202   for(int i = 0;i<2;++i){
 203     fh2MultSpdChips[i] = new TH2F(Form("fh2ChipsSpdMult_%d",i),"mult SPD vs chips;chips;Mult",201,-1.5,199.5,201,-1.5,199.5);
 204   }
 205   fh2VtxTpcSpd  = new TH2F("fh2VtxTpcSpd","SPD vtx vs TPC;tpc-z;spd-z",200,-20.,20.,200,-20,20.);
 206
 207   TString labelsVtx[kTotalVtx];
 208   labelsVtx[kTPC] = "TPC";
 209   labelsVtx[kSPD] = "SPD";
 210   for(int i = 0; i < kTotalVtx; ++i){
 211     fh2VertexCorrelation[i] = new TH2F(Form("fh2VertexCorrelation_%d",i),Form("VertexCorrelation %s;MC z-vtx;ESD z-vtx",labelsVtx[i].Data()), 120, -30, 30, 120, -30, 30);
 212     fh2VertexCorrelationShift[i] = new TH2F(Form("fh2VertexCorrelationShift_%d",i), Form("VertexCorrelationShift %s;MC z-vtx;MC z-vtx - ESD z-vtx",labelsVtx[i].Data()), 120, -30, 30, 100, -1, 1);
 213     fh1VertexShift[i] =  new TH1F(Form("fh1VertexShift_%d",i), Form("VertexShift %s;(MC z-vtx - ESD z-vtx);Entries",labelsVtx[i].Data()), 201, -2, 2);
 214     fh1VertexShiftNorm[i] = new TH1F(Form("fh1VertexShiftNorm_%d",i), Form("VertexShiftNorm %s;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries",labelsVtx[i].Data()), 200, -100, 100);
 215   }
 216
 217
 218   fh2EtaCorrelation = new TH2F("fh2EtaCorrelation", "EtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
 219   fh2EtaCorrelationShift = new TH2F("fh2EtaCorrelationShift", "EtaCorrelationShift;ESD #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);
 220   fh2PhiCorrelation = new TH2F("fh2PhiCorrelation", "PhiCorrelation;MC #phi;ESD #phi", 120, 0, 360, 120, 0, 360);
 221   fh2PhiCorrelationShift = new TH2F("fh2PhiCorrelationShift", "PhiCorrelationShift;ESD #phi;MC #phi - ESD #phi", 120, 0, 360, 100, -10, 10);
 222   fh2PtCorrelation = new TH2F("fh2PtCorrelation", "PtCorrelation;MC p_T;ESD p_T", 200, 0., 200., 200, 0, 200);
 223   fh2PtCorrelationShift = new TH2F("fh2PtCorrelationShift", "PtCorrelationShift;ESD p_T;MC p_T - ESD #p_T", 120, 0, 10, 100, -2, 2);
 224
 225   fHistDeltaphiprimaries  =new TH1F("fHistDeltaphiprimaries", " #Delta#phi distribution: primaries",     314, -0.2,0.2);
 226   fHistDeltaphisecondaries=new TH1F("fHistDeltaphisecondaries", "#Delta#phi distribution: secondaries", 314, -0.2,0.2);
 227   fHistDeltaphireject     =new TH1F("fHistDeltaphireject", " #Delta#phi distribution: reject trackled",  314, -0.2,0.2);
 228
 229   fHistTPCITSdeltax =new TH1F("fHistTPCITSdeltax", "TPC-ITS matching:#Delta x", 100,-20,20);
 230   fHistTPCITSdeltay =new TH1F("fHistTPCITSdeltay", "TPC-ITS matching:#Delta y", 100,-20,20);
 231   fHistTPCITSdeltaz =new TH1F("fHistTPCITSdeltaz", "TPC-ITS matching:#Delta z", 100,-20,20);
 232   fHistTPCITSdeltar =new TH1F("fHistTPCITSdeltar", "TPC-ITS matching:#Delta r", 100,0,20);
 233
 234   fHistTPCTRDdeltax   = new TH1F("fHistTPCTRDdeltax", "TPC-TRD matching:#Delta x", 400,-400, 400);
 235   fHistTPCTRDdeltay   = new TH1F("fHistTPCTRDdeltay", "TPC-TRD matching:#Delta y", 400,-400, 400);
 236   fHistTPCTRDdeltaz   = new TH1F("fHistTPCTRDdeltaz", "TPC-TRD matching:#Delta z", 400,-400, 400);
 237   fHistTPCTRDdeltar   = new TH1F("fHistTPCTRDdeltar", "TPC-TRD matching:#Delta r", 100,0,20);
 238   fHistTPCTRDdeltaphi = new TH1F("fHistTPCTRDdeltarphi", "TPC-TRD matching:#Delta #phi", 128,-3.14 , 3.14);
 239   // Pulls
 240   fPtPullsPos = new TProfile("fPtPullsPos", "Pt Pulls for pos. primaries", 50 , 0., 10., -5., 5.,"S");
 241   fPtPullsNeg = new TProfile("fPtPullsNeg", "Pt Pulls for neg. primaries", 50 , 0., 10., -5., 5.,"S");
 242   fPtShiftPos = new TProfile("fPtShiftPos", "Pt Shift for pos. primaries", 50 , 0., 10., -0.2, 0.2,"S");
 243   fPtShiftNeg = new TProfile("fPtShiftNeg", "Pt Shift for neg. primaries", 50 , 0., 10., -0.2, 0.2,"S");
 244
 245   fEtaMultiFluc = new TProfile("fEtaMultiFluc", "eta multiplicity fluctuations", 10, -2., 2., 0., 600., "S");
 246   fEtaMulti     = new TH1F("fEtaMulti",     "eta multiplicity fluctuations", 10, -2., 2.);
 247   fEtaMultiH    = new TH1F("fEtaMultiH",    "eta multiplicity fluctuations", 600,  0., 600.);
 248
 249   fPhiMultiFluc = new TProfile("fPhiMultiFluc", "phi multiplicity fluctuations", 64, 0., 6.4, 0., 100., "S");
 250   fPhiMulti     = new TH1F("fPhiMulti",         "phi multiplicity fluctuations", 64, 0., 6.4);
 251   fPhiMultiH    = new TH1F("fPhiMultiH",        "phi multiplicity fluctuations", 100, 0., 100.);
 252
 253   // SPD Vertex
 254   fVertexRvsZ  = new TH2F("fVertexRvsZ", "SPD Vertex R vs Z", 200, -1., 1., 200, -1., 1.);
 255   fVertexRvsZC = new TH2F("fVertexRvsZC", "SPD Vertex R vs Z", 200, -1., 1., 200, -1., 1.);
 256
 257   TString charge[8]={"Deu", "AntiDeu", "Tri", "AntiTri", "He3", "AntiHe3", "He4", "AntiHe4"};
 258   for(Int_t i=0;i<8;i++){
 259     fHistRECptCharge[i]   = new TH1F(Form("fHistRECptCharge%s",charge[i].Data()),
 260                                      "p_{T} distribution: all reconstructed",
 261                                      100, 0., 20.);
 262     fHistMCptCharge[i]    = new TH1F(Form("fHistMCptCharge%s",charge[i].Data()),
 263                                      "p_{T} distribution: all MC",
 264                                      100, 0., 20.);
 265     fHistMCNRptCharge[i]  = new TH1F(Form("fHistMCNRptCharge%s",charge[i].Data()),
 266                                      "p_{T} distribution: all not-reconstructed MC",
 267                                      100, 0., 20.);
 268     fHistFAKEptCharge[i]  = new TH1F( Form("fHistFAKEptCharge%s",charge[i].Data()),
 269                                      "p_{T} distribution: all Fake",
 270                                      100, 0., 20.);
 271
 272     fHistRECetaCharge[i]   = new TH1F(Form("fHistRECetaCharge%s",charge[i].Data()),
 273                                      "p_{T} distribution: all reconstructed",
 274                                      100, 0., 20.);
 275     fHistMCetaCharge[i]    = new TH1F(Form("fHistMCetaCharge%s",charge[i].Data()),
 276                                      "p_{T} distribution: all MC",
 277                                      100, 0., 20.);
 278     fHistMCNRetaCharge[i]  = new TH1F(Form("fHistMCNRetaCharge%s",charge[i].Data()),
 279                                      "p_{T} distribution: all not-reconstructed MC",
 280                                      100, 0., 20.);
 281     fHistFAKEetaCharge[i]  = new TH1F( Form("fHistFAKEetaCharge%s",charge[i].Data()),
 282                                      "p_{T} distribution: all Fake",
 283                                      100, 0., 20.);
 284
 285     fHistRECphiCharge[i]   = new TH1F(Form("fHistRECphiCharge%s",charge[i].Data()),
 286                                      "p_{T} distribution: all reconstructed",
 287                                      100, 0., 20.);
 288     fHistMCphiCharge[i]    = new TH1F(Form("fHistMCphiCharge%s",charge[i].Data()),
 289                                      "p_{T} distribution: all MC",
 290                                      100, 0., 20.);
 291     fHistMCNRphiCharge[i]  = new TH1F(Form("fHistMCNRphiCharge%s",charge[i].Data()),
 292                                      "p_{T} distribution: all not-reconstructed MC",
 293                                      100, 0., 20.);
 294     fHistFAKEphiCharge[i]  = new TH1F( Form("fHistFAKEphiCharge%s",charge[i].Data()),
 295                                      "p_{T} distribution: all Fake",
 296                                      100, 0., 20.);
 297
 298   }
 299   // BKC
 300
 301
 302   fHists->SetOwner();
 303
 304   fHists->Add(fHistRECpt);
 305   fHists->Add(fHistMCpt);
 306   fHists->Add(fHistMCNRpt);
 307   fHists->Add(fHistFAKEpt);
 308   fHists->Add(fHistMULTpt);
 309
 310   fHists->Add(fHistRECeta);
 311   fHists->Add(fHistMCeta);
 312   fHists->Add(fHistMCNReta);
 313   fHists->Add(fHistFAKEeta);
 314   fHists->Add(fHistMULTeta);
 315
 316   fHists->Add(fHistRECphi);
 317   fHists->Add(fHistMCphi);
 318   fHists->Add(fHistMCNRphi);
 319   fHists->Add(fHistFAKEphi);
 320   fHists->Add(fHistMULTphi);
 321
 322   fHists->Add(fHistRECHPTeta);
 323   fHists->Add(fHistMCHPTeta);
 324   fHists->Add(fHistMCNRHPTeta);
 325   fHists->Add(fHistFAKEHPTeta);
 326
 327   fHists->Add(fHistRECHPTphi);
 328   fHists->Add(fHistMCHPTphi);
 329   fHists->Add(fHistMCNRHPTphi);
 330   fHists->Add(fHistFAKEHPTphi);
 331
 332   // CKB
 333
 334   fHists->Add(fh1VtxEff);
 335   for(int i = 0;i < kTotalAC;++i){
 336     fHists->Add(fh2PhiPadRow[i]);
 337     fHists->Add(fh2PhiLayer[i]);
 338   }
 339   for(int i = 0;i<2;++i){
 340     fHists->Add(fh2MultSpdChips[i]);
 341   }
 342   fHists->Add(fh2VtxTpcSpd);
 343
 344   for(int i = 0; i < kTotalVtx; ++i){
 345     fHists->Add(fh2VertexCorrelation[i]);
 346     fHists->Add(fh2VertexCorrelationShift[i]);
 347     fHists->Add(fh1VertexShift[i]);
 348     fHists->Add(fh1VertexShiftNorm[i]);
 349   }
 350
 351
 352   fHists->Add(fh2EtaCorrelation);
 353   fHists->Add(fh2EtaCorrelationShift);
 354   fHists->Add(fh2PhiCorrelation);
 355   fHists->Add(fh2PhiCorrelationShift);
 356   fHists->Add(fh2PtCorrelation);
 357   fHists->Add(fh2PtCorrelationShift);
 358
 359   fHists->Add(fHistDeltaphiprimaries);
 360   fHists->Add(fHistDeltaphisecondaries);
 361   fHists->Add(fHistDeltaphireject);
 362
 363   fHists->Add(fHistTPCITSdeltax);
 364   fHists->Add(fHistTPCITSdeltay);
 365   fHists->Add(fHistTPCITSdeltaz);
 366   fHists->Add(fHistTPCITSdeltar);
 367
 368   fHists->Add(fHistTPCTRDdeltax);
 369   fHists->Add(fHistTPCTRDdeltay);
 370   fHists->Add(fHistTPCTRDdeltaz);
 371   fHists->Add(fHistTPCTRDdeltar);
 372   fHists->Add(fHistTPCTRDdeltaphi);
 373
 374   fHists->Add(fHistRecMult);
 375   fHists->Add(fHistNCluster);
 376
 377
 378   fHists->Add(fPtPullsPos);
 379   fHists->Add(fPtPullsNeg);
 380   fHists->Add(fPtShiftPos);
 381   fHists->Add(fPtShiftNeg);
 382
 383   fHists->Add(fVertexRvsZ);
 384   fHists->Add(fVertexRvsZC);
 385   fHists->Add(fEtaMultiFluc);
 386   fHists->Add(fEtaMulti);
 387   fHists->Add(fEtaMultiH);
 388   fHists->Add(fPhiMultiFluc);
 389   fHists->Add(fPhiMulti);
 390   fHists->Add(fPhiMultiH);
 391
 392   for(Int_t i=0;i<8;i++){
 393     fHists->Add(fHistRECptCharge[i]);
 394     fHists->Add(fHistMCptCharge[i]);
 395     fHists->Add(fHistMCNRptCharge[i]);
 396     fHists->Add(fHistFAKEptCharge[i]);
 397
 398     fHists->Add(fHistRECetaCharge[i]);
 399     fHists->Add(fHistMCetaCharge[i]);
 400     fHists->Add(fHistMCNRetaCharge[i]);
 401     fHists->Add(fHistFAKEetaCharge[i]);
 402
 403     fHists->Add(fHistRECphiCharge[i]);
 404     fHists->Add(fHistMCphiCharge[i]);
 405     fHists->Add(fHistMCNRphiCharge[i]);
 406     fHists->Add(fHistFAKEphiCharge[i]);
 407   }
 408
 409   for (Int_t i=0; i<fHists->GetEntries(); ++i) {
 410     TH1 *h1 = dynamic_cast<TH1*>(fHists->At(i));
 411     if (h1){
 412       h1->Sumw2();
 413     }
 414   }
 415   // BKC
 416
 417
 418   // Post output data.
 419   PostData(1, fHists);
 420
 421 }
 422
 423 //________________________________________________________________________
 424 void AliAnalysisTaskEfficiency::UserExec(Option_t *)
 425 {
 426   // Event loop
 427   // Analysis of MC true particles and reconstructed tracks
 428   // Different track types (Global, TPC, ITS) can be selected
 429
 430   if (!fInputEvent) {
 431     Printf("ERROR: fESD not available");
 432     return;
 433   }
 434   static Int_t nfc = 0;
 435
 436   //AliESDInputHandler* esdI = (AliESDInputHandler*)
 437   //(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler();
 438   //AliESDEvent* hltEvent = esdI->GetHLTEvent();
 439
 440   AliStack* stack = MCEvent()->Stack();
 441   AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
 442
 443   // CKB
 444   // Fetch the SPD vertex:
 445   Double_t vtxSPD[3];
 446   const AliESDVertex* spdVertex = esdE->GetPrimaryVertexSPD();
 447   if (spdVertex) {
 448     if(spdVertex->GetNContributors()<=0){
 449       spdVertex = 0;
 450     }
 451     else{
 452       spdVertex->GetXYZ(vtxSPD);
 453     }
 454   }
 455
 456   // Fetch the TPC vertex
 457   Double_t vtxTPC[3];
 458   const AliESDVertex* tpcVertex = esdE->GetPrimaryVertexTPC();
 459   if (tpcVertex) {
 460     if(tpcVertex->GetNContributors()<=0){
 461       tpcVertex = 0;
 462     }
 463     else{
 464       tpcVertex->GetXYZ(vtxTPC);
 465     }
 466   }
 467   // SPD Vertex
 468   if (spdVertex) {
 469     Double_t x,y,z;
 470     x = spdVertex->GetX() + 0.07;
 471     y = spdVertex->GetY() - 0.25;
 472     z = spdVertex->GetZ();
 473     if (TMath::Abs(z) < 10.) {
 474       fVertexRvsZ->Fill(x,y);
 475       if (TMath::Sqrt(x*x + y*y) > 5. * 0.028)
 476         fVertexRvsZC->Fill(x,y);
 477     }
 478   }
 479   // BKC
 480   //Printf("%s:%d %5d",(char*)__FILE__,__LINE__, esdE->GetNumberOfTracks());
 481   TArrayI labels(esdE->GetNumberOfTracks());
 482   Int_t igood = 0;
 483   // Track loop to fill a pT spectrum
 484   //printf("ESD track loop \n");
 485
 486   AliESDtrack *tpcP = 0x0;
 487
 488   for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
 489
 490     //prevent mem leak for TPConly track
 491     if(fTrackType==2&&tpcP){
 492       delete tpcP;
 493       tpcP = 0;
 494     }
 495
 496     AliVParticle *track = esdE->GetTrack(iTracks);
 497     AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
 498     esdtrack->PropagateToDCA(esdE->GetPrimaryVertex(),
 499                              esdE->GetMagneticField(), 10000.);
 500
 501     if (!track) {
 502       Printf("ERROR: Could not receive track %d", iTracks);
 503       continue;
 504     }
 505     //__________
 506     // run Task for global tracks or ITS tracks or TPC tracks
 507     const AliExternalTrackParam *tpcPin = 0x0;
 508     //Double_t phiIn=0.;
 509     Float_t phi = 0.;
 510     Float_t eta = 0.;
 511     Float_t pt  = 0.;
 512     Int_t indexAC = GetIndexAC(esdtrack);
 513
 514     //select in the steering macro, which track type should be analysed.
 515     // Four track types are selectable:
 516     // 0. Global tracks
 517     // 1. ITS tracks (SA or Pure SA)
 518     // 2. TPC only tracks
 519     // Track selection copied from PWGPP/AliAnalysisTaskQAsym
 520
 521     if(fTrackType==0){
 522       //Fill all histograms with global tracks
 523       tpcP = esdtrack;
 524       if (!tpcP) continue;
 525       if (!fCuts->AcceptTrack(tpcP)) continue;
 526       phi = tpcP->Phi();
 527       eta = tpcP->Eta();
 528       pt  = tpcP->Pt();
 529     }
 530     else if(fTrackType==1){
 531       //Fill all histograms with ITS tracks
 532       tpcP = esdtrack;
 533       if (!tpcP) continue;
 534       if (!fCuts->AcceptTrack(tpcP)) continue;
 535       phi = tpcP->Phi();
 536       eta = tpcP->Eta();
 537       pt  = tpcP->Pt();
 538     }
 539     else if(fTrackType==2){
 540       //Fill all histograms with TPC track information
 541       tpcPin = esdtrack->GetInnerParam();
 542       if (!tpcPin) continue;
 543       tpcP = AliESDtrackCuts::
 544         GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdE),
 545                         esdtrack->GetID());
 546       if (!tpcP) continue;
 547       if (!fCuts->AcceptTrack(tpcP)) continue;
 548       if(tpcP->GetNcls(1)>160)continue;//jacek's track cut
 549       if(tpcP->GetConstrainedChi2TPC()<0)continue; // jacek's track cut
 550       phi=tpcPin->Phi();
 551       eta=tpcPin->Eta();
 552       pt=tpcPin->Pt();
 553     }
 554     else{
 555       Printf("ERROR: wrong track type \n");
 556       continue;
 557     }
 558     //___________
 559     //
 560
 561
 562     labels.AddAt(-1, iTracks);
 563
 564
 565       // Selected
 566       if (TMath::Abs(eta) > 0.9) {
 567       } else {
 568
 569         //Int_t charge=track->Charge();
 570         // Within acceptance
 571         Int_t ind = TMath::Abs(esdtrack->GetLabel());
 572         AliMCParticle* mcPtest = (AliMCParticle*) MCEvent()->GetTrack(ind);
 573         if (stack->IsPhysicalPrimary(mcPtest->Label())){
 574           if(TMath::Abs(mcPtest->PdgCode())>=1000020030){//all helium (+-1000020030,+-1000020040)
 575             pt*=2;// reconstruction takes charge=1 -> for particles with charge=2, pT,rec need to be corrected
 576           }
 577
 578           Int_t index=ConvertePDG(mcPtest->PdgCode());
 579           if(fDebug>1)Printf("PDG=%d, index=%d", mcPtest->PdgCode(), index);
 580           //fill tracks comming from d,t,3He, 4He
 581           if(index<8){
 582             fHistRECptCharge [index] ->Fill(pt);
 583             fHistRECetaCharge[index]->Fill(eta);
 584             fHistRECphiCharge[index]->Fill(phi);
 585           }
 586
 587         }
 588         fHistRECpt ->Fill(pt);
 589         fHistRECeta->Fill(eta);
 590         fHistRECphi->Fill(phi);
 591
 592         if (pt > 2.) {
 593           fHistRECHPTeta->Fill(eta);
 594           fHistRECHPTphi->Fill(phi);
 595         }
 596
 597         //      Int_t ind = TMath::Abs(esdtrack->GetLabel());
 598         AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(ind);
 599         if (!(stack->IsPhysicalPrimary(mcP->Label()))) {
 600           fHistFAKEpt ->Fill(pt);
 601           fHistFAKEeta->Fill(eta);
 602           fHistFAKEphi->Fill(phi);
 603
 604           if (pt > 2.) {
 605             fHistFAKEHPTeta->Fill(eta);
 606             fHistFAKEHPTphi->Fill(phi);
 607           }
 608
 609         }
 610
 611         labels.AddAt(TMath::Abs(esdtrack->GetLabel()), iTracks);
 612         igood++;
 613
 614
 615         Float_t phiDeg = TMath::RadToDeg()*track->Phi();
 616
 617         //TPC-ITS matching
 618         //TPC-TRD matching
 619
 620         if (tpcP){
 621           Int_t labeltpcits = esdtrack->GetLabel();
 622           AliExternalTrackParam * tpcPCopy = new AliExternalTrackParam(*tpcP);
 623           Double_t xk = 43.6;  // cm
 624           Double_t bz = 5.0;   // kG
 625           if(tpcPCopy->PropagateTo(xk,bz)){
 626             Double_t alpha=tpcPCopy->GetAlpha();
 627             // if(tpcPCopy->Rotate(0.)){
 628             Float_t xtpc=tpcPCopy->GetX();
 629             Float_t ytpc=tpcPCopy->GetY();
 630             Float_t ztpc=tpcPCopy->GetZ();
 631             Float_t xpr,ypr ;
 632             xpr=xtpc*TMath::Cos(alpha)-ytpc*TMath::Sin(alpha);
 633             ypr=xtpc*TMath::Sin(alpha)+ytpc*TMath::Cos(alpha);
 634
 635             AliMCParticle* mcPart = (AliMCParticle*) MCEvent()->GetTrack(abs(labeltpcits));
 636             Int_t ntref = mcPart->GetNumberOfTrackReferences();
 637
 638             for (Int_t k = 0; k < ntref; k++) {
 639               AliTrackReference* ref = mcPart->GetTrackReference(k);
 640               Float_t xhits = ref->X();
 641               Float_t yhits = ref->Y();
 642               Float_t radio = TMath::Sqrt(xhits*xhits+yhits*yhits);
 643
 644               if(ref->DetectorId() == 0 && (radio > 42.8 && radio < 43.7)) {
 645
 646                 Float_t xits = ref->X();
 647                 Float_t yits = ref->Y();
 648                 Float_t zits = ref->Z();
 649
 650                 Float_t difx=(xits-xpr);
 651                 fHistTPCITSdeltax->Fill(difx);
 652                 Float_t dify=(yits-ypr);
 653                 fHistTPCITSdeltay->Fill(dify);
 654                 Float_t difz=(zits-ztpc);
 655                 fHistTPCITSdeltaz->Fill(difz);
 656                 Float_t deltar = TMath::Sqrt(difx * difx + dify * dify + difz * difz);
 657                 fHistTPCITSdeltar->Fill(deltar);
 658                 break;
 659               }
 660               // }
 661             } // trackrefs
 662           }
 663         } //ITS-TPC maching
 664           //TPC-TRD maching
 665
 666         const AliExternalTrackParam *trd = esdtrack->GetOuterParam();
 667
 668         if (trd){Int_t labeltpctrd = track->GetLabel();
 669
 670           AliExternalTrackParam * trdCopy = new AliExternalTrackParam(*trd);
 671           Double_t xktrd=296.0;
 672           Double_t bztrd=5.0;
 673           if(trdCopy->PropagateTo(xktrd,bztrd)){
 674             Float_t xtpc2=trdCopy->GetX();
 675             Float_t ytpc2=trdCopy->GetY();
 676             Float_t ztpc2=trdCopy->GetZ();
 677             Double_t alpha=trdCopy->GetAlpha();
 678             Float_t xpr,ypr ;
 679             xpr=xtpc2*TMath::Cos(alpha)-ytpc2*TMath::Sin(alpha);
 680             ypr=xtpc2*TMath::Sin(alpha)+ytpc2*TMath::Cos(alpha);
 681             AliMCParticle* mcPart = (AliMCParticle*) MCEvent()->GetTrack(abs(labeltpctrd));
 682             Int_t ntreftrd = mcPart->GetNumberOfTrackReferences();
 683
 684             for (Int_t k = 0; k < ntreftrd; k++) {
 685               AliTrackReference* ref = mcPart->GetTrackReference(k);
 686               if(ref->DetectorId() == 3 && ref->Label() == labeltpctrd){
 687
 688                 Float_t xtrd = ref->X();
 689                 Float_t ytrd = ref->Y();
 690                 Float_t ztrd = ref->Z();
 691
 692                 Float_t difx=(xtrd-xpr);
 693                 fHistTPCTRDdeltax->Fill(difx);
 694                 Float_t dify=(ytrd-ypr);
 695                 fHistTPCTRDdeltay->Fill(dify);
 696                 Float_t difz=(ztrd-ztpc2);
 697                 fHistTPCTRDdeltaz->Fill(difz);
 698                 Float_t deltar = TMath::Sqrt(difx * difx + dify * dify + difz * difz);
 699                 fHistTPCTRDdeltar->Fill(deltar);
 700                 Float_t phi_tpc = TMath::ATan2(ypr, xpr);
 701                 Float_t phi_trd = TMath::ATan2(ytrd, xtrd);
 702                 Float_t dphi = phi_trd - phi_tpc;
 703                 if (dphi >   TMath::Pi()) dphi -= 2. * TMath::Pi();
 704                 if (dphi < - TMath::Pi()) dphi += 2. * TMath::Pi();
 705
 706                 fHistTPCTRDdeltaphi->Fill(dphi);
 707                 break;
 708               }
 709
 710             }
 711           }
 712
 713         }//TRD-TPC maching
 714
 715         // CKB
 716         if(pt>2.&&fFieldOn){
 717           TBits bits(esdtrack->GetTPCClusterMap());
 718           for(unsigned int ip = 0;ip < bits.GetNbits();++ip){
 719             if(bits[ip]){
 720               fh2PhiPadRow[indexAC]->Fill(phiDeg,ip);
 721             }
 722           }
 723           for(int i = 0;i < 6;++i){
 724             if(esdtrack->HasPointOnITSLayer(i))fh2PhiLayer[indexAC]->Fill(phiDeg,i);
 725           }
 726         }
 727         else if(!fFieldOn){ // field not on all momenta are set to MPV
 728           TBits bits(esdtrack->GetTPCClusterMap());
 729           for(unsigned int ip = 0;ip < bits.GetNbits();++ip){
 730             if(bits[ip]){
 731               fh2PhiPadRow[indexAC]->Fill(phiDeg,ip);
 732             }
 733             for(int i = 0;i < 6;++i){
 734               if(esdtrack->HasPointOnITSLayer(i))fh2PhiLayer[indexAC]->Fill(phiDeg,i);
 735             }
 736           }
 737
 738         }
 739         // Fill the correlation
 740         if(MCEvent()){
 741           TParticle *part = MCEvent()->Stack()->Particle(TMath::Abs(esdtrack->GetLabel()));
 742           if(part){
 743             Float_t mcEta = part->Eta();
 744             Float_t mcPhi = TMath::RadToDeg()*part->Phi();
 745             Float_t mcPt  = part->Pt();
 746             //             if (pt - mcPt > 20.) {
 747             fh2EtaCorrelation->Fill(mcEta,eta);
 748             fh2EtaCorrelationShift->Fill(eta,mcEta-eta);
 749             fh2PhiCorrelation->Fill(mcPhi,phiDeg);
 750             fh2PhiCorrelationShift->Fill(phiDeg,mcPhi-phiDeg);
 751             fh2PtCorrelation->Fill(mcPt,pt);
 752             fh2PtCorrelationShift->Fill(pt,mcPt-pt);
 753             //}
 754             Double_t sigmaPt = TMath::Sqrt(esdtrack->GetSigma1Pt2());
 755             if (track->Charge() > 0.) {
 756               fPtPullsPos->Fill(mcPt, (1./pt - 1./mcPt) / sigmaPt);
 757               fPtShiftPos->Fill(mcPt, (1./pt - 1./mcPt));
 758             } else {
 759               fPtPullsNeg->Fill(mcPt, (1./pt - 1./mcPt) / sigmaPt);
 760               fPtShiftNeg->Fill(mcPt, (1./pt - 1./mcPt));
 761             }
 762           }
 763           // BKC
 764         }
 765       }
 766
 767   } //track loop
 768
 769   //prevent mem leak for TPConly track
 770   if(fTrackType==2&&tpcP){
 771     delete tpcP;
 772     tpcP = 0;
 773   }
 774
 775   //Printf("%s:%d",(char*)__FILE__,__LINE__);
 776   // CKB
 777   // Vertex eff
 778   if(tpcVertex){
 779     fh1VtxEff->Fill("TPC VTX",1);
 780   }
 781   else{
 782     fh1VtxEff->Fill("NO TPC VTX",1);
 783   }
 784
 785   if(spdVertex){
 786     fh1VtxEff->Fill("SPD VTX",1);
 787   }
 788   else{
 789     fh1VtxEff->Fill("NO SPD VTX",1);
 790   }
 791
 792
 793   // Vertex correlation SPD vs. TPC
 794   if(tpcVertex&&spdVertex){
 795     fh2VtxTpcSpd->Fill(vtxTPC[2],vtxSPD[2]);
 796   }
 797   //  Printf("%s:%d",(char*)__FILE__,__LINE__);
 798   // Multiplicity checks in the SPD
 799   const AliMultiplicity *spdMult = esdE->GetMultiplicity();
 800   // Multiplicity Analysis
 801   Int_t nt = spdMult->GetNumberOfTracklets();
 802   nfc++;
 803   if (nfc == 520) {
 804     nfc = 0;
 805     for (Int_t ib = 1; ib <= 10; ib++) {
 806       Int_t ic = Int_t(fEtaMulti->GetBinContent(ib));
 807       Float_t eta = -1.8 + Float_t(ib-1) * 0.4;
 808       fEtaMultiFluc->Fill(eta, Float_t(ic));
 809       if (ib == 5) fEtaMultiH->Fill(Float_t(ic));
 810     }
 811
 812     for (Int_t ib = 1; ib <= 64; ib++) {
 813       Int_t ic = Int_t(fPhiMulti->GetBinContent(ib));
 814       Float_t phi = 0.05 + Float_t(ib-1) * 0.1;
 815       fPhiMultiFluc->Fill(phi, Float_t(ic));
 816       if (ib  == 2) fPhiMultiH->Fill(Float_t(ic));
 817     }
 818
 819     fEtaMulti->Reset();
 820     fPhiMulti->Reset();
 821   }
 822
 823   for (Int_t j = 0; j < nt; j++) {
 824     fEtaMulti->Fill(spdMult->GetEta(j));
 825     fPhiMulti->Fill(spdMult->GetPhi(j));
 826   }
 827
 828   Short_t chips0 = spdMult->GetNumberOfFiredChips(0);
 829   Short_t chips1 = spdMult->GetNumberOfFiredChips(1);
 830   //Printf("%s:%d",(char*)__FILE__,__LINE__);
 831   Int_t inputCount = 0;
 832
 833   if(spdVertex){
 834     // get multiplicity from ITS tracklets
 835     for (Int_t i=0; i<spdMult->GetNumberOfTracklets(); ++i){
 836       Float_t deltaPhi = spdMult->GetDeltaPhi(i);
 837       // prevent values to be shifted by 2 Pi()
 838       if (deltaPhi < -TMath::Pi())
 839         deltaPhi += TMath::Pi() * 2;
 840       if (deltaPhi > TMath::Pi())
 841         deltaPhi -= TMath::Pi() * 2;
 842       if (TMath::Abs(deltaPhi) > 1)
 843         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, spdMult->GetTheta(i), spdMult->GetPhi(i), deltaPhi);
 844       //trackled Deltaphi
 845       Int_t label1=spdMult->GetLabel(i,0);
 846       Int_t label2=spdMult->GetLabel(i,1);
 847       if ( label1==label2){
 848         if(stack->IsPhysicalPrimary(label1) == kTRUE)
 849           fHistDeltaphiprimaries->Fill(deltaPhi);
 850         else
 851           fHistDeltaphisecondaries->Fill(deltaPhi);
 852       }
 853       else{
 854         fHistDeltaphireject->Fill(deltaPhi);
 855       }
 856       ++inputCount;
 857     }
 858     // Printf("%s:%d",(char*)__FILE__,__LINE__);
 859     fh2MultSpdChips[0]->Fill(chips0,inputCount);
 860     fh2MultSpdChips[1]->Fill(chips1,inputCount);
 861     //Printf("%s:%d",(char*)__FILE__,__LINE__);
 862   }
 863   // BKC
 864
 865   // MC analysis
 866   Int_t nmc = MCEvent()->GetNumberOfTracks();
 867   AliGenEventHeader*       header = MCEvent()->GenEventHeader();
 868   // some test
 869   AliGenCocktailEventHeader* cH = dynamic_cast<AliGenCocktailEventHeader*> (header);
 870   AliGenPythiaEventHeader* pH;
 871   if (cH == 0)
 872     {
 873       header->Dump();
 874       pH = dynamic_cast<AliGenPythiaEventHeader*>(header);
 875     } else {
 876     TObject* entry = cH->GetHeaders()->FindObject("Pythia");
 877     pH = dynamic_cast<AliGenPythiaEventHeader*>(entry);
 878   }
 879   Int_t iproc = -1;
 880
 881   if (pH) iproc = pH->ProcessType();
 882
 883   TArrayF mcV;
 884   header->PrimaryVertex(mcV);
 885   Float_t vzMC = mcV[2];
 886
 887   // CKB
 888   // Fill the correlation with MC
 889   if(spdVertex&&MCEvent()){
 890     Float_t diff = vzMC-vtxSPD[2];
 891     fh2VertexCorrelation[kSPD]->Fill(vzMC,vtxSPD[2]);
 892     fh2VertexCorrelationShift[kSPD]->Fill(vzMC,diff);
 893     fh1VertexShift[kSPD]->Fill(diff);
 894     if(spdVertex->GetZRes()>0)fh1VertexShiftNorm[kSPD]->Fill(diff/spdVertex->GetZRes());
 895   }
 896   if(tpcVertex&&MCEvent()){
 897     Float_t diff = vzMC-vtxTPC[2];
 898     fh2VertexCorrelation[kTPC]->Fill(vzMC,vtxTPC[2]);
 899     fh2VertexCorrelationShift[kTPC]->Fill(vzMC,diff);
 900     fh1VertexShift[kTPC]->Fill(diff);
 901     if(tpcVertex->GetZRes()>0)fh1VertexShiftNorm[kTPC]->Fill(diff/tpcVertex->GetZRes());
 902   }
 903   // BKC
 904
 905
 906   // MC event loop
 907   //printf("MC particle loop \n");
 908   for (Int_t i = 0; i < nmc; i++) {
 909     AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(i);
 910     //printf("MC particle loop %5d \n", i);
 911     // Primaries only
 912     if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
 913     if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
 914     // Int_t charge=Int_t(mcP->Particle()->GetPDG()->Charge());
 915     Float_t phi = mcP->Phi();
 916     Float_t eta = mcP->Eta();
 917     Float_t pt  = mcP->Pt();
 918     if (TMath::Abs(eta) > 0.9) continue;
 919     Int_t ntr = mcP->GetNumberOfTrackReferences();
 920     Int_t nITS = 0;
 921     Int_t nTPC = 0;
 922     Int_t nFRA = 0;
 923     Float_t  x = 0.;
 924     Float_t  y = 0.;
 925
 926     for (Int_t j = 0; j < ntr; j++) {
 927
 928       AliTrackReference* ref = mcP->GetTrackReference(j);
 929
 930       if (ref->DetectorId() == 0) nITS++;
 931       if (ref->DetectorId() == 1) nTPC++;
 932       if (ref->DetectorId() == 2) nFRA++;
 933       if (nTPC == 1) {
 934         x = ref->X();
 935         y = ref->Y();
 936         break;
 937       }
 938     }
 939
 940     fHistMCpt ->Fill(pt);
 941     fHistMCeta->Fill(eta);
 942     fHistMCphi->Fill(phi);
 943     Int_t index=ConvertePDG(mcP->PdgCode());
 944     if(index<8){
 945       fHistMCptCharge [index] ->Fill(pt);
 946       fHistMCetaCharge[index]->Fill(eta);
 947       fHistMCphiCharge[index]->Fill(phi);
 948     }
 949
 950
 951     if (pt > 2.) {
 952       fHistMCHPTeta->Fill(eta);
 953       fHistMCHPTphi->Fill(phi);
 954     }
 955
 956     Bool_t reco = kFALSE;
 957     Int_t multR =  0;
 958     Int_t jold  = -1;
 959
 960     for (Int_t j = 0; j < esdE->GetNumberOfTracks(); j++) {
 961       if (i == labels.At(j)) {
 962         reco = kTRUE;
 963         multR++;
 964         //AliESDtrack* test = esdE->GetTrack(j);
 965         if (multR > 1) {
 966           Int_t nclus = 0;
 967           AliESDtrack* track = esdE->GetTrack(jold);
 968           nclus = track->GetTPCclusters(0);
 969           fHistNCluster->Fill(nclus);
 970
 971
 972           track = esdE->GetTrack(j);
 973           nclus = track->GetTPCclusters(0);
 974           fHistNCluster->Fill(nclus);
 975
 976
 977         }
 978         jold = j;
 979       }
 980     }
 981
 982     fHistRecMult->Fill(Float_t(multR), 1.);
 983     if (multR == 0) {
 984       fHistMCNRpt ->Fill(pt);
 985       fHistMCNReta->Fill(eta);
 986       fHistMCNRphi->Fill(phi);
 987       if(index<8){
 988         fHistMCNRptCharge [index] ->Fill(pt);
 989         fHistMCNRetaCharge[index]->Fill(eta);
 990         fHistMCNRphiCharge[index]->Fill(phi);
 991       }
 992
 993
 994       if (pt > 2.) {
 995         fHistMCNRHPTeta->Fill(eta);
 996         fHistMCNRHPTphi->Fill(phi);
 997       }
 998     }
 999     else if (multR > 1)
1000       {
1001         fHistMULTpt ->Fill(pt);
1002         fHistMULTeta->Fill(eta);
1003         fHistMULTphi->Fill(phi);
1004       }
1005   } // MC particle loop
1006   //printf("End of MC particle loop \n");
1007
1008
1009 }
1010
1011 //________________________________________________________________________
1012 void AliAnalysisTaskEfficiency::Terminate(Option_t *)
1013 {
1014
1015 }
1016
1017
1018 Bool_t AliAnalysisTaskEfficiency::SelectJouri(AliESDtrack* track)
1019 {
1020
1021   Bool_t selected = kTRUE;
1022   AliESDEvent* esdE = (AliESDEvent*) fInputEvent;
1023   // > 50 Clusters
1024   if (track->GetTPCclusters(0) < 50) selected = kFALSE;
1025
1026   const AliESDVertex *vtx=esdE->GetPrimaryVertexSPD();
1027   if (!vtx->GetStatus()) selected = kFALSE;
1028
1029   Double_t zv=vtx->GetZ();
1030
1031
1032   const AliExternalTrackParam *ct=track->GetTPCInnerParam();
1033   if (!ct)  selected = kFALSE;
1034
1035   Double_t xyz[3];
1036   ct->GetXYZ(xyz);
1037   Float_t rv = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
1038   if (rv > 3.0)                     selected = kFALSE;
1039   if (TMath::Abs(xyz[2] - zv)>2.5) selected = kFALSE;
1040   return (selected);
1041
1042 }
1043
1044 Int_t AliAnalysisTaskEfficiency::GetIndexAC(AliESDtrack *track){
1045   if(!track)return -1;
1046
1047   // Crude selection for A and C side
1048   // just with eta
1049   if (track->Eta() > 0) {
1050     if (track->Charge() > 0)
1051       return kPosA;
1052     else
1053       return kNegA;
1054   }
1055   else { // C side
1056     if (track->Charge() > 0)
1057       return kPosC;
1058     else
1059       return kNegC;
1060   }
1061
1062   return -1;
1063
1064 }
1065
1066
1067 Int_t AliAnalysisTaskEfficiency::ConvertePDG(Int_t pdg)
1068 {
1069
1070   // function converts the pdg values for nuclei (d, t, 3He,
1071   // 4He and their antiparticles) into indices for histo-array
1072   // filled in UserExec
1073
1074   if      ( pdg ==  1000010020 )return 0;
1075   else if ( pdg == -1000010020 )return 1;
1076   else if ( pdg ==  1000010030 )return 2;
1077   else if ( pdg == -1000010030 )return 3;
1078   else if ( pdg ==  1000020030 )return 4;
1079   else if ( pdg == -1000020030 )return 5;
1080   else if ( pdg ==  1000020040 )return 6;
1081   else if ( pdg == -1000020040 )return 7;
1082   else return 9;
1083
1084
1085 }