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

add Cleopatra from digios

Browse Source
This commit is contained in:
Ryan Tang 2023-04-03 16:03:48 -04:00
parent e1446834e3
commit eb2d06e67b
42 changed files with 19122 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
.vscode/settings.json vendored
View File

@ -75,7 +75,23 @@
"ChainMonitors.C": "cpp",
"GeneralSortAgent.C": "cpp",
"readRawTrace.C": "cpp",
"readTrace.C": "cpp"
"readTrace.C": "cpp",
"Cleopatra.C": "cpp",
"alpha.C": "cpp",
"DWBA_compare.C": "cpp",
"DWBARatio.C": "cpp",
"ExtractXSec.C": "cpp",
"ExtractXSecFromText.C": "cpp",
"FindThetaCM.C": "cpp",
"InFileCreator.C": "cpp",
"IsotopeShort.C": "cpp",
"knockout.C": "cpp",
"PlotSimulation.C": "cpp",
"PlotTGraphTObjArray.C": "cpp",
"transfer_test.C": "cpp",
"Transfer.C": "cpp",
"Simulation_Helper.C": "cpp",
"Check_Simulation.C": "cpp"
},
"better-comments.multilineComments": true,

646
Cleopatra/Check_Simulation.C Normal file
View File

@ -0,0 +1,646 @@
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TROOT.h>
#include <TObjArray.h>
#include <TStyle.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TF1.h>
#include <TArc.h>
#include <TMath.h>
#include <TLine.h>
#include <TSpectrum.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TLatex.h>
#include <TMacro.h>
#include <TObjArray.h>
#include <fstream>
#include <TCutG.h>
#include "../armory/AnalysisLib.h"
#include "../Cleopatra/Isotope.h"
double * FindRange(TString branch, TString gate, TTree * tree, double output[2]);
double ExtractNumber(int index, TMacro * macro);
TString ExtractString(int index, TMacro * macro);
vector<TString> StringToVector(TString str);
vector<int> intConvertor(vector<TString> arr);
vector<double> doubleConvertor(vector<TString> arr);
enum plotID { pEZ, /// 0
pRecoilXY, /// 1
pRecoilXY1, /// 2
pRecoilXY2, /// 3
pRecoilRZ, /// 4
pRecoilRTR, /// 5
pTDiffZ, /// 6
pThetaCM, /// 7
pThetaCM_Z, /// 8
pExCal, /// 9
pRecoilRThetaCM, /// 10
pArrayXY, /// 11
pInfo, /// 12
pHitID, /// 13
pElum1XY, /// 14
pEElum1R, /// 15
pElum1RThetaCM, /// 16
pEmpty }; /// 17
plotID StringToPlotID(TString str);
void Check_Simulation(TString filename = "transfer.root",
TString configFile = "../Armory/Check_Simulation_Config.txt",
Int_t padSize = 500,
bool outputCanvas = false){
printf("=========================== Check_Simulation.C\n");
TMacro * config = new TMacro(configFile);
int numLine = config->GetListOfLines()->GetSize();
int startLineNum = 0;
for( int i = 0; i < numLine ; i++){
TString haha = config->GetListOfLines()->At(i)->GetName();
haha.Remove(4);
if( haha != "////" ) {
startLineNum = i;
break;
}
}
TString gate = ExtractString(startLineNum+1, config);
double elumRange = ExtractNumber(startLineNum+2, config);
vector<double> thetaCMRange = doubleConvertor( StringToVector( ExtractString(startLineNum+3,config) ));
bool shownKELines = (ExtractString(startLineNum+4, config).Remove(4) == "true" ? true : false);
bool isOverRideEx = (ExtractString(startLineNum+5, config).Remove(4) == "true" ? true : false);
vector<double> oExRange = doubleConvertor( StringToVector ( ExtractString(startLineNum+6, config )));
printf("%s \n", gate.Data());
///==== config Canvas
vector<TString> plotConfig = StringToVector( ExtractString(startLineNum, config));
vector<plotID> canvas;
int colCount = 0;
int colCount_new = 0;
int rowCount = 1;
for( int i = 0; i < (int) plotConfig.size(); i++){
if( plotConfig[i] == "break" ) {
rowCount ++;
if( colCount_new > colCount ) colCount = colCount_new;
colCount_new = 0;
continue;
}
canvas.push_back( StringToPlotID(plotConfig[i]));
colCount_new ++;
}
if( colCount == 0 ) colCount = colCount_new;
///printf("plot row: %d, col: %d \n", rowCount, colCount);
vector<int> Div = {colCount, rowCount};
TFile * file = new TFile(filename, "read");
TTree * tree = (TTree*) file->Get("tree");
TObjArray * fxList = (TObjArray *) file->FindObjectAny("fxList");
TObjArray * txList = (TObjArray *) file->FindObjectAny("txList");
//================== reactionConfig
TMacro * reactionConfigTxt = (TMacro *) file->FindObjectAny("reactionConfig");
TString Reaction=reactionConfigTxt->GetName();
ReactionConfig reactionConfig = LoadReactionConfig(reactionConfigTxt);
int nEvent = reactionConfig.numEvents;
printf("number of events generated : %d \n", nEvent);
double xBeam = reactionConfig.beamX;
double yBeam = reactionConfig.beamY;
printf(" beam position : (%5.2f, %5.2f) mm \n", xBeam, yBeam);
gStyle->SetOptStat("");
gStyle->SetStatY(0.9);
gStyle->SetStatX(0.9);
gStyle->SetStatW(0.4);
gStyle->SetStatH(0.2);
gStyle->SetLabelSize(0.05, "XY");
gStyle->SetTitleFontSize(0.1);
double eRange[2] = {0, 10};
double zRange[3] = {400, -1000, 1000}; /// zRange[0] = nBin
double recoilERange[2];
vector<double> exList;
double ExRange[2];
//================================== detetcor Geometry
printf("=================================\n");
printf(" loading detector Geometry.\n");
TMacro * detGeoTxt = (TMacro *) file->FindObjectAny("detGeo");
DetGeo detGeo = LoadDetectorGeo(detGeoTxt);
double field = detGeo.Bfield;
TString fdmsg = field > 0 ? "out of plan" : "into plan";
TString msg2;
msg2.Form("field = %.2f T, %s", field, fdmsg.Data());
double prepDist = detGeo.detPerpDist;
double length = detGeo.detLength;
double posRecoil = detGeo.recoilPos;
double rhoRecoilIn = detGeo.recoilInnerRadius;
double rhoRecoilOut = detGeo.recoilOuterRadius;
double posRecoil1 = detGeo.recoilPos1;
double posRecoil2 = detGeo.recoilPos2;
vector<double> pos = detGeo.detPos;
float firstPos = detGeo.firstPos;
int rDet = detGeo.nDet;
int cDet = detGeo.mDet;
double elum1 = detGeo.elumPos1;
printf("number of row-Det : %d \n", rDet);
printf("number of col-Det : %d \n", cDet);
for(int i = 0; i < rDet ; i++){
if( firstPos > 0 ){
printf("%d, %10.2f mm - %10.2f mm \n", i, pos[i], pos[i] + length);
}else{
printf("%d, %10.2f mm - %10.2f mm \n", i, pos[i] - length , pos[i]);
}
}
printf("=================================\n");
int numDet = rDet * cDet;
zRange[1] = detGeo.zMin - 50;
zRange[2] = detGeo.zMax + 50;
printf(" zRange : %f - %f \n", zRange[1], zRange[2]);
printf("=================================\n");
//========================================= Ex List;
printf(" loading Ex list\n");
TMacro * exListMacro = (TMacro *) file->FindObjectAny("ExList");
int numEx = exListMacro->GetListOfLines()->GetSize() - 1 ;
for(int i = 1; i <= numEx ; i++){
string temp = exListMacro->GetListOfLines()->At(i)->GetName();
if( temp[0] == '#' ) break;
if( temp[0] == '/' ) continue;
vector<string> tempStr = SplitStr(temp, " ");
printf("%d | %s \n", i, tempStr[0].c_str());
exList.push_back(atof(tempStr[0].c_str()));
}
double exSpan = exList.back() - exList[0];
const int nExID = exList.size();
printf("========= number of excited states : %d \n", nExID);
ExRange[0] = exList[0] - exSpan * 0.2;
ExRange[1] = exList.back() + exSpan * 0.2;
if( isOverRideEx ) {
ExRange[0] = oExRange[0];
ExRange[1] = oExRange[1];
}
printf("=================================\n");
//========================================= reaction parameters
printf(" loading reaction parameters \n");
TMacro * reactionData = (TMacro *) file->FindObjectAny("reactionData");
double mass = ExtractNumber(0, reactionData);
double q = ExtractNumber(1, reactionData);
double beta = ExtractNumber(2, reactionData);
double Et = ExtractNumber(3, reactionData);
double massB = ExtractNumber(4, reactionData);
double alpha = ExtractNumber(5, reactionData);
double gamm = 1./TMath::Sqrt(1-beta*beta);
double slope = alpha * beta;
printf("\tmass-b : %f MeV/c2 \n", mass);
printf("\tcharge-b : %f \n", q);
printf("\tE-total : %f MeV \n", Et);
printf("\tmass-B : %f MeV/c2 \n", massB);
printf("\tbeta : %f \n", beta);
printf("\tslope : %f MeV/mm \n", slope);
printf("=================================\n");
//=================================== calculate Ranges
//eRange by zRange and exList
double QQ = (Et*Et + mass*mass - (massB-exList[0])*(massB-exList[0]))/2/Et;
double intercept = QQ/gamm - mass;
eRange[1] = intercept + zRange[2] * slope;
///printf("intercept of 0 MeV : %f MeV \n", intercept);
///printf("eRange 0 MeV : %f MeV \n", eRange[1]);
//thetaCMRange
///double momentum = sqrt(( Et*Et - pow(mass + massB - exList[0],2)) * ( Et*Et - pow(mass - massB + exList[0],2)))/2/Et;
///double thetaMax = acos( (beta * QQ- alpha / gamm * zRange[2])/momentum) * TMath::RadToDeg();
///thetaCMRange[1] = (int) TMath::Ceil(thetaMax/10.)*10;
///printf(" momentum : %f \n", momentum);
///printf(" thetaCM Max : %f \n", thetaMax);
///printf(" thetaCM Range : %d \n", thetaCMRange[1]);
//===================================================
printf("============================== Gate\n");
printf("gate : %s\n", gate.Data());
printf("====================================\n");
Int_t size[2] = {padSize,padSize}; ///x,y, single Canvas size
TCanvas * cCheck = new TCanvas("cCheck", "Check For Simulation", 0, 0, size[0]*Div[0], size[1]*Div[1]);
if(cCheck->GetShowEditor() )cCheck->ToggleEditor();
if(cCheck->GetShowToolBar() )cCheck->ToggleToolBar();
cCheck->Divide(Div[0],Div[1]);
for( int i = 1; i <= Div[0]*Div[1] ; i++){
cCheck->cd(i);
cCheck->cd(i)->SetGrid();
if( canvas[i-1] == pThetaCM ) {
cCheck->cd(i)->SetGrid(0,0);
cCheck->cd(i)->SetLogy();
}
if( canvas[i-1] == pHitID ){
cCheck->cd(i)->SetLogy();
}
plotID pID = canvas[i-1];
///########################################
if( pID == pEZ){
TH2F * hez = new TH2F("hez", Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", firstPos), zRange[0], zRange[1], zRange[2], 400, eRange[0], eRange[1]);
tree->Draw("e:z>>hez", gate, "colz");
if( shownKELines){
for( int i = 0; i < nExID ; i++){
fxList->At(i)->Draw("same");
}
}
}
if( pID == pRecoilXY ){
TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", posRecoil ), 400, -rhoRecoilOut, rhoRecoilOut, 400,-rhoRecoilOut, rhoRecoilOut);
tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz");
TArc * detArc1 = new TArc(0,0, rhoRecoilOut);
detArc1->SetLineColor(kBlue-8);
detArc1->SetFillStyle(0);
detArc1->Draw("same");
TArc * detArc2 = new TArc(0,0, rhoRecoilIn);
detArc2->SetLineColor(kBlue-8);
detArc2->SetFillStyle(0);
detArc2->Draw("same");
if( xBeam != 0. || yBeam != 0. ){
TArc * arc = new TArc(xBeam, yBeam, 1);
arc->SetLineColor(2);
detArc1->SetFillStyle(0);
arc->Draw("same");
}
}
if( pID == pRecoilXY1 ){
TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", posRecoil1 ), 400, -rhoRecoilOut, rhoRecoilOut, 400,-rhoRecoilOut, rhoRecoilOut);
tree->Draw("yRecoil1:xRecoil1>>hRecoilXY1", gate, "colz");
}
if( pID == pRecoilXY2 ){
TH2F * hRecoilXY2 = new TH2F("hRecoilXY2", Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", posRecoil2 ), 400, -rhoRecoilOut, rhoRecoilOut, 400,-rhoRecoilOut, rhoRecoilOut);
tree->Draw("yRecoil2:xRecoil2>>hRecoilXY2", gate, "colz");
}
if( pID == pRecoilRZ ){
TH2F * hRecoilRZ = new TH2F("hRecoilRZ", "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", zRange[0], zRange[1], zRange[2], 400,0, rhoRecoilOut);
tree->Draw("rhoRecoil:z>>hRecoilRZ", gate, "colz");
}
if( pID == pRecoilRTR ){
FindRange("TB", gate, tree, recoilERange);
TH2F * hRecoilRTR = new TH2F("hRecoilRTR", "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", 500, recoilERange[0], recoilERange[1], 500, 0, rhoRecoilOut);
tree->Draw("rhoRecoil:TB>>hRecoilRTR", gate, "colz");
}
if( pID == pTDiffZ ){
double tDiffRange [2];
FindRange("t-tB", gate, tree, tDiffRange);
TH2F * hTDiffZ = new TH2F("hTDiffZ", "time(Array) - time(Recoil) vs Z [gated]; z [mm]; time diff [ns]", zRange[0], zRange[1], zRange[2], 500, tDiffRange[0], tDiffRange[1]);
tree->Draw("t - tB : z >> hTDiffZ", gate, "colz");
}
if( pID == pThetaCM ){
TH1F * hThetaCM[nExID];
TLegend * legend = new TLegend(0.8,0.2,0.99,0.8);
double maxCount = 0;
int startID = 0; // set the start ExID
for( int i = startID; i < nExID; i++){
hThetaCM[i] = new TH1F(Form("hThetaCM%d", i), Form("thetaCM [gated] (ExID=%d); thetaCM [deg]; count", i), 200, thetaCMRange[0], thetaCMRange[1]);
hThetaCM[i]->SetLineColor(i+1-startID);
hThetaCM[i]->SetFillColor(i+1-startID);
hThetaCM[i]->SetFillStyle(3000+i-startID);
tree->Draw(Form("thetaCM>>hThetaCM%d", i), gate + Form("&& ExID==%d", i), "");
legend->AddEntry(hThetaCM[i], Form("Ex=%5.1f MeV", exList[i]));
double max = hThetaCM[i]->GetMaximum();
if( max > maxCount ) maxCount = max;
}
for( int i = startID; i < nExID; i++){
hThetaCM[i]->GetYaxis()->SetRangeUser(1, maxCount * 1.2);
if( i == startID ) {
hThetaCM[i]->Draw();
}else{
hThetaCM[i]->Draw("same");
}
}
legend->Draw();
}
if( pID == pThetaCM_Z ){
TH2F *hThetaCM_Z = new TH2F("hThetaCM_Z","ThetaCM vs Z ; Z [mm]; thetaCM [deg]",zRange[0], zRange[1], zRange[2], 200, thetaCMRange[0], thetaCMRange[1]);
tree->Draw("thetaCM:z>>hThetaCM_Z",gate,"col");
if( shownKELines){
for( int i = 0; i < nExID ; i++){
txList->At(i)->Draw("same");
}
}
}
if( pID == pExCal ){
TH1F * hExCal = new TH1F("hExCal", Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV", (ExRange[1]-ExRange[0])/400.*1000), 400, ExRange[0], ExRange[1]);
tree->Draw("ExCal>>hExCal", gate, "");
Isotope hRecoil(reactionConfig.recoilHeavyA, reactionConfig.recoilHeavyZ);
double Sn = hRecoil.CalSp(0,1);
double Sp = hRecoil.CalSp(1,0);
double Sa = hRecoil.CalSp2(4,2);
double S2n = hRecoil.CalSp(0, 2);
printf("Heavy recoil: %s \n", hRecoil.Name.c_str());
printf("Sn : %f MeV/u \n", Sn);
printf("Sp : %f MeV/u \n", Sp);
printf("Sa : %f MeV/u \n", Sa);
printf("S2n : %f MeV/u \n", S2n);
double yMax = hExCal->GetMaximum();
TLine * lineSn = new TLine(Sn, 0, Sn, yMax); lineSn->SetLineColor(2); lineSn->Draw("");
TLine * lineSp = new TLine(Sp, 0, Sp, yMax); lineSp->SetLineColor(4); lineSp->Draw("same");
TLine * lineSa = new TLine(Sa, 0, Sa, yMax); lineSa->SetLineColor(6); lineSa->Draw("same");
TLine * lineS2n = new TLine(S2n, 0, S2n, yMax); lineS2n->SetLineColor(8); lineS2n->Draw("same");
TLatex * text = new TLatex();
text->SetTextFont(82);
text->SetTextSize(0.06);
text->SetTextColor(2); text->DrawLatex(Sn, yMax*0.9, "S_{n}");
text->SetTextColor(4); text->DrawLatex(Sp, yMax*0.9, "S_{p}");
text->SetTextColor(6); text->DrawLatex(Sa, yMax*0.9, "S_{a}");
text->SetTextColor(8); text->DrawLatex(S2n, yMax*0.9, "S_{2n}");
}
if( pID == pRecoilRThetaCM ){
TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, rhoRecoilOut);
tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz");
}
if( pID == pArrayXY ){
TH2F * hArrayXY = new TH2F("hArrayXY", "Array-XY [gated]; X [mm]; Y [mm]", 400, -prepDist*1.5, prepDist*1.5, 400, -prepDist*1.5, prepDist*1.5);
tree->Draw("yArray:xArray>>hArrayXY", gate, "colz");
}
if( pID == pInfo ){
TLatex text;
text.SetNDC();
text.SetTextFont(82);
text.SetTextSize(0.06);
text.SetTextColor(2);
text.DrawLatex(0., 0.9, Reaction);
text.DrawLatex(0., 0.8, msg2);
text.SetTextColor(1);
text.DrawLatex(0., 0.7, "gate:");
text.SetTextColor(2);
//check gate text length, if > 30, break by "&&"
int ll = gate.Length();
if( ll > 30 ) {
vector<string> strList = SplitStr( (string) gate.Data(), "&&");
for( int i = 0; i < strList.size(); i++){
text.DrawLatex(0., 0.6 - 0.05*i, (TString) strList[i]);
}
}else{
text.DrawLatex(0., 0.6, gate);
}
if( xBeam != 0.0 || yBeam != 0.0 ){
text.DrawLatex(0.0, 0.1, Form("Bema pos: (%4.1f, %4.1f) mm", xBeam, yBeam));
}
}
if( pID == pElum1XY ){
TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", elum1), 400, -elumRange, elumRange, 400, -elumRange, elumRange);
tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz");
double count = hElum1XY->GetEntries();
if( count < 2000. ) {
hElum1XY->SetMarkerStyle(7);
if( count < 500. ) hElum1XY->SetMarkerStyle(3);
hElum1XY->Draw("scat");
}
}
if( pID == pEElum1R ){
TH2F * hEElum1Rho = new TH2F("hEElum1Rho", "Elum-1 E-R [gated]; R[mm]; Energy[MeV]", 400, 0, elumRange, 400, eRange[0], eRange[1]);
tree->Draw("Tb:rhoElum1>>hEElum1Rho", gate, "colz");
}
if( pID == pElum1RThetaCM){
int angBin = 400;
TH2F * hElum1RThetaCM = new TH2F("hElum1RThetaCM", "Elum-1 rho vs ThetaCM [gated]; thatCM [deg]; Elum- rho [mm]", angBin, thetaCMRange[0], thetaCMRange[1], 400, 0, elumRange);
tree->Draw("rhoElum1:thetaCM>>hElum1RThetaCM", gate, "colz");
TH1F * htemp = (TH1F *) hElum1RThetaCM->ProjectionX("htemp");
double rel = (thetaCMRange[1] - thetaCMRange[0])*1.0/angBin;
printf("angular resolution : %f deg \n", rel);
vector<double> xList;
double old_y = 0;
for( int i = 1; i <= angBin; i++){
double y = htemp->GetBinContent(i);
if( old_y == 0 && y > 0) xList.push_back(htemp->GetBinCenter(i));
if( old_y > 0 && y == 0 ) xList.push_back(htemp->GetBinCenter(i));
old_y = y;
}
printf("list of gaps :\n");
for( int i = 0; i < (int) xList.size(); i+=2){
printf("%d | %.3f - %.3f deg\n", i, xList[i], xList[i+1]);
}
TF1 f1("f1", "sin(x)");
double acceptance = 0;
double err1 = 0;
double err2 = 0;
for( int i = 0; i < (int) xList.size(); i += 2 ){
acceptance += f1.Integral(xList[i] * TMath::DegToRad(), xList[i+1] * TMath::DegToRad() ) * TMath::TwoPi();
err1 += f1.Integral((xList[i]-rel) * TMath::DegToRad(), (xList[i+1] + rel) * TMath::DegToRad() ) * TMath::TwoPi();
err2 += f1.Integral((xList[i]+rel) * TMath::DegToRad(), (xList[i+1] - rel) * TMath::DegToRad() ) * TMath::TwoPi();
}
printf("acceptance = %f sr +- %f \n", acceptance, (err1-err2)/2);
TLatex text;
text.SetTextFont(82);
text.SetTextSize(0.06);
text.SetTextColor(2);
text.SetTextAngle(90);
for( int i = 0; i < (int) xList.size(); i++){
text.DrawLatex(xList[i], elumRange/2, Form("%.2f", xList[i]));
}
text.SetNDC();
text.SetTextAngle(0);
text.DrawLatex(0.15, 0.15, Form("accp. = %.2f(%.2f) msr", acceptance * 1000., (err1-err2)*1000./2));
}
if( pID == pHitID ){
printf("=======================meaning of Hit ID\n");
printf(" 1 = light recoil hit array & heavy recoil hit recoil\n");
printf(" 0 = no detector\n");
printf(" -1 = light recoil go opposite side of array\n");
printf(" -2 = light recoil hit > det width\n");
printf(" -3 = light recoil hit > array \n");
printf(" -4 = light recoil hit blocker \n");
printf(" -10 = light recoil orbit radius too big \n");
printf(" -11 = light recoil orbit radius too small\n");
printf(" -12 = when reocol at the same side of array, light recoil blocked by recoil detector\n");
printf(" -13 = more than 3 loops\n");
printf(" -14 = heavy recoil did not hit recoil \n");
printf(" -15 = cannot find hit on array\n");
printf(" -20 = unknown\n");
printf("===========================================\n");
TH1F * hHit = new TH1F("hHit", "hit; hit-ID; count", 13, -11, 2);
tree->Draw("hit>>hHit", "", "");
}
///#######################################################
}
cCheck->Modified();
cCheck->Update();
if( outputCanvas ){
TDatime dateTime;
TString outPNGName = Form("Sim_%d%02d%02d_%06d.png", dateTime.GetYear(), dateTime.GetMonth(), dateTime.GetDay(), dateTime.GetTime());
cCheck->SaveAs(outPNGName);
printf("%s\n", outPNGName.Data());
gROOT->ProcessLine(".q");
}
}
///============================================================
///============================================================
double * FindRange(TString branch, TString gate, TTree * tree, double output[2]){
tree->Draw(Form("%s>>temp1", branch.Data()), gate);
TH1F * temp1 = (TH1F *) gROOT->FindObjectAny("temp1");
output[1] = temp1->GetXaxis()->GetXmax();
output[0] = temp1->GetXaxis()->GetXmin();
delete temp1;
return output;
}
double ExtractNumber(int index, TMacro * macro){
TString field = macro->GetListOfLines()->At(index)->GetName();
int pos = field.First('/');
if( pos >= 0 ) field.Remove(pos);
return field.Atof();
}
TString ExtractString(int index, TMacro * macro){
TString field = macro->GetListOfLines()->At(index)->GetName();
int pos = field.First('/');
if( pos >= 0 ) field.Remove(pos);
return field;
}
vector<TString> StringToVector(TString str){
vector<TString> temp;
bool startFlag = false;
bool endFlag = false;
string jaja="";
for(int i = 0; i < str.Length(); i++){
if( str[i] == '{' ) {
startFlag = true;
continue;
}
if( str[i] == ' '){
continue;
}
if( startFlag && !endFlag){
if( str[i] == ',' ){
temp.push_back(jaja);
jaja="";
continue;
}
if( str[i] == '}') {
temp.push_back(jaja);
endFlag = true;
continue;
}
jaja += str[i];
}
}
return temp;
}
vector<int> intConvertor(vector<TString> arr){
vector<int> out ;
for( int i = 0 ; i < (int) arr.size(); i++){
out.push_back( arr[i].Atoi());
}
return out;
}
vector<double> doubleConvertor(vector<TString> arr){
vector<double> out ;
for( int i = 0 ; i < (int) arr.size(); i++){
out.push_back( arr[i].Atof());
}
return out;
}
plotID StringToPlotID(TString str){
if( str == "pEZ") return plotID::pEZ; ///0
if( str == "pRecoilXY") return plotID::pRecoilXY; /// 1
if( str == "pRecoilXY1" ) return plotID::pRecoilXY1; /// 2
if( str == "pRecoilXY2" ) return plotID::pRecoilXY2; /// 3
if( str == "pRecoilRZ" ) return plotID::pRecoilRZ; /// 4
if( str == "pRecoilRTR" ) return plotID::pRecoilRTR; /// 5
if( str == "pTDiffZ" ) return plotID::pTDiffZ; /// 6
if( str == "pThetaCM" ) return plotID::pThetaCM; /// 7
if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z; /// 8
if( str == "pExCal" ) return plotID::pExCal; /// 9
if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM; /// 10
if( str == "pArrayXY" ) return plotID::pArrayXY; /// 11
if( str == "pInfo" ) return plotID::pInfo; /// 12
if( str == "pHitID" ) return plotID::pHitID; /// 13
if( str == "pElum1XY" ) return plotID::pElum1XY; /// 14
if( str == "pEElum1R" ) return plotID::pEElum1R; /// 14
if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM; /// 15
if( str == "pEmpty" ) return plotID::pEmpty ; /// 16
return plotID::pEmpty;
}

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);
}

