Changeset 95a917c in git for readers/DelphesHepMC3.cpp

Timestamp:

May 28, 2021, 6:07:25 PM (3 years ago)

Author:

Pavel Demin <pavel.demin@…>

Branches:

master

Children:

302624f

Parents:

91eef4d

Message:

add HepMC3 library

File:

• readers/DelphesHepMC3.cpp (modified) (10 diffs)

readers/DelphesHepMC3.cpp

@@ -18 +18 @@
 
 #include <iostream>
+#include <fstream>
+#include <memory>
 #include <sstream>
 #include <stdexcept>
 
+#include <map>
+
 #include <signal.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 #include "TApplication.h"
@@ -35 +41 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
-#include "classes/DelphesHepMC3Reader.h"
 #include "modules/Delphes.h"
 
@@ -42 +47 @@
 #include "ExRootAnalysis/ExRootTreeWriter.h"
 
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/GenCrossSection.h"
+#include "HepMC3/GenPdfInfo.h"
+#include "HepMC3/GenParticle.h"
+#include "HepMC3/GenVertex.h"
+#include "HepMC3/Units.h"
+
 using namespace std;
+using namespace HepMC3;
+
+map<Int_t, pair<Int_t, Int_t> > gMotherMap;
+map<Int_t, pair<Int_t, Int_t> > gDaughterMap;
+
+//---------------------------------------------------------------------------
+
+void AnalyzeParticle(Bool_t in, Int_t counter,
+  Double_t momentumCoefficient,
+  Double_t positionCoefficient,
+  shared_ptr<HepMC3::GenVertex> vertex,
+  shared_ptr<HepMC3::GenParticle> particle,
+  DelphesFactory *factory,
+  TObjArray *allParticleOutputArray,
+  TObjArray *stableParticleOutputArray,
+  TObjArray *partonOutputArray)
+{
+  map<Int_t, pair<Int_t, Int_t> >::iterator itMotherMap;
+  map<Int_t, pair<Int_t, Int_t> >::iterator itDaughterMap;
+
+  Candidate *candidate;
+  TDatabasePDG *pdg;
+  TParticlePDG *pdgParticle;
+  Int_t pdgCode;
+
+  Int_t pid, status, inVertexCode, outVertexCode;
+  Double_t px, py, pz, e, mass;
+  Double_t x, y, z, t;
+
+  pdg = TDatabasePDG::Instance();
+
+  candidate = factory->NewCandidate();
+
+  pid = particle->pid();
+  px = particle->momentum().px();
+  py = particle->momentum().py();
+  pz = particle->momentum().pz();
+  e = particle->momentum().e();
+  mass = particle->generated_mass();
+  x = vertex->position().x();
+  y = vertex->position().y();
+  z = vertex->position().z();
+  t = vertex->position().t();
+  status = particle->status();
+
+  outVertexCode = vertex->id();
+  inVertexCode = particle->end_vertex() ? particle->end_vertex()->id() : 0;
+
+  candidate->PID = pid;
+  pdgCode = TMath::Abs(pid);
+
+  candidate->Status = status;
+
+  pdgParticle = pdg->GetParticle(pid);
+  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
+  candidate->Mass = mass;
+
+  candidate->Momentum.SetPxPyPzE(px, py, pz, e);
+  if(momentumCoefficient != 1.0)
+  {
+    candidate->Momentum *= momentumCoefficient;
+  }
+
+  candidate->M2 = 1;
+  candidate->D2 = 1;
+
+  if(in)
+  {
+    candidate->M1 = 1;
+    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
+  }
+  else
+  {
+    candidate->M1 = outVertexCode;
+    candidate->Position.SetXYZT(x, y, z, t);
+    if(positionCoefficient != 1.0)
+    {
+      candidate->Position *= positionCoefficient;
+    }
+
+    itDaughterMap = gDaughterMap.find(outVertexCode);
+    if(itDaughterMap == gDaughterMap.end())
+    {
+      gDaughterMap[outVertexCode] = make_pair(counter, counter);
+    }
+    else
+    {
+      itDaughterMap->second.second = counter;
+    }
+  }
+
+  if(inVertexCode < 0)
+  {
+    candidate->D1 = inVertexCode;
+
+    itMotherMap = gMotherMap.find(inVertexCode);
+    if(itMotherMap == gMotherMap.end())
+    {
+      gMotherMap[inVertexCode] = make_pair(counter, -1);
+    }
+    else
+    {
+      itMotherMap->second.second = counter;
+    }
+  }
+  else
+  {
+    candidate->D1 = 1;
+  }
+
+  allParticleOutputArray->Add(candidate);
+
+  if(!pdgParticle) return;
+
+  if(status == 1)
+  {
+    stableParticleOutputArray->Add(candidate);
+  }
+  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
+  {
+    partonOutputArray->Add(candidate);
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void AnalyzeEvent(GenEvent &event, ExRootTreeBranch *branchEvent, DelphesFactory *factory,
+  TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray,
+  TObjArray *partonOutputArray, TStopwatch *readStopWatch, TStopwatch *procStopWatch)
+{
+  Int_t i, counter;
+  map<Int_t, pair<Int_t, Int_t> >::iterator itMotherMap;
+  map<Int_t, pair<Int_t, Int_t> >::iterator itDaughterMap;
+
+  HepMCEvent *element;
+  Candidate *candidate;
+  Double_t momentumCoefficient, positionCoefficient;
+
+  shared_ptr<IntAttribute> processID = event.attribute<IntAttribute>("signal_process_id");
+  shared_ptr<IntAttribute> mpi = event.attribute<IntAttribute>("mpi");
+  shared_ptr<DoubleAttribute> scale = event.attribute<DoubleAttribute>("event_scale");
+  shared_ptr<DoubleAttribute> alphaQED = event.attribute<DoubleAttribute>("alphaQED");
+  shared_ptr<DoubleAttribute> alphaQCD = event.attribute<DoubleAttribute>("alphaQCD");
+
+  shared_ptr<GenCrossSection> cs = event.attribute<GenCrossSection>("GenCrossSection");
+  shared_ptr<GenPdfInfo> pdf = event.attribute<GenPdfInfo>("GenPdfInfo");
+
+  element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
+
+  element->Number = event.event_number();
+
+  element->ProcessID = processID ? processID->value() : 0;
+  element->MPI = mpi ? mpi->value() : 0;
+  element->Weight = event.weights().size() > 0 ? event.weights()[0] : 1.0;
+  element->Scale = scale ? scale->value() : 0.0;
+  element->AlphaQED = alphaQED ? alphaQED->value() : 0.0;
+  element->AlphaQCD = alphaQCD ? alphaQCD->value() : 0.0;
+
+  if(cs)
+  {
+    element->CrossSection = cs->xsec();
+    element->CrossSectionError = cs->xsec_err();
+  }
+  else
+  {
+    element->CrossSection = 0.0;
+    element->CrossSectionError = 0.0;;
+  }
+
+  if(pdf)
+  {
+    element->ID1 = pdf->parton_id[0];
+    element->ID2 = pdf->parton_id[1];
+    element->X1 = pdf->x[0];
+    element->X2 = pdf->x[1];
+    element->ScalePDF = pdf->scale;
+    element->PDF1 = pdf->xf[0];
+    element->PDF2 = pdf->xf[1];
+  }
+  else
+  {
+    element->ID1 = 0;
+    element->ID2 = 0;
+    element->X1 = 0.0;
+    element->X2 = 0.0;
+    element->ScalePDF = 0.0;
+    element->PDF1 = 0.0;
+    element->PDF2 = 0.0;
+  }
+
+  element->ReadTime = readStopWatch->RealTime();
+  element->ProcTime = procStopWatch->RealTime();
+
+  momentumCoefficient = (event.momentum_unit() == Units::MEV) ? 0.001 : 1.0;
+  positionCoefficient = (event.length_unit() == Units::MM) ? 1.0 : 10.0;
+
+  counter = 0;
+  for(auto vertex: event.vertices())
+  {
+    for(auto particle: vertex->particles_in())
+    {
+      if(!particle->production_vertex() || particle->production_vertex()->id() == 0)
+      {
+        AnalyzeParticle(kTRUE, counter, momentumCoefficient, positionCoefficient, vertex, particle,
+          factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
+        ++counter;
+      }
+    }
+    for(auto particle: vertex->particles_out())
+    {
+      AnalyzeParticle(kFALSE, counter, momentumCoefficient, positionCoefficient, vertex, particle,
+        factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
+      ++counter;
+    }
+  }
+
+  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
+  {
+    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
+
+    if(candidate->M1 > 0)
+    {
+      candidate->M1 = -1;
+      candidate->M2 = -1;
+    }
+    else
+    {
+      itMotherMap = gMotherMap.find(candidate->M1);
+      if(itMotherMap == gMotherMap.end())
+      {
+        candidate->M1 = -1;
+        candidate->M2 = -1;
+      }
+      else
+      {
+        candidate->M1 = itMotherMap->second.first;
+        candidate->M2 = itMotherMap->second.second;
+      }
+    }
+    if(candidate->D1 > 0)
+    {
+      candidate->D1 = -1;
+      candidate->D2 = -1;
+    }
+    else
+    {
+      itDaughterMap = gDaughterMap.find(candidate->D1);
+      if(itDaughterMap == gDaughterMap.end())
+      {
+        candidate->D1 = -1;
+        candidate->D2 = -1;
+      }
+      else
+      {
+        candidate->D1 = itDaughterMap->second.first;
+        candidate->D2 = itDaughterMap->second.second;
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void AnalyzeWeight(GenEvent &event, ExRootTreeBranch *branchWeight)
+{
+  Weight *element;
+
+  for(auto weight: event.weights())
+  {
+    element = static_cast<Weight *>(branchWeight->NewEntry());
+
+    element->Weight = weight;
+  }
+}
 
 //---------------------------------------------------------------------------
@@ -59 +346 @@
   char appName[] = "DelphesHepMC3";
   stringstream message;
-  FILE *inputFile = 0;
+  ifstream inputFile;
+  filebuf *inputBuffer;
   TFile *outputFile = 0;
   TStopwatch readStopWatch, procStopWatch;
@@ -67 +355 @@
   Delphes *modularDelphes = 0;
   DelphesFactory *factory = 0;
-  TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
-  DelphesHepMC3Reader *reader = 0;
+  TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
+  ReaderAscii *reader = 0;
+  GenEvent event;
   Int_t i, maxEvents, skipEvents;
   Long64_t length, eventCounter;
@@ -132 +421 @@
     partonOutputArray = modularDelphes->ExportArray("partons");
 
-    reader = new DelphesHepMC3Reader;
+    inputBuffer = inputFile.rdbuf();
+    reader = new ReaderAscii(inputFile);
 
     modularDelphes->InitTask();
@@ -144 +434 @@
       {
         cout << "** Reading standard input" << endl;
-        inputFile = stdin;
+        inputFile.ios::rdbuf(cin.rdbuf());
         length = -1;
       }
@@ -150 +440 @@
       {
         cout << "** Reading " << argv[i] << endl;
-        inputFile = fopen(argv[i], "r");
-
-        if(inputFile == NULL)
+        inputFile.ios::rdbuf(inputBuffer);
+        inputFile.open(argv[i]);
+
+        if(inputFile.fail())
         {
           message << "can't open " << argv[i];
@@ -158 +449 @@
         }
 
-        fseek(inputFile, 0L, SEEK_END);
-        length = ftello(inputFile);
-        fseek(inputFile, 0L, SEEK_SET);
+        inputFile.seekg(0, ios::end);
+        length = inputFile.tellg();
+        inputFile.seekg(0, ios::beg);
 
         if(length <= 0)
         {
-          fclose(inputFile);
+          inputFile.close();
+          inputFile.clear();
           ++i;
           continue;
@@ -170 +462 @@
       }
 
-      reader->SetInputFile(inputFile);
-
       ExRootProgressBar progressBar(length);
 
-      // Loop over all objects
-      eventCounter = 0;
+      reader->skip(skipEvents);
+      progressBar.Update(inputFile.tellg(), skipEvents);
+
+      eventCounter = skipEvents;
+      modularDelphes->Clear();
       treeWriter->Clear();
-      modularDelphes->Clear();
-      reader->Clear();
+      event.clear();
+      gMotherMap.clear();
+      gDaughterMap.clear();
       readStopWatch.Start();
-      while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && reader->ReadBlock(factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray) && !interrupted)
-      {
-        if(reader->EventReady())
-        {
-          ++eventCounter;
-
-          readStopWatch.Stop();
-
-          if(eventCounter > skipEvents)
-          {
-            procStopWatch.Start();
-            modularDelphes->ProcessTask();
-            procStopWatch.Stop();
-
-            reader->AnalyzeEvent(branchEvent, eventCounter, &readStopWatch, &procStopWatch);
-            reader->AnalyzeWeight(branchWeight);
-
-            treeWriter->Fill();
-
-            treeWriter->Clear();
-          }
-
-          modularDelphes->Clear();
-          reader->Clear();
-
-          readStopWatch.Start();
-        }
-        progressBar.Update(ftello(inputFile), eventCounter);
-      }
-
-      fseek(inputFile, 0L, SEEK_END);
-      progressBar.Update(ftello(inputFile), eventCounter, kTRUE);
+      reader->read_event(event);
+      while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && !reader->failed() && !interrupted)
+      {
+        readStopWatch.Stop();
+
+        ++eventCounter;
+
+        procStopWatch.Start();
+        AnalyzeEvent(event, branchEvent, factory,
+          allParticleOutputArray, stableParticleOutputArray,
+          partonOutputArray, &readStopWatch, &procStopWatch);
+        modularDelphes->ProcessTask();
+        AnalyzeWeight(event, branchWeight);
+        procStopWatch.Stop();
+
+        treeWriter->Fill();
+
+        modularDelphes->Clear();
+        treeWriter->Clear();
+        event.clear();
+        gMotherMap.clear();
+        gDaughterMap.clear();
+
+        progressBar.Update(inputFile.tellg(), eventCounter);
+
+        readStopWatch.Start();
+        reader->read_event(event);
+      }
+
+      progressBar.Update(length, eventCounter, kTRUE);
       progressBar.Finish();
 
-      if(inputFile != stdin) fclose(inputFile);
+      if(length > 0) inputFile.close();
 
       ++i;