[01f457a] | 1 | /*
|
---|
| 2 | * Delphes: a framework for fast simulation of a generic collider experiment
|
---|
| 3 | * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
|
---|
[1fa50c2] | 4 | *
|
---|
[01f457a] | 5 | * This program is free software: you can redistribute it and/or modify
|
---|
| 6 | * it under the terms of the GNU General Public License as published by
|
---|
| 7 | * the Free Software Foundation, either version 3 of the License, or
|
---|
| 8 | * (at your option) any later version.
|
---|
[1fa50c2] | 9 | *
|
---|
[01f457a] | 10 | * This program is distributed in the hope that it will be useful,
|
---|
| 11 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 12 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 13 | * GNU General Public License for more details.
|
---|
[1fa50c2] | 14 | *
|
---|
[01f457a] | 15 | * You should have received a copy of the GNU General Public License
|
---|
| 16 | * along with this program. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 17 | */
|
---|
| 18 |
|
---|
[d7d2da3] | 19 |
|
---|
| 20 | /** \class FastJetFinder
|
---|
| 21 | *
|
---|
| 22 | * Finds jets using FastJet library.
|
---|
| 23 | *
|
---|
| 24 | * \author P. Demin - UCL, Louvain-la-Neuve
|
---|
| 25 | *
|
---|
| 26 | */
|
---|
| 27 |
|
---|
| 28 | #include "modules/FastJetFinder.h"
|
---|
| 29 |
|
---|
| 30 | #include "classes/DelphesClasses.h"
|
---|
| 31 | #include "classes/DelphesFactory.h"
|
---|
| 32 | #include "classes/DelphesFormula.h"
|
---|
| 33 |
|
---|
| 34 | #include "ExRootAnalysis/ExRootResult.h"
|
---|
| 35 | #include "ExRootAnalysis/ExRootFilter.h"
|
---|
| 36 | #include "ExRootAnalysis/ExRootClassifier.h"
|
---|
| 37 |
|
---|
| 38 | #include "TMath.h"
|
---|
| 39 | #include "TString.h"
|
---|
| 40 | #include "TFormula.h"
|
---|
| 41 | #include "TRandom3.h"
|
---|
| 42 | #include "TObjArray.h"
|
---|
| 43 | #include "TDatabasePDG.h"
|
---|
| 44 | #include "TLorentzVector.h"
|
---|
| 45 |
|
---|
[8336b6e] | 46 | #include <algorithm>
|
---|
[d7d2da3] | 47 | #include <stdexcept>
|
---|
| 48 | #include <iostream>
|
---|
| 49 | #include <sstream>
|
---|
| 50 | #include <vector>
|
---|
| 51 |
|
---|
| 52 | #include "fastjet/PseudoJet.hh"
|
---|
| 53 | #include "fastjet/JetDefinition.hh"
|
---|
| 54 | #include "fastjet/ClusterSequence.hh"
|
---|
| 55 | #include "fastjet/Selector.hh"
|
---|
| 56 | #include "fastjet/ClusterSequenceArea.hh"
|
---|
| 57 | #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
|
---|
| 58 |
|
---|
| 59 | #include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
|
---|
| 60 | #include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
|
---|
| 61 | #include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
|
---|
| 62 |
|
---|
[e4c3fef] | 63 | #include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
|
---|
[9687203] | 64 | #include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
|
---|
| 65 | #include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
|
---|
[973b92a] | 66 | #include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
|
---|
[9687203] | 67 |
|
---|
[f319c1d] | 68 | #include "fastjet/contribs/ValenciaPlugin/ValenciaPlugin.hh"
|
---|
| 69 |
|
---|
[de6d698] | 70 | #include "fastjet/tools/Filter.hh"
|
---|
| 71 | #include "fastjet/tools/Pruner.hh"
|
---|
| 72 | #include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
|
---|
| 73 |
|
---|
[d7d2da3] | 74 | using namespace std;
|
---|
| 75 | using namespace fastjet;
|
---|
[9687203] | 76 | using namespace fastjet::contrib;
|
---|
| 77 |
|
---|
[d7d2da3] | 78 |
|
---|
| 79 | //------------------------------------------------------------------------------
|
---|
| 80 |
|
---|
| 81 | FastJetFinder::FastJetFinder() :
|
---|
[e9c0d73] | 82 | fPlugin(0), fRecomb(0), fAxesDef(0), fMeasureDef(0), fNjettinessPlugin(0), fValenciaPlugin(0),
|
---|
| 83 | fDefinition(0), fAreaDefinition(0), fItInputArray(0)
|
---|
[d7d2da3] | 84 | {
|
---|
| 85 |
|
---|
| 86 | }
|
---|
| 87 |
|
---|
| 88 | //------------------------------------------------------------------------------
|
---|
| 89 |
|
---|
| 90 | FastJetFinder::~FastJetFinder()
|
---|
| 91 | {
|
---|
| 92 |
|
---|
| 93 | }
|
---|
| 94 |
|
---|
| 95 | //------------------------------------------------------------------------------
|
---|
| 96 |
|
---|
| 97 | void FastJetFinder::Init()
|
---|
| 98 | {
|
---|
[c6321ad] | 99 | JetDefinition::Plugin *plugin = 0;
|
---|
| 100 | JetDefinition::Recombiner *recomb = 0;
|
---|
[7278220] | 101 | ExRootConfParam param;
|
---|
[8336b6e] | 102 | Long_t i, size;
|
---|
[7278220] | 103 | Double_t etaMin, etaMax;
|
---|
| 104 | TEstimatorStruct estimatorStruct;
|
---|
[8336b6e] | 105 |
|
---|
[d7d2da3] | 106 | // define algorithm
|
---|
| 107 |
|
---|
| 108 | fJetAlgorithm = GetInt("JetAlgorithm", 6);
|
---|
| 109 | fParameterR = GetDouble("ParameterR", 0.5);
|
---|
| 110 |
|
---|
| 111 | fConeRadius = GetDouble("ConeRadius", 0.5);
|
---|
| 112 | fSeedThreshold = GetDouble("SeedThreshold", 1.0);
|
---|
| 113 | fConeAreaFraction = GetDouble("ConeAreaFraction", 1.0);
|
---|
| 114 | fMaxIterations = GetInt("MaxIterations", 100);
|
---|
| 115 | fMaxPairSize = GetInt("MaxPairSize", 2);
|
---|
| 116 | fIratch = GetInt("Iratch", 1);
|
---|
[3ccc8586] | 117 | fAdjacencyCut = GetInt("AdjacencyCut", 2);
|
---|
[d7d2da3] | 118 | fOverlapThreshold = GetDouble("OverlapThreshold", 0.75);
|
---|
| 119 |
|
---|
| 120 | fJetPTMin = GetDouble("JetPTMin", 10.0);
|
---|
| 121 |
|
---|
[9687203] | 122 | //-- N(sub)jettiness parameters --
|
---|
[e4c3fef] | 123 |
|
---|
[9687203] | 124 | fComputeNsubjettiness = GetBool("ComputeNsubjettiness", false);
|
---|
| 125 | fBeta = GetDouble("Beta", 1.0);
|
---|
| 126 | fAxisMode = GetInt("AxisMode", 1);
|
---|
[e4c3fef] | 127 | fRcutOff = GetDouble("RcutOff", 0.8); // used only if Njettiness is used as jet clustering algo (case 8)
|
---|
| 128 | fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
|
---|
[f319c1d] | 129 |
|
---|
| 130 | //-- Exclusive clustering for e+e- collisions --
|
---|
[3b3d0ad] | 131 |
|
---|
[f319c1d] | 132 | fNJets = GetInt("NJets",2);
|
---|
| 133 | fExclusiveClustering = GetBool("ExclusiveClustering", false);
|
---|
| 134 |
|
---|
| 135 | //-- Valencia Linear Collider algorithm
|
---|
[3b3d0ad] | 136 |
|
---|
[f319c1d] | 137 | fGamma = GetDouble("Gamma", 1.0);
|
---|
| 138 | //fBeta parameter see above
|
---|
[3b3d0ad] | 139 |
|
---|
[919f6a6] | 140 | fMeasureDef = new NormalizedMeasure(fBeta, fParameterR);
|
---|
[3b3d0ad] | 141 |
|
---|
[8abab33] | 142 | switch(fAxisMode)
|
---|
| 143 | {
|
---|
| 144 | default:
|
---|
[3b3d0ad] | 145 | case 1:
|
---|
| 146 | fAxesDef = new WTA_KT_Axes();
|
---|
| 147 | break;
|
---|
| 148 | case 2:
|
---|
| 149 | fAxesDef = new OnePass_WTA_KT_Axes();
|
---|
| 150 | break;
|
---|
| 151 | case 3:
|
---|
| 152 | fAxesDef = new KT_Axes();
|
---|
| 153 | break;
|
---|
| 154 | case 4:
|
---|
| 155 | fAxesDef = new OnePass_KT_Axes();
|
---|
| 156 | }
|
---|
[e4c3fef] | 157 |
|
---|
[de6d698] | 158 | //-- Trimming parameters --
|
---|
[3b3d0ad] | 159 |
|
---|
[de6d698] | 160 | fComputeTrimming = GetBool("ComputeTrimming", false);
|
---|
| 161 | fRTrim = GetDouble("RTrim", 0.2);
|
---|
| 162 | fPtFracTrim = GetDouble("PtFracTrim", 0.05);
|
---|
[3b3d0ad] | 163 |
|
---|
[de6d698] | 164 |
|
---|
| 165 | //-- Pruning parameters --
|
---|
[3b3d0ad] | 166 |
|
---|
[de6d698] | 167 | fComputePruning = GetBool("ComputePruning", false);
|
---|
| 168 | fZcutPrun = GetDouble("ZcutPrun", 0.1);
|
---|
| 169 | fRcutPrun = GetDouble("RcutPrun", 0.5);
|
---|
| 170 | fRPrun = GetDouble("RPrun", 0.8);
|
---|
[3b3d0ad] | 171 |
|
---|
[de6d698] | 172 | //-- SoftDrop parameters --
|
---|
[3b3d0ad] | 173 |
|
---|
[de6d698] | 174 | fComputeSoftDrop = GetBool("ComputeSoftDrop", false);
|
---|
| 175 | fBetaSoftDrop = GetDouble("BetaSoftDrop", 0.0);
|
---|
| 176 | fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
|
---|
| 177 | fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
|
---|
| 178 |
|
---|
[d7d2da3] | 179 | // --- Jet Area Parameters ---
|
---|
[3b3d0ad] | 180 |
|
---|
[d7d2da3] | 181 | fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
|
---|
| 182 | fComputeRho = GetBool("ComputeRho", false);
|
---|
[e4c3fef] | 183 |
|
---|
[d7d2da3] | 184 | // - ghost based areas -
|
---|
| 185 | fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
|
---|
| 186 | fRepeat = GetInt("Repeat", 1);
|
---|
| 187 | fGhostArea = GetDouble("GhostArea", 0.01);
|
---|
| 188 | fGridScatter = GetDouble("GridScatter", 1.0);
|
---|
| 189 | fPtScatter = GetDouble("PtScatter", 0.1);
|
---|
| 190 | fMeanGhostPt = GetDouble("MeanGhostPt", 1.0E-100);
|
---|
[e4c3fef] | 191 |
|
---|
[d7d2da3] | 192 | // - voronoi based areas -
|
---|
| 193 | fEffectiveRfact = GetDouble("EffectiveRfact", 1.0);
|
---|
| 194 |
|
---|
| 195 | switch(fAreaAlgorithm)
|
---|
| 196 | {
|
---|
[3b3d0ad] | 197 | default:
|
---|
| 198 | case 0:
|
---|
| 199 | fAreaDefinition = 0;
|
---|
| 200 | break;
|
---|
[d7d2da3] | 201 | case 1:
|
---|
[c6321ad] | 202 | fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
[d7d2da3] | 203 | break;
|
---|
| 204 | case 2:
|
---|
[c6321ad] | 205 | fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
[d7d2da3] | 206 | break;
|
---|
| 207 | case 3:
|
---|
[c6321ad] | 208 | fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
[d7d2da3] | 209 | break;
|
---|
| 210 | case 4:
|
---|
[c6321ad] | 211 | fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
|
---|
[d7d2da3] | 212 | break;
|
---|
| 213 | case 5:
|
---|
[c6321ad] | 214 | fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
|
---|
[d7d2da3] | 215 | break;
|
---|
| 216 | }
|
---|
| 217 |
|
---|
| 218 | switch(fJetAlgorithm)
|
---|
| 219 | {
|
---|
[8336b6e] | 220 | case 1:
|
---|
[c6321ad] | 221 | plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
|
---|
| 222 | fDefinition = new JetDefinition(plugin);
|
---|
[d7d2da3] | 223 | break;
|
---|
| 224 | case 2:
|
---|
[c6321ad] | 225 | plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
|
---|
| 226 | fDefinition = new JetDefinition(plugin);
|
---|
[d7d2da3] | 227 | break;
|
---|
| 228 | case 3:
|
---|
[c6321ad] | 229 | plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
|
---|
| 230 | fDefinition = new JetDefinition(plugin);
|
---|
[d7d2da3] | 231 | break;
|
---|
| 232 | case 4:
|
---|
[c6321ad] | 233 | fDefinition = new JetDefinition(kt_algorithm, fParameterR);
|
---|
[d7d2da3] | 234 | break;
|
---|
| 235 | case 5:
|
---|
[c6321ad] | 236 | fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
|
---|
[d7d2da3] | 237 | break;
|
---|
| 238 | default:
|
---|
| 239 | case 6:
|
---|
[c6321ad] | 240 | fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
|
---|
[d7d2da3] | 241 | break;
|
---|
[9687203] | 242 | case 7:
|
---|
[c6321ad] | 243 | recomb = new WinnerTakeAllRecombiner();
|
---|
| 244 | fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
|
---|
[9687203] | 245 | break;
|
---|
| 246 | case 8:
|
---|
[7278220] | 247 | fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
|
---|
| 248 | fDefinition = new JetDefinition(fNjettinessPlugin);
|
---|
[9687203] | 249 | break;
|
---|
[3b3d0ad] | 250 | case 9:
|
---|
[f319c1d] | 251 | fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
|
---|
| 252 | fDefinition = new JetDefinition(fValenciaPlugin);
|
---|
| 253 | break;
|
---|
| 254 |
|
---|
[d7d2da3] | 255 | }
|
---|
[8336b6e] | 256 |
|
---|
[3b3d0ad] | 257 |
|
---|
[8abab33] | 258 |
|
---|
[d7d2da3] | 259 | fPlugin = plugin;
|
---|
[9687203] | 260 | fRecomb = recomb;
|
---|
[e4c3fef] | 261 |
|
---|
[d7d2da3] | 262 | ClusterSequence::print_banner();
|
---|
[8336b6e] | 263 |
|
---|
[7278220] | 264 | if(fComputeRho && fAreaDefinition)
|
---|
| 265 | {
|
---|
| 266 | // read eta ranges
|
---|
| 267 |
|
---|
| 268 | param = GetParam("RhoEtaRange");
|
---|
| 269 | size = param.GetSize();
|
---|
| 270 |
|
---|
| 271 | fEstimators.clear();
|
---|
| 272 | for(i = 0; i < size/2; ++i)
|
---|
| 273 | {
|
---|
| 274 | etaMin = param[i*2].GetDouble();
|
---|
| 275 | etaMax = param[i*2 + 1].GetDouble();
|
---|
[58dcd2e] | 276 | estimatorStruct.estimator = new JetMedianBackgroundEstimator(SelectorRapRange(etaMin, etaMax), *fDefinition, *fAreaDefinition);
|
---|
[7278220] | 277 | estimatorStruct.etaMin = etaMin;
|
---|
| 278 | estimatorStruct.etaMax = etaMax;
|
---|
| 279 | fEstimators.push_back(estimatorStruct);
|
---|
| 280 | }
|
---|
| 281 | }
|
---|
| 282 |
|
---|
[d7d2da3] | 283 | // import input array
|
---|
| 284 |
|
---|
| 285 | fInputArray = ImportArray(GetString("InputArray", "Calorimeter/towers"));
|
---|
| 286 | fItInputArray = fInputArray->MakeIterator();
|
---|
| 287 |
|
---|
| 288 | // create output arrays
|
---|
| 289 |
|
---|
| 290 | fOutputArray = ExportArray(GetString("OutputArray", "jets"));
|
---|
| 291 | fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
|
---|
[70bb4cb] | 292 | fConstituentsOutputArray = ExportArray(GetString("ConstituentsOutputArray", "constituents"));
|
---|
[d7d2da3] | 293 | }
|
---|
| 294 |
|
---|
| 295 | //------------------------------------------------------------------------------
|
---|
| 296 |
|
---|
| 297 | void FastJetFinder::Finish()
|
---|
| 298 | {
|
---|
[7278220] | 299 | vector< TEstimatorStruct >::iterator itEstimators;
|
---|
| 300 |
|
---|
| 301 | for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
|
---|
| 302 | {
|
---|
| 303 | if(itEstimators->estimator) delete itEstimators->estimator;
|
---|
| 304 | }
|
---|
| 305 |
|
---|
[d7d2da3] | 306 | if(fItInputArray) delete fItInputArray;
|
---|
| 307 | if(fDefinition) delete fDefinition;
|
---|
| 308 | if(fAreaDefinition) delete fAreaDefinition;
|
---|
| 309 | if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
|
---|
[9687203] | 310 | if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
|
---|
| 311 | if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
|
---|
[8abab33] | 312 | if(fAxesDef) delete fAxesDef;
|
---|
| 313 | if(fMeasureDef) delete fMeasureDef;
|
---|
[f319c1d] | 314 | if(fValenciaPlugin) delete static_cast<JetDefinition::Plugin*>(fValenciaPlugin);
|
---|
| 315 |
|
---|
[d7d2da3] | 316 | }
|
---|
| 317 |
|
---|
| 318 | //------------------------------------------------------------------------------
|
---|
| 319 |
|
---|
| 320 | void FastJetFinder::Process()
|
---|
| 321 | {
|
---|
| 322 | Candidate *candidate, *constituent;
|
---|
| 323 | TLorentzVector momentum;
|
---|
[df04eb1] | 324 |
|
---|
[d7d2da3] | 325 | Double_t deta, dphi, detaMax, dphiMax;
|
---|
[df04eb1] | 326 | Double_t time, timeWeight;
|
---|
[ce4feac] | 327 | Int_t number, ncharged, nneutrals;
|
---|
[3b3d0ad] | 328 | Int_t charge;
|
---|
[df04eb1] | 329 | Double_t rho = 0.0;
|
---|
[d7d2da3] | 330 | PseudoJet jet, area;
|
---|
| 331 | ClusterSequence *sequence;
|
---|
[de6d698] | 332 | vector< PseudoJet > inputList, outputList, subjets;
|
---|
[7278220] | 333 | vector< PseudoJet >::iterator itInputList, itOutputList;
|
---|
| 334 | vector< TEstimatorStruct >::iterator itEstimators;
|
---|
[e9c0d73] | 335 | Double_t excl_ymerge23 = 0.0;
|
---|
| 336 | Double_t excl_ymerge34 = 0.0;
|
---|
| 337 | Double_t excl_ymerge45 = 0.0;
|
---|
| 338 | Double_t excl_ymerge56 = 0.0;
|
---|
[3b3d0ad] | 339 |
|
---|
[d7d2da3] | 340 | DelphesFactory *factory = GetFactory();
|
---|
| 341 |
|
---|
| 342 | inputList.clear();
|
---|
| 343 |
|
---|
| 344 | // loop over input objects
|
---|
| 345 | fItInputArray->Reset();
|
---|
| 346 | number = 0;
|
---|
| 347 | while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
|
---|
[8336b6e] | 348 | {
|
---|
[d7d2da3] | 349 | momentum = candidate->Momentum;
|
---|
| 350 | jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
|
---|
| 351 | jet.set_user_index(number);
|
---|
| 352 | inputList.push_back(jet);
|
---|
| 353 | ++number;
|
---|
| 354 | }
|
---|
[8336b6e] | 355 |
|
---|
| 356 | // construct jets
|
---|
[d7d2da3] | 357 | if(fAreaDefinition)
|
---|
| 358 | {
|
---|
| 359 | sequence = new ClusterSequenceArea(inputList, *fDefinition, *fAreaDefinition);
|
---|
| 360 | }
|
---|
| 361 | else
|
---|
| 362 | {
|
---|
| 363 | sequence = new ClusterSequence(inputList, *fDefinition);
|
---|
[8336b6e] | 364 | }
|
---|
[d7d2da3] | 365 |
|
---|
| 366 | // compute rho and store it
|
---|
| 367 | if(fComputeRho && fAreaDefinition)
|
---|
| 368 | {
|
---|
[7278220] | 369 | for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
|
---|
[8336b6e] | 370 | {
|
---|
[7278220] | 371 | itEstimators->estimator->set_particles(inputList);
|
---|
| 372 | rho = itEstimators->estimator->rho();
|
---|
[8336b6e] | 373 |
|
---|
| 374 | candidate = factory->NewCandidate();
|
---|
| 375 | candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
|
---|
[7278220] | 376 | candidate->Edges[0] = itEstimators->etaMin;
|
---|
| 377 | candidate->Edges[1] = itEstimators->etaMax;
|
---|
[8336b6e] | 378 | fRhoOutputArray->Add(candidate);
|
---|
| 379 | }
|
---|
[d7d2da3] | 380 | }
|
---|
[8336b6e] | 381 |
|
---|
[d7d2da3] | 382 | outputList.clear();
|
---|
| 383 |
|
---|
[f319c1d] | 384 | if(fExclusiveClustering)
|
---|
| 385 | {
|
---|
[35807af] | 386 | try{
|
---|
| 387 | outputList = sorted_by_pt(sequence->exclusive_jets( fNJets ));
|
---|
| 388 | }
|
---|
| 389 | catch(fastjet::Error)
|
---|
| 390 | {
|
---|
| 391 | outputList.clear();
|
---|
| 392 | }
|
---|
[e9c0d73] | 393 |
|
---|
| 394 | excl_ymerge23 = sequence->exclusive_ymerge( 2 );
|
---|
| 395 | excl_ymerge34 = sequence->exclusive_ymerge( 3 );
|
---|
| 396 | excl_ymerge45 = sequence->exclusive_ymerge( 4 );
|
---|
| 397 | excl_ymerge56 = sequence->exclusive_ymerge( 5 );
|
---|
[f319c1d] | 398 | }
|
---|
| 399 | else
|
---|
| 400 | {
|
---|
| 401 | outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
|
---|
| 402 | }
|
---|
[9687203] | 403 |
|
---|
[d7d2da3] | 404 | // loop over all jets and export them
|
---|
| 405 | detaMax = 0.0;
|
---|
| 406 | dphiMax = 0.0;
|
---|
[3b3d0ad] | 407 |
|
---|
[d7d2da3] | 408 | for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
|
---|
| 409 | {
|
---|
[d244bc9] | 410 | jet = *itOutputList;
|
---|
| 411 | if(fJetAlgorithm == 7) jet = join(jet.constituents());
|
---|
[df04eb1] | 412 |
|
---|
[d244bc9] | 413 | momentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
|
---|
[df04eb1] | 414 |
|
---|
[d7d2da3] | 415 | area.reset(0.0, 0.0, 0.0, 0.0);
|
---|
| 416 | if(fAreaDefinition) area = itOutputList->area_4vector();
|
---|
| 417 |
|
---|
| 418 | candidate = factory->NewCandidate();
|
---|
| 419 |
|
---|
[df04eb1] | 420 | time = 0.0;
|
---|
| 421 | timeWeight = 0.0;
|
---|
[22dc7fd] | 422 |
|
---|
[5a44a72] | 423 | charge = 0;
|
---|
| 424 |
|
---|
[ce4feac] | 425 | ncharged = 0;
|
---|
| 426 | nneutrals = 0;
|
---|
| 427 |
|
---|
[d7d2da3] | 428 | inputList.clear();
|
---|
| 429 | inputList = sequence->constituents(*itOutputList);
|
---|
[e4c3fef] | 430 |
|
---|
[d7d2da3] | 431 | for(itInputList = inputList.begin(); itInputList != inputList.end(); ++itInputList)
|
---|
| 432 | {
|
---|
[29f8a06] | 433 | if(itInputList->user_index() < 0) continue;
|
---|
[d7d2da3] | 434 | constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
|
---|
| 435 |
|
---|
| 436 | deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
|
---|
| 437 | dphi = TMath::Abs(momentum.DeltaPhi(constituent->Momentum));
|
---|
| 438 | if(deta > detaMax) detaMax = deta;
|
---|
| 439 | if(dphi > dphiMax) dphiMax = dphi;
|
---|
[e4c3fef] | 440 |
|
---|
[ce4feac] | 441 | if(constituent->Charge == 0) nneutrals++;
|
---|
| 442 | else ncharged++;
|
---|
| 443 |
|
---|
[22dc7fd] | 444 | time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
|
---|
[df04eb1] | 445 | timeWeight += TMath::Sqrt(constituent->Momentum.E());
|
---|
[e4c3fef] | 446 |
|
---|
[5a44a72] | 447 | charge += constituent->Charge;
|
---|
| 448 |
|
---|
[70bb4cb] | 449 | fConstituentsOutputArray->Add(constituent);
|
---|
[d7d2da3] | 450 | candidate->AddCandidate(constituent);
|
---|
| 451 | }
|
---|
[e4c3fef] | 452 |
|
---|
[d7d2da3] | 453 | candidate->Momentum = momentum;
|
---|
[df04eb1] | 454 | candidate->Position.SetT(time/timeWeight);
|
---|
[d7d2da3] | 455 | candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
|
---|
| 456 |
|
---|
| 457 | candidate->DeltaEta = detaMax;
|
---|
| 458 | candidate->DeltaPhi = dphiMax;
|
---|
[3b3d0ad] | 459 | candidate->Charge = charge;
|
---|
[ce4feac] | 460 | candidate->NNeutrals = nneutrals;
|
---|
| 461 | candidate->NCharged = ncharged;
|
---|
[e9c0d73] | 462 |
|
---|
| 463 |
|
---|
| 464 | //for exclusive clustering, access y_n,n+1 as exclusive_ymerge (fNJets);
|
---|
| 465 | candidate->ExclYmerge23 = excl_ymerge23;
|
---|
| 466 | candidate->ExclYmerge34 = excl_ymerge34;
|
---|
| 467 | candidate->ExclYmerge45 = excl_ymerge45;
|
---|
| 468 | candidate->ExclYmerge56 = excl_ymerge56;
|
---|
[3b3d0ad] | 469 |
|
---|
[de6d698] | 470 | //------------------------------------
|
---|
| 471 | // Trimming
|
---|
| 472 | //------------------------------------
|
---|
[d7d2da3] | 473 |
|
---|
[8abab33] | 474 | if(fComputeTrimming)
|
---|
[de6d698] | 475 | {
|
---|
| 476 |
|
---|
| 477 | fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
|
---|
| 478 | fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
|
---|
[3b3d0ad] | 479 |
|
---|
[de6d698] | 480 | trimmed_jet = join(trimmed_jet.constituents());
|
---|
[3b3d0ad] | 481 |
|
---|
[de6d698] | 482 | candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
|
---|
[3b3d0ad] | 483 |
|
---|
| 484 | // four hardest subjets
|
---|
[de6d698] | 485 | subjets.clear();
|
---|
| 486 | subjets = trimmed_jet.pieces();
|
---|
| 487 | subjets = sorted_by_pt(subjets);
|
---|
[3b3d0ad] | 488 |
|
---|
[de6d698] | 489 | candidate->NSubJetsTrimmed = subjets.size();
|
---|
| 490 |
|
---|
[8abab33] | 491 | for (size_t i = 0; i < subjets.size() and i < 4; i++)
|
---|
| 492 | {
|
---|
[3b3d0ad] | 493 | if(subjets.at(i).pt() < 0) continue ;
|
---|
[8abab33] | 494 | candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
|
---|
[de6d698] | 495 | }
|
---|
| 496 | }
|
---|
[3b3d0ad] | 497 |
|
---|
| 498 |
|
---|
[de6d698] | 499 | //------------------------------------
|
---|
| 500 | // Pruning
|
---|
| 501 | //------------------------------------
|
---|
[3b3d0ad] | 502 |
|
---|
| 503 |
|
---|
[de6d698] | 504 | if(fComputePruning)
|
---|
| 505 | {
|
---|
| 506 |
|
---|
| 507 | fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun);
|
---|
| 508 | fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
|
---|
| 509 |
|
---|
| 510 | candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
|
---|
[3b3d0ad] | 511 |
|
---|
| 512 | // four hardest subjet
|
---|
[de6d698] | 513 | subjets.clear();
|
---|
| 514 | subjets = pruned_jet.pieces();
|
---|
| 515 | subjets = sorted_by_pt(subjets);
|
---|
[3b3d0ad] | 516 |
|
---|
[de6d698] | 517 | candidate->NSubJetsPruned = subjets.size();
|
---|
| 518 |
|
---|
[8abab33] | 519 | for (size_t i = 0; i < subjets.size() and i < 4; i++)
|
---|
| 520 | {
|
---|
[3b3d0ad] | 521 | if(subjets.at(i).pt() < 0) continue ;
|
---|
[8abab33] | 522 | candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
|
---|
[de6d698] | 523 | }
|
---|
| 524 |
|
---|
[3b3d0ad] | 525 | }
|
---|
| 526 |
|
---|
[de6d698] | 527 | //------------------------------------
|
---|
| 528 | // SoftDrop
|
---|
| 529 | //------------------------------------
|
---|
[3b3d0ad] | 530 |
|
---|
[de6d698] | 531 | if(fComputeSoftDrop)
|
---|
| 532 | {
|
---|
[3b3d0ad] | 533 |
|
---|
[8abab33] | 534 | contrib::SoftDrop softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
|
---|
[de6d698] | 535 | fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
|
---|
[3b3d0ad] | 536 |
|
---|
[de6d698] | 537 | candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
|
---|
[3b3d0ad] | 538 |
|
---|
| 539 | // four hardest subjet
|
---|
| 540 |
|
---|
[de6d698] | 541 | subjets.clear();
|
---|
| 542 | subjets = softdrop_jet.pieces();
|
---|
| 543 | subjets = sorted_by_pt(subjets);
|
---|
| 544 | candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
|
---|
| 545 |
|
---|
[ba75867] | 546 | candidate->SoftDroppedJet = candidate->SoftDroppedP4[0];
|
---|
| 547 |
|
---|
[8abab33] | 548 | for (size_t i = 0; i < subjets.size() and i < 4; i++)
|
---|
| 549 | {
|
---|
[3b3d0ad] | 550 | if(subjets.at(i).pt() < 0) continue ;
|
---|
[8abab33] | 551 | candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
|
---|
[ba75867] | 552 | if(i==0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i+1];
|
---|
| 553 | if(i==1) candidate->SoftDroppedSubJet2 = candidate->SoftDroppedP4[i+1];
|
---|
[de6d698] | 554 | }
|
---|
| 555 | }
|
---|
[3b3d0ad] | 556 |
|
---|
[9687203] | 557 | // --- compute N-subjettiness with N = 1,2,3,4,5 ----
|
---|
[e4c3fef] | 558 |
|
---|
[9687203] | 559 | if(fComputeNsubjettiness)
|
---|
| 560 | {
|
---|
[3b3d0ad] | 561 |
|
---|
[8abab33] | 562 | Nsubjettiness nSub1(1, *fAxesDef, *fMeasureDef);
|
---|
| 563 | Nsubjettiness nSub2(2, *fAxesDef, *fMeasureDef);
|
---|
| 564 | Nsubjettiness nSub3(3, *fAxesDef, *fMeasureDef);
|
---|
| 565 | Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
|
---|
| 566 | Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
|
---|
[3b3d0ad] | 567 |
|
---|
[e4c3fef] | 568 | candidate->Tau[0] = nSub1(*itOutputList);
|
---|
| 569 | candidate->Tau[1] = nSub2(*itOutputList);
|
---|
| 570 | candidate->Tau[2] = nSub3(*itOutputList);
|
---|
| 571 | candidate->Tau[3] = nSub4(*itOutputList);
|
---|
| 572 | candidate->Tau[4] = nSub5(*itOutputList);
|
---|
[3b3d0ad] | 573 |
|
---|
[9687203] | 574 | }
|
---|
| 575 |
|
---|
[d7d2da3] | 576 | fOutputArray->Add(candidate);
|
---|
| 577 | }
|
---|
| 578 | delete sequence;
|
---|
| 579 | }
|
---|