552
Cleopatra/ExtractXSec.h Normal file
View File

@ -0,0 +1,552 @@
/***********************************************************************
*
* 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 <stdlib.h>
#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/AnalysisLib.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);
}

58
Cleopatra/FindThetaCM.C Normal file
View File

@ -0,0 +1,58 @@
/***********************************************************************
*
* This is FindThetaCM.C, To calculate the thetaCM convrage for each detector
*
* This required two inout files: basicReactionConfig.txt
* detectorGeo.txt
*
*-------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h> /* atof */
#include <cmath>
#include <vector>
#include "FindThetaCM.h"
int main(int argc, char *argv[]){
printf("=================================================================\n");
printf("=== Find ThetaCM convrage for each detector at Ex ====\n");
printf("=================================================================\n");
if(argc < 2 || argc > 6) {
printf("Usage: ./FindThetaCM Ex\n");
printf("Usage: ./FindThetaCM Ex nDiv\n");
printf("Usage: ./FindThetaCM Ex nDiv X-Ratio\n");
printf("Usage: ./FindThetaCM Ex nDiv X-Ratio reactionTxt detGeoTxt\n");
exit(0);
}
double Ex = 0;
double xRatio = 0.95;
int nDiv = 1;
string reactionTxt = "reactionConfig.txt";
string detGeoTxt = "detectorGeo.txt";
if ( argc >= 2 ){
Ex = atof(argv[1]);
}
if ( argc >= 3 ){
nDiv = atoi(argv[2]);
}
if ( argc >= 4 ){
xRatio = atof(argv[3]);
}
if ( argc >= 5 ){
reactionTxt = argv[4];
}
if ( argc >= 6 ){
detGeoTxt = argv[5];
}
FindThetaCM(Ex, nDiv, xRatio, reactionTxt, detGeoTxt);
return 0;
}

214
Cleopatra/FindThetaCM.h Normal file
View File

@ -0,0 +1,214 @@
/***********************************************************************
*
* This is FindThetaCM.h, To calculate the thetaCM convrage for each detector
*
* This required two inout files: basicReactionConfig.txt
* detectorGeo.txt
*
*-------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Nov-18, 2018
* email: goluckyryan@gmail.com
* ********************************************************************/
#include "TFile.h"
#include "TTree.h"
#include "TMacro.h"
#include "TObjArray.h"
#include "TGraph.h"
#include "../Cleopatra/HELIOS_LIB.h"
void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95,
string basicConfig="reactionConfig.txt",
string detGeoFileName = "detectorGeo.txt"){
//---- reaction
int AA, zA; //beam
int Aa, za; //target
int Ab, zb; //recoil-1
double ExA;
//---- beam
double KEAmean, KEAsigma; // MeV/u , assume Guassian
double thetaMean, thetaSigma; // mrad , assume Guassian due to small angle
double xBeam, yBeam; // mm
/**///========================================================= load files
ReactionConfig reactionConfig;
DetGeo detGeo;
TMacro * haha = new TMacro();
if( haha->ReadFile(basicConfig.c_str()) > 0 ){
reactionConfig = LoadReactionConfig(haha);
PrintReactionConfig(reactionConfig);
KEAmean = reactionConfig.beamEnergy;
KEAsigma = reactionConfig.beamEnergySigma;
thetaMean = reactionConfig.beamAngle;
thetaSigma = reactionConfig.beamAngleSigma;
xBeam = reactionConfig.beamX;
yBeam = reactionConfig.beamY;
AA = reactionConfig.beamA; zA = reactionConfig.beamZ;
Aa = reactionConfig.targetA; za = reactionConfig.targetZ;
Ab = reactionConfig.recoilLightA; zb = reactionConfig.recoilLightZ;
ExA = reactionConfig.beamEx[0];
}else{
printf("cannot load %s \n", basicConfig.c_str());
return;
}
vector<double> pos;
double a = 11.5;
double length = 50.5;
double firstPos = 0;
int iDet = 6;
int jDet = 4;
double BField = 0;
//=============================================================
//=============================================================
//=============================================================
//===== Set Reaction
TransferReaction reaction;
int AB = AA+Aa-Ab, zB = zA+za-zb;
reaction.SetA(AA,zA);
reaction.Seta(Aa,za);
reaction.Setb(Ab,zb);
reaction.SetB(AB,zB);
reaction.SetIncidentEnergyAngle(KEAmean, 0, 0);
reaction.SetExB(Ex);
reaction.SetExA(ExA);
reaction.CalReactionConstant();
printf("===================================================\n");
printf("=========== %27s ===========\n", reaction.GetReactionName().Data());
printf("===================================================\n");
printf("----- loading reaction from : %s. \n", basicConfig.c_str());
printf(" Ex A: %7.3f MeV\n", ExA);
printf(" KE: %7.4f \n", KEAmean);
printf(" theta: %7.4f \n", thetaMean);
printf("offset(x,y): %7.4f, %7.4f mm \n", xBeam, yBeam);
printf(" Q-value: %7.4f MeV \n", reaction.GetQValue() );
printf(" Max Ex: %7.4f MeV \n", reaction.GetMaxExB() );
printf("===================================================\n");
printf("----- loading detector geometery : %s.", detGeoFileName.c_str());
TMacro * kaka = new TMacro();
if( kaka->ReadFile(detGeoFileName.c_str()) > 0 ){
detGeo = LoadDetectorGeo(kaka);
pos = detGeo.detPos;
a = detGeo.detPerpDist;
length = detGeo.detLength;
firstPos = detGeo.firstPos;
iDet = detGeo.nDet;
jDet = detGeo.mDet;
BField = detGeo.Bfield;
printf("... done.\n");
PrintDetGeo(detGeo);
}else{
printf("... fail\n");
return;
}
vector<double> midPos;
for(int i = 0; i < iDet; i++){
if( firstPos > 0 ){
midPos.push_back(pos[i]+length/2.);
}else{
midPos.push_back(pos[i]-length/2.);
}
}
//calculate a TGraph for thetaCM vs z
double px[100];
double py[100];
double mb = reaction.GetMass_b();
double kCM = reaction.GetMomentumbCM();
double q = TMath::Sqrt(mb*mb + kCM * kCM );
double beta = reaction.GetReactionBeta() ;
double slope = 299.792458 * zb * abs(BField) / TMath::TwoPi() * beta / 1000.; // MeV/mm
double gamma = reaction.GetReactionGamma();
for(int i = 0; i < 100; i++){
double thetacm = (i + 5.) * TMath::DegToRad();
double temp = TMath::TwoPi() * slope / beta / kCM * a / TMath::Sin(thetacm);
px[i] = beta /slope * (gamma * beta * q - gamma * kCM * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi());
py[i] = thetacm * TMath::RadToDeg();
}
//find minimum z position
TGraph * xt = new TGraph(100, py, px);
xt->SetName("xt");
///double zMin0 = xt->Eval(0);
///printf("z for thetaCM = 0 : %f mm \n", zMin0);
///xt->Draw("AC*");
/// find the minimum z position and the corresponding theta
double zMin0 = 0;
double tMin0 = 0;
for( double ttt = 3 ; ttt < 20 ; ttt += 0.1 ){
double zzz = xt->Eval(ttt);
if( zzz < zMin0 ) {
zMin0 = zzz;
tMin0 = ttt;
}
}
printf(" z min %f mm at thetaCM %f deg \n", zMin0, tMin0);
TGraph * tx = new TGraph(100, px, py);
tx->SetName(Form("tx"));
tx->SetLineColor(4);
//tx->Draw("AC*");
/**///========================================================= result
printf("==== ThetaCM in degree =================\n");
printf("========================= x-ratio : %f, number of division : %d \n", XRATION, nDivision);
printf("\n");
for( int j = 0; j < nDivision + 1; j++) printf("%5.2f ", -XRATION + 2*XRATION/nDivision*j);
printf(" <<-- in X \n");
for( int j = 0; j < nDivision + 1; j++) printf("%5s ", " | ");
printf("\n");
for( int j = 0; j < nDivision + 1; j++) printf("%5.2f ", length/2 -length*XRATION/2 + length*XRATION/nDivision*j);
printf(" <<-- in cm \n\n");
printf("========================= Ex : %6.4f MeV\n", Ex);
printf(" %6s - %6s | %6s, %6s, %6s\n", "Min", "Max", "Mean", "Dt", "sin(x)dx * 180/pi");
printf("-------------------------------------------------\n");
for( int i = 0; i < iDet; i++){
double zMin = midPos[i]-length*XRATION/2.;
double zMax = midPos[i]+length*XRATION/2.;
double zLength = zMax - zMin;
double zStep = zLength/(nDivision);
for( int j = 0 ; j < nDivision ; j++){
double tMin = (zMin + j*zStep > zMin0) ? tx->Eval(zMin + j*zStep) : TMath::QuietNaN();
double tMax = (zMin + (j+1)*zStep > zMin0) ? tx->Eval(zMin + (j+1)*zStep) : TMath::QuietNaN();
double tMean = (tMax + tMin)/2.;
double dt = (tMax - tMin);
double sintdt = TMath::Sin(tMean * TMath::DegToRad()) * dt ;
printf(" det-%d[%d]: %6.2f - %6.2f | %6.2f, %6.2f, %6.4f\n", i, j, tMin, tMax, tMean, dt, sintdt);
}
if( nDivision > 0 ) printf("--------------\n");
}
printf("================================================= \n");
}

1791
Cleopatra/HELIOS_LIB.h Normal file
View File

File diff suppressed because it is too large Load Diff

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);
}

404
Cleopatra/InFileCreator.h Normal file
View File

