Changes in modules/PileUpJetID.cc [84edab9:6cdc544] in git

File:

• modules/PileUpJetID.cc (modified) (8 diffs)

modules/PileUpJetID.cc

@@ +1 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
 
 /** \class PileUpJetID
  *
- *  CMS PileUp Jet ID Variables
- *
- *  \author S. Zenz
+ *  CMS PileUp Jet ID Variables, based on http://cds.cern.ch/record/1581583
+ *
+ *  \author S. Zenz, December 2013
  *
  */
@@ -23 +41 @@
 #include "TRandom3.h"
 #include "TObjArray.h"
-//#include "TDatabasePDG.h"
+#include "TDatabasePDG.h"
 #include "TLorentzVector.h"
 
@@ -36 +54 @@
 
 PileUpJetID::PileUpJetID() :
-  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
+  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
 {
 
@@ -52 +70 @@
 void PileUpJetID::Init()
 {
-
   fJetPTMin = GetDouble("JetPTMin", 20.0);
   fParameterR = GetDouble("ParameterR", 0.5);
   fUseConstituents = GetInt("UseConstituents", 0);
 
-  fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
-  fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
-  fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
-  fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
-  fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
-  fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
-  fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
-
   fAverageEachTower = false; // for timing
 
@@ -72 +81 @@
   fItJetInputArray = fJetInputArray->MakeIterator();
 
-  fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/eflowTracks"));
   fItTrackInputArray = fTrackInputArray->MakeIterator();
 
-  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
+  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
   fItNeutralInputArray = fNeutralInputArray->MakeIterator();
 
+  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
+  fItVertexInputArray = fVertexInputArray->MakeIterator();
+
+  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
+
   // create output array(s)
 
   fOutputArray = ExportArray(GetString("OutputArray", "jets"));
-
-  fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
-
 }
 
@@ -90 +101 @@
 void PileUpJetID::Finish()
 {
-  //  cout << "In finish" << endl;
-
   if(fItJetInputArray) delete fItJetInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
   if(fItNeutralInputArray) delete fItNeutralInputArray;
-
+  if(fItVertexInputArray) delete fItVertexInputArray;
 }
 
@@ -104 +113 @@
   Candidate *candidate, *constituent;
   TLorentzVector momentum, area;
-
-  Candidate *trk;
+  Int_t i, nc, nn;
+  Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
+  Double_t dr, pt, pt_ann[5];
+  Double_t zvtx = 0.0;
+
+  Candidate *track;
+
+  // find z position of primary vertex
+
+  fItVertexInputArray->Reset();
+  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
+  {
+    if(!candidate->IsPU)
+    {
+      zvtx = candidate->Position.Z();
+      break;
+    }
+  }
 
   // loop over all input candidates
@@ -114 +139 @@
     area = candidate->Area;
 
-    float sumT0 = 0.;
-    float sumT1 = 0.;
-    float sumT10 = 0.;
-    float sumT20 = 0.;
-    float sumT30 = 0.;
-    float sumT40 = 0.;
-    float sumWeightsForT = 0.;
-    candidate->NTimeHits = 0;
-
-    float sumpt = 0.;
-    float sumptch = 0.;
-    float sumptchpv = 0.;
-    float sumptchpu = 0.;
-    float sumdrsqptsq = 0.;
-    float sumptsq = 0.;
-    int nc = 0;
-    int nn = 0;
-    float pt_ann[5];
-
-    for (int i = 0 ; i < 5 ; i++) {
-      pt_ann[i] = 0.;
-    }
-
-    if (fUseConstituents) {
+    sumpt = 0.0;
+    sumptch = 0.0;
+    sumptchpv = 0.0;
+    sumptchpu = 0.0;
+    sumdrsqptsq = 0.0;
+    sumptsq = 0.0;
+    nc = 0;
+    nn = 0;
+
+    for(i = 0; i < 5; ++i)
+    {
+      pt_ann[i] = 0.0;
+    }
+
+    if(fUseConstituents)
+    {
       TIter itConstituents(candidate->GetCandidates());
-      while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
-        float pt = constituent->Momentum.Pt();
-        float dr = candidate->Momentum.DeltaR(constituent->Momentum);
-        //      cout << " There exists a constituent with dr=" << dr << endl;
-        sumpt += pt;
-        sumdrsqptsq += dr*dr*pt*pt;
-        sumptsq += pt*pt;
-        if (constituent->Charge == 0) {
-          nn++;
-        } else {
-          if (constituent->IsRecoPU) {
-            sumptchpu += pt;
-          } else {
-            sumptchpv += pt;
-          }
-          sumptch += pt;
-          nc++;
-        }
-        for (int i = 0 ; i < 5 ; i++) {
-          if (dr > 0.1*i && dr < 0.1*(i+1)) {
-            pt_ann[i] += pt;
-          }
-        }
-        float tow_sumT = 0;
-        float tow_sumW = 0;
-        for (int i = 0 ; i < constituent->ECalEnergyTimePairs.size() ; i++) {
-          float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first);
-          if (fAverageEachTower) {
-            tow_sumT += w*constituent->ECalEnergyTimePairs[i].second;
-            tow_sumW += w;
-          } else {
-            sumT0 += w*constituent->ECalEnergyTimePairs[i].second;
-            sumT1 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.001);
-            sumT10 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.010);
-            sumT20 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.020);
-            sumT30 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.030);
-            sumT40 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.040);
-            sumWeightsForT += w;
-            candidate->NTimeHits++;
-          }
-        }
-        if (fAverageEachTower && tow_sumW > 0.) {
-          sumT0 += tow_sumT;
-          sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
-          sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
-          sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
-          sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
-          sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
-          sumWeightsForT += tow_sumW;
-          candidate->NTimeHits++;
-        }
-      }
-    } else {
+      while((constituent = static_cast<Candidate*>(itConstituents.Next())))
+      {
+        pt = constituent->Momentum.Pt();
+        dr = candidate->Momentum.DeltaR(constituent->Momentum);
+        sumpt += pt;
+        sumdrsqptsq += dr*dr*pt*pt;
+        sumptsq += pt*pt;
+        if(constituent->Charge == 0)
+        {
+          // neutrals
+          ++nn;
+        }
+        else
+        {
+          // charged
+          if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
+          {
+            sumptchpu += pt;
+          }
+          else
+          {
+            sumptchpv += pt;
+          }
+          sumptch += pt;
+          ++nc;
+        }
+        for(i = 0; i < 5; ++i)
+        {
+          if(dr > 0.1*i && dr < 0.1*(i + 1))
+          {
+            pt_ann[i] += pt;
+          }
+        }
+      }
+    }
+    else
+    {
       // Not using constituents, using dr
       fItTrackInputArray->Reset();
-      while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
-        if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
-          float pt = trk->Momentum.Pt();
-          sumpt += pt;
-          sumptch += pt;
-          if (trk->IsRecoPU) {
-            sumptchpu += pt;
-          } else {
-            sumptchpv += pt;
-          }
-          float dr = candidate->Momentum.DeltaR(trk->Momentum);
-          sumdrsqptsq += dr*dr*pt*pt;
-          sumptsq += pt*pt;
-          nc++;
-          for (int i = 0 ; i < 5 ; i++) {
-            if (dr > 0.1*i && dr < 0.1*(i+1)) {
-              pt_ann[i] += pt;
-            }
-          }
-        }
-      }
+      while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
+      {
+        if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
+        {
+          pt = track->Momentum.Pt();
+          sumpt += pt;
+          sumptch += pt;
+          if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
+          {
+            sumptchpu += pt;
+          }
+          else
+          {
+            sumptchpv += pt;
+          }
+          dr = candidate->Momentum.DeltaR(track->Momentum);
+          sumdrsqptsq += dr*dr*pt*pt;
+          sumptsq += pt*pt;
+          nc++;
+          for(i = 0; i < 5; ++i)
+          {
+            if(dr > 0.1*i && dr < 0.1*(i + 1))
+            {
+              pt_ann[i] += pt;
+            }
+          }
+        }
+      }
+
       fItNeutralInputArray->Reset();
