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

replace SimTransfer2 to SimTransfer. Already able to use N-arrays and reactions

Browse Source
This commit is contained in:
Ryan Tang 2024-07-04 13:06:26 -04:00
parent 83177e57c7
commit e944682888
13 changed files with 630 additions and 605 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -27,7 +27,7 @@ Cleopatra/PlotSimulation
Cleopatra/PlotTGraphTObjArray Cleopatra/PlotTGraphTObjArray
Cleopatra/SimAlpha Cleopatra/SimAlpha
Cleopatra/SimTransfer Cleopatra/SimTransfer
Cleopatra/SimTransfer2 Cleopatra/SimTransfer_single
Cleopatra/Cleopatra Cleopatra/Cleopatra
__pycache__ __pycache__

3
.vscode/settings.json vendored
View File

@ -130,7 +130,8 @@
"httpaccess.C": "cpp", "httpaccess.C": "cpp",
"httpcontrol.C": "cpp", "httpcontrol.C": "cpp",
"SimTransfer2.C": "cpp", "SimTransfer2.C": "cpp",
"haha.C": "cpp" "haha.C": "cpp",
"SimTransfer_2.C": "cpp"
}, },
"better-comments.multilineComments": true, "better-comments.multilineComments": true,

2
Armory/GeneralSort.C
View File

@ -20,8 +20,6 @@ Bool_t GeneralSort::Process(Long64_t entry){
if( entry < 1 ) printf("============================== start processing data\n"); if( entry < 1 ) printf("============================== start processing data\n");
if( entry > 5000) return true;
///initialization ///initialization
for( int i = 0; i < mapping::nDetType; i++){ for( int i = 0; i < mapping::nDetType; i++){
if( mapping::detNum[i] == 0 ) continue; if( mapping::detNum[i] == 0 ) continue;

2
Armory/Process_Run
View File

@ -34,7 +34,7 @@ fi;
RUN=$1 RUN=$1
runNum=$1 runNum=$1
EventBld=2 EventBld=1
nWorker=1 nWorker=1
TraceMethod=0 TraceMethod=0
isMonitor=1 isMonitor=1

6
Cleopatra/ClassHelios.h
View File

@ -77,7 +77,7 @@ class HELIOS{
public: public:
HELIOS(); HELIOS();
HELIOS(std::string filename, unsigned short ID); HELIOS(std::string detGeoFile, unsigned short ID);
~HELIOS(); ~HELIOS();
void SetCoincidentWithRecoil(bool TorF){ this->isCoincidentWithRecoil = TorF;} void SetCoincidentWithRecoil(bool TorF){ this->isCoincidentWithRecoil = TorF;}
@ -176,9 +176,9 @@ HELIOS::HELIOS(){
Clear(); Clear();
} }
HELIOS::HELIOS(std::string filename, unsigned short ID){ HELIOS::HELIOS(std::string detGeoFile, unsigned short ID){
Clear(); Clear();
SetDetectorGeometry(filename, ID); SetDetectorGeometry(detGeoFile, ID);
} }
HELIOS::~HELIOS(){ HELIOS::~HELIOS(){

0
Cleopatra/SimHelper.py → Cleopatra/PySimHelper.py
View File

581
Cleopatra/SimTransfer.C
View File

@ -19,9 +19,7 @@
#include "ClassTransfer.h" #include "ClassTransfer.h"
#include "ClassHelios.h" #include "ClassHelios.h"
double exDistFunc(Double_t *x, Double_t * par){
return par[(int) x[0]];
}
void PrintEZPlotPara(TransferReaction tran, HELIOS helios){ void PrintEZPlotPara(TransferReaction tran, HELIOS helios){
@ -33,121 +31,100 @@ void PrintEZPlotPara(TransferReaction tran, HELIOS helios){
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
double slope = tran.GetEZSlope(helios.GetBField()); /// MeV/mm double slope = tran.GetEZSlope(helios.GetBField()); /// MeV/mm
printf(" e-z slope : %f MeV/mm\n", slope); printf(" e-z slope : %f MeV/mm\n", slope);
double intercept = q/gamma - mb; // MeV // double intercept = q/gamma - mb; // MeV
printf(" e-z intercept (ground state) : %f MeV\n", intercept); // printf(" e-z intercept (ground state) : %f MeV\n", intercept);
} }
void Transfer( void Transfer(
std::string basicConfig = "reactionConfig.txt", std::string basicConfig = "reactionConfig.txt",
std::string heliosDetGeoFile = "detectorGeo.txt", std::string detGeoFile = "detectorGeo.txt",
unsigned short ID = 0, // this is the ID for the array TString ptolemyRoot = "DWBA.root",
TString ptolemyRoot = "DWBA.root", TString saveFileName = "transfer.root"){
TString saveFileName = "transfer.root"){
//*############################################# Set Reaction //*############################################# Set Reaction
TransferReaction transfer;
HELIOS helios;
Decay decay;
std::vector<double> kbCM; /// momentum of b in CM frame // std::vector<double> kbCM; /// momentum of b in CM frame
TF1 * exDist = nullptr; // TF1 * exDistribution = nullptr;
transfer.SetReactionFromFile(basicConfig, ID); DetGeo detGeoConfig;
helios.SetDetectorGeometry(heliosDetGeoFile, ID); ReactionConfig reactionConfig;
detGeoConfig.LoadDetectorGeo(detGeoFile, false);
reactionConfig.LoadReactionConfig(basicConfig);
const unsigned short numDetGeo = detGeoConfig.array.size();
const unsigned short numReact = reactionConfig.recoil.size();
if( numDetGeo != numReact ){
printf("\e[31m !!!!!! number of array is not equal to number of reaction.!!! \e[0m\n");
printf("Abort\n");
return;
}
unsigned short numTransfer = 0;
for( int i = 0; i < std::min(numDetGeo, numReact); i++){
if( detGeoConfig.array[i].enable ) numTransfer ++;
}
TransferReaction * transfer = new TransferReaction[numTransfer];
Decay * decay = new Decay[numTransfer];
HELIOS * helios = new HELIOS[numTransfer];
int count = 0;
for( unsigned short i = 0 ; i < numDetGeo; i++){
if( detGeoConfig.array[i].enable ){
transfer[count].SetReactionFromFile(basicConfig, i);
if(transfer[count].GetRecoil().isDecay) {
decay[count].SetMotherDaugther(transfer[count].GetRecoil());
}
helios[count].SetDetectorGeometry(detGeoFile, i);
count ++;
}
}
printf("*****************************************************************\n");
printf("*\e[1m\e[33m %27s \e[0m*\n", transfer.GetReactionName().Data());
printf("*****************************************************************\n");
printf("----- loading reaction setting from %s. \n", basicConfig.c_str()); printf("----- loading reaction setting from %s. \n", basicConfig.c_str());
printf("----- loading geometry setting from %s. \n", heliosDetGeoFile.c_str()); printf("----- loading geometry setting from %s. \n", detGeoFile.c_str());
printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n"); printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n");
transfer.PrintReaction(false);
if(transfer.GetRecoil().isDecay) {
decay.SetMotherDaugther(transfer.GetRecoil());
}
helios.PrintGeometry();
PrintEZPlotPara(transfer, helios);
DetGeo detGeo = helios.GetDetectorGeometry();
Array array = helios.GetArrayGeometry();
Auxillary aux = helios.GetAuxGeometry();
ReactionConfig reactConfig = transfer.GetRectionConfig();
Recoil recoil = transfer.GetRecoil();
//*############################################# 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", transfer.GetMass_b(), "mass_b");
// fprintf(keyParaOut, "%-15d //%s\n", reactConfig.recoilLightZ, "charge_b");
// fprintf(keyParaOut, "%-15.8f //%s\n", transfer.GetReactionBeta(), "betaCM");
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetCMTotalEnergy(), "Ecm");
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.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 = reactConfig.isTargetScattering;
// float density = reactConfig.targetDensity;
// float targetThickness = reactConfig.targetThickness;
// if(isTargetScattering) printf("\e[32m#################################### Target Scattering\e[0m\n");
// TargetScattering msA;
// TargetScattering msB;
// TargetScattering msb;
// if(reactConfig.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( reactConfig.isTargetScattering ){
// msA.LoadStoppingPower(reactConfig.beamStoppingPowerFile);
// msb.LoadStoppingPower(reactConfig.recoilLightStoppingPowerFile);
// msB.LoadStoppingPower(reactConfig.recoilHeavyStoppingPowerFile);
// }
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"); TFile * distFile = new TFile(ptolemyRoot, "read");
TObjArray * distList = nullptr; TObjArray * distList = nullptr;
TMacro * dwbaExList = nullptr; TMacro * dwbaExList = nullptr;
TMacro * dwbaReactList = nullptr; TMacro * dwbaReactList = nullptr;
TMacro dwbaExList_Used; TMacro dwbaExList_Used;
TMacro dwbaReactList_Used;
if( distFile->IsOpen() ) { if( distFile->IsOpen() ) {
printf("\e[32m#################################### Load DWBA input : %s \e[0m\n", ptolemyRoot.Data());
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 = (TMacro *) distFile->FindObjectAny("ExList");
dwbaExList_Used.AddLine(dwbaExList->GetListOfLines()->At(0)->GetName()); dwbaExList_Used.AddLine(dwbaExList->GetListOfLines()->At(0)->GetName());
dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList"); dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList");
exList.Clear();
int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ; int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ;
for(int i = 1; i <= numEx ; i++){
for( int i = 0; i < numTransfer; i++){ transfer[i].GetExList()->Clear(); }
for( int i = 1; i <= numEx ; i++){
//Check DWBA reaction is same as transfer setting
std::string reactionName = dwbaReactList->GetListOfLines()->At(i-1)->GetName(); 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(); for( int j = 0; j < numTransfer; j++){
dwbaExList_Used.AddLine(temp.c_str()); if( reactionName.find( transfer[j].GetReactionName().Data() ) != std::string::npos) {
if( temp[0] == '/' ) continue; std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName();
std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " "); dwbaReactList_Used.AddLine((reactionName + " | " + std::to_string(j)).c_str());
exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00); dwbaExList_Used.AddLine(temp.c_str());
if( temp[0] == '/' ) continue;
std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " ");
transfer[j].GetExList()->Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
}
} }
} }
@ -155,25 +132,20 @@ void Transfer(
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());
} }
std::vector<bool> listOfTransfer(numTransfer, false);
printf("------------------------------ Heavy Recoil excitation\n"); for( int i = 0; i < numTransfer; i++){
printf("Energy[MeV] Rel.Xsec SF sigma\n"); if( transfer[i].GetExList()->ExList.size() > 0 ){
int numEx = exList.ExList.size(); listOfTransfer[i] = true;
for( int j = 0; j < numEx; j++){ transfer[i].PrintReaction(false);
double ex = exList.ExList[j].Ex; transfer[i].GetExList()->Print();
kbCM.push_back(transfer.CalkCM(ex)); helios[i].PrintGeometry();
int decayID = decay.CalDecay(TLorentzVector (0,0,0,0), ex, 0); transfer[i].CreateExDistribution();
exList.ExList[j].Print(decayID == 1 ? "-->Decay" : "\n"); // PrintEZPlotPara(transfer[i], helios[i]);
} }else{
printf(" Reaction : %s has no excited energy. Skipped. \n", transfer[i].GetReactionName().Data());
//---- create Ex-distribution
if( exList.ExList.size() > 1 ) {
printf("---- creating Ex-distribution \n");
exDist = new TF1("exDist", exDistFunc, 0, numEx, numEx);
for(int q = 0; q < numEx; q++){
exDist->SetParameter(q, exList.ExList[q].xsec*exList.ExList[q].SF);
} }
} }
@ -183,23 +155,22 @@ void Transfer(
TTree * tree = new TTree("tree", "tree"); TTree * tree = new TTree("tree", "tree");
TMacro config(basicConfig.c_str()); TMacro config(basicConfig.c_str());
TMacro detGeoTxt(heliosDetGeoFile.c_str()); TMacro detGeoTxt(detGeoFile.c_str());
config.SetName(transfer.GetReactionName_Latex().Data()); config.SetName("ReactionConfig");
config.Write("reactionConfig"); config.Write("reactionConfig");
detGeoTxt.Write("detGeo"); detGeoTxt.Write("detGeo");
if( distList != NULL ) distList->Write("DWBA", 1); if( distList != NULL ) distList->Write("DWBA", 1);
if( dwbaExList != NULL ) dwbaExList_Used.Write("DWBA_ExList", 1); if( dwbaExList != NULL ) {
dwbaExList_Used.Write("DWBA_ExList", 1);
dwbaReactList_Used.Write("DWBA_ReactionList", 1);
}
TMacro idMacro;
idMacro.AddLine(Form("%d", ID));
idMacro.Write("detGeoID");
TMacro hitMeaning; TMacro hitMeaning;
hitMeaning.AddLine("======================= meaning of Hit\n"); hitMeaning.AddLine("======================= meaning of Hit\n");
for( int code = -15 ; code <= 1; code ++ ){ for( int code = -15 ; code <= 1; code ++ ){
hitMeaning.AddLine( Form( "%4d = %s", code, helios.AcceptanceCodeToMsg(code).Data() )); hitMeaning.AddLine( Form( "%4d = %s", code, helios[0].AcceptanceCodeToMsg(code).Data() ));
} }
hitMeaning.AddLine(" other = unknown\n"); hitMeaning.AddLine(" other = unknown\n");
hitMeaning.AddLine("===========================================\n"); hitMeaning.AddLine("===========================================\n");
@ -208,6 +179,9 @@ void Transfer(
int hit; /// the output of Helios.CalHit int hit; /// the output of Helios.CalHit
tree->Branch("hit", &hit, "hit/I"); tree->Branch("hit", &hit, "hit/I");
int rID; /// reaction ID
tree->Branch("rID", &rID, "reactionID/I");
double thetab, phib, Tb; double thetab, phib, Tb;
double thetaB, phiB, TB; double thetaB, phiB, TB;
tree->Branch("thetab", &thetab, "thetab/D"); tree->Branch("thetab", &thetab, "thetab/D");
@ -259,7 +233,14 @@ void Transfer(
double decayTheta; /// the change of thetaB due to decay double decayTheta; /// the change of thetaB due to decay
double xRecoil_d, yRecoil_d, rhoRecoil_d, Td; double xRecoil_d, yRecoil_d, rhoRecoil_d, Td;
if( recoil.isDecay ) {
bool isAnyDecay = false;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyDecay |= transfer[i].GetRecoil().isDecay;
}
if( isAnyDecay ) {
tree->Branch("decayTheta", &decayTheta, "decayTheta/D"); tree->Branch("decayTheta", &decayTheta, "decayTheta/D");
tree->Branch("xRecoil_d", &xRecoil_d, "xRecoil_d/D"); tree->Branch("xRecoil_d", &xRecoil_d, "xRecoil_d/D");
tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D"); tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D");
@ -280,14 +261,24 @@ void Transfer(
///in case need ELUM ///in case need ELUM
double xElum1, yElum1, rhoElum1; double xElum1, yElum1, rhoElum1;
if( aux.elumPos1 != 0 ) { bool isAnyElum1 = false;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyElum1 |= (helios[i].GetAuxGeometry().elumPos1 != 0);
}
if( isAnyElum1 ) {
tree->Branch("xElum1", &xElum1, "xElum1/D"); tree->Branch("xElum1", &xElum1, "xElum1/D");
tree->Branch("yElum1", &yElum1, "yElum1/D"); tree->Branch("yElum1", &yElum1, "yElum1/D");
tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D"); tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D");
} }
double xElum2, yElum2, rhoElum2; double xElum2, yElum2, rhoElum2;
if( aux.elumPos2 != 0 ) { bool isAnyElum2 = false;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyElum2 |= (helios[i].GetAuxGeometry().elumPos2 != 0);
}
if( isAnyElum2 ) {
tree->Branch("xElum2", &xElum2, "xElum2/D"); tree->Branch("xElum2", &xElum2, "xElum2/D");
tree->Branch("yElum2", &yElum2, "yElum2/D"); tree->Branch("yElum2", &yElum2, "yElum2/D");
tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D"); tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D");
@ -295,31 +286,52 @@ void Transfer(
///in case need other recoil detector. ///in case need other recoil detector.
double xRecoil1, yRecoil1, rhoRecoil1; double xRecoil1, yRecoil1, rhoRecoil1;
if( aux.detPos1 != 0 ){ bool isAnyRecoil1 = false;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyRecoil1 |= (helios[i].GetAuxGeometry().detPos1 != 0);
}
if( isAnyRecoil1 != 0 ){
tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D"); tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D");
tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D"); tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D");
tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D"); tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D");
} }
double xRecoil2, yRecoil2, rhoRecoil2; double xRecoil2, yRecoil2, rhoRecoil2;
if( aux.detPos2 != 0 ){ bool isAnyRecoil2 = false;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyRecoil2 |= (helios[i].GetAuxGeometry().detPos2 != 0);
}
if( isAnyRecoil2 != 0 ){
tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D"); tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D");
tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D"); tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D");
tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D"); tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D");
} }
//======= list of reaction used.
TMacro listOfReaction;
for( int i = 0; i < numTransfer ; i++){
if( !listOfTransfer[i] ) continue;
listOfReaction.AddLine(Form("%2d | %s", i, transfer[i].GetReactionName_Latex().Data()));
}
listOfReaction.Write("ListOfReactions");
//======= function for e-z plot for ideal case //======= function for e-z plot for ideal case
printf("++++ generate functions\n"); printf("++++ generate functions\n");
TObjArray * gList = new TObjArray(); TObjArray * gList = new TObjArray();
gList->SetName("Constant thetaCM lines"); gList->SetName("Constant thetaCM = 0 lines");
const int gxSize = 50; const int gxSize = numTransfer;
TF1 ** gx = new TF1*[gxSize]; TF1 ** gx = new TF1*[gxSize];
TString name; TString name;
double mb = transfer.GetMass_b();
double betaRect = transfer.GetReactionBeta();
double gamma = transfer.GetReactionGamma();
double slope = transfer.GetEZSlope(helios.GetBField()); /// MeV/mm
for( int i = 0; i < gxSize; i++){ for( int i = 0; i < gxSize; i++){
double mb = transfer[i].GetMass_b();
double betaRect = transfer[i].GetReactionBeta();
double gamma = transfer[i].GetReactionGamma();
double slope = transfer[i].GetEZSlope(helios[0].GetBField()); /// MeV/mm
name.Form("g%d", i); name.Form("g%d", i);
gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000); gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000);
double thetacm = i * TMath::DegToRad(); double thetacm = i * TMath::DegToRad();
@ -338,67 +350,96 @@ void Transfer(
gList->Write("EZ_thetaCM", TObject::kSingleKey); gList->Write("EZ_thetaCM", TObject::kSingleKey);
printf(" %d constant thetaCM functions\n", gxSize); printf(" %d constant thetaCM functions\n", gxSize);
for( int i = 0; i < gxSize; i++){
delete gx[i];
}
delete [] gx;
delete gList;
//--- cal modified f //--- cal modified f
int numEx = 0;
for( int i = 0; i < numTransfer; i++){
if( !listOfTransfer[i] ) continue;
numEx += transfer[i].GetExList()->ExList.size();
}
TObjArray * fxList = new TObjArray(); TObjArray * fxList = new TObjArray();
TGraph ** fx = new TGraph*[numEx]; TGraph ** fx = new TGraph*[numEx];
std::vector<double> px, py; std::vector<double> px, py;
int countfx = 0; int countfx = 0;
for( int j = 0 ; j < numEx; j++){
double a = helios.GetDetRadius();
double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[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 / kbCM[j] * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
double pyTemp = gamma * q - mb - gamma * betaRect * kbCM[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]); count = 0;
name.Form("fx%d", j); for( int i = 0; i < numTransfer; i++ ){
fx[j]->SetName(name); if( !listOfTransfer[i] ) continue;
fx[j]->SetLineColor(4); double mb = transfer[i].GetMass_b();
fxList->Add(fx[j]); double betaRect = transfer[i].GetReactionBeta();
printf(","); double gamma = transfer[i].GetReactionGamma();
double slope = transfer[i].GetEZSlope(helios[0].GetBField()); /// MeV/mm
for( size_t j = 0 ; j < transfer[i].GetExList()->ExList.size(); j++){
double Ex = transfer[i].GetExList()->ExList[j].Ex;
double kbCM = transfer[i].CalkCM(Ex);
double a = helios[i].GetDetRadius();
double q = TMath::Sqrt(mb*mb + kbCM * kbCM );
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 / kbCM * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
double pyTemp = gamma * q - mb - gamma * betaRect * kbCM * TMath::Cos(thetacm);
if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
px.push_back(pxTemp);
py.push_back(pyTemp);
countfx ++;
}
fx[count] = new TGraph(countfx, &px[0], &py[0]);
name.Form("fx%d_%ld", i, j);
fx[count]->SetName(name);
fx[count]->SetLineColor(4);
fxList->Add(fx[count]);
printf(",");
count ++;
}
} }
fxList->Write("EZCurve", TObject::kSingleKey); fxList->Write("EZCurve", TObject::kSingleKey);
printf(" %d e-z finite-size detector functions\n", numEx); printf(" %d (%d) e-z finite-size detector functions\n", numEx, count);
//--- cal modified thetaCM vs z for( int i = 0 ; i < numEx; i++) delete fx[i];
TObjArray * txList = new TObjArray(); delete [] fx;
TGraph ** tx = new TGraph*[numEx]; delete fxList;
for( int j = 0 ; j < numEx; j++){
double a = helios.GetDetRadius();
double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[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 / kbCM[j] * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[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]); // //--- cal modified thetaCM vs z
name.Form("tx%d", j); // TObjArray * txList = new TObjArray();
tx[j]->SetName(name); // TGraph ** tx = new TGraph*[numEx];
tx[j]->SetLineColor(4); // for( int j = 0 ; j < numEx; j++){
txList->Add(tx[j]); // double a = helios.GetDetRadius();
printf("*"); // double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[j] );
} // px.clear();
txList->Write("thetaCM_Z", TObject::kSingleKey); // py.clear();
printf(" %d thetaCM-z for finite-size detector functions\n", numEx); // countfx = 0;
// for(int i = 0; i < 100; i++){
// double thetacm = (i + 8.) * TMath::DegToRad();
// double temp = TMath::TwoPi() * slope / betaRect / kbCM[j] * a / TMath::Sin(thetacm);
// double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[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("thetaCM_Z", TObject::kSingleKey);
// printf(" %d thetaCM-z for finite-size detector functions\n", numEx);
//========timer //========timer
TBenchmark clock; TBenchmark clock;
@ -408,7 +449,7 @@ void Transfer(
shown = false; shown = false;
//change the number of event into human easy-to-read form //change the number of event into human easy-to-read form
int numEvent = reactConfig.numEvents; int numEvent = reactionConfig.numEvents;
int digitLen = TMath::Floor(TMath::Log10(numEvent)); int digitLen = TMath::Floor(TMath::Log10(numEvent));
TString numEventStr; TString numEventStr;
if( 3 <= digitLen && digitLen < 6 ){ if( 3 <= digitLen && digitLen < 6 ){
@ -420,38 +461,40 @@ void Transfer(
} }
printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data()); printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data());
double KEA = reactConfig.beamEnergy; double KEA = reactionConfig.beamEnergy;
double theta = reactConfig.beamTheta; double theta = reactionConfig.beamTheta;
double phi = 0.0; double phi = 0.0;
TF1 * angDist = nullptr;
//*====================================================== calculate event //*====================================================== calculate event
int count = 0; count = 0;
for( int i = 0; i < numEvent; i++){ for( int i = 0; i < numEvent; i++){
bool redoFlag = true; bool redoFlag = true;
if( !reactConfig.isRedo ) redoFlag = false; if( !reactionConfig.isRedo ) redoFlag = false;
do{ do{
rID = gRandom->Integer( numTransfer );
if( !listOfTransfer[rID] ) continue;
//==== Set Ex of B //==== Set Ex of B
if( numEx == 1 ) { ExID = transfer[rID].GetRandomExID();
ExID = 0; double sigma = transfer[rID].GetExList()->ExList[ExID].sigma;
Ex = exList.ExList[0].Ex + (exList.ExList[0].sigma == 0 ? 0 : gRandom->Gaus(0, exList.ExList[0].sigma)); Ex = transfer[rID].GetExList()->ExList[ExID].Ex + gRandom->Gaus(0, sigma);
}else{
ExID = exDist->GetRandom(); transfer[rID].SetExB(Ex);
Ex = exList.ExList[ExID].Ex + (exList.ExList[ExID].sigma == 0 ? 0 : gRandom->Gaus(0, exList.ExList[ExID].sigma));
}
transfer.SetExB(Ex);
//==== Set incident beam //==== Set incident beam
if( reactConfig.beamEnergySigma != 0 ){ if( reactionConfig.beamEnergySigma != 0 ){
KEA = gRandom->Gaus(reactConfig.beamEnergy, reactConfig.beamEnergySigma); KEA = gRandom->Gaus(reactionConfig.beamEnergy, reactionConfig.beamEnergySigma);
} }
if( reactConfig.beamThetaSigma != 0 ){ if( reactionConfig.beamThetaSigma != 0 ){
theta = gRandom->Gaus(reactConfig.beamTheta, reactConfig.beamThetaSigma); theta = gRandom->Gaus(reactionConfig.beamTheta, reactionConfig.beamThetaSigma);
} }
//==== for taregt scattering //==== for taregt scattering
transfer.SetIncidentEnergyAngle(KEA, theta, 0.); transfer[rID].SetIncidentEnergyAngle(KEA, theta, 0.);
transfer.CalReactionConstant(); transfer[rID].CalReactionConstant();
// TLorentzVector PA = transfer.GetPA(); // TLorentzVector PA = transfer.GetPA();
//depth = 0; //depth = 0;
@ -468,8 +511,8 @@ void Transfer(
//==== Calculate thetaCM, phiCM //==== Calculate thetaCM, phiCM
if( distFile->IsOpen()){ if( distFile->IsOpen()){
dist = (TF1 *) distList->At(ExID); angDist = (TF1 *) distList->At(ExID);
thetaCM = dist->GetRandom() / 180. * TMath::Pi(); thetaCM = angDist->GetRandom() / 180. * TMath::Pi();
}else{ }else{
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
} }
@ -477,9 +520,9 @@ void Transfer(
double phiCM = TMath::TwoPi() * gRandom->Rndm(); double phiCM = TMath::TwoPi() * gRandom->Rndm();
//==== Calculate reaction //==== Calculate reaction
transfer.Event(thetaCM, phiCM); transfer[rID].Event(thetaCM, phiCM);
TLorentzVector Pb = transfer.GetPb(); TLorentzVector Pb = transfer[rID].GetPb();
TLorentzVector PB = transfer.GetPB(); TLorentzVector PB = transfer[rID].GetPB();
// //==== Calculate energy loss of scattered and recoil in target // //==== Calculate energy loss of scattered and recoil in target
// if( isTargetScattering ){ // if( isTargetScattering ){
@ -498,14 +541,14 @@ void Transfer(
//======= Decay of particle-B //======= Decay of particle-B
int decayID = 0; int decayID = 0;
if( recoil.isDecay){ if( transfer[rID].GetRecoil().isDecay){
decayID = decay.CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state decayID = decay[rID].CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state
if( decayID == 1 ){ if( decayID == 1 ){
PB = decay.GetDaugther_D(); PB = decay[rID].GetDaugther_D();
//decayTheta = decay.GetAngleChange(); //decayTheta = decay.GetAngleChange();
decayTheta = decay.GetThetaCM(); decayTheta = decay[rID].GetThetaCM();
PB.SetUniqueID(recoil.decayZ); PB.SetUniqueID(transfer[rID].GetRecoil().decayZ);
}else{ }else{
decayTheta = TMath::QuietNaN(); decayTheta = TMath::QuietNaN();
} }
@ -528,19 +571,19 @@ void Transfer(
if( Tb > 0 || TB > 0 ){ if( Tb > 0 || TB > 0 ){
helios.CalArrayHit(Pb); helios[rID].CalArrayHit(Pb);
helios.CalRecoilHit(PB); helios[rID].CalRecoilHit(PB);
hit = 2; hit = 2;
while( hit > 1 ){ hit = helios.CheckDetAcceptance(); } /// while hit > 1, goto next loop; while( hit > 1 ){ hit = helios[rID].CheckDetAcceptance(); } /// while hit > 1, goto next loop;
trajectory orb_b = helios.GetTrajectory_b(); trajectory orb_b = helios[rID].GetTrajectory_b();
trajectory orb_B = helios.GetTrajectory_B(); trajectory orb_B = helios[rID].GetTrajectory_B();
e = helios.GetEnergy() + gRandom->Gaus(0, array.eSigma ); e = helios[rID].GetEnergy() + gRandom->Gaus(0, helios[rID].GetArrayGeometry().eSigma );
double ranX = gRandom->Gaus(0, array.zSigma); double ranX = gRandom->Gaus(0, helios[rID].GetArrayGeometry().zSigma);
z = orb_b.z + ranX; z = orb_b.z + ranX;
detX = helios.GetDetX() + ranX; detX = helios[rID].GetDetX() + ranX;
z0 = orb_b.z0; z0 = orb_b.z0;
t = orb_b.t; t = orb_b.t;
@ -553,17 +596,18 @@ void Transfer(
xArray = orb_b.x; xArray = orb_b.x;
yArray = orb_b.y; yArray = orb_b.y;
//ELUM //ELUM
if( aux.elumPos1 != 0 ){ double elumPos1 = helios[rID].GetAuxGeometry().elumPos1;
xElum1 = helios.GetXPos(aux.elumPos1); if( elumPos1 != 0 ){
yElum1 = helios.GetYPos(aux.elumPos1); xElum1 = helios[rID].GetXPos(elumPos1);
rhoElum1 = helios.GetR(aux.elumPos1); yElum1 = helios[rID].GetYPos(elumPos1);
rhoElum1 = helios[rID].GetR (elumPos1);
} }
if( aux.elumPos2 != 0 ){ double elumPos2 = helios[rID].GetAuxGeometry().elumPos2;
xElum2 = helios.GetXPos(aux.elumPos2); if( elumPos2 ){
yElum2 = helios.GetYPos(aux.elumPos2); xElum2 = helios[rID].GetXPos(elumPos2);
rhoElum2 = helios.GetR(aux.elumPos2); yElum2 = helios[rID].GetYPos(elumPos2);
rhoElum2 = helios[rID].GetR (elumPos2);
} }
//Recoil //Recoil
@ -574,18 +618,20 @@ void Transfer(
rhoB = orb_B.rho; rhoB = orb_B.rho;
//other recoil detectors //other recoil detectors
if ( aux.detPos1 != 0 ){ double recoilPos1 = helios[rID].GetAuxGeometry().detPos1;
xRecoil1 = helios.GetRecoilXPos(aux.detPos1); if ( recoilPos1 != 0 ){
yRecoil1 = helios.GetRecoilYPos(aux.detPos1); xRecoil1 = helios[rID].GetRecoilXPos(recoilPos1);
rhoRecoil1 = helios.GetRecoilR(aux.detPos1); yRecoil1 = helios[rID].GetRecoilYPos(recoilPos1);
rhoRecoil1 = helios[rID].GetRecoilR (recoilPos1);
} }
if ( aux.detPos2 != 0 ){ double recoilPos2 = helios[rID].GetAuxGeometry().detPos2;
xRecoil2 = helios.GetRecoilXPos(aux.detPos2); if ( recoilPos2 != 0 ){
yRecoil2 = helios.GetRecoilYPos(aux.detPos2); xRecoil2 = helios[rID].GetRecoilXPos(recoilPos2);
rhoRecoil2 = helios.GetRecoilR(aux.detPos2); yRecoil2 = helios[rID].GetRecoilYPos(recoilPos2);
rhoRecoil2 = helios[rID].GetRecoilR (recoilPos2);
} }
std::pair<double,double> ExThetaCM = transfer.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius()); std::pair<double,double> ExThetaCM = transfer[rID].CalExThetaCM(e, z, helios[rID].GetBField(), helios[rID].GetDetRadius());
ExCal = ExThetaCM.first; ExCal = ExThetaCM.first;
thetaCMCal = ExThetaCM.second; thetaCMCal = ExThetaCM.second;
@ -593,15 +639,15 @@ void Transfer(
thetaCM = thetaCM * TMath::RadToDeg(); thetaCM = thetaCM * TMath::RadToDeg();
//if decay, get the light decay particle on the recoil; //if decay, get the light decay particle on the recoil;
if( recoil.isDecay ){ if( transfer[rID].GetRecoil().isDecay ){
if( decayID == 1 ){ if( decayID == 1 ){
TLorentzVector Pd = decay.GetDaugther_d(); TLorentzVector Pd = decay[rID].GetDaugther_d();
Td = Pd.E() - Pd.M(); Td = Pd.E() - Pd.M();
helios.CalRecoilHit(Pd); helios[rID].CalRecoilHit(Pd);
trajectory orb_d = helios.GetTrajectory_B(); trajectory orb_d = helios[rID].GetTrajectory_B();
rhoRecoil_d = orb_d.R; rhoRecoil_d = orb_d.R;
xRecoil_d = orb_d.x; xRecoil_d = orb_d.x;
yRecoil_d = orb_d.y; yRecoil_d = orb_d.y;
@ -619,7 +665,7 @@ void Transfer(
if( hit == 1) count ++; if( hit == 1) count ++;
if( reactConfig.isRedo ){ if( reactionConfig.isRedo ){
if( hit == 1) { if( hit == 1) {
redoFlag = false; redoFlag = false;
}else{ }else{
@ -655,11 +701,21 @@ void Transfer(
saveFile->Close(); saveFile->Close();
distFile->Close(); distFile->Close();
delete exDist; delete angDist;
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count); printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
//gROOT->ProcessLine(".q"); //gROOT->ProcessLine(".q");
delete [] transfer;
delete [] decay;
delete [] helios;
distFile->Close();
delete distFile;
delete distList;
delete dwbaExList;
delete dwbaReactList;
return; return;
} }
@ -671,50 +727,49 @@ int main (int argc, char *argv[]) {
printf("========== Simulate Transfer reaction in HELIOS ==========\n"); printf("========== Simulate Transfer reaction in HELIOS ==========\n");
printf("=================================================================\n"); printf("=================================================================\n");
if(argc == 2 || argc > 7) { if(argc == 2 || argc > 5 ) {
printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n"); printf("Usage: ./Transfer [1] [2] [3] [4]\n");
printf(" default file name \n"); printf(" default file name \n");
printf(" [1] reactionConfig.txt (input) reaction Setting \n"); printf(" [1] reactionConfig.txt (input) reaction Setting \n");
printf(" [2] detectorGeo.txt (input) detector Setting \n"); printf(" [2] detectorGeo.txt (input) detector Setting \n");
printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n"); printf(" [3] DWBA.root (input) thetaCM distribution from DWBA \n");
printf(" [4] DWBA.root (input) thetaCM distribution from DWBA \n"); printf(" [4] transfer.root (output) rootFile name for output \n");
printf(" [5] transfer.root (output) rootFile name for output \n");
printf(" [6] plot (input) will it plot stuffs [1/0] \n");
printf("------------------------------------------------------\n"); printf("-----------------------------------------------------------------\n");
printf(" When DWBA.root provided.\n");
printf(" The excitation energies from the DWBA are used.\n");
printf(" And the excitation energies in reactionConfig.txt will be ignored\n");
printf("=================================================================\n");
return 0 ; return 0 ;
} }
std::string basicConfig = "reactionConfig.txt"; std::string basicConfig = "reactionConfig.txt";
std::string heliosDetGeoFile = "detectorGeo.txt"; std::string detGeoFile = "detectorGeo.txt";
int ID = 0;
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root"; TString saveFileName = "transfer.root";
bool isPlot = false; bool isPlot = false;
if( argc >= 2) basicConfig = argv[1]; if( argc >= 2) basicConfig = argv[1];
if( argc >= 3) heliosDetGeoFile = argv[2]; if( argc >= 3) detGeoFile = argv[2];
if( argc >= 4) ID = atoi(argv[3]); if( argc >= 4) ptolemyRoot = argv[3];
if( argc >= 5) ptolemyRoot = argv[4]; if( argc >= 5) saveFileName = argv[4];
if( argc >= 6) saveFileName = argv[5];
if( argc >= 7) isPlot = atoi(argv[7]);
Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName); Transfer( basicConfig, detGeoFile, ptolemyRoot, saveFileName);
//run Armory/Check_Simulation //run Armory/Check_Simulation
if( isPlot ){ // if( isPlot ){
std::ifstream file_in; // std::ifstream file_in;
file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in); // file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in);
if( file_in){ // if( file_in){
printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data()); // printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data());
TString cmd; // TString cmd;
cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data()); // cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data());
system(cmd.Data()); // system(cmd.Data());
}else{ // }else{
printf("cannot find ../Cleopatra/Check_Simulation.C \n"); // printf("cannot find ../Cleopatra/Check_Simulation.C \n");
} // }
} // }
} }

579
Cleopatra/SimTransfer2.C → Cleopatra/SimTransfer_single.C
View File

@ -19,7 +19,9 @@
#include "ClassTransfer.h" #include "ClassTransfer.h"
#include "ClassHelios.h" #include "ClassHelios.h"
double exDistFunc(Double_t *x, Double_t * par){
return par[(int) x[0]];
}
void PrintEZPlotPara(TransferReaction tran, HELIOS helios){ void PrintEZPlotPara(TransferReaction tran, HELIOS helios){
@ -31,100 +33,121 @@ void PrintEZPlotPara(TransferReaction tran, HELIOS helios){
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
double slope = tran.GetEZSlope(helios.GetBField()); /// MeV/mm double slope = tran.GetEZSlope(helios.GetBField()); /// MeV/mm
printf(" e-z slope : %f MeV/mm\n", slope); printf(" e-z slope : %f MeV/mm\n", slope);
// double intercept = q/gamma - mb; // MeV double intercept = q/gamma - mb; // MeV
// printf(" e-z intercept (ground state) : %f MeV\n", intercept); printf(" e-z intercept (ground state) : %f MeV\n", intercept);
} }
void Transfer( void Transfer(
std::string basicConfig = "reactionConfig.txt", std::string basicConfig = "reactionConfig.txt",
std::string detGeoFile = "detectorGeo.txt", std::string heliosDetGeoFile = "detectorGeo.txt",
TString ptolemyRoot = "DWBA.root", unsigned short ID = 0, // this is the ID for the array
TString saveFileName = "transfer.root"){ TString ptolemyRoot = "DWBA.root",
TString saveFileName = "transfer.root"){
//*############################################# Set Reaction //*############################################# Set Reaction
TransferReaction transfer;
HELIOS helios;
Decay decay;
// std::vector<double> kbCM; /// momentum of b in CM frame std::vector<double> kbCM; /// momentum of b in CM frame
// TF1 * exDistribution = nullptr; TF1 * exDist = nullptr;
DetGeo detGeoConfig; transfer.SetReactionFromFile(basicConfig, ID);
ReactionConfig reactionConfig; helios.SetDetectorGeometry(heliosDetGeoFile, ID);
detGeoConfig.LoadDetectorGeo(detGeoFile, false);
reactionConfig.LoadReactionConfig(basicConfig);
const unsigned short numDetGeo = detGeoConfig.array.size();
const unsigned short numReact = reactionConfig.recoil.size();
if( numDetGeo != numReact ){
printf("\e[31m !!!!!! number of array is not equal to number of reaction.!!! \e[0m\n");
printf("Abort\n");
return;
}
unsigned short numTransfer = 0;
for( int i = 0; i < std::min(numDetGeo, numReact); i++){
if( detGeoConfig.array[i].enable ) numTransfer ++;
}
TransferReaction * transfer = new TransferReaction[numTransfer];
Decay * decay = new Decay[numTransfer];
HELIOS * helios = new HELIOS[numTransfer];
int count = 0;
for( unsigned short i = 0 ; i < numDetGeo; i++){
if( detGeoConfig.array[i].enable ){
transfer[count].SetReactionFromFile(basicConfig, i);
if(transfer[count].GetRecoil().isDecay) {
decay[count].SetMotherDaugther(transfer[count].GetRecoil());
}
helios[count].SetDetectorGeometry(detGeoFile, i);
count ++;
}
}
printf("*****************************************************************\n");
printf("*\e[1m\e[33m %27s \e[0m*\n", transfer.GetReactionName().Data());
printf("*****************************************************************\n");
printf("----- loading reaction setting from %s. \n", basicConfig.c_str()); printf("----- loading reaction setting from %s. \n", basicConfig.c_str());
printf("----- loading geometry setting from %s. \n", detGeoFile.c_str()); printf("----- loading geometry setting from %s. \n", heliosDetGeoFile.c_str());
printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n"); printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n");
transfer.PrintReaction(false);
if(transfer.GetRecoil().isDecay) {
decay.SetMotherDaugther(transfer.GetRecoil());
}
helios.PrintGeometry();
PrintEZPlotPara(transfer, helios);
DetGeo detGeo = helios.GetDetectorGeometry();
Array array = helios.GetArrayGeometry();
Auxillary aux = helios.GetAuxGeometry();
ReactionConfig reactConfig = transfer.GetRectionConfig();
Recoil recoil = transfer.GetRecoil();
//*############################################# 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", transfer.GetMass_b(), "mass_b");
// fprintf(keyParaOut, "%-15d //%s\n", reactConfig.recoilLightZ, "charge_b");
// fprintf(keyParaOut, "%-15.8f //%s\n", transfer.GetReactionBeta(), "betaCM");
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.GetCMTotalEnergy(), "Ecm");
// fprintf(keyParaOut, "%-15.4f //%s\n", transfer.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 = reactConfig.isTargetScattering;
// float density = reactConfig.targetDensity;
// float targetThickness = reactConfig.targetThickness;
// if(isTargetScattering) printf("\e[32m#################################### Target Scattering\e[0m\n");
// TargetScattering msA;
// TargetScattering msB;
// TargetScattering msb;
// if(reactConfig.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( reactConfig.isTargetScattering ){
// msA.LoadStoppingPower(reactConfig.beamStoppingPowerFile);
// msb.LoadStoppingPower(reactConfig.recoilLightStoppingPowerFile);
// msB.LoadStoppingPower(reactConfig.recoilHeavyStoppingPowerFile);
// }
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"); TFile * distFile = new TFile(ptolemyRoot, "read");
TObjArray * distList = nullptr; TObjArray * distList = nullptr;
TMacro * dwbaExList = nullptr; TMacro * dwbaExList = nullptr;
TMacro * dwbaReactList = nullptr; TMacro * dwbaReactList = nullptr;
TMacro dwbaExList_Used; TMacro dwbaExList_Used;
TMacro dwbaReactList_Used;
if( distFile->IsOpen() ) { if( distFile->IsOpen() ) {
printf("\e[32m#################################### Load DWBA input : %s \e[0m\n", ptolemyRoot.Data());
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 = (TMacro *) distFile->FindObjectAny("ExList");
dwbaExList_Used.AddLine(dwbaExList->GetListOfLines()->At(0)->GetName()); dwbaExList_Used.AddLine(dwbaExList->GetListOfLines()->At(0)->GetName());
dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList"); dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList");
exList.Clear();
int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ; int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ;
for(int i = 1; i <= numEx ; i++){
for( int i = 0; i < numTransfer; i++){ transfer[i].GetExList()->Clear(); }
for( int i = 1; i <= numEx ; i++){
//Check DWBA reaction is same as transfer setting
std::string reactionName = dwbaReactList->GetListOfLines()->At(i-1)->GetName(); std::string reactionName = dwbaReactList->GetListOfLines()->At(i-1)->GetName();
if( reactionName.find( transfer.GetReactionName().Data() ) != std::string::npos) {
for( int j = 0; j < numTransfer; j++){ std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName();
if( reactionName.find( transfer[j].GetReactionName().Data() ) != std::string::npos) { dwbaExList_Used.AddLine(temp.c_str());
std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName(); if( temp[0] == '/' ) continue;
dwbaReactList_Used.AddLine((reactionName + " | " + std::to_string(j)).c_str()); std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " ");
dwbaExList_Used.AddLine(temp.c_str()); exList.Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
if( temp[0] == '/' ) continue;
std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " ");
transfer[j].GetExList()->Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
}
} }
} }
@ -132,20 +155,25 @@ void Transfer(
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());
} }
std::vector<bool> listOfTransfer(numTransfer, false);
for( int i = 0; i < numTransfer; i++){ printf("------------------------------ Heavy Recoil excitation\n");
if( transfer[i].GetExList()->ExList.size() > 0 ){ printf("Energy[MeV] Rel.Xsec SF sigma\n");
listOfTransfer[i] = true; int numEx = exList.ExList.size();
transfer[i].PrintReaction(false); for( int j = 0; j < numEx; j++){
transfer[i].GetExList()->Print(); double ex = exList.ExList[j].Ex;
helios[i].PrintGeometry(); kbCM.push_back(transfer.CalkCM(ex));
transfer[i].CreateExDistribution(); int decayID = decay.CalDecay(TLorentzVector (0,0,0,0), ex, 0);
// PrintEZPlotPara(transfer[i], helios[i]); exList.ExList[j].Print(decayID == 1 ? "-->Decay" : "\n");
}else{ }
printf(" Reaction : %s has no excited energy. Skipped. \n", transfer[i].GetReactionName().Data());
//---- create Ex-distribution
if( exList.ExList.size() > 1 ) {
printf("---- creating Ex-distribution \n");
exDist = new TF1("exDist", exDistFunc, 0, numEx, numEx);
for(int q = 0; q < numEx; q++){
exDist->SetParameter(q, exList.ExList[q].xsec*exList.ExList[q].SF);
} }
} }
@ -155,22 +183,23 @@ void Transfer(
TTree * tree = new TTree("tree", "tree"); TTree * tree = new TTree("tree", "tree");
TMacro config(basicConfig.c_str()); TMacro config(basicConfig.c_str());
TMacro detGeoTxt(detGeoFile.c_str()); TMacro detGeoTxt(heliosDetGeoFile.c_str());
config.SetName("ReactionConfig"); config.SetName(transfer.GetReactionName_Latex().Data());
config.Write("reactionConfig"); config.Write("reactionConfig");
detGeoTxt.Write("detGeo"); detGeoTxt.Write("detGeo");
if( distList != NULL ) distList->Write("DWBA", 1); if( distList != NULL ) distList->Write("DWBA", 1);
if( dwbaExList != NULL ) { if( dwbaExList != NULL ) dwbaExList_Used.Write("DWBA_ExList", 1);
dwbaExList_Used.Write("DWBA_ExList", 1);
dwbaReactList_Used.Write("DWBA_ReactionList", 1);
}
TMacro idMacro;
idMacro.AddLine(Form("%d", ID));
idMacro.Write("detGeoID");
TMacro hitMeaning; TMacro hitMeaning;
hitMeaning.AddLine("======================= meaning of Hit\n"); hitMeaning.AddLine("======================= meaning of Hit\n");
for( int code = -15 ; code <= 1; code ++ ){ for( int code = -15 ; code <= 1; code ++ ){
hitMeaning.AddLine( Form( "%4d = %s", code, helios[0].AcceptanceCodeToMsg(code).Data() )); hitMeaning.AddLine( Form( "%4d = %s", code, helios.AcceptanceCodeToMsg(code).Data() ));
} }
hitMeaning.AddLine(" other = unknown\n"); hitMeaning.AddLine(" other = unknown\n");
hitMeaning.AddLine("===========================================\n"); hitMeaning.AddLine("===========================================\n");
@ -179,9 +208,6 @@ void Transfer(
int hit; /// the output of Helios.CalHit int hit; /// the output of Helios.CalHit
tree->Branch("hit", &hit, "hit/I"); tree->Branch("hit", &hit, "hit/I");
int rID; /// reaction ID
tree->Branch("rID", &rID, "reactionID/I");
double thetab, phib, Tb; double thetab, phib, Tb;
double thetaB, phiB, TB; double thetaB, phiB, TB;
tree->Branch("thetab", &thetab, "thetab/D"); tree->Branch("thetab", &thetab, "thetab/D");
@ -233,14 +259,7 @@ void Transfer(
double decayTheta; /// the change of thetaB due to decay double decayTheta; /// the change of thetaB due to decay
double xRecoil_d, yRecoil_d, rhoRecoil_d, Td; double xRecoil_d, yRecoil_d, rhoRecoil_d, Td;
if( recoil.isDecay ) {
bool isAnyDecay = false;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyDecay |= transfer[i].GetRecoil().isDecay;
}
if( isAnyDecay ) {
tree->Branch("decayTheta", &decayTheta, "decayTheta/D"); tree->Branch("decayTheta", &decayTheta, "decayTheta/D");
tree->Branch("xRecoil_d", &xRecoil_d, "xRecoil_d/D"); tree->Branch("xRecoil_d", &xRecoil_d, "xRecoil_d/D");
tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D"); tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D");
@ -261,24 +280,14 @@ void Transfer(
///in case need ELUM ///in case need ELUM
double xElum1, yElum1, rhoElum1; double xElum1, yElum1, rhoElum1;
bool isAnyElum1 = false; if( aux.elumPos1 != 0 ) {
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyElum1 |= (helios[i].GetAuxGeometry().elumPos1 != 0);
}
if( isAnyElum1 ) {
tree->Branch("xElum1", &xElum1, "xElum1/D"); tree->Branch("xElum1", &xElum1, "xElum1/D");
tree->Branch("yElum1", &yElum1, "yElum1/D"); tree->Branch("yElum1", &yElum1, "yElum1/D");
tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D"); tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D");
} }
double xElum2, yElum2, rhoElum2; double xElum2, yElum2, rhoElum2;
bool isAnyElum2 = false; if( aux.elumPos2 != 0 ) {
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyElum2 |= (helios[i].GetAuxGeometry().elumPos2 != 0);
}
if( isAnyElum2 ) {
tree->Branch("xElum2", &xElum2, "xElum2/D"); tree->Branch("xElum2", &xElum2, "xElum2/D");
tree->Branch("yElum2", &yElum2, "yElum2/D"); tree->Branch("yElum2", &yElum2, "yElum2/D");
tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D"); tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D");
@ -286,52 +295,31 @@ void Transfer(
///in case need other recoil detector. ///in case need other recoil detector.
double xRecoil1, yRecoil1, rhoRecoil1; double xRecoil1, yRecoil1, rhoRecoil1;
bool isAnyRecoil1 = false; if( aux.detPos1 != 0 ){
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyRecoil1 |= (helios[i].GetAuxGeometry().detPos1 != 0);
}
if( isAnyRecoil1 != 0 ){
tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D"); tree->Branch("xRecoil1", &xRecoil1, "xRecoil1/D");
tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D"); tree->Branch("yRecoil1", &yRecoil1, "yRecoil1/D");
tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D"); tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D");
} }
double xRecoil2, yRecoil2, rhoRecoil2; double xRecoil2, yRecoil2, rhoRecoil2;
bool isAnyRecoil2 = false; if( aux.detPos2 != 0 ){
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
isAnyRecoil2 |= (helios[i].GetAuxGeometry().detPos2 != 0);
}
if( isAnyRecoil2 != 0 ){
tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D"); tree->Branch("xRecoil2", &xRecoil2, "xRecoil2/D");
tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D"); tree->Branch("yRecoil2", &yRecoil2, "yRecoil2/D");
tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D"); tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D");
} }
//======= list of reaction used.
TMacro listOfReaction;
for( int i = 0; i < numTransfer ; i++){
if( !listOfTransfer[i] ) continue;
listOfReaction.AddLine(Form("%2d | %s", i, transfer[i].GetReactionName_Latex().Data()));
}
listOfReaction.Write("ListOfReactions");
//======= function for e-z plot for ideal case //======= function for e-z plot for ideal case
printf("++++ generate functions\n"); printf("++++ generate functions\n");
TObjArray * gList = new TObjArray(); TObjArray * gList = new TObjArray();
gList->SetName("Constant thetaCM = 0 lines"); gList->SetName("Constant thetaCM lines");
const int gxSize = numTransfer; const int gxSize = 50;
TF1 ** gx = new TF1*[gxSize]; TF1 ** gx = new TF1*[gxSize];
TString name; TString name;
for( int i = 0; i < gxSize; i++){ double mb = transfer.GetMass_b();
double mb = transfer[i].GetMass_b(); double betaRect = transfer.GetReactionBeta();
double betaRect = transfer[i].GetReactionBeta(); double gamma = transfer.GetReactionGamma();
double gamma = transfer[i].GetReactionGamma(); double slope = transfer.GetEZSlope(helios.GetBField()); /// MeV/mm
double slope = transfer[i].GetEZSlope(helios[0].GetBField()); /// MeV/mm
for( int i = 0; i < gxSize; i++){
name.Form("g%d", i); name.Form("g%d", i);
gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000); gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000);
double thetacm = i * TMath::DegToRad(); double thetacm = i * TMath::DegToRad();
@ -350,96 +338,67 @@ void Transfer(
gList->Write("EZ_thetaCM", TObject::kSingleKey); gList->Write("EZ_thetaCM", TObject::kSingleKey);
printf(" %d constant thetaCM functions\n", gxSize); printf(" %d constant thetaCM functions\n", gxSize);
for( int i = 0; i < gxSize; i++){
delete gx[i];
}
delete [] gx;
delete gList;
//--- cal modified f //--- cal modified f
int numEx = 0;
for( int i = 0; i < numTransfer; i++){
if( !listOfTransfer[i] ) continue;
numEx += transfer[i].GetExList()->ExList.size();
}
TObjArray * fxList = new TObjArray(); TObjArray * fxList = new TObjArray();
TGraph ** fx = new TGraph*[numEx]; TGraph ** fx = new TGraph*[numEx];
std::vector<double> px, py; std::vector<double> px, py;
int countfx = 0; int countfx = 0;
for( int j = 0 ; j < numEx; j++){
count = 0; double a = helios.GetDetRadius();
for( int i = 0; i < numTransfer; i++ ){ double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[j] );
if( !listOfTransfer[i] ) continue; px.clear();
double mb = transfer[i].GetMass_b(); py.clear();
double betaRect = transfer[i].GetReactionBeta(); countfx = 0;
double gamma = transfer[i].GetReactionGamma(); for(int i = 0; i < 100; i++){
double slope = transfer[i].GetEZSlope(helios[0].GetBField()); /// MeV/mm 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 / kbCM[j] * a / TMath::Sin(thetacm);
for( size_t j = 0 ; j < transfer[i].GetExList()->ExList.size(); j++){ double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[j] * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
double Ex = transfer[i].GetExList()->ExList[j].Ex; double pyTemp = gamma * q - mb - gamma * betaRect * kbCM[j] * TMath::Cos(thetacm);
double kbCM = transfer[i].CalkCM(Ex); if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
double a = helios[i].GetDetRadius(); px.push_back(pxTemp);
double q = TMath::Sqrt(mb*mb + kbCM * kbCM ); py.push_back(pyTemp);
px.clear(); countfx ++;
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 / kbCM * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi()) ;
double pyTemp = gamma * q - mb - gamma * betaRect * kbCM * TMath::Cos(thetacm);
if( TMath::IsNaN(pxTemp) || TMath::IsNaN(pyTemp) ) continue;
px.push_back(pxTemp);
py.push_back(pyTemp);
countfx ++;
}
fx[count] = new TGraph(countfx, &px[0], &py[0]);
name.Form("fx%d_%ld", i, j);
fx[count]->SetName(name);
fx[count]->SetLineColor(4);
fxList->Add(fx[count]);
printf(",");
count ++;
} }
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("EZCurve", TObject::kSingleKey); fxList->Write("EZCurve", TObject::kSingleKey);
printf(" %d (%d) e-z finite-size detector functions\n", numEx, count); printf(" %d e-z finite-size detector functions\n", numEx);
for( int i = 0 ; i < numEx; i++) delete fx[i]; //--- cal modified thetaCM vs z
delete [] fx; TObjArray * txList = new TObjArray();
delete fxList; TGraph ** tx = new TGraph*[numEx];
for( int j = 0 ; j < numEx; j++){
double a = helios.GetDetRadius();
double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[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 / kbCM[j] * a / TMath::Sin(thetacm);
double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[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 ++;
}
// //--- cal modified thetaCM vs z tx[j] = new TGraph(countfx, &px[0], &py[0]);
// TObjArray * txList = new TObjArray(); name.Form("tx%d", j);
// TGraph ** tx = new TGraph*[numEx]; tx[j]->SetName(name);
// for( int j = 0 ; j < numEx; j++){ tx[j]->SetLineColor(4);
// double a = helios.GetDetRadius(); txList->Add(tx[j]);
// double q = TMath::Sqrt(mb*mb + kbCM[j] * kbCM[j] ); printf("*");
// px.clear(); }
// py.clear(); txList->Write("thetaCM_Z", TObject::kSingleKey);
// countfx = 0; printf(" %d thetaCM-z for finite-size detector functions\n", numEx);
// for(int i = 0; i < 100; i++){
// double thetacm = (i + 8.) * TMath::DegToRad();
// double temp = TMath::TwoPi() * slope / betaRect / kbCM[j] * a / TMath::Sin(thetacm);
// double pxTemp = betaRect /slope * (gamma * betaRect * q - gamma * kbCM[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("thetaCM_Z", TObject::kSingleKey);
// printf(" %d thetaCM-z for finite-size detector functions\n", numEx);
//========timer //========timer
TBenchmark clock; TBenchmark clock;
@ -449,7 +408,7 @@ void Transfer(
shown = false; shown = false;
//change the number of event into human easy-to-read form //change the number of event into human easy-to-read form
int numEvent = reactionConfig.numEvents; int numEvent = reactConfig.numEvents;
int digitLen = TMath::Floor(TMath::Log10(numEvent)); int digitLen = TMath::Floor(TMath::Log10(numEvent));
TString numEventStr; TString numEventStr;
if( 3 <= digitLen && digitLen < 6 ){ if( 3 <= digitLen && digitLen < 6 ){
@ -461,40 +420,38 @@ void Transfer(
} }
printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data()); printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data());
double KEA = reactionConfig.beamEnergy; double KEA = reactConfig.beamEnergy;
double theta = reactionConfig.beamTheta; double theta = reactConfig.beamTheta;
double phi = 0.0; double phi = 0.0;
TF1 * angDist = nullptr;
//*====================================================== calculate event //*====================================================== calculate event
count = 0; int count = 0;
for( int i = 0; i < numEvent; i++){ for( int i = 0; i < numEvent; i++){
bool redoFlag = true; bool redoFlag = true;
if( !reactionConfig.isRedo ) redoFlag = false; if( !reactConfig.isRedo ) redoFlag = false;
do{ do{
rID = gRandom->Integer( numTransfer );
if( !listOfTransfer[rID] ) continue;
//==== Set Ex of B //==== Set Ex of B
ExID = transfer[rID].GetRandomExID(); if( numEx == 1 ) {
double sigma = transfer[rID].GetExList()->ExList[ExID].sigma; ExID = 0;
Ex = transfer[rID].GetExList()->ExList[ExID].Ex + gRandom->Gaus(0, sigma); Ex = exList.ExList[0].Ex + (exList.ExList[0].sigma == 0 ? 0 : gRandom->Gaus(0, exList.ExList[0].sigma));
}else{
transfer[rID].SetExB(Ex); ExID = exDist->GetRandom();
Ex = exList.ExList[ExID].Ex + (exList.ExList[ExID].sigma == 0 ? 0 : gRandom->Gaus(0, exList.ExList[ExID].sigma));
}
transfer.SetExB(Ex);
//==== Set incident beam //==== Set incident beam
if( reactionConfig.beamEnergySigma != 0 ){ if( reactConfig.beamEnergySigma != 0 ){
KEA = gRandom->Gaus(reactionConfig.beamEnergy, reactionConfig.beamEnergySigma); KEA = gRandom->Gaus(reactConfig.beamEnergy, reactConfig.beamEnergySigma);
} }
if( reactionConfig.beamThetaSigma != 0 ){ if( reactConfig.beamThetaSigma != 0 ){
theta = gRandom->Gaus(reactionConfig.beamTheta, reactionConfig.beamThetaSigma); theta = gRandom->Gaus(reactConfig.beamTheta, reactConfig.beamThetaSigma);
} }
//==== for taregt scattering //==== for taregt scattering
transfer[rID].SetIncidentEnergyAngle(KEA, theta, 0.); transfer.SetIncidentEnergyAngle(KEA, theta, 0.);
transfer[rID].CalReactionConstant(); transfer.CalReactionConstant();
// TLorentzVector PA = transfer.GetPA(); // TLorentzVector PA = transfer.GetPA();
//depth = 0; //depth = 0;
@ -511,8 +468,8 @@ void Transfer(
//==== Calculate thetaCM, phiCM //==== Calculate thetaCM, phiCM
if( distFile->IsOpen()){ if( distFile->IsOpen()){
angDist = (TF1 *) distList->At(ExID); dist = (TF1 *) distList->At(ExID);
thetaCM = angDist->GetRandom() / 180. * TMath::Pi(); thetaCM = dist->GetRandom() / 180. * TMath::Pi();
}else{ }else{
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
} }
@ -520,9 +477,9 @@ void Transfer(
double phiCM = TMath::TwoPi() * gRandom->Rndm(); double phiCM = TMath::TwoPi() * gRandom->Rndm();
//==== Calculate reaction //==== Calculate reaction
transfer[rID].Event(thetaCM, phiCM); transfer.Event(thetaCM, phiCM);
TLorentzVector Pb = transfer[rID].GetPb(); TLorentzVector Pb = transfer.GetPb();
TLorentzVector PB = transfer[rID].GetPB(); TLorentzVector PB = transfer.GetPB();
// //==== Calculate energy loss of scattered and recoil in target // //==== Calculate energy loss of scattered and recoil in target
// if( isTargetScattering ){ // if( isTargetScattering ){
@ -541,14 +498,14 @@ void Transfer(
//======= Decay of particle-B //======= Decay of particle-B
int decayID = 0; int decayID = 0;
if( transfer[rID].GetRecoil().isDecay){ if( recoil.isDecay){
decayID = decay[rID].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 ){ if( decayID == 1 ){
PB = decay[rID].GetDaugther_D(); PB = decay.GetDaugther_D();
//decayTheta = decay.GetAngleChange(); //decayTheta = decay.GetAngleChange();
decayTheta = decay[rID].GetThetaCM(); decayTheta = decay.GetThetaCM();
PB.SetUniqueID(transfer[rID].GetRecoil().decayZ); PB.SetUniqueID(recoil.decayZ);
}else{ }else{
decayTheta = TMath::QuietNaN(); decayTheta = TMath::QuietNaN();
} }
@ -571,19 +528,19 @@ void Transfer(
if( Tb > 0 || TB > 0 ){ if( Tb > 0 || TB > 0 ){
helios[rID].CalArrayHit(Pb); helios.CalArrayHit(Pb);
helios[rID].CalRecoilHit(PB); helios.CalRecoilHit(PB);
hit = 2; hit = 2;
while( hit > 1 ){ hit = helios[rID].CheckDetAcceptance(); } /// while hit > 1, goto next loop; while( hit > 1 ){ hit = helios.CheckDetAcceptance(); } /// while hit > 1, goto next loop;
trajectory orb_b = helios[rID].GetTrajectory_b(); trajectory orb_b = helios.GetTrajectory_b();
trajectory orb_B = helios[rID].GetTrajectory_B(); trajectory orb_B = helios.GetTrajectory_B();
e = helios[rID].GetEnergy() + gRandom->Gaus(0, helios[rID].GetArrayGeometry().eSigma ); e = helios.GetEnergy() + gRandom->Gaus(0, array.eSigma );
double ranX = gRandom->Gaus(0, helios[rID].GetArrayGeometry().zSigma); double ranX = gRandom->Gaus(0, array.zSigma);
z = orb_b.z + ranX; z = orb_b.z + ranX;
detX = helios[rID].GetDetX() + ranX; detX = helios.GetDetX() + ranX;
z0 = orb_b.z0; z0 = orb_b.z0;
t = orb_b.t; t = orb_b.t;
@ -596,18 +553,17 @@ void Transfer(
xArray = orb_b.x; xArray = orb_b.x;
yArray = orb_b.y; yArray = orb_b.y;
//ELUM //ELUM
double elumPos1 = helios[rID].GetAuxGeometry().elumPos1; if( aux.elumPos1 != 0 ){
if( elumPos1 != 0 ){ xElum1 = helios.GetXPos(aux.elumPos1);
xElum1 = helios[rID].GetXPos(elumPos1); yElum1 = helios.GetYPos(aux.elumPos1);
yElum1 = helios[rID].GetYPos(elumPos1); rhoElum1 = helios.GetR(aux.elumPos1);
rhoElum1 = helios[rID].GetR (elumPos1);
} }
double elumPos2 = helios[rID].GetAuxGeometry().elumPos2; if( aux.elumPos2 != 0 ){
if( elumPos2 ){ xElum2 = helios.GetXPos(aux.elumPos2);
xElum2 = helios[rID].GetXPos(elumPos2); yElum2 = helios.GetYPos(aux.elumPos2);
yElum2 = helios[rID].GetYPos(elumPos2); rhoElum2 = helios.GetR(aux.elumPos2);
rhoElum2 = helios[rID].GetR (elumPos2);
} }
//Recoil //Recoil
@ -618,20 +574,18 @@ void Transfer(
rhoB = orb_B.rho; rhoB = orb_B.rho;
//other recoil detectors //other recoil detectors
double recoilPos1 = helios[rID].GetAuxGeometry().detPos1; if ( aux.detPos1 != 0 ){
if ( recoilPos1 != 0 ){ xRecoil1 = helios.GetRecoilXPos(aux.detPos1);
xRecoil1 = helios[rID].GetRecoilXPos(recoilPos1); yRecoil1 = helios.GetRecoilYPos(aux.detPos1);
yRecoil1 = helios[rID].GetRecoilYPos(recoilPos1); rhoRecoil1 = helios.GetRecoilR(aux.detPos1);
rhoRecoil1 = helios[rID].GetRecoilR (recoilPos1);
} }
double recoilPos2 = helios[rID].GetAuxGeometry().detPos2; if ( aux.detPos2 != 0 ){
if ( recoilPos2 != 0 ){ xRecoil2 = helios.GetRecoilXPos(aux.detPos2);
xRecoil2 = helios[rID].GetRecoilXPos(recoilPos2); yRecoil2 = helios.GetRecoilYPos(aux.detPos2);
yRecoil2 = helios[rID].GetRecoilYPos(recoilPos2); rhoRecoil2 = helios.GetRecoilR(aux.detPos2);
rhoRecoil2 = helios[rID].GetRecoilR (recoilPos2);
} }
std::pair<double,double> ExThetaCM = transfer[rID].CalExThetaCM(e, z, helios[rID].GetBField(), helios[rID].GetDetRadius()); std::pair<double,double> ExThetaCM = transfer.CalExThetaCM(e, z, helios.GetBField(), helios.GetDetRadius());
ExCal = ExThetaCM.first; ExCal = ExThetaCM.first;
thetaCMCal = ExThetaCM.second; thetaCMCal = ExThetaCM.second;
@ -639,15 +593,15 @@ void Transfer(
thetaCM = thetaCM * TMath::RadToDeg(); thetaCM = thetaCM * TMath::RadToDeg();
//if decay, get the light decay particle on the recoil; //if decay, get the light decay particle on the recoil;
if( transfer[rID].GetRecoil().isDecay ){ if( recoil.isDecay ){
if( decayID == 1 ){ if( decayID == 1 ){
TLorentzVector Pd = decay[rID].GetDaugther_d(); TLorentzVector Pd = decay.GetDaugther_d();
Td = Pd.E() - Pd.M(); Td = Pd.E() - Pd.M();
helios[rID].CalRecoilHit(Pd); helios.CalRecoilHit(Pd);
trajectory orb_d = helios[rID].GetTrajectory_B(); trajectory orb_d = helios.GetTrajectory_B();
rhoRecoil_d = orb_d.R; rhoRecoil_d = orb_d.R;
xRecoil_d = orb_d.x; xRecoil_d = orb_d.x;
yRecoil_d = orb_d.y; yRecoil_d = orb_d.y;
@ -665,7 +619,7 @@ void Transfer(
if( hit == 1) count ++; if( hit == 1) count ++;
if( reactionConfig.isRedo ){ if( reactConfig.isRedo ){
if( hit == 1) { if( hit == 1) {
redoFlag = false; redoFlag = false;
}else{ }else{
@ -701,21 +655,11 @@ void Transfer(
saveFile->Close(); saveFile->Close();
distFile->Close(); distFile->Close();
delete angDist; delete exDist;
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count); printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
//gROOT->ProcessLine(".q"); //gROOT->ProcessLine(".q");
delete [] transfer;
delete [] decay;
delete [] helios;
distFile->Close();
delete distFile;
delete distList;
delete dwbaExList;
delete dwbaReactList;
return; return;
} }
@ -727,45 +671,52 @@ int main (int argc, char *argv[]) {
printf("========== Simulate Transfer reaction in HELIOS ==========\n"); printf("========== Simulate Transfer reaction in HELIOS ==========\n");
printf("=================================================================\n"); printf("=================================================================\n");
if(argc == 2 || argc > 5) { if(argc == 2 || argc > 7) {
printf("Usage: ./Transfer [1] [2] [3] [4]\n"); printf("Usage: ./Transfer [1] [2] [3] [4] [5] [6]\n");
printf(" default file name \n"); printf(" default file name \n");
printf(" [1] reactionConfig.txt (input) reaction Setting \n"); printf(" [1] reactionConfig.txt (input) reaction Setting \n");
printf(" [2] detectorGeo.txt (input) detector Setting \n"); printf(" [2] detectorGeo.txt (input) detector Setting \n");
printf(" [3] DWBA.root (input) thetaCM distribution from DWBA \n"); printf(" [3] ID (input) detector & reaction ID (default = 0 ) \n");
printf(" [4] transfer.root (output) rootFile name for output \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"); printf("------------------------------------------------------\n");
return 0 ; return 0 ;
} }
std::string basicConfig = "reactionConfig.txt"; std::string basicConfig = "reactionConfig.txt";
std::string detGeoFile = "detectorGeo.txt"; std::string heliosDetGeoFile = "detectorGeo.txt";
int ID = 0;
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root"; TString saveFileName; // format based on ID;
bool isPlot = false; bool isPlot = false;
if( argc >= 2) basicConfig = argv[1]; if( argc >= 2) basicConfig = argv[1];
if( argc >= 3) detGeoFile = argv[2]; if( argc >= 3) heliosDetGeoFile = argv[2];
if( argc >= 4) ptolemyRoot = argv[3]; if( argc >= 4) ID = atoi(argv[3]);
if( argc >= 5) saveFileName = argv[4]; if( argc >= 5) ptolemyRoot = argv[4];
if( argc >= 6) saveFileName = argv[5];
if( argc >= 7) isPlot = atoi(argv[7]);
Transfer( basicConfig, detGeoFile, ptolemyRoot, saveFileName); saveFileName = Form("transfer_%d.root", ID);
Transfer( basicConfig, heliosDetGeoFile, ID, ptolemyRoot, saveFileName);
//run Armory/Check_Simulation //run Armory/Check_Simulation
// if( isPlot ){ if( isPlot ){
// std::ifstream file_in; std::ifstream file_in;
// file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in); file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in);
// if( file_in){ if( file_in){
// printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data()); printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data());
// TString cmd; TString cmd;
// cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data()); cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data());
// system(cmd.Data()); system(cmd.Data());
// }else{ }else{
// printf("cannot find ../Cleopatra/Check_Simulation.C \n"); printf("cannot find ../Cleopatra/Check_Simulation.C \n");
// } }
// } }
} }

6
Cleopatra/makefile
View File

@ -1,6 +1,6 @@
CC=g++ CC=g++
ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimTransfer2 SimAlpha ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimTransfer_single SimAlpha
all: $(ALL) all: $(ALL)
@ -28,8 +28,8 @@ 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 ../Armory/ClassReactionConfig.h ../Armory/ClassDetGeo.h SimTransfer_single: SimTransfer_single.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) SimTransfer_single.C -o SimTransfer_single `root-config --cflags --glibs`
SimAlpha: SimAlpha.C ClassHelios.h SimAlpha: SimAlpha.C ClassHelios.h
$(CC) SimAlpha.C -o SimAlpha `root-config --cflags --glibs` $(CC) SimAlpha.C -o SimAlpha `root-config --cflags --glibs`

17
Cleopatra/test.C Normal file
View File

@ -0,0 +1,17 @@
#include "ClassHelios.h"
void test() {
// HELIOS helios("../working/detectorGeo.txt", 1);
// helios.GetDetectorGeometry().Print(true);
std::vector<unsigned short> ID;
{
DetGeo temp("../working/detectorGeo.txt");
for( size_t i = 0; i < temp.array.size(); i++ ){
if( temp.array[i].enable ) ID.push_back(i);
}
}
const unsigned short numID = ID.size();
}

10
SOLARIS.sh
View File

@ -16,14 +16,18 @@ export SOLARISANADIR
echo "####### set global variable SOLARISANADIR = ${SOLARISANADIR}" echo "####### set global variable SOLARISANADIR = ${SOLARISANADIR}"
export PATH=$PATH:$SOLARISANADIR/Armory export PATH=$PATH:$SOLARISANADIR/Armory:$SOLARISANADIR/Cleopatra
echo "####### add ${SOLARISANADIR}/Armory into PATH" echo "####### add ${SOLARISANADIR}/Armory into PATH"
echo "####### add ${SOLARISANADIR}/Cleopatra into PATH"
########################### ###########################
echo "####### Define BASH Alias / Functions for SOLARIS" echo "####### Define BASH Alias and Functions for SOLARIS"
echo " 2Working = goto the working directory"
echo " ShowRunTimeStamp = show Run Timestamp"
echo " ShowRunSize = show Run Size"
alias 2Working='cd ${SOLARISANADIR}/working' alias 2Working='cd ${SOLARISANADIR}/working'
alias ShowRunTimeStamp='cat $SOLARISANADIR/data_raw/RunTimeStamp.dat' alias ShowRunTimeStamp='cat $SOLARISANADIR/data_raw/RunTimeStamp.dat'
@ -46,5 +50,3 @@ function ShowRunSize {
fi fi
du -hc $SOLARISANADIR/data_raw/${expName}_${RUN}_*.sol du -hc $SOLARISANADIR/data_raw/${expName}_${RUN}_*.sol
} }
2Working

3
working/Monitor.C
View File

@ -177,7 +177,8 @@ void Monitor::Begin(TTree *tree){
printf("###########################################################\n"); printf("###########################################################\n");
//===================================================== loading parameter //===================================================== loading parameter
corr->LoadDetGeoAndReactionConfigFile();
// corr->LoadDetGeoAndReactionConfigFile();
corr->LoadXNCorr(); corr->LoadXNCorr();
corr->LoadXFXN2ECorr(); corr->LoadXFXN2ECorr();
corr->LoadXScaleCorr(); corr->LoadXScaleCorr();

2
working/detectorGeo.txt
View File

@ -28,7 +28,7 @@ Out //detector_facing_Out_or_In
294.0 294.0
#===============2nd_Array + Recoil #===============2nd_Array + Recoil
true //is_this_array_exist_or_use_for_Simulation false //is_this_array_exist_or_use_for_Simulation
1000 //recoil_position_+_for_downstream_[mm] 1000 //recoil_position_+_for_downstream_[mm]
10.0 //inner_radius_of_recoil_detector_[mm] 10.0 //inner_radius_of_recoil_detector_[mm]
40.2 //outter_radius_of_recoil_detector_[mm] 40.2 //outter_radius_of_recoil_detector_[mm]