@ -0,0 +1,404 @@
/***********************************************************************
*
* 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 <stdlib.h> /* atof */
#include <vector>
#include "../Cleopatra/Isotope.h" // for geting Z
#include "potentials.h"
#include "../armory/AnalysisLib.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];
int lenStr20 = str2[0].length();
size_t posTemp1 = str2[0].find(",");
string mass_a = str2[0].substr(1, posTemp1-1);
size_t posTemp2 = str2[0].find(")");
string mass_b = str2[0].substr(posTemp1+1, posTemp2-posTemp1-1);
///printf(" mass_a : |%s| , mass_b : |%s| \n", mass_a.c_str(), mass_b.c_str());
Isotope iso_a(mass_a);
Isotope iso_b(mass_b);
// Check is the Reaction supported
bool isReactionSupported = false;
bool isTransferReaction = true;
if( iso_a.A <= 4 && iso_a.Z <= 2 && iso_b.A <=4 && iso_b.Z <=2 ) isReactionSupported = true;
///======= elastics-ish scattering
if( iso_a.Mass == iso_b.Mass ) isTransferReaction = false;
///======= p/n-exchange is not supported
if( iso_a.A == iso_b.A && iso_a.Z != iso_b.Z ) isReactionSupported = false;
///======= 3-nucleons transfer is not supported. e.g. (n,a), (p,a), (a,n), (a,p)
int numNucleonsTransfer = iso_a.A - iso_b.A;
if( abs(numNucleonsTransfer) >= 3 ) isReactionSupported = false;
if( isReactionSupported == false ){
printf(" ===> Ignored. Reaction type not supported. \n");
continue;
}
// Continues to decode the input string
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 iso_A(str1[0]);
Isotope iso_B(str2[1]);
/// check is iso_A or iso_B exist in the mass table
if( iso_A.Mass == -404 || iso_B.Mass == -404 ){
printf(" ===> Error! mass does not found. \n");
continue;
}
/// check reaction valid by balancing the A and Z number;
if( iso_A.A + iso_a.A != iso_B.A + iso_b.A || iso_A.Z + iso_a.Z != iso_B.Z + iso_b.Z ) {
printf("====> ERROR! A-number or Z-number not balanced. \n");
Isotope isotopeK(iso_A.A + iso_a.A - iso_b.A, iso_A.Z + iso_a.Z - iso_b.Z);
printf(" try : %s(%s,%s)%s ??\n", iso_A.Name.c_str(), iso_a.Name.c_str(), iso_b.Name.c_str(), isotopeK.Name.c_str());
continue;
}
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(iso_a.A - iso_b.A));
node = orbital.substr(0,1);
// single nucleon transfer
if( abs(iso_a.A - iso_b.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(iso_a.A - iso_b.A) == 2 ){
size_t posEq = orbital.find('=');
lValue = orbital.substr(posEq+1,1);
spdf=atoi(lValue.c_str());
}
if( abs(iso_a.A - iso_b.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);
int factor = 1;
if( unit == "MeV/u") factor = iso_a.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 = iso_a.Mass + iso_A.Mass - iso_b.Mass - iso_B.Mass;
///printf("Q-Value = %f MeV\n", Qvalue);
//##########################################################
//============ write ptolmey infile
numOfReaction ++ ;
//================ elastic-ish transfer
if( isTransferReaction == false ){
if ( atof(Ex.c_str()) == 0.0 ) {
fprintf(file_out, "$============================================ ELab=%5.2f(%s+%s)%s\n",
totalBeamEnergy, mass_a.c_str(), isoA.c_str(), potential.c_str());
fprintf(file_out, "reset\n");
fprintf(file_out, "CHANNEL %s + %s\n", mass_a.c_str(), isoA.c_str());
fprintf(file_out, "r0target\n");
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, iso_A.A, iso_A.Z, totalBeamEnergy, iso_a.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(), mass_a.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, iso_A.A, iso_A.Z, totalBeamEnergy, iso_a.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, iso_A.A, iso_A.Z, totalBeamEnergy - atof(Ex.c_str()), iso_a.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");
}
}
//================ Transfer reaction
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);
//-------- Projectile (the light particle)
if( abs(numNucleonsTransfer) == 1 ){
if( iso_a.A <= 2 && iso_a.Z <= 1 && iso_b.A <=2 && iso_b.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 <= iso_a.A && iso_a.A <= 4) || (3 <= iso_b.A && iso_b.A <= 4) ){
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( iso_a.Z + iso_b.Z == 2){ // (t,d) or (d,t)
fprintf(file_out, "nodes=0 l=0 jp=1/2 spfacp=1.30 v=172.88 r=0.56 a=0.69 param1=0.64 param2=1.15 rc=2.0\n");
}
if( iso_a.Z + iso_b.Z == 3){ // (3He,d) or (d, 3He)
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( iso_b.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");
}
}
}else if( abs(numNucleonsTransfer) == 2 ){ // 2 nucleons 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, "wavefunction phiffer\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 = iso_B.CalSp(0,1);
double pThreshold = iso_B.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(iso_a.A-iso_b.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(iso_a.A-iso_b.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, iso_A.A, iso_A.Z, totalBeamEnergy, iso_a.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, iso_B.A, iso_B.Z, eBeam, iso_b.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");
}
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();
}

535
Cleopatra/Isotope.h Normal file
View File

@ -0,0 +1,535 @@
/***********************************************************************
*
* 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="/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(){
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

54
Cleopatra/IsotopeShort.C Normal file
View File

@ -0,0 +1,54 @@
/***********************************************************************
*
* This is IsotopeShort.C for isotope mass-related quatilies
* for python web
*
* -----------------------------------------------------
* To compile
*
* g++ IsotopeShort.C -o IsotopesShort
*
* ------------------------------------------------------
* 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 << "./NuclearData Sym" << endl;
cout << "./NuclearData A Z" << endl;
exit(0);
}
int main (int argc, char *argv[]) {
if ( argc != 2 && argc != 3 && argc != 6) Usage();
Isotope iso1, iso2;
int Z, A, Za, Aa;
if (argc == 2){
iso1.SetIsoByName(argv[1]);
}
if (argc == 3){
A= atoi(argv[1]);
Z= atoi(argv[2]);
iso1.SetIso(A, Z);
}
//iso1.Print();
printf("A:%3d\n", iso1.A);
printf("Z:%3d\n", iso1.Z);
printf("N:%3d\n", iso1.A-iso1.Z);
printf("Name:%s|\n",iso1.Name.c_str());
printf("Mass:%13.4f\n", iso1.Mass);
printf("Sn:%13.4f\n", iso1.CalSp(0,1));
printf("Sp:%13.4f\n", iso1.CalSp(1,0));
printf("Sa:%13.4f\n", iso1.CalSp2(4,2));
printf("S2n:%13.4f\n", iso1.CalSp(0,2));
}

31
Cleopatra/PlotSimulation.C Normal file
View File

@ -0,0 +1,31 @@
#include <fstream>
#include <stdlib.h>
#include "../Armory/Check_Simulation.C"
using namespace std;
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("=================== Plot Simulation Canvas ======================\n");
printf("=================================================================\n");
if(argc < 2 ) {
printf("Usage: ./PlotSimulation input_root_file [config]\n");
exit(0);
}else{
printf("ROOT file : %s \n", argv[1]);
}
string rootFile = argv[1];
string config = "../Armory/Check_Simulation_Config.txt";
if( argc >= 3 ) config = argv[2];
printf("Config File : %s \n", config.c_str());
Int_t padSize = 500;
Check_Simulation(rootFile, config, padSize, true);
}

58
Cleopatra/PlotTGraphTObjArray.C Normal file
View File

@ -0,0 +1,58 @@
/***********************************************************************
*
* 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 [savePNG]\n");
printf(" savePNG : 1 or 0, default is 0\n");
exit(0);
}else{
printf("From file : %s \n", argv[1]);
printf("=========== Press Ctrl+C to end.\n");
}
string readFile = argv[1];
bool isSavePNG = false;
if( argc >= 3) isSavePNG = atoi(argv[2]);
if( isSavePNG ){
PlotTGraphTObjArray(readFile, true);
}else{
TApplication app ("app", &argc, argv);
PlotTGraphTObjArray(readFile);
app.Run(); //anything after this line is not running
}
return 0;
}

98
Cleopatra/PlotTGraphTObjArray.h Normal file
View File

@ -0,0 +1,98 @@
/***********************************************************************
*
* 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, bool isSavePNG = false){
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.6, 0.9, 0.9); //x1, y1, x2, y2
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();
if( isSavePNG ){
TDatime dateTime;
TString outPNGName = Form("Xsec_%d%02d%02d_%06d.png", dateTime.GetYear(), dateTime.GetMonth(), dateTime.GetDay(), dateTime.GetTime());
cPlots->SaveAs(outPNGName);
printf("%s\n", outPNGName.Data());
gROOT->ProcessLine(".q");
}
}

591
Cleopatra/Simulation_Helper.C Normal file
View File

@ -0,0 +1,591 @@
#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 "../Cleopatra/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
TString isoFileName;
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;
TGTextEntry * txtName ;
TGTextEntry * txtEx ;
public:
MyMainFrame(const TGWindow *p,UInt_t w,UInt_t h);
virtual ~MyMainFrame();
void Command(int);
void OpenFile(int);
void GetData();
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 | kLHintsExpandX | kLHintsExpandY, 2,2,2,2));
TGVerticalFrame *hframe1 = new TGVerticalFrame(fMain,600,600 );
hframe->AddFrame(hframe1);
TGVerticalFrame *hframe2 = new TGVerticalFrame(fMain,600,800 );
hframe->AddFrame(hframe2,new TGLayoutHints( kLHintsExpandX | kLHintsExpandY, 2,2,2,2));
fileName = "reactionConfig.txt";
TGHorizontalFrame *hframe00 = new TGHorizontalFrame(hframe2,600,600 );
hframe2->AddFrame(hframe00, new TGLayoutHints(kLHintsCenterX | kLHintsExpandX , 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, new TGLayoutHints(kLHintsCenterX | kLHintsExpandX | kLHintsExpandY, 2,2,2,2));
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 | kLHintsExpandX, 2,2,2,2));
{//================= Simulation group
TGGroupFrame * simFrame = new TGGroupFrame(hframe1, "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(hframe1, "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));
}
{//====================== Nuclear data API
TGGroupFrame * dataFrame = new TGGroupFrame(hframe1, "Nuclear Data", kVerticalFrame);
hframe1->AddFrame(dataFrame, new TGLayoutHints(kLHintsCenterX, 5,5,3,4));
TGHorizontalFrame * hfData = new TGHorizontalFrame(dataFrame); dataFrame->AddFrame(hfData, new TGLayoutHints(kLHintsNormal, 0, 0, 5, 0));
TGVerticalFrame * vfLabel = new TGVerticalFrame(hfData, 200); hfData->AddFrame(vfLabel );
TGVerticalFrame * vfTxt = new TGVerticalFrame(hfData); hfData->AddFrame(vfTxt);
TGLayoutHints * haha = new TGLayoutHints(kLHintsRight | kLHintsCenterY, 5,5,5,2);
TGLayoutHints * kaka = new TGLayoutHints(kLHintsLeft | kLHintsCenterY, 5,5,0,0);
TGLabel * lb1 = new TGLabel(vfLabel, "Nuclear Name :"); vfLabel->AddFrame(lb1, haha);
TGLabel * lb2 = new TGLabel(vfLabel, "Max Ex [MeV] :"); vfLabel->AddFrame(lb2, haha);
txtName = new TGTextEntry(vfTxt, "25F"); vfTxt->AddFrame(txtName, kaka); txtName->Resize(50, 20);
txtEx = new TGTextEntry(vfTxt, "0"); vfTxt->AddFrame(txtEx, kaka); txtEx->Resize(50, 20);
TGTextButton *GetData = new TGTextButton(dataFrame, "Get Data");
GetData->SetWidth(150);
GetData->SetHeight(40);
GetData->ChangeOptions( GetData->GetOptions() | kFixedSize );
GetData->Connect("Clicked()","MyMainFrame",this,"GetData()");
dataFrame->AddFrame(GetData,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 | kLHintsExpandX, 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 = "../working/detectorGeo.txt";
if ( ID == 1 ) fileName = "../working/reactionConfig.txt";
if ( ID == 2 ) fileName = "../working/Ex.txt";
if ( ID == 3 ) fileName = "../working/DWBA";
if ( ID == 4 ) fileName = "../working/Check_Simulation_Config.txt";
if ( ID == 5 ) fileName = "../working/DWBA.in";
if ( ID == 6 ) fileName = "../working/DWBA.out";
if ( ID == 7 ) fileName = "../working/DWBA.Xsec.txt";
if ( ID == 8 ) fileName = isoFileName;
//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::GetData(){
TString name = txtName->GetText();
TString maxEx = txtEx->GetText();
TString cmd = "../Cleopatra/nuclear_data.py " + name + " " + maxEx;
system(cmd.Data());
statusLabel->SetText("Check termial.");
//isoFileName = name + ".txt";
//OpenFile(8);
}
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("../Cleopatra/ptolemy <DWBA.in> DWBA.out");
}else{
output = system("../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->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 Simulation_Helper() {
new MyMainFrame(gClient->GetRoot(),800,600);
}

75
Cleopatra/Transfer.C Normal file
View File

@ -0,0 +1,75 @@
/***********************************************************************
*
* This is Transfer.C for simulation of transfer reaction.
*
* -----------------------------------------------------
* This program will call the root library and compile in g++
* compilation:
* g++ Transfer.C -o Transfer `root-config --cflags --glibs`
*
* ------------------------------------------------------
* created by Ryan (Tsz Leung) Tang, Feb-4, 2019
* email: goluckyryan@gmail.com
* ********************************************************************/
#include <fstream>
#include <stdlib.h>
#include "Transfer.h"
using namespace std;
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("========== Simulate Transfer reaction in HELIOS ==========\n");
printf("=================================================================\n");
if(argc == 2 || argc > 8) {
printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6] [7]\n");
printf(" default file name \n");
printf(" [1] reactionConfig.txt (input) reaction Setting \n");
printf(" [2] detectorGeo.txt (input) detector Setting \n");
printf(" [3] Ex.txt (input) Excitation energies \n");
printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n");
printf(" [5] transfer.root (output) rootFile name for output \n");
printf(" [6] reaction.dat (output) Key reaction parameters \n");
printf(" [7] plot (input) will it plot stuffs [1/0] \n");
printf("------------------------------------------------------\n");
return 0 ;
}
string basicConfig = "reactionConfig.txt";
string heliosDetGeoFile = "detectorGeo.txt";
string excitationFile = "Ex.txt"; //when no file, only ground state
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root";
TString filename = "reaction.dat"; //when no file, no output
bool isPlot = false;
if( argc >= 2) basicConfig = argv[1];
if( argc >= 3) heliosDetGeoFile = argv[2];
if( argc >= 4) excitationFile = argv[3];
if( argc >= 5) ptolemyRoot = argv[4];
if( argc >= 6) saveFileName = argv[5];
if( argc >= 7) filename = argv[6];
if( argc >= 8) isPlot = atoi(argv[7]);
Transfer( basicConfig, heliosDetGeoFile, excitationFile, ptolemyRoot, saveFileName, filename);
//run Armory/Check_Simulation
if( isPlot ){
ifstream file_in;
file_in.open("../Armory/Check_Simulation.C", ios::in);
if( file_in){
printf("---- running ../Armory/Check_Simulation.C on %s \n", saveFileName.Data());
TString cmd;
cmd.Form("root -l '../Armory/Check_Simulation.C(\"%s\", 500)'", saveFileName.Data());
system(cmd.Data());
}else{
printf("cannot find ../Armory/Check_Simulation.C \n");
}
}
}

722
Cleopatra/Transfer.h Normal file
View File

