diff --git a/Inputs/DetectorConfiguration/27Al.target b/Inputs/DetectorConfiguration/27Al.target new file mode 100644 index 0000000000000000000000000000000000000000..08f2199265cb8df825706729838e9529cbd3e883 --- /dev/null +++ b/Inputs/DetectorConfiguration/27Al.target @@ -0,0 +1,10 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +GeneralTarget +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +Target + THICKNESS= 0.3231 + RADIUS= 7.5 + MATERIAL= Al + X= 0 + Y= 0 + Z= 0 diff --git a/Inputs/DetectorConfiguration/W1.detector b/Inputs/DetectorConfiguration/W1.detector new file mode 100644 index 0000000000000000000000000000000000000000..7b2cf7f9cb21106020f482f190de30929f78e99a --- /dev/null +++ b/Inputs/DetectorConfiguration/W1.detector @@ -0,0 +1,7 @@ +%%%%%%% Telescope 1 %%%%%%% +W1 + THETA= 0 + PHI= 0 + R= 150 + BETA= 0 0 0 + VIS= all diff --git a/Inputs/EventGenerator/27Alppp.reaction b/Inputs/EventGenerator/27Alppp.reaction new file mode 100644 index 0000000000000000000000000000000000000000..3b5ebe181917689702972e8d16446e39c45524c0 --- /dev/null +++ b/Inputs/EventGenerator/27Alppp.reaction @@ -0,0 +1,25 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%% Reaction file for 11Li(d,3He)10He reaction %%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%Beam energy given in MeV ; Excitation in MeV +TransfertToResonance + Beam= 1H + Target= 27Al + Light= 1H + Heavy= 27Al + ExcitationEnergy= 13.0 + BeamEnergy= 18 + BeamEnergySpread= 0 + SigmaThetaX= 0 + SigmaPhiY= 0 + SigmaX= 0 + SigmaY= 0 + ResonanceWidth= 0 + ResonanceDecayZ= 12 + ResonanceDecayA= 26 + CrossSectionPath= flat.txt + ShootLight= 1 + ShootHeavy= 0 + ShootDecayProduct= 1 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source index 31cee696d6e75bd1f8e852af71fe4cb1b80f68fb..f840fa369c2183eaf8d988c6457c4dc48d487736 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= 0 + EnergyLow= 5 EnergyHigh= 5 - HalfOpenAngleMin= 158 - HalfOpenAngleMax= 180 + HalfOpenAngleMin= 0 + HalfOpenAngleMax= 20 x0= 0 y0= 0 z0= 0 diff --git a/NPSimulation/include/W1.hh b/NPSimulation/include/W1.hh new file mode 100644 index 0000000000000000000000000000000000000000..57eaf8dec747a1268d1a41af8bdcb2bdeefcb58d --- /dev/null +++ b/NPSimulation/include/W1.hh @@ -0,0 +1,177 @@ +/***************************************************************************** + * Copyright (C) 2009-2010 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : 12/01/11 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Define the W1 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef W1_h +#define W1_h 1 + +// C++ headers +#include <vector> + +// NPTool header +#include "VDetector.hh" + +// NPTool - ROOT headers +#include "TW1Data.h" + +using namespace std; + + + +class W1 : public VDetector +{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// +public: + W1(); + virtual ~W1(); + + //////////////////////////////////////////////////// + //////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// +public: + // Detector positionning + // By Position Method + void AddDetector(G4ThreeVector TL, G4ThreeVector BL, + G4ThreeVector BR, G4ThreeVector CT); + // By Angle Method + void AddDetector(G4double R, G4double Theta, G4double Phi, + G4double beta_u, G4double beta_v, G4double beta_w); + + // Effectively construct Volume + // Avoid to have two time same code for Angle and Point definition + void VolumeMaker(G4int DetecNumber, + G4ThreeVector pos, + G4RotationMatrix* rot, + G4LogicalVolume* world); + + + ///////////////////////////////////////// + //// Inherite from VDetector class ///// + ///////////////////////////////////////// +public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + // Called in DetecorConstruction::ReadDetextorConfiguration Method + void ReadConfiguration(string Path); + + // Construct detector and inialise sensitive part. + // Called After DetecorConstruction::AddDetector Method + void ConstructDetector(G4LogicalVolume* world); + + // Add Detector branch to the EventTree. + // Called After DetecorConstruction::AddDetector Method + void InitializeRootOutput(); + + // Initialize all scorers necessary for the detector + void InitializeScorers(); + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event); + + + //////////////////////////////////////////////////// + ///////////Event class to store Data//////////////// + //////////////////////////////////////////////////// +private: + TW1Data* m_Event; + + + //////////////////////////////////////////////////// + //////////////////// Scorers /////////////////////// + //////////////////////////////////////////////////// +private: + G4MultiFunctionalDetector* m_Scorer; + + + //////////////////////////////////////////////////// + //////////////////// Material ////////////////////// + //////////////////////////////////////////////////// +private: + // Declare all material used by theW1 detector + void InitializeMaterials(); + + // Vacuum + G4Material* m_MaterialVacuum; + // Si + G4Material* m_MaterialSilicon; + // Al + G4Material* m_MaterialAluminium; + // Iron + G4Material* m_MaterialIron; + + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// +private: + // True if Define by Position, False is Define by angle + vector<bool> m_DefinitionType; + + // Used for "By Point Definition" + vector<G4ThreeVector> m_X1_Y1; // Top Left Corner Position Vector + vector<G4ThreeVector> m_X1_Y16; // Bottom Left Corner Position Vector + vector<G4ThreeVector> m_X16_Y1; // Bottom Right Corner Position Vector + vector<G4ThreeVector> m_X16_Y16; // Center Corner Position Vector + + // Used for "By Angle Definition" + vector<G4double> m_R; // | + vector<G4double> m_Theta; // > Spherical coordinate of Strips Silicium Plate + vector<G4double> m_Phi; // | + + vector<G4double> m_beta_u; // | + vector<G4double> m_beta_v; // > Tilt angle of the Telescope + vector<G4double> m_beta_w; // | + + // Set to true if you want to see Telescope Frame in your visualisation + bool m_non_sensitive_part_visiualisation; +}; + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +namespace W1SQUARE +{ + // Energy/Time resolutions for the different layers + const G4double EnergyResolution = 0; // = 52keV of Resolution // Unit is MeV/2.35 +// const G4double ResoFirstStage = 0.0106 ;// = 52keV of Resolution // Unit is MeV/2.35 + const G4double TimeResolution = 0.638; // 1.5 ns (FWHM) + + // Geometry + const G4double FaceFront = 50*mm; + const G4double Length = 1*mm; + + // First stage +// const G4double AluStripThickness = 0.00000001*micrometer; + const G4double AluStripThickness = 0.4*micrometer; + const G4double SiliconThickness = 500*micrometer; + const G4double SiliconFace = 49.6*mm; + + // Characteristics + const G4int NbStrips = 16; + + // Starting at the front and going in direction of third stage + const G4double AluStripFront_PosZ = Length* -0.5 + 0.5*AluStripThickness; + const G4double Silicon_PosZ = AluStripFront_PosZ + 0.5*AluStripThickness + 0.5*SiliconThickness; + const G4double AluStripBack_PosZ = Silicon_PosZ + 0.5*SiliconThickness + 0.5*AluStripThickness; +} + +#endif diff --git a/NPSimulation/include/W1Scorers.hh b/NPSimulation/include/W1Scorers.hh new file mode 100644 index 0000000000000000000000000000000000000000..564dc95ac68b4630f34aaf3cc936496b236f7f29 --- /dev/null +++ b/NPSimulation/include/W1Scorers.hh @@ -0,0 +1,79 @@ +/***************************************************************************** + * Copyright (C) 2009-2010 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : 12/01/11 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class holds all the scorers needed by the W1 object * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef W1Scorer_h +#define W1Scorer_h 1 + +#include "G4VPrimitiveScorer.hh" +#include "G4THitsMap.hh" + +namespace W1SCORERS +{ + +class W1ScorerFrontStripNumber : public G4VPrimitiveScorer +{ +public: // with description + W1ScorerFrontStripNumber(G4String name, G4int depth = 0, G4int NumberOfStrip = 16); + virtual ~W1ScorerFrontStripNumber(); + +protected: // with description + virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*); + +public: + virtual void Initialize(G4HCofThisEvent*); + virtual void EndOfEvent(G4HCofThisEvent*); + virtual void Clear(); + virtual void DrawAll(); + virtual void PrintAll(); + +private: + G4int m_NumberOfStrip; + G4int HCID; + G4THitsMap<G4double>* EvtMap; +}; + +class W1ScorerBackStripNumber : public G4VPrimitiveScorer +{ +public: // with description + W1ScorerBackStripNumber(G4String name, G4int depth = 0, G4int NumberOfStrip = 16); + virtual ~W1ScorerBackStripNumber(); + +protected: // with description + virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*); + +public: + virtual void Initialize(G4HCofThisEvent*); + virtual void EndOfEvent(G4HCofThisEvent*); + virtual void Clear(); + virtual void DrawAll(); + virtual void PrintAll(); + +private: + G4int m_NumberOfStrip; + G4int HCID; + G4THitsMap<G4double>* EvtMap; +}; + + +} + +using namespace W1SCORERS; +#endif diff --git a/NPSimulation/src/DetectorConstruction.cc b/NPSimulation/src/DetectorConstruction.cc index 8ca7279ed95935fcf913c79e8160ee61b014da84..a35092c4a986621050ae0db4afad9c31be1276c9 100644 --- a/NPSimulation/src/DetectorConstruction.cc +++ b/NPSimulation/src/DetectorConstruction.cc @@ -50,6 +50,7 @@ #include "Plastic.hh" #include "Paris.hh" #include "Shield.hh" +#include "W1.hh" //Not G4 #include <cstdlib> @@ -167,8 +168,8 @@ void DetectorConstruction::ReadConfigurationFile(string Path) bool cDummy = false; bool cParis = false; // Paris Calorimeter bool cShield = false; // Paris Shield CsI + bool cW1 = false; // W1 Micron DSSD ////////////////////////////////////////////////////////////////////////////////////////// - // added by Nicolas [07/05/09] string GlobalPath = getenv("NPTOOL"); string StandardPath = GlobalPath + "/Inputs/DetectorConfiguration/" + Path; ifstream ConfigFile; @@ -293,6 +294,25 @@ void DetectorConstruction::ReadConfigurationFile(string Path) AddDetector(myDetector) ; } + //////////////////////////////////////////// + ///// Search for S1 Annular detector ////// + //////////////////////////////////////////// + else if (LineBuffer.compare(0, 2, "W1") == 0 && cW1 == false) { + cW1 = true ; + G4cout << "//////// W1 Square detector ////////" << G4endl << G4endl; + + // Instantiate the new array as a VDetector Object + VDetector* myDetector = new W1(); + + // Read Position of Telescope + ConfigFile.close(); + myDetector->ReadConfiguration(Path); + ConfigFile.open(Path.c_str()); + + // Add array to the VDetector Vector + AddDetector(myDetector); + } + //////////////////////////////////////////// //////// Search for MUST2 Array //////// //////////////////////////////////////////// diff --git a/NPSimulation/src/W1.cc b/NPSimulation/src/W1.cc new file mode 100644 index 0000000000000000000000000000000000000000..65858527a51d85ce0b2bed1b290223614ca804fd --- /dev/null +++ b/NPSimulation/src/W1.cc @@ -0,0 +1,776 @@ +/***************************************************************************** + * Copyright (C) 2009-2010 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : 12/01/11 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Define the W1 detector from Micron * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <sstream> +#include <string> +#include <cmath> + +// G4 Geometry headers +#include "G4Box.hh" +#include "G4Tubs.hh" + +// G4 various headers +#include "G4MaterialTable.hh" +#include "G4Element.hh" +#include "G4ElementTable.hh" +#include "G4VisAttributes.hh" +#include "G4Colour.hh" +#include "G4RotationMatrix.hh" +#include "G4Transform3D.hh" +#include "G4PVPlacement.hh" +#include "G4PVDivision.hh" + +// G4 sensitive +#include "G4SDManager.hh" +#include "G4MultiFunctionalDetector.hh" + +// NPTool headers +#include "GeneralScorers.hh" +#include "W1.hh" +#include "W1Scorers.hh" +#include "TW1Data.h" +#include "RootOutput.h" + +// CLHEP +#include "CLHEP/Random/RandGauss.h" + +using namespace std; +using namespace CLHEP; +using namespace W1SQUARE; + + + +W1::W1() + : m_Event(new TW1Data) +{ + InitializeMaterials(); +} + + + +W1::~W1() +{ + delete m_Event; + delete m_MaterialAluminium; + delete m_MaterialSilicon; + delete m_MaterialIron; + delete m_MaterialVacuum ; +} + + + +void W1::AddDetector(G4ThreeVector X1_Y1, G4ThreeVector X16_Y1, + G4ThreeVector X1_Y16, G4ThreeVector X16_Y16) +{ + m_DefinitionType.push_back(true); + + m_X1_Y1.push_back(X1_Y1); + m_X16_Y1.push_back(X16_Y1); + m_X1_Y16.push_back(X1_Y16); + m_X16_Y16.push_back(X16_Y16); + + m_R.push_back(0); + m_Theta.push_back(0); + m_Phi.push_back(0); + m_beta_u.push_back(0); + m_beta_v.push_back(0); + m_beta_w.push_back(0); +} + + + +void W1::AddDetector(G4double R, G4double Theta, G4double Phi, + G4double beta_u, G4double beta_v, G4double beta_w) +{ + G4ThreeVector empty = G4ThreeVector(0, 0, 0); + + m_DefinitionType.push_back(false); + + m_R.push_back(R); + m_Theta.push_back(Theta); + m_Phi.push_back(Phi); + m_beta_u.push_back(beta_u); + m_beta_v.push_back(beta_v); + m_beta_w.push_back(beta_w); + + m_X1_Y1.push_back(empty); + m_X16_Y1.push_back(empty); + m_X1_Y16.push_back(empty); + m_X16_Y16.push_back(empty); + +} + + + +void W1::VolumeMaker(G4int DetecNumber, + G4ThreeVector position, + G4RotationMatrix* rotation, + G4LogicalVolume* world) +{ + G4double NbrTelescopes = DetecNumber; + G4String DetectorNumber; + ostringstream Number; + Number << NbrTelescopes; + DetectorNumber = Number.str(); + + //////////////////////////////////////////////////////////////// + ////////////// Starting Volume Definition ////////////////////// + //////////////////////////////////////////////////////////////// + // Little trick to avoid warning in compilation: Use a PVPlacement "buffer". + // If don't you will have a Warning unused variable 'myPVP' + G4PVPlacement* PVPBuffer ; + G4String Name = "W1Square" + DetectorNumber; + + // Definition of the volume containing the sensitive detector + G4Box* solidW1 = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length); + G4LogicalVolume* logicW1 = new G4LogicalVolume(solidW1, m_MaterialVacuum, Name, 0, 0, 0); + + PVPBuffer = new G4PVPlacement(G4Transform3D(*rotation, position), logicW1, Name, world, false, 0); + + logicW1->SetVisAttributes(G4VisAttributes::Invisible); + if (m_non_sensitive_part_visiualisation) logicW1->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90))); + + // Aluminium dead layers + G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, AluStripFront_PosZ); + G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripBack_PosZ); + + G4Box* solidAluStrip = new G4Box("AluBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*AluStripThickness); +// G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0); + G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialVacuum, "logicAluStrip", 0, 0, 0); + + PVPBuffer = new G4PVPlacement(0, positionAluStripFront, logicAluStrip, Name + "_AluStripFront", logicW1, false, 0); + PVPBuffer = new G4PVPlacement(0, positionAluStripBack, logicAluStrip, Name + "_AluStripBack", logicW1, false, 0); + + logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible); + + // Silicon detector itself + G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ); + + G4Box* solidSilicon = new G4Box("solidSilicon", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconThickness); + G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, m_MaterialSilicon, "logicSilicon", 0, 0, 0); + + PVPBuffer = new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicW1, false, 0); + + // Set Silicon strip sensible + logicSilicon->SetSensitiveDetector(m_Scorer); + + ///Visualisation of Silicon Strip + G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)); + logicSilicon->SetVisAttributes(SiliconVisAtt); +} + + + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void W1::ReadConfiguration(string Path) +{ + ifstream ConfigFile; + ConfigFile.open(Path.c_str()); + string LineBuffer, DataBuffer; + + // A:X1_Y1 --> X:1 Y:1 + // B:X16_Y1 --> X:16 Y:1 + // C:X1_Y16 --> X:1 Y:16 + // D:X16_Y16 --> X:16 Y:16 + + G4double Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Az, Bz, Cz, Dz; + G4ThreeVector A, B, C, D; + G4double Theta = 0, Phi = 0, R = 0, beta_u = 0, beta_v = 0, beta_w = 0; + + bool ReadingStatus = false; + + bool check_A = false; + bool check_C = false; + bool check_B = false; + bool check_D = false; + bool check_R = false; + bool check_Theta = false; + bool check_Phi = false; + bool check_beta = false; + bool checkVis = false ; + + while (!ConfigFile.eof()) { + getline(ConfigFile, LineBuffer); + + if (LineBuffer.compare(0, 2, "W1") == 0) { + G4cout << "///" << G4endl; + G4cout << "W1 element found: " << G4endl; + ReadingStatus = true; + } + else ReadingStatus = false; + + while (ReadingStatus) { + ConfigFile >> DataBuffer; + + // Search for comment Symbol % + if (DataBuffer.compare(0, 1, "%") == 0) { + ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' ); + } + + // Position method + else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) { + check_A = true; + ConfigFile >> DataBuffer; + Ax = atof(DataBuffer.c_str()); + Ax = Ax * mm; + ConfigFile >> DataBuffer; + Ay = atof(DataBuffer.c_str()); + Ay = Ay * mm; + ConfigFile >> DataBuffer; + Az = atof(DataBuffer.c_str()); + Az = Az * mm; + + A = G4ThreeVector(Ax, Ay, Az); + cout << "X1 Y1 corner position : " << A << endl; + } + + else if (DataBuffer.compare(0, 7, "X16_Y1=") == 0) { + check_B = true; + ConfigFile >> DataBuffer; + Bx = atof(DataBuffer.c_str()); + Bx = Bx * mm; + ConfigFile >> DataBuffer; + By = atof(DataBuffer.c_str()); + By = By * mm; + ConfigFile >> DataBuffer; + Bz = atof(DataBuffer.c_str()); + Bz = Bz * mm ; + + B = G4ThreeVector(Bx, By, Bz); + cout << "X16 Y1 corner position : " << B << endl; + } + + else if (DataBuffer.compare(0, 7, "X1_Y16=") == 0) { + check_C = true; + ConfigFile >> DataBuffer; + Cx = atof(DataBuffer.c_str()); + Cx = Cx * mm; + ConfigFile >> DataBuffer; + Cy = atof(DataBuffer.c_str()); + Cy = Cy * mm; + ConfigFile >> DataBuffer; + Cz = atof(DataBuffer.c_str()); + Cz = Cz * mm; + + C = G4ThreeVector(Cx, Cy, Cz); + cout << "X1 Y16 corner position : " << C << endl; + } + + else if (DataBuffer.compare(0, 8, "X16_Y16=") == 0) { + check_D = true; + ConfigFile >> DataBuffer; + Dx = atof(DataBuffer.c_str()); + Dx = Dx * mm; + ConfigFile >> DataBuffer; + Dy = atof(DataBuffer.c_str()); + Dy = Dy * mm; + ConfigFile >> DataBuffer; + Dz = atof(DataBuffer.c_str()); + Dz = Dz * mm; + + D = G4ThreeVector(Dx, Dy, Dz); + cout << "X16 Y16 corner position : " << D << endl; + } + + + // Angle method + else if (DataBuffer.compare(0, 6, "THETA=") == 0) { + check_Theta = true; + ConfigFile >> DataBuffer ; + Theta = atof(DataBuffer.c_str()) ; + Theta = Theta * deg; + cout << "Theta: " << Theta / deg << endl; + } + + else if (DataBuffer.compare(0, 4, "PHI=") == 0) { + check_Phi = true; + ConfigFile >> DataBuffer ; + Phi = atof(DataBuffer.c_str()) ; + Phi = Phi * deg; + cout << "Phi: " << Phi / deg << endl; + } + + else if (DataBuffer.compare(0, 2, "R=") == 0) { + check_R = true; + ConfigFile >> DataBuffer ; + R = atof(DataBuffer.c_str()) ; + R = R * mm; + cout << "R: " << R / mm << endl; + } + + else if (DataBuffer.compare(0, 5, "BETA=") == 0) { + check_beta = true; + ConfigFile >> DataBuffer; + beta_u = atof(DataBuffer.c_str()); + beta_u = beta_u * deg; + ConfigFile >> DataBuffer; + beta_v = atof(DataBuffer.c_str()); + beta_v = beta_v * deg; + ConfigFile >> DataBuffer; + beta_w = atof(DataBuffer.c_str()); + beta_w = beta_w * deg; + G4cout << "Beta: " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl ; + } + + else if (DataBuffer.compare(0, 4, "VIS=") == 0) { + checkVis = true; + ConfigFile >> DataBuffer; + if (DataBuffer.compare(0, 3, "all") == 0) m_non_sensitive_part_visiualisation = true; + } + + else { + /////////////////////////////////////////////////// + // If no Detector Token and no comment, toggle out + ReadingStatus = false; + G4cout << "Wrong Token Sequence: Getting out " << DataBuffer << G4endl; + } + + // Add The previously define detector + // With position method + if ((check_A && check_B && check_C && check_D && checkVis) && !(check_Theta && check_Phi && check_R)) { + ReadingStatus = false; + check_A = false; + check_C = false; + check_B = false; + check_D = false; + checkVis = false; + + AddDetector(A, B, C, D); + } + + // with angle method + if ((check_Theta && check_Phi && check_R && checkVis) && !(check_A && check_B && check_C && check_D)) { + ReadingStatus = false; + check_R = false; + check_Theta = false; + check_Phi = false; + check_beta = false; + checkVis = false; + + AddDetector(R, Theta, Phi, beta_u, beta_v, beta_w); + } + } + } +} + + + +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void W1::ConstructDetector(G4LogicalVolume* world) +{ + G4RotationMatrix* W1rot = NULL; + G4ThreeVector W1pos = G4ThreeVector(0, 0, 0) ; + G4ThreeVector W1u = G4ThreeVector(0, 0, 0) ; + G4ThreeVector W1v = G4ThreeVector(0, 0, 0) ; + G4ThreeVector W1w = G4ThreeVector(0, 0, 0) ; + G4ThreeVector W1Center = G4ThreeVector(0, 0, 0) ; + + G4int NumberOfDetector = m_DefinitionType.size() ; + for (G4int i = 0; i < NumberOfDetector; i++) { + // By Point + if (m_DefinitionType[i]) { + // (u,v,w) unitary vector associated to telescope referencial + // (u,v) // to silicon plan + // w perpendicular to (u,v) plan and pointing ThirdStage + W1u = m_X16_Y1[i] - m_X1_Y1[i]; + W1u = W1u.unit(); + + W1v = m_X1_Y16[i] - m_X1_Y1[i]; + W1v = W1v.unit(); + + W1w = W1u.cross(W1v); + W1w = W1w.unit(); + + W1Center = (m_X1_Y1[i] + m_X1_Y16[i] + m_X16_Y1[i] + m_X16_Y16[i]) / 4; + + // Passage Matrix from Lab Referential to Telescope Referential + W1rot = new G4RotationMatrix(W1u, W1v, W1w); + // translation to place Telescope + W1pos = W1w * Length * 0.5 + W1Center; + } + + // By Angle + else { + G4double Theta = m_Theta[i]; + G4double Phi = m_Phi[i]; + + // (u,v,w) unitary vector associated to telescope referencial + // (u,v) // to silicon plan + // w perpendicular to (u,v) plan and pointing ThirdStage + // Phi is angle between X axis and projection in (X,Y) plan + // Theta is angle between position vector and z axis + G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad); + G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad); + G4double wZ = m_R[i] * cos(Theta / rad); + W1w = G4ThreeVector(wX, wY, wZ); + + // vector corresponding to the center of the module + W1Center = W1w; + + // vector parallel to one axis of silicon plane + G4double ii = cos(Theta / rad) * cos(Phi / rad); + G4double jj = cos(Theta / rad) * sin(Phi / rad); + G4double kk = -sin(Theta / rad); + G4ThreeVector Y = G4ThreeVector(ii, jj, kk); + + W1w = W1w.unit(); + W1u = W1w.cross(Y); + W1v = W1w.cross(W1u); + W1v = W1v.unit(); + W1u = W1u.unit(); + + // Passage Matrix from Lab Referential to Telescope Referential + // MUST2 + W1rot = new G4RotationMatrix(W1u, W1v, W1w); + // Telescope is rotate of Beta angle around W1v axis. + W1rot->rotate(m_beta_u[i], W1u); + W1rot->rotate(m_beta_v[i], W1v); + W1rot->rotate(m_beta_w[i], W1w); + // translation to place Telescope + W1pos = W1w * Length * 0.5 + W1Center; + } + + VolumeMaker(i + 1, W1pos, W1rot, world); + } + + delete W1rot; +} + + + +// Connect the GaspardTrackingData class to the output TTree +// of the simulation +void W1::InitializeRootOutput() +{ + RootOutput *pAnalysis = RootOutput::getInstance(); + TTree *pTree = pAnalysis->GetTree(); + pTree->Branch("W1", "TW1Data", &m_Event); +} + + + +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAvtion +void W1::ReadSensitive(const G4Event* event) +{ + // Clear ROOT objects + m_Event->Clear(); + + ////////////////////////////////////////////////////////////////////////////////// + /////////////// Variables Used to Read Event Map of detector //////////////////// + ////////////////////////////////////////////////////////////////////////////////// + // First Stage + std::map<G4int, G4int*>::iterator DetectorNumber_itr; + std::map<G4int, G4double*>::iterator Energy_itr; + std::map<G4int, G4double*>::iterator Time_itr; + std::map<G4int, G4double*>::iterator X_itr; + std::map<G4int, G4double*>::iterator Y_itr; + std::map<G4int, G4double*>::iterator Pos_X_itr; + std::map<G4int, G4double*>::iterator Pos_Y_itr; + std::map<G4int, G4double*>::iterator Pos_Z_itr; + std::map<G4int, G4double*>::iterator Ang_Theta_itr; + std::map<G4int, G4double*>::iterator Ang_Phi_itr; + + G4THitsMap<G4int>* DetectorNumberHitMap; + G4THitsMap<G4double>* EnergyHitMap; + G4THitsMap<G4double>* TimeHitMap; + G4THitsMap<G4double>* XHitMap; + G4THitsMap<G4double>* YHitMap; + G4THitsMap<G4double>* PosXHitMap; + G4THitsMap<G4double>* PosYHitMap; + G4THitsMap<G4double>* PosZHitMap; + G4THitsMap<G4double>* AngThetaHitMap; + G4THitsMap<G4double>* AngPhiHitMap; + + ///////////////////////////////////////////////////////////////////////////////// + ////////////////// Read scorers //////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////////// + // Read the Scorer associated to the first Stage + //Detector Number + G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/DetectorNumber"); + DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)); + DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin(); + + //Energy + G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/StripEnergy") ; + EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ; + Energy_itr = EnergyHitMap->GetMap()->begin() ; + + //Time of Flight + G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/StripTime") ; + TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ; + Time_itr = TimeHitMap->GetMap()->begin() ; + + //Strip Number X + G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/FrontStripNumber"); + XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)); + X_itr = XHitMap->GetMap()->begin(); + + //Strip Number Y + G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/BackStripNumber"); + YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ; + Y_itr = YHitMap->GetMap()->begin() ; + + //Interaction Coordinate X + G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/InterCoordX") ; + PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ; + Pos_X_itr = PosXHitMap->GetMap()->begin() ; + + //Interaction Coordinate Y + G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/InterCoordY") ; + PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ; + Pos_Y_itr = PosYHitMap->GetMap()->begin() ; + + //Interaction Coordinate Z + G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/InterCoordZ") ; + PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ; + Pos_Z_itr = PosXHitMap->GetMap()->begin() ; + + //Interaction Coordinate Angle Theta + G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/InterCoordAngTheta") ; + AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ; + Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ; + + //Interaction Coordinate Angle Phi + G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ScorerW1/InterCoordAngPhi"); + AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)); + Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); + + // Check the size of different map + G4int sizeN = DetectorNumberHitMap->entries(); + G4int sizeE = EnergyHitMap->entries(); + G4int sizeT = TimeHitMap->entries(); + G4int sizeX = XHitMap->entries(); + G4int sizeY = YHitMap->entries(); + + if (sizeE != sizeT || sizeT != sizeX || sizeX != sizeY) { + G4cout << "No match size Si Event Map: sE:" + << sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << endl ; + return; + } + + // Loop on W1 number + for (G4int l = 0; l < sizeN; l++) { + G4double N = *(DetectorNumber_itr->second); + G4int NTrackID = DetectorNumber_itr->first - N; + + if (N > 0) { + // Fill detector number + m_Event->SetW1FrontEDetectorNbr(N); + m_Event->SetW1FrontTDetectorNbr(N); + m_Event->SetW1BackEDetectorNbr(N); + m_Event->SetW1BackTDetectorNbr(N); + + // Energy + Energy_itr = EnergyHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeE ; h++) { + G4int ETrackID = Energy_itr->first - N; + G4double E = *(Energy_itr->second); + if (ETrackID == NTrackID) { + m_Event->SetW1FrontEEnergy(RandGauss::shoot(E, EnergyResolution)); + m_Event->SetW1BackEEnergy(RandGauss::shoot(E, EnergyResolution)); + } + Energy_itr++; + } + + // Time + Time_itr = TimeHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeT ; h++) { + G4int TTrackID = Time_itr->first - N; + G4double T = *(Time_itr->second); + if (TTrackID == NTrackID) { + m_Event->SetW1FrontTTime(RandGauss::shoot(T, TimeResolution)); + m_Event->SetW1BackTTime(RandGauss::shoot(T, TimeResolution)); + } + Time_itr++; + } + + // strip front + X_itr = XHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeX ; h++) { + G4int XTrackID = X_itr->first - N; + G4double X = *(X_itr->second); + if (XTrackID == NTrackID) { + m_Event->SetW1FrontEStripNbr(X) ; + m_Event->SetW1FrontTStripNbr(X) ; + } + X_itr++; + } + + // strip back + Y_itr = YHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeY ; h++) { + G4int YTrackID = Y_itr->first - N; + G4double Y = *(Y_itr->second); + if (YTrackID == NTrackID) { + m_Event->SetW1BackEStripNbr(Y); + m_Event->SetW1BackTStripNbr(Y); + } + Y_itr++; + } + + // Pos X + Pos_X_itr = PosXHitMap->GetMap()->begin(); + for (G4int h = 0; h < PosXHitMap->entries(); h++) { + G4int PosXTrackID = Pos_X_itr->first - N; + G4double PosX = *(Pos_X_itr->second); + if (PosXTrackID == NTrackID) { + ms_InterCoord->SetDetectedPositionX(PosX); + } + Pos_X_itr++; + } + + // Pos Y + Pos_Y_itr = PosYHitMap->GetMap()->begin(); + for (G4int h = 0; h < PosYHitMap->entries(); h++) { + G4int PosYTrackID = Pos_Y_itr->first - N; + G4double PosY = *(Pos_Y_itr->second); + if (PosYTrackID == NTrackID) { + ms_InterCoord->SetDetectedPositionY(PosY); + } + Pos_Y_itr++; + } + + // Pos Z + Pos_Z_itr = PosZHitMap->GetMap()->begin(); + for (G4int h = 0; h < PosZHitMap->entries(); h++) { + G4int PosZTrackID = Pos_Z_itr->first - N; + G4double PosZ = *(Pos_Z_itr->second); + if (PosZTrackID == NTrackID) { + ms_InterCoord->SetDetectedPositionZ(PosZ); + } + Pos_Z_itr++; + } + + // Angle Theta + Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); + for (G4int h = 0; h < AngThetaHitMap->entries(); h++) { + G4int AngThetaTrackID = Ang_Theta_itr->first - N; + G4double AngTheta = *(Ang_Theta_itr->second); + if (AngThetaTrackID == NTrackID) { + ms_InterCoord->SetDetectedAngleTheta(AngTheta); + } + Ang_Theta_itr++; + } + + // Angle Phi + Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); + for (G4int h = 0; h < AngPhiHitMap->entries(); h++) { + G4int AngPhiTrackID = Ang_Phi_itr->first - N; + G4double AngPhi = *(Ang_Phi_itr->second); + if (AngPhiTrackID == NTrackID) { + ms_InterCoord->SetDetectedAnglePhi(AngPhi); + } + Ang_Phi_itr++; + } + + } // end if number of detector > 0 + + DetectorNumber_itr++; + } // end loop on detector multiplicity + + // clear map for next event + DetectorNumberHitMap -> clear(); + EnergyHitMap -> clear(); + TimeHitMap -> clear(); + XHitMap -> clear(); + YHitMap -> clear(); + PosXHitMap -> clear(); + PosYHitMap -> clear(); + PosZHitMap -> clear(); + AngThetaHitMap -> clear(); + AngPhiHitMap -> clear(); +} + + + +void W1::InitializeMaterials() +{ + //////////////////////////////////////////////////////////////// + /////////////////Element Definition /////////////////////////// + //////////////////////////////////////////////////////////////// + G4String symbol; + G4double density = 0, a = 0, z = 0; + G4int ncomponents = 0; + + G4Element* N = new G4Element("Nitrogen" , symbol = "N" , z = 7 , a = 14.01 * g / mole); + G4Element* O = new G4Element("Oxigen" , symbol = "O" , z = 8 , a = 16.00 * g / mole); + + //////////////////////////////////////////////////////////////// + /////////////////Material Definition /////////////////////////// + //////////////////////////////////////////////////////////////// + // Si + a = 28.0855 * g / mole; + density = 2.321 * g / cm3; + m_MaterialSilicon = new G4Material("Si", z = 14., a, density); + + // Al + density = 2.702 * g / cm3; + a = 26.98 * g / mole; + m_MaterialAluminium = new G4Material("Aluminium", z = 13., a, density); + + // Iron + density = 7.874 * g / cm3; + a = 55.847 * g / mole; + m_MaterialIron = new G4Material("Iron", z = 26., a, density); + + // Vacuum + density = 0.000000001 * mg / cm3; + m_MaterialVacuum = new G4Material("Vacuum", density, ncomponents = 2); + m_MaterialVacuum->AddElement(N, .7); + m_MaterialVacuum->AddElement(O, .3); +} + + + +void W1::InitializeScorers() +{ + // Associate Scorer + m_Scorer = new G4MultiFunctionalDetector("ScorerW1"); + G4VPrimitiveScorer* DetNbr = new GENERALSCORERS::PSDetectorNumber("DetectorNumber", "W1Square", 0); + G4VPrimitiveScorer* Energy = new GENERALSCORERS::PSEnergy("StripEnergy", "W1Square", 0); + G4VPrimitiveScorer* TOF = new GENERALSCORERS::PSTOF("StripTime", "W1Square", 0); + G4VPrimitiveScorer* InteractionCoordinatesX = new GENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","W1Square", 0); + G4VPrimitiveScorer* InteractionCoordinatesY = new GENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","W1Square", 0); + G4VPrimitiveScorer* InteractionCoordinatesZ = new GENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","W1Square", 0); + G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new GENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","W1Square", 0); + G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new GENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","W1Square", 0); + G4VPrimitiveScorer* ThetaStripPosition = new W1ScorerFrontStripNumber("FrontStripNumber", 0, NbStrips); + G4VPrimitiveScorer* PhiStripPosition = new W1ScorerBackStripNumber("BackStripNumber", 0, NbStrips); + + //and register it to the multifunctionnal detector + m_Scorer->RegisterPrimitive(DetNbr); + m_Scorer->RegisterPrimitive(Energy); + m_Scorer->RegisterPrimitive(TOF); + m_Scorer->RegisterPrimitive(ThetaStripPosition); + m_Scorer->RegisterPrimitive(PhiStripPosition); + m_Scorer->RegisterPrimitive(InteractionCoordinatesX); + m_Scorer->RegisterPrimitive(InteractionCoordinatesY); + m_Scorer->RegisterPrimitive(InteractionCoordinatesZ); + m_Scorer->RegisterPrimitive(InteractionCoordinatesAngleTheta); + m_Scorer->RegisterPrimitive(InteractionCoordinatesAnglePhi); + + // Add All Scorer to the Global Scorer Manager + G4SDManager::GetSDMpointer()->AddNewDetector(m_Scorer); +} diff --git a/NPSimulation/src/W1Scorers.cc b/NPSimulation/src/W1Scorers.cc new file mode 100644 index 0000000000000000000000000000000000000000..5d4d0fd13236996d660cfdba5f9b77db11b1033b --- /dev/null +++ b/NPSimulation/src/W1Scorers.cc @@ -0,0 +1,172 @@ +/***************************************************************************** + * Copyright (C) 2009-2010 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr * + * * + * Creation Date : 12/01/11 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class holds all the scorers needed by the W1 object * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// G4 headers +#include "G4UnitsTable.hh" + +// NPTool headers +#include "GeneralScorers.hh" +#include "W1Scorers.hh" + +#include "W1.hh" +using namespace W1SQUARE; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Added by Adrien MATTA: +// Those Scorer use TrackID as map index. This way ones can rebuild energy deposit, +// time of flight or position,... particle by particle for each event. Because standard +// scorer provide by G4 don't work this way but using a global ID for each event you should +// not use those scorer with some G4 provided ones or being very carefull doing so. + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Front Strip position Scorer +W1ScorerFrontStripNumber::W1ScorerFrontStripNumber(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip ; +} + +W1ScorerFrontStripNumber::~W1ScorerFrontStripNumber() +{ +} + +G4bool W1ScorerFrontStripNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "W1Square"); + + // get front strip number + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); +// G4cout << "POS world: " << POS << G4endl; + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); +// G4cout << "POS local: " << POS << G4endl; + + G4double StripPitch = W1SQUARE::SiliconFace / m_NumberOfStrip; + + G4double temp = (POS(0) + W1SQUARE::SiliconFace / 2.) / StripPitch; + G4double X = int(temp) + 1 ; +// G4cout << "strip X: " << X << G4endl; + + //Rare case where particle is close to edge of silicon plan + if (X == m_NumberOfStrip+1) X = m_NumberOfStrip; + G4double edep = aStep->GetTotalEnergyDeposit(); + if (edep < 100*keV) return FALSE; + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void W1ScorerFrontStripNumber::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void W1ScorerFrontStripNumber::EndOfEvent(G4HCofThisEvent*) +{ +} + +void W1ScorerFrontStripNumber::Clear() +{ + EvtMap->clear(); +} + +void W1ScorerFrontStripNumber::DrawAll() +{ +} + +void W1ScorerFrontStripNumber::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Back Strip position Scorer +W1ScorerBackStripNumber::W1ScorerBackStripNumber(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip ; +} + +W1ScorerBackStripNumber::~W1ScorerBackStripNumber() +{ +} + +G4bool W1ScorerBackStripNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "W1Square"); + + // get back strip number + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = W1SQUARE::SiliconFace / m_NumberOfStrip; + + G4double temp = (POS(1) + W1SQUARE::SiliconFace / 2.) / StripPitch ; + G4double X = int(temp) + 1 ; +// G4cout << "strip Y: " << X << G4endl; + + //Rare case where particle is close to edge of silicon plan + if (X == m_NumberOfStrip+1) X = m_NumberOfStrip; + G4double edep = aStep->GetTotalEnergyDeposit(); + if (edep < 100*keV) return FALSE; + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void W1ScorerBackStripNumber::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void W1ScorerBackStripNumber::EndOfEvent(G4HCofThisEvent*) +{ +} + +void W1ScorerBackStripNumber::Clear() +{ + EvtMap->clear(); +} + +void W1ScorerBackStripNumber::DrawAll() +{ +} + +void W1ScorerBackStripNumber::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +}