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