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 fUseGlobalSelection(kFALSE),
94 fUseBackgroundCalc(kFALSE),
96 fTrackTypeRec(kTrackUndef),
97 fTrackTypeGen(kTrackUndef),
103 fTrackEtaWindow(0.9),
106 fJetOutputMinPt(0.150),
111 fBackgroundBranch(""),
114 fAlgorithm(fastjet::kt_algorithm),
115 fStrategy(fastjet::Best),
116 fRecombScheme(fastjet::BIpt_scheme),
117 fAreaType(fastjet::active_area),
119 fActiveAreaRepeats(1),
123 fTCARandomConesOut(0x0),
124 fTCARandomConesOutRan(0x0),
125 fAODJetBackgroundOut(0x0),
131 fh1PtHardTrials(0x0),
134 fh1NConstLeadingRec(0x0),
136 fh1PtJetsLeadingRecIn(0x0),
137 fh1PtJetConstRec(0x0),
138 fh1PtJetConstLeadingRec(0x0),
139 fh1PtTracksRecIn(0x0),
140 fh1PtTracksLeadingRecIn(0x0),
142 fh1NConstRecRan(0x0),
143 fh1PtJetsLeadingRecInRan(0x0),
144 fh1NConstLeadingRecRan(0x0),
145 fh1PtJetsRecInRan(0x0),
146 fh1PtTracksGenIn(0x0),
148 fh1CentralityPhySel(0x0),
150 fh1CentralitySelect(0x0),
155 fh2NRecTracksPt(0x0),
157 fh2NConstLeadingPt(0x0),
159 fh2LeadingJetPhiEta(0x0),
161 fh2LeadingJetEtaPt(0x0),
163 fh2LeadingTrackEtaPt(0x0),
164 fh2JetsLeadingPhiEta(0x0),
165 fh2JetsLeadingPhiPt(0x0),
166 fh2TracksLeadingPhiEta(0x0),
167 fh2TracksLeadingPhiPt(0x0),
168 fh2TracksLeadingJetPhiPt(0x0),
169 fh2JetsLeadingPhiPtW(0x0),
170 fh2TracksLeadingPhiPtW(0x0),
171 fh2TracksLeadingJetPhiPtW(0x0),
172 fh2NRecJetsPtRan(0x0),
174 fh2NConstLeadingPtRan(0x0),
179 fh2TracksLeadingJetPhiPtRan(0x0),
180 fh2TracksLeadingJetPhiPtWRan(0x0),
183 for(int i = 0;i<3;i++){
184 fh1BiARandomCones[i] = 0;
185 fh1BiARandomConesRan[i] = 0;
187 for(int i = 0;i<kMaxCent;i++){
188 fh2JetsLeadingPhiPtC[i] = 0;
189 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
190 fh2TracksLeadingJetPhiPtC[i] = 0;
191 fh2TracksLeadingJetPhiPtWC[i] = 0;
195 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
196 AliAnalysisTaskSE(name),
200 fUseAODTrackInput(kFALSE),
201 fUseAODMCInput(kFALSE),
202 fUseGlobalSelection(kFALSE),
203 fUseBackgroundCalc(kFALSE),
205 fTrackTypeRec(kTrackUndef),
206 fTrackTypeGen(kTrackUndef),
208 fNSkipLeadingCone(0),
212 fTrackEtaWindow(0.9),
215 fJetOutputMinPt(0.150),
220 fBackgroundBranch(""),
223 fAlgorithm(fastjet::kt_algorithm),
224 fStrategy(fastjet::Best),
225 fRecombScheme(fastjet::BIpt_scheme),
226 fAreaType(fastjet::active_area),
228 fActiveAreaRepeats(1),
232 fTCARandomConesOut(0x0),
233 fTCARandomConesOutRan(0x0),
234 fAODJetBackgroundOut(0x0),
240 fh1PtHardTrials(0x0),
243 fh1NConstLeadingRec(0x0),
245 fh1PtJetsLeadingRecIn(0x0),
246 fh1PtJetConstRec(0x0),
247 fh1PtJetConstLeadingRec(0x0),
248 fh1PtTracksRecIn(0x0),
249 fh1PtTracksLeadingRecIn(0x0),
251 fh1NConstRecRan(0x0),
252 fh1PtJetsLeadingRecInRan(0x0),
253 fh1NConstLeadingRecRan(0x0),
254 fh1PtJetsRecInRan(0x0),
255 fh1PtTracksGenIn(0x0),
257 fh1CentralityPhySel(0x0),
259 fh1CentralitySelect(0x0),
264 fh2NRecTracksPt(0x0),
266 fh2NConstLeadingPt(0x0),
268 fh2LeadingJetPhiEta(0x0),
270 fh2LeadingJetEtaPt(0x0),
272 fh2LeadingTrackEtaPt(0x0),
273 fh2JetsLeadingPhiEta(0x0),
274 fh2JetsLeadingPhiPt(0x0),
275 fh2TracksLeadingPhiEta(0x0),
276 fh2TracksLeadingPhiPt(0x0),
277 fh2TracksLeadingJetPhiPt(0x0),
278 fh2JetsLeadingPhiPtW(0x0),
279 fh2TracksLeadingPhiPtW(0x0),
280 fh2TracksLeadingJetPhiPtW(0x0),
281 fh2NRecJetsPtRan(0x0),
283 fh2NConstLeadingPtRan(0x0),
288 fh2TracksLeadingJetPhiPtRan(0x0),
289 fh2TracksLeadingJetPhiPtWRan(0x0),
292 for(int i = 0;i<3;i++){
293 fh1BiARandomCones[i] = 0;
294 fh1BiARandomConesRan[i] = 0;
296 for(int i = 0;i<kMaxCent;i++){
297 fh2JetsLeadingPhiPtC[i] = 0;
298 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
299 fh2TracksLeadingJetPhiPtC[i] = 0;
300 fh2TracksLeadingJetPhiPtWC[i] = 0;
302 DefineOutput(1, TList::Class());
307 Bool_t AliAnalysisTaskJetCluster::Notify()
310 // Implemented Notify() to read the cross sections
311 // and number of trials from pyxsec.root
316 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
320 // Create the output container
323 fRandom = new TRandom3(0);
329 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
333 if(fNonStdBranch.Length()!=0)
335 // only create the output branch if we have a name
336 // Create a new branch for jets...
337 // -> cleared in the UserExec....
338 // here we can also have the case that the brnaches are written to a separate file
340 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
341 fTCAJetsOut->SetName(fNonStdBranch.Data());
342 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
346 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
347 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
348 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
350 if(fUseBackgroundCalc){
351 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
352 fAODJetBackgroundOut = new AliAODJetEventBackground();
353 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
354 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
357 // create the branch for the random cones with the same R
358 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
361 if(!AODEvent()->FindListObject(cName.Data())){
362 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
363 fTCARandomConesOut->SetName(cName.Data());
364 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
366 // create the branch with the random for the random cones on the random event
367 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
368 if(!AODEvent()->FindListObject(cName.Data())){
369 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
370 fTCARandomConesOutRan->SetName(cName.Data());
371 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
375 if(fNonStdFile.Length()!=0){
377 // case that we have an AOD extension we need to fetch the jets from the extended output
378 // we identify the extension aod event by looking for the branchname
379 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
380 // case that we have an AOD extension we need can fetch the background maybe from the extended output
381 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
386 if(!fHistList)fHistList = new TList();
387 fHistList->SetOwner();
388 PostData(1, fHistList); // post data in any case once
390 Bool_t oldStatus = TH1::AddDirectoryStatus();
391 TH1::AddDirectory(kFALSE);
396 const Int_t nBinPt = 100;
397 Double_t binLimitsPt[nBinPt+1];
398 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
400 binLimitsPt[iPt] = 0.0;
403 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
407 const Int_t nBinPhi = 90;
408 Double_t binLimitsPhi[nBinPhi+1];
409 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
411 binLimitsPhi[iPhi] = -1.*TMath::Pi();
414 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
420 const Int_t nBinEta = 40;
421 Double_t binLimitsEta[nBinEta+1];
422 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
424 binLimitsEta[iEta] = -2.0;
427 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
431 const int nChMax = 5000;
433 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
434 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
436 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
437 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
440 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
441 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
443 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
444 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
445 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
446 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
449 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
450 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
451 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
453 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
454 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
455 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
456 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
457 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
458 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
459 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
460 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
461 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
462 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
464 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
465 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
466 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
468 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
469 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
470 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
472 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
473 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
474 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
478 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
479 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
480 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
481 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
483 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
484 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);
485 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);
486 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);
490 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
491 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
492 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
493 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
495 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
496 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
497 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
498 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
500 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
501 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
502 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
503 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
507 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
508 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
509 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
510 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
511 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
512 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
513 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
514 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
515 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
516 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
517 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
518 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
520 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
521 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
522 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
523 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
525 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
526 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
527 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
528 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
531 if(fNRandomCones>0&&fUseBackgroundCalc){
532 for(int i = 0;i<3;i++){
533 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
534 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
538 for(int i = 0;i < kMaxCent;i++){
539 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
540 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
541 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
542 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
545 const Int_t saveLevel = 3; // large save level more histos
547 fHistList->Add(fh1Xsec);
548 fHistList->Add(fh1Trials);
550 fHistList->Add(fh1NJetsRec);
551 fHistList->Add(fh1NConstRec);
552 fHistList->Add(fh1NConstLeadingRec);
553 fHistList->Add(fh1PtJetsRecIn);
554 fHistList->Add(fh1PtJetsLeadingRecIn);
555 fHistList->Add(fh1PtTracksRecIn);
556 fHistList->Add(fh1PtTracksLeadingRecIn);
557 fHistList->Add(fh1PtJetConstRec);
558 fHistList->Add(fh1PtJetConstLeadingRec);
559 fHistList->Add(fh1NJetsRecRan);
560 fHistList->Add(fh1NConstRecRan);
561 fHistList->Add(fh1PtJetsLeadingRecInRan);
562 fHistList->Add(fh1NConstLeadingRecRan);
563 fHistList->Add(fh1PtJetsRecInRan);
564 fHistList->Add(fh1Nch);
565 fHistList->Add(fh1Centrality);
566 fHistList->Add(fh1CentralitySelect);
567 fHistList->Add(fh1CentralityPhySel);
568 fHistList->Add(fh1Z);
569 fHistList->Add(fh1ZSelect);
570 fHistList->Add(fh1ZPhySel);
571 if(fNRandomCones>0&&fUseBackgroundCalc){
572 for(int i = 0;i<3;i++){
573 fHistList->Add(fh1BiARandomCones[i]);
574 fHistList->Add(fh1BiARandomConesRan[i]);
577 for(int i = 0;i < kMaxCent;i++){
578 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
579 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
580 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
581 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
584 fHistList->Add(fh2NRecJetsPt);
585 fHistList->Add(fh2NRecTracksPt);
586 fHistList->Add(fh2NConstPt);
587 fHistList->Add(fh2NConstLeadingPt);
588 fHistList->Add(fh2PtNch);
589 fHistList->Add(fh2PtNchRan);
590 fHistList->Add(fh2PtNchN);
591 fHistList->Add(fh2PtNchNRan);
592 fHistList->Add(fh2JetPhiEta);
593 fHistList->Add(fh2LeadingJetPhiEta);
594 fHistList->Add(fh2JetEtaPt);
595 fHistList->Add(fh2LeadingJetEtaPt);
596 fHistList->Add(fh2TrackEtaPt);
597 fHistList->Add(fh2LeadingTrackEtaPt);
598 fHistList->Add(fh2JetsLeadingPhiEta );
599 fHistList->Add(fh2JetsLeadingPhiPt);
600 fHistList->Add(fh2TracksLeadingPhiEta);
601 fHistList->Add(fh2TracksLeadingPhiPt);
602 fHistList->Add(fh2TracksLeadingJetPhiPt);
603 fHistList->Add(fh2JetsLeadingPhiPtW);
604 fHistList->Add(fh2TracksLeadingPhiPtW);
605 fHistList->Add(fh2TracksLeadingJetPhiPtW);
606 fHistList->Add(fh2NRecJetsPtRan);
607 fHistList->Add(fh2NConstPtRan);
608 fHistList->Add(fh2NConstLeadingPtRan);
609 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
610 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
613 // =========== Switch on Sumw2 for all histos ===========
614 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
615 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
620 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
623 TH1::AddDirectory(oldStatus);
626 void AliAnalysisTaskJetCluster::Init()
632 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
636 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
639 if(fUseGlobalSelection){
640 // no selection by the service task, we continue
641 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
642 PostData(1, fHistList);
648 // handle and reset the output jet branch
650 if(fTCAJetsOut)fTCAJetsOut->Delete();
651 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
652 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
653 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
654 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
656 AliAODJetEventBackground* externalBackground = 0;
657 if(!externalBackground&&fBackgroundBranch.Length()){
658 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
659 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
660 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
663 // Execute analysis for current event
665 AliESDEvent *fESD = 0;
666 if(fUseAODTrackInput){
667 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
669 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
675 // assume that the AOD is in the general output...
678 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
682 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
686 Bool_t selectEvent = false;
687 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
693 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
694 TString vtxTitle(vtxAOD->GetTitle());
695 zVtx = vtxAOD->GetZ();
697 cent = fAOD->GetHeader()->GetCentrality();
698 if(cent<10)cenClass = 0;
699 else if(cent<30)cenClass = 1;
700 else if(cent<50)cenClass = 2;
701 else if(cent<80)cenClass = 3;
702 if(physicsSelection){
703 fh1CentralityPhySel->Fill(cent);
704 fh1ZPhySel->Fill(zVtx);
708 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
709 Float_t yvtx = vtxAOD->GetY();
710 Float_t xvtx = vtxAOD->GetX();
711 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
712 if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
713 if(physicsSelection){
719 if(cent<fCentCutLo||cent>fCentCutUp){
726 PostData(1, fHistList);
729 fh1Centrality->Fill(cent);
731 fh1Trials->Fill("#sum{ntrials}",1);
734 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
736 // ==== General variables needed
740 // we simply fetch the tracks/mc particles as a list of AliVParticles
745 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
746 Float_t nCh = recParticles.GetEntries();
748 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
749 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
750 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
754 vector<fastjet::PseudoJet> inputParticlesRec;
755 vector<fastjet::PseudoJet> inputParticlesRecRan;
757 // Generate the random cones
759 AliAODJet vTmpRan(1,0,0,1);
760 for(int i = 0; i < recParticles.GetEntries(); i++){
761 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
762 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
763 // we take total momentum here
764 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
765 jInp.set_user_index(i);
766 inputParticlesRec.push_back(jInp);
768 // the randomized input changes eta and phi, but keeps the p_T
769 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
770 Double_t pT = vp->Pt();
771 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
772 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
774 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
775 Double_t pZ = pT/TMath::Tan(theta);
777 Double_t pX = pT * TMath::Cos(phi);
778 Double_t pY = pT * TMath::Sin(phi);
779 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
780 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
782 jInpRan.set_user_index(i);
783 inputParticlesRecRan.push_back(jInpRan);
784 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
787 // fill the tref array, only needed when we write out jets
790 fRef->Delete(); // make sure to delete before placement new...
791 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
797 if(inputParticlesRec.size()==0){
798 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
799 PostData(1, fHistList);
804 // employ setters for these...
807 // now create the object that holds info about ghosts
808 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
809 // reduce CPU time...
811 fActiveAreaRepeats = 0;
814 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
815 fastjet::AreaType areaType = fastjet::active_area;
816 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
817 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
818 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
820 //range where to compute background
821 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
823 phiMax = 2*TMath::Pi();
824 rapMax = fGhostEtamax - fRparam;
825 rapMin = - fGhostEtamax + fRparam;
826 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
829 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
830 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
833 fh1NJetsRec->Fill(sortedJets.size());
835 // loop over all jets an fill information, first on is the leading jet
837 Int_t nRecOver = inclusiveJets.size();
838 Int_t nRec = inclusiveJets.size();
839 if(inclusiveJets.size()>0){
840 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
841 Double_t area = clustSeq.area(sortedJets[0]);
842 leadingJet.SetEffArea(area,0);
843 Float_t pt = leadingJet.Pt();
844 Int_t nAodOutJets = 0;
845 Int_t nAodOutTracks = 0;
846 AliAODJet *aodOutJet = 0;
849 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
850 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
851 while(pt<ptCut&&iCount<nRec){
855 pt = sortedJets[iCount].perp();
858 if(nRecOver<=0)break;
859 fh2NRecJetsPt->Fill(ptCut,nRecOver);
861 Float_t phi = leadingJet.Phi();
862 if(phi<0)phi+=TMath::Pi()*2.;
863 Float_t eta = leadingJet.Eta();
865 if(externalBackground){
866 // carefull has to be filled in a task before
867 // todo, ReArrange to the botom
868 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
870 pt = leadingJet.Pt() - pTback;
871 // correlation of leading jet with tracks
872 TIterator *recIter = recParticles.MakeIterator();
874 AliVParticle *tmpRecTrack = 0;
875 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
876 Float_t tmpPt = tmpRecTrack->Pt();
878 Float_t tmpPhi = tmpRecTrack->Phi();
879 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
880 Float_t dPhi = phi - tmpPhi;
881 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
882 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
883 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
884 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
886 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
887 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
893 TLorentzVector vecareab;
894 for(int j = 0; j < nRec;j++){
895 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
898 Float_t tmpPt = tmpRec.Pt();
900 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
901 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
902 aodOutJet->GetRefTracks()->Clear();
903 Double_t area1 = clustSeq.area(sortedJets[j]);
904 aodOutJet->SetEffArea(area1,0);
905 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
906 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
907 aodOutJet->SetVectorAreaCharged(&vecareab);
911 Float_t tmpPtBack = 0;
912 if(externalBackground){
913 // carefull has to be filled in a task before
914 // todo, ReArrange to the botom
915 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
917 tmpPt = tmpPt - tmpPtBack;
918 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
920 fh1PtJetsRecIn->Fill(tmpPt);
921 // Fill Spectra with constituents
922 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
924 fh1NConstRec->Fill(constituents.size());
925 fh2PtNch->Fill(nCh,tmpPt);
926 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
927 fh2NConstPt->Fill(tmpPt,constituents.size());
928 // loop over constiutents and fill spectrum
930 for(unsigned int ic = 0; ic < constituents.size();ic++){
931 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
932 fh1PtJetConstRec->Fill(part->Pt());
934 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
936 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
940 Float_t tmpPhi = tmpRec.Phi();
941 Float_t tmpEta = tmpRec.Eta();
942 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
944 fh1PtJetsLeadingRecIn->Fill(tmpPt);
945 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
946 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
947 fh1NConstLeadingRec->Fill(constituents.size());
948 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
951 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
952 fh2JetEtaPt->Fill(tmpEta,tmpPt);
953 Float_t dPhi = phi - tmpPhi;
954 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
955 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
956 Float_t dEta = eta - tmpRec.Eta();
957 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
958 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
960 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
961 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
963 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
964 }// loop over reconstructed jets
969 // Add the random cones...
970 if(fNRandomCones>0&&fTCARandomConesOut){
971 // create a random jet within the acceptance
972 Double_t etaMax = fTrackEtaWindow - fRparam;
975 Double_t pTC = 1; // small number
976 for(int ir = 0;ir < fNRandomCones;ir++){
977 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
978 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
980 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
981 Double_t pZC = pTC/TMath::Tan(thetaC);
982 Double_t pXC = pTC * TMath::Cos(phiC);
983 Double_t pYC = pTC * TMath::Sin(phiC);
984 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
985 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
987 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
988 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
989 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
994 // test for overlap with previous cones to avoid double counting
995 for(int iic = 0;iic<ir;iic++){
996 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
998 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1005 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1006 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1007 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1008 }// loop over random cones creation
1011 // loop over the reconstructed particles and add up the pT in the random cones
1012 // maybe better to loop over randomized particles not in the real jets...
1013 // but this by definition brings dow average energy in the whole event
1014 AliAODJet vTmpRanR(1,0,0,1);
1015 for(int i = 0; i < recParticles.GetEntries(); i++){
1016 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1017 if(fTCARandomConesOut){
1018 for(int ir = 0;ir < fNRandomCones;ir++){
1019 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1020 if(jC&&jC->DeltaR(vp)<fRparam){
1021 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1024 }// add up energy in cone
1026 // the randomized input changes eta and phi, but keeps the p_T
1027 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1028 Double_t pTR = vp->Pt();
1029 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1030 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1032 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1033 Double_t pZR = pTR/TMath::Tan(thetaR);
1035 Double_t pXR = pTR * TMath::Cos(phiR);
1036 Double_t pYR = pTR * TMath::Sin(phiR);
1037 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1038 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1039 if(fTCARandomConesOutRan){
1040 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1041 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1042 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1043 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1048 }// loop over recparticles
1050 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1051 if(fTCARandomConesOut){
1052 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1053 // rescale the momntum vectors for the random cones
1055 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1057 Double_t etaC = rC->Eta();
1058 Double_t phiC = rC->Phi();
1059 // massless jet, unit vector
1060 pTC = rC->ChargedBgEnergy();
1061 if(pTC<=0)pTC = 0.001; // for almost empty events
1062 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1063 Double_t pZC = pTC/TMath::Tan(thetaC);
1064 Double_t pXC = pTC * TMath::Cos(phiC);
1065 Double_t pYC = pTC * TMath::Sin(phiC);
1066 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1067 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1068 rC->SetBgEnergy(0,0);
1069 rC->SetEffArea(jetArea,0);
1073 if(fTCARandomConesOutRan){
1074 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1075 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1078 Double_t etaC = rC->Eta();
1079 Double_t phiC = rC->Phi();
1080 // massless jet, unit vector
1081 pTC = rC->ChargedBgEnergy();
1082 if(pTC<=0)pTC = 0.001;// for almost empty events
1083 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1084 Double_t pZC = pTC/TMath::Tan(thetaC);
1085 Double_t pXC = pTC * TMath::Cos(phiC);
1086 Double_t pYC = pTC * TMath::Sin(phiC);
1087 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1088 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1089 rC->SetBgEnergy(0,0);
1090 rC->SetEffArea(jetArea,0);
1094 }// if(fNRandomCones
1096 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1102 if(fAODJetBackgroundOut){
1103 vector<fastjet::PseudoJet> jets2=sortedJets;
1104 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1107 Double_t meanarea1=0.;
1110 Double_t meanarea2=0.;
1112 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1113 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1115 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1116 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1118 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1119 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1120 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1121 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1130 // fill track information
1131 Int_t nTrackOver = recParticles.GetSize();
1132 // do the same for tracks and jets
1135 TIterator *recIter = recParticles.MakeIterator();
1136 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1137 Float_t pt = tmpRec->Pt();
1139 // Printf("Leading track p_t %3.3E",pt);
1140 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1141 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1142 while(pt<ptCut&&tmpRec){
1144 tmpRec = (AliVParticle*)(recIter->Next());
1149 if(nTrackOver<=0)break;
1150 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1154 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1155 Float_t phi = leading->Phi();
1156 if(phi<0)phi+=TMath::Pi()*2.;
1157 Float_t eta = leading->Eta();
1159 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1160 Float_t tmpPt = tmpRec->Pt();
1161 Float_t tmpEta = tmpRec->Eta();
1162 fh1PtTracksRecIn->Fill(tmpPt);
1163 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1164 if(tmpRec==leading){
1165 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1166 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1170 Float_t tmpPhi = tmpRec->Phi();
1172 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1173 Float_t dPhi = phi - tmpPhi;
1174 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1175 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1176 Float_t dEta = eta - tmpRec->Eta();
1177 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1178 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1179 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1184 // find the random jets
1186 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1188 // fill the jet information from random track
1189 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1190 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1192 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1194 // loop over all jets an fill information, first on is the leading jet
1196 Int_t nRecOverRan = inclusiveJetsRan.size();
1197 Int_t nRecRan = inclusiveJetsRan.size();
1199 if(inclusiveJetsRan.size()>0){
1200 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1201 Float_t pt = leadingJet.Pt();
1204 TLorentzVector vecarearanb;
1206 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1207 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1208 while(pt<ptCut&&iCount<nRecRan){
1212 pt = sortedJetsRan[iCount].perp();
1215 if(nRecOverRan<=0)break;
1216 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1218 Float_t phi = leadingJet.Phi();
1219 if(phi<0)phi+=TMath::Pi()*2.;
1220 pt = leadingJet.Pt();
1222 // correlation of leading jet with random tracks
1224 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1226 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1228 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1229 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1230 Float_t dPhi = phi - tmpPhi;
1231 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1232 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1233 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1234 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1237 Int_t nAodOutJetsRan = 0;
1238 AliAODJet *aodOutJetRan = 0;
1239 for(int j = 0; j < nRecRan;j++){
1240 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1241 Float_t tmpPt = tmpRec.Pt();
1242 fh1PtJetsRecInRan->Fill(tmpPt);
1243 // Fill Spectra with constituents
1244 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1245 fh1NConstRecRan->Fill(constituents.size());
1246 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1247 fh2PtNchRan->Fill(nCh,tmpPt);
1248 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1251 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1252 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1253 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1254 aodOutJetRan->GetRefTracks()->Clear();
1255 aodOutJetRan->SetEffArea(arearan,0);
1256 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1257 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1258 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1263 Float_t tmpPhi = tmpRec.Phi();
1264 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1267 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1268 fh1NConstLeadingRecRan->Fill(constituents.size());
1269 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1275 if(fAODJetBackgroundOut){
1278 Double_t meanarea3=0.;
1279 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1280 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1281 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1283 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1284 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1293 // do the event selection if activated
1294 if(fJetTriggerPtCut>0){
1295 bool select = false;
1296 Float_t minPt = fJetTriggerPtCut;
1298 // hard coded for now ...
1299 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1300 if(cent<10)minPt = 50;
1301 else if(cent<30)minPt = 42;
1302 else if(cent<50)minPt = 28;
1303 else if(cent<80)minPt = 18;
1306 if(externalBackground)rho = externalBackground->GetBackground(2);
1308 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1309 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1310 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1319 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1320 fh1CentralitySelect->Fill(cent);
1321 fh1ZSelect->Fill(zVtx);
1322 aodH->SetFillAOD(kTRUE);
1325 if (fDebug > 10)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1326 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1327 PostData(1, fHistList);
1330 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1332 // Terminate analysis
1334 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1338 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1340 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1343 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1344 if(type!=kTrackAODextraonly) {
1345 AliAODEvent *aod = 0;
1346 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1347 else aod = AODEvent();
1349 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1352 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1353 AliAODTrack *tr = aod->GetTrack(it);
1354 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))){
1355 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1358 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1359 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1362 if(tr->Pt()<fTrackPtCut){
1363 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1366 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1371 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1372 AliAODEvent *aod = 0;
1373 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1374 else aod = AODEvent();
1379 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1380 if(!aodExtraTracks)return iCount;
1381 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1382 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1383 if (!track) continue;
1385 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1386 if(!trackAOD)continue;
1387 if(trackAOD->Pt()<fTrackPtCut) continue;
1388 list->Add(trackAOD);
1393 else if (type == kTrackKineAll||type == kTrackKineCharged){
1394 AliMCEvent* mcEvent = MCEvent();
1395 if(!mcEvent)return iCount;
1396 // we want to have alivpartilces so use get track
1397 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1398 if(!mcEvent->IsPhysicalPrimary(it))continue;
1399 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1400 if(type == kTrackKineAll){
1401 if(part->Pt()<fTrackPtCut)continue;
1405 else if(type == kTrackKineCharged){
1406 if(part->Particle()->GetPDG()->Charge()==0)continue;
1407 if(part->Pt()<fTrackPtCut)continue;
1413 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1414 AliAODEvent *aod = 0;
1415 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1416 else aod = AODEvent();
1417 if(!aod)return iCount;
1418 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1419 if(!tca)return iCount;
1420 for(int it = 0;it < tca->GetEntriesFast();++it){
1421 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1422 if(!part->IsPhysicalPrimary())continue;
1423 if(type == kTrackAODMCAll){
1424 if(part->Pt()<fTrackPtCut)continue;
1428 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1429 if(part->Charge()==0)continue;
1430 if(part->Pt()<fTrackPtCut)continue;
1431 if(kTrackAODMCCharged){
1435 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1447 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1448 for(int i = 0; i < particles.GetEntries(); i++){
1449 AliVParticle *vp = (AliVParticle*)particles.At(i);
1450 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1451 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1452 jInp.set_user_index(i);
1453 inputParticles.push_back(jInp);