Fork me on GitHub

Ticket #44: HepMCConverter.cc

File HepMCConverter.cc, 12.1 KB (added by muanza, 13 years ago)

Delphes_V_2.0.2/src/HepMCConverter.cc

Line 
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//-------------------------------------------------------------------------
50int 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//--------------------------------------------------------------------------
57void 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//-------------------------------------------------------------------------
108void 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
169void 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
193void 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
260HepMCConverter::~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
277HepMCConverter::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