X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FEMCALJetTasks%2FAliEmcalJetTask.cxx;h=17f55eaf846106e263f90e469ef76a7f76995d08;hb=9196956a6ca4a867f0aabe4a44be192c90e1b1dc;hp=c0ef3b551d7595a1000a6f1afcdad84c8063a8b6;hpb=7071e730d31f5772cb9422fc3a5860cd473be6cf;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx index c0ef3b551d7..17f55eaf846 100644 --- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx +++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx @@ -1,8 +1,7 @@ -// $Id$ // // Emcal jet finder task. // -// Authors: C.Loizides, S.Aiola +// Authors: C.Loizides, S.Aiola, M. Verweij #include #include "AliEmcalJetTask.h" @@ -25,6 +24,7 @@ #include "AliVCluster.h" #include "AliVEvent.h" #include "AliVParticle.h" +#include "AliRhoParameter.h" using std::cout; using std::endl; using std::cerr; @@ -37,6 +37,7 @@ AliEmcalJetTask::AliEmcalJetTask() : fTracksName("Tracks"), fCaloName("CaloClusters"), fJetsName("Jets"), + fJetsSubName(""), fJetType(kNone), fConstSel(0), fMCConstSel(0), @@ -56,16 +57,31 @@ AliEmcalJetTask::AliEmcalJetTask() : fJetEtaMax(+1), fGhostArea(0.005), fMinMCLabel(0), - fRecombScheme(-1), + fRecombScheme(fastjet::pt_scheme), fTrackEfficiency(1.), + fMCFlag(0), + fGeneratorIndex(-1), fIsInit(0), + fLocked(0), fIsPSelSet(0), fIsMcPart(0), fIsEmcPart(0), + fLegacyMode(kFALSE), + fCodeDebug(kFALSE), + fDoGenericSubtraction(kFALSE), + fDoConstituentSubtraction(kFALSE), + fUseExternalBkg(kFALSE), + fRhoName(""), + fRhomName(""), + fRho(0), + fRhom(0), fJets(0), + fJetsSub(0), fEvent(0), fTracks(0), - fClus(0) + fClus(0), + fRhoParam(0), + fRhomParam(0) { // Default constructor. } @@ -76,6 +92,7 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) : fTracksName("Tracks"), fCaloName("CaloClusters"), fJetsName("Jets"), + fJetsSubName(""), fJetType(kAKT|kFullJet|kRX1Jet), fConstSel(0), fMCConstSel(0), @@ -95,16 +112,31 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) : fJetEtaMax(+1), fGhostArea(0.005), fMinMCLabel(0), - fRecombScheme(-1), + fRecombScheme(fastjet::pt_scheme), fTrackEfficiency(1.), + fMCFlag(0), + fGeneratorIndex(-1), fIsInit(0), + fLocked(0), fIsPSelSet(0), fIsMcPart(0), fIsEmcPart(0), + fLegacyMode(kFALSE), + fCodeDebug(kFALSE), + fDoGenericSubtraction(kFALSE), + fDoConstituentSubtraction(kFALSE), + fUseExternalBkg(kFALSE), + fRhoName(""), + fRhomName(""), + fRho(0), + fRhom(0), fJets(0), + fJetsSub(0), fEvent(0), fTracks(0), - fClus(0) + fClus(0), + fRhoParam(0), + fRhomParam(0) { // Standard constructor. @@ -124,6 +156,11 @@ void AliEmcalJetTask::UserCreateOutputObjects() fJets = new TClonesArray("AliEmcalJet"); fJets->SetName(fJetsName); + + if(!fJetsSubName.IsNull()) { + fJetsSub = new TClonesArray("AliEmcalJet"); + fJetsSub->SetName(fJetsSubName); + } } //________________________________________________________________________ @@ -138,6 +175,7 @@ void AliEmcalJetTask::UserExec(Option_t *) // clear the jet array (normally a null operation) fJets->Delete(); + if(fJetsSub) fJetsSub->Delete(); FindJets(); } @@ -157,6 +195,11 @@ void AliEmcalJetTask::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) { @@ -181,14 +224,13 @@ void AliEmcalJetTask::FindJets() fjw.SetGhostArea(fGhostArea); fjw.SetR(fRadius); fjw.SetAlgorithm(jalgo); - if(fRecombScheme>0) - fjw.SetRecombScheme(static_cast(fRecombScheme)); + fjw.SetRecombScheme(static_cast(fRecombScheme)); fjw.SetMaxRap(fEtaMax); fjw.Clear(); // get primary vertex Double_t vertex[3] = {0, 0, 0}; - fEvent->GetPrimaryVertex()->GetXYZ(vertex); + if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex); AliDebug(2,Form("Jet type = %d", fJetType)); @@ -225,7 +267,15 @@ void AliEmcalJetTask::FindJets() 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(); @@ -245,7 +295,7 @@ 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->P(), iTracks + 100); + fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100); } } @@ -334,10 +384,27 @@ void AliEmcalJetTask::FindJets() fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100); } } - + + // setting legacy mode + if (fLegacyMode) { + fjw.SetLegacyMode(kTRUE); + } + // run jet finder fjw.Run(); + //run generic subtractor + if(fDoGenericSubtraction) { + 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) { @@ -366,6 +433,21 @@ void AliEmcalJetTask::FindJets() 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 + if(fDoGenericSubtraction) { +#ifdef FASTJET_CONTRIB + std::vector jetMassInfo = fjw.GetGenSubtractorInfoJetMass(); + UInt_t n = (UInt_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()); + } +#endif + } // loop over constituents std::vector constituents(fjw.GetJetConstituents(ij)); @@ -523,6 +605,30 @@ void AliEmcalJetTask::FindJets() jetCount++; } //fJets->Sort(); + + //run constituent subtractor if requested + if(fDoConstituentSubtraction) { +#ifdef FASTJET_CONTRIB + 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); + 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 } //________________________________________________________________________ @@ -573,6 +679,9 @@ Bool_t AliEmcalJetTask::DoInit() return 0; } + if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub) + fEvent->AddObject(fJetsSub); + if (fTracksName == "Tracks") am->LoadBranch("Tracks"); if (!fTracks && !fTracksName.IsNull()) { @@ -603,5 +712,20 @@ Bool_t AliEmcalJetTask::DoInit() fIsEmcPart = 1; } + 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 %s!", GetName(), fRhomName.Data())); + return 0; + } + } + return 1; }