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

bug fix on ClassHelios

Browse Source
This commit is contained in:
Ryan Tang 2024-02-14 16:30:13 -05:00
parent 871a4d0b26
commit fca3602769
6 changed files with 197 additions and 156 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Armory/ClassDetGeo.h
View File

@ -30,13 +30,14 @@ struct Array{
double zMin, zMax; double zMin, zMax;
void DeduceAbsolutePos(){ void DeduceAbsolutePos(){
nDet = pos.size(); nDet = pos.size();
detPos.clear(); detPos.clear();
for(int id = 0; id < nDet; id++){ for(int id = 0; id < nDet; id++){
if( firstPos > 0 ) detPos.push_back(firstPos + pos[id]); if( firstPos > 0 ) detPos.push_back(firstPos + pos[id]);
if( firstPos < 0 ) detPos.push_back(firstPos - pos[nDet - 1 - id]); if( firstPos < 0 ) detPos.push_back(firstPos - pos[nDet - 1 - id]);
///printf("%d | %f, %f \n", id, pos[id], detPos[id]); // printf("%d | %f, %f \n", id, pos[id], detPos[id]);
} }
zMin = TMath::Min(detPos.front(), detPos.back()) - (firstPos < 0 ? detLength : 0); zMin = TMath::Min(detPos.front(), detPos.back()) - (firstPos < 0 ? detLength : 0);
@ -65,8 +66,10 @@ struct Array{
class DetGeo { class DetGeo {
public: public:
DetGeo(); DetGeo(){};
~DetGeo(); DetGeo(TString detGeoTxt){ LoadDetectorGeo(detGeoTxt, false);}
DetGeo(TMacro * macro){ LoadDetectorGeo(macro, false);}
~DetGeo(){};
double Bfield; /// T double Bfield; /// T
int BfieldSign ; /// sign of B-field int BfieldSign ; /// sign of B-field
@ -100,24 +103,19 @@ private:
}; };
inline DetGeo::DetGeo(){
}
inline DetGeo::~DetGeo(){
}
inline bool DetGeo::LoadDetectorGeo(TString fileName, bool verbose){ inline bool DetGeo::LoadDetectorGeo(TString fileName, bool verbose){
TMacro * haha = new TMacro(); TMacro * haha = new TMacro();
if( haha->ReadFile(fileName) > 0 ) { if( haha->ReadFile(fileName) > 0 ) {
if( LoadDetectorGeo(haha, verbose) ){ if( LoadDetectorGeo(haha, verbose) ){
delete haha;
return true; return true;
}else{ }else{
delete haha;
return false; return false;
} }
}else{ }else{
delete haha;
return false; return false;
} }
} }
@ -142,7 +140,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
std::vector<std::string> str = AnalysisLib::SplitStr(macro->GetListOfLines()->At(i)->GetName(), " "); std::vector<std::string> str = AnalysisLib::SplitStr(macro->GetListOfLines()->At(i)->GetName(), " ");
//printf("%3d | %s\n", i, str[0].c_str()); // printf("%3d | %s\n", i, str[0].c_str());
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 ) {

5
Armory/ClassReactionConfig.h
View File

@ -17,6 +17,8 @@ class ReactionConfig {
public: public:
ReactionConfig(){} ReactionConfig(){}
ReactionConfig(TString configTxt){ LoadReactionConfig(configTxt);}
ReactionConfig(TMacro * macro){ LoadReactionConfig(macro);}
~ReactionConfig(){} ~ReactionConfig(){}
int beamA, beamZ; int beamA, beamZ;
@ -83,11 +85,14 @@ inline bool ReactionConfig::LoadReactionConfig(TString fileName){
TMacro * haha = new TMacro(); TMacro * haha = new TMacro();
if( haha->ReadFile(fileName) > 0 ) { if( haha->ReadFile(fileName) > 0 ) {
if( LoadReactionConfig(haha) ){ if( LoadReactionConfig(haha) ){
delete haha;
return true; return true;
}else{ }else{
delete haha;
return false; return false;
} }
}else{ }else{
delete haha;
return false; return false;
} }
} }

152
Cleopatra/ClassHelios.h
View File

@ -35,26 +35,27 @@ struct trajectory{
int detID, detRowID; int detID, detRowID;
int loop; int loop;
double effLoop; double effLoop;
void PrintTrajectory(){
printf("=====================\n");
printf(" theta : %f deg\n", theta*TMath::RadToDeg());
printf(" phi : %f deg\n", phi*TMath::RadToDeg());
printf(" vt : %f mm/ns\n", vt);
printf(" vp : %f mm/ns\n", vp);
printf(" rho : %f mm\n", rho);
printf(" z0 : %f mm\n", z0);
printf(" t0 : %f ns\n", t0);
printf("(x, y, z) : (%f, %f. %f) mm\n", x, y, z);
printf(" R : %f mm\n", R);
printf(" t : %f ns\n", t);
printf(" effLoop : %f cycle\n", effLoop);
printf(" Loop : %d cycle\n", loop);
printf(" detRowID : %d \n", detRowID);
printf(" detID : %d \n", detID);
}
}; };
void PrintTrajectory(trajectory a){
printf("=====================\n");
printf(" theta : %f deg\n", a.theta*TMath::RadToDeg());
printf(" phi : %f deg\n", a.phi*TMath::RadToDeg());
printf(" vt : %f mm/ns\n", a.vt);
printf(" vp : %f mm/ns\n", a.vp);
printf(" rho : %f mm\n", a.rho);
printf(" z0 : %f mm\n", a.z0);
printf(" t0 : %f ns\n", a.t0);
printf("(x, y, z) : (%f, %f. %f) mm\n", a.x, a.y, a.z);
printf(" R : %f mm\n", a.R);
printf(" t : %f ns\n", a.t);
printf(" effLoop : %f cycle\n", a.effLoop);
printf(" Loop : %d cycle\n", a.loop);
printf(" detRowID : %d \n", a.detRowID);
printf(" detID : %d \n", a.detID);
}
class HELIOS{ class HELIOS{
public: public:
@ -86,11 +87,11 @@ public:
} }
int DetAcceptance(); int DetAcceptance();
int CalArrayHit(TLorentzVector Pb, int Zb); int CalArrayHit(TLorentzVector Pb, int Zb, bool debug = false);
int CalRecoilHit(TLorentzVector PB, int ZB); int CalRecoilHit(TLorentzVector PB, int ZB);
//int CalHit(TLorentzVector Pb, int Zb, TLorentzVector PB, int ZB, double xOff = 0, double yOff = 0 ); // return 0 for no hit, 1 for hit //int CalHit(TLorentzVector Pb, int Zb, TLorentzVector PB, int ZB, double xOff = 0, double yOff = 0 ); // return 0 for no hit, 1 for hit
void CalTrajectoryPara(TLorentzVector P, int Z, int id); // id = 0 for Pb, id = 1 for PB. void CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil);
int GetNumberOfDetectorsInSamePos(){return array.mDet;} int GetNumberOfDetectorsInSamePos(){return array.mDet;}
double GetEnergy(){return e;} double GetEnergy(){return e;}
@ -98,25 +99,28 @@ public:
/// clockwise rotation for B-field along the z-axis, sign = 1. /// clockwise rotation for B-field along the z-axis, sign = 1.
double XPos(double Zpos, double theta, double phi, double rho, int sign){ double XPos(double Zpos, double theta, double phi, double rho, int sign){
if( TMath::IsNaN(Zpos) ) return TMath::QuietNaN();
return rho * ( TMath::Sin( TMath::Tan(theta) * Zpos / rho - sign * phi ) + sign * TMath::Sin(phi) ) + xOff; return rho * ( TMath::Sin( TMath::Tan(theta) * Zpos / rho - sign * phi ) + sign * TMath::Sin(phi) ) + xOff;
} }
double YPos(double Zpos, double theta, double phi, double rho, int sign){ double YPos(double Zpos, double theta, double phi, double rho, int sign){
if( TMath::IsNaN(Zpos) ) return TMath::QuietNaN();
return rho * sign * (TMath::Cos( TMath::Tan(theta) * Zpos / rho - sign * phi ) - TMath::Cos(phi)) + yOff; return rho * sign * (TMath::Cos( TMath::Tan(theta) * Zpos / rho - sign * phi ) - TMath::Cos(phi)) + yOff;
} }
double RPos(double Zpos, double theta, double phi, double rho, int sign){ double RPos(double Zpos, double theta, double phi, double rho, int sign){
double x = XPos(Zpos, theta, phi, rho, sign) ; if( TMath::IsNaN(Zpos) ) return TMath::QuietNaN();
double y = YPos(Zpos, theta, phi, rho, sign) ; double x = XPos(Zpos, theta, phi, rho, sign) ;
return sqrt(x*x+y*y); double y = YPos(Zpos, theta, phi, rho, sign) ;
return sqrt(x*x+y*y);
} }
double GetXPos(double ZPos){ return XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } double GetXPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : XPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); }
double GetYPos(double ZPos){ return YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } double GetYPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : YPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); }
double GetR(double ZPos) { return RPos( ZPos, orbitb.theta, orbitb.phi, orbitb.rho, detGeo.BfieldSign); } double GetR(double ZPos) { return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : 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, detGeo.BfieldSign); } double GetRecoilXPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : XPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); }
double GetRecoilYPos(double ZPos){ return YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } double GetRecoilYPos(double ZPos){ return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : YPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); }
double GetRecoilR(double ZPos) { return RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); } double GetRecoilR(double ZPos) { return TMath::IsNaN(ZPos) ? TMath::QuietNaN() : RPos( ZPos, orbitB.theta, orbitB.phi, orbitB.rho, detGeo.BfieldSign); }
double GetBField() {return detGeo.Bfield;} double GetBField() {return detGeo.Bfield;}
double GetDetRadius() {return array.detPerpDist;} double GetDetRadius() {return array.detPerpDist;}
@ -125,6 +129,9 @@ public:
trajectory GetTrajectory_B() {return orbitB;} trajectory GetTrajectory_B() {return orbitB;}
DetGeo GetDetectorGeometry() {return detGeo;} DetGeo GetDetectorGeometry() {return detGeo;}
TString GetHitMessage() const {return hitMessage;}
TString GetAcceptanceMessage() const {return accMessage;}
private: private:
@ -156,10 +163,12 @@ private:
double eB; ///energy of heavy recoil double eB; ///energy of heavy recoil
bool isDetReady; bool isDetReady;
TString hitMessage;
TString accMessage; //acceptance check
double xOff, yOff; // beam position double xOff, yOff; // beam position
bool overrideDetDistance; bool overrideDetDistance;
bool overrideFirstPos; bool overrideFirstPos;
bool isCoincidentWithRecoil; bool isCoincidentWithRecoil;
@ -182,6 +191,9 @@ HELIOS::HELIOS(){
isDetReady = false; isDetReady = false;
hitMessage = "";
accMessage = "";
overrideDetDistance = false; overrideDetDistance = false;
overrideFirstPos = false; overrideFirstPos = false;
isCoincidentWithRecoil = false; isCoincidentWithRecoil = false;
@ -220,40 +232,43 @@ 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( array.firstPos < 0 && orbitb.z > 0 ) return -1; if( array.firstPos < 0 && orbitb.z > 0 ) {accMessage = "array at upstream, z is downstream."; return -1;}
if( array.firstPos > 0 && orbitb.z < 0 ) return -1; if( array.firstPos > 0 && orbitb.z < 0 ) {accMessage = "array at downstream, z is upstream."; return -1;}
// -11 ======== rho is too small // -11 ======== rho is too small
if( 2 * orbitb.rho < array.detPerpDist ) return -11; if( 2 * orbitb.rho < array.detPerpDist ) { accMessage = "rho is too small"; 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 ) {accMessage = "det Row ID == -1"; return -15;}
// -10 =========== when rho is too big . // -10 =========== when rho is too big .
if( detGeo.bore < 2 * orbitb.rho) return -10; if( detGeo.bore < 2 * orbitb.rho) { accMessage = "rho is too big"; return -10;}
// -14 ========== check particle-B hit radius on recoil dectector // -14 ========== check particle-B hit radius on recoil dectector
if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) return -14; if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) {
accMessage = "heavy recoil does not hit recoil detector";
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(detGeo.recoilPos); rhoHit = GetR(detGeo.recoilPos);
if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) return -12; if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { accMessage = "light recoil blocked by recoil detector"; return -12;}
if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) return -12; if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { accMessage = "light recoil blocked by recoil detector"; return -12;}
// -13 ========= not more than 3 loops // -13 ========= not more than 3 loops
if( orbitb.loop > 3 ) return -13; if( orbitb.loop > 3 ) {accMessage = "more than 3 loops."; return -13;}
// -2 ========= calculate the "y"-distance from detector center // -2 ========= calculate the "y"-distance from detector center
if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) return -2; if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) { accMessage = "hit at the XY gap."; 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( array.firstPos < 0 && orbitb.z < array.detPos[0] - array.detLength ) return -3; if( array.firstPos < 0 && orbitb.z < array.detPos[0] - array.detLength ) { accMessage = "hit more upstream than the array length"; return -3; }
if( array.firstPos > 0 && orbitb.z > array.detPos[array.nDet-1] + array.detLength ) return -3; if( array.firstPos > 0 && orbitb.z > array.detPos[array.nDet-1] + array.detLength ) { accMessage = "hit more downstream than the array length"; return -3;}
// -4 ======== Hit on blacker // -4 ======== Hit on blacker
if( array.blocker != 0 && array.firstPos > 0 && array.detPos[0] - array.blocker < orbitb.z && orbitb.z < array.detPos[0] ) return -4; if( array.blocker != 0 && array.firstPos > 0 && array.detPos[0] - array.blocker < orbitb.z && orbitb.z < array.detPos[0] ) { accMessage = "hit 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; if( array.blocker != 0 && array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z && orbitb.z < array.detPos[array.nDet-1] + array.blocker ) { accMessage = "hit 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;
@ -264,6 +279,7 @@ int HELIOS::DetAcceptance(){
orbitb.effLoop += 1.0; orbitb.effLoop += 1.0;
orbitb.loop += 1; orbitb.loop += 1;
orbitb.t = orbitb.t0 * orbitb.effLoop; orbitb.t = orbitb.t0 * orbitb.effLoop;
accMessage = " hit less than the nearest array. increase loop ";
return 2; return 2;
} }
@ -273,6 +289,7 @@ int HELIOS::DetAcceptance(){
if( array.detPos[i] - array.detLength <= orbitb.z && orbitb.z <= array.detPos[i]) { if( array.detPos[i] - array.detLength <= orbitb.z && orbitb.z <= array.detPos[i]) {
orbitb.detID = i; orbitb.detID = i;
detX = ( orbitb.z - (array.detPos[i] + array.detLength/2 ))/ array.detLength * 2 ;// range from -1 , 1 detX = ( orbitb.z - (array.detPos[i] + array.detLength/2 ))/ array.detLength * 2 ;// range from -1 , 1
accMessage = "hit array";
return 1; return 1;
} }
} }
@ -282,6 +299,7 @@ int HELIOS::DetAcceptance(){
///printf(" %d | %f < z = %f < %f \n", i, array.detPos[i], orbitb.z, array.detPos[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 - (array.detPos[i] - array.detLength/2 ))/ array.detLength*2 ;// range from -1 , 1 detX = ( orbitb.z - (array.detPos[i] - array.detLength/2 ))/ array.detLength*2 ;// range from -1 , 1
accMessage = "hit array";
return 1; return 1;
} }
} }
@ -291,11 +309,11 @@ int HELIOS::DetAcceptance(){
// -5 ======== check hit array gap // -5 ======== check hit array gap
if( array.firstPos < 0 ){ if( array.firstPos < 0 ){
for( int i = 0; i < array.nDet-1 ; i++){ for( int i = 0; i < array.nDet-1 ; i++){
if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) return -5; //increaseLoopFlag = 3; if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) { accMessage = "hit array Z-gap"; return -5; }//increaseLoopFlag = 3;
} }
}else{ }else{
for( int i = 0; i < array.nDet-1 ; i++){ for( int i = 0; i < array.nDet-1 ; i++){
if( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) return -5; //increaseLoopFlag = 3; if( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) { accMessage = "hit array Z-gap"; return -5; }//increaseLoopFlag = 3;
} }
} }
if (increaseLoopFlag == 3 ) { if (increaseLoopFlag == 3 ) {
@ -303,16 +321,17 @@ int HELIOS::DetAcceptance(){
orbitb.effLoop += 1.0; orbitb.effLoop += 1.0;
orbitb.loop += 1; orbitb.loop += 1;
orbitb.t = orbitb.t0 * orbitb.effLoop; orbitb.t = orbitb.t0 * orbitb.effLoop;
accMessage = " try one more loop. ";
return 3; return 3;
} }
accMessage = " unknown reason ";
return -20; // for unknown reason return -20; // for unknown reason
} }
void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, bool isLightRecoil){
if( id == 0 ){ if( isLightRecoil ){
orbitb.theta = P.Theta(); orbitb.theta = P.Theta();
orbitb.phi = P.Phi(); orbitb.phi = P.Phi();
orbitb.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm orbitb.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm
@ -323,10 +342,8 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){
orbitb.detID = -1; orbitb.detID = -1;
orbitb.detRowID = -1; orbitb.detRowID = -1;
} }else{
if( id == 1 ){
orbitB.theta = P.Theta(); orbitB.theta = P.Theta();
orbitB.phi = P.Phi(); orbitB.phi = P.Phi();
orbitB.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm orbitB.rho = P.Pt() / abs(detGeo.Bfield) / Z / c * 1000; //mm
@ -337,29 +354,27 @@ void HELIOS::CalTrajectoryPara(TLorentzVector P, int Z, int id){
orbitB.detID = -1; orbitB.detID = -1;
orbitB.detRowID = -1; orbitB.detRowID = -1;
} }
} }
int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb, bool debug){
e = Pb.E() - Pb.M(); e = Pb.E() - Pb.M();
detX = TMath::QuietNaN(); detX = TMath::QuietNaN();
rhoHit = TMath::QuietNaN(); rhoHit = TMath::QuietNaN();
CalTrajectoryPara(Pb, Zb, 0); CalTrajectoryPara(Pb, Zb, true);
int targetLoop = 1; int targetLoop = 1;
int inOut = array.detFaceOut == 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;
if( debug ) { if( debug ) {
printf("===================================\n"); printf("===================================\n");
printf("theta : %f deg, phi : %f deg \n", orbitb.theta * TMath::RadToDeg(), orbitb.phi * TMath::RadToDeg()); printf("theta : %f deg, phi : %f deg \n", orbitb.theta * TMath::RadToDeg(), orbitb.phi * TMath::RadToDeg());
printf("z0: %f mm, rho : %f mm \n", orbitb.z0, orbitb.rho); printf("z0: %f mm, rho : %f mm \n", orbitb.z0, orbitb.rho);
printf(" inOut : %d = %s \n", inOut, inOut == 1 ? "Out" : "in"); printf(" inOut : %d = %s \n", inOut, inOut == 1 ? "Out" : "in");
printf(" z range : %.2f - %.2f \n", detGeo.zMin, detGeo.zMax); printf(" z range : %.2f - %.2f \n", detGeo.zMin, detGeo.zMax);
printf(" B-field sign : %d\n", detGeo.BfieldSign);
printf("-----------------------------------\n"); printf("-----------------------------------\n");
} }
@ -373,7 +388,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
double phiD = TMath::TwoPi()/array.mDet * i ; double phiD = TMath::TwoPi()/array.mDet * i ;
double dphi = orbitb.phi - phiD; double dphi = orbitb.phi - phiD;
double aEff = array.detPerpDist - (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 - detGeo.BfieldTheta * sin(dphi)); double hahaha = asin( aEff/ orbitb.rho - detGeo.BfieldSign * sin(dphi));
int n = 2*targetLoop + inOut; int n = 2*targetLoop + inOut;
@ -390,22 +405,12 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
} }
///Selection ///Selection
if( !TMath::IsNaN(zP) && 0< zP/orbitb.z0 && TMath::Max(0, targetLoop-1) < zP/orbitb.z0 && zP/orbitb.z0 < targetLoop ) { if( !TMath::IsNaN(zP) && 0 < zP/orbitb.z0 && TMath::Max(0, targetLoop-1) < zP/orbitb.z0 && zP/orbitb.z0 < targetLoop ) {
zPossible.push_back(zP); zPossible.push_back(zP);
dID.push_back(i); dID.push_back(i);
} }
} }
/*
if( zPossible.size() == 0 ){ // will not happen
zHit = TMath::QuietNaN();
xPos = TMath::QuietNaN();
yPos = TMath::QuietNaN();
loop = -1;
detID = -1;
detRowID = -1;
return -1 ;
}*/
if( debug ) printf("-----------------------------------\n"); if( debug ) printf("-----------------------------------\n");
double dMin = 1; double dMin = 1;
for( int i = 0; i < (int) zPossible.size(); i++){ for( int i = 0; i < (int) zPossible.size(); i++){
@ -434,6 +439,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
orbitb.z = TMath::QuietNaN(); orbitb.z = TMath::QuietNaN();
orbitb.loop = -1; orbitb.loop = -1;
orbitb.detRowID = -1; orbitb.detRowID = -1;
hitMessage = "wrong direction.";
return - 2; return - 2;
} }
@ -447,6 +453,7 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
orbitb.z = TMath::QuietNaN(); orbitb.z = TMath::QuietNaN();
orbitb.loop = -1; orbitb.loop = -1;
orbitb.detRowID = -1; orbitb.detRowID = -1;
hitMessage = "not on the detector plan.";
return -3; return -3;
} }
} }
@ -459,12 +466,13 @@ int HELIOS::CalArrayHit(TLorentzVector Pb, int Zb){
orbitb.y = GetYPos(orbitb.z) ; orbitb.y = GetYPos(orbitb.z) ;
orbitb.R = GetR(orbitb.z); orbitb.R = GetR(orbitb.z);
hitMessage = "successful hit.";
return 1; // return 1 when OK return 1; // return 1 when OK
} }
int HELIOS::CalRecoilHit(TLorentzVector PB, int ZB){ int HELIOS::CalRecoilHit(TLorentzVector PB, int ZB){
CalTrajectoryPara(PB, ZB, 1); CalTrajectoryPara(PB, ZB, false);
orbitB.z = detGeo.recoilPos; orbitB.z = detGeo.recoilPos;
orbitB.x = GetRecoilXPos(detGeo.recoilPos) ; orbitB.x = GetRecoilXPos(detGeo.recoilPos) ;

94
Cleopatra/ClassTransfer.h
View File

@ -21,7 +21,20 @@
//======================================================= //=======================================================
class TransferReaction { class TransferReaction {
public: public:
TransferReaction(); TransferReaction(){Inititization();};
TransferReaction(int beamA, int beamZ,
int targetA, int targetZ,
int recoilA, int recoilZ, float beamEnergy_AMeV){
Inititization();
SetReactionSimple(beamA, beamZ,
targetA, targetZ,
recoilA, recoilZ, beamEnergy_AMeV);
}
TransferReaction(string reactionConfigFile){
Inititization();
SetReactionFromFile(reactionConfigFile);
}
~TransferReaction(); ~TransferReaction();
void SetA(int A, int Z, double Ex); void SetA(int A, int Z, double Ex);
@ -41,15 +54,13 @@ public:
TString GetReactionName(); TString GetReactionName();
TString GetReactionName_Latex(); TString GetReactionName_Latex();
ReactionConfig GetRectionConfig() { return reaction;} ReactionConfig GetRectionConfig() { return reactionConfig;}
double GetMass_A() const {return mA + ExA;} double GetMass_A() const {return mA + ExA;}
double GetMass_a() const {return ma;} double GetMass_a() const {return ma;}
double GetMass_b() const {return mb;} double GetMass_b() const {return mb;}
double GetMass_B() const {return mB + ExB;} double GetMass_B() const {return mB + ExB;}
int GetCharge_b() const {return reaction.recoilLightZ;}
double GetCMTotalKE() {return Etot - mA - ma;} double GetCMTotalKE() {return Etot - mA - ma;}
double GetQValue() {return mA + ExA + ma - mb - mB - ExB;} double GetQValue() {return mA + ExA + ma - mb - mB - ExB;}
double GetMaxExB() {return Etot - mb - mB;} double GetMaxExB() {return Etot - mb - mB;}
@ -59,6 +70,8 @@ public:
TLorentzVector GetPb(){return Pb;} TLorentzVector GetPb(){return Pb;}
TLorentzVector GetPB(){return PB;} TLorentzVector GetPB(){return PB;}
void PrintFourVectors();
void CalReactionConstant(); void CalReactionConstant();
void Event(double thetaCM_rad, double phiCM_rad); void Event(double thetaCM_rad, double phiCM_rad);
@ -72,7 +85,7 @@ public:
private: private:
ReactionConfig reaction; ReactionConfig reactionConfig;
string nameA, namea, nameb, nameB; string nameA, namea, nameb, nameB;
double thetaIN, phiIN; double thetaIN, phiIN;
@ -92,10 +105,12 @@ private:
TLorentzVector PA, Pa, Pb, PB; TLorentzVector PA, Pa, Pb, PB;
TString format(TString name); TString format(TString name);
void Inititization();
}; };
TransferReaction::TransferReaction(){ void TransferReaction::Inititization(){
thetaIN = 0.; thetaIN = 0.;
phiIN = 0.; phiIN = 0.;
@ -104,7 +119,7 @@ TransferReaction::TransferReaction(){
Setb(1,1); Setb(1,1);
SetB(13,6); SetB(13,6);
TA = 6; TA = 6;
T = TA * reaction.beamA; T = TA * reactionConfig.beamA;
ExA = 0; ExA = 0;
ExB = 0; ExB = 0;
@ -126,8 +141,8 @@ TransferReaction::~TransferReaction(){
void TransferReaction::SetA(int A, int Z, double Ex = 0){ void TransferReaction::SetA(int A, int Z, double Ex = 0){
Isotope temp (A, Z); Isotope temp (A, Z);
mA = temp.Mass; mA = temp.Mass;
reaction.beamA = A; reactionConfig.beamA = A;
reaction.beamZ = Z; reactionConfig.beamZ = Z;
ExA = Ex; ExA = Ex;
nameA = temp.Name; nameA = temp.Name;
isReady = false; isReady = false;
@ -137,8 +152,8 @@ void TransferReaction::SetA(int A, int Z, double Ex = 0){
void TransferReaction::Seta(int A, int Z){ void TransferReaction::Seta(int A, int Z){
Isotope temp (A, Z); Isotope temp (A, Z);
ma = temp.Mass; ma = temp.Mass;
reaction.targetA = A; reactionConfig.targetA = A;
reaction.targetZ = Z; reactionConfig.targetZ = Z;
namea = temp.Name; namea = temp.Name;
isReady = false; isReady = false;
isBSet = false; isBSet = false;
@ -146,8 +161,8 @@ void TransferReaction::Seta(int A, int Z){
void TransferReaction::Setb(int A, int Z){ void TransferReaction::Setb(int A, int Z){
Isotope temp (A, Z); Isotope temp (A, Z);
mb = temp.Mass; mb = temp.Mass;
reaction.recoilLightA = A; reactionConfig.recoilLightA = A;
reaction.recoilLightZ = Z; reactionConfig.recoilLightZ = Z;
nameb = temp.Name; nameb = temp.Name;
isReady = false; isReady = false;
isBSet = false; isBSet = false;
@ -155,8 +170,8 @@ void TransferReaction::Setb(int A, int Z){
void TransferReaction::SetB(int A, int Z){ void TransferReaction::SetB(int A, int Z){
Isotope temp (A, Z); Isotope temp (A, Z);
mB = temp.Mass; mB = temp.Mass;
reaction.recoilHeavyA = A; reactionConfig.recoilHeavyA = A;
reaction.recoilHeavyZ = Z; reactionConfig.recoilHeavyZ = Z;
nameB = temp.Name; nameB = temp.Name;
isReady = false; isReady = false;
isBSet = true; isBSet = true;
@ -164,7 +179,7 @@ void TransferReaction::SetB(int A, int Z){
void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){ void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){
this->TA = KEA; this->TA = KEA;
this->T = TA * reaction.beamA; this->T = TA * reactionConfig.beamA;
this->thetaIN = theta; this->thetaIN = theta;
this->phiIN = phi; this->phiIN = phi;
isReady = false; isReady = false;
@ -174,15 +189,15 @@ void TransferReaction::SetReactionSimple(int beamA, int beamZ,
int targetA, int targetZ, int targetA, int targetZ,
int recoilA, int recoilZ, float beamEnergy_AMeV){ int recoilA, int recoilZ, float beamEnergy_AMeV){
reaction.SetReactionSimple(beamA, beamZ, reactionConfig.SetReactionSimple(beamA, beamZ,
targetA, targetZ, targetA, targetZ,
recoilA, recoilZ, beamEnergy_AMeV); recoilA, recoilZ, beamEnergy_AMeV);
SetA(reaction.beamA, reaction.beamZ); SetA(reactionConfig.beamA, reactionConfig.beamZ);
Seta(reaction.targetA, reaction.targetZ); Seta(reactionConfig.targetA, reactionConfig.targetZ);
Setb(reaction.recoilLightA, reaction.recoilLightZ); Setb(reactionConfig.recoilLightA, reactionConfig.recoilLightZ);
SetB(reaction.recoilHeavyA, reaction.recoilHeavyZ); SetB(reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ);
SetIncidentEnergyAngle(reaction.beamEnergy, 0, 0); SetIncidentEnergyAngle(reactionConfig.beamEnergy, 0, 0);
CalReactionConstant(); CalReactionConstant();
@ -200,13 +215,13 @@ void TransferReaction::SetExB(double Ex){
void TransferReaction::SetReactionFromFile(string reactionConfigFile){ void TransferReaction::SetReactionFromFile(string reactionConfigFile){
if( reaction.LoadReactionConfig(reactionConfigFile) ){ if( reactionConfig.LoadReactionConfig(reactionConfigFile) ){
SetA(reaction.beamA, reaction.beamZ); SetA(reactionConfig.beamA, reactionConfig.beamZ);
Seta(reaction.targetA, reaction.targetZ); Seta(reactionConfig.targetA, reactionConfig.targetZ);
Setb(reaction.recoilLightA, reaction.recoilLightZ); Setb(reactionConfig.recoilLightA, reactionConfig.recoilLightZ);
SetB(reaction.recoilHeavyA, reaction.recoilHeavyZ); SetB(reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ);
SetIncidentEnergyAngle(reaction.beamEnergy, 0, 0); SetIncidentEnergyAngle(reactionConfig.beamEnergy, 0, 0);
CalReactionConstant(); CalReactionConstant();
@ -245,9 +260,9 @@ TString TransferReaction::GetReactionName_Latex(){
void TransferReaction::CalReactionConstant(){ void TransferReaction::CalReactionConstant(){
if( !isBSet){ if( !isBSet){
reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA; reactionConfig.recoilHeavyA = reactionConfig.beamA + reactionConfig.targetA - reactionConfig.recoilLightA;
reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ; reactionConfig.recoilHeavyZ = reactionConfig.beamZ + reactionConfig.targetZ - reactionConfig.recoilLightZ;
Isotope temp (reaction.recoilHeavyA, reaction.recoilHeavyZ); Isotope temp (reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ);
mB = temp.Mass; mB = temp.Mass;
isBSet = true; isBSet = true;
} }
@ -267,6 +282,21 @@ void TransferReaction::CalReactionConstant(){
isReady = true; isReady = true;
} }
void TransferReaction::PrintFourVectors(){
printf("A : %10.2f %10.2f %10.2f %10.2f\n", PA.E(), PA.Px(), PA.Py(), PA.Pz());
printf("a : %10.2f %10.2f %10.2f %10.2f\n", Pa.E(), Pa.Px(), Pa.Py(), Pa.Pz());
printf("b : %10.2f %10.2f %10.2f %10.2f\n", Pb.E(), Pb.Px(), Pb.Py(), Pb.Pz());
printf("B : %10.2f %10.2f %10.2f %10.2f\n", PB.E(), PB.Px(), PB.Py(), PB.Pz());
printf("-------------------------------------------------------\n");
printf("B : %10.2f %10.2f %10.2f %10.2f\n",
PA.E() + Pa.E() - Pb.E() - PB.E(),
PA.Px() + Pa.Px() - Pb.Px() - PB.Px(),
PA.Py() + Pa.Py() - Pb.Py() - PB.Py(),
PA.Pz() + Pa.Pz() - Pb.Pz() - PB.Pz());
}
void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){ void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){
if( isReady == false ){ if( isReady == false ){
@ -325,7 +355,7 @@ std::pair<double, double> TransferReaction::CalExThetaCM(double e, double z, dou
double mass = mb; double mass = mb;
double massB = mB; double massB = mB;
double y = e + mass; double y = e + mass;
double slope = 299.792458 * reaction.recoilLightZ * abs(Bfield) / TMath::TwoPi() * beta / 1000.; // MeV/mm; double slope = 299.792458 * reactionConfig.recoilLightZ * abs(Bfield) / TMath::TwoPi() * beta / 1000.; // MeV/mm;
double alpha = slope/beta; double alpha = slope/beta;
double G = alpha * gamma * beta * perpDist ; double G = alpha * gamma * beta * perpDist ;
double Z = alpha * gamma * beta * z; double Z = alpha * gamma * beta * z;

76
Cleopatra/Transfer.h
View File

@ -32,16 +32,16 @@ void Transfer(
TString filename = "reaction.dat"){ /// when no file, no output. TString filename = "reaction.dat"){ /// when no file, no output.
//############################################# Set Reaction //############################################# Set Reaction
TransferReaction reaction; TransferReaction transfer;
reaction.SetReactionFromFile(basicConfig); transfer.SetReactionFromFile(basicConfig);
printf("*****************************************************************\n"); printf("*****************************************************************\n");
printf("*\e[1m\e[33m %27s \e[0m*\n", reaction.GetReactionName().Data()); printf("*\e[1m\e[33m %27s \e[0m*\n", transfer.GetReactionName().Data());
printf("*****************************************************************\n"); printf("*****************************************************************\n");
printf("----- loading reaction setting from %s. \n", basicConfig.c_str()); printf("----- loading reaction setting from %s. \n", basicConfig.c_str());
printf("\e[32m#################################### Beam \e[0m\n"); printf("\e[32m#################################### Beam \e[0m\n");
const ReactionConfig reactionConfig = reaction.GetRectionConfig(); const ReactionConfig reactionConfig = transfer.GetRectionConfig();
reactionConfig.PrintReactionConfig(); reactionConfig.PrintReactionConfig();
@ -56,12 +56,12 @@ void Transfer(
const DetGeo detGeo = helios.GetDetectorGeometry(); const DetGeo detGeo = helios.GetDetectorGeometry();
printf("==================================== E-Z plot slope\n"); printf("==================================== E-Z plot slope\n");
double betaRect = reaction.GetReactionBeta() ; double betaRect = transfer.GetReactionBeta() ;
double gamma = reaction.GetReactionGamma(); double gamma = transfer.GetReactionGamma();
double mb = reaction.GetMass_b(); double mb = transfer.GetMass_b();
double pCM = reaction.GetMomentumbCM(); double pCM = transfer.GetMomentumbCM();
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
double slope = 299.792458 * reaction.GetCharge_b() * abs(helios.GetBField()) / TMath::TwoPi() * betaRect / 1000.; /// MeV/mm double slope = 299.792458 * reactionConfig.recoilLightZ * abs(helios.GetBField()) / TMath::TwoPi() * betaRect / 1000.; /// MeV/mm
printf(" e-z slope : %f MeV/mm\n", slope); printf(" e-z slope : %f MeV/mm\n", slope);
double intercept = q/gamma - mb; // MeV double intercept = q/gamma - mb; // MeV
printf(" e-z intercept (ground state) : %f MeV\n", intercept); printf(" e-z intercept (ground state) : %f MeV\n", intercept);
@ -72,11 +72,11 @@ void Transfer(
keyParaOut = fopen (filename.Data(), "w+"); keyParaOut = fopen (filename.Data(), "w+");
printf("=========== save key reaction constants to %s \n", filename.Data()); printf("=========== save key reaction constants to %s \n", filename.Data());
fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetMass_b(), "mass_b"); fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_b(), "mass_b");
fprintf(keyParaOut, "%-15d //%s\n", reaction.GetCharge_b(), "charge_b"); fprintf(keyParaOut, "%-15d //%s\n", reactionConfig.recoilLightZ, "charge_b");
fprintf(keyParaOut, "%-15.8f //%s\n", reaction.GetReactionBeta(), "betaCM"); fprintf(keyParaOut, "%-15.8f //%s\n", transfer.GetReactionBeta(), "betaCM");
fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetCMTotalEnergy(), "Ecm"); fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetCMTotalEnergy(), "Ecm");
fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetMass_B(), "mass_B"); fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_B(), "mass_B");
fprintf(keyParaOut, "%-15.4f //%s\n", slope/betaRect, "alpha=slope/betaRect"); fprintf(keyParaOut, "%-15.4f //%s\n", slope/betaRect, "alpha=slope/betaRect");
fflush(keyParaOut); fflush(keyParaOut);
@ -147,9 +147,9 @@ void Transfer(
printf("%3s | %7s | %5s | %3s | %10s | %5s \n", "", "Ex[MeV]", "Xsec", "SF", "sigma[MeV]", "y0[MeV]"); printf("%3s | %7s | %5s | %3s | %10s | %5s \n", "", "Ex[MeV]", "Xsec", "SF", "sigma[MeV]", "y0[MeV]");
printf("----+---------+------+-----+------------+--------\n"); printf("----+---------+------+-----+------------+--------\n");
for(int i = 0; i < n ; i++){ for(int i = 0; i < n ; i++){
reaction.SetExB(ExKnown[i]); transfer.SetExB(ExKnown[i]);
reaction.CalReactionConstant(); transfer.CalReactionConstant();
kCM.push_back(reaction.GetMomentumbCM()); kCM.push_back(transfer.GetMomentumbCM());
y0.push_back(TMath::Sqrt(mb*mb + kCM[i]*kCM[i])/gamma - mb); y0.push_back(TMath::Sqrt(mb*mb + kCM[i]*kCM[i])/gamma - mb);
if( reactionConfig.isDecay ) { if( reactionConfig.isDecay ) {
TLorentzVector temp(0,0,0,0); TLorentzVector temp(0,0,0,0);
@ -169,9 +169,9 @@ void Transfer(
ExKnown.push_back(0.0); ExKnown.push_back(0.0);
ExStrength.push_back(1.0); ExStrength.push_back(1.0);
ExWidth.push_back(0.0); ExWidth.push_back(0.0);
reaction.SetExB(ExKnown[0]); transfer.SetExB(ExKnown[0]);
reaction.CalReactionConstant(); transfer.CalReactionConstant();
kCM.push_back(reaction.GetMomentumbCM()); kCM.push_back(transfer.GetMomentumbCM());
y0.push_back(TMath::Sqrt(mb*mb + kCM[0]*kCM[0])/gamma - mb); y0.push_back(TMath::Sqrt(mb*mb + kCM[0]*kCM[0])/gamma - mb);
} }
@ -215,7 +215,7 @@ void Transfer(
TMacro reactionData(filename.Data()); TMacro reactionData(filename.Data());
double KEAmean = reactionConfig.beamEnergy; double KEAmean = reactionConfig.beamEnergy;
TString str; TString str;
str.Form("%s @ %.2f MeV/u", reaction.GetReactionName_Latex().Data(), KEAmean); str.Form("%s @ %.2f MeV/u", transfer.GetReactionName_Latex().Data(), KEAmean);
config.SetName(str.Data()); config.SetName(str.Data());
config.Write("reactionConfig"); config.Write("reactionConfig");
detGeoTxt.Write("detGeo"); detGeoTxt.Write("detGeo");
@ -486,7 +486,7 @@ void Transfer(
//==== Set Ex of A //==== Set Ex of A
ExAID = gRandom->Integer(nExA); ExAID = gRandom->Integer(nExA);
ExA = ExAList[ExAID]; ExA = ExAList[ExAID];
reaction.SetExA(ExA); transfer.SetExA(ExA);
//==== Set Ex of B //==== Set Ex of B
if( ExKnown.size() == 1 ) { if( ExKnown.size() == 1 ) {
@ -496,7 +496,7 @@ void Transfer(
ExID = exDist->GetRandom(); ExID = exDist->GetRandom();
Ex = ExKnown[ExID]+ (ExWidth[ExID] == 0 ? 0 : gRandom->Gaus(0, ExWidth[ExID])); Ex = ExKnown[ExID]+ (ExWidth[ExID] == 0 ? 0 : gRandom->Gaus(0, ExWidth[ExID]));
} }
reaction.SetExB(Ex); transfer.SetExB(Ex);
//==== Set incident beam //==== Set incident beam
KEA = reactionConfig.beamEnergy; KEA = reactionConfig.beamEnergy;
@ -514,9 +514,9 @@ void Transfer(
phi = 0.0; phi = 0.0;
//==== for taregt scattering //==== for taregt scattering
reaction.SetIncidentEnergyAngle(KEA, theta, 0.); transfer.SetIncidentEnergyAngle(KEA, theta, 0.);
reaction.CalReactionConstant(); transfer.CalReactionConstant();
TLorentzVector PA = reaction.GetPA(); TLorentzVector PA = transfer.GetPA();
//depth = 0; //depth = 0;
if( isTargetScattering ){ if( isTargetScattering ){
@ -525,9 +525,9 @@ void Transfer(
msA.SetTarget(density, depth); msA.SetTarget(density, depth);
TLorentzVector PAnew = msA.Scattering(PA); TLorentzVector PAnew = msA.Scattering(PA);
KEAnew = msA.GetKE()/reactionConfig.beamA; KEAnew = msA.GetKE()/reactionConfig.beamA;
reaction.SetIncidentEnergyAngle(KEAnew, theta, phi); transfer.SetIncidentEnergyAngle(KEAnew, theta, phi);
reaction.CalReactionConstant(); transfer.CalReactionConstant();
Ecm = reaction.GetCMTotalKE(); Ecm = transfer.GetCMTotalKE();
} }
//==== Calculate thetaCM, phiCM //==== Calculate thetaCM, phiCM
@ -541,9 +541,9 @@ void Transfer(
double phiCM = TMath::TwoPi() * gRandom->Rndm(); double phiCM = TMath::TwoPi() * gRandom->Rndm();
//==== Calculate reaction //==== Calculate reaction
reaction.Event(thetaCM, phiCM); transfer.Event(thetaCM, phiCM);
TLorentzVector Pb = reaction.GetPb(); TLorentzVector Pb = transfer.GetPb();
TLorentzVector PB = reaction.GetPB(); TLorentzVector PB = transfer.GetPB();
//==== Calculate energy loss of scattered and recoil in target //==== Calculate energy loss of scattered and recoil in target
if( isTargetScattering ){ if( isTargetScattering ){
@ -593,7 +593,7 @@ void Transfer(
///printf(" thetaCM : %f \n", thetaCM * TMath::RadToDeg()); ///printf(" thetaCM : %f \n", thetaCM * TMath::RadToDeg());
if( Tb > 0 || TB > 0 ){ if( Tb > 0 || TB > 0 ){
helios.CalArrayHit(Pb, reaction.GetCharge_b()); helios.CalArrayHit(Pb, reactionConfig.recoilLightZ);
helios.CalRecoilHit(PB, new_zB); helios.CalRecoilHit(PB, new_zB);
hit = 2; hit = 2;
while( hit > 1 ){ hit = helios.DetAcceptance(); } /// while hit > 1, goto next loop; while( hit > 1 ){ hit = helios.DetAcceptance(); } /// while hit > 1, goto next loop;
@ -640,17 +640,17 @@ void Transfer(
//other recoil detectors //other recoil detectors
if ( detGeo.recoilPos1 != 0 ){ if ( detGeo.recoilPos1 != 0 ){
xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1); xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1);
yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1); yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1);
rhoRecoil1 = helios.GetRecoilR(detGeo.recoilPos1); rhoRecoil1 = helios.GetRecoilR(detGeo.recoilPos1);
} }
if ( detGeo.recoilPos2 != 0 ){ if ( detGeo.recoilPos2 != 0 ){
xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2); xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2);
yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2); yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2);
rhoRecoil2 = helios.GetRecoilR(detGeo.recoilPos2); rhoRecoil2 = helios.GetRecoilR(detGeo.recoilPos2);
} }
std::pair<double,double> ExThetaCM = reaction.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius()); std::pair<double,double> ExThetaCM = transfer.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius());
ExCal = ExThetaCM.first; ExCal = ExThetaCM.first;
thetaCMCal = ExThetaCM.second; thetaCMCal = ExThetaCM.second;

4
README.md
View File

@ -12,9 +12,9 @@ Analysis
├── SOLARIS.sh // bash script to define some env variable and functions ├── SOLARIS.sh // bash script to define some env variable and functions
├── armory // analysis codes, independent from experiment. ├── Armory // analysis codes, independent from experiment.
├── Cleopatra // Swaper for DWBA code Ptolomey ├── Cleopatra // Swaper for DWBA code Ptolomey and simulation
├── data_raw // should be the symbolic link to the raw data, created by SetUpNewExp ├── data_raw // should be the symbolic link to the raw data, created by SetUpNewExp