ReaderSAD Class Reference

Class to read files that have been created by SAD and store their content in a MCParticle graph. More...

#include <ReaderSAD.h>

Classes
struct	bendingElement
	Sensitive Element More...
struct	straightElement
	Calculates the transformation matrix from local SAD to global geant4 space. More...

Public Types
enum	AcceleratorRings {   c_HER = 0 ,   c_LER = 1 }
	The both accelerator rings. More...

Public Member Functions
	BELLE2_DEFINE_EXCEPTION (SADCouldNotOpenFileError, "Could not open file %1% !")
	Exception is thrown if the SAD file could not be opened.
	BELLE2_DEFINE_EXCEPTION (SADEndOfFile, "End of the SAD file.")
	Exception is thrown if the end of the SAD file has been reached.
	ReaderSAD ()
	Constructor of the ReaderSAD class.
	~ReaderSAD ()
	Destructor.
void	initialize (TGeoHMatrix *transMatrix, double sRange, AcceleratorRings accRing, double readoutTime)
	Initializes the reader, sets the particle parameters and calculates important values.
void	open (const std::string &filename)
	Opens a root file and prepares it for reading.
void	setMomentumRes (double pxRes, double pyRes)
	Sets the resolution of the momentum for the real particles.
double	getSADParticle (MCParticleGraph &graph)
	Reads one SAD particle from the file and creates one event per SAD particle.
bool	getRealParticle (MCParticleGraph &graph)
	Reads one SAD particle from the file, calculates the number of real particles which are represented by the SAD particle and creates one event per real particle.
void	addAllSADParticles (MCParticleGraph &graph)
	Reads all SAD particles from the file into the MCParticles collection which are inside the specified s range.
TGeoHMatrix	SADtoGeant (ReaderSAD::AcceleratorRings accRing, double s)
	Transformation matrix.

Protected Attributes
TFile *	m_file
	The input root file.
TTree *	m_tree
	The input root tree.
TGeoHMatrix *	m_transMatrix
	Transformation matrix from local SAD to global geant4 space.
double	m_sRange
	The +- range for the s value for which particles are loaded.
AcceleratorRings	m_accRing
	The accelerator ring from which the particles originate.
double	m_pxRes
	The resolution for the x momentum component of the SAD real particle.
double	m_pyRes
	The resolution for the y momentum component of the SAD real particle.
double	m_SADToRealFactor
	The factor to calculate the number of real particles from a SAD particle.
double	m_readoutTime
	The readout time.
unsigned int	m_realPartNum
	The current number of the created real particles.
unsigned int	m_realPartEntry
	The current number of the created real particles.
int	m_readEntry
	The current number of the SAD entry that is read.
double	m_lostX
	x at lost position [m].
double	m_lostY
	y at lost position [m].
double	m_lostS
	lost position [m] along ring.
double	m_lostPx
	x momentum at lost position [m].
double	m_lostPy
	y momentum at lost position [m].
double	m_lostRate
	loss rate [Hz]>
double	m_lostE
	energy at lost position [m].
double	m_inputSAD_ssraw
	scattered position [m]
double	m_inputSAD_sraw
	lost position [m
double	m_inputSAD_ss
	scattered position (|s|<Ltot/2) [m]
double	m_inputSAD_s
	lost position (|s|<Ltot/2) [m]
double	m_inputSAD_Lss
	length of element in which scattered [m]
int	m_inputSAD_nturn
	number of turns from scattered to lost
double	m_inputSAD_x
	x at lost position [m].
double	m_inputSAD_y
	y at lost position [m].
double	m_inputSAD_px
	x momentum at lost position [m].
double	m_inputSAD_py
	y momentum at lost position [m].
double	m_inputSAD_xraw
	x at lost position [m] before matching onto beam pipe inner surface
double	m_inputSAD_yraw
	y at lost position [m] before matching onto beam pipe inner surface
double	m_inputSAD_r
	sqrt(x*x+y*y) [m]
double	m_inputSAD_rr
	sqrt(x*x+y*y) [m] before matching onto beam pipe inner surface
double	m_inputSAD_dp_over_p0
	dp_over_p0
double	m_inputSAD_E
	energy at lost position [m].
double	m_inputSAD_rate
	loss rate [Hz]
double	m_inputSAD_watt
	loss wattage [W]

Private Member Functions
void	convertParamsToSADUnits ()
	Convert the parameters from the SAD units to the basf2 units.
void	addParticleToMCParticles (MCParticleGraph &graph, bool gaussSmearing=false)
	Adds the current particle described by the member variables to the MCParticles collection.
int	calculateRealParticleNumber (double rate)
	Calculates the number of real particles for a SAD particle.

Detailed Description

Class to read files that have been created by SAD and store their content in a MCParticle graph.

The input data is stored in a root file and contains the particles together with their lost rate. The reader reads one particle from the file, calculates the number of 'real' particles and creates a new event for each of them.

Definition at line 35 of file ReaderSAD.h.

Member Enumeration Documentation

AcceleratorRings

enum AcceleratorRings

The both accelerator rings.

Enumerator
c_HER	High Energy Ring (electrons)
c_LER	Low Energy Ring (positrons)

Definition at line 46 of file ReaderSAD.h.

                          {
      c_HER = 0, 
      c_LER = 1  
    };

Constructor & Destructor Documentation

ReaderSAD()

ReaderSAD

(

)

Constructor of the ReaderSAD class.

Definition at line 28 of file ReaderSAD.cc.

                    : m_file(NULL), m_tree(NULL), m_transMatrix(NULL),
  m_sRange(3000.0), m_accRing(ReaderSAD::c_LER), m_pxRes(0.01), m_pyRes(0.01),
  m_SADToRealFactor(5.76e6), m_readoutTime(20.0), m_realPartNum(0),
  m_realPartEntry(0), m_readEntry(0)
{
  m_lostX = 0.0;
  m_lostY = 0.0;
  m_lostS = 0.0;
  m_lostPx = 0.0;
  m_lostPy = 0.0;
  m_lostRate = 0.0;
  m_lostE = 0.0;
 
  //Start my addition
  m_inputSAD_ssraw =  0;
  m_inputSAD_sraw = 0;
  m_inputSAD_ss = 0;
  m_inputSAD_s = 0;
  m_inputSAD_Lss = 0;
  m_inputSAD_nturn = 0;
  m_inputSAD_x = 0;
  m_inputSAD_y = 0;
  m_inputSAD_px = 0;
  m_inputSAD_py = 0;
  m_inputSAD_xraw = 0;
  m_inputSAD_yraw = 0;
  m_inputSAD_r = 0;
  m_inputSAD_rr = 0;
  m_inputSAD_dp_over_p0 = 0;
  m_inputSAD_E = 0;
  m_inputSAD_rate = 0;
  m_inputSAD_watt = 0;
  //End my addition
}

~ReaderSAD()

~ReaderSAD

(

)

Destructor.

Definition at line 64 of file ReaderSAD.cc.

{
  // comment out the following line for the march 16 build (TF)
  //  if (m_file != NULL) m_file->Close();
}

Member Function Documentation

addAllSADParticles()

void addAllSADParticles

(

MCParticleGraph &

graph

)

Reads all SAD particles from the file into the MCParticles collection which are inside the specified s range.

Parameters

graph

Reference to the graph which should be filled with the information from the SAD file.

Definition at line 212 of file ReaderSAD.cc.

{
  if (m_tree == NULL) {
    B2ERROR("The SAD tree doesn't exist !");
    return ;
  }
 
  int nPart = m_tree->GetEntries();
 
  for (int iPart = 0; iPart < nPart; ++iPart) {
    m_tree->GetEntry(iPart);
    convertParamsToSADUnits();
    if (fabs(m_lostS) <= m_sRange) addParticleToMCParticles(graph);
 
    //std::cout << " in sad " <<  m_inputSAD_sraw << std::endl;
 
  }
}

addParticleToMCParticles()

void addParticleToMCParticles ( MCParticleGraph &  graph, bool  gaussSmearing = false  )

private

Adds the current particle described by the member variables to the MCParticles collection.

Parameters

graph	Reference to the graph which should be filled with the information from the SAD file.
gaussSmearing	If set to true the particle momentum is smeared using a Gaussian.

Definition at line 266 of file ReaderSAD.cc.

{
  double particlePosSAD[3] = {0.0, 0.0, 0.0};
  double particlePosSADfar[3] = {0.0, 0.0, 0.0};
  double particlePosGeant4[3] = {0.0, 0.0, 0.0};
  double particleMomSAD[3] = {0.0, 0.0, 0.0};
  double particleMomGeant4[3] = {0.0, 0.0, 0.0};
 
  //Add particle to MCParticle collection
  MCParticleGraph::GraphParticle& particle = graph.addParticle();
  particle.setStatus(MCParticle::c_PrimaryParticle);
 
  switch (m_accRing) {
    case c_HER: particle.setPDG(11); //electrons
      break;
    case c_LER: particle.setPDG(-11); //positrons
      break;
  }
 
  particle.setMassFromPDG();
 
  //Convert the position of the particle from local SAD space to global geant4 space.
  //Flip the sign for the y and z component to go from the accelerator to the detector coordinate system.
  particlePosSAD[0] = m_lostX;
  particlePosSAD[1] = -m_lostY;
  particlePosSAD[2] = -m_lostS;
 
  particlePosSADfar[0] = m_lostX;
  particlePosSADfar[1] = -m_lostY;
  particlePosSADfar[2] = 0;
 
 
  static GearDir content = Gearbox::getInstance().getDetectorComponent("FarBeamLine");
  if (!content)
    B2FATAL("You need FarBeamLine.xml to run SADInput module. Please include FarBeamLine.xml in Belle2.xml. You also need to change 'length' in Belle2.xml to be 40m.");
 
  TGeoHMatrix* m_transMatrix2 = new TGeoHMatrix(SADtoGeant(m_accRing, m_lostS)); //overwrite m_transMatrix given by initialize()
 
  if (abs(m_lostS) <= 400.) { //4m
    m_transMatrix->LocalToMaster(particlePosSAD, particlePosGeant4);
  } else {
    m_transMatrix2->LocalToMaster(particlePosSADfar, particlePosGeant4);
  }
 
  //Convert the momentum of the particle from local SAD space to global geant4 space.
  //Flip the sign for the y and z component to go from the accelerator to the detector coordinate system.
  //Calculate the missing pz by using the energy of the particle at the position where it has been lost.
  //double totalMomSqr = (m_lostE * m_lostE - (particle.getMass() * particle.getMass()));
  double totalMomSqr = m_lostE * m_lostE; // SAD output "E" is "P", in fact.
 
  if (gaussSmearing) {
    particleMomSAD[0] = gRandom->Gaus(m_lostPx, m_pxRes * m_lostPx); //1% px resolution
    particleMomSAD[1] = -1.0 * gRandom->Gaus(m_lostPy, m_pyRes * m_lostPy);
  } else {
    particleMomSAD[0] = m_lostPx;
    particleMomSAD[1] = -m_lostPy;
  }
 
  double zMom = sqrt(totalMomSqr - (particleMomSAD[0] *  particleMomSAD[0]) - (particleMomSAD[1] *  particleMomSAD[1]));
 
  switch (m_accRing) {
    case c_HER: particleMomSAD[2] = zMom;
      break;
    case c_LER: particleMomSAD[2] = -zMom;
      break;
  }
 
  int ring = 0;
 
  switch (m_accRing) {
    case c_HER: ring = 1;
      break;
    case c_LER: ring = 2;
      break;
  }
  /*
  int ler_section[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
  int ler_inf_section[12] = {0., 254.74, 498.59, 754.14, 1009.66, 1253.52, 1488.7, 1764.1, 2007.05, 2262.14, 2517.2, 2760.15};
  int her_inf_section[12] = {254.32, 494.31, 750.25, 1009.24, 1249.24, 1488.43, 1763.68, 2002.77, 2258.25, 2516.78, 2755.87, 3011.33};
  int her_section[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
  int ler_sup_section[12] = {0., 254.61, 496.52, 727.23, 1007.7, 1249.3, 1479.82, 1761.31, 2003.88, 2239.89, 2516.31, 2758.89};
  int her_sup_section[12] = {254.3, 495.32, 726.92, 1007.39, 1248.41, 1479.51, 1761, 2002.99, 2239.58, 2516, 2758, 3011.87};
 
  int section = -1;
 
  if (ring == 1) {
    for (int i = 0; i < 12; i++) {
      if (her_inf_section[i] <= (m_inputSAD_ssraw + 1500.) && (m_inputSAD_ssraw + 1500.) <= her_sup_section[i])
  section = her_section[i] - 1;
    }
  } else if  (ring == 2) {
    for (int i = 0; i < 12; i++) {
      if (ler_inf_section[i] <= (m_inputSAD_ssraw + 1500.) && (m_inputSAD_ssraw + 1500.) <= ler_sup_section[i])
  section = ler_section[i] - 1;
    }
  }
  */
  //each rings have 12 section of ~250m
  //the 1st section D01, the second section is D12, followed by D11, D10 .... for both rings
  const int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
  //int ring_section = section_ordering[(int)((m_inputSAD_ssraw + 1500.) / 12.)];
  double ssraw = 0.;
  if (ring == 1) {
    if (m_inputSAD_ssraw >= 0) ssraw = m_inputSAD_ssraw / 100.;
    else  ssraw = 3000. + m_inputSAD_ssraw / 100.;
  } else if (ring == 2) {
    //if (m_inputSAD_ssraw >= 0) ssraw = 3000. - m_inputSAD_ssraw / 100.;
    //else if (m_inputSAD_ssraw < 0) ssraw = -m_inputSAD_ssraw / 100.;
    if (m_inputSAD_ssraw >= 0) ssraw = m_inputSAD_ssraw / 100.;
    else  ssraw = 3000. + m_inputSAD_ssraw / 100.;
  }
  int ring_section = section_ordering[(int)((ssraw) / 250.)];
 
  if (abs(m_lostS) <= 400.) {
    m_transMatrix->LocalToMasterVect(particleMomSAD, particleMomGeant4);
  } else {
    m_transMatrix2->LocalToMasterVect(particleMomSAD, particleMomGeant4);
  }
 
  //Set missing particle information
  particle.setMomentum(ROOT::Math::XYZVector(particleMomGeant4[0], particleMomGeant4[1], particleMomGeant4[2]));
  particle.setProductionVertex(ROOT::Math::XYZVector(particlePosGeant4[0], particlePosGeant4[1], particlePosGeant4[2]));
  particle.setProductionTime(0.0);
  particle.setEnergy(sqrt(m_lostE * m_lostE + particle.getMass()*particle.getMass()));
  particle.setValidVertex(true);
 
  //Start my addition
  StoreArray<SADMetaHit> SADMetaHits;
  SADMetaHits.appendNew(SADMetaHit(m_inputSAD_ssraw, m_inputSAD_sraw, m_inputSAD_ss, m_lostS,
                                   m_inputSAD_Lss, m_inputSAD_nturn,
                                   m_lostX, m_lostY, m_lostPx, m_lostPy, m_inputSAD_xraw, m_inputSAD_yraw,
                                   m_inputSAD_r, m_inputSAD_rr, m_inputSAD_dp_over_p0, m_lostE, m_lostRate,
                                   m_inputSAD_watt, ring, ring_section));
  //End my addition
}

calculateRealParticleNumber()

int calculateRealParticleNumber ( double  rate )

private

Calculates the number of real particles for a SAD particle.

Parameters

rate	The loss rate of the SAD particle.

Returns

The number of real particles for the given loss rate.

Definition at line 403 of file ReaderSAD.cc.

{
  double numPart = rate * m_readoutTime / Unit::s; //readoutTime [ns]
 
  int numPart_int = static_cast<int>(floor(numPart));
  double numPart_dec = numPart  - numPart_int;
 
  double rnd = gRandom->Uniform(); //returns a random number in the interval ]0, 1]
  if (rnd < numPart_dec) numPart_int += 1;
 
  return numPart_int;
}

convertParamsToSADUnits()

void convertParamsToSADUnits ( )

private

Convert the parameters from the SAD units to the basf2 units.

Definition at line 236 of file ReaderSAD.cc.

{
  m_lostX = m_lostX * Unit::m;
  m_lostY = m_lostY * Unit::m;
  m_lostS = m_lostS * Unit::m;
  m_lostPx = m_lostPx * Unit::GeV;
  m_lostPy = m_lostPy * Unit::GeV;
 
  //Start my addition
  m_inputSAD_ssraw =  m_inputSAD_ssraw * Unit::m;
  m_inputSAD_sraw = m_inputSAD_sraw * Unit::m;
  m_inputSAD_ss = m_inputSAD_ss * Unit::m;
  m_inputSAD_s = m_inputSAD_s * Unit::m;
  //m_inputSAD_Lss
  //m_inputSAD_nturn
  m_inputSAD_x = m_inputSAD_x * Unit::m;
  m_inputSAD_y = m_inputSAD_y * Unit::m;
  m_inputSAD_px = m_inputSAD_px * Unit::GeV;
  m_inputSAD_py =  m_inputSAD_py * Unit::GeV;
  m_inputSAD_xraw = m_inputSAD_xraw * Unit::m;
  m_inputSAD_yraw = m_inputSAD_yraw * Unit::m;
  //m_inputSAD_r
  //m_inputSAD_rr
  //m_inputSAD_dp_over_p0
  //m_inputSAD_E
  //m_inputSAD_rate
  //End my addition
}

getRealParticle()

bool getRealParticle

(

MCParticleGraph &

graph

)

Reads one SAD particle from the file, calculates the number of real particles which are represented by the SAD particle and creates one event per real particle.

Parameters

graph

Reference to the graph which should be filled with the information from the SAD file.

Returns

True if the particle could be read.

Definition at line 165 of file ReaderSAD.cc.

{
  if (m_tree == NULL) {
    B2ERROR("The SAD tree doesn't exist !");
    return false;
  }
 
  //Check for end of file
  if ((m_readEntry >= m_tree->GetEntries()) || (m_tree->GetEntries() == 0)) throw SADEndOfFile();
 
  //Check if the number of the real particles is reached
  if (m_realPartEntry >= m_realPartNum) {
    m_realPartEntry = 0;
    m_realPartNum   = 0;
  }
 
  //Check if a new SAD particle has to be read from the file
  if (m_realPartNum == 0) {
    //Read only SAD particles which are inside the chosen sRange
    do {
      m_readEntry++;
      if (m_readEntry >= m_tree->GetEntries()) throw SADEndOfFile();
 
      m_tree->GetEntry(m_readEntry);
      convertParamsToSADUnits();
 
      B2DEBUG(10, "> Read particle " << m_readEntry + 1 << "/" << m_tree->GetEntries() << " with s = " << m_lostS << " cm" <<
              " and rate = " << m_lostRate << " Hz");
    } while ((fabs(m_lostS) > m_sRange) && (m_readEntry < m_tree->GetEntries()));
 
    if (fabs(m_lostS) <= m_sRange)
      m_realPartNum = calculateRealParticleNumber(m_lostRate);
  }
 
  //Create a new real particle from the SAD particle
  if ((fabs(m_lostS) <= m_sRange) && (m_realPartNum > 0)) {
    addParticleToMCParticles(graph);
    B2DEBUG(10, "* Created real particle " << m_realPartEntry + 1 << "/" << m_realPartNum << " for SAD particle " << m_readEntry + 1 <<
            "/" << m_tree->GetEntries());
  }
 
  m_realPartEntry++;
 
  return true;
}

getSADParticle()

double getSADParticle

(

MCParticleGraph &

graph

)

Reads one SAD particle from the file and creates one event per SAD particle.

The loss rate of the SAD particle is stored in the weight attribute of the event meta info.

Parameters

graph

Reference to the graph which should be filled with the information from the SAD file.

Returns

The loss rate of the SAD particle which was read. Returns -1 if an error occurred.

Definition at line 122 of file ReaderSAD.cc.

{
  if (m_tree == NULL) {
    B2ERROR("The SAD tree doesn't exist !");
    return -1;
  }
 
  do {
    m_readEntry++;
 
    //Check for end of file
    if (m_readEntry >= m_tree->GetEntries()) throw SADEndOfFile();
 
    //Load the SAD particle
    m_tree->GetEntry(m_readEntry);
    convertParamsToSADUnits();
    B2DEBUG(10, "> Read particle " << m_readEntry + 1 << "/" << m_tree->GetEntries() << " with s = " << m_lostS << " cm" <<
            " and rate = " << m_lostRate << " Hz");
 
    //printf("Read particle %d / %d with s= %f [m]\n", m_readEntry + 1 , (int)m_tree->GetEntries(), m_lostS / 100.);
 
    double zMom2 = m_lostE * m_lostE - m_lostPx * m_lostPx - m_lostPy * m_lostPy ;
    if (zMom2 < 0) printf("zMom2= %f is negative. Skipped!\n", zMom2);
 
  } while ((fabs(m_lostS) > m_sRange) or (m_lostE * m_lostE - m_lostPx * m_lostPx - m_lostPy * m_lostPy < 0));
 
//  do {
//    //Check for end of file
//    m_readEntry++;
//    if (m_readEntry >= m_tree->GetEntries()) throw SADEndOfFile();
//
//    //Load the SAD particle
//    m_tree->GetEntry(m_readEntry);
//    convertParamsToSADUnits();
//
//    B2DEBUG(10, "> Read particle " << m_readEntry + 1 << "/" << m_tree->GetEntries() << " with s = " << m_lostS << " cm" << " and rate = " << m_lostRate << " Hz" )
//  } while (fabs(m_lostS) > m_sRange);
 
  addParticleToMCParticles(graph);
  return m_lostRate;
}

initialize()

void initialize	(	TGeoHMatrix *	transMatrix,
		double	sRange,
		ReaderSAD::AcceleratorRings	accRing,
		double	readoutTime
	)

Initializes the reader, sets the particle parameters and calculates important values.

Parameters

transMatrix	Pointer to the matrix which transforms the particles from the local SAD to the global geant4 coordinate system.
sRange	The +- range for the s value for which particles are loaded.
accRing	The accelerator ring from which the particles originate.
readoutTime	The readout time of the detector in [ns].

Definition at line 71 of file ReaderSAD.cc.

{
  m_transMatrix = transMatrix;
  m_sRange = sRange;
  m_accRing = accRing;
  m_readoutTime = readoutTime;
}

open()

void open

(

const std::string &

filename

)

Opens a root file and prepares it for reading.

Parameters

filename

The filename of the SAD root file which should be read.

Definition at line 80 of file ReaderSAD.cc.

{
  if (m_file != NULL) {
    m_file->Close();
    delete m_file;
  }
 
  m_file = new TFile(filename.c_str(), "READ");
  if (m_file == NULL) throw(SADCouldNotOpenFileError() << filename);
 
  m_file->cd("");
  m_tree = dynamic_cast<TTree*>(m_file->Get("tp"));
  if (m_tree == NULL) throw(SADCouldNotOpenFileError() << filename);
 
  m_tree->SetBranchAddress("x", &m_lostX);
  m_tree->SetBranchAddress("y", &m_lostY);
  m_tree->SetBranchAddress("s", &m_lostS);
  m_tree->SetBranchAddress("px", &m_lostPx);
  m_tree->SetBranchAddress("py", &m_lostPy);
  m_tree->SetBranchAddress("rate", &m_lostRate);
  m_tree->SetBranchAddress("E", &m_lostE);
 
  m_tree->SetBranchAddress("ssraw", &m_inputSAD_ssraw);
  m_tree->SetBranchAddress("sraw", &m_inputSAD_sraw);
  m_tree->SetBranchAddress("ss", &m_inputSAD_ss);
  m_tree->SetBranchAddress("Lss", &m_inputSAD_Lss);
  m_tree->SetBranchAddress("nturn", &m_inputSAD_nturn);
  m_tree->SetBranchAddress("xraw", &m_inputSAD_xraw);
  m_tree->SetBranchAddress("yraw", &m_inputSAD_yraw);
  m_tree->SetBranchAddress("r", &m_inputSAD_r);
  m_tree->SetBranchAddress("rr", &m_inputSAD_rr);
  m_tree->SetBranchAddress("dp_over_p0", &m_inputSAD_dp_over_p0);
  m_tree->SetBranchAddress("watt", &m_inputSAD_watt);
 
  m_readEntry = -1;
 
  StoreArray<SADMetaHit> hit;
  hit.registerInDataStore();
 
}

SADtoGeant()

TGeoHMatrix SADtoGeant	(	ReaderSAD::AcceleratorRings	accRing,
		double	s
	)

Transformation matrix.

Parameters

accRing	The accelerator ring from which the particles originate.
s	s value.

Definition at line 416 of file ReaderSAD.cc.

{
  // 0<sraw<3016.3145 m
  // -1500<s<1500 m
 
  //static double max_s_her = 3016.3145 * Unit::m;
  //static double max_s_ler = 3016.3026 * Unit::m;
 
  //get parameters from .xml file
  static GearDir content = Gearbox::getInstance().getDetectorComponent("FarBeamLine");
 
  map<string, straightElement> straights;
  map<string, bendingElement> bendings;
  for (const GearDir& element : content.getNodes("Straight")) {
 
    string name = element.getString("@name");
    string type = element.getString("@type");
 
    if (type != "pipe") continue;
 
    straightElement straight;
 
    for (const GearDir& slot : element.getNodes("sec")) {
      string nameSec = slot.getString("@name");
      if (nameSec.find("X0")  != std::string::npos) {straight.x0  = slot.getLength();}
      if (nameSec.find("Z0")  != std::string::npos) {straight.z0  = slot.getLength();}
      if (nameSec.find("L")   != std::string::npos) {straight.l   = slot.getLength();}
      if (nameSec.find("PHI") != std::string::npos) {straight.phi = slot.getAngle();}
    }
 
    //straight.x0 = element.getLength("X0");
    //straight.z0 = element.getLength("Z0");
    //straight.l = element.getLength("L");
    //straight.phi = element.getLength("PHI");
 
    straights[name] = straight;
  }
 
  string str_checklist[] = {"LHR1", "LHR2", "LLR1", "LLR2", "LLR3", "LLR4", "LLR5", "LLR6", "LHL1", "LHL2", "LLL1", "LLL2", "LLL3", "LLL4", "LLL5"};
  for (const string& str : str_checklist) {
    if (straights.count(str) == 0)
      B2FATAL("You need FarBeamLine.xml to run SADInput module. Please include FarBeamLine.xml in Belle2.xml. You also need to change 'length' in Belle2.xml to be 40m.");
  }
 
  for (const GearDir& element : content.getNodes("Bending")) {
 
    string name = element.getString("@name");
    string type = element.getString("@type");
 
    if (type != "pipe") continue;
 
    bendingElement bending;
 
    for (const GearDir& slot : element.getNodes("sec")) {
      string nameSec = slot.getString("@name");
      if (nameSec.find("RT")   != std::string::npos) {bending.rt   = slot.getLength();}
      if (nameSec.find("X0")   != std::string::npos) {bending.x0   = slot.getLength();}
      if (nameSec.find("Z0")   != std::string::npos) {bending.z0   = slot.getLength();}
      if (nameSec.find("SPHI") != std::string::npos) {bending.sphi = slot.getAngle();}
      if (nameSec.find("DPHI") != std::string::npos) {bending.dphi = slot.getAngle();}
    }
 
    //bending.rt = element.getLength("RT");
    //bending.x0 = element.getLength("X0");
    //bending.z0 = element.getLength("Z0");
    //bending.sphi = element.getLength("SPHI");
    //bending.dphi = element.getLength("DPHI");
 
    bendings[name] = bending;
  }
 
  string bend_checklist[] = {"BLC2RE", "BC1RP", "BLCWRP", "BLC1RP", "BLC2RP", "BLC2RP.2", "BLY2RP.2", "BLC1LE", "BC1LP", "BLC1LP1", "BLC1LP2", "BLC2LP"};
  for (const string& bnd : bend_checklist) {
    if (bendings.count(bnd) == 0)
      B2FATAL("You need FarBeamLine.xml to run SADInput module. Please include FarBeamLine.xml in Belle2.xml. You also need to change 'length' in Belle2.xml to be 40m.");
  }
 
  static double her_breakpoints[6];
  static double ler_breakpoints[21];
 
  // positive s
  her_breakpoints[0] = straights["LHL1"].l;
  her_breakpoints[1] = her_breakpoints[0] + bendings["BLC1LE"].rt * bendings["BLC1LE"].dphi;
  her_breakpoints[2] = her_breakpoints[1] + straights["LHL2"].l;
 
  // negative s
  her_breakpoints[3] = -straights["LHR1"].l;
  her_breakpoints[4] = her_breakpoints[3] - bendings["BLC2RE"].rt * bendings["BLC2RE"].dphi;
  her_breakpoints[5] = her_breakpoints[4] - straights["LHR2"].l;
 
  // positive s
  ler_breakpoints[0] = straights["LLL1"].l;
  ler_breakpoints[1] = ler_breakpoints[0] + bendings["BC1LP"].rt * bendings["BC1LP"].dphi;
  ler_breakpoints[2] = ler_breakpoints[1] + straights["LLL2"].l;
  ler_breakpoints[3] = ler_breakpoints[2] + bendings["BLC1LP1"].rt * bendings["BLC1LP1"].dphi;
  ler_breakpoints[4] = ler_breakpoints[3] + straights["LLL3"].l;
  ler_breakpoints[5] = ler_breakpoints[4] + bendings["BLC1LP2"].rt * bendings["BLC1LP2"].dphi;
  ler_breakpoints[6] = ler_breakpoints[5] + straights["LLL4"].l;
  ler_breakpoints[7] = ler_breakpoints[6] + bendings["BLC2LP"].rt * bendings["BLC2LP"].dphi;
  ler_breakpoints[8] = ler_breakpoints[7] + straights["LLL5"].l;
 
  // negative s
  ler_breakpoints[9] = -straights["LLR1"].l;
  ler_breakpoints[10] = ler_breakpoints[9] - bendings["BC1RP"].rt * bendings["BC1RP"].dphi;
  ler_breakpoints[11] = ler_breakpoints[10] - straights["LLR2"].l;
  ler_breakpoints[12] = ler_breakpoints[11] - bendings["BLCWRP"].rt * bendings["BLCWRP"].dphi;
  ler_breakpoints[13] = ler_breakpoints[12] - straights["LLR3"].l;
  ler_breakpoints[14] = ler_breakpoints[13] - bendings["BLC1RP"].rt * bendings["BLC1RP"].dphi;
  ler_breakpoints[15] = ler_breakpoints[14] - straights["LLR4"].l;
  ler_breakpoints[16] = ler_breakpoints[15] - bendings["BLC2RP"].rt * bendings["BLC2RP"].dphi;
  ler_breakpoints[17] = ler_breakpoints[16] - bendings["BLC2RP.2"].rt * bendings["BLC2RP.2"].dphi;
  ler_breakpoints[18] = ler_breakpoints[17] - straights["LLR5"].l;
  ler_breakpoints[19] = ler_breakpoints[18] - bendings["BLY2RP.2"].rt * bendings["BLY2RP.2"].dphi;
  ler_breakpoints[20] = ler_breakpoints[20] - straights["LLR6"].l;
 
 
  double dx = 0;
  double dz = 0;
  double phi = 0;
  if (accRing == c_LER) {
    // LER
    // positive s
    if (400.0 * Unit::cm < s) {
      if (s < ler_breakpoints[0]) {
        phi = straights["LLL1"].phi;
        dx = straights["LLL1"].x0 + s * sin(phi);
        dz = straights["LLL1"].z0 + s * cos(phi);
      } else if (s < ler_breakpoints[1]) {
        double sloc = s - ler_breakpoints[0];
        phi = bendings["BC1LP"].sphi + sloc / bendings["BC1LP"].rt;
        // Torus is created in x-y plain.
        // It is then rotated to x-z plain,
        // and its direction changes to reversed,
        // thus phi_real=-phi_xml
        phi = -phi;
        dx = bendings["BC1LP"].x0 + bendings["BC1LP"].rt * cos(-phi);
        dz = bendings["BC1LP"].z0 + bendings["BC1LP"].rt * sin(-phi);
      } else if (s < ler_breakpoints[2]) {
        double sloc = s - ler_breakpoints[1];
        phi = straights["LLL2"].phi;
        dx = straights["LLL2"].x0 + sloc * sin(phi);
        dz = straights["LLL2"].z0 + sloc * cos(phi);
      } else if (s < ler_breakpoints[3]) {
        double sloc = s - ler_breakpoints[2];
        phi = bendings["BLC1LP1"].sphi + sloc / bendings["BLC1LP1"].rt;
        phi = -phi;
        dx = bendings["BLC1LP1"].x0 + bendings["BLC1LP1"].rt * cos(-phi);
        dz = bendings["BLC1LP1"].z0 + bendings["BLC1LP1"].rt * sin(-phi);
      } else if (s < ler_breakpoints[4]) {
        double sloc = s - ler_breakpoints[3];
        phi = straights["LLL3"].phi;
        dx = straights["LLL3"].x0 + sloc * sin(phi);
        dz = straights["LLL3"].z0 + sloc * cos(phi);
      } else if (s < ler_breakpoints[5]) {
        double sloc = s - ler_breakpoints[4];
        // Torus dphi may be only positive,
        // while direction of increasing |s| is sometimes negative,
        // and we need to use -s and not change phi.
        // Since we add pi to phi later,
        // we subtract it now for this element.
        phi = bendings["BLC1LP2"].sphi + bendings["BLC1LP2"].dphi - sloc / bendings["BLC1LP2"].rt;
        phi = -phi;
        dx = bendings["BLC1LP2"].x0 + bendings["BLC1LP2"].rt * cos(-phi);
        dz = bendings["BLC1LP2"].z0 + bendings["BLC1LP2"].rt * sin(-phi);
        phi -= M_PI;
      } else if (s < ler_breakpoints[6]) {
        double sloc = s - ler_breakpoints[5];
        phi = straights["LLL4"].phi;
        dx = straights["LLL4"].x0 + sloc * sin(phi);
        dz = straights["LLL4"].z0 + sloc * cos(phi);
      } else if (s < ler_breakpoints[8]) {
        double sloc = s - ler_breakpoints[7];
        phi = straights["LLL5"].phi;
        dx = straights["LLL5"].x0 + sloc * sin(phi);
        dz = straights["LLL5"].z0 + sloc * cos(phi);
      }
      // For this direction rotation angle of elements changes to negative,
      // while SAD coordinates keep orientation.
      // We need to compensate.
      phi += M_PI;
    }
    // negative s
    else if (s < -400.0 * Unit::cm) {
      if (s > ler_breakpoints[9]) {
        double sloc = -s;
        phi = straights["LLR1"].phi;
        dx = straights["LLR1"].x0 + sloc * sin(phi);
        dz = straights["LLR1"].z0 + sloc * cos(phi);
      } else if (s > ler_breakpoints[10]) {
        double sloc = ler_breakpoints[9] - s;
        phi = bendings["BC1RP"].sphi + bendings["BC1RP"].dphi - sloc / bendings["BC1RP"].rt;
        phi = -phi;
        dx = bendings["BC1RP"].x0 + bendings["BC1RP"].rt * cos(-phi);
        dz = bendings["BC1RP"].z0 + bendings["BC1RP"].rt * sin(-phi);
        phi += M_PI;
      } else if (s > ler_breakpoints[11]) {
        double sloc = ler_breakpoints[10] - s;
        phi = straights["LLR2"].phi;
        dx = straights["LLR2"].x0 + sloc * sin(phi);
        dz = straights["LLR2"].z0 + sloc * cos(phi);
      } else if (s > ler_breakpoints[12]) {
        double sloc = ler_breakpoints[11] - s;
        phi = bendings["BLCWRP"].sphi + bendings["BLCWRP"].dphi - sloc / bendings["BLCWRP"].rt;
        phi = -phi;
        dx = bendings["BLCWRP"].x0 + bendings["BLCWRP"].rt * cos(-phi);
        dz = bendings["BLCWRP"].z0 + bendings["BLCWRP"].rt * sin(-phi);
        phi += M_PI;
      } else if (s > ler_breakpoints[13]) {
        double sloc = ler_breakpoints[12] - s;
        phi = straights["LLR3"].phi;
        dx = straights["LLR3"].x0 + sloc * sin(phi);
        dz = straights["LLR3"].z0 + sloc * cos(phi);
      } else if (s > ler_breakpoints[14]) {
        double sloc = ler_breakpoints[13] - s;
        phi = bendings["BLC1RP"].sphi + bendings["BLC1RP"].dphi - sloc / bendings["BLC1RP"].rt;
        phi = -phi;
        dx = bendings["BLC1RP"].x0 + bendings["BLC1RP"].rt * cos(-phi);
        dz = bendings["BLC1RP"].z0 + bendings["BLC1RP"].rt * sin(-phi);
        phi += M_PI;
      } else if (s > ler_breakpoints[15]) {
        double sloc = ler_breakpoints[14] - s;
        phi = straights["LLR4"].phi;
        dx = straights["LLR4"].x0 + sloc * sin(phi);
        dz = straights["LLR4"].z0 + sloc * cos(phi);
      } else if (s > ler_breakpoints[16]) {
        double sloc = ler_breakpoints[15] - s;
        phi = bendings["BLC2RP"].sphi + sloc / bendings["BLC2RP"].rt;
        phi = -phi;
        dx = bendings["BLC2RP"].x0 + bendings["BLC2RP"].rt * cos(-phi);
        dz = bendings["BLC2RP"].z0 + bendings["BLC2RP"].rt * sin(-phi);
      } else if (s > ler_breakpoints[17]) {
        double sloc = ler_breakpoints[16] - s;
        phi = bendings["BLC2RP.2"].sphi + sloc / bendings["BLC2RP.2"].rt;
        phi = -phi;
        dx = bendings["BLC2RP.2"].x0 + bendings["BLC2RP.2"].rt * cos(-phi);
        dz = bendings["BLC2RP.2"].z0 + bendings["BLC2RP.2"].rt * sin(-phi);
      } else if (s > ler_breakpoints[18]) {
        double sloc = ler_breakpoints[17] - s;
        phi = straights["LLR5"].phi;
        dx = straights["LLR5"].x0 + sloc * sin(phi);
        dz = straights["LLR5"].z0 + sloc * cos(phi);
      } else if (s > ler_breakpoints[19]) {
        double sloc = ler_breakpoints[18] - s;
        phi = bendings["BLY2RP.2"].sphi + sloc / bendings["BLY2RP.2"].rt;
        phi = -phi;
        dx = bendings["BLY2RP.2"].x0 + bendings["BLY2RP.2"].rt * cos(-phi);
        dz = bendings["BLY2RP.2"].z0 + bendings["BLY2RP.2"].rt * sin(-phi);
      } else if (s > ler_breakpoints[20]) {
        double sloc = ler_breakpoints[19] - s;
        phi = straights["LLR6"].phi;
        dx = straights["LLR6"].x0 + sloc * sin(phi);
        dz = straights["LLR6"].z0 + sloc * cos(phi);
      }
    }
  }
  if (accRing == c_HER) {
    // HER
    // positive s
    if (400.0 * Unit::cm < s) {
      if (s < her_breakpoints[0]) {
        phi = straights["LHL1"].phi;
        dx = straights["LHL1"].x0 + s * sin(phi);
        dz = straights["LHL1"].z0 + s * cos(phi);
      } else if (s < her_breakpoints[1]) {
        double sloc = s - her_breakpoints[0];
        phi = bendings["BLC1LE"].sphi + sloc / bendings["BLC1LE"].rt;
        phi = -phi;
        dx = bendings["BLC1LE"].x0 + bendings["BLC1LE"].rt * cos(-phi);
        dz = bendings["BLC1LE"].z0 + bendings["BLC1LE"].rt * sin(-phi);
      } else if (s < her_breakpoints[2]) {
        double sloc = s - her_breakpoints[1];
        phi = straights["LHL2"].phi;
        dx = straights["LHL2"].x0 + sloc * sin(phi);
        dz = straights["LHL2"].z0 + sloc * cos(phi);
      }
      phi += M_PI;
    }
    // negative s
    else if (s < -400.0 * Unit::cm) {
      if (s > her_breakpoints[3]) {
        double sloc = -s;
        phi = straights["LHR1"].phi;
        dx = straights["LHR1"].x0 + sloc * sin(phi);
        dz = straights["LHR1"].z0 + sloc * cos(phi);
      } else if (s > her_breakpoints[4]) {
        double sloc = her_breakpoints[3] - s;
        phi = bendings["BLC2RE"].sphi + sloc / bendings["BLC2RE"].rt;
        phi = -phi;
        dx = bendings["BLC2RE"].x0 + bendings["BLC2RE"].rt * cos(-phi);
        dz = bendings["BLC2RE"].z0 + bendings["BLC2RE"].rt * sin(-phi);
      } else if (s > her_breakpoints[5]) {
        double sloc = her_breakpoints[4] - s;
        phi = straights["LHR2"].phi;
        dx = straights["LHR2"].x0 + sloc * sin(phi);
        dz = straights["LHR2"].z0 + sloc * cos(phi);
      }
    }
  }
 
  TGeoHMatrix matrix("SADTrafo");
  matrix.RotateY(phi / Unit::deg);
  matrix.SetDx(dx);
  matrix.SetDz(dz);
  return matrix;
}

setMomentumRes()

void setMomentumRes ( double  pxRes, double  pyRes  )

inline

Sets the resolution of the momentum for the real particles.

Allows each real particle to be smeared according to the specified momentum resolution.

Parameters

pxRes	The resolution for the x momentum component.
pyRes	The resolution for the y momentum component.

Definition at line 82 of file ReaderSAD.h.

{ m_pxRes = pxRes; m_pyRes = pyRes; }

Member Data Documentation

m_accRing

AcceleratorRings m_accRing

protected

The accelerator ring from which the particles originate.

Definition at line 125 of file ReaderSAD.h.

m_file

TFile* m_file

protected

The input root file.

Definition at line 120 of file ReaderSAD.h.

m_inputSAD_dp_over_p0

double m_inputSAD_dp_over_p0

protected

dp_over_p0

Definition at line 158 of file ReaderSAD.h.

m_inputSAD_E

double m_inputSAD_E

protected

energy at lost position [m].

Definition at line 159 of file ReaderSAD.h.

m_inputSAD_Lss

double m_inputSAD_Lss

protected

length of element in which scattered [m]

Definition at line 148 of file ReaderSAD.h.

m_inputSAD_nturn

int m_inputSAD_nturn

protected

number of turns from scattered to lost

Definition at line 149 of file ReaderSAD.h.

m_inputSAD_px

double m_inputSAD_px

protected

x momentum at lost position [m].

Definition at line 152 of file ReaderSAD.h.

m_inputSAD_py

double m_inputSAD_py

protected

y momentum at lost position [m].

Definition at line 153 of file ReaderSAD.h.

m_inputSAD_r

double m_inputSAD_r

protected

sqrt(x*x+y*y) [m]

Definition at line 156 of file ReaderSAD.h.

m_inputSAD_rate

double m_inputSAD_rate

protected

loss rate [Hz]

Definition at line 160 of file ReaderSAD.h.

m_inputSAD_rr

double m_inputSAD_rr

protected

sqrt(x*x+y*y) [m] before matching onto beam pipe inner surface

Definition at line 157 of file ReaderSAD.h.

m_inputSAD_s

double m_inputSAD_s

protected

lost position (|s|<Ltot/2) [m]

Definition at line 147 of file ReaderSAD.h.

m_inputSAD_sraw

double m_inputSAD_sraw

protected

lost position [m

Definition at line 145 of file ReaderSAD.h.

m_inputSAD_ss

double m_inputSAD_ss

protected

scattered position (|s|<Ltot/2) [m]

Definition at line 146 of file ReaderSAD.h.

m_inputSAD_ssraw

double m_inputSAD_ssraw

protected

scattered position [m]

Definition at line 144 of file ReaderSAD.h.

m_inputSAD_watt

double m_inputSAD_watt

protected

loss wattage [W]

Definition at line 161 of file ReaderSAD.h.

m_inputSAD_x

double m_inputSAD_x

protected

x at lost position [m].

Definition at line 150 of file ReaderSAD.h.

m_inputSAD_xraw

double m_inputSAD_xraw

protected

x at lost position [m] before matching onto beam pipe inner surface

Definition at line 154 of file ReaderSAD.h.

m_inputSAD_y

double m_inputSAD_y

protected

y at lost position [m].

Definition at line 151 of file ReaderSAD.h.

m_inputSAD_yraw

double m_inputSAD_yraw

protected

y at lost position [m] before matching onto beam pipe inner surface

Definition at line 155 of file ReaderSAD.h.

m_lostE

double m_lostE

protected

energy at lost position [m].

Definition at line 142 of file ReaderSAD.h.

m_lostPx

double m_lostPx

protected

x momentum at lost position [m].

Definition at line 139 of file ReaderSAD.h.

m_lostPy

double m_lostPy

protected

y momentum at lost position [m].

Definition at line 140 of file ReaderSAD.h.

m_lostRate

double m_lostRate

protected

loss rate [Hz]>

Definition at line 141 of file ReaderSAD.h.

m_lostS

double m_lostS

protected

lost position [m] along ring.

0 is the IP. Range goes from -L/2 to L/2, where L is the ring circumference.

Definition at line 138 of file ReaderSAD.h.

m_lostX

double m_lostX

protected

x at lost position [m].

Definition at line 136 of file ReaderSAD.h.

m_lostY

double m_lostY

protected

y at lost position [m].

Definition at line 137 of file ReaderSAD.h.

m_pxRes

double m_pxRes

protected

The resolution for the x momentum component of the SAD real particle.

Definition at line 126 of file ReaderSAD.h.

m_pyRes

double m_pyRes

protected

The resolution for the y momentum component of the SAD real particle.

Definition at line 127 of file ReaderSAD.h.

m_readEntry

int m_readEntry

protected

The current number of the SAD entry that is read.

Definition at line 134 of file ReaderSAD.h.

m_readoutTime

double m_readoutTime

protected

The readout time.

Definition at line 130 of file ReaderSAD.h.

m_realPartEntry

unsigned int m_realPartEntry

protected

The current number of the created real particles.

Definition at line 133 of file ReaderSAD.h.

m_realPartNum

unsigned int m_realPartNum

protected

The current number of the created real particles.

Definition at line 132 of file ReaderSAD.h.

m_SADToRealFactor

double m_SADToRealFactor

protected

The factor to calculate the number of real particles from a SAD particle.

Definition at line 129 of file ReaderSAD.h.

m_sRange

double m_sRange

protected

The +- range for the s value for which particles are loaded.

Definition at line 124 of file ReaderSAD.h.

m_transMatrix

TGeoHMatrix* m_transMatrix

protected

Transformation matrix from local SAD to global geant4 space.

Definition at line 123 of file ReaderSAD.h.

m_tree

TTree* m_tree

protected

The input root tree.

Definition at line 121 of file ReaderSAD.h.

---

The documentation for this class was generated from the following files:

• generators/SAD/include/ReaderSAD.h
• generators/SAD/src/ReaderSAD.cc