fix alpha simulation and FindThetaCM

.gitignore

@@ -24,4 +24,5 @@ Cleopatra/Isotope
 Cleopatra/IsotopeShort
 Cleopatra/PlotSimulation
 Cleopatra/PlotTGraphTObjArray
+Cleopatra/SimAlpha
 Cleopatra/Transfer

Cleopatra/ClassHelios.h

@@ -77,6 +77,7 @@ class HELIOS{
 public:
 
   HELIOS();
+  HELIOS(std::string filename, unsigned short ID);
   ~HELIOS();
   
   void SetCoincidentWithRecoil(bool TorF){ this->isCoincidentWithRecoil = TorF;}
@@ -166,10 +167,24 @@ private:
   
   const double c = 299.792458; //mm/ns
 
+  void Clear();
 
 };
 
 HELIOS::HELIOS(){
+  Clear();
+}
+
+HELIOS::HELIOS(std::string filename, unsigned short ID){
+  Clear();
+  SetDetectorGeometry(filename, ID);
+}
+
+HELIOS::~HELIOS(){
+  
+}
+
+void HELIOS::Clear(){
 
   orbitb.Clear();
   orbitB.Clear();
@@ -191,11 +206,6 @@ HELIOS::HELIOS(){
   overrideDetDistance = false;
   overrideFirstPos = false;
   isCoincidentWithRecoil = false;
-   
-}
-
-HELIOS::~HELIOS(){
-  
 }
 
 void HELIOS::OverrideMagneticField(double BField){ 

Cleopatra/ClassTransfer.h

@@ -43,7 +43,7 @@ public:
   void SetExA(double Ex);
   void SetExB(double Ex);
 
-  TString GetReactionName();
+  TString GetReactionName() const;
   TString GetReactionName_Latex();
 
   ReactionConfig  GetRectionConfig() { return config;}
@@ -65,7 +65,7 @@ public:
   TLorentzVector GetPB() const {return PB;}
   
   void PrintFourVectors() const;
-  void PrintReaction() const;
+  void PrintReaction(bool withEx = true) const;
 
   double CalkCM(double ExB); //momentum at CM frame
   void CalReactionConstant();
@@ -252,7 +252,7 @@ void TransferReaction::SetReactionFromFile(string configFile, unsigned short ID)
 
 }
 
-TString TransferReaction::GetReactionName(){
+TString TransferReaction::GetReactionName() const{
   TString rName;
   rName.Form("%s(%s,%s)%s", nameA.c_str(), namea.c_str(), nameb.c_str(), nameB.c_str()); 
   return rName;
@@ -272,6 +272,7 @@ TString TransferReaction::format(TString name){
   }
   return "^{"+temp2+"}"+temp;
 }
+
 TString TransferReaction::GetReactionName_Latex(){
   TString rName;
   rName.Form("%s(%s,%s)%s @ %.2f MeV/u", format(nameA).Data(), format(namea).Data(), format(nameb).Data(), format(nameB).Data(), config.beamEnergy); 
@@ -327,8 +328,10 @@ void TransferReaction::PrintFourVectors() const {
 
 }
 
-void TransferReaction::PrintReaction() const {
+void TransferReaction::PrintReaction(bool withEx) const {
 
+  printf("=====================================================\n");  
+  printf("\e[1m\e[33m %s \e[0m\n", GetReactionName().Data());
   printf("=====================================================\n");  
   printf("------------------------------ Beam\n");
   printf("   beam : A = %3d, Z = %2d, Ex = %.2f MeV\n", config.beamA, config.beamZ, config.beamEx);
@@ -343,9 +346,11 @@ void TransferReaction::PrintReaction() const {
   printf("  light : A = %3d, Z = %2d \n", recoil.lightA, recoil.lightZ);
   printf("  heavy : A = %3d, Z = %2d \n", recoil.heavyA, recoil.heavyZ);
   printf("=====================================================\n"); 
+
+  if( withEx ) { 
     exList.Print();
     printf("=====================================================\n");  
-
+  }
 }
 
 void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){
@@ -395,7 +400,6 @@ void TransferReaction::Event(double thetaCM_rad, double phiCM_rad){
    
 }
 
-
 std::pair<double, double> TransferReaction::CalExThetaCM(double e, double z, double Bfield, double perpDist){
 
   double Ex = TMath::QuietNaN();

Cleopatra/FindThetaCM.C

@@ -22,11 +22,12 @@ int main(int argc, char *argv[]){
   printf("===     Find ThetaCM convrage for each detector at Ex        ====\n");
   printf("=================================================================\n");
 
-  if(argc < 2 || argc > 6) { 
+  if(argc < 2 || argc > 7) { 
     printf("Usage: ./FindThetaCM Ex\n");
     printf("Usage: ./FindThetaCM Ex nDiv\n");
     printf("Usage: ./FindThetaCM Ex nDiv X-Ratio\n");
-    printf("Usage: ./FindThetaCM Ex nDiv X-Ratio reactionTxt detGeoTxt\n");
+    printf("Usage: ./FindThetaCM Ex nDiv X-Ratio reactionTxt detGeoTxt ID\n");
+    printf("   * default is the first settings from reaction and detGeo.\n");
     exit(0);
   }  
   
@@ -35,24 +36,21 @@ int main(int argc, char *argv[]){
   int nDiv = 1;
   string reactionTxt = "reactionConfig.txt";
   string detGeoTxt = "detectorGeo.txt";
+  int ID = 0;
 
-  if ( argc >= 2 ){
+  if ( argc >= 2 ) Ex = atof(argv[1]);
-    Ex = atof(argv[1]);
+  if ( argc >= 3 ) nDiv = atoi(argv[2]);
-  }
+  if ( argc >= 4 ) xRatio = atof(argv[3]);
-  if ( argc >= 3 ){
+  if ( argc >= 5 ) reactionTxt = argv[4];
-    nDiv = atoi(argv[2]);
+  if ( argc >= 6 ) detGeoTxt = argv[5];
-  }
+  if ( argc >= 6 ) ID = atoi(argv[6]);
-  if ( argc >= 4 ){
+
-    xRatio = atof(argv[3]);
+  if( nDiv < 1 ) {
-  }
+    printf(" nDiv must be >= 1 \n");
-  if ( argc >= 5 ){
+    return -1;
-    reactionTxt = argv[4];
-  }
-  if ( argc >= 6 ){
-    detGeoTxt = argv[5];
   }
   
-  FindThetaCM(Ex, nDiv, xRatio, reactionTxt, detGeoTxt);
+  FindThetaCM(Ex, nDiv, xRatio, reactionTxt, detGeoTxt, ID);
   
   return 0;
 }

Cleopatra/FindThetaCM.h

@@ -20,122 +20,41 @@
 #include "../Cleopatra/ClassTransfer.h"
 
 void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95, 
-            std::string basicConfig="reactionConfig.txt",  
+            std::string reactionConfigFileName="reactionConfig.txt",  
             std::string detGeoFileName = "detectorGeo.txt", unsigned short ID = 0){
 
-  //---- reaction
-  int AA, zA; //beam
-  int Aa, za; //target
-  int Ab, zb; //recoil-1
-  double ExA; 
 
-  //---- beam
+  ///========================================================= load files
-  double KEAmean, KEAsigma; // MeV/u , assume Guassian
-  double thetaMean, thetaSigma; // mrad , assume Guassian due to small angle
-  double xBeam, yBeam; // mm
 
-   /**///========================================================= load files
+  TransferReaction reaction(reactionConfigFileName, ID);
-  ReactionConfig reConfig;
-  DetGeo detGeo;
-  if( reConfig.LoadReactionConfig(basicConfig) ){
-
-    KEAmean = reConfig.beamEnergy;
-    KEAsigma = reConfig.beamEnergySigma;
-
-    thetaMean = reConfig.beamTheta;
-    thetaSigma = reConfig.beamThetaSigma;
-
-    xBeam = reConfig.beamX;
-    yBeam = reConfig.beamY;
-
-    AA = reConfig.beamA;        zA = reConfig.beamZ;
-    Aa = reConfig.targetA;      za = reConfig.targetZ;
-    Ab = reConfig.recoil[ID].lightA; zb = reConfig.recoil[ID].lightZ;
-
-    ExA = reConfig.beamEx;
-    
-  }else{
-    printf("cannot load %s \n", basicConfig.c_str());
-    return;
-  }
-  
-   std::vector<double> pos;
-   double a = 11.5;
-   double length = 50.5;
-   double firstPos = 0;
-   int iDet = 6;
-   int jDet = 4;
-   double BField = 0;
-   
-   //=============================================================
-   //=============================================================
-   //=============================================================
-   //===== Set Reaction
-   TransferReaction reaction;
-   int AB = AA+Aa-Ab, zB = zA+za-zb;
-   reaction.SetA(AA,zA);
-   reaction.Seta(Aa,za);
-   reaction.Setb(Ab,zb);
-   reaction.SetB(AB,zB);
-   reaction.SetIncidentEnergyAngle(KEAmean, 0, 0);
   reaction.SetExB(Ex);  
-   reaction.SetExA(ExA);
   reaction.CalReactionConstant();
+  reaction.PrintReaction(false);
+  ReactionConfig reConfig = reaction.GetRectionConfig();
+  Recoil recoil = reaction.GetRecoil();
 
-   printf("===================================================\n");
+  HELIOS helios(detGeoFileName, ID);
-   printf("=========== %27s ===========\n", reaction.GetReactionName().Data());
+  helios.PrintGeometry();
-   printf("===================================================\n");
+  DetGeo detGeo = helios.GetDetectorGeometry();
-   printf("----- loading reaction from : %s. \n", basicConfig.c_str());
+  Array array = helios.GetArrayGeometry();
-   printf("       Ex A: %7.3f MeV\n", ExA);
-   printf("         KE: %7.4f \n", KEAmean);
-   printf("      theta: %7.4f \n", thetaMean);
-   printf("offset(x,y): %7.4f, %7.4f mm \n", xBeam, yBeam);
-   printf("    Q-value: %7.4f MeV \n", reaction.GetQValue() );
-   printf("     Max Ex: %7.4f MeV \n", reaction.GetMaxExB() );
-   printf("===================================================\n");
-   
-   
-   printf("----- loading detector geometery : %s.", detGeoFileName.c_str());
-   if(detGeo.LoadDetectorGeo(detGeoFileName) ){
-     pos = detGeo.array[ID].detPos;
-     a = detGeo.array[ID].detPerpDist;
-     length = detGeo.array[ID].detLength;
-     firstPos = detGeo.array[ID].firstPos;
-     iDet = detGeo.array[ID].nDet;
-     jDet = detGeo.array[ID].mDet;
-     BField = detGeo.Bfield;
-     
-     printf("... done.\n");
-     
-   }else{
-     printf("... fail\n");
-     return;
-   }
-    
-   
-   vector<double> midPos;
-   
-   for(int i = 0; i < iDet; i++){
-      if( firstPos > 0 ){
-         midPos.push_back(pos[i]+length/2.);
-      }else{
-         midPos.push_back(pos[i]-length/2.);
-      }
-   }
    
   //calculate a TGraph for thetaCM vs z
-   double px[100];
+  const int nData = 170;
-   double py[100];
+  double px[nData];
+  double py[nData];
   
   double mb = reaction.GetMass_b();
   double kCM = reaction.GetMomentumbCM();
   double q = TMath::Sqrt(mb*mb + kCM * kCM );
   double beta = reaction.GetReactionBeta() ;
-   double slope = 299.792458 * zb * abs(BField) / TMath::TwoPi() * beta / 1000.; // MeV/mm
+  double BField = detGeo.Bfield;
+  double slope = reaction.GetEZSlope(BField);
   double gamma = reaction.GetReactionGamma();
-   for(int i = 0; i < 100; i++){
+  double perpDist = array.detPerpDist;
+
+  for(int i = 0; i < nData; i++){
     double thetacm = (i + 5.) * TMath::DegToRad();
-      double temp = TMath::TwoPi() * slope / beta / kCM * a / TMath::Sin(thetacm); 
+    double temp = TMath::TwoPi() * slope / beta / kCM * perpDist / TMath::Sin(thetacm); 
     px[i] = beta /slope * (gamma * beta * q - gamma * kCM * TMath::Cos(thetacm)) * (1 - TMath::ASin(temp)/TMath::TwoPi());
     py[i] = thetacm * TMath::RadToDeg();   
   }
@@ -149,8 +68,8 @@ void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95,
   ///xt->Draw("AC*");
   
   /// find the minimum z position and the corresponding theta
-   double zMin0 = 0;
+  double zMin0 = 99999999;
-   double tMin0 = 0;
+  double tMin0 = 99999999;
   for( double ttt = 3 ; ttt < 20 ; ttt += 0.1 ){
     double zzz = xt->Eval(ttt);
     if( zzz < zMin0 ) {
@@ -160,24 +79,41 @@ void FindThetaCM(double Ex, int nDivision=1, double XRATION = 0.95,
   }   
   printf(" z min %f mm at thetaCM %f deg \n", zMin0, tMin0);
 
-
+  TGraph * tx = new TGraph(nData, px, py);
-   TGraph * tx = new TGraph(100, px, py);
   tx->SetName(Form("tx"));
   tx->SetLineColor(4);
 
-   //tx->Draw("AC*");
+  //Remove nan data
+  for( int i = tx->GetN() -1 ;  i >= 0 ; i--){
+    if( TMath::IsNaN(tx->GetPointX(i)) ) tx->RemovePoint(i);
+  }
   
-   /**///========================================================= result
+  // tx->Draw("AC");
+
+  ///========================================================= result
+    
+  int iDet = array.nDet;
+  double length = array.detLength;
+  vector<double> midPos;
+   
+  for(int i = 0; i < iDet; i++){
+    if( array.firstPos > 0 ){
+      midPos.push_back(array.detPos[i]+length/2.);
+    }else{
+      midPos.push_back(array.detPos[i]-length/2.);
+    }
+    // printf("%2d | %f \n", i, midPos.back());
+  }
 
   printf("==== ThetaCM in degree =================\n");
-   printf("========================= x-ratio : %f, number of division : %d \n", XRATION, nDivision);
+  printf(" x-ratio : %f, number of division : %d \n", XRATION, nDivision);
   printf("\n");
   for( int j = 0; j < nDivision + 1; j++) printf("%5.2f   ", -XRATION + 2*XRATION/nDivision*j);
   printf(" <<-- in X \n");
   for( int j = 0; j < nDivision + 1; j++) printf("%5s   ", "  |  ");
   printf("\n");
   for( int j = 0; j < nDivision + 1; j++) printf("%5.2f   ", length/2 -length*XRATION/2 + length*XRATION/nDivision*j);
-   printf(" <<-- in cm \n\n");
+  printf(" <<-- in mm \n\n");
   printf("=========================      Ex : %6.4f MeV\n", Ex);
   printf("            %6s - %6s |  %6s, %6s, %6s\n", "Min",  "Max", "Mean", "Dt", "sin(x)dx * 180/pi");
   printf("-------------------------------------------------\n");

Cleopatra/alpha.C

@@ -1,4 +1,4 @@
-#include "ClassHelios.h"
+#include "../Cleopatra/ClassHelios.h"
 #include "TROOT.h"
 #include "TBenchmark.h"
 #include "TLorentzVector.h"
@@ -18,25 +18,24 @@
 
 const double ma = 3727.3792; // alpha mass
 
-void alpha(){
+void alpha( int numEvent = 100000){
 
   //================================================= User Setting
-   const int numEnergy = 4;
+  std::string heliosDetGeoFile = "detectorGeo.txt";
-   double energy [numEnergy] = {3.18, 5.16, 5.49, 5.81};
+  int geoID = 0;
  
-   int numEvent = 1000000;
+  //std::vector<double> energy = {3.18, 5.16, 5.49, 5.81};
+  std::vector<double> energy = {5.34, 5.42, 5.68, 6.05, 6.23, 6.77, 8.78}; //228Th
  
-   //---- HELIOS detector geometry
+  //---- Over-ride HELIOS detector geometry
-   //string heliosDetGeoFile = "detectorGeo.txt";
+  // double BField = 2.5; // T
-   string heliosDetGeoFile = "";
+  // double BFieldTheta = 0.; // direction of B-field
-   double BField = 2.5; // T
+  // bool isCoincidentWithRecoil = false; 
-   double BFieldTheta = 0.; // direction of B-field
+  // double eSigma = 0.040 ; // detector energy sigma MeV
-   bool isCoincidentWithRecoil = false; 
+  // double zSigma = 0.500 ; // detector position sigma mm
-   double eSigma = 0.040 ; // detector energy sigma MeV
-   double zSigma = 0.500 ; // detector position sigma mm
   
   //---- save root file name
-   TString saveFileName = "alpha.root";
+  TString saveFileName = "SimAlpha.root";
   
   //=============================================================
   //=============================================================
@@ -45,6 +44,7 @@ void alpha(){
   printf("=============  Alpha source in HELIOS  ============\n");
   printf("===================================================\n");
   
+  int numEnergy = energy.size();
   printf("========= Alpha Enegry : \n");
   for( int i = 0; i < numEnergy ; i++){
     printf("%2d | %6.2f MeV\n", i, energy[i]);
@@ -54,17 +54,19 @@ void alpha(){
   //======== Set HELIOS
   printf("############################################## HELIOS configuration\n");   
   HELIOS helios;
-   helios.OverrideMagneticFieldDirection(BFieldTheta);
+  // helios.OverrideMagneticFieldDirection(BFieldTheta);
-   helios.OverrideFirstPos(-700);
+  // helios.OverrideFirstPos(-700);
   //helios.OverrideDetectorDistance(5);
-   bool sethelios = helios.SetDetectorGeometry(heliosDetGeoFile);
+  // bool sethelios = helios.SetDetectorGeometry(heliosDetGeoFile, geoID);
-   if( !sethelios){
+  // if( !sethelios){
-		helios.OverrideMagneticField(BField);
+  //   helios.OverrideMagneticField(BField);
-		printf("======== B-field : %5.2f T, Theta : %6.2f deg\n", BField, BFieldTheta);
+  //   printf("======== B-field : %5.2f T, Theta : %6.2f deg\n", BField, BFieldTheta);
-	}
+  // }
-   helios.SetCoincidentWithRecoil(isCoincidentWithRecoil);
+  // helios.SetCoincidentWithRecoil(isCoincidentWithRecoil);
-   printf("========== energy resol.: %f MeV\n", eSigma);
+  // printf("========== energy resol.: %f MeV\n", eSigma);
-   printf("=========== pos-Z resol.: %f mm \n", zSigma);
+  // printf("=========== pos-Z resol.: %f mm \n", zSigma);
+  helios.SetDetectorGeometry(heliosDetGeoFile, geoID);
+  helios.PrintGeometry();
   
   //====================== build tree
   TFile * saveFile = new TFile(saveFileName, "recreate");
@@ -87,7 +89,6 @@ void alpha(){
   tree->Branch("energyID", &energyID, "energyID/I");
   
   tree->Branch("e", &e, "e/D");
-   tree->Branch("x", &x, "x/D");
   tree->Branch("z", &z, "z/D");
   tree->Branch("t", &t, "t/D");
   tree->Branch("detID", &detID, "detID/I");
@@ -122,24 +123,26 @@ void alpha(){
     v.SetMagThetaPhi(p, theta, phi);
     
     P.SetVectM(v, ma);
+    P.SetUniqueID(2); //alpha particle has charge 2
     
     //################################### tree branches
     
     //==== Helios
-      hit = helios.CalHit(P, 2, P, 2);
+    helios.CalArrayHit(P);
+    hit = helios.CheckDetAcceptance();
     
-      e = helios.GetEnergy() + gRandom->Gaus(0, eSigma);
+    e = helios.GetEnergy() + gRandom->Gaus(0, helios.GetDetectorGeometry().array[geoID].eSigma);
-      z = helios.GetZ() ; 
-      x = helios.GetX() + gRandom->Gaus(0, zSigma);
-      t = helios.GetTime();
-      loop = helios.GetLoop();
-      detID = helios.GetDetID();
-      dphi = helios.GetdPhi();
-      rho = helios.GetRho();
-      xHit = helios.GetXPos(z);
-      yHit = helios.GetYPos(z);
-      z += gRandom->Gaus(0, zSigma);
 
+    trajectory orb = helios.GetTrajectory_b();
+
+    z = orb.z + gRandom->Gaus(0,  helios.GetDetectorGeometry().array[geoID].zSigma);
+    t = orb.t;
+    loop = orb.effLoop;
+    detID = orb.detID;
+    dphi = orb.phi;
+    rho = orb.rho;
+    xHit = orb.x;
+    yHit = orb.y;
     
     if( hit == 1) {
       count ++;
@@ -170,3 +173,9 @@ void alpha(){
   printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
   gROOT->ProcessLine(".q");
 }
+
+int main(){
+
+  alpha();
+  return 0;
+}

Cleopatra/makefile

@@ -1,6 +1,6 @@
 CC=g++
 
-ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM Transfer
+ALL = Isotope InFileCreator ExtractXSec ExtractXSecFromText PlotTGraphTObjArray Cleopatra FindThetaCM Transfer SimAlpha
 
 all: $(ALL)
 
@@ -22,11 +22,14 @@ PlotTGraphTObjArray: PlotTGraphTObjArray.C PlotTGraphTObjArray.h
 Cleopatra: Cleopatra.C
 	$(CC) Cleopatra.C -o Cleopatra `root-config --cflags --glibs`
 
-FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
+FindThetaCM: FindThetaCM.C FindThetaCM.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h
 	$(CC) FindThetaCM.C -o FindThetaCM `root-config --cflags --glibs`
 
 Transfer: Transfer.C Transfer.h ../Cleopatra/ClassTransfer.h ../Cleopatra/ClassHelios.h ../Cleopatra/ClassIsotope.h ../Cleopatra/constant.h 
 	$(CC) Transfer.C -o Transfer `root-config --cflags --glibs`
 
+SimAlpha: alpha.C ../Cleopatra/ClassHelios.h 
+	$(CC) alpha.C -o SimAlpha `root-config --cflags --glibs`
+
 clean:
 	/bin/rm -f $(ALL)