diff --git a/NPSimulation/Detectors/Actar/Actar.cc b/NPSimulation/Detectors/Actar/Actar.cc index 7f197a0fb5c944725f35e7b037b112fdf7fcf99e..548ab89a415d369f7f12b2844358f64709762208 100644 --- a/NPSimulation/Detectors/Actar/Actar.cc +++ b/NPSimulation/Detectors/Actar/Actar.cc @@ -1,18 +1,18 @@ /***************************************************************************** - * Copyright (C) 2009-2017 this file is part of the NPTool Project * + * Copyright (C) 2009-2017 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: Pierre Morfouace contact address: morfouac@nscl.msu.edu * + * Original Author: Pierre Morfouace contact address: morfouac@nscl.msu.edu * * * * Creation Date : September 2017 * * Last update : * *---------------------------------------------------------------------------* * Decription: * - * This class describe Actar simulation * + * This class describe Actar simulation * * * *---------------------------------------------------------------------------* * Comment: * @@ -52,7 +52,7 @@ // NPTool header #include "Actar.hh" -#include "SiliconScorers.hh" +#include "DSSDScorers.hh" #include "TPCScorers.hh" #include "RootOutput.h" #include "MaterialManager.hh" @@ -694,23 +694,17 @@ void Actar::ReadSensitive(const G4Event* event){ // Silicon // if(m_build_Silicon){ - ReducedData DataReduced; - - NPS::HitsMap<G4double*>* SiHitMap; - std::map<G4int, G4double**>::iterator Si_itr; - - G4int SiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID(m_SiliconCollectionID); - SiHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(SiCollectionID)); - - // Loop on the ThinSi map - for (Si_itr = SiHitMap->GetMap()->begin() ; Si_itr != SiHitMap->GetMap()->end() ; Si_itr++){ - G4double* Info = *(Si_itr->second); - double E_Si = RandGauss::shoot(Info[0],Actar_NS::ResoSilicon); + DSSDScorers::PS_Rectangle* SiScorer= (DSSDScorers::PS_Rectangle*) m_SiliconScorer->GetPrimitive(0); + unsigned int sizeSi = SiScorer->GetLengthMult(); + // Loop on the ThinSi map + for(unsigned int i = 0 ; i < sizeSi ; i++){ + DataReduced.clear(); + double E_Si = RandGauss::shoot(SiScorer->GetEnergyLength(i),Actar_NS::ResoSilicon); - int co=31; - int as=0; - int ag=0; - int ch=Info[7]; + int co = 31; + int as = 0; + int ag = 0; + int ch = SiScorer->GetStripLength(i); if(E_Si>Actar_NS::EnergyThreshold){ DataReduced.globalchannelid = ch+(ag<<7)+(as<<9)+(co<<11); @@ -720,7 +714,7 @@ void Actar::ReadSensitive(const G4Event* event){ m_EventReduced->CoboAsad.push_back(DataReduced); } // Clear Map for next event - SiHitMap->clear(); + SiScorer->clear(); } // CsI // @@ -733,8 +727,8 @@ void Actar::ReadSensitive(const G4Event* event){ // Loop on the CsI map for (CsI_itr = CsIHitMap->GetMap()->begin() ; CsI_itr !=CsIHitMap->GetMap()->end() ; CsI_itr++){ - G4double* Info = *(CsI_itr->second); - double E_CsI = RandGauss::shoot(Info[0],Actar_NS::ResoCsI); + //G4double* Info = *(CsI_itr->second); + //double E_CsI = RandGauss::shoot(Info[0],Actar_NS::ResoCsI); /*if(E_CsI>Actar_NS::EnergyThreshold){ m_Event->SetCsIEnergy(E_CsI); @@ -796,7 +790,7 @@ void Actar::InitializeScorers() { if(already_exist) return; - G4VPrimitiveScorer* SiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("SiliconScorer",0,Actar_NS::SiliconHeight,Actar_NS::SiliconWidth,1,1); + G4VPrimitiveScorer* SiScorer = new DSSDScorers::PS_Rectangle("SiliconScorer",0,Actar_NS::SiliconHeight,Actar_NS::SiliconWidth,1,1); m_SiliconScorer->RegisterPrimitive(SiScorer); /*G4VPrimitiveScorer* VamosSiScorer = new SILICONSCORERS::PS_Silicon_Rectangle("VamosSiliconScorer",0,Actar_NS::VamosSiliconHeight,Actar_NS::VamosSiliconWidth,1,1); diff --git a/NPSimulation/Detectors/AnnularS1/AnnularS1.cc b/NPSimulation/Detectors/AnnularS1/AnnularS1.cc index 5cbae4e7642fee8b2dcd1e1477a1c9f1a9b98085..fbad5cdbefd6463ae400508b5f716b6c4ccc799b 100644 --- a/NPSimulation/Detectors/AnnularS1/AnnularS1.cc +++ b/NPSimulation/Detectors/AnnularS1/AnnularS1.cc @@ -44,7 +44,8 @@ #include "MaterialManager.hh" #include "NPSDetectorFactory.hh" #include "AnnularS1.hh" -#include "SiliconScorers.hh" +#include "DSSDScorers.hh" +#include "InteractionScorers.hh" #include "TS1Data.h" #include "RootOutput.h" #include "NPOptionManager.h" @@ -299,28 +300,32 @@ void AnnularS1::InitializeRootOutput(){ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... // Read sensitive part and fill the Root tree. // Called at in the EventAction::EndOfEventAvtion -void AnnularS1::ReadSensitive(const G4Event* event){ +void AnnularS1::ReadSensitive(const G4Event*){ // Clear ROOT objects m_Event->Clear(); - NPS::HitsMap<G4double*>* SiliconHitMap; - std::map<G4int, G4double**>::iterator Silicon_itr; - - G4int SiliconCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("AnnularS1_Scorer/AnnularS1_Scorer"); - SiliconHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(SiliconCollectionID)); - - // Loop on the Silicon map - for (Silicon_itr = SiliconHitMap->GetMap()->begin() ; Silicon_itr != SiliconHitMap->GetMap()->end() ; Silicon_itr++){ - G4double* Info = *(Silicon_itr->second); - - double Energy = Info[0]; + DSSDScorers::PS_Annular* Scorer= (DSSDScorers::PS_Annular*) m_Scorer->GetPrimitive(0); + + // Loop on Silicon Ring Hit + unsigned int sizeRing = Scorer->GetRingMult(); + for(unsigned int i = 0 ; i < sizeRing ; i++){ + double Energy = Scorer->GetEnergyRing(i); if(Energy>EnergyThreshold){ - double Time = Info[1]; - int DetNbr = (int) Info[7]; - int StripTheta = (int) Info[8]; - int StripPhi = (int) Info[9]; - int StripQuadrant = (int) Info[10] - 1; + double Time = Scorer->GetTimeRing(i); + unsigned int DetNbr = Scorer->GetDetectorRing(i);; + unsigned int StripTheta = Scorer->GetStripRing(i); + + + // Check for associated Quadrant strip + int StripQuadrant = 0; + unsigned int sizeQ = Scorer->GetQuadrantMult(); + for(unsigned int q = 0 ; q < sizeQ ; q++){ + if(Scorer->GetDetectorQuadrant(q)==DetNbr){ + StripQuadrant = Scorer->GetStripQuadrant(q)-1; + break; + } + } m_Event->SetS1ThetaEDetectorNbr(DetNbr); m_Event->SetS1ThetaEStripNbr(StripTheta+StripQuadrant*NbrRingStrips); @@ -330,6 +335,20 @@ void AnnularS1::ReadSensitive(const G4Event* event){ m_Event->SetS1ThetaTStripNbr(StripTheta+StripQuadrant*NbrRingStrips); m_Event->SetS1ThetaTTime(RandGauss::shoot(Time, ResoTime)); + } + + } + + // Loop on Silicon Sector Hit + unsigned int sizeSector = Scorer->GetSectorMult(); + for(unsigned int i = 0 ; i < sizeSector ; i++){ + double Energy = Scorer->GetEnergyRing(i); + + if(Energy>EnergyThreshold){ + double Time = Scorer->GetTimeRing(i); + unsigned int DetNbr = Scorer->GetDetectorRing(i);; + unsigned int StripPhi = Scorer->GetStripSector(i); + m_Event->SetS1PhiEDetectorNbr(DetNbr); m_Event->SetS1PhiEStripNbr(StripPhi); m_Event->SetS1PhiEEnergy(RandGauss::shoot(Energy, ResoEnergy)); @@ -337,18 +356,8 @@ void AnnularS1::ReadSensitive(const G4Event* event){ m_Event->SetS1PhiTDetectorNbr(DetNbr); m_Event->SetS1PhiTStripNbr(StripPhi); m_Event->SetS1PhiTTime(RandGauss::shoot(Time, ResoTime)); - - // Interraction Coordinates - ms_InterCoord->SetDetectedPositionX(Info[2]) ; - ms_InterCoord->SetDetectedPositionY(Info[3]) ; - ms_InterCoord->SetDetectedPositionZ(Info[4]) ; - ms_InterCoord->SetDetectedAngleTheta(Info[5]/deg) ; - ms_InterCoord->SetDetectedAnglePhi(Info[6]/deg) ; - } } - // clear map for next event - SiliconHitMap->clear(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -360,7 +369,7 @@ void AnnularS1::InitializeScorers(){ if(already_exist) return; G4VPrimitiveScorer* AnnularScorer = - new SILICONSCORERS::PS_Silicon_Annular("AnnularS1_Scorer", + new DSSDScorers::PS_Annular("AnnularS1_Scorer", 2, ActiveWaferInnerRadius, ActiveWaferOutterRadius, @@ -371,7 +380,8 @@ void AnnularS1::InitializeScorers(){ NbQuadrant); m_Scorer->RegisterPrimitive(AnnularScorer); - + G4VPrimitiveScorer* InteractionScorer = new InteractionScorers::PS_Interactions("InteractionS1",ms_InterCoord,2); + m_Scorer->RegisterPrimitive(InteractionScorer); // Add All Scorer to the Global Scorer Manager G4SDManager::GetSDMpointer()->AddNewDetector(m_Scorer); } diff --git a/NPSimulation/Detectors/Helios2/Helios2.cc b/NPSimulation/Detectors/Helios2/Helios2.cc index c2e5f35b0faadae524726ae11441ed8023506dbe..aa2558c62d355add0ce918604f1ec05448417cc1 100644 --- a/NPSimulation/Detectors/Helios2/Helios2.cc +++ b/NPSimulation/Detectors/Helios2/Helios2.cc @@ -47,7 +47,8 @@ // NPTool header #include "Helios2.hh" -#include "SiliconScorers.hh" +#include "DSSDScorers.hh" +#include "InteractionScorers.hh" #include "RootOutput.h" #include "MaterialManager.hh" #include "NPSDetectorFactory.hh" @@ -117,9 +118,9 @@ G4LogicalVolume* Helios2::BuildSquareTube(){ Helios2_NS::SquareTubeSide*0.5,0.5*(Helios2_NS::SquareTubeExcess+Helios2_NS::WaferLength)); G4Tubs* tubs = new G4Tubs("Helios2_Box",0,Helios2_NS::SquareTubeRadius, - (Helios2_NS::SquareTubeExcess+Helios2_NS::WaferLength),0,360*deg); - - + (Helios2_NS::SquareTubeExcess+Helios2_NS::WaferLength),0,360*deg); + + G4RotationMatrix* R = new G4RotationMatrix(); G4ThreeVector P(0,0,0); G4SubtractionSolid* sub = new G4SubtractionSolid("Helios2_Sub",box,tubs,R,P); @@ -149,15 +150,15 @@ G4LogicalVolume* Helios2::BuildSiliconWafer(){ G4ThreeVector AWPos(0,0,0); G4RotationMatrix* AWRot = new G4RotationMatrix(); new G4PVPlacement(G4Transform3D(*AWRot,AWPos),m_ActiveWafer, - "Helios2_ActiveWafer",m_SiliconWafer, true, 0); + "Helios2_ActiveWafer",m_SiliconWafer, true, 0); m_ActiveWafer->SetSensitiveDetector(m_Helios2Scorer); m_SiliconWafer->SetVisAttributes(m_VisPassiveSilicon); m_ActiveWafer->SetVisAttributes(m_VisSilicon); - - } - + } + + return m_SiliconWafer; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -165,20 +166,20 @@ G4LogicalVolume* Helios2::BuildMagnet(){ if(!m_Magnet){ G4Tubs* tubs1 = new G4Tubs("Helios2_MainFull",0, Helios2_NS::MagnetOutterRadius,Helios2_NS::MagnetLength*0.5,0,360*deg); - + // Inner part of the Soleinoid minus the Target it self (placed in the world) G4SubtractionSolid* tubs = new G4SubtractionSolid("Helios_Main", - tubs1, Target::GetTarget()->GetTargetSolid(), new G4RotationMatrix() ,Target::GetTarget()->GetTargetPosition()); + tubs1, Target::GetTarget()->GetTargetSolid(), new G4RotationMatrix() ,Target::GetTarget()->GetTargetPosition()); G4Tubs* tubs2 = new G4Tubs("Helios2_Mag",Helios2_NS::MagnetInnerRadius, Helios2_NS::MagnetOutterRadius,Helios2_NS::MagnetLength*0.5,0,360*deg); G4Material* Fe= MaterialManager::getInstance()->GetMaterialFromLibrary("Fe"); G4Material* Vc= MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); - + m_Magnet= new G4LogicalVolume(tubs,Vc,"logic_Helios2_Main",0,0,0); G4LogicalVolume* Mag = new G4LogicalVolume(tubs2,Fe,"logic_Helios2_Magnet",0,0,0); - + Mag->SetVisAttributes(m_VisMagnet); m_Magnet->SetVisAttributes(G4VisAttributes::Invisible); // Place the Solenoid @@ -188,7 +189,7 @@ G4LogicalVolume* Helios2::BuildMagnet(){ new G4PVPlacement(G4Transform3D(*MagRot,MagPos), Mag, "Helios2_Magnet",m_Magnet,false,0); - + } return m_Magnet; } @@ -252,7 +253,7 @@ void Helios2::ConstructDetector(G4LogicalVolume* world){ fieldMgr->SetMinimumEpsilonStep( 1*mm); fieldMgr->SetMaximumEpsilonStep( 10*m ); fieldMgr->SetDeltaOneStep( 1 * mm ); - + // Place detectors and support inside it for (unsigned short i = 0 ; i < m_Z.size() ; i++) { G4ThreeVector DetPos; @@ -316,50 +317,44 @@ void Helios2::InitializeRootOutput(){ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... // Read sensitive part and fill the Root tree. // Called at in the EventAction::EndOfEventAvtion -void Helios2::ReadSensitive(const G4Event* event){ - m_Event->Clear(); +void Helios2::ReadSensitive(const G4Event* ){ + m_Event->Clear(); /////////// // Resistiverimeter scorer - NPS::HitsMap<G4double*>* ResistiveHitMap; - std::map<G4int, G4double**>::iterator Resistive_itr; - - G4int ResistiveCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("Helios2Scorer/Resistive"); - ResistiveHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(ResistiveCollectionID)); - - // Loop on the Resistive map - for (Resistive_itr = ResistiveHitMap->GetMap()->begin() ; Resistive_itr != ResistiveHitMap->GetMap()->end() ; Resistive_itr++){ - G4double* Info = *(Resistive_itr->second); - double EBack = RandGauss::shoot(Info[0]+Info[1],Helios2_NS::ResoEnergyBack); - double TBack = RandGauss::shoot(Info[2],Helios2_NS::ResoTime); - double EUp = RandGauss::shoot(Info[1],Helios2_NS::ResoEnergyFront); - double TUp = RandGauss::shoot(Info[2],Helios2_NS::ResoTime); - double EDw = RandGauss::shoot(Info[0],Helios2_NS::ResoEnergyFront); - double TDw = RandGauss::shoot(Info[2],Helios2_NS::ResoTime); - - if(EBack>Helios2_NS::EnergyThreshold){ - m_Event->SetEBack(Info[3],EBack); - m_Event->SetTBack(Info[3],TBack); - } - - if(EUp>Helios2_NS::EnergyThreshold){ - m_Event->SetEUp(Info[3],EUp); - m_Event->SetTUp(Info[3],TUp); - } - - if(EUp>Helios2_NS::EnergyThreshold){ - m_Event->SetEUp(Info[3],EUp); - m_Event->SetTUp(Info[3],TUp); + DSSDScorers::PS_Resistive* Scorer= (DSSDScorers::PS_Resistive*) m_Helios2Scorer->GetPrimitive(0); + + // Loop on the Back + unsigned int sizeBack = Scorer->GetBackMult(); + for(unsigned int i = 0 ; i < sizeBack ; i++){ + double EBack = RandGauss::shoot(Scorer->GetEnergyBack(i),Helios2_NS::ResoEnergyBack); + double TBack = RandGauss::shoot(Scorer->GetTimeBack(i),Helios2_NS::ResoTime); + if(EBack>Helios2_NS::EnergyThreshold){ + m_Event->SetEBack(Scorer->GetDetectorBack(i),EBack); + m_Event->SetTBack(Scorer->GetDetectorBack(i),TBack); + } } - - if(EDw>Helios2_NS::EnergyThreshold){ - m_Event->SetEDw(Info[3],EDw); - m_Event->SetTDw(Info[3],TDw); + // Loop on the Up + unsigned int sizeUp = Scorer->GetUpMult(); + for(unsigned int i = 0 ; i < sizeUp ; i++){ + double EUp = RandGauss::shoot(Scorer->GetEnergyUp(i),Helios2_NS::ResoEnergyFront); + double TUp = RandGauss::shoot(Scorer->GetTimeUp(i),Helios2_NS::ResoTime); + if(EUp>Helios2_NS::EnergyThreshold){ + m_Event->SetEUp(Scorer->GetDetectorUp(i),EUp); + m_Event->SetTUp(Scorer->GetDetectorUp(i),TUp); + } } - + + // Loop on the Down + unsigned int sizeDown = Scorer->GetDownMult(); + for(unsigned int i = 0 ; i < sizeDown ; i++){ + double EDw = RandGauss::shoot(Scorer->GetEnergyDown(i),Helios2_NS::ResoEnergyFront); + double TDw = RandGauss::shoot(Scorer->GetTimeDown(i),Helios2_NS::ResoTime); + if(EDw>Helios2_NS::EnergyThreshold){ + m_Event->SetEDw(Scorer->GetDetectorDown(i),EDw); + m_Event->SetTDw(Scorer->GetDetectorDown(i),TDw); + } } - // clear map for next event - ResistiveHitMap->clear(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -373,9 +368,12 @@ void Helios2::InitializeScorers() { return ; // Otherwise the scorer is initialised - G4VPrimitiveScorer* Resistive= new SILICONSCORERS::PS_Silicon_Resistive("Resistive",1,Helios2_NS::WaferLength,Helios2_NS::WaferWidth,1) ; + G4VPrimitiveScorer* Resistive= new DSSDScorers::PS_Resistive("Resistive",1,Helios2_NS::WaferLength,Helios2_NS::WaferWidth,1) ; //and register it to the multifunctionnal detector m_Helios2Scorer->RegisterPrimitive(Resistive); + G4VPrimitiveScorer* Inter = new InteractionScorers::PS_Interactions("Resistive",ms_InterCoord,1) ; + m_Helios2Scorer->RegisterPrimitive(Inter); + G4SDManager::GetSDMpointer()->AddNewDetector(m_Helios2Scorer) ; } diff --git a/NPSimulation/Scorers/DSSDScorers.cc b/NPSimulation/Scorers/DSSDScorers.cc index 30daf2e7f440a671a9f2f9a6a128d0c0269cb0af..dec8a8e7c7e670ee6e22ae3e9fe80e2aab3f57f7 100644 --- a/NPSimulation/Scorers/DSSDScorers.cc +++ b/NPSimulation/Scorers/DSSDScorers.cc @@ -327,17 +327,16 @@ void PS_Annular::PrintAll(){ } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... PS_Resistive::PS_Resistive(G4String name,G4int Level, G4double StripPlaneLength, G4double StripPlaneWidth, G4int NumberOfStripWidth,G4int depth) - :G4VPrimitiveScorer(name, depth),HCID(-1){ + :G4VPrimitiveScorer(name, depth){ m_StripPlaneLength = StripPlaneLength; m_StripPlaneWidth = StripPlaneWidth; m_NumberOfStripWidth = NumberOfStripWidth; m_StripPitchWidth = m_StripPlaneWidth / m_NumberOfStripWidth; m_Level = Level; - m_Position = G4ThreeVector(-1000,-1000,-1000); - m_DetectorNumber = -1; - m_StripWidthNumber = -1; - m_Index = -1 ; + t_Position = G4ThreeVector(-1000,-1000,-1000); + t_DetectorNumber = -1; + t_StripWidthNumber = -1; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -348,61 +347,58 @@ PS_Resistive::~PS_Resistive(){ G4bool PS_Resistive::ProcessHits(G4Step* aStep, G4TouchableHistory*){ // contain Energy Total, E1, E2, Time, DetNbr, and StripWidth - G4double* EnergyAndTime = new G4double[10]; - - EnergyAndTime[2] = aStep->GetPreStepPoint()->GetGlobalTime(); - - m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); - m_Position = aStep->GetPreStepPoint()->GetPosition(); - - // Interaction coordinates (used to fill the InteractionCoordinates branch) - EnergyAndTime[5] = m_Position.x(); - EnergyAndTime[6] = m_Position.y(); - EnergyAndTime[7] = m_Position.z(); - EnergyAndTime[8] = m_Position.theta(); - EnergyAndTime[9] = m_Position.phi(); - - m_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(m_Position); + t_Energy = aStep->GetTotalEnergyDeposit(); + t_Time = aStep->GetPreStepPoint()->GetGlobalTime(); + + t_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); + t_Position = aStep->GetPreStepPoint()->GetPosition(); + t_Position = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(t_Position); - m_StripWidthNumber = (int)((m_Position.x() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ; - // m_StripWidthNumber = m_NumberOfStripWidth - m_StripWidthNumber + 1 ; + t_StripWidthNumber = (int)((t_Position.x() + m_StripPlaneWidth / 2.) / m_StripPitchWidth ) + 1 ; // The energy is divided in two depending on the position // position along the resistive strip - double P = (m_Position.z())/(0.5*m_StripPlaneLength); - // Downstream Energy - EnergyAndTime[0] = aStep->GetTotalEnergyDeposit()*(1+P)*0.5; // downsteam collects MORE energy when P is positive - - // Upstream Energy - EnergyAndTime[1] = aStep->GetTotalEnergyDeposit()-EnergyAndTime[0]; - - EnergyAndTime[3] = m_DetectorNumber; - EnergyAndTime[4] = m_StripWidthNumber; + double P = (t_Position.z())/(0.5*m_StripPlaneLength); + // Energy + t_EnergyUp = aStep->GetTotalEnergyDeposit()*(1+P)*0.5; + t_EnergyDown = t_Energy-t_EnergyUp; //Rare case where particle is close to edge of silicon plan - if (m_StripWidthNumber > m_NumberOfStripWidth) m_StripWidthNumber = m_NumberOfStripWidth; + if (t_StripWidthNumber > m_NumberOfStripWidth) t_StripWidthNumber = m_NumberOfStripWidth; - m_Index = m_DetectorNumber * 1e3 + m_StripWidthNumber * 1e6; - // Check if the particle has interact before, if yes, add up the energies. - map<G4int, G4double**>::iterator it; - it= EvtMap->GetMap()->find(m_Index); - if(it!=EvtMap->GetMap()->end()){ - G4double* dummy = *(it->second); - EnergyAndTime[0]+=dummy[0]; - EnergyAndTime[1]+=dummy[1]; - } - EvtMap->set(m_Index, EnergyAndTime); + // Up + vector<DSSDData>::iterator it; + it = m_HitUp.find(DSSDData::CalculateIndex(t_DetectorNumber,t_StripWidthNumber)); + if(it!=m_HitUp.end()) + it->Add(t_EnergyUp); + else + m_HitUp.Set(t_EnergyUp,t_Time,t_DetectorNumber,t_StripWidthNumber); + + // Down + it = m_HitDown.find(DSSDData::CalculateIndex(t_DetectorNumber,t_StripWidthNumber)); + if(it!=m_HitDown.end()) + it->Add(t_EnergyDown); + else + m_HitDown.Set(t_EnergyDown,t_Time,t_DetectorNumber,t_StripWidthNumber); + + // Back + it = m_HitBack.find(DSSDData::CalculateIndex(t_DetectorNumber,t_StripWidthNumber)); + if(it!=m_HitBack.end()) + it->Add(t_Energy); + else + m_HitBack.Set(t_Energy,t_Time,t_DetectorNumber,1); + + + return TRUE; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void PS_Resistive::Initialize(G4HCofThisEvent* HCE){ - EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName()); - if (HCID < 0) { - HCID = GetCollectionID(0); - } - HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); + m_HitUp.clear(); + m_HitDown.clear(); + m_HitBack.clear(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -411,12 +407,9 @@ void PS_Resistive::EndOfEvent(G4HCofThisEvent*){ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void PS_Resistive::clear(){ - std::map<G4int, G4double**>::iterator MapIterator; - for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){ - delete *(MapIterator->second); - } - - EvtMap->clear(); + m_HitUp.clear(); + m_HitDown.clear(); + m_HitBack.clear(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -426,7 +419,4 @@ void PS_Resistive::DrawAll(){ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void PS_Resistive::PrintAll(){ - G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; - G4cout << " PrimitiveScorer " << GetName() << G4endl ; - G4cout << " Number of entries " << EvtMap->entries() << G4endl ; } diff --git a/NPSimulation/Scorers/DSSDScorers.hh b/NPSimulation/Scorers/DSSDScorers.hh index ec50ed784fc249fa5d79b4a78032f48ae4982083..a6722ee5a91b8556bac726ac74d6a3bc7d21eec3 100644 --- a/NPSimulation/Scorers/DSSDScorers.hh +++ b/NPSimulation/Scorers/DSSDScorers.hh @@ -291,21 +291,7 @@ namespace DSSDScorers { }; - //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - struct ResistiveOutput { - G4double upstreamEnergy; - G4double downstreamEnergy; - G4double globalTime; - G4double detectorNumber; - G4double stripWidthNumber; - G4double x; - G4double y; - G4double z; - G4double theta; - G4double phi; - }; - class PS_Resistive : public G4VPrimitiveScorer{ public: // with description @@ -326,22 +312,45 @@ namespace DSSDScorers { void PrintAll(); private: // Geometry of the detector - G4double m_StripPlaneLength; - G4double m_StripPlaneWidth; - G4int m_NumberOfStripWidth; - G4double m_StripPitchWidth; + double m_StripPlaneLength; + double m_StripPlaneWidth; + unsigned int m_NumberOfStripWidth; + double m_StripPitchWidth; // Level at which to find the copy number linked to the detector number G4int m_Level; - private: // inherited from G4VPrimitiveScorer - G4int HCID; - NPS::HitsMap<G4double*>* EvtMap; - + private: + // Up and Down are each extremities of the resistive strip + DSSDDataVector m_HitUp; + DSSDDataVector m_HitDown; + DSSDDataVector m_HitBack; + private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit) - G4ThreeVector m_Position ; - G4int m_DetectorNumber ; - G4int m_StripWidthNumber ; - G4long m_Index ; + G4ThreeVector t_Position; + double t_Energy; + double t_EnergyUp; + double t_EnergyDown; + double t_Time; + unsigned int t_DetectorNumber; + unsigned int t_StripWidthNumber; + unsigned int t_Index; + public: + unsigned int GetUpMult() {return m_HitUp.size();}; + unsigned int GetStripUp(const unsigned int& i){return m_HitUp[i]->GetStrip();}; + unsigned int GetDetectorUp(const unsigned int& i){return m_HitUp[i]->GetDetector();}; + double GetEnergyUp(const unsigned int& i){return m_HitUp[i]->GetEnergy();}; + double GetTimeUp(const unsigned int& i){return m_HitUp[i]->GetTime();}; + unsigned int GetDownMult() {return m_HitDown.size();}; + unsigned int GetStripDown(const unsigned int& i){return m_HitDown[i]->GetStrip();}; + unsigned int GetDetectorDown(const unsigned int& i){return m_HitDown[i]->GetDetector();}; + double GetEnergyDown(const unsigned int& i){return m_HitDown[i]->GetEnergy();}; + double GetTimeDown(const unsigned int& i){return m_HitDown[i]->GetTime();}; + unsigned int GetBackMult() {return m_HitBack.size();}; + unsigned int GetStripBack(const unsigned int& i){return m_HitBack[i]->GetStrip();}; + unsigned int GetDetectorBack(const unsigned int& i){return m_HitBack[i]->GetDetector();}; + double GetEnergyBack(const unsigned int& i){return m_HitBack[i]->GetEnergy();}; + double GetTimeBack(const unsigned int& i){return m_HitBack[i]->GetTime();}; + }; diff --git a/Outputs/Analysis/.gitignore b/Outputs/Analysis/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/Outputs/Simulation/.gitignore b/Outputs/Simulation/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391