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

init

Browse Source
This commit is contained in:
Ryan Tang 2022-07-13 14:08:02 -04:00
commit 36d8ffeedc
23 changed files with 11578 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

496
Armory/AnalysisLibrary.h Normal file
View File

@ -0,0 +1,496 @@
#ifndef Analysis_Library_h
#define Analysis_Library_h
#include <TF1.h>
#include <TGraph.h>
#include <TSpectrum.h>
#include <TMath.h>
#include <TMacro.h>
#include <TList.h>
#include <iostream>
#include <vector>
#include <string>
struct DetGeo{
double Bfield; /// T
int BfieldSign ; /// sign of B-field
double BfieldTheta; /// rad, 0 = z-axis, pi/2 = y axis, pi = -z axis
double bore; /// bore , mm
double detPerpDist; /// distance from axis
double detWidth; /// width
double detLength; /// length
double recoilPos; /// recoil, downstream
double recoilInnerRadius; /// radius recoil inner
double recoilOuterRadius; /// radius recoil outter
double recoilPos1, recoilPos2; /// imaginary recoils
double elumPos1, elumPos2; /// imaginary elum, only sensitive to light recoil
double blocker;
double firstPos; /// meter
double eSigma; /// intrinsic energy resolution MeV
double zSigma; /// intrinsic position resolution mm
std::vector<double> pos; /// near position in meter
int nDet, mDet; /// nDet = number of different pos, mDet, number of same pos
double zMin, zMax; /// range of detectors
std::vector<double> detPos; ///absolute position of detector
bool isCoincidentWithRecoil;
bool detFaceOut; ///detector_facing_Out_or_In
};
struct ReactionConfig{
int beamA;
int beamZ;
int targetA;
int targetZ;
int recoilLightA;
int recoilLightZ;
int recoilHeavyA;
int 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]
};
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;
}
///Using TMacro to load the detectorGeo frist,
///this indrect method is good for loading detectorGeo from TMacro in root file
DetGeo LoadDetectorGeo(TMacro * macro){
DetGeo detGeo;
if( macro == NULL ) return detGeo;
TList * haha = macro->GetListOfLines();
int numLine = (haha)->GetSize();
detGeo.pos.clear();
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 ) {
detGeo.Bfield = atof(str[0].c_str());
detGeo.BfieldSign = detGeo.Bfield > 0 ? 1: -1;
}
if ( i == 1 ) detGeo.BfieldTheta = atof(str[0].c_str());
if ( i == 2 ) detGeo.bore = atof(str[0].c_str());
if ( i == 3 ) detGeo.detPerpDist = atof(str[0].c_str());
if ( i == 4 ) detGeo.detWidth = atof(str[0].c_str());
if ( i == 5 ) detGeo.detLength = atof(str[0].c_str());
if ( i == 6 ) detGeo.recoilPos = atof(str[0].c_str());
if ( i == 7 ) detGeo.recoilInnerRadius = atof(str[0].c_str());
if ( i == 8 ) detGeo.recoilOuterRadius = atof(str[0].c_str());
if ( i == 9 ) detGeo.isCoincidentWithRecoil = str[0] == "false" ? false: true;
if ( i == 10 ) detGeo.recoilPos1 = atof(str[0].c_str());
if ( i == 11 ) detGeo.recoilPos2 = atof(str[0].c_str());
if ( i == 12 ) detGeo.elumPos1 = atof(str[0].c_str());
if ( i == 13 ) detGeo.elumPos2 = atof(str[0].c_str());
if ( i == 14 ) detGeo.blocker = atof(str[0].c_str());
if ( i == 15 ) detGeo.firstPos = atof(str[0].c_str());
if ( i == 16 ) detGeo.eSigma = atof(str[0].c_str());
if ( i == 17 ) detGeo.zSigma = atof(str[0].c_str());
if ( i == 18 ) detGeo.detFaceOut = str[0] == "Out" ? true : false;
if ( i == 19 ) detGeo.mDet = atoi(str[0].c_str());
if ( i >= 20 ) (detGeo.pos).push_back(atof(str[0].c_str()));
}
detGeo.nDet = (detGeo.pos).size();
(detGeo.detPos).clear();
for(int id = 0; id < detGeo.nDet; id++){
if( detGeo.firstPos > 0 ) detGeo.detPos.push_back(detGeo.firstPos + detGeo.pos[id]);
if( detGeo.firstPos < 0 ) detGeo.detPos.push_back(detGeo.firstPos - detGeo.pos[detGeo.nDet - 1 - id]);
///printf("%d | %f, %f \n", id, detGeo.pos[id], detGeo.detPos[id]);
}
detGeo.zMin = TMath::Min(detGeo.detPos.front(), detGeo.detPos.back()) - (detGeo.firstPos < 0 ? detGeo.detLength : 0);
detGeo.zMax = TMath::Max(detGeo.detPos.front(), detGeo.detPos.back()) + (detGeo.firstPos > 0 ? detGeo.detLength : 0);
return detGeo;
}
void PrintDetGeo(DetGeo detGeo){
printf("=====================================================\n");
printf(" B-field: %8.2f T, Theta : %6.2f deg \n", detGeo.Bfield, detGeo.BfieldTheta);
if( detGeo.BfieldTheta != 0.0 ) {
printf(" +---- field angle != 0 is not supported!!! \n");
}
printf(" Recoil detector pos: %8.2f mm, radius: %6.2f - %6.2f mm \n", detGeo.recoilPos, detGeo.recoilInnerRadius, detGeo.recoilOuterRadius);
printf(" Blocker Position: %8.2f mm \n", detGeo.firstPos > 0 ? detGeo.firstPos - detGeo.blocker : detGeo.firstPos + detGeo.blocker );
printf(" First Position: %8.2f mm \n", detGeo.firstPos);
printf("------------------------------------- Detector Position \n");
for(int i = 0; i < detGeo.nDet ; i++){
if( detGeo.firstPos > 0 ){
printf("%d, %8.2f mm - %8.2f mm \n", i, detGeo.detPos[i], detGeo.detPos[i] + detGeo.detLength);
}else{
printf("%d, %8.2f mm - %8.2f mm \n", i, detGeo.detPos[i] - detGeo.detLength , detGeo.detPos[i]);
}
}
printf(" number of det : %d x %d \n", detGeo.mDet, detGeo.nDet);
printf(" detector facing : %s\n", detGeo.detFaceOut ? "Out" : "In");
printf(" energy resol.: %f MeV\n", detGeo.eSigma);
printf(" pos-Z resol.: %f mm \n", detGeo.zSigma);
if( detGeo.elumPos1 != 0 || detGeo.elumPos2 != 0 || detGeo.recoilPos1 != 0 || detGeo.recoilPos2 != 0){
printf("=================================== Auxillary/Imaginary Detectors\n");
}
if( detGeo.elumPos1 != 0 ) printf(" Elum 1 pos.: %f mm \n", detGeo.elumPos1);
if( detGeo.elumPos2 != 0 ) printf(" Elum 2 pos.: %f mm \n", detGeo.elumPos2);
if( detGeo.recoilPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", detGeo.recoilPos1);
if( detGeo.recoilPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", detGeo.recoilPos2);
printf("=====================================================\n");
}
ReactionConfig LoadReactionConfig(TMacro * macro){
ReactionConfig reaction;
if( macro == NULL ) return reaction;
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 ) reaction.beamA = atoi(str[0].c_str());
if( i == 1 ) reaction.beamZ = atoi(str[0].c_str());
if( i == 2 ) reaction.targetA = atoi(str[0].c_str());
if( i == 3 ) reaction.targetZ = atoi(str[0].c_str());
if( i == 4 ) reaction.recoilLightA = atoi(str[0].c_str());
if( i == 5 ) reaction.recoilLightZ = atoi(str[0].c_str());
if( i == 6 ) reaction.beamEnergy = atof(str[0].c_str());
if( i == 7 ) reaction.beamEnergySigma = atof(str[0].c_str());
if( i == 8 ) reaction.beamAngle = atof(str[0].c_str());
if( i == 9 ) reaction.beamAngleSigma = atof(str[0].c_str());
if( i == 10 ) reaction.beamX = atof(str[0].c_str());
if( i == 11 ) reaction.beamY = atof(str[0].c_str());
if( i == 12 ) reaction.numEvents = atoi(str[0].c_str());
if( i == 13 ) {
if( str[0].compare("false") == 0 ) reaction.isTargetScattering = false;
if( str[0].compare("true") == 0 ) reaction.isTargetScattering = true;
}
if( i == 14 ) reaction.targetDensity = atof(str[0].c_str());
if( i == 15 ) reaction.targetThickness = atof(str[0].c_str());
if( i == 16 ) reaction.beamStoppingPowerFile = str[0];
if( i == 17 ) reaction.recoilLightStoppingPowerFile = str[0];
if( i == 18 ) reaction.recoilHeavyStoppingPowerFile = str[0];
if( i == 19 ) {
if( str[0].compare("false") == 0 ) reaction.isDecay = false;
if( str[0].compare("true") == 0 ) reaction.isDecay = true;
}
if( i == 20 ) reaction.heavyDecayA = atoi(str[0].c_str());
if( i == 21 ) reaction.heavyDecayZ = atoi(str[0].c_str());
if( i == 22 ) {
if( str[0].compare("false") == 0 ) reaction.isRedo = false;
if( str[0].compare("true" ) == 0 ) reaction.isRedo = true;
}
if( i >= 23) {
reaction.beamEx.push_back( atof(str[0].c_str()) );
}
}
reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA;
reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ;
return reaction;
}
void PrintReactionConfig(ReactionConfig reaction){
printf("=====================================================\n");
printf(" beam : A = %3d, Z = %2d \n", reaction.beamA, reaction.beamZ);
printf(" target : A = %3d, Z = %2d \n", reaction.targetA, reaction.targetZ);
printf(" light : A = %3d, Z = %2d \n", reaction.recoilLightA, reaction.recoilLightZ);
printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", reaction.beamEnergy, reaction.beamEnergySigma, reaction.beamEnergySigma/reaction.beamEnergy);
printf(" Angle : %.2f +- %.2f mrad\n", reaction.beamAngle, reaction.beamAngleSigma);
printf(" offset : (x,y) = (%.2f, %.2f) mmm \n", reaction.beamX, reaction.beamY);
printf("##### number of Simulation Events : %d \n", reaction.numEvents);
printf(" is target scattering : %s \n", reaction.isTargetScattering ? "Yes" : "No");
if(reaction.isTargetScattering){
printf(" target density : %.f g/cm3\n", reaction.targetDensity);
printf(" thickness : %.f cm\n", reaction.targetThickness);
printf(" beam stopping file : %s \n", reaction.beamStoppingPowerFile.c_str());
printf(" recoil light stopping file : %s \n", reaction.recoilLightStoppingPowerFile.c_str());
printf(" recoil heavy stopping file : %s \n", reaction.recoilHeavyStoppingPowerFile.c_str());
}
printf(" is simulate decay : %s \n", reaction.isDecay ? "Yes" : "No");
if( reaction.isDecay ){
printf(" heavy decay : A = %d, Z = %d \n", reaction.heavyDecayA, reaction.heavyDecayZ);
}
printf(" is Redo until hit array : %s \n", reaction.isRedo ? "Yes" : "No");
printf(" beam Ex : %.2f MeV \n", reaction.beamEx[0]);
for( int i = 1; i < (int) reaction.beamEx.size(); i++){
printf(" %.2f MeV \n", reaction.beamEx[i]);
}
printf("=====================================================\n");
}
std::vector<std::vector<double>> combination(std::vector<double> arr, int r){
std::vector<std::vector<double>> output;
int n = arr.size();
std::vector<int> v(n);
std::fill(v.begin(), v.begin()+r, 1);
do {
//for( int i = 0; i < n; i++) { printf("%d ", v[i]); }; printf("\n");
std::vector<double> temp;
for (int i = 0; i < n; ++i) {
if (v[i]) {
//printf("%.1f, ", arr[i]);
temp.push_back(arr[i]);
}
}
//printf("\n");
output.push_back(temp);
} while (std::prev_permutation(v.begin(), v.end()));
return output;
}
double* sumMeanVar(std::vector<double> data){
int n = data.size();
double sum = 0;
for( int k = 0; k < n; k++) sum += data[k];
double mean = sum/n;
double var = 0;
for( int k = 0; k < n; k++) var += pow(data[k] - mean,2);
static double output[3];
output[0] = sum;
output[1] = mean;
output[2] = var;
return output;
}
double* fitSlopeIntercept(std::vector<double> dataX, std::vector<double> dataY){
double * smvY = sumMeanVar(dataY);
double sumY = smvY[0];
double meanY = smvY[1];
double * smvX = sumMeanVar(dataX);
double sumX = smvX[0];
double meanX = smvX[1];
double varX = smvX[2];
int n = dataX.size();
double sumXY = 0;
for( int j = 0; j < n; j++) sumXY += dataX[j] * dataY[j];
double slope = ( sumXY - sumX * sumY/n ) / varX;
double intercept = meanY - slope * meanX;
static double output[2];
output[0] = slope;
output[1] = intercept;
return output;
}
std::vector<std::vector<double>> FindMatchingPair(std::vector<double> enX, std::vector<double> enY){
//output[0] = fitEnergy;
//output[1] = refEnergy;
int nX = enX.size();
int nY = enY.size();
std::vector<double> fitEnergy;
std::vector<double> refEnergy;
if( nX > nY ){
std::vector<std::vector<double>> output = combination(enX, nY);
double * smvY = sumMeanVar(enY);
double sumY = smvY[0];
double meanY = smvY[1];
double varY = smvY[2];
double optRSquared = 0;
double absRSqMinusOne = 1;
int maxID = 0;
for( int k = 0; k < (int) output.size(); k++){
double * smvX = sumMeanVar(output[k]);
double sumX = smvX[0];
double meanX = smvX[1];
double varX = smvX[2];
double sumXY = 0;
for( int j = 0; j < nY; j++) sumXY += output[k][j] * enY[j];
double rSq = abs(sumXY - sumX*sumY/nY)/sqrt(varX*varY);
//for( int j = 0; j < nY ; j++){ printf("%.1f, ", output[k][j]); }; printf("| %.10f\n", rSq);
if( abs(rSq-1) < absRSqMinusOne ) {
absRSqMinusOne = abs(rSq-1);
optRSquared = rSq;
maxID = k;
}
}
fitEnergy = output[maxID];
refEnergy = enY;
printf(" R^2 : %.20f\n", optRSquared);
//calculation fitting coefficient
//double * si = fitSlopeIntercept(fitEnergy, refEnergy);
//printf( " y = %.4f x + %.4f\n", si[0], si[1]);
}else if( nX < nY ){
std::vector<std::vector<double>> output = combination(enY, nX);
double * smvX = sumMeanVar(enX);
double sumX = smvX[0];
double meanX = smvX[1];
double varX = smvX[2];
double optRSquared = 0;
double absRSqMinusOne = 1;
int maxID = 0;
for( int k = 0; k < (int) output.size(); k++){
double * smvY = sumMeanVar(output[k]);
double sumY = smvY[0];
double meanY = smvY[1];
double varY = smvY[2];
double sumXY = 0;
for( int j = 0; j < nX; j++) sumXY += output[k][j] * enX[j];
double rSq = abs(sumXY - sumX*sumY/nX)/sqrt(varX*varY);
//for( int j = 0; j < nX ; j++){ printf("%.1f, ", output[k][j]); }; printf("| %.10f\n", rSq);
if( abs(rSq-1) < absRSqMinusOne ) {
absRSqMinusOne = abs(rSq-1);
optRSquared = rSq;
maxID = k;
}
}
fitEnergy = enX;
refEnergy = output[maxID];
printf(" R^2 : %.20f\n", optRSquared);
}else{
fitEnergy = enX;
refEnergy = enY;
//if nX == nY, ther could be cases that only partial enX and enY are matched.
}
printf("fitEnergy = ");for( int k = 0; k < (int) fitEnergy.size() ; k++){ printf("%7.2f, ", fitEnergy[k]); }; printf("\n");
printf("refEnergy = ");for( int k = 0; k < (int) refEnergy.size() ; k++){ printf("%7.2f, ", refEnergy[k]); }; printf("\n");
std::vector<std::vector<double>> haha;
haha.push_back(fitEnergy);
haha.push_back(refEnergy);
return haha;
}
#endif

