update the Check_Simulation.C

This commit is contained in:
Ryan Tang 2024-02-20 18:15:16 -05:00
parent 364530f73c
commit 7a5107998a
10 changed files with 572 additions and 603 deletions

2
.gitignore vendored
View File

@ -25,4 +25,4 @@ Cleopatra/IsotopeShort
Cleopatra/PlotSimulation
Cleopatra/PlotTGraphTObjArray
Cleopatra/SimAlpha
Cleopatra/Transfer
Cleopatra/SimTransfer

View File

@ -95,7 +95,9 @@ public:
bool LoadDetectorGeo(TString fileName, bool verbose = true);
bool LoadDetectorGeo(TMacro * macro, bool verbose = true);
void Print( bool printAll = true) const;
void PrintWithoutArray() ;
void Print( bool printAll = true) ;
private:
@ -203,24 +205,14 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
}
inline void DetGeo::Print(bool printAll) const{
inline void DetGeo::PrintWithoutArray(){
printf("=====================================================\n");
printf(" B-field: %8.2f T, Theta : %6.2f deg \n", Bfield, BfieldTheta);
if( BfieldTheta != 0.0 ) {
printf(" +---- field angle != 0 is not supported!!! \n");
}
printf(" B-field: %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan");
printf(" B-field Theta : %6.2f deg \n", BfieldTheta);
if( 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", recoilPos, recoilInnerRadius, recoilOuterRadius);
for( int i = 0; i < 2 ; i++){
if( printAll || array[i].enable ) {
printf("-----------------------------------%d-th Detector Position \n", i);
array[i].PrintArray();
}
}
if( elumPos1 != 0 || elumPos2 != 0 || recoilPos1 != 0 || recoilPos2 != 0){
printf("=================================== Auxillary/Imaginary Detectors\n");
}
@ -229,6 +221,21 @@ inline void DetGeo::Print(bool printAll) const{
if( recoilPos1 != 0 ) printf(" Recoil 1 pos.: %f mm \n", recoilPos1);
if( recoilPos2 != 0 ) printf(" Recoil 2 pos.: %f mm \n", recoilPos2);
printf("=====================================================\n");
}
inline void DetGeo::Print(bool printAll){
PrintWithoutArray();
for( int i = 0; i < 2 ; i++){
if( printAll || array[i].enable ) {
printf("-----------------------------------%d-th Detector Position \n", i);
array[i].PrintArray();
}
}
printf("=====================================================\n");
}

View File

@ -114,7 +114,7 @@ public:
bool LoadReactionConfig(TString fileName);
bool LoadReactionConfig(TMacro * macro);
void Print() const;
void Print(int ID = -1, bool withEx = true) const;
private:
@ -237,7 +237,7 @@ inline bool ReactionConfig::LoadReactionConfig(TMacro * macro){
return true;
}
inline void ReactionConfig::Print() const{
inline void ReactionConfig::Print(int ID, bool withEx) const{
printf("=====================================================\n");
@ -261,8 +261,10 @@ inline void ReactionConfig::Print() const{
for( int i = 0; i < 2; i ++ ){
printf("------------------------------ Recoil-%d\n", i);
if( ID == i || ID < 0 ){
recoil[i].Print();
exList[i].Print();
if( withEx ) exList[i].Print();
}
}

View File

@ -22,6 +22,7 @@
#include "../Armory/ClassDetGeo.h"
#include "../Armory/ClassReactionConfig.h"
#include "../Cleopatra/ClassIsotope.h"
#include "../Cleopatra/ClassTransfer.h"
double * FindRange(TString branch, TString gate, TTree * tree, double output[2]);
double ExtractNumber(int index, TMacro * macro);
@ -49,7 +50,9 @@ enum plotID { pEZ, /// 0
pElum1RThetaCM, /// 16
pEmpty }; /// 17
plotID StringToPlotID(TString str);
void Check_Simulation(TString filename = "transfer1.root",
//*=====================================================
void Check_Simulation(TString filename = "transfer.root",
TString configFile = "../working/Check_Simulation_Config.txt",
Int_t padSize = 500,
bool outputCanvas = false){
@ -77,7 +80,6 @@ void Check_Simulation(TString filename = "transfer1.root",
printf("%s \n", gate.Data());
///==== config Canvas
vector<TString> plotConfig = StringToVector( ExtractString(startLineNum, config));
vector<plotID> canvas;
@ -96,30 +98,15 @@ void Check_Simulation(TString filename = "transfer1.root",
}
if( colCount == 0 ) colCount = colCount_new;
///printf("plot row: %d, col: %d \n", rowCount, colCount);
//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;
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);
TObjArray * fxList = (TObjArray *) file->FindObjectAny("EZCure");
TObjArray * txList = (TObjArray *) file->FindObjectAny("EZ_thetaCM");
gStyle->SetOptStat("");
gStyle->SetStatY(0.9);
@ -129,138 +116,121 @@ void Check_Simulation(TString filename = "transfer1.root",
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];
//*================== detGeoID
TMacro * detGeoIDTxt = (TMacro *) file->FindObjectAny("detGeoID");
int detGeoID = atoi(detGeoIDTxt->GetListOfLines()->At(0)->GetName());
//================================== detetcor Geometry
//*================== reactionConfig
TMacro * reactionConfigTxt = (TMacro *) file->FindObjectAny("reactionConfig");
TString Reaction = reactionConfigTxt->GetName();
ReactionConfig reactionConfig(reactionConfigTxt);
Recoil recoil = reactionConfig.recoil[detGeoID];
int nEvent = reactionConfig.numEvents;
printf("number of events generated : %d \n", nEvent);
reactionConfig.Print(detGeoID, false);
//*================================== detetcor Geometry
printf("=================================\n");
printf(" loading detector Geometry.\n");
TMacro * detGeoTxt = (TMacro *) file->FindObjectAny("detGeo");
DetGeo detGeo;
detGeo.LoadDetectorGeo(detGeoTxt);
DetGeo detGeo(detGeoTxt);
Array array = detGeo.array[detGeoID];
Array array;
if( detGeo.use2ndArray){
array = detGeo.array2;
}else{
array = detGeo.array1;
}
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 = array.detPerpDist;
double length = array.detLength;
double posRecoil = detGeo.recoilPos;
double rhoRecoilIn = detGeo.recoilInnerRadius;
double rhoRecoilOut = detGeo.recoilOuterRadius;
double posRecoil1 = detGeo.recoilPos1;
double posRecoil2 = detGeo.recoilPos2;
vector<double> pos = array.detPos;
float firstPos = array.firstPos;
int rDet = array.nDet;
int cDet = array.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]);
}
}
detGeo.PrintWithoutArray();
array.PrintArray();
printf("=================================\n");
int numDet = rDet * cDet;
int numDet = array.nDet * array.mDet ;
double zRange[3] = {400, -1000, 1000}; /// zRange[0] = nBin
zRange[1] = array.zMin - 50;
zRange[2] = array.zMax + 50;
printf(" zRange : %f - %f \n", zRange[1], zRange[2]);
printf("=================================\n");
//========================================= Ex List;
printf(" loading Ex list\n");
//*========================================= Ex List;
double ExRange[2];
int numEx = 0;
ExcitedEnergies exList;
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 = AnalysisLib::SplitStr(temp, " ");
printf("%d | %s \n", i, tempStr[0].c_str());
exList.push_back(atof(tempStr[0].c_str()));
// if DEBA_ExList exist, use this, else use the recoil ExList
TMacro * exListTxt = (TMacro *) file->FindObjectAny("DWBA_ExList");
ExRange[0] = 9999999;
ExRange[1] = -9999999;
if( exListTxt == nullptr ){
exList = reactionConfig.exList[detGeoID];
numEx = exList.ExList.size();
for( size_t i = 0; i < numEx; i++ ){
double ex = exList.ExList[i].Ex;
if( ex < ExRange[0] ) ExRange[0] = ex;
if( ex > ExRange[1] ) ExRange[1] = ex;
}
double exSpan = exList.back() - exList[0];
}else{
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];
numEx = exListTxt->GetListOfLines()->GetSize();
for( int i = 0 ; i < numEx ; i++){
double ex = atof(exListTxt->GetListOfLines()->At(i)->GetName());
if( ex < ExRange[0] ) ExRange[0] = ex;
if( ex > ExRange[1] ) ExRange[1] = ex;
exList.Add(ex, 0, 0, 0);
}
printf("=================================\n");
}
//========================================= reaction parameters
double dExRange = ExRange[1] - ExRange[0];
ExRange[0] = ExRange[0] - dExRange * 0.1;
ExRange[1] = ExRange[1] + dExRange * 0.1;
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
printf("Number of Ex states = %d \n", numEx);
//*=================================== 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]);
TransferReaction transfer;
transfer.SetReactionSimple( reactionConfig.beamA, reactionConfig.beamZ, reactionConfig.targetA, reactionConfig.targetZ, recoil.lightA, recoil.lightZ, reactionConfig.beamEnergy);
double QQ = transfer.GetCMTotalEnergy();
double gamm = transfer.GetReactionGamma();
double mass = transfer.GetMass_b();
double slope = transfer.GetEZSlope( detGeo.Bfield);
double eRange[2] = {0, 10};
// double intercept = QQ/gamm - mass;
eRange[1] = zRange[2] * slope;
// printf("intercept of 0 MeV : %f MeV \n", intercept);
printf("eRange 0 MeV : %f MeV \n", eRange[1]);
double dERange = eRange[1] - eRange[0];
eRange[0] = eRange[0] - dERange * 0.1;
eRange[1] = eRange[1] + dERange * 0.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;
double momentum = transfer.GetMomentumbCM();
double beta = transfer.GetReactionBeta();
double alpha = slope / beta;
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]);
double recoilERange[2] = {0, 100};
//===================================================
printf("============================== Gate\n");
printf("gate : %s\n", gate.Data());
@ -275,83 +245,98 @@ void Check_Simulation(TString filename = "transfer1.root",
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]);
TH2F * hez = new TH2F("hez", Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", array.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++){
for( int i = 0; i < numEx ; 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);
TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius,
400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius);
tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz");
TArc * detArc1 = new TArc(0,0, rhoRecoilOut);
TArc * detArc1 = new TArc(0,0, detGeo.recoilOuterRadius);
detArc1->SetLineColor(kBlue-8);
detArc1->SetFillStyle(0);
detArc1->Draw("same");
TArc * detArc2 = new TArc(0,0, rhoRecoilIn);
TArc * detArc2 = new TArc(0,0, detGeo.recoilInnerRadius);
detArc2->SetLineColor(kBlue-8);
detArc2->SetFillStyle(0);
detArc2->Draw("same");
if( xBeam != 0. || yBeam != 0. ){
TArc * arc = new TArc(xBeam, yBeam, 1);
if( reactionConfig.beamX != 0. || reactionConfig.beamY != 0. ){
TArc * arc = new TArc(reactionConfig.beamX, reactionConfig.beamY, 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);
TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos1 ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius,
400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius);
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);
TH2F * hRecoilXY2 = new TH2F("hRecoilXY2", Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos2 ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius,
400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius);
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);
TH2F * hRecoilRZ = new TH2F("hRecoilRZ", "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", zRange[0], zRange[1], zRange[2], 400,0, detGeo.recoilOuterRadius);
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);
TH2F * hRecoilRTR = new TH2F("hRecoilRTR", "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", 500, recoilERange[0], recoilERange[1], 500, 0, detGeo.recoilOuterRadius);
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];
TH1F * hThetaCM[numEx];
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++){
for( int i = startID; i < numEx; 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]));
legend->AddEntry(hThetaCM[i], Form("Ex=%5.1f MeV", exList.ExList[i].Ex));
double max = hThetaCM[i]->GetMaximum();
if( max > maxCount ) maxCount = max;
}
for( int i = startID; i < nExID; i++){
for( int i = startID; i < numEx; i++){
hThetaCM[i]->GetYaxis()->SetRangeUser(1, maxCount * 1.2);
if( i == startID ) {
hThetaCM[i]->Draw();
@ -361,20 +346,21 @@ void Check_Simulation(TString filename = "transfer1.root",
}
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++){
for( int i = 0; i < numEx ; 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);
Isotope hRecoil(recoil.heavyA, recoil.heavyZ);
double Sn = hRecoil.CalSp(0,1);
double Sp = hRecoil.CalSp(1,0);
double Sa = hRecoil.CalSp2(4,2);
@ -401,14 +387,17 @@ void Check_Simulation(TString filename = "transfer1.root",
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);
TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, detGeo.recoilOuterRadius);
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);
TH2F * hArrayXY = new TH2F("hArrayXY", "Array-XY [gated]; X [mm]; Y [mm]", 400, -array.detPerpDist*1.5, array.detPerpDist*1.5, 400, -array.detPerpDist*1.5, array.detPerpDist*1.5);
tree->Draw("yArray:xArray>>hArrayXY", gate, "colz");
}
if( pID == pInfo ){
TLatex text;
text.SetNDC();
@ -417,7 +406,7 @@ void Check_Simulation(TString filename = "transfer1.root",
text.SetTextColor(2);
text.DrawLatex(0., 0.9, Reaction);
text.DrawLatex(0., 0.8, msg2);
text.DrawLatex(0., 0.8, detGeo.Bfield > 0 ? "out of plan" : "into plan");
text.SetTextColor(1);
text.DrawLatex(0., 0.7, "gate:");
@ -433,14 +422,13 @@ void Check_Simulation(TString filename = "transfer1.root",
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( reactionConfig.beamX != 0.0 || reactionConfig.beamY != 0.0 ){
text.DrawLatex(0.0, 0.1, Form("Bema pos: (%4.1f, %4.1f) mm", reactionConfig.beamX, reactionConfig.beamY));
}
}
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);
TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.elumPos1), 400, -elumRange, elumRange, 400, -elumRange, elumRange);
tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz");
double count = hElum1XY->GetEntries();
@ -529,7 +517,6 @@ void Check_Simulation(TString filename = "transfer1.root",
TH1F * hHit = new TH1F("hHit", "hit; hit-ID; count", 13, -11, 2);
tree->Draw("hit>>hHit", "", "");
}
///#######################################################
}

View File

@ -45,8 +45,7 @@ void Transfer(
TString ptolemyRoot = "DWBA.root",
TString saveFileName = "transfer.root"){
//############################################# Set Reaction
//*############################################# Set Reaction
TransferReaction transfer;
HELIOS helios;
Decay decay;
@ -65,7 +64,7 @@ void Transfer(
printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n");
transfer.PrintReaction();
transfer.PrintReaction(false);
if(transfer.GetRecoil().isDecay) {
decay.SetMotherDaugther(transfer.GetRecoil());
@ -118,19 +117,9 @@ void Transfer(
// msB.LoadStoppingPower(reactConfig.recoilHeavyStoppingPowerFile);
// }
//*############################################# Decay of particle-B
// Decay decay[2];
// if(reactConfig.isDecay) {
// printf("\e[32m#################################### Decay\e[0m\n");
// decay.SetMotherDaugther(reactConfig.recoilHeavyA,
// reactConfig.recoilHeavyZ,
// reactConfig.heavyDecayA,
// reactConfig.heavyDecayZ);
// }
ExcitedEnergies exList = transfer.GetRectionConfig().exList[ID];
//############################################# Load DWBAroot for thetaCM distribution
//*############################################# 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");
@ -178,7 +167,7 @@ void Transfer(
}
}
//############################################# build tree
//*############################################# 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");
@ -192,6 +181,9 @@ void Transfer(
if( distList != NULL ) distList->Write("DWBA", 1);
if( dwbaExList != NULL ) dwbaExList->Write("DWBA_ExList", 1);
TMacro idMacro;
idMacro.AddLine(Form("%d", ID));
idMacro.Write("detGeoID");
TMacro hitMeaning;
hitMeaning.AddLine("======================= meaning of Hit\n");
@ -421,7 +413,7 @@ void Transfer(
double theta = reactConfig.beamTheta;
double phi = 0.0;
//====================================================== calculate event
//*====================================================== calculate event
int count = 0;
for( int i = 0; i < numEvent; i++){
bool redoFlag = true;
@ -497,7 +489,6 @@ void Transfer(
int decayID = 0;
if( recoil.isDecay){
//decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2.); // decay to ground state
decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state
if( decayID == 1 ){
PB = decay.GetDaugther_D();
@ -661,3 +652,58 @@ void Transfer(
return;
}
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("========== Simulate Transfer reaction in HELIOS ==========\n");
printf("=================================================================\n");
if(argc == 2 || argc > 7) {
printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n");
printf(" default file name \n");
printf(" [1] reactionConfig.txt (input) reaction Setting \n");
printf(" [2] detectorGeo.txt (input) detector Setting \n");
printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n");
printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n");
printf(" [5] transfer.root (output) rootFile name for output \n");
printf(" [6] plot (input) will it plot stuffs [1/0] \n");
printf("------------------------------------------------------\n");
return 0 ;
}
std::string basicConfig = "reactionConfig.txt";
std::string heliosDetGeoFile = "detectorGeo.txt";
int ID = 0;
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root";
bool isPlot = false;
if( argc >= 2) basicConfig = argv[1];
if( argc >= 3) heliosDetGeoFile = argv[2];
if( argc >= 4) ID = atoi(argv[3]);
if( argc >= 5) ptolemyRoot = argv[4];
if( argc >= 6) saveFileName = argv[5];
if( argc >= 7) isPlot = atoi(argv[7]);
Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName);
//run Armory/Check_Simulation
if( isPlot ){
ifstream file_in;
file_in.open("../Cleopatra/Check_Simulation.C", ios::in);
if( file_in){
printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data());
TString cmd;
cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data());
system(cmd.Data());
}else{
printf("cannot find ../Cleopatra/Check_Simulation.C \n");
}
}
}

View File

@ -12,7 +12,7 @@
#include <RQ_OBJECT.h>
#include "../Cleopatra/Transfer.h"
#include "../Cleopatra/SimTransfer.C"
#include "../Cleopatra/InFileCreator.h"
#include "../Cleopatra/ExtractXSec.h"
#include "../Cleopatra/PlotTGraphTObjArray.h"
@ -58,7 +58,6 @@ private:
TGComboBox * extractFlag;
TGTextEntry * txtName ;
TGTextEntry * txtEx ;

View File

@ -1,72 +0,0 @@
/***********************************************************************
*
* 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 > 7) {
printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n");
printf(" default file name \n");
printf(" [1] reactionConfig.txt (input) reaction Setting \n");
printf(" [2] detectorGeo.txt (input) detector Setting \n");
printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n");
printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n");
printf(" [5] transfer.root (output) rootFile name for output \n");
printf(" [6] plot (input) will it plot stuffs [1/0] \n");
printf("------------------------------------------------------\n");
return 0 ;
}
string basicConfig = "reactionConfig.txt";
string heliosDetGeoFile = "detectorGeo.txt";
int ID = 0;
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root";
bool isPlot = false;
if( argc >= 2) basicConfig = argv[1];
if( argc >= 3) heliosDetGeoFile = argv[2];
if( argc >= 4) ID = atoi(argv[3]);
if( argc >= 5) ptolemyRoot = argv[4];
if( argc >= 6) saveFileName = argv[5];
if( argc >= 7) isPlot = atoi(argv[7]);
Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName);
//run Armory/Check_Simulation
if( isPlot ){
ifstream file_in;
file_in.open("../Cleopatra/Check_Simulation.C", ios::in);
if( file_in){
printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data());
TString cmd;
cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data());
system(cmd.Data());
}else{
printf("cannot find ../Cleopatra/Check_Simulation.C \n");
}
}
}

View File

@ -1,6 +1,6 @@
CC=g++
ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM Transfer SimAlpha
ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimAlpha
all: $(ALL)
@ -25,11 +25,11 @@ Cleopatra: Cleopatra.C
FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
$(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs`
Transfer: Transfer.C Transfer.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
$(CC) Transfer.C -o Transfer `root-config --cflags --glibs`
SimTransfer: SimTransfer.C ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
$(CC) SimTransfer.C -o SimTransfer `root-config --cflags --glibs`
SimAlpha: alpha.C ../Cleopatra/ClassHelios.h
$(CC) alpha.C -o SimAlpha `root-config --cflags --glibs`
SimAlpha: SimAlpha.C ../Cleopatra/ClassHelios.h
$(CC) SimAlpha.C -o SimAlpha `root-config --cflags --glibs`
clean:
/bin/rm -f $(ALL)