//
// Emcal jet finder task.
//
-// Authors: C.Loizides, S.Aiola
+// Authors: C.Loizides, S.Aiola, M. Verweij
#include <vector>
#include "AliEmcalJetTask.h"
#include "AliVCluster.h"
#include "AliVEvent.h"
#include "AliVParticle.h"
+#include "AliRhoParameter.h"
using std::cout;
using std::endl;
using std::cerr;
fTracksName("Tracks"),
fCaloName("CaloClusters"),
fJetsName("Jets"),
+ fJetsSubName(""),
fJetType(kNone),
fConstSel(0),
fMCConstSel(0),
fMinMCLabel(0),
fRecombScheme(fastjet::pt_scheme),
fTrackEfficiency(1.),
+ fMCFlag(0),
+ fGeneratorIndex(-1),
fIsInit(0),
+ fLocked(0),
fIsPSelSet(0),
fIsMcPart(0),
fIsEmcPart(0),
fLegacyMode(kFALSE),
fCodeDebug(kFALSE),
+ fDoGenericSubtractionJetMass(kFALSE),
+ fDoGenericSubtractionGR(kFALSE),
+ fDoConstituentSubtraction(kFALSE),
+ fUseExternalBkg(kFALSE),
+ fRhoName(""),
+ fRhomName(""),
+ fRho(0),
+ fRhom(0),
+ fRMax(0.4),
+ fDRStep(0.04),
+ fPtMinGR(40.),
fJets(0),
+ fJetsSub(0),
fEvent(0),
fTracks(0),
- fClus(0)
+ fClus(0),
+ fRhoParam(0),
+ fRhomParam(0)
{
// Default constructor.
}
fTracksName("Tracks"),
fCaloName("CaloClusters"),
fJetsName("Jets"),
+ fJetsSubName(""),
fJetType(kAKT|kFullJet|kRX1Jet),
fConstSel(0),
fMCConstSel(0),
fMinMCLabel(0),
fRecombScheme(fastjet::pt_scheme),
fTrackEfficiency(1.),
+ fMCFlag(0),
+ fGeneratorIndex(-1),
fIsInit(0),
+ fLocked(0),
fIsPSelSet(0),
fIsMcPart(0),
fIsEmcPart(0),
fLegacyMode(kFALSE),
fCodeDebug(kFALSE),
+ fDoGenericSubtractionJetMass(kFALSE),
+ fDoGenericSubtractionGR(kFALSE),
+ fDoConstituentSubtraction(kFALSE),
+ fUseExternalBkg(kFALSE),
+ fRhoName(""),
+ fRhomName(""),
+ fRho(0),
+ fRhom(0),
+ fRMax(0.4),
+ fDRStep(0.04),
+ fPtMinGR(40.),
fJets(0),
+ fJetsSub(0),
fEvent(0),
fTracks(0),
- fClus(0)
+ fClus(0),
+ fRhoParam(0),
+ fRhomParam(0)
{
// Standard constructor.
fJets = new TClonesArray("AliEmcalJet");
fJets->SetName(fJetsName);
+
+ if(!fJetsSubName.IsNull()) {
+ fJetsSub = new TClonesArray("AliEmcalJet");
+ fJetsSub->SetName(fJetsSubName);
+ }
}
//________________________________________________________________________
// clear the jet array (normally a null operation)
fJets->Delete();
+ if(fJetsSub) fJetsSub->Delete();
FindJets();
}
return;
}
+ if (fRhoParam)
+ fRho = fRhoParam->GetVal();
+ if (fRhomParam)
+ fRhom = fRhomParam->GetVal();
+
TString name("kt");
fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
if ((fJetType & kAKT) != 0) {
else {
AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
}
- }
+ }
+ if ((t->GetFlag() & fMCFlag) != fMCFlag) {
+ AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
+ continue;
+ }
+ if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) {
+ AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
+ continue;
+ }
if (t->Pt() < fMinJetTrackPt)
continue;
Double_t eta = t->Eta();
// run jet finder
fjw.Run();
+ //run generic subtractor
+ if(fDoGenericSubtractionJetMass) {
+ fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+ fjw.DoGenericSubtractionJetMass();
+ }
+
+ //run constituent subtractor
+ if(fDoConstituentSubtraction) {
+ fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+ fjw.DoConstituentSubtraction();
+ }
+
// get geometry
AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
if (!geom) {
AliEmcalJet *jet = new ((*fJets)[jetCount])
AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
+ jet->SetLabel(ij);
+
+ //do generic subtraction if requested
+#ifdef FASTJET_VERSION
+ if(fDoGenericSubtractionJetMass) {
+ std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
+ Int_t n = (Int_t)jetMassInfo.size();
+ if(n>ij && n>0) {
+ jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
+ jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
+ jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
+ jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
+ }
+ }
+ //here do generic subtraction for angular structure function
+ Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
+ fRMax = fRadius+0.2;
+ if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
+ fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+ fjw.SetRMaxAndStep(fRMax,fDRStep);
+ fjw.DoGenericSubtractionGR(ij);
+ std::vector<double> num = fjw.GetGRNumerator();
+ std::vector<double> den = fjw.GetGRDenominator();
+ std::vector<double> nums = fjw.GetGRNumeratorSub();
+ std::vector<double> dens = fjw.GetGRDenominatorSub();
+ //pass this to AliEmcalJet
+ jet->SetGRNumSize(num.size());
+ jet->SetGRDenSize(den.size());
+ jet->SetGRNumSubSize(nums.size());
+ jet->SetGRDenSubSize(dens.size());
+ Int_t nsize = (Int_t)num.size();
+ for(Int_t g = 0; g<nsize; ++g) {
+ jet->AddGRNumAt(num[g],g);
+ jet->AddGRNumSubAt(nums[g],g);
+ }
+ Int_t dsize = (Int_t)den.size();
+ for(Int_t g = 0; g<dsize; ++g) {
+ jet->AddGRDenAt(den[g],g);
+ jet->AddGRDenSubAt(dens[g],g);
+ }
+ }
+#endif
// loop over constituents
std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
jetCount++;
}
//fJets->Sort();
+
+ //run constituent subtractor if requested
+ if(fDoConstituentSubtraction) {
+#ifdef FASTJET_VERSION
+ if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
+ else {
+ std::vector<fastjet::PseudoJet> jets_sub;
+ jets_sub = fjw.GetConstituentSubtrJets();
+ AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
+ for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
+ //Only storing 4-vector and jet area of unsubtracted jet
+ AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount])
+ AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
+ jet_sub->SetLabel(ijet);
+ fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
+ jet_sub->SetArea(area.perp());
+ jet_sub->SetAreaEta(area.eta());
+ jet_sub->SetAreaPhi(area.phi());
+ jet_sub->SetAreaEmc(area.perp());
+ jetCount++;
+ }
+ }
+#endif
+ } //constituent subtraction
}
//________________________________________________________________________
return 0;
}
+ if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
+ fEvent->AddObject(fJetsSub);
+
if (fTracksName == "Tracks")
am->LoadBranch("Tracks");
if (!fTracks && !fTracksName.IsNull()) {
fIsEmcPart = 1;
}
+ if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
+ fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
+ if (!fRhoParam) {
+ AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
+ return 0;
+ }
+ }
+ if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
+ fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
+ if (!fRhomParam) {
+ AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
+ return 0;
+ }
+ }
+
return 1;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff