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>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include "TDatabasePDG.h"
41 #include "AliAnalysisTaskJetCluster.h"
42 #include "AliAnalysisManager.h"
43 #include "AliJetFinder.h"
44 #include "AliJetHeader.h"
45 #include "AliJetReader.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliAODHandler.h"
49 #include "AliAODExtension.h"
50 #include "AliAODTrack.h"
51 #include "AliAODJet.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliJetKineReaderHeader.h"
58 #include "AliGenCocktailEventHeader.h"
59 #include "AliInputEventHandler.h"
60 #include "AliAODJetEventBackground.h"
62 #include "fastjet/PseudoJet.hh"
63 #include "fastjet/ClusterSequenceArea.hh"
64 #include "fastjet/AreaDefinition.hh"
65 #include "fastjet/JetDefinition.hh"
66 // get info on how fastjet was configured
67 #include "fastjet/config.h"
70 ClassImp(AliAnalysisTaskJetCluster)
72 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
80 if(fTCAJetsOut)fTCAJetsOut->Delete();
83 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
84 delete fTCAJetsOutRan;
86 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
87 delete fTCARandomConesOut;
89 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
90 delete fTCARandomConesOutRan;
92 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
93 delete fAODJetBackgroundOut;
96 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
101 fUseAODTrackInput(kFALSE),
102 fUseAODMCInput(kFALSE),
103 fUseBackgroundCalc(kFALSE),
104 fEventSelection(kFALSE),
108 fTrackTypeRec(kTrackUndef),
109 fTrackTypeGen(kTrackUndef),
111 fNSkipLeadingCone(0),
115 fTrackEtaWindow(0.9),
118 fJetOutputMinPt(0.150),
119 fMaxTrackPtInJet(100.),
126 fBackgroundBranch(""),
137 fUseTrMomentumSmearing(kFALSE),
138 fUseDiceEfficiency(kFALSE),
140 fAlgorithm(fastjet::kt_algorithm),
141 fStrategy(fastjet::Best),
142 fRecombScheme(fastjet::BIpt_scheme),
143 fAreaType(fastjet::active_area),
145 fActiveAreaRepeats(1),
149 fTCARandomConesOut(0x0),
150 fTCARandomConesOutRan(0x0),
151 fAODJetBackgroundOut(0x0),
157 fh1PtHardTrials(0x0),
160 fh1NConstLeadingRec(0x0),
162 fh1PtJetsLeadingRecIn(0x0),
163 fh1PtJetConstRec(0x0),
164 fh1PtJetConstLeadingRec(0x0),
165 fh1PtTracksRecIn(0x0),
166 fh1PtTracksLeadingRecIn(0x0),
168 fh1NConstRecRan(0x0),
169 fh1PtJetsLeadingRecInRan(0x0),
170 fh1NConstLeadingRecRan(0x0),
171 fh1PtJetsRecInRan(0x0),
172 fh1PtTracksGenIn(0x0),
174 fh1CentralityPhySel(0x0),
176 fh1CentralitySelect(0x0),
181 fh2NRecTracksPt(0x0),
183 fh2NConstLeadingPt(0x0),
185 fh2LeadingJetPhiEta(0x0),
187 fh2LeadingJetEtaPt(0x0),
189 fh2LeadingTrackEtaPt(0x0),
190 fh2JetsLeadingPhiEta(0x0),
191 fh2JetsLeadingPhiPt(0x0),
192 fh2TracksLeadingPhiEta(0x0),
193 fh2TracksLeadingPhiPt(0x0),
194 fh2TracksLeadingJetPhiPt(0x0),
195 fh2JetsLeadingPhiPtW(0x0),
196 fh2TracksLeadingPhiPtW(0x0),
197 fh2TracksLeadingJetPhiPtW(0x0),
198 fh2NRecJetsPtRan(0x0),
200 fh2NConstLeadingPtRan(0x0),
205 fh2TracksLeadingJetPhiPtRan(0x0),
206 fh2TracksLeadingJetPhiPtWRan(0x0),
207 fh2PtGenPtSmeared(0x0),
209 fp1PtResolution(0x0),
216 for(int i = 0;i<3;i++){
217 fh1BiARandomCones[i] = 0;
218 fh1BiARandomConesRan[i] = 0;
220 for(int i = 0;i<kMaxCent;i++){
221 fh2JetsLeadingPhiPtC[i] = 0;
222 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
223 fh2TracksLeadingJetPhiPtC[i] = 0;
224 fh2TracksLeadingJetPhiPtWC[i] = 0;
228 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
229 AliAnalysisTaskSE(name),
233 fUseAODTrackInput(kFALSE),
234 fUseAODMCInput(kFALSE),
235 fUseBackgroundCalc(kFALSE),
236 fEventSelection(kFALSE),
240 fTrackTypeRec(kTrackUndef),
241 fTrackTypeGen(kTrackUndef),
243 fNSkipLeadingCone(0),
247 fTrackEtaWindow(0.9),
250 fJetOutputMinPt(0.150),
251 fMaxTrackPtInJet(100.),
258 fBackgroundBranch(""),
269 fUseTrMomentumSmearing(kFALSE),
270 fUseDiceEfficiency(kFALSE),
272 fAlgorithm(fastjet::kt_algorithm),
273 fStrategy(fastjet::Best),
274 fRecombScheme(fastjet::BIpt_scheme),
275 fAreaType(fastjet::active_area),
277 fActiveAreaRepeats(1),
281 fTCARandomConesOut(0x0),
282 fTCARandomConesOutRan(0x0),
283 fAODJetBackgroundOut(0x0),
289 fh1PtHardTrials(0x0),
292 fh1NConstLeadingRec(0x0),
294 fh1PtJetsLeadingRecIn(0x0),
295 fh1PtJetConstRec(0x0),
296 fh1PtJetConstLeadingRec(0x0),
297 fh1PtTracksRecIn(0x0),
298 fh1PtTracksLeadingRecIn(0x0),
300 fh1NConstRecRan(0x0),
301 fh1PtJetsLeadingRecInRan(0x0),
302 fh1NConstLeadingRecRan(0x0),
303 fh1PtJetsRecInRan(0x0),
304 fh1PtTracksGenIn(0x0),
306 fh1CentralityPhySel(0x0),
308 fh1CentralitySelect(0x0),
313 fh2NRecTracksPt(0x0),
315 fh2NConstLeadingPt(0x0),
317 fh2LeadingJetPhiEta(0x0),
319 fh2LeadingJetEtaPt(0x0),
321 fh2LeadingTrackEtaPt(0x0),
322 fh2JetsLeadingPhiEta(0x0),
323 fh2JetsLeadingPhiPt(0x0),
324 fh2TracksLeadingPhiEta(0x0),
325 fh2TracksLeadingPhiPt(0x0),
326 fh2TracksLeadingJetPhiPt(0x0),
327 fh2JetsLeadingPhiPtW(0x0),
328 fh2TracksLeadingPhiPtW(0x0),
329 fh2TracksLeadingJetPhiPtW(0x0),
330 fh2NRecJetsPtRan(0x0),
332 fh2NConstLeadingPtRan(0x0),
337 fh2TracksLeadingJetPhiPtRan(0x0),
338 fh2TracksLeadingJetPhiPtWRan(0x0),
339 fh2PtGenPtSmeared(0x0),
341 fp1PtResolution(0x0),
348 for(int i = 0;i<3;i++){
349 fh1BiARandomCones[i] = 0;
350 fh1BiARandomConesRan[i] = 0;
352 for(int i = 0;i<kMaxCent;i++){
353 fh2JetsLeadingPhiPtC[i] = 0;
354 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
355 fh2TracksLeadingJetPhiPtC[i] = 0;
356 fh2TracksLeadingJetPhiPtWC[i] = 0;
358 DefineOutput(1, TList::Class());
363 Bool_t AliAnalysisTaskJetCluster::Notify()
366 // Implemented Notify() to read the cross sections
367 // and number of trials from pyxsec.root
372 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
376 // Create the output container
379 fRandom = new TRandom3(0);
385 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
389 if(fNonStdBranch.Length()!=0)
391 // only create the output branch if we have a name
392 // Create a new branch for jets...
393 // -> cleared in the UserExec....
394 // here we can also have the case that the brnaches are written to a separate file
397 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
398 fTCAJetsOut->SetName(fNonStdBranch.Data());
399 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
402 if(fJetTypes&kJetRan){
403 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
404 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
405 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
406 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%d",fNonStdBranch.Data(),"random",fUseTrMomentumSmearing,fUseDiceEfficiency));
407 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
410 if(fUseBackgroundCalc){
411 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
412 fAODJetBackgroundOut = new AliAODJetEventBackground();
413 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
414 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
415 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency));
417 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
420 // create the branch for the random cones with the same R
421 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
422 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
423 cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
427 if(!AODEvent()->FindListObject(cName.Data())){
428 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
429 fTCARandomConesOut->SetName(cName.Data());
430 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
433 // create the branch with the random for the random cones on the random event
434 if(fJetTypes&kRCRan){
435 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
436 if(!AODEvent()->FindListObject(cName.Data())){
437 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
438 fTCARandomConesOutRan->SetName(cName.Data());
439 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
444 if(fNonStdFile.Length()!=0){
446 // case that we have an AOD extension we need to fetch the jets from the extended output
447 // we identify the extension aod event by looking for the branchname
448 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
449 // case that we have an AOD extension we need can fetch the background maybe from the extended output
450 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
454 // FitMomentumResolution();
457 if(!fHistList)fHistList = new TList();
458 fHistList->SetOwner();
459 PostData(1, fHistList); // post data in any case once
461 Bool_t oldStatus = TH1::AddDirectoryStatus();
462 TH1::AddDirectory(kFALSE);
467 const Int_t nBinPt = 100;
468 Double_t binLimitsPt[nBinPt+1];
469 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
471 binLimitsPt[iPt] = 0.0;
474 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
478 const Int_t nBinPhi = 90;
479 Double_t binLimitsPhi[nBinPhi+1];
480 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
482 binLimitsPhi[iPhi] = -1.*TMath::Pi();
485 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
491 const Int_t nBinEta = 40;
492 Double_t binLimitsEta[nBinEta+1];
493 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
495 binLimitsEta[iEta] = -2.0;
498 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
502 const int nChMax = 5000;
504 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
505 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
507 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
508 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
511 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
512 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
514 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
515 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
516 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
517 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
520 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
521 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
522 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
524 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
525 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
526 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
527 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
528 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
529 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
530 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
531 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
532 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
533 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
535 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
536 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
537 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
539 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
540 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
541 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
543 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
544 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
545 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
549 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
550 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
551 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
552 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
554 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
555 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);
556 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);
557 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);
561 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
562 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
563 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
564 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
566 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
567 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
568 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
569 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
571 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
572 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
573 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
574 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
578 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
579 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
580 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
581 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
582 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
583 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
584 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
585 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
586 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
587 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
588 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
589 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
591 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
592 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
593 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
594 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
596 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
597 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
598 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
599 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
601 //Detector level effects histos
602 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
604 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
605 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
607 fHistList->Add(fh2PtGenPtSmeared);
608 fHistList->Add(fp1Efficiency);
609 fHistList->Add(fp1PtResolution);
611 if(fNRandomCones>0&&fUseBackgroundCalc){
612 for(int i = 0;i<3;i++){
613 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
614 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
618 for(int i = 0;i < kMaxCent;i++){
619 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
620 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
621 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
622 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
625 const Int_t saveLevel = 3; // large save level more histos
627 fHistList->Add(fh1Xsec);
628 fHistList->Add(fh1Trials);
630 fHistList->Add(fh1NJetsRec);
631 fHistList->Add(fh1NConstRec);
632 fHistList->Add(fh1NConstLeadingRec);
633 fHistList->Add(fh1PtJetsRecIn);
634 fHistList->Add(fh1PtJetsLeadingRecIn);
635 fHistList->Add(fh1PtTracksRecIn);
636 fHistList->Add(fh1PtTracksLeadingRecIn);
637 fHistList->Add(fh1PtJetConstRec);
638 fHistList->Add(fh1PtJetConstLeadingRec);
639 fHistList->Add(fh1NJetsRecRan);
640 fHistList->Add(fh1NConstRecRan);
641 fHistList->Add(fh1PtJetsLeadingRecInRan);
642 fHistList->Add(fh1NConstLeadingRecRan);
643 fHistList->Add(fh1PtJetsRecInRan);
644 fHistList->Add(fh1Nch);
645 fHistList->Add(fh1Centrality);
646 fHistList->Add(fh1CentralitySelect);
647 fHistList->Add(fh1CentralityPhySel);
648 fHistList->Add(fh1Z);
649 fHistList->Add(fh1ZSelect);
650 fHistList->Add(fh1ZPhySel);
651 if(fNRandomCones>0&&fUseBackgroundCalc){
652 for(int i = 0;i<3;i++){
653 fHistList->Add(fh1BiARandomCones[i]);
654 fHistList->Add(fh1BiARandomConesRan[i]);
657 for(int i = 0;i < kMaxCent;i++){
658 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
659 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
660 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
661 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
664 fHistList->Add(fh2NRecJetsPt);
665 fHistList->Add(fh2NRecTracksPt);
666 fHistList->Add(fh2NConstPt);
667 fHistList->Add(fh2NConstLeadingPt);
668 fHistList->Add(fh2PtNch);
669 fHistList->Add(fh2PtNchRan);
670 fHistList->Add(fh2PtNchN);
671 fHistList->Add(fh2PtNchNRan);
672 fHistList->Add(fh2JetPhiEta);
673 fHistList->Add(fh2LeadingJetPhiEta);
674 fHistList->Add(fh2JetEtaPt);
675 fHistList->Add(fh2LeadingJetEtaPt);
676 fHistList->Add(fh2TrackEtaPt);
677 fHistList->Add(fh2LeadingTrackEtaPt);
678 fHistList->Add(fh2JetsLeadingPhiEta );
679 fHistList->Add(fh2JetsLeadingPhiPt);
680 fHistList->Add(fh2TracksLeadingPhiEta);
681 fHistList->Add(fh2TracksLeadingPhiPt);
682 fHistList->Add(fh2TracksLeadingJetPhiPt);
683 fHistList->Add(fh2JetsLeadingPhiPtW);
684 fHistList->Add(fh2TracksLeadingPhiPtW);
685 fHistList->Add(fh2TracksLeadingJetPhiPtW);
686 fHistList->Add(fh2NRecJetsPtRan);
687 fHistList->Add(fh2NConstPtRan);
688 fHistList->Add(fh2NConstLeadingPtRan);
689 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
690 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
693 // =========== Switch on Sumw2 for all histos ===========
694 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
695 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
700 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
703 TH1::AddDirectory(oldStatus);
706 void AliAnalysisTaskJetCluster::Init()
712 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
714 FitMomentumResolution();
718 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
721 // handle and reset the output jet branch
723 if(fTCAJetsOut)fTCAJetsOut->Delete();
724 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
725 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
726 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
727 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
729 AliAODJetEventBackground* externalBackground = 0;
730 if(!externalBackground&&fBackgroundBranch.Length()){
731 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
732 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
733 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
736 // Execute analysis for current event
738 AliESDEvent *fESD = 0;
739 if(fUseAODTrackInput){
740 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
742 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
748 // assume that the AOD is in the general output...
751 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
755 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
759 //Check if information is provided detector level effects
760 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
761 if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
763 Bool_t selectEvent = false;
764 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
770 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
771 TString vtxTitle(vtxAOD->GetTitle());
772 zVtx = vtxAOD->GetZ();
774 cent = fAOD->GetHeader()->GetCentrality();
775 if(cent<10)cenClass = 0;
776 else if(cent<30)cenClass = 1;
777 else if(cent<50)cenClass = 2;
778 else if(cent<80)cenClass = 3;
779 if(physicsSelection){
780 fh1CentralityPhySel->Fill(cent);
781 fh1ZPhySel->Fill(zVtx);
785 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
786 Float_t yvtx = vtxAOD->GetY();
787 Float_t xvtx = vtxAOD->GetX();
788 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
789 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
790 if(physicsSelection){
796 if(cent<fCentCutLo||cent>fCentCutUp){
807 PostData(1, fHistList);
810 fh1Centrality->Fill(cent);
812 fh1Trials->Fill("#sum{ntrials}",1);
815 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
817 // ==== General variables needed
821 // we simply fetch the tracks/mc particles as a list of AliVParticles
826 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
827 Float_t nCh = recParticles.GetEntries();
829 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
830 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
831 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
835 vector<fastjet::PseudoJet> inputParticlesRec;
836 vector<fastjet::PseudoJet> inputParticlesRecRan;
838 // Generate the random cones
840 AliAODJet vTmpRan(1,0,0,1);
841 for(int i = 0; i < recParticles.GetEntries(); i++){
842 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
844 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
845 // we take total momentum here
847 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
848 //Add particles to fastjet in case we are not running toy model
849 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
850 jInp.set_user_index(i);
851 inputParticlesRec.push_back(jInp);
853 else if(fUseDiceEfficiency) {
855 // Dice to decide if particle is kept or not - toy model for efficiency
857 Double_t rnd = fRandom->Uniform(1.);
858 Double_t pT = vp->Pt();
859 Double_t eff[3] = {0.};
861 if(pT>10.) pTtmp = 10.;
862 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
863 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
864 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
866 //Sort efficiencies from large to small
867 if(eff1>eff2 && eff1>eff3) {
882 else if(eff2>eff1 && eff2>eff3) {
897 else if(eff3>eff1 && eff3>eff2) {
913 Double_t sumEff = eff[0]+eff[1]+eff[2];
914 fp1Efficiency->Fill(vp->Pt(),sumEff);
915 if(rnd>sumEff) continue;
917 if(fUseTrMomentumSmearing) {
918 //Smear momentum of generated particle
920 //Select hybrid track category
922 smear = GetMomentumSmearing(cat[2],pT);
923 else if(rnd<=(eff[2]+eff[1]))
924 smear = GetMomentumSmearing(cat[1],pT);
926 smear = GetMomentumSmearing(cat[0],pT);
928 fp1PtResolution->Fill(vp->Pt(),smear);
930 Double_t sigma = vp->Pt()*smear;
931 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
933 Double_t phi = vp->Phi();
934 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
935 Double_t pX = pTrec * TMath::Cos(phi);
936 Double_t pY = pTrec * TMath::Sin(phi);
937 Double_t pZ = pTrec/TMath::Tan(theta);
938 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
940 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
942 fastjet::PseudoJet jInp(pX,pY,pZ,p);
943 jInp.set_user_index(i);
944 inputParticlesRec.push_back(jInp);
948 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
949 jInp.set_user_index(i);
950 inputParticlesRec.push_back(jInp);
956 // the randomized input changes eta and phi, but keeps the p_T
957 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
958 Double_t pT = vp->Pt();
959 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
960 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
962 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
963 Double_t pZ = pT/TMath::Tan(theta);
965 Double_t pX = pT * TMath::Cos(phi);
966 Double_t pY = pT * TMath::Sin(phi);
967 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
968 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
970 jInpRan.set_user_index(i);
971 inputParticlesRecRan.push_back(jInpRan);
972 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
975 // fill the tref array, only needed when we write out jets
978 fRef->Delete(); // make sure to delete before placement new...
979 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
980 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
983 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
987 if(inputParticlesRec.size()==0){
988 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
989 PostData(1, fHistList);
994 // employ setters for these...
997 // now create the object that holds info about ghosts
999 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1000 // reduce CPU time...
1002 fActiveAreaRepeats = 0;
1006 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1007 fastjet::AreaType areaType = fastjet::active_area;
1008 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1009 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1010 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1012 //range where to compute background
1013 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1015 phiMax = 2*TMath::Pi();
1016 rapMax = fGhostEtamax - fRparam;
1017 rapMin = - fGhostEtamax + fRparam;
1018 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1021 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1022 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1025 fh1NJetsRec->Fill(sortedJets.size());
1027 // loop over all jets an fill information, first on is the leading jet
1029 Int_t nRecOver = inclusiveJets.size();
1030 Int_t nRec = inclusiveJets.size();
1031 if(inclusiveJets.size()>0){
1032 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1033 Double_t area = clustSeq.area(sortedJets[0]);
1034 leadingJet.SetEffArea(area,0);
1035 Float_t pt = leadingJet.Pt();
1036 Int_t nAodOutJets = 0;
1037 Int_t nAodOutTracks = 0;
1038 AliAODJet *aodOutJet = 0;
1041 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1042 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1043 while(pt<ptCut&&iCount<nRec){
1047 pt = sortedJets[iCount].perp();
1050 if(nRecOver<=0)break;
1051 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1053 Float_t phi = leadingJet.Phi();
1054 if(phi<0)phi+=TMath::Pi()*2.;
1055 Float_t eta = leadingJet.Eta();
1057 if(externalBackground){
1058 // carefull has to be filled in a task before
1059 // todo, ReArrange to the botom
1060 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1062 pt = leadingJet.Pt() - pTback;
1063 // correlation of leading jet with tracks
1064 TIterator *recIter = recParticles.MakeIterator();
1066 AliVParticle *tmpRecTrack = 0;
1067 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1068 Float_t tmpPt = tmpRecTrack->Pt();
1070 Float_t tmpPhi = tmpRecTrack->Phi();
1071 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1072 Float_t dPhi = phi - tmpPhi;
1073 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1074 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1075 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1076 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1078 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1079 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1085 TLorentzVector vecareab;
1086 for(int j = 0; j < nRec;j++){
1087 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1090 Float_t tmpPt = tmpRec.Pt();
1092 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1093 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1094 aodOutJet->GetRefTracks()->Clear();
1095 Double_t area1 = clustSeq.area(sortedJets[j]);
1096 aodOutJet->SetEffArea(area1,0);
1097 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1098 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1099 aodOutJet->SetVectorAreaCharged(&vecareab);
1103 Float_t tmpPtBack = 0;
1104 if(externalBackground){
1105 // carefull has to be filled in a task before
1106 // todo, ReArrange to the botom
1107 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1109 tmpPt = tmpPt - tmpPtBack;
1110 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1112 fh1PtJetsRecIn->Fill(tmpPt);
1113 // Fill Spectra with constituentsemacs
1114 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1116 fh1NConstRec->Fill(constituents.size());
1117 fh2PtNch->Fill(nCh,tmpPt);
1118 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1119 fh2NConstPt->Fill(tmpPt,constituents.size());
1120 // loop over constiutents and fill spectrum
1122 for(unsigned int ic = 0; ic < constituents.size();ic++){
1123 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1126 fh1PtJetConstRec->Fill(part->Pt());
1128 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1129 if(part->Pt()>fMaxTrackPtInJet)aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1131 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1135 Float_t tmpPhi = tmpRec.Phi();
1136 Float_t tmpEta = tmpRec.Eta();
1137 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1139 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1140 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1141 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1142 fh1NConstLeadingRec->Fill(constituents.size());
1143 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1146 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1147 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1148 Float_t dPhi = phi - tmpPhi;
1149 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1150 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1151 Float_t dEta = eta - tmpRec.Eta();
1152 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1153 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1155 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1156 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1158 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1159 }// loop over reconstructed jets
1164 // Add the random cones...
1165 if(fNRandomCones>0&&fTCARandomConesOut){
1166 // create a random jet within the acceptance
1167 Double_t etaMax = fTrackEtaWindow - fRparam;
1170 Double_t pTC = 1; // small number
1171 for(int ir = 0;ir < fNRandomCones;ir++){
1172 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1173 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1175 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1176 Double_t pZC = pTC/TMath::Tan(thetaC);
1177 Double_t pXC = pTC * TMath::Cos(phiC);
1178 Double_t pYC = pTC * TMath::Sin(phiC);
1179 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1180 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1182 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1183 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1184 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1189 // test for overlap with previous cones to avoid double counting
1190 for(int iic = 0;iic<ir;iic++){
1191 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1193 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1200 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1201 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1202 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1203 }// loop over random cones creation
1206 // loop over the reconstructed particles and add up the pT in the random cones
1207 // maybe better to loop over randomized particles not in the real jets...
1208 // but this by definition brings dow average energy in the whole event
1209 AliAODJet vTmpRanR(1,0,0,1);
1210 for(int i = 0; i < recParticles.GetEntries(); i++){
1211 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1212 if(fTCARandomConesOut){
1213 for(int ir = 0;ir < fNRandomCones;ir++){
1214 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1215 if(jC&&jC->DeltaR(vp)<fRparam){
1216 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1217 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1220 }// add up energy in cone
1222 // the randomized input changes eta and phi, but keeps the p_T
1223 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1224 Double_t pTR = vp->Pt();
1225 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1226 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1228 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1229 Double_t pZR = pTR/TMath::Tan(thetaR);
1231 Double_t pXR = pTR * TMath::Cos(phiR);
1232 Double_t pYR = pTR * TMath::Sin(phiR);
1233 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1234 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1235 if(fTCARandomConesOutRan){
1236 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1237 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1238 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1239 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1240 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1245 }// loop over recparticles
1247 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1248 if(fTCARandomConesOut){
1249 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1250 // rescale the momntum vectors for the random cones
1252 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1254 Double_t etaC = rC->Eta();
1255 Double_t phiC = rC->Phi();
1256 // massless jet, unit vector
1257 pTC = rC->ChargedBgEnergy();
1258 if(pTC<=0)pTC = 0.001; // for almost empty events
1259 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1260 Double_t pZC = pTC/TMath::Tan(thetaC);
1261 Double_t pXC = pTC * TMath::Cos(phiC);
1262 Double_t pYC = pTC * TMath::Sin(phiC);
1263 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1264 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1265 rC->SetBgEnergy(0,0);
1266 rC->SetEffArea(jetArea,0);
1270 if(fTCARandomConesOutRan){
1271 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1272 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1275 Double_t etaC = rC->Eta();
1276 Double_t phiC = rC->Phi();
1277 // massless jet, unit vector
1278 pTC = rC->ChargedBgEnergy();
1279 if(pTC<=0)pTC = 0.001;// for almost empty events
1280 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1281 Double_t pZC = pTC/TMath::Tan(thetaC);
1282 Double_t pXC = pTC * TMath::Cos(phiC);
1283 Double_t pYC = pTC * TMath::Sin(phiC);
1284 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1285 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1286 rC->SetBgEnergy(0,0);
1287 rC->SetEffArea(jetArea,0);
1291 }// if(fNRandomCones
1293 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1299 if(fAODJetBackgroundOut){
1300 vector<fastjet::PseudoJet> jets2=sortedJets;
1301 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1304 Double_t meanarea1=0.;
1307 Double_t meanarea2=0.;
1309 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1310 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1312 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1313 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1315 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1316 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1317 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1318 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1327 // fill track information
1328 Int_t nTrackOver = recParticles.GetSize();
1329 // do the same for tracks and jets
1332 TIterator *recIter = recParticles.MakeIterator();
1333 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1334 Float_t pt = tmpRec->Pt();
1336 // Printf("Leading track p_t %3.3E",pt);
1337 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1338 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1339 while(pt<ptCut&&tmpRec){
1341 tmpRec = (AliVParticle*)(recIter->Next());
1346 if(nTrackOver<=0)break;
1347 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1351 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1352 Float_t phi = leading->Phi();
1353 if(phi<0)phi+=TMath::Pi()*2.;
1354 Float_t eta = leading->Eta();
1356 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1357 Float_t tmpPt = tmpRec->Pt();
1358 Float_t tmpEta = tmpRec->Eta();
1359 fh1PtTracksRecIn->Fill(tmpPt);
1360 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1361 if(tmpRec==leading){
1362 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1363 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1367 Float_t tmpPhi = tmpRec->Phi();
1369 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1370 Float_t dPhi = phi - tmpPhi;
1371 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1372 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1373 Float_t dEta = eta - tmpRec->Eta();
1374 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1375 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1376 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1381 // find the random jets
1383 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1385 // fill the jet information from random track
1386 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1387 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1389 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1391 // loop over all jets an fill information, first on is the leading jet
1393 Int_t nRecOverRan = inclusiveJetsRan.size();
1394 Int_t nRecRan = inclusiveJetsRan.size();
1396 if(inclusiveJetsRan.size()>0){
1397 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1398 Float_t pt = leadingJet.Pt();
1401 TLorentzVector vecarearanb;
1403 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1404 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1405 while(pt<ptCut&&iCount<nRecRan){
1409 pt = sortedJetsRan[iCount].perp();
1412 if(nRecOverRan<=0)break;
1413 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1415 Float_t phi = leadingJet.Phi();
1416 if(phi<0)phi+=TMath::Pi()*2.;
1417 pt = leadingJet.Pt();
1419 // correlation of leading jet with random tracks
1421 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1423 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1425 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1426 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1427 Float_t dPhi = phi - tmpPhi;
1428 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1429 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1430 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1431 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1434 Int_t nAodOutJetsRan = 0;
1435 AliAODJet *aodOutJetRan = 0;
1436 for(int j = 0; j < nRecRan;j++){
1437 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1438 Float_t tmpPt = tmpRec.Pt();
1439 fh1PtJetsRecInRan->Fill(tmpPt);
1440 // Fill Spectra with constituents
1441 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1442 fh1NConstRecRan->Fill(constituents.size());
1443 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1444 fh2PtNchRan->Fill(nCh,tmpPt);
1445 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1448 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1449 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1450 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1451 aodOutJetRan->GetRefTracks()->Clear();
1452 aodOutJetRan->SetEffArea(arearan,0);
1453 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1454 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1455 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1460 Float_t tmpPhi = tmpRec.Phi();
1461 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1464 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1465 fh1NConstLeadingRecRan->Fill(constituents.size());
1466 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1472 if(fAODJetBackgroundOut){
1475 Double_t meanarea3=0.;
1476 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1477 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1478 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1480 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1481 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1490 // do the event selection if activated
1491 if(fJetTriggerPtCut>0){
1492 bool select = false;
1493 Float_t minPt = fJetTriggerPtCut;
1495 // hard coded for now ...
1496 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1497 if(cent<10)minPt = 50;
1498 else if(cent<30)minPt = 42;
1499 else if(cent<50)minPt = 28;
1500 else if(cent<80)minPt = 18;
1503 if(externalBackground)rho = externalBackground->GetBackground(2);
1505 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1506 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1507 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1516 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1517 fh1CentralitySelect->Fill(cent);
1518 fh1ZSelect->Fill(zVtx);
1519 aodH->SetFillAOD(kTRUE);
1523 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1524 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1525 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1526 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1528 PostData(1, fHistList);
1531 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1534 // Terminate analysis
1536 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1538 if(fMomResH1Fit) delete fMomResH1Fit;
1539 if(fMomResH2Fit) delete fMomResH2Fit;
1540 if(fMomResH3Fit) delete fMomResH3Fit;
1545 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1548 // get list of tracks/particles for different types
1551 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1554 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1555 if(type!=kTrackAODextraonly) {
1556 AliAODEvent *aod = 0;
1557 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1558 else aod = AODEvent();
1560 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1563 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1564 AliAODTrack *tr = aod->GetTrack(it);
1565 Bool_t bGood = false;
1566 if(fFilterType == 0)bGood = true;
1567 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1568 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1569 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1570 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1573 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1574 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1577 if(tr->Pt()<fTrackPtCut){
1578 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1581 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1586 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1587 AliAODEvent *aod = 0;
1588 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1589 else aod = AODEvent();
1594 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1595 if(!aodExtraTracks)return iCount;
1596 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1597 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1598 if (!track) continue;
1600 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1601 if(!trackAOD)continue;
1602 Bool_t bGood = false;
1603 if(fFilterType == 0)bGood = true;
1604 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1605 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1606 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1607 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1608 if(trackAOD->Pt()<fTrackPtCut) continue;
1609 list->Add(trackAOD);
1614 else if (type == kTrackKineAll||type == kTrackKineCharged){
1615 AliMCEvent* mcEvent = MCEvent();
1616 if(!mcEvent)return iCount;
1617 // we want to have alivpartilces so use get track
1618 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1619 if(!mcEvent->IsPhysicalPrimary(it))continue;
1620 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1621 if(type == kTrackKineAll){
1622 if(part->Pt()<fTrackPtCut)continue;
1626 else if(type == kTrackKineCharged){
1627 if(part->Particle()->GetPDG()->Charge()==0)continue;
1628 if(part->Pt()<fTrackPtCut)continue;
1634 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1635 AliAODEvent *aod = 0;
1636 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1637 else aod = AODEvent();
1638 if(!aod)return iCount;
1639 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1640 if(!tca)return iCount;
1641 for(int it = 0;it < tca->GetEntriesFast();++it){
1642 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1643 if(!part->IsPhysicalPrimary())continue;
1644 if(type == kTrackAODMCAll){
1645 if(part->Pt()<fTrackPtCut)continue;
1649 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1650 if(part->Charge()==0)continue;
1651 if(part->Pt()<fTrackPtCut)continue;
1652 if(kTrackAODMCCharged){
1656 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1667 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1670 // set mom res profiles
1673 fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
1674 fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
1675 fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
1678 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1680 // set tracking efficiency histos
1683 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1684 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1685 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1688 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1691 // Get smearing on generated momentum
1694 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1696 TProfile *fMomRes = 0x0;
1697 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1698 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1699 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1706 Double_t smear = 0.;
1709 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1710 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1711 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1715 Int_t bin = fMomRes->FindBin(pt);
1717 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1721 if(fMomRes) delete fMomRes;
1726 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1728 // Fit linear function on momentum resolution at high pT
1731 if(!fMomResH1Fit && fMomResH1) {
1732 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1733 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1734 fMomResH1Fit ->SetRange(5.,100.);
1737 if(!fMomResH2Fit && fMomResH2) {
1738 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1739 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1740 fMomResH2Fit ->SetRange(5.,100.);
1743 if(!fMomResH3Fit && fMomResH3) {
1744 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1745 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1746 fMomResH3Fit ->SetRange(5.,100.);
1752 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1753 for(int i = 0; i < particles.GetEntries(); i++){
1754 AliVParticle *vp = (AliVParticle*)particles.At(i);
1755 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1756 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1757 jInp.set_user_index(i);
1758 inputParticles.push_back(jInp);