From 85c423bce787c5ed38c1c6bc535fa6c084593cf4 Mon Sep 17 00:00:00 2001 From: saiola Date: Fri, 12 Dec 2014 16:03:25 -0500 Subject: [PATCH] Reorganization of the EMCal jet finder task: the manin task only performs the basic jet finding; other features (e.g. FJ contribs) are implemented in external utilities such as AliEmcalJetUtilityGenSubtractor that can be attached to AliEmcalJetTask (similar approach as in AliTender). See this twiki for more details: https://twiki.cern.ch/twiki/bin/view/ALICE/EMCalJetFinderTask --- PWG/EMCAL/AliEmcalParticle.cxx | 1 - PWG/EMCAL/AliEmcalParticle.h | 6 +- PWGJE/CMakelibPWGJEEMCALJetTasks.pkg | 3 + PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx | 1015 ++++++----------- PWGJE/EMCALJetTasks/AliEmcalJetTask.h | 249 ++-- PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx | 42 + PWGJE/EMCALJetTasks/AliEmcalJetUtility.h | 36 + .../AliEmcalJetUtilityConstSubtractor.cxx | 195 ++++ .../AliEmcalJetUtilityConstSubtractor.h | 53 + .../AliEmcalJetUtilityGenSubtractor.cxx | 256 +++++ .../AliEmcalJetUtilityGenSubtractor.h | 57 + PWGJE/EMCALJetTasks/AliFJWrapper.h | 4 +- PWGJE/PWGJEEMCALJetTasksLinkDef.h | 3 + 13 files changed, 1123 insertions(+), 797 deletions(-) create mode 100644 PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx create mode 100644 PWGJE/EMCALJetTasks/AliEmcalJetUtility.h create mode 100644 PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx create mode 100644 PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h create mode 100644 PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx create mode 100644 PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h diff --git a/PWG/EMCAL/AliEmcalParticle.cxx b/PWG/EMCAL/AliEmcalParticle.cxx index f6c70a90119..80a793964f3 100644 --- a/PWG/EMCAL/AliEmcalParticle.cxx +++ b/PWG/EMCAL/AliEmcalParticle.cxx @@ -1,4 +1,3 @@ -// $Id$ // // Emcal particle class, which can contain either an AliVTrack or an AliVCluster // diff --git a/PWG/EMCAL/AliEmcalParticle.h b/PWG/EMCAL/AliEmcalParticle.h index 3b2803fec5e..689bed84379 100644 --- a/PWG/EMCAL/AliEmcalParticle.h +++ b/PWG/EMCAL/AliEmcalParticle.h @@ -1,8 +1,6 @@ #ifndef ALIEMCALPARTICLE_H #define ALIEMCALPARTICLE_H -// $Id$ - #include #include #include @@ -69,6 +67,10 @@ class AliEmcalParticle: public AliVParticle { void SetMatchedObj(Int_t id, Double_t d) { ResetMatchedObjects(); fMatchedIds[0] = id; fMatchedDist[0] = d; fNMatched = 1; } void SetMatchedPtr(TObjArray *arr) { fMatchedPtr = arr; } + void SetPt(Double_t pt) { fPt = pt ; } + void SetPhi(Double_t phi) { fPhi = phi; } + void SetEta(Double_t eta) { fEta = eta; } + protected: TLorentzVector &GetLorentzVector(const Double_t *vertex = 0) const; diff --git a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg index 8c8576846d2..a891508816e 100644 --- a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg +++ b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg @@ -126,6 +126,9 @@ set ( SRCS if (FASTJET_FOUND) LIST(APPEND SRCS EMCALJetTasks/AliFJWrapper.cxx + EMCALJetTasks/AliEmcalJetUtility.cxx + EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx + EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx EMCALJetTasks/AliEmcalJetTask.cxx EMCALJetTasks/AliEmcalJetFinder.cxx EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx index 856d5ff7a5c..620ab9d6594 100644 --- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx +++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx @@ -12,7 +12,6 @@ #include #include #include -#include #include "AliAnalysisManager.h" #include "AliCentrality.h" @@ -23,12 +22,9 @@ #include "AliFJWrapper.h" #include "AliMCEvent.h" #include "AliVCluster.h" -#include "AliAODCaloCluster.h" -#include "AliESDCaloCluster.h" #include "AliVEvent.h" #include "AliVParticle.h" -#include "AliRhoParameter.h" -#include "AliPicoTrack.h" +#include "AliEmcalJetUtility.h" using std::cout; using std::endl; @@ -42,13 +38,12 @@ AliEmcalJetTask::AliEmcalJetTask() : fTracksName("Tracks"), fCaloName("CaloClusters"), fJetsName("Jets"), - fJetsSubName(""), - fTracksSubName(""), - fClusSubName(""), - fJetType(kNone), - fConstSel(0), - fMCConstSel(0), - fMarkConst(kFALSE), + fJetType(kAKT|kFullJet|kRX1Jet), + fMinLabelTracks(-kMaxInt), + fMaxLabelTracks(kMaxInt), + fMinLabelClusters(-kMaxInt), + fMaxLabelClusters(kMaxInt), + fMinMCLabel(0), fRadius(0.4), fMinJetTrackPt(0.15), fMinJetClusPt(0.15), @@ -56,51 +51,35 @@ AliEmcalJetTask::AliEmcalJetTask() : fPhiMax(TMath::TwoPi()), fEtaMin(-0.9), fEtaMax(+0.9), - fMinJetArea(0.005), + fMinJetArea(0.001), fMinJetPt(1.0), fJetPhiMin(-10), fJetPhiMax(+10), fJetEtaMin(-1), fJetEtaMax(+1), fGhostArea(0.005), - fMinMCLabel(0), fRecombScheme(fastjet::pt_scheme), fTrackEfficiency(1.), fMCFlag(0), fGeneratorIndex(-1), - fIsInit(0), + fUtilities(0), fLocked(0), + fIsInit(0), fIsPSelSet(0), - fIsMcPart(0), fIsEmcPart(0), - fIsPicoTrack(0), - fIsEsdClus(0), fLegacyMode(kFALSE), - fCodeDebug(kFALSE), - fPionMassClusters(kFALSE), - fDoGenericSubtractionJetMass(kFALSE), - fDoGenericSubtractionGR(kFALSE), - fDoGenericSubtractionExtraJetShapes(kFALSE), - fDoConstituentSubtraction(kFALSE), - fUseExternalBkg(kFALSE), - fRhoName(""), - fRhomName(""), - fRho(0), - fRhom(0), - fRMax(0.4), - fDRStep(0.04), - fPtMinGR(40.), + fGeom(0), fJets(0), - fJetsSub(0), - fTracksSub(0), - fClusSub(0), fEvent(0), fTracks(0), fClus(0), - fRhoParam(0), - fRhomParam(0) + fFastJetWrapper("AliEmcalJetTask","AliEmcalJetTask") { // Default constructor. + + fVertex[0] = 0.; + fVertex[1] = 0.; + fVertex[2] = 0.; } //________________________________________________________________________ @@ -109,13 +88,12 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) : fTracksName("Tracks"), fCaloName("CaloClusters"), fJetsName("Jets"), - fJetsSubName(""), - fTracksSubName(""), - fClusSubName(""), fJetType(kAKT|kFullJet|kRX1Jet), - fConstSel(0), - fMCConstSel(0), - fMarkConst(kFALSE), + fMinLabelTracks(-kMaxInt), + fMaxLabelTracks(kMaxInt), + fMinLabelClusters(-kMaxInt), + fMaxLabelClusters(kMaxInt), + fMinMCLabel(0), fRadius(0.4), fMinJetTrackPt(0.15), fMinJetClusPt(0.15), @@ -130,46 +108,30 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) : fJetEtaMin(-1), fJetEtaMax(+1), fGhostArea(0.005), - fMinMCLabel(0), fRecombScheme(fastjet::pt_scheme), fTrackEfficiency(1.), fMCFlag(0), fGeneratorIndex(-1), - fIsInit(0), + fUtilities(0), fLocked(0), + fIsInit(0), fIsPSelSet(0), - fIsMcPart(0), fIsEmcPart(0), - fIsPicoTrack(0), - fIsEsdClus(0), fLegacyMode(kFALSE), - fCodeDebug(kFALSE), - fPionMassClusters(kFALSE), - fDoGenericSubtractionJetMass(kFALSE), - fDoGenericSubtractionGR(kFALSE), - fDoGenericSubtractionExtraJetShapes(kFALSE), - fDoConstituentSubtraction(kFALSE), - fUseExternalBkg(kFALSE), - fRhoName(""), - fRhomName(""), - fRho(0), - fRhom(0), - fRMax(0.4), - fDRStep(0.04), - fPtMinGR(40.), + fGeom(0), fJets(0), - fJetsSub(0), - fTracksSub(0), - fClusSub(0), fEvent(0), fTracks(0), fClus(0), - fRhoParam(0), - fRhomParam(0) + fFastJetWrapper(name,name) { // Standard constructor. fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex."; + + fVertex[0] = 0.; + fVertex[1] = 0.; + fVertex[2] = 0.; } //________________________________________________________________________ @@ -179,28 +141,57 @@ AliEmcalJetTask::~AliEmcalJetTask() } //________________________________________________________________________ -void AliEmcalJetTask::UserCreateOutputObjects() +AliEmcalJetUtility* AliEmcalJetTask::AddUtility(AliEmcalJetUtility* utility) { - // Create user objects. + // Addition of utilities. - fJets = new TClonesArray("AliEmcalJet"); - fJets->SetName(fJetsName); + if (!fUtilities) fUtilities = new TObjArray(); + if (fUtilities->FindObject(utility)) { + Error("AddSupply", "Jet utility %s already connected.", utility->GetName()); + return utility; + } + fUtilities->Add(utility); + utility->SetJetTask(this); - if(!fJetsSubName.IsNull()) { - fJetsSub = new TClonesArray("AliEmcalJet"); - fJetsSub->SetName(fJetsSubName); - } + return utility; +} - if (fDoConstituentSubtraction && !fTracksSubName.IsNull()) { - fTracksSub = new TClonesArray(); - fTracksSub->SetName(fTracksSubName); - } +//________________________________________________________________________ +void AliEmcalJetTask::InitUtilities() +{ + TIter next(fUtilities); + AliEmcalJetUtility *utility = 0; + while ((utility=static_cast(next()))) utility->Init(); +} - if (fDoConstituentSubtraction && !fClusSubName.IsNull()) { - fClusSub = new TClonesArray(); - fClusSub->SetName(fClusSubName); - } - +//________________________________________________________________________ +void AliEmcalJetTask::PrepareUtilities() +{ + TIter next(fUtilities); + AliEmcalJetUtility *utility = 0; + while ((utility=static_cast(next()))) utility->Prepare(fFastJetWrapper); +} + +//________________________________________________________________________ +void AliEmcalJetTask::ExecuteUtilities(AliEmcalJet* jet, Int_t ij) +{ + TIter next(fUtilities); + AliEmcalJetUtility *utility = 0; + while ((utility=static_cast(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper); +} + +//________________________________________________________________________ +void AliEmcalJetTask::TerminateUtilities() +{ + TIter next(fUtilities); + AliEmcalJetUtility *utility = 0; + while ((utility=static_cast(next()))) utility->Terminate(fFastJetWrapper); +} + +//________________________________________________________________________ +void AliEmcalJetTask::UserCreateOutputObjects() +{ + // Create user objects. } //________________________________________________________________________ @@ -215,117 +206,68 @@ void AliEmcalJetTask::UserExec(Option_t *) // clear the jet array (normally a null operation) fJets->Delete(); - if(fJetsSub) fJetsSub->Delete(); - if(fTracksSub) fTracksSub->Delete(); - if(fClusSub) fClusSub->Delete(); - //Run jet finding + // get primary vertex + if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(fVertex); + FindJets(); -} -//________________________________________________________________________ -void AliEmcalJetTask::Terminate(Option_t *) -{ - // Called once at the end of the analysis. + FillJetBranch(); } //________________________________________________________________________ void AliEmcalJetTask::FindJets() { // Find jets. + if (!fTracks && !fClus){ - AliInfo("WARNING NO TRACKS OR CLUSTERS"); + AliError("No tracks or clusters, returning."); return; } - if (fRhoParam) - fRho = fRhoParam->GetVal(); - if (fRhomParam) - fRhom = fRhomParam->GetVal(); - - TString name("kt"); - fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm); - if ((fJetType & kAKT) != 0) { - name = "antikt"; - jalgo = fastjet::antikt_algorithm; - AliDebug(1,"Using AKT algorithm"); - } - else { - AliDebug(1,"Using KT algorithm"); - } - - if ((fJetType & kR020Jet) != 0) - fRadius = 0.2; - else if ((fJetType & kR030Jet) != 0) - fRadius = 0.3; - else if ((fJetType & kR040Jet) != 0) - fRadius = 0.4; - - // setup fj wrapper - AliFJWrapper fjw(name, name); - fjw.SetAreaType(fastjet::active_area_explicit_ghosts); - fjw.SetGhostArea(fGhostArea); - fjw.SetR(fRadius); - fjw.SetAlgorithm(jalgo); - fjw.SetRecombScheme(static_cast(fRecombScheme)); - fjw.SetMaxRap(fEtaMax); - fjw.Clear(); - - // get primary vertex - Double_t vertex[3] = {0, 0, 0}; - if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex); + fFastJetWrapper.Clear(); AliDebug(2,Form("Jet type = %d", fJetType)); - if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) { + if (fTracks) { const Int_t Ntracks = fTracks->GetEntries(); for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { AliVParticle *t = static_cast(fTracks->At(iTracks)); - if (!t) + + if (!t) continue; + + if (t->Pt() < fMinJetTrackPt) continue; + Double_t eta = t->Eta(); + Double_t phi = t->Phi(); + if ((etafEtaMax) || + (phifPhiMax)) + continue; + + if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) { + AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks)); continue; - if (fIsMcPart) { - if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) { - AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks)); - continue; - } - if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) { - AliDebug(2,Form("Skipping track %d because it is charged.", iTracks)); - continue; - } } - if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) { - if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) { - AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask))); - continue; - } - else { - AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask))); - } + + if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) { + AliDebug(2,Form("Skipping track %d because it is charged.", iTracks)); + continue; } - else { - if (t->TestBits(fConstSel) != (Int_t)fConstSel) { - AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask))); - continue; - } - else { - AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask))); - } + + Int_t lab = TMath::Abs(t->GetLabel()); + if (lab < fMinLabelTracks || lab > fMaxLabelTracks) { + AliDebug(2,Form("Skipping track %d because label %d is not in range (%d, %d)", iTracks, lab, fMinLabelTracks, fMaxLabelTracks)); + continue; } + 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(); - Double_t phi = t->Phi(); - if ((etafEtaMax) || - (phifPhiMax)) - continue; // artificial inefficiency if (fTrackEfficiency < 1.) { @@ -337,12 +279,12 @@ void AliEmcalJetTask::FindJets() } // offset of 100 for consistency with cluster ids - AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt())); - fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100); + AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, lab, t->Pt())); + fFastJetWrapper.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100); } } - if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) { + if (fClus) { const Int_t Nclus = fClus->GetEntries(); for (Int_t iClus = 0; iClus < Nclus; ++iClus) { AliVCluster *c = 0; @@ -350,31 +292,9 @@ void AliEmcalJetTask::FindJets() Double_t cPx=0,cPy=0,cPz=0; if (fIsEmcPart) { AliEmcalParticle *ep = static_cast(fClus->At(iClus)); - if (!ep) - continue; - + if (!ep) continue; c = ep->GetCluster(); - if (!c) - continue; - - if (c->GetLabel() > fMinMCLabel) { - if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) { - AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask))); - continue; - } - else { - AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask))); - } - } - else { - if (ep->TestBits(fConstSel) != (Int_t)fConstSel) { - AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask))); - continue; - } - else { - AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask))); - } - } + if (!c) continue; cEta = ep->Eta(); cPhi = ep->Phi(); @@ -382,32 +302,13 @@ void AliEmcalJetTask::FindJets() cPx = ep->Px(); cPy = ep->Py(); cPz = ep->Pz(); - } else { + } + else { c = static_cast(fClus->At(iClus)); - if (!c) - continue; - - if (c->GetLabel() > fMinMCLabel) { - if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) { - AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask))); - continue; - } - else { - AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask))); - } - } - else { - if (c->TestBits(fConstSel) != (Int_t)fConstSel) { - AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask))); - continue; - } - else { - AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask))); - } - } + if (!c) continue; TLorentzVector nP; - c->GetMomentum(nP, vertex); + c->GetMomentum(nP, fVertex); cEta = nP.Eta(); cPhi = nP.Phi(); cPt = nP.Pt(); @@ -415,278 +316,83 @@ void AliEmcalJetTask::FindJets() cPy = nP.Py(); cPz = nP.Pz(); } - if (!c->IsEMCAL()) - continue; - if (cPt < fMinJetClusPt) - continue; + + if (!c->IsEMCAL()) continue; + if (cPt < fMinJetClusPt) continue; if ((cEtafEtaMax) || (cPhifPhiMax)) continue; + + Int_t lab = TMath::Abs(c->GetLabel()); + if (lab < fMinLabelClusters || lab > fMaxLabelClusters) { + AliDebug(2,Form("Skipping cluster %d because label %d is not in range (%d, %d)", iClus, lab, fMinLabelClusters, fMaxLabelClusters)); + continue; + } + // offset of 100 to skip ghost particles uid = -1 AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel())); Double_t e = TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz); - if(fPionMassClusters) e = TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz + 0.13957*0.13957); //MV: dirty, need better solution - fjw.AddInputVector(cPx, cPy, cPz, e, -iClus - 100); - // fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100); + fFastJetWrapper.AddInputVector(cPx, cPy, cPz, e, -iClus - 100); } } - // setting legacy mode - if (fLegacyMode) { - fjw.SetLegacyMode(kTRUE); - } - // run jet finder - fjw.Run(); - - //run generic subtractor - if(fDoGenericSubtractionJetMass) { - fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom); - fjw.DoGenericSubtractionJetMass(); - } - - if(fDoGenericSubtractionExtraJetShapes) { - fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom); - fjw.DoGenericSubtractionJetAngularity(); - fjw.DoGenericSubtractionJetpTD(); - fjw.DoGenericSubtractionJetCircularity(); - fjw.DoGenericSubtractionJetSigma2(); - fjw.DoGenericSubtractionJetConstituent(); - fjw.DoGenericSubtractionJetLeSub(); - } + fFastJetWrapper.Run(); +} - //run constituent subtractor - if(fDoConstituentSubtraction) { - fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom); - fjw.DoConstituentSubtraction(); - } +//________________________________________________________________________ +void AliEmcalJetTask::FillJetBranch() +{ + // Fill jet branch with fastjet jets. - // get geometry - AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); - if (!geom) { - AliFatal(Form("%s: Can not create geometry", GetName())); - return; - } + PrepareUtilities(); // loop over fastjet jets - std::vector jets_incl = fjw.GetInclusiveJets(); + std::vector jets_incl = fFastJetWrapper.GetInclusiveJets(); // sort jets according to jet pt static Int_t indexes[9999] = {-1}; GetSortedArray(indexes, jets_incl); AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size())); - for (UInt_t ijet=0, jetCount=0; ijetfJetEtaMax) || - (jets_incl[ij].phi()fJetPhiMax)) + if (jets_incl[ij].perp() < fMinJetPt) continue; + if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue; + if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) || + (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax)) continue; 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 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 num = fjw.GetGRNumerator(); - std::vector den = fjw.GetGRDenominator(); - std::vector nums = fjw.GetGRNumeratorSub(); - std::vector 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; gAddGRNumAt(num[g],g); - jet->AddGRNumSubAt(nums[g],g); - } - Int_t dsize = (Int_t)den.size(); - for(Int_t g = 0; gAddGRDenAt(den[g],g); - jet->AddGRDenSubAt(dens[g],g); - } - } - - if(fDoGenericSubtractionExtraJetShapes){ - std::vector jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity(); - Int_t na = (Int_t)jetAngularityInfo.size(); - if(na>ij && na>0) { - jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative()); - jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative()); - jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted()); - jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted()); - } - - std::vector jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD(); - Int_t np = (Int_t)jetpTDInfo.size(); - if(np>ij && np>0) { - jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative()); - jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative()); - jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted()); - jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted()); - } - - std::vector jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity(); - Int_t nc = (Int_t)jetCircularityInfo.size(); - if(nc>ij && nc>0) { - jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative()); - jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative()); - jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted()); - jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted()); - } - - std::vector jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2(); - Int_t ns = (Int_t)jetSigma2Info.size(); - if(ns>ij && ns>0) { - jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative()); - jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative()); - jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted()); - jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted()); - } - - - std::vector jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent(); - Int_t nco= (Int_t)jetConstituentInfo.size(); - if(nco>ij && nco>0) { - jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative()); - jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative()); - jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted()); - jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted()); - } - - std::vector jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub(); - Int_t nlsub= (Int_t)jetLeSubInfo.size(); - if(nlsub>ij && nlsub>0) { - jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative()); - jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative()); - jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted()); - jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted()); - } - } -#endif - - // Fill constituent info - std::vector constituents(fjw.GetJetConstituents(ij)); - jet->SetNumberOfTracks(constituents.size()); - jet->SetNumberOfClusters(constituents.size()); - - 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; - - FillJetConstituents(constituents,jet,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents,0); - jet->SetNumberOfTracks(nt); - jet->SetNumberOfClusters(nc); - jet->SortConstituents(); - jet->SetMaxNeutralPt(maxNe); - jet->SetMaxChargedPt(maxCh); - jet->SetNEF(neutralE / jet->E()); - fastjet::PseudoJet area(fjw.GetJetAreaVector(ij)); + fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij)); jet->SetArea(area.perp()); jet->SetAreaEta(area.eta()); jet->SetAreaPhi(area.phi()); jet->SetAreaE(area.E()); - 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); - - AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size())); + + // Fill constituent info + std::vector constituents(fFastJetWrapper.GetJetConstituents(ij)); + FillJetConstituents(jet, constituents, fTracks, fClus, constituents); + + if (fGeom) { + if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) && + (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) && + (jet->Eta() > fGeom->GetArm1EtaMin()) && + (jet->Eta() < fGeom->GetArm1EtaMax())) + jet->SetAxisInEmcal(kTRUE); + } + + ExecuteUtilities(jet, ij); + + AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents())); 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 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; ijetSetLabel(ijet); - // Fill constituent info - std::vector constituents_unsub(fjw.GetJetConstituents(ijet)); - std::vector constituents_sub = jets_sub[ijet].constituents(); - jet_sub->SetNumberOfTracks(constituents_sub.size()); - jet_sub->SetNumberOfClusters(constituents_sub.size()); - 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; - - FillJetConstituents(constituents_sub,jet_sub,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents_unsub,1); - jet_sub->SetNumberOfTracks(nt); - jet_sub->SetNumberOfClusters(nc); - jet_sub->SortConstituents(); - - 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 + TerminateUtilities(); } //________________________________________________________________________ @@ -719,6 +425,12 @@ Bool_t AliEmcalJetTask::DoInit() gRandom = new TRandom3(0); } + // get geometry + fGeom = AliEMCALGeometry::GetInstance(); + if (!fGeom) { + AliWarning(Form("%s: Can not create EMCal geometry, some features will not work...", GetName())); + } + // get input collections AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); @@ -738,15 +450,6 @@ Bool_t AliEmcalJetTask::DoInit() return 0; } } - if (fTracks) { - TClass cls(fTracks->GetClass()->GetName()); - if (cls.InheritsFrom("AliPicoTrack")) - fIsPicoTrack = 1; - else if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle")) - fIsMcPart = 1; - else if (cls.InheritsFrom("AliEmcalParticle")) - fIsEmcPart = 1; - } if (fCaloName == "CaloClusters") am->LoadBranch("CaloClusters"); @@ -763,249 +466,225 @@ Bool_t AliEmcalJetTask::DoInit() fIsEmcPart = 1; } - //Add constituent subtracted jets to event - if(!fJetsSubName.IsNull()) { - if (!(fEvent->FindListObject(fJetsSubName)) ) { - if(fJetsSub) fEvent->AddObject(fJetsSub); - } else { - AliError(Form("%s: Object for subtracted jet branch with name %s already in event! Returning", GetName(), fJetsSubName.Data())); - return 0; - } + // add jets to event if not yet there + if (!(fEvent->FindListObject(fJetsName))) { + fJets = new TClonesArray("AliEmcalJet"); + fJets->SetName(fJetsName); + fEvent->AddObject(fJets); } - - //Add tracks from constituent subtracted jets to event - if(!fTracksSubName.IsNull()) { - if (!(fEvent->FindListObject(fTracksSubName))) { - if(fTracksSub) fEvent->AddObject(fTracksSub); - } else { - AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fTracksSubName.Data())); - return 0; - } + else { + AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data())); + return 0; } - if(fDoConstituentSubtraction && fTracksSub) { - if(fIsEmcPart==1) { - AliFatal("Constituent subtraction method not implemented for AliEmcalParticle as input tracks. Use AliPicoTracks instead. Aborting!"); - return 0; - } - if(fIsPicoTrack==1) fTracksSub->SetClass("AliPicoTrack"); - else if(fIsMcPart==1) fTracksSub->SetClass("AliMCParticle"); - else { - AliError(Form("%s Object type of subtracted track branch %s not known",GetName(),fTracksSubName.Data())); - return 0; - } - fTracksSub->SetName(fTracksSubName); + TString name("kt"); + fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm); + if ((fJetType & kAKT) != 0) { + name = "antikt"; + jalgo = fastjet::antikt_algorithm; + AliDebug(1,"Using AKT algorithm"); } - - //Add clusters from constituent subtracted jets to event - if(!fClusSubName.IsNull()) { - if (!(fEvent->FindListObject(fClusSubName))) { - if(fClusSub) fEvent->AddObject(fClusSub); - } else { - AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fClusSubName.Data())); - return 0; - } + else { + AliDebug(1,"Using KT algorithm"); } - if(fDoConstituentSubtraction && fClusSub) { - TClass cls(fClus->GetClass()->GetName()); - if (cls.InheritsFrom("AliESDCaloCluster")) { - fClusSub->SetClass("AliESDCaloCluster"); - fIsEsdClus = kTRUE; - } - else if (cls.InheritsFrom("AliAODCaloCluster")) - fClusSub->SetClass("AliAODCaloCluster"); - else if (cls.InheritsFrom("AliEmcalParticle")) { - AliFatal("Constituent subtraction method not implemented for AliEmcalParticle as input clusters. Use AliESD(AOD)CaloCluster instead. Aborting!"); - return 0; - } else { - AliError(Form("%s Object type of subtracted cluster branch %s not known",GetName(),fClusSubName.Data())); - return 0; - } - fClusSub->SetName(fClusSubName); - } + if ((fJetType & kR020Jet) != 0) + fRadius = 0.2; + else if ((fJetType & kR030Jet) != 0) + fRadius = 0.3; + else if ((fJetType & kR040Jet) != 0) + fRadius = 0.4; - if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event - fRhoParam = dynamic_cast(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(InputEvent()->FindListObject(fRhomName)); - if (!fRhomParam) { - AliError(Form("%s: Could not retrieve rho_m %s!", GetName(), fRhomName.Data())); - return 0; - } - } + // setup fj wrapper + fFastJetWrapper.SetName(name); + fFastJetWrapper.SetTitle(name); + fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts); + fFastJetWrapper.SetGhostArea(fGhostArea); + fFastJetWrapper.SetR(fRadius); + fFastJetWrapper.SetAlgorithm(jalgo); + fFastJetWrapper.SetRecombScheme(static_cast(fRecombScheme)); + fFastJetWrapper.SetMaxRap(fEtaMax); - // add jets to event if not yet there - if (!(fEvent->FindListObject(fJetsName))) - fEvent->AddObject(fJets); - else { - AliError(Form("%s: Object for jet branch with name %s already in event! Returning", GetName(), fJetsName.Data())); - return 0; + // setting legacy mode + if (fLegacyMode) { + fFastJetWrapper.SetLegacyMode(kTRUE); } - return 1; + InitUtilities(); + + return kTRUE; } //___________________________________________________________________________________________________________________ -void AliEmcalJetTask::FillJetConstituents(std::vector& constituents,AliEmcalJet *jet,Double_t vertex[3],UInt_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc,std::vector& constituentsunsub,Int_t flagsub){ - nt = 0; - nc = 0; - neutralE = 0; - maxCh = 0; - maxNe = 0; - gall = 0; - gemc = 0; - cemc = 0; - ncharged = 0; - nneutral = 0; - mcpt = 0; - emcpt = 0; - Int_t uid = -1; - AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(); - for(UInt_t ic = 0; ic < constituents.size(); ++ic) { - if(flagsub==0) uid = constituents[ic].user_index(); - else uid = GetIndexSub(constituents[ic].phi(),constituentsunsub); - AliDebug(3,Form("Processing constituent %d", uid)); - if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle - ++gall; +void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector& constituents, TClonesArray *tracks, TClonesArray *clusters, + std::vector& constituents_unsub, Int_t flag, TClonesArray *particles_sub) +{ + 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.; + + Int_t uid = -1; + + jet->SetNumberOfTracks(constituents.size()); + jet->SetNumberOfClusters(constituents.size()); + + for (UInt_t ic = 0; ic < constituents.size(); ++ic) { + + if (flag == 0) uid = constituents[ic].user_index(); + else uid = GetIndexSub(constituents[ic].phi(), constituents_unsub); + + AliDebug(3,Form("Processing constituent %d", uid)); + if (uid == -1) { //ghost particle + ++gall; + if (fGeom) { Double_t gphi = constituents[ic].phi(); - if (gphi<0) - gphi += TMath::TwoPi(); + 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())) + if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) && + (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax())) ++gemc; - } else if ((uid > 0) && fTracks) { // track constituent - Int_t tid = uid - 100; - AliVParticle *t = static_cast(fTracks->At(tid)); - if (!t) { - AliError(Form("Could not find track %d",tid)); - continue; - } - if (jetCount < fMarkConst) { - AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType)); - t->SetBit(fJetType); - } - 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; - } - if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle - mcpt += cPt; + } + } + else if ((uid >= 100) && tracks) { // track constituent + Int_t tid = uid - 100; + AliVParticle *t = static_cast(tracks->At(tid)); + if (!t) { + AliError(Form("Could not find track %d",tid)); + continue; + } - if (cPhi<0) - cPhi += TMath::TwoPi(); + 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; + } + + // check if MC particle + if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt; + + if (fGeom) { + if (cPhi < 0) cPhi += TMath::TwoPi(); cPhi *= TMath::RadToDeg(); - if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) && - (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) { + if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) && + (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) { emcpt += cPt; ++cemc; } + } - if(flagsub!=0 && fTracksSub) { //fill the subtracted track branch - if(fIsPicoTrack) - new ((*fTracksSub)[tid]) AliPicoTrack(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),0,0,0,0,0,0,0,constituents[ic].m()); - else if(fIsMcPart) { - AliMCParticle *tmc = static_cast(fTracks->At(tid)); - TParticle *mcPart = new TParticle(tmc->PdgCode(),tmc->Particle()->GetStatusCode(),tmc->GetMother(),0,tmc->GetFirstDaughter(),tmc->GetLastDaughter(),constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),constituents[ic].e(),0.,0.,0.,0.); - new ((*fTracksSub)[tid]) AliMCParticle(mcPart,0x0,tid); - } - } + if (flag == 0 || particles_sub == 0) { 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(fClus->At(cid)); - if (!ep) - continue; - c = ep->GetCluster(); - if (!c) - continue; - if (jetCount < fMarkConst) - ep->SetBit(fJetType); - cEta = ep->Eta(); - cPhi = ep->Phi(); - cPt = ep->Pt(); - cP = ep->P(); - } else { - c = static_cast(fClus->At(cid)); - if (!c) - continue; - if (jetCount < fMarkConst) - c->SetBit(fJetType); - TLorentzVector nP; - c->GetMomentum(nP, vertex); - cEta = nP.Eta(); - cPhi = nP.Phi(); - cPt = nP.Pt(); - cP = nP.P(); - } + } + else { + Int_t part_sub_id = particles_sub->GetEntriesFast(); + AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(t); + part_sub->SetPt(constituents[ic].perp()); + jet->AddTrackAt(part_sub_id, nt); + } - neutralE += cP; - if (cPt > maxNe) - maxNe = cPt; + ++nt; + } + else if ((uid <= -100) && clusters) { // 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(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(fClus->At(cid)); + if (!c) continue; + + TLorentzVector nP; + c->GetMomentum(nP, fVertex); + cEta = nP.Eta(); + cPhi = nP.Phi(); + cPt = nP.Pt(); + cP = nP.P(); + } + + neutralE += cP; + if (cPt > maxNe) maxNe = cPt; - if (c->GetLabel() > fMinMCLabel) // MC particle - mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt; + // MC particle + if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt; - if (cPhi<0) - cPhi += TMath::TwoPi(); + if (fGeom) { + if (cPhi<0) cPhi += TMath::TwoPi(); cPhi *= TMath::RadToDeg(); - if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) && - (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) { + if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) && + (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) { emcpt += cPt; ++cemc; } + } - //Add constituent subtracted clusters - if(flagsub!=0 && fClusSub && !fIsEmcPart) { //fill the subtracted track branch - AliVCluster *oc; - if(fIsEsdClus) { - AliESDCaloCluster *ec = dynamic_cast(c); - if (!ec) continue; - oc = new ((*fClusSub)[cid]) AliESDCaloCluster(*ec); - } else { - AliAODCaloCluster *ac = dynamic_cast(c); - if (!ac) continue; - oc = new ((*fClusSub)[cid]) AliAODCaloCluster(*ac); - } - oc->SetE(constituents[ic].E()); - } - + if (flag == 0 || particles_sub == 0) { jet->AddClusterAt(cid, nc); - ++nc; - ++nneutral; - } else { - AliError(Form("%s: No logical way to end up here.", GetName())); - continue; } + else { + Int_t part_sub_id = particles_sub->GetEntriesFast(); + AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c); + part_sub->SetPt(constituents[ic].perp()); + jet->AddTrackAt(part_sub_id, nt); + } + + ++nc; + ++nneutral; + } + else { + AliError(Form("%s: No logical way to end up here.", GetName())); + continue; } + } + + jet->SetNumberOfTracks(nt); + jet->SetNumberOfClusters(nc); + jet->SetNEF(neutralE / jet->E()); + jet->SetMaxChargedPt(maxCh); + jet->SetMaxNeutralPt(maxNe); + if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall); + else jet->SetAreaEmc(-1); + jet->SetNEmc(cemc); + jet->SetNumberOfCharged(ncharged); + jet->SetNumberOfNeutrals(nneutral); + jet->SetMCPt(mcpt); + jet->SetPtEmc(emcpt); + jet->SortConstituents(); } + //______________________________________________________________________________________ -Int_t AliEmcalJetTask::GetIndexSub(Double_t phisub,std::vector& constituentsunsub) { - for(UInt_t ii=0;ii& constituents_unsub) +{ + for (UInt_t ii = 0; ii < constituents_unsub.size(); ii++) { + Double_t phi_unsub = constituents_unsub[ii].phi(); + if (phi_sub == phi_unsub) return constituents_unsub[ii].user_index(); + } + return -1; } -//______________________________________________________________________________________ diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.h b/PWGJE/EMCALJetTasks/AliEmcalJetTask.h index 69b2a58de42..56e5fd0d0c3 100644 --- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.h +++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.h @@ -2,8 +2,11 @@ #define ALIEMCALJETTASK_H class TClonesArray; +class TObjArray; class AliVEvent; -class AliRhoParameter; +class AliEMCALGeometry; +class AliEmcalJetUtility; +class AliFJWrapper; namespace fastjet { class PseudoJet; @@ -37,40 +40,83 @@ class AliEmcalJetTask : public AliAnalysisTaskSE { AliEmcalJetTask(const char *name); virtual ~AliEmcalJetTask(); - void UserCreateOutputObjects(); - void UserExec(Option_t *option); - void Terminate(Option_t *option); - - Bool_t IsLocked() { if(fLocked) {AliFatal("Jet finder task is locked! Changing properties is not allowed."); return kTRUE;} else return kFALSE;}; - void SetLocked() { fLocked = kTRUE;} - void SelectConstituents(UInt_t constSel, UInt_t MCconstSel) { fConstSel = constSel; fMCConstSel = MCconstSel; }; - void SetAlgo(Int_t a) { if(IsLocked()) return; if (a==0) fJetType |= kKT; else fJetType |= kAKT; } // for backward compatibility only - void SetClusName(const char *n) { if(IsLocked()) return; fCaloName = n ; } - void SetJetsName(const char *n) { if(IsLocked()) return; fJetsName = n ; } - void SetJetsSubName(const char *n) { fJetsSubName = n ; } - void SetTracksSubName(const char *n) { fTracksSubName = n ; } - void SetClusSubName(const char *n) { fClusSubName = n ; } - void SetJetType(UInt_t t) { if(IsLocked()) return; fJetType = t ; } - void SetMarkConstituents(UInt_t m) { if(IsLocked()) return; fMarkConst = m ; } - void SetMinJetArea(Double_t a) { if(IsLocked()) return; fMinJetArea = a ; } - void SetMinJetClusPt(Double_t min) { if(IsLocked()) return; fMinJetClusPt = min ; } - void SetMinJetPt(Double_t j) { if(IsLocked()) return; fMinJetPt = j ; } - void SetMinJetTrackPt(Double_t min) { if(IsLocked()) return; fMinJetTrackPt = min ; } - void SetRadius(Double_t r) { if(IsLocked()) return; fRadius = r ; if ((fJetType & (kRX1Jet|kRX2Jet|kRX3Jet)) == 0) AliWarning("Radius value will be ignored if jet type is not set to a user defined radius (kRX1Jet,kRX2Jet,kRX3Jet)."); } - void SetTrackEfficiency(Double_t t) { if(IsLocked()) return; fTrackEfficiency = t ; } - void SetTracksName(const char *n) { if(IsLocked()) return; fTracksName = n ; } - void SetType(Int_t t) { if(IsLocked()) return; - if (t==0) fJetType |= kFullJet; - else if (t==1) fJetType |= kChargedJet; - else if (t==2) fJetType |= kNeutralJet; } // for backward compatibility only - void SetEtaRange(Double_t emi, Double_t ema) {if(IsLocked()) return; fEtaMin = emi; fEtaMax = ema; } - void SetPhiRange(Double_t pmi, Double_t pma) {if(IsLocked()) return; fPhiMin = pmi; fPhiMax = pma; } - void SetJetEtaRange(Double_t emi, Double_t ema) {if(IsLocked()) return; fJetEtaMin = emi; fJetEtaMax = ema; } - void SetJetPhiRange(Double_t pmi, Double_t pma) {if(IsLocked()) return; fJetPhiMin = pmi; fJetPhiMax = pma; } - void SetGhostArea(Double_t gharea) { if(IsLocked()) return; fGhostArea = gharea; } - void SetMinMCLabel(Int_t s) { if(IsLocked()) return; fMinMCLabel = s ; } - void SetRecombScheme(Int_t scheme) { if(IsLocked()) return; fRecombScheme = scheme; } - void SelectCollisionCandidates(UInt_t offlineTriggerMask = AliVEvent::kMB) + void UserCreateOutputObjects(); + void UserExec(Option_t* option); + void Terminate(Option_t* /*option*/) {;} + + void SetAlgo(Int_t a) { if (IsLocked()) return; if (a==0) fJetType |= kKT; else fJetType |= kAKT; } // for backward compatibility only + void SetClusName(const char *n) { if (IsLocked()) return; fCaloName = n ; } + void SetClusLabelRange(Int_t min, Int_t max) { if (IsLocked()) return; fMinLabelClusters = min ; fMaxLabelClusters = max; } + void SetEtaRange(Double_t emi, Double_t ema) { if (IsLocked()) return; fEtaMin = emi ; fEtaMax = ema ; } + void SetGhostArea(Double_t gharea) { if (IsLocked()) return; fGhostArea = gharea; } + void SetJetsName(const char *n) { if (IsLocked()) return; fJetsName = n ; } + void SetJetEtaRange(Double_t emi, Double_t ema) { if (IsLocked()) return; fJetEtaMin = emi ; fJetEtaMax = ema; } + void SetJetPhiRange(Double_t pmi, Double_t pma) { if (IsLocked()) return; fJetPhiMin = pmi ; fJetPhiMax = pma; } + void SetJetType(UInt_t t) { if (IsLocked()) return; fJetType = t ; } + void SetLocked() { fLocked = kTRUE;} + void SetMinJetArea(Double_t a) { if (IsLocked()) return; fMinJetArea = a ; } + void SetMinJetClusPt(Double_t min) { if (IsLocked()) return; fMinJetClusPt = min ; } + void SetMinJetPt(Double_t j) { if (IsLocked()) return; fMinJetPt = j ; } + void SetMinJetTrackPt(Double_t min) { if (IsLocked()) return; fMinJetTrackPt = min ; } + void SetMinMCLabel(Int_t s) { if (IsLocked()) return; fMinMCLabel = s ; } + void SetPhiRange(Double_t pmi, Double_t pma) { if (IsLocked()) return; fPhiMin = pmi ; fPhiMax = pma; } + void SetRecombScheme(Int_t scheme) { if (IsLocked()) return; fRecombScheme = scheme; } + void SetTracksName(const char *n) { if (IsLocked()) return; fTracksName = n ; } + void SetTrackEfficiency(Double_t t) { if (IsLocked()) return; fTrackEfficiency = t ; } + void SetTrackLabelRange(Int_t min, Int_t max) { if (IsLocked()) return; fMinLabelTracks = min ; fMaxLabelTracks = max; } + void SetLegacyMode(Bool_t mode) { if (IsLocked()) return; fLegacyMode = mode ; } + void SetMCFlag(UInt_t m) { if (IsLocked()) return; fMCFlag = m ; } + void SelectHIJING(Bool_t s) { if (IsLocked()) return; if (s) fGeneratorIndex = 0; else fGeneratorIndex = -1; } + void SetGeneratorIndex(Short_t i) { if (IsLocked()) return; fGeneratorIndex = i ; } + + AliEmcalJetUtility* AddUtility(AliEmcalJetUtility* utility); + + const char* GetClusName() { return fCaloName.Data() ; } + Double_t GetEtaMin() { return fEtaMin ; } + Double_t GetEtaMax() { return fEtaMax ; } + Double_t GetGhostArea() { return fGhostArea ; } + const char* GetJetsName() { return fJetsName.Data() ; } + Double_t GetJetEtaMin() { return fJetEtaMin ; } + Double_t GetJetEtaMax() { return fJetEtaMax ; } + Double_t GetJetPhiMin() { return fJetPhiMin ; } + Double_t GetJetPhiMax() { return fJetPhiMax ; } + UInt_t GetJetType() { return fJetType ; } + Bool_t GetLegacyMode() { return fLegacyMode ; } + Double_t GetMinJetArea() { return fMinJetArea ; } + Double_t GetMinJetClusPt() { return fMinJetClusPt ; } + Double_t GetMinJetPt() { return fMinJetPt ; } + Double_t GetMinJetTrackPt() { return fMinJetTrackPt ; } + Int_t GetMinMCLabel() { return fMinMCLabel ; } + Double_t GetPhiMin() { return fPhiMin ; } + Double_t GetPhiMax() { return fPhiMax ; } + Double_t GetRadius() { return fRadius ; } + Int_t GetRecombScheme() { return fRecombScheme ; } + const char* GetTracksName() { return fTracksName.Data() ; } + Double_t GetTrackEfficiency() { return fTrackEfficiency ; } + + AliVEvent* GetEvent() { return fEvent ; } + TClonesArray* GetClusters() { return fClus ; } + TClonesArray* GetTracks() { return fTracks ; } + TClonesArray* GetJets() { return fJets ; } + TObjArray* GetUtilities() { return fUtilities ; } + + void FillJetConstituents(AliEmcalJet *jet, std::vector& constituents, TClonesArray *tracks, TClonesArray *clusters, + std::vector& constituents_sub, Int_t flag = 0, TClonesArray *particles_sub = 0); + + Int_t GetIndexSub(Double_t phi_sub, std::vector& constituents_unsub); + + Bool_t IsLocked() + { + if (fLocked) { + AliFatal("Jet finder task is locked! Changing properties is not allowed."); + return kTRUE; + } + else { + return kFALSE; + } + } + + void SelectCollisionCandidates(UInt_t offlineTriggerMask = AliVEvent::kMB) { if(!fIsPSelSet) { @@ -83,78 +129,52 @@ class AliEmcalJetTask : public AliAnalysisTaskSE { fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask; } } - void SetLegacyMode(Bool_t mode) { if(IsLocked()) return; fLegacyMode = mode; } - void SetMCFlag(UInt_t m) { if(IsLocked()) return; fMCFlag = m ; } - void SelectHIJING(Bool_t s) { if(IsLocked()) return; if (s) fGeneratorIndex = 0; else fGeneratorIndex = -1; } - void SetGeneratorIndex(Short_t i) { if(IsLocked()) return; fGeneratorIndex = i; } - void SelectPhysicalPrimaries(Bool_t s) { if(IsLocked()) return; - if (s) fMCFlag |= AliAODMCParticle::kPhysicalPrim ; - else fMCFlag &= ~AliAODMCParticle::kPhysicalPrim ; } - - void SetCodeDebug(Bool_t val) { fCodeDebug = val ; } - void SetForceIsMcPart(Bool_t b) { fIsMcPart = b ; } - void SetPionMassForClusters(Bool_t b) { fPionMassClusters = b ; } - - void SetRhoName(const char *n) { fRhoName = n ; } - void SetRhomName(const char *n) { fRhomName = n ; } - void SetGenericSubtractionJetMass(Bool_t b) { fDoGenericSubtractionJetMass = b; } - void SetGenericSubtractionGR(Bool_t b, Double_t rmax = 2., Double_t dr = 0.04, Double_t ptmin = 0.) { fDoGenericSubtractionGR = b; fRMax=rmax; fDRStep=dr; fPtMinGR=ptmin;} - void SetConstituentSubtraction(Bool_t b) { fDoConstituentSubtraction = b; } - void SetGenericSubtractionExtraJetShapes(Bool_t b) { fDoGenericSubtractionExtraJetShapes =b;} - void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;} - - UInt_t GetJetType() { return fJetType; } - Bool_t GetLegacyMode() { return fLegacyMode; } - TString GetJetsSubName() { return fJetsSubName; } - - const char* GetJetsName() { return fJetsName.Data(); } - Double_t GetRadius() { return fRadius; } - const char* GetTracksName() { return fTracksName.Data(); } - const char* GetClusName() { return fCaloName.Data(); } - const char* GetRhoName() { return fRhoName.Data(); } - UInt_t GetMarkConstituents() { return fMarkConst; } - Double_t GetMinJetArea() { return fMinJetArea; } - Double_t GetMinJetClusPt() { return fMinJetClusPt; } - Double_t GetMinJetPt() { return fMinJetPt; } - Double_t GetMinJetTrackPt() { return fMinJetTrackPt; } - Double_t GetTrackEfficiency() { return fTrackEfficiency; } - Double_t GetGhostArea() { return fGhostArea; } - Int_t GetMinMCLabel() { return fMinMCLabel; } - Int_t GetRecombScheme() { return fRecombScheme; } - Double_t GetEtaMin() { return fEtaMin; } - Double_t GetEtaMax() { return fEtaMax; } - Double_t GetPhiMin() { return fPhiMin; } - Double_t GetPhiMax() { return fPhiMax; } - Double_t GetJetEtaMin() { return fJetEtaMin; } - Double_t GetJetEtaMax() { return fJetEtaMax; } - Double_t GetJetPhiMin() { return fJetPhiMin; } - Double_t GetJetPhiMax() { return fJetPhiMax; } - - AliVEvent* GetEvent() { return fEvent;} - TClonesArray* GetClusters() { return fClus;} - TClonesArray* GetTracks() { return fTracks;} - TClonesArray* GetJets() { return fJets;} + + void SetRadius(Double_t r) + { + if (IsLocked()) return; + fRadius = r; + if ((fJetType & (kRX1Jet|kRX2Jet|kRX3Jet)) == 0) + AliWarning("Radius value will be ignored if jet type is not set to a user defined radius (kRX1Jet,kRX2Jet,kRX3Jet)."); + } + + void SetType(Int_t t) + { + if(IsLocked()) return; + if (t==0) fJetType |= kFullJet; + else if (t==1) fJetType |= kChargedJet; + else if (t==2) fJetType |= kNeutralJet; + } // for backward compatibility only + + void SelectPhysicalPrimaries(Bool_t s) + { + if(IsLocked()) return; + if (s) fMCFlag |= AliAODMCParticle::kPhysicalPrim; + else fMCFlag &= ~AliAODMCParticle::kPhysicalPrim; + } protected: + void FindJets(); + void FillJetBranch(); Bool_t DoInit(); + void InitUtilities(); + void PrepareUtilities(); + void ExecuteUtilities(AliEmcalJet* jet, Int_t ij); + void TerminateUtilities(); Bool_t GetSortedArray(Int_t indexes[], std::vector array) const; - void FillJetConstituents(std::vector& constituents,AliEmcalJet *jet,Double_t vertex[3],UInt_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc,std::vector& constituentsunsub,Int_t flag); - Int_t GetIndexSub(Double_t phisub,std::vector& constituentsunsub); - TString fTracksName; // name of track collection TString fCaloName; // name of calo cluster collection TString fJetsName; // name of jet collection - TString fJetsSubName; // name of subtracted jet collection - TString fTracksSubName; // name of subtracted track collection - TString fClusSubName; // name of subtracted clusters collection UInt_t fJetType; // jet type (algorithm, radius, constituents) - UInt_t fConstSel; // select constituents from a previous jet finding - UInt_t fMCConstSel; // select MC constituents (label!=0) from a previous jet finding - UInt_t fMarkConst; // constituents are marked (via TObject::SetBit) as belonging to the # leading jets + Int_t fMinLabelTracks; // select track constituents with a minimum label index (track->GetLabel()) + Int_t fMaxLabelTracks; // select track constituents with a maximum label index (track->GetLabel()) + Int_t fMinLabelClusters; // select cluster constituents with a minimum label index (cluster->GetLabel()) + Int_t fMaxLabelClusters; // select cluster constituents with a maximum label index (cluster->GetLabel()) + Int_t fMinMCLabel; // minimum MC label value for the tracks/clusters being considered MC particles Double_t fRadius; // jet radius - Double_t fMinJetTrackPt; // min jet track momentum (applied before clustering) + Double_t fMinJetTrackPt; // min jet track momentum (applied before clustering) Double_t fMinJetClusPt; // min jet cluster momentum (applied before clustering) Double_t fPhiMin; // minimum phi for constituents (applied before clustering) Double_t fPhiMax; // maximum phi for constituents (applied before clustering) @@ -167,44 +187,25 @@ class AliEmcalJetTask : public AliAnalysisTaskSE { Double_t fJetEtaMin; // minimum eta to keep jet in output Double_t fJetEtaMax; // maximum eta to keep jet in output Double_t fGhostArea; // ghost area - Int_t fMinMCLabel; // minimum MC label value for the tracks/clusters being considered MC particles Int_t fRecombScheme; // recombination scheme used by fastjet Double_t fTrackEfficiency; // artificial tracking inefficiency (0...1) - UInt_t fMCFlag; // select MC particles with flags - Short_t fGeneratorIndex; // select MC particles with generator index (default = -1 = switch off selection) - Bool_t fIsInit; //!=true if already initialized + UInt_t fMCFlag; // select MC particles with flags (e.g. to select physical primaries) + Short_t fGeneratorIndex; // select MC particles with generator index (default = -1 to switch off selection) + TObjArray *fUtilities; // jet utilities (gen subtractor, constituent subtractor etc.) Bool_t fLocked; // true if lock is set + + Bool_t fIsInit; //!=true if already initialized Bool_t fIsPSelSet; //!=true if physics selection was set - Bool_t fIsMcPart; //!=true if MC particles are given as input Bool_t fIsEmcPart; //!=true if emcal particles are given as input (for clusters) - Bool_t fIsPicoTrack; //!=true if pico tracks are given as input - Bool_t fIsEsdClus; //!=true if AliESDCaloCluster is given as input - Bool_t fLegacyMode; //! if true, enable FJ 2.x behavior - Bool_t fCodeDebug; // use nontested code changes - Bool_t fPionMassClusters; // assume pion mass for clusters - Bool_t fDoGenericSubtractionJetMass; // calculate generic subtraction - Bool_t fDoGenericSubtractionGR; // calculate generic subtraction for angular structure function GR - Bool_t fDoGenericSubtractionExtraJetShapes; // calculate generic subtraction for other jet shapes like radialmoment,pTD etc - Bool_t fDoConstituentSubtraction; // calculate constituent subtraction - Bool_t fUseExternalBkg; // use external background for generic subtractor - TString fRhoName; // name of rho - TString fRhomName; // name of rhom - Double_t fRho; // pT background density - Double_t fRhom; // mT background density - Double_t fRMax; // R max for GR calculation - Double_t fDRStep; // step width for GR calculation - Double_t fPtMinGR; // min pT for GR calculation - + Bool_t fLegacyMode; //!=true to enable FJ 2.x behavior + Double_t fVertex[3]; //!vertex of the current event + AliEMCALGeometry *fGeom; //!emcal geometry object TClonesArray *fJets; //!jet collection - TClonesArray *fJetsSub; //!subtracted jet collection - TClonesArray *fTracksSub; //!subtracted track collection - TClonesArray *fClusSub; //!subtracted cluster collection AliVEvent *fEvent; //!current event TClonesArray *fTracks; //!tracks collection TClonesArray *fClus; //!cluster collection - AliRhoParameter *fRhoParam; //!event rho - AliRhoParameter *fRhomParam; //!event rhom + AliFJWrapper fFastJetWrapper; //!fastjet wrapper private: AliEmcalJetTask(const AliEmcalJetTask&); // not implemented diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx new file mode 100644 index 00000000000..35e53e52880 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx @@ -0,0 +1,42 @@ +#include "AliEmcalJetUtility.h" + +ClassImp(AliEmcalJetUtility) + +//______________________________________________________________________________ +AliEmcalJetUtility::AliEmcalJetUtility() : +TNamed(), + fJetTask(0), + fInit(kFALSE) +{ + // Dummy constructor. + +} + +//______________________________________________________________________________ +AliEmcalJetUtility::AliEmcalJetUtility(const char* name) : + TNamed(name, name), + fJetTask(0), + fInit(kFALSE) +{ + // Default constructor. +} + +//______________________________________________________________________________ +AliEmcalJetUtility::AliEmcalJetUtility(const AliEmcalJetUtility &other) : + TNamed(other), + fJetTask(other.fJetTask), + fInit(other.fInit) +{ + // Copy constructor. +} + +//______________________________________________________________________________ +AliEmcalJetUtility& AliEmcalJetUtility::operator=(const AliEmcalJetUtility &other) +{ + // Assignment + if (&other == this) return *this; + TNamed::operator=(other); + fJetTask = other.fJetTask; + fInit = other.fInit; + return *this; +} diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtility.h b/PWGJE/EMCALJetTasks/AliEmcalJetUtility.h new file mode 100644 index 00000000000..3283afd1c96 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliEmcalJetUtility.h @@ -0,0 +1,36 @@ +#ifndef ALIEMCALJETUTILITY_H +#define ALIEMCALJETUTILITY_H + +#include + +#include "AliFJWrapper.h" + +class AliEmcalJetTask; +class AliEmcalJet; +class AliFJWrapper; + +class AliEmcalJetUtility : public TNamed +{ + public: + + AliEmcalJetUtility(); + AliEmcalJetUtility(const char* name); + AliEmcalJetUtility(const AliEmcalJetUtility &jet); + AliEmcalJetUtility& operator=(const AliEmcalJetUtility &jet); + ~AliEmcalJetUtility() {;} + + virtual void Init() = 0; // Executed only once at the end of AliEmcalJetTask::DoInit() + virtual void Prepare(AliFJWrapper& fjw) = 0; // Executed for each event at the beginning of AliEmcalJetTask::FillJetBranch() + virtual void ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw) = 0; // Executed for each jet in the loop in AliEmcalJetTask::FillJetBranch() + virtual void Terminate(AliFJWrapper& fjw) = 0; // Executed for each event at the end of AliEmcalJetTask::FillJetBranch() + + void SetJetTask(AliEmcalJetTask* jetTask) { fJetTask = jetTask; } + + protected: + + AliEmcalJetTask *fJetTask ; // pointer to the main jet task + Bool_t fInit ; //! whether or not the utility has been initialized + + ClassDef(AliEmcalJetUtility, 1) // Abstract Emcal jet utility class +}; +#endif diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx new file mode 100644 index 00000000000..e46568ffdf0 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx @@ -0,0 +1,195 @@ +#include "AliEmcalJetUtilityConstSubtractor.h" +#include "AliEmcalJet.h" +#include "AliRhoParameter.h" +#include "AliEmcalJetTask.h" + +ClassImp(AliEmcalJetUtilityConstSubtractor) + +//______________________________________________________________________________ +AliEmcalJetUtilityConstSubtractor::AliEmcalJetUtilityConstSubtractor() : + AliEmcalJetUtility(), + fJetsSubName(""), + fParticlesSubName(""), + fUseExternalBkg(kFALSE), + fRhoName(""), + fRhomName(""), + fRho(0), + fRhom(0), + fJetsSub(0x0), + fParticlesSub(0x0), + fRhoParam(0), + fRhomParam(0) +{ + // Dummy constructor. + +} + +//______________________________________________________________________________ +AliEmcalJetUtilityConstSubtractor::AliEmcalJetUtilityConstSubtractor(const char* name) : + AliEmcalJetUtility(name), + fJetsSubName(""), + fParticlesSubName(""), + fUseExternalBkg(kFALSE), + fRhoName(""), + fRhomName(""), + fRho(0), + fRhom(0), + fJetsSub(0x0), + fParticlesSub(0x0), + fRhoParam(0), + fRhomParam(0) +{ + // Default constructor. +} + +//______________________________________________________________________________ +AliEmcalJetUtilityConstSubtractor::AliEmcalJetUtilityConstSubtractor(const AliEmcalJetUtilityConstSubtractor &other) : + AliEmcalJetUtility(other), + fJetsSubName(other.fJetsSubName), + fParticlesSubName(other.fParticlesSubName), + fUseExternalBkg(other.fUseExternalBkg), + fRhoName(other.fRhoName), + fRhomName(other.fRhomName), + fRho(other.fRho), + fRhom(other.fRhom), + fJetsSub(other.fJetsSub), + fParticlesSub(other.fParticlesSub), + fRhoParam(other.fRhoParam), + fRhomParam(other.fRhomParam) +{ + // Copy constructor. +} + +//______________________________________________________________________________ +AliEmcalJetUtilityConstSubtractor& AliEmcalJetUtilityConstSubtractor::operator=(const AliEmcalJetUtilityConstSubtractor &other) +{ + // Assignment. + + if (&other == this) return *this; + AliEmcalJetUtility::operator=(other); + fJetsSubName = other.fJetsSubName; + fParticlesSubName = other.fParticlesSubName; + fUseExternalBkg = other.fUseExternalBkg; + fRhoName = other.fRhoName; + fRhomName = other.fRhomName; + fRho = other.fRho; + fRhom = other.fRhom; + fJetsSub = other.fJetsSub; + fParticlesSub = other.fParticlesSub; + fRhoParam = other.fRhoParam; + fRhomParam = other.fRhomParam; + return *this; +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityConstSubtractor::Init() +{ + // Initialize the utility. + + // Add constituent subtracted jets to event + if (!fJetsSubName.IsNull()) { + if (!(fJetTask->GetEvent()->FindListObject(fJetsSubName)) ) { + fJetsSub = new TClonesArray("AliEmcalJet"); + fJetsSub->SetName(fJetsSubName); + fJetTask->GetEvent()->AddObject(fJetsSub); + } + else { + AliError(Form("%s: Object for subtracted jet branch with name %s already in event! Returning", GetName(), fJetsSubName.Data())); + return; + } + } + + // Add tracks from constituent subtracted jets to event + if (!fParticlesSubName.IsNull()) { + if (!(fJetTask->GetEvent()->FindListObject(fParticlesSubName))) { + fParticlesSub = new TClonesArray("AliEmcalParticle"); + fParticlesSub->SetName(fParticlesSubName); + fJetTask->GetEvent()->AddObject(fParticlesSub); + } else { + AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fParticlesSubName.Data())); + return; + } + } + + if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event + fRhoParam = dynamic_cast(fJetTask->GetEvent()->FindListObject(fRhoName)); + if (!fRhoParam) { + AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data())); + return; + } + } + + if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event + fRhomParam = dynamic_cast(fJetTask->GetEvent()->FindListObject(fRhomName)); + if (!fRhomParam) { + AliError(Form("%s: Could not retrieve rho_m %s!", GetName(), fRhomName.Data())); + return; + } + } + + fInit = kTRUE; +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityConstSubtractor::Prepare(AliFJWrapper& fjw) +{ + // Prepare the utility. + + if (!fInit) return; + + if (fRhoParam) fRho = fRhoParam->GetVal(); + if (fRhomParam) fRhom = fRhomParam->GetVal(); + + if (fJetsSub) fJetsSub->Delete(); + + fjw.SetUseExternalBkg(fUseExternalBkg, fRho, fRhom); + fjw.DoConstituentSubtraction(); +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityConstSubtractor::ProcessJet(AliEmcalJet* /*jet*/, Int_t /*ij*/, AliFJWrapper& /*fjw*/) +{ + // Proceess each jet. +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityConstSubtractor::Terminate(AliFJWrapper& fjw) +{ + // Run termination of the utility (after each event). + + if (!fInit) return; + + if (!fJetsSub) { + AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s", fJetsSubName.Data())); + return; + } + +#ifdef FASTJET_VERSION + std::vector 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()); + + // Fill constituent info + std::vector constituents_unsub(fjw.GetJetConstituents(ijet)); + std::vector constituents_sub = jets_sub[ijet].constituents(); + jet_sub->SetNumberOfTracks(constituents_sub.size()); + fJetTask->FillJetConstituents(jet_sub, constituents_sub, fJetTask->GetTracks(), fJetTask->GetClusters(), constituents_unsub, 1, fParticlesSub); + + jetCount++; + } + +#endif + +} + diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h new file mode 100644 index 00000000000..2d7dc774552 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h @@ -0,0 +1,53 @@ +#ifndef ALIEMCALJETUTILITYCONSTSUBTRACTOR_H +#define ALIEMCALJETUTILITYCONSTSUBTRACTOR_H + +#include + +#include "AliEmcalJetUtility.h" +#include "AliFJWrapper.h" + +class AliEmcalJetTask; +class AliEmcalJet; +class AliFJWrapper; +class AliRhoParameter; + +class AliEmcalJetUtilityConstSubtractor : public AliEmcalJetUtility +{ + public: + + AliEmcalJetUtilityConstSubtractor(); + AliEmcalJetUtilityConstSubtractor(const char* name); + AliEmcalJetUtilityConstSubtractor(const AliEmcalJetUtilityConstSubtractor &jet); + AliEmcalJetUtilityConstSubtractor& operator=(const AliEmcalJetUtilityConstSubtractor &jet); + ~AliEmcalJetUtilityConstSubtractor() {;} + + void SetRhoName(const char *n) { fRhoName = n ; } + void SetRhomName(const char *n) { fRhomName = n ; } + void SetUseExternalBkg(Bool_t b) { fUseExternalBkg = b ; } + + void SetJetsSubName(const char *n) { fJetsSubName = n ; } + void SetParticlesSubName(const char *n) { fParticlesSubName = n ; } + + void Init(); + void Prepare(AliFJWrapper& fjw); + void ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw); + void Terminate(AliFJWrapper& fjw); + + protected: + + TString fJetsSubName; // name of subtracted jet collection + TString fParticlesSubName; // name of subtracted particle collection + Bool_t fUseExternalBkg; // use external background for generic subtractor + TString fRhoName; // name of rho + TString fRhomName; // name of rhom + Double_t fRho; // pT background density + Double_t fRhom; // mT background density + + TClonesArray *fJetsSub; //!subtracted jet collection + TClonesArray *fParticlesSub; //!subtracted particle collection + AliRhoParameter *fRhoParam; //!event rho + AliRhoParameter *fRhomParam; //!event rhom + + ClassDef(AliEmcalJetUtilityConstSubtractor, 1) // Emcal jet utility that implements the constituent subtractor form the fastjet contrib +}; +#endif diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx new file mode 100644 index 00000000000..68f9adf64e6 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx @@ -0,0 +1,256 @@ +#include "AliEmcalJetUtilityGenSubtractor.h" +#include "AliEmcalJet.h" +#include "AliRhoParameter.h" +#include "AliEmcalJetTask.h" + +ClassImp(AliEmcalJetUtilityGenSubtractor) + +//______________________________________________________________________________ +AliEmcalJetUtilityGenSubtractor::AliEmcalJetUtilityGenSubtractor() : +AliEmcalJetUtility(), + fDoGenericSubtractionJetMass(kFALSE), + fDoGenericSubtractionGR(kFALSE), + fDoGenericSubtractionExtraJetShapes(kFALSE), + fUseExternalBkg(kFALSE), + fRhoName(""), + fRhomName(""), + fRho(0), + fRhom(0), + fRMax(0.4), + fDRStep(0.04), + fPtMinGR(40.), + fRhoParam(0), + fRhomParam(0) +{ + // Dummy constructor. + +} + +//______________________________________________________________________________ +AliEmcalJetUtilityGenSubtractor::AliEmcalJetUtilityGenSubtractor(const char* name) : + AliEmcalJetUtility(name), + fDoGenericSubtractionJetMass(kFALSE), + fDoGenericSubtractionGR(kFALSE), + fDoGenericSubtractionExtraJetShapes(kFALSE), + fUseExternalBkg(kFALSE), + fRhoName(""), + fRhomName(""), + fRho(0), + fRhom(0), + fRMax(0.4), + fDRStep(0.04), + fPtMinGR(40.), + fRhoParam(0), + fRhomParam(0) +{ + // Default constructor. +} + +//______________________________________________________________________________ +AliEmcalJetUtilityGenSubtractor::AliEmcalJetUtilityGenSubtractor(const AliEmcalJetUtilityGenSubtractor &other) : + AliEmcalJetUtility(other), + fDoGenericSubtractionJetMass(other.fDoGenericSubtractionJetMass), + fDoGenericSubtractionGR(other.fDoGenericSubtractionGR), + fDoGenericSubtractionExtraJetShapes(other.fDoGenericSubtractionExtraJetShapes), + fUseExternalBkg(other.fUseExternalBkg), + fRhoName(other.fRhoName), + fRhomName(other.fRhomName), + fRho(other.fRho), + fRhom(other.fRhom), + fRMax(other.fRMax), + fDRStep(other.fDRStep), + fPtMinGR(other.fPtMinGR), + fRhoParam(other.fRhoParam), + fRhomParam(other.fRhomParam) +{ + // Copy constructor. +} + +//______________________________________________________________________________ +AliEmcalJetUtilityGenSubtractor& AliEmcalJetUtilityGenSubtractor::operator=(const AliEmcalJetUtilityGenSubtractor &other) +{ + // Assignment. + + if (&other == this) return *this; + AliEmcalJetUtility::operator=(other); + fDoGenericSubtractionJetMass = other.fDoGenericSubtractionJetMass; + fDoGenericSubtractionGR = other.fDoGenericSubtractionGR; + fDoGenericSubtractionExtraJetShapes = other.fDoGenericSubtractionExtraJetShapes; + fUseExternalBkg = other.fUseExternalBkg; + fRhoName = other.fRhoName; + fRhomName = other.fRhomName; + fRho = other.fRho; + fRhom = other.fRhom; + fRMax = other.fRMax; + fDRStep = other.fDRStep; + fPtMinGR = other.fPtMinGR; + fRhoParam = other.fRhoParam; + fRhomParam = other.fRhomParam; + return *this; +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityGenSubtractor::Init() +{ + // Initialize the utility. + + if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event + if(!fJetTask) return; + fRhoParam = dynamic_cast(fJetTask->GetEvent()->FindListObject(fRhoName)); + if (!fRhoParam) { + AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data())); + return; + } + } + + if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event + fRhomParam = dynamic_cast(fJetTask->GetEvent()->FindListObject(fRhomName)); + if (!fRhomParam) { + AliError(Form("%s: Could not retrieve rho_m %s!", GetName(), fRhomName.Data())); + return; + } + } + + fInit = kTRUE; +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityGenSubtractor::Prepare(AliFJWrapper& fjw) +{ + // Prepare the utility. + + if (!fInit) return; + + if (fRhoParam) fRho = fRhoParam->GetVal(); + if (fRhomParam) fRhom = fRhomParam->GetVal(); + + //run generic subtractor + if (fDoGenericSubtractionJetMass) { + fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom); + fjw.DoGenericSubtractionJetMass(); + } + + if (fDoGenericSubtractionExtraJetShapes) { + fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom); + fjw.DoGenericSubtractionJetAngularity(); + fjw.DoGenericSubtractionJetpTD(); + fjw.DoGenericSubtractionJetCircularity(); + fjw.DoGenericSubtractionJetSigma2(); + fjw.DoGenericSubtractionJetConstituent(); + fjw.DoGenericSubtractionJetLeSub(); + } +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityGenSubtractor::ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw) +{ + // Proceess each jet. + + if (!fInit) return; + +#ifdef FASTJET_VERSION + + if (fDoGenericSubtractionJetMass) { + std::vector 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 = jet->Pt()-fjw.GetJetArea(ij)*fRho; + if (fDoGenericSubtractionGR && ptcorr>fPtMinGR) { + fjw.SetUseExternalBkg(fUseExternalBkg, fRho, fRhom); + fRMax = fJetTask->GetRadius()+0.2; + fjw.SetRMaxAndStep(fRMax, fDRStep); + fjw.DoGenericSubtractionGR(ij); + std::vector num = fjw.GetGRNumerator(); + std::vector den = fjw.GetGRDenominator(); + std::vector nums = fjw.GetGRNumeratorSub(); + std::vector 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); + } + } + + if (fDoGenericSubtractionExtraJetShapes) { + std::vector jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity(); + Int_t na = (Int_t)jetAngularityInfo.size(); + if(na > ij && na > 0) { + jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative()); + jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative()); + jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted()); + jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted()); + } + + std::vector jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD(); + Int_t np = (Int_t)jetpTDInfo.size(); + if(np > ij && np > 0) { + jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative()); + jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative()); + jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted()); + jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted()); + } + + std::vector jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity(); + Int_t nc = (Int_t)jetCircularityInfo.size(); + if(nc > ij && nc > 0) { + jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative()); + jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative()); + jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted()); + jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted()); + } + + std::vector jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2(); + Int_t ns = (Int_t)jetSigma2Info.size(); + if (ns > ij && ns > 0) { + jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative()); + jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative()); + jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted()); + jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted()); + } + + + std::vector jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent(); + Int_t nco = (Int_t)jetConstituentInfo.size(); + if(nco > ij && nco > 0) { + jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative()); + jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative()); + jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted()); + jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted()); + } + + std::vector jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub(); + Int_t nlsub = (Int_t)jetLeSubInfo.size(); + if(nlsub > ij && nlsub > 0) { + jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative()); + jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative()); + jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted()); + jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted()); + } + } + +#endif +} + +//______________________________________________________________________________ +void AliEmcalJetUtilityGenSubtractor::Terminate(AliFJWrapper& /*fjw*/) +{ + // Run termination of the utility (after each event). +} diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h new file mode 100644 index 00000000000..f66084b0679 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h @@ -0,0 +1,57 @@ +#ifndef ALIEMCALJETUTILITYGENSUBTRACTOR_H +#define ALIEMCALJETUTILITYGENSUBTRACTOR_H + +#include + +#include "AliEmcalJetUtility.h" +#include "AliFJWrapper.h" + +class AliEmcalJetTask; +class AliEmcalJet; +class AliFJWrapper; +class AliRhoParameter; + +class AliEmcalJetUtilityGenSubtractor : public AliEmcalJetUtility +{ + public: + + AliEmcalJetUtilityGenSubtractor(); + AliEmcalJetUtilityGenSubtractor(const char* name); + AliEmcalJetUtilityGenSubtractor(const AliEmcalJetUtilityGenSubtractor &jet); + AliEmcalJetUtilityGenSubtractor& operator=(const AliEmcalJetUtilityGenSubtractor &jet); + ~AliEmcalJetUtilityGenSubtractor() {;} + + void Init(); + void Prepare(AliFJWrapper& fjw); + void ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw); + void Terminate(AliFJWrapper& fjw); + + void SetRhoName(const char *n) { fRhoName = n ; } + void SetRhomName(const char *n) { fRhomName = n ; } + + void SetGenericSubtractionJetMass(Bool_t b) { fDoGenericSubtractionJetMass = b; } + void SetGenericSubtractionGR(Bool_t b, Double_t rmax = 2., + Double_t dr = 0.04, Double_t ptmin = 0.) { fDoGenericSubtractionGR = b; fRMax = rmax; fDRStep = dr; fPtMinGR = ptmin;} + void SetGenericSubtractionExtraJetShapes(Bool_t b) { fDoGenericSubtractionExtraJetShapes = b; } + void SetUseExternalBkg(Bool_t b) { fUseExternalBkg = b; } + + protected: + + Bool_t fDoGenericSubtractionJetMass; // calculate generic subtraction + Bool_t fDoGenericSubtractionGR; // calculate generic subtraction for angular structure function GR + Bool_t fDoGenericSubtractionExtraJetShapes; // calculate generic subtraction for other jet shapes like radialmoment,pTD etc + Bool_t fUseExternalBkg; // use external background for generic subtractor + TString fRhoName; // name of rho + TString fRhomName; // name of rhom + Double_t fRho; // pT background density + Double_t fRhom; // mT background density + Double_t fRMax; // R max for GR calculation + Double_t fDRStep; // step width for GR calculation + Double_t fPtMinGR; // min pT for GR calculation + + AliRhoParameter *fRhoParam; //!event rho + AliRhoParameter *fRhomParam; //!event rhom + + ClassDef(AliEmcalJetUtilityGenSubtractor, 1) // Emcal jet utility that implements generic subtractors form the fastjet contrib +}; +#endif diff --git a/PWGJE/EMCALJetTasks/AliFJWrapper.h b/PWGJE/EMCALJetTasks/AliFJWrapper.h index 1d4fd71a7f1..c919da4bdbc 100644 --- a/PWGJE/EMCALJetTasks/AliFJWrapper.h +++ b/PWGJE/EMCALJetTasks/AliFJWrapper.h @@ -1,8 +1,6 @@ #ifndef AliFJWrapper_H #define AliFJWrapper_H -// $Id$ - #if !defined(__CINT__) #include @@ -62,6 +60,8 @@ class AliFJWrapper virtual Int_t DoGenericSubtractionJetLeSub(); virtual Int_t DoConstituentSubtraction(); + void SetName(const char* name) { fName = name; } + void SetTitle(const char* title) { fTitle = title; } void SetStrategy(const fastjet::Strategy &strat) { fStrategy = strat; } void SetAlgorithm(const fastjet::JetAlgorithm &algor) { fAlgor = algor; } void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme = scheme; } diff --git a/PWGJE/PWGJEEMCALJetTasksLinkDef.h b/PWGJE/PWGJEEMCALJetTasksLinkDef.h index 2c4dc3a0c7f..baf5143c3e5 100644 --- a/PWGJE/PWGJEEMCALJetTasksLinkDef.h +++ b/PWGJE/PWGJEEMCALJetTasksLinkDef.h @@ -111,6 +111,9 @@ #pragma link C++ class AliAnalysisTaskEmcalTriggerTreeWriter+; #ifdef HAVE_FASTJET +#pragma link C++ class AliEmcalJetUtility+; +#pragma link C++ class AliEmcalJetUtilityGenSubtractor+; +#pragma link C++ class AliEmcalJetUtilityConstSubtractor+; #pragma link C++ class AliEmcalJetTask+; #pragma link C++ class AliEmcalJetFinder+; #pragma link C++ class AliJetEmbeddingFromAODTask+; -- 2.43.0