Fork me on GitHub

source: svn/trunk/src/HepMCConverter.cc@ 368

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

bug fix if no mother written in the input data. small code-cleaning
File size: 10.8 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \----------------------------------------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
12** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
13** FROG: [hep-ex/0901.2718v1] **
14** **
15** ------------------------------------------------------------------ **
16** **
17** Main authors: **
18** ------------- **
19** **
20** Severine Ovyn Xavier Rouby **
21** severine.ovyn@uclouvain.be xavier.rouby@cern **
22** **
23** Center for Particle Physics and Phenomenology (CP3) **
24** Universite catholique de Louvain (UCL) **
25** Louvain-la-Neuve, Belgium **
26** **
27** Copyright (C) 2008-2009, **
28** All rights reserved. **
29** **
30***********************************************************************/
31
32
33#include <iostream>
34#include <fstream>
35#include "TLorentzVector.h"
36#include "BlockClasses.h"
37
38#include "ExRootTreeWriter.h"
39#include "ExRootTreeBranch.h"
40#include "HepMCConverter.h"
41
42#include "GenParticle.h"
43#include "GenVertex.h"
44#include "IO_Ascii.h"
45#include "IO_GenEvent.h"
46
47//-------------------------------------------------------------------------
48int HepMCConverter::find_in_map( const std::map<HepMC::GenParticle*,int>& m, HepMC::GenParticle *p) const
49{
50 std::map<HepMC::GenParticle*,int>::const_iterator iter = m.find(p);
51 return (iter == m.end()) ? 0 : iter->second;
52}
53
54//--------------------------------------------------------------------------
55void HepMCConverter::ReadStats()
56{
57
58 unsigned int particle_counter=0;
59 index_to_particle.reserve(evt->particles_size());
60 index_to_particle[0] = 0;
61 HepMC::GenEvent::vertex_const_iterator v;
62 for (v = evt->vertices_begin(); v != evt->vertices_end(); ++v )
63 {
64 // making a list of incoming particles of the vertices
65 // so that the mother indices in HEPEVT can be filled properly
66 HepMC::GenVertex::particles_out_const_iterator p1;
67 for (p1 = (*v)->particles_in_const_begin();p1 != (*v)->particles_in_const_end(); ++p1 )
68 {
69
70 ++particle_counter;
71 //particle_counter can be very large for heavy ions
72 if(particle_counter >= index_to_particle.size() )
73 {
74 //make it large enough to hold up to this index
75 index_to_particle.resize(particle_counter+1);
76 }
77 index_to_particle[particle_counter] = *p1;
78 particle_to_index[*p1] = particle_counter;
79 }
80 // daughters are entered only if they aren't a mother of
81 // another vertex
82 HepMC::GenVertex::particles_out_const_iterator p2;
83 for (p2 = (*v)->particles_out_const_begin();p2 != (*v)->particles_out_const_end(); ++p2)
84 {
85 if (!(*p2)->end_vertex())
86 {
87 ++particle_counter;
88 //particle_counter can be very large for heavy ions
89 if(particle_counter >= index_to_particle.size() )
90 {
91 //make it large enough to hold up to this index
92 index_to_particle.resize(particle_counter+1);
93 }
94 index_to_particle[particle_counter] = *p2;
95 particle_to_index[*p2] = particle_counter;
96 }
97 }
98
99 }
100}
101
102
103//-------------------------------------------------------------------------
104void HepMCConverter::getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2, int &status, int &pid, int j) const
105{
106 if (!evt)
107 {
108 cout << "HepMCFileReader: Got no event :-( Game over already ?" <<endl;
109 }
110 else
111 {
112 status = index_to_particle[j]->status();
113 pid = index_to_particle[j]->pdg_id();
114 if ( index_to_particle[j]->production_vertex() )
115 {
116 int num_mothers = index_to_particle[j]->production_vertex()->particles_in_size();
117 if (num_mothers ==0) {
118 mo1 = 0;
119 mo2 = 0;
120 }
121 else {
122 int first_mother = find_in_map( particle_to_index,*(index_to_particle[j]->production_vertex()->particles_in_const_begin()));
123 int last_mother = first_mother + num_mothers - 1;
124 if ( first_mother == 0 ) last_mother = 0;
125 mo1=first_mother;
126 mo2=last_mother;
127 } // if num_mothers !=0
128 }
129 else // no data on production_vertex
130 {
131 mo1 =0;
132 mo2 =0;
133 }
134 if (index_to_particle[j]->end_vertex())
135 {
136 //find # of 1. daughter
137 int first_daughter = find_in_map( particle_to_index,*(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
138 //cout <<"first_daughter "<< first_daughter << "num_daughters " << num_daughters << endl;
139 HepMC::GenVertex::particle_iterator ic;
140 int last_daughter=0;
141 //find # of last daughter
142 for (ic = index_to_particle[j]->end_vertex()->particles_begin(HepMC::children);ic != index_to_particle[j]->end_vertex()->particles_end(HepMC::children); ++ic)
143 last_daughter= find_in_map( particle_to_index,*ic);
144
145 if (first_daughter== 0) last_daughter = 0;
146 da1=first_daughter;
147 da2=last_daughter;
148 }
149 else
150 {
151 da1=0;
152 da2=0;
153 }
154 }
155}
156
157
158
159//---------------------------------------------------------------------------
160
161void HepMCConverter::AnalyseEvent(ExRootTreeBranch *branch,HepMC::GenEvent& evt,const Long64_t eventNumber)
162{
163 TRootLHEFEvent *element;
164
165 element = static_cast<TRootLHEFEvent*>(branch->NewEntry());
166 element->Number = eventNumber;
167 element->Nparticles = evt.particles_size();
168 element->ProcessID = evt.signal_process_id();
169 // element->Weight = hepeup.XWGTUP;
170 element->ScalePDF = evt.event_scale();
171 element->CouplingQED = evt.alphaQED();
172 element->CouplingQCD = evt.alphaQCD();
173
174}
175
176//---------------------------------------------------------------------------
177
178void HepMCConverter::AnalyseParticles(ExRootTreeBranch *branch, const HepMC::GenEvent& evt)
179{
180 TRootC::GenParticle *element;
181
182 TLorentzVector momentum;
183 Double_t signPz;
184
185 ReadStats();
186 for(int n=1; n<=evt.particles_size(); n++)
187 {
188 getStatsFromTuple( mo1,mo2,da1,da2,status,pid,n);
189
190 element = static_cast<TRootC::GenParticle*>(branch->NewEntry());
191
192 element->PID = pid;
193 element->Status = status;
194 element->M1 = mo1;
195 element->M2 = mo2;
196 element->D1 = da1;
197 element->D2 = da2;
198 element->Charge = pid;
199
200 element->E = index_to_particle[n]->momentum().e();
201 element->Px = index_to_particle[n]->momentum().px();
202 element->Py = index_to_particle[n]->momentum().py();
203 element->Pz = index_to_particle[n]->momentum().pz();
204 element->M = index_to_particle[n]->momentum().m();
205
206 float PT = sqrt(pow(element->Px,2)+pow(element->Py,2));
207 element->PT = PT;
208
209 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
210 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
211 //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
212 element->Eta = element->PT < 1e-6 ? signPz*999.9 : momentum.Eta();
213 element->Phi = index_to_particle[n]->momentum().phi();
214
215 //In particle at vertex
216 HepMC::GenVertex* vrtI = (index_to_particle[n])->production_vertex();
217 HepMC::GenVertex::particles_in_const_iterator partI;
218
219 if(vrtI)
220 {
221 element->T = vrtI->position().t();
222 element->X = vrtI->position().x();
223 element->Y = vrtI->position().y();
224 element->Z = vrtI->position().z();
225 }
226 else
227 {
228 element->T = 0.;
229 element->X = 0.;
230 element->Y = 0.;
231 element->Z = 0.;
232 }
233 }
234}
235
236//
237
238//------------------------------------------------------------------------------
239
240HepMCConverter::~HepMCConverter()
241{
242 delete evt;
243
244 /*cout << "delete index to particle" << endl;
245 std::vector<HepMC::GenParticle*>::iterator i;
246 / for (i=index_to_particle.begin();i != index_to_particle.end(); i++) delete (*i);
247 cout << "done" << endl;
248
249 std::map<HepMC::GenParticle*,int>::iterator j;
250 for (j=particle_to_index.begin(); j != particle_to_index.end(); j++) delete j->first;
251 */
252
253}
254
255//------------------------------------------------------------------------------
256
257HepMCConverter::HepMCConverter(const string& inputFileList, const string& outputFileName, const int& Nevents) : DataConverter(Nevents)
258{
259
260 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
261 // information about generated event
262 ExRootTreeBranch *branchGenEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
263 // generated particles from HEPEVT
264 ExRootTreeBranch *branchGenParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());
265
266 // Open a stream connected to an event file:
267 ifstream infile(inputFileList.c_str());
268 string filename;
269 if(!infile.is_open()) {
270 cerr << left << setw(30) <<"** ERROR: Can't open "<<""
271 << left << setw(20) << inputFileList <<""
272 << right << setw(19) <<"for input **"<<""<<endl;
273 exit(1);
274 }
275
276 Long64_t entry = 0;
277 while(1)
278 {
279 infile >> filename;
280 if(!infile.good()) break;
281 ifstream checking_the_file(filename.c_str());
282 if(!checking_the_file.good())
283 {
284 cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
285 << left << setw(20) << filename <<""
286 << right << setw(19) <<"for input **"<<""<<endl;
287 continue;
288 }
289 else checking_the_file.close();
290
291 //Open the file
292 HepMC::IO_GenEvent ascii_in(filename,std::ios::in);
293 // get the first event
294 evt = ascii_in.read_next_event();
295
296 while ( evt )
297 {
298 treeWriter->Clear();
299 AnalyseEvent(branchGenEvent, *evt,entry+1);
300 AnalyseParticles(branchGenParticle, *evt);
301 delete evt;
302 // read the next event
303 ascii_in >> evt;
304 treeWriter->Fill();
305 ++entry;
306 }
307 }
308
309 treeWriter->Write();
310 delete treeWriter;
311
312}
313
314
315
316
Note: See TracBrowser for help on using the repository browser.