Updating sofia simulation parameters

NPLib/Detectors/Sofia/GladFieldMap.cxx

@@ -158,7 +158,6 @@ double GladFieldMap::GetFlightPath(TVector3 vStart, double Brho, TVector3 Pos, T
   track = Propagate(Brho, Pos, Dir, true);
 
   FlightPath += (Pos - vStart).Mag();
-
   unsigned int vsize = track.size();
   for(unsigned int i=0; i<vsize-1; i++){
     TVector3 point1 = track[i];
@@ -404,7 +403,8 @@ void GladFieldMap::LoadMap(string filename) {
     }
   }
 
-  m_R_max = m_z_max;
+  m_R_max = m_x_max;
+  //m_R_max = m_z_max;
   cout << endl;
   cout << "///////// ASCII file loaded"<< endl;
   cout << "m_field size= " << m_By.size() << endl;

Projects/Sofia/Analysis.cxx

@@ -41,9 +41,7 @@ void Analysis::Init(){
   InitialConditions = new TInitialConditions();
   InteractionCoordinates = new TInteractionCoordinates();
   
-  /*if(beam_ana==false)
-    FissionConditions = new TFissionConditions();
-*/
+  FissionConditions = new TFissionConditions();
 
   InitInputBranch();
   InitOutputBranch();
@@ -69,16 +67,14 @@ void Analysis::TreatEvent(){
 
   ReInit();
 
-  if(beam_ana==true)
-    BeamAnalysis();
-  else 
-    FFAnalysis();
+  //BeamAnalysis();
+  FFAnalysis();
 }
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::FFAnalysis(){
 
-  /*int mult = InteractionCoordinates->GetDetectedMultiplicity();
+  int mult = InteractionCoordinates->GetDetectedMultiplicity();
   if(mult==2){
     int A_CN = FissionConditions->GetA_CN();
     int Z_CN = FissionConditions->GetZ_CN();
@@ -92,10 +88,10 @@ void Analysis::FFAnalysis(){
     double beta = beam->GetBeta();
     double BeamTimeOffset = v1.Mag()*mm/(beta*NPUNITS::c_light);
     for(int i=0; i<mult; i++){
-      double Time = ran.Gaus(InteractionCoordinates->GetTime(i) - BeamTimeOffset,0.);
-      m_XE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionX(i),0));
+      double Time = ran.Gaus(InteractionCoordinates->GetTime(i) - BeamTimeOffset,0.02);
+      m_XE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionX(i),1));
       m_YE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionY(i),0));
-      m_ZE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(i),0));
+      m_ZE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(i),1));
       TVector3 vE = TVector3(m_XE[i],m_YE[i],m_ZE[i]);
       //TVector3 vE = TVector3(m_XE[i],0,m_ZE[i]);
       TVector3 vA = TVector3(0,0,z_target);
@@ -104,8 +100,8 @@ void Analysis::FFAnalysis(){
       //int Z = FissionConditions->GetFragmentZ(i);
       int Z = InteractionCoordinates->GetZ(i);
       double ThetaIn = FissionConditions->GetFragmentTheta(i);
-      if(FissionConditions->GetFragmentMomentumX(i)<0)
-        ThetaIn = -ThetaIn;
+      //if(FissionConditions->GetFragmentMomentumX(i)<0)
+        //ThetaIn = -ThetaIn;
       double Phi = FissionConditions->GetFragmentPhi(i);
       double Brho = FissionConditions->GetFragmentBrho(i);
  
@@ -114,6 +110,7 @@ void Analysis::FFAnalysis(){
       m_TOF.push_back(Time);
       m_Adet.push_back(Adet);
       m_Asim.push_back(Asim);
+      m_Zsim.push_back(Z);
       m_Brho.push_back(Brho);
       m_ThetaIn.push_back(ThetaIn);
  
@@ -124,14 +121,13 @@ void Analysis::FFAnalysis(){
       double beta = velocity * m/ns / NPUNITS::c_light;
       double gamma = 1. / sqrt(1 - beta*beta);
       double AoQ = Brho_calc / (3.107 * beta * gamma);
-      double A_calc = AoQ*Z;
+      double A_calc = AoQ*(double)Z;
 
       m_Brho_calc.push_back(Brho_calc);
       m_FlightPath.push_back(PathLength);
       m_A_calc.push_back(A_calc);
     }
-  }*/
-    
+  }
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -192,11 +188,9 @@ void Analysis::InitInputBranch(){
   RootInput::getInstance()->GetChain()->SetBranchStatus("fDetecetd_*",true);
   RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&InteractionCoordinates);
 
-  /*if(beam_ana==false){
-    RootInput::getInstance()->GetChain()->SetBranchStatus("FissionConditions",true);
-    RootInput::getInstance()->GetChain()->SetBranchStatus("fFC_*",true);
-    RootInput::getInstance()->GetChain()->SetBranchAddress("FissionConditions",&FissionConditions);
-  }*/
+  RootInput::getInstance()->GetChain()->SetBranchStatus("FissionConditions",true);
+  RootInput::getInstance()->GetChain()->SetBranchStatus("fFC_*",true);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("FissionConditions",&FissionConditions);
 
 }
 
@@ -206,6 +200,7 @@ void Analysis::InitOutputBranch(){
   RootOutput::getInstance()->GetTree()->Branch("ThetaIn",&m_ThetaIn);
   RootOutput::getInstance()->GetTree()->Branch("Brho",&m_Brho);
   RootOutput::getInstance()->GetTree()->Branch("Asim",&m_Asim);
+  RootOutput::getInstance()->GetTree()->Branch("Zsim",&m_Zsim);
   RootOutput::getInstance()->GetTree()->Branch("Adet",&m_Adet);
   RootOutput::getInstance()->GetTree()->Branch("Brho_calc",&m_Brho_calc);
   RootOutput::getInstance()->GetTree()->Branch("A_calc",&m_A_calc);
@@ -221,6 +216,7 @@ void Analysis::ReInit(){
   m_TOF.clear();
   m_Brho.clear();
   m_Asim.clear();
+  m_Zsim.clear();
   m_Adet.clear();
   m_ThetaIn.clear();
   m_XE.clear();

Projects/Sofia/Analysis.h

@@ -50,9 +50,10 @@ class Analysis: public NPL::VAnalysis{
     vector<double> m_ThetaIn;
     vector<double> m_Brho;
     vector<int> m_Asim;
+    vector<int> m_Zsim;
     vector<int> m_Adet;
     vector<double> m_Brho_calc;
-    vector<int> m_A_calc;
+    vector<double> m_A_calc;
     vector<double> m_FlightPath;
     vector<double> m_XE;
     vector<double> m_YE;
@@ -61,7 +62,7 @@ class Analysis: public NPL::VAnalysis{
   private:
     TInitialConditions* InitialConditions;
     TInteractionCoordinates* InteractionCoordinates;
-    //TFissionConditions* FissionConditions;
+    TFissionConditions* FissionConditions;
 
     TRandom3 ran;
     GladFieldMap* m_GladField;

Projects/Sofia/sofia.detector

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%Target
+Target
  THICKNESS= 0.03 mm
  RADIUS= 100 mm
  MATERIAL= CD2

Projects/Sofia/sofia.reaction

@@ -6,18 +6,18 @@ Beam
  Particle= 190Pb
  Energy= 100276
  ExcitationEnergy= 0 MeV
- SigmaEnergy= 10000
+ SigmaEnergy= 0
  SigmaX= 0 mm
  SigmaY= 0 mm
- SigmaThetaX= 0 
- SigmaPhiY= 0
+ SigmaThetaX= 0 deg
+ SigmaPhiY= 0 deg
  MeanThetaX= 0
  MeanPhiY= 0 deg
  MeanX= 0 mm
  MeanY= 0 mm
  ZEmission= -5000 mm
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%TwoBodyReaction
+TwoBodyReaction
  Beam= 190Pb
  Target= 208Pb
  Light= 190Pb
@@ -29,7 +29,7 @@ Beam
  ShootHeavy= 1
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 FissionDecay 
- CompoundNucleus= 191Pb
+ CompoundNucleus= 190Pb
  FissionModel= GEF
  VamosChargeStates= 0
  Shoot_FF= 1