Fork me on GitHub

Changeset 95a917c in git for readers


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • readers/DelphesHepMC3.cpp

    r91eef4d r95a917c  
    1818
    1919#include <iostream>
     20#include <fstream>
     21#include <memory>
    2022#include <sstream>
    2123#include <stdexcept>
    2224
     25#include <map>
     26
    2327#include <signal.h>
     28#include <stdio.h>
     29#include <stdlib.h>
    2430
    2531#include "TApplication.h"
     
    3541#include "classes/DelphesClasses.h"
    3642#include "classes/DelphesFactory.h"
    37 #include "classes/DelphesHepMC3Reader.h"
    3843#include "modules/Delphes.h"
    3944
     
    4247#include "ExRootAnalysis/ExRootTreeWriter.h"
    4348
     49#include "HepMC3/ReaderAscii.h"
     50#include "HepMC3/GenEvent.h"
     51#include "HepMC3/GenCrossSection.h"
     52#include "HepMC3/GenPdfInfo.h"
     53#include "HepMC3/GenParticle.h"
     54#include "HepMC3/GenVertex.h"
     55#include "HepMC3/Units.h"
     56
    4457using namespace std;
     58using namespace HepMC3;
     59
     60map<Int_t, pair<Int_t, Int_t> > gMotherMap;
     61map<Int_t, pair<Int_t, Int_t> > gDaughterMap;
     62
     63//---------------------------------------------------------------------------
     64
     65void AnalyzeParticle(Bool_t in, Int_t counter,
     66  Double_t momentumCoefficient,
     67  Double_t positionCoefficient,
     68  shared_ptr<HepMC3::GenVertex> vertex,
     69  shared_ptr<HepMC3::GenParticle> particle,
     70  DelphesFactory *factory,
     71  TObjArray *allParticleOutputArray,
     72  TObjArray *stableParticleOutputArray,
     73  TObjArray *partonOutputArray)
     74{
     75  map<Int_t, pair<Int_t, Int_t> >::iterator itMotherMap;
     76  map<Int_t, pair<Int_t, Int_t> >::iterator itDaughterMap;
     77
     78  Candidate *candidate;
     79  TDatabasePDG *pdg;
     80  TParticlePDG *pdgParticle;
     81  Int_t pdgCode;
     82
     83  Int_t pid, status, inVertexCode, outVertexCode;
     84  Double_t px, py, pz, e, mass;
     85  Double_t x, y, z, t;
     86
     87  pdg = TDatabasePDG::Instance();
     88
     89  candidate = factory->NewCandidate();
     90
     91  pid = particle->pid();
     92  px = particle->momentum().px();
     93  py = particle->momentum().py();
     94  pz = particle->momentum().pz();
     95  e = particle->momentum().e();
     96  mass = particle->generated_mass();
     97  x = vertex->position().x();
     98  y = vertex->position().y();
     99  z = vertex->position().z();
     100  t = vertex->position().t();
     101  status = particle->status();
     102
     103  outVertexCode = vertex->id();
     104  inVertexCode = particle->end_vertex() ? particle->end_vertex()->id() : 0;
     105
     106  candidate->PID = pid;
     107  pdgCode = TMath::Abs(pid);
     108
     109  candidate->Status = status;
     110
     111  pdgParticle = pdg->GetParticle(pid);
     112  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
     113  candidate->Mass = mass;
     114
     115  candidate->Momentum.SetPxPyPzE(px, py, pz, e);
     116  if(momentumCoefficient != 1.0)
     117  {
     118    candidate->Momentum *= momentumCoefficient;
     119  }
     120
     121  candidate->M2 = 1;
     122  candidate->D2 = 1;
     123
     124  if(in)
     125  {
     126    candidate->M1 = 1;
     127    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
     128  }
     129  else
     130  {
     131    candidate->M1 = outVertexCode;
     132    candidate->Position.SetXYZT(x, y, z, t);
     133    if(positionCoefficient != 1.0)
     134    {
     135      candidate->Position *= positionCoefficient;
     136    }
     137
     138    itDaughterMap = gDaughterMap.find(outVertexCode);
     139    if(itDaughterMap == gDaughterMap.end())
     140    {
     141      gDaughterMap[outVertexCode] = make_pair(counter, counter);
     142    }
     143    else
     144    {
     145      itDaughterMap->second.second = counter;
     146    }
     147  }
     148
     149  if(inVertexCode < 0)
     150  {
     151    candidate->D1 = inVertexCode;
     152
     153    itMotherMap = gMotherMap.find(inVertexCode);
     154    if(itMotherMap == gMotherMap.end())
     155    {
     156      gMotherMap[inVertexCode] = make_pair(counter, -1);
     157    }
     158    else
     159    {
     160      itMotherMap->second.second = counter;
     161    }
     162  }
     163  else
     164  {
     165    candidate->D1 = 1;
     166  }
     167
     168  allParticleOutputArray->Add(candidate);
     169
     170  if(!pdgParticle) return;
     171
     172  if(status == 1)
     173  {
     174    stableParticleOutputArray->Add(candidate);
     175  }
     176  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
     177  {
     178    partonOutputArray->Add(candidate);
     179  }
     180}
     181
     182//---------------------------------------------------------------------------
     183
     184void AnalyzeEvent(GenEvent &event, ExRootTreeBranch *branchEvent, DelphesFactory *factory,
     185  TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray,
     186  TObjArray *partonOutputArray, TStopwatch *readStopWatch, TStopwatch *procStopWatch)
     187{
     188  Int_t i, counter;
     189  map<Int_t, pair<Int_t, Int_t> >::iterator itMotherMap;
     190  map<Int_t, pair<Int_t, Int_t> >::iterator itDaughterMap;
     191
     192  HepMCEvent *element;
     193  Candidate *candidate;
     194  Double_t momentumCoefficient, positionCoefficient;
     195
     196  shared_ptr<IntAttribute> processID = event.attribute<IntAttribute>("signal_process_id");
     197  shared_ptr<IntAttribute> mpi = event.attribute<IntAttribute>("mpi");
     198  shared_ptr<DoubleAttribute> scale = event.attribute<DoubleAttribute>("event_scale");
     199  shared_ptr<DoubleAttribute> alphaQED = event.attribute<DoubleAttribute>("alphaQED");
     200  shared_ptr<DoubleAttribute> alphaQCD = event.attribute<DoubleAttribute>("alphaQCD");
     201
     202  shared_ptr<GenCrossSection> cs = event.attribute<GenCrossSection>("GenCrossSection");
     203  shared_ptr<GenPdfInfo> pdf = event.attribute<GenPdfInfo>("GenPdfInfo");
     204
     205  element = static_cast<HepMCEvent *>(branchEvent->NewEntry());
     206
     207  element->Number = event.event_number();
     208
     209  element->ProcessID = processID ? processID->value() : 0;
     210  element->MPI = mpi ? mpi->value() : 0;
     211  element->Weight = event.weights().size() > 0 ? event.weights()[0] : 1.0;
     212  element->Scale = scale ? scale->value() : 0.0;
     213  element->AlphaQED = alphaQED ? alphaQED->value() : 0.0;
     214  element->AlphaQCD = alphaQCD ? alphaQCD->value() : 0.0;
     215
     216  if(cs)
     217  {
     218    element->CrossSection = cs->xsec();
     219    element->CrossSectionError = cs->xsec_err();
     220  }
     221  else
     222  {
     223    element->CrossSection = 0.0;
     224    element->CrossSectionError = 0.0;;
     225  }
     226
     227  if(pdf)
     228  {
     229    element->ID1 = pdf->parton_id[0];
     230    element->ID2 = pdf->parton_id[1];
     231    element->X1 = pdf->x[0];
     232    element->X2 = pdf->x[1];
     233    element->ScalePDF = pdf->scale;
     234    element->PDF1 = pdf->xf[0];
     235    element->PDF2 = pdf->xf[1];
     236  }
     237  else
     238  {
     239    element->ID1 = 0;
     240    element->ID2 = 0;
     241    element->X1 = 0.0;
     242    element->X2 = 0.0;
     243    element->ScalePDF = 0.0;
     244    element->PDF1 = 0.0;
     245    element->PDF2 = 0.0;
     246  }
     247
     248  element->ReadTime = readStopWatch->RealTime();
     249  element->ProcTime = procStopWatch->RealTime();
     250
     251  momentumCoefficient = (event.momentum_unit() == Units::MEV) ? 0.001 : 1.0;
     252  positionCoefficient = (event.length_unit() == Units::MM) ? 1.0 : 10.0;
     253
     254  counter = 0;
     255  for(auto vertex: event.vertices())
     256  {
     257    for(auto particle: vertex->particles_in())
     258    {
     259      if(!particle->production_vertex() || particle->production_vertex()->id() == 0)
     260      {
     261        AnalyzeParticle(kTRUE, counter, momentumCoefficient, positionCoefficient, vertex, particle,
     262          factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
     263        ++counter;
     264      }
     265    }
     266    for(auto particle: vertex->particles_out())
     267    {
     268      AnalyzeParticle(kFALSE, counter, momentumCoefficient, positionCoefficient, vertex, particle,
     269        factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
     270      ++counter;
     271    }
     272  }
     273
     274  for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
     275  {
     276    candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
     277
     278    if(candidate->M1 > 0)
     279    {
     280      candidate->M1 = -1;
     281      candidate->M2 = -1;
     282    }
     283    else
     284    {
     285      itMotherMap = gMotherMap.find(candidate->M1);
     286      if(itMotherMap == gMotherMap.end())
     287      {
     288        candidate->M1 = -1;
     289        candidate->M2 = -1;
     290      }
     291      else
     292      {
     293        candidate->M1 = itMotherMap->second.first;
     294        candidate->M2 = itMotherMap->second.second;
     295      }
     296    }
     297    if(candidate->D1 > 0)
     298    {
     299      candidate->D1 = -1;
     300      candidate->D2 = -1;
     301    }
     302    else
     303    {
     304      itDaughterMap = gDaughterMap.find(candidate->D1);
     305      if(itDaughterMap == gDaughterMap.end())
     306      {
     307        candidate->D1 = -1;
     308        candidate->D2 = -1;
     309      }
     310      else
     311      {
     312        candidate->D1 = itDaughterMap->second.first;
     313        candidate->D2 = itDaughterMap->second.second;
     314      }
     315    }
     316  }
     317}
     318
     319//---------------------------------------------------------------------------
     320
     321void AnalyzeWeight(GenEvent &event, ExRootTreeBranch *branchWeight)
     322{
     323  Weight *element;
     324
     325  for(auto weight: event.weights())
     326  {
     327    element = static_cast<Weight *>(branchWeight->NewEntry());
     328
     329    element->Weight = weight;
     330  }
     331}
    45332
    46333//---------------------------------------------------------------------------
     
    59346  char appName[] = "DelphesHepMC3";
    60347  stringstream message;
    61   FILE *inputFile = 0;
     348  ifstream inputFile;
     349  filebuf *inputBuffer;
    62350  TFile *outputFile = 0;
    63351  TStopwatch readStopWatch, procStopWatch;
     
    67355  Delphes *modularDelphes = 0;
    68356  DelphesFactory *factory = 0;
    69   TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
    70   DelphesHepMC3Reader *reader = 0;
     357  TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
     358  ReaderAscii *reader = 0;
     359  GenEvent event;
    71360  Int_t i, maxEvents, skipEvents;
    72361  Long64_t length, eventCounter;
     
    132421    partonOutputArray = modularDelphes->ExportArray("partons");
    133422
    134     reader = new DelphesHepMC3Reader;
     423    inputBuffer = inputFile.rdbuf();
     424    reader = new ReaderAscii(inputFile);
    135425
    136426    modularDelphes->InitTask();
     
    144434      {
    145435        cout << "** Reading standard input" << endl;
    146         inputFile = stdin;
     436        inputFile.ios::rdbuf(cin.rdbuf());
    147437        length = -1;
    148438      }
     
    150440      {
    151441        cout << "** Reading " << argv[i] << endl;
    152         inputFile = fopen(argv[i], "r");
    153 
    154         if(inputFile == NULL)
     442        inputFile.ios::rdbuf(inputBuffer);
     443        inputFile.open(argv[i]);
     444
     445        if(inputFile.fail())
    155446        {
    156447          message << "can't open " << argv[i];
     
    158449        }
    159450
    160         fseek(inputFile, 0L, SEEK_END);
    161         length = ftello(inputFile);
    162         fseek(inputFile, 0L, SEEK_SET);
     451        inputFile.seekg(0, ios::end);
     452        length = inputFile.tellg();
     453        inputFile.seekg(0, ios::beg);
    163454
    164455        if(length <= 0)
    165456        {
    166           fclose(inputFile);
     457          inputFile.close();
     458          inputFile.clear();
    167459          ++i;
    168460          continue;
     
    170462      }
    171463
    172       reader->SetInputFile(inputFile);
    173 
    174464      ExRootProgressBar progressBar(length);
    175465
    176       // Loop over all objects
    177       eventCounter = 0;
     466      reader->skip(skipEvents);
     467      progressBar.Update(inputFile.tellg(), skipEvents);
     468
     469      eventCounter = skipEvents;
     470      modularDelphes->Clear();
    178471      treeWriter->Clear();
    179       modularDelphes->Clear();
    180       reader->Clear();
     472      event.clear();
     473      gMotherMap.clear();
     474      gDaughterMap.clear();
    181475      readStopWatch.Start();
    182       while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && reader->ReadBlock(factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray) && !interrupted)
    183       {
    184         if(reader->EventReady())
    185         {
    186           ++eventCounter;
    187 
    188           readStopWatch.Stop();
    189 
    190           if(eventCounter > skipEvents)
    191           {
    192             procStopWatch.Start();
    193             modularDelphes->ProcessTask();
    194             procStopWatch.Stop();
    195 
    196             reader->AnalyzeEvent(branchEvent, eventCounter, &readStopWatch, &procStopWatch);
    197             reader->AnalyzeWeight(branchWeight);
    198 
    199             treeWriter->Fill();
    200 
    201             treeWriter->Clear();
    202           }
    203 
    204           modularDelphes->Clear();
    205           reader->Clear();
    206 
    207           readStopWatch.Start();
    208         }
    209         progressBar.Update(ftello(inputFile), eventCounter);
    210       }
    211 
    212       fseek(inputFile, 0L, SEEK_END);
    213       progressBar.Update(ftello(inputFile), eventCounter, kTRUE);
     476      reader->read_event(event);
     477      while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && !reader->failed() && !interrupted)
     478      {
     479        readStopWatch.Stop();
     480
     481        ++eventCounter;
     482
     483        procStopWatch.Start();
     484        AnalyzeEvent(event, branchEvent, factory,
     485          allParticleOutputArray, stableParticleOutputArray,
     486          partonOutputArray, &readStopWatch, &procStopWatch);
     487        modularDelphes->ProcessTask();
     488        AnalyzeWeight(event, branchWeight);
     489        procStopWatch.Stop();
     490
     491        treeWriter->Fill();
     492
     493        modularDelphes->Clear();
     494        treeWriter->Clear();
     495        event.clear();
     496        gMotherMap.clear();
     497        gDaughterMap.clear();
     498
     499        progressBar.Update(inputFile.tellg(), eventCounter);
     500
     501        readStopWatch.Start();
     502        reader->read_event(event);
     503      }
     504
     505      progressBar.Update(length, eventCounter, kTRUE);
    214506      progressBar.Finish();
    215507
    216       if(inputFile != stdin) fclose(inputFile);
     508      if(length > 0) inputFile.close();
    217509
    218510      ++i;
Note: See TracChangeset for help on using the changeset viewer.