CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HepMCFileReader.cc
Go to the documentation of this file.
1 
9 #include <iostream>
10 #include <iomanip>
11 #include <string>
12 
13 #include "HepMC/GenEvent.h"
14 #include "HepMC/GenParticle.h"
15 #include "HepMC/IO_GenEvent.h"
16 
17 #include "CLHEP/Units/GlobalPhysicalConstants.h"
18 
22 
23 using namespace std;
24 
25 
27 
28 
29 //-------------------------------------------------------------------------
31 {
32  // Implement a HepMCFileReader singleton.
33 
34  if (instance_== nullptr) {
35  instance_ = new HepMCFileReader();
36  }
37  return instance_;
38 }
39 
40 
41 //-------------------------------------------------------------------------
43  evt_(nullptr), input_(nullptr)
44 {
45  // Default constructor.
46  if (instance_ == nullptr) {
47  instance_ = this;
48  } else {
49  edm::LogError("HepMCFileReader") << "Constructing a new instance";
50  }
51 }
52 
53 
54 //-------------------------------------------------------------------------
56 {
57  edm::LogInfo("HepMCFileReader") << "Destructing HepMCFileReader";
58 
59  instance_=nullptr;
60  delete input_.get();
61 }
62 
63 
64 //-------------------------------------------------------------------------
66 {
67  if (isInitialized()) {
68  edm::LogError("HepMCFileReader") << "Was already initialized... reinitializing";
69  delete input_.get();
70  }
71 
72  edm::LogInfo("HepMCFileReader") << "Opening file" << filename << "using HepMC::IO_GenEvent";
73  input_ = new HepMC::IO_GenEvent(filename.c_str(), std::ios::in);
74 
75  if (rdstate() == std::ios::failbit) {
76  throw cms::Exception("FileNotFound", "HepMCFileReader::initialize()")
77  << "File " << filename << " was not found.\n";
78  }
79 }
80 
81 
82 //-------------------------------------------------------------------------
84 {
85  // work around a HepMC IO_ inheritence shortfall
86 
87  HepMC::IO_GenEvent const* p = dynamic_cast<HepMC::IO_GenEvent const*>(input());
88  if (p) return p->rdstate();
89 
90  return std::ios::failbit;
91 }
92 
93 
94 //-------------------------------------------------------------------------
96 {
97  evt_ = input_->read_next_event();
98  if (evt_) {
99  edm::LogInfo("HepMCFileReader") << "| --- Event Nr. " << evt_->event_number()
100  << " with " << evt_->particles_size() << " particles --- !";
101  ReadStats();
102  // printHepMcEvent();
103  // printEvent();
104  } else {
105  edm::LogInfo("HepMCFileReader") << "Got no event" <<endl;
106  }
107 
108  return evt_ != nullptr;
109 }
110 
111 
112 //-------------------------------------------------------------------------
114 {
115  return true;
116 }
117 
118 
119 //-------------------------------------------------------------------------
121 {
122  if (evt_ != nullptr) evt_->print();
123  return true;
124 }
125 
126 
127 //-------------------------------------------------------------------------
129 {
131  return evt_;
132 }
133 
134 
135 //-------------------------------------------------------------------------
136 // Print out in old CMKIN style for comparisons
138  int mo1=0,mo2=0,da1=0,da2=0,status=0,pid=0;
139  if (evt_ != nullptr) {
140  cout << "---#-------pid--st---Mo1---Mo2---Da1---Da2------px------py------pz-------E-";
141  cout << "------m---------x---------y---------z---------t-";
142  cout << endl;
143  cout.setf(ios::right, ios::adjustfield);
144  for(int n=1; n<=evt_->particles_size(); n++) {
146  getStatsFromTuple( mo1,mo2,da1,da2,status,pid,n);
147  cout << setw(4) << n
148  << setw(8) << pid
149  << setw(5) << status
150  << setw(6) << mo1
151  << setw(6) << mo2
152  << setw(6) << da1
153  << setw(6) << da2;
154  cout.setf(ios::fixed, ios::floatfield);
155  cout.setf(ios::right, ios::adjustfield);
156  cout << setw(10) << setprecision(2) << g->momentum().x();
157  cout << setw(8) << setprecision(2) << g->momentum().y();
158  cout << setw(10) << setprecision(2) << g->momentum().z();
159  cout << setw(8) << setprecision(2) << g->momentum().t();
160  cout << setw(8) << setprecision(2) << g->generatedMass();
161  // tau=L/(gamma*beta*c)
162  if (g->production_vertex() != nullptr && g->end_vertex() != nullptr && status == 2) {
163  cout << setw(10) << setprecision(2) <<g->production_vertex()->position().x();
164  cout << setw(10) << setprecision(2) <<g->production_vertex()->position().y();
165  cout << setw(10) << setprecision(2) <<g->production_vertex()->position().z();
166 
167  double xm = g->production_vertex()->position().x();
168  double ym = g->production_vertex()->position().y();
169  double zm = g->production_vertex()->position().z();
170  double xd = g->end_vertex()->position().x();
171  double yd = g->end_vertex()->position().y();
172  double zd = g->end_vertex()->position().z();
173  double decl = sqrt((xd-xm)*(xd-xm)+(yd-ym)*(yd-ym)+(zd-zm)*(zd-zm));
174  double labTime = decl/c_light;
175  // convert lab time to proper time
176  double properTime = labTime/g->momentum().rho()*(g->generatedMass() );
177  // set the proper time in nanoseconds
178  cout << setw(8) << setprecision(2) <<properTime;
179  }
180  else{
181  cout << setw(10) << setprecision(2) << 0.0;
182  cout << setw(10) << setprecision(2) << 0.0;
183  cout << setw(10) << setprecision(2) << 0.0;
184  cout << setw(8) << setprecision(2) << 0.0;
185  }
186  cout << endl;
187  }
188  } else {
189  cout << " HepMCFileReader: No event available !" << endl;
190  }
191 }
192 
193 
194 //-------------------------------------------------------------------------
196 {
197  unsigned int particle_counter=0;
198  index_to_particle.reserve(evt_->particles_size()+1);
199  index_to_particle[0] = nullptr;
200  for (HepMC::GenEvent::vertex_const_iterator v = evt_->vertices_begin();
201  v != evt_->vertices_end(); ++v ) {
202  // making a list of incoming particles of the vertices
203  // so that the mother indices in HEPEVT can be filled properly
204  for (HepMC::GenVertex::particles_in_const_iterator p1 = (*v)->particles_in_const_begin();
205  p1 != (*v)->particles_in_const_end(); ++p1 ) {
206  ++particle_counter;
207  //particle_counter can be very large for heavy ions
208  if(particle_counter >= index_to_particle.size() ) {
209  //make it large enough to hold up to this index
210  index_to_particle.resize(particle_counter+1);
211  }
212  index_to_particle[particle_counter] = *p1;
213  particle_to_index[*p1] = particle_counter;
214  }
215  // daughters are entered only if they aren't a mother of
216  // another vertex
217  for (HepMC::GenVertex::particles_out_const_iterator p2 = (*v)->particles_out_const_begin();
218  p2 != (*v)->particles_out_const_end(); ++p2) {
219  if (!(*p2)->end_vertex()) {
220  ++particle_counter;
221  //particle_counter can be very large for heavy ions
222  if(particle_counter >= index_to_particle.size() ) {
223  //make it large enough to hold up to this index
224  index_to_particle.resize(particle_counter+1);
225  }
226  index_to_particle[particle_counter] = *p2;
227  particle_to_index[*p2] = particle_counter;
228  }
229  }
230  }
231 }
232 
233 
234 //-------------------------------------------------------------------------
235 void HepMCFileReader::getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2,
236  int &status, int &pid, int j) const
237 {
238  if (!evt_) {
239  cout << "HepMCFileReader: Got no event :-( Game over already ?" <<endl;
240  } else {
241  status = index_to_particle[j]->status();
242  pid = index_to_particle[j]->pdg_id();
243  if ( index_to_particle[j]->production_vertex() ) {
244 
245  //HepLorentzVector p = index_to_particle[j]->
246  //production_vertex()->position();
247 
248  int num_mothers = index_to_particle[j]->production_vertex()->
249  particles_in_size();
250  int first_mother = find_in_map( particle_to_index,
251  *(index_to_particle[j]->
252  production_vertex()->
253  particles_in_const_begin()));
254  int last_mother = first_mother + num_mothers - 1;
255  if ( first_mother == 0 ) last_mother = 0;
256  mo1=first_mother;
257  mo2=last_mother;
258  } else {
259  mo1 =0;
260  mo2 =0;
261  }
262  if (!index_to_particle[j]->end_vertex()) {
263  //find # of 1. daughter
264  int first_daughter = find_in_map( particle_to_index,
265  *(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
266  //cout <<"first_daughter "<< first_daughter << "num_daughters " << num_daughters << endl;
267  HepMC::GenVertex::particle_iterator ic;
268  int last_daughter=0;
269  //find # of last daughter
270  for (ic = index_to_particle[j]->end_vertex()->particles_begin(HepMC::children);
271  ic != index_to_particle[j]->end_vertex()->particles_end(HepMC::children); ++ic)
272  last_daughter= find_in_map( particle_to_index,*ic);
273 
274  if (first_daughter== 0) last_daughter = 0;
275  da1=first_daughter;
276  da2=last_daughter;
277  } else {
278  da1=0;
279  da2=0;
280  }
281  }
282 }
283 
284 
285 //-------------------------------------------------------------------------
286 int HepMCFileReader::find_in_map( const std::map<HepMC::GenParticle*,int>& m,
287  HepMC::GenParticle *p) const
288 {
289  std::map<HepMC::GenParticle*,int>::const_iterator iter = m.find(p);
290  return (iter == m.end()) ? 0 : iter->second;
291 }
292 
void printEvent() const
static HepMCFileReader * instance_
virtual void ReadStats()
std::map< HepMC::GenParticle *, int > particle_to_index
int rdstate() const
virtual void initialize(const std::string &filename)
#define nullptr
virtual bool setEvent(int event)
HepMC::IO_BaseClass const * input() const
std::vector< HepMC::GenParticle * > index_to_particle
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
virtual ~HepMCFileReader()
static HepMCFileReader * instance()
virtual void getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2, int &status, int &pid, int j) const
int find_in_map(const std::map< HepMC::GenParticle *, int > &m, HepMC::GenParticle *p) const
T sqrt(T t)
Definition: SSEVec.h:18
int j
Definition: DBlmapReader.cc:9
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
bool isInitialized() const
double p2[4]
Definition: TauolaWrapper.h:90
element_type const * get() const
virtual bool printHepMcEvent() const
virtual bool readCurrentEvent()
HepMC::GenEvent * fillCurrentEventData()
tuple pid
Definition: sysUtil.py:22
string const
Definition: compareJSON.py:14
edm::propagate_const< HepMC::IO_BaseClass * > input_
double p1[4]
Definition: TauolaWrapper.h:89
tuple filename
Definition: lut2db_cfg.py:20
tuple cout
Definition: gather_cfg.py:145
tuple status
Definition: ntuplemaker.py:245
edm::propagate_const< HepMC::GenEvent * > evt_