Advancing on strasse missing mass

NPLib/Physics/NPQFS.cxx

@@ -215,7 +215,7 @@ void QFS::CalculateVariables(){
     //cout<<"---- COMPUTE ------"<<endl;
    // cout<<"--CM--"<<endl; 
 
-    mA = fNucleiA.Mass();            // Beam mass in MeV
+    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);
@@ -307,7 +307,7 @@ void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEn
                                         pCM_2*sin(thetaCM_2)*sin(phiCM_2),
                                         pCM_2*cos(thetaCM_2),
                                         ECM_2);
-
+  
     fEnergyImpulsionCM_1	= fTotalEnergyImpulsionCM - fEnergyImpulsionCM_2;
 
     //-- Boost in the direction of the moving cluster "a" --//
@@ -346,7 +346,6 @@ void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEn
     // 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;

NPLib/Physics/TReactionConditions.h

@@ -104,7 +104,12 @@ public:
     void SetMomentumDirectionX (const double & Momentum_Direction_X)  {fRC_Momentum_Direction_X.push_back(Momentum_Direction_X);}//!
     void SetMomentumDirectionY (const double & Momentum_Direction_Y)  {fRC_Momentum_Direction_Y.push_back(Momentum_Direction_Y);}//!
     void SetMomentumDirectionZ (const double & Momentum_Direction_Z)  {fRC_Momentum_Direction_Z.push_back(Momentum_Direction_Z);}//!
-    void SetMomentum (const TVector3 & Momentum)  {fRC_Momentum.push_back(Momentum);}//!
+    void SetMomentum (const TVector3 & Momentum)  {
+      Momentum.Unit();
+      SetMomentumDirectionX(Momentum.X());
+      SetMomentumDirectionY(Momentum.Y());
+      SetMomentumDirectionZ(Momentum.Z());
+    }//!
     
     /////////////////////           GETTERS           ////////////////////////
     // Beam parameter
@@ -133,7 +138,7 @@ public:
     double GetMomentumDirectionX  (const int &i) const {return fRC_Momentum_Direction_X[i];}//!
     double GetMomentumDirectionY  (const int &i) const {return fRC_Momentum_Direction_Y[i];}//!
     double GetMomentumDirectionZ  (const int &i) const {return fRC_Momentum_Direction_Z[i];}//!
-    TVector3 GetParticleMomentum  (const int &i) const {return fRC_Momentum[i];}//!
+    TVector3 GetParticleMomentum  (const int &i) const {return TVector3(fRC_Momentum_Direction_X[i],fRC_Momentum_Direction_Y[i],fRC_Momentum_Direction_Z[i]).Unit();}//!
 
     TVector3 GetBeamDirection         () const ;
     TVector3 GetParticleDirection     (const int i) const ; 

NPSimulation/Process/BeamReaction.cc

@@ -482,10 +482,10 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
 
     TLorentzVector* P_A = m_QFS.GetEnergyImpulsionLab_A();
     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();
-    TKEB = m_QFS.GetEnergyImpulsionLab_B()->Energy() - m_QFS.GetNucleusB()->Mass();
+    TKEB = P_B->Energy() - m_QFS.GetNucleusB()->Mass();
     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
@@ -549,15 +549,15 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
     m_ReactionConditions->SetInternalMomentum(m_QFS.GetInternalMomentum());
     //m_ReactionConditions->SetExcitationEnergy3(m_QFS.GetExcitation3());
     //m_ReactionConditions->SetExcitationEnergy4(m_QFS.GetExcitation4());
-    // Momuntum X 3 and 4 //
+    // Momuntum X 1,2 and B //
     m_ReactionConditions->SetMomentumDirectionX(momentum_kine1_world.x());
     m_ReactionConditions->SetMomentumDirectionX(momentum_kine2_world.x());
     m_ReactionConditions->SetMomentumDirectionX(momentum_kineB_world.x());
-    // Momuntum Y 3 and 4 //
+    // Momuntum Y 1,2 and B //
     m_ReactionConditions->SetMomentumDirectionY(momentum_kine1_world.y());
     m_ReactionConditions->SetMomentumDirectionY(momentum_kine2_world.y());
     m_ReactionConditions->SetMomentumDirectionY(momentum_kineB_world.y());
-    // Momuntum Z 3 and 4 //
+    // Momuntum Z 1,2 and B //
     m_ReactionConditions->SetMomentumDirectionZ(momentum_kine1_world.z());
     m_ReactionConditions->SetMomentumDirectionZ(momentum_kine2_world.z());
     m_ReactionConditions->SetMomentumDirectionZ(momentum_kineB_world.z());

Projects/Strasse/Analysis.cxx

@@ -27,6 +27,7 @@ using namespace std;
 #include"NPOptionManager.h"
 #include"NPFunction.h"
 #include"NPTrackingUtility.h"
