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