@ -0,0 +1,722 @@
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TF1.h"
#include "TTree.h"
#include "TRandom.h"
#include "TGraph.h"
#include "TMacro.h"
#include <stdlib.h>
#include <vector>
#include <fstream>
#include <TObjArray.h>
#include "HELIOS_LIB.h"
double exDistFunc(Double_t *x, Double_t * par){
return par[(int) x[0]];
}
void Transfer(
string basicConfig = "reactionConfig.txt",
string heliosDetGeoFile = "detectorGeo.txt",
string excitationFile = "Ex.txt", ///when no file, only ground state
TString ptolemyRoot = "DWBA.root", /// when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root",
TString filename = "reaction.dat"){ /// when no file, no output.
//############################################# Set Reaction
TransferReaction reaction;
reaction.SetReactionFromFile(basicConfig);
printf("*****************************************************************\n");
printf("*\e[1m\e[33m %27s \e[0m*\n", reaction.GetReactionName().Data());
printf("*****************************************************************\n");
printf("----- loading reaction setting from %s. \n", basicConfig.c_str());
printf("\e[32m#################################### Beam \e[0m\n");
const ReactionConfig reactionConfig = reaction.GetRectionConfig();
PrintReactionConfig(reactionConfig);
vector<float> ExAList = reactionConfig.beamEx;
int nExA = (int) ExAList.size();
//############################################# Set HELIOS
printf("\e[32m#################################### HELIOS configuration\e[0m\n");
HELIOS helios;
helios.SetDetectorGeometry(heliosDetGeoFile);
const DetGeo detGeo = helios.GetDetectorGeometry();
printf("==================================== E-Z plot slope\n");
double betaRect = reaction.GetReactionBeta() ;
double gamma = reaction.GetReactionGamma();
double mb = reaction.GetMass_b();
double pCM = reaction.GetMomentumbCM();
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
double slope = 299.792458 * reaction.GetCharge_b() * abs(helios.GetBField()) / TMath::TwoPi() * betaRect / 1000.; /// MeV/mm
printf(" e-z slope : %f MeV/mm\n", slope);
double intercept = q/gamma - mb; // MeV
printf(" e-z intercept (ground state) : %f MeV\n", intercept);
//############################################# save reaction.dat
if( filename != "" ) {
FILE * keyParaOut;
keyParaOut = fopen (filename.Data(), "w+");
printf("=========== save key reaction constants to %s \n", filename.Data());
fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetMass_b(), "mass_b");
fprintf(keyParaOut, "%-15d //%s\n", reaction.GetCharge_b(), "charge_b");
fprintf(keyParaOut, "%-15.8f //%s\n", reaction.GetReactionBeta(), "betaCM");
fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetCMTotalEnergy(), "Ecm");
fprintf(keyParaOut, "%-15.4f //%s\n", reaction.GetMass_B(), "mass_B");
fprintf(keyParaOut, "%-15.4f //%s\n", slope/betaRect, "alpha=slope/betaRect");
fflush(keyParaOut);
fclose(keyParaOut);
}
//############################################# Target scattering, only energy loss
bool isTargetScattering = reactionConfig.isTargetScattering;
float density = reactionConfig.targetDensity;
float targetThickness = reactionConfig.targetThickness;
if(isTargetScattering) printf("\e[32m#################################### Target Scattering\e[0m\n");
TargetScattering msA;
TargetScattering msB;
TargetScattering msb;
if(reactionConfig.isTargetScattering) printf("======== Target : (thickness : %6.2f um) x (density : %6.2f g/cm3) = %6.2f ug/cm2\n",
targetThickness * 1e+4,
density,
targetThickness * density * 1e+6);
if( reactionConfig.isTargetScattering ){
msA.LoadStoppingPower(reactionConfig.beamStoppingPowerFile);
msb.LoadStoppingPower(reactionConfig.recoilLightStoppingPowerFile);
msB.LoadStoppingPower(reactionConfig.recoilHeavyStoppingPowerFile);
}
//############################################# Decay of particle-B
Decay decay;
if(reactionConfig.isDecay) {
printf("\e[32m#################################### Decay\e[0m\n");
decay.SetMotherDaugther(reactionConfig.recoilHeavyA,
reactionConfig.recoilHeavyZ,
reactionConfig.heavyDecayA,
reactionConfig.heavyDecayZ);
}
//############################################# loading excitation energy
printf("\e[32m#################################### excitation energies\e[0m\n");
vector<double> ExKnown;
vector<double> ExStrength;
vector<double> ExWidth;
vector<double> SF;
vector<double> y0; /// intercept of e-z plot
vector<double> kCM; /// momentum of b in CM frame
printf("----- loading excitation energy levels (%s).", excitationFile.c_str());
ifstream file;
file.open(excitationFile.c_str());
string isotopeName;
if( file.is_open() ){
string line;
while( getline(file, line) ){
///printf("%s \n", line.c_str());
if( line.substr(0,2) == "//" ) continue;
if( line.substr(0,2) == "#=" ) break;
vector<string> str = SplitStr(line, " ");
ExKnown.push_back(atof(str[0].c_str()));
ExStrength.push_back(atof(str[1].c_str()));
SF.push_back(atof(str[2].c_str()));
ExWidth.push_back(atof(str[3].c_str()));
}
file.close();
printf("... done.\n");
int n = (int) ExKnown.size();
printf("%3s | %7s | %5s | %3s | %10s | %5s \n", "", "Ex[MeV]", "Xsec", "SF", "sigma[MeV]", "y0[MeV]");
printf("----+---------+------+-----+------------+--------\n");
for(int i = 0; i < n ; i++){
reaction.SetExB(ExKnown[i]);
reaction.CalReactionConstant();
kCM.push_back(reaction.GetMomentumbCM());
y0.push_back(TMath::Sqrt(mb*mb + kCM[i]*kCM[i])/gamma - mb);
if( reactionConfig.isDecay ) {
TLorentzVector temp(0,0,0,0);
int decayID = decay.CalDecay(temp, ExKnown[i], 0);
if( decayID == 1) {
printf("%3d | %7.2f | %5.2f | %3.1f | %5.3f | %5.2f --> Decay. \n", i, ExKnown[i], ExStrength[i], SF[i], ExWidth[i], y0[i]);
}else{
printf("%3d | %7.2f | %5.2f | %3.1f | %5.3f | %5.2f \n", i, ExKnown[i], ExStrength[i], SF[i], ExWidth[i], y0[i]);
}
}else{
printf("%3d | %7.2f | %5.2f | %3.1f | %5.3f | %5.2f \n", i, ExKnown[i], ExStrength[i], SF[i], ExWidth[i], y0[i]);
}
}
printf("----+---------+-------+-----+------------+--------\n");
}else{
printf("... fail ------> only ground state.\n");
ExKnown.push_back(0.0);
ExStrength.push_back(1.0);
ExWidth.push_back(0.0);
reaction.SetExB(ExKnown[0]);
reaction.CalReactionConstant();
kCM.push_back(reaction.GetMomentumbCM());
y0.push_back(TMath::Sqrt(mb*mb + kCM[0]*kCM[0])/gamma - mb);
}
//---- create Ex-distribution
TF1 * exDist = NULL;
if( ExKnown.size() > 1 ) {
printf("---- creating Ex-distribution \n");
int exSize = ExKnown.size();
exDist = new TF1("exDist", exDistFunc, 0, exSize, exSize);
for(int i = 0; i < exSize; i++){
exDist->SetParameter(i, ExStrength[i]*SF[i]);
}
}
//############################################# Load DWBAroot for thetaCM distribution
printf("\e[32m#################################### Load DWBA input : %s \e[0m\n", ptolemyRoot.Data());
TF1 * dist = NULL;
TFile * distFile = new TFile(ptolemyRoot, "read");
TObjArray * distList = NULL;
if( distFile->IsOpen() ) {
distList = (TObjArray *) distFile->FindObjectAny("pList"); // the function List
int distSize = distList->GetLast() + 1;
if( distSize != ExKnown.size() ) {
printf(" The number of distribution from Ptolmey Calculation is not equal to number of Ex input \n");
printf(" --> the Ptolmey calculation is probably not matched with Ex input.\n");
printf(" .... not use DWBA input. \n");
distFile->Close();
}
}else{
printf("------- no DWBA input. \n");
}
//############################################# build tree
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
TFile * saveFile = new TFile(saveFileName, "recreate");
TTree * tree = new TTree("tree", "tree");
TMacro config(basicConfig.c_str());
TMacro detGeoTxt(heliosDetGeoFile.c_str());
TMacro exList(excitationFile.c_str());
TMacro reactionData(filename.Data());
double KEAmean = reactionConfig.beamEnergy;
TString str;
str.Form("%s @ %.2f MeV/u", reaction.GetReactionName_Latex().Data(), KEAmean);
config.SetName(str.Data());
config.Write("reactionConfig");
detGeoTxt.Write("detGeo");
exList.Write("ExList");
reactionData.Write("reactionData");
if( distList != NULL ) distList->Write("DWBA", 1);
TMacro hitMeaning;
str = "=======================meaning of Hit ID\n"; hitMeaning.AddLine(str.Data());
str = " 1 = light recoil hit array & heavy recoil hit recoil\n"; hitMeaning.AddLine(str.Data());
str = " 0 = no detector\n"; hitMeaning.AddLine(str.Data());
str = " -1 = light recoil go opposite side of array\n"; hitMeaning.AddLine(str.Data());
str = " -2 = light recoil hit > det width\n"; hitMeaning.AddLine(str.Data());
str = " -3 = light recoil hit > array \n"; hitMeaning.AddLine(str.Data());
str = " -4 = light recoil hit blocker \n"; hitMeaning.AddLine(str.Data());
str = " -10 = light recoil orbit radius too big \n"; hitMeaning.AddLine(str.Data());
str = " -11 = light recoil orbit radius too small\n"; hitMeaning.AddLine(str.Data());
str = " -12 = when reocol at the same side of array, light recoil blocked by recoil detector\n"; hitMeaning.AddLine(str.Data());
str = " -13 = more than 3 loops\n"; hitMeaning.AddLine(str.Data());
str = " -14 = heavy recoil did not hit recoil \n"; hitMeaning.AddLine(str.Data());
str = " -15 = cannot find hit on array\n"; hitMeaning.AddLine(str.Data());
str = " -20 = unknown\n"; hitMeaning.AddLine(str.Data());
str = "===========================================\n"; hitMeaning.AddLine(str.Data());
hitMeaning.Write("hitMeaning");
int hit; /// the output of Helios.CalHit
tree->Branch("hit", &hit, "hit/I");
double thetab, phib, Tb;
double thetaB, phiB, TB;
tree->Branch("thetab", &thetab, "thetab/D");
tree->Branch("phib", &phib, "phib/D");
tree->Branch("Tb", &Tb, "Tb/D");
tree->Branch("thetaB", &thetaB, "thetaB/D");
tree->Branch("phiB", &phiB, "phiB/D");
tree->Branch("TB", &TB, "TB/D");
double thetaCM;
tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
double e, z, detX, t, z0, tB;
tree->Branch("e", &e, "energy_light/D");
tree->Branch("x", &detX, "detector_x/D");
tree->Branch("z", &z, "array_hit_z/D");
tree->Branch("z0", &z0, "z-cycle/D");
tree->Branch("t", &t, "cycle_time_light/D");
tree->Branch("tB", &tB, "recoil_hit_time/D"); /// hit time for recoil on the recoil detector
int loop, detID, detRowID;
tree->Branch("detID", &detID, "detID/I");
tree->Branch("detRowID", &detRowID, "detRowID/I");
tree->Branch("loop", &loop, "loop/I");
double rho, rhoB; ///orbit radius
tree->Branch("rho", &rho, "orbit_radius_light/D");
tree->Branch("rhoB", &rhoB, "orbit_radius_heavy/D");
int ExAID;
double ExA;
tree->Branch("ExAID", &ExAID, "ExAID/I");
tree->Branch("ExA", &ExA, "ExA/D");
int ExID;
double Ex;
tree->Branch("ExID", &ExID, "ExID/I");
tree->Branch("Ex", &Ex, "Ex/D");
double ExCal, thetaCMCal;
tree->Branch("ExCal", &ExCal, "ExCal/D");
tree->Branch("thetaCMCal", &thetaCMCal, "thetaCMCal/D");
double KEA, theta, phi;
tree->Branch("beamTheta", &theta, "beamTheta/D");
tree->Branch("beamPhi", &phi, "beamPhi/D");
tree->Branch("beamKEA", &KEA, "beamKEA/D");
double TbLoss; /// energy loss of particle-b from target scattering
double KEAnew; ///beam energy after target scattering
double depth; /// reaction depth;
double Ecm;
if( reactionConfig.isTargetScattering ){
tree->Branch("depth", &depth, "depth/D");
tree->Branch("TbLoss", &TbLoss, "TbLoss/D");
tree->Branch("KEAnew", &KEAnew, "KEAnew/D");
tree->Branch("Ecm", &Ecm, "Ecm/D");
}
double decayTheta; /// the change of thetaB due to decay
double xRecoil_d, yRecoil_d, rhoRecoil_d, Td;
if( reactionConfig.isDecay ) {
tree->Branch("decayTheta", &decayTheta, "decayTheta/D");
tree->Branch("xRecoil_d", &xRecoil_d, "xRecoil_d/D");
tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D");
tree->Branch("rhoRecoil_d", &rhoRecoil_d, "rhoRecoil_d/D");
tree->Branch("Td", &Td, "Td/D");
}
double xArray, yArray, rhoArray; ///x, y, rho positon of particle-b on PSD
tree->Branch("xArray", &xArray, "xArray/D");
tree->Branch("yArray", &yArray, "yArray/D");
tree->Branch("rhoArray", &rhoArray, "rhoArray/D");
double xRecoil, yRecoil, rhoRecoil; /// x, y, rho position of particle-B on recoil-detector
tree->Branch("xRecoil", &xRecoil, "xRecoil/D");
tree->Branch("yRecoil", &yRecoil, "yRecoil/D");
tree->Branch("rhoRecoil", &rhoRecoil, "rhoRecoil/D");
///in case need ELUM
double xElum1, yElum1, rhoElum1;
if( detGeo.elumPos1 != 0 ) {
tree->Branch("xElum1", &xElum1, "xElum1/D");
tree->Branch("yElum1", &yElum1, "yElum1/D");
tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D");
}
double xElum2, yElum2, rhoElum2;
if( detGeo.elumPos2 != 0 ) {
tree->Branch("xElum2", &xElum2, "xElum2/D");
tree->Branch("yElum2", &yElum2, "yElum2/D");
tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D");
}
///in case need other recoil detector.
double xRecoil1, yRecoil1, rhoRecoil1;
if( detGeo.recoilPos1 != 0 ){
tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D");
tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D");
tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D");
}
double xRecoil2, yRecoil2, rhoRecoil2;
if( detGeo.recoilPos2 != 0 ){
tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D");
tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D");
tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D");
}
//======= function for e-z plot for ideal case
printf("++++ generate functions\n");
TObjArray * gList = new TObjArray();
gList->SetName("Constant thetaCM lines");
const int gxSize = 50;
TF1 ** gx = new TF1*[gxSize];
TString name;
for( int i = 0; i < gxSize; i++){
name.Form("g%d", i);
gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000);
double thetacm = i * TMath::DegToRad();
double gS2 = TMath::Power(TMath::Sin(thetacm)*gamma,2);
gx[i]->SetParameter(0, TMath::Cos(thetacm));
gx[i]->SetParameter(1, mb*mb*(1-gS2));
gx[i]->SetParameter(2, TMath::Power(slope/betaRect,2));
gx[i]->SetParameter(3, 1-gS2);
gx[i]->SetParameter(4, mb);
gx[i]->SetParameter(5, -gS2*slope);
gx[i]->SetNpx(1000);
gList->Add(gx[i]);
printf("/");
if( i > 1 && i % 40 == 0 ) printf("\n");
}
gList->Write("gList", TObject::kSingleKey);
printf(" %d constant thetaCM functions\n", gxSize);
int n = ExKnown.size();
TObjArray * fList = new TObjArray();
TF1** f = new TF1*[n];
for( int i = 0; i< n ; i++){
name.Form("f%d", i);
f[i] = new TF1(name, "[0] + [1] * x", -1000, 1000);
f[i]->SetParameter(0, y0[i]);
f[i]->SetParameter(1, slope);
f[i]->SetNpx(1000);
fList->Add(f[i]);
printf(".");
}
fList->Write("fList", TObject::kSingleKey);
printf(" %d e-z infinte-small detector functions\n", n);
//--- cal modified f
TObjArray * fxList = new TObjArray();
TGraph ** fx = new TGraph*[n];
vector<double> px, py;
int countfx = 0;
for( int j = 0 ; j < n; j++){
double a = helios.GetDetRadius();
double q = TMath::Sqrt(mb*mb + kCM[j] * kCM[j] );
px.clear();
py.clear();
countfx = 0;
for(int i = 0; i < 100; i++){
double thetacm = TMath::Pi()/TMath::Log(100) * (TMath::Log(100) - TMath::Log(100-i)) ;//using log scale, for more point in small angle.
double temp = TMath::TwoPi() * slope / betaRect / kCM[j] * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
double pyTemp = gamma * q - mb - gamma * betaRect * kCM[j] * TMath::Cos(thetacm);
if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
px.push_back(pxTemp);
py.push_back(pyTemp);
countfx ++;
}
fx[j] = new TGraph(countfx, &px[0], &py[0]);
name.Form("fx%d", j);
fx[j]->SetName(name);
fx[j]->SetLineColor(4);
fxList->Add(fx[j]);
printf(",");
}
fxList->Write("fxList", TObject::kSingleKey);
printf(" %d e-z finite-size detector functions\n", n);
//--- cal modified thetaCM vs z
TObjArray * txList = new TObjArray();
TGraph ** tx = new TGraph*[n];
for( int j = 0 ; j < n; j++){
double a = helios.GetDetRadius();
double q = TMath::Sqrt(mb*mb + kCM[j] * kCM[j] );
px.clear();
py.clear();
countfx = 0;
for(int i = 0; i < 100; i++){
double thetacm = (i + 8.) * TMath::DegToRad();
double temp = TMath::TwoPi() * slope / betaRect / kCM[j] * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi());
double pyTemp = thetacm * TMath::RadToDeg();
if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
px.push_back(pxTemp);
py.push_back(pyTemp);
countfx ++;
}
tx[j] = new TGraph(countfx, &px[0], &py[0]);
name.Form("tx%d", j);
tx[j]->SetName(name);
tx[j]->SetLineColor(4);
txList->Add(tx[j]);
printf("*");
}
txList->Write("txList", TObject::kSingleKey);
printf(" %d thetaCM-z for finite-size detector functions\n", n);
//========timer
TBenchmark clock;
bool shown ;
clock.Reset();
clock.Start("timer");
shown = false;
//change the number of event into human easy-to-read form
int numEvent = reactionConfig.numEvents;
int digitLen = TMath::Floor(TMath::Log10(numEvent));
TString numEventStr;
if( 3 <= digitLen && digitLen < 6 ){
numEventStr.Form("%5.1f kilo", numEvent/1000.);
}else if ( 6<= digitLen && digitLen < 9 ){
numEventStr.Form("%6.2f million", numEvent/1e6);
}else if ( 9<= digitLen ){
numEventStr.Form("%6.2f billion", numEvent/1e9);
}
printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data());
//====================================================== calculate event
int count = 0;
for( int i = 0; i < numEvent; i++){
bool redoFlag = true;
if( !reactionConfig.isRedo ) redoFlag = false;
do{
//==== Set Ex of A
ExAID = gRandom->Integer(nExA);
ExA = ExAList[ExAID];
reaction.SetExA(ExA);
//==== Set Ex of B
if( ExKnown.size() == 1 ) {
ExID = 0;
Ex = ExKnown[0] + (ExWidth[0] == 0 ? 0 : gRandom->Gaus(0, ExWidth[0]));
}else{
ExID = exDist->GetRandom();
Ex = ExKnown[ExID]+ (ExWidth[ExID] == 0 ? 0 : gRandom->Gaus(0, ExWidth[ExID]));
}
reaction.SetExB(Ex);
//==== Set incident beam
KEA = reactionConfig.beamEnergy;
if( reactionConfig.beamEnergySigma == 0 ){
KEA = reactionConfig.beamEnergy;
}else{
KEA = gRandom->Gaus(reactionConfig.beamEnergy, reactionConfig.beamEnergySigma);
}
theta = 0.0;
if( reactionConfig.beamAngleSigma == 0 ){
theta = reactionConfig.beamAngle;
}else{
theta = gRandom->Gaus(reactionConfig.beamAngle, reactionConfig.beamAngleSigma);
}
phi = 0.0;
//==== for taregt scattering
reaction.SetIncidentEnergyAngle(KEA, theta, 0.);
reaction.CalReactionConstant();
TLorentzVector PA = reaction.GetPA();
//depth = 0;
if( isTargetScattering ){
//==== Target scattering, only energy loss
depth = targetThickness * gRandom->Rndm();
msA.SetTarget(density, depth);
TLorentzVector PAnew = msA.Scattering(PA);
KEAnew = msA.GetKE()/reactionConfig.beamA;
reaction.SetIncidentEnergyAngle(KEAnew, theta, phi);
reaction.CalReactionConstant();
Ecm = reaction.GetCMTotalKE();
}
//==== Calculate thetaCM, phiCM
if( distFile->IsOpen()){
dist = (TF1 *) distList->At(ExID);
thetaCM = dist->GetRandom() / 180. * TMath::Pi();
}else{
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
}
double phiCM = TMath::TwoPi() * gRandom->Rndm();
//==== Calculate reaction
TLorentzVector * output = reaction.Event(thetaCM, phiCM);
TLorentzVector Pb = output[2];
TLorentzVector PB = output[3];
//==== Calculate energy loss of scattered and recoil in target
if( isTargetScattering ){
if( Pb.Theta() < TMath::PiOver2() ){
msb.SetTarget(density, targetThickness - depth);
}else{
msb.SetTarget(density, depth);
}
Pb = msb.Scattering(Pb);
TbLoss = msb.GetKELoss();
msB.SetTarget(density, targetThickness - depth);
PB = msB.Scattering(PB);
}else{
TbLoss = 0;
}
//======= Decay of particle-B
int decayID = 0;
int new_zB = reactionConfig.recoilHeavyZ;
if( reactionConfig.isDecay){
//decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2.); // decay to ground state
decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state
if( decayID == 1 ){
PB = decay.GetDaugther_D();
//decayTheta = decay.GetAngleChange();
decayTheta = decay.GetThetaCM();
new_zB = reactionConfig.heavyDecayZ;
}else{
decayTheta = TMath::QuietNaN();
}
}
//################################### tree branches
//===== reaction
thetab = Pb.Theta() * TMath::RadToDeg();
thetaB = PB.Theta() * TMath::RadToDeg();
Tb = Pb.E() - Pb.M();
TB = PB.E() - PB.M();
phib = Pb.Phi() * TMath::RadToDeg();
phiB = PB.Phi() * TMath::RadToDeg();
//==== Helios
///printf(" thetaCM : %f \n", thetaCM * TMath::RadToDeg());
if( Tb > 0 || TB > 0 ){
helios.CalArrayHit(Pb, reaction.GetCharge_b());
helios.CalRecoilHit(PB, new_zB);
hit = 2;
while( hit > 1 ){ hit = helios.DetAcceptance(); } /// while hit > 1, goto next loop;
trajectory orb_b = helios.GetTrajectory_b();
trajectory orb_B = helios.GetTrajectory_B();
e = helios.GetEnergy() + gRandom->Gaus(0, detGeo.eSigma);
double ranX = gRandom->Gaus(0, detGeo.zSigma);
z = orb_b.z + ranX;
detX = helios.GetDetX() + ranX;
z0 = orb_b.z0;
t = orb_b.t;
loop = orb_b.loop;
detID = orb_b.detID;
detRowID = orb_b.detRowID;
rho = orb_b.rho;
rhoArray = orb_b.R;
xArray = orb_b.x;
yArray = orb_b.y;
//ELUM
if( detGeo.elumPos1 != 0 ){
xElum1 = helios.GetXPos(detGeo.elumPos1);
yElum1 = helios.GetYPos(detGeo.elumPos1);
rhoElum1 = helios.GetR(detGeo.elumPos1);
}
if( detGeo.elumPos2 != 0 ){
xElum2 = helios.GetXPos(detGeo.elumPos2);
yElum2 = helios.GetYPos(detGeo.elumPos2);
rhoElum2 = helios.GetR(detGeo.elumPos2);
}
//Recoil
rhoRecoil = orb_B.R;
tB = orb_B.t;
xRecoil = orb_B.x;
yRecoil = orb_B.y;
rhoB = orb_B.rho;
//other recoil detectors
if ( detGeo.recoilPos1 != 0 ){
xRecoil1 = helios.GetRecoilXPos(detGeo.recoilPos1);
yRecoil1 = helios.GetRecoilYPos(detGeo.recoilPos1);
rhoRecoil1 = helios.GetRecoilR(detGeo.recoilPos1);
}
if ( detGeo.recoilPos2 != 0 ){
xRecoil2 = helios.GetRecoilXPos(detGeo.recoilPos2);
yRecoil2 = helios.GetRecoilYPos(detGeo.recoilPos2);
rhoRecoil2 = helios.GetRecoilR(detGeo.recoilPos2);
}
reaction.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius());
ExCal = reaction.GetEx();
thetaCMCal = reaction.GetThetaCM();
//change thetaCM into deg
thetaCM = thetaCM * TMath::RadToDeg();
//if decay, get the light decay particle on the recoil;
if( reactionConfig.isDecay ){
if( decayID == 1 ){
TLorentzVector Pd = decay.GetDaugther_d();
Td = Pd.E() - Pd.M();
helios.CalRecoilHit(Pd, reactionConfig.heavyDecayZ);
trajectory orb_d = helios.GetTrajectory_B();
rhoRecoil_d = orb_d.R;
xRecoil_d = orb_d.x;
yRecoil_d = orb_d.y;
}else{
rhoRecoil_d = TMath::QuietNaN();
xRecoil_d = TMath::QuietNaN();
yRecoil_d = TMath::QuietNaN();
}
}
}else{
hit = -404;
}
if( hit == 1) count ++;
if( reactionConfig.isRedo ){
if( hit == 1) {
redoFlag = false;
}else{
redoFlag = true;
//printf("%d, %2d, thetaCM : %f, theta : %f, z0: %f \n", i, hit, thetaCM * TMath::RadToDeg(), thetab, helios.GetZ0());
}
}else{
redoFlag = false;
}
}while( redoFlag );
tree->Fill();
//#################################################################### Timer
clock.Stop("timer");
Double_t time = clock.GetRealTime("timer");
clock.Start("timer");
if ( !shown ) {
if (fmod(time, 10) < 1 ){
printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
shown = 1;
}
}else{
if (fmod(time, 10) > 9 ){
shown = 0;
}
}
}
saveFile->Write();
saveFile->Close();
distFile->Close();
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
//gROOT->ProcessLine(".q");
}

