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

modified SimTransfer for the modified DetGeo and ReactionConfig

Browse Source
This commit is contained in:
Ryan Tang 2024-04-02 14:30:50 -04:00
parent 0905fe73dc
commit 3f22531698
2 changed files with 37 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
Cleopatra/ClassHelios.h
View File

@ -86,7 +86,6 @@ 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);
@ -135,6 +134,7 @@ 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,6 +145,7 @@ private:
DetGeo detGeo; DetGeo detGeo;
Array array; Array array;
Auxillary aux;
trajectory orbitb, orbitB; trajectory orbitb, orbitB;
@ -213,10 +214,6 @@ 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);
@ -239,7 +236,8 @@ 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];
isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil; aux = detGeo.aux[ID];
isCoincidentWithRecoil = detGeo.aux[ID].isCoincident;
isDetReady = true; isDetReady = true;
}else{ }else{
@ -253,25 +251,13 @@ 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, Theta : %6.2f deg \n", detGeo.Bfield, detGeo.BfieldTheta); printf(" B-field : %8.2f T, %s\n", detGeo.Bfield, detGeo.Bfield > 0 ? "out of plan" : "into plan");
if( detGeo.BfieldTheta != 0.0 ) { printf(" Bore : %8.2f mm\n", detGeo.bore);
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.PrintArray(); array.Print();
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 ){
@ -325,14 +311,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 > detGeo.recoilOuterRadius ) {acceptanceCode = -14; return acceptanceCode;} if( isCoincidentWithRecoil && orbitB.R > aux.outerRadius ) {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(detGeo.recoilPos); rhoHit = GetR(aux.detPos);
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;} if( orbitb.z < 0 && aux.detPos < 0 && orbitb.z < aux.detPos && rhoHit < aux.outerRadius ) { 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;}
@ -552,10 +538,10 @@ int HELIOS::CalRecoilHit(TLorentzVector PB){
CalTrajectoryPara(PB, false); CalTrajectoryPara(PB, false);
orbitB.z = detGeo.recoilPos; orbitB.z = aux.detPos;
orbitB.x = GetRecoilXPos(detGeo.recoilPos) ; orbitB.x = GetRecoilXPos(aux.detPos) ;
orbitB.y = GetRecoilYPos(detGeo.recoilPos) ; orbitB.y = GetRecoilYPos(aux.detPos) ;
orbitB.R = GetRecoilR(detGeo.recoilPos); orbitB.R = GetRecoilR(aux.detPos);
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,6 +76,7 @@ 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();
@ -269,14 +270,14 @@ void Transfer(
///in case need ELUM ///in case need ELUM
double xElum1, yElum1, rhoElum1; double xElum1, yElum1, rhoElum1;
if( detGeo.elumPos1 != 0 ) { if( aux.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( detGeo.elumPos2 != 0 ) { if( aux.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");
@ -284,13 +285,13 @@ void Transfer(
///in case need other recoil detector. ///in case need other recoil detector.
double xRecoil1, yRecoil1, rhoRecoil1; double xRecoil1, yRecoil1, rhoRecoil1;
if( detGeo.recoilPos1 != 0 ){ if( aux.detPos1 != 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( detGeo.recoilPos2 != 0 ){ if( aux.detPos2 != 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");
@ -544,15 +545,15 @@ void Transfer(
//ELUM //ELUM
if( detGeo.elumPos1 != 0 ){ if( aux.elumPos1 != 0 ){
xElum1 = helios.GetXPos(detGeo.elumPos1); xElum1 = helios.GetXPos(aux.elumPos1);
yElum1 = helios.GetYPos(detGeo.elumPos1); yElum1 = helios.GetYPos(aux.elumPos1);
rhoElum1 = helios.GetR(detGeo.elumPos1); rhoElum1 = helios.GetR(aux.elumPos1);
} }
if( detGeo.elumPos2 != 0 ){ if( aux.elumPos2 != 0 ){
xElum2 = helios.GetXPos(detGeo.elumPos2); xElum2 = helios.GetXPos(aux.elumPos2);
yElum2 = helios.GetYPos(detGeo.elumPos2); yElum2 = helios.GetYPos(aux.elumPos2);
rhoElum2 = helios.GetR(detGeo.elumPos2); rhoElum2 = helios.GetR(aux.elumPos2);
} }
//Recoil //Recoil
@ -563,15 +564,15 @@ void Transfer(
rhoB = orb_B.rho; rhoB = orb_B.rho;
//other recoil detectors //other recoil detectors
if ( detGeo.recoilPos1 != 0 ){ if ( aux.detPos1 != 0 ){
xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1); xRecoil1 = helios.GetRecoilXPos(aux.detPos1);
yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1); yRecoil1 = helios.GetRecoilYPos(aux.detPos1);
rhoRecoil1 = helios.GetRecoilR(detGeo.recoilPos1); rhoRecoil1 = helios.GetRecoilR(aux.detPos1);
} }
if ( detGeo.recoilPos2 != 0 ){ if ( aux.detPos2 != 0 ){
xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2); xRecoil2 = helios.GetRecoilXPos(aux.detPos2);
yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2); yRecoil2 = helios.GetRecoilYPos(aux.detPos2);
rhoRecoil2 = helios.GetRecoilR(detGeo.recoilPos2); rhoRecoil2 = helios.GetRecoilR(aux.detPos2);
} }
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());