diff --git a/NPAnalysis/10He_Riken/Analysis b/NPAnalysis/10He_Riken/Analysis
index f5b004c15c19c002bffe3614626ee9f880d82d44..73acf47f2d6e340653baac089ce92f0d1716e695 100755
Binary files a/NPAnalysis/10He_Riken/Analysis and b/NPAnalysis/10He_Riken/Analysis differ
diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh
index f42981842a9a40aebd4476f7087cedfc4e020c5d..8cef0cdee67777ad8b16301617c9a42042b4adc7 100644
--- a/NPAnalysis/10He_Riken/include/ObjectManager.hh
+++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh
@@ -31,6 +31,7 @@
#include "NPReaction.h"
#include "RootInput.h"
#include "RootOutput.h"
+#include "TInitialConditions.h"
// Use CLHEP System of unit and Physical Constant
#include "CLHEP/Units/GlobalSystemOfUnits.h"
@@ -104,10 +105,20 @@ namespace ENERGYLOSS
{
// Declare your Energy loss here :
- EnergyLoss He3Target = EnergyLoss ( "3He_Mylar.txt" ,
- 100 ,
- 1 ,
- 3 );
+ EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" ,
+ 1000 ,
+ 1 ,
+ 3 );
+
+ EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2gaz_1b_26K.txt" ,
+ 100 ,
+ 1 ,
+ 3 );
+
+ EnergyLoss He3StripAl = EnergyLoss ( "3He_Al.txt" ,
+ 100 ,
+ 1 ,
+ 3 );
}
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index e23f787f3ca8e0b569b3066483ca515d8356bce2..a88cb3a92c65fd92962699030df914937fb7478d 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -31,50 +31,131 @@ int main(int argc,char** argv)
myDetector -> ReadConfigurationFile(detectorfileName) ;
// Attach more branch to the output
- double Ex = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
+ double ThinSi=-1 ;double Ex = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ; double ResolTheta=0;
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&Ex,"Ex/D") ;
RootOutput::getInstance()->GetTree()->Branch("E",&EE,"EE/D") ;
RootOutput::getInstance()->GetTree()->Branch("A",&TT,"TT/D") ;
RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ;
RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ThinSi_E",&ThinSi,"ThinSi/D") ;
+
+ // Get the formed Chained Tree and Treat it
+ TChain* Chain = RootInput:: getInstance() -> GetChain() ;
// Open the ThinSi Branch
- RootInput::getInstance() -> GetChain()->SetBranchStatus("ThinSiEnergy",true) ;
- double ThinSi=-1 ;
- RootInput::getInstance() -> GetChain()->SetBranchAddress("ThinSiEnergy",&ThinSi) ;
+ Chain->SetBranchStatus("ThinSiEnergy",true) ;
+ Chain->SetBranchStatus("InitialConditions",true) ;
+ Chain->SetBranchStatus("fIC_*",true) ;
+
+ TInitialConditions* Init = new TInitialConditions();
+ Chain->SetBranchAddress("ThinSiEnergy" ,&ThinSi );
+ Chain->SetBranchAddress("InitialConditions" ,&Init );
+
+ double TargetX=0 ; double TargetY=0; double BeamTheta = 0 ; double BeamPhi = 0 ;
+double TrueE=0 ; double TrueTheta=0 ;
+
// Get Must2 Pointer:
MUST2Array* M2 = (MUST2Array*) myDetector -> m_Detector["MUST2"] ;
- // Get the formed Chained Tree and Treat it
- TChain* Chain = RootInput:: getInstance() -> GetChain() ;
+
int i;
for ( i = 0 ; i < Chain -> GetEntries() ; i ++ )
{
if( i%10000 == 0 && i!=0) cout << i << " Event annalysed " << endl ;
Chain -> GetEntry(i);
-
myDetector -> ClearEventPhysics() ;
myDetector -> BuildPhysicalEvent() ;
-
double E = M2 -> GetEnergyDeposit();
- TVector3 A = M2 -> GetPositionOfInteraction();
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction() - TVector3(Init->GetICPositionX(0),Init->GetICPositionY(0),0);
- if(ThinSi > 0) E = E + ThinSi ;
-
- double Theta = ThetaCalculation ( A , TVector3(0,0,1) ) ;
-
- E= He3Target.EvaluateInitialEnergy( E , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- Theta );
+ BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ; BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
+
+// TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ;
+ TVector3 BeamDirection = TVector3(0,0,1) ; BeamDirection.SetTheta(BeamTheta) ; BeamDirection.SetPhi(BeamPhi) ;
+ double Theta = ThetaCalculation ( HitDirection , BeamDirection ) ;
- if(E>-1000) Ex = myReaction -> ReconstructRelativistic( E , Theta ) ;
- else Ex = -100 ;
+ if(E>-1000 && ThinSi>0 )
+ {
+
+ E=E+ThinSi;
+
+// E= He3StripAl.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 2*2*0.4*micrometer , // One for ThinSi and one for Must
+// 0 );
+
+// E= He3StripAl.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 2*0.4*micrometer , // Target Thickness at 0 degree
+// 0 );
+//
+// E= He3TargetWind.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 2*15*micrometer , // Target Thickness at 0 degree
+// Theta );
+
+// E= He3TargetGaz.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 3*mm , // Target Thickness at 0 degree
+// Theta );
+
+
+
+
+ E=E+ThinSi;
+
+ Ex = myReaction -> ReconstructRelativistic( E , Theta ) ;
+ X = HitDirection . X();
+ Y = HitDirection . Y();
+ }
+
+ else if(ThinSi>0)
+ {
+
+// ThinSi= He3StripAl.EvaluateInitialEnergy( ThinSi , // Energy of the detected particle
+// 2*0.4*micrometer , // Target Thickness at 0 degree
+// 0 );
+//
+// ThinSi= He3TargetWind.EvaluateInitialEnergy( ThinSi , // Energy of the detected particle
+// 2*15*micrometer , // Target Thickness at 0 degree
+// Theta );
+//
+// ThinSi= He3TargetGaz.EvaluateInitialEnergy( ThinSi , // Energy of the detected particle
+// 3*mm , // Target Thickness at 0 degree
+// Theta );
+
+ E= ThinSi;
+
+ Ex = myReaction -> ReconstructRelativistic( E , Theta ) ;
+ X = HitDirection . X();
+ Y = HitDirection . Y();
+
+ }
+ if(E>-1000 )
+ {
+
+
+
+// E= He3StripAl.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 2*0.4*micrometer ,
+// 0 );
+//
+// E= He3TargetWind.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 2*15*micrometer , // Target Thickness at 0 degree
+// Theta );
+//
+// E= He3TargetGaz.EvaluateInitialEnergy( E , // Energy of the detected particle
+// 3*mm , // Target Thickness at 0 degree
+// Theta );
+
+ ResolTheta = Theta/deg - Init->GetICEmittedAngleThetaLabWorldFrame(0);
+ Ex = myReaction -> ReconstructRelativistic( E, Theta ) ;
+ X = HitDirection . X();
+ Y = HitDirection . Y();
+
+ }
+
+ else {Ex=-100 ; X = -100 ; Y = -100 ;}
+
EE = E ; TT = Theta/deg ;
- if(E>-1000){
- X = A . X();
- Y = A . Y();}
- else{X = -1000 ; Y = -1000;}
RootOutput::getInstance()->GetTree()->Fill() ;
ThinSi = -1 ;
@@ -89,6 +170,6 @@ int main(int argc,char** argv)
double ThetaCalculation (TVector3 A , TVector3 B)
{
double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ;
- return Theta ;
+ return Theta*rad ;
}
diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h
index 86260d82f9b683a8ffc75c4e7866ad700236d5fb..86db26ad31b7f5faaea8198d0e00481bb2e471bc 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -94,7 +94,7 @@ public:
// Vertex of interaction
Double_t GetICPositionX(Int_t i) {return fIC_Position_X.at(i);}
Double_t GetICPositionY(Int_t i) {return fIC_Position_Y.at(i);}
- Double_t GetICPositionZ(Int_t i) {return fIC_Position_Z.at(i);}
+ Double_t GetICPositionZ(Int_t i) {return fIC_Position_Z.at(i);}
// Incident particle angles
Double_t GetICIncidentAngleTheta(Int_t i){return fIC_Incident_Angle_Theta.at(i);}
Double_t GetICIncidentAnglePhi(Int_t i) {return fIC_Incident_Angle_Phi.at(i);}
diff --git a/NPLib/MUST2/Must2Array.cxx b/NPLib/MUST2/Must2Array.cxx
index 26f51bae47a771c256e11d79998a6ac3247b0616..8041b3ae25d1a641577b3752687fe79941fbe360 100644
--- a/NPLib/MUST2/Must2Array.cxx
+++ b/NPLib/MUST2/Must2Array.cxx
@@ -373,10 +373,10 @@ void MUST2Array::AddTelescope( TVector3 C_X1_Y1 ,
NumberOfTelescope++;
// Vector U on Telescope Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
- TVector3 U = C_X128_Y1 - C_X1_Y1 ;
- U = -U.Unit() ;
+ TVector3 U = C_X1_Y1 - C_X128_Y1 ;
+ U = U.Unit() ;
// Vector V on Telescope Face (parallele to X Strip)
- TVector3 V = C_X1_Y128 - C_X1_Y1 ;
+ TVector3 V = C_X128_Y128 - C_X128_Y1 ;
V = V.Unit() ;
// Position Vector of Strip Center
diff --git a/NPLib/Tools/NPEnergyLoss.cxx b/NPLib/Tools/NPEnergyLoss.cxx
index d478a388fb6903425f7c49b5668cde500fa7e964..4e14530371ac83ed11e306f0f4a6447fa25177e4 100644
--- a/NPLib/Tools/NPEnergyLoss.cxx
+++ b/NPLib/Tools/NPEnergyLoss.cxx
@@ -236,7 +236,7 @@ double EnergyLoss::EvaluateInitialEnergy( double Energy , // Energy of the de
Energy = Energy / (double) fNumberOfMass ;
if (Angle > halfpi) Angle = pi-Angle ;
- TargetThickness = TargetThickness / ( 2*cos(Angle) ) ;
+ TargetThickness = TargetThickness / ( cos(Angle) ) ;
double SliceThickness = TargetThickness / (double)fNumberOfSlice ;
Interpolator* s = new Interpolator( fEnergy , fdEdX_Total ) ;
diff --git a/NPLib/Tools/NPEnergyLoss.h b/NPLib/Tools/NPEnergyLoss.h
index 7ec6bb53f5e46923fa0ec4e27ed913b099246877..6388a12221827576ba9253265bdf0ad41b8ed48f 100644
--- a/NPLib/Tools/NPEnergyLoss.h
+++ b/NPLib/Tools/NPEnergyLoss.h
@@ -53,7 +53,7 @@ namespace NPL
const;
// Evaluate Initial Energy of particle before crossing material knowing Angle, final Energy
- // and Target Thickness. Interaction is supposed to be in the middle of Target.
+ // and Target Thickness.
double EvaluateInitialEnergy( double energy , // Energy of the detected particle
double TargetThickness , // Target Thickness at 0 degree
double Angle ) // Particle Angle
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 018a695f4620ebdf9485585f55559c5e7730f1db..85b6a3fe84a204d0cae5a699a65c2596aeb1ff35 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -334,8 +334,10 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
G4double Ydir = cos( pi/2. - Beam_phiY ) ;
G4double Zdir = sin( pi/2. - Beam_thetaX ) + sin( pi/2. - Beam_phiY) ;
+
G4double Beam_theta = acos ( Zdir / sqrt( Xdir*Xdir + Ydir*Ydir + Zdir*Zdir ) );
- G4double Beam_phi = atan2( Ydir , Xdir );
+
+ G4double Beam_phi = atan2( Ydir , Xdir ) ;
// write angles to ROOT file
m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
@@ -394,8 +396,7 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa
// write angles/energy to ROOT file
m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(ThetaLight / deg);
- m_InitConditions->SetICEmittedEnergy(EnergyLight);
-
+ m_InitConditions->SetICEmittedEnergy(EnergyLight/MeV);
//must shoot inside the target.
G4double z0 = (-m_TargetThickness / 2 + RandFlat::shoot() * m_TargetThickness);
diff --git a/NPSimulation/src/ThinSi.cc b/NPSimulation/src/ThinSi.cc
index 432792a656dd312a37724d483a5eaecf7cc7859f..88e236f3fbb5029005967fb0747d278519f864d1 100644
--- a/NPSimulation/src/ThinSi.cc
+++ b/NPSimulation/src/ThinSi.cc
@@ -46,7 +46,7 @@ namespace THINSI
// Geometry
const G4double DetectorSize = 68*mm ;
const G4double SiliconThickness = 20*micrometer ;
- const G4double FrameThickness = 1*mm ;
+ const G4double FrameThickness = 3*mm ;
const G4double SiliconSize = 50*mm ;
const G4double AluThickness = 0.4*micrometer ;
const G4int NumberOfStrip = 32 ;
@@ -589,7 +589,7 @@ void ThinSi::ReadSensitive(const G4Event* event)
{
G4String DetectorNumber ;
bool checkSi = false ;
-
+ m_Energy = 0 ;
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////// Used to Read Event Map of detector //////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
@@ -598,6 +598,8 @@ void ThinSi::ReadSensitive(const G4Event* event)
std::map<G4int, G4double*>::iterator Energy_itr ;
G4THitsMap<G4double>* EnergyHitMap ;
+
+
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
G4int HitNumber = 0;