11 #include <generators/SAD/ReaderSAD.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Unit.h>
15 #include <framework/gearbox/GearDir.h>
16 #include <mdst/dataobjects/MCParticle.h>
18 #include <generators/SAD/dataobjects/SADMetaHit.h>
21 #include <framework/datastore/StoreArray.h>
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)
88 m_file =
new TFile(filename.c_str(),
"READ");
89 if (
m_file == NULL)
throw(SADCouldNotOpenFileError() << filename);
93 if (
m_tree == NULL)
throw(SADCouldNotOpenFileError() << filename);
118 hit.registerInDataStore();
126 B2ERROR(
"The SAD tree doesn't exist !");
139 B2DEBUG(10,
"> Read particle " <<
m_readEntry + 1 <<
"/" <<
m_tree->GetEntries() <<
" with s = " <<
m_lostS <<
" cm" <<
145 if (zMom2 < 0) printf(
"zMom2= %f is negative. Skipped!\n", zMom2);
169 B2ERROR(
"The SAD tree doesn't exist !");
192 B2DEBUG(10,
"> Read particle " <<
m_readEntry + 1 <<
"/" <<
m_tree->GetEntries() <<
" with s = " <<
m_lostS <<
" cm" <<
194 }
while ((fabs(
m_lostS) >
m_sRange) && (m_readEntry < m_tree->GetEntries()));
204 "/" <<
m_tree->GetEntries());
216 B2ERROR(
"The SAD tree doesn't exist !");
220 int nPart =
m_tree->GetEntries();
222 for (
int iPart = 0; iPart < nPart; ++iPart) {
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};
280 case c_HER: particle.setPDG(11);
282 case c_LER: particle.setPDG(-11);
286 particle.setMassFromPDG();
294 particlePosSADfar[0] =
m_lostX;
295 particlePosSADfar[1] = -
m_lostY;
296 particlePosSADfar[2] = 0;
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.");
306 m_transMatrix->LocalToMaster(particlePosSAD, particlePosGeant4);
308 m_transMatrix2->LocalToMaster(particlePosSADfar, particlePosGeant4);
325 double zMom = sqrt(totalMomSqr - (particleMomSAD[0] * particleMomSAD[0]) - (particleMomSAD[1] * particleMomSAD[1]));
328 case c_HER: particleMomSAD[2] = zMom;
330 case c_LER: particleMomSAD[2] = -zMom;
337 case c_HER: ring = 1;
339 case c_LER: ring = 2;
366 int section_ordering[12] = {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2};
372 }
else if (ring == 2) {
378 int ring_section = section_ordering[(int)((ssraw) / 250.)];
381 m_transMatrix->LocalToMasterVect(particleMomSAD, particleMomGeant4);
383 m_transMatrix2->LocalToMasterVect(particleMomSAD, particleMomGeant4);
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);
408 int numPart_int =
static_cast<int>(floor(numPart));
409 double numPart_dec = numPart - numPart_int;
411 double rnd = gRandom->Uniform();
412 if (rnd < numPart_dec) numPart_int += 1;
428 map<string, straightElement> straights;
429 map<string, bendingElement> bendings;
430 for (
const GearDir& element : content.getNodes(
"Straight")) {
432 string name = element.getString(
"@name");
433 string type = element.getString(
"@type");
435 if (type !=
"pipe")
continue;
439 straight.
x0 = element.getLength(
"X0");
440 straight.
z0 = element.getLength(
"Z0");
441 straight.
l = element.getLength(
"L");
442 straight.
phi = element.getLength(
"PHI");
444 straights[name] = straight;
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.");
453 for (
const GearDir& element : content.getNodes(
"Bending")) {
455 string name = element.getString(
"@name");
456 string type = element.getString(
"@type");
458 if (type !=
"pipe")
continue;
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");
468 bendings[name] = bending;
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.");
477 static double her_breakpoints[6];
478 static double ler_breakpoints[16];
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;
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;
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;
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;
513 if (accRing ==
c_LER) {
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;
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;
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];
554 phi = bendings[
"BLC2LP"].sphi + bendings[
"BLC2LP"].dphi - sloc / bendings[
"BLC2LP"].rt;
556 dx = bendings[
"BLC2LP"].x0 + bendings[
"BLC2LP"].rt * cos(-phi);
557 dz = bendings[
"BLC2LP"].z0 + bendings[
"BLC2LP"].rt * sin(-phi);
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);
572 if (s > ler_breakpoints[7]) {
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;
581 dx = bendings[
"BC1RP"].x0 + bendings[
"BC1RP"].rt * cos(-phi);
582 dz = bendings[
"BC1RP"].z0 + bendings[
"BC1RP"].rt * sin(-phi);
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;
593 dx = bendings[
"BLCWRP"].x0 + bendings[
"BLCWRP"].rt * cos(-phi);
594 dz = bendings[
"BLCWRP"].z0 + bendings[
"BLCWRP"].rt * sin(-phi);
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;
605 dx = bendings[
"BLC1RP"].x0 + bendings[
"BLC1RP"].rt * cos(-phi);
606 dz = bendings[
"BLC1RP"].z0 + bendings[
"BLC1RP"].rt * sin(-phi);
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;
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);
627 if (accRing ==
c_HER) {
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;
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);
651 if (s > her_breakpoints[3]) {
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;
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);
671 TGeoHMatrix matrix(
"SADTrafo");