From 406247ef1661aa1ea1535d39a975efa24c03b035 Mon Sep 17 00:00:00 2001
From: Morfouace <pierre.morfouace@cea.fr>
Date: Wed, 25 Jan 2023 17:01:44 +0100
Subject: [PATCH] Updating sofia simulation
---
NPLib/Detectors/Sofia/GladFieldMap.cxx | 18 +-
NPSimulation/Detectors/Sofia/Glad.cc | 51 +-
NPSimulation/Detectors/Sofia/Glad.hh | 9 +-
.../Detectors/Sofia/GladFieldPropagation.cc | 3 +-
.../Detectors/Sofia/GladFieldPropagation.hh | 4 +
NPSimulation/Process/FissionDecay.cc | 5 +-
NPSimulation/Scorers/InteractionScorers.cc | 3 +-
Projects/Sofia/Analysis.cxx | 135 +-
Projects/Sofia/Analysis.h | 13 +-
Projects/Sofia/GladFieldMap_50mm.dat | 401845 +++++++++++++++
Projects/Sofia/sim.sh | 7 +
Projects/Sofia/sim_exo.sh | 7 +
Projects/Sofia/sofia.detector | 10 +-
Projects/Sofia/sofia.reaction | 37 +
14 files changed, 402064 insertions(+), 83 deletions(-)
create mode 100644 Projects/Sofia/GladFieldMap_50mm.dat
create mode 100755 Projects/Sofia/sim.sh
create mode 100755 Projects/Sofia/sim_exo.sh
create mode 100644 Projects/Sofia/sofia.reaction
diff --git a/NPLib/Detectors/Sofia/GladFieldMap.cxx b/NPLib/Detectors/Sofia/GladFieldMap.cxx
index 0566d486b..df7f7ce7a 100644
--- a/NPLib/Detectors/Sofia/GladFieldMap.cxx
+++ b/NPLib/Detectors/Sofia/GladFieldMap.cxx
@@ -101,7 +101,7 @@ double GladFieldMap::Delta(const double* parameter){
double GladFieldMap::FindBrho(TVector3 Pos_init, TVector3 Dir_init, TVector3 Pos_final){
if(!m_BrhoScan)
- BrhoScan(6,11,0.1,TVector3(0,0,-2500), TVector3(0,0,1));
+ BrhoScan(6,11,0.1,TVector3(0,0,-5000), TVector3(0,0,1));
m_InitPos = Pos_init;
m_InitDir = Dir_init;
@@ -177,21 +177,21 @@ TVector3 GladFieldMap::PropagateToMWPC(TVector3 pos, TVector3 dir){
pos.RotateY(-m_CentralTheta);
dir.RotateY(-m_CentralTheta);
- double deltaZ = m_MWPC3_POS.Z() - pos.Z();
+ double deltaZ = m_MWPC3_POS.Mag() - pos.Z();
dir*=deltaZ/dir.Z();
pos+=dir;
pos.SetX(pos.X());
pos.RotateY(m_CentralTheta);
-
+
return pos;
}
//////////////////////////////////////////////////////////////////////
vector<TVector3> GladFieldMap::Propagate(double Brho, TVector3 Pos, TVector3 Dir, bool store){
+
Pos.RotateY(m_Tilt);
Dir.RotateY(m_Tilt);
Dir = Dir.Unit();
-
static NPL::Particle N("90Zr");
N.SetBrho(Brho);
@@ -256,7 +256,6 @@ vector<TVector3> GladFieldMap::Propagate(double Brho, TVector3 Pos, TVector3 Dir
Pos += (m_dt/6)*(xk1 + 2*xk2 + 2*xk3 + xk4);
Imp += (m_dt/6)*(pk1 + 2*pk2 + 2*pk3 + pk4);
-
if(store){
Pos.RotateY(-m_Tilt);
track.push_back(Pos);
@@ -304,6 +303,7 @@ void GladFieldMap::func(NPL::Particle& N, TVector3 Pos, TVector3 Imp, TVector3&
By = GetB(Pos,"Y");
Bz = GetB(Pos,"Z");
+
q = N.GetZ()*eplus;
pk.SetX(q*(vy*Bz - vz*By));
pk.SetY(q*(vz*Bx - vx*Bz));
@@ -404,7 +404,7 @@ void GladFieldMap::LoadMap(string filename) {
}
}
- 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;
@@ -544,10 +544,10 @@ vector<double> GladFieldMap::InterpolateB(const vector<double>& pos)
//////////////////////////////////
double GladFieldMap::GetB(TVector3 localpoint, string field_component)
{
- TVector3 vtrans(0,0,-m_Glad_Entrance.Z());
-
- localpoint = localpoint + vtrans;
+ TVector3 vtrans(m_Glad_Entrance.X(),m_Glad_Entrance.Y(),m_Glad_Entrance.Z());
+ localpoint = localpoint - vtrans;
+
static int ix, iy, iz;
static double dx, dy, dz;
double val = 0;
diff --git a/NPSimulation/Detectors/Sofia/Glad.cc b/NPSimulation/Detectors/Sofia/Glad.cc
index eb194a33f..812b65810 100644
--- a/NPSimulation/Detectors/Sofia/Glad.cc
+++ b/NPSimulation/Detectors/Sofia/Glad.cc
@@ -88,6 +88,8 @@ Glad::Glad(){
m_GLAD_Y = 0;
m_GLAD_Z = -1113.5*mm;
+ m_MW3_R = 4440*mm;
+
// RGB Color + Transparency
m_VisSquare = new G4VisAttributes(G4Colour(0.53, 0.81, 0.98, 1));
m_VisGLAD = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.9));
@@ -99,28 +101,18 @@ Glad::Glad(){
Glad::~Glad(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Glad::AddMagnet(G4ThreeVector POS, double Tilt_Angle, string fieldmap){
+void Glad::AddMagnet(G4ThreeVector POS, string fieldmap){
// Convert the POS value to R theta Phi as Spherical coordinate is easier in G4
+ m_X = POS.x();
+ m_Y = POS.y();
+ m_Z = POS.z();
m_R = POS.mag();
m_Theta = POS.theta();
m_Phi = POS.phi();
- m_GLAD_TiltAngle = Tilt_Angle;
- m_FieldMapFile = fieldmap;
-
- return;
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Glad::AddMagnet(double R, double Theta, double Phi, double Tilt_Angle, string fieldmap){
- m_R = R;
- m_Theta = Theta;
- m_Phi = Phi;
-
- m_GLAD_TiltAngle= Tilt_Angle;
m_FieldMapFile = fieldmap;
+ m_MW3_POS = TVector3(-m_MW3_R*sin(20*deg),0,m_MW3_R*cos(20*deg));
return;
}
@@ -174,7 +166,6 @@ void Glad::ReadConfiguration(NPL::InputParser parser){
cout << "//// GLAD Magnet found ////" << endl;
vector<string> cart = {"POS","TILT_ANGLE","FIELDMAP","CURRENT","STEPSIZE"};
- vector<string> sphe = {"R","Theta","Phi","TILT_ANGLE","FIELDMAP","CURRENT","STEPSIZE"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
@@ -186,22 +177,9 @@ void Glad::ReadConfiguration(NPL::InputParser parser){
string fieldmap = blocks[i]->GetString("FIELDMAP");
m_Current = blocks[i]->GetDouble("CURRENT","A");
m_StepSize = blocks[i]->GetDouble("STEPSIZE","mm");
+ m_MW3_R = blocks[i]->GetDouble("MW3_DISTANCE","mm");
- AddMagnet(Pos,m_GLAD_TiltAngle,fieldmap);
- }
- else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Glad information: " << endl;
-
- double R = blocks[i]->GetDouble("R","mm");
- double Theta = blocks[i]->GetDouble("Theta","deg");
- double Phi = blocks[i]->GetDouble("Phi","deg");
- m_GLAD_TiltAngle = blocks[i]->GetDouble("TILT_ANGLE","deg");
- string fieldmap = blocks[i]->GetString("FIELDMAP");
- m_Current = blocks[i]->GetDouble("CURRENT","A");
- m_StepSize = blocks[i]->GetDouble("STEPSIZE","mm");
-
- AddMagnet(R,Theta,Phi,m_GLAD_TiltAngle,fieldmap);
+ AddMagnet(Pos,fieldmap);
}
else{
cout << "ERROR: check your GLAD input file formatting " << endl;
@@ -217,9 +195,12 @@ void Glad::ReadConfiguration(NPL::InputParser parser){
// Called After DetecorConstruction::AddDetector Method
void Glad::ConstructDetector(G4LogicalVolume* world){
- G4double wX = m_R * sin(m_Theta) * cos(m_Phi);
- G4double wY = m_R * sin(m_Theta) * sin(m_Phi);
- G4double wZ = m_R * cos(m_Theta);
+ //G4double wX = m_R * sin(m_Theta) * cos(m_Phi);
+ //G4double wY = m_R * sin(m_Theta) * sin(m_Phi);
+ //G4double wZ = m_R * cos(m_Theta);
+ G4double wX = m_X - 1500*sin(m_GLAD_TiltAngle);
+ G4double wY = m_Y;
+ G4double wZ = m_Z - 1500*(1-cos(m_GLAD_TiltAngle));
G4ThreeVector Mag_pos = G4ThreeVector(wX,wY,wZ);
G4ThreeVector u (cos(m_GLAD_TiltAngle), 0, -sin(m_GLAD_TiltAngle));
@@ -265,9 +246,11 @@ void Glad::SetPropagationRegion()
G4VFastSimulationModel* fsm;
fsm = new NPS::GladFieldPropagation("GladFieldPropagation", m_PropagationRegion);
((NPS::GladFieldPropagation*) fsm)->SetGladEntrance(m_GLAD_X, m_GLAD_Y, m_GLAD_Z);
+ ((NPS::GladFieldPropagation*) fsm)->SetGladTiltAngle(m_GLAD_TiltAngle);
((NPS::GladFieldPropagation*) fsm)->SetCurrent(m_Current);
((NPS::GladFieldPropagation*) fsm)->SetStepSize(m_StepSize);
((NPS::GladFieldPropagation*) fsm)->SetFieldMap(m_FieldMapFile);
+ ((NPS::GladFieldPropagation*) fsm)->Set_MWPC3_Position(m_MW3_POS.X(),m_MW3_POS.Y(),m_MW3_POS.Z());
/*if(m_Method == NPS::RungeKutta){
double r_max = sqrt(
Glad_NS::GLAD_Width * Glad_NS::GLAD_Width /4. +
diff --git a/NPSimulation/Detectors/Sofia/Glad.hh b/NPSimulation/Detectors/Sofia/Glad.hh
index 1c3f215f9..2129322fb 100644
--- a/NPSimulation/Detectors/Sofia/Glad.hh
+++ b/NPSimulation/Detectors/Sofia/Glad.hh
@@ -54,9 +54,7 @@ class Glad : public NPS::VDetector{
////////////////////////////////////////////////////
public:
// Cartesian
- void AddMagnet(G4ThreeVector POS, double Tilt_Angle, string fieldmap);
- // Spherical
- void AddMagnet(double R,double Theta,double Phi, double Tilt_Angle, string fieldmap);
+ void AddMagnet(G4ThreeVector POS, string fieldmap);
G4LogicalVolume* BuildGLADFromSTL();
G4LogicalVolume* BuildMagnet();
@@ -106,6 +104,9 @@ class Glad : public NPS::VDetector{
////////////////////////////////////////////////////
private: // Geometry
// Detector Coordinate
+ double m_X;
+ double m_Y;
+ double m_Z;
double m_R;
double m_Theta;
double m_Phi;
@@ -118,6 +119,8 @@ class Glad : public NPS::VDetector{
double m_Current;
int m_StepSize;
string m_FieldMapFile;
+ TVector3 m_MW3_POS;
+ double m_MW3_R;
// Visualisation Attribute
G4VisAttributes* m_VisSquare;
diff --git a/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc b/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc
index 55d366fd6..74897da78 100644
--- a/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc
+++ b/NPSimulation/Detectors/Sofia/GladFieldPropagation.cc
@@ -93,7 +93,8 @@ void NPS::GladFieldPropagation::DoIt(const G4FastTrack& fastTrack, G4FastStep& f
m_Map->SetLimit( 10.*m / m_StepSize );
m_Map->SetPropagationTimeInterval(0.8);
m_Map->SetGladEntrance(m_GladEntrance.X(), m_GladEntrance.Y(), m_GladEntrance.Z());
-
+ m_Map->SetGladTiltAngle(m_GladTiltAngle);
+ m_Map->Set_MWPC3_Position(m_MW3_POS.X(), m_MW3_POS.Y(), m_MW3_POS.Z());
m_Map->LoadMap(m_FieldMap);
m_Initialized = true;
}
diff --git a/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh b/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh
index 7e833a2d8..777c19b6c 100644
--- a/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh
+++ b/NPSimulation/Detectors/Sofia/GladFieldPropagation.hh
@@ -57,7 +57,9 @@ namespace NPS{
double m_StepSize;
double m_Rmax;
double m_Current;
+ double m_GladTiltAngle;
TVector3 m_GladEntrance;
+ TVector3 m_MW3_POS;
string m_FieldMap;
GladFieldMap* m_Map;
@@ -70,7 +72,9 @@ namespace NPS{
void SetFieldMap(string fieldmap){m_FieldMap=fieldmap;};
void SetRmax(double r_max){m_Rmax=r_max;};
void SetCurrent(double val){m_Current=val;}
+ void SetGladTiltAngle(double val){m_GladTiltAngle=val;}
void SetGladEntrance(double x, double y, double z){m_GladEntrance=TVector3(x,y,z);}
+ void Set_MWPC3_Position(double x, double y, double z){m_MW3_POS=TVector3(x,y,z);}
};
}
diff --git a/NPSimulation/Process/FissionDecay.cc b/NPSimulation/Process/FissionDecay.cc
index e87765164..cf89b91f8 100644
--- a/NPSimulation/Process/FissionDecay.cc
+++ b/NPSimulation/Process/FissionDecay.cc
@@ -183,12 +183,13 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
if(Zsum == m_CompoundParticle.GetZ())
m_FissionConditions->SetNeutronMultiplicity(m_CompoundParticle.GetA()-Asum);
- for(int i=size-1; i>=0; --i){
+ //for(int i=size-1; i>=0; --i){
+ for(unsigned int i=0; i<size; i++){
// Get the decaying particle
int FFZ = FissionFragment[i].GetZ();
int FFA = FissionFragment[i].GetA();
FissionFragmentDef=NULL;
-
+
// Set the momentum direction
G4ThreeVector Momentum (DPx[i],DPy[i],DPz[i]);
Momentum=Momentum.unit();
diff --git a/NPSimulation/Scorers/InteractionScorers.cc b/NPSimulation/Scorers/InteractionScorers.cc
index 4f0966282..6eddde1bb 100644
--- a/NPSimulation/Scorers/InteractionScorers.cc
+++ b/NPSimulation/Scorers/InteractionScorers.cc
@@ -81,7 +81,8 @@ void PS_Interactions::Initialize(G4HCofThisEvent*){
void PS_Interactions::EndOfEvent(G4HCofThisEvent*){
unsigned int size = m_DataVector.size();
- for(unsigned int i = 0 ; i < size ; i++){
+ //for(unsigned int i = 0 ; i < size ; i++){
+ for(int i = size-1 ; i>=0 ; --i){
int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
int Index = t_Index + eventID;
m_InterractionCoordinates->SetInteraction(Index,m_DataVector[i]->GetEnergy(),m_DataVector[i]->GetTime(),
diff --git a/Projects/Sofia/Analysis.cxx b/Projects/Sofia/Analysis.cxx
index f2cb83078..f01f92f7b 100644
--- a/Projects/Sofia/Analysis.cxx
+++ b/Projects/Sofia/Analysis.cxx
@@ -37,24 +37,30 @@ Analysis::~Analysis(){
void Analysis::Init(){
//Sofia= (TSofiaPhysics*) m_DetectorManager->GetDetector("Sofia");
+ beam_ana = true;
InitialConditions = new TInitialConditions();
InteractionCoordinates = new TInteractionCoordinates();
- FissionConditions = new TFissionConditions();
+
+ /*if(beam_ana==false)
+ FissionConditions = new TFissionConditions();
+*/
InitInputBranch();
InitOutputBranch();
m_GladField = new GladFieldMap();
m_GladField->SetCurrent(2185);
- m_GladField->SetGladEntrance(0,0,-1.1135*m);
- m_GladField->SetGladTiltAngle(-14*deg);
+ m_GladField->SetGladEntrance(0,0,-1113.5*mm);
+ m_GladField->SetGladTiltAngle(-0*deg);
m_GladField->LoadMap("GladFieldMap_50mm.dat");
m_GladField->SetBin(50);
m_GladField->SetTimeStep(0.8);
m_GladField->SetCentralTheta(-20*deg);
- double Z_MW3 = 4440 * cos(20*deg);
- double X_MW3 = Z_MW3 * tan(20*deg);
+ double R_MW3 = 4440*mm;
+ double Z_MW3 = R_MW3 * cos(20*deg);
+ double X_MW3 = -R_MW3 * sin(20*deg);
+
m_GladField->Set_MWPC3_Position(X_MW3,0,Z_MW3);
}
@@ -63,7 +69,16 @@ void Analysis::TreatEvent(){
ReInit();
- int mult = InteractionCoordinates->GetDetectedMultiplicity();
+ if(beam_ana==true)
+ BeamAnalysis();
+ else
+ FFAnalysis();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::FFAnalysis(){
+
+ /*int mult = InteractionCoordinates->GetDetectedMultiplicity();
if(mult==2){
int A_CN = FissionConditions->GetA_CN();
int Z_CN = FissionConditions->GetZ_CN();
@@ -77,25 +92,34 @@ void Analysis::TreatEvent(){
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.02);
- double XE = ran.Gaus(InteractionCoordinates->GetDetectedPositionX(i),0.1);
- double YE = 0;//ran.Gaus(InteractionCoordinates->GetDetectedPositionY(i),0.1);
- double ZE = ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(i),0.1);
- TVector3 vE = TVector3(XE,YE,ZE);
+ double Time = ran.Gaus(InteractionCoordinates->GetTime(i) - BeamTimeOffset,0.);
+ m_XE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionX(i),0));
+ m_YE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionY(i),0));
+ m_ZE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(i),0));
+ 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);
- int A = FissionConditions->GetFragmentA(i);
- int Z = FissionConditions->GetFragmentZ(i);
- double Theta = FissionConditions->GetFragmentTheta(i);
+ int Asim = FissionConditions->GetFragmentA(i);
+ int Adet = InteractionCoordinates->GetA(i);
+ //int Z = FissionConditions->GetFragmentZ(i);
+ int Z = InteractionCoordinates->GetZ(i);
+ double ThetaIn = FissionConditions->GetFragmentTheta(i);
+ if(FissionConditions->GetFragmentMomentumX(i)<0)
+ ThetaIn = -ThetaIn;
double Phi = FissionConditions->GetFragmentPhi(i);
double Brho = FissionConditions->GetFragmentBrho(i);
- TVector3 dir = TVector3(sin(Theta*deg)*cos(Phi*deg), sin(Theta*deg)*sin(Phi*deg), cos(Theta*deg));
+ TVector3 dir = TVector3(sin(ThetaIn*deg)*cos(Phi*deg), sin(ThetaIn*deg)*sin(Phi*deg), cos(ThetaIn*deg));
+ //TVector3 dir = TVector3(sin(ThetaIn*deg), 0, cos(ThetaIn*deg));
m_TOF.push_back(Time);
- m_A.push_back(A);
+ m_Adet.push_back(Adet);
+ m_Asim.push_back(Asim);
m_Brho.push_back(Brho);
-
+ m_ThetaIn.push_back(ThetaIn);
+
+
double Brho_calc = m_GladField->FindBrho(vA,dir,vE);
- double PathLength = m_GladField->GetFlightPath(vA,Brho_calc,vA,dir);
+ double PathLength = m_GladField->GetFlightPath(vA,Brho_calc,vA,dir)/1000;
double velocity = PathLength/Time;
double beta = velocity * m/ns / NPUNITS::c_light;
double gamma = 1. / sqrt(1 - beta*beta);
@@ -104,12 +128,60 @@ void Analysis::TreatEvent(){
m_Brho_calc.push_back(Brho_calc);
m_FlightPath.push_back(PathLength);
- m_A.push_back(A_calc);
+ m_A_calc.push_back(A_calc);
}
- }
+ }*/
}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::BeamAnalysis(){
+
+ int mult = InteractionCoordinates->GetDetectedMultiplicity();
+ if(mult==1){
+ double Elab_CN = InitialConditions->GetIncidentInitialKineticEnergy();
+ double x_start_beam = InitialConditions->GetIncidentPositionX();
+ double y_start_beam = InitialConditions->GetIncidentPositionY();
+ double z_start_beam = InitialConditions->GetIncidentPositionZ();
+
+ int Adet = InteractionCoordinates->GetA(0);
+ int Zdet = InteractionCoordinates->GetZ(0);
+ NPL::Particle* beam = new NPL::Particle(Zdet,Adet);
+ beam->SetKineticEnergy(Elab_CN);
+ //double beta = beam->GetBeta();
+ double Brho = beam->GetBrho();
+ double Time = ran.Gaus(InteractionCoordinates->GetTime(0),0.);
+ m_XE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionX(0),0));
+ m_YE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionY(0),0));
+ m_ZE.push_back(ran.Gaus(InteractionCoordinates->GetDetectedPositionZ(0),0));
+ TVector3 vE = TVector3(m_XE[0],m_YE[0],m_ZE[0]);
+ TVector3 vA = TVector3(x_start_beam,y_start_beam,z_start_beam);
+
+ double ThetaIn = InitialConditions->GetIncidentEmittanceTheta();
+ double Phi = InitialConditions->GetIncidentEmittancePhi();
+
+ TVector3 dir = TVector3(sin(ThetaIn*deg)*cos(Phi*deg), sin(ThetaIn*deg)*sin(Phi*deg), cos(ThetaIn*deg));
+ m_TOF.push_back(Time);
+ m_Adet.push_back(Adet);
+ m_Brho.push_back(Brho);
+ m_ThetaIn.push_back(ThetaIn);
+
+ double Brho_calc = m_GladField->FindBrho(vA,dir,vE);
+ double PathLength = m_GladField->GetFlightPath(vA,Brho_calc,vA,dir)/1000;
+ double velocity = PathLength/Time;
+ 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*Zdet;
+
+ m_Brho_calc.push_back(Brho_calc);
+ m_FlightPath.push_back(PathLength);
+ m_A_calc.push_back(A_calc);
+ }
+
+}
+
+
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
@@ -120,21 +192,27 @@ void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchStatus("fDetecetd_*",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&InteractionCoordinates);
- RootInput::getInstance()->GetChain()->SetBranchStatus("FissionConditions",true);
- RootInput::getInstance()->GetChain()->SetBranchStatus("fFC_*",true);
- RootInput::getInstance()->GetChain()->SetBranchAddress("FissionConditions",&FissionConditions);
-
+ /*if(beam_ana==false){
+ RootInput::getInstance()->GetChain()->SetBranchStatus("FissionConditions",true);
+ RootInput::getInstance()->GetChain()->SetBranchStatus("fFC_*",true);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("FissionConditions",&FissionConditions);
+ }*/
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("TOF",&m_TOF);
+ RootOutput::getInstance()->GetTree()->Branch("ThetaIn",&m_ThetaIn);
RootOutput::getInstance()->GetTree()->Branch("Brho",&m_Brho);
- RootOutput::getInstance()->GetTree()->Branch("A",&m_A);
+ RootOutput::getInstance()->GetTree()->Branch("Asim",&m_Asim);
+ 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);
RootOutput::getInstance()->GetTree()->Branch("FlightPath",&m_FlightPath);
+ RootOutput::getInstance()->GetTree()->Branch("XE",&m_XE);
+ RootOutput::getInstance()->GetTree()->Branch("YE",&m_YE);
+ RootOutput::getInstance()->GetTree()->Branch("ZE",&m_ZE);
}
@@ -142,7 +220,12 @@ void Analysis::InitOutputBranch(){
void Analysis::ReInit(){
m_TOF.clear();
m_Brho.clear();
- m_A.clear();
+ m_Asim.clear();
+ m_Adet.clear();
+ m_ThetaIn.clear();
+ m_XE.clear();
+ m_YE.clear();
+ m_ZE.clear();
m_Brho_calc.clear();
m_A_calc.clear();
diff --git a/Projects/Sofia/Analysis.h b/Projects/Sofia/Analysis.h
index 636524f5d..50fc980f5 100644
--- a/Projects/Sofia/Analysis.h
+++ b/Projects/Sofia/Analysis.h
@@ -40,21 +40,28 @@ class Analysis: public NPL::VAnalysis{
void TreatEvent();
void ReInit();
void End();
+ void FFAnalysis();
+ void BeamAnalysis();
static NPL::VAnalysis* Construct();
private:
vector<double> m_TOF;
+ vector<double> m_ThetaIn;
vector<double> m_Brho;
- vector<int> m_A;
+ vector<int> m_Asim;
+ vector<int> m_Adet;
vector<double> m_Brho_calc;
vector<int> m_A_calc;
vector<double> m_FlightPath;
-
+ vector<double> m_XE;
+ vector<double> m_YE;
+ vector<double>