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(): AliAnalysisTaskSE(),
91 fUseAODTrackInput(kFALSE),
92 fUseAODMCInput(kFALSE),
93 fUseBackgroundCalc(kFALSE),
94 fEventSelection(kFALSE),
96 fTrackTypeRec(kTrackUndef),
97 fTrackTypeGen(kTrackUndef),
103 fTrackEtaWindow(0.9),
106 fJetOutputMinPt(0.150),
113 fBackgroundBranch(""),
116 fAlgorithm(fastjet::kt_algorithm),
117 fStrategy(fastjet::Best),
118 fRecombScheme(fastjet::BIpt_scheme),
119 fAreaType(fastjet::active_area),
121 fActiveAreaRepeats(1),
125 fTCARandomConesOut(0x0),
126 fTCARandomConesOutRan(0x0),
127 fAODJetBackgroundOut(0x0),
133 fh1PtHardTrials(0x0),
136 fh1NConstLeadingRec(0x0),
138 fh1PtJetsLeadingRecIn(0x0),
139 fh1PtJetConstRec(0x0),
140 fh1PtJetConstLeadingRec(0x0),
141 fh1PtTracksRecIn(0x0),
142 fh1PtTracksLeadingRecIn(0x0),
144 fh1NConstRecRan(0x0),
145 fh1PtJetsLeadingRecInRan(0x0),
146 fh1NConstLeadingRecRan(0x0),
147 fh1PtJetsRecInRan(0x0),
148 fh1PtTracksGenIn(0x0),
150 fh1CentralityPhySel(0x0),
152 fh1CentralitySelect(0x0),
157 fh2NRecTracksPt(0x0),
159 fh2NConstLeadingPt(0x0),
161 fh2LeadingJetPhiEta(0x0),
163 fh2LeadingJetEtaPt(0x0),
165 fh2LeadingTrackEtaPt(0x0),
166 fh2JetsLeadingPhiEta(0x0),
167 fh2JetsLeadingPhiPt(0x0),
168 fh2TracksLeadingPhiEta(0x0),
169 fh2TracksLeadingPhiPt(0x0),
170 fh2TracksLeadingJetPhiPt(0x0),
171 fh2JetsLeadingPhiPtW(0x0),
172 fh2TracksLeadingPhiPtW(0x0),
173 fh2TracksLeadingJetPhiPtW(0x0),
174 fh2NRecJetsPtRan(0x0),
176 fh2NConstLeadingPtRan(0x0),
181 fh2TracksLeadingJetPhiPtRan(0x0),
182 fh2TracksLeadingJetPhiPtWRan(0x0),
185 for(int i = 0;i<3;i++){
186 fh1BiARandomCones[i] = 0;
187 fh1BiARandomConesRan[i] = 0;
189 for(int i = 0;i<kMaxCent;i++){
190 fh2JetsLeadingPhiPtC[i] = 0;
191 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
192 fh2TracksLeadingJetPhiPtC[i] = 0;
193 fh2TracksLeadingJetPhiPtWC[i] = 0;
197 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
198 AliAnalysisTaskSE(name),
202 fUseAODTrackInput(kFALSE),
203 fUseAODMCInput(kFALSE),
204 fUseBackgroundCalc(kFALSE),
205 fEventSelection(kFALSE),
207 fTrackTypeRec(kTrackUndef),
208 fTrackTypeGen(kTrackUndef),
210 fNSkipLeadingCone(0),
214 fTrackEtaWindow(0.9),
217 fJetOutputMinPt(0.150),
224 fBackgroundBranch(""),
227 fAlgorithm(fastjet::kt_algorithm),
228 fStrategy(fastjet::Best),
229 fRecombScheme(fastjet::BIpt_scheme),
230 fAreaType(fastjet::active_area),
232 fActiveAreaRepeats(1),
236 fTCARandomConesOut(0x0),
237 fTCARandomConesOutRan(0x0),
238 fAODJetBackgroundOut(0x0),
244 fh1PtHardTrials(0x0),
247 fh1NConstLeadingRec(0x0),
249 fh1PtJetsLeadingRecIn(0x0),
250 fh1PtJetConstRec(0x0),
251 fh1PtJetConstLeadingRec(0x0),
252 fh1PtTracksRecIn(0x0),
253 fh1PtTracksLeadingRecIn(0x0),
255 fh1NConstRecRan(0x0),
256 fh1PtJetsLeadingRecInRan(0x0),
257 fh1NConstLeadingRecRan(0x0),
258 fh1PtJetsRecInRan(0x0),
259 fh1PtTracksGenIn(0x0),
261 fh1CentralityPhySel(0x0),
263 fh1CentralitySelect(0x0),
268 fh2NRecTracksPt(0x0),
270 fh2NConstLeadingPt(0x0),
272 fh2LeadingJetPhiEta(0x0),
274 fh2LeadingJetEtaPt(0x0),
276 fh2LeadingTrackEtaPt(0x0),
277 fh2JetsLeadingPhiEta(0x0),
278 fh2JetsLeadingPhiPt(0x0),
279 fh2TracksLeadingPhiEta(0x0),
280 fh2TracksLeadingPhiPt(0x0),
281 fh2TracksLeadingJetPhiPt(0x0),
282 fh2JetsLeadingPhiPtW(0x0),
283 fh2TracksLeadingPhiPtW(0x0),
284 fh2TracksLeadingJetPhiPtW(0x0),
285 fh2NRecJetsPtRan(0x0),
287 fh2NConstLeadingPtRan(0x0),
292 fh2TracksLeadingJetPhiPtRan(0x0),
293 fh2TracksLeadingJetPhiPtWRan(0x0),
296 for(int i = 0;i<3;i++){
297 fh1BiARandomCones[i] = 0;
298 fh1BiARandomConesRan[i] = 0;
300 for(int i = 0;i<kMaxCent;i++){
301 fh2JetsLeadingPhiPtC[i] = 0;
302 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
303 fh2TracksLeadingJetPhiPtC[i] = 0;
304 fh2TracksLeadingJetPhiPtWC[i] = 0;
306 DefineOutput(1, TList::Class());
311 Bool_t AliAnalysisTaskJetCluster::Notify()
314 // Implemented Notify() to read the cross sections
315 // and number of trials from pyxsec.root
320 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
324 // Create the output container
327 fRandom = new TRandom3(0);
333 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
337 if(fNonStdBranch.Length()!=0)
339 // only create the output branch if we have a name
340 // Create a new branch for jets...
341 // -> cleared in the UserExec....
342 // here we can also have the case that the brnaches are written to a separate file
344 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
345 fTCAJetsOut->SetName(fNonStdBranch.Data());
346 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
350 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
351 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
352 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
354 if(fUseBackgroundCalc){
355 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
356 fAODJetBackgroundOut = new AliAODJetEventBackground();
357 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
358 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
361 // create the branch for the random cones with the same R
362 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
365 if(!AODEvent()->FindListObject(cName.Data())){
366 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
367 fTCARandomConesOut->SetName(cName.Data());
368 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
370 // create the branch with the random for the random cones on the random event
371 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
372 if(!AODEvent()->FindListObject(cName.Data())){
373 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
374 fTCARandomConesOutRan->SetName(cName.Data());
375 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
379 if(fNonStdFile.Length()!=0){
381 // case that we have an AOD extension we need to fetch the jets from the extended output
382 // we identify the extension aod event by looking for the branchname
383 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
384 // case that we have an AOD extension we need can fetch the background maybe from the extended output
385 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
390 if(!fHistList)fHistList = new TList();
391 fHistList->SetOwner();
392 PostData(1, fHistList); // post data in any case once
394 Bool_t oldStatus = TH1::AddDirectoryStatus();
395 TH1::AddDirectory(kFALSE);
400 const Int_t nBinPt = 100;
401 Double_t binLimitsPt[nBinPt+1];
402 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
404 binLimitsPt[iPt] = 0.0;
407 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
411 const Int_t nBinPhi = 90;
412 Double_t binLimitsPhi[nBinPhi+1];
413 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
415 binLimitsPhi[iPhi] = -1.*TMath::Pi();
418 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
424 const Int_t nBinEta = 40;
425 Double_t binLimitsEta[nBinEta+1];
426 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
428 binLimitsEta[iEta] = -2.0;
431 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
435 const int nChMax = 5000;
437 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
438 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
440 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
441 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
444 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
445 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
447 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
448 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
449 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
450 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
453 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
454 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
455 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
457 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
458 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
459 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
460 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
461 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
462 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
463 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
464 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
465 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
466 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
468 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
469 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
470 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
472 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
473 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
474 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
476 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
477 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
478 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
482 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
483 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
484 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
485 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
487 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
488 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);
489 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);
490 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);
494 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
495 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
496 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
497 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
499 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
500 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
501 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
502 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
504 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
505 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
506 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
507 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
511 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
512 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
513 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
514 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
515 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
516 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
517 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
518 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
519 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
520 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
521 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
522 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
524 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
525 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
526 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
527 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
529 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
530 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
531 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
532 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
535 if(fNRandomCones>0&&fUseBackgroundCalc){
536 for(int i = 0;i<3;i++){
537 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
538 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
542 for(int i = 0;i < kMaxCent;i++){
543 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
544 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
545 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
546 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
549 const Int_t saveLevel = 3; // large save level more histos
551 fHistList->Add(fh1Xsec);
552 fHistList->Add(fh1Trials);
554 fHistList->Add(fh1NJetsRec);
555 fHistList->Add(fh1NConstRec);
556 fHistList->Add(fh1NConstLeadingRec);
557 fHistList->Add(fh1PtJetsRecIn);
558 fHistList->Add(fh1PtJetsLeadingRecIn);
559 fHistList->Add(fh1PtTracksRecIn);
560 fHistList->Add(fh1PtTracksLeadingRecIn);
561 fHistList->Add(fh1PtJetConstRec);
562 fHistList->Add(fh1PtJetConstLeadingRec);
563 fHistList->Add(fh1NJetsRecRan);
564 fHistList->Add(fh1NConstRecRan);
565 fHistList->Add(fh1PtJetsLeadingRecInRan);
566 fHistList->Add(fh1NConstLeadingRecRan);
567 fHistList->Add(fh1PtJetsRecInRan);
568 fHistList->Add(fh1Nch);
569 fHistList->Add(fh1Centrality);
570 fHistList->Add(fh1CentralitySelect);
571 fHistList->Add(fh1CentralityPhySel);
572 fHistList->Add(fh1Z);
573 fHistList->Add(fh1ZSelect);
574 fHistList->Add(fh1ZPhySel);
575 if(fNRandomCones>0&&fUseBackgroundCalc){
576 for(int i = 0;i<3;i++){
577 fHistList->Add(fh1BiARandomCones[i]);
578 fHistList->Add(fh1BiARandomConesRan[i]);
581 for(int i = 0;i < kMaxCent;i++){
582 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
583 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
584 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
585 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
588 fHistList->Add(fh2NRecJetsPt);
589 fHistList->Add(fh2NRecTracksPt);
590 fHistList->Add(fh2NConstPt);
591 fHistList->Add(fh2NConstLeadingPt);
592 fHistList->Add(fh2PtNch);
593 fHistList->Add(fh2PtNchRan);
594 fHistList->Add(fh2PtNchN);
595 fHistList->Add(fh2PtNchNRan);
596 fHistList->Add(fh2JetPhiEta);
597 fHistList->Add(fh2LeadingJetPhiEta);
598 fHistList->Add(fh2JetEtaPt);
599 fHistList->Add(fh2LeadingJetEtaPt);
600 fHistList->Add(fh2TrackEtaPt);
601 fHistList->Add(fh2LeadingTrackEtaPt);
602 fHistList->Add(fh2JetsLeadingPhiEta );
603 fHistList->Add(fh2JetsLeadingPhiPt);
604 fHistList->Add(fh2TracksLeadingPhiEta);
605 fHistList->Add(fh2TracksLeadingPhiPt);
606 fHistList->Add(fh2TracksLeadingJetPhiPt);
607 fHistList->Add(fh2JetsLeadingPhiPtW);
608 fHistList->Add(fh2TracksLeadingPhiPtW);
609 fHistList->Add(fh2TracksLeadingJetPhiPtW);
610 fHistList->Add(fh2NRecJetsPtRan);
611 fHistList->Add(fh2NConstPtRan);
612 fHistList->Add(fh2NConstLeadingPtRan);
613 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
614 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
617 // =========== Switch on Sumw2 for all histos ===========
618 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
619 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
624 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
627 TH1::AddDirectory(oldStatus);
630 void AliAnalysisTaskJetCluster::Init()
636 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
640 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
643 // handle and reset the output jet branch
645 if(fTCAJetsOut)fTCAJetsOut->Delete();
646 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
647 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
648 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
649 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
651 AliAODJetEventBackground* externalBackground = 0;
652 if(!externalBackground&&fBackgroundBranch.Length()){
653 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
654 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
655 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
658 // Execute analysis for current event
660 AliESDEvent *fESD = 0;
661 if(fUseAODTrackInput){
662 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
664 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
670 // assume that the AOD is in the general output...
673 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
677 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
681 Bool_t selectEvent = false;
682 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
688 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
689 TString vtxTitle(vtxAOD->GetTitle());
690 zVtx = vtxAOD->GetZ();
692 cent = fAOD->GetHeader()->GetCentrality();
693 if(cent<10)cenClass = 0;
694 else if(cent<30)cenClass = 1;
695 else if(cent<50)cenClass = 2;
696 else if(cent<80)cenClass = 3;
697 if(physicsSelection){
698 fh1CentralityPhySel->Fill(cent);
699 fh1ZPhySel->Fill(zVtx);
703 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
704 Float_t yvtx = vtxAOD->GetY();
705 Float_t xvtx = vtxAOD->GetX();
706 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
707 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
708 if(physicsSelection){
714 if(cent<fCentCutLo||cent>fCentCutUp){
725 PostData(1, fHistList);
728 fh1Centrality->Fill(cent);
730 fh1Trials->Fill("#sum{ntrials}",1);
733 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
735 // ==== General variables needed
739 // we simply fetch the tracks/mc particles as a list of AliVParticles
744 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
745 Float_t nCh = recParticles.GetEntries();
747 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
748 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
749 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
753 vector<fastjet::PseudoJet> inputParticlesRec;
754 vector<fastjet::PseudoJet> inputParticlesRecRan;
756 // Generate the random cones
758 AliAODJet vTmpRan(1,0,0,1);
759 for(int i = 0; i < recParticles.GetEntries(); i++){
760 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
761 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
762 // we take total momentum here
763 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
764 jInp.set_user_index(i);
765 inputParticlesRec.push_back(jInp);
767 // the randomized input changes eta and phi, but keeps the p_T
768 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
769 Double_t pT = vp->Pt();
770 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
771 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
773 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
774 Double_t pZ = pT/TMath::Tan(theta);
776 Double_t pX = pT * TMath::Cos(phi);
777 Double_t pY = pT * TMath::Sin(phi);
778 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
779 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
781 jInpRan.set_user_index(i);
782 inputParticlesRecRan.push_back(jInpRan);
783 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
786 // fill the tref array, only needed when we write out jets
789 fRef->Delete(); // make sure to delete before placement new...
790 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
796 if(inputParticlesRec.size()==0){
797 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
798 PostData(1, fHistList);
803 // employ setters for these...
806 // now create the object that holds info about ghosts
807 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
808 // reduce CPU time...
810 fActiveAreaRepeats = 0;
813 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
814 fastjet::AreaType areaType = fastjet::active_area;
815 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
816 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
817 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
819 //range where to compute background
820 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
822 phiMax = 2*TMath::Pi();
823 rapMax = fGhostEtamax - fRparam;
824 rapMin = - fGhostEtamax + fRparam;
825 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
828 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
829 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
832 fh1NJetsRec->Fill(sortedJets.size());
834 // loop over all jets an fill information, first on is the leading jet
836 Int_t nRecOver = inclusiveJets.size();
837 Int_t nRec = inclusiveJets.size();
838 if(inclusiveJets.size()>0){
839 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
840 Double_t area = clustSeq.area(sortedJets[0]);
841 leadingJet.SetEffArea(area,0);
842 Float_t pt = leadingJet.Pt();
843 Int_t nAodOutJets = 0;
844 Int_t nAodOutTracks = 0;
845 AliAODJet *aodOutJet = 0;
848 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
849 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
850 while(pt<ptCut&&iCount<nRec){
854 pt = sortedJets[iCount].perp();
857 if(nRecOver<=0)break;
858 fh2NRecJetsPt->Fill(ptCut,nRecOver);
860 Float_t phi = leadingJet.Phi();
861 if(phi<0)phi+=TMath::Pi()*2.;
862 Float_t eta = leadingJet.Eta();
864 if(externalBackground){
865 // carefull has to be filled in a task before
866 // todo, ReArrange to the botom
867 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
869 pt = leadingJet.Pt() - pTback;
870 // correlation of leading jet with tracks
871 TIterator *recIter = recParticles.MakeIterator();
873 AliVParticle *tmpRecTrack = 0;
874 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
875 Float_t tmpPt = tmpRecTrack->Pt();
877 Float_t tmpPhi = tmpRecTrack->Phi();
878 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
879 Float_t dPhi = phi - tmpPhi;
880 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
881 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
882 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
883 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
885 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
886 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
892 TLorentzVector vecareab;
893 for(int j = 0; j < nRec;j++){
894 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
897 Float_t tmpPt = tmpRec.Pt();
899 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
900 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
901 aodOutJet->GetRefTracks()->Clear();
902 Double_t area1 = clustSeq.area(sortedJets[j]);
903 aodOutJet->SetEffArea(area1,0);
904 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
905 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
906 aodOutJet->SetVectorAreaCharged(&vecareab);
910 Float_t tmpPtBack = 0;
911 if(externalBackground){
912 // carefull has to be filled in a task before
913 // todo, ReArrange to the botom
914 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
916 tmpPt = tmpPt - tmpPtBack;
917 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
919 fh1PtJetsRecIn->Fill(tmpPt);
920 // Fill Spectra with constituentsemacs
921 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
923 fh1NConstRec->Fill(constituents.size());
924 fh2PtNch->Fill(nCh,tmpPt);
925 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
926 fh2NConstPt->Fill(tmpPt,constituents.size());
927 // loop over constiutents and fill spectrum
929 for(unsigned int ic = 0; ic < constituents.size();ic++){
930 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
931 fh1PtJetConstRec->Fill(part->Pt());
933 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
935 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
939 Float_t tmpPhi = tmpRec.Phi();
940 Float_t tmpEta = tmpRec.Eta();
941 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
943 fh1PtJetsLeadingRecIn->Fill(tmpPt);
944 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
945 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
946 fh1NConstLeadingRec->Fill(constituents.size());
947 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
950 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
951 fh2JetEtaPt->Fill(tmpEta,tmpPt);
952 Float_t dPhi = phi - tmpPhi;
953 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
954 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
955 Float_t dEta = eta - tmpRec.Eta();
956 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
957 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
959 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
960 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
962 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
963 }// loop over reconstructed jets
968 // Add the random cones...
969 if(fNRandomCones>0&&fTCARandomConesOut){
970 // create a random jet within the acceptance
971 Double_t etaMax = fTrackEtaWindow - fRparam;
974 Double_t pTC = 1; // small number
975 for(int ir = 0;ir < fNRandomCones;ir++){
976 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
977 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
979 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
980 Double_t pZC = pTC/TMath::Tan(thetaC);
981 Double_t pXC = pTC * TMath::Cos(phiC);
982 Double_t pYC = pTC * TMath::Sin(phiC);
983 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
984 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
986 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
987 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
988 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
993 // test for overlap with previous cones to avoid double counting
994 for(int iic = 0;iic<ir;iic++){
995 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
997 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1004 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1005 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1006 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1007 }// loop over random cones creation
1010 // loop over the reconstructed particles and add up the pT in the random cones
1011 // maybe better to loop over randomized particles not in the real jets...
1012 // but this by definition brings dow average energy in the whole event
1013 AliAODJet vTmpRanR(1,0,0,1);
1014 for(int i = 0; i < recParticles.GetEntries(); i++){
1015 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1016 if(fTCARandomConesOut){
1017 for(int ir = 0;ir < fNRandomCones;ir++){
1018 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1019 if(jC&&jC->DeltaR(vp)<fRparam){
1020 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1023 }// add up energy in cone
1025 // the randomized input changes eta and phi, but keeps the p_T
1026 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1027 Double_t pTR = vp->Pt();
1028 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1029 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1031 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1032 Double_t pZR = pTR/TMath::Tan(thetaR);
1034 Double_t pXR = pTR * TMath::Cos(phiR);
1035 Double_t pYR = pTR * TMath::Sin(phiR);
1036 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1037 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1038 if(fTCARandomConesOutRan){
1039 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1040 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1041 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1042 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1047 }// loop over recparticles
1049 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1050 if(fTCARandomConesOut){
1051 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1052 // rescale the momntum vectors for the random cones
1054 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1056 Double_t etaC = rC->Eta();
1057 Double_t phiC = rC->Phi();
1058 // massless jet, unit vector
1059 pTC = rC->ChargedBgEnergy();
1060 if(pTC<=0)pTC = 0.001; // for almost empty events
1061 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1062 Double_t pZC = pTC/TMath::Tan(thetaC);
1063 Double_t pXC = pTC * TMath::Cos(phiC);
1064 Double_t pYC = pTC * TMath::Sin(phiC);
1065 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1066 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1067 rC->SetBgEnergy(0,0);
1068 rC->SetEffArea(jetArea,0);
1072 if(fTCARandomConesOutRan){
1073 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1074 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1077 Double_t etaC = rC->Eta();
1078 Double_t phiC = rC->Phi();
1079 // massless jet, unit vector
1080 pTC = rC->ChargedBgEnergy();
1081 if(pTC<=0)pTC = 0.001;// for almost empty events
1082 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1083 Double_t pZC = pTC/TMath::Tan(thetaC);
1084 Double_t pXC = pTC * TMath::Cos(phiC);
1085 Double_t pYC = pTC * TMath::Sin(phiC);
1086 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1087 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1088 rC->SetBgEnergy(0,0);
1089 rC->SetEffArea(jetArea,0);
1093 }// if(fNRandomCones
1095 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1101 if(fAODJetBackgroundOut){
1102 vector<fastjet::PseudoJet> jets2=sortedJets;
1103 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1106 Double_t meanarea1=0.;
1109 Double_t meanarea2=0.;
1111 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1112 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1114 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1115 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1117 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1118 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1119 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1120 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1129 // fill track information
1130 Int_t nTrackOver = recParticles.GetSize();
1131 // do the same for tracks and jets
1134 TIterator *recIter = recParticles.MakeIterator();
1135 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1136 Float_t pt = tmpRec->Pt();
1138 // Printf("Leading track p_t %3.3E",pt);
1139 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1140 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1141 while(pt<ptCut&&tmpRec){
1143 tmpRec = (AliVParticle*)(recIter->Next());
1148 if(nTrackOver<=0)break;
1149 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1153 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1154 Float_t phi = leading->Phi();
1155 if(phi<0)phi+=TMath::Pi()*2.;
1156 Float_t eta = leading->Eta();
1158 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1159 Float_t tmpPt = tmpRec->Pt();
1160 Float_t tmpEta = tmpRec->Eta();
1161 fh1PtTracksRecIn->Fill(tmpPt);
1162 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1163 if(tmpRec==leading){
1164 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1165 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1169 Float_t tmpPhi = tmpRec->Phi();
1171 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1172 Float_t dPhi = phi - tmpPhi;
1173 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1174 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1175 Float_t dEta = eta - tmpRec->Eta();
1176 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1177 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1178 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1183 // find the random jets
1185 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1187 // fill the jet information from random track
1188 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1189 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1191 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1193 // loop over all jets an fill information, first on is the leading jet
1195 Int_t nRecOverRan = inclusiveJetsRan.size();
1196 Int_t nRecRan = inclusiveJetsRan.size();
1198 if(inclusiveJetsRan.size()>0){
1199 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1200 Float_t pt = leadingJet.Pt();
1203 TLorentzVector vecarearanb;
1205 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1206 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1207 while(pt<ptCut&&iCount<nRecRan){
1211 pt = sortedJetsRan[iCount].perp();
1214 if(nRecOverRan<=0)break;
1215 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1217 Float_t phi = leadingJet.Phi();
1218 if(phi<0)phi+=TMath::Pi()*2.;
1219 pt = leadingJet.Pt();
1221 // correlation of leading jet with random tracks
1223 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1225 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1227 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1228 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1229 Float_t dPhi = phi - tmpPhi;
1230 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1231 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1232 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1233 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1236 Int_t nAodOutJetsRan = 0;
1237 AliAODJet *aodOutJetRan = 0;
1238 for(int j = 0; j < nRecRan;j++){
1239 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1240 Float_t tmpPt = tmpRec.Pt();
1241 fh1PtJetsRecInRan->Fill(tmpPt);
1242 // Fill Spectra with constituents
1243 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1244 fh1NConstRecRan->Fill(constituents.size());
1245 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1246 fh2PtNchRan->Fill(nCh,tmpPt);
1247 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1250 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1251 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1252 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1253 aodOutJetRan->GetRefTracks()->Clear();
1254 aodOutJetRan->SetEffArea(arearan,0);
1255 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1256 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1257 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1262 Float_t tmpPhi = tmpRec.Phi();
1263 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1266 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1267 fh1NConstLeadingRecRan->Fill(constituents.size());
1268 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1274 if(fAODJetBackgroundOut){
1277 Double_t meanarea3=0.;
1278 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1279 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1280 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1282 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1283 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1292 // do the event selection if activated
1293 if(fJetTriggerPtCut>0){
1294 bool select = false;
1295 Float_t minPt = fJetTriggerPtCut;
1297 // hard coded for now ...
1298 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1299 if(cent<10)minPt = 50;
1300 else if(cent<30)minPt = 42;
1301 else if(cent<50)minPt = 28;
1302 else if(cent<80)minPt = 18;
1305 if(externalBackground)rho = externalBackground->GetBackground(2);
1307 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1308 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1309 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1318 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1319 fh1CentralitySelect->Fill(cent);
1320 fh1ZSelect->Fill(zVtx);
1321 aodH->SetFillAOD(kTRUE);
1325 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1326 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1327 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1328 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1330 PostData(1, fHistList);
1333 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1335 // Terminate analysis
1337 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1341 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1343 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1346 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1347 if(type!=kTrackAODextraonly) {
1348 AliAODEvent *aod = 0;
1349 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1350 else aod = AODEvent();
1352 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1355 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1356 AliAODTrack *tr = aod->GetTrack(it);
1357 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))){
1358 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1361 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1362 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1365 if(tr->Pt()<fTrackPtCut){
1366 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1369 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1374 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1375 AliAODEvent *aod = 0;
1376 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1377 else aod = AODEvent();
1382 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1383 if(!aodExtraTracks)return iCount;
1384 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1385 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1386 if (!track) continue;
1388 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1389 if(!trackAOD)continue;
1390 if((fFilterMask>0)&&!(trackAOD->TestFilterBit(fFilterMask))) continue;
1391 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1392 if(trackAOD->Pt()<fTrackPtCut) continue;
1393 list->Add(trackAOD);
1398 else if (type == kTrackKineAll||type == kTrackKineCharged){
1399 AliMCEvent* mcEvent = MCEvent();
1400 if(!mcEvent)return iCount;
1401 // we want to have alivpartilces so use get track
1402 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1403 if(!mcEvent->IsPhysicalPrimary(it))continue;
1404 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1405 if(type == kTrackKineAll){
1406 if(part->Pt()<fTrackPtCut)continue;
1410 else if(type == kTrackKineCharged){
1411 if(part->Particle()->GetPDG()->Charge()==0)continue;
1412 if(part->Pt()<fTrackPtCut)continue;
1418 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1419 AliAODEvent *aod = 0;
1420 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1421 else aod = AODEvent();
1422 if(!aod)return iCount;
1423 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1424 if(!tca)return iCount;
1425 for(int it = 0;it < tca->GetEntriesFast();++it){
1426 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1427 if(!part->IsPhysicalPrimary())continue;
1428 if(type == kTrackAODMCAll){
1429 if(part->Pt()<fTrackPtCut)continue;
1433 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1434 if(part->Charge()==0)continue;
1435 if(part->Pt()<fTrackPtCut)continue;
1436 if(kTrackAODMCCharged){
1440 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1452 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1453 for(int i = 0; i < particles.GetEntries(); i++){
1454 AliVParticle *vp = (AliVParticle*)particles.At(i);
1455 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1456 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1457 jInp.set_user_index(i);
1458 inputParticles.push_back(jInp);