[Major] overhaul the Cleopaatra due to the change of reactionConfig and detectorGeo
This commit is contained in:
parent
6fad708ee2
commit
b0b37ce950
4
.gitignore
vendored
4
.gitignore
vendored
|
|
@ -12,6 +12,10 @@ data
|
|||
data_raw
|
||||
root_data
|
||||
|
||||
*.in
|
||||
*.out
|
||||
*.txt
|
||||
|
||||
Cleopatra/ExtractXSec
|
||||
Cleopatra/ExtractXSecFromText
|
||||
Cleopatra/FindThetaCM
|
||||
|
|
|
|||
|
|
@ -39,6 +39,44 @@ struct Recoil {
|
|||
|
||||
};
|
||||
|
||||
struct EnergyLevel{
|
||||
|
||||
float Ex, xsec, SF, sigma;
|
||||
|
||||
EnergyLevel(float Ex, float xsec, float SF, float sigma) {
|
||||
this->Ex = Ex;
|
||||
this->xsec = xsec;
|
||||
this->SF = SF;
|
||||
this->sigma = sigma;
|
||||
}
|
||||
|
||||
void Print(std::string str) const {
|
||||
printf("%11.3f %8.1f %5.1f %5.3f%s", Ex, xsec, SF, sigma, str.c_str() );
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
struct ExcitedEnergies {
|
||||
|
||||
std::vector<EnergyLevel> ExList;
|
||||
|
||||
void Clear(){
|
||||
ExList.clear();
|
||||
}
|
||||
|
||||
void Add(float Ex, float xsec, float SF, float sigma){
|
||||
ExList.push_back( EnergyLevel(Ex, xsec, SF, sigma));
|
||||
}
|
||||
|
||||
void Print() const {
|
||||
printf("Energy[MeV] Rel.Xsec SF sigma\n");
|
||||
for( size_t i = 0; i < ExList.size(); i++){
|
||||
ExList[i].Print("\n");
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
class ReactionConfig {
|
||||
|
||||
|
|
@ -52,8 +90,8 @@ public:
|
|||
float beamEx; ///excitation_energy_of_A[MeV]
|
||||
float beamEnergy; ///MeV/u
|
||||
float beamEnergySigma; ///beam-energy_sigma_in_MeV/u
|
||||
float beamAngle; ///beam-angle_in_mrad
|
||||
float beamAngleSigma; ///beam-emittance_in_mrad
|
||||
float beamTheta; ///beam-angle_in_mrad
|
||||
float beamThetaSigma; ///beam-emittance_in_mrad
|
||||
float beamX; ///x_offset_of_Beam_in_mm
|
||||
float beamY; ///y_offset_of_Beam_in_mm
|
||||
|
||||
|
|
@ -64,6 +102,7 @@ public:
|
|||
std::string beamStoppingPowerFile; ///stopping_power_for_beam
|
||||
|
||||
Recoil recoil[2];
|
||||
ExcitedEnergies exList[2];
|
||||
|
||||
int numEvents; ///number_of_Event_being_generated
|
||||
bool isRedo; ///isReDo
|
||||
|
|
@ -97,8 +136,8 @@ inline void ReactionConfig::SetReactionSimple(int beamA, int beamZ,
|
|||
|
||||
beamEnergy = beamEnergy_AMeV;
|
||||
beamEnergySigma = 0;
|
||||
beamAngle = 0;
|
||||
beamAngleSigma = 0;
|
||||
beamTheta = 0;
|
||||
beamThetaSigma = 0;
|
||||
beamX = 0;
|
||||
beamY = 0;
|
||||
|
||||
|
|
@ -124,6 +163,9 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
|
|||
|
||||
if( macro == NULL ) return false;
|
||||
|
||||
exList[0].Clear();
|
||||
exList[1].Clear();
|
||||
|
||||
int recoilFlag = 0;
|
||||
int recoilLine = 0;
|
||||
|
||||
|
|
@ -138,6 +180,7 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
|
|||
// printf("%d |%s|%d|%d\n", i, str[0].c_str(), recoilFlag, recoilLine);
|
||||
|
||||
if( str[0].find("####") != std::string::npos ) break;
|
||||
if( str[0].find("#---") != std::string::npos ) continue;
|
||||
if( str[0].find("#===") != std::string::npos ) {
|
||||
recoilFlag ++;
|
||||
recoilLine = 0;
|
||||
|
|
@ -151,8 +194,8 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
|
|||
|
||||
if( recoilLine == 3 ) beamEnergy = atof(str[0].c_str());
|
||||
if( recoilLine == 4 ) beamEnergySigma = atof(str[0].c_str());
|
||||
if( recoilLine == 5 ) beamAngle = atof(str[0].c_str());
|
||||
if( recoilLine == 6 ) beamAngleSigma = atof(str[0].c_str());
|
||||
if( recoilLine == 5 ) beamTheta = atof(str[0].c_str());
|
||||
if( recoilLine == 6 ) beamThetaSigma = atof(str[0].c_str());
|
||||
if( recoilLine == 7 ) beamX = atof(str[0].c_str());
|
||||
if( recoilLine == 8 ) beamY = atof(str[0].c_str());
|
||||
|
||||
|
|
@ -179,6 +222,8 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
|
|||
if( recoilLine == 5 ) recoil[ID].decayA = atoi(str[0].c_str());
|
||||
if( recoilLine == 6 ) recoil[ID].decayZ = atoi(str[0].c_str());
|
||||
|
||||
if( recoilLine > 6 && str.size() == 4) exList[ID].Add( atoi(str[0].c_str()), atoi(str[1].c_str()), atoi(str[2].c_str()), atoi(str[3].c_str()));
|
||||
|
||||
}
|
||||
|
||||
recoilLine ++;
|
||||
|
|
@ -202,7 +247,7 @@ inline void ReactionConfig::Print() const{
|
|||
printf("------------------------------ Beam\n");
|
||||
printf(" beam : A = %3d, Z = %2d, Ex = %.2f MeV\n", beamA, beamZ, beamEx);
|
||||
printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", beamEnergy, beamEnergySigma, beamEnergySigma/beamEnergy);
|
||||
printf(" Angle : %.2f +- %.2f mrad\n", beamAngle, beamAngleSigma);
|
||||
printf(" Angle : %.2f +- %.2f mrad\n", beamTheta, beamThetaSigma);
|
||||
printf(" offset : (x,y) = (%.2f, %.2f) mmm \n", beamX, beamY);
|
||||
|
||||
printf("------------------------------ Target\n");
|
||||
|
|
@ -215,7 +260,9 @@ inline void ReactionConfig::Print() const{
|
|||
}
|
||||
|
||||
for( int i = 0; i < 2; i ++ ){
|
||||
printf("------------------------------ Recoil-%d\n", i); recoil[i].Print();
|
||||
printf("------------------------------ Recoil-%d\n", i);
|
||||
recoil[i].Print();
|
||||
exList[i].Print();
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@
|
|||
#include "TVector3.h"
|
||||
|
||||
#include "../Cleopatra/ClassIsotope.h"
|
||||
#include "../Armory/ClassReactionConfig.h"
|
||||
|
||||
//=======================================================
|
||||
//#######################################################
|
||||
|
|
@ -34,6 +35,28 @@ public:
|
|||
TLorentzVector GetDaugther_d() {return Pd;}
|
||||
TLorentzVector GetDaugther_D() {return PD;}
|
||||
|
||||
void SetMotherDaugther(Recoil recoil){
|
||||
|
||||
Isotope Mother(recoil.heavyA, recoil.heavyZ);
|
||||
Isotope Daugther_D(recoil.decayA, recoil.decayZ);
|
||||
Isotope Daugther_d(recoil.heavyA - recoil.decayA, recoil.heavyZ - recoil.decayZ);
|
||||
|
||||
zB = recoil.heavyZ;
|
||||
zD = recoil.decayZ;
|
||||
zd = recoil.heavyZ - recoil.decayZ;
|
||||
|
||||
mB = Mother.Mass;
|
||||
mD = Daugther_D.Mass;
|
||||
md = Daugther_d.Mass;
|
||||
|
||||
double Q = mB - mD - md;
|
||||
|
||||
printf("====== decay mode : %s --> %s + %s, Q = %.3f MeV \n", Mother.Name.c_str(), Daugther_d.Name.c_str(), Daugther_D.Name.c_str(), Q);
|
||||
|
||||
isMotherSet = true;
|
||||
|
||||
}
|
||||
|
||||
void SetMotherDaugther(int AB, int zB, int AD, int zD){
|
||||
Isotope Mother(AB, zB);
|
||||
Isotope Daugther_D(AD, zD);
|
||||
|
|
@ -71,6 +94,9 @@ public:
|
|||
PD = temp;
|
||||
Pd = temp;
|
||||
|
||||
PD.SetUniqueID(zD);
|
||||
Pd.SetUniqueID(zd);
|
||||
|
||||
k = TMath::Sqrt((MB+MD+md)*(MB+MD-md)*(MB-MD+md)*(MB-MD-md))/2./MB;
|
||||
|
||||
//in mother's frame, assume isotropic decay
|
||||
|
|
@ -109,6 +135,7 @@ private:
|
|||
TLorentzVector PB, Pd, PD;
|
||||
|
||||
double mB, mD, md;
|
||||
double zB, zD, zd;
|
||||
double theta;
|
||||
|
||||
TF1 * f1;
|
||||
|
|
@ -128,6 +155,11 @@ Decay::Decay(){
|
|||
mB = TMath::QuietNaN();
|
||||
mD = TMath::QuietNaN();
|
||||
md = TMath::QuietNaN();
|
||||
|
||||
zB = 0;
|
||||
zD = 0;
|
||||
zd = 0;
|
||||
|
||||
theta = TMath::QuietNaN();
|
||||
|
||||
k = TMath::QuietNaN();
|
||||
|
|
|
|||
|
|
@ -133,9 +133,12 @@ public:
|
|||
trajectory GetTrajectory_B() const {return orbitB;}
|
||||
|
||||
DetGeo GetDetectorGeometry() const {return detGeo;}
|
||||
Array GetArrayGeometry() const {return array;}
|
||||
|
||||
TString GetHitMessage() const {return hitMessage;}
|
||||
TString GetAcceptanceMessage() const {return accMessage;}
|
||||
TString GetHitMessage() {return hitMessage;}
|
||||
TString GetAcceptanceMessage() { AcceptanceCodeToMsg(acceptanceCode); return acceptanceMsg;}
|
||||
|
||||
TString AcceptanceCodeToMsg(short code );
|
||||
|
||||
private:
|
||||
|
||||
|
|
@ -152,7 +155,8 @@ private:
|
|||
bool isDetReady;
|
||||
|
||||
TString hitMessage;
|
||||
TString accMessage; //acceptance check
|
||||
TString acceptanceMsg; //acceptance check
|
||||
short acceptanceCode;
|
||||
|
||||
double xOff, yOff; // beam position
|
||||
|
||||
|
|
@ -161,6 +165,8 @@ private:
|
|||
bool isCoincidentWithRecoil;
|
||||
|
||||
const double c = 299.792458; //mm/ns
|
||||
|
||||
|
||||
};
|
||||
|
||||
HELIOS::HELIOS(){
|
||||
|
|
@ -179,7 +185,8 @@ HELIOS::HELIOS(){
|
|||
isDetReady = false;
|
||||
|
||||
hitMessage = "";
|
||||
accMessage = "";
|
||||
acceptanceMsg = "";
|
||||
acceptanceCode = 0;
|
||||
|
||||
overrideDetDistance = false;
|
||||
overrideFirstPos = false;
|
||||
|
|
@ -257,50 +264,79 @@ void HELIOS::PrintGeometry() const{
|
|||
|
||||
}
|
||||
|
||||
TString HELIOS::AcceptanceCodeToMsg(short code ){
|
||||
|
||||
switch(code){
|
||||
case 3 : acceptanceMsg = "try one more loop"; break;
|
||||
case 2 : acceptanceMsg = "hit less than the nearest array. increase loop"; break;
|
||||
case 1 : acceptanceMsg = "GOOD!! hit Array"; break;
|
||||
|
||||
case 0 : acceptanceMsg = "detector geometry incomplete."; break;
|
||||
case -1 : acceptanceMsg = "array at upstream, z is downstream."; break;
|
||||
case -2 : acceptanceMsg = "array at downstream, z is upstream."; break;
|
||||
case -3 : acceptanceMsg = "hit at the XY gap."; break;
|
||||
case -4 : acceptanceMsg = "hit more upstream than the array length"; break;
|
||||
case -5 : acceptanceMsg = "hit more downstream than the array length"; break;
|
||||
case -6 : acceptanceMsg = "hit blocker"; break;
|
||||
case -7 : acceptanceMsg = "hit array Z-gap"; break;
|
||||
|
||||
case -10 : acceptanceMsg = "rho is too big"; break;
|
||||
case -11 : acceptanceMsg = "rho is too small"; break;
|
||||
case -12 : acceptanceMsg = "light recoil blocked by recoil detector"; break;
|
||||
case -13 : acceptanceMsg = "more than 3 loops."; break;
|
||||
case -14 : acceptanceMsg = "heavy recoil does not hit recoil detector"; break;
|
||||
case -15 : acceptanceMsg = "det Row ID == -1"; break;
|
||||
default : acceptanceMsg = "unknown error."; break;
|
||||
}
|
||||
|
||||
return acceptanceMsg;
|
||||
|
||||
}
|
||||
|
||||
int HELIOS::CheckDetAcceptance(){
|
||||
|
||||
//CalArrayHit and CalRecoilHit must be done before.
|
||||
|
||||
if( isDetReady == false ) return 0;
|
||||
if( isDetReady == false ) { acceptanceCode = 0; return acceptanceCode; }
|
||||
|
||||
// -1 ========= when recoil direction is not same side of array
|
||||
if( array.firstPos < 0 && orbitb.z > 0 ) {accMessage = "array at upstream, z is downstream."; return -1;}
|
||||
if( array.firstPos > 0 && orbitb.z < 0 ) {accMessage = "array at downstream, z is upstream."; return -1;}
|
||||
if( array.firstPos < 0 && orbitb.z > 0 ) {acceptanceCode = -1; return acceptanceCode;}
|
||||
|
||||
// -2 ========= when recoil direction is not same side of array
|
||||
if( array.firstPos > 0 && orbitb.z < 0 ) {acceptanceCode = -2; return acceptanceCode;}
|
||||
|
||||
// -11 ======== rho is too small
|
||||
if( 2 * orbitb.rho < array.detPerpDist ) { accMessage = "rho is too small"; return -11;}
|
||||
if( 2 * orbitb.rho < array.detPerpDist ) { acceptanceCode = -11; return acceptanceCode;}
|
||||
|
||||
// -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist)
|
||||
if( orbitb.detRowID == -1 ) {accMessage = "det Row ID == -1"; return -15;}
|
||||
if( orbitb.detRowID == -1 ) {acceptanceCode = -15; return acceptanceCode;}
|
||||
|
||||
// -10 =========== when rho is too big .
|
||||
if( detGeo.bore < 2 * orbitb.rho) { accMessage = "rho is too big"; return -10;}
|
||||
if( detGeo.bore < 2 * orbitb.rho) { acceptanceCode = -10; return acceptanceCode;}
|
||||
|
||||
// -14 ========== check particle-B hit radius on recoil dectector
|
||||
if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) {
|
||||
accMessage = "heavy recoil does not hit recoil detector";
|
||||
return -14;
|
||||
}
|
||||
if( isCoincidentWithRecoil && orbitB.R > detGeo.recoilOuterRadius ) {acceptanceCode = -14; return acceptanceCode;}
|
||||
|
||||
//if( isCoincidentWithRecoil && (orbitB.R > rhoRecoilout || orbitB.R < rhoRecoilin) ) return -14;
|
||||
|
||||
// -12 ========= check is particle-b was blocked by recoil detector
|
||||
rhoHit = GetR(detGeo.recoilPos);
|
||||
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 ) { accMessage = "light recoil blocked by recoil detector"; return -12;}
|
||||
if( orbitb.z > 0 && detGeo.recoilPos > 0 && orbitb.z > detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { acceptanceCode = -12; return acceptanceCode;}
|
||||
if( orbitb.z < 0 && detGeo.recoilPos < 0 && orbitb.z < detGeo.recoilPos && rhoHit < detGeo.recoilOuterRadius ) { acceptanceCode = -12; return acceptanceCode;}
|
||||
|
||||
// -13 ========= not more than 3 loops
|
||||
if( orbitb.loop > 3 ) {accMessage = "more than 3 loops."; return -13;}
|
||||
if( orbitb.loop > 3 ) {acceptanceCode = -13; return acceptanceCode;}
|
||||
|
||||
// -2 ========= calculate the "y"-distance from detector center
|
||||
if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) { accMessage = "hit at the XY gap."; return -2;}
|
||||
// -3 ========= calculate the "y"-distance from detector center
|
||||
if( sqrt(orbitb.R*orbitb.R - array.detPerpDist * array.detPerpDist)> array.detWidth/2 ) { acceptanceCode = -3; return acceptanceCode;}
|
||||
|
||||
// -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 ) { accMessage = "hit more upstream than the array length"; 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, -5 ==== when zPos further the range of whole array, more loop would not save
|
||||
if( array.firstPos < 0 && orbitb.z < array.detPos[0] - array.detLength ) { acceptanceCode = -4; return acceptanceCode;}
|
||||
if( array.firstPos > 0 && orbitb.z > array.detPos[array.nDet-1] + array.detLength ) { acceptanceCode = -5; return acceptanceCode;}
|
||||
|
||||
// -4 ======== Hit on blacker
|
||||
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 ) { accMessage = "hit blocker"; return -4;}
|
||||
// -6 ======== Hit on blacker
|
||||
if( array.blocker != 0 && array.firstPos > 0 && array.detPos[0] - array.blocker < orbitb.z && orbitb.z < array.detPos[0] ) {acceptanceCode = -6; return acceptanceCode;}
|
||||
if( array.blocker != 0 && array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z && orbitb.z < array.detPos[array.nDet-1] + array.blocker ) { acceptanceCode = -6; return acceptanceCode;}
|
||||
|
||||
// 2 ====== when zPos less then the nearest position, more loop may hit
|
||||
int increaseLoopFlag = 0;
|
||||
|
|
@ -311,8 +347,8 @@ int HELIOS::CheckDetAcceptance(){
|
|||
orbitb.effLoop += 1.0;
|
||||
orbitb.loop += 1;
|
||||
orbitb.t = orbitb.t0 * orbitb.effLoop;
|
||||
accMessage = " hit less than the nearest array. increase loop ";
|
||||
return 2;
|
||||
acceptanceCode = 2;
|
||||
return acceptanceCode;
|
||||
}
|
||||
|
||||
// 1 ======= check hit array z- position
|
||||
|
|
@ -321,8 +357,8 @@ int HELIOS::CheckDetAcceptance(){
|
|||
if( array.detPos[i] - array.detLength <= orbitb.z && orbitb.z <= array.detPos[i]) {
|
||||
orbitb.detID = i;
|
||||
detX = ( orbitb.z - (array.detPos[i] + array.detLength/2 ))/ array.detLength * 2 ;// range from -1 , 1
|
||||
accMessage = "hit array";
|
||||
return 1;
|
||||
acceptanceCode = 1;
|
||||
return acceptanceCode;
|
||||
}
|
||||
}
|
||||
}else{
|
||||
|
|
@ -331,21 +367,21 @@ int HELIOS::CheckDetAcceptance(){
|
|||
///printf(" %d | %f < z = %f < %f \n", i, array.detPos[i], orbitb.z, array.detPos[i]+length);
|
||||
orbitb.detID = i;
|
||||
detX = ( orbitb.z - (array.detPos[i] - array.detLength/2 ))/ array.detLength*2 ;// range from -1 , 1
|
||||
accMessage = "hit array";
|
||||
return 1;
|
||||
acceptanceCode = 1;
|
||||
return acceptanceCode;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// -5 ======== check hit array gap
|
||||
// -7 ======== check hit array gap
|
||||
if( array.firstPos < 0 ){
|
||||
for( int i = 0; i < array.nDet-1 ; i++){
|
||||
if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) { accMessage = "hit array Z-gap"; return -5; }//increaseLoopFlag = 3;
|
||||
if( array.detPos[i] < orbitb.z && orbitb.z < array.detPos[i+1] - array.detLength ) { acceptanceCode = -7; return acceptanceCode; }//increaseLoopFlag = 3;
|
||||
}
|
||||
}else{
|
||||
for( int i = 0; i < array.nDet-1 ; i++){
|
||||
if( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) { accMessage = "hit array Z-gap"; return -5; }//increaseLoopFlag = 3;
|
||||
if( array.detPos[i] + array.detLength < orbitb.z && orbitb.z < array.detPos[i+1] ) { acceptanceCode = -7; return acceptanceCode; }//increaseLoopFlag = 3;
|
||||
}
|
||||
}
|
||||
if (increaseLoopFlag == 3 ) {
|
||||
|
|
@ -353,12 +389,12 @@ int HELIOS::CheckDetAcceptance(){
|
|||
orbitb.effLoop += 1.0;
|
||||
orbitb.loop += 1;
|
||||
orbitb.t = orbitb.t0 * orbitb.effLoop;
|
||||
accMessage = " try one more loop. ";
|
||||
return 3;
|
||||
acceptanceCode = 3;
|
||||
return acceptanceCode;
|
||||
}
|
||||
|
||||
accMessage = " unknown reason ";
|
||||
return -20; // for unknown reason
|
||||
acceptanceCode = -20 ;
|
||||
return acceptanceCode; // for unknown reason
|
||||
}
|
||||
|
||||
void HELIOS::CalTrajectoryPara(TLorentzVector P, bool isLightRecoil){
|
||||
|
|
|
|||
|
|
@ -22,18 +22,10 @@
|
|||
class TransferReaction {
|
||||
public:
|
||||
TransferReaction(){Inititization();};
|
||||
TransferReaction(string configFile, unsigned short ID = 0);
|
||||
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 configFile, unsigned short ID = 0){
|
||||
Inititization();
|
||||
SetReactionFromFile(configFile, ID);
|
||||
}
|
||||
int recoilA, int recoilZ, float beamEnergy_AMeV);
|
||||
|
||||
~TransferReaction();
|
||||
|
||||
|
|
@ -43,18 +35,20 @@ public:
|
|||
void SetB(int A, int Z);
|
||||
void SetIncidentEnergyAngle(double KEA, double theta, double phi);
|
||||
|
||||
void SetReactionFromFile(std::string configFile, unsigned short ID = 0);
|
||||
void SetReactionSimple(int beamA, int beamZ,
|
||||
int targetA, int targetZ,
|
||||
int recoilA, int recoilZ, float beamEnergy_AMeV);
|
||||
|
||||
void SetExA(double Ex);
|
||||
void SetExB(double Ex);
|
||||
void SetReactionFromFile(string configFile, unsigned short ID = 0);
|
||||
|
||||
TString GetReactionName();
|
||||
TString GetReactionName_Latex();
|
||||
|
||||
ReactionConfig GetRectionConfig() { return config;}
|
||||
Recoil GetRecoil() { return recoil;}
|
||||
ExcitedEnergies GetExList() { return exList;}
|
||||
|
||||
double GetMass_A() const {return mA + ExA;}
|
||||
double GetMass_a() const {return ma;}
|
||||
|
|
@ -65,27 +59,28 @@ public:
|
|||
double GetQValue() {return mA + ExA + ma - mb - mB - ExB;}
|
||||
double GetMaxExB() {return Etot - mb - mB;}
|
||||
|
||||
TLorentzVector GetPA(){return PA;}
|
||||
TLorentzVector GetPa(){return Pa;}
|
||||
TLorentzVector GetPb(){return Pb;}
|
||||
TLorentzVector GetPB(){return PB;}
|
||||
TLorentzVector GetPA() const {return PA;}
|
||||
TLorentzVector GetPa() const {return Pa;}
|
||||
TLorentzVector GetPb() const {return Pb;}
|
||||
TLorentzVector GetPB() const {return PB;}
|
||||
|
||||
void PrintFourVectors() const;
|
||||
void PrintReaction() const;
|
||||
|
||||
double CalkCM(double ExB); //momentum at CM frame
|
||||
void CalReactionConstant();
|
||||
|
||||
std::pair<double, double> CalExThetaCM(double e, double z, double Bfield, double a);
|
||||
void Event(double thetaCM_rad, double phiCM_rad);
|
||||
|
||||
double GetMomentumbCM() {return p;}
|
||||
double GetReactionBeta() {return beta;}
|
||||
double GetReactionGamma() {return gamma;}
|
||||
double GetCMTotalEnergy() {return Etot;}
|
||||
|
||||
std::pair<double, double> CalExThetaCM(double e, double z, double Bfield, double a);
|
||||
double GetEZSlope(double BField) {return 299.792458 * recoil.lightZ * abs(BField) / TMath::TwoPi() * beta / 1000.;} // MeV/mm
|
||||
|
||||
private:
|
||||
|
||||
ExcitedEnergies exList;
|
||||
Recoil recoil;
|
||||
ReactionConfig config;
|
||||
|
||||
|
|
@ -103,6 +98,7 @@ private:
|
|||
double beta, gamma; //CM boost beta
|
||||
double Etot;
|
||||
double p; // CM frame momentum of b, B
|
||||
double slope; // slope of the
|
||||
|
||||
TLorentzVector PA, Pa, Pb, PB;
|
||||
|
||||
|
|
@ -112,6 +108,20 @@ private:
|
|||
|
||||
};
|
||||
|
||||
TransferReaction::TransferReaction(string configFile, unsigned short ID){
|
||||
Inititization();
|
||||
SetReactionFromFile(configFile, ID);
|
||||
}
|
||||
|
||||
TransferReaction::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);
|
||||
}
|
||||
|
||||
void TransferReaction::Inititization(){
|
||||
|
||||
thetaIN = 0.;
|
||||
|
|
@ -224,7 +234,10 @@ void TransferReaction::SetReactionFromFile(string configFile, unsigned short ID)
|
|||
SetA(config.beamA, config.beamZ);
|
||||
Seta(config.targetA, config.targetZ);
|
||||
|
||||
SetExA(config.beamEx);
|
||||
|
||||
recoil = config.recoil[ID];
|
||||
exList = config.exList[ID];
|
||||
|
||||
Setb(recoil.lightA, recoil.lightZ);
|
||||
SetB(recoil.heavyA, recoil.heavyZ);
|
||||
|
|
@ -261,10 +274,15 @@ TString TransferReaction::format(TString name){
|
|||
}
|
||||
TString TransferReaction::GetReactionName_Latex(){
|
||||
TString rName;
|
||||
rName.Form("%s(%s,%s)%s", format(nameA).Data(), format(namea).Data(), format(nameb).Data(), format(nameB).Data());
|
||||
rName.Form("%s(%s,%s)%s @ %.2f MeV/u", format(nameA).Data(), format(namea).Data(), format(nameb).Data(), format(nameB).Data(), config.beamEnergy);
|
||||
return rName;
|
||||
}
|
||||
|
||||
double TransferReaction::CalkCM(double ExB){
|
||||
if( !isBSet || !isReady) return TMath::QuietNaN();
|
||||
return TMath::Sqrt( (Etot*Etot - TMath::Power(mb + mB + ExB,2)) * (Etot*Etot - TMath::Power(mb - mB - ExB,2)) ) / 2 / Etot;
|
||||
}
|
||||
|
||||
void TransferReaction::CalReactionConstant(){
|
||||
if( !isBSet){
|
||||
recoil.heavyA = config.beamA + config.targetA - recoil.lightA;
|
||||
|
|
@ -278,12 +296,10 @@ void TransferReaction::CalReactionConstant(){
|
|||
beta = k / (mA + ExA + ma + T);
|
||||
gamma = 1 / TMath::Sqrt(1- beta * beta);
|
||||
Etot = TMath::Sqrt(TMath::Power(mA + ExA + ma + T,2) - k * k);
|
||||
p = TMath::Sqrt( (Etot*Etot - TMath::Power(mb + mB + ExB,2)) * (Etot*Etot - TMath::Power(mb - mB - ExB,2)) ) / 2 / Etot;
|
||||
|
||||
PA.SetXYZM(0, 0, k, mA + ExA);
|
||||
PA.RotateY(thetaIN);
|
||||
PA.RotateZ(phiIN);
|
||||
|
||||
Pa.SetXYZM(0,0,0,ma);
|
||||
|
||||
PA.SetUniqueID(config.beamZ);
|
||||
|
|
@ -292,6 +308,8 @@ void TransferReaction::CalReactionConstant(){
|
|||
PB.SetUniqueID(recoil.heavyZ);
|
||||
|
||||
isReady = true;
|
||||
p = CalkCM(ExB);
|
||||
|
||||
}
|
||||
|
||||
void TransferReaction::PrintFourVectors() const {
|
||||
|
|
@ -315,7 +333,7 @@ void TransferReaction::PrintReaction() const {
|
|||
printf("------------------------------ Beam\n");
|
||||
printf(" beam : A = %3d, Z = %2d, Ex = %.2f MeV\n", config.beamA, config.beamZ, config.beamEx);
|
||||
printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", config.beamEnergy, config.beamEnergySigma, config.beamEnergySigma/config.beamEnergy);
|
||||
printf(" Angle : %.2f +- %.2f mrad\n", config.beamAngle, config.beamAngleSigma);
|
||||
printf(" Angle : %.2f +- %.2f mrad\n", config.beamTheta, config.beamThetaSigma);
|
||||
printf(" offset : (x,y) = (%.2f, %.2f) mmm \n", config.beamX, config.beamY);
|
||||
|
||||
printf("------------------------------ Target\n");
|
||||
|
|
@ -325,14 +343,14 @@ void TransferReaction::PrintReaction() const {
|
|||
printf(" light : A = %3d, Z = %2d \n", recoil.lightA, recoil.lightZ);
|
||||
printf(" heavy : A = %3d, Z = %2d \n", recoil.heavyA, recoil.heavyZ);
|
||||
printf("=====================================================\n");
|
||||
exList.Print();
|
||||
printf("=====================================================\n");
|
||||
|
||||
}
|
||||
|
||||
void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){
|
||||
|
||||
if( isReady == false ){
|
||||
CalReactionConstant();
|
||||
}
|
||||
if( !isReady ) CalReactionConstant();
|
||||
|
||||
//---- to CM frame
|
||||
TLorentzVector Pc = PA + Pa;
|
||||
|
|
|
|||
BIN
Cleopatra/Cleopatra
Executable file
BIN
Cleopatra/Cleopatra
Executable file
Binary file not shown.
|
|
@ -78,11 +78,10 @@ int main (int argc, char *argv[]) { //TODO add angle range
|
|||
InFileCreator( readFile, ptolemyInFileName, angMin, angMax, angStep);
|
||||
|
||||
//================= run ptolemy
|
||||
|
||||
char command[200];
|
||||
string ptolemyOutFileName = argv[1];
|
||||
ptolemyOutFileName += ".out";
|
||||
sprintf(command, "./ptolemy <%s> %s", ptolemyInFileName.c_str(), ptolemyOutFileName.c_str());
|
||||
sprintf(command, "../Cleopatra/ptolemy <%s> %s", ptolemyInFileName.c_str(), ptolemyOutFileName.c_str());
|
||||
printf("=================== Run Ptolemy\n");
|
||||
printf("%s \n", command);
|
||||
system(command);
|
||||
|
|
|
|||
|
|
@ -23,6 +23,7 @@
|
|||
#include <TString.h>
|
||||
#include <TMath.h>
|
||||
#include <TGraph.h>
|
||||
#include <TMacro.h>
|
||||
#include <TF1.h>
|
||||
#include <TObjArray.h>
|
||||
#include "../Armory/AnalysisLib.h"
|
||||
|
|
@ -293,26 +294,9 @@ int ExtractXSec (string readFile, int indexForElastic=1) {
|
|||
}
|
||||
printf("---------------------------------------------------\n");
|
||||
|
||||
//================================== save *.Ex.txt
|
||||
string saveExName = readFile;
|
||||
int len = saveExName.length();
|
||||
saveExName = saveExName.substr(0, len - 4);
|
||||
saveExName += ".Ex.txt";
|
||||
printf("Output : %s \n", saveExName.c_str());
|
||||
FILE * file_Ex;
|
||||
file_Ex = fopen(saveExName.c_str(), "w+");
|
||||
|
||||
fprintf(file_Ex, "//generated_by_ExtractXSec.h____Ex____Xsec(4pi)____SF____sigma\n");
|
||||
|
||||
for( int i = 0; i < numCal ; i++){
|
||||
fprintf(file_Ex, "%9.5f %9.5f 1.0 0.000\n", Ex[i], partialXsec[i]);
|
||||
}
|
||||
fprintf(file_Ex, "#=====End_of_File\n");
|
||||
fclose(file_Ex);
|
||||
|
||||
//================================== save file.Xsec.txt
|
||||
string saveFileName = readFile;
|
||||
len = saveFileName.length();
|
||||
int len = saveFileName.length();
|
||||
saveFileName = saveFileName.substr(0, len - 4);
|
||||
saveFileName += ".Xsec.txt";
|
||||
printf("Output : %s \n", saveFileName.c_str());
|
||||
|
|
@ -324,7 +308,7 @@ int ExtractXSec (string readFile, int indexForElastic=1) {
|
|||
}
|
||||
|
||||
int space = 19;
|
||||
fprintf(file_out, "%8s\t", "Angel");
|
||||
fprintf(file_out, "%8s\t", "Angle");
|
||||
for( int i = 0; i < numCal ; i++){
|
||||
fprintf(file_out, "%*s", space, title[i].c_str());
|
||||
}
|
||||
|
|
@ -339,14 +323,22 @@ int ExtractXSec (string readFile, int indexForElastic=1) {
|
|||
}
|
||||
fclose(file_out);
|
||||
|
||||
//================================== Make TMacro for ExList
|
||||
|
||||
TMacro ExList;
|
||||
ExList.AddLine("#---Ex relative_xsec SF sigma_in_MeV");
|
||||
for( int i = 0; i < numCal ; i++){
|
||||
ExList.AddLine(Form("%9.5f %9.5f 1.0 0.000", Ex[i], partialXsec[i]));
|
||||
}
|
||||
|
||||
//================================== Save in ROOT
|
||||
len = saveFileName.length();
|
||||
saveFileName = saveFileName.substr(0, len - 9);
|
||||
TString fileName = saveFileName;
|
||||
fileName += ".root";
|
||||
printf("Output : %s \n", fileName.Data());
|
||||
TFile * fileOut = new TFile(fileName, "RECREATE" );
|
||||
|
||||
TFile * fileOut = new TFile(fileName, "RECREATE" );
|
||||
gList = new TObjArray(); ///no SetTitle() method for TObjArray
|
||||
gList->SetName("TGraph of d.s.c");
|
||||
TObjArray * fList = new TObjArray();
|
||||
|
|
@ -372,12 +364,11 @@ int ExtractXSec (string readFile, int indexForElastic=1) {
|
|||
|
||||
fList->Add(dist[i]);
|
||||
|
||||
//delete tempFunc;
|
||||
|
||||
}
|
||||
gList->Write("qList", 1);
|
||||
fList->Write("pList", 1);
|
||||
gList->Write("thetaCM_TGraph", 1);
|
||||
fList->Write("thetaCM_TF1", 1);
|
||||
|
||||
ExList.Write("ExList");
|
||||
|
||||
fileOut->Write();
|
||||
fileOut->Close();
|
||||
|
|
|
|||
|
|
@ -21,7 +21,7 @@
|
|||
|
||||
void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95,
|
||||
std::string basicConfig="reactionConfig.txt",
|
||||
std::string detGeoFileName = "detectorGeo.txt"){
|
||||
std::string detGeoFileName = "detectorGeo.txt", unsigned short ID = 0){
|
||||
|
||||
//---- reaction
|
||||
int AA, zA; //beam
|
||||
|
|
@ -35,24 +35,24 @@ void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95,
|
|||
double xBeam, yBeam; // mm
|
||||
|
||||
/**///========================================================= load files
|
||||
ReactionConfig reactionConfig;
|
||||
ReactionConfig reConfig;
|
||||
DetGeo detGeo;
|
||||
if( reactionConfig.LoadReactionConfig(basicConfig) ){
|
||||
if( reConfig.LoadReactionConfig(basicConfig) ){
|
||||
|
||||
KEAmean = reactionConfig.beamEnergy;
|
||||
KEAsigma = reactionConfig.beamEnergySigma;
|
||||
KEAmean = reConfig.beamEnergy;
|
||||
KEAsigma = reConfig.beamEnergySigma;
|
||||
|
||||
thetaMean = reactionConfig.beamAngle;
|
||||
thetaSigma = reactionConfig.beamAngleSigma;
|
||||
thetaMean = reConfig.beamTheta;
|
||||
thetaSigma = reConfig.beamThetaSigma;
|
||||
|
||||
xBeam = reactionConfig.beamX;
|
||||
yBeam = reactionConfig.beamY;
|
||||
xBeam = reConfig.beamX;
|
||||
yBeam = reConfig.beamY;
|
||||
|
||||
AA = reactionConfig.beamA; zA = reactionConfig.beamZ;
|
||||
Aa = reactionConfig.targetA; za = reactionConfig.targetZ;
|
||||
Ab = reactionConfig.recoilLightA; zb = reactionConfig.recoilLightZ;
|
||||
AA = reConfig.beamA; zA = reConfig.beamZ;
|
||||
Aa = reConfig.targetA; za = reConfig.targetZ;
|
||||
Ab = reConfig.recoil[ID].lightA; zb = reConfig.recoil[ID].lightZ;
|
||||
|
||||
ExA = reactionConfig.beamEx[0];
|
||||
ExA = reConfig.beamEx;
|
||||
|
||||
}else{
|
||||
printf("cannot load %s \n", basicConfig.c_str());
|
||||
|
|
@ -97,12 +97,12 @@ void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95,
|
|||
|
||||
printf("----- loading detector geometery : %s.", detGeoFileName.c_str());
|
||||
if(detGeo.LoadDetectorGeo(detGeoFileName) ){
|
||||
pos = detGeo.array1.detPos;
|
||||
a = detGeo.array1.detPerpDist;
|
||||
length = detGeo.array1.detLength;
|
||||
firstPos = detGeo.array1.firstPos;
|
||||
iDet = detGeo.array1.nDet;
|
||||
jDet = detGeo.array1.mDet;
|
||||
pos = detGeo.array[ID].detPos;
|
||||
a = detGeo.array[ID].detPerpDist;
|
||||
length = detGeo.array[ID].detLength;
|
||||
firstPos = detGeo.array[ID].firstPos;
|
||||
iDet = detGeo.array[ID].nDet;
|
||||
jDet = detGeo.array[ID].mDet;
|
||||
BField = detGeo.Bfield;
|
||||
|
||||
printf("... done.\n");
|
||||
|
|
|
|||
|
|
@ -1,31 +0,0 @@
|
|||
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include "Check_Simulation.C"
|
||||
|
||||
using namespace std;
|
||||
|
||||
int main (int argc, char *argv[]) {
|
||||
|
||||
printf("=================================================================\n");
|
||||
printf("=================== Plot Simulation Canvas ======================\n");
|
||||
printf("=================================================================\n");
|
||||
|
||||
if(argc < 2 ) {
|
||||
printf("Usage: ./PlotSimulation input_root_file [config]\n");
|
||||
exit(0);
|
||||
}else{
|
||||
printf("ROOT file : %s \n", argv[1]);
|
||||
}
|
||||
|
||||
string rootFile = argv[1];
|
||||
string config = "../Armory/Check_Simulation_Config.txt";
|
||||
if( argc >= 3 ) config = argv[2];
|
||||
|
||||
printf("Config File : %s \n", config.c_str());
|
||||
|
||||
Int_t padSize = 500;
|
||||
|
||||
Check_Simulation(rootFile, config, padSize, true);
|
||||
|
||||
}
|
||||
|
|
@ -30,7 +30,7 @@ void PlotTGraphTObjArray(TString rootFileName, bool isSavePNG = false){
|
|||
|
||||
TFile * file = new TFile(rootFileName, "READ");
|
||||
|
||||
TObjArray * gList = (TObjArray *) file->FindObjectAny("qList");
|
||||
TObjArray * gList = (TObjArray *) file->FindObjectAny("thetaCM_TGraph");
|
||||
|
||||
if( gList == NULL ) {
|
||||
printf("No Result was found.\n");
|
||||
|
|
|
|||
|
|
@ -24,16 +24,15 @@ int main (int argc, char *argv[]) {
|
|||
printf("========== Simulate Transfer reaction in HELIOS ==========\n");
|
||||
printf("=================================================================\n");
|
||||
|
||||
if(argc == 2 || argc > 8) {
|
||||
printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6] [7]\n");
|
||||
if(argc == 2 || argc > 7) {
|
||||
printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n");
|
||||
printf(" default file name \n");
|
||||
printf(" [1] reactionConfig.txt (input) reaction Setting \n");
|
||||
printf(" [2] detectorGeo.txt (input) detector Setting \n");
|
||||
printf(" [3] Ex.txt (input) Excitation energies \n");
|
||||
printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n");
|
||||
printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n");
|
||||
printf(" [5] transfer.root (output) rootFile name for output \n");
|
||||
printf(" [6] reaction.dat (output) Key reaction parameters \n");
|
||||
printf(" [7] plot (input) will it plot stuffs [1/0] \n");
|
||||
printf(" [6] plot (input) will it plot stuffs [1/0] \n");
|
||||
|
||||
printf("------------------------------------------------------\n");
|
||||
return 0 ;
|
||||
|
|
@ -41,21 +40,19 @@ int main (int argc, char *argv[]) {
|
|||
|
||||
string basicConfig = "reactionConfig.txt";
|
||||
string heliosDetGeoFile = "detectorGeo.txt";
|
||||
string excitationFile = "Ex.txt"; //when no file, only ground state
|
||||
int ID = 0;
|
||||
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
|
||||
TString saveFileName = "transfer.root";
|
||||
TString filename = "reaction.dat"; //when no file, no output
|
||||
bool isPlot = false;
|
||||
|
||||
if( argc >= 2) basicConfig = argv[1];
|
||||
if( argc >= 3) heliosDetGeoFile = argv[2];
|
||||
if( argc >= 4) excitationFile = argv[3];
|
||||
if( argc >= 4) ID = atoi(argv[3]);
|
||||
if( argc >= 5) ptolemyRoot = argv[4];
|
||||
if( argc >= 6) saveFileName = argv[5];
|
||||
if( argc >= 7) filename = argv[6];
|
||||
if( argc >= 8) isPlot = atoi(argv[7]);
|
||||
if( argc >= 7) isPlot = atoi(argv[7]);
|
||||
|
||||
Transfer( basicConfig, heliosDetGeoFile, excitationFile, ptolemyRoot, saveFileName, filename);
|
||||
Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName);
|
||||
|
||||
//run Armory/Check_Simulation
|
||||
if( isPlot ){
|
||||
|
|
|
|||
|
|
@ -23,185 +23,159 @@ double exDistFunc(Double_t *x, Double_t * par){
|
|||
return par[(int) x[0]];
|
||||
}
|
||||
|
||||
void PrintEZPlotPara(TransferReaction tran, HELIOS helios){
|
||||
|
||||
printf("==================================== E-Z plot slope\n");
|
||||
double betaRect = tran.GetReactionBeta() ;
|
||||
double gamma = tran.GetReactionGamma();
|
||||
double mb = tran.GetMass_b();
|
||||
double pCM = tran.GetMomentumbCM();
|
||||
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
|
||||
double slope = tran.GetEZSlope(helios.GetBField()); /// MeV/mm
|
||||
printf(" e-z slope : %f MeV/mm\n", slope);
|
||||
double intercept = q/gamma - mb; // MeV
|
||||
printf(" e-z intercept (ground state) : %f MeV\n", intercept);
|
||||
|
||||
}
|
||||
|
||||
void Transfer(
|
||||
string basicConfig = "reactionConfig.txt",
|
||||
string heliosDetGeoFile = "detectorGeo.txt",
|
||||
string excitationFile = "Ex.txt", ///when no file, only ground state
|
||||
TString ptolemyRoot = "DWBA.root", /// when no file, use isotropic distribution of thetaCM
|
||||
TString saveFileName = "transfer.root",
|
||||
TString filename = "reaction.dat"){ /// when no file, no output.
|
||||
std::string basicConfig = "reactionConfig.txt",
|
||||
std::string heliosDetGeoFile = "detectorGeo.txt",
|
||||
unsigned short ID = 0, // this is the ID for the array
|
||||
TString ptolemyRoot = "DWBA.root",
|
||||
TString saveFileName = "transfer.root"){
|
||||
|
||||
//############################################# Set Reaction
|
||||
|
||||
TransferReaction transfer;
|
||||
transfer.SetReactionFromFile(basicConfig);
|
||||
HELIOS helios;
|
||||
Decay decay;
|
||||
|
||||
std::vector<double> kbCM; /// momentum of b in CM frame
|
||||
TF1 * exDist = nullptr;
|
||||
|
||||
transfer.SetReactionFromFile(basicConfig, ID);
|
||||
helios.SetDetectorGeometry(heliosDetGeoFile, ID);
|
||||
|
||||
printf("*****************************************************************\n");
|
||||
printf("*\e[1m\e[33m %27s \e[0m*\n", transfer.GetReactionName().Data());
|
||||
printf("*****************************************************************\n");
|
||||
printf("----- loading reaction setting from %s. \n", basicConfig.c_str());
|
||||
printf("\e[32m#################################### Beam \e[0m\n");
|
||||
printf("----- loading geometry setting from %s. \n", heliosDetGeoFile.c_str());
|
||||
|
||||
const ReactionConfig reactionConfig = transfer.GetRectionConfig();
|
||||
printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n");
|
||||
|
||||
reactionConfig.PrintReactionConfig();
|
||||
transfer.PrintReaction();
|
||||
|
||||
vector<float> ExAList = reactionConfig.beamEx;
|
||||
int nExA = (int) ExAList.size();
|
||||
|
||||
//############################################# Set HELIOS
|
||||
printf("\e[32m#################################### HELIOS configuration\e[0m\n");
|
||||
HELIOS helios;
|
||||
helios.SetDetectorGeometry(heliosDetGeoFile);
|
||||
|
||||
const DetGeo detGeo = helios.GetDetectorGeometry();
|
||||
|
||||
printf("==================================== E-Z plot slope\n");
|
||||
double betaRect = transfer.GetReactionBeta() ;
|
||||
double gamma = transfer.GetReactionGamma();
|
||||
double mb = transfer.GetMass_b();
|
||||
double pCM = transfer.GetMomentumbCM();
|
||||
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
|
||||
double slope = 299.792458 * reactionConfig.recoilLightZ * abs(helios.GetBField()) / TMath::TwoPi() * betaRect / 1000.; /// MeV/mm
|
||||
printf(" e-z slope : %f MeV/mm\n", slope);
|
||||
double intercept = q/gamma - mb; // MeV
|
||||
printf(" e-z intercept (ground state) : %f MeV\n", intercept);
|
||||
|
||||
//############################################# save reaction.dat
|
||||
if( filename != "" ) {
|
||||
FILE * keyParaOut;
|
||||
keyParaOut = fopen (filename.Data(), "w+");
|
||||
|
||||
printf("=========== save key reaction constants to %s \n", filename.Data());
|
||||
fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_b(), "mass_b");
|
||||
fprintf(keyParaOut, "%-15d //%s\n", reactionConfig.recoilLightZ, "charge_b");
|
||||
fprintf(keyParaOut, "%-15.8f //%s\n", transfer.GetReactionBeta(), "betaCM");
|
||||
fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetCMTotalEnergy(), "Ecm");
|
||||
fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_B(), "mass_B");
|
||||
fprintf(keyParaOut, "%-15.4f //%s\n", slope/betaRect, "alpha=slope/betaRect");
|
||||
|
||||
fflush(keyParaOut);
|
||||
fclose(keyParaOut);
|
||||
if(transfer.GetRecoil().isDecay) {
|
||||
decay.SetMotherDaugther(transfer.GetRecoil());
|
||||
}
|
||||
|
||||
//############################################# Target scattering, only energy loss
|
||||
bool isTargetScattering = reactionConfig.isTargetScattering;
|
||||
float density = reactionConfig.targetDensity;
|
||||
float targetThickness = reactionConfig.targetThickness;
|
||||
helios.PrintGeometry();
|
||||
PrintEZPlotPara(transfer, helios);
|
||||
|
||||
if(isTargetScattering) printf("\e[32m#################################### Target Scattering\e[0m\n");
|
||||
TargetScattering msA;
|
||||
TargetScattering msB;
|
||||
TargetScattering msb;
|
||||
|
||||
if(reactionConfig.isTargetScattering) printf("======== Target : (thickness : %6.2f um) x (density : %6.2f g/cm3) = %6.2f ug/cm2\n",
|
||||
targetThickness * 1e+4,
|
||||
density,
|
||||
targetThickness * density * 1e+6);
|
||||
DetGeo detGeo = helios.GetDetectorGeometry();
|
||||
Array array = helios.GetArrayGeometry();
|
||||
ReactionConfig reactConfig = transfer.GetRectionConfig();
|
||||
Recoil recoil = transfer.GetRecoil();
|
||||
|
||||
if( reactionConfig.isTargetScattering ){
|
||||
msA.LoadStoppingPower(reactionConfig.beamStoppingPowerFile);
|
||||
msb.LoadStoppingPower(reactionConfig.recoilLightStoppingPowerFile);
|
||||
msB.LoadStoppingPower(reactionConfig.recoilHeavyStoppingPowerFile);
|
||||
}
|
||||
//*############################################# save reaction.dat
|
||||
// if( filename != "" ) {
|
||||
// FILE * keyParaOut;
|
||||
// keyParaOut = fopen (filename.Data(), "w+");
|
||||
|
||||
//############################################# Decay of particle-B
|
||||
Decay decay;
|
||||
if(reactionConfig.isDecay) {
|
||||
printf("\e[32m#################################### Decay\e[0m\n");
|
||||
decay.SetMotherDaugther(reactionConfig.recoilHeavyA,
|
||||
reactionConfig.recoilHeavyZ,
|
||||
reactionConfig.heavyDecayA,
|
||||
reactionConfig.heavyDecayZ);
|
||||
}
|
||||
//############################################# loading excitation energy
|
||||
printf("\e[32m#################################### excitation energies\e[0m\n");
|
||||
vector<double> ExKnown;
|
||||
vector<double> ExStrength;
|
||||
vector<double> ExWidth;
|
||||
vector<double> SF;
|
||||
vector<double> y0; /// intercept of e-z plot
|
||||
vector<double> kCM; /// momentum of b in CM frame
|
||||
printf("----- loading excitation energy levels (%s).", excitationFile.c_str());
|
||||
ifstream file;
|
||||
file.open(excitationFile.c_str());
|
||||
string isotopeName;
|
||||
if( file.is_open() ){
|
||||
string line;
|
||||
while( getline(file, line) ){
|
||||
///printf("%s \n", line.c_str());
|
||||
if( line.substr(0,2) == "//" ) continue;
|
||||
if( line.substr(0,2) == "#=" ) break;
|
||||
// printf("=========== save key reaction constants to %s \n", filename.Data());
|
||||
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_b(), "mass_b");
|
||||
// fprintf(keyParaOut, "%-15d //%s\n", reactConfig.recoilLightZ, "charge_b");
|
||||
// fprintf(keyParaOut, "%-15.8f //%s\n", transfer.GetReactionBeta(), "betaCM");
|
||||
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetCMTotalEnergy(), "Ecm");
|
||||
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetMass_B(), "mass_B");
|
||||
// fprintf(keyParaOut, "%-15.4f //%s\n", slope/betaRect, "alpha=slope/betaRect");
|
||||
|
||||
vector<string> str = AnalysisLib::SplitStr(line, " ");
|
||||
// fflush(keyParaOut);
|
||||
// fclose(keyParaOut);
|
||||
// }
|
||||
|
||||
ExKnown.push_back(atof(str[0].c_str()));
|
||||
ExStrength.push_back(atof(str[1].c_str()));
|
||||
SF.push_back(atof(str[2].c_str()));
|
||||
ExWidth.push_back(atof(str[3].c_str()));
|
||||
//*############################################# Target scattering, only energy loss
|
||||
// bool isTargetScattering = reactConfig.isTargetScattering;
|
||||
// float density = reactConfig.targetDensity;
|
||||
// float targetThickness = reactConfig.targetThickness;
|
||||
|
||||
}
|
||||
file.close();
|
||||
printf("... done.\n");
|
||||
int n = (int) ExKnown.size();
|
||||
// if(isTargetScattering) printf("\e[32m#################################### Target Scattering\e[0m\n");
|
||||
// TargetScattering msA;
|
||||
// TargetScattering msB;
|
||||
// TargetScattering msb;
|
||||
|
||||
printf("%3s | %7s | %5s | %3s | %10s | %5s \n", "", "Ex[MeV]", "Xsec", "SF", "sigma[MeV]", "y0[MeV]");
|
||||
printf("----+---------+------+-----+------------+--------\n");
|
||||
for(int i = 0; i < n ; i++){
|
||||
transfer.SetExB(ExKnown[i]);
|
||||
transfer.CalReactionConstant();
|
||||
kCM.push_back(transfer.GetMomentumbCM());
|
||||
y0.push_back(TMath::Sqrt(mb*mb + kCM[i]*kCM[i])/gamma - mb);
|
||||
if( reactionConfig.isDecay ) {
|
||||
TLorentzVector temp(0,0,0,0);
|
||||
int decayID = decay.CalDecay(temp, ExKnown[i], 0);
|
||||
if( decayID == 1) {
|
||||
printf("%3d | %7.2f | %5.2f | %3.1f | %5.3f | %5.2f --> Decay. \n", i, ExKnown[i], ExStrength[i], SF[i], ExWidth[i], y0[i]);
|
||||
}else{
|
||||
printf("%3d | %7.2f | %5.2f | %3.1f | %5.3f | %5.2f \n", i, ExKnown[i], ExStrength[i], SF[i], ExWidth[i], y0[i]);
|
||||
}
|
||||
}else{
|
||||
printf("%3d | %7.2f | %5.2f | %3.1f | %5.3f | %5.2f \n", i, ExKnown[i], ExStrength[i], SF[i], ExWidth[i], y0[i]);
|
||||
}
|
||||
}
|
||||
printf("----+---------+-------+-----+------------+--------\n");
|
||||
}else{
|
||||
printf("... fail ------> only ground state.\n");
|
||||
ExKnown.push_back(0.0);
|
||||
ExStrength.push_back(1.0);
|
||||
ExWidth.push_back(0.0);
|
||||
transfer.SetExB(ExKnown[0]);
|
||||
transfer.CalReactionConstant();
|
||||
kCM.push_back(transfer.GetMomentumbCM());
|
||||
y0.push_back(TMath::Sqrt(mb*mb + kCM[0]*kCM[0])/gamma - mb);
|
||||
}
|
||||
// if(reactConfig.isTargetScattering) printf("======== Target : (thickness : %6.2f um) x (density : %6.2f g/cm3) = %6.2f ug/cm2\n",
|
||||
// targetThickness * 1e+4,
|
||||
// density,
|
||||
// targetThickness * density * 1e+6);
|
||||
|
||||
//---- create Ex-distribution
|
||||
TF1 * exDist = NULL;
|
||||
if( ExKnown.size() > 1 ) {
|
||||
printf("---- creating Ex-distribution \n");
|
||||
int exSize = ExKnown.size();
|
||||
exDist = new TF1("exDist", exDistFunc, 0, exSize, exSize);
|
||||
for(int i = 0; i < exSize; i++){
|
||||
exDist->SetParameter(i, ExStrength[i]*SF[i]);
|
||||
}
|
||||
}
|
||||
// if( reactConfig.isTargetScattering ){
|
||||
// msA.LoadStoppingPower(reactConfig.beamStoppingPowerFile);
|
||||
// msb.LoadStoppingPower(reactConfig.recoilLightStoppingPowerFile);
|
||||
// msB.LoadStoppingPower(reactConfig.recoilHeavyStoppingPowerFile);
|
||||
// }
|
||||
|
||||
//*############################################# Decay of particle-B
|
||||
// Decay decay[2];
|
||||
// if(reactConfig.isDecay) {
|
||||
// printf("\e[32m#################################### Decay\e[0m\n");
|
||||
// decay.SetMotherDaugther(reactConfig.recoilHeavyA,
|
||||
// reactConfig.recoilHeavyZ,
|
||||
// reactConfig.heavyDecayA,
|
||||
// reactConfig.heavyDecayZ);
|
||||
// }
|
||||
|
||||
ExcitedEnergies exList = transfer.GetRectionConfig().exList[ID];
|
||||
|
||||
//############################################# Load DWBAroot for thetaCM distribution
|
||||
printf("\e[32m#################################### Load DWBA input : %s \e[0m\n", ptolemyRoot.Data());
|
||||
TF1 * dist = NULL;
|
||||
TFile * distFile = new TFile(ptolemyRoot, "read");
|
||||
TObjArray * distList = NULL;
|
||||
TObjArray * distList = nullptr;
|
||||
TMacro * dwbaExList = nullptr;
|
||||
if( distFile->IsOpen() ) {
|
||||
distList = (TObjArray *) distFile->FindObjectAny("pList"); // the function List
|
||||
int distSize = distList->GetLast() + 1;
|
||||
if( distSize != ExKnown.size() ) {
|
||||
printf(" The number of distribution from Ptolmey Calculation is not equal to number of Ex input \n");
|
||||
printf(" --> the Ptolmey calculation is probably not matched with Ex input.\n");
|
||||
printf(" .... not use DWBA input. \n");
|
||||
distFile->Close();
|
||||
printf("--------- Found DWBA thetaCM distributions. Use the ExList from DWBA.\n");
|
||||
|
||||
distList = (TObjArray *) distFile->FindObjectAny("thetaCM_TF1"); // the function List
|
||||
|
||||
exList.Clear();
|
||||
|
||||
dwbaExList = (TMacro *) distFile->FindObjectAny("ExList");
|
||||
int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ;
|
||||
for(int i = 1; i <= numEx ; i++){
|
||||
string temp = dwbaExList->GetListOfLines()->At(i)->GetName();
|
||||
if( temp[0] == '/' ) continue;
|
||||
vector<string> tempStr = AnalysisLib::SplitStr(temp, " ");
|
||||
exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
|
||||
}
|
||||
|
||||
}else{
|
||||
printf("------- no DWBA input. \n");
|
||||
printf("------- no DWBA input. Use the ExList from %s\n", basicConfig.c_str());
|
||||
}
|
||||
|
||||
|
||||
printf("------------------------------ Heavy Recoil excitation\n");
|
||||
printf("Energy[MeV] Rel.Xsec SF sigma\n");
|
||||
|
||||
int numEx = exList.ExList.size();
|
||||
|
||||
for( int j = 0; j < numEx; j++){
|
||||
double ex = exList.ExList[j].Ex;
|
||||
kbCM.push_back(transfer.CalkCM(ex));
|
||||
int decayID = decay.CalDecay(TLorentzVector (0,0,0,0), ex, 0);
|
||||
exList.ExList[j].Print(decayID == 1 ? "-->Decay" : "\n");
|
||||
}
|
||||
|
||||
//---- create Ex-distribution
|
||||
if( exList.ExList.size() > 1 ) {
|
||||
printf("---- creating Ex-distribution \n");
|
||||
exDist = new TF1("exDist", exDistFunc, 0, numEx, numEx);
|
||||
for(int q = 0; q < numEx; q++){
|
||||
exDist->SetParameter(q, exList.ExList[q].xsec*exList.ExList[q].SF);
|
||||
}
|
||||
}
|
||||
|
||||
//############################################# build tree
|
||||
|
|
@ -211,36 +185,21 @@ void Transfer(
|
|||
|
||||
TMacro config(basicConfig.c_str());
|
||||
TMacro detGeoTxt(heliosDetGeoFile.c_str());
|
||||
TMacro exList(excitationFile.c_str());
|
||||
TMacro reactionData(filename.Data());
|
||||
double KEAmean = reactionConfig.beamEnergy;
|
||||
TString str;
|
||||
str.Form("%s @ %.2f MeV/u", transfer.GetReactionName_Latex().Data(), KEAmean);
|
||||
config.SetName(str.Data());
|
||||
config.SetName(transfer.GetReactionName_Latex().Data());
|
||||
config.Write("reactionConfig");
|
||||
detGeoTxt.Write("detGeo");
|
||||
exList.Write("ExList");
|
||||
reactionData.Write("reactionData");
|
||||
|
||||
if( distList != NULL ) distList->Write("DWBA", 1);
|
||||
if( dwbaExList != NULL ) dwbaExList->Write("DWBA_ExList", 1);
|
||||
|
||||
|
||||
TMacro hitMeaning;
|
||||
str = "=======================meaning of Hit ID\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " 1 = light recoil hit array & heavy recoil hit recoil\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " 0 = no detector\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -1 = light recoil go opposite side of array\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -2 = light recoil hit > det width\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -3 = light recoil hit > array \n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -4 = light recoil hit blocker \n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -10 = light recoil orbit radius too big \n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -11 = light recoil orbit radius too small\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -12 = when reocol at the same side of array, light recoil blocked by recoil detector\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -13 = more than 3 loops\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -14 = heavy recoil did not hit recoil \n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -15 = cannot find hit on array\n"; hitMeaning.AddLine(str.Data());
|
||||
str = " -20 = unknown\n"; hitMeaning.AddLine(str.Data());
|
||||
str = "===========================================\n"; hitMeaning.AddLine(str.Data());
|
||||
|
||||
hitMeaning.AddLine("======================= meaning of Hit\n");
|
||||
for( int code = -15 ; code <= 1; code ++ ){
|
||||
hitMeaning.AddLine( Form( "%4d = %s", code, helios.AcceptanceCodeToMsg(code).Data() ));
|
||||
}
|
||||
hitMeaning.AddLine(" other = unknown\n");
|
||||
hitMeaning.AddLine("===========================================\n");
|
||||
hitMeaning.Write("hitMeaning");
|
||||
|
||||
int hit; /// the output of Helios.CalHit
|
||||
|
|
@ -275,11 +234,6 @@ void Transfer(
|
|||
tree->Branch("rho", &rho, "orbit_radius_light/D");
|
||||
tree->Branch("rhoB", &rhoB, "orbit_radius_heavy/D");
|
||||
|
||||
int ExAID;
|
||||
double ExA;
|
||||
tree->Branch("ExAID", &ExAID, "ExAID/I");
|
||||
tree->Branch("ExA", &ExA, "ExA/D");
|
||||
|
||||
int ExID;
|
||||
double Ex;
|
||||
tree->Branch("ExID", &ExID, "ExID/I");
|
||||
|
|
@ -289,25 +243,20 @@ void Transfer(
|
|||
tree->Branch("ExCal", &ExCal, "ExCal/D");
|
||||
tree->Branch("thetaCMCal", &thetaCMCal, "thetaCMCal/D");
|
||||
|
||||
double KEA, theta, phi;
|
||||
tree->Branch("beamTheta", &theta, "beamTheta/D");
|
||||
tree->Branch("beamPhi", &phi, "beamPhi/D");
|
||||
tree->Branch("beamKEA", &KEA, "beamKEA/D");
|
||||
|
||||
double TbLoss; /// energy loss of particle-b from target scattering
|
||||
double KEAnew; ///beam energy after target scattering
|
||||
double depth; /// reaction depth;
|
||||
double Ecm;
|
||||
if( reactionConfig.isTargetScattering ){
|
||||
tree->Branch("depth", &depth, "depth/D");
|
||||
tree->Branch("TbLoss", &TbLoss, "TbLoss/D");
|
||||
tree->Branch("KEAnew", &KEAnew, "KEAnew/D");
|
||||
tree->Branch("Ecm", &Ecm, "Ecm/D");
|
||||
}
|
||||
// double TbLoss; /// energy loss of particle-b from target scattering
|
||||
// double KEAnew; ///beam energy after target scattering
|
||||
// double depth; /// reaction depth;
|
||||
// double Ecm;
|
||||
// if( reactConfig.isTargetScattering ){
|
||||
// tree->Branch("depth", &depth, "depth/D");
|
||||
// tree->Branch("TbLoss", &TbLoss, "TbLoss/D");
|
||||
// tree->Branch("KEAnew", &KEAnew, "KEAnew/D");
|
||||
// tree->Branch("Ecm", &Ecm, "Ecm/D");
|
||||
// }
|
||||
|
||||
double decayTheta; /// the change of thetaB due to decay
|
||||
double xRecoil_d, yRecoil_d, rhoRecoil_d, Td;
|
||||
if( reactionConfig.isDecay ) {
|
||||
if( recoil.isDecay ) {
|
||||
tree->Branch("decayTheta", &decayTheta, "decayTheta/D");
|
||||
tree->Branch("xRecoil_d", &xRecoil_d, "xRecoil_d/D");
|
||||
tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D");
|
||||
|
|
@ -325,6 +274,7 @@ void Transfer(
|
|||
tree->Branch("yRecoil", &yRecoil, "yRecoil/D");
|
||||
tree->Branch("rhoRecoil", &rhoRecoil, "rhoRecoil/D");
|
||||
|
||||
|
||||
///in case need ELUM
|
||||
double xElum1, yElum1, rhoElum1;
|
||||
if( detGeo.elumPos1 != 0 ) {
|
||||
|
|
@ -360,6 +310,12 @@ void Transfer(
|
|||
const int gxSize = 50;
|
||||
TF1 ** gx = new TF1*[gxSize];
|
||||
TString name;
|
||||
|
||||
double mb = transfer.GetMass_b();
|
||||
double betaRect = transfer.GetReactionBeta();
|
||||
double gamma = transfer.GetReactionGamma();
|
||||
double slope = transfer.GetEZSlope(helios.GetBField()); /// MeV/mm
|
||||
|
||||
for( int i = 0; i < gxSize; i++){
|
||||
name.Form("g%d", i);
|
||||
gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000);
|
||||
|
|
@ -376,40 +332,25 @@ void Transfer(
|
|||
printf("/");
|
||||
if( i > 1 && i % 40 == 0 ) printf("\n");
|
||||
}
|
||||
gList->Write("gList", TObject::kSingleKey);
|
||||
gList->Write("EZ_thetaCM", TObject::kSingleKey);
|
||||
printf(" %d constant thetaCM functions\n", gxSize);
|
||||
|
||||
int n = ExKnown.size();
|
||||
TObjArray * fList = new TObjArray();
|
||||
TF1** f = new TF1*[n];
|
||||
for( int i = 0; i< n ; i++){
|
||||
name.Form("f%d", i);
|
||||
f[i] = new TF1(name, "[0] + [1] * x", -1000, 1000);
|
||||
f[i]->SetParameter(0, y0[i]);
|
||||
f[i]->SetParameter(1, slope);
|
||||
f[i]->SetNpx(1000);
|
||||
fList->Add(f[i]);
|
||||
printf(".");
|
||||
}
|
||||
fList->Write("fList", TObject::kSingleKey);
|
||||
printf(" %d e-z infinte-small detector functions\n", n);
|
||||
|
||||
//--- cal modified f
|
||||
TObjArray * fxList = new TObjArray();
|
||||
TGraph ** fx = new TGraph*[n];
|
||||
TGraph ** fx = new TGraph*[numEx];
|
||||
vector<double> px, py;
|
||||
int countfx = 0;
|
||||
for( int j = 0 ; j < n; j++){
|
||||
for( int j = 0 ; j < numEx; j++){
|
||||
double a = helios.GetDetRadius();
|
||||
double q = TMath::Sqrt(mb*mb + kCM[j] * kCM[j] );
|
||||
double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[j] );
|
||||
px.clear();
|
||||
py.clear();
|
||||
countfx = 0;
|
||||
for(int i = 0; i < 100; i++){
|
||||
double thetacm = TMath::Pi()/TMath::Log(100) * (TMath::Log(100) - TMath::Log(100-i)) ;//using log scale, for more point in small angle.
|
||||
double temp = TMath::TwoPi() * slope / betaRect / kCM[j] * a / TMath::Sin(thetacm);
|
||||
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
|
||||
double pyTemp = gamma * q - mb - gamma * betaRect * kCM[j] * TMath::Cos(thetacm);
|
||||
double temp = TMath::TwoPi() * slope / betaRect / kbCM[j] * a / TMath::Sin(thetacm);
|
||||
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
|
||||
double pyTemp = gamma * q - mb - gamma * betaRect * kbCM[j] * TMath::Cos(thetacm);
|
||||
if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
|
||||
px.push_back(pxTemp);
|
||||
py.push_back(pyTemp);
|
||||
|
|
@ -423,22 +364,22 @@ void Transfer(
|
|||
fxList->Add(fx[j]);
|
||||
printf(",");
|
||||
}
|
||||
fxList->Write("fxList", TObject::kSingleKey);
|
||||
printf(" %d e-z finite-size detector functions\n", n);
|
||||
fxList->Write("EZCurve", TObject::kSingleKey);
|
||||
printf(" %d e-z finite-size detector functions\n", numEx);
|
||||
|
||||
//--- cal modified thetaCM vs z
|
||||
TObjArray * txList = new TObjArray();
|
||||
TGraph ** tx = new TGraph*[n];
|
||||
for( int j = 0 ; j < n; j++){
|
||||
TGraph ** tx = new TGraph*[numEx];
|
||||
for( int j = 0 ; j < numEx; j++){
|
||||
double a = helios.GetDetRadius();
|
||||
double q = TMath::Sqrt(mb*mb + kCM[j] * kCM[j] );
|
||||
double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[j] );
|
||||
px.clear();
|
||||
py.clear();
|
||||
countfx = 0;
|
||||
for(int i = 0; i < 100; i++){
|
||||
double thetacm = (i + 8.) * TMath::DegToRad();
|
||||
double temp = TMath::TwoPi() * slope / betaRect / kCM[j] * a / TMath::Sin(thetacm);
|
||||
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi());
|
||||
double temp = TMath::TwoPi() * slope / betaRect / kbCM[j] * a / TMath::Sin(thetacm);
|
||||
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi());
|
||||
double pyTemp = thetacm * TMath::RadToDeg();
|
||||
if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
|
||||
px.push_back(pxTemp);
|
||||
|
|
@ -453,8 +394,8 @@ void Transfer(
|
|||
txList->Add(tx[j]);
|
||||
printf("*");
|
||||
}
|
||||
txList->Write("txList", TObject::kSingleKey);
|
||||
printf(" %d thetaCM-z for finite-size detector functions\n", n);
|
||||
txList->Write("thetaCM_Z", TObject::kSingleKey);
|
||||
printf(" %d thetaCM-z for finite-size detector functions\n", numEx);
|
||||
|
||||
//========timer
|
||||
TBenchmark clock;
|
||||
|
|
@ -464,7 +405,7 @@ void Transfer(
|
|||
shown = false;
|
||||
|
||||
//change the number of event into human easy-to-read form
|
||||
int numEvent = reactionConfig.numEvents;
|
||||
int numEvent = reactConfig.numEvents;
|
||||
int digitLen = TMath::Floor(TMath::Log10(numEvent));
|
||||
TString numEventStr;
|
||||
if( 3 <= digitLen && digitLen < 6 ){
|
||||
|
|
@ -476,59 +417,51 @@ void Transfer(
|
|||
}
|
||||
printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data());
|
||||
|
||||
double KEA = reactConfig.beamEnergy;
|
||||
double theta = reactConfig.beamTheta;
|
||||
double phi = 0.0;
|
||||
|
||||
//====================================================== calculate event
|
||||
int count = 0;
|
||||
for( int i = 0; i < numEvent; i++){
|
||||
bool redoFlag = true;
|
||||
if( !reactionConfig.isRedo ) redoFlag = false;
|
||||
if( !reactConfig.isRedo ) redoFlag = false;
|
||||
do{
|
||||
|
||||
//==== Set Ex of A
|
||||
ExAID = gRandom->Integer(nExA);
|
||||
ExA = ExAList[ExAID];
|
||||
transfer.SetExA(ExA);
|
||||
|
||||
//==== Set Ex of B
|
||||
if( ExKnown.size() == 1 ) {
|
||||
if( numEx == 1 ) {
|
||||
ExID = 0;
|
||||
Ex = ExKnown[0] + (ExWidth[0] == 0 ? 0 : gRandom->Gaus(0, ExWidth[0]));
|
||||
Ex = exList.ExList[0].Ex + (exList.ExList[0].sigma == 0 ? 0 : gRandom->Gaus(0, exList.ExList[0].sigma));
|
||||
}else{
|
||||
ExID = exDist->GetRandom();
|
||||
Ex = ExKnown[ExID]+ (ExWidth[ExID] == 0 ? 0 : gRandom->Gaus(0, ExWidth[ExID]));
|
||||
Ex = exList.ExList[ExID].Ex + (exList.ExList[ExID].sigma == 0 ? 0 : gRandom->Gaus(0, exList.ExList[ExID].sigma));
|
||||
}
|
||||
transfer.SetExB(Ex);
|
||||
|
||||
//==== Set incident beam
|
||||
KEA = reactionConfig.beamEnergy;
|
||||
if( reactionConfig.beamEnergySigma == 0 ){
|
||||
KEA = reactionConfig.beamEnergy;
|
||||
}else{
|
||||
KEA = gRandom->Gaus(reactionConfig.beamEnergy, reactionConfig.beamEnergySigma);
|
||||
if( reactConfig.beamEnergySigma != 0 ){
|
||||
KEA = gRandom->Gaus(reactConfig.beamEnergy, reactConfig.beamEnergySigma);
|
||||
}
|
||||
theta = 0.0;
|
||||
if( reactionConfig.beamAngleSigma == 0 ){
|
||||
theta = reactionConfig.beamAngle;
|
||||
}else{
|
||||
theta = gRandom->Gaus(reactionConfig.beamAngle, reactionConfig.beamAngleSigma);
|
||||
if( reactConfig.beamThetaSigma != 0 ){
|
||||
theta = gRandom->Gaus(reactConfig.beamTheta, reactConfig.beamThetaSigma);
|
||||
}
|
||||
phi = 0.0;
|
||||
|
||||
//==== for taregt scattering
|
||||
transfer.SetIncidentEnergyAngle(KEA, theta, 0.);
|
||||
transfer.CalReactionConstant();
|
||||
TLorentzVector PA = transfer.GetPA();
|
||||
|
||||
// TLorentzVector PA = transfer.GetPA();
|
||||
//depth = 0;
|
||||
if( isTargetScattering ){
|
||||
//==== Target scattering, only energy loss
|
||||
depth = targetThickness * gRandom->Rndm();
|
||||
msA.SetTarget(density, depth);
|
||||
TLorentzVector PAnew = msA.Scattering(PA);
|
||||
KEAnew = msA.GetKE()/reactionConfig.beamA;
|
||||
transfer.SetIncidentEnergyAngle(KEAnew, theta, phi);
|
||||
transfer.CalReactionConstant();
|
||||
Ecm = transfer.GetCMTotalKE();
|
||||
}
|
||||
// if( isTargetScattering ){
|
||||
// //==== Target scattering, only energy loss
|
||||
// depth = targetThickness * gRandom->Rndm();
|
||||
// msA.SetTarget(density, depth);
|
||||
// TLorentzVector PAnew = msA.Scattering(PA);
|
||||
// KEAnew = msA.GetKE()/reactConfig.beamA;
|
||||
// transfer.SetIncidentEnergyAngle(KEAnew, theta, phi);
|
||||
// transfer.CalReactionConstant();
|
||||
// Ecm = transfer.GetCMTotalKE();
|
||||
// }
|
||||
|
||||
//==== Calculate thetaCM, phiCM
|
||||
if( distFile->IsOpen()){
|
||||
|
|
@ -545,25 +478,24 @@ void Transfer(
|
|||
TLorentzVector Pb = transfer.GetPb();
|
||||
TLorentzVector PB = transfer.GetPB();
|
||||
|
||||
//==== Calculate energy loss of scattered and recoil in target
|
||||
if( isTargetScattering ){
|
||||
if( Pb.Theta() < TMath::PiOver2() ){
|
||||
msb.SetTarget(density, targetThickness - depth);
|
||||
}else{
|
||||
msb.SetTarget(density, depth);
|
||||
}
|
||||
Pb = msb.Scattering(Pb);
|
||||
TbLoss = msb.GetKELoss();
|
||||
msB.SetTarget(density, targetThickness - depth);
|
||||
PB = msB.Scattering(PB);
|
||||
}else{
|
||||
TbLoss = 0;
|
||||
}
|
||||
// //==== Calculate energy loss of scattered and recoil in target
|
||||
// if( isTargetScattering ){
|
||||
// if( Pb.Theta() < TMath::PiOver2() ){
|
||||
// msb.SetTarget(density, targetThickness - depth);
|
||||
// }else{
|
||||
// msb.SetTarget(density, depth);
|
||||
// }
|
||||
// Pb = msb.Scattering(Pb);
|
||||
// TbLoss = msb.GetKELoss();
|
||||
// msB.SetTarget(density, targetThickness - depth);
|
||||
// PB = msB.Scattering(PB);
|
||||
// }else{
|
||||
// TbLoss = 0;
|
||||
// }
|
||||
|
||||
//======= Decay of particle-B
|
||||
int decayID = 0;
|
||||
int new_zB = reactionConfig.recoilHeavyZ;
|
||||
if( reactionConfig.isDecay){
|
||||
if( recoil.isDecay){
|
||||
|
||||
//decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2.); // decay to ground state
|
||||
decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state
|
||||
|
|
@ -571,7 +503,7 @@ void Transfer(
|
|||
PB = decay.GetDaugther_D();
|
||||
//decayTheta = decay.GetAngleChange();
|
||||
decayTheta = decay.GetThetaCM();
|
||||
new_zB = reactionConfig.heavyDecayZ;
|
||||
PB.SetUniqueID(recoil.decayZ);
|
||||
}else{
|
||||
decayTheta = TMath::QuietNaN();
|
||||
}
|
||||
|
|
@ -590,20 +522,21 @@ void Transfer(
|
|||
|
||||
//==== Helios
|
||||
|
||||
///printf(" thetaCM : %f \n", thetaCM * TMath::RadToDeg());
|
||||
// printf(" thetaCM : %f, Tb : %f\n", thetaCM * TMath::RadToDeg(), Pb.M());
|
||||
|
||||
if( Tb > 0 || TB > 0 ){
|
||||
helios.CalArrayHit(Pb, reactionConfig.recoilLightZ);
|
||||
helios.CalRecoilHit(PB, new_zB);
|
||||
|
||||
helios.CalArrayHit(Pb);
|
||||
helios.CalRecoilHit(PB);
|
||||
hit = 2;
|
||||
while( hit > 1 ){ hit = helios.DetAcceptance(); } /// while hit > 1, goto next loop;
|
||||
while( hit > 1 ){ hit = helios.CheckDetAcceptance(); } /// while hit > 1, goto next loop;
|
||||
|
||||
trajectory orb_b = helios.GetTrajectory_b();
|
||||
trajectory orb_B = helios.GetTrajectory_B();
|
||||
|
||||
e = helios.GetEnergy() + gRandom->Gaus(0, detGeo.array1.eSigma);
|
||||
e = helios.GetEnergy() + gRandom->Gaus(0, array.eSigma );
|
||||
|
||||
double ranX = gRandom->Gaus(0, detGeo.array1.zSigma);
|
||||
double ranX = gRandom->Gaus(0, array.zSigma);
|
||||
z = orb_b.z + ranX;
|
||||
detX = helios.GetDetX() + ranX;
|
||||
|
||||
|
|
@ -658,13 +591,13 @@ void Transfer(
|
|||
thetaCM = thetaCM * TMath::RadToDeg();
|
||||
|
||||
//if decay, get the light decay particle on the recoil;
|
||||
if( reactionConfig.isDecay ){
|
||||
if( recoil.isDecay ){
|
||||
if( decayID == 1 ){
|
||||
TLorentzVector Pd = decay.GetDaugther_d();
|
||||
|
||||
Td = Pd.E() - Pd.M();
|
||||
|
||||
helios.CalRecoilHit(Pd, reactionConfig.heavyDecayZ);
|
||||
helios.CalRecoilHit(Pd);
|
||||
|
||||
trajectory orb_d = helios.GetTrajectory_B();
|
||||
rhoRecoil_d = orb_d.R;
|
||||
|
|
@ -684,7 +617,7 @@ void Transfer(
|
|||
|
||||
if( hit == 1) count ++;
|
||||
|
||||
if( reactionConfig.isRedo ){
|
||||
if( reactConfig.isRedo ){
|
||||
if( hit == 1) {
|
||||
redoFlag = false;
|
||||
}else{
|
||||
|
|
@ -720,7 +653,11 @@ void Transfer(
|
|||
saveFile->Close();
|
||||
|
||||
distFile->Close();
|
||||
delete exDist;
|
||||
|
||||
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
|
||||
//gROOT->ProcessLine(".q");
|
||||
|
||||
return;
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,9 +1,12 @@
|
|||
CC=g++
|
||||
|
||||
ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray FindThetaCM Transfer PlotSimulation
|
||||
ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM Transfer
|
||||
|
||||
all: $(ALL)
|
||||
|
||||
Isotope: ../Cleopatra/ClassIsotope.h ../Cleopatra/Isotope.C
|
||||
$(CC) Isotope.C -o Isotope
|
||||
|
||||
InFileCreator: InFileCreator.C InFileCreator.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h potentials.h
|
||||
$(CC) InFileCreator.C -o InFileCreator `root-config --cflags --glibs`
|
||||
|
||||
|
|
@ -16,17 +19,14 @@ ExtractXSecFromText: ExtractXSecFromText.C ExtractXSec.h
|
|||
PlotTGraphTObjArray: PlotTGraphTObjArray.C PlotTGraphTObjArray.h
|
||||
$(CC) PlotTGraphTObjArray.C -o PlotTGraphTObjArray `root-config --cflags --glibs`
|
||||
|
||||
Cleopatra: Cleopatra.C
|
||||
$(CC) Cleopatra.C -o Cleopatra `root-config --cflags --glibs`
|
||||
|
||||
FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
|
||||
$(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs`
|
||||
|
||||
Transfer: Transfer.C Transfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
|
||||
Transfer: Transfer.C Transfer.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
|
||||
$(CC) Transfer.C -o Transfer `root-config --cflags --glibs`
|
||||
|
||||
PlotSimulation: PlotSimulation.C Check_Simulation.C
|
||||
$(CC) PlotSimulation.C -o PlotSimulation `root-config --cflags --glibs`
|
||||
|
||||
Isotope: ../Cleopatra/ClassIsotope.h ../Cleopatra/Isotope.C
|
||||
$(CC) Isotope.C -o Isotope
|
||||
|
||||
clean:
|
||||
/bin/rm -f $(ALL)
|
||||
30
README.md
30
README.md
|
|
@ -5,23 +5,37 @@ This is the analysis package for the SOLARIS DAQ. It is supposed to be the analy
|
|||
The folder struture is
|
||||
|
||||
Analysis
|
||||
|
||||
├── README.md
|
||||
|
||||
├── SetupNewExp // bash script to create new branch and raw data folder
|
||||
|
||||
├── SOLARIS.sh // bash script to define some env variable and functions
|
||||
|
||||
├── Armory // analysis codes, independent from experiment.
|
||||
|
||||
├── Cleopatra // Swaper for DWBA code Ptolomey and simulation
|
||||
|
||||
├── data_raw // should be the symbolic link to the raw data, created by SetUpNewExp
|
||||
|
||||
├── root_data // symbolic link to converted root file, created by SetUpNewExp
|
||||
|
||||
└── working // working directory, depends on experiment.
|
||||
|
||||
# Analysis & Simulation
|
||||
|
||||
The Armory/AnalysisLib.h constains many small but handy functions.
|
||||
|
||||
All class headers are started with Class*.h
|
||||
|
||||
The classes **DetGeo**** and **ReactionConfig** are fundamental for loading the detectorGeo.txt and reactionConfig.txt.
|
||||
|
||||
Both txt file support empty lines, and up to 2 settings. The reason for that is for dual-array configuration. It has potentail to extend and include more settings. But it is two now, one for upstream array (reaction) and downstream array (reaction).
|
||||
|
||||
The **TransferReaction** class is only use one of the reaction from the reactionConfig.txt.
|
||||
|
||||
```C++
|
||||
TransferReaction::SetReactionFromFile("reactionConfig.txt", ID); // ID = 0 or 1
|
||||
```
|
||||
Same for the **Helios** class
|
||||
|
||||
```C++
|
||||
HELIOS::SetDetectorGeometry("detectorGeo.txt", ID); // ID = 0 or 1
|
||||
```
|
||||
|
||||
|
||||
# Event Builder
|
||||
|
||||
The EventBuilder is at the armory. It depends on the Hit.h and SolReader.h.
|
||||
|
|
|
|||
|
|
@ -44,3 +44,6 @@
|
|||
#32Si(t,p)34Si 0 0L=0 0+ 0.000 8MeV/u lA #two-nucleon_transfer
|
||||
#36Ar(d,a)34Cl 0 4L=2 3+ 0.000 8MeV/u As # (d,a) reaction
|
||||
|
||||
|
||||
20F(d,t)19F 2 1s1/2 3/2+ 0.000 10MeV/u Vl
|
||||
20F(d,t)19F 2 0d5/2 5/2+ 0.197 10MeV/u Vl
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
//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
|
||||
|
|
@ -1,3 +1,5 @@
|
|||
#--- '#---' comment line identifier
|
||||
#--- beam
|
||||
32 //beam_A
|
||||
14 //beam_Z
|
||||
0.0 //excitation_energy_of_beam[MeV]
|
||||
|
|
@ -8,6 +10,7 @@
|
|||
0.00 //x_offset_of_Beam_in_mm
|
||||
0.00 //y_offset_of_Beam_in_mm
|
||||
|
||||
#--- target
|
||||
2 //target_A
|
||||
1 //target_Z
|
||||
false //isTargetScattering
|
||||
|
|
@ -15,7 +18,8 @@ false //isTargetScattering
|
|||
2.2e-4 //targetThickness_in_cm
|
||||
../SRIM/20F_in_CD2.txt //stopping_power_for_beam
|
||||
|
||||
100000 //number_of_Event_being_generated
|
||||
#--- Monte Carlo Setting
|
||||
1000000 //number_of_Event_being_generated
|
||||
false //Redo_until_hit_array=all_events_hit_array
|
||||
|
||||
#=====reaction_for_1st_Array
|
||||
|
|
@ -26,6 +30,11 @@ false //Redo_until_hit_array=all_events_hit_array
|
|||
false //isDacay
|
||||
32 //decayNucleus_A
|
||||
14 //decayNucleus_Z
|
||||
#--- List of Ex of heavy recoil
|
||||
#---Ex relative_xsec SF sigma_in_MeV
|
||||
0.000 1.0 1.0 0.01
|
||||
1.000 1.0 1.0 0.01
|
||||
2.000 1.0 1.0 0.01
|
||||
|
||||
#=====_reaction_for_2nd_Array_use_ony_when_2nd_array_used
|
||||
3 //recoil_light_A
|
||||
|
|
@ -36,4 +45,12 @@ false //isDacay
|
|||
32 //decayNucleus_A
|
||||
14 //decayNucleus_Z
|
||||
|
||||
#--- List of Ex of heavy recoil
|
||||
#---Ex relative_xsec SF sigma_in_MeV
|
||||
|
||||
0.000 1.0 1.0 0.01
|
||||
1.000 1.0 1.0 0.01
|
||||
2.000 1.0 1.0 0.01
|
||||
|
||||
|
||||
################## end of file
|
||||
Loading…
Reference in New Issue
Block a user