SimTransfer plot result by defautl

.vscode/c_cpp_properties.json

@@ -29,7 +29,7 @@
             "name": "WinLinux",
             "includePath": [
                 "${workspaceFolder}/**",
-                "/home/ryan/Downloads/root/include/**"
+                "/home/ryan/root_v6.30.06/include/**"
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",

.vscode/settings.json

@@ -1,139 +1,13 @@
 {
   "files.associations": {
-    "script.C": "cpp",
-    "SOLARIS_Qt6_DAQ.pro": "makefile",
-    "qlineseries": "cpp",
-    "cctype": "cpp",
-    "clocale": "cpp",
-    "cmath": "cpp",
-    "cstdarg": "cpp",
-    "cstddef": "cpp",
-    "cstdio": "cpp",
-    "cstdlib": "cpp",
-    "cstring": "cpp",
-    "ctime": "cpp",
-    "cwchar": "cpp",
-    "cwctype": "cpp",
-    "array": "cpp",
-    "atomic": "cpp",
-    "bit": "cpp",
-    "*.tcc": "cpp",
-    "chrono": "cpp",
-    "codecvt": "cpp",
-    "compare": "cpp",
-    "concepts": "cpp",
-    "condition_variable": "cpp",
-    "cstdint": "cpp",
-    "deque": "cpp",
-    "list": "cpp",
-    "map": "cpp",
-    "string": "cpp",
-    "unordered_map": "cpp",
-    "unordered_set": "cpp",
-    "vector": "cpp",
-    "exception": "cpp",
-    "algorithm": "cpp",
-    "functional": "cpp",
-    "iterator": "cpp",
-    "memory": "cpp",
-    "memory_resource": "cpp",
-    "numeric": "cpp",
-    "optional": "cpp",
-    "random": "cpp",
-    "ratio": "cpp",
-    "string_view": "cpp",
-    "system_error": "cpp",
-    "tuple": "cpp",
-    "type_traits": "cpp",
-    "utility": "cpp",
-    "future": "cpp",
-    "initializer_list": "cpp",
-    "iomanip": "cpp",
-    "iosfwd": "cpp",
-    "iostream": "cpp",
-    "istream": "cpp",
-    "limits": "cpp",
-    "mutex": "cpp",
-    "new": "cpp",
-    "numbers": "cpp",
-    "ostream": "cpp",
-    "semaphore": "cpp",
-    "span": "cpp",
-    "sstream": "cpp",
-    "stdexcept": "cpp",
-    "stop_token": "cpp",
-    "streambuf": "cpp",
-    "thread": "cpp",
-    "cinttypes": "cpp",
-    "typeinfo": "cpp",
-    "variant": "cpp",
-    "qdatetime": "cpp",
-    "fstream": "cpp",
-    "allocator": "cpp",
-    "qsignalmapper": "cpp",
-    "GeneralSort.C": "cpp",
-    "ChainMonitors.C": "cpp",
-    "GeneralSortAgent.C": "cpp",
-    "readRawTrace.C": "cpp",
-    "readTrace.C": "cpp",
-    "Cleopatra.C": "cpp",
-    "alpha.C": "cpp",
-    "DWBA_compare.C": "cpp",
-    "DWBARatio.C": "cpp",
-    "ExtractXSec.C": "cpp",
-    "ExtractXSecFromText.C": "cpp",
-    "FindThetaCM.C": "cpp",
-    "InFileCreator.C": "cpp",
-    "IsotopeShort.C": "cpp",
-    "knockout.C": "cpp",
-    "PlotSimulation.C": "cpp",
-    "PlotTGraphTObjArray.C": "cpp",
-    "transfer_test.C": "cpp",
-    "Transfer.C": "cpp",
-    "Simulation_Helper.C": "cpp",
-    "Check_Simulation.C": "cpp",
-    "AutoFit.C": "cpp",
-    "__hash_table": "cpp",
-    "__split_buffer": "cpp",
-    "__tree": "cpp",
-    "bitset": "cpp",
-    "stack": "cpp",
-    "Monitor.C": "cpp",
-    "DataHoSei.C": "cpp",
-    "Monitors.C": "cpp",
-    "__bit_reference": "cpp",
-    "__bits": "cpp",
-    "__config": "cpp",
-    "__debug": "cpp",
-    "__errc": "cpp",
-    "__locale": "cpp",
-    "__mutex_base": "cpp",
-    "__node_handle": "cpp",
-    "__nullptr": "cpp",
-    "__string": "cpp",
-    "__threading_support": "cpp",
-    "__tuple": "cpp",
-    "complex": "cpp",
-    "forward_list": "cpp",
-    "ios": "cpp",
-    "locale": "cpp",
-    "__verbose_abort": "cpp",
-    "Monitors_Util.C": "cpp",
-    "Monitor_Util.C": "cpp",
-    "script_single.C": "cpp",
-    "script_multi.C": "cpp",
-    "Isotope.C": "cpp",
-    "classisotope.h": "c",
-    "knockoutSim.C": "cpp",
-    "test.C": "cpp",
-    "SimTransfer.C": "cpp",
-    "httpaccess.C": "cpp",
-    "httpcontrol.C": "cpp",
-    "SimTransfer2.C": "cpp",
-    "haha.C": "cpp",
-    "SimTransfer_2.C": "cpp"
+    "detectorGeo.txt": "shellscript",
+    "reactionConfig.txt": "shellscript",
+    "DWBA": "shellscript",
+    "SimChecker_Config.txt": "fsharp",
+    "*.C": "cpp"
   },
   
+  "C_Cpp.autoAddFileAssociations": false,
   "better-comments.multilineComments": true,
   "better-comments.tags" : [
     {
@@ -146,10 +20,10 @@
     },
     {
       "tag": "?",
-      "color": "#0076FF",
+      "color": "#1FE0CB",
       "strikethrough": false,
       "backgroundColor": "transparent",
-      "bold": false,
+      "bold": true,
       "italic": false
     },
     {

Cleopatra/Check_Simulation_obsolete.C

@@ -252,7 +252,6 @@ void Check_Simulation(TString filename = "transfer.root",
    
   for( int i = 1; i <= Div[0]*Div[1] ; i++){
     cCheck->cd(i);
-    cCheck->cd(i)->SetGrid();
 
     if( canvas[i-1] == pThetaCM ) {
       cCheck->cd(i)->SetGrid(0,0);
@@ -633,23 +632,24 @@ vector<double> doubleConvertor(vector<TString> arr){
 
 plotID StringToPlotID(TString str){
    
-  if( str == "pEZ") return plotID::pEZ;                      ///0
-  if( str == "pRecoilXY") return plotID::pRecoilXY;                /// 1
-  if( str == "pRecoilXY1" ) return plotID::pRecoilXY1;       /// 2
-  if( str == "pRecoilXY2" ) return plotID::pRecoilXY2;       /// 3
-  if( str == "pRecoilRZ" ) return plotID::pRecoilRZ;         /// 4
-  if( str == "pRecoilRTR" ) return plotID::pRecoilRTR;       /// 5
-  if( str == "pTDiffZ" ) return plotID::pTDiffZ;             /// 6
-  if( str == "pThetaCM" ) return plotID::pThetaCM;           /// 7
-  if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z;       /// 8
-  if( str == "pExCal" ) return plotID::pExCal;               /// 9
-  if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM;   /// 10
-  if( str == "pArrayXY" ) return plotID::pArrayXY;           /// 11
-  if( str == "pInfo" ) return plotID::pInfo;                 /// 12
+  if( str == "pEZ") return plotID::pEZ;                      /// 0
+  if( str == "pRecoilXY") return plotID::pRecoilXY;          /// 1
+  if( str == "pThetaCM" ) return plotID::pThetaCM;           /// 2
+  if( str == "pExCal" ) return plotID::pExCal;               /// 2
+  if( str == "pArrayXY" ) return plotID::pArrayXY;           /// 3
+  if( str == "pInfo" ) return plotID::pInfo;                 /// 4
+  if( str == "pElum1XY" ) return plotID::pElum1XY;           /// 5
+  if( str == "pRecoilXY1" ) return plotID::pRecoilXY1;       /// 6
+  if( str == "pRecoilXY2" ) return plotID::pRecoilXY2;       /// 7
+  if( str == "pTDiffZ" ) return plotID::pTDiffZ;             /// 8
+  if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM;   /// 9
+  if( str == "pRecoilRZ" ) return plotID::pRecoilRZ;         /// 10
+  if( str == "pEElum1R" ) return plotID::pEElum1R;           /// 11
+  if( str == "pRecoilRTR" ) return plotID::pRecoilRTR;       /// 12
+  if( str == "pThetaCM_Z" ) return plotID::pThetaCM_Z;       /// 13
+  if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM;    /// 14
+
   if( str == "pHitID" ) return plotID::pHitID;               /// 13
-  if( str == "pElum1XY" ) return plotID::pElum1XY;           /// 14
-  if( str == "pEElum1R" ) return plotID::pEElum1R;           /// 14
-  if( str == "pElum1RThetaCM" ) return plotID::pElum1RThetaCM;    /// 15
   if( str == "pEmpty" ) return plotID::pEmpty ;              /// 16
   
   return plotID::pEmpty;

Cleopatra/ClassSimPlotter.h

@@ -0,0 +1,501 @@
+#ifndef SimPlotter_h
+#define SimPlotter_h
+
+#include <TTree.h>
+#include <TH2F.h>
+#include <TH1F.h>
+#include <TString.h>
+#include <TCanvas.h>
+#include <TArc.h>
+#include <TLegend.h>
+#include <TLine.h>
+#include <TLatex.h>
+#include "../Armory/ClassDetGeo.h"
+#include "../Armory/ClassReactionConfig.h"
+#include "../Cleopatra/ClassIsotope.h"
+
+enum plotID { pEZ,               /// 0
+              pRecoilXY,         /// 1
+              pRecoilXY1,        /// 2
+              pRecoilXY2,        /// 3
+              pRecoilRZ,         /// 4
+              pRecoilRTR,        /// 5
+              pTDiffZ,           /// 6
+              pThetaCM,          /// 7
+              pThetaCM_Z,        /// 8
+              pExCal,            /// 9
+              pRecoilRThetaCM,   /// 10
+              pArrayXY,          /// 11
+              pInfo,             /// 12
+              pElum1XY,          /// 13
+              pEElum1R,          /// 14
+              pElum1RThetaCM,    /// 15
+              pEmpty };          /// 16
+
+class Plotter{
+public:
+  Plotter(TTree * tree, short reactionID, std::string reactionName, DetGeo detGeo, ExcitedEnergies exList, TString gate, std::vector<plotID> padPlotID);
+  ~Plotter();
+
+  void SetRanges(double zMin, double zMax, double eMax, double elumMax, double thetaCM_Max){
+    this->zRange[0] = zMin;
+    this->zRange[1] = zMax;
+    this->eMax = eMax;
+    this->elumMax = elumMax;
+    this->thetaCM_Max = thetaCM_Max;
+  } 
+  void SetCanvas( int colSize, int rowSize , int padSize, std::vector<plotID> padPlotID ){    
+    cCheck = new TCanvas(Form("SimChecker-%d", rID), 
+                         "Simulation Checker", 
+                         200 * rID, 
+                         200 * rID, 
+                         padSize*colSize, 
+                         padSize*rowSize);
+
+    if(cCheck->GetShowEditor() ) cCheck->ToggleEditor();
+    if(cCheck->GetShowToolBar() )cCheck->ToggleToolBar();
+    cCheck->Divide(colSize,rowSize);
+
+    numPad = colSize * rowSize;
+  }
+
+  void Plot();
+
+private:
+
+  TTree * tree;
+  DetGeo detGeo;
+  std::string reactionName;
+  short rID;
+  ExcitedEnergies exList;
+  int numEx;
+  TString gate;
+  TString gateTrue;
+  std::vector<plotID> padPlotID;
+
+  double zRange[2];
+  double eMax;
+  double recoilOut;
+  double recoilIn;
+  double elumMax;
+  double thetaCM_Max;
+
+  int numPad;
+  TCanvas * cCheck;
+  TH2F * hez;
+  TH2F * hArrayXY;
+  TH2F * heavyXY;
+  TH2F * hElum1XY;
+  TH1F ** hThetaCM;
+  TH1F * hExCal;
+  TH1F * hHit;
+  TH2F * hRecoilXY1;
+  TH2F * hRecoilXY2;
+  TH2F * hTDiffZ;
+  TH2F * hRecoilRThetaCM;
+  TH2F * hRecoilRZ;
+  TH2F * hEElum1Rho;
+  TH2F * hRecoilRTR;
+  TH2F * hElum1RThetaCM;
+  TH2F * hThetaCM_Z;
+  
+  Isotope heavy;
+
+  std::vector<double> FindRange(TString branch);
+
+}; 
+
+Plotter::Plotter(TTree * tree, short reactionID, std::string reactionName, DetGeo detGeo, ExcitedEnergies exList, TString gate, std::vector<plotID> padPlotID){
+  this->tree = tree;
+  this->rID = reactionID;
+  this->reactionName = reactionName;
+  this->exList = exList;
+  numEx = this->exList.ExList.size();
+
+  this->detGeo = detGeo;
+
+  this->recoilOut = detGeo.aux[rID].outerRadius;
+  this->recoilIn = detGeo.aux[rID].innerRadius;
+
+  this->gate = gate;
+  this->gateTrue = gate + Form(" && rID == %d", rID);
+
+  this->padPlotID = padPlotID;
+
+  cCheck = nullptr;
+  hez = nullptr;
+  hArrayXY = nullptr;
+  heavyXY = nullptr;
+  hElum1XY = nullptr;
+  hThetaCM = nullptr;
+  hExCal = nullptr;
+  hHit = nullptr;
+  hRecoilXY1 = nullptr;
+  hRecoilXY2 = nullptr;
+  hTDiffZ = nullptr;
+  hRecoilRThetaCM = nullptr;
+  hRecoilRZ = nullptr;
+  hEElum1Rho = nullptr;
+  hRecoilRTR = nullptr;
+  hElum1RThetaCM = nullptr;
+  hThetaCM_Z = nullptr;
+
+  std::size_t pos1 = reactionName.find(")") + 1;
+  std::size_t pos2 = reactionName.find(" ", pos1);
+  std::string result = reactionName.substr(pos1, pos2 - pos1);
+
+  pos1 = result.find("{") + 1;
+  pos2 = result.find("}", pos1);
+
+  std::string massHeavy = result.substr(pos1, pos2-pos1);
+  std::string symHeavy = result.substr(pos2+1);
+
+  heavy.SetIsoByName(massHeavy + symHeavy);
+
+}
+
+Plotter::~Plotter(){
+
+  if( cCheck ) delete cCheck;
+  if( hez ) delete hez;
+  if( hArrayXY ) delete hArrayXY;
+  if( heavyXY ) delete heavyXY;
+  if( hElum1XY ) delete hElum1XY;
+  if( hExCal ) delete hExCal;
+  if( hRecoilXY1 ) delete hRecoilXY1;
+  if( hRecoilXY2 ) delete hRecoilXY2;
+  if( hTDiffZ ) delete hTDiffZ;
+  if( hRecoilRThetaCM ) delete hRecoilRThetaCM;
+  if( hRecoilRZ ) delete hRecoilRZ;
+  if( hEElum1Rho ) delete hEElum1Rho;
+  if( hRecoilRTR ) delete hRecoilRTR;
+  if( hElum1RThetaCM ) delete hElum1RThetaCM;
+  if( hThetaCM_Z ) delete hThetaCM_Z;
+
+  if(hThetaCM ){
+    for( int i = 0; i < numEx; i++) delete [] hThetaCM[i];
+    delete [] hThetaCM;
+  }
+
+}
+
+inline void Plotter::Plot(){
+
+  for( int i = 1; i <= numPad; i++ ){
+    cCheck->cd(i);
+    plotID  pID = padPlotID[i-1];
+
+    //&======================================= 0 
+    if( pID == pEZ){
+      hez = new TH2F(Form("hez%d",rID), //need unqie name in root.
+                     Form("e-z [gated] @ %5.0f mm; z [mm]; e [MeV]", detGeo.array[rID].firstPos), 
+                     400, zRange[0], zRange[1], 
+                     400, 0, eMax);
+
+      tree->Draw(Form("e:z>>hez%d", rID), gateTrue , "colz");
+    }
+
+    //&======================================= 1
+    if( pID == pRecoilXY  ){
+      heavyXY = new TH2F(Form("heavyXY%d",rID), 
+                           Form("RecoilXY [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[rID].detPos ), 
+                           400, -recoilOut, recoilOut, 
+                           400, -recoilOut, recoilOut);
+      tree->Draw(Form("yRecoil:xRecoil>>heavyXY%d", rID), gateTrue, "colz");
+      
+      // printf("%f, %f\n", recoilOut, recoilIn);
+      TArc * detArc1 = new TArc(0, 0, recoilOut);
+      detArc1->SetLineColor(kBlue-8);
+      detArc1->SetFillStyle(0);
+      detArc1->Draw("same");  
+      TArc * detArc2 = new TArc(0, 0, recoilIn);
+      detArc2->SetLineColor(kBlue-8);
+      detArc2->SetFillStyle(0);
+      detArc2->Draw("same");  
+        
+      // if( reactionConfig.beamX != 0. || reactionConfig.beamY != 0. ){
+      //   TArc * arc = new TArc(reactionConfig.beamX, reactionConfig.beamY, 1);
+      //   arc->SetLineColor(2);
+      //   detArc1->SetFillStyle(0);
+      //   arc->Draw("same");
+      // }
+    }
+
+    //&======================================= 2
+    if( pID == pThetaCM   ){
+      hThetaCM = new TH1F *[numEx]; 
+      TLegend * legend = new TLegend(0.8,0.2,0.99,0.8);
+      double maxCount = 0;
+      for( int h = 0; h < numEx; h++){
+        hThetaCM[h] = new TH1F(Form("hThetaCM%d-%d", rID, h), Form("thetaCM [gated] (ExID=%d); thetaCM [deg]; count", h), 200, 0, thetaCM_Max);
+        hThetaCM[h]->SetLineColor(h+1);
+        hThetaCM[h]->SetFillColor(h+1);
+        hThetaCM[h]->SetFillStyle(3000+h);
+        tree->Draw(Form("thetaCM>>hThetaCM%d-%d", rID, h), gateTrue + Form("&& ExID==%d", h), "");
+        legend->AddEntry(hThetaCM[h], Form("Ex=%5.1f MeV", exList.ExList[h].Ex));
+        double max = hThetaCM[h]->GetMaximum();
+        if( max > maxCount ) maxCount = max;
+      }
+
+      for( int h = 0; h < numEx; h++){
+        hThetaCM[h]->GetYaxis()->SetRangeUser(1, maxCount * 1.2);
+        if( h == 0 ) {
+          hThetaCM[h]->Draw();
+        }else{
+          hThetaCM[h]->Draw("same");
+        }
+      }   
+      legend->Draw();
+    }
+
+    //&======================================= 3
+    if( pID == pExCal  ){
+      double exMin = 9999;
+      double exMax = -1;
+      for( size_t k = 0 ; k < exList.ExList.size(); k++){
+        double kuku = exList.ExList[k].Ex;
+        if( kuku > exMax ) exMax = kuku;
+        if( kuku < exMin ) exMin = kuku;
+      }
+
+      double exPadding = (exMax - exMin) < 1 ? 1 : (exMax - exMin) * 0.3;
+
+      exMin = exMin - exPadding ;
+      exMax = exMax + exPadding ;
+
+      hExCal = new TH1F(Form("hExCal%d",rID), 
+                            Form("calculated Ex [gated]; Ex [MeV]; count / %.2f keV",  (exMax-exMin)/400.*1000),  
+                            400, exMin, exMax);
+      tree->Draw(Form("ExCal>>hExCal%d", rID), gateTrue, "");
+
+
+      double Sn = heavy.CalSp(0,1);
+      double Sp = heavy.CalSp(1,0);
+      double Sa = heavy.CalSp2(4,2);
+      double S2n = heavy.CalSp(0, 2);
+
+      printf("======= Heavy recoil: %s \n", heavy.Name.c_str());
+      printf("Sn  : %6.2f MeV/u \n", Sn);
+      printf("Sp  : %6.2f MeV/u \n", Sp);
+      printf("Sa  : %6.2f MeV/u \n", Sa);
+      printf("S2n : %6.2f MeV/u \n", S2n);
+
+      double yMax = hExCal->GetMaximum();
+      TLatex * text = new TLatex();
+      text->SetTextFont(82);
+      text->SetTextSize(0.06);
+      
+      if(  Sn < exMax ) {TLine * lineSn  = new TLine(Sn,  0,  Sn, yMax); lineSn->SetLineColor(2);  lineSn->Draw("");      text->SetTextColor(2); text->DrawLatex(Sn,  yMax*0.9, "S_{n}");}
+      if(  Sp < exMax ) {TLine * lineSp  = new TLine(Sp,  0,  Sp, yMax); lineSp->SetLineColor(4);  lineSp->Draw("same");  text->SetTextColor(4); text->DrawLatex(Sp,  yMax*0.9, "S_{p}");}
+      if(  Sa < exMax ) {TLine * lineSa  = new TLine(Sa,  0,  Sa, yMax); lineSa->SetLineColor(6);  lineSa->Draw("same");  text->SetTextColor(6); text->DrawLatex(Sa,  yMax*0.9, "S_{a}");}
+      if( S2n < exMax ) {TLine * lineS2n = new TLine(S2n, 0, S2n, yMax); lineS2n->SetLineColor(8); lineS2n->Draw("same"); text->SetTextColor(8); text->DrawLatex(S2n, yMax*0.9, "S_{2n}");}
+      
+    }
+
+    //&======================================= 4
+    if( pID == pArrayXY ){
+
+      hArrayXY = new TH2F(Form("hArrayXY%d", rID), 
+                              "Array-XY [gated]; X [mm]; Y [mm]", 
+                              400, -detGeo.array[rID].detPerpDist*1.5, detGeo.array[rID].detPerpDist*1.5, 
+                              400, -detGeo.array[rID].detPerpDist*1.5, detGeo.array[rID].detPerpDist*1.5);
+      tree->Draw(Form("yArray:xArray>>hArrayXY%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 5
+    if( pID == pInfo ){
+      TLatex text;
+      text.SetNDC();
+      text.SetTextFont(82);
+      text.SetTextSize(0.06);
+      text.SetTextColor(2);
+
+      text.DrawLatex(0., 0.9, reactionName.c_str());
+      text.DrawLatex(0., 0.8, detGeo.Bfield > 0 ? "out of plan" : "into plan");
+      text.SetTextColor(1);
+      text.DrawLatex(0., 0.7, "gate:");
+      
+      text.SetTextColor(2);
+      //check gate text length, if > 30, break by "&&" 
+      int ll = gate.Length();
+      if( ll > 30 ) {
+        std::vector<std::string> strList = AnalysisLib::SplitStr( (std::string) gate.Data(), "&&");
+        for( int i = 0; i < strList.size(); i++){
+          text.DrawLatex(0., 0.6 - 0.05*i, (TString) strList[i]);
+        }
+      }else{
+        text.DrawLatex(0., 0.6, gate);
+      }
+
+      // if( reactionConfig.beamX != 0.0 || reactionConfig.beamY != 0.0 ){
+      //   text.DrawLatex(0.0, 0.1, Form("Bema pos: (%4.1f, %4.1f) mm", reactionConfig.beamX, reactionConfig.beamY));
+      // }
+    }
+
+    //&======================================= 6
+    if( pID == pElum1XY ){
+      hElum1XY = new TH2F(Form("hElum1XY%d", rID), 
+                                Form("Elum-1 XY [gated] @ %.0f mm ; X [mm]; Y [mm]", detGeo.aux[rID].elumPos1),  
+                                400, -elumMax, elumMax, 
+                                400, -elumMax, elumMax);
+      tree->Draw(Form("yElum1:xElum1>>hElum1XY%d", rID), gateTrue, "colz");
+        
+      double count = hElum1XY->GetEntries();        
+      if( count < 2000. ) {
+        hElum1XY->SetMarkerStyle(7);
+        if( count < 500. ) hElum1XY->SetMarkerStyle(3);
+        hElum1XY->Draw("scat");
+      }
+    }
+
+    //&======================================= 7
+    if( pID == pRecoilXY1       ){
+      TH2F * hRecoilXY1 = new TH2F(Form("hRecoilXY1-%d", rID), 
+                                   Form("RecoilXY-1 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[rID].detPos1 ), 
+                                   400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius, 
+                                   400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius);
+      tree->Draw(Form("yRecoil1:xRecoil1>>hRecoilXY1-%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 8
+    if( pID == pRecoilXY2       ){
+      hRecoilXY2 = new TH2F(Form("hRecoilXY2=%d", rID),
+                                   Form("RecoilXY-2 [gated] @ %4.0f mm; X [mm]; Y [mm]", detGeo.aux[rID].detPos2 ), 
+                                  400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius, 
+                                  400, -detGeo.aux[rID].outerRadius, detGeo.aux[rID].outerRadius);
+      tree->Draw(Form("yRecoil2:xRecoil2>>hRecoilXY2-%d", rID), gate, "colz");
+    }
+
+    //&======================================= 9
+    if( pID == pTDiffZ         ){
+      std::vector<double> tDiffRange = FindRange("t-tB");
+      hTDiffZ = new TH2F(Form("hTDiffZ%d", rID), 
+                         "time(Array) - time(Recoil) vs Z [gated]; z [mm]; time diff [ns]", 
+                         400, zRange[0], zRange[1],  
+                         400, tDiffRange[0], tDiffRange[1]);
+      tree->Draw(Form("t - tB : z >> hTDiffZ%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 10
+    if( pID == pRecoilRThetaCM ){
+      hRecoilRThetaCM = new TH2F(Form("hRecoilRThetaCM%d", rID), 
+                                 "RecoilR - thetaCM [gated]; thetaCM [deg]; RecoilR [mm]", 
+                                 400, 0, 60, 
+                                 400, 0, detGeo.aux[rID].outerRadius);
+      tree->Draw(Form("rhoRecoil:thetaCM>>hRecoilRThetaCM%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 11
+    if( pID == pRecoilRZ       ){
+      hRecoilRZ = new TH2F(Form("hRecoilRZ%d", rID), 
+                           "RecoilR - Z [gated]; z [mm]; RecoilR [mm]",  
+                           400, zRange[0], zRange[1], 
+                           400,0, detGeo.aux[rID].outerRadius);
+      tree->Draw(Form("rhoRecoil:z>>hRecoilRZ%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 12
+    if( pID == pEElum1R ){
+      hEElum1Rho = new TH2F(Form("hEElum1Rho%d", rID), 
+                            "Elum-1 E-R [gated]; R[mm]; Energy[MeV]",  
+                            400, 0, elumMax, 
+                            400, 0, eMax);
+      tree->Draw(Form("Tb:rhoElum1>>hEElum1Rho%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 13
+    if( pID == pRecoilRTR      ){
+      std::vector<double> recoilERange =  FindRange("TB");
+      hRecoilRTR = new TH2F(Form("hRecoilRTR%d", rID), 
+                            "RecoilR - recoilE [gated]; recoil Energy [MeV]; RecoilR [mm]", 
+                            500, recoilERange[0], recoilERange[1], 
+                            500, 0, detGeo.aux[rID].outerRadius);
+      tree->Draw(Form("rhoRecoil:TB>>hRecoilRTR%d", rID), gateTrue, "colz");
+    }
+
+    //&======================================= 14
+    if( pID == pThetaCM_Z      ){
+      hThetaCM_Z = new TH2F(Form("hThetaCM_Z%d", rID),
+                            "ThetaCM vs Z ; Z [mm]; thetaCM [deg]",
+                            400, zRange[0], zRange[1], 
+                            400, 0, thetaCM_Max);
+      tree->Draw(Form("thetaCM:z>>hThetaCM_Z%d", rID), gateTrue,"col");
+    }
+
+    //&======================================= 15
+    if( pID == pElum1RThetaCM){
+      int angBin = 400;
+      
+      hElum1RThetaCM = new TH2F(Form("hElum1RThetaCM%d", rID), 
+                                "Elum-1 rho vs ThetaCM [gated]; thatCM [deg]; Elum- rho [mm]", 
+                                angBin, 0, thetaCM_Max,  
+                                400, 0, elumMax);
+      tree->Draw(Form("rhoElum1:thetaCM>>hElum1RThetaCM%d", rID), gateTrue, "colz");   
+      
+      TH1F * htemp = (TH1F *) hElum1RThetaCM->ProjectionX("htemp");  
+      
+      double rel = (thetaCM_Max)*1.0/angBin;
+      printf("angular resolution : %f deg \n", rel);
+      
+      std::vector<double> xList;
+      double old_y = 0;
+      for( int i = 1; i <= angBin; i++){
+        double y = htemp->GetBinContent(i);
+        if( old_y == 0 && y > 0) xList.push_back(htemp->GetBinCenter(i));
+        if( old_y > 0 && y == 0 ) xList.push_back(htemp->GetBinCenter(i));
+        old_y = y;
+      }
+              
+      printf("list of gaps :\n");
+      for( int i = 0; i < (int) xList.size(); i+=2){
+        printf("%d | %.3f - %.3f deg\n", i, xList[i], xList[i+1]);
+      }
+      
+      TF1 f1("f1", "sin(x)");
+      double acceptance = 0;
+      double err1 = 0;
+      double err2 = 0;
+      for( int i = 0; i < (int) xList.size(); i += 2 ){
+        acceptance += f1.Integral(xList[i] * TMath::DegToRad(), xList[i+1] * TMath::DegToRad() ) * TMath::TwoPi();
+        err1 += f1.Integral((xList[i]-rel) * TMath::DegToRad(), (xList[i+1] + rel) * TMath::DegToRad() ) * TMath::TwoPi();
+        err2 += f1.Integral((xList[i]+rel) * TMath::DegToRad(), (xList[i+1] - rel) * TMath::DegToRad() ) * TMath::TwoPi();
+      }
+      printf("acceptance = %f sr +- %f \n", acceptance, (err1-err2)/2);
+      
+      TLatex text;
+      text.SetTextFont(82);
+      text.SetTextSize(0.06);
+      text.SetTextColor(2);
+      text.SetTextAngle(90);
+
+      for( int i = 0; i < (int) xList.size(); i++){
+        text.DrawLatex(xList[i], elumMax/2, Form("%.2f", xList[i]));
+      }
+      
+      text.SetNDC();
+      text.SetTextAngle(0);
+      text.DrawLatex(0.15, 0.15, Form("accp. = %.2f(%.2f) msr", acceptance * 1000., (err1-err2)*1000./2));         
+
+      delete htemp;
+
+    }
+
+  }///===== end of pad loop
+
+}
+
+std::vector<double> Plotter::FindRange(TString branch){
+
+  tree->Draw(Form("%s>>temp1", branch.Data()), gateTrue);
+  TH1F * temp1 = (TH1F *) gROOT->FindObjectAny("temp1");
+
+  std::vector<double> output;
+  output.push_back( temp1->GetXaxis()->GetXmax() );
+  output.push_back( temp1->GetXaxis()->GetXmin() );
+
+  delete temp1;
+  return output;
+}
+
+
+#endif

Cleopatra/ClassTransfer.h

@@ -10,6 +10,7 @@
 #include "TLorentzVector.h"
 #include "TMath.h"
 #include "TF1.h"
+#include "TMacro.h"
 
 //=======================================================
 //#######################################################
@@ -25,6 +26,7 @@ public:
   TransferReaction(){Inititization();};
   TransferReaction(ReactionConfig config, unsigned short ID = 0);
   TransferReaction(std::string configFile, unsigned short ID = 0);
+  TransferReaction(TMacro configMarco, unsigned short ID = 0);
   TransferReaction(int beamA,   int beamZ,
                    int targetA, int targetZ,
                    int recoilA, int recoilZ, float beamEnergy_AMeV);
@@ -38,6 +40,7 @@ public:
   void SetIncidentEnergyAngle(double KEA, double theta, double phi);
 
   void SetReactionFromFile(std::string configFile, unsigned short ID = 0);
+  void SetReactionFromTMacro(TMacro configMacro, unsigned short ID = 0);
   void SetReactionSimple(int beamA, int beamZ,
                          int targetA, int targetZ,
                          int recoilA, int recoilZ, float beamEnergy_AMeV);
@@ -156,6 +159,11 @@ TransferReaction::TransferReaction(std::string configFile, unsigned short ID){
   SetReactionFromFile(configFile, ID);
 }
 
+TransferReaction::TransferReaction(TMacro configMarco, unsigned short ID){
+  Inititization();
+  SetReactionFromTMacro(configMarco, ID);
+}
+
 TransferReaction::TransferReaction(int beamA, int beamZ,
                                    int targetA, int targetZ,
                                    int recoilA, int recoilZ, 
@@ -305,6 +313,36 @@ void TransferReaction::SetReactionFromFile(std::string configFile, unsigned shor
 
 }
 
+void TransferReaction::SetReactionFromTMacro(TMacro configMacro, unsigned short ID){
+
+  if( config.LoadReactionConfig(&configMacro) ){
+
+    SetA(config.beamA, config.beamZ);
+    Seta(config.targetA, config.targetZ);
+
+    SetExA(config.beamEx);
+
+    if( ID > config.recoil.size() ){
+      printf("Reaction Config only has %zu recoil settings. input is %u. Abort.\n", config.recoil.size(), ID);
+      return;
+    }
+
+    recoil = config.recoil[ID];
+    exList = config.exList[ID];
+
+    Setb(recoil.lightA, recoil.lightZ);
+    SetB(recoil.heavyA, recoil.heavyZ);
+    SetIncidentEnergyAngle(config.beamEnergy, 0, 0);
+
+    CalReactionConstant();
+
+  }else{     
+    printf("cannot read TMacro %s.\n", configMacro.GetName());
+    isReady = false;
+  }
+
+}
+
 TString TransferReaction::GetReactionName() const{
   TString rName;
   rName.Form("%s(%s,%s)%s", nameA.c_str(), namea.c_str(), nameb.c_str(), nameB.c_str()); 

Cleopatra/SimChecker.C

@@ -0,0 +1,235 @@
+#include <TFile.h>
+#include <TTree.h>
+#include <TCanvas.h>
+#include <TROOT.h>
+#include <TObjArray.h>
+#include <TStyle.h>
+#include <TH2F.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TArc.h>
+#include <TMath.h>
+#include <TLine.h>
+#include <TSpectrum.h>
+#include <TGraph.h>
+#include <TLegend.h>
+#include <TLatex.h>
+#include <TMacro.h>
+#include <TObjArray.h>
+#include <fstream>
+#include <TCutG.h>
+#include "../Armory/AnalysisLib.h"
+#include "../Armory/ClassDetGeo.h"
+#include "../Armory/ClassReactionConfig.h"
+#include "../Cleopatra/ClassIsotope.h"
+#include "../Cleopatra/ClassTransfer.h"
+#include "../Cleopatra/ClassSimPlotter.h"
+
+plotID StringToPlotID(TString str);
+
+void SimChecker(TString filename = "transfer.root",
+                TString configFile = "../working/SimChecker_Config.txt",
+                Int_t padSize = 500,
+                bool outputCanvas = false){
+
+  printf("=================================================================\n");
+  printf("====================       Simulate Checker    ==================\n");
+  printf("=================================================================\n");   
+
+  TFile * file = new TFile(filename, "read");
+  TTree * tree = (TTree*) file->Get("tree");
+
+  //*================= Get reactions and Ex
+  TMacro * ListOfReactions = (TMacro *) file->FindObjectAny("ListOfReactions");
+  const short numReact = ListOfReactions->GetListOfLines()->GetEntries();
+
+  printf(">>>>> %d reactions found.\n", numReact);
+
+  std::string reactionList[numReact];
+  for( int i = 0; i < numReact; i++ ){
+    std::string haha = ListOfReactions->GetListOfLines()->At(i)->GetName();
+    std::vector<std::string> kaka = AnalysisLib::SplitStr(haha, "|");
+    reactionList[i]= kaka[1];
+  }
+  
+  ExcitedEnergies reactEx[numReact]; //2-D array [i][j] = i-reaction, j-Ex
+  TMacro * AllExList = (TMacro *) file->FindObjectAny("AllExList");
+
+  TMacro * ExID_ReactID_List = (TMacro *) file->FindObjectAny("ExID_ReactID_List");
+  const short numEx = ExID_ReactID_List->GetListOfLines()->GetEntries()-1;
+  
+  for( int i = 1; i <= numEx; i++){
+    std::string haha = ExID_ReactID_List->GetListOfLines()->At(i)->GetName();
+    std::vector<std::string> kaka = AnalysisLib::SplitStr(haha, " ");
+    
+    std::string dudu = AllExList->GetListOfLines()->At(i)->GetName();
+    std::vector<std::string> dada = AnalysisLib::SplitStr(dudu, " ");
+    
+    short rID = atoi(kaka[1].c_str());
+    reactEx[rID].Add(atof(dada[0].c_str()), 
+                     atof(dada[1].c_str()), 
+                     atof(dada[2].c_str()),
+                     atof(dada[3].c_str()));
+  }
+
+  for( int i = 0; i < numReact; i++ ){
+    printf("=========== %s\n", reactionList[i].c_str());
+    reactEx[i].Print();
+  }
+
+  //*================== detGeoID
+  TMacro * detGeotxt = (TMacro *) file->FindObjectAny("detGeo");
+  DetGeo detGeo(detGeotxt);
+  detGeo.Print(true);
+
+  //*================== Get EZ-curve
+  TObjArray * ezList = (TObjArray *) file->FindObjectAny("EZCurve");
+
+  //*================== Get thetaCM = 0
+  TObjArray * thetaCM0List = (TObjArray *) file->FindObjectAny("thetaCM_Z");
+
+  //^################################################ Find e-range, z-range
+
+  double zRange[numReact][2];
+  double eMax[numReact];
+  
+  int count = 0;
+  for( int i = 0; i < numReact; i++ ){
+    zRange[i][0] = detGeo.array[i].zMin-50;
+    zRange[i][1] = detGeo.array[i].zMax+50;
+
+    eMax[i] = -1;
+    for( size_t j = 0; j < reactEx[i].ExList.size() ; j ++){
+      TGraph * func = (TGraph *) ezList->At(count);
+      double aaa = func->Eval(zRange[i][1]);
+      // printf(" xxxxxxxxxxxx %d, %d | %d %.3f\n", i, j, count, aaa);
+      if( aaa > eMax[i] ) eMax[i] = aaa; 
+      count++;
+    }
+
+    eMax[i] = TMath::Ceil( eMax[i] * 1.1 );
+  }
+
+  // for( int i = 0; i < numReact; i++ ){
+  //   printf(" %d | eMax : %.2f, zRange : %.2f, %.2f \n", i, eMax[i], zRange[i][0], zRange[i][1]);
+  // }
+
+  // //^################################################
+  TMacro * config = new TMacro(configFile);
+  int numLine = config->GetListOfLines()->GetSize();
+
+  TString gate;
+  std::vector<plotID> padPlotID;
+
+  float elumMax = 60;
+  float thetaCMMax = 60; //TODO add thetaCM curves in transfer, so that it can be determinated automatically
+
+  int rowCount = 0;
+  int colCount = 0;
+
+  bool startCanvasConfig = false;
+  bool startGateConfig = false;
+  bool startExtra = false;
+  for( int i = 0; i < numLine ; i++){
+    std::string haha = config->GetListOfLines()->At(i)->GetName();
+    std::vector<std::string> dudu = AnalysisLib::SplitStr(haha, ",");
+    TString lala = haha;
+    lala.Remove(3);
+    if( (lala == "   " || lala == "// " || lala == "//=") && dudu.size() == 0) continue;
+    if( lala == "//#" ) break;
+    if( lala == "//*" ) {
+      startCanvasConfig = true;
+      // rowCount ++;
+      continue;
+    }
+    if( lala == "//^" ) {
+      startCanvasConfig = false;
+      startGateConfig = true;
+      continue;
+    }
+    if( lala == "//@" ) {
+      startGateConfig = false;
+      startExtra = true;
+    }
+
+    if( startCanvasConfig ){
+      rowCount ++;
+      // printf("|%s|\n", haha.c_str());
+      if( dudu.size() > colCount ) colCount = dudu.size();
+
+      for( size_t k = 0; k < dudu.size() ; k++){
+        padPlotID.push_back(StringToPlotID(dudu[k]));
+      }
+    }
+
+    if( startGateConfig ){
+      gate = haha;
+    }
+
+    if( startExtra ){
+      if( dudu[0] == "elum_Max" ) elumMax = atof(dudu[2].c_str());
+      if( dudu[0] == "thetaCM_Max" ) thetaCMMax = atof(dudu[2].c_str());
+    }
+
+  }
+
+  gStyle->SetOptStat("");
+  gStyle->SetStatY(0.9);
+  gStyle->SetStatX(0.9);
+  gStyle->SetStatW(0.4);
+  gStyle->SetStatH(0.2);
+  gStyle->SetLabelSize(0.05, "XY");
+  gStyle->SetTitleFontSize(0.1);
+
+  printf(" Canvas division | col : %d, row : %d \n", colCount, rowCount);
+  count = 0;
+  for( int i = 0; i < rowCount; i++){
+    for( int j = 0; j < colCount; j++){
+      printf("%6d", padPlotID[count]);
+      count++;
+    }
+    printf("\n");
+  }
+
+  printf("Gate : %s \n", gate.Data());
+  printf("elum Max : %.2f mm\n", elumMax);
+  printf("thetaCM Max : %.2f deg\n", thetaCMMax);
+
+  printf("#####################################################\n");
+
+  Plotter * plotter[numReact];
+  
+  for( int i = 0; i < numReact; i++){
+    plotter[i] = new Plotter(tree, i, reactionList[i], detGeo, reactEx[i], gate, padPlotID);
+    plotter[i]->SetRanges(zRange[i][0], zRange[i][1], eMax[i], elumMax, thetaCMMax);
+    plotter[i]->SetCanvas(colCount, rowCount, 500, padPlotID);
+    plotter[i]->Plot();
+  }
+
+  return;
+
+}
+
+
+plotID StringToPlotID(TString str){
+   
+  if( str == "pEZ")              return plotID::pEZ;              /// 0
+  if( str == "pRecoilXY")        return plotID::pRecoilXY;        /// 1
+  if( str == "pThetaCM" )        return plotID::pThetaCM;         /// 2
+  if( str == "pExCal" )          return plotID::pExCal;           /// 3
+  if( str == "pArrayXY" )        return plotID::pArrayXY;         /// 4
+  if( str == "pInfo" )           return plotID::pInfo;            /// 5
+  if( str == "pElum1XY" )        return plotID::pElum1XY;         /// 6
+  if( str == "pRecoilXY1" )      return plotID::pRecoilXY1;       /// 7
+  if( str == "pRecoilXY2" )      return plotID::pRecoilXY2;       /// 8
+  if( str == "pTDiffZ" )         return plotID::pTDiffZ;          /// 9
+  if( str == "pRecoilRThetaCM" ) return plotID::pRecoilRThetaCM;  /// 10
+  if( str == "pRecoilRZ" )       return plotID::pRecoilRZ;        /// 11
+  if( str == "pEElum1R" )        return plotID::pEElum1R;         /// 12
+  if( str == "pRecoilRTR" )      return plotID::pRecoilRTR;       /// 13
+  if( str == "pThetaCM_Z" )      return plotID::pThetaCM_Z;       /// 14
+  if( str == "pElum1RThetaCM" )  return plotID::pElum1RThetaCM;   /// 15
+  if( str == "pEmpty" )          return plotID::pEmpty ;          /// 16
+  
+  return plotID::pEmpty;
+}

Cleopatra/SimTransfer.C

@@ -96,9 +96,13 @@ void Transfer(
   TMacro dwbaExList_Used;
   TMacro dwbaReactList_Used;
 
+  bool useDWBA[numTransfer];
+  for( int i = 0; i < numTransfer; i++ ) useDWBA[i] = false;
+  
   if( distFile->IsOpen() ) {
     printf("\e[32m#################################### Load DWBA input : %s  \e[0m\n", ptolemyRoot.Data());
-    printf("--------- Found DWBA thetaCM distributions. Use the ExList from DWBA.\n"); 
+    printf("--------- Found DWBA thetaCM distributions.\n"); 
+    printf("          Checking DWBA matches withe %s.\n",  basicConfig.c_str()); 
 
     distList = (TObjArray *) distFile->FindObjectAny("thetaCM_TF1"); // the function List
     dwbaExList = (TMacro *) distFile->FindObjectAny("ExList");   
@@ -107,22 +111,23 @@ void Transfer(
 
     int numEx = dwbaExList->GetListOfLines()->GetSize() - 1 ;
     
-    // for( int i = 0; i < numTransfer; i++){ transfer[i].GetExList()->Clear(); }
     ExcitedEnergies dwbaExTemp[numTransfer];
 
     for( int i = 1; i <= numEx ; i++){
-      //Check DWBA reaction is same as transfer setting
       std::string reactionName = dwbaReactList->GetListOfLines()->At(i-1)->GetName();
-
+      printf(" %d | Checking %s from DWBA \n", i, reactionName.c_str());
       for( int j = 0; j < numTransfer; j++){
+        //Check DWBA reaction is same as transfer setting
         if( reactionName.find( transfer[j].GetReactionName().Data() ) != std::string::npos) {
+          printf("    >>> found %s in %s\n", transfer[j].GetReactionName().Data(), basicConfig.c_str());
           std::string temp = dwbaExList->GetListOfLines()->At(i)->GetName();
           dwbaReactList_Used.AddLine((reactionName + " | " + std::to_string(j)).c_str());
           dwbaExList_Used.AddLine(temp.c_str());
           if( temp[0] == '/' ) continue;
           std::vector<std::string> tempStr = AnalysisLib::SplitStr(temp, " ");
-          // transfer[j].GetExList()->Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
           dwbaExTemp[j].Add( atof(tempStr[0].c_str()), atof(tempStr[1].c_str()), 1.0, 0.00);
+        }else{
+          printf("    XXX Not found\n");
         }
       }
     }
@@ -130,11 +135,13 @@ void Transfer(
     for( int i = 0; i < numTransfer; i++ ){
       if( dwbaExTemp[i].ExList.size() > 0 ) {
         transfer[i].GetExList()->Clear();
-        for( size_t j = 0 ; dwbaExTemp[i].ExList.size(); j ++ ){
+        for( size_t j = 0 ; j < dwbaExTemp[i].ExList.size(); j ++ ){
           transfer[i].GetExList()->Add( dwbaExTemp[i].ExList[j].Ex, dwbaExTemp[i].ExList[j].xsec, 1.0, 0.00);
         }
+        useDWBA[i] = true;
       }else{
         printf("Cannot match %s with DWBA, use Reaction Ex List\n", transfer[i].GetReactionName().Data());
+        useDWBA[i] = false;
       }
     }
 
@@ -176,6 +183,20 @@ void Transfer(
     dwbaReactList_Used.Write("DWBA_ReactionList", 1);
   }
 
+  TMacro allExList;
+  allExList.AddLine("#---Ex   relative_xsec   SF   sigma_in_MeV");
+  TMacro exIDReactIDList; //list of all ex and corresponding Reaction ID
+  exIDReactIDList.AddLine("#-- ExID  ReactionID");
+  for( int i = 0; i < numTransfer; i++){
+    std::vector<EnergyLevel> tempExList = transfer[i].GetExList()->ExList;
+    for( size_t j = 0; j < tempExList.size(); j++){
+      allExList.AddLine(Form("%9.5f     %9.5f  %3.1f  %5.3f", tempExList[j].Ex, tempExList[j].xsec, tempExList[j].SF, tempExList[j].sigma));
+      exIDReactIDList.AddLine(Form("%ld   %d", j, i));
+    }
+  }
+
+  allExList.Write("AllExList");
+  exIDReactIDList.Write("ExID_ReactID_List");
 
   TMacro hitMeaning;
   hitMeaning.AddLine("======================= meaning of Hit\n"); 
@@ -366,7 +387,7 @@ void Transfer(
   delete [] gx;
   delete gList;
 
-  //--- cal modified f
+  //--- cal E-Z curve with finite detector correction
   int numEx = 0;
   for( int i = 0; i < numTransfer; i++){
     if( !listOfTransfer[i] ) continue;
@@ -520,7 +541,7 @@ void Transfer(
       // }
 
       //==== Calculate thetaCM, phiCM
-      if( distFile->IsOpen()){
+      if( distFile->IsOpen() && useDWBA[rID] ){
         angDist = (TF1 *) distList->At(ExID);
         thetaCM = angDist->GetRandom() / 180. * TMath::Pi();
       }else{
@@ -737,13 +758,14 @@ int main (int argc, char *argv[]) {
   printf("==========     Simulate Transfer reaction in HELIOS    ==========\n");
   printf("=================================================================\n");
   
-  if(argc == 2 || argc > 5 ) { 
-    printf("Usage: ./Transfer [1] [2] [3] [4]\n");
+  if(argc == 2 || argc > 6 ) { 
+    printf("Usage: ./Transfer [1] [2] [3] [4] [5]\n");
     printf("       default file name \n");
     printf("   [1] reactionConfig.txt (input) reaction Setting \n");
     printf("   [2] detectorGeo.txt    (input) detector Setting \n");
     printf("   [3] DWBA.root          (input) thetaCM distribution from DWBA \n");  
     printf("   [4] transfer.root      (output) rootFile name for output \n");
+    printf("   [5] 1                  (input) 0 = no plot, 1 = plot \n");
 
     printf("-----------------------------------------------------------------\n");
     printf(" When DWBA.root provided.\n");
@@ -757,29 +779,32 @@ int main (int argc, char *argv[]) {
   std::string   detGeoFile = "detectorGeo.txt";
   TString      ptolemyRoot = "DWBA.root"; // when no file, use isotropic distribution of thetaCM
   TString     saveFileName = "transfer.root";
-  bool           isPlot = false;
+  bool           isPlot = true;
    
   if( argc >= 2) basicConfig = argv[1];  
   if( argc >= 3) detGeoFile = argv[2];  
   if( argc >= 4) ptolemyRoot = argv[3];
   if( argc >= 5) saveFileName = argv[4];
+  if( argc >= 6) isPlot = atoi(argv[5]);
   
   Transfer( basicConfig, detGeoFile, ptolemyRoot, saveFileName);
 
-  //run Armory/Check_Simulation
-  // if( isPlot ){
-  //   std::ifstream file_in;
-  //   file_in.open("../Cleopatra/Check_Simulation.C", std::ios::in);
-  //   if( file_in){
-  //     printf("---- running ../Cleopatra/Check_Simulation.C on %s \n", saveFileName.Data());
-  //     TString cmd;
-  //     cmd.Form("root -l '../Cleopatra/Check_Simulation.C(\"%s\")'", saveFileName.Data());
+  //run Cleopatra/SimChecker.C
+  if( isPlot ){
+    std::ifstream file_in;
+    file_in.open("../Cleopatra/SimChecker.C", std::ios::in);
+    if( file_in){
+      printf("---- running ../Cleopatra/SimChecker.C on %s \n", saveFileName.Data());
+      TString cmd;
+      cmd.Form("root -l '../Cleopatra/SimChecker.C(\"%s\")'", saveFileName.Data());
       
-  //     system(cmd.Data());
-  //   }else{
-  //     printf("cannot find ../Cleopatra/Check_Simulation.C \n");
-  //   }
-  // }
+      system(cmd.Data());
+    }else{
+      printf("cannot find ../Cleopatra/SimChecker.C \n");
+    }
+  }
+
+  return 0;
   
 }
 

Cleopatra/makefile

@@ -3,7 +3,7 @@ CFLAG = -O2
 
 depend = ClassTransfer.h ClassHelios.h ClassIsotope.h ClassDecay.h constant.h potentials.h
 
-ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimTransfer_single SimAlpha
+ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM SimTransfer SimAlpha
 
 all: $(ALL)
 

README.md

@@ -14,10 +14,36 @@ Analysis
 ├── root_data // symbolic link to converted root file, created by SetUpNewExp  
 └── working // working directory, depends on experiment.
 
+# SOLARIS.sh
+
+this batch shell script adds few enviroment variables and functions. Add the Armory and Cleopatra into the system PATH.
+
+```sh
+>source SOLARIS.sh
+```
+
+# Event Builder
+
+Please download the SOLARIS_DAQ, under the Aux directory, make, and link the EventBuilder to Armory.
+
+The reason for having EventBuilder in the DAQ code is the Hit.h is original from the DAQ code.
+
 # ROOT issue
 
 We are still using TProof for parallel calculation. TProof is not pre-compiled since 6.32+. And 6.30 only precompiled for Ubuntu 22.04. So, for system using Ubuntu 24.04, user must precompiled to root in order to work.
 
+## Compilation
+
+We can manual download the git repository of the root following the instruction. in the cmake 
+
+```sh
+root_install="/opt/root_v6.32.00_compiled"
+root_src="root_src"
+root_build="root_build"
+cmake -DCMAKE_INSTALL_PREFIX=$root_install $root_src -Dproof=ON -Dmathmore=ON
+sudo cmake --build . --target install -j <number_of_threads>
+```
+
 # Analysis & Simulation
 
 The Armory/AnalysisLib.h constains many small but handy functions.
@@ -26,28 +52,30 @@ All class headers are started with Class*.h
 
 The classes **DetGeo**** and **ReactionConfig** are fundamental for loading the detectorGeo.txt and reactionConfig.txt. 
 
-Both txt file support empty lines, and up to 2 settings. The reason for that is for dual-array configuration. It has potentail to extend and include more settings. But it is two now, one for upstream array (reaction) and downstream array (reaction).
+Both txt file support empty lines, can have multiple settings. The reason for that is for many-array configuration. 
 
-The **TransferReaction** class is only use one of the reaction from the reactionConfig.txt. 
+The **TransferReaction** class is only use one of the reaction from the reactionConfig.txt. This class generates the TLorentzVector for ligh and heavy recoils.
 
 ```C++
   TransferReaction::SetReactionFromFile("reactionConfig.txt", ID); // ID = 0 or 1
 ```
-Same for the **Helios** class
+Same for the **Helios** class, **Helios** class use the detectorGeo.txt. It takes TLorentzVector and calculate does it be detected by the array or recoil detector.
 
 ```C++
   HELIOS::SetDetectorGeometry("detectorGeo.txt", ID); // ID = 0 or 1
 ```
 
+## Simulation
 
-# Event Builder
+Simply run
+```sh
+>SimTransfer
+```
 
-The EventBuilder is at the armory. It depends on the Hit.h and SolReader.h.
+it will digest the detectorGeo.txt, reactionConfig.txt, if DWBA.root exist, find the reactions. 
 
-## Hit.h
+* it does not have TargetScattering (yet)
+* for multiple reactions, it will randomly use any and disregard the total Xsec of different reactions. The Xsec only takes effect within same reaction.
+* the decay of heavy recoil only have isotropic decay.
 
-The Hit class stores a hit (or a data block)
 
-## SolReader.h
-
-The SolReader class read the sol file. It can be loaded in CERN ROOT alone. 

working/Check_Simulation_Config.txt

@@ -1,30 +0,0 @@
-/////enum plotID { pEZ,               /// 0
-/////              pRecoilXY,         /// 1
-/////              pRecoilXY1,        /// 2
-/////              pRecoilXY2,        /// 3
-/////              pRecoilRZ,         /// 4
-/////              pRecoilRTR,        /// 5
-/////              pTDiffZ,           /// 6
-/////              pThetaCM,          /// 7
-/////              pThetaCM_Z,        /// 8
-/////              pExCal,            /// 9
-/////              pRecoilRThetaCM,   /// 10
-/////              pArrayXY,          /// 11
-/////              pInfo,             /// 12
-/////              pHitID,            /// 13
-/////              pElum1XY,          /// 14
-/////              pEElum1R,          /// 15
-/////              pElum1RThetaCM,    /// 16
-/////              pEmpty };          /// 17
-/////=============================================== User Config
-{pEZ, pExCal, pThetaCM, pRecoilRZ, break, pThetaCM_Z, pRecoilXY, pInfo, pArrayXY}   //Canvas config
-hit == 1 && loop <= 1 && thetaCM > 10
-60              //elum range
-{0,50}        //thetaCM range
-true          //shownKELines
-false          //isOverRideEx
-{-0.5, 4.0}    // over-rdied Ex range
-///============================== example of gate
-hit == 1 && loop <= 1 && thetaCM > 10 && detRowID == 0
-hit == 1 && loop <= 1
-15 < rhoElum1 && rhoElum1 < 50 && rhoElum2 > 60

working/detectorGeo.txt

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