Fork me on GitHub

Changeset fc6bce2 in git for readers/DelphesHepMC3.cpp


Ignore:
Timestamp:
Jun 1, 2021, 8:13:13 PM (3 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
master
Children:
bfc6722
Parents:
4006893
Message:

remove HepMC3 library

File:
1 edited

Legend:

Unmodified
Added
Removed

  • readers/DelphesHepMC3.cpp

    r4006893 rfc6bce2  
    1818
    1919#include <iostream>
    20 #include <fstream>
    21 #include <memory>
    2220#include <sstream>
    2321#include <stdexcept>
    2422
    25 #include <map>
    26 
    2723#include <signal.h>
    28 #include <stdio.h>
    29 #include <stdlib.h>
    3024
    3125#include "TApplication.h"
    … …  
    4135#include "classes/DelphesClasses.h"
    4236#include "classes/DelphesFactory.h"
     37#include "classes/DelphesHepMC3Reader.h"
    4338#include "modules/Delphes.h"
    4439
    … …  
    4742#include "ExRootAnalysis/ExRootTreeWriter.h"
    4843
    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 
    5744using namespace std;
    58 using namespace HepMC3;
    59 
    60 map<Int_t, pair<Int_t, Int_t> > gMotherMap;
    61 map<Int_t, pair<Int_t, Int_t> > gDaughterMap;
    62 
    63 //---------------------------------------------------------------------------
    64 
    65 void 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 
    184 void 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 
    321 void 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 }
    33245
    33346//---------------------------------------------------------------------------
    … …  
    34659  char appName[] = "DelphesHepMC3";
    34760  stringstream message;
    348   ifstream inputFile;
    349   filebuf *inputBuffer;
     61  FILE *inputFile = 0;
    35062  TFile *outputFile = 0;
    35163  TStopwatch readStopWatch, procStopWatch;
    … …  
    35567  Delphes *modularDelphes = 0;
    35668  DelphesFactory *factory = 0;
    357   TObjArray *allParticleOutputArray = 0, *stableParticleOutputArray = 0, *partonOutputArray = 0;
    358   ReaderAscii *reader = 0;
    359   GenEvent event;
     69  TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
     70  DelphesHepMC3Reader *reader = 0;
    36071  Int_t i, maxEvents, skipEvents;
    36172  Long64_t length, eventCounter;
    … …  
    421132    partonOutputArray = modularDelphes->ExportArray("partons");
    422133
    423     inputBuffer = inputFile.rdbuf();
    424     reader = new ReaderAscii(inputFile);
     134    reader = new DelphesHepMC3Reader;
    425135
    426136    modularDelphes->InitTask();
    … …  
    434144      {
    435145        cout << "** Reading standard input" << endl;
    436         inputFile.ios::rdbuf(cin.rdbuf());
     146        inputFile = stdin;
    437147        length = -1;
    438148      }
    … …  
    440150      {
    441151        cout << "** Reading " << argv[i] << endl;
    442         inputFile.ios::rdbuf(inputBuffer);
    443         inputFile.open(argv[i]);
    444 
    445         if(inputFile.fail())
     152        inputFile = fopen(argv[i], "r");
     153
     154        if(inputFile == NULL)
    446155        {
    447156          message << "can't open " << argv[i];
    … …  
    449158        }
    450159
    451         inputFile.seekg(0, ios::end);
    452         length = inputFile.tellg();
    453         inputFile.seekg(0, ios::beg);
     160        fseek(inputFile, 0L, SEEK_END);
     161        length = ftello(inputFile);
     162        fseek(inputFile, 0L, SEEK_SET);
    454163
    455164        if(length <= 0)
    456165        {
    457           inputFile.close();
    458           inputFile.clear();
     166          fclose(inputFile);
    459167          ++i;
    460168          continue;
    … …  
    462170      }
    463171
     172      reader->SetInputFile(inputFile);
     173
    464174      ExRootProgressBar progressBar(length);
    465175
    466       reader->skip(skipEvents);
    467       progressBar.Update(inputFile.tellg(), skipEvents);
    468 
    469       eventCounter = skipEvents;
     176      // Loop over all objects
     177      eventCounter = 0;
     178      treeWriter->Clear();
    470179      modularDelphes->Clear();
    471       treeWriter->Clear();
    472       event.clear();
    473       gMotherMap.clear();
    474       gDaughterMap.clear();
     180      reader->Clear();
    475181      readStopWatch.Start();
    476       reader->read_event(event);
    477       while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && !reader->failed() && !interrupted)
     182      while((maxEvents <= 0 || eventCounter - skipEvents < maxEvents) && reader->ReadBlock(factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray) && !interrupted)
    478183      {
    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);
     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);
    503210      }
    504211
    505       progressBar.Update(length, eventCounter, kTRUE);
     212      fseek(inputFile, 0L, SEEK_END);
     213      progressBar.Update(ftello(inputFile), eventCounter, kTRUE);
    506214      progressBar.Finish();
    507215
    508       if(length > 0) inputFile.close();
     216      if(inputFile != stdin) fclose(inputFile);
    509217
    510218      ++i;
Note: See TracChangeset for help on using the changeset viewer.