diff --git a/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.cxx b/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.cxx
index c8030d807ff8107d408a56fb52c3b39f22185247..4d0982a3a705a6fc82c329f8123b88fa67aa2c95 100644
--- a/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.cxx
+++ b/NPLib/Detectors/ComptonTelescope/TComptonTelescopePhysics.cxx
@@ -89,7 +89,7 @@ void TComptonTelescopePhysics::BuildSimplePhysicalEvent()
//// DSSSD analysis ////
// Check event type
- Int_t evtType = CheckEvent();
+ //Int_t evtType = CheckEvent();
//cout << "event type = " << evtType << endl;
// check possible interstrip
@@ -192,7 +192,7 @@ void TComptonTelescopePhysics::BuildSimplePhysicalEvent()
//cout << "couple " << i << " CT" << Tower << " D" << N << " SF" << Front << " SB" << Back << " EF " << Front_E << " EB " << Back_E << " TF " << Front_T << " TB " << Back_T << endl;
// Fill TComptonTelescopePhysics members
- EventType.push_back(evtType);
+ //EventType.push_back(evtType);
TowerNumber.push_back(Tower);
DetectorNumber.push_back(N);
Strip_Front.push_back(Front);
diff --git a/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx b/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
index 9bfc53f313adb71d0996453765fcdf35b017331e..42ad1e3903d2e389c28407b30d4d211703acbf57 100644
--- a/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
+++ b/NPLib/Detectors/Samurai/SamuraiFieldMap.cxx
@@ -47,16 +47,19 @@ SamuraiFieldMap::SamuraiFieldMap(){
}
////////////////////////////////////////////////////////////////////////////////
double SamuraiFieldMap::Delta(const double* parameter){
- static vector<TVector3>pos ;
- static TVector3 diff;
- pos =Propagate(parameter[0],m_FitPosFDC0,m_FitDirFDC0,false);
- // Move the fdc2 pos from lab frame to fdc2 frame
- pos.back().RotateY(-m_fdc2angle+m_angle);
-
- //double d = (pos.back().X()-m_FitPosFDC2.X())*(pos.back().X()-m_FitPosFDC2.X());
- // return d;
- diff = pos.back()-m_FitPosFDC2;
- return diff.Mag2();
+ static vector<TVector3>pos ;
+ static TVector3 diff;
+ static int i = 0;
+ //pos =Propagate(parameter[0],m_FitPosFDC0,m_FitDirFDC0,false);
+ pos =Propagate(parameter[0],m_FitPosFDC0,m_FitDirFDC0,true);
+ // Move the fdc2 pos from lab frame to fdc2 frame
+ pos.back().RotateY(-m_fdc2angle+m_angle);
+
+
+ //double d = (pos.back().X()-m_FitPosFDC2.X())*(pos.back().X()-m_FitPosFDC2.X());
+ // return d;
+ diff = pos.back()-m_FitPosFDC2;
+ return diff.Mag2();
}
////////////////////////////////////////////////////////////////////////////////
@@ -69,7 +72,9 @@ double SamuraiFieldMap::FindBrho(TVector3 p_fdc0,TVector3 d_fdc0,TVector3 p_fdc2
if(!m_BrhoScan)
BrhoScan(1,10,0.1);
// do a first guess based on fdc2 pos
- double b0[1] ={m_BrhoScan->Eval(p_fdc2.X())};
+ double b0[1] ={m_BrhoScan->Eval(p_fdc2.X())};
+ //cout << "First guess Brho " << b0[0] << " "; //endl;
+
m_min->Clear();
m_min->SetPrecision(1e-6);
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index f7c7551413e47cd51ee2955f29fe38101eca41a9..f7277bdd7df1fd7588b7913a89114206cfe01dec 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -233,19 +233,11 @@ void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, doubl
fXThetaXHist->GetRandom2(X,ThetaX);
fYPhiYHist->GetRandom2(Y,PhiY);
}
- // Direction
- double Xdir = sin(ThetaX); // cos(90-x) = sin(x)
- double Ydir = sin(PhiY);
- double Zdir = sqrt(1-Xdir*Xdir-Ydir*Ydir); // alpha^2 + beta^2 + gamma^2 = 1
- // Stretch factor to extend unitary vector from ZEmission to ZProfile
- double S = fZProfile-fZEmission ;
- TVector3 BeamDir(Xdir*S/Zdir,Ydir*S/Zdir,S);
-
- Xdir = BeamDir.X();
- Ydir = BeamDir.Y();
- // POSITION/DIRECTION AT Z EMISSION//
- X = X0 - Xdir;
- Y = Y0 - Ydir;
+ // Backward propagtion from ZProfile (in general on target)
+ // to ZEmission (position where beam is generated in simulation.)
+ double dZ = fZProfile-fZEmission ; //(Z1-Z0)
+ X= X0-tan(ThetaX)*dZ;
+ Y= Y0-tan(PhiY)*dZ;
Z = fZEmission;
}
diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index b19851ca6e57f04a99f0d9bd7eae04f764df1235..587919787ea6cbb1859a616b0989b37aabe9ad31 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -78,7 +78,11 @@ QFS::QFS(){
fshootB=false;
fshoot1=true;
fshoot2=true;
- fisotropic = true;
+ fIsotropic = true;
+
+ fCondEcmPos = false;
+ fCondOffshellMassPos = false;
+ fIsAllowed = false;
fUseExInGeant4=true;
@@ -220,6 +224,9 @@ Particle QFS::GetParticle(string name, NPL::InputParser parser){
////////////////////////////////////////////////////////////////////////////////////////////
void QFS::CalculateVariables(){
+ fCondEcmPos = false;
+ fCondOffshellMassPos = false;
+
if(fBeamEnergy < 0)
fBeamEnergy = 0 ;
@@ -250,7 +257,8 @@ void QFS::CalculateVariables(){
// 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;}
+ if(buffer>0) {fCondOffshellMassPos = true;}
+ else{/*cout<<"ERROR off shell mass ma_off=\t"<<buffer<<endl; Dump();*/ return;}
ma_off = sqrt(buffer);
//deduced total energies of "a" and "B" in restframe of A
@@ -279,7 +287,8 @@ void QFS::CalculateVariables(){
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;Dump(); return;}
+ if(fEcm>0) {fCondEcmPos = true;}
+ if(fEcm<=0) { /*cout<<"ERROR Ecm negative =\t"<<fEcm<<endl;Dump();*/ 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));
@@ -298,8 +307,11 @@ void QFS::CalculateVariables(){
void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEnergyLab1, double& ThetaLab2, double& PhiLab2, double& KineticEnergyLab2){
+ fIsAllowed=false;
CalculateVariables();
-
+ if(fCondOffshellMassPos && fCondEcmPos) fIsAllowed= true;
+ else return;
+
double thetaCM_1 = fThetaCM;
double thetaCM_2 = M_PI - thetaCM_1;
double phiCM_2 = fPhiCM;
@@ -501,17 +513,6 @@ TGraph* QFS::GetPhi2VsPhi1(double AngleStep_CM){
return(fPhi2VsPhi1);
}
-///////////////////////////////////////////////////////////////////////////////
-// Check whenever the reaction is allowed at a given energy
-bool QFS::IsAllowed(){//double Energy){
- //double AvailableEnergy = Energy + fParticle1.Mass() + fParticle2.Mass();
- //double RequiredEnergy = fParticle3.Mass() + fParticle4.Mass();
-
- //if(AvailableEnergy>RequiredEnergy)
- return true;
- //else
- // return false;
-}
////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index 0afde4ccda55ae51395e2116f19b584d1c115b59..c5954ef76ee1cd8eb3dfde21efaf66651ef3c4a1 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -83,7 +83,11 @@ namespace NPL{
TVector3 fInternalMomentum; // Internal momentum of the removed cluster
TH1F* fPerpMomentumHist; // Perpendicular momentum distribution in beam at rest frame
TH1F* fParMomentumHist; // Parallel momentum distribution in beam at rest frame
- double fisotropic;
+ double fIsotropic;
+
+ bool fCondEcmPos;
+ bool fCondOffshellMassPos;
+ bool fIsAllowed;
TGraph* fTheta2VsTheta1;
TGraph* fPhi2VsPhi1;
@@ -103,7 +107,6 @@ namespace NPL{
void TestR3B();
void KineRelativistic(double &ThetaLab1, double &PhiLab1, double &KineticEnergyLab1, double &ThetaLab2, double &PhiLab2, double &KineticEnergyLab2);
TVector3 ShootInternalMomentum();
- bool IsAllowed();
void Dump();
private: // intern precompute variable
@@ -169,6 +172,7 @@ namespace NPL{
void SetPhiCM(const double& angle) {fPhiCM = angle;}
void SetInternalMomentum(const TVector3& mom) {fInternalMomentum = mom;}
void SetMomentumSigma(const double& sigma) {fMomentumSigma = sigma;}
+ void SetIsAllowed(const bool& val) {fIsAllowed = val;}
void SetPerpMomentumHist (TH1F* PerpMomentumHist)
{delete fPerpMomentumHist; fPerpMomentumHist = PerpMomentumHist;}
void SetParMomentumHist (TH1F* ParMomentumHist)
@@ -187,7 +191,8 @@ namespace NPL{
double GetThetaCM() const {return fThetaCM;}
double GetPhiCM() const {return fPhiCM;}
double GetMomentumSigma() const {return fMomentumSigma;}
- bool GetUseExInGeant4() const { return fUseExInGeant4; }
+ bool IsAllowed() const {return fIsAllowed;}
+ bool GetUseExInGeant4() const { return fUseExInGeant4; }
double GetExcitationA() const {return fExcitationA;}
double GetExcitationB() const {return fExcitationB;}
TVector3 GetInternalMomentum() const {return fInternalMomentum;}
diff --git a/NPLib/Physics/TReactionConditions.cxx b/NPLib/Physics/TReactionConditions.cxx
index 460075917dcdfd0c6dab8a20379546905751a64c..f4882241d3a74dc507eeebbe8ba27d3594ff1905 100644
--- a/NPLib/Physics/TReactionConditions.cxx
+++ b/NPLib/Physics/TReactionConditions.cxx
@@ -72,7 +72,8 @@ void TReactionConditions::Dump() const{
cout << "\t Phi_Y: " << fRC_Beam_Emittance_PhiY << endl;
cout << "\t Theta: " << fRC_Beam_Emittance_Theta << endl;
cout << "\t Phi: " << fRC_Beam_Emittance_Phi << endl;
-
+ cout << "\t Direction: " ;
+ GetBeamDirection().Print();
// Beam status at the initial interaction point
cout << "\t ---- Interaction Point ---- " << endl;
@@ -87,17 +88,23 @@ void TReactionConditions::Dump() const{
<< fRC_Internal_Momentum.Y() << " ; "
<< fRC_Internal_Momentum.Z() << ")" << endl;
-
+ TVector3 *emitted= new TVector3();
// emmitted particle
unsigned int size = fRC_Particle_Name.size();
for(unsigned int i = 0 ; i < size; i ++){
cout << "\t ---- Particle " << i << " ---- " << endl;
cout << "\t Particle Name: " << fRC_Particle_Name[i] << endl;
- cout << "\t Kinetic Energy: " << fRC_Kinetic_Energy[i] << endl;
- cout << "\t Momentum Direction: ( "
+ cout << "\t Kinetic Energy: " << fRC_Kinetic_Energy[i] << endl;
+ cout << "\t Angles in beam frame (along Z)"<< endl;
+ cout << "\t Theta: " << fRC_Theta[i] << endl;
+ cout << "\t Phi: " << fRC_Phi[i] << endl;
+ cout << "\t Momentum Direction in world frame: ( "
<< fRC_Momentum_Direction_X[i] << " ; "
<< fRC_Momentum_Direction_Y[i] << " ; "
<< fRC_Momentum_Direction_Z[i] << ")" << endl;
+ emitted->SetXYZ(fRC_Momentum_Direction_X[i],fRC_Momentum_Direction_Y[i],fRC_Momentum_Direction_Z[i]);
+ cout << "\t ThetaWorld: " << emitted->Theta()*180./pi<< endl;
+ cout << "\t PhiWorld: " << emitted->Phi()*180./pi<< endl;
}
diff --git a/NPLib/Physics/TReactionConditions.h b/NPLib/Physics/TReactionConditions.h
index 1a85cda8c9fb448ba6c1f33e5490b542ee6f91c2..ade0078b866ab31e50e5d22d1626030c07a4bba8 100644
--- a/NPLib/Physics/TReactionConditions.h
+++ b/NPLib/Physics/TReactionConditions.h
@@ -48,27 +48,33 @@ private:
// Beam beam parameter
string fRC_Beam_Particle_Name;
- double fRC_Beam_Emittance_ThetaX;
- double fRC_Beam_Emittance_PhiY;
- double fRC_Beam_Emittance_Theta;
- double fRC_Beam_Emittance_Phi;
- double fRC_Beam_Reaction_Energy;
-
+ double fRC_Beam_Emittance_ThetaX; //beam angle between Pxz and Z axis
+ double fRC_Beam_Emittance_PhiY; //beam angle between Pyz and Z axis
+ double fRC_Beam_Emittance_Theta; //spher. theta (betw. beam dir. and Z axis)
+ double fRC_Beam_Emittance_Phi; //spher. phi (betw. Pyz and X axis)
+ double fRC_Beam_Reaction_Energy; //beam kinetic energy at vertex
+
+ //Reaction vertex coordinates
double fRC_Vertex_Position_X;
double fRC_Vertex_Position_Y;
double fRC_Vertex_Position_Z;
+ //Center of mass angle for the reaction
+ double fRC_ThetaCM;
+
+ //Two-Body reaction: Exc. energy of the two products
double fRC_ExcitationEnergy3;
double fRC_ExcitationEnergy4;
- double fRC_ThetaCM;
+ //QFS reacion only: Internal mom. of the removed particle
TVector3 fRC_Internal_Momentum;
- // emmitted particles
+
+ // Emitted reaction products properties
vector<string> fRC_Particle_Name;
- vector<double> fRC_Theta;
- vector<double> fRC_Phi;
+ vector<double> fRC_Theta; //in the frame with beam on Z axis
+ vector<double> fRC_Phi; //in the frame with beam on Z axis
vector<double> fRC_Kinetic_Energy;
- vector<double> fRC_Momentum_Direction_X;
- vector<double> fRC_Momentum_Direction_Y;
- vector<double> fRC_Momentum_Direction_Z;
+ vector<double> fRC_Momentum_Direction_X; //in the world frame
+ vector<double> fRC_Momentum_Direction_Y; //in the world frame
+ vector<double> fRC_Momentum_Direction_Z; //in the world frame
public:
TReactionConditions();
@@ -139,31 +145,29 @@ 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 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 ;
-
- double GetThetaLab_WorldFrame(const int i) const{
+ double GetTheta_WorldFrame(const int i) const{
return (GetParticleDirection(i).Theta())/deg;
}
- double GetThetaLab_BeamFrame (const int i) const{
- return (GetParticleDirection(i).Angle(GetBeamDirection()))/deg;
- }
-
- double GetPhiLab_WorldFrame (const int i) const {
+ double GetPhi_WorldFrame (const int i) const {
return (M_PI + GetParticleDirection(i).Phi())/deg;
// to have Phi in [0,2pi]] and not [-pi,pi]]
}
- double GetPhiLab_BeamFrame (const int i) const{
+
+ //Two following methods should return angles identical as GetTheta() and GetPhi()
+ //spherical angles when beam axis along Z
+ // Only used for consistency checks
+ double GetTheta_BeamFrame (const int i) const{
+ return (GetParticleDirection(i).Angle(GetBeamDirection()))/deg;
+ }
+ double GetPhi_BeamFrame (const int i) const{
TVector3 rot = GetParticleDirection(i);
rot.RotateUz(GetBeamDirection());
return (M_PI + rot.Phi())/deg;
}
-
unsigned int GetEmittedMult() const {return fRC_Particle_Name.size();}
ClassDef(TReactionConditions, 1) // TReactionConditions structure
diff --git a/NPSimulation/Detectors/Samurai/CMakeLists.txt b/NPSimulation/Detectors/Samurai/CMakeLists.txt
index 96e0eaedf049e6bfd4b49a3c52b0e9812b1634cd..c77d1a47210d739baf592595c3fd6c6187abcf33 100644
--- a/NPSimulation/Detectors/Samurai/CMakeLists.txt
+++ b/NPSimulation/Detectors/Samurai/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSSamurai SHARED Samurai.cc SamuraiFieldPropagation.cc SamuraiFDC2.cc)
+add_library(NPSSamurai SHARED Samurai.cc SamuraiFieldPropagation.cc SamuraiFDC2.cc SamuraiFDC1.cc)
target_link_libraries(NPSSamurai NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPSamurai)
diff --git a/NPSimulation/Detectors/Samurai/Samurai.cc b/NPSimulation/Detectors/Samurai/Samurai.cc
index 2a2731ff586bae503b12e519cf5db4c34bb2ee57..92da9b578022c6d2181756c6c7612959ce89451f 100644
--- a/NPSimulation/Detectors/Samurai/Samurai.cc
+++ b/NPSimulation/Detectors/Samurai/Samurai.cc
@@ -47,7 +47,7 @@
// NPTool header
#include "Samurai.hh"
#include "SamuraiFieldPropagation.hh"
-#include "CalorimeterScorers.hh"
+//#include "CalorimeterScorers.hh"
#include "InteractionScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
@@ -69,7 +69,8 @@ namespace Samurai_NS{
const double Magnet_Width = 6700*mm; //(x)
const double Magnet_Height = 4640*mm;//(y)
const double Magnet_Depth = 6000*mm;//(z)
- const string Magnet_Material = "G4_Galactic";
+ const double Magnet_Depth_eps = 0.2*mm; //(delta-z)
+ const string Magnet_Material = "G4_Galactic"; //G4_Galactic
//Gap between top and bottom yoke slabs
const double Magnet_Gap = 880*mm;//(y)
@@ -99,6 +100,10 @@ namespace Samurai_NS{
const double pv_main_Height = Magnet_Gap;
const double pv_main_Depth = Magnet_Depth;
const string Propvol_Material = "G4_Galactic";
+
+ //Ideal Data saving location
+ const short int Magnet_DetectorNumber = -1;
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -122,6 +127,12 @@ Samurai::Samurai(){
//Propagation region
m_PropagationRegion = NULL;
+
+ //Scorer
+ m_SamuraiScorer = NULL;
+
+ //Data
+ m_Event = new TSamuraiIdealData;
}
Samurai::~Samurai(){
@@ -191,7 +202,7 @@ G4LogicalVolume* Samurai::BuildMagnet(){
if(!m_Magnet){
//Shape - G4Box
G4Box* box = new G4Box("Samurai_Box",Samurai_NS::Magnet_Width*0.5,
- Samurai_NS::Magnet_Height*0.5,Samurai_NS::Magnet_Depth*0.5);
+ Samurai_NS::Magnet_Height*0.5,Samurai_NS::Magnet_Depth*0.5 + Samurai_NS::Magnet_Depth_eps*0.5);
//Material - vacuum
G4Material* VacuumMaterial = MaterialManager::getInstance()
@@ -200,6 +211,8 @@ G4LogicalVolume* Samurai::BuildMagnet(){
//Logical Volume
m_Magnet = new G4LogicalVolume(box, VacuumMaterial, "logic_Samurai_box",0,0,0);
m_Magnet->SetVisAttributes(m_VisMagnet);
+ m_Magnet->SetSensitiveDetector(m_SamuraiScorer);
+ cout << "ATTACHING SCORER\n";
}
return m_Magnet;
}
@@ -420,70 +433,63 @@ void Samurai::SetPropagationRegion(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
-//FIXME
-/*
void Samurai::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("Samurai")){
- pTree->Branch("Samurai", "TSamuraiData", &m_Event) ;
+ if(!pTree->FindBranch("IdealDataMagnet")){
+ pTree->Branch("IdealDataMagnet", "TSamuraiIdealData", &m_Event) ;
}
- pTree->SetBranchAddress("Samurai", &m_Event) ;
+ pTree->SetBranchAddress("IdealDataMagnet", &m_Event) ;
}
-*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
-//FIXME
-void Samurai::ReadSensitive(const G4Event* ){
-}
-
-/*
-//FIXME
void Samurai::ReadSensitive(const G4Event* ){
- m_Event->Clear();
- ///////////
- // Calorimeter scorer
- CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_SamuraiScorer->GetPrimitive(0);
+ //cout << "Read Sensitive\n";
+ m_Event->Clear();
+ //Interaction Scorer
+ InteractionScorers::PS_Interactions* Scorer= (InteractionScorers::PS_Interactions*) m_SamuraiScorer->GetPrimitive(0);
unsigned int size = Scorer->GetMult();
+ //cout << "size " << size << endl;
for(unsigned int i = 0 ; i < size ; i++){
- vector<unsigned int> level = Scorer->GetLevel(i);
- double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Samurai_NS::ResoEnergy);
- if(Energy>Samurai_NS::EnergyThreshold){
- double Time = RandGauss::shoot(Scorer->GetTime(i),Samurai_NS::ResoTime);
- int DetectorNbr = level[0];
- m_Event->SetEnergy(DetectorNbr,Energy);
- m_Event->SetTime(DetectorNbr,Time);
- }
+ double energy = Scorer->GetEnergy(i);
+ double brho = Scorer->GetBrho(i);
+ double posx = Scorer->GetPositionX(i);
+ double posy = Scorer->GetPositionY(i);
+ double posz = Scorer->GetPositionZ(i);
+ double mom_mag = brho*Scorer->GetCharge(i);
+ double theta = Scorer->GetTheta(i);
+ double phi = Scorer->GetPhi(i);
+ m_Event->SetData(Samurai_NS::Magnet_DetectorNumber, energy, posx, posy, posz, mom_mag, theta, phi, brho);
+ //cout << Samurai_NS::Magnet_DetectorNumber << endl;
}
+
}
-*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////
-/*
-//FIXME
-void Samurai::InitializeScorers() {
+
+void Samurai::InitializeScorers() {
+ cout << "I'M IN InitializeScorers" << endl;
// This check is necessary in case the geometry is reloaded
bool already_exist = false;
m_SamuraiScorer = CheckScorer("SamuraiScorer",already_exist) ;
+ if (already_exist) cout << "//////////////// already exist ///////// " << endl;
+ else cout << "//////////////// does not already exist ///////// " << endl;
if(already_exist)
return ;
// Otherwise the scorer is initialised
vector<int> level; level.push_back(0);
- G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",level, 0) ;
G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ;
//and register it to the multifunctionnal detector
- m_SamuraiScorer->RegisterPrimitive(Calorimeter);
m_SamuraiScorer->RegisterPrimitive(Interaction);
G4SDManager::GetSDMpointer()->AddNewDetector(m_SamuraiScorer) ;
}
-*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Detectors/Samurai/Samurai.hh b/NPSimulation/Detectors/Samurai/Samurai.hh
index dabe269d1e0d99a7f38561063025445b60f3dfa2..83fc8ec51a07f2254a8203f5088f6f37f7809054 100644
--- a/NPSimulation/Detectors/Samurai/Samurai.hh
+++ b/NPSimulation/Detectors/Samurai/Samurai.hh
@@ -39,6 +39,7 @@ using namespace std;
#include "NPInputParser.h"
#include "SamuraiFieldPropagation.hh"
+#include "TSamuraiIdealData.h"
class Samurai : public NPS::VDetector{
@@ -92,13 +93,13 @@ class Samurai : public NPS::VDetector{
// Called in DetecorConstruction::ReadDetectorConfiguration Method
void ReadConfiguration(NPL::InputParser) ;
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
+ // Construct Magnet and inialise sensitive part.
+ // Called After DetectorConstruction::AddDetector Method
void ConstructDetector(G4LogicalVolume* world) ;
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- //void InitializeRootOutput() ;
+ // Add Magnet branch to the EventTree.
+ // Called After DetectorConstruction::AddDetector Method
+ void InitializeRootOutput() ;
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAction
@@ -106,18 +107,20 @@ class Samurai : public NPS::VDetector{
public:
// Scorer
- // Initialize all Scorer used by the MUST2Array
- //void InitializeScorers() ;
+ // Initialize all Scorer used by the MUST2Array
+ void InitializeScorers() ;
// Set region were magnetic field is active:
void SetPropagationRegion();
- // Associated Scorer
- //G4MultiFunctionalDetector* m_Sa//(x)muraiScorer ;
+ // Associated Scorer
+ G4MultiFunctionalDetector* m_SamuraiScorer ;
////////////////////////////////////////////////////
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
private:
+ // Store data
+ TSamuraiIdealData* m_Event;
// Region were magnetic field is active:
G4Region* m_PropagationRegion;
diff --git a/NPSimulation/Detectors/Samurai/SamuraiFDC1.cc b/NPSimulation/Detectors/Samurai/SamuraiFDC1.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3fc0ffbe90e36d3d52d7f2b1a6adb72c18e98ad7
--- /dev/null
+++ b/NPSimulation/Detectors/Samurai/SamuraiFDC1.cc
@@ -0,0 +1,266 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2021 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: Omar Nasr contact address: omar.h.nasr@outlook.com *
+ * *
+ * Creation Date : septembre 2021 *
+ * Last update : septembre 2021 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe Samurai simulation *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// C++ headers
+#include <sstream>
+#include <cmath>
+//G4 Geometry object
+#include "G4Box.hh"
+
+//G4 sensitive
+#include "G4SDManager.hh"
+#include "G4MultiFunctionalDetector.hh"
+
+//G4 various object
+#include "G4Material.hh"
+#include "G4Transform3D.hh"
+#include "G4PVPlacement.hh"
+#include "G4VisAttributes.hh"
+#include "G4Colour.hh"
+#include "G4RegionStore.hh"
+
+// NPTool header
+#include "SamuraiFDC1.hh"
+//#include "CalorimeterScorers.hh"
+#include "InteractionScorers.hh"
+#include "RootOutput.h"
+#include "MaterialManager.hh"
+#include "NPSDetectorFactory.hh"
+#include "NPOptionManager.h"
+#include "NPSHitsMap.hh"
+// CLHEP header
+#include "CLHEP/Random/RandGauss.h"
+
+using namespace std;
+using namespace CLHEP;
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+namespace SamuraiFDC1_NS{
+ // Samurai magnet construction paramethers
+
+ //Main outer box
+ const double FDC1_Width = 1000*mm; //(x)
+ const double FDC1_Height = 696*mm;//(y)
+ const double FDC1_Depth = 336*mm;//(z)
+ const string FDC1_Material = "G4_Galactic";
+
+ //Detector Number
+ const short int FDC1_DetectorNumber = 1;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Samurai Specific Method
+SamuraiFDC1::SamuraiFDC1(){
+
+ //Visualization attributes
+ m_VisFDC1 = new G4VisAttributes(G4Colour(1,1,0,0.5));
+ //Logical volumes
+ m_FDC1 = NULL;
+ //Scorer
+ m_FDC1Scorer = NULL;
+ //Ideal Data event
+ m_Event = new TSamuraiIdealData;
+
+}
+
+SamuraiFDC1::~SamuraiFDC1(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void SamuraiFDC1::AddDetector(G4ThreeVector Mag_Pos, G4ThreeVector Offset){
+
+ m_Pos = Mag_Pos + Offset;
+
+ return;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4LogicalVolume* SamuraiFDC1::BuildFDC1(){
+ if(!m_FDC1){
+ //Shape - G4Box
+ G4Box* box = new G4Box("FDC1_Box",SamuraiFDC1_NS::FDC1_Width*0.5,
+ SamuraiFDC1_NS::FDC1_Height*0.5,SamuraiFDC1_NS::FDC1_Depth*0.5);
+
+ //Material - vacuum
+ G4Material* VacuumMaterial = MaterialManager::getInstance()
+ ->GetMaterialFromLibrary(SamuraiFDC1_NS::FDC1_Material);
+
+ //Logical Volume
+ m_FDC1 = new G4LogicalVolume(box, VacuumMaterial, "logic_SamuraiFDC1_box",0,0,0);
+ m_FDC1->SetVisAttributes(m_VisFDC1);
+ m_FDC1->SetSensitiveDetector(m_FDC1Scorer);
+ }
+ return m_FDC1;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Virtual Method of NPS::VDetector class
+
+// Read stream at Configfile to pick-up parameters of detector (Position,...)
+// (Called in DetecorConstruction::ReadDetectorConfiguration Method)
+void SamuraiFDC1::ReadConfiguration(NPL::InputParser parser){
+
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Samurai");
+ vector<NPL::InputBlock*> blocks2 = parser.GetAllBlocksWithToken("SAMURAIFDC1");
+
+ if(blocks.size()==1 && blocks2.size()==1){
+ if(NPOptionManager::getInstance()->GetVerboseLevel()) {
+ cout << "/////// Samurai FDC1 found ///////" << endl;
+ }
+ vector<string> cart = {"POS","ANGLE"};
+ vector<string> sphe = {"R","Theta","Phi","ANGLE"};
+
+ G4ThreeVector Mag_Pos;
+
+ if(blocks[0]->HasTokenList(cart)){
+ Mag_Pos = NPS::ConvertVector(blocks[0]->GetTVector3("POS", "cm"));
+ }
+ else if(blocks[0]->HasTokenList(sphe)){
+ double R = blocks[0]->GetDouble("R","mm");
+ double Theta = blocks[0]->GetDouble("Theta","deg");
+ double Phi = blocks[0]->GetDouble("Phi","deg");
+ Mag_Pos.setMag(R);
+ Mag_Pos.setTheta(Theta);
+ Mag_Pos.setPhi(Phi);
+ }
+
+ G4ThreeVector Offset = NPS::ConvertVector(blocks2[0]->GetTVector3("Offset", "mm"));
+
+ //string xml = blocks2[0]->GetString("XML");
+ //bool invert_x = blocks2[0]->GetBool("InvertX");
+ //bool invert_y = blocks2[0]->GetBool("InvertY");
+ //bool invert_z = blocks2[0]->GetBool("InvertD");
+
+ AddDetector(Mag_Pos, Offset);
+
+ }
+ else{
+ cout << "ERROR: there should be only one Samurai magnet, check your input file" << endl;
+ exit(1);
+ }
+
+
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Construct detector and inialise sensitive part.
+// (Called After DetectorConstruction::AddDetector Method)
+void SamuraiFDC1::ConstructDetector(G4LogicalVolume* world){
+
+ new G4PVPlacement(0, m_Pos,
+ BuildFDC1(), "SamuraiFDC1", world, false, 0);
+
+ return;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Add Detector branch to the EventTree.
+// Called After DetecorConstruction::AddDetector Method
+void SamuraiFDC1::InitializeRootOutput(){
+ RootOutput *pAnalysis = RootOutput::getInstance();
+ TTree *pTree = pAnalysis->GetTree();
+ if(!pTree->FindBranch("IdealDataFDC1")){
+ pTree->Branch("IdealDataFDC1", "TSamuraiIdealData", &m_Event) ;
+ }
+ pTree->SetBranchAddress("IdealDataFDC1", &m_Event);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+// Read sensitive part and fill the Root tree.
+// (Called at in the EventAction::EndOfEventAvtion)
+
+void SamuraiFDC1::ReadSensitive(const G4Event* event){
+
+ m_Event->Clear();
+ //Interaction Scorer
+ InteractionScorers::PS_Interactions* Scorer= (InteractionScorers::PS_Interactions*) m_FDC1Scorer->GetPrimitive(0);
+
+ unsigned int size = Scorer->GetMult();
+ for(unsigned int i = 0 ; i < size ; i++){
+ //vector<unsigned int> level = Scorer->GetLevel(i);
+ //double Energy = RandGauss::shoot(Scorer->GetEnergy(i),SamuraiFDC1_NS::ResoEnergy);
+ //double Energy = Scorer->GetEnergy(i);
+
+ double energy = Scorer->GetEnergy(i);
+ double brho = Scorer->GetBrho(i);
+ double posx = Scorer->GetPositionX(i);
+ double posy = Scorer->GetPositionY(i);
+ double posz = Scorer->GetPositionZ(i);
+ double mom_mag = brho*Scorer->GetCharge(i);
+ double theta = Scorer->GetTheta(i);
+ double phi = Scorer->GetPhi(i);
+ m_Event->SetData(SamuraiFDC1_NS::FDC1_DetectorNumber, energy, posx, posy, posz, mom_mag, theta, phi, brho);
+ }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void SamuraiFDC1::InitializeScorers() {
+ // This check is necessary in case the geometry is reloaded
+ bool already_exist = false;
+ m_FDC1Scorer = CheckScorer("FDC1Scorer",already_exist) ;
+
+ if(already_exist)
+ return ;
+
+ // Otherwise the scorer is initialised
+ vector<int> level; level.push_back(0);
+ //Interaction
+ G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ;
+ m_FDC1Scorer->RegisterPrimitive(Interaction);
+ G4SDManager::GetSDMpointer()->AddNewDetector(m_FDC1Scorer) ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+
+NPS::VDetector* SamuraiFDC1::Construct(){
+ return (NPS::VDetector*) new SamuraiFDC1();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern"C" {
+ class proxy_nps_samuraiFDC1{
+ public:
+ proxy_nps_samuraiFDC1(){
+ NPS::DetectorFactory::getInstance()->AddToken("SAMURAIFDC1","SAMURAIFDC1");
+ NPS::DetectorFactory::getInstance()->AddDetector("SAMURAIFDC1",SamuraiFDC1::Construct);
+ }
+ };
+
+ proxy_nps_samuraiFDC1 p_nps_samuraiFDC1;
+}
+
+
diff --git a/NPSimulation/Detectors/Samurai/SamuraiFDC1.hh b/NPSimulation/Detectors/Samurai/SamuraiFDC1.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e18be19416240c61c3e789be4a0ed7f59284ac6d
--- /dev/null
+++ b/NPSimulation/Detectors/Samurai/SamuraiFDC1.hh
@@ -0,0 +1,127 @@
+#ifndef Samurai_h
+#define Samurai_h 1
+/*****************************************************************************
+ * Copyright (C) 2009-2021 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: Omar Nasr contact address: omar.h.nasr@outlook.com *
+ * *
+ * Creation Date : septembre 2021 *
+ * Last update : septembre 2021 *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class describe Samurai simulation *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+// C++ header
+#include <string>
+#include <vector>
+using namespace std;
+
+// G4 headers
+#include "G4ThreeVector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4LogicalVolume.hh"
+#include "G4MultiFunctionalDetector.hh"
+//#include "G4VFastSimulationModel.hh"
+#include "G4VSolid.hh"
+
+// NPTool header
+#include "NPSVDetector.hh"
+#include "NPInputParser.h"
+
+//#include "SamuraiFieldPropagation.hh"
+#include "TSamuraiIdealData.h"
+
+
+class SamuraiFDC1 : public NPS::VDetector{
+ ////////////////////////////////////////////////////
+ /////// Default Constructor and Destructor /////////
+ ////////////////////////////////////////////////////
+ public:
+ SamuraiFDC1() ;
+ virtual ~SamuraiFDC1() ;
+
+ ////////////////////////////////////////////////////
+ /////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+ public:
+
+ // Cartezian FDC1
+ void AddDetector(G4ThreeVector Mag_Pos, G4ThreeVector Offset);
+
+ G4LogicalVolume* BuildFDC1();
+ private:
+
+ //Logical Volume
+ G4LogicalVolume* m_FDC1;
+
+ // Visualisation Attributes
+ G4VisAttributes* m_VisFDC1;
+
+
+ ////////////////////////////////////////////////////
+ ////// Inherite from NPS::VDetector class /////////
+ ////////////////////////////////////////////////////
+ public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetectorConfiguration Method
+ void ReadConfiguration(NPL::InputParser) ;
+
+ // Construct detector and initialise sensitive part.
+ // (Called After DetecorConstruction::AddDetector Method)
+ void ConstructDetector(G4LogicalVolume* world) ;
+
+ // Add Detector branch to the EventTree.
+ // (Called After DetecorConstruction::AddDetector Method)
+ void InitializeRootOutput() ;
+
+ // Read sensitive part and fill the Root tree.
+ // (Called at in the EventAction::EndOfEventAction)
+ void ReadSensitive(const G4Event* event) ;
+
+ public:
+ // Scorer
+ // Initialize the scorer(s) used by the FDC1 detector
+ void InitializeScorers() ;
+
+
+ // Associated Scorer
+ G4MultiFunctionalDetector* m_FDC1Scorer ;
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+ private:
+ //TSamuraiFDC1Data* m_Event;
+ //////////////////////////////////////////////////////////////////
+
+ TSamuraiIdealData* m_Event;
+
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+ private:
+ //Detector coordinates
+ G4ThreeVector m_Pos;
+
+
+ // Needed for dynamic loading of the library
+ public:
+ static NPS::VDetector* Construct();
+};
+#endif
+
+
+
+
+
+
+
diff --git a/NPSimulation/Detectors/Samurai/SamuraiFDC2.cc b/NPSimulation/Detectors/Samurai/SamuraiFDC2.cc
index d898c2ab4de76d78c7d449c3d15dc2ca77ea2cd8..e403e631735fbbcd135259176de575dc2a8fa28d 100644
--- a/NPSimulation/Detectors/Samurai/SamuraiFDC2.cc
+++ b/NPSimulation/Detectors/Samurai/SamuraiFDC2.cc
@@ -22,7 +22,6 @@
// C++ headers
#include <sstream>
#include <cmath>
-#include <limits>
//G4 Geometry object
#include "G4Box.hh"
@@ -37,12 +36,10 @@
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4RegionStore.hh"
-#include "G4FastSimulationManager.hh"
-#include "G4UserLimits.hh"
// NPTool header
#include "SamuraiFDC2.hh"
-#include "CalorimeterScorers.hh"
+//#include "CalorimeterScorers.hh"
#include "InteractionScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
@@ -61,16 +58,18 @@ namespace SamuraiFDC2_NS{
// Samurai magnet construction paramethers
//Main outer box
- const double FDC2_Width = 6700*mm; //(x)
- const double FDC2_Height = 4640*mm;//(y)
- const double FDC2_Depth = 500*mm;//(z)
- const string FDC2_Material = "G4_Fe";
+ const double FDC2_Width = 2616*mm; //(x)
+ const double FDC2_Height = 1156*mm;//(y)
+ const double FDC2_Depth = 876*mm;//(z)
+ const string FDC2_Material = "G4_Galactic";
+
+ //Detector Number
+ const short int FDC2_DetectorNumber = 2;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
// Samurai Specific Method
SamuraiFDC2::SamuraiFDC2(){
@@ -80,7 +79,7 @@ SamuraiFDC2::SamuraiFDC2(){
m_FDC2 = NULL;
//Scorer
m_FDC2Scorer = NULL;
-
+ //Ideal Data event
m_Event = new TSamuraiIdealData;
}
@@ -98,40 +97,9 @@ void SamuraiFDC2::AddDetector(G4ThreeVector Mag_Pos, double Mag_Angle, G4ThreeVe
return;
}
-/*
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void SamuraiFDC2::AddFDC2(double r, double theta, double phi,
- double R, double Theta, double Phi, double Angle){
-
- m_R_d = r;
- m_Theta_d = theta;
- m_Phi_d = phi;
-
- m_R = R;
- m_Theta = Theta;
- m_Phi = Phi;
- m_Angle = Angle;
-
- return;
-}
-//Spherical Magnet
-void SamuraiFDC2::AddMagnet(double R, double Theta, double Phi, double Angle){
-
-
- return;
-}
-//Cartezian Magnet
-void SamuraiFDC2::AddMagnet(G4ThreeVector POS, double Angle){
-
- m_R = POS.mag();
- m_Theta = POS.theta();
- m_Phi = POS.phi();
- m_Angle = Angle;
-
- return;
-}*/
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
G4LogicalVolume* SamuraiFDC2::BuildFDC2(){
if(!m_FDC2){
//Shape - G4Box
@@ -150,15 +118,12 @@ G4LogicalVolume* SamuraiFDC2::BuildFDC2(){
return m_FDC2;
}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
// Virtual Method of NPS::VDetector class
// Read stream at Configfile to pick-up parameters of detector (Position,...)
-// Called in DetecorConstruction::ReadDetectorConfiguration Method
+// (Called in DetecorConstruction::ReadDetectorConfiguration Method)
void SamuraiFDC2::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Samurai");
@@ -211,7 +176,7 @@ void SamuraiFDC2::ReadConfiguration(NPL::InputParser parser){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Construct detector and inialise sensitive part.
-// Called After DetectorConstruction::AddDetector Method
+// (Called After DetectorConstruction::AddDetector Method)
void SamuraiFDC2::ConstructDetector(G4LogicalVolume* world){
G4RotationMatrix* Rot = new G4RotationMatrix();
@@ -224,26 +189,22 @@ void SamuraiFDC2::ConstructDetector(G4LogicalVolume* world){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
-//FIXME
-
void SamuraiFDC2::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("IdealData")){
- pTree->Branch("IdealData", "TSamuraiIdealData", &m_Event) ; //WATCH OUT !!!!!!
+ if(!pTree->FindBranch("IdealDataFDC2")){
+ pTree->Branch("IdealDataFDC2", "TSamuraiIdealData", &m_Event) ;
}
- pTree->SetBranchAddress("IdealData", &m_Event) ; //WATCH OUT !!!!!!
+ pTree->SetBranchAddress("IdealDataFDC2", &m_Event);
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
// Read sensitive part and fill the Root tree.
-// Called at in the EventAction::EndOfEventAvtion
+// (Called at in the EventAction::EndOfEventAvtion)
-//FIXME
void SamuraiFDC2::ReadSensitive(const G4Event* event){
m_Event->Clear();
@@ -255,7 +216,7 @@ void SamuraiFDC2::ReadSensitive(const G4Event* event){
//vector<unsigned int> level = Scorer->GetLevel(i);
//double Energy = RandGauss::shoot(Scorer->GetEnergy(i),SamuraiFDC2_NS::ResoEnergy);
//double Energy = Scorer->GetEnergy(i);
- short int detector = 2;
+
double energy = Scorer->GetEnergy(i);
double brho = Scorer->GetBrho(i);
double posx = Scorer->GetPositionX(i);
@@ -264,38 +225,12 @@ void SamuraiFDC2::ReadSensitive(const G4Event* event){
double mom_mag = brho*Scorer->GetCharge(i);
double theta = Scorer->GetTheta(i);
double phi = Scorer->GetPhi(i);
- m_Event->SetData(detector, energy, posx, posy, posz, mom_mag, theta, phi, brho);
+ m_Event->SetData(SamuraiFDC2_NS::FDC2_DetectorNumber, energy, posx, posy, posz, mom_mag, theta, phi, brho);
}
- /*
- m_Event->Clear();
-
- ///////////
- // Calorimeter scorer
- //CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_FDC2Scorer->GetPrimitive(0);
-
- //Interaction Scorer
- InteractionScorers::PS_Interactions* Scorer= (InteractionScorers::PS_Interactions*) m_FDC2Scorer->GetPrimitive(0);
-
- unsigned int size = Scorer->GetMult();
- for(unsigned int i = 0 ; i < size ; i++){
- vector<unsigned int> level = Scorer->GetLevel(i);
- //double Energy = RandGauss::shoot(Scorer->GetEnergy(i),SamuraiFDC2_NS::ResoEnergy);
- double Energy = Scorer->GetEnergy();
-
- if(Energy>SamuraiFDC2_NS::EnergyThreshold){
- double Time = RandGauss::shoot(Scorer->GetTime(i),SamuraiFDC2_NS::ResoTime);
- int DetectorNbr = level[0];
- m_Event->SetEnergy(DetectorNbr,Energy);
- m_Event->SetTime(DetectorNbr,Time);
- }
- }*/
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-////////////////////////////////////////////////////////////////
-
-//FIXME
void SamuraiFDC2::InitializeScorers() {
// This check is necessary in case the geometry is reloaded
bool already_exist = false;
@@ -306,21 +241,17 @@ void SamuraiFDC2::InitializeScorers() {
// Otherwise the scorer is initialised
vector<int> level; level.push_back(0);
- //G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter",level, 0) ;
- G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ; //WATCH OUT
- //and register it to the multifunctionnal detector
- //m_FDC2Scorer->RegisterPrimitive(Calorimeter);
+ //Interaction
+ G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ;
m_FDC2Scorer->RegisterPrimitive(Interaction);
G4SDManager::GetSDMpointer()->AddNewDetector(m_FDC2Scorer) ;
}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
+
NPS::VDetector* SamuraiFDC2::Construct(){
return (NPS::VDetector*) new SamuraiFDC2();
}
@@ -329,7 +260,7 @@ NPS::VDetector* SamuraiFDC2::Construct(){
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
-extern"C" {////// WHAT IS THIS??
+extern"C" {
class proxy_nps_samuraiFDC2{
public:
proxy_nps_samuraiFDC2(){
diff --git a/NPSimulation/Detectors/Samurai/SamuraiFDC2.hh b/NPSimulation/Detectors/Samurai/SamuraiFDC2.hh
index 79afe5eede1a14ed349b3b16ffdcf85006706fed..1dc277fc1347124b30879263f3a70a0e16cee23f 100644
--- a/NPSimulation/Detectors/Samurai/SamuraiFDC2.hh
+++ b/NPSimulation/Detectors/Samurai/SamuraiFDC2.hh
@@ -31,14 +31,14 @@ using namespace std;
#include "G4RotationMatrix.hh"
#include "G4LogicalVolume.hh"
#include "G4MultiFunctionalDetector.hh"
-#include "G4VFastSimulationModel.hh"
+//#include "G4VFastSimulationModel.hh"
#include "G4VSolid.hh"
// NPTool header
#include "NPSVDetector.hh"
#include "NPInputParser.h"
-#include "SamuraiFieldPropagation.hh"
+//#include "SamuraiFieldPropagation.hh"
#include "TSamuraiIdealData.h"
@@ -55,18 +55,13 @@ class SamuraiFDC2 : public NPS::VDetector{
////////////////////////////////////////////////////
public:
- // Cartesian fdc0
+ // Cartezian FDC2
void AddDetector(G4ThreeVector Mag_Pos, double Mag_Angle, G4ThreeVector Offset, double Off_Angle);
- // Spherical fdc0
- /*void AddFDC2(double R, double Theta, double Phi, double Angle,
- double r, double theta, double phi);*/
- //G4VSolid* BuildRYokeSolid();
G4LogicalVolume* BuildFDC2();
private:
- //Solids used for the magnet construction
- //G4VSolid* m_RYokeSolid;
+ //Logical Volume
G4LogicalVolume* m_FDC2;
// Visualisation Attributes
@@ -81,21 +76,21 @@ class SamuraiFDC2 : public NPS::VDetector{
// Called in DetecorConstruction::ReadDetectorConfiguration Method
void ReadConfiguration(NPL::InputParser) ;
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
+ // Construct detector and initialise sensitive part.
+ // (Called After DetecorConstruction::AddDetector Method)
void ConstructDetector(G4LogicalVolume* world) ;
// Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
+ // (Called After DetecorConstruction::AddDetector Method)
void InitializeRootOutput() ;
// Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAction
+ // (Called at in the EventAction::EndOfEventAction)
void ReadSensitive(const G4Event* event) ;
public:
// Scorer
- // Initialize all Scorer used by the MUST2Array
+ // Initialize the scorer(s) used by the FDC2 detector
void InitializeScorers() ;
diff --git a/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.cc b/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.cc
index a21467a3095760d2c97e6f1ae800fe6b8cdc12ba..6571c9964ad050c1bd9581f853b43f20b0abd8ef 100644
--- a/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.cc
+++ b/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.cc
@@ -49,11 +49,11 @@
#include "RootOutput.h"
#include "TLorentzVector.h"
-bool single_particle = false;
-ofstream outEO ("EliaOmar.txt" );
-ofstream outRK ("RungeKutta.txt" );
-ofstream outEOsp ("EliaOmarSP.txt" , ios::app );
-ofstream outRKsp ("RungeKuttaSP.txt" , ios::app );
+//bool single_particle = false;
+//ofstream outEO ("EliaOmar.txt" );
+//ofstream outRK ("RungeKutta.txt" );
+//ofstream outEOsp ("EliaOmarSP.txt" , ios::app );
+//ofstream outRKsp ("RungeKuttaSP.txt" , ios::app );
////////////////////////////////////////////////////////////////////////////////
NPS::SamuraiFieldPropagation::SamuraiFieldPropagation(G4String modelName, G4Region* envelope)
: G4VFastSimulationModel(modelName, envelope) {
@@ -251,8 +251,8 @@ void NPS::SamuraiFieldPropagation::EliaOmarPropagation (const G4FastTrack& fastT
newMomentum = localMomentum;
//benchmark
- if(single_particle) PrintData(m_StepSize, newPosition, newMomentum, outEOsp);
- else PrintData(m_StepSize, newPosition, newMomentum, outEO);
+ //if(single_particle) PrintData(m_StepSize, newPosition, newMomentum, outEOsp);
+ //else PrintData(m_StepSize, newPosition, newMomentum, outEO);
}
double time = PrimaryTrack->GetGlobalTime()+(newPosition - localPosition).mag()/speed;
@@ -287,7 +287,7 @@ void NPS::SamuraiFieldPropagation::RungeKuttaPropagation (const G4FastTrack& fas
double charge = PrimaryTrack->GetParticleDefinition()->GetPDGCharge() / coulomb;
double ConF_p = 1 / (joule * c_light / (m/s)) ; // MeV/c to kg*m/s (SI units)
- //Initially inside is false
+ //Initially inside is false -> calculate trajectory
if (!inside){
count = 2;//skip first two positions as they are the same as the current position
double Brho = localMomentum.mag() * ConF_p / charge;
@@ -304,8 +304,8 @@ void NPS::SamuraiFieldPropagation::RungeKuttaPropagation (const G4FastTrack& fas
G4ThreeVector newPosition (trajectory[count].x(), trajectory[count].y(), trajectory[count].z());
//benchmark
- G4ThreeVector newDir = (newPosition - localPosition).unit();
- G4ThreeVector newMomentum = newDir * localMomentum.mag();
+ //G4ThreeVector newDir = (newPosition - localPosition).unit();
+ //G4ThreeVector newMomentum = newDir * localMomentum.mag();
//Check if newPosition is not inside
if (solid->Inside(newPosition) != kInside){
@@ -313,18 +313,18 @@ void NPS::SamuraiFieldPropagation::RungeKuttaPropagation (const G4FastTrack& fas
G4ThreeVector toOut = solid->DistanceToOut(localPosition, localDir) * localDir;
newPosition = localPosition + toOut;
//benchmark
- newDir = (newPosition - localPosition).unit();
- newMomentum = newDir * localMomentum.mag();
+ //newDir = (newPosition - localPosition).unit();
+ //newMomentum = newDir * localMomentum.mag();
//benchmark
- if(single_particle) PrintData(m_StepSize, newPosition, newMomentum, outRKsp);
- else PrintData(m_StepSize, newPosition, newMomentum, outRK);
+ //if(single_particle) PrintData(m_StepSize, newPosition, newMomentum, outRKsp);
+ //else PrintData(m_StepSize, newPosition, newMomentum, outRK);
//counter++;
}
//benchmark
- //G4ThreeVector newDir = (newPosition - localPosition).unit();
- //G4ThreeVector newMomentum = newDir * localMomentum.mag();
+ G4ThreeVector newDir = (newPosition - localPosition).unit();
+ G4ThreeVector newMomentum = newDir * localMomentum.mag();
double time = PrimaryTrack->GetGlobalTime()+(newPosition - localPosition).mag()/speed;
diff --git a/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.hh b/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.hh
index 8914d824d0fdd1664b6fd5e0c6893e5929353a37..8cfc6c13a0c8bf7b50a1ab4e2860cf43f51e0d52 100644
--- a/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.hh
+++ b/NPSimulation/Detectors/Samurai/SamuraiFieldPropagation.hh
@@ -77,12 +77,12 @@ namespace NPS{
return to_string(V.getX()) + " " + to_string(V.getY()) + " " + to_string(V.getZ());
}
- void PrintData(double stepsize, G4ThreeVector pos, G4ThreeVector mom, ofstream& out){
- out << setprecision(17) << stepsize/micrometer << "\t";
- out << setprecision(17) << pos.x() << "\t" << pos.y() << "\t" << pos.z() << "\t";
- out << setprecision(17) << mom.x() << "\t" << mom.y() << "\t" << mom.z() << "\t";
- out << setprecision(17) << mom.getR() << "\t" << mom.getPhi() << "\t" << mom.getTheta() << endl;
- }
+ //void PrintData(double stepsize, G4ThreeVector pos, G4ThreeVector mom, ofstream& out){
+ // out << setprecision(17) << stepsize/micrometer << "\t";
+ // out << setprecision(17) << pos.x() << "\t" << pos.y() << "\t" << pos.z() << "\t";
+ // out << setprecision(17) << mom.x() << "\t" << mom.y() << "\t" << mom.z() << "\t";
+ // out << setprecision(17) << mom.getR() << "\t" << mom.getPhi() << "\t" << mom.getTheta() << endl;
+ //}
bool m_Initialized; //field map initialized
diff --git a/NPSimulation/EventGenerator/EventGeneratorBeam.cc b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
index d169d083067a2a4f8f3c1217c4e8a81957db65c4..4e33ce3cf296bcfc2ae0d443f2409589355d3ae9 100644
--- a/NPSimulation/EventGenerator/EventGeneratorBeam.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
@@ -97,12 +97,13 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
double InitialBeamEnergy, x0, y0, z0, Beam_thetaX, Beam_phiY;
m_Beam->GenerateRandomEvent(InitialBeamEnergy, x0, y0, z0, Beam_thetaX, Beam_phiY);
- //Set the direction cosines: alpha=90-Beam_thetaX, beta=90-Beam_phiY
- double Xdir = sin(Beam_thetaX); // cos(90-x) = sin(x)
- double Ydir = sin(Beam_phiY);
- double Zdir = sqrt(1-Xdir*Xdir-Ydir*Ydir); // alpha^2 + beta^2 + gamma^2 = 1
+ double Xdir = tan(Beam_thetaX); //tan(thetax)= px/pz
+ double Ydir = tan(Beam_phiY); //tan(phiy)= py/pz
+ double Zdir = 1; // fix pz=1 arbitrarily
G4ThreeVector BeamDir(Xdir,Ydir,Zdir);
+ BeamDir = BeamDir.unit();
G4ThreeVector BeamPos(x0,y0,z0);
+
///////////////////////////////////////////////////////
///// Add the Beam particle to the particle Stack /////
///////////////////////////////////////////////////////
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 189d397dc3e91dbeb312c14bdf742b88424c6673..1a8ba7f3edb150c20c14cf7562837c173ced425d 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -183,13 +183,13 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
// If the condition is met, the event is generated
if (m_shoot && m_length < m_StepSize) {
if(m_ReactionType==QFS){
- if ( m_QFS.IsAllowed() ) {
+ //if ( m_QFS.IsAllowed() ) {
return true;
- }
- else{
- m_shoot=false;
- std::cout << "QFS not allowed" << std::endl;
- }
+ //}
+ //else{
+ // m_shoot=false;
+ // std::cout << "QFS not allowed" << std::endl;
+ //}
}
else if(m_ReactionType==TwoBody){
@@ -475,14 +475,6 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
// Set the Energy of the reaction
m_QFS.SetBeamEnergy(reac_energy);
- double Beam_theta = pdirection.theta();
- double Beam_phi = pdirection.phi();
-
-
- /////////////////////////////////////////////////////////////////
- ///// Angles for emitted particles following Cross Section //////
- ///// Angles are in the beam frame //////
- /////////////////////////////////////////////////////////////////
// Shoot and Set a Random ThetaCM
double costheta = G4RandFlat::shoot() * 2 - 1;
@@ -491,45 +483,55 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_QFS.SetThetaCM(theta);
m_QFS.SetPhiCM(phi);
- // Shoot and set internal momentum for the removed cluster
- // if a momentum Sigma is given then shoot in 3 indep. Gaussians
- // if input files are given for distributions use them instead
- m_QFS.ShootInternalMomentum();
-
- //////////////////////////////////////////////////
- ///// Momentum and angles from kinematics /////
- //////////////////////////////////////////////////
- // Variable where to store results
+ // Lab frame variables where to store results
double Theta1, Phi1, TKE1, Theta2, Phi2, TKE2, ThetaB, PhiB, TKEB;
- // Compute Kinematic using previously defined ThetaCM
- m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
+ int j=0;
+ m_QFS.SetIsAllowed(false);
+ while(!m_QFS.IsAllowed()){
+ // Shoot internal momentum for the removed cluster
+ // if a momentum Sigma is given then shoot in 3 indep. Gaussians
+ // if input files are given for distributions use them instead
+ m_QFS.ShootInternalMomentum();
- G4ThreeVector U(1, 0, 0);
- G4ThreeVector V(0, 1, 0);
- G4ThreeVector ZZ(0, 0, 1);
+ // Go from CM to Lab
+ m_QFS.KineRelativistic(Theta1, Phi1, TKE1, Theta2, Phi2, TKE2);
+
+ j++;
+ if(j>100) cout<< "ERROR: too many iteration and QFS kinematical conditions not allowed"<<endl;
+ }
+
+ //---------------------------------------------------------
+ //Rotations to switch from frame with beam along Z to world
+ //---------------------------------------------------------
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+ G4ThreeVector ux(1, 0, 0);
+ G4ThreeVector uy(0, 1, 0);
+ G4ThreeVector uz(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_kine1_world.rotate(Beam_theta, uy); // rotation of Beam_theta around Y axis
+ momentum_kine1_world.rotate(Beam_phi, uz); // rotation of Beam_phi around 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_kine2_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ momentum_kine2_world.rotate(Beam_phi, uz); // 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
+ //momentum_kineB_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ //momentum_kineB_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
TLorentzVector* P_A = m_QFS.GetEnergyImpulsionLab_A();
TLorentzVector* P_B = m_QFS.GetEnergyImpulsionLab_B();
@@ -538,12 +540,13 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
momentum_kineB_beam = momentum_kineB_beam.unit();
TKEB = P_B->Energy() - m_QFS.GetParticleB()->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
+ momentum_kineB_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ momentum_kineB_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
ThetaB = P_B->Angle(P_A->Vect());
- if (ThetaB < 0) ThetaB += M_PI;
- PhiB = M_PI + P_B->Vect().Phi();
+ //if (ThetaB < 0) ThetaB += M_PI;
+ //PhiB = M_PI + P_B->Vect().Phi();
+ PhiB = P_B->Vect().Phi();
if (fabs(PhiB) < 1e-6) PhiB = 0;
@@ -602,44 +605,45 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
m_ReactionConditions->SetVertexPositionX(localPosition.x());
m_ReactionConditions->SetVertexPositionY(localPosition.y());
m_ReactionConditions->SetVertexPositionZ(localPosition.z());
+
m_ReactionConditions->SetBeamEmittanceTheta(
PrimaryTrack->GetMomentumDirection().theta() / deg);
m_ReactionConditions->SetBeamEmittancePhi(
PrimaryTrack->GetMomentumDirection().phi() / deg);
m_ReactionConditions->SetBeamEmittanceThetaX(
- PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+ PrimaryTrack->GetMomentumDirection().perpPart(uy).angle(uz) / deg);
m_ReactionConditions->SetBeamEmittancePhiY(
- PrimaryTrack->GetMomentumDirection().angle(V) / deg);
+ PrimaryTrack->GetMomentumDirection().perpPart(ux).angle(uz) / deg);
// Names 1,2 and B//
m_ReactionConditions->SetParticleName(Light1Name->GetParticleName());
m_ReactionConditions->SetParticleName(Light2Name->GetParticleName());
m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
- // Angle 1,2 and B //
+ // Angle 1,2 and B in fram with beam along Z axis//
m_ReactionConditions->SetTheta(Theta1 / deg);
m_ReactionConditions->SetTheta(Theta2 / deg);
m_ReactionConditions->SetTheta(ThetaB / deg);
m_ReactionConditions->SetPhi(Phi1 / deg);
m_ReactionConditions->SetPhi(Phi2 / deg);
m_ReactionConditions->SetPhi(PhiB / deg);
- // Energy 1,2 and B //
+ // Total Kinetic Energy 1,2 and B //
m_ReactionConditions->SetKineticEnergy(TKE1);
m_ReactionConditions->SetKineticEnergy(TKE2);
m_ReactionConditions->SetKineticEnergy(TKEB);
// ThetaCM, Ex and Internal Momentum of removed cluster//
m_ReactionConditions->SetThetaCM(m_QFS.GetThetaCM() / deg);
m_ReactionConditions->SetInternalMomentum(m_QFS.GetInternalMomentum());
- //m_ReactionConditions->SetExcitationEnergy3(m_QFS.GetExcitation3());
+ //Excitation energy of reaction product B
m_ReactionConditions->SetExcitationEnergy4(m_QFS.GetExcitationB());
- // Momuntum X 1,2 and B //
+ // Momentum Dir. X 1,2 and B in world frame //
m_ReactionConditions->SetMomentumDirectionX(momentum_kine1_world.x());
m_ReactionConditions->SetMomentumDirectionX(momentum_kine2_world.x());
m_ReactionConditions->SetMomentumDirectionX(momentum_kineB_world.x());
- // Momuntum Y 1,2 and B //
+ // Momentum Dir. Y 1,2 and B in world frame//
m_ReactionConditions->SetMomentumDirectionY(momentum_kine1_world.y());
m_ReactionConditions->SetMomentumDirectionY(momentum_kine2_world.y());
m_ReactionConditions->SetMomentumDirectionY(momentum_kineB_world.y());
- // Momuntum Z 1,2 and B //
+ // Momentum Dir. Z 1,2 and B in world frame//
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine1_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kine2_world.z());
m_ReactionConditions->SetMomentumDirectionZ(momentum_kineB_world.z());
diff --git a/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C b/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
index 8f06d8d0e0563e9205a691939be705739935b411..5d40ac68940a6866ff8ea84c7e63fdc5858f546a 100644
--- a/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
+++ b/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
@@ -162,6 +162,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
: Form("h_D%d_BACK_E%d",DETECTOR_ID,n+1);
auto h = (TH1F*) f->Get(hname.c_str());
h->Sumw2();
+ //h->GetXaxis()->SetLimits(10,1024);
+ h->SetAxisRange(1,1024,"X");
// draw histogram
can->cd(1);
auto pad = (TPad*) can->FindObject("can_1");
diff --git a/Projects/ComptonTelescope/online/Calibration_DSSSD/ExtractRawHisto_DSSSD_E.C b/Projects/ComptonTelescope/online/Calibration_DSSSD/ExtractRawHisto_DSSSD_E.C
index f4b1f1a5308bac09a63502e5214fc96694163130..3f1ee3a6cd732ba63cc142adb9ce6a3d04bb716b 100644
--- a/Projects/ComptonTelescope/online/Calibration_DSSSD/ExtractRawHisto_DSSSD_E.C
+++ b/Projects/ComptonTelescope/online/Calibration_DSSSD/ExtractRawHisto_DSSSD_E.C
@@ -91,14 +91,19 @@ void ExtractRawHisto_DSSSD_E(const char* filename = "20200128_10h44_bi207_conv")
switch (i) {
case 0: //Assuming 0 is front - to be checked
for (int k = 0; k < NBSTRIPS; k++) { // strips
- hFrontEnergy[j][k]->Fill(event->sample[i][j][k]);
- hFront[j]->Fill(k, event->sample[i][j][k]);
+ //if (event->sample[i][j][k] != 0) {
+ hFrontEnergy[j][k]->Fill(event->sample[i][j][k]);
+ hFront[j]->Fill(k, event->sample[i][j][k]);
+ //cout << "E = " << event->sample[i][j][k] << endl;
+ //}
}
break;
case 1://Assuming 1 is back - to be checked
- for (int k = 0; k < NBSTRIPS; k++) { // strips
- hBackEnergy[j][k]->Fill(event->sample[i][j][k]);
- hBack[j]->Fill(k, event->sample[i][j][k]);
+ for (int k = 0; k < NBSTRIPS; k++) { // strips
+ //if (event->sample[i][j][k] != 0) {
+ hBackEnergy[j][k]->Fill(event->sample[i][j][k]);
+ hBack[j]->Fill(k, event->sample[i][j][k]);
+ //}
}
}
}
diff --git a/Projects/ComptonTelescope/online/calibrations.txt b/Projects/ComptonTelescope/online/calibrations.txt
index 5f15faf94b9435da5397855df503bd79eed893c1..271959706716dc3fef83baba282d8696ee58047d 100644
--- a/Projects/ComptonTelescope/online/calibrations.txt
+++ b/Projects/ComptonTelescope/online/calibrations.txt
@@ -1,7 +1,8 @@
CalibrationFilePath
% ./calibrations/DSSSD_D1_Calibration_run10_newCal.txt
% ./calibrations/DSSSD_calibration_withPed_pol3_ok.txt
- ./calibrations/DSSSD_D1_Calibration_run13.txt
+% ./calibrations/DSSSD_D1_Calibration_run13.txt
+ ./calibrations/DSSSD_D1_Calibration_210923.txt
./calibrations/CeBr3_PED.txt
./calibrations/CeBr3_PED3.txt
./calibrations/CeBr3_PED4.txt
diff --git a/Projects/ComptonTelescope/online/configs/ConfigComptonTelescope.dat b/Projects/ComptonTelescope/online/configs/ConfigComptonTelescope.dat
index 8fb526fa62a3ae3dfd5a4dd0ba166048c78777ef..5cb483e94975e30fc13e849a8e10d3254d3920ef 100644
--- a/Projects/ComptonTelescope/online/configs/ConfigComptonTelescope.dat
+++ b/Projects/ComptonTelescope/online/configs/ConfigComptonTelescope.dat
@@ -3,15 +3,11 @@ ConfigComptonTelescope
ONLY_GOOD_COUPLE_WITH_MULTIPLICITY_ONE ON
FRONT_BACK_ENERGY_MATCHING_SIGMA 6
FRONT_BACK_ENERGY_MATCHING_NUMBER_OF_SIGMA 2
- DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_FRONT25_E
+ DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_FRONT1_E
DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_FRONT26_E
- DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_FRONT30_E
DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK0_E
- DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK1_E
DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK4_E
- DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK12_E
DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK13_E
- DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK14_E
DISABLE_CHANNEL COMPTONTELESCOPE_D1_STRIP_BACK26_E
STRIP_FRONT_E_RAW_THRESHOLD 0
STRIP_BACK_E_RAW_THRESHOLD 0
diff --git a/Projects/ComptonTelescope/online/src/online_coinc.cpp b/Projects/ComptonTelescope/online/src/online_coinc.cpp
index 8a08af820097286adc477b7c4c276d66a49f2a4e..0658d5c58c4576c0b8bc1cc213c172905a6ae29e 100644
--- a/Projects/ComptonTelescope/online/src/online_coinc.cpp
+++ b/Projects/ComptonTelescope/online/src/online_coinc.cpp
@@ -16,16 +16,16 @@
#include "DecodeD.h"
#include "DecodeT.h"
-#define __RUN__ 7
+#define __RUN__ 3
#define __OUTPUT_ALL__
#undef __OUTPUT_ALL__
#define __1DET__
-#undef __1DET__
+//#undef __1DET__
#define __TEST_ZONE__
-//#undef __TEST_ZONE__
+#undef __TEST_ZONE__
#define __CIRCULAR_TREE__
#undef __CIRCULAR_TREE__
@@ -120,7 +120,9 @@ int main(int argc, char** argv)
#if __RUN__ == 3 || __RUN__ == 7
ifstream is;
#if __RUN__ == 3
- is.open("/disk/proto-data/data/20210510_Bi207/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-05-10_07_41_50.raw");
+ //is.open("/disk/proto-data/data/20210510_Bi207/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-05-10_07_41_50.raw");
+ //is.open("/projet/astronuc/data/ccam/data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001");
+ is.open("/projet/astronuc/data/ccam/data/banc_104/211109/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-11-09_16_06_21.raw.0001");
#elif __RUN__ == 7
#endif
is.seekg(0, ios::end);
@@ -146,7 +148,8 @@ int main(int argc, char** argv)
m_NPDetectorManager->ClearEventPhysics();
m_NPDetectorManager->ClearEventData();
#if __RUN__ == 3 || __RUN__ == 7
- ccamData -> SetCTCalorimeter(1, 4, DR->getPixelNumber(), DR->getTime(), DR->getData(), 64);
+ ccamData -> SetCTCalorimeter(1, 4, Decoder->getPixelNumber(), Decoder->getTime(), Decoder->getData(), 64);
+ //ccamData -> SetCTCalorimeter(1, 4, DR->getPixelNumber(), DR->getTime(), DR->getData(), 64);
#elif __RUN__ == 0
setCTTracker(ccamData, Decoder);
//ccamData -> Dump();
diff --git a/Projects/ComptonTelescope/online/src/online_dsssd.cpp b/Projects/ComptonTelescope/online/src/online_dsssd.cpp
index 228f3e43e5044da909f1d33ce0e890100c4dc316..7bb12f7d2d1747e10765b7d363439b1d1a27796f 100644
--- a/Projects/ComptonTelescope/online/src/online_dsssd.cpp
+++ b/Projects/ComptonTelescope/online/src/online_dsssd.cpp
@@ -30,7 +30,8 @@ int main(int argc, char *argv[])
NPOptionManager::getInstance(arg);
// open ROOT output file
- RootOutput::getInstance("OnlineTree_DSSSD.root", "OnlineTree");
+ RootOutput::getInstance("OnlineTree.root", "OnlineTree");
+ //RootOutput::getInstance("OnlineTree_DSSSD.root", "OnlineTree");
// get tree pointer
auto m_OutputTree = RootOutput::getInstance()->GetTree();
diff --git a/Projects/Samurai/Samurai.detector b/Projects/Samurai/Samurai.detector
index 5ce30b423c92a9667738f7eca3d44143e16b4a93..53e981171413777e66ca006fa8570d48c60a3baf 100644
--- a/Projects/Samurai/Samurai.detector
+++ b/Projects/Samurai/Samurai.detector
@@ -9,7 +9,7 @@ Target
Z= -10000 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Samurai
- POS= 0 0 0 mm
+ POS= 0 0 10000 mm
ANGLE= 30 deg
%
METHOD= 1
@@ -19,9 +19,10 @@ Samurai
% fieldmap path
STEPS_PER_METER= 1000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%SAMURAIFDC0
+SAMURAIFDC1
%XML= db/SAMURAIFDC0_20200109.xml
%Offset= -0.00666226 0.102191 -3370.01 mm
+ Offset= -0.00666226 0.102191 -4370.01 mm
%InvertX= 1
%InvertY= 0
%InvertD= 1
@@ -29,7 +30,7 @@ Samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMURAIFDC2
%XML= db/SAMURAIFDC2.xml
- Offset= 1000 1000 5000 mm
+ Offset= 1000 200 5000 mm
OffAngle= 59.930 deg
%-252.416 -0.228477 4122.57 mm
%InvertX= 0
diff --git a/Projects/Samurai/benchmark/bm2/bm2ReconstructBrho.C b/Projects/Samurai/benchmark/bm2/bm2ReconstructBrho.C
new file mode 100644
index 0000000000000000000000000000000000000000..a5f55efde71179dd35041e88559a2fb45d20f3ba
--- /dev/null
+++ b/Projects/Samurai/benchmark/bm2/bm2ReconstructBrho.C
@@ -0,0 +1,181 @@
+
+double mean(vector<double>V){
+ double s = 0;
+ for (auto x:V) s += x;
+ return s / V.size();
+}
+double var(vector<double>V){
+ double s = 0, m=mean(V);
+ for (auto x:V) s += (x-m)*(x-m);
+ return s / (V.size()-1);
+}
+
+
+void Brhos(string path_to_file);
+void findbrho(string path_to_file, int time_int, SamuraiFieldMap* Map, double mag_angle, double fdc2_angle,
+ const vector<double>& Brho, const vector<TVector3>& Pos_FDC1,
+ const vector<TVector3>& Pos_FDC2, const vector<TVector3>& Dir_FDC1,
+ const vector<TVector3>& Dir_FDC2);
+
+
+void bm2ReconstructBrho (){
+
+ vector<string> paths ={"spm1000/1cm_E0/", "spm1000/1cm_E480/", "spm1000/5cm_E480/",
+ "spm10000/1cm_E0/", "spm10000/1cm_E480/", "spm10000/5cm_E480/"};
+ for(auto p:paths) Brhos(p);
+
+ /*
+ for (auto i =0; i<Br_rec.size(); i++) if(Br_rec[i] < 6) cout << i << " " << Br_rec[i] << endl;
+ TH1D* hist1d = new TH1D("brho", "Brho", 200, 6, 7);
+ for (int i = 0; i < Br_rec.size(); i++) hist1d->Fill(Br_rec[i]);
+ TCanvas* tc = new TCanvas();
+ tc->SetLogy();
+ hist1d->Draw();*/
+/*
+ vector<TVector3> track1 = Map->Propagate(6.26709, Pos_FDC1[7839], Dir_FDC1[7839], true);
+ vector<TVector3> track2 = Map->Propagate(Br_rec[839], Pos_FDC1[7839], Dir_FDC1[7839], true);
+ vector<TVector3> track3 = Map->Propagate(1, Pos_FDC1[7839], Dir_FDC1[7839], true);
+
+ vector<double> x1, z1, x2, z2, x3, z3;
+ for (int i=0; i<track1.size(); i++){
+ x1.push_back(-track1[i].X());
+ z1.push_back(track1[i].Z());
+ }
+ for (int i=0; i<track2.size(); i++){
+ x2.push_back(-track2[i].X());
+ z2.push_back(track2[i].Z());
+ }
+ for (int i=0; i<track3.size(); i++){
+ x3.push_back(-track3[i].X());
+ z3.push_back(track3[i].Z());
+ }
+
+ TGraph* tg1 = new TGraph(x1.size(), &x1[0], &z1[0]);
+ TGraph* tg2 = new TGraph(x2.size(), &x2[0], &z2[0]);
+ TGraph* tg3 = new TGraph(x3.size(), &x3[0], &z3[0]);
+ TMarker* mk = new TMarker(-Pos_FDC2[7839].X(), Pos_FDC2[7839].Z(), 29);
+ tg1->SetMarkerStyle(2);
+ tg2->SetMarkerStyle(3);
+ tg3->SetMarkerStyle(4);
+
+ tg1->Draw("APL");
+ tg2->Draw("PLsame");
+ tg3->Draw("PLsame");
+ mk->Draw();
+ */
+
+ // cout << "Real position FDC1: " << "\t" << Pos_FDC1[0].X() << " " << Pos_FDC1[0].Y() << " " << Pos_FDC1[0].Z() << endl;
+ // cout << "Real direction FDC1: " << "\t" << Dir_FDC1[0].X() << " " << Dir_FDC1[0].Y() << " " << Dir_FDC1[0].Z() << endl;
+ // cout << "Real position FDC2: " << "\t" << Pos_FDC2[0].X() << " " << Pos_FDC2[0].Y() << " " << Pos_FDC2[0].Z() << endl;
+ // cout << "Real direction FDC2: " << "\t" << Dir_FDC2[0].X() << " " << Dir_FDC2[0].Y() << " " << Dir_FDC2[0].Z() << endl;
+ // cout << "I found Brho = " << br << endl;
+ // cout << "Real Brho = " << Brho[0] << endl;
+
+
+
+ //TH1D* h1d = new TH1D("h1d", "Brho at FDC2", 100, 5, 8);
+ //for (auto i=0; i<Brho.size(); i++) h1d->Fill(Brho[i]);
+ //h1d->Draw();
+ //tc->Print((path+set+"histoBrho_FDC2.pdf").c_str());
+
+ // cout << "Brho\n";
+ // cout << "mean = " << mean(Brho) << endl;
+ // cout << "std dev = " << sqrt(var(Brho)) << endl;
+ //
+ // cout << "Brho Reconstructed\n";
+ // cout << "mean = " << mean(Br_rec) << endl;
+ // cout << "std dev = " << sqrt(var(Br_rec)) << endl;
+
+}
+
+void Brhos(string path_to_file){
+ /*string path = "";
+ string set = "spm1000/1cm_E0/";
+ string filename = path + set + "testanalysis.root";*/
+
+ //cout << "Read from " << path + set << endl;
+ cout <<"Read from "<<path_to_file << endl;
+
+ string filename = path_to_file + "testanalysis.root";
+
+
+ string fieldmap_file = "../../field_map/3T.table.bin";
+ //string fieldmap_file = path_to_fieldmap + "3T.table.bin";
+ TFile* file = new TFile(filename.c_str(), "read");
+ TTree* tree = (TTree*)file->Get("SimulatedTree");
+
+ SamuraiFieldMap* Map = new SamuraiFieldMap();
+ TSamuraiIdealData* FDC1 = new TSamuraiIdealData();
+ TSamuraiIdealData* FDC2 = new TSamuraiIdealData();
+ TSamuraiIdealData* Magnet = new TSamuraiIdealData();
+
+ tree->SetBranchAddress("IdealDataFDC1", &FDC1);
+ tree->SetBranchAddress("IdealDataFDC2", &FDC2);
+ tree->SetBranchAddress("IdealDataMagnet", &Magnet);
+
+ vector <double> Brho;
+ vector <TVector3> Pos_FDC1, Dir_FDC1, Pos_FDC2, Dir_FDC2;
+
+ int entries = tree->GetEntries();
+
+ TVector3 Mag_Pos (0,0,10000*mm);
+ TVector3 dir1, dir2;
+ for (int i=0; i<entries; i++){
+ tree->GetEntry(i);
+ if (FDC1->GetMult() == FDC2->GetMult()){
+ for (auto j=0;j<FDC2->GetMult(); j++){
+ Brho.push_back(FDC2->GetBrho(j));
+ Pos_FDC1.push_back(TVector3(FDC1->GetPosX(j), FDC1->GetPosY(j), FDC1->GetPosZ(j) ) - Mag_Pos);
+ Pos_FDC2.push_back(TVector3(FDC2->GetPosX(j), FDC2->GetPosY(j), FDC2->GetPosZ(j) ) - Mag_Pos);
+ dir1.SetMagThetaPhi(FDC1->GetMomMag(j), FDC1->GetMomTheta(j), FDC1->GetMomPhi(j));
+ Dir_FDC1.push_back(dir1.Unit());
+ dir2.SetMagThetaPhi(FDC2->GetMomMag(i), FDC2->GetMomTheta(j), FDC2->GetMomPhi(j));
+ Dir_FDC2.push_back(dir2.Unit());
+ }
+ }
+ }
+ //int time_int = 1000;
+ double mag_angle = 30*deg, fdc2_angle=-(59.930-30)*deg;
+ Map->LoadMap(mag_angle, fieldmap_file, 10);
+ Map->SetFDC2R((5000-438)*mm); //z pos - half thickness of FDC2
+ Map->SetFDC2Angle(fdc2_angle);
+ Map->SetStepLimit(1e5);
+
+ vector <int> time {10000,3000,1000,300,100};
+ //vector <int> time {10000,3000};
+ for(auto t:time) findbrho(path_to_file, t, Map, mag_angle,
+ fdc2_angle, Brho, Pos_FDC1, Pos_FDC2, Dir_FDC1, Dir_FDC2);
+////////////////////////////////////////////////////////////
+
+
+}
+
+void findbrho(string path_to_file, int time_int, SamuraiFieldMap* Map, double mag_angle, double fdc2_angle,
+ const vector<double>& Brho, const vector<TVector3>& Pos_FDC1,
+ const vector<TVector3>& Pos_FDC2, const vector<TVector3>& Dir_FDC1,
+ const vector<TVector3>& Dir_FDC2)
+{
+ Map->SetTimeIntervalSize(time_int*picosecond);
+
+ TVector3 posF2, dirF2;
+ vector <double> Br_rec;
+ for (auto i=0; i<Brho.size();i++){
+ cout<<i<<" ";
+ posF2 = Pos_FDC2[i];
+ posF2.RotateY(mag_angle-fdc2_angle);
+ dirF2 = Dir_FDC2[i];
+ dirF2.RotateY(mag_angle-fdc2_angle);
+ double brs = Map->FindBrho(Pos_FDC1[i], Dir_FDC1[i], posF2, dirF2);
+ Br_rec.push_back(brs);
+ }
+ double brho;
+ TFile* file2 = new TFile((path_to_file + "brho_" + to_string(time_int) + "ps.root").c_str(), "recreate");
+ TTree* brhos = new TTree("tree", "brhos");
+ brhos->Branch("Brho", &brho, "brho/D");
+ for(auto i=0; i<Br_rec.size(); i++) {
+ brho = Br_rec[i];
+ brhos->Fill();
+ }
+ file2->Write();
+ file2->Close();
+}
\ No newline at end of file
diff --git a/Projects/Samurai/benchmarkMacro.C b/Projects/Samurai/benchmarkMacro.C
new file mode 100644
index 0000000000000000000000000000000000000000..87db8f462b7134c397fc4d90ea016a7b4cc5e838
--- /dev/null
+++ b/Projects/Samurai/benchmarkMacro.C
@@ -0,0 +1,129 @@
+
+void xy_graph(vector<double> x, vector <double> y);
+
+double mean(vector<double>V){
+ double s = 0;
+ for (auto x:V) s += x;
+ return s / V.size();
+}
+double var(vector<double>V){
+ double s = 0, m=mean(V);
+ for (auto x:V) s += (x-m)*(x-m);
+ return s / (V.size()-1);
+}
+
+
+void benchmarkMacro (){
+ int n_file = 20;
+ string path = "root/simulation/bm1/";
+ string set = "EO/spm3/";
+ cout << "Read from " << path + set << endl;
+ TChain* chain = new TChain("SimulatedTree");
+ for (auto i=0; i<n_file; i++){
+ string filename = path + set + "benchmark" + to_string(i+1) + ".root";
+ chain->Add(filename.c_str());
+ }
+
+ TSamuraiIdealData* FDC1 = new TSamuraiIdealData();
+ TSamuraiIdealData* FDC2 = new TSamuraiIdealData();
+ TSamuraiIdealData* Magnet = new TSamuraiIdealData();
+
+ chain->SetBranchAddress("IdealDataFDC1", &FDC1);
+ chain->SetBranchAddress("IdealDataFDC2", &FDC2);
+ chain->SetBranchAddress("IdealDataMagnet", &Magnet);
+
+ vector <double> X, Y;
+
+ int entries = chain->GetEntries();
+ for (int i=0; i<entries; i++){
+ chain->GetEntry(i);
+ for (auto j=0;j<FDC2->GetMult(); j++){
+ X.push_back(FDC2->GetPosX(j));
+ Y.push_back(FDC2->GetPosY(j));
+ //Y.push_back(FDC2->GetPosZ(j));
+ //Y.push_back(FDC2->GetBrho(j));
+ }
+
+ }
+
+ TCanvas* tc = new TCanvas();
+
+ TH1D* h1d = new TH1D("h1d", "X at FDC2", 100, -3812, -3794);
+ for (auto i=0; i<X.size(); i++) h1d->Fill(X[i]);
+ h1d->Draw();
+ tc->Print((path+set+"histoX_FDC2.pdf").c_str());
+
+ TH2D* h2d = new TH2D("h2d", "X:Y", 100, -3812, -3794, 100, -20, 20);
+ for (int i=0; i<X.size(); i++) h2d->Fill(X[i], Y[i]);
+ h2d->Draw("colz");
+ tc->Print((path+set+"histoXY_FDC2.pdf").c_str());
+
+ cout << "X\n";
+ cout << "mean = " << mean(X) << endl;
+ cout << "std dev = " << sqrt(var(X)) << endl;
+ cout << "Y\n";
+ cout << "mean = " << mean(Y) << endl;
+ cout << "std dev = " << sqrt(var(Y)) << endl;
+
+
+}
+
+//-------------------------------------------------------------
+
+
+//-------------------------------------------------------------
+void t_hist(vector<double> t){
+
+ TH1D* hist = new TH1D("h1", "h1 title", 50, 0, 1);
+ for (unsigned int i = 0; i < t.size(); i++) hist->Fill(t[i]);
+
+ TF1* tf = new TF1("tf", "[0]", 0, 1);
+ tf->SetParameters(1, 1);
+ tf->SetParNames("A", "B");
+ hist->Fit(tf);
+
+ cout << "chi square : " << tf->GetChisquare() << endl;
+ cout << "prob : " << tf->GetProb() << endl;
+
+ TCanvas* tc = new TCanvas();
+ hist->GetYaxis()->SetRangeUser(0, 1.5*t.size()/50);
+ hist->Draw();
+ tc->Print("hist.pdf");
+
+ return;
+}
+
+void xy_graph(vector <double> x, vector <double> y){
+
+ TGraph* tg = new TGraph(x.size(), &x[0], &y[0]);
+
+ //TF1* tf = new TF1("tf", "pol2", 0, 1);
+ //tf->SetParameters(1, 1, 1);
+ //tf->SetParNames("A", "B", "C");
+ //tg->Fit(tf);
+
+ tg->SetMarkerStyle('+'); //draws x
+
+ //TCanvas* tc = new TCanvas();
+ tg->Draw("AP");
+ //tc->Print("graph.pdf");
+
+ return;
+}
+
+void tx_graperrors(vector <double> t, vector <double> x, vector <double> et, vector <double> ex){
+
+ TGraphErrors* tg = new TGraphErrors(t.size(), &t[0], &x[0], &et[0], &ex[0]);
+
+ TF1* tf = new TF1("tf", "pol2", 0, 1);
+ tf->SetParameters(1, 1, 1);
+ tf->SetParNames("A", "B", "C");
+ tg->Fit(tf);
+
+ cout << "chi square : " << tf->GetChisquare() << endl;
+ cout << "prob : " << tf->GetProb() << endl;
+
+ TCanvas* tc = new TCanvas();
+ tg->Draw("AP");
+ tc->Print("grapherrors.pdf");
+}
diff --git a/Projects/Samurai/run1.mac b/Projects/Samurai/run1.mac
index 0334083983567840763276104873862dd9f35d2c..e944da9d7a48206b50f7233fae95554919cfa3ec 100644
--- a/Projects/Samurai/run1.mac
+++ b/Projects/Samurai/run1.mac
@@ -1,3 +1,3 @@
-/run/beamOn 1
-/vis/scene/add/axes 0 0 0 5 m
-/vis/viewer/set/viewpointThetaPhi 90 90
+/run/beamOn 10000
+#/vis/scene/add/axes 0 0 0 5 m
+#/vis/viewer/set/viewpointThetaPhi 90 90