102
Cleopatra/Cleopatra.C Normal file
View File

@ -0,0 +1,102 @@
/***********************************************************************
*
* This is Cleopatra, a sucessor of Ptolemy
* only for (d,p), (d,p), (d,d), or (p,p)
*
* 1) it read a simple infile.in from reaction_setting file
* 2) use Ptolemy to calculate the the creation
* 3) extract the cross sec distribution into txt and root file
*
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ cleopatra.C -o cleopatra `root-config --cflags --glibs`
*
*------------------------------------------------------
* The reaction_setting file is simple like:
*
* 206Hg(d,p)207Hg(1s1/2 0.000) 10MeV/u AK
*
* the first is similar to usual reaction setting, the word AK is a
* short name for Optical Potential, user can put as many line as
* they like, Cleopatra can create the suitable infile.in for Ptolomy
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h> /* atof */
#include <cmath>
#include <vector>
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <string>
#include <TROOT.h>
#include <TFile.h>
#include <TString.h>
#include "InFileCreator.h"
#include "ExtractXSec.h"
#include <TApplication.h>
#include "PlotTGraphTObjArray.h"
using namespace std;
int main (int argc, char *argv[]) { //TODO add angle range
printf("=================================================================\n");
printf("===== Cleopatra, Ptolemy for (d,p),(p,d), (p,p) and (d,d) =====\n");
printf("=================================================================\n");
if(argc < 2 || argc > 5) {
printf("Usage: ./Cleopatra input_file (angMin = 0 deg, angMax = 180 deg, angStep = 1 deg)\n");
printf("Usage: ./Cleopatra input_file angMin angMax (angStep = 1 deg)\n");
printf("Usage: ./Cleopatra input_file angMin angMax angStep\n");
exit(0);
}else{
printf("From file : %s \n", argv[1]);
}
//================= read infile. extract the reactions, write pptolemy infile for each reaction
string readFile = argv[1];
double angMin = 0.;
double angMax = 180.;
double angStep = 1.;
if( argc >= 4 ){
angMin = atof(argv[2]);
angMax = atof(argv[3]);
}
if( argc == 5 ){
angStep = atof(argv[4]);
}
string ptolemyInFileName = argv[1];
ptolemyInFileName += ".in";
printf("=================== Create InFile\n");
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());
printf("=================== Run Ptolemy\n");
printf("%s \n", command);
system(command);
//================= extract the Xsec and save as txt and root
printf("=================== Extract Cross-section\n");
ExtractXSec(ptolemyOutFileName);
//================= Call root to plot the d.s.c.
printf("=================== Plot Result using ROOT.\n");
string rootFileName = argv[1];
rootFileName += ".root";
TApplication app ("app", &argc, argv);
PlotTGraphTObjArray(rootFileName);
app.Run(); //nothing run after
return 0;
}

205
Cleopatra/Cleopatra.sh Executable file
View File

