Belle II Software  release-05-01-25
TouschekReaderTURTLE.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/touschek/TouschekReaderTURTLE.h>
12 
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Unit.h>
15 #include <mdst/dataobjects/MCParticle.h>
16 
17 #include <boost/foreach.hpp>
18 #include <boost/lexical_cast.hpp>
19 #include <boost/algorithm/string.hpp>
20 #include <boost/tokenizer.hpp>
21 
22 using namespace std;
23 using namespace Belle2;
24 
25 
26 TouschekReaderTURTLE::TouschekReaderTURTLE(TGeoHMatrix* transMatrix, int pdg): m_transMatrix(transMatrix), m_pdg(pdg),
27  m_lineNum(0)
28 {
29 }
30 
31 
33 {
34  if (m_input) m_input.close();
35 }
36 
37 
38 void TouschekReaderTURTLE::open(const string& filename)
39 {
40  m_lineNum = 0;
41  m_input.open(filename.c_str());
42  if (!m_input) throw(TouschekCouldNotOpenFileError() << filename);
43 }
44 
45 
47 {
48  if (m_transMatrix == NULL) {
49  B2ERROR("The transformation matrix is NULL !");
50  }
51 
52  int numParticles = 0;
53  string currLine;
54  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
55  boost::char_separator<char> sep(",; \t");
56  double fields[7];
57  double particlePosTouschek[3] = {0.0, 0.0, 0.0};
58  double particlePosGeant4[3] = {0.0, 0.0, 0.0};
59  double particleMomTouschek[3] = {0.0, 0.0, 0.0};
60  double particleMomGeant4[3] = {0.0, 0.0, 0.0};
61 
62  do {
63  getline(m_input, currLine);
64  m_lineNum++;
65 
66  boost::trim(currLine);
67  if (currLine == "") continue;
68 
69  tokenizer tokens(currLine, sep);
70 
71  int index = 0;
72  BOOST_FOREACH(const string & tok, tokens) {
73  try {
74  if (index >= 7) {
75  B2WARNING("This Touschek file has more than 7 fields ! Only the first 7 fields were read.");
76  break;
77  }
78  fields[index] = boost::lexical_cast<double>(tok);
79  index++;
80  } catch (boost::bad_lexical_cast&) {
81  throw (TouschekConvertFieldError() << m_lineNum << index << tok);
82  }
83  }
84 
85  //Convert the position of the particle from local Touschek plane space to global geant4 space.
86  //Flip the sign for the y and z component to go from the accelerator to the detector coordinate system
87  particlePosTouschek[0] = fields[2] * Unit::m;
88  particlePosTouschek[1] = -fields[3] * Unit::m;
89  if (m_pdg == -11) particlePosTouschek[2] = 400.0; //forLER
90  else if (m_pdg == 11) particlePosTouschek[2] = -400.0; //forHER
91  else { B2WARNING("m_pdg is not 11/-11"); break;}
92 
93  m_transMatrix->LocalToMaster(particlePosTouschek, particlePosGeant4);
94 
95  //Convert the momentum of the particle from local Touschek plane space to global geant4 space.
96  //Flip the sign for the y and z component to go from the accelerator to the detector coordinate system
97  particleMomTouschek[0] = fields[4];
98  particleMomTouschek[1] = -fields[5];
99  if (m_pdg == -11) particleMomTouschek[2] = -fields[6]; //forLER
100  else if (m_pdg == 11) particleMomTouschek[2] = fields[6]; //forHER
101  else { B2WARNING("m_pdg is not 11/-11"); break;}
102 
103  m_transMatrix->LocalToMasterVect(particleMomTouschek, particleMomGeant4);
104 
105  double totalMom2 = particleMomGeant4[0] * particleMomGeant4[0];
106  totalMom2 += particleMomGeant4[1] * particleMomGeant4[1];
107  totalMom2 += particleMomGeant4[2] * particleMomGeant4[2];
108 
109  //Add particles to MCParticle collection
110  MCParticleGraph::GraphParticle& particle = graph.addParticle();
111  particle.setStatus(MCParticle::c_PrimaryParticle);
112  particle.setPDG(m_pdg);
113  particle.setMassFromPDG();
114  particle.setMomentum(TVector3(particleMomGeant4));
115  particle.setProductionVertex(TVector3(particlePosGeant4));
116  particle.setEnergy(sqrt(totalMom2 + particle.getMass()*particle.getMass()));
117  particle.setProductionTime(0.0);
118  particle.setValidVertex(true);
119 
120  numParticles++;
121  } while ((!m_input.eof()) && ((number > -1) && (numParticles < number)));
122 
123  return numParticles;
124 }
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::TouschekReaderTURTLE::open
void open(const std::string &filename)
Opens an ascii file and prepares it for reading.
Definition: TouschekReaderTURTLE.cc:38
Belle2::TouschekReaderTURTLE::m_lineNum
int m_lineNum
The line number in the ascii file of the last particle which was read.
Definition: TouschekReaderTURTLE.h:89
Belle2::TouschekReaderTURTLE::m_pdg
int m_pdg
The pdg value of the type of particle that is read (e.g.
Definition: TouschekReaderTURTLE.h:87
Belle2::TouschekReaderTURTLE::m_transMatrix
TGeoHMatrix * m_transMatrix
Transformation matrix from local Touschek to global geant4 space.
Definition: TouschekReaderTURTLE.h:86
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::Unit::m
static const double m
[meters]
Definition: Unit.h:79
Belle2::TouschekReaderTURTLE::getParticles
int getParticles(int number, MCParticleGraph &graph)
Reads the specified number of particles from the file and stores the result in the given MCParticle g...
Definition: TouschekReaderTURTLE.cc:46
Belle2::TouschekReaderTURTLE::~TouschekReaderTURTLE
~TouschekReaderTURTLE()
Destructor.
Definition: TouschekReaderTURTLE.cc:32
Belle2::TouschekReaderTURTLE::m_input
std::ifstream m_input
The input stream of the ascii file.
Definition: TouschekReaderTURTLE.h:88
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58