From 40e268fb37f2635d8d898e7035f9dacf66657eff Mon Sep 17 00:00:00 2001 From: matta <matta@npt> Date: Fri, 5 Mar 2010 10:05:47 +0000 Subject: [PATCH] * Adding new option to Plastic classes in NPS - Now can deal with squarred and cylindrical shape - Now can deal with spherical and cartesian coordinate --- .../DetectorConfiguration/Riken_65mm.detector | 28 +- Inputs/EventGenerator/10He.reaction | 2 +- Inputs/EventGenerator/alpha.source | 2 +- NPAnalysis/10He_Riken/src/Analysis.cc | 5 +- NPSimulation/include/Plastic.hh | 36 +- NPSimulation/src/EventGeneratorPhaseSpace.cc | 4 +- .../src/EventGeneratorTransfertToResonance.cc | 6 +- NPSimulation/src/Plastic.cc | 315 +++++++++++++----- 8 files changed, 297 insertions(+), 101 deletions(-) diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector index af4c56547..1931bbd73 100644 --- a/Inputs/DetectorConfiguration/Riken_65mm.detector +++ b/Inputs/DetectorConfiguration/Riken_65mm.detector @@ -17,12 +17,12 @@ GeneralTarget Z= 0 CryoTarget - THICKNESS= 3000 + THICKNESS= 4400 RADIUS= 45 TEMPERATURE= 26 PRESSURE= 1 MATERIAL= D2 - WINDOWSTHICKNESS= 10 + WINDOWSTHICKNESS= 12 WINDOWSMATERIAL= Mylar X= 0 Y= 0 @@ -149,6 +149,28 @@ D= 17.61 -9.85 104.11 ScintillatorPlastic %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Plastic + X= 0 + Y= -15 + Z= 318 + Shape= Square + Height= 70 + Width= 40 + Thickness= 20 + Scintillator= BC400 + LeadThickness= 0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Plastic + X= 0 + Y= -15 + Z= 343 + Thickness= 30 + Shape= Square + Height= 70 + Width= 40 + Scintillator= BC400 + LeadThickness= 2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Plastic THETA= 0 PHI= 0 R= 318 @@ -157,7 +179,7 @@ Plastic Scintillator= BC400 LeadThickness= 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Plastic +%Plastic THETA= 0 PHI= 0 R= 343 diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction index 138ab2f5c..2cbb5f098 100644 --- a/Inputs/EventGenerator/10He.reaction +++ b/Inputs/EventGenerator/10He.reaction @@ -19,7 +19,7 @@ TransfertToResonance ResonanceDecayA= 8 CrossSectionPath= 11Li(d,3He)10He.txt ShootLight= 1 - ShootHeavy= 0 + ShootHeavy= 1 ShootDecayProduct= 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source index 66b3898df..343524fe6 100644 --- a/Inputs/EventGenerator/alpha.source +++ b/Inputs/EventGenerator/alpha.source @@ -5,7 +5,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic EnergyLow= 0 - EnergyHigh= 300 + EnergyHigh= 25 HalfOpenAngleMin= 0 HalfOpenAngleMax= 90 x0= 0 diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc index 36770271e..857bb639e 100644 --- a/NPAnalysis/10He_Riken/src/Analysis.cc +++ b/NPAnalysis/10He_Riken/src/Analysis.cc @@ -169,7 +169,7 @@ RootOutput::getInstance()->GetList()->Add(myHist1D); if(M2 -> GetPositionOfInteraction(hit).Z() > 0) { - if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0 ) + if( M2 -> CsI_E[hit] < 0 && M2 -> Si_E[hit] > 0 ) { ELab[hit] = M2 -> Si_E[hit] ; @@ -197,8 +197,7 @@ RootOutput::getInstance()->GetList()->Add(myHist1D); ThetaN ); ThetaCM[hit] = He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ; -// ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) ; - ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( Init->GetICEmittedEnergy(hit) , Init->GetICEmittedAngleThetaLabIncidentFrame(hit) ) ; + ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) ; if(ThinSi -> Energy.size() > 0 ) if(cut3He_M2_SSSD->IsInside(ThinSi -> Energy[hit], M2 -> Si_E[hit]) ) diff --git a/NPSimulation/include/Plastic.hh b/NPSimulation/include/Plastic.hh index 9987fea5f..02f2020db 100644 --- a/NPSimulation/include/Plastic.hh +++ b/NPSimulation/include/Plastic.hh @@ -54,14 +54,24 @@ public: //////// Specific Function of this Class /////////// //////////////////////////////////////////////////// public: - // By Angle Method - void AddPlastic( G4double R , - G4double Theta , - G4double Phi , - G4double PlasticThickness , - G4double PlasticRadius , - G4String Scintillator , - G4double LeadThickness ); + // Cylindric plastic + void AddPlastic( G4double R , + G4double Theta , + G4double Phi , + G4double PlasticThickness , + G4double PlasticRadius , + G4String Scintillator , + G4double LeadThickness ); + + // Squared Plastic + void AddPlastic( G4double R , + G4double Theta , + G4double Phi , + G4double Height , + G4double Width , + G4double Thickness , + G4String Scintillator , + G4double LeadThickness ); void VolumeMaker(G4ThreeVector Det_pos, int DetNumber,G4LogicalVolume* world) ; //////////////////////////////////////////////////// @@ -111,12 +121,16 @@ private: //////////////////////////////////////////////////// private: - // if true a Lead plate is added in front or back of the detector + + + // Lead plate is added in front or back of the detector vector<double> m_LeadThickness ; vector<double> m_PlasticThickness ; - vector<double> m_PlasticRadius ; - + vector<double> m_PlasticRadius ; // cylindrical shape + vector<double> m_PlasticHeight ; // squared shape + vector<double> m_PlasticWidth ; // squared shape + // Used for "By Angle Definition" vector<G4double> m_R ; // | vector<G4double> m_Theta ; // > Spherical coordinate plastic volume center diff --git a/NPSimulation/src/EventGeneratorPhaseSpace.cc b/NPSimulation/src/EventGeneratorPhaseSpace.cc index 7f8e13b1f..ed612b9ad 100644 --- a/NPSimulation/src/EventGeneratorPhaseSpace.cc +++ b/NPSimulation/src/EventGeneratorPhaseSpace.cc @@ -247,13 +247,15 @@ while(ReadingStatus){ ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); else if ( A==1 && Z==1) ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("proton"); +// else if ( A==2 && Z==0) +// ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("neutron"); +// ParticleBuffer->SetPDGMass( ParticleBuffer->GetPDGMass*2); else ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ; ReactionProducts.push_back(ParticleBuffer); G4cout << "Decay Product: A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl; - G4cout << "Super " <<ParticleBuffer->GetPDGMass() << G4endl ; } else if (DataBuffer.compare(0, 21, "EndOfDecayProductList") == 0) { diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc index 17141d8f2..eb96fa7ec 100644 --- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc +++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc @@ -595,16 +595,16 @@ void EventGeneratorTransfertToResonance::ResonanceDecay( G4double EnergyHeavy G4ParticleDefinition* proton = G4ParticleTable::GetParticleTable()->FindParticle("proton"); - G4double M = parent -> GetPDGMass() + m_Reaction->GetExcitation() ; + G4double M = parent -> GetPDGMass() + m_Reaction->GetExcitation()*MeV ; G4double md = daughter -> GetPDGMass() ; G4double mn = neutron -> GetPDGMass() ; G4double mp = proton -> GetPDGMass() ; // Check that we are above threshold: // If the Resonnance go below the threshold, decay is forced at thereshold - if (M < md + NumberOfNeutrons*mm + NumberOfProtons*mp) + if (M < md + NumberOfNeutrons*mn + NumberOfProtons*mp) { - double NewExx = M - md - NumberOfNeutrons*mm - NumberOfProtons*mp ; + double NewExx = parent -> GetPDGMass() - md - NumberOfNeutrons*mn - NumberOfProtons*mp ; M = parent -> GetPDGMass() + NewExx ; } diff --git a/NPSimulation/src/Plastic.cc b/NPSimulation/src/Plastic.cc index ee6ef4b47..a8ea86188 100644 --- a/NPSimulation/src/Plastic.cc +++ b/NPSimulation/src/Plastic.cc @@ -27,6 +27,7 @@ #include <limits> //G4 Geometry object #include "G4Tubs.hh" +#include "G4Box.hh" //G4 sensitive #include "G4SDManager.hh" @@ -91,16 +92,40 @@ void Plastic::AddPlastic( G4double R , G4double PlasticRadius , G4String Scintillator , G4double LeadThickness ) -{ + { + + m_R.push_back(R) ; + m_Theta.push_back(Theta) ; + m_Phi.push_back(Phi) ; + m_PlasticThickness.push_back(PlasticThickness) ; + m_LeadThickness.push_back(LeadThickness) ; + m_Scintillator.push_back(Scintillator) ; + m_PlasticRadius.push_back(PlasticRadius) ; // cylindrical shape + m_PlasticHeight.push_back(-1) ; // squared shape + m_PlasticWidth.push_back(-1) ; // squared shape + } + +void Plastic::AddPlastic( G4double R , + G4double Theta , + G4double Phi , + G4double Height , + G4double Width , + G4double PlasticThickness , + G4String Scintillator , + G4double LeadThickness ) + { + m_R.push_back(R) ; + m_Theta.push_back(Theta) ; + m_Phi.push_back(Phi) ; + m_PlasticThickness.push_back(PlasticThickness) ; + m_LeadThickness.push_back(LeadThickness) ; + m_Scintillator.push_back(Scintillator) ; + m_PlasticRadius.push_back(-1) ; // cylindrical shape + m_PlasticHeight.push_back(Height) ; // squared shape + m_PlasticWidth.push_back(Width) ; // squared shape + + } - m_R.push_back(R) ; - m_Theta.push_back(Theta) ; - m_Phi.push_back(Phi) ; - m_PlasticThickness.push_back(PlasticThickness) ; - m_PlasticRadius.push_back(PlasticRadius) ; - m_LeadThickness.push_back(LeadThickness) ; - m_Scintillator.push_back(Scintillator) ; -} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -117,18 +142,24 @@ void Plastic::ReadConfiguration(string Path) string LineBuffer ; string DataBuffer ; - G4double Theta = 0 , Phi = 0 , R = 0 , Thickness = 0 , Radius = 0 , LeadThickness = 0; - G4String Scintillator ; + G4double Theta = 0 , Phi = 0 , R = 0 , Thickness = 0 , Radius = 0 , LeadThickness = 0, X = 0 , Y = 0 , Z = 0 , Width = 0 , Height = 0 ; + G4String Scintillator, Shape ; - bool check_Theta = false ; - bool check_Phi = false ; - bool check_R = false ; + bool check_Theta = false ; + bool check_Phi = false ; + bool check_R = false ; bool check_Thickness = false ; - bool check_Radius = false ; - bool check_LeadThickness = false ; + bool check_Radius = false ; + bool check_LeadThickness = false ; bool check_Scintillator = false ; - bool ReadingStatus = false ; - + bool check_Height = false ; + bool check_Width = false ; + bool check_Shape = false ; + bool check_X = false ; + bool check_Y = false ; + bool check_Z = false ; + bool ReadingStatus = false ; + while (!ConfigFile.eof()) { @@ -187,6 +218,41 @@ void Plastic::ReadConfiguration(string Path) cout << "R: " << R/mm << endl; } + //Position method + else if (DataBuffer.compare(0, 2, "X=") == 0) { + check_X = true; + ConfigFile >> DataBuffer ; + X = atof(DataBuffer.c_str()) ; + X = X * mm; + cout << "X: " << X / mm << endl; + } + + else if (DataBuffer.compare(0, 2, "Y=") == 0) { + check_Y = true; + ConfigFile >> DataBuffer ; + Y = atof(DataBuffer.c_str()) ; + Y = Y * mm; + cout << "Y: " << Y / mm << endl; + } + + else if (DataBuffer.compare(0, 2, "Z=") == 0) { + check_Z = true; + ConfigFile >> DataBuffer ; + Z = atof(DataBuffer.c_str()) ; + Z = Z * mm; + cout << "Z: " << Z / mm << endl; + } + + + //General + else if (DataBuffer.compare(0, 6, "Shape=") == 0) { + check_Shape = true; + ConfigFile >> DataBuffer ; + Shape = DataBuffer ; + cout << "Shape: " << Shape << endl; + } + + // Cylindrical shape else if (DataBuffer.compare(0, 7, "Radius=") == 0) { check_Radius = true; ConfigFile >> DataBuffer ; @@ -195,6 +261,24 @@ void Plastic::ReadConfiguration(string Path) cout << "Plastic Radius: " << Radius/mm << endl; } + // Squared shape + else if (DataBuffer.compare(0, 7, "Width=") == 0) { + check_Width = true; + ConfigFile >> DataBuffer ; + Width = atof(DataBuffer.c_str()) ; + Width = Width * mm; + cout << "Plastic Width: " << Width/mm << endl; + } + + else if (DataBuffer.compare(0, 7, "Height=") == 0) { + check_Height = true; + ConfigFile >> DataBuffer ; + Height = atof(DataBuffer.c_str()) ; + Height = Height * mm; + cout << "Plastic Height: " << Height/mm << endl; + } + + // Common else if (DataBuffer.compare(0, 10, "Thickness=") == 0) { check_Thickness = true; ConfigFile >> DataBuffer ; @@ -226,27 +310,57 @@ void Plastic::ReadConfiguration(string Path) ///////////////////////////////////////////////// // If All necessary information there, toggle out - if ( check_Theta && check_Phi && check_R && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator) + if ( ( check_Theta && check_Phi && check_R && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator && check_Shape) // Cylindrical case + || ( check_X && check_Y && check_Z && check_Thickness && check_Radius && check_LeadThickness && check_Scintillator ) + + || ( check_Theta && check_Phi && check_R && check_Thickness && check_Width && check_Height && check_LeadThickness && check_Scintillator && check_Shape ) // Squared case + || ( check_X && check_Y && check_Z && check_Thickness && check_Width && check_Height && check_LeadThickness && check_Scintillator ) + ) { - AddPlastic( R , - Theta , - Phi , - Thickness , - Radius , - Scintillator , - LeadThickness ); + + if (check_X && check_Y && check_Z) + { + R = sqrt (X*X+Y*Y+Z*Z) ; + Theta = acos(Z / (R) ) ; + Phi = atan2(Y,X) ; + } + + if (Shape == "Cylinder") + AddPlastic( R , + Theta , + Phi , + Thickness , + Radius , + Scintillator , + LeadThickness ); + + else if (Shape == "Square") + AddPlastic( R , + Theta , + Phi , + Height , + Width , + Thickness , + Scintillator , + LeadThickness ); // Reinitialisation of Check Boolean check_Theta = false ; - check_Phi = false ; + check_Phi = false ; check_R = false ; check_Thickness = false ; check_Radius = false ; - check_LeadThickness = false ; - check_Scintillator = false ; + check_LeadThickness = false ; + check_Scintillator = false ; + check_Height = false ; + check_Width = false ; + check_Shape = false ; + check_X = false ; + check_Y = false ; + check_Z = false ; ReadingStatus = false ; - cout << "///"<< endl ; + cout << "///"<< endl ; } } @@ -303,58 +417,103 @@ void Plastic::VolumeMaker(G4ThreeVector Det_pos, int DetNumber, G4LogicalVolume* // Definition of the volume containing the sensitive detector - if(m_PlasticThickness[i]>0 && m_PlasticRadius[i]>0) - { - G4Tubs* solidPlastic = new G4Tubs( Name , - 0 , - m_PlasticRadius[i] , - m_PlasticThickness[i]/2 , - 0*deg , - 360*deg ); - - G4LogicalVolume* logicPlastic = new G4LogicalVolume(solidPlastic, PlasticMaterial, Name+ "_Scintillator", 0, 0, 0); - logicPlastic->SetSensitiveDetector(m_PlasticScorer); + + + // Cylindrical Case + if(m_PlasticRadius[i]!=-1) + { + if(m_PlasticThickness[i]>0 && m_PlasticRadius[i]>0) + { + G4Tubs* solidPlastic = new G4Tubs( Name , + 0 , + m_PlasticRadius[i] , + m_PlasticThickness[i]/2 , + 0*deg , + 360*deg ); + + G4LogicalVolume* logicPlastic = new G4LogicalVolume(solidPlastic, PlasticMaterial, Name+ "_Scintillator", 0, 0, 0); + logicPlastic->SetSensitiveDetector(m_PlasticScorer); + + G4VisAttributes* PlastVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ; + logicPlastic->SetVisAttributes(PlastVisAtt) ; + + + + PVPBuffer = new G4PVPlacement( 0 , + Det_pos , + logicPlastic , + Name + "_Scintillator" , + world , + false , + 0 ); + } - G4VisAttributes* PlastVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ; - logicPlastic->SetVisAttributes(PlastVisAtt) ; - - - - PVPBuffer = new G4PVPlacement( 0 , - Det_pos , - logicPlastic , - Name + "_Scintillator" , - world , - false , - 0 ); - - + if(m_LeadThickness[i]>0&& m_PlasticRadius[i]>0) + { + G4Tubs* solidLead = new G4Tubs( Name+"_Lead" , + 0 , + m_PlasticRadius[i] , + m_LeadThickness[i]/2 , + 0*deg , + 360*deg ); + + G4LogicalVolume* logicLead = new G4LogicalVolume(solidLead, m_MaterialLead, Name+"_Lead", 0, 0, 0); + G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ; + logicLead->SetVisAttributes(LeadVisAtt) ; + + PVPBuffer = new G4PVPlacement( 0 , + Det_pos+(m_PlasticThickness[i]/2+m_LeadThickness[i]/2)*Det_pos.unit() , + logicLead , + Name+"_Lead" , + world , + false , + 0 ); + + } + } + + // Squared case + if(m_PlasticHeight[i]!=-1) + { + + if(m_PlasticThickness[i]>0 && m_PlasticHeight[i]>0 && m_PlasticWidth[i]>0) + { + G4Box* solidPlastic = new G4Box(Name, 0.5*m_PlasticWidth[i], 0.5*m_PlasticHeight[i], 0.5*m_PlasticThickness[i]); + G4LogicalVolume* logicPlastic = new G4LogicalVolume(solidPlastic, PlasticMaterial, Name+ "_Scintillator", 0, 0, 0); + logicPlastic->SetSensitiveDetector(m_PlasticScorer); + + G4VisAttributes* PlastVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)) ; + logicPlastic->SetVisAttributes(PlastVisAtt) ; + + PVPBuffer = new G4PVPlacement( 0 , + Det_pos , + logicPlastic , + Name + "_Scintillator" , + world , + false , + 0 ); + } + + if(m_LeadThickness[i]>0&& m_PlasticHeight[i]>0 && m_PlasticWidth[i]>0) + { + G4Box* solidLead = new G4Box(Name+"_Lead", 0.5*m_PlasticWidth[i], 0.5*m_PlasticHeight[i], 0.5*m_LeadThickness[i]); + + G4LogicalVolume* logicLead = new G4LogicalVolume(solidLead, m_MaterialLead, Name+"_Lead", 0, 0, 0); + G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ; + logicLead->SetVisAttributes(LeadVisAtt) ; + + PVPBuffer = new G4PVPlacement( 0 , + Det_pos+(m_PlasticThickness[i]/2+m_LeadThickness[i]/2)*Det_pos.unit() , + logicLead , + Name+"_Lead" , + world , + false , + 0 ); + + } } - - if(m_LeadThickness[i]>0&& m_PlasticRadius[i]>0) - { - G4Tubs* solidLead = new G4Tubs( Name+"_Lead" , - 0 , - m_PlasticRadius[i] , - m_LeadThickness[i]/2 , - 0*deg , - 360*deg ); - - G4LogicalVolume* logicLead = new G4LogicalVolume(solidLead, m_MaterialLead, Name+"_Lead", 0, 0, 0); - G4VisAttributes* LeadVisAtt = new G4VisAttributes(G4Colour(0.1, 0.1, 0.1)) ; - logicLead->SetVisAttributes(LeadVisAtt) ; - - PVPBuffer = new G4PVPlacement( 0 , - Det_pos+(m_PlasticThickness[i]/2+m_LeadThickness[i]/2)*Det_pos.unit() , - logicLead , - Name+"_Lead" , - world , - false , - 0 ); - - } } // Add Detector branch to the EventTree. -- GitLab