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 "fOutRhoName".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"
18 #include "AliJetContainer.h"
20 ClassImp(AliAnalysisTaskRhoSparse)
22 //________________________________________________________________________
23 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse() :
24 AliAnalysisTaskRhoBase("AliAnalysisTaskRhoSparse"),
32 //________________________________________________________________________
33 AliAnalysisTaskRhoSparse::AliAnalysisTaskRhoSparse(const char *name, Bool_t histo) :
34 AliAnalysisTaskRhoBase(name, histo),
42 //________________________________________________________________________
43 void AliAnalysisTaskRhoSparse::UserCreateOutputObjects()
45 if (!fCreateHisto) return;
47 AliAnalysisTaskRhoBase::UserCreateOutputObjects();
49 fHistOccCorrvsCent = new TH2F("OccCorrvsCent", "OccCorrvsCent", 101, -1, 100, 2000, 0 , 2);
50 fOutput->Add(fHistOccCorrvsCent);
53 //________________________________________________________________________
54 Bool_t AliAnalysisTaskRhoSparse::IsJetOverlapping(AliEmcalJet* jet1, AliEmcalJet* jet2)
56 for (Int_t i = 0; i < jet1->GetNumberOfTracks(); ++i)
58 Int_t jet1Track = jet1->TrackAt(i);
59 for (Int_t j = 0; j < jet2->GetNumberOfTracks(); ++j)
61 Int_t jet2Track = jet2->TrackAt(j);
62 if (jet1Track == jet2Track)
69 //________________________________________________________________________
70 Bool_t AliAnalysisTaskRhoSparse::IsJetSignal(AliEmcalJet* jet)
80 //________________________________________________________________________
81 Bool_t AliAnalysisTaskRhoSparse::Run()
87 fOutRhoScaled->SetVal(0);
92 const Int_t Njets = fJets->GetEntries();
94 AliJetContainer *sigjets = static_cast<AliJetContainer*>(fJetCollArray.At(1));
97 if (sigjets) NjetsSig = sigjets->GetNJets();
99 Int_t maxJetIds[] = {-1, -1};
100 Float_t maxJetPts[] = { 0, 0};
102 if (fNExclLeadJets > 0) {
103 for (Int_t ij = 0; ij < Njets; ++ij) {
104 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
106 AliError(Form("%s: Could not receive jet %d", GetName(), ij));
113 if (jet->Pt() > maxJetPts[0]) {
114 maxJetPts[1] = maxJetPts[0];
115 maxJetIds[1] = maxJetIds[0];
116 maxJetPts[0] = jet->Pt();
118 } else if (jet->Pt() > maxJetPts[1]) {
119 maxJetPts[1] = jet->Pt();
123 if (fNExclLeadJets < 2) {
129 static Double_t rhovec[999];
131 Double_t TotaljetArea=0;
132 Double_t TotaljetAreaPhys=0;
134 // push all jets within selected acceptance into stack
135 for (Int_t iJets = 0; iJets < Njets; ++iJets) {
137 // exlcuding lead jets
138 if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
141 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
143 AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
147 TotaljetArea+=jet->Area();
150 TotaljetAreaPhys+=jet->Area();
156 // Search for overlap with signal jets
157 Bool_t isOverlapping = kFALSE;
159 for(Int_t j=0;j<NjetsSig;j++)
161 AliEmcalJet* signalJet = sigjets->GetAcceptJet(j);
164 if(!IsJetSignal(signalJet))
167 if(IsJetOverlapping(signalJet, jet))
169 isOverlapping = kTRUE;
179 rhovec[NjetAcc] = jet->Pt() / jet->Area();
184 Double_t OccCorr=0.0;
185 if(TotaljetArea>0) OccCorr=TotaljetAreaPhys/TotaljetArea;
187 fHistOccCorrvsCent->Fill(fCent, OccCorr);
191 Double_t rho = TMath::Median(NjetAcc, rhovec);
197 fOutRho->SetVal(rho);
200 Double_t rhoScaled = rho * GetScaleFactor(fCent);
201 fOutRhoScaled->SetVal(rhoScaled);