Fork me on GitHub

source: svn/trunk/src/LHEFConverter.cc@ 219

Last change on this file since 219 was 219, checked in by Xavier Rouby, 16 years ago

JetUtils.cc

File size: 3.9 KB
Line 
1#include <iostream>
2#include <fstream>
3
4#include "TLorentzVector.h"
5
6#include "ExRootTreeWriter.h"
7#include "BlockClasses.h"
8#include "LHEFConverter.h"
9
10using namespace std;
11
12
13//------------------------------------------------------------------------------
14
15void LHEFConverter::AnalyseEvent(LHEF::Reader *reader, ExRootTreeBranch *branch, const Long64_t eventNumber)
16{
17 const LHEF::HEPEUP &hepeup = reader->hepeup;
18
19 TRootLHEFEvent *element;
20
21 element = (TRootLHEFEvent*) branch->NewEntry();
22
23 element->Number = eventNumber;
24 element->Nparticles = hepeup.NUP;
25 element->ProcessID = hepeup.IDPRUP;
26 element->Weight = hepeup.XWGTUP;
27 element->ScalePDF = hepeup.SCALUP;
28 element->CouplingQED = hepeup.AQEDUP;
29 element->CouplingQCD = hepeup.AQCDUP;
30}
31
32//---------------------------------------------------------------------------
33
34void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
35{
36 const LHEF::HEPEUP &hepeup = reader->hepeup;
37 Double_t signPz;
38 TLorentzVector momentum;
39 TRootLHEFParticle *element;
40
41 for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
42 {
43 element = (TRootLHEFParticle*) branch->NewEntry();
44
45 element->PID = hepeup.IDUP[particle];
46 element->Status = hepeup.ISTUP[particle];
47 element->Mother1 = hepeup.MOTHUP[particle].first;
48 element->Mother2 = hepeup.MOTHUP[particle].second;
49 element->ColorLine1 = hepeup.ICOLUP[particle].first;
50 element->ColorLine2 = hepeup.ICOLUP[particle].second;
51 element->Px = hepeup.PUP[particle][0];
52 element->Py = hepeup.PUP[particle][1];
53 element->Pz = hepeup.PUP[particle][2];
54 element->E = hepeup.PUP[particle][3];
55 element->M = hepeup.PUP[particle][4];
56
57 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
58 element->PT = momentum.Perp();
59 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
60 element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta();
61 element->Phi = momentum.Phi();
62 element->Rapidity = element->PT == 0.0 ? signPz*999.9 : momentum.Rapidity();
63 element->LifeTime = hepeup.VTIMUP[particle];
64 element->Spin = hepeup.SPINUP[particle];
65 }
66}
67
68LHEFConverter::~LHEFConverter()
69{
70}
71
72//------------------------------------------------------------------------------
73
74LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName) {
75
76 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
77
78 // generated event from LHEF
79 ExRootTreeBranch *branchEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
80
81 // generated partons from LHEF
82 ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootLHEFParticle::Class());
83
84 // Open a stream connected to an event file:
85 ifstream infile(inputFileList.c_str());
86 string filename;
87 if(!infile.is_open()) {
88 cerr << "** ERROR: Can't open '" << inputFileList << "' for input" << endl;
89 exit(1);
90 }
91
92 while(1) { // parsing the list of files
93
94 infile >> filename;
95 if(!infile.good()) break;
96 ifstream checking_the_file(filename.c_str());
97 if(!checking_the_file.good()) { cout << filename << ": file not found\n"; continue;}
98 else checking_the_file.close();
99
100 // Create the Reader object:
101 LHEF::Reader *inputReader = new LHEF::Reader(filename);
102
103 cout << "** Calculating number of events to process. Please wait..." << endl;
104 Long64_t allEntries = inputReader->getNumberOfEvents();
105 cout << "** Input file contains " << allEntries << " events" << endl;
106
107 if(allEntries > 0) {
108 // Loop over all events
109 Long64_t entry = 0;
110 while(inputReader->readEvent()) {
111 treeWriter->Clear();
112 AnalyseEvent(inputReader, branchEvent, entry + 1);
113 AnalyseParticles(inputReader, branchParticle);
114 treeWriter->Fill();
115 ++entry;
116 }
117 }
118 }
119 treeWriter->Write();
120
121 cout << "** Exiting..." << endl;
122
123 delete treeWriter;
124 //delete inputReader;
125}
126
Note: See TracBrowser for help on using the repository browser.