Fork me on GitHub

Changeset 341014c in git for modules/FastJetFinder.cc


Ignore:
Timestamp:
Feb 12, 2019, 9:29:17 PM (6 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, llp, master
Children:
6455202
Parents:
45e58be
Message:

apply .clang-format to all .h, .cc and .cpp files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/FastJetFinder.cc

    r45e58be r341014c  
    1717 */
    1818
    19 
    2019/** \class FastJetFinder
    2120 *
    … …  
    3231#include "classes/DelphesFormula.h"
    3332
     33#include "ExRootAnalysis/ExRootClassifier.h"
     34#include "ExRootAnalysis/ExRootFilter.h"
    3435#include "ExRootAnalysis/ExRootResult.h"
    35 #include "ExRootAnalysis/ExRootFilter.h"
    36 #include "ExRootAnalysis/ExRootClassifier.h"
    37 
     36
     37#include "TDatabasePDG.h"
     38#include "TFormula.h"
     39#include "TLorentzVector.h"
    3840#include "TMath.h"
     41#include "TObjArray.h"
     42#include "TRandom3.h"
    3943#include "TString.h"
    40 #include "TFormula.h"
    41 #include "TRandom3.h"
    42 #include "TObjArray.h"
    43 #include "TDatabasePDG.h"
    44 #include "TLorentzVector.h"
    4544
    4645#include <algorithm>
    47 #include <stdexcept>
    4846#include <iostream>
    4947#include <sstream>
     48#include <stdexcept>
    5049#include <vector>
    5150
     51#include "fastjet/ClusterSequence.hh"
     52#include "fastjet/ClusterSequenceArea.hh"
     53#include "fastjet/JetDefinition.hh"
    5254#include "fastjet/PseudoJet.hh"
    53 #include "fastjet/JetDefinition.hh"
    54 #include "fastjet/ClusterSequence.hh"
    5555#include "fastjet/Selector.hh"
    56 #include "fastjet/ClusterSequenceArea.hh"
    5756#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
    5857
     58#include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
     59#include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
    5960#include "fastjet/plugins/SISCone/fastjet/SISConePlugin.hh"
    60 #include "fastjet/plugins/CDFCones/fastjet/CDFMidPointPlugin.hh"
    61 #include "fastjet/plugins/CDFCones/fastjet/CDFJetCluPlugin.hh"
    62 
    63 #include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
     61
     62#include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
    6463#include "fastjet/contribs/Nsubjettiness/Njettiness.hh"
    6564#include "fastjet/contribs/Nsubjettiness/NjettinessPlugin.hh"
    66 #include "fastjet/contribs/Nsubjettiness/ExtraRecombiners.hh"
     65#include "fastjet/contribs/Nsubjettiness/Nsubjettiness.hh"
    6766
    6867#include "fastjet/contribs/ValenciaPlugin/ValenciaPlugin.hh"
    6968
     69#include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
    7070#include "fastjet/tools/Filter.hh"
    7171#include "fastjet/tools/Pruner.hh"
    72 #include "fastjet/contribs/RecursiveTools/SoftDrop.hh"
    7372
    7473using namespace std;
    7574using namespace fastjet;
    7675using namespace fastjet::contrib;
    77 
    7876
    7977//------------------------------------------------------------------------------
    … …  
    8381  fDefinition(0), fAreaDefinition(0), fItInputArray(0)
    8482{
    85 
    8683}
    8784
    … …  
    9087FastJetFinder::~FastJetFinder()
    9188{
    92 
    9389}
    9490
    … …  
    126122  fAxisMode = GetInt("AxisMode", 1);
    127123  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)
     124  fN = GetInt("N", 2); // used only if Njettiness is used as jet clustering algo (case 8)
    129125
    130126  //-- Exclusive clustering for e+e- collisions --
    131127
    132   fNJets = GetInt("NJets",2);
     128  fNJets = GetInt("NJets", 2);
    133129  fExclusiveClustering = GetBool("ExclusiveClustering", false);
    134130
    … …  
    142138  switch(fAxisMode)
    143139  {
    144     default:
    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();
     140  default:
     141  case 1:
     142    fAxesDef = new WTA_KT_Axes();
     143    break;
     144  case 2:
     145    fAxesDef = new OnePass_WTA_KT_Axes();
     146    break;
     147  case 3:
     148    fAxesDef = new KT_Axes();
     149    break;
     150  case 4:
     151    fAxesDef = new OnePass_KT_Axes();
    156152  }
    157153
    … …  
    161157  fRTrim = GetDouble("RTrim", 0.2);
    162158  fPtFracTrim = GetDouble("PtFracTrim", 0.05);
    163 
    164159
    165160  //-- Pruning parameters --
    … …  
    172167  //-- SoftDrop parameters --
    173168
    174   fComputeSoftDrop     = GetBool("ComputeSoftDrop", false);
    175   fBetaSoftDrop        = GetDouble("BetaSoftDrop", 0.0);
     169  fComputeSoftDrop = GetBool("ComputeSoftDrop", false);
     170  fBetaSoftDrop = GetDouble("BetaSoftDrop", 0.0);
    176171  fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
    177   fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
     172  fR0SoftDrop = GetDouble("R0SoftDrop=", 0.8);
    178173
    179174  // ---  Jet Area Parameters ---
    … …  
    195190  switch(fAreaAlgorithm)
    196191  {
    197     default:
    198     case 0:
    199       fAreaDefinition = 0;
    200       break;
    201     case 1:
    202       fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    203       break;
    204     case 2:
    205       fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    206       break;
    207     case 3:
    208       fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    209       break;
    210     case 4:
    211       fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
    212       break;
    213     case 5:
    214       fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    215       break;
     192  default:
     193  case 0:
     194    fAreaDefinition = 0;
     195    break;
     196  case 1:
     197    fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
     198    break;
     199  case 2:
     200    fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
     201    break;
     202  case 3:
     203    fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
     204    break;
     205  case 4:
     206    fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
     207    break;
     208  case 5:
     209    fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
     210    break;
    216211  }
    217212
    218213  switch(fJetAlgorithm)
    219214  {
    220     case 1:
    221       plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
    222       fDefinition = new JetDefinition(plugin);
    223       break;
    224     case 2:
    225       plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
    226       fDefinition = new JetDefinition(plugin);
    227       break;
    228     case 3:
    229       plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
    230       fDefinition = new JetDefinition(plugin);
    231       break;
    232     case 4:
    233       fDefinition = new JetDefinition(kt_algorithm, fParameterR);
    234       break;
    235     case 5:
    236       fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
    237       break;
    238     default:
    239     case 6:
    240       fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
    241       break;
    242     case 7:
    243       recomb = new WinnerTakeAllRecombiner();
    244       fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
    245       break;
    246     case 8:
    247       fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
    248       fDefinition = new JetDefinition(fNjettinessPlugin);
    249       break;
    250     case 9:
    251       fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
    252       fDefinition = new JetDefinition(fValenciaPlugin);
    253       break;
    254 
    255   }
    256 
    257 
     215  case 1:
     216    plugin = new CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
     217    fDefinition = new JetDefinition(plugin);
     218    break;
     219  case 2:
     220    plugin = new CDFMidPointPlugin(fSeedThreshold, fConeRadius, fConeAreaFraction, fMaxPairSize, fMaxIterations, fOverlapThreshold);
     221    fDefinition = new JetDefinition(plugin);
     222    break;
     223  case 3:
     224    plugin = new SISConePlugin(fConeRadius, fOverlapThreshold, fMaxIterations, fJetPTMin);
     225    fDefinition = new JetDefinition(plugin);
     226    break;
     227  case 4:
     228    fDefinition = new JetDefinition(kt_algorithm, fParameterR);
     229    break;
     230  case 5:
     231    fDefinition = new JetDefinition(cambridge_algorithm, fParameterR);
     232    break;
     233  default:
     234  case 6:
     235    fDefinition = new JetDefinition(antikt_algorithm, fParameterR);
     236    break;
     237  case 7:
     238    recomb = new WinnerTakeAllRecombiner();
     239    fDefinition = new JetDefinition(antikt_algorithm, fParameterR, recomb, Best);
     240    break;
     241  case 8:
     242    fNjettinessPlugin = new NjettinessPlugin(fN, Njettiness::wta_kt_axes, Njettiness::unnormalized_cutoff_measure, fBeta, fRcutOff);
     243    fDefinition = new JetDefinition(fNjettinessPlugin);
     244    break;
     245  case 9:
     246    fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
     247    fDefinition = new JetDefinition(fValenciaPlugin);
     248    break;
     249  }
    258250
    259251  fPlugin = plugin;
    … …  
    270262
    271263    fEstimators.clear();
    272     for(i = 0; i < size/2; ++i)
    273     {
    274       etaMin = param[i*2].GetDouble();
    275       etaMax = param[i*2 + 1].GetDouble();
     264    for(i = 0; i < size / 2; ++i)
     265    {
     266      etaMin = param[i * 2].GetDouble();
     267      etaMax = param[i * 2 + 1].GetDouble();
    276268      estimatorStruct.estimator = new JetMedianBackgroundEstimator(SelectorRapRange(etaMin, etaMax), *fDefinition, *fAreaDefinition);
    277269      estimatorStruct.etaMin = etaMin;
    … …  
    297289void FastJetFinder::Finish()
    298290{
    299   vector< TEstimatorStruct >::iterator itEstimators;
     291  vector<TEstimatorStruct>::iterator itEstimators;
    300292
    301293  for(itEstimators = fEstimators.begin(); itEstimators != fEstimators.end(); ++itEstimators)
    … …  
    307299  if(fDefinition) delete fDefinition;
    308300  if(fAreaDefinition) delete fAreaDefinition;
    309   if(fPlugin) delete static_cast<JetDefinition::Plugin*>(fPlugin);
    310   if(fRecomb) delete static_cast<JetDefinition::Recombiner*>(fRecomb);
    311   if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin*>(fNjettinessPlugin);
     301  if(fPlugin) delete static_cast<JetDefinition::Plugin *>(fPlugin);
     302  if(fRecomb) delete static_cast<JetDefinition::Recombiner *>(fRecomb);
     303  if(fNjettinessPlugin) delete static_cast<JetDefinition::Plugin *>(fNjettinessPlugin);
    312304  if(fAxesDef) delete fAxesDef;
    313305  if(fMeasureDef) delete fMeasureDef;
    314   if(fValenciaPlugin) delete static_cast<JetDefinition::Plugin*>(fValenciaPlugin);
    315 
     306  if(fValenciaPlugin) delete static_cast<JetDefinition::Plugin *>(fValenciaPlugin);
    316307}
    317308
    … …  
    330321  PseudoJet jet, area;
    331322  ClusterSequence *sequence;
    332   vector< PseudoJet > inputList, outputList, subjets;
    333   vector< PseudoJet >::iterator itInputList, itOutputList;
    334   vector< TEstimatorStruct >::iterator itEstimators;
     323  vector<PseudoJet> inputList, outputList, subjets;
     324  vector<PseudoJet>::iterator itInputList, itOutputList;
     325  vector<TEstimatorStruct>::iterator itEstimators;
    335326  Double_t excl_ymerge23 = 0.0;
    336327  Double_t excl_ymerge34 = 0.0;
    … …  
    345336  fItInputArray->Reset();
    346337  number = 0;
    347   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
     338  while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
    348339  {
    349340    momentum = candidate->Momentum;
    … …  
    383374
    384375  if(fExclusiveClustering)
    385     {
    386    try{
    387      outputList = sorted_by_pt(sequence->exclusive_jets( fNJets ));
    388    }
    389    catch(fastjet::Error)
    390      {
    391        outputList.clear();
    392      }
    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 );
    398     }
     376  {
     377    try
     378    {
     379      outputList = sorted_by_pt(sequence->exclusive_jets(fNJets));
     380    }
     381    catch(fastjet::Error)
     382    {
     383      outputList.clear();
     384    }
     385
     386    excl_ymerge23 = sequence->exclusive_ymerge(2);
     387    excl_ymerge34 = sequence->exclusive_ymerge(3);
     388    excl_ymerge45 = sequence->exclusive_ymerge(4);
     389    excl_ymerge56 = sequence->exclusive_ymerge(5);
     390  }
    399391  else
    400     {
    401       outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
    402     }
     392  {
     393    outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));
     394  }
    403395
    404396  // loop over all jets and export them
    … …  
    432424    {
    433425      if(itInputList->user_index() < 0) continue;
    434       constituent = static_cast<Candidate*>(fInputArray->At(itInputList->user_index()));
     426      constituent = static_cast<Candidate *>(fInputArray->At(itInputList->user_index()));
    435427
    436428      deta = TMath::Abs(momentum.Eta() - constituent->Momentum.Eta());
    … …  
    439431      if(dphi > dphiMax) dphiMax = dphi;
    440432
    441       if(constituent->Charge == 0) nneutrals++;
    442       else ncharged++;
    443 
    444       time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T());
     433      if(constituent->Charge == 0)
     434        nneutrals++;
     435      else
     436        ncharged++;
     437
     438      time += TMath::Sqrt(constituent->Momentum.E()) * (constituent->Position.T());
    445439      timeWeight += TMath::Sqrt(constituent->Momentum.E());
    446440
    … …  
    452446
    453447    candidate->Momentum = momentum;
    454     candidate->Position.SetT(time/timeWeight);
     448    candidate->Position.SetT(time / timeWeight);
    455449    candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E());
    456450
    … …  
    461455    candidate->NCharged = ncharged;
    462456
    463 
    464457    //for exclusive clustering, access y_n,n+1 as exclusive_ymerge (fNJets);
    465458    candidate->ExclYmerge23 = excl_ymerge23;
    … …  
    475468    {
    476469
    477       fastjet::Filter    trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
     470      fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, fRTrim), fastjet::SelectorPtFractionMin(fPtFracTrim));
    478471      fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
    479472
    … …  
    489482      candidate->NSubJetsTrimmed = subjets.size();
    490483
    491       for (size_t i = 0; i < subjets.size() and i < 4; i++)
     484      for(size_t i = 0; i < subjets.size() and i < 4; i++)
    492485      {
    493             if(subjets.at(i).pt() < 0) continue ;
    494             candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
     486        if(subjets.at(i).pt() < 0) continue;
     487        candidate->TrimmedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    495488      }
    496489    }
    497 
    498490
    499491    //------------------------------------
    … …  
    501493    //------------------------------------
    502494
    503 
    504495    if(fComputePruning)
    505496    {
    506497
    507       fastjet::Pruner    pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm,fRPrun),fZcutPrun,fRcutPrun);
     498      fastjet::Pruner pruner(fastjet::JetDefinition(fastjet::cambridge_algorithm, fRPrun), fZcutPrun, fRcutPrun);
    508499      fastjet::PseudoJet pruned_jet = pruner(*itOutputList);
    509500
    … …  
    517508      candidate->NSubJetsPruned = subjets.size();
    518509
    519       for (size_t i = 0; i < subjets.size() and i < 4; i++)
     510      for(size_t i = 0; i < subjets.size() and i < 4; i++)
    520511      {
    521             if(subjets.at(i).pt() < 0) continue ;
    522             candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
     512        if(subjets.at(i).pt() < 0) continue;
     513        candidate->PrunedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    523514      }
    524 
    525515    }
    526516
    … …  
    532522    {
    533523
    534       contrib::SoftDrop softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
     524      contrib::SoftDrop softDrop(fBetaSoftDrop, fSymmetryCutSoftDrop, fR0SoftDrop);
    535525      fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
    536526
    … …  
    540530
    541531      subjets.clear();
    542       subjets    = softdrop_jet.pieces();
    543       subjets    = sorted_by_pt(subjets);
     532      subjets = softdrop_jet.pieces();
     533      subjets = sorted_by_pt(subjets);
    544534      candidate->NSubJetsSoftDropped = softdrop_jet.pieces().size();
    545535
    546536      candidate->SoftDroppedJet = candidate->SoftDroppedP4[0];
    547537
    548       for (size_t i = 0; i < subjets.size() and i < 4; i++)
     538      for(size_t i = 0; i < subjets.size() and i < 4; i++)
    549539      {
    550             if(subjets.at(i).pt() < 0) continue ;
    551             candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    552             if(i==0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i+1];
    553             if(i==1) candidate->SoftDroppedSubJet2 = candidate->SoftDroppedP4[i+1];
     540        if(subjets.at(i).pt() < 0) continue;
     541        candidate->SoftDroppedP4[i + 1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
     542        if(i == 0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i + 1];
     543        if(i == 1) candidate->SoftDroppedSubJet2 = candidate->SoftDroppedP4[i + 1];
    554544      }
    555545    }
    … …  
    571561      candidate->Tau[3] = nSub4(*itOutputList);
    572562      candidate->Tau[4] = nSub5(*itOutputList);
    573 
    574563    }
    575564
Note: See TracChangeset for help on using the changeset viewer.