diff --git a/Benchmarks/physics/physics.cxx b/Benchmarks/physics/physics.cxx
index 10db7db4439e808982a9b0320179fc3ac045100e..feae0da77b57589985fac93780315983805def27 100644
--- a/Benchmarks/physics/physics.cxx
+++ b/Benchmarks/physics/physics.cxx
@@ -31,12 +31,13 @@ void EnergyLoss(){
}
////////////////////////////////////////////////////////////////////////////////
void Reaction() {
- unsigned int cycles = 1000000;
+ unsigned int cycles = 10000;
// test Random Ex generation
NPL::Reaction r2;
r2.ReadConfigurationFile("test.reaction");
TH1F* h2 = new TH1F("hE2","hE2",1000,-1,1);
double E4,T4,E3,T3,Ex,E;
+ r2.ShootRandomExcitationEnergy();
clock_t begin = clock();
for (unsigned int i = 0 ; i < cycles ; i++){
r2.ShootRandomExcitationEnergy();
diff --git a/Inputs/EventGenerator/11Li.beam b/Inputs/EventGenerator/11Li.beam
index 361a7b898c8924347fa6e013af88ea4032d4bbee..cd92cf20ee125da3baa228a1304b6f3f139cea6a 100644
--- a/Inputs/EventGenerator/11Li.beam
+++ b/Inputs/EventGenerator/11Li.beam
@@ -2,7 +2,7 @@
%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
- Particle= 11Li
+ Particle= 10He
Energy= 553
SigmaEnergy= 20
SigmaThetaX= 0.6138
diff --git a/Inputs/EventGenerator/Example1.reaction b/Inputs/EventGenerator/Example1.reaction
index 1b69a917b0a96e05b4c804ebb4f4362a8f2b75a9..24c22817ff15ddffa0fe5384736a43bcdcb0f49e 100644
--- a/Inputs/EventGenerator/Example1.reaction
+++ b/Inputs/EventGenerator/Example1.reaction
@@ -30,9 +30,12 @@ TwoBodyReaction
ShootHeavy= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%ParticleDecay 10He
-% Daughter= 8He n n
-% ExcitationEnergy= 0 0 0 MeV
-% shoot= 1 1 1
-
+Decay 10He
+ Daughter= 8He n n
+ ExcitationEnergy= 0 0 0 MeV
+ Threshold= 0 MeV
+ BranchingRatio= 0.5
+ LifeTime= 0 ns
+ Shoot= 1 1 1
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index ea19e4bc82f69f9e898e26f1beadd9a70658a35b..2f871a6b62f1646ac4dca1715f4a8f6fa5289520 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 5.5 MeV
- EnergyHigh= 5.5 MeV
+ EnergyLow= 3.5 MeV
+ EnergyHigh= 3.5 MeV
HalfOpenAngleMin= 0 deg
- HalfOpenAngleMax= 90 deg
+ HalfOpenAngleMax= 180 deg
x0= 0 mm
y0= 0 mm
z0= -110 mm
diff --git a/NPLib/Core/NPInputParser.cxx b/NPLib/Core/NPInputParser.cxx
index 577983975a9a8249375ecabbd9f3ded97dcc1502..5f1b3fa52ae586d104fdd29275462011af779186 100644
--- a/NPLib/Core/NPInputParser.cxx
+++ b/NPLib/Core/NPInputParser.cxx
@@ -331,6 +331,18 @@ std::vector<std::string> NPL::InputParser::GetAllBlocksToken(){
return token;
}
+////////////////////////////////////////////////////////////////////////////////
+std::vector<std::string> NPL::InputParser::GetAllBlocksValues(std::string token){
+ std::vector<std::string> value;
+ std::vector<NPL::InputBlock*> blocks = GetAllBlocksWithToken(token);
+ unsigned int size = blocks.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ value.push_back(blocks[i]->GetMainValue());
+ }
+
+ return value;
+}
+
////////////////////////////////////////////////////////////////////////////////
void NPL::InputParser::TreatAliases(){
diff --git a/NPLib/Core/NPInputParser.h b/NPLib/Core/NPInputParser.h
index cb808f59138af230b54f8567e7a54d28e8148d65..0791fe5c74a8c47f83b8bfd0390d79066d4dbc61 100644
--- a/NPLib/Core/NPInputParser.h
+++ b/NPLib/Core/NPInputParser.h
@@ -110,7 +110,7 @@ namespace NPL{
void Clear();
std::vector<InputBlock*> GetAllBlocksWithToken(std::string Token);
std::vector<InputBlock*> GetAllBlocksWithTokenAndValue(std::string Token,std::string Value);
-
+ std::vector<std::string> GetAllBlocksValues(std::string);
std::vector<std::string> GetAllBlocksToken();
diff --git a/NPLib/Detectors/CATS/TCATSPhysics.cxx b/NPLib/Detectors/CATS/TCATSPhysics.cxx
index 7cfc5aa9ed6bf9e24b269668bfdd66663fe3d4de..916cb16b978888725d54ed8f57efb4739f35c2b8 100644
--- a/NPLib/Detectors/CATS/TCATSPhysics.cxx
+++ b/NPLib/Detectors/CATS/TCATSPhysics.cxx
@@ -49,6 +49,7 @@ ClassImp(TCATSPhysics)
m_EventPhysics = this ;
m_NumberOfCATS = 0 ;
m_Spectra = NULL ;
+ m_Zproj = 0 ;
}
///////////////////////////////////////////////////////////////////////////
@@ -108,7 +109,11 @@ void TCATSPhysics::BuildSimplePhysicalEvent(){
//////////////////////////////////////////////////////////////////////////////
void TCATSPhysics::BuildPhysicalEvent(){
+
+
+
PreTreat();
+
// Look how many CATS were fired
// use a set to identify which detector has been hit
set<int> DetectorHitX; // X only
@@ -164,7 +169,7 @@ void TCATSPhysics::BuildPhysicalEvent(){
if(NX == DetMaxX[j] ){
Buffer_X_Q[j][StrX-1]= CATS_X_Q;
QsumX[j]+= CATS_X_Q;
- if(CATS_X_Q > Buffer_X_Q[j][StripMaxX[j] -1]){
+ if(CATS_X_Q > Buffer_X_Q[j][StripMaxX[j]-1]){
StripMaxX[j] = StrX ;
ChargeMaxX[j]= CATS_X_Q;
}
@@ -185,7 +190,7 @@ void TCATSPhysics::BuildPhysicalEvent(){
if(NY == DetMaxY[j] ){
Buffer_Y_Q[j][StrY-1]= CATS_Y_Q;
QsumY[j]+= CATS_Y_Q;
- if(CATS_Y_Q > Buffer_Y_Q[j][StripMaxY[j] -1]){
+ if(CATS_Y_Q > Buffer_Y_Q[j][StripMaxY[j]-1]){
StripMaxY[j] = StrY ;
ChargeMaxY[j]= CATS_Y_Q;
}
@@ -208,8 +213,8 @@ void TCATSPhysics::BuildPhysicalEvent(){
// a shift - -1 is made to have PosX in between -0.5 and 27.5
// for the following calculation of the position in the lab.
- PosX = PosX -1;
- PosY = PosY -1;
+ PosX = PosX;
+ PosY = PosY;
// sx and sy are the X and Y strip number between which the PosX and PosY are
int sx0 = (int) PosX;
@@ -227,6 +232,9 @@ void TCATSPhysics::BuildPhysicalEvent(){
PositionX.push_back(px0+(px1-px0)*(PosX-sx0));
PositionY.push_back(py0+(py1-py0)*(PosY-sy0));
+ //PositionX.push_back(2.54*(PosX-14));
+ //PositionY.push_back(2.54*(PosY-14));
+
PositionZ.push_back(StripPositionZ[DetMaxX[i]-1]);
}
@@ -235,14 +243,23 @@ void TCATSPhysics::BuildPhysicalEvent(){
// At least two CATS need to gave back position in order to reconstruct on Target
if(PositionX.size()>1){
if(DetMaxX[0]<DetMaxX[1]){
- double t = -PositionZ[1]/(PositionZ[1]-PositionZ[0]);
+
+ // cout << "Test " << m_Zproj << endl ;
+ // cout << "Test2 " << PositionZ[0]<< endl ;
+ // cout << "Test3 " << PositionZ[1]<< endl ;
+ double t = (m_Zproj-PositionZ[1])/(PositionZ[1]-PositionZ[0]);
+ // cout << "t " << t << endl ;
+ // cout << "X1 " << PositionX[0] << endl ;
+ // cout << "X2 " << PositionX[1] << endl ;
PositionOnTargetX= PositionX[1] + (PositionX[1]-PositionX[0])*t;
PositionOnTargetY= PositionY[1] + (PositionY[1]-PositionY[0])*t;
+
+ //cout << "X3 " << PositionOnTargetX << endl ;
BeamDirection = GetBeamDirection();
}
else{
- double t = -PositionZ[0]/(PositionZ[0]-PositionZ[1]);
+ double t = (m_Zproj-PositionZ[1])/(PositionZ[0]-PositionZ[1]);
PositionOnTargetX= PositionX[0] + (PositionX[0]-PositionX[1])*t;
PositionOnTargetY= PositionY[0] + (PositionY[0]-PositionY[1])*t;
BeamDirection = GetBeamDirection();
@@ -564,6 +581,12 @@ void TCATSPhysics::ReadAnalysisConfig(){
SetReconstructionMethod(CATSNumber,XorY,DataBuffer);
}
+ else if (whatToDo == "ZPROJ") {
+ AnalysisConfigFile >> DataBuffer;
+ cout << whatToDo << " " << DataBuffer << endl ;
+ m_Zproj = atoi(DataBuffer.c_str());
+ }
+
else {ReadingStatus = false;}
}
@@ -852,7 +875,8 @@ namespace CATS_LOCAL{
name+= NPL::itoa( m_EventData->GetCATSStripX(i) ) ;
name+= "_Q";
return CalibrationManager::getInstance()->ApplyCalibration( name,
- m_EventData->GetCATSChargeX(i) + gRandom->Rndm() - fCATS_Ped_X(m_EventData, i) );
+ m_EventData->GetCATSChargeX(i) + gRandom->Rndm() );
+ //m_EventData->GetCATSChargeX(i) + gRandom->Rndm() - fCATS_Ped_X(m_EventData, i) );
}
////////////////////////////////////////////////////////////////////////
double fCATS_Y_Q(const TCATSData* m_EventData , const int i){
@@ -863,7 +887,8 @@ namespace CATS_LOCAL{
name+= NPL::itoa( m_EventData->GetCATSStripY(i) ) ;
name+= "_Q";
return CalibrationManager::getInstance()->ApplyCalibration( name ,
- m_EventData->GetCATSChargeY(i) + gRandom->Rndm() - fCATS_Ped_Y(m_EventData, i) );
+ m_EventData->GetCATSChargeY(i) + gRandom->Rndm() );
+ //m_EventData->GetCATSChargeY(i) + gRandom->Rndm() - fCATS_Ped_Y(m_EventData, i) );
}
////////////////////////////////////////////////////////////////////////
bool fCATS_Threshold_X(const TCATSData* m_EventData , const int i){
diff --git a/NPLib/Detectors/CATS/TCATSPhysics.h b/NPLib/Detectors/CATS/TCATSPhysics.h
index 64fd27f7c909b0489f51e7004b4aa302e30deddf..ab5b2860b69d3eca78dcec6548b99e53e60c6211 100644
--- a/NPLib/Detectors/CATS/TCATSPhysics.h
+++ b/NPLib/Detectors/CATS/TCATSPhysics.h
@@ -84,6 +84,7 @@ class TCATSPhysics : public TObject, public NPL::VDetector
vector<double> QsumY;
double PositionOnTargetX;
double PositionOnTargetY;
+ double m_Zproj;;
TVector3 BeamDirection;//!
@@ -167,7 +168,7 @@ class TCATSPhysics : public TObject, public NPL::VDetector
void Clear(const Option_t*) {};
// Give and external TCATSData object to TCATSPhysics, needed for online analysis
- void SetRawDataPointer(TCATSData* rawDataPointer) {m_EventData = rawDataPointer;}
+ void SetRawDataPointer(void* rawDataPointer) {m_EventData = (TCATSData*)rawDataPointer;}
// Return false if the channel is disabled by user
bool IsValidChannel(const string DetectorType, const int Detector , const int channel);
diff --git a/NPLib/Detectors/MUST2/TMust2Physics.cxx b/NPLib/Detectors/MUST2/TMust2Physics.cxx
index 8466864423ec2a06227af41cc302a594f73c27b2..7c5f7d735f685f378ddb4ca5662658a6d62a9f26 100644
--- a/NPLib/Detectors/MUST2/TMust2Physics.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Physics.cxx
@@ -76,7 +76,7 @@ ClassImp(TMust2Physics)
for(int i = 0 ; i < 16 ; ++i){
m_SiLi_MatchingX[0]=112;
- m_SiLi_MatchingY[0]=112;
+ m_SiLi_MatchingY[1]=112;
m_SiLi_MatchingX[1]=112;
m_SiLi_MatchingY[1]=80;
@@ -128,52 +128,52 @@ ClassImp(TMust2Physics)
m_CsI_MatchingX.resize(16,0);
m_CsI_MatchingY.resize(16,0);
for(int i = 0 ; i < 16 ; ++i){
- m_CsI_MatchingX[0]=80;
- m_CsI_MatchingY[0]=48;
+ m_CsI_MatchingX[0]=112;
+ m_CsI_MatchingY[0]=112;
m_CsI_MatchingX[1]=112;
- m_CsI_MatchingY[1]=48;
+ m_CsI_MatchingY[1]=80;
- m_CsI_MatchingX[2]=80;
- m_CsI_MatchingY[2]=16;
+ m_CsI_MatchingX[2]=112;
+ m_CsI_MatchingY[2]=48;
m_CsI_MatchingX[3]=112;
m_CsI_MatchingY[3]=16;
//
- m_CsI_MatchingX[4]=48;
- m_CsI_MatchingY[4]=80;
+ m_CsI_MatchingX[4]=80;
+ m_CsI_MatchingY[4]=16;
- m_CsI_MatchingX[5]=16;
- m_CsI_MatchingY[5]=80;
+ m_CsI_MatchingX[5]=80;
+ m_CsI_MatchingY[5]=48;
- m_CsI_MatchingX[6]=48;
- m_CsI_MatchingY[6]=112;
+ m_CsI_MatchingX[6]=80;
+ m_CsI_MatchingY[6]=80;
- m_CsI_MatchingX[7]=16;
+ m_CsI_MatchingX[7]=80;
m_CsI_MatchingY[7]=112;
//
m_CsI_MatchingX[8]=48;
m_CsI_MatchingY[8]=48;
- m_CsI_MatchingX[9]=16;
+ m_CsI_MatchingX[9]=48;
m_CsI_MatchingY[9]=48;
m_CsI_MatchingX[10]=48;
- m_CsI_MatchingY[10]=16;
+ m_CsI_MatchingY[10]=80;
- m_CsI_MatchingX[11]=16;
- m_CsI_MatchingY[11]=16;
+ m_CsI_MatchingX[11]=48;
+ m_CsI_MatchingY[11]=112;
//
- m_CsI_MatchingX[12]=80;
- m_CsI_MatchingY[12]=80;
+ m_CsI_MatchingX[12]=16;
+ m_CsI_MatchingY[12]=16;
- m_CsI_MatchingX[13]=112;
- m_CsI_MatchingY[13]=80;
+ m_CsI_MatchingX[13]=16;
+ m_CsI_MatchingY[13]=48;
- m_CsI_MatchingX[14]=80;
- m_CsI_MatchingY[14]=112;
+ m_CsI_MatchingX[14]=16;
+ m_CsI_MatchingY[14]=80;
- m_CsI_MatchingX[15]=112;
+ m_CsI_MatchingX[15]=16;
m_CsI_MatchingY[15]=112;
}
@@ -191,7 +191,9 @@ void TMust2Physics::BuildSimplePhysicalEvent(){
///////////////////////////////////////////////////////////////////////////
void TMust2Physics::BuildPhysicalEvent(){
+
PreTreat();
+
bool check_SILI = false ;
bool check_CSI = false ;
@@ -203,8 +205,9 @@ void TMust2Physics::BuildPhysicalEvent(){
m_SiLiTMult = m_PreTreatedData->GetMMSiLiTMult();
m_CsIEMult = m_PreTreatedData->GetMMCsIEMult();
m_CsITMult = m_PreTreatedData->GetMMCsITMult();
- if( CheckEvent() == 1 ){
+ if( 1 /*CheckEvent() == 1*/ ){
vector< TVector2 > couple = Match_X_Y() ;
+
EventMultiplicity = couple.size();
for(unsigned int i = 0 ; i < couple.size() ; ++i){
check_SILI = false ;
@@ -324,7 +327,7 @@ void TMust2Physics::PreTreat(){
m_SiLiTMult = m_EventData->GetMMSiLiTMult();
m_CsIEMult = m_EventData->GetMMCsIEMult();
m_CsITMult = m_EventData->GetMMCsITMult();
- // X
+ // X
// E
for(unsigned int i = 0 ; i < m_StripXEMult ; ++i){
if( m_EventData->GetMMStripXEEnergy(i)>m_Si_X_E_RAW_Threshold && IsValidChannel(0, m_EventData->GetMMStripXEDetectorNbr(i), m_EventData->GetMMStripXEStripNbr(i)) ){
@@ -421,14 +424,15 @@ bool TMust2Physics :: ResolvePseudoEvent(){
vector < TVector2 > TMust2Physics :: Match_X_Y(){
vector < TVector2 > ArrayOfGoodCouple ;
m_StripXEMult = m_PreTreatedData->GetMMStripXEMult();
- m_StripYEMult = m_PreTreatedData->GetMMStripXEMult();
-
+ m_StripYEMult = m_PreTreatedData->GetMMStripYEMult();
+
// Prevent code from treating very high multiplicity Event
// Those event are not physical anyway and that improve speed.
- if( m_StripXEMult > m_MaximumStripMultiplicityAllowed || m_StripYEMult > m_MaximumStripMultiplicityAllowed )
+ if( m_StripXEMult > m_MaximumStripMultiplicityAllowed || m_StripYEMult > m_MaximumStripMultiplicityAllowed ){
return ArrayOfGoodCouple;
+ }
- for(unsigned int i = 0 ; i < m_StripXEMult ; ++i){
+ for(unsigned int i = 0 ; i < m_StripXEMult ; i++){
for(unsigned int j = 0 ; j < m_StripYEMult ; j++){
// if same detector check energy
if ( m_PreTreatedData->GetMMStripXEDetectorNbr(i) == m_PreTreatedData->GetMMStripYEDetectorNbr(j) ){
@@ -449,7 +453,7 @@ vector < TVector2 > TMust2Physics :: Match_X_Y(){
}
// Special Option, if the event is between two SiLi pad , it is rejected.
- if(m_Ignore_not_matching_SiLi){
+ else if(m_Ignore_not_matching_SiLi){
bool check_validity=false;
for (unsigned int hh = 0 ; hh<16 ; ++hh ){
if( Match_Si_SiLi(m_PreTreatedData->GetMMStripXEStripNbr(i), m_PreTreatedData->GetMMStripYEStripNbr(j) , hh+1) )
@@ -461,7 +465,9 @@ vector < TVector2 > TMust2Physics :: Match_X_Y(){
}
// Regular case, keep the event
- else ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
+ else {
+ ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
+ }
}
}
}
@@ -850,7 +856,7 @@ void TMust2Physics::ReadConfiguration(NPL::InputParser parser){
InitializeStandardParameter();
ReadAnalysisConfig();
}
-///////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
void TMust2Physics::InitSpectra(){
m_Spectra = new TMust2Spectra(m_NumberOfTelescope);
}
@@ -968,18 +974,13 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1,
m_NumberOfTelescope++;
- // Geometry Parameter
- double Face = 97.3; //mm
- double NumberOfStrip = 128;
- double StripPitch = Face/NumberOfStrip ; //mm
-
// 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 ;
- double Ushift = (U.Mag()-Face)/2.;
+ double Ushift = (U.Mag()-98)/2.;
U = U.Unit();
// Vector V on Telescope Face (parallele to X Strip)
TVector3 V = C_X1_Y128 - C_X1_Y1 ;
- double Vshift = (V.Mag() -Face)/2. ;
+ double Vshift = (V.Mag() -98)/2. ;
V = V.Unit() ;
// Position Vector of Strip Center
@@ -987,6 +988,10 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1,
// Position Vector of X=1 Y=1 Strip
TVector3 Strip_1_1;
+ // Geometry Parameter
+ double Face = 98; //mm
+ double NumberOfStrip = 128;
+ double StripPitch = Face/NumberOfStrip ; //mm
// Buffer object to fill Position Array
vector<double> lineX ; vector<double> lineY ; vector<double> lineZ ;
@@ -1095,7 +1100,7 @@ void TMust2Physics::AddTelescope( double theta,
U.Rotate( beta_w * Pi/180. , W ) ;
V.Rotate( beta_w * Pi/180. , W ) ;
- double Face = 97.3 ; //mm
+ double Face = 98 ; //mm
double NumberOfStrip = 128 ;
double StripPitch = Face/NumberOfStrip ; //mm
@@ -1144,7 +1149,9 @@ TVector3 TMust2Physics::GetPositionOfInteraction(const int i) const{
TVector3 Position = TVector3 ( GetStripPositionX( TelescopeNumber[i] , Si_X[i] , Si_Y[i] ) ,
GetStripPositionY( TelescopeNumber[i] , Si_X[i] , Si_Y[i] ) ,
GetStripPositionZ( TelescopeNumber[i] , Si_X[i] , Si_Y[i] ) ) ;
+
return(Position) ;
+
}
TVector3 TMust2Physics::GetTelescopeNormal( const int i) const{
diff --git a/NPLib/Detectors/MUST2/TMust2Spectra.cxx b/NPLib/Detectors/MUST2/TMust2Spectra.cxx
index 0f29265db6591da7676a0b3bcd7e31bda9909086..aa5a33bf65be58d01a2856b75566b79511ed61fd 100644
--- a/NPLib/Detectors/MUST2/TMust2Spectra.cxx
+++ b/NPLib/Detectors/MUST2/TMust2Spectra.cxx
@@ -146,11 +146,11 @@ void TMust2Spectra::InitPreTreatedSpectra()
// STRX_T_CAL
name = "MM"+NPL::itoa(i+1)+"_STRX_T_CAL";
- AddHisto2D(name, name, fStripX, 1, fStripX+1, 500, 0, 500, "MUST2/CAL/STRXT");
+ AddHisto2D(name, name, fStripX, 1, fStripX+1, 1000, 0, 1000, "MUST2/CAL/STRXT");
// STRY_T_CAL
name = "MM"+NPL::itoa(i+1)+"_STRY_T_CAL";
- AddHisto2D(name, name, fStripY, 1, fStripY+1, 500, 0, 500, "MUST2/CAL/STRYT");
+ AddHisto2D(name, name, fStripY, 1, fStripY+1, 1000, 0, 1000, "MUST2/CAL/STRYT");
// SILI_E_CAL
name = "MM"+NPL::itoa(i+1)+"_SILI_E_CAL";
diff --git a/NPLib/Detectors/ModularLeaf/TModularLeafPhysics.h b/NPLib/Detectors/ModularLeaf/TModularLeafPhysics.h
index 3bcddfac2611567cd1a2d45e4b5a07fe73077153..83cf334201e0771947b678ae72a5ccdb9130c358 100644
--- a/NPLib/Detectors/ModularLeaf/TModularLeafPhysics.h
+++ b/NPLib/Detectors/ModularLeaf/TModularLeafPhysics.h
@@ -90,7 +90,11 @@ class TModularLeafPhysics : public TObject, public NPL::VDetector{
//TModularLeafData* EventData ;//!
TModularLeafPhysics* EventPhysics ;//!
- public: // Static constructor to be passed to the Detector Factory
+
+ public:
+ inline short GetRawValue(std::string label){return m_RawData[label];};
+ inline double GetCalibratedValue(std::string label){return m_CalibratedData[label];};
+ public: // Static constructor to be passed to the Detector Factory
static NPL::VDetector* Construct();
ClassDef(TModularLeafPhysics,1) // ModularLeafPhysics structure
};
diff --git a/NPLib/Physics/CMakeLists.txt b/NPLib/Physics/CMakeLists.txt
index 06d6713640464d23bb452541b91deaa5921414f8..e8863cec47184e174223181e68e9a87f532ac21d 100644
--- a/NPLib/Physics/CMakeLists.txt
+++ b/NPLib/Physics/CMakeLists.txt
@@ -3,7 +3,7 @@ add_custom_command(OUTPUT NPEnergyLossDict.cxx COMMAND ../scripts/build_dict.sh
add_custom_command(OUTPUT TInitialConditionsDict.cxx COMMAND ../scripts/build_dict.sh TInitialConditions.h TInitialConditionsDict.cxx TInitialConditions.rootmap libNPInitialConditions.so DEPENDS TInitialConditions.h)
add_custom_command(OUTPUT TInteractionCoordinatesDict.cxx COMMAND ../scripts/build_dict.sh TInteractionCoordinates.h TInteractionCoordinatesDict.cxx TInteractionCoordinates.rootmap libNPInteractionCoordinates.so DEPENDS TInteractionCoordinates.h)
-add_library(NPPhysics SHARED NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPNucleus.cxx NPReaction.cxx NPReactionDict.cxx NPEnergyLossDict.cxx )
+add_library(NPPhysics SHARED NPDecay.cxx NPBeam.cxx NPEnergyLoss.cxx NPFunction.cxx NPNucleus.cxx NPReaction.cxx NPReactionDict.cxx NPEnergyLossDict.cxx )
target_link_libraries(NPPhysics ${ROOT_LIBRARIES} Physics NPCore)
add_library(NPInitialConditions SHARED TInitialConditions.cxx TInitialConditionsDict.cxx )
@@ -12,4 +12,4 @@ target_link_libraries(NPInitialConditions ${ROOT_LIBRARIES} )
add_library(NPInteractionCoordinates SHARED TInteractionCoordinates.cxx TInteractionCoordinatesDict.cxx)
target_link_libraries(NPInteractionCoordinates ${ROOT_LIBRARIES} )
-install(FILES NPBeam.h NPEnergyLoss.h NPFunction.h NPNucleus.h NPReaction.h TInitialConditions.h TInteractionCoordinates.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES NPDecay.h NPBeam.h NPEnergyLoss.h NPFunction.h NPNucleus.h NPReaction.h TInitialConditions.h TInteractionCoordinates.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index 454f4f894dd7f3ac32ce35fdc92e11d85354c002..ab7007695d8d5b7d590bda7cba496613ec426e9a 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -66,6 +66,8 @@ Beam::Beam(){
fEffectiveTargetThickness = 0 ;
fTargetAngle = 0 ;
fTargetZ = 0 ;
+ fZEmission=-1*NPUNITS::m;
+ fZProfile=0;
fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
// case of user given distribution
@@ -98,6 +100,8 @@ Beam::Beam(string isotope){
fEffectiveTargetThickness = 0 ;
fTargetAngle = 0 ;
fTargetZ = 0 ;
+ fZEmission=-1*NPUNITS::m;
+ fZProfile=0;
fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
// case of user given distribution
@@ -141,6 +145,12 @@ void Beam::ReadConfigurationFile(NPL::InputParser parser){
if(blocks[i]->HasToken("ExcitationEnergy"))
fExcitationEnergy = blocks[i]->GetDouble("ExcitationEnergy","MeV");
+ if(blocks[i]->HasToken("ZEmission"))
+ fZEmission = blocks[i]->GetDouble("ZEmission","mm");
+
+ if(blocks[i]->HasToken("ZProfile"))
+ fZProfile = blocks[i]->GetDouble("ZProfile","mm");
+
// Energy analytic
if(blocks[i]->HasTokenList(energyA)){
fEnergy = blocks[i]->GetDouble("Energy","MeV");
@@ -175,7 +185,7 @@ void Beam::ReadConfigurationFile(NPL::InputParser parser){
SetXThetaXHist( Read2DProfile(XThetaX[0], XThetaX[1]));
vector<string> YPhiY= blocks[i]->GetVectorString("YPhiYProfilePath");
SetYPhiYHist( Read2DProfile(YPhiY[0], YPhiY[1]));
-
+
}
else{
@@ -195,28 +205,43 @@ void Beam::ReadConfigurationFile(NPL::InputParser parser){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY ){
- X = Y = 1000000*cm;
+ double X0,Y0;
+ X0 = Y0 = 1*km;
+ // ENERGY //
+ // Gaussian energy distribution
if(fSigmaEnergy!=-1)
E = gRandom->Gaus(fEnergy,fSigmaEnergy);
+ // User Profile
else
E = fEnergyHist->GetRandom();
+ // POSITION/DIRECTION AT Z PROFILE//
+ // Gaussian Distribution
if(fSigmaX!=-1){
- // Shoot within the target unless target size is null (no limit)
- while(sqrt(X*X+Y*Y)>fEffectiveTargetSize || fEffectiveTargetSize == 0){
- NPL::RandomGaussian2D(fMeanX, fMeanThetaX, fSigmaX, fSigmaThetaX, X, ThetaX);
- NPL::RandomGaussian2D(fMeanY, fMeanPhiY, fSigmaY, fSigmaPhiY, Y, PhiY);
- }
+ NPL::RandomGaussian2D(fMeanX, fMeanThetaX, fSigmaX, fSigmaThetaX, X0, ThetaX);
+ NPL::RandomGaussian2D(fMeanY, fMeanPhiY, fSigmaY, fSigmaPhiY, Y0, PhiY);
}
+ // Profile
else{
- while(sqrt(X*X+Y*Y)>fEffectiveTargetSize || fEffectiveTargetSize == 0){
fXThetaXHist->GetRandom2(X,ThetaX);
fYPhiYHist->GetRandom2(Y,PhiY);
- }
}
- Z = fTargetZ + fEffectiveTargetThickness*(gRandom->Uniform()-0.5);
+ // 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;
+ Z = fZEmission;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index 574937e1109eec1ac5980586a1033db183214b1b..38d2dcd3883f6af8248d788efe173ae3666276e5 100644
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -1,5 +1,5 @@
-#ifndef __Beam__
-#define __Beam__
+#ifndef NPBEAM
+#define NPBEAM
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
@@ -40,102 +40,109 @@ using namespace std;
using namespace NPL;
namespace NPL{
-
+
class Beam:public NPL::Nucleus{
-
- public: // Constructors and Destructors
- Beam();
- Beam(string);
- ~Beam();
-
- public: // Various Method
- void ReadConfigurationFile(string Path);
- void ReadConfigurationFile(NPL::InputParser);
-
- private:
- int fVerboseLevel;
-
- private:
- //Nucleus* fBeamNucleus;
- double fEnergy;
- double fExcitationEnergy;
- double fSigmaEnergy;
- double fMeanX;
- double fMeanY;
- double fSigmaX;
- double fSigmaY;
- double fMeanThetaX;
- double fMeanPhiY;
- double fSigmaThetaX;
- double fSigmaPhiY;
-
- // case of user given distribution
- TH1F* fEnergyHist;
- TH2F* fXThetaXHist;
- TH2F* fYPhiYHist;
-
- public:
- // Getters and Setters
- // Set
- // void SetBeamNucleus (Nucleus* BeamNucleus) {delete fBeamNucleus ; fBeamNucleus = new Nucleus(BeamNucleus->GetZ(),BeamNucleus->GetA());}
- void SetEnergy (double Energy) {fEnergy=Energy;}
- void SetSigmaEnergy (double SigmaEnergy) {fSigmaEnergy=SigmaEnergy;}
- void SetMeanX (double MeanX) {fMeanX=MeanX;}
- void SetMeanY (double MeanY) {fMeanY=MeanY;}
- void SetSigmaX (double SigmaX) {fSigmaX=SigmaX;}
- void SetSigmaY (double SigmaY) {fSigmaY=SigmaY;}
- void SetMeanThetaX (double MeanThetaX) {fMeanThetaX=MeanThetaX;}
- void SetMeanPhiY (double MeanPhiY) {fMeanPhiY=MeanPhiY;}
- void SetSigmaThetaX (double SigmaThetaX) {fSigmaThetaX=SigmaThetaX;}
- void SetSigmaPhiY (double SigmaPhiY) {fSigmaPhiY=SigmaPhiY;}
- void SetEnergyHist (TH1F* EnergyHist) {delete fEnergyHist; fEnergyHist = EnergyHist;}
- void SetXThetaXHist (TH2F* XThetaXHist) {delete fXThetaXHist; fXThetaXHist = XThetaXHist;}
- void SetYPhiYHist (TH2F* YPhiYHist) {delete fYPhiYHist; fYPhiYHist = YPhiYHist;}
- void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
-
- // Get
- // Nucleus* GetNucleus () const {return fBeamNucleus;}
- double GetEnergy () const {return fEnergy;}
- double GetExcitationEnergy() const {return fExcitationEnergy;}
- double GetSigmaEnergy () const {return fSigmaEnergy;}
- double GetMeanX () const {return fMeanX;}
- double GetMeanY () const {return fMeanY;}
- double GetSigmaX () const {return fSigmaX;}
- double GetSigmaY () const {return fSigmaY;}
- double GetMeanThetaX () const {return fMeanThetaX;}
- double GetMeanPhiY () const {return fMeanPhiY;}
- double GetSigmaThetaX () const {return fSigmaThetaX;}
- double GetSigmaPhiY () const {return fSigmaPhiY;}
- TH1F* GetEnergyHist () const {return fEnergyHist;}
- TH2F* GetXThetaXHist () const {return fXThetaXHist;}
- TH2F* GetYPhiYHist () const {return fYPhiYHist;}
- int GetVerboseLevel() const {return fVerboseLevel;}
-
- private: // Event Generation private variable
- double fTargetSize;
- double fEffectiveTargetSize; // target size has seen from the beam axis
- double fTargetThickness;
- double fEffectiveTargetThickness; // target thickness has seen by the beam
- double fTargetAngle;
- double fTargetZ;
-
- public:
- void SetTargetSize(double TargetSize);
- void SetTargetThickness(double TargetThickness);
- void SetTargetAngle(double TargetAngle);
- void SetTargetZ(double TargetZ) {fTargetZ = TargetZ;}
- double GetTargetSize(){return fTargetSize;}
- double GetTargetThickness(){return fTargetThickness;}
- double GetTargetAngle(){return fTargetAngle;}
- double GetTargetZ() {return fTargetZ;}
- double GetTargetEffectiveThickness(){return fEffectiveTargetThickness;}
- double GetTargetEffectiveTargetSize(){return fEffectiveTargetSize;}
-
-
- public: // Event Generation
- void GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY );
- public: // Print private paremeter
- void Print() const;
+
+ public: // Constructors and Destructors
+ Beam();
+ Beam(string);
+ ~Beam();
+
+ public: // Various Method
+ void ReadConfigurationFile(string Path);
+ void ReadConfigurationFile(NPL::InputParser);
+
+ private:
+ int fVerboseLevel;
+
+ private:
+ //Nucleus* fBeamNucleus;
+ double fEnergy;
+ double fExcitationEnergy;
+ double fSigmaEnergy;
+ double fMeanX;
+ double fMeanY;
+ double fSigmaX;
+ double fSigmaY;
+ double fMeanThetaX;
+ double fMeanPhiY;
+ double fSigmaThetaX;
+ double fSigmaPhiY;
+ double fZProfile;
+ double fZEmission;
+
+ // case of user given distribution
+ TH1F* fEnergyHist;
+ TH2F* fXThetaXHist;
+ TH2F* fYPhiYHist;
+
+ public:
+ // Getters and Setters
+ // Set
+ // void SetBeamNucleus (Nucleus* BeamNucleus) {delete fBeamNucleus ; fBeamNucleus = new Nucleus(BeamNucleus->GetZ(),BeamNucleus->GetA());}
+ inline void SetEnergy (double& Energy) {fEnergy=Energy;}
+ inline void SetSigmaEnergy (double& SigmaEnergy) {fSigmaEnergy=SigmaEnergy;}
+ inline void SetMeanX (double& MeanX) {fMeanX=MeanX;}
+ inline void SetMeanY (double& MeanY) {fMeanY=MeanY;}
+ inline void SetSigmaX (double& SigmaX) {fSigmaX=SigmaX;}
+ inline void SetSigmaY (double& SigmaY) {fSigmaY=SigmaY;}
+ inline void SetMeanThetaX (double& MeanThetaX) {fMeanThetaX=MeanThetaX;}
+ inline void SetMeanPhiY (double& MeanPhiY) {fMeanPhiY=MeanPhiY;}
+ inline void SetSigmaThetaX (double& SigmaThetaX) {fSigmaThetaX=SigmaThetaX;}
+ inline void SetSigmaPhiY (double& SigmaPhiY) {fSigmaPhiY=SigmaPhiY;}
+ inline void SetEnergyHist (TH1F* EnergyHist) {delete fEnergyHist; fEnergyHist = EnergyHist;}
+ inline void SetXThetaXHist (TH2F* XThetaXHist) {delete fXThetaXHist; fXThetaXHist = XThetaXHist;}
+ inline void SetYPhiYHist (TH2F* YPhiYHist) {delete fYPhiYHist; fYPhiYHist = YPhiYHist;}
+ inline void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
+ inline void SetZProfile (double& ZProfile) {fZProfile=ZProfile;}
+ inline void SetZEmission (double& ZEmission) {fZEmission=ZEmission;}
+
+ // Get
+ // Nucleus* GetNucleus () const {return fBeamNucleus;}
+ inline double GetEnergy () const {return fEnergy;}
+ inline double GetExcitationEnergy() const {return fExcitationEnergy;}
+ inline double GetSigmaEnergy () const {return fSigmaEnergy;}
+ inline double GetMeanX () const {return fMeanX;}
+ inline double GetMeanY () const {return fMeanY;}
+ inline double GetSigmaX () const {return fSigmaX;}
+ inline double GetSigmaY () const {return fSigmaY;}
+ inline double GetMeanThetaX () const {return fMeanThetaX;}
+ inline double GetMeanPhiY () const {return fMeanPhiY;}
+ inline double GetSigmaThetaX () const {return fSigmaThetaX;}
+ inline double GetSigmaPhiY () const {return fSigmaPhiY;}
+ inline TH1F* GetEnergyHist () const {return fEnergyHist;}
+ inline TH2F* GetXThetaXHist () const {return fXThetaXHist;}
+ inline TH2F* GetYPhiYHist () const {return fYPhiYHist;}
+ inline int GetVerboseLevel() const {return fVerboseLevel;}
+ inline double GetZProfile () const {return fZProfile;}
+ inline double GetZEmission () const {return fZEmission;}
+
+
+ private: // Event Generation private variable
+ double fTargetSize;
+ double fEffectiveTargetSize; // target size has seen from the beam axis
+ double fTargetThickness;
+ double fEffectiveTargetThickness; // target thickness has seen by the beam
+ double fTargetAngle;
+ double fTargetZ;
+
+ public:
+ void SetTargetSize(double TargetSize);
+ void SetTargetThickness(double TargetThickness);
+ void SetTargetAngle(double TargetAngle);
+ void SetTargetZ(double TargetZ) {fTargetZ = TargetZ;}
+ double GetTargetSize(){return fTargetSize;}
+ double GetTargetThickness(){return fTargetThickness;}
+ double GetTargetAngle(){return fTargetAngle;}
+ double GetTargetZ() {return fTargetZ;}
+ double GetTargetEffectiveThickness(){return fEffectiveTargetThickness;}
+ double GetTargetEffectiveTargetSize(){return fEffectiveTargetSize;}
+
+
+ public: // Event Generation
+ void GenerateRandomEvent(double& E, double& X, double& Y, double& Z, double& ThetaX, double& PhiY );
+ public: // Print private paremeter
+ void Print() const;
};
}
#endif
diff --git a/NPLib/Physics/NPDecay.cxx b/NPLib/Physics/NPDecay.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..68a6eb6dbed984cecc7421e6ed593c3898fa5644
--- /dev/null
+++ b/NPLib/Physics/NPDecay.cxx
@@ -0,0 +1,438 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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 : Adrien Matta contact: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This Class hold data for all decay scheme of a given nuclei *
+ * *
+ * *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+#include "NPDecay.h"
+#include <iostream>
+#include "NPOptionManager.h"
+#include "NPFunction.h"
+#include "NPCore.h"
+#include "TF1.h"
+#include "TRandom.h"
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////// Single Decay //////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+NPL::SingleDecay::SingleDecay(std::string mother, double threshold, std::vector<std::string> daughter,std::vector<double> Ex, TH1F* CrossSection){
+
+ m_MotherName = mother;
+ m_DaughterName = daughter;
+ // Checking if there is a mixt of Gamma and particle
+ bool hasGamma=false;
+ bool hasOther=false;
+ unsigned int sizeD = daughter.size();
+ for(unsigned int d = 0 ; d < sizeD ; d++){
+ if(daughter[d]=="Gamma"|| daughter[d]=="gamma")
+ hasGamma=true;
+ else if (daughter[d]!=m_MotherName)
+ hasOther=true;
+ }
+
+ // If the same channel contain both, issue an error and exit
+ if(hasOther&&hasGamma){
+ NPL::SendErrorAndExit("NPL::SingleDecay","Cannot have particle/gamma mixt decay channel");
+ }
+ else if (hasGamma)
+ m_GammaOnly=true;
+ else
+ m_GammaOnly=false;
+
+ m_ExcitationEnergies = Ex;
+ m_CrossSectionHist = CrossSection;
+ m_Mother = NPL::Nucleus(m_MotherName);
+ m_MotherMass = m_Mother.Mass();
+ m_Threshold = threshold;
+ unsigned int size = m_DaughterName.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ m_Daughter.push_back(NPL::Nucleus(m_DaughterName[i])) ;
+
+ m_Daughter[i].SetExcitationEnergy(Ex[i]);
+ if (m_DaughterName[i]=="Gamma" || m_DaughterName[i]=="gamma")
+ m_DaughterMasses.push_back(0);
+ else
+ m_DaughterMasses.push_back(m_Daughter[i].Mass()/GeV);
+
+ }
+
+}
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::SingleDecay::GenerateEvent(double MEx, double MEK, double MPX, double MPY, double MPZ, std::vector<double>& DEK, std::vector<double>& DPx, std::vector<double>& DPy, std::vector<double>& DPz){
+ // Clear the output object
+ DPx.clear();
+ DPy.clear();
+ DPz.clear();
+ DEK.clear();
+ if(MEx < m_Threshold)
+ return false;
+
+ double Effective_Mass= m_MotherMass+MEx;
+ double NucleiEnergy= MEK+Effective_Mass;
+
+ double NucleiMomentum=
+ sqrt(NucleiEnergy*NucleiEnergy
+ - Effective_Mass*Effective_Mass);
+ TVector3 Momentum(MPX,MPY,MPZ);
+ Momentum = Momentum.Unit();
+
+ if(m_CrossSectionHist){
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.X(),
+ NucleiMomentum*Momentum.Y(),
+ NucleiMomentum*Momentum.Z(),
+ NucleiEnergy);
+ // Shoot the angle in Center of Mass (CM) frame
+ double ThetaCM = m_CrossSectionHist->GetRandom()* deg;
+ double phi = gRandom->Uniform()*2.*pi;
+
+ // Build daughter particule CM LV
+ // Pre compute variable for the decay
+ double m1 = m_DaughterMasses[0]*GeV;
+ double m2 = m_DaughterMasses[1]*GeV;
+ if(Effective_Mass<(m1+m2))
+ return false;
+
+ else {
+ double Energy = ( 1./(2.*Effective_Mass) )*(Effective_Mass*Effective_Mass + m1*m1 - m2*m2);
+ double Momentum1 = sqrt(Energy*Energy - m1*m1);
+
+ TVector3 FirstDaughterMomentum = Momentum1 * TVector3( sin(ThetaCM) * cos(phi),
+ sin(ThetaCM) * sin(phi),
+ cos(ThetaCM));
+
+ TLorentzVector FirstDaughterLV(FirstDaughterMomentum,Energy);
+
+ FirstDaughterLV.Boost( NucleiLV.BoostVector() );
+ TLorentzVector SecondDaughterLV = NucleiLV - FirstDaughterLV;
+
+ DEK.push_back(FirstDaughterLV.E()-m1);
+ DPx.push_back(FirstDaughterLV.X());
+ DPy.push_back(FirstDaughterLV.Y());
+ DPz.push_back(FirstDaughterLV.Z());
+
+ DEK.push_back(SecondDaughterLV.E()-m2);
+ DPx.push_back(SecondDaughterLV.X());
+ DPy.push_back(SecondDaughterLV.Y());
+ DPz.push_back(SecondDaughterLV.Z());
+ return true;
+ }
+ }
+
+ // Case of a TGenPhaseSpace
+ else if (!m_GammaOnly){
+ // TGenPhaseSpace require to input Energy and Momentum in GeV
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.x()/GeV,
+ NucleiMomentum*Momentum.y()/GeV,
+ NucleiMomentum*Momentum.z()/GeV,
+ NucleiEnergy/GeV);
+
+ if( !m_TGenPhaseSpace.SetDecay(NucleiLV, m_Daughter.size(), &m_DaughterMasses[0]) )
+ return false;
+
+ else{
+ // FIXME : weight has to be returned for normalisation
+ double weight = m_TGenPhaseSpace.Generate();
+
+ TLorentzVector* DaughterLV ;
+ double KineticEnergy;
+
+ for (unsigned int i = 0 ; i < m_Daughter.size(); i++) {
+ DaughterLV = m_TGenPhaseSpace.GetDecay(i);
+ DEK.push_back(GeV*DaughterLV->E()-GeV*m_DaughterMasses[i]);
+ DPx.push_back(GeV*DaughterLV->X());
+ DPy.push_back(GeV*DaughterLV->Y());
+ DPz.push_back(GeV*DaughterLV->Z());
+ }
+ return true;
+ }
+ }
+
+ // Case of a Gamma Cascade
+ else{
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.X(),
+ NucleiMomentum*Momentum.Y(),
+ NucleiMomentum*Momentum.Z(),
+ NucleiEnergy);
+
+ TLorentzVector TotalGamma(0,0,0,0);
+ unsigned int pos;
+
+ unsigned int sizeM = m_DaughterMasses.size();
+ for (unsigned int i = 0; i < sizeM; i++) {
+ // One of the Gamma
+ if(m_DaughterName[i]=="gamma"|| m_DaughterName[i]=="Gamma"){
+ // Shoot flat in cos(theta) and Phi to have isotropic emission
+ double cos_theta = -1+2*gRandom->Uniform();
+ double theta = acos(cos_theta);
+ double phi = gRandom->Uniform()*2.*pi;
+ double gammaEnergy = m_ExcitationEnergies[i];
+ TLorentzVector GammaLV( gammaEnergy*cos(phi)*sin(theta),
+ gammaEnergy*sin(phi)*sin(theta),
+ gammaEnergy*cos(theta),
+ gammaEnergy);
+
+ GammaLV.Boost(NucleiLV.BoostVector());
+ TotalGamma+=GammaLV;
+ DEK.push_back(GammaLV.E());
+ DPx.push_back(GammaLV.X());
+ DPy.push_back(GammaLV.Y());
+ DPz.push_back(GammaLV.Z());
+
+ }
+ else{
+ pos = i;
+ DEK.push_back(0);
+ DPx.push_back(0);
+ DPy.push_back(0);
+ DPz.push_back(0);
+ }
+ }
+
+ // The original nuclei get what is left
+ NucleiLV -= TotalGamma;
+ DEK[pos]=NucleiLV.E()-m_MotherMass;
+ DPx[pos]=NucleiLV.X();
+ DPy[pos]=NucleiLV.Y();
+ DPz[pos]=NucleiLV.Z();
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////// Decay ///////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+NPL::Decay::Decay(std::string MotherName,std::string path){
+ ReadConfiguration(MotherName,path);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+NPL::Decay::Decay(std::string MotherName, NPL::InputParser parser){
+ ReadConfiguration(MotherName,parser);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPL::Decay::ReadConfiguration(std::string MotherName,std::string path){
+ NPL::InputParser parser(path);
+ ReadConfiguration(MotherName,parser);
+}
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::Decay::AnyAboveThreshold(double MEx){
+ unsigned int size = m_SingleDecay.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ if(MEx >= m_SingleDecay[i].GetThreshold())
+ return true;
+ }
+ return false;
+}
+////////////////////////////////////////////////////////////////////////////////
+void NPL::Decay::ReadConfiguration(std::string MotherName, NPL::InputParser parser){
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("Decay",MotherName);
+ m_BRTotal=0;
+ m_MotherName=MotherName;
+
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ std::cout << "//// " << blocks.size() << " decay path found for " << MotherName << std::endl;
+
+ std::vector<std::string> token =
+ {"Daughter","Threshold","ExcitationEnergy","Shoot","LifeTime","BranchingRatio"};
+ std::vector<std::string> CStoken =
+ {"DifferentialCrossSection"};
+
+ unsigned int size = blocks.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ if(blocks[i]->HasTokenList(token)){
+ if(NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// "<< MotherName << " Decay " << i+1 << endl;
+
+ std::vector<std::string> Daughter = blocks[i]->GetVectorString("Daughter");
+
+ double Threshold = blocks[i]->GetDouble("Threshold","MeV");
+ std::vector<double> Ex = blocks[i]->GetVectorDouble("ExcitationEnergy","MeV");
+ double BR = blocks[i]->GetDouble("BranchingRatio","void");
+ m_BranchingRatio.push_back(BR);
+ m_BRTotal+= BR;
+ double LT = blocks[i]->GetDouble("LifeTime","ns");
+ m_Shoot = blocks[i]->GetVectorInt("Shoot");
+
+ TH1F* h = NULL;
+ std::vector<std::string> cs ;
+
+
+ if(blocks[i]->HasTokenList(CStoken)){
+ if(Ex.size()!=2){
+ std::cout << "ERROR : Differential cross section is only used for two body decay" << std::endl;
+ exit(1);
+ }
+ cs = blocks[i]->GetVectorString("DifferentialCrossSection");
+ }
+
+ // Load the Cross section
+ if(cs.size()==2){
+ h = NPL::Read1DProfile(cs[0],cs[1]);
+ TF1* sinus = new TF1("sinus","1/sin(x*3.141592653589793/180.)",0,180);
+ h->Divide(sinus,1);
+ delete sinus;
+ }
+ else if (cs.size()!=0){
+ std::cout << "ERROR : You should provide 2 arguments for the cross section : path and name" << std::endl;
+ exit(1);
+ }
+
+ m_SingleDecay.push_back(SingleDecay(m_MotherName,Threshold,Daughter,Ex,h));
+ }
+
+ else{
+ cout << "ERROR: check your input file formatting" << endl;
+ exit(1);
+ }
+
+ }
+ // Shift the Branching ratio for faster shooting during event generation
+ for (unsigned int i = 1; i < m_BranchingRatio.size(); i++) {
+ m_BranchingRatio[i] = m_BranchingRatio[i-1]+m_BranchingRatio[i];
+ }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::Decay::GenerateEvent(double MEx,double MEK,double MPx,double MPy,double MPz,
+ std::vector<NPL::Nucleus>& Daughter, std::vector<double>& Ex,
+ std::vector<double>& DEK,
+ std::vector<double>& DPx,std::vector<double>& DPy,std::vector<double>& DPz){
+ // Choose one of the Decay
+ // if the decay is not allowed, then try again
+ // FIXME: does take into account case where the decay is never possible
+ bool worked = false;
+ // Limit the number of attempt
+ unsigned int count =0;
+ while (!worked){
+ if(count++ > 1000)
+ break;
+ double rand = m_BRTotal*gRandom->Uniform();
+ unsigned int size = m_BranchingRatio.size();
+ unsigned int ChoosenDecay = 0;
+ for (unsigned int i = 1; i < size; i++) {
+ if(rand > m_BranchingRatio[i-1] && rand< m_BranchingRatio[i])
+ ChoosenDecay = i;
+ }
+
+ // Generate the event
+ worked = m_SingleDecay[ChoosenDecay].GenerateEvent(MEx,MEK,MPx,MPy,MPz,DEK,DPx,DPy,DPz);
+ if(worked){
+ Daughter = m_SingleDecay[ChoosenDecay].GetDaughter();
+ Ex = m_SingleDecay[ChoosenDecay].GetExcitationEnergies();
+ }
+ }
+
+ return worked;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Decay Store ////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+NPL::DecayStore::DecayStore(std::string path){
+ ReadConfiguration(path);
+}
+////////////////////////////////////////////////////////////////////////////////
+NPL::DecayStore::DecayStore(NPL::InputParser parser){
+ ReadConfiguration(parser);
+}
+////////////////////////////////////////////////////////////////////////////////
+void NPL::DecayStore::ReadConfiguration(std::string path){
+ NPL::InputParser parser(path);
+ ReadConfiguration(parser);
+}
+////////////////////////////////////////////////////////////////////////////////
+void NPL::DecayStore::ReadConfiguration(NPL::InputParser parser){
+ // Get all the Nucleus concern by a decay
+ std::vector<std::string> nuclei = parser.GetAllBlocksValues("Decay");
+ m_Store.clear();
+ unsigned int size = nuclei.size();
+ for(unsigned int i = 0 ; i < size ; i++){
+ if(m_Store.find(nuclei[i])==m_Store.end()){
+ m_Store[nuclei[i]] = Decay(nuclei[i],parser);
+ }
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void NPL::DecayStore::GenerateEvent(std::string MotherName, double MEx,
+ double MEK,double MPx,double MPy,double MPz,
+ std::vector<NPL::Nucleus>& Daughter, std::vector<double>& Ex,
+ std::vector<double>& DEK,
+ std::vector<double>& DPx,std::vector<double>& DPy,std::vector<double>& DPz){
+
+ if(m_Store.find(MotherName)!=m_Store.end())
+ m_Store[MotherName].GenerateEvent(MEx,MEK,MPx,MPy,MPz,Daughter,Ex,DEK,DPx,DPy,DPz);
+ else
+ cout << "Warning: Decay for " << MotherName << " requested but not found" << endl;
+}
+////////////////////////////////////////////////////////////////////////////////
+bool NPL::DecayStore::AnyAboveThreshold(std::string MotherName, double MEx){
+ if(m_Store.find(MotherName)!=m_Store.end())
+ m_Store[MotherName].AnyAboveThreshold(MEx);
+ else
+ return false;
+}
+////////////////////////////////////////////////////////////////////////////////
+std::set<std::string> NPL::DecayStore::GetAllMotherName(){
+ std::map<std::string,NPL::Decay>::iterator it;
+ std::set<std::string> result;
+ for(it = m_Store.begin(); it != m_Store.end() ; it++)
+ result.insert(it->first);
+
+ return result;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/NPLib/Physics/NPDecay.h b/NPLib/Physics/NPDecay.h
new file mode 100644
index 0000000000000000000000000000000000000000..9478a22d94d48c4eaa299468a2855e174e1145df
--- /dev/null
+++ b/NPLib/Physics/NPDecay.h
@@ -0,0 +1,133 @@
+#ifndef NPDECAY_H
+#define NPDECAY_H
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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 : Adrien Matta contact: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This Class hold data for all decay scheme of a given nuclei *
+ * *
+ * *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ * *
+ * *
+ *****************************************************************************/
+// C++ header
+#include <string>
+#include <vector>
+#include <set>
+
+// NPL Header
+#include "NPInputParser.h"
+#include "NPNucleus.h"
+
+// Root
+#include "TH1D.h"
+#include "TGenPhaseSpace.h"
+namespace NPL{
+ // A given decay Path
+ class SingleDecay{
+ public:
+ SingleDecay(){};
+ SingleDecay(std::string mother,double threshold, std::vector<std::string> daughter,
+ std::vector<double> Ex, TH1F* CrossSection);
+ ~SingleDecay(){};
+
+ private:
+ std::string m_MotherName;
+ NPL::Nucleus m_Mother;
+ double m_MotherMass;
+ std::vector<std::string> m_DaughterName;
+ std::vector<NPL::Nucleus> m_Daughter;
+ bool m_GammaOnly;
+ std::vector<double> m_DaughterMasses;
+ std::vector<double> m_ExcitationEnergies;
+ double m_Threshold;
+ TH1F* m_CrossSectionHist;
+ TGenPhaseSpace m_TGenPhaseSpace;
+
+ public:
+ // Given Energy and Momentum direction of Mother,
+ // Send back Momemtum and Energy of Daugther
+ // Return false if the decay is not allowed
+ // true other wise
+ bool GenerateEvent(double MEx,double MEK,double MPx, double MPy,double MPz,
+ std::vector<double>& DEK,
+ std::vector<double>& DPx,std::vector<double>& DPy,std::vector<double>& DPz);
+
+ public:// Getter
+ inline std::vector<std::string> GetDaughterName() {return m_DaughterName;};
+ inline std::vector<double> GetExcitationEnergies() {return m_ExcitationEnergies;};
+ inline std::vector<NPL::Nucleus> GetDaughter() {return m_Daughter;};
+ inline NPL::Nucleus GetMother() {return m_Mother;};
+ inline double GetThreshold(){return m_Threshold;};
+ };
+
+ // All decay Path for a given nuclei, with Branching Ratio
+ class Decay{
+ public:
+ Decay(){};
+ Decay(std::string MotherName,std::string path);
+ Decay(std::string MotherName,NPL::InputParser parser);
+ ~Decay(){};
+
+ public:
+ void ReadConfiguration(std::string MotherName, std::string path);
+ void ReadConfiguration(std::string MotherName, NPL::InputParser parser);
+ bool GenerateEvent(double MEx,double MEK,double MPx,double MPy,double MPz,
+ std::vector<NPL::Nucleus>& Daughter, std::vector<double>& Ex,
+ std::vector<double>& DEK,
+ std::vector<double>& DPx,std::vector<double>& DPy,std::vector<double>& DPz);
+ bool AnyAboveThreshold(double MEx);
+
+
+ private:
+ std::string m_MotherName;
+ std::vector<int> m_Shoot;
+ std::vector<NPL::SingleDecay> m_SingleDecay;
+ std::vector<double> m_BranchingRatio;
+ double m_BRTotal;
+ };
+
+ // All decay read from files
+ class DecayStore{
+ public:
+ DecayStore(){};
+ DecayStore(std::string path);
+ DecayStore(NPL::InputParser parser);
+ ~DecayStore(){};
+
+ public:
+ void ReadConfiguration(std::string path);
+ void ReadConfiguration(NPL::InputParser parser);
+ std::set<std::string> GetAllMotherName();
+
+ private:
+ std::map<std::string, NPL::Decay> m_Store;
+
+ public:
+ void GenerateEvent(std::string MotherName, double MEx,
+ double MEK,double MPx,double MPy,double MPz,
+ std::vector<NPL::Nucleus>& Daughter, std::vector<double>& Ex,
+ std::vector<double>& DEK,
+ std::vector<double>& DPx,std::vector<double>& DPy,std::vector<double>& DPz);
+ bool AnyAboveThreshold(std::string MotherName, double MEx);
+ };
+ }
+#endif
diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 71602f164667e633ce7966235034c47924788d89..c48465391db83de68e6fc3ab39ee7da70204b038 100644
--- a/NPLib/Physics/NPFunction.cxx
+++ b/NPLib/Physics/NPFunction.cxx
@@ -28,409 +28,394 @@ using namespace NPL;
#include "TDirectory.h"
namespace NPL{
-// Check the type of Filename (root or ASCII) and extract build/extract a 1D histogramm
-TH1F* Read1DProfile(string filename,string HistName)
-{
- ifstream ASCII;
- TFile ROOT;
- TH1F* h;
-
- // test whether file format is ASCII or ROOT
- bool type = OpenASCIIorROOTFile(filename, ASCII , ROOT);
-
- // ASCII File case
- if (type) {
- string LineBuffer;
-
- // storing vector
- vector <double> x, w;
-
- // variable buffer
- double xb, wb;
-
- // read the file
- Double_t xmin = 200;
- Double_t xmax = -200;
- int mysize = 0;
- while (getline(ASCII, LineBuffer)) {
- stringstream iss(LineBuffer);
- if (!(iss >> xb >> wb)) {continue;} // skip comment lines
- // fill vectors
- x.push_back(xb);
- w.push_back(wb);
-// cout << xb << "\t" << wb << endl;
- // compute xmin / xmax / size of x array
- if (xb > xmax) xmax = xb;
- if (xb < xmin) xmin = xb;
- mysize++;
- }
- Double_t dx = (xmax - xmin) / (mysize);
-// cout << xmin << "\t" << xmax << "\t" << mysize << "\t" << dx << endl;
-
- // fill histo
- h = new TH1F(HistName.c_str(), HistName.c_str(), mysize, xmin, xmax+dx);
- for (unsigned int i = 0; i < mysize; i++) {
- int bin = h->FindBin(x[i]);
- h->SetBinContent(bin,w[i]);
-// cout << i << "\t" << x[i] << "\t" << bin << "\t" << w[i] << "\t" << h->GetBinContent(bin) << endl;
- }
- }
-
- // ROOT File case
- else{
- h = (TH1F*) gDirectory->FindObjectAny(HistName.c_str());
- if(!h){
- cout << "Error: Histogram " << HistName << " not found in file " << filename << endl;
- exit(1);
- }
- }
-
- return h;
-}
+ // Check the type of Filename (root or ASCII) and extract build/extract a 1D histogramm
+ TH1F* Read1DProfile(string filename,string HistName) {
+ ifstream ASCII;
+ TFile ROOT;
+ TH1F* h;
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Check the type of Filename (root or ASCII) and extract build/extract a 2D histogramm
-TH2F* Read2DProfile(string filename,string HistName){
-
- ifstream ASCII;
- TFile ROOT;
- TH2F* h;
-
- bool type = OpenASCIIorROOTFile(filename, ASCII , ROOT);
-
-
- // ASCII File case
- if(type){
- string LineBuffer;
-
- // storing vector
- vector <double> x,y,w;
-
- // variable buffer
- double xb,yb,wb;
-
- // Read the file
- while(!ASCII.eof()){
- getline(ASCII,LineBuffer);
- stringstream LineStream(LineBuffer);
- // ignore comment line
- if (LineBuffer.compare(0,1,"%")!=0){
- LineStream >> xb >> yb >> wb ;
+ // test whether file format is ASCII or ROOT
+ bool type = OpenASCIIorROOTFile(filename, ASCII , ROOT);
+
+ // ASCII File case
+ if (type) {
+ string LineBuffer;
+
+ // storing vector
+ vector <double> x, w;
+
+ // variable buffer
+ double xb, wb;
+
+ // read the file
+ Double_t xmin = 200;
+ Double_t xmax = -200;
+ int mysize = 0;
+ while (getline(ASCII, LineBuffer)) {
+ stringstream iss(LineBuffer);
+ if (!(iss >> xb >> wb)) {continue;} // skip comment lines
+ // fill vectors
x.push_back(xb);
- y.push_back(yb);
w.push_back(wb);
+ // cout << xb << "\t" << wb << endl;
+ // compute xmin / xmax / size of x array
+ if (xb > xmax) xmax = xb;
+ if (xb < xmin) xmin = xb;
+ mysize++;
+ }
+ Double_t dx = (xmax - xmin) / (mysize);
+ // cout << xmin << "\t" << xmax << "\t" << mysize << "\t" << dx << endl;
+
+ // fill histo
+ h = new TH1F(HistName.c_str(), HistName.c_str(), mysize, xmin, xmax+dx);
+ for (unsigned int i = 0; i < mysize; i++) {
+ int bin = h->FindBin(x[i]);
+ h->SetBinContent(bin,w[i]);
+ // cout << i << "\t" << x[i] << "\t" << bin << "\t" << w[i] << "\t" << h->GetBinContent(bin) << endl;
}
}
-
- // Look for the step size, min and max of the distribution
- double xmin = 0;
- double xmax = 0;
- unsigned int xsize = x.size();
-
- double ymin = 0;
- double ymax = 0;
- unsigned int ysize = y.size();
-
- if(xsize > 0){
- xmin = x[0] ;
- xmax = x[0] ;
- }
-
- if(ysize > 0){
- ymin = y[0] ;
- ymax = y[0] ;
- }
-
- for(unsigned int i = 0 ; i < xsize ; i++){
- if(x[i] > xmax) xmax = x[i] ;
- if(x[i] < xmin) xmin = x[i] ;
- }
-
- for(unsigned int i = 0 ; i < ysize ; i++){
- if(y[i] > ymax) ymax = y[i] ;
- if(y[i] < ymin) ymin = y[i] ;
- }
-
- h = new TH2F(HistName.c_str(),HistName.c_str(),xsize,xmin,xmax,ysize,ymin,ymax);
-
- for(unsigned int i = 0 ; i < xsize ; i++){
- int bin = h->FindBin(x[i],y[i]);
- h->SetBinContent(bin,w[i]);
+
+ // ROOT File case
+ else{
+ h = (TH1F*) gDirectory->FindObjectAny(HistName.c_str());
+ if(!h){
+ cout << "Error: Histogram " << HistName << " not found in file " << filename << endl;
+ exit(1);
+ }
}
+
+ return h;
}
-
- // ROOT File case
- else{
- h = (TH2F*) gDirectory->FindObjectAny(HistName.c_str());
- if(!h){
- cout << "Error: Histogramm " << HistName << " not found in file " << filename << endl;
- exit(1);
+
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ // Check the type of Filename (root or ASCII) and extract build/extract a 2D histogramm
+ TH2F* Read2DProfile(string filename,string HistName){
+
+ ifstream ASCII;
+ TFile ROOT;
+ TH2F* h;
+
+ bool type = OpenASCIIorROOTFile(filename, ASCII , ROOT);
+
+
+ // ASCII File case
+ if(type){
+ string LineBuffer;
+
+ // storing vector
+ vector <double> x,y,w;
+
+ // variable buffer
+ double xb,yb,wb;
+
+ // Read the file
+ while(!ASCII.eof()){
+ getline(ASCII,LineBuffer);
+ stringstream LineStream(LineBuffer);
+ // ignore comment line
+ if (LineBuffer.compare(0,1,"%")!=0){
+ LineStream >> xb >> yb >> wb ;
+ x.push_back(xb);
+ y.push_back(yb);
+ w.push_back(wb);
+ }
+ }
+
+ // Look for the step size, min and max of the distribution
+ double xmin = 0;
+ double xmax = 0;
+ unsigned int xsize = x.size();
+
+ double ymin = 0;
+ double ymax = 0;
+ unsigned int ysize = y.size();
+
+ if(xsize > 0){
+ xmin = x[0] ;
+ xmax = x[0] ;
+ }
+
+ if(ysize > 0){
+ ymin = y[0] ;
+ ymax = y[0] ;
+ }
+
+ for(unsigned int i = 0 ; i < xsize ; i++){
+ if(x[i] > xmax) xmax = x[i] ;
+ if(x[i] < xmin) xmin = x[i] ;
+ }
+
+ for(unsigned int i = 0 ; i < ysize ; i++){
+ if(y[i] > ymax) ymax = y[i] ;
+ if(y[i] < ymin) ymin = y[i] ;
+ }
+
+ h = new TH2F(HistName.c_str(),HistName.c_str(),xsize,xmin,xmax,ysize,ymin,ymax);
+
+ for(unsigned int i = 0 ; i < xsize ; i++){
+ int bin = h->FindBin(x[i],y[i]);
+ h->SetBinContent(bin,w[i]);
+ }
}
- }
- return h;
-}
-// Open a file at Filename after checking the type of file it is
-// true for a ASCII file
-// False for a Root file
-bool OpenASCIIorROOTFile(string filename, ifstream &ASCII , TFile &ROOT){
-
- // look for .root extension
- size_t pos = filename.find(".root");
-
- string GlobalPath = getenv("NPTOOL");
- string StandardPath = GlobalPath + "/Inputs/CrossSection/" + filename;
-
- // extension not found, file is assume to be a ASCII file
- if(pos > filename.size()) {
- ASCII.open(filename.c_str());
-
- if(!ASCII.is_open()){
- ASCII.open(StandardPath.c_str());
- if(!ASCII.is_open()){
- cout << "Error, file " << filename << " not found " << endl ;
+ // ROOT File case
+ else{
+ h = (TH2F*) gDirectory->FindObjectAny(HistName.c_str());
+ if(!h){
+ cout << "Error: Histogramm " << HistName << " not found in file " << filename << endl;
exit(1);
}
}
-
- return true;
+ return h;
}
-
- else {
- if(ROOT.Open(filename.c_str(),"READ")){
- return false;
+
+ // Open a file at Filename after checking the type of file it is
+ // true for a ASCII file
+ // False for a Root file
+ bool OpenASCIIorROOTFile(string filename, ifstream &ASCII , TFile &ROOT){
+
+ // look for .root extension
+ size_t pos = filename.find(".root");
+
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/CrossSection/" + filename;
+
+ // extension not found, file is assume to be a ASCII file
+ if(pos > filename.size()) {
+ ASCII.open(filename.c_str());
+
+ if(!ASCII.is_open()){
+ ASCII.open(StandardPath.c_str());
+ if(!ASCII.is_open()){
+ cout << "Error, file " << filename << " not found " << endl ;
+ exit(1);
+ }
+ }
+
+ return true;
}
-
- else{
-
- if(ROOT.Open(StandardPath.c_str(),"READ")){
+
+ else {
+ if(ROOT.Open(filename.c_str(),"READ")){
return false;
}
+
else{
- cout << "Error, file " << StandardPath << " not found " << endl ;
- exit(1);
+
+ if(ROOT.Open(StandardPath.c_str(),"READ")){
+ return false;
+ }
+ else{
+ cout << "Error, file " << StandardPath << " not found " << endl ;
+ exit(1);
+ }
}
+
}
-
- }
-
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y, double NumberOfSigma)
-{
- if (SigmaX != 0 && SigmaY!=0) {
- X = 2 * NumberOfSigma*SigmaX;
- while (X > NumberOfSigma*SigmaX) X = gRandom->Gaus(MeanX, SigmaX);
-
- double a = NumberOfSigma * SigmaX/2;
- double b = NumberOfSigma * SigmaY/2;
- double SigmaYPrim = b * sqrt(1 - X*X / (a*a));
-
- SigmaYPrim = 2*SigmaYPrim / NumberOfSigma;
- Y = gRandom->Gaus(MeanY, SigmaYPrim);
- }
-
- else if(SigmaX == 0 && SigmaY!=0) {
- X = MeanX;
- Y = gRandom->Gaus(MeanY, SigmaY);
}
-
- else if(SigmaX != 0 && SigmaY==0) {
- Y = MeanY;
- X = gRandom->Gaus(MeanX, SigmaX);
- }
-
- else {
- X = MeanX;
- Y = MeanY;
+
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ void RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y){
+ if (SigmaX != 0 && SigmaY!=0) {
+ X = gRandom->Gaus(MeanX, SigmaX);
+
+ double a = SigmaX/2.;
+ double b = SigmaY/2.;
+ double SigmaYPrim = b * sqrt(1 - X*X / (a*a));
+ Y = gRandom->Gaus(MeanY, SigmaY);
+ }
+
+ else if(SigmaX == 0 && SigmaY!=0) {
+ X = MeanX;
+ Y = gRandom->Gaus(MeanY, SigmaY);
+ }
+
+ else if(SigmaX != 0 && SigmaY==0) {
+ Y = MeanY;
+ X = gRandom->Gaus(MeanX, SigmaX);
+ }
+
+ else {
+ X = MeanX;
+ Y = MeanY;
+ }
}
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-string ChangeNameToG4Standard(string OriginalName){
- string NumberOfMass ;
- string Nucleid;
-
- for (unsigned int i = 0; i < OriginalName.length(); i++) {
- ostringstream character;
- character << OriginalName[i];
- if (character.str()=="0") NumberOfMass+="0";
- else if (character.str()=="1") NumberOfMass+="1";
- else if (character.str()=="2") NumberOfMass+="2";
- else if (character.str()=="3") NumberOfMass+="3";
- else if (character.str()=="4") NumberOfMass+="4";
- else if (character.str()=="5") NumberOfMass+="5";
- else if (character.str()=="6") NumberOfMass+="6";
- else if (character.str()=="7") NumberOfMass+="7";
- else if (character.str()=="8") NumberOfMass+="8";
- else if (character.str()=="9") NumberOfMass+="9";
-
- else if (character.str()=="A") Nucleid+="A";
- else if (character.str()=="B") Nucleid+="B";
- else if (character.str()=="C") Nucleid+="C";
- else if (character.str()=="D") Nucleid+="D";
- else if (character.str()=="E") Nucleid+="E";
- else if (character.str()=="F") Nucleid+="F";
- else if (character.str()=="G") Nucleid+="G";
- else if (character.str()=="H") Nucleid+="H";
- else if (character.str()=="I") Nucleid+="I";
- else if (character.str()=="J") Nucleid+="J";
- else if (character.str()=="K") Nucleid+="K";
- else if (character.str()=="L") Nucleid+="L";
- else if (character.str()=="M") Nucleid+="M";
- else if (character.str()=="N") Nucleid+="N";
- else if (character.str()=="O") Nucleid+="O";
- else if (character.str()=="P") Nucleid+="P";
- else if (character.str()=="Q") Nucleid+="Q";
- else if (character.str()=="R") Nucleid+="R";
- else if (character.str()=="S") Nucleid+="S";
- else if (character.str()=="T") Nucleid+="T";
- else if (character.str()=="U") Nucleid+="U";
- else if (character.str()=="V") Nucleid+="V";
- else if (character.str()=="W") Nucleid+="W";
- else if (character.str()=="X") Nucleid+="X";
- else if (character.str()=="Y") Nucleid+="Y";
- else if (character.str()=="Z") Nucleid+="Z";
-
- else if (character.str()=="a") Nucleid+="a";
- else if (character.str()=="b") Nucleid+="b";
- else if (character.str()=="c") Nucleid+="c";
- else if (character.str()=="d") Nucleid+="d";
- else if (character.str()=="e") Nucleid+="e";
- else if (character.str()=="f") Nucleid+="f";
- else if (character.str()=="g") Nucleid+="g";
- else if (character.str()=="h") Nucleid+="h";
- else if (character.str()=="i") Nucleid+="i";
- else if (character.str()=="j") Nucleid+="j";
- else if (character.str()=="k") Nucleid+="k";
- else if (character.str()=="l") Nucleid+="l";
- else if (character.str()=="m") Nucleid+="m";
- else if (character.str()=="n") Nucleid+="n";
- else if (character.str()=="o") Nucleid+="o";
- else if (character.str()=="p") Nucleid+="p";
- else if (character.str()=="q") Nucleid+="q";
- else if (character.str()=="r") Nucleid+="r";
- else if (character.str()=="s") Nucleid+="s";
- else if (character.str()=="t") Nucleid+="t";
- else if (character.str()=="u") Nucleid+="u";
- else if (character.str()=="v") Nucleid+="v";
- else if (character.str()=="w") Nucleid+="w";
- else if (character.str()=="x") Nucleid+="x";
- else if (character.str()=="y") Nucleid+="y";
- else if (character.str()=="z") Nucleid+="z";
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ string ChangeNameToG4Standard(string OriginalName){
+ string NumberOfMass ;
+ string Nucleid;
+
+ for (unsigned int i = 0; i < OriginalName.length(); i++) {
+ ostringstream character;
+ character << OriginalName[i];
+ if (character.str()=="0") NumberOfMass+="0";
+ else if (character.str()=="1") NumberOfMass+="1";
+ else if (character.str()=="2") NumberOfMass+="2";
+ else if (character.str()=="3") NumberOfMass+="3";
+ else if (character.str()=="4") NumberOfMass+="4";
+ else if (character.str()=="5") NumberOfMass+="5";
+ else if (character.str()=="6") NumberOfMass+="6";
+ else if (character.str()=="7") NumberOfMass+="7";
+ else if (character.str()=="8") NumberOfMass+="8";
+ else if (character.str()=="9") NumberOfMass+="9";
+
+ else if (character.str()=="A") Nucleid+="A";
+ else if (character.str()=="B") Nucleid+="B";
+ else if (character.str()=="C") Nucleid+="C";
+ else if (character.str()=="D") Nucleid+="D";
+ else if (character.str()=="E") Nucleid+="E";
+ else if (character.str()=="F") Nucleid+="F";
+ else if (character.str()=="G") Nucleid+="G";
+ else if (character.str()=="H") Nucleid+="H";
+ else if (character.str()=="I") Nucleid+="I";
+ else if (character.str()=="J") Nucleid+="J";
+ else if (character.str()=="K") Nucleid+="K";
+ else if (character.str()=="L") Nucleid+="L";
+ else if (character.str()=="M") Nucleid+="M";
+ else if (character.str()=="N") Nucleid+="N";
+ else if (character.str()=="O") Nucleid+="O";
+ else if (character.str()=="P") Nucleid+="P";
+ else if (character.str()=="Q") Nucleid+="Q";
+ else if (character.str()=="R") Nucleid+="R";
+ else if (character.str()=="S") Nucleid+="S";
+ else if (character.str()=="T") Nucleid+="T";
+ else if (character.str()=="U") Nucleid+="U";
+ else if (character.str()=="V") Nucleid+="V";
+ else if (character.str()=="W") Nucleid+="W";
+ else if (character.str()=="X") Nucleid+="X";
+ else if (character.str()=="Y") Nucleid+="Y";
+ else if (character.str()=="Z") Nucleid+="Z";
+
+ else if (character.str()=="a") Nucleid+="a";
+ else if (character.str()=="b") Nucleid+="b";
+ else if (character.str()=="c") Nucleid+="c";
+ else if (character.str()=="d") Nucleid+="d";
+ else if (character.str()=="e") Nucleid+="e";
+ else if (character.str()=="f") Nucleid+="f";
+ else if (character.str()=="g") Nucleid+="g";
+ else if (character.str()=="h") Nucleid+="h";
+ else if (character.str()=="i") Nucleid+="i";
+ else if (character.str()=="j") Nucleid+="j";
+ else if (character.str()=="k") Nucleid+="k";
+ else if (character.str()=="l") Nucleid+="l";
+ else if (character.str()=="m") Nucleid+="m";
+ else if (character.str()=="n") Nucleid+="n";
+ else if (character.str()=="o") Nucleid+="o";
+ else if (character.str()=="p") Nucleid+="p";
+ else if (character.str()=="q") Nucleid+="q";
+ else if (character.str()=="r") Nucleid+="r";
+ else if (character.str()=="s") Nucleid+="s";
+ else if (character.str()=="t") Nucleid+="t";
+ else if (character.str()=="u") Nucleid+="u";
+ else if (character.str()=="v") Nucleid+="v";
+ else if (character.str()=="w") Nucleid+="w";
+ else if (character.str()=="x") Nucleid+="x";
+ else if (character.str()=="y") Nucleid+="y";
+ else if (character.str()=="z") Nucleid+="z";
+ }
+
+ // Special case for light particles
+ string FinalName=Nucleid+NumberOfMass;
+ if (FinalName=="H1") FinalName="proton";
+ else if (FinalName=="H2") FinalName="deuteron";
+ else if (FinalName=="H3") FinalName="triton";
+ else if (FinalName=="He4") FinalName="alpha";
+ else if (FinalName=="p") FinalName="proton";
+ else if (FinalName=="d") FinalName="deuteron";
+ else if (FinalName=="t") FinalName="triton";
+ else if (FinalName=="a") FinalName="alpha";
+ else if (FinalName=="n") FinalName="neutron";
+ return(FinalName);
}
-
- // Special case for light particles
- string FinalName=Nucleid+NumberOfMass;
- if (FinalName=="H1") FinalName="proton";
- else if (FinalName=="H2") FinalName="deuteron";
- else if (FinalName=="H3") FinalName="triton";
- else if (FinalName=="He4") FinalName="alpha";
- else if (FinalName=="p") FinalName="proton";
- else if (FinalName=="d") FinalName="deuteron";
- else if (FinalName=="t") FinalName="triton";
- else if (FinalName=="a") FinalName="alpha";
- else if (FinalName=="proton") FinalName="proton";
- else if (FinalName=="deuteron") FinalName="deuteron";
- else if (FinalName=="triton") FinalName="triton";
- else if (FinalName=="alpha") FinalName="alpha";
- else if (FinalName=="n") FinalName="neutron";
- else if (FinalName=="neutron") FinalName="neutron";
- return(FinalName);
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-string ChangeNameFromG4Standard(string OriginalName){
- string NumberOfMass ;
- string Nucleid;
-
- for (unsigned int i = 0; i < OriginalName.length(); i++) {
- ostringstream character;
- character << OriginalName[i];
- if (character.str()=="0") NumberOfMass+="0";
- else if (character.str()=="1") NumberOfMass+="1";
- else if (character.str()=="2") NumberOfMass+="2";
- else if (character.str()=="3") NumberOfMass+="3";
- else if (character.str()=="4") NumberOfMass+="4";
- else if (character.str()=="5") NumberOfMass+="5";
- else if (character.str()=="6") NumberOfMass+="6";
- else if (character.str()=="7") NumberOfMass+="7";
- else if (character.str()=="8") NumberOfMass+="8";
- else if (character.str()=="9") NumberOfMass+="9";
-
- else if (character.str()=="A") Nucleid+="A";
- else if (character.str()=="B") Nucleid+="B";
- else if (character.str()=="C") Nucleid+="C";
- else if (character.str()=="D") Nucleid+="D";
- else if (character.str()=="E") Nucleid+="E";
- else if (character.str()=="F") Nucleid+="F";
- else if (character.str()=="G") Nucleid+="G";
- else if (character.str()=="H") Nucleid+="H";
- else if (character.str()=="I") Nucleid+="I";
- else if (character.str()=="J") Nucleid+="J";
- else if (character.str()=="K") Nucleid+="K";
- else if (character.str()=="L") Nucleid+="L";
- else if (character.str()=="M") Nucleid+="M";
- else if (character.str()=="N") Nucleid+="N";
- else if (character.str()=="O") Nucleid+="O";
- else if (character.str()=="P") Nucleid+="P";
- else if (character.str()=="Q") Nucleid+="Q";
- else if (character.str()=="R") Nucleid+="R";
- else if (character.str()=="S") Nucleid+="S";
- else if (character.str()=="T") Nucleid+="T";
- else if (character.str()=="U") Nucleid+="U";
- else if (character.str()=="V") Nucleid+="V";
- else if (character.str()=="W") Nucleid+="W";
- else if (character.str()=="X") Nucleid+="X";
- else if (character.str()=="Y") Nucleid+="Y";
- else if (character.str()=="Z") Nucleid+="Z";
-
- else if (character.str()=="a") Nucleid+="a";
- else if (character.str()=="b") Nucleid+="b";
- else if (character.str()=="c") Nucleid+="c";
- else if (character.str()=="d") Nucleid+="d";
- else if (character.str()=="e") Nucleid+="e";
- else if (character.str()=="f") Nucleid+="f";
- else if (character.str()=="g") Nucleid+="g";
- else if (character.str()=="h") Nucleid+="h";
- else if (character.str()=="i") Nucleid+="i";
- else if (character.str()=="j") Nucleid+="j";
- else if (character.str()=="k") Nucleid+="k";
- else if (character.str()=="l") Nucleid+="l";
- else if (character.str()=="m") Nucleid+="m";
- else if (character.str()=="n") Nucleid+="n";
- else if (character.str()=="o") Nucleid+="o";
- else if (character.str()=="p") Nucleid+="p";
- else if (character.str()=="q") Nucleid+="q";
- else if (character.str()=="r") Nucleid+="r";
- else if (character.str()=="s") Nucleid+="s";
- else if (character.str()=="t") Nucleid+="t";
- else if (character.str()=="u") Nucleid+="u";
- else if (character.str()=="v") Nucleid+="v";
- else if (character.str()=="w") Nucleid+="w";
- else if (character.str()=="x") Nucleid+="x";
- else if (character.str()=="y") Nucleid+="y";
- else if (character.str()=="z") Nucleid+="z";
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ string ChangeNameFromG4Standard(string OriginalName){
+ string NumberOfMass ;
+ string Nucleid;
+
+ for (unsigned int i = 0; i < OriginalName.length(); i++) {
+ ostringstream character;
+ character << OriginalName[i];
+ if (character.str()=="0") NumberOfMass+="0";
+ else if (character.str()=="1") NumberOfMass+="1";
+ else if (character.str()=="2") NumberOfMass+="2";
+ else if (character.str()=="3") NumberOfMass+="3";
+ else if (character.str()=="4") NumberOfMass+="4";
+ else if (character.str()=="5") NumberOfMass+="5";
+ else if (character.str()=="6") NumberOfMass+="6";
+ else if (character.str()=="7") NumberOfMass+="7";
+ else if (character.str()=="8") NumberOfMass+="8";
+ else if (character.str()=="9") NumberOfMass+="9";
+
+ else if (character.str()=="A") Nucleid+="A";
+ else if (character.str()=="B") Nucleid+="B";
+ else if (character.str()=="C") Nucleid+="C";
+ else if (character.str()=="D") Nucleid+="D";
+ else if (character.str()=="E") Nucleid+="E";
+ else if (character.str()=="F") Nucleid+="F";
+ else if (character.str()=="G") Nucleid+="G";
+ else if (character.str()=="H") Nucleid+="H";
+ else if (character.str()=="I") Nucleid+="I";
+ else if (character.str()=="J") Nucleid+="J";
+ else if (character.str()=="K") Nucleid+="K";
+ else if (character.str()=="L") Nucleid+="L";
+ else if (character.str()=="M") Nucleid+="M";
+ else if (character.str()=="N") Nucleid+="N";
+ else if (character.str()=="O") Nucleid+="O";
+ else if (character.str()=="P") Nucleid+="P";
+ else if (character.str()=="Q") Nucleid+="Q";
+ else if (character.str()=="R") Nucleid+="R";
+ else if (character.str()=="S") Nucleid+="S";
+ else if (character.str()=="T") Nucleid+="T";
+ else if (character.str()=="U") Nucleid+="U";
+ else if (character.str()=="V") Nucleid+="V";
+ else if (character.str()=="W") Nucleid+="W";
+ else if (character.str()=="X") Nucleid+="X";
+ else if (character.str()=="Y") Nucleid+="Y";
+ else if (character.str()=="Z") Nucleid+="Z";
+
+ else if (character.str()=="a") Nucleid+="a";
+ else if (character.str()=="b") Nucleid+="b";
+ else if (character.str()=="c") Nucleid+="c";
+ else if (character.str()=="d") Nucleid+="d";
+ else if (character.str()=="e") Nucleid+="e";
+ else if (character.str()=="f") Nucleid+="f";
+ else if (character.str()=="g") Nucleid+="g";
+ else if (character.str()=="h") Nucleid+="h";
+ else if (character.str()=="i") Nucleid+="i";
+ else if (character.str()=="j") Nucleid+="j";
+ else if (character.str()=="k") Nucleid+="k";
+ else if (character.str()=="l") Nucleid+="l";
+ else if (character.str()=="m") Nucleid+="m";
+ else if (character.str()=="n") Nucleid+="n";
+ else if (character.str()=="o") Nucleid+="o";
+ else if (character.str()=="p") Nucleid+="p";
+ else if (character.str()=="q") Nucleid+="q";
+ else if (character.str()=="r") Nucleid+="r";
+ else if (character.str()=="s") Nucleid+="s";
+ else if (character.str()=="t") Nucleid+="t";
+ else if (character.str()=="u") Nucleid+="u";
+ else if (character.str()=="v") Nucleid+="v";
+ else if (character.str()=="w") Nucleid+="w";
+ else if (character.str()=="x") Nucleid+="x";
+ else if (character.str()=="y") Nucleid+="y";
+ else if (character.str()=="z") Nucleid+="z";
+ }
+
+ // Special case for light particles
+ string FinalName=NumberOfMass+Nucleid;
+ if (FinalName=="H1") FinalName="proton";
+ else if (FinalName=="H2") FinalName="deuteron";
+ else if (FinalName=="H3") FinalName="triton";
+ else if (FinalName=="He4") FinalName="alpha";
+ else if (FinalName=="p") FinalName="proton";
+ else if (FinalName=="d") FinalName="deuteron";
+ else if (FinalName=="t") FinalName="triton";
+ else if (FinalName=="a") FinalName="alpha";
+ else if (FinalName=="n") FinalName="neutron";
+ return(FinalName);
}
-
- // Special case for light particles
- string FinalName=NumberOfMass+Nucleid;
- if (FinalName=="H1") FinalName="proton";
- else if (FinalName=="H2") FinalName="deuteron";
- else if (FinalName=="H3") FinalName="triton";
- else if (FinalName=="He4") FinalName="alpha";
- else if (FinalName=="p") FinalName="proton";
- else if (FinalName=="d") FinalName="deuteron";
- else if (FinalName=="t") FinalName="triton";
- else if (FinalName=="a") FinalName="alpha";
- else if (FinalName=="proton") FinalName="proton";
- else if (FinalName=="deuteron") FinalName="deuteron";
- else if (FinalName=="triton") FinalName="triton";
- else if (FinalName=="alpha") FinalName="alpha";
- else if (FinalName=="n") FinalName="neutron";
- else if (FinalName=="neutron") FinalName="neutron";
- return(FinalName);
-}
}
diff --git a/NPLib/Physics/NPFunction.h b/NPLib/Physics/NPFunction.h
index a9fc43fde32def47b053120998a9a99173fc8f48..64858fd864032fdb9ac60ebaa5bc755def98e6bf 100644
--- a/NPLib/Physics/NPFunction.h
+++ b/NPLib/Physics/NPFunction.h
@@ -54,7 +54,7 @@ namespace NPL{
// False for a Root file
bool OpenASCIIorROOTFile(string filename, ifstream &ASCII , TFile &ROOT);
- void RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y, double NumberOfSigma= 10000);
+ void RandomGaussian2D(double MeanX, double MeanY, double SigmaX, double SigmaY, double &X, double &Y);
// Change nucleus name from G4 standard to Physics standard (11Li vs Li11)
string ChangeNameToG4Standard(string name);
diff --git a/NPLib/Physics/NPNucleus.cxx b/NPLib/Physics/NPNucleus.cxx
index 44038da618164558bd4df79f4ba6b84fbde418c5..4761b080bec6ad2c35ab2f50856c6b2e3969d42c 100644
--- a/NPLib/Physics/NPNucleus.cxx
+++ b/NPLib/Physics/NPNucleus.cxx
@@ -12,7 +12,7 @@
* Last update : may 2012 morfouac@ipno.in2p3.fr *
*---------------------------------------------------------------------------*
* Decription: *
- * This class manage a nucleus, data are read in the nubtab03.asc file *
+ * This class manage a nucleus, data are read in the nubtab12.asc file *
* *
*---------------------------------------------------------------------------*
* Comment: *
@@ -24,16 +24,7 @@
// NPTOOL headers
#include "NPNucleus.h"
using namespace NPL;
-//#include "Constant.h"
-// Use CLHEP System of unit and Physical Constant
-//#include "CLHEP/Units/GlobalSystemOfUnits.h"
-//#include "CLHEP/Units/PhysicalConstants.h"
-#include "NPGlobalSystemOfUnits.h"
-#include "NPPhysicalConstants.h"
-#ifdef NP_SYSTEM_OF_UNITS_H
-using namespace NPUNITS;
-#endif
// C++ headers
#include <iostream>
@@ -54,6 +45,7 @@ Nucleus::Nucleus(){
fCharge= 0;
fAtomicWeight= 0;
fMassExcess= 0;
+ fExcitationEnergy=0;
fSpinParity= "";
fSpin= 0;
fParity= "";
@@ -66,18 +58,23 @@ Nucleus::Nucleus(string isotope){
}
void Nucleus::SetUp(string isotope){
- //----------- Constructor Using nubtab03.asc by name----------
+ //----------- Constructor Using nubtab12.asc by name----------
// open the file to read and check if it is open
-
- // Replace the p,d,t,a by there standard name:
+ fExcitationEnergy=0;
+ // Replace the n,p,d,t,a by there standard name:
if(isotope=="p") isotope="1H";
else if(isotope=="d") isotope="2H";
else if(isotope=="t") isotope="3H";
else if(isotope=="a") isotope="4He";
+ else if(isotope=="n") isotope="1n";
+ else if(isotope=="neutron") isotope="1n";
+ else if(isotope=="g") isotope="gamma";
+ else if(isotope=="gamma") isotope="gamma";
+
ifstream inFile;
string Path = getenv("NPTOOL") ;
- string FileName = Path + "/NPLib/Physics/nubtab03.asc";
+ string FileName = Path + "/NPLib/Physics/nubtab12.asc";
inFile.open(FileName.c_str());
// reading the file
@@ -103,12 +100,12 @@ void Nucleus::SetUp(string isotope){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Nucleus::Nucleus(int Z, int A)
{
- //----------- Constructor Using nubtab03.asc by Z and A----------
+ //----------- Constructor Using nubtab12.asc by Z and A----------
// open the file to read and check if it is open
ifstream inFile;
string Path = getenv("NPTOOL") ;
- string FileName = Path + "/NPLib/Physics/nubtab03.asc";
+ string FileName = Path + "/NPLib/Physics/nubtab12.asc";
inFile.open(FileName.c_str());
// reading the file
@@ -209,9 +206,6 @@ void Nucleus::Extract(string line){
fSpinParity = s_spinparity.data();
size_t found_p = s_spinparity.find("+");
size_t found_m = s_spinparity.find("-");
- // size_t found_( = s_jpi.find("(");
- // size_t found_) = s_jpi.find(")");
- // size_t found_# = s_jpi.find("#");
// parity
if (found_p != string::npos) fParity = "+";
if (found_m != string::npos) fParity = "-";
diff --git a/NPLib/Physics/NPNucleus.h b/NPLib/Physics/NPNucleus.h
index 18bf35d0ef4bfc00ecf01047aa65309474899946..58373c56b48923109d4f57a26e8acbe14c2f0bc6 100644
--- a/NPLib/Physics/NPNucleus.h
+++ b/NPLib/Physics/NPNucleus.h
@@ -28,14 +28,20 @@
// NPTOOL headers
#include "NPGlobalSystemOfUnits.h"
#include "NPPhysicalConstants.h"
+
#ifdef NP_SYSTEM_OF_UNITS_H
using namespace NPUNITS;
#endif
+#ifdef HEP_PHYSICAL_CONSTANTS_H
+using namespace CLHEP;
+#endif
+
// C++ headers
#include <string>
using namespace std;
+#include <iostream>
namespace NPL {
class Nucleus {
@@ -50,16 +56,16 @@ namespace NPL {
private :
//intrinsic properties
- string fName; // Nucleus name
- string fNucleusName;
- int fCharge; // Nucleus charge
- int fAtomicWeight; // Nucleus atomic weight
- double fMassExcess; // Nucleus mass excess in keV
- string fSpinParity; // Nucleus spin and parity
- double fSpin; // Nucleus spin
- string fParity; // Nucleus parity
- double fLifeTime; // life time
-
+ string fName; // Nucleus name
+ string fNucleusName;
+ int fCharge; // Nucleus charge
+ int fAtomicWeight; // Nucleus atomic weight
+ double fMassExcess; // Nucleus mass excess in keV
+ string fSpinParity; // Nucleus spin and parity
+ double fSpin; // Nucleus spin
+ string fParity; // Nucleus parity
+ double fLifeTime; // life time
+ double fExcitationEnergy; // Excitation Energy
//kinematic properties
double fKineticEnergy;
double fBeta;
@@ -102,18 +108,19 @@ namespace NPL {
double GetBeta() const {return fBeta;}
double GetGamma() const {return fGamma;}
double GetVelocity() const {return fVelocity;}
- TLorentzVector GetEnergyImpulsion() const {return fEnergyImpulsion;}
- void SetName(const char* name) {fName = name;}
- void SetZ(int charge) {fCharge = charge;}
- void SetA(int mass) {fAtomicWeight = mass;}
- void SetMassExcess(double massexcess) {fMassExcess = massexcess;}
+ TLorentzVector GetEnergyImpulsion() const {return fEnergyImpulsion;}
+ double GetExcitationEnergy() const {return fExcitationEnergy;}
+ void SetName(const char* name) {fName = name;}
+ void SetZ(int charge) {fCharge = charge;}
+ void SetA(int mass) {fAtomicWeight = mass;}
+ void SetMassExcess(double massexcess) {fMassExcess = massexcess;}
void SetSpinParity(const char* spinparity) {fSpinParity = spinparity;}
- void SetSpin(double spin) {fSpin = spin;}
- void SetParity(const char* parity) {fParity = parity;}
+ void SetSpin(double spin) {fSpin = spin;}
+ void SetParity(const char* parity) {fParity = parity;}
void SetLifeTime(double LifeTime) {fLifeTime=LifeTime;}
- void SetKineticEnergy(double energy) {fKineticEnergy = energy; EnergyToBrho(); EnergyToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
- void SetBrho(double brho) {fBrho = brho; BrhoToEnergy(); BrhoToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
- void SetTimeOfFlight(double tof) {fTimeOfFlight = tof; TofToEnergy(); TofToBrho(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
+ void SetKineticEnergy(double energy) {fKineticEnergy = energy; EnergyToBrho(); EnergyToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
+ void SetBrho(double brho) {fBrho = brho; BrhoToEnergy(); BrhoToTof(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
+ void SetTimeOfFlight(double tof) {fTimeOfFlight = tof; TofToEnergy(); TofToBrho(); EnergyToBeta(); BetaToGamma();BetaToVelocity();}
void SetEnergyImpulsion(TLorentzVector P) {fEnergyImpulsion = P;
fKineticEnergy = P.E() - Mass();
EnergyToBrho();
@@ -121,12 +128,13 @@ namespace NPL {
EnergyToBeta();
BetaToGamma();
BetaToVelocity();}
+ void SetExcitationEnergy(double Ex) {fExcitationEnergy=Ex;}
void SetBeta(double beta) {fBeta = beta; BetaToGamma(); BetaToEnergy(); EnergyToTof(); EnergyToBrho();BetaToVelocity();}
double GetEnergyCM(double EnergyLab, double ThetaLab, double PhiLab, double relativisticboost);
double GetThetaCM(double EnergyLab, double ThetaLab, double PhiLab, double relativisticboost);
// Nuclear mass in MeV
- double Mass() const {return (fAtomicWeight*amu_c2 + fMassExcess/1000. - fCharge*electron_mass_c2);}
+ double Mass() const {return (fAtomicWeight*amu_c2 + fMassExcess/1000. - fCharge*electron_mass_c2+fExcitationEnergy);}
double GetBindingEnergy() const {return (fCharge*proton_mass_c2 + (fAtomicWeight-fCharge)*neutron_mass_c2 + fCharge*electron_mass_c2 - fAtomicWeight*amu_c2 - fMassExcess/1000);}
void Print() const ;
};
diff --git a/NPLib/Physics/NPPhysicsLinkDef.h b/NPLib/Physics/NPPhysicsLinkDef.h
index 7f7bb77816c1424f53d2fb4cd94b851cbd8b9975..4b4ace982dd68c055180e2d0a7e3b23a36145313 100644
--- a/NPLib/Physics/NPPhysicsLinkDef.h
+++ b/NPLib/Physics/NPPhysicsLinkDef.h
@@ -3,4 +3,5 @@
#pragma link C++ defined_in "./NPBeam.h";
#pragma link C++ defined_in "./NPNucleus.h";
#pragma link C++ defined_in "./NPEnergyLoss.h";
+#pragma link C++ defined_in "./NPDecay.h";
#endif
diff --git a/NPLib/Physics/nubtab12.asc b/NPLib/Physics/nubtab12.asc
index 86e7023a03acef84f7cf1f02a056afcda3a6d4ae..b16be589f1a75ebbdda755400aa8487cc2c3cebe 100644
--- a/NPLib/Physics/nubtab12.asc
+++ b/NPLib/Physics/nubtab12.asc
@@ -1,3 +1,4 @@
+000 0000 gamma 0.0 0.0000 stbl 0+ 06 1932 IS=0.000000 00
001 0000 1n 8071.3171 0.0005 613.9 s 0.6 1/2+ 06 1932 B-=100
001 0010 1H 7288.9705 0.0001 stbl 1/2+ 06 11Be53d 1920 IS=99.9885 70
002 0010 2H 13135.7217 0.0001 stbl 1+ 03 1932 IS=0.0115 70
@@ -5508,4 +5509,4 @@
293 1180 293Ei 198930# 730# RN 1# ms 1/2+# 00 02Ni10i A ?
294 1170 294Eh 196040# 690# 290 ms 230 10 10Og01td 2010 A=100
294 1180 294Ei 199270# 660# 1.4 ms 0.7 0+ 05 06Og05td 2006 A=100
-295 1180 295Ei 201430# 640# 10# ms 04Og05td A ?
\ No newline at end of file
+295 1180 295Ei 201430# 640# 10# ms 04Og05td A ?
diff --git a/NPLib/scripts/Style_nponline.C b/NPLib/scripts/Style_nponline.C
index 8f5ac365611e969a108a9e62cd10cac9514b9b27..aeb7eb4a4debd4eefc56f4ac033c5693e3407732 100644
--- a/NPLib/scripts/Style_nponline.C
+++ b/NPLib/scripts/Style_nponline.C
@@ -81,7 +81,7 @@ style->SetFrameFillColor(kGray+3);
style->SetTitlePS("nptool");
- const UInt_t Number = 2;
+/* const UInt_t Number = 2;
Double_t Red[Number] = { 0,0 };
Double_t Green[Number] = { 0,0.8 };
Double_t Blue[Number] = { 0,1.00 };
@@ -89,5 +89,7 @@ style->SetFrameFillColor(kGray+3);
Double_t Length[Number] = { 0,1.00 };
Int_t nb=255;
TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
- style->SetNumberContours(99);
+*/
+ style->SetNumberContours(99);
+ style->SetPalette();
}
diff --git a/NPLib/scripts/Style_nptool.C b/NPLib/scripts/Style_nptool.C
index 6026c490c028391d99680b2427621c72944dc2ed..60c8f49926f3e4e2d0831939d74a16a88076b15b 100644
--- a/NPLib/scripts/Style_nptool.C
+++ b/NPLib/scripts/Style_nptool.C
@@ -80,13 +80,15 @@ void Style_nptool(){
style->SetFuncWidth(2);
// Create the color gradiant
- const UInt_t Number = 2;
- Double_t Red[Number] = { 0,0 };
- Double_t Green[Number] = { 0,0.8 };
- Double_t Blue[Number] = { 0,1.00 };
+// const UInt_t Number = 2;
+// Double_t Red[Number] = { 0,0 };
+// Double_t Green[Number] = { 0,0.8 };
+// Double_t Blue[Number] = { 0,1.00 };
- Double_t Length[Number] = { 0,1.00 };
- Int_t nb=255;
- TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
+// Double_t Length[Number] = { 0,1.00 };
+// Int_t nb=255;
+// TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
style->SetNumberContours(99);
+ style->SetPalette(0);
+
}
diff --git a/NPSimulation/Core/DetectorConstruction.cc b/NPSimulation/Core/DetectorConstruction.cc
index a5f9cb379cd9dcf47e9fda11e2becad57f88ec85..16b31980cd7b6cf68490c14574f809808b258827 100644
--- a/NPSimulation/Core/DetectorConstruction.cc
+++ b/NPSimulation/Core/DetectorConstruction.cc
@@ -56,11 +56,10 @@
#include "NPSDetectorFactory.hh"
#include "MaterialManager.hh"
#include "DetectorMessenger.hh"
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
DetectorConstruction::DetectorConstruction(): world_log(0), world_phys(0){
- m_Target = 0;
- m_Chamber = 0 ;
+ m_Target = NULL;
+ m_Chamber = NULL ;
m_Messenger = new DetectorMessenger(this);
m_ReadSensitivePtr = &NPS::VDetector::ReadSensitive;
}
@@ -74,8 +73,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct(){
return ReadConfigurationFile();
}
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void DetectorConstruction::AddDetector(NPS::VDetector* NewDetector){
// Add new detector to vector
@@ -92,7 +89,6 @@ void DetectorConstruction::AddDetector(NPS::VDetector* NewDetector){
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4VPhysicalVolume* DetectorConstruction::ReadConfigurationFile(){
@@ -110,7 +106,7 @@ G4VPhysicalVolume* DetectorConstruction::ReadConfigurationFile(){
world_log = new G4LogicalVolume(world_box, Vacuum, "world_log", 0, 0, 0);
world_phys = new G4PVPlacement(0, G4ThreeVector(), world_log, "world", 0, false, 0);
- G4VisAttributes* VisAtt = new G4VisAttributes(G4VisAttributes::Invisible);
+ G4VisAttributes* VisAtt = new G4VisAttributes(G4VisAttributes::Invisible);
world_log->SetVisAttributes(VisAtt);
//------------------------------------------------------------------
diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index 236eab375ff5d25d7c3c03017806295546578e9f..c8d6641f9c7a01c458b11adbac240899c8dea4d7 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -106,18 +106,19 @@ void ParticleStack::AddBeamParticleToStack(Particle& particle){
m_ParticleStack.push_back(particle);
// Incident beam parameter
m_InitialConditions-> SetIncidentParticleName (particle.GetParticleDefinition()->GetParticleName());
- m_InitialConditions-> SetIncidentFinalKineticEnergy (particle. GetParticleKineticEnergy());
+ //m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleThetaCM());
+ m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleKineticEnergy());
G4ThreeVector U(1,0,0);
G4ThreeVector V(0,1,0);
m_InitialConditions-> SetIncidentEmittanceThetaX (particle.GetParticleMomentumDirection().angle(U)/deg);
m_InitialConditions-> SetIncidentEmittancePhiY (particle.GetParticleMomentumDirection().angle(V)/deg);
- m_InitialConditions-> SetIncidentEmittanceTheta (particle.GetParticleMomentumDirection().theta()/deg);
- m_InitialConditions-> SetIncidentEmittancePhi (particle.GetParticleMomentumDirection().phi()/deg);
+ m_InitialConditions-> SetIncidentEmittanceTheta (particle.GetParticleMomentumDirection().theta()/deg);
+ m_InitialConditions-> SetIncidentEmittancePhi (particle.GetParticleMomentumDirection().phi()/deg);
// Beam status at the initial interaction point
- m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleThetaCM());
+ m_InitialConditions-> SetIncidentFinalKineticEnergy (particle. GetParticleKineticEnergy());
m_InitialConditions-> SetIncidentPositionX (particle. GetParticlePosition().x());
m_InitialConditions-> SetIncidentPositionY (particle. GetParticlePosition().y());
m_InitialConditions-> SetIncidentPositionZ (particle. GetParticlePosition().z());
diff --git a/NPSimulation/Core/PrimaryGeneratorAction.cc b/NPSimulation/Core/PrimaryGeneratorAction.cc
index 361a95fe196c7878b9b5bd53098263659423f761..4c761b3c19d9231c4edb6e5bfe57f828385323e9 100644
--- a/NPSimulation/Core/PrimaryGeneratorAction.cc
+++ b/NPSimulation/Core/PrimaryGeneratorAction.cc
@@ -51,7 +51,8 @@
PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* det): m_detector(det){
m_Messenger = new PrimaryGeneratorActionMessenger(this);
m_GenerateEvent = &NPS::VEventGenerator::GenerateEvent;
-}
+ m_Target=NULL;
+ }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PrimaryGeneratorAction::~PrimaryGeneratorAction(){
@@ -111,7 +112,7 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
myEventGenerator->SetTarget(m_detector->GetTarget());
m_EventGenerator.push_back(myEventGenerator);
}
- blocks.clear();
+/* blocks.clear();
blocks = parser.GetAllBlocksWithToken("TwoBodyReaction");
if (blocks.size()>0) {
NPS::VEventGenerator* myEventGenerator = new EventGeneratorTwoBodyReaction();
@@ -138,6 +139,10 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
myEventGenerator->SetTarget(m_detector->GetTarget());
m_EventGenerator.push_back(myEventGenerator);
}
+*/
+ m_Target=m_detector->GetTarget();
+ m_Target->SetReactionRegion();
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::ClearEventGenerator(){
@@ -152,7 +157,7 @@ void PrimaryGeneratorAction::ClearEventGenerator(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::SetTarget(){
for (unsigned int i = 0 ; i < m_EventGenerator.size(); i++) {
- m_EventGenerator[i]->SetTarget(m_detector->GetTarget());
+ m_EventGenerator[i]->SetTarget(m_Target);
}
}
diff --git a/NPSimulation/Core/PrimaryGeneratorAction.hh b/NPSimulation/Core/PrimaryGeneratorAction.hh
index 7a52d1b5a36468447c4367770e1a3f19b4361069..3a09784f4d160335f149712cf8e4d5b0b9201462 100644
--- a/NPSimulation/Core/PrimaryGeneratorAction.hh
+++ b/NPSimulation/Core/PrimaryGeneratorAction.hh
@@ -62,7 +62,7 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction{
DetectorConstruction* m_detector;
vector<NPS::VEventGenerator*> m_EventGenerator;
PrimaryGeneratorActionMessenger* m_Messenger;
-
+ Target* m_Target;
};
#endif
diff --git a/NPSimulation/Core/PrimaryGeneratorActionMessenger.hh b/NPSimulation/Core/PrimaryGeneratorActionMessenger.hh
index 37c9d314f582f0ae34cf9e6746b60c900095d513..0b9e6af09b15118d854ad7de361659260fa1c8a8 100644
--- a/NPSimulation/Core/PrimaryGeneratorActionMessenger.hh
+++ b/NPSimulation/Core/PrimaryGeneratorActionMessenger.hh
@@ -14,7 +14,7 @@
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
- * This class describe the PrimaryGeneratorAction Messenger *
+ * This class describe the PrimaryGeneratorAction Messenger *
* *
*---------------------------------------------------------------------------*
* Comment: *
diff --git a/NPSimulation/Core/Target.cc b/NPSimulation/Core/Target.cc
index ecb81aab77764c8ae2203f2f1c2c99c420777d45..535abe37a3b3bb62a464aae86f3b0922a735f898 100644
--- a/NPSimulation/Core/Target.cc
+++ b/NPSimulation/Core/Target.cc
@@ -45,13 +45,16 @@
#include "G4EmCalculator.hh"
#include "G4ParticleDefinition.hh"
#include "G4ParticleTable.hh"
+#include "G4UserLimits.hh"
#include "Randomize.hh"
+#include "BeamReaction.hh"
+#include "G4FastSimulationManager.hh"
using namespace CLHEP ;
// NPS
-#include"Target.hh"
-#include"MaterialManager.hh"
-
+#include "Target.hh"
+#include "MaterialManager.hh"
+#include "Decay.hh"
// NPL
#include "NPOptionManager.h"
#include "NPInputParser.h"
@@ -71,6 +74,7 @@ Target::Target(){
m_TargetNbLayers = 5; // Number of steps by default
// Set the global pointer
TargetInstance = this;
+ m_ReactionRegion=0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -122,7 +126,7 @@ void Target::ReadConfiguration(NPL::InputParser parser){
m_WindowsThickness= ctarget[0]->GetDouble("WindowsThickness","micrometer");
m_WindowsMaterial= GetMaterialFromLibrary(ctarget[0]->GetString("WindowsMaterial"));
m_TargetRadius=ctarget[0]->GetDouble("Radius","mm");
-
+
m_TargetX=ctarget[0]->GetDouble("X","mm");
m_TargetY=ctarget[0]->GetDouble("Y","mm");
m_TargetZ =ctarget[0]->GetDouble("Z","mm");
@@ -131,7 +135,7 @@ void Target::ReadConfiguration(NPL::InputParser parser){
cout << "ERROR: Target token list incomplete, check your input file" << endl;
exit(1);
}
-
+
if(ctarget[0]->HasToken("NBLAYERS"))
m_TargetNbLayers = ctarget[0]->GetInt("NBLAYERS");
@@ -175,14 +179,14 @@ void Target::ConstructDetector(G4LogicalVolume* world){
new G4Tubs("solidTarget", 0, m_TargetRadius,
0.5*m_TargetThickness+m_WindowsThickness, 0*deg, 360*deg);
- m_TargetLogic =
+ m_TargetLogic =
new G4LogicalVolume(m_TargetSolid,
- GetMaterialFromLibrary("Vacuum")
- , "logicTarget");
+ GetMaterialFromLibrary("Vacuum")
+ , "logicTarget");
m_TargetLogic->SetVisAttributes(G4VisAttributes::Invisible);
G4Tubs* solidTarget =
- new G4Tubs("solidTarget", 0, m_TargetRadius,
+ new G4Tubs("solidTarget", 0, m_TargetRadius,
0.5*m_TargetThickness, 0*deg, 360*deg);
G4LogicalVolume* logicTarget =
@@ -191,7 +195,7 @@ void Target::ConstructDetector(G4LogicalVolume* world){
new G4PVPlacement(0, G4ThreeVector(0, 0, 0),
logicTarget, "Target", m_TargetLogic, false, 0);
- new G4PVPlacement(0, G4ThreeVector(m_TargetX, m_TargetY, m_TargetZ),
+ new G4PVPlacement(0, G4ThreeVector(m_TargetX, m_TargetY, m_TargetZ),
m_TargetLogic, "Target", world, false, 0);
G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour(0., 0., 1.));
@@ -231,7 +235,30 @@ void Target::ConstructDetector(G4LogicalVolume* world){
}
}
+ //SetReactionRegion();
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Target::SetReactionRegion(){
+ if(!m_ReactionRegion){
+ m_ReactionRegion= new G4Region("NPSimulationProcess");
+ m_ReactionRegion->AddRootLogicalVolume(m_TargetLogic);
+ m_ReactionRegion->SetUserLimits(new G4UserLimits(m_TargetThickness/10.));
+ }
+
+ G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
+
+ unsigned int size = m_ReactionModel.size();
+ for(unsigned int i = 0 ; i < size ; i++)
+ mng->RemoveFastSimulationModel(m_ReactionModel[i]);
+
+ m_ReactionModel.clear();
+ G4VFastSimulationModel* fsm;
+ fsm = new NPS::BeamReaction("BeamReaction",m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
+ fsm = new NPS::Decay("Decay",m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
+
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
diff --git a/NPSimulation/Core/Target.hh b/NPSimulation/Core/Target.hh
index 833243887047c596dc17dc4e0c70219bf8826f97..d2b04f9895414d386236b799a22565b5ff462f65 100644
--- a/NPSimulation/Core/Target.hh
+++ b/NPSimulation/Core/Target.hh
@@ -38,7 +38,7 @@
#include "G4Material.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
-
+#include "G4VFastSimulationModel.hh"
// NPTool headers
#include "NPSVDetector.hh"
#include "NPInputParser.h"
@@ -121,6 +121,14 @@ private:
G4double m_TargetY;
G4double m_TargetZ;
+
+public: // Region were reaction can occure
+ void SetReactionRegion();
+
+private:
+ // Region were reaction can occure:
+ G4Region* m_ReactionRegion;
+ vector<G4VFastSimulationModel*> m_ReactionModel;
// Pointer to the last target generated (0 if none)
private:
static Target* TargetInstance ;
diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc
index a02d0cdee388aa7bf2cb7180891552ee83eccdc0..d809e2960fdf4bca8a38a0b668a98dc87512422a 100644
--- a/NPSimulation/Detectors/Chio/Chio.cc
+++ b/NPSimulation/Detectors/Chio/Chio.cc
@@ -148,7 +148,9 @@ G4LogicalVolume* Chio::BuildDetector(){
G4Material* Fe= MaterialManager::getInstance()->GetMaterialFromLibrary("Fe");
G4Material* Al= MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin);
+ //G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("CF4",0.0693276*bar,273.15*kelvin);
+ G4Material* CF4= MaterialManager::getInstance()->GetGasFromLibrary("iC4H10",8.*bar/1000.,273.15*kelvin);
+
G4Material* Mylar= MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
diff --git a/NPSimulation/Detectors/TRex/TRex.cc b/NPSimulation/Detectors/TRex/TRex.cc
index 5908162cf760b6dcf4fb71b7462b19698ba0003b..f8c0376c0a0a0ffaae3b1d31804af06b5920b676 100644
--- a/NPSimulation/Detectors/TRex/TRex.cc
+++ b/NPSimulation/Detectors/TRex/TRex.cc
@@ -134,7 +134,7 @@ G4LogicalVolume* TRex::BuildChamber(){
cout << "TRex geometry is based on Munich Group Simulation exported in GDML"<< endl;
string basepath = getenv("NPTOOL");
string path=basepath+"/NPSimulation/Detectors/TRex/TRex_Miniball.gdml";
- m_gdmlparser.Read(path);
+ m_gdmlparser.Read(path,false);
}
m_Chamber= m_gdmlparser.GetVolume("chamber_log");
}
diff --git a/NPSimulation/EventGenerator/EventGeneratorBeam.cc b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
index c6b67dc1f74bc6c516f6ab2e6516f7d81168f5ea..748fa663713a2de716f2c4edf0685262abf71866 100644
--- a/NPSimulation/EventGenerator/EventGeneratorBeam.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
@@ -55,9 +55,8 @@ EventGeneratorBeam::~EventGeneratorBeam(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorBeam::SetTarget(Target* Target){
- if(Target!=0){
+ if(Target!=0)
m_Target = Target;
- }
// Set the target parameter for the internal event generator of m_Beam
m_Beam->SetTargetSize(m_Target->GetTargetRadius());
@@ -75,12 +74,12 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
if( anEvent->GetEventID()==0){
// Define the particle to be shoot
- if(m_Beam->GetZ()==0 && m_Beam->GetA()==1){
+ if(m_Beam->GetZ()==0 && m_Beam->GetA()==1)
m_particle = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
- }
else
- m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
+ m_particle =
+ G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
}
@@ -89,28 +88,23 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
///// of interaction in target and Energy Loss of the beam within /////
///// the target. /////
///////////////////////////////////////////////////////////////////////
- G4ThreeVector InterCoord;
+ //G4ThreeVector InterCoord;
- G4double Beam_theta, Beam_phi, FinalBeamEnergy, InitialBeamEnergy, x0, y0, z0, Beam_thetaX, Beam_phiY;
+ 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
- G4double Xdir = sin(Beam_thetaX); // cos(90-x) = sin(x)
- G4double Ydir = sin(Beam_phiY);
- G4double Zdir = sqrt(1-Xdir*Xdir-Ydir*Ydir); // alpha^2 + beta^2 + gamma^2 = 1
+ 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
G4ThreeVector BeamDir(Xdir,Ydir,Zdir);
G4ThreeVector BeamPos(x0,y0,z0);
- Beam_theta = BeamDir.theta() ;
- Beam_phi = BeamDir.phi() ; Beam_phi *= 1;
- FinalBeamEnergy = m_Target->SlowDownBeam(m_particle, InitialBeamEnergy,z0-m_Beam->GetTargetZ(),Beam_theta);
- if(FinalBeamEnergy<0 || FinalBeamEnergy!=FinalBeamEnergy)
- FinalBeamEnergy=0;
///////////////////////////////////////////////////////
///// Add the Beam particle to the particle Stack /////
///////////////////////////////////////////////////////
Particle BeamParticle( m_particle,
InitialBeamEnergy,
- FinalBeamEnergy,
+ InitialBeamEnergy,
BeamDir.unit(),
BeamPos,
1);
diff --git a/NPSimulation/EventGenerator/EventGeneratorGammaDecay.cc b/NPSimulation/EventGenerator/EventGeneratorGammaDecay.cc
index d7d7471ce6c8352707168fc2d483e16e6a9e7901..af815fab7e4fc5eded5ac4fa081b5e88b1e457fa 100644
--- a/NPSimulation/EventGenerator/EventGeneratorGammaDecay.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorGammaDecay.cc
@@ -89,7 +89,8 @@ void EventGeneratorGammaDecay::ReadConfiguration(NPL::InputParser parser){
else
AddCascade(E, BranchingRatio);
}
-PrepareCascade();
+
+ PrepareCascade();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a55b2b36d852d2ab867d8588cf249c742707ab6d
--- /dev/null
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -0,0 +1,232 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#include <iostream>
+#include <string>
+#include "NPFunction.h"
+#include "BeamReaction.hh"
+#include "NPOptionManager.h"
+#include "NPInputParser.h"
+#include "G4VPhysicalVolume.hh"
+#include "G4Electron.hh"
+#include "G4Gamma.hh"
+#include "G4SystemOfUnits.hh"
+////////////////////////////////////////////////////////////////////////////////
+NPS::BeamReaction::BeamReaction(G4String modelName,G4Region* envelope) :
+ G4VFastSimulationModel(modelName, envelope) {
+ ReadConfiguration();
+ m_PreviousEnergy=0 ;
+ m_PreviousLength=0 ;
+ }
+
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::BeamReaction::BeamReaction(G4String modelName) :
+ G4VFastSimulationModel(modelName) {
+ }
+
+////////////////////////////////////////////////////////////////////////////////
+NPS::BeamReaction::~BeamReaction() {
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::BeamReaction::ReadConfiguration(){
+ NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
+ m_Reaction.ReadConfigurationFile(input);
+ m_BeamName=NPL::ChangeNameToG4Standard(m_Reaction.GetNucleus1().GetName());
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool NPS::BeamReaction::IsApplicable( const G4ParticleDefinition& particleType) {
+ std::string particleName = particleType.GetParticleName();
+ if (particleName.find(m_BeamName)!=std::string::npos) {
+ return true;
+ }
+ return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
+ static bool shoot = false;
+ static double rand = 0;
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector V = PrimaryTrack->GetMomentum();
+ G4ThreeVector P = fastTrack.GetPrimaryTrackLocalPosition();
+ G4VSolid* solid =
+ fastTrack.GetPrimaryTrack()->GetVolume()->GetLogicalVolume()->GetSolid();
+ double in = solid->DistanceToOut(P,V);
+ double out = solid->DistanceToOut(P,-V);
+ double ratio = in / (out+in) ;
+ // Generate a random for this event
+ if(!shoot){
+ rand = G4RandFlat::shoot();
+ shoot = true;
+ }
+
+ // If the condition is met, the event is generated
+ // if PreviousLength is NULL then the particle has taken no step yet in the
+ // region, so reaction could not happen
+ if(ratio<rand && m_PreviousLength!=0){
+ // Reset the static for next event
+ shoot = false;
+ return true;
+ }
+
+ // Record the situation of the current step
+ // so it can be used in the next one
+ if(!PrimaryTrack->GetStep()->IsLastStepInVolume()){
+ m_PreviousLength = PrimaryTrack->GetStep()->GetStepLength();
+ m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
+ }
+
+ return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep) {
+ // Get the track info
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
+ G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
+
+ G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
+ G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
+
+ double energy = PrimaryTrack->GetKineticEnergy();
+ double time = PrimaryTrack->GetGlobalTime();
+
+ // Randomize within the step
+ // Assume energy loss is linear within the step
+ // Assume no scattering
+ double rand = G4RandFlat::shoot();
+ double length = rand*(m_PreviousLength);
+ energy += (1-rand)*(m_PreviousEnergy-energy);
+ G4ThreeVector ldir = pdirection;
+ ldir*=length;
+ localPosition = localPosition - ldir;
+ // Set the end of the step conditions
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0,pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetTotalEnergyDeposited(0);
+ fastStep.SetPrimaryTrackFinalTime (time);
+ fastStep.KillPrimaryTrack();
+ fastStep.SetPrimaryTrackPathLength(0.0);
+
+ //////////////////////////////////////////////////
+ //////Define the kind of particle to shoot////////
+ //////////////////////////////////////////////////
+
+ // Nucleus 3
+ int LightZ = m_Reaction.GetNucleus3().GetZ() ;
+ int LightA = m_Reaction.GetNucleus3().GetA() ;
+ static G4IonTable* IonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+ G4ParticleDefinition* LightName
+ = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3()*MeV);
+
+ // Nucleus 4
+ G4int HeavyZ = m_Reaction.GetNucleus4().GetZ() ;
+ G4int HeavyA = m_Reaction.GetNucleus4().GetA() ;
+
+ // Generate the excitation energy if a distribution is given
+ m_Reaction.ShootRandomExcitationEnergy();
+
+ // Use to clean up the IonTable in case of the Ex changing at every event
+ static G4ParticleDefinition* previousHeavy=0;
+ if(previousHeavy)
+ delete previousHeavy;
+ // IonTable->Remove(previousHeavy);
+
+ G4ParticleDefinition* HeavyName
+ = IonTable->GetIon(HeavyZ, HeavyA, m_Reaction.GetExcitation4()*MeV);
+ //cout << IonTable->Entries()<< endl;
+ // Set the Energy of the reaction
+ m_Reaction.SetBeamEnergy(energy);
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+ G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ cos(Beam_theta) * sin(Beam_phi),
+ -sin(Beam_theta));
+ G4ThreeVector col2(-sin(Beam_phi),
+ cos(Beam_phi),
+ 0);
+ G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ sin(Beam_theta) * sin(Beam_phi),
+ cos(Beam_theta));
+ G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+
+ // Angles
+ // Shoot and Set a Random ThetaCM
+ m_Reaction.ShootRandomThetaCM();
+ double phi = RandFlat::shoot() * 2. * pi;
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ double Theta3, Energy3, Theta4, Energy4;
+
+ // Compute Kinematic using previously defined ThetaCM
+ m_Reaction.KineRelativistic(Theta3, Energy3, Theta4, Energy4);
+
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kine3_beam(sin(Theta3) * cos(phi),
+ sin(Theta3) * sin(phi),
+ cos(Theta3));
+ // Momentum in World frame
+ G4ThreeVector momentum_kine3_world = BeamToWorld * momentum_kine3_beam;
+
+
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kine4_beam(sin(Theta4) * cos(phi+pi),
+ sin(Theta4) * sin(phi+pi),
+ cos(Theta4));
+ // Momentum in World frame
+ G4ThreeVector momentum_kine4_world = BeamToWorld * momentum_kine4_beam;
+
+
+ // Emitt secondary
+ G4DynamicParticle particle3(LightName,momentum_kine3_world,Energy3);
+ fastStep.CreateSecondaryTrack(particle3, localPosition, time);
+ G4DynamicParticle particle4(HeavyName,momentum_kine4_world,Energy4);
+ fastStep.CreateSecondaryTrack(particle4, localPosition, time);
+
+ // Reinit for next event
+ m_PreviousEnergy=0 ;
+ m_PreviousLength=0 ;
+ //previousHeavy=HeavyName;
+ //IonTable->fIonList->begin()->second->Mass();
+ //IonTable->fIonList->erase(IonTable->fIonList->begin());
+
+}
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3ec04679e5ed141163773b290d23f9bb67a5711a
--- /dev/null
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -0,0 +1,53 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#ifndef BeamReaction_h
+#define BeamReaction_h
+
+#include "G4VFastSimulationModel.hh"
+#include "PhysicsList.hh"
+#include "NPReaction.h"
+class G4VPhysicalVolume;
+namespace NPS{
+ class BeamReaction : public G4VFastSimulationModel{
+ public:
+ BeamReaction (G4String, G4Region*);
+ BeamReaction (G4String);
+ ~BeamReaction ();
+
+ public:
+ void ReadConfiguration();
+ virtual G4bool IsApplicable(const G4ParticleDefinition&);
+ virtual G4bool ModelTrigger(const G4FastTrack &);
+ virtual void DoIt(const G4FastTrack&, G4FastStep&);
+
+ private:
+ NPL::Reaction m_Reaction;
+ string m_BeamName;
+ double m_PreviousEnergy;
+ double m_PreviousLength;
+ };
+}
+
+
+#endif
diff --git a/NPSimulation/Process/CMakeLists.txt b/NPSimulation/Process/CMakeLists.txt
index c5e28e4e2f7c60cba71dbd956e9d1550ef8fcd9d..769c3a6c874cb5986645a8c90e9ae94e7376e5a8 100644
--- a/NPSimulation/Process/CMakeLists.txt
+++ b/NPSimulation/Process/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSProcess SHARED FastDriftElectron.cc NPIonIonInelasticPhysic.cc PhysicsList.cc G4DriftElectron.cc G4IonizationWithDE.cc G4DriftElectronPhysics.cc G4DEAbsorption.cc G4DEAmplification.cc G4DETransport.cc )
+add_library(NPSProcess SHARED Decay.cc BeamReaction.cc FastDriftElectron.cc NPIonIonInelasticPhysic.cc PhysicsList.cc G4DriftElectron.cc G4IonizationWithDE.cc G4DriftElectronPhysics.cc G4DEAbsorption.cc G4DEAmplification.cc G4DETransport.cc )
target_link_libraries(NPSProcess ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates)
diff --git a/NPSimulation/Process/Decay.cc b/NPSimulation/Process/Decay.cc
new file mode 100644
index 0000000000000000000000000000000000000000..10040f64cc6c3868900cdc3a22e78b0fb8aefa75
--- /dev/null
+++ b/NPSimulation/Process/Decay.cc
@@ -0,0 +1,173 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+ * *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#include <iostream>
+#include <string>
+#include <set>
+#include "NPFunction.h"
+#include "Decay.hh"
+#include "NPOptionManager.h"
+#include "NPInputParser.h"
+#include "G4VPhysicalVolume.hh"
+#include "G4Electron.hh"
+#include "G4Gamma.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
+
+
+using namespace NPS;
+////////////////////////////////////////////////////////////////////////////////
+Decay::Decay(G4String modelName,G4Region* envelope) :
+ G4VFastSimulationModel(modelName, envelope) {
+ ReadConfiguration();
+ m_PreviousEnergy=0 ;
+ m_PreviousLength=0 ;
+ }
+
+
+////////////////////////////////////////////////////////////////////////////////
+Decay::Decay(G4String modelName) :
+ G4VFastSimulationModel(modelName) {
+ }
+
+////////////////////////////////////////////////////////////////////////////////
+Decay::~Decay() {
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Decay::ReadConfiguration(){
+ NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
+ m_Decay.ReadConfiguration(input);
+ std::set<std::string> Mother = m_Decay.GetAllMotherName();
+ std::set<std::string>::iterator it ;
+ for(it = Mother.begin() ; it != Mother.end() ; it++)
+ m_MotherName.insert(NPL::ChangeNameToG4Standard(*it));
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool Decay::IsApplicable( const G4ParticleDefinition& particleType) {
+ m_CurrentName = particleType.GetParticleName();
+ // Extract Ex from name
+ if(m_CurrentName.find("[")!=std::string::npos)
+ m_ExcitationEnergy = atof(m_CurrentName.substr(m_CurrentName.find("[")+1,m_CurrentName.find("]")-1).c_str())*keV;
+ else
+ m_ExcitationEnergy=0;
+
+
+ // Strip name from excitation energy
+ m_CurrentName = m_CurrentName.substr(0,m_CurrentName.find("["));
+ // If the decay exist
+ if (m_MotherName.find(m_CurrentName)!=m_MotherName.end()) {
+ return true;
+ }
+ return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+G4bool Decay::ModelTrigger(const G4FastTrack& fastTrack) {
+ //FIXME: Solve the issue of long lived decay
+ // Check that a decay is possible:
+ return m_Decay.AnyAboveThreshold(NPL::ChangeNameFromG4Standard(m_CurrentName),m_ExcitationEnergy);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Decay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
+ // Get the track info
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
+ G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
+
+ G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
+ G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
+
+ double energy = PrimaryTrack->GetKineticEnergy();
+ double time = PrimaryTrack->GetGlobalTime();
+
+ // Randomize within the step
+ // Assume energy loss is linear within the step
+ // Assume no scattering
+ double rand = G4RandFlat::shoot();
+ double length = rand*(m_PreviousLength);
+ energy += (1-rand)*(m_PreviousEnergy-energy);
+ G4ThreeVector ldir = pdirection;
+ ldir*=length;
+ localPosition = localPosition - ldir;
+ //////////////////////////////////////////////////
+ //////Define the kind of particle to shoot////////
+ //////////////////////////////////////////////////
+ std::vector<NPL::Nucleus> Daughter;
+ std::vector<double> Ex;
+ std::vector<double> DEK;
+ std::vector<double> DPx;
+ std::vector<double> DPy;
+ std::vector<double> DPz;
+
+ m_Decay.GenerateEvent(NPL::ChangeNameFromG4Standard(m_CurrentName),m_ExcitationEnergy,energy,
+ pdirection.x(),pdirection.y(),pdirection.z(),
+ Daughter, Ex,DEK,DPx,DPy,DPz);
+
+
+ G4ParticleDefinition* DaughterDef;
+ unsigned int size = Daughter.size();
+ if(size == 0)
+ return;
+ for(unsigned int i = 0 ; i < size ; i++){
+ // Get the decaying particle
+ int DaughterZ = Daughter[i].GetZ();
+ int DaughterA = Daughter[i].GetA();
+
+ // neutral particle
+ if(DaughterZ==0){
+ if(DaughterA==1)
+ DaughterDef=G4ParticleTable::GetParticleTable()->FindParticle("neutron");
+
+ else if(DaughterA==0){
+ DaughterDef=G4ParticleTable::GetParticleTable()->FindParticle("gamma");
+ }
+
+ }
+ // proton
+ else if (DaughterZ==1 && DaughterA==1 )
+ DaughterDef=G4ParticleTable::GetParticleTable()->FindParticle("proton");
+ // the rest
+ else
+ DaughterDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(DaughterZ, DaughterA, Ex[i]);
+
+ // Set the momentum direction
+ G4ThreeVector Momentum (DPx[i],DPy[i],DPz[i]);
+ Momentum=Momentum.unit();
+
+ G4DynamicParticle DynamicDaughter(DaughterDef,Momentum,DEK[i]);
+ fastStep.CreateSecondaryTrack(DynamicDaughter, localPosition, time);
+ }
+ if(size){
+ // Set the end of the step conditions
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0,pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetTotalEnergyDeposited(0);
+ fastStep.SetPrimaryTrackFinalTime (time);
+ fastStep.KillPrimaryTrack();
+ fastStep.SetPrimaryTrackPathLength(0.0);
+ }
+}
diff --git a/NPSimulation/Process/Decay.hh b/NPSimulation/Process/Decay.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ae6190811e2b287151770ef2e3ed6a7fd7381203
--- /dev/null
+++ b/NPSimulation/Process/Decay.hh
@@ -0,0 +1,58 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * *
+ * Creation Date : Octobre 2017 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Use to kill the beam track and replace it with the reaction product *
+ * *
+* *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ *****************************************************************************/
+
+#ifndef Decay_h
+#define Decay_h
+
+#include "G4VFastSimulationModel.hh"
+#include "PhysicsList.hh"
+#include "NPDecay.h"
+class G4VPhysicalVolume;
+
+////////////////////////////////////////////////////////////////////////////////
+namespace NPS{
+
+class Decay : public G4VFastSimulationModel{
+ public:
+ Decay (G4String, G4Region*);
+ Decay (G4String);
+ ~Decay();
+
+ public:
+ void ReadConfiguration();
+ virtual G4bool IsApplicable(const G4ParticleDefinition&);
+ virtual G4bool ModelTrigger(const G4FastTrack &);
+ virtual void DoIt(const G4FastTrack&, G4FastStep&);
+
+ private:
+ NPL::DecayStore m_Decay;
+ std::set<std::string> m_MotherName;
+ std::string m_CurrentName;
+ double m_ExcitationEnergy;
+ double m_PreviousEnergy;
+ double m_PreviousLength;
+
+};
+}
+
+#endif
diff --git a/NPSimulation/Process/PhysicsList.cc b/NPSimulation/Process/PhysicsList.cc
index cf08ee4881d969ba626ae8c5bbe1c2e33d2e191d..29bf7065ab635184d02f4354f2fd39d89381ae23 100644
--- a/NPSimulation/Process/PhysicsList.cc
+++ b/NPSimulation/Process/PhysicsList.cc
@@ -27,7 +27,7 @@
// Local physic directly implemented from the Hadronthrapy example
// Physic dedicated to the ion-ion inelastic processes
#include "NPIonIonInelasticPhysic.hh"
-
+#include "Decay.hh"
// G4
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
@@ -37,6 +37,7 @@
#include "G4UnitsTable.hh"
#include "G4ProcessManager.hh"
#include "G4FastSimulationManagerProcess.hh"
+#include "G4StepLimiter.hh"
/////////////////////////////////////////////////////////////////////////////
PhysicsList::PhysicsList() : G4VUserPhysicsList(){
m_EmList = "Option4";
@@ -283,7 +284,7 @@ void PhysicsList::ConstructProcess(){
em_option.SetAuger(true);
AddParametrisation();
-
+
return;
}
/////////////////////////////////////////////////////////////////////////////
@@ -291,9 +292,8 @@ void PhysicsList::AddStepMax(){
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddParametrisation(){
-/*
- G4FastSimulationManagerProcess* drift =
- new G4FastSimulationManagerProcess("DriftElectron");
+ G4FastSimulationManagerProcess* BeamReaction =
+ new G4FastSimulationManagerProcess("NPSProcess");
// For 10.3 and higher
#ifndef theParticleIterator
@@ -304,10 +304,12 @@ void PhysicsList::AddParametrisation(){
while ((*theParticleIterator)()){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
-
- if(particle->GetParticleName()=="e-")
- pmanager->AddDiscreteProcess(drift);
- }*/
+ std::string name = particle->GetParticleName();
+ pmanager->AddDiscreteProcess(BeamReaction);
+ // Add a Step limiter to the beam particle.
+ // This will be used to limit the step of the beam in the target
+ pmanager->AddProcess(new G4StepLimiter,-1,-1,5);
+ }
}
diff --git a/NPSimulation/Process/PhysicsList.hh b/NPSimulation/Process/PhysicsList.hh
index c171bc6271a3f9eb9897f189f1f4c42b9edcb0fb..66ceaa685f2f212011a250413130a2cfd56076a1 100644
--- a/NPSimulation/Process/PhysicsList.hh
+++ b/NPSimulation/Process/PhysicsList.hh
@@ -84,7 +84,7 @@ class PhysicsList: public G4VUserPhysicsList{
void AddParametrisation();
void AddPackage(const G4String& name);
void BiasCrossSectionByFactor(double factor);
-
+
private:
std::map<std::string,G4VPhysicsConstructor*> m_PhysList;
G4EmConfigurator em_config;
@@ -96,7 +96,7 @@ class PhysicsList: public G4VUserPhysicsList{
G4VPhysicsConstructor* emPhysicsList;
G4VPhysicsConstructor* decay_List;
G4VPhysicsConstructor* radioactiveDecay_List;
-
+
private: // Physics option
std::string m_EmList;
double m_IonBinaryCascadePhysics;
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index df6c6886f7565858cd0f9720f3f284a441484d53..b15420bc908d392d02336750a1a3d417ed9d7c72 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -74,7 +74,7 @@ int main(int argc, char** argv){
PhysicsList* physicsList = new PhysicsList();
runManager->SetUserInitialization(physicsList);
PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector);
-
+
// Initialize Geant4 kernel
runManager->Initialize();
diff --git a/Projects/TRex_Miniball/Analysis.cxx b/Projects/TRex_Miniball/Analysis.cxx
index 8bd30237c82fed822b2dcc135d0458255f865263..fbe9c96f2dbf88d67870ca7500f1a29366f5fb95 100644
--- a/Projects/TRex_Miniball/Analysis.cxx
+++ b/Projects/TRex_Miniball/Analysis.cxx
@@ -42,12 +42,12 @@ void Analysis::Init(){
myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
string ligth,heavy;
- if(myReaction->GetNucleus3()->GetZ()==1 && myReaction->GetNucleus3()->GetA()==1)
+ if(myReaction->GetNucleus3().GetZ()==1 && myReaction->GetNucleus3().GetA()==1)
ligth = "proton";
else
ligth = "deuteron";
- heavy = "N17";
+ heavy = "Ga80";
Sharc = (TSharcPhysics*) m_DetectorManager -> GetDetector("Sharc");
diff --git a/Projects/TRex_Miniball/Reaction/17Ndp.reaction b/Projects/TRex_Miniball/Reaction/17Ndp.reaction
index f1fbd52ea4d0693e056d9e1c71651c5fad902e71..9287cda0c85d7ba896ecf69a7ce486a310004ddf 100644
--- a/Projects/TRex_Miniball/Reaction/17Ndp.reaction
+++ b/Projects/TRex_Miniball/Reaction/17Ndp.reaction
@@ -1,30 +1,30 @@
%%%%%%%%%%%%%%%%%%%%%% IS591 at ISOLDE %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
- Particle= 17N
- ExcitationEnergy= 0
- Energy= 93
- SigmaEnergy= 0.448
- SigmaThetaX= 0.01
- SigmaPhiY= 0.01
- SigmaX= 0.5
- SigmaY= 0.5
- MeanThetaX= 0
- MeanPhiY= 0
- MeanX= 0
- MeanY= 0
- %EnergyProfilePath=
- %XThetaXProfilePath=
- %YPhiYProfilePath=
+ Particle= 17N
+ ExcitationEnergy= 0
+ Energy= 93
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TwoBodyReaction
- Beam= 17N
- Target= 2H
- Light= 1H
- Heavy= 18N
- ExcitationEnergyLight= 0.0
- ExcitationEnergyHeavy= 0.0
- CrossSectionPath= 18Ngs.txt CSR
- ShootLight= 1
- ShootHeavy= 1
+ Beam= 17N
+ Target= 2H
+ Light= 1H
+ Heavy= 18N
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 0.0
+ CrossSectionPath= flat.txt CSR
+ ShootLight= 1
+ ShootHeavy= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/TRex_Miniball/TRex_Miniball.detector b/Projects/TRex_Miniball/TRex_Miniball.detector
index 210bb5c22c38f25706688592b59a375800706db3..62f238675ef5d32ee0793a91ba3d8eed73e9e2cd 100644
--- a/Projects/TRex_Miniball/TRex_Miniball.detector
+++ b/Projects/TRex_Miniball/TRex_Miniball.detector
@@ -1,135 +1,132 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-GeneralTarget
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
Target
- THICKNESS= 5
- RADIUS= 5
- MATERIAL= CD2
- ANGLE= 0
- X= 0
- Y= 0
- Z= 0
+ THICKNESS= 5 micrometer
+ RADIUS= 5 mm
+ MATERIAL= CD2
+ ANGLE= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- Chamber= 1
- X= 0
- Y= 29
- Z= 33
- Shape= Square
+ Chamber= 1
+ X= 0 mm
+ Y= 29 mm
+ Z= 33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= 0
- Y= -29
- Z= 33
- Shape= Square
+ X= 0 mm
+ Y= -29 mm
+ Z= 33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= 29
- Y= 0
- Z= 33
- Shape= Square
+ X= 29 mm
+ Y= 0 mm
+ Z= 33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= -29
- Y= 0
- Z= 33
- Shape= Square
+ X= -29 mm
+ Y= 0 mm
+ Z= 33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= 0
- Y= 29
- Z= -33
- Shape= Square
+ X= 0 mm
+ Y= 29 mm
+ Z= -33 mm
+Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= 0
- Y= -29
- Z= -33
- Shape= Square
+ X= 0 mm
+ Y= -29 mm
+ Z= -33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= 29
- Y= 0
- Z= -33
- Shape= Square
+ X= 29 mm
+ Y= 0 mm
+ Z= -33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TRex
- X= -29
- Y= 0
- Z= -33
- Shape= Square
+ X= -29 mm
+ Y= 0 mm
+ Z= -33 mm
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Sharc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %Upstream CD
- SharcQQQ
- Z= -70
- R= 0
- Phi= 2
- ThicknessDector= 500
- SharcQQQ
- Z= -70
- R= 0
- Phi= 92
- ThicknessDector= 500
- SharcQQQ
- Z= -70
- R= 0
- Phi= 182
- ThicknessDector= 500
- SharcQQQ
- Z= -70
- R= 0
- Phi= 272
- ThicknessDector= 500
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Miniball
- R= 110
- THETA= 110
- PHI= 320
- Shape= Square
+%Upstream CD
+Sharc QQQ
+ Z= -70 mm
+ R= 0 mm
+ Phi= 2 deg
+ ThicknessDetector= 500
+Sharc QQQ
+ Z= -70 mm
+ R= 0 mm
+ Phi= 92 deg
+ ThicknessDetector= 500
+Sharc QQQ
+ Z= -70 mm
+ R= 0 mm
+ Phi= 182 deg
+ ThicknessDetector= 500
+Sharc QQQ
+ Z= -70 mm
+ R= 0 mm
+ Phi= 272 deg
+ ThicknessDetector= 500
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 110
- PHI= 220
- Shape= Square
+ R= 110 mm
+ THETA= 110 deg
+ PHI= 320 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 67
- PHI= 320
- Shape= Square
+ R= 110 mm
+ THETA= 110 deg
+ PHI= 220 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 67
- PHI= 220
- Shape= Square
+ R= 110 mm
+ THETA= 67 deg
+ PHI= 320 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 67
- PHI= 40
- Shape= Square
+ R= 110 mm
+ THETA= 67 deg
+ PHI= 220 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 67
- PHI= 140
- Shape= Square
+ R= 110 mm
+ THETA= 67 deg
+ PHI= 40 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 110
- PHI= 40
- Shape= Square
+ R= 110 mm
+ THETA= 67 deg
+ PHI= 140 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Miniball
- R= 110
- THETA= 125
- PHI= 140
- Shape= Square
+ R= 110 mm
+ THETA= 110 deg
+ PHI= 40 deg
+ Shape= Square
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Miniball
+ R= 110 mm
+ THETA= 125 deg
+ PHI= 140 deg
+ Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/e748/Analysis.cxx b/Projects/e748/Analysis.cxx
index a66d024f8f42a54f16bfe368c0f97b21132ede4e..376c1dd901c721db230647806bffb4ec1d9645ff 100644
--- a/Projects/e748/Analysis.cxx
+++ b/Projects/e748/Analysis.cxx
@@ -36,11 +36,9 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){
- M2= (TMust2Physics*) m_DetectorManager->GetDetector("M2Telescope");
- CATS= ( TCATSPhysics*) m_DetectorManager->GetDetector("CATSDetector");
-
- //Initial=new TInitialConditions();
-
+ M2 = (TMust2Physics*) m_DetectorManager->GetDetector("M2Telescope");
+ CATS = (TCATSPhysics*) m_DetectorManager->GetDetector("CATSDetector");
+ ModularLeaf = (TModularLeafPhysics*) m_DetectorManager->GetDetector("ModularLeaf");
InitOutputBranch();
InitInputBranch();
/* Rand = TRandom3(); */
@@ -62,17 +60,29 @@ void Analysis::Init(){
TimeCorr=0;
TargetThickness = m_DetectorManager->GetTargetThickness();
// Energy loss table: the G4Table are generated by the simulation
- He3CD2 = EnergyLoss("Example/He3_CD2.G4table","G4Table",101 );
+ He3CD2 = EnergyLoss("Example/He3_CD2.G4table","G4Table",10 );
He3Al = EnergyLoss("Example/He3_Al.G4table","G4Table",10);
He3Si = EnergyLoss("Example/He3_Si.G4table","G4Table",10);
- //Li11CD2 = EnergyLoss("ExampleLi11_CD2.G4table","G4Table",100);
+ BeamCD2 = EnergyLoss("Be12_CD2.G4table","G4Table",10);
+ BeamMylar = EnergyLoss("Be12_Mylar.G4table","G4Table",10);
+ BeamIsobutane = EnergyLoss("Be12_iC4H10.G4table","G4Table",10);
+
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
// Reinitiate calculated variable
ReInitValue();
- // Get the Init information on beam position and energy
+
+ // Calculate Run Number
+ static string filename ;
+ filename = RootInput::getInstance()->GetChain()->GetFile()->GetName();
+ size_t minor_pos = filename.rfind("_");
+ run_minor = atoi(filename.substr(minor_pos+1,1).c_str());
+ filename = filename.substr(0,minor_pos);
+ size_t major_pos = filename.rfind("_");
+ run_major = atoi(filename.substr(major_pos+1,4).c_str());
+ // Get the Init information on beam position and energy
// and apply by hand the experimental resolution
// This is because the beam diagnosis are not simulated
// PPAC position resolution on target is assumed to be 1mm
@@ -82,9 +92,6 @@ void Analysis::TreatEvent(){
TVector3 BeamDirection = CATS->GetBeamDirection();
// Beam energy is measured using F3 and F2 plastic TOF
//double BeamEnergy= myReaction-> GetBeamEnergy()* MeV;
- double BeamEnergy= 30*12* MeV;
- //BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
- myReaction->SetBeamEnergy(BeamEnergy);
//////////////////////////// LOOP on MUST2 Hit //////////////////
for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
/* //Part 0 : Get the usefull Data */
@@ -95,37 +102,70 @@ void Analysis::TreatEvent(){
Si_X_M2 = X ;
Si_Y_M2 = Y ;
- /* // Forward Telescope Only */
- // if(TelescopeNumber<6){
- /* // All Telescopes */
if(TelescopeNumber<9){
- // if(TelescopeNumber<6){
DetectorNumber = TelescopeNumber ;
- /* // Part 1 : Impact Angle */
+ /* // Part 1 : Impact Angle */
ThetaM2Surface = 0;
ThetaNormalTarget = 0;
if(XTarget>-1000 && YTarget>-1000){
TVector3 BeamImpact(XTarget,YTarget,0);
- //TVector3 BeamImpact(0,0,0);
-
TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
ThetaLab = HitDirection.Angle( BeamDirection );
-
-
ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
- static double BeamSpeed = 74.2;// mm per nano
- static double ZCats1 = 1458; //mm
- BeamLength = ZCats1/(1-sin(BeamDirection.Angle(TVector3(0,0,1))));
- ParticleLength = HitDirection.Mag();
- TimeCorr=BeamLength/BeamSpeed ;
-
+
+ // Beam Energy from Cav Time of Flight //
+
+ // Beam speed from Beam Energy
+
+ // double BeamSpeed = 10.8727 + ModularLeaf->GetCalibratedValue("T_CATS1_CAV")*0.276825; // mm/ns
+ //double BeamSpeed = 5.17952 + ModularLeaf->GetCalibratedValue("T_CATS1_CAV")*0.305315; // mm/ns
+ double BeamSpeed = 11.0476 + ModularLeaf->GetCalibratedValue("T_CATS1_CAV")*0.278917; // mm/ns
+ // Beam Energy before CATS1
+ static double c2 = 299.792458*299.792458;// mm/ns
+ double gamma = 1./sqrt(1-BeamSpeed*BeamSpeed/c2);
+ BeamEnergy= 11200.962140*(gamma-1);
+ double BeamAngle= BeamDirection.Angle(TVector3(0,0,1));
+
+ // Beam Energy and speed after CATS1
+ double BeamEnergyC1 = BeamMylar.Slow(BeamEnergy,1.2*micrometer,BeamAngle);
+ BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3.,BeamAngle);
+ BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,0.9*micrometer,BeamAngle);
+ BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3,BeamAngle);
+ BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,0.9*micrometer,BeamAngle);
+ BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3.,BeamAngle);
+ BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,1.2*micrometer,BeamAngle);
+
+ double gammaC1 = (BeamEnergyC1+11200.962140) / 11200.962140 ;
+ double BeamSpeedC1 = sqrt(c2*(1-1/(gammaC1*gammaC1)));
+ TVector3 C1toC2 = TVector3(CATS->PositionX[1],CATS->PositionY[1],CATS->PositionZ[1])
+ - TVector3(CATS->PositionX[0],CATS->PositionY[0],CATS->PositionZ[0]) ;
+ TimeCorr = C1toC2.Mag()/BeamSpeedC1;
+
+ // Beam Energy and speed after CATS2
+ double BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC1,1.2*micrometer,BeamAngle);
+ BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3.,BeamAngle);
+ BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,0.9*micrometer,BeamAngle);
+ BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3,BeamAngle);
+ BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,0.9*micrometer,BeamAngle);
+ BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3.,BeamAngle);
+ BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,1.2*micrometer,BeamAngle);
+
+ double gammaC2 = (BeamEnergyC2+11200.962140) / 11200.962140;
+ double BeamSpeedC2 = sqrt(c2*(1-1/(gammaC2*gammaC2)));
+ TVector3 C2toTarget = BeamImpact-TVector3(CATS->PositionX[1],CATS->PositionY[1],CATS->PositionZ[1]);
+ TimeCorr += C2toTarget.Mag()/BeamSpeedC2;
+
+ // slow down beam inside the target
+ BeamEnergy = BeamCD2.Slow(BeamEnergyC2,TargetThickness*0.5,BeamDirection.Angle(TVector3(0,0,1)));
+ myReaction->SetBeamEnergy(BeamEnergy);
+
}
@@ -135,8 +175,7 @@ void Analysis::TreatEvent(){
ThetaNormalTarget = -1000 ;
}
-
-/* // Part 2 : Impact Energy */
+/* // Part 2 : Impact Energy */
Energy = ELab = E_M2 = 0;
Si_E_M2 = M2->Si_E[countMust2];
CsI_E_M2= M2->CsI_E[countMust2];
@@ -194,8 +233,13 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("TimeCorr",&TimeCorr,"TimeCorr/D");
RootOutput::getInstance()->GetTree()->Branch("BeamLength",&BeamLength,"BeamLength/D");
RootOutput::getInstance()->GetTree()->Branch("ParticleLength",&ParticleLength,"ParticleLength/D");
+ RootOutput::getInstance()->GetTree()->Branch("BeamEnergy",&BeamEnergy,"BeamEnergy/D");
+
RootOutput::getInstance()->GetTree()->Branch("GATCONF",&vGATCONF,"GATCONF/s");
+ RootOutput::getInstance()->GetTree()->Branch("RunMajor",&run_major,"RunMajor/I");
+ RootOutput::getInstance()->GetTree()->Branch("RunMinor",&run_minor,"RunMinor/I");
+
/* RootOutput::getInstance()->GetTree()->Branch("ADC_CHIO_V15",&vADC_CHIO_V15,"ADC_CHIO_V15/s");
RootOutput::getInstance()->GetTree()->Branch("ADC_VOIE_29",&vADC_VOIE_29,"ADC_VOIE_29/s");
RootOutput::getInstance()->GetTree()->Branch("CHIO",&vCHIO,"CHIO/s");
@@ -276,6 +320,9 @@ void Analysis::ReInitValue(){
TimeCorr=-1000;
BeamLength=-1000;
ParticleLength=-1000;
+ run_minor=-1000;
+ run_major=-1000;
+ BeamEnergy=-1000;
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/e748/Analysis.h b/Projects/e748/Analysis.h
index c614c7ae4ce79b8ba8ee573c5ee273d4305e32d5..6bac47cdcb550985ba913f5ed81119a7689000a6 100644
--- a/Projects/e748/Analysis.h
+++ b/Projects/e748/Analysis.h
@@ -24,7 +24,8 @@
#include"NPVAnalysis.h"
#include "TMust2Physics.h"
#include "TCATSPhysics.h"
- #include "TInitialConditions.h"
+#include "TModularLeafPhysics.h"
+#include "TInitialConditions.h"
#include "NPEnergyLoss.h"
#include "NPReaction.h"
#include "TRandom3.h"
@@ -70,13 +71,21 @@ class Analysis: public NPL::VAnalysis{
double TargetThickness;
double OriginalThetaLab;
double OriginalELab;
+ double BeamEnergy;
+ int run_major;
+ int run_minor;
NPL::EnergyLoss He3CD2 ;
NPL::EnergyLoss He3Al ;
NPL::EnergyLoss He3Si ;
- NPL::EnergyLoss Li11CD2 ;
+ NPL::EnergyLoss BeamCD2 ;
+ NPL::EnergyLoss BeamMylar ;
+ NPL::EnergyLoss BeamIsobutane ;
+
+
TMust2Physics* M2;
TCATSPhysics* CATS;
+ TModularLeafPhysics* ModularLeaf;
TInitialConditions* Initial;
//other variables
Short_t vADC_CHIO_V15;
diff --git a/Projects/e748/Calibration/Energy/ExtractRawHisto_E.C b/Projects/e748/Calibration/Energy/ExtractRawHisto_E.C
index 29d419a57ea364a41dc422ea959b4d7dd03326e1..839018e011adacaf55aecd51fb21c77dbb090082 100644
--- a/Projects/e748/Calibration/Energy/ExtractRawHisto_E.C
+++ b/Projects/e748/Calibration/Energy/ExtractRawHisto_E.C
@@ -19,11 +19,11 @@
#include "/home/muvi/e748/nptool/NPLib/Detectors/MUST2/TMust2Data.h"
#include "/home/muvi/e748/nptool/NPLib/include/RootInput.h"
-#define NBTELESCOPE 8
+#define NBTELESCOPE 4
#define NBSTRIPS 128
#define NBSILI 16
-void ExtractMust2Histos(const char* fname = "run_0019")
+void ExtractMust2Histos(const char* fname = "run_1031")
{
RootInput* Input = RootInput::getInstance("RunToTreat.txt");
diff --git a/Projects/e748/configs/ConfigMust2.dat b/Projects/e748/configs/ConfigMust2.dat
index 247888b9066a495290646e88bdb2f6ae7c65bd4c..97ac0c84fa8b1939aaffe20273ef57a681982097 100644
--- a/Projects/e748/configs/ConfigMust2.dat
+++ b/Projects/e748/configs/ConfigMust2.dat
@@ -69,6 +69,6 @@ ConfigMust2
SI_X_E_RAW_THRESHOLD 8270
SI_Y_E_RAW_THRESHOLD 8120
CSI_E_RAW_THRESHOLD 8250
- CSI_SIZE 40
- TAKE_T_Y
- TAKE_E_Y
+ CSI_SIZE 34
+ TAKE_T_X
+ TAKE_E_X
diff --git a/Projects/e748/e748.detector b/Projects/e748/e748.detector
index 302a7fe4f7ac6751f3ca066a0e30103a20055f31..90336a06ac1a2a26bb9f421b822cd4081433e7eb 100644
--- a/Projects/e748/e748.detector
+++ b/Projects/e748/e748.detector
@@ -100,18 +100,24 @@ M2Telescope
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CATSDetector
- X1_Y1= -35.86 -34.76 -1457 mm
- X28_Y1= 35.26 -34.76 -1457 mm
- X1_Y28= -35.86 36.36 -1457 mm
- X28_Y28= 35.26 36.36 -1457 mm
+ X1_Y1= -35.86 -35.26 -1457 mm
+ X28_Y1= 35.26 -35.26 -1457 mm
+ X1_Y28= -35.86 35.86 -1457 mm
+ X28_Y28= 35.26 35.86 -1457 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CATSDetector
- X1_Y1= 34.86 -35.36 -949 mm
- X28_Y1= -36.26 -35.36 -949 mm
- X1_Y28= 34.66 35.76 -949 mm
- X28_Y28= -36.26 35.76 -949 mm
+ X1_Y1= 34.86 -35.86 -949 mm
+ X28_Y1= -36.26 -35.86 -949 mm
+ X1_Y28= 34.86 35.26 -949 mm
+ X28_Y28= -36.26 35.26 -949 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+ModularLeaf
+ DefaultValue= -1000
+ Leafs= QPlast QCaviar CHIO T_CATS1_CAV T_CATS1_2 T_MUVI_CATS1 T_PL_CATS2 EXL_HF T_PL_CHIO T_PL_CATS1
+
Chio_dig
Pos= -30. -30. 600. mm
@@ -125,21 +131,18 @@ Chio_an
%% warning all parameters need to be there for routine to exit
%Plastic
- THETA= 0 deg
- PHI= 0 deg
- R= 0 mm
- Thickness= 20 mm
- Shape= Square
- Height= 30 mm
- Width= 30 mm
- Scintillator= BC400
- LeadThickness= 0 mm
-
+% THETA= 0 deg
+ % PHI= 0 deg
+ % R= 0 mm
+ % Thickness= 20 mm
+ %Shape= Square
+ %Height= 30 mm
+ %Width= 30 mm
+ %Scintillator= BC400
+ %LeadThickness= 0 mm
+ %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-EXL
- POS= -20 30 -20 cm
+%EXL
+ % POS= -20 30 -20 cm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-ModularLeaf
- DefaultValue= -1000
- Leafs= T_CATS1_CAV T_CATS1_2 T_MUVI_CATS1 T_PL_CATS2 EXL_HF T_PL_CHIO T_PL_CATS1
diff --git a/Projects/macros/GeometricalEfficiency.C b/Projects/macros/GeometricalEfficiency.C
index 99ba200fc583cb1b1dbf12b851e083302ced89d6..9655b7afddfbe9ddacc42377bb0394f7a002d5d7 100644
--- a/Projects/macros/GeometricalEfficiency.C
+++ b/Projects/macros/GeometricalEfficiency.C
@@ -48,8 +48,7 @@
using namespace std;
-void GeometricalEfficiency(const char * fname = "46Ar_pd_gs"){
- gROOT->SetStyle("pierre_style");
+void GeometricalEfficiency(const char * fname = "myResult"){
// Open output ROOT file from NPTool simulation run
TString path = gSystem->Getenv("NPTOOL");
path += "/Outputs/Simulation/";
@@ -78,6 +77,7 @@ void GeometricalEfficiency(const char * fname = "46Ar_pd_gs"){
int nentries = tree->GetEntries();
for (int i = 0; i < nentries; i++) {
tree->GetEntry(i);
+
// Fill histos
hEmittTheta->Fill(initCond->GetThetaLab_WorldFrame(0));
hEmittThetaCM->Fill(initCond->GetThetaCM(0));
@@ -113,13 +113,13 @@ void GeometricalEfficiency(const char * fname = "46Ar_pd_gs"){
c4->SetRightMargin(0.03);
TH1F* SolidACM = new TH1F(*hDetecThetaCM);
SolidACM->Sumw2();
- TF1* C = new TF1("C",Form("1./(2*%f*sin(x*%f/180.)*1*%f/180)",M_PI,M_PI,M_PI),0,90);
+ TF1* C = new TF1("C",Form("1./(2*%f*sin(x*%f/180.)*1*%f/180)",M_PI,M_PI,M_PI),0,180);
SolidACM->Divide(hEmittThetaCM);
SolidACM->Divide(C,1);
SolidACM->Draw();
SolidACM->GetXaxis()->SetTitle("#theta_{CM} (deg)");
SolidACM->GetYaxis()->SetTitle("d#Omega (sr) ");
- TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *1*%f/180.",M_PI,M_PI,M_PI),0,90);
+ TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *1*%f/180.",M_PI,M_PI,M_PI),0,180);
f->Draw("SAME");
c4->Update();