+#include"NPPhysicalConstants.h"
 ////////////////////////////////////////////////////////////////////////////////
 Analysis::Analysis(){
 }
@@ -40,19 +41,18 @@ void Analysis::Init(){
   DC= new TInteractionCoordinates;
   RC= new TReactionConditions;
   
-
   InitOutputBranch();
   InitInputBranch();
   
   Strasse = (TStrassePhysics*)  m_DetectorManager -> GetDetector("Strasse");
   Catana = (TCatanaPhysics*)  m_DetectorManager -> GetDetector("Catana");
   // reaction properties
-  myQFS = new NPL::QFS();
-  myQFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+  m_QFS = new NPL::QFS();
+  m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
   // reaction properties
   myBeam = new NPL::Beam();
   myBeam->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
-  InitialBeamEnergy = myBeam->GetEnergy() * myBeam->GetA();
+  InitialBeamEnergy = myBeam->GetEnergy()/* * myBeam->GetA()*/;
   // target thickness
   TargetThickness = m_DetectorManager->GetTargetThickness();
   string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -62,8 +62,7 @@ void Analysis::Init(){
   protonTarget = NPL::EnergyLoss("proton_"+TargetMaterial+".G4table","G4Table",100);
   protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100);
   protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100);
-
-
+  LV_T.SetVectM(TVector3(0,0,0),NPUNITS::proton_mass_c2);
 } 
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -109,20 +108,28 @@ void Analysis::TreatEvent(){
     if(i1!=i2){
       E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]); 
       E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]);
+      double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),0);
+      //double TA = RC->GetBeamEnergy();
+      // setting up Lorentz Vector from measured trajectories and energies
+      TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetNucleusA()->Mass()))); // for like there is no BDC
+      Proton1=E1*Proton1.Unit();
+      Proton2=E2*Proton2.Unit();
+      
+      LV_A.SetVectM(PA,m_QFS->GetNucleusA()->Mass());
+      double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2));
+      double P2= sqrt(E2*(E2+2*NPUNITS::proton_mass_c2));
+
+      LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2); 
+      LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2); 
+
+      // computing Ex from Missing Mass
+      LV_B = LV_A + LV_T - LV_p1 - LV_p2;
+      //LV_B = RC->GetParticleMomentum(2);
+      Ex = LV_B.M() - m_QFS->GetNucleusB()->Mass();
+      cout << Ex << " " << m_QFS->GetNucleusB()->Mass() << endl;
     }
+    
   }
-    //double thickness_before = 0;
-    //double EA_vertex = BeamTarget.Slow(InitialBeamEnergy,thickness_before,0);
-
-    // setting up Lorentz Vector from measured trajectories and energies
-    //LV_A.SetVect(PA); LV_p1.SetE(EA_vertex); 
-    //LV_p1.SetVect(P1); LV_p1.SetE(E1); 
-    //LV_p2.SetVect(P2); LV_p1.SetE(E2); 
-
-    // computing Ex from Missing Mass
-    //double EB = LV_A.E() + LV_T.E() - LV_p1.E() - LV_p2.E();   
-    //TVector3 PB = LV_A.Vect() + LV_p1.Vect() - LV_p2.Vect();   
-    //Ex = TMath::Sqrt( EB*EB - PB.Mag2() ) - myQFS->GetNucleusB()->Mass();
 }
 
 ////////////////////////////////////////////////////////////////////////////////

Projects/Strasse/Analysis.h

@@ -73,7 +73,7 @@ class Analysis: public NPL::VAnalysis{
     TLorentzVector LV_p2;
 
     NPL::Beam* myBeam;
-    NPL::QFS* myQFS;
+    NPL::QFS* m_QFS;
     //	Energy loss table: the G4Table are generated by the simulation
     EnergyLoss BeamTarget;
     EnergyLoss protonTarget;

Projects/Strasse/reaction/C12_p2p.reaction

@@ -4,10 +4,10 @@ Beam
   Particle= 12C
   Energy= 4800 MeV
   SigmaEnergy= 0 MeV
-  SigmaThetaX= 0.1 deg
-  SigmaPhiY= 0.1 deg
-  SigmaX= 5 mm
-  SigmaY= 5 mm
+  SigmaThetaX= 0 deg
+  SigmaPhiY= 0 deg
+  SigmaX= 0 mm
+  SigmaY= 0 mm
   MeanThetaX= 0 deg
   MeanPhiY= 0 deg
   MeanX= 0 mm
@@ -22,7 +22,7 @@ QFSReaction
  Heavy= 11B
  ExcitationEnergyBeam= 0.0 MeV
  ExcitationEnergyHeavy= 0.0 MeV
- MomentumSigma= 50.0 
+ MomentumSigma= 0.0 
  ShootHeavy= 1
  ShootLight= 1