Possibility to select the Lc resonant channel in the MC analysis (Rossella)
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
CommitLineData
dbf053f0 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
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>
d7ca0e14 24#include <TRandom3.h>
dbf053f0 25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
54424110 28#include <TRefArray.h>
dbf053f0 29#include <TFile.h>
30#include <TKey.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TH3F.h>
34#include <TProfile.h>
006b2a30 35#include <TF1.h>
dbf053f0 36#include <TList.h>
37#include <TLorentzVector.h>
38#include <TClonesArray.h>
39#include "TDatabasePDG.h"
40
41#include "AliAnalysisTaskJetCluster.h"
42#include "AliAnalysisManager.h"
43#include "AliJetFinder.h"
44#include "AliJetHeader.h"
45#include "AliJetReader.h"
dbf053f0 46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
1240edf5 49#include "AliAODExtension.h"
dbf053f0 50#include "AliAODTrack.h"
51#include "AliAODJet.h"
52#include "AliAODMCParticle.h"
53#include "AliMCEventHandler.h"
54#include "AliMCEvent.h"
55#include "AliStack.h"
56#include "AliGenPythiaEventHeader.h"
57#include "AliJetKineReaderHeader.h"
58#include "AliGenCocktailEventHeader.h"
59#include "AliInputEventHandler.h"
d6e66a82 60#include "AliAODJetEventBackground.h"
dbf053f0 61
62#include "fastjet/PseudoJet.hh"
63#include "fastjet/ClusterSequenceArea.hh"
64#include "fastjet/AreaDefinition.hh"
65#include "fastjet/JetDefinition.hh"
66// get info on how fastjet was configured
67#include "fastjet/config.h"
68
69
70ClassImp(AliAnalysisTaskJetCluster)
71
54424110 72AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
f0659f11 73 //
74 // Destructor
75 //
006b2a30 76
54424110 77 delete fRef;
d7ca0e14 78 delete fRandom;
8f85b773 79
80 if(fTCAJetsOut)fTCAJetsOut->Delete();
81 delete fTCAJetsOut;
006b2a30 82
8f85b773 83 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
84 delete fTCAJetsOutRan;
006b2a30 85
8f85b773 86 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
87 delete fTCARandomConesOut;
006b2a30 88
8f85b773 89 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
90 delete fTCARandomConesOutRan;
006b2a30 91
8f85b773 92 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
93 delete fAODJetBackgroundOut;
54424110 94}
95
b2c1e955 96AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
97 AliAnalysisTaskSE(),
dbf053f0 98 fAOD(0x0),
54424110 99 fAODExtension(0x0),
100 fRef(new TRefArray),
dbf053f0 101 fUseAODTrackInput(kFALSE),
102 fUseAODMCInput(kFALSE),
9fd9f55b 103 fUseBackgroundCalc(kFALSE),
b2c1e955 104 fEventSelection(kFALSE),
dbf053f0 105 fFilterMask(0),
d7396b04 106 fFilterMaskBestPt(0),
b2c1e955 107 fFilterType(0),
bea24d2d 108 fJetTypes(kJet),
dbf053f0 109 fTrackTypeRec(kTrackUndef),
bd80a748 110 fTrackTypeGen(kTrackUndef),
111 fNSkipLeadingRan(0),
af399bc2 112 fNSkipLeadingCone(0),
c4856527 113 fNRandomCones(0),
bd80a748 114 fAvgTrials(1),
8a320ab9 115 fExternalWeight(1),
116 fTrackEtaWindow(0.9),
dbf053f0 117 fRecEtaWindow(0.5),
54424110 118 fTrackPtCut(0.),
617932d0 119 fJetOutputMinPt(0.150),
65c43de8 120 fMaxTrackPtInJet(100.),
0341d0d8 121 fJetTriggerPtCut(0),
ca6d12f9 122 fVtxZCut(8),
123 fVtxR2Cut(1),
aae7993b 124 fCentCutUp(0),
125 fCentCutLo(0),
54424110 126 fNonStdBranch(""),
0341d0d8 127 fBackgroundBranch(""),
54424110 128 fNonStdFile(""),
006b2a30 129 fMomResH1(0x0),
130 fMomResH2(0x0),
131 fMomResH3(0x0),
132 fMomResH1Fit(0x0),
133 fMomResH2Fit(0x0),
134 fMomResH3Fit(0x0),
135 fhEffH1(0x0),
136 fhEffH2(0x0),
137 fhEffH3(0x0),
138 fUseTrMomentumSmearing(kFALSE),
139 fUseDiceEfficiency(kFALSE),
dbf053f0 140 fRparam(1.0),
141 fAlgorithm(fastjet::kt_algorithm),
142 fStrategy(fastjet::Best),
143 fRecombScheme(fastjet::BIpt_scheme),
144 fAreaType(fastjet::active_area),
d6e66a82 145 fGhostArea(0.01),
146 fActiveAreaRepeats(1),
147 fGhostEtamax(1.5),
8f85b773 148 fTCAJetsOut(0x0),
149 fTCAJetsOutRan(0x0),
150 fTCARandomConesOut(0x0),
151 fTCARandomConesOutRan(0x0),
152 fAODJetBackgroundOut(0x0),
d7ca0e14 153 fRandom(0),
dbf053f0 154 fh1Xsec(0x0),
155 fh1Trials(0x0),
156 fh1PtHard(0x0),
157 fh1PtHardNoW(0x0),
158 fh1PtHardTrials(0x0),
159 fh1NJetsRec(0x0),
160 fh1NConstRec(0x0),
161 fh1NConstLeadingRec(0x0),
162 fh1PtJetsRecIn(0x0),
163 fh1PtJetsLeadingRecIn(0x0),
164 fh1PtJetConstRec(0x0),
165 fh1PtJetConstLeadingRec(0x0),
166 fh1PtTracksRecIn(0x0),
167 fh1PtTracksLeadingRecIn(0x0),
168 fh1NJetsRecRan(0x0),
169 fh1NConstRecRan(0x0),
170 fh1PtJetsLeadingRecInRan(0x0),
171 fh1NConstLeadingRecRan(0x0),
172 fh1PtJetsRecInRan(0x0),
173 fh1PtTracksGenIn(0x0),
a0d08b6d 174 fh1Nch(0x0),
aae7993b 175 fh1CentralityPhySel(0x0),
0341d0d8 176 fh1Centrality(0x0),
177 fh1CentralitySelect(0x0),
e2ca7519 178 fh1ZPhySel(0x0),
179 fh1Z(0x0),
180 fh1ZSelect(0x0),
dbf053f0 181 fh2NRecJetsPt(0x0),
182 fh2NRecTracksPt(0x0),
183 fh2NConstPt(0x0),
184 fh2NConstLeadingPt(0x0),
185 fh2JetPhiEta(0x0),
186 fh2LeadingJetPhiEta(0x0),
187 fh2JetEtaPt(0x0),
188 fh2LeadingJetEtaPt(0x0),
189 fh2TrackEtaPt(0x0),
190 fh2LeadingTrackEtaPt(0x0),
191 fh2JetsLeadingPhiEta(0x0),
192 fh2JetsLeadingPhiPt(0x0),
193 fh2TracksLeadingPhiEta(0x0),
194 fh2TracksLeadingPhiPt(0x0),
195 fh2TracksLeadingJetPhiPt(0x0),
196 fh2JetsLeadingPhiPtW(0x0),
197 fh2TracksLeadingPhiPtW(0x0),
198 fh2TracksLeadingJetPhiPtW(0x0),
199 fh2NRecJetsPtRan(0x0),
200 fh2NConstPtRan(0x0),
201 fh2NConstLeadingPtRan(0x0),
a0d08b6d 202 fh2PtNch(0x0),
203 fh2PtNchRan(0x0),
204 fh2PtNchN(0x0),
205 fh2PtNchNRan(0x0),
dbf053f0 206 fh2TracksLeadingJetPhiPtRan(0x0),
207 fh2TracksLeadingJetPhiPtWRan(0x0),
006b2a30 208 fh2PtGenPtSmeared(0x0),
209 fp1Efficiency(0x0),
210 fp1PtResolution(0x0),
dbf053f0 211 fHistList(0x0)
212{
f0659f11 213 //
214 // Constructor
215 //
216
82a35b1a 217 for(int i = 0;i<3;i++){
218 fh1BiARandomCones[i] = 0;
219 fh1BiARandomConesRan[i] = 0;
220 }
aae7993b 221 for(int i = 0;i<kMaxCent;i++){
222 fh2JetsLeadingPhiPtC[i] = 0;
223 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
224 fh2TracksLeadingJetPhiPtC[i] = 0;
225 fh2TracksLeadingJetPhiPtWC[i] = 0;
226 }
dbf053f0 227}
228
229AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
230 AliAnalysisTaskSE(name),
231 fAOD(0x0),
54424110 232 fAODExtension(0x0),
233 fRef(new TRefArray),
dbf053f0 234 fUseAODTrackInput(kFALSE),
235 fUseAODMCInput(kFALSE),
9fd9f55b 236 fUseBackgroundCalc(kFALSE),
d7396b04 237 fEventSelection(kFALSE), fFilterMask(0),
238 fFilterMaskBestPt(0),
b2c1e955 239 fFilterType(0),
bea24d2d 240 fJetTypes(kJet),
dbf053f0 241 fTrackTypeRec(kTrackUndef),
242 fTrackTypeGen(kTrackUndef),
bd80a748 243 fNSkipLeadingRan(0),
af399bc2 244 fNSkipLeadingCone(0),
c4856527 245 fNRandomCones(0),
dbf053f0 246 fAvgTrials(1),
247 fExternalWeight(1),
8a320ab9 248 fTrackEtaWindow(0.9),
dbf053f0 249 fRecEtaWindow(0.5),
54424110 250 fTrackPtCut(0.),
617932d0 251 fJetOutputMinPt(0.150),
65c43de8 252 fMaxTrackPtInJet(100.),
0341d0d8 253 fJetTriggerPtCut(0),
ca6d12f9 254 fVtxZCut(8),
255 fVtxR2Cut(1),
aae7993b 256 fCentCutUp(0),
257 fCentCutLo(0),
54424110 258 fNonStdBranch(""),
0341d0d8 259 fBackgroundBranch(""),
54424110 260 fNonStdFile(""),
006b2a30 261 fMomResH1(0x0),
262 fMomResH2(0x0),
263 fMomResH3(0x0),
264 fMomResH1Fit(0x0),
265 fMomResH2Fit(0x0),
266 fMomResH3Fit(0x0),
267 fhEffH1(0x0),
268 fhEffH2(0x0),
269 fhEffH3(0x0),
270 fUseTrMomentumSmearing(kFALSE),
271 fUseDiceEfficiency(kFALSE),
dbf053f0 272 fRparam(1.0),
273 fAlgorithm(fastjet::kt_algorithm),
274 fStrategy(fastjet::Best),
275 fRecombScheme(fastjet::BIpt_scheme),
276 fAreaType(fastjet::active_area),
d6e66a82 277 fGhostArea(0.01),
278 fActiveAreaRepeats(1),
279 fGhostEtamax(1.5),
8f85b773 280 fTCAJetsOut(0x0),
281 fTCAJetsOutRan(0x0),
282 fTCARandomConesOut(0x0),
283 fTCARandomConesOutRan(0x0),
284 fAODJetBackgroundOut(0x0),
d7ca0e14 285 fRandom(0),
dbf053f0 286 fh1Xsec(0x0),
287 fh1Trials(0x0),
288 fh1PtHard(0x0),
289 fh1PtHardNoW(0x0),
290 fh1PtHardTrials(0x0),
291 fh1NJetsRec(0x0),
292 fh1NConstRec(0x0),
293 fh1NConstLeadingRec(0x0),
294 fh1PtJetsRecIn(0x0),
295 fh1PtJetsLeadingRecIn(0x0),
296 fh1PtJetConstRec(0x0),
297 fh1PtJetConstLeadingRec(0x0),
298 fh1PtTracksRecIn(0x0),
299 fh1PtTracksLeadingRecIn(0x0),
300 fh1NJetsRecRan(0x0),
301 fh1NConstRecRan(0x0),
302 fh1PtJetsLeadingRecInRan(0x0),
303 fh1NConstLeadingRecRan(0x0),
304 fh1PtJetsRecInRan(0x0),
305 fh1PtTracksGenIn(0x0),
a0d08b6d 306 fh1Nch(0x0),
aae7993b 307 fh1CentralityPhySel(0x0),
0341d0d8 308 fh1Centrality(0x0),
309 fh1CentralitySelect(0x0),
e2ca7519 310 fh1ZPhySel(0x0),
311 fh1Z(0x0),
312 fh1ZSelect(0x0),
dbf053f0 313 fh2NRecJetsPt(0x0),
314 fh2NRecTracksPt(0x0),
315 fh2NConstPt(0x0),
316 fh2NConstLeadingPt(0x0),
317 fh2JetPhiEta(0x0),
318 fh2LeadingJetPhiEta(0x0),
319 fh2JetEtaPt(0x0),
320 fh2LeadingJetEtaPt(0x0),
321 fh2TrackEtaPt(0x0),
322 fh2LeadingTrackEtaPt(0x0),
323 fh2JetsLeadingPhiEta(0x0),
324 fh2JetsLeadingPhiPt(0x0),
325 fh2TracksLeadingPhiEta(0x0),
326 fh2TracksLeadingPhiPt(0x0),
327 fh2TracksLeadingJetPhiPt(0x0),
328 fh2JetsLeadingPhiPtW(0x0),
329 fh2TracksLeadingPhiPtW(0x0),
330 fh2TracksLeadingJetPhiPtW(0x0),
331 fh2NRecJetsPtRan(0x0),
332 fh2NConstPtRan(0x0),
333 fh2NConstLeadingPtRan(0x0),
a0d08b6d 334 fh2PtNch(0x0),
335 fh2PtNchRan(0x0),
336 fh2PtNchN(0x0),
337 fh2PtNchNRan(0x0),
dbf053f0 338 fh2TracksLeadingJetPhiPtRan(0x0),
339 fh2TracksLeadingJetPhiPtWRan(0x0),
006b2a30 340 fh2PtGenPtSmeared(0x0),
341 fp1Efficiency(0x0),
342 fp1PtResolution(0x0),
dbf053f0 343 fHistList(0x0)
344{
f0659f11 345 //
346 // named ctor
347 //
006b2a30 348
82a35b1a 349 for(int i = 0;i<3;i++){
350 fh1BiARandomCones[i] = 0;
351 fh1BiARandomConesRan[i] = 0;
352 }
aae7993b 353 for(int i = 0;i<kMaxCent;i++){
354 fh2JetsLeadingPhiPtC[i] = 0;
355 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
356 fh2TracksLeadingJetPhiPtC[i] = 0;
357 fh2TracksLeadingJetPhiPtWC[i] = 0;
358 }
82a35b1a 359 DefineOutput(1, TList::Class());
dbf053f0 360}
361
362
363
364Bool_t AliAnalysisTaskJetCluster::Notify()
365{
366 //
367 // Implemented Notify() to read the cross sections
368 // and number of trials from pyxsec.root
369 //
370 return kTRUE;
371}
372
373void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
374{
375
376 //
377 // Create the output container
378 //
379
d7ca0e14 380 fRandom = new TRandom3(0);
54424110 381
382
dbf053f0 383 // Connect the AOD
384
385
386 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
387
d6e66a82 388
54424110 389
390 if(fNonStdBranch.Length()!=0)
391 {
392 // only create the output branch if we have a name
393 // Create a new branch for jets...
394 // -> cleared in the UserExec....
395 // here we can also have the case that the brnaches are written to a separate file
d508795b 396
bea24d2d 397 if(fJetTypes&kJet){
398 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
399 fTCAJetsOut->SetName(fNonStdBranch.Data());
400 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
401 }
006b2a30 402
bea24d2d 403 if(fJetTypes&kJetRan){
404 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
405 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
006b2a30 406 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
407 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%d",fNonStdBranch.Data(),"random",fUseTrMomentumSmearing,fUseDiceEfficiency));
bea24d2d 408 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
409 }
410
9fd9f55b 411 if(fUseBackgroundCalc){
412 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
8f85b773 413 fAODJetBackgroundOut = new AliAODJetEventBackground();
414 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
006b2a30 415 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
416 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency));
417
8f85b773 418 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
9fd9f55b 419 }
c4856527 420 }
421 // create the branch for the random cones with the same R
af399bc2 422 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
006b2a30 423 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
424 cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
425
c4856527 426 if(fNRandomCones>0){
bea24d2d 427 if(fJetTypes&kRC){
428 if(!AODEvent()->FindListObject(cName.Data())){
429 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
430 fTCARandomConesOut->SetName(cName.Data());
431 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
432 }
d7ca0e14 433 }
d7ca0e14 434 // create the branch with the random for the random cones on the random event
bea24d2d 435 if(fJetTypes&kRCRan){
436 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
437 if(!AODEvent()->FindListObject(cName.Data())){
438 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
439 fTCARandomConesOutRan->SetName(cName.Data());
440 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
441 }
d7ca0e14 442 }
9fd9f55b 443 }
c4856527 444
54424110 445 if(fNonStdFile.Length()!=0){
446 //
447 // case that we have an AOD extension we need to fetch the jets from the extended output
f840d64c 448 // we identify the extension aod event by looking for the branchname
54424110 449 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
8f85b773 450 // case that we have an AOD extension we need can fetch the background maybe from the extended output
451 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
54424110 452 }
453 }
454
006b2a30 455 // FitMomentumResolution();
456
54424110 457
dbf053f0 458 if(!fHistList)fHistList = new TList();
9c6bd749 459 fHistList->SetOwner();
750e22e8 460 PostData(1, fHistList); // post data in any case once
461
dbf053f0 462 Bool_t oldStatus = TH1::AddDirectoryStatus();
463 TH1::AddDirectory(kFALSE);
464
465 //
466 // Histogram
467
fe80c230 468 const Int_t nBinPt = 100;
dbf053f0 469 Double_t binLimitsPt[nBinPt+1];
470 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
471 if(iPt == 0){
472 binLimitsPt[iPt] = 0.0;
473 }
474 else {// 1.0
fe80c230 475 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
dbf053f0 476 }
477 }
478
479 const Int_t nBinPhi = 90;
480 Double_t binLimitsPhi[nBinPhi+1];
481 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
482 if(iPhi==0){
483 binLimitsPhi[iPhi] = -1.*TMath::Pi();
484 }
485 else{
486 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
487 }
488 }
489
490
491
492 const Int_t nBinEta = 40;
493 Double_t binLimitsEta[nBinEta+1];
494 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
495 if(iEta==0){
496 binLimitsEta[iEta] = -2.0;
497 }
498 else{
499 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
500 }
501 }
502
fe80c230 503 const int nChMax = 5000;
dbf053f0 504
505 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
506 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
507
508 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
509 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
510
511
512 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
513 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
514
515 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
516 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
517 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
518 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
519
520
521 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
522 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
523 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
524
525 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
526 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
527 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
528 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
529 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
530 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
8a320ab9 531 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
532 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
533 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
a0d08b6d 534 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
dbf053f0 535
0341d0d8 536 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
537 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
aae7993b 538 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
0341d0d8 539
e2ca7519 540 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
541 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
542 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
543
dbf053f0 544 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
545 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
546 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
547 //
548
549
550 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
551 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
552 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
553 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
554
a0d08b6d 555 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
556 fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
557 fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
558 fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
559
560
561
dbf053f0 562 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
563 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
564 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
565 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
566
567 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
568 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
569 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
570 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
571
572 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
573 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
574 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
575 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
576
577
578
579 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
580 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
581 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
582 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
583 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
584 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
585 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
586 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
587 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
588 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
589 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
590 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
591
592 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
593 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
594 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
595 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
596
597 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
598 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
599 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
600 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
601
006b2a30 602 //Detector level effects histos
603 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
604
605 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
606 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
607
608 fHistList->Add(fh2PtGenPtSmeared);
609 fHistList->Add(fp1Efficiency);
610 fHistList->Add(fp1PtResolution);
dbf053f0 611
c4856527 612 if(fNRandomCones>0&&fUseBackgroundCalc){
d7ca0e14 613 for(int i = 0;i<3;i++){
614 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
615 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
616 }
617 }
c4856527 618
aae7993b 619 for(int i = 0;i < kMaxCent;i++){
620 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
621 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
622 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
623 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
624 }
dbf053f0 625
626 const Int_t saveLevel = 3; // large save level more histos
627 if(saveLevel>0){
628 fHistList->Add(fh1Xsec);
629 fHistList->Add(fh1Trials);
630
631 fHistList->Add(fh1NJetsRec);
632 fHistList->Add(fh1NConstRec);
633 fHistList->Add(fh1NConstLeadingRec);
634 fHistList->Add(fh1PtJetsRecIn);
635 fHistList->Add(fh1PtJetsLeadingRecIn);
636 fHistList->Add(fh1PtTracksRecIn);
637 fHistList->Add(fh1PtTracksLeadingRecIn);
638 fHistList->Add(fh1PtJetConstRec);
639 fHistList->Add(fh1PtJetConstLeadingRec);
640 fHistList->Add(fh1NJetsRecRan);
641 fHistList->Add(fh1NConstRecRan);
642 fHistList->Add(fh1PtJetsLeadingRecInRan);
643 fHistList->Add(fh1NConstLeadingRecRan);
644 fHistList->Add(fh1PtJetsRecInRan);
a0d08b6d 645 fHistList->Add(fh1Nch);
0341d0d8 646 fHistList->Add(fh1Centrality);
647 fHistList->Add(fh1CentralitySelect);
aae7993b 648 fHistList->Add(fh1CentralityPhySel);
e2ca7519 649 fHistList->Add(fh1Z);
650 fHistList->Add(fh1ZSelect);
651 fHistList->Add(fh1ZPhySel);
8f85b773 652 if(fNRandomCones>0&&fUseBackgroundCalc){
d7ca0e14 653 for(int i = 0;i<3;i++){
654 fHistList->Add(fh1BiARandomCones[i]);
655 fHistList->Add(fh1BiARandomConesRan[i]);
656 }
82a35b1a 657 }
aae7993b 658 for(int i = 0;i < kMaxCent;i++){
659 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
660 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
661 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
662 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
663 }
664
dbf053f0 665 fHistList->Add(fh2NRecJetsPt);
666 fHistList->Add(fh2NRecTracksPt);
667 fHistList->Add(fh2NConstPt);
668 fHistList->Add(fh2NConstLeadingPt);
a0d08b6d 669 fHistList->Add(fh2PtNch);
670 fHistList->Add(fh2PtNchRan);
671 fHistList->Add(fh2PtNchN);
672 fHistList->Add(fh2PtNchNRan);
dbf053f0 673 fHistList->Add(fh2JetPhiEta);
674 fHistList->Add(fh2LeadingJetPhiEta);
675 fHistList->Add(fh2JetEtaPt);
676 fHistList->Add(fh2LeadingJetEtaPt);
677 fHistList->Add(fh2TrackEtaPt);
678 fHistList->Add(fh2LeadingTrackEtaPt);
679 fHistList->Add(fh2JetsLeadingPhiEta );
680 fHistList->Add(fh2JetsLeadingPhiPt);
681 fHistList->Add(fh2TracksLeadingPhiEta);
682 fHistList->Add(fh2TracksLeadingPhiPt);
683 fHistList->Add(fh2TracksLeadingJetPhiPt);
684 fHistList->Add(fh2JetsLeadingPhiPtW);
685 fHistList->Add(fh2TracksLeadingPhiPtW);
686 fHistList->Add(fh2TracksLeadingJetPhiPtW);
687 fHistList->Add(fh2NRecJetsPtRan);
688 fHistList->Add(fh2NConstPtRan);
689 fHistList->Add(fh2NConstLeadingPtRan);
690 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
691 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
85587ff2 692 }
dbf053f0 693
694 // =========== Switch on Sumw2 for all histos ===========
695 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
696 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
697 if (h1){
698 h1->Sumw2();
699 continue;
700 }
701 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
702 if(hn)hn->Sumw2();
703 }
704 TH1::AddDirectory(oldStatus);
705}
706
707void AliAnalysisTaskJetCluster::Init()
708{
709 //
710 // Initialization
711 //
712
dbf053f0 713 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
714
006b2a30 715 FitMomentumResolution();
716
dbf053f0 717}
718
719void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
720{
721
54424110 722 // handle and reset the output jet branch
c4856527 723
8f85b773 724 if(fTCAJetsOut)fTCAJetsOut->Delete();
725 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
726 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
727 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
728 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
54424110 729
a18eedbb 730 AliAODJetEventBackground* externalBackground = 0;
0341d0d8 731 if(!externalBackground&&fBackgroundBranch.Length()){
732 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
750e22e8 733 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
8f85b773 734 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
0341d0d8 735 }
dbf053f0 736 //
737 // Execute analysis for current event
738 //
739 AliESDEvent *fESD = 0;
740 if(fUseAODTrackInput){
741 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
742 if(!fAOD){
743 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
744 return;
745 }
006b2a30 746 // fetch the header
dbf053f0 747 }
748 else{
749 // assume that the AOD is in the general output...
750 fAOD = AODEvent();
751 if(!fAOD){
752 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
753 return;
754 }
755 if(fDebug>0){
756 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
757 }
758 }
006b2a30 759
760 //Check if information is provided detector level effects
761 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
762 if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
dbf053f0 763
764 Bool_t selectEvent = false;
54424110 765 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
dbf053f0 766
aae7993b 767 Float_t cent = 0;
e2ca7519 768 Float_t zVtx = 0;
aae7993b 769 Int_t cenClass = -1;
dbf053f0 770 if(fAOD){
771 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
54424110 772 TString vtxTitle(vtxAOD->GetTitle());
e2ca7519 773 zVtx = vtxAOD->GetZ();
774
aae7993b 775 cent = fAOD->GetHeader()->GetCentrality();
776 if(cent<10)cenClass = 0;
777 else if(cent<30)cenClass = 1;
778 else if(cent<50)cenClass = 2;
779 else if(cent<80)cenClass = 3;
e2ca7519 780 if(physicsSelection){
781 fh1CentralityPhySel->Fill(cent);
782 fh1ZPhySel->Fill(zVtx);
783 }
784
ca6d12f9 785 if(fEventSelection){
786 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
dbf053f0 787 Float_t yvtx = vtxAOD->GetY();
788 Float_t xvtx = vtxAOD->GetX();
789 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
ca6d12f9 790 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
e2ca7519 791 if(physicsSelection){
792 selectEvent = true;
793 }
dbf053f0 794 }
aae7993b 795 }
ca6d12f9 796 if(fCentCutUp>0){
797 if(cent<fCentCutLo||cent>fCentCutUp){
798 selectEvent = false;
799 }
800 }
801 }else{
802 selectEvent = true;
aae7993b 803 }
dbf053f0 804 }
ca6d12f9 805
806
dbf053f0 807 if(!selectEvent){
808 PostData(1, fHistList);
809 return;
810 }
aae7993b 811 fh1Centrality->Fill(cent);
e2ca7519 812 fh1Z->Fill(zVtx);
dbf053f0 813 fh1Trials->Fill("#sum{ntrials}",1);
814
815
816 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
817
818 // ==== General variables needed
819
820
821
822 // we simply fetch the tracks/mc particles as a list of AliVParticles
823
824 TList recParticles;
825 TList genParticles;
826
827 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
a0d08b6d 828 Float_t nCh = recParticles.GetEntries();
829 fh1Nch->Fill(nCh);
dbf053f0 830 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
831 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
832 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
833
834 // find the jets....
835
836 vector<fastjet::PseudoJet> inputParticlesRec;
837 vector<fastjet::PseudoJet> inputParticlesRecRan;
d7ca0e14 838
839 // Generate the random cones
840
841 AliAODJet vTmpRan(1,0,0,1);
dbf053f0 842 for(int i = 0; i < recParticles.GetEntries(); i++){
843 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
006b2a30 844
dbf053f0 845 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
d7ca0e14 846 // we take total momentum here
006b2a30 847
848 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
849 //Add particles to fastjet in case we are not running toy model
850 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
851 jInp.set_user_index(i);
852 inputParticlesRec.push_back(jInp);
853 }
854 else if(fUseDiceEfficiency) {
855
856 // Dice to decide if particle is kept or not - toy model for efficiency
857 //
858 Double_t rnd = fRandom->Uniform(1.);
859 Double_t pT = vp->Pt();
860 Double_t eff[3] = {0.};
861 Double_t pTtmp = pT;
862 if(pT>10.) pTtmp = 10.;
863 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
864 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
865 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
866 Int_t cat[3] = {0};
867 //Sort efficiencies from large to small
868 if(eff1>eff2 && eff1>eff3) {
869 eff[0] = eff1;
870 cat[0] = 1;
871 if(eff2>eff3) {
872 eff[1] = eff2;
873 eff[2] = eff3;
874 cat[1] = 2;
875 cat[2] = 3;
876 } else {
877 eff[1] = eff3;
878 eff[2] = eff2;
879 cat[1] = 3;
880 cat[2] = 2;
881 }
882 }
883 else if(eff2>eff1 && eff2>eff3) {
884 eff[0] = eff2;
885 cat[0] = 2;
886 if(eff1>eff3) {
887 eff[1] = eff1;
888 eff[2] = eff3;
889 cat[1] = 1;
890 cat[2] = 3;
891 } else {
892 eff[1] = eff3;
893 eff[2] = eff1;
894 cat[1] = 3;
895 cat[2] = 1;
896 }
897 }
898 else if(eff3>eff1 && eff3>eff2) {
899 eff[0] = eff3;
900 cat[0] = 3;
901 if(eff1>eff2) {
902 eff[1] = eff1;
903 eff[2] = eff2;
904 cat[1] = 1;
905 cat[2] = 2;
906 } else {
907 eff[1] = eff2;
908 eff[2] = eff1;
909 cat[1] = 2;
910 cat[2] = 1;
911 }
912 }
913
914 Double_t sumEff = eff[0]+eff[1]+eff[2];
915 fp1Efficiency->Fill(vp->Pt(),sumEff);
916 if(rnd>sumEff) continue;
917
918 if(fUseTrMomentumSmearing) {
919 //Smear momentum of generated particle
920 Double_t smear = 1.;
921 //Select hybrid track category
922 if(rnd<=eff[2])
923 smear = GetMomentumSmearing(cat[2],pT);
924 else if(rnd<=(eff[2]+eff[1]))
925 smear = GetMomentumSmearing(cat[1],pT);
926 else
927 smear = GetMomentumSmearing(cat[0],pT);
928
929 fp1PtResolution->Fill(vp->Pt(),smear);
930
931 Double_t sigma = vp->Pt()*smear;
932 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
933
934 Double_t phi = vp->Phi();
935 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
936 Double_t pX = pTrec * TMath::Cos(phi);
937 Double_t pY = pTrec * TMath::Sin(phi);
938 Double_t pZ = pTrec/TMath::Tan(theta);
939 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
940
941 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
942
943 fastjet::PseudoJet jInp(pX,pY,pZ,p);
944 jInp.set_user_index(i);
945 inputParticlesRec.push_back(jInp);
946
947 }
948 else {
949 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
950 jInp.set_user_index(i);
951 inputParticlesRec.push_back(jInp);
952
953 }
954
955 }
dbf053f0 956
957 // the randomized input changes eta and phi, but keeps the p_T
bd80a748 958 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
959 Double_t pT = vp->Pt();
8a320ab9 960 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
d7ca0e14 961 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
bd80a748 962
d7ca0e14 963 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
bd80a748 964 Double_t pZ = pT/TMath::Tan(theta);
965
966 Double_t pX = pT * TMath::Cos(phi);
967 Double_t pY = pT * TMath::Sin(phi);
968 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
969 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
d7ca0e14 970
bd80a748 971 jInpRan.set_user_index(i);
972 inputParticlesRecRan.push_back(jInpRan);
d7ca0e14 973 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
bd80a748 974 }
54424110 975
976 // fill the tref array, only needed when we write out jets
8f85b773 977 if(fTCAJetsOut){
54424110 978 if(i == 0){
979 fRef->Delete(); // make sure to delete before placement new...
006b2a30 980 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
981 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
982 }
54424110 983 }
006b2a30 984 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
54424110 985 }
d7ca0e14 986 }// recparticles
d7ca0e14 987
dbf053f0 988 if(inputParticlesRec.size()==0){
989 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
990 PostData(1, fHistList);
991 return;
992 }
993
994 // run fast jet
4dd6159d 995 // employ setters for these...
996
d6e66a82 997
9fd9f55b 998 // now create the object that holds info about ghosts
bea24d2d 999 /*
9fd9f55b 1000 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1001 // reduce CPU time...
1002 fGhostArea = 0.5;
1003 fActiveAreaRepeats = 0;
1004 }
bea24d2d 1005 */
1006
1007 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
4dd6159d 1008 fastjet::AreaType areaType = fastjet::active_area;
1009 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
82a35b1a 1010 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
4dd6159d 1011 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
dbf053f0 1012
d6e66a82 1013 //range where to compute background
1014 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1015 phiMin = 0;
1016 phiMax = 2*TMath::Pi();
1017 rapMax = fGhostEtamax - fRparam;
1018 rapMin = - fGhostEtamax + fRparam;
1019 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1020
1021
fe80c230 1022 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1023 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
d6e66a82 1024
dbf053f0 1025
dbf053f0 1026 fh1NJetsRec->Fill(sortedJets.size());
1027
1028 // loop over all jets an fill information, first on is the leading jet
1029
54424110 1030 Int_t nRecOver = inclusiveJets.size();
1031 Int_t nRec = inclusiveJets.size();
1032 if(inclusiveJets.size()>0){
1033 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
aae7993b 1034 Double_t area = clustSeq.area(sortedJets[0]);
0341d0d8 1035 leadingJet.SetEffArea(area,0);
dbf053f0 1036 Float_t pt = leadingJet.Pt();
54424110 1037 Int_t nAodOutJets = 0;
1038 Int_t nAodOutTracks = 0;
1039 AliAODJet *aodOutJet = 0;
dbf053f0 1040
1041 Int_t iCount = 0;
1042 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1043 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1044 while(pt<ptCut&&iCount<nRec){
1045 nRecOver--;
1046 iCount++;
1047 if(iCount<nRec){
1048 pt = sortedJets[iCount].perp();
1049 }
1050 }
1051 if(nRecOver<=0)break;
1052 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1053 }
1054 Float_t phi = leadingJet.Phi();
1055 if(phi<0)phi+=TMath::Pi()*2.;
1056 Float_t eta = leadingJet.Eta();
0341d0d8 1057 Float_t pTback = 0;
1058 if(externalBackground){
1059 // carefull has to be filled in a task before
006b2a30 1060 // todo, ReArrange to the botom
83decb5f 1061 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
0341d0d8 1062 }
1063 pt = leadingJet.Pt() - pTback;
dbf053f0 1064 // correlation of leading jet with tracks
1065 TIterator *recIter = recParticles.MakeIterator();
dbf053f0 1066 recIter->Reset();
225f0094 1067 AliVParticle *tmpRecTrack = 0;
dbf053f0 1068 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1069 Float_t tmpPt = tmpRecTrack->Pt();
1070 // correlation
1071 Float_t tmpPhi = tmpRecTrack->Phi();
1072 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1073 Float_t dPhi = phi - tmpPhi;
1074 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1075 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1076 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1077 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
aae7993b 1078 if(cenClass>=0){
1079 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1080 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1081 }
1082
54424110 1083 }
1084
1085
d73e1768 1086 TLorentzVector vecareab;
4096256e 1087 for(int j = 0; j < nRec;j++){
1088 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1089 aodOutJet = 0;
1090 nAodOutTracks = 0;
1091 Float_t tmpPt = tmpRec.Pt();
fb10c4b8 1092
8f85b773 1093 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1094 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
d508795b 1095 aodOutJet->GetRefTracks()->Clear();
fb10c4b8 1096 Double_t area1 = clustSeq.area(sortedJets[j]);
1097 aodOutJet->SetEffArea(area1,0);
d73e1768 1098 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1099 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1100 aodOutJet->SetVectorAreaCharged(&vecareab);
fb10c4b8 1101 }
1102
1103
4096256e 1104 Float_t tmpPtBack = 0;
1105 if(externalBackground){
1106 // carefull has to be filled in a task before
2b018f7a 1107 // todo, ReArrange to the botom
4096256e 1108 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1109 }
1110 tmpPt = tmpPt - tmpPtBack;
1111 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1112
1113 fh1PtJetsRecIn->Fill(tmpPt);
ea950747 1114 // Fill Spectra with constituentsemacs
fe80c230 1115 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
fb10c4b8 1116
4096256e 1117 fh1NConstRec->Fill(constituents.size());
1118 fh2PtNch->Fill(nCh,tmpPt);
1119 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1120 fh2NConstPt->Fill(tmpPt,constituents.size());
1121 // loop over constiutents and fill spectrum
d508795b 1122
d7396b04 1123 AliVParticle *partLead = 0;
1124 Float_t ptLead = -1;
1125
d508795b 1126 for(unsigned int ic = 0; ic < constituents.size();ic++){
1127 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
006b2a30 1128 if(!part) continue;
d508795b 1129 fh1PtJetConstRec->Fill(part->Pt());
1130 if(aodOutJet){
006b2a30 1131 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
d7396b04 1132 if(part->Pt()>fMaxTrackPtInJet){
1133 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1134 }
1135 }
1136 if(part->Pt()>ptLead){
1137 partLead = part;
d508795b 1138 }
1139 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1140 }
17fdf3ee 1141 //set pT of leading constituent of jet
bc4b0d28 1142 aodOutJet->SetPtLeading(partLead->Pt());
d7396b04 1143
1144 AliAODTrack *aodT = 0;
1145 if(partLead){
b01d3ab7 1146 aodT = dynamic_cast<AliAODTrack*>(partLead);
1147 if(aodT){
d7396b04 1148 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1149 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1150 }
1151 }
1152 }
d508795b 1153
54424110 1154 // correlation
1155 Float_t tmpPhi = tmpRec.Phi();
1156 Float_t tmpEta = tmpRec.Eta();
c4856527 1157 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
54424110 1158 if(j==0){
1159 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1160 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1161 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1162 fh1NConstLeadingRec->Fill(constituents.size());
1163 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1164 continue;
1165 }
1166 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1167 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1168 Float_t dPhi = phi - tmpPhi;
1169 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1170 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1171 Float_t dEta = eta - tmpRec.Eta();
1172 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1173 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
aae7993b 1174 if(cenClass>=0){
1175 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1176 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1177 }
54424110 1178 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
8f85b773 1179 }// loop over reconstructed jets
54424110 1180 delete recIter;
8f85b773 1181
1182
1183
1184 // Add the random cones...
1185 if(fNRandomCones>0&&fTCARandomConesOut){
1186 // create a random jet within the acceptance
46465e39 1187 Double_t etaMax = fTrackEtaWindow - fRparam;
8f85b773 1188 Int_t nCone = 0;
1189 Int_t nConeRan = 0;
1190 Double_t pTC = 1; // small number
1191 for(int ir = 0;ir < fNRandomCones;ir++){
1192 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1193 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1194 // massless jet
1195 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1196 Double_t pZC = pTC/TMath::Tan(thetaC);
1197 Double_t pXC = pTC * TMath::Cos(phiC);
1198 Double_t pYC = pTC * TMath::Sin(phiC);
1199 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1200 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1201 bool skip = false;
af399bc2 1202 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
8f85b773 1203 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1204 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1205 skip = true;
1206 break;
1207 }
1208 }
1209 // test for overlap with previous cones to avoid double counting
1210 for(int iic = 0;iic<ir;iic++){
1211 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1212 if(iicone){
1213 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1214 skip = true;
1215 break;
1216 }
1217 }
1218 }
1219 if(skip)continue;
1220 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
82ac956f 1221 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1222 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
8f85b773 1223 }// loop over random cones creation
1224
1225
1226 // loop over the reconstructed particles and add up the pT in the random cones
1227 // maybe better to loop over randomized particles not in the real jets...
1228 // but this by definition brings dow average energy in the whole event
1229 AliAODJet vTmpRanR(1,0,0,1);
1230 for(int i = 0; i < recParticles.GetEntries(); i++){
1231 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1232 if(fTCARandomConesOut){
1233 for(int ir = 0;ir < fNRandomCones;ir++){
1234 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1235 if(jC&&jC->DeltaR(vp)<fRparam){
65c43de8 1236 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
8f85b773 1237 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1238 }
1239 }
1240 }// add up energy in cone
1241
1242 // the randomized input changes eta and phi, but keeps the p_T
1243 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1244 Double_t pTR = vp->Pt();
8a320ab9 1245 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
8f85b773 1246 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1247
1248 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1249 Double_t pZR = pTR/TMath::Tan(thetaR);
1250
1251 Double_t pXR = pTR * TMath::Cos(phiR);
1252 Double_t pYR = pTR * TMath::Sin(phiR);
1253 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1254 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1255 if(fTCARandomConesOutRan){
1256 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1257 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1258 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
65c43de8 1259 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
8f85b773 1260 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1261 }
1262 }
1263 }
1264 }
1265 }// loop over recparticles
1266
1267 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1268 if(fTCARandomConesOut){
1269 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1270 // rescale the momntum vectors for the random cones
1271
1272 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1273 if(rC){
1274 Double_t etaC = rC->Eta();
1275 Double_t phiC = rC->Phi();
1276 // massless jet, unit vector
1277 pTC = rC->ChargedBgEnergy();
c48ad116 1278 if(pTC<=0)pTC = 0.001; // for almost empty events
8f85b773 1279 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1280 Double_t pZC = pTC/TMath::Tan(thetaC);
1281 Double_t pXC = pTC * TMath::Cos(phiC);
1282 Double_t pYC = pTC * TMath::Sin(phiC);
1283 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1284 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1285 rC->SetBgEnergy(0,0);
1286 rC->SetEffArea(jetArea,0);
1287 }
1288 }
1289 }
82ac956f 1290 if(fTCARandomConesOutRan){
8f85b773 1291 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1292 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1293 // same wit random
1294 if(rC){
1295 Double_t etaC = rC->Eta();
1296 Double_t phiC = rC->Phi();
1297 // massless jet, unit vector
1298 pTC = rC->ChargedBgEnergy();
c48ad116 1299 if(pTC<=0)pTC = 0.001;// for almost empty events
8f85b773 1300 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1301 Double_t pZC = pTC/TMath::Tan(thetaC);
1302 Double_t pXC = pTC * TMath::Cos(phiC);
1303 Double_t pYC = pTC * TMath::Sin(phiC);
1304 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1305 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1306 rC->SetBgEnergy(0,0);
1307 rC->SetEffArea(jetArea,0);
1308 }
1309 }
1310 }
1311 }// if(fNRandomCones
d7ca0e14 1312
1313 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
d6e66a82 1314
4096256e 1315
8f85b773 1316
1317
1318
1319 if(fAODJetBackgroundOut){
d6e66a82 1320 vector<fastjet::PseudoJet> jets2=sortedJets;
1321 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1322 Double_t bkg1=0;
1323 Double_t sigma1=0.;
1324 Double_t meanarea1=0.;
1325 Double_t bkg2=0;
1326 Double_t sigma2=0.;
1327 Double_t meanarea2=0.;
82a35b1a 1328
d73e1768 1329 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
8f85b773 1330 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 1331
1332 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1333 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 1334
4a33f620 1335 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
8f85b773 1336 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
f974a156 1337 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1338 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1339
d6e66a82 1340 }
54424110 1341 }
d7ca0e14 1342
82a35b1a 1343
1344
1345
dbf053f0 1346
82a35b1a 1347 // fill track information
1348 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1349 // do the same for tracks and jets
82a35b1a 1350
1351 if(nTrackOver>0){
dbf053f0 1352 TIterator *recIter = recParticles.MakeIterator();
1353 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1354 Float_t pt = tmpRec->Pt();
82a35b1a 1355
dbf053f0 1356 // Printf("Leading track p_t %3.3E",pt);
1357 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1358 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1359 while(pt<ptCut&&tmpRec){
1360 nTrackOver--;
1361 tmpRec = (AliVParticle*)(recIter->Next());
1362 if(tmpRec){
1363 pt = tmpRec->Pt();
1364 }
1365 }
1366 if(nTrackOver<=0)break;
1367 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1368 }
1369
1370 recIter->Reset();
1371 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1372 Float_t phi = leading->Phi();
1373 if(phi<0)phi+=TMath::Pi()*2.;
1374 Float_t eta = leading->Eta();
1375 pt = leading->Pt();
1376 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1377 Float_t tmpPt = tmpRec->Pt();
1378 Float_t tmpEta = tmpRec->Eta();
1379 fh1PtTracksRecIn->Fill(tmpPt);
1380 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1381 if(tmpRec==leading){
1382 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1383 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1384 continue;
1385 }
1386 // correlation
1387 Float_t tmpPhi = tmpRec->Phi();
1388
1389 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1390 Float_t dPhi = phi - tmpPhi;
1391 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1392 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1393 Float_t dEta = eta - tmpRec->Eta();
1394 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1395 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1396 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1397 }
1398 delete recIter;
1399 }
1400
1401 // find the random jets
fe80c230 1402
d6e66a82 1403 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1404
dbf053f0 1405 // fill the jet information from random track
fe80c230 1406 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1407 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
dbf053f0 1408
1409 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1410
1411 // loop over all jets an fill information, first on is the leading jet
1412
1413 Int_t nRecOverRan = inclusiveJetsRan.size();
1414 Int_t nRecRan = inclusiveJetsRan.size();
d73e1768 1415
dbf053f0 1416 if(inclusiveJetsRan.size()>0){
1417 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1418 Float_t pt = leadingJet.Pt();
1419
1420 Int_t iCount = 0;
d73e1768 1421 TLorentzVector vecarearanb;
1422
dbf053f0 1423 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1424 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1425 while(pt<ptCut&&iCount<nRecRan){
1426 nRecOverRan--;
1427 iCount++;
1428 if(iCount<nRecRan){
1429 pt = sortedJetsRan[iCount].perp();
1430 }
1431 }
1432 if(nRecOverRan<=0)break;
1433 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1434 }
1435 Float_t phi = leadingJet.Phi();
1436 if(phi<0)phi+=TMath::Pi()*2.;
1437 pt = leadingJet.Pt();
1438
1439 // correlation of leading jet with random tracks
1440
1441 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1442 {
1443 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1444 // correlation
1445 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1446 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1447 Float_t dPhi = phi - tmpPhi;
1448 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1449 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1450 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1451 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1452 }
1453
d6e66a82 1454 Int_t nAodOutJetsRan = 0;
1455 AliAODJet *aodOutJetRan = 0;
dbf053f0 1456 for(int j = 0; j < nRecRan;j++){
1457 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1458 Float_t tmpPt = tmpRec.Pt();
1459 fh1PtJetsRecInRan->Fill(tmpPt);
1460 // Fill Spectra with constituents
fe80c230 1461 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
dbf053f0 1462 fh1NConstRecRan->Fill(constituents.size());
1463 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1464 fh2PtNchRan->Fill(nCh,tmpPt);
1465 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1466
1467
8f85b773 1468 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1469 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
d6e66a82 1470 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
d508795b 1471 aodOutJetRan->GetRefTracks()->Clear();
d73e1768 1472 aodOutJetRan->SetEffArea(arearan,0);
1473 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1474 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1475 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1476
1477 }
d6e66a82 1478
dbf053f0 1479 // correlation
1480 Float_t tmpPhi = tmpRec.Phi();
1481 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1482
1483 if(j==0){
1484 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1485 fh1NConstLeadingRecRan->Fill(constituents.size());
1486 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1487 continue;
1488 }
1489 }
d6e66a82 1490
1491
8f85b773 1492 if(fAODJetBackgroundOut){
d6e66a82 1493 Double_t bkg3=0.;
1494 Double_t sigma3=0.;
1495 Double_t meanarea3=0.;
4a33f620 1496 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
8f85b773 1497 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1498 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1499 /*
82a35b1a 1500 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1501 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1502 */
d6e66a82 1503 }
1504
1505
1506
dbf053f0 1507 }
0341d0d8 1508
1509
1510 // do the event selection if activated
1511 if(fJetTriggerPtCut>0){
1512 bool select = false;
2b018f7a 1513 Float_t minPt = fJetTriggerPtCut;
1514 /*
0341d0d8 1515 // hard coded for now ...
2b018f7a 1516 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
0341d0d8 1517 if(cent<10)minPt = 50;
aae7993b 1518 else if(cent<30)minPt = 42;
0341d0d8 1519 else if(cent<50)minPt = 28;
1520 else if(cent<80)minPt = 18;
2b018f7a 1521 */
0341d0d8 1522 float rho = 0;
1523 if(externalBackground)rho = externalBackground->GetBackground(2);
8f85b773 1524 if(fTCAJetsOut){
1525 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1526 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
0341d0d8 1527 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1528 if(ptSub>=minPt){
1529 select = true;
1530 break;
1531 }
1532 }
1533 }
ea950747 1534
0341d0d8 1535 if(select){
1536 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1537 fh1CentralitySelect->Fill(cent);
e2ca7519 1538 fh1ZSelect->Fill(zVtx);
0341d0d8 1539 aodH->SetFillAOD(kTRUE);
1540 }
1541 }
ea950747 1542 if (fDebug > 2){
1543 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1544 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1545 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1546 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1547 }
dbf053f0 1548 PostData(1, fHistList);
1549}
1550
1551void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1552{
f0659f11 1553 //
1554 // Terminate analysis
1555 //
006b2a30 1556 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1557
1558 if(fMomResH1Fit) delete fMomResH1Fit;
1559 if(fMomResH2Fit) delete fMomResH2Fit;
1560 if(fMomResH3Fit) delete fMomResH3Fit;
1561
dbf053f0 1562}
1563
1564
1565Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1566
f0659f11 1567 //
1568 // get list of tracks/particles for different types
1569 //
1570
dbf053f0 1571 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1572
1573 Int_t iCount = 0;
f840d64c 1574 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1575 if(type!=kTrackAODextraonly) {
1576 AliAODEvent *aod = 0;
1577 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1578 else aod = AODEvent();
1579 if(!aod){
617932d0 1580 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
f840d64c 1581 return iCount;
1582 }
d7396b04 1583
f840d64c 1584 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1585 AliAODTrack *tr = aod->GetTrack(it);
b2c1e955 1586 Bool_t bGood = false;
1587 if(fFilterType == 0)bGood = true;
8a0f6112 1588 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1589 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
b2c1e955 1590 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
617932d0 1591 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1592 continue;
1593 }
1594 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1595 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1596 continue;
1597 }
1598 if(tr->Pt()<fTrackPtCut){
1599 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1600 continue;
1601 }
1602 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
f840d64c 1603 list->Add(tr);
1604 iCount++;
1605 }
dbf053f0 1606 }
f840d64c 1607 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1608 AliAODEvent *aod = 0;
1609 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1610 else aod = AODEvent();
1611
1612 if(!aod){
1613 return iCount;
1614 }
1615 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
8ddcf605 1616 if(!aodExtraTracks)return iCount;
f840d64c 1617 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1618 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1619 if (!track) continue;
1620
1621 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
8ddcf605 1622 if(!trackAOD)continue;
b2c1e955 1623 Bool_t bGood = false;
1624 if(fFilterType == 0)bGood = true;
8a0f6112 1625 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1626 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
b2c1e955 1627 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
ca6d12f9 1628 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
f840d64c 1629 if(trackAOD->Pt()<fTrackPtCut) continue;
f840d64c 1630 list->Add(trackAOD);
1631 iCount++;
1632 }
dbf053f0 1633 }
1634 }
1635 else if (type == kTrackKineAll||type == kTrackKineCharged){
1636 AliMCEvent* mcEvent = MCEvent();
1637 if(!mcEvent)return iCount;
1638 // we want to have alivpartilces so use get track
1639 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1640 if(!mcEvent->IsPhysicalPrimary(it))continue;
1641 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1642 if(type == kTrackKineAll){
0341d0d8 1643 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1644 list->Add(part);
1645 iCount++;
1646 }
1647 else if(type == kTrackKineCharged){
1648 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 1649 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1650 list->Add(part);
1651 iCount++;
1652 }
1653 }
1654 }
1655 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1656 AliAODEvent *aod = 0;
959508f4 1657 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1658 else aod = AODEvent();
1659 if(!aod)return iCount;
1660 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1661 if(!tca)return iCount;
1662 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1663 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1664 if(!part->IsPhysicalPrimary())continue;
1665 if(type == kTrackAODMCAll){
0341d0d8 1666 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1667 list->Add(part);
1668 iCount++;
1669 }
1670 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1671 if(part->Charge()==0)continue;
0341d0d8 1672 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1673 if(kTrackAODMCCharged){
1674 list->Add(part);
1675 }
1676 else {
8a320ab9 1677 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
dbf053f0 1678 list->Add(part);
1679 }
1680 iCount++;
1681 }
1682 }
1683 }// AODMCparticle
1684 list->Sort();
1685 return iCount;
dbf053f0 1686}
1687
006b2a30 1688void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1689
1690 //
1691 // set mom res profiles
1692 //
1693
1694 fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
1695 fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
1696 fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
1697}
1698
1699void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1700 //
1701 // set tracking efficiency histos
1702 //
1703
1704 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1705 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1706 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1707}
1708
1709Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1710
1711 //
1712 // Get smearing on generated momentum
1713 //
1714
1715 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1716
1717 TProfile *fMomRes = 0x0;
1718 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1719 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1720 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1721
1722 if(!fMomRes) {
1723 return 0.;
1724 }
1725
1726
1727 Double_t smear = 0.;
1728
1729 if(pt>20.) {
1730 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1731 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1732 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1733 }
1734 else {
1735
1736 Int_t bin = fMomRes->FindBin(pt);
1737
1738 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1739
1740 }
1741
1742 if(fMomRes) delete fMomRes;
1743
1744 return smear;
1745}
1746
1747void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1748 //
1749 // Fit linear function on momentum resolution at high pT
1750 //
1751
1752 if(!fMomResH1Fit && fMomResH1) {
1753 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1754 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1755 fMomResH1Fit ->SetRange(5.,100.);
1756 }
1757
1758 if(!fMomResH2Fit && fMomResH2) {
1759 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1760 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1761 fMomResH2Fit ->SetRange(5.,100.);
1762 }
1763
1764 if(!fMomResH3Fit && fMomResH3) {
1765 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1766 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1767 fMomResH3Fit ->SetRange(5.,100.);
1768 }
1769
1770}
1771
dbf053f0 1772/*
1773Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1774 for(int i = 0; i < particles.GetEntries(); i++){
1775 AliVParticle *vp = (AliVParticle*)particles.At(i);
1776 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1777 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1778 jInp.set_user_index(i);
1779 inputParticles.push_back(jInp);
1780 }
1781
1782 return 0;
1783
1784}
1785*/