Fork me on GitHub

Changeset 8336b6e in git for modules/FastJetFinder.cc


Ignore:
Timestamp:
Nov 4, 2013, 11:59:27 AM (11 years ago)
Author:
pavel <pavel@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
e9065e7
Parents:
498cda4
Message:

add eta ranges for rho

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/FastJetFinder.cc

    r498cda4 r8336b6e  
    3030#include "TLorentzVector.h"
    3131
    32 #include <algorithm> 
     32#include <algorithm>
    3333#include <stdexcept>
    3434#include <iostream>
     
    7070{
    7171  JetDefinition::Plugin *plugin = NULL;
     72
     73  // read eta ranges
     74
     75  ExRootConfParam param = GetParam("RhoEtaRange");
     76  Long_t i, size;
     77
     78  fEtaRangeMap.clear();
     79  size = param.GetSize();
     80  for(i = 0; i < size/2; ++i)
     81  {
     82    fEtaRangeMap[param[i*2].GetDouble()] = param[i*2 + 1].GetDouble();
     83  }
    7284
    7385  // define algorithm
     
    90102  fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
    91103  fComputeRho = GetBool("ComputeRho", false);
    92   fRhoEtaMax = GetDouble("RhoEtaMax", 5.0);
    93104  // - ghost based areas -
    94105  fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
     
    126137  switch(fJetAlgorithm)
    127138  {
    128     case 1: 
     139    case 1:
    129140      plugin = new fastjet::CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
    130141      fDefinition = new fastjet::JetDefinition(plugin);
     
    149160      break;
    150161  }
    151  
     162
    152163  fPlugin = plugin;
    153164
    154165  ClusterSequence::print_banner();
    155  
     166
    156167  // import input array
    157168
     
    187198  vector<PseudoJet> inputList, outputList;
    188199  ClusterSequence *sequence;
     200  map< Double_t, Double_t >::iterator itEtaRangeMap;
    189201
    190202  DelphesFactory *factory = GetFactory();
     
    196208  number = 0;
    197209  while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
    198   {   
     210  {
    199211    momentum = candidate->Momentum;
    200212    jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
     
    203215    ++number;
    204216  }
    205    
    206   // construct jets 
     217
     218  // construct jets
    207219  if(fAreaDefinition)
    208220  {
     
    212224  {
    213225    sequence = new ClusterSequence(inputList, *fDefinition);
    214   } 
     226  }
    215227
    216228  // compute rho and store it
    217229  if(fComputeRho && fAreaDefinition)
    218230  {
    219     Selector select_rapidity = SelectorAbsRapMax(fRhoEtaMax);
    220     JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
    221     estimator.set_particles(inputList);
    222     rho = estimator.rho();
    223 
    224     candidate = factory->NewCandidate();
    225     candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
    226     fRhoOutputArray->Add(candidate);
    227   }
    228  
     231    for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
     232    {
     233      Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
     234      JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
     235      estimator.set_particles(inputList);
     236      rho = estimator.rho();
     237
     238      candidate = factory->NewCandidate();
     239      candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
     240      candidate->Edges[0] = itEtaRangeMap->first;
     241      candidate->Edges[1] = itEtaRangeMap->second;
     242      fRhoOutputArray->Add(candidate);
     243    }
     244  }
     245
    229246  outputList.clear();
    230247  outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
Note: See TracChangeset for help on using the changeset viewer.