PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.cxx

   1 #include <vector>\r
   2 #include "TChain.h"\r
   3 #include "TList.h"\r
   4 #include "TCanvas.h"\r
   5 #include "TLorentzVector.h"\r
   6 #include "TGraphErrors.h"\r
   7 #include "TH1F.h"\r
   8 #include "TH2F.h"\r
   9 #include "TArrayF.h"\r
  10 #include "TF1.h"\r
  11 #include "TRandom.h"\r
  12 \r
  13 #include "AliLog.h"\r
  14 \r
  15 #include "AliAnalysisTaskSE.h"\r
  16 #include "AliAnalysisManager.h"\r
  17 \r
  18 #include "AliESDVertex.h"\r
  19 #include "AliESDEvent.h"\r
  20 #include "AliESDInputHandler.h"\r
  21 #include "AliAODEvent.h"\r
  22 #include "AliAODTrack.h"\r
  23 #include "AliAODInputHandler.h"\r
  24 #include "AliGenEventHeader.h"\r
  25 #include "AliGenHijingEventHeader.h"\r
  26 #include "AliMCEventHandler.h"\r
  27 #include "AliMCEvent.h"\r
  28 #include "AliMixInputEventHandler.h"\r
  29 #include "AliStack.h"\r
  30 \r
  31 #include "TH2D.h"    \r
  32 #include "AliTHn.h"             \r
  33 \r
  34 #include "AliEventPoolManager.h" \r
  35 \r
  36 #include "AliAnalysisTaskTriggeredBF.h"\r
  37 #include "AliBalanceTriggered.h"\r
  38 \r
  39 \r
  40 // Analysis task for the TriggeredBF code\r
  41 // Authors: Panos.Christakoglou@nikhef.nl, m.weber@cern.ch\r
  42 \r
  43 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+\r
  44 //\r
  45 // For the V0 part:\r
  46 // --> AliAnalysisTaskExtractV0AOD (by david.chinellato@gmail.com)\r
  47 //\r
  48 // +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+\r
  49 \r
  50 using std::cout;\r
  51 using std::endl;\r
  52 using std::vector;\r
  53 \r
  54 ClassImp(AliAnalysisTaskTriggeredBF)\r
  55 \r
  56 //________________________________________________________________________\r
  57 AliAnalysisTaskTriggeredBF::AliAnalysisTaskTriggeredBF(const char *name) \r
  58 : AliAnalysisTaskSE(name), \r
  59   fBalance(0),\r
  60   fRunShuffling(kFALSE),\r
  61   fShuffledBalance(0),\r
  62   fRunMixing(kFALSE),\r
  63   fMixingTracks(50000),\r
  64   fMixedBalance(0),\r
  65   fPoolMgr(0),\r
  66   fRunV0(kFALSE),\r
  67   fPIDResponse(0),\r
  68   fPIDCombined(0),\r
  69   fList(0),\r
  70   fListTriggeredBF(0),\r
  71   fListTriggeredBFS(0),\r
  72   fListTriggeredBFM(0),\r
  73   fHistListPIDQA(0),\r
  74   fHistListV0(0),\r
  75   fHistEventStats(0),\r
  76   fHistCentStats(0),\r
  77   fHistTriggerStats(0),\r
  78   fHistTrackStats(0),\r
  79   fHistVx(0),\r
  80   fHistVy(0),\r
  81   fHistVz(0),\r
  82   fHistClus(0),\r
  83   fHistDCA(0),\r
  84   fHistChi2(0),\r
  85   fHistPt(0),\r
  86   fHistEta(0),\r
  87   fHistPhi(0),\r
  88   fHistPhiBefore(0),\r
  89   fHistPhiAfter(0),\r
  90   fHistV0M(0),\r
  91   fHistRefTracks(0),\r
  92   fHistV0MultiplicityBeforeTrigSel(0),\r
  93   fHistV0MultiplicityForTrigEvt(0),\r
  94   fHistV0MultiplicityForSelEvt(0),\r
  95   fHistV0MultiplicityForSelEvtNoTPCOnly(0),\r
  96   fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup(0),\r
  97   fHistMultiplicityBeforeTrigSel(0),\r
  98   fHistMultiplicityForTrigEvt(0),\r
  99   fHistMultiplicity(0),\r
 100   fHistMultiplicityNoTPCOnly(0),\r
 101   fHistMultiplicityNoTPCOnlyNoPileup(0),\r
 102   fHistV0InvMassK0(0),\r
 103   fHistV0InvMassLambda(0),\r
 104   fHistV0InvMassAntiLambda(0),\r
 105   fHistV0Armenteros(0),\r
 106   fHistV0SelInvMassK0(0),\r
 107   fHistV0SelInvMassLambda(0),\r
 108   fHistV0SelInvMassAntiLambda(0),\r
 109   fHistV0SelArmenteros(0),\r
 110   fCentralityEstimator("V0M"),\r
 111   fUseCentrality(kFALSE),\r
 112   fCentralityPercentileMin(0.), \r
 113   fCentralityPercentileMax(5.),\r
 114   fImpactParameterMin(0.),\r
 115   fImpactParameterMax(20.),\r
 116   fUseMultiplicity(kFALSE),\r
 117   fNumberOfAcceptedTracksMin(0),\r
 118   fNumberOfAcceptedTracksMax(10000),\r
 119   fHistNumberOfAcceptedTracks(0),\r
 120   fUseOfflineTrigger(kFALSE),\r
 121   fVxMax(0.3),\r
 122   fVyMax(0.3),\r
 123   fVzMax(10.),\r
 124   nAODtrackCutBit(128),\r
 125   fPtMin(0.3),\r
 126   fPtMax(1.5),\r
 127   fEtaMin(-0.8),\r
 128   fEtaMax(-0.8),\r
 129   fDCAxyCut(-1),\r
 130   fDCAzCut(-1),\r
 131   fTPCchi2Cut(-1),\r
 132   fNClustersTPCCut(-1)\r
 133 {\r
 134   // Constructor\r
 135   // Define input and output slots here\r
 136   // Input slot #0 works with a TChain\r
 137   DefineInput(0, TChain::Class());\r
 138   // Output slot #0 writes into a TH1 container\r
 139   DefineOutput(1, TList::Class());\r
 140   DefineOutput(2, TList::Class());\r
 141   DefineOutput(3, TList::Class());\r
 142   DefineOutput(4, TList::Class());\r
 143   DefineOutput(5, TList::Class());\r
 144 }\r
 145 \r
 146 //________________________________________________________________________\r
 147 AliAnalysisTaskTriggeredBF::~AliAnalysisTaskTriggeredBF() {\r
 148 \r
 149   // Destructor\r
 150 \r
 151 }\r
 152 \r
 153 //________________________________________________________________________\r
 154 void AliAnalysisTaskTriggeredBF::UserCreateOutputObjects() {\r
 155   // Create histograms\r
 156   // Called once\r
 157 \r
 158   // global switch disabling the reference \r
 159   // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
 160   Bool_t oldStatus = TH1::AddDirectoryStatus();\r
 161   TH1::AddDirectory(kFALSE);\r
 162 \r
 163   if(!fBalance) {\r
 164     fBalance = new AliBalanceTriggered();\r
 165     fBalance->SetAnalysisLevel("AOD");\r
 166   }\r
 167   if(fRunShuffling) {\r
 168     if(!fShuffledBalance) {\r
 169       fShuffledBalance = new AliBalanceTriggered();\r
 170       fShuffledBalance->SetAnalysisLevel("AOD");\r
 171     }\r
 172   }\r
 173   if(fRunMixing) {\r
 174     if(!fMixedBalance) {\r
 175       fMixedBalance = new AliBalanceTriggered();\r
 176       fMixedBalance->SetAnalysisLevel("AOD");\r
 177     }\r
 178   }\r
 179 \r
 180   //QA list\r
 181   fList = new TList();\r
 182   fList->SetName("listQA");\r
 183   fList->SetOwner();\r
 184 \r
 185   //Balance Function list\r
 186   fListTriggeredBF = new TList();\r
 187   fListTriggeredBF->SetName("listTriggeredBF");\r
 188   fListTriggeredBF->SetOwner();\r
 189 \r
 190   if(fRunShuffling) {\r
 191     fListTriggeredBFS = new TList();\r
 192     fListTriggeredBFS->SetName("listTriggeredBFShuffled");\r
 193     fListTriggeredBFS->SetOwner();\r
 194   }\r
 195   if(fRunMixing) {\r
 196     fListTriggeredBFM = new TList();\r
 197     fListTriggeredBFM->SetName("listTriggeredBFMixed");\r
 198     fListTriggeredBFM->SetOwner();\r
 199   }\r
 200   \r
 201   \r
 202   //Event stats.\r
 203   TString gCutName[4] = {"Total","Offline trigger",\r
 204                          "Vertex","Analyzed"};\r
 205   fHistEventStats = new TH1F("fHistEventStats",\r
 206                              "Event statistics;;N_{events}",\r
 207                              4,0.5,4.5);\r
 208   for(Int_t i = 1; i <= 4; i++)\r
 209     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 210   fList->Add(fHistEventStats);\r
 211   \r
 212   TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
 213   fHistCentStats = new TH2F("fHistCentStats",\r
 214                             "Centrality statistics;;Cent percentile",\r
 215                             9,-0.5,8.5,220,-5,105);\r
 216   for(Int_t i = 1; i <= 9; i++)\r
 217     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
 218   fList->Add(fHistCentStats);\r
 219   \r
 220   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
 221   fList->Add(fHistTriggerStats);\r
 222   \r
 223   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
 224   fList->Add(fHistTrackStats);\r
 225 \r
 226   fHistNumberOfAcceptedTracks = new TH1F("fHistNumberOfAcceptedTracks",";N_{acc.};Entries",4001,-0.5,4000.5);\r
 227   fList->Add(fHistNumberOfAcceptedTracks);\r
 228 \r
 229   // Vertex distributions\r
 230   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
 231   fList->Add(fHistVx);\r
 232   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
 233   fList->Add(fHistVy);\r
 234   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
 235   fList->Add(fHistVz);\r
 236 \r
 237   // QA histograms\r
 238   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
 239   fList->Add(fHistClus);\r
 240   fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
 241   fList->Add(fHistChi2);\r
 242   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
 243   fList->Add(fHistDCA);\r
 244   fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
 245   fList->Add(fHistPt);\r
 246   fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
 247   fList->Add(fHistEta);\r
 248   fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
 249   fList->Add(fHistPhi);\r
 250   fHistPhiBefore  = new TH1F("fHistPhiBefore","#phi distribution",200,0.,2*TMath::Pi());\r
 251   fList->Add(fHistPhiBefore);\r
 252   fHistPhiAfter  = new TH1F("fHistPhiAfter","#phi distribution",200,0.,2*TMath::Pi());\r
 253   fList->Add(fHistPhiAfter);\r
 254   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
 255   fList->Add(fHistV0M);\r
 256   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
 257   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
 258   for(Int_t i = 1; i <= 6; i++)\r
 259     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
 260   fList->Add(fHistRefTracks);\r
 261 \r
 262   //------------------------------------------------\r
 263   // V0 Multiplicity Histograms\r
 264   //------------------------------------------------\r
 265   if(fRunV0){\r
 266     fHistListV0 = new TList();\r
 267     fHistListV0->SetOwner();  // See http://root.cern.ch/root/html/TCollection.html#TCollection:SetOwner\r
 268     \r
 269     if(! fHistV0MultiplicityBeforeTrigSel) {\r
 270       fHistV0MultiplicityBeforeTrigSel = new TH1F("fHistV0MultiplicityBeforeTrigSel", \r
 271                                                   "V0s per event (before Trig. Sel.);Nbr of V0s/Evt;Events", \r
 272                                                   25, 0, 25);\r
 273       fHistListV0->Add(fHistV0MultiplicityBeforeTrigSel);\r
 274     }\r
 275     \r
 276     if(! fHistV0MultiplicityForTrigEvt) {\r
 277       fHistV0MultiplicityForTrigEvt = new TH1F("fHistV0MultiplicityForTrigEvt", \r
 278                                                "V0s per event (for triggered evt);Nbr of V0s/Evt;Events", \r
 279                                                25, 0, 25);\r
 280       fHistListV0->Add(fHistV0MultiplicityForTrigEvt);\r
 281     }\r
 282     \r
 283     if(! fHistV0MultiplicityForSelEvt) {\r
 284       fHistV0MultiplicityForSelEvt = new TH1F("fHistV0MultiplicityForSelEvt", \r
 285                                               "V0s per event;Nbr of V0s/Evt;Events", \r
 286                                               25, 0, 25);\r
 287       fHistListV0->Add(fHistV0MultiplicityForSelEvt);\r
 288     }\r
 289     \r
 290     if(! fHistV0MultiplicityForSelEvtNoTPCOnly) {\r
 291     fHistV0MultiplicityForSelEvtNoTPCOnly = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnly", \r
 292                                                      "V0s per event;Nbr of V0s/Evt;Events", \r
 293                                                      25, 0, 25);\r
 294     fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnly);\r
 295     }\r
 296     \r
 297     if(! fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup) {\r
 298       fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup = new TH1F("fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup", \r
 299                                                              "V0s per event;Nbr of V0s/Evt;Events", \r
 300                                                                25, 0, 25);\r
 301       fHistListV0->Add(fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup);\r
 302     }\r
 303     \r
 304     //------------------------------------------------\r
 305     // Track Multiplicity Histograms\r
 306     //------------------------------------------------\r
 307     \r
 308     if(! fHistMultiplicityBeforeTrigSel) {\r
 309       fHistMultiplicityBeforeTrigSel = new TH1F("fHistMultiplicityBeforeTrigSel", \r
 310                                               "Tracks per event;Nbr of Tracks;Events", \r
 311                                                 200, 0, 200);           \r
 312       fHistListV0->Add(fHistMultiplicityBeforeTrigSel);\r
 313     }\r
 314     if(! fHistMultiplicityForTrigEvt) {\r
 315       fHistMultiplicityForTrigEvt = new TH1F("fHistMultiplicityForTrigEvt", \r
 316                                              "Tracks per event;Nbr of Tracks;Events", \r
 317                                              200, 0, 200);              \r
 318       fHistListV0->Add(fHistMultiplicityForTrigEvt);\r
 319     }\r
 320     if(! fHistMultiplicity) {\r
 321       fHistMultiplicity = new TH1F("fHistMultiplicity", \r
 322                                    "Tracks per event;Nbr of Tracks;Events", \r
 323                                    200, 0, 200);                \r
 324       fHistListV0->Add(fHistMultiplicity);\r
 325     }\r
 326     if(! fHistMultiplicityNoTPCOnly) {\r
 327       fHistMultiplicityNoTPCOnly = new TH1F("fHistMultiplicityNoTPCOnly", \r
 328                                             "Tracks per event;Nbr of Tracks;Events", \r
 329                                             200, 0, 200);               \r
 330     fHistListV0->Add(fHistMultiplicityNoTPCOnly);\r
 331     }\r
 332     if(! fHistMultiplicityNoTPCOnlyNoPileup) {\r
 333       fHistMultiplicityNoTPCOnlyNoPileup = new TH1F("fHistMultiplicityNoTPCOnlyNoPileup", \r
 334                                                     "Tracks per event;Nbr of Tracks;Events", \r
 335                                                     200, 0, 200);               \r
 336       fHistListV0->Add(fHistMultiplicityNoTPCOnlyNoPileup);\r
 337     }\r
 338 \r
 339     //------------------------------------------------\r
 340     // V0 selection Histograms (before)\r
 341     //------------------------------------------------\r
 342     if(!fHistV0InvMassK0) {\r
 343       fHistV0InvMassK0 = new TH1F("fHistV0InvMassK0",\r
 344                                   "Invariant Mass for K0;Mass (GeV/c^{2});Events",\r
 345                                   200,0,2);\r
 346       fHistListV0->Add(fHistV0InvMassK0);\r
 347     }\r
 348     if(!fHistV0InvMassLambda) {\r
 349       fHistV0InvMassLambda = new TH1F("fHistV0InvMassLambda",\r
 350                                   "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",\r
 351                                   200,0,2);\r
 352       fHistListV0->Add(fHistV0InvMassLambda);\r
 353     }\r
 354     if(!fHistV0InvMassAntiLambda) {\r
 355       fHistV0InvMassAntiLambda = new TH1F("fHistV0InvMassAntiLambda",\r
 356                                   "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",\r
 357                                   200,0,2);\r
 358       fHistListV0->Add(fHistV0InvMassAntiLambda);\r
 359     }\r
 360     if(!fHistV0Armenteros) {\r
 361       fHistV0Armenteros = new TH2F("fHistV0Armenteros",\r
 362                                   "Armenteros plot;#alpha;q_{t}",\r
 363                                    200,-1,1,200,0,0.5);\r
 364       fHistListV0->Add(fHistV0Armenteros);\r
 365     }\r
 366     \r
 367     //------------------------------------------------\r
 368     // V0 selection Histograms (after)\r
 369     //------------------------------------------------\r
 370     if(!fHistV0SelInvMassK0) {\r
 371       fHistV0SelInvMassK0 = new TH1F("fHistV0SelInvMassK0",\r
 372                                   "Invariant Mass for K0;Mass (GeV/c^{2});Events",\r
 373                                   200,0,2);\r
 374       fHistListV0->Add(fHistV0SelInvMassK0);\r
 375     }\r
 376     if(!fHistV0SelInvMassLambda) {\r
 377       fHistV0SelInvMassLambda = new TH1F("fHistV0SelInvMassLambda",\r
 378                                   "Invariant Mass for Lambda;Mass (GeV/c^{2});Events",\r
 379                                   200,0,2);\r
 380       fHistListV0->Add(fHistV0SelInvMassLambda);\r
 381     }\r
 382     if(!fHistV0SelInvMassAntiLambda) {\r
 383       fHistV0SelInvMassAntiLambda = new TH1F("fHistV0SelInvMassAntiLambda",\r
 384                                   "Invariant Mass for AntiLambda;Mass (GeV/c^{2});Events",\r
 385                                   200,0,2);\r
 386       fHistListV0->Add(fHistV0SelInvMassAntiLambda);\r
 387     }\r
 388     if(!fHistV0SelArmenteros) {\r
 389       fHistV0SelArmenteros = new TH2F("fHistV0SelArmenteros",\r
 390                                   "Armenteros plot;#alpha;q_{t}",\r
 391                                    200,-1,1,200,0,0.5);\r
 392       fHistListV0->Add(fHistV0SelArmenteros);\r
 393     }\r
 394   }//V0\r
 395     \r
 396   // Balance function histograms\r
 397   // Initialize histograms if not done yet\r
 398   if(!fBalance->GetHistNp()){\r
 399     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 400     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 401     fBalance->InitHistograms();\r
 402   }\r
 403 \r
 404   if(fRunShuffling) {\r
 405     if(!fShuffledBalance->GetHistNp()) {\r
 406       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 407       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 408       fShuffledBalance->InitHistograms();\r
 409     }\r
 410   }\r
 411 \r
 412   if(fRunMixing) {\r
 413     if(!fMixedBalance->GetHistNp()) {\r
 414       AliWarning("Histograms (mixing) not yet initialized! --> Will be done now");\r
 415       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 416       fMixedBalance->InitHistograms();\r
 417     }\r
 418   }\r
 419 \r
 420   fListTriggeredBF->Add(fBalance->GetHistNp());\r
 421   fListTriggeredBF->Add(fBalance->GetHistNn());\r
 422   fListTriggeredBF->Add(fBalance->GetHistNpn());\r
 423   fListTriggeredBF->Add(fBalance->GetHistNnn());\r
 424   fListTriggeredBF->Add(fBalance->GetHistNpp());\r
 425   fListTriggeredBF->Add(fBalance->GetHistNnp());\r
 426   \r
 427   if(fRunShuffling) {\r
 428     fListTriggeredBFS->Add(fShuffledBalance->GetHistNp());\r
 429     fListTriggeredBFS->Add(fShuffledBalance->GetHistNn());\r
 430     fListTriggeredBFS->Add(fShuffledBalance->GetHistNpn());\r
 431     fListTriggeredBFS->Add(fShuffledBalance->GetHistNnn());\r
 432     fListTriggeredBFS->Add(fShuffledBalance->GetHistNpp());\r
 433     fListTriggeredBFS->Add(fShuffledBalance->GetHistNnp());\r
 434   }  \r
 435 \r
 436   if(fRunMixing) {\r
 437     fListTriggeredBFM->Add(fMixedBalance->GetHistNp());\r
 438     fListTriggeredBFM->Add(fMixedBalance->GetHistNn());\r
 439     fListTriggeredBFM->Add(fMixedBalance->GetHistNpn());\r
 440     fListTriggeredBFM->Add(fMixedBalance->GetHistNnn());\r
 441     fListTriggeredBFM->Add(fMixedBalance->GetHistNpp());\r
 442     fListTriggeredBFM->Add(fMixedBalance->GetHistNnp());\r
 443   }  \r
 444 \r
 445   // PID Response task active?\r
 446   if(fRunV0) {\r
 447     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
 448     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 449   }\r
 450 \r
 451   // Event Mixing\r
 452   Int_t trackDepth = fMixingTracks; \r
 453   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
 454    \r
 455   Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,90.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 456   Double_t* centbins        = centralityBins;\r
 457   Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
 458   \r
 459   // bins for second buffer are shifted by 100 cm\r
 460   Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 461   Double_t* vtxbins     = vertexBins;\r
 462   Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
 463 \r
 464   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 465 \r
 466 \r
 467   // Post output data.\r
 468   PostData(1, fList);\r
 469   PostData(2, fListTriggeredBF);\r
 470   if(fRunShuffling) PostData(3, fListTriggeredBFS);\r
 471   if(fRunMixing) PostData(4, fListTriggeredBFM);\r
 472   if(fRunV0) PostData(5,fHistListV0);\r
 473 \r
 474   TH1::AddDirectory(oldStatus);\r
 475 \r
 476 }\r
 477 \r
 478 //________________________________________________________________________\r
 479 void AliAnalysisTaskTriggeredBF::UserExec(Option_t *) {\r
 480   // Main loop\r
 481   // Called for each event\r
 482 \r
 483   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 484   Float_t fCentrality = -1.;  \r
 485 \r
 486   // -------------------------------------------------------------                   \r
 487   // AOD analysis (vertex and track cuts also here!!!!)\r
 488   if(gAnalysisLevel == "AOD") {\r
 489     AliVEvent* eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
 490     if(!eventMain) {\r
 491       AliError("eventMain not available");\r
 492       return;\r
 493     }\r
 494 \r
 495     // check event cuts and fill event histograms\r
 496     if((fCentrality = IsEventAccepted(eventMain)) < 0){\r
 497       return;\r
 498     }\r
 499     \r
 500     // get the accepted tracks in main event\r
 501     TObjArray *tracksMain = NULL;\r
 502     if(fRunV0) tracksMain = GetAcceptedV0s(eventMain);\r
 503     else       tracksMain = GetAcceptedTracks(eventMain);\r
 504 \r
 505     // store charges of all accepted tracks, shuffle and reassign (two extra loops!)\r
 506     TObjArray* tracksShuffled = NULL;\r
 507     if(fRunShuffling){\r
 508       tracksShuffled = GetShuffledTracks(tracksMain);\r
 509     }\r
 510     \r
 511     // Event mixing --> UPDATE POOL IS MISSING!!!\r
 512     if (fRunMixing)\r
 513       {\r
 514         // 1. First get an event pool corresponding in mult (cent) and\r
 515         //    zvertex to the current event. Once initialized, the pool\r
 516         //    should contain nMix (reduced) events. This routine does not\r
 517         //    pre-scan the chain. The first several events of every chain\r
 518         //    will be skipped until the needed pools are filled to the\r
 519         //    specified depth. If the pool categories are not too rare, this\r
 520         //    should not be a problem. If they are rare, you could lose`\r
 521         //    statistics.\r
 522         \r
 523         // 2. Collect the whole pool's content of tracks into one TObjArray\r
 524         //    (bgTracks), which is effectively a single background super-event.\r
 525         \r
 526         // 3. The reduced and bgTracks arrays must both be passed into\r
 527         //    FillCorrelations(). Also nMix should be passed in, so a weight\r
 528         //    of 1./nMix can be applied.\r
 529         \r
 530         AliEventPool* pool = fPoolMgr->GetEventPool(fCentrality, eventMain->GetPrimaryVertex()->GetZ());\r
 531         \r
 532         if (!pool){\r
 533           AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, eventMain->GetPrimaryVertex()->GetZ()));\r
 534         }\r
 535         else{\r
 536 \r
 537         //pool->SetDebug(1);\r
 538         \r
 539           if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
 540             \r
 541             \r
 542             Int_t nMix = pool->GetCurrentNEvents();\r
 543             //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
 544             \r
 545             //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
 546             //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
 547             //if (pool->IsReady())\r
 548             //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
 549             \r
 550             // Fill mixed-event histos here  \r
 551             for (Int_t jMix=0; jMix<nMix; jMix++) \r
 552               {\r
 553                 TObjArray* tracksMixed = pool->GetEvent(jMix);\r
 554                 fMixedBalance->FillBalance(fCentrality,tracksMain,tracksMixed); \r
 555               }\r
 556           }\r
 557           \r
 558           // Update the Event pool\r
 559           pool->UpdatePool(tracksMain);\r
 560           //pool->PrintInfo();\r
 561           \r
 562         }//pool NULL check  \r
 563       }//run mixing\r
 564     \r
 565     // calculate balance function\r
 566     fBalance->FillBalance(fCentrality,tracksMain,NULL);\r
 567     \r
 568     // calculate shuffled balance function\r
 569     if(fRunShuffling && tracksShuffled != NULL) {\r
 570       fShuffledBalance->FillBalance(fCentrality,tracksShuffled,NULL);\r
 571     }\r
 572     \r
 573   }//AOD analysis\r
 574   else{\r
 575     AliError("Triggered Balance Function analysis only for AODs!");\r
 576   }\r
 577 }     \r
 578 \r
 579 //________________________________________________________________________\r
 580 Float_t AliAnalysisTaskTriggeredBF::IsEventAccepted(AliVEvent *event){\r
 581   // Checks the Event cuts\r
 582   // Fills Event statistics histograms\r
 583   \r
 584   // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
 585   fHistEventStats->Fill(1); //all events\r
 586 \r
 587   Bool_t isSelectedMain = kTRUE;\r
 588   Float_t fCentrality = -1.;\r
 589   Int_t nV0s          = event->GetNumberOfV0s();\r
 590   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 591   \r
 592   if(fUseOfflineTrigger)\r
 593     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 594   \r
 595   //V0 QA histograms (before trigger selection)\r
 596   if(fRunV0){\r
 597     fHistMultiplicityBeforeTrigSel->Fill ( -1 );\r
 598     fHistV0MultiplicityBeforeTrigSel->Fill ( nV0s );\r
 599   }\r
 600   \r
 601   if(isSelectedMain) {\r
 602     fHistEventStats->Fill(2); //triggered events\r
 603     \r
 604     //Centrality stuff \r
 605     if(fUseCentrality) {\r
 606       if(gAnalysisLevel == "AOD") { //centrality in AOD header\r
 607         AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 608         fCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 609 \r
 610         // QA for centrality estimators\r
 611         fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 612         fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 613         fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 614         fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("TKL"));\r
 615         fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 616         fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 617         fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 618         fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 619         fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 620         \r
 621         // centrality QA (V0M)\r
 622         fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 623         \r
 624         // centrality QA (reference tracks)\r
 625         fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
 626         fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
 627         fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
 628         fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
 629         fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
 630         fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
 631         fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
 632         fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
 633         fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
 634 \r
 635         //V0 QA histograms (after trigger selection)\r
 636         if(fRunV0){\r
 637           fHistMultiplicityForTrigEvt->Fill ( fCentrality );\r
 638           fHistV0MultiplicityForTrigEvt->Fill ( nV0s );\r
 639         }\r
 640       }\r
 641     }\r
 642     \r
 643     \r
 644     const AliVVertex *vertex = event->GetPrimaryVertex();\r
 645     \r
 646     if(vertex) {\r
 647       Double32_t fCov[6];\r
 648       vertex->GetCovarianceMatrix(fCov);\r
 649       if(vertex->GetNContributors() > 0) {\r
 650         if(fCov[5] != 0) {\r
 651           fHistEventStats->Fill(3); //events with a proper vertex\r
 652           if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 653             if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 654               if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 655                 fHistEventStats->Fill(4); //analyzed events\r
 656                 fHistVx->Fill(vertex->GetX());\r
 657                 fHistVy->Fill(vertex->GetY());\r
 658                 fHistVz->Fill(vertex->GetZ());\r
 659 \r
 660 \r
 661 \r
 662                 //V0 QA histograms (vertex Z check)\r
 663                 if(fRunV0){\r
 664                   fHistV0MultiplicityForSelEvt ->Fill( nV0s );\r
 665                   fHistMultiplicity->Fill(fCentrality);\r
 666 \r
 667                   //V0 QA histograms (Only look at events with well-established PV)\r
 668                   const AliAODVertex *lPrimarySPDVtx = ((AliAODEvent*)event)->GetPrimaryVertexSPD();\r
 669                   if(lPrimarySPDVtx){\r
 670                     fHistMultiplicityNoTPCOnly->Fill ( fCentrality );\r
 671                     fHistV0MultiplicityForSelEvtNoTPCOnly->Fill ( nV0s );\r
 672                     \r
 673                     //V0 QA histograms (Pileup Rejection)\r
 674                     // FIXME : quality selection regarding pile-up rejection \r
 675                     fHistMultiplicityNoTPCOnlyNoPileup->Fill(fCentrality);\r
 676                     fHistV0MultiplicityForSelEvtNoTPCOnlyNoPileup ->Fill( nV0s );\r
 677 \r
 678                   }\r
 679                   else{\r
 680                     return -1;\r
 681                   }\r
 682                 }\r
 683                 \r
 684                 \r
 685                 // take only events inside centrality class\r
 686                 if((fCentrality > fCentralityPercentileMin) && (fCentrality < fCentralityPercentileMax)){\r
 687                   return fCentrality;           \r
 688                 }//centrality class\r
 689               }//Vz cut\r
 690             }//Vy cut\r
 691           }//Vx cut\r
 692         }//proper vertex resolution\r
 693       }//proper number of contributors\r
 694     }//vertex object valid\r
 695   }//triggered event \r
 696   \r
 697   // in all other cases return -1 (event not accepted)\r
 698   return -1;\r
 699 }\r
 700 \r
 701 //________________________________________________________________________\r
 702 TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedTracks(AliVEvent *event){\r
 703   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
 704   // Fills QA histograms\r
 705 \r
 706   //output TObjArray holding all good tracks\r
 707   TObjArray* tracksAccepted = new TObjArray;\r
 708   tracksAccepted->SetOwner(kTRUE);\r
 709 \r
 710   Short_t vCharge = 0;\r
 711   Double_t vEta    = 0.;\r
 712   Double_t vPhi    = 0.;\r
 713   Double_t vPt     = 0.;\r
 714   \r
 715   // Loop over tracks in event\r
 716   for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
 717     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
 718     if (!aodTrack) {\r
 719       AliError(Form("Could not receive track %d", iTracks));\r
 720       continue;\r
 721     }\r
 722     \r
 723     // AOD track cuts\r
 724     \r
 725     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
 726     // take only TPC only tracks \r
 727     fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
 728     if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
 729     \r
 730     vCharge = aodTrack->Charge();\r
 731     vEta    = aodTrack->Eta();\r
 732     vPhi    = aodTrack->Phi() * TMath::RadToDeg();\r
 733     vPt     = aodTrack->Pt();\r
 734     \r
 735     Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
 736     Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
 737     \r
 738     \r
 739     // Kinematics cuts from ESD track cuts\r
 740     if( vPt < fPtMin || vPt > fPtMax)      continue;\r
 741     if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
 742     \r
 743     // Extra DCA cuts (for systematic studies [!= -1])\r
 744     if( fDCAxyCut != -1 && fDCAzCut != -1){\r
 745       if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
 746         continue;  // 2D cut\r
 747       }\r
 748     }\r
 749     \r
 750     // Extra TPC cuts (for systematic studies [!= -1])\r
 751     if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
 752       continue;\r
 753     }\r
 754     if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
 755       continue;\r
 756     }\r
 757     \r
 758     // fill QA histograms\r
 759     fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
 760     fHistDCA->Fill(dcaZ,dcaXY);\r
 761     fHistChi2->Fill(aodTrack->Chi2perNDF());\r
 762     fHistPt->Fill(vPt);\r
 763     fHistEta->Fill(vEta);\r
 764     fHistPhi->Fill(vPhi);\r
 765     \r
 766     // add the track to the TObjArray\r
 767     tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge,1.));\r
 768   }\r
 769 \r
 770   return tracksAccepted;\r
 771 }\r
 772 \r
 773 //________________________________________________________________________\r
 774 TObjArray* AliAnalysisTaskTriggeredBF::GetAcceptedV0s(AliVEvent *event){\r
 775   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
 776   // Fills QA histograms\r
 777 \r
 778   //output TObjArray holding all good tracks\r
 779   TObjArray* tracksAccepted = new TObjArray;\r
 780   tracksAccepted->SetOwner(kTRUE);\r
 781 \r
 782   Short_t vCharge = 0;\r
 783   Double_t vEta    = 0.;\r
 784   Double_t vPhi    = 0.;\r
 785   Double_t vPt     = 0.;\r
 786   \r
 787   //------------------------------------------------\r
 788   // MAIN LAMBDA LOOP STARTS HERE (basically a copy of AliAnalysisTaskExtractV0AOD)\r
 789   //------------------------------------------------\r
 790 \r
 791   // parameters (for the time being hard coded here) --> from David for EbyE Lambdas\r
 792   Bool_t fkUseOnTheFly = kFALSE;\r
 793   Double_t fRapidityBoundary  = 0.5; \r
 794   Double_t fCutDaughterEta    = 0.8;\r
 795   Double_t fCutV0Radius       = 0.9;\r
 796   Double_t fCutDCANegToPV     = 0.1;\r
 797   Double_t fCutDCAPosToPV     = 0.1;\r
 798   Double_t fCutDCAV0Daughters = 1.0;\r
 799   Double_t fCutV0CosPA        = 0.9995;\r
 800   Double_t fMassLambda        = 1.115683;\r
 801   Double_t fCutMassLambda     = 0.007;\r
 802   Double_t fCutProperLifetime = 3*7.9;\r
 803   Double_t fCutLeastNumberOfCrossedRows = 70;\r
 804   Double_t fCutLeastNumberOfCrossedRowsOverFindable = 0.8;\r
 805   Double_t fCutTPCPIDNSigmasProton  = 3.0;\r
 806   Double_t fCutTPCPIDNSigmasPion    = 5.0;\r
 807 \r
 808 \r
 809   //Variable definition\r
 810   Int_t    lOnFlyStatus = 0;// nv0sOn = 0, nv0sOff = 0;\r
 811   Double_t lChi2V0 = 0;\r
 812   Double_t lDcaV0Daughters = 0, lDcaV0ToPrimVertex = 0;\r
 813   Double_t lDcaPosToPrimVertex = 0, lDcaNegToPrimVertex = 0;\r
 814   Double_t lV0CosineOfPointingAngle = 0;\r
 815   Double_t lV0Radius = 0, lPt = 0;\r
 816   Double_t lEta = 0, lPhi = 0;\r
 817   Double_t lRap = 0, lRapK0Short = 0, lRapLambda = 0;\r
 818   Double_t lInvMassK0s = 0, lInvMassLambda = 0, lInvMassAntiLambda = 0;\r
 819   Double_t lAlphaV0 = 0, lPtArmV0 = 0;\r
 820   \r
 821   Double_t fMinV0Pt = 0; \r
 822   Double_t fMaxV0Pt = 100; \r
 823   \r
 824 \r
 825   \r
 826   // some event observables\r
 827   Int_t nv0s = event->GetNumberOfV0s();\r
 828   Double_t tPrimaryVtxPosition[3];\r
 829   const AliVVertex *primaryVtx = event->GetPrimaryVertex();\r
 830   tPrimaryVtxPosition[0] = primaryVtx->GetX();\r
 831   tPrimaryVtxPosition[1] = primaryVtx->GetY();\r
 832   tPrimaryVtxPosition[2] = primaryVtx->GetZ();\r
 833 \r
 834 \r
 835   //loop over V0s  \r
 836   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) \r
 837     {// This is the begining of the V0 loop\r
 838       AliAODv0 *v0 = ((AliAODEvent*)event)->GetV0(iV0);\r
 839       if (!v0) continue;\r
 840 \r
 841       //Obsolete at AOD level... \r
 842       //---> Fix On-the-Fly candidates, count how many swapped\r
 843       //if( v0->GetParamN()->Charge() > 0 && v0->GetParamP()->Charge() < 0 ){\r
 844       //  fHistSwappedV0Counter -> Fill( 1 );\r
 845       //}else{\r
 846       //  fHistSwappedV0Counter -> Fill( 0 ); \r
 847       //}\r
 848       //if ( fkUseOnTheFly ) CheckChargeV0(v0); \r
 849       \r
 850       Double_t tDecayVertexV0[3]; v0->GetXYZ(tDecayVertexV0); \r
 851       Double_t tV0mom[3];\r
 852       v0->GetPxPyPz( tV0mom ); \r
 853       Double_t lV0TotalMomentum = TMath::Sqrt(\r
 854                                               tV0mom[0]*tV0mom[0]+tV0mom[1]*tV0mom[1]+tV0mom[2]*tV0mom[2] );\r
 855       \r
 856       lV0Radius = TMath::Sqrt(tDecayVertexV0[0]*tDecayVertexV0[0]+tDecayVertexV0[1]*tDecayVertexV0[1]);\r
 857       lPt = v0->Pt();\r
 858       lEta = v0->Eta();\r
 859       lPhi = v0->Phi()*TMath::RadToDeg();\r
 860       lRapK0Short = v0->RapK0Short();\r
 861       lRapLambda  = v0->RapLambda();\r
 862       lRap        = lRapLambda;//v0->Y(); //FIXME!!!\r
 863       if ((lPt<fMinV0Pt)||(fMaxV0Pt<lPt)) continue;\r
 864       \r
 865       //UInt_t lKeyPos = (UInt_t)TMath::Abs(v0->GetPosID());\r
 866       //UInt_t lKeyNeg = (UInt_t)TMath::Abs(v0->GetPosID());\r
 867 \r
 868       Double_t lMomPos[3]; //v0->GetPPxPyPz(lMomPos[0],lMomPos[1],lMomPos[2]);\r
 869       Double_t lMomNeg[3]; //v0->GetNPxPyPz(lMomNeg[0],lMomNeg[1],lMomNeg[2]);\r
 870       lMomPos[0] = v0->MomPosX();\r
 871       lMomPos[1] = v0->MomPosY();\r
 872       lMomPos[2] = v0->MomPosZ();\r
 873       lMomNeg[0] = v0->MomNegX();\r
 874       lMomNeg[1] = v0->MomNegY();\r
 875       lMomNeg[2] = v0->MomNegZ();\r
 876       \r
 877       AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter\r
 878       AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter\r
 879       if (!pTrack || !nTrack) {\r
 880         AliError("ERROR: Could not retreive one of the daughter track");\r
 881         continue;\r
 882       }\r
 883 \r
 884       //Daughter Eta for Eta selection, afterwards\r
 885       Double_t lNegEta = nTrack->Eta();\r
 886       Double_t lPosEta = pTrack->Eta();\r
 887       \r
 888       // Filter like-sign V0 (next: add counter and distribution)\r
 889       if ( pTrack->Charge() == nTrack->Charge()){\r
 890         continue;\r
 891       } \r
 892       \r
 893       //Quick test this far! \r
 894       \r
 895 \r
 896       //________________________________________________________________________\r
 897       // Track quality cuts \r
 898       Float_t lPosTrackCrossedRows = pTrack->GetTPCClusterInfo(2,1);\r
 899       Float_t lNegTrackCrossedRows = nTrack->GetTPCClusterInfo(2,1);\r
 900       Float_t lLeastNbrCrossedRows =  (lPosTrackCrossedRows>lNegTrackCrossedRows) ? lNegTrackCrossedRows : lPosTrackCrossedRows;\r
 901 \r
 902       // TPC refit condition (done during reconstruction for Offline but not for On-the-fly)\r
 903       if( !(pTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;\r
 904       if( !(nTrack->GetStatus() & AliESDtrack::kTPCrefit)) continue;\r
 905       \r
 906       if ( ( ( pTrack->GetTPCClusterInfo(2,1) ) < 70 ) || ( ( nTrack->GetTPCClusterInfo(2,1) ) < 70 ) ) continue;\r
 907       \r
 908       //Findable clusters > 0 condition\r
 909       if( pTrack->GetTPCNclsF()<=0 || nTrack->GetTPCNclsF()<=0 ) continue;\r
 910       \r
 911       //Compute ratio Crossed Rows / Findable clusters\r
 912       //Note: above test avoids division by zero! \r
 913       Float_t lPosTrackCrossedRowsOverFindable = lPosTrackCrossedRows / ((double)(pTrack->GetTPCNclsF())); \r
 914       Float_t lNegTrackCrossedRowsOverFindable = lNegTrackCrossedRows / ((double)(nTrack->GetTPCNclsF())); \r
 915       Float_t lLeastNbrCrossedRowsOverFindable = (lPosTrackCrossedRowsOverFindable>lNegTrackCrossedRowsOverFindable) ? lNegTrackCrossedRowsOverFindable : lPosTrackCrossedRowsOverFindable;\r
 916 \r
 917       //Lowest Cut Level for Ratio Crossed Rows / Findable = 0.8, set here\r
 918       if ( lLeastNbrCrossedRowsOverFindable < 0.8) continue;\r
 919       \r
 920       //End track Quality Cuts\r
 921       //________________________________________________________________________\r
 922       \r
 923       \r
 924       lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();\r
 925       lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();\r
 926           \r
 927       lOnFlyStatus = v0->GetOnFlyStatus();\r
 928       lChi2V0 = v0->Chi2V0();\r
 929       lDcaV0Daughters = v0->DcaV0Daughters();\r
 930       lDcaV0ToPrimVertex = v0->DcaV0ToPrimVertex();\r
 931       lV0CosineOfPointingAngle = v0->CosPointingAngle(tPrimaryVtxPosition);\r
 932       \r
 933       // Distance over total momentum\r
 934       Double_t lDistOverTotMom = TMath::Sqrt(\r
 935                                     TMath::Power( tDecayVertexV0[0] - tPrimaryVtxPosition[0] , 2) +\r
 936                                     TMath::Power( tDecayVertexV0[1] - tPrimaryVtxPosition[1] , 2) +\r
 937                                     TMath::Power( tDecayVertexV0[2] - tPrimaryVtxPosition[2] , 2)\r
 938                                     );\r
 939       lDistOverTotMom /= (lV0TotalMomentum+1e-10); //avoid division by zero, to be sure\r
 940       \r
 941       \r
 942       // Getting invariant mass infos directly from ESD\r
 943       lInvMassK0s        = v0->MassK0Short();\r
 944       lInvMassLambda     = v0->MassLambda();\r
 945       lInvMassAntiLambda = v0->MassAntiLambda();\r
 946       lAlphaV0 = v0->AlphaV0();\r
 947       lPtArmV0 = v0->PtArmV0();\r
 948 \r
 949       //Official means of acquiring N-sigmas \r
 950       Double_t lNSigmasPosProton = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kProton );\r
 951       Double_t lNSigmasPosPion   = fPIDResponse->NumberOfSigmasTPC( pTrack, AliPID::kPion );\r
 952       Double_t lNSigmasNegProton = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kProton );\r
 953       Double_t lNSigmasNegPion   = fPIDResponse->NumberOfSigmasTPC( nTrack, AliPID::kPion );\r
 954 \r
 955       //V0 QA histograms (before V0 selection)\r
 956       fHistV0InvMassK0->Fill(lInvMassK0s);\r
 957       fHistV0InvMassLambda->Fill(lInvMassLambda);\r
 958       fHistV0InvMassAntiLambda->Fill(lInvMassAntiLambda);\r
 959       fHistV0Armenteros->Fill(lAlphaV0,lPtArmV0);\r
 960       \r
 961       \r
 962       //First Selection: Reject OnFly\r
 963       if( (lOnFlyStatus == 0 && fkUseOnTheFly == kFALSE) || (lOnFlyStatus != 0 && fkUseOnTheFly == kTRUE ) ){\r
 964         \r
 965 \r
 966         //Second Selection: rough 20-sigma band, parametric.    \r
 967         //K0Short: Enough to parametrize peak broadening with linear function.    \r
 968         Double_t lUpperLimitK0Short = (5.63707e-01) + (1.14979e-02)*lPt; \r
 969         Double_t lLowerLimitK0Short = (4.30006e-01) - (1.10029e-02)*lPt;\r
 970         \r
 971         //Lambda: Linear (for higher pt) plus exponential (for low-pt broadening)\r
 972         //[0]+[1]*x+[2]*TMath::Exp(-[3]*x)\r
 973         Double_t lUpperLimitLambda = (1.13688e+00) + (5.27838e-03)*lPt + (8.42220e-02)*TMath::Exp(-(3.80595e+00)*lPt); \r
 974         Double_t lLowerLimitLambda = (1.09501e+00) - (5.23272e-03)*lPt - (7.52690e-02)*TMath::Exp(-(3.46339e+00)*lPt);\r
 975         \r
 976         //Do Selection      \r
 977         if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda     ) || \r
 978             (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda     ) || \r
 979             (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short    ) ){\r
 980 \r
 981 \r
 982       //          //Pre-selection in case this is AA...\r
 983       //          //if( fkIsNuclear == kFALSE ) fTree->Fill();\r
 984       //          //if( fkIsNuclear == kTRUE){ \r
 985       //          //If this is a nuclear collision___________________\r
 986       //          // ... pre-filter with TPC, daughter eta selection\r
 987 \r
 988           \r
 989           if( (lInvMassLambda     < lUpperLimitLambda  && lInvMassLambda     > lLowerLimitLambda \r
 990                && TMath::Abs(lNSigmasPosProton) < 6.0 && TMath::Abs(lNSigmasNegPion) < 6.0 ) || \r
 991               (lInvMassAntiLambda < lUpperLimitLambda  && lInvMassAntiLambda > lLowerLimitLambda \r
 992                && TMath::Abs(lNSigmasNegProton) < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ||  \r
 993               (lInvMassK0s        < lUpperLimitK0Short && lInvMassK0s        > lLowerLimitK0Short \r
 994                && TMath::Abs(lNSigmasNegPion)   < 6.0 && TMath::Abs(lNSigmasPosPion) < 6.0 ) ){\r
 995             \r
 996             //insane test\r
 997             if ( TMath::Abs(lNegEta)<0.8 && TMath::Abs(lPosEta)<0.8 ){\r
 998 \r
 999               // start the fine selection (usually done in post processing, but we don't have time to waste) --> Lambdas!\r
1000               if(\r
1001                  TMath::Abs(lRap)<fRapidityBoundary &&\r
1002                  TMath::Abs(lNegEta)       <= fCutDaughterEta               &&                   \r
1003                  TMath::Abs(lPosEta)       <= fCutDaughterEta               &&\r
1004                  lV0Radius                 >= fCutV0Radius                  &&\r
1005                  lDcaNegToPrimVertex       >= fCutDCANegToPV                &&\r
1006                  lDcaPosToPrimVertex       >= fCutDCAPosToPV                &&\r
1007                  lDcaV0Daughters           <= fCutDCAV0Daughters            &&\r
1008                  lV0CosineOfPointingAngle  >= fCutV0CosPA                   && \r
1009                  fMassLambda*lDistOverTotMom    <= fCutProperLifetime       &&\r
1010                  lLeastNbrCrossedRows             >= fCutLeastNumberOfCrossedRows             &&\r
1011                  lLeastNbrCrossedRowsOverFindable >= fCutLeastNumberOfCrossedRowsOverFindable &&\r
1012                  lPtArmV0 * 5 < TMath::Abs(lAlphaV0)                        && \r
1013                  ((TMath::Abs(lNSigmasNegPion)   <= fCutTPCPIDNSigmasPion     &&\r
1014                   TMath::Abs(lNSigmasPosProton) <= fCutTPCPIDNSigmasProton) ||\r
1015                   (TMath::Abs(lNSigmasPosPion)   <= fCutTPCPIDNSigmasPion     &&\r
1016                    TMath::Abs(lNSigmasNegProton) <= fCutTPCPIDNSigmasProton))            \r
1017                  )\r
1018                 {\r
1019 \r
1020                   //V0 QA histograms (after V0 selection)\r
1021                   fHistV0SelInvMassK0->Fill(lInvMassK0s);\r
1022                   fHistV0SelInvMassLambda->Fill(lInvMassLambda);\r
1023                   fHistV0SelInvMassAntiLambda->Fill(lInvMassAntiLambda);\r
1024 \r
1025                   // this means a V0 candidate is found\r
1026                   if(TMath::Abs(lInvMassLambda-fMassLambda) < fCutMassLambda ||\r
1027                      TMath::Abs(lInvMassAntiLambda-fMassLambda) < fCutMassLambda){\r
1028 \r
1029                     fHistV0SelArmenteros->Fill(lAlphaV0,lPtArmV0);                \r
1030 \r
1031                     vEta    = lEta;\r
1032                     vPhi    = lPhi;\r
1033                     vPt     = lPt;\r
1034                     if(lAlphaV0 > 0) vCharge = 1;\r
1035                     if(lAlphaV0 < 0) vCharge = -1;\r
1036 \r
1037                     // fill QA histograms\r
1038                     fHistPt->Fill(vPt);\r
1039                     fHistEta->Fill(vEta);\r
1040                     fHistPhi->Fill(vPhi);\r
1041                     \r
1042                     // add the track to the TObjArray\r
1043                     tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge,1.));\r
1044                   }\r
1045                 }\r
1046               }\r
1047           }\r
1048           //}//end nuclear_____________________________________\r
1049         }\r
1050       }\r
1051     }//V0 loop\r
1052   \r
1053   return tracksAccepted;\r
1054 }\r
1055 \r
1056 //________________________________________________________________________\r
1057 TObjArray* AliAnalysisTaskTriggeredBF::GetShuffledTracks(TObjArray *tracks){\r
1058   // Clones TObjArray and returns it with tracks after shuffling the charges\r
1059 \r
1060   TObjArray* tracksShuffled = new TObjArray;\r
1061   tracksShuffled->SetOwner(kTRUE);\r
1062 \r
1063   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
1064 \r
1065   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1066   {\r
1067     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1068     chargeVector->push_back(track->Charge());\r
1069   }  \r
1070  \r
1071   random_shuffle(chargeVector->begin(), chargeVector->end());\r
1072   \r
1073   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1074     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1075     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i),1.));\r
1076   }\r
1077 \r
1078   delete chargeVector;\r
1079    \r
1080   return tracksShuffled;\r
1081 }\r
1082 \r
1083 //________________________________________________________________________\r
1084 void  AliAnalysisTaskTriggeredBF::FinishTaskOutput(){\r
1085   //checks if Balance Function objects are there (needed to write the histograms)\r
1086   if (!fBalance) {\r
1087     AliError("fBalance not available");\r
1088     return;\r
1089   }  \r
1090   if(fRunShuffling) {\r
1091     if (!fShuffledBalance) {\r
1092       AliError("fShuffledBalance not available");\r
1093       return;\r
1094     }\r
1095   }\r
1096 \r
1097 }\r
1098 \r
1099 //________________________________________________________________________\r
1100 void AliAnalysisTaskTriggeredBF::Terminate(Option_t *) {\r
1101   // Called once at the end of the query\r
1102 \r
1103   // not implemented ...\r
1104 \r
1105 }\r
1106 \r
1107 void AliAnalysisTaskTriggeredBF::UserExecMix(Option_t *)\r
1108 {\r
1109 \r
1110   // not yet done for event mixing!\r
1111   return;\r
1112 \r
1113 }\r
1114 \r