source: HiSusy/trunk/Delphes-3.0.0/readers/DelphesHepMC.cpp @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 13.4 KB
Line 
1#include <stdexcept>
2#include <iostream>
3#include <sstream>
4
5#include <map>
6
7#include <stdlib.h>
8#include <signal.h>
9#include <stdio.h>
10
11#include "TROOT.h"
12#include "TApplication.h"
13
14#include "TFile.h"
15#include "TObjArray.h"
16#include "TStopwatch.h"
17#include "TDatabasePDG.h"
18#include "TParticlePDG.h"
19#include "TLorentzVector.h"
20
21#include "modules/Delphes.h"
22#include "classes/DelphesStream.h"
23#include "classes/DelphesClasses.h"
24#include "classes/DelphesFactory.h"
25
26#include "ExRootAnalysis/ExRootTreeWriter.h"
27#include "ExRootAnalysis/ExRootTreeBranch.h"
28#include "ExRootAnalysis/ExRootProgressBar.h"
29
30using namespace std;
31
32static const int kBufferSize  = 1024;
33
34class HepMCConverter
35{
36public:
37
38  HepMCConverter(DelphesFactory *factory, TObjArray *particleOutputArray,
39                 TObjArray *candidateOutputArray, TObjArray *partonOutputArray);
40  ~HepMCConverter();
41
42  void Clear();
43  Bool_t EventReady();
44
45  Bool_t ReadLine(FILE *inputFile);
46
47  void AnalyzeEvent(ExRootTreeBranch *branch,
48                    TStopwatch *readStopWatch, TStopwatch *procStopWatch);
49  void AnalyzeParticle();
50  void FinalizeParticles();
51 
52private:
53
54  Int_t fEventNumber, fMPI, fProcessID, fSignalCode, fVertexCounter;
55  Double_t fScale, fAlphaQCD, fAlphaQED;
56
57  Int_t fID1, fID2;
58  Double_t fX1, fX2, fScalePDF, fPDF1, fPDF2;
59
60  Int_t fOutVertexCode, fVertexID, fInCounter, fOutCounter;
61  Double_t fX, fY, fZ, fT;
62 
63  Int_t fParticleCode, fPID, fStatus, fInVertexCode;
64  Double_t fPx, fPy, fPz, fE, fMass, fTheta, fPhi;
65 
66  Int_t fParticleCounter;
67
68  char *fBuffer;
69
70  TDatabasePDG *fPDG;
71  DelphesFactory *fFactory;
72 
73  TIterator *fItParticleOutputArray;
74 
75  TObjArray *fParticleOutputArray;
76  TObjArray *fCandidateOutputArray;
77  TObjArray *fPartonOutputArray;
78
79  map< Int_t, pair < Int_t, Int_t > > fMotherMap;   
80  map< Int_t, pair < Int_t, Int_t > > fDaughterMap;   
81};
82
83//---------------------------------------------------------------------------
84
85void HepMCConverter::Clear()
86{
87  fVertexCounter = -1;
88  fInCounter = -1;
89  fOutCounter = -1;
90  fMotherMap.clear();
91  fParticleCounter = 0;
92}
93
94//---------------------------------------------------------------------------
95
96Bool_t HepMCConverter::EventReady()
97{
98  return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
99}
100
101//---------------------------------------------------------------------------
102
103Bool_t HepMCConverter::ReadLine(FILE *inputFile)
104{
105  map< Int_t, pair< Int_t, Int_t > >::iterator itMotherMap;
106  map< Int_t, pair< Int_t, Int_t > >::iterator itDaughterMap;
107  char key;
108  int rc;
109
110  if(!fgets(fBuffer, kBufferSize, inputFile)) return kFALSE;
111
112  DelphesStream bufferStream(fBuffer + 1);
113
114  key = fBuffer[0];
115
116  if(key == 'E')
117  {
118    Clear();
119
120    rc = bufferStream.ReadInt(fEventNumber)
121      && bufferStream.ReadInt(fMPI)
122      && bufferStream.ReadDbl(fScale)
123      && bufferStream.ReadDbl(fAlphaQCD)
124      && bufferStream.ReadDbl(fAlphaQED)
125      && bufferStream.ReadInt(fProcessID)
126      && bufferStream.ReadInt(fSignalCode)
127      && bufferStream.ReadInt(fVertexCounter);
128
129    if(!rc)
130    {
131      cerr << "** ERROR: " << "invalid event format" << endl;
132      return kFALSE;
133    }
134  }
135  else if(key == 'F')
136  {
137    rc = bufferStream.ReadInt(fID1)
138      && bufferStream.ReadInt(fID2)
139      && bufferStream.ReadDbl(fX1)
140      && bufferStream.ReadDbl(fX2)
141      && bufferStream.ReadDbl(fScalePDF)
142      && bufferStream.ReadDbl(fPDF1)
143      && bufferStream.ReadDbl(fPDF2);
144
145    if(!rc)
146    {
147      cerr << "** ERROR: " << "invalid PDF format" << endl;
148      return kFALSE;
149    }
150  }
151  else if(key == 'V' && fVertexCounter > 0)
152  {
153    rc = bufferStream.ReadInt(fOutVertexCode)
154      && bufferStream.ReadInt(fVertexID)
155      && bufferStream.ReadDbl(fX)
156      && bufferStream.ReadDbl(fY)
157      && bufferStream.ReadDbl(fZ)
158      && bufferStream.ReadDbl(fT)
159      && bufferStream.ReadInt(fInCounter)
160      && bufferStream.ReadInt(fOutCounter);
161
162    if(!rc)
163    {
164      cerr << "** ERROR: " << "invalid vertex format" << endl;
165      return kFALSE;
166    }
167    --fVertexCounter;
168  }
169  else if(key == 'P' && fOutCounter > 0)
170  {
171    rc = bufferStream.ReadInt(fParticleCode)
172      && bufferStream.ReadInt(fPID)
173      && bufferStream.ReadDbl(fPx)
174      && bufferStream.ReadDbl(fPy)
175      && bufferStream.ReadDbl(fPz)
176      && bufferStream.ReadDbl(fE)
177      && bufferStream.ReadDbl(fMass)
178      && bufferStream.ReadInt(fStatus)
179      && bufferStream.ReadDbl(fTheta)
180      && bufferStream.ReadDbl(fPhi)
181      && bufferStream.ReadInt(fInVertexCode);
182
183    if(!rc)
184    {
185      cerr << "** ERROR: " << "invalid particle format" << endl;
186      return kFALSE;
187    }
188
189    if(fInVertexCode < 0)
190    {
191      itMotherMap = fMotherMap.find(fInVertexCode);
192      if(itMotherMap == fMotherMap.end())
193      {
194        fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
195      }
196      else
197      {
198        itMotherMap->second.second = fParticleCounter;
199      }
200    }
201
202    if(fInCounter <= 0)
203    {
204      itDaughterMap = fDaughterMap.find(fOutVertexCode);
205      if(itDaughterMap == fDaughterMap.end())
206      {
207        fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
208      }
209      else
210      {
211        itDaughterMap->second.second = fParticleCounter;
212      }
213    }
214
215    AnalyzeParticle();
216
217    if(fInCounter > 0)
218    {
219      --fInCounter;
220    }
221    else
222    {
223      --fOutCounter;
224    }
225
226    ++fParticleCounter;
227  }
228 
229  return kTRUE;
230 
231}
232
233//---------------------------------------------------------------------------
234
235HepMCConverter::HepMCConverter(DelphesFactory *factory, TObjArray *particleOutputArray,
236                               TObjArray *candidateOutputArray, TObjArray *partonOutputArray) :
237  fBuffer(0), fPDG(0),
238  fFactory(factory),
239  fItParticleOutputArray(0),
240  fParticleOutputArray(particleOutputArray),
241  fCandidateOutputArray(candidateOutputArray),
242  fPartonOutputArray(partonOutputArray)
243{
244  fBuffer = new char[kBufferSize];
245   
246  fPDG = TDatabasePDG::Instance();
247  fItParticleOutputArray = fParticleOutputArray->MakeIterator();
248}
249
250//---------------------------------------------------------------------------
251
252HepMCConverter::~HepMCConverter()
253{
254  if(fItParticleOutputArray) delete fItParticleOutputArray;
255  if(fBuffer) delete[] fBuffer;
256}
257
258//---------------------------------------------------------------------------
259
260void HepMCConverter::AnalyzeEvent(ExRootTreeBranch *branch,
261                                  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
262{
263  HepMCEvent *element;
264 
265  element = static_cast<HepMCEvent *>(branch->NewEntry());
266  element->Number = fEventNumber;
267
268  element->ProcessID = fProcessID;
269  element->MPI = fMPI;
270  element->Scale = fScale;
271  element->AlphaQED = fAlphaQED;
272  element->AlphaQCD = fAlphaQCD;
273
274  element->ID1 = fID1;
275  element->ID2 = fID2;
276  element->X1 = fX1;
277  element->X2 = fX2;
278  element->ScalePDF = fScalePDF;
279  element->PDF1 = fPDF1;
280  element->PDF2 = fPDF2;
281 
282  element->ReadTime = readStopWatch->RealTime();
283  element->ProcTime = procStopWatch->RealTime();
284}
285
286//---------------------------------------------------------------------------
287
288void HepMCConverter::AnalyzeParticle()
289{
290  Candidate *candidate;
291  TParticlePDG *pdgParticle;
292
293  candidate = fFactory->NewCandidate();
294 
295  candidate->PID = fPID;
296
297  candidate->Status = fStatus;
298
299  pdgParticle = fPDG->GetParticle(fPID);
300  candidate->Charge = pdgParticle ? Int_t(pdgParticle->Charge()/3.0) : -999;
301  candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
302 
303  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
304
305  candidate->M2 = 1;
306  candidate->D2 = 1;
307  if(fInCounter > 0)
308  {
309    candidate->M1 = 1;
310    candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
311  }
312  else
313  {
314    candidate->M1 = fOutVertexCode;
315    candidate->Position.SetXYZT(fX, fY, fZ, fT);
316  }
317  if(fInVertexCode < 0)
318  {
319    candidate->D1 = fInVertexCode;
320  }
321  else
322  {
323    candidate->D1 = 1;
324  }
325
326  fParticleOutputArray->Add(candidate);
327 
328  if(!pdgParticle) return;
329
330  if(fStatus == 1)
331  {
332    fCandidateOutputArray->Add(candidate);
333  }
334  else if(fStatus == 2)
335  {
336    fPartonOutputArray->Add(candidate);
337  }
338}
339
340//---------------------------------------------------------------------------
341
342void HepMCConverter::FinalizeParticles()
343{
344  Candidate *candidate;
345  map< Int_t, pair< Int_t, Int_t > >::iterator itMotherMap;
346  map< Int_t, pair< Int_t, Int_t > >::iterator itDaughterMap;
347 
348  fItParticleOutputArray->Reset();
349  while((candidate = static_cast<Candidate *>(fItParticleOutputArray->Next())))
350  {
351    if(candidate->M1 > 0)
352    {
353      candidate->M1 = -1;
354      candidate->M2 = -1;
355    }
356    else
357    {
358      itMotherMap = fMotherMap.find(candidate->M1);
359      if(itMotherMap == fMotherMap.end())
360      {
361        candidate->M1 = -1;
362        candidate->M2 = -1;
363      }
364      else
365      {
366        candidate->M1 = itMotherMap->second.first;
367        candidate->M2 = itMotherMap->second.second;
368      }
369    }
370    if(candidate->D1 > 0)
371    {
372      candidate->D1 = -1;
373      candidate->D2 = -1;
374    }
375    else
376    {
377      itDaughterMap = fDaughterMap.find(candidate->D1);
378      if(itDaughterMap == fDaughterMap.end())
379      {
380        candidate->D1 = -1;
381        candidate->D2 = -1;
382      }
383      else
384      {
385        candidate->D1 = itDaughterMap->second.first;
386        candidate->D2 = itDaughterMap->second.second;
387      }
388    }
389  }
390}
391
392//---------------------------------------------------------------------------
393
394static bool interrupted = false;
395
396void SignalHandler(int sig)
397{
398        interrupted = true;
399}
400
401//---------------------------------------------------------------------------
402
403int main(int argc, char *argv[])
404{
405  char appName[] = "DelphesHepMC";
406  stringstream message;
407  FILE *inputFile = 0;
408  TFile *outputFile = 0;
409  TStopwatch readStopWatch, procStopWatch;
410  ExRootTreeWriter *treeWriter = 0;
411  ExRootTreeBranch *branchEvent = 0;
412  ExRootConfReader *confReader = 0;
413  Delphes *modularDelphes = 0;
414  DelphesFactory *factory = 0;
415  TObjArray *candidateOutputArray = 0, *particleOutputArray = 0, *partonOutputArray = 0;
416  HepMCConverter *converter = 0;
417  Int_t i;
418  Long64_t length, eventCounter;
419
420  if(argc < 3)
421  {
422    cout << " Usage: " << appName << " config_file" << " output_file" << " [input_file(s)]" << endl;
423    cout << " config_file - configuration file in Tcl format," << endl;
424    cout << " output_file - output file in ROOT format," << endl;
425    cout << " input_file(s) - input file(s) in HepMC format," << endl;
426    cout << " with no input_file, or when input_file is -, read standard input." << endl;
427    return 1;
428  }
429
430  signal(SIGINT, SignalHandler);
431
432  gROOT->SetBatch();
433
434  int appargc = 1;
435  char *appargv[] = {appName};
436  TApplication app(appName, &appargc, appargv);
437       
438  try
439  {
440    outputFile = TFile::Open(argv[2], "RECREATE");
441
442    if(outputFile == NULL)
443    {
444      message << "can't open " << argv[2];
445      throw runtime_error(message.str());
446    }
447
448    treeWriter = new ExRootTreeWriter(outputFile, "Delphes");
449
450    branchEvent = treeWriter->NewBranch("Event", HepMCEvent::Class());
451
452    confReader = new ExRootConfReader;
453    confReader->ReadFile(argv[1]);
454   
455    modularDelphes = new Delphes("Delphes");
456    modularDelphes->SetConfReader(confReader);
457    modularDelphes->SetTreeWriter(treeWriter);
458   
459    factory = modularDelphes->GetFactory();
460    particleOutputArray = modularDelphes->ExportArray("particles");
461    candidateOutputArray = modularDelphes->ExportArray("candidates");
462    partonOutputArray = modularDelphes->ExportArray("partons");
463
464    converter = new HepMCConverter(factory, particleOutputArray, candidateOutputArray, partonOutputArray);
465
466    modularDelphes->InitTask();
467 
468    i = 3;
469    do
470    {
471      if(interrupted) break;
472
473      if(i == argc || strcmp(argv[i], "-") == 0)
474      {
475        cout << "** Reading standard input" << endl;
476        inputFile = stdin;
477        length = -1;
478      }
479      else
480      {
481        cout << "** Reading " << argv[i] << endl;
482        inputFile = fopen(argv[i], "r");
483
484        if(inputFile == NULL)
485        {
486          message << "can't open " << argv[i];
487          throw runtime_error(message.str());
488        }
489
490        fseek(inputFile, 0L, SEEK_END);
491        length = ftell(inputFile);
492        fseek(inputFile, 0L, SEEK_SET);
493
494        if(length <= 0) continue;
495      }
496
497      ExRootProgressBar progressBar(length);
498 
499      // Loop over all objects
500      eventCounter = 0;
501      treeWriter->Clear();
502      modularDelphes->Clear();
503      converter->Clear();
504      readStopWatch.Start();
505      while(converter->ReadLine(inputFile) && !interrupted)
506      {
507        if(converter->EventReady())
508        {
509          ++eventCounter;
510
511          converter->FinalizeParticles();
512
513          readStopWatch.Stop();
514          procStopWatch.Start();
515          modularDelphes->ProcessTask();
516          procStopWatch.Stop();
517
518          converter->AnalyzeEvent(branchEvent, &readStopWatch, &procStopWatch);
519   
520          treeWriter->Fill();
521
522          treeWriter->Clear();
523          modularDelphes->Clear();
524          converter->Clear();
525
526          readStopWatch.Start();
527        }
528        progressBar.Update(ftell(inputFile), eventCounter);
529      }
530
531      fseek(inputFile, 0L, SEEK_END);
532      progressBar.Update(ftell(inputFile), eventCounter, kTRUE);
533      progressBar.Finish();
534
535      if(inputFile != stdin) fclose(inputFile);
536
537      ++i;
538    }
539    while(i < argc);
540
541    modularDelphes->FinishTask();
542    treeWriter->Write();
543
544    cout << "** Exiting..." << endl;
545
546    delete converter;
547    delete modularDelphes;
548    delete confReader;
549    delete treeWriter;
550    delete outputFile;
551   
552    return 0;
553  }
554  catch(runtime_error &e)
555  {
556    if(treeWriter) delete treeWriter;
557    if(outputFile) delete outputFile;
558    cerr << "** ERROR: " << e.what() << endl;
559    return 1;
560  }
561}
562
Note: See TracBrowser for help on using the repository browser.