diff --git a/.gitignore b/.gitignore index 73ab6e2..71270d2 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,10 @@ data data_raw root_data +*.in +*.out +*.txt + Cleopatra/ExtractXSec Cleopatra/ExtractXSecFromText Cleopatra/FindThetaCM diff --git a/Armory/ClassReactionConfig.h b/Armory/ClassReactionConfig.h index 145b204..b24c02d 100644 --- a/Armory/ClassReactionConfig.h +++ b/Armory/ClassReactionConfig.h @@ -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 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(); } diff --git a/Cleopatra/ClassDecay.h b/Cleopatra/ClassDecay.h index b313ae3..47bb6dd 100644 --- a/Cleopatra/ClassDecay.h +++ b/Cleopatra/ClassDecay.h @@ -4,6 +4,7 @@ #include "TVector3.h" #include "../Cleopatra/ClassIsotope.h" +#include "../Armory/ClassReactionConfig.h" //======================================================= //####################################################### @@ -14,133 +15,164 @@ //======================================================= class Decay{ public: - Decay(); - ~Decay(); - - double GetQValue() { return Q;} - - double GetAngleChange(){ - TVector3 vD = PD.Vect(); - TVector3 vB = PB.Vect(); - vD.SetMag(1); - vB.SetMag(1); - double dot = vD.Dot(vB); - return TMath::ACos(dot)*TMath::RadToDeg() ; - } - - double GetThetaCM() { return theta * TMath::RadToDeg();} + Decay(); + ~Decay(); + + double GetQValue() { return Q;} + + double GetAngleChange(){ + TVector3 vD = PD.Vect(); + TVector3 vB = PB.Vect(); + vD.SetMag(1); + vB.SetMag(1); + double dot = vD.Dot(vB); + return TMath::ACos(dot)*TMath::RadToDeg() ; + } + + double GetThetaCM() { return theta * TMath::RadToDeg();} + + double GetCMMomentum(){ return k;} + TLorentzVector GetDaugther_d() {return Pd;} + TLorentzVector GetDaugther_D() {return PD;} + + void SetMotherDaugther(Recoil recoil){ - double GetCMMomentum(){ return k;} - TLorentzVector GetDaugther_d() {return Pd;} - TLorentzVector GetDaugther_D() {return PD;} - - void SetMotherDaugther(int AB, int zB, int AD, int zD){ - Isotope Mother(AB, zB); - Isotope Daugther_D(AD, zD); - Isotope Daugther_d(AB-AD, zB-zD); + Isotope Mother(recoil.heavyA, recoil.heavyZ); + Isotope Daugther_D(recoil.decayA, recoil.decayZ); + Isotope Daugther_d(recoil.heavyA - recoil.decayA, 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; - } - - int CalDecay(TLorentzVector P_mother, double ExB, double ExD, double normOfReactionPlane = 0){ - if( !isMotherSet ) { - return -1; - } - this->PB = P_mother; - - double MB = mB + ExB; ///mother - double MD = mD + ExD; ///Big_Daugther - Q = MB - MD - md; - if( Q < 0 ) { - this->PD = this->PB; - dTheta = TMath::QuietNaN(); - k = TMath::QuietNaN(); - return -2; - } - - //clear - TLorentzVector temp(0,0,0,0); - PD = temp; - Pd = temp; - - k = TMath::Sqrt((MB+MD+md)*(MB+MD-md)*(MB-MD+md)*(MB-MD-md))/2./MB; - - //in mother's frame, assume isotropic decay - theta = TMath::ACos(2 * gRandom->Rndm() - 1) ; - - //for non isotropic decay, edit f1. - //theta = TMath::ACos(f1->GetRandom()); - - double phi = TMath::TwoPi() * gRandom->Rndm(); - PD.SetE(TMath::Sqrt(mD * mD + k * k )); - PD.SetPz(k); - PD.SetTheta(theta); - PD.SetPhi(phi); - - Pd.SetE(TMath::Sqrt(md * md + k * k )); - Pd.SetPz(k); - Pd.SetTheta(theta + TMath::Pi()); - Pd.SetPhi(phi + TMath::Pi()); - - PD.RotateY(TMath::Pi()/2.); - PD.RotateZ(normOfReactionPlane); + zB = recoil.heavyZ; + zD = recoil.decayZ; + zd = recoil.heavyZ - recoil.decayZ; - Pd.RotateY(TMath::Pi()/2.); - Pd.RotateZ(normOfReactionPlane); - - //Transform to Lab frame; - TVector3 boost = PB.BoostVector(); - - PD.Boost(boost); - Pd.Boost(boost); - - return 1; - } - + 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); + Isotope Daugther_d(AB-AD, zB-zD); + + 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; + } + + int CalDecay(TLorentzVector P_mother, double ExB, double ExD, double normOfReactionPlane = 0){ + if( !isMotherSet ) { + return -1; + } + this->PB = P_mother; + + double MB = mB + ExB; ///mother + double MD = mD + ExD; ///Big_Daugther + Q = MB - MD - md; + if( Q < 0 ) { + this->PD = this->PB; + dTheta = TMath::QuietNaN(); + k = TMath::QuietNaN(); + return -2; + } + + //clear + TLorentzVector temp(0,0,0,0); + 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 + theta = TMath::ACos(2 * gRandom->Rndm() - 1) ; + + //for non isotropic decay, edit f1. + //theta = TMath::ACos(f1->GetRandom()); + + double phi = TMath::TwoPi() * gRandom->Rndm(); + PD.SetE(TMath::Sqrt(mD * mD + k * k )); + PD.SetPz(k); + PD.SetTheta(theta); + PD.SetPhi(phi); + + Pd.SetE(TMath::Sqrt(md * md + k * k )); + Pd.SetPz(k); + Pd.SetTheta(theta + TMath::Pi()); + Pd.SetPhi(phi + TMath::Pi()); + + PD.RotateY(TMath::Pi()/2.); + PD.RotateZ(normOfReactionPlane); + + Pd.RotateY(TMath::Pi()/2.); + Pd.RotateZ(normOfReactionPlane); + + //Transform to Lab frame; + TVector3 boost = PB.BoostVector(); + + PD.Boost(boost); + Pd.Boost(boost); + + return 1; + } + private: - TLorentzVector PB, Pd, PD; - - double mB, mD, md; - double theta; - - TF1 * f1; - - bool isMotherSet; - double Q; - double k; // momentum in B-frame - double dTheta; // change of angle + TLorentzVector PB, Pd, PD; + + double mB, mD, md; + double zB, zD, zd; + double theta; + + TF1 * f1; + + bool isMotherSet; + double Q; + double k; // momentum in B-frame + double dTheta; // change of angle }; Decay::Decay(){ - TLorentzVector temp(0,0,0,0); - PB = temp; - Pd = temp; - PD = temp; + TLorentzVector temp(0,0,0,0); + PB = temp; + Pd = temp; + PD = temp; + + mB = TMath::QuietNaN(); + mD = TMath::QuietNaN(); + md = TMath::QuietNaN(); + + zB = 0; + zD = 0; + zd = 0; + + theta = TMath::QuietNaN(); - mB = TMath::QuietNaN(); - mD = TMath::QuietNaN(); - md = TMath::QuietNaN(); - theta = TMath::QuietNaN(); - - k = TMath::QuietNaN(); - - Q = TMath::QuietNaN(); - dTheta = TMath::QuietNaN(); - isMotherSet = false; - - f1 = new TF1("f1", "(1+ROOT::Math::legendre(2,x))/2.", -1, 1); + k = TMath::QuietNaN(); + + Q = TMath::QuietNaN(); + dTheta = TMath::QuietNaN(); + isMotherSet = false; + + f1 = new TF1("f1", "(1+ROOT::Math::legendre(2,x))/2.", -1, 1); } Decay::~Decay(){ - delete f1; + delete f1; } diff --git a/Cleopatra/ClassHelios.h b/Cleopatra/ClassHelios.h index 30ea62b..ed1a908 100644 --- a/Cleopatra/ClassHelios.h +++ b/Cleopatra/ClassHelios.h @@ -133,34 +133,40 @@ 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: - DetGeo detGeo; - Array array; - - trajectory orbitb, orbitB; - - double e,detX ; ///energy of light recoil, position X - double rhoHit; /// radius of particle-b hit on recoil detector - - double eB; ///energy of heavy recoil - - bool isDetReady; + DetGeo detGeo; + Array array; + + trajectory orbitb, orbitB; + + double e,detX ; ///energy of light recoil, position X + double rhoHit; /// radius of particle-b hit on recoil detector + + double eB; ///energy of heavy recoil + + bool isDetReady; + + TString hitMessage; + TString acceptanceMsg; //acceptance check + short acceptanceCode; + + double xOff, yOff; // beam position + + bool overrideDetDistance; + bool overrideFirstPos; + bool isCoincidentWithRecoil; + + const double c = 299.792458; //mm/ns - TString hitMessage; - TString accMessage; //acceptance check - - double xOff, yOff; // beam position - bool overrideDetDistance; - bool overrideFirstPos; - 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,108 +264,137 @@ 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; + //CalArrayHit and CalRecoilHit must be done before. + + 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;} + // -1 ========= when recoil direction is not same side of array + 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;} - - // -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist) - if( orbitb.detRowID == -1 ) {accMessage = "det Row ID == -1"; return -15;} - - // -10 =========== when rho is too big . - if( detGeo.bore < 2 * orbitb.rho) { accMessage = "rho is too big"; return -10;} - - // -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 > 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;} - - // -13 ========= not more than 3 loops - if( orbitb.loop > 3 ) {accMessage = "more than 3 loops."; return -13;} - - // -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 ==== 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;} + // -11 ======== rho is too small + if( 2 * orbitb.rho < array.detPerpDist ) { acceptanceCode = -11; return acceptanceCode;} + + // -15 ========= if detRowID == -1, should be (2 * orbitb.rho < perpDist) + if( orbitb.detRowID == -1 ) {acceptanceCode = -15; return acceptanceCode;} + + // -10 =========== when rho is too big . + 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 ) {acceptanceCode = -14; 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;} + //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 ) { 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 ) {acceptanceCode = -13; return acceptanceCode;} + + // -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;} + + // -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;} - // 2 ====== when zPos less then the nearest position, more loop may hit - int increaseLoopFlag = 0; - if( array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z ) increaseLoopFlag = 2; - if( array.firstPos > 0 && array.detPos[0] > orbitb.z ) increaseLoopFlag = 2; - if (increaseLoopFlag == 2 ) { - orbitb.z += orbitb.z0; - orbitb.effLoop += 1.0; - orbitb.loop += 1; - orbitb.t = orbitb.t0 * orbitb.effLoop; - accMessage = " hit less than the nearest array. increase loop "; - return 2; - } + // -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; + if( array.firstPos < 0 && array.detPos[array.nDet-1] < orbitb.z ) increaseLoopFlag = 2; + if( array.firstPos > 0 && array.detPos[0] > orbitb.z ) increaseLoopFlag = 2; + if (increaseLoopFlag == 2 ) { + orbitb.z += orbitb.z0; + orbitb.effLoop += 1.0; + orbitb.loop += 1; + orbitb.t = orbitb.t0 * orbitb.effLoop; + acceptanceCode = 2; + return acceptanceCode; + } - // 1 ======= check hit array z- position - if( array.firstPos < 0 ){ - for( int i = 0; i < array.nDet; i++){ - 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; - } - } - }else{ - for( int i = 0; i < array.nDet ; i++){ - if( array.detPos[i] <= orbitb.z && orbitb.z <= array.detPos[i] + array.detLength) { - ///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; - } + // 1 ======= check hit array z- position + if( array.firstPos < 0 ){ + for( int i = 0; i < array.nDet; i++){ + 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 + acceptanceCode = 1; + return acceptanceCode; } - } + } + }else{ + for( int i = 0; i < array.nDet ; i++){ + if( array.detPos[i] <= orbitb.z && orbitb.z <= array.detPos[i] + array.detLength) { + ///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 + acceptanceCode = 1; + return acceptanceCode; + } + } + } - // -5 ======== 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; - } - }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 (increaseLoopFlag == 3 ) { - orbitb.z += orbitb.z0; - orbitb.effLoop += 1.0; - orbitb.loop += 1; - orbitb.t = orbitb.t0 * orbitb.effLoop; - accMessage = " try one more loop. "; - return 3; - } + // -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 ) { 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] ) { acceptanceCode = -7; return acceptanceCode; }//increaseLoopFlag = 3; + } + } + if (increaseLoopFlag == 3 ) { + orbitb.z += orbitb.z0; + orbitb.effLoop += 1.0; + orbitb.loop += 1; + orbitb.t = orbitb.t0 * orbitb.effLoop; + 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){ diff --git a/Cleopatra/ClassTransfer.h b/Cleopatra/ClassTransfer.h index c4c38d5..3578e66 100644 --- a/Cleopatra/ClassTransfer.h +++ b/Cleopatra/ClassTransfer.h @@ -22,18 +22,10 @@ class TransferReaction { public: 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 configFile, unsigned short ID = 0){ - Inititization(); - SetReactionFromFile(configFile, ID); - } + TransferReaction(string configFile, unsigned short ID = 0); + TransferReaction(int beamA, int beamZ, + int targetA, int targetZ, + 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;} + 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 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 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.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; diff --git a/Cleopatra/Cleopatra b/Cleopatra/Cleopatra new file mode 100755 index 0000000..7dea56f Binary files /dev/null and b/Cleopatra/Cleopatra differ diff --git a/Cleopatra/Cleopatra.C b/Cleopatra/Cleopatra.C index 40133ae..f20b83b 100644 --- a/Cleopatra/Cleopatra.C +++ b/Cleopatra/Cleopatra.C @@ -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); diff --git a/Cleopatra/ExtractXSec.h b/Cleopatra/ExtractXSec.h index 668de9b..ae85cf9 100644 --- a/Cleopatra/ExtractXSec.h +++ b/Cleopatra/ExtractXSec.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #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()); } @@ -338,6 +322,14 @@ int ExtractXSec (string readFile, int indexForElastic=1) { fprintf(file_out, "\n"); } 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(); @@ -345,8 +337,8 @@ int ExtractXSec (string readFile, int indexForElastic=1) { TString fileName = saveFileName; fileName += ".root"; printf("Output : %s \n", fileName.Data()); + 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(); diff --git a/Cleopatra/FindThetaCM.h b/Cleopatra/FindThetaCM.h index 3e124df..6d74cd4 100644 --- a/Cleopatra/FindThetaCM.h +++ b/Cleopatra/FindThetaCM.h @@ -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"); diff --git a/Cleopatra/PlotSimulation.C b/Cleopatra/PlotSimulation.C deleted file mode 100644 index e3080c1..0000000 --- a/Cleopatra/PlotSimulation.C +++ /dev/null @@ -1,31 +0,0 @@ - -#include -#include -#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); - -} diff --git a/Cleopatra/PlotTGraphTObjArray.h b/Cleopatra/PlotTGraphTObjArray.h index f6d3403..d646f5b 100644 --- a/Cleopatra/PlotTGraphTObjArray.h +++ b/Cleopatra/PlotTGraphTObjArray.h @@ -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"); diff --git a/Cleopatra/Transfer.C b/Cleopatra/Transfer.C index ffe99cb..3f851e3 100644 --- a/Cleopatra/Transfer.C +++ b/Cleopatra/Transfer.C @@ -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 ){ diff --git a/Cleopatra/Transfer.h b/Cleopatra/Transfer.h index 861b15b..ef98d46 100644 --- a/Cleopatra/Transfer.h +++ b/Cleopatra/Transfer.h @@ -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 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()); + + printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n"); - const ReactionConfig reactionConfig = transfer.GetRectionConfig(); + transfer.PrintReaction(); - reactionConfig.PrintReactionConfig(); + if(transfer.GetRecoil().isDecay) { + decay.SetMotherDaugther(transfer.GetRecoil()); + } - vector ExAList = reactionConfig.beamEx; - int nExA = (int) ExAList.size(); + helios.PrintGeometry(); + PrintEZPlotPara(transfer, helios); - //############################################# Set HELIOS - printf("\e[32m#################################### HELIOS configuration\e[0m\n"); - HELIOS helios; - helios.SetDetectorGeometry(heliosDetGeoFile); - const DetGeo detGeo = helios.GetDetectorGeometry(); + DetGeo detGeo = helios.GetDetectorGeometry(); + Array array = helios.GetArrayGeometry(); + ReactionConfig reactConfig = transfer.GetRectionConfig(); + Recoil recoil = transfer.GetRecoil(); - 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+"); - //############################################# 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", 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"); - 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); - } + // fflush(keyParaOut); + // fclose(keyParaOut); + // } - //############################################# Target scattering, only energy loss - bool isTargetScattering = reactionConfig.isTargetScattering; - float density = reactionConfig.targetDensity; - float targetThickness = reactionConfig.targetThickness; + //*############################################# Target scattering, only energy loss + // bool isTargetScattering = reactConfig.isTargetScattering; + // float density = reactConfig.targetDensity; + // float targetThickness = reactConfig.targetThickness; - if(isTargetScattering) printf("\e[32m#################################### Target Scattering\e[0m\n"); - TargetScattering msA; - TargetScattering msB; - TargetScattering msb; + // 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); + // 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); - if( reactionConfig.isTargetScattering ){ - msA.LoadStoppingPower(reactionConfig.beamStoppingPowerFile); - msb.LoadStoppingPower(reactionConfig.recoilLightStoppingPowerFile); - msB.LoadStoppingPower(reactionConfig.recoilHeavyStoppingPowerFile); - } + // if( reactConfig.isTargetScattering ){ + // msA.LoadStoppingPower(reactConfig.beamStoppingPowerFile); + // msb.LoadStoppingPower(reactConfig.recoilLightStoppingPowerFile); + // msB.LoadStoppingPower(reactConfig.recoilHeavyStoppingPowerFile); + // } - //############################################# 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 ExKnown; - vector ExStrength; - vector ExWidth; - vector SF; - vector y0; /// intercept of e-z plot - vector 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; + //*############################################# 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); + // } - vector str = AnalysisLib::SplitStr(line, " "); - - 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())); - - } - file.close(); - printf("... done.\n"); - int n = (int) ExKnown.size(); - - 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); - } - - //---- 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]); - } - } + 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 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 @@ -274,11 +233,6 @@ void Transfer( double rho, rhoB; ///orbit radius 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; @@ -288,26 +242,21 @@ void Transfer( double ExCal, thetaCMCal; 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 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; + } diff --git a/Cleopatra/makefile b/Cleopatra/makefile index 26ce40c..f253b1e 100644 --- a/Cleopatra/makefile +++ b/Cleopatra/makefile @@ -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) \ No newline at end of file diff --git a/README.md b/README.md index 6ee3403..5207adf 100644 --- a/README.md +++ b/README.md @@ -4,24 +4,38 @@ 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 - +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. diff --git a/working/DWBA b/working/DWBA index c94977f..f380e7a 100644 --- a/working/DWBA +++ b/working/DWBA @@ -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 \ No newline at end of file diff --git a/working/Ex.txt b/working/Ex.txt deleted file mode 100644 index 160e6c1..0000000 --- a/working/Ex.txt +++ /dev/null @@ -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 diff --git a/working/reactionConfig.txt b/working/reactionConfig.txt index 5f0662d..cd9fa6c 100644 --- a/working/reactionConfig.txt +++ b/working/reactionConfig.txt @@ -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 \ No newline at end of file