diff --git a/NPSimulation/Detectors/Vendeta/Vendeta.cc b/NPSimulation/Detectors/Vendeta/Vendeta.cc index d492843039fbafb1644295102113e2846cd9799c..369206b1b7371a47c8dc9f4b0b4e5e891eee9660 100644 --- a/NPSimulation/Detectors/Vendeta/Vendeta.cc +++ b/NPSimulation/Detectors/Vendeta/Vendeta.cc @@ -60,12 +60,13 @@ using namespace CLHEP; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... namespace Vendeta_NS{ // Energy and time Resolution - const double EnergyThreshold = 0.01*MeV; - const double ResoTime = 0.454*ns ; // reso Cf - /* const double ResoTime = 0.61*ns ; // reso U8 */ + const double EnergyThreshold = 0.02*MeV; + /* const double ResoTime = 0.454*ns ; // reso Cf */ + double ResoTime = 0.61*ns ; // reso U8 /* const double ResoTime = 0.*ns ; // reso U8 */ const double ResoEnergyLG = 0.43*MeV ; - const double ResoEnergyHG = 0.255*MeV ; + /* const double ResoEnergyHG = 0.255*MeV ; */ + const double ResoEnergyHG = 0.2*MeV ; const double Thickness = 51.*mm ; const double Radius = 127./2*mm ; @@ -83,7 +84,7 @@ Vendeta::Vendeta(){ m_VendetaDetector = 0; m_SensitiveCell = 0; m_MecanicalStructure = 0; - m_Build_MecanicalStructure = 0; + m_Build_MecanicalStructure = 1; m_LeadShield = 0; m_BuildLeadShield = 1; @@ -278,7 +279,7 @@ void Vendeta::ReadConfiguration(NPL::InputParser parser){ vector<string> cart = {"POS"}; vector<string> sphe = {"R","Theta","Phi"}; - + Vendeta_NS::ResoTime = blocks[0]->GetDouble("SigmaT","ns"); for(unsigned int i = 0 ; i < blocks.size() ; i++){ if(blocks[i]->HasTokenList(cart)){ if(NPOptionManager::getInstance()->GetVerboseLevel()) @@ -294,6 +295,7 @@ void Vendeta::ReadConfiguration(NPL::InputParser parser){ double Theta = blocks[i]->GetDouble("Theta","deg"); double Phi = blocks[i]->GetDouble("Phi","deg"); AddDetector(R,Theta,Phi); + } else{ cout << "ERROR: check your input file formatting " << endl; @@ -308,6 +310,14 @@ void Vendeta::ReadConfiguration(NPL::InputParser parser){ // Construct detector and inialise sensitive part. // Called After DetecorConstruction::AddDetector Method void Vendeta::ConstructDetector(G4LogicalVolume* world){ + + // Air density at Los Alamos (temp ~ 25 degres) + double density = 0.90*mg/cm3; + cout <<"units "<< mg << " " << cm3 << endl; + G4Material* Air = MaterialManager::getInstance()->GetMaterialFromLibrary("Air"); + + world->SetMaterial(Air); + for (unsigned short i = 0 ; i < m_R.size() ; i++) { G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ; @@ -338,16 +348,8 @@ void Vendeta::ConstructDetector(G4LogicalVolume* world){ BuildMecanicalStructure()->MakeImprint(world,PosMeca,RotMeca); } - - G4double platformWidth = 3.0 * m; - G4double platformLength = 3.0 * m; - G4double platformThickness = 1.2 * cm; - G4ThreeVector *platformPosition = new G4ThreeVector(0,0,-1000); - G4Box* platformSolid = new G4Box("PlatformSolid", platformWidth / 2.0, platformThickness / 2.0, platformLength / 2.0); - G4LogicalVolume* platformLogical = new G4LogicalVolume(platformSolid, m_Al, "PlatformLogical"); - G4PVPlacement* platformPhysical = new G4PVPlacement(0, G4ThreeVector(0,-1000,0), platformLogical, "PlatformPhysical",world,0,false); - } + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... // Add Detector branch to the EventTree. // Called After DetecorConstruction::AddDetector Method @@ -374,28 +376,28 @@ void Vendeta::ReadSensitive(const G4Event* ){ for(unsigned int i = 0 ; i < size ; i++){ vector<unsigned int> level = Scorer->GetLevel(i); /* double Energy =Scorer->GetEnergy(i); */ - double EnergyLG = RandGauss::shoot(Scorer->GetEnergy(i),Vendeta_NS::ResoEnergyLG); - double EnergyHG = RandGauss::shoot(Scorer->GetEnergy(i),Vendeta_NS::ResoEnergyHG); + double sigmaE = Scorer->GetEnergy(i) * 0.1897 + 0.04; + + double EnergyLG = RandGauss::shoot(Scorer->GetEnergy(i), sigmaE); + double EnergyHG = RandGauss::shoot(Scorer->GetEnergy(i), sigmaE); + //Convert enrgy deposit to Light Output thanks the F. Pino et al. formula double LightOutHG = 0.62*EnergyHG-1.3*(1-exp(-0.39*pow(EnergyHG,0.97))); // F. Pino et al double LightOutLG = 0.62*EnergyLG-1.3*(1-exp(-0.39*pow(EnergyLG,0.97))); // F. Pino et al // Apply Resolution on Charge measured on Vendeta data - /* LightOut = RandGauss::shoot(LightOut, Vendeta_NS::ResoEnergy); */ - if( EnergyHG > Vendeta_NS::EnergyThreshold){ + if( LightOutHG > Vendeta_NS::EnergyThreshold){ double Time = RandGauss::shoot(Scorer->GetTime(i),Vendeta_NS::ResoTime); int DetectorNbr = level[0]-1; - // Filling HG m_Event->SetHGDetectorNbr(DetectorNbr); - m_Event->SetHGQ1(LightOutHG*50000); - m_Event->SetHGQ2(LightOutHG*50000); + m_Event->SetHGQ1(LightOutHG); + m_Event->SetHGQ2(LightOutHG); m_Event->SetHGTime(Time); m_Event->SetHGQmax(0); - // Filling LG m_Event->SetLGDetectorNbr(DetectorNbr); - m_Event->SetLGQ1(LightOutLG*20000); - m_Event->SetLGQ2(LightOutLG*20000); + m_Event->SetLGQ1(LightOutLG); + m_Event->SetLGQ2(LightOutLG); m_Event->SetLGTime(Time); m_Event->SetLGQmax(0);