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),
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),
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 = 200;
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] + 1.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 = 4000;
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 vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
819 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
821 //range where to compute background
822 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
824 phiMax = 2*TMath::Pi();
825 rapMax = fGhostEtamax - fRparam;
826 rapMin = - fGhostEtamax + fRparam;
827 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
831 inclusiveJets = clustSeq.inclusive_jets();
832 sortedJets = sorted_by_pt(inclusiveJets);
834 fh1NJetsRec->Fill(sortedJets.size());
836 // loop over all jets an fill information, first on is the leading jet
838 Int_t nRecOver = inclusiveJets.size();
839 Int_t nRec = inclusiveJets.size();
840 if(inclusiveJets.size()>0){
841 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
842 Double_t area = clustSeq.area(sortedJets[0]);
843 leadingJet.SetEffArea(area,0);
844 Float_t pt = leadingJet.Pt();
845 Int_t nAodOutJets = 0;
846 Int_t nAodOutTracks = 0;
847 AliAODJet *aodOutJet = 0;
850 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
851 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
852 while(pt<ptCut&&iCount<nRec){
856 pt = sortedJets[iCount].perp();
859 if(nRecOver<=0)break;
860 fh2NRecJetsPt->Fill(ptCut,nRecOver);
862 Float_t phi = leadingJet.Phi();
863 if(phi<0)phi+=TMath::Pi()*2.;
864 Float_t eta = leadingJet.Eta();
866 if(externalBackground){
867 // carefull has to be filled in a task before
868 // todo, ReArrange to the botom
869 pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
871 pt = leadingJet.Pt() - pTback;
872 // correlation of leading jet with tracks
873 TIterator *recIter = recParticles.MakeIterator();
875 AliVParticle *tmpRecTrack = 0;
876 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
877 Float_t tmpPt = tmpRecTrack->Pt();
879 Float_t tmpPhi = tmpRecTrack->Phi();
880 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
881 Float_t dPhi = phi - tmpPhi;
882 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
883 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
884 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
885 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
887 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
888 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
894 TLorentzVector vecareab;
895 for(int j = 0; j < nRec;j++){
896 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
899 Float_t tmpPt = tmpRec.Pt();
901 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
902 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
903 aodOutJet->GetRefTracks()->Clear();
904 Double_t area1 = clustSeq.area(sortedJets[j]);
905 aodOutJet->SetEffArea(area1,0);
906 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
907 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
908 aodOutJet->SetVectorAreaCharged(&vecareab);
912 Float_t tmpPtBack = 0;
913 if(externalBackground){
914 // carefull has to be filled in a task before
915 // todo, ReArrange to the botom
916 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
918 tmpPt = tmpPt - tmpPtBack;
919 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
921 fh1PtJetsRecIn->Fill(tmpPt);
922 // Fill Spectra with constituents
923 vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
925 fh1NConstRec->Fill(constituents.size());
926 fh2PtNch->Fill(nCh,tmpPt);
927 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
928 fh2NConstPt->Fill(tmpPt,constituents.size());
929 // loop over constiutents and fill spectrum
931 for(unsigned int ic = 0; ic < constituents.size();ic++){
932 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
933 fh1PtJetConstRec->Fill(part->Pt());
935 aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
937 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
941 Float_t tmpPhi = tmpRec.Phi();
942 Float_t tmpEta = tmpRec.Eta();
943 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
945 fh1PtJetsLeadingRecIn->Fill(tmpPt);
946 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
947 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
948 fh1NConstLeadingRec->Fill(constituents.size());
949 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
952 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
953 fh2JetEtaPt->Fill(tmpEta,tmpPt);
954 Float_t dPhi = phi - tmpPhi;
955 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
956 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
957 Float_t dEta = eta - tmpRec.Eta();
958 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
959 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
961 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
962 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
964 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
965 }// loop over reconstructed jets
970 // Add the random cones...
971 if(fNRandomCones>0&&fTCARandomConesOut){
972 // create a random jet within the acceptance
973 Double_t etaMax = fTrackEtaWindow - fRparam;
976 Double_t pTC = 1; // small number
977 for(int ir = 0;ir < fNRandomCones;ir++){
978 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
979 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
981 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
982 Double_t pZC = pTC/TMath::Tan(thetaC);
983 Double_t pXC = pTC * TMath::Cos(phiC);
984 Double_t pYC = pTC * TMath::Sin(phiC);
985 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
986 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
988 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
989 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
990 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
995 // test for overlap with previous cones to avoid double counting
996 for(int iic = 0;iic<ir;iic++){
997 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
999 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1006 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1007 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1008 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1009 }// loop over random cones creation
1012 // loop over the reconstructed particles and add up the pT in the random cones
1013 // maybe better to loop over randomized particles not in the real jets...
1014 // but this by definition brings dow average energy in the whole event
1015 AliAODJet vTmpRanR(1,0,0,1);
1016 for(int i = 0; i < recParticles.GetEntries(); i++){
1017 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1018 if(fTCARandomConesOut){
1019 for(int ir = 0;ir < fNRandomCones;ir++){
1020 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1021 if(jC&&jC->DeltaR(vp)<fRparam){
1022 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1025 }// add up energy in cone
1027 // the randomized input changes eta and phi, but keeps the p_T
1028 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1029 Double_t pTR = vp->Pt();
1030 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1031 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1033 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1034 Double_t pZR = pTR/TMath::Tan(thetaR);
1036 Double_t pXR = pTR * TMath::Cos(phiR);
1037 Double_t pYR = pTR * TMath::Sin(phiR);
1038 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1039 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1040 if(fTCARandomConesOutRan){
1041 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1042 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1043 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1044 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1049 }// loop over recparticles
1051 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1052 if(fTCARandomConesOut){
1053 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1054 // rescale the momntum vectors for the random cones
1056 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1058 Double_t etaC = rC->Eta();
1059 Double_t phiC = rC->Phi();
1060 // massless jet, unit vector
1061 pTC = rC->ChargedBgEnergy();
1062 if(pTC<=0)pTC = 0.1; // for almost empty events
1063 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1064 Double_t pZC = pTC/TMath::Tan(thetaC);
1065 Double_t pXC = pTC * TMath::Cos(phiC);
1066 Double_t pYC = pTC * TMath::Sin(phiC);
1067 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1068 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1069 rC->SetBgEnergy(0,0);
1070 rC->SetEffArea(jetArea,0);
1074 if(fTCARandomConesOutRan){
1075 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1076 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1079 Double_t etaC = rC->Eta();
1080 Double_t phiC = rC->Phi();
1081 // massless jet, unit vector
1082 pTC = rC->ChargedBgEnergy();
1083 if(pTC<=0)pTC = 0.1;// for almost empty events
1084 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1085 Double_t pZC = pTC/TMath::Tan(thetaC);
1086 Double_t pXC = pTC * TMath::Cos(phiC);
1087 Double_t pYC = pTC * TMath::Sin(phiC);
1088 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1089 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1090 rC->SetBgEnergy(0,0);
1091 rC->SetEffArea(jetArea,0);
1095 }// if(fNRandomCones
1097 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1103 if(fAODJetBackgroundOut){
1104 vector<fastjet::PseudoJet> jets2=sortedJets;
1105 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1108 Double_t meanarea1=0.;
1111 Double_t meanarea2=0.;
1113 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1114 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1116 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1117 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1119 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1120 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1121 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1122 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1131 // fill track information
1132 Int_t nTrackOver = recParticles.GetSize();
1133 // do the same for tracks and jets
1136 TIterator *recIter = recParticles.MakeIterator();
1137 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1138 Float_t pt = tmpRec->Pt();
1140 // Printf("Leading track p_t %3.3E",pt);
1141 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1142 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1143 while(pt<ptCut&&tmpRec){
1145 tmpRec = (AliVParticle*)(recIter->Next());
1150 if(nTrackOver<=0)break;
1151 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1155 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1156 Float_t phi = leading->Phi();
1157 if(phi<0)phi+=TMath::Pi()*2.;
1158 Float_t eta = leading->Eta();
1160 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1161 Float_t tmpPt = tmpRec->Pt();
1162 Float_t tmpEta = tmpRec->Eta();
1163 fh1PtTracksRecIn->Fill(tmpPt);
1164 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1165 if(tmpRec==leading){
1166 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1167 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1171 Float_t tmpPhi = tmpRec->Phi();
1173 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1174 Float_t dPhi = phi - tmpPhi;
1175 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1176 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1177 Float_t dEta = eta - tmpRec->Eta();
1178 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1179 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1180 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1185 // find the random jets
1186 vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
1187 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1189 inclusiveJetsRan = clustSeqRan.inclusive_jets();
1190 sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1192 // fill the jet information from random track
1194 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1196 // loop over all jets an fill information, first on is the leading jet
1198 Int_t nRecOverRan = inclusiveJetsRan.size();
1199 Int_t nRecRan = inclusiveJetsRan.size();
1201 if(inclusiveJetsRan.size()>0){
1202 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1203 Float_t pt = leadingJet.Pt();
1206 TLorentzVector vecarearanb;
1208 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1209 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1210 while(pt<ptCut&&iCount<nRecRan){
1214 pt = sortedJetsRan[iCount].perp();
1217 if(nRecOverRan<=0)break;
1218 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1220 Float_t phi = leadingJet.Phi();
1221 if(phi<0)phi+=TMath::Pi()*2.;
1222 pt = leadingJet.Pt();
1224 // correlation of leading jet with random tracks
1226 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1228 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1230 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1231 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1232 Float_t dPhi = phi - tmpPhi;
1233 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1234 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1235 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1236 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1239 Int_t nAodOutJetsRan = 0;
1240 AliAODJet *aodOutJetRan = 0;
1241 for(int j = 0; j < nRecRan;j++){
1242 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1243 Float_t tmpPt = tmpRec.Pt();
1244 fh1PtJetsRecInRan->Fill(tmpPt);
1245 // Fill Spectra with constituents
1246 vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1247 fh1NConstRecRan->Fill(constituents.size());
1248 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1249 fh2PtNchRan->Fill(nCh,tmpPt);
1250 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1253 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1254 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1255 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1256 aodOutJetRan->GetRefTracks()->Clear();
1257 aodOutJetRan->SetEffArea(arearan,0);
1258 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1259 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1260 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1265 Float_t tmpPhi = tmpRec.Phi();
1266 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1269 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1270 fh1NConstLeadingRecRan->Fill(constituents.size());
1271 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1277 if(fAODJetBackgroundOut){
1280 Double_t meanarea3=0.;
1281 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1282 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1283 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1285 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1286 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1295 // do the event selection if activated
1296 if(fJetTriggerPtCut>0){
1297 bool select = false;
1298 Float_t minPt = fJetTriggerPtCut;
1300 // hard coded for now ...
1301 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1302 if(cent<10)minPt = 50;
1303 else if(cent<30)minPt = 42;
1304 else if(cent<50)minPt = 28;
1305 else if(cent<80)minPt = 18;
1308 if(externalBackground)rho = externalBackground->GetBackground(2);
1310 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1311 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1312 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1321 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1322 fh1CentralitySelect->Fill(cent);
1323 fh1ZSelect->Fill(zVtx);
1324 aodH->SetFillAOD(kTRUE);
1328 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1329 PostData(1, fHistList);
1332 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1334 // Terminate analysis
1336 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1340 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1342 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1345 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1346 if(type!=kTrackAODextraonly) {
1347 AliAODEvent *aod = 0;
1348 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1349 else aod = AODEvent();
1353 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1354 AliAODTrack *tr = aod->GetTrack(it);
1355 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1356 if(TMath::Abs(tr->Eta())>fTrackEtaWindow)continue;
1357 if(tr->Pt()<fTrackPtCut)continue;
1362 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1363 AliAODEvent *aod = 0;
1364 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1365 else aod = AODEvent();
1370 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1371 if(!aodExtraTracks)return iCount;
1372 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1373 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1374 if (!track) continue;
1376 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1377 if(!trackAOD)continue;
1378 if(trackAOD->Pt()<fTrackPtCut) continue;
1379 list->Add(trackAOD);
1384 else if (type == kTrackKineAll||type == kTrackKineCharged){
1385 AliMCEvent* mcEvent = MCEvent();
1386 if(!mcEvent)return iCount;
1387 // we want to have alivpartilces so use get track
1388 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1389 if(!mcEvent->IsPhysicalPrimary(it))continue;
1390 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1391 if(type == kTrackKineAll){
1392 if(part->Pt()<fTrackPtCut)continue;
1396 else if(type == kTrackKineCharged){
1397 if(part->Particle()->GetPDG()->Charge()==0)continue;
1398 if(part->Pt()<fTrackPtCut)continue;
1404 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1405 AliAODEvent *aod = 0;
1406 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1407 else aod = AODEvent();
1408 if(!aod)return iCount;
1409 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1410 if(!tca)return iCount;
1411 for(int it = 0;it < tca->GetEntriesFast();++it){
1412 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1413 if(!part->IsPhysicalPrimary())continue;
1414 if(type == kTrackAODMCAll){
1415 if(part->Pt()<fTrackPtCut)continue;
1419 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1420 if(part->Charge()==0)continue;
1421 if(part->Pt()<fTrackPtCut)continue;
1422 if(kTrackAODMCCharged){
1426 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1438 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1439 for(int i = 0; i < particles.GetEntries(); i++){
1440 AliVParticle *vp = (AliVParticle*)particles.At(i);
1441 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1442 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1443 jInp.set_user_index(i);
1444 inputParticles.push_back(jInp);