-      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
-        if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
-          float pt = constituent->Momentum.Pt();
-          sumpt += pt;
-          float dr = candidate->Momentum.DeltaR(constituent->Momentum);
-          sumdrsqptsq += dr*dr*pt*pt;
-          sumptsq += pt*pt;
-          nn++;
-          for (int i = 0 ; i < 5 ; i++) {
-            if (dr > 0.1*i && dr < 0.1*(i+1)) {
-            pt_ann[i] += pt;
-            }
-          }
-        }
-      }
-    }
-
-    if (sumptch > 0.) {
-      candidate->Beta = sumptchpv/sumptch;
-      candidate->BetaStar = sumptchpu/sumptch;
-    } else {
-      candidate->Beta = -999.;
-      candidate->BetaStar = -999.;
-    }
-    if (sumptsq > 0.) {
+      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
+      {
+        if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
+        {
+          pt = constituent->Momentum.Pt();
+          sumpt += pt;
+          dr = candidate->Momentum.DeltaR(constituent->Momentum);
+          sumdrsqptsq += dr*dr*pt*pt;
+          sumptsq += pt*pt;
+          nn++;
+          for(i = 0; i < 5; ++i)
+          {
+            if(dr > 0.1*i && dr < 0.1*(i + 1))
+            {
+              pt_ann[i] += pt;
+            }
+          }
+        }
+      }
+    }
+
+    if(sumptch > 0.0)
+    {
+      candidate->Beta = sumptchpu/sumptch;
+      candidate->BetaStar = sumptchpv/sumptch;
+    }
+    else
+    {
+      candidate->Beta = -999.0;
+      candidate->BetaStar = -999.0;
+    }
+    if(sumptsq > 0.0)
+    {
       candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
-    } else {
-      candidate->MeanSqDeltaR = -999.;
+    }
+    else
+    {
+      candidate->MeanSqDeltaR = -999.0;
     }
     candidate->NCharged = nc;
     candidate->NNeutrals = nn;
