Utilities/src/HepMCFileReader.cc

// $Id: HepMCFileReader.cc,v 1.1 2012/06/07 19:10:12 antoniov Exp $

/**  
*  See header file for a description of this class.
*
*
*  $Date: 2012/06/07 19:10:12 $
*  $Revision: 1.1 $
*  \author Jo. Weng  - CERN, Ph Division & Uni Karlsruhe
*/

#include <iostream>
#include <iomanip>
#include <string>

#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/IO_GenEvent.h"

#include "CLHEP/Units/GlobalPhysicalConstants.h"

#include "ForwardAnalysis/Utilities/interface/HepMCFileReader.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"

using namespace std;


HepMCFileReader *HepMCFileReader::instance_=0;


//-------------------------------------------------------------------------
HepMCFileReader *HepMCFileReader::instance()
{
  // Implement a HepMCFileReader singleton.

  if (instance_== 0) {
    instance_ = new HepMCFileReader();
  }
  return instance_;
}


//-------------------------------------------------------------------------
HepMCFileReader::HepMCFileReader() :
  evt_(0), input_(0)
{ 
  // Default constructor.
  if (instance_ == 0) {
    instance_ = this;  
  } else {
    edm::LogError("HepMCFileReader") << "Constructing a new instance";
  }
} 


//-------------------------------------------------------------------------
HepMCFileReader::~HepMCFileReader()
{
  edm::LogInfo("HepMCFileReader") << "Destructing HepMCFileReader";
  
  instance_=0;    
  delete input_;
}


//-------------------------------------------------------------------------
void HepMCFileReader::initialize(const string &filename)
{
  if (isInitialized()) {
    edm::LogError("HepMCFileReader") << "Was already initialized... reinitializing";
    delete input_;
  }

  edm::LogInfo("HepMCFileReader") << "Opening file" << filename << "using HepMC::IO_GenEvent";
  input_ = new HepMC::IO_GenEvent(filename.c_str(), std::ios::in);

  if (rdstate() == std::ios::failbit) {
    throw cms::Exception("FileNotFound", "HepMCFileReader::initialize()")
      << "File " << filename << " was not found.\n";
  }
}


//-------------------------------------------------------------------------
int HepMCFileReader::rdstate() const
{
  // work around a HepMC IO_ inheritence shortfall

  HepMC::IO_GenEvent *p = dynamic_cast<HepMC::IO_GenEvent*>(input_);
  if (p) return p->rdstate();

  return std::ios::failbit;
}


//-------------------------------------------------------------------------
bool HepMCFileReader::readCurrentEvent()
{
  evt_ = input_->read_next_event();
  if (evt_) { 
    edm::LogInfo("HepMCFileReader") << "| --- Event Nr. "  << evt_->event_number()
         << " with " << evt_->particles_size() << " particles --- !";  
    ReadStats();
    //  printHepMcEvent();
    //          printEvent();
  } else {
    edm::LogInfo("HepMCFileReader") << "Got no event" <<endl;
  }

  return evt_ != 0;
}


//-------------------------------------------------------------------------
bool HepMCFileReader::setEvent(int event)
{
  return true;
}


//-------------------------------------------------------------------------
bool HepMCFileReader::printHepMcEvent() const
{
  if (evt_ != 0) evt_->print(); 
  return true;
}


//-------------------------------------------------------------------------
HepMC::GenEvent * HepMCFileReader::fillCurrentEventData()
{
  readCurrentEvent();
  return evt_;
}


//-------------------------------------------------------------------------
// Print out in old CMKIN style for comparisons
void HepMCFileReader::printEvent() const {
  int mo1=0,mo2=0,da1=0,da2=0,status=0,pid=0;   
  if (evt_ != 0) {   
    cout << "---#-------pid--st---Mo1---Mo2---Da1---Da2------px------py------pz-------E-";
    cout << "------m---------x---------y---------z---------t-";
    cout << endl;
    cout.setf(ios::right, ios::adjustfield);
    for(int n=1; n<=evt_->particles_size(); n++) {  
      HepMC::GenParticle *g = index_to_particle[n];  
      getStatsFromTuple( mo1,mo2,da1,da2,status,pid,n);
      cout << setw(4) << n
      << setw(8) << pid
      << setw(5) << status
      << setw(6) << mo1    
      << setw(6) << mo2
      << setw(6) << da1
      << setw(6) << da2;
      cout.setf(ios::fixed, ios::floatfield);
      cout.setf(ios::right, ios::adjustfield);
      cout << setw(10) << setprecision(2) << g->momentum().x();
      cout << setw(8) << setprecision(2) << g->momentum().y();
      cout << setw(10) << setprecision(2) << g->momentum().z();
      cout << setw(8) << setprecision(2) << g->momentum().t();
      cout << setw(8) << setprecision(2) << g->generatedMass();       
      // tau=L/(gamma*beta*c) 
      if (g->production_vertex() != 0 && g->end_vertex() != 0 && status == 2) {
        cout << setw(10) << setprecision(2) <<g->production_vertex()->position().x();
        cout << setw(10) << setprecision(2) <<g->production_vertex()->position().y();
        cout << setw(10) << setprecision(2) <<g->production_vertex()->position().z();
        
        double xm = g->production_vertex()->position().x();
        double ym = g->production_vertex()->position().y();
        double zm = g->production_vertex()->position().z();
        double xd = g->end_vertex()->position().x();
        double yd = g->end_vertex()->position().y();
        double zd = g->end_vertex()->position().z();
        double decl = sqrt((xd-xm)*(xd-xm)+(yd-ym)*(yd-ym)+(zd-zm)*(zd-zm));
        double labTime = decl/c_light;
        // convert lab time to proper time
        double properTime =  labTime/g->momentum().rho()*(g->generatedMass() );
        // set the proper time in nanoseconds
        cout << setw(8) << setprecision(2) <<properTime;
      }
      else{
        cout << setw(10) << setprecision(2) << 0.0;
        cout << setw(10) << setprecision(2) << 0.0;
        cout << setw(10) << setprecision(2) << 0.0;
        cout << setw(8) << setprecision(2) << 0.0;  
      }
      cout << endl; 
    }
  } else {
    cout << " HepMCFileReader: No event available !" << endl;
  }
}


//-------------------------------------------------------------------------
void HepMCFileReader::ReadStats()
{
  unsigned int particle_counter=0;
  index_to_particle.reserve(evt_->particles_size()+1); 
  index_to_particle[0] = 0; 
  for (HepMC::GenEvent::vertex_const_iterator v = evt_->vertices_begin();
       v != evt_->vertices_end(); ++v ) {
    // making a list of incoming particles of the vertices
    // so that the mother indices in HEPEVT can be filled properly
    for (HepMC::GenVertex::particles_in_const_iterator p1 = (*v)->particles_in_const_begin();
         p1 != (*v)->particles_in_const_end(); ++p1 ) {
      ++particle_counter;
      //particle_counter can be very large for heavy ions
      if(particle_counter >= index_to_particle.size() ) {
        //make it large enough to hold up to this index
         index_to_particle.resize(particle_counter+1);
      }                 
      index_to_particle[particle_counter] = *p1;
      particle_to_index[*p1] = particle_counter;
    }   
    // daughters are entered only if they aren't a mother of
    // another vertex
    for (HepMC::GenVertex::particles_out_const_iterator p2 = (*v)->particles_out_const_begin();
         p2 != (*v)->particles_out_const_end(); ++p2) {
      if (!(*p2)->end_vertex()) {
        ++particle_counter;
        //particle_counter can be very large for heavy ions
        if(particle_counter  >= index_to_particle.size() ) {        
          //make it large enough to hold up to this index
          index_to_particle.resize(particle_counter+1);
        }                   
        index_to_particle[particle_counter] = *p2;
        particle_to_index[*p2] = particle_counter;
      }
    }
  } 
}


