Changes in / [a0f2226:a7c9002] in git

File:

• readers/DelphesCMSFWLite.cpp (modified) (12 diffs)

readers/DelphesCMSFWLite.cpp

@@ -53 +53 @@
 #include "DataFormats/FWLite/interface/Handle.h"
 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/PatCandidates/interface/PackedGenParticle.h"
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
@@ -63 +64 @@
 
 void ConvertInput(fwlite::Event &event, Long64_t eventCounter,
-  ExRootTreeBranch *branchEvent, ExRootTreeBranch *branchRwgt,
-  DelphesFactory *factory, TObjArray *allParticleOutputArray,
-  TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray, Bool_t firstEvent)
+    ExRootTreeBranch *branchEvent, ExRootTreeBranch *branchRwgt,
+    DelphesFactory *factory, TObjArray *allParticleOutputArray,
+    TObjArray *stableParticleOutputArray, TObjArray *partonOutputArray, Bool_t firstEvent)
 {
 
@@ -71 +72 @@
   fwlite::Handle< LHEEventProduct > handleLHEEvent;
   fwlite::Handle< vector< reco::GenParticle > > handleParticle;
+  fwlite::Handle< vector< pat::PackedGenParticle > > handlePackedParticle;
 
   vector< reco::GenParticle >::const_iterator itParticle;
+  vector< pat::PackedGenParticle >::const_iterator itPackedParticle;
 
   vector< const reco::Candidate * > vectorCandidate;
@@ -96 +99 @@
     handleParticle.getByLabel(event, "genParticles");
   }
-  else if (!((handleParticle.getBranchNameFor(event, "prunedGenParticles")).empty()))
+  else if (!((handlePackedParticle.getBranchNameFor(event, "packedGenParticles")).empty()) && !((handleParticle.getBranchNameFor(event,"prunedGenParticles")).empty()))
   {
     handleParticle.getByLabel(event, "prunedGenParticles");
+    handlePackedParticle.getByLabel(event, "packedGenParticles");
   }
   else
@@ -107 +111 @@
 
   Bool_t foundLHE = !((handleLHEEvent.getBranchNameFor(event, "source")).empty()) || !((handleLHEEvent.getBranchNameFor(event, "externalLHEProducer")).empty());
+  Bool_t isMiniAOD = !((handlePackedParticle.getBranchNameFor(event, "packedGenParticles")).empty()) && ((handleParticle.getBranchNameFor(event, "genParticles")).empty()) ;
 
   HepMCEvent *element;
@@ -146 +151 @@
     const vector< gen::WeightsInfo > &vectorWeightsInfo = handleLHEEvent->weights();
     vector< gen::WeightsInfo >::const_iterator itWeightsInfo;
-    
+
     for(itWeightsInfo = vectorWeightsInfo.begin(); itWeightsInfo != vectorWeightsInfo.end(); ++itWeightsInfo)
     {
@@ -158 +163 @@
   for(itParticle = handleParticle->begin(); itParticle != handleParticle->end(); ++itParticle)
   {
-    vectorCandidate.push_back(&*itParticle);
+    const reco::GenParticle &particle = *itParticle;
+    if( !isMiniAOD || particle.status() != 1 ) vectorCandidate.push_back(&*itParticle);
   }
 
@@ -164 +170 @@
   {
     const reco::GenParticle &particle = *itParticle;
+
+    pid = particle.pdgId();
+    status = particle.status();
+    if( isMiniAOD && particle.status() == 1 ) continue;
+    px = particle.px(); py = particle.py(); pz = particle.pz(); e = particle.energy(); mass = particle.mass();
+    x = particle.vx(); y = particle.vy(); z = particle.vz();
+
+    candidate = factory->NewCandidate();
+
+    candidate->PID = pid;
+    pdgCode = TMath::Abs(candidate->PID);
+
+    candidate->Status = status;
+
+    if(particle.mother())
+    {
+      itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother());
+      if(itCandidate != vectorCandidate.end()) candidate->M1 = distance(vectorCandidate.begin(), itCandidate);
+    }
+
+    itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(0));
+    if(itCandidate != vectorCandidate.end()) candidate->D1 = distance(vectorCandidate.begin(), itCandidate);
+
+    itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.daughter(particle.numberOfDaughters() - 1));
+    if(itCandidate != vectorCandidate.end()) candidate->D2 = distance(vectorCandidate.begin(), itCandidate);
+
+    pdgParticle = pdg->GetParticle(pid);
+    candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
+    candidate->Mass = mass;
+
+    candidate->Momentum.SetPxPyPzE(px, py, pz, e);
+
+    candidate->Position.SetXYZT(x*10.0, y*10.0, z*10.0, 0.0);
+
+    allParticleOutputArray->Add(candidate);
+
+    if(!pdgParticle) continue;
+
+    if( status == 1)
+    {
+      // Prevent duplicated particle.
+      if ( !isMiniAOD ) stableParticleOutputArray->Add(candidate);
+    }
+    else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
+    {
+      partonOutputArray->Add(candidate);
+    }
+  }
+  if ( !isMiniAOD) return ;
+  // For MiniAOD sample,
+  // Only status==1 particles are saved to packedGenParticles.
+  for(itPackedParticle = handlePackedParticle->begin(); itPackedParticle != handlePackedParticle->end(); ++itPackedParticle)
+  {
+    vectorCandidate.push_back(&*itPackedParticle);
+  }
+
+  for(itPackedParticle = handlePackedParticle->begin(); itPackedParticle != handlePackedParticle->end(); ++itPackedParticle)
+  {
+    const pat::PackedGenParticle &particle = *itPackedParticle;
 
     pid = particle.pdgId();
@@ -177 +242 @@
     candidate->Status = status;
 
-    if(particle.mother())
-    {
-      itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother());
+    if(particle.mother(0))
+    {
+      itCandidate = find(vectorCandidate.begin(), vectorCandidate.end(), particle.mother(0));
       if(itCandidate != vectorCandidate.end()) candidate->M1 = distance(vectorCandidate.begin(), itCandidate);
     }
@@ -204 +269 @@
     {
       stableParticleOutputArray->Add(candidate);
-    }
-    else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
-    {
-      partonOutputArray->Add(candidate);
     }
   }
@@ -258 +319 @@
 
   FWLiteEnabler::enable();
-  
+
   try
   {
@@ -317 +378 @@
       {
         ConvertInput(event, eventCounter, branchEvent, branchRwgt, factory,
-          allParticleOutputArray, stableParticleOutputArray, partonOutputArray, firstEvent);
+            allParticleOutputArray, stableParticleOutputArray, partonOutputArray, firstEvent);
         modularDelphes->ProcessTask();
-          
+
         firstEvent = kFALSE;