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 "TLorentzVector.h"
|
---|
37 | #include "BlockClasses.h"
|
---|
38 |
|
---|
39 | #include "PdgParticle.h"
|
---|
40 | #include "ExRootTreeWriter.h"
|
---|
41 | #include "ExRootTreeBranch.h"
|
---|
42 | #include "HepMCConverter.h"
|
---|
43 |
|
---|
44 | #include "GenParticle.h"
|
---|
45 | #include "GenVertex.h"
|
---|
46 | #include "IO_AsciiParticles.h"
|
---|
47 | #include "IO_GenEvent.h"
|
---|
48 |
|
---|
49 | //-------------------------------------------------------------------------
|
---|
50 | int HepMCConverter::find_in_map( const std::map<HepMC::GenParticle*,int>& m, HepMC::GenParticle *p) const
|
---|
51 | {
|
---|
52 | std::map<HepMC::GenParticle*,int>::const_iterator iter = m.find(p);
|
---|
53 | return (iter == m.end()) ? 0 : iter->second;
|
---|
54 | }
|
---|
55 |
|
---|
56 | //--------------------------------------------------------------------------
|
---|
57 | void HepMCConverter::ReadStats()
|
---|
58 | {
|
---|
59 |
|
---|
60 | unsigned int particle_counter=0;
|
---|
61 | index_to_particle.clear();
|
---|
62 | particle_to_index.clear();
|
---|
63 | index_to_particle.reserve(evt->particles_size());
|
---|
64 | index_to_particle[0] = 0;
|
---|
65 | HepMC::GenEvent::vertex_const_iterator v;
|
---|
66 | for (v = evt->vertices_begin(); v != evt->vertices_end(); ++v )
|
---|
67 | {
|
---|
68 | // making a list of incoming particles of the vertices
|
---|
69 | // so that the mother indices in HEPEVT can be filled properly
|
---|
70 | HepMC::GenVertex::particles_out_const_iterator p1;
|
---|
71 | for (p1 = (*v)->particles_in_const_begin();p1 != (*v)->particles_in_const_end(); ++p1 )
|
---|
72 | {
|
---|
73 |
|
---|
74 | ++particle_counter;
|
---|
75 | //particle_counter can be very large for heavy ions
|
---|
76 | if(particle_counter >= index_to_particle.size() )
|
---|
77 | {
|
---|
78 | //make it large enough to hold up to this index
|
---|
79 | index_to_particle.resize(particle_counter+1);
|
---|
80 | }
|
---|
81 | index_to_particle[particle_counter] = *p1;
|
---|
82 | particle_to_index[*p1] = particle_counter;
|
---|
83 | }
|
---|
84 | // daughters are entered only if they aren't a mother of
|
---|
85 | // another vertex
|
---|
86 | HepMC::GenVertex::particles_out_const_iterator p2;
|
---|
87 | for (p2 = (*v)->particles_out_const_begin();p2 != (*v)->particles_out_const_end(); ++p2)
|
---|
88 | {
|
---|
89 | if (!(*p2)->end_vertex())
|
---|
90 | {
|
---|
91 | ++particle_counter;
|
---|
92 | //particle_counter can be very large for heavy ions
|
---|
93 | if(particle_counter >= index_to_particle.size() )
|
---|
94 | {
|
---|
95 | //make it large enough to hold up to this index
|
---|
96 | index_to_particle.resize(particle_counter+1);
|
---|
97 | }
|
---|
98 | index_to_particle[particle_counter] = *p2;
|
---|
99 | particle_to_index[*p2] = particle_counter;
|
---|
100 | }
|
---|
101 | }
|
---|
102 |
|
---|
103 | }
|
---|
104 | }
|
---|
105 |
|
---|
106 |
|
---|
107 | //-------------------------------------------------------------------------
|
---|
108 | void HepMCConverter::getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2, int &status, int &pid, int j) const
|
---|
109 | {
|
---|
110 | if (!evt)
|
---|
111 | {
|
---|
112 | cout << "HepMCFileReader: Got no event :-( Game over already ?" <<endl;
|
---|
113 | }
|
---|
114 | else
|
---|
115 | {
|
---|
116 | status = index_to_particle[j]->status();
|
---|
117 | pid = index_to_particle[j]->pdg_id();
|
---|
118 | if ( index_to_particle[j]->production_vertex() )
|
---|
119 | {
|
---|
120 | int num_mothers = index_to_particle[j]->production_vertex()->particles_in_size();
|
---|
121 | if (num_mothers ==0) {
|
---|
122 | mo1 = 0;
|
---|
123 | mo2 = 0;
|
---|
124 | }
|
---|
125 | else {
|
---|
126 | int first_mother = find_in_map( particle_to_index,*(index_to_particle[j]->production_vertex()->particles_in_const_begin()));
|
---|
127 | int last_mother = first_mother + num_mothers - 1;
|
---|
128 | if ( first_mother == 0 ) last_mother = 0;
|
---|
129 | mo1=first_mother;
|
---|
130 | mo2=last_mother;
|
---|
131 | } // if num_mothers !=0
|
---|
132 | }
|
---|
133 | else // no data on production_vertex
|
---|
134 | {
|
---|
135 | mo1 =0;
|
---|
136 | mo2 =0;
|
---|
137 | }
|
---|
138 | if (index_to_particle[j]->end_vertex())
|
---|
139 | {
|
---|
140 | // make sure first and last daughter are indeed the first and last
|
---|
141 | int first_daughter = find_in_map( particle_to_index,*(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
|
---|
142 | int last_daughter = 0;
|
---|
143 |
|
---|
144 | HepMC::GenVertex::particle_iterator ic;
|
---|
145 | for (ic = index_to_particle[j]->end_vertex()->particles_begin(HepMC::children);ic != index_to_particle[j]->end_vertex()->particles_end(HepMC::children); ++ic) {
|
---|
146 | int current_daughter = find_in_map( particle_to_index,*ic) ;
|
---|
147 | if (current_daughter < first_daughter)
|
---|
148 | first_daughter = current_daughter;
|
---|
149 | if (current_daughter > last_daughter)
|
---|
150 | last_daughter = current_daughter;
|
---|
151 | }
|
---|
152 | if (first_daughter== 0) last_daughter = 0;
|
---|
153 | da1=first_daughter;
|
---|
154 | da2=last_daughter;
|
---|
155 |
|
---|
156 | }
|
---|
157 | else
|
---|
158 | {
|
---|
159 | da1=0;
|
---|
160 | da2=0;
|
---|
161 | }
|
---|
162 | }
|
---|
163 | }
|
---|
164 |
|
---|
165 |
|
---|
166 |
|
---|
167 | //---------------------------------------------------------------------------
|
---|
168 |
|
---|
169 | void HepMCConverter::AnalyseEvent(ExRootTreeBranch *branch,HepMC::GenEvent& evt,const Long64_t eventNumber)
|
---|
170 | {
|
---|
171 | TRootLHEFEvent *element;
|
---|
172 |
|
---|
173 | element = static_cast<TRootLHEFEvent*>(branch->NewEntry());
|
---|
174 | element->Number = eventNumber;
|
---|
175 | element->Nparticles = evt.particles_size();
|
---|
176 | element->ProcessID = evt.signal_process_id();
|
---|
177 | // element->Weight = hepeup.XWGTUP;
|
---|
178 | element->ScalePDF = evt.event_scale();
|
---|
179 | element->CouplingQED = evt.alphaQED();
|
---|
180 | element->CouplingQCD = evt.alphaQCD();
|
---|
181 | //
|
---|
182 | element->ID1 = evt.pdf_info()->id1();
|
---|
183 | element->ID2 = evt.pdf_info()->id2();
|
---|
184 | element->x1 = evt.pdf_info()->x1();
|
---|
185 | element->x2 = evt.pdf_info()->x2();
|
---|
186 | element->Q_PDF = evt.pdf_info()->scalePDF();
|
---|
187 | element->xpdf1 = evt.pdf_info()->pdf1();
|
---|
188 | element->xpdf2 = evt.pdf_info()->pdf2();
|
---|
189 | }
|
---|
190 |
|
---|
191 | //---------------------------------------------------------------------------
|
---|
192 |
|
---|
193 | void HepMCConverter::AnalyseParticles(ExRootTreeBranch *branch, const HepMC::GenEvent& evt)
|
---|
194 | {
|
---|
195 | TRootC::GenParticle *element;
|
---|
196 |
|
---|
197 | TLorentzVector momentum;
|
---|
198 | Double_t signPz;
|
---|
199 |
|
---|
200 | ReadStats();
|
---|
201 | for(int n=1; n<=evt.particles_size(); n++)
|
---|
202 | {
|
---|
203 | getStatsFromTuple( mo1,mo2,da1,da2,status,pid,n);
|
---|
204 |
|
---|
205 | element = static_cast<TRootC::GenParticle*>(branch->NewEntry());
|
---|
206 |
|
---|
207 | element->PID = pid;
|
---|
208 | element->Status = status;
|
---|
209 | element->M1 = mo1 - 1; // added -1 as the numbering in the tree starts from 0
|
---|
210 | element->M2 = mo2 - 1;
|
---|
211 | element->D1 = da1 - 1;
|
---|
212 | element->D2 = da2 - 1;
|
---|
213 |
|
---|
214 | element->E = index_to_particle[n]->momentum().e();
|
---|
215 | element->Px = index_to_particle[n]->momentum().px();
|
---|
216 | element->Py = index_to_particle[n]->momentum().py();
|
---|
217 | element->Pz = index_to_particle[n]->momentum().pz();
|
---|
218 |
|
---|
219 |
|
---|
220 | //cout << "element->PID = " << pid << "\t";
|
---|
221 | //PdgParticle pdg_part(PdgID[pid]);
|
---|
222 | //element->M = index_to_particle[n]->momentum().m(); // this is the particle virtuality, not its rest mass
|
---|
223 | //element->M = pdg_part.mass();
|
---|
224 | //element->Charge = pdg_part.charge();
|
---|
225 | //cout << "element->M = " << element->M << " \t element->Charge = " << element->Charge << endl;
|
---|
226 |
|
---|
227 | element->PT = sqrt(pow(element->Px,2)+pow(element->Py,2));
|
---|
228 |
|
---|
229 | momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
|
---|
230 | signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
|
---|
231 | //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
|
---|
232 | element->Eta = element->PT < 1e-6 ? signPz*999.9 : momentum.Eta();
|
---|
233 | element->Phi = index_to_particle[n]->momentum().phi();
|
---|
234 |
|
---|
235 | //In particle at vertex
|
---|
236 | HepMC::GenVertex* vrtI = (index_to_particle[n])->production_vertex();
|
---|
237 | HepMC::GenVertex::particles_in_const_iterator partI;
|
---|
238 |
|
---|
239 | if(vrtI)
|
---|
240 | {
|
---|
241 | element->T = vrtI->position().t();
|
---|
242 | element->X = vrtI->position().x();
|
---|
243 | element->Y = vrtI->position().y();
|
---|
244 | element->Z = vrtI->position().z();
|
---|
245 | }
|
---|
246 | else
|
---|
247 | {
|
---|
248 | element->T = 0.;
|
---|
249 | element->X = 0.;
|
---|
250 | element->Y = 0.;
|
---|
251 | element->Z = 0.;
|
---|
252 | }
|
---|
253 | }
|
---|
254 | }
|
---|
255 |
|
---|
256 | //
|
---|
257 |
|
---|
258 | //------------------------------------------------------------------------------
|
---|
259 |
|
---|
260 | HepMCConverter::~HepMCConverter()
|
---|
261 | {
|
---|
262 | delete evt;
|
---|
263 |
|
---|
264 | /*cout << "delete index to particle" << endl;
|
---|
265 | std::vector<HepMC::GenParticle*>::iterator i;
|
---|
266 | / for (i=index_to_particle.begin();i != index_to_particle.end(); i++) delete (*i);
|
---|
267 | cout << "done" << endl;
|
---|
268 |
|
---|
269 | std::map<HepMC::GenParticle*,int>::iterator j;
|
---|
270 | for (j=particle_to_index.begin(); j != particle_to_index.end(); j++) delete j->first;
|
---|
271 | */
|
---|
272 |
|
---|
273 | }
|
---|
274 |
|
---|
275 | //------------------------------------------------------------------------------
|
---|
276 |
|
---|
277 | HepMCConverter::HepMCConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg,Nevents)
|
---|
278 | {
|
---|
279 |
|
---|
280 | ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
|
---|
281 | // information about generated event
|
---|
282 | ExRootTreeBranch *branchGenEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
|
---|
283 | // generated particles from HEPEVT
|
---|
284 | ExRootTreeBranch *branchGenParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());
|
---|
285 |
|
---|
286 | int nevt_already_processed=0;
|
---|
287 |
|
---|
288 | // Open a stream connected to an event file:
|
---|
289 | ifstream infile(inputFileList.c_str());
|
---|
290 | string filename;
|
---|
291 | if(!infile.is_open()) {
|
---|
292 | cerr << left << setw(30) <<"** ERROR: Can't open "<<""
|
---|
293 | << left << setw(20) << inputFileList <<""
|
---|
294 | << right << setw(19) <<"for input **"<<""<<endl;
|
---|
295 | exit(1);
|
---|
296 | }
|
---|
297 |
|
---|
298 | Long64_t entry = 0;
|
---|
299 | while(1)
|
---|
300 | {
|
---|
301 | infile >> filename;
|
---|
302 | if(!infile.good()) break; // end of listfile reached
|
---|
303 | if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
|
---|
304 |
|
---|
305 | ifstream checking_the_file(filename.c_str());
|
---|
306 | if(!checking_the_file.good())
|
---|
307 | {
|
---|
308 | cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
|
---|
309 | << left << setw(36) << filename <<""
|
---|
310 | << right << setw(3) <<" **"<<""<<endl;
|
---|
311 | continue;
|
---|
312 | }
|
---|
313 | else checking_the_file.close();
|
---|
314 |
|
---|
315 | //Open the file
|
---|
316 | HepMC::IO_GenEvent ascii_in(filename,std::ios::in);
|
---|
317 | // get the first event
|
---|
318 |
|
---|
319 | evt = ascii_in.read_next_event();
|
---|
320 |
|
---|
321 | while ( evt )
|
---|
322 | {
|
---|
323 | if(Nevt>0 && entry>=Nevt) break;
|
---|
324 | if(Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
|
---|
325 |
|
---|
326 | treeWriter->Clear();
|
---|
327 | AnalyseEvent(branchGenEvent, *evt,entry+1);
|
---|
328 | AnalyseParticles(branchGenParticle, *evt);
|
---|
329 | delete evt;
|
---|
330 | // read the next event
|
---|
331 | ascii_in >> evt;
|
---|
332 | treeWriter->Fill();
|
---|
333 | ++entry;
|
---|
334 | ++nevt_already_processed;
|
---|
335 | }
|
---|
336 | }
|
---|
337 |
|
---|
338 | treeWriter->Write();
|
---|
339 | delete treeWriter;
|
---|
340 |
|
---|
341 | }
|
---|
342 |
|
---|