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