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