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

modified simulation to adopt 2nd array

Browse Source
This commit is contained in:
Ryan Tang 2023-04-10 15:46:24 -04:00
parent 72ce53279d
commit effdbe7617
11 changed files with 348 additions and 445 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
Cleopatra/Check_Simulation.C
View File

@ -47,7 +47,7 @@ enum plotID { pEZ, /// 0
pElum1RThetaCM, /// 16 pElum1RThetaCM, /// 16
pEmpty }; /// 17 pEmpty }; /// 17
plotID StringToPlotID(TString str); plotID StringToPlotID(TString str);
void Check_Simulation(TString filename = "transfer.root", void Check_Simulation(TString filename = "transfer1.root",
TString configFile = "../working/Check_Simulation_Config.txt", TString configFile = "../working/Check_Simulation_Config.txt",
Int_t padSize = 500, Int_t padSize = 500,
bool outputCanvas = false){ bool outputCanvas = false){
@ -139,25 +139,33 @@ void Check_Simulation(TString filename = "transfer.root",
AnalysisLib::DetGeo detGeo = AnalysisLib::LoadDetectorGeo(detGeoTxt); AnalysisLib::DetGeo detGeo = AnalysisLib::LoadDetectorGeo(detGeoTxt);
AnalysisLib::Array array;
if( detGeo.use2ndArray){
array = detGeo.array2;
}else{
array = detGeo.array1;
}
double field = detGeo.Bfield; double field = detGeo.Bfield;
TString fdmsg = field > 0 ? "out of plan" : "into plan"; TString fdmsg = field > 0 ? "out of plan" : "into plan";
TString msg2; TString msg2;
msg2.Form("field = %.2f T, %s", field, fdmsg.Data()); msg2.Form("field = %.2f T, %s", field, fdmsg.Data());
double prepDist = detGeo.array1.detPerpDist; double prepDist = array.detPerpDist;
double length = detGeo.array1.detLength; double length = array.detLength;
double posRecoil = detGeo.recoilPos;
double posRecoil = detGeo.recoilPos;
double rhoRecoilIn = detGeo.recoilInnerRadius; double rhoRecoilIn = detGeo.recoilInnerRadius;
double rhoRecoilOut = detGeo.recoilOuterRadius; double rhoRecoilOut = detGeo.recoilOuterRadius;
double posRecoil1 = detGeo.recoilPos1; double posRecoil1 = detGeo.recoilPos1;
double posRecoil2 = detGeo.recoilPos2; double posRecoil2 = detGeo.recoilPos2;
vector<double> pos = detGeo.array1.detPos; vector<double> pos = array.detPos;
float firstPos = detGeo.array1.firstPos; float firstPos = array.firstPos;
int rDet = detGeo.array1.nDet; int rDet = array.nDet;
int cDet = detGeo.array1.mDet; int cDet = array.mDet;
double elum1 = detGeo.elumPos1; double elum1 = detGeo.elumPos1;
@ -176,8 +184,8 @@ void Check_Simulation(TString filename = "transfer.root",
int numDet = rDet * cDet; int numDet = rDet * cDet;
zRange[1] = detGeo.zMin - 50; zRange[1] = array.zMin - 50;
zRange[2] = detGeo.zMax + 50; zRange[2] = array.zMax + 50;
printf(" zRange : %f - %f \n", zRange[1], zRange[2]); printf(" zRange : %f - %f \n", zRange[1], zRange[2]);
printf("=================================\n"); printf("=================================\n");

425
Cleopatra/HELIOS_LIB.h
View File

@ -444,22 +444,22 @@ public:
bool SetDetectorGeometry(string filename); bool SetDetectorGeometry(string filename);
void SetBeamPosition(double x, double y) { xOff = x; yOff = y;} void SetBeamPosition(double x, double y) { xOff = x; yOff = y;}
void OverrideMagneticField(double BField){ this->Bfield = BField; this->sign = BField > 0 ? 1: -1;} void OverrideMagneticField(double BField){ this->detGeo.Bfield = BField; this->detGeo.BfieldSign = BField > 0 ? 1: -1;}
void OverrideMagneticFieldDirection(double BfieldThetaInDeg){ this->BfieldTheta = BfieldThetaInDeg;} void OverrideMagneticFieldDirection(double BfieldThetaInDeg){ this->detGeo.BfieldTheta = BfieldThetaInDeg;}
void OverrideFirstPos(double firstPos){ void 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);
this->firstPos = firstPos; this->array.firstPos = firstPos;
} }
void OverrideDetectorDistance(double perpDist){ void OverrideDetectorDistance(double perpDist){
overrideDetDistance = true; overrideDetDistance = true;
printf("------ Overriding Detector Distance to : %8.2f mm \n", perpDist); printf("------ Overriding Detector Distance to : %8.2f mm \n", perpDist);
this->perpDist = perpDist; this->array.detPerpDist = perpDist;
} }
void SetDetectorOutside(bool isOutside){ void SetDetectorOutside(bool isOutside){
this->isFromOutSide = isOutside; this->array.detFaceOut = isOutside;
printf(" Detectors are facing %s\n", isFromOutSide ? "outside": "inside" ); printf(" Detectors are facing %s\n", array.detFaceOut ? "outside": "inside" );
} }
int DetAcceptance(); int DetAcceptance();
@ -469,7 +469,7 @@ public:
void CalTrajectoryPara(TLorentzVector P, int Z, int id); // id = 0 for Pb, id = 1 for PB. void CalTrajectoryPara(TLorentzVector P, int Z, int id); // id = 0 for Pb, id = 1 for PB.
int GetNumberOfDetectorsInSamePos(){return mDet;} int GetNumberOfDetectorsInSamePos(){return array.mDet;}
double GetEnergy(){return e;} double GetEnergy(){return e;}
double GetDetX(){return detX;} // position in each detector, range from -1, 1 double GetDetX(){return detX;} // position in each detector, range from -1, 1
@ -486,17 +486,17 @@ public:
return sqrt(x*x+y*y); return sqrt(x*x+y*y);
} }
double GetXPos(double ZPos){ return XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, sign); } double GetXPos(double ZPos){ return XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); }
double GetYPos(double ZPos){ return YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, sign); } double GetYPos(double ZPos){ return YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); }
double GetR(double ZPos) { return RPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, sign); } double GetR(double ZPos) { return RPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); }
double GetRecoilEnergy(){return eB;} double GetRecoilEnergy(){return eB;}
double GetRecoilXPos(double ZPos){ return XPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, sign); } double GetRecoilXPos(double ZPos){ return XPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); }
double GetRecoilYPos(double ZPos){ return YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, sign); } double GetRecoilYPos(double ZPos){ return YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); }
double GetRecoilR(double ZPos) { return RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, sign); } double GetRecoilR(double ZPos) { return RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); }
double GetBField() {return Bfield;} double GetBField() {return detGeo.Bfield;}
double GetDetRadius() {return perpDist;} double GetDetRadius() {return array.detPerpDist;}
trajectory GetTrajectory_b() {return orbitb;} trajectory GetTrajectory_b() {return orbitb;}
trajectory GetTrajectory_B() {return orbitB;} trajectory GetTrajectory_B() {return orbitB;}
@ -523,6 +523,7 @@ private:
} }
AnalysisLib::DetGeo detGeo; AnalysisLib::DetGeo detGeo;
AnalysisLib::Array array;
trajectory orbitb, orbitB; trajectory orbitb, orbitB;
@ -536,22 +537,21 @@ private:
double xOff, yOff; // beam position double xOff, yOff; // beam position
//detetcor Geometry //detetcor Geometry
double Bfield; // T //int sign ; // sign of B-field
int sign ; // sign of B-field // double Bfield; // T
double BfieldTheta; // rad, 0 = z-axis, pi/2 = y axis, pi = -z axis // double BfieldTheta; // rad, 0 = z-axis, pi/2 = y axis, pi = -z axis
double bore; // bore , mm // double bore; // bore , mm
double perpDist; // distance from axis // double perpDist; // distance from axis
double width; // width // double width; // width
double posRecoil; // recoil, downstream // double posRecoil; // recoil, downstream
double rhoRecoilin; // radius recoil inner // double rhoRecoilin; // radius recoil inner
double rhoRecoilout; // radius recoil outter // double rhoRecoilout; // radius recoil outter
double length; // length // double length; // length
double blocker; // double blocker;
double firstPos; // m // double firstPos; // m
vector<double> pos; // near position in m // vector<double> pos; // near position in m
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
// bool isFromOutSide;
bool isFromOutSide;
bool overrideDetDistance; bool overrideDetDistance;
bool overrideFirstPos; bool overrideFirstPos;
@ -575,24 +575,6 @@ HELIOS::HELIOS(){
isDetReady = false; isDetReady = false;
Bfield = 0;
BfieldTheta = 0;
sign = 1;
bore = 1000;
perpDist = 10;
width = 10;
posRecoil = 0;
rhoRecoilin = 0;
rhoRecoilout = 0;
length = 0;
blocker = 0;
firstPos = 0;
pos.clear();
nDet = 0;
mDet = 0;
isFromOutSide = true; //default is facing outside
overrideDetDistance = false; overrideDetDistance = false;
overrideFirstPos = false; overrideFirstPos = false;
isCoincidentWithRecoil = false; isCoincidentWithRecoil = false;
@ -609,25 +591,13 @@ bool HELIOS::SetDetectorGeometry(string filename){
if( haha->ReadFile(filename.c_str()) > 0 ) { if( haha->ReadFile(filename.c_str()) > 0 ) {
detGeo = AnalysisLib::LoadDetectorGeo(haha); detGeo = AnalysisLib::LoadDetectorGeo(haha);
PrintDetGeo(detGeo); if( detGeo.use2ndArray ){
array = detGeo.array2;
}else{
array = detGeo.array1;
}
Bfield = detGeo.Bfield; PrintDetGeo(detGeo, false);
BfieldTheta = detGeo.BfieldTheta;
sign = detGeo.BfieldSign;
bore = detGeo.bore;
posRecoil = detGeo.recoilPos;
rhoRecoilin = detGeo.recoilInnerRadius;
rhoRecoilout = detGeo.recoilOuterRadius;
length = detGeo.array1.detLength;
width = detGeo.array1.detWidth;
perpDist = detGeo.array1.detPerpDist;
blocker = detGeo.array1.blocker;
firstPos = detGeo.array1.firstPos;
pos = detGeo.array1.detPos;
nDet = detGeo.array1.nDet;
mDet = detGeo.array1.mDet;
isFromOutSide = detGeo.array1.detFaceOut;
isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil; isCoincidentWithRecoil = detGeo.isCoincidentWithRecoil;
@ -647,45 +617,45 @@ int HELIOS::DetAcceptance(){
if( isDetReady == false ) return 0; if( isDetReady == false ) return 0;
// -1 ========= when recoil direction is not same side of array // -1 ========= when recoil direction is not same side of array
if( firstPos < 0 && orbitb.z > 0 ) return -1; if( array.firstPos < 0 && orbitb.z > 0 ) return -1;
if( firstPos > 0 && orbitb.z < 0 ) return -1; if( array.firstPos > 0 && orbitb.z < 0 ) return -1;
// -11 ======== rho is too small // -11 ======== rho is too small
if( 2 * orbitb.rho < perpDist ) return -11; if( 2 * orbitb.rho < array.detPerpDist ) return -11;
// -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist) // -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist)
if( orbitb.detRowID == -1 ) return -15; if( orbitb.detRowID == -1 ) return -15;
// -10 =========== when rho is too big . // -10 =========== when rho is too big .
if( bore < 2 * orbitb.rho) return -10; if( detGeo.bore < 2 * orbitb.rho) return -10;
// -14 ========== check particle-B hit radius on recoil dectector // -14 ========== check particle-B hit radius on recoil dectector
if( isCoincidentWithRecoil && orbitB.R > rhoRecoilout ) return -14; if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) return -14;
//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(posRecoil); rhoHit = GetR(detGeo.recoilPos);
if( orbitb.z > 0 && posRecoil > 0 && orbitb.z > posRecoil && rhoHit < rhoRecoilout ) return -12; if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) return -12;
if( orbitb.z < 0 && posRecoil < 0 && orbitb.z < posRecoil && rhoHit < rhoRecoilout ) return -12; if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) return -12;
// -13 ========= not more than 3 loops // -13 ========= not more than 3 loops
if( orbitb.loop > 3 ) return -13; if( orbitb.loop > 3 ) return -13;
// -2 ========= calculate the "y"-distance from detector center // -2 ========= calculate the "y"-distance from detector center
if( sqrt(orbitb.R*orbitb.R - perpDist*perpDist)> width/2 ) return -2; if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) return -2;
// -3 ==== when zPos further the range of whole array, more loop would not save // -3 ==== when zPos further the range of whole array, more loop would not save
if( firstPos < 0 && orbitb.z < pos[0] - length ) return -3; if( array.firstPos < 0 && orbitb.z < array.detPos[0] - array.detLength ) return -3;
if( firstPos > 0 && orbitb.z > pos[nDet-1] + length ) return -3; if( array.firstPos > 0 && orbitb.z > array.detPos[array.nDet-1] + array.detLength ) return -3;
// -4 ======== Hit on blacker // -4 ======== Hit on blacker
if( blocker != 0 && firstPos > 0 && pos[0] - blocker < orbitb.z && orbitb.z < pos[0] ) return -4; if( array.blocker != 0 && array.firstPos > 0 && array.detPos[0] - array.blocker < orbitb.z && orbitb.z < array.detPos[0] ) return -4;
if( blocker != 0 && firstPos < 0 && pos[nDet-1] < orbitb.z && orbitb.z < pos[nDet-1] + blocker ) return -4; if( array.blocker != 0 && array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z && orbitb.z < array.detPos[array.nDet-1] + array.blocker ) return -4;
// 2 ====== when zPos less then the nearest position, more loop may hit // 2 ====== when zPos less then the nearest position, more loop may hit
int increaseLoopFlag = 0; int increaseLoopFlag = 0;
if( firstPos < 0 && pos[nDet-1] < orbitb.z ) increaseLoopFlag = 2; if( array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z ) increaseLoopFlag = 2;
if( firstPos > 0 && pos[0] > orbitb.z ) increaseLoopFlag = 2; if( array.firstPos > 0 && array.detPos[0] > orbitb.z ) increaseLoopFlag = 2;
if (increaseLoopFlag == 2 ) { if (increaseLoopFlag == 2 ) {
orbitb.z += orbitb.z0; orbitb.z += orbitb.z0;
orbitb.effLoop += 1.0; orbitb.effLoop += 1.0;
@ -695,33 +665,34 @@ int HELIOS::DetAcceptance(){
} }
// 1 ======= check hit array z- position // 1 ======= check hit array z- position
if( firstPos < 0 ){ if( array.firstPos < 0 ){
for( int i = 0; i < nDet; i++){ for( int i = 0; i < array.nDet; i++){
if( pos[i]-length <= orbitb.z && orbitb.z <= pos[i]) { if( array.detPos[i] - array.detLength <= orbitb.z && orbitb.z <= array.detPos[i]) {
orbitb.detID = i; orbitb.detID = i;
detX = ( orbitb.z - (pos[i] + length/2 ))/ length*2 ;// range from -1 , 1 detX = ( orbitb.z - (array.detPos[i] + array.detLength/2 ))/ array.detLength * 2 ;// range from -1 , 1
return 1; return 1;
} }
} }
}else{ }else{
for( int i = 0; i < nDet ; i++){ for( int i = 0; i < array.nDet ; i++){
if( pos[i] <= orbitb.z && orbitb.z <= pos[i] + length) { if( array.detPos[i] <= orbitb.z && orbitb.z <= array.detPos[i] + array.detLength) {
///printf(" %d | %f < z = %f < %f \n", i, pos[i], orbitb.z, pos[i]+length); ///printf(" %d | %f < z = %f < %f \n", i, array.detPos[i], orbitb.z, array.detPos[i]+length);
orbitb.detID = i; orbitb.detID = i;
detX = ( orbitb.z - (pos[i] - length/2 ))/ length*2 ;// range from -1 , 1 detX = ( orbitb.z - (array.detPos[i] - array.detLength/2 ))/ array.detLength*2 ;// range from -1 , 1
return 1; return 1;
} }
} }
} }
// -5 ======== check hit array gap // -5 ======== check hit array gap
if( firstPos < 0 ){ if( array.firstPos < 0 ){
for( int i = 0; i < nDet-1 ; i++){ for( int i = 0; i < array.nDet-1 ; i++){
if( pos[i] < orbitb.z && orbitb.z < pos[i+1] - length ) return -5; //increaseLoopFlag = 3; if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) return -5; //increaseLoopFlag = 3;
} }
}else{ }else{
for( int i = 0; i < nDet-1 ; i++){ for( int i = 0; i < array.nDet-1 ; i++){
if( pos[i] + length < orbitb.z && orbitb.z < pos[i+1] ) return -5; //increaseLoopFlag = 3; if( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) return -5; //increaseLoopFlag = 3;
} }
} }
if (increaseLoopFlag == 3 ) { if (increaseLoopFlag == 3 ) {
@ -741,7 +712,7 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){
if( id == 0 ){ if( id == 0 ){
orbitb.theta = P.Theta(); orbitb.theta = P.Theta();
orbitb.phi = P.Phi(); orbitb.phi = P.Phi();
orbitb.rho = P.Pt() / abs(Bfield) / Z / c * 1000; //mm orbitb.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm
orbitb.vt = P.Beta() * TMath::Sin(P.Theta()) * c ; // mm / nano-second orbitb.vt = P.Beta() * TMath::Sin(P.Theta()) * c ; // mm / nano-second
orbitb.vp = P.Beta() * TMath::Cos(P.Theta()) * c ; // mm / nano-second orbitb.vp = P.Beta() * TMath::Cos(P.Theta()) * c ; // mm / nano-second
orbitb.t0 = TMath::TwoPi() * orbitb.rho / orbitb.vt; // nano-second orbitb.t0 = TMath::TwoPi() * orbitb.rho / orbitb.vt; // nano-second
@ -755,7 +726,7 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){
if( id == 1 ){ if( id == 1 ){
orbitB.theta = P.Theta(); orbitB.theta = P.Theta();
orbitB.phi = P.Phi(); orbitB.phi = P.Phi();
orbitB.rho = P.Pt() / abs(Bfield) / Z / c * 1000; //mm orbitB.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm
orbitB.vt = P.Beta() * TMath::Sin(P.Theta()) * c ; // mm / nano-second orbitB.vt = P.Beta() * TMath::Sin(P.Theta()) * c ; // mm / nano-second
orbitB.vp = P.Beta() * TMath::Cos(P.Theta()) * c ; // mm / nano-second orbitB.vp = P.Beta() * TMath::Cos(P.Theta()) * c ; // mm / nano-second
orbitB.t0 = TMath::TwoPi() * orbitB.rho / orbitB.vt; // nano-second orbitB.t0 = TMath::TwoPi() * orbitB.rho / orbitB.vt; // nano-second
@ -776,7 +747,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
CalTrajectoryPara(Pb, Zb, 0); CalTrajectoryPara(Pb, Zb, 0);
int targetLoop = 1; int targetLoop = 1;
int inOut = isFromOutSide == true ? 1: 0; //1 = from Outside, 0 = from inside int inOut = array.detFaceOut == true ? 1: 0; //1 = from Outside, 0 = from inside
bool debug = false; bool debug = false;
@ -792,18 +763,18 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
vector<double> zPossible; vector<double> zPossible;
vector<int> dID; //detRowID vector<int> dID; //detRowID
int iStart = ( sign == 1 ? 0 : -mDet ); int iStart = ( detGeo.BfieldSign == 1 ? 0 : -array.mDet );
int iEnd = ( sign == 1 ? 2*mDet : mDet ); int iEnd = ( detGeo.BfieldSign == 1 ? 2 * array.mDet : array.mDet );
for( int i = iStart; i < iEnd ; i++){ for( int i = iStart; i < iEnd ; i++){
double phiD = TMath::TwoPi()/mDet * i ; double phiD = TMath::TwoPi()/array.mDet * i ;
double dphi = orbitb.phi - phiD; double dphi = orbitb.phi - phiD;
double aEff = perpDist - (xOff * TMath::Cos(phiD) + yOff * TMath::Sin(phiD)); double aEff = array.detPerpDist - (xOff * TMath::Cos(phiD) + yOff * TMath::Sin(phiD));
double hahaha = asin( aEff/ orbitb.rho - sign * sin(dphi)); double hahaha = asin( aEff/ orbitb.rho - detGeo.BfieldTheta * sin(dphi));
int n = 2*targetLoop + inOut; int n = 2*targetLoop + inOut;
double zP = orbitb.z0 /TMath::TwoPi() * ( sign * dphi + n * TMath::Pi() + pow(-1, n) * hahaha ); double zP = orbitb.z0 /TMath::TwoPi() * ( detGeo.BfieldSign * dphi + n * TMath::Pi() + pow(-1, n) * hahaha );
if( debug ) { if( debug ) {
double xP = GetXPos(zP) ; double xP = GetXPos(zP) ;
@ -848,15 +819,15 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
orbitb.detRowID = (12+dID[i])%4; orbitb.detRowID = (12+dID[i])%4;
orbitb.t = orbitb.t0 * orbitb.effLoop; orbitb.t = orbitb.t0 * orbitb.effLoop;
double phiD = TMath::TwoPi()/mDet * dID[i] ; double phiD = TMath::TwoPi()/array.mDet * dID[i] ;
double dphi = orbitb.phi - phiD ; double dphi = orbitb.phi - phiD ;
if( debug ) { if( debug ) {
// Check is in or out // Check is in or out
double hitDir = cos( orbitb.z/orbitb.z0 * TMath::TwoPi() - sign * dphi ); double hitDir = cos( orbitb.z/orbitb.z0 * TMath::TwoPi() - detGeo.BfieldSign * dphi );
printf(" hitDir : %4.1f ", hitDir); printf(" hitDir : %4.1f ", hitDir);
if( ( inOut == 1 && hitDir > 0 ) || (inOut == 0 && hitDir < 0 ) ) { if( ( inOut == 1 && hitDir > 0 ) || (inOut == 0 && hitDir < 0 ) ) {
printf(" != %f ", perpDist); printf(" != %f ", array.detPerpDist);
orbitb.z = TMath::QuietNaN(); orbitb.z = TMath::QuietNaN();
orbitb.loop = -1; orbitb.loop = -1;
orbitb.detRowID = -1; orbitb.detRowID = -1;
@ -868,8 +839,8 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
double yPos = GetYPos(orbitb.z ) ; double yPos = GetYPos(orbitb.z ) ;
double a = xPos * cos(phiD) + yPos * sin(phiD); double a = xPos * cos(phiD) + yPos * sin(phiD);
printf(" a : %f ", a); printf(" a : %f ", a);
if( abs(a - perpDist) > 0.01) { if( abs(a - array.detPerpDist) > 0.01) {
printf(" != %f ", perpDist); printf(" != %f ", array.detPerpDist);
orbitb.z = TMath::QuietNaN(); orbitb.z = TMath::QuietNaN();
orbitb.loop = -1; orbitb.loop = -1;
orbitb.detRowID = -1; orbitb.detRowID = -1;
@ -892,244 +863,16 @@ int HELIOS::CalRecoilHit(TLorentzVector PB, int ZB){
CalTrajectoryPara(PB, ZB, 1); CalTrajectoryPara(PB, ZB, 1);
orbitB.z = posRecoil; orbitB.z = detGeo.recoilPos;
orbitB.x = GetRecoilXPos(posRecoil) ; orbitB.x = GetRecoilXPos(detGeo.recoilPos) ;
orbitB.y = GetRecoilYPos(posRecoil) ; orbitB.y = GetRecoilYPos(detGeo.recoilPos) ;
orbitB.R = GetRecoilR(posRecoil); 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 ;
return 1; return 1;
} }
/*
int HELIOS::CalHit(TLorentzVector Pb, int Zb, TLorentzVector PB, int ZB, double xOff, double yOff){
//initialization
int hit = 0;
const double c = 299.792458; // mm/ns, standard speed for MeV conversion
theta = TMath::QuietNaN();
e = TMath::QuietNaN();
z = TMath::QuietNaN();
detX = TMath::QuietNaN();
t = TMath::QuietNaN();
rho = TMath::QuietNaN();
dphi = TMath::QuietNaN();
detID = -1;
loop = -1;
detRowID = -1;
eB = TMath::QuietNaN();
zB = TMath::QuietNaN();
tB = TMath::QuietNaN();
rhoB = TMath::QuietNaN();
rxHit = TMath::QuietNaN();
ryHit = TMath::QuietNaN();
phiB = TMath::QuietNaN();
//range of detector azimuth angle, for beam at center
double azimu = TMath::Pi()/ mDet;
double azimuDet = TMath::ATan2(width/2., perpDist);
//rotate Pb and PB to B-Field
//Pb.RotateX(BfieldTheta);
//PB.RotateX(BfieldTheta);
//====================== X-Y plane, light particle
rho = Pb.Pt() / abs(Bfield) / Zb / c * 1000; //mm
theta = Pb.Theta();
phi = Pb.Phi();
//====================== recoil detector
thetaB = PB.Theta();
phiB = PB.Phi();
rhoB = PB.Pt() / abs(Bfield) / ZB / c * 1000; //mm
vt0B = PB.Beta() * TMath::Sin(thetaB) * c ; // mm / nano-second
vp0B = PB.Beta() * TMath::Cos(thetaB) * c ; // mm / nano-second
//tB = TMath::TwoPi() * rhoB / vt0B; // nano-second
tB = posRecoil / vp0B; // nano-second
eB = PB.E() - PB.M();
zB = vp0B * tB;
rhoBHit = GetRecoilR(posRecoil);
rxHit = GetRecoilXPos(posRecoil);
ryHit = GetRecoilYPos(posRecoil);
if( isDetReady == false ) {
//====================== infinite small detector
vt0 = Pb.Beta() * TMath::Sin(theta) * c ; // mm / nano-second
vp0 = Pb.Beta() * TMath::Cos(theta) * c ; // mm / nano-second
t = TMath::TwoPi() * rho / vt0; // nano-second
z = vp0 * t; // mm
e = Pb.E() - Pb.M();
dphi = TMath::TwoPi();
return 0;
}
if( bore > 2 * rho && 2*rho > perpDist && ((firstPos > 0 && theta < TMath::PiOver2()) || (firstPos < 0 && theta > TMath::PiOver2())) ){
//====================== infinite small detector
vt0 = Pb.Beta() * TMath::Sin(theta) * c ; // mm / nano-second
vp0 = Pb.Beta() * TMath::Cos(theta) * c ; // mm / nano-second
t0 = TMath::TwoPi() * rho / vt0; // nano-second
z0 = vp0 * t0; // mm
//========== check particle-B hit radius on recoil dectector
if(isCoincidentWithRecoil && rhoBHit > rhoRecoilout ) return -2; // when particle-B miss the recoil detector
//========= check is particle-b was blocked by recoil detector
rhoHit = GetR(posRecoil) ;// radius of light particle b at recoil detector
if( z0 > 0 && posRecoil > 0 && z0 > posRecoil && rhoHit < rhoRecoilout) {
return -1 ; // when particle-b blocked by recoil detector
}
if( z0 < 0 && posRecoil < 0 && z0 < posRecoil && rhoHit < rhoRecoilout) {
return -1 ;
}
//================ Calulate z hit
double zHit = TMath::QuietNaN();
bool isHit = false;
bool isHitFromOutside = false;
bool isReachArrayCoverage = false;
loop = 0;
int startJ = (int) fmod(TMath::Ceil(mDet*phi/TMath::TwoPi() - 0.5) ,mDet) ;
//printf("======================= T : %5.2f, theta : %7.2f , phi : %7.2f \n", Pb.E() - Pb.M(), theta*TMath::RadToDeg(), phi*TMath::RadToDeg());
// loop until reach the detector position covrage.
do{
loop += 1;
//int n = 2*loop + sign;
int n = 2*loop -1;
if( blocker != 0.0 && abs(firstPos/blocker) < loop ) return -6;
//======= check Block By blocker using z0, speed thing up
if( firstPos > 0 ){
if( pos[0] - blocker < z0 * loop && z0 * loop < pos[0] ) return -6; // blocked by blocker
}else{
if( pos[nDet-1] < z0 * loop && z0 * loop < pos[nDet-1] + blocker ) return -6;
}
if( loop > 10 ) {
return -3; // when loop > 10
break; // maximum 10 loops
}
for( int j = startJ ; j < startJ + mDet; j++){
double phiDet = TMath::TwoPi() / mDet * (j); // detector plane angle
isReachArrayCoverage = false;
isHitFromOutside = false;
//========== calculate zHit by solving the cycltron orbit and the detector planes.
double aEff = perpDist - (xOff * TMath::Cos(phiDet) + yOff * TMath::Sin(phiDet));
double dphi = phi - phiDet;
double z0 = TMath::TwoPi() * rho / TMath::Tan(theta); // the cycle
//with sign is the correct formula, but somehow, probably the ASin? give incorrect result for sign = 1;
//zHit = z0 / TMath::TwoPi() * ( sign * dphi + TMath::Power(-1, n) * TMath::ASin(aEff/rho - sign * TMath::Sin(dphi)) + TMath::Pi() * n );
zHit = z0 / TMath::TwoPi() * ( -1 * dphi + TMath::Power(-1, n) * TMath::ASin(aEff/rho + TMath::Sin(dphi)) + TMath::Pi() * n );
e = Pb.E() - Pb.M();
z = zHit;
//======= calculate the distance from middle of detector
double xHit = GetXPos(zHit) + xOff; // resotre the beam to be at the center
double yHit = GetYPos(zHit) + yOff;
double sHit = TMath::Sqrt(xHit*xHit + yHit*yHit - perpDist*perpDist); // is it hit on the detector
///if( zHit > z0 ) continue;
///printf("==== %7.2f, %7.2f | n : %d, row : %2d, phiD : %4.0f, rho : %9.4f, z0 : %9.4f, zHit : %9.4f, xHit : %9.4f, yHit : %9.4f \n",
/// theta*TMath::RadToDeg(), phi*TMath::RadToDeg(), n, j, phiDet*TMath::RadToDeg(), rho, z0, zHit, xHit, yHit);
if( sHit > width /2.) continue; // if the sHit is large, it does not hit on detector, go to next mDet
//======== this is the particel direction (normalized) dot normal vector of the detector plane
//double dir = TMath::Cos(zHit/z0 * TMath::TwoPi() - sign * dphi);
double dir = TMath::Cos(zHit/z0 * TMath::TwoPi() + dphi);
if( dir < 0) {// when dir == 0, no solution
isHitFromOutside = true;
}else{
return -5 ; // hit from inside.
}
//### after this, the particle is hit on detector, from outside.
//======= check Block By blocker using z0 and zHit; z0 has to be inside the nearst array, but zHit is in the blocker.
if( firstPos > 0 && blocker > 0.0){
if( pos[0] < z0*loop && pos[0] - blocker < zHit && zHit < pos[0] ) return -6;
}else{
if( z0*loop < pos[nDet-1] && pos[nDet-1] < zHit && zHit < pos[nDet-1] + blocker ) return -6;
}
//### after this, the particle is hit on detector, from outside, and not blocked by blocker
//======= check within the detector range
if( firstPos > 0 ){
if( zHit < pos[0] ) continue; // goto next mDet, after progress of all side, still not OK, then next loop
if( zHit > pos[nDet-1] + length) return -4; // since the zHit is mono-increse, when zHit shoot over the detector
}else{
if( zHit < pos[0] - length ) return -4;
if( zHit > pos[nDet-1]) continue;
}
//====== check hit
if( !isReachArrayCoverage && isHitFromOutside && sHit < width/2.){
isHit = true;
isReachArrayCoverage = true;
detRowID = (j+mDet) % mDet;
break; // if isHit, break, don't calculate for the rest of the detector
}else{
isReachArrayCoverage = false;
}
}
}while(!isReachArrayCoverage);
if( !isHit ) return -7; // zHit in the gap of detector
//===== final calculation for light particle
e = Pb.E() - Pb.M();
z = zHit;
t = zHit / vp0;
dphi = t * vt0 / rho;
//sometimes, dphi = 2pi + something, i.e. loop a bit more than a single loop.
if( dphi / TMath::TwoPi() > loop ) return -8;
double xHit = GetXPos(zHit) + xOff; // resotre the beam to be at the center
double yHit = GetYPos(zHit) + yOff;
//======= Check inside the detector
double eta = TMath::Abs(TMath::ATan2(yHit,xHit));
double etaMod = TMath::Abs(fmod(eta + azimu, 2* azimu) - azimu);
if ( etaMod > azimuDet ) return -8;
//=========== check hit on detector gap
for( int i = 0 ; i < nDet ; i++){
if( firstPos > 0 ){
if( pos[i] < z && z < pos[i] + length ){
detID = i;
detX = ( z - (pos[i] + length/2 ))/ length*2 ;// range from -1 , 1
}
}else{
if( pos[i] - length < z && z < pos[i] ){
detID = i;
detX = ( z - (pos[i] - length/2 ))/ length*2 ;// range from -1 , 1
}
}
}
if( detID >=0 ){
hit = 1;
}else{
hit = -9; // particle-b hit the gap
}
}else{
hit = -10; // hit helio wall or (detector upstream, particle go downstream, via versa)
}
return hit;
}
*/
//======================================================= //=======================================================
//####################################################### //#######################################################
// Class for multi-scattering of the beam inside target // Class for multi-scattering of the beam inside target

107
Cleopatra/Simulation_Helper.C
View File

@ -17,6 +17,7 @@
#include "../Cleopatra/ExtractXSec.h" #include "../Cleopatra/ExtractXSec.h"
#include "../Cleopatra/PlotTGraphTObjArray.h" #include "../Cleopatra/PlotTGraphTObjArray.h"
#include "../armory/AutoFit.C" #include "../armory/AutoFit.C"
#include "../armory/AnalysisLib.h"
#include "../Cleopatra/Check_Simulation.C" #include "../Cleopatra/Check_Simulation.C"
#include <iostream> #include <iostream>
@ -60,7 +61,8 @@ private:
TGTextEntry * txtName ; TGTextEntry * txtName ;
TGTextEntry * txtEx ; TGTextEntry * txtEx ;
bool isUse2ndArray;
public: public:
MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h); MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h);
@ -69,6 +71,8 @@ public:
void OpenFile(int); void OpenFile(int);
void GetData(); void GetData();
bool IsFileExist(TString filename); bool IsFileExist(TString filename);
void CheckIsUse2ndArray();
}; };
@ -85,7 +89,7 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) {
TGVerticalFrame *hframe2 = new TGVerticalFrame(fMain,600,800 ); TGVerticalFrame *hframe2 = new TGVerticalFrame(fMain,600,800 );
hframe->AddFrame(hframe2,new TGLayoutHints( kLHintsExpandX | kLHintsExpandY, 2,2,2,2)); hframe->AddFrame(hframe2,new TGLayoutHints( kLHintsExpandX | kLHintsExpandY, 2,2,2,2));
fileName = "../working/reactionConfig.txt"; fileName = "../working/detectorGeo.txt";
TGHorizontalFrame *hframe00 = new TGHorizontalFrame(hframe2,600,600 ); TGHorizontalFrame *hframe00 = new TGHorizontalFrame(hframe2,600,600 );
hframe2->AddFrame(hframe00, new TGLayoutHints(kLHintsCenterX | kLHintsExpandX , 2,2,2,2)); hframe2->AddFrame(hframe00, new TGLayoutHints(kLHintsCenterX | kLHintsExpandX , 2,2,2,2));
@ -127,13 +131,6 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) {
{//================= Simulation group {//================= Simulation group
TGGroupFrame * simFrame = new TGGroupFrame(hframe1, "Kinematics Simulation", kVerticalFrame); TGGroupFrame * simFrame = new TGGroupFrame(hframe1, "Kinematics Simulation", kVerticalFrame);
hframe1->AddFrame(simFrame, new TGLayoutHints(kLHintsCenterX, 5,5,3,4)); hframe1->AddFrame(simFrame, new TGLayoutHints(kLHintsCenterX, 5,5,3,4));
TGTextButton *openRec = new TGTextButton(simFrame, "reaction Config");
openRec->SetWidth(150);
openRec->SetHeight(20);
openRec->ChangeOptions( openRec->GetOptions() | kFixedSize );
openRec->Connect("Clicked()","MyMainFrame",this, "OpenFile(=1)");
simFrame->AddFrame(openRec,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openDet = new TGTextButton(simFrame, "detector Geo."); TGTextButton *openDet = new TGTextButton(simFrame, "detector Geo.");
openDet->SetWidth(150); openDet->SetWidth(150);
@ -142,6 +139,13 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) {
openDet->Connect("Clicked()","MyMainFrame",this, "OpenFile(=0)"); openDet->Connect("Clicked()","MyMainFrame",this, "OpenFile(=0)");
simFrame->AddFrame(openDet,new TGLayoutHints(kLHintsRight, 5,5,3,4)); simFrame->AddFrame(openDet,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openRec = new TGTextButton(simFrame, "reaction Config");
openRec->SetWidth(150);
openRec->SetHeight(20);
openRec->ChangeOptions( openRec->GetOptions() | kFixedSize );
openRec->Connect("Clicked()","MyMainFrame",this, "OpenFile(=1)");
simFrame->AddFrame(openRec,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openEx = new TGTextButton(simFrame, "Ex List"); TGTextButton *openEx = new TGTextButton(simFrame, "Ex List");
openEx->SetWidth(150); openEx->SetWidth(150);
openEx->SetHeight(20); openEx->SetHeight(20);
@ -168,16 +172,16 @@ MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) {
openSimChk->Connect("Clicked()","MyMainFrame",this, "OpenFile(=4)"); openSimChk->Connect("Clicked()","MyMainFrame",this, "OpenFile(=4)");
simFrame->AddFrame(openSimChk,new TGLayoutHints(kLHintsRight, 5,5,3,4)); simFrame->AddFrame(openSimChk,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *SimChk = new TGTextButton(simFrame,"Plot Simulation"); TGTextButton *SimChk = new TGTextButton(simFrame,"Re-Plot Simulation");
SimChk->SetWidth(150); SimChk->SetWidth(150);
SimChk->SetHeight(40); //SimChk->SetHeight(40);
SimChk->ChangeOptions( SimChk->GetOptions() | kFixedSize ); SimChk->ChangeOptions( SimChk->GetOptions() | kFixedSize );
SimChk->Connect("Clicked()","MyMainFrame",this,"Command(=2)"); SimChk->Connect("Clicked()","MyMainFrame",this,"Command(=2)");
simFrame->AddFrame(SimChk, new TGLayoutHints(kLHintsRight,5,5,3,4)); simFrame->AddFrame(SimChk, new TGLayoutHints(kLHintsRight,5,5,3,4));
TGTextButton *autoFit = new TGTextButton(simFrame,"AutoFit ExCal"); TGTextButton *autoFit = new TGTextButton(simFrame,"AutoFit ExCal");
autoFit->SetWidth(150); autoFit->SetWidth(150);
autoFit->SetHeight(40); //autoFit->SetHeight(40);
autoFit->ChangeOptions( autoFit->GetOptions() | kFixedSize ); autoFit->ChangeOptions( autoFit->GetOptions() | kFixedSize );
autoFit->Connect("Clicked()","MyMainFrame",this,"Command(=5)"); autoFit->Connect("Clicked()","MyMainFrame",this,"Command(=5)");
simFrame->AddFrame(autoFit, new TGLayoutHints(kLHintsRight,5,5,3,4)); simFrame->AddFrame(autoFit, new TGLayoutHints(kLHintsRight,5,5,3,4));
@ -352,20 +356,40 @@ bool MyMainFrame::IsFileExist(TString filename){
return file.is_open(); return file.is_open();
} }
void MyMainFrame::CheckIsUse2ndArray(){
TMacro * haha = new TMacro("../working/detectorGeo.txt");
AnalysisLib::DetGeo detGeo = AnalysisLib::LoadDetectorGeo(haha);
delete haha;
isUse2ndArray = detGeo.use2ndArray;
}
void MyMainFrame::OpenFile(int ID){ void MyMainFrame::OpenFile(int ID){
editor->SaveFile(fileName); editor->SaveFile(fileName);
TString oldFileName = fileName; TString oldFileName = fileName;
if ( ID == 0 ) fileName = "../working/detectorGeo.txt"; if ( ID == 0 ) fileName = "../working/detectorGeo.txt";
if ( ID == 1 ) fileName = "../working/reactionConfig.txt";
if ( ID == 2 ) fileName = "../working/Ex.txt"; CheckIsUse2ndArray();
if ( ID == 3 ) fileName = "../working/DWBA"; if( isUse2ndArray){
if ( ID == 1 ) fileName = "../working/reactionConfig2.txt";
if ( ID == 2 ) fileName = "../working/Ex2.txt";
if ( ID == 3 ) fileName = "../working/DWBA2";
if ( ID == 5 ) fileName = "../working/DWBA2.in";
if ( ID == 6 ) fileName = "../working/DWBA2.out";
if ( ID == 7 ) fileName = "../working/DWBA2.Xsec.txt";
}else{
if ( ID == 1 ) fileName = "../working/reactionConfig1.txt";
if ( ID == 2 ) fileName = "../working/Ex1.txt";
if ( ID == 3 ) fileName = "../working/DWBA1";
if ( ID == 5 ) fileName = "../working/DWBA1.in";
if ( ID == 6 ) fileName = "../working/DWBA1.out";
if ( ID == 7 ) fileName = "../working/DWBA1.Xsec.txt";
}
if ( ID == 4 ) fileName = "../working/Check_Simulation_Config.txt"; if ( ID == 4 ) fileName = "../working/Check_Simulation_Config.txt";
if ( ID == 5 ) fileName = "../working/DWBA.in";
if ( ID == 6 ) fileName = "../working/DWBA.out";
if ( ID == 7 ) fileName = "../working/DWBA.Xsec.txt";
if ( ID == 8 ) fileName = isoFileName; if ( ID == 8 ) fileName = isoFileName;
//test if file exist //test if file exist
@ -417,6 +441,8 @@ void MyMainFrame::GetData(){
void MyMainFrame::Command(int ID) { void MyMainFrame::Command(int ID) {
editor->SaveFile(fileName); editor->SaveFile(fileName);
CheckIsUse2ndArray();
if( ID == 0 ){ if( ID == 0 ){
@ -468,29 +494,54 @@ void MyMainFrame::Command(int ID) {
} }
if( ID == 1 ){ if( ID == 1 ){
string basicConfig = "reactionConfig.txt";
string basicConfig = "reactionConfig1.txt";
string heliosDetGeoFile = "detectorGeo.txt"; string heliosDetGeoFile = "detectorGeo.txt";
string excitationFile = "Ex.txt"; //when no file, only ground state string excitationFile = "Ex1.txt"; //when no file, only ground state
TString ptolemyRoot = ""; // when no file, use isotropic distribution of thetaCM TString ptolemyRoot = ""; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root"; TString saveFileName = "transfer1.root";
TString filename = "reaction.dat"; //when no file, no output TString filename = "reaction1.dat"; //when no file, no output
if( withDWBA->GetState() ) { if( withDWBA->GetState() ) {
ptolemyRoot = "DWBA.root"; ptolemyRoot = "DWBA1.root";
excitationFile = "DWBA.Ex.txt"; excitationFile = "DWBA1.Ex.txt";
} }
if( isUse2ndArray ){
basicConfig = "reactionConfig2.txt";
heliosDetGeoFile = "detectorGeo.txt";
excitationFile = "Ex2.txt"; //when no file, only ground state
ptolemyRoot = ""; // when no file, use isotropic distribution of thetaCM
saveFileName = "transfer2.root";
filename = "reaction2.dat"; //when no file, no output
if( withDWBA->GetState() ) {
ptolemyRoot = "DWBA2.root";
excitationFile = "DWBA2.Ex.txt";
}
}
statusLabel->SetText("Running simulation......."); statusLabel->SetText("Running simulation.......");
Transfer( basicConfig, heliosDetGeoFile, excitationFile, ptolemyRoot, saveFileName, filename); Transfer( basicConfig, heliosDetGeoFile, excitationFile, ptolemyRoot, saveFileName, filename);
statusLabel->SetText("Plotting simulation......."); statusLabel->SetText("Plotting simulation.......");
Check_Simulation();
if( isUse2ndArray ){
Check_Simulation("transfer2.root");
}else{
Check_Simulation("transfer1.root");
}
statusLabel->SetText("Plotted Simulation result"); statusLabel->SetText("Plotted Simulation result");
} }
if( ID == 2 ){ if( ID == 2 ){
Check_Simulation(); if( isUse2ndArray ){
statusLabel->SetText(" Run Simulation first."); Check_Simulation("transfer2.root");
}else{
Check_Simulation("transfer1.root");
}
statusLabel->SetText(" Run Simulation first.");
} }
if( ID == 3 ){ if( ID == 3 ){

183
armory/AnalysisLib.h
View File

@ -47,7 +47,7 @@ std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, in
}while(pos != std::string::npos ); }while(pos != std::string::npos );
return output; return output;
} };
struct Array{ struct Array{
@ -117,7 +117,7 @@ struct DetGeo{
Array array1; Array array1;
//==================2nd array //==================2nd array
bool is2ndArrayExist; bool use2ndArray;
Array array2; Array array2;
double zMin, zMax; /// range of detectors double zMin, zMax; /// range of detectors
@ -169,7 +169,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){
detGeo.array1.pos.clear(); detGeo.array1.pos.clear();
detGeo.array2.pos.clear(); detGeo.array2.pos.clear();
detGeo.is2ndArrayExist = false; detGeo.use2ndArray = false;
int detFlag = 0; int detFlag = 0;
int detLine = 0; int detLine = 0;
@ -219,7 +219,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){
} }
if( detFlag == 2){ if( detFlag == 2){
if ( detLine == 0 ) detGeo.is2ndArrayExist = str[0] == "false" ? false: true; if ( detLine == 0 ) detGeo.use2ndArray = str[0] == "true" ? true : false;
if ( detLine == 1 ) detGeo.array2.detPerpDist = atof(str[0].c_str()); if ( detLine == 1 ) detGeo.array2.detPerpDist = atof(str[0].c_str());
if ( detLine == 2 ) detGeo.array2.detWidth = atof(str[0].c_str()); if ( detLine == 2 ) detGeo.array2.detWidth = atof(str[0].c_str());
if ( detLine == 3 ) detGeo.array2.detLength = atof(str[0].c_str()); if ( detLine == 3 ) detGeo.array2.detLength = atof(str[0].c_str());
@ -240,7 +240,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){
detGeo.zMin = detGeo.array1.zMin; detGeo.zMin = detGeo.array1.zMin;
detGeo.zMax = detGeo.array1.zMax; detGeo.zMax = detGeo.array1.zMax;
if( detGeo.is2ndArrayExist) { if( detGeo.use2ndArray) {
detGeo.array2.DeduceAbsolutePos(); detGeo.array2.DeduceAbsolutePos();
detGeo.zMin = TMath::Min(detGeo.array1.zMin, detGeo.array2.zMin); detGeo.zMin = TMath::Min(detGeo.array1.zMin, detGeo.array2.zMin);
@ -250,7 +250,7 @@ DetGeo LoadDetectorGeo(TMacro * macro){
return detGeo; return detGeo;
} }
void PrintDetGeo(DetGeo detGeo){ void PrintDetGeo(DetGeo detGeo, bool printAll = true){
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, Theta : %6.2f deg \n", detGeo.Bfield, detGeo.BfieldTheta);
@ -258,15 +258,24 @@ void PrintDetGeo(DetGeo detGeo){
printf(" +---- field angle != 0 is not supported!!! \n"); 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(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", detGeo.recoilPos, detGeo.recoilInnerRadius, detGeo.recoilOuterRadius);
printf("------------------------------------- Detector Position \n"); if( printAll ){
detGeo.array1.PrintArray(); printf("------------------------------------- Detector Position \n");
detGeo.array1.PrintArray();
if( detGeo.is2ndArrayExist){ if( detGeo.use2ndArray){
printf("--------------------------------- 2nd Detector Position \n"); printf("--------------------------------- 2nd Detector Position \n");
detGeo.array2.PrintArray(); detGeo.array2.PrintArray();
}
}else{
if( detGeo.use2ndArray){
printf("--------------------------------- 2nd Detector Position \n");
detGeo.array2.PrintArray();
}else{
printf("------------------------------------- Detector Position \n");
detGeo.array1.PrintArray();
}
} }
if( detGeo.elumPos1 != 0 || detGeo.elumPos2 != 0 || detGeo.recoilPos1 != 0 || detGeo.recoilPos2 != 0){ if( detGeo.elumPos1 != 0 || detGeo.elumPos2 != 0 || detGeo.recoilPos1 != 0 || detGeo.recoilPos2 != 0){
printf("=================================== Auxillary/Imaginary Detectors\n"); printf("=================================== Auxillary/Imaginary Detectors\n");
} }
@ -334,7 +343,7 @@ ReactionConfig LoadReactionConfig(TMacro * macro){
reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA; reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA;
reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ; reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ;
return reaction; return reaction;
} }
@ -378,9 +387,12 @@ void PrintReactionConfig(ReactionConfig reaction){
} }
DetGeo detGeo; DetGeo detGeo;
ReactionConfig reactionConfig; ReactionConfig reactionConfig1;
ReactionConfig reactionConfig2;
void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt", std::string reactionConfigFileName = "reactionConfig.txt"){ void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt",
std::string reactionConfigFileName1 = "reactionConfig1.txt",
std::string reactionConfigFileName2 = "reactionConfig2.txt"){
printf("=====================================================\n"); printf("=====================================================\n");
printf(" loading detector geometery : %s.", detGeoFileName.c_str()); printf(" loading detector geometery : %s.", detGeoFileName.c_str());
TMacro * haha = new TMacro(); TMacro * haha = new TMacro();
@ -391,21 +403,34 @@ void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.t
}else{ }else{
printf("... fail\n"); printf("... fail\n");
} }
delete haha;
printf("=====================================================\n"); printf("=====================================================\n");
printf(" loading reaction config : %s.", reactionConfigFileName.c_str()); printf(" loading reaction1 config : %s.", reactionConfigFileName1.c_str());
TMacro * kaka = new TMacro(); TMacro * kaka = new TMacro();
if( kaka->ReadFile(reactionConfigFileName.c_str()) > 0 ) { if( kaka->ReadFile(reactionConfigFileName1.c_str()) > 0 ) {
reactionConfig = AnalysisLib::LoadReactionConfig(kaka); reactionConfig1 = AnalysisLib::LoadReactionConfig(kaka);
printf("..... done.\n"); printf("..... done.\n");
AnalysisLib::PrintReactionConfig(reactionConfig); AnalysisLib::PrintReactionConfig(reactionConfig1);
}else{ }else{
printf("..... fail\n"); printf("..... fail\n");
} }
delete haha;
delete kaka; delete kaka;
if( detGeo.use2ndArray){
printf("=====================================================\n");
printf(" loading reaction2 config : %s.", reactionConfigFileName2.c_str());
TMacro * jaja = new TMacro();
if( jaja->ReadFile(reactionConfigFileName2.c_str()) > 0 ) {
reactionConfig2 = AnalysisLib::LoadReactionConfig(kaka);
printf("..... done.\n");
AnalysisLib::PrintReactionConfig(reactionConfig2);
}else{
printf("..... fail\n");
}
delete jaja;
}
} }
//************************************** Correction parameters; //************************************** Correction parameters;
@ -502,11 +527,28 @@ void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.d
if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]); if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]);
} }
double q, alpha, Et, betRel, gamm, G, massB, mass; //variables for Ex calculation
bool hasReactionPara = false; struct ReactionParas{
double Et; // total energy in CM frame
double beta; // Lorentz beta from Lab to CM
double gamma; // Lorentz gamma from Lab to CM
double alpha; // E-Z slope / beta
double G; //The G-coefficient....
double massB; // heavy mass
double q; // charge of light particle
double mass; //light mass
bool hasReactionPara;
double detPrepDist;
};
ReactionParas reactParas1;
ReactionParas reactParas2;
//~========================================= reaction parameters //~========================================= reaction parameters
void LoadReactionParasForArray1(bool verbose = false){ void LoadReactionParas(int arrayID, bool verbose = false){
//check is the transfer.root is using the latest reactionConfig.txt //check is the transfer.root is using the latest reactionConfig.txt
//sicne reaction.dat is generated as a by-product of transfer.root //sicne reaction.dat is generated as a by-product of transfer.root
@ -527,59 +569,86 @@ void LoadReactionParasForArray1(bool verbose = false){
// system("../Cleopatra/Transfer"); // system("../Cleopatra/Transfer");
// printf("########################## transfer.root updated\n"); // printf("########################## transfer.root updated\n");
//} //}
ReactionParas * reactParas = nullptr;
std::string fileName;
if( arrayID == 1){
reactParas = &AnalysisLib::reactParas1;
fileName = "reaction.dat";
}else if( arrayID == 2){
reactParas = &AnalysisLib::reactParas2;
fileName = "reaction2.dat";
}else{
printf("arrayID must be either 1 or 2. Abort.\n");
return;
}
reactParas->detPrepDist = AnalysisLib::detGeo.array1.detPerpDist;
printf(" loading reaction parameters"); printf(" loading reaction parameters");
std::ifstream file; std::ifstream file;
file.open("reaction.dat"); file.open(fileName.c_str());
hasReactionPara = false; reactParas->hasReactionPara = false;
if( file.is_open() ){ if( file.is_open() ){
std::string x; std::string x;
int i = 0; int i = 0;
while( file >> x ){ while( file >> x ){
if( x.substr(0,2) == "//" ) continue; if( x.substr(0,2) == "//" ) continue;
if( i == 0 ) mass = atof(x.c_str()); if( i == 0 ) reactParas->mass = atof(x.c_str());
if( i == 1 ) q = atof(x.c_str()); if( i == 1 ) reactParas->q = atof(x.c_str());
if( i == 2 ) betRel = atof(x.c_str()); if( i == 2 ) reactParas->beta = atof(x.c_str());
if( i == 3 ) Et = atof(x.c_str()); if( i == 3 ) reactParas->Et = atof(x.c_str());
if( i == 4 ) massB = atof(x.c_str()); if( i == 4 ) reactParas->massB = atof(x.c_str());
i = i + 1; i = i + 1;
} }
printf("........ done.\n"); printf("........ done.\n");
hasReactionPara = true; reactParas->hasReactionPara = true;
alpha = 299.792458 * abs(detGeo.Bfield) * q / TMath::TwoPi()/1000.; //MeV/mm reactParas->alpha = 299.792458 * abs(detGeo.Bfield) * reactParas->q / TMath::TwoPi()/1000.; //MeV/mm
gamm = 1./TMath::Sqrt(1-betRel*betRel); reactParas->gamma = 1./TMath::Sqrt(1-reactParas->beta * reactParas->beta);
G = alpha * gamm * betRel * detGeo.array1.detPerpDist ; reactParas->G = reactParas->alpha * reactParas->gamma * reactParas->beta * reactParas->detPrepDist ;
if( verbose ){ if( verbose ){
printf("\tmass-b : %f MeV/c2 \n", mass); printf("\tmass-b : %f MeV/c2 \n", reactParas->mass);
printf("\tcharge-b : %f \n", q); printf("\tcharge-b : %f \n", reactParas->q);
printf("\tE-total : %f MeV \n", Et); printf("\tE-total : %f MeV \n", reactParas->Et);
printf("\tmass-B : %f MeV/c2 \n", massB); printf("\tmass-B : %f MeV/c2 \n", reactParas->massB);
printf("\tbeta : %f \n", betRel); printf("\tbeta : %f \n", reactParas->beta);
printf("\tB-field : %f T \n", detGeo.Bfield); printf("\tB-field : %f T \n", detGeo.Bfield);
printf("\tslope : %f MeV/mm \n", alpha * betRel); printf("\tslope : %f MeV/mm \n", reactParas->alpha * reactParas->beta);
printf("\tdet radius: %f mm \n", detGeo.array1.detPerpDist); printf("\tdet radius: %f mm \n", reactParas->detPrepDist);
printf("\tG-coeff : %f MeV \n", G); printf("\tG-coeff : %f MeV \n", reactParas->G);
printf("=====================================================\n"); printf("=====================================================\n");
} }
}else{ }else{
printf("........ fail.\n"); printf("........ fail.\n");
hasReactionPara = false;
} }
file.close(); file.close();
} }
std::vector<double> CalExTheta(double e, double z){ std::vector<double> CalExTheta(double e, double z){
if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
ReactionParas * reactParas = nullptr;
if( detGeo.array1.zMin <= z && z <= detGeo.array1.zMax ){
reactParas = &reactParas1;
if( !reactParas->hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
}
if( detGeo.array2.zMin <= z && z <= detGeo.array2.zMax ){
reactParas = &reactParas2;
if( !reactParas->hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
}
double Ex = TMath::QuietNaN(); double Ex = TMath::QuietNaN();
double thetaCM = TMath::QuietNaN(); double thetaCM = TMath::QuietNaN();
double y = e + mass; // to give the KE + mass of proton; double y = e + reactParas->mass; // to give the KE + mass of proton;
double Z = alpha * gamm * betRel * z; double Z = reactParas->alpha * reactParas->gamma * reactParas->beta * z;
double H = TMath::Sqrt(TMath::Power(gamm * betRel,2) * (y*y - mass * mass) ) ; double H = TMath::Sqrt(TMath::Power(reactParas->gamma * reactParas->beta,2) * (y*y - reactParas->mass * reactParas->mass) ) ;
if( TMath::Abs(Z) < H ) { if( TMath::Abs(Z) < H ) {
///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi) ///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi)
@ -590,26 +659,26 @@ std::vector<double> CalExTheta(double e, double z){
int iter = 0; int iter = 0;
do{ do{
phi = nPhi; phi = nPhi;
nPhi = phi - (H * TMath::Sin(phi) - G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - G /TMath::Power( TMath::Cos(phi), 2)); nPhi = phi - (H * TMath::Sin(phi) - reactParas->G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - reactParas->G /TMath::Power( TMath::Cos(phi), 2));
iter ++; iter ++;
if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break; if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break;
}while( TMath::Abs(phi - nPhi ) > tolerrence); }while( TMath::Abs(phi - nPhi ) > tolerrence);
phi = nPhi; phi = nPhi;
/// check f'(phi) > 0 /// check f'(phi) > 0
double Df = H * TMath::Cos(phi) - G / TMath::Power( TMath::Cos(phi),2); double Df = H * TMath::Cos(phi) - reactParas->G / TMath::Power( TMath::Cos(phi),2);
if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){ if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){
double K = H * TMath::Sin(phi); double K = H * TMath::Sin(phi);
double x = TMath::ACos( mass / ( y * gamm - K)); double x = TMath::ACos( reactParas->mass / ( y * reactParas->gamma - K));
double momt = mass * TMath::Tan(x); /// momentum of particel b or B in CM frame double momt = reactParas->mass * TMath::Tan(x); /// momentum of particel b or B in CM frame
double EB = TMath::Sqrt(mass*mass + Et*Et - 2*Et*TMath::Sqrt(momt*momt + mass * mass)); double EB = TMath::Sqrt(reactParas->mass * reactParas->mass + reactParas->Et * reactParas->Et - 2 * reactParas->Et * TMath::Sqrt(momt*momt + reactParas->mass * reactParas->mass));
Ex = EB - massB; Ex = EB - reactParas->massB;
double hahaha1 = gamm* TMath::Sqrt(mass * mass + momt * momt) - y; double hahaha1 = reactParas->gamma * TMath::Sqrt(reactParas->mass * reactParas->mass + momt * momt) - y;
double hahaha2 = gamm* betRel * momt; double hahaha2 = reactParas->gamma * reactParas->beta * momt;
thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg(); thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg();
}
}
} }
return {Ex, thetaCM}; return {Ex, thetaCM};

0
working/DWBA → working/DWBA2
View File

0
working/Ex.txt → working/Ex1.txt
View File

6
working/Ex2.txt Normal file
View File

@ -0,0 +1,6 @@
//Ex relative_xsec SF sigma_in_MeV
//<--- use "//" for line comment
0.000 1.0 1.0 0.0100
//4.400 1.0 1.0 0.0100
//4.600 1.0 1.0 0.0100
#============_End_of_file

3
working/Monitor.C
View File

@ -183,7 +183,8 @@ void Monitor::Begin(TTree *tree){
AnalysisLib::LoadXScaleCorr(); AnalysisLib::LoadXScaleCorr();
AnalysisLib::LoadECorr(); AnalysisLib::LoadECorr();
AnalysisLib::LoadRDTCorr(); AnalysisLib::LoadRDTCorr();
AnalysisLib::LoadReactionParasForArray1(true); AnalysisLib::LoadReactionParas(1, true);
if( AnalysisLib::detGeo.use2ndArray ) AnalysisLib::LoadReactionParas(2, true);
if( (int) AnalysisLib::xnCorr.size() < mapping::NARRAY ) { isXNCorrOK = false; printf(" !!!!!!!! size of xnCorr < NARRAY .\n"); } if( (int) AnalysisLib::xnCorr.size() < mapping::NARRAY ) { isXNCorrOK = false; printf(" !!!!!!!! size of xnCorr < NARRAY .\n"); }
if( (int) AnalysisLib::xfxneCorr.size() < mapping::NARRAY ) { isXFXNCorrOK = false; printf(" !!!!!!!! size of xfxneCorr < NARRAY .\n"); } if( (int) AnalysisLib::xfxneCorr.size() < mapping::NARRAY ) { isXFXNCorrOK = false; printf(" !!!!!!!! size of xfxneCorr < NARRAY .\n"); }

2
working/detectorGeo.txt
View File

@ -26,7 +26,7 @@ Out //detector_facing_Out_or_In
235.8 //5th_det 235.8 //5th_det
290.0 290.0
#===============2nd_Array #===============2nd_Array
true //is_2nd_detctor_exist_false=false_other=true false //is_2nd_array_exist_is_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]

14
working/reactionConfig.txt → working/reactionConfig1.txt
View File

@ -1,16 +1,16 @@
32 //beam_A 32 //beam_A
14 //beam_Z 14 //beam_Z
2 //target_A 2 //target_A
1 //target_Z 1 //target_Z
1 //recoil_light_A 1 //recoil_light_A
1 //recoil-light_Z 1 //recoil-light_Z
8.8 //beam-energy_in_MeV/u 8.8 //beam-energy_in_MeV/u
0.000 //beam-energy_sigma_in_MeV/u 0.000 //beam-energy_sigma_in_MeV/u
0.000 //beam-angle_in_mrad 0.000 //beam-angle_in_mrad
0.000 //beam-emittance_in_mrad 0.000 //beam-emittance_in_mrad
0.00 //x_offset_of_Beam_in_mm 0.00 //x_offset_of_Beam_in_mm
0.00 //y_offset_of_Beam_in_mm 0.00 //y_offset_of_Beam_in_mm
100000 //number_of_Event_being_generated 10000 //number_of_Event_being_generated
false //isTargetScattering false //isTargetScattering
0.913 //target_density_in_g/cm3 0.913 //target_density_in_g/cm3
2.2e-4 //targetThickness_in_cm 2.2e-4 //targetThickness_in_cm
@ -18,8 +18,8 @@ false //isTargetScattering
../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil ../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil
../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil ../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil
false //isDacay false //isDacay
32 //decayNucleus_A 32 //decayNucleus_A
14 //decayNucleus_Z 14 //decayNucleus_Z
false //isReDo false //isReDo
0.0 //excitation_energy_of_A[MeV] 0.0 //List_of_excitation_energy_of_A[MeV]
#===== end of file #===== end of file

25
working/reactionConfig2.txt Normal file
View File

@ -0,0 +1,25 @@
16 //beam_A
8 //beam_Z
2 //target_A
1 //target_Z
1 //recoil_light_A
1 //recoil-light_Z
10.0 //beam-energy_in_MeV/u
0.000 //beam-energy_sigma_in_MeV/u
0.000 //beam-angle_in_mrad
0.000 //beam-emittance_in_mrad
0.00 //x_offset_of_Beam_in_mm
0.00 //y_offset_of_Beam_in_mm
100000 //number_of_Event_being_generated
false //isTargetScattering
0.913 //target_density_in_g/cm3
2.2e-4 //targetThickness_in_cm
../SRIM/20F_in_CD2.txt //stopping_power_for_beam
../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil
../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil
false //isDacay
32 //decayNucleus_A
14 //decayNucleus_Z
false //isReDo
0.0 //List_of_excitation_energy_of_A[MeV]
#===== end of file