]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add delta area to THnSparse for delta pT output
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
CommitLineData
3b7ffecf 1// **************************************
3170a3f8 2// used for the correction of determiantion of reconstructed jet spectra
3b7ffecf 3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
21
22
23#include <TROOT.h>
24#include <TRandom.h>
742ee86c 25#include <TRandom3.h>
3b7ffecf 26#include <TSystem.h>
27#include <TInterpreter.h>
28#include <TChain.h>
29#include <TFile.h>
30#include <TKey.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TH3F.h>
34#include <TProfile.h>
35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
37eb26ea 38#include <TRefArray.h>
3b7ffecf 39#include "TDatabasePDG.h"
40
41#include "AliAnalysisTaskJetSpectrum2.h"
42#include "AliAnalysisManager.h"
43#include "AliJetFinder.h"
44#include "AliJetHeader.h"
45#include "AliJetReader.h"
46#include "AliJetReaderHeader.h"
47#include "AliUA1JetHeaderV1.h"
48#include "AliESDEvent.h"
49#include "AliAODEvent.h"
50#include "AliAODHandler.h"
51#include "AliAODTrack.h"
52#include "AliAODJet.h"
c2785065 53#include "AliAODJetEventBackground.h"
3b7ffecf 54#include "AliAODMCParticle.h"
55#include "AliMCEventHandler.h"
56#include "AliMCEvent.h"
57#include "AliStack.h"
58#include "AliGenPythiaEventHeader.h"
59#include "AliJetKineReaderHeader.h"
60#include "AliGenCocktailEventHeader.h"
61#include "AliInputEventHandler.h"
62
63
64#include "AliAnalysisHelperJetTasks.h"
65
66ClassImp(AliAnalysisTaskJetSpectrum2)
67
37eb26ea 68AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
69 AliAnalysisTaskSE(),
a31b8a87 70 fJetHeaderRec(0x0),
71 fJetHeaderGen(0x0),
72 fAODIn(0x0),
73 fAODOut(0x0),
74 fAODExtension(0x0),
75 fhnCorrelation(0x0),
c8eabe24 76 f1PtScale(0x0),
a31b8a87 77 fBranchRec("jets"),
78 fBranchGen(""),
37eb26ea 79 fBranchBkgRec(""),
80 fBranchBkgGen(""),
a31b8a87 81 fNonStdFile(""),
742ee86c 82 fRandomizer(0x0),
a31b8a87 83 fUseAODJetInput(kFALSE),
84 fUseAODTrackInput(kFALSE),
85 fUseAODMCInput(kFALSE),
86 fUseGlobalSelection(kFALSE),
87 fUseExternalWeightOnly(kFALSE),
88 fLimitGenJetEta(kFALSE),
a31b8a87 89 fNMatchJets(5),
7b822bc5 90 fNRPBins(3),
a31b8a87 91 fFilterMask(0),
f2dd0695 92 fEventSelectionMask(0),
a31b8a87 93 fAnalysisType(0),
3b7ffecf 94 fTrackTypeRec(kTrackUndef),
95 fTrackTypeGen(kTrackUndef),
f4132e7d 96 fEventClass(0),
742ee86c 97 fRPSubeventMethod(0),
3b7ffecf 98 fAvgTrials(1),
99 fExternalWeight(1),
a31b8a87 100 fJetRecEtaWindow(0.5),
101 fTrackRecEtaWindow(0.5),
9280dfa6 102 fMinJetPt(0),
a31b8a87 103 fMinTrackPt(0.15),
104 fDeltaPhiWindow(90./180.*TMath::Pi()),
742ee86c 105 fCentrality(100),
7b822bc5 106 fRPAngle(0),
37eb26ea 107 fMultRec(0),
108 fMultGen(0),
3b7ffecf 109 fh1Xsec(0x0),
110 fh1Trials(0x0),
111 fh1PtHard(0x0),
112 fh1PtHardNoW(0x0),
113 fh1PtHardTrials(0x0),
57ca1193 114 fh1ZVtx(0x0),
7b822bc5 115 fh1RP(0x0),
742ee86c 116 fh1Centrality(0x0),
c8eabe24 117 fh1TmpRho(0x0),
37eb26ea 118 fh2MultRec(0x0),
119 fh2MultGen(0x0),
742ee86c 120 fh2RPSubevents(0x0),
121 fh2RPCentrality(0x0),
122 fh2RPDeltaRP(0x0),
123 fh2RPQxQy(0x0),
124 fh2RPCosDeltaRP(0x0),
a31b8a87 125 fh2PtFGen(0x0),
9a83d4af 126 fh2RelPtFGen(0x0),
c08c3ad2 127 fh3PhiWeights(0x0),
742ee86c 128 fh3RPPhiTracks(0x0),
3b7ffecf 129 fHistList(0x0)
130{
c17f867b 131
132 fFlatA[0] = fFlatA[1] = 0;
133 fFlatB[0] = fFlatB[1] = 0;
134 fDeltaQxy[0] = fDeltaQxy[1] = 0;
135
3b7ffecf 136 for(int i = 0;i < kMaxStep*2;++i){
137 fhnJetContainer[i] = 0;
138 }
a31b8a87 139
140 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 141 fFlagJetType[ij] = 1; // default = on
a31b8a87 142 fh1NJets[ij] = 0;
143 fh1SumPtTrack[ij] = 0;
144 fh1PtJetsIn[ij] = 0;
145 fh1PtTracksIn[ij] = 0;
cb883243 146 fh1PtTracksInLow[ij] = 0;
a31b8a87 147 fh1PtTracksLeadingIn[ij] = 0;
148 fh2NJetsPt[ij] = 0;
149 fh2NTracksPt[ij] = 0;
3170a3f8 150 fh2TrackEtaPt[ij] = 0;
7b822bc5 151 fh3MultTrackPtRP[ij] = 0;
d95fc15a 152 fh3MultTrackPtLowRP[ij] = 0;
a31b8a87 153 fh2LeadingTrackPtTrackPhi[ij] = 0;
37eb26ea 154 for(int i = 0;i <= kMaxJets;++i){
7b822bc5 155 fh3MultPtRP[ij][i] = 0;
f12de05f 156 fh2PhiPt[ij][i] = 0;
157 fh2EtaPt[ij][i] = 0;
37eb26ea 158 fh2AreaPt[ij][i] = 0;
159 fh2EtaArea[ij][i] = 0;
f12de05f 160 fh2PhiEta[ij][i] = 0;
d7117cbd 161 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 162 fh1PtIn[ij][i] = 0;
163 fh2RhoPt[ij][i] = 0;
164 fh2PsiPt[ij][i] = 0;
165 }
166
167 fh1DijetMinv[ij] = 0;
168 fh2DijetDeltaPhiPt[ij] = 0;
169 fh2DijetAsymPt[ij] = 0;
170 fh2DijetPt2vsPt1[ij] = 0;
171 fh2DijetDifvsSum[ij] = 0;
172
173 }
3b7ffecf 174}
175
176AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
177 AliAnalysisTaskSE(name),
178 fJetHeaderRec(0x0),
179 fJetHeaderGen(0x0),
a31b8a87 180 fAODIn(0x0),
181 fAODOut(0x0),
182 fAODExtension(0x0),
3b7ffecf 183 fhnCorrelation(0x0),
c8eabe24 184 f1PtScale(0x0),
3b7ffecf 185 fBranchRec("jets"),
186 fBranchGen(""),
37eb26ea 187 fBranchBkgRec(""),
188 fBranchBkgGen(""),
a31b8a87 189 fNonStdFile(""),
742ee86c 190 fRandomizer(0x0),
565584e8 191 fUseAODJetInput(kFALSE),
192 fUseAODTrackInput(kFALSE),
193 fUseAODMCInput(kFALSE),
b5cc0c6d 194 fUseGlobalSelection(kFALSE),
3b7ffecf 195 fUseExternalWeightOnly(kFALSE),
196 fLimitGenJetEta(kFALSE),
a31b8a87 197 fNMatchJets(5),
7b822bc5 198 fNRPBins(3),
3b7ffecf 199 fFilterMask(0),
f2dd0695 200 fEventSelectionMask(0),
3b7ffecf 201 fAnalysisType(0),
202 fTrackTypeRec(kTrackUndef),
203 fTrackTypeGen(kTrackUndef),
f4132e7d 204 fEventClass(0),
742ee86c 205 fRPSubeventMethod(0),
3b7ffecf 206 fAvgTrials(1),
207 fExternalWeight(1),
a31b8a87 208 fJetRecEtaWindow(0.5),
209 fTrackRecEtaWindow(0.5),
9280dfa6 210 fMinJetPt(0),
a31b8a87 211 fMinTrackPt(0.15),
212 fDeltaPhiWindow(90./180.*TMath::Pi()),
742ee86c 213 fCentrality(100),
7b822bc5 214 fRPAngle(0),
37eb26ea 215 fMultRec(0),
216 fMultGen(0),
3b7ffecf 217 fh1Xsec(0x0),
218 fh1Trials(0x0),
219 fh1PtHard(0x0),
220 fh1PtHardNoW(0x0),
221 fh1PtHardTrials(0x0),
57ca1193 222 fh1ZVtx(0x0),
7b822bc5 223 fh1RP(0x0),
742ee86c 224 fh1Centrality(0x0),
c8eabe24 225 fh1TmpRho(0x0),
37eb26ea 226 fh2MultRec(0x0),
227 fh2MultGen(0x0),
742ee86c 228 fh2RPSubevents(0x0),
229 fh2RPCentrality(0x0),
230 fh2RPDeltaRP(0x0),
231 fh2RPQxQy(0x0),
232 fh2RPCosDeltaRP(0x0),
a31b8a87 233 fh2PtFGen(0x0),
9a83d4af 234 fh2RelPtFGen(0x0),
c08c3ad2 235 fh3PhiWeights(0x0),
742ee86c 236 fh3RPPhiTracks(0x0),
3b7ffecf 237 fHistList(0x0)
238{
239
c17f867b 240 fFlatA[0] = fFlatA[1] = 0;
241 fFlatB[0] = fFlatB[1] = 0;
242 fDeltaQxy[0] = fDeltaQxy[1] = 0;
243
3b7ffecf 244 for(int i = 0;i < kMaxStep*2;++i){
245 fhnJetContainer[i] = 0;
246 }
a31b8a87 247
248 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 249 fFlagJetType[ij] = 1; // default = on
a31b8a87 250 fh1NJets[ij] = 0;
251 fh1SumPtTrack[ij] = 0;
252 fh1PtJetsIn[ij] = 0;
253 fh1PtTracksIn[ij] = 0;
cb883243 254 fh1PtTracksInLow[ij] = 0;
a31b8a87 255 fh1PtTracksLeadingIn[ij] = 0;
256 fh2NJetsPt[ij] = 0;
257 fh2NTracksPt[ij] = 0;
3170a3f8 258 fh2TrackEtaPt[ij] = 0;
7b822bc5 259 fh3MultTrackPtRP[ij] = 0;
d95fc15a 260 fh3MultTrackPtLowRP[ij] = 0;
a31b8a87 261 fh2LeadingTrackPtTrackPhi[ij] = 0;
37eb26ea 262 for(int i = 0;i <= kMaxJets;++i){
7b822bc5 263 fh3MultPtRP[ij][i] = 0;
f12de05f 264 fh2PhiPt[ij][i] = 0;
265 fh2EtaPt[ij][i] = 0;
37eb26ea 266 fh2AreaPt[ij][i] = 0;
267 fh2EtaArea[ij][i] = 0;
f12de05f 268 fh2PhiEta[ij][i] = 0;
d7117cbd 269 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 270 fh1PtIn[ij][i] = 0;
271 fh2RhoPt[ij][i] = 0;
272 fh2PsiPt[ij][i] = 0;
273 }
274
275 fh1DijetMinv[ij] = 0;
276 fh2DijetDeltaPhiPt[ij] = 0;
277 fh2DijetAsymPt[ij] = 0;
278 fh2DijetPt2vsPt1[ij] = 0;
279 fh2DijetDifvsSum[ij] = 0;
280 }
281
37eb26ea 282 DefineOutput(1, TList::Class());
3b7ffecf 283}
284
285
286
287Bool_t AliAnalysisTaskJetSpectrum2::Notify()
288{
a31b8a87 289
290
291
3b7ffecf 292 //
293 // Implemented Notify() to read the cross sections
294 // and number of trials from pyxsec.root
295 //
a31b8a87 296
297 // Fetch the aod also from the input in,
298 // have todo it in notify
299
300
301 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
302 // assume that the AOD is in the general output...
303 fAODOut = AODEvent();
304
305 if(fNonStdFile.Length()!=0){
306 // case that we have an AOD extension we need can fetch the jets from the extended output
307 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
e855f5c5 308 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
a31b8a87 309 if(!fAODExtension){
310 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
311 }
312 }
313
3b7ffecf 314
315 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
316 Float_t xsection = 0;
317 Float_t ftrials = 1;
318
319 fAvgTrials = 1;
320 if(tree){
321 TFile *curfile = tree->GetCurrentFile();
322 if (!curfile) {
323 Error("Notify","No current file");
324 return kFALSE;
325 }
326 if(!fh1Xsec||!fh1Trials){
327 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
328 return kFALSE;
329 }
330 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
331 fh1Xsec->Fill("<#sigma>",xsection);
332 // construct a poor man average trials
333 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 334 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
3b7ffecf 335 }
a31b8a87 336
d7117cbd 337 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
a31b8a87 338
3b7ffecf 339 return kTRUE;
340}
341
342void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
343{
344
742ee86c 345
3b7ffecf 346 // Connect the AOD
347
3b7ffecf 348 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
3b7ffecf 349 OpenFile(1);
37eb26ea 350 if(!fHistList)fHistList = new TList();
40440b28 351 PostData(1, fHistList); // post data in any case once
352
742ee86c 353 if(!fRandomizer)fRandomizer = new TRandom3(0);
354
5cb0f438 355 fHistList->SetOwner(kTRUE);
37eb26ea 356 Bool_t oldStatus = TH1::AddDirectoryStatus();
3b7ffecf 357 TH1::AddDirectory(kFALSE);
358
a31b8a87 359 MakeJetContainer();
360 fHistList->Add(fhnCorrelation);
b5d90baf 361 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
a31b8a87 362
3b7ffecf 363 //
364 // Histogram
365
d95fc15a 366 const Int_t nBinPt = 125;
3b7ffecf 367 Double_t binLimitsPt[nBinPt+1];
368 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
369 if(iPt == 0){
370 binLimitsPt[iPt] = 0.0;
371 }
372 else {// 1.0
742ee86c 373 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
3b7ffecf 374 }
375 }
edfbe476 376 const Int_t nBinPhi = 90;
3b7ffecf 377 Double_t binLimitsPhi[nBinPhi+1];
378 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
379 if(iPhi==0){
37eb26ea 380 binLimitsPhi[iPhi] = 0;
3b7ffecf 381 }
382 else{
383 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
384 }
385 }
386
edfbe476 387 const Int_t nBinEta = 40;
388 Double_t binLimitsEta[nBinEta+1];
389 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
390 if(iEta==0){
391 binLimitsEta[iEta] = -2.0;
392 }
393 else{
cf049e94 394 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 395 }
396 }
397
a31b8a87 398
3b7ffecf 399 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
400 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
a31b8a87 401 fHistList->Add(fh1Xsec);
3b7ffecf 402 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
403 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
a31b8a87 404 fHistList->Add(fh1Trials);
3b7ffecf 405 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 406 fHistList->Add(fh1PtHard);
3b7ffecf 407 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 408 fHistList->Add(fh1PtHardNoW);
3b7ffecf 409 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 410 fHistList->Add(fh1PtHardTrials);
a31b8a87 411
57ca1193 412 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
a31b8a87 413 fHistList->Add(fh1ZVtx);
7b822bc5 414
9939142f 415 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
7b822bc5 416 fHistList->Add(fh1RP);
417
742ee86c 418 fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",101,-0.5,100.5);
419 fHistList->Add(fh1Centrality);
420
3170a3f8 421 fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
37eb26ea 422 fHistList->Add(fh2MultRec);
3170a3f8 423 fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
37eb26ea 424 fHistList->Add(fh2MultGen);
a31b8a87 425
742ee86c 426 fh2RPSubevents = new TH2F("fh2RPSubevents" ,"Reaction Plane Angle" , 180, 0, TMath::Pi(), 180, 0, TMath::Pi());
427 fHistList->Add( fh2RPSubevents);
a31b8a87 428
742ee86c 429 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
430 fHistList->Add(fh2RPCentrality);
a31b8a87 431
742ee86c 432 fh2RPDeltaRP = new TH2F("fh2DeltaRP" ,"Delta Reaction Plane Angle" , 100, -TMath::Pi()/2, TMath::Pi()/2,20,0.,100.0);
433 fHistList->Add(fh2RPDeltaRP);
a31b8a87 434
742ee86c 435 fh2RPQxQy = new TH2F("fh2RPQxQy" ,"" , 100, -100,100,100,-100,100);
436 fHistList->Add(fh2RPQxQy);
a31b8a87 437
742ee86c 438 fh2RPCosDeltaRP = new TH2F("fh2RPCosDeltaRP" ,"" , 20, 0.001,100.001,100,-1,1);
439 fHistList->Add(fh2RPCosDeltaRP);
cb883243 440
a31b8a87 441
742ee86c 442 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
443 fHistList->Add(fh2PtFGen);
a31b8a87 444
742ee86c 445 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
446 fHistList->Add(fh2RelPtFGen);
a31b8a87 447
d95fc15a 448 fh3RPPhiTracks = new TH3F("fh3RPPhiTracks","Phi Tracks Pt Centrality", 10, 0.,100.,20,-5,5,200, 0, 2*TMath::Pi());
742ee86c 449 fHistList->Add(fh3RPPhiTracks);
450
451 for(int ij = 0;ij <kJetTypes;++ij){
452 TString cAdd = "";
453 TString cJetBranch = "";
454 if(ij==kJetRec){
455 cAdd = "Rec";
456 cJetBranch = fBranchRec.Data();
457 }
458 else if (ij==kJetGen){
459 cAdd = "Gen";
460 cJetBranch = fBranchGen.Data();
461 }
462 else if (ij==kJetRecFull){
463 cAdd = "RecFull";
464 cJetBranch = fBranchRec.Data();
465 }
466 else if (ij==kJetGenFull){
467 cAdd = "GenFull";
468 cJetBranch = fBranchGen.Data();
469 }
a31b8a87 470
742ee86c 471 if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
472 if(!fFlagJetType[ij])continue;
7b822bc5 473
742ee86c 474 fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
475 fHistList->Add(fh1NJets[ij]);
476
477 fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
478 fHistList->Add(fh1PtJetsIn[ij]);
479
480 fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
481 fHistList->Add(fh1PtTracksIn[ij]);
482
d95fc15a 483 fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),50,0.,5.);
742ee86c 484 fHistList->Add(fh1PtTracksInLow[ij]);
485
486 fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
487 fHistList->Add(fh1PtTracksLeadingIn[ij]);
488
3170a3f8 489 fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
742ee86c 490 fHistList->Add(fh1SumPtTrack[ij]);
491
492 fh2NJetsPt[ij] = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
493 fHistList->Add(fh2NJetsPt[ij]);
494
3170a3f8 495 fh2NTracksPt[ij] = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,5000);
742ee86c 496 fHistList->Add(fh2NTracksPt[ij]);
497
3170a3f8 498
499 fh2TrackEtaPt[ij] = new TH2F(Form("fh2TrackEtaPt%s",cAdd.Data()),";#eta;p_{T,track}",nBinPt,binLimitsPt,50,-1.,1.);
500 fHistList->Add(fh2TrackEtaPt[ij]);
501
742ee86c 502 fh3MultTrackPtRP[ij] = new TH3F(Form("fh3MultTrackPtRP%s",
d95fc15a 503 cAdd.Data()),Form(";# tracks;%s track p_T;RP Bin",cAdd.Data()),250,0,5000,50,0.,50.,(Int_t)fNRPBins,-0.5,fNRPBins-0.5);
742ee86c 504 fHistList->Add(fh3MultTrackPtRP[ij]);
d95fc15a 505
506 fh3MultTrackPtLowRP[ij] = new TH3F(Form("fh3MultTrackPtLowRP%s",
507 cAdd.Data()),Form(";# tracks;%s track p_T;RP Bin",cAdd.Data()),250,0,5000,50,0.,5.,(Int_t)fNRPBins,-0.5,fNRPBins-0.5);
508 fHistList->Add(fh3MultTrackPtLowRP[ij]);
742ee86c 509
510 fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
bd2402ef 511 nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
a31b8a87 512 fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
a31b8a87 513
37eb26ea 514 for(int i = 0;i <= kMaxJets;++i){
a31b8a87 515 fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
516 fHistList->Add(fh1PtIn[ij][i]);
517
518 fh2RhoPt[ij][i] = new TH2F(Form("fh2RhoPt%s_j%d",cAdd.Data(),i),Form("jet shape rho for %s jets;r;p_{T};",cAdd.Data()),
519 40,0.,2.,nBinPt,binLimitsPt);
520 fHistList->Add(fh2RhoPt[ij][i]);
521 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
a31b8a87 522 fh2PsiPt[ij][i] = new TH2F(Form("fh2PsiPt%s_j%d",cAdd.Data(),i),Form("jet shape #psi for %s jets;r;p_{T};",cAdd.Data()),
523 40,0.,2.,nBinPt,binLimitsPt);
524 fHistList->Add(fh2PsiPt[ij][i]);
37eb26ea 525 fh2PhiPt[ij][i] = new TH2F(Form("fh2PhiPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;#phi;p_{T};",cAdd.Data()),
f12de05f 526 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
7b822bc5 527
d95fc15a 528 fh3MultPtRP[ij][i] = new TH3F(Form("fh3MultPtRP%s_j%d",cAdd.Data(),i),Form("%s jets p_T;# tracks;;p_{T} (GeV/c)",cAdd.Data()),250,0,5000,nBinPt,0,binLimitsPt[nBinPt-1],(Int_t)fNRPBins,-0.5,fNRPBins-0.5);
7b822bc5 529 fHistList->Add(fh3MultPtRP[ij][i]);
530
531
532
f12de05f 533 fHistList->Add(fh2PhiPt[ij][i]);
37eb26ea 534 fh2EtaPt[ij][i] = new TH2F(Form("fh2EtaPt%s_j%d",cAdd.Data(),i),Form("pt vs eta %s;#eta;p_{T};",cAdd.Data()),
535 50,-1.,1.,nBinPt,binLimitsPt);
f12de05f 536 fHistList->Add(fh2EtaPt[ij][i]);
37eb26ea 537 fh2AreaPt[ij][i] = new TH2F(Form("fh2AreaPt%s_j%d",cAdd.Data(),i),
538 Form("pt vs area %s;area;p_{T};",
539 cAdd.Data()),
540 50,0.,1.,nBinPt,binLimitsPt);
541 fHistList->Add(fh2AreaPt[ij][i]);
542 fh2EtaArea[ij][i] = new TH2F(Form("fh2EtaArea%s_j%d",cAdd.Data(),i),
543 Form("area vs eta %s;#eta;area;",
544 cAdd.Data()),
545 50,-1.,1.,50,0,1.);
546 fHistList->Add(fh2EtaArea[ij][i]);
547
f12de05f 548 fh2PhiEta[ij][i] = new TH2F(Form("fh2PhiEta%s_j%d",cAdd.Data(),i),Form("phi vs eta %s ;phi;p_{T};",cAdd.Data()),
549 nBinPhi,binLimitsPhi,20,-1.,1.);
550 fHistList->Add(fh2PhiEta[ij][i]);
d7117cbd 551
552 fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
553 Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
554 cAdd.Data()),
555 200,0.,200.,nBinPt,binLimitsPt);
556 fHistList->Add(fh2LTrackPtJetPt[ij][i]);
a31b8a87 557 }
558
559
560 fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
561 fHistList->Add(fh1DijetMinv[ij]);
562
b5d90baf 563 fh2DijetDeltaPhiPt[ij] = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
a31b8a87 564 fHistList->Add(fh2DijetDeltaPhiPt[ij]);
565
566 fh2DijetAsymPt[ij] = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
567 fHistList->Add(fh2DijetAsymPt[ij]);
568
569 fh2DijetPt2vsPt1[ij] = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
570 fHistList->Add(fh2DijetPt2vsPt1[ij]);
a31b8a87 571 fh2DijetDifvsSum[ij] = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
572 fHistList->Add( fh2DijetDifvsSum[ij]);
573 }
3b7ffecf 574 // =========== Switch on Sumw2 for all histos ===========
575 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
576 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
577 if (h1){
578 h1->Sumw2();
579 continue;
580 }
581 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
582 if(hn)hn->Sumw2();
583 }
584 TH1::AddDirectory(oldStatus);
585}
586
587void AliAnalysisTaskJetSpectrum2::Init()
588{
589 //
590 // Initialization
591 //
592
3b7ffecf 593 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
594
595}
596
a31b8a87 597void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
598
b5cc0c6d 599
f2dd0695 600 Bool_t selected = kTRUE;
601
602 if(fUseGlobalSelection&&fEventSelectionMask==0){
603 selected = AliAnalysisHelperJetTasks::Selected();
604 }
605 else if(fUseGlobalSelection&&fEventSelectionMask>0){
45a11aef 606 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
f2dd0695 607 }
608
f4132e7d 609 if(fEventClass>0){
610 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
611 }
612
f2dd0695 613 if(!selected){
b5cc0c6d 614 // no selection by the service task, we continue
f4132e7d 615 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
b5cc0c6d 616 PostData(1, fHistList);
617 return;
618 }
619
f2dd0695 620
a31b8a87 621 static AliAODEvent* aod = 0;
622
623 // take all other information from the aod we take the tracks from
624 if(!aod){
625 if(fUseAODTrackInput)aod = fAODIn;
626 else aod = fAODOut;
627 }
628
629
3b7ffecf 630 //
631 // Execute analysis for current event
632 //
a31b8a87 633 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
634 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
635 if(!aodH){
636 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
637 return;
638 }
639
640 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
641 TClonesArray *aodRecJets = 0;
642
643 if(fAODOut&&!aodRecJets){
644 aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
645 }
646 if(fAODExtension&&!aodRecJets){
647 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
648 }
649 if(fAODIn&&!aodRecJets){
650 aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
651 }
652
653
654
655 if(!aodRecJets){
d95fc15a 656 if(fDebug){
657
658 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
659 if(fAODIn){
660 Printf("Input AOD >>>>");
661 fAODIn->Print();
662 }
663 if(fAODExtension){
664 Printf("AOD Extension >>>>");
665 fAODExtension->Print();
666 }
667 if(fAODOut){
668 Printf("Output AOD >>>>");
669 fAODOut->Print();
670 }
671 }
672 return;
a31b8a87 673 }
674
675 TClonesArray *aodGenJets = 0;
676 if(fBranchGen.Length()>0){
677 if(fAODOut&&!aodGenJets){
678 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
679 }
680 if(fAODExtension&&!aodGenJets){
681 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
682 }
683 if(fAODIn&&!aodGenJets){
684 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
685 }
686
687 if(!aodGenJets){
688 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
3b7ffecf 689 return;
690 }
3b7ffecf 691 }
a31b8a87 692
37eb26ea 693 TClonesArray *aodBackRecJets = 0;
694 if(fBranchBkgRec.Length()>0){
695 if(fAODOut&&!aodBackRecJets){
696 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
697 }
698 if(fAODExtension&&!aodBackRecJets){
699 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
700 }
701 if(fAODIn&&!aodBackRecJets){
702 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
703 }
704
705 if(!aodBackRecJets){
706 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
707 return;
708 }
709 }
710
711
712 TClonesArray *aodBackGenJets = 0;
713
714 if(fBranchBkgGen.Length()>0){
715 if(fAODOut&&!aodBackGenJets){
716 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
717 }
718 if(fAODExtension&&!aodBackGenJets){
719 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
720 }
721 if(fAODIn&&!aodBackGenJets){
722 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
723 }
724
725 if(!aodBackGenJets){
726 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
727 return;
728 }
729 }
730
a31b8a87 731
732 // new Scheme
733 // first fill all the pure histograms, i.e. full jets
734 // in case of the correaltion limit to the n-leading jets
735
736 // reconstructed
737
738
739 // generated
740
741
742 // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
743
744
742ee86c 745
a31b8a87 746 TList genJetsList; // full acceptance
747 TList genJetsListCut; // acceptance cut
748 TList recJetsList; // full acceptance
749 TList recJetsListCut; // acceptance cut
750
751
752 GetListOfJets(&recJetsList,aodRecJets,0);
753 GetListOfJets(&recJetsListCut,aodRecJets,1);
754
755 if(fBranchGen.Length()>0){
756 GetListOfJets(&genJetsList,aodGenJets,0);
757 GetListOfJets(&genJetsListCut,aodGenJets,1);
758 }
759
760 Double_t eventW = 1;
761 Double_t ptHard = 0;
762 Double_t nTrials = 1; // Trials for MC trigger
763 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
764
742ee86c 765 // Getting some global properties
766 fCentrality = GetCentrality();
767 fh1Centrality->Fill(fCentrality);
768
769
a31b8a87 770 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
771 // this is the part we only use when we have MC information
772 AliMCEvent* mcEvent = MCEvent();
773 // AliStack *pStack = 0;
774 if(!mcEvent){
775 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
3b7ffecf 776 return;
777 }
a31b8a87 778 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
779 if(pythiaGenHeader){
780 nTrials = pythiaGenHeader->Trials();
781 ptHard = pythiaGenHeader->GetPtHard();
782 int iProcessType = pythiaGenHeader->ProcessType();
783 // 11 f+f -> f+f
784 // 12 f+barf -> f+barf
785 // 13 f+barf -> g+g
786 // 28 f+g -> f+g
787 // 53 g+g -> f+barf
788 // 68 g+g -> g+g
789 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
790 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
791 }
792 }// (fAnalysisType&kMCESD)==kMCESD)
793
794
795 // we simply fetch the tracks/mc particles as a list of AliVParticles
796
797 TList recParticles;
798 TList genParticles;
799
800 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
801 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
802 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
803 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
804
742ee86c 805 CalculateReactionPlaneAngle(&recParticles);
806
a31b8a87 807 // Event control and counting ...
808 // MC
809 fh1PtHard->Fill(ptHard,eventW);
810 fh1PtHardNoW->Fill(ptHard,1);
811 fh1PtHardTrials->Fill(ptHard,nTrials);
812
813 // Real
784b1747 814 if(aod->GetPrimaryVertex()){// No vtx for pure MC
815 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
816 }
a31b8a87 817
a2522483 818
37eb26ea 819 Int_t recMult1 = recParticles.GetEntries();
820 Int_t genMult1 = genParticles.GetEntries();
821
742ee86c 822 Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
823 Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
37eb26ea 824
825 fh2MultRec->Fill(recMult1,recMult2);
826 fh2MultGen->Fill(genMult1,genMult2);
827 fMultRec = recMult1;
828 if(fMultRec<=0)fMultRec = recMult2;
829 fMultGen = genMult1;
830 if(fMultGen<=0)fMultGen = genMult2;
a2522483 831
a31b8a87 832 // the loops for rec and gen should be indentical... pass it to a separate
833 // function ...
834 // Jet Loop
835 // Track Loop
836 // Jet Jet Loop
837 // Jet Track loop
838
839 FillJetHistos(recJetsListCut,recParticles,kJetRec);
a2522483 840 FillJetHistos(recJetsList,recParticles,kJetRecFull);
a31b8a87 841 FillTrackHistos(recParticles,kJetRec);
842
843 FillJetHistos(genJetsListCut,genParticles,kJetGen);
a2522483 844 FillJetHistos(genJetsList,genParticles,kJetGenFull);
a31b8a87 845 FillTrackHistos(genParticles,kJetGen);
846
847 // Here follows the jet matching part
848 // delegate to a separated method?
849
850 FillMatchHistos(recJetsList,genJetsList);
d95fc15a 851
a31b8a87 852 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
853 PostData(1, fHistList);
d95fc15a 854 }
a31b8a87 855
856void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
857
858 if(iType>=kJetTypes){
859 return;
3b7ffecf 860 }
742ee86c 861 if(!fFlagJetType[iType])return;
a31b8a87 862
37eb26ea 863 Int_t refMult = fMultRec;
864 if(iType==kJetGen||iType==kJetGenFull){
865 refMult = fMultGen;
866 }
867
a31b8a87 868 Int_t nJets = jetsList.GetEntries();
869 fh1NJets[iType]->Fill(nJets);
870
871 if(nJets<=0)return;
3b7ffecf 872
a31b8a87 873 Float_t ptOld = 1.E+32;
874 Float_t pT0 = 0;
875 Float_t pT1 = 0;
876 Float_t phi0 = 0;
877 Float_t phi1 = 0;
878 Int_t ij0 = -1;
879 Int_t ij1 = -1;
880
881
882 for(int ij = 0;ij < nJets;ij++){
883 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
884 Float_t ptJet = jet->Pt();
885 fh1PtJetsIn[iType]->Fill(ptJet);
886 if(ptJet>ptOld){
887 Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
888 }
889 ptOld = ptJet;
890
b5d90baf 891 // find the dijets assume sorting and acceptance cut...
a31b8a87 892 if(ij==0){
893 ij0 = ij;
894 pT0 = ptJet;
895 phi0 = jet->Phi();
896 if(phi0<0)phi0+=TMath::Pi()*2.;
897 }
898 else if(ptJet>pT1){
899 // take only the backward hemisphere??
900 phi1 = jet->Phi();
b5d90baf 901 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 902 Float_t dPhi = phi1 - phi0;
903 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
904 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
905 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
906 ij1 = ij;
907 pT1 = ptJet;
908 }
909 }
910 // fill jet histos for kmax jets
37eb26ea 911
a31b8a87 912 Float_t phiJet = jet->Phi();
f12de05f 913 Float_t etaJet = jet->Eta();
a31b8a87 914 if(phiJet<0)phiJet+=TMath::Pi()*2.;
915 fh1TmpRho->Reset();
37eb26ea 916 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
917 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
a31b8a87 918 // fill leading jets...
d7117cbd 919 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
920 // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
742ee86c 921 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
37eb26ea 922 if(ptJet>10){
923 if(ij<kMaxJets){
924 fh2PhiEta[iType][ij]->Fill(phiJet,etaJet);
925 fh2AreaPt[iType][ij]->Fill(jet->EffectiveAreaCharged(),ptJet);
926 fh2EtaArea[iType][ij]->Fill(etaJet,jet->EffectiveAreaCharged());
a31b8a87 927 }
37eb26ea 928 fh2PhiEta[iType][kMaxJets]->Fill(phiJet,etaJet);
929 fh2AreaPt[iType][kMaxJets]->Fill(jet->EffectiveAreaCharged(),ptJet);
930 fh2EtaArea[iType][kMaxJets]->Fill(etaJet,jet->EffectiveAreaCharged());
931 }
932 if(ij<kMaxJets){
742ee86c 933 fh3MultPtRP[iType][ij]->Fill(refMult,ptJet,phiBin);
37eb26ea 934 fh2PhiPt[iType][ij]->Fill(phiJet,ptJet);
935 fh2EtaPt[iType][ij]->Fill(etaJet,ptJet);
d7117cbd 936 if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
a31b8a87 937 }
742ee86c 938 fh3MultPtRP[iType][kMaxJets]->Fill(refMult,ptJet,phiBin);
37eb26ea 939 fh2PhiPt[iType][kMaxJets]->Fill(phiJet,ptJet);
940 fh2EtaPt[iType][kMaxJets]->Fill(etaJet,ptJet);
d7117cbd 941 if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
a31b8a87 942
37eb26ea 943 if(particlesList.GetSize()&&ij<kMaxJets){
a31b8a87 944 // Particles... correlated with jets...
945 for(int it = 0;it<particlesList.GetEntries();++it){
946 AliVParticle *part = (AliVParticle*)particlesList.At(it);
947 Float_t deltaR = jet->DeltaR(part);
948 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
949 }
950 // fill the jet shapes
951 Float_t rhoSum = 0;
952 for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
953 Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
954 Float_t rho = fh1TmpRho->GetBinContent(ibx);
955 rhoSum += rho;
956 fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
957 fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
958 }
959 }// if we have particles
a31b8a87 960 }// Jet Loop
961
962
963 // do something with dijets...
964 if(ij0>=0&&ij1>0){
965 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
966 Double_t ptJet0 = jet0->Pt();
967 Double_t phiJet0 = jet0->Phi();
968 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
969
970 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
971 Double_t ptJet1 = jet1->Pt();
972 Double_t phiJet1 = jet1->Phi();
973 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
974
975 Float_t deltaPhi = phiJet0 - phiJet1;
976 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
977 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
978 deltaPhi = TMath::Abs(deltaPhi);
b5d90baf 979 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 980
981 Float_t asym = 9999;
982 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
983 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
984 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
985 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
986 Float_t minv = 2.*(jet0->P()*jet1->P()-
987 jet0->Px()*jet1->Px()-
988 jet0->Py()*jet1->Py()-
989 jet0->Pz()*jet1->Pz()); // assume mass == 0;
990 if(minv<0)minv=0; // prevent numerical instabilities
991 minv = TMath::Sqrt(minv);
992 fh1DijetMinv[iType]->Fill(minv);
993 }
994
995
996
997 // count down the jets above thrueshold
998 Int_t nOver = nJets;
999 if(nOver>0){
1000 TIterator *jetIter = jetsList.MakeIterator();
1001 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
a31b8a87 1002 if(tmpJet){
e855f5c5 1003 Float_t pt = tmpJet->Pt();
a31b8a87 1004 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1005 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1006 while(pt<ptCut&&tmpJet){
1007 nOver--;
1008 tmpJet = (AliAODJet*)(jetIter->Next());
1009 if(tmpJet){
1010 pt = tmpJet->Pt();
1011 }
1012 }
1013 if(nOver<=0)break;
1014 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1015 }
1016 }
1017 delete jetIter;
1018 }
1019}
1020
1021void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1022
742ee86c 1023 if(!fFlagJetType[iType])return;
7b822bc5 1024 Int_t refMult = fMultRec;
1025 if(iType==kJetGen||iType==kJetGenFull){
1026 refMult = fMultGen;
742ee86c 1027
7b822bc5 1028 }
1029
a31b8a87 1030 Int_t nTrackOver = particlesList.GetSize();
1031 // do the same for tracks and jets
1032 if(nTrackOver>0){
1033 TIterator *trackIter = particlesList.MakeIterator();
1034 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1035 Float_t pt = tmpTrack->Pt();
1036 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1037 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1038 while(pt<ptCut&&tmpTrack){
1039 nTrackOver--;
1040 tmpTrack = (AliVParticle*)(trackIter->Next());
1041 if(tmpTrack){
1042 pt = tmpTrack->Pt();
1043 }
1044 }
1045 if(nTrackOver<=0)break;
1046 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1047 }
1048
1049
1050 trackIter->Reset();
1051 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1052 Float_t sumPt = 0;
1053
1054 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1055 Float_t tmpPt = tmpTrack->Pt();
1056 fh1PtTracksIn[iType]->Fill(tmpPt);
cb883243 1057 fh1PtTracksInLow[iType]->Fill(tmpPt);
742ee86c 1058
a31b8a87 1059 sumPt += tmpPt;
1060 Float_t tmpPhi = tmpTrack->Phi();
1061 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
742ee86c 1062 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1063 fh3MultTrackPtRP[iType]->Fill(refMult,tmpPt,phiBin);
d95fc15a 1064 fh3MultTrackPtLowRP[iType]->Fill(refMult,tmpPt,phiBin);
3170a3f8 1065 fh2TrackEtaPt[iType]->Fill(tmpTrack->Eta(),tmpPt);
a31b8a87 1066 if(tmpTrack==leading){
1067 fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
1068 fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
1069 continue;
1070 }
1071 }
1072 fh1SumPtTrack[iType]->Fill(sumPt);
1073 delete trackIter;
1074 }
1075
1076}
1077
1078
1079void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1080
1081
1082 // Fill al the matching histos
1083 // It is important that the acceptances for the mathing are as large as possible
1084 // to avoid false matches on the edge of acceptance
1085 // therefore we add some extra matching jets as overhead
1086
1087 static TArrayI aGenIndex(recJetsList.GetEntries());
1088 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1089
a31b8a87 1090 static TArrayI aRecIndex(genJetsList.GetEntries());
1091 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1092
a31b8a87 1093 if(fDebug){
1094 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1095 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1096 }
1097 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1098 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1099 aGenIndex,aRecIndex,fDebug);
1100
1101 if(fDebug){
1102 for(int i = 0;i< aGenIndex.GetSize();++i){
1103 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1104 }
1105 for(int i = 0;i< aRecIndex.GetSize();++i){
1106 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1107 }
1108 }
1109
1110 Double_t container[6];
a31b8a87 1111
1112 // loop over generated jets
1113 // consider the
1114 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1115 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1116 Double_t ptGen = genJet->Pt();
1117 Double_t phiGen = genJet->Phi();
1118 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1119 Double_t etaGen = genJet->Eta();
1120 container[3] = ptGen;
1121 container[4] = etaGen;
1122 container[5] = phiGen;
1123 fhnJetContainer[kStep0]->Fill(&container[3]);
1124 if(JetSelected(genJet)){
1125 fhnJetContainer[kStep1]->Fill(&container[3]);
1126 Int_t ir = aRecIndex[ig];
1127 if(ir>=0&&ir<recJetsList.GetEntries()){
1128 fhnJetContainer[kStep2]->Fill(&container[3]);
1129 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1130 if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
1131 if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1132 }
1133 }
1134 }// loop over generated jets used for matching...
1135
1136
1137
1138 // loop over reconstructed jets
1139 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1140 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1141 Double_t etaRec = recJet->Eta();
1142 Double_t ptRec = recJet->Pt();
1143 Double_t phiRec = recJet->Phi();
1144 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1145 // do something with dijets...
1146
1147 container[0] = ptRec;
1148 container[1] = etaRec;
1149 container[2] = phiRec;
a31b8a87 1150
1151 fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1152 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1153
1154 if(JetSelected(recJet)){
1155 fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1156 // Fill Correlation
1157 Int_t ig = aGenIndex[ir];
1158 if(ig>=0 && ig<genJetsList.GetEntries()){
1159 fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1160 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1161 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1162 Double_t ptGen = genJet->Pt();
1163 Double_t phiGen = genJet->Phi();
1164 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1165 Double_t etaGen = genJet->Eta();
1166
1167 container[3] = ptGen;
1168 container[4] = etaGen;
1169 container[5] = phiGen;
a31b8a87 1170 //
1171 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1172 //
1173 if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1174 fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1175 fhnCorrelation->Fill(container);
1176 if(ptGen>0){
1177 Float_t delta = (ptRec-ptGen)/ptGen;
1178 fh2RelPtFGen->Fill(ptGen,delta);
1179 fh2PtFGen->Fill(ptGen,ptRec);
1180 }
a31b8a87 1181 }
a31b8a87 1182 }// loop over reconstructed jets
1183 }
1184 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1185}
1186
1187
3b7ffecf 1188void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1189 //
1190 // Create the particle container for the correction framework manager and
1191 // link it
1192 //
1193 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
d95fc15a 1194 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
3b7ffecf 1195 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1196 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1197
1198 // can we neglect migration in eta and phi?
1199 // phi should be no problem since we cover full phi and are phi symmetric
1200 // eta migration is more difficult due to needed acceptance correction
1201 // in limited eta range
1202
3b7ffecf 1203 //arrays for the number of bins in each dimension
1204 Int_t iBin[kNvar];
d95fc15a 1205 iBin[0] = 125; //bins in pt
3b7ffecf 1206 iBin[1] = 1; //bins in eta
1207 iBin[2] = 1; // bins in phi
1208
1209 //arrays for lower bounds :
1210 Double_t* binEdges[kNvar];
1211 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1212 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1213
1214 //values for bin lower bounds
1215 // for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
a923bd34 1216 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
3b7ffecf 1217 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1218 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1219
1220
1221 for(int i = 0;i < kMaxStep*2;++i){
f51451be 1222 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
3b7ffecf 1223 for (int k=0; k<kNvar; k++) {
1224 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1225 }
1226 }
1227 //create correlation matrix for unfolding
1228 Int_t thnDim[2*kNvar];
1229 for (int k=0; k<kNvar; k++) {
1230 //first half : reconstructed
1231 //second half : MC
1232 thnDim[k] = iBin[k];
1233 thnDim[k+kNvar] = iBin[k];
1234 }
1235
1236 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1237 for (int k=0; k<kNvar; k++) {
1238 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1239 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1240 }
1241 fhnCorrelation->Sumw2();
1242
742ee86c 1243
1244 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1245 delete [] binEdges[ivar];
1246
1247
3b7ffecf 1248}
1249
1250void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1251{
1252// Terminate analysis
1253//
1254 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1255}
1256
1257
1258Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1259
565584e8 1260 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1261
1262 Int_t iCount = 0;
565584e8 1263 if(type==kTrackAOD){
3b7ffecf 1264 AliAODEvent *aod = 0;
565584e8 1265 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1266 else aod = AODEvent();
3b7ffecf 1267 if(!aod){
1268 return iCount;
1269 }
1270 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1271 AliAODTrack *tr = aod->GetTrack(it);
1272 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 1273 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1274 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 1275 if(fDebug>0){
1276 if(tr->Pt()>20){
1277 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 1278 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1279 tr->Print();
1280 // tr->Dump();
1281 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1282 if(fESD){
a2522483 1283 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
38ecb6a5 1284 esdTr->Print("");
bb3d13aa 1285 // esdTr->Dump();
38ecb6a5 1286 }
1287 }
1288 }
3b7ffecf 1289 list->Add(tr);
1290 iCount++;
1291 }
1292 }
1293 else if (type == kTrackKineAll||type == kTrackKineCharged){
1294 AliMCEvent* mcEvent = MCEvent();
1295 if(!mcEvent)return iCount;
1296 // we want to have alivpartilces so use get track
1297 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1298 if(!mcEvent->IsPhysicalPrimary(it))continue;
1299 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 1300 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1301 if(type == kTrackKineAll){
1302 list->Add(part);
1303 iCount++;
1304 }
1305 else if(type == kTrackKineCharged){
1306 if(part->Particle()->GetPDG()->Charge()==0)continue;
1307 list->Add(part);
1308 iCount++;
1309 }
1310 }
1311 }
565584e8 1312 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1313 AliAODEvent *aod = 0;
5301f754 1314 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1315 else aod = AODEvent();
5010a3f7 1316 if(!aod)return iCount;
3b7ffecf 1317 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1318 if(!tca)return iCount;
1319 for(int it = 0;it < tca->GetEntriesFast();++it){
1320 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
e855f5c5 1321 if(!part)continue;
a31b8a87 1322 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1323 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1324 if(type == kTrackAODMCAll){
3b7ffecf 1325 list->Add(part);
1326 iCount++;
1327 }
565584e8 1328 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1329 if(part->Charge()==0)continue;
565584e8 1330 if(kTrackAODMCCharged){
1331 list->Add(part);
1332 }
1333 else {
a31b8a87 1334 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 1335 list->Add(part);
1336 }
3b7ffecf 1337 iCount++;
1338 }
1339 }
1340 }// AODMCparticle
cc0649e4 1341 list->Sort();
c4631cdb 1342 return iCount;
3b7ffecf 1343
1344}
a31b8a87 1345
1346
742ee86c 1347Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1348 AliAODEvent *aod = 0;
1349 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1350 else aod = AODEvent();
1351 if(!aod){
3170a3f8 1352 return 101;
742ee86c 1353 }
1354 return aod->GetHeader()->GetCentrality();
1355}
1356
1357
a31b8a87 1358
1359Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1360 Bool_t selected = false;
1361
1362 if(!jet)return selected;
1363
1364 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1365 selected = kTRUE;
1366 }
1367 return selected;
1368
1369}
1370
1371Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1372
1373 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1374 Int_t iCount = 0;
1375
1376 if(!jarray){
1377 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1378 return iCount;
1379 }
1380
1381
1382 for(int ij=0;ij<jarray->GetEntries();ij++){
1383 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1384 if(!jet)continue;
1385 if(type==0){
1386 // No cut at all, main purpose here is sorting
1387 list->Add(jet);
1388 iCount++;
1389 }
1390 else if(type == 1){
1391 // eta cut
1392 if(JetSelected(jet)){
1393 list->Add(jet);
1394 iCount++;
1395 }
1396 }
1397 }
1398
1399 list->Sort();
1400 return iCount;
1401
1402}
1403
1404
37eb26ea 1405Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1406 if(!jets)return 0;
1407
1408 Int_t refMult = 0;
1409 for(int ij = 0;ij < jets->GetEntries();++ij){
1410 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1411 if(!jet)continue;
1412 TRefArray *refs = jet->GetRefTracks();
1413 if(!refs)continue;
1414 refMult += refs->GetEntries();
1415 }
1416 return refMult;
1417
1418}
d7117cbd 1419
1420
1421AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1422 if(!jet)return 0;
1423 TRefArray *refs = jet->GetRefTracks();
1424 if(!refs) return 0;
1425 AliVParticle *leading = 0;
1426 Float_t fMaxPt = 0;
1427 for(int i = 0;i<refs->GetEntries();i++){
1428 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1429 if(!tmp)continue;
1430 if(tmp->Pt()>fMaxPt){
1431 leading = tmp;
1432 fMaxPt = tmp->Pt();
1433 }
1434 }
1435 return leading;
1436}
1437
1438
1439AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1440 if(!jet)return 0;
1441 if(!list) return 0;
1442 AliVParticle *leading = 0;
1443 Float_t fMaxPt = 0;
1444 for(int i = 0;i<list->GetEntries();i++){
1445 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1446 if(!tmp)continue;
1447 if(jet->DeltaR(tmp)>r)continue;
1448 if(tmp->Pt()>fMaxPt){
1449 leading = tmp;
1450 fMaxPt = tmp->Pt();
1451 }
1452 }
1453 return leading;
1454}
742ee86c 1455
1456Bool_t AliAnalysisTaskJetSpectrum2::CalculateReactionPlaneAngle(const TList *trackList)
1457{
1458
1459 if(!trackList)return kFALSE;
1460 fRPAngle=0;
1461
1462 // need to get this info from elsewhere??
742ee86c 1463
1464 Double_t fPsiRP =0,fDeltaPsiRP = 0;
1465
1466
1467
1468 TVector2 mQ,mQ1,mQ2;
c17f867b 1469 Float_t mQx= fDeltaQxy[0], mQy=fDeltaQxy[1];
742ee86c 1470
c17f867b 1471 Float_t mQx1=fDeltaQxy[0], mQy1=fDeltaQxy[1];
1472 Float_t mQx2=fDeltaQxy[0], mQy2=fDeltaQxy[1];
742ee86c 1473
1474 AliVParticle *track=0x0;
1475 Int_t count[3]={0,0,0};
1476
1477
1478 for (Int_t iter=0;iter<trackList->GetEntries();iter++){
1479
1480 track=(AliVParticle*)trackList->At(iter);
1481
1482 //cuts already applied before
1483 // Comment DCA not correctly implemented yet for AOD tracks
1484
1485 Double_t momentum;
1486 if(track->Charge()>0){momentum=track->Pt();}
1487 else{momentum=-track->Pt();}
1488
1489
1490
1491 // For Weighting
1492 fh3RPPhiTracks->Fill(fCentrality,momentum,track->Phi());
1493 count[0]++;
1494
c08c3ad2 1495 Double_t phiweight=GetPhiWeight(track->Phi(),momentum);
1496 // Double_t phiweight=1;
742ee86c 1497 Double_t weight=2;
1498 if(track->Pt()<2){weight=track->Pt();}
1499
1500
1501 mQx += (cos(2*track->Phi()))*weight*phiweight;
1502 mQy += (sin(2*track->Phi()))*weight*phiweight;
1503
1504 // Make random Subevents
1505
1506 if(fRPSubeventMethod==0){
1507 if(fRandomizer->Binomial(1,0.5)){
1508 mQx1 += (cos(2*track->Phi()))*weight*phiweight;
1509 mQy1 += (sin(2*track->Phi()))*weight*phiweight;
1510 count[1]++;}
1511 else{
1512 mQx2 += (cos(2*track->Phi()))*weight*phiweight;
1513 mQy2 += (sin(2*track->Phi()))*weight*phiweight;
1514 count[2]++;}
1515 }
1516 else if(fRPSubeventMethod==1){
1517 // Make eta dependent subevents
1518 if(track->Eta()>0){
1519 mQx1 += (cos(2*track->Phi()))*weight*phiweight;
1520 mQy1 += (sin(2*track->Phi()))*weight*phiweight;
1521 count[1]++;}
1522 else{
1523 mQx2 += (cos(2*track->Phi()))*weight*phiweight;
1524 mQy2 += (sin(2*track->Phi()))*weight*phiweight;
1525 count[2]++;}
1526 }
1527
1528 }
1529
1530
1531
1532 //If no track passes the cuts, the ,Q.Phi() will return Pi and a peak at Pi/2 in the RP Angular Distribution will appear
1533 if(count[0]==0||count[1]==0||count[2]==0){
1534 return kFALSE;
1535 }
1536
1537 mQ.Set(mQx,mQy);
1538 mQ1.Set(mQx1,mQy1);
1539 mQ2.Set(mQx2,mQy2);
1540
1541 // cout<<"MQ"<<mQx<<" " <<mQy<<" psi"<<endl;
1542
1543 fPsiRP=mQ.Phi()/2;
1544
1545 //Correction
1546 fPsiRP+=fFlatA[0]*TMath::Cos(2*fPsiRP)+fFlatB[0]*TMath::Sin(2*fPsiRP)+fFlatA[1]*TMath::Cos(4*fPsiRP)+fFlatB[1]*TMath::Sin(4*fPsiRP);
1547
1548 Double_t fPsiRP1=mQ1.Phi()/2;
1549 fPsiRP1+=fFlatA[0]*TMath::Cos(2*fPsiRP1)+fFlatB[0]*TMath::Sin(2*fPsiRP1)+fFlatA[1]*TMath::Cos(4*fPsiRP1)+fFlatB[1]*TMath::Sin(4*fPsiRP1);
1550 Double_t fPsiRP2=mQ2.Phi()/2;
1551 fPsiRP2+=fFlatA[0]*TMath::Cos(2*fPsiRP2)+fFlatB[0]*TMath::Sin(2*fPsiRP2)+fFlatA[1]*TMath::Cos(4*fPsiRP2)+fFlatB[1]*TMath::Sin(4*fPsiRP2);
1552 fDeltaPsiRP=fPsiRP1-fPsiRP2;
1553
1554 if(fPsiRP>TMath::Pi()){fPsiRP-=TMath::Pi();}
1555 if(fPsiRP<0){fPsiRP+=TMath::Pi();}
1556
1557 // reactionplaneangle + Pi() is the same angle
1558 if(TMath::Abs(fDeltaPsiRP)>TMath::Pi()/2){
1559 if(fDeltaPsiRP>0)fDeltaPsiRP-=TMath::Pi();
1560 else fDeltaPsiRP+=TMath::Pi();
1561 }
1562
1563 Double_t cos2deltaRP=TMath::Cos(2*fDeltaPsiRP);
1564
1565 // FillHistograms
1566 fh2RPSubevents->Fill(fPsiRP1,fPsiRP2);
1567 fh1RP->Fill(fPsiRP);
1568 fh2RPCentrality->Fill(fCentrality,fPsiRP);
1569 fh2RPDeltaRP->Fill(fDeltaPsiRP,fCentrality);
1570 fh2RPQxQy->Fill(mQx,mQy);
1571 fh2RPCosDeltaRP->Fill(fCentrality,cos2deltaRP);
1572
1573 fRPAngle=fPsiRP;
1574
1575 return kTRUE;
1576}
1577
1578Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1579{
1580 Int_t phibin=-1;
1581 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1582 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1583 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1584 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1585 return phibin;
1586}
1587
c08c3ad2 1588
1589
1590Double_t AliAnalysisTaskJetSpectrum2::GetPhiWeight(Double_t phi,Double_t signedpt){
1591 if(!fh3PhiWeights)return 1;
1592 else return fh3PhiWeights->GetBinContent(fh3PhiWeights->GetXaxis()->FindBin(fCentrality),fh3PhiWeights->GetYaxis()->FindBin(signedpt),fh3PhiWeights->GetZaxis()->FindBin(phi));
1593}
1594
742ee86c 1595 //________________________________________________________________________