Fork me on GitHub

Changeset 341014c in git for readers/DelphesProIO.cpp


Ignore:
Timestamp:
Feb 12, 2019, 9:29:17 PM (6 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, llp, master
Children:
6455202
Parents:
45e58be
Message:

apply .clang-format to all .h, .cc and .cpp files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • readers/DelphesProIO.cpp

    r45e58be r341014c  
    1818 */
    1919
     20#include <iostream>
     21#include <memory>
     22#include <sstream>
    2023#include <stdexcept>
    21 #include <iostream>
    22 #include <sstream>
    23 #include <memory>
    2424
    2525#include <map>
    2626
    27 #include <stdlib.h>
    2827#include <signal.h>
    2928#include <stdio.h>
    30 
     29#include <stdlib.h>
     30
     31#include "TApplication.h"
    3132#include "TROOT.h"
    32 #include "TApplication.h"
    33 
     33
     34#include "TDatabasePDG.h"
    3435#include "TFile.h"
     36#include "TLorentzVector.h"
    3537#include "TObjArray.h"
     38#include "TParticlePDG.h"
    3639#include "TStopwatch.h"
    37 #include "TDatabasePDG.h"
    38 #include "TParticlePDG.h"
    39 #include "TLorentzVector.h"
    40 
    41 #include "modules/Delphes.h"
    42 #include "classes/DelphesStream.h"
     40
    4341#include "classes/DelphesClasses.h"
    4442#include "classes/DelphesFactory.h"
    45 
     43#include "classes/DelphesStream.h"
     44#include "modules/Delphes.h"
     45
     46#include "ExRootAnalysis/ExRootProgressBar.h"
     47#include "ExRootAnalysis/ExRootTreeBranch.h"
    4648#include "ExRootAnalysis/ExRootTreeWriter.h"
    47 #include "ExRootAnalysis/ExRootTreeBranch.h"
    48 #include "ExRootAnalysis/ExRootProgressBar.h"
    49 
    50 #include <proio/reader.h>
     49
    5150#include <proio/event.h>
    5251#include <proio/model/mc.pb.h>
    53 namespace model=proio::model::mc;
    54 
     52#include <proio/reader.h>
     53namespace model = proio::model::mc;
    5554
    5655using namespace std;
    … …  
    5958// This method dynamically checks the message type (varint or not) depending on
    6059// non-zero value of units momentumUnit and positionUnit.
    61 void ConvertInput(proio::Event *event, double momentumUnit, double positionUnit, 
     60void ConvertInput(proio::Event *event, double momentumUnit, double positionUnit,
    6261  ExRootTreeBranch *branch, DelphesFactory *factory,
    6362  TObjArray *allParticleOutputArray, TObjArray *stableParticleOutputArray,
    … …  
    8079  element = static_cast<HepMCEvent *>(branch->NewEntry());
    8180
    82  int nID=0;
    83  double weight=0;
    84  int process_id=0;
    85  auto mcentries = event->TaggedEntries("MCParameters");
    86  for (uint64_t mcentryID : mcentries) {
    87    auto mcpar = dynamic_cast<proio::model::mc::MCParameters *>(event->GetEntry(mcentryID));
    88    nID=mcpar->number();
    89    weight=mcpar->weight();
    90    process_id=mcpar->processid();
    91    break; // consider only most generic from 1st entry
    92  };
     81  int nID = 0;
     82  double weight = 0;
     83  int process_id = 0;
     84  auto mcentries = event->TaggedEntries("MCParameters");
     85  for(uint64_t mcentryID : mcentries)
     86  {
     87    auto mcpar = dynamic_cast<proio::model::mc::MCParameters *>(event->GetEntry(mcentryID));
     88    nID = mcpar->number();
     89    weight = mcpar->weight();
     90    process_id = mcpar->processid();
     91    break; // consider only most generic from 1st entry
     92  };
    9393
    9494  element->Number = nID;
    95   element->ProcessID =process_id;
     95  element->ProcessID = process_id;
    9696  element->Weight = weight;
    9797
    98 /*
     98  /*
    9999  // Pythia8 specific
    100100  element->MPI = mutableEvent->mpi();
    … …  
    114114  element->ProcTime = procStopWatch->RealTime();
    115115
    116 
    117 
    118  if ( momentumUnit >0 && positionUnit>0) {
    119 
    120 
    121   auto entries = event->TaggedEntries("VarintPackedParticles");
    122 
    123   for (uint64_t entryID : entries) {
    124 
    125     auto mutableParticles = dynamic_cast<model::VarintPackedParticles *>(event->GetEntry(entryID));
    126 
    127     for(int i = 0; i < mutableParticles->pdg_size(); ++i)
    128    {
    129     pid = mutableParticles->pdg(i);
    130     status = mutableParticles->status(i);
    131     px = mutableParticles->px(i)/momentumUnit;
    132     py = mutableParticles->py(i)/momentumUnit;
    133     pz = mutableParticles->pz(i)/momentumUnit;
    134     mass = mutableParticles->mass(i)/momentumUnit;
    135     x = mutableParticles->x(i)/positionUnit;
    136     y = mutableParticles->y(i)/positionUnit;
    137     z = mutableParticles->z(i)/positionUnit;
    138     t = mutableParticles->t(i)/positionUnit;
    139    
    140     candidate = factory->NewCandidate();
    141     candidate->PID = pid;
    142     pdgCode = TMath::Abs(candidate->PID);
    143     candidate->Status = status;
    144     candidate->M1 = mutableParticles->parent1(i);
    145     candidate->M2 = mutableParticles->parent2(i);
    146     candidate->D1 = mutableParticles->child1(i);
    147     candidate->D2 = mutableParticles->child2(i);
    148     pdgParticle = pdg->GetParticle(pid);
    149     candidate->Charge = mutableParticles->charge(i)/3.0;
    150     //candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
    151     candidate->Mass = mass;
    152     candidate->Momentum.SetXYZM(px, py, pz, mass);
    153     candidate->Position.SetXYZT(x, y, z, t);
    154     allParticleOutputArray->Add(candidate);
    155     if(!pdgParticle) continue;
    156 
    157     if(status == 1)
     116  if(momentumUnit > 0 && positionUnit > 0)
     117  {
     118
     119    auto entries = event->TaggedEntries("VarintPackedParticles");
     120
     121    for(uint64_t entryID : entries)
    158122    {
    159       stableParticleOutputArray->Add(candidate);
    160     }
    161     else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
    162     {
    163       partonOutputArray->Add(candidate);
     123
     124      auto mutableParticles = dynamic_cast<model::VarintPackedParticles *>(event->GetEntry(entryID));
     125
     126      for(int i = 0; i < mutableParticles->pdg_size(); ++i)
     127      {
     128        pid = mutableParticles->pdg(i);
     129        status = mutableParticles->status(i);
     130        px = mutableParticles->px(i) / momentumUnit;
     131        py = mutableParticles->py(i) / momentumUnit;
     132        pz = mutableParticles->pz(i) / momentumUnit;
     133        mass = mutableParticles->mass(i) / momentumUnit;
     134        x = mutableParticles->x(i) / positionUnit;
     135        y = mutableParticles->y(i) / positionUnit;
     136        z = mutableParticles->z(i) / positionUnit;
     137        t = mutableParticles->t(i) / positionUnit;
     138
     139        candidate = factory->NewCandidate();
     140        candidate->PID = pid;
     141        pdgCode = TMath::Abs(candidate->PID);
     142        candidate->Status = status;
     143        candidate->M1 = mutableParticles->parent1(i);
     144        candidate->M2 = mutableParticles->parent2(i);
     145        candidate->D1 = mutableParticles->child1(i);
     146        candidate->D2 = mutableParticles->child2(i);
     147        pdgParticle = pdg->GetParticle(pid);
     148        candidate->Charge = mutableParticles->charge(i) / 3.0;
     149        //candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
     150        candidate->Mass = mass;
     151        candidate->Momentum.SetXYZM(px, py, pz, mass);
     152        candidate->Position.SetXYZT(x, y, z, t);
     153        allParticleOutputArray->Add(candidate);
     154        if(!pdgParticle) continue;
     155
     156        if(status == 1)
     157        {
     158          stableParticleOutputArray->Add(candidate);
     159        }
     160        else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
     161        {
     162          partonOutputArray->Add(candidate);
     163        }
     164      }
    164165    }
    165166  }
    166 
    167   }
    168 
    169  } else {
    170 
    171 
    172  auto entries = event->TaggedEntries("Particle");
    173 
    174  for (uint64_t entryID : entries) {
    175     auto part = dynamic_cast<model::Particle *>(event->GetEntry(entryID));
    176     pid = part->pdg();
    177     status = part->status();
    178     model::XYZF pvector=part->p();
    179     px=pvector.x();
    180     py=pvector.y();
    181     pz=pvector.z();
    182     mass = part->mass();
    183     auto v =part->vertex();
    184     x=v.x();
    185     y=v.y();
    186     z=v.z();
    187     t=v.t();
    188 
    189     candidate = factory->NewCandidate();
    190     candidate->PID = pid;
    191     pdgCode = TMath::Abs(candidate->PID);
    192     candidate->Status = status;
    193 
    194     int M1=0;
    195     int M2=0;
    196     for (int k1=0; k1<part->parent_size(); k1++){
    197        if (k1==0) {
    198                    auto mother = dynamic_cast<model::Particle *>(event->GetEntry(part->parent(0)));
    199                    M1=mother->barcode();
    200        }
    201        if (k1==1) {
    202                    auto mother = dynamic_cast<model::Particle *>(event->GetEntry(part->parent(1)));
    203                    M2=mother->barcode();
    204        }
     167  else
     168  {
     169
     170    auto entries = event->TaggedEntries("Particle");
     171
     172    for(uint64_t entryID : entries)
     173    {
     174      auto part = dynamic_cast<model::Particle *>(event->GetEntry(entryID));
     175      pid = part->pdg();
     176      status = part->status();
     177      model::XYZF pvector = part->p();
     178      px = pvector.x();
     179      py = pvector.y();
     180      pz = pvector.z();
     181      mass = part->mass();
     182      auto v = part->vertex();
     183      x = v.x();
     184      y = v.y();
     185      z = v.z();
     186      t = v.t();
     187
     188      candidate = factory->NewCandidate();
     189      candidate->PID = pid;
     190      pdgCode = TMath::Abs(candidate->PID);
     191      candidate->Status = status;
     192
     193      int M1 = 0;
     194      int M2 = 0;
     195      for(int k1 = 0; k1 < part->parent_size(); k1++)
     196      {
     197        if(k1 == 0)
     198        {
     199          auto mother = dynamic_cast<model::Particle *>(event->GetEntry(part->parent(0)));
     200          M1 = mother->barcode();
     201        }
     202        if(k1 == 1)
     203        {
     204          auto mother = dynamic_cast<model::Particle *>(event->GetEntry(part->parent(1)));
     205          M2 = mother->barcode();
     206        }
     207      }
     208
     209      int D1 = 0;
     210      int D2 = 0;
     211      for(int k1 = 0; k1 < part->child_size(); k1++)
     212      {
     213        if(k1 == 0)
     214        {
     215          auto child = dynamic_cast<model::Particle *>(event->GetEntry(part->child(0)));
     216          D1 = child->barcode();
     217        }
     218
     219        if(k1 == 1)
     220        {
     221          auto child = dynamic_cast<model::Particle *>(event->GetEntry(part->child(1)));
     222          D2 = child->barcode();
     223        };
     224      }
     225
     226      candidate->M1 = M1;
     227      candidate->M2 = M2;
     228      candidate->D1 = D1;
     229      candidate->D2 = D2;
     230
     231      pdgParticle = pdg->GetParticle(pid);
     232      candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge() / 3.0) : -999;
     233      candidate->Mass = mass;
     234
     235      candidate->Momentum.SetXYZM(px, py, pz, mass);
     236
     237      candidate->Position.SetXYZT(x, y, z, t);
     238
     239      allParticleOutputArray->Add(candidate);
     240
     241      if(!pdgParticle) continue;
     242
     243      if(status == 1)
     244      {
     245        stableParticleOutputArray->Add(candidate);
     246      }
     247      else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
     248      {
     249        partonOutputArray->Add(candidate);
     250      }
    205251    }
    206252
    207 
    208     int D1=0;
    209     int D2=0;
    210     for (int k1=0; k1<part->child_size(); k1++){
    211        if (k1==0) {
    212                    auto child = dynamic_cast<model::Particle *>(event->GetEntry(part->child(0)));
    213                    D1=child->barcode();
    214        }
    215 
    216       if (k1==1) {
    217                    auto child = dynamic_cast<model::Particle *>(event->GetEntry(part->child(1)));
    218                    D2=child->barcode();
    219          };
    220       }
    221 
    222 
    223     candidate->M1 = M1;
    224     candidate->M2 = M2;
    225     candidate->D1 = D1;
    226     candidate->D2 = D2;
    227 
    228     pdgParticle = pdg->GetParticle(pid);
    229     candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
    230     candidate->Mass = mass;
    231 
    232     candidate->Momentum.SetXYZM(px, py, pz, mass);
    233 
    234     candidate->Position.SetXYZT(x, y, z, t);
    235 
    236     allParticleOutputArray->Add(candidate);
    237 
    238     if(!pdgParticle) continue;
    239 
    240     if(status == 1)
    241     {
    242       stableParticleOutputArray->Add(candidate);
    243     }
    244     else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
    245     {
    246       partonOutputArray->Add(candidate);
    247     }
    248 
    249   }
    250 
    251 
    252   } // end particle type
    253 
     253  } // end particle type
    254254}
    255255
    … …  
    283283  if(argc < 4)
    284284  {
    285     cout << " Usage: " << appName << " config_file" << " output_file" << " input_file(s)" << endl;
     285    cout << " Usage: " << appName << " config_file"
     286         << " output_file"
     287         << " input_file(s)" << endl;
    286288    cout << " config_file - configuration file in Tcl format," << endl;
    287289    cout << " output_file - output file in ROOT format," << endl;
    … …  
    330332      cout << "** INFO: Reading " << argv[i] << endl;
    331333
    332       inputFile=new proio::Reader( argv[i] );
     334      inputFile = new proio::Reader(argv[i]);
    333335
    334336      if(inputFile == NULL)
    … …  
    338340      }
    339341
    340 
    341 /*
     342      /*
    342343// this is slow method, but general
    343344      inputFile->SeekToStart();
    … …  
    350351*/
    351352
    352 
    353     auto event = new proio::Event();
    354 
    355 
    356    double varint_energy=0;
    357    double varint_length=0;
    358 
    359 
    360     auto max_n_events = std::numeric_limits<uint64_t>::max();
    361     auto nn = inputFile->Skip(max_n_events);
    362     cout << "** INFO: " << nn-1 << " events found in ProIO file" << endl;
    363     inputFile->SeekToStart();
    364     numberOfEvents = nn-1; // last event has only metadata (no particle record)
    365  
    366     if(numberOfEvents <= 0) continue;
    367 
    368     ExRootProgressBar progressBar(numberOfEvents - 1);
    369 
     353      auto event = new proio::Event();
     354
     355      double varint_energy = 0;
     356      double varint_length = 0;
     357
     358      auto max_n_events = std::numeric_limits<uint64_t>::max();
     359      auto nn = inputFile->Skip(max_n_events);
     360      cout << "** INFO: " << nn - 1 << " events found in ProIO file" << endl;
     361      inputFile->SeekToStart();
     362      numberOfEvents = nn - 1; // last event has only metadata (no particle record)
     363
     364      if(numberOfEvents <= 0) continue;
     365
     366      ExRootProgressBar progressBar(numberOfEvents - 1);
    370367
    371368      // Loop over all objects
    … …  
    376373      for(eventCounter = 0; eventCounter < numberOfEvents && !interrupted; ++eventCounter)
    377374      {
    378         inputFile->Next(event); 
     375        inputFile->Next(event);
    379376        if(event == 0) continue;
    380377
    381        // get metadata
    382        if (eventCounter == 0) {
    383        auto metadata = event->Metadata();
    384        std::cout << "** INFO: ProIO file metadata:" << std::endl;
    385        for (auto element : metadata) {
    386         string key=(string)element.first;
    387         string value=(string)(*element.second);
    388         std::cout << "** INFO:   " << key << " = " << value << std::endl;
    389         if (key=="info:varint_energy") varint_energy=std::stod(value);
    390         if (key=="info:varint_length") varint_length=std::stod(value);
    391         }
    392        }
    393 
     378        // get metadata
     379        if(eventCounter == 0)
     380        {
     381          auto metadata = event->Metadata();
     382          std::cout << "** INFO: ProIO file metadata:" << std::endl;
     383          for(auto element : metadata)
     384          {
     385            string key = (string)element.first;
     386            string value = (string)(*element.second);
     387            std::cout << "** INFO:   " << key << " = " << value << std::endl;
     388            if(key == "info:varint_energy") varint_energy = std::stod(value);
     389            if(key == "info:varint_length") varint_length = std::stod(value);
     390          }
     391        }
    394392
    395393        readStopWatch.Stop();
    … …  
    397395        procStopWatch.Start();
    398396
    399         ConvertInput(event, varint_energy, varint_length, 
     397        ConvertInput(event, varint_energy, varint_length,
    400398          branchEvent, factory,
    401399          allParticleOutputArray, stableParticleOutputArray,
    402400          partonOutputArray, &readStopWatch, &procStopWatch);
    403  
     401
    404402        modularDelphes->ProcessTask();
    405403        procStopWatch.Stop();
    … …  
    414412      }
    415413
    416 
    417414      progressBar.Update(eventCounter, eventCounter, kTRUE);
    418415      progressBar.Finish();
    … …  
    420417      delete inputFile;
    421418    }
    422 
    423419
    424420    modularDelphes->FinishTask();
Note: See TracChangeset for help on using the changeset viewer.