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"
71 ClassImp(AliAnalysisTaskJetCluster)
73 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
81 if(fTCAJetsOut)fTCAJetsOut->Delete();
84 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
85 delete fTCAJetsOutRan;
87 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
88 delete fTCARandomConesOut;
90 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
91 delete fTCARandomConesOutRan;
93 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
94 delete fAODJetBackgroundOut;
97 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
102 fUseAODTrackInput(kFALSE),
103 fUseAODMCInput(kFALSE),
104 fUseBackgroundCalc(kFALSE),
105 fEventSelection(kFALSE),
107 fFilterMaskBestPt(0),
110 fTrackTypeRec(kTrackUndef),
111 fTrackTypeGen(kTrackUndef),
113 fNSkipLeadingCone(0),
117 fTrackEtaWindow(0.9),
120 fJetOutputMinPt(0.150),
121 fMaxTrackPtInJet(100.),
128 fBackgroundBranch(""),
139 fUseTrMomentumSmearing(kFALSE),
140 fUseDiceEfficiency(kFALSE),
142 fAlgorithm(fastjet::kt_algorithm),
143 fStrategy(fastjet::Best),
144 fRecombScheme(fastjet::BIpt_scheme),
145 fAreaType(fastjet::active_area),
147 fActiveAreaRepeats(1),
151 fTCARandomConesOut(0x0),
152 fTCARandomConesOutRan(0x0),
153 fAODJetBackgroundOut(0x0),
159 fh1PtHardTrials(0x0),
162 fh1NConstLeadingRec(0x0),
164 fh1PtJetsLeadingRecIn(0x0),
165 fh1PtJetConstRec(0x0),
166 fh1PtJetConstLeadingRec(0x0),
167 fh1PtTracksRecIn(0x0),
168 fh1PtTracksLeadingRecIn(0x0),
170 fh1NConstRecRan(0x0),
171 fh1PtJetsLeadingRecInRan(0x0),
172 fh1NConstLeadingRecRan(0x0),
173 fh1PtJetsRecInRan(0x0),
174 fh1PtTracksGenIn(0x0),
176 fh1CentralityPhySel(0x0),
178 fh1CentralitySelect(0x0),
183 fh2NRecTracksPt(0x0),
185 fh2NConstLeadingPt(0x0),
187 fh2LeadingJetPhiEta(0x0),
189 fh2LeadingJetEtaPt(0x0),
191 fh2LeadingTrackEtaPt(0x0),
192 fh2JetsLeadingPhiEta(0x0),
193 fh2JetsLeadingPhiPt(0x0),
194 fh2TracksLeadingPhiEta(0x0),
195 fh2TracksLeadingPhiPt(0x0),
196 fh2TracksLeadingJetPhiPt(0x0),
197 fh2JetsLeadingPhiPtW(0x0),
198 fh2TracksLeadingPhiPtW(0x0),
199 fh2TracksLeadingJetPhiPtW(0x0),
200 fh2NRecJetsPtRan(0x0),
202 fh2NConstLeadingPtRan(0x0),
207 fh2TracksLeadingJetPhiPtRan(0x0),
208 fh2TracksLeadingJetPhiPtWRan(0x0),
209 fh2PtGenPtSmeared(0x0),
211 fp1PtResolution(0x0),
218 for(int i = 0;i<3;i++){
219 fh1BiARandomCones[i] = 0;
220 fh1BiARandomConesRan[i] = 0;
222 for(int i = 0;i<kMaxCent;i++){
223 fh2JetsLeadingPhiPtC[i] = 0;
224 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
225 fh2TracksLeadingJetPhiPtC[i] = 0;
226 fh2TracksLeadingJetPhiPtWC[i] = 0;
230 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
231 AliAnalysisTaskSE(name),
235 fUseAODTrackInput(kFALSE),
236 fUseAODMCInput(kFALSE),
237 fUseBackgroundCalc(kFALSE),
238 fEventSelection(kFALSE), fFilterMask(0),
239 fFilterMaskBestPt(0),
242 fTrackTypeRec(kTrackUndef),
243 fTrackTypeGen(kTrackUndef),
245 fNSkipLeadingCone(0),
249 fTrackEtaWindow(0.9),
252 fJetOutputMinPt(0.150),
253 fMaxTrackPtInJet(100.),
260 fBackgroundBranch(""),
271 fUseTrMomentumSmearing(kFALSE),
272 fUseDiceEfficiency(kFALSE),
274 fAlgorithm(fastjet::kt_algorithm),
275 fStrategy(fastjet::Best),
276 fRecombScheme(fastjet::BIpt_scheme),
277 fAreaType(fastjet::active_area),
279 fActiveAreaRepeats(1),
283 fTCARandomConesOut(0x0),
284 fTCARandomConesOutRan(0x0),
285 fAODJetBackgroundOut(0x0),
291 fh1PtHardTrials(0x0),
294 fh1NConstLeadingRec(0x0),
296 fh1PtJetsLeadingRecIn(0x0),
297 fh1PtJetConstRec(0x0),
298 fh1PtJetConstLeadingRec(0x0),
299 fh1PtTracksRecIn(0x0),
300 fh1PtTracksLeadingRecIn(0x0),
302 fh1NConstRecRan(0x0),
303 fh1PtJetsLeadingRecInRan(0x0),
304 fh1NConstLeadingRecRan(0x0),
305 fh1PtJetsRecInRan(0x0),
306 fh1PtTracksGenIn(0x0),
308 fh1CentralityPhySel(0x0),
310 fh1CentralitySelect(0x0),
315 fh2NRecTracksPt(0x0),
317 fh2NConstLeadingPt(0x0),
319 fh2LeadingJetPhiEta(0x0),
321 fh2LeadingJetEtaPt(0x0),
323 fh2LeadingTrackEtaPt(0x0),
324 fh2JetsLeadingPhiEta(0x0),
325 fh2JetsLeadingPhiPt(0x0),
326 fh2TracksLeadingPhiEta(0x0),
327 fh2TracksLeadingPhiPt(0x0),
328 fh2TracksLeadingJetPhiPt(0x0),
329 fh2JetsLeadingPhiPtW(0x0),
330 fh2TracksLeadingPhiPtW(0x0),
331 fh2TracksLeadingJetPhiPtW(0x0),
332 fh2NRecJetsPtRan(0x0),
334 fh2NConstLeadingPtRan(0x0),
339 fh2TracksLeadingJetPhiPtRan(0x0),
340 fh2TracksLeadingJetPhiPtWRan(0x0),
341 fh2PtGenPtSmeared(0x0),
343 fp1PtResolution(0x0),
350 for(int i = 0;i<3;i++){
351 fh1BiARandomCones[i] = 0;
352 fh1BiARandomConesRan[i] = 0;
354 for(int i = 0;i<kMaxCent;i++){
355 fh2JetsLeadingPhiPtC[i] = 0;
356 fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
357 fh2TracksLeadingJetPhiPtC[i] = 0;
358 fh2TracksLeadingJetPhiPtWC[i] = 0;
360 DefineOutput(1, TList::Class());
365 Bool_t AliAnalysisTaskJetCluster::Notify()
368 // Implemented Notify() to read the cross sections
369 // and number of trials from pyxsec.root
374 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
378 // Create the output container
381 fRandom = new TRandom3(0);
387 if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
391 if(fNonStdBranch.Length()!=0)
393 // only create the output branch if we have a name
394 // Create a new branch for jets...
395 // -> cleared in the UserExec....
396 // here we can also have the case that the brnaches are written to a separate file
399 fTCAJetsOut = new TClonesArray("AliAODJet", 0);
400 fTCAJetsOut->SetName(fNonStdBranch.Data());
401 AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
404 if(fJetTypes&kJetRan){
405 fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
406 fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
407 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
408 fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%d",fNonStdBranch.Data(),"random",fUseTrMomentumSmearing,fUseDiceEfficiency));
409 AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
412 if(fUseBackgroundCalc){
413 if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
414 fAODJetBackgroundOut = new AliAODJetEventBackground();
415 fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
416 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
417 fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency));
419 AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
422 // create the branch for the random cones with the same R
423 TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
424 if(fUseDiceEfficiency || fUseTrMomentumSmearing)
425 cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
429 if(!AODEvent()->FindListObject(cName.Data())){
430 fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
431 fTCARandomConesOut->SetName(cName.Data());
432 AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
435 // create the branch with the random for the random cones on the random event
436 if(fJetTypes&kRCRan){
437 cName = Form("%sRandomCone_random",fNonStdBranch.Data());
438 if(!AODEvent()->FindListObject(cName.Data())){
439 fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
440 fTCARandomConesOutRan->SetName(cName.Data());
441 AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
446 if(fNonStdFile.Length()!=0){
448 // case that we have an AOD extension we need to fetch the jets from the extended output
449 // we identify the extension aod event by looking for the branchname
450 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
451 // case that we have an AOD extension we need can fetch the background maybe from the extended output
452 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
456 // FitMomentumResolution();
459 if(!fHistList)fHistList = new TList();
460 fHistList->SetOwner();
461 PostData(1, fHistList); // post data in any case once
463 Bool_t oldStatus = TH1::AddDirectoryStatus();
464 TH1::AddDirectory(kFALSE);
469 const Int_t nBinPt = 100;
470 Double_t binLimitsPt[nBinPt+1];
471 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
473 binLimitsPt[iPt] = 0.0;
476 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
480 const Int_t nBinPhi = 90;
481 Double_t binLimitsPhi[nBinPhi+1];
482 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
484 binLimitsPhi[iPhi] = -1.*TMath::Pi();
487 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
493 const Int_t nBinEta = 40;
494 Double_t binLimitsEta[nBinEta+1];
495 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
497 binLimitsEta[iEta] = -2.0;
500 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
504 const int nChMax = 5000;
506 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
507 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
509 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
510 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
513 fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
514 fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
516 fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
517 fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
518 fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
519 fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
522 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
523 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
524 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
526 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
527 fh1PtJetsRecInRan = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
528 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
529 fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
530 fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
531 fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
532 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
533 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
534 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
535 fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
537 fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
538 fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
539 fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
541 fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
542 fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
543 fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
545 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
546 fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
547 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
551 fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
552 fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
553 fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
554 fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
556 fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
557 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);
558 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);
559 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);
563 fh2JetPhiEta = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
564 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
565 fh2LeadingJetPhiEta = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
566 nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
568 fh2JetEtaPt = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
569 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
570 fh2LeadingJetEtaPt = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
571 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
573 fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
574 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
575 fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
576 nBinEta,binLimitsEta,nBinPt,binLimitsPt);
580 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
581 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
582 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
583 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
584 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
585 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
586 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
587 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
588 fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
589 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
590 fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
591 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
593 fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
594 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
595 fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
596 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
598 fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
599 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
600 fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
601 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
603 //Detector level effects histos
604 fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
606 fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
607 fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
609 fHistList->Add(fh2PtGenPtSmeared);
610 fHistList->Add(fp1Efficiency);
611 fHistList->Add(fp1PtResolution);
613 if(fNRandomCones>0&&fUseBackgroundCalc){
614 for(int i = 0;i<3;i++){
615 fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
616 fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
620 for(int i = 0;i < kMaxCent;i++){
621 fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
622 fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
623 fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
624 fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
627 const Int_t saveLevel = 3; // large save level more histos
629 fHistList->Add(fh1Xsec);
630 fHistList->Add(fh1Trials);
632 fHistList->Add(fh1NJetsRec);
633 fHistList->Add(fh1NConstRec);
634 fHistList->Add(fh1NConstLeadingRec);
635 fHistList->Add(fh1PtJetsRecIn);
636 fHistList->Add(fh1PtJetsLeadingRecIn);
637 fHistList->Add(fh1PtTracksRecIn);
638 fHistList->Add(fh1PtTracksLeadingRecIn);
639 fHistList->Add(fh1PtJetConstRec);
640 fHistList->Add(fh1PtJetConstLeadingRec);
641 fHistList->Add(fh1NJetsRecRan);
642 fHistList->Add(fh1NConstRecRan);
643 fHistList->Add(fh1PtJetsLeadingRecInRan);
644 fHistList->Add(fh1NConstLeadingRecRan);
645 fHistList->Add(fh1PtJetsRecInRan);
646 fHistList->Add(fh1Nch);
647 fHistList->Add(fh1Centrality);
648 fHistList->Add(fh1CentralitySelect);
649 fHistList->Add(fh1CentralityPhySel);
650 fHistList->Add(fh1Z);
651 fHistList->Add(fh1ZSelect);
652 fHistList->Add(fh1ZPhySel);
653 if(fNRandomCones>0&&fUseBackgroundCalc){
654 for(int i = 0;i<3;i++){
655 fHistList->Add(fh1BiARandomCones[i]);
656 fHistList->Add(fh1BiARandomConesRan[i]);
659 for(int i = 0;i < kMaxCent;i++){
660 fHistList->Add(fh2JetsLeadingPhiPtC[i]);
661 fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
662 fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
663 fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
666 fHistList->Add(fh2NRecJetsPt);
667 fHistList->Add(fh2NRecTracksPt);
668 fHistList->Add(fh2NConstPt);
669 fHistList->Add(fh2NConstLeadingPt);
670 fHistList->Add(fh2PtNch);
671 fHistList->Add(fh2PtNchRan);
672 fHistList->Add(fh2PtNchN);
673 fHistList->Add(fh2PtNchNRan);
674 fHistList->Add(fh2JetPhiEta);
675 fHistList->Add(fh2LeadingJetPhiEta);
676 fHistList->Add(fh2JetEtaPt);
677 fHistList->Add(fh2LeadingJetEtaPt);
678 fHistList->Add(fh2TrackEtaPt);
679 fHistList->Add(fh2LeadingTrackEtaPt);
680 fHistList->Add(fh2JetsLeadingPhiEta );
681 fHistList->Add(fh2JetsLeadingPhiPt);
682 fHistList->Add(fh2TracksLeadingPhiEta);
683 fHistList->Add(fh2TracksLeadingPhiPt);
684 fHistList->Add(fh2TracksLeadingJetPhiPt);
685 fHistList->Add(fh2JetsLeadingPhiPtW);
686 fHistList->Add(fh2TracksLeadingPhiPtW);
687 fHistList->Add(fh2TracksLeadingJetPhiPtW);
688 fHistList->Add(fh2NRecJetsPtRan);
689 fHistList->Add(fh2NConstPtRan);
690 fHistList->Add(fh2NConstLeadingPtRan);
691 fHistList->Add(fh2TracksLeadingJetPhiPtRan);
692 fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
695 // =========== Switch on Sumw2 for all histos ===========
696 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
697 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
702 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
705 TH1::AddDirectory(oldStatus);
708 void AliAnalysisTaskJetCluster::Init()
714 if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
716 FitMomentumResolution();
720 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
723 // handle and reset the output jet branch
725 if(fTCAJetsOut)fTCAJetsOut->Delete();
726 if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
727 if(fTCARandomConesOut)fTCARandomConesOut->Delete();
728 if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
729 if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
731 AliAODJetEventBackground* externalBackground = 0;
732 if(!externalBackground&&fBackgroundBranch.Length()){
733 externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
734 if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
735 if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
738 // Execute analysis for current event
740 AliESDEvent *fESD = 0;
741 if(fUseAODTrackInput){
742 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
744 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
750 // assume that the AOD is in the general output...
753 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
757 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
761 //Check if information is provided detector level effects
762 if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
763 if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
765 Bool_t selectEvent = false;
766 Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
772 const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
773 TString vtxTitle(vtxAOD->GetTitle());
774 zVtx = vtxAOD->GetZ();
776 cent = fAOD->GetHeader()->GetCentrality();
777 if(cent<10)cenClass = 0;
778 else if(cent<30)cenClass = 1;
779 else if(cent<50)cenClass = 2;
780 else if(cent<80)cenClass = 3;
781 if(physicsSelection){
782 fh1CentralityPhySel->Fill(cent);
783 fh1ZPhySel->Fill(zVtx);
787 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
788 Float_t yvtx = vtxAOD->GetY();
789 Float_t xvtx = vtxAOD->GetX();
790 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
791 if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
792 if(physicsSelection){
798 if(cent<fCentCutLo||cent>fCentCutUp){
809 PostData(1, fHistList);
812 fh1Centrality->Fill(cent);
814 fh1Trials->Fill("#sum{ntrials}",1);
817 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
819 // ==== General variables needed
823 // we simply fetch the tracks/mc particles as a list of AliVParticles
828 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
829 Float_t nCh = recParticles.GetEntries();
831 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
832 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
833 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
837 vector<fastjet::PseudoJet> inputParticlesRec;
838 vector<fastjet::PseudoJet> inputParticlesRecRan;
840 // Generate the random cones
842 AliAODJet vTmpRan(1,0,0,1);
843 for(int i = 0; i < recParticles.GetEntries(); i++){
844 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
846 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
847 // we take total momentum here
849 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
850 //Add particles to fastjet in case we are not running toy model
851 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
852 jInp.set_user_index(i);
853 inputParticlesRec.push_back(jInp);
855 else if(fUseDiceEfficiency) {
857 // Dice to decide if particle is kept or not - toy model for efficiency
859 Double_t rnd = fRandom->Uniform(1.);
860 Double_t pT = vp->Pt();
861 Double_t eff[3] = {0.};
863 if(pT>10.) pTtmp = 10.;
864 Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
865 Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
866 Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
868 //Sort efficiencies from large to small
869 if(eff1>eff2 && eff1>eff3) {
884 else if(eff2>eff1 && eff2>eff3) {
899 else if(eff3>eff1 && eff3>eff2) {
915 Double_t sumEff = eff[0]+eff[1]+eff[2];
916 fp1Efficiency->Fill(vp->Pt(),sumEff);
917 if(rnd>sumEff) continue;
919 if(fUseTrMomentumSmearing) {
920 //Smear momentum of generated particle
922 //Select hybrid track category
924 smear = GetMomentumSmearing(cat[2],pT);
925 else if(rnd<=(eff[2]+eff[1]))
926 smear = GetMomentumSmearing(cat[1],pT);
928 smear = GetMomentumSmearing(cat[0],pT);
930 fp1PtResolution->Fill(vp->Pt(),smear);
932 Double_t sigma = vp->Pt()*smear;
933 Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
935 Double_t phi = vp->Phi();
936 Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
937 Double_t pX = pTrec * TMath::Cos(phi);
938 Double_t pY = pTrec * TMath::Sin(phi);
939 Double_t pZ = pTrec/TMath::Tan(theta);
940 Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
942 fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
944 fastjet::PseudoJet jInp(pX,pY,pZ,p);
945 jInp.set_user_index(i);
946 inputParticlesRec.push_back(jInp);
950 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
951 jInp.set_user_index(i);
952 inputParticlesRec.push_back(jInp);
958 // the randomized input changes eta and phi, but keeps the p_T
959 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
960 Double_t pT = vp->Pt();
961 Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
962 Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
964 Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
965 Double_t pZ = pT/TMath::Tan(theta);
967 Double_t pX = pT * TMath::Cos(phi);
968 Double_t pY = pT * TMath::Sin(phi);
969 Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
970 fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
972 jInpRan.set_user_index(i);
973 inputParticlesRecRan.push_back(jInpRan);
974 vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
977 // fill the tref array, only needed when we write out jets
980 fRef->Delete(); // make sure to delete before placement new...
981 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
982 new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
985 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
989 if(inputParticlesRec.size()==0){
990 if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
991 PostData(1, fHistList);
996 // employ setters for these...
999 // now create the object that holds info about ghosts
1001 if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1002 // reduce CPU time...
1004 fActiveAreaRepeats = 0;
1008 fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1009 fastjet::AreaType areaType = fastjet::active_area;
1010 fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1011 fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1012 fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1014 //range where to compute background
1015 Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1017 phiMax = 2*TMath::Pi();
1018 rapMax = fGhostEtamax - fRparam;
1019 rapMin = - fGhostEtamax + fRparam;
1020 fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1023 const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1024 const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1027 fh1NJetsRec->Fill(sortedJets.size());
1029 // loop over all jets an fill information, first on is the leading jet
1031 Int_t nRecOver = inclusiveJets.size();
1032 Int_t nRec = inclusiveJets.size();
1033 if(inclusiveJets.size()>0){
1034 AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1035 Double_t area = clustSeq.area(sortedJets[0]);
1036 leadingJet.SetEffArea(area,0);
1037 Float_t pt = leadingJet.Pt();
1038 Int_t nAodOutJets = 0;
1039 Int_t nAodOutTracks = 0;
1040 AliAODJet *aodOutJet = 0;
1043 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1044 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1045 while(pt<ptCut&&iCount<nRec){
1049 pt = sortedJets[iCount].perp();
1052 if(nRecOver<=0)break;
1053 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1055 Float_t phi = leadingJet.Phi();
1056 if(phi<0)phi+=TMath::Pi()*2.;
1057 Float_t eta = leadingJet.Eta();
1059 if(externalBackground){
1060 // carefull has to be filled in a task before
1061 // todo, ReArrange to the botom
1062 pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1064 pt = leadingJet.Pt() - pTback;
1065 // correlation of leading jet with tracks
1066 TIterator *recIter = recParticles.MakeIterator();
1068 AliVParticle *tmpRecTrack = 0;
1069 while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1070 Float_t tmpPt = tmpRecTrack->Pt();
1072 Float_t tmpPhi = tmpRecTrack->Phi();
1073 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1074 Float_t dPhi = phi - tmpPhi;
1075 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1076 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1077 fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1078 fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1080 fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1081 fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1087 TLorentzVector vecareab;
1088 for(int j = 0; j < nRec;j++){
1089 AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1092 Float_t tmpPt = tmpRec.Pt();
1094 if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1095 aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1096 aodOutJet->GetRefTracks()->Clear();
1097 Double_t area1 = clustSeq.area(sortedJets[j]);
1098 aodOutJet->SetEffArea(area1,0);
1099 fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
1100 vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
1101 aodOutJet->SetVectorAreaCharged(&vecareab);
1105 Float_t tmpPtBack = 0;
1106 if(externalBackground){
1107 // carefull has to be filled in a task before
1108 // todo, ReArrange to the botom
1109 tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1111 tmpPt = tmpPt - tmpPtBack;
1112 if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1114 fh1PtJetsRecIn->Fill(tmpPt);
1115 // Fill Spectra with constituentsemacs
1116 const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1118 fh1NConstRec->Fill(constituents.size());
1119 fh2PtNch->Fill(nCh,tmpPt);
1120 fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1121 fh2NConstPt->Fill(tmpPt,constituents.size());
1122 // loop over constiutents and fill spectrum
1124 AliVParticle *partLead = 0;
1125 Float_t ptLead = -1;
1127 for(unsigned int ic = 0; ic < constituents.size();ic++){
1128 AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1130 fh1PtJetConstRec->Fill(part->Pt());
1132 if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1133 if(part->Pt()>fMaxTrackPtInJet){
1134 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1137 if(part->Pt()>ptLead){
1138 ptLead = part->Pt();
1141 if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1143 //set pT of leading constituent of jet
1144 aodOutJet->SetPtLeading(partLead->Pt());
1146 AliAODTrack *aodT = 0;
1148 aodT = dynamic_cast<AliAODTrack*>(partLead);
1150 if(aodT->TestFilterBit(fFilterMaskBestPt)){
1151 aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1157 Float_t tmpPhi = tmpRec.Phi();
1158 Float_t tmpEta = tmpRec.Eta();
1159 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1161 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1162 fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1163 fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1164 fh1NConstLeadingRec->Fill(constituents.size());
1165 fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1168 fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1169 fh2JetEtaPt->Fill(tmpEta,tmpPt);
1170 Float_t dPhi = phi - tmpPhi;
1171 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1172 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1173 Float_t dEta = eta - tmpRec.Eta();
1174 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1175 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1177 fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1178 fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1180 fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1181 }// loop over reconstructed jets
1186 // Add the random cones...
1187 if(fNRandomCones>0&&fTCARandomConesOut){
1188 // create a random jet within the acceptance
1189 Double_t etaMax = fTrackEtaWindow - fRparam;
1192 Double_t pTC = 1; // small number
1193 for(int ir = 0;ir < fNRandomCones;ir++){
1194 Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1195 Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1197 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1198 Double_t pZC = pTC/TMath::Tan(thetaC);
1199 Double_t pXC = pTC * TMath::Cos(phiC);
1200 Double_t pYC = pTC * TMath::Sin(phiC);
1201 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1202 AliAODJet tmpRecC (pXC,pYC,pZC, pC);
1204 for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1205 AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1206 if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1211 // test for overlap with previous cones to avoid double counting
1212 for(int iic = 0;iic<ir;iic++){
1213 AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1215 if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1222 tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1223 tmpRecC.SetPtLeading(-1.);
1224 if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1225 if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1226 }// loop over random cones creation
1229 // loop over the reconstructed particles and add up the pT in the random cones
1230 // maybe better to loop over randomized particles not in the real jets...
1231 // but this by definition brings dow average energy in the whole event
1232 AliAODJet vTmpRanR(1,0,0,1);
1233 for(int i = 0; i < recParticles.GetEntries(); i++){
1234 AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1235 if(fTCARandomConesOut){
1236 for(int ir = 0;ir < fNRandomCones;ir++){
1237 AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
1238 if(jC&&jC->DeltaR(vp)<fRparam){
1239 if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1240 jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1241 if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1244 }// add up energy in cone
1246 // the randomized input changes eta and phi, but keeps the p_T
1247 if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1248 Double_t pTR = vp->Pt();
1249 Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1250 Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1252 Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
1253 Double_t pZR = pTR/TMath::Tan(thetaR);
1255 Double_t pXR = pTR * TMath::Cos(phiR);
1256 Double_t pYR = pTR * TMath::Sin(phiR);
1257 Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
1258 vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1259 if(fTCARandomConesOutRan){
1260 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1261 AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1262 if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1263 if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1264 jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1265 if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1270 }// loop over recparticles
1272 Float_t jetArea = fRparam*fRparam*TMath::Pi();
1273 if(fTCARandomConesOut){
1274 for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1275 // rescale the momentum vectors for the random cones
1277 AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1279 Double_t etaC = rC->Eta();
1280 Double_t phiC = rC->Phi();
1281 // massless jet, unit vector
1282 pTC = rC->ChargedBgEnergy();
1283 if(pTC<=0)pTC = 0.001; // for almost empty events
1284 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1285 Double_t pZC = pTC/TMath::Tan(thetaC);
1286 Double_t pXC = pTC * TMath::Cos(phiC);
1287 Double_t pYC = pTC * TMath::Sin(phiC);
1288 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1289 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1290 rC->SetBgEnergy(0,0);
1291 rC->SetEffArea(jetArea,0);
1295 if(fTCARandomConesOutRan){
1296 for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1297 AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1300 Double_t etaC = rC->Eta();
1301 Double_t phiC = rC->Phi();
1302 // massless jet, unit vector
1303 pTC = rC->ChargedBgEnergy();
1304 if(pTC<=0)pTC = 0.001;// for almost empty events
1305 Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
1306 Double_t pZC = pTC/TMath::Tan(thetaC);
1307 Double_t pXC = pTC * TMath::Cos(phiC);
1308 Double_t pYC = pTC * TMath::Sin(phiC);
1309 Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
1310 rC->SetPxPyPzE(pXC,pYC,pZC, pC);
1311 rC->SetBgEnergy(0,0);
1312 rC->SetEffArea(jetArea,0);
1316 }// if(fNRandomCones
1318 //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1324 if(fAODJetBackgroundOut){
1325 vector<fastjet::PseudoJet> jets2=sortedJets;
1326 if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
1329 Double_t meanarea1=0.;
1332 Double_t meanarea2=0.;
1334 clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1335 fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1337 // fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));
1338 // fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
1340 clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1341 fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1342 // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
1343 // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
1352 // fill track information
1353 Int_t nTrackOver = recParticles.GetSize();
1354 // do the same for tracks and jets
1357 TIterator *recIter = recParticles.MakeIterator();
1358 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1359 Float_t pt = tmpRec->Pt();
1361 // Printf("Leading track p_t %3.3E",pt);
1362 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1363 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1364 while(pt<ptCut&&tmpRec){
1366 tmpRec = (AliVParticle*)(recIter->Next());
1371 if(nTrackOver<=0)break;
1372 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1376 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1377 Float_t phi = leading->Phi();
1378 if(phi<0)phi+=TMath::Pi()*2.;
1379 Float_t eta = leading->Eta();
1381 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1382 Float_t tmpPt = tmpRec->Pt();
1383 Float_t tmpEta = tmpRec->Eta();
1384 fh1PtTracksRecIn->Fill(tmpPt);
1385 fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1386 if(tmpRec==leading){
1387 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1388 fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1392 Float_t tmpPhi = tmpRec->Phi();
1394 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1395 Float_t dPhi = phi - tmpPhi;
1396 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1397 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1398 Float_t dEta = eta - tmpRec->Eta();
1399 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1400 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1401 fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1406 // find the random jets
1408 fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1410 // fill the jet information from random track
1411 const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1412 const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
1414 fh1NJetsRecRan->Fill(sortedJetsRan.size());
1416 // loop over all jets an fill information, first on is the leading jet
1418 Int_t nRecOverRan = inclusiveJetsRan.size();
1419 Int_t nRecRan = inclusiveJetsRan.size();
1421 if(inclusiveJetsRan.size()>0){
1422 AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1423 Float_t pt = leadingJet.Pt();
1426 TLorentzVector vecarearanb;
1428 for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1429 Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1430 while(pt<ptCut&&iCount<nRecRan){
1434 pt = sortedJetsRan[iCount].perp();
1437 if(nRecOverRan<=0)break;
1438 fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1440 Float_t phi = leadingJet.Phi();
1441 if(phi<0)phi+=TMath::Pi()*2.;
1442 pt = leadingJet.Pt();
1444 // correlation of leading jet with random tracks
1446 for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1448 Float_t tmpPt = inputParticlesRecRan[ip].perp();
1450 Float_t tmpPhi = inputParticlesRecRan[ip].phi();
1451 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1452 Float_t dPhi = phi - tmpPhi;
1453 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1454 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1455 fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1456 fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1459 Int_t nAodOutJetsRan = 0;
1460 AliAODJet *aodOutJetRan = 0;
1461 for(int j = 0; j < nRecRan;j++){
1462 AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1463 Float_t tmpPt = tmpRec.Pt();
1464 fh1PtJetsRecInRan->Fill(tmpPt);
1465 // Fill Spectra with constituents
1466 const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1467 fh1NConstRecRan->Fill(constituents.size());
1468 fh2NConstPtRan->Fill(tmpPt,constituents.size());
1469 fh2PtNchRan->Fill(nCh,tmpPt);
1470 fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1473 if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1474 aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1475 Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1476 aodOutJetRan->GetRefTracks()->Clear();
1477 aodOutJetRan->SetEffArea(arearan,0);
1478 fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
1479 vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1480 aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1485 Float_t tmpPhi = tmpRec.Phi();
1486 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1489 fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1490 fh1NConstLeadingRecRan->Fill(constituents.size());
1491 fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1497 if(fAODJetBackgroundOut){
1500 Double_t meanarea3=0.;
1501 clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1502 fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1503 // float areaRandomCone = rRandomCone2 *TMath::Pi();
1505 fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
1506 fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
1515 // do the event selection if activated
1516 if(fJetTriggerPtCut>0){
1517 bool select = false;
1518 Float_t minPt = fJetTriggerPtCut;
1520 // hard coded for now ...
1521 // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1522 if(cent<10)minPt = 50;
1523 else if(cent<30)minPt = 42;
1524 else if(cent<50)minPt = 28;
1525 else if(cent<80)minPt = 18;
1528 if(externalBackground)rho = externalBackground->GetBackground(2);
1530 for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1531 AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1532 Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1541 static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1542 fh1CentralitySelect->Fill(cent);
1543 fh1ZSelect->Fill(zVtx);
1544 aodH->SetFillAOD(kTRUE);
1548 if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1549 if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1550 if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1551 if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1553 PostData(1, fHistList);
1556 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1559 // Terminate analysis
1561 if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1563 if(fMomResH1Fit) delete fMomResH1Fit;
1564 if(fMomResH2Fit) delete fMomResH2Fit;
1565 if(fMomResH3Fit) delete fMomResH3Fit;
1570 Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1573 // get list of tracks/particles for different types
1576 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1579 if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1580 if(type!=kTrackAODextraonly) {
1581 AliAODEvent *aod = 0;
1582 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1583 else aod = AODEvent();
1585 if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1589 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1590 AliAODTrack *tr = aod->GetTrack(it);
1591 Bool_t bGood = false;
1592 if(fFilterType == 0)bGood = true;
1593 else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1594 else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1595 if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1596 if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
1599 if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1600 if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1603 if(tr->Pt()<fTrackPtCut){
1604 if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1607 if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
1612 if(type==kTrackAODextra || type==kTrackAODextraonly) {
1613 AliAODEvent *aod = 0;
1614 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1615 else aod = AODEvent();
1620 TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1621 if(!aodExtraTracks)return iCount;
1622 for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1623 AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1624 if (!track) continue;
1626 AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1627 if(!trackAOD)continue;
1628 Bool_t bGood = false;
1629 if(fFilterType == 0)bGood = true;
1630 else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1631 else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1632 if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1633 if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1634 if(trackAOD->Pt()<fTrackPtCut) continue;
1635 list->Add(trackAOD);
1640 else if (type == kTrackKineAll||type == kTrackKineCharged){
1641 AliMCEvent* mcEvent = MCEvent();
1642 if(!mcEvent)return iCount;
1643 // we want to have alivpartilces so use get track
1644 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1645 if(!mcEvent->IsPhysicalPrimary(it))continue;
1646 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1647 if(type == kTrackKineAll){
1648 if(part->Pt()<fTrackPtCut)continue;
1652 else if(type == kTrackKineCharged){
1653 if(part->Particle()->GetPDG()->Charge()==0)continue;
1654 if(part->Pt()<fTrackPtCut)continue;
1660 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1661 AliAODEvent *aod = 0;
1662 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1663 else aod = AODEvent();
1664 if(!aod)return iCount;
1665 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1666 if(!tca)return iCount;
1667 for(int it = 0;it < tca->GetEntriesFast();++it){
1668 AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1669 if(!part->IsPhysicalPrimary())continue;
1670 if(type == kTrackAODMCAll){
1671 if(part->Pt()<fTrackPtCut)continue;
1675 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1676 if(part->Charge()==0)continue;
1677 if(part->Pt()<fTrackPtCut)continue;
1678 if(kTrackAODMCCharged){
1682 if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1693 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
1696 // set mom res profiles
1699 fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
1700 fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
1701 fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
1704 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
1706 // set tracking efficiency histos
1709 fhEffH1 = (TH1*)h1->Clone("fhEffH1");
1710 fhEffH2 = (TH1*)h2->Clone("fhEffH2");
1711 fhEffH3 = (TH1*)h3->Clone("fhEffH3");
1714 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
1717 // Get smearing on generated momentum
1720 //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
1722 TProfile *fMomRes = 0x0;
1723 if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
1724 if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
1725 if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
1732 Double_t smear = 0.;
1735 if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
1736 if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
1737 if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
1741 Int_t bin = fMomRes->FindBin(pt);
1743 smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
1747 if(fMomRes) delete fMomRes;
1752 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
1754 // Fit linear function on momentum resolution at high pT
1757 if(!fMomResH1Fit && fMomResH1) {
1758 fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
1759 fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
1760 fMomResH1Fit ->SetRange(5.,100.);
1763 if(!fMomResH2Fit && fMomResH2) {
1764 fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
1765 fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
1766 fMomResH2Fit ->SetRange(5.,100.);
1769 if(!fMomResH3Fit && fMomResH3) {
1770 fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
1771 fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
1772 fMomResH3Fit ->SetRange(5.,100.);
1778 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
1779 for(int i = 0; i < particles.GetEntries(); i++){
1780 AliVParticle *vp = (AliVParticle*)particles.At(i);
1781 // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1782 fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
1783 jInp.set_user_index(i);
1784 inputParticles.push_back(jInp);