//-------------------------------------------------------------------------
void HepMCFileReader::getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2,
                                        int &status, int &pid, int j) const
{
  if (!evt_) {
    cout <<  "HepMCFileReader: Got no event :-(  Game over already  ?" <<endl; 
  } else {
    status =  index_to_particle[j]->status();
    pid = index_to_particle[j]->pdg_id();
    if ( index_to_particle[j]->production_vertex() ) {
      
          //HepLorentzVector p = index_to_particle[j]->
          //production_vertex()->position();
      
      int num_mothers = index_to_particle[j]->production_vertex()->
      particles_in_size();
      int first_mother = find_in_map( particle_to_index,
      *(index_to_particle[j]->
      production_vertex()->
      particles_in_const_begin()));
      int last_mother = first_mother + num_mothers - 1;
      if ( first_mother == 0 ) last_mother = 0;
      mo1=first_mother;
      mo2=last_mother;
    } else {
      mo1 =0;
      mo2 =0;
    }   
    if (!index_to_particle[j]->end_vertex()) {
      //find # of 1. daughter
      int first_daughter = find_in_map( particle_to_index,
      *(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
      //cout <<"first_daughter "<< first_daughter <<  "num_daughters " << num_daughters << endl; 
      HepMC::GenVertex::particle_iterator ic;
      int last_daughter=0;
      //find # of last daughter
      for (ic = index_to_particle[j]->end_vertex()->particles_begin(HepMC::children);
      ic != index_to_particle[j]->end_vertex()->particles_end(HepMC::children); ++ic) 
      last_daughter= find_in_map( particle_to_index,*ic);       
      
      if (first_daughter== 0) last_daughter = 0;
      da1=first_daughter;
      da2=last_daughter;
    } else {
      da1=0;
      da2=0;
    }
  } 
} 


//-------------------------------------------------------------------------
int HepMCFileReader::find_in_map( const std::map<HepMC::GenParticle*,int>& m, 
                                  HepMC::GenParticle *p) const
{
  std::map<HepMC::GenParticle*,int>::const_iterator iter = m.find(p);
  return (iter == m.end()) ? 0 : iter->second;
}

Revision:	1.1
Committed:	Mon Jun 11 12:59:57 2012 UTC (12 years, 10 months ago) by antoniov
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-01-01, V01-01-00, antoniov-forwardAnalysis-09Jul2012-v1, antoniov-forwardAnalysis-29Jun2012-v1, V01-00-00, antoniov-utilities-11Jun2012-v1, HEAD
Log Message:	update
#	User	Rev	Content
1	antoniov	1.1	// $Id: HepMCFileReader.cc,v 1.1 2012/06/07 19:10:12 antoniov Exp $
2
3			/**
4			* See header file for a description of this class.
5			*
6			*
7			* $Date: 2012/06/07 19:10:12 $
8			* $Revision: 1.1 $
9			* \author Jo. Weng - CERN, Ph Division & Uni Karlsruhe
10			*/
11
12			#include <iostream>
13			#include <iomanip>
14			#include <string>
15
16			#include "HepMC/GenEvent.h"
17			#include "HepMC/GenParticle.h"
18			#include "HepMC/IO_GenEvent.h"
19
20			#include "CLHEP/Units/GlobalPhysicalConstants.h"
21
22			#include "ForwardAnalysis/Utilities/interface/HepMCFileReader.h"
23			#include "FWCore/MessageLogger/interface/MessageLogger.h"
24			#include "FWCore/Utilities/interface/Exception.h"
25
26			using namespace std;
27
28
29			HepMCFileReader *HepMCFileReader::instance_=0;
30
31
32			//-------------------------------------------------------------------------
33			HepMCFileReader *HepMCFileReader::instance()
34			{
35			// Implement a HepMCFileReader singleton.
36
37			if (instance_== 0) {
38			instance_ = new HepMCFileReader();
39			}
40			return instance_;
41			}
42
43
44			//-------------------------------------------------------------------------
45			HepMCFileReader::HepMCFileReader() :
46			evt_(0), input_(0)
47			{
48			// Default constructor.
49			if (instance_ == 0) {
50			instance_ = this;
51			} else {
52			edm::LogError("HepMCFileReader") << "Constructing a new instance";
53			}
54			}
55
56
57			//-------------------------------------------------------------------------
58			HepMCFileReader::~HepMCFileReader()
59			{
60			edm::LogInfo("HepMCFileReader") << "Destructing HepMCFileReader";
61
62			instance_=0;
63			delete input_;
64			}
65
66
67			//-------------------------------------------------------------------------
68			void HepMCFileReader::initialize(const string &filename)
69			{
70			if (isInitialized()) {
71			edm::LogError("HepMCFileReader") << "Was already initialized... reinitializing";
72			delete input_;
73			}
74
75			edm::LogInfo("HepMCFileReader") << "Opening file" << filename << "using HepMC::IO_GenEvent";
76			input_ = new HepMC::IO_GenEvent(filename.c_str(), std::ios::in);
77
78			if (rdstate() == std::ios::failbit) {
79			throw cms::Exception("FileNotFound", "HepMCFileReader::initialize()")
80			<< "File " << filename << " was not found.\n";
81			}
82			}
83
84
85			//-------------------------------------------------------------------------
86			int HepMCFileReader::rdstate() const
87			{
88			// work around a HepMC IO_ inheritence shortfall
89
90			HepMC::IO_GenEvent p = dynamic_cast<HepMC::IO_GenEvent>(input_);
91			if (p) return p->rdstate();
92
93			return std::ios::failbit;
94			}
95
96
97			//-------------------------------------------------------------------------
98			bool HepMCFileReader::readCurrentEvent()
99			{
100			evt_ = input_->read_next_event();
101			if (evt_) {
102			edm::LogInfo("HepMCFileReader") << "\| --- Event Nr. " << evt_->event_number()
103			<< " with " << evt_->particles_size() << " particles --- !";
104			ReadStats();
105			// printHepMcEvent();
106			// printEvent();
107			} else {
108			edm::LogInfo("HepMCFileReader") << "Got no event" <<endl;
109			}
110
111			return evt_ != 0;
112			}
113
114
115			//-------------------------------------------------------------------------
116			bool HepMCFileReader::setEvent(int event)
117			{
118			return true;
119			}
120
121
122			//-------------------------------------------------------------------------
123			bool HepMCFileReader::printHepMcEvent() const
124			{
125			if (evt_ != 0) evt_->print();
126			return true;
127			}
128
129
130			//-------------------------------------------------------------------------
131			HepMC::GenEvent * HepMCFileReader::fillCurrentEventData()
132			{
133			readCurrentEvent();
134			return evt_;
135			}
136
137
138			//-------------------------------------------------------------------------
139			// Print out in old CMKIN style for comparisons
140			void HepMCFileReader::printEvent() const {
141			int mo1=0,mo2=0,da1=0,da2=0,status=0,pid=0;
142			if (evt_ != 0) {
143			cout << "---#-------pid--st---Mo1---Mo2---Da1---Da2------px------py------pz-------E-";
144			cout << "------m---------x---------y---------z---------t-";
145			cout << endl;
146			cout.setf(ios::right, ios::adjustfield);
147			for(int n=1; n<=evt_->particles_size(); n++) {
148			HepMC::GenParticle *g = index_to_particle[n];
149			getStatsFromTuple( mo1,mo2,da1,da2,status,pid,n);
150			cout << setw(4) << n
151			<< setw(8) << pid
152			<< setw(5) << status
153			<< setw(6) << mo1
154			<< setw(6) << mo2
155			<< setw(6) << da1
156			<< setw(6) << da2;
157			cout.setf(ios::fixed, ios::floatfield);
158			cout.setf(ios::right, ios::adjustfield);
159			cout << setw(10) << setprecision(2) << g->momentum().x();
160			cout << setw(8) << setprecision(2) << g->momentum().y();
161			cout << setw(10) << setprecision(2) << g->momentum().z();
162			cout << setw(8) << setprecision(2) << g->momentum().t();
163			cout << setw(8) << setprecision(2) << g->generatedMass();
164			// tau=L/(gammabetac)
165			if (g->production_vertex() != 0 && g->end_vertex() != 0 && status == 2) {
166			cout << setw(10) << setprecision(2) <<g->production_vertex()->position().x();
167			cout << setw(10) << setprecision(2) <<g->production_vertex()->position().y();
168			cout << setw(10) << setprecision(2) <<g->production_vertex()->position().z();
169
170			double xm = g->production_vertex()->position().x();
171			double ym = g->production_vertex()->position().y();
172			double zm = g->production_vertex()->position().z();
173			double xd = g->end_vertex()->position().x();
174			double yd = g->end_vertex()->position().y();
175			double zd = g->end_vertex()->position().z();
176			double decl = sqrt((xd-xm)(xd-xm)+(yd-ym)(yd-ym)+(zd-zm)*(zd-zm));
177			double labTime = decl/c_light;
178			// convert lab time to proper time
179			double properTime = labTime/g->momentum().rho()*(g->generatedMass() );
180			// set the proper time in nanoseconds
181			cout << setw(8) << setprecision(2) <<properTime;
182			}
183			else{
184			cout << setw(10) << setprecision(2) << 0.0;
185			cout << setw(10) << setprecision(2) << 0.0;
186			cout << setw(10) << setprecision(2) << 0.0;
187			cout << setw(8) << setprecision(2) << 0.0;
188			}
189			cout << endl;
190			}
191			} else {
192			cout << " HepMCFileReader: No event available !" << endl;
193			}
194			}
195
196
197			//-------------------------------------------------------------------------
198			void HepMCFileReader::ReadStats()
199			{
200			unsigned int particle_counter=0;
201			index_to_particle.reserve(evt_->particles_size()+1);
202			index_to_particle[0] = 0;
203			for (HepMC::GenEvent::vertex_const_iterator v = evt_->vertices_begin();
204			v != evt_->vertices_end(); ++v ) {
205			// making a list of incoming particles of the vertices
206			// so that the mother indices in HEPEVT can be filled properly
207			for (HepMC::GenVertex::particles_in_const_iterator p1 = (*v)->particles_in_const_begin();
208			p1 != (*v)->particles_in_const_end(); ++p1 ) {
209			++particle_counter;
210			//particle_counter can be very large for heavy ions
211			if(particle_counter >= index_to_particle.size() ) {
212			//make it large enough to hold up to this index
213			index_to_particle.resize(particle_counter+1);
214			}
215			index_to_particle[particle_counter] = *p1;
216			particle_to_index[*p1] = particle_counter;
217			}
218			// daughters are entered only if they aren't a mother of
219			// another vertex
220			for (HepMC::GenVertex::particles_out_const_iterator p2 = (*v)->particles_out_const_begin();
221			p2 != (*v)->particles_out_const_end(); ++p2) {
222			if (!(*p2)->end_vertex()) {
223			++particle_counter;
224			//particle_counter can be very large for heavy ions
225			if(particle_counter >= index_to_particle.size() ) {
226			//make it large enough to hold up to this index
227			index_to_particle.resize(particle_counter+1);
228			}
229			index_to_particle[particle_counter] = *p2;
230			particle_to_index[*p2] = particle_counter;
231			}
232			}
233			}
234			}
235
236
237			//-------------------------------------------------------------------------
238			void HepMCFileReader::getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2,
239			int &status, int &pid, int j) const
240			{
241			if (!evt_) {
242			cout << "HepMCFileReader: Got no event :-( Game over already ?" <<endl;
243			} else {
244			status = index_to_particle[j]->status();
245			pid = index_to_particle[j]->pdg_id();
246			if ( index_to_particle[j]->production_vertex() ) {
247
248			//HepLorentzVector p = index_to_particle[j]->
249			//production_vertex()->position();
250
251			int num_mothers = index_to_particle[j]->production_vertex()->
252			particles_in_size();
253			int first_mother = find_in_map( particle_to_index,
254			*(index_to_particle[j]->
255			production_vertex()->
256			particles_in_const_begin()));
257			int last_mother = first_mother + num_mothers - 1;
258			if ( first_mother == 0 ) last_mother = 0;
259			mo1=first_mother;
260			mo2=last_mother;
261			} else {
262			mo1 =0;
263			mo2 =0;
264			}
265			if (!index_to_particle[j]->end_vertex()) {
266			//find # of 1. daughter
267			int first_daughter = find_in_map( particle_to_index,
268			*(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
269			//cout <<"first_daughter "<< first_daughter << "num_daughters " << num_daughters << endl;
270			HepMC::GenVertex::particle_iterator ic;
271			int last_daughter=0;
272			//find # of last daughter
273			for (ic = index_to_particle[j]->end_vertex()->particles_begin(HepMC::children);
274			ic != index_to_particle[j]->end_vertex()->particles_end(HepMC::children); ++ic)
275			last_daughter= find_in_map( particle_to_index,*ic);
276
277			if (first_daughter== 0) last_daughter = 0;
278			da1=first_daughter;
279			da2=last_daughter;
280			} else {
281			da1=0;
282			da2=0;
283			}
284			}
285			}
286
287
288			//-------------------------------------------------------------------------
289			int HepMCFileReader::find_in_map( const std::map<HepMC::GenParticle*,int>& m,
290			HepMC::GenParticle *p) const
291			{
292			std::map<HepMC::GenParticle*,int>::const_iterator iter = m.find(p);
293			return (iter == m.end()) ? 0 : iter->second;
294			}
295