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