PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskEfficiencyBF.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TH3D.h"
   6 #include "TCanvas.h"
   7 #include "TParticle.h"
   8 #include "TObjArray.h"
   9
  10 #include "AliAnalysisTask.h"
  11 #include "AliAnalysisManager.h"
  12
  13 #include "AliStack.h"
  14 #include "AliMCEvent.h"
  15 #include "AliESDEvent.h"
  16 #include "AliESDInputHandler.h"
  17 #include "AliESDtrackCuts.h"
  18 #include "AliCentrality.h"
  19
  20 #include "AliAnalysisTaskEfficiencyBF.h"
  21
  22 // ---------------------------------------------------------------------
  23 //
  24 // Task for calculating the efficiency of the Balance Function
  25 // for single particles and pairs
  26 //
  27 // Authors: Panos Christakoglou, Michael Weber
  28 //
  29 // ---------------------------------------------------------------------
  30
  31 ClassImp(AliAnalysisTaskEfficiencyBF)
  32
  33 //________________________________________________________________________
  34 AliAnalysisTaskEfficiencyBF::AliAnalysisTaskEfficiencyBF(const char *name)
  35   : AliAnalysisTaskSE(name), fESD(0), fQAList(0), fOutputList(0),
  36   fHistEventStats(0), fHistCentrality(0), fHistNMult(0),
  37   fHistGeneratedEtaPtPhiPlus(0), fHistFindableEtaPtPhiPlus(0),
  38   fHistReconstructedEtaPtPhiPlus(0), fHistSurvivedEtaPtPhiPlus(0),
  39   fHistGeneratedEtaPtPhiMinus(0), fHistFindableEtaPtPhiMinus(0),
  40   fHistReconstructedEtaPtPhiMinus(0), fHistSurvivedEtaPtPhiMinus(0),
  41   fHistGeneratedEtaPtPlusControl(0), fHistFindableEtaPtPlusControl(0),
  42   fHistReconstructedEtaPtPlusControl(0), fHistSurvivedEtaPtPlusControl(0),
  43   fHistGeneratedEtaPtMinusControl(0), fHistFindableEtaPtMinusControl(0),
  44   fHistReconstructedEtaPtMinusControl(0), fHistSurvivedEtaPtMinusControl(0),
  45   fHistGeneratedEtaPtPlusPlus(0), fHistFindableEtaPtPlusPlus(0),
  46   fHistReconstructedEtaPtPlusPlus(0), fHistSurvivedEtaPtPlusPlus(0),
  47   fHistGeneratedEtaPtMinusMinus(0), fHistFindableEtaPtMinusMinus(0),
  48   fHistReconstructedEtaPtMinusMinus(0), fHistSurvivedEtaPtMinusMinus(0),
  49   fHistGeneratedEtaPtPlusMinus(0), fHistFindableEtaPtPlusMinus(0),
  50   fHistReconstructedEtaPtPlusMinus(0), fHistSurvivedEtaPtPlusMinus(0),
  51   fHistGeneratedPhiEtaPlusPlus(0), fHistFindablePhiEtaPlusPlus(0),
  52   fHistReconstructedPhiEtaPlusPlus(0), fHistSurvivedPhiEtaPlusPlus(0),
  53   fHistGeneratedPhiEtaMinusMinus(0), fHistFindablePhiEtaMinusMinus(0),
  54   fHistReconstructedPhiEtaMinusMinus(0), fHistSurvivedPhiEtaMinusMinus(0),
  55   fHistGeneratedPhiEtaPlusMinus(0), fHistFindablePhiEtaPlusMinus(0),
  56   fHistReconstructedPhiEtaPlusMinus(0), fHistSurvivedPhiEtaPlusMinus(0),
  57   fESDtrackCuts(0), fAnalysisMode(0),
  58   fCentralityEstimator("V0M"), fCentralityPercentileMin(0.0), fCentralityPercentileMax(5.0),
  59   fVxMax(3.0), fVyMax(3.0), fVzMax(10.),
  60   fMinNumberOfTPCClusters(80), fMaxChi2PerTPCCluster(4.0), fMaxDCAxy(3.0), fMaxDCAz(3.0),
  61     fMinPt(0.3), fMaxPt(1.5), fMinEta(-0.8), fMaxEta(0.8),fEtaRangeMin(0.0), fEtaRangeMax(1.6), fPtRangeMin(0.3), fPtRangeMax(1.5), fPhiRangeMin(0.0),fPhiRangeMax(360.),fdPhiRangeMax(180.), fEtaBin(100),fdEtaBin(64),fPtBin(49),fPhiBin(100),fdPhiBin(90)   {
  62   // Define input and output slots here
  63   // Input slot #0 works with a TChain
  64   DefineInput(0, TChain::Class());
  65   // Output slot #0 id reserved by the base class for AOD
  66   // Output slot #1 writes into a TH1 container
  67   DefineOutput(1, TList::Class());
  68   DefineOutput(2, TList::Class());
  69 }
  70
  71 //________________________________________________________________________
  72 void AliAnalysisTaskEfficiencyBF::UserCreateOutputObjects() {
  73   // Create histograms
  74   // Called once
  75
  76   fQAList = new TList();
  77   fQAList->SetName("QAList");
  78   fQAList->SetOwner();
  79
  80   fOutputList = new TList();
  81   fOutputList->SetName("OutputList");
  82   fOutputList->SetOwner();
  83
  84   //Event stats.
  85   TString gCutName[4] = {"Total","Offline trigger",
  86                          "Vertex","Analyzed"};
  87   fHistEventStats = new TH1F("fHistEventStats",
  88                              "Event statistics;;N_{events}",
  89                              4,0.5,4.5);
  90   for(Int_t i = 1; i <= 4; i++)
  91     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
  92   fQAList->Add(fHistEventStats);
  93
  94   //ESD analysis
  95   fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
  96                              20,0.5,20.5);
  97   fQAList->Add(fHistCentrality);
  98
  99   //multiplicity (good MC tracks)
 100   TString histName;
 101   histName = "fHistNMult";
 102   fHistNMult = new TH1F(histName.Data(),
 103                         ";N_{mult.}",
 104                         200,0,20000);
 105   fQAList->Add(fHistNMult);
 106
 107   //eta vs pt for MC positives
 108   fHistGeneratedEtaPtPhiPlus = new TH3D("fHistGeneratedEtaPtPhiPlus",
 109                                      "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
 110                                         fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 111   fOutputList->Add(fHistGeneratedEtaPtPhiPlus);
 112   fHistFindableEtaPtPhiPlus = new TH3D("fHistFindableEtaPtPhiPlus",
 113                                      "Findable positive primaries;#eta;p_{T} (GeV/c);#phi",
 114                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 115   fOutputList->Add(fHistFindableEtaPtPhiPlus);
 116   fHistReconstructedEtaPtPhiPlus = new TH3D("fHistReconstructedEtaPtPhiPlus",
 117                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c);#phi",
 118                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 119   fOutputList->Add(fHistReconstructedEtaPtPhiPlus);
 120   fHistSurvivedEtaPtPhiPlus = new TH3D("fHistSurvivedEtaPtPhiPlus",
 121                                      "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
 122                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 123   fOutputList->Add(fHistSurvivedEtaPtPhiPlus);
 124
 125   //eta vs pt for MC negatives
 126   fHistGeneratedEtaPtPhiMinus = new TH3D("fHistGeneratedEtaPtPhiMinus",
 127                                      "Generated positive primaries;#eta;p_{T} (GeV/c);#phi",
 128                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 129   fOutputList->Add(fHistGeneratedEtaPtPhiMinus);
 130   fHistFindableEtaPtPhiMinus = new TH3D("fHistFindableEtaPtPhiMinus",
 131                                      "Findable positive primaries;#eta;p_{T} (GeV/c);#phi",
 132                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 133   fOutputList->Add(fHistFindableEtaPtPhiMinus);
 134   fHistReconstructedEtaPtPhiMinus = new TH3D("fHistReconstructedEtaPtPhiMinus",
 135                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c);#phi",
 136                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 137   fOutputList->Add(fHistReconstructedEtaPtPhiMinus);
 138   fHistSurvivedEtaPtPhiMinus = new TH3D("fHistSurvivedEtaPtPhiMinus",
 139                                      "Survived positive primaries;#eta;p_{T} (GeV/c);#phi",
 140                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax,fPhiBin,fPhiRangeMin,fPhiRangeMax);
 141   fOutputList->Add(fHistSurvivedEtaPtPhiMinus);
 142
 143   //eta vs pt for MC positives (control)
 144   fHistGeneratedEtaPtPlusControl = new TH2F("fHistGeneratedEtaPtPlusControl",
 145                                      "Generated positive primaries;#eta;p_{T} (GeV/c)",
 146                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 147   fOutputList->Add(fHistGeneratedEtaPtPlusControl);
 148   fHistFindableEtaPtPlusControl = new TH2F("fHistFindableEtaPtPlusControl",
 149                                      "Findable positive primaries;#eta;p_{T} (GeV/c)",
 150                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 151   fOutputList->Add(fHistFindableEtaPtPlusControl);
 152   fHistReconstructedEtaPtPlusControl = new TH2F("fHistReconstructedEtaPtPlusControl",
 153                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
 154                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 155   fOutputList->Add(fHistReconstructedEtaPtPlusControl);
 156   fHistSurvivedEtaPtPlusControl = new TH2F("fHistSurvivedEtaPtPlusControl",
 157                                      "Survived positive primaries;#eta;p_{T} (GeV/c)",
 158                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 159   fOutputList->Add(fHistSurvivedEtaPtPlusControl);
 160
 161   //eta vs pt for MC negatives (control)
 162   fHistGeneratedEtaPtMinusControl = new TH2F("fHistGeneratedEtaPtMinusControl",
 163                                      "Generated positive primaries;#eta;p_{T} (GeV/c)",
 164                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 165   fOutputList->Add(fHistGeneratedEtaPtMinusControl);
 166   fHistFindableEtaPtMinusControl = new TH2F("fHistFindableEtaPtMinusControl",
 167                                      "Findable positive primaries;#eta;p_{T} (GeV/c)",
 168                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 169   fOutputList->Add(fHistFindableEtaPtMinusControl);
 170   fHistReconstructedEtaPtMinusControl = new TH2F("fHistReconstructedEtaPtMinusControl",
 171                                      "Reconstructed positive primaries;#eta;p_{T} (GeV/c)",
 172                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 173   fOutputList->Add(fHistReconstructedEtaPtMinusControl);
 174   fHistSurvivedEtaPtMinusControl = new TH2F("fHistSurvivedEtaPtMinusControl",
 175                                      "Survived positive primaries;#eta;p_{T} (GeV/c)",
 176                                      fEtaBin,fMinEta,fMaxEta,fPtBin,fPtRangeMin,fPtRangeMax);
 177   fOutputList->Add(fHistSurvivedEtaPtMinusControl);
 178
 179   //eta vs pt for MC ++
 180   fHistGeneratedEtaPtPlusPlus = new TH2F("fHistGeneratedEtaPtPlusPlus",
 181                                      "Generated ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 182                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 183   fOutputList->Add(fHistGeneratedEtaPtPlusPlus);
 184   fHistFindableEtaPtPlusPlus = new TH2F("fHistFindableEtaPtPlusPlus",
 185                                      "Findable ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 186                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 187   fOutputList->Add(fHistFindableEtaPtPlusPlus);
 188   fHistReconstructedEtaPtPlusPlus = new TH2F("fHistReconstructedEtaPtPlusPlus",
 189                                      "Reconstructed ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 190                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 191   fOutputList->Add(fHistReconstructedEtaPtPlusPlus);
 192   fHistSurvivedEtaPtPlusPlus = new TH2F("fHistSurvivedEtaPtPlusPlus",
 193                                      "Survived ++ primaries;#Delta#eta;p_{T} (GeV/c)",
 194                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 195   fOutputList->Add(fHistSurvivedEtaPtPlusPlus);
 196
 197   //eta vs pt for MC --
 198   fHistGeneratedEtaPtMinusMinus = new TH2F("fHistGeneratedEtaPtMinusMinus",
 199                                      "Generated -- primaries;#Delta#eta;p_{T} (GeV/c)",
 200                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 201   fOutputList->Add(fHistGeneratedEtaPtMinusMinus);
 202   fHistFindableEtaPtMinusMinus = new TH2F("fHistFindableEtaPtMinusMinus",
 203                                      "Findable -- primaries;#Delta#eta;p_{T} (GeV/c)",
 204                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 205   fOutputList->Add(fHistFindableEtaPtMinusMinus);
 206   fHistReconstructedEtaPtMinusMinus = new TH2F("fHistReconstructedEtaPtMinusMinus",
 207                                      "Reconstructed -- primaries;#Delta#eta;p_{T} (GeV/c)",
 208                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 209   fOutputList->Add(fHistReconstructedEtaPtMinusMinus);
 210   fHistSurvivedEtaPtMinusMinus = new TH2F("fHistSurvivedEtaPtMinusMinus",
 211                                      "Survived -- primaries;#Delta#eta;p_{T} (GeV/c)",
 212                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 213   fOutputList->Add(fHistSurvivedEtaPtMinusMinus);
 214
 215   //eta vs pt for MC +-
 216   fHistGeneratedEtaPtPlusMinus = new TH2F("fHistGeneratedEtaPtPlusMinus",
 217                                      "Generated +- primaries;#Delta#eta;p_{T} (GeV/c)",
 218                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 219   fOutputList->Add(fHistGeneratedEtaPtPlusMinus);
 220   fHistFindableEtaPtPlusMinus = new TH2F("fHistFindableEtaPtPlusMinus",
 221                                      "Findable +- primaries;#Delta#eta;p_{T} (GeV/c)",
 222                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 223   fOutputList->Add(fHistFindableEtaPtPlusMinus);
 224   fHistReconstructedEtaPtPlusMinus = new TH2F("fHistReconstructedEtaPtPlusMinus",
 225                                      "Reconstructed +- primaries;#Delta#eta;p_{T} (GeV/c)",
 226                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 227   fOutputList->Add(fHistReconstructedEtaPtPlusMinus);
 228   fHistSurvivedEtaPtPlusMinus = new TH2F("fHistSurvivedEtaPtPlusMinus",
 229                                      "Survived +- primaries;#Delta#eta;p_{T} (GeV/c)",
 230                                      fdEtaBin,fEtaRangeMin,fEtaRangeMax,fPtBin,fPtRangeMin,fPtRangeMax);
 231   fOutputList->Add(fHistSurvivedEtaPtPlusMinus);
 232
 233   //=============================//
 234    //phi vs eta for MC ++
 235   fHistGeneratedPhiEtaPlusPlus = new TH2F("fHistGeneratedPhiEtaPlusPlus",
 236                                      "Generated ++ primaries;#Delta#phi",
 237                                           fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 238   fOutputList->Add(fHistGeneratedPhiEtaPlusPlus);
 239   fHistFindablePhiEtaPlusPlus = new TH2F("fHistFindablePhiEtaPlusPlus",
 240                                      "Findable ++ primaries;#Delta#phi;#Delta#eta",
 241                                          fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 242   fOutputList->Add(fHistFindablePhiEtaPlusPlus);
 243   fHistReconstructedPhiEtaPlusPlus = new TH2F("fHistReconstructedPhiEtaPlusPlus",
 244                                      "Reconstructed ++ primaries;#Delta#phi;#Delta#eta",
 245                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 246   fOutputList->Add(fHistReconstructedPhiEtaPlusPlus);
 247   fHistSurvivedPhiEtaPlusPlus = new TH2F("fHistSurvivedPhiEtaPlusPlus",
 248                                      "Survived ++ primaries;#Delta#phi;#Delta#eta",
 249                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 250   fOutputList->Add(fHistSurvivedPhiEtaPlusPlus);
 251
 252   //phi vs eta for MC --
 253   fHistGeneratedPhiEtaMinusMinus = new TH2F("fHistGeneratedPhiEtaMinusMinus",
 254                                      "Generated -- primaries;#Delta#phi;#Delta#eta",
 255                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 256   fOutputList->Add(fHistGeneratedPhiEtaMinusMinus);
 257   fHistFindablePhiEtaMinusMinus = new TH2F("fHistFindablePhiEtaMinusMinus",
 258                                      "Findable -- primaries;#Delta#phi;#Delta#eta",
 259                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 260   fOutputList->Add(fHistFindablePhiEtaMinusMinus);
 261   fHistReconstructedPhiEtaMinusMinus = new TH2F("fHistReconstructedPhiEtaMinusMinus",
 262                                      "Reconstructed -- primaries;#Delta#phi;#Delta#eta",
 263                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 264   fOutputList->Add(fHistReconstructedPhiEtaMinusMinus);
 265   fHistSurvivedPhiEtaMinusMinus = new TH2F("fHistSurvivedPhiEtaMinusMinus",
 266                                      "Survived -- primaries;#Delta#phi;#Delta#eta",
 267                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 268   fOutputList->Add(fHistSurvivedPhiEtaMinusMinus);
 269
 270   //phi vs eta for MC +-
 271   fHistGeneratedPhiEtaPlusMinus = new TH2F("fHistGeneratedPhiEtaPlusMinus",
 272                                      "Generated +- primaries;#Delta#phi;#Delta#eta",
 273                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 274   fOutputList->Add(fHistGeneratedPhiEtaPlusMinus);
 275   fHistFindablePhiEtaPlusMinus = new TH2F("fHistFindablePhiEtaPlusMinus",
 276                                      "Findable +- primaries;#Delta#phi;#Delta#eta",
 277                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 278   fOutputList->Add(fHistFindablePhiEtaPlusMinus);
 279   fHistReconstructedPhiEtaPlusMinus = new TH2F("fHistReconstructedPhiEtaPlusMinus",
 280                                      "Reconstructed +- primaries;#Delta#phi;#Delta#eta",
 281                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 282   fOutputList->Add(fHistReconstructedPhiEtaPlusMinus);
 283   fHistSurvivedPhiEtaPlusMinus = new TH2F("fHistSurvivedPhiEtaPlusMinus",
 284                                      "Survived +- primaries;#Delta#phi;#Delta#eta",
 285                                      fdPhiBin,fPhiRangeMin,fdPhiRangeMax,fdEtaBin,fEtaRangeMin,fEtaRangeMax);
 286   fOutputList->Add(fHistSurvivedPhiEtaPlusMinus);
 287   //=============================//
 288
 289   fQAList->Print();
 290   fOutputList->Print();
 291
 292   PostData(1, fQAList);
 293   PostData(2, fOutputList);
 294 }
 295
 296 //________________________________________________________________________
 297 void AliAnalysisTaskEfficiencyBF::UserExec(Option_t *) {
 298   // Main loop
 299   // Called for each event
 300
 301   // Post output data.
 302   //ESD analysis
 303   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 304   if (!fESD) {
 305     printf("ERROR: fESD not available\n");
 306     return;
 307   }
 308
 309   AliMCEvent* mcEvent = MCEvent();
 310   if (!mcEvent) {
 311     Printf("ERROR: Could not retrieve MC event");
 312     return;
 313   }
 314   AliStack* stack = mcEvent->Stack();
 315   if (!stack) {
 316     Printf("ERROR: Could not retrieve MC stack");
 317     return;
 318   }
 319
 320   // arrays for 2 particle histograms
 321   Int_t nMCLabelCounter         = 0;
 322   const Int_t maxMCLabelCounter = 20000;
 323
 324   Double_t eta[maxMCLabelCounter];
 325   Double_t pt[maxMCLabelCounter];
 326   Double_t phi[maxMCLabelCounter];
 327   Int_t level[maxMCLabelCounter];
 328   Int_t charge[maxMCLabelCounter];
 329
 330   //AliInfo(Form("%d %d",mcEvent->GetNumberOfTracks(),fESD->GetNumberOfTracks()));
 331   fHistEventStats->Fill(1); //all events
 332
 333   //Centrality stuff
 334   AliCentrality *centrality = fESD->GetCentrality();
 335   Int_t nCentrality = 0;
 336   nCentrality = (Int_t)(centrality->GetCentralityPercentile(fCentralityEstimator.Data())/10.);
 337
 338   //Printf("Centrality: %lf",centrality->GetCentralityPercentile(fCentralityEstimator.Data()));
 339
 340   if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
 341                                           fCentralityPercentileMax,
 342                                           fCentralityEstimator.Data())) {
 343     fHistEventStats->Fill(2); //triggered + centrality
 344     fHistCentrality->Fill(nCentrality+1);
 345
 346     //Printf("Centrality selection: %lf - %lf",fCentralityPercentileMin,fCentralityPercentileMax);
 347
 348     if(fAnalysisMode.CompareTo("TPC") == 0 ) {
 349       const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
 350       if(vertex) {
 351         if(vertex->GetNContributors() > 0) {
 352           if(vertex->GetZRes() != 0) {
 353             fHistEventStats->Fill(3); //events with a proper vertex
 354             if(TMath::Abs(vertex->GetXv()) < fVxMax) {
 355               if(TMath::Abs(vertex->GetYv()) < fVyMax) {
 356                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
 357                   fHistEventStats->Fill(4); //analyzed events
 358
 359                   Int_t nMCParticles = mcEvent->GetNumberOfTracks();
 360                   TArrayI labelMCArray(nMCParticles);
 361
 362                   for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
 363                     AliMCParticle *mcTrack = (AliMCParticle*) mcEvent->GetTrack(iTracks);
 364                     if (!mcTrack) {
 365                       Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
 366                       continue;
 367                     }
 368
 369                     //exclude particles generated out of the acceptance
 370                     Double_t vz = mcTrack->Zv();
 371                     if (TMath::Abs(vz) > 50.) continue;
 372
 373                     //acceptance
 374                     if(TMath::Abs(mcTrack->Eta()) > fMaxEta)
 375                       continue;
 376                     if((mcTrack->Pt() > fMaxPt)||(mcTrack->Pt() < fMinPt))
 377                       continue;
 378                     if((mcTrack->Phi() > fPhiRangeMax)||(mcTrack->Phi() < fPhiRangeMin))
 379                       continue;
 380
 381                     TParticle* particle = mcTrack->Particle();
 382                     if(!particle) continue;
 383                     if(!stack->IsPhysicalPrimary(iTracks)) continue;
 384
 385                     if(iTracks <= stack->GetNprimary()) {
 386                       Short_t gMCCharge = mcTrack->Charge();
 387                       Double_t phiRad = particle->Phi();
 388                       Double_t phiDeg = phiRad*TMath::RadToDeg();
 389
 390                       if(gMCCharge > 0)
 391                         fHistGeneratedEtaPtPhiPlus->Fill(particle->Eta(),
 392                                                          particle->Pt(),
 393                                                          phiDeg);
 394                       else if(gMCCharge < 0)
 395                         fHistGeneratedEtaPtPhiMinus->Fill(particle->Eta(),
 396                                                           particle->Pt(),
 397                                                           phiDeg);
 398
 399
 400                       // findable tracks --> DOES NOT WORK????
 401                       // Loop over Track References
 402                       Bool_t labelTPC = kTRUE;
 403                       AliTrackReference* trackRef = 0;
 404
 405                       for (Int_t iTrackRef = 0; iTrackRef < mcTrack->GetNumberOfTrackReferences(); iTrackRef++) {
 406                         trackRef = mcTrack->GetTrackReference(iTrackRef);
 407                         if(trackRef) {
 408                           Int_t detectorId = trackRef->DetectorId();
 409                           if (detectorId == AliTrackReference::kTPC) {
 410                             labelTPC = kTRUE;
 411                             break;
 412                           }
 413                         }
 414                       }//loop over track references
 415
 416                       if(labelTPC) {
 417                         labelMCArray.AddAt(iTracks,nMCLabelCounter);
 418
 419                         if(nMCLabelCounter >= maxMCLabelCounter){
 420                           AliWarning(Form("MC Label Counter > Limit (%d) --> stop loop here",maxMCLabelCounter));
 421                           break;
 422                         }
 423
 424                         //fill the arrays for 2 particle analysis
 425                         eta[nMCLabelCounter]    = particle->Eta();
 426                         pt[nMCLabelCounter]     = particle->Pt();
 427                         phi[nMCLabelCounter]     = particle->Phi()*TMath::RadToDeg();
 428                         charge[nMCLabelCounter] = gMCCharge;
 429
 430                         // findable = generated in this case!
 431
 432                         level[nMCLabelCounter]  = 1;
 433                         nMCLabelCounter += 1;
 434
 435
 436                         if(gMCCharge > 0)
 437                           fHistFindableEtaPtPhiPlus->Fill(particle->Eta(),
 438                                                           particle->Pt(),
 439                                                           phiDeg);
 440                         else if(gMCCharge < 0)
 441                           fHistFindableEtaPtPhiMinus->Fill(particle->Eta(),
 442                                                            particle->Pt(),
 443                                                            phiDeg);
 444                       }
 445                     }//primaries
 446                   }//loop over MC particles
 447
 448                   fHistNMult->Fill(nMCLabelCounter);
 449
 450                   // not used so far
 451                   //Float_t dcaXY = 0.0, dcaZ = 0.0;
 452
 453                   //ESD track loop
 454                   Int_t nGoodTracks = fESD->GetNumberOfTracks();
 455
 456                   TArrayI labelArray(nGoodTracks);
 457                   Int_t labelCounter = 0;
 458                   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 459                     AliESDtrack* track = fESD->GetTrack(iTracks);
 460                     //AliESDtrack* track = fESDtrackCuts->GetTPCOnlyTrack(fESD,iTracks);
 461                     if(!track) continue;
 462
 463                     AliESDtrack *tpcOnlyTrack = new AliESDtrack();
 464
 465                     if (!track->FillTPCOnlyTrack(*tpcOnlyTrack)) {
 466                       delete tpcOnlyTrack;
 467                       continue;
 468                     }
 469
 470                     Int_t label = TMath::Abs(track->GetTPCLabel());
 471                     if(IsLabelUsed(labelArray,label)) continue;
 472                     labelArray.AddAt(label,labelCounter);
 473                     labelCounter += 1;
 474
 475                     Bool_t iFound = kFALSE;
 476                     Int_t mcGoods = nMCLabelCounter;
 477                     for (Int_t k = 0; k < mcGoods; k++) {
 478                       Int_t mcLabel = labelMCArray.At(k);
 479                       iFound = kFALSE;
 480
 481                       if (mcLabel != TMath::Abs(label)) continue;
 482                       if(mcLabel != label) continue;
 483                       if(label > stack->GetNtrack()) continue;
 484
 485                       TParticle *particle = stack->Particle(label);
 486                       if(!particle) continue;
 487
 488                       //acceptance
 489                       if(TMath::Abs(particle->Eta()) > fMaxEta)
 490                         continue;
 491                       if((particle->Pt() > fMaxPt)||(particle->Pt() <  fMinPt))
 492                         continue;
 493                       if((particle->Phi() > fPhiRangeMax)||(particle->Phi() < fPhiRangeMin))
 494                         continue;
 495
 496                       if(!stack->IsPhysicalPrimary(label)) continue;
 497
 498                       if(label <= stack->GetNprimary()) {
 499
 500                         // reconstructed
 501                         level[k]  = 2;
 502
 503                         Short_t gCharge = track->Charge();
 504                         Double_t phiRad = particle->Phi();
 505                         Double_t phiDeg = phiRad*TMath::RadToDeg();
 506
 507                         if(gCharge > 0)
 508                           fHistReconstructedEtaPtPhiPlus->Fill(particle->Eta(),
 509                                                                particle->Pt(),
 510                                                                phiDeg);
 511                         else if(gCharge < 0)
 512                           fHistReconstructedEtaPtPhiMinus->Fill(particle->Eta(),
 513                                                                 particle->Pt(),
 514                                                                 phiDeg);
 515
 516                         // track cuts + analysis kinematic cuts
 517                         if(fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta()) < fMaxEta && track->Pt() > fMinPt && track->Pt() < fMaxPt ){
 518
 519                           // survived
 520                           level[k]  = 3;
 521
 522                           if(gCharge > 0)
 523                             fHistSurvivedEtaPtPhiPlus->Fill(particle->Eta(),
 524                                                             particle->Pt(),
 525                                                             phiDeg);
 526                           else if(gCharge < 0)
 527                             fHistSurvivedEtaPtPhiMinus->Fill(particle->Eta(),
 528                                                              particle->Pt(),
 529                                                              phiDeg);
 530
 531                         }//track cuts
 532                       }//primary particles
 533                     }//findable track loop
 534                   }//ESD track loop
 535
 536                   labelMCArray.Reset();
 537                   labelArray.Reset();
 538
 539
 540                 }//Vz cut
 541               }//Vy cut
 542             }//Vx cut
 543           }//Vz resolution
 544         }//number of contributors
 545       }//valid vertex
 546     }//TPC analysis mode
 547   }//centrality
 548
 549
 550
 551   // Here comes the 2 particle analysis
 552   // loop over all good MC particles
 553   for (Int_t i = 0; i < nMCLabelCounter ; i++) {
 554
 555     // control 1D histograms (charge might be different?)
 556     if(charge[i] > 0){
 557       if(level[i] > 0) fHistGeneratedEtaPtPlusControl->Fill(eta[i],pt[i]);
 558       if(level[i] > 1) fHistReconstructedEtaPtPlusControl->Fill(eta[i],pt[i]);
 559       if(level[i] > 2) fHistSurvivedEtaPtPlusControl->Fill(eta[i],pt[i]);
 560     }
 561     else if(charge[i] < 0){
 562       if(level[i] > 0) fHistGeneratedEtaPtMinusControl->Fill(eta[i],pt[i]);
 563       if(level[i] > 1) fHistReconstructedEtaPtMinusControl->Fill(eta[i],pt[i]);
 564       if(level[i] > 2) fHistSurvivedEtaPtMinusControl->Fill(eta[i],pt[i]);
 565     }
 566
 567
 568     for (Int_t j = i+1; j < nMCLabelCounter ; j++) {
 569
 570       if(charge[i] > 0 && charge[j] > 0 ){
 571         if(level[i] > 0 && level[j] > 0) {
 572           fHistGeneratedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 573           if (TMath::Abs(phi[i]-phi[j]) < 180)
 574           fHistGeneratedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 575         }
 576         if(level[i] > 1 && level[j] > 1) {
 577           fHistReconstructedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 578           if (TMath::Abs(phi[i]-phi[j]) < 180)
 579           fHistReconstructedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 580         }
 581         if(level[i] > 2 && level[j] > 2) {
 582           fHistSurvivedEtaPtPlusPlus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 583           if (TMath::Abs(phi[i]-phi[j]) < 180)
 584           fHistSurvivedPhiEtaPlusPlus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 585         }
 586       }
 587
 588       else if(charge[i] < 0 && charge[j] < 0 ){
 589         if(level[i] > 0 && level[j] > 0) {
 590           fHistGeneratedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 591           if (TMath::Abs(phi[i]-phi[j]) < 180)
 592           fHistGeneratedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 593         }
 594         if(level[i] > 1 && level[j] > 1) {
 595           fHistReconstructedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 596           fHistReconstructedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 597         }
 598         if(level[i] > 2 && level[j] > 2) {
 599           fHistSurvivedEtaPtMinusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 600           if (TMath::Abs(phi[i]-phi[j]) < 180)
 601           fHistSurvivedPhiEtaMinusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 602         }
 603       }
 604
 605       else if((charge[i] > 0 && charge[j] < 0)||(charge[i] < 0 && charge[j] > 0)){
 606         if(level[i] > 0 && level[j] > 0) {
 607           fHistGeneratedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 608           if (TMath::Abs(phi[i]-phi[j]) < 180)
 609           fHistGeneratedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 610         }
 611         if(level[i] > 1 && level[j] > 1) {
 612           fHistReconstructedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 613           fHistReconstructedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 614         }
 615         if(level[i] > 2 && level[j] > 2) {
 616           fHistSurvivedEtaPtPlusMinus->Fill(TMath::Abs(eta[i]-eta[j]),pt[i]);
 617           if (TMath::Abs(phi[i]-phi[j]) < 180)
 618           fHistSurvivedPhiEtaPlusMinus->Fill(TMath::Abs(phi[i]-phi[j]),TMath::Abs(eta[i]-eta[j]));
 619         }
 620       }
 621     }
 622   }
 623
 624   //Checking Entries for generated and survived
 625   /*TH1F* GeneratedEtaPt = (TH1F*)fHistGeneratedEtaPtPlusMinus->ProjectionX("GeneratedEtaPt",5,15);
 626     Int_t xGeneratedPt = fHistGeneratedEtaPtPlusMinus->GetNbinsX();
 627     for (Int_t h=1; h < xGeneratedPt+1; h++){
 628     Double_t binEntriesGenerated = GeneratedEtaPt->GetBinContent(h);
 629     Printf("binEntriesGenerated: %lf - xGeneratedPt: %d",binEntriesGenerated,h);
 630     }
 631     TH1F* GeneratedPhiEta = (TH1F*)fHistGeneratedPhiEtaPlusMinus->ProjectionY("GeneratedPhiEta",5,15);
 632     Int_t yGeneratedPhi = fHistGeneratedPhiEtaPlusMinus->GetNbinsY();
 633     for (Int_t h=1; h < yGeneratedPhi+1; h++){
 634     Double_t binEntriesGenerated = GeneratedPhiEta->GetBinContent(h);
 635     Printf("binEntriesGenerated: %lf - yGeneratedPhi: %d",binEntriesGenerated,h);
 636     }*/
 637
 638   /*TH1F* SurvivedEtaPt = (TH1F*)fHistSurvivedEtaPtPlusMinus->ProjectionX("SurvivedEtaPt",5,15);
 639     Int_t xSurvivedPt = fHistSurvivedEtaPtPlusMinus->GetNbinsX();
 640     for (Int_t h=1; h < xSurvivedPt+1; h++){
 641     Double_t binEntriesSurvived = SurvivedEtaPt->GetBinContent(h);
 642     Printf("binEntriesSurvived: %lf - xSurvivedPt: %d",binEntriesSurvived,h);
 643     }
 644     TH1F* SurvivedPhiEta = (TH1F*)fHistSurvivedPhiEtaPlusMinus->ProjectionY("SurvivedPhiEta",5,15);
 645     Int_t ySurvivedPhi = fHistSurvivedPhiEtaPlusMinus->GetNbinsY();
 646     for (Int_t h=1; h < ySurvivedPhi+1; h++){
 647     Double_t binEntriesSurvived = SurvivedPhiEta->GetBinContent(h);
 648     Printf("binEntriesSurvived: %lf - ySurvivedPhi: %d",binEntriesSurvived,h);
 649     }*/
 650 }
 651   //________________________________________________________________________
 652 void AliAnalysisTaskEfficiencyBF::Terminate(Option_t *) {
 653   // Draw result to the screen
 654   // Called once at the end of the query
 655
 656   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 657   if (!fOutputList) {
 658     printf("ERROR: Output list not available\n");
 659     return;
 660   }
 661 }
 662
 663 //____________________________________________________________________//
 664 Bool_t AliAnalysisTaskEfficiencyBF::IsLabelUsed(TArrayI labelArray, Int_t label) {
 665   //Checks if the label is used already
 666   Bool_t status = kFALSE;
 667   for(Int_t i = 0; i < labelArray.GetSize(); i++) {
 668     if(labelArray.At(i) == label)
 669       status = kTRUE;
 670   }
 671
 672   return status;
 673 }