Fork me on GitHub

Changes in / [dd64cff:8a11cdc] in git


Ignore:
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • examples/CMakeLists.txt

    rdd64cff r8a11cdc  
    11include_directories(
    22  ${CMAKE_SOURCE_DIR}
    3   ${DelphesExternals_INCLUDE_DIR} 
     3  ${DelphesExternals_INCLUDE_DIR}
    44)
    55
    … …  
    1717# take all other relevant files and put them into examples/
    1818install(FILES ${macros} DESTINATION examples)
     19install(DIRECTORY ExternalFastJet DESTINATION examples)
     20if(PYTHIA8_FOUND)
     21  install(DIRECTORY Pythia8 DESTINATION examples)
     22endif()

  • modules/CMakeLists.txt

    rdd64cff r8a11cdc  
    2828if (NOT ${ROOT_VERSION} VERSION_LESS "6.0.0")
    2929  install(FILES
    30       ${PROJECT_BINARY_DIR}/modules/libModulesDict_rdict.pcm
    31       ${PROJECT_BINARY_DIR}/modules/libFastJetDict_rdict.pcm
     30    ${PROJECT_BINARY_DIR}/modules/libModulesDict_rdict.pcm
     31    ${PROJECT_BINARY_DIR}/modules/libFastJetDict_rdict.pcm
     32    DESTINATION lib)
     33  if(PYTHIA8_FOUND)
     34    install(FILES
     35      ${PROJECT_BINARY_DIR}/modules/libPythia8Dict_rdict.pcm
    3236      DESTINATION lib)
     37  endif()
    3338endif()

  • modules/FastJetFinder.cc

    rdd64cff r8a11cdc  
    120120  fJetPTMin = GetDouble("JetPTMin", 10.0);
    121121
    122  
    123122  //-- N(sub)jettiness parameters --
    124123
    … …  
    130129
    131130  //-- Exclusive clustering for e+e- collisions --
    132  
     131
    133132  fNJets = GetInt("NJets",2);
    134133  fExclusiveClustering = GetBool("ExclusiveClustering", false);
    135134
    136135  //-- Valencia Linear Collider algorithm
     136
    137137  fGamma = GetDouble("Gamma", 1.0);
    138138  //fBeta parameter see above
    139  
     139
    140140  fMeasureDef = new NormalizedMeasure(fBeta, fParameterR);
    141    
     141
    142142  switch(fAxisMode)
    143143  {
    144144    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();
    156    }
     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  }
    157157
    158158  //-- Trimming parameters --
    159  
     159
    160160  fComputeTrimming = GetBool("ComputeTrimming", false);
    161161  fRTrim = GetDouble("RTrim", 0.2);
    162162  fPtFracTrim = GetDouble("PtFracTrim", 0.05);
    163  
     163
    164164
    165165  //-- Pruning parameters --
    166  
     166
    167167  fComputePruning = GetBool("ComputePruning", false);
    168168  fZcutPrun = GetDouble("ZcutPrun", 0.1);
    169169  fRcutPrun = GetDouble("RcutPrun", 0.5);
    170170  fRPrun = GetDouble("RPrun", 0.8);
    171  
     171
    172172  //-- SoftDrop parameters --
    173  
     173
    174174  fComputeSoftDrop     = GetBool("ComputeSoftDrop", false);
    175175  fBetaSoftDrop        = GetDouble("BetaSoftDrop", 0.0);
    176176  fSymmetryCutSoftDrop = GetDouble("SymmetryCutSoftDrop", 0.1);
    177177  fR0SoftDrop= GetDouble("R0SoftDrop=", 0.8);
    178  
    179178
    180179  // ---  Jet Area Parameters ---
     180
    181181  fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
    182182  fComputeRho = GetBool("ComputeRho", false);
    … …  
    195195  switch(fAreaAlgorithm)
    196196  {
    197     case 1:
    198       fAreaDefinition = new AreaDefinition(active_area_explicit_ghosts, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    199       break;
    200     case 2:
    201       fAreaDefinition = new AreaDefinition(one_ghost_passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    202       break;
    203     case 3:
    204       fAreaDefinition = new AreaDefinition(passive_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    205       break;
    206     case 4:
    207       fAreaDefinition = new AreaDefinition(VoronoiAreaSpec(fEffectiveRfact));
    208       break;
    209     case 5:
    210       fAreaDefinition = new AreaDefinition(active_area, GhostedAreaSpec(fGhostEtaMax, fRepeat, fGhostArea, fGridScatter, fPtScatter, fMeanGhostPt));
    211       break;
    212197    default:
    213198    case 0:
    214199      fAreaDefinition = 0;
    215200      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;
    216216  }
    217217
    … …  
    248248      fDefinition = new JetDefinition(fNjettinessPlugin);
    249249      break;
    250   case 9:
     250    case 9:
    251251      fValenciaPlugin = new ValenciaPlugin(fParameterR, fBeta, fGamma);
    252252      fDefinition = new JetDefinition(fValenciaPlugin);
    … …  
    255255  }
    256256
    257    
     257
    258258
    259259  fPlugin = plugin;
    … …  
    290290  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
    291291  fRhoOutputArray = ExportArray(GetString("RhoOutputArray", "rho"));
     292  fConstituentsOutputArray = ExportArray(GetString("ConstituentsOutputArray", "constituents"));
    292293}
    293294
    … …  
    325326  Double_t time, timeWeight;
    326327  Int_t number, ncharged, nneutrals;
    327   Int_t charge; 
     328  Int_t charge;
    328329  Double_t rho = 0.0;
    329330  PseudoJet jet, area;
    … …  
    336337  Double_t excl_ymerge45 = 0.0;
    337338  Double_t excl_ymerge56 = 0.0;
    338  
     339
    339340  DelphesFactory *factory = GetFactory();
    340341
    … …  
    381382  outputList.clear();
    382383
    383  
    384384  if(fExclusiveClustering)
    385385    {
    … …  
    405405  detaMax = 0.0;
    406406  dphiMax = 0.0;
    407  
     407
    408408  for(itOutputList = outputList.begin(); itOutputList != outputList.end(); ++itOutputList)
    409409  {
    … …  
    416416    if(fAreaDefinition) area = itOutputList->area_4vector();
    417417
    418 
    419    
    420418    candidate = factory->NewCandidate();
    421419
    … …  
    449447      charge += constituent->Charge;
    450448
     449      fConstituentsOutputArray->Add(constituent);
    451450      candidate->AddCandidate(constituent);
    452451    }
    … …  
    458457    candidate->DeltaEta = detaMax;
    459458    candidate->DeltaPhi = dphiMax;
    460     candidate->Charge = charge; 
     459    candidate->Charge = charge;
    461460    candidate->NNeutrals = nneutrals;
    462461    candidate->NCharged = ncharged;
    … …  
    468467    candidate->ExclYmerge45 = excl_ymerge45;
    469468    candidate->ExclYmerge56 = excl_ymerge56;
    470    
     469
    471470    //------------------------------------
    472471    // Trimming
    … …  
    478477      fastjet::Filter    trimmer(fastjet::JetDefinition(fastjet::kt_algorithm,fRTrim),fastjet::SelectorPtFractionMin(fPtFracTrim));
    479478      fastjet::PseudoJet trimmed_jet = trimmer(*itOutputList);
    480      
     479
    481480      trimmed_jet = join(trimmed_jet.constituents());
    482      
     481
    483482      candidate->TrimmedP4[0].SetPtEtaPhiM(trimmed_jet.pt(), trimmed_jet.eta(), trimmed_jet.phi(), trimmed_jet.m());
    484        
    485       // four hardest subjets 
     483
     484      // four hardest subjets
    486485      subjets.clear();
    487486      subjets = trimmed_jet.pieces();
    488487      subjets = sorted_by_pt(subjets);
    489      
     488
    490489      candidate->NSubJetsTrimmed = subjets.size();
    491490
    492491      for (size_t i = 0; i < subjets.size() and i < 4; i++)
    493492      {
    494             if(subjets.at(i).pt() < 0) continue ; 
     493            if(subjets.at(i).pt() < 0) continue ;
    495494            candidate->TrimmedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    496495      }
    497496    }
    498    
    499    
     497
     498
    500499    //------------------------------------
    501500    // Pruning
    502501    //------------------------------------
    503    
    504    
     502
     503
    505504    if(fComputePruning)
    506505    {
    … …  
    510509
    511510      candidate->PrunedP4[0].SetPtEtaPhiM(pruned_jet.pt(), pruned_jet.eta(), pruned_jet.phi(), pruned_jet.m());
    512          
    513       // four hardest subjet 
     511
     512      // four hardest subjet
    514513      subjets.clear();
    515514      subjets = pruned_jet.pieces();
    516515      subjets = sorted_by_pt(subjets);
    517      
     516
    518517      candidate->NSubJetsPruned = subjets.size();
    519518
    520519      for (size_t i = 0; i < subjets.size() and i < 4; i++)
    521520      {
    522             if(subjets.at(i).pt() < 0) continue ; 
     521            if(subjets.at(i).pt() < 0) continue ;
    523522            candidate->PrunedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    524523      }
    525524
    526     } 
    527      
     525    }
     526
    528527    //------------------------------------
    529528    // SoftDrop
    530529    //------------------------------------
    531    
     530
    532531    if(fComputeSoftDrop)
    533532    {
    534    
     533
    535534      contrib::SoftDrop softDrop(fBetaSoftDrop,fSymmetryCutSoftDrop,fR0SoftDrop);
    536535      fastjet::PseudoJet softdrop_jet = softDrop(*itOutputList);
    537      
     536
    538537      candidate->SoftDroppedP4[0].SetPtEtaPhiM(softdrop_jet.pt(), softdrop_jet.eta(), softdrop_jet.phi(), softdrop_jet.m());
    539        
    540       // four hardest subjet 
    541      
     538
     539      // four hardest subjet
     540
    542541      subjets.clear();
    543542      subjets    = softdrop_jet.pieces();
    … …  
    549548      for (size_t i = 0; i < subjets.size()  and i < 4; i++)
    550549      {
    551             if(subjets.at(i).pt() < 0) continue ; 
     550            if(subjets.at(i).pt() < 0) continue ;
    552551            candidate->SoftDroppedP4[i+1].SetPtEtaPhiM(subjets.at(i).pt(), subjets.at(i).eta(), subjets.at(i).phi(), subjets.at(i).m());
    553552            if(i==0) candidate->SoftDroppedSubJet1 = candidate->SoftDroppedP4[i+1];
    … …  
    555554      }
    556555    }
    557  
     556
    558557    // --- compute N-subjettiness with N = 1,2,3,4,5 ----
    559558
    560559    if(fComputeNsubjettiness)
    561560    {
    562      
     561
    563562      Nsubjettiness nSub1(1, *fAxesDef, *fMeasureDef);
    564563      Nsubjettiness nSub2(2, *fAxesDef, *fMeasureDef);
    … …  
    566565      Nsubjettiness nSub4(4, *fAxesDef, *fMeasureDef);
    567566      Nsubjettiness nSub5(5, *fAxesDef, *fMeasureDef);
    568      
     567
    569568      candidate->Tau[0] = nSub1(*itOutputList);
    570569      candidate->Tau[1] = nSub2(*itOutputList);
    … …  
    572571      candidate->Tau[3] = nSub4(*itOutputList);
    573572      candidate->Tau[4] = nSub5(*itOutputList);
    574          
     573
    575574    }
    576575

  • modules/FastJetFinder.h

    rdd64cff r8a11cdc  
    4242    class NjettinessPlugin;
    4343    class ValenciaPlugin;
    44     class AxesDefinition; 
    45     class MeasureDefinition;   
     44    class AxesDefinition;
     45    class MeasureDefinition;
    4646  }
    4747}
    … …  
    6363  void *fRecomb; //!
    6464
    65   fastjet::contrib::AxesDefinition *fAxesDef; 
    66   fastjet::contrib::MeasureDefinition *fMeasureDef;   
     65  fastjet::contrib::AxesDefinition *fAxesDef;
     66  fastjet::contrib::MeasureDefinition *fMeasureDef;
    6767
    6868  fastjet::contrib::NjettinessPlugin *fNjettinessPlugin; //!
    … …  
    9090  //-- Valencia Linear Collider algorithm
    9191  Double_t fGamma;
    92  
     92
    9393  //-- N (sub)jettiness parameters --
    9494
    … …  
    100100
    101101  //-- Trimming parameters --
    102  
     102
    103103  Bool_t fComputeTrimming;
    104104  Double_t fRTrim;
    105105  Double_t fPtFracTrim;
    106  
     106
    107107  //-- Pruning parameters --
    108108
    … …  
    152152  TObjArray *fOutputArray; //!
    153153  TObjArray *fRhoOutputArray; //!
     154  TObjArray *fConstituentsOutputArray; //!
    154155
    155156  ClassDef(FastJetFinder, 1)

  • modules/PdgCodeFilter.cc

    rdd64cff r8a11cdc  
    7676  fInvert = GetBool("Invert", false);
    7777
     78  // no pileup
     79  fRequireNotPileup = GetBool("RequireNotPileup", false);
     80
    7881  fRequireStatus = GetBool("RequireStatus", false);
    7982  fStatus = GetInt("Status", 1);
    … …  
    127130    if(fRequireStatus && (candidate->Status != fStatus)) continue;
    128131    if(fRequireCharge && (candidate->Charge != fCharge)) continue;
     132    if(fRequireNotPileup && (candidate->IsPU   >0     )) continue;
    129133
    130134    pass = kTRUE;

  • modules/PdgCodeFilter.h

    rdd64cff r8a11cdc  
    5555  Bool_t fRequireCharge; //!
    5656  Int_t fCharge; //!
    57 
     57  Bool_t fRequireNotPileup; //!
    5858
    5959  std::vector<Int_t> fPdgCodes;

  • modules/StatusPidFilter.cc

    rdd64cff r8a11cdc  
    153153{
    154154  // PT threshold
    155 
    156155  fPTMin = GetDouble("PTMin", 0.5);
     156
     157  // keep or remove pileup particles
     158  fRequireNotPileup = GetBool("RequireNotPileup", false);
    157159
    158160  // import input array
    … …  
    227229    if(!pass || (candidate->Momentum.Pt() < fPTMin && !(is_b_hadron || is_b_quark || is_tau_daughter || is_W_daughter)) ) continue;
    228230
     231    // not pileup particles
     232    if(fRequireNotPileup && (candidate->IsPU >0)) continue;
     233
    229234    fOutputArray->Add(candidate);
    230235  }

  • modules/StatusPidFilter.h

    rdd64cff r8a11cdc  
    5151  Double_t fPTMin; //!
    5252
     53  Bool_t fRequireNotPileup; //!
     54
    5355  TIterator *fItInputArray; //!
    5456

  • modules/UniqueObjectFinder.cc

    rdd64cff r8a11cdc  
    6767void UniqueObjectFinder::Init()
    6868{
     69  // use GetUniqueID algorithm to find unique objects (faster than the default Overlaps method)
     70  fUseUniqueID = GetBool("UseUniqueID", false);
     71
    6972  // import arrays with output from other modules
    7073
    … …  
    146149    while((previousCandidate = static_cast<Candidate*>(iterator.Next())))
    147150    {
    148       if(candidate->Overlaps(previousCandidate))
     151      if(fUseUniqueID)
    149152      {
    150         return kFALSE;
     153        if(candidate->GetUniqueID() == previousCandidate->GetUniqueID())
     154        {
     155          return kFALSE;
     156        }
     157      }
     158      else
     159      {
     160        if(candidate->Overlaps(previousCandidate))
     161        {
     162          return kFALSE;
     163        }
    151164      }
    152165    }

  • modules/UniqueObjectFinder.h

    rdd64cff r8a11cdc  
    5050private:
    5151
     52  Bool_t fUseUniqueID;
     53
    5254  Bool_t Unique(Candidate *candidate, std::vector< std::pair< TIterator *, TObjArray * > >::iterator itInputMap);
    5355

  • readers/DelphesProMC.cpp

    rdd64cff r8a11cdc  
    125125    candidate->Status = status;
    126126
     127    candidate->IsPU=0;
     128    if (mutableParticles->barcode(i)>0) candidate->IsPU=1; // pileup particle
     129
    127130    candidate->M1 = mutableParticles->mother1(i);
    128131    candidate->M2 = mutableParticles->mother2(i);
    … …  
    231234      cout << "** Reading " << argv[i] << endl;
    232235
     236      // use 64 bit
     237      //inputFile = new ProMCBook(argv[i], "r", true);
     238
     239      //use 32 bit zip (faster but limitted to 64k events)
    233240      inputFile = new ProMCBook(argv[i], "r");
    234241
Note: See TracChangeset for help on using the changeset viewer.