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

added SimTransfer2, this can simulate many arrays many reactions at once

Browse Source
This commit is contained in:
Ryan Tang 2024-04-03 16:10:26 -04:00
parent 3f22531698
commit f921d03a39
11 changed files with 837 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

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

3
Armory/ClassDetGeo.h
View File

@ -112,6 +112,7 @@ public:
int BfieldSign ; /// sign of B-field int BfieldSign ; /// sign of B-field
double bore; /// bore , mm double bore; /// bore , mm
unsigned short numEnableGeo;
double zMin, zMax; /// total range span of all arrays double zMin, zMax; /// total range span of all arrays
bool LoadDetectorGeo(TString fileName, bool verbose = true); bool LoadDetectorGeo(TString fileName, bool verbose = true);
@ -214,10 +215,12 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
zMin = 99999; zMin = 99999;
zMax = -99999; zMax = -99999;
numEnableGeo = 0;
for( int i = 0; i < detFlag; i ++ ){ for( int i = 0; i < detFlag; i ++ ){
array[i].DeduceAbsolutePos(); array[i].DeduceAbsolutePos();
if (array[i].enable ) { if (array[i].enable ) {
numEnableGeo ++;
double zmax = TMath::Max(array[i].zMin, array[i].zMax); double zmax = TMath::Max(array[i].zMin, array[i].zMax);
double zmin = TMath::Min(array[i].zMin, array[i].zMax); double zmin = TMath::Min(array[i].zMin, array[i].zMax);
if( zmax > zMax ) zMax = zmax; if( zmax > zMax ) zMax = zmax;

1
Armory/ClassReactionConfig.h
View File

@ -70,7 +70,6 @@ struct ExcitedEnergies {
} }
void Print() const { void Print() const {
printf("...................................\n");
printf("Energy[MeV] Rel.Xsec SF sigma\n"); printf("Energy[MeV] Rel.Xsec SF sigma\n");
for( size_t i = 0; i < ExList.size(); i++){ for( size_t i = 0; i < ExList.size(); i++){
ExList[i].Print("\n"); ExList[i].Print("\n");

67
Cleopatra/ClassTransfer.h
View File

@ -9,6 +9,7 @@
#include "TLorentzVector.h" #include "TLorentzVector.h"
#include "TMath.h" #include "TMath.h"
#include "TF1.h"
//======================================================= //=======================================================
//####################################################### //#######################################################
@ -22,6 +23,7 @@
class TransferReaction { class TransferReaction {
public: public:
TransferReaction(){Inititization();}; TransferReaction(){Inititization();};
TransferReaction(ReactionConfig config, unsigned short ID = 0);
TransferReaction(std::string configFile, unsigned short ID = 0); TransferReaction(std::string configFile, unsigned short ID = 0);
TransferReaction(int beamA, int beamZ, TransferReaction(int beamA, int beamZ,
int targetA, int targetZ, int targetA, int targetZ,
@ -29,7 +31,7 @@ public:
~TransferReaction(); ~TransferReaction();
void SetA(int A, int Z, double Ex); void SetA(int A, int Z, double Ex = 0);
void Seta(int A, int Z); void Seta(int A, int Z);
void Setb(int A, int Z); void Setb(int A, int Z);
void SetB(int A, int Z); void SetB(int A, int Z);
@ -46,9 +48,9 @@ public:
TString GetReactionName() const; TString GetReactionName() const;
TString GetReactionName_Latex(); TString GetReactionName_Latex();
ReactionConfig GetRectionConfig() { return config;} ReactionConfig GetRectionConfig() { return config;}
Recoil GetRecoil() { return recoil;} Recoil GetRecoil() { return recoil;}
ExcitedEnergies GetExList() { return exList;} ExcitedEnergies * GetExList() { return &exList;}
double GetMass_A() const {return mA + ExA;} double GetMass_A() const {return mA + ExA;}
double GetMass_a() const {return ma;} double GetMass_a() const {return ma;}
@ -78,6 +80,22 @@ public:
double GetCMTotalEnergy() {return Etot;} double GetCMTotalEnergy() {return Etot;}
double GetEZSlope(double BField) {return 299.792458 * recoil.lightZ * abs(BField) / TMath::TwoPi() * beta / 1000.;} // MeV/mm double GetEZSlope(double BField) {return 299.792458 * recoil.lightZ * abs(BField) / TMath::TwoPi() * beta / 1000.;} // MeV/mm
void CreateExDistribution();
int GetRandomExID(){
if( exDistribution ) {
return exDistribution->GetRandom();
}
return 0;
}
double GetRandomEx(){
if( exDistribution ) {
int exID = exDistribution->GetRandom();
return exList.ExList[exID].Ex;
}
return TMath::QuietNaN();
}
private: private:
ExcitedEnergies exList; ExcitedEnergies exList;
@ -106,8 +124,33 @@ private:
void Inititization(); void Inititization();
TF1 * exDistribution;
static double exDistFunc(Double_t *x, Double_t * par){
return par[(int) x[0]];
}
}; };
TransferReaction::TransferReaction(ReactionConfig config, unsigned short ID){
Inititization();
SetA(config.beamA, config.beamZ);
Seta(config.targetA, config.targetZ);
SetExA(config.beamEx);
recoil = config.recoil[ID];
exList = config.exList[ID];
Setb(recoil.lightA, recoil.lightZ);
SetB(recoil.heavyA, recoil.heavyZ);
SetIncidentEnergyAngle(config.beamEnergy, 0, 0);
CalReactionConstant();
}
TransferReaction::TransferReaction(std::string configFile, unsigned short ID){ TransferReaction::TransferReaction(std::string configFile, unsigned short ID){
Inititization(); Inititization();
SetReactionFromFile(configFile, ID); SetReactionFromFile(configFile, ID);
@ -133,6 +176,8 @@ void TransferReaction::Inititization(){
TA = 6; TA = 6;
T = TA * config.beamA; T = TA * config.beamA;
exDistribution = nullptr;
ExA = 0; ExA = 0;
ExB = 0; ExB = 0;
@ -147,10 +192,10 @@ void TransferReaction::Inititization(){
} }
TransferReaction::~TransferReaction(){ TransferReaction::~TransferReaction(){
delete exDistribution;
} }
void TransferReaction::SetA(int A, int Z, double Ex = 0){ void TransferReaction::SetA(int A, int Z, double Ex){
Isotope temp (A, Z); Isotope temp (A, Z);
mA = temp.Mass; mA = temp.Mass;
config.beamA = A; config.beamA = A;
@ -463,5 +508,15 @@ std::pair<double, double> TransferReaction::CalExThetaCM(double e, double z, dou
return std::make_pair(Ex, thetaCM); return std::make_pair(Ex, thetaCM);
} }
void TransferReaction::CreateExDistribution(){
int numEx = exList.ExList.size();
exDistribution = new TF1("exDistribution", TransferReaction::exDistFunc, 0, numEx, numEx);
for(int q = 0; q < numEx; q++){
exDistribution->SetParameter(q, exList.ExList[q].xsec*exList.ExList[q].SF);
}
}
#endif #endif

BIN
Cleopatra/Cleopatra
View File

Binary file not shown.

73
Cleopatra/DWBARatio.C
View File

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

16
Cleopatra/DWBA_compare.C
View File

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

3
Cleopatra/ExtractXSec.h
View File

@ -326,9 +326,11 @@ int ExtractXSec (string readFile, int indexForElastic=1) {
//================================== Make TMacro for ExList //================================== Make TMacro for ExList
TMacro ExList; TMacro ExList;
TMacro ReactList;
ExList.AddLine("#---Ex relative_xsec SF sigma_in_MeV"); ExList.AddLine("#---Ex relative_xsec SF sigma_in_MeV");
for( int i = 0; i < numCal ; i++){ for( int i = 0; i < numCal ; i++){
ExList.AddLine(Form("%9.5f %9.5f 1.0 0.000", Ex[i], partialXsec[i])); ExList.AddLine(Form("%9.5f %9.5f 1.0 0.000", Ex[i], partialXsec[i]));
ReactList.AddLine(reaction[i].c_str());
} }
//================================== Save in ROOT //================================== Save in ROOT
@ -369,6 +371,7 @@ int ExtractXSec (string readFile, int indexForElastic=1) {
fList->Write("thetaCM_TF1", 1); fList->Write("thetaCM_TF1", 1);
ExList.Write("ExList"); ExList.Write("ExList");
ReactList.Write("ReactionList");
fileOut->Write(); fileOut->Write();
fileOut->Close(); fileOut->Close();

755
Cleopatra/SimTransfer2.C Normal file
View File

@ -0,0 +1,755 @@
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TF1.h"
#include "TTree.h"
#include "TRandom.h"
#include "TGraph.h"
#include "TMacro.h"
#include <stdlib.h>
#include <vector>
#include <fstream>
#include <TObjArray.h>
#include "../Armory/ClassDetGeo.h"
#include "ClassTargetScattering.h"
#include "ClassDecay.h"
#include "ClassTransfer.h"
#include "ClassHelios.h"
void PrintEZPlotPara(TransferReaction tran, HELIOS helios){
printf("==================================== E-Z plot slope\n");
double betaRect = tran.GetReactionBeta() ;
double gamma = tran.GetReactionGamma();
double mb = tran.GetMass_b();
double pCM = tran.GetMomentumbCM();
double q = TMath::Sqrt(mb*mb + pCM*pCM); ///energy of light recoil in center of mass
double slope = tran.GetEZSlope(helios.GetBField()); /// MeV/mm
printf(" e-z slope : %f MeV/mm\n", slope);
// double intercept = q/gamma - mb; // MeV
// printf(" e-z intercept (ground state) : %f MeV\n", intercept);
}
void Transfer(
std::string basicConfig = "reactionConfig.txt",
std::string detGeoFile = "detectorGeo.txt",
TString ptolemyRoot = "DWBA.root",
TString saveFileName = "transfer.root"){
//*############################################# Set Reaction
// std::vector<double> kbCM; /// momentum of b in CM frame
// TF1 * exDistribution = nullptr;
DetGeo detGeoConfig;
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("----- loading reaction setting from %s. \n", basicConfig.c_str());
printf("----- loading geometry setting from %s. \n", detGeoFile.c_str());
printf("\e[32m#################################### Reaction & HELIOS configuration\e[0m\n");
//*############################################# Load DWBAroot for thetaCM distribution
TFile * distFile = new TFile(ptolemyRoot, "read");
TObjArray * distList = nullptr;
TMacro * dwbaExList = nullptr;
TMacro * dwbaReactList = nullptr;
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");
distList = (TObjArray *) distFile->FindObjectAny("thetaCM_TF1"); // the function List
dwbaExList = (TMacro *) distFile->FindObjectAny("ExList");
dwbaReactList = (TMacro *) distFile->FindObjectAny("ReactionList");
int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ;
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();
for( int j = 0; j < numTransfer; j++){
if( reactionName.find( transfer[j].GetReactionName().Data() ) != std::string::npos) {
std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName();
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);
}
}
}
}else{
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++){
if( transfer[i].GetExList()->ExList.size() > 0 ){
listOfTransfer[i] = true;
transfer[i].PrintReaction(false);
transfer[i].GetExList()->Print();
helios[i].PrintGeometry();
transfer[i].CreateExDistribution();
// PrintEZPlotPara(transfer[i], helios[i]);
}else{
printf(" Reaction : %s has no excited energy. Skipped. \n", transfer[i].GetReactionName().Data());
}
}
//*############################################# build tree
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
TFile * saveFile = new TFile(saveFileName, "recreate");
TTree * tree = new TTree("tree", "tree");
TMacro config(basicConfig.c_str());
TMacro detGeoTxt(detGeoFile.c_str());
config.SetName("ReactionConfig");
config.Write("reactionConfig");
detGeoTxt.Write("detGeo");
if( distList != NULL ) distList->Write("DWBA", 1);
if( dwbaExList != NULL ) dwbaExList->Write("DWBA_ExList", 1);
TMacro hitMeaning;
hitMeaning.AddLine("======================= meaning of Hit\n");
for( int code = -15 ; code <= 1; code ++ ){
hitMeaning.AddLine( Form( "%4d = %s", code, helios[0].AcceptanceCodeToMsg(code).Data() ));
}
hitMeaning.AddLine(" other = unknown\n");
hitMeaning.AddLine("===========================================\n");
hitMeaning.Write("hitMeaning");
int hit; /// the output of Helios.CalHit
tree->Branch("hit", &hit, "hit/I");
int rID; /// reaction ID
tree->Branch("rID", &rID, "reactionID/I");
double thetab, phib, Tb;
double thetaB, phiB, TB;
tree->Branch("thetab", &thetab, "thetab/D");
tree->Branch("phib", &phib, "phib/D");
tree->Branch("Tb", &Tb, "Tb/D");
tree->Branch("thetaB", &thetaB, "thetaB/D");
tree->Branch("phiB", &phiB, "phiB/D");
tree->Branch("TB", &TB, "TB/D");
double thetaCM;
tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
double e, z, detX, t, z0, tB;
tree->Branch("e", &e, "energy_light/D");
tree->Branch("x", &detX, "detector_x/D");
tree->Branch("z", &z, "array_hit_z/D");
tree->Branch("z0", &z0, "z-cycle/D");
tree->Branch("t", &t, "cycle_time_light/D");
tree->Branch("tB", &tB, "recoil_hit_time/D"); /// hit time for recoil on the recoil detector
int loop, detID, detRowID;
tree->Branch("detID", &detID, "detID/I");
tree->Branch("detRowID", &detRowID, "detRowID/I");
tree->Branch("loop", &loop, "loop/I");
double rho, rhoB; ///orbit radius
tree->Branch("rho", &rho, "orbit_radius_light/D");
tree->Branch("rhoB", &rhoB, "orbit_radius_heavy/D");
int ExID;
double Ex;
tree->Branch("ExID", &ExID, "ExID/I");
tree->Branch("Ex", &Ex, "Ex/D");
double ExCal, thetaCMCal;
tree->Branch("ExCal", &ExCal, "ExCal/D");
tree->Branch("thetaCMCal", &thetaCMCal, "thetaCMCal/D");
// double TbLoss; /// energy loss of particle-b from target scattering
// double KEAnew; ///beam energy after target scattering
// double depth; /// reaction depth;
// double Ecm;
// if( reactConfig.isTargetScattering ){
// tree->Branch("depth", &depth, "depth/D");
// tree->Branch("TbLoss", &TbLoss, "TbLoss/D");
// tree->Branch("KEAnew", &KEAnew, "KEAnew/D");
// tree->Branch("Ecm", &Ecm, "Ecm/D");
// }
double decayTheta; /// the change of thetaB due to decay
double xRecoil_d, yRecoil_d, rhoRecoil_d, Td;
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("xRecoil_d", &xRecoil_d, "xRecoil_d/D");
tree->Branch("yRecoil_d", &yRecoil_d, "yRecoil_d/D");
tree->Branch("rhoRecoil_d", &rhoRecoil_d, "rhoRecoil_d/D");
tree->Branch("Td", &Td, "Td/D");
}
double xArray, yArray, rhoArray; ///x, y, rho positon of particle-b on PSD
tree->Branch("xArray", &xArray, "xArray/D");
tree->Branch("yArray", &yArray, "yArray/D");
tree->Branch("rhoArray", &rhoArray, "rhoArray/D");
double xRecoil, yRecoil, rhoRecoil; /// x, y, rho position of particle-B on recoil-detector
tree->Branch("xRecoil", &xRecoil, "xRecoil/D");
tree->Branch("yRecoil", &yRecoil, "yRecoil/D");
tree->Branch("rhoRecoil", &rhoRecoil, "rhoRecoil/D");
///in case need ELUM
double xElum1, yElum1, rhoElum1;
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("yElum1", &yElum1, "yElum1/D");
tree->Branch("rhoElum1", &rhoElum1, "rhoElum1/D");
}
double xElum2, yElum2, rhoElum2;
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("yElum2", &yElum2, "yElum2/D");
tree->Branch("rhoElum2", &rhoElum2, "rhoElum2/D");
}
///in case need other recoil detector.
double xRecoil1, yRecoil1, rhoRecoil1;
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("yRecoil1", &yRecoil1, "yRecoil1/D");
tree->Branch("rhoRecoil1", &rhoRecoil1, "rhoRecoil1/D");
}
double xRecoil2, yRecoil2, rhoRecoil2;
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("yRecoil2", &yRecoil2, "yRecoil2/D");
tree->Branch("rhoRecoil2", &rhoRecoil2, "rhoRecoil2/D");
}
//======= function for e-z plot for ideal case
printf("++++ generate functions\n");
TObjArray * gList = new TObjArray();
gList->SetName("Constant thetaCM = 0 lines");
const int gxSize = numTransfer;
TF1 ** gx = new TF1*[gxSize];
TString name;
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);
gx[i] = new TF1(name, "([0]*TMath::Sqrt([1]+[2]*x*x)+[5]*x)/([3]) - [4]", -1000, 1000);
double thetacm = i * TMath::DegToRad();
double gS2 = TMath::Power(TMath::Sin(thetacm)*gamma,2);
gx[i]->SetParameter(0, TMath::Cos(thetacm));
gx[i]->SetParameter(1, mb*mb*(1-gS2));
gx[i]->SetParameter(2, TMath::Power(slope/betaRect,2));
gx[i]->SetParameter(3, 1-gS2);
gx[i]->SetParameter(4, mb);
gx[i]->SetParameter(5, -gS2*slope);
gx[i]->SetNpx(1000);
gList->Add(gx[i]);
printf("/");
if( i > 1 && i % 40 == 0 ) printf("\n");
}
gList->Write("EZ_thetaCM", TObject::kSingleKey);
printf(" %d constant thetaCM functions\n", gxSize);
for( int i = 0; i < gxSize; i++){
delete gx[i];
}
delete [] gx;
delete gList;
//--- 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();
TGraph ** fx = new TGraph*[numEx];
std::vector<double> px, py;
int countfx = 0;
count = 0;
for( int i = 0; i < numTransfer; i++ ){
if( !listOfTransfer[i] ) continue;
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
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);
printf(" %d (%d) e-z finite-size detector functions\n", numEx, count);
for( int i = 0 ; i < numEx; i++) delete fx[i];
delete [] fx;
delete fxList;
printf("============================== dasjdlasdj\n");
// //--- cal modified thetaCM vs z
// TObjArray * txList = new TObjArray();
// 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 ++;
// }
// 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
TBenchmark clock;
bool shown ;
clock.Reset();
clock.Start("timer");
shown = false;
//change the number of event into human easy-to-read form
int numEvent = reactionConfig.numEvents;
int digitLen = TMath::Floor(TMath::Log10(numEvent));
TString numEventStr;
if( 3 <= digitLen && digitLen < 6 ){
numEventStr.Form("%5.1f kilo", numEvent/1000.);
}else if ( 6<= digitLen && digitLen < 9 ){
numEventStr.Form("%6.2f million", numEvent/1e6);
}else if ( 9<= digitLen ){
numEventStr.Form("%6.2f billion", numEvent/1e9);
}
printf("\e[32m#################################### generating %s events \e[0m\n", numEventStr.Data());
double KEA = reactionConfig.beamEnergy;
double theta = reactionConfig.beamTheta;
double phi = 0.0;
TF1 * angDist = nullptr;
//*====================================================== calculate event
count = 0;
for( int i = 0; i < numEvent; i++){
bool redoFlag = true;
if( !reactionConfig.isRedo ) redoFlag = false;
do{
rID = gRandom->Integer( numTransfer );
if( !listOfTransfer[rID] ) continue;
//==== Set Ex of B
ExID = transfer[rID].GetRandomExID();
double sigma = transfer[rID].GetExList()->ExList[ExID].sigma;
Ex = transfer[rID].GetExList()->ExList[ExID].Ex + gRandom->Gaus(0, sigma);
transfer[rID].SetExB(Ex);
//==== Set incident beam
if( reactionConfig.beamEnergySigma != 0 ){
KEA = gRandom->Gaus(reactionConfig.beamEnergy, reactionConfig.beamEnergySigma);
}
if( reactionConfig.beamThetaSigma != 0 ){
theta = gRandom->Gaus(reactionConfig.beamTheta, reactionConfig.beamThetaSigma);
}
//==== for taregt scattering
transfer[rID].SetIncidentEnergyAngle(KEA, theta, 0.);
transfer[rID].CalReactionConstant();
// TLorentzVector PA = transfer.GetPA();
//depth = 0;
// if( isTargetScattering ){
// //==== Target scattering, only energy loss
// depth = targetThickness * gRandom->Rndm();
// msA.SetTarget(density, depth);
// TLorentzVector PAnew = msA.Scattering(PA);
// KEAnew = msA.GetKE()/reactConfig.beamA;
// transfer.SetIncidentEnergyAngle(KEAnew, theta, phi);
// transfer.CalReactionConstant();
// Ecm = transfer.GetCMTotalKE();
// }
//==== Calculate thetaCM, phiCM
if( distFile->IsOpen()){
angDist = (TF1 *) distList->At(ExID);
thetaCM = angDist->GetRandom() / 180. * TMath::Pi();
}else{
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
}
double phiCM = TMath::TwoPi() * gRandom->Rndm();
//==== Calculate reaction
transfer[rID].Event(thetaCM, phiCM);
TLorentzVector Pb = transfer[rID].GetPb();
TLorentzVector PB = transfer[rID].GetPB();
// //==== Calculate energy loss of scattered and recoil in target
// if( isTargetScattering ){
// if( Pb.Theta() < TMath::PiOver2() ){
// msb.SetTarget(density, targetThickness - depth);
// }else{
// msb.SetTarget(density, depth);
// }
// Pb = msb.Scattering(Pb);
// TbLoss = msb.GetKELoss();
// msB.SetTarget(density, targetThickness - depth);
// PB = msB.Scattering(PB);
// }else{
// TbLoss = 0;
// }
//======= Decay of particle-B
int decayID = 0;
if( transfer[rID].GetRecoil().isDecay){
decayID = decay[rID].CalDecay(PB, Ex, 0, phiCM + TMath::Pi()/2); // decay to ground state
if( decayID == 1 ){
PB = decay[rID].GetDaugther_D();
//decayTheta = decay.GetAngleChange();
decayTheta = decay[rID].GetThetaCM();
PB.SetUniqueID(transfer[rID].GetRecoil().decayZ);
}else{
decayTheta = TMath::QuietNaN();
}
}
//################################### tree branches
//===== reaction
thetab = Pb.Theta() * TMath::RadToDeg();
thetaB = PB.Theta() * TMath::RadToDeg();
Tb = Pb.E() - Pb.M();
TB = PB.E() - PB.M();
phib = Pb.Phi() * TMath::RadToDeg();
phiB = PB.Phi() * TMath::RadToDeg();
//==== Helios
// printf(" thetaCM : %f, Tb : %f\n", thetaCM * TMath::RadToDeg(), Pb.M());
if( Tb > 0 || TB > 0 ){
helios[rID].CalArrayHit(Pb);
helios[rID].CalRecoilHit(PB);
hit = 2;
while( hit > 1 ){ hit = helios[rID].CheckDetAcceptance(); } /// while hit > 1, goto next loop;
trajectory orb_b = helios[rID].GetTrajectory_b();
trajectory orb_B = helios[rID].GetTrajectory_B();
e = helios[rID].GetEnergy() + gRandom->Gaus(0, helios[rID].GetArrayGeometry().eSigma );
double ranX = gRandom->Gaus(0, helios[rID].GetArrayGeometry().zSigma);
z = orb_b.z + ranX;
detX = helios[rID].GetDetX() + ranX;
z0 = orb_b.z0;
t = orb_b.t;
loop = orb_b.loop;
detID = orb_b.detID;
detRowID = orb_b.detRowID;
rho = orb_b.rho;
rhoArray = orb_b.R;
xArray = orb_b.x;
yArray = orb_b.y;
//ELUM
double elumPos1 = helios[rID].GetAuxGeometry().elumPos1;
if( elumPos1 != 0 ){
xElum1 = helios[rID].GetXPos(elumPos1);
yElum1 = helios[rID].GetYPos(elumPos1);
rhoElum1 = helios[rID].GetR (elumPos1);
}
double elumPos2 = helios[rID].GetAuxGeometry().elumPos2;
if( elumPos2 ){
xElum2 = helios[rID].GetXPos(elumPos2);
yElum2 = helios[rID].GetYPos(elumPos2);
rhoElum2 = helios[rID].GetR (elumPos2);
}
//Recoil
rhoRecoil = orb_B.R;
tB = orb_B.t;
xRecoil = orb_B.x;
yRecoil = orb_B.y;
rhoB = orb_B.rho;
//other recoil detectors
double recoilPos1 = helios[rID].GetAuxGeometry().detPos1;
if ( recoilPos1 != 0 ){
xRecoil1 = helios[rID].GetRecoilXPos(recoilPos1);
yRecoil1 = helios[rID].GetRecoilYPos(recoilPos1);
rhoRecoil1 = helios[rID].GetRecoilR (recoilPos1);
}
double recoilPos2 = helios[rID].GetAuxGeometry().detPos2;
if ( recoilPos2 != 0 ){
xRecoil2 = helios[rID].GetRecoilXPos(recoilPos2);
yRecoil2 = helios[rID].GetRecoilYPos(recoilPos2);
rhoRecoil2 = helios[rID].GetRecoilR (recoilPos2);
}
std::pair<double,double> ExThetaCM = transfer[rID].CalExThetaCM(e, z, helios[rID].GetBField(), helios[rID].GetDetRadius());
ExCal = ExThetaCM.first;
thetaCMCal = ExThetaCM.second;
//change thetaCM into deg
thetaCM = thetaCM * TMath::RadToDeg();
//if decay, get the light decay particle on the recoil;
if( transfer[rID].GetRecoil().isDecay ){
if( decayID == 1 ){
TLorentzVector Pd = decay[rID].GetDaugther_d();
Td = Pd.E() - Pd.M();
helios[rID].CalRecoilHit(Pd);
trajectory orb_d = helios[rID].GetTrajectory_B();
rhoRecoil_d = orb_d.R;
xRecoil_d = orb_d.x;
yRecoil_d = orb_d.y;
}else{
rhoRecoil_d = TMath::QuietNaN();
xRecoil_d = TMath::QuietNaN();
yRecoil_d = TMath::QuietNaN();
}
}
}else{
hit = -404;
}
if( hit == 1) count ++;
if( reactionConfig.isRedo ){
if( hit == 1) {
redoFlag = false;
}else{
redoFlag = true;
//printf("%d, %2d, thetaCM : %f, theta : %f, z0: %f \n", i, hit, thetaCM * TMath::RadToDeg(), thetab, helios.GetZ0());
}
}else{
redoFlag = false;
}
}while( redoFlag );
tree->Fill();
//#################################################################### Timer
clock.Stop("timer");
Double_t time = clock.GetRealTime("timer");
clock.Start("timer");
if ( !shown ) {
if (fmod(time, 10) < 1 ){
printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
shown = 1;
}
}else{
if (fmod(time, 10) > 9 ){
shown = 0;
}
}
}
saveFile->Write();
saveFile->Close();
distFile->Close();
delete angDist;
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
//gROOT->ProcessLine(".q");
delete [] transfer;
delete [] decay;
delete [] helios;
distFile->Close();
delete distFile;
delete distList;
delete dwbaExList;
delete dwbaReactList;
return;
}
int main (int argc, char *argv[]) {
printf("=================================================================\n");
printf("========== Simulate Transfer reaction in HELIOS ==========\n");
printf("=================================================================\n");
if(argc == 2 || argc > 5) {
printf("Usage: ./Transfer [1] [2] [3] [4]\n");
printf(" default file name \n");
printf(" [1] reactionConfig.txt (input) reaction Setting \n");
printf(" [2] detectorGeo.txt (input) detector Setting \n");
printf(" [3] DWBA.root (input) thetaCM distribution from DWBA \n");
printf(" [4] transfer.root (output) rootFile name for output \n");
printf("------------------------------------------------------\n");
return 0 ;
}
std::string basicConfig = "reactionConfig.txt";
std::string detGeoFile = "detectorGeo.txt";
TString ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
TString saveFileName = "transfer.root";
bool isPlot = false;
if( argc >= 2) basicConfig = argv[1];
if( argc >= 3) detGeoFile = argv[2];
if( argc >= 4) ptolemyRoot = argv[3];
if( argc >= 5) saveFileName = argv[4];
Transfer( basicConfig, detGeoFile, ptolemyRoot, saveFileName);
//run Armory/Check_Simulation
// if( isPlot ){
// std::ifstream file_in;
// file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in);
// if( file_in){
// printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data());
// TString cmd;
// cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data());
// system(cmd.Data());
// }else{
// printf("cannot find ../Cleopatra/Check_Simulation.C \n");
// }
// }
}

17
Cleopatra/makefile
View File

@ -1,13 +1,13 @@
CC=g++ CC=g++
ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimAlpha ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimTransfer2 SimAlpha
all: $(ALL) all: $(ALL)
Isotope: ../Cleopatra/ClassIsotope.h ../Cleopatra/Isotope.C Isotope: ClassIsotope.h Isotope.C
$(CC) Isotope.C -o Isotope $(CC) Isotope.C -o Isotope
InFileCreator: InFileCreator.C InFileCreator.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h potentials.h InFileCreator: InFileCreator.C InFileCreator.h ClassIsotope.h constant.h potentials.h
$(CC) InFileCreator.C -o InFileCreator `root-config --cflags --glibs` $(CC) InFileCreator.C -o InFileCreator `root-config --cflags --glibs`
ExtractXSec: ExtractXSec.C ExtractXSec.h ExtractXSec: ExtractXSec.C ExtractXSec.h
@ -19,16 +19,19 @@ ExtractXSecFromText: ExtractXSecFromText.C ExtractXSec.h
PlotTGraphTObjArray: PlotTGraphTObjArray.C PlotTGraphTObjArray.h PlotTGraphTObjArray: PlotTGraphTObjArray.C PlotTGraphTObjArray.h
$(CC) PlotTGraphTObjArray.C -o PlotTGraphTObjArray `root-config --cflags --glibs` $(CC) PlotTGraphTObjArray.C -o PlotTGraphTObjArray `root-config --cflags --glibs`
Cleopatra: Cleopatra.C Cleopatra: Cleopatra.C InFileCreator.h ExtractXSec.h
$(CC) Cleopatra.C -o Cleopatra `root-config --cflags --glibs` $(CC) Cleopatra.C -o Cleopatra `root-config --cflags --glibs`
FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h FindThetaCM: FindThetaCM.C FindThetaCM.h ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h
$(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs` $(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs`
SimTransfer: SimTransfer.C ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/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`
SimAlpha: SimAlpha.C ../Cleopatra/ClassHelios.h SimTransfer2: SimTransfer2.C ClassTransfer.h ClassHelios.h ClassIsotope.h constant.h
$(CC) SimTransfer2.C -o SimTransfer2 `root-config --cflags --glibs`
SimAlpha: SimAlpha.C ClassHelios.h
$(CC) SimAlpha.C -o SimAlpha `root-config --cflags --glibs` $(CC) SimAlpha.C -o SimAlpha `root-config --cflags --glibs`
clean: clean:

5
working/DWBA
View File

@ -45,5 +45,6 @@
#36Ar(d,a)34Cl 0 4L=2 3+ 0.000 8MeV/u As # (d,a) reaction #36Ar(d,a)34Cl 0 4L=2 3+ 0.000 8MeV/u As # (d,a) reaction
32Si(d,p)33Si 0 1s1/2 1/2+ 0.000 10MeV/u Vl 32Si(d,p)33Si 0 1s1/2 1/2+ 0.000 10MeV/u AK
32Si(d,p)33Si 0 0d5/2 5/2+ 0.197 10MeV/u Vl 32Si(d,p)33Si 0 0d5/2 5/2+ 0.197 10MeV/u AK
32Si(d,3He)31Al 0 0d5/2 5/2+ 0.000 10MeV/u Ax