PWGCF/FLOW/Attic/AliAnalysisTaskJetFlow.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 /*
  17  * AliAnaysisTaskJetFlow
  18  * author: Redmer Alexander Bertens
  19  * rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
  20  *
  21  * Interface task between EMCal jet framework and flow package
  22  *
  23  * This task expects POI's in a TClonesArray  (e.g. from running it after
  24  * $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation() )
  25  * and connects them into an AliFlowEvent which is filled with either VZERO tracks or
  26  * TPC trakcs.
  27  *
  28  * POI's can be supplied as AliEmcalJets or as AliVTracks, note the different behavior
  29  * with respect to the Pt value: for AliEmcalJets subtracted Pt is used by default, for VParticles
  30  * Pt is taken directly.
  31  *
  32  * */
  33
  34 // root includes
  35 #include <TChain.h>
  36 #include <TH1F.h>
  37 #include <TArrayD.h>
  38 #include <TProfile.h>
  39 #include <TString.h>
  40 #include <TList.h>
  41 #include <TClonesArray.h>
  42 // aliroot includes
  43 #include <AliAODEvent.h>
  44 #include <AliESDEvent.h>
  45 #include <AliEmcalJet.h>
  46 #include <AliAnalysisTask.h>
  47 #include <AliAnalysisManager.h>
  48 #include <AliFlowEventSimple.h>
  49 #include <AliFlowEvent.h>
  50 #include <AliFlowTrack.h>
  51 #include <AliFlowTrackCuts.h>
  52 #include <AliFlowEventCuts.h>
  53 #include <AliFlowCommonConstants.h>
  54 // local includes
  55 #include "AliAnalysisTaskJetFlow.h"
  56 // EMCAL jet framework includes
  57 #include <AliRhoParameter.h>
  58 #include <AliLocalRhoParameter.h>
  59 #include <AliPicoTrack.h>
  60 #include <AliAnalysisTaskRhoVnModulation.h>
  61
  62 class AliAnalysisTaskJetFlow;
  63
  64 using namespace std;
  65
  66 ClassImp(AliAnalysisTaskJetFlow)
  67
  68 AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow() : AliAnalysisTaskSE(),
  69     fDebug(-1), fJetsName(0), fJetRadius(0.3), fTracksName(0), fLocalRhoName(0), fPois(0x0), fRPs(0x0), fLocalRho(0x0), fOutputList(0), fDataType(kESD), fVParticleAnalysis(kFALSE), fMinimizeDiffBins(kTRUE), fDoVZEROFlowAnalysis(kTRUE), fDoGappedQC2Analysis(kTRUE), fDoQC2FlowAnalysis(kTRUE), fDoQC4FlowAnalysis(kFALSE), fDoQCFPAnalysis(kFALSE), fDoSPFPAnalysis(kFALSE), fDoMultWeight(kTRUE), fDoPtWeight(0), fInitialized(kFALSE), fUsePtWeight(kFALSE), fCCMinPt(1), fCCMaxPt(150), fCCBinsInPt(50), fCentralityMin(-1), fCentralityMax(-1), fPtBins(0), fCutsRP_VZERO(0), fCutsNull(0), fCutsEvent(0), fFlowEvent_TPC(0), fFlowEvent_VZERO(0), fRhoVn(0), fHistAnalysisSummary(0), fCentralitySelection(0), fVZEROAEP(0), fVZEROCEP(0), fv2VZEROA(0), fv2VZEROC(0), fRefCumulants(0), fDiffCumlantsV2(0), fDiffCumlantsV3(0), fQC2v2(0), fQC2v3(0), fTempA(0), fTempC(0)
  70 { /* default constructor for ROOT IO */ }
  71 //_____________________________________________________________________________
  72 AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow(
  73         const char* name,
  74         AliAnalysisTaskRhoVnModulation* rhoTask,
  75         Bool_t VPart,
  76         Bool_t VZEROEP,
  77         Bool_t GQC2,
  78         Bool_t QC2,
  79         Bool_t QC4,
  80         Bool_t FlowPackageSP,
  81         Bool_t FlowPackageQC
  82         ) : AliAnalysisTaskSE(name),
  83     fDebug(-1), fJetsName(0), fJetRadius(0.3), fTracksName(0), fLocalRhoName(0), fPois(0x0), fRPs(0x0), fLocalRho(0x0), fOutputList(0), fDataType(kESD), fVParticleAnalysis(VPart), fMinimizeDiffBins(kTRUE), fDoVZEROFlowAnalysis(VZEROEP), fDoGappedQC2Analysis(GQC2), fDoQC2FlowAnalysis(QC2), fDoQC4FlowAnalysis(QC4), fDoQCFPAnalysis(FlowPackageQC), fDoSPFPAnalysis(FlowPackageSP), fDoMultWeight(kTRUE), fDoPtWeight(0), fInitialized(kFALSE), fUsePtWeight(kFALSE), fCCMinPt(1), fCCMaxPt(150), fCCBinsInPt(50), fCentralityMin(-1), fCentralityMax(-1), fPtBins(0), fCutsRP_VZERO(0x0), fCutsNull(0), fCutsEvent(0), fFlowEvent_TPC(0), fFlowEvent_VZERO(0), fRhoVn(rhoTask), fHistAnalysisSummary(0), fCentralitySelection(0), fVZEROAEP(0), fVZEROCEP(0), fv2VZEROA(0), fv2VZEROC(0), fRefCumulants(0), fDiffCumlantsV2(0), fDiffCumlantsV3(0), fQC2v2(0), fQC2v3(0), fTempA(0), fTempC(0)
  84 {
  85     // constructor
  86     DefineInput(0, TChain::Class());
  87     DefineOutput(1, TList::Class());
  88     fJetsName = rhoTask->GetJetsName();
  89     fTracksName = rhoTask->GetTracksName();
  90     if(GQC2 && QC2) {
  91         printf(" > Warning, QC2 and gapped QC2 method are both called <\n   will only run gapped QC2 !");
  92         fDoQC2FlowAnalysis = kFALSE;
  93     }
  94     if(FlowPackageSP || FlowPackageQC)    DefineOutput(2, AliFlowEventSimple::Class());
  95     if(FlowPackageSP && FlowPackageQC)    DefineOutput(3, AliFlowEventSimple::Class());
  96 }
  97 //_____________________________________________________________________________
  98 AliAnalysisTaskJetFlow::~AliAnalysisTaskJetFlow()
  99 {
 100     // destructor
 101     if(fOutputList)             delete fOutputList;
 102     if(fFlowEvent_TPC)          delete fFlowEvent_TPC;
 103     if(fFlowEvent_VZERO)        delete fFlowEvent_VZERO;
 104     if(fCutsEvent)              delete fCutsEvent;
 105     if(fCutsRP_VZERO)           delete fCutsRP_VZERO;
 106     if(fCutsNull)               delete fCutsNull;
 107     if(fPtBins)                 delete fPtBins;
 108     if(fTempA)                  delete fTempA;
 109     if(fTempC)                  delete fTempC;
 110 }
 111 //_____________________________________________________________________________
 112 void AliAnalysisTaskJetFlow::UserCreateOutputObjects()
 113 {
 114     // create output objects
 115     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 116     // create the cut objects
 117     if(fDoSPFPAnalysis || fDoQCFPAnalysis) {
 118         fCutsEvent = new AliFlowEventCuts();
 119         fCutsEvent->SetRefMultMethod(AliESDtrackCuts::kTrackletsITSTPC);
 120         fCutsNull = new AliFlowTrackCuts("CutsNull");
 121         fCutsNull->SetParamType(AliFlowTrackCuts::kGlobal);
 122         fCutsNull->SetEtaRange(+1, -1);
 123         fCutsNull->SetPtRange(+1, -1);
 124         fCutsRP_VZERO = new AliFlowTrackCuts();
 125         fCutsRP_VZERO = fCutsRP_VZERO->GetStandardVZEROOnlyTrackCuts();
 126         if(fDoSPFPAnalysis) {
 127             (fVParticleAnalysis) ? fFlowEvent_VZERO = new AliFlowEvent(10000) : fFlowEvent_VZERO = new AliFlowEvent(100);
 128         }
 129         if(fDoQCFPAnalysis) fFlowEvent_TPC = new AliFlowEvent(10000);
 130     }
 131     fOutputList = new TList();
 132     fOutputList->SetOwner(kTRUE);
 133     // histograms
 134     fHistAnalysisSummary = new TH1F("fHistAnalysisSummary", "fHistAnalysisSummary", 4, -.5, 3.5);
 135     fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fDataType");
 136     fHistAnalysisSummary->SetBinContent(1, (int)fDataType);
 137     fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fCentralityMin");
 138     fHistAnalysisSummary->SetBinContent(2, fCentralityMin);
 139     fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fCentralityMax");
 140     fHistAnalysisSummary->SetBinContent(3, fCentralityMax);
 141     fOutputList->Add(fHistAnalysisSummary);
 142     if(fVParticleAnalysis && ! fPtBins) { // FIXME inherits from flow package now
 143         Double_t pt[fCCBinsInPt+1];
 144         Double_t width = (fCCMaxPt - fCCMinPt ) / (float)fCCBinsInPt;
 145         for(Int_t i(0); i < fCCBinsInPt+1; i++) pt[i] = fCCMinPt+width*i;
 146         fPtBins = new TArrayD(fCCBinsInPt+1, pt);  // assume they will be charged particles
 147     } else if (!fPtBins) {
 148         Double_t pt[] = {0., 10., 20., 30., 40., 50., 80., 110., 140., 170., 200.};
 149         fPtBins = new TArrayD(sizeof(pt)/sizeof(pt[0]), pt);     // assuming jets
 150     }
 151     if(fDoVZEROFlowAnalysis) {
 152         fv2VZEROA = new TProfile("Differential v_{2}^{obs} VZEROA", "Differential v_{2}^{obs} VZEROA", fPtBins->GetSize()-1, fPtBins->GetArray());
 153         fv2VZEROC = new TProfile("Differential v_{2}^{obs} VZEROC", "Differential v_{2}^{obs} VZEROC", fPtBins->GetSize()-1, fPtBins->GetArray());
 154         fv2VZEROA->GetXaxis()->SetTitle("Pt [GeV/c]");
 155         fv2VZEROA->GetYaxis()->SetTitle("v_{2}^{obs}");
 156         fv2VZEROC->GetXaxis()->SetTitle("Pt [GeV/c]");
 157         fv2VZEROC->GetYaxis()->SetTitle("v_{2}^{obs}");
 158         fOutputList->Add(fv2VZEROA);
 159         fOutputList->Add(fv2VZEROC);
 160         fVZEROAEP = new TH1F("V0A_EP", "V0A_EP", 100, -TMath::Pi()/2., TMath::Pi()/2.);
 161         fVZEROCEP = new TH1F("V0C_EP", "V0C_EP", 100, -TMath::Pi()/2., TMath::Pi()/2.);
 162         fOutputList->Add(fVZEROAEP);
 163         fOutputList->Add(fVZEROCEP);
 164         // bookkeeping histo's, will not be stored
 165         fTempA = (TProfile*)fv2VZEROA->Clone("temp_a");
 166         fTempC = (TProfile*)fv2VZEROC->Clone("temp_c");
 167     }
 168     if(fDoQC2FlowAnalysis || fDoGappedQC2Analysis) {
 169         fRefCumulants = new TProfile("Reference cumulants", "Reference cumulants", 2, -0.5, 1.5);
 170         fRefCumulants->GetXaxis()->SetBinLabel(1, "c_{2}[2]");
 171         fRefCumulants->GetXaxis()->SetBinLabel(2, "c_{3}[2]");
 172         fOutputList->Add(fRefCumulants);
 173         fDiffCumlantsV2 = new TProfile("Differential cumulants v2", "Differential cumulants v2", fPtBins->GetSize()-1, fPtBins->GetArray());
 174         fDiffCumlantsV2->GetXaxis()->SetTitle("Pt [GeV/c]");
 175         fDiffCumlantsV2->GetYaxis()->SetTitle("v_{2}[2]");
 176         fOutputList->Add(fDiffCumlantsV2);
 177         fDiffCumlantsV3 = new TProfile("Differential cumulants v3", "Differential cumulants v3", fPtBins->GetSize()-1, fPtBins->GetArray());
 178         fDiffCumlantsV3->GetXaxis()->SetTitle("Pt [GeV/c]");
 179         fDiffCumlantsV3->GetYaxis()->SetTitle("v_{3}[2]");
 180         fOutputList->Add(fDiffCumlantsV3);
 181         fQC2v2 = new TH1F("Differential v_{2}[2]", "Differential v_{2}[2]", fPtBins->GetSize()-1, fPtBins->GetArray());
 182         fQC2v2->Sumw2();
 183         fQC2v3 = new TH1F("Differential v_{3}[2]", "Differential v_{3}[2]", fPtBins->GetSize()-1, fPtBins->GetArray());
 184         fQC2v3->Sumw2();
 185         fOutputList->Add(fQC2v2);
 186         fOutputList->Add(fQC2v3);
 187     }
 188     // qa
 189     fCentralitySelection = new TH1F("fCentralitySelection", "fCentralitySelection", 100, 0, 100);
 190     fOutputList->Add(fCentralitySelection);
 191     PostData(1, fOutputList);
 192     // create the flow event and configure the static cc object
 193     Bool_t slotTwoFilled(kFALSE);
 194     if(fFlowEvent_VZERO) {
 195         PostData(2, fFlowEvent_VZERO);
 196         slotTwoFilled = kTRUE;
 197     }
 198     if(fFlowEvent_TPC) {
 199         (slotTwoFilled) ? PostData(3, fFlowEvent_TPC) : PostData(2, fFlowEvent_TPC);
 200     }
 201     if(fFlowEvent_VZERO || fFlowEvent_TPC) {
 202         AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 203         cc->SetPtMin(fCCMinPt);
 204         cc->SetPtMax(fCCMaxPt);
 205         cc->SetNbinsPt(fCCBinsInPt);
 206         if(fMinimizeDiffBins) { // minimize differential bins to reduce the risk of numerical instability
 207             cc->SetNbinsMult(1);    // only reduces output size
 208             cc->SetNbinsPhi(1);     // only reduces output size
 209             cc->SetNbinsEta(1);     // reduces instability
 210             cc->SetNbinsQ(1);       // only reduces output size
 211             cc->SetNbinsMass(1);    // reduces instability but should be one in any case ...
 212         }
 213     }
 214 }
 215 //_____________________________________________________________________________
 216 void AliAnalysisTaskJetFlow::UserExec(Option_t *)
 217 {
 218     // user exec
 219     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 220     if(!InputEvent()) return; // coverity (and sanity)
 221     // see if the analysis is initialized
 222     if(!fInitialized) {
 223         if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
 224         else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
 225         (fVParticleAnalysis) ? fPois = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName.Data())) : fPois = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName.Data()));
 226         fRPs = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName.Data()));
 227         if(!fPois || !fRPs) return; // couldn't get expected input data
 228         fInitialized = kTRUE;
 229         fLocalRho = InputEvent()->FindListObject(fLocalRhoName.Data());
 230         if(!fLocalRho && !fVParticleAnalysis) {
 231             AliFatal(Form("Couldn't find %s, aborting!", fLocalRhoName.Data()));
 232         }
 233     }
 234     if(!PassesCuts()) return; // event quality cuts and centrality determination
 235     // execute the requested flow methods
 236     if(fDoVZEROFlowAnalysis)            DoVZEROFlowAnalysis();
 237     if(fDoGappedQC2Analysis)            DoGappedQC2Analysis();
 238     if(fDoQC2FlowAnalysis)              DoQC2FlowAnalysis();
 239     if(fDoQC4FlowAnalysis)              DoQC4FlowAnalysis();
 240     Bool_t post(0x0);   // post only when analysis succeeded
 241     if(fFlowEvent_TPC || fFlowEvent_VZERO) DoFlowPackageFlowAnalysis();
 242     // push the output for the different analyses to file
 243     PostData(1, fOutputList);
 244     if(fFlowEvent_VZERO) {
 245         PostData(2, fFlowEvent_VZERO);
 246         post = kTRUE;
 247     }
 248     if(fFlowEvent_TPC) {
 249         (post) ? PostData(3, fFlowEvent_TPC) : PostData(2, fFlowEvent_TPC);
 250     }
 251 }
 252 //_____________________________________________________________________________
 253 Bool_t AliAnalysisTaskJetFlow::PassesCuts()
 254 {
 255     // passes event cuts
 256     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 257     if(!fRhoVn->PassesCuts(InputEvent())) return kFALSE;          // event quality cuts
 258     if(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M") < fCentralityMin || InputEvent()->GetCentrality()->GetCentralityPercentile("V0M") > fCentralityMax) return kFALSE;  // cutting on centrality quality is done in event quality cuts
 259     fCentralitySelection->Fill(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 260     return kTRUE;
 261 }
 262 //_____________________________________________________________________________
 263 void AliAnalysisTaskJetFlow::DoVZEROFlowAnalysis()
 264 {
 265     // flow with the VZERO event plane method
 266     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 267     Double_t vzero[2][2];
 268     fRhoVn->CalculateEventPlaneVZERO(vzero);
 269     Double_t Q2a(vzero[0][0]), Q2c(vzero[1][0]);        // just for readability
 270     fTempA->Reset();            fTempC->Reset();        // clear the containers for a new iteration
 271     fVZEROAEP->Fill(Q2a);       fVZEROCEP->Fill(Q2c);
 272     Int_t iPois(fPois->GetEntriesFast());
 273     if(fVParticleAnalysis) {
 274         for(Int_t i(0); i < iPois; i++) {
 275             AliVTrack* poi = static_cast<AliVTrack*>(fPois->At(i));
 276             if(fRhoVn->PassesCuts(poi)) {
 277                 if(!fDoMultWeight) {
 278                     fTempA->Fill(poi->Pt(), TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2a), 2)));
 279                     fTempC->Fill(poi->Pt(), TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2c), 2)));
 280                 } else {
 281                     fv2VZEROA->Fill(poi->Pt(), TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2a), 2)));
 282                     fv2VZEROC->Fill(poi->Pt(), TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2c), 2)));
 283                 }
 284             }
 285         }
 286     } else {
 287         for(Int_t i(0); i < iPois; i++) {
 288             AliEmcalJet* poi = static_cast<AliEmcalJet*>(fPois->At(i));
 289             if(fRhoVn->PassesCuts(poi)) {
 290                 AliLocalRhoParameter* localRho = static_cast<AliLocalRhoParameter*>(fLocalRho);
 291                 if(!localRho) break;
 292                 Double_t rho(localRho->GetLocalVal(poi->Phi(), fJetRadius, localRho->GetVal()));
 293                 Double_t pt(poi->Pt() - poi->Area() * rho);
 294                 if(!fDoMultWeight) {
 295                     fTempA->Fill(pt, TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2a), 2)));
 296                     fTempC->Fill(pt, TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2c), 2)));
 297                 } else {
 298                     fv2VZEROA->Fill(pt, TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2a), 2)));
 299                     fv2VZEROC->Fill(pt, TMath::Cos(2.*fRhoVn->PhaseShift((poi->Phi()-Q2c), 2)));
 300                 }
 301             }
 302         }
 303     }
 304     if(!fDoMultWeight) {
 305         for(Int_t i(0); i < fv2VZEROA->GetXaxis()->GetNbins(); i++) {
 306             fv2VZEROA->Fill(fPtBins->At(i)+(fPtBins->At(i)+fPtBins->At(1+i))/2., fTempA->GetBinContent(i+1));
 307             fv2VZEROC->Fill(fPtBins->At(i)+(fPtBins->At(i)+fPtBins->At(1+i))/2., fTempC->GetBinContent(i+1));
 308         }
 309     }
 310 }
 311  //_____________________________________________________________________________
 312 void AliAnalysisTaskJetFlow::DoGappedQC2Analysis()
 313 {
 314     // do q-cumulant analysis with eta gaps (avoiding autocorrelation of rps and jet constituents)
 315     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 316     // first step, get lhs tpc rp's
 317     fRhoVn->SetTrackEtaLimits(-0.9, -0.7);
 318     // get LHS rp's multiplicity and q-vector
 319     Double_t LHSreQn(0), LHSimQn(0), LHSmQ(0);
 320     (fDoPtWeight) ? fRhoVn->QCnQnk(2, 1, LHSreQn, LHSimQn) : fRhoVn->QCnQnk(2, 0, LHSreQn, LHSimQn);
 321     (fDoPtWeight) ? LHSmQ = fRhoVn->QCnM11() : LHSmQ = fRhoVn->QCnM();
 322     // get the RHS rp's multiplicity and q-vector
 323     fRhoVn->SetTrackEtaLimits(0.7, 0.9);
 324     Double_t RHSreQn(0), RHSimQn(0), RHSmQ(0);
 325     (fDoPtWeight) ? fRhoVn->QCnQnk(2, 1, RHSreQn, RHSimQn) : fRhoVn->QCnQnk(2, 0, RHSreQn, RHSimQn);
 326     (fDoPtWeight) ? RHSmQ = fRhoVn->QCnM11() : RHSmQ = fRhoVn->QCnM();
 327     // differential flow vectors
 328     Double_t repn[fPtBins->GetSize()-1];      // real part of q-vector of all poi's
 329     Double_t impn[fPtBins->GetSize()-1];      // im part of q-vector of all poi's
 330     Double_t mp[fPtBins->GetSize()-1];        // poi multiplicity
 331     Double_t reqn[fPtBins->GetSize()-1];      // real part of q-vectors of poi's labeled as rp
 332     Double_t imqn[fPtBins->GetSize()-1];      // im part of q-vectors of poi's labeled as rp
 333     Double_t mq[fPtBins->GetSize()-1];        // multiplicity of poi's labeled as rp
 334     for(Int_t i(0); i < fPtBins->GetSize(); i++) {
 335         repn[i] = 0;
 336         impn[i] = 0;
 337         mp[i] = 0;
 338         reqn[i] = 0;
 339         imqn[i] = 0;
 340         mq[i] = 0;
 341     }
 342     // calculate differential q-vectors and fill the profile with cumulants
 343     (fVParticleAnalysis) ? fRhoVn->SetTrackEtaLimits(-0.7, 0.7) : fRhoVn->SetLocalJetMinMaxEta(fJetRadius+.2);       // avoid overlap in poi and rp region
 344     QCnDifferentialFlowVectors(repn, impn, mp, reqn, imqn, mq, 2);
 345     // do the calculation
 346     if(RHSmQ*LHSmQ < 1) {
 347         fRhoVn->SetTrackEtaLimits(-0.9, 0.9);
 348         fRhoVn->SetLocalJetMinMaxEta(fJetRadius);
 349         return;
 350     }
 351     fRefCumulants->Fill(0., (LHSreQn*RHSreQn+LHSimQn*RHSimQn)/(RHSmQ*LHSmQ), RHSmQ*LHSmQ);
 352     for(Int_t i(0); i < fPtBins->GetSize(); i++) {
 353         if(LHSmQ*mp[i] < 1. ) continue;        // avoid division by zero
 354         Double_t atPt(fPtBins->At(i)+0.5*(fPtBins->At(i+1)-fPtBins->At(i)));      // pt value
 355         Double_t diffC((repn[i]*LHSreQn+impn[i]*LHSimQn)/(LHSmQ*mp[i]));
 356         Double_t eventW(mp[i]*LHSmQ);
 357         fDiffCumlantsV2->Fill(atPt, diffC, eventW);
 358     }
 359     // last step: roll back the eta cuts of the EmcalTask
 360     fRhoVn->SetTrackEtaLimits(-0.9, 0.9);
 361     fRhoVn->SetLocalJetMinMaxEta(fJetRadius);
 362 }
 363 //_____________________________________________________________________________
 364 void AliAnalysisTaskJetFlow::DoQC2FlowAnalysis()
 365 {
 366     // flow analysis with the qc2 method
 367     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 368     // reference flow is taken from the pico track selection and evaluated in
 369     // the AliAnalysisTaskiRhoVnModulation class
 370     for(Int_t harm(2); harm < 4; harm++) {        // loop over harmonics
 371         Double_t reQn(0), imQn(0), mQ(0);         // total q-vector
 372         // get the total q-vectors and reference cumulants
 373         (fDoPtWeight) ? fRhoVn->QCnQnk(harm, 1, reQn, imQn) : fRhoVn->QCnQnk(harm, 0, reQn, imQn);
 374         (fDoPtWeight) ? mQ = fRhoVn->QCnM11() : mQ = fRhoVn->QCnM();
 375         if(mQ < 2) continue;    // avoid division by zero
 376         fRefCumulants->Fill((double)(harm-2), ((reQn*reQn+imQn*imQn)-mQ)/(mQ*(mQ-1)), mQ*(mQ-1));
 377         // differential flow vectors
 378         Double_t repn[fPtBins->GetSize()-1];      // real part of q-vector of all poi's
 379         Double_t impn[fPtBins->GetSize()-1];      // im part of q-vector of all poi's
 380         Double_t mp[fPtBins->GetSize()-1];        // poi multiplicity
 381         Double_t reqn[fPtBins->GetSize()-1];      // real part of q-vectors of poi's labeled as rp
 382         Double_t imqn[fPtBins->GetSize()-1];      // im part of q-vectors of poi's labeled as rp
 383         Double_t mq[fPtBins->GetSize()-1];        // multiplicity of poi's labeled as rp
 384         for(Int_t i(0); i < fPtBins->GetSize(); i++) {
 385             repn[i] = 0;
 386             impn[i] = 0;
 387             mp[i] = 0;
 388             reqn[i] = 0;
 389             imqn[i] = 0;
 390             mq[i] = 0;
 391         }
 392         // calculate differential q-vectors and fill the profile with cumulants
 393         QCnDifferentialFlowVectors(repn, impn, mp, reqn, imqn, mq, harm);
 394         // FIXME differential evnet weights
 395         for(Int_t i(0); i < fPtBins->GetSize(); i++) {
 396             if(mp[i]*mQ - mq[i] <= 0 ) continue;        // avoid division by zero
 397             Double_t atPt(fPtBins->At(i)+0.5*(fPtBins->At(i+1)-fPtBins->At(i)));      // pt value
 398             Double_t diffC(((repn[i]*reQn+impn[i]*imQn)-mq[i])/(mp[i]*mQ-mq[i]));
 399             Double_t eventW(mp[i]*mQ-mq[i]);
 400             (harm == 2 ) ? fDiffCumlantsV2->Fill(atPt, diffC, eventW) : fDiffCumlantsV3->Fill(atPt, diffC, eventW);
 401         }
 402     }
 403 }
 404 //_____________________________________________________________________________
 405 void AliAnalysisTaskJetFlow::DoQC4FlowAnalysis()
 406 {
 407     // flow analysis with the qc4 method - see comments at qc2 FIXME not implemented yet
 408     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 409 }
 410 //_____________________________________________________________________________
 411 Bool_t AliAnalysisTaskJetFlow::DoFlowPackageFlowAnalysis()
 412 {
 413     // name's a bit misleading: this function does anlaysis using FlowAnalysis classes in /PWG/FLOW/Base/
 414     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 415     // get the jet array, which is added as an extension to the AliVEvent by the jetfinder
 416     Int_t nAcceptedJets(0), iPois(fPois->GetEntriesFast());
 417     if(iPois <= 0) return kFALSE;
 418     if(fFlowEvent_VZERO) {
 419         fCutsNull->SetEvent(InputEvent(), MCEvent());
 420         fCutsRP_VZERO->SetEvent(InputEvent(), MCEvent());
 421         fFlowEvent_VZERO->ClearFast();
 422         // the event is filled with rp's only, poi's will be added manually
 423         fFlowEvent_VZERO->Fill(fCutsRP_VZERO, fCutsNull);
 424         fFlowEvent_VZERO->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
 425     }
 426     if(fFlowEvent_TPC) {
 427         fCutsNull->SetEvent(InputEvent(), MCEvent());
 428         // in this case, both poi's and rp's will be added manually
 429         fFlowEvent_TPC->ClearFast();
 430         fFlowEvent_TPC->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
 431     }
 432     // loop over jets and inject them as POI's
 433     if(fVParticleAnalysis) {
 434         for(Int_t i(0); i < iPois; i++) {
 435             AliVTrack* poi = static_cast<AliVTrack*>(fPois->At(i));
 436             if(fRhoVn->PassesCuts(poi)) {
 437                 nAcceptedJets++;
 438                 // AliFlowTracks are created on the stack
 439                 AliFlowTrack flowTrack = AliFlowTrack(poi);
 440                 flowTrack.SetForPOISelection(kTRUE);
 441                 if(fFlowEvent_TPC) {
 442                     fFlowEvent_TPC->SetNumberOfRPs(fFlowEvent_TPC->GetNumberOfRPs()+1);
 443                     fFlowEvent_TPC->SetNumberOfPOIs(fFlowEvent_TPC->GetNumberOfPOIs()+1);
 444                     flowTrack.SetForRPSelection(kTRUE);
 445                     fFlowEvent_TPC->InsertTrack(&flowTrack);
 446                 }
 447                 if(fFlowEvent_VZERO) {
 448                     flowTrack.SetForRPSelection(kFALSE);
 449                     fFlowEvent_VZERO->InsertTrack(&flowTrack);
 450                     fFlowEvent_VZERO->SetNumberOfPOIs(fFlowEvent_VZERO->GetNumberOfPOIs()+1);
 451                 }
 452             }
 453         }
 454     } else {
 455         // add the jets as pois
 456         for(Int_t i(0); i < iPois; i++) {
 457             AliEmcalJet* poi = static_cast<AliEmcalJet*>(fPois->At(i));
 458             if(fRhoVn->PassesCuts(poi)) {
 459                 nAcceptedJets++;
 460                 AliFlowTrack flowTrack = AliFlowTrack(poi);
 461                 AliLocalRhoParameter* localRho = static_cast<AliLocalRhoParameter*>(fLocalRho);
 462                 if(!localRho) break;
 463                 Double_t rho(localRho->GetLocalVal(poi->Phi(), fJetRadius, localRho->GetVal()));
 464                 flowTrack.SetPt(poi->Pt() - poi->Area() * rho);
 465                 flowTrack.SetForPOISelection(kTRUE);
 466                 flowTrack.SetForRPSelection(kFALSE);
 467                 if(fFlowEvent_TPC) {
 468                     fFlowEvent_TPC->InsertTrack(&flowTrack);
 469                     fFlowEvent_TPC->SetNumberOfPOIs(fFlowEvent_TPC->GetNumberOfPOIs()+1);
 470                 }
 471                 if(fFlowEvent_VZERO) {
 472                     fFlowEvent_VZERO->InsertTrack(&flowTrack);
 473                     fFlowEvent_VZERO->SetNumberOfPOIs(fFlowEvent_VZERO->GetNumberOfPOIs()+1);
 474                 }
 475             }
 476         }
 477         // then add the reference section only for the TPC reference case
 478         for(Int_t i(0); i < fRPs->GetEntriesFast(); i++) {
 479             AliVTrack* rp = static_cast<AliVTrack*>(fRPs->At(i));
 480             if(fRhoVn->PassesCuts(rp) && rp->Pt() >= .15 && rp->Pt() <= 5.) {
 481                 AliFlowTrack flowTrack = AliFlowTrack(rp);
 482                 flowTrack.SetForPOISelection(kFALSE);
 483                 flowTrack.SetForRPSelection(kTRUE);
 484                 if(fFlowEvent_TPC) {
 485                     fFlowEvent_TPC->SetNumberOfRPs(fFlowEvent_TPC->GetNumberOfRPs()+1);
 486                     fFlowEvent_TPC->InsertTrack(&flowTrack);
 487                 }
 488             }
 489         }
 490     }
 491     if(fFlowEvent_TPC)      fFlowEvent_TPC->TagSubeventsInEta(-10, -1, 1, 10);
 492     if(fFlowEvent_VZERO)    fFlowEvent_VZERO->TagSubeventsInEta(-10, -1, 1, 10);
 493     return (nAcceptedJets < 1) ? kFALSE : kTRUE;
 494 }
 495 //_____________________________________________________________________________
 496 void AliAnalysisTaskJetFlow::QCnDifferentialFlowVectors(Double_t* repn, Double_t* impn, Double_t *mp, Double_t *reqn, Double_t *imqn, Double_t* mq, Int_t n)
 497 {
 498     // get (for now) unweighted differential flow vectors
 499     // FIXME move (part of) this code to AliAnalysisTaskRhoVnModulation
 500     Int_t iPois(fPois->GetEntriesFast());
 501     if(fVParticleAnalysis) {
 502         for(Int_t i(0); i < iPois; i++) {
 503             for(Int_t ptBin(0); ptBin < fPtBins->GetSize()-1; ptBin++) {
 504                 AliVTrack* poi = static_cast<AliVTrack*>(fPois->At(i));
 505                 if(fRhoVn->PassesCuts(poi)) {       // inherit cuts from mother task
 506                     if(poi->Pt() >= fPtBins->At(ptBin) && poi->Pt() < fPtBins->At(ptBin+1)) {
 507                             // fill the flow vectors assuming that all poi's are in the rp selection (true by design)
 508                             repn[ptBin]+=TMath::Cos(((Double_t)n)*poi->Phi());
 509                             impn[ptBin]+=TMath::Sin(((Double_t)n)*poi->Phi());
 510                             mp[ptBin]++;
 511                             reqn[ptBin]+=TMath::Cos(((Double_t)n)*poi->Phi());
 512                             imqn[ptBin]+=TMath::Sin(((Double_t)n)*poi->Phi());
 513                             mq[ptBin]++;
 514                     }
 515                 }
 516             }
 517         }
 518     } else {
 519         for(Int_t i(0); i < iPois; i++) {
 520             for(Int_t ptBin(0); ptBin < fPtBins->GetSize()-1; ptBin++) {
 521                 AliEmcalJet* poi = static_cast<AliEmcalJet*>(fPois->At(i));
 522                 AliLocalRhoParameter* localRho = static_cast<AliLocalRhoParameter*>(fLocalRho);
 523                 if(!localRho) break;
 524                 Double_t rho(localRho->GetLocalVal(poi->Phi(), fJetRadius, localRho->GetVal()));
 525                 Double_t pt(poi->Pt() - poi->Area() * rho);
 526                 if(fRhoVn->PassesCuts(poi)) {
 527                     if(pt >= fPtBins->At(ptBin) && pt < fPtBins->At(ptBin+1)) {
 528                             // fill the flow vectors assuming that all poi's are in the rp selection (true by design)
 529                             repn[ptBin]+=TMath::Cos(((Double_t)n)*poi->Phi());
 530                             impn[ptBin]+=TMath::Sin(((Double_t)n)*poi->Phi());
 531                             mp[ptBin]++;        // qn isn't filled, no overlap between poi's and rp's
 532                     }
 533                 }
 534             }
 535         }
 536     }
 537 }
 538 //_____________________________________________________________________________
 539 void AliAnalysisTaskJetFlow::Terminate(Option_t *)
 540 {
 541     // terminate
 542     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 543 }
 544 //_____________________________________________________________________________