1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets
4 // *******************************************
7 /**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
26 #include <TInterpreter.h>
28 #include <TRefArray.h>
36 #include <TLorentzVector.h>
37 #include <TClonesArray.h>
38 #include "TDatabasePDG.h"
40 #include "AliAnalysisTaskJetCluster.h"
41 #include "AliAnalysisManager.h"
42 #include "AliJetFinder.h"
43 #include "AliJetHeader.h"
44 #include "AliJetReader.h"
45 #include "AliESDEvent.h"
46 #include "AliAODEvent.h"
47 #include "AliAODHandler.h"
48 #include "AliAODExtension.h"
49 #include "AliAODTrack.h"
50 #include "AliAODJet.h"
51 #include "AliAODMCParticle.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliJetKineReaderHeader.h"
57 #include "AliGenCocktailEventHeader.h"
58 #include "AliInputEventHandler.h"
59 #include "AliAODJetEventBackground.h"
61 #include "fastjet/PseudoJet.hh"
62 #include "fastjet/ClusterSequenceArea.hh"
63 #include "fastjet/AreaDefinition.hh"
64 #include "fastjet/JetDefinition.hh"
65 // get info on how fastjet was configured
66 #include "fastjet/config.h"
69 ClassImp(AliAnalysisTaskJetCluster)
71 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
78 if(fTCAJetsOut)fTCAJetsOut->Delete();
80 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
81 delete fTCAJetsOutRan;
82 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
83 delete fTCARandomConesOut;
84 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
85 delete fTCARandomConesOutRan;
86 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
87 delete fAODJetBackgroundOut;
90 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
95 fUseAODTrackInput(kFALSE),
96 fUseAODMCInput(kFALSE),
97 fUseBackgroundCalc(kFALSE),
98 fEventSelection(kFALSE),
102 fTrackTypeRec(kTrackUndef),
103 fTrackTypeGen(kTrackUndef),
105 fNSkipLeadingCone(0),
109 fTrackEtaWindow(0.9),
112 fJetOutputMinPt(0.150),
113 fMaxTrackPtInJet(100.),
120 fBackgroundBranch(""),
123 fAlgorithm(fastjet::kt_algorithm),
124 fStrategy(fastjet::Best),
125 fRecombScheme(fastjet::BIpt_scheme),
126 fAreaType(fastjet::active_area),
128 fActiveAreaRepeats(1),
132 fTCARandomConesOut(0x0),
133 fTCARandomConesOutRan(0x0),
134 fAODJetBackgroundOut(0x0),
140 fh1PtHardTrials(0x0),
143 fh1NConstLeadingRec(0x0),
145 fh1PtJetsLeadingRecIn(0x0),
146 fh1PtJetConstRec(0x0),
147 fh1PtJetConstLeadingRec(0x0),
148 fh1PtTracksRecIn(0x0),
149 fh1PtTracksLeadingRecIn(0x0),
151 fh1NConstRecRan(0x0),
152 fh1PtJetsLeadingRecInRan(0x0),
153 fh1NConstLeadingRecRan(0x0),
154 fh1PtJetsRecInRan(0x0),
155 fh1PtTracksGenIn(0x0),
157 fh1CentralityPhySel(0x0),
159 fh1CentralitySelect(0x0),
164 fh2NRecTracksPt(0x0),
166 fh2NConstLeadingPt(0x0),
168 fh2LeadingJetPhiEta(0x0),
170 fh2LeadingJetEtaPt(0x0),
172 fh2LeadingTrackEtaPt(0x0),
173 fh2JetsLeadingPhiEta(0x0),
174 fh2JetsLeadingPhiPt(0x0),
175 fh2TracksLeadingPhiEta(0x0),
176 fh2TracksLeadingPhiPt(0x0),
177 fh2TracksLeadingJetPhiPt(0x0),
178 fh2JetsLeadingPhiPtW(0x0),
179 fh2TracksLeadingPhiPtW(0x0),
180 fh2TracksLeadingJetPhiPtW(0x0),
181 fh2NRecJetsPtRan(0x0),
183 fh2NConstLeadingPtRan(0x0),
188 fh2TracksLeadingJetPhiPtRan(0x0),
189 fh2TracksLeadingJetPhiPtWRan(0x0),
196 for(int i = 0;i<3;i++){
197 fh1BiARandomCones[i] = 0;
198 fh1BiARandomConesRan[i] = 0;
200 for(int i = 0;i<kMaxCent;i++){
201 fh2JetsLeadingPhiPtC[i] = 0;
202 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
203 fh2TracksLeadingJetPhiPtC[i] = 0;
204 fh2TracksLeadingJetPhiPtWC[i] = 0;
208 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
209 AliAnalysisTaskSE(name),
213 fUseAODTrackInput(kFALSE),
214 fUseAODMCInput(kFALSE),
215 fUseBackgroundCalc(kFALSE),
216 fEventSelection(kFALSE),
220 fTrackTypeRec(kTrackUndef),
221 fTrackTypeGen(kTrackUndef),
223 fNSkipLeadingCone(0),
227 fTrackEtaWindow(0.9),
230 fJetOutputMinPt(0.150),
231 fMaxTrackPtInJet(100.),
238 fBackgroundBranch(""),
241 fAlgorithm(fastjet::kt_algorithm),
242 fStrategy(fastjet::Best),
243 fRecombScheme(fastjet::BIpt_scheme),
244 fAreaType(fastjet::active_area),
246 fActiveAreaRepeats(1),
250 fTCARandomConesOut(0x0),
251 fTCARandomConesOutRan(0x0),
252 fAODJetBackgroundOut(0x0),
258 fh1PtHardTrials(0x0),
261 fh1NConstLeadingRec(0x0),
263 fh1PtJetsLeadingRecIn(0x0),
264 fh1PtJetConstRec(0x0),
265 fh1PtJetConstLeadingRec(0x0),
266 fh1PtTracksRecIn(0x0),
267 fh1PtTracksLeadingRecIn(0x0),
269 fh1NConstRecRan(0x0),
270 fh1PtJetsLeadingRecInRan(0x0),
271 fh1NConstLeadingRecRan(0x0),
272 fh1PtJetsRecInRan(0x0),
273 fh1PtTracksGenIn(0x0),
275 fh1CentralityPhySel(0x0),
277 fh1CentralitySelect(0x0),
282 fh2NRecTracksPt(0x0),
284 fh2NConstLeadingPt(0x0),
286 fh2LeadingJetPhiEta(0x0),
288 fh2LeadingJetEtaPt(0x0),
290 fh2LeadingTrackEtaPt(0x0),
291 fh2JetsLeadingPhiEta(0x0),
292 fh2JetsLeadingPhiPt(0x0),
293 fh2TracksLeadingPhiEta(0x0),
294 fh2TracksLeadingPhiPt(0x0),
295 fh2TracksLeadingJetPhiPt(0x0),
296 fh2JetsLeadingPhiPtW(0x0),
297 fh2TracksLeadingPhiPtW(0x0),
298 fh2TracksLeadingJetPhiPtW(0x0),
299 fh2NRecJetsPtRan(0x0),
301 fh2NConstLeadingPtRan(0x0),
306 fh2TracksLeadingJetPhiPtRan(0x0),
307 fh2TracksLeadingJetPhiPtWRan(0x0),
313 for(int i = 0;i<3;i++){
314 fh1BiARandomCones[i] = 0;
315 fh1BiARandomConesRan[i] = 0;
317 for(int i = 0;i<kMaxCent;i++){
318 fh2JetsLeadingPhiPtC[i] = 0;
319 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
320 fh2TracksLeadingJetPhiPtC[i] = 0;
321 fh2TracksLeadingJetPhiPtWC[i] = 0;
323 DefineOutput(1, TList::Class());
328 Bool_t AliAnalysisTaskJetCluster::Notify()
331 // Implemented Notify() to read the cross sections
332 // and number of trials from pyxsec.root
337 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
341 // Create the output container
344 fRandom = new TRandom3(0);
350 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
354 if(fNonStdBranch.Length()!=0)
356 // only create the output branch if we have a name
357 // Create a new branch for jets...
358 // -> cleared in the UserExec....
359 // here we can also have the case that the brnaches are written to a separate file
362 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
363 fTCAJetsOut->SetName(fNonStdBranch.Data());
364 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
367 if(fJetTypes&kJetRan){
368 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
369 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
370 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
373 if(fUseBackgroundCalc){
374 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
375 fAODJetBackgroundOut = new AliAODJetEventBackground();
376 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
377 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
380 // create the branch for the random cones with the same R
381 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
385 if(!AODEvent()->FindListObject(cName.Data())){
386 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
387 fTCARandomConesOut->SetName(cName.Data());
388 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
391 // create the branch with the random for the random cones on the random event
392 if(fJetTypes&kRCRan){
393 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
394 if(!AODEvent()->FindListObject(cName.Data())){
395 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
396 fTCARandomConesOutRan->SetName(cName.Data());
397 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
402 if(fNonStdFile.Length()!=0){
404 // case that we have an AOD extension we need to fetch the jets from the extended output
405 // we identify the extension aod event by looking for the branchname
406 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
407 // case that we have an AOD extension we need can fetch the background maybe from the extended output
408 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
413 if(!fHistList)fHistList = new TList();
414 fHistList->SetOwner();
415 PostData(1, fHistList); // post data in any case once
417 Bool_t oldStatus = TH1::AddDirectoryStatus();
418 TH1::AddDirectory(kFALSE);
423 const Int_t nBinPt = 100;
424 Double_t binLimitsPt[nBinPt+1];
425 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
427 binLimitsPt[iPt] = 0.0;
430 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
434 const Int_t nBinPhi = 90;
435 Double_t binLimitsPhi[nBinPhi+1];
436 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
438 binLimitsPhi[iPhi] = -1.*TMath::Pi();
441 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
447 const Int_t nBinEta = 40;
448 Double_t binLimitsEta[nBinEta+1];
449 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
451 binLimitsEta[iEta] = -2.0;
454 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
458 const int nChMax = 5000;
460 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
461 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
463 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
464 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
467 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
468 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
470 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
471 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
472 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
473 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
476 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
477 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
478 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
480 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
481 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
482 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
483 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
484 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
485 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
486 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
487 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
488 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
489 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
491 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
492 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
493 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
495 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
496 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
497 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
499 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
500 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
501 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
505 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
506 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
507 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
508 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
510 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
511 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);
512 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);
513 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);
517 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
518 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
519 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
520 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
522 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
523 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
524 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
525 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
527 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
528 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
529 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
530 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
534 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
535 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
536 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
537 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
538 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
539 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
540 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
541 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
542 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
543 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
544 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
545 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
547 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
548 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
549 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
550 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
552 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
553 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
554 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
555 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
558 if(fNRandomCones>0&&fUseBackgroundCalc){
559 for(int i = 0;i<3;i++){
560 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
561 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
565 for(int i = 0;i < kMaxCent;i++){
566 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
567 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
568 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
569 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
572 const Int_t saveLevel = 3; // large save level more histos
574 fHistList->Add(fh1Xsec);
575 fHistList->Add(fh1Trials);
577 fHistList->Add(fh1NJetsRec);
578 fHistList->Add(fh1NConstRec);
579 fHistList->Add(fh1NConstLeadingRec);
580 fHistList->Add(fh1PtJetsRecIn);
581 fHistList->Add(fh1PtJetsLeadingRecIn);
582 fHistList->Add(fh1PtTracksRecIn);
583 fHistList->Add(fh1PtTracksLeadingRecIn);
584 fHistList->Add(fh1PtJetConstRec);
585 fHistList->Add(fh1PtJetConstLeadingRec);
586 fHistList->Add(fh1NJetsRecRan);
587 fHistList->Add(fh1NConstRecRan);
588 fHistList->Add(fh1PtJetsLeadingRecInRan);
589 fHistList->Add(fh1NConstLeadingRecRan);
590 fHistList->Add(fh1PtJetsRecInRan);
591 fHistList->Add(fh1Nch);
592 fHistList->Add(fh1Centrality);
593 fHistList->Add(fh1CentralitySelect);
594 fHistList->Add(fh1CentralityPhySel);
595 fHistList->Add(fh1Z);
596 fHistList->Add(fh1ZSelect);
597 fHistList->Add(fh1ZPhySel);
598 if(fNRandomCones>0&&fUseBackgroundCalc){
599 for(int i = 0;i<3;i++){
600 fHistList->Add(fh1BiARandomCones[i]);
601 fHistList->Add(fh1BiARandomConesRan[i]);
604 for(int i = 0;i < kMaxCent;i++){
605 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
606 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
607 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
608 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
611 fHistList->Add(fh2NRecJetsPt);
612 fHistList->Add(fh2NRecTracksPt);
613 fHistList->Add(fh2NConstPt);
614 fHistList->Add(fh2NConstLeadingPt);
615 fHistList->Add(fh2PtNch);
616 fHistList->Add(fh2PtNchRan);
617 fHistList->Add(fh2PtNchN);
618 fHistList->Add(fh2PtNchNRan);
619 fHistList->Add(fh2JetPhiEta);
620 fHistList->Add(fh2LeadingJetPhiEta);
621 fHistList->Add(fh2JetEtaPt);
622 fHistList->Add(fh2LeadingJetEtaPt);
623 fHistList->Add(fh2TrackEtaPt);
624 fHistList->Add(fh2LeadingTrackEtaPt);
625 fHistList->Add(fh2JetsLeadingPhiEta );
626 fHistList->Add(fh2JetsLeadingPhiPt);
627 fHistList->Add(fh2TracksLeadingPhiEta);
628 fHistList->Add(fh2TracksLeadingPhiPt);
629 fHistList->Add(fh2TracksLeadingJetPhiPt);
630 fHistList->Add(fh2JetsLeadingPhiPtW);
631 fHistList->Add(fh2TracksLeadingPhiPtW);
632 fHistList->Add(fh2TracksLeadingJetPhiPtW);
633 fHistList->Add(fh2NRecJetsPtRan);
634 fHistList->Add(fh2NConstPtRan);
635 fHistList->Add(fh2NConstLeadingPtRan);
636 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
637 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
640 // =========== Switch on Sumw2 for all histos ===========
641 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
642 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
647 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
650 TH1::AddDirectory(oldStatus);
653 void AliAnalysisTaskJetCluster::Init()
659 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
663 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
666 // handle and reset the output jet branch
668 if(fTCAJetsOut)fTCAJetsOut->Delete();
669 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
670 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
671 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
672 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
674 AliAODJetEventBackground* externalBackground = 0;
675 if(!externalBackground&&fBackgroundBranch.Length()){
676 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
677 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
678 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
681 // Execute analysis for current event
683 AliESDEvent *fESD = 0;
684 if(fUseAODTrackInput){
685 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
687 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
693 // assume that the AOD is in the general output...
696 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
700 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
704 Bool_t selectEvent = false;
705 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
711 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
712 TString vtxTitle(vtxAOD->GetTitle());
713 zVtx = vtxAOD->GetZ();
715 cent = fAOD->GetHeader()->GetCentrality();
716 if(cent<10)cenClass = 0;
717 else if(cent<30)cenClass = 1;
718 else if(cent<50)cenClass = 2;
719 else if(cent<80)cenClass = 3;
720 if(physicsSelection){
721 fh1CentralityPhySel->Fill(cent);
722 fh1ZPhySel->Fill(zVtx);
726 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
727 Float_t yvtx = vtxAOD->GetY();
728 Float_t xvtx = vtxAOD->GetX();
729 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
730 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
731 if(physicsSelection){
737 if(cent<fCentCutLo||cent>fCentCutUp){
748 PostData(1, fHistList);
751 fh1Centrality->Fill(cent);
753 fh1Trials->Fill("#sum{ntrials}",1);
756 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
758 // ==== General variables needed
762 // we simply fetch the tracks/mc particles as a list of AliVParticles
767 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
768 Float_t nCh = recParticles.GetEntries();
770 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
771 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
772 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
776 vector<fastjet::PseudoJet> inputParticlesRec;
777 vector<fastjet::PseudoJet> inputParticlesRecRan;
779 // Generate the random cones
781 AliAODJet vTmpRan(1,0,0,1);
782 for(int i = 0; i < recParticles.GetEntries(); i++){
783 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
784 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
785 // we take total momentum here
786 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
787 jInp.set_user_index(i);
788 inputParticlesRec.push_back(jInp);
790 // the randomized input changes eta and phi, but keeps the p_T
791 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
792 Double_t pT = vp->Pt();
793 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
794 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
796 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
797 Double_t pZ = pT/TMath::Tan(theta);
799 Double_t pX = pT * TMath::Cos(phi);
800 Double_t pY = pT * TMath::Sin(phi);
801 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
802 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
804 jInpRan.set_user_index(i);
805 inputParticlesRecRan.push_back(jInpRan);
806 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
809 // fill the tref array, only needed when we write out jets
812 fRef->Delete(); // make sure to delete before placement new...
813 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
819 if(inputParticlesRec.size()==0){
820 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
821 PostData(1, fHistList);
826 // employ setters for these...
829 // now create the object that holds info about ghosts
831 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
832 // reduce CPU time...
834 fActiveAreaRepeats = 0;
838 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
839 fastjet::AreaType areaType = fastjet::active_area;
840 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
841 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
842 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
844 //range where to compute background
845 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
847 phiMax = 2*TMath::Pi();
848 rapMax = fGhostEtamax - fRparam;
849 rapMin = - fGhostEtamax + fRparam;
850 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
853 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
854 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
857 fh1NJetsRec->Fill(sortedJets.size());
859 // loop over all jets an fill information, first on is the leading jet
861 Int_t nRecOver = inclusiveJets.size();
862 Int_t nRec = inclusiveJets.size();
863 if(inclusiveJets.size()>0){
864 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
865 Double_t area = clustSeq.area(sortedJets[0]);
866 leadingJet.SetEffArea(area,0);
867 Float_t pt = leadingJet.Pt();
868 Int_t nAodOutJets = 0;
869 Int_t nAodOutTracks = 0;
870 AliAODJet *aodOutJet = 0;
873 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
874 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
875 while(pt<ptCut&&iCount<nRec){
879 pt = sortedJets[iCount].perp();
882 if(nRecOver<=0)break;
883 fh2NRecJetsPt->Fill(ptCut,nRecOver);
885 Float_t phi = leadingJet.Phi();
886 if(phi<0)phi+=TMath::Pi()*2.;
887 Float_t eta = leadingJet.Eta();
889 if(externalBackground){
890 // carefull has to be filled in a task before
891 // todo, ReArrange to the botom
892 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
894 pt = leadingJet.Pt() - pTback;
895 // correlation of leading jet with tracks
896 TIterator *recIter = recParticles.MakeIterator();
898 AliVParticle *tmpRecTrack = 0;
899 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
900 Float_t tmpPt = tmpRecTrack->Pt();
902 Float_t tmpPhi = tmpRecTrack->Phi();
903 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
904 Float_t dPhi = phi - tmpPhi;
905 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
906 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
907 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
908 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
910 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
911 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
917 TLorentzVector vecareab;
918 for(int j = 0; j < nRec;j++){
919 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
922 Float_t tmpPt = tmpRec.Pt();
924 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
925 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
926 aodOutJet->GetRefTracks()->Clear();
927 Double_t area1 = clustSeq.area(sortedJets[j]);
928 aodOutJet->SetEffArea(area1,0);
929 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
930 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
931 aodOutJet->SetVectorAreaCharged(&vecareab);
935 Float_t tmpPtBack = 0;
936 if(externalBackground){
937 // carefull has to be filled in a task before
938 // todo, ReArrange to the botom
939 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
941 tmpPt = tmpPt - tmpPtBack;
942 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
944 fh1PtJetsRecIn->Fill(tmpPt);
945 // Fill Spectra with constituentsemacs
946 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
948 fh1NConstRec->Fill(constituents.size());
949 fh2PtNch->Fill(nCh,tmpPt);
950 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
951 fh2NConstPt->Fill(tmpPt,constituents.size());
952 // loop over constiutents and fill spectrum
954 for(unsigned int ic = 0; ic < constituents.size();ic++){
955 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
956 fh1PtJetConstRec->Fill(part->Pt());
958 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
959 if(part->Pt()>fMaxTrackPtInJet)aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
961 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
965 Float_t tmpPhi = tmpRec.Phi();
966 Float_t tmpEta = tmpRec.Eta();
967 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
969 fh1PtJetsLeadingRecIn->Fill(tmpPt);
970 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
971 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
972 fh1NConstLeadingRec->Fill(constituents.size());
973 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
976 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
977 fh2JetEtaPt->Fill(tmpEta,tmpPt);
978 Float_t dPhi = phi - tmpPhi;
979 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
980 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
981 Float_t dEta = eta - tmpRec.Eta();
982 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
983 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
985 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
986 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
988 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
989 }// loop over reconstructed jets
994 // Add the random cones...
995 if(fNRandomCones>0&&fTCARandomConesOut){
996 // create a random jet within the acceptance
997 Double_t etaMax = fTrackEtaWindow - fRparam;
1000 Double_t pTC = 1; // small number
1001 for(int ir = 0;ir < fNRandomCones;ir++){
1002 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1003 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1005 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1006 Double_t pZC = pTC/TMath::Tan(thetaC);
1007 Double_t pXC = pTC * TMath::Cos(phiC);
1008 Double_t pYC = pTC * TMath::Sin(phiC);
1009 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1010 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1012 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1013 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1014 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1019 // test for overlap with previous cones to avoid double counting
1020 for(int iic = 0;iic<ir;iic++){
1021 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1023 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1030 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1031 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1032 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1033 }// loop over random cones creation
1036 // loop over the reconstructed particles and add up the pT in the random cones
1037 // maybe better to loop over randomized particles not in the real jets...
1038 // but this by definition brings dow average energy in the whole event
1039 AliAODJet vTmpRanR(1,0,0,1);
1040 for(int i = 0; i < recParticles.GetEntries(); i++){
1041 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1042 if(fTCARandomConesOut){
1043 for(int ir = 0;ir < fNRandomCones;ir++){
1044 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1045 if(jC&&jC->DeltaR(vp)<fRparam){
1046 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1047 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1050 }// add up energy in cone
1052 // the randomized input changes eta and phi, but keeps the p_T
1053 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1054 Double_t pTR = vp->Pt();
1055 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1056 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1058 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1059 Double_t pZR = pTR/TMath::Tan(thetaR);
1061 Double_t pXR = pTR * TMath::Cos(phiR);
1062 Double_t pYR = pTR * TMath::Sin(phiR);
1063 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1064 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1065 if(fTCARandomConesOutRan){
1066 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1067 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1068 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1069 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1070 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1075 }// loop over recparticles
1077 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1078 if(fTCARandomConesOut){
1079 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1080 // rescale the momntum vectors for the random cones
1082 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1084 Double_t etaC = rC->Eta();
1085 Double_t phiC = rC->Phi();
1086 // massless jet, unit vector
1087 pTC = rC->ChargedBgEnergy();
1088 if(pTC<=0)pTC = 0.001; // for almost empty events
1089 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1090 Double_t pZC = pTC/TMath::Tan(thetaC);
1091 Double_t pXC = pTC * TMath::Cos(phiC);
1092 Double_t pYC = pTC * TMath::Sin(phiC);
1093 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1094 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1095 rC->SetBgEnergy(0,0);
1096 rC->SetEffArea(jetArea,0);
1100 if(fTCARandomConesOutRan){
1101 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1102 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1105 Double_t etaC = rC->Eta();
1106 Double_t phiC = rC->Phi();
1107 // massless jet, unit vector
1108 pTC = rC->ChargedBgEnergy();
1109 if(pTC<=0)pTC = 0.001;// for almost empty events
1110 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1111 Double_t pZC = pTC/TMath::Tan(thetaC);
1112 Double_t pXC = pTC * TMath::Cos(phiC);
1113 Double_t pYC = pTC * TMath::Sin(phiC);
1114 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1115 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1116 rC->SetBgEnergy(0,0);
1117 rC->SetEffArea(jetArea,0);
1121 }// if(fNRandomCones
1123 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1129 if(fAODJetBackgroundOut){
1130 vector<fastjet::PseudoJet> jets2=sortedJets;
1131 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1134 Double_t meanarea1=0.;
1137 Double_t meanarea2=0.;
1139 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1140 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1142 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1143 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1145 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1146 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1147 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1148 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1157 // fill track information
1158 Int_t nTrackOver = recParticles.GetSize();
1159 // do the same for tracks and jets
1162 TIterator *recIter = recParticles.MakeIterator();
1163 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1164 Float_t pt = tmpRec->Pt();
1166 // Printf("Leading track p_t %3.3E",pt);
1167 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1168 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1169 while(pt<ptCut&&tmpRec){
1171 tmpRec = (AliVParticle*)(recIter->Next());
1176 if(nTrackOver<=0)break;
1177 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1181 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1182 Float_t phi = leading->Phi();
1183 if(phi<0)phi+=TMath::Pi()*2.;
1184 Float_t eta = leading->Eta();
1186 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1187 Float_t tmpPt = tmpRec->Pt();
1188 Float_t tmpEta = tmpRec->Eta();
1189 fh1PtTracksRecIn->Fill(tmpPt);
1190 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1191 if(tmpRec==leading){
1192 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1193 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1197 Float_t tmpPhi = tmpRec->Phi();
1199 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1200 Float_t dPhi = phi - tmpPhi;
1201 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1202 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1203 Float_t dEta = eta - tmpRec->Eta();
1204 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1205 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1206 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1211 // find the random jets
1213 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1215 // fill the jet information from random track
1216 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1217 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1219 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1221 // loop over all jets an fill information, first on is the leading jet
1223 Int_t nRecOverRan = inclusiveJetsRan.size();
1224 Int_t nRecRan = inclusiveJetsRan.size();
1226 if(inclusiveJetsRan.size()>0){
1227 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1228 Float_t pt = leadingJet.Pt();
1231 TLorentzVector vecarearanb;
1233 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1234 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1235 while(pt<ptCut&&iCount<nRecRan){
1239 pt = sortedJetsRan[iCount].perp();
1242 if(nRecOverRan<=0)break;
1243 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1245 Float_t phi = leadingJet.Phi();
1246 if(phi<0)phi+=TMath::Pi()*2.;
1247 pt = leadingJet.Pt();
1249 // correlation of leading jet with random tracks
1251 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1253 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1255 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1256 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1257 Float_t dPhi = phi - tmpPhi;
1258 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1259 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1260 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1261 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1264 Int_t nAodOutJetsRan = 0;
1265 AliAODJet *aodOutJetRan = 0;
1266 for(int j = 0; j < nRecRan;j++){
1267 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1268 Float_t tmpPt = tmpRec.Pt();
1269 fh1PtJetsRecInRan->Fill(tmpPt);
1270 // Fill Spectra with constituents
1271 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1272 fh1NConstRecRan->Fill(constituents.size());
1273 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1274 fh2PtNchRan->Fill(nCh,tmpPt);
1275 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1278 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1279 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1280 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1281 aodOutJetRan->GetRefTracks()->Clear();
1282 aodOutJetRan->SetEffArea(arearan,0);
1283 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1284 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1285 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1290 Float_t tmpPhi = tmpRec.Phi();
1291 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1294 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1295 fh1NConstLeadingRecRan->Fill(constituents.size());
1296 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1302 if(fAODJetBackgroundOut){
1305 Double_t meanarea3=0.;
1306 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1307 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1308 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1310 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1311 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1320 // do the event selection if activated
1321 if(fJetTriggerPtCut>0){
1322 bool select = false;
1323 Float_t minPt = fJetTriggerPtCut;
1325 // hard coded for now ...
1326 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1327 if(cent<10)minPt = 50;
1328 else if(cent<30)minPt = 42;
1329 else if(cent<50)minPt = 28;
1330 else if(cent<80)minPt = 18;
1333 if(externalBackground)rho = externalBackground->GetBackground(2);
1335 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1336 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1337 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1346 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1347 fh1CentralitySelect->Fill(cent);
1348 fh1ZSelect->Fill(zVtx);
1349 aodH->SetFillAOD(kTRUE);
1353 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1354 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1355 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1356 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1358 PostData(1, fHistList);
1361 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1364 // Terminate analysis
1366 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1370 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1373 // get list of tracks/particles for different types
1376 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1379 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1380 if(type!=kTrackAODextraonly) {
1381 AliAODEvent *aod = 0;
1382 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1383 else aod = AODEvent();
1385 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1388 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1389 AliAODTrack *tr = aod->GetTrack(it);
1390 Bool_t bGood = false;
1391 if(fFilterType == 0)bGood = true;
1392 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1393 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1394 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1395 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1398 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1399 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1402 if(tr->Pt()<fTrackPtCut){
1403 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1406 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1411 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1412 AliAODEvent *aod = 0;
1413 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1414 else aod = AODEvent();
1419 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1420 if(!aodExtraTracks)return iCount;
1421 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1422 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1423 if (!track) continue;
1425 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1426 if(!trackAOD)continue;
1427 Bool_t bGood = false;
1428 if(fFilterType == 0)bGood = true;
1429 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1430 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1431 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1432 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1433 if(trackAOD->Pt()<fTrackPtCut) continue;
1434 list->Add(trackAOD);
1439 else if (type == kTrackKineAll||type == kTrackKineCharged){
1440 AliMCEvent* mcEvent = MCEvent();
1441 if(!mcEvent)return iCount;
1442 // we want to have alivpartilces so use get track
1443 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1444 if(!mcEvent->IsPhysicalPrimary(it))continue;
1445 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1446 if(type == kTrackKineAll){
1447 if(part->Pt()<fTrackPtCut)continue;
1451 else if(type == kTrackKineCharged){
1452 if(part->Particle()->GetPDG()->Charge()==0)continue;
1453 if(part->Pt()<fTrackPtCut)continue;
1459 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1460 AliAODEvent *aod = 0;
1461 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1462 else aod = AODEvent();
1463 if(!aod)return iCount;
1464 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1465 if(!tca)return iCount;
1466 for(int it = 0;it < tca->GetEntriesFast();++it){
1467 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1468 if(!part->IsPhysicalPrimary())continue;
1469 if(type == kTrackAODMCAll){
1470 if(part->Pt()<fTrackPtCut)continue;
1474 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1475 if(part->Charge()==0)continue;
1476 if(part->Pt()<fTrackPtCut)continue;
1477 if(kTrackAODMCCharged){
1481 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1493 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1494 for(int i = 0; i < particles.GetEntries(); i++){
1495 AliVParticle *vp = (AliVParticle*)particles.At(i);
1496 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1497 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1498 jInp.set_user_index(i);
1499 inputParticles.push_back(jInp);