172
Cleopatra/alpha.C Normal file
View File

@ -0,0 +1,172 @@
#include "HELIOS_LIB.h"
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TF1.h"
#include "TTree.h"
#include "TRandom.h"
#include <stdlib.h>
#include <vector>
#include <fstream>
#include <TObjArray.h>
//----------- usage
// $root transfer.C+ | tee output.txt
// this will same the massage to output.txt
const double ma = 3727.3792; // alpha mass
void alpha(){
//================================================= User Setting
const int numEnergy = 4;
double energy [numEnergy] = {3.18, 5.16, 5.49, 5.81};
int numEvent = 1000000;
//---- HELIOS detector geometry
//string heliosDetGeoFile = "detectorGeo.txt";
string heliosDetGeoFile = "";
double BField = 2.5; // T
double BFieldTheta = 0.; // direction of B-field
bool isCoincidentWithRecoil = false;
double eSigma = 0.040 ; // detector energy sigma MeV
double zSigma = 0.500 ; // detector position sigma mm
//---- save root file name
TString saveFileName = "alpha.root";
//=============================================================
//=============================================================
printf("===================================================\n");
printf("============= Alpha source in HELIOS ============\n");
printf("===================================================\n");
printf("========= Alpha Enegry : \n");
for( int i = 0; i < numEnergy ; i++){
printf("%2d | %6.2f MeV\n", i, energy[i]);
}
//======== Set HELIOS
printf("############################################## HELIOS configuration\n");
HELIOS helios;
helios.OverrideMagneticFieldDirection(BFieldTheta);
helios.OverrideFirstPos(-700);
//helios.OverrideDetectorDistance(5);
bool sethelios = helios.SetDetectorGeometry(heliosDetGeoFile);
if( !sethelios){
helios.OverrideMagneticField(BField);
printf("======== B-field : %5.2f T, Theta : %6.2f deg\n", BField, BFieldTheta);
}
helios.SetCoincidentWithRecoil(isCoincidentWithRecoil);
printf("========== energy resol.: %f MeV\n", eSigma);
printf("=========== pos-Z resol.: %f mm \n", zSigma);
//====================== build tree
TFile * saveFile = new TFile(saveFileName, "recreate");
TTree * tree = new TTree("tree", "tree");
double theta, phi, T;
int hit; // the output of Helios.CalHit
double e, z, x, t;
int loop, detID;
double dphi, rho; //rad of rotation, and radius
int energyID;
double xHit, yHit;
tree->Branch("hit", &hit, "hit/I");
tree->Branch("theta", &theta, "theta/D");
tree->Branch("phi", &phi, "phi/D");
tree->Branch("T", &T, "T/D");
tree->Branch("energy", &energy, "energy/D");
tree->Branch("energyID", &energyID, "energyID/I");
tree->Branch("e", &e, "e/D");
tree->Branch("x", &x, "x/D");
tree->Branch("z", &z, "z/D");
tree->Branch("t", &t, "t/D");
tree->Branch("detID", &detID, "detID/I");
tree->Branch("loop", &loop, "loop/I");
tree->Branch("dphi", &dphi, "dphi/D");
tree->Branch("rho", &rho, "rho/D");
tree->Branch("xHit", &xHit, "xHit/D");
tree->Branch("yHit", &yHit, "yHit/D");
//========timer
TBenchmark clock;
bool shown ;
clock.Reset();
clock.Start("timer");
shown = false;
printf("############################################## generating %d events \n", numEvent);
//====================================================== calculate
int count = 0;
TLorentzVector P;
TVector3 v;
for( int i = 0; i < numEvent; i++){
//==== generate alpha
theta = TMath::ACos(2 * gRandom->Rndm() - 1) ;
phi = TMath::TwoPi() * gRandom->Rndm();
energyID = gRandom->Integer(numEnergy);
T = energy[energyID];
double p = TMath::Sqrt( ( ma + T )*(ma + T) - ma* ma);
v.SetMagThetaPhi(p, theta, phi);
P.SetVectM(v, ma);
//################################### tree branches
//==== Helios
hit = helios.CalHit(P, 2, P, 2);
e = helios.GetEnergy() + gRandom->Gaus(0, eSigma);
z = helios.GetZ() ;
x = helios.GetX() + gRandom->Gaus(0, zSigma);
t = helios.GetTime();
loop = helios.GetLoop();
detID = helios.GetDetID();
dphi = helios.GetdPhi();
rho = helios.GetRho();
xHit = helios.GetXPos(z);
yHit = helios.GetYPos(z);
z += gRandom->Gaus(0, zSigma);
if( hit == 1) {
count ++;
}
tree->Fill();
//#################################################################### Timer
clock.Stop("timer");
Double_t time = clock.GetRealTime("timer");
clock.Start("timer");
if ( !shown ) {
if (fmod(time, 10) < 1 ){
printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
shown = 1;
}
}else{
if (fmod(time, 10) > 9 ){
shown = 0;
}
}
}
saveFile->Write();
saveFile->Close();
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
gROOT->ProcessLine(".q");
}

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

37
Cleopatra/install_ptolmey_mac.txt Normal file
View File

@ -0,0 +1,37 @@
Installation instructions
- Install Docker Desktop for Mac.
https://download.docker.com/mac/stable/Docker.dmg
Start the app. If it is running, the docker icon (a whale stacked with
container boxes) should be present on the top menu bar on the right hand
side.
- Open a terminal window and type the following at the command prompt ::
Kens-MBP:~ $ docker pull tehatanlgov/ptolemy:latest
This downloads the ptolemy container from docker hub. You need do this only
once. List downloaded containers with ::
Kens-MBP:~ $ docker images
- Make the ptolemy script executable by typing ::
Kens-MBP:~ $ chmod 755 /path/to/ptolemy
where /path/to/ptolemy is the full path to where you saved the file. Eg, if
you saved it to your home directory, then it would be ~/ptolemy. To run
ptolemy, type ::
Kens-MBP:~ $ /path/to/ptolemy < something.in > something.out
Tip: I suggest you save the ptolemy script to your home bin directory. Make
one if you don't have one and add the bin directory path to your PATH
environment variable. Ask if you need instructions on how to do this. After
you do this, then all you need to say is ::
Kens-MBP:~ $ ptolemy < another.in > another.out

507
Cleopatra/knockout.C Normal file
View File

