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