@ -0,0 +1,205 @@
#!/bin/bash
########################################################################
#
# This is Cleopatra.sh, a scripted version for Cleopatra
#
# Using bash script provide flexibility that user can add difference
# compoenents during the calculation
#
# A full package includes fellowing:
# 1) create a in-file for ptolemy
# 2) run ptolemy from that in-file and output an out-file
# 3) extract cross-section distribution from the out-file
# save as txt or root TGraph format
# 4) call ROOT to draw the TGraph
# 5) load possible experimental Xsec and fit with Ptolemy calulation
#
# User can easily select/comment-out different component
# to suit their needs
#-------------------------------------------------------
# created by Ryan (Tsz Leung) Tang, Nov-18, 2018
# email: goluckyryan@gmail.com
########################################################################
#===== Call thisroot.h
ROOTPATH=$(which root)
len=${#ROOTPATH}
ROOTSOURCE=${ROOTPATH:0:$len-4}"thisroot.sh"
echo $ROOTSOURCE
source $ROOTSOURCE
#===== go to Cleopatra and make
cd ../Cleopatra
make
cd ../working
#================================ User Defualt Control
CreateInFile=0 # 0 = false, 1 = true
RunPtolemy=0
IsExtractXSec=0
PlotResult=0
SimTransfer=0
#============================================ USER don't need to change thing below
if [ $# -eq 0 ] ; then
echo "$./Cleopatra in-file X X X X X angMin angMax angStep"
echo " | | | | |"
echo " | | | | Simulate Transfer reaction? (1/0)"
echo " | | | |"
echo " | | | PlotResult? (1/0)"
echo " | | Extract cross-section? (2/1/0)"
echo " | | if 1 = extract Ratio to Rutherford for (d,d), (p,p)"
echo " | | if 2 = extract total Xsec for (d,d), (p,p), (n,n)"
echo " | | if 3 = extract Rutherford"
echo " | Run Ptolemy? (1/0)"
echo " Create infile? (1/0)"
echo "default angMin = 0, angMax = 50, angStep = 0.5"
exit 1
fi;
loadfile=$1
infile=$1".in"
outfile=$1".out"
rootfile=$1".root"
exFile=$1".Ex.txt"
if [ $# -eq 2 ]; then
CreateInFile=$2
fi;
if [ $# -eq 3 ]; then
CreateInFile=$2
RunPtolemy=$3
fi;
if [ $# -eq 4 ]; then
CreateInFile=$2
RunPtolemy=$3
IsExtractXSec=$4
fi;
if [ $# -eq 5 ]; then
CreateInFile=$2
RunPtolemy=$3
IsExtractXSec=$4
PlotResult=$5
fi;
if [ $# -eq 6 ]; then
CreateInFile=$2
RunPtolemy=$3
IsExtractXSec=$4
PlotResult=$5
SimTransfer=$6
fi;
if [ $# -eq 9 ]; then
CreateInFile=$2
RunPtolemy=$3
IsExtractXSec=$4
PlotResult=$5
SimTransfer=$6
angMin=$7
angMax=$8
angStep=$9
fi;
ExtractXsecMsg=""
if [ $IsExtractXSec -eq 1 ]; then
ExtractXsecMsg=", for (d,d)(p,p), extract Ratio to Rutherford"
elif [ $IsExtractXSec -eq 2 ]; then
ExtractXsecMsg=", for (d,d)(p,p), extract Total Xsec"
fi;
if [ $SimTransfer -eq 1 ]; then
angMin=0
angMax=180
angStep=0.5
fi
echo "#################################################################"
echo "## @@@@ @@ @@@@ @@@@ @@@@@ @@@@ @@@@@@ @@@@@ @@@@ ##"
echo "## @@ @@ @@ @@ @@ @@ @@ @@ @@ @@ @@ @@ @@ @@ ##"
echo "## @@ @@ @@@@ @@ @@ @@@@@ @@@@@@ @@ @@@@@ @@@@@@ ##"
echo "## @@ @@ @@ @@ @@ @@ @@ @@ @@ @@ @ @@ @@ ##"
echo "## @@@@ @@@@@ @@@@ @@@@ @@ @@ @@ @@ @@ @ @@ @@ ##"
echo "#################################################################"
echo "##### Cleopatra, Ptolemy for p,d,t,3He, a #####"
echo "#################################################################"
echo "Supported reactions:"
echo "elastics/inelastics : (p,p), (d,d)"
echo "1-neutron adding : (d,p), (t,d) | (a,3He)?"
echo "1-proton adding : (3He,d)"
echo "2-neutrons adding : (t,p)"
echo "1-neutron removal : (p,d), (d,t) "
echo "1-proton removal : (d,3He), (t,a)"
echo "================================================================="
echo "USER OPTION:"
echo " --- Is Create Ptolemy infile ? " ${CreateInFile}
echo " --- Is Run Ptolemy ? " ${RunPtolemy}
echo " --- Is Extract Cross-Section ? " ${IsExtractXSec} ${ExtractXsecMsg}
echo " --- Is Plot Results ? " ${PlotResult}
echo " ----Is Simulation Transfer ? " ${SimTransfer}
echo "================================================================="
#if [ ${CreateInFile} -eq 1 ] ; then
# echo "infile ----> "${loadfile}
#fi;
#
#if [ ${RunPtolemy} -eq 1 ] ; then
# echo "Ptolemy infile ----> "${infile}
# echo "Ptolemy outfile ----> "${outfile}
#fi;
if [ ${CreateInFile} -eq 1 ] ; then
if [ $# -eq 9 ]; then
../Cleopatra/InFileCreator ${loadfile} $angMin $angMax $angStep
else
../Cleopatra/InFileCreator ${loadfile} 0.0 50.0 0.5
fi
fi;
if [ ${RunPtolemy} -eq 1 ] ; then
echo "================================================================="
echo "===== Ptolemy Calcualtion ==================================="
echo "================================================================="
#check is linux or Mac
arch=$(uname)
if [ ${arch} == "Darwin" ] ; then
../Cleopatra/ptolemy_mac <${infile}> ${outfile}
ptolemyOUTPUT=$?
else
../Cleopatra/ptolemy <${infile}> ${outfile}
ptolemyOUTPUT=$?
fi
echo "ptolmey exit code : " $ptolemyOUTPUT
if [ ${ptolemyOUTPUT} -eq 0 ] ; then
echo "Ptolmey finished without problem. "
else
echo "Ptolemy has error, check ${infile} or ${outfile}"
exit 1;
fi
fi;
#===== Extracting XSec and save into *txt and *root
if [ ${IsExtractXSec} -ge 1 ] ; then
../Cleopatra/ExtractXSec ${outfile} ${IsExtractXSec}
fi;
if [ ${PlotResult} -eq 1 ] ; then
#===== Plot the result from the *.root
#./PlotTGraphTObjArray ${rootfile}
#--- other way within ROOT
echo "================================================================="
echo "===== Plot Result from ${rootfile}"
echo "================================================================="
com='../Cleopatra/PlotTGraphTObjArray.h("'${rootfile}'")'
echo ${com}
root -l ${com}
fi;
if [ ${SimTransfer} -eq 1 ] ; then
../Cleopatra/Transfer reactionConfig.txt detectorGeo.txt ${exFile} ${rootfile} transfer.root reaction.dat
fi;

73
Cleopatra/DWBARatio.C Normal file
View File

@ -0,0 +1,73 @@
TGraph * DWBARatio(int id1, int id2, TString rootFile="DWBA.root", bool isPlot = true){
TGraph * gR = NULL;
TFile * file = new TFile(rootFile, "READ");
if( file != NULL ){
printf("----- Opening %s\n", rootFile.Data());
}else{
printf("----- Opening %s .... Fail.\n", rootFile.Data());
return gR;
}
///get the TGraph of the distribution.
TObjArray * gList = (TObjArray *) file->Get("qList");
int size = gList->GetLast()+1;
printf("----- found %d d.s.c\n", size);
if( id1 > size || id2 > size ) {
printf(" id1 > %d || id2 > %d \n", size, size);
return gR;
}
TGraph * g1 = (TGraph *) gList->At(id1);
TGraph * g2 = (TGraph *) gList->At(id2);
double g1MaxY = g1->GetHistogram()->GetMaximum();
double g2MaxY = g2->GetHistogram()->GetMaximum();
g2->SetLineColor(2);
TCanvas * cDWBA = NULL ;
if( isPlot ){
cDWBA= new TCanvas("cDWBA", "DWBA Ratio", 1000, 500);
cDWBA->Divide(2,1);
cDWBA->cd(1);
cDWBA->cd(1)->SetLogy();
if( g1MaxY > g2MaxY ) {
g1->Draw();
g2->Draw("same");
}else{
g2->Draw();
g1->Draw("same");
}
TLegend * legend = new TLegend( 0.1, 0.9, 0.9, 1.0);
legend->AddEntry(g1, g1->GetName());
legend->AddEntry(g2, g2->GetName());
legend->Draw();
cDWBA->cd(2);
}
gR = new TGraph();
double x, y1, y2;
for( int i = 0 ; i < g1->GetN(); i++){
g1->GetPoint(i, x, y1);
g2->GetPoint(i, x, y2);
gR->SetPoint(i, x, y1/y2);
}
if( isPlot) gR->Draw();
return gR;
}

16
Cleopatra/DWBA_compare.C Normal file
View File

@ -0,0 +1,16 @@
#include "DWBARatio.C"
void DWBA_compare(){
TGraph * w[12];
for( int i = 0 ; i < 12; i++){
w[i] = DWBARatio(i, i+12, "DWBA.root", false);
w[i]->SetLineColor(i+1);
i == 0 ? w[i]->Draw("Al") : w[i]->Draw("same");
}
}

48
Cleopatra/ExtractXSec.C Normal file
View File

@ -0,0 +1,48 @@
/***********************************************************************
*
* This is ExtractXSec for *.out for Ptolemy
*
* This program will extract cross section distribution from Ptolmey output.
* save as *.Xsec.txt and *.root for distribution
*
* save ExPtolemy.txt for excitation energies and total X-section
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ ExtractXSec.C -o ExtractXSec `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h>
#include "ExtractXSec.h"
using namespace std;
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("========== Extract Cross-Section From Ptolemy out file ==========\n");
printf("=================================================================\n");
if(argc < 2 || argc > 3) {
printf("Usage: ./ExtractXSec input_file <ElasticFlag>\n");
printf(" ElasticFlag = 1 , default, extarct Ratio to Rutherford\n");
printf(" ElasticFlag = 2 , extarct Total Xsec\n");
printf(" ElasticFlag = 3 , Rutherford\n");
exit(0);
}else{
printf("From file : %s \n", argv[1]);
}
string readFile = argv[1];
int ElasticFlag = 1;
if( argc == 3 ){
ElasticFlag = atoi(argv[2]);
}
ExtractXSec(readFile, ElasticFlag);
}

554
Cleopatra/ExtractXSec.h Normal file
View File

@ -0,0 +1,554 @@
/***********************************************************************
*
* This is ExtractXSec for *.out for Ptolemy
*
* This program will extract cross section distribution from Ptolmey output.
* save as *.Xsec.txt and *.root for distribution
*
* save ExPtolemy.txt for excitation energies and total X-section
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ ExtractXSec.C -o ExtractXSec `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h>
#include <cmath>
#include <vector>
#include <TROOT.h>
#include <TFile.h>
#include <TString.h>
#include <TMath.h>
#include <TGraph.h>
#include <TF1.h>
#include <TObjArray.h>
#include "../Armory/AnalysisLibrary.h"
using namespace std;
TObjArray * gList = NULL;
double distfunct(double *x, double *par){
TGraph * gTemp = (TGraph *) gList->At(par[0]);
return (gTemp->Eval(x[0]))* sin(x[0] * TMath::DegToRad());
}
bool isFloat( string str ) {
int length = str.length();
for( int i = 0; i < length; i++){
string it = str.substr(i,1);
if( it == " " || it == "."|| it == "1"|| it == "2"||
it == "3"|| it == "4"|| it == "5"|| it == "6"|| it == "7"|| it == "8"|| it == "9"|| it == "0"){
continue;
}else{
return false;
}
}
return true;
}
int ExtractXSec (string readFile, int indexForElastic=1) {
//indexForElastic = 1 ; for Ratio
//indexForElastic = 2 ; for Total
//indexForElastic = 3 ; for Rutherford
//--- open file.out
ifstream file_in;
file_in.open(readFile.c_str(), ios::in);
if( !file_in){
printf("Unable to open %s \n", readFile.c_str());
exit(1);
}
//---- variable for Xsec
vector<string> title;
vector<vector<double>> dataMatrix;
vector<double> angle;
vector<double> Ex;
vector<double> totalXsec;
vector<string> reaction;
double angleMin = 0;
double angleMax = 0;
double angleStep = -1;
//================================== read file.out
string line;
int lineNum = 0;
vector<double> dataXsec;
bool startExtract = false;
bool angleFilled = false;
int numCal = 0;
int reactionFlag = 0;
bool preFindForElastic = false;
printf("======================================================\n");
while(getline(file_in, line)){
lineNum ++;
//printf("%d , %s\n", lineNum, line.c_str());
size_t pos;
string findStr;
int findLen;
//---- find Title
findStr = "$============================================";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos){
title.push_back(line.substr(pos + findLen+1));
pos = title.back().find("(");
string ex = title.back().substr(3, pos-3);
Ex.push_back( atof(ex.c_str()) );
//clear dataXsex for new data set
dataXsec.clear();
numCal ++;
continue;
}
//---- find Reaction
findStr = "0INPUT... CHANNEL";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos ) {
reaction.push_back( line.substr(pos + findLen + 1) );
reactionFlag = 1; // 1 for (d,d) or (p,p)
///printf("%d |----------- (d,d), %d\n", lineNum, reactionFlag);
continue; // nextline
}
findStr = "REACTION:";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos ) {
reaction.push_back( line.substr(pos + findLen + 1) );
reactionFlag = 2; // 2 for (d,p) or (p,d)
///printf("%d |----------- (d,p), %d\n", lineNum, reactionFlag);
continue; // nextline
}
//----- pre find for Elastic scattering
findStr = "0 OPTICAL MODEL SCATTERING FOR THE OUTGOING CHANNEL";
pos = line.find(findStr);
if( pos != string::npos ) {
preFindForElastic = true;
///printf("%d |----------- pre Elastics Flag : %d\n", lineNum, preFindForElastic);
continue;
}
//----- find angle stetting when not known
if( angleStep == -1 ){
findStr = "anglemin=";
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos){
size_t pos1 = line.find("anglemax=");
angleMin = atof( line.substr(pos+findLen, pos1 - pos - findLen).c_str() );
pos = pos1;
pos1 = line.find("anglestep=");
angleMax = atof( line.substr(pos+findLen, pos1 - pos - findLen).c_str() );
angleStep = atof( line.substr(pos1 + findLen+1).c_str() );
//printf("%d |---- angle range found.\n", lineNum);
}
continue; //nextline
}
//----- check if start extracting Xsec or not
findStr = "dumpdumpdump";
if( reactionFlag == 1 && preFindForElastic ){
findStr = "C.M. LAB RUTHERFORD";
}
if( reactionFlag == 2 ){
findStr = "0 C.M. REACTION REACTION LOW L HIGH L % FROM";
}
pos = line.find(findStr);
if( pos != string::npos ) {
startExtract = true;
//printf("%d |----- start extraction \n", lineNum);
continue;
}
//----- start extracting Xsec
if( startExtract ){
if( line.length() < 20 ) continue;
///printf("%d |%s \n", line.length(), line.c_str());
double numAngle, numXsec;
if( reactionFlag == 1 ){ // Elastics (d,d) or (p,p)
numAngle = atof( line.substr(0, 7).c_str());
if( indexForElastic == 1 ){
numXsec = atof( line.substr(15, 15).c_str());
}else if( indexForElastic == 2 ){
if ( line.length() > 60 ) {
numXsec = atof( line.substr(30, 13).c_str());
}else{
numXsec = -404;
}
}else{
if ( line.length() > 60 ) {
numXsec = atof( line.substr(57, 14).c_str());
}else{
numXsec = -404;
}
}
}
if( reactionFlag == 2 ){ // InElastics (d,p) or (p,d)
if( isFloat( line.substr(0, 7) ) ){
numAngle = atof( line.substr(0, 7).c_str());
numXsec = atof( line.substr(7, 19).c_str());
}else{
numAngle = -1.0;
numXsec = -1.0;
}
}
if( numAngle != numXsec && numXsec > 0. ){
if( !angleFilled ) angle.push_back(numAngle);
dataXsec.push_back(numXsec);
///printf(" %f , %f \n", angle.back(), dataXsec.back());
}
}
//------ find total Xsec, if found stop extraction
findStr = "dumpdumpdump";
if( reactionFlag == 1 && preFindForElastic){
findStr = "0TOTAL REACTION CROSS SECTION =";
}
if( reactionFlag == 2){
findStr = "0TOTAL:";
}
findLen = findStr.length();
pos = line.find(findStr);
if( pos != string::npos ) {
totalXsec.push_back( atof(line.substr(pos + findLen).c_str()) );
printf("%2d | %s | total Xsec(4pi): %f mb \n", numCal, title.back().c_str(), totalXsec.back());
//push back dataXsec to dataMatrix
dataMatrix.push_back(dataXsec);
//printf("%d |----- end extraction \n", lineNum);
angleFilled = true;
startExtract = false;
reactionFlag = 0;
preFindForElastic = false;
continue;
}
}
file_in.close();
//================================== summary
printf("====================== Total number of line read : %d \n", lineNum);
printf("Angle : %5.2f, %5.2f | step : %5.2f \n", angleMin, angleMax, angleStep);
printf("Number of Angle read : %lu \n", angle.size());
printf("Number of data read : %lu \n", dataXsec.size());
printf("Number of Reaction : %lu \n", reaction.size());
printf("----------------------------- list of Calculation \n");
//... find suitable lenght for displaying reaction string
int reactionStrLen = 0;
for( int i = 0; i < numCal ; i++){
int len = reaction[i].length();
if( len > reactionStrLen ) reactionStrLen = len;
}
vector<double> partialXsec;
partialXsec.clear();
for( int i = 0; i < numCal ; i++){
double partialSumXsec = 0.0;
for( int j = 0; j < (dataMatrix[i]).size() ; j++ ){
//double theta = (angle[j] + angleStep/2.) * TMath::DegToRad();
double theta = (angle[j]) * TMath::DegToRad();
double dTheta = angleStep * TMath::DegToRad();
double phi = TMath::TwoPi();
partialSumXsec += dataMatrix[i][j] * sin( theta ) * dTheta * phi ;
}
partialXsec.push_back(partialSumXsec);
size_t pos = title[i].find(")");
printf("%*s| %s | Xsec(%3.0f-%3.0f deg) : %f mb\n", reactionStrLen + 3, reaction[i].c_str(), title[i].substr(pos+1).c_str(), angleMin, angleMax, partialSumXsec);
}
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();
saveFileName = saveFileName.substr(0, len - 4);
saveFileName += ".Xsec.txt";
printf("Output : %s \n", saveFileName.c_str());
FILE * file_out;
file_out = fopen(saveFileName.c_str(), "w+");
for( int i = 0; i < numCal ; i++){
fprintf(file_out, "#%14s\n", reaction[i].c_str());
}
int space = 19;
fprintf(file_out, "%8s\t", "Angel");
for( int i = 0; i < numCal ; i++){
fprintf(file_out, "%*s", space, title[i].c_str());
}
fprintf(file_out, "\n");
for( int i = 0; i < angle.size() ; i ++){
fprintf(file_out, "%8.3f\t", angle[i]);
for( int j = 0; j < numCal ; j++){
fprintf(file_out, "%*f", space, dataMatrix[j][i]);
}
fprintf(file_out, "\n");
}
fclose(file_out);
//================================== Save in ROOT
len = saveFileName.length();
saveFileName = saveFileName.substr(0, len - 9);
TString fileName = saveFileName;
fileName += ".root";
printf("Output : %s \n", fileName.Data());
TFile * fileOut = new TFile(fileName, "RECREATE" );
gList = new TObjArray(); ///no SetTitle() method for TObjArray
gList->SetName("TGraph of d.s.c");
TObjArray * fList = new TObjArray();
fList->SetName("TF1 of distributions = d.s.c. * sin()");
TGraph ** gGraph = new TGraph *[numCal];
TF1 ** dist = new TF1*[numCal];
for( int i = 0; i < numCal ; i++){
gGraph[i] = new TGraph();
TString name = reaction[i];
name += "|";
name += title[i];
gGraph[i]->SetName(name);
for( int j = 0; j < angle.size() ; j++){
gGraph[i]->SetPoint(j, angle[j], dataMatrix[i][j]);
}
gList->Add(gGraph[i]);
name.Form("dist%d", i);
dist[i] = new TF1(name, distfunct, angleMin, angleMax, 1 );
dist[i]->SetParameter(0, i);
fList->Add(dist[i]);
//delete tempFunc;
}
gList->Write("qList", 1);
fList->Write("pList", 1);
fileOut->Write();
fileOut->Close();
printf("---------------------------------------------------\n");
return 0;
}
int ExtractXSecFromText(string readFile){
//read file
ifstream file_in;
file_in.open(readFile.c_str(), ios::in);
if( !file_in){
printf("Unable to open %s \n", readFile.c_str());
exit(1);
}
//---- variable for Xsec
vector<vector<double>> dataMatrix;
vector<double> angle;
vector<double> Ex;
//================================== read file.out
string line;
int lineNum = 0;
vector<double> dataXsec;
int numCal = 0;
printf("======================================================\n");
bool headerDone = false;
while(getline(file_in, line)){
lineNum ++;
if( line.substr(0,1) == "#" ) continue;
//printf("%d | %s\n", lineNum, line.c_str());
//after the comment line, the next line must be column name
vector<string> header= SplitStr(line, " ");
//printf("---%lu #", header.size());
//for( int i = 0; i < header.size(); i++){
// printf("%s|", header[i].c_str());
//}
//printf("\n");
if ( !headerDone ) {
for( int i = 1 ; i < header.size(); i++){
string energy = header[i].substr(3, header[i].length());
Ex.push_back(atof(energy.c_str()));
//printf("%f---\n", Ex.back());
}
headerDone = true;
}else{
vector<double> temp;
for( int i = 0; i < header.size(); i++){
if( i == 0 ) angle.push_back(atof(header[i].c_str()));
if( i > 0 ) temp.push_back(atof(header[i].c_str()));
}
dataMatrix.push_back(temp);
}
}
file_in.close();
double angleMin = angle.front();
double angleMax = angle.back();
double angleStep = (angleMax - angleMin)/(angle.size()-1);
//------ print out read data
printf("====================== data read as:\n");
printf("%5s, ", "Ex");
for(int i = 0; i < Ex.size(); i++){
printf("%6.3f|", Ex[i]);
}
printf("\n");
for(int i = 0; i < dataMatrix.size(); i++){
printf("%5.2f, ", angle[i]);
for( int j = 0; j < dataMatrix[i].size(); j++) printf("%6.3f|", dataMatrix[i][j]);
printf("\n");
}
printf("---------------------------------\n");
printf("angle min :%f\n", angleMin);
printf("angle max :%f\n", angleMax);
printf("angle step:%f\n", angleStep);
printf("---------------------------------\n");
//------- calculate summed cross section
vector<double> partialXsec(Ex.size());
for( int i = 0; i < dataMatrix.size() ; i++){
for( int j = 0; j < (dataMatrix[i]).size() ; j++ ){
//double theta = (angle[j] + angleStep/2.) * TMath::DegToRad();
double theta = (angle[i]) * TMath::DegToRad();
double dTheta = angleStep * TMath::DegToRad();
double phi = TMath::TwoPi();
partialXsec[j] += dataMatrix[i][j] * sin( theta ) * dTheta * phi ;
}
}
for(int i = 0; i < Ex.size(); i++){
printf("Ex=%6.3f| Xsec(%3.0f-%3.0f deg) : %f mb\n", Ex[i] , angleMin, angleMax, partialXsec[i]);
}
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___SF____sigma\n");
for( int i = 0; i < Ex.size() ; 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 in ROOT
string saveFileName = readFile;
len = saveFileName.length();
saveFileName = saveFileName.substr(0, len - 4);
TString fileName = saveFileName;
fileName += ".root";
printf("Output : %s \n", fileName.Data());
TFile * fileOut = new TFile(fileName, "RECREATE" );
gList = new TObjArray();
gList->SetName("TGraph of d.s.c");
TObjArray * fList = new TObjArray();
fList->SetName("TF1 of distributions = d.s.c. * sin()");
TGraph ** gGraph = new TGraph *[numCal];
TF1 ** dist = new TF1*[numCal];
for( int i = 0; i < Ex.size() ; i++){
gGraph[i] = new TGraph();
TString name ;
name.Form("Ex=%6.3f MeV", Ex[i]);
gGraph[i]->SetName(name);
for( int j = 0; j < angle.size() ; j++){
gGraph[i]->SetPoint(j, angle[j], dataMatrix[j][i]);
}
gList->Add(gGraph[i]);
name.Form("dist%d", i);
dist[i] = new TF1(name, distfunct, angleMin, angleMax, 1 );
dist[i]->SetParameter(0, i);
fList->Add(dist[i]);
//delete tempFunc;
}
gList->Write("qList", 1);
fList->Write("pList", 1);
fileOut->Write();
fileOut->Close();
printf("---------------------------------------------------\n");
return 0;
}

41
Cleopatra/ExtractXSecFromText.C Normal file
View File

@ -0,0 +1,41 @@
/***********************************************************************
*
* This is ExtractXSecFromText for *.txt file
*
* This program will extract cross section distribution
* save as *.root for distribution
*
* save *.Ex.txt for excitation energies and total X-section
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ ExtractXSecFromText.C -o ExtractXSecFromText `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Dec-29, 2019
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h>
#include "ExtractXSec.h"
using namespace std;
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("========== Extract Cross-Section From text file ==========\n");
printf("=================================================================\n");
if(argc < 2 ) {
printf("Usage: ./ExtractXSecFromText input_file\n");
exit(0);
}else{
printf("From file : %s \n", argv[1]);
}
string readFile = argv[1];
ExtractXSecFromText(readFile);
}

68
Cleopatra/InFileCreator.C Normal file
View File

@ -0,0 +1,68 @@
/***********************************************************************
*
* This is InFileCreator, To creator the in-file for Ptolemy
* only for (d,p), (d,p), (d,d), or (p,p)
*
* It read a simple infile.in from reaction_setting file
*
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ InFileCreator.C -o InFileCreator `root-config --cflags --glibs`
*
*------------------------------------------------------
* The reaction_setting file is simple like:
*
* 206Hg(d,p)207Hg(1s1/2 0.000) 10MeV/u AK
*
* the first is similar to usual reaction setting, the word AK is a
* short name for Optical Potential, user can put as many line as
* they like, Cleopatra can create the suitable infile.in for Ptolomy
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h> /* atof */
#include <cmath>
#include <vector>
#include "InFileCreator.h"
using namespace std;
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("=== InfileCreator, Ptolemy for p,d,t,3He ====\n");
printf("=================================================================\n");
if(argc < 2 || argc > 5) {
printf("Usage: ./InfileCreator input_file (angMin = 0 deg, angMax = 180 deg, angStep = 1 deg)\n");
printf("Usage: ./InfileCreator input_file angMin angMax (angStep = 1 deg)\n");
printf("Usage: ./InfileCreator input_file angMin angMax angStep\n");
exit(0);
}else{
printf("From file : %s \n", argv[1]);
}
//================= read infile. extract the reactions, write pptolemy infile for each reaction
string readFile = argv[1];
double angMin = 0.;
double angMax = 180.;
double angStep = 1.;
if( argc >= 4 ){
angMin = atof(argv[2]);
angMax = atof(argv[3]);
}
if( argc == 5 ){
angStep = atof(argv[4]);
}
string ptolemyInFileName = argv[1];
ptolemyInFileName += ".in";
InFileCreator( readFile, ptolemyInFileName, angMin, angMax, angStep);
}

393
Cleopatra/InFileCreator.h Normal file
View File

@ -0,0 +1,393 @@
/***********************************************************************
*
* This is InFileCreator, To creator the in-file for Ptolemy
* only for (x,y), x or y = n, p, d, t, 3He
*
* It read a simple infile.in from reaction_setting file
*
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ InFileCreator.C -o InFileCreator `root-config --cflags --glibs`
*
*------------------------------------------------------
* The reaction_setting file is simple like:
*
* 206Hg(d,p)207Hg(1s1/2 0.000) 10MeV/u AK
*
* the first is similar to usual reaction setting, the word AK is a
* short name for Optical Potential, user can put as many line as
* they like, Cleopatra can create the suitable infile.in for Ptolomy
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h> /* atof */
#include <cmath>
#include <vector>
#include "../Cleopatra/Isotope.h" // for geting Z
#include "potentials.h"
#include "../Armory/AnalysisLibrary.h"
using namespace std;
int GetLValue(string spdf){
if( spdf == "s" ) return 0;
if( spdf == "p" ) return 1;
if( spdf == "d" ) return 2;
if( spdf == "f" ) return 3;
if( spdf == "g" ) return 4;
if( spdf == "h" ) return 5;
if( spdf == "i" ) return 6;
if( spdf == "j" ) return 7;
return -1;
}
int InFileCreator(string readFile, string infile, double angMin, double angMax, double angStep) {
//================= read infile. extract the reactions, write pptolemy infile for each reaction
ifstream file_in;
file_in.open(readFile.c_str(), ios::in);
if( !file_in ){
printf(" cannot read file. \n");
return 0 ;
}
printf("Save to infile : %s \n", infile.c_str());
FILE * file_out;
file_out = fopen (infile.c_str(), "w+");
printf("Angle setting (%5.2f, %5.2f) deg | Step : %5.2f deg\n", angMin, angMax, angStep);
printf("---------------------------\n");
//extract information
int numOfReaction = 0;
while( file_in.good() ) {
string tempLine;
getline(file_in, tempLine );
if( tempLine.substr(0, 1) == "#" ) continue;
if( tempLine.size() < 5 ) continue;
//split line using space
vector<string> str0 = SplitStr(tempLine, " ");
if ( str0.size() == 0 ) continue;
printf(" %s\n", tempLine.c_str());
///for( int i = 0 ; i < str0.size() ; i++){
/// printf(" str0[%d] %s \n", i, str0[i].c_str());
///}
vector<string> str1 = SplitStr(str0[0], "(", 0);
vector<string> str2 = SplitStr(str1[1], ")", 1);
str2[0] = "(" + str2[0];
//TODO use mass table for d, p, t, 3He
int lenStr20 = str2[0].length();
size_t posTemp1 = str2[0].find(",");
string ma = str2[0].substr(1, posTemp1-1);
size_t posTemp2 = str2[0].find(")");
string mb = str2[0].substr(posTemp1+1, posTemp2-posTemp1-1);
///printf(" ma : |%s| , mb : |%s| \n", ma.c_str(), mb.c_str());
Isotope isotopea(ma);
Isotope isotopeb(mb);
bool isReactionSupported = false;
if( isotopea.A <= 4 && isotopea.Z <= 2 && isotopeb.A <=4 && isotopeb.Z <=2 ) isReactionSupported = true;
///if( ma == "p" || ma == "d" || ma == "t" || ma == "3He" || mb == "n"){
/// if( mb == "p" || mb == "d" || mb == "t" || mb == "3He" || mb == "n" ) isReactionSupported = true;
///}
if( isReactionSupported == false ){
printf(" ===> Ignored. Reaction type not supported. \n");
continue;
}
string gsSpinA = str0[1];
string orbital = str0[2];
string spinParity = str0[3];
int lenSpinParity = spinParity.length();
string spin = spinParity.substr(0, lenSpinParity-1);
string parity = spinParity.substr(lenSpinParity-1);
string Ex = str0[4];
string reactionEnergy = str0[5];
string potential = str0[6];
string isoA = str1[0];
string isoB = str2[1];
string reactionType = str2[0];
Isotope isotopeA(str1[0]);
Isotope isotopeB(str2[1]);
///check reaction valid by balancing the A and Z number;
///printf("A: %d + %d = %d + %d \n", isotopeA.A, isotopea.A, isotopeB.A, isotopeb.A);
///printf("Z: %d + %d = %d + %d \n", isotopeA.Z, isotopea.Z, isotopeB.Z, isotopeb.Z);
if( isotopeA.A + isotopea.A != isotopeB.A + isotopeb.A || isotopeA.Z + isotopea.Z != isotopeB.Z + isotopeb.Z ) {
printf(" ====> ERROR! A-number or Z-number not balanced. \n");
Isotope isotopeK(isotopeA.A + isotopea.A - isotopeb.A, isotopeA.Z + isotopea.Z - isotopeb.Z);
printf(" try : %s ??\n", isotopeK.Name.c_str());
continue;
}
bool isTransferReaction = true;
if( ma == mb ) isTransferReaction = false;
if( isTransferReaction && potential.length() != 2 ){
printf(" ===> ERROR! Potential input should be 2 charaters! skipped. \n");
continue;
}
string node ;
string jValue ;
string lValue ;
int spdf = -1;
if( isTransferReaction ) {
///printf("------------ %d nucleon(s) transfer \n", abs(isotopea.A - isotopeb.A));
node = orbital.substr(0,1);
// single nucleon transfer
if( abs(isotopea.A - isotopeb.A) == 1 ){
lValue = orbital.substr(1,1);
jValue = orbital.substr(2);
///printf(" l : %s, j : %s \n", lValue.c_str(), jValue.c_str());
spdf = GetLValue(lValue);
}
// two-nucleons transfer
if( abs(isotopea.A - isotopeb.A) == 2 ){
size_t posEq = orbital.find('=');
lValue = orbital.substr(posEq+1,1);
spdf=atoi(lValue.c_str());
}
if( abs(isotopea.A - isotopeb.A) == 0 ){
printf(" ===? skipped. p-n exchange reaction does not support. \n");
}
if( spdf == -1 ){
printf(" ===> skipped. Not reconginzed orbital-label. (user input : l=%s | %s) \n", lValue.c_str(), orbital.c_str());
continue;
}
}
//get Beam energy, distingusih MeV or MeV/u
int pos = reactionEnergy.length() - 1;
for( int i = pos; i >= 0 ; i--){
if( isdigit(reactionEnergy[i]) ) {
pos = i;
break;
}
}
string unit = reactionEnergy.substr(pos+1);
//string beamParticle = reactionType.substr(1,1);
int factor = 1;
if( unit == "MeV/u") factor = isotopea.A;
double totalBeamEnergy = atof(reactionEnergy.substr(0, pos+1).c_str()) * factor;
///printf("unit : |%s| , %f\n", unit.c_str(), totalBeamEnergy);
///printf(" target nucleus : %s \n", isoA.c_str());
///printf(" reaction : %s \n", reactionType.c_str());
///printf(" remain : %s \n", isoB.c_str());
///printf(" reaction energy : %s \n", reactionEnergy.c_str());
///printf(" Potential : %s \n", potential.c_str());
///printf(" orbital : %s \n", orbital.c_str());
///printf(" Ex [MeV] : %s \n", Ex.c_str());
double Qvalue = isotopea.Mass + isotopeA.Mass - isotopeb.Mass - isotopeB.Mass;
///printf("Q-Value = %f MeV\n", Qvalue);
///printf("A: %d, Z: %d, mass: %f MeV/c2 \n", isotopeA.A, isotopeA.Z, isotopeA.Mass);
///printf("A: %d, Z: %d, mass: %f MeV/c2 \n", isotopeB.A, isotopeB.Z, isotopeB.Mass);
if( isotopeA.Mass == -404 || isotopeB.Mass == -404 ){
printf(" ===> Error! mass does not found. \n");
continue;
}
//write ptolmey infile
numOfReaction ++ ;
if( isTransferReaction ){
fprintf(file_out, "$============================================ Ex=%s(%s)%s\n", Ex.c_str(), orbital.c_str(), potential.c_str());
fprintf(file_out, "reset\n");
fprintf(file_out, "REACTION: %s%s%s(%s%s %s) ELAB=%7.3f\n",
isoA.c_str(), reactionType.c_str(), isoB.c_str(), spin.c_str(), parity.c_str(), Ex.c_str(), totalBeamEnergy);
if( isotopea.A <= 2 && isotopea.Z <= 1){ // incoming d or p
fprintf(file_out, "PARAMETERSET dpsb r0target \n");
fprintf(file_out, "$lstep=1 lmin=0 lmax=30 maxlextrap=0 asymptopia=50 \n");
fprintf(file_out, "\n");
fprintf(file_out, "PROJECTILE \n");
fprintf(file_out, "wavefunction av18 \n");
fprintf(file_out, "r0=1 a=0.5 l=0 rc0=1.2\n");
}
if( 3 <= isotopea.A && isotopea.A <= 4 && abs(isotopea.A-isotopeb.A) == 1){
fprintf(file_out, "PARAMETERSET alpha3 r0target \n");
fprintf(file_out, "$lstep=1 lmin=0 lmax=30 maxlextrap=0 asymptopia=50 \n");
fprintf(file_out, "\n");
fprintf(file_out, "PROJECTILE \n");
fprintf(file_out, "wavefunction phiffer \n");
if( isotopeb.A <= 3 ){
fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.31 v=179.94 r=0.54 a=0.68 param1=0.64 param2=1.13 rc=2.0\n");
}
if( isotopeb.A == 4 ){
fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.61 v=202.21 r=.93 a=.66 param1=.81 param2=.87 rc=2.0 $ rc=2 is a quirk\n");
}
}
if( abs(isotopea.A-isotopeb.A) == 2 ){ // 2 neutron transfer
fprintf(file_out, "PARAMETERSET alpha3 r0target\n");
fprintf(file_out, "lstep=1 lmin=0 lmax=30 maxlextrap=0 ASYMPTOPIA=40\n");
fprintf(file_out, "\n");
fprintf(file_out, "PROJECTILE\n");
fprintf(file_out, "L = 0 NODES=0 R0 = 1.25 A = .65 RC0 = 1.25\n");
}
fprintf(file_out, ";\n");
//===== TARGET
fprintf(file_out, "TARGET\n");
///check Ex is above particle threshold
double nThreshold = isotopeB.CalSp(0,1);
double pThreshold = isotopeB.CalSp(1,0);
bool isAboveThreshold = false;
double ExEnergy = atof(Ex.c_str());
if( ExEnergy > nThreshold || ExEnergy > pThreshold ) {
isAboveThreshold = true;
printf(" Ex = %.3f MeV is above thresholds; Sp = %.3f MeV, Sn = %.3f MeV\n", ExEnergy, pThreshold, nThreshold);
}
if( abs(isotopea.A-isotopeb.A) == 1 ){
fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str());
if( isAboveThreshold ) {
fprintf(file_out, "nodes=%s l=%d jp=%s E=-.2 $node is n-1, set binding 200 keV \n", node.c_str(), spdf, jValue.c_str());
}else{
fprintf(file_out, "nodes=%s l=%d jp=%s $node is n-1 \n", node.c_str(), spdf, jValue.c_str());
}
fprintf(file_out, "r0=1.25 a=.65 \n");
fprintf(file_out, "vso=6 rso0=1.10 aso=.65 \n");
fprintf(file_out, "rc0=1.3 \n");
}
if( abs(isotopea.A-isotopeb.A) == 2 ){
fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str());
if( isAboveThreshold ){
fprintf(file_out, "nodes=%s L=%d E=-.2 $node is n-1, binding by 200 keV \n", node.c_str(), spdf);
}else{
fprintf(file_out, "nodes=%s L=%d $node is n-1 \n", node.c_str(), spdf);
}
}
fprintf(file_out, ";\n");
//===== POTENTIAL
string pot1Name = potential.substr(0,1);
string pot1Ref = potentialRef(pot1Name);
fprintf(file_out, "INCOMING $%s\n", pot1Ref.c_str());
CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy, isotopea.Z);
fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a);
fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai);
fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi);
fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso);
fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0);
fprintf(file_out, ";\n");
string pot2Name = potential.substr(1,1);
string pot2Ref = potentialRef(pot2Name);
fprintf(file_out, "OUTGOING $%s\n", pot2Ref.c_str());
//printf(" total Beam Energy : %f | Qvalue : %f | Ex : %f \n", totalBeamEnergy, Qvalue, atof(Ex.c_str()));
double eBeam = totalBeamEnergy + Qvalue - atof(Ex.c_str());
CallPotential(pot2Name, isotopeB.A, isotopeB.Z, eBeam, isotopeb.Z);
fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a);
fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai);
fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi);
fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso);
fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0);
fprintf(file_out, ";\n");
}
if( isTransferReaction == false ){
if ( atof(Ex.c_str()) == 0.0 ) {
fprintf(file_out, "$============================================ ELab=%5.2f(%s+%s)%s\n",
totalBeamEnergy, ma.c_str(), isoA.c_str(), potential.c_str());
fprintf(file_out, "reset\n");
fprintf(file_out, "CHANNEL %s + %s\n", ma.c_str(), isoA.c_str());
fprintf(file_out, "ELAB = %f\n", totalBeamEnergy);
fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str());
string pot1Name = potential.substr(0,1);
string pot1Ref = potentialRef(pot1Name);
fprintf(file_out, "$%s\n", pot1Ref.c_str());
CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy, isotopea.Z);
fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a);
fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai);
fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi);
fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso);
fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0);
fprintf(file_out, "ELASTIC SCATTERING\n");
fprintf(file_out, ";\n");
}else{
fprintf(file_out, "$============================================ Ex=%s(%s+%s|%s%s)%s,ELab=%5.2f\n",
Ex.c_str(), ma.c_str(), isoA.c_str(), spin.c_str(), parity.c_str(), potential.c_str(),totalBeamEnergy);
fprintf(file_out, "reset\n");
fprintf(file_out, "REACTION: %s%s%s(%s%s %s) ELAB=%7.3f\n",
isoA.c_str(), reactionType.c_str(), isoB.c_str(), spin.c_str(), parity.c_str(), Ex.c_str(), totalBeamEnergy);
fprintf(file_out, "PARAMETERSET ineloca2 r0target\n");
fprintf(file_out, "JBIGA=%s\n", gsSpinA.c_str());
if( str0.size() >= 8 ){
fprintf(file_out, "BETA=%s\n", str0[7].c_str()); //deformation length
}
string pot1Name = potential.substr(0,1);
string pot1Ref = potentialRef(pot1Name);
fprintf(file_out, "$%s\n", pot1Ref.c_str());
CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy, isotopea.Z);
fprintf(file_out, "INCOMING\n");
fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a);
fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai);
fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f rc0 = %7.3f\n", vsi, rsi0, asi, rc0);
//fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso);
//fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0);
fprintf(file_out, ";\n");
fprintf(file_out, "OUTGOING\n");
fprintf(file_out, "$%s\n", pot1Ref.c_str());
CallPotential(pot1Name, isotopeA.A, isotopeA.Z, totalBeamEnergy - atof(Ex.c_str()), isotopea.Z);
fprintf(file_out, "v = %7.3f r0 = %7.3f a = %7.3f\n", v, r0, a);
fprintf(file_out, "vi = %7.3f ri0 = %7.3f ai = %7.3f\n", vi, ri0, ai);
fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f rc0 = %7.3f\n", vsi, rsi0, asi, rc0);
//fprintf(file_out, "vsi = %7.3f rsi0 = %7.3f asi = %7.3f\n", vsi, rsi0, asi);
//fprintf(file_out, "vso = %7.3f rso0 = %7.3f aso = %7.3f\n", vso, rso0, aso);
//fprintf(file_out, "vsoi = %7.3f rsoi0 = %7.3f asoi = %7.3f rc0 = %7.3f\n", vsoi, rsoi0, asoi, rc0);
fprintf(file_out, ";\n");
}
}
fprintf(file_out, "anglemin=%f anglemax=%f anglestep=%f\n", angMin, angMax, angStep);
fprintf(file_out, ";\n");
}
printf("================= end of input. Number of Reaction : %d \n", numOfReaction);
fprintf(file_out, "end $================================== end of input\n");
file_in.close();
fclose(file_out);
return 1;
}