@ -0,0 +1,507 @@
#include "HELIOS_LIB.h"
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TF1.h"
#include "TTree.h"
#include "TRandom.h"
#include "TClonesArray.h"
#include <vector>
#include <fstream>
//----------- usage
// $root sim.C+ | tee output.txt
// this will same the massage to output.txt
void knockout(){
//================================================= User Setting
//---- reaction
int AA = 23, ZA = 9;
int Aa = 1, Za = 1;
int A2 = 1, Z2 = 1;
bool isNormalKinematics = false;
bool isOverRideExNegative = true;
double maxkb = 200.;
//---- beam
double KEAmean = 100; // MeV/u
const int nKEA = 1;
double KEAList[nKEA] = {300};
double KEAsigma = 0; //KEAmean*0.001; // MeV/u , assume Guassian
double thetaMean = 0.; // mrad
double thetaSigma = 0.; // mrad , assume Guassian due to small angle
int numEvent = 100000;
//---- HELIOS detector geometry
string heliosDetGeoFile = "";//"detectorGeo_upstream.txt";
double BField = 4.0; // if not detector, must set B-field, else, this value is not used.
double BFieldTheta = 0.; // direction of B-field
double eSigma = 0.0001 ; // detector energy sigma MeV
double zSigma = 0.1 ; // detector position sigma mm
//---- excitation of Beam
int nExA = 1;
double ExAList[nExA];
ExAList[0] = 0.000; // MeV
//ExAList[1] = 1.567;
//---- Separation energy
string separationFile = "separation_energies.txt";
//---- save root file name
TString saveFileName = "knockout.root";
//---- Auxiliary setting
bool isTargetScattering = false;
bool isDecay = false;
bool isReDo = false; // redo calculation until detected.
//---- target
double density = 0.913; // 0.913 g/cm3
double targetThickness = 2./2. * 2.2e-4; // 2.2 um = 201 ug/cm2
string stoppingPowerForA = "208Pb_in_CD2.txt"; // generated by SRIM
string stoppingPowerForb = "1H_in_CD2.txt";
string stoppingPowerForB = "209Pb_in_CD2.txt";
//=============================================================
//=============================================================
//=============================================================
//===== Set Reaction
Knockout reaction;
int AB = AA-A2, ZB = ZA-Z2;
int A1 = Aa , Z1 = Za;
reaction.SetA(AA,ZA);
reaction.Seta(Aa,Za);
reaction.Set2(A2,Z2);
reaction.SetIncidentEnergyAngle(KEAmean, 0, 0);
reaction.OverRideExNegative(isOverRideExNegative);
printf("===================================================\n");
printf("=========== %s ===========\n", reaction.GetReactionName().Data());
printf("=========== KE: %9.4f +- %5.4f MeV/u, dp/p = %5.2f %% \n", KEAmean, KEAsigma, KEAsigma/KEAmean * 50.);
printf("======== theta: %9.4f +- %5.4f MeV/u \n", thetaMean, thetaSigma);
printf("===================================================\n");
//======== Set HELIOS
printf("############################################## HELIOS configuration\n");
HELIOS helios1; // for particle-1
HELIOS helios2; // for particle-2
helios1.SetMagneticFieldDirection(BFieldTheta);
helios2.SetMagneticFieldDirection(BFieldTheta);
bool sethelios1 = helios1.SetDetectorGeometry(heliosDetGeoFile);
bool sethelios2 = helios2.SetDetectorGeometry(heliosDetGeoFile);
if( sethelios1 && sethelios2 ) {
int mDet = helios1.GetNumberOfDetectorsInSamePos();
printf("========== energy resol.: %f MeV\n", eSigma);
printf("=========== pos-Z resol.: %f mm \n", zSigma);
}else{
helios1.SetMagneticField(BField);
helios2.SetMagneticField(BField);
printf("======== B-field : %5.2f T, Theta: %5.2f deg\n", BField, BFieldTheta);
}
//==== Target scattering, only energy loss
if(isTargetScattering) printf("############################################## Target Scattering\n");
TargetScattering msA;
TargetScattering msB;
TargetScattering msb;
if(isTargetScattering) printf("======== Target : (thickness : %6.2f um) x (density : %6.2f g/cm3) = %6.2f ug/cm2\n",
targetThickness * 1e+4,
density,
targetThickness * density * 1e+6);
if( isTargetScattering ){
msA.LoadStoppingPower(stoppingPowerForA);
msb.LoadStoppingPower(stoppingPowerForb);
msB.LoadStoppingPower(stoppingPowerForB);
}
//======= Decay of particle-B
Decay decay;
decay.SetMotherDaugther(AB, ZB, AB-1,ZB); //neutron decay
//======= loading excitation energy
printf("############################################## excitation energies\n");
vector<double> SpList;
printf("----- loading separation energies.");
ifstream file;
file.open(separationFile.c_str());
string isotopeName;
if( file.is_open() ){
string line;
int i = 0;
while( file >> line){
//printf("%d, %s \n", i, line.c_str());
if( line.substr(0,2) == "//" ) continue;
if( i == 0 ) isotopeName = line;
if ( i >= 1 ){
SpList.push_back(atof(line.c_str()));
}
i = i + 1;
}
file.close();
printf("... done.\n");
printf("========== %s\n", isotopeName.c_str());
int n = SpList.size();
for(int i = 0; i < n ; i++){
if( isDecay ) {
TLorentzVector temp(0,0,0,0);
int decayID = decay.CalDecay(temp, SpList[i], 0);
if( decayID == 1) {
printf("%d, Sp: %6.2f MeV --> Decay. \n", i, SpList[i]);
}else{
printf("%d, Sp: %6.2f MeV\n", i, SpList[i]);
}
}else{
printf("%d, Sp: %6.2f MeV \n", i, SpList[i]);
}
}
}else{
printf("... fail\n");
return;
}
//====================== build tree
TFile * saveFile = new TFile(saveFileName, "recreate");
TTree * tree = new TTree("tree", "tree");
double thetaNN, phiNN;
double theta1, phi1, T1;
double theta2, phi2, T2;
double thetaB, TB;
double Sp, kb, thetab, phib;
int SpID;
double ExA;
int ExAID;
double KEA, KEAscattered, theta, phi;
double mB,mb;
tree->Branch("theta1", &theta1, "theta1/D");
tree->Branch("phi1", &phi1, "phi1/D");
tree->Branch("T1", &T1, "T1/D");
tree->Branch("theta2", &theta2, "theta2/D");
tree->Branch("phi2", &phi2, "phi2/D");
tree->Branch("T2", &T2, "T2/D");
tree->Branch("thetaB", &thetaB, "thetaB/D");
tree->Branch("TB", &TB, "TB/D");
tree->Branch("thetaNN", &thetaNN, "thetaNN/D");
tree->Branch("phiNN", &phiNN, "phiNN/D");
tree->Branch("Sp", &Sp, "Sp/D");
tree->Branch("kb", &kb, "kb/D");
tree->Branch("thetab", &thetab, "thetab/D");
tree->Branch("phib", &phib, "phib/D");
tree->Branch("SpID", &SpID, "SpID/I");
tree->Branch("ExAID", &ExAID, "ExAID/I");
tree->Branch("KEA", &KEA, "KEA/D");
if(isTargetScattering) tree->Branch("KEAscattered", &KEAscattered, "KEAscattered/D");
tree->Branch("theta", &theta, "theta/D");
tree->Branch("phi", &phi, "phi/D");
tree->Branch("mB", &mB, "mB/D");
tree->Branch("mb", &mb, "mb/D");
tree->Branch("Bfield", &BField, "Bfield/D");
TClonesArray * arr = new TClonesArray("TLorentzVector");
tree->Branch("fV", arr, 256000);
arr->BypassStreamer();
TClonesArray * arrN = new TClonesArray("TLorentzVector");
tree->Branch("fVN", arrN, 256000);
arrN->BypassStreamer();
TLorentzVector* fourVector = NULL;
// the output of Helios.CalHit
double e1, z1, t1, rho1, x1h, y1h, r1h;
double e2, z2, t2, rho2, x2h, y2h, r2h;
tree->Branch("e1", &e1, "e1/D");
tree->Branch("z1", &z1, "z1/D");
tree->Branch("t1", &t1, "t1/D");
tree->Branch("x1h", &x1h, "x1h/D"); //at 200 mm downstream
tree->Branch("y1h", &y1h, "y1h/D");
tree->Branch("r1h", &r1h, "r1h/D");
tree->Branch("rho1", &rho1, "rho1/D");
tree->Branch("e2", &e2, "e2/D");
tree->Branch("z2", &z2, "z2/D");
tree->Branch("t2", &t2, "t2/D");
tree->Branch("x2h", &x2h, "x2h/D"); //at 200 mm downstream
tree->Branch("y2h", &y2h, "y2h/D");
tree->Branch("r2h", &r2h, "r2h/D");
tree->Branch("rho2", &rho2, "rho2/D");
double Sp2;
tree->Branch("Sp2", &Sp2, "Sp2/D");
//different coordinate
double a0, a1, a2; // in k1-k2 coordinate
tree->Branch("a0", &a0, "a0/D");
tree->Branch("a1", &a1, "a1/D");
tree->Branch("a2", &a2, "a2/D");
double b; // in THREEDEE coordinate
tree->Branch("b", &b, "b/D");
//========timer
TBenchmark clock;
bool shown ;
clock.Reset();
clock.Start("timer");
shown = false;
printf("############################################## generating %d events \n", numEvent);
//====================================================== calculate
int count = 0;
for( int i = 0; i < numEvent; i++){
bool redoFlag = true;
if( !isReDo ) redoFlag = false;
do{
//==== Set Ex of A
ExAID = gRandom->Integer(nExA);
ExA = ExAList[ExAID];
reaction.SetExA(ExA);
//==== Set Sp of B
SpID = gRandom->Integer(SpList.size());
Sp = SpList[SpID];
//==== Set incident beam
if( KEAsigma == 0 ){
KEA = KEAmean;
}else{
KEA = gRandom->Gaus(KEAmean, KEAsigma);
}
if( thetaSigma == 0 ){
theta = thetaMean;
}else{
theta = gRandom->Gaus(thetaMean, thetaSigma);
}
int rKEA = gRandom->Integer(nKEA);
KEA = KEAList[rKEA];
/*
KEA = 300*gRandom->Rndm();
BField = 4*gRandom->Rndm();
helios1.SetMagneticField(BField);
helios2.SetMagneticField(BField);
*/
reaction.SetIncidentEnergyAngle(KEA, theta, 0.);
/*
//For target scattering
reaction.CalIncidentChannel(isNormalKinematics); // but only need is PA
TLorentzVector PA = reaction.GetPA();
double depth = 0;
if( isTargetScattering ){
//==== Target scattering, only energy loss
depth = targetThickness * gRandom->Rndm();
msA.SetTarget(density, depth);
TLorentzVector PAnew = msA.Scattering(PA);
KEAscattered = msA.GetKE()/AA;
reaction.SetIncidentEnergyAngle(KEAscattered, theta, phi);
}*/
//==== Calculate reaction
thetaNN = TMath::ACos(2 * gRandom->Rndm() - 1) ;
phiNN = TMath::TwoPi() * gRandom->Rndm();
kb = maxkb * gRandom->Rndm();
thetab = TMath::ACos(2 * gRandom->Rndm() - 1);
phib = TMath::TwoPi() * gRandom->Rndm();
reaction.SetBSpk(Sp, kb, thetab, phib);
reaction.CalReactionConstant(isNormalKinematics);
reaction.Event(thetaNN, phiNN);
TLorentzVector PA = reaction.GetPA();
TLorentzVector Pa = reaction.GetPa();
TLorentzVector P1 = reaction.GetP1();
TLorentzVector P2 = reaction.GetP2();
TLorentzVector Pb = reaction.GetPb();
TLorentzVector PB = reaction.GetPB();
/*
//==== Calculate energy loss of scattered and recoil in target
if( isTargetScattering ){
if( Pb.Theta() < TMath::PiOver2() ){
msb.SetTarget(density, targetThickness - depth);
}else{
msb.SetTarget(density, depth);
}
Pb = msb.Scattering(Pb);
TbLoss = msb.GetKELoss();
msB.SetTarget(density, targetThickness - depth);
PB = msB.Scattering(PB);
}else{
TbLoss = 0;
}
//======= Decay of particle-B
if( isDecay){
int decayID = decay.CalDecay(PB, Sp, 0); // decay to ground state
if( decayID == 1 ){
PB = decay.GetDaugther_D();
decayTheta = decay.GetAngleChange();
}else{
decayTheta = TMath::QuietNaN();
}
}
*/
//################################### tree branches
//===== reaction
theta1 = P1.Theta() * TMath::RadToDeg();
theta2 = P2.Theta() * TMath::RadToDeg();
thetaB = PB.Theta() * TMath::RadToDeg();
T1 = P1.E() - P1.M();
T2 = P2.E() - P2.M();
TB = PB.E() - PB.M();
phi1 = P1.Phi() * TMath::RadToDeg();
phi2 = P2.Phi() * TMath::RadToDeg();
//--------- diff coordinate
b = TMath::ASin( TMath::Sin(P2.Theta()) * TMath::Sin(P1.Phi() - P2.Phi() )) * TMath::RadToDeg();
TVector3 kA = (PA.Vect()).Unit();
TVector3 k1 = (P1.Vect()).Unit();
TVector3 k2 = (P2.Vect()).Unit();
TVector3 n = (k1.Cross(k2)).Unit();
TVector3 j = (kA - (kA.Dot(n))*n).Unit();
a0 = TMath::ASin(n.Dot(kA)) * TMath::RadToDeg();
a1 = TMath::ACos(k1.Dot(j)) * TMath::RadToDeg();
a2 = TMath::ACos(k2.Dot(j)) * TMath::RadToDeg();
//----------- mass
mB = PB.M();
mb = Pb.M();
TVector3 bA = PA.BoostVector();
for(int i = 0; i < 6 ; i++){
TLorentzVector temp;
double xyzt[4];
switch(i){
case 0: temp = PA; break;
case 1: temp = Pa; break;
case 2: temp = P1; break;
case 3: temp = P2; break;
case 4: temp = PB; break;
case 5: temp = Pb; break;
}
temp.GetXYZT(xyzt);
fourVector = (TLorentzVector*) arr->ConstructedAt(i);
fourVector->SetXYZT(xyzt[0], xyzt[1], xyzt[2], xyzt[3]);
//into normal kinematic
temp.Boost(-bA);
temp.GetXYZT(xyzt);
fourVector = (TLorentzVector*) arrN->ConstructedAt(i);
fourVector->SetXYZT(xyzt[0], xyzt[1], xyzt[2], xyzt[3]);
}
//==== Helios
int hit1 = helios1.CalHit(P1, Z1, PB, ZB);
int hit2 = helios2.CalHit(P2, Z2, PB, ZB);
double recoilZ = 1000;
e1 = helios1.GetEnergy() + gRandom->Gaus(0, eSigma);
z1 = helios1.GetZ() ;
t1 = helios1.GetTime();
rho1 = helios1.GetRho();
x1h = helios1.GetXPos(recoilZ);
y1h = helios1.GetYPos(recoilZ);
r1h = helios1.GetR(recoilZ);
e2 = helios2.GetEnergy() + gRandom->Gaus(0, eSigma);
z2 = helios2.GetZ() ;
t2 = helios2.GetTime();
rho2 = helios2.GetRho();
x2h = helios2.GetXPos(recoilZ);
y2h = helios2.GetYPos(recoilZ);
r2h = helios2.GetR(recoilZ);
double gammaA = PA.Gamma();
double betaA = PA.Beta();
Sp2 = (1-gammaA)* 938.272 - gammaA*(e1+e2) + betaA*gammaA*(P1.Pz() + P2.Pz()) ;
//printf("%f, %f | %f, %f \n", e1, z1, e2, z2);
//change thetaNN into deg
thetaNN = thetaNN * TMath::RadToDeg();
if( hit1 == 1) {
count ++;
}
if( isReDo ){
if( hit1 == 1) {
redoFlag = false;
}else{
redoFlag = true;
//printf("%d, %2d, thetaNN : %f, theta : %f, z0: %f \n", i, hit, thetaNN * TMath::RadToDeg(), thetab, helios.GetZ0());
}
}else{
redoFlag = false;
}
}while( redoFlag );
tree->Fill();
//#################################################################### Timer
clock.Stop("timer");
Double_t time = clock.GetRealTime("timer");
clock.Start("timer");
if ( !shown ) {
if (fmod(time, 10) < 1 ){
printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
shown = 1;
}
}else{
if (fmod(time, 10) > 9 ){
shown = 0;
}
}
}
saveFile->Write();
saveFile->Close();
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
gROOT->ProcessLine(".q");
/**/
}

33
Cleopatra/makefile Normal file
View File

@ -0,0 +1,33 @@
CC=g++
all: Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray FindThetaCM Transfer PlotSimulation IsotopeShort
#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`
PlotSimulation: PlotSimulation.C ../Armory/Check_Simulation.C
$(CC) PlotSimulation.C -o PlotSimulation `root-config --cflags --glibs`
Isotope: ../Cleopatra/Isotope.h ../Cleopatra/Isotope.C
$(CC) Isotope.C -o Isotope
IsotopeShort: ../Cleopatra/Isotope.h ../Cleopatra/IsotopeShort.C
$(CC) IsotopeShort.C -o IsotopeShort

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