-    if (sumpt > 0.) {
+    if(sumpt > 0.0)
+    {
       candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
-      for (int i = 0 ; i < 5 ; i++) {
+      for(i = 0; i < 5; ++i)
+      {
         candidate->FracPt[i] = pt_ann[i]/sumpt;
       }
-    } else {
-      candidate->PTD = -999.;
-      for (int i = 0 ; i < 5 ; i++) {
-        candidate->FracPt[i] = -999.;
+    }
+    else
+    {
+      candidate->PTD = -999.0;
+      for(i = 0; i < 5; ++i)
+      {
+        candidate->FracPt[i] = -999.0;
       }
     }
 
     fOutputArray->Add(candidate);
-
-    // New stuff
-    /*
-    fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
-    fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
-    fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
-    fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
-    fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
-    */
-
-    bool passId = false;
-    if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
-      if (fabs(candidate->Momentum.Eta())<1.5) {
-        passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
-      } else if (fabs(candidate->Momentum.Eta())<4.0) {
-        passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
-      } else {
-        passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
-      }
-    }
-
-    //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt() 
-    //   << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
-
-    if (passId) {
-      if (fUseConstituents) {
-        TIter itConstituents(candidate->GetCandidates());
-        while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
-          if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
-            fNeutralsInPassingJets->Add(constituent);
-            //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
-          }
-        }
-      } else { // use DeltaR
-        fItNeutralInputArray->Reset();
-        while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
-          if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
-            fNeutralsInPassingJets->Add(constituent);
-            //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
-          }
-        }
-      }
-    }
-
-
   }
 }
 
 //------------------------------------------------------------------------------
+