64
Cleopatra/Isotope.C Normal file
View File

@ -0,0 +1,64 @@
/***********************************************************************
*
* This is Isotope.C for isotope mass-related quatilies
*
* -----------------------------------------------------
* To compile
*
* g++ Isotope.C -o Isotope
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Feb-20, 2021
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h>
#include "Isotope.h"
using namespace std;
void Usage(){
cout << "./Isotope Sym" << endl;
cout << "./Isotope A Z" << endl;
///cout << "./Isotope A Z A' Z'" << endl;
///cout << " Opt = rkm, e.g. 001 only calculate mass" << endl;
///cout << " |||_ mass " << endl;
///cout << " ||__ kinematics " << endl;
///cout << " |___ reaction kinematics " << endl;
exit(0);
}
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("========== Isotope mass-related quantities ==========\n");
printf("=================================================================\n");
if ( argc != 2 && argc != 3 && argc != 6) Usage();
Isotope iso1, iso2;
int Z, A, Za, Aa;
//string Opt = argv[1];
if (argc == 2){
iso1.SetIsoByName(argv[1]);
}
if (argc == 3){
A= atoi(argv[1]);
Z= atoi(argv[2]);
iso1.SetIso(A, Z);
//}else if ( argc == 6){
// A= atoi(argv[2]);
// Z= atoi(argv[3]);
// Aa= atoi(argv[4]);
// Za= atoi(argv[5]);
// iso1.SetIso(A, Z);
// iso2.SetIso(Aa,Za);
}
iso1.Print();
iso1.ListShell();
}

536
Cleopatra/Isotope.h Normal file
View File

@ -0,0 +1,536 @@
/***********************************************************************
*
* 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="/opt/Ptolemy/Cleopatra/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(){findHeliosPath();};
Isotope(int a, int z){ findHeliosPath();SetIso(a,z); };
Isotope(string name){ findHeliosPath(); 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 findHeliosPath(){
dataSource = massData;
// heliosPath = getenv("HELIOSSYS");
// if( heliosPath ){
// dataSource = heliosPath;
// dataSource += "/analysis" + massData;
// }else{
// dataSource = ".." + massData;
// }
}
};
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){
findHeliosPath();
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

50
Cleopatra/PlotTGraphTObjArray.C Normal file
View File

@ -0,0 +1,50 @@
/***********************************************************************
*
* This is PlotResultInRoot.C for ExtractXSec *.root output
*
* The Xsec are stored in (TObjArray *) gList
*
* This program is simple get plot all the member in the gList
*
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ PlotResultInROOT.C -o PlotResultInROOT `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h>
#include <TApplication.h>
#include "PlotTGraphTObjArray.h"
using namespace std;
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("==================== Plot Results in ROOT =======================\n");
printf("=================================================================\n");
if(argc != 2) {
printf("Usage: ./PlotTGraphTObjArray root_file\n");
exit(0);
}else{
printf("From file : %s \n", argv[1]);
printf("=========== Press Ctrl+C to end.\n");
}
string readFile = argv[1];
TApplication app ("app", &argc, argv);
PlotTGraphTObjArray(readFile);
app.Run(); //anything after this line is not running
return 0;
}

88
Cleopatra/PlotTGraphTObjArray.h Normal file
View File

@ -0,0 +1,88 @@
/***********************************************************************
*
* This is PlotResultInRoot.C for ExtractXSec *.root output
*
* The Xsec are stored in (TObjArray *) gList
*
* This program is simple get plot all the member in the gList
*
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ PlotResultInROOT.C -o PlotResultInROOT `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <TROOT.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TObjArray.h>
#include <TGraph.h>
#include <TF1.h>
#include <TAxis.h>
#include <TH1F.h>
#include <TLegend.h>
void PlotTGraphTObjArray(TString rootFileName){
TFile * file = new TFile(rootFileName, "READ");
TObjArray * gList = (TObjArray *) file->FindObjectAny("qList");
if( gList == NULL ) {
printf("No Result was found.\n");
return;
}
TCanvas * cPlots = new TCanvas("cPlots", "Ptolemy Results", 0, 0, 800, 600);
cPlots->SetLogy();
TLegend * legend = new TLegend( 0.6, 0.2, 0.9, 0.4);
const int n = gList->GetLast() + 1 ;
TGraph * gr[n];
//Get minimum, maximum Y
double maxY = 0;
double minY = 10000000;
for ( int i = 0; i < n ; i++){
gr[i] = (TGraph *) gList->At(i);
gr[i]->SetLineColor(i+1);
gr[i]->GetXaxis()->SetTitle("#theta_{CM} [deg]");
gr[i]->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/sr]");
TString name = gr[i]->GetName();
int pos = name.First("|");
name.Remove(0, pos+1);
legend->AddEntry(gr[i], name);
double miny = gr[i]->GetHistogram()->GetMinimum();
double maxy = gr[i]->GetHistogram()->GetMaximum();
if( miny < minY ) minY = miny;
if( maxy > maxY ) maxY = maxy;
}
for ( int i = 0; i < n ; i++){
gr[i]->Draw("same");
if( i == 0 ){
gr[i]->Draw();
gr[i]->GetYaxis()->SetRangeUser(minY * 0.8, maxY * 1.2);
}else{
gr[i]->Draw("same");
}
}
legend->Draw();
cPlots->Update();
cPlots->Draw();
}

71
Cleopatra/constant.h Normal file
View File

@ -0,0 +1,71 @@
/***********************************************************************
*
* 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;
#endif

27
Cleopatra/makefile Normal file
View File

@ -0,0 +1,27 @@
CC=g++
all: Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray FindThetaCM Transfer
#Cleopatra: Cleopatra.C ../Simulation/Isotope.h ../Simulation/constant.h potentials.h InFileCreator.h ExtractXSec.h PlotTGraphTObjArray.h
# $(CC) Cleopatra.C -o Cleopatra `root-config --cflags --glibs`
InFileCreator: InFileCreator.C InFileCreator.h ../Cleopatra/Isotope.h ../Cleopatra/constant.h potentials.h
$(CC) InFileCreator.C -o InFileCreator `root-config --cflags --glibs`
ExtractXSec: ExtractXSec.C ExtractXSec.h
$(CC) ExtractXSec.C -o ExtractXSec `root-config --cflags --glibs`
ExtractXSecFromText: ExtractXSecFromText.C ExtractXSec.h
$(CC) ExtractXSecFromText.C -o ExtractXSecFromText `root-config --cflags --glibs`
PlotTGraphTObjArray: PlotTGraphTObjArray.C PlotTGraphTObjArray.h
$(CC) PlotTGraphTObjArray.C -o PlotTGraphTObjArray `root-config --cflags --glibs`
FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/HELIOS_LIB.h ../Cleopatra/Isotope.h ../Cleopatra/constant.h
$(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs`
Transfer: Transfer.C Transfer.h ../Cleopatra/HELIOS_LIB.h ../Cleopatra/Isotope.h ../Cleopatra/constant.h
$(CC) Transfer.C -o Transfer `root-config --cflags --glibs`
Isotope: ../Cleopatra/Isotope.h ../Cleopatra/Isotope.C
$(CC) Isotope.C -o Isotope

3475
Cleopatra/mass16.txt Normal file
View File

File diff suppressed because it is too large Load Diff

3594
Cleopatra/mass20.txt Normal file
View File

File diff suppressed because it is too large Load Diff

1073
Cleopatra/potentials.h Normal file
View File

File diff suppressed because it is too large Load Diff

BIN
Cleopatra/ptolemy Executable file
View File

Binary file not shown.

51
DWBA Normal file
View File

@ -0,0 +1,51 @@
#========= Input for Cleopatra
#===== # for comment line, must be at the beginning of line
#===== the potential contain two words
# one for incoming
# one for outgoing
#================================================= Potenital abberation
#========================= deuteron
# A = An, Cai, 2006 | E < 183 | 12 < A < 238 | http://dx.doi.org/10.1103/PhysRevC.73.054605
# H = Han, Shi, Shen, 2006 | E < 200 | 12 < A < 209 | http://dx.doi.org/10.1103/PhysRevC.74.044615
# B = Bojowald et al., 1988 | 50 < E < 80 | 27 < A < 208 | http://dx.doi.org/10.1103/PhysRevC.38.1153
# D = Daehnick, Childs, Vrcelj, 1980 | 11.8 < E < 80 | 27 < A < 238 (REL) | http://dx.doi.org/10.1103/PhysRevC.21.2253
# C = Daehnick, Childs, Vrcelj, 1980 | 11.8 < E < 80 | 27 < A < 238 (NON-REL) | http://dx.doi.org/10.1103/PhysRevC.21.2253 // not impletmented yet
# L = Lohr and Haeberli, 1974 | 9 < E < 13 | 40 < A | http://dx.doi.org/10.1016/0375-9474(74)90627-7
# Q = Perey and Perey, 1963 | 12 < E < 25 | 40 < A | http://dx.doi.org/10.1016/0370-1573(91)90039-O
# Z = Zhang, Pang, Lou, 2016 | 5 < E < 170 | A < 18, spe 6-7Li | https://doi.org/10.1103/PhysRevC.94.014619
#========================= proton
# K = Koning and Delaroche, 2009 | 0.001 < E < 200 | 24 < A < 209 | Iso. Dep. | http://dx.doi.org/10.1016/S0375-9474(02)01321-0
# V = Varner et al., (CH89), 1991 | 16 < E < 65 | 4 < A < 209 | http://dx.doi.org/10.1016/0370-1573(91)90039-O
# M = Menet et al., 1971 | 30 < E < 60 | 40 < A | http://dx/doi.org/10.1016/0092-640X(76)90007-3
# G = Becchetti and Greenlees, 1969 | E < 50 | 40 < A | http://dx.doi.org/10.1103/PhysRev.182.1190
# P = Perey, 1963 | E < 20 | 30 < A < 100 | http://dx/doi.org/10.1016/0092-640X(76)90007-3
#========================= A=3
# x = Xu, Guo, Han, Shen, 2011 | E < 250 | 20 < A < 209 | 3He | http://dx.doi.org/10.1007/s11433-011-4488-5
# l = Liang, Li, Cai, 2009 | E < 270 | All masses | http://dx.doi.org/10.1088/0954-3899/36/8/085104
# p = Pang et al., 2009 | All E | All masses | Iso. Dep. | http://dx.doi.org/10.1103/PhysRevC.79.024615
# c = Li, Liang, Cai, 2007 | E < 40 | 48 < A < 232 | Tritons | http://dx.doi.org/10.1016/j.nuclphysa.2007.03.004
# t = Trost et al., 1987 | 10 < E < 220 | 10 < A < 208 | http://dx.doi.org/10.1016/0375-9474(87)90551-3
# h = Hyakutake et al., 1980 | 90 < E < 120 | About 58 < A < 92 | http://dx.doi.org/10.1016/0375-9474(80)90013-5
# b = Becchetti and Greenlees, 1971 | E < 40 | 40 < A | Iso. Dep.
#========================= alpha
# s = Su and Han, 2015 | E < 398 | 20 < A < 209 | http://dx.doi/org/10.1142/S0218301315500925
# a = Avrigeanu et al., 2009 | E ??? | A ??? | http://dx.doi/org/10.1016/j.adt.2009.02.001
# f = Bassani and Picard, 1969(FIXED)| 24 < E < 31 | A = 90 | https://doi.org/10.1016/0375-9474(69)90601-0
#=======================================================================
#reaction gs-spin orbital spin-pi(Ex) Ex ELab Potentials
#206Hg(d,d)206Hg 0 none 9/2+ 0.000 7.39MeV/u AA #elastic
#206Hg(d,d)206Hg 0 none 9/2+ 1.000 7.39MeV/u AA 0.12 #inelastics_0.12=beta
#206Hg(d,p)207Hg 0 1g9/2 9/2+ 0.000 7.39MeV/u AK
#20F(d,t)19F 2 0d5/2 5/2+ 0.197 10MeV/u Vl
#16N(d,3He)15C 2 0p1/2 5/2+ 0.74 12MeV/u Ax
#10Be(t,p)12Be 0 1L=0 0+ 0.000 5MeV/u lA #two-nucleon_transfer
#32Si(t,p)34Si 0 0L=0 0+ 0.000 8MeV/u lA #two-nucleon_transfer
#133Sb(t,3He)133Sn 7/2 0g7/2 0+ 0.000 8.5MeV/u Ax .... cannot cal
#82Kr(d,p)83Kr 0 0g9/2 9/2+ 0.000 10MeV/u AK
#82Kr(d,p)83Kr 0 1p1/2 1/2- 0.000 10MeV/u AK
#82Kr(d,p)83Kr 0 1d3/2 3/2+ 1.220 10MeV/u AK
#82Kr(d,p)83Kr 0 2s1/2 1/2+ 1.220 10MeV/u AK
10Be(d,d)10Be 0 none 0+ 0.000 9MeV/u ZZ
10Be(d,d)10Be 0 none 0+ 0.000 9MeV/u AA

9
working/PtolemyGUI Executable file
View File

@ -0,0 +1,9 @@
#!/usr/bin/bash
if [[ -f DWBA ]]; then
dummy=10
else
cp /opt/Ptolemy/DWBA .
fi
root -l /opt/Ptolemy/working/PtolemyGUI.C

544
working/PtolemyGUI.C Normal file
View File

@ -0,0 +1,544 @@
#include <TGClient.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TRandom.h>
#include <TGButton.h>
#include <TGLabel.h>
#include <TGFrame.h>
#include <TGTextEditor.h>
#include <TGNumberEntry.h>
#include <TGComboBox.h>
#include <TRootEmbeddedCanvas.h>
#include <RQ_OBJECT.h>
//#include "../Cleopatra/Transfer.h"
#include "../Cleopatra/InFileCreator.h"
#include "../Cleopatra/ExtractXSec.h"
#include "../Cleopatra/PlotTGraphTObjArray.h"
//#include "../Armory/AutoFit.C"
//#include "../Armory/Check_Simulation.C"
#include <iostream>
#include <stdexcept>
#include <stdio.h>
#include <string>
#ifdef __linux__
#define OS_Type 1
#elif __APPLE__
#define OS_Type 0
#endif
class MyMainFrame {
RQ_OBJECT("MyMainFrame")
private:
TGMainFrame *fMain;
TGTextEdit * editor;
TString fileName;
TGLabel * fileLabel;
TGLabel * statusLabel;
TGNumberEntry * angMin;
TGNumberEntry * angMax;
TGNumberEntry * angStep;
TGCheckButton * withDWBA;
TGCheckButton * isInFile;
TGCheckButton * isRun;
TGCheckButton * isExtract;
TGCheckButton * isPlot;
TGComboBox * extractFlag;
public:
MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h);
virtual ~MyMainFrame();
void Command(int);
void OpenFile(int);
bool IsFileExist(TString filename);
};
MyMainFrame::MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h) {
// Create a main frame
fMain = new TGMainFrame(p,w,h);
TGHorizontalFrame *hframe = new TGHorizontalFrame(fMain,600,600 );
fMain->AddFrame(hframe, new TGLayoutHints(kLHintsCenterX, 2,2,2,2));
TGVerticalFrame *hframe1 = new TGVerticalFrame(fMain,600,600 );
hframe->AddFrame(hframe1);
TGVerticalFrame *hframe2 = new TGVerticalFrame(fMain,600,800 );
hframe->AddFrame(hframe2);
fileName = "DWBA";
TGHorizontalFrame *hframe00 = new TGHorizontalFrame(hframe2,600,600 );
hframe2->AddFrame(hframe00, new TGLayoutHints(kLHintsCenterX, 2,2,2,2));
fileLabel = new TGLabel(hframe00, "");
fileLabel->SetWidth(370);
fileLabel->SetHeight(20);
fileLabel->SetTextColor(kRed);
fileLabel->ChangeOptions(kFixedSize | kSunkenFrame);
fileLabel->SetText(fileName);
hframe00->AddFrame(fileLabel, new TGLayoutHints(kLHintsLeft, 2,2,2,2));
TGTextButton *save = new TGTextButton(hframe00,"Save");
save->SetWidth(100);
save->SetHeight(20);
save->ChangeOptions( save->GetOptions() | kFixedSize );
save->Connect("Clicked()","MyMainFrame",this,"Command(=3)");
hframe00->AddFrame(save, new TGLayoutHints(kLHintsLeft,5,5,3,4));
TGTextButton *help = new TGTextButton(hframe00, "Help");
help->SetWidth(100);
help->SetHeight(20);
help->ChangeOptions( help->GetOptions() | kFixedSize );
help->Connect("Clicked()","MyMainFrame",this,"Command(=4)");
hframe00->AddFrame(help,new TGLayoutHints(kLHintsLeft, 5,5,3,4));
editor = new TGTextEdit(hframe2, 600, 700);
editor->LoadFile(fileName);
hframe2->AddFrame(editor);
statusLabel = new TGLabel(hframe2, "");
statusLabel->SetWidth(600);
statusLabel->SetHeight(20);
statusLabel->SetTextJustify(kTextLeft);
statusLabel->SetTextColor(1);
statusLabel->ChangeOptions(kFixedSize | kSunkenFrame);
hframe2->AddFrame(statusLabel, new TGLayoutHints(kLHintsLeft, 2,2,2,2));
//================= Simulation group
/*
TGGroupFrame * simFrame = new TGGroupFrame(fMain, "Kinematics Simulation", kVerticalFrame);
hframe1->AddFrame(simFrame, new TGLayoutHints(kLHintsCenterX, 5,5,3,4));
TGTextButton *openRec = new TGTextButton(simFrame, "reaction Config");
openRec->SetWidth(150);
openRec->SetHeight(20);
openRec->ChangeOptions( openRec->GetOptions() | kFixedSize );
openRec->Connect("Clicked()","MyMainFrame",this, "OpenFile(=1)");
simFrame->AddFrame(openRec,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openDet = new TGTextButton(simFrame, "detector Geo.");
openDet->SetWidth(150);
openDet->SetHeight(20);
openDet->ChangeOptions( openDet->GetOptions() | kFixedSize );
openDet->Connect("Clicked()","MyMainFrame",this, "OpenFile(=0)");
simFrame->AddFrame(openDet,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openEx = new TGTextButton(simFrame, "Ex List");
openEx->SetWidth(150);
openEx->SetHeight(20);
openEx->ChangeOptions( openEx->GetOptions() | kFixedSize );
openEx->Connect("Clicked()","MyMainFrame",this, "OpenFile(=2)");
simFrame->AddFrame(openEx,new TGLayoutHints(kLHintsRight, 5,5,3,4));
withDWBA = new TGCheckButton(simFrame, "Sim with DWBA\n+DWBA.root\n+DWBA.Ex.txt");
withDWBA->SetWidth(140);
withDWBA->ChangeOptions(kFixedSize );
simFrame->AddFrame(withDWBA, new TGLayoutHints(kLHintsLeft, 5,5,3,4));
TGTextButton *Sim = new TGTextButton(simFrame,"Simulate");
Sim->SetWidth(150);
Sim->SetHeight(40);
Sim->ChangeOptions( Sim->GetOptions() | kFixedSize );
Sim->Connect("Clicked()","MyMainFrame",this,"Command(=1)");
simFrame->AddFrame(Sim, new TGLayoutHints(kLHintsRight,5,5,3,4));
TGTextButton *openSimChk = new TGTextButton(simFrame, "Config Simulation Plot");
openSimChk->SetWidth(150);
openSimChk->SetHeight(20);
openSimChk->ChangeOptions( openSimChk->GetOptions() | kFixedSize );
openSimChk->Connect("Clicked()","MyMainFrame",this, "OpenFile(=4)");
simFrame->AddFrame(openSimChk,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *SimChk = new TGTextButton(simFrame,"Plot Simulation");
SimChk->SetWidth(150);
SimChk->SetHeight(40);
SimChk->ChangeOptions( SimChk->GetOptions() | kFixedSize );
SimChk->Connect("Clicked()","MyMainFrame",this,"Command(=2)");
simFrame->AddFrame(SimChk, new TGLayoutHints(kLHintsRight,5,5,3,4));
TGTextButton *autoFit = new TGTextButton(simFrame,"AutoFit ExCal");
autoFit->SetWidth(150);
autoFit->SetHeight(40);
autoFit->ChangeOptions( autoFit->GetOptions() | kFixedSize );
autoFit->Connect("Clicked()","MyMainFrame",this,"Command(=5)");
simFrame->AddFrame(autoFit, new TGLayoutHints(kLHintsRight,5,5,3,4));
*/
//================= DWBA group
TGGroupFrame * DWBAFrame = new TGGroupFrame(fMain, "DWBA calculation", kVerticalFrame);
hframe1->AddFrame(DWBAFrame, new TGLayoutHints(kLHintsCenterX, 5,5,3,4));
TGTextButton *openDWBA = new TGTextButton(DWBAFrame, "DWBA setting");
openDWBA->SetWidth(150);
openDWBA->SetHeight(20);
openDWBA->ChangeOptions( openDWBA->GetOptions() | kFixedSize );
openDWBA->Connect("Clicked()","MyMainFrame",this, "OpenFile(=3)");
DWBAFrame->AddFrame(openDWBA,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openInFile = new TGTextButton(DWBAFrame, "InFile");
openInFile->SetWidth(150);
openInFile->SetHeight(20);
openInFile->ChangeOptions( openDWBA->GetOptions() | kFixedSize );
openInFile->Connect("Clicked()","MyMainFrame",this, "OpenFile(=5)");
DWBAFrame->AddFrame(openInFile,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *openOutFile = new TGTextButton(DWBAFrame, "OutFile");
openOutFile->SetWidth(150);
openOutFile->SetHeight(20);
openOutFile->ChangeOptions( openDWBA->GetOptions() | kFixedSize );
openOutFile->Connect("Clicked()","MyMainFrame",this, "OpenFile(=6)");
DWBAFrame->AddFrame(openOutFile,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *xsecFile = new TGTextButton(DWBAFrame, "X-Sec");
xsecFile->SetWidth(150);
xsecFile->SetHeight(20);
xsecFile->ChangeOptions( openDWBA->GetOptions() | kFixedSize );
xsecFile->Connect("Clicked()","MyMainFrame",this, "OpenFile(=7)");
DWBAFrame->AddFrame(xsecFile,new TGLayoutHints(kLHintsRight, 5,5,3,4));
//-------- angle setting
TGHorizontalFrame * hframe000 = new TGHorizontalFrame(DWBAFrame, 150, 30, kFixedSize);
DWBAFrame->AddFrame(hframe000);
TGLabel * lb1 = new TGLabel(hframe000, "angMin");
lb1->SetWidth(50); lb1->ChangeOptions( kFixedSize);
hframe000->AddFrame(lb1, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 5, 5, 0, 0));
TGLabel * lb2 = new TGLabel(hframe000, "angMax");
lb2->SetWidth(50); lb2->ChangeOptions( kFixedSize);
hframe000->AddFrame(lb2, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 5, 5, 0, 0));
TGLabel * lb3 = new TGLabel(hframe000, "angStep");
lb3->SetWidth(50); lb3->ChangeOptions( kFixedSize);
hframe000->AddFrame(lb3, new TGLayoutHints(kLHintsCenterX | kLHintsCenterY, 5, 5, 0, 0));
TGHorizontalFrame * hframe001 = new TGHorizontalFrame(DWBAFrame, 150, 30, kFixedSize);
DWBAFrame->AddFrame(hframe001);
angMin = new TGNumberEntry(hframe001, 0, 0, 0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative);
angMin->SetWidth(50);
angMin->SetLimits(TGNumberFormat::kNELLimitMinMax, 0, 180);
hframe001->AddFrame(angMin, new TGLayoutHints(kLHintsCenterX , 5, 5, 0, 0));
angMax = new TGNumberEntry(hframe001, 60, 0, 0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative);
angMax->SetWidth(50);
angMax->SetLimits(TGNumberFormat::kNELLimitMinMax, 0, 180);
hframe001->AddFrame(angMax, new TGLayoutHints(kLHintsCenterX , 5, 5, 0, 0));
angStep = new TGNumberEntry(hframe001, 1, 0, 0, TGNumberFormat::kNESRealOne, TGNumberFormat::kNEAPositive);
angStep->SetWidth(50);
angStep->SetLimits(TGNumberFormat::kNELLimitMinMax, 0, 30);
hframe001->AddFrame(angStep, new TGLayoutHints(kLHintsCenterX, 5, 5, 0, 0));
//------- Check Boxes
isInFile = new TGCheckButton(DWBAFrame, "Create inFile");
isInFile->SetWidth(100);
isInFile->ChangeOptions(kFixedSize );
isInFile->SetState(kButtonDown);
DWBAFrame->AddFrame(isInFile, new TGLayoutHints(kLHintsLeft, 5,5,3,4));
isRun = new TGCheckButton(DWBAFrame, "Run Ptolemy");
isRun->SetWidth(100);
isRun->ChangeOptions(kFixedSize );
isRun->SetState(kButtonDown);
DWBAFrame->AddFrame(isRun, new TGLayoutHints(kLHintsLeft, 5,5,3,4));
isExtract = new TGCheckButton(DWBAFrame, "Extract Xsec");
isExtract->SetWidth(100);
isExtract->ChangeOptions(kFixedSize );
isExtract->SetState(kButtonDown);
DWBAFrame->AddFrame(isExtract, new TGLayoutHints(kLHintsLeft, 5,5,3,4));
isPlot = new TGCheckButton(DWBAFrame, "Plot");
isPlot->SetWidth(100);
isPlot->ChangeOptions(kFixedSize );
isPlot->SetState(kButtonDown);
DWBAFrame->AddFrame(isPlot, new TGLayoutHints(kLHintsLeft, 5,5,3,4));
extractFlag = new TGComboBox(DWBAFrame, 100);
extractFlag->SetWidth(130);
extractFlag->SetHeight(30);
extractFlag->AddEntry("Xsec.", 2);
extractFlag->AddEntry("Ratio to Ruth.", 1);
extractFlag->AddEntry("Rutherford", 3);
extractFlag->Select(2);
DWBAFrame->AddFrame(extractFlag, new TGLayoutHints(kLHintsLeft, 5,5,3,4));
TGTextButton *DWBA = new TGTextButton(DWBAFrame, "DWBA");
DWBA->SetWidth(150);
DWBA->SetHeight(40);
DWBA->ChangeOptions( DWBA->GetOptions() | kFixedSize );
DWBA->Connect("Clicked()","MyMainFrame",this,"Command(=0)");
DWBAFrame->AddFrame(DWBA,new TGLayoutHints(kLHintsRight, 5,5,3,4));
TGTextButton *exit = new TGTextButton(hframe1,"Exit", "gApplication->Terminate(0)");
exit->SetWidth(150);
exit->SetHeight(40);
exit->ChangeOptions( exit->GetOptions() | kFixedSize );
hframe1->AddFrame(exit, new TGLayoutHints(kLHintsCenterX, 5,5,3,4));
// Set a name to the main frame
fMain->SetWindowName("Simulation Helper");
// Map all subwindows of main frame
fMain->MapSubwindows();
// Initialize the layout algorithm
fMain->Resize(fMain->GetDefaultSize());
// Map main frame
fMain->MapWindow();
int versionInt = gROOT->GetVersionInt();
if( versionInt < 62600 ) {
statusLabel->SetText(Form("Root version : %s. Please Update Root to v6.26/00",gROOT->GetVersion()));
}else{
statusLabel->SetText(Form("Root version : %s",gROOT->GetVersion()));
}
}
bool MyMainFrame::IsFileExist(TString filename){
ifstream file (filename.Data());
return file.is_open();
}
void MyMainFrame::OpenFile(int ID){
editor->SaveFile(fileName);
TString oldFileName = fileName;
if ( ID == 0 ) fileName = "detectorGeo.txt";
if ( ID == 1 ) fileName = "reactionConfig.txt";
if ( ID == 2 ) fileName = "Ex.txt";
if ( ID == 3 ) fileName = "DWBA";
if ( ID == 4 ) fileName = "../Armory/Check_Simulation_Config.txt";
if ( ID == 5 ) fileName = "DWBA.in";
if ( ID == 6 ) fileName = "DWBA.out";
if ( ID == 7 ) fileName = "DWBA.Xsec.txt";
//test if file exist
if ( IsFileExist(fileName) ){
fileLabel->SetText(fileName);
editor->LoadFile(fileName);
if( ID >= 6 ) {
editor->SetReadOnly(true);
}else{
editor->SetReadOnly(false);
}
editor->ShowBottom();
if( ID < 6){
statusLabel->SetText(fileName + " opened.");
}else{
statusLabel->SetText(fileName + " opened. (READ ONLY)");
}
}else{
statusLabel->SetText(fileName + " not exist.");
fileName = oldFileName;
}
}
void MyMainFrame::Command(int ID) {
editor->SaveFile(fileName);
if( ID == 0 ){
if( isInFile->GetState()) {
double aMin = angMin->GetNumber();
double aMax = angMax->GetNumber();
double aStep = angStep->GetNumber();
statusLabel->SetText("Creating DWBA.in.....");
InFileCreator("DWBA", "DWBA.in", aMin, aMax, aStep);
statusLabel->SetText("in-file created.");
}
bool isRunOK = true;
if( isRun->GetState() && IsFileExist("DWBA.in") ) {
//printf("run ptolemy...........\n");
statusLabel->SetText("Running Ptolemy.....");
int output = 1;
if( OS_Type == 1 ){
output = system("/opt/Ptolemy/Cleopatra/ptolemy <DWBA.in> DWBA.out");
}else{
output = system("/opt/Ptolemy/Cleopatra/ptolemy_mac <DWBA.in> DWBA.out");
}
statusLabel->SetText("Check terminal, if no massage, Ptolemy run well.");
printf("Ptolemy exist code : %d\n", output);
if( output == 0 ) {
isRunOK = true;
}else{
isRunOK = false;
statusLabel->SetText("Ptolemy exist with problems.");
}
}
if( isRunOK && isExtract->GetState() && IsFileExist("DWBA.out")){
int ElasticFlag = 0; // 1 for ratio to Rutherford, 2 for total Xsec, 3 for (n,n) Xsec
ElasticFlag = extractFlag->GetSelected();
statusLabel->SetText("Extracting X-sec.....");
ExtractXSec("DWBA.out", ElasticFlag);
statusLabel->SetText("X-sec Extracted.");
}
if( isRunOK && isPlot->GetState() && IsFileExist("DWBA.root")){
statusLabel->SetText("Plot X-sec.....");
PlotTGraphTObjArray("DWBA.root");
statusLabel->SetText("Plotted X-sec.");
}
}
/*
if( ID == 1 ){
string basicConfig = "reactionConfig.txt";
string heliosDetGeoFile = "detectorGeo.txt";
string excitationFile = "Ex.txt"; //when no file, only ground state
TString ptolemyRoot = ""; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root";
TString filename = "reaction.dat"; //when no file, no output
if( withDWBA->GetState() ) {
ptolemyRoot = "DWBA.root";
excitationFile = "DWBA.Ex.txt";
}
statusLabel->SetText("Running simulation.......");
Transfer( basicConfig, heliosDetGeoFile, excitationFile, ptolemyRoot, saveFileName, filename);
statusLabel->SetText("Plotting simulation.......");
Check_Simulation();
statusLabel->SetText("Plotted Simulation result");
}
if( ID == 2 ){
Check_Simulation();
statusLabel->SetText(" Run Simulation first.");
}
*/
if( ID == 3 ){
if( fileName != "" ){
statusLabel->SetText(fileName + " saved.");
}else{
statusLabel->SetText("cannot save HELP page.");
}
}
if( ID == 4 ){
fileName = "";
statusLabel->SetText("Help Page.");
/*
editor->LoadBuffer("==================== For Simulation");
editor->AddLine("");
editor->AddLine("1) Make sure you check :");
editor->AddLine(" a) reaction Config");
editor->AddLine(" b) detector Geo.");
editor->AddLine(" c) Ex List");
editor->AddLine("");
editor->AddLine("2) Not need to save file, fiel save when any button (except the Exit) is pressed.");
editor->AddLine("");
editor->AddLine("3) There is a checkbox for simulation with DWBA");
editor->AddLine(" This requires the existance of DWBA.root and DWBA.Ex.txt");
editor->AddLine(" These files can be generated by DWBA calculation.");
editor->AddLine(" Please change the angMin = 0 and angMax = 180.");
editor->AddLine("");
editor->AddLine("4) After simulation, it will plot the result.");
editor->AddLine(" To change the plotting, Click on the Config Simulation Plot.");
editor->AddLine("");
editor->AddLine("5) If the transfer.root is already here, simply Plot Simulation.");
editor->AddLine("");
*/
editor->LoadBuffer("==================== For DWBA");
///editor->AddLine("========================= For DWBA ");
editor->AddLine("");
editor->AddLine("1) Only need to change the DWBA setting.");
editor->AddLine("");
editor->AddLine("2) The GUI offer a view on the infile and outfile.");
editor->AddLine("");
editor->AddLine("3) For elastics scattering, there is a checkbox for plotting the ratio to RutherFord.");
editor->AddLine("");
editor->AddLine("4) The flow of the DWBA calculation is like this:");
editor->AddLine(" a) read the DWBA file and convert to DWBA.in");
editor->AddLine(" b) run Ptolemy from DWBA.in, and the output is DWBA.out");
editor->AddLine(" c) extract the cross section from the DWBA.out, and save :");
editor->AddLine(" * DWBA.Xsec.txt");
editor->AddLine(" * DWBA.Ex.txt");
editor->AddLine(" * DWBA.root");
editor->AddLine(" d) Plot the cross section from the DWBA.root.");
editor->AddLine("");
editor->AddLine("================ Tips for using the editor, both MAC or LINUX");
editor->AddLine("");
editor->AddLine("Ctrl+U | Delete current line. ");
editor->AddLine("Ctrl+C | Copy ");
editor->AddLine("Ctrl+V | Paste ");
editor->AddLine("=================================================== eof");
TString osTypeStr;
osTypeStr.Form("OS type is %s", (OS_Type == 0 ? "Mac" : "Linux"));
editor->AddLine(osTypeStr);
editor->AddLine(Form("Root version : %s",gROOT->GetVersion()));
int versionInt = gROOT->GetVersionInt();
if( versionInt < 62600 ) {
editor->AddLine("Please Update Root to v6.26/00");
}
}
/*
if( ID == 5) {
TH1F * temp = (TH1F*) gROOT->FindObjectAny("hExCal");
if( temp != NULL ){
fitAuto(temp, -1);
statusLabel->SetText("Auto Fit hExCal");
}else{
statusLabel->SetText("Cannot find historgram hExCal. Please Run Plot Simulation first.");
}
//gROOT->ProcessLine("fitAuto(hExCal, -1)");
}
*/
}
MyMainFrame::~MyMainFrame() {
// Clean up used widgets: frames, buttons, layout hints
fMain->Cleanup();
delete fMain;
}
void PtolemyGUI() {
new MyMainFrame(gClient->GetRoot(),800,600);
}