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

modified working code for the recent changes

Browse Source
This commit is contained in:
Ryan Tang 2024-02-13 19:25:33 -05:00
parent 96e7d67bbd
commit 871a4d0b26
14 changed files with 411 additions and 411 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.vscode/settings.json vendored
View File

@ -123,7 +123,8 @@
"script_single.C": "cpp",
"script_multi.C": "cpp",
"Isotope.C": "cpp",
"classisotope.h": "c"
"classisotope.h": "c",
"knockoutSim.C": "cpp"
},
"better-comments.multilineComments": true,

39
Armory/AnalysisLib.h
View File

@ -12,9 +12,47 @@
#include <TMath.h>
#include <TObjArray.h>
#include <TCutG.h>
#include <TGraph.h>
namespace AnalysisLib {
//*######################################### TRAPEZOID
TGraph * TrapezoidFilter(TGraph * trace, int baseLineEnd = 80, int riseTime = 10, int flatTop = 20, float decayTime = 2000){
///Trapezoid filter https://doi.org/10.1016/0168-9002(94)91652-7
TGraph * trapezoid = new TGraph();
trapezoid->Clear();
///find baseline;
double baseline = 0;
for( int i = 0; i < baseLineEnd; i++){
baseline += trace->Eval(i);
}
baseline = baseline*1./baseLineEnd;
int length = trace->GetN();
double pn = 0.;
double sn = 0.;
for( int i = 0; i < length ; i++){
double dlk = trace->Eval(i) - baseline;
if( i - riseTime >= 0 ) dlk -= trace->Eval(i - riseTime) - baseline;
if( i - flatTop - riseTime >= 0 ) dlk -= trace->Eval(i - flatTop - riseTime) - baseline;
if( i - flatTop - 2*riseTime >= 0) dlk += trace->Eval(i - flatTop - 2*riseTime) - baseline;
if( i == 0 ){
pn = dlk;
sn = pn + dlk*decayTime;
}else{
pn = pn + dlk;
sn = sn + pn + dlk*decayTime;
}
trapezoid->SetPoint(i, i, sn / decayTime / riseTime);
}
return trapezoid;
}
std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
std::vector<std::string> output;
@ -49,7 +87,6 @@ std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, in
return output;
};
//************************************** TCutG
TObjArray * LoadListOfTCut(TString fileName, TString cutName = "cutList"){

160
Armory/ClassCorrParas.h Normal file
View File

@ -0,0 +1,160 @@
#ifndef Parameters_H
#define Parameters_H
// #include "ClassDetGeo.h"
// #include "ClassReactionConfig.h"
// DetGeo detGeo;
// ReactionConfig reactionConfig1;
// ReactionConfig reactionConfig2;
// void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt",
// std::string reactionConfigFileName1 = "reactionConfig1.txt",
// std::string reactionConfigFileName2 = "reactionConfig2.txt"){
// printf("=====================================================\n");
// printf(" loading detector geometery : %s.", detGeoFileName.c_str());
// TMacro * haha = new TMacro();
// if( haha->ReadFile(detGeoFileName.c_str()) > 0 ) {
// detGeo = AnalysisLib::LoadDetectorGeo(haha);
// printf("... done.\n");
// AnalysisLib::PrintDetGeo(detGeo);
// }else{
// printf("... fail\n");
// }
// delete haha;
// printf("=====================================================\n");
// printf(" loading reaction1 config : %s.", reactionConfigFileName1.c_str());
// TMacro * kaka = new TMacro();
// if( kaka->ReadFile(reactionConfigFileName1.c_str()) > 0 ) {
// reactionConfig1 = AnalysisLib::LoadReactionConfig(kaka);
// printf("..... done.\n");
// AnalysisLib::PrintReactionConfig(reactionConfig1);
// }else{
// printf("..... fail\n");
// }
// delete kaka;
// if( detGeo.use2ndArray){
// printf("=====================================================\n");
// printf(" loading reaction2 config : %s.", reactionConfigFileName2.c_str());
// TMacro * jaja = new TMacro();
// if( jaja->ReadFile(reactionConfigFileName2.c_str()) > 0 ) {
// reactionConfig2 = AnalysisLib::LoadReactionConfig(kaka);
// printf("..... done.\n");
// AnalysisLib::PrintReactionConfig(reactionConfig2);
// }else{
// printf("..... fail\n");
// }
// delete jaja;
// }
// }
//************************************** Correction parameters;
class CorrParas {
public:
CorrParas(){
xnCorr.clear();
xScale.clear();
xfxneCorr.clear();
eCorr.clear();
rdtCorr.clear();
};
std::vector<float> xnCorr; //correction of xn to match xf
std::vector<float> xScale; // correction of x to be (0,1)
std::vector<std::vector<float>> xfxneCorr; //correction of xn and xf to match e
std::vector<std::vector<float>> eCorr; // correction to e, ch -> MeV
std::vector<std::vector<float>> rdtCorr; // correction of rdt, ch -> MeV
//~========================================= xf = xn correction
void LoadXNCorr(bool verbose = false, const char * fileName = "correction_xf_xn.dat"){
printf(" loading xf-xn correction.");
xnCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a;
while( file >> a ) xnCorr.push_back(a);
printf(".......... done.\n");
}else{
printf(".......... fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xnCorr.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
}
//~========================================= X-Scale correction
void LoadXScaleCorr(bool verbose = false, const char * fileName = "correction_scaleX.dat"){
printf(" loading x-Scale correction.");
xScale.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a ) xScale.push_back(a);
printf("........ done.\n");
}else{
printf("........ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xScale.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
}
//~========================================= e = xf + xn correction
void LoadXFXN2ECorr(bool verbose = false, const char * fileName = "correction_xfxn_e.dat"){
printf(" loading xf/xn-e correction.");
xfxneCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a >> b) xfxneCorr.push_back({a, b});
printf("........ done.\n");
}else{
printf("........ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xfxneCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, xfxneCorr[i][0], xfxneCorr[i][1]);
}
//~========================================= e correction
void LoadECorr(bool verbose = false, const char * fileName = "correction_e.dat"){
printf(" loading e correction.");
eCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a >> b) eCorr.push_back( {a, b} ); // 1/a1, a0 , e' = e * a1 + a0
printf(".............. done.\n");
}else{
printf(".............. fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) eCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, eCorr[i][0], eCorr[i][1]);
}
//~========================================= rdt correction
void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.dat"){
printf(" loading rdt correction.");
rdtCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a >> b) rdtCorr.push_back({a, b});
printf("............ done.\n");
}else{
printf("............ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]);
}
};
#endif

12
Armory/ClassDetGeo.h
View File

@ -91,8 +91,8 @@ public:
double zMin, zMax; /// range of detectors
bool isCoincidentWithRecoil;
bool LoadDetectorGeo(TString fileName);
bool LoadDetectorGeo(TMacro * macro);
bool LoadDetectorGeo(TString fileName, bool verbose = true);
bool LoadDetectorGeo(TMacro * macro, bool verbose = true);
void PrintDetGeo( bool printAll = true) const;
@ -108,11 +108,11 @@ inline DetGeo::~DetGeo(){
}
inline bool DetGeo::LoadDetectorGeo(TString fileName){
inline bool DetGeo::LoadDetectorGeo(TString fileName, bool verbose){
TMacro * haha = new TMacro();
if( haha->ReadFile(fileName) > 0 ) {
if( LoadDetectorGeo(haha) ){
if( LoadDetectorGeo(haha, verbose) ){
return true;
}else{
return false;
@ -124,7 +124,7 @@ inline bool DetGeo::LoadDetectorGeo(TString fileName){
///Using TMacro to load the detectorGeo frist,
///this indrect method is good for loading detectorGeo from TMacro in root file
inline bool DetGeo::LoadDetectorGeo(TMacro * macro){
inline bool DetGeo::LoadDetectorGeo(TMacro * macro, bool verbose){
if( macro == NULL ) return false;
@ -211,7 +211,7 @@ inline bool DetGeo::LoadDetectorGeo(TMacro * macro){
zMin = TMath::Min(array1.zMin, array2.zMin);
}
PrintDetGeo(false);
if( verbose ) PrintDetGeo(false);
return true;

158
Armory/ClassReactionParas_tobeKill.h Normal file
View File

@ -0,0 +1,158 @@
#ifndef ReactionParameters_H
#define ReactionParameters_H
#include "ClassDetGeo.h"
class ReactionParas{
public:
ReactionParas();
double Et; // total energy in CM frame
double beta; // Lorentz beta from Lab to CM
double gamma; // Lorentz gamma from Lab to CM
double alpha; // E-Z slope / beta
double G; //The G-coefficient....
double massB; // heavy mass
double q; // charge of light particle
double mass; //light mass
bool hasReactionPara;
double detPrepDist;
void LoadReactionParas(bool verbose = false);
std::pair<double, double> CalExTheta(double e, double z)
};
ReactionParas::ReactionParas(){
}
//~========================================= reaction parameters
inline void ReactionParas::LoadReactionParas(bool verbose = false){
//check is the transfer.root is using the latest reactionConfig.txt
//sicne reaction.dat is generated as a by-product of transfer.root
//TFile * transfer = new TFile("transfer.root");
//TString aaa1 = "";
//TString aaa2 = "";
//if( transfer->IsOpen() ){
// TMacro * reactionConfig = (TMacro *) transfer->FindObjectAny("reactionConfig");
// TMacro presentReactionConfig ("reactionConfig.txt");
// aaa1 = ((TMD5*) reactionConfig->Checksum())->AsString();
// aaa2 = ((TMD5*) presentReactionConfig.Checksum())->AsString();
//}
//printf("%s\n", aaa1.Data());
//printf("%s\n", aaa2.Data());
//if( aaa1 != aaa2 ) {
// printf("########################## recalculate transfer.root \n");
// system("../Cleopatra/Transfer");
// printf("########################## transfer.root updated\n");
//}
std::string fileName;
detPrepDist = Array::detPerpDist;
printf(" loading reaction parameters");
std::ifstream file;
file.open(fileName.c_str());
hasReactionPara = false;
if( file.is_open() ){
std::string x;
int i = 0;
while( file >> x ){
if( x.substr(0,2) == "//" ) continue;
if( i == 0 ) mass = atof(x.c_str());
if( i == 1 ) q = atof(x.c_str());
if( i == 2 ) beta = atof(x.c_str());
if( i == 3 ) Et = atof(x.c_str());
if( i == 4 ) massB = atof(x.c_str());
i = i + 1;
}
printf("........ done.\n");
hasReactionPara = true;
alpha = 299.792458 * abs(detGeo.Bfield) * q / TMath::TwoPi()/1000.; //MeV/mm
gamma = 1./TMath::Sqrt(1-beta * beta);
G = alpha * gamma * beta * detPrepDist ;
if( verbose ){
printf("\tmass-b : %f MeV/c2 \n", mass);
printf("\tcharge-b : %f \n", q);
printf("\tE-total : %f MeV \n", Et);
printf("\tmass-B : %f MeV/c2 \n", massB);
printf("\tbeta : %f \n", beta);
printf("\tB-field : %f T \n", detGeo.Bfield);
printf("\tslope : %f MeV/mm \n", alpha * beta);
printf("\tdet radius: %f mm \n", detPrepDist);
printf("\tG-coeff : %f MeV \n", G);
printf("=====================================================\n");
}
}else{
printf("........ fail.\n");
}
file.close();
}
inline std::pair<double, double> ReactionParas::CalExTheta(double e, double z){
ReactionParas * reactParas = nullptr;
if( detGeo.array1.zMin <= z && z <= detGeo.array1.zMax ){
reactParas = &reactParas1;
if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
}
if( detGeo.array2.zMin <= z && z <= detGeo.array2.zMax ){
reactParas = &reactParas2;
if( !hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
}
double Ex = TMath::QuietNaN();
double thetaCM = TMath::QuietNaN();
double y = e + mass; // to give the KE + mass of proton;
double Z = alpha * gamma * beta * z;
double H = TMath::Sqrt(TMath::Power(gamma * beta,2) * (y*y - mass * mass) ) ;
if( TMath::Abs(Z) < H ) {
///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi)
double tolerrence = 0.001;
double phi = 0; ///initial phi = 0 -> ensure the solution has f'(phi) > 0
double nPhi = 0; /// new phi
int iter = 0;
do{
phi = nPhi;
nPhi = phi - (H * TMath::Sin(phi) - G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - G /TMath::Power( TMath::Cos(phi), 2));
iter ++;
if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break;
}while( TMath::Abs(phi - nPhi ) > tolerrence);
phi = nPhi;
/// check f'(phi) > 0
double Df = H * TMath::Cos(phi) - G / TMath::Power( TMath::Cos(phi),2);
if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){
double K = H * TMath::Sin(phi);
double x = TMath::ACos( mass / ( y * gamma - K));
double momt = mass * TMath::Tan(x); /// momentum of particel b or B in CM frame
double EB = TMath::Sqrt(mass * mass + Et * Et - 2 * Et * TMath::Sqrt(momt*momt + mass * mass));
Ex = EB - massB;
double hahaha1 = gamma * TMath::Sqrt(mass * mass + momt * momt) - y;
double hahaha2 = gamma * beta * momt;
thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg();
}
}
return std::make_pair(Ex, thetaCM);
}
#endif

43
Armory/GeneralSort.h
View File

@ -50,49 +50,6 @@ double fitFunc(double * x, double * par){
return par[3] + par[0] * (1 - TMath::Exp(- (x[0] - par[1]) / par[2]) ) * TMath::Exp(- (x[0] - par[1]) / par[4]);
}
//^######################################### TRAPEZOID
TGraph * TrapezoidFilter(TGraph * trace){
///Trapezoid filter https://doi.org/10.1016/0168-9002(94)91652-7
//TODO how to not hard code?
int baseLineEnd = 80;
int riseTime = 10; //ch
int flatTop = 20;
float decayTime = 2000;
TGraph * trapezoid = new TGraph();
trapezoid->Clear();
///find baseline;
double baseline = 0;
for( int i = 0; i < baseLineEnd; i++){
baseline += trace->Eval(i);
}
baseline = baseline*1./baseLineEnd;
int length = trace->GetN();
double pn = 0.;
double sn = 0.;
for( int i = 0; i < length ; i++){
double dlk = trace->Eval(i) - baseline;
if( i - riseTime >= 0 ) dlk -= trace->Eval(i - riseTime) - baseline;
if( i - flatTop - riseTime >= 0 ) dlk -= trace->Eval(i - flatTop - riseTime) - baseline;
if( i - flatTop - 2*riseTime >= 0) dlk += trace->Eval(i - flatTop - 2*riseTime) - baseline;
if( i == 0 ){
pn = dlk;
sn = pn + dlk*decayTime;
}else{
pn = pn + dlk;
sn = sn + pn + dlk*decayTime;
}
trapezoid->SetPoint(i, i, sn / decayTime / riseTime);
}
return trapezoid;
}
TStopwatch stpWatch;
//^######################################### Class definition

306
Armory/Parameters.h
View File

@ -1,306 +0,0 @@
#ifndef Parameters_H
#define Parameters_H
#include "ClassDetGeo.h"
#include "ClassReactionConfig.h"
struct ReactionParas{
double Et; // total energy in CM frame
double beta; // Lorentz beta from Lab to CM
double gamma; // Lorentz gamma from Lab to CM
double alpha; // E-Z slope / beta
double G; //The G-coefficient....
double massB; // heavy mass
double q; // charge of light particle
double mass; //light mass
bool hasReactionPara;
double detPrepDist;
};
ReactionParas reactParas1;
ReactionParas reactParas2;
//~========================================= reaction parameters
void LoadReactionParas(int arrayID, bool verbose = false){
//check is the transfer.root is using the latest reactionConfig.txt
//sicne reaction.dat is generated as a by-product of transfer.root
//TFile * transfer = new TFile("transfer.root");
//TString aaa1 = "";
//TString aaa2 = "";
//if( transfer->IsOpen() ){
// TMacro * reactionConfig = (TMacro *) transfer->FindObjectAny("reactionConfig");
// TMacro presentReactionConfig ("reactionConfig.txt");
// aaa1 = ((TMD5*) reactionConfig->Checksum())->AsString();
// aaa2 = ((TMD5*) presentReactionConfig.Checksum())->AsString();
//}
//printf("%s\n", aaa1.Data());
//printf("%s\n", aaa2.Data());
//if( aaa1 != aaa2 ) {
// printf("########################## recalculate transfer.root \n");
// system("../Cleopatra/Transfer");
// printf("########################## transfer.root updated\n");
//}
ReactionParas * reactParas = nullptr;
std::string fileName;
if( arrayID == 1){
reactParas = &AnalysisLib::reactParas1;
fileName = "reaction.dat";
}else if( arrayID == 2){
reactParas = &AnalysisLib::reactParas2;
fileName = "reaction2.dat";
}else{
printf("arrayID must be either 1 or 2. Abort.\n");
return;
}
reactParas->detPrepDist = AnalysisLib::detGeo.array1.detPerpDist;
printf(" loading reaction parameters");
std::ifstream file;
file.open(fileName.c_str());
reactParas->hasReactionPara = false;
if( file.is_open() ){
std::string x;
int i = 0;
while( file >> x ){
if( x.substr(0,2) == "//" ) continue;
if( i == 0 ) reactParas->mass = atof(x.c_str());
if( i == 1 ) reactParas->q = atof(x.c_str());
if( i == 2 ) reactParas->beta = atof(x.c_str());
if( i == 3 ) reactParas->Et = atof(x.c_str());
if( i == 4 ) reactParas->massB = atof(x.c_str());
i = i + 1;
}
printf("........ done.\n");
reactParas->hasReactionPara = true;
reactParas->alpha = 299.792458 * abs(detGeo.Bfield) * reactParas->q / TMath::TwoPi()/1000.; //MeV/mm
reactParas->gamma = 1./TMath::Sqrt(1-reactParas->beta * reactParas->beta);
reactParas->G = reactParas->alpha * reactParas->gamma * reactParas->beta * reactParas->detPrepDist ;
if( verbose ){
printf("\tmass-b : %f MeV/c2 \n", reactParas->mass);
printf("\tcharge-b : %f \n", reactParas->q);
printf("\tE-total : %f MeV \n", reactParas->Et);
printf("\tmass-B : %f MeV/c2 \n", reactParas->massB);
printf("\tbeta : %f \n", reactParas->beta);
printf("\tB-field : %f T \n", detGeo.Bfield);
printf("\tslope : %f MeV/mm \n", reactParas->alpha * reactParas->beta);
printf("\tdet radius: %f mm \n", reactParas->detPrepDist);
printf("\tG-coeff : %f MeV \n", reactParas->G);
printf("=====================================================\n");
}
}else{
printf("........ fail.\n");
}
file.close();
}
std::vector<double> CalExTheta(double e, double z){
ReactionParas * reactParas = nullptr;
if( detGeo.array1.zMin <= z && z <= detGeo.array1.zMax ){
reactParas = &reactParas1;
if( !reactParas->hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
}
if( detGeo.array2.zMin <= z && z <= detGeo.array2.zMax ){
reactParas = &reactParas2;
if( !reactParas->hasReactionPara) return {TMath::QuietNaN(), TMath::QuietNaN()};
}
double Ex = TMath::QuietNaN();
double thetaCM = TMath::QuietNaN();
double y = e + reactParas->mass; // to give the KE + mass of proton;
double Z = reactParas->alpha * reactParas->gamma * reactParas->beta * z;
double H = TMath::Sqrt(TMath::Power(reactParas->gamma * reactParas->beta,2) * (y*y - reactParas->mass * reactParas->mass) ) ;
if( TMath::Abs(Z) < H ) {
///using Newton's method to solve 0 == H * sin(phi) - G * tan(phi) - Z = f(phi)
double tolerrence = 0.001;
double phi = 0; ///initial phi = 0 -> ensure the solution has f'(phi) > 0
double nPhi = 0; /// new phi
int iter = 0;
do{
phi = nPhi;
nPhi = phi - (H * TMath::Sin(phi) - reactParas->G * TMath::Tan(phi) - Z) / (H * TMath::Cos(phi) - reactParas->G /TMath::Power( TMath::Cos(phi), 2));
iter ++;
if( iter > 10 || TMath::Abs(nPhi) > TMath::PiOver2()) break;
}while( TMath::Abs(phi - nPhi ) > tolerrence);
phi = nPhi;
/// check f'(phi) > 0
double Df = H * TMath::Cos(phi) - reactParas->G / TMath::Power( TMath::Cos(phi),2);
if( Df > 0 && TMath::Abs(phi) < TMath::PiOver2() ){
double K = H * TMath::Sin(phi);
double x = TMath::ACos( reactParas->mass / ( y * reactParas->gamma - K));
double momt = reactParas->mass * TMath::Tan(x); /// momentum of particel b or B in CM frame
double EB = TMath::Sqrt(reactParas->mass * reactParas->mass + reactParas->Et * reactParas->Et - 2 * reactParas->Et * TMath::Sqrt(momt*momt + reactParas->mass * reactParas->mass));
Ex = EB - reactParas->massB;
double hahaha1 = reactParas->gamma * TMath::Sqrt(reactParas->mass * reactParas->mass + momt * momt) - y;
double hahaha2 = reactParas->gamma * reactParas->beta * momt;
thetaCM = TMath::ACos(hahaha1/hahaha2) * TMath::RadToDeg();
}
}
return {Ex, thetaCM};
}
DetGeo detGeo;
ReactionConfig reactionConfig1;
ReactionConfig reactionConfig2;
void LoadDetGeoAndReactionConfigFile(std::string detGeoFileName = "detectorGeo.txt",
std::string reactionConfigFileName1 = "reactionConfig1.txt",
std::string reactionConfigFileName2 = "reactionConfig2.txt"){
printf("=====================================================\n");
printf(" loading detector geometery : %s.", detGeoFileName.c_str());
TMacro * haha = new TMacro();
if( haha->ReadFile(detGeoFileName.c_str()) > 0 ) {
detGeo = AnalysisLib::LoadDetectorGeo(haha);
printf("... done.\n");
AnalysisLib::PrintDetGeo(detGeo);
}else{
printf("... fail\n");
}
delete haha;
printf("=====================================================\n");
printf(" loading reaction1 config : %s.", reactionConfigFileName1.c_str());
TMacro * kaka = new TMacro();
if( kaka->ReadFile(reactionConfigFileName1.c_str()) > 0 ) {
reactionConfig1 = AnalysisLib::LoadReactionConfig(kaka);
printf("..... done.\n");
AnalysisLib::PrintReactionConfig(reactionConfig1);
}else{
printf("..... fail\n");
}
delete kaka;
if( detGeo.use2ndArray){
printf("=====================================================\n");
printf(" loading reaction2 config : %s.", reactionConfigFileName2.c_str());
TMacro * jaja = new TMacro();
if( jaja->ReadFile(reactionConfigFileName2.c_str()) > 0 ) {
reactionConfig2 = AnalysisLib::LoadReactionConfig(kaka);
printf("..... done.\n");
AnalysisLib::PrintReactionConfig(reactionConfig2);
}else{
printf("..... fail\n");
}
delete jaja;
}
}
//************************************** Correction parameters;
std::vector<float> xnCorr; //correction of xn to match xf
std::vector<float> xScale; // correction of x to be (0,1)
std::vector<std::vector<float>> xfxneCorr; //correction of xn and xf to match e
std::vector<std::vector<float>> eCorr; // correction to e, ch -> MeV
std::vector<std::vector<float>> rdtCorr; // correction of rdt, ch -> MeV
//~========================================= xf = xn correction
void LoadXNCorr(bool verbose = false, const char * fileName = "correction_xf_xn.dat"){
printf(" loading xf-xn correction.");
xnCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a;
while( file >> a ) xnCorr.push_back(a);
printf(".......... done.\n");
}else{
printf(".......... fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xnCorr.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
}
//~========================================= X-Scale correction
void LoadXScaleCorr(bool verbose = false, const char * fileName = "correction_scaleX.dat"){
printf(" loading x-Scale correction.");
xScale.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a ) xScale.push_back(a);
printf("........ done.\n");
}else{
printf("........ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xScale.size(); i++) printf("%2d | %10.3f\n", i, xnCorr[i]);
}
//~========================================= e = xf + xn correction
void LoadXFXN2ECorr(bool verbose = false, const char * fileName = "correction_xfxn_e.dat"){
printf(" loading xf/xn-e correction.");
xfxneCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a >> b) xfxneCorr.push_back({a, b});
printf("........ done.\n");
}else{
printf("........ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) xfxneCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, xfxneCorr[i][0], xfxneCorr[i][1]);
}
//~========================================= e correction
void LoadECorr(bool verbose = false, const char * fileName = "correction_e.dat"){
printf(" loading e correction.");
eCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a >> b) eCorr.push_back( {a, b} ); // 1/a1, a0 , e' = e * a1 + a0
printf(".............. done.\n");
}else{
printf(".............. fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) eCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, eCorr[i][0], eCorr[i][1]);
}
//~========================================= rdt correction
void LoadRDTCorr(bool verbose = false, const char * fileName = "correction_rdt.dat"){
printf(" loading rdt correction.");
rdtCorr.clear();
std::ifstream file;
file.open(fileName);
if( file.is_open() ){
float a, b;
while( file >> a >> b) rdtCorr.push_back({a, b});
printf("............ done.\n");
}else{
printf("............ fail.\n");
}
file.close();
if( verbose ) for(int i = 0; i < (int) rdtCorr.size(); i++) printf("%2d | %10.3f, %10.3f\n", i, rdtCorr[i][0], rdtCorr[i][1]);
}
#endif

8
Cleopatra/ClassTransfer.h
View File

@ -36,7 +36,7 @@ public:
void SetExA(double Ex);
void SetExB(double Ex);
void SetReactionFromFile(string settingFile);
void SetReactionFromFile(string reactionConfigFile);
TString GetReactionName();
TString GetReactionName_Latex();
@ -198,9 +198,9 @@ void TransferReaction::SetExB(double Ex){
isReady = false;
}
void TransferReaction::SetReactionFromFile(string settingFile){
void TransferReaction::SetReactionFromFile(string reactionConfigFile){
if( reaction.LoadReactionConfig(settingFile) ){
if( reaction.LoadReactionConfig(reactionConfigFile) ){
SetA(reaction.beamA, reaction.beamZ);
Seta(reaction.targetA, reaction.targetZ);
@ -211,7 +211,7 @@ void TransferReaction::SetReactionFromFile(string settingFile){
CalReactionConstant();
}else{
printf("cannot read file %s.\n", settingFile.c_str());
printf("cannot read file %s.\n", reactionConfigFile.c_str());
isReady = false;
}

9
Cleopatra/Simulation_Helper.C
View File

@ -16,8 +16,8 @@
#include "../Cleopatra/InFileCreator.h"
#include "../Cleopatra/ExtractXSec.h"
#include "../Cleopatra/PlotTGraphTObjArray.h"
#include "../armory/AutoFit.C"
#include "../armory/AnalysisLib.h"
#include "../Armory/AutoFit.C"
#include "../Armory/AnalysisLib.h"
#include "../Cleopatra/Check_Simulation.C"
#include <iostream>
@ -358,9 +358,8 @@ bool MyMainFrame::IsFileExist(TString filename){
void MyMainFrame::CheckIsUse2ndArray(){
TMacro * haha = new TMacro("../working/detectorGeo.txt");
AnalysisLib::DetGeo detGeo = AnalysisLib::LoadDetectorGeo(haha);
delete haha;
DetGeo detGeo;
detGeo.LoadDetectorGeo("../working/detectorGeo.txt", false);
isUse2ndArray = detGeo.use2ndArray;
}

2
Cleopatra/alpha.C
View File

@ -1,4 +1,4 @@
#include "HELIOS_LIB.h"
#include "ClassHelios.h"
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"

2
Cleopatra/knockout.C → Cleopatra/knockoutSim.C
View File

@ -1,4 +1,4 @@
#include "HELIOS_LIB.h"
#include "ClassHelio.h"
#include "TROOT.h"
#include "TBenchmark.h"
#include "TLorentzVector.h"

0
working/DWBA2 → working/DWBA
View File

51
working/Monitor.C
View File

@ -21,7 +21,7 @@
#include <TObjArray.h>
#include <fstream>
#include <vector>
#include "../Cleopatra/Isotope.h"
#include "../Cleopatra/ClassIsotope.h"
#include "Mapping.h"
#define tick2ns 8. // 1clock tick = 8 ns
@ -177,27 +177,25 @@ void Monitor::Begin(TTree *tree){
printf("###########################################################\n");
//===================================================== loading parameter
AnalysisLib::LoadDetGeoAndReactionConfigFile();
AnalysisLib::LoadXNCorr();
AnalysisLib::LoadXFXN2ECorr();
AnalysisLib::LoadXScaleCorr();
AnalysisLib::LoadECorr();
AnalysisLib::LoadRDTCorr();
AnalysisLib::LoadReactionParas(1, true);
if( AnalysisLib::detGeo.use2ndArray ) AnalysisLib::LoadReactionParas(2, true);
corr->LoadDetGeoAndReactionConfigFile();
corr->LoadXNCorr();
corr->LoadXFXN2ECorr();
corr->LoadXScaleCorr();
corr->LoadECorr();
corr->LoadRDTCorr();
if( (int) AnalysisLib::xnCorr.size() < mapping::NARRAY ) { isXNCorrOK = false; printf(" !!!!!!!! size of xnCorr < NARRAY .\n"); }
if( (int) AnalysisLib::xfxneCorr.size() < mapping::NARRAY ) { isXFXNCorrOK = false; printf(" !!!!!!!! size of xfxneCorr < NARRAY .\n"); }
if( (int) AnalysisLib::eCorr.size() < mapping::NARRAY ) { isXScaleCorrOK = false; printf(" !!!!!!!! size of eCorr < NARRAY .\n"); }
if( (int) AnalysisLib::xScale.size() < mapping::NARRAY ) { isECorrOK = false; printf(" !!!!!!!! size of xScale < NARRAY .\n"); }
if( (int) AnalysisLib::rdtCorr.size() < mapping::NRDT ) { isRDTCorrOK = false; printf(" !!!!!!!! size of rdtCorr < NRDT .\n"); }
if( (int) corr->xnCorr.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of xnCorr < NARRAY .\n"); }
if( (int) corr->xfxneCorr.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of xfxneCorr < NARRAY .\n"); }
if( (int) corr->eCorr.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of eCorr < NARRAY .\n"); }
if( (int) corr->xScale.size() < mapping::NARRAY ) { printf(" !!!!!!!! size of xScale < NARRAY .\n"); }
if( (int) corr->rdtCorr.size() < mapping::NRDT ) { printf(" !!!!!!!! size of rdtCorr < NRDT .\n"); }
numRow = AnalysisLib::detGeo.array1.nDet;
numRow = detGeo->use2ndArray ? detGeo->array2.nDet : detGeo->array1.nDet;
numCol = mapping::NARRAY/numRow;
numDet = mapping::NARRAY;
zRange[0] = AnalysisLib::detGeo.zMax - 50;
zRange[1] = AnalysisLib::detGeo.zMax + 50;
zRange[0] = detGeo->zMax - 50;
zRange[1] = detGeo->zMax + 50;
printf("=====================================================\n");
printf(" z Range : %5.0f - %5.0f mm\n", zRange[0], zRange[1]);
@ -213,7 +211,6 @@ void Monitor::Begin(TTree *tree){
//================ Get EZ cuts;
EZCut = AnalysisLib::LoadSingleTCut(ezCutFile);
//========================= Generate all of the histograms needed for drawing later on
printf("============================================ Histograms declaration\n");
@ -447,7 +444,8 @@ Bool_t Monitor::Process(Long64_t entry){
if( isXScaleCorrOK ) xCal[detID] = (xCal[detID]-0.5)/AnalysisLib::xScale[detID] + 0.5; /// if include this scale, need to also inclused in Cali_littleTree
if( abs(xCal[detID] - 0.5) > xGate/2. ) continue;
//TODO two arrays?
//@==================== calculate Z
if( AnalysisLib::detGeo.array1.firstPos > 0 ) {
z[detID] = AnalysisLib::detGeo.array1.detLength*(1.0-xCal[detID]) + AnalysisLib::detGeo.array1.detPos[detID%numCol];
@ -569,18 +567,9 @@ Bool_t Monitor::Process(Long64_t entry){
if( eCal[detID] < eCalCut[0] ) continue ;
if( eCal[detID] > eCalCut[1] ) continue ;
double Ex, thetaCM;
if( AnalysisLib::hasReactionPara ){
std::vector<double> ExThetaCM = AnalysisLib::CalExTheta(eCal[detID], x[detID]);
Ex = ExThetaCM[0];
thetaCM = ExThetaCM[1];
}else{
Ex = TMath::QuietNaN();
thetaCM = TMath::QuietNaN();
}
std::pair<double, double> ExThetaCM = transfer->CalExThetaCM(eCal[detID], x[detID], detGeo->Bfield, detGeo->array1.detPerpDist);
double Ex = ExThetaCM.first;
double thetaCM = ExThetaCM.second;
if( thetaCM > thetaCMGate ) {

29
working/Monitor.h
View File

@ -12,7 +12,11 @@
#include <TCutG.h>
#include "Mapping.h"
#include "../armory/AnalysisLib.h"
#include "../Armory/AnalysisLib.h"
#include "../Armory/ClassDetGeo.h"
#include "../Armory/ClassReactionConfig.h"
#include "../Armory/ClassCorrParas.h"
#include "../Cleopatra/ClassTransfer.h"
class Monitor : public TSelector {
public :
@ -58,11 +62,10 @@ public :
bool isRDTExist;
bool isXNCorrOK;
bool isXFXNCorrOK;
bool isXScaleCorrOK;
bool isECorrOK;
bool isRDTCorrOK;
CorrParas * corr; //!
DetGeo * detGeo; //!
TransferReaction * transfer; //!
//==== global variable
float * x, * z;
@ -95,12 +98,6 @@ public :
xnCal = new float [mapping::NARRAY];
eCal = new float [mapping::NARRAY];
isXNCorrOK = true;
isXFXNCorrOK = true;
isXScaleCorrOK = true;
isECorrOK = true;
isRDTCorrOK = true;
padID = 0;
timeRangeInMin[0] = 0;
@ -112,6 +109,12 @@ public :
baseTimeStamp = 0;
treeID = -1;
corr = new CorrParas();
detGeo = new DetGeo();
detGeo->LoadDetectorGeo("detectorGeo.txt");
transfer = new TransferReaction();
transfer->SetReactionFromFile("reactionConfig1.txt");
}
virtual ~Monitor() {
@ -131,6 +134,8 @@ public :
delete xnCal;
delete eCal;
delete corr;
}
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);