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