]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/DEV/AliAnalysisTaskJetCluster.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Compatibility with the ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / DEV / AliAnalysisTaskJetCluster.cxx
CommitLineData
d89b8229 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>
24#include <TRandom3.h>
25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
28#include <TRefArray.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>
006b2a30 35#include <TF1.h>
d89b8229 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"
46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
49#include "AliAODExtension.h"
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"
60#include "AliAODJetEventBackground.h"
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
72AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
f0659f11 73 //
74 // Destructor
75 //
006b2a30 76
d89b8229 77 delete fRef;
78 delete fRandom;
79
80 if(fTCAJetsOut)fTCAJetsOut->Delete();
81 delete fTCAJetsOut;
006b2a30 82
d89b8229 83 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
84 delete fTCAJetsOutRan;
006b2a30 85
d89b8229 86 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
87 delete fTCARandomConesOut;
006b2a30 88
d89b8229 89 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
90 delete fTCARandomConesOutRan;
006b2a30 91
d89b8229 92 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
93 delete fAODJetBackgroundOut;
94}
95
96AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
97 AliAnalysisTaskSE(),
98 fAOD(0x0),
99 fAODExtension(0x0),
100 fRef(new TRefArray),
101 fUseAODTrackInput(kFALSE),
102 fUseAODMCInput(kFALSE),
103 fUseBackgroundCalc(kFALSE),
104 fEventSelection(kFALSE),
105 fFilterMask(0),
d7396b04 106 fFilterMaskBestPt(0),
d89b8229 107 fFilterType(0),
108 fJetTypes(kJet),
109 fTrackTypeRec(kTrackUndef),
110 fTrackTypeGen(kTrackUndef),
111 fNSkipLeadingRan(0),
112 fNSkipLeadingCone(0),
113 fNRandomCones(0),
114 fAvgTrials(1),
115 fExternalWeight(1),
116 fTrackEtaWindow(0.9),
117 fRecEtaWindow(0.5),
118 fTrackPtCut(0.),
119 fJetOutputMinPt(0.150),
120 fMaxTrackPtInJet(100.),
121 fJetTriggerPtCut(0),
122 fVtxZCut(8),
123 fVtxR2Cut(1),
124 fCentCutUp(0),
125 fCentCutLo(0),
126 fNonStdBranch(""),
127 fBackgroundBranch(""),
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),
d89b8229 140 fRparam(1.0),
141 fAlgorithm(fastjet::kt_algorithm),
142 fStrategy(fastjet::Best),
143 fRecombScheme(fastjet::BIpt_scheme),
144 fAreaType(fastjet::active_area),
145 fGhostArea(0.01),
146 fActiveAreaRepeats(1),
147 fGhostEtamax(1.5),
148 fTCAJetsOut(0x0),
149 fTCAJetsOutRan(0x0),
150 fTCARandomConesOut(0x0),
151 fTCARandomConesOutRan(0x0),
152 fAODJetBackgroundOut(0x0),
153 fRandom(0),
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),
174 fh1Nch(0x0),
175 fh1CentralityPhySel(0x0),
176 fh1Centrality(0x0),
177 fh1CentralitySelect(0x0),
178 fh1ZPhySel(0x0),
179 fh1Z(0x0),
180 fh1ZSelect(0x0),
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),
202 fh2PtNch(0x0),
203 fh2PtNchRan(0x0),
204 fh2PtNchN(0x0),
205 fh2PtNchNRan(0x0),
206 fh2TracksLeadingJetPhiPtRan(0x0),
207 fh2TracksLeadingJetPhiPtWRan(0x0),
006b2a30 208 fh2PtGenPtSmeared(0x0),
209 fp1Efficiency(0x0),
210 fp1PtResolution(0x0),
d89b8229 211 fHistList(0x0)
212{
f0659f11 213 //
214 // Constructor
215 //
216
d89b8229 217 for(int i = 0;i<3;i++){
218 fh1BiARandomCones[i] = 0;
219 fh1BiARandomConesRan[i] = 0;
220 }
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 }
227}
228
229AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
230 AliAnalysisTaskSE(name),
231 fAOD(0x0),
232 fAODExtension(0x0),
233 fRef(new TRefArray),
234 fUseAODTrackInput(kFALSE),
235 fUseAODMCInput(kFALSE),
236 fUseBackgroundCalc(kFALSE),
d7396b04 237 fEventSelection(kFALSE), fFilterMask(0),
238 fFilterMaskBestPt(0),
d89b8229 239 fFilterType(0),
240 fJetTypes(kJet),
241 fTrackTypeRec(kTrackUndef),
242 fTrackTypeGen(kTrackUndef),
243 fNSkipLeadingRan(0),
244 fNSkipLeadingCone(0),
245 fNRandomCones(0),
246 fAvgTrials(1),
247 fExternalWeight(1),
248 fTrackEtaWindow(0.9),
249 fRecEtaWindow(0.5),
250 fTrackPtCut(0.),
251 fJetOutputMinPt(0.150),
252 fMaxTrackPtInJet(100.),
253 fJetTriggerPtCut(0),
254 fVtxZCut(8),
255 fVtxR2Cut(1),
256 fCentCutUp(0),
257 fCentCutLo(0),
258 fNonStdBranch(""),
259 fBackgroundBranch(""),
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),
d89b8229 272 fRparam(1.0),
273 fAlgorithm(fastjet::kt_algorithm),
274 fStrategy(fastjet::Best),
275 fRecombScheme(fastjet::BIpt_scheme),
276 fAreaType(fastjet::active_area),
277 fGhostArea(0.01),
278 fActiveAreaRepeats(1),
279 fGhostEtamax(1.5),
280 fTCAJetsOut(0x0),
281 fTCAJetsOutRan(0x0),
282 fTCARandomConesOut(0x0),
283 fTCARandomConesOutRan(0x0),
284 fAODJetBackgroundOut(0x0),
285 fRandom(0),
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),
306 fh1Nch(0x0),
307 fh1CentralityPhySel(0x0),
308 fh1Centrality(0x0),
309 fh1CentralitySelect(0x0),
310 fh1ZPhySel(0x0),
311 fh1Z(0x0),
312 fh1ZSelect(0x0),
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),
334 fh2PtNch(0x0),
335 fh2PtNchRan(0x0),
336 fh2PtNchN(0x0),
337 fh2PtNchNRan(0x0),
338 fh2TracksLeadingJetPhiPtRan(0x0),
339 fh2TracksLeadingJetPhiPtWRan(0x0),
006b2a30 340 fh2PtGenPtSmeared(0x0),
341 fp1Efficiency(0x0),
342 fp1PtResolution(0x0),
d89b8229 343 fHistList(0x0)
344{
f0659f11 345 //
346 // named ctor
347 //
006b2a30 348
d89b8229 349 for(int i = 0;i<3;i++){
350 fh1BiARandomCones[i] = 0;
351 fh1BiARandomConesRan[i] = 0;
352 }
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 }
359 DefineOutput(1, TList::Class());
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
380 fRandom = new TRandom3(0);
381
382
383 // Connect the AOD
384
385
386 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
387
388
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
396
397 if(fJetTypes&kJet){
398 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
399 fTCAJetsOut->SetName(fNonStdBranch.Data());
400 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
401 }
006b2a30 402
d89b8229 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));
d89b8229 408 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
409 }
410
411 if(fUseBackgroundCalc){
412 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
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
d89b8229 418 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
419 }
420 }
421 // create the branch for the random cones with the same R
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
d89b8229 426 if(fNRandomCones>0){
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 }
433 }
434 // create the branch with the random for the random cones on the random event
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 }
442 }
443 }
444
445 if(fNonStdFile.Length()!=0){
446 //
447 // case that we have an AOD extension we need to fetch the jets from the extended output
448 // we identify the extension aod event by looking for the branchname
449 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
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);
452 }
453 }
454
006b2a30 455 // FitMomentumResolution();
456
d89b8229 457
458 if(!fHistList)fHistList = new TList();
459 fHistList->SetOwner();
460 PostData(1, fHistList); // post data in any case once
461
462 Bool_t oldStatus = TH1::AddDirectoryStatus();
463 TH1::AddDirectory(kFALSE);
464
465 //
466 // Histogram
467
468 const Int_t nBinPt = 100;
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
475 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
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
503 const int nChMax = 5000;
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);
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);
534 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
535
536 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
537 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
538 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
539
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
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
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
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);
d89b8229 611
612 if(fNRandomCones>0&&fUseBackgroundCalc){
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 }
618
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 }
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);
645 fHistList->Add(fh1Nch);
646 fHistList->Add(fh1Centrality);
647 fHistList->Add(fh1CentralitySelect);
648 fHistList->Add(fh1CentralityPhySel);
649 fHistList->Add(fh1Z);
650 fHistList->Add(fh1ZSelect);
651 fHistList->Add(fh1ZPhySel);
652 if(fNRandomCones>0&&fUseBackgroundCalc){
653 for(int i = 0;i<3;i++){
654 fHistList->Add(fh1BiARandomCones[i]);
655 fHistList->Add(fh1BiARandomConesRan[i]);
656 }
657 }
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
665 fHistList->Add(fh2NRecJetsPt);
666 fHistList->Add(fh2NRecTracksPt);
667 fHistList->Add(fh2NConstPt);
668 fHistList->Add(fh2NConstLeadingPt);
669 fHistList->Add(fh2PtNch);
670 fHistList->Add(fh2PtNchRan);
671 fHistList->Add(fh2PtNchN);
672 fHistList->Add(fh2PtNchNRan);
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);
692 }
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
713 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
714
006b2a30 715 FitMomentumResolution();
716
d89b8229 717}
718
719void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
720{
721
722 // handle and reset the output jet branch
723
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();
729
730 AliAODJetEventBackground* externalBackground = 0;
731 if(!externalBackground&&fBackgroundBranch.Length()){
732 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
733 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
734 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
735 }
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
d89b8229 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;
d89b8229 763
764 Bool_t selectEvent = false;
765 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
766
767 Float_t cent = 0;
768 Float_t zVtx = 0;
769 Int_t cenClass = -1;
770 if(fAOD){
771 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
772 TString vtxTitle(vtxAOD->GetTitle());
773 zVtx = vtxAOD->GetZ();
774
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;
780 if(physicsSelection){
781 fh1CentralityPhySel->Fill(cent);
782 fh1ZPhySel->Fill(zVtx);
783 }
784
785 if(fEventSelection){
786 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
787 Float_t yvtx = vtxAOD->GetY();
788 Float_t xvtx = vtxAOD->GetX();
789 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
790 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
791 if(physicsSelection){
792 selectEvent = true;
793 }
794 }
795 }
796 if(fCentCutUp>0){
797 if(cent<fCentCutLo||cent>fCentCutUp){
798 selectEvent = false;
799 }
800 }
801 }else{
802 selectEvent = true;
803 }
804 }
805
806
807 if(!selectEvent){
808 PostData(1, fHistList);
809 return;
810 }
811 fh1Centrality->Fill(cent);
812 fh1Z->Fill(zVtx);
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);
828 Float_t nCh = recParticles.GetEntries();
829 fh1Nch->Fill(nCh);
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;
838
839 // Generate the random cones
840
841 AliAODJet vTmpRan(1,0,0,1);
842 for(int i = 0; i < recParticles.GetEntries(); i++){
843 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
006b2a30 844
d89b8229 845 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
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 }
d89b8229 956
957 // the randomized input changes eta and phi, but keeps the p_T
958 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
959 Double_t pT = vp->Pt();
960 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
961 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
962
963 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
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);
970
971 jInpRan.set_user_index(i);
972 inputParticlesRecRan.push_back(jInpRan);
973 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
974 }
975
976 // fill the tref array, only needed when we write out jets
977 if(fTCAJetsOut){
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 }
d89b8229 983 }
006b2a30 984 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
d89b8229 985 }
986 }// recparticles
987
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
995 // employ setters for these...
996
997
998 // now create the object that holds info about ghosts
999 /*
1000 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1001 // reduce CPU time...
1002 fGhostArea = 0.5;
1003 fActiveAreaRepeats = 0;
1004 }
1005 */
1006
1007 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1008 fastjet::AreaType areaType = fastjet::active_area;
1009 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1010 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1011 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1012
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
1022 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1023 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1024
1025
1026 fh1NJetsRec->Fill(sortedJets.size());
1027
1028 // loop over all jets an fill information, first on is the leading jet
1029
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());
1034 Double_t area = clustSeq.area(sortedJets[0]);
1035 leadingJet.SetEffArea(area,0);
1036 Float_t pt = leadingJet.Pt();
1037 Int_t nAodOutJets = 0;
1038 Int_t nAodOutTracks = 0;
1039 AliAODJet *aodOutJet = 0;
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();
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();
d89b8229 1062 }
1063 pt = leadingJet.Pt() - pTback;
1064 // correlation of leading jet with tracks
1065 TIterator *recIter = recParticles.MakeIterator();
1066 recIter->Reset();
1067 AliVParticle *tmpRecTrack = 0;
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();
1076 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1077 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1078 if(cenClass>=0){
1079 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1080 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1081 }
1082
1083 }
1084
1085
1086 TLorentzVector vecareab;
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();
1092
1093 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1094 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1095 aodOutJet->GetRefTracks()->Clear();
1096 Double_t area1 = clustSeq.area(sortedJets[j]);
1097 aodOutJet->SetEffArea(area1,0);
1098 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1099 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1100 aodOutJet->SetVectorAreaCharged(&vecareab);
1101 }
1102
1103
1104 Float_t tmpPtBack = 0;
1105 if(externalBackground){
1106 // carefull has to be filled in a task before
1107 // todo, ReArrange to the botom
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);
1114 // Fill Spectra with constituentsemacs
1115 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1116
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
1122
d7396b04 1123 AliVParticle *partLead = 0;
1124 Float_t ptLead = -1;
1125
d89b8229 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;
d89b8229 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;
d89b8229 1138 }
1139 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1140 }
d7396b04 1141
1142 AliAODTrack *aodT = 0;
1143 if(partLead){
b01d3ab7 1144 aodT = dynamic_cast<AliAODTrack*>(partLead);
1145 if(aodT){
d7396b04 1146 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1147 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1148 }
1149 }
1150 }
d89b8229 1151
1152 // correlation
1153 Float_t tmpPhi = tmpRec.Phi();
1154 Float_t tmpEta = tmpRec.Eta();
1155 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1156 if(j==0){
1157 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1158 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1159 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1160 fh1NConstLeadingRec->Fill(constituents.size());
1161 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1162 continue;
1163 }
1164 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1165 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1166 Float_t dPhi = phi - tmpPhi;
1167 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1168 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1169 Float_t dEta = eta - tmpRec.Eta();
1170 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1171 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1172 if(cenClass>=0){
1173 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1174 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1175 }
1176 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1177 }// loop over reconstructed jets
1178 delete recIter;
1179
1180
1181
1182 // Add the random cones...
1183 if(fNRandomCones>0&&fTCARandomConesOut){
1184 // create a random jet within the acceptance
1185 Double_t etaMax = fTrackEtaWindow - fRparam;
1186 Int_t nCone = 0;
1187 Int_t nConeRan = 0;
1188 Double_t pTC = 1; // small number
1189 for(int ir = 0;ir < fNRandomCones;ir++){
1190 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1191 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1192 // massless jet
1193 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1194 Double_t pZC = pTC/TMath::Tan(thetaC);
1195 Double_t pXC = pTC * TMath::Cos(phiC);
1196 Double_t pYC = pTC * TMath::Sin(phiC);
1197 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1198 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1199 bool skip = false;
1200 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1201 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1202 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1203 skip = true;
1204 break;
1205 }
1206 }
1207 // test for overlap with previous cones to avoid double counting
1208 for(int iic = 0;iic<ir;iic++){
1209 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1210 if(iicone){
1211 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1212 skip = true;
1213 break;
1214 }
1215 }
1216 }
1217 if(skip)continue;
1218 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1219 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1220 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1221 }// loop over random cones creation
1222
1223
1224 // loop over the reconstructed particles and add up the pT in the random cones
1225 // maybe better to loop over randomized particles not in the real jets...
1226 // but this by definition brings dow average energy in the whole event
1227 AliAODJet vTmpRanR(1,0,0,1);
1228 for(int i = 0; i < recParticles.GetEntries(); i++){
1229 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1230 if(fTCARandomConesOut){
1231 for(int ir = 0;ir < fNRandomCones;ir++){
1232 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1233 if(jC&&jC->DeltaR(vp)<fRparam){
1234 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1235 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1236 }
1237 }
1238 }// add up energy in cone
1239
1240 // the randomized input changes eta and phi, but keeps the p_T
1241 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1242 Double_t pTR = vp->Pt();
1243 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1244 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1245
1246 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1247 Double_t pZR = pTR/TMath::Tan(thetaR);
1248
1249 Double_t pXR = pTR * TMath::Cos(phiR);
1250 Double_t pYR = pTR * TMath::Sin(phiR);
1251 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1252 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1253 if(fTCARandomConesOutRan){
1254 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1255 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1256 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1257 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1258 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1259 }
1260 }
1261 }
1262 }
1263 }// loop over recparticles
1264
1265 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1266 if(fTCARandomConesOut){
1267 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1268 // rescale the momntum vectors for the random cones
1269
1270 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1271 if(rC){
1272 Double_t etaC = rC->Eta();
1273 Double_t phiC = rC->Phi();
1274 // massless jet, unit vector
1275 pTC = rC->ChargedBgEnergy();
1276 if(pTC<=0)pTC = 0.001; // for almost empty events
1277 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1278 Double_t pZC = pTC/TMath::Tan(thetaC);
1279 Double_t pXC = pTC * TMath::Cos(phiC);
1280 Double_t pYC = pTC * TMath::Sin(phiC);
1281 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1282 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1283 rC->SetBgEnergy(0,0);
1284 rC->SetEffArea(jetArea,0);
1285 }
1286 }
1287 }
1288 if(fTCARandomConesOutRan){
1289 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1290 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1291 // same wit random
1292 if(rC){
1293 Double_t etaC = rC->Eta();
1294 Double_t phiC = rC->Phi();
1295 // massless jet, unit vector
1296 pTC = rC->ChargedBgEnergy();
1297 if(pTC<=0)pTC = 0.001;// for almost empty events
1298 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1299 Double_t pZC = pTC/TMath::Tan(thetaC);
1300 Double_t pXC = pTC * TMath::Cos(phiC);
1301 Double_t pYC = pTC * TMath::Sin(phiC);
1302 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1303 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1304 rC->SetBgEnergy(0,0);
1305 rC->SetEffArea(jetArea,0);
1306 }
1307 }
1308 }
1309 }// if(fNRandomCones
1310
1311 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1312
1313
1314
1315
1316
1317 if(fAODJetBackgroundOut){
1318 vector<fastjet::PseudoJet> jets2=sortedJets;
1319 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1320 Double_t bkg1=0;
1321 Double_t sigma1=0.;
1322 Double_t meanarea1=0.;
1323 Double_t bkg2=0;
1324 Double_t sigma2=0.;
1325 Double_t meanarea2=0.;
1326
1327 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1328 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1329
1330 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1331 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1332
1333 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1334 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1335 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1336 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1337
1338 }
1339 }
1340
1341
1342
1343
1344
1345 // fill track information
1346 Int_t nTrackOver = recParticles.GetSize();
1347 // do the same for tracks and jets
1348
1349 if(nTrackOver>0){
1350 TIterator *recIter = recParticles.MakeIterator();
1351 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1352 Float_t pt = tmpRec->Pt();
1353
1354 // Printf("Leading track p_t %3.3E",pt);
1355 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1356 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1357 while(pt<ptCut&&tmpRec){
1358 nTrackOver--;
1359 tmpRec = (AliVParticle*)(recIter->Next());
1360 if(tmpRec){
1361 pt = tmpRec->Pt();
1362 }
1363 }
1364 if(nTrackOver<=0)break;
1365 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1366 }
1367
1368 recIter->Reset();
1369 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1370 Float_t phi = leading->Phi();
1371 if(phi<0)phi+=TMath::Pi()*2.;
1372 Float_t eta = leading->Eta();
1373 pt = leading->Pt();
1374 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1375 Float_t tmpPt = tmpRec->Pt();
1376 Float_t tmpEta = tmpRec->Eta();
1377 fh1PtTracksRecIn->Fill(tmpPt);
1378 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1379 if(tmpRec==leading){
1380 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1381 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1382 continue;
1383 }
1384 // correlation
1385 Float_t tmpPhi = tmpRec->Phi();
1386
1387 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1388 Float_t dPhi = phi - tmpPhi;
1389 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1390 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1391 Float_t dEta = eta - tmpRec->Eta();
1392 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1393 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1394 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1395 }
1396 delete recIter;
1397 }
1398
1399 // find the random jets
1400
1401 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1402
1403 // fill the jet information from random track
1404 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1405 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1406
1407 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1408
1409 // loop over all jets an fill information, first on is the leading jet
1410
1411 Int_t nRecOverRan = inclusiveJetsRan.size();
1412 Int_t nRecRan = inclusiveJetsRan.size();
1413
1414 if(inclusiveJetsRan.size()>0){
1415 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1416 Float_t pt = leadingJet.Pt();
1417
1418 Int_t iCount = 0;
1419 TLorentzVector vecarearanb;
1420
1421 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1422 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1423 while(pt<ptCut&&iCount<nRecRan){
1424 nRecOverRan--;
1425 iCount++;
1426 if(iCount<nRecRan){
1427 pt = sortedJetsRan[iCount].perp();
1428 }
1429 }
1430 if(nRecOverRan<=0)break;
1431 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1432 }
1433 Float_t phi = leadingJet.Phi();
1434 if(phi<0)phi+=TMath::Pi()*2.;
1435 pt = leadingJet.Pt();
1436
1437 // correlation of leading jet with random tracks
1438
1439 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1440 {
1441 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1442 // correlation
1443 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1444 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1445 Float_t dPhi = phi - tmpPhi;
1446 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1447 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1448 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1449 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1450 }
1451
1452 Int_t nAodOutJetsRan = 0;
1453 AliAODJet *aodOutJetRan = 0;
1454 for(int j = 0; j < nRecRan;j++){
1455 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1456 Float_t tmpPt = tmpRec.Pt();
1457 fh1PtJetsRecInRan->Fill(tmpPt);
1458 // Fill Spectra with constituents
1459 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1460 fh1NConstRecRan->Fill(constituents.size());
1461 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1462 fh2PtNchRan->Fill(nCh,tmpPt);
1463 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1464
1465
1466 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1467 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1468 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1469 aodOutJetRan->GetRefTracks()->Clear();
1470 aodOutJetRan->SetEffArea(arearan,0);
1471 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1472 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1473 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1474
1475 }
1476
1477 // correlation
1478 Float_t tmpPhi = tmpRec.Phi();
1479 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1480
1481 if(j==0){
1482 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1483 fh1NConstLeadingRecRan->Fill(constituents.size());
1484 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1485 continue;
1486 }
1487 }
1488
1489
1490 if(fAODJetBackgroundOut){
1491 Double_t bkg3=0.;
1492 Double_t sigma3=0.;
1493 Double_t meanarea3=0.;
1494 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1495 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1496 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1497 /*
1498 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1499 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1500 */
1501 }
1502
1503
1504
1505 }
1506
1507
1508 // do the event selection if activated
1509 if(fJetTriggerPtCut>0){
1510 bool select = false;
1511 Float_t minPt = fJetTriggerPtCut;
1512 /*
1513 // hard coded for now ...
1514 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1515 if(cent<10)minPt = 50;
1516 else if(cent<30)minPt = 42;
1517 else if(cent<50)minPt = 28;
1518 else if(cent<80)minPt = 18;
1519 */
1520 float rho = 0;
1521 if(externalBackground)rho = externalBackground->GetBackground(2);
1522 if(fTCAJetsOut){
1523 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1524 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1525 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1526 if(ptSub>=minPt){
1527 select = true;
1528 break;
1529 }
1530 }
1531 }
1532
1533 if(select){
1534 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1535 fh1CentralitySelect->Fill(cent);
1536 fh1ZSelect->Fill(zVtx);
1537 aodH->SetFillAOD(kTRUE);
1538 }
1539 }
1540 if (fDebug > 2){
1541 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1542 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1543 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1544 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1545 }
1546 PostData(1, fHistList);
1547}
1548
1549void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1550{
f0659f11 1551 //
1552 // Terminate analysis
1553 //
006b2a30 1554 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1555
1556 if(fMomResH1Fit) delete fMomResH1Fit;
1557 if(fMomResH2Fit) delete fMomResH2Fit;
1558 if(fMomResH3Fit) delete fMomResH3Fit;
1559
d89b8229 1560}
1561
1562
1563Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1564
f0659f11 1565 //
1566 // get list of tracks/particles for different types
1567 //
1568
d89b8229 1569 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1570
1571 Int_t iCount = 0;
1572 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1573 if(type!=kTrackAODextraonly) {
1574 AliAODEvent *aod = 0;
1575 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1576 else aod = AODEvent();
1577 if(!aod){
1578 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1579 return iCount;
1580 }
d7396b04 1581
d89b8229 1582 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1583 AliAODTrack *tr = aod->GetTrack(it);
1584 Bool_t bGood = false;
1585 if(fFilterType == 0)bGood = true;
1586 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1587 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1588 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1589 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1590 continue;
1591 }
1592 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1593 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1594 continue;
1595 }
1596 if(tr->Pt()<fTrackPtCut){
1597 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1598 continue;
1599 }
1600 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1601 list->Add(tr);
1602 iCount++;
1603 }
1604 }
1605 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1606 AliAODEvent *aod = 0;
1607 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1608 else aod = AODEvent();
1609
1610 if(!aod){
1611 return iCount;
1612 }
1613 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1614 if(!aodExtraTracks)return iCount;
1615 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1616 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1617 if (!track) continue;
1618
1619 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1620 if(!trackAOD)continue;
1621 Bool_t bGood = false;
1622 if(fFilterType == 0)bGood = true;
1623 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1624 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1625 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1626 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1627 if(trackAOD->Pt()<fTrackPtCut) continue;
1628 list->Add(trackAOD);
1629 iCount++;
1630 }
1631 }
1632 }
1633 else if (type == kTrackKineAll||type == kTrackKineCharged){
1634 AliMCEvent* mcEvent = MCEvent();
1635 if(!mcEvent)return iCount;
1636 // we want to have alivpartilces so use get track
1637 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1638 if(!mcEvent->IsPhysicalPrimary(it))continue;
1639 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1640 if(type == kTrackKineAll){
1641 if(part->Pt()<fTrackPtCut)continue;
1642 list->Add(part);
1643 iCount++;
1644 }
1645 else if(type == kTrackKineCharged){
1646 if(part->Particle()->GetPDG()->Charge()==0)continue;
1647 if(part->Pt()<fTrackPtCut)continue;
1648 list->Add(part);
1649 iCount++;
1650 }
1651 }
1652 }
1653 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1654 AliAODEvent *aod = 0;
1655 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1656 else aod = AODEvent();
1657 if(!aod)return iCount;
1658 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1659 if(!tca)return iCount;
1660 for(int it = 0;it < tca->GetEntriesFast();++it){
1661 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1662 if(!part->IsPhysicalPrimary())continue;
1663 if(type == kTrackAODMCAll){
1664 if(part->Pt()<fTrackPtCut)continue;
1665 list->Add(part);
1666 iCount++;
1667 }
1668 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1669 if(part->Charge()==0)continue;
1670 if(part->Pt()<fTrackPtCut)continue;
1671 if(kTrackAODMCCharged){
1672 list->Add(part);
1673 }
1674 else {
1675 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1676 list->Add(part);
1677 }
1678 iCount++;
1679 }
1680 }
1681 }// AODMCparticle
1682 list->Sort();
1683 return iCount;
1684}
1685
006b2a30 1686void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1687
1688 //
1689 // set mom res profiles
1690 //
1691
1692 fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
1693 fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
1694 fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
1695}
1696
1697void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1698 //
1699 // set tracking efficiency histos
1700 //
1701
1702 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1703 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1704 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1705}
1706
1707Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1708
1709 //
1710 // Get smearing on generated momentum
1711 //
1712
1713 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1714
1715 TProfile *fMomRes = 0x0;
1716 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1717 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1718 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1719
1720 if(!fMomRes) {
1721 return 0.;
1722 }
1723
1724
1725 Double_t smear = 0.;
1726
1727 if(pt>20.) {
1728 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1729 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1730 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1731 }
1732 else {
1733
1734 Int_t bin = fMomRes->FindBin(pt);
1735
1736 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1737
1738 }
1739
1740 if(fMomRes) delete fMomRes;
1741
1742 return smear;
1743}
1744
1745void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1746 //
1747 // Fit linear function on momentum resolution at high pT
1748 //
1749
1750 if(!fMomResH1Fit && fMomResH1) {
1751 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1752 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1753 fMomResH1Fit ->SetRange(5.,100.);
1754 }
1755
1756 if(!fMomResH2Fit && fMomResH2) {
1757 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1758 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1759 fMomResH2Fit ->SetRange(5.,100.);
1760 }
1761
1762 if(!fMomResH3Fit && fMomResH3) {
1763 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1764 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1765 fMomResH3Fit ->SetRange(5.,100.);
1766 }
1767
1768}
1769
d89b8229 1770/*
1771Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1772 for(int i = 0; i < particles.GetEntries(); i++){
1773 AliVParticle *vp = (AliVParticle*)particles.At(i);
1774 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1775 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1776 jInp.set_user_index(i);
1777 inputParticles.push_back(jInp);
1778 }
1779
1780 return 0;
1781
1782}
1783*/