Geant4_10
G4HEPEvtInterface.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4HEPEvtInterface.cc 70215 2013-05-27 07:27:49Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 
32 #include "G4HEPEvtInterface.hh"
33 
34 #include "G4Types.hh"
35 #include "G4SystemOfUnits.hh"
36 
37 #include "G4ios.hh"
38 #include "G4PrimaryVertex.hh"
39 #include "G4PrimaryParticle.hh"
40 #include "G4HEPEvtParticle.hh"
41 #include "G4Event.hh"
42 
44 :vLevel(vl)
45 {
46  inputFile.open((char*)evfile);
47  if (inputFile.is_open()) {
48  fileName = evfile;
49  if(vl>0) G4cout << "G4HEPEvtInterface - " << fileName << " is open." << G4endl;
50  }
51  else {
52  G4Exception("G4HEPEvtInterface::G4HEPEvtInterface","Event0201",FatalException,
53  "G4HEPEvtInterface:: cannot open file.");
54  }
55  G4ThreeVector zero;
56  particle_position = zero;
57  particle_time = 0.0;
58 
59 }
60 /*************************************
61 G4HEPEvtInterface::G4HEPEvtInterface(G4String evfile)
62 {
63  const char* fn = evfile.data();
64  inputFile.open((char*)fn);
65  if (inputFile.is_open()) {
66  fileName = evfile;
67  }
68  else {
69  G4Exception("G4HEPEvtInterface::G4HEPEvtInterface","Event0201",FatalException,
70  "G4HEPEvtInterface:: cannot open file.");
71  }
72  G4ThreeVector zero;
73  particle_position = zero;
74  particle_time = 0.0;
75 }
76 **************************************/
78 {;}
79 
81 {
82  G4int NHEP = 0; // number of entries
83  if (inputFile.is_open()) {
84  inputFile >> NHEP;
85  }
86  else {
87  G4Exception("G4HEPEvtInterface::G4HEPEvtInterface","Event0201",FatalException,
88  "G4HEPEvtInterface:: cannot open file.");
89  }
90  if( inputFile.eof() )
91  {
92  G4Exception("G4HEPEvtInterface::GeneratePrimaryVertex","Event0202",
93  JustWarning,"End-Of-File : HEPEvt input file");
94  return;
95  }
96 
97  if(vLevel > 0)
98  {
99  G4cout << "G4HEPEvtInterface - reading " << NHEP << " HEPEvt particles from "
100  << fileName << "." << G4endl;
101  }
102  for( G4int IHEP=0; IHEP<NHEP; IHEP++ )
103  {
104  G4int ISTHEP; // status code
105  G4int IDHEP; // PDG code
106  G4int JDAHEP1; // first daughter
107  G4int JDAHEP2; // last daughter
108  G4double PHEP1; // px in GeV
109  G4double PHEP2; // py in GeV
110  G4double PHEP3; // pz in GeV
111  G4double PHEP5; // mass in GeV
112 
113  inputFile >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
114  >> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5;
115  if( inputFile.eof() )
116  {
117  G4Exception("G4HEPEvtInterface::GeneratePrimaryVertex","Event0202",
118  FatalException,"Unexpected End-Of-File : HEPEvt input file");
119  }
120  if(vLevel > 1)
121  {
122  G4cout << " " << ISTHEP << " " << IDHEP << " " << JDAHEP1 << " " << JDAHEP2
123  << " " << PHEP1 << " " << PHEP2 << " " << PHEP3 << " " << PHEP5
124  << G4endl;
125  }
126 
127  // create G4PrimaryParticle object
128  G4PrimaryParticle* particle
129  = new G4PrimaryParticle( IDHEP );
130  particle->SetMass( PHEP5*GeV );
131  particle->SetMomentum(PHEP1*GeV, PHEP2*GeV, PHEP3*GeV );
132 
133  // create G4HEPEvtParticle object
134  G4HEPEvtParticle* hepParticle
135  = new G4HEPEvtParticle( particle, ISTHEP, JDAHEP1, JDAHEP2 );
136 
137  // Store
138  HPlist.push_back( hepParticle );
139  }
140 
141  // check if there is at least one particle
142  if( HPlist.size() == 0 ) return;
143 
144  // make connection between daughter particles decayed from
145  // the same mother
146  for( size_t i=0; i<HPlist.size(); i++ )
147  {
148  if( HPlist[i]->GetJDAHEP1() > 0 ) // it has daughters
149  {
150  G4int jda1 = HPlist[i]->GetJDAHEP1()-1; // FORTRAN index starts from 1
151  G4int jda2 = HPlist[i]->GetJDAHEP2()-1; // but C++ starts from 0.
152  G4PrimaryParticle* mother = HPlist[i]->GetTheParticle();
153  for( G4int j=jda1; j<=jda2; j++ )
154  {
155  G4PrimaryParticle* daughter = HPlist[j]->GetTheParticle();
156  if(HPlist[j]->GetISTHEP()>0)
157  {
158  mother->SetDaughter( daughter );
159  HPlist[j]->Done();
160  }
161  }
162  }
163  }
164 
165  // create G4PrimaryVertex object
167 
168  // put initial particles to the vertex
169  for( size_t ii=0; ii<HPlist.size(); ii++ )
170  {
171  if( HPlist[ii]->GetISTHEP() > 0 ) // ISTHEP of daughters had been
172  // set to negative
173  {
174  G4PrimaryParticle* initialParticle = HPlist[ii]->GetTheParticle();
175  vertex->SetPrimary( initialParticle );
176  }
177  }
178 
179  // clear G4HEPEvtParticles
180  //HPlist.clearAndDestroy();
181  for(size_t iii=0;iii<HPlist.size();iii++)
182  { delete HPlist[iii]; }
183  HPlist.clear();
184 
185  // Put the vertex to G4Event object
186  evt->AddPrimaryVertex( vertex );
187 }
188 
void AddPrimaryVertex(G4PrimaryVertex *aPrimaryVertex)
Definition: G4Event.hh:143
void GeneratePrimaryVertex(G4Event *evt)
int G4int
Definition: G4Types.hh:78
G4ThreeVector particle_position
G4HEPEvtInterface(const char *evfile, G4int vl=0)
G4GLOB_DLL std::ostream G4cout
void SetMass(G4double mas)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void SetDaughter(G4PrimaryParticle *np)
void SetPrimary(G4PrimaryParticle *pp)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetMomentum(G4double px, G4double py, G4double pz)