Changeset fc6bce2 in git for readers

Timestamp:

Jun 1, 2021, 8:13:13 PM (4 years ago)

Author:

Pavel Demin <pavel.demin@…>

Branches:

master

Children:

bfc6722

Parents:

4006893

Message:

remove 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"
@@ -41 +35 @@
 #include "classes/DelphesClasses.h"
 #include "classes/DelphesFactory.h"
+#include "classes/DelphesHepMC3Reader.h"
 #include "modules/Delphes.h"
 
@@ -47 +42 @@
 #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;
-  }
-}
 
 //---------------------------------------------------------------------------
@@ -346 +59 @@
   char appName[] = "DelphesHepMC3";
   stringstream message;
-  ifstream inputFile;
-  filebuf *inputBuffer;
+  FILE *inputFile = 0;
   TFile *outputFile = 0;
   TStopwatch readStopWatch, procStopWatch;
@@ -355 +67 @@
   Delphes *modularDelphes = 0;
   DelphesFactory *factory = 0;
-  TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
-  ReaderAscii *reader = 0;
-  GenEvent event;
+  TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
+  DelphesHepMC3Reader *reader = 0;
   Int_t i, maxEvents, skipEvents;
   Long64_t length, eventCounter;
@@ -421 +132 @@
     partonOutputArray = modularDelphes->ExportArray("partons");
 
-    inputBuffer = inputFile.rdbuf();
-    reader = new ReaderAscii(inputFile);
+    reader = new DelphesHepMC3Reader;
 
     modularDelphes->InitTask();
@@ -434 +144 @@
       {
         cout << "** Reading standard input" << endl;
-        inputFile.ios::rdbuf(cin.rdbuf());
+        inputFile = stdin;
         length = -1;
       }
@@ -440 +150 @@
       {
         cout << "** Reading " << argv[i] << endl;
-        inputFile.ios::rdbuf(inputBuffer);
-        inputFile.open(argv[i]);
-
-        if(inputFile.fail())
+        inputFile = fopen(argv[i], "r");
+
+        if(inputFile == NULL)
         {
           message << "can't open " << argv[i];
@@ -449 +158 @@
         }
 
-        inputFile.seekg(0, ios::end);
-        length = inputFile.tellg();
-        inputFile.seekg(0, ios::beg);
+        fseek(inputFile, 0L, SEEK_END);
+        length = ftello(inputFile);
+        fseek(inputFile, 0L, SEEK_SET);
 
         if(length <= 0)
         {
-          inputFile.close();
-          inputFile.clear();
+          fclose(inputFile);
           ++i;
           continue;
@@ -462 +170 @@
       }
 
+      reader->SetInputFile(inputFile);
+
       ExRootProgressBar progressBar(length);
 
-      reader->skip(skipEvents);
-      progressBar.Update(inputFile.tellg(), skipEvents);
-
-      eventCounter = skipEvents;
+      // Loop over all objects
+      eventCounter = 0;
+      treeWriter->Clear();
       modularDelphes->Clear();
-      treeWriter->Clear();
-      event.clear();
-      gMotherMap.clear();
-      gDaughterMap.clear();
+      reader->Clear();
       readStopWatch.Start();
-      reader->read_event(event);
-      while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && !reader->failed() && !interrupted)
+      while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && reader->ReadBlock(factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray) && !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);
+        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);
       }
 
-      progressBar.Update(length, eventCounter, kTRUE);
+      fseek(inputFile, 0L, SEEK_END);
+      progressBar.Update(ftello(inputFile), eventCounter, kTRUE);
       progressBar.Finish();
 
-      if(length > 0) inputFile.close();
+      if(inputFile != stdin) fclose(inputFile);
 
       ++i;