153
Cleopatra/nuclear_data.py Executable file
View File

@ -0,0 +1,153 @@
#!/usr/local/bin/python3
################################################
import sys
if len(sys.argv) == 1 :
print( "usage :" )
print( "%s AZ [maxEx]" % sys.argv[0] )
print( "e.g. %s 17O --> 17O ground state data" % sys.argv[0])
print( "e.g. %s 20O 5 --> 20O ground state data and excited state < 5 MeV" % sys.argv[0])
exit()
################################################
import pandas as pd
# the API webpage
# https://www-nds.iaea.org/relnsd/vcharthtml/api_v0_guide.html#examples
# the service URL
livechart = "https://nds.iaea.org/relnsd/v0/data?"
import urllib.request
def lc_read_csv(url):
req = urllib.request.Request(url)
req.add_header('User-Agent', 'Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:77.0) Gecko/20100101 Firefox/77.0')
return pd.read_csv(urllib.request.urlopen(req))
haha = lc_read_csv(livechart + 'fields=ground_states&nuclides=all')
mp = 938.27208816; #MeV/c^2
mn = 939.56542052;
def FindZ(AZ):
query = livechart + "fields=ground_states&nuclides=" + AZ
temp = lc_read_csv(query);
try :
return temp['z']
except :
return 'na'
def FindSym(Z):
try:
return (haha['symbol'][haha['z']==Z]).iloc[0]
except:
return 'na'
def Mass(A, Z):
try :
BEA = float(haha['binding'][haha['z']==Z][haha['n']==(A-Z)])/1000
return (A-Z)*mn + Z*mp - A * BEA
except :
return -404
def MassSym(AZ):
query = livechart + "fields=ground_states&nuclides=" + AZ
temp = lc_read_csv(query);
Z = temp['z']
N = temp['n']
try :
return Z*mp + N*mn - (Z+N)*temp['binding']/1000
except:
return -404
def Sp(A,Z,a,z):
mA = Mass(A,Z)
mB = Mass(A-a, Z-z)
if z == 0 :
mb = a * mn
elif a == z :
mb = a * mp
else :
mb = Mass(a,z)
if (mB == -404 or mb == -404 or mA == -404) :
return -404
else:
return mB + mb - mA
def Ex(AZ, maxMeV):
query = livechart + "fields=levels&nuclides=" + AZ
tempEx = lc_read_csv(query);
try :
return tempEx[['energy', 'jp']][tempEx['energy']<= maxMeV * 1000]
except:
return -404
def Info(AZ):
query = livechart + "fields=ground_states&nuclides=" + AZ
temp = lc_read_csv(query);
print("============================================== " + AZ )
try :
Z = temp['z'][0]
N = temp['n'][0]
mass = Z*mp + N*mn - (Z+N)*temp['binding']/1000
halfLife = temp['half_life_sec'][0]
print(" A : %3d, Z : %3d, N : %3d, Mass : %.4f MeV" % (Z+N, Z, N, mass))
print("Jpi : %3s, half-live : %s sec" % (temp['jp'][0], halfLife))
print("Sn : %8.3f MeV, Sp : %8.3f MeV" % (Sp(Z+N,Z, 1, 0), Sp(Z+N,Z, 1, 1)))
print("S2n : %8.3f MeV, S2p : %8.3f MeV, Sd : %8.3f MeV" % (Sp(Z+N,Z, 2, 0), Sp(Z+N,Z, 2, 2), Sp(Z+N, Z, 2, 1)))
print("S3n : %8.3f MeV, S3p : %8.3f MeV, St : %8.3f MeV, S(3He) : %8.3f MeV" % (Sp(Z+N,Z, 3, 0), Sp(Z+N,Z, 3, 3), Sp(Z+N, Z, 3, 1), Sp(Z+N, Z, 3, 2)))
print("S4n : %8.3f MeV, S4p : %8.3f MeV, Sa : %8.3f MeV" % (Sp(Z+N,Z, 4, 0), Sp(Z+N,Z, 4, 4), Sp(Z+N, Z, 4, 2)))
print(" magnetic dipole : " + temp['magnetic_dipole'][0] + " mu.N")
print("electric quadruple : " + temp['electric_quadrupole'][0] + " barn")
if halfLife > 0 :
print('------------ decay mode:')
for i in range(1, 4) :
print("%5s %s %%" % (temp["decay_%d" % i][0], temp["decay_%d_%%" % i][0]))
print('--------------------------')
except :
dummy = 1
print("====================================================\n")
#outfile = open("%s.txt" % AZ, "w")
#outfile.write("============================================== " + AZ + "\n")
#try :
# Z = temp['z'][0]
# N = temp['n'][0]
# mass = Z*mp + N*mn - (Z+N)*temp['binding']/1000
# halfLife = temp['half_life_sec'][0]
# outfile.write(" A : %3d, Z : %3d, N : %3d, Mass : %.4f MeV \n" % (Z+N, Z, N, mass))
# outfile.write("Jpi : %3s, half-live : %s sec \n" % (temp['jp'][0], halfLife))
# outfile.write("Sn : %8.3f MeV, Sp : %8.3f MeV \n" % (Sp(Z+N,Z, 1, 0), Sp(Z+N,Z, 1, 1)))
# outfile.write("S2n : %8.3f MeV, S2p : %8.3f MeV, Sd : %8.3f MeV \n" % (Sp(Z+N,Z, 2, 0), Sp(Z+N,Z, 2, 2), Sp(Z+N, Z, 2, 1)))
# outfile.write("S3n : %8.3f MeV, S3p : %8.3f MeV, St : %8.3f MeV, S(3He) : %8.3f MeV \n" % (Sp(Z+N,Z, 3, 0), Sp(Z+N,Z, 3, 3), Sp(Z+N, Z, 3, 1), Sp(Z+N, Z, 3, 2)))
# outfile.write("S4n : %8.3f MeV, S4p : %8.3f MeV, Sa : %8.3f MeV \n" % (Sp(Z+N,Z, 4, 0), Sp(Z+N,Z, 4, 4), Sp(Z+N, Z, 4, 2)))
# outfile.write(" magnetic dipole : " + temp['magnetic_dipole'][0] + " mu.N\n")
# outfile.write("electric quadruple : " + temp['electric_quadrupole'][0] + " barn \n")
# if halfLife > 0 :
# outfile.write('------------ decay mode: \n')
# for i in range(1, 4) :
# outfile.write("%5s %s %% \n" % (temp["decay_%d" % i][0], temp["decay_%d_%%" % i][0]))
# outfile.write('-------------------------\n')
#except :
# dummy = 1
#outfile.write("====================================================\n")
#outfile.close();
#
#infile = open("%s.txt" % AZ, "r")
#Lines = infile.readlines()
#
#for line in Lines:
# print(line.strip("\n"))
#
#infile.close();
################################################
AZ = sys.argv[1];
maxEx = 0;
if len(sys.argv) > 2 :
maxEx = float(sys.argv[2]);
Info(AZ)
if maxEx > 0 :
print( Ex(AZ, float(maxEx)) )

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.

2
Cleopatra/ptolemy_mac Executable file
View File

@ -0,0 +1,2 @@
#!/bin/bash
docker run --rm -i -w $PWD tehatanlgov/ptolemy

77
Cleopatra/transfer_test.C Normal file
View File

@ -0,0 +1,77 @@
#include "HELIOS_LIB.h"
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TF1.h"
#include "TTree.h"
#include "TRandom.h"
#include "TGraph.h"
#include "TMacro.h"
#include <stdlib.h>
#include <vector>
#include <fstream>
#include <TObjArray.h>
void transfer_test(double t, double p, double bField, bool fromOutSide){
TransferReaction reaction;
reaction.SetA(14, 6);
reaction.Seta( 2, 1);
reaction.Setb( 1, 1);
reaction.SetB(15, 6);
reaction.SetExB(0);
reaction.SetIncidentEnergyAngle(10, 0, 0);
reaction.CalReactionConstant();
HELIOS helios;
helios.SetDetectorGeometry("../working/detectorGeo.txt");
helios.OverrideMagneticField(bField);
helios.SetDetectorOutside(fromOutSide);
double beta = reaction.GetReactionBeta() ;
double slope = 299.792458 * abs(helios.GetBField()) / TMath::TwoPi() * beta / 1000.; // MeV/mm
double alpha = slope / beta;
printf("===================================\n");
printf("Mass A : %8.2f MeV/c2\n", reaction.GetMass_A());
printf("Mass a : %8.2f MeV/c2\n", reaction.GetMass_a());
printf("Mass b : %8.2f MeV/c2\n", reaction.GetMass_b());
printf("Mass B : %8.2f MeV/c2\n", reaction.GetMass_B());
printf("CM Mass : %8.2f MeV\n", reaction.GetCMTotalEnergy());
printf("CM beta : %8.6f \n", beta);
printf("slope : %8.6f MeV\n", alpha * beta);
printf("alpha : %8.6f MeV\n", alpha);
double thetaCM = t * TMath::DegToRad();
double phiCM = - p * TMath::DegToRad();
TLorentzVector * output = reaction.Event(thetaCM, phiCM);
TLorentzVector Pb = output[2];
TLorentzVector PB = output[3];
Pb.Print();
PB.Print();
helios.CalArrayHit(Pb, 1);
helios.CalRecoilHit(PB, 6);
printf("+++++++++++++++++++++++++++++++++++++\n");
int hitID = 2;
while( hitID > 1 ){
printf("==================== check accp.\n");
hitID = helios.DetAcceptance();
printf("-------------------- hitID %d\n", hitID);
}
PrintTrajectory(helios.GetTrajectory_b());
}

483
armory/AnalysisLib.h
View File

@ -6,6 +6,9 @@
#include <string>
#include <TMacro.h>
int FindDetType(int detID, std::vector<int> detMaxID){
for( int k = 0; k < (int) detMaxID.size(); k++){
int low = (k == 0 ? 0 : detMaxID[k-1]);
@ -61,6 +64,486 @@ void PrintMapping(std::vector<std::vector<int>> mapping, std::vector<std::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

2848
armory/AutoFit.C Normal file
View File

File diff suppressed because it is too large Load Diff

30
working/Check_Simulation_Config.txt Normal file
View File

@ -0,0 +1,30 @@
/////enum plotID { pEZ, /// 0
///// pRecoilXY, /// 1
///// pRecoilXY1, /// 2
///// pRecoilXY2, /// 3
///// pRecoilRZ, /// 4
///// pRecoilRTR, /// 5
///// pTDiffZ, /// 6
///// pThetaCM, /// 7
///// pThetaCM_Z, /// 8
///// pExCal, /// 9
///// pRecoilRThetaCM, /// 10
///// pArrayXY, /// 11
///// pInfo, /// 12
///// pHitID, /// 13
///// pElum1XY, /// 14
///// pEElum1R, /// 15
///// pElum1RThetaCM, /// 16
///// pEmpty }; /// 17
/////=============================================== User Config
{pEZ, pExCal, pThetaCM, pRecoilRZ, break, pThetaCM_Z, pRecoilXY, pInfo, pRecoilRThetaCM} //Canvas config
hit == 1 && loop <= 1 && thetaCM > 10
60 //elum range
{0,50} //thetaCM range
false //shownKELines
false //isOverRideEx
{-0.5, 4.0} // over-rdied Ex range
///============================== example of gate
hit == 1 && loop <= 1 && thetaCM > 10 && detRowID == 0
hit == 1 && loop <= 1
15 < rhoElum1 && rhoElum1 < 50 && rhoElum2 > 60

46
working/DWBA Normal file
View File

@ -0,0 +1,46 @@
#========= Input for Cleopatra
#===== # for comment line, must be at the beginning of line
# This is a simplified input for creating Ptolmey in-file.
# Ptolemy has many control parameters for specific situation. Use with caution or ask Ben Kay.
#===== 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
#36Ar(d,a)34Cl 0 4L=2 3+ 0.000 8MeV/u As # (d,a) reaction

6
working/Ex.txt Normal file
View File

@ -0,0 +1,6 @@
//Ex relative_xsec SF sigma_in_MeV
//<--- use "//" for line comment
0.000 1.0 1.0 0.0100
//4.400 1.0 1.0 0.0100
//4.600 1.0 1.0 0.0100
#============_End_of_file

27
working/detectorGeo.txt Normal file
View File

@ -0,0 +1,27 @@
-3.00 //Bfield_[T]
0.00 //Bfield_direction_to_z-axis_[deg]_should_not_use
462.5 //bore_[mm]
11.5 //distance_from_axis_[mm]
10.0 //width_of_detector_[mm]
50 //length_of_detector_[mm]
1000 //recoil_position_+_for_downstream_[mm]
10.0 //inner_radius_of_recoil_detector_[mm]
40.2 //outter_radius_of_recoil_detector_[mm]
false //is_coincident_with_recoil
0 //Recoil_1_position_[mm]_when_0_disable_tree_branch
0 //Recoil_2_position_[mm]
0.00 //Elum_1_position_[mm]_(just_another_recoil_detector_but_for_light_recoil)
0.00 //Elum_2_position_[mm]_when_Elum=0_disable_tree_branch
0 //support_length_[mm]
-121 //first_position_-_for_upstream_[mm]
0.03 //energy_resolution_of_PSD_array_[MeV]
1.00 //position_resolution_of_PSD_array_[mm]
Out //detector_facing_Out_or_In
4 //number_of_detector_as_same_side
0.00 //1st_detector_near_position_in_reference_to_det6_[mm]
58.6 //2nd_det
117.9
176.8
235.8 //5th_det
290.0
#============= end of file

25
working/reactionConfig.txt Normal file
View File

@ -0,0 +1,25 @@
32 //beam_A
14 //beam_Z
2 //target_A
1 //target_Z
1 //recoil_light_A
1 //recoil-light_Z
8.8 //beam-energy_in_MeV/u
0.000 //beam-energy_sigma_in_MeV/u
0.000 //beam-angle_in_mrad
0.000 //beam-emittance_in_mrad
0.00 //x_offset_of_Beam_in_mm
0.00 //y_offset_of_Beam_in_mm
100000 //number_of_Event_being_generated
false //isTargetScattering
0.913 //target_density_in_g/cm3
2.2e-4 //targetThickness_in_cm
../SRIM/20F_in_CD2.txt //stopping_power_for_beam
../SRIM/3H_in_CD2.txt //stopping_power_for_light_recoil
../SRIM19F_in_CD2.txt //stopping_power_for_heavy_recoil
false //isDacay
32 //decayNucleus_A
14 //decayNucleus_Z
false //isReDo
0.0 //excitation_energy_of_A[MeV]
#===== end of file