PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // AliAnalysisTaskPhiFlow:
  17 // author: Redmer Alexander Bertens (rbertens@cern.ch)
  18 // analyis task for phi-meson reconstruction and determination of v_n
  19 // AliPhiMesonHelperTrack provides a lightweight helper track for reconstruction
  20 // new in this version (use with caution): vzero event plane, event mixing
  21
  22 #include "TChain.h"
  23 #include "TH1.h"
  24 #include "TH1F.h"
  25 #include "TH2F.h"
  26 #include "TMath.h"
  27 #include "TObjArray.h"
  28 #include "AliAnalysisTaskSE.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliAODEvent.h"
  31 #include "AliAODInputHandler.h"
  32 #include "AliCentrality.h"
  33 #include "AliAnalysisTaskPhiFlow.h"
  34 #include "AliFlowBayesianPID.h"
  35 #include "AliPIDCombined.h"
  36 #include "AliMCEvent.h"
  37 #include "TProfile.h"
  38 #include "AliFlowCandidateTrack.h"
  39 #include "AliFlowTrackCuts.h"
  40 #include "AliFlowEventSimple.h"
  41 #include "AliFlowTrackSimple.h"
  42 #include "AliFlowCommonConstants.h"
  43 #include "AliFlowEvent.h"
  44 #include "TVector3.h"
  45 #include "AliAODVZERO.h"
  46 #include "AliPIDResponse.h"
  47 #include "AliAODMCParticle.h"
  48 #include "AliAnalysisTaskVnV0.h"
  49 #include "AliEventPoolManager.h"
  50
  51 class AliFlowTrackCuts;
  52
  53 using std::cout;
  54 using std::endl;
  55
  56 ClassImp(AliAnalysisTaskPhiFlow)
  57 ClassImp(AliPhiMesonHelperTrack)
  58
  59 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
  60    fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCandidateYCut(kFALSE), fCandidateMinY(-.5), fCandidateMaxY(.5), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0),fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0)
  61 {
  62    // Default constructor
  63    for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
  64    for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
  65    for(Int_t i(0); i < 20; i++) {
  66        fVertexMixingBins[i] = 0;
  67        fCentralityMixingBins[i] = 0;
  68    }
  69    fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
  70    for(Int_t i(0); i < 18; i++) {
  71        for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
  72        fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
  73    }
  74 }
  75 //_____________________________________________________________________________
  76 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
  77    fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCandidateYCut(kFALSE), fCandidateMinY(-.5), fCandidateMaxY(.5), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0)
  78 {
  79    // Constructor
  80    for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
  81    for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
  82    for(Int_t i(0); i < 20; i++) {
  83        fVertexMixingBins[i] = 0;
  84        fCentralityMixingBins[i] = 0;
  85    }
  86    fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
  87    for(Int_t i(0); i < 18; i++) {
  88        for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
  89        fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
  90    }
  91    DefineInput(0, TChain::Class());
  92    DefineOutput(1, TList::Class());
  93    DefineOutput(2, AliFlowEventSimple::Class());
  94    if(fDebug > 0) cout << " === Created instance of AliAnaysisTaskPhiFlow === " << endl;
  95 }
  96 //_____________________________________________________________________________
  97 AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
  98 {
  99    // Destructor
 100    if (fNullCuts) delete fNullCuts;
 101    if (fOutputList) delete fOutputList;
 102    if (fCandidates) delete fCandidates;
 103    if (fFlowEvent) delete fFlowEvent;
 104    if (fBayesianResponse) delete fBayesianResponse;
 105    if (fEventMixing) delete fPoolManager;
 106    if (fPIDCombined) delete fPIDCombined;
 107    if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTaskPhiFlow === " << endl;
 108 }
 109 //_____________________________________________________________________________
 110 TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
 111 {
 112    // Return a pointer to a TH1 with predefined binning
 113    if(fDebug > 0) cout << " *** BookHistogram() *** " << name << endl;
 114    TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, .99, 1.092);
 115    hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
 116    hist->GetYaxis()->SetTitle("No. of pairs");
 117    hist->SetMarkerStyle(kFullCircle);
 118    hist->Sumw2();
 119    fOutputList->Add(hist);
 120    return hist;
 121 }
 122 //_____________________________________________________________________________
 123 TH2F* AliAnalysisTaskPhiFlow::BookPIDHistogram(const char* name, Bool_t TPC)
 124 {
 125    // Return a pointer to a TH2 with predefined binning
 126    if(fDebug > 0) cout << " *** BookPIDHisotgram() *** " << endl;
 127    TH2F *hist = 0x0;
 128    if(TPC) {
 129        hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
 130        hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
 131    }
 132    if(!TPC) {
 133        hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
 134        hist->GetYaxis()->SetTitle("#beta");
 135    }
 136    hist->GetXaxis()->SetTitle("P (GeV / c)");
 137    fOutputList->Add(hist);
 138    return hist;
 139 }
 140 //_____________________________________________________________________________
 141 TH1F* AliAnalysisTaskPhiFlow::InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
 142 {
 143    // intialize p_t histograms for each p_t bin
 144    if(fDebug > 0) cout << " *** InitPtSpectraHistograms() *** " << endl;
 145    TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 60, nmin, nmax);
 146    hist->GetXaxis()->SetTitle("p_{T} GeV / c");
 147    fOutputList->Add(hist);
 148    return hist;
 149 }
 150 //_____________________________________________________________________________
 151 TH1F* AliAnalysisTaskPhiFlow::BookPtHistogram(const char* name)
 152 {
 153    // Return a pointer to a p_T spectrum histogram
 154    if(fDebug > 0) cout << " *** BookPtHistogram() *** " << endl;
 155    TH1F* ratio = new TH1F(name, name, 100, 0, 7);
 156    ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
 157    ratio->GetYaxis()->SetTitle("No. of events");
 158    ratio->Sumw2();
 159    fOutputList->Add(ratio);
 160    return ratio;
 161 }
 162 //_____________________________________________________________________________
 163 void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
 164 {
 165    // Add Phi Identification Output Objects
 166    if(fDebug > 0) cout << " ** AddPhiIdentificationOutputObjects() *** " << endl;
 167    if(fQA) {
 168        fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
 169        fEventStats->GetXaxis()->SetTitle("No. of events");
 170        fEventStats->GetYaxis()->SetTitle("Statistics");
 171        fOutputList->Add(fEventStats);
 172    }
 173    fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
 174    fOutputList->Add(fCentralityPass);
 175    if(fQA) {
 176        fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
 177        fOutputList->Add(fCentralityNoPass);
 178        fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
 179        fPIDk = BookPIDHistogram("TPC signal, kaons", kTRUE);
 180        fNOPIDTOF = BookPIDHistogram("TOF signal, all particles", kFALSE);
 181        fPIDTOF = BookPIDHistogram("TOF signal, kaons", kFALSE);
 182    }
 183    for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
 184        fInvMNP[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
 185        fInvMNN[ptbin] = BookHistogram(Form("NN, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
 186        fInvMPP[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
 187    }
 188    if(fQA) {
 189        for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
 190        fPtP = BookPtHistogram("i^{+}");
 191        fPtN = BookPtHistogram("i^{-}");
 192        fPtKP = BookPtHistogram("K^{+}");
 193        fPtKN = BookPtHistogram("K^{-}");
 194        fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
 195        fOutputList->Add(fPhi);
 196        fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
 197        fOutputList->Add(fPt);
 198        fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
 199        fOutputList->Add(fEta);
 200        fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
 201        fOutputList->Add(fVZEROA);
 202        fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
 203        fOutputList->Add(fVZEROC);
 204        fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
 205        fOutputList->Add(fTPCM);
 206        fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
 207        fOutputList->Add(fDCAXYQA);
 208        fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
 209        fOutputList->Add(fDCAZQA);
 210        if(fCentralityCut2010 || fCentralityCut2011) {
 211            fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
 212            fOutputList->Add(fMultCorAfterCuts);
 213            fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
 214            fOutputList->Add(fMultvsCentr);
 215        }
 216    }
 217    if(fIsMC || fQA) {
 218        fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
 219        fOutputList->Add(fDCAAll);
 220        fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
 221        fOutputList->Add(fDCAPrim);
 222        fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
 223        fOutputList->Add(fDCASecondaryWeak);
 224        fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
 225        fOutputList->Add(fDCAMaterial);
 226    }
 227    if(fV0) {
 228        fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 5, 0, 5);
 229        fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<#Psi_{a} - #Psi_{b}>");
 230        fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<#Psi_{a} - #Psi_{c}>");
 231        fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<#Psi_{b} - #Psi_{c}>");
 232        fSubEventDPhiv2->GetXaxis()->SetBinLabel(4, "#sqrt{#frac{<#Psi_{a} - #Psi_{b}><#Psi_{a} - #Psi_{c}>}{<#Psi_{b} - #Psi_{c}>}}");
 233        fSubEventDPhiv2->GetXaxis()->SetBinLabel(5, "#sqrt{#frac{<#Psi_{a} - #Psi_{c}><#Psi_{b} - #Psi_{c}>}{<#Psi_{a} - #Psi_{b}>}}");
 234        fOutputList->Add(fSubEventDPhiv2);
 235        const char* V0[] = {"V0A", "V0C"};
 236        for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
 237            for(Int_t iV0(0); iV0 < 2; iV0++) {
 238                    fV0Data[ptbin][iV0] = new TProfile(Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), fMassBins, fMinMass, fMaxMass);
 239                    fOutputList->Add(fV0Data[ptbin][iV0]);
 240            }
 241    }
 242 }
 243 //_____________________________________________________________________________
 244 void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
 245 {
 246    // Create user defined output objects
 247    if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
 248    fNullCuts = new AliFlowTrackCuts("null_cuts");
 249    fBayesianResponse = new AliFlowBayesianPID();
 250       // combined pid
 251    fPIDCombined = new AliPIDCombined;
 252    fPIDCombined->SetDefaultTPCPriors();
 253    fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 254
 255    // flag to mc
 256    if(fSkipEventSelection || fIsMC) fBayesianResponse->SetMC(kTRUE);
 257    Double_t t(0);
 258    for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
 259    if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
 260    if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
 261    AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 262    cc->SetNbinsQ(500);           cc->SetNbinsPhi(180);           cc->SetNbinsMult(10000);
 263    cc->SetQMin(0.0);             cc->SetPhiMin(0.0);             cc->SetMultMin(0);
 264    cc->SetQMax(3.0);             cc->SetPhiMax(TMath::TwoPi());  cc->SetMultMax(10000);
 265    cc->SetNbinsMass(fMassBins);  cc->SetNbinsEta(200);           (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
 266    cc->SetMassMin(fMinMass);     cc->SetEtaMin(-5.0);            cc->SetPtMin(0);
 267    cc->SetMassMax(fMaxMass);     cc->SetEtaMax(+5.0);            (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
 268    fBayesianResponse->SetNewTrackParam();
 269    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
 270    if (man) {
 271       AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 272       if (inputHandler)   fPIDResponse = inputHandler->GetPIDResponse();
 273    }
 274    // Create all output objects and store them to a list
 275    fOutputList = new TList();
 276    fOutputList->SetOwner(kTRUE);
 277    // Create and post the Phi identification output objects
 278    AddPhiIdentificationOutputObjects();
 279    PostData(1, fOutputList);
 280    // create candidate array
 281    fCandidates = new TObjArray(1000);
 282    fCandidates->SetOwner(kTRUE);
 283    // create and post flowevent
 284    fFlowEvent = new AliFlowEvent(10000);
 285    PostData(2, fFlowEvent);
 286    if(fEventMixing) fPoolManager = InitializeEventMixing();
 287 }
 288 //_____________________________________________________________________________
 289 AliEventPoolManager* AliAnalysisTaskPhiFlow::InitializeEventMixing()
 290 {
 291    // initialize event mixing
 292   if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
 293   Int_t _c(0), _v(0);
 294   for(Int_t i(0); i < 19; i++) {
 295       if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
 296       else _c = 19;
 297   }
 298   for(Int_t i(0); i < 19; i++) {
 299       if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
 300       else _v = 19;
 301   }
 302   if(fDebug > 0 ) cout << Form("   --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
 303   Double_t centralityBins[_c];
 304   Double_t vertexBins[_v];
 305   for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
 306   for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
 307   return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
 308 }
 309 //_____________________________________________________________________________
 310 template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
 311 {
 312    // Return the invariant mass of two tracks, assuming both tracks are kaons
 313    if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
 314    if ((!track2) || (!track1)) return 0.;
 315    Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
 316    Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
 317    Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
 318    Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
 319    Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
 320    Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
 321    Double_t es = TMath::Power((e1 + e2), 2);
 322    if ((es - (pxs + pys + pzs)) < 0) return 0.;
 323    return TMath::Sqrt((es - (pxs + pys + pzs)));
 324 }
 325 //_____________________________________________________________________________
 326 /*
 327 template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
 328 {
 329    // Calculate the delta dip angle between two particles
 330    if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
 331    if (track1->P()*track2->P() == 0) return 999;
 332    return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
 333 }
 334 //_____________________________________________________________________________
 335 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T* track1, const T* track2) const
 336 {
 337    // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
 338    if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
 339    if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
 340    return kTRUE;
 341 }
 342 */
 343 //_____________________________________________________________________________
 344 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
 345 {
 346    // Check if pair passes eta and pt cut
 347    if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
 348    if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
 349    TVector3 a(track1->Px(), track1->Py(), track1->Pz());
 350    TVector3 b(track2->Px(), track2->Py(), track2->Pz());
 351    TVector3 c = a + b;
 352    if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
 353    return kTRUE;
 354 }
 355 //_____________________________________________________________________________
 356 template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
 357 {
 358    // Impose event cuts
 359    if(fDebug > 0) cout << " *** EventCut() *** " << endl;
 360    if (!event) return kFALSE;
 361    if (fSkipEventSelection) return kTRUE;
 362    if (!CheckVertex(event)) return kFALSE;
 363    if (!CheckCentrality(event)) return kFALSE;
 364    if(fQA) PlotMultiplcities(event);
 365    return kTRUE;
 366 }
 367 //_____________________________________________________________________________
 368 template <typename T> void AliAnalysisTaskPhiFlow::PlotMultiplcities(const T* event) const
 369 {
 370    // QA multiplicity plots
 371    if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
 372    fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
 373    fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
 374    fTPCM->Fill(event->GetNumberOfTracks());
 375 }
 376 //_____________________________________________________________________________
 377 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const T* event)
 378 {
 379    // Check if event vertex is within given range
 380    if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
 381    if (!event->GetPrimaryVertex()) return 0x0;
 382    fVertex = event->GetPrimaryVertex()->GetZ();
 383    if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
 384    return kTRUE;
 385 }
 386 //_____________________________________________________________________________
 387 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
 388 {
 389    // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
 390    if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
 391    if (!fkCentralityMethodA) AliFatal("No centrality method set! FATAL ERROR!");
 392    fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
 393    Double_t cenB(-999);
 394    // if a second centrality estimator is requited, set it
 395    (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
 396    if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
 397       if(fQA) fCentralityNoPass->Fill(fCentrality) ;
 398       return kFALSE;
 399    }
 400    const Int_t nGoodTracks = event->GetNumberOfTracks();
 401    if(fCentralityCut2010) { // cut on outliers
 402       Float_t multTPC(0.); // tpc mult estimate
 403       Float_t multGlob(0.); // global multiplicity
 404       for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
 405           AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
 406           if(!trackAOD) AliFatal("Not a standard AOD");
 407           if (!trackAOD) continue;
 408           if (!(trackAOD->TestFilterBit(1))) continue;
 409           if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
 410           multTPC++;
 411       }
 412       for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
 413           AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
 414           if(!trackAOD) AliFatal("Not a standard AOD");
 415           if (!trackAOD) continue;
 416           if (!(trackAOD->TestFilterBit(16))) continue;
 417           if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
 418           Double_t b[2] = {-99., -99.};
 419           Double_t bCov[3] = {-99., -99., -99.};
 420           AliAODTrack copy(*trackAOD);
 421           if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
 422           if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
 423           multGlob++;
 424       } //track loop
 425  //     printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
 426       if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
 427       if(fQA) {
 428           fMultCorAfterCuts->Fill(multGlob, multTPC);
 429           fMultvsCentr->Fill(fCentrality, multTPC);
 430       }
 431    }
 432
 433    if(fCentralityCut2011) { // cut on outliers
 434       Float_t multTPC(0.); // tpc mult estimate
 435       Float_t multGlob(0.); // global multiplicity
 436       for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
 437           AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
 438           if(!trackAOD) AliFatal("Not a standard AOD");
 439           if (!trackAOD) continue;
 440           if (!(trackAOD->TestFilterBit(1))) continue;
 441           if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
 442           multTPC++;
 443       }
 444       for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
 445           AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
 446           if(!trackAOD) AliFatal("Not a standard AOD");
 447           if (!trackAOD) continue;
 448           if (!(trackAOD->TestFilterBit(16))) continue;
 449           if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
 450           Double_t b[2] = {-99., -99.};
 451           Double_t bCov[3] = {-99., -99., -99.};
 452           AliAODTrack copy(*trackAOD);
 453           if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
 454           if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
 455           multGlob++;
 456       } //track loop
 457       //printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
 458       if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
 459       if(fQA) {
 460           fMultCorAfterCuts->Fill(multGlob, multTPC);
 461           fMultvsCentr->Fill(fCentrality, multTPC);
 462       }
 463    }
 464
 465    fCentralityPass->Fill(fCentrality);
 466    return kTRUE;
 467 }
 468 //_____________________________________________________________________________
 469 void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
 470 {
 471    // Initialize the Bayesian PID object for AOD
 472    if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
 473    fBayesianResponse->SetDetResponse(event, fCentrality);
 474 }
 475 //_____________________________________________________________________________
 476 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
 477 {
 478    // Check if the particle passes the TPC TOF bayesian cut.
 479    if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
 480    fBayesianResponse->ComputeProb(track);
 481    if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
 482    Float_t *probabilities = fBayesianResponse->GetProb();
 483    if (probabilities[3] > fPIDConfig[6]) {
 484       if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
 485       return kTRUE;
 486    }
 487    return kFALSE;
 488 }
 489 //_____________________________________________________________________________
 490 Bool_t AliAnalysisTaskPhiFlow::PassesDCACut(AliAODTrack* track) const
 491 {
 492     // check if track passes dca cut according to dca routine
 493     // setup the routine as follows:
 494     // fDCAConfig[0] < -1 no pt dependence
 495     // fDCAConfig[0] =  0 do nothing
 496     // fDCAConfig[0] >  1 pt dependent dca cut
 497     if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
 498     if(fIsMC) return kTRUE;
 499     AliAODTrack copy(*track);
 500     if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
 501         if(fQA) {
 502             Double_t b[2] = { -99., -99.};
 503             Double_t bCov[3] = { -99., -99., -99.};
 504             if(copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
 505                 fDCAXYQA->Fill(b[0]);
 506                 fDCAZQA->Fill(b[1]);
 507                 fDCAPrim->Fill(track->Pt(), b[0]);
 508             }
 509         }
 510         return kTRUE;
 511     }
 512     Double_t b[2] = { -99., -99.};
 513     Double_t bCov[3] = { -99., -99., -99.};
 514     if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
 515     if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
 516     if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
 517     if(fDCAConfig[0] > .9) {
 518         if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
 519         Double_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
 520         if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
 521         if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
 522     }
 523     if(fQA) {
 524         fDCAXYQA->Fill(b[0]);
 525         fDCAZQA->Fill(b[1]);
 526         fDCAPrim->Fill(track->Pt(), b[0]);
 527     }
 528     return kTRUE;
 529 }
 530 //_____________________________________________________________________________
 531 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
 532 {
 533    // Kaon identification routine, based on multiple detectors and approaches
 534    if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
 535    if(fUsePidResponse) {
 536        Double_t prob[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 537        fPIDCombined->ComputeProbabilities(track, fPIDResponse, prob);
 538        if(prob[3] > fPIDConfig[6])  return kTRUE;
 539    }
 540    if(!PassesDCACut(track)) return kFALSE;
 541    if(fQA) {fNOPID->Fill(track->P(), track->GetTPCsignal());fNOPIDTOF->Fill(track->P(), track->GetTOFsignal());}
 542    if(track->Pt() < fPIDConfig[1]) {
 543        if(fDebug>1) cout << " ITS received track with p_t " << track->Pt() << endl;
 544        // if tpc control is disabled, pure its pid
 545        if(fPIDConfig[2] < 0.) {
 546            if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
 547            return kFALSE;
 548        }
 549        // else, switch to ITS pid with TPC rejection of protons and pions
 550        if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
 551        else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
 552        else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) {
 553            if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
 554            return kTRUE;
 555        }
 556        return kFALSE;
 557    }
 558    if((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4])) {
 559        if(fDebug>1) cout << " TPC received track with p_t " << track->Pt() << endl;
 560        // if its control is disabled, pure tpc pid
 561        if(fPIDConfig[5] < 0.) {
 562            if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
 563            return kFALSE;
 564        }
 565        // else, switch to TPC pid with ITS rejection of protons and pions
 566        if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
 567        else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
 568        else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) {
 569            if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
 570            return kTRUE;
 571        }
 572        return kFALSE;
 573    }
 574    if(fDebug>1) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
 575    // switch to bayesian PID
 576    if (PassesTPCbayesianCut(track)) {
 577        if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal());fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
 578        return kTRUE;
 579    }
 580    return kFALSE;
 581 }
 582 //_____________________________________________________________________________
 583 template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
 584 {
 585    // return p_t of track pair
 586    TVector3 a(track1->Px(), track1->Py(), track1->Pz());
 587    TVector3 b(track2->Px(), track2->Py(), track2->Pz());
 588    TVector3 c = a + b;
 589    return c.Pt();
 590 }
 591 //_____________________________________________________________________________
 592 template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
 593 {
 594    // plot m_inv spectra of like- and unlike-sign kaon pairs for each pt bin
 595    Double_t pt = PhiPt(track1, track2);
 596    if (tracktype == 0) {
 597        for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
 598            if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
 599                fInvMNP[ptbin]->Fill(InvariantMass(track1, track2));
 600                if(fQA) fPtSpectra[ptbin]->Fill(pt);
 601            }
 602        }
 603    }
 604    if (tracktype == 1) {
 605        for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
 606            if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
 607                fInvMPP[ptbin]->Fill(InvariantMass(track1, track2));
 608            }
 609        }
 610    }
 611    if (tracktype == 2) {
 612        for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
 613            if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
 614                fInvMNN[ptbin]->Fill(InvariantMass(track1, track2));
 615            }
 616        }
 617    }
 618 }
 619 //_____________________________________________________________________________
 620 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
 621 {
 622    // check if track meets quality cuts
 623    if(!track) return kFALSE;
 624    return fPOICuts->IsSelected(track);
 625 }
 626 //_____________________________________________________________________________
 627 template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
 628 {
 629    // Set null cuts
 630    if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
 631    fCutsRP->SetEvent(event, MCEvent());
 632    fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
 633    fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
 634    fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
 635    fNullCuts->SetEvent(event, MCEvent());
 636 }
 637 //_____________________________________________________________________________
 638 void AliAnalysisTaskPhiFlow::PrepareFlowEvent(Int_t iMulti)
 639 {
 640    // Prepare flow events
 641    if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
 642    fFlowEvent->ClearFast();
 643    fFlowEvent->Fill(fCutsRP, fNullCuts);
 644    fFlowEvent->SetReferenceMultiplicity(iMulti);
 645    fFlowEvent->DefineDeadZone(0, 0, 0, 0);
 646 }
 647 //_____________________________________________________________________________
 648 void AliAnalysisTaskPhiFlow::VZEROSubEventAnalysis()
 649 {
 650     // vzero event plane analysis using three subevents
 651     if(fDebug > 0) cout << " ** VZEROSubEventAnalysis() *** " << endl;
 652     if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
 653         if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
 654         return;
 655     }
 656     // retrieve data
 657     Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
 658     for(Int_t i(0); i < 3; i++) if(abcPsi2[i] == 0)  {
 659         if(fDebug > 0) cout << " Warning: I've encountered 0 in one of the symmetry panes (TPC, VZERO) - skipping event !" << endl;
 660             return;
 661     }
 662     // save info for resolution calculations
 663     Float_t qaqb = TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1])); // vzeroa - tpc
 664     Float_t qaqc = TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2])); // vzeroa - vzeroc
 665     Float_t qbqc = TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2])); // tpc - vzeroc
 666     fSubEventDPhiv2->Fill(0.5, qaqb);
 667     fSubEventDPhiv2->Fill(1.5, qaqc);
 668     fSubEventDPhiv2->Fill(2.5, qbqc);
 669     Float_t minv[31];
 670     Float_t _dphi[30][fNPtBins][2]; //minv, pt, v0a-c
 671     Int_t occurence[30][fNPtBins]; //minv, pt
 672     for(Int_t mb(0); mb < 31; mb++) minv[mb] = 0.99 + mb * 0.0034;
 673     for(Int_t i(0); i < 30; i++) for (Int_t j(0); j < fNPtBins; j++) for(Int_t k(0); k < 2; k ++) {
 674         _dphi[i][j][k] = 0;
 675         if(k==0) occurence[i][j] = 0;
 676     }
 677     for(Int_t iCand(0); iCand < fCandidates->GetEntriesFast(); iCand++) { // loop over unlike sign kaon pairs
 678         AliFlowTrackSimple *track = dynamic_cast<AliFlowTrackSimple*>(fCandidates->At(iCand));
 679         if(!track) {
 680             if(fDebug > 1) cout << " --> dynamic_cast returned null-pointer, skipping candidate <-- " << endl;
 681             continue;
 682         }
 683         for(Int_t mb(0); mb < 30; mb++) { // loop over massbands
 684             if(track->Mass() <= minv[mb] || track->Mass() > minv[mb+1]) continue;
 685             for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
 686                 if(track->Pt() <= fPtBins[ptbin] || track->Pt() > fPtBins[ptbin+1]) continue;
 687                 _dphi[mb][ptbin][0]+=TMath::Cos(2.*(track->Phi() - abcPsi2[0]));
 688                 _dphi[mb][ptbin][1]+=TMath::Cos(2.*(track->Phi() - abcPsi2[2]));
 689                 occurence[mb][ptbin]+=1;
 690             }
 691         }
 692         for(Int_t mb(0); mb < 30; mb++) // write vn values to tprofiles
 693             for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
 694                 if(occurence[mb][ptbin]==0) continue;
 695                 fV0Data[ptbin][0]->Fill(mb*0.0034+0.99+0.0017, _dphi[mb][ptbin][0]/(Float_t)occurence[mb][ptbin]);
 696                 fV0Data[ptbin][1]->Fill(mb*0.0034+0.99+0.0017, _dphi[mb][ptbin][1]/(Float_t)occurence[mb][ptbin]);
 697             }
 698     }
 699 }
 700 //_____________________________________________________________________________
 701 void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
 702 {
 703    // UserExec: called for each event. Commented where necessary
 704    if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
 705    TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
 706    if(fEventMixing) {
 707        MixingCandidates = new TObjArray();
 708        MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
 709    }
 710    if (!fPIDResponse) {
 711       if(fDebug > 0 ) cout << " Could not get PID response " << endl;
 712       return;
 713    }
 714    fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
 715    if (fAOD) {
 716       if (!EventCut(fAOD)) return; // check for event cuts
 717       InitializeBayesianPID(fAOD); // init event objects
 718       SetNullCuts(fAOD);
 719       PrepareFlowEvent(fAOD->GetNumberOfTracks());
 720       fCandidates->SetLast(-1);
 721       if(fIsMC) IsMC(); // launch mc stuff
 722       if(fQA) fEventStats->Fill(0);
 723       Int_t unTracks = fAOD->GetNumberOfTracks();
 724       AliAODTrack* un[unTracks];
 725       AliAODTrack* up[unTracks];
 726       Int_t unp(0);
 727       Int_t unn(0);
 728       for (Int_t iTracks = 0; iTracks < unTracks; iTracks++) { // select analysis candidates
 729          AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 730          if(!track) AliFatal("Not a standard AOD");
 731          if (!PhiTrack(track)) continue;
 732          if (fQA) {
 733             if(track->Charge() > 0) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
 734             if(track->Charge() < 0) {fEventStats->Fill(2); fPtN->Fill(track->Pt());}
 735          }
 736          if (IsKaon(track)) {
 737             if(fEventMixing) MixingCandidates->Add(new AliPhiMesonHelperTrack(track->Eta(), track->Phi(), track->P(),
 738                                                                               track->Px(), track->Py(), track->Pz(),
 739                                                                               track->Pt(), track->Charge()));
 740             if (track->Charge() > 0) {
 741                up[unp] = track;
 742                unp++;
 743                if(fQA) {fEventStats->Fill(3);fPtKP->Fill(track->Pt());}
 744             }
 745             else if (track->Charge() < 0) {
 746                un[unn] = track;
 747                unn++;
 748                if(fQA) {fEventStats->Fill(4); fPtKN->Fill(track->Pt());}
 749             }
 750          }
 751       }
 752       for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
 753          for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
 754 //            if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
 755             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
 756             PtSelector(0, up[pTracks], un[nTracks]);
 757             Double_t pt = PhiPt(up[pTracks], un[nTracks]);
 758             Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
 759             TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
 760             TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
 761             TVector3 c = a + b;
 762             Double_t phi = c.Phi();
 763             Double_t eta = c.Eta();
 764             Double_t p = TMath::Sqrt(c.Px()*c.Px()+c.Py()*c.Py()+c.Pz()*c.Pz());
 765             Int_t nIDs[2];
 766             nIDs[0] = up[pTracks]->GetID();
 767             nIDs[1] = un[nTracks]->GetID();
 768             MakeTrack(mass, pt, phi, eta, 2, nIDs, p, c.Pz());
 769          }
 770       }
 771       if(fV0) VZEROSubEventAnalysis();
 772       if (fDebug > 0)  printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
 773       for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
 774          AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
 775          if (!cand) continue;
 776          if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
 777          for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
 778             if (fDebug>1) printf("     *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
 779             for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
 780                AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
 781                if (!iRP) continue;
 782                if (!iRP->InRPSelection()) continue;
 783                if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
 784                   if (fDebug > 1) printf("      was in RP set");
 785                   iRP->SetForRPSelection(kFALSE);
 786                   fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
 787                }
 788             }
 789             if (fDebug > 1) printf("\n");
 790          }
 791          cand->SetForPOISelection(kTRUE);
 792          fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
 793       }
 794       if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
 795       if(!fEventMixing) { // combinatorial background
 796           for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
 797              for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
 798 //                if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
 799                 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
 800                 PtSelector(1, up[pTracks], up[nTracks]);
 801              }
 802           }
 803           for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
 804              for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
 805 //                if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
 806                 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
 807                 PtSelector(2, un[nTracks], un[pTracks]);
 808              }
 809           }
 810       }
 811       if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
 812       PostData(1, fOutputList);
 813       PostData(2, fFlowEvent);
 814    }
 815 }
 816 //_____________________________________________________________________________
 817 void AliAnalysisTaskPhiFlow::Exec(Option_t* c)
 818 {
 819     // skip the event selection for SE task (e.g. for MC productions)
 820     if(fSkipEventSelection) AliAnalysisTaskPhiFlow::UserExec(c);
 821     // exec of task se will do event selection and call UserExec
 822     else AliAnalysisTaskSE::Exec(c);
 823 }
 824 //_____________________________________________________________________________
 825 void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
 826 {
 827     // perform phi reconstruction with event mixing
 828     if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
 829     AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
 830     if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
 831     if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
 832         Int_t nEvents = pool->GetCurrentNEvents();
 833         if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
 834         for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
 835             TObjArray* mixed_candidates = pool->GetEvent(iEvent);
 836             if(!mixed_candidates) continue; // this should NEVER happen
 837             Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
 838             Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
 839             if(fDebug > 0) cout << Form("   - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
 840             AliPhiMesonHelperTrack* buffer_un[bufferTracks]; // set values for buffered candidates
 841             AliPhiMesonHelperTrack* buffer_up[bufferTracks];
 842             Int_t buffer_unp(0);
 843             Int_t buffer_unn(0);
 844             AliPhiMesonHelperTrack* mix_un[candidates];// set values for mixing candidates
 845             AliPhiMesonHelperTrack* mix_up[candidates];
 846             Int_t mix_unp(0);
 847             Int_t mix_unn(0);
 848             for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { // distinguish between k+ and k- for mixing candidates
 849                 AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)MixingCandidates->At(iTracks);
 850                 if(!track) continue;
 851                 if (track->Charge() > 0) {
 852                     mix_up[mix_unp] = track;
 853                     mix_unp++;
 854                 }
 855                 else if (track->Charge() < 0 ) {
 856                    mix_un[mix_unn] = track;
 857                    mix_unn++;
 858                 }
 859             }
 860             for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { // distinguish between k+ and k- for buffered candidates
 861                 AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)mixed_candidates->At(iTracks);
 862                 if(!track) continue;
 863                 if (track->Charge() > 0) {
 864                     buffer_up[buffer_unp] = track;
 865                     buffer_unp++;
 866                 }
 867                 else if (track->Charge() < 0 ) {
 868                    buffer_un[buffer_unn] = track;
 869                    buffer_unn++;
 870                 }
 871             }
 872             for (Int_t pMix = 0; pMix < mix_unp; pMix++) { // mix k+ (current event) k+ (buffer)
 873                 if(fDebug > 1 ) cout << Form("mixing current k+ track %d with", pMix);
 874                 if(!fTypeMixing) { // mix with like-sign kaons
 875                     for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
 876                         if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
 877 //                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
 878                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
 879                         PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
 880                     }
 881                 }
 882                 else if(fTypeMixing) { // mix with unlike-sign kaons
 883                     for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
 884                         if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
 885 //                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
 886                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
 887                         PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
 888                     }
 889                 }
 890             }
 891             for (Int_t nMix = 0; nMix < mix_unn; nMix++) { // mix k- (current event) k- (buffer)
 892                 if(fDebug > 1 ) cout << Form("mixing current k- track %d with", nMix);
 893                 if(!fTypeMixing) { // mix with like-sign kaons
 894                     for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
 895                         if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
 896 //                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
 897                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
 898                         PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
 899                     }
 900                 }
 901                 else if(fTypeMixing) { // mix with unlike-sign kaons
 902                     for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
 903                         if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
 904 //                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
 905                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
 906                         PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
 907                     }
 908                 }
 909             }
 910         } // end of mixed events loop
 911     } // end of checking to see whether pool is filled correctly
 912     pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
 913     if(fDebug > 0 ) pool->PrintInfo();
 914 }
 915 //_____________________________________________________________________________
 916 void AliAnalysisTaskPhiFlow::Terminate(Option_t *)
 917 {
 918     // terminate
 919     if(fDebug > 0) cout << " *** Terminate() *** " << endl;
 920 }
 921 //______________________________________________________________________________
 922 void  AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass, Double_t pt, Double_t phi, Double_t eta, Int_t nDau, Int_t iID[], Double_t p, Double_t pz) const
 923 {
 924    // Construct Flow Candidate Track from two selected candidates
 925    if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
 926    // if requested, check rapidity at this point
 927    if(fCandidateYCut) {
 928        Double_t y = 0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
 929        if (y > fCandidateMaxY || y < fCandidateMinY) return;
 930    }
 931    Bool_t overwrite = kTRUE;
 932    AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
 933    if (!sTrack) {
 934       sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
 935       overwrite = kFALSE;
 936    }
 937    else sTrack->ClearMe();
 938    sTrack->SetMass(mass);
 939    sTrack->SetPt(pt);
 940    sTrack->SetPhi(phi);
 941    sTrack->SetEta(eta);
 942    for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
 943    sTrack->SetForPOISelection(kTRUE);
 944    sTrack->SetForRPSelection(kFALSE);
 945    if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
 946    else fCandidates->AddLast(sTrack);
 947    return;
 948 }
 949 //_____________________________________________________________________________
 950 void AliAnalysisTaskPhiFlow::IsMC()
 951 {
 952     // Fill QA histos for MC analysis
 953    TClonesArray *arrayMC = 0;
 954    if(fDebug > 0) cout << " -> Switching to MC mode <- " << endl;
 955    // fill array with mc tracks
 956    arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 957    if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
 958    for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
 959      AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
 960      if(!track) AliFatal("Not a standard AOD");
 961      if(!PhiTrack(track) || !IsKaon(track)) { // check for kaons
 962          if(fDebug>1) cout << " Rejected track" << endl;
 963          continue;
 964      }
 965      if (fDebug>1) cout << " Received MC kaon " << endl;
 966      Double_t b[2] = { -99., -99.};
 967      Double_t bCov[3] = { -99., -99., -99.};
 968      AliAODTrack copy(*track);
 969      if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
 970      // find corresponding mc particle
 971      AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
 972      if (!partMC) {
 973          if(fDebug > 1) cout << "Cannot get MC particle" << endl;
 974          continue;
 975      }
 976      // Check if it is primary, secondary from material or secondary from weak decay
 977      Bool_t isPrimary           = partMC->IsPhysicalPrimary();
 978      Bool_t isSecondaryMaterial = kFALSE;
 979      Bool_t isSecondaryWeak     = kFALSE;
 980      if (!isPrimary) {
 981          Int_t mfl = -999, codemoth = -999;
 982          Int_t indexMoth = partMC->GetMother();
 983          if (indexMoth >= 0) { // is not fake
 984             AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
 985             codemoth = TMath::Abs(moth->GetPdgCode());
 986             mfl = Int_t(codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
 987          }
 988          if (mfl == 3) isSecondaryWeak     = kTRUE;
 989          else       isSecondaryMaterial = kTRUE;
 990       }
 991       if (isPrimary) {
 992           fDCAPrim->Fill(track->Pt(), b[0]);
 993           fDCAXYQA->Fill(b[0]);
 994           fDCAZQA->Fill(b[1]);
 995       }
 996       if (isSecondaryWeak)  fDCASecondaryWeak->Fill(track->Pt(), b[0]);
 997       if (isSecondaryMaterial) fDCAMaterial->Fill(track->Pt(), b[0]);
 998    }
 999 }
1000 //_____________________________________________________________________________