1 // $Id: AliAnalysisTaskRho.cxx 58408 2012-09-03 07:00:58Z loizides $
3 // Calculation of rho from a collection of jets.
4 // If scale function is given the scaled rho will be exported
5 // with the name as "fRhoName".apppend("_Scaled").
7 // Authors: R.Reed, S.Aiola, M.Connors
9 #include "AliAnalysisTaskRhoSparse.h"
11 #include <TClonesArray.h>
14 #include "AliAnalysisManager.h"
15 #include "AliEmcalJet.h"
17 #include "AliRhoParameter.h"
19 ClassImp(AliAnalysisTaskRhoSparse)
21 //________________________________________________________________________
22 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse() :
23 AliAnalysisTaskRhoBase("AliAnalysisTaskRhoSparse"),
24 fHistOccCorrvsCent(0),
32 //________________________________________________________________________
33 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse(const char *name, Bool_t histo) :
34 AliAnalysisTaskRhoBase(name, histo),
35 fHistOccCorrvsCent(0),
43 //________________________________________________________________________
44 void AliAnalysisTaskRhoSparse::UserCreateOutputObjects()
49 AliAnalysisTaskRhoBase::UserCreateOutputObjects();
51 fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
52 fOutput->Add(fHistOccCorrvsCent);
55 //________________________________________________________________________
56 Bool_t AliAnalysisTaskRhoSparse::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2)
58 for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
60 Int_t jet1Track = jet1->TrackAt(i);
61 for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
63 Int_t jet2Track = jet2->TrackAt(j);
64 if (jet1Track == jet2Track)
71 //________________________________________________________________________
72 Bool_t AliAnalysisTaskRhoSparse::IsJetSignal(AliEmcalJet* jet)
82 //________________________________________________________________________
83 Bool_t AliAnalysisTaskRhoSparse::Run()
89 fRhoScaled->SetVal(0);
94 const Int_t Njets = fJets->GetEntries();
96 TClonesArray *sigjets = 0;
97 sigjets= dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fSigJetsName));
100 if(sigjets) NjetsSig = sigjets->GetEntries();
103 Int_t maxJetIds[] = {-1, -1};
104 Float_t maxJetPts[] = { 0, 0};
106 if (fNExclLeadJets > 0) {
107 for (Int_t ij = 0; ij < Njets; ++ij) {
108 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
110 AliError(Form("%s: Could not receive jet %d", GetName(), ij));
117 if (jet->Pt() > maxJetPts[0]) {
118 maxJetPts[1] = maxJetPts[0];
119 maxJetIds[1] = maxJetIds[0];
120 maxJetPts[0] = jet->Pt();
122 } else if (jet->Pt() > maxJetPts[1]) {
123 maxJetPts[1] = jet->Pt();
127 if (fNExclLeadJets < 2) {
133 static Double_t rhovec[999];
135 Double_t TotaljetArea=0;
136 Double_t TotaljetAreaPhys=0;
138 // push all jets within selected acceptance into stack
139 for (Int_t iJets = 0; iJets < Njets; ++iJets) {
141 // exlcuding lead jets
142 if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
145 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
147 AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
151 TotaljetArea+=jet->Area();
154 TotaljetAreaPhys+=jet->Area();
162 // Search for overlap with signal jets
163 Bool_t isOverlapping = kFALSE;
165 for(Int_t j=0;j<NjetsSig;j++)
167 AliEmcalJet* signalJet = static_cast<AliEmcalJet*>(sigjets->At(j));
168 if(!AcceptJet(signalJet))
170 if(!IsJetSignal(signalJet))
173 if(IsJetOverlapping(signalJet, jet))
175 isOverlapping = kTRUE;
185 rhovec[NjetAcc] = jet->Pt() / jet->Area();
192 Double_t OccCorr=0.0;
193 if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
195 fHistOccCorrvsCent->Fill(fCent, OccCorr);
200 Double_t rho = TMath::Median(NjetAcc, rhovec);
210 Double_t rhoScaled = rho * GetScaleFactor(fCent);
211 fRhoScaled->SetVal(rhoScaled);