From 871a4d0b26d62073fc7eae6fcf9aaff7aa0586bc Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Tue, 13 Feb 2024 19:25:33 -0500 Subject: [PATCH] modified working code for the recent changes --- .vscode/settings.json | 3 +- Armory/AnalysisLib.h | 39 ++- Armory/ClassCorrParas.h | 160 +++++++++++++ Armory/ClassDetGeo.h | 12 +- Armory/ClassReactionParas_tobeKill.h | 158 ++++++++++++ Armory/GeneralSort.h | 43 ---- Armory/Parameters.h | 306 ------------------------ Cleopatra/ClassTransfer.h | 8 +- Cleopatra/Simulation_Helper.C | 9 +- Cleopatra/alpha.C | 2 +- Cleopatra/{knockout.C => knockoutSim.C} | 2 +- working/{DWBA2 => DWBA} | 0 working/Monitor.C | 51 ++-- working/Monitor.h | 29 ++- 14 files changed, 411 insertions(+), 411 deletions(-) create mode 100644 Armory/ClassCorrParas.h create mode 100644 Armory/ClassReactionParas_tobeKill.h delete mode 100644 Armory/Parameters.h rename Cleopatra/{knockout.C => knockoutSim.C} (99%) rename working/{DWBA2 => DWBA} (100%) diff --git a/.vscode/settings.json b/.vscode/settings.json index 415af4f..162a9a6 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -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, diff --git a/Armory/AnalysisLib.h b/Armory/AnalysisLib.h index 4ff335c..c689b0c 100644 --- a/Armory/AnalysisLib.h +++ b/Armory/AnalysisLib.h @@ -12,9 +12,47 @@ #include #include #include +#include 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 SplitStr(std::string tempLine, std::string splitter, int shift = 0){ std::vector output; @@ -49,7 +87,6 @@ std::vector SplitStr(std::string tempLine, std::string splitter, in return output; }; - //************************************** TCutG TObjArray * LoadListOfTCut(TString fileName, TString cutName = "cutList"){ diff --git a/Armory/ClassCorrParas.h b/Armory/ClassCorrParas.h new file mode 100644 index 0000000..b5c40de --- /dev/null +++ b/Armory/ClassCorrParas.h @@ -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 xnCorr; //correction of xn to match xf + std::vector xScale; // correction of x to be (0,1) + std::vector> xfxneCorr; //correction of xn and xf to match e + std::vector> eCorr; // correction to e, ch -> MeV + std::vector> 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 \ No newline at end of file diff --git a/Armory/ClassDetGeo.h b/Armory/ClassDetGeo.h index 1bae2a1..a50c76f 100644 --- a/Armory/ClassDetGeo.h +++ b/Armory/ClassDetGeo.h @@ -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; diff --git a/Armory/ClassReactionParas_tobeKill.h b/Armory/ClassReactionParas_tobeKill.h new file mode 100644 index 0000000..7abdf74 --- /dev/null +++ b/Armory/ClassReactionParas_tobeKill.h @@ -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 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 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 \ No newline at end of file diff --git a/Armory/GeneralSort.h b/Armory/GeneralSort.h index 632ed7f..b708e94 100644 --- a/Armory/GeneralSort.h +++ b/Armory/GeneralSort.h @@ -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 diff --git a/Armory/Parameters.h b/Armory/Parameters.h deleted file mode 100644 index 660d63b..0000000 --- a/Armory/Parameters.h +++ /dev/null @@ -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 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 xnCorr; //correction of xn to match xf -std::vector xScale; // correction of x to be (0,1) -std::vector> xfxneCorr; //correction of xn and xf to match e -std::vector> eCorr; // correction to e, ch -> MeV -std::vector> 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 \ No newline at end of file diff --git a/Cleopatra/ClassTransfer.h b/Cleopatra/ClassTransfer.h index 6ae8aba..8bf4de8 100644 --- a/Cleopatra/ClassTransfer.h +++ b/Cleopatra/ClassTransfer.h @@ -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; } diff --git a/Cleopatra/Simulation_Helper.C b/Cleopatra/Simulation_Helper.C index 41b20b7..e14204e 100644 --- a/Cleopatra/Simulation_Helper.C +++ b/Cleopatra/Simulation_Helper.C @@ -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 @@ -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; } diff --git a/Cleopatra/alpha.C b/Cleopatra/alpha.C index 623960c..cd68622 100644 --- a/Cleopatra/alpha.C +++ b/Cleopatra/alpha.C @@ -1,4 +1,4 @@ -#include "HELIOS_LIB.h" +#include "ClassHelios.h" #include "TROOT.h" #include "TBenchmark.h" #include "TLorentzVector.h" diff --git a/Cleopatra/knockout.C b/Cleopatra/knockoutSim.C similarity index 99% rename from Cleopatra/knockout.C rename to Cleopatra/knockoutSim.C index ceb0d38..4c05593 100644 --- a/Cleopatra/knockout.C +++ b/Cleopatra/knockoutSim.C @@ -1,4 +1,4 @@ -#include "HELIOS_LIB.h" +#include "ClassHelio.h" #include "TROOT.h" #include "TBenchmark.h" #include "TLorentzVector.h" diff --git a/working/DWBA2 b/working/DWBA similarity index 100% rename from working/DWBA2 rename to working/DWBA diff --git a/working/Monitor.C b/working/Monitor.C index a7d4b62..0393219 100644 --- a/working/Monitor.C +++ b/working/Monitor.C @@ -21,7 +21,7 @@ #include #include #include -#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 ExThetaCM = AnalysisLib::CalExTheta(eCal[detID], x[detID]); - Ex = ExThetaCM[0]; - thetaCM = ExThetaCM[1]; - - }else{ - Ex = TMath::QuietNaN(); - thetaCM = TMath::QuietNaN(); - } + std::pair ExThetaCM = transfer->CalExThetaCM(eCal[detID], x[detID], detGeo->Bfield, detGeo->array1.detPerpDist); + double Ex = ExThetaCM.first; + double thetaCM = ExThetaCM.second; if( thetaCM > thetaCMGate ) { diff --git a/working/Monitor.h b/working/Monitor.h index 4708ae4..fd7df6e 100644 --- a/working/Monitor.h +++ b/working/Monitor.h @@ -12,7 +12,11 @@ #include #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);