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