]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliAnalysisTaskJetFlow.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' into dev
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskJetFlow.cxx
CommitLineData
e1dc3f11 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/*
6f5e34fa 17 * AliAnaysisTaskJetFlow
e1dc3f11 18 * author: Redmer Alexander Bertens
19 * rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
20 *
6f5e34fa 21 * Interface task between EMCal jet framework and flow package
e1dc3f11 22 *
6f5e34fa 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 *
7f232b07 28 * POI's can be supplied as AliEmcalJets or as AliVTracks, note the different behavior
6f5e34fa 29 * with respect to the Pt value: for AliEmcalJets subtracted Pt is used by default, for VParticles
30 * Pt is taken directly.
e1dc3f11 31 *
32 * */
33
6f5e34fa 34// root includes
e1dc3f11 35#include <TChain.h>
36#include <TH1F.h>
6f5e34fa 37#include <TArrayD.h>
38#include <TProfile.h>
39#include <TString.h>
40#include <TList.h>
7f232b07 41#include <TClonesArray.h>
6f5e34fa 42// aliroot includes
43#include <AliAODEvent.h>
44#include <AliESDEvent.h>
45#include <AliEmcalJet.h>
e1dc3f11 46#include <AliAnalysisTask.h>
47#include <AliAnalysisManager.h>
48#include <AliFlowEventSimple.h>
6f5e34fa 49#include <AliFlowEvent.h>
e1dc3f11 50#include <AliFlowTrack.h>
51#include <AliFlowTrackCuts.h>
2045ca8b 52#include <AliFlowEventCuts.h>
e1dc3f11 53#include <AliFlowCommonConstants.h>
6f5e34fa 54// local includes
e1dc3f11 55#include "AliAnalysisTaskJetFlow.h"
7f232b07 56// EMCAL jet framework includes
57#include <AliRhoParameter.h>
ba9aeb09 58#include <AliLocalRhoParameter.h>
7f232b07 59#include <AliPicoTrack.h>
60#include <AliAnalysisTaskRhoVnModulation.h>
e1dc3f11 61
62class AliAnalysisTaskJetFlow;
63
64using namespace std;
65
66ClassImp(AliAnalysisTaskJetFlow)
67
68AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow() : AliAnalysisTaskSE(),
ba9aeb09 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)
6f5e34fa 70{ /* default constructor for ROOT IO */ }
e1dc3f11 71//_____________________________________________________________________________
06e0f8e5 72AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow(
73 const char* name,
7f232b07 74 AliAnalysisTaskRhoVnModulation* rhoTask,
75 Bool_t VPart,
ba9aeb09 76 Bool_t VZEROEP,
77 Bool_t GQC2,
7f232b07 78 Bool_t QC2,
79 Bool_t QC4,
80 Bool_t FlowPackageSP,
81 Bool_t FlowPackageQC
82 ) : AliAnalysisTaskSE(name),
ad0edaf5 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)
e1dc3f11 84{
85 // constructor
86 DefineInput(0, TChain::Class());
87 DefineOutput(1, TList::Class());
7f232b07 88 fJetsName = rhoTask->GetJetsName();
89 fTracksName = rhoTask->GetTracksName();
ba9aeb09 90 if(GQC2 && QC2) {
91 printf(" > Warning, QC2 and gapped QC2 method are both called <\n will only run gapped QC2 !");
a4b3e9e4 92 fDoQC2FlowAnalysis = kFALSE;
ba9aeb09 93 }
7f232b07 94 if(FlowPackageSP || FlowPackageQC) DefineOutput(2, AliFlowEventSimple::Class());
95 if(FlowPackageSP && FlowPackageQC) DefineOutput(3, AliFlowEventSimple::Class());
e1dc3f11 96}
97//_____________________________________________________________________________
98AliAnalysisTaskJetFlow::~AliAnalysisTaskJetFlow()
99{
100 // destructor
6f5e34fa 101 if(fOutputList) delete fOutputList;
102 if(fFlowEvent_TPC) delete fFlowEvent_TPC;
103 if(fFlowEvent_VZERO) delete fFlowEvent_VZERO;
104 if(fCutsEvent) delete fCutsEvent;
5a4d96c7 105 if(fCutsRP_VZERO) delete fCutsRP_VZERO;
5a4d96c7 106 if(fCutsNull) delete fCutsNull;
5a4d96c7 107 if(fPtBins) delete fPtBins;
108 if(fTempA) delete fTempA;
109 if(fTempC) delete fTempC;
e1dc3f11 110}
111//_____________________________________________________________________________
112void AliAnalysisTaskJetFlow::UserCreateOutputObjects()
113{
114 // create output objects
115 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2b7447f0 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 }
51300ca6 129 if(fDoQCFPAnalysis) fFlowEvent_TPC = new AliFlowEvent(10000);
2b7447f0 130 }
e1dc3f11 131 fOutputList = new TList();
132 fOutputList->SetOwner(kTRUE);
98bd4878 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);
98bd4878 141 fOutputList->Add(fHistAnalysisSummary);
2b7447f0 142 if(fVParticleAnalysis && ! fPtBins) { // FIXME inherits from flow package now
7f232b07 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());
6f5e34fa 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);
06e0f8e5 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);
5a4d96c7 164 // bookkeeping histo's, will not be stored
165 fTempA = (TProfile*)fv2VZEROA->Clone("temp_a");
166 fTempC = (TProfile*)fv2VZEROC->Clone("temp_c");
6f5e34fa 167 }
ba9aeb09 168 if(fDoQC2FlowAnalysis || fDoGappedQC2Analysis) {
7f232b07 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);
2b7447f0 173 fDiffCumlantsV2 = new TProfile("Differential cumulants v2", "Differential cumulants v2", fPtBins->GetSize()-1, fPtBins->GetArray());
7f232b07 174 fDiffCumlantsV2->GetXaxis()->SetTitle("Pt [GeV/c]");
175 fDiffCumlantsV2->GetYaxis()->SetTitle("v_{2}[2]");
176 fOutputList->Add(fDiffCumlantsV2);
2b7447f0 177 fDiffCumlantsV3 = new TProfile("Differential cumulants v3", "Differential cumulants v3", fPtBins->GetSize()-1, fPtBins->GetArray());
7f232b07 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 }
a4a06f0e 188 // qa
189 fCentralitySelection = new TH1F("fCentralitySelection", "fCentralitySelection", 100, 0, 100);
190 fOutputList->Add(fCentralitySelection);
e1dc3f11 191 PostData(1, fOutputList);
192 // create the flow event and configure the static cc object
f0438bd7 193 Bool_t slotTwoFilled(kFALSE);
7f232b07 194 if(fFlowEvent_VZERO) {
f0438bd7 195 PostData(2, fFlowEvent_VZERO);
196 slotTwoFilled = kTRUE;
197 }
7f232b07 198 if(fFlowEvent_TPC) {
f0438bd7 199 (slotTwoFilled) ? PostData(3, fFlowEvent_TPC) : PostData(2, fFlowEvent_TPC);
200 }
201 if(fFlowEvent_VZERO || fFlowEvent_TPC) {
202 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
7f232b07 203 cc->SetPtMin(fCCMinPt);
204 cc->SetPtMax(fCCMaxPt);
f0438bd7 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 }
06e0f8e5 213 }
e1dc3f11 214}
215//_____________________________________________________________________________
216void AliAnalysisTaskJetFlow::UserExec(Option_t *)
217{
218 // user exec
219 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
f0438bd7 220 if(!InputEvent()) return; // coverity (and sanity)
7f232b07 221 // see if the analysis is initialized
a4a06f0e 222 if(!fInitialized) {
223 if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
224 else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
7f232b07 225 (fVParticleAnalysis) ? fPois = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName.Data())) : fPois = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName.Data()));
ea4212a0 226 fRPs = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName.Data()));
227 if(!fPois || !fRPs) return; // couldn't get expected input data
a4a06f0e 228 fInitialized = kTRUE;
ba9aeb09 229 fLocalRho = InputEvent()->FindListObject(fLocalRhoName.Data());
230 if(!fLocalRho && !fVParticleAnalysis) {
231 AliFatal(Form("Couldn't find %s, aborting!", fLocalRhoName.Data()));
232 }
a4a06f0e 233 }
7f232b07 234 if(!PassesCuts()) return; // event quality cuts and centrality determination
235 // execute the requested flow methods
236 if(fDoVZEROFlowAnalysis) DoVZEROFlowAnalysis();
ba9aeb09 237 if(fDoGappedQC2Analysis) DoGappedQC2Analysis();
7f232b07 238 if(fDoQC2FlowAnalysis) DoQC2FlowAnalysis();
239 if(fDoQC4FlowAnalysis) DoQC4FlowAnalysis();
240 Bool_t post(0x0); // post only when analysis succeeded
2b7447f0 241 if(fFlowEvent_TPC || fFlowEvent_VZERO) DoFlowPackageFlowAnalysis();
7f232b07 242 // push the output for the different analyses to file
243 PostData(1, fOutputList);
2b7447f0 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);
7f232b07 250 }
251}
252//_____________________________________________________________________________
253Bool_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//_____________________________________________________________________________
263void 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)));
98bd4878 283 }
2db2a6d0 284 }
e1dc3f11 285 }
7f232b07 286 } else {
287 for(Int_t i(0); i < iPois; i++) {
288 AliEmcalJet* poi = static_cast<AliEmcalJet*>(fPois->At(i));
ea4212a0 289 if(fRhoVn->PassesCuts(poi)) {
ba9aeb09 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);
7f232b07 294 if(!fDoMultWeight) {
ba9aeb09 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)));
7f232b07 297 } else {
ba9aeb09 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)));
7f232b07 300 }
301 }
f0438bd7 302 }
f0438bd7 303 }
7f232b07 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 }
f0438bd7 309 }
ba9aeb09 310}
311 //_____________________________________________________________________________
312void 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
a4b3e9e4 343 (fVParticleAnalysis) ? fRhoVn->SetTrackEtaLimits(-0.7, 0.7) : fRhoVn->SetLocalJetMinMaxEta(fJetRadius+.2); // avoid overlap in poi and rp region
ba9aeb09 344 QCnDifferentialFlowVectors(repn, impn, mp, reqn, imqn, mq, 2);
345 // do the calculation
ad0edaf5 346 if(RHSmQ*LHSmQ < 1) {
347 fRhoVn->SetTrackEtaLimits(-0.9, 0.9);
348 fRhoVn->SetLocalJetMinMaxEta(fJetRadius);
349 return;
350 }
ba9aeb09 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);
e1dc3f11 362}
363//_____________________________________________________________________________
7f232b07 364void AliAnalysisTaskJetFlow::DoQC2FlowAnalysis()
98bd4878 365{
7f232b07 366 // flow analysis with the qc2 method
98bd4878 367 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
7f232b07 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
ea4212a0 393 QCnDifferentialFlowVectors(repn, impn, mp, reqn, imqn, mq, harm);
7f232b07 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]);
ea4212a0 400 (harm == 2 ) ? fDiffCumlantsV2->Fill(atPt, diffC, eventW) : fDiffCumlantsV3->Fill(atPt, diffC, eventW);
7f232b07 401 }
402 }
98bd4878 403}
404//_____________________________________________________________________________
7f232b07 405void AliAnalysisTaskJetFlow::DoQC4FlowAnalysis()
5a4d96c7 406{
7f232b07 407 // flow analysis with the qc4 method - see comments at qc2 FIXME not implemented yet
5a4d96c7 408 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
5a4d96c7 409}
410//_____________________________________________________________________________
7f232b07 411Bool_t AliAnalysisTaskJetFlow::DoFlowPackageFlowAnalysis()
6f5e34fa 412{
7f232b07 413 // name's a bit misleading: this function does anlaysis using FlowAnalysis classes in /PWG/FLOW/Base/
6f5e34fa 414 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
7f232b07 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);
ea4212a0 450 fFlowEvent_VZERO->SetNumberOfPOIs(fFlowEvent_VZERO->GetNumberOfPOIs()+1);
7f232b07 451 }
452 }
453 }
454 } else {
ea4212a0 455 // add the jets as pois
7f232b07 456 for(Int_t i(0); i < iPois; i++) {
457 AliEmcalJet* poi = static_cast<AliEmcalJet*>(fPois->At(i));
ea4212a0 458 if(fRhoVn->PassesCuts(poi)) {
7f232b07 459 nAcceptedJets++;
460 AliFlowTrack flowTrack = AliFlowTrack(poi);
ba9aeb09 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);
7f232b07 465 flowTrack.SetForPOISelection(kTRUE);
466 flowTrack.SetForRPSelection(kFALSE);
ea4212a0 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 }
6f5e34fa 488 }
7f232b07 489 }
490 }
51300ca6 491 if(fFlowEvent_TPC) fFlowEvent_TPC->TagSubeventsInEta(-10, -1, 1, 10);
492 if(fFlowEvent_VZERO) fFlowEvent_VZERO->TagSubeventsInEta(-10, -1, 1, 10);
7f232b07 493 return (nAcceptedJets < 1) ? kFALSE : kTRUE;
494}
495//_____________________________________________________________________________
ea4212a0 496void AliAnalysisTaskJetFlow::QCnDifferentialFlowVectors(Double_t* repn, Double_t* impn, Double_t *mp, Double_t *reqn, Double_t *imqn, Double_t* mq, Int_t n)
7f232b07 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]++;
06e0f8e5 514 }
6f5e34fa 515 }
516 }
517 }
7f232b07 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));
ba9aeb09 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)) {
7f232b07 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 }
06e0f8e5 534 }
6f5e34fa 535 }
6f5e34fa 536 }
537}
538//_____________________________________________________________________________
e1dc3f11 539void 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//_____________________________________________________________________________