//
// Emcal jet finder task.
//
-// Author: C.Loizides
+// Authors: C.Loizides, S.Aiola
+#include <vector>
#include "AliEmcalJetTask.h"
+
#include <TChain.h>
#include <TClonesArray.h>
-#include <TH1F.h>
-#include <TH2F.h>
#include <TList.h>
#include <TLorentzVector.h>
-#include <TParticle.h>
+
#include "AliAnalysisManager.h"
#include "AliCentrality.h"
+#include "AliEMCALGeometry.h"
+#include "AliESDEvent.h"
+#include "AliEmcalJet.h"
+#include "AliEmcalParticle.h"
#include "AliFJWrapper.h"
+#include "AliMCEvent.h"
#include "AliVCluster.h"
-#include "AliVTrack.h"
-#include "AliEmcalJet.h"
+#include "AliVEvent.h"
+#include "AliVParticle.h"
ClassImp(AliEmcalJetTask)
fType(0),
fMinJetTrackPt(0.15),
fMinJetClusPt(0.15),
- fJets(0)
+ fPhiMin(0),
+ fPhiMax(TMath::TwoPi()),
+ fEtaMin(-1),
+ fEtaMax(1),
+ fMinJetArea(0.01),
+ fMinJetPt(1.0),
+ fIsInit(0),
+ fIsMcPart(0),
+ fIsEmcPart(0),
+ fJets(0),
+ fEvent(0),
+ fTracks(0),
+ fClus(0)
{
// Default constructor.
}
fType(0),
fMinJetTrackPt(0.15),
fMinJetClusPt(0.15),
- fJets(0)
+ fPhiMin(0),
+ fPhiMax(TMath::TwoPi()),
+ fEtaMin(-1),
+ fEtaMax(1),
+ fMinJetArea(0.01),
+ fMinJetPt(1.0),
+ fIsInit(0),
+ fIsMcPart(0),
+ fIsEmcPart(0),
+ fJets(0),
+ fEvent(0),
+ fTracks(0),
+ fClus(0)
{
// Standard constructor.
{
// Main loop, called for each event.
- // add jets to event if not yet there
- if (!(InputEvent()->FindListObject(fJetsName)))
- InputEvent()->AddObject(fJets);
-
- // delete jet output
- fJets->Delete();
-
-
- // get input collections
- TClonesArray *tracks = 0;
- TClonesArray *clus = 0;
- TList *l = InputEvent()->GetList();
- AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
-
- if ((fType==0)||(fType==1)) {
- if (fTracksName == "Tracks")
- am->LoadBranch("Tracks");
- tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
- if (!tracks) {
- AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
+ if (!fIsInit) {
+ if (!DoInit())
return;
- }
+ fIsInit = kTRUE;
}
- if ((fType==0)||(fType==2)) {
- if (fCaloName == "CaloClusters")
- am->LoadBranch("CaloClusters");
- clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
- if (!clus) {
- AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
- return;
- }
- }
-
- // get centrality
- Float_t cent = -1;
- AliCentrality *centrality = InputEvent()->GetCentrality();
- if (centrality)
- cent = centrality->GetCentralityPercentile("V0M");
- else
- cent=99; // probably pp data
- if (cent<0) {
- AliError(Form("Centrality negative: %f", cent));
- return;
- }
+ // delete jet output
+ fJets->Delete();
- FindJets(tracks, clus, fAlgo, fRadius, cent);
+ FindJets();
}
//________________________________________________________________________
}
//________________________________________________________________________
-void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
+void AliEmcalJetTask::FindJets()
{
// Find jets.
- if (!tracks && !clus)
+ if (!fTracks && !fClus)
return;
TString name("kt");
fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
- if (algo>=1) {
+ if (fAlgo>=1) {
name = "antikt";
jalgo = fastjet::antikt_algorithm;
}
+ // setup fj wrapper
AliFJWrapper fjw(name, name);
- fjw.SetR(radius);
+ fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
+ fjw.SetR(fRadius);
fjw.SetAlgorithm(jalgo);
- fjw.SetMaxRap(0.9);
+ fjw.SetMaxRap(1);
fjw.Clear();
- if (tracks) {
- const Int_t Ntracks = tracks->GetEntries();
+ // get primary vertex
+ Double_t vertex[3] = {0, 0, 0};
+ fEvent->GetPrimaryVertex()->GetXYZ(vertex);
+
+ if (fTracks) {
+ const Int_t Ntracks = fTracks->GetEntries();
for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
- AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
+ AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
if (!t)
continue;
- if (t->Pt()<fMinJetTrackPt)
+ if (t->Pt() < fMinJetTrackPt)
+ continue;
+ Double_t eta = t->Eta();
+ Double_t phi = t->Phi();
+ if ((eta<fEtaMin) && (eta>fEtaMax) &&
+ (phi<fPhiMin) && (phi<fPhiMax))
continue;
- fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
+ // offset of 100 for consistency with cluster ids
+ fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
}
}
- if (clus) {
- Double_t vertex[3] = {0, 0, 0};
- InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
- const Int_t Nclus = clus->GetEntries();
- for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
- AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
- if (!c)
- continue;
- if (!c->IsEMCAL())
- continue;
- TLorentzVector nPart;
- c->GetMomentum(nPart, vertex);
- Double_t et = nPart.Pt();
- if (et<fMinJetClusPt)
- continue;
- fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-1);
+ if (fClus) {
+ if (fIsEmcPart) {
+ const Int_t Nclus = fClus->GetEntries();
+ for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
+ AliVCluster *c = 0;
+ Double_t cEta=0,cPhi=0,cPt=0;
+ Double_t cPx=0,cPy=0,cPz=0;
+ if (fIsEmcPart) {
+ AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
+ if (!ep)
+ continue;
+ c = ep->GetCluster();
+ if (!c)
+ continue;
+ cEta = ep->Eta();
+ cPhi = ep->Phi();
+ cPt = ep->Pt();
+ cPx = ep->Px();
+ cPy = ep->Py();
+ cPz = ep->Pz();
+ } else {
+ c = static_cast<AliVCluster*>(fClus->At(iClus));
+ if (!c)
+ continue;
+ TLorentzVector nP;
+ c->GetMomentum(nP, vertex);
+ cEta = nP.Eta();
+ cPhi = nP.Phi();
+ cPt = nP.Pt();
+ cPx = nP.Px();
+ cPy = nP.Py();
+ cPz = nP.Pz();
+ }
+ if (!c->IsEMCAL())
+ continue;
+ if (cPt < fMinJetClusPt)
+ continue;
+ if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
+ (cPhi<fPhiMin) && (cPhi<fPhiMax))
+ continue;
+ // offset of 100 to skip ghost particles uid = -1
+ fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
+ }
+ } else { /*CaloClusters given as input*/
+ const Int_t Nclus = fClus->GetEntries();
+ for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
+ AliVCluster *c = static_cast<AliVCluster*>(fClus->At(iClus));
+ if (!c)
+ continue;
+ if (!c->IsEMCAL())
+ continue;
+ TLorentzVector nPart;
+ c->GetMomentum(nPart, vertex);
+ if (nPart.Pt() < fMinJetClusPt)
+ continue;
+ // offset of 100 to skip ghost particles uid = -1
+ fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);
+ }
}
}
// run jet finder
fjw.Run();
-
+
+ // get geometry
+ AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+ if (!geom) {
+ AliFatal(Form("%s: Can not create geometry", GetName()));
+ return;
+ }
+
+ // loop over fastjet jets
std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
- for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
- if (jets_incl[ij].perp()<1)
+ for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
+ if (jets_incl[ij].perp()<fMinJetPt)
+ continue;
+ if (fjw.GetJetArea(ij)<fMinJetArea)
continue;
+
AliEmcalJet *jet = new ((*fJets)[jetCount])
AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
- jet->SetArea(fjw.GetJetArea(ij));
- Double_t vertex[3] = {0, 0, 0};
- InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
- vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
- Double_t neutralE = 0;Double_t maxTrack = 0;Double_t maxCluster=0;
+
+ // loop over constituents
+ std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
jet->SetNumberOfTracks(constituents.size());
jet->SetNumberOfClusters(constituents.size());
- Int_t nt = 0;
- Int_t nc = 0;
- for(UInt_t ic=0; ic<constituents.size(); ++ic) {
+
+ Int_t nt = 0;
+ Int_t nc = 0;
+ Double_t neutralE = 0;
+ Double_t maxCh = 0;
+ Double_t maxNe = 0;
+ Int_t gall = 0;
+ Int_t gemc = 0;
+ Int_t cemc = 0;
+ Int_t ncharged = 0;
+ Int_t nneutral = 0;
+ Double_t mcpt = 0;
+ Double_t emcpt = 0;
+
+ for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
Int_t uid = constituents[ic].user_index();
- if (uid>=0){
- jet->AddTrackAt(uid, nt);
- AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
- if (t) {
- if (t->Pt()>maxTrack)
- maxTrack=t->Pt();
- nt++;
+
+ if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
+ ++gall;
+ Double_t gphi = constituents[ic].phi();
+ if (gphi<0)
+ gphi += TMath::TwoPi();
+ gphi *= TMath::RadToDeg();
+ Double_t geta = constituents[ic].eta();
+ if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
+ (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
+ ++gemc;
+ } else if ((uid > 0) && fTracks) { // track constituent
+ Int_t tid = uid - 100;
+ AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
+ if (!t)
+ continue;
+ Double_t cEta = t->Eta();
+ Double_t cPhi = t->Phi();
+ Double_t cPt = t->Pt();
+ Double_t cP = t->P();
+ if (t->Charge() == 0) {
+ neutralE += cP;
+ ++nneutral;
+ if (cPt > maxNe)
+ maxNe = cPt;
+ } else {
+ ++ncharged;
+ if (cPt > maxCh)
+ maxCh = cPt;
}
- } else {
- jet->AddClusterAt(-(uid+1),nc);
- AliVCluster *c = static_cast<AliVCluster*>(clus->At(-(uid+1)));
- if (c) {
+
+ if (fIsMcPart || t->GetLabel() == 100) // check if MC particle
+ mcpt += cPt;
+
+ if (cPhi<0)
+ cPhi += TMath::TwoPi();
+ cPhi *= TMath::RadToDeg();
+ if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
+ (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+ emcpt += cPt;
+ ++cemc;
+ }
+
+ jet->AddTrackAt(tid, nt);
+ ++nt;
+ } else if (fClus) { // cluster constituent
+ Int_t cid = -uid - 100;
+ AliVCluster *c = 0;
+ Double_t cEta=0,cPhi=0,cPt=0,cP=0;
+ if (fIsEmcPart) {
+ AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
+ if (!ep)
+ continue;
+ c = ep->GetCluster();
+ if (!c)
+ continue;
+ cEta = ep->Eta();
+ cPhi = ep->Phi();
+ cPt = ep->Pt();
+ cP = ep->P();
+ } else {
+ c = static_cast<AliVCluster*>(fClus->At(cid));
+ if (!c)
+ continue;
TLorentzVector nP;
c->GetMomentum(nP, vertex);
- neutralE += nP.P();
- if (nP.Pt()>maxCluster)
- maxCluster=nP.Pt();
+ cEta = nP.Eta();
+ cPhi = nP.Phi();
+ cPt = nP.Pt();
+ cP = nP.P();
+ }
+
+ neutralE += cP;
+ if (cPt > maxNe)
+ maxNe = cPt;
+
+ if (c->Chi2() == 100) // MC particle
+ mcpt += cPt;
+
+ if (cPhi<0)
+ cPhi += TMath::TwoPi();
+ cPhi *= TMath::RadToDeg();
+ if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
+ (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+ emcpt += cPt;
+ ++cemc;
}
- nc++;
+
+ jet->AddClusterAt(cid, nc);
+ ++nc;
+ ++nneutral;
+ } else {
+ AliError(Form("%s: No logical way to end up here.", GetName()));
+ continue;
}
- }
+ } /* end of constituent loop */
+
jet->SetNumberOfTracks(nt);
jet->SetNumberOfClusters(nc);
- jet->SetMaxTrackPt(maxTrack);
- jet->SetMaxClusterPt(maxCluster);
- jet->SetNEF(neutralE/jet->E());
+ jet->SortConstituents();
+ jet->SetMaxNeutralPt(maxNe);
+ jet->SetMaxChargedPt(maxCh);
+ jet->SetNEF(neutralE / jet->E());
+ jet->SetArea(fjw.GetJetArea(ij));
+ jet->SetNumberOfCharged(ncharged);
+ jet->SetNumberOfNeutrals(nneutral);
+ jet->SetMCPt(mcpt);
+ jet->SetNEmc(cemc);
+ jet->SetPtEmc(emcpt);
+
+ if (gall > 0)
+ jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
+ else
+ jet->SetAreaEmc(-1);
+ if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
+ (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
+ (jet->Eta() > geom->GetArm1EtaMin()) &&
+ (jet->Eta() < geom->GetArm1EtaMax()))
+ jet->SetAxisInEmcal(kTRUE);
jetCount++;
}
}
+
+//________________________________________________________________________
+Bool_t AliEmcalJetTask::DoInit()
+{
+ // Init. Return true if successful.
+
+ // get input collections
+ AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+
+ // get the event
+ fEvent = InputEvent();
+ if (!fEvent) {
+ AliError(Form("%s: Could not retrieve event! Returning", GetName()));
+ return 0;
+ }
+
+ // add jets to event if not yet there
+ if (!(fEvent->FindListObject(fJetsName)))
+ fEvent->AddObject(fJets);
+ else {
+ AliError(Form("%s: Object with name %s already in event! Returning", fJetsName.Data(), GetName()));
+ return 0;
+ }
+
+ if ((fType == 0) || (fType == 1)) {
+ if (fTracksName == "Tracks")
+ am->LoadBranch("Tracks");
+ if (!fTracks && !fTracksName.IsNull()) {
+ fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
+ if (!fTracks) {
+ AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
+ return 0;
+ }
+ }
+ TString objname(fTracks->GetClass()->GetName());
+ TClass cls(objname);
+ if (cls.InheritsFrom("AliMCParticle"))
+ fIsMcPart = 1;
+ }
+
+ if ((fType == 0) || (fType == 2)) {
+ if (fCaloName == "CaloClusters")
+ am->LoadBranch("CaloClusters");
+ if (!fClus && !fCaloName.IsNull()) {
+ fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
+ if (!fClus) {
+ AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
+ return 0;
+ }
+ }
+ TString objname(fClus->GetClass()->GetName());
+ TClass cls(objname);
+ if (cls.InheritsFrom("AliEmcalParticle"))
+ fIsEmcPart = 1;
+ }
+
+ return 1;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff