]> 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
update from megan
[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){
346cfe4a 327 Printf("%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
e5ee1c19 363 const Int_t nBinsSparse0 = 4;
364 const Int_t nBins0[nBinsSparse0] = { 100, 500,fNTrigger,125};
365 const Double_t xmin0[nBinsSparse0] = { 0, 0, -0.5,-2};
366 const Double_t xmax0[nBinsSparse0] = { 100,5000,fNTrigger-0.5,248};
90a006c2 367
368
e5ee1c19 369 fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0,
90a006c2 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;
5ff21917 548 Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 60, 8, 18, 72, 10,fNTrigger,fNAcceptance+0.5,20};
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
e5ee1c19 898 Double_t var0[4] = {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;
e5ee1c19 904 var0[3] = GetRho(recJetsList);
3d1dba29 905 fhnEvent->Fill(var0);
906 }
907 }
a31b8a87 908 // the loops for rec and gen should be indentical... pass it to a separate
909 // function ...
910 // Jet Loop
911 // Track Loop
912 // Jet Jet Loop
913 // Jet Track loop
914
915 FillJetHistos(recJetsListCut,recParticles,kJetRec);
a2522483 916 FillJetHistos(recJetsList,recParticles,kJetRecFull);
a31b8a87 917 FillTrackHistos(recParticles,kJetRec);
918
919 FillJetHistos(genJetsListCut,genParticles,kJetGen);
a2522483 920 FillJetHistos(genJetsList,genParticles,kJetGenFull);
a31b8a87 921 FillTrackHistos(genParticles,kJetGen);
922
923 // Here follows the jet matching part
924 // delegate to a separated method?
925
b99018c6 926 if(fDoMatching){
927 FillMatchHistos(recJetsList,genJetsList);
928 }
929
a31b8a87 930 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
931 PostData(1, fHistList);
90a006c2 932}
a31b8a87 933
934void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
935
936 if(iType>=kJetTypes){
937 return;
3b7ffecf 938 }
742ee86c 939 if(!fFlagJetType[iType])return;
a31b8a87 940
37eb26ea 941 Int_t refMult = fMultRec;
942 if(iType==kJetGen||iType==kJetGenFull){
943 refMult = fMultGen;
944 }
945
a31b8a87 946 Int_t nJets = jetsList.GetEntries();
947 fh1NJets[iType]->Fill(nJets);
948
949 if(nJets<=0)return;
3b7ffecf 950
a31b8a87 951 Float_t ptOld = 1.E+32;
952 Float_t pT0 = 0;
953 Float_t pT1 = 0;
954 Float_t phi0 = 0;
955 Float_t phi1 = 0;
956 Int_t ij0 = -1;
957 Int_t ij1 = -1;
958
db745d29 959 Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
17f4943e 960 Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
7bee4353 961
90a006c2 962 var1[2] = fCentrality;
17f4943e 963 var1b[1] = fCentrality;
90a006c2 964 var1[3] = refMult;
965
17f4943e 966
967
5662640e 968 Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
90a006c2 969 var2[2] = fCentrality;
a31b8a87 970
971 for(int ij = 0;ij < nJets;ij++){
972 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
973 Float_t ptJet = jet->Pt();
17f4943e 974
975
4a828dd2 976 if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
17f4943e 977
978 var1b[0] = ptJet;
d7396b04 979 if(jet->Trigger()&fJetTriggerBestMask){
980 fh1PtJetsInBest[iType]->Fill(ptJet);
17f4943e 981 for(int it = 0;it <fNTrigger;it++){
982 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
983 var1b[2] = it;
984 fhnJetPtBest[iType]->Fill(var1b);
985 }
986 }
d7396b04 987 }
65c43de8 988 if(jet->Trigger()&fJetTriggerExcludeMask){
989 fh1PtJetsInRej[iType]->Fill(ptJet);
17f4943e 990 for(int it = 0;it <fNTrigger;it++){
991 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
992 var1b[2] = it;
993 fhnJetPtRej[iType]->Fill(var1b);
994 }
995 }
65c43de8 996 continue;
997 }
a31b8a87 998 fh1PtJetsIn[iType]->Fill(ptJet);
a31b8a87 999 ptOld = ptJet;
1000
b5d90baf 1001 // find the dijets assume sorting and acceptance cut...
a31b8a87 1002 if(ij==0){
1003 ij0 = ij;
1004 pT0 = ptJet;
1005 phi0 = jet->Phi();
1006 if(phi0<0)phi0+=TMath::Pi()*2.;
1007 }
1008 else if(ptJet>pT1){
1009 // take only the backward hemisphere??
1010 phi1 = jet->Phi();
b5d90baf 1011 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 1012 Float_t dPhi = phi1 - phi0;
1013 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1014 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1015 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1016 ij1 = ij;
1017 pT1 = ptJet;
1018 }
1019 }
1020 // fill jet histos for kmax jets
37eb26ea 1021
adf3920d 1022 Float_t phiJet = jet->Phi();
1023 Float_t etaJet = jet->Eta();
1024 if(phiJet<0)phiJet+=TMath::Pi()*2.;
1025 fh1TmpRho->Reset();
1026 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1027
1028 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1029 // fill leading jets...
1030 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1031 // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1032 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1033 Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1034
1035 var1[1] = ptJet;
1036 var1[4] = phiBin;
1037 var1[5] = jet->EffectiveAreaCharged();
1038 var1[7] = ptLead;
1039 var1[8] = CheckAcceptance(phiJet,etaJet);
1040
1041 var2[1] = ptJet;
1042 var2[3] = etaJet;
1043 var2[4] = phiJet;
1044 var2[5] = jet->EffectiveAreaCharged();
1045 var2[7] = var1[8];
1046 var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1047
1048 if(ij<kMaxJets){
1049 fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1050 var1[0] = ij;
1051 var2[0] = ij;
3d1dba29 1052 for(int it = 0;it <fNTrigger;it++){
1053 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1054 var1[6] = it;
adf3920d 1055 var2[6] = it;
3d1dba29 1056 fhnJetPt[iType]->Fill(var1);
1d478650 1057 fhnJetPtQA[iType]->Fill(var2);
3d1dba29 1058 }
1059 }
adf3920d 1060 }
1061 var1[0] = kMaxJets;// fill for all jets
1062 var2[0] = kMaxJets;// fill for all jets
1063 for(int it = 0;it <fNTrigger;it++){
1064 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1065 var1[6] = it;
94c3d5e4 1066 var2[6] = it;
adf3920d 1067 fhnJetPt[iType]->Fill(var1);
1068 fhnJetPtQA[iType]->Fill(var2);
1069 }
1070 }
3d1dba29 1071
adf3920d 1072 fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
a31b8a87 1073
adf3920d 1074 if(particlesList.GetSize()&&ij<kMaxJets){
1075 // Particles... correlated with jets...
1076 for(int it = 0;it<particlesList.GetEntries();++it){
1077 AliVParticle *part = (AliVParticle*)particlesList.At(it);
1078 Float_t deltaR = jet->DeltaR(part);
1079 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1080 }
1081 // fill the jet shapes
1082 }// if we have particles
a31b8a87 1083 }// Jet Loop
1084
1085
1086 // do something with dijets...
1087 if(ij0>=0&&ij1>0){
1088 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1089 Double_t ptJet0 = jet0->Pt();
1090 Double_t phiJet0 = jet0->Phi();
1091 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1092
1093 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1094 Double_t ptJet1 = jet1->Pt();
1095 Double_t phiJet1 = jet1->Phi();
1096 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1097
1098 Float_t deltaPhi = phiJet0 - phiJet1;
1099 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1100 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1101 deltaPhi = TMath::Abs(deltaPhi);
b5d90baf 1102 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 1103
1104 Float_t asym = 9999;
1105 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
adf3920d 1106 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1107 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1108 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1109 Float_t minv = 2.*(jet0->P()*jet1->P()-
1110 jet0->Px()*jet1->Px()-
1111 jet0->Py()*jet1->Py()-
1112 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1113 if(minv<0)minv=0; // prevent numerical instabilities
1114 minv = TMath::Sqrt(minv);
1115 fh1DijetMinv[iType]->Fill(minv);
a31b8a87 1116 }
1117
1118
1119
1120 // count down the jets above thrueshold
1121 Int_t nOver = nJets;
1122 if(nOver>0){
1123 TIterator *jetIter = jetsList.MakeIterator();
1124 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
a31b8a87 1125 if(tmpJet){
e855f5c5 1126 Float_t pt = tmpJet->Pt();
a31b8a87 1127 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1128 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1129 while(pt<ptCut&&tmpJet){
1130 nOver--;
1131 tmpJet = (AliAODJet*)(jetIter->Next());
1132 if(tmpJet){
1133 pt = tmpJet->Pt();
1134 }
1135 }
1136 if(nOver<=0)break;
1137 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1138 }
1139 }
1140 delete jetIter;
1141 }
1142}
1143
1144void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1145
90a006c2 1146 if(fFlagJetType[iType]<=0)return;
7b822bc5 1147 Int_t refMult = fMultRec;
1148 if(iType==kJetGen||iType==kJetGenFull){
1149 refMult = fMultGen;
742ee86c 1150
7b822bc5 1151 }
1152
90a006c2 1153 //
17f4943e 1154 Double_t var3[6]; // track number;p_{T};cent;#tracks;RP
90a006c2 1155 var3[2] = fCentrality;
1156 var3[3] = refMult;
5662640e 1157 Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
90a006c2 1158 var4[2] = fCentrality;
a31b8a87 1159 Int_t nTrackOver = particlesList.GetSize();
1160 // do the same for tracks and jets
1161 if(nTrackOver>0){
1162 TIterator *trackIter = particlesList.MakeIterator();
1163 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1164 Float_t pt = tmpTrack->Pt();
1165 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1166 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1167 while(pt<ptCut&&tmpTrack){
1168 nTrackOver--;
1169 tmpTrack = (AliVParticle*)(trackIter->Next());
1170 if(tmpTrack){
1171 pt = tmpTrack->Pt();
1172 }
1173 }
1174 if(nTrackOver<=0)break;
1175 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1176 }
1177
1178
1179 trackIter->Reset();
1180 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1181 Float_t sumPt = 0;
1182
1183 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1184 Float_t tmpPt = tmpTrack->Pt();
1185 fh1PtTracksIn[iType]->Fill(tmpPt);
cb883243 1186 fh1PtTracksInLow[iType]->Fill(tmpPt);
742ee86c 1187
a31b8a87 1188 sumPt += tmpPt;
1189 Float_t tmpPhi = tmpTrack->Phi();
1190 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
2eded560 1191
1192
1193 Float_t phiRP = tmpPhi-fRPAngle;
1194 if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1195 if(phiRP<0)phiRP += TMath::Pi();
1196 if(phiRP<0)phiRP += TMath::Pi();
1197 const float allPhi = -1./180.*TMath::Pi();
1198
dae7dd67 1199 if(tmpPt<100){
1200 fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1201 fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1202
1203 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1204 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1205 }
742ee86c 1206 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
90a006c2 1207 var3[0] = 1;
1208 var3[1] = tmpPt;
1209 var3[4] = phiBin;
1210
1211 var4[0] = 1;
1212 var4[1] = tmpPt;
1213 var4[3] = tmpTrack->Eta();
1214 var4[4] = tmpPhi;
5662640e 1215 var4[5] = tmpTrack->Charge();
90a006c2 1216
17f4943e 1217
1218 for(int it = 0;it <fNTrigger;it++){
1219 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1220 var3[0] = 1;
1221 var3[5] = it;
1222 fhnTrackPt[iType]->Fill(var3);
1223 if(tmpTrack==leading){
1224 var3[0] = 0;
1225 fhnTrackPt[iType]->Fill(var3);
1226 }
1227 }
1228 }
1229
1230
1231
1232
1233
90a006c2 1234 fhnTrackPtQA[iType]->Fill(var4);
1235
a31b8a87 1236 if(tmpTrack==leading){
90a006c2 1237 var4[0] = 0;
90a006c2 1238 fhnTrackPtQA[iType]->Fill(var4);
a31b8a87 1239 continue;
1240 }
1241 }
1242 fh1SumPtTrack[iType]->Fill(sumPt);
1243 delete trackIter;
1244 }
1245
1246}
1247
1248
1249void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1250
1251
1252 // Fill al the matching histos
1253 // It is important that the acceptances for the mathing are as large as possible
1254 // to avoid false matches on the edge of acceptance
1255 // therefore we add some extra matching jets as overhead
1256
1257 static TArrayI aGenIndex(recJetsList.GetEntries());
1258 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1259
a31b8a87 1260 static TArrayI aRecIndex(genJetsList.GetEntries());
1261 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1262
a31b8a87 1263 if(fDebug){
1264 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1265 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1266 }
1267 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1268 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1269 aGenIndex,aRecIndex,fDebug);
1270
1271 if(fDebug){
1272 for(int i = 0;i< aGenIndex.GetSize();++i){
1273 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1274 }
1275 for(int i = 0;i< aRecIndex.GetSize();++i){
1276 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1277 }
1278 }
1279
152a0744 1280 Double_t container[8];
1281 Double_t containerRec[4];
1282 Double_t containerGen[4];
a31b8a87 1283
1284 // loop over generated jets
1285 // consider the
1286 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1287 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1288 Double_t ptGen = genJet->Pt();
1289 Double_t phiGen = genJet->Phi();
1290 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1291 Double_t etaGen = genJet->Eta();
152a0744 1292 containerGen[0] = ptGen;
1293 containerGen[1] = etaGen;
1294 containerGen[2] = phiGen;
1295 containerGen[3] = genJet->GetPtLeading();
1296
1297 fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1298
1299 Int_t ir = aRecIndex[ig];
1300 if(ir>=0&&ir<recJetsList.GetEntries()){
1301 fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1302
1303 if(JetSelected(genJet)){
1304 fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1305
a31b8a87 1306 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
152a0744 1307
1308 fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1309 if(JetSelected(recJet)) {
1310 fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1311
1312 containerGen[3] = recJet->GetPtLeading();
1313 fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner in eta window with leading track on reconstructed level
1314 }
a31b8a87 1315 }
1316 }
1317 }// loop over generated jets used for matching...
1318
1319
1320
1321 // loop over reconstructed jets
1322 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1323 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1324 Double_t etaRec = recJet->Eta();
1325 Double_t ptRec = recJet->Pt();
1326 Double_t phiRec = recJet->Phi();
1327 if(phiRec<0)phiRec+=TMath::Pi()*2.;
a31b8a87 1328
152a0744 1329 containerRec[0] = ptRec;
1330 containerRec[1] = etaRec;
1331 containerRec[2] = phiRec;
1332 containerRec[3] = recJet->GetPtLeading();
1333
a31b8a87 1334 container[0] = ptRec;
1335 container[1] = etaRec;
1336 container[2] = phiRec;
152a0744 1337 container[3] = recJet->GetPtLeading();
a31b8a87 1338
a31b8a87 1339 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1340
1341 if(JetSelected(recJet)){
152a0744 1342 fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
a31b8a87 1343 // Fill Correlation
1344 Int_t ig = aGenIndex[ir];
1345 if(ig>=0 && ig<genJetsList.GetEntries()){
a31b8a87 1346 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1347 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1348 Double_t ptGen = genJet->Pt();
1349 Double_t phiGen = genJet->Phi();
1350 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1351 Double_t etaGen = genJet->Eta();
152a0744 1352
1353 containerGen[0] = ptGen;
1354 containerGen[1] = etaGen;
1355 containerGen[2] = phiGen;
1356 containerGen[3] = genJet->GetPtLeading();
1357
1358 container[4] = ptGen;
1359 container[5] = etaGen;
1360 container[6] = phiGen;
1361 container[7] = genJet->GetPtLeading();
1362
1363 fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1364
1365 fhnCorrelation->Fill(container);
a31b8a87 1366 if(ptGen>0){
1367 Float_t delta = (ptRec-ptGen)/ptGen;
1368 fh2RelPtFGen->Fill(ptGen,delta);
1369 fh2PtFGen->Fill(ptGen,ptRec);
1370 }
a31b8a87 1371 }
a31b8a87 1372 }// loop over reconstructed jets
1373 }
1374 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1375}
1376
1377
3b7ffecf 1378void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1379 //
1380 // Create the particle container for the correction framework manager and
1381 // link it
1382 //
152a0744 1383 const Int_t kNvar = 4 ; //number of variables on the grid:pt,eta, phi
d95fc15a 1384 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
152a0744 1385 const Double_t kEtamin = -1.0, kEtamax = 1.0;
3b7ffecf 1386 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
152a0744 1387 const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
3b7ffecf 1388
1389 // can we neglect migration in eta and phi?
1390 // phi should be no problem since we cover full phi and are phi symmetric
1391 // eta migration is more difficult due to needed acceptance correction
1392 // in limited eta range
1393
3b7ffecf 1394 //arrays for the number of bins in each dimension
1395 Int_t iBin[kNvar];
d95fc15a 1396 iBin[0] = 125; //bins in pt
152a0744 1397 iBin[1] = 4; //bins in eta
1398 iBin[2] = 1; // bins in phi
1399 iBin[3] = 10; // bins in leading track Pt
3b7ffecf 1400
1401 //arrays for lower bounds :
1402 Double_t* binEdges[kNvar];
1403 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1404 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1405
1406 //values for bin lower bounds
1407 // 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 1408 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 1409 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1410 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
152a0744 1411 binEdges[3][0]= kPtLeadingTrackPtMin;
1412 binEdges[3][1]= 1.;
1413 binEdges[3][2]= 2.;
1414 binEdges[3][3]= 3.;
1415 binEdges[3][4]= 4.;
1416 binEdges[3][5]= 5.;
1417 binEdges[3][6]= 6.;
1418 binEdges[3][7]= 8.;
1419 binEdges[3][8]= 10.;
1420 binEdges[3][9]= 12.;
1421 binEdges[3][10]= kPtLeadingTrackPtMax;
1422
1423 fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
b99018c6 1424 for (int k=0; k<kNvar; k++) {
1425 fhnJetContainer->SetBinLimits(k,binEdges[k]);
3b7ffecf 1426 }
b99018c6 1427
3b7ffecf 1428 //create correlation matrix for unfolding
1429 Int_t thnDim[2*kNvar];
1430 for (int k=0; k<kNvar; k++) {
1431 //first half : reconstructed
1432 //second half : MC
1433 thnDim[k] = iBin[k];
1434 thnDim[k+kNvar] = iBin[k];
1435 }
1436
152a0744 1437 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
3b7ffecf 1438 for (int k=0; k<kNvar; k++) {
152a0744 1439 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1440 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
3b7ffecf 1441 }
742ee86c 1442
1443 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1444 delete [] binEdges[ivar];
1445
1446
3b7ffecf 1447}
1448
1449void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1450{
adf3920d 1451 // Terminate analysis
1452 //
1453 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
3b7ffecf 1454}
1455
1456
1457Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1458
565584e8 1459 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1460
1461 Int_t iCount = 0;
565584e8 1462 if(type==kTrackAOD){
3b7ffecf 1463 AliAODEvent *aod = 0;
565584e8 1464 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1465 else aod = AODEvent();
3b7ffecf 1466 if(!aod){
1467 return iCount;
1468 }
1469 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1470 AliAODTrack *tr = aod->GetTrack(it);
1471 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 1472 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1473 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 1474 if(fDebug>0){
1475 if(tr->Pt()>20){
1476 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 1477 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1478 tr->Print();
1479 // tr->Dump();
1480 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1481 if(fESD){
a2522483 1482 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
38ecb6a5 1483 esdTr->Print("");
bb3d13aa 1484 // esdTr->Dump();
38ecb6a5 1485 }
1486 }
1487 }
3b7ffecf 1488 list->Add(tr);
1489 iCount++;
1490 }
1491 }
1492 else if (type == kTrackKineAll||type == kTrackKineCharged){
1493 AliMCEvent* mcEvent = MCEvent();
1494 if(!mcEvent)return iCount;
1495 // we want to have alivpartilces so use get track
1496 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1497 if(!mcEvent->IsPhysicalPrimary(it))continue;
1498 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 1499 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1500 if(type == kTrackKineAll){
1501 list->Add(part);
1502 iCount++;
1503 }
1504 else if(type == kTrackKineCharged){
1505 if(part->Particle()->GetPDG()->Charge()==0)continue;
1506 list->Add(part);
1507 iCount++;
1508 }
1509 }
1510 }
565584e8 1511 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1512 AliAODEvent *aod = 0;
5301f754 1513 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1514 else aod = AODEvent();
5010a3f7 1515 if(!aod)return iCount;
3b7ffecf 1516 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1517 if(!tca)return iCount;
1518 for(int it = 0;it < tca->GetEntriesFast();++it){
1519 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
e855f5c5 1520 if(!part)continue;
a31b8a87 1521 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1522 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1523 if(type == kTrackAODMCAll){
3b7ffecf 1524 list->Add(part);
1525 iCount++;
1526 }
565584e8 1527 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1528 if(part->Charge()==0)continue;
565584e8 1529 if(kTrackAODMCCharged){
1530 list->Add(part);
1531 }
1532 else {
a31b8a87 1533 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 1534 list->Add(part);
1535 }
3b7ffecf 1536 iCount++;
1537 }
1538 }
1539 }// AODMCparticle
cc0649e4 1540 list->Sort();
c4631cdb 1541 return iCount;
3b7ffecf 1542
1543}
a31b8a87 1544
1545
742ee86c 1546Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
adf3920d 1547 AliAODEvent *aod = 0;
1548 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1549 else aod = AODEvent();
1550 if(!aod){
1551 return 101;
1552 }
1553 return aod->GetHeader()->GetCentrality();
742ee86c 1554}
1555
1556
a31b8a87 1557
1558Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1559 Bool_t selected = false;
1560
1561 if(!jet)return selected;
1562
1563 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1564 selected = kTRUE;
1565 }
1566 return selected;
1567
1568}
1569
1570Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1571
1572 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1573 Int_t iCount = 0;
1574
1575 if(!jarray){
1576 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1577 return iCount;
1578 }
1579
1580
1581 for(int ij=0;ij<jarray->GetEntries();ij++){
1582 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1583 if(!jet)continue;
1584 if(type==0){
1585 // No cut at all, main purpose here is sorting
1586 list->Add(jet);
1587 iCount++;
1588 }
1589 else if(type == 1){
1590 // eta cut
1591 if(JetSelected(jet)){
1592 list->Add(jet);
1593 iCount++;
1594 }
1595 }
1596 }
1597
1598 list->Sort();
1599 return iCount;
1600
1601}
1602
1603
37eb26ea 1604Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1605 if(!jets)return 0;
1606
1607 Int_t refMult = 0;
1608 for(int ij = 0;ij < jets->GetEntries();++ij){
1609 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1610 if(!jet)continue;
1611 TRefArray *refs = jet->GetRefTracks();
1612 if(!refs)continue;
1613 refMult += refs->GetEntries();
1614 }
1615 return refMult;
1616
1617}
d7117cbd 1618
1619
1620AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1621 if(!jet)return 0;
1622 TRefArray *refs = jet->GetRefTracks();
1623 if(!refs) return 0;
1624 AliVParticle *leading = 0;
1625 Float_t fMaxPt = 0;
1626 for(int i = 0;i<refs->GetEntries();i++){
1627 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1628 if(!tmp)continue;
1629 if(tmp->Pt()>fMaxPt){
1630 leading = tmp;
1631 fMaxPt = tmp->Pt();
1632 }
1633 }
1634 return leading;
1635}
1636
1637
1638AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1639 if(!jet)return 0;
1640 if(!list) return 0;
1641 AliVParticle *leading = 0;
1642 Float_t fMaxPt = 0;
1643 for(int i = 0;i<list->GetEntries();i++){
1644 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1645 if(!tmp)continue;
1646 if(jet->DeltaR(tmp)>r)continue;
1647 if(tmp->Pt()>fMaxPt){
1648 leading = tmp;
1649 fMaxPt = tmp->Pt();
1650 }
1651 }
1652 return leading;
1653}
742ee86c 1654
742ee86c 1655Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1656{
adf3920d 1657 Int_t phibin=-1;
1658 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1659 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1660 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1661 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1662 return phibin;
742ee86c 1663}
1664
3d1dba29 1665void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
db745d29 1666 //
1667 // set number of trigger bins
1668 //
3d1dba29 1669 if(n>0){
adf3920d 1670 fNTrigger = n;
1671 delete [] fTriggerName;
1672 fTriggerName = new TString [fNTrigger];
1673 delete [] fTriggerBit;fTriggerBit = 0;
1674 fTriggerBit = new UInt_t [fNTrigger];
3d1dba29 1675 }
1676 else{
1677 fNTrigger = 0;
1678 }
1679}
1680
db745d29 1681
3d1dba29 1682void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
db745d29 1683 //
1684 // set trigger bin
1685 //
3d1dba29 1686 if(i<fNTrigger){
1687 fTriggerBit[i] = it;
1688 fTriggerName[i] = c;
1689 }
1690}
1691
db745d29 1692
1693void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1694 //
1695 // Set number of acceptance bins
1696 //
1697
1698 if(n>0){
adf3920d 1699 fNAcceptance = n;
1700 delete [] fAcceptancePhiMin;
1701 delete [] fAcceptancePhiMax;
1702 delete [] fAcceptanceEtaMin;
1703 delete [] fAcceptanceEtaMax;
1704
1705 fAcceptancePhiMin = new Float_t[fNAcceptance];
1706 fAcceptancePhiMax = new Float_t[fNAcceptance];
1707 fAcceptanceEtaMin = new Float_t[fNAcceptance];
1708 fAcceptanceEtaMax = new Float_t[fNAcceptance];
db745d29 1709 }
1710 else{
1711 fNTrigger = 0;
1712 }
1713}
1714
1715void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1716 //
1717 // Set acceptance bins
1718 //
1719
1720 if(i<fNAcceptance){
24cc3482 1721 Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1722
1723 fAcceptancePhiMin[i] = phiMin;
1724 fAcceptancePhiMax[i] = phiMax;
1725 fAcceptanceEtaMin[i] = etaMin;
1726 fAcceptanceEtaMax[i] = etaMax;
db745d29 1727 }
1728 else{
1729 fNTrigger = 0;
1730 }
1731}
1732
1733
1734Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1735
1736 //
1737 // return aceptnace bin do no use ovelapping bins
1738 //
1739
24cc3482 1740 for(int i = 0;i<fNAcceptance;i++){
db745d29 1741 if(phi<fAcceptancePhiMin[i])continue;
1742 if(phi>fAcceptancePhiMax[i])continue;
1743 if(eta<fAcceptanceEtaMin[i])continue;
1744 if(eta>fAcceptanceEtaMax[i])continue;
1745 return i;
1746 }
1747 // catch the rest
1748 return fNAcceptance;
1749}
1750
e5ee1c19 1751Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
db745d29 1752
e5ee1c19 1753 // invert the correction
1754 AliAODJet *jet = (AliAODJet*)list.At(0); // highest pt jet
1755 if(!jet)return -1;
1756 if(jet->EffectiveAreaCharged()<=0)return -1;
1757 Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1758 return rho;
1759
1760}
db745d29 1761
1762
1763
3d1dba29 1764AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1765 //
1766 delete [] fTriggerBit;
1767 delete [] fTriggerName;
db745d29 1768
1769 delete [] fAcceptancePhiMin;
1770 delete [] fAcceptancePhiMax;
1771 delete [] fAcceptanceEtaMin;
1772 delete [] fAcceptanceEtaMax;
1773
1774
1775
3d1dba29 1776}