diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 391ffa6c0571bf27e68b88a1156b365ab328e086..1b20bfbde3c40df7b0cb83bac4de05360e7e40d0 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -7,29 +7,32 @@
/*****************************************************************************
* *
- * Original Author : F. Flavigny contact address: flavigny@ipno.in2p3.fr *
+ * Original Author : F. Flavigny contact address: flavigny@lpccaen.in2p3.fr *
* *
* Creation Date : April 2019 *
- * Last update : April 2019 *
+ * Last update : Nov 2019 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with Quasi Free Scattering Reaction in which a cluster *
* or a nucleon is removed from a projectile by interaction with a target *
* nucleon (proton target in general) *
* *
- * Labeling is: *
- * A --> i ==> B + (i -> a) => B + 1 + 2 *
+ * First step (dissociation): A -> B + c *
+ * Second step (scattering) : c + T -> 1 + 2 *
+ * Labeling is: *
+ * *
+ * A --> T ==> B + (c -> T) => B + 1 + 2 *
* *
* where: *
* + A is the beam nucleus *
- * + i is the target nucleon (proton) *
+ * + T is the target nucleon (proton) *
* *
* + B is the residual fragment (beam-like) *
- * + 1 is the scattered target nucleon (former i) *
- * + 2 is the knocked-out cluster/nucleon (noted a) in the intermediate *
+ * + 1 is the scattered target nucleon (former T) *
+ * + 2 is the knocked-out cluster/nucleon (noted c) in the intermediate *
*---------------------------------------------------------------------------*
* Comment: *
- * + *
+ * + Adapted from original event generator from V. Panin (R3B collab) *
* *
*****************************************************************************/
@@ -55,16 +58,625 @@ using namespace NPUNITS;
// ROOT
#include"TF1.h"
+#include"TRandom3.h"
ClassImp(QFS)
-/*
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+////////////////////////////////////////////////////////////////////////////////
+
QFS::QFS(){
- fVerboseLevel=0;
+
+ //------------- Default Constructor -------------
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ fBeamEnergy = 0;
+ fThetaCM = 0;
+ fPhiCM = 0;
+
+ fExcitationA = 0;
+ fExcitationB = 0;
+ fMomentumSigma = 0;
+ fshootB=false;
+ fshoot1=true;
+ fshoot2=true;
+ fisotropic = true;
+
+ fTheta2VsTheta1 = 0;
+ fPhi2VsPhi1 = 0;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+////////////////////////////////////////////////////////////////////////////////
+
QFS::~QFS(){
+
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void QFS::ReadConfigurationFile(string Path){
+ ifstream ReactionFile;
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
+ ReactionFile.open(Path.c_str());
+ if (!ReactionFile.is_open()) {
+ ReactionFile.open(StandardPath.c_str());
+ if(ReactionFile.is_open()) {
+ Path = StandardPath;
+ }
+ else {cout << "QFS File " << Path << " not found" << endl;exit(1);}
+ }
+ NPL::InputParser parser(Path);
+ ReadConfigurationFile(parser);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+void QFS::ReadConfigurationFile(NPL::InputParser parser){
+
+ cout << " In QFS ReadConfiguration " << endl;
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("QFSReaction");
+ if(blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "\033[1;35m//// QFS reaction found " << endl;
+
+ vector<string> token1 = {"Beam","Target","Scattered","KnockedOut","Heavy"};
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(token1)){
+ int v = NPOptionManager::getInstance()->GetVerboseLevel();
+ NPOptionManager::getInstance()->SetVerboseLevel(0);
+ fNucleiA.ReadConfigurationFile(parser);
+ NPOptionManager::getInstance()->SetVerboseLevel(v);
+
+ fBeamEnergy= fNucleiA.GetEnergy();
+ GetNucleus(blocks[i]->GetString("Beam"),parser);
+ fNucleiT = GetNucleus(blocks[i]->GetString("Target"),parser);
+ fNucleiB = GetNucleus(blocks[i]->GetString("Heavy"),parser);
+ fNuclei1 = GetNucleus(blocks[i]->GetString("Scattered"),parser);
+ fNuclei2 = GetNucleus(blocks[i]->GetString("KnockedOut"),parser);
+ }
+ else{
+ cout << "ERROR: check your input file formatting \033[0m" << endl;
+ exit(1);
+ }
+ if(blocks[i]->HasToken("ExcitationEnergyBeam")){
+ fExcitationA = blocks[i]->GetDouble("ExcitationEnergyBeam","MeV");
+ }
+ if(blocks[i]->HasToken("ExcitationEnergyHeavy")){
+ fExcitationB = blocks[i]->GetDouble("ExcitationEnergyHeavy","MeV");
+ }
+ if(blocks[i]->HasToken("MomentumSigma")){
+ fMomentumSigma = blocks[i]->GetDouble("MomentumSigma","MeV");
+ }
+ if(blocks[i]->HasToken("ShootHeavy")){
+ fshootB = blocks[i]->GetInt("ShootHeavy");
+ }
+ if(blocks[i]->HasToken("ShootLight")){
+ fshoot1 = blocks[i]->GetInt("ShootLight");
+ fshoot2 = blocks[i]->GetInt("ShootLight");
+ }
+ }
+
+ cout << "\033[0m" ;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+Nucleus QFS::GetNucleus(string name, NPL::InputParser parser){
+
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("DefineNucleus",name);
+ unsigned int size = blocks.size();
+ if(size==0)
+ return NPL::Nucleus(name);
+ else if(size==1){
+ cout << " -- User defined nucleus " << name << " -- " << endl;
+ vector<string> token = {"SubPart","BindingEnergy"};
+ if(blocks[0]->HasTokenList(token)){
+ NPL::Nucleus 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::QFS","Too many nuclei define with the same name");
+ }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void QFS::CalculateVariables(){
+
+ if(fBeamEnergy < 0)
+ fBeamEnergy = 0 ;
+
+ //cout<<"---- COMPUTE ------"<<endl;
+ // cout<<"--CM--"<<endl;
+
+ mA = fNucleiA.Mass(); // Beam mass in MeV
+ mT = fNucleiT.Mass(); // Target mass in MeV
+ mB = fNucleiB.Mass(); // Heavy residual mass in MeV
+ m1 = mT; // scattered target nucleon (same mass);
+ m2 = fNuclei2.Mass(); // knocked out cluster mass in MeV
+ ma = m2; // intermediate cluster mass in MeV (same);
+
+ double TA = fBeamEnergy; // Beam kinetic energy
+ double PA = sqrt(TA*(TA+2*mA)); // Beam momentum (norm)
+ double EA = sqrt(mA*mA + PA*PA); // Beam total energy
+ fEnergyImpulsionLab_A = TLorentzVector(0.,0.,PA,EA);
+
+ //Internal momentum of removed cluster/nucleon
+ static TRandom3 r;
+ //r.SetSeed(0);
+ //double mom_sigma = 0; // MeV/c
+ Pa.SetX(r.Gaus(0.,fMomentumSigma));
+ Pa.SetY(r.Gaus(0.,fMomentumSigma));
+ Pa.SetZ(r.Gaus(0.,fMomentumSigma));
+
+ //Internal momentum of heavy recoil after removal
+ PB.SetXYZ( (-Pa.X()) , (-Pa.Y()) , (-Pa.Z()) );
+
+ //cout<<"Pa_cm=\t("<<Pa.Px()<<","<<Pa.Py()<<","<<Pa.Pz()<<")"<<endl;
+ //cout<<"||Pa||^2=\t("<<Pa.Mag2()<<endl;
+ //cout<<"mA^2=\t"<<mA*mA<<endl;
+ //cout<<"mB^2=\t"<<mB*mB<<endl;
+
+
+ // Off-shell mass of the bound nucleon from E conservation
+ // in virtual dissociation of A -> B + a
+ double buffer = mA*mA + mB*mB - 2*mA*sqrt(mB*mB+Pa.Mag2()) ;
+ if(buffer<=0) { cout<<"ERROR off shell mass ma_off=\t"<<buffer<<endl; return;}
+ ma_off = sqrt(buffer);
+
+ //deduced total energies of "a" and "B" in restframe of A
+ double Ea = sqrt(ma_off*ma_off + Pa.Mag2());
+ double EB = sqrt(mB*mB + PB.Mag2());
+
+ //cout<<"ma_off^2=\t"<<buffer<<endl;
+ //cout<<"Ea=\t"<<Ea<<endl;
+ //cout<<"EB=\t"<<EB<<endl;
+
+ fEnergyImpulsionLab_a = TLorentzVector(Pa,Ea);
+ fEnergyImpulsionLab_B = TLorentzVector(PB,EB);
+ fEnergyImpulsionLab_a.Boost(0,0,fEnergyImpulsionLab_A.Beta());
+ fEnergyImpulsionLab_B.Boost(0,0,fEnergyImpulsionLab_A.Beta());
+ Ea_lab = fEnergyImpulsionLab_a.E();
+ EB_lab = fEnergyImpulsionLab_B.E();
+ Pa = fEnergyImpulsionLab_a.Vect();
+ PB = fEnergyImpulsionLab_B.Vect();
+
+ // Scattering part (2-body kinematics)
+ // virtual cluster of mass "ma_off" scattering on target T
+ // to give scattered cluster with real mass (ma=m2)
+ // and scattered target (mT=m1)
+
+ fQValue =ma_off+mT-m1-m2;
+
+ s = ma_off*ma_off + mT*mT + 2*mT*Ea_lab ;
+ fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
+ fEcm = sqrt(s) - m1 -m2;
+ if(fEcm<=0) { cout<<"ERROR Ecm negative =\t"<<fEcm<<endl; return;}
+
+ ECM_a = (s + ma_off*ma_off - mT*mT)/(2*sqrt(s));
+ ECM_T = (s + mT*mT - ma_off*ma_off)/(2*sqrt(s));
+ ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
+ ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
+
+ pCM_a = sqrt(ECM_a*ECM_a - ma_off*ma_off);
+ pCM_T = sqrt(ECM_T*ECM_T - mT*mT);
+ pCM_1 = sqrt(ECM_1*ECM_1 - m1*m1);
+ pCM_2 = sqrt(ECM_2*ECM_2 - m2*m2);
+
+ BetaCM = Pa.Mag() / (Ea_lab + mT);
+
+ //if(ECM_a<0 || ECM_T<0 || ECM_1<0 || ECM_2<0 ) {
+ /*
+ if(pCM_a<0 || pCM_T<0 || pCM_1<0 || pCM_2<0 ) {
+ cout<<"--LAB after Boost--"<<endl;
+ cout<<"Pa_lab=\t("<<Pa.Px()<<","<<Pa.Py()<<","<<Pa.Pz()<<")"<<endl;
+ cout<<"PB_lab=\t("<<PB.Px()<<","<<PB.Py()<<","<<PB.Pz()<<")"<<endl;
+ cout<<"Ea_lab=\t"<<Ea_lab<<endl;
+ cout<<"EB_lab=\t"<<EB_lab<<endl;
+ cout<<"--Back to CM--"<<endl;
+ cout<<"fQValue=\t"<<fQValue<<endl;
+ cout<<"s=\t"<<s<<endl;
+ cout<<"Ecm=\t"<<fEcm<<endl;
+ cout<<"ea*=\t"<<ECM_a<<endl;
+ cout<<"eT*=\t"<<ECM_T<<endl;
+ cout<<"e1*=\t"<<ECM_1<<endl;
+ cout<<"p1*=\t"<<pCM_1<<endl;
+ cout<<"e2*=\t"<<ECM_2<<endl;
+ cout<<"p2*=\t"<<pCM_2<<endl;
+ cout<<"beta_cm=\t"<<BetaCM<<endl;
+ }
+ */
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEnergyLab1, double& ThetaLab2, double& PhiLab2, double& KineticEnergyLab2){
+
+ CalculateVariables();
+
+ //double thetaCM_1 = fThetaCM*TMath::Pi()/180.;
+ double thetaCM_1 = fThetaCM;
+ double thetaCM_2 = M_PI - thetaCM_1;
+ //double phiCM_1 = fPhiCM*TMath::Pi()/180.;
+ double phiCM_1 = fPhiCM;
+ double phiCM_2 = 2*M_PI - phiCM_1;
+
+ TVector3 z_axis(0.,0.,1.);
+
+ fEnergyImpulsionCM_2 = TLorentzVector(
+ pCM_2*sin(thetaCM_2)*cos(phiCM_2),
+ pCM_2*cos(thetaCM_2)*sin(phiCM_2),
+ pCM_2*cos(thetaCM_2),
+ ECM_2);
+ fEnergyImpulsionCM_1 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_2;
+
+ fEnergyImpulsionLab_1 = fEnergyImpulsionCM_1;
+ fEnergyImpulsionLab_1.Boost(0,0,BetaCM);
+ fEnergyImpulsionLab_2 = fEnergyImpulsionCM_2;
+ fEnergyImpulsionLab_2.Boost(0,0,BetaCM);
+
+ // Angle in the lab frame
+ //ThetaLab1 = fEnergyImpulsionLab_1.Angle(fEnergyImpulsionLab_A.Vect());
+ ThetaLab1 = fEnergyImpulsionLab_1.Angle(z_axis);
+ if (ThetaLab1 < 0) ThetaLab1 += M_PI;
+ //ThetaLab2 = fEnergyImpulsionLab_4.Angle(fEnergyImpulsionLab_A.Vect());
+ ThetaLab2 = fEnergyImpulsionLab_2.Angle(z_axis);
+ if (fabs(ThetaLab1) < 1e-6) ThetaLab1 = 0;
+ ThetaLab2 = fabs(ThetaLab2);
+ if (fabs(ThetaLab2) < 1e-6) ThetaLab2 = 0;
+
+ PhiLab1 = M_PI + fEnergyImpulsionLab_1.Vect().Phi();
+ if (fabs(PhiLab1) < 1e-6) PhiLab1 = 0;
+ PhiLab2 = M_PI + fEnergyImpulsionLab_2.Vect().Phi();
+ if (fabs(PhiLab2) < 1e-6) PhiLab2 = 0;
+
+ // Kinetic Energy in the lab frame
+ KineticEnergyLab1 = fEnergyImpulsionLab_1.E() - m1;
+ KineticEnergyLab2 = fEnergyImpulsionLab_2.E() - m2;
+
+ // test for total energy conversion
+ //if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_1.E()+fEnergyImpulsionLab_2.E())) > 1e-6)
+ // cout << "Problem for energy conservation" << endl;
+
+/*
+ cout<<"--KINE RELATIVISTIC--"<<endl;
+ cout<<"theta_cm:"<<fThetaCM*180./TMath::Pi()<<endl;
+ cout<<"phi_cm:"<<fPhiCM*180./TMath::Pi()<<endl;
+ cout<<"P1_CM=\t("<<fEnergyImpulsionCM_1.Px()<<","<<fEnergyImpulsionCM_1.Py()<<","<<fEnergyImpulsionCM_1.Pz()<<")"<<endl;
+ cout<<"P2_CM=\t("<<fEnergyImpulsionCM_2.Px()<<","<<fEnergyImpulsionCM_2.Py()<<","<<fEnergyImpulsionCM_2.Pz()<<")"<<endl;
+ cout<<"P1_lab=\t("<<fEnergyImpulsionLab_1.Px()<<","<<fEnergyImpulsionLab_1.Py()<<","<<fEnergyImpulsionLab_1.Pz()<<")"<<endl;
+ cout<<"P2_lab=\t("<<fEnergyImpulsionLab_2.Px()<<","<<fEnergyImpulsionLab_2.Py()<<","<<fEnergyImpulsionLab_2.Pz()<<")"<<endl;
+ cout<<"Theta1:\t"<<fEnergyImpulsionLab_1.Vect().Theta()*180./TMath::Pi()<<endl;
+ cout<<"Theta2:\t"<<fEnergyImpulsionLab_2.Vect().Theta()*180./TMath::Pi()<<endl;
+ cout<<"Phi1:\t"<<M_PI+fEnergyImpulsionLab_1.Vect().Phi()*180./TMath::Pi()<<endl;
+ cout<<"Phi2:\t"<<M_PI+fEnergyImpulsionLab_2.Vect().Phi()*180./TMath::Pi()<<endl;
+*/
+
}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+double QFS::ShootRandomThetaCM(){
+
+ double theta; // CM
+/*
+ if(fDoubleDifferentialCrossSectionHist){
+ // Take a slice in energy
+ TAxis* Y = fDoubleDifferentialCrossSectionHist->GetYaxis();
+ int binY;
+
+ // Those test are there for the tail event of the energy distribution
+ // In case the energy is outside the range of the 2D histo we take the
+ // closest available CS
+ if(Y->FindBin(fBeamEnergy) > Y->GetLast())
+ binY = Y->GetLast()-1;
+
+ else if(Y->FindBin(fBeamEnergy) < Y->GetFirst())
+ binY = Y->GetFirst()+1;
+
+ else
+ binY = Y->FindBin(fBeamEnergy);
+
+ TH1D* Proj = fDoubleDifferentialCrossSectionHist->ProjectionX("proj",binY,binY);
+ SetThetaCM( theta=Proj->GetRandom()*deg );
+ }
+ else if (fLabCrossSection){
+ double thetalab=-1;
+ double energylab=-1;
+ while(energylab<0){
+ thetalab=fCrossSectionHist->GetRandom()*deg; //shoot in lab
+ energylab=EnergyLabFromThetaLab(thetalab); //get corresponding energy
+ }
+ theta = EnergyLabToThetaCM(energylab, thetalab); //transform to theta CM
+ SetThetaCM( theta );
+ }
+ else{
+ // When root perform a Spline interpolation to shoot random number out of
+ // the distribution, it can over shoot and output a number larger that 180
+ // this lead to an additional signal at 0-4 deg Lab, especially when using a
+ // flat distribution.
+ // This fix it.
+ theta=181;
+ if(theta/deg>180)
+ theta=fCrossSectionHist->GetRandom();
+ //cout << " Shooting Random ThetaCM " << theta << endl;
+ SetThetaCM( theta*deg );
+ }
*/
+ return theta;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+TGraph* QFS::GetTheta2VsTheta1(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
+ double theta1,phi1,E1,theta2,phi2,E2;
+
+ for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*TMath::Pi()/180.);
+ KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
+
+ vx.push_back(theta1*180./M_PI);
+ vy.push_back(theta2*180./M_PI);
+ }
+ fTheta2VsTheta1 = new TGraph(vx.size(),&vx[0],&vy[0]);
+
+ return(fTheta2VsTheta1);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+TGraph* QFS::GetPhi2VsPhi1(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
+ double theta1,phi1,E1,theta2,phi2,E2;
+
+ for (double theta=0 ; theta <= 180 ; theta+=AngleStep_CM){
+ for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
+ SetThetaCM(theta*TMath::Pi()/180.);
+ SetPhiCM(angle*TMath::Pi()/180.);
+ KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
+ vx.push_back(phi1*180./M_PI);
+ vy.push_back(phi2*180./M_PI);
+ }
+ }
+ fPhi2VsPhi1 = new TGraph(vx.size(),&vx[0],&vy[0]);
+
+ return(fPhi2VsPhi1);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Check whenever the reaction is allowed at a given energy
+bool QFS::IsAllowed(){//double Energy){
+ //double AvailableEnergy = Energy + fNuclei1.Mass() + fNuclei2.Mass();
+ //double RequiredEnergy = fNuclei3.Mass() + fNuclei4.Mass();
+
+ //if(AvailableEnergy>RequiredEnergy)
+ return true;
+ //else
+ // return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////
+///////////////// Old R3B method not using TLorentz Vector ////////////////////////////////
+/////////////////// (used as a reference)///////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////
+/*
+void QFS::CalculateVariablesOld(){
+
+ initializePrecomputeVariable();
+
+ vector<double> theta1;
+ vector<double> theta2;
+ vector<double> phi1;
+ vector<double> phi2;
+
+ //for(int i=0; i<=180; i++){
+ int i = 29;
+ KineR3B(s, ma_off, mT, ma, (double)i);
+ if(!good) { cout<<"ERROR CM calculations!!!"<<endl; return;}
+
+ TVector3 P1_cm(0.,0.,1.), P2_cm(0.,0.,1.);
+ P2_cm.SetMag(p_clust);
+ P2_cm.SetTheta(theta_clust);
+ TRandom3 ra;
+ ra.SetSeed(0);
+ //P2_cm.SetPhi(ra.Uniform(-1*TMath::Pi(),+1*TMath::Pi()));
+ P2_cm.SetPhi(0.);
+ P1_cm.SetX(-P2_cm.X());
+ P1_cm.SetY(-P2_cm.Y());
+ P1_cm.SetZ(-P2_cm.Z());
+ cout<<"FFFFFFFFFFFFFFFFF"<<endl;
+ cout<<"P1_CM=\t("<<P1_cm.X()<<","<<P1_cm.Y()<<","<<P1_cm.Z()<<")"<<endl;
+ cout<<"P2_CM=\t("<<P2_cm.X()<<","<<P2_cm.Y()<<","<<P2_cm.Z()<<")"<<endl;
+
+
+ // Calculate relative to direction of quasi-particle (cluster)
+
+ double beta_cm = -Pa.Mag() / (Ea_lab + mT);
+ double gamma_cm = 1/sqrt(1-beta_cm*beta_cm);
+
+ pair<double,double> lor_a1 = Lorentz(gamma_cm,beta_cm,e_scat,P1_cm.Z());
+ pair<double,double> lor_a2 = Lorentz(gamma_cm,beta_cm,e_clust,P2_cm.Z());
+
+ P1_cm.SetZ(lor_a1.second);
+ P2_cm.SetZ(lor_a2.second);
+
+ //Rotating back to beam direction
+ //TVector3 P1_L = Rotations(P1_cm, Pa);
+ //TVector3 P2_L = Rotations(P2_cm, Pa);
+
+ TVector3 P1_L = P1_cm;
+ TVector3 P2_L = P2_cm;
+ TVector3 direction = Pa.Unit();
+ //P1_L.RotateUz(direction);
+ //P1_L.RotateUz(direction);
+
+ cout<<"----Calculate variables output------"<<endl;
+ cout<<"--CM--"<<endl;
+ cout<<"theta1*=\t"<<theta_scat*180/TMath::Pi()<<endl;
+ cout<<"theta2*=\t"<<theta_clust*180/TMath::Pi()<<endl;
+ cout<<"e1*=\t"<<e_scat<<endl;
+ cout<<"p1*=\t"<<p_scat<<endl;
+ cout<<"e2*=\t"<<e_clust<<endl;
+ cout<<"p2*=\t"<<p_clust<<endl;
+ cout<<"T=\t"<<T<<endl;
+ cout<<"beta_cm=\t"<<beta_cm<<endl;
+ cout<<"gamma_cm=\t"<<gamma_cm<<endl;
+
+ cout<<"--LAB (cluster dir)--"<<endl;
+ cout<<"P1_lab=\t("<<P1_cm.X()<<","<<P1_cm.Y()<<","<<P1_cm.Z()<<")"<<endl;
+ cout<<"P2_lab=\t("<<P2_cm.X()<<","<<P2_cm.Y()<<","<<P2_cm.Z()<<")"<<endl;
+ cout<<"Theta1:\t"<<P1_cm.Theta()*180./TMath::Pi()<<endl;
+ cout<<"Theta2:\t"<<P2_cm.Theta()*180./TMath::Pi()<<endl;
+
+ cout<<"--LAB--"<<endl;
+ cout<<"Theta1L:\t"<<P1_L.Theta()*180./TMath::Pi()<<endl;
+ cout<<"Theta2L:\t"<<P2_L.Theta()*180./TMath::Pi()<<endl;
+ cout<<"Phi1L:\t"<<P1_L.Phi()*180./TMath::Pi()<<endl;
+ cout<<"Phi2L:\t"<<P2_L.Phi()*180./TMath::Pi()<<endl;
+
+ //cout<<P1_cm.Theta()*180./TMath::Pi()<<"\t"<<P2_cm.Theta()*180./TMath::Pi()<<endl;
+ //cout<<P1_L.Phi()*180./TMath::Pi()<<"\t"<<P2_L.Phi()*180./TMath::Pi()<<endl;
+
+ theta1.push_back(P1_L.Theta()*180./TMath::Pi());
+ theta2.push_back(P2_L.Theta()*180./TMath::Pi());
+
+ double temp_phi1 = P1_L.Phi();
+ double temp_phi2 = P2_L.Phi();
+ phi1.push_back(180. + P1_L.Phi()*180./TMath::Pi());
+ phi2.push_back(180. + P2_L.Phi()*180./TMath::Pi());
+
+ //}
+ TGraph* fTheta2VsTheta1 = new TGraph(theta1.size(),&theta1[0],&theta2[0]);
+ TGraph* fPhi2VsPhi1 = new TGraph(phi1.size(),&phi1[0],&phi2[0]);
+ fTheta2VsTheta1->SetName("Theta2VsTheta1");
+ fPhi2VsPhi1->SetName("Phi2VsPhi1");
+ TFile* f = new TFile("graphs.root","RECREATE");
+ fTheta2VsTheta1->Write();
+ fPhi2VsPhi1->Write();
+ f->Close();
+
+
+}
+
+// Calculate elastic scattering kinematics in CM-system (1-target proton, 2-cluster)
+void QFS::KineR3B(double s,double m2off,double m1,double m2,double thetacm)
+{
+ if(thetacm>180 || thetacm<0){
+ cout << "\nERROR! ThetaCM (in deg) should be between 0 and 180"<<endl;
+ return;
+ }
+
+ e_clust = 0;
+ p_clust = 0;
+ theta_clust = 0;
+ e_scat = 0;
+ p_scat = 0;
+ theta_scat = 0;
+ T = 0;
+ good = false;
+
+
+ double X = s;
+ double Y = m2off*m2off;
+ double Z = m1*m1;
+ double sqrt_s = sqrt(s);
+
+ // Kinematics before the scattering process
+ // (with one off-shell mass)
+ double p2_off = sqrt(function(X,Y,Z))/2/sqrt_s;
+ double p1_off = p2_off;
+ // CM energies
+ double e1_off = (s+Z-Y)/2/sqrt_s;
+ double e2_off = (s+Y-Z)/2/sqrt_s;
+
+ // Now take the real masses (after scattering)
+ Y = m2*m2; Z = m1*m1;
+ //And check whether the kinematical function is ok
+ //for this specific kinematical case
+ double ERROR_CI = function(X,Y,Z);
+ if(ERROR_CI <= 0.){
+ cout << "\nERROR!!! Kinematical function is negative!";
+ return;
+ }
+
+ // Kinematics after the scattering process
+ // (with all real masses)
+ double p2 = sqrt(function(X,Y,Z))/2/sqrt_s;
+ double p1 = p2;
+ double e1 = (s+Z-Y)/2/sqrt_s;
+ double e2 = (s+Y-Z)/2/sqrt_s;
+
+ // Let's consider momentum transfer <t> from the
+ // target particle 1 to the cluster 2
+ double t = 2*(m1*m1 - e1_off*e1 + p1_off*p1*cos(thetacm*TMath::Pi()/180.));
+
+ //CM scattering angles
+ double theta1 = thetacm*TMath::Pi()/180.;
+ double theta2 = TMath::Pi() - theta1;
+
+ e_clust = e2;
+ p_clust = p2;
+ theta_clust = theta2;
+
+ e_scat = e1;
+ p_scat = p1;
+ theta_scat = theta1;
+
+ T = t;
+ good = true;
+
+
+}
+
+
+
+double QFS::function(double x,double y,double z)
+{
+ double lambda = x*x + y*y + z*z - 2*x*y - 2*x*z - 2*y*z;
+ return lambda;
+}
+
+//---- Two consecutive rotations
+//first around Z on <phi>, then around new X' on <theta> (1=Pcm, 2=Pa in lab)
+TVector3 QFS::Rotations(TVector3 v1,TVector3 v2)
+{
+ double CT = v2.Z()/v2.Mag(); // cos(theta) of v2 wrt. Z-axis
+ double ST = sqrt(1-CT*CT); // sin(theta)
+ double CF = v2.X()/v2.Mag()/ST;
+ double SF = v2.Y()/v2.Mag()/ST;
+
+ TVector3 v3;
+ double _v3x = v1.X()*CT*CF - v1.Y()*SF + v1.Z()*ST*CF;
+ double _v3y = v1.X()*CT*SF + v1.Y()*CF + v1.Z()*ST*SF;
+ double _v3z = -v1.X()*ST + v1.Z()*CT;
+ v3.SetXYZ(_v3x,_v3y,_v3z);
+ return v3;
+}
+
+
+
+pair<double, double> QFS::Lorentz(double gamma,double beta,double e,double p)
+{
+ double eL = gamma*e - gamma*beta*p;
+ double pL = gamma*p - gamma*beta*e;
+ return make_pair(eL, pL);
+}
+*/
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index b0107b6c85c3def61911987017ee613fa918ed36..0d7f63c753157247acd905f10497bcd11f7eb4a7 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -1,6 +1,5 @@
#ifndef NPQFS_h
#define NPQFS_h
-/*****************************************************************************/
/*****************************************************************************
* Copyright (C) 2009-2019 this file is part of the NPTool Project *
* *
@@ -10,29 +9,32 @@
/*****************************************************************************
* *
- * Original Author : F. Flavigny contact address: flavigny@ipno.in2p3.fr *
+ * Original Author : F. Flavigny contact address: flavigny@lpccaen.in2p3.fr *
* *
* Creation Date : April 2019 *
- * Last update : April 2019 *
+ * Last update : Nov 2019 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with Quasi Free Scattering Reaction in which a cluster *
* or a nucleon is removed from a projectile by interaction with a target *
* nucleon (proton target in general) *
* *
- * Labeling is: *
- * A --> i ==> B + (i -> a) => B + 1 + 2 *
+ * First step (dissociation): A -> B + c *
+ * Second step (scattering) : c + T -> 1 + 2 *
+ * Labeling is: *
+ * *
+ * A --> T ==> B + (c -> T) => B + 1 + 2 *
* *
* where: *
* + A is the beam nucleus *
- * + i is the target nucleon (proton) *
+ * + T is the target nucleon (proton) *
* *
* + B is the residual fragment (beam-like) *
- * + 1 is the scattered target nucleon (former i) *
- * + 2 is the knocked-out cluster/nucleon (noted a) in the intermediate *
+ * + 1 is the scattered target nucleon (former T) *
+ * + 2 is the knocked-out cluster/nucleon (noted c) in the intermediate *
*---------------------------------------------------------------------------*
* Comment: *
- * + *
+ * + Adapted from original event generator from V. Panin (R3B collab) *
* *
*****************************************************************************/
@@ -41,6 +43,7 @@
// NPL
#include "NPNucleus.h"
+#include "NPBeam.h"
#include "NPInputParser.h"
using namespace NPL;
@@ -60,12 +63,138 @@ namespace NPL{
class QFS{
public: // Constructors and Destructors
- QFS(){};
- ~QFS(){};
+ QFS();
+ ~QFS();
private:
int fVerboseLevel;
+ Beam fNucleiA; // Beam (A)
+ Nucleus fNucleiT; // Target (T)
+ Nucleus fNucleiB; // Beam-like ejectile (B)
+ Nucleus fNuclei1; // Target-like ejectile (1)
+ Nucleus fNuclei2; // Knocked-out nucleon/cluster (2)
+ double fQValue; // Q-value in MeV
+ double fEcm; // Ecm in MeV
+ double fThetaCM; // Center-of-mass theta angle in radian
+ double fPhiCM; // Center-of-mass Phi angle in radian
+ double fBeamEnergy; // Beam energy in MeV
+ double fExcitationA; // Excitation energy in MeV of the beam, useful for isomers
+ double fExcitationB; // Excitation energy in MeV of beam-like ejectile
+ double fMomentumSigma; // Width of the momentum distribution of the removed cluster (sigma)
+ double fisotropic;
+
+ TGraph* fTheta2VsTheta1;
+ TGraph* fPhi2VsPhi1;
+
+ // used for MC simulations
+ bool fshootB; // shoot beam-like ejectile
+ bool fshoot1; // shoot light ejectile &
+ bool fshoot2; // shoot light ejectile 2
+
+ public:
+ Nucleus GetNucleus(string name, NPL::InputParser parser);
+ void ReadConfigurationFile(string Path);
+ void ReadConfigurationFile(NPL::InputParser);
+ void CalculateVariables();
+ void KineRelativistic(double &ThetaLab1, double &PhiLab1, double &KineticEnergyLab1, double &ThetaLab2, double &PhiLab2, double &KineticEnergyLab2);
+ double ShootRandomThetaCM();
+ bool IsAllowed();
+
+ private: // intern precompute variable
+ double mA;
+ double mB;
+ double ma_off;
+ double ma;
+ double mT;
+ double m1;
+ double m2;
+
+ TVector3 Pa, PB;
+ double Ea_lab,EB_lab;
+
+ // Lorentz Vector
+ TLorentzVector fEnergyImpulsionLab_A;
+ TLorentzVector fEnergyImpulsionLab_B;
+ TLorentzVector fEnergyImpulsionLab_a;
+ TLorentzVector fEnergyImpulsionLab_T;
+ TLorentzVector fEnergyImpulsionLab_1;
+ TLorentzVector fEnergyImpulsionLab_2;
+ TLorentzVector fTotalEnergyImpulsionLab;
+
+ TLorentzVector fEnergyImpulsionCM_a;
+ TLorentzVector fEnergyImpulsionCM_T;
+ TLorentzVector fEnergyImpulsionCM_1;
+ TLorentzVector fEnergyImpulsionCM_2;
+ TLorentzVector fTotalEnergyImpulsionCM;
+
+ // Impulsion Vector3
+ TVector3 fImpulsionLab_a;
+ TVector3 fImpulsionLab_T;
+ TVector3 fImpulsionLab_1;
+ TVector3 fImpulsionLab_2;
+
+ TVector3 fImpulsionCM_a;
+ TVector3 fImpulsionCM_T;
+ TVector3 fImpulsionCM_1;
+ TVector3 fImpulsionCM_2;
+
+ // CM Energy composante & CM impulsion norme
+ Double_t ECM_a;
+ Double_t ECM_T;
+ Double_t ECM_1;
+ Double_t ECM_2;
+ Double_t pCM_a;
+ Double_t pCM_T;
+ Double_t pCM_1;
+ Double_t pCM_2;
+
+ // Mandelstam variable
+ Double_t s;
+
+ // Center of Mass Kinematic
+ Double_t BetaCM;
+
+ public:
+ //SETTERS
+ void SetBeamEnergy(const double& eBeam) {fBeamEnergy = eBeam;}
+ void SetThetaCM(const double& angle) {fThetaCM = angle;}
+ void SetPhiCM(const double& angle) {fPhiCM = angle;}
+
+ //GETTERS
+ Nucleus* GetNucleusA() {return &fNucleiA;}
+ Nucleus* GetNucleusT() {return &fNucleiT;}
+ Nucleus* GetNucleusB() {return &fNucleiB;}
+ Nucleus* GetNucleus1() {return &fNuclei1;}
+ Nucleus* GetNucleus2() {return &fNuclei2;}
+ bool GetShoot1() const {return fshoot1;}
+ bool GetShoot2() const {return fshoot2;}
+ bool GetShootB() const {return fshootB;}
+ TLorentzVector* GetEnergyImpulsionLab_B() {return &fEnergyImpulsionLab_B;}
+ TGraph* GetTheta2VsTheta1(double AngleStep_CM=1);
+ TGraph* GetPhi2VsPhi1(double AngleStep_CM=1);
+
+
+
+ /*
+ //TO REMOVE AT SOME POINT WHEN CLASS IS ROBUSTLY TESTED
+ private:
+ // R3B Methods and Variables used as a starting point for this class (useful for checks)
+ void CalculateVariablesOld();
+ pair<double,double> Lorentz(double, double, double, double);
+ double function(double, double, double);
+ void KineR3B(double, double, double, double, double);
+ TVector3 Rotations(TVector3,TVector3);
+ double e_clust;
+ double p_clust;
+ double theta_clust;
+ double e_scat;
+ double p_scat;
+ double theta_scat;
+ bool good;
+ double T;
+ */
+
ClassDef(QFS,0)
};
}
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 5a94c3020a0e63b004ac278a9c5469f4e5780139..07f42fcf5a43f69cea71093467eedfa6b5a1896a 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -58,19 +58,37 @@ void NPS::BeamReaction::AttachReactionConditions() {
////////////////////////////////////////////////////////////////////////////////
void NPS::BeamReaction::ReadConfiguration() {
NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
- m_Reaction.ReadConfigurationFile(input);
- m_BeamName = NPL::ChangeNameToG4Standard(m_Reaction.GetNucleus1()->GetName());
-
- if (m_Reaction.GetNucleus3()->GetName() != "") {
- m_active = true;
- m_ReactionConditions = new TReactionConditions();
- AttachReactionConditions();
- if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
- RootOutput::getInstance()->GetTree()->Branch(
- "ReactionConditions", "TReactionConditions", &m_ReactionConditions);
-
- } else {
- m_active = false;
+
+ //vector<NPL::InputBlock*> blocks;
+ //blocks = input.GetAllBlocksWithToken("TwoBodyReaction");
+ //if(blocks.size()>0) m_ReactionType ="TwoBodyReaction";
+ //
+ //blocks = input.GetAllBlocksWithToken("QFSReaction");
+ //if(blocks.size()>0) m_ReactionType ="QFSReaction";
+
+ if(input.GetAllBlocksWithToken("TwoBodyReaction").size()>0) m_ReactionType ="TwoBodyReaction";
+ if(input.GetAllBlocksWithToken("QFSReaction").size()>0) m_ReactionType ="QFSReaction";
+
+
+ if (m_ReactionType=="TwoBodyReaction" ) {
+ m_Reaction.ReadConfigurationFile(input);
+ m_BeamName = NPL::ChangeNameToG4Standard(m_Reaction.GetNucleus1()->GetName());
+ if(m_Reaction.GetNucleus3()->GetName() != ""){
+ m_active = true;
+ m_ReactionConditions = new TReactionConditions();
+ AttachReactionConditions();
+ if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+ RootOutput::getInstance()->GetTree()->Branch(
+ "ReactionConditions", "TReactionConditions", &m_ReactionConditions);
+ }
+ } else if (m_ReactionType=="QFSReaction") {
+ m_QFS.ReadConfigurationFile(input);
+ m_BeamName = NPL::ChangeNameToG4Standard(m_QFS.GetNucleusA()->GetName());
+ m_active = true;
+ m_ReactionConditions = new TReactionConditions();
+ }
+ else {
+ m_active = false;
}
}
@@ -90,6 +108,8 @@ NPS::BeamReaction::IsApplicable(const G4ParticleDefinition& particleType) {
////////////////////////////////////////////////////////////////////////////////
G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
+
+ //cout<< "MODEL TRIG"<<endl;
static bool shoot = false;
static double rand = 0;
const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
@@ -103,9 +123,13 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
double out = solid->DistanceToOut(P, -V);
double ratio = in / (out + in);
+ //cout<< "in:"<<in<<std::scientific<<endl;
+ //cout<< "ou:"<<out<<std::scientific<<endl;
+ //cout<< "ratio:"<<ratio<<std::scientific<<endl;
// m_StepSize = m_StepSize/100000.;
if (out == 0) { // first step of current event
+ //cout<< "FIRST"<<endl;
rand = G4RandFlat::shoot();
m_PreviousLength = m_StepSize;
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
@@ -113,21 +137,35 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
m_ReactionConditions->Clear();
shoot = true;
}
-
- else if (in == 0) { // last step
+ else if ((in-m_StepSize) <= 1E-9) { // last step
+ //cout<< "LAST"<<endl;
return true;
}
+ //cout.precision(17);
+ //cout<< "rand:"<<rand<<std::scientific<<endl;
+
// If the condition is met, the event is generated
if (ratio < rand) {
+
// Reset the static for next event
// shoot = false;
- if (shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy())) {
- shoot = false;
- return true;
-
- } else {
- return false;
+ if(m_ReactionType=="QFSReaction"){
+ if ( shoot && m_QFS.IsAllowed() ) {
+ shoot = false;
+ return true;
+ } else {
+ return false;
+ }
+ }else if(m_ReactionType=="TwoBodyReaction"){
+ if ( shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
+ shoot = false;
+ return true;
+ } else {
+ return false;
+ }
+ }else{
+ return false;
}
}
@@ -135,6 +173,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
// so it can be used in the next one
if (!PrimaryTrack->GetStep()->IsLastStepInVolume()) {
m_PreviousLength = PrimaryTrack->GetStep()->GetStepLength();
+ //cout<< "PreviousLength="<<m_PreviousLength<<endl;
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
}
@@ -143,207 +182,337 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
////////////////////////////////////////////////////////////////////////////////
void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
- G4FastStep& fastStep) {
-
- // Get the track info
- const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
- G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
- G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
-
- G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
- G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
-
- double energy = PrimaryTrack->GetKineticEnergy();
- double time = PrimaryTrack->GetGlobalTime();
-
- // Randomize within the step
- // Assume energy loss is linear within the step
- // Assume no scattering
- double rand = G4RandFlat::shoot();
- double length = rand * (m_PreviousLength);
- energy -= (1 - rand) * (m_PreviousEnergy - energy);
- G4ThreeVector ldir = pdirection;
- ldir *= length;
- localPosition = localPosition - ldir;
- // Set the end of the step conditions
- fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
- fastStep.SetPrimaryTrackFinalPosition(worldPosition);
- fastStep.SetTotalEnergyDeposited(0);
- fastStep.SetPrimaryTrackFinalTime(time);
- fastStep.KillPrimaryTrack();
- fastStep.SetPrimaryTrackPathLength(0.0);
-
- /////////////////////////////////////////////////
- //////Define the kind of particle to shoot////////
- //////////////////////////////////////////////////
-
- // Nucleus 3
- int LightZ = m_Reaction.GetNucleus3()->GetZ();
- int LightA = m_Reaction.GetNucleus3()->GetA();
- static G4IonTable* IonTable
- = G4ParticleTable::GetParticleTable()->GetIonTable();
-
- G4ParticleDefinition* LightName;
-
- if (LightZ == 0 && LightA == 1) // neutron is special case
- {
- LightName = G4Neutron::Definition();
- } else {
- if (m_Reaction.GetUseExInGeant4())
- LightName
- = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3() * MeV);
- else
- LightName = IonTable->GetIon(LightZ, LightA);
- }
+ G4FastStep& fastStep) {
+
+ // Get the track info
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
+ G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
+
+ G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
+ G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
+
+ double energy = PrimaryTrack->GetKineticEnergy();
+ double time = PrimaryTrack->GetGlobalTime();
+
+ // Randomize within the step
+ // Assume energy loss is linear within the step
+ // Assume no scattering
+ double rand = G4RandFlat::shoot();
+ double length = rand * (m_PreviousLength);
+ energy -= (1 - rand) * (m_PreviousEnergy - energy);
+ G4ThreeVector ldir = pdirection;
+ ldir *= length;
+ localPosition = localPosition - ldir;
+ // Set the end of the step conditions
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetTotalEnergyDeposited(0);
+ fastStep.SetPrimaryTrackFinalTime(time);
+ fastStep.KillPrimaryTrack();
+ fastStep.SetPrimaryTrackPathLength(0.0);
+
+
+ if (m_ReactionType=="TwoBodyReaction" ) {
+
+ ///////////////////////////////
+ // Two-Body Reaction Case /////
+ ///////////////////////////////
+
+ //////Define the kind of particle to shoot////////
+ // Nucleus 3
+ int LightZ = m_Reaction.GetNucleus3()->GetZ();
+ int LightA = m_Reaction.GetNucleus3()->GetA();
+ static G4IonTable* IonTable
+ = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+ G4ParticleDefinition* LightName;
+
+ if (LightZ == 0 && LightA == 1) // neutron is special case
+ {
+ LightName = G4Neutron::Definition();
+ } else {
+ if (m_Reaction.GetUseExInGeant4())
+ LightName
+ = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3() * MeV);
+ else
+ LightName = IonTable->GetIon(LightZ, LightA);
+ }
+
+ // Nucleus 4
+ G4int HeavyZ = m_Reaction.GetNucleus4()->GetZ();
+ G4int HeavyA = m_Reaction.GetNucleus4()->GetA();
+
+ // Generate the excitation energy if a distribution is given
+ m_Reaction.ShootRandomExcitationEnergy();
+
+ // Use to clean up the IonTable in case of the Ex changing at every event
+ G4ParticleDefinition* HeavyName;
+
+ if (m_Reaction.GetUseExInGeant4())
+ HeavyName
+ = IonTable->GetIon(HeavyZ, HeavyA, m_Reaction.GetExcitation4() * MeV);
+ else
+ HeavyName = IonTable->GetIon(HeavyZ, HeavyA);
+
+ // Set the Energy of the reaction
+ m_Reaction.SetBeamEnergy(energy);
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+
+ ///////////////////////////
+ ///// Beam Parameters /////
+ ///////////////////////////
+ m_ReactionConditions->SetBeamParticleName(
+ PrimaryTrack->GetParticleDefinition()->GetParticleName());
+
+ m_ReactionConditions->SetBeamReactionEnergy(energy);
+ m_ReactionConditions->SetVertexPositionX(localPosition.x());
+ m_ReactionConditions->SetVertexPositionY(localPosition.y());
+ m_ReactionConditions->SetVertexPositionZ(localPosition.z());
+
+ G4ThreeVector U(1, 0, 0);
+ G4ThreeVector V(0, 1, 0);
+ G4ThreeVector ZZ(0, 0, 1);
+ m_ReactionConditions->SetBeamEmittanceTheta(
+ PrimaryTrack->GetMomentumDirection().theta() / deg);
+ m_ReactionConditions->SetBeamEmittancePhi(
+ PrimaryTrack->GetMomentumDirection().phi() / deg);
+ m_ReactionConditions->SetBeamEmittanceThetaX(
+ PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+ m_ReactionConditions->SetBeamEmittancePhiY(
+ PrimaryTrack->GetMomentumDirection().angle(V) / deg);
+
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+
+
+ // G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ // cos(Beam_theta) * sin(Beam_phi),
+ // -sin(Beam_theta));
+ // G4ThreeVector col2(-sin(Beam_phi),
+ // cos(Beam_phi),
+ // 0);
+ // G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ // sin(Beam_theta) * sin(Beam_phi),
+ // cos(Beam_theta));
+ // G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+
+ // Angles
+ // Shoot and Set a Random ThetaCM
+ m_Reaction.ShootRandomThetaCM();
+ double phi = RandFlat::shoot() * 2. * pi;
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ double Theta3, Energy3, Theta4, Energy4;
+
+ // Compute Kinematic using previously defined ThetaCM
+ m_Reaction.KineRelativistic(Theta3, Energy3, Theta4, Energy4);
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kine3_beam(sin(Theta3) * cos(phi),
+ sin(Theta3) * sin(phi), cos(Theta3));
+ // Momentum in World frame //to go from the incident beam frame to the "world"
+ // frame
+ G4ThreeVector momentum_kine3_world = momentum_kine3_beam;
+ momentum_kine3_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kine3_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kine4_beam(sin(Theta4) * cos(phi + pi),
+ sin(Theta4) * sin(phi + pi), cos(Theta4));
+ // Momentum in World frame
+ G4ThreeVector momentum_kine4_world = momentum_kine4_beam;
+ momentum_kine4_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kine4_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Emitt secondary
+ if (m_Reaction.GetShoot3()) {
+ G4DynamicParticle particle3(LightName, momentum_kine3_world, Energy3);
+ fastStep.CreateSecondaryTrack(particle3, localPosition, time);
+ }
+
+ if (m_Reaction.GetShoot4()) {
+ G4DynamicParticle particle4(HeavyName, momentum_kine4_world, Energy4);
+ fastStep.CreateSecondaryTrack(particle4, localPosition, time);
+ }
+
+ // Reinit for next event
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+
+ ///////////////////////////////////////
+ ///// Emitted particle Parameters /////
+ ///////////////////////////////////////
+ // Names 3 and 4//
+ m_ReactionConditions->SetParticleName(LightName->GetParticleName());
+ m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
+ // Angle 3 and 4 //
+ m_ReactionConditions->SetTheta(Theta3 / deg);
+ m_ReactionConditions->SetTheta(Theta4 / deg);
+
+ m_ReactionConditions->SetPhi(phi / deg);
+ if ((phi + pi) / deg > 360)
+ m_ReactionConditions->SetPhi((phi - pi) / deg);
+ else
+ m_ReactionConditions->SetPhi((phi + pi) / deg);
+
+ // Energy 3 and 4 //
+ m_ReactionConditions->SetKineticEnergy(Energy3);
+ m_ReactionConditions->SetKineticEnergy(Energy4);
+ // ThetaCM and Ex//
+ m_ReactionConditions->SetThetaCM(m_Reaction.GetThetaCM() / deg);
+ m_ReactionConditions->SetExcitationEnergy3(m_Reaction.GetExcitation3());
+ m_ReactionConditions->SetExcitationEnergy4(m_Reaction.GetExcitation4());
+ // Momuntum X 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kine3_world.x());
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kine4_world.x());
+ // Momuntum Y 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kine3_world.y());
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kine4_world.y());
+ // Momuntum Z 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kine3_world.z());
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kine4_world.z());
+
+ }// end if TwoBodyReaction
+
+ else if (m_ReactionType=="QFSReaction" ) {
+
+ //////Define the kind of particle to shoot////////
+ // A --> T ==> B + (c -> T) => B + 1 + 2
+
+ int Light1_Z = m_QFS.GetNucleus1()->GetZ();
+ int Light1_A = m_QFS.GetNucleus1()->GetA();
+ int Light2_Z = m_QFS.GetNucleus2()->GetZ();
+ int Light2_A = m_QFS.GetNucleus2()->GetA();
+
+ static G4IonTable* IonTable
+ = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+ G4ParticleDefinition* Light1Name;
+ G4ParticleDefinition* Light2Name;
+
+ if (Light1_Z == 0 && Light1_A == 1) // neutron is special case
+ {
+ Light1Name = G4Neutron::Definition();
+ } else {
+ Light1Name = IonTable->GetIon(Light1_Z, Light1_A);
+ }
+
+ if (Light2_Z == 0 && Light2_A == 1) // neutron is special case
+ {
+ Light2Name = G4Neutron::Definition();
+ } else {
+ Light2Name = IonTable->GetIon(Light2_Z, Light2_A);
+ }
+
+ // Nucleus B
+ G4int Heavy_Z = m_QFS.GetNucleusB()->GetZ();
+ G4int Heavy_A = m_QFS.GetNucleusB()->GetA();
+
+ G4ParticleDefinition* HeavyName;
+ HeavyName = IonTable->GetIon(Heavy_Z, Heavy_A);
+
+ // Set the Energy of the reaction
+ m_QFS.SetBeamEnergy(energy);
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+
+ // Angles
+ // Shoot and Set a Random ThetaCM
+ //m_QFS.ShootRandomThetaCM();
+ //m_QFS.ShootRandomPhiCM();
+ double theta = RandFlat::shoot() * pi;
+ double phi = RandFlat::shoot() * 2. * pi;
+
+ m_QFS.SetThetaCM(theta);
+ m_QFS.SetPhiCM(phi);
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ double Theta1, Phi1, Energy1, Theta2, Phi2, Energy2;
+
+ // Compute Kinematic using previously defined ThetaCM
+ m_QFS.KineRelativistic(Theta1, Phi1, Energy1, Theta2, Phi2, Energy2);
+
+ G4ThreeVector U(1, 0, 0);
+ G4ThreeVector V(0, 1, 0);
+ G4ThreeVector ZZ(0, 0, 1);
+
+ // Momentum in beam and world frame for light particle 1
+ G4ThreeVector momentum_kine1_beam(sin(Theta1) * cos(Phi1),
+ sin(Theta1) * sin(Phi1), cos(Theta1));
+ G4ThreeVector momentum_kine1_world = momentum_kine1_beam;
+ momentum_kine1_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ momentum_kine1_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam and world frame for light particle 2
+ G4ThreeVector momentum_kine2_beam(sin(Theta2) * cos(Phi2),
+ sin(Theta2) * sin(Phi2), cos(Theta2));
+ G4ThreeVector momentum_kine2_world = momentum_kine2_beam;
+ momentum_kine2_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ momentum_kine2_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam and world frame for heavy residual
+ //
+ //G4ThreeVector momentum_kineB_beam(sin(ThetaB) * cos(PhiB + pi),
+ // sin(ThetaB) * sin(PhiB + pi), cos(ThetaB));
+ //G4ThreeVector momentum_kineB_world = momentum_kineB_beam;
+ //momentum_kineB_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ //momentum_kineB_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ TLorentzVector* P_B = m_QFS.GetEnergyImpulsionLab_B();
+
+ G4ThreeVector momentum_kineB_beam( P_B->Px(), P_B->Py(), P_B->Pz() );
+ momentum_kineB_beam = momentum_kineB_beam.unit();
+ double EnergyB = m_QFS.GetEnergyImpulsionLab_B()->Energy();
+ G4ThreeVector momentum_kineB_world = momentum_kineB_beam;
+ momentum_kineB_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ momentum_kineB_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
- // Nucleus 4
- G4int HeavyZ = m_Reaction.GetNucleus4()->GetZ();
- G4int HeavyA = m_Reaction.GetNucleus4()->GetA();
-
- // Generate the excitation energy if a distribution is given
- m_Reaction.ShootRandomExcitationEnergy();
-
- // Use to clean up the IonTable in case of the Ex changing at every event
- G4ParticleDefinition* HeavyName;
-
- if (m_Reaction.GetUseExInGeant4())
- HeavyName
- = IonTable->GetIon(HeavyZ, HeavyA, m_Reaction.GetExcitation4() * MeV);
- else
- HeavyName = IonTable->GetIon(HeavyZ, HeavyA);
-
- // Set the Energy of the reaction
- m_Reaction.SetBeamEnergy(energy);
-
- double Beam_theta = pdirection.theta();
- double Beam_phi = pdirection.phi();
-
- ///////////////////////////
- ///// Beam Parameters /////
- ///////////////////////////
- m_ReactionConditions->SetBeamParticleName(
- PrimaryTrack->GetParticleDefinition()->GetParticleName());
-
- m_ReactionConditions->SetBeamReactionEnergy(energy);
- m_ReactionConditions->SetVertexPositionX(localPosition.x());
- m_ReactionConditions->SetVertexPositionY(localPosition.y());
- m_ReactionConditions->SetVertexPositionZ(localPosition.z());
-
- G4ThreeVector U(1, 0, 0);
- G4ThreeVector V(0, 1, 0);
- G4ThreeVector ZZ(0, 0, 1);
- m_ReactionConditions->SetBeamEmittanceTheta(
- PrimaryTrack->GetMomentumDirection().theta() / deg);
- m_ReactionConditions->SetBeamEmittancePhi(
- PrimaryTrack->GetMomentumDirection().phi() / deg);
- m_ReactionConditions->SetBeamEmittanceThetaX(
- PrimaryTrack->GetMomentumDirection().angle(U) / deg);
- m_ReactionConditions->SetBeamEmittancePhiY(
- PrimaryTrack->GetMomentumDirection().angle(V) / deg);
-
- //////////////////////////////////////////////////////////
- ///// Build rotation matrix to go from the incident //////
- ///// beam frame to the "world" frame //////
- //////////////////////////////////////////////////////////
-
- /*
-
- G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
- cos(Beam_theta) * sin(Beam_phi),
- -sin(Beam_theta));
- G4ThreeVector col2(-sin(Beam_phi),
- cos(Beam_phi),
- 0);
- G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
- sin(Beam_theta) * sin(Beam_phi),
- cos(Beam_theta));
- G4RotationMatrix BeamToWorld(col1, col2, col3);
-
- */
-
- /////////////////////////////////////////////////////////////////
- ///// Angles for emitted particles following Cross Section //////
- ///// Angles are in the beam frame //////
- /////////////////////////////////////////////////////////////////
-
- // Angles
- // Shoot and Set a Random ThetaCM
- m_Reaction.ShootRandomThetaCM();
- double phi = RandFlat::shoot() * 2. * pi;
-
- //////////////////////////////////////////////////
- ///// Momentum and angles from kinematics /////
- ///// Angles are in the beam frame /////
- //////////////////////////////////////////////////
- // Variable where to store results
- double Theta3, Energy3, Theta4, Energy4;
-
- // Compute Kinematic using previously defined ThetaCM
- m_Reaction.KineRelativistic(Theta3, Energy3, Theta4, Energy4);
- // Momentum in beam frame for light particle
- G4ThreeVector momentum_kine3_beam(sin(Theta3) * cos(phi),
- sin(Theta3) * sin(phi), cos(Theta3));
- // Momentum in World frame //to go from the incident beam frame to the "world"
- // frame
- G4ThreeVector momentum_kine3_world = /*BeamToWorld */ momentum_kine3_beam;
- momentum_kine3_world.rotate(Beam_theta,
- V); // rotation of Beam_theta on Y axis
- momentum_kine3_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
-
- // Momentum in beam frame for heavy particle
- G4ThreeVector momentum_kine4_beam(sin(Theta4) * cos(phi + pi),
- sin(Theta4) * sin(phi + pi), cos(Theta4));
- // Momentum in World frame
- G4ThreeVector momentum_kine4_world = /*BeamToWorld */ momentum_kine4_beam;
- momentum_kine4_world.rotate(Beam_theta,
- V); // rotation of Beam_theta on Y axis
- momentum_kine4_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
-
- // Emitt secondary
- if (m_Reaction.GetShoot3()) {
- G4DynamicParticle particle3(LightName, momentum_kine3_world, Energy3);
- fastStep.CreateSecondaryTrack(particle3, localPosition, time);
- }
+ // Emitt secondary
+ if (m_QFS.GetShoot1()) {
+ G4DynamicParticle particle1(Light1Name, momentum_kine1_world, Energy1);
+ fastStep.CreateSecondaryTrack(particle1, localPosition, time);
+ }
- if (m_Reaction.GetShoot4()) {
- G4DynamicParticle particle4(HeavyName, momentum_kine4_world, Energy4);
- fastStep.CreateSecondaryTrack(particle4, localPosition, time);
- }
+ if (m_QFS.GetShoot2()) {
+ G4DynamicParticle particle2(Light2Name, momentum_kine2_world, Energy2);
+ fastStep.CreateSecondaryTrack(particle2, localPosition, time);
+ }
+ if (m_QFS.GetShootB()) {
+ G4DynamicParticle particleB(HeavyName, momentum_kineB_world, EnergyB);
+ fastStep.CreateSecondaryTrack(particleB, localPosition, time);
+ }
- // Reinit for next event
- m_PreviousEnergy = 0;
- m_PreviousLength = 0;
- ///////////////////////////////////////
- ///// Emitted particle Parameters /////
- ///////////////////////////////////////
- // Names 3 and 4//
- m_ReactionConditions->SetParticleName(LightName->GetParticleName());
- m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
- // Angle 3 and 4 //
- m_ReactionConditions->SetTheta(Theta3 / deg);
- m_ReactionConditions->SetTheta(Theta4 / deg);
-
- m_ReactionConditions->SetPhi(phi / deg);
- if ((phi + pi) / deg > 360)
- m_ReactionConditions->SetPhi((phi - pi) / deg);
- else
- m_ReactionConditions->SetPhi((phi + pi) / deg);
-
- // Energy 3 and 4 //
- m_ReactionConditions->SetKineticEnergy(Energy3);
- m_ReactionConditions->SetKineticEnergy(Energy4);
- // ThetaCM and Ex//
- m_ReactionConditions->SetThetaCM(m_Reaction.GetThetaCM() / deg);
- m_ReactionConditions->SetExcitationEnergy3(m_Reaction.GetExcitation3());
- m_ReactionConditions->SetExcitationEnergy4(m_Reaction.GetExcitation4());
- // Momuntum X 3 and 4 //
- m_ReactionConditions->SetMomentumDirectionX(momentum_kine3_world.x());
- m_ReactionConditions->SetMomentumDirectionX(momentum_kine4_world.x());
- // Momuntum Y 3 and 4 //
- m_ReactionConditions->SetMomentumDirectionY(momentum_kine3_world.y());
- m_ReactionConditions->SetMomentumDirectionY(momentum_kine4_world.y());
- // Momuntum Z 3 and 4 //
- m_ReactionConditions->SetMomentumDirectionZ(momentum_kine3_world.z());
- m_ReactionConditions->SetMomentumDirectionZ(momentum_kine4_world.z());
+ // Reinit for next event
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+
+
+
+ }
}
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index cb81602f18eb28cde9f4fcbd6af5d146f307e288..b77c424d8885112d5152c7b1f2bf32b37999e16d 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -27,6 +27,7 @@
#include "G4VFastSimulationModel.hh"
#include "PhysicsList.hh"
#include "NPReaction.h"
+#include "NPQFS.h"
#include "TReactionConditions.h"
class G4VPhysicalVolume;
namespace NPS{
@@ -44,7 +45,9 @@ namespace NPS{
private:
NPL::Reaction m_Reaction;
+ NPL::QFS m_QFS;
string m_BeamName;
+ string m_ReactionType;
double m_PreviousEnergy;
double m_PreviousLength;
bool m_active;// is the process active