update Check_Simulation.C for the updated ClassDetGeo and ClassReactionConfig.h, some improvement of SimTransfer

This commit is contained in:
Ryan Tang 2024-04-03 17:15:24 -04:00
parent f921d03a39
commit e9f397cb16
7 changed files with 89 additions and 56 deletions

View File

@ -128,7 +128,8 @@
"test.C": "cpp", "test.C": "cpp",
"SimTransfer.C": "cpp", "SimTransfer.C": "cpp",
"httpaccess.C": "cpp", "httpaccess.C": "cpp",
"httpcontrol.C": "cpp" "httpcontrol.C": "cpp",
"SimTransfer2.C": "cpp"
}, },
"better-comments.multilineComments": true, "better-comments.multilineComments": true,

View File

@ -122,7 +122,7 @@ public:
std::vector<Array> array; std::vector<Array> array;
std::vector<Auxillary> aux; std::vector<Auxillary> aux;
void Print( bool printAll = false) ; void Print( bool printArray = false) ;
private: private:
@ -234,7 +234,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
} }
inline void DetGeo::Print(bool printAll){ inline void DetGeo::Print(bool printArray){
printf("#####################################################\n"); printf("#####################################################\n");
printf(" B-field : %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan"); printf(" B-field : %8.2f T, %s\n", Bfield, Bfield > 0 ? "out of plan" : "into plan");
@ -243,9 +243,9 @@ inline void DetGeo::Print(bool printAll){
printf(" z-Min : %8.2f mm\n", zMin); printf(" z-Min : %8.2f mm\n", zMin);
printf(" z-Max : %8.2f mm\n", zMax); printf(" z-Max : %8.2f mm\n", zMax);
if( printArray ) {
for( size_t i = 0; i < array.size() ; i++){ for( size_t i = 0; i < array.size() ; i++){
printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled"); printf("================================= %zu-th Detector Info (%s)\n", i, array[i].enable ? "enabled" : "disabled");
if( printAll || array[i].enable ) {
array[i].Print(); array[i].Print();
aux[i].Print(); aux[i].Print();
} }

View File

@ -109,7 +109,8 @@ public:
void SetReactionSimple(int beamA, int beamZ, void SetReactionSimple(int beamA, int beamZ,
int targetA, int targetZ, int targetA, int targetZ,
int recoilA, int recoilZ, float beamEnergy_AMeV, unsigned short ID); int recoilA, int recoilZ,
float beamEnergy_AMeV);
bool LoadReactionConfig(TString fileName); bool LoadReactionConfig(TString fileName);
bool LoadReactionConfig(TMacro * macro); bool LoadReactionConfig(TMacro * macro);
@ -122,17 +123,20 @@ private:
inline void ReactionConfig::SetReactionSimple(int beamA, int beamZ, inline void ReactionConfig::SetReactionSimple(int beamA, int beamZ,
int targetA, int targetZ, int targetA, int targetZ,
int recoilA, int recoilZ, float beamEnergy_AMeV, unsigned short ID){ int recoilA, int recoilZ,
float beamEnergy_AMeV){
this->beamA = beamA; this->beamA = beamA;
this->beamZ = beamZ; this->beamZ = beamZ;
this->targetA = targetA; this->targetA = targetA;
this->targetZ = targetZ; this->targetZ = targetZ;
this->recoil[ID].lightA = recoilA; this->recoil.push_back(Recoil());
this->recoil[ID].lightZ = recoilZ;
recoil[ID].heavyA = this->beamA + this->targetA - recoil[ID].lightA; this->recoil.back().lightA = recoilA;
recoil[ID].heavyZ = this->beamZ + this->targetZ - recoil[ID].lightZ; this->recoil.back().lightZ = recoilZ;
recoil.back().heavyA = this->beamA + this->targetA - recoil.back().lightA;
recoil.back().heavyZ = this->beamZ + this->targetZ - recoil.back().lightZ;
beamEnergy = beamEnergy_AMeV; beamEnergy = beamEnergy_AMeV;
beamEnergySigma = 0; beamEnergySigma = 0;
@ -260,8 +264,8 @@ inline void ReactionConfig::Print(int ID, bool withEx) const{
printf("================================= Number of recoil reactions : %zu\n", recoil.size()); printf("================================= Number of recoil reactions : %zu\n", recoil.size());
for( size_t i = 0; i < recoil.size(); i ++ ){ for( size_t i = 0; i < recoil.size(); i ++ ){
printf("------------------------------------------ Recoil-%zu\n", i);
if( ID == i || ID < 0 ){ if( ID == i || ID < 0 ){
printf("------------------------------------------ Recoil-%zu\n", i);
recoil[i].Print(); recoil[i].Print();
if( withEx ) exList[i].Print(); if( withEx ) exList[i].Print();
} }

View File

@ -73,16 +73,16 @@ void Check_Simulation(TString filename = "transfer.root",
TString gate = ExtractString(startLineNum+1, config); TString gate = ExtractString(startLineNum+1, config);
double elumRange = ExtractNumber(startLineNum+2, config); double elumRange = ExtractNumber(startLineNum+2, config);
vector<double> thetaCMRange = doubleConvertor( StringToVector( ExtractString(startLineNum+3,config) )); std::vector<double> thetaCMRange = doubleConvertor( StringToVector( ExtractString(startLineNum+3,config) ));
bool shownKELines = (ExtractString(startLineNum+4, config).Remove(4) == "true" ? true : false); bool shownKELines = (ExtractString(startLineNum+4, config).Remove(4) == "true" ? true : false);
bool isOverRideEx = (ExtractString(startLineNum+5, 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 ))); std::vector<double> oExRange = doubleConvertor( StringToVector ( ExtractString(startLineNum+6, config )));
printf("%s \n", gate.Data()); printf("%s \n", gate.Data());
///==== config Canvas ///==== config Canvas
vector<TString> plotConfig = StringToVector( ExtractString(startLineNum, config)); std::vector<TString> plotConfig = StringToVector( ExtractString(startLineNum, config));
vector<plotID> canvas; std::vector<plotID> canvas;
int colCount = 0; int colCount = 0;
int colCount_new = 0; int colCount_new = 0;
int rowCount = 1; int rowCount = 1;
@ -100,7 +100,7 @@ void Check_Simulation(TString filename = "transfer.root",
if( colCount == 0 ) colCount = colCount_new; 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}; std::vector<int> Div = {colCount, rowCount};
TFile * file = new TFile(filename, "read"); TFile * file = new TFile(filename, "read");
TTree * tree = (TTree*) file->Get("tree"); TTree * tree = (TTree*) file->Get("tree");
@ -140,8 +140,8 @@ void Check_Simulation(TString filename = "transfer.root",
DetGeo detGeo(detGeoTxt); DetGeo detGeo(detGeoTxt);
Array array = detGeo.array[detGeoID]; Array array = detGeo.array[detGeoID];
detGeo.PrintWithoutArray(); detGeo.Print();
array.PrintArray(); array.Print();
printf("=================================\n"); printf("=================================\n");
@ -179,11 +179,12 @@ void Check_Simulation(TString filename = "transfer.root",
}else{ }else{
numEx = exListTxt->GetListOfLines()->GetSize()-1; numEx = exListTxt->GetListOfLines()->GetSize()-1;
for( int i = 0 ; i < numEx ; i++){ for( int i = 1 ; i <= numEx ; i++){
double ex = atof(exListTxt->GetListOfLines()->At(i)->GetName()); std::vector<std::string> tempStr = AnalysisLib::SplitStr(exListTxt->GetListOfLines()->At(i)->GetName(), " ");
double ex = atof(tempStr[0].c_str());
if( ex < ExRange[0] ) ExRange[0] = ex; if( ex < ExRange[0] ) ExRange[0] = ex;
if( ex > ExRange[1] ) ExRange[1] = ex; if( ex > ExRange[1] ) ExRange[1] = ex;
exList.Add(ex, 0, 0, 0); exList.Add(ex, atof(tempStr[1].c_str()), 1.0, 0.00);
} }
} }
@ -191,8 +192,8 @@ void Check_Simulation(TString filename = "transfer.root",
exList.Print(); exList.Print();
double dExRange = ExRange[1] - ExRange[0]; double dExRange = ExRange[1] - ExRange[0];
ExRange[0] = ExRange[0] - dExRange * 0.1; ExRange[0] = ExRange[0] - 0.3 - dExRange * 0.1;
ExRange[1] = ExRange[1] + dExRange * 0.1; ExRange[1] = ExRange[1] + 0.3 + dExRange * 0.1;
printf("Number of Ex states = %d \n", numEx); printf("Number of Ex states = %d \n", numEx);
@ -200,7 +201,13 @@ void Check_Simulation(TString filename = "transfer.root",
//eRange by zRange and exList //eRange by zRange and exList
TransferReaction transfer; TransferReaction transfer;
transfer.SetReactionSimple( reactionConfig.beamA, reactionConfig.beamZ, reactionConfig.targetA, reactionConfig.targetZ, recoil.lightA, recoil.lightZ, reactionConfig.beamEnergy); transfer.SetReactionSimple( reactionConfig.beamA,
reactionConfig.beamZ,
reactionConfig.targetA,
reactionConfig.targetZ,
recoil.lightA,
recoil.lightZ,
reactionConfig.beamEnergy);
double QQ = transfer.GetCMTotalEnergy(); double QQ = transfer.GetCMTotalEnergy();
double gamm = transfer.GetReactionGamma(); double gamm = transfer.GetReactionGamma();
@ -208,7 +215,6 @@ void Check_Simulation(TString filename = "transfer.root",
double slope = transfer.GetEZSlope( detGeo.Bfield); double slope = transfer.GetEZSlope( detGeo.Bfield);
double eRange[2] = {0, 10}; double eRange[2] = {0, 10};
// double intercept = QQ/gamm - mass;
eRange[1] = zRange[2] * slope; eRange[1] = zRange[2] * slope;
// printf("intercept of 0 MeV : %f MeV \n", intercept); // printf("intercept of 0 MeV : %f MeV \n", intercept);
@ -272,14 +278,15 @@ void Check_Simulation(TString filename = "transfer.root",
} }
if( pID == pRecoilXY ){ if( pID == pRecoilXY ){
TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius, TH2F * hRecoilXY = new TH2F("hRecoilXY", Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[detGeoID].detPos ),
400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius); 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius,
400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius);
tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz"); tree->Draw("yRecoil:xRecoil>>hRecoilXY", gate, "colz");
TArc * detArc1 = new TArc(0,0, detGeo.recoilOuterRadius); TArc * detArc1 = new TArc(0,0, detGeo.aux[detGeoID].outerRadius);
detArc1->SetLineColor(kBlue-8); detArc1->SetLineColor(kBlue-8);
detArc1->SetFillStyle(0); detArc1->SetFillStyle(0);
detArc1->Draw("same"); detArc1->Draw("same");
TArc * detArc2 = new TArc(0,0, detGeo.recoilInnerRadius); TArc * detArc2 = new TArc(0,0, detGeo.aux[detGeoID].innerRadius);
detArc2->SetLineColor(kBlue-8); detArc2->SetLineColor(kBlue-8);
detArc2->SetFillStyle(0); detArc2->SetFillStyle(0);
detArc2->Draw("same"); detArc2->Draw("same");
@ -293,25 +300,27 @@ void Check_Simulation(TString filename = "transfer.root",
} }
if( pID == pRecoilXY1 ){ if( pID == pRecoilXY1 ){
TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos1 ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius, TH2F * hRecoilXY1 = new TH2F("hRecoilXY1", Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[detGeoID].detPos1 ),
400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius); 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius,
400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius);
tree->Draw("yRecoil1:xRecoil1>>hRecoilXY1", gate, "colz"); tree->Draw("yRecoil1:xRecoil1>>hRecoilXY1", gate, "colz");
} }
if( pID == pRecoilXY2 ){ if( pID == pRecoilXY2 ){
TH2F * hRecoilXY2 = new TH2F("hRecoilXY2", Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.recoilPos2 ), 400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius, TH2F * hRecoilXY2 = new TH2F("hRecoilXY2", Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[detGeoID].detPos2 ),
400, -detGeo.recoilOuterRadius, detGeo.recoilOuterRadius); 400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius,
400, -detGeo.aux[detGeoID].outerRadius, detGeo.aux[detGeoID].outerRadius);
tree->Draw("yRecoil2:xRecoil2>>hRecoilXY2", gate, "colz"); tree->Draw("yRecoil2:xRecoil2>>hRecoilXY2", gate, "colz");
} }
if( pID == pRecoilRZ ){ if( pID == pRecoilRZ ){
TH2F * hRecoilRZ = new TH2F("hRecoilRZ", "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", zRange[0], zRange[1], zRange[2], 400,0, detGeo.recoilOuterRadius); TH2F * hRecoilRZ = new TH2F("hRecoilRZ", "RecoilR - Z [gated]; z [mm]; RecoilR [mm]", zRange[0], zRange[1], zRange[2], 400,0, detGeo.aux[detGeoID].outerRadius);
tree->Draw("rhoRecoil:z>>hRecoilRZ", gate, "colz"); tree->Draw("rhoRecoil:z>>hRecoilRZ", gate, "colz");
} }
if( pID == pRecoilRTR ){ if( pID == pRecoilRTR ){
FindRange("TB", gate, tree, recoilERange); 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, detGeo.recoilOuterRadius); TH2F * hRecoilRTR = new TH2F("hRecoilRTR", "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", 500, recoilERange[0], recoilERange[1], 500, 0, detGeo.aux[detGeoID].outerRadius);
tree->Draw("rhoRecoil:TB>>hRecoilRTR", gate, "colz"); tree->Draw("rhoRecoil:TB>>hRecoilRTR", gate, "colz");
} }
@ -391,7 +400,7 @@ void Check_Simulation(TString filename = "transfer.root",
} }
if( pID == pRecoilRThetaCM ){ if( pID == pRecoilRThetaCM ){
TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, detGeo.recoilOuterRadius); TH2F * hRecoilRThetaCM = new TH2F("hRecoilRThetaCM", "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 400, 0, 60, 400,0, detGeo.aux[detGeoID].outerRadius);
tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz"); tree->Draw("rhoRecoil:thetaCM>>hRecoilRThetaCM", gate, "colz");
} }
@ -430,7 +439,7 @@ void Check_Simulation(TString filename = "transfer.root",
} }
if( pID == pElum1XY ){ if( pID == pElum1XY ){
TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.elumPos1), 400, -elumRange, elumRange, 400, -elumRange, elumRange); TH2F * hElum1XY = new TH2F("hElum1XY", Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.aux[detGeoID].elumPos1), 400, -elumRange, elumRange, 400, -elumRange, elumRange);
tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz"); tree->Draw("yElum1:xElum1>>hElum1XY", gate, "colz");
double count = hElum1XY->GetEntries(); double count = hElum1XY->GetEntries();

View File

@ -158,11 +158,13 @@ TransferReaction::TransferReaction(std::string configFile, unsigned short ID){
TransferReaction::TransferReaction(int beamA, int beamZ, TransferReaction::TransferReaction(int beamA, int beamZ,
int targetA, int targetZ, int targetA, int targetZ,
int recoilA, int recoilZ, float beamEnergy_AMeV){ int recoilA, int recoilZ,
float beamEnergy_AMeV){
Inititization(); Inititization();
SetReactionSimple(beamA, beamZ, SetReactionSimple(beamA, beamZ,
targetA, targetZ, targetA, targetZ,
recoilA, recoilZ, beamEnergy_AMeV); recoilA, recoilZ,
beamEnergy_AMeV);
} }
void TransferReaction::Inititization(){ void TransferReaction::Inititization(){
@ -244,11 +246,12 @@ void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double p
void TransferReaction::SetReactionSimple(int beamA, int beamZ, void TransferReaction::SetReactionSimple(int beamA, int beamZ,
int targetA, int targetZ, int targetA, int targetZ,
int recoilA, int recoilZ, float beamEnergy_AMeV){ int recoilA, int recoilZ,
float beamEnergy_AMeV){
config.SetReactionSimple(beamA, beamZ, config.SetReactionSimple(beamA, beamZ,
targetA, targetZ, targetA, targetZ,
recoilA, recoilZ, beamEnergy_AMeV, 0); recoilA, recoilZ, beamEnergy_AMeV);
recoil = config.recoil[0]; recoil = config.recoil[0];
@ -281,6 +284,11 @@ void TransferReaction::SetReactionFromFile(std::string configFile, unsigned shor
SetExA(config.beamEx); SetExA(config.beamEx);
if( ID > config.recoil.size() ){
printf("Reaction Config only has %zu recoil settings. input is %u. Abort.\n", config.recoil.size(), ID);
return;
}
recoil = config.recoil[ID]; recoil = config.recoil[ID];
exList = config.exList[ID]; exList = config.exList[ID];

View File

@ -126,21 +126,30 @@ void Transfer(
TFile * distFile = new TFile(ptolemyRoot, "read"); TFile * distFile = new TFile(ptolemyRoot, "read");
TObjArray * distList = nullptr; TObjArray * distList = nullptr;
TMacro * dwbaExList = nullptr; TMacro * dwbaExList = nullptr;
TMacro * dwbaReactList = nullptr;
TMacro dwbaExList_Used;
if( distFile->IsOpen() ) { if( distFile->IsOpen() ) {
printf("--------- Found DWBA thetaCM distributions. Use the ExList from DWBA.\n"); printf("--------- Found DWBA thetaCM distributions. Use the ExList from DWBA.\n");
distList = (TObjArray *) distFile->FindObjectAny("thetaCM_TF1"); // the function List distList = (TObjArray *) distFile->FindObjectAny("thetaCM_TF1"); // the function List
dwbaExList = (TMacro *) distFile->FindObjectAny("ExList");
dwbaExList_Used.AddLine(dwbaExList->GetListOfLines()->At(0)->GetName());
dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList");
exList.Clear(); exList.Clear();
dwbaExList = (TMacro *) distFile->FindObjectAny("ExList");
int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ; int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ;
for(int i = 1; i <= numEx ; i++){ for(int i = 1; i <= numEx ; i++){
std::string reactionName = dwbaReactList->GetListOfLines()->At(i-1)->GetName();
if( reactionName.find( transfer.GetReactionName().Data() ) != std::string::npos) {
std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName();
dwbaExList_Used.AddLine(temp.c_str());
if( temp[0] == '/' ) continue; if( temp[0] == '/' ) continue;
std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " "); std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " ");
exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00); exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
} }
}
}else{ }else{
printf("------- no DWBA input. Use the ExList from %s\n", basicConfig.c_str()); printf("------- no DWBA input. Use the ExList from %s\n", basicConfig.c_str());
@ -180,7 +189,9 @@ void Transfer(
detGeoTxt.Write("detGeo"); detGeoTxt.Write("detGeo");
if( distList != NULL ) distList->Write("DWBA", 1); if( distList != NULL ) distList->Write("DWBA", 1);
if( dwbaExList != NULL ) dwbaExList->Write("DWBA_ExList", 1); if( dwbaExList != NULL ) {
dwbaExList_Used.Write("DWBA_ExList", 1);
}
TMacro idMacro; TMacro idMacro;
idMacro.AddLine(Form("%d", ID)); idMacro.AddLine(Form("%d", ID));

View File

@ -28,7 +28,7 @@ FindThetaCM: FindThetaCM.C FindThetaCM.h ClassTransfer.h ClassHelios.h ClassIsot
SimTransfer: SimTransfer.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h SimTransfer: SimTransfer.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h
$(CC) SimTransfer.C -o SimTransfer `root-config --cflags --glibs` $(CC) SimTransfer.C -o SimTransfer `root-config --cflags --glibs`
SimTransfer2: SimTransfer2.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h SimTransfer2: SimTransfer2.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h ../Armory/ClassReactionConfig.h ../Armory/ClassDetGeo.h
$(CC) SimTransfer2.C -o SimTransfer2 `root-config --cflags --glibs` $(CC) SimTransfer2.C -o SimTransfer2 `root-config --cflags --glibs`
SimAlpha: SimAlpha.C ClassHelios.h SimAlpha: SimAlpha.C ClassHelios.h