Belle II Software  release-05-01-25
ReaderSAD.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2016 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll, Hiroyuki Nakayama, Igal Jaegle *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/SAD/ReaderSAD.h>
12 
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Unit.h>
15 #include <framework/gearbox/GearDir.h>
16 #include <mdst/dataobjects/MCParticle.h>
17 
18 #include <generators/SAD/dataobjects/SADMetaHit.h>
19 
20 // framework - DataStore
21 #include <framework/datastore/StoreArray.h>
22 
23 #include <TRandom3.h>
24 
25 using namespace std;
26 using namespace Belle2;
27 
28 
29 ReaderSAD::ReaderSAD(): m_file(NULL), m_tree(NULL), m_transMatrix(NULL),
30  m_sRange(3000.0), m_accRing(ReaderSAD::c_LER), m_pxRes(0.01), m_pyRes(0.01),
31  m_SADToRealFactor(5.76e6), m_readoutTime(20.0), m_realPartNum(0),
32  m_realPartEntry(0), m_readEntry(0)
33 {
34  m_lostX = 0.0;
35  m_lostY = 0.0;
36  m_lostS = 0.0;
37  m_lostPx = 0.0;
38  m_lostPy = 0.0;
39  m_lostRate = 0.0;
40  m_lostE = 0.0;
41 
42  //Start my addition
43  m_inputSAD_ssraw = 0;
44  m_inputSAD_sraw = 0;
45  m_inputSAD_ss = 0;
46  m_inputSAD_s = 0;
47  m_inputSAD_Lss = 0;
48  m_inputSAD_nturn = 0;
49  m_inputSAD_x = 0;
50  m_inputSAD_y = 0;
51  m_inputSAD_px = 0;
52  m_inputSAD_py = 0;
53  m_inputSAD_xraw = 0;
54  m_inputSAD_yraw = 0;
55  m_inputSAD_r = 0;
56  m_inputSAD_rr = 0;
58  m_inputSAD_E = 0;
59  m_inputSAD_rate = 0;
60  m_inputSAD_watt = 0;
61  //End my addition
62 }
63 
64 
66 {
67  // comment out the following line for the march 16 build (TF)
68  // if (m_file != NULL) m_file->Close();
69 }
70 
71 
72 void ReaderSAD::initialize(TGeoHMatrix* transMatrix, double sRange, ReaderSAD::AcceleratorRings accRing, double readoutTime)
73 {
74  m_transMatrix = transMatrix;
75  m_sRange = sRange;
76  m_accRing = accRing;
77  m_readoutTime = readoutTime;
78 }
79 
80 
81 void ReaderSAD::open(const string& filename)
82 {
83  if (m_file != NULL) {
84  m_file->Close();
85  delete m_file;
86  }
87 
88  m_file = new TFile(filename.c_str(), "READ");
89  if (m_file == NULL) throw(SADCouldNotOpenFileError() << filename);
90 
91  m_file->cd("");
92  m_tree = dynamic_cast<TTree*>(m_file->Get("tp"));
93  if (m_tree == NULL) throw(SADCouldNotOpenFileError() << filename);
94 
95  m_tree->SetBranchAddress("x", &m_lostX);
96  m_tree->SetBranchAddress("y", &m_lostY);
97  m_tree->SetBranchAddress("s", &m_lostS);
98  m_tree->SetBranchAddress("px", &m_lostPx);
99  m_tree->SetBranchAddress("py", &m_lostPy);
100  m_tree->SetBranchAddress("rate", &m_lostRate);
101  m_tree->SetBranchAddress("E", &m_lostE);
102 
103  m_tree->SetBranchAddress("ssraw", &m_inputSAD_ssraw);
104  m_tree->SetBranchAddress("sraw", &m_inputSAD_sraw);
105  m_tree->SetBranchAddress("ss", &m_inputSAD_ss);
106  m_tree->SetBranchAddress("Lss", &m_inputSAD_Lss);
107  m_tree->SetBranchAddress("nturn", &m_inputSAD_nturn);
108  m_tree->SetBranchAddress("xraw", &m_inputSAD_xraw);
109  m_tree->SetBranchAddress("yraw", &m_inputSAD_yraw);
110  m_tree->SetBranchAddress("r", &m_inputSAD_r);
111  m_tree->SetBranchAddress("rr", &m_inputSAD_rr);
112  m_tree->SetBranchAddress("dp_over_p0", &m_inputSAD_dp_over_p0);
113  m_tree->SetBranchAddress("watt", &m_inputSAD_watt);
114 
115  m_readEntry = -1;
116 
118  hit.registerInDataStore();
119 
120 }
121 
122 
124 {
125  if (m_tree == NULL) {
126  B2ERROR("The SAD tree doesn't exist !");
127  return -1;
128  }
129 
130  do {
131  m_readEntry++;
132 
133  //Check for end of file
134  if (m_readEntry >= m_tree->GetEntries()) throw SADEndOfFile();
135 
136  //Load the SAD particle
137  m_tree->GetEntry(m_readEntry);
139  B2DEBUG(10, "> Read particle " << m_readEntry + 1 << "/" << m_tree->GetEntries() << " with s = " << m_lostS << " cm" <<
140  " and rate = " << m_lostRate << " Hz");
141 
142  //printf("Read particle %d / %d with s= %f [m]\n", m_readEntry + 1 , (int)m_tree->GetEntries(), m_lostS / 100.);
143 
144  double zMom2 = m_lostE * m_lostE - m_lostPx * m_lostPx - m_lostPy * m_lostPy ;
145  if (zMom2 < 0) printf("zMom2= %f is negative. Skipped!\n", zMom2);
146 
147  } while ((fabs(m_lostS) > m_sRange) or (m_lostE * m_lostE - m_lostPx * m_lostPx - m_lostPy * m_lostPy < 0));
148 
149 // do {
150 // //Check for end of file
151 // m_readEntry++;
152 // if (m_readEntry >= m_tree->GetEntries()) throw SADEndOfFile();
153 //
154 // //Load the SAD particle
155 // m_tree->GetEntry(m_readEntry);
156 // convertParamsToSADUnits();
157 //
158 // B2DEBUG(10, "> Read particle " << m_readEntry + 1 << "/" << m_tree->GetEntries() << " with s = " << m_lostS << " cm" << " and rate = " << m_lostRate << " Hz" )
159 // } while (fabs(m_lostS) > m_sRange);
160 
162  return m_lostRate;
163 }
164 
165 
167 {
168  if (m_tree == NULL) {
169  B2ERROR("The SAD tree doesn't exist !");
170  return false;
171  }
172 
173  //Check for end of file
174  if ((m_readEntry >= m_tree->GetEntries()) || (m_tree->GetEntries() == 0)) throw SADEndOfFile();
175 
176  //Check if the number of the real particles is reached
178  m_realPartEntry = 0;
179  m_realPartNum = 0;
180  }
181 
182  //Check if a new SAD particle has to be read from the file
183  if (m_realPartNum == 0) {
184  //Read only SAD particles which are inside the chosen sRange
185  do {
186  m_readEntry++;
187  if (m_readEntry >= m_tree->GetEntries()) throw SADEndOfFile();
188 
189  m_tree->GetEntry(m_readEntry);
191 
192  B2DEBUG(10, "> Read particle " << m_readEntry + 1 << "/" << m_tree->GetEntries() << " with s = " << m_lostS << " cm" <<
193  " and rate = " << m_lostRate << " Hz");
194  } while ((fabs(m_lostS) > m_sRange) && (m_readEntry < m_tree->GetEntries()));
195 
196  if (fabs(m_lostS) <= m_sRange)
198  }
199 
200  //Create a new real particle from the SAD particle
201  if ((fabs(m_lostS) <= m_sRange) && (m_realPartNum > 0)) {
203  B2DEBUG(10, "* Created real particle " << m_realPartEntry + 1 << "/" << m_realPartNum << " for SAD particle " << m_readEntry + 1 <<
204  "/" << m_tree->GetEntries());
205  }
206 
207  m_realPartEntry++;
208 
209  return true;
210 }
211 
212 
214 {
215  if (m_tree == NULL) {
216  B2ERROR("The SAD tree doesn't exist !");
217  return ;
218  }
219 
220  int nPart = m_tree->GetEntries();
221 
222  for (int iPart = 0; iPart < nPart; ++iPart) {
223  m_tree->GetEntry(iPart);
225  if (fabs(m_lostS) <= m_sRange) addParticleToMCParticles(graph);
226 
227  //std::cout << " in sad " << m_inputSAD_sraw << std::endl;
228 
229  }
230 }
231 
232 
233 //======================================================================
234 // Private methods
235 //======================================================================
236 
238 {
239  m_lostX = m_lostX * Unit::m;
240  m_lostY = m_lostY * Unit::m;
241  m_lostS = m_lostS * Unit::m;
244 
245  //Start my addition
250  //m_inputSAD_Lss
251  //m_inputSAD_nturn
258  //m_inputSAD_r
259  //m_inputSAD_rr
260  //m_inputSAD_dp_over_p0
261  //m_inputSAD_E
262  //m_inputSAD_rate
263  //End my addition
264 }
265 
266 
267 void ReaderSAD::addParticleToMCParticles(MCParticleGraph& graph, bool gaussSmearing)
268 {
269  double particlePosSAD[3] = {0.0, 0.0, 0.0};
270  double particlePosSADfar[3] = {0.0, 0.0, 0.0};
271  double particlePosGeant4[3] = {0.0, 0.0, 0.0};
272  double particleMomSAD[3] = {0.0, 0.0, 0.0};
273  double particleMomGeant4[3] = {0.0, 0.0, 0.0};
274 
275  //Add particle to MCParticle collection
276  MCParticleGraph::GraphParticle& particle = graph.addParticle();
277  particle.setStatus(MCParticle::c_PrimaryParticle);
278 
279  switch (m_accRing) {
280  case c_HER: particle.setPDG(11); //electrons
281  break;
282  case c_LER: particle.setPDG(-11); //positrons
283  break;
284  }
285 
286  particle.setMassFromPDG();
287 
288  //Convert the position of the particle from local SAD space to global geant4 space.
289  //Flip the sign for the y and z component to go from the accelerator to the detector coordinate system.
290  particlePosSAD[0] = m_lostX;
291  particlePosSAD[1] = -m_lostY;
292  particlePosSAD[2] = -m_lostS;
293 
294  particlePosSADfar[0] = m_lostX;
295  particlePosSADfar[1] = -m_lostY;
296  particlePosSADfar[2] = 0;
297 
298 
299  static GearDir content = Gearbox::getInstance().getDetectorComponent("FarBeamLine");
300  if (!content)
301  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.");
302 
303  TGeoHMatrix* m_transMatrix2 = new TGeoHMatrix(SADtoGeant(m_accRing, m_lostS)); //overwrite m_transMatrix given by initialize()
304 
305  if (abs(m_lostS) <= 400.) { //4m
306  m_transMatrix->LocalToMaster(particlePosSAD, particlePosGeant4);
307  } else {
308  m_transMatrix2->LocalToMaster(particlePosSADfar, particlePosGeant4);
309  }
310 
311  //Convert the momentum of the particle from local SAD space to global geant4 space.
312  //Flip the sign for the y and z component to go from the accelerator to the detector coordinate system.
313  //Calculate the missing pz by using the energy of the particle at the position where it has been lost.
314  //double totalMomSqr = (m_lostE * m_lostE - (particle.getMass() * particle.getMass()));
315  double totalMomSqr = m_lostE * m_lostE; // SAD output "E" is "P", in fact.
316 
317  if (gaussSmearing) {
318  particleMomSAD[0] = gRandom->Gaus(m_lostPx, m_pxRes * m_lostPx); //1% px resolution
319  particleMomSAD[1] = -1.0 * gRandom->Gaus(m_lostPy, m_pyRes * m_lostPy);
320  } else {
321  particleMomSAD[0] = m_lostPx;
322  particleMomSAD[1] = -m_lostPy;
323  }
324 
325  double zMom = sqrt(totalMomSqr - (particleMomSAD[0] * particleMomSAD[0]) - (particleMomSAD[1] * particleMomSAD[1]));
326 
327  switch (m_accRing) {
328  case c_HER: particleMomSAD[2] = zMom;
329  break;
330  case c_LER: particleMomSAD[2] = -zMom;
331  break;
332  }
333 
334  int ring = 0;
335 
336  switch (m_accRing) {
337  case c_HER: ring = 1;
338  break;
339  case c_LER: ring = 2;
340  break;
341  }
342  /*
343  int ler_section[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
344  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};
345  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};
346  int her_section[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
347  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};
348  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};
349 
350  int section = -1;
351 
352  if (ring == 1) {
353  for (int i = 0; i < 12; i++) {
354  if (her_inf_section[i] <= (m_inputSAD_ssraw + 1500.) && (m_inputSAD_ssraw + 1500.) <= her_sup_section[i])
355  section = her_section[i] - 1;
356  }
357  } else if (ring == 2) {
358  for (int i = 0; i < 12; i++) {
359  if (ler_inf_section[i] <= (m_inputSAD_ssraw + 1500.) && (m_inputSAD_ssraw + 1500.) <= ler_sup_section[i])
360  section = ler_section[i] - 1;
361  }
362  }
363  */
364  //each rings have 12 section of ~250m
365  //the 1st section D01, the second section is D12, followed by D11, D10 .... for both rings
366  int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
367  //int ring_section = section_ordering[(int)((m_inputSAD_ssraw + 1500.) / 12.)];
368  double ssraw = 0.;
369  if (ring == 1) {
370  if (m_inputSAD_ssraw >= 0) ssraw = m_inputSAD_ssraw / 100.;
371  else ssraw = 3000. + m_inputSAD_ssraw / 100.;
372  } else if (ring == 2) {
373  //if (m_inputSAD_ssraw >= 0) ssraw = 3000. - m_inputSAD_ssraw / 100.;
374  //else if (m_inputSAD_ssraw < 0) ssraw = -m_inputSAD_ssraw / 100.;
375  if (m_inputSAD_ssraw >= 0) ssraw = m_inputSAD_ssraw / 100.;
376  else ssraw = 3000. + m_inputSAD_ssraw / 100.;
377  }
378  int ring_section = section_ordering[(int)((ssraw) / 250.)];
379 
380  if (abs(m_lostS) <= 400.) {
381  m_transMatrix->LocalToMasterVect(particleMomSAD, particleMomGeant4);
382  } else {
383  m_transMatrix2->LocalToMasterVect(particleMomSAD, particleMomGeant4);
384  }
385 
386  //Set missing particle information
387  particle.setMomentum(TVector3(particleMomGeant4));
388  particle.setProductionVertex(TVector3(particlePosGeant4));
389  particle.setProductionTime(0.0);
390  particle.setEnergy(sqrt(m_lostE * m_lostE + particle.getMass()*particle.getMass()));
391  particle.setValidVertex(true);
392 
393  //Start my addition
394  StoreArray<SADMetaHit> SADMetaHits;
399  m_inputSAD_watt, ring, ring_section));
400  //End my addition
401 }
402 
403 
405 {
406  double numPart = rate * m_readoutTime / Unit::s; //readoutTime [ns]
407 
408  int numPart_int = static_cast<int>(floor(numPart));
409  double numPart_dec = numPart - numPart_int;
410 
411  double rnd = gRandom->Uniform(); //returns a random number in the interval ]0, 1]
412  if (rnd < numPart_dec) numPart_int += 1;
413 
414  return numPart_int;
415 }
416 
418 {
419  // 0<sraw<3016.3145 m
420  // -1500<s<1500 m
421 
422  //static double max_s_her = 3016.3145 * Unit::m;
423  //static double max_s_ler = 3016.3026 * Unit::m;
424 
425  //get parameters from .xml file
426  static GearDir content = Gearbox::getInstance().getDetectorComponent("FarBeamLine");
427 
428  map<string, straightElement> straights;
429  map<string, bendingElement> bendings;
430  for (const GearDir& element : content.getNodes("Straight")) {
431 
432  string name = element.getString("@name");
433  string type = element.getString("@type");
434 
435  if (type != "pipe") continue;
436 
437  straightElement straight;
438 
439  straight.x0 = element.getLength("X0");
440  straight.z0 = element.getLength("Z0");
441  straight.l = element.getLength("L");
442  straight.phi = element.getLength("PHI");
443 
444  straights[name] = straight;
445  }
446 
447  string str_checklist[] = {"LHR1", "LHR2", "LLR1", "LLR2", "LLR3", "LLR4", "LLR5", "LHL1", "LHL2", "LLL1", "LLL2", "LLL3", "LLL4"};
448  for (const string& str : str_checklist) {
449  if (straights.count(str) == 0)
450  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.");
451  }
452 
453  for (const GearDir& element : content.getNodes("Bending")) {
454 
455  string name = element.getString("@name");
456  string type = element.getString("@type");
457 
458  if (type != "pipe") continue;
459 
460  bendingElement bending;
461 
462  bending.rt = element.getLength("RT");
463  bending.x0 = element.getLength("X0");
464  bending.z0 = element.getLength("Z0");
465  bending.sphi = element.getLength("SPHI");
466  bending.dphi = element.getLength("DPHI");
467 
468  bendings[name] = bending;
469  }
470 
471  string bend_checklist[] = {"BLC2RE", "BC1RP", "BLCWRP", "BLC1RP", "BLC2RP", "BLC1LE", "BC1LP", "BLC1LP", "BLC2LP"};
472  for (const string& bnd : bend_checklist) {
473  if (bendings.count(bnd) == 0)
474  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.");
475  }
476 
477  static double her_breakpoints[6];
478  static double ler_breakpoints[16];
479 
480  // positive s
481  her_breakpoints[0] = straights["LHL1"].l;
482  her_breakpoints[1] = her_breakpoints[0] + bendings["BLC1LE"].rt * bendings["BLC1LE"].dphi;
483  her_breakpoints[2] = her_breakpoints[1] + straights["LHL2"].l;
484 
485  // negative s
486  her_breakpoints[3] = -straights["LHR1"].l;
487  her_breakpoints[4] = her_breakpoints[3] - bendings["BLC2RE"].rt * bendings["BLC2RE"].dphi;
488  her_breakpoints[5] = her_breakpoints[4] - straights["LHR2"].l;
489 
490  // positive s
491  ler_breakpoints[0] = straights["LLL1"].l;
492  ler_breakpoints[1] = ler_breakpoints[0] + bendings["BC1LP"].rt * bendings["BC1LP"].dphi;
493  ler_breakpoints[2] = ler_breakpoints[1] + straights["LLL2"].l;
494  ler_breakpoints[3] = ler_breakpoints[2] + bendings["BLC1LP"].rt * bendings["BLC1LP"].dphi;
495  ler_breakpoints[4] = ler_breakpoints[3] + straights["LLL3"].l;
496  ler_breakpoints[5] = ler_breakpoints[4] + bendings["BLC2LP"].rt * bendings["BLC2LP"].dphi;
497  ler_breakpoints[6] = ler_breakpoints[5] + straights["LLL4"].l;
498 
499  // negative s
500  ler_breakpoints[7] = -straights["LLR1"].l;
501  ler_breakpoints[8] = ler_breakpoints[7] - bendings["BC1RP"].rt * bendings["BC1RP"].dphi;
502  ler_breakpoints[9] = ler_breakpoints[8] - straights["LLR2"].l;
503  ler_breakpoints[10] = ler_breakpoints[9] - bendings["BLCWRP"].rt * bendings["BLCWRP"].dphi;
504  ler_breakpoints[11] = ler_breakpoints[10] - straights["LLR3"].l;
505  ler_breakpoints[12] = ler_breakpoints[11] - bendings["BLC1RP"].rt * bendings["BLC1RP"].dphi;
506  ler_breakpoints[13] = ler_breakpoints[12] - straights["LLR4"].l;
507  ler_breakpoints[14] = ler_breakpoints[13] - bendings["BLC2RP"].rt * bendings["BLC2RP"].dphi;
508  ler_breakpoints[15] = ler_breakpoints[14] - straights["LLR5"].l;
509 
510  double dx = 0;
511  double dz = 0;
512  double phi = 0;
513  if (accRing == c_LER) {
514  // LER
515  // positive s
516  if (400.0 * Unit::cm < s) {
517  if (s < ler_breakpoints[0]) {
518  phi = straights["LLL1"].phi;
519  dx = straights["LLL1"].x0 + s * sin(phi);
520  dz = straights["LLL1"].z0 + s * cos(phi);
521  } else if (s < ler_breakpoints[1]) {
522  double sloc = s - ler_breakpoints[0];
523  phi = bendings["BC1LP"].sphi + sloc / bendings["BC1LP"].rt;
524  // Torus is created in x-y plain.
525  // It is then rotated to x-z plain,
526  // and its direction changes to reversed,
527  // thus phi_real=-phi_xml
528  phi = -phi;
529  dx = bendings["BC1LP"].x0 + bendings["BC1LP"].rt * cos(-phi);
530  dz = bendings["BC1LP"].z0 + bendings["BC1LP"].rt * sin(-phi);
531  } else if (s < ler_breakpoints[2]) {
532  double sloc = s - ler_breakpoints[1];
533  phi = straights["LLL2"].phi;
534  dx = straights["LLL2"].x0 + sloc * sin(phi);
535  dz = straights["LLL2"].z0 + sloc * cos(phi);
536  } else if (s < ler_breakpoints[3]) {
537  double sloc = s - ler_breakpoints[2];
538  phi = bendings["BLC1LP"].sphi + sloc / bendings["BLC1LP"].rt;
539  phi = -phi;
540  dx = bendings["BLC1LP"].x0 + bendings["BLC1LP"].rt * cos(-phi);
541  dz = bendings["BLC1LP"].z0 + bendings["BLC1LP"].rt * sin(-phi);
542  } else if (s < ler_breakpoints[4]) {
543  double sloc = s - ler_breakpoints[3];
544  phi = straights["LLL3"].phi;
545  dx = straights["LLL3"].x0 + sloc * sin(phi);
546  dz = straights["LLL3"].z0 + sloc * cos(phi);
547  } else if (s < ler_breakpoints[5]) {
548  double sloc = s - ler_breakpoints[4];
549  // Torus dphi may be only positive,
550  // while direction of increasing |s| is sometimes negative,
551  // and we need to use -s and not change phi.
552  // Since we add pi to phi later,
553  // we subtract it now for this element.
554  phi = bendings["BLC2LP"].sphi + bendings["BLC2LP"].dphi - sloc / bendings["BLC2LP"].rt;
555  phi = -phi;
556  dx = bendings["BLC2LP"].x0 + bendings["BLC2LP"].rt * cos(-phi);
557  dz = bendings["BLC2LP"].z0 + bendings["BLC2LP"].rt * sin(-phi);
558  phi -= M_PI;
559  } else if (s < ler_breakpoints[6]) {
560  double sloc = s - ler_breakpoints[5];
561  phi = straights["LLL4"].phi;
562  dx = straights["LLL4"].x0 + sloc * sin(phi);
563  dz = straights["LLL4"].z0 + sloc * cos(phi);
564  }
565  // For this direction rotation angle of elements changes to negative,
566  // while SAD coordinates keep orientation.
567  // We need to compensate.
568  phi += M_PI;
569  }
570  // negative s
571  else if (s < -400.0 * Unit::cm) {
572  if (s > ler_breakpoints[7]) {
573  double sloc = -s;
574  phi = straights["LLR1"].phi;
575  dx = straights["LLR1"].x0 + sloc * sin(phi);
576  dz = straights["LLR1"].z0 + sloc * cos(phi);
577  } else if (s > ler_breakpoints[8]) {
578  double sloc = ler_breakpoints[7] - s;
579  phi = bendings["BC1RP"].sphi + bendings["BC1RP"].dphi - sloc / bendings["BC1RP"].rt;
580  phi = -phi;
581  dx = bendings["BC1RP"].x0 + bendings["BC1RP"].rt * cos(-phi);
582  dz = bendings["BC1RP"].z0 + bendings["BC1RP"].rt * sin(-phi);
583  phi += M_PI;
584  } else if (s > ler_breakpoints[9]) {
585  double sloc = ler_breakpoints[8] - s;
586  phi = straights["LLR2"].phi;
587  dx = straights["LLR2"].x0 + sloc * sin(phi);
588  dz = straights["LLR2"].z0 + sloc * cos(phi);
589  } else if (s > ler_breakpoints[10]) {
590  double sloc = ler_breakpoints[9] - s;
591  phi = bendings["BLCWRP"].sphi + bendings["BLCWRP"].dphi - sloc / bendings["BLCWRP"].rt;
592  phi = -phi;
593  dx = bendings["BLCWRP"].x0 + bendings["BLCWRP"].rt * cos(-phi);
594  dz = bendings["BLCWRP"].z0 + bendings["BLCWRP"].rt * sin(-phi);
595  phi += M_PI;
596  } else if (s > ler_breakpoints[11]) {
597  double sloc = ler_breakpoints[10] - s;
598  phi = straights["LLR3"].phi;
599  dx = straights["LLR3"].x0 + sloc * sin(phi);
600  dz = straights["LLR3"].z0 + sloc * cos(phi);
601  } else if (s > ler_breakpoints[12]) {
602  double sloc = ler_breakpoints[11] - s;
603  phi = bendings["BLC1RP"].sphi + bendings["BLC1RP"].dphi - sloc / bendings["BLC1RP"].rt;
604  phi = -phi;
605  dx = bendings["BLC1RP"].x0 + bendings["BLC1RP"].rt * cos(-phi);
606  dz = bendings["BLC1RP"].z0 + bendings["BLC1RP"].rt * sin(-phi);
607  phi += M_PI;
608  } else if (s > ler_breakpoints[13]) {
609  double sloc = ler_breakpoints[12] - s;
610  phi = straights["LLR4"].phi;
611  dx = straights["LLR4"].x0 + sloc * sin(phi);
612  dz = straights["LLR4"].z0 + sloc * cos(phi);
613  } else if (s > ler_breakpoints[14]) {
614  double sloc = ler_breakpoints[13] - s;
615  phi = bendings["BLC2RP"].sphi + sloc / bendings["BLC2RP"].rt;
616  phi = -phi;
617  dx = bendings["BLC2RP"].x0 + bendings["BLC2RP"].rt * cos(-phi);
618  dz = bendings["BLC2RP"].z0 + bendings["BLC2RP"].rt * sin(-phi);
619  } else if (s > ler_breakpoints[15]) {
620  double sloc = ler_breakpoints[14] - s;
621  phi = straights["LLR5"].phi;
622  dx = straights["LLR5"].x0 + sloc * sin(phi);
623  dz = straights["LLR5"].z0 + sloc * cos(phi);
624  }
625  }
626  }
627  if (accRing == c_HER) {
628  // HER
629  // positive s
630  if (400.0 * Unit::cm < s) {
631  if (s < her_breakpoints[0]) {
632  phi = straights["LHL1"].phi;
633  dx = straights["LHL1"].x0 + s * sin(phi);
634  dz = straights["LHL1"].z0 + s * cos(phi);
635  } else if (s < her_breakpoints[1]) {
636  double sloc = s - her_breakpoints[0];
637  phi = bendings["BLC1LE"].sphi + sloc / bendings["BLC1LE"].rt;
638  phi = -phi;
639  dx = bendings["BLC1LE"].x0 + bendings["BLC1LE"].rt * cos(-phi);
640  dz = bendings["BLC1LE"].z0 + bendings["BLC1LE"].rt * sin(-phi);
641  } else if (s < her_breakpoints[2]) {
642  double sloc = s - her_breakpoints[1];
643  phi = straights["LHL2"].phi;
644  dx = straights["LHL2"].x0 + sloc * sin(phi);
645  dz = straights["LHL2"].z0 + sloc * cos(phi);
646  }
647  phi += M_PI;
648  }
649  // negative s
650  else if (s < -400.0 * Unit::cm) {
651  if (s > her_breakpoints[3]) {
652  double sloc = -s;
653  phi = straights["LHR1"].phi;
654  dx = straights["LHR1"].x0 + sloc * sin(phi);
655  dz = straights["LHR1"].z0 + sloc * cos(phi);
656  } else if (s > her_breakpoints[4]) {
657  double sloc = her_breakpoints[3] - s;
658  phi = bendings["BLC2RE"].sphi + sloc / bendings["BLC2RE"].rt;
659  phi = -phi;
660  dx = bendings["BLC2RE"].x0 + bendings["BLC2RE"].rt * cos(-phi);
661  dz = bendings["BLC2RE"].z0 + bendings["BLC2RE"].rt * sin(-phi);
662  } else if (s > her_breakpoints[5]) {
663  double sloc = her_breakpoints[4] - s;
664  phi = straights["LHR2"].phi;
665  dx = straights["LHR2"].x0 + sloc * sin(phi);
666  dz = straights["LHR2"].z0 + sloc * cos(phi);
667  }
668  }
669  }
670 
671  TGeoHMatrix matrix("SADTrafo");
672  matrix.RotateY(phi / Unit::deg);
673  matrix.SetDx(dx);
674  matrix.SetDz(dz);
675  return matrix;
676 }
677 
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::ReaderSAD::m_lostY
double m_lostY
y at lost position [m].
Definition: ReaderSAD.h:147
Belle2::ReaderSAD::m_lostX
double m_lostX
x at lost position [m].
Definition: ReaderSAD.h:146
Belle2::Unit::s
static const double s
[second]
Definition: Unit.h:105
Belle2::ReaderSAD::m_inputSAD_py
double m_inputSAD_py
y momentum at lost position [m].
Definition: ReaderSAD.h:163
Belle2::ReaderSAD::m_inputSAD_rate
double m_inputSAD_rate
loss rate [Hz]
Definition: ReaderSAD.h:170
Belle2::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:57
Belle2::Gearbox::getDetectorComponent
GearDir getDetectorComponent(const std::string &component)
Return GearDir representing a given DetectorComponent.
Definition: Gearbox.cc:307
Belle2::ReaderSAD::m_accRing
AcceleratorRings m_accRing
The accelerator ring from which the particles originate.
Definition: ReaderSAD.h:135
Belle2::ReaderSAD::m_lostE
double m_lostE
energy at lost position [m].
Definition: ReaderSAD.h:152
Belle2::ReaderSAD::m_sRange
double m_sRange
The +- range for the s value for which particles are loaded.
Definition: ReaderSAD.h:134
Belle2::ReaderSAD::m_pxRes
double m_pxRes
The resolution for the x momentum component of the SAD real particle.
Definition: ReaderSAD.h:136
Belle2::ReaderSAD::m_lostS
double m_lostS
lost position [m] along ring.
Definition: ReaderSAD.h:148
Belle2::Gearbox::getInstance
static Gearbox & getInstance()
Return reference to the Gearbox instance.
Definition: Gearbox.cc:74
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::ReaderSAD::m_pyRes
double m_pyRes
The resolution for the y momentum component of the SAD real particle.
Definition: ReaderSAD.h:137
Belle2::ReaderSAD::~ReaderSAD
~ReaderSAD()
Destructor.
Definition: ReaderSAD.cc:65
Belle2::ReaderSAD::c_LER
@ c_LER
Low Energy Ring (positrons)
Definition: ReaderSAD.h:58
Belle2::ReaderSAD::c_HER
@ c_HER
High Energy Ring (electrons)
Definition: ReaderSAD.h:57
Belle2::SADMetaHit
ClassSADMetaHit - digitization simulated metahit for the SAD detector.
Definition: SADMetaHit.h:35
Belle2::ReaderSAD::straightElement::phi
double phi
Phi angle.
Definition: ReaderSAD.h:199
Belle2::ReaderSAD::getSADParticle
double getSADParticle(MCParticleGraph &graph)
Reads one SAD particle from the file and creates one event per SAD particle.
Definition: ReaderSAD.cc:123
Belle2::ReaderSAD::m_inputSAD_x
double m_inputSAD_x
x at lost position [m].
Definition: ReaderSAD.h:160
Belle2::ReaderSAD::m_tree
TTree * m_tree
The input root tree.
Definition: ReaderSAD.h:131
Belle2::ReaderSAD::m_inputSAD_nturn
int m_inputSAD_nturn
number of turns from scattered to lost
Definition: ReaderSAD.h:159
Belle2::ReaderSAD::open
void open(const std::string &filename)
Opens a root file and prepares it for reading.
Definition: ReaderSAD.cc:81
Belle2::ReaderSAD::m_readEntry
int m_readEntry
The current number of the SAD entry that is read.
Definition: ReaderSAD.h:144
Belle2::ReaderSAD::m_inputSAD_rr
double m_inputSAD_rr
sqrt(x*x+y*y) [m] before matching onto beam pipe inner surface
Definition: ReaderSAD.h:167
Belle2::ReaderSAD::m_realPartEntry
unsigned int m_realPartEntry
The current number of the created real particles.
Definition: ReaderSAD.h:143
Belle2::ReaderSAD::SADtoGeant
TGeoHMatrix SADtoGeant(ReaderSAD::AcceleratorRings accRing, double s)
Transformation matrix.
Definition: ReaderSAD.cc:417
Belle2::ReaderSAD::getRealParticle
bool getRealParticle(MCParticleGraph &graph)
Reads one SAD particle from the file, calculates the number of real particles which are represented b...
Definition: ReaderSAD.cc:166
Belle2::ReaderSAD::m_inputSAD_watt
double m_inputSAD_watt
loss wattage [W]
Definition: ReaderSAD.h:171
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ReaderSAD::m_readoutTime
double m_readoutTime
The readout time.
Definition: ReaderSAD.h:140
Belle2::ReaderSAD::m_inputSAD_sraw
double m_inputSAD_sraw
lost position [m
Definition: ReaderSAD.h:155
Belle2::ReaderSAD::bendingElement::rt
double rt
Parameter.
Definition: ReaderSAD.h:204
Belle2::ReaderSAD::addParticleToMCParticles
void addParticleToMCParticles(MCParticleGraph &graph, bool gaussSmearing=false)
Adds the current particle described by the member variables to the MCParticles collection.
Definition: ReaderSAD.cc:267
Belle2::ReaderSAD::m_inputSAD_px
double m_inputSAD_px
x momentum at lost position [m].
Definition: ReaderSAD.h:162
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::ReaderSAD::convertParamsToSADUnits
void convertParamsToSADUnits()
Convert the parameters from the SAD units to the basf2 units.
Definition: ReaderSAD.cc:237
Belle2::ReaderSAD::bendingElement::z0
double z0
Initial position in Z
Definition: ReaderSAD.h:206
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::ReaderSAD::m_inputSAD_ss
double m_inputSAD_ss
scattered position (|s|<Ltot/2) [m]
Definition: ReaderSAD.h:156
Belle2::ReaderSAD::m_transMatrix
TGeoHMatrix * m_transMatrix
Transformation matrix from local SAD to global geant4 space.
Definition: ReaderSAD.h:133
Belle2::ReaderSAD
Class to read files that have been created by SAD and store their content in a MCParticle graph.
Definition: ReaderSAD.h:45
Belle2::ReaderSAD::bendingElement
Sensitive Element
Definition: ReaderSAD.h:203
Belle2::ReaderSAD::m_inputSAD_y
double m_inputSAD_y
y at lost position [m].
Definition: ReaderSAD.h:161
Belle2::ReaderSAD::initialize
void initialize(TGeoHMatrix *transMatrix, double sRange, AcceleratorRings accRing, double readoutTime)
Initializes the reader, sets the particle parameters and calculates important values.
Definition: ReaderSAD.cc:72
Belle2::Unit::deg
static const double deg
degree to radians
Definition: Unit.h:119
Belle2::Unit::m
static const double m
[meters]
Definition: Unit.h:79
Belle2::ReaderSAD::m_lostPy
double m_lostPy
y momentum at lost position [m].
Definition: ReaderSAD.h:150
Belle2::ReaderSAD::AcceleratorRings
AcceleratorRings
The both accelerator rings.
Definition: ReaderSAD.h:56
Belle2::ReaderSAD::bendingElement::x0
double x0
Initial position in X
Definition: ReaderSAD.h:205
Belle2::ReaderSAD::m_inputSAD_Lss
double m_inputSAD_Lss
length of element in which scattered [m]
Definition: ReaderSAD.h:158
Belle2::ReaderSAD::m_file
TFile * m_file
The input root file.
Definition: ReaderSAD.h:130
Belle2::ReaderSAD::bendingElement::dphi
double dphi
Bending Angle for torus in phi.
Definition: ReaderSAD.h:208
Belle2::ReaderSAD::m_inputSAD_r
double m_inputSAD_r
sqrt(x*x+y*y) [m]
Definition: ReaderSAD.h:166
Belle2::ReaderSAD::m_inputSAD_dp_over_p0
double m_inputSAD_dp_over_p0
dp_over_p0
Definition: ReaderSAD.h:168
Belle2::ReaderSAD::m_inputSAD_s
double m_inputSAD_s
lost position (|s|<Ltot/2) [m]
Definition: ReaderSAD.h:157
Belle2::ReaderSAD::calculateRealParticleNumber
int calculateRealParticleNumber(double rate)
Calculates the number of real particles for a SAD particle.
Definition: ReaderSAD.cc:404
Belle2::ReaderSAD::bendingElement::sphi
double sphi
Bending Angle.
Definition: ReaderSAD.h:207
Belle2::ReaderSAD::m_inputSAD_xraw
double m_inputSAD_xraw
x at lost position [m] before matching onto beam pipe inner surface
Definition: ReaderSAD.h:164
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::ReaderSAD::m_inputSAD_ssraw
double m_inputSAD_ssraw
scattered position [m]
Definition: ReaderSAD.h:154
Belle2::ReaderSAD::m_inputSAD_E
double m_inputSAD_E
energy at lost position [m].
Definition: ReaderSAD.h:169
Belle2::ReaderSAD::straightElement::l
double l
Parameter.
Definition: ReaderSAD.h:198
Belle2::ReaderSAD::straightElement
Calculates the transformation matrix from local SAD to global geant4 space.
Definition: ReaderSAD.h:195
Belle2::Unit::GeV
static const double GeV
Standard of [energy, momentum, mass].
Definition: Unit.h:61
Belle2::ReaderSAD::addAllSADParticles
void addAllSADParticles(MCParticleGraph &graph)
Reads all SAD particles from the file into the MCParticles collection which are inside the specified ...
Definition: ReaderSAD.cc:213
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::ReaderSAD::straightElement::x0
double x0
Initial position in X.
Definition: ReaderSAD.h:196
Belle2::ReaderSAD::m_lostPx
double m_lostPx
x momentum at lost position [m].
Definition: ReaderSAD.h:149
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58
Belle2::ReaderSAD::m_lostRate
double m_lostRate
loss rate [Hz]>
Definition: ReaderSAD.h:151
Belle2::ReaderSAD::m_inputSAD_yraw
double m_inputSAD_yraw
y at lost position [m] before matching onto beam pipe inner surface
Definition: ReaderSAD.h:165
Belle2::ReaderSAD::m_realPartNum
unsigned int m_realPartNum
The current number of the created real particles.
Definition: ReaderSAD.h:142
Belle2::ReaderSAD::straightElement::z0
double z0
Initial position in Z.
Definition: ReaderSAD.h:197