diff --git a/Inputs/DetectorConfiguration/Vamos.detector b/Inputs/DetectorConfiguration/Vamos.detector new file mode 100644 index 0000000000000000000000000000000000000000..188e5613668a2c99fb99518b73ec6f4ee70e2f5b --- /dev/null +++ b/Inputs/DetectorConfiguration/Vamos.detector @@ -0,0 +1,98 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Target + THICKNESS= 5 micrometer + RADIUS= 20 mm + MATERIAL= CD2 + ANGLE= 0 deg + X= 0 mm + Y= 0 mm + Z= 0 mm + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Vamos Configuration +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Vamos + R= 1 m + Theta= 0 deg + +Vamos DC + Z= 60 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos DC + Z= 200 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos BeamCatcher + Material= BC400 + Width= 30 mm + Length= 80 mm + Thickness= 30 mm + Pos= 0 0 600 mm + +Vamos MWPPAC + Z= 750 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos DC + Z= 940 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + +Vamos DC + Z= 1060 mm + Gas= iC4H10 + Pressure= 0.01 bar + Temperature= 295 kelvin + + +Vamos IC + Z= 1200 mm + Thickness= 50 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1250 mm + Thickness= 50 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1300 mm + Thickness= 50 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1375 mm + Thickness= 100 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1525 mm + Thickness= 200 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin + +Vamos IC + Z= 1725 mm + Thickness= 200 mm + Gas= CF4 + Pressure= 0.2 bar + Temperature= 295 kelvin +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx index 14d18e037da057d2ce4c55d77bed100503d20b65..f990f16e2af1daf53500025a093d541f07b8516b 100644 --- a/NPLib/Detectors/MUST2/TMust2Physics.cxx +++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx @@ -42,6 +42,7 @@ using namespace NPUNITS; /////////////////////////////////////////////////////////////////////////// ClassImp(TMust2Physics) + /////////////////////////////////////////////////////////////////////////// TMust2Physics::TMust2Physics() { EventMultiplicity = 0; @@ -70,6 +71,10 @@ ClassImp(TMust2Physics) m_Take_E_Y = false; m_Take_T_Y = true; + ////////////////// + // SiLi matching // + ////////////////// + m_SiLi_Size = 32; m_SiLi_MatchingX.resize(16, 0); m_SiLi_MatchingY.resize(16, 0); @@ -122,9 +127,59 @@ ClassImp(TMust2Physics) m_SiLi_MatchingX[15] = 16; m_SiLi_MatchingY[15] = 16; + ////////////////// + // CsI matching // + ////////////////// + m_CsI_Size = 32; m_CsI_MatchingX.resize(16, 0); m_CsI_MatchingY.resize(16, 0); + + m_CsI_MatchingX[0] = 112; + m_CsI_MatchingY[0] = 112; + + m_CsI_MatchingX[1] = 112; + m_CsI_MatchingY[1] = 80; + + m_CsI_MatchingX[2] = 112; + m_CsI_MatchingY[2] = 48; + + m_CsI_MatchingX[3] = 112; + m_CsI_MatchingY[3] = 16; + // + m_CsI_MatchingX[4] = 80; + m_CsI_MatchingY[4] = 16; + + m_CsI_MatchingX[5] = 80; + m_CsI_MatchingY[5] = 48; + + m_CsI_MatchingX[6] = 80; + m_CsI_MatchingY[6] = 80; + + m_CsI_MatchingX[7] = 80; + m_CsI_MatchingY[7] = 112; + // + m_CsI_MatchingX[8] = 48; + m_CsI_MatchingY[8] = 16; + + m_CsI_MatchingX[9] = 48; + m_CsI_MatchingY[9] = 48; + + m_CsI_MatchingX[10] = 48; + m_CsI_MatchingY[10] = 80; + + m_CsI_MatchingX[11] = 48; + m_CsI_MatchingY[11] = 112; + // + m_CsI_MatchingX[12] = 16; + m_CsI_MatchingY[12] = 16; + + m_CsI_MatchingX[13] = 16; + m_CsI_MatchingY[13] = 48; + + m_CsI_MatchingX[14] = 16; + m_CsI_MatchingY[14] = 80; + m_CsI_MatchingX[0] = 112; m_CsI_MatchingY[0] = 112; @@ -185,6 +240,9 @@ void TMust2Physics::BuildPhysicalEvent() { PreTreat(); + // FIXME has to be set in the configuration file + bool InterstripTreatment = false; + bool check_SILI = false; bool check_CSI = false; @@ -197,48 +255,120 @@ void TMust2Physics::BuildPhysicalEvent() { m_CsIEMult = m_PreTreatedData->GetMMCsIEMult(); m_CsITMult = m_PreTreatedData->GetMMCsITMult(); - if (1 /*CheckEvent() == 1*/) { - vector<TVector2> couple = Match_X_Y(); + ///////////////////////////////////////////////// + vector<TVector2> couple = Match_X_Y(); + ///////////////////////////////////////////////// + + EventMultiplicity = std::min(m_StripXEMult, m_StripYEMult); + + unsigned int couple_size = couple.size(); - EventMultiplicity = couple.size(); - for (unsigned int i = 0; i < couple.size(); ++i) { - check_SILI = false; - check_CSI = false; + for (unsigned int i = 0; i < couple_size; ++i) { - int N = m_PreTreatedData->GetMMStripXEDetectorNbr(couple[i].X()); + check_SILI = false; + check_CSI = false; - int X = m_PreTreatedData->GetMMStripXEStripNbr(couple[i].X()); - int Y = m_PreTreatedData->GetMMStripYEStripNbr(couple[i].Y()); + int N = m_PreTreatedData->GetMMStripXEDetectorNbr(couple[i].X()); + + int X = m_PreTreatedData->GetMMStripXEStripNbr(couple[i].X()); + int Y = m_PreTreatedData->GetMMStripYEStripNbr(couple[i].Y()); + + double Si_X_E = m_PreTreatedData->GetMMStripXEEnergy(couple[i].X()); + double Si_Y_E = m_PreTreatedData->GetMMStripYEEnergy(couple[i].Y()); + + // Search for associate Time + double Si_X_T = -1000; + for (unsigned int t = 0; t < m_StripXTMult; ++t) { + if (m_PreTreatedData->GetMMStripXTStripNbr(couple[i].X()) + == m_PreTreatedData->GetMMStripXTStripNbr(t)) { + Si_X_T = m_PreTreatedData->GetMMStripXTTime(t); + break; + } + } - double Si_X_E = m_PreTreatedData->GetMMStripXEEnergy(couple[i].X()); - double Si_Y_E = m_PreTreatedData->GetMMStripYEEnergy(couple[i].Y()); + double Si_Y_T = -1000; + for (unsigned int t = 0; t < m_StripYTMult; ++t) { + if (m_PreTreatedData->GetMMStripYTStripNbr(couple[i].Y()) + == m_PreTreatedData->GetMMStripYTStripNbr(t)) { + Si_Y_T = m_PreTreatedData->GetMMStripYTTime(t); + break; + } + } - // Search for associate Time - double Si_X_T = -1000; - for (unsigned int t = 0; t < m_StripXTMult; ++t) { - if (m_PreTreatedData->GetMMStripXTStripNbr(couple[i].X()) - == m_PreTreatedData->GetMMStripXTStripNbr(t) - && m_PreTreatedData->GetMMStripXTDetectorNbr(couple[i].X()) - == m_PreTreatedData->GetMMStripXTDetectorNbr(t)) { - Si_X_T = m_PreTreatedData->GetMMStripXTTime(t); - break; + for (unsigned int j = 0; j < m_SiLiEMult; ++j) { + if (m_PreTreatedData->GetMMSiLiEDetectorNbr(j) == N) { + // pad vs strip number match + if (Match_Si_SiLi(X, Y, m_PreTreatedData->GetMMSiLiEPadNbr(j))) { + SiLi_N.push_back(m_PreTreatedData->GetMMSiLiEPadNbr(j)); + SiLi_E.push_back(m_PreTreatedData->GetMMSiLiEEnergy(j)); + SiLi_T.push_back(-1000); + // Look for associate time + for (unsigned int k = 0; k < m_SiLiTMult; ++k) { + // Same Pad, same Detector + if (m_PreTreatedData->GetMMSiLiEPadNbr(j) + == m_PreTreatedData->GetMMSiLiTPadNbr(k)) { + SiLi_T[SiLi_T.size() - 1] = m_PreTreatedData->GetMMSiLiTTime(k); + break; + } + } + check_SILI = true; } } + } - double Si_Y_T = -1000; - for (unsigned int t = 0; t < m_StripYTMult; ++t) { - if (m_PreTreatedData->GetMMStripYTStripNbr(couple[i].Y()) - == m_PreTreatedData->GetMMStripYTStripNbr(t) - && m_PreTreatedData->GetMMStripYTDetectorNbr(couple[i].Y()) - == m_PreTreatedData->GetMMStripYTDetectorNbr(t)) { - Si_Y_T = m_PreTreatedData->GetMMStripYTTime(t); - break; + for (unsigned int j = 0; j < m_CsIEMult; ++j) { + if (m_PreTreatedData->GetMMCsIEDetectorNbr(j) == N) { + if (Match_Si_CsI(X, Y, m_PreTreatedData->GetMMCsIECristalNbr(j))) { + CsI_N.push_back(m_PreTreatedData->GetMMCsIECristalNbr(j)); + CsI_E.push_back(m_PreTreatedData->GetMMCsIEEnergy(j)); + CsI_T.push_back(-1000); + // Look for associate Time + for (unsigned int k = 0; k < m_CsITMult; ++k) { + // Same Cristal, Same Detector + if (m_PreTreatedData->GetMMCsIECristalNbr(j) + == m_PreTreatedData->GetMMCsITCristalNbr(k)) { + CsI_T[CsI_T.size() - 1] = m_PreTreatedData->GetMMCsITTime(j); + break; + } + } + check_CSI = true; } } + } - Si_X.push_back(X); - Si_Y.push_back(Y); - TelescopeNumber.push_back(N); + TelescopeNumber.push_back(N); + Si_X.push_back(X); + Si_Y.push_back(Y); + + // Store the other value for checking purpose + Si_EX.push_back(Si_X_E); + Si_TX.push_back(Si_X_T); + + Si_EY.push_back(Si_Y_E); + Si_TY.push_back(Si_Y_T); + + ///////////////////////////////////////////////// + CheckEvent(N); + ///////////////////////////////////////////////// + + if (m_OrderMatch == 2) { + Clear(); + } else { + // This test retrieves only Y interstrips + // see Xavier Mougeot's PhD. In the case where + // if (InterstripTreatment == true) { + // EventType.push_back(2); + // Si_E.push_back(std::max(Si_X_E, Si_Y_E)); + // } else if (InterstripTreatment == false) { + // EventType.push_back(2); + // Si_E.push_back(-1000); + // } else { + // EventType.push_back(1); + // if (m_Take_E_Y) + // Si_E.push_back(Si_Y_E); + // else + // Si_E.push_back(Si_X_E); + // } if (m_Take_E_Y) Si_E.push_back(Si_Y_E); @@ -250,60 +380,6 @@ void TMust2Physics::BuildPhysicalEvent() { else Si_T.push_back(Si_X_T); - // Store the other value for checking purpose - Si_EX.push_back(Si_X_E); - Si_TX.push_back(Si_X_T); - - Si_EY.push_back(Si_Y_E); - Si_TY.push_back(Si_Y_T); - - for (unsigned int j = 0; j < m_SiLiEMult; ++j) { - if (m_PreTreatedData->GetMMSiLiEDetectorNbr(j) == N) { - // pad vs strip number match - if (Match_Si_SiLi(X, Y, m_PreTreatedData->GetMMSiLiEPadNbr(j))) { - SiLi_N.push_back(m_PreTreatedData->GetMMSiLiEPadNbr(j)); - SiLi_E.push_back(m_PreTreatedData->GetMMSiLiEEnergy(j)); - SiLi_T.push_back(-1000); - // Look for associate time - for (unsigned int k = 0; k < m_SiLiTMult; ++k) { - // Same Pad, same Detector - if (m_PreTreatedData->GetMMSiLiEPadNbr(j) - == m_PreTreatedData->GetMMSiLiTPadNbr(k) - && m_PreTreatedData->GetMMSiLiEDetectorNbr(j) - == m_PreTreatedData->GetMMSiLiTDetectorNbr(k)) { - SiLi_T[SiLi_T.size() - 1] = m_PreTreatedData->GetMMSiLiTTime(k); - break; - } - } - - check_SILI = true; - } - } - } - - for (unsigned int j = 0; j < m_CsIEMult; ++j) { - if (m_PreTreatedData->GetMMCsIEDetectorNbr(j) == N) { - if (Match_Si_CsI(X, Y, m_PreTreatedData->GetMMCsIECristalNbr(j))) { - CsI_N.push_back(m_PreTreatedData->GetMMCsIECristalNbr(j)); - CsI_E.push_back(m_PreTreatedData->GetMMCsIEEnergy(j)); - CsI_T.push_back(-1000); - // Look for associate Time - for (unsigned int k = 0; k < m_CsITMult; ++k) { - // Same Cristal, Same Detector - if (m_PreTreatedData->GetMMCsIECristalNbr(j) - == m_PreTreatedData->GetMMCsITCristalNbr(k) - && m_PreTreatedData->GetMMCsIEDetectorNbr(j) - == m_PreTreatedData->GetMMCsITDetectorNbr(k)) { - CsI_T[CsI_T.size() - 1] = m_PreTreatedData->GetMMCsITTime(j); - break; - } - } - - check_CSI = true; - } - } - } - if (!check_SILI) { SiLi_N.push_back(0); SiLi_E.push_back(-1000); @@ -316,7 +392,9 @@ void TMust2Physics::BuildPhysicalEvent() { CsI_T.push_back(-1000); } } - } + } // loop on event multiplicity + + // Check CsI_Size return; } @@ -333,12 +411,10 @@ void TMust2Physics::PreTreat() { m_CsITMult = m_EventData->GetMMCsITMult(); map<int, int> hitX; map<int, int> hitY; + // X // E for (unsigned int i = 0; i < m_StripXEMult; ++i) { - if (m_EventData->GetMMStripXEDetectorNbr(i) == 5) - hitX[m_EventData->GetMMStripXEStripNbr(i)]++; - if (m_EventData->GetMMStripXEEnergy(i) > m_Si_X_E_RAW_Threshold && IsValidChannel(0, m_EventData->GetMMStripXEDetectorNbr(i), m_EventData->GetMMStripXEStripNbr(i))) { @@ -348,12 +424,7 @@ void TMust2Physics::PreTreat() { m_EventData->GetMMStripXEStripNbr(i), EX); } } - map<int, int>::iterator it = hitX.begin(); - for (; it != hitX.end(); it++) { - if (it->second > 1) - cout << "Raw Strip X " << it->first << " hit " << it->second << " times" - << endl; - } + // T for (unsigned int i = 0; i < m_StripXTMult; ++i) { if (IsValidChannel(0, m_EventData->GetMMStripXTDetectorNbr(i), @@ -366,9 +437,6 @@ void TMust2Physics::PreTreat() { // Y // E for (unsigned int i = 0; i < m_StripYEMult; ++i) { - if (m_EventData->GetMMStripYEDetectorNbr(i) == 5) - hitY[m_EventData->GetMMStripYEStripNbr(i)]++; - if (m_EventData->GetMMStripYEEnergy(i) < m_Si_Y_E_RAW_Threshold && IsValidChannel(1, m_EventData->GetMMStripYEDetectorNbr(i), m_EventData->GetMMStripYEStripNbr(i))) { @@ -379,13 +447,6 @@ void TMust2Physics::PreTreat() { } } - map<int, int>::iterator it2 = hitY.begin(); - for (; it2 != hitY.end(); it2++) { - if (it2->second > 1) - cout << "Raw Strip Y " << it2->first << " hit " << it2->second << " times" - << endl; - } - // T for (unsigned int i = 0; i < m_StripYTMult; ++i) { if (IsValidChannel(1, m_EventData->GetMMStripYTDetectorNbr(i), @@ -442,33 +503,21 @@ void TMust2Physics::PreTreat() { return; } -/////////////////////////////////////////////////////////////////////////// -int TMust2Physics::CheckEvent() { - // Check the size of the different elements - if (m_PreTreatedData->GetMMStripXEMult() - == m_PreTreatedData->GetMMStripYEMult()) - return 1; // Regular Event - - // INterstrip management is not coded, so waste of time to make this test - /* else if( m_PreTreatedData->GetMMStripXEMult() == - m_PreTreatedData->GetMMStripYEMult()+1 - || m_PreTreatedData->GetMMStripXEMult() == - m_PreTreatedData->GetMMStripYEMult()-1 ) - return 2 ; // Pseudo Event, potentially interstrip*/ - - else - return -1; // Rejected Event -} - /////////////////////////////////////////////////////////////////////////// bool TMust2Physics::ResolvePseudoEvent() { return false; } /////////////////////////////////////////////////////////////////////////// + vector<TVector2> TMust2Physics::Match_X_Y() { vector<TVector2> ArrayOfGoodCouple; + ArrayOfGoodCouple.clear(); + m_StripXEMult = m_PreTreatedData->GetMMStripXEMult(); m_StripYEMult = m_PreTreatedData->GetMMStripYEMult(); + double matchSigma = m_StripEnergyMatchingSigma; + double NmatchSigma = m_StripEnergyMatchingNumberOfSigma; + // Prevent code from treating very high multiplicity Event // Those event are not physical anyway and that improve speed. if (m_StripXEMult > m_MaximumStripMultiplicityAllowed @@ -476,25 +525,42 @@ vector<TVector2> TMust2Physics::Match_X_Y() { return ArrayOfGoodCouple; } + // Get Detector multiplicity + // FIXME Variable names + for (unsigned int i = 0; i < m_StripXEMult; i++) { + int N = m_PreTreatedData->GetMMStripXEDetectorNbr(i); + m_StripXMultDet[N] += 1; + } + + for (unsigned int j = 0; j < m_StripYEMult; j++) { + int N = m_PreTreatedData->GetMMStripYEDetectorNbr(j); + m_StripYMultDet[N] += 1; + } + for (unsigned int i = 0; i < m_StripXEMult; i++) { for (unsigned int j = 0; j < m_StripYEMult; j++) { // Declaration of variable for clarity - double StripXDetNbr = m_PreTreatedData->GetMMStripXEDetectorNbr(i); - double StripYDetNbr = m_PreTreatedData->GetMMStripYEDetectorNbr(j); + int StripXDetNbr = m_PreTreatedData->GetMMStripXEDetectorNbr(i); + int StripYDetNbr = m_PreTreatedData->GetMMStripYEDetectorNbr(j); // if same detector check energy if (StripXDetNbr == StripYDetNbr) { + int DetNbr = StripXDetNbr; + // Declaration of variable for clarity double StripXEnergy = m_PreTreatedData->GetMMStripXEEnergy(i); double StripXNbr = m_PreTreatedData->GetMMStripXEStripNbr(i); + double StripYEnergy = m_PreTreatedData->GetMMStripYEEnergy(j); double StripYNbr = m_PreTreatedData->GetMMStripYEStripNbr(j); - double matchSigma = m_StripEnergyMatchingSigma; - double NmatchSigma = m_StripEnergyMatchingNumberOfSigma; // Look if energy match + // FIXME Should be proportional to the energy loss in the DSSDs + // if (abs(StripXEnergy - StripYEnergy) + // < 0.09 * (std::max(StripXEnergy, StripYEnergy))) { + // negligible correction according to Adrien if (abs((StripXEnergy - StripYEnergy) / 2.) < NmatchSigma * matchSigma) { @@ -512,54 +578,38 @@ vector<TVector2> TMust2Physics::Match_X_Y() { } } - // Special Option, if the event is between two SiLi pad , it is - // rejected. + // Special Option, if the event is between two SiLi pad , + // it is rejected. else if (m_Ignore_not_matching_SiLi) { bool check_validity = false; for (unsigned int hh = 0; hh < 16; ++hh) { if (Match_Si_SiLi(StripXNbr, StripYNbr, hh + 1)) check_validity = true; } - if (check_validity) ArrayOfGoodCouple.push_back(TVector2(i, j)); } - // Regular case, keep the event else { - // Gives a unique ID for every telescope and strip combination - int IDX = m_NumberOfTelescope * StripXNbr + StripXDetNbr; - int IDY = m_NumberOfTelescope * StripYNbr + StripYDetNbr; - - m_HitStripX[IDX]++; - m_HitStripY[IDY]++; - ArrayOfGoodCouple.push_back(TVector2(i, j)); + m_NMatchDet[DetNbr] += 1; } } - } - } - } + } // if same detector + } // loop on StripY Mult + } // loop on StripX Mult - // Prevent to treat event with ambiguous matching beetween X and Y - map<int, int>::iterator itX = m_HitStripX.begin(); - for (; itX != m_HitStripX.end(); itX++) { - if (itX->second > 1) { - ArrayOfGoodCouple.clear(); - } - } + return ArrayOfGoodCouple; +} - map<int, int>::iterator itY = m_HitStripY.begin(); - for (; itY != m_HitStripY.end(); itY++) { - if (itY->second > 1) { - ArrayOfGoodCouple.clear(); - } +//////////////////////////////////////////////////////////////////////////// +void TMust2Physics::CheckEvent(int N) { + if (m_NMatchDet[N] > m_StripXMultDet[N] + || m_NMatchDet[N] > m_StripYMultDet[N]) { + m_OrderMatch = 2; + } else { + m_OrderMatch = 1; } - - m_HitStripX.clear(); - m_HitStripY.clear(); - - return ArrayOfGoodCouple; } //////////////////////////////////////////////////////////////////////////// @@ -593,7 +643,8 @@ void TMust2Physics::ReadAnalysisConfig() { AnalysisConfigFile.open(FileName.c_str()); if (!AnalysisConfigFile.is_open()) { - cout << " No ConfigMust2.dat found: Default parameters loaded for Analysis " + cout << " No ConfigMust2.dat found: Default parameters loaded for " + "Analysis " << FileName << endl; return; } @@ -806,8 +857,9 @@ bool TMust2Physics::Match_Si_SiLi(int X, int Y, int PadNbr) { bool TMust2Physics::Match_Si_CsI(int X, int Y, int CristalNbr) { if (abs(m_CsI_MatchingX[CristalNbr - 1] - X) < (double)m_CsI_Size / 2. - && abs(m_CsI_MatchingY[CristalNbr - 1] - Y) < (double)m_CsI_Size / 2.) + && abs(m_CsI_MatchingY[CristalNbr - 1] - Y) < (double)m_CsI_Size / 2.) { return true; + } else return false; @@ -817,6 +869,12 @@ bool TMust2Physics::Match_Si_CsI(int X, int Y, int CristalNbr) { void TMust2Physics::Clear() { EventMultiplicity = 0; + m_OrderMatch = 0; + + m_StripXMultDet.clear(); + m_StripYMultDet.clear(); + m_NMatchDet.clear(); + TelescopeNumber.clear(); EventType.clear(); TotalEnergy.clear(); @@ -920,6 +978,8 @@ void TMust2Physics::ReadConfiguration(NPL::InputParser parser) { TVector3 C = blocks[i]->GetTVector3("X1_Y128", "mm"); TVector3 D = blocks[i]->GetTVector3("X128_Y128", "mm"); AddTelescope(A, B, C, D); + m_CsIPresent[i + 1] = blocks[i]->GetInt("CSI"); + m_SiLiPresent[i + 1] = blocks[i]->GetInt("SILI"); } else if (blocks[i]->HasTokenList(sphe)) { @@ -931,16 +991,22 @@ void TMust2Physics::ReadConfiguration(NPL::InputParser parser) { double R = blocks[i]->GetDouble("R", "mm"); vector<double> beta = blocks[i]->GetVectorDouble("BETA", "deg"); AddTelescope(Theta, Phi, R, beta[0], beta[1], beta[2]); + + m_CsIPresent[i + 1] = blocks[i]->GetInt("CSI"); + m_SiLiPresent[i + 1] = blocks[i]->GetInt("SILI"); } else { - cout << "ERROR: Missing token for M2Telescope blocks, check your input " + cout << "ERROR: Missing token for M2Telescope blocks, check your " + "input " "file" << endl; exit(1); } } + // m_MultEvt = new TMust2MultTelescope[m_NumberOfTelescope]; + InitializeStandardParameter(); ReadAnalysisConfig(); } @@ -1085,8 +1151,8 @@ void TMust2Physics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1, m_NumberOfTelescope++; - // Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y - // strip are allong X axis) + // Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that + // Y strip are allong X axis) TVector3 U = C_X128_Y1 - C_X1_Y1; double Ushift = (U.Mag() - 98) / 2.; U = U.Unit(); @@ -1376,21 +1442,21 @@ double fCsI_T(const TMust2Data* m_EventData, const int& i) { return CalibrationManager::getInstance()->ApplyCalibration( name, m_EventData->GetMMCsITTime(i)); } -} +} // namespace MUST2_LOCAL //////////////////////////////////////////////////////////////////////////////// -// Construct Method to be pass to the DetectorFactory // +// Construct Method to be pass to the DetectorFactory // //////////////////////////////////////////////////////////////////////////////// NPL::VDetector* TMust2Physics::Construct() { return (NPL::VDetector*)new TMust2Physics(); } //////////////////////////////////////////////////////////////////////////////// -// Registering the construct method to the factory // +// Registering the construct method to the factory // //////////////////////////////////////////////////////////////////////////////// extern "C" { class proxy_must2 { - public: +public: proxy_must2() { NPL::DetectorFactory::getInstance()->AddToken("M2Telescope", "MUST2"); NPL::DetectorFactory::getInstance()->AddDetector("M2Telescope", diff --git a/NPLib/Detectors/MUST2/TMust2Physics.h b/NPLib/Detectors/MUST2/TMust2Physics.h index 3705b691812997690d5f73052bca730514827ad5..32338b205f9c847917f1dda04a6c9ab932dc5f05 100644 --- a/NPLib/Detectors/MUST2/TMust2Physics.h +++ b/NPLib/Detectors/MUST2/TMust2Physics.h @@ -31,6 +31,7 @@ #include "NPVDetector.h" #include "TMust2Data.h" #include "TMust2Spectra.h" + // ROOT #include "TH1.h" #include "TObject.h" @@ -43,24 +44,25 @@ using namespace std; class TMust2Spectra; class TMust2Physics : public TObject, public NPL::VDetector { - public: +public: TMust2Physics(); ~TMust2Physics(); - public: +public: void Clear(); void Clear(const Option_t*){}; - public: +public: vector<TVector2> Match_X_Y(); - bool Match_Si_CsI(int X, int Y, int CristalNbr); - bool Match_Si_SiLi(int X, int Y, int PadNbr); - bool ResolvePseudoEvent(); - int CheckEvent(); + void CheckEvent(int N); + bool Match_Si_CsI(int X, int Y, int CristalNbr); + bool Match_Si_SiLi(int X, int Y, int PadNbr); + bool ResolvePseudoEvent(); - public: +public: // Provide Physical Multiplicity - Int_t EventMultiplicity; + // Int_t EventMultiplicity; + int EventMultiplicity; // Provide a Classification of Event vector<int> EventType; @@ -95,7 +97,7 @@ class TMust2Physics : public TObject, public NPL::VDetector { // Physical Value vector<double> TotalEnergy; - public: // Innherited from VDetector Class +public: // Innherited from VDetector Class // Read stream at ConfigFile to pick-up parameters of detector // (Position,...) using Token void ReadConfiguration(NPL::InputParser parser); @@ -149,7 +151,7 @@ class TMust2Physics : public TObject, public NPL::VDetector { // Used for Online only, clear all the spectra hold by the Spectra class void ClearSpectra(); - public: // Specific to MUST2 Array +public: // Specific to MUST2 Array // Clear The PreTeated object void ClearPreTreatedData() { m_PreTreatedData->Clear(); } @@ -211,7 +213,7 @@ class TMust2Physics : public TObject, public NPL::VDetector { TVector3 GetPositionOfInteraction(const int i) const; TVector3 GetTelescopeNormal(const int i) const; - private: // Parameter used in the analysis +private: // Parameter used in the analysis // By default take EX and TY. bool m_Take_E_Y; //! bool m_Take_T_Y; //! @@ -268,39 +270,45 @@ class TMust2Physics : public TObject, public NPL::VDetector { bool m_Ignore_not_matching_SiLi; //! bool m_Ignore_not_matching_CsI; //! - private: // Root Input and Output tree classes +private: // Root Input and Output tree classes TMust2Data* m_EventData; //! TMust2Data* m_PreTreatedData; //! TMust2Physics* m_EventPhysics; //! - private: // Map of activated channel +private: // Map of activated channel map<int, vector<bool>> m_XChannelStatus; //! map<int, vector<bool>> m_YChannelStatus; //! map<int, vector<bool>> m_SiLiChannelStatus; //! map<int, vector<bool>> m_CsIChannelStatus; //! - private - : // Spatial Position of Strip Calculated on bases of detector position +private: int m_NumberOfTelescope; //! vector<vector<vector<double>>> m_StripPositionX; //! vector<vector<vector<double>>> m_StripPositionY; //! vector<vector<vector<double>>> m_StripPositionZ; //! - private: - map<int, int> m_HitStripX; //! - map<int, int> m_HitStripY; //! +public: + // Prevent to treat event with ambiguous matching beetween X and Y + int m_OrderMatch; //! + map<int, int> m_NMatchDet; //! + map<int, int> m_StripXMultDet; //! + map<int, int> m_StripYMultDet; //! + +private: + map<int, bool> m_CsIPresent; //! + map<int, bool> m_SiLiPresent; //! - private: // Spectra Class +private: // Spectra Class TMust2Spectra* m_Spectra; //! - public: +public: void WriteSpectra(); //! - public: // Spectra Getter +public: // Spectra Getter map<string, TH1*> GetSpectra(); - public: // Static constructor to be passed to the Detector Factory +public: // Static constructor to be passed to the Detector Factory static NPL::VDetector* Construct(); ClassDef(TMust2Physics, 1) // Must2Physics structure }; @@ -322,6 +330,6 @@ double fSiLi_T(const TMust2Data* Data, const int& i); // CsI double fCsI_E(const TMust2Data* Data, const int& i); double fCsI_T(const TMust2Data* Data, const int& i); -} +} // namespace MUST2_LOCAL #endif diff --git a/NPLib/Detectors/ModularLeaf/CMakeLists.txt b/NPLib/Detectors/ModularLeaf/CMakeLists.txt index 6a32ba0c2d498087542b9e429c5047ddb1fec192..104b310bb4d476a71a8a09537ded520680292c05 100644 --- a/NPLib/Detectors/ModularLeaf/CMakeLists.txt +++ b/NPLib/Detectors/ModularLeaf/CMakeLists.txt @@ -1,5 +1,5 @@ add_custom_command(OUTPUT TModularLeafPhysicsDict.cxx COMMAND ../../scripts/build_dict.sh TModularLeafPhysics.h TModularLeafPhysicsDict.cxx TModularLeafPhysics.rootmap libNPModularLeaf.dylib DEPENDS TModularLeafPhysics.h) -add_library(NPModularLeaf SHARED TModularLeafPhysics.cxx TModularLeafPhysicsDict.cxx) +add_library(NPModularLeaf SHARED TModularLeafPhysics.cxx TModularLeafPhysicsDict.cxx TModularLeafSpectra.cxx) target_link_libraries(NPModularLeaf ${ROOT_LIBRARIES} NPCore) -install(FILES TModularLeafPhysics.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) +install(FILES TModularLeafPhysics.h TModularLeafSpectra.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY}) diff --git a/NPLib/Detectors/Mugast/TMugastPhysics.cxx b/NPLib/Detectors/Mugast/TMugastPhysics.cxx index 270a284d0365d3d3a1edbb1d62de10997fbbb7f0..467811c598eecb45fbd68359d0d702818281b6db 100644 --- a/NPLib/Detectors/Mugast/TMugastPhysics.cxx +++ b/NPLib/Detectors/Mugast/TMugastPhysics.cxx @@ -524,6 +524,8 @@ map<string, TH1*> TMugastPhysics::GetSpectra() { map<string, TH1*> empty; return empty; }*/ + map<string, TH1*> empty; + return empty; } /////////////////////////////////////////////////////////////////////////// diff --git a/NPLib/Detectors/Vamos/TVamosPhysics.cxx b/NPLib/Detectors/Vamos/TVamosPhysics.cxx index fa6baa2e388ed2246b4db739e5c98b0c3e0b16bc..0ea9b62bc94a575fbdbe0992cc84d0ef1269054d 100644 --- a/NPLib/Detectors/Vamos/TVamosPhysics.cxx +++ b/NPLib/Detectors/Vamos/TVamosPhysics.cxx @@ -101,10 +101,8 @@ void TVamosPhysics::BuildPhysicalEvent() { Drift_X.push_back(m_EventData->Get_X(d)); X_Event.push_back(m_EventData->Get_X(d)); - DriftTime = m_EventData->Get_DriftTime(d); + DriftTime = m_EventData->Get_DriftTime(d)*microsecond; - /* cout<<7.5*centimeter/DriftSpeed<<endl; */ - // 7.5 centimeter is half times the height of Chamber // to separate negative Y from positive Y if (DriftTime > 7.5*centimeter/DriftSpeed) { diff --git a/NPSimulation/Detectors/Dali/Dali.cc b/NPSimulation/Detectors/Dali/Dali.cc index bbf2fa3c32561259e383a8472cd7c41818f2d491..ec87cd71ce39d7defd5d078fb027d5a8316ea706 100644 --- a/NPSimulation/Detectors/Dali/Dali.cc +++ b/NPSimulation/Detectors/Dali/Dali.cc @@ -73,7 +73,7 @@ namespace Dali_NS{ // Energy and time Resolution const double EnergyThreshold = 0*MeV; const double ResoTime = 0.0*ns; //4.5*ns ; - const double ResoEnergy = 1.36*MeV ; // mean Resolution(FWHM) 1.7% of 80MeV from slides 20170214-SAMURAI34-setup-DALI.pdf // 0.001*MeV ; + const double ResoEnergy = 0.51*MeV ; // mean Resolution(FWHM) 1.7% of 80MeV from slides 20170214-SAMURAI34-setup-DALI.pdf if 1.7% of 30MeV = 0.51 MeV // 0.001*MeV ; const double Radius = 50*mm ; const double Width = 49.76*mm ; const double Hight = 84.81*mm ; @@ -96,8 +96,8 @@ Dali::Dali(){ Logic_ArrayDali_1 =0; // RGB Color + Transparency - m_VisSquare = new G4VisAttributes(G4Colour(0, 1, 1, 0.3)); - m_VisCylinder = new G4VisAttributes(G4Colour(0, 0, 1, 0.3)); + m_VisSquare = new G4VisAttributes(G4Colour(0, 1, 1/*, 0.3*/)); + m_VisCylinder = new G4VisAttributes(G4Colour(0, 0, 1/*, 0.3*/)); } @@ -181,7 +181,7 @@ G4LogicalVolume* Dali::BuildSquareDetector(){ G4Material* Aria = MaterialManager::getInstance()->GetMaterialFromLibrary("Air"); Logic_ArrayDali_1 = new G4LogicalVolume(box_3can,Aria,"logic_ArrayDali",0,0,0); - Logic_ArrayDali_1->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1, 0.01))); + Logic_ArrayDali_1->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1, 0))); G4Box* box_can = new G4Box("Dali_BoxCan", Dali_NS::Hight*0.5, @@ -225,9 +225,9 @@ G4LogicalVolume* Dali::BuildSquareDetector(){ G4VisAttributes* Can_Attributes = new G4VisAttributes(G4Colour(0.5,0.5,0.5, .3)); m_SquareDetector_Can->SetVisAttributes(Can_Attributes); - m_Square2Detector_Can->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0.1))); + m_Square2Detector_Can->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0))); - //Extrudedbox_can->SetVisAttributes(Can_Attributes); + AriaExtrude->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0))); lAlPMT->SetVisAttributes(Can_Attributes); lMuPMT->SetVisAttributes(Can_Attributes); lTopPlatePMT->SetVisAttributes(Can_Attributes); @@ -308,7 +308,6 @@ G4LogicalVolume* Dali::BuildSquareDetector(){ 0); G4VisAttributes* MgO_Attributes = new G4VisAttributes(G4Colour(1,1,1, .3)); m_SquareDetector_CanMgO->SetVisAttributes(MgO_Attributes); - AriaExtrude->SetVisAttributes(MgO_Attributes); // NaI Volume - diff --git a/NPSimulation/Detectors/MUST2/MUST2Array.cc b/NPSimulation/Detectors/MUST2/MUST2Array.cc index 0935d0d23dd352df9d450b4b7ca979107bf69836..94860cbe3740461838caa8a64bf95949e7ea3039 100644 --- a/NPSimulation/Detectors/MUST2/MUST2Array.cc +++ b/NPSimulation/Detectors/MUST2/MUST2Array.cc @@ -22,42 +22,42 @@ * - 16 Si(Li) pad * * - 16 CsI scintillator Crystal * *****************************************************************************/ -#include"MUST2Array.hh" +#include "MUST2Array.hh" -//Geant4 -#include "G4Trd.hh" +// Geant4 #include "G4Box.hh" -#include "G4Trap.hh" -#include "G4MaterialTable.hh" +#include "G4Colour.hh" #include "G4Element.hh" #include "G4ElementTable.hh" +#include "G4MaterialTable.hh" +#include "G4PVDivision.hh" +#include "G4PVPlacement.hh" #include "G4RotationMatrix.hh" #include "G4Transform3D.hh" -#include "G4PVPlacement.hh" +#include "G4Trap.hh" +#include "G4Trd.hh" #include "G4VisAttributes.hh" -#include "G4Colour.hh" -#include "G4PVDivision.hh" #include "Randomize.hh" // NPS -#include "MaterialManager.hh" -#include "NPSDetectorFactory.hh" -#include "InteractionScorers.hh" -#include "DSSDScorers.hh" #include "CalorimeterScorers.hh" +#include "DSSDScorers.hh" +#include "InteractionScorers.hh" +#include "MaterialManager.hh" #include "NPOptionManager.h" +#include "NPSDetectorFactory.hh" // NPL #include "NPCore.h" -//ROOT +// ROOT #include "RootOutput.h" // CLHEP #include "CLHEP/Random/RandGauss.h" // STL -#include <sstream> -#include <string> #include <cmath> #include <set> +#include <sstream> +#include <string> using namespace std; using namespace CLHEP; using namespace MUST2; @@ -65,33 +65,26 @@ using namespace MUST2; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... // MUST2Array Specific Method -MUST2Array::MUST2Array(){ - m_Event = new TMust2Data() ; +MUST2Array::MUST2Array() { + m_Event = new TMust2Data(); InitializeMaterial(); - m_StripScorer=0; - m_SiLiScorer=0; - m_CsIScorer=0; -} - -MUST2Array::~MUST2Array(){ + m_StripScorer = 0; + m_SiLiScorer = 0; + m_CsIScorer = 0; } +MUST2Array::~MUST2Array() {} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void MUST2Array::AddTelescope( G4ThreeVector X1_Y1, - G4ThreeVector X128_Y1, - G4ThreeVector X1_Y128, - G4ThreeVector X128_Y128, - bool wSi, - bool wSiLi, - bool wCsI) -{ +void MUST2Array::AddTelescope(G4ThreeVector X1_Y1, G4ThreeVector X128_Y1, + G4ThreeVector X1_Y128, G4ThreeVector X128_Y128, + bool wSi, bool wSiLi, bool wCsI) { m_DefinitionType.push_back(true); - m_X1_Y1 .push_back(X1_Y1); - m_X128_Y1 .push_back(X128_Y1); - m_X1_Y128 .push_back(X1_Y128); - m_X128_Y128 .push_back(X128_Y128); + m_X1_Y1.push_back(X1_Y1); + m_X128_Y1.push_back(X128_Y1); + m_X1_Y128.push_back(X1_Y128); + m_X128_Y128.push_back(X128_Y128); m_wSi.push_back(wSi); m_wSiLi.push_back(wSiLi); m_wCsI.push_back(wCsI); @@ -104,16 +97,9 @@ void MUST2Array::AddTelescope( G4ThreeVector X1_Y1, m_beta_w.push_back(0); } -void MUST2Array::AddTelescope( G4double R, - G4double Theta, - G4double Phi, - G4double beta_u, - G4double beta_v, - G4double beta_w, - bool wSi, - bool wSiLi, - bool wCsI) -{ +void MUST2Array::AddTelescope(G4double R, G4double Theta, G4double Phi, + G4double beta_u, G4double beta_v, G4double beta_w, + bool wSi, bool wSiLi, bool wCsI) { G4ThreeVector empty = G4ThreeVector(0, 0, 0); m_DefinitionType.push_back(false); @@ -128,90 +114,95 @@ void MUST2Array::AddTelescope( G4double R, m_wSiLi.push_back(wSiLi); m_wCsI.push_back(wCsI); - m_X1_Y1 .push_back(empty) ; - m_X128_Y1 .push_back(empty) ; - m_X1_Y128 .push_back(empty) ; - m_X128_Y128 .push_back(empty) ; - + m_X1_Y1.push_back(empty); + m_X128_Y1.push_back(empty); + m_X1_Y128.push_back(empty); + m_X128_Y128.push_back(empty); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void MUST2Array::VolumeMaker( G4int TelescopeNumber, - G4ThreeVector MMpos, - G4RotationMatrix* MMrot, - bool wSi, - bool wSiLi, - bool wCsI, - G4LogicalVolume* world) -{ - G4double NbrTelescopes = TelescopeNumber; - G4String DetectorNumber; +void MUST2Array::VolumeMaker(G4int TelescopeNumber, G4ThreeVector MMpos, + G4RotationMatrix* MMrot, bool wSi, bool wSiLi, + bool wCsI, G4LogicalVolume* world) { + G4double NbrTelescopes = TelescopeNumber; + G4String DetectorNumber; std::ostringstream Number; - Number << NbrTelescopes ; + Number << NbrTelescopes; DetectorNumber = Number.str(); //////////////////////////////////////////////////////////////// ////////////// Starting Volume Definition ////////////////////// //////////////////////////////////////////////////////////////// - G4Trd* solidMM = new G4Trd("MUST2Telescope" + DetectorNumber, 0.5*FaceFront, 0.5*FaceBack, 0.5*FaceFront, 0.5*FaceBack, 0.5*Length); - G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, m_MaterialIron, "MUST2Telescope" + DetectorNumber, 0, 0, 0); - G4String Name = "MUST2Telescope" + DetectorNumber ; - - new G4PVPlacement( G4Transform3D(*MMrot, MMpos), - logicMM , - Name , - world , - false , - TelescopeNumber ); - - if (m_non_sensitive_part_visiualisation){ - G4VisAttributes* FrameVisAtt = new G4VisAttributes(G4Colour(0.80, 0.80, 0.80)); + G4Trd* solidMM = new G4Trd("MUST2Telescope" + DetectorNumber, 0.5 * FaceFront, + 0.5 * FaceBack, 0.5 * FaceFront, 0.5 * FaceBack, + 0.5 * Length); + G4LogicalVolume* logicMM = new G4LogicalVolume( + solidMM, m_MaterialIron, "MUST2Telescope" + DetectorNumber, 0, 0, 0); + G4String Name = "MUST2Telescope" + DetectorNumber; + + new G4PVPlacement(G4Transform3D(*MMrot, MMpos), logicMM, Name, world, false, + TelescopeNumber); + + if (m_non_sensitive_part_visiualisation) { + G4VisAttributes* FrameVisAtt + = new G4VisAttributes(G4Colour(0.80, 0.80, 0.80)); FrameVisAtt->SetForceWireframe(true); - logicMM->SetVisAttributes(FrameVisAtt) ; - } - else logicMM->SetVisAttributes(G4VisAttributes::Invisible) ; + logicMM->SetVisAttributes(FrameVisAtt); + } else + logicMM->SetVisAttributes(G4VisAttributes::Invisible); G4ThreeVector positionVacBox = G4ThreeVector(0, 0, VacBox_PosZ); - G4Trd* solidVacBox = new G4Trd("solidVacBox", 0.5*SiliconFace, 0.5*CsIFaceFront, 0.5*SiliconFace, 0.5*CsIFaceFront, 0.5*VacBoxThickness); - G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0); + G4Trd* solidVacBox + = new G4Trd("solidVacBox", 0.5 * SiliconFace, 0.5 * CsIFaceFront, + 0.5 * SiliconFace, 0.5 * CsIFaceFront, 0.5 * VacBoxThickness); + G4LogicalVolume* logicVacBox = new G4LogicalVolume( + solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0); - - new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, TelescopeNumber ); + new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, + false, TelescopeNumber); logicVacBox->SetVisAttributes(G4VisAttributes::Invisible); - G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)) ; + G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)); //////////////////////////////////////////////////////////////// /////////////////Si Strip Construction////////////////////////// //////////////////////////////////////////////////////////////// if (wSi) { - G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, AluStripFront_PosZ); - G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripBack_PosZ); + G4ThreeVector positionAluStripFront + = G4ThreeVector(0, 0, AluStripFront_PosZ); + G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripBack_PosZ); - G4Box* solidAluStrip = new G4Box("AluBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*AluStripThickness); - G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0); + G4Box* solidAluStrip + = new G4Box("AluBox", 0.5 * SiliconFace, 0.5 * SiliconFace, + 0.5 * AluStripThickness); + G4LogicalVolume* logicAluStrip = new G4LogicalVolume( + solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0); - new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicMM, false, TelescopeNumber ); - new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicMM, false, TelescopeNumber ); + new G4PVPlacement(0, positionAluStripFront, logicAluStrip, + Name + "_AluStripFront", logicMM, false, TelescopeNumber); + new G4PVPlacement(0, positionAluStripBack, logicAluStrip, + Name + "_AluStripBack", logicMM, false, TelescopeNumber); logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible); - G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ); - - G4Box* solidSilicon = new G4Box("solidSilicon", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconThickness); - G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0); + G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ); - new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicMM, false, TelescopeNumber ); + G4Box* solidSilicon = new G4Box("solidSilicon", 0.5 * SiliconFace, + 0.5 * SiliconFace, 0.5 * SiliconThickness); + G4LogicalVolume* logicSilicon = new G4LogicalVolume( + solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0); + new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", + logicMM, false, TelescopeNumber); - ///Set Silicon strip sensible + /// Set Silicon strip sensible logicSilicon->SetSensitiveDetector(m_StripScorer); - ///Visualisation of Silicon Strip - logicSilicon->SetVisAttributes(SiliconVisAtt) ; + /// Visualisation of Silicon Strip + logicSilicon->SetVisAttributes(SiliconVisAtt); } //////////////////////////////////////////////////////////////// @@ -219,185 +210,207 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, //////////////////////////////////////////////////////////////// if (wSiLi) { - G4double SiLiSpace = 8 * mm; - G4RotationMatrix* rotSiLi = new G4RotationMatrix(0,0,0); - G4Box* solidSiLi = new G4Box("SiLi", 0.5*SiliconFace+0.5*SiLiSpace, 0.5*SiliconFace, 0.5*SiLiThickness); - G4LogicalVolume* logicSiLi = new G4LogicalVolume(solidSiLi, m_MaterialAluminium, Name + "_SiLi" , 0, 0, 0); + G4double SiLiSpace = 8 * mm; + G4RotationMatrix* rotSiLi = new G4RotationMatrix(0, 0, 0); + G4Box* solidSiLi = new G4Box("SiLi", 0.5 * SiliconFace + 0.5 * SiLiSpace, + 0.5 * SiliconFace, 0.5 * SiLiThickness); + G4LogicalVolume* logicSiLi = new G4LogicalVolume( + solidSiLi, m_MaterialAluminium, Name + "_SiLi", 0, 0, 0); logicSiLi->SetVisAttributes(G4VisAttributes::Invisible); - new G4PVPlacement( G4Transform3D(*rotSiLi, G4ThreeVector(0,0,0) ) , - logicSiLi , - Name + "_SiLi", - logicVacBox , - false , - 0); + new G4PVPlacement(G4Transform3D(*rotSiLi, G4ThreeVector(0, 0, 0)), + logicSiLi, Name + "_SiLi", logicVacBox, false, 0); // SiLi are placed inside of the VacBox... // Left/Right define when looking to detector from Si to CsI - G4double SiLi_HighY_Upper = 19.86 * mm; + G4double SiLi_HighY_Upper = 19.86 * mm; G4double SiLi_HighY_Center = 25.39 * mm; G4double SiLi_WidthX_Left = 22.85 * mm; - G4double SiLi_WidthX_Right = 24.9 * mm; - G4double SiLi_ShiftX = 0.775 * mm; + G4double SiLi_WidthX_Right = 24.9 * mm; + G4double SiLi_ShiftX = 0.775 * mm; // SiLi are organized by two group of 8 Up(9 to 15) and Down(1 to 8). - G4ThreeVector ShiftSiLiUp = G4ThreeVector(-0.25 * SiliconFace - 0.5 * SiLiSpace, 0, 0) ; - G4ThreeVector ShiftSiLiDown = G4ThreeVector(0.25 * SiliconFace + 0.5 * SiLiSpace, 0, 0) ; + G4ThreeVector ShiftSiLiUp + = G4ThreeVector(-0.25 * SiliconFace - 0.5 * SiLiSpace, 0, 0); + G4ThreeVector ShiftSiLiDown + = G4ThreeVector(0.25 * SiliconFace + 0.5 * SiLiSpace, 0, 0); // SiLi : left side of SiLi detector - G4Box* solidSiLi_LT = new G4Box("SiLi_LT", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); - G4Box* solidSiLi_RT = new G4Box("SiLi_RT", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); - G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); - G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); - G4Box* solidSiLi_LB = new G4Box("SiLi_LB", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); - G4Box* solidSiLi_RB = new G4Box("SiLi_RB", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Upper, 0.5*SiLiThickness); - G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2", 0.5*SiLi_WidthX_Left, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); - G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2", 0.5*SiLi_WidthX_Right, 0.5*SiLi_HighY_Center, 0.5*SiLiThickness); - - G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume(solidSiLi_LT , m_MaterialSilicon, "SiLi_LT" , 0, 0, 0); - G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume(solidSiLi_RT , m_MaterialSilicon, "SiLi_RT" , 0, 0, 0); - G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume(solidSiLi_LC1, m_MaterialSilicon, "SiLi_LC1", 0, 0, 0); - G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume(solidSiLi_RC1, m_MaterialSilicon, "SiLi_RC1", 0, 0, 0); - G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume(solidSiLi_LB , m_MaterialSilicon, "SiLi_LB" , 0, 0, 0); - G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume(solidSiLi_RB , m_MaterialSilicon, "SiLi_RB" , 0, 0, 0); - G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume(solidSiLi_LC2, m_MaterialSilicon, "SiLi_LC2", 0, 0, 0); - G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume(solidSiLi_RC2, m_MaterialSilicon, "SiLi_RC2", 0, 0, 0); + G4Box* solidSiLi_LT + = new G4Box("SiLi_LT", 0.5 * SiLi_WidthX_Left, 0.5 * SiLi_HighY_Upper, + 0.5 * SiLiThickness); + G4Box* solidSiLi_RT + = new G4Box("SiLi_RT", 0.5 * SiLi_WidthX_Right, 0.5 * SiLi_HighY_Upper, + 0.5 * SiLiThickness); + G4Box* solidSiLi_LC1 + = new G4Box("SiLi_LC1", 0.5 * SiLi_WidthX_Left, 0.5 * SiLi_HighY_Center, + 0.5 * SiLiThickness); + G4Box* solidSiLi_RC1 + = new G4Box("SiLi_RC1", 0.5 * SiLi_WidthX_Right, + 0.5 * SiLi_HighY_Center, 0.5 * SiLiThickness); + G4Box* solidSiLi_LB + = new G4Box("SiLi_LB", 0.5 * SiLi_WidthX_Left, 0.5 * SiLi_HighY_Upper, + 0.5 * SiLiThickness); + G4Box* solidSiLi_RB + = new G4Box("SiLi_RB", 0.5 * SiLi_WidthX_Right, 0.5 * SiLi_HighY_Upper, + 0.5 * SiLiThickness); + G4Box* solidSiLi_LC2 + = new G4Box("SiLi_LC2", 0.5 * SiLi_WidthX_Left, 0.5 * SiLi_HighY_Center, + 0.5 * SiLiThickness); + G4Box* solidSiLi_RC2 + = new G4Box("SiLi_RC2", 0.5 * SiLi_WidthX_Right, + 0.5 * SiLi_HighY_Center, 0.5 * SiLiThickness); + + G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume( + solidSiLi_LT, m_MaterialSilicon, "SiLi_LT", 0, 0, 0); + G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume( + solidSiLi_RT, m_MaterialSilicon, "SiLi_RT", 0, 0, 0); + G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume( + solidSiLi_LC1, m_MaterialSilicon, "SiLi_LC1", 0, 0, 0); + G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume( + solidSiLi_RC1, m_MaterialSilicon, "SiLi_RC1", 0, 0, 0); + G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume( + solidSiLi_LB, m_MaterialSilicon, "SiLi_LB", 0, 0, 0); + G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume( + solidSiLi_RB, m_MaterialSilicon, "SiLi_RB", 0, 0, 0); + G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume( + solidSiLi_LC2, m_MaterialSilicon, "SiLi_LC2", 0, 0, 0); + G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume( + solidSiLi_RC2, m_MaterialSilicon, "SiLi_RC2", 0, 0, 0); G4double interSiLi = 0.5 * mm; // Top - G4ThreeVector positionSiLi_LT_up = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LT_up = G4ThreeVector( + -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, 0); - positionSiLi_LT_up += ShiftSiLiUp ; + positionSiLi_LT_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_RT_up = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RT_up = G4ThreeVector( + 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, 0); - positionSiLi_RT_up += ShiftSiLiUp ; + positionSiLi_RT_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_LC1_up = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LC1_up + = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, 0); - positionSiLi_LC1_up += ShiftSiLiUp ; + positionSiLi_LC1_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_RC1_up = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_RC1_up + = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, 0); - positionSiLi_RC1_up += ShiftSiLiUp ; + positionSiLi_RC1_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_LB_up = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , - 0); + G4ThreeVector positionSiLi_LB_up = G4ThreeVector( + -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, 0); - positionSiLi_LB_up += ShiftSiLiUp ; + positionSiLi_LB_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_RB_up = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX , - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , - 0); + G4ThreeVector positionSiLi_RB_up = G4ThreeVector( + 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, 0); - positionSiLi_RB_up += ShiftSiLiUp ; + positionSiLi_RB_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_LC2_up = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, - 0); + G4ThreeVector positionSiLi_LC2_up + = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, 0); - positionSiLi_LC2_up += ShiftSiLiUp ; + positionSiLi_LC2_up += ShiftSiLiUp; - G4ThreeVector positionSiLi_RC2_up = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX , - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi , - 0); - - positionSiLi_RC2_up += ShiftSiLiUp ; + G4ThreeVector positionSiLi_RC2_up + = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, 0); + positionSiLi_RC2_up += ShiftSiLiUp; // Down - G4ThreeVector positionSiLi_LT_down = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); - - positionSiLi_LT_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_RT_down = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, - 0); - - positionSiLi_RT_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_LC1_down = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, - 0); - - positionSiLi_LC1_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_RC1_down = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - 0.5 * SiLi_HighY_Center + 0.5 * interSiLi , - 0); - - positionSiLi_RC1_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_LB_down = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, - 0); - - positionSiLi_LB_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_RB_down = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, - 0); - - positionSiLi_RB_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_LC2_down = - G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, - 0); - - positionSiLi_LC2_down += ShiftSiLiDown ; - - G4ThreeVector positionSiLi_RC2_down = - G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, - -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, - 0); - - positionSiLi_RC2_down += ShiftSiLiDown ; - - new G4PVPlacement(0 , positionSiLi_RT_down , logicSiLi_RT , Name + "_SiLi_Pad1" , logicSiLi , false ,1); - new G4PVPlacement(0 , positionSiLi_LT_down , logicSiLi_LT , Name + "_SiLi_Pad2" , logicSiLi , false ,2); - new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" , logicSiLi , false ,3); - new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" , logicSiLi , false ,4); - new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" , logicSiLi , false ,5); - new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" , logicSiLi , false ,6); - new G4PVPlacement(0 , positionSiLi_LB_down , logicSiLi_LB , Name + "_SiLi_Pad7" , logicSiLi , false ,7); - new G4PVPlacement(0 , positionSiLi_RB_down , logicSiLi_RB , Name + "_SiLi_Pad8" , logicSiLi , false ,8); - new G4PVPlacement(0 , positionSiLi_LT_up , logicSiLi_LT , Name + "_SiLi_Pad9" , logicSiLi , false ,9); - new G4PVPlacement(0 , positionSiLi_RT_up , logicSiLi_RT , Name + "_SiLi_Pad10" , logicSiLi , false ,10); - new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" , logicSiLi , false ,11); - new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" , logicSiLi , false ,12); - new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" , logicSiLi , false ,13); - new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" , logicSiLi , false ,14); - new G4PVPlacement(0 , positionSiLi_RB_up , logicSiLi_RB , Name + "_SiLi_Pad15" , logicSiLi , false ,15); - new G4PVPlacement(0 , positionSiLi_LB_up , logicSiLi_LB , Name + "_SiLi_Pad16" , logicSiLi , false ,16); - + G4ThreeVector positionSiLi_LT_down = G4ThreeVector( + -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, 0); + + positionSiLi_LT_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_RT_down = G4ThreeVector( + 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi, 0); + + positionSiLi_RT_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_LC1_down + = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, 0); + + positionSiLi_LC1_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_RC1_down + = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi, 0); + + positionSiLi_RC1_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_LB_down = G4ThreeVector( + -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, 0); + + positionSiLi_LB_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_RB_down = G4ThreeVector( + 0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi, 0); + + positionSiLi_RB_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_LC2_down + = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX, + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, 0); + + positionSiLi_LC2_down += ShiftSiLiDown; + + G4ThreeVector positionSiLi_RC2_down + = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX, + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi, 0); + + positionSiLi_RC2_down += ShiftSiLiDown; + + new G4PVPlacement(0, positionSiLi_RT_down, logicSiLi_RT, + Name + "_SiLi_Pad1", logicSiLi, false, 1); + new G4PVPlacement(0, positionSiLi_LT_down, logicSiLi_LT, + Name + "_SiLi_Pad2", logicSiLi, false, 2); + new G4PVPlacement(0, positionSiLi_RC1_down, logicSiLi_RC1, + Name + "_SiLi_Pad3", logicSiLi, false, 3); + new G4PVPlacement(0, positionSiLi_LC1_down, logicSiLi_LC1, + Name + "_SiLi_Pad4", logicSiLi, false, 4); + new G4PVPlacement(0, positionSiLi_LC2_down, logicSiLi_LC2, + Name + "_SiLi_Pad5", logicSiLi, false, 5); + new G4PVPlacement(0, positionSiLi_RC2_down, logicSiLi_RC2, + Name + "_SiLi_Pad6", logicSiLi, false, 6); + new G4PVPlacement(0, positionSiLi_LB_down, logicSiLi_LB, + Name + "_SiLi_Pad7", logicSiLi, false, 7); + new G4PVPlacement(0, positionSiLi_RB_down, logicSiLi_RB, + Name + "_SiLi_Pad8", logicSiLi, false, 8); + new G4PVPlacement(0, positionSiLi_LT_up, logicSiLi_LT, Name + "_SiLi_Pad9", + logicSiLi, false, 9); + new G4PVPlacement(0, positionSiLi_RT_up, logicSiLi_RT, Name + "_SiLi_Pad10", + logicSiLi, false, 10); + new G4PVPlacement(0, positionSiLi_LC1_up, logicSiLi_LC1, + Name + "_SiLi_Pad11", logicSiLi, false, 11); + new G4PVPlacement(0, positionSiLi_RC1_up, logicSiLi_RC1, + Name + "_SiLi_Pad12", logicSiLi, false, 12); + new G4PVPlacement(0, positionSiLi_RC2_up, logicSiLi_RC2, + Name + "_SiLi_Pad13", logicSiLi, false, 13); + new G4PVPlacement(0, positionSiLi_LC2_up, logicSiLi_LC2, + Name + "_SiLi_Pad14", logicSiLi, false, 14); + new G4PVPlacement(0, positionSiLi_RB_up, logicSiLi_RB, Name + "_SiLi_Pad15", + logicSiLi, false, 15); + new G4PVPlacement(0, positionSiLi_LB_up, logicSiLi_LB, Name + "_SiLi_Pad16", + logicSiLi, false, 16); // Set SiLi sensible logicSiLi_LT->SetSensitiveDetector(m_SiLiScorer); @@ -423,7 +436,6 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, logicSiLi_LC2->SetVisAttributes(SiLiVisAtt); logicSiLi_RC2->SetVisAttributes(SiLiVisAtt); - delete rotSiLi; } @@ -432,63 +444,104 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, //////////////////////////////////////////////////////////////// if (wCsI) { - m_MaterialMyl = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar"); - m_MaterialCsI = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI"); + m_MaterialMyl + = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar"); + m_MaterialCsI + = MaterialManager::getInstance()->GetMaterialFromLibrary("CsI"); G4ThreeVector positionCsI = G4ThreeVector(0, 0, CsI_PosZ); - G4Trd* solidCsI = new G4Trd("csI", 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIThickness); - - G4LogicalVolume* logicCsI = new G4LogicalVolume(solidCsI, m_MaterialAluminium, Name + "_CsI_Mylar", 0, 0, 0); - new G4PVPlacement(0, positionCsI, logicCsI, Name + "_CsI_Mylar", logicMM, false, 0); + G4Trd* solidCsI + = new G4Trd("csI", 0.5 * CsIFaceFront, 0.5 * CsIFaceBack, + 0.5 * CsIFaceFront, 0.5 * CsIFaceBack, 0.5 * CsIThickness); - G4ThreeVector positionMylarCsI = G4ThreeVector(0, 0, MylarCsIThickness * 0.5 - CsIThickness * 0.5); + G4LogicalVolume* logicCsI = new G4LogicalVolume( + solidCsI, m_MaterialAluminium, Name + "_CsI_Mylar", 0, 0, 0); + new G4PVPlacement(0, positionCsI, logicCsI, Name + "_CsI_Mylar", logicMM, + false, 0); - G4Box* solidMylarCsI = new G4Box("MylarCsIBox", 0.5*CsIFaceFront, 0.5*CsIFaceFront, 0.5*MylarCsIThickness); - G4LogicalVolume* logicMylarCsI = new G4LogicalVolume(solidMylarCsI, m_MaterialMyl, Name + "_CsI_Mylar", 0, 0, 0); + G4ThreeVector positionMylarCsI + = G4ThreeVector(0, 0, MylarCsIThickness * 0.5 - CsIThickness * 0.5); - new G4PVPlacement(0, positionMylarCsI, logicMylarCsI, Name + "_CsI_Mylar", logicCsI, false, 0); + G4Box* solidMylarCsI + = new G4Box("MylarCsIBox", 0.5 * CsIFaceFront, 0.5 * CsIFaceFront, + 0.5 * MylarCsIThickness); + G4LogicalVolume* logicMylarCsI = new G4LogicalVolume( + solidMylarCsI, m_MaterialMyl, Name + "_CsI_Mylar", 0, 0, 0); + new G4PVPlacement(0, positionMylarCsI, logicMylarCsI, Name + "_CsI_Mylar", + logicCsI, false, 0); logicCsI->SetVisAttributes(G4VisAttributes::Invisible); logicMylarCsI->SetVisAttributes(G4VisAttributes::Invisible); // Cristal1 - G4Trap* solidCristal1 = new G4Trap("Cristal1", 40.*mm / 2., 6.693896*deg, 41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal1 = new G4LogicalVolume(solidCristal1, m_MaterialCsI, Name + "_CsI_Cristal1", 0, 0, 0); + G4Trap* solidCristal1 = new G4Trap( + "Cristal1", 40. * mm / 2., 6.693896 * deg, 41.97814 * deg, + 33.1 * mm / 2., 37.39 * mm / 2., 37.39 * mm / 2., 0. * deg, + 26.9 * mm / 2., 30.41 * mm / 2., 30.41 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal1 = new G4LogicalVolume( + solidCristal1, m_MaterialCsI, Name + "_CsI_Cristal1", 0, 0, 0); // Cristal2 - G4Trap* solidCristal2 = new G4Trap("Cristal2", 40.*mm / 2., 17.8836*deg, (74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal2 = new G4LogicalVolume(solidCristal2, m_MaterialCsI, Name + "_CsI_Cristal2", 0, 0, 0); + G4Trap* solidCristal2 = new G4Trap( + "Cristal2", 40. * mm / 2., 17.8836 * deg, (74.3122 + 180) * deg, + 43.49 * mm / 2., 37.39 * mm / 2., 37.39 * mm / 2., 0. * deg, + 31.0377 * mm / 2., 30.41 * mm / 2., 30.41 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal2 = new G4LogicalVolume( + solidCristal2, m_MaterialCsI, Name + "_CsI_Cristal2", 0, 0, 0); // Cristal3 - G4Trap* solidCristal3 = new G4Trap("Cristal3", 40.*mm / 2., 18.243*deg, 13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal3 = new G4LogicalVolume(solidCristal3, m_MaterialCsI, Name + "_CsI_Cristal3", 0, 0, 0); + G4Trap* solidCristal3 = new G4Trap( + "Cristal3", 40. * mm / 2., 18.243 * deg, 13.5988 * deg, 33.11 * mm / 2., + 39.25 * mm / 2., 39.25 * mm / 2., 0. * deg, 26.91 * mm / 2., + 27.58 * mm / 2., 27.58 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal3 = new G4LogicalVolume( + solidCristal3, m_MaterialCsI, Name + "_CsI_Cristal3", 0, 0, 0); // Cristal4 - G4Trap* solidCristal4 = new G4Trap("Cristal4", 40.*mm / 2., 24.0482*deg, 44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal4 = new G4LogicalVolume(solidCristal4, m_MaterialCsI, Name + "_CsI_Cristal4", 0, 0, 0); - + G4Trap* solidCristal4 = new G4Trap( + "Cristal4", 40. * mm / 2., 24.0482 * deg, 44.1148 * deg, + 43.49 * mm / 2., 39.19 * mm / 2., 39.19 * mm / 2., 0. * deg, + 31.04 * mm / 2., 27.52 * mm / 2., 27.52 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal4 = new G4LogicalVolume( + solidCristal4, m_MaterialCsI, Name + "_CsI_Cristal4", 0, 0, 0); // Cristal1s - G4Trap* solidCristal1s = new G4Trap("Cristal1s", 40.*mm / 2., 6.693896*deg, -41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal1s = new G4LogicalVolume(solidCristal1s, m_MaterialCsI, Name + "_CsI_Cristal1s", 0, 0, 0); + G4Trap* solidCristal1s = new G4Trap( + "Cristal1s", 40. * mm / 2., 6.693896 * deg, -41.97814 * deg, + 33.1 * mm / 2., 37.39 * mm / 2., 37.39 * mm / 2., 0. * deg, + 26.9 * mm / 2., 30.41 * mm / 2., 30.41 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal1s = new G4LogicalVolume( + solidCristal1s, m_MaterialCsI, Name + "_CsI_Cristal1s", 0, 0, 0); // Cristal2s - G4Trap* solidCristal2s = new G4Trap("Cristal2s", 40.*mm / 2., 17.8836*deg, -(74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal2s = new G4LogicalVolume(solidCristal2s, m_MaterialCsI, Name + "_CsI_Cristal2s", 0, 0, 0); + G4Trap* solidCristal2s = new G4Trap( + "Cristal2s", 40. * mm / 2., 17.8836 * deg, -(74.3122 + 180) * deg, + 43.49 * mm / 2., 37.39 * mm / 2., 37.39 * mm / 2., 0. * deg, + 31.0377 * mm / 2., 30.41 * mm / 2., 30.41 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal2s = new G4LogicalVolume( + solidCristal2s, m_MaterialCsI, Name + "_CsI_Cristal2s", 0, 0, 0); // Cristal3s - G4Trap* solidCristal3s = new G4Trap("Cristal3s", 40.*mm / 2., 18.243*deg, -13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal3s = new G4LogicalVolume(solidCristal3s, m_MaterialCsI, Name + "_CsI_Cristal3s", 0, 0, 0); + G4Trap* solidCristal3s = new G4Trap( + "Cristal3s", 40. * mm / 2., 18.243 * deg, -13.5988 * deg, + 33.11 * mm / 2., 39.25 * mm / 2., 39.25 * mm / 2., 0. * deg, + 26.91 * mm / 2., 27.58 * mm / 2., 27.58 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal3s = new G4LogicalVolume( + solidCristal3s, m_MaterialCsI, Name + "_CsI_Cristal3s", 0, 0, 0); // Cristal4s - G4Trap* solidCristal4s = new G4Trap("Cristal4s", 40.*mm / 2., 24.0482*deg, -44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg); - G4LogicalVolume* logicCristal4s = new G4LogicalVolume(solidCristal4s, m_MaterialCsI, Name + "_CsI_Cristal4s", 0, 0, 0); + G4Trap* solidCristal4s = new G4Trap( + "Cristal4s", 40. * mm / 2., 24.0482 * deg, -44.1148 * deg, + 43.49 * mm / 2., 39.19 * mm / 2., 39.19 * mm / 2., 0. * deg, + 31.04 * mm / 2., 27.52 * mm / 2., 27.52 * mm / 2., 0. * deg); + G4LogicalVolume* logicCristal4s = new G4LogicalVolume( + solidCristal4s, m_MaterialCsI, Name + "_CsI_Cristal4s", 0, 0, 0); G4double XEdge1 = 16.96 * mm + DistInterCsI * 0.5; G4double YEdge1 = 15.01 * mm + DistInterCsI * 0.5; @@ -500,43 +553,65 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, G4ThreeVector positionCristal3 = G4ThreeVector(-XEdge2, YEdge1, 0); G4ThreeVector positionCristal4 = G4ThreeVector(-XEdge2, YEdge2, 0); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1, logicCristal1, Name + "_CsI_Cristal1", logicCsI, false, 1); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2, logicCristal2, Name + "_CsI_Cristal2", logicCsI, false, 2); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3, logicCristal3, Name + "_CsI_Cristal3", logicCsI, false, 3); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4, logicCristal4, Name + "_CsI_Cristal4", logicCsI, false, 4); - + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1, logicCristal1, + Name + "_CsI_Cristal1", logicCsI, false, 6); // 1 + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2, logicCristal2, + Name + "_CsI_Cristal2", logicCsI, false, 3); // 2 + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3, logicCristal3, + Name + "_CsI_Cristal3", logicCsI, false, 5); // 3 + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4, logicCristal4, + Name + "_CsI_Cristal4", logicCsI, false, 4); // 4 G4ThreeVector positionCristal1b = G4ThreeVector(XEdge1, -YEdge1, 0 * mm); G4ThreeVector positionCristal2b = G4ThreeVector(XEdge1, -YEdge2, 0); G4ThreeVector positionCristal3b = G4ThreeVector(XEdge2, -YEdge1, 0); G4ThreeVector positionCristal4b = G4ThreeVector(XEdge2, -YEdge2, 0); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1b, logicCristal1, Name + "_CsI_Cristal5", logicCsI, false, 5); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal2b, logicCristal2, Name + "_CsI_Cristal6", logicCsI, false, 6); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3b, logicCristal3, Name + "_CsI_Cristal7", logicCsI, false, 7); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4b, logicCristal4, Name + "_CsI_Cristal8", logicCsI, false, 8); + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1b, + logicCristal1, Name + "_CsI_Cristal5", logicCsI, false, + 11); // 5 + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal2b, logicCristal2, + Name + "_CsI_Cristal6", logicCsI, false, 15); // 6 + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3b, + logicCristal3, Name + "_CsI_Cristal7", logicCsI, false, + 12); // 7 + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4b, + logicCristal4, Name + "_CsI_Cristal8", logicCsI, false, + 16); // 8 G4ThreeVector positionCristal1s = G4ThreeVector(-XEdge1, -YEdge1, 0 * mm); G4ThreeVector positionCristal2s = G4ThreeVector(-XEdge1, -YEdge2, 0); G4ThreeVector positionCristal3s = G4ThreeVector(-XEdge2, -YEdge1, 0); G4ThreeVector positionCristal4s = G4ThreeVector(-XEdge2, -YEdge2, 0); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1s, logicCristal1s, Name + "_CsI_Cristal9", logicCsI, false, 9); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2s, logicCristal2s, Name + "_CsI_Cristal10", logicCsI, false, 10); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3s, logicCristal3s, Name + "_CsI_Cristal11", logicCsI, false, 11); - new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4s, logicCristal4s, Name + "_CsI_Cristal12", logicCsI, false, 12); + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1s, logicCristal1s, + Name + "_CsI_Cristal9", logicCsI, false, 10); // 9 + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2s, + logicCristal2s, Name + "_CsI_Cristal10", logicCsI, false, + 14); // 10 + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3s, logicCristal3s, + Name + "_CsI_Cristal11", logicCsI, false, 9); // 11 + new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4s, logicCristal4s, + Name + "_CsI_Cristal12", logicCsI, false, 13); // 12 G4ThreeVector positionCristal1sb = G4ThreeVector(XEdge1, YEdge1, 0 * mm); G4ThreeVector positionCristal2sb = G4ThreeVector(XEdge1, YEdge2, 0); G4ThreeVector positionCristal3sb = G4ThreeVector(XEdge2, YEdge1, 0); G4ThreeVector positionCristal4sb = G4ThreeVector(XEdge2, YEdge2, 0); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1sb, logicCristal1s, Name + "_CsI_Cristal13", logicCsI, false, 13); - new G4PVPlacement(Rotation(180, 0, 0), positionCristal2sb, logicCristal2s, Name + "_CsI_Cristal14", logicCsI, false, 14); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3sb, logicCristal3s, Name + "_CsI_Cristal15", logicCsI, false, 15); - new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4sb, logicCristal4s, Name + "_CsI_Cristal16", logicCsI, false, 16); - - ///Set CsI sensible + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1sb, + logicCristal1s, Name + "_CsI_Cristal13", logicCsI, false, + 7); // 13 + new G4PVPlacement(Rotation(180, 0, 0), positionCristal2sb, logicCristal2s, + Name + "_CsI_Cristal14", logicCsI, false, 2); // 14 + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3sb, + logicCristal3s, Name + "_CsI_Cristal15", logicCsI, false, + 8); // 15 + new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4sb, + logicCristal4s, Name + "_CsI_Cristal16", logicCsI, false, + 1); // 16 + + /// Set CsI sensible logicCristal1->SetSensitiveDetector(m_CsIScorer); logicCristal2->SetSensitiveDetector(m_CsIScorer); logicCristal3->SetSensitiveDetector(m_CsIScorer); @@ -547,17 +622,16 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, logicCristal3s->SetSensitiveDetector(m_CsIScorer); logicCristal4s->SetSensitiveDetector(m_CsIScorer); - //Mark blue a CsI corner crystal to see telescope orientation + // Mark blue a CsI corner crystal to see telescope orientation G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1, 0.5, 0)); - logicCristal1 ->SetVisAttributes(CsIVisAtt); - logicCristal2 ->SetVisAttributes(CsIVisAtt); - logicCristal3 ->SetVisAttributes(CsIVisAtt); - logicCristal4 ->SetVisAttributes(CsIVisAtt); - logicCristal1s ->SetVisAttributes(CsIVisAtt); - logicCristal2s ->SetVisAttributes(CsIVisAtt); - logicCristal3s ->SetVisAttributes(CsIVisAtt); - logicCristal4s ->SetVisAttributes(CsIVisAtt); - + logicCristal1->SetVisAttributes(CsIVisAtt); + logicCristal2->SetVisAttributes(CsIVisAtt); + logicCristal3->SetVisAttributes(CsIVisAtt); + logicCristal4->SetVisAttributes(CsIVisAtt); + logicCristal1s->SetVisAttributes(CsIVisAtt); + logicCristal2s->SetVisAttributes(CsIVisAtt); + logicCristal3s->SetVisAttributes(CsIVisAtt); + logicCristal4s->SetVisAttributes(CsIVisAtt); } } @@ -568,114 +642,121 @@ void MUST2Array::VolumeMaker( G4int TelescopeNumber, // Read stream at Configfile to pick-up parameters of detector (Position,...) // Called in DetecorConstruction::ReadDetextorConfiguration Method -void MUST2Array::ReadConfiguration(NPL::InputParser parser){ +void MUST2Array::ReadConfiguration(NPL::InputParser parser) { vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("M2Telescope"); - if(NPOptionManager::getInstance()->GetVerboseLevel()) + if (NPOptionManager::getInstance()->GetVerboseLevel()) cout << "//// " << blocks.size() << " telescope found" << endl; - for(unsigned int i = 0 ; i < blocks.size() ; i++){ - if(NPOptionManager::getInstance()->GetVerboseLevel()) - cout << endl << "//// Must 2 Telecope " << i+1 << endl; + for (unsigned int i = 0; i < blocks.size(); i++) { + if (NPOptionManager::getInstance()->GetVerboseLevel()) + cout << endl << "//// Must 2 Telecope " << i + 1 << endl; // Cartesian Case - vector<string> cart = {"X1_Y1","X1_Y128","X128_Y1","X128_Y128","SI","SILI","CSI"}; + vector<string> cart + = {"X1_Y1", "X1_Y128", "X128_Y1", "X128_Y128", "SI", "SILI", "CSI"}; // Spherical Case - vector<string> sphe= {"R","THETA","PHI","BETA","SI","SILI","CSI"}; - - if(blocks[i]->HasTokenList(cart)){ - G4ThreeVector A = NPS::ConvertVector( blocks[i]->GetTVector3("X1_Y1","mm")); - G4ThreeVector B = NPS::ConvertVector( blocks[i]->GetTVector3("X128_Y1","mm")); - G4ThreeVector C = NPS::ConvertVector( blocks[i]->GetTVector3("X1_Y128","mm")); - G4ThreeVector D = NPS::ConvertVector( blocks[i]->GetTVector3("X128_Y128","mm")); - int SI = blocks[i]->GetInt("SI"); + vector<string> sphe = {"R", "THETA", "PHI", "BETA", "SI", "SILI", "CSI"}; + + if (blocks[i]->HasTokenList(cart)) { + G4ThreeVector A + = NPS::ConvertVector(blocks[i]->GetTVector3("X1_Y1", "mm")); + G4ThreeVector B + = NPS::ConvertVector(blocks[i]->GetTVector3("X128_Y1", "mm")); + G4ThreeVector C + = NPS::ConvertVector(blocks[i]->GetTVector3("X1_Y128", "mm")); + G4ThreeVector D + = NPS::ConvertVector(blocks[i]->GetTVector3("X128_Y128", "mm")); + int SI = blocks[i]->GetInt("SI"); int SILI = blocks[i]->GetInt("SILI"); - int CSI = blocks[i]->GetInt("CSI"); - AddTelescope(A,B,C,D,SI==1,SILI==1,CSI==1) ; + int CSI = blocks[i]->GetInt("CSI"); + AddTelescope(A, B, C, D, SI == 1, SILI == 1, CSI == 1); } - else if(blocks[i]->HasTokenList(sphe)){ - - double Theta = blocks[i]->GetDouble("THETA","deg"); - double Phi= blocks[i]->GetDouble("PHI","deg"); - double R = blocks[i]->GetDouble("R","mm"); - vector<double> beta = blocks[i]->GetVectorDouble("BETA","deg"); - int SI = blocks[i]->GetInt("SI"); - int SILI = blocks[i]->GetInt("SILI"); - int CSI = blocks[i]->GetInt("CSI"); - AddTelescope( R,Theta,Phi,beta[0],beta[1],beta[2],SI==1,SILI==1,CSI==1); + else if (blocks[i]->HasTokenList(sphe)) { + + double Theta = blocks[i]->GetDouble("THETA", "deg"); + double Phi = blocks[i]->GetDouble("PHI", "deg"); + double R = blocks[i]->GetDouble("R", "mm"); + vector<double> beta = blocks[i]->GetVectorDouble("BETA", "deg"); + int SI = blocks[i]->GetInt("SI"); + int SILI = blocks[i]->GetInt("SILI"); + int CSI = blocks[i]->GetInt("CSI"); + AddTelescope(R, Theta, Phi, beta[0], beta[1], beta[2], SI == 1, SILI == 1, + CSI == 1); } - else{ - cout << "WARNING: Missing token for M2Telescope blocks, check your input file" << endl; + else { + cout << "WARNING: Missing token for M2Telescope blocks, check your input " + "file" + << endl; exit(1); } - if(blocks[i]->GetString("VIS")=="all") + if (blocks[i]->GetString("VIS") == "all") m_non_sensitive_part_visiualisation = true; - } } - // Construct detector and inialise sensitive part. // Called After DetecorConstruction::AddDetector Method -void MUST2Array::ConstructDetector(G4LogicalVolume* world){ - G4RotationMatrix* MMrot = NULL ; - G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ; - bool Si = true ; - bool SiLi = true ; - bool CsI = true ; - - G4int NumberOfTelescope = m_DefinitionType.size() ; - - for (G4int i = 0 ; i < NumberOfTelescope ; i++) { +void MUST2Array::ConstructDetector(G4LogicalVolume* world) { + G4RotationMatrix* MMrot = NULL; + G4ThreeVector MMpos = G4ThreeVector(0, 0, 0); + G4ThreeVector MMu = G4ThreeVector(0, 0, 0); + G4ThreeVector MMv = G4ThreeVector(0, 0, 0); + G4ThreeVector MMw = G4ThreeVector(0, 0, 0); + G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0); + bool Si = true; + bool SiLi = true; + bool CsI = true; + + G4int NumberOfTelescope = m_DefinitionType.size(); + + for (G4int i = 0; i < NumberOfTelescope; i++) { // By Point if (m_DefinitionType[i]) { // (u,v,w) unitary vector associated to telescope referencial // (u,v) // to silicon plan // w perpendicular to (u,v) plan and pointing CsI - MMu = m_X128_Y1[i] - m_X1_Y1[i] ; - MMu = MMu.unit() ; + MMu = m_X128_Y1[i] - m_X1_Y1[i]; + MMu = MMu.unit(); - MMv = m_X1_Y128[i] - m_X1_Y1[i] ; - MMv = MMv.unit() ; + MMv = m_X1_Y128[i] - m_X1_Y1[i]; + MMv = MMv.unit(); - MMw = MMv.cross(MMu) ; + MMw = MMv.cross(MMu); // if (MMw.z() > 0)MMw = MMv.cross(MMu) ; - MMw = MMw.unit() ; + MMw = MMw.unit(); - MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4 ; + MMCenter + = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4; // Passage Matrix from Lab Referential to Telescope Referential MMrot = new G4RotationMatrix(MMv, MMu, MMw); - MMpos = MMw * Length * 0.5 + MMCenter ; + MMpos = MMw * Length * 0.5 + MMCenter; } // By Angle else { - G4double Theta = m_Theta[i] ; - G4double Phi = m_Phi[i] ; + G4double Theta = m_Theta[i]; + G4double Phi = m_Phi[i]; // (u,v,w) unitary vector associated to telescope referencial // (u,v) // to silicon plan // w perpendicular to (u,v) plan and pointing ThirdStage // Phi is angle between X axis and projection in (X,Y) plan // Theta is angle between position vector and z axis - G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad) ; - G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad) ; + G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad); + G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad); G4double wZ = m_R[i] * cos(Theta / rad); - MMw = G4ThreeVector(wX, wY, wZ) ; + MMw = G4ThreeVector(wX, wY, wZ); // vector corresponding to the center of the module G4ThreeVector CT = MMw; // vector parallel to one axis of silicon plane - G4double ii = cos(Theta / rad) * cos(Phi / rad); - G4double jj = cos(Theta / rad) * sin(Phi / rad); - G4double kk = -sin(Theta / rad); - G4ThreeVector Y = G4ThreeVector(ii, jj, kk); + G4double ii = cos(Theta / rad) * cos(Phi / rad); + G4double jj = cos(Theta / rad) * sin(Phi / rad); + G4double kk = -sin(Theta / rad); + G4ThreeVector Y = G4ThreeVector(ii, jj, kk); MMw = MMw.unit(); MMu = MMw.cross(Y); @@ -691,215 +772,232 @@ void MUST2Array::ConstructDetector(G4LogicalVolume* world){ MMrot->rotate(m_beta_v[i], MMv); MMrot->rotate(m_beta_w[i], MMw); // translation to place Telescope - MMpos = MMw * Length * 0.5 + CT ; + MMpos = MMw * Length * 0.5 + CT; } - Si = m_wSi[i] ; - SiLi = m_wSiLi[i] ; - CsI = m_wCsI[i] ; + Si = m_wSi[i]; + SiLi = m_wSiLi[i]; + CsI = m_wCsI[i]; - VolumeMaker(i + 1 , MMpos , MMrot , Si , SiLi , CsI , world); + VolumeMaker(i + 1, MMpos, MMrot, Si, SiLi, CsI, world); } - delete MMrot ; + delete MMrot; } // Add Detector branch to the EventTree. // Called After DetecorConstruction::AddDetector Method -void MUST2Array::InitializeRootOutput(){ - RootOutput *pAnalysis = RootOutput::getInstance(); - TTree *pTree = pAnalysis->GetTree(); - if(!pTree->FindBranch("MUST2")){ - pTree->Branch("MUST2", "TMust2Data", &m_Event) ; - } - pTree->SetBranchAddress("MUST2", &m_Event) ; +void MUST2Array::InitializeRootOutput() { + RootOutput* pAnalysis = RootOutput::getInstance(); + TTree* pTree = pAnalysis->GetTree(); + if (!pTree->FindBranch("MUST2")) { + pTree->Branch("MUST2", "TMust2Data", &m_Event); + } + pTree->SetBranchAddress("MUST2", &m_Event); } // Read sensitive part and fill the Root tree. // Called at in the EventAction::EndOfEventAvtion -void MUST2Array::ReadSensitive(const G4Event*){ +void MUST2Array::ReadSensitive(const G4Event*) { m_Event->Clear(); ////////////////////////////////////////////////////////////////////////////////////// - //////////////////////// Used to Read Event Map of detector ////////////////////////// + //////////////////////// Used to Read Event Map of detector + ///////////////////////////// ////////////////////////////////////////////////////////////////////////////////////// ///////////////////// // Read the Scorer associate to the Silicon Strip - DSSDScorers::PS_Images* SiScorer = (DSSDScorers::PS_Images*) m_StripScorer->GetPrimitive(0); + DSSDScorers::PS_Images* SiScorer + = (DSSDScorers::PS_Images*)m_StripScorer->GetPrimitive(0); - bool SiScoredHit; // flag true if first stage scores a hit above threshold + bool SiScoredHit; // flag true if first stage scores a hit above threshold set<int> trig; // list of telescope that got a Si trigger unsigned int sizeFront = SiScorer->GetFrontMult(); unsigned int sizeBack = SiScorer->GetBackMult(); - for(unsigned int i = 0 ; i < sizeFront ; i++){ - double energy = SiScorer->GetEnergyFront(i); - double energyX = RandGauss::shoot(energy, ResoStrip); - int detectorNbr = SiScorer->GetDetectorFront(i); - double time = SiScorer->GetTimeFront(i); + for (unsigned int i = 0; i < sizeFront; i++) { + double energy = SiScorer->GetEnergyFront(i); + double energyX = RandGauss::shoot(energy, ResoStrip); + int detectorNbr = SiScorer->GetDetectorFront(i); + double time = SiScorer->GetTimeFront(i); // X - if(energyX>0.1*keV){ // above threshold - SiScoredHit=true; + if (energyX > 0.1 * keV) { // above threshold + SiScoredHit = true; // Pixel value at interaction point - unsigned int a,r,g,b; + unsigned int a, r, g, b; // pixel - SiScorer->GetARGBFront(i,a,r,g,b); - b=b+2; - g=g+2; - if(r==0){ + SiScorer->GetARGBFront(i, a, r, g, b); + b = b + 2; + g = g + 2; + if (r == 0) { trig.insert(detectorNbr); - // Energy - m_Event->SetStripXE(detectorNbr,b+1,NPL::EnergyToADC(energyX,0,63,8192,16384)); - // Time + // Energy + m_Event->SetStripXE(detectorNbr, b + 1, + NPL::EnergyToADC(energyX, 0, 63, 8192, 16384)); + // Time double timeX = TimeOffset - RandGauss::shoot(time, ResoTimeMust); - m_Event->SetStripXT(detectorNbr,b+1,NPL::EnergyToADC(timeX,0,1000,8192,16384)); - } - else{ // Interstrip X, keep maximum shared energy + m_Event->SetStripXT(detectorNbr, b + 1, + NPL::EnergyToADC(timeX, 0, 1000, 8192, 16384)); + } else { // Interstrip X, keep maximum shared energy double rand = G4UniformRand(); - if(rand>0.5){ - energyX = rand*energyX; - if(energyX>0.1*keV){ + if (rand > 0.5) { + energyX = rand * energyX; + if (energyX > 0.1 * keV) { trig.insert(detectorNbr); - // Energy - m_Event->SetStripXE(detectorNbr,b+1,NPL::EnergyToADC(energyX,0,63,8192,16384)) ; - // Time + // Energy + m_Event->SetStripXE(detectorNbr, b + 1, + NPL::EnergyToADC(energyX, 0, 63, 8192, 16384)); + // Time double timeX = TimeOffset - RandGauss::shoot(time, ResoTimeMust); - m_Event->SetStripXT(detectorNbr,b+1,NPL::EnergyToADC(timeX,0,1000,8192,16384)); - } - } - else{ - energyX = (1-rand)*energyX; - if(energyX>0.1*keV){ + m_Event->SetStripXT(detectorNbr, b + 1, + NPL::EnergyToADC(timeX, 0, 1000, 8192, 16384)); + } + } else { + energyX = (1 - rand) * energyX; + if (energyX > 0.1 * keV) { trig.insert(detectorNbr); - // Energy - m_Event->SetStripXE(detectorNbr,g+1,NPL::EnergyToADC(energyX,0,63,8192,16384)) ; - // Time + // Energy + m_Event->SetStripXE(detectorNbr, g + 1, + NPL::EnergyToADC(energyX, 0, 63, 8192, 16384)); + // Time double timeX = TimeOffset - RandGauss::shoot(time, ResoTimeMust); - m_Event->SetStripXT(detectorNbr,g+1,NPL::EnergyToADC(timeX,0,1000,8192,16384)); + m_Event->SetStripXT(detectorNbr, g + 1, + NPL::EnergyToADC(timeX, 0, 1000, 8192, 16384)); } } } } }; - for(unsigned int i = 0 ; i < sizeBack ; i++){ - double energy = SiScorer->GetEnergyBack(i); - double energyY = RandGauss::shoot(energy, ResoStrip); - int detectorNbr = SiScorer->GetDetectorBack(i); - double time = SiScorer->GetTimeBack(i); + for (unsigned int i = 0; i < sizeBack; i++) { + double energy = SiScorer->GetEnergyBack(i); + double energyY = RandGauss::shoot(energy, ResoStrip); + int detectorNbr = SiScorer->GetDetectorBack(i); + double time = SiScorer->GetTimeBack(i); // Y - if(energyY>0.1*keV){ // above threshold - SiScoredHit=true; + if (energyY > 0.1 * keV) { // above threshold + SiScoredHit = true; // Pixel value at interaction point - unsigned int a,r,g,b; + unsigned int a, r, g, b; // pixel - SiScorer->GetARGBBack(i,a,r,g,b); - b=b+2; - g=g+2; - if(r==0){ + SiScorer->GetARGBBack(i, a, r, g, b); + b = b + 2; + g = g + 2; + if (r == 0) { trig.insert(detectorNbr); - // Energy - m_Event->SetStripYE(detectorNbr,b+1,NPL::EnergyToADC(energyY,0,63,8192,0)); - // Time + // Energy + m_Event->SetStripYE(detectorNbr, b + 1, + NPL::EnergyToADC(energyY, 0, 63, 8192, 0)); + // Time double timeY = TimeOffset - RandGauss::shoot(time, ResoTimeMust); - m_Event->SetStripYT(detectorNbr,b+1,NPL::EnergyToADC(timeY,0,1000,8192,16384)); - } - else{ // Interstrip Y, keep both strip with shared energy - double rand = G4UniformRand(); - double energyY1 = rand*energyY; - if(energyY1>0.1*keV){ + m_Event->SetStripYT(detectorNbr, b + 1, + NPL::EnergyToADC(timeY, 0, 1000, 8192, 16384)); + } else { // Interstrip Y, keep both strip with shared energy + double rand = G4UniformRand(); + double energyY1 = rand * energyY; + if (energyY1 > 0.1 * keV) { trig.insert(detectorNbr); - // Energy - m_Event->SetStripYE(detectorNbr,b+1,NPL::EnergyToADC(energyY1,0,63,8192,0)); + // Energy + m_Event->SetStripYE(detectorNbr, b + 1, + NPL::EnergyToADC(energyY1, 0, 63, 8192, 0)); - // Time + // Time double timeY = TimeOffset - RandGauss::shoot(time, ResoTimeMust); - m_Event->SetStripYT(detectorNbr,b+1,NPL::EnergyToADC(timeY,0,1000,8192,16384)); - } + m_Event->SetStripYT(detectorNbr, b + 1, + NPL::EnergyToADC(timeY, 0, 1000, 8192, 16384)); + } - if(energyY1>0.1*keV){ + if (energyY1 > 0.1 * keV) { trig.insert(detectorNbr); - double energyY2 = (1-rand)*energyY; - // Energy - m_Event->SetStripYE(detectorNbr,g+1,NPL::EnergyToADC(energyY2,0,63,8192,0)); - // Time + double energyY2 = (1 - rand) * energyY; + // Energy + m_Event->SetStripYE(detectorNbr, g + 1, + NPL::EnergyToADC(energyY2, 0, 63, 8192, 0)); + // Time double timeY = TimeOffset - RandGauss::shoot(time, ResoTimeMust); - m_Event->SetStripYT(detectorNbr,g+1,NPL::EnergyToADC(timeY,0,1000,8192,16384)); - } + m_Event->SetStripYT(detectorNbr, g + 1, + NPL::EnergyToADC(timeY, 0, 1000, 8192, 16384)); + } } } } // Look for 2nd and 3rd stage only if 1st stage is hit - if(SiScoredHit){ + if (SiScoredHit) { // SiLi // - CalorimeterScorers::PS_Calorimeter* SiLiScorer= (CalorimeterScorers::PS_Calorimeter*) m_SiLiScorer->GetPrimitive(0); - - unsigned int sizeSiLi = SiLiScorer->GetMult(); - for(unsigned int i = 0 ; i < sizeSiLi ; i++){ - double ESiLi = RandGauss::shoot(SiLiScorer->GetEnergy(i),ResoSiLi); - vector<unsigned int> level = SiLiScorer->GetLevel(i); - m_Event->SetSiLiE(level[0],level[1],NPL::EnergyToADC(ESiLi,0,250,8192,16384)); - double timeSiLi = RandGauss::shoot(SiLiScorer->GetTime(i),ResoTimeMust); - m_Event->SetSiLiT(level[0],level[1],NPL::EnergyToADC(timeSiLi,0,1000,16384,8192)); - - } - + CalorimeterScorers::PS_Calorimeter* SiLiScorer + = (CalorimeterScorers::PS_Calorimeter*)m_SiLiScorer->GetPrimitive(0); + + unsigned int sizeSiLi = SiLiScorer->GetMult(); + for (unsigned int i = 0; i < sizeSiLi; i++) { + double ESiLi = RandGauss::shoot(SiLiScorer->GetEnergy(i), ResoSiLi); + vector<unsigned int> level = SiLiScorer->GetLevel(i); + m_Event->SetSiLiE(level[0], level[1], + NPL::EnergyToADC(ESiLi, 0, 250, 8192, 16384)); + double timeSiLi = RandGauss::shoot(SiLiScorer->GetTime(i), ResoTimeMust); + m_Event->SetSiLiT(level[0], level[1], + NPL::EnergyToADC(timeSiLi, 0, 1000, 16384, 8192)); + } + // CsI // - CalorimeterScorers::PS_Calorimeter* CsIScorer= (CalorimeterScorers::PS_Calorimeter*) m_CsIScorer->GetPrimitive(0); - - unsigned int sizeCsI = CsIScorer->GetMult(); - for(unsigned int i = 0 ; i < sizeCsI ; i++){ - double ECsI = RandGauss::shoot(CsIScorer->GetEnergy(i),ResoCsI); - vector<unsigned int> level = CsIScorer->GetLevel(i); - m_Event->SetCsIE(level[0],level[1],NPL::EnergyToADC(ECsI,0,250,8192,16384)); - double timeCsI = RandGauss::shoot(CsIScorer->GetTime(i),ResoTimeMust); - m_Event->SetCsIT(level[0],level[1],NPL::EnergyToADC(timeCsI,0,1000,16384,8192)); - } + CalorimeterScorers::PS_Calorimeter* CsIScorer + = (CalorimeterScorers::PS_Calorimeter*)m_CsIScorer->GetPrimitive(0); + + unsigned int sizeCsI = CsIScorer->GetMult(); + for (unsigned int i = 0; i < sizeCsI; i++) { + double ECsI = RandGauss::shoot(CsIScorer->GetEnergy(i), ResoCsI); + vector<unsigned int> level = CsIScorer->GetLevel(i); + m_Event->SetCsIE(level[0], level[1], + NPL::EnergyToADC(ECsI, 0, 250, 8192, 16384)); + double timeCsI = RandGauss::shoot(CsIScorer->GetTime(i), ResoTimeMust); + m_Event->SetCsIT(level[0], level[1], + NPL::EnergyToADC(timeCsI, 0, 1000, 16384, 8192)); + } } } - //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void MUST2Array::InitializeScorers() { +void MUST2Array::InitializeScorers() { // Silicon Associate Scorer - bool already_exist = false; - m_StripScorer = CheckScorer("MUST2_StripScorer",already_exist); - m_SiLiScorer = CheckScorer("MUST2_SiLiScorer",already_exist); - m_CsIScorer = CheckScorer("MUST2_CsIScorer",already_exist); + bool already_exist = false; + m_StripScorer = CheckScorer("MUST2_StripScorer", already_exist); + m_SiLiScorer = CheckScorer("MUST2_SiLiScorer", already_exist); + m_CsIScorer = CheckScorer("MUST2_CsIScorer", already_exist); // if the scorer were created previously nothing else need to be made - if(already_exist) return; + if (already_exist) + return; - string nptool = getenv("NPTOOL"); - G4VPrimitiveScorer* SiScorer = - new DSSDScorers::PS_Images("SiScorer",nptool+"/NPLib/Detectors/MUST2/ressources/maskFront.png",nptool+"/NPLib/Detectors/MUST2/ressources/maskBack.png",0.01,0.01,0,0,0xffff0000,0); + string nptool = getenv("NPTOOL"); + G4VPrimitiveScorer* SiScorer = new DSSDScorers::PS_Images( + "SiScorer", nptool + "/NPLib/Detectors/MUST2/ressources/maskFront.png", + nptool + "/NPLib/Detectors/MUST2/ressources/maskBack.png", 0.01, 0.01, 0, + 0, 0xffff0000, 0); - G4VPrimitiveScorer* InterScorer = new InteractionScorers::PS_Interactions("SiScorer",ms_InterCoord,0); + G4VPrimitiveScorer* InterScorer + = new InteractionScorers::PS_Interactions("SiScorer", ms_InterCoord, 0); - - //and register it to the multifunctionnal detector + // and register it to the multifunctionnal detector m_StripScorer->RegisterPrimitive(SiScorer); m_StripScorer->RegisterPrimitive(InterScorer); - // SiLi Associate Scorer - vector<int> SiLi_nesting={3,0}; - G4VPrimitiveScorer* SiLiScorer= - new CalorimeterScorers::PS_Calorimeter("SiLiScorer",SiLi_nesting) ; + vector<int> SiLi_nesting = {3, 0}; + G4VPrimitiveScorer* SiLiScorer + = new CalorimeterScorers::PS_Calorimeter("SiLiScorer", SiLi_nesting); m_SiLiScorer->RegisterPrimitive(SiLiScorer); - // CsI Associate Scorer - vector<int> CsI_nesting = {2,0}; - G4VPrimitiveScorer* CsIScorer= - new CalorimeterScorers::PS_Calorimeter("CsIScorer",CsI_nesting, 0) ; - m_CsIScorer->RegisterPrimitive(CsIScorer) ; + // CsI Associate Scorer + vector<int> CsI_nesting = {2, 0}; + G4VPrimitiveScorer* CsIScorer + = new CalorimeterScorers::PS_Calorimeter("CsIScorer", CsI_nesting, 0); + m_CsIScorer->RegisterPrimitive(CsIScorer); // Add All Scorer to the Global Scorer Manager G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer); @@ -908,53 +1006,57 @@ void MUST2Array::InitializeScorers() { } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -void MUST2Array::InitializeMaterial(){ +void MUST2Array::InitializeMaterial() { - m_MaterialSilicon = MaterialManager::getInstance()->GetMaterialFromLibrary("Si"); - m_MaterialAluminium = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); + m_MaterialSilicon + = MaterialManager::getInstance()->GetMaterialFromLibrary("Si"); + m_MaterialAluminium + = MaterialManager::getInstance()->GetMaterialFromLibrary("Al"); m_MaterialIron = MaterialManager::getInstance()->GetMaterialFromLibrary("Fe"); - m_MaterialVacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); + m_MaterialVacuum + = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -G4RotationMatrix* Rotation(double tetaX, double tetaY, double tetaZ){ - double PI = 3.141592653589793238; +G4RotationMatrix* Rotation(double tetaX, double tetaY, double tetaZ) { + double PI = 3.141592653589793238; double radX = tetaX * PI / 180.; double radY = tetaY * PI / 180.; double radZ = tetaZ * PI / 180.; G4ThreeVector col1 = G4ThreeVector(cos(radZ) * cos(radY), - -sin(radZ) * cos(radY), - -sin(radY)); - G4ThreeVector col2 = G4ThreeVector(sin(radZ) * cos(radX) - cos(radZ) * sin(radY) * sin(radX), - cos(radZ) * cos(radX) + sin(radZ) * sin(radY) * sin(radX), - -cos(radY) * sin(radX)); - G4ThreeVector col3 = G4ThreeVector(sin(radZ) * sin(radX) + cos(radZ) * sin(radY) * sin(radX), - cos(radZ) * sin(radX) - sin(radZ) * sin(radY) * cos(radX), - cos(radY) * cos(radX)); + -sin(radZ) * cos(radY), -sin(radY)); + G4ThreeVector col2 + = G4ThreeVector(sin(radZ) * cos(radX) - cos(radZ) * sin(radY) * sin(radX), + cos(radZ) * cos(radX) + sin(radZ) * sin(radY) * sin(radX), + -cos(radY) * sin(radX)); + G4ThreeVector col3 + = G4ThreeVector(sin(radZ) * sin(radX) + cos(radZ) * sin(radY) * sin(radX), + cos(radZ) * sin(radX) - sin(radZ) * sin(radY) * cos(radX), + cos(radY) * cos(radX)); return (new G4RotationMatrix(col1, col2, col3)); } - //////////////////////////////////////////////////////////////////////////////// // Construct Method to be pass to the DetectorFactory // //////////////////////////////////////////////////////////////////////////////// -NPS::VDetector* MUST2Array::Construct(){ - return (NPS::VDetector*) new MUST2Array(); +NPS::VDetector* MUST2Array::Construct() { + return (NPS::VDetector*)new MUST2Array(); } //////////////////////////////////////////////////////////////////////////////// // Registering the construct method to the factory // //////////////////////////////////////////////////////////////////////////////// -extern"C" { - class proxy_nps_must2{ - public: - proxy_nps_must2(){ - NPS::DetectorFactory::getInstance()->AddToken("M2Telescope","MUST2"); - NPS::DetectorFactory::getInstance()->AddDetector("M2Telescope",MUST2Array::Construct); - } - }; +extern "C" { +class proxy_nps_must2 { +public: + proxy_nps_must2() { + NPS::DetectorFactory::getInstance()->AddToken("M2Telescope", "MUST2"); + NPS::DetectorFactory::getInstance()->AddDetector("M2Telescope", + MUST2Array::Construct); + } +}; - proxy_nps_must2 p_nps_must2; +proxy_nps_must2 p_nps_must2; } diff --git a/NPSimulation/EventGenerator/EventGeneratorCosmic.cc b/NPSimulation/EventGenerator/EventGeneratorCosmic.cc index af76a85dcd29508ab8c040edd7069f8112d32864..fa19c0ef8a39fee0ea66a125e5468844495a074e 100644 --- a/NPSimulation/EventGenerator/EventGeneratorCosmic.cc +++ b/NPSimulation/EventGenerator/EventGeneratorCosmic.cc @@ -136,8 +136,8 @@ void EventGeneratorCosmic::ReadConfiguration(NPL::InputParser parser){ G4double CosmicAngle; - G4double H = 1500*mm; - G4double R = 500*mm; + G4double H = 1500.*mm; + G4double R = 500.*mm; G4double x0 =0; G4double m_y0 = 0; G4double z0 = 0; @@ -148,7 +148,7 @@ void EventGeneratorCosmic::ReadConfiguration(NPL::InputParser parser){ G4double angle = 0; //<<<<<<< HEAD - TF1* cosSq= new TF1("cosSq", "TMath::Power(cos(x),2)", 0, (TMath::Pi())/2); + TF1* cosSq= new TF1("cosSq", "TMath::Power(cos(x),2)", 0, (TMath::Pi())/2.); //TH1F* DistribCosmicAngle= new TH1F("DistribCosmicAngle","DistribCosmicAngle",100, 0, (TMath::Pi())/2); //TH2F* DistribMomZMomX= new TH2F("DistribMomZMomX","Horizontals components of Cosmic rays Momentums",50, -1, 1, 50, -1, 1); //TCanvas* DisCanva = new TCanvas("DisCanva","Distribution"); @@ -185,7 +185,7 @@ void EventGeneratorCosmic::GenerateEvent(G4Event*){ - angle = RandFlat::shoot()*2*pi; + angle = RandFlat::shoot()*2.*pi; CosmicAngle = cosSq->GetRandom(); randomize1 = RandFlat::shoot() ; randomize2 = RandFlat::shoot() ; @@ -211,16 +211,16 @@ void EventGeneratorCosmic::GenerateEvent(G4Event*){ // Begin Constrain to pass in a square with L = 3* R - x0 = R*(randomize1-0.5)*3; - z0 = R*(randomize2-0.5)*3; + x0 = R*(randomize1-0.5)*3.; + z0 = R*(randomize2-0.5)*3.; momentum_y = cos(CosmicAngle); momentum_x = cos(angle)*sin(CosmicAngle); momentum_z = sin(angle)*sin(CosmicAngle); - x0 = x0-momentum_x*( H/2 / momentum_y);// *( H/2 / momentum_y) this is to have the origin always with par.m_y0 = H/2; - z0 = z0-momentum_z*( H/2 / momentum_y); - par.m_y0 = H/2; // momentum_y*( H/2 / momentum_y); + x0 = x0-momentum_x*( H/2. / momentum_y);// *( H/2 / momentum_y) this is to have the origin always with par.m_y0 = H/2; + z0 = z0-momentum_z*( H/2. / momentum_y); + par.m_y0 = H/2.; // momentum_y*( H/2 / momentum_y); // End Constrain to pass in a square with L = 3* R diff --git a/Projects/Dali/Dali.detector b/Projects/Dali/Dali.detector index faab7e259bd3400da50a7e6c33143b4879f9da45..08bb567942278e6606bc93264f792f765c6de864 100644 --- a/Projects/Dali/Dali.detector +++ b/Projects/Dali/Dali.detector @@ -1,12 +1,12 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Target - THICKNESS= 20 mm - RADIUS= 20 mm - MATERIAL= CD2 - ANGLE= 0 deg - X= 0 mm - Y= 0 mm - Z= 220 mm +%Target +% THICKNESS= 20 mm +% RADIUS= 20 mm +% MATERIAL= CD2 +% ANGLE= 0 deg +% X= 0 mm +% Y= 0 mm +% Z= 220 mm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Dali R = 212.4 mm diff --git a/Projects/Dali/protonSingle.reaction b/Projects/Dali/protonSingle.reaction index cf1f63939d811b52d92832c94f2104af8d755f2e..6e5ead539b6db5101109689fc2322131c53be365 100644 --- a/Projects/Dali/protonSingle.reaction +++ b/Projects/Dali/protonSingle.reaction @@ -11,7 +11,7 @@ Beam %0.5 SigmaY= 0.01 %0.5 - MeanThetaX= 5 + MeanThetaX= 0 MeanPhiY= 0 MeanX= 0 MeanY= 0