//
// Emcal jet finder task.
//
-// Author: C.Loizides
+// Authors: C.Loizides, S.Aiola
-#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 "AliEmcalJet.h"
#include "AliFJWrapper.h"
#include "AliVCluster.h"
-#include "AliVTrack.h"
-#include "AliEmcalJet.h"
+#include "AliVParticle.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+
+#include "AliEmcalJetTask.h"
ClassImp(AliEmcalJetTask)
fType(0),
fMinJetTrackPt(0.15),
fMinJetClusPt(0.15),
- fJets(0)
+ fMinJetArea(0.01),
+ fMinJetPt(1.0),
+ fJets(0),
+ fEvent(0)
{
// Default constructor.
}
fType(0),
fMinJetTrackPt(0.15),
fMinJetClusPt(0.15),
- fJets(0)
+ fMinJetArea(0.01),
+ fMinJetPt(1.0),
+ fJets(0),
+ fEvent(0)
{
// Standard constructor.
fJets->SetName(fJetsName);
}
+//_____________________________________________________
+TString AliEmcalJetTask::GetBeamType()
+{
+ // Get beam type : pp-AA-pA
+ // ESDs have it directly, AODs get it from hardcoded run number ranges
+
+ TString beamType;
+
+ AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fEvent);
+ if (esd) {
+ const AliESDRun *run = esd->GetESDRun();
+ beamType = run->GetBeamType();
+ }
+ else
+ {
+ Int_t runNumber = fEvent->GetRunNumber();
+ if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
+ (runNumber >= 166529 && runNumber <= 170593)) // LHC11h
+ {
+ beamType = "A-A";
+ }
+ else
+ {
+ beamType = "p-p";
+ }
+ }
+
+ return beamType;
+}
+
+
//________________________________________________________________________
void AliEmcalJetTask::UserExec(Option_t *)
{
// Main loop, called for each event.
+ // get the event
+ fEvent = InputEvent();
+
+ if (!fEvent) {
+ AliError("Could not retrieve event! Returning");
+ return;
+ }
+
// add jets to event if not yet there
- if (!(InputEvent()->FindListObject(fJetsName)))
- InputEvent()->AddObject(fJets);
+ if (!(fEvent->FindListObject(fJetsName)))
+ fEvent->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 ((fType == 0) || (fType == 1)) {
if (fTracksName == "Tracks")
am->LoadBranch("Tracks");
- tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
+ tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
if (!tracks) {
- AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
+ AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
return;
}
}
- if ((fType==0)||(fType==2)) {
+ if ((fType == 0) || (fType == 2)) {
if (fCaloName == "CaloClusters")
am->LoadBranch("CaloClusters");
- clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
+ clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
if (!clus) {
- AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
+ 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;
+
+ // get centrality
+ Double_t cent = 99;
+
+ if (GetBeamType() == "A-A") {
+ AliCentrality *centrality = InputEvent()->GetCentrality();
+
+ if (centrality)
+ cent = centrality->GetCentralityPercentile("V0M");
+ else
+ cent = 99; // probably pp data
+
+ if (cent < 0) {
+ AliWarning(Form("Centrality negative: %f, assuming 99", cent));
+ cent = 99;
+ }
}
-
+
FindJets(tracks, clus, fAlgo, fRadius, cent);
}
}
AliFJWrapper fjw(name, name);
+ fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
fjw.SetR(radius);
fjw.SetAlgorithm(jalgo);
fjw.SetMaxRap(0.9);
if (tracks) {
const Int_t Ntracks = tracks->GetEntries();
for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
- AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
+ AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks));
if (!t)
continue;
- if (t->Pt()<fMinJetTrackPt)
+
+ if (t->Pt() < fMinJetTrackPt)
continue;
- fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
+ fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100); // offset of 100 for consistency with cluster ids
}
}
if (clus) {
Double_t vertex[3] = {0, 0, 0};
- InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+ fEvent->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));
continue;
TLorentzVector nPart;
c->GetMomentum(nPart, vertex);
- Double_t et = nPart.Pt();
- if (et<fMinJetClusPt)
+ if (nPart.Pt() < fMinJetClusPt)
continue;
- fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-1);
+ fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100); // offset of 100 to skip ghost particles uid = -1
}
}
// run jet finder
fjw.Run();
+
+ // get geometry
+ AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
+ if (!geom) {
+ AliFatal("Can not create geometry");
+ return;
+ }
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);
+ fEvent->GetPrimaryVertex()->GetXYZ(vertex);
vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
- Double_t neutralE = 0;Double_t maxTrack = 0;Double_t maxCluster=0;
+ 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 ncharged = 0;
+ Int_t nneutral = 0;
+ Double_t MCpt = 0;
+
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) {
+ 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 ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
+ ++gall;
+ Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
+ Double_t geta = constituents[ic].eta();
+ if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
+ (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
+ ++gemc;
+ continue;
+ } else if (uid > 0) {
+ Int_t tid = uid - 100;
+ AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
if (t) {
- if (t->Pt()>maxTrack)
- maxTrack=t->Pt();
- nt++;
+ if (t->Charge() == 0) {
+ neutralE += t->P();
+ ++nneutral;
+ if (t->Pt() > maxNe)
+ maxNe = t->Pt();
+ } else {
+ ++ncharged;
+ if (t->Pt() > maxCh)
+ maxCh = t->Pt();
+ }
+
+ if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
+ MCpt += t->Pt();
+
+ jet->AddTrackAt(tid, nt);
+ ++nt;
}
} else {
- jet->AddClusterAt(-(uid+1),nc);
- AliVCluster *c = static_cast<AliVCluster*>(clus->At(-(uid+1)));
- if (c) {
- TLorentzVector nP;
- c->GetMomentum(nP, vertex);
- neutralE += nP.P();
- if (nP.Pt()>maxCluster)
- maxCluster=nP.Pt();
- }
- nc++;
+ Int_t cid = -uid - 100;
+ AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
+ if (c) {
+ TLorentzVector nP;
+ c->GetMomentum(nP, vertex);
+ neutralE += nP.P();
+ if (nP.Pt() > maxNe)
+ maxNe = nP.Pt();
+
+ if (c->Chi2() == 100) // MC particle
+ MCpt += nP.Pt();
+
+ jet->AddClusterAt(cid, nc);
+ ++nc;
+ ++nneutral;
+ }
}
}
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);
+ 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++;
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff