added Monte Carlo simulation
This commit is contained in:
parent
18fe2ed62c
commit
caff146185
12
.vscode/c_cpp_properties.json
vendored
12
.vscode/c_cpp_properties.json
vendored
|
|
@ -24,6 +24,18 @@
|
|||
"cStandard": "c17",
|
||||
"cppStandard": "gnu++17",
|
||||
"intelliSenseMode": "linux-gcc-x64"
|
||||
},
|
||||
{
|
||||
"name": "RyanUbuntu",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/opt/root/**"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
"cStandard": "c17",
|
||||
"cppStandard": "gnu++17",
|
||||
"intelliSenseMode": "linux-gcc-x64"
|
||||
}
|
||||
],
|
||||
"version": 4
|
||||
|
|
|
|||
|
|
@ -69,15 +69,15 @@ private:
|
|||
|
||||
const int nWire = 24;
|
||||
const int wireShift = 3;
|
||||
const int zLen = 380; //mm
|
||||
const int radiusA = 37;
|
||||
const int radiusC = 43;
|
||||
const float zLen = 380; //mm
|
||||
const float radiusA = 37;
|
||||
const float radiusC = 43;
|
||||
|
||||
std::vector<std::pair<TVector3,TVector3>> P1; // the anode wire position vector in space
|
||||
std::vector<std::pair<TVector3,TVector3>> Q1; // the cathode wire position vector in space
|
||||
|
||||
std::vector<std::pair<TVector3,TVector3>> S1; // coners of the SX3 0-11, z = mid point
|
||||
std::vector<std::pair<TVector3,TVector3>> S2; // coners of the SX3 12-23, z = mid point
|
||||
std::vector<std::pair<TVector3,TVector3>> SD; // coners of the SX3 0-11, z = mid point
|
||||
std::vector<std::pair<TVector3,TVector3>> SU; // coners of the SX3 12-23, z = mid point
|
||||
std::vector<TVector3> SNorml; // normal of the SX3 (outward)
|
||||
|
||||
void CalGeometry();
|
||||
|
|
@ -92,14 +92,14 @@ private:
|
|||
TVector3 trackVecErrorC; // error vector prependicular to the Cathode-Pos plan
|
||||
|
||||
const int nSX3 = 12;
|
||||
const int sx3Radius = 88;
|
||||
const int sx3Width = 40;
|
||||
const int sx3Length = 75;
|
||||
const int sx3Gap = 46;
|
||||
const float sx3Radius = 88;
|
||||
const float sx3Width = 40;
|
||||
const float sx3Length = 75;
|
||||
const float sx3Gap = 46;
|
||||
|
||||
const int qqqR1 = 50;
|
||||
const int qqqR2 = 100;
|
||||
const int qqqZPos = sx3Gap/2 + sx3Length + 30;
|
||||
const float qqqR1 = 50;
|
||||
const float qqqR2 = 100;
|
||||
const float qqqZPos = sx3Gap/2 + sx3Length + 30;
|
||||
|
||||
// int geomID;
|
||||
TGeoManager *geom;
|
||||
|
|
@ -168,7 +168,7 @@ inline void ANASEN::CalGeometry(){
|
|||
|
||||
sa.RotateZ( TMath::TwoPi() / nSX3 * (i + 0.5) );
|
||||
sb.RotateZ( TMath::TwoPi() / nSX3 * (i + 0.5) );
|
||||
S1.push_back(std::pair(sa,sb));
|
||||
SD.push_back(std::pair(sa,sb));
|
||||
|
||||
|
||||
sc.SetXYZ( sx3Radius, -sx3Width/2, sx3Gap/2 );
|
||||
|
|
@ -182,7 +182,7 @@ inline void ANASEN::CalGeometry(){
|
|||
|
||||
sa.RotateZ( TMath::TwoPi() / nSX3 * (i + 0.5) );
|
||||
sb.RotateZ( TMath::TwoPi() / nSX3 * (i + 0.5) );
|
||||
S2.push_back(std::pair(sa,sb));
|
||||
SU.push_back(std::pair(sa,sb));
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -344,22 +344,47 @@ inline std::pair<int, int> ANASEN::FindWireID(TVector3 pos, TVector3 direction,
|
|||
|
||||
double phi = direction.Phi();
|
||||
|
||||
|
||||
for( int i = 0; i < nWire; i++){
|
||||
|
||||
double disA = 99999999;
|
||||
double disC = 99999999;
|
||||
|
||||
if(P1[i].first.Phi()-TMath::PiOver2() < phi && phi < P1[i].second.Phi()+TMath::PiOver2()) {
|
||||
disA = Distance( pos, pos + direction, P1[i].first, P1[i].second);
|
||||
double phiS = P1[i].first.Phi() - TMath::PiOver4();
|
||||
double phiL = P1[i].second.Phi() + TMath::PiOver4();
|
||||
|
||||
// printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg());
|
||||
|
||||
if( phi > 0 && phiS > phiL ) {
|
||||
phiL = phiL + TMath::TwoPi();
|
||||
// printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
|
||||
}
|
||||
if( phi < 0 && phiS > phiL ) {
|
||||
phiS = phiS - TMath::TwoPi();
|
||||
// printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
|
||||
}
|
||||
|
||||
if( phiS < phi && phi < phiL) {
|
||||
disA = Distance( pos, pos + direction, P1[i].first, P1[i].second);
|
||||
if( disA < minAnodeDis ){
|
||||
minAnodeDis = disA;
|
||||
anodeID = i;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if(Q1[i].second.Phi() - TMath::PiOver2() < phi && phi < Q1[i].first.Phi() + TMath::PiOver2()) {
|
||||
phiS = Q1[i].second.Phi()- TMath::PiOver4();
|
||||
phiL = Q1[i].first.Phi() + TMath::PiOver4();
|
||||
// printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
|
||||
if( phi > 0 && phiS > phiL ) {
|
||||
phiL = phiL + TMath::TwoPi();
|
||||
// printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
|
||||
}
|
||||
if( phi < 0 && phiS > phiL ) {
|
||||
phiS = phiS - TMath::TwoPi();
|
||||
// printf("------ %f %f\n", phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
|
||||
}
|
||||
|
||||
if(phiS < phi && phi < phiL) {
|
||||
disC = Distance( pos, pos + direction, Q1[i].first, Q1[i].second);
|
||||
|
||||
if( disC < minCathodeDis ){
|
||||
|
|
@ -383,7 +408,7 @@ inline SX3 ANASEN::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){
|
|||
for( int i = 0 ; i < nSX3; i++){
|
||||
|
||||
if(verbose) printf(" %d ", i);
|
||||
std::pair<double, double> frac = Intersect( pos, pos + direction, S1[i].first, S1[i].second, verbose);
|
||||
std::pair<double, double> frac = Intersect( pos, pos + direction, SD[i].first, SD[i].second, verbose);
|
||||
|
||||
|
||||
if( frac.second < 0 || frac.second > 1 ) continue;
|
||||
|
|
@ -394,7 +419,7 @@ inline SX3 ANASEN::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){
|
|||
if(verbose) {
|
||||
printf("reduced distance : %f\n", dis);
|
||||
printf(" %d*", (i+1)%nSX3);
|
||||
Intersect( pos, pos + direction, S1[(i+1)%nSX3].first, S1[(i+1)%nSX3].second, verbose);
|
||||
Intersect( pos, pos + direction, SD[(i+1)%nSX3].first, SD[(i+1)%nSX3].second, verbose);
|
||||
}
|
||||
|
||||
if( TMath::Abs(dis - sx3Radius) > 0.1 ) continue;
|
||||
|
|
@ -407,9 +432,9 @@ inline SX3 ANASEN::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){
|
|||
|
||||
haha.id = zPos > 0 ? i : i + 12;
|
||||
|
||||
haha.zFrac = zPos > 0 ? (zPos - sx3Gap/2 - sx3Length/2)/sx3Length : (zPos - ( - sx3Gap/2 - sx3Length/2) )/sx3Length ;
|
||||
haha.zFrac = zPos > 0 ? (zPos - sx3Gap/2. - sx3Length/2.)/sx3Length : (zPos - ( - sx3Gap/2. - sx3Length/2.) )/sx3Length ;
|
||||
|
||||
haha.chBack = TMath::Floor( haha.zFrac * 4 ) + 8;
|
||||
haha.chBack = TMath::Floor( (haha.zFrac + 0.5) * 4 ) + 8;
|
||||
|
||||
if( verbose) haha.Print();
|
||||
|
||||
|
|
@ -443,7 +468,8 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima
|
|||
int a1 = id.first - 1; if( a1 < 0 ) a1 += nWire;
|
||||
int b1 = id.second - 1; if( b1 < 0 ) b1 += nWire;
|
||||
|
||||
Construct3DModel(a1, id.first+1, b1, id.second+1, false);
|
||||
//Construct3DModel(a1, id.first+1, b1, id.second+1, false);
|
||||
Construct3DModel(id.first, id.first, id.second, id.second, false);
|
||||
|
||||
double theta = direction.Theta() * TMath::RadToDeg();
|
||||
double phi = direction.Phi() * TMath::RadToDeg();
|
||||
|
|
@ -534,13 +560,13 @@ inline TVector3 ANASEN::CalSX3Pos(unsigned short ID, unsigned short chUp, unsign
|
|||
|
||||
if( ID < nSX3 ){ //down
|
||||
|
||||
sa = S1[reducedID].first;
|
||||
sb = S1[reducedID].second;
|
||||
sa = SD[reducedID].first;
|
||||
sb = SD[reducedID].second;
|
||||
|
||||
}else{
|
||||
|
||||
sa = S2[reducedID].first;
|
||||
sb = S2[reducedID].second;
|
||||
sa = SU[reducedID].first;
|
||||
sb = SU[reducedID].second;
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
468
ClassTransfer.h
Normal file
468
ClassTransfer.h
Normal file
|
|
@ -0,0 +1,468 @@
|
|||
#ifndef ClassTransfer_h
|
||||
#define ClassTransfer_h
|
||||
|
||||
#include "TBenchmark.h"
|
||||
#include "TLorentzVector.h"
|
||||
#include "TVector3.h"
|
||||
#include "TMath.h"
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
#include "TRandom.h"
|
||||
#include "TMacro.h"
|
||||
#include "TGraph.h"
|
||||
#include <vector>
|
||||
#include <fstream>
|
||||
|
||||
#include "Isotope.h"
|
||||
|
||||
class ReactionConfig{
|
||||
|
||||
public:
|
||||
|
||||
ReactionConfig(){}
|
||||
~ReactionConfig(){}
|
||||
|
||||
int beamA, beamZ;
|
||||
int targetA, targetZ;
|
||||
int recoilLightA, recoilLightZ;
|
||||
int recoilHeavyA, recoilHeavyZ;
|
||||
|
||||
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 beamX; ///x_offset_of_Beam_in_mm
|
||||
float beamY; ///y_offset_of_Beam_in_mm
|
||||
int numEvents; ///number_of_Event_being_generated
|
||||
bool isTargetScattering; ///isTargetScattering
|
||||
float targetDensity; ///target_density_in_g/cm3
|
||||
float targetThickness; ///targetThickness_in_cm
|
||||
std::string beamStoppingPowerFile; ///stopping_power_for_beam
|
||||
std::string recoilLightStoppingPowerFile; ///stopping_power_for_light_recoil
|
||||
std::string recoilHeavyStoppingPowerFile; ///stopping_power_for_heavy_recoil
|
||||
bool isDecay; ///isDacay
|
||||
int heavyDecayA; ///decayNucleus_A
|
||||
int heavyDecayZ; ///decayNucleus_Z
|
||||
bool isRedo; ///isReDo
|
||||
std::vector<float> beamEx; ///excitation_energy_of_A[MeV]
|
||||
|
||||
|
||||
void SetReaction(int beamA, int beamZ,
|
||||
int targetA, int targetZ,
|
||||
int recoilA, int recoilZ, float beamEnergy_AMeV){
|
||||
this->beamA = beamA;
|
||||
this->beamZ = beamZ;
|
||||
this->targetA = targetA;
|
||||
this->targetZ = targetZ;
|
||||
this->recoilLightA = recoilA;
|
||||
this->recoilLightZ = recoilZ;
|
||||
|
||||
recoilHeavyA = this->beamA + this->targetA - recoilLightA;
|
||||
recoilHeavyZ = this->beamZ + this->targetZ - recoilLightZ;
|
||||
|
||||
}
|
||||
|
||||
void LoadReactionConfig(TMacro * macro){
|
||||
|
||||
if( macro == NULL ) return ;
|
||||
|
||||
int numLine = macro->GetListOfLines()->GetSize();
|
||||
|
||||
for( int i = 0; i < numLine; i ++){
|
||||
|
||||
std::vector<std::string> str = SplitStr(macro->GetListOfLines()->At(i)->GetName(), " ");
|
||||
|
||||
///printf("%d | %s\n", i, str[0].c_str());
|
||||
|
||||
if( str[0].find_first_of("#") == 0 ) break;
|
||||
|
||||
if( i == 0 ) beamA = atoi(str[0].c_str());
|
||||
if( i == 1 ) beamZ = atoi(str[0].c_str());
|
||||
if( i == 2 ) targetA = atoi(str[0].c_str());
|
||||
if( i == 3 ) targetZ = atoi(str[0].c_str());
|
||||
if( i == 4 ) recoilLightA = atoi(str[0].c_str());
|
||||
if( i == 5 ) recoilLightZ = atoi(str[0].c_str());
|
||||
if( i == 6 ) beamEnergy = atof(str[0].c_str());
|
||||
if( i == 7 ) beamEnergySigma = atof(str[0].c_str());
|
||||
if( i == 8 ) beamAngle = atof(str[0].c_str());
|
||||
if( i == 9 ) beamAngleSigma = atof(str[0].c_str());
|
||||
if( i == 10 ) beamX = atof(str[0].c_str());
|
||||
if( i == 11 ) beamY = atof(str[0].c_str());
|
||||
if( i == 12 ) numEvents = atoi(str[0].c_str());
|
||||
if( i == 13 ) {
|
||||
if( str[0].compare("false") == 0 ) isTargetScattering = false;
|
||||
if( str[0].compare("true") == 0 ) isTargetScattering = true;
|
||||
}
|
||||
if( i == 14 ) targetDensity = atof(str[0].c_str());
|
||||
if( i == 15 ) targetThickness = atof(str[0].c_str());
|
||||
if( i == 16 ) beamStoppingPowerFile = str[0];
|
||||
if( i == 17 ) recoilLightStoppingPowerFile = str[0];
|
||||
if( i == 18 ) recoilHeavyStoppingPowerFile = str[0];
|
||||
if( i == 19 ) {
|
||||
if( str[0].compare("false") == 0 ) isDecay = false;
|
||||
if( str[0].compare("true") == 0 ) isDecay = true;
|
||||
}
|
||||
if( i == 20 ) heavyDecayA = atoi(str[0].c_str());
|
||||
if( i == 21 ) heavyDecayZ = atoi(str[0].c_str());
|
||||
if( i == 22 ) {
|
||||
if( str[0].compare("false") == 0 ) isRedo = false;
|
||||
if( str[0].compare("true" ) == 0 ) isRedo = true;
|
||||
}
|
||||
|
||||
if( i >= 23) {
|
||||
beamEx.push_back( atof(str[0].c_str()) );
|
||||
}
|
||||
}
|
||||
|
||||
recoilHeavyA = beamA + targetA - recoilLightA;
|
||||
recoilHeavyZ = beamZ + targetZ - recoilLightZ;
|
||||
|
||||
}
|
||||
|
||||
void PrintReactionConfig(){
|
||||
|
||||
printf("=====================================================\n");
|
||||
printf(" beam : A = %3d, Z = %2d \n", beamA, beamZ);
|
||||
printf(" target : A = %3d, Z = %2d \n", targetA, targetZ);
|
||||
printf(" light : A = %3d, Z = %2d \n", recoilLightA, recoilLightZ);
|
||||
|
||||
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(" offset : (x,y) = (%.2f, %.2f) mmm \n", beamX, beamY);
|
||||
|
||||
printf("##### number of Simulation Events : %d \n", numEvents);
|
||||
|
||||
printf(" is target scattering : %s \n", isTargetScattering ? "Yes" : "No");
|
||||
|
||||
if(isTargetScattering){
|
||||
printf(" target density : %.f g/cm3\n", targetDensity);
|
||||
printf(" thickness : %.f cm\n", targetThickness);
|
||||
printf(" beam stopping file : %s \n", beamStoppingPowerFile.c_str());
|
||||
printf(" recoil light stopping file : %s \n", recoilLightStoppingPowerFile.c_str());
|
||||
printf(" recoil heavy stopping file : %s \n", recoilHeavyStoppingPowerFile.c_str());
|
||||
}
|
||||
|
||||
printf(" is simulate decay : %s \n", isDecay ? "Yes" : "No");
|
||||
if( isDecay ){
|
||||
printf(" heavy decay : A = %d, Z = %d \n", heavyDecayA, heavyDecayZ);
|
||||
}
|
||||
printf(" is Redo until hit array : %s \n", isRedo ? "Yes" : "No");
|
||||
|
||||
printf(" beam Ex : %.2f MeV \n", beamEx[0]);
|
||||
for( int i = 1; i < (int) beamEx.size(); i++){
|
||||
printf(" %.2f MeV \n", beamEx[i]);
|
||||
}
|
||||
|
||||
printf("=====================================================\n");
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
//=======================================================
|
||||
//#######################################################
|
||||
// Class for Transfer Reaction
|
||||
// reaction notation A(a,b)B
|
||||
// A = incident particle
|
||||
// a = target
|
||||
// b = light scattered particle
|
||||
// B = heavy scattered particle
|
||||
//=======================================================
|
||||
class TransferReaction {
|
||||
public:
|
||||
TransferReaction();
|
||||
~TransferReaction();
|
||||
|
||||
void SetA(int A, int Z, double Ex);
|
||||
void Seta(int A, int Z);
|
||||
void Setb(int A, int Z);
|
||||
void SetB(int A, int Z);
|
||||
void SetIncidentEnergyAngle(double KEA, double theta, double phi);
|
||||
void SetExA(double Ex);
|
||||
void SetExB(double Ex);
|
||||
void SetReactionFromFile(string settingFile);
|
||||
|
||||
TString GetReactionName();
|
||||
TString GetReactionName_Latex();
|
||||
|
||||
ReactionConfig GetRectionConfig() { return reaction;}
|
||||
|
||||
double GetMass_A(){return mA + ExA;}
|
||||
double GetMass_a(){return ma;}
|
||||
double GetMass_b(){return mb;}
|
||||
double GetMass_B(){return mB + ExB;}
|
||||
|
||||
double GetCMTotalKE() {return Etot - mA - ma;}
|
||||
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;}
|
||||
|
||||
void CalReactionConstant();
|
||||
|
||||
TLorentzVector * Event(double thetaCM, double phiCM);
|
||||
|
||||
double GetEx(){return Ex;}
|
||||
double GetThetaCM(){return thetaCM;}
|
||||
|
||||
double GetMomentumbCM() {return p;}
|
||||
double GetReactionBeta() {return beta;}
|
||||
double GetReactionGamma() {return gamma;}
|
||||
double GetCMTotalEnergy() {return Etot;}
|
||||
|
||||
private:
|
||||
|
||||
ReactionConfig reaction;
|
||||
|
||||
string nameA, namea, nameb, nameB;
|
||||
double thetaIN, phiIN;
|
||||
double mA, ma, mb, mB;
|
||||
|
||||
double TA, T; // TA = KE of A pre u, T = total energy
|
||||
double ExA, ExB;
|
||||
double Ex, thetaCM; //calculated Ex using inverse mapping from e and z to thetaCM
|
||||
|
||||
bool isReady;
|
||||
bool isBSet;
|
||||
|
||||
double k; // CM Boost momentum
|
||||
double beta, gamma; //CM boost beta
|
||||
double Etot;
|
||||
double p; // CM frame momentum of b, B
|
||||
|
||||
TLorentzVector PA, Pa, Pb, PB;
|
||||
|
||||
TString format(TString name);
|
||||
|
||||
};
|
||||
|
||||
TransferReaction::TransferReaction(){
|
||||
|
||||
thetaIN = 0.;
|
||||
phiIN = 0.;
|
||||
SetA(12, 6, 0);
|
||||
Seta(2,1);
|
||||
Setb(1,1);
|
||||
SetB(13,6);
|
||||
TA = 6;
|
||||
T = TA * reaction.beamA;
|
||||
|
||||
ExA = 0;
|
||||
ExB = 0;
|
||||
|
||||
Ex = TMath::QuietNaN();
|
||||
thetaCM = TMath::QuietNaN();
|
||||
|
||||
CalReactionConstant();
|
||||
|
||||
TLorentzVector temp (0,0,0,0);
|
||||
PA = temp;
|
||||
Pa = temp;
|
||||
Pb = temp;
|
||||
PB = temp;
|
||||
|
||||
}
|
||||
|
||||
TransferReaction::~TransferReaction(){
|
||||
|
||||
}
|
||||
|
||||
void TransferReaction::SetA(int A, int Z, double Ex = 0){
|
||||
Isotope temp (A, Z);
|
||||
mA = temp.Mass;
|
||||
reaction.beamA = A;
|
||||
reaction.beamZ = Z;
|
||||
ExA = Ex;
|
||||
nameA = temp.Name;
|
||||
isReady = false;
|
||||
isBSet = true;
|
||||
|
||||
}
|
||||
|
||||
void TransferReaction::Seta(int A, int Z){
|
||||
Isotope temp (A, Z);
|
||||
ma = temp.Mass;
|
||||
reaction.targetA = A;
|
||||
reaction.targetZ = Z;
|
||||
namea = temp.Name;
|
||||
isReady = false;
|
||||
isBSet = false;
|
||||
}
|
||||
|
||||
void TransferReaction::Setb(int A, int Z){
|
||||
Isotope temp (A, Z);
|
||||
mb = temp.Mass;
|
||||
reaction.recoilLightA = A;
|
||||
reaction.recoilLightZ = Z;
|
||||
nameb = temp.Name;
|
||||
isReady = false;
|
||||
isBSet = false;
|
||||
}
|
||||
void TransferReaction::SetB(int A, int Z){
|
||||
Isotope temp (A, Z);
|
||||
mB = temp.Mass;
|
||||
reaction.recoilHeavyA = A;
|
||||
reaction.recoilHeavyZ = Z;
|
||||
nameB = temp.Name;
|
||||
isReady = false;
|
||||
isBSet = true;
|
||||
}
|
||||
|
||||
void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){
|
||||
this->TA = KEA;
|
||||
this->T = TA * reaction.beamA;
|
||||
this->thetaIN = theta;
|
||||
this->phiIN = phi;
|
||||
isReady = false;
|
||||
}
|
||||
|
||||
void TransferReaction::SetExA(double Ex){
|
||||
this->ExA = Ex;
|
||||
isReady = false;
|
||||
}
|
||||
|
||||
void TransferReaction::SetExB(double Ex){
|
||||
this->ExB = Ex;
|
||||
isReady = false;
|
||||
}
|
||||
|
||||
void TransferReaction::SetReactionFromFile(string settingFile){
|
||||
|
||||
TMacro * haha = new TMacro();
|
||||
if( haha->ReadFile(settingFile.c_str()) > 0 ) {
|
||||
reaction.LoadReactionConfig(haha);
|
||||
|
||||
SetA(reaction.beamA, reaction.beamZ);
|
||||
Seta(reaction.targetA, reaction.targetZ);
|
||||
Setb(reaction.recoilLightA, reaction.recoilLightZ);
|
||||
SetB(reaction.recoilHeavyA, reaction.recoilHeavyZ);
|
||||
|
||||
SetIncidentEnergyAngle(reaction.beamEnergy, 0, 0);
|
||||
CalReactionConstant();
|
||||
}else{
|
||||
|
||||
printf("cannot read file %s.\n", settingFile.c_str());
|
||||
isReady = false;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
TString TransferReaction::GetReactionName(){
|
||||
TString rName;
|
||||
rName.Form("%s(%s,%s)%s", nameA.c_str(), namea.c_str(), nameb.c_str(), nameB.c_str());
|
||||
return rName;
|
||||
}
|
||||
|
||||
TString TransferReaction::format(TString name){
|
||||
if( name.IsAlpha() ) return name;
|
||||
int len = name.Length();
|
||||
TString temp = name;
|
||||
TString temp2 = name;
|
||||
if( temp.Remove(0, len-2).IsAlpha()){
|
||||
temp2.Remove(len-2);
|
||||
}else{
|
||||
temp = name;
|
||||
temp.Remove(0, len-1);
|
||||
temp2.Remove(len-1);
|
||||
}
|
||||
return "^{"+temp2+"}"+temp;
|
||||
}
|
||||
|
||||
TString TransferReaction::GetReactionName_Latex(){
|
||||
TString rName;
|
||||
rName.Form("%s(%s,%s)%s", format(nameA).Data(), format(namea).Data(), format(nameb).Data(), format(nameB).Data());
|
||||
return rName;
|
||||
}
|
||||
|
||||
void TransferReaction::CalReactionConstant(){
|
||||
if( !isBSet){
|
||||
reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA;
|
||||
reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ;
|
||||
Isotope temp (reaction.recoilHeavyA, reaction.recoilHeavyZ);
|
||||
mB = temp.Mass;
|
||||
isBSet = true;
|
||||
}
|
||||
|
||||
k = TMath::Sqrt(TMath::Power(mA + ExA + T, 2) - (mA + ExA) * (mA + ExA));
|
||||
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);
|
||||
|
||||
isReady = true;
|
||||
}
|
||||
|
||||
TLorentzVector * TransferReaction::Event(double thetaCM, double phiCM)
|
||||
{
|
||||
if( isReady == false ){
|
||||
CalReactionConstant();
|
||||
}
|
||||
|
||||
//TLorentzVector Pa(0, 0, 0, ma);
|
||||
|
||||
//---- to CM frame
|
||||
TLorentzVector Pc = PA + Pa;
|
||||
TVector3 b = Pc.BoostVector();
|
||||
|
||||
TVector3 vb(0,0,0);
|
||||
|
||||
if( b.Mag() > 0 ){
|
||||
TVector3 v0 (0,0,0);
|
||||
TVector3 nb = v0 - b;
|
||||
|
||||
TLorentzVector PAc = PA;
|
||||
PAc.Boost(nb);
|
||||
TVector3 vA = PAc.Vect();
|
||||
|
||||
TLorentzVector Pac = Pa;
|
||||
Pac.Boost(nb);
|
||||
TVector3 va = Pac.Vect();
|
||||
|
||||
//--- construct vb
|
||||
vb = va;
|
||||
vb.SetMag(p);
|
||||
|
||||
TVector3 ub = vb.Orthogonal();
|
||||
vb.Rotate(thetaCM, ub);
|
||||
vb.Rotate(phiCM + TMath::PiOver2(), va); // somehow, the calculation turn the vector 90 degree.
|
||||
//vb.Rotate(phiCM , va); // somehow, the calculation turn the vector 90 degree.
|
||||
}
|
||||
|
||||
//--- from Pb
|
||||
TLorentzVector Pbc;
|
||||
Pbc.SetVectM(vb, mb);
|
||||
|
||||
//--- from PB
|
||||
TLorentzVector PBc;
|
||||
//PBc.SetVectM(vB, mB + ExB);
|
||||
PBc.SetVectM(-vb, mB + ExB);
|
||||
|
||||
//---- to Lab Frame
|
||||
TLorentzVector Pb = Pbc;
|
||||
Pb.Boost(b);
|
||||
TLorentzVector PB = PBc;
|
||||
PB.Boost(b);
|
||||
|
||||
TLorentzVector * output = new TLorentzVector[4];
|
||||
output[0] = PA;
|
||||
output[1] = Pa;
|
||||
output[2] = Pb;
|
||||
output[3] = PB;
|
||||
|
||||
this->Pb = Pb;
|
||||
this->PB = PB;
|
||||
|
||||
|
||||
return output;
|
||||
}
|
||||
|
||||
#endif
|
||||
523
Isotope.h
Normal file
523
Isotope.h
Normal file
|
|
@ -0,0 +1,523 @@
|
|||
/***********************************************************************
|
||||
*
|
||||
* This is Isotope.h, To extract the isotope mass from massXX.txt
|
||||
*
|
||||
*-------------------------------------------------------
|
||||
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
|
||||
* email: goluckyryan@gmail.com
|
||||
* ********************************************************************/
|
||||
|
||||
|
||||
#ifndef ISOTOPE_H
|
||||
#define ISOTOPE_H
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <string>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include "constant.h" // amu
|
||||
#include <stdlib.h> //atoi
|
||||
#include <algorithm>
|
||||
using namespace std;
|
||||
|
||||
string massData="mass20.txt";
|
||||
|
||||
// about the mass**.txt
|
||||
// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu
|
||||
// mass excess uncertaintly
|
||||
// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons
|
||||
// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons
|
||||
|
||||
class Isotope {
|
||||
public:
|
||||
int A, Z;
|
||||
double Mass, MassError, BEA;
|
||||
string Name, Symbol;
|
||||
string dataSource;
|
||||
|
||||
Isotope(){ dataSource = massData; };
|
||||
Isotope(int a, int z){ dataSource = massData; SetIso(a,z); };
|
||||
Isotope(string name){ dataSource = massData; SetIsoByName(name); };
|
||||
|
||||
void SetIso(int a, int z);
|
||||
void SetIsoByName(string name);
|
||||
|
||||
double CalSp(int Np, int Nn); // this for the Np-proton, Nn-neutron removal
|
||||
double CalSp2(int a, int z); // this is for (a,z) nucleus removal
|
||||
|
||||
double CalBeta(double T){
|
||||
double Etot = Mass + T;
|
||||
double gamma = 1 + T/Mass;
|
||||
double beta = sqrt(1 - 1 / gamma / gamma ) ;
|
||||
return beta;
|
||||
}
|
||||
|
||||
void Print();
|
||||
void ListShell();
|
||||
|
||||
private:
|
||||
void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol;
|
||||
void FindMassByName(string name); // give Z, mass, massError, BEA;
|
||||
|
||||
int TwoJ(int nShell);
|
||||
string Orbital(int nShell);
|
||||
int magic(int i){
|
||||
switch (i){
|
||||
case 0: return 2; break;
|
||||
case 1: return 8; break;
|
||||
case 2: return 20; break;
|
||||
case 3: return 28; break;
|
||||
case 4: return 40; break;
|
||||
case 5: return 50; break;
|
||||
case 6: return 82; break;
|
||||
case 7: return 128; break;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
int magicShellID(int i){
|
||||
switch (i){
|
||||
case 0: return 0; break;
|
||||
case 1: return 2; break;
|
||||
case 2: return 5; break;
|
||||
case 3: return 6; break;
|
||||
case 4: return 9; break;
|
||||
case 5: return 10; break;
|
||||
case 6: return 15; break;
|
||||
case 7: return 21; break;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
int fileStartLine;
|
||||
int fileEndLine;
|
||||
int lineMass050_099;
|
||||
int lineMass100_149;
|
||||
int lineMass150_199;
|
||||
int lineMass200;
|
||||
|
||||
|
||||
void setFileLines(){
|
||||
fileStartLine = 37;
|
||||
fileEndLine = 3594;
|
||||
|
||||
lineMass050_099 = 466;
|
||||
lineMass100_149 = 1160;
|
||||
lineMass150_199 = 1994;
|
||||
lineMass200 = 2774;
|
||||
}
|
||||
|
||||
char * heliosPath;
|
||||
bool isFindOnce;
|
||||
|
||||
|
||||
};
|
||||
|
||||
void Isotope::SetIso(int a, int z){
|
||||
this->A = a;
|
||||
this->Z = z;
|
||||
FindMassByAZ(a,z);
|
||||
}
|
||||
|
||||
void Isotope::SetIsoByName(string name){
|
||||
FindMassByName(name);
|
||||
}
|
||||
|
||||
void Isotope::FindMassByAZ(int A, int Z){
|
||||
string line;
|
||||
int lineNum=0;
|
||||
int list_A, list_Z;
|
||||
|
||||
ifstream myfile;
|
||||
int flag=0;
|
||||
|
||||
setFileLines();
|
||||
|
||||
int numLineStart = fileStartLine;
|
||||
int numLineEnd = fileEndLine;
|
||||
|
||||
if ( A >= 50 && A < 100) numLineStart = lineMass050_099;
|
||||
if ( A >=100 && A < 150) numLineStart = lineMass100_149;
|
||||
if ( A >=150 && A < 200) numLineStart = lineMass150_199;
|
||||
if ( A >=200 ) numLineStart = lineMass200;
|
||||
|
||||
myfile.open(dataSource.c_str());
|
||||
|
||||
if (myfile.is_open()) {
|
||||
while (/*! myfile.eof() &&*/ flag == 0 && lineNum <numLineEnd){
|
||||
lineNum ++ ;
|
||||
//printf("%3d ",lineNum);
|
||||
getline (myfile,line);
|
||||
|
||||
if (lineNum >= numLineStart ){
|
||||
list_Z = atoi((line.substr(10,5)).c_str());
|
||||
list_A = atoi((line.substr(15,5)).c_str());
|
||||
|
||||
if ( A == list_A && Z == list_Z) {
|
||||
this->BEA = atof((line.substr(54,11)).c_str());
|
||||
this->Mass = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A;
|
||||
this->MassError = atof((line.substr(65,7)).c_str());
|
||||
string str = line.substr(20,2);
|
||||
str.erase(remove(str.begin(), str.end(), ' '), str.end());
|
||||
this->Symbol = str;
|
||||
|
||||
ostringstream ss;
|
||||
ss << A << this->Symbol;
|
||||
this->Name = ss.str();
|
||||
flag = 1;
|
||||
}else if ( list_A > A) {
|
||||
this->BEA = -404;
|
||||
this->Mass = -404;
|
||||
this->MassError = -404;
|
||||
this->Symbol = "non";
|
||||
this->Name = "non";
|
||||
break;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if( this->Name == "1H" ) this->Name = "p";
|
||||
if( this->Name == "2H" ) this->Name = "d";
|
||||
if( this->Name == "3H" ) this->Name = "t";
|
||||
if( this->Name == "4He" ) this->Name = "a";
|
||||
|
||||
myfile.close();
|
||||
}else {
|
||||
printf("Unable to open %s\n", dataSource.c_str());
|
||||
}
|
||||
}
|
||||
|
||||
void Isotope::FindMassByName(string name){
|
||||
|
||||
// done seperate the Mass number and the name
|
||||
if( name == "n" ) {
|
||||
this->Name = "1n";
|
||||
this->BEA = 0;
|
||||
this->Mass = mn;
|
||||
this->MassError = 0;
|
||||
this->Name = "n";
|
||||
this->A = 1;
|
||||
this->Z = 0;
|
||||
return;
|
||||
}
|
||||
if( name == "p" ) name = "1H";
|
||||
if( name == "d" ) name = "2H";
|
||||
if( name == "t" ) name = "3H";
|
||||
if( name == "a" ) name = "4He";
|
||||
|
||||
string temp = name;
|
||||
int lastDigit = 0;
|
||||
|
||||
for(int i=0; temp[i]; i++){
|
||||
if(temp[i] == '0') lastDigit = i;
|
||||
if(temp[i] == '1') lastDigit = i;
|
||||
if(temp[i] == '2') lastDigit = i;
|
||||
if(temp[i] == '3') lastDigit = i;
|
||||
if(temp[i] == '4') lastDigit = i;
|
||||
if(temp[i] == '5') lastDigit = i;
|
||||
if(temp[i] == '6') lastDigit = i;
|
||||
if(temp[i] == '7') lastDigit = i;
|
||||
if(temp[i] == '8') lastDigit = i;
|
||||
if(temp[i] == '9') lastDigit = i;
|
||||
}
|
||||
|
||||
this->Symbol = temp.erase(0, lastDigit +1);
|
||||
//check is Symbol is 2 charaters, if not, add " " at the end
|
||||
if( this->Symbol.length() == 1 ){
|
||||
this->Symbol = this->Symbol + " ";
|
||||
}
|
||||
|
||||
|
||||
temp = name;
|
||||
int len = temp.length();
|
||||
temp = temp.erase(lastDigit+1, len);
|
||||
|
||||
this->A = atoi(temp.c_str());
|
||||
//printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A);
|
||||
|
||||
// find the nucleus in the data
|
||||
string line;
|
||||
int lineNum=0;
|
||||
int list_A;
|
||||
string list_symbol;
|
||||
|
||||
ifstream myfile;
|
||||
int flag=0;
|
||||
|
||||
setFileLines();
|
||||
|
||||
int numLineStart = fileStartLine;
|
||||
int numLineEnd = fileEndLine;
|
||||
|
||||
if ( A >= 50 && A < 100) numLineStart = lineMass050_099;
|
||||
if ( A >=100 && A < 150) numLineStart = lineMass100_149;
|
||||
if ( A >=150 && A < 200) numLineStart = lineMass150_199;
|
||||
if ( A >=200 ) numLineStart = lineMass200;
|
||||
|
||||
myfile.open(dataSource.c_str());
|
||||
|
||||
if (myfile.is_open()) {
|
||||
while (/*! myfile.eof() &&*/ flag == 0 && lineNum <numLineEnd){
|
||||
lineNum ++ ;
|
||||
//printf("%3d ",lineNum);
|
||||
getline (myfile,line);
|
||||
|
||||
if (lineNum >= numLineStart ){
|
||||
list_symbol = line.substr(20,2);
|
||||
list_A = atoi((line.substr(15,5)).c_str());
|
||||
|
||||
//printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str());
|
||||
|
||||
if ( this->A == list_A && this->Symbol == list_symbol) {
|
||||
this->Z = atoi((line.substr(10,5)).c_str());
|
||||
this->BEA = atof((line.substr(54,11)).c_str());
|
||||
this->Mass = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A;
|
||||
this->MassError = atof((line.substr(65,7)).c_str());
|
||||
|
||||
string str = line.substr(20,2);
|
||||
str.erase(remove(str.begin(), str.end(), ' '), str.end());
|
||||
this->Symbol = str;
|
||||
|
||||
ostringstream ss;
|
||||
ss << this->A << this->Symbol;
|
||||
this->Name = ss.str();
|
||||
flag = 1;
|
||||
}else if ( list_A > this->A) {
|
||||
this->BEA = -404;
|
||||
this->Mass = -404;
|
||||
this->MassError = -404;
|
||||
this->Symbol = "non";
|
||||
this->Name = "non";
|
||||
break;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
myfile.close();
|
||||
}else {
|
||||
printf("Unable to open %s\n", dataSource.c_str());
|
||||
}
|
||||
}
|
||||
|
||||
double Isotope::CalSp(int Np, int Nn){
|
||||
Isotope nucleusD(A - Np - Nn, Z - Np);
|
||||
|
||||
if( nucleusD.Mass != -404){
|
||||
return nucleusD.Mass + Nn*mn + Np*mp - this->Mass;
|
||||
}else{
|
||||
return -404;
|
||||
}
|
||||
}
|
||||
|
||||
double Isotope::CalSp2(int a, int z){
|
||||
Isotope nucleusD(A - a , Z - z);
|
||||
Isotope nucleusS(a,z);
|
||||
|
||||
if( nucleusD.Mass != -404 && nucleusS.Mass != -404){
|
||||
return nucleusD.Mass + nucleusS.Mass - this->Mass;
|
||||
}else{
|
||||
return -404;
|
||||
}
|
||||
}
|
||||
|
||||
int Isotope::TwoJ(int nShell){
|
||||
|
||||
switch(nShell){
|
||||
case 0: return 1; break; // 0s1/2
|
||||
case 1: return 3; break; // 0p3/2
|
||||
case 2: return 1; break; // 0p1/2 -- 8
|
||||
case 3: return 5; break; // 0d5/2
|
||||
case 4: return 1; break; // 1s1/2
|
||||
case 5: return 3; break; // 0d3/2 -- 20
|
||||
case 6: return 7; break; // 0f7/2 -- 28
|
||||
case 7: return 3; break; // 1p3/2
|
||||
case 8: return 1; break; // 1p1/2
|
||||
case 9: return 5; break; // 0f5/2 -- 40
|
||||
case 10: return 9; break; // 0g9/2 -- 50
|
||||
case 11: return 7; break; // 0g7/2
|
||||
case 12: return 5; break; // 1d5/2
|
||||
case 13: return 11; break; // 0h11/2
|
||||
case 14: return 3; break; // 1d3/2
|
||||
case 15: return 1; break; // 2s1/2 -- 82
|
||||
case 16: return 9; break; // 0h9/2
|
||||
case 17: return 7; break; // 1f7/2
|
||||
case 18: return 13; break; // 0i13/2
|
||||
case 19: return 3; break; // 2p3/2
|
||||
case 20: return 5; break; // 1f5/2
|
||||
case 21: return 1; break; // 1p1/2 -- 126
|
||||
case 22: return 9; break; // 1g9/2
|
||||
case 23: return 11; break; // 0i11/2
|
||||
case 24: return 15; break; // 0j15/2
|
||||
case 25: return 5; break; // 2d5/2
|
||||
case 26: return 1; break; // 3s1/2
|
||||
case 27: return 3; break; // 2d3/2
|
||||
case 28: return 7; break; // 1g7/2
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
string Isotope::Orbital(int nShell){
|
||||
|
||||
switch(nShell){
|
||||
case 0: return "0s1 "; break; //
|
||||
case 1: return "0p3 "; break; //
|
||||
case 2: return "0p1 "; break; //-- 8
|
||||
case 3: return "0d5 "; break; //
|
||||
case 4: return "1s1 "; break; //
|
||||
case 5: return "0d3 "; break; //-- 20
|
||||
case 6: return "0f7 "; break; //-- 28
|
||||
case 7: return "1p3 "; break; //
|
||||
case 8: return "1p1 "; break; //
|
||||
case 9: return "0f5 "; break; //-- 40
|
||||
case 10: return "0g9 "; break; //-- 50
|
||||
case 11: return "0g7 "; break; //
|
||||
case 12: return "1d5 "; break; //
|
||||
case 13: return "0h11"; break; //
|
||||
case 14: return "1d3 "; break; //
|
||||
case 15: return "2s1 "; break; //-- 82
|
||||
case 16: return "0h9 "; break; //
|
||||
case 17: return "1f7 "; break; //
|
||||
case 18: return "0i13"; break; //
|
||||
case 19: return "2p3 "; break; //
|
||||
case 20: return "1f5 "; break; //
|
||||
case 21: return "1p1 "; break; //-- 126
|
||||
case 22: return "1g9 "; break; //
|
||||
case 23: return "0i11"; break; //
|
||||
case 24: return "0j15"; break; //
|
||||
case 25: return "2d5 "; break; //
|
||||
case 26: return "3s1 "; break; //
|
||||
case 27: return "2d3 "; break; //
|
||||
case 28: return "1g7 "; break; //
|
||||
}
|
||||
|
||||
return "nan";
|
||||
}
|
||||
|
||||
void Isotope::ListShell(){
|
||||
|
||||
if( Mass < 0 ) return;
|
||||
|
||||
int n = A-Z;
|
||||
int p = Z;
|
||||
|
||||
int k = min(n,p);
|
||||
int nMagic = 0;
|
||||
for( int i = 0; i < 7; i++){
|
||||
if( magic(i) < k && k <= magic(i+1) ){
|
||||
nMagic = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
int coreShell = magicShellID(nMagic-1);
|
||||
int coreZ1 = magic(nMagic-1);
|
||||
int coreZ2 = magic(nMagic);
|
||||
|
||||
Isotope core1( 2*coreZ1, coreZ1);
|
||||
Isotope core2( 2*coreZ2, coreZ2);
|
||||
|
||||
printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str());
|
||||
printf(" || ");
|
||||
int t = max(n,p);
|
||||
int nShell = 0;
|
||||
do{
|
||||
int occ = TwoJ(nShell)+1;
|
||||
if( nShell > coreShell) {
|
||||
printf("%4s", Orbital(nShell).c_str());
|
||||
if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){
|
||||
printf("|");
|
||||
}else{
|
||||
printf(",");
|
||||
}
|
||||
}
|
||||
t = t - occ;
|
||||
nShell++;
|
||||
}while( t > 0 && nShell < 29);
|
||||
for( int i = 1; i <= 6; i++){
|
||||
if (nShell < 28) {
|
||||
printf("%4s,", Orbital(nShell).c_str());
|
||||
}else if( nShell == 28 ) {
|
||||
printf("%4s", Orbital(nShell).c_str());
|
||||
}
|
||||
nShell ++;
|
||||
}
|
||||
if (nShell < 29) printf("%4s", Orbital(nShell).c_str());
|
||||
printf("\n");
|
||||
|
||||
|
||||
printf(" Z = %3d || ", p);
|
||||
nShell = 0;
|
||||
do{
|
||||
int occ = TwoJ(nShell)+1;
|
||||
if( nShell > coreShell ){
|
||||
if( p > occ ) {
|
||||
printf("%-4d", occ);
|
||||
if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){
|
||||
printf("|");
|
||||
}else{
|
||||
printf(",");
|
||||
}
|
||||
}else{
|
||||
printf("%-4d", p);
|
||||
}
|
||||
}
|
||||
p = p - occ;
|
||||
nShell++;
|
||||
}while( p > 0 && nShell < 29);
|
||||
printf("\n");
|
||||
|
||||
printf(" N = %3d || ", n);
|
||||
nShell = 0;
|
||||
do{
|
||||
int occ = TwoJ(nShell)+1;
|
||||
if ( nShell > coreShell ){
|
||||
if( n > occ ) {
|
||||
printf("%-4d", occ);
|
||||
if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9 || nShell == 10 || nShell == 15 || nShell == 21){
|
||||
printf("|");
|
||||
}else{
|
||||
printf(",");
|
||||
}
|
||||
}else{
|
||||
printf("%-4d", n);
|
||||
}
|
||||
}
|
||||
n = n - occ;
|
||||
nShell++;
|
||||
}while( n > 0 && nShell < 29);
|
||||
printf("\n");
|
||||
|
||||
printf("------------------ \n");
|
||||
}
|
||||
|
||||
|
||||
|
||||
void Isotope::Print(){
|
||||
|
||||
if (Mass > 0){
|
||||
|
||||
dataSource = massData;
|
||||
|
||||
printf(" using mass data : %s \n", dataSource.c_str());
|
||||
printf(" mass of \e[47m\e[31m%s\e[m nucleus (Z,A)=(%3d,%3d) is \e[47m\e[31m%12.5f\e[m MeV, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.);
|
||||
printf(" total BE : %12.5f MeV\n",BEA*A/1000.);
|
||||
printf(" mass in amu : %12.5f u\n",Mass/amu);
|
||||
printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu);
|
||||
printf("-------------- Seperation energy \n");
|
||||
printf(" S1p: %8.4f| S1n: %8.4f| S(2H ): %8.4f| S1p1n : %8.4f\n", CalSp(1, 0), CalSp(0, 1), CalSp2(2, 1), CalSp(1, 1));
|
||||
printf(" S2p: %8.4f| S2n: %8.4f| S(3He): %8.4f| S(3H) : %8.4f\n", CalSp(2, 0), CalSp(0, 2), CalSp2(3, 2), CalSp2(3, 1));
|
||||
printf(" S3p: %8.4f| S3n: %8.4f| S(4He): %8.4f|\n", CalSp(3, 0), CalSp(0, 3), CalSp2(4, 2));
|
||||
printf(" S4p: %8.4f| S4n: %8.4f| \n", CalSp(4, 0), CalSp(0, 4));
|
||||
|
||||
}else{
|
||||
printf("Error %6.0f, no nucleus with (Z,A) = (%3d,%3d). \n", Mass, Z, A);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
208
anasenMS.C
208
anasenMS.C
|
|
@ -1,5 +1,3 @@
|
|||
#include "ClassAnasen.h"
|
||||
|
||||
#include "TRandom.h"
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
|
|
@ -7,26 +5,218 @@
|
|||
#include "TH2.h"
|
||||
#include "TStyle.h"
|
||||
#include "TCanvas.h"
|
||||
#include "TBenchmark.h"
|
||||
|
||||
#include "ClassTransfer.h"
|
||||
#include "ClassAnasen.h"
|
||||
|
||||
//======== Gerneate light particle based on reaction
|
||||
// find out the CalTrack and the real track
|
||||
// find out the Q-value uncertaintly
|
||||
|
||||
std::pair<TVector3, TVector3> Transfer(){
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
void anasenMS(){
|
||||
|
||||
const int numEvent = 100000;
|
||||
|
||||
//Reaction
|
||||
TransferReaction transfer;
|
||||
|
||||
transfer.SetA(12, 6, 0);
|
||||
transfer.SetIncidentEnergyAngle(10, 0, 0);
|
||||
|
||||
transfer.Seta( 2, 1);
|
||||
transfer.Setb( 1, 1);
|
||||
|
||||
std::vector<float> ExAList = {0};
|
||||
std::vector<float> ExList = {0, 1, 2};
|
||||
|
||||
double vertexXRange[2] = { 0, 0};
|
||||
double vertexYRange[2] = { 0, 0};
|
||||
double vertexZRange[2] = { 0, 0};
|
||||
|
||||
//###################################################
|
||||
transfer.CalReactionConstant();
|
||||
|
||||
int nExA = ExAList.size();
|
||||
int nEx = ExList.size();
|
||||
|
||||
ANASEN anasen;
|
||||
|
||||
TString saveFileName = "msAnasen.root";
|
||||
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
|
||||
TFile * saveFile = new TFile(saveFileName, "recreate");
|
||||
TTree * tree = new TTree("tree", "tree");
|
||||
|
||||
double KEA;
|
||||
tree->Branch("beamKEA", &KEA, "beamKEA/D");
|
||||
|
||||
double thetaCM, phiCM;
|
||||
tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
|
||||
tree->Branch("phiCM", &phiCM, "phiCM/D");
|
||||
|
||||
double thetab, phib, Tb;
|
||||
double thetaB, phiB, TB;
|
||||
tree->Branch("thetab", &thetab, "thetab/D");
|
||||
tree->Branch("phib", &phib, "phib/D");
|
||||
tree->Branch("Tb", &Tb, "Tb/D");
|
||||
tree->Branch("thetaB", &thetaB, "thetaB/D");
|
||||
tree->Branch("phiB", &phiB, "phiB/D");
|
||||
tree->Branch("TB", &TB, "TB/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");
|
||||
tree->Branch("Ex", &Ex, "Ex/D");
|
||||
|
||||
double vertexX, vertexY, vertexZ;
|
||||
tree->Branch("vX", &vertexX, "VertexX/D");
|
||||
tree->Branch("vY", &vertexY, "VertexY/D");
|
||||
tree->Branch("vZ", &vertexZ, "VertexZ/D");
|
||||
|
||||
double sx3X, sx3Y, sx3Z;
|
||||
tree->Branch("sx3X", &sx3X, "sx3X/D");
|
||||
tree->Branch("sx3Y", &sx3Y, "sx3Y/D");
|
||||
tree->Branch("sx3Z", &sx3Z, "sx3Z/D");
|
||||
|
||||
int anodeID, cathodeID;
|
||||
tree->Branch("aID", &anodeID, "anodeID/I");
|
||||
tree->Branch("cID", &cathodeID, "cathodeID/I");
|
||||
|
||||
int sx3ID, sx3Up, sx3Down, sx3Back;
|
||||
double sx3ZFrac;
|
||||
tree->Branch("sx3ID", &sx3ID, "sx3ID/I");
|
||||
tree->Branch("sx3Up", &sx3Up, "sx3Up/I");
|
||||
tree->Branch("sx3Down", &sx3Down, "sx3Down/I");
|
||||
tree->Branch("sx3Back", &sx3Back, "sx3Back/I");
|
||||
tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
|
||||
|
||||
double reTheta, rePhi;
|
||||
tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
|
||||
tree->Branch("rePhi", &rePhi, "reconstucted_phi/D");
|
||||
|
||||
|
||||
//========timer
|
||||
TBenchmark clock;
|
||||
bool shown ;
|
||||
clock.Reset();
|
||||
clock.Start("timer");
|
||||
shown = false;
|
||||
|
||||
//================================= Calculate event
|
||||
for( int i = 0; i < numEvent ; i++){
|
||||
|
||||
ExAID = gRandom->Integer(nExA);
|
||||
ExA = ExAList[ExAID];
|
||||
transfer.SetExA(ExA);
|
||||
|
||||
ExID = gRandom->Integer(nEx);
|
||||
Ex = ExList[ExID];
|
||||
transfer.SetExB(Ex);
|
||||
|
||||
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
|
||||
phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
|
||||
|
||||
//==== Calculate reaction
|
||||
TLorentzVector * output = transfer.Event(thetaCM, phiCM);
|
||||
TLorentzVector Pb = output[2];
|
||||
TLorentzVector PB = output[3];
|
||||
|
||||
thetab = Pb.Theta() * TMath::RadToDeg();
|
||||
thetaB = PB.Theta() * TMath::RadToDeg();
|
||||
|
||||
Tb = Pb.E() - Pb.M();
|
||||
TB = PB.E() - PB.M();
|
||||
|
||||
phib = Pb.Phi() * TMath::RadToDeg();
|
||||
phiB = PB.Phi() * TMath::RadToDeg();
|
||||
|
||||
vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() + vertexXRange[0];
|
||||
vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() + vertexYRange[0];
|
||||
vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0];
|
||||
|
||||
TVector3 vertex(vertexX, vertexY, vertexZ);
|
||||
|
||||
TVector3 dir(1, 0, 0);
|
||||
dir.SetTheta(thetab * TMath::DegToRad());
|
||||
dir.SetPhi(phib * TMath::DegToRad());
|
||||
|
||||
std::pair<int, int> wireID = anasen.FindWireID(vertex, dir, false);
|
||||
SX3 sx3 = anasen.FindSX3Pos(vertex, dir, false);
|
||||
|
||||
anodeID = wireID.first;
|
||||
cathodeID = wireID.second;
|
||||
|
||||
sx3ID = sx3.id;
|
||||
if( sx3.id >= 0 ){
|
||||
sx3Up = sx3.chUp;
|
||||
sx3Down = sx3.chDown;
|
||||
sx3Back = sx3.chBack;
|
||||
sx3ZFrac = sx3.zFrac;
|
||||
|
||||
sx3X = sx3.hitPos.X();
|
||||
sx3Y = sx3.hitPos.Y();
|
||||
sx3Z = sx3.hitPos.Z();
|
||||
|
||||
// for( int i = 0; i < 12; i++){
|
||||
// sx3Index[i] = -1;
|
||||
// if( i == sx3Up ) sx3Index[i] = sx3ID * 12 + sx3Up;
|
||||
// if( i == sx3Down ) sx3Index[i] = sx3ID * 12 + sx3Down;
|
||||
// if( i == sx3Back ) sx3Index[i] = sx3ID * 12 + sx3Back;
|
||||
// }
|
||||
anasen.CalTrack(sx3.hitPos, wireID.first, wireID.second, false);
|
||||
|
||||
reTheta = anasen.GetTrackTheta() * TMath::RadToDeg();
|
||||
rePhi = anasen.GetTrackPhi() * TMath::RadToDeg();
|
||||
|
||||
}else{
|
||||
sx3Up = -1;
|
||||
sx3Down = -1;
|
||||
sx3Back = -1;
|
||||
sx3ZFrac = TMath::QuietNaN();
|
||||
|
||||
sx3X = TMath::QuietNaN();
|
||||
sx3Y = TMath::QuietNaN();
|
||||
sx3Z = TMath::QuietNaN();
|
||||
|
||||
// for( int i = 0; i < 12; i++){
|
||||
// sx3Index[i] = -1;
|
||||
// }
|
||||
|
||||
reTheta = TMath::QuietNaN();
|
||||
rePhi = TMath::QuietNaN();
|
||||
}
|
||||
|
||||
|
||||
|
||||
tree->Fill();
|
||||
|
||||
//#################################################################### Timer
|
||||
clock.Stop("timer");
|
||||
Double_t time = clock.GetRealTime("timer");
|
||||
clock.Start("timer");
|
||||
|
||||
if ( !shown ) {
|
||||
if (fmod(time, 10) < 1 ){
|
||||
printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
|
||||
shown = 1;
|
||||
}
|
||||
}else{
|
||||
if (fmod(time, 10) > 9 ){
|
||||
shown = 0;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
tree->Write();
|
||||
int count = tree->GetEntries();
|
||||
saveFile->Close();
|
||||
|
||||
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
|
||||
|
||||
}
|
||||
106
constant.h
Normal file
106
constant.h
Normal file
|
|
@ -0,0 +1,106 @@
|
|||
/***********************************************************************
|
||||
*
|
||||
* This is constant.h, to provide various physical constants.
|
||||
*
|
||||
*-------------------------------------------------------
|
||||
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
|
||||
* email: goluckyryan@gmail.com
|
||||
* ********************************************************************/
|
||||
|
||||
|
||||
#ifndef constant
|
||||
#define constant
|
||||
#include <cmath>
|
||||
|
||||
const double pi = acos(-1.0);
|
||||
const double E = 2.718281828459 ;
|
||||
const double hbar_SI = 1.054571628e-34; //Js
|
||||
const double kB = 1.3806504e-23; //JK^-1
|
||||
const double e = 1.602176487e-19; //C
|
||||
const double c_SI = 299792458; //ms^-1
|
||||
const double me_SI = 9.10938215e-31 ; //kg
|
||||
const double mp_SI = 1.672621637e-27 ; //kg
|
||||
const double mn_SI = 1.67492729e-27 ; //kg
|
||||
const double NA = 6.022141e+23 ; //mol^-1
|
||||
|
||||
const double deg2rad = pi/180 ;
|
||||
const double rad2deg = 180/pi ;
|
||||
|
||||
//======================================================================
|
||||
const double amu = 931.49432; // MeV/c^2
|
||||
const double hbarc = 197.326979; // MeV fm;
|
||||
const double c = 299.792458; // mm/ns;
|
||||
const double ee = 1.439964454; // MeV.fm
|
||||
|
||||
//======================================================================
|
||||
double kg2MeV(double m){
|
||||
return m*c_SI*c_SI/e/1e6;
|
||||
}
|
||||
|
||||
double T2Brho(double mass, int Z, int A, double T){
|
||||
//mass in MeV
|
||||
// Z in e
|
||||
// T in MeV/A
|
||||
double gamma = (T*A + mass)/mass;
|
||||
double beta = sqrt(1-1/gamma/gamma);
|
||||
return mass*beta*gamma/Z/c;
|
||||
}
|
||||
|
||||
double Brho2T(double mass, int Z, int A, double Brho){
|
||||
//mass in MeV
|
||||
// Z in e
|
||||
return (sqrt(pow(Brho*Z*c,2)+mass*mass)-mass)/A;
|
||||
}
|
||||
|
||||
double T2beta(double mass, int A, double T){
|
||||
double gamma = 1.0 + T*A/mass;
|
||||
return sqrt(1-1/gamma/gamma);
|
||||
}
|
||||
|
||||
double ev2nm(double eV){
|
||||
// photon energy to nm
|
||||
return hbarc/2/pi/eV;
|
||||
}
|
||||
|
||||
//======================================================================
|
||||
const double mp = kg2MeV(mp_SI);
|
||||
const double mn = kg2MeV(mn_SI);
|
||||
const double hbar = 197.326979;
|
||||
|
||||
|
||||
//======================================================================
|
||||
std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
|
||||
|
||||
std::vector<std::string> output;
|
||||
|
||||
size_t pos;
|
||||
do{
|
||||
pos = tempLine.find(splitter); /// fine splitter
|
||||
if( pos == 0 ){ ///check if it is splitter again
|
||||
tempLine = tempLine.substr(pos+1);
|
||||
continue;
|
||||
}
|
||||
|
||||
std::string secStr;
|
||||
if( pos == std::string::npos ){
|
||||
secStr = tempLine;
|
||||
}else{
|
||||
secStr = tempLine.substr(0, pos+shift);
|
||||
tempLine = tempLine.substr(pos+shift);
|
||||
}
|
||||
|
||||
///check if secStr is begin with space
|
||||
while( secStr.substr(0, 1) == " ") secStr = secStr.substr(1);
|
||||
|
||||
///check if secStr is end with space
|
||||
while( secStr.back() == ' ') secStr = secStr.substr(0, secStr.size()-1);
|
||||
|
||||
output.push_back(secStr);
|
||||
///printf(" |%s---\n", secStr.c_str());
|
||||
|
||||
}while(pos != std::string::npos );
|
||||
|
||||
return output;
|
||||
}
|
||||
|
||||
#endif
|
||||
3594
mass20.txt
Normal file
3594
mass20.txt
Normal file
File diff suppressed because it is too large
Load Diff
2
script.C
2
script.C
|
|
@ -253,7 +253,7 @@ void script(TString fileName = "", int maxEvent = -1){
|
|||
int xRan = gRandom->Integer(20) - 10;
|
||||
int yRan = gRandom->Integer(20) - 10;
|
||||
int zRan = gRandom->Integer(20) - 10;
|
||||
int pRan = gRandom->Integer(360); // phi deg
|
||||
int pRan = gRandom->Integer(10) + 175; // phi deg
|
||||
|
||||
int tRan = 0 ;
|
||||
do{
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user