1 | /***********************************************************************
|
---|
2 | ** **
|
---|
3 | ** /----------------------------------------------\ **
|
---|
4 | ** | Delphes, a framework for the fast simulation | **
|
---|
5 | ** | of a generic collider experiment | **
|
---|
6 | ** \------------- arXiv:0903.2225v1 ------------/ **
|
---|
7 | ** **
|
---|
8 | ** **
|
---|
9 | ** This package uses: **
|
---|
10 | ** ------------------ **
|
---|
11 | ** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
|
---|
12 | ** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
|
---|
13 | ** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
|
---|
14 | ** FROG: [hep-ex/0901.2718v1] **
|
---|
15 | ** HepMC: Comput. Phys. Commun.134 (2001) 41 **
|
---|
16 | ** **
|
---|
17 | ** ------------------------------------------------------------------ **
|
---|
18 | ** **
|
---|
19 | ** Main authors: **
|
---|
20 | ** ------------- **
|
---|
21 | ** **
|
---|
22 | ** Severine Ovyn Xavier Rouby **
|
---|
23 | ** severine.ovyn@uclouvain.be xavier.rouby@cern **
|
---|
24 | ** **
|
---|
25 | ** Center for Particle Physics and Phenomenology (CP3) **
|
---|
26 | ** Universite catholique de Louvain (UCL) **
|
---|
27 | ** Louvain-la-Neuve, Belgium **
|
---|
28 | ** **
|
---|
29 | ** Copyright (C) 2008-2009, **
|
---|
30 | ** All rights reserved. **
|
---|
31 | ** **
|
---|
32 | ***********************************************************************/
|
---|
33 |
|
---|
34 | #include <iostream>
|
---|
35 | #include <fstream>
|
---|
36 | #include <cstdlib> // for exit()
|
---|
37 |
|
---|
38 | #include "TLorentzVector.h"
|
---|
39 |
|
---|
40 | #include "ExRootTreeWriter.h"
|
---|
41 | #include "BlockClasses.h"
|
---|
42 | #include "LHEFConverter.h"
|
---|
43 | #include "SmearUtil.h"
|
---|
44 | #include "PdgParticle.h"
|
---|
45 |
|
---|
46 | using namespace std;
|
---|
47 |
|
---|
48 |
|
---|
49 | //------------------------------------------------------------------------------
|
---|
50 |
|
---|
51 | void LHEFConverter::AnalyseEvent(LHEF::Reader *reader, ExRootTreeBranch *branch, const Long64_t eventNumber)
|
---|
52 | {
|
---|
53 | const LHEF::HEPEUP &hepeup = reader->hepeup;
|
---|
54 |
|
---|
55 | TRootLHEFEvent *element;
|
---|
56 |
|
---|
57 | element = (TRootLHEFEvent*) branch->NewEntry();
|
---|
58 |
|
---|
59 | element->Number = eventNumber;
|
---|
60 | element->Nparticles = hepeup.NUP;
|
---|
61 | element->ProcessID = hepeup.IDPRUP;
|
---|
62 | element->Weight = hepeup.XWGTUP;
|
---|
63 | element->ScalePDF = hepeup.SCALUP;
|
---|
64 | element->CouplingQED = hepeup.AQEDUP;
|
---|
65 | element->CouplingQCD = hepeup.AQCDUP;
|
---|
66 | }
|
---|
67 |
|
---|
68 | //---------------------------------------------------------------------------
|
---|
69 |
|
---|
70 | void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
|
---|
71 | {
|
---|
72 | const LHEF::HEPEUP &hepeup = reader->hepeup;
|
---|
73 | Double_t signPz;
|
---|
74 | TLorentzVector momentum;
|
---|
75 | TRootLHEFParticle *element;
|
---|
76 |
|
---|
77 | for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
|
---|
78 | {
|
---|
79 | element = (TRootLHEFParticle*) branch->NewEntry();
|
---|
80 |
|
---|
81 | element->PID = hepeup.IDUP[particle];
|
---|
82 | element->Status = hepeup.ISTUP[particle];
|
---|
83 | element->M1 = hepeup.MOTHUP[particle].first;
|
---|
84 | element->M2 = hepeup.MOTHUP[particle].second;
|
---|
85 | element->D1 = 0;
|
---|
86 | element->D2 = 0;
|
---|
87 | element->ColorLine1 = hepeup.ICOLUP[particle].first;
|
---|
88 | element->ColorLine2 = hepeup.ICOLUP[particle].second;
|
---|
89 | element->Px = hepeup.PUP[particle][0];
|
---|
90 | element->Py = hepeup.PUP[particle][1];
|
---|
91 | element->Pz = hepeup.PUP[particle][2];
|
---|
92 | element->E = hepeup.PUP[particle][3];
|
---|
93 | //element->M = hepeup.PUP[particle][4];
|
---|
94 | //element->Charge = ChargeVal(element->PID);
|
---|
95 |
|
---|
96 | PdgParticle pdg_part = pdg_table[element->PID];
|
---|
97 | element->Charge = pdg_part.charge();
|
---|
98 | element->M = pdg_part.mass();
|
---|
99 |
|
---|
100 | momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
|
---|
101 | element->PT = momentum.Perp();
|
---|
102 | signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
|
---|
103 | //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
|
---|
104 | element->Eta = element->PT < 1E-6 ? signPz*999.9 : momentum.Eta();
|
---|
105 | element->Phi = momentum.Phi();
|
---|
106 | element->Rapidity = element->PT < 1E-6 ? signPz*999.9 : momentum.Rapidity();
|
---|
107 | element->LifeTime = hepeup.VTIMUP[particle];
|
---|
108 | element->Spin = hepeup.SPINUP[particle];
|
---|
109 | }
|
---|
110 | }
|
---|
111 |
|
---|
112 |
|
---|
113 | //---------------------------------------------------------------------------
|
---|
114 | /*
|
---|
115 | void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
|
---|
116 | {
|
---|
117 | const LHEF::HEPEUP &hepeup = reader->hepeup;
|
---|
118 | Double_t signPz;
|
---|
119 | TLorentzVector momentum;
|
---|
120 | TRootC::GenParticle *element;
|
---|
121 |
|
---|
122 | for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
|
---|
123 | {
|
---|
124 | element = (TRootC::GenParticle*) branch->NewEntry();
|
---|
125 |
|
---|
126 | element->PID = hepeup.IDUP[particle];
|
---|
127 | element->Status = hepeup.ISTUP[particle];
|
---|
128 | element->M1 = hepeup.MOTHUP[particle].first;
|
---|
129 | element->M2 = hepeup.MOTHUP[particle].second;
|
---|
130 | //element->ColorLine1 = hepeup.ICOLUP[particle].first;
|
---|
131 | //element->ColorLine2 = hepeup.ICOLUP[particle].second;
|
---|
132 | element->D1 = 0;
|
---|
133 | element->D2 = 0;
|
---|
134 |
|
---|
135 | PdgParticle pdg_part = pdg_table[element->PID];
|
---|
136 | element->Charge = pdg_part.charge();
|
---|
137 | element->M = pdg_part.mass();
|
---|
138 | //element->Charge = ChargeVal(element->PID);
|
---|
139 | //element->M = hepeup.PUP[particle][4];
|
---|
140 |
|
---|
141 | element->E = hepeup.PUP[particle][3];
|
---|
142 | element->Px = hepeup.PUP[particle][0];
|
---|
143 | element->Py = hepeup.PUP[particle][1];
|
---|
144 | element->Pz = hepeup.PUP[particle][2];
|
---|
145 |
|
---|
146 | momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
|
---|
147 | element->PT = momentum.Perp();
|
---|
148 | signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
|
---|
149 | //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
|
---|
150 | element->Eta = element->PT < 1E-6 ? signPz*999.9 : momentum.Eta();
|
---|
151 | element->Phi = momentum.Phi();
|
---|
152 | // element->Rapidity = element->PT < 1E-6 ? signPz*999.9 : momentum.Rapidity();
|
---|
153 | // element->LifeTime = hepeup.VTIMUP[particle];
|
---|
154 | // element->Spin = hepeup.SPINUP[particle];
|
---|
155 | }
|
---|
156 | }
|
---|
157 | */
|
---|
158 |
|
---|
159 | LHEFConverter::~LHEFConverter()
|
---|
160 | {
|
---|
161 | }
|
---|
162 |
|
---|
163 | //------------------------------------------------------------------------------
|
---|
164 | LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg,Nevents) {
|
---|
165 |
|
---|
166 | ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
|
---|
167 |
|
---|
168 | // generated event from LHEF
|
---|
169 | ExRootTreeBranch *branchEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
|
---|
170 |
|
---|
171 | // generated partons from LHEF
|
---|
172 | ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootLHEFParticle::Class());
|
---|
173 | //ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());
|
---|
174 |
|
---|
175 | int nevt_already_processed=0;
|
---|
176 |
|
---|
177 | // Open a stream connected to an event file:
|
---|
178 | ifstream infile(inputFileList.c_str());
|
---|
179 | string filename;
|
---|
180 | if(!infile.is_open()) {
|
---|
181 | cerr << left << setw(30) <<"** ERROR: Can't open "<<""
|
---|
182 | << left << setw(20) << inputFileList <<""
|
---|
183 | << right << setw(19) <<"for input **"<<""<<endl;
|
---|
184 | exit(1);
|
---|
185 | }
|
---|
186 |
|
---|
187 | while(1) { // parsing the list of files
|
---|
188 |
|
---|
189 | infile >> filename;
|
---|
190 | if(!infile.good()) break;
|
---|
191 | if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
|
---|
192 |
|
---|
193 | ifstream checking_the_file(filename.c_str());
|
---|
194 | if(!checking_the_file.good())
|
---|
195 | {
|
---|
196 | cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
|
---|
197 | << left << setw(36) << filename <<""
|
---|
198 | << right << setw(3) <<" **"<<""<<endl;
|
---|
199 | continue;
|
---|
200 | }
|
---|
201 | else checking_the_file.close();
|
---|
202 |
|
---|
203 | // Create the Reader object:
|
---|
204 | LHEF::Reader *inputReader = new LHEF::Reader(filename);
|
---|
205 |
|
---|
206 | Long64_t allEntries = inputReader->getNumberOfEvents();
|
---|
207 | allEntries = (Nevt<0)?allEntries:min((int)allEntries,Nevt); // do not miss the "+1"
|
---|
208 |
|
---|
209 | if(allEntries > 0) {
|
---|
210 | // Loop over all events
|
---|
211 | Long64_t entry = 0;
|
---|
212 | while(inputReader->readEvent()) {
|
---|
213 | if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
|
---|
214 | treeWriter->Clear();
|
---|
215 | AnalyseEvent(inputReader, branchEvent, entry + 1);
|
---|
216 | AnalyseParticles(inputReader, branchParticle);
|
---|
217 | treeWriter->Fill();
|
---|
218 | ++entry;
|
---|
219 | ++nevt_already_processed;
|
---|
220 | if(allEntries<entry+1) break;
|
---|
221 | }
|
---|
222 | }
|
---|
223 | }
|
---|
224 | treeWriter->Write();
|
---|
225 | delete treeWriter;
|
---|
226 | //delete inputReader;
|
---|
227 | }
|
---|
228 |
|
---|