rtang/SOLARIS_Analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Compare commits

base: rtang:3f225316980f5e7e51a78e14d07b0bb6cda120cf
Branches Tags
rtang:Master
rtang:NUMEN
..
compare: rtang:806a7124a4f3ed5c260fbea21e0f1a77af66902f
Branches Tags
rtang:Master
rtang:NUMEN

No commits in common. "3f225316980f5e7e51a78e14d07b0bb6cda120cf" and "806a7124a4f3ed5c260fbea21e0f1a77af66902f" have entirely different histories.
3f22531698 ... 806a7124a4

6 changed files with 138 additions and 145 deletions
Show Stats Expand all files Collapse all files

144
Armory/ClassDetGeo.h
View File

@ -25,7 +25,7 @@ struct Array{
double eSigma; /// intrinsic energy resolution MeV double eSigma; /// intrinsic energy resolution MeV
double zSigma; /// intrinsic position resolution mm double zSigma; /// intrinsic position resolution mm
bool detFaceOut; ///detector_facing_Out_or_In bool detFaceOut; ///detector_facing_Out_or_In
std::vector<double> pos; /// realtive position in meter std::vector<double> pos; /// near position in meter
int nDet, mDet; /// nDet = number of different pos, mDet, number of same pos int nDet, mDet; /// nDet = number of different pos, mDet, number of same pos
std::vector<double> detPos; ///absolute position of detector std::vector<double> detPos; ///absolute position of detector
@ -46,10 +46,7 @@ struct Array{
zMax = TMath::Max(detPos.front(), detPos.back()) + (firstPos > 0 ? detLength : 0); zMax = TMath::Max(detPos.front(), detPos.back()) + (firstPos > 0 ? detLength : 0);
} }
void Print() const{ void PrintArray() const{
printf("------------------------------- Array\n");
for(int i = 0; i < nDet ; i++){ for(int i = 0; i < nDet ; i++){
if( firstPos > 0 ){ if( firstPos > 0 ){
printf("%d, %8.2f mm - %8.2f mm \n", i, detPos[i], detPos[i] + detLength); printf("%d, %8.2f mm - %8.2f mm \n", i, detPos[i], detPos[i] + detLength);
@ -60,7 +57,7 @@ struct Array{
printf(" Blocker Position: %8.2f mm \n", firstPos > 0 ? firstPos - blocker : firstPos + blocker ); printf(" Blocker Position: %8.2f mm \n", firstPos > 0 ? firstPos - blocker : firstPos + blocker );
printf(" First Position: %8.2f mm \n", firstPos); printf(" First Position: %8.2f mm \n", firstPos);
printf(" number of det : %d x %d (side x col) \n", mDet, nDet); printf(" number of det : %d x %d \n", mDet, nDet);
printf(" detector facing : %s\n", detFaceOut ? "Out" : "In"); printf(" detector facing : %s\n", detFaceOut ? "Out" : "In");
printf(" energy resol.: %f MeV\n", eSigma); printf(" energy resol.: %f MeV\n", eSigma);
printf(" pos-Z resol.: %f mm \n", zSigma); printf(" pos-Z resol.: %f mm \n", zSigma);
@ -68,38 +65,6 @@ struct Array{
}; };
struct Auxillary{
//========== recoil
double detPos; //
double outerRadius;
double innerRadius;
bool isCoincident;
double detPos1 = 0; // virtual recoil
double detPos2 = 0; // virtual recoil
//========== enum
double elumPos1 = 0;
double elumPos2 = 0;
void Print() const {
printf("------------------------------- Auxillary\n");
printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", detPos, innerRadius, outerRadius);
if( elumPos1 != 0 || elumPos2 != 0 || detPos1 != 0 || detPos2 != 0){
printf("=================================== Virtual Detectors\n");
if( elumPos1 != 0 ) printf(" Elum 1 pos.: %f mm \n", elumPos1);
if( elumPos2 != 0 ) printf(" Elum 2 pos.: %f mm \n", elumPos2);
if( detPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", detPos1);
if( detPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", detPos2);
printf("=====================================================\n");
}
}
};
class DetGeo { class DetGeo {
public: public:
@ -110,18 +75,28 @@ public:
double Bfield; /// T double Bfield; /// T
int BfieldSign ; /// sign of B-field int BfieldSign ; /// sign of B-field
double BfieldTheta; /// rad, 0 = z-axis, pi/2 = y axis, pi = -z axis
double bore; /// bore , mm double bore; /// bore , mm
double zMin, zMax; /// total range span of all arrays double recoilPos; /// recoil, downstream
double recoilInnerRadius; /// radius recoil inner
double recoilOuterRadius; /// radius recoil outter
double recoilPos1, recoilPos2; /// imaginary recoils
double elumPos1, elumPos2; /// imaginary elum, only sensitive to light recoil
//===================1st array
Array array[2];
double zMin, zMax; /// range of detectors
bool isCoincidentWithRecoil;
bool LoadDetectorGeo(TString fileName, bool verbose = true); bool LoadDetectorGeo(TString fileName, bool verbose = true);
bool LoadDetectorGeo(TMacro * macro, bool verbose = true); bool LoadDetectorGeo(TMacro * macro, bool verbose = true);
//=================== array void PrintWithoutArray() ;
std::vector<Array> array; void Print( bool printAll = true) ;
std::vector<Auxillary> aux;
void Print( bool printAll = false) ;
private: private:
@ -154,7 +129,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
TList * haha = macro->GetListOfLines(); TList * haha = macro->GetListOfLines();
int numLine = (haha)->GetSize(); int numLine = (haha)->GetSize();
// for( int i = 0; i < 2 ; i++) array[i].pos.clear(); for( int i = 0; i < 2 ; i++) array[i].pos.clear();
int detFlag = 0; int detFlag = 0;
int detLine = 0; int detLine = 0;
@ -172,8 +147,6 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
if( str[0].find("####") != std::string::npos ) break; if( str[0].find("####") != std::string::npos ) break;
if( str[0].find("#===") != std::string::npos ) { if( str[0].find("#===") != std::string::npos ) {
detFlag ++; detFlag ++;
array.push_back(Array());
aux.push_back(Auxillary());
detLine = 0; detLine = 0;
continue;; continue;;
} }
@ -183,30 +156,31 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
Bfield = atof(str[0].c_str()); Bfield = atof(str[0].c_str());
BfieldSign = Bfield > 0 ? 1: -1; BfieldSign = Bfield > 0 ? 1: -1;
} }
if ( detLine == 1 ) bore = atof(str[0].c_str()); if ( detLine == 1 ) BfieldTheta = atof(str[0].c_str());
if ( detLine == 2 ) bore = atof(str[0].c_str());
if ( detLine == 3 ) recoilPos = atof(str[0].c_str());
if ( detLine == 4 ) recoilInnerRadius = atof(str[0].c_str());
if ( detLine == 5 ) recoilOuterRadius = atof(str[0].c_str());
if ( detLine == 6 ) isCoincidentWithRecoil = str[0] == "false" ? false: true;
if ( detLine == 7 ) recoilPos1 = atof(str[0].c_str());
if ( detLine == 8 ) recoilPos2 = atof(str[0].c_str());
if ( detLine == 9 ) elumPos1 = atof(str[0].c_str());
if ( detLine == 10 ) elumPos2 = atof(str[0].c_str());
} }
if( detFlag > 0){ if( detFlag > 0){
unsigned short ID = detFlag - 1; unsigned short ID = detFlag - 1;
if ( detLine == 0 ) array[ID].enable = str[0] == "true" ? true : false; if ( detLine == 0 ) array[ID].enable = str[0] == "true" ? true : false;
if ( detLine == 1 ) aux[ID].detPos = atof(str[0].c_str()); if ( detLine == 1 ) array[ID].detPerpDist = atof(str[0].c_str());
if ( detLine == 2 ) aux[ID].innerRadius = atof(str[0].c_str()); if ( detLine == 2 ) array[ID].detWidth = atof(str[0].c_str());
if ( detLine == 3 ) aux[ID].outerRadius = atof(str[0].c_str()); if ( detLine == 3 ) array[ID].detLength = atof(str[0].c_str());
if ( detLine == 4 ) aux[ID].isCoincident = str[0] == "true" ? true : false; if ( detLine == 4 ) array[ID].blocker = atof(str[0].c_str());
if ( detLine == 5 ) aux[ID].detPos1 = atof(str[0].c_str()); if ( detLine == 5 ) array[ID].firstPos = atof(str[0].c_str());
if ( detLine == 6 ) aux[ID].detPos2 = atof(str[0].c_str()); if ( detLine == 6 ) array[ID].eSigma = atof(str[0].c_str());
if ( detLine == 7 ) aux[ID].elumPos1 = atof(str[0].c_str()); if ( detLine == 7 ) array[ID].zSigma = atof(str[0].c_str());
if ( detLine == 8 ) aux[ID].elumPos2 = atof(str[0].c_str()); if ( detLine == 8 ) array[ID].detFaceOut = str[0] == "Out" ? true : false;
if ( detLine == 9 ) array[ID].detPerpDist = atof(str[0].c_str()); if ( detLine == 9 ) array[ID].mDet = atoi(str[0].c_str());
if ( detLine == 10 ) array[ID].detWidth = atof(str[0].c_str()); if ( detLine >= 10 ) array[ID].pos.push_back(atof(str[0].c_str()));
if ( detLine == 11 ) array[ID].detLength = atof(str[0].c_str());
if ( detLine == 12 ) array[ID].blocker = atof(str[0].c_str());
if ( detLine == 13 ) array[ID].firstPos = atof(str[0].c_str());
if ( detLine == 14 ) array[ID].eSigma = atof(str[0].c_str());
if ( detLine == 15 ) array[ID].zSigma = atof(str[0].c_str());
if ( detLine == 16 ) array[ID].detFaceOut = str[0] == "Out" ? true : false;
if ( detLine == 17 ) array[ID].mDet = atoi(str[0].c_str());
if ( detLine >= 18 ) array[ID].pos.push_back(atof(str[0].c_str()));
} }
detLine ++; detLine ++;
@ -215,7 +189,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
zMin = 99999; zMin = 99999;
zMax = -99999; zMax = -99999;
for( int i = 0; i < detFlag; i ++ ){ for( int i = 0; i < 2; i ++ ){
array[i].DeduceAbsolutePos(); array[i].DeduceAbsolutePos();
if (array[i].enable ) { if (array[i].enable ) {
double zmax = TMath::Max(array[i].zMin, array[i].zMax); double zmax = TMath::Max(array[i].zMin, array[i].zMax);
@ -231,23 +205,37 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
} }
inline void DetGeo::PrintWithoutArray(){
printf("=====================================================\n");
printf(" B-field: %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan");
printf(" B-field Theta : %6.2f deg \n", BfieldTheta);
if( BfieldTheta != 0.0 ) printf(" +---- field angle != 0 is not supported!!! \n");
printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", recoilPos, recoilInnerRadius, recoilOuterRadius);
if( elumPos1 != 0 || elumPos2 != 0 || recoilPos1 != 0 || recoilPos2 != 0){
printf("=================================== Auxillary/Imaginary Detectors\n");
}
if( elumPos1 != 0 ) printf(" Elum 1 pos.: %f mm \n", elumPos1);
if( elumPos2 != 0 ) printf(" Elum 2 pos.: %f mm \n", elumPos2);
if( recoilPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", recoilPos1);
if( recoilPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", recoilPos2);
printf("=====================================================\n");
}
inline void DetGeo::Print(bool printAll){ inline void DetGeo::Print(bool printAll){
printf("#####################################################\n"); PrintWithoutArray();
printf(" B-field : %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan");
printf(" Bore : %8.2f mm\n", bore); for( int i = 0; i < 2 ; i++){
printf(" No. of det. Set. : %zu \n", array.size());
printf(" z-Min : %8.2f mm\n", zMin);
printf(" z-Max : %8.2f mm\n", zMax);
for( size_t i = 0; i < array.size() ; i++){
printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled");
if( printAll || array[i].enable ) { if( printAll || array[i].enable ) {
array[i].Print(); printf("-----------------------------------%d-th Detector Position \n", i);
aux[i].Print(); array[i].PrintArray();
} }
} }
printf("#####################################################\n"); printf("=====================================================\n");
} }

27
Armory/ClassReactionConfig.h
View File

@ -70,7 +70,6 @@ struct ExcitedEnergies {
} }
void Print() const { void Print() const {
printf("...................................\n");
printf("Energy[MeV] Rel.Xsec SF sigma\n"); printf("Energy[MeV] Rel.Xsec SF sigma\n");
for( size_t i = 0; i < ExList.size(); i++){ for( size_t i = 0; i < ExList.size(); i++){
ExList[i].Print("\n"); ExList[i].Print("\n");
@ -102,8 +101,8 @@ public:
float targetThickness; ///targetThickness_in_cm float targetThickness; ///targetThickness_in_cm
std::string beamStoppingPowerFile; ///stopping_power_for_beam std::string beamStoppingPowerFile; ///stopping_power_for_beam
std::vector<Recoil> recoil; Recoil recoil[2];
std::vector<ExcitedEnergies> exList; ExcitedEnergies exList[2];
int numEvents; ///number_of_Event_being_generated int numEvents; ///number_of_Event_being_generated
bool isRedo; ///isReDo bool isRedo; ///isReDo
@ -164,6 +163,9 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
if( macro == NULL ) return false; if( macro == NULL ) return false;
exList[0].Clear();
exList[1].Clear();
int recoilFlag = 0; int recoilFlag = 0;
int recoilLine = 0; int recoilLine = 0;
@ -182,8 +184,6 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
if( str[0].find("#===") != std::string::npos ) { if( str[0].find("#===") != std::string::npos ) {
recoilFlag ++; recoilFlag ++;
recoilLine = 0; recoilLine = 0;
recoil.push_back(Recoil());
exList.push_back(ExcitedEnergies());
continue; continue;
} }
@ -230,7 +230,7 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
} }
for( size_t i = 0; i < recoil.size(); i++){ for( int i = 0; i < 2; i++){
recoil[i].heavyA = beamA + targetA - recoil[i].lightA; recoil[i].heavyA = beamA + targetA - recoil[i].lightA;
recoil[i].heavyZ = beamZ + targetZ - recoil[i].lightZ; recoil[i].heavyZ = beamZ + targetZ - recoil[i].lightZ;
} }
@ -239,29 +239,28 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
inline void ReactionConfig::Print(int ID, bool withEx) const{ inline void ReactionConfig::Print(int ID, bool withEx) const{
printf("#####################################################\n"); printf("=====================================================\n");
printf("number of Simulation Events : %d \n", numEvents); printf("number of Simulation Events : %d \n", numEvents);
printf(" is Redo until hit array : %s \n", isRedo ? "Yes" : "No"); printf(" is Redo until hit array : %s \n", isRedo ? "Yes" : "No");
printf("================================= Beam\n"); printf("------------------------------ Beam\n");
printf(" beam : A = %3d, Z = %2d, Ex = %.2f MeV\n", beamA, beamZ, beamEx); printf(" beam : A = %3d, Z = %2d, Ex = %.2f MeV\n", beamA, beamZ, beamEx);
printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", beamEnergy, beamEnergySigma, beamEnergySigma/beamEnergy); printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", beamEnergy, beamEnergySigma, beamEnergySigma/beamEnergy);
printf(" Angle : %.2f +- %.2f mrad\n", beamTheta, beamThetaSigma); printf(" Angle : %.2f +- %.2f mrad\n", beamTheta, beamThetaSigma);
printf(" offset : (x,y) = (%.2f, %.2f) mmm \n", beamX, beamY); printf(" offset : (x,y) = (%.2f, %.2f) mmm \n", beamX, beamY);
printf("================================= Target\n"); printf("------------------------------ Target\n");
printf(" target : A = %3d, Z = %2d \n", targetA, targetZ); printf(" target : A = %3d, Z = %2d \n", targetA, targetZ);
printf(" enable scattering : %s \n", isTargetScattering ? "Yes" : "No"); printf(" is target scattering : %s \n", isTargetScattering ? "Yes" : "No");
if(isTargetScattering){ if(isTargetScattering){
printf(" target density : %.f g/cm3\n", targetDensity); printf(" target density : %.f g/cm3\n", targetDensity);
printf(" thickness : %.f cm\n", targetThickness); printf(" thickness : %.f cm\n", targetThickness);
printf(" beam stopping file : %s \n", beamStoppingPowerFile.c_str()); printf(" beam stopping file : %s \n", beamStoppingPowerFile.c_str());
} }
printf("================================= Number of recoil reactions : %zu\n", recoil.size()); for( int i = 0; i < 2; i ++ ){
for( size_t i = 0; i < recoil.size(); i ++ ){ printf("------------------------------ Recoil-%d\n", i);
printf("------------------------------------------ Recoil-%zu\n", i);
if( ID == i || ID < 0 ){ if( ID == i || ID < 0 ){
recoil[i].Print(); recoil[i].Print();
if( withEx ) exList[i].Print(); if( withEx ) exList[i].Print();
@ -269,7 +268,7 @@ inline void ReactionConfig::Print(int ID, bool withEx) const{
} }
printf("#####################################################\n"); printf("=====================================================\n");
} }
#endif #endif

46
Cleopatra/ClassHelios.h
View File

@ -86,6 +86,7 @@ public:
void SetBeamPosition(double x, double y) { xOff = x; yOff = y;} void SetBeamPosition(double x, double y) { xOff = x; yOff = y;}
void OverrideMagneticField(double BField); void OverrideMagneticField(double BField);
void OverrideMagneticFieldDirection(double BfieldThetaInDeg);
void OverrideFirstPos(double firstPos); void OverrideFirstPos(double firstPos);
void OverrideDetectorDistance(double perpDist); void OverrideDetectorDistance(double perpDist);
void OverrideDetectorFacing(bool isOutside); void OverrideDetectorFacing(bool isOutside);
@ -134,7 +135,6 @@ public:
DetGeo GetDetectorGeometry() const {return detGeo;} DetGeo GetDetectorGeometry() const {return detGeo;}
Array GetArrayGeometry() const {return array;} Array GetArrayGeometry() const {return array;}
Auxillary GetAuxGeometry() const {return aux;}
TString GetHitMessage() {return hitMessage;} TString GetHitMessage() {return hitMessage;}
TString GetAcceptanceMessage() { AcceptanceCodeToMsg(acceptanceCode); return acceptanceMsg;} TString GetAcceptanceMessage() { AcceptanceCodeToMsg(acceptanceCode); return acceptanceMsg;}
@ -145,7 +145,6 @@ private:
DetGeo detGeo; DetGeo detGeo;
Array array; Array array;
Auxillary aux;
trajectory orbitb, orbitB; trajectory orbitb, orbitB;
@ -214,6 +213,10 @@ void HELIOS::OverrideMagneticField(double BField){
this->detGeo.BfieldSign = BField > 0 ? 1: -1; this->detGeo.BfieldSign = BField > 0 ? 1: -1;
} }
void HELIOS::OverrideMagneticFieldDirection(double BfieldThetaInDeg){
this->detGeo.BfieldTheta = BfieldThetaInDeg;
}
void HELIOS::OverrideFirstPos(double firstPos){ void HELIOS::OverrideFirstPos(double firstPos){
overrideFirstPos = true; overrideFirstPos = true;
printf("------ Overriding FirstPosition to : %8.2f mm \n", firstPos); printf("------ Overriding FirstPosition to : %8.2f mm \n", firstPos);
@ -236,8 +239,7 @@ bool HELIOS::SetDetectorGeometry(std::string filename, unsigned short ID){
if( detGeo.LoadDetectorGeo(filename, false)) { if( detGeo.LoadDetectorGeo(filename, false)) {
array = detGeo.array[ID]; array = detGeo.array[ID];
aux = detGeo.aux[ID]; isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil;
isCoincidentWithRecoil = detGeo.aux[ID].isCoincident;
isDetReady = true; isDetReady = true;
}else{ }else{
@ -251,13 +253,25 @@ bool HELIOS::SetDetectorGeometry(std::string filename, unsigned short ID){
void HELIOS::PrintGeometry() const{ void HELIOS::PrintGeometry() const{
printf("=====================================================\n"); printf("=====================================================\n");
printf(" B-field : %8.2f T, %s\n", detGeo.Bfield, detGeo.Bfield > 0 ? "out of plan" : "into plan"); printf(" B-field: %8.2f T, Theta : %6.2f deg \n", detGeo.Bfield, detGeo.BfieldTheta);
printf(" Bore : %8.2f mm\n", detGeo.bore); if( detGeo.BfieldTheta != 0.0 ) {
printf(" +---- field angle != 0 is not supported!!! \n");
}
printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", detGeo.recoilPos, detGeo.recoilInnerRadius, detGeo.recoilOuterRadius);
printf("----------------------------------- Detector Position \n"); printf("----------------------------------- Detector Position \n");
array.Print(); array.PrintArray();
aux.Print();
if( detGeo.elumPos1 != 0 || detGeo.elumPos2 != 0 || detGeo.recoilPos1 != 0 || detGeo.recoilPos2 != 0){
printf("=================================== Auxillary/Imaginary Detectors\n");
}
if( detGeo.elumPos1 != 0 ) printf(" Elum 1 pos.: %f mm \n", detGeo.elumPos1);
if( detGeo.elumPos2 != 0 ) printf(" Elum 2 pos.: %f mm \n", detGeo.elumPos2);
if( detGeo.recoilPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", detGeo.recoilPos1);
if( detGeo.recoilPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", detGeo.recoilPos2);
printf("=====================================================\n"); printf("=====================================================\n");
} }
TString HELIOS::AcceptanceCodeToMsg(short code ){ TString HELIOS::AcceptanceCodeToMsg(short code ){
@ -311,14 +325,14 @@ int HELIOS::CheckDetAcceptance(){
if( detGeo.bore < 2 * orbitb.rho) { acceptanceCode = -10; return acceptanceCode;} if( detGeo.bore < 2 * orbitb.rho) { acceptanceCode = -10; return acceptanceCode;}
// -14 ========== check particle-B hit radius on recoil dectector // -14 ========== check particle-B hit radius on recoil dectector
if( isCoincidentWithRecoil && orbitB.R > aux.outerRadius ) {acceptanceCode = -14; return acceptanceCode;} if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) {acceptanceCode = -14; return acceptanceCode;}
//if( isCoincidentWithRecoil && (orbitB.R > rhoRecoilout || orbitB.R < rhoRecoilin) ) return -14; //if( isCoincidentWithRecoil && (orbitB.R > rhoRecoilout || orbitB.R < rhoRecoilin) ) return -14;
// -12 ========= check is particle-b was blocked by recoil detector // -12 ========= check is particle-b was blocked by recoil detector
rhoHit = GetR(aux.detPos); rhoHit = GetR(detGeo.recoilPos);
if( orbitb.z > 0 && aux.detPos > 0 && orbitb.z > aux.detPos && rhoHit < aux.outerRadius ) { acceptanceCode = -12; return acceptanceCode;} if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { acceptanceCode = -12; return acceptanceCode;}
if( orbitb.z < 0 && aux.detPos < 0 && orbitb.z < aux.detPos && rhoHit < aux.outerRadius ) { acceptanceCode = -12; return acceptanceCode;} if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { acceptanceCode = -12; return acceptanceCode;}
// -13 ========= not more than 3 loops // -13 ========= not more than 3 loops
if( orbitb.loop > 3 ) {acceptanceCode = -13; return acceptanceCode;} if( orbitb.loop > 3 ) {acceptanceCode = -13; return acceptanceCode;}
@ -538,10 +552,10 @@ int HELIOS::CalRecoilHit(TLorentzVector PB){
CalTrajectoryPara(PB, false); CalTrajectoryPara(PB, false);
orbitB.z = aux.detPos; orbitB.z = detGeo.recoilPos;
orbitB.x = GetRecoilXPos(aux.detPos) ; orbitB.x = GetRecoilXPos(detGeo.recoilPos) ;
orbitB.y = GetRecoilYPos(aux.detPos) ; orbitB.y = GetRecoilYPos(detGeo.recoilPos) ;
orbitB.R = GetRecoilR(aux.detPos); orbitB.R = GetRecoilR(detGeo.recoilPos);
orbitB.effLoop = orbitB.z/orbitB.z0; orbitB.effLoop = orbitB.z/orbitB.z0;
orbitB.t = orbitB.t0 * orbitB.effLoop ; orbitB.t = orbitB.t0 * orbitB.effLoop ;

41
Cleopatra/SimTransfer.C
View File

@ -76,7 +76,6 @@ void Transfer(
DetGeo detGeo = helios.GetDetectorGeometry(); DetGeo detGeo = helios.GetDetectorGeometry();
Array array = helios.GetArrayGeometry(); Array array = helios.GetArrayGeometry();
Auxillary aux = helios.GetAuxGeometry();
ReactionConfig reactConfig = transfer.GetRectionConfig(); ReactionConfig reactConfig = transfer.GetRectionConfig();
Recoil recoil = transfer.GetRecoil(); Recoil recoil = transfer.GetRecoil();
@ -270,14 +269,14 @@ void Transfer(
///in case need ELUM ///in case need ELUM
double xElum1, yElum1, rhoElum1; double xElum1, yElum1, rhoElum1;
if( aux.elumPos1 != 0 ) { if( detGeo.elumPos1 != 0 ) {
tree->Branch("xElum1", &xElum1, "xElum1/D"); tree->Branch("xElum1", &xElum1, "xElum1/D");
tree->Branch("yElum1", &yElum1, "yElum1/D"); tree->Branch("yElum1", &yElum1, "yElum1/D");
tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D"); tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D");
} }
double xElum2, yElum2, rhoElum2; double xElum2, yElum2, rhoElum2;
if( aux.elumPos2 != 0 ) { if( detGeo.elumPos2 != 0 ) {
tree->Branch("xElum2", &xElum2, "xElum2/D"); tree->Branch("xElum2", &xElum2, "xElum2/D");
tree->Branch("yElum2", &yElum2, "yElum2/D"); tree->Branch("yElum2", &yElum2, "yElum2/D");
tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D"); tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D");
@ -285,13 +284,13 @@ void Transfer(
///in case need other recoil detector. ///in case need other recoil detector.
double xRecoil1, yRecoil1, rhoRecoil1; double xRecoil1, yRecoil1, rhoRecoil1;
if( aux.detPos1 != 0 ){ if( detGeo.recoilPos1 != 0 ){
tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D"); tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D");
tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D"); tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D");
tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D"); tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D");
} }
double xRecoil2, yRecoil2, rhoRecoil2; double xRecoil2, yRecoil2, rhoRecoil2;
if( aux.detPos2 != 0 ){ if( detGeo.recoilPos2 != 0 ){
tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D"); tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D");
tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D"); tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D");
tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D"); tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D");
@ -545,15 +544,15 @@ void Transfer(
//ELUM //ELUM
if( aux.elumPos1 != 0 ){ if( detGeo.elumPos1 != 0 ){
xElum1 = helios.GetXPos(aux.elumPos1); xElum1 = helios.GetXPos(detGeo.elumPos1);
yElum1 = helios.GetYPos(aux.elumPos1); yElum1 = helios.GetYPos(detGeo.elumPos1);
rhoElum1 = helios.GetR(aux.elumPos1); rhoElum1 = helios.GetR(detGeo.elumPos1);
} }
if( aux.elumPos2 != 0 ){ if( detGeo.elumPos2 != 0 ){
xElum2 = helios.GetXPos(aux.elumPos2); xElum2 = helios.GetXPos(detGeo.elumPos2);
yElum2 = helios.GetYPos(aux.elumPos2); yElum2 = helios.GetYPos(detGeo.elumPos2);
rhoElum2 = helios.GetR(aux.elumPos2); rhoElum2 = helios.GetR(detGeo.elumPos2);
} }
//Recoil //Recoil
@ -564,15 +563,15 @@ void Transfer(
rhoB = orb_B.rho; rhoB = orb_B.rho;
//other recoil detectors //other recoil detectors
if ( aux.detPos1 != 0 ){ if ( detGeo.recoilPos1 != 0 ){
xRecoil1 = helios.GetRecoilXPos(aux.detPos1); xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1);
yRecoil1 = helios.GetRecoilYPos(aux.detPos1); yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1);
rhoRecoil1 = helios.GetRecoilR(aux.detPos1); rhoRecoil1 = helios.GetRecoilR(detGeo.recoilPos1);
} }
if ( aux.detPos2 != 0 ){ if ( detGeo.recoilPos2 != 0 ){
xRecoil2 = helios.GetRecoilXPos(aux.detPos2); xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2);
yRecoil2 = helios.GetRecoilYPos(aux.detPos2); yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2);
rhoRecoil2 = helios.GetRecoilR(aux.detPos2); rhoRecoil2 = helios.GetRecoilR(detGeo.recoilPos2);
} }
std::pair<double,double> ExThetaCM = transfer.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius()); std::pair<double,double> ExThetaCM = transfer.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius());

17
working/detectorGeo.txt
View File

@ -1,8 +1,8 @@
-3.00 //Bfield_[T] -3.00 //Bfield_[T]
0.00 //Bfield_direction_to_z-axis_[deg]_should_not_use
462.5 //bore_[mm] 462.5 //bore_[mm]
#===============1st_Array + Recoil
true ////is_this_array_exist_or_use_for_Simulation
1000 //recoil_position_+_for_downstream_[mm] 1000 //recoil_position_+_for_downstream_[mm]
10.0 //inner_radius_of_recoil_detector_[mm] 10.0 //inner_radius_of_recoil_detector_[mm]
40.2 //outter_radius_of_recoil_detector_[mm] 40.2 //outter_radius_of_recoil_detector_[mm]
@ -11,6 +11,9 @@ false //is_coincident_with_recoil
0 //Recoil_2_position_[mm] 0 //Recoil_2_position_[mm]
0.00 //Elum_1_position_[mm]_(just_another_recoil_detector_but_for_light_recoil) 0.00 //Elum_1_position_[mm]_(just_another_recoil_detector_but_for_light_recoil)
0.00 //Elum_2_position_[mm]_when_Elum=0_disable_tree_branch 0.00 //Elum_2_position_[mm]_when_Elum=0_disable_tree_branch
#===============1st_Array
true ////is_this_array_exist_or_use_for_Simulation
11.5 //distance_from_axis_[mm] 11.5 //distance_from_axis_[mm]
10.0 //width_of_detector_[mm] 10.0 //width_of_detector_[mm]
50 //length_of_detector_[mm] 50 //length_of_detector_[mm]
@ -27,16 +30,8 @@ Out //detector_facing_Out_or_In
235.8 //5th_det 235.8 //5th_det
294.0 294.0
#===============2nd_Array + Recoil #===============2nd_Array
true //is_this_array_exist_or_use_for_Simulation true //is_this_array_exist_or_use_for_Simulation
1000 //recoil_position_+_for_downstream_[mm]
10.0 //inner_radius_of_recoil_detector_[mm]
40.2 //outter_radius_of_recoil_detector_[mm]
false //is_coincident_with_recoil
0 //Recoil_1_position_[mm]_when_0_disable_tree_branch
0 //Recoil_2_position_[mm]
0.00 //Elum_1_position_[mm]_(just_another_recoil_detector_but_for_light_recoil)
0.00 //Elum_2_position_[mm]_when_Elum=0_disable_tree_branch
11.5 //distance_from_axis_[mm] 11.5 //distance_from_axis_[mm]
10.0 //width_of_detector_[mm] 10.0 //width_of_detector_[mm]
50 //length_of_detector_[mm] 50 //length_of_detector_[mm]

4
working/reactionConfig.txt
View File

@ -22,9 +22,6 @@ false //isTargetScattering
1000000 //number_of_Event_being_generated 1000000 //number_of_Event_being_generated
false //Redo_until_hit_array=all_events_hit_array false //Redo_until_hit_array=all_events_hit_array
#---- Number of Reaction Should equal to number of array in detectorGeo.txt
#---- If not, somethings funny would happens.
#=================reaction_for_1st_Array #=================reaction_for_1st_Array
1 //recoil_light_A 1 //recoil_light_A
1 //recoil-light_Z 1 //recoil-light_Z
@ -53,4 +50,5 @@ false //isDacay
1.000 1.0 1.0 0.01 1.000 1.0 1.0 0.01
2.000 1.0 1.0 0.01 2.000 1.0 1.0 0.01
################## end of file ################## end of file