Fork me on GitHub

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

Last change on this file since 362 was 362, checked in by severine ovyn, 16 years ago

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