Merge branch 'NPTool.2.dev' of https://gitlab.in2p3.fr/np/nptool into NPTool.2.dev

NPLib/Physics/CMakeLists.txt

@@ -2,6 +2,8 @@ add_custom_command(OUTPUT NPParticleDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts
 
 add_custom_command(OUTPUT NPReactionDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPReaction.h NPReactionDict.cxx NPReaction.rootmap libNPPhysics.so NPReactionLinkDef.h DEPENDS NPReaction.h NPReactionLinkDef.h) 
 
+add_custom_command(OUTPUT NPInelasticBreakupDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPInelasticBreakup.h NPInelasticBreakupDict.cxx NPInelasticBreakup.rootmap libNPPhysics.so NPInelasticBreakupLinkDef.h DEPENDS NPInelasticBreakup.h NPInelasticBreakupLinkDef.h) 
+
 add_custom_command(OUTPUT NPPhaseSpaceDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPPhaseSpace.h NPPhaseSpaceDict.cxx NPPhaseSpace.rootmap libNPPhysics.so NPPhaseSpaceLinkDef.h DEPENDS NPPhaseSpace.h NPPhaseSpaceLinkDef.h) 
 
 add_custom_command(OUTPUT NPQFSDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPQFS.h NPQFSDict.cxx NPQFS.rootmap libNPPhysics.so NPQFSLinkDef.h DEPENDS NPQFS.h NPQFSLinkDef.h) 
@@ -26,7 +28,7 @@ add_custom_command(OUTPUT TFissionConditionsDict.cxx COMMAND ${CMAKE_BINARY_DIR}
 
 add_custom_command(OUTPUT NPTimeStampDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh NPTimeStamp.h NPTimeStampDict.cxx NPTimeStamp.rootmap libNPPhysics.so NPTimeStampLinkDef.h DEPENDS NPTimeStamp.h NPTimeStampLinkDef.h) 
 
-add_library(NPPhysics SHARED GEF.cxx NPFissionDecay.cxx NPDecay.cxx NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPParticle.cxx NPReaction.cxx NPTimeStamp.cxx NPPhaseSpace.cxx NPQFS.cxx NPParticleDict.cxx NPReactionDict.cxx NPTimeStampDict.cxx NPPhaseSpaceDict.cxx NPQFSDict.cxx NPEnergyLossDict.cxx )
+add_library(NPPhysics SHARED GEF.cxx NPInelasticBreakup.cxx NPInelasticBreakupDict.cxx NPFissionDecay.cxx NPDecay.cxx NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPParticle.cxx NPReaction.cxx NPTimeStamp.cxx NPPhaseSpace.cxx NPQFS.cxx NPParticleDict.cxx NPReactionDict.cxx NPTimeStampDict.cxx NPPhaseSpaceDict.cxx NPQFSDict.cxx NPEnergyLossDict.cxx )
 target_link_libraries(NPPhysics ${ROOT_LIBRARIES} Physics NPCore) 
 
 add_library(NPInitialConditions  SHARED  TInitialConditions.cxx TInitialConditionsDict.cxx )

NPLib/Physics/NPInelasticBreakup.cxx

+/*****************************************************************************
+ * Copyright (C) 2009-2016    this file is part of the NPTool Project        *
+ *                                                                           *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
+ * For the list of contributors see $NPTOOL/Licence/Contributors             *
+ *****************************************************************************/
+
+/*****************************************************************************
+ *                                                                           *
+ * Original Author : Adrien MAtta contact: matta@lpccaen.in2p3.fr            *
+ *                                                                           *
+ * Creation Date   : April 2024                                              *
+ *---------------------------------------------------------------------------*
+ * Decription:                                                               *
+ *                                                                           *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ * Comment:                                                                  *
+ *                                                                           *
+ *                                                                           *
+ *****************************************************************************/
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "NPCore.h"
+#include "NPFunction.h"
+#include "NPInelasticBreakup.h"
+#include "NPOptionManager.h"
+
+// Use CLHEP System of unit and Physical Constant
+#include "NPGlobalSystemOfUnits.h"
+#include "NPPhysicalConstants.h"
+#ifdef NP_SYSTEM_OF_UNITS_H
+using namespace NPUNITS;
+#endif
+
+// ROOT
+#include "TF1.h"
+
+ClassImp(InelasticBreakup)
+
+    //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+    InelasticBreakup::InelasticBreakup() {
+  //------------- Default Constructor -------------
+
+  //
+  fBeamEnergy = 0;
+  fThetaCM = 0;
+  fQValue = 0;
+  //  fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+  fVerboseLevel = 1; // NPOptionManager::getInstance()->GetVerboseLevel();
+
+  fCrossSectionHist = NULL;
+  fDoubleDifferentialCrossSectionHist = NULL;
+
+  fLabCrossSection = false; // flag if the provided cross-section is in the lab or not
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+InelasticBreakup::~InelasticBreakup() {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void InelasticBreakup::GenerateEvent(double& ThetaLight, double& PhiLight, double& KineticEnergyLight,
+                                     double& ThetaLabTarget, double& KineticEnergyLabTarget, double& ThetaLabHeavy,
+                                     double& KineticEnergyLabHeavy) {
+  // 2-body relativistic kinematics: direct + inverse
+  // EnergieLabTarget,4 : lab energy in MeV of the 2 ejectiles
+  // ThetaLabTarget,4   : angles in rad
+  // case of inverse kinematics
+
+  double theta = fThetaCM;
+  if (fParticle1.Mass() > fParticle2.Mass()) {
+    theta = M_PI - fThetaCM;
+    // fThetaCM = M_PI - fThetaCM;
+  }
+
+  // shoot light first
+  ThetaLight = -1;
+  KineticEnergyLight = -1;
+  while (ThetaLight < 0 || KineticEnergyLight < 0) {
+    ThetaLight = gRandom->Gaus(fMeanAngleLight, fSigmaAngleLight);
+    KineticEnergyLight = gRandom->Gaus(fMeanEnergyLight, fSigmaEnergyLight);
+  }
+
+  PhiLight = gRandom->Uniform() * 2. * M_PI;
+  double totalELight = fParticle3.Mass() + KineticEnergyLight;
+  double momentumLight = sqrt(KineticEnergyLight * KineticEnergyLight + 2 * fParticle3.Mass() * KineticEnergyLight);
+  TLorentzVector lightLV(momentumLight * sin(ThetaLight) * cos(PhiLight),
+                         momentumLight * sin(ThetaLight) * sin(PhiLight), momentumLight * cos(ThetaLight), totalELight);
+
+  // Substract light LV to Beam
+  double momentumBeam = sqrt(fBeamEnergy * fBeamEnergy + 2 * fBeam.Mass() * fBeamEnergy);
+  TLorentzVector beamLV(0, 0, momentumBeam, fBeamEnergy + fBeam.Mass());
+
+  auto newLV = beamLV - lightLV;
+
+  // Two body kine elastic scattering with heavy
+  fReaction.SetBeamEnergy(newLV.E() - newLV.Mag());
+  fReaction.ShootRandomThetaCM();
+  fReaction.KineRelativistic(ThetaLabTarget, KineticEnergyLabTarget, ThetaLabHeavy, KineticEnergyLabHeavy);
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::ReadConfigurationFile(string Path) {
+  ifstream InelasticBreakupFile;
+  string GlobalPath = getenv("NPTOOL");
+  string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
+  InelasticBreakupFile.open(Path.c_str());
+  if (!InelasticBreakupFile.is_open()) {
+    InelasticBreakupFile.open(StandardPath.c_str());
+    if (InelasticBreakupFile.is_open()) {
+      Path = StandardPath;
+    }
+    else {
+      cout << "InelasticBreakup File " << Path << " not found" << endl;
+      exit(1);
+    }
+  }
+  NPL::InputParser parser(Path);
+  ReadConfigurationFile(parser);
+}
+////////////////////////////////////////////////////////////////////////////////
+Particle InelasticBreakup::GetParticle(string name, NPL::InputParser parser) {
+  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("DefineParticle", name);
+  unsigned int size = blocks.size();
+  if (size == 0)
+    return NPL::Particle(name);
+  else if (size == 1) {
+    cout << " -- User defined nucleus " << name << " -- " << endl;
+    vector<string> token = {"SubPart", "BindingEnergy"};
+    if (blocks[0]->HasTokenList(token)) {
+      NPL::Particle N(name, blocks[0]->GetVectorString("SubPart"), blocks[0]->GetDouble("BindingEnergy", "MeV"));
+      if (blocks[0]->HasToken("ExcitationEnergy"))
+        N.SetExcitationEnergy(blocks[0]->GetDouble("ExcitationEnergy", "MeV"));
+      if (blocks[0]->HasToken("SpinParity"))
+        N.SetSpinParity(blocks[0]->GetString("SpinParity").c_str());
+      if (blocks[0]->HasToken("Spin"))
+        N.SetSpin(blocks[0]->GetDouble("Spin", ""));
+      if (blocks[0]->HasToken("Parity"))
+        N.SetParity(blocks[0]->GetString("Parity").c_str());
+      if (blocks[0]->HasToken("LifeTime"))
+        N.SetLifeTime(blocks[0]->GetDouble("LifeTime", "s"));
+
+      cout << " -- -- -- -- -- -- -- -- -- -- --" << endl;
+      return N;
+    }
+  }
+  else {
+    NPL::SendErrorAndExit("NPL::InelasticBreakup", "Too many nuclei define with the same name");
+  }
+
+  return (NPL::Particle());
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::ReadConfigurationFile(NPL::InputParser parser) {
+
+  vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("InelasticBreakup");
+  if (blocks.size() > 0 && NPOptionManager::getInstance()->GetVerboseLevel())
+    cout << endl << "\033[1;35m//// Inelastic Breakup reaction found " << endl;
+
+  vector<string> token1 = {"Beam",           "Target",          "Light",
+                           "Heavy",          "MeanEnergyLight", "SigmaEnergyLight",
+                           "MeanAngleLight", "SigmaAngleLight", "CrossSectionPath"};
+  double CSHalfOpenAngleMin = 0 * deg;
+  double CSHalfOpenAngleMax = 180 * deg;
+  for (unsigned int i = 0; i < blocks.size(); i++) {
+    if (blocks[i]->HasTokenList(token1)) {
+      int v = NPOptionManager::getInstance()->GetVerboseLevel();
+      NPOptionManager::getInstance()->SetVerboseLevel(0);
+      // Read the beam block
+      fBeam.ReadConfigurationFile(parser);
+      NPOptionManager::getInstance()->SetVerboseLevel(v);
+
+      fBeamEnergy = fBeam.GetEnergy();
+      // set the particle
+      fReaction.SetParticle1(fBeam);
+      fReaction.GetParticle1()->SetUp(blocks[i]->GetString("Heavy"));
+      fReaction.SetParticle2(GetParticle(blocks[i]->GetString("Target"), parser));
+      fReaction.SetParticle3(GetParticle(blocks[i]->GetString("Target"), parser));
+      fReaction.SetParticle4(GetParticle(blocks[i]->GetString("Heavy"), parser));
+      fReaction.SetBeamEnergy(fBeamEnergy);
+      fReaction.initializePrecomputeVariable();
+
+      fParticle3 = GetParticle(blocks[i]->GetString("Light"), parser);
+      fMeanEnergyLight = blocks[i]->GetDouble("MeanEnergyLight", "MeV");
+      fSigmaEnergyLight = blocks[i]->GetDouble("SigmaEnergyLight", "MeV");
+      fMeanAngleLight = blocks[i]->GetDouble("MeanAngleLight", "deg");
+      fSigmaAngleLight = blocks[i]->GetDouble("SigmaAngleLight", "deg");
+    }
+    else {
+      cout << "ERROR: check your input file formatting \033[0m" << endl;
+      exit(1);
+    }
+
+    if (blocks[i]->HasToken("ExcitationEnergyLight"))
+      fExcitationLight = blocks[i]->GetDouble("ExcitationEnergyLight", "MeV");
+
+    if (blocks[i]->HasToken("ExcitationEnergyHeavy"))
+      fExcitationHeavy = blocks[i]->GetDouble("ExcitationEnergyHeavy", "MeV");
+
+    if (blocks[i]->HasToken("CrossSectionPath")) {
+      vector<string> file = blocks[i]->GetVectorString("CrossSectionPath");
+      TH1F* CStemp = Read1DProfile(file[0], file[1]);
+
+      // multiply CStemp by sin(theta)
+      TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
+      CStemp->Divide(fsin, 1);
+      fReaction.SetCrossSectionHist(CStemp);
+      delete fsin;
+    }
+
+    if (blocks[i]->HasToken("LabCrossSectionPath")) {
+      fLabCrossSection = true;
+
+      vector<string> file = blocks[i]->GetVectorString("LabCrossSectionPath");
+      TH1F* CStemp = Read1DProfile(file[0], file[1]);
+
+      // multiply CStemp by sin(theta)
+      TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
+      CStemp->Divide(fsin, 1);
+      SetCrossSectionHist(CStemp);
+      delete fsin;
+    }
+
+    if (blocks[i]->HasToken("DoubleDifferentialCrossSectionPath")) {
+      vector<string> file = blocks[i]->GetVectorString("DoubleDifferentialCrossSectionPath");
+      TH2F* CStemp = Read2DProfile(file[0], file[1]);
+
+      // multiply CStemp by sin(theta)
+      // X axis is theta CM
+      // Y axis is beam energy
+      // Division affect only X axis
+      TF1* fsin = new TF1("sin", Form("1/(sin(x*%f/180.))", M_PI), 0, 180);
+      CStemp->Divide(fsin, 1);
+
+      SetDoubleDifferentialCrossSectionHist(CStemp);
+      delete fsin;
+    }
+
+    if (blocks[i]->HasToken("HalfOpenAngleMin")) {
+      CSHalfOpenAngleMin = blocks[i]->GetDouble("HalfOpenAngleMin", "deg");
+    }
+    if (blocks[i]->HasToken("HalfOpenAngleMax")) {
+      CSHalfOpenAngleMax = blocks[i]->GetDouble("HalfOpenAngleMax", "deg");
+    }
+  }
+  SetCSAngle(CSHalfOpenAngleMin / deg, CSHalfOpenAngleMax / deg);
+  cout << "\033[0m";
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void InelasticBreakup::SetCSAngle(double CSHalfOpenAngleMin, double CSHalfOpenAngleMax) {
+  if (fCrossSectionHist) {
+    for (int i = 0; i < fCrossSectionHist->GetNbinsX(); i++) {
+      if (fCrossSectionHist->GetBinCenter(i) > CSHalfOpenAngleMax ||
+          fCrossSectionHist->GetBinCenter(i) < CSHalfOpenAngleMin) {
+        fCrossSectionHist->SetBinContent(i, 0);
+      }
+    }
+  }
+}
+

NPLib/Physics/NPInelasticBreakup.h

+#ifndef NPInelasticBreakUp_h
+#define NPInelasticBreakUp_h
+/*****************************************************************************
+ * Copyright (C) 2009-2016    this file is part of the NPTool Project        *
+ *                                                                           *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
+ * For the list of contributors see $NPTOOL/Licence/Contributors             *
+ *****************************************************************************/
+
+/*****************************************************************************
+ *                                                                           *
+ * Original Author : Adrien MAtta contact: matta@lpccaen.in2p3.fr            *
+ *                                                                           *
+ * Creation Date   : April 2024                                              *
+ *---------------------------------------------------------------------------*
+ * Decription:                                                               *
+ *                                                                           *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ * Comment:                                                                  *
+ *                                                                           *
+ *                                                                           *
+ *****************************************************************************/
+// C++ header
+#include <string>
+
+// NPL
+#include "NPBeam.h"
+#include "NPInputParser.h"
+#include "NPParticle.h"
+#include "NPReaction.h"
+using namespace NPL;
+
+// ROOT header
+#include "NPReaction.h"
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "TH1F.h"
+#include "TLorentzRotation.h"
+#include "TLorentzVector.h"
+#include "TRandom.h"
+#include "TVector3.h"
+
+using namespace std;
+
+namespace NPL {
+  class InelasticBreakup {
+
+   public: // Constructors and Destructors
+    InelasticBreakup();
+    ~InelasticBreakup();
+
+   public: // Various Method
+    Particle GetParticle(string name, NPL::InputParser parser);
+    void ReadConfigurationFile(string Path);
+    void ReadConfigurationFile(NPL::InputParser);
+
+   private:
+    TRandom* RandGen;
+
+   private:
+    int fVerboseLevel;
+
+   private: // use for Monte Carlo simulation
+    bool fLabCrossSection;
+
+   private:
+    NPL::Reaction fReaction; // two body part of the kinematic
+    Beam fBeam;              // Beam
+    Beam fParticle1;         // Beam
+    Particle fParticle2;     // Target
+    Particle fParticle3;     // Light ejectile
+    Particle fParticle4;     // Heavy ejectile
+    double fMeanEnergyLight;
+    double fSigmaEnergyLight;
+    double fMeanAngleLight;
+    double fSigmaAngleLight;
+    double fExcitationHeavy;
+    double fExcitationLight;
+
+    double fQValue;                            // Q-value in MeV
+    double fBeamEnergy;                        // Beam energy in MeV
+    double fThetaCM;                           // Center-of-mass angle in radian
+    TH1F* fCrossSectionHist;                   // Differential cross section in CM frame
+    TH2F* fDoubleDifferentialCrossSectionHist; // Diff. CS CM frame vs Beam E
+    TH1F* fExcitationEnergyHist;               // Distribution of Excitation energy
+   public:
+    // Getters and Setters
+    void SetBeamEnergy(const double& eBeam) { fBeamEnergy = eBeam; }
+    void SetThetaCM(const double& angle) { fThetaCM = angle; }
+
+    void SetExcitationLight(const double& exci) { fExcitationLight = exci; }
+    void SetExcitationHeavy(const double& exci) { fExcitationHeavy = exci; }
+    void SetVerboseLevel(const int& verbose) { fVerboseLevel = verbose; }
+    void SetCrossSectionHist(TH1F* CrossSectionHist) {
+      delete fCrossSectionHist;
+      fCrossSectionHist = CrossSectionHist;
+    }
+
+    void SetDoubleDifferentialCrossSectionHist(TH2F* CrossSectionHist) {
+      fDoubleDifferentialCrossSectionHist = CrossSectionHist;
+    }
+    double GetBeamEnergy() const { return fBeamEnergy; }
+    double GetThetaCM() const { return fThetaCM; }
+    double GetQValue() const { return fQValue; }
+    Particle* GetParticleBeam() { return &fBeam; }
+    Particle* GetParticle1() { return &fParticle1; }
+    Particle* GetParticle2() { return &fParticle2; }
+    Particle* GetParticle3() { return &fParticle3; }
+    Particle* GetParticle4() { return &fParticle4; }
+    Particle* GetNucleus1() { return GetParticle1(); }
+    Particle* GetNucleus2() { return GetParticle2(); }
+    Particle* GetNucleus3() { return GetParticle3(); }
+    Particle* GetNucleus4() { return GetParticle4(); }
+
+    TH1F* GetCrossSectionHist() const { return fCrossSectionHist; }
+    int GetVerboseLevel() const { return fVerboseLevel; }
+
+   public:
+    // Modify the CS histo so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
+    void SetCSAngle(double CSHalfOpenAngleMin, double CSHalfOpenAngleMax);
+
+   public: // Kinematics
+    // Check that the reaction is alowed
+    bool IsLabCrossSection() { return fLabCrossSection; };
+
+    // Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
+    double ShootRandomThetaCM();
+
+    // Compute ThetaLab and EnergyLab for product of reaction
+    void GenerateEvent(double& ThetaLight, double& PhiLight, double& KineticEnergyLight, double& ThetaLab3,
+                       double& KineticEnergyLab3, double& ThetaLab4, double& KineticEnergyLab4);
+
+    void SetParticle3(double EnergyLab, double ThetaLab);
+
+    void SetParticle3(NPL::Particle p) { fParticle3 = p; };
+
+    ClassDef(InelasticBreakup, 0)
+  };
+} // namespace NPL
+#endif

NPLib/Physics/NPInelasticBreakupLinkDef.h

+#ifdef __CINT__
+#pragma link C++ class InelasticBreakup+;
+#endif

NPLib/Physics/NPReaction.h

@@ -173,9 +173,9 @@ namespace NPL {
    public:
     // Modify the CS histo so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
     void SetCSAngle(double CSHalfOpenAngleMin, double CSHalfOpenAngleMax);
+    void initializePrecomputeVariable();
 
    private: // intern precompute variable
-    void initializePrecomputeVariable();
     double m1;
     double m2;
     double m3;
@@ -243,7 +243,7 @@ namespace NPL {
 
     // Return Excitation Energy
     double ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab = 0);
-    
+
     // Return Lorentz Vector of Heavy Recoil after reaction
     TLorentzVector LorentzAfterReaction(double EnergyLab, double ThetaLab, double PhiLab = 0);
 
@@ -261,6 +261,11 @@ namespace NPL {
 
     void SetParticle3(double EnergyLab, double ThetaLab);
 
+    void SetParticle1(Beam p) { fParticle1 = p; };
+    void SetParticle2(Particle p) { fParticle2 = p; };
+    void SetParticle3(Particle p) { fParticle3 = p; };
+    void SetParticle4(Particle p) { fParticle4 = p; };
+
     TGraph* GetKinematicLine3(double AngleStep_CM = 1);
     TGraph* GetKinematicLine4(double AngleStep_CM = 1);
     TGraph* GetBrhoLine3(double AngleStep_CM = 1);