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 <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"){
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<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
\ 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<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
\ 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<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
\ 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 <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;
 
 }
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 <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 ) {
 
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 <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);