- Timestamp:
- Jun 1, 2021, 8:13:13 PM (3 years ago)
- Branches:
- master
- Children:
- bfc6722
- Parents:
- 4006893
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
readers/DelphesHepMC3.cpp
r4006893 rfc6bce2 18 18 19 19 #include <iostream> 20 #include <fstream>21 #include <memory>22 20 #include <sstream> 23 21 #include <stdexcept> 24 22 25 #include <map>26 27 23 #include <signal.h> 28 #include <stdio.h>29 #include <stdlib.h>30 24 31 25 #include "TApplication.h" … … 41 35 #include "classes/DelphesClasses.h" 42 36 #include "classes/DelphesFactory.h" 37 #include "classes/DelphesHepMC3Reader.h" 43 38 #include "modules/Delphes.h" 44 39 … … 47 42 #include "ExRootAnalysis/ExRootTreeWriter.h" 48 43 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 57 44 using 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 else130 {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 else144 {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 else159 {160 itMotherMap->second.second = counter;161 }162 }163 else164 {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 else222 {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 else238 {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 else284 {285 itMotherMap = gMotherMap.find(candidate->M1);286 if(itMotherMap == gMotherMap.end())287 {288 candidate->M1 = -1;289 candidate->M2 = -1;290 }291 else292 {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 else303 {304 itDaughterMap = gDaughterMap.find(candidate->D1);305 if(itDaughterMap == gDaughterMap.end())306 {307 candidate->D1 = -1;308 candidate->D2 = -1;309 }310 else311 {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 }332 45 333 46 //--------------------------------------------------------------------------- … … 346 59 char appName[] = "DelphesHepMC3"; 347 60 stringstream message; 348 ifstream inputFile; 349 filebuf *inputBuffer; 61 FILE *inputFile = 0; 350 62 TFile *outputFile = 0; 351 63 TStopwatch readStopWatch, procStopWatch; … … 355 67 Delphes *modularDelphes = 0; 356 68 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; 360 71 Int_t i, maxEvents, skipEvents; 361 72 Long64_t length, eventCounter; … … 421 132 partonOutputArray = modularDelphes->ExportArray("partons"); 422 133 423 inputBuffer = inputFile.rdbuf(); 424 reader = new ReaderAscii(inputFile); 134 reader = new DelphesHepMC3Reader; 425 135 426 136 modularDelphes->InitTask(); … … 434 144 { 435 145 cout << "** Reading standard input" << endl; 436 inputFile .ios::rdbuf(cin.rdbuf());146 inputFile = stdin; 437 147 length = -1; 438 148 } … … 440 150 { 441 151 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) 446 155 { 447 156 message << "can't open " << argv[i]; … … 449 158 } 450 159 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); 454 163 455 164 if(length <= 0) 456 165 { 457 inputFile.close(); 458 inputFile.clear(); 166 fclose(inputFile); 459 167 ++i; 460 168 continue; … … 462 170 } 463 171 172 reader->SetInputFile(inputFile); 173 464 174 ExRootProgressBar progressBar(length); 465 175 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(); 470 179 modularDelphes->Clear(); 471 treeWriter->Clear(); 472 event.clear(); 473 gMotherMap.clear(); 474 gDaughterMap.clear(); 180 reader->Clear(); 475 181 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) 478 183 { 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); 503 210 } 504 211 505 progressBar.Update(length, eventCounter, kTRUE); 212 fseek(inputFile, 0L, SEEK_END); 213 progressBar.Update(ftello(inputFile), eventCounter, kTRUE); 506 214 progressBar.Finish(); 507 215 508 if( length > 0) inputFile.close();216 if(inputFile != stdin) fclose(inputFile); 509 217 510 218 ++i;
Note:
See TracChangeset
for help on using the changeset viewer.