-// $Id: $
+// $Id$
//
-// Calculation of rho
+// Calculation of rho from a collection of jets.
+// If scale function is given the scaled rho will be exported
+// with the name as "fRhoName".apppend("_Scaled").
//
// Authors: R.Reed, S.Aiola
+#include "AliAnalysisTaskRho.h"
-#include <TList.h>
-#include <TH1F.h>
-#include <TH2F.h>
#include <TClonesArray.h>
#include <TF1.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TList.h>
+#include <TMath.h>
-#include "AliLog.h"
#include "AliAnalysisManager.h"
-#include "AliVEventHandler.h"
#include "AliCentrality.h"
#include "AliEmcalJet.h"
+#include "AliLog.h"
+#include "AliRhoParameter.h"
#include "AliVCluster.h"
-
-#include "AliAnalysisTaskRho.h"
+#include "AliVEventHandler.h"
ClassImp(AliAnalysisTaskRho)
//________________________________________________________________________
AliAnalysisTaskRho::AliAnalysisTaskRho() :
AliAnalysisTaskRhoBase(),
- fTracksName("tracks"),
- fJetsName("KtJets"),
- fClustersName("caloClusters"),
+ fTracksName(),
+ fJetsName(),
fRhoScaledName(""),
fPhiMin(0),
fPhiMax(0),
fEtaMin(0),
fEtaMax(0),
fAreaCut(0),
+ fAreaEmcCut(0),
+ fNExclLeadJets(0),
fScaleFunction(0),
+ fCreateHisto(kFALSE),
+ fTracks(0),
+ fJets(0),
fOutputList(0),
fHistCentrality(0),
fHistJetPt(0),
fHistJetPtvsNtrack(0),
fHistJetAreavsNtrack(0),
fHistNjetvsNtrack(0),
- fRhoScaled(0),
- fNewRhoFunction(0)
+ fRhoScaled(0)
{
- // Constructor
+ // Constructor.
}
+
//________________________________________________________________________
-AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name) :
+AliAnalysisTaskRho::AliAnalysisTaskRho(const char *name, Bool_t histo) :
AliAnalysisTaskRhoBase(name),
fTracksName("tracks"),
fJetsName("KtJets"),
- fClustersName("caloClusters"),
fRhoScaledName(""),
fPhiMin(0),
- fPhiMax(0),
- fEtaMin(0),
- fEtaMax(0),
- fAreaCut(0),
+ fPhiMax(TMath::TwoPi()),
+ fEtaMin(-0.5),
+ fEtaMax(+0.5),
+ fAreaCut(0.01),
+ fAreaEmcCut(0),
+ fNExclLeadJets(1),
fScaleFunction(0),
+ fCreateHisto(histo),
+ fTracks(0),
+ fJets(0),
fOutputList(0),
fHistCentrality(0),
fHistJetPt(0),
fHistJetPtvsNtrack(0),
fHistJetAreavsNtrack(0),
fHistNjetvsNtrack(0),
- fRhoScaled(0),
- fNewRhoFunction(0)
+ fRhoScaled(0)
{
- // Constructor
+ // Constructor.
- DefineOutput(1, TList::Class());
+ if (fCreateHisto)
+ DefineOutput(1, TList::Class());
}
//________________________________________________________________________
void AliAnalysisTaskRho::UserCreateOutputObjects()
{
+ // User create output objects, called at the beginning of the analysis.
+
AliAnalysisTaskRhoBase::UserCreateOutputObjects();
- fRhoScaledName = fRhoName;
- fRhoScaledName += "_Scaled";
- fRhoScaled = new TParameter<Double_t>(fRhoScaledName, 0);
+ if (fScaleFunction) {
+ fRhoScaledName = fRhoName;
+ fRhoScaledName += "_Scaled";
+ fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);
+ }
+
+ if (!fCreateHisto)
+ return;
OpenFile(1);
fOutputList = new TList();
fHistJetPt = new TH1F("JetPt", "Jet Pt", 100, 0, 250);
fHistJetArea = new TH1F("JetArea", "Jet Area", 100, 0.0, 1.0);
-
- fNewRhoFunction = new TF1("rfunc","pol4", -1, 100);
fOutputList->Add(fHistCentrality);
fOutputList->Add(fHistRhovsCent);
fOutputList->Add(fHistJetPtvsNtrack);
fOutputList->Add(fHistJetAreavsNtrack);
fOutputList->Add(fHistNjetvsNtrack);
-
- fOutputList->Add(fNewRhoFunction);
PostData(1, fOutputList);
-
-}
-
-//________________________________________________________________________
-Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
-{
- Double_t scale = 1;
- if (fScaleFunction)
- scale = fScaleFunction->Eval(cent);
- return scale;
}
//________________________________________________________________________
void AliAnalysisTaskRho::UserExec(Option_t *)
{
// Main loop, called for each event.
-
- AliAnalysisTaskRhoBase::UserExec("");
- // add rho to event if not yet there
- if (!(InputEvent()->FindListObject(fRhoScaledName))) {
- new(fRhoScaled) TParameter<Double_t>(fRhoScaledName, 0);
- InputEvent()->AddObject(fRhoScaled);
+ if (!fIsInit) {
+ ExecOnce();
+ fIsInit = 1;
}
- // optimization in case autobranch loading is off
- AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
- if (fTracksName == "Tracks")
- am->LoadBranch("Tracks");
-
- TClonesArray *jets = 0;
- TClonesArray *tracks = 0;
-
- tracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
- if (!tracks) {
- AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
- return;
- }
-
- jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
- if (!jets) {
- AliError(Form("Pointer to jets %s == 0", fTracksName.Data() ));
+ if (!fJets)
return;
- }
- fHistCentrality->Fill(fCent);
+ fRho->SetVal(-1);
+ if (fRhoScaled)
+ fRhoScaled->SetVal(-1);
+
+ DetermineCent();
+
+ const Int_t Njets = fJets->GetEntries();
+
+ Int_t maxJetIds[] = {-1, -1};
+ Float_t maxJetPts[] = { 0, 0};
+
+ if (fNExclLeadJets > 0) {
+ for (Int_t ij = 0; ij < Njets; ++ij) {
+ AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
+ if (!jet) {
+ AliError(Form("%s: Could not receive jet %d", GetName(), ij));
+ continue;
+ }
+ if (jet->Pt() > maxJetPts[0]) {
+ maxJetPts[1] = maxJetPts[0];
+ maxJetIds[1] = maxJetIds[0];
+ maxJetPts[0] = jet->Pt();
+ maxJetIds[0] = ij;
+ } else if (jet->Pt() > maxJetPts[1]) {
+ maxJetPts[1] = jet->Pt();
+ maxJetIds[1] = ij;
+ }
+ }
+ if (fNExclLeadJets < 2) {
+ maxJetIds[1] = -1;
+ maxJetPts[1] = -1;
+ }
+ }
- const Int_t Ntracks = tracks->GetEntries();
- const Int_t Njets = jets->GetEntries();
+ static Double_t rhovec[999];
Int_t NjetAcc = 0;
- vector<Double_t> rhovec;
+ // push all jets within selected acceptance into stack
for (Int_t iJets = 0; iJets < Njets; ++iJets) {
- //push all jets within selected acceptance into stack
- AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
+
+ // exlcuding lead jets
+ if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
+ continue;
+
+ AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
if (!jet)
continue;
+
+ // applying some other cuts
if (jet->Area() < fAreaCut)
continue;
+ if (jet->AreaEmc()<fAreaEmcCut)
+ continue;
if ((jet->Phi() < fPhiMin) || (jet->Phi() > fPhiMax))
continue;
if ((jet->Eta() < fEtaMin) || (jet->Eta() > fEtaMax))
continue;
- if (jet->Area() == 0)
- continue;
- NjetAcc++;
- rhovec.push_back(jet->Pt() / jet->Area());
- fHistJetPt->Fill(jet->Pt());
- fHistJetArea->Fill(jet->Area());
- fHistJetPtvsCent->Fill(fCent, jet->Pt());
- fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
- fHistJetAreavsCent->Fill(fCent, jet->Area());
- fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
- }
-
- fHistNjetvsCent->Fill(fCent,NjetAcc);
- fHistNjetvsNtrack->Fill(Ntracks,NjetAcc);
- Double_t scale = GetScaleFactor(fCent);
- Double_t rhochem = GetRhoFactor(fCent);
+ rhovec[NjetAcc] = jet->Pt() / jet->Area();
- Double_t rho0 = -1;
-
- if (rhovec.size() > 0){
- //find median value
- Sort(rhovec);
- rho0 = GetMedian(rhovec,0);
- fHistRhovsCent->Fill(fCent, rho0);
- fHistDeltaRhovsCent->Fill(fCent, rho0 - rhochem);
- fHistDeltaRhoScalevsCent->Fill(fCent, rho0 * scale - rhochem);
- fHistRhovsNtrack->Fill(Ntracks, rho0);
- fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rhochem);
- fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rho0 * scale - rhochem);
+ ++NjetAcc;
+
+ if (fCreateHisto) {
+ // filling histograms
+ const Int_t Ntracks = fTracks->GetEntriesFast();
+ fHistJetPt->Fill(jet->Pt());
+ fHistJetArea->Fill(jet->Area());
+ fHistJetPtvsCent->Fill(fCent, jet->Pt());
+ fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
+ fHistJetAreavsCent->Fill(fCent, jet->Area());
+ fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
+ }
}
- fRho->SetVal(rho0);
+ if (NjetAcc > 0) {
+ //find median value
+ Double_t rho0 = TMath::Median(NjetAcc, rhovec);
+ fRho->SetVal(rho0);
+ Double_t rhoScaled = rho0;
+ if (fRhoScaled) {
+ rhoScaled *= GetScaleFactor(fCent);
+ fRhoScaled->SetVal(rhoScaled);
+ }
- Double_t rho_scaled = rho0 * scale;
+ if (fCreateHisto) {
+ // filling other histograms
+ Double_t rho = GetRhoFactor(fCent);
+ const Int_t Ntracks = fTracks->GetEntriesFast();
+ fHistRhovsCent->Fill(fCent, rho0);
+ fHistDeltaRhovsCent->Fill(fCent, rho0 - rho);
+ fHistDeltaRhoScalevsCent->Fill(fCent, rhoScaled - rho);
+ fHistRhovsNtrack->Fill(Ntracks, rho0);
+ fHistDeltaRhovsNtrack->Fill(Ntracks, rho0 - rho);
+ fHistDeltaRhoScalevsNtrack->Fill(Ntracks, rhoScaled - rho);
+ }
+ }
- fRhoScaled->SetVal(rho_scaled);
-
- PostData(1, fOutputList);
+ if (fCreateHisto) {
+ const Int_t Ntracks = fTracks->GetEntriesFast();
+ fHistCentrality->Fill(fCent);
+ fHistNjetvsCent->Fill(fCent, NjetAcc);
+ fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
+ PostData(1, fOutputList);
+ }
}
//________________________________________________________________________
-void AliAnalysisTaskRho::Sort(vector<Double_t>& v)
+void AliAnalysisTaskRho::ExecOnce()
{
- vector<Double_t> temp;
- temp.push_back(v[0]);
- for (UInt_t i = 1; i < v.size(); i++) {
- Bool_t insert = kFALSE;
- for (vector<Double_t>::iterator j = temp.begin(); j < temp.end(); j++){
- if (v[i]> * j) {
- temp.insert(j, v[i]);
- insert = kTRUE;
- j = temp.end();
- }
+ // Initialize some settings that need to be determined in UserExec.
+
+ AliAnalysisTaskRhoBase::ExecOnce();
+
+ if (fRhoScaled) {
+ // add rho to event if not yet there
+ if (!(InputEvent()->FindListObject(fRhoScaledName))) {
+ InputEvent()->AddObject(fRhoScaled);
+ } else {
+ AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
+ return;
}
- if (!insert)
- temp.insert(temp.end(), v[i]);
}
- v = temp;
- return;
-}
-//________________________________________________________________________
-Double_t AliAnalysisTaskRho::GetMedian(vector<Double_t> v, Int_t c)
-{
- if (v.size() == 0)
- return -1;
- if ((v.size() < 2) && (c == 1))
- return -1;
- if ((v.size() < 3) && ( c== 2))
- return -1;
-
- if (c == 1)
- v.erase(v.begin());
-
- if (c == 2) {
- v.erase(v.begin());
- v.erase(v.begin());
+ if (fCreateHisto) {
+ fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
+ if (!fTracks) {
+ AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data() ));
+ return;
+ }
}
- Float_t middle = (Float_t)v.size() / 2.0;
- if (middle != ceil(middle)) {
- //odd number
- return v[floor(middle)];
+ fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
+ if (!fJets) {
+ AliError(Form("%s: Pointer to jets %s == 0", GetName(), fJetsName.Data() ));
+ return;
}
- else{
- //even
- return (v[middle] + v[middle-1]) / 2.0;
- }
}
//________________________________________________________________________
-void AliAnalysisTaskRho::Terminate(Option_t *)
+Double_t AliAnalysisTaskRho::GetScaleFactor(Double_t cent)
{
- fHistRhovsCent->Fit(fNewRhoFunction, "NO");
- PostData(1, fOutputList);
+ // Get scale factor.
+
+ Double_t scale = 1;
+ if (fScaleFunction)
+ scale = fScaleFunction->Eval(cent);
+ return scale;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff