Changeset 8336b6e in git

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

Location:

modules

Files:

• FastJetFinder.cc (modified) (9 diffs)
• FastJetFinder.h (modified) (4 diffs)
• JetPileUpSubtractor.cc (modified) (5 diffs)
• JetPileUpSubtractor.h (modified) (1 diff)

modules/FastJetFinder.cc

@@ -30 +30 @@
 #include "TLorentzVector.h"
 
-#include <algorithm> 
+#include <algorithm>
 #include <stdexcept>
 #include <iostream>
@@ -70 +70 @@
 {
   JetDefinition::Plugin *plugin = NULL;
+
+  // read eta ranges
+
+  ExRootConfParam param = GetParam("RhoEtaRange");
+  Long_t i, size;
+
+  fEtaRangeMap.clear();
+  size = param.GetSize();
+  for(i = 0; i < size/2; ++i)
+  {
+    fEtaRangeMap[param[i*2].GetDouble()] = param[i*2 + 1].GetDouble();
+  }
 
   // define algorithm
@@ -90 +102 @@
   fAreaAlgorithm = GetInt("AreaAlgorithm", 0);
   fComputeRho = GetBool("ComputeRho", false);
-  fRhoEtaMax = GetDouble("RhoEtaMax", 5.0);
   // - ghost based areas -
   fGhostEtaMax = GetDouble("GhostEtaMax", 5.0);
@@ -126 +137 @@
   switch(fJetAlgorithm)
   {
-    case 1: 
+    case 1:
       plugin = new fastjet::CDFJetCluPlugin(fSeedThreshold, fConeRadius, fAdjacencyCut, fMaxIterations, fIratch, fOverlapThreshold);
       fDefinition = new fastjet::JetDefinition(plugin);
@@ -149 +160 @@
       break;
   }
-  
+
   fPlugin = plugin;
 
   ClusterSequence::print_banner();
-  
+
   // import input array
 
@@ -187 +198 @@
   vector<PseudoJet> inputList, outputList;
   ClusterSequence *sequence;
+  map< Double_t, Double_t >::iterator itEtaRangeMap;
 
   DelphesFactory *factory = GetFactory();
@@ -196 +208 @@
   number = 0;
   while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
-  {   
+  {
     momentum = candidate->Momentum;
     jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
@@ -203 +215 @@
     ++number;
   }
-   
-  // construct jets 
+
+  // construct jets
   if(fAreaDefinition)
   {
@@ -212 +224 @@
   {
     sequence = new ClusterSequence(inputList, *fDefinition);
-  } 
+  }
 
   // compute rho and store it
   if(fComputeRho && fAreaDefinition)
   {
-    Selector select_rapidity = SelectorAbsRapMax(fRhoEtaMax);
-    JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
-    estimator.set_particles(inputList);
-    rho = estimator.rho();
-
-    candidate = factory->NewCandidate();
-    candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho); 
-    fRhoOutputArray->Add(candidate);
-  }
-  
+    for(itEtaRangeMap = fEtaRangeMap.begin(); itEtaRangeMap != fEtaRangeMap.end(); ++itEtaRangeMap)
+    {
+      Selector select_rapidity = SelectorAbsRapRange(itEtaRangeMap->first, itEtaRangeMap->second);
+      JetMedianBackgroundEstimator estimator(select_rapidity, *fDefinition, *fAreaDefinition);
+      estimator.set_particles(inputList);
+      rho = estimator.rho();
+
+      candidate = factory->NewCandidate();
+      candidate->Momentum.SetPtEtaPhiE(rho, 0.0, 0.0, rho);
+      candidate->Edges[0] = itEtaRangeMap->first;
+      candidate->Edges[1] = itEtaRangeMap->second;
+      fRhoOutputArray->Add(candidate);
+    }
+  }
+
   outputList.clear();
   outputList = sorted_by_pt(sequence->inclusive_jets(fJetPTMin));

modules/FastJetFinder.h

@@ -16 +16 @@
 #include "classes/DelphesModule.h"
 
-#include <vector>
+#include <map>
 
 class TObjArray;
@@ -37 +37 @@
   void Process();
   void Finish();
-  
+
 private:
 
   void *fPlugin; //!
   fastjet::JetDefinition *fDefinition; //!
-   
+
   Int_t fJetAlgorithm;
   Double_t fParameterR;
@@ -54 +54 @@
   Double_t fAdjacencyCut;
   Double_t fOverlapThreshold;
-  
+
   // --- FastJet Area method --------
-  
+
   fastjet::AreaDefinition *fAreaDefinition;
-  Int_t fAreaAlgorithm;  
-  Bool_t  fComputeRho; 
-  Double_t fRhoEtaMax;
-  
+  Int_t fAreaAlgorithm;
+  Bool_t  fComputeRho;
+
   // -- ghost based areas --
   Double_t fGhostEtaMax;
@@ -68 +67 @@
   Double_t fGridScatter;
   Double_t fPtScatter;
-  Double_t fMeanGhostPt; 
-  
+  Double_t fMeanGhostPt;
+
   // -- voronoi areas --
   Double_t fEffectiveRfact;
-  
-  
+
+  std::map< Double_t, Double_t > fEtaRangeMap; //!
+
   TIterator *fItInputArray; //!
 

modules/JetPileUpSubtractor.cc

@@ -39 +39 @@
 
 JetPileUpSubtractor::JetPileUpSubtractor() :
-  fItJetInputArray(0)
+  fItJetInputArray(0), fItRhoInputArray(0)
 {
 
@@ -63 +63 @@
 
   fRhoInputArray = ImportArray(GetString("RhoInputArray", "Rho/rho"));
+  fItRhoInputArray = fRhoInputArray->MakeIterator();
 
   // create output array(s)
@@ -74 +75 @@
 void JetPileUpSubtractor::Finish()
 {
+  if(fItRhoInputArray) delete fItRhoInputArray;
   if(fItJetInputArray) delete fItJetInputArray;
 }
@@ -83 +85 @@
   Candidate *candidate;
   TLorentzVector momentum, area;
+  Double_t eta = 0.0;
   Double_t rho = 0.0;
 
-  if(fRhoInputArray && fRhoInputArray->GetEntriesFast() > 0)
-  {
-    candidate = static_cast<Candidate*>(fRhoInputArray->At(0));
-    rho = candidate->Momentum.Pt();
-  }
+  if(!fRhoInputArray) return;
 
   // loop over all input candidates
@@ -97 +96 @@
     momentum = candidate->Momentum;
     area = candidate->Area;
+    eta = TMath::Abs(momentum.Eta());
+
+    // find rho
+    rho = 0.0;
+    while((candidate = static_cast<Candidate*>(fItRhoInputArray->Next())))
+    {
+      if(eta >= candidate->Edges[0] && eta < candidate->Edges[1])
+      {
+        rho = candidate->Momentum.Pt();
+      }
+    }  
 
     // apply pile-up correction

modules/JetPileUpSubtractor.h

@@ -36 +36 @@
 
   TIterator *fItJetInputArray; //!
+  TIterator *fItRhoInputArray; //!
 
   const TObjArray *fJetInputArray; //!