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(){
75 if(fTCAJetsOut)fTCAJetsOut->Delete();
77 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
78 delete fTCAJetsOutRan;
79 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
80 delete fTCARandomConesOut;
81 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
82 delete fTCARandomConesOutRan;
83 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
84 delete fAODJetBackgroundOut;
87 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
92 fUseAODTrackInput(kFALSE),
93 fUseAODMCInput(kFALSE),
94 fUseBackgroundCalc(kFALSE),
95 fEventSelection(kFALSE),
99 fTrackTypeRec(kTrackUndef),
100 fTrackTypeGen(kTrackUndef),
102 fNSkipLeadingCone(0),
106 fTrackEtaWindow(0.9),
109 fJetOutputMinPt(0.150),
110 fMaxTrackPtInJet(100.),
117 fBackgroundBranch(""),
120 fAlgorithm(fastjet::kt_algorithm),
121 fStrategy(fastjet::Best),
122 fRecombScheme(fastjet::BIpt_scheme),
123 fAreaType(fastjet::active_area),
125 fActiveAreaRepeats(1),
129 fTCARandomConesOut(0x0),
130 fTCARandomConesOutRan(0x0),
131 fAODJetBackgroundOut(0x0),
137 fh1PtHardTrials(0x0),
140 fh1NConstLeadingRec(0x0),
142 fh1PtJetsLeadingRecIn(0x0),
143 fh1PtJetConstRec(0x0),
144 fh1PtJetConstLeadingRec(0x0),
145 fh1PtTracksRecIn(0x0),
146 fh1PtTracksLeadingRecIn(0x0),
148 fh1NConstRecRan(0x0),
149 fh1PtJetsLeadingRecInRan(0x0),
150 fh1NConstLeadingRecRan(0x0),
151 fh1PtJetsRecInRan(0x0),
152 fh1PtTracksGenIn(0x0),
154 fh1CentralityPhySel(0x0),
156 fh1CentralitySelect(0x0),
161 fh2NRecTracksPt(0x0),
163 fh2NConstLeadingPt(0x0),
165 fh2LeadingJetPhiEta(0x0),
167 fh2LeadingJetEtaPt(0x0),
169 fh2LeadingTrackEtaPt(0x0),
170 fh2JetsLeadingPhiEta(0x0),
171 fh2JetsLeadingPhiPt(0x0),
172 fh2TracksLeadingPhiEta(0x0),
173 fh2TracksLeadingPhiPt(0x0),
174 fh2TracksLeadingJetPhiPt(0x0),
175 fh2JetsLeadingPhiPtW(0x0),
176 fh2TracksLeadingPhiPtW(0x0),
177 fh2TracksLeadingJetPhiPtW(0x0),
178 fh2NRecJetsPtRan(0x0),
180 fh2NConstLeadingPtRan(0x0),
185 fh2TracksLeadingJetPhiPtRan(0x0),
186 fh2TracksLeadingJetPhiPtWRan(0x0),
189 for(int i = 0;i<3;i++){
190 fh1BiARandomCones[i] = 0;
191 fh1BiARandomConesRan[i] = 0;
193 for(int i = 0;i<kMaxCent;i++){
194 fh2JetsLeadingPhiPtC[i] = 0;
195 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
196 fh2TracksLeadingJetPhiPtC[i] = 0;
197 fh2TracksLeadingJetPhiPtWC[i] = 0;
201 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
202 AliAnalysisTaskSE(name),
206 fUseAODTrackInput(kFALSE),
207 fUseAODMCInput(kFALSE),
208 fUseBackgroundCalc(kFALSE),
209 fEventSelection(kFALSE),
213 fTrackTypeRec(kTrackUndef),
214 fTrackTypeGen(kTrackUndef),
216 fNSkipLeadingCone(0),
220 fTrackEtaWindow(0.9),
223 fJetOutputMinPt(0.150),
224 fMaxTrackPtInJet(100.),
231 fBackgroundBranch(""),
234 fAlgorithm(fastjet::kt_algorithm),
235 fStrategy(fastjet::Best),
236 fRecombScheme(fastjet::BIpt_scheme),
237 fAreaType(fastjet::active_area),
239 fActiveAreaRepeats(1),
243 fTCARandomConesOut(0x0),
244 fTCARandomConesOutRan(0x0),
245 fAODJetBackgroundOut(0x0),
251 fh1PtHardTrials(0x0),
254 fh1NConstLeadingRec(0x0),
256 fh1PtJetsLeadingRecIn(0x0),
257 fh1PtJetConstRec(0x0),
258 fh1PtJetConstLeadingRec(0x0),
259 fh1PtTracksRecIn(0x0),
260 fh1PtTracksLeadingRecIn(0x0),
262 fh1NConstRecRan(0x0),
263 fh1PtJetsLeadingRecInRan(0x0),
264 fh1NConstLeadingRecRan(0x0),
265 fh1PtJetsRecInRan(0x0),
266 fh1PtTracksGenIn(0x0),
268 fh1CentralityPhySel(0x0),
270 fh1CentralitySelect(0x0),
275 fh2NRecTracksPt(0x0),
277 fh2NConstLeadingPt(0x0),
279 fh2LeadingJetPhiEta(0x0),
281 fh2LeadingJetEtaPt(0x0),
283 fh2LeadingTrackEtaPt(0x0),
284 fh2JetsLeadingPhiEta(0x0),
285 fh2JetsLeadingPhiPt(0x0),
286 fh2TracksLeadingPhiEta(0x0),
287 fh2TracksLeadingPhiPt(0x0),
288 fh2TracksLeadingJetPhiPt(0x0),
289 fh2JetsLeadingPhiPtW(0x0),
290 fh2TracksLeadingPhiPtW(0x0),
291 fh2TracksLeadingJetPhiPtW(0x0),
292 fh2NRecJetsPtRan(0x0),
294 fh2NConstLeadingPtRan(0x0),
299 fh2TracksLeadingJetPhiPtRan(0x0),
300 fh2TracksLeadingJetPhiPtWRan(0x0),
303 for(int i = 0;i<3;i++){
304 fh1BiARandomCones[i] = 0;
305 fh1BiARandomConesRan[i] = 0;
307 for(int i = 0;i<kMaxCent;i++){
308 fh2JetsLeadingPhiPtC[i] = 0;
309 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
310 fh2TracksLeadingJetPhiPtC[i] = 0;
311 fh2TracksLeadingJetPhiPtWC[i] = 0;
313 DefineOutput(1, TList::Class());
318 Bool_t AliAnalysisTaskJetCluster::Notify()
321 // Implemented Notify() to read the cross sections
322 // and number of trials from pyxsec.root
327 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
331 // Create the output container
334 fRandom = new TRandom3(0);
340 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
344 if(fNonStdBranch.Length()!=0)
346 // only create the output branch if we have a name
347 // Create a new branch for jets...
348 // -> cleared in the UserExec....
349 // here we can also have the case that the brnaches are written to a separate file
352 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
353 fTCAJetsOut->SetName(fNonStdBranch.Data());
354 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
357 if(fJetTypes&kJetRan){
358 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
359 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
360 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
363 if(fUseBackgroundCalc){
364 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
365 fAODJetBackgroundOut = new AliAODJetEventBackground();
366 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
367 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
370 // create the branch for the random cones with the same R
371 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
375 if(!AODEvent()->FindListObject(cName.Data())){
376 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
377 fTCARandomConesOut->SetName(cName.Data());
378 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
381 // create the branch with the random for the random cones on the random event
382 if(fJetTypes&kRCRan){
383 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
384 if(!AODEvent()->FindListObject(cName.Data())){
385 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
386 fTCARandomConesOutRan->SetName(cName.Data());
387 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
392 if(fNonStdFile.Length()!=0){
394 // case that we have an AOD extension we need to fetch the jets from the extended output
395 // we identify the extension aod event by looking for the branchname
396 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
397 // case that we have an AOD extension we need can fetch the background maybe from the extended output
398 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
403 if(!fHistList)fHistList = new TList();
404 fHistList->SetOwner();
405 PostData(1, fHistList); // post data in any case once
407 Bool_t oldStatus = TH1::AddDirectoryStatus();
408 TH1::AddDirectory(kFALSE);
413 const Int_t nBinPt = 100;
414 Double_t binLimitsPt[nBinPt+1];
415 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
417 binLimitsPt[iPt] = 0.0;
420 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
424 const Int_t nBinPhi = 90;
425 Double_t binLimitsPhi[nBinPhi+1];
426 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
428 binLimitsPhi[iPhi] = -1.*TMath::Pi();
431 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
437 const Int_t nBinEta = 40;
438 Double_t binLimitsEta[nBinEta+1];
439 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
441 binLimitsEta[iEta] = -2.0;
444 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
448 const int nChMax = 5000;
450 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
451 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
453 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
454 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
457 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
458 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
460 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
461 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
462 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
463 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
466 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
467 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
468 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
470 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
471 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
472 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
473 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
474 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
475 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
476 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
477 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
478 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
479 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
481 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
482 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
483 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
485 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
486 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
487 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
489 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
490 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
491 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
495 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
496 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
497 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
498 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
500 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
501 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);
502 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);
503 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);
507 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
508 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
509 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
510 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
512 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
513 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
514 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
515 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
517 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
518 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
519 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
520 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
524 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
525 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
526 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
527 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
528 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
529 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
530 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
531 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
532 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
533 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
534 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
535 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
537 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
538 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
539 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
540 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
542 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
543 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
544 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
545 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
548 if(fNRandomCones>0&&fUseBackgroundCalc){
549 for(int i = 0;i<3;i++){
550 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
551 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
555 for(int i = 0;i < kMaxCent;i++){
556 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
557 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
558 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
559 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
562 const Int_t saveLevel = 3; // large save level more histos
564 fHistList->Add(fh1Xsec);
565 fHistList->Add(fh1Trials);
567 fHistList->Add(fh1NJetsRec);
568 fHistList->Add(fh1NConstRec);
569 fHistList->Add(fh1NConstLeadingRec);
570 fHistList->Add(fh1PtJetsRecIn);
571 fHistList->Add(fh1PtJetsLeadingRecIn);
572 fHistList->Add(fh1PtTracksRecIn);
573 fHistList->Add(fh1PtTracksLeadingRecIn);
574 fHistList->Add(fh1PtJetConstRec);
575 fHistList->Add(fh1PtJetConstLeadingRec);
576 fHistList->Add(fh1NJetsRecRan);
577 fHistList->Add(fh1NConstRecRan);
578 fHistList->Add(fh1PtJetsLeadingRecInRan);
579 fHistList->Add(fh1NConstLeadingRecRan);
580 fHistList->Add(fh1PtJetsRecInRan);
581 fHistList->Add(fh1Nch);
582 fHistList->Add(fh1Centrality);
583 fHistList->Add(fh1CentralitySelect);
584 fHistList->Add(fh1CentralityPhySel);
585 fHistList->Add(fh1Z);
586 fHistList->Add(fh1ZSelect);
587 fHistList->Add(fh1ZPhySel);
588 if(fNRandomCones>0&&fUseBackgroundCalc){
589 for(int i = 0;i<3;i++){
590 fHistList->Add(fh1BiARandomCones[i]);
591 fHistList->Add(fh1BiARandomConesRan[i]);
594 for(int i = 0;i < kMaxCent;i++){
595 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
596 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
597 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
598 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
601 fHistList->Add(fh2NRecJetsPt);
602 fHistList->Add(fh2NRecTracksPt);
603 fHistList->Add(fh2NConstPt);
604 fHistList->Add(fh2NConstLeadingPt);
605 fHistList->Add(fh2PtNch);
606 fHistList->Add(fh2PtNchRan);
607 fHistList->Add(fh2PtNchN);
608 fHistList->Add(fh2PtNchNRan);
609 fHistList->Add(fh2JetPhiEta);
610 fHistList->Add(fh2LeadingJetPhiEta);
611 fHistList->Add(fh2JetEtaPt);
612 fHistList->Add(fh2LeadingJetEtaPt);
613 fHistList->Add(fh2TrackEtaPt);
614 fHistList->Add(fh2LeadingTrackEtaPt);
615 fHistList->Add(fh2JetsLeadingPhiEta );
616 fHistList->Add(fh2JetsLeadingPhiPt);
617 fHistList->Add(fh2TracksLeadingPhiEta);
618 fHistList->Add(fh2TracksLeadingPhiPt);
619 fHistList->Add(fh2TracksLeadingJetPhiPt);
620 fHistList->Add(fh2JetsLeadingPhiPtW);
621 fHistList->Add(fh2TracksLeadingPhiPtW);
622 fHistList->Add(fh2TracksLeadingJetPhiPtW);
623 fHistList->Add(fh2NRecJetsPtRan);
624 fHistList->Add(fh2NConstPtRan);
625 fHistList->Add(fh2NConstLeadingPtRan);
626 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
627 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
630 // =========== Switch on Sumw2 for all histos ===========
631 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
632 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
637 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
640 TH1::AddDirectory(oldStatus);
643 void AliAnalysisTaskJetCluster::Init()
649 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
653 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
656 // handle and reset the output jet branch
658 if(fTCAJetsOut)fTCAJetsOut->Delete();
659 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
660 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
661 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
662 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
664 AliAODJetEventBackground* externalBackground = 0;
665 if(!externalBackground&&fBackgroundBranch.Length()){
666 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
667 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
668 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
671 // Execute analysis for current event
673 AliESDEvent *fESD = 0;
674 if(fUseAODTrackInput){
675 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
677 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
683 // assume that the AOD is in the general output...
686 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
690 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
694 Bool_t selectEvent = false;
695 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
701 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
702 TString vtxTitle(vtxAOD->GetTitle());
703 zVtx = vtxAOD->GetZ();
705 cent = fAOD->GetHeader()->GetCentrality();
706 if(cent<10)cenClass = 0;
707 else if(cent<30)cenClass = 1;
708 else if(cent<50)cenClass = 2;
709 else if(cent<80)cenClass = 3;
710 if(physicsSelection){
711 fh1CentralityPhySel->Fill(cent);
712 fh1ZPhySel->Fill(zVtx);
716 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
717 Float_t yvtx = vtxAOD->GetY();
718 Float_t xvtx = vtxAOD->GetX();
719 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
720 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
721 if(physicsSelection){
727 if(cent<fCentCutLo||cent>fCentCutUp){
738 PostData(1, fHistList);
741 fh1Centrality->Fill(cent);
743 fh1Trials->Fill("#sum{ntrials}",1);
746 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
748 // ==== General variables needed
752 // we simply fetch the tracks/mc particles as a list of AliVParticles
757 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
758 Float_t nCh = recParticles.GetEntries();
760 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
761 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
762 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
766 vector<fastjet::PseudoJet> inputParticlesRec;
767 vector<fastjet::PseudoJet> inputParticlesRecRan;
769 // Generate the random cones
771 AliAODJet vTmpRan(1,0,0,1);
772 for(int i = 0; i < recParticles.GetEntries(); i++){
773 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
774 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
775 // we take total momentum here
776 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
777 jInp.set_user_index(i);
778 inputParticlesRec.push_back(jInp);
780 // the randomized input changes eta and phi, but keeps the p_T
781 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
782 Double_t pT = vp->Pt();
783 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
784 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
786 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
787 Double_t pZ = pT/TMath::Tan(theta);
789 Double_t pX = pT * TMath::Cos(phi);
790 Double_t pY = pT * TMath::Sin(phi);
791 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
792 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
794 jInpRan.set_user_index(i);
795 inputParticlesRecRan.push_back(jInpRan);
796 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
799 // fill the tref array, only needed when we write out jets
802 fRef->Delete(); // make sure to delete before placement new...
803 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
809 if(inputParticlesRec.size()==0){
810 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
811 PostData(1, fHistList);
816 // employ setters for these...
819 // now create the object that holds info about ghosts
821 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
822 // reduce CPU time...
824 fActiveAreaRepeats = 0;
828 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
829 fastjet::AreaType areaType = fastjet::active_area;
830 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
831 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
832 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
834 //range where to compute background
835 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
837 phiMax = 2*TMath::Pi();
838 rapMax = fGhostEtamax - fRparam;
839 rapMin = - fGhostEtamax + fRparam;
840 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
843 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
844 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
847 fh1NJetsRec->Fill(sortedJets.size());
849 // loop over all jets an fill information, first on is the leading jet
851 Int_t nRecOver = inclusiveJets.size();
852 Int_t nRec = inclusiveJets.size();
853 if(inclusiveJets.size()>0){
854 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
855 Double_t area = clustSeq.area(sortedJets[0]);
856 leadingJet.SetEffArea(area,0);
857 Float_t pt = leadingJet.Pt();
858 Int_t nAodOutJets = 0;
859 Int_t nAodOutTracks = 0;
860 AliAODJet *aodOutJet = 0;
863 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
864 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
865 while(pt<ptCut&&iCount<nRec){
869 pt = sortedJets[iCount].perp();
872 if(nRecOver<=0)break;
873 fh2NRecJetsPt->Fill(ptCut,nRecOver);
875 Float_t phi = leadingJet.Phi();
876 if(phi<0)phi+=TMath::Pi()*2.;
877 Float_t eta = leadingJet.Eta();
879 if(externalBackground){
880 // carefull has to be filled in a task before
881 // todo, ReArrange to the botom
882 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
884 pt = leadingJet.Pt() - pTback;
885 // correlation of leading jet with tracks
886 TIterator *recIter = recParticles.MakeIterator();
888 AliVParticle *tmpRecTrack = 0;
889 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
890 Float_t tmpPt = tmpRecTrack->Pt();
892 Float_t tmpPhi = tmpRecTrack->Phi();
893 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
894 Float_t dPhi = phi - tmpPhi;
895 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
896 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
897 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
898 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
900 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
901 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
907 TLorentzVector vecareab;
908 for(int j = 0; j < nRec;j++){
909 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
912 Float_t tmpPt = tmpRec.Pt();
914 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
915 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
916 aodOutJet->GetRefTracks()->Clear();
917 Double_t area1 = clustSeq.area(sortedJets[j]);
918 aodOutJet->SetEffArea(area1,0);
919 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
920 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
921 aodOutJet->SetVectorAreaCharged(&vecareab);
925 Float_t tmpPtBack = 0;
926 if(externalBackground){
927 // carefull has to be filled in a task before
928 // todo, ReArrange to the botom
929 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
931 tmpPt = tmpPt - tmpPtBack;
932 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
934 fh1PtJetsRecIn->Fill(tmpPt);
935 // Fill Spectra with constituentsemacs
936 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
938 fh1NConstRec->Fill(constituents.size());
939 fh2PtNch->Fill(nCh,tmpPt);
940 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
941 fh2NConstPt->Fill(tmpPt,constituents.size());
942 // loop over constiutents and fill spectrum
944 for(unsigned int ic = 0; ic < constituents.size();ic++){
945 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
946 fh1PtJetConstRec->Fill(part->Pt());
948 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
949 if(part->Pt()>fMaxTrackPtInJet)aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
951 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
955 Float_t tmpPhi = tmpRec.Phi();
956 Float_t tmpEta = tmpRec.Eta();
957 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
959 fh1PtJetsLeadingRecIn->Fill(tmpPt);
960 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
961 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
962 fh1NConstLeadingRec->Fill(constituents.size());
963 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
966 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
967 fh2JetEtaPt->Fill(tmpEta,tmpPt);
968 Float_t dPhi = phi - tmpPhi;
969 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
970 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
971 Float_t dEta = eta - tmpRec.Eta();
972 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
973 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
975 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
976 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
978 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
979 }// loop over reconstructed jets
984 // Add the random cones...
985 if(fNRandomCones>0&&fTCARandomConesOut){
986 // create a random jet within the acceptance
987 Double_t etaMax = fTrackEtaWindow - fRparam;
990 Double_t pTC = 1; // small number
991 for(int ir = 0;ir < fNRandomCones;ir++){
992 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
993 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
995 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
996 Double_t pZC = pTC/TMath::Tan(thetaC);
997 Double_t pXC = pTC * TMath::Cos(phiC);
998 Double_t pYC = pTC * TMath::Sin(phiC);
999 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1000 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1002 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1003 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1004 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1009 // test for overlap with previous cones to avoid double counting
1010 for(int iic = 0;iic<ir;iic++){
1011 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1013 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1020 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1021 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1022 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1023 }// loop over random cones creation
1026 // loop over the reconstructed particles and add up the pT in the random cones
1027 // maybe better to loop over randomized particles not in the real jets...
1028 // but this by definition brings dow average energy in the whole event
1029 AliAODJet vTmpRanR(1,0,0,1);
1030 for(int i = 0; i < recParticles.GetEntries(); i++){
1031 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1032 if(fTCARandomConesOut){
1033 for(int ir = 0;ir < fNRandomCones;ir++){
1034 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1035 if(jC&&jC->DeltaR(vp)<fRparam){
1036 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1037 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1040 }// add up energy in cone
1042 // the randomized input changes eta and phi, but keeps the p_T
1043 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1044 Double_t pTR = vp->Pt();
1045 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1046 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1048 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1049 Double_t pZR = pTR/TMath::Tan(thetaR);
1051 Double_t pXR = pTR * TMath::Cos(phiR);
1052 Double_t pYR = pTR * TMath::Sin(phiR);
1053 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1054 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1055 if(fTCARandomConesOutRan){
1056 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1057 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1058 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1059 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1060 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1065 }// loop over recparticles
1067 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1068 if(fTCARandomConesOut){
1069 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1070 // rescale the momntum vectors for the random cones
1072 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1074 Double_t etaC = rC->Eta();
1075 Double_t phiC = rC->Phi();
1076 // massless jet, unit vector
1077 pTC = rC->ChargedBgEnergy();
1078 if(pTC<=0)pTC = 0.001; // for almost empty events
1079 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1080 Double_t pZC = pTC/TMath::Tan(thetaC);
1081 Double_t pXC = pTC * TMath::Cos(phiC);
1082 Double_t pYC = pTC * TMath::Sin(phiC);
1083 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1084 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1085 rC->SetBgEnergy(0,0);
1086 rC->SetEffArea(jetArea,0);
1090 if(fTCARandomConesOutRan){
1091 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1092 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1095 Double_t etaC = rC->Eta();
1096 Double_t phiC = rC->Phi();
1097 // massless jet, unit vector
1098 pTC = rC->ChargedBgEnergy();
1099 if(pTC<=0)pTC = 0.001;// for almost empty events
1100 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1101 Double_t pZC = pTC/TMath::Tan(thetaC);
1102 Double_t pXC = pTC * TMath::Cos(phiC);
1103 Double_t pYC = pTC * TMath::Sin(phiC);
1104 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1105 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1106 rC->SetBgEnergy(0,0);
1107 rC->SetEffArea(jetArea,0);
1111 }// if(fNRandomCones
1113 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1119 if(fAODJetBackgroundOut){
1120 vector<fastjet::PseudoJet> jets2=sortedJets;
1121 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1124 Double_t meanarea1=0.;
1127 Double_t meanarea2=0.;
1129 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1130 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1132 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1133 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1135 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1136 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1137 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1138 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1147 // fill track information
1148 Int_t nTrackOver = recParticles.GetSize();
1149 // do the same for tracks and jets
1152 TIterator *recIter = recParticles.MakeIterator();
1153 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1154 Float_t pt = tmpRec->Pt();
1156 // Printf("Leading track p_t %3.3E",pt);
1157 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1158 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1159 while(pt<ptCut&&tmpRec){
1161 tmpRec = (AliVParticle*)(recIter->Next());
1166 if(nTrackOver<=0)break;
1167 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1171 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1172 Float_t phi = leading->Phi();
1173 if(phi<0)phi+=TMath::Pi()*2.;
1174 Float_t eta = leading->Eta();
1176 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1177 Float_t tmpPt = tmpRec->Pt();
1178 Float_t tmpEta = tmpRec->Eta();
1179 fh1PtTracksRecIn->Fill(tmpPt);
1180 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1181 if(tmpRec==leading){
1182 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1183 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1187 Float_t tmpPhi = tmpRec->Phi();
1189 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1190 Float_t dPhi = phi - tmpPhi;
1191 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1192 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1193 Float_t dEta = eta - tmpRec->Eta();
1194 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1195 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1196 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1201 // find the random jets
1203 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1205 // fill the jet information from random track
1206 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1207 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1209 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1211 // loop over all jets an fill information, first on is the leading jet
1213 Int_t nRecOverRan = inclusiveJetsRan.size();
1214 Int_t nRecRan = inclusiveJetsRan.size();
1216 if(inclusiveJetsRan.size()>0){
1217 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1218 Float_t pt = leadingJet.Pt();
1221 TLorentzVector vecarearanb;
1223 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1224 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1225 while(pt<ptCut&&iCount<nRecRan){
1229 pt = sortedJetsRan[iCount].perp();
1232 if(nRecOverRan<=0)break;
1233 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1235 Float_t phi = leadingJet.Phi();
1236 if(phi<0)phi+=TMath::Pi()*2.;
1237 pt = leadingJet.Pt();
1239 // correlation of leading jet with random tracks
1241 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1243 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1245 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1246 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1247 Float_t dPhi = phi - tmpPhi;
1248 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1249 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1250 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1251 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1254 Int_t nAodOutJetsRan = 0;
1255 AliAODJet *aodOutJetRan = 0;
1256 for(int j = 0; j < nRecRan;j++){
1257 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1258 Float_t tmpPt = tmpRec.Pt();
1259 fh1PtJetsRecInRan->Fill(tmpPt);
1260 // Fill Spectra with constituents
1261 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1262 fh1NConstRecRan->Fill(constituents.size());
1263 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1264 fh2PtNchRan->Fill(nCh,tmpPt);
1265 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1268 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1269 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1270 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1271 aodOutJetRan->GetRefTracks()->Clear();
1272 aodOutJetRan->SetEffArea(arearan,0);
1273 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1274 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1275 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1280 Float_t tmpPhi = tmpRec.Phi();
1281 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1284 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1285 fh1NConstLeadingRecRan->Fill(constituents.size());
1286 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1292 if(fAODJetBackgroundOut){
1295 Double_t meanarea3=0.;
1296 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1297 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1298 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1300 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1301 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1310 // do the event selection if activated
1311 if(fJetTriggerPtCut>0){
1312 bool select = false;
1313 Float_t minPt = fJetTriggerPtCut;
1315 // hard coded for now ...
1316 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1317 if(cent<10)minPt = 50;
1318 else if(cent<30)minPt = 42;
1319 else if(cent<50)minPt = 28;
1320 else if(cent<80)minPt = 18;
1323 if(externalBackground)rho = externalBackground->GetBackground(2);
1325 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1326 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1327 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1336 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1337 fh1CentralitySelect->Fill(cent);
1338 fh1ZSelect->Fill(zVtx);
1339 aodH->SetFillAOD(kTRUE);
1343 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1344 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1345 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1346 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1348 PostData(1, fHistList);
1351 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1353 // Terminate analysis
1355 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1359 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1361 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1364 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1365 if(type!=kTrackAODextraonly) {
1366 AliAODEvent *aod = 0;
1367 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1368 else aod = AODEvent();
1370 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1373 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1374 AliAODTrack *tr = aod->GetTrack(it);
1375 Bool_t bGood = false;
1376 if(fFilterType == 0)bGood = true;
1377 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1378 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1379 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1380 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1383 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1384 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1387 if(tr->Pt()<fTrackPtCut){
1388 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1391 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1396 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1397 AliAODEvent *aod = 0;
1398 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1399 else aod = AODEvent();
1404 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1405 if(!aodExtraTracks)return iCount;
1406 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1407 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1408 if (!track) continue;
1410 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1411 if(!trackAOD)continue;
1412 Bool_t bGood = false;
1413 if(fFilterType == 0)bGood = true;
1414 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1415 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1416 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1417 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1418 if(trackAOD->Pt()<fTrackPtCut) continue;
1419 list->Add(trackAOD);
1424 else if (type == kTrackKineAll||type == kTrackKineCharged){
1425 AliMCEvent* mcEvent = MCEvent();
1426 if(!mcEvent)return iCount;
1427 // we want to have alivpartilces so use get track
1428 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1429 if(!mcEvent->IsPhysicalPrimary(it))continue;
1430 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1431 if(type == kTrackKineAll){
1432 if(part->Pt()<fTrackPtCut)continue;
1436 else if(type == kTrackKineCharged){
1437 if(part->Particle()->GetPDG()->Charge()==0)continue;
1438 if(part->Pt()<fTrackPtCut)continue;
1444 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1445 AliAODEvent *aod = 0;
1446 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1447 else aod = AODEvent();
1448 if(!aod)return iCount;
1449 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1450 if(!tca)return iCount;
1451 for(int it = 0;it < tca->GetEntriesFast();++it){
1452 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1453 if(!part->IsPhysicalPrimary())continue;
1454 if(type == kTrackAODMCAll){
1455 if(part->Pt()<fTrackPtCut)continue;
1459 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1460 if(part->Charge()==0)continue;
1461 if(part->Pt()<fTrackPtCut)continue;
1462 if(kTrackAODMCCharged){
1466 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1478 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1479 for(int i = 0; i < particles.GetEntries(); i++){
1480 AliVParticle *vp = (AliVParticle*)particles.At(i);
1481 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1482 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1483 jInp.set_user_index(i);
1484 inputParticles.push_back(jInp);