PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskTriggeredBF.cxx

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