enable tagging of jets with high qualitiy leading track
[u/mrichter/AliRoot.git] / PWGJE / 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>
2eded560 35#include <TProfile2D.h>
3b7ffecf 36#include <TList.h>
37#include <TLorentzVector.h>
38#include <TClonesArray.h>
37eb26ea 39#include <TRefArray.h>
3b7ffecf 40#include "TDatabasePDG.h"
41
42#include "AliAnalysisTaskJetSpectrum2.h"
43#include "AliAnalysisManager.h"
44#include "AliJetFinder.h"
b99018c6 45#include "AliTHn.h"
3b7ffecf 46#include "AliJetHeader.h"
47#include "AliJetReader.h"
48#include "AliJetReaderHeader.h"
49#include "AliUA1JetHeaderV1.h"
50#include "AliESDEvent.h"
51#include "AliAODEvent.h"
52#include "AliAODHandler.h"
53#include "AliAODTrack.h"
54#include "AliAODJet.h"
c2785065 55#include "AliAODJetEventBackground.h"
3b7ffecf 56#include "AliAODMCParticle.h"
57#include "AliMCEventHandler.h"
58#include "AliMCEvent.h"
59#include "AliStack.h"
60#include "AliGenPythiaEventHeader.h"
61#include "AliJetKineReaderHeader.h"
62#include "AliGenCocktailEventHeader.h"
63#include "AliInputEventHandler.h"
64
65
66#include "AliAnalysisHelperJetTasks.h"
67
68ClassImp(AliAnalysisTaskJetSpectrum2)
69
37eb26ea 70AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
71 AliAnalysisTaskSE(),
a31b8a87 72 fJetHeaderRec(0x0),
73 fJetHeaderGen(0x0),
74 fAODIn(0x0),
75 fAODOut(0x0),
76 fAODExtension(0x0),
b99018c6 77 fhnJetContainer(0x0),
a31b8a87 78 fhnCorrelation(0x0),
90a006c2 79 fhnEvent(0x0),
c8eabe24 80 f1PtScale(0x0),
a31b8a87 81 fBranchRec("jets"),
82 fBranchGen(""),
37eb26ea 83 fBranchBkgRec(""),
84 fBranchBkgGen(""),
a31b8a87 85 fNonStdFile(""),
742ee86c 86 fRandomizer(0x0),
a31b8a87 87 fUseAODJetInput(kFALSE),
88 fUseAODTrackInput(kFALSE),
89 fUseAODMCInput(kFALSE),
90 fUseGlobalSelection(kFALSE),
91 fUseExternalWeightOnly(kFALSE),
92 fLimitGenJetEta(kFALSE),
b99018c6 93 fDoMatching(kFALSE),
a31b8a87 94 fNMatchJets(5),
7b822bc5 95 fNRPBins(3),
65c43de8 96 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
d7396b04 97 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
a31b8a87 98 fFilterMask(0),
f2dd0695 99 fEventSelectionMask(0),
3d1dba29 100 fNTrigger(0),
101 fTriggerBit(0x0),
a31b8a87 102 fAnalysisType(0),
3b7ffecf 103 fTrackTypeRec(kTrackUndef),
104 fTrackTypeGen(kTrackUndef),
f4132e7d 105 fEventClass(0),
2eded560 106 fRPMethod(0),
3b7ffecf 107 fAvgTrials(1),
108 fExternalWeight(1),
a31b8a87 109 fJetRecEtaWindow(0.5),
110 fTrackRecEtaWindow(0.5),
9280dfa6 111 fMinJetPt(0),
a31b8a87 112 fMinTrackPt(0.15),
113 fDeltaPhiWindow(90./180.*TMath::Pi()),
742ee86c 114 fCentrality(100),
7b822bc5 115 fRPAngle(0),
37eb26ea 116 fMultRec(0),
117 fMultGen(0),
3d1dba29 118 fTriggerName(0x0),
3b7ffecf 119 fh1Xsec(0x0),
120 fh1Trials(0x0),
121 fh1PtHard(0x0),
122 fh1PtHardNoW(0x0),
123 fh1PtHardTrials(0x0),
57ca1193 124 fh1ZVtx(0x0),
7b822bc5 125 fh1RP(0x0),
742ee86c 126 fh1Centrality(0x0),
c8eabe24 127 fh1TmpRho(0x0),
37eb26ea 128 fh2MultRec(0x0),
129 fh2MultGen(0x0),
742ee86c 130 fh2RPCentrality(0x0),
a31b8a87 131 fh2PtFGen(0x0),
9a83d4af 132 fh2RelPtFGen(0x0),
3b7ffecf 133 fHistList(0x0)
134{
c17f867b 135
a31b8a87 136 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 137 fFlagJetType[ij] = 1; // default = on
a31b8a87 138 fh1NJets[ij] = 0;
139 fh1SumPtTrack[ij] = 0;
140 fh1PtJetsIn[ij] = 0;
65c43de8 141 fh1PtJetsInRej[ij] = 0;
a31b8a87 142 fh1PtTracksIn[ij] = 0;
cb883243 143 fh1PtTracksInLow[ij] = 0;
a31b8a87 144 fh2NJetsPt[ij] = 0;
145 fh2NTracksPt[ij] = 0;
2eded560 146 fp2MultRPPhiTrackPt[ij] = 0;
147 fp2CentRPPhiTrackPt[ij] = 0;
90a006c2 148 fhnJetPt[ij] = 0;
149 fhnJetPtQA[ij] = 0;
150 fhnTrackPt[ij] = 0;
151 fhnTrackPtQA[ij] = 0;
37eb26ea 152 for(int i = 0;i <= kMaxJets;++i){
d7117cbd 153 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 154 fh1PtIn[ij][i] = 0;
a31b8a87 155 }
156
157 fh1DijetMinv[ij] = 0;
158 fh2DijetDeltaPhiPt[ij] = 0;
159 fh2DijetAsymPt[ij] = 0;
160 fh2DijetPt2vsPt1[ij] = 0;
161 fh2DijetDifvsSum[ij] = 0;
162
163 }
3b7ffecf 164}
165
166AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
167 AliAnalysisTaskSE(name),
168 fJetHeaderRec(0x0),
169 fJetHeaderGen(0x0),
a31b8a87 170 fAODIn(0x0),
171 fAODOut(0x0),
172 fAODExtension(0x0),
b99018c6 173 fhnJetContainer(0x0),
3b7ffecf 174 fhnCorrelation(0x0),
90a006c2 175 fhnEvent(0x0),
c8eabe24 176 f1PtScale(0x0),
3b7ffecf 177 fBranchRec("jets"),
178 fBranchGen(""),
37eb26ea 179 fBranchBkgRec(""),
180 fBranchBkgGen(""),
a31b8a87 181 fNonStdFile(""),
742ee86c 182 fRandomizer(0x0),
565584e8 183 fUseAODJetInput(kFALSE),
184 fUseAODTrackInput(kFALSE),
185 fUseAODMCInput(kFALSE),
b5cc0c6d 186 fUseGlobalSelection(kFALSE),
3b7ffecf 187 fUseExternalWeightOnly(kFALSE),
188 fLimitGenJetEta(kFALSE),
b99018c6 189 fDoMatching(kFALSE),
a31b8a87 190 fNMatchJets(5),
7b822bc5 191 fNRPBins(3),
65c43de8 192 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
d7396b04 193 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
3b7ffecf 194 fFilterMask(0),
f2dd0695 195 fEventSelectionMask(0),
3d1dba29 196 fNTrigger(0),
197 fTriggerBit(0x0),
3b7ffecf 198 fAnalysisType(0),
199 fTrackTypeRec(kTrackUndef),
200 fTrackTypeGen(kTrackUndef),
f4132e7d 201 fEventClass(0),
2eded560 202 fRPMethod(0),
3b7ffecf 203 fAvgTrials(1),
204 fExternalWeight(1),
a31b8a87 205 fJetRecEtaWindow(0.5),
206 fTrackRecEtaWindow(0.5),
9280dfa6 207 fMinJetPt(0),
a31b8a87 208 fMinTrackPt(0.15),
209 fDeltaPhiWindow(90./180.*TMath::Pi()),
742ee86c 210 fCentrality(100),
7b822bc5 211 fRPAngle(0),
37eb26ea 212 fMultRec(0),
213 fMultGen(0),
3d1dba29 214 fTriggerName(0x0),
3b7ffecf 215 fh1Xsec(0x0),
216 fh1Trials(0x0),
217 fh1PtHard(0x0),
218 fh1PtHardNoW(0x0),
219 fh1PtHardTrials(0x0),
57ca1193 220 fh1ZVtx(0x0),
7b822bc5 221 fh1RP(0x0),
742ee86c 222 fh1Centrality(0x0),
c8eabe24 223 fh1TmpRho(0x0),
37eb26ea 224 fh2MultRec(0x0),
225 fh2MultGen(0x0),
742ee86c 226 fh2RPCentrality(0x0),
a31b8a87 227 fh2PtFGen(0x0),
9a83d4af 228 fh2RelPtFGen(0x0),
3b7ffecf 229 fHistList(0x0)
230{
231
a31b8a87 232 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 233 fFlagJetType[ij] = 1; // default = on
a31b8a87 234 fh1NJets[ij] = 0;
235 fh1SumPtTrack[ij] = 0;
236 fh1PtJetsIn[ij] = 0;
65c43de8 237 fh1PtJetsInRej[ij] = 0;
a31b8a87 238 fh1PtTracksIn[ij] = 0;
cb883243 239 fh1PtTracksInLow[ij] = 0;
a31b8a87 240 fh2NJetsPt[ij] = 0;
241 fh2NTracksPt[ij] = 0;
2eded560 242 fp2MultRPPhiTrackPt[ij] = 0;
243 fp2CentRPPhiTrackPt[ij] = 0;
90a006c2 244 fhnJetPt[ij] = 0;
245 fhnJetPtQA[ij] = 0;
246 fhnTrackPt[ij] = 0;
247 fhnTrackPtQA[ij] = 0;
37eb26ea 248 for(int i = 0;i <= kMaxJets;++i){
d7117cbd 249 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 250 fh1PtIn[ij][i] = 0;
a31b8a87 251 }
252
253 fh1DijetMinv[ij] = 0;
254 fh2DijetDeltaPhiPt[ij] = 0;
255 fh2DijetAsymPt[ij] = 0;
256 fh2DijetPt2vsPt1[ij] = 0;
257 fh2DijetDifvsSum[ij] = 0;
258 }
259
37eb26ea 260 DefineOutput(1, TList::Class());
3b7ffecf 261}
262
263
264
265Bool_t AliAnalysisTaskJetSpectrum2::Notify()
266{
a31b8a87 267
268
269
3b7ffecf 270 //
271 // Implemented Notify() to read the cross sections
272 // and number of trials from pyxsec.root
273 //
a31b8a87 274
275 // Fetch the aod also from the input in,
276 // have todo it in notify
277
278
279 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
280 // assume that the AOD is in the general output...
281 fAODOut = AODEvent();
282
283 if(fNonStdFile.Length()!=0){
284 // case that we have an AOD extension we need can fetch the jets from the extended output
285 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
e855f5c5 286 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
a31b8a87 287 if(!fAODExtension){
288 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
289 }
290 }
291
3b7ffecf 292
293 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
294 Float_t xsection = 0;
295 Float_t ftrials = 1;
296
297 fAvgTrials = 1;
298 if(tree){
299 TFile *curfile = tree->GetCurrentFile();
300 if (!curfile) {
301 Error("Notify","No current file");
302 return kFALSE;
303 }
304 if(!fh1Xsec||!fh1Trials){
305 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
306 return kFALSE;
307 }
308 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
309 fh1Xsec->Fill("<#sigma>",xsection);
310 // construct a poor man average trials
311 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 312 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
3b7ffecf 313 }
a31b8a87 314
d7117cbd 315 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
a31b8a87 316
3b7ffecf 317 return kTRUE;
318}
319
320void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
321{
322
742ee86c 323
3b7ffecf 324 // Connect the AOD
325
3b7ffecf 326 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
3b7ffecf 327 OpenFile(1);
37eb26ea 328 if(!fHistList)fHistList = new TList();
40440b28 329 PostData(1, fHistList); // post data in any case once
330
742ee86c 331 if(!fRandomizer)fRandomizer = new TRandom3(0);
332
5cb0f438 333 fHistList->SetOwner(kTRUE);
37eb26ea 334 Bool_t oldStatus = TH1::AddDirectoryStatus();
3b7ffecf 335 TH1::AddDirectory(kFALSE);
336
90a006c2 337
338
339 // event npsparse cent, mult
3d1dba29 340 const Int_t nBinsSparse0 = 3;
341 const Int_t nBins0[nBinsSparse0] = { 100, 500,fNTrigger};
342 const Double_t xmin0[nBinsSparse0] = { 0, 0, -0.5};
343 const Double_t xmax0[nBinsSparse0] = { 100,5000,fNTrigger-0.5};
90a006c2 344
345
346 fhnEvent = new THnSparseF("fhnEvent",";cent;mult",nBinsSparse0,
347 nBins0,xmin0,xmax0);
348 fHistList->Add(fhnEvent);
349
b99018c6 350 if(fDoMatching){
351 MakeJetContainer();
352 fHistList->Add(fhnCorrelation);
353 fHistList->Add(fhnJetContainer);
354 }
3b7ffecf 355 //
356 // Histogram
357
90a006c2 358
359
ed5c9f3f 360 const Int_t nBinPt = 120;
3b7ffecf 361 Double_t binLimitsPt[nBinPt+1];
362 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
363 if(iPt == 0){
ed5c9f3f 364 binLimitsPt[iPt] = -50.0;
3b7ffecf 365 }
366 else {// 1.0
e24d63c0 367 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
3b7ffecf 368 }
369 }
edfbe476 370 const Int_t nBinPhi = 90;
3b7ffecf 371 Double_t binLimitsPhi[nBinPhi+1];
372 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
373 if(iPhi==0){
37eb26ea 374 binLimitsPhi[iPhi] = 0;
3b7ffecf 375 }
376 else{
377 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
378 }
379 }
380
edfbe476 381 const Int_t nBinEta = 40;
382 Double_t binLimitsEta[nBinEta+1];
383 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
384 if(iEta==0){
385 binLimitsEta[iEta] = -2.0;
386 }
387 else{
cf049e94 388 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 389 }
390 }
391
a31b8a87 392
3b7ffecf 393 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
394 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
90a006c2 395 fHistList->Add(fh1Xsec);
3b7ffecf 396 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
397 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
a31b8a87 398 fHistList->Add(fh1Trials);
3b7ffecf 399 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 400 fHistList->Add(fh1PtHard);
3b7ffecf 401 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 402 fHistList->Add(fh1PtHardNoW);
3b7ffecf 403 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 404 fHistList->Add(fh1PtHardTrials);
a31b8a87 405
57ca1193 406 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
a31b8a87 407 fHistList->Add(fh1ZVtx);
7b822bc5 408
2eded560 409
9939142f 410 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
7b822bc5 411 fHistList->Add(fh1RP);
412
eb1739b7 413 fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
742ee86c 414 fHistList->Add(fh1Centrality);
415
3170a3f8 416 fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
37eb26ea 417 fHistList->Add(fh2MultRec);
3170a3f8 418 fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
37eb26ea 419 fHistList->Add(fh2MultGen);
a31b8a87 420
a31b8a87 421
742ee86c 422 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
423 fHistList->Add(fh2RPCentrality);
a31b8a87 424
742ee86c 425 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
426 fHistList->Add(fh2PtFGen);
a31b8a87 427
742ee86c 428 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
429 fHistList->Add(fh2RelPtFGen);
742ee86c 430
431 for(int ij = 0;ij <kJetTypes;++ij){
432 TString cAdd = "";
433 TString cJetBranch = "";
434 if(ij==kJetRec){
435 cAdd = "Rec";
436 cJetBranch = fBranchRec.Data();
437 }
438 else if (ij==kJetGen){
439 cAdd = "Gen";
440 cJetBranch = fBranchGen.Data();
441 }
442 else if (ij==kJetRecFull){
443 cAdd = "RecFull";
444 cJetBranch = fBranchRec.Data();
445 }
446 else if (ij==kJetGenFull){
447 cAdd = "GenFull";
448 cJetBranch = fBranchGen.Data();
449 }
a31b8a87 450
742ee86c 451 if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
452 if(!fFlagJetType[ij])continue;
7b822bc5 453
742ee86c 454 fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
455 fHistList->Add(fh1NJets[ij]);
456
457 fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
458 fHistList->Add(fh1PtJetsIn[ij]);
65c43de8 459
460 fh1PtJetsInRej[ij] = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
461 fHistList->Add(fh1PtJetsInRej[ij]);
d7396b04 462
463 fh1PtJetsInBest[ij] = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
464 fHistList->Add(fh1PtJetsInBest[ij]);
742ee86c 465
466 fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
467 fHistList->Add(fh1PtTracksIn[ij]);
468
2eded560 469 fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
742ee86c 470 fHistList->Add(fh1PtTracksInLow[ij]);
471
3170a3f8 472 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 473 fHistList->Add(fh1SumPtTrack[ij]);
474
475 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);
476 fHistList->Add(fh2NJetsPt[ij]);
477
3170a3f8 478 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 479 fHistList->Add(fh2NTracksPt[ij]);
90a006c2 480
2eded560 481
482 fp2MultRPPhiTrackPt[ij] = new TProfile2D(Form("fp2MultRPPhiTrackPt%s",cAdd.Data()),"RP phi vs Number of tracks;# tracks;#Delta#phi_{RP}; <p_{T}>",20,0,4000,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
483 fHistList->Add(fp2MultRPPhiTrackPt[ij]);
484 fp2CentRPPhiTrackPt[ij] = new TProfile2D(Form("fp2CentRPPhiTrackPt%s",cAdd.Data()),"RP phi vs cent;# cent;#Delta#phi_{RP}; <p_{T}>",10,0,100,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
485 fHistList->Add(fp2CentRPPhiTrackPt[ij]);
486
3d1dba29 487 // Bins: Jet number: pTJet, cent, mult, RP, Area, trigger total bins = 4.5M
488 const Int_t nBinsSparse1 = 7;
489 Int_t nBins1[nBinsSparse1] = { kMaxJets+1,120, 10, 25, fNRPBins, 10,fNTrigger};
d01931c3 490 if(cJetBranch.Contains("RandomCone")){
491 nBins1[1] = 600;
492 nBins1[5] = 1;
493 }
3d1dba29 494 const Double_t xmin1[nBinsSparse1] = { -0.5,-50, 0, 0, -0.5, 0.,-0.5};
495 const Double_t xmax1[nBinsSparse1] = {kMaxJets+0.5,250,100,5000,fNRPBins-0.5,1.0,fNTrigger-0.5};
2eded560 496
3d1dba29 497 fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger",nBinsSparse1,nBins1,xmin1,xmax1);
2eded560 498 fHistList->Add(fhnJetPt[ij]);
499
500 // Bins: Jet number: pTJet, cent, eta, phi, Area. total bins = 9.72 M
501 const Int_t nBinsSparse2 = 6;
d01931c3 502 Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 25, 5, 18, 360, 10};
503 if(cJetBranch.Contains("RandomCone")){
504 nBins2[5] = 1;
505 }
506 const Double_t xmin2[nBinsSparse2] = { -0.5, 0, 0,-0.9, 0, 0.};
507 const Double_t xmax2[nBinsSparse2] = {kMaxJets+0.5,250, 100, 0.9, 2.*TMath::Pi(),1.0};
508 fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area",nBinsSparse2,nBins2,xmin2,xmax2);
509 fHistList->Add(fhnJetPtQA[ij]);
510
511 // Bins:track number pTtrack, cent, mult, RP. total bins = 224 k
512 const Int_t nBinsSparse3 = 5;
513 const Int_t nBins3[nBinsSparse3] = { 2, 100, 10, 20, fNRPBins};
514 const Double_t xmin3[nBinsSparse3] = { -0.5, 0, 0, 0, -0.5};
515 const Double_t xmax3[nBinsSparse3] = { 1.5, 200, 100, 4000,fNRPBins-0.5};
516
90a006c2 517 // change the binning ot the pT axis:
d01931c3 518 Double_t *xPt3 = new Double_t[nBins3[1]+1];
519 xPt3[0] = 0.;
520 for(int i = 1; i<=nBins3[1];i++){
521 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
522 else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
523 else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
524 else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 1.; // 62 - 72
2eded560 525 else if(xPt3[i-1]<30)xPt3[i] = xPt3[i-1] + 2.5; // 74 - 78
526 else xPt3[i] = xPt3[i-1] + 5.; // 78 - 100 = 140
d01931c3 527 }
528
90a006c2 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;
30f922ad 536 const Int_t nBins4[nBinsSparse4] = { 2, 50, 10, 20, 360};
537 const Double_t xmin4[nBinsSparse4] = { -0.5, 0, 0, -1.0, 0.};
538 const Double_t xmax4[nBinsSparse4] = { 1.5,150, 100, 1.0,2.*TMath::Pi()};
90a006c2 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
d3a3f33d 815 // CalculateReactionPlaneAngle(&recParticles);
2eded560 816 fRPAngle = 0;
817
818 if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane();
819 else if(fRPMethod==1||fRPMethod==2){
820 fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod);
821 }
822 fh1RP->Fill(fRPAngle);
823 fh2RPCentrality->Fill(fCentrality,fRPAngle);
a31b8a87 824 // Event control and counting ...
825 // MC
826 fh1PtHard->Fill(ptHard,eventW);
827 fh1PtHardNoW->Fill(ptHard,1);
828 fh1PtHardTrials->Fill(ptHard,nTrials);
829
830 // Real
784b1747 831 if(aod->GetPrimaryVertex()){// No vtx for pure MC
832 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
833 }
a31b8a87 834
a2522483 835
37eb26ea 836 Int_t recMult1 = recParticles.GetEntries();
837 Int_t genMult1 = genParticles.GetEntries();
838
742ee86c 839 Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
840 Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
37eb26ea 841
842 fh2MultRec->Fill(recMult1,recMult2);
843 fh2MultGen->Fill(genMult1,genMult2);
844 fMultRec = recMult1;
845 if(fMultRec<=0)fMultRec = recMult2;
846 fMultGen = genMult1;
847 if(fMultGen<=0)fMultGen = genMult2;
a2522483 848
441afcb9 849 Double_t var0[3] = {0,};
90a006c2 850 var0[0] = fCentrality;
851 var0[1] = fMultRec;
3d1dba29 852 for(int it=0;it<fNTrigger;it++){
853 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
854 var0[2] = it;
855 fhnEvent->Fill(var0);
856 }
857 }
a31b8a87 858 // the loops for rec and gen should be indentical... pass it to a separate
859 // function ...
860 // Jet Loop
861 // Track Loop
862 // Jet Jet Loop
863 // Jet Track loop
864
865 FillJetHistos(recJetsListCut,recParticles,kJetRec);
a2522483 866 FillJetHistos(recJetsList,recParticles,kJetRecFull);
a31b8a87 867 FillTrackHistos(recParticles,kJetRec);
868
869 FillJetHistos(genJetsListCut,genParticles,kJetGen);
a2522483 870 FillJetHistos(genJetsList,genParticles,kJetGenFull);
a31b8a87 871 FillTrackHistos(genParticles,kJetGen);
872
873 // Here follows the jet matching part
874 // delegate to a separated method?
875
b99018c6 876 if(fDoMatching){
877 FillMatchHistos(recJetsList,genJetsList);
878 }
879
a31b8a87 880 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
881 PostData(1, fHistList);
90a006c2 882}
a31b8a87 883
884void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
885
886 if(iType>=kJetTypes){
887 return;
3b7ffecf 888 }
742ee86c 889 if(!fFlagJetType[iType])return;
a31b8a87 890
37eb26ea 891 Int_t refMult = fMultRec;
892 if(iType==kJetGen||iType==kJetGenFull){
893 refMult = fMultGen;
894 }
895
a31b8a87 896 Int_t nJets = jetsList.GetEntries();
897 fh1NJets[iType]->Fill(nJets);
898
899 if(nJets<=0)return;
3b7ffecf 900
a31b8a87 901 Float_t ptOld = 1.E+32;
902 Float_t pT0 = 0;
903 Float_t pT1 = 0;
904 Float_t phi0 = 0;
905 Float_t phi1 = 0;
906 Int_t ij0 = -1;
907 Int_t ij1 = -1;
908
441afcb9 909 Double_t var1[7] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area
90a006c2 910 var1[2] = fCentrality;
911 var1[3] = refMult;
912
913 Double_t var2[6] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area
914 var2[2] = fCentrality;
a31b8a87 915
916 for(int ij = 0;ij < nJets;ij++){
917 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
918 Float_t ptJet = jet->Pt();
4a828dd2 919 if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
d7396b04 920 if(jet->Trigger()&fJetTriggerBestMask){
921 fh1PtJetsInBest[iType]->Fill(ptJet);
922 }
65c43de8 923 if(jet->Trigger()&fJetTriggerExcludeMask){
924 fh1PtJetsInRej[iType]->Fill(ptJet);
925 continue;
926 }
a31b8a87 927 fh1PtJetsIn[iType]->Fill(ptJet);
a31b8a87 928 ptOld = ptJet;
929
b5d90baf 930 // find the dijets assume sorting and acceptance cut...
a31b8a87 931 if(ij==0){
932 ij0 = ij;
933 pT0 = ptJet;
934 phi0 = jet->Phi();
935 if(phi0<0)phi0+=TMath::Pi()*2.;
936 }
937 else if(ptJet>pT1){
938 // take only the backward hemisphere??
939 phi1 = jet->Phi();
b5d90baf 940 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 941 Float_t dPhi = phi1 - phi0;
942 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
943 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
944 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
945 ij1 = ij;
946 pT1 = ptJet;
947 }
948 }
949 // fill jet histos for kmax jets
37eb26ea 950
a31b8a87 951 Float_t phiJet = jet->Phi();
f12de05f 952 Float_t etaJet = jet->Eta();
a31b8a87 953 if(phiJet<0)phiJet+=TMath::Pi()*2.;
954 fh1TmpRho->Reset();
37eb26ea 955 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
90a006c2 956
37eb26ea 957 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
a31b8a87 958 // fill leading jets...
d7117cbd 959 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
960 // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
742ee86c 961 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
90a006c2 962
963 var1[1] = ptJet;
964 var1[4] = phiBin;
965 var1[5] = jet->EffectiveAreaCharged();
966
967 var2[1] = ptJet;
968 var2[3] = etaJet;
969 var2[4] = phiJet;
970 var2[5] = jet->EffectiveAreaCharged();
37eb26ea 971 if(ij<kMaxJets){
d7117cbd 972 if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
90a006c2 973 var1[0] = ij;
974 var2[0] = ij;
3d1dba29 975 for(int it = 0;it <fNTrigger;it++){
976 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
977 var1[6] = it;
978 fhnJetPt[iType]->Fill(var1);
979 }
980 }
90a006c2 981 fhnJetPtQA[iType]->Fill(var2);
a31b8a87 982 }
90a006c2 983 var1[0] = kMaxJets;// fill for all jets
984 var2[0] = kMaxJets;// fill for all jets
3d1dba29 985 for(int it = 0;it <fNTrigger;it++){
986 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
987 var1[6] = it;
988 fhnJetPt[iType]->Fill(var1);
989 }
990 }
991
90a006c2 992 fhnJetPtQA[iType]->Fill(var2);
d7117cbd 993 if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
a31b8a87 994
37eb26ea 995 if(particlesList.GetSize()&&ij<kMaxJets){
a31b8a87 996 // Particles... correlated with jets...
997 for(int it = 0;it<particlesList.GetEntries();++it){
998 AliVParticle *part = (AliVParticle*)particlesList.At(it);
999 Float_t deltaR = jet->DeltaR(part);
1000 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1001 }
1002 // fill the jet shapes
a31b8a87 1003 }// if we have particles
a31b8a87 1004 }// Jet Loop
1005
1006
1007 // do something with dijets...
1008 if(ij0>=0&&ij1>0){
1009 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1010 Double_t ptJet0 = jet0->Pt();
1011 Double_t phiJet0 = jet0->Phi();
1012 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1013
1014 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1015 Double_t ptJet1 = jet1->Pt();
1016 Double_t phiJet1 = jet1->Phi();
1017 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1018
1019 Float_t deltaPhi = phiJet0 - phiJet1;
1020 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1021 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1022 deltaPhi = TMath::Abs(deltaPhi);
b5d90baf 1023 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 1024
1025 Float_t asym = 9999;
1026 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1027 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1028 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1029 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1030 Float_t minv = 2.*(jet0->P()*jet1->P()-
1031 jet0->Px()*jet1->Px()-
1032 jet0->Py()*jet1->Py()-
1033 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1034 if(minv<0)minv=0; // prevent numerical instabilities
1035 minv = TMath::Sqrt(minv);
1036 fh1DijetMinv[iType]->Fill(minv);
1037 }
1038
1039
1040
1041 // count down the jets above thrueshold
1042 Int_t nOver = nJets;
1043 if(nOver>0){
1044 TIterator *jetIter = jetsList.MakeIterator();
1045 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
a31b8a87 1046 if(tmpJet){
e855f5c5 1047 Float_t pt = tmpJet->Pt();
a31b8a87 1048 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1049 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1050 while(pt<ptCut&&tmpJet){
1051 nOver--;
1052 tmpJet = (AliAODJet*)(jetIter->Next());
1053 if(tmpJet){
1054 pt = tmpJet->Pt();
1055 }
1056 }
1057 if(nOver<=0)break;
1058 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1059 }
1060 }
1061 delete jetIter;
1062 }
1063}
1064
1065void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1066
90a006c2 1067 if(fFlagJetType[iType]<=0)return;
7b822bc5 1068 Int_t refMult = fMultRec;
1069 if(iType==kJetGen||iType==kJetGenFull){
1070 refMult = fMultGen;
742ee86c 1071
7b822bc5 1072 }
1073
90a006c2 1074 //
1075 Double_t var3[5]; // track number;p_{T};cent;#tracks;RP
1076 var3[2] = fCentrality;
1077 var3[3] = refMult;
1078 Double_t var4[5]; // track number;p_{T};cent;#eta;#phi
1079 var4[2] = fCentrality;
a31b8a87 1080 Int_t nTrackOver = particlesList.GetSize();
1081 // do the same for tracks and jets
1082 if(nTrackOver>0){
1083 TIterator *trackIter = particlesList.MakeIterator();
1084 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1085 Float_t pt = tmpTrack->Pt();
1086 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1087 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1088 while(pt<ptCut&&tmpTrack){
1089 nTrackOver--;
1090 tmpTrack = (AliVParticle*)(trackIter->Next());
1091 if(tmpTrack){
1092 pt = tmpTrack->Pt();
1093 }
1094 }
1095 if(nTrackOver<=0)break;
1096 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1097 }
1098
1099
1100 trackIter->Reset();
1101 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1102 Float_t sumPt = 0;
1103
1104 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1105 Float_t tmpPt = tmpTrack->Pt();
1106 fh1PtTracksIn[iType]->Fill(tmpPt);
cb883243 1107 fh1PtTracksInLow[iType]->Fill(tmpPt);
742ee86c 1108
a31b8a87 1109 sumPt += tmpPt;
1110 Float_t tmpPhi = tmpTrack->Phi();
1111 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
2eded560 1112
1113
1114 Float_t phiRP = tmpPhi-fRPAngle;
1115 if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1116 if(phiRP<0)phiRP += TMath::Pi();
1117 if(phiRP<0)phiRP += TMath::Pi();
1118 const float allPhi = -1./180.*TMath::Pi();
1119
dae7dd67 1120 if(tmpPt<100){
1121 fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1122 fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1123
1124 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1125 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1126 }
742ee86c 1127 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
90a006c2 1128 var3[0] = 1;
1129 var3[1] = tmpPt;
1130 var3[4] = phiBin;
1131
1132 var4[0] = 1;
1133 var4[1] = tmpPt;
1134 var4[3] = tmpTrack->Eta();
1135 var4[4] = tmpPhi;
1136
1137
1138 fhnTrackPt[iType]->Fill(var3);
1139 fhnTrackPtQA[iType]->Fill(var4);
1140
a31b8a87 1141 if(tmpTrack==leading){
90a006c2 1142 var3[0] = 0;
1143 var4[0] = 0;
1144 fhnTrackPt[iType]->Fill(var3);
1145 fhnTrackPtQA[iType]->Fill(var4);
a31b8a87 1146 continue;
1147 }
1148 }
1149 fh1SumPtTrack[iType]->Fill(sumPt);
1150 delete trackIter;
1151 }
1152
1153}
1154
1155
1156void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1157
1158
1159 // Fill al the matching histos
1160 // It is important that the acceptances for the mathing are as large as possible
1161 // to avoid false matches on the edge of acceptance
1162 // therefore we add some extra matching jets as overhead
1163
1164 static TArrayI aGenIndex(recJetsList.GetEntries());
1165 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1166
a31b8a87 1167 static TArrayI aRecIndex(genJetsList.GetEntries());
1168 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1169
a31b8a87 1170 if(fDebug){
1171 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1172 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1173 }
1174 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1175 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1176 aGenIndex,aRecIndex,fDebug);
1177
1178 if(fDebug){
1179 for(int i = 0;i< aGenIndex.GetSize();++i){
1180 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1181 }
1182 for(int i = 0;i< aRecIndex.GetSize();++i){
1183 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1184 }
1185 }
1186
1187 Double_t container[6];
a31b8a87 1188
1189 // loop over generated jets
1190 // consider the
1191 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1192 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1193 Double_t ptGen = genJet->Pt();
1194 Double_t phiGen = genJet->Phi();
1195 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1196 Double_t etaGen = genJet->Eta();
1197 container[3] = ptGen;
1198 container[4] = etaGen;
1199 container[5] = phiGen;
b99018c6 1200 fhnJetContainer->Fill(&container[3],kStep0);
a31b8a87 1201 if(JetSelected(genJet)){
b99018c6 1202 fhnJetContainer->Fill(&container[3],kStep1);
a31b8a87 1203 Int_t ir = aRecIndex[ig];
1204 if(ir>=0&&ir<recJetsList.GetEntries()){
b99018c6 1205 fhnJetContainer->Fill(&container[3],kStep2);
a31b8a87 1206 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
b99018c6 1207 if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4);
1208 if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3);
a31b8a87 1209 }
1210 }
1211 }// loop over generated jets used for matching...
1212
1213
1214
1215 // loop over reconstructed jets
1216 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1217 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1218 Double_t etaRec = recJet->Eta();
1219 Double_t ptRec = recJet->Pt();
1220 Double_t phiRec = recJet->Phi();
1221 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1222 // do something with dijets...
1223
1224 container[0] = ptRec;
1225 container[1] = etaRec;
1226 container[2] = phiRec;
a31b8a87 1227
b99018c6 1228 fhnJetContainer->Fill(container,kStep0+kMaxStep);
a31b8a87 1229 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1230
1231 if(JetSelected(recJet)){
b99018c6 1232 fhnJetContainer->Fill(container,kStep1+kMaxStep);
a31b8a87 1233 // Fill Correlation
1234 Int_t ig = aGenIndex[ir];
1235 if(ig>=0 && ig<genJetsList.GetEntries()){
b99018c6 1236 fhnJetContainer->Fill(container,kStep2+kMaxStep);
a31b8a87 1237 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1238 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1239 Double_t ptGen = genJet->Pt();
1240 Double_t phiGen = genJet->Phi();
1241 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1242 Double_t etaGen = genJet->Eta();
1243
1244 container[3] = ptGen;
1245 container[4] = etaGen;
1246 container[5] = phiGen;
a31b8a87 1247 //
1248 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1249 //
b99018c6 1250 if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep);
1251 fhnJetContainer->Fill(container,kStep3+kMaxStep);
1252 fhnCorrelation->Fill(container,0);
a31b8a87 1253 if(ptGen>0){
1254 Float_t delta = (ptRec-ptGen)/ptGen;
1255 fh2RelPtFGen->Fill(ptGen,delta);
1256 fh2PtFGen->Fill(ptGen,ptRec);
1257 }
a31b8a87 1258 }
a31b8a87 1259 }// loop over reconstructed jets
1260 }
1261 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1262}
1263
1264
3b7ffecf 1265void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1266 //
1267 // Create the particle container for the correction framework manager and
1268 // link it
1269 //
1270 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
d95fc15a 1271 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
3b7ffecf 1272 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1273 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1274
1275 // can we neglect migration in eta and phi?
1276 // phi should be no problem since we cover full phi and are phi symmetric
1277 // eta migration is more difficult due to needed acceptance correction
1278 // in limited eta range
1279
3b7ffecf 1280 //arrays for the number of bins in each dimension
1281 Int_t iBin[kNvar];
d95fc15a 1282 iBin[0] = 125; //bins in pt
3b7ffecf 1283 iBin[1] = 1; //bins in eta
1284 iBin[2] = 1; // bins in phi
1285
1286 //arrays for lower bounds :
1287 Double_t* binEdges[kNvar];
1288 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1289 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1290
1291 //values for bin lower bounds
1292 // 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 1293 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 1294 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1295 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1296
1297
b99018c6 1298 fhnJetContainer = new AliTHn(Form("fahnJetContainer"),Form("AliTHn jet info"),kMaxStep*2,kNvar,iBin);
1299 for (int k=0; k<kNvar; k++) {
1300 fhnJetContainer->SetBinLimits(k,binEdges[k]);
3b7ffecf 1301 }
b99018c6 1302
3b7ffecf 1303 //create correlation matrix for unfolding
1304 Int_t thnDim[2*kNvar];
1305 for (int k=0; k<kNvar; k++) {
1306 //first half : reconstructed
1307 //second half : MC
1308 thnDim[k] = iBin[k];
1309 thnDim[k+kNvar] = iBin[k];
1310 }
1311
b99018c6 1312 fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim);
3b7ffecf 1313 for (int k=0; k<kNvar; k++) {
b99018c6 1314 fhnCorrelation->SetBinLimits(k,binEdges[k]);
1315 fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]);
3b7ffecf 1316 }
742ee86c 1317
1318 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1319 delete [] binEdges[ivar];
1320
1321
3b7ffecf 1322}
1323
1324void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1325{
1326// Terminate analysis
1327//
1328 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1329}
1330
1331
1332Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1333
565584e8 1334 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1335
1336 Int_t iCount = 0;
565584e8 1337 if(type==kTrackAOD){
3b7ffecf 1338 AliAODEvent *aod = 0;
565584e8 1339 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1340 else aod = AODEvent();
3b7ffecf 1341 if(!aod){
1342 return iCount;
1343 }
1344 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1345 AliAODTrack *tr = aod->GetTrack(it);
1346 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 1347 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1348 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 1349 if(fDebug>0){
1350 if(tr->Pt()>20){
1351 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 1352 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1353 tr->Print();
1354 // tr->Dump();
1355 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1356 if(fESD){
a2522483 1357 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
38ecb6a5 1358 esdTr->Print("");
bb3d13aa 1359 // esdTr->Dump();
38ecb6a5 1360 }
1361 }
1362 }
3b7ffecf 1363 list->Add(tr);
1364 iCount++;
1365 }
1366 }
1367 else if (type == kTrackKineAll||type == kTrackKineCharged){
1368 AliMCEvent* mcEvent = MCEvent();
1369 if(!mcEvent)return iCount;
1370 // we want to have alivpartilces so use get track
1371 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1372 if(!mcEvent->IsPhysicalPrimary(it))continue;
1373 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 1374 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1375 if(type == kTrackKineAll){
1376 list->Add(part);
1377 iCount++;
1378 }
1379 else if(type == kTrackKineCharged){
1380 if(part->Particle()->GetPDG()->Charge()==0)continue;
1381 list->Add(part);
1382 iCount++;
1383 }
1384 }
1385 }
565584e8 1386 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1387 AliAODEvent *aod = 0;
5301f754 1388 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1389 else aod = AODEvent();
5010a3f7 1390 if(!aod)return iCount;
3b7ffecf 1391 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1392 if(!tca)return iCount;
1393 for(int it = 0;it < tca->GetEntriesFast();++it){
1394 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
e855f5c5 1395 if(!part)continue;
a31b8a87 1396 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1397 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1398 if(type == kTrackAODMCAll){
3b7ffecf 1399 list->Add(part);
1400 iCount++;
1401 }
565584e8 1402 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1403 if(part->Charge()==0)continue;
565584e8 1404 if(kTrackAODMCCharged){
1405 list->Add(part);
1406 }
1407 else {
a31b8a87 1408 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 1409 list->Add(part);
1410 }
3b7ffecf 1411 iCount++;
1412 }
1413 }
1414 }// AODMCparticle
cc0649e4 1415 list->Sort();
c4631cdb 1416 return iCount;
3b7ffecf 1417
1418}
a31b8a87 1419
1420
742ee86c 1421Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1422 AliAODEvent *aod = 0;
1423 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1424 else aod = AODEvent();
1425 if(!aod){
3170a3f8 1426 return 101;
742ee86c 1427 }
1428 return aod->GetHeader()->GetCentrality();
1429}
1430
1431
a31b8a87 1432
1433Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1434 Bool_t selected = false;
1435
1436 if(!jet)return selected;
1437
1438 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1439 selected = kTRUE;
1440 }
1441 return selected;
1442
1443}
1444
1445Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1446
1447 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1448 Int_t iCount = 0;
1449
1450 if(!jarray){
1451 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1452 return iCount;
1453 }
1454
1455
1456 for(int ij=0;ij<jarray->GetEntries();ij++){
1457 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1458 if(!jet)continue;
1459 if(type==0){
1460 // No cut at all, main purpose here is sorting
1461 list->Add(jet);
1462 iCount++;
1463 }
1464 else if(type == 1){
1465 // eta cut
1466 if(JetSelected(jet)){
1467 list->Add(jet);
1468 iCount++;
1469 }
1470 }
1471 }
1472
1473 list->Sort();
1474 return iCount;
1475
1476}
1477
1478
37eb26ea 1479Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1480 if(!jets)return 0;
1481
1482 Int_t refMult = 0;
1483 for(int ij = 0;ij < jets->GetEntries();++ij){
1484 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1485 if(!jet)continue;
1486 TRefArray *refs = jet->GetRefTracks();
1487 if(!refs)continue;
1488 refMult += refs->GetEntries();
1489 }
1490 return refMult;
1491
1492}
d7117cbd 1493
1494
1495AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1496 if(!jet)return 0;
1497 TRefArray *refs = jet->GetRefTracks();
1498 if(!refs) return 0;
1499 AliVParticle *leading = 0;
1500 Float_t fMaxPt = 0;
1501 for(int i = 0;i<refs->GetEntries();i++){
1502 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1503 if(!tmp)continue;
1504 if(tmp->Pt()>fMaxPt){
1505 leading = tmp;
1506 fMaxPt = tmp->Pt();
1507 }
1508 }
1509 return leading;
1510}
1511
1512
1513AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1514 if(!jet)return 0;
1515 if(!list) return 0;
1516 AliVParticle *leading = 0;
1517 Float_t fMaxPt = 0;
1518 for(int i = 0;i<list->GetEntries();i++){
1519 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1520 if(!tmp)continue;
1521 if(jet->DeltaR(tmp)>r)continue;
1522 if(tmp->Pt()>fMaxPt){
1523 leading = tmp;
1524 fMaxPt = tmp->Pt();
1525 }
1526 }
1527 return leading;
1528}
742ee86c 1529
742ee86c 1530Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1531{
1532 Int_t phibin=-1;
1533 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1534 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1535 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1536 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1537 return phibin;
1538}
1539
3d1dba29 1540void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1541 if(n>0){
1542 fNTrigger = n;
1543 delete [] fTriggerName;
1544 fTriggerName = new TString [fNTrigger];
1545 delete [] fTriggerBit;fTriggerBit = 0;
1546 fTriggerBit = new UInt_t [fNTrigger];
1547 }
1548 else{
1549 fNTrigger = 0;
1550 }
1551}
1552
1553void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1554 if(i<fNTrigger){
1555 fTriggerBit[i] = it;
1556 fTriggerName[i] = c;
1557 }
1558}
1559
1560AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1561 //
1562 delete [] fTriggerBit;
1563 delete [] fTriggerName;
1564}