Bug fix for case of missing first chamber (Jochen)
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
CommitLineData
dbf053f0 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
21
22
23#include <TROOT.h>
d7ca0e14 24#include <TRandom3.h>
dbf053f0 25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
54424110 28#include <TRefArray.h>
dbf053f0 29#include <TFile.h>
30#include <TKey.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TH3F.h>
34#include <TProfile.h>
35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
38#include "TDatabasePDG.h"
39
40#include "AliAnalysisTaskJetCluster.h"
41#include "AliAnalysisManager.h"
42#include "AliJetFinder.h"
43#include "AliJetHeader.h"
44#include "AliJetReader.h"
dbf053f0 45#include "AliESDEvent.h"
46#include "AliAODEvent.h"
47#include "AliAODHandler.h"
48#include "AliAODTrack.h"
49#include "AliAODJet.h"
50#include "AliAODMCParticle.h"
51#include "AliMCEventHandler.h"
52#include "AliMCEvent.h"
53#include "AliStack.h"
54#include "AliGenPythiaEventHeader.h"
55#include "AliJetKineReaderHeader.h"
56#include "AliGenCocktailEventHeader.h"
57#include "AliInputEventHandler.h"
d6e66a82 58#include "AliAODJetEventBackground.h"
dbf053f0 59
60#include "fastjet/PseudoJet.hh"
61#include "fastjet/ClusterSequenceArea.hh"
62#include "fastjet/AreaDefinition.hh"
63#include "fastjet/JetDefinition.hh"
64// get info on how fastjet was configured
65#include "fastjet/config.h"
66
67
68ClassImp(AliAnalysisTaskJetCluster)
69
54424110 70AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
71 delete fRef;
d7ca0e14 72 delete fRandom;
54424110 73}
74
dbf053f0 75AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
76 fAOD(0x0),
54424110 77 fAODExtension(0x0),
78 fRef(new TRefArray),
dbf053f0 79 fUseAODTrackInput(kFALSE),
80 fUseAODMCInput(kFALSE),
81 fUseGlobalSelection(kFALSE),
9fd9f55b 82 fUseBackgroundCalc(kFALSE),
dbf053f0 83 fFilterMask(0),
84 fTrackTypeRec(kTrackUndef),
bd80a748 85 fTrackTypeGen(kTrackUndef),
86 fNSkipLeadingRan(0),
d7ca0e14 87 fNRandomCones(5),
bd80a748 88 fAvgTrials(1),
dbf053f0 89 fExternalWeight(1),
90 fRecEtaWindow(0.5),
54424110 91 fTrackPtCut(0.),
92 fJetOutputMinPt(1),
0341d0d8 93 fJetTriggerPtCut(0),
54424110 94 fNonStdBranch(""),
0341d0d8 95 fBackgroundBranch(""),
54424110 96 fNonStdFile(""),
dbf053f0 97 fRparam(1.0),
98 fAlgorithm(fastjet::kt_algorithm),
99 fStrategy(fastjet::Best),
100 fRecombScheme(fastjet::BIpt_scheme),
101 fAreaType(fastjet::active_area),
d6e66a82 102 fGhostArea(0.01),
103 fActiveAreaRepeats(1),
104 fGhostEtamax(1.5),
d7ca0e14 105 fRandom(0),
dbf053f0 106 fh1Xsec(0x0),
107 fh1Trials(0x0),
108 fh1PtHard(0x0),
109 fh1PtHardNoW(0x0),
110 fh1PtHardTrials(0x0),
111 fh1NJetsRec(0x0),
112 fh1NConstRec(0x0),
113 fh1NConstLeadingRec(0x0),
114 fh1PtJetsRecIn(0x0),
115 fh1PtJetsLeadingRecIn(0x0),
116 fh1PtJetConstRec(0x0),
117 fh1PtJetConstLeadingRec(0x0),
118 fh1PtTracksRecIn(0x0),
119 fh1PtTracksLeadingRecIn(0x0),
120 fh1NJetsRecRan(0x0),
121 fh1NConstRecRan(0x0),
122 fh1PtJetsLeadingRecInRan(0x0),
123 fh1NConstLeadingRecRan(0x0),
124 fh1PtJetsRecInRan(0x0),
125 fh1PtTracksGenIn(0x0),
a0d08b6d 126 fh1Nch(0x0),
0341d0d8 127 fh1Centrality(0x0),
128 fh1CentralitySelect(0x0),
dbf053f0 129 fh2NRecJetsPt(0x0),
130 fh2NRecTracksPt(0x0),
131 fh2NConstPt(0x0),
132 fh2NConstLeadingPt(0x0),
133 fh2JetPhiEta(0x0),
134 fh2LeadingJetPhiEta(0x0),
135 fh2JetEtaPt(0x0),
136 fh2LeadingJetEtaPt(0x0),
137 fh2TrackEtaPt(0x0),
138 fh2LeadingTrackEtaPt(0x0),
139 fh2JetsLeadingPhiEta(0x0),
140 fh2JetsLeadingPhiPt(0x0),
141 fh2TracksLeadingPhiEta(0x0),
142 fh2TracksLeadingPhiPt(0x0),
143 fh2TracksLeadingJetPhiPt(0x0),
144 fh2JetsLeadingPhiPtW(0x0),
145 fh2TracksLeadingPhiPtW(0x0),
146 fh2TracksLeadingJetPhiPtW(0x0),
147 fh2NRecJetsPtRan(0x0),
148 fh2NConstPtRan(0x0),
149 fh2NConstLeadingPtRan(0x0),
a0d08b6d 150 fh2PtNch(0x0),
151 fh2PtNchRan(0x0),
152 fh2PtNchN(0x0),
153 fh2PtNchNRan(0x0),
dbf053f0 154 fh2TracksLeadingJetPhiPtRan(0x0),
155 fh2TracksLeadingJetPhiPtWRan(0x0),
156 fHistList(0x0)
157{
82a35b1a 158 for(int i = 0;i<3;i++){
159 fh1BiARandomCones[i] = 0;
160 fh1BiARandomConesRan[i] = 0;
161 }
dbf053f0 162}
163
164AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
165 AliAnalysisTaskSE(name),
166 fAOD(0x0),
54424110 167 fAODExtension(0x0),
168 fRef(new TRefArray),
dbf053f0 169 fUseAODTrackInput(kFALSE),
170 fUseAODMCInput(kFALSE),
171 fUseGlobalSelection(kFALSE),
9fd9f55b 172 fUseBackgroundCalc(kFALSE),
dbf053f0 173 fFilterMask(0),
174 fTrackTypeRec(kTrackUndef),
175 fTrackTypeGen(kTrackUndef),
bd80a748 176 fNSkipLeadingRan(0),
d7ca0e14 177 fNRandomCones(5),
dbf053f0 178 fAvgTrials(1),
179 fExternalWeight(1),
180 fRecEtaWindow(0.5),
54424110 181 fTrackPtCut(0.),
182 fJetOutputMinPt(1),
0341d0d8 183 fJetTriggerPtCut(0),
54424110 184 fNonStdBranch(""),
0341d0d8 185 fBackgroundBranch(""),
54424110 186 fNonStdFile(""),
dbf053f0 187 fRparam(1.0),
188 fAlgorithm(fastjet::kt_algorithm),
189 fStrategy(fastjet::Best),
190 fRecombScheme(fastjet::BIpt_scheme),
191 fAreaType(fastjet::active_area),
d6e66a82 192 fGhostArea(0.01),
193 fActiveAreaRepeats(1),
194 fGhostEtamax(1.5),
d7ca0e14 195 fRandom(0),
dbf053f0 196 fh1Xsec(0x0),
197 fh1Trials(0x0),
198 fh1PtHard(0x0),
199 fh1PtHardNoW(0x0),
200 fh1PtHardTrials(0x0),
201 fh1NJetsRec(0x0),
202 fh1NConstRec(0x0),
203 fh1NConstLeadingRec(0x0),
204 fh1PtJetsRecIn(0x0),
205 fh1PtJetsLeadingRecIn(0x0),
206 fh1PtJetConstRec(0x0),
207 fh1PtJetConstLeadingRec(0x0),
208 fh1PtTracksRecIn(0x0),
209 fh1PtTracksLeadingRecIn(0x0),
210 fh1NJetsRecRan(0x0),
211 fh1NConstRecRan(0x0),
212 fh1PtJetsLeadingRecInRan(0x0),
213 fh1NConstLeadingRecRan(0x0),
214 fh1PtJetsRecInRan(0x0),
215 fh1PtTracksGenIn(0x0),
a0d08b6d 216 fh1Nch(0x0),
0341d0d8 217 fh1Centrality(0x0),
218 fh1CentralitySelect(0x0),
dbf053f0 219 fh2NRecJetsPt(0x0),
220 fh2NRecTracksPt(0x0),
221 fh2NConstPt(0x0),
222 fh2NConstLeadingPt(0x0),
223 fh2JetPhiEta(0x0),
224 fh2LeadingJetPhiEta(0x0),
225 fh2JetEtaPt(0x0),
226 fh2LeadingJetEtaPt(0x0),
227 fh2TrackEtaPt(0x0),
228 fh2LeadingTrackEtaPt(0x0),
229 fh2JetsLeadingPhiEta(0x0),
230 fh2JetsLeadingPhiPt(0x0),
231 fh2TracksLeadingPhiEta(0x0),
232 fh2TracksLeadingPhiPt(0x0),
233 fh2TracksLeadingJetPhiPt(0x0),
234 fh2JetsLeadingPhiPtW(0x0),
235 fh2TracksLeadingPhiPtW(0x0),
236 fh2TracksLeadingJetPhiPtW(0x0),
237 fh2NRecJetsPtRan(0x0),
238 fh2NConstPtRan(0x0),
239 fh2NConstLeadingPtRan(0x0),
a0d08b6d 240 fh2PtNch(0x0),
241 fh2PtNchRan(0x0),
242 fh2PtNchN(0x0),
243 fh2PtNchNRan(0x0),
dbf053f0 244 fh2TracksLeadingJetPhiPtRan(0x0),
245 fh2TracksLeadingJetPhiPtWRan(0x0),
246 fHistList(0x0)
247{
82a35b1a 248 for(int i = 0;i<3;i++){
249 fh1BiARandomCones[i] = 0;
250 fh1BiARandomConesRan[i] = 0;
251 }
252 DefineOutput(1, TList::Class());
dbf053f0 253}
254
255
256
257Bool_t AliAnalysisTaskJetCluster::Notify()
258{
259 //
260 // Implemented Notify() to read the cross sections
261 // and number of trials from pyxsec.root
262 //
263 return kTRUE;
264}
265
266void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
267{
268
269 //
270 // Create the output container
271 //
272
d7ca0e14 273 fRandom = new TRandom3(0);
54424110 274
275
dbf053f0 276 // Connect the AOD
277
278
279 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
280
d6e66a82 281
54424110 282
283 if(fNonStdBranch.Length()!=0)
284 {
285 // only create the output branch if we have a name
286 // Create a new branch for jets...
287 // -> cleared in the UserExec....
288 // here we can also have the case that the brnaches are written to a separate file
289
290 TClonesArray *tca = new TClonesArray("AliAODJet", 0);
291 tca->SetName(fNonStdBranch.Data());
292 AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
d6e66a82 293
294
295 TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
296 tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
297 AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
298
299
9fd9f55b 300 if(fUseBackgroundCalc){
301 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
302 AliAODJetEventBackground* evBkg = new AliAODJetEventBackground();
303 evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
304 AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data());
305 }
d7ca0e14 306 // create the branch for the random cones with the same R
307 TString cName = Form("%sRandomCone",fNonStdBranch.Data());
308 if(!AODEvent()->FindListObject(cName.Data())){
309 TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
310 tcaC->SetName(cName.Data());
311 AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
312 }
313
314 // create the branch with the random for the random cones on the random event
315 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
316 if(!AODEvent()->FindListObject(cName.Data())){
317 TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0);
318 tcaCran->SetName(cName.Data());
319 AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data());
320 }
9fd9f55b 321 }
d6e66a82 322
54424110 323 if(fNonStdFile.Length()!=0){
324 //
325 // case that we have an AOD extension we need to fetch the jets from the extended output
326 // we identifay the extension aod event by looking for the branchname
327 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
328 TObjArray* extArray = aodH->GetExtensions();
329 if (extArray) {
330 TIter next(extArray);
331 while ((fAODExtension=(AliAODExtension*)next())){
332 TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
333 if(fDebug>10){
334 Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
335 fAODExtension->GetAOD()->Dump();
336 }
337 if(obj){
338 if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
339 break;
340 }
341 fAODExtension = 0;
342 }
343 }
344 }
345 }
346
347
dbf053f0 348 if(!fHistList)fHistList = new TList();
9c6bd749 349 fHistList->SetOwner();
dbf053f0 350
351 Bool_t oldStatus = TH1::AddDirectoryStatus();
352 TH1::AddDirectory(kFALSE);
353
354 //
355 // Histogram
356
357 const Int_t nBinPt = 200;
358 Double_t binLimitsPt[nBinPt+1];
359 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
360 if(iPt == 0){
361 binLimitsPt[iPt] = 0.0;
362 }
363 else {// 1.0
d7ca0e14 364 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
dbf053f0 365 }
366 }
367
368 const Int_t nBinPhi = 90;
369 Double_t binLimitsPhi[nBinPhi+1];
370 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
371 if(iPhi==0){
372 binLimitsPhi[iPhi] = -1.*TMath::Pi();
373 }
374 else{
375 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
376 }
377 }
378
379
380
381 const Int_t nBinEta = 40;
382 Double_t binLimitsEta[nBinEta+1];
383 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
384 if(iEta==0){
385 binLimitsEta[iEta] = -2.0;
386 }
387 else{
388 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
389 }
390 }
391
d7ca0e14 392 const int nChMax = 4000;
dbf053f0 393
394 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
395 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
396
397 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
398 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
399
400
401 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
402 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
403
404 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
405 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
406 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
407 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
408
409
410 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
411 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
412 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
413
414 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
415 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
416 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
417 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
418 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
419 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
420 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
421 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
422 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
a0d08b6d 423 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
dbf053f0 424
0341d0d8 425 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
426 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
427
dbf053f0 428 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
429 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
430 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
431 //
432
433
434 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
435 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
436 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
437 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
438
a0d08b6d 439 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
440 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);
441 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);
442 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);
443
444
445
dbf053f0 446 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
447 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
448 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
449 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
450
451 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
452 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
453 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
454 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
455
456 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
457 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
458 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
459 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
460
461
462
463 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
464 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
465 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
466 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
467 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
468 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
469 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
470 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
471 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
472 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
473 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
474 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
475
476 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
477 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
478 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
479 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
480
481 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
482 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
483 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
484 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
485
486
d7ca0e14 487 if(fUseBackgroundCalc){
488 for(int i = 0;i<3;i++){
489 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
490 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
491 }
492 }
dbf053f0 493
494 const Int_t saveLevel = 3; // large save level more histos
495 if(saveLevel>0){
496 fHistList->Add(fh1Xsec);
497 fHistList->Add(fh1Trials);
498
499 fHistList->Add(fh1NJetsRec);
500 fHistList->Add(fh1NConstRec);
501 fHistList->Add(fh1NConstLeadingRec);
502 fHistList->Add(fh1PtJetsRecIn);
503 fHistList->Add(fh1PtJetsLeadingRecIn);
504 fHistList->Add(fh1PtTracksRecIn);
505 fHistList->Add(fh1PtTracksLeadingRecIn);
506 fHistList->Add(fh1PtJetConstRec);
507 fHistList->Add(fh1PtJetConstLeadingRec);
508 fHistList->Add(fh1NJetsRecRan);
509 fHistList->Add(fh1NConstRecRan);
510 fHistList->Add(fh1PtJetsLeadingRecInRan);
511 fHistList->Add(fh1NConstLeadingRecRan);
512 fHistList->Add(fh1PtJetsRecInRan);
a0d08b6d 513 fHistList->Add(fh1Nch);
0341d0d8 514 fHistList->Add(fh1Centrality);
515 fHistList->Add(fh1CentralitySelect);
d7ca0e14 516 if(fUseBackgroundCalc){
517 for(int i = 0;i<3;i++){
518 fHistList->Add(fh1BiARandomCones[i]);
519 fHistList->Add(fh1BiARandomConesRan[i]);
520 }
82a35b1a 521 }
dbf053f0 522 fHistList->Add(fh2NRecJetsPt);
523 fHistList->Add(fh2NRecTracksPt);
524 fHistList->Add(fh2NConstPt);
525 fHistList->Add(fh2NConstLeadingPt);
a0d08b6d 526 fHistList->Add(fh2PtNch);
527 fHistList->Add(fh2PtNchRan);
528 fHistList->Add(fh2PtNchN);
529 fHistList->Add(fh2PtNchNRan);
dbf053f0 530 fHistList->Add(fh2JetPhiEta);
531 fHistList->Add(fh2LeadingJetPhiEta);
532 fHistList->Add(fh2JetEtaPt);
533 fHistList->Add(fh2LeadingJetEtaPt);
534 fHistList->Add(fh2TrackEtaPt);
535 fHistList->Add(fh2LeadingTrackEtaPt);
536 fHistList->Add(fh2JetsLeadingPhiEta );
537 fHistList->Add(fh2JetsLeadingPhiPt);
538 fHistList->Add(fh2TracksLeadingPhiEta);
539 fHistList->Add(fh2TracksLeadingPhiPt);
540 fHistList->Add(fh2TracksLeadingJetPhiPt);
541 fHistList->Add(fh2JetsLeadingPhiPtW);
542 fHistList->Add(fh2TracksLeadingPhiPtW);
543 fHistList->Add(fh2TracksLeadingJetPhiPtW);
544 fHistList->Add(fh2NRecJetsPtRan);
545 fHistList->Add(fh2NConstPtRan);
546 fHistList->Add(fh2NConstLeadingPtRan);
547 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
548 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
549 }
550
551 // =========== Switch on Sumw2 for all histos ===========
552 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
553 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
554 if (h1){
555 h1->Sumw2();
556 continue;
557 }
558 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
559 if(hn)hn->Sumw2();
560 }
561 TH1::AddDirectory(oldStatus);
562}
563
564void AliAnalysisTaskJetCluster::Init()
565{
566 //
567 // Initialization
568 //
569
dbf053f0 570 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
571
572}
573
574void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
575{
576
577 if(fUseGlobalSelection){
578 // no selection by the service task, we continue
579 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
580 PostData(1, fHistList);
581 return;
582 }
583
0341d0d8 584
585
54424110 586 // handle and reset the output jet branch
587 // only need this once
588 TClonesArray* jarray = 0;
d6e66a82 589 TClonesArray* jarrayran = 0;
d7ca0e14 590 TClonesArray* rConeArray = 0;
591 TClonesArray* rConeArrayRan = 0;
d6e66a82 592 AliAODJetEventBackground* evBkg = 0;
54424110 593 if(fNonStdBranch.Length()!=0) {
594 if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
595 if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
596 if(jarray)jarray->Delete(); // this is our responsibility, clear before filling again
d6e66a82 597 if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
598 if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random")));
599 if(jarrayran)jarrayran->Delete(); // this is our responsibility, clear before filling again
d7ca0e14 600
9fd9f55b 601 if(fUseBackgroundCalc){
602 if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
603 if(!evBkg) evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())));
604 if(evBkg)evBkg->Reset();
d7ca0e14 605
606 TString cName = Form("%sRandomCone",fNonStdBranch.Data());
607 if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
608 if(!rConeArray)rConeArray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
609 if(rConeArray)rConeArray->Delete();
610
611 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
612 if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
613 if(!rConeArrayRan)rConeArrayRan = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
614 if(rConeArrayRan)rConeArrayRan->Delete();
9fd9f55b 615 }
54424110 616 }
617
0341d0d8 618 static AliAODJetEventBackground* externalBackground = 0;
619 if(!externalBackground&&fBackgroundBranch.Length()){
620 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
621 }
dbf053f0 622 //
623 // Execute analysis for current event
624 //
625 AliESDEvent *fESD = 0;
626 if(fUseAODTrackInput){
627 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
628 if(!fAOD){
629 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
630 return;
631 }
632 // fethc the header
633 }
634 else{
635 // assume that the AOD is in the general output...
636 fAOD = AODEvent();
637 if(!fAOD){
638 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
639 return;
640 }
641 if(fDebug>0){
642 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
643 }
644 }
645
646 Bool_t selectEvent = false;
54424110 647 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
dbf053f0 648
649 if(fAOD){
650 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
54424110 651 TString vtxTitle(vtxAOD->GetTitle());
652
653 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
dbf053f0 654 Float_t zvtx = vtxAOD->GetZ();
655 Float_t yvtx = vtxAOD->GetY();
656 Float_t xvtx = vtxAOD->GetX();
657 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
0341d0d8 658 if(TMath::Abs(zvtx)<20.&&r2<1.){ // apply vertex cut later on
dbf053f0 659 if(physicsSelection)selectEvent = true;
660 }
661 }
662 }
663 if(!selectEvent){
664 PostData(1, fHistList);
665 return;
666 }
667
668 fh1Trials->Fill("#sum{ntrials}",1);
669
670
671 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
672
673 // ==== General variables needed
674
675
676
677 // we simply fetch the tracks/mc particles as a list of AliVParticles
678
679 TList recParticles;
680 TList genParticles;
681
682 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
a0d08b6d 683 Float_t nCh = recParticles.GetEntries();
684 fh1Nch->Fill(nCh);
dbf053f0 685 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
686 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
687 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
688
689 // find the jets....
690
691 vector<fastjet::PseudoJet> inputParticlesRec;
692 vector<fastjet::PseudoJet> inputParticlesRecRan;
82a35b1a 693
694
695 // create a random jet within the acceptance
82a35b1a 696
d7ca0e14 697 if(fUseBackgroundCalc){
698 Double_t etaMax = 0.9 - fRparam;
699 Int_t nCone = 0;
700 Int_t nConeRan = 0;
701 Double_t pT = 1; // small number
702 for(int ir = 0;ir < fNRandomCones;ir++){
703 Double_t eta = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
704 Double_t phi = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
705 // massless jet
706 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
707 Double_t pZ = pT/TMath::Tan(theta);
708 Double_t pX = pT * TMath::Cos(phi);
709 Double_t pY = pT * TMath::Sin(phi);
710 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
711 AliAODJet tmpRec (pX,pY,pZ, p);
712 tmpRec.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
713 if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRec);
714 if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRec);
715 }
716 }
717
718
719 // Generate the random cones
720
721 AliAODJet vTmpRan(1,0,0,1);
dbf053f0 722 for(int i = 0; i < recParticles.GetEntries(); i++){
723 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
724 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
d7ca0e14 725 // we take total momentum here
dbf053f0 726 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
727 jInp.set_user_index(i);
728 inputParticlesRec.push_back(jInp);
729
d7ca0e14 730 if(fUseBackgroundCalc&&rConeArray){
731 for(int ir = 0;ir < fNRandomCones;ir++){
732 AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
733 if(jC&&jC->DeltaR(vp)<fRparam){
734 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
735 }
736 }
82a35b1a 737 }
738
dbf053f0 739 // the randomized input changes eta and phi, but keeps the p_T
bd80a748 740 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
741 Double_t pT = vp->Pt();
d7ca0e14 742 Double_t eta = 1.8 * fRandom->Rndm() - 0.9;
743 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
bd80a748 744
d7ca0e14 745 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
bd80a748 746 Double_t pZ = pT/TMath::Tan(theta);
747
748 Double_t pX = pT * TMath::Cos(phi);
749 Double_t pY = pT * TMath::Sin(phi);
750 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
751 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
d7ca0e14 752
bd80a748 753 jInpRan.set_user_index(i);
754 inputParticlesRecRan.push_back(jInpRan);
d7ca0e14 755 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
82a35b1a 756
d7ca0e14 757 if(fUseBackgroundCalc&&rConeArrayRan){
758 for(int ir = 0;ir < fNRandomCones;ir++){
759 AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
760 if(jC&&jC->DeltaR(&vTmpRan)<fRparam){
761 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRan.Pt(),0);
762 }
763 }
82a35b1a 764 }
bd80a748 765 }
54424110 766
767 // fill the tref array, only needed when we write out jets
768 if(jarray){
769 if(i == 0){
770 fRef->Delete(); // make sure to delete before placement new...
771 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
772 }
773 fRef->Add(vp);
774 }
d7ca0e14 775 }// recparticles
776 if(fUseBackgroundCalc){
0341d0d8 777 Float_t jetArea = fRparam*fRparam*TMath::Pi();
d7ca0e14 778 for(int ir = 0;ir < fNRandomCones;ir++){
779 // rescale the momntum vectors for the random cones
780 if(!rConeArray)continue;
781 AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
782 if(rC){
783 Double_t eta = rC->Eta();
784 Double_t phi = rC->Phi();
785 // massless jet, unit vector
786 Double_t pT = rC->ChargedBgEnergy();
787 if(pT<=0)pT = 0.1; // for almost empty events
788 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
789 Double_t pZ = pT/TMath::Tan(theta);
790 Double_t pX = pT * TMath::Cos(phi);
791 Double_t pY = pT * TMath::Sin(phi);
792 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
793 rC->SetPxPyPzE(pX,pY,pZ, p);
794 rC->SetBgEnergy(0,0);
0341d0d8 795 rC->SetEffArea(jetArea,0);
d7ca0e14 796 }
797 rC = (AliAODJet*)rConeArrayRan->At(ir);
798 // same wit random
799 if(rC){
800 Double_t eta = rC->Eta();
801 Double_t phi = rC->Phi();
802 // massless jet, unit vector
803 Double_t pT = rC->ChargedBgEnergy();
804 if(pT<=0)pT = 0.1;// for almost empty events
805 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
806 Double_t pZ = pT/TMath::Tan(theta);
807 Double_t pX = pT * TMath::Cos(phi);
808 Double_t pY = pT * TMath::Sin(phi);
809 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
810 rC->SetPxPyPzE(pX,pY,pZ, p);
811 rC->SetBgEnergy(0,0);
0341d0d8 812 rC->SetEffArea(jetArea,0);
d7ca0e14 813 }
814 }
dbf053f0 815 }
816
d7ca0e14 817
dbf053f0 818 if(inputParticlesRec.size()==0){
819 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
820 PostData(1, fHistList);
821 return;
822 }
823
824 // run fast jet
4dd6159d 825 // employ setters for these...
826
d6e66a82 827
9fd9f55b 828 // now create the object that holds info about ghosts
829 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
830 // reduce CPU time...
831 fGhostArea = 0.5;
832 fActiveAreaRepeats = 0;
833 }
d7ca0e14 834
835 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
4dd6159d 836 fastjet::AreaType areaType = fastjet::active_area;
837 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
82a35b1a 838 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
dbf053f0 839 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
4dd6159d 840 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
dbf053f0 841
d6e66a82 842 //range where to compute background
843 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
844 phiMin = 0;
845 phiMax = 2*TMath::Pi();
846 rapMax = fGhostEtamax - fRparam;
847 rapMin = - fGhostEtamax + fRparam;
848 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
849
850
851
dbf053f0 852 inclusiveJets = clustSeq.inclusive_jets();
853 sortedJets = sorted_by_pt(inclusiveJets);
854
dbf053f0 855 fh1NJetsRec->Fill(sortedJets.size());
856
857 // loop over all jets an fill information, first on is the leading jet
858
54424110 859 Int_t nRecOver = inclusiveJets.size();
860 Int_t nRec = inclusiveJets.size();
861 if(inclusiveJets.size()>0){
862 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
0341d0d8 863 Double_t area=clustSeq.area(sortedJets[0]);
864 leadingJet.SetEffArea(area,0);
dbf053f0 865 Float_t pt = leadingJet.Pt();
54424110 866 Int_t nAodOutJets = 0;
867 Int_t nAodOutTracks = 0;
868 AliAODJet *aodOutJet = 0;
dbf053f0 869
870 Int_t iCount = 0;
871 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
872 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
873 while(pt<ptCut&&iCount<nRec){
874 nRecOver--;
875 iCount++;
876 if(iCount<nRec){
877 pt = sortedJets[iCount].perp();
878 }
879 }
880 if(nRecOver<=0)break;
881 fh2NRecJetsPt->Fill(ptCut,nRecOver);
882 }
883 Float_t phi = leadingJet.Phi();
884 if(phi<0)phi+=TMath::Pi()*2.;
885 Float_t eta = leadingJet.Eta();
0341d0d8 886 Float_t pTback = 0;
887 if(externalBackground){
888 // carefull has to be filled in a task before
889 // todo, ReArrange to the botom
890 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
891 }
892 pt = leadingJet.Pt() - pTback;
dbf053f0 893 // correlation of leading jet with tracks
894 TIterator *recIter = recParticles.MakeIterator();
dbf053f0 895 recIter->Reset();
225f0094 896 AliVParticle *tmpRecTrack = 0;
dbf053f0 897 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
898 Float_t tmpPt = tmpRecTrack->Pt();
899 // correlation
900 Float_t tmpPhi = tmpRecTrack->Phi();
901 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
902 Float_t dPhi = phi - tmpPhi;
903 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
904 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 905 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
906 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
54424110 907 }
908
909
910
911 for(int j = 0; j < nRec;j++){
912 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
913 aodOutJet = 0;
914 nAodOutTracks = 0;
915 Float_t tmpPt = tmpRec.Pt();
916 fh1PtJetsRecIn->Fill(tmpPt);
917 // Fill Spectra with constituents
918 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
919 fh1NConstRec->Fill(constituents.size());
920 fh2PtNch->Fill(nCh,tmpPt);
921 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
922 fh2NConstPt->Fill(tmpPt,constituents.size());
923 // loop over constiutents and fill spectrum
924
925 // Add the jet information and the track references to the output AOD
926
927 if(tmpPt>fJetOutputMinPt&&jarray){
928 aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
42bfdee2 929 Double_t area=clustSeq.area(sortedJets[j]);
930 aodOutJet->SetEffArea(area,0);
54424110 931 }
932
933
934 for(unsigned int ic = 0; ic < constituents.size();ic++){
935 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
936 fh1PtJetConstRec->Fill(part->Pt());
937 if(aodOutJet){
938 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
939 }
940 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
941 }
54424110 942
943 // correlation
944 Float_t tmpPhi = tmpRec.Phi();
945 Float_t tmpEta = tmpRec.Eta();
946 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
947
0341d0d8 948 Float_t tmpPtBack = 0;
949 if(externalBackground){
950 // carefull has to be filled in a task before
951 // todo, ReArrange to the botom
952 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
953 }
954 tmpPt = tmpPt - tmpPtBack;
955 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
956
54424110 957 if(j==0){
958 fh1PtJetsLeadingRecIn->Fill(tmpPt);
959 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
960 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
961 fh1NConstLeadingRec->Fill(constituents.size());
962 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
963 continue;
964 }
965 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
966 fh2JetEtaPt->Fill(tmpEta,tmpPt);
967 Float_t dPhi = phi - tmpPhi;
968 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
969 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
970 Float_t dEta = eta - tmpRec.Eta();
971 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
972 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
973 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
974 }
975 delete recIter;
d7ca0e14 976
977 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
d6e66a82 978
979 if(evBkg){
980 vector<fastjet::PseudoJet> jets2=sortedJets;
981 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
982 Double_t bkg1=0;
983 Double_t sigma1=0.;
984 Double_t meanarea1=0.;
985 Double_t bkg2=0;
986 Double_t sigma2=0.;
987 Double_t meanarea2=0.;
82a35b1a 988
4a33f620 989 clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
d6e66a82 990 evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
f974a156 991
992 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
993 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
82a35b1a 994
4a33f620 995 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
d6e66a82 996 evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
f974a156 997 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
998 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
999
d6e66a82 1000 }
54424110 1001 }
d7ca0e14 1002
82a35b1a 1003
1004
1005
dbf053f0 1006
82a35b1a 1007 // fill track information
1008 Int_t nTrackOver = recParticles.GetSize();
dbf053f0 1009 // do the same for tracks and jets
82a35b1a 1010
1011 if(nTrackOver>0){
dbf053f0 1012 TIterator *recIter = recParticles.MakeIterator();
1013 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1014 Float_t pt = tmpRec->Pt();
82a35b1a 1015
dbf053f0 1016 // Printf("Leading track p_t %3.3E",pt);
1017 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1018 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1019 while(pt<ptCut&&tmpRec){
1020 nTrackOver--;
1021 tmpRec = (AliVParticle*)(recIter->Next());
1022 if(tmpRec){
1023 pt = tmpRec->Pt();
1024 }
1025 }
1026 if(nTrackOver<=0)break;
1027 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1028 }
1029
1030 recIter->Reset();
1031 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1032 Float_t phi = leading->Phi();
1033 if(phi<0)phi+=TMath::Pi()*2.;
1034 Float_t eta = leading->Eta();
1035 pt = leading->Pt();
1036 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1037 Float_t tmpPt = tmpRec->Pt();
1038 Float_t tmpEta = tmpRec->Eta();
1039 fh1PtTracksRecIn->Fill(tmpPt);
1040 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1041 if(tmpRec==leading){
1042 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1043 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1044 continue;
1045 }
1046 // correlation
1047 Float_t tmpPhi = tmpRec->Phi();
1048
1049 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1050 Float_t dPhi = phi - tmpPhi;
1051 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1052 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1053 Float_t dEta = eta - tmpRec->Eta();
1054 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1055 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1056 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1057 }
1058 delete recIter;
1059 }
1060
1061 // find the random jets
1062 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
d6e66a82 1063 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
dbf053f0 1064
1065 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1066 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1067
1068 // fill the jet information from random track
1069
1070 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1071
1072 // loop over all jets an fill information, first on is the leading jet
1073
1074 Int_t nRecOverRan = inclusiveJetsRan.size();
1075 Int_t nRecRan = inclusiveJetsRan.size();
1076 if(inclusiveJetsRan.size()>0){
1077 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1078 Float_t pt = leadingJet.Pt();
1079
1080 Int_t iCount = 0;
1081 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1082 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1083 while(pt<ptCut&&iCount<nRecRan){
1084 nRecOverRan--;
1085 iCount++;
1086 if(iCount<nRecRan){
1087 pt = sortedJetsRan[iCount].perp();
1088 }
1089 }
1090 if(nRecOverRan<=0)break;
1091 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1092 }
1093 Float_t phi = leadingJet.Phi();
1094 if(phi<0)phi+=TMath::Pi()*2.;
1095 pt = leadingJet.Pt();
1096
1097 // correlation of leading jet with random tracks
1098
1099 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1100 {
1101 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1102 // correlation
1103 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1104 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1105 Float_t dPhi = phi - tmpPhi;
1106 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1107 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
dbf053f0 1108 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1109 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1110 }
1111
d6e66a82 1112 Int_t nAodOutJetsRan = 0;
1113 AliAODJet *aodOutJetRan = 0;
dbf053f0 1114 for(int j = 0; j < nRecRan;j++){
1115 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1116 Float_t tmpPt = tmpRec.Pt();
1117 fh1PtJetsRecInRan->Fill(tmpPt);
1118 // Fill Spectra with constituents
1119 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1120 fh1NConstRecRan->Fill(constituents.size());
1121 fh2NConstPtRan->Fill(tmpPt,constituents.size());
a0d08b6d 1122 fh2PtNchRan->Fill(nCh,tmpPt);
1123 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
d6e66a82 1124
1125
1126 if(tmpPt>fJetOutputMinPt&&jarrayran){
1127 aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1128 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1129
1130 aodOutJetRan->SetEffArea(arearan,0); }
1131
1132
1133
1134
1135 for(unsigned int ic = 0; ic < constituents.size();ic++){
1136 // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1137 // fh1PtJetConstRec->Fill(part->Pt());
1138 if(aodOutJetRan){
1139 aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
d7ca0e14 1140 }
1141 }
d6e66a82 1142
1143
1144
1145
dbf053f0 1146 // correlation
1147 Float_t tmpPhi = tmpRec.Phi();
1148 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1149
1150 if(j==0){
1151 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1152 fh1NConstLeadingRecRan->Fill(constituents.size());
1153 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1154 continue;
1155 }
1156 }
d6e66a82 1157
1158
1159 if(evBkg){
1160 Double_t bkg3=0.;
1161 Double_t sigma3=0.;
1162 Double_t meanarea3=0.;
4a33f620 1163 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
d6e66a82 1164 evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
d7ca0e14 1165 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1166 /*
82a35b1a 1167 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
9c6bd749 1168 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
d7ca0e14 1169 */
d6e66a82 1170 }
1171
1172
1173
dbf053f0 1174 }
0341d0d8 1175
1176
1177 // do the event selection if activated
1178 if(fJetTriggerPtCut>0){
1179 bool select = false;
1180 Float_t cent = 100;
1181 if(fAOD)cent = fAOD->GetHeader()->GetCentrality();
1182 Float_t minPt = 9999999;
1183 // hard coded for now ...
1184 // 54.50 44.5 29.5 18.5 for anti-kt rejection 10E-3
1185 if(cent<10)minPt = 50;
1186 else if(cent<20)minPt = 42;
1187 else if(cent<50)minPt = 28;
1188 else if(cent<80)minPt = 18;
1189 fh1Centrality->Fill(cent);
1190
1191 float rho = 0;
1192 if(externalBackground)rho = externalBackground->GetBackground(2);
1193 if(jarray&&cent<80){
1194 for(int i = 0;i < jarray->GetEntriesFast();i++){
1195 AliAODJet *jet = (AliAODJet*)jarray->At(i);
1196 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1197 if(ptSub>=minPt){
1198 select = true;
1199 break;
1200 }
1201 }
1202 }
1203
1204 if(select){
1205 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1206 fh1CentralitySelect->Fill(cent);
1207 aodH->SetFillAOD(kTRUE);
1208 }
1209 }
dbf053f0 1210
1211 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1212 PostData(1, fHistList);
1213}
1214
1215void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1216{
1217// Terminate analysis
1218//
1219 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1220}
1221
1222
1223Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1224
1225 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1226
1227 Int_t iCount = 0;
1228 if(type==kTrackAOD){
1229 AliAODEvent *aod = 0;
1230 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1231 else aod = AODEvent();
1232 if(!aod){
1233 return iCount;
1234 }
1235 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1236 AliAODTrack *tr = aod->GetTrack(it);
1237 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1238 if(TMath::Abs(tr->Eta())>0.9)continue;
0341d0d8 1239 if(tr->Pt()<fTrackPtCut)continue;
dbf053f0 1240 list->Add(tr);
1241 iCount++;
1242 }
1243 }
1244 else if (type == kTrackKineAll||type == kTrackKineCharged){
1245 AliMCEvent* mcEvent = MCEvent();
1246 if(!mcEvent)return iCount;
1247 // we want to have alivpartilces so use get track
1248 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1249 if(!mcEvent->IsPhysicalPrimary(it))continue;
1250 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1251 if(type == kTrackKineAll){
0341d0d8 1252 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1253 list->Add(part);
1254 iCount++;
1255 }
1256 else if(type == kTrackKineCharged){
1257 if(part->Particle()->GetPDG()->Charge()==0)continue;
0341d0d8 1258 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1259 list->Add(part);
1260 iCount++;
1261 }
1262 }
1263 }
1264 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1265 AliAODEvent *aod = 0;
959508f4 1266 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
dbf053f0 1267 else aod = AODEvent();
1268 if(!aod)return iCount;
1269 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1270 if(!tca)return iCount;
1271 for(int it = 0;it < tca->GetEntriesFast();++it){
225f0094 1272 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
dbf053f0 1273 if(!part->IsPhysicalPrimary())continue;
1274 if(type == kTrackAODMCAll){
0341d0d8 1275 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1276 list->Add(part);
1277 iCount++;
1278 }
1279 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1280 if(part->Charge()==0)continue;
0341d0d8 1281 if(part->Pt()<fTrackPtCut)continue;
dbf053f0 1282 if(kTrackAODMCCharged){
1283 list->Add(part);
1284 }
1285 else {
1286 if(TMath::Abs(part->Eta())>0.9)continue;
1287 list->Add(part);
1288 }
1289 iCount++;
1290 }
1291 }
1292 }// AODMCparticle
1293 list->Sort();
1294 return iCount;
dbf053f0 1295}
1296
1297/*
1298Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1299 for(int i = 0; i < particles.GetEntries(); i++){
1300 AliVParticle *vp = (AliVParticle*)particles.At(i);
1301 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1302 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1303 jInp.set_user_index(i);
1304 inputParticles.push_back(jInp);
1305 }
1306
1307 return 0;
1308
1309}
1310*/