From fc1a8cf4d09161306958ed5a3ff4f8e0a8cec447 Mon Sep 17 00:00:00 2001 From: matta <matta@npt> Date: Thu, 7 Feb 2013 16:29:25 +0000 Subject: [PATCH] * Porting Helios form Marc Labiche to new NPTool version --- Inputs/DetectorConfiguration/helios.detector | 426 +++++++++ NPLib/Helios/Makefile | 2 +- NPSimulation/include/Helios.hh | 73 ++ NPSimulation/include/HeliosDetDummyShape.hh | 170 ++++ NPSimulation/include/HeliosModule.hh | 95 ++ NPSimulation/include/HeliosScorers.hh | 249 ++++++ NPSimulation/include/MyMagneticField.hh | 50 ++ NPSimulation/src/DetectorConstruction.cc | 128 ++- NPSimulation/src/Helios.cc | 194 ++++ NPSimulation/src/HeliosDetDummyShape.cc | 886 +++++++++++++++++++ NPSimulation/src/HeliosModule.cc | 63 ++ NPSimulation/src/HeliosScorers.cc | 642 ++++++++++++++ NPSimulation/src/MyMagneticField.cc | 71 ++ 13 files changed, 3047 insertions(+), 2 deletions(-) create mode 100644 Inputs/DetectorConfiguration/helios.detector create mode 100644 NPSimulation/include/Helios.hh create mode 100644 NPSimulation/include/HeliosDetDummyShape.hh create mode 100644 NPSimulation/include/HeliosModule.hh create mode 100644 NPSimulation/include/HeliosScorers.hh create mode 100644 NPSimulation/include/MyMagneticField.hh create mode 100644 NPSimulation/src/Helios.cc create mode 100644 NPSimulation/src/HeliosDetDummyShape.cc create mode 100644 NPSimulation/src/HeliosModule.cc create mode 100644 NPSimulation/src/HeliosScorers.cc create mode 100644 NPSimulation/src/MyMagneticField.cc diff --git a/Inputs/DetectorConfiguration/helios.detector b/Inputs/DetectorConfiguration/helios.detector new file mode 100644 index 000000000..1337a3707 --- /dev/null +++ b/Inputs/DetectorConfiguration/helios.detector @@ -0,0 +1,426 @@ +%%%%%%%%%Target%%%%%%%%%%%%%%%%%%%%%% +%% +%% Thickness in micrometer +%% Radius in mm +%% Position in mm +%% Magnetic field in tesla +%% +%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%% +%% +%% Position and distance given in mm +%% Angle given in degree +%% using the data from the experimental mesurement +%% special care is given for the X Y direction +%% NOTATTION USED IN THE FOLLOWING: +%% +%% X1_Y1 --> X:1 Y:1 +%% X128_Y1 --> X:128 Y:1 +%% X1_Y128 --> X:1 Y:128 +%% X128_Y128 --> X:128 Y:128 +%%Option: 0,1 for Si SiLi and CSI +%%Option: all or sensible for VISualisation +% 1.2 12.371134021 +% 2 20.618556701 +% 30 309.278350515 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +GeneralTarget +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +Target + THICKNESS= 1. +% THICKNESS= 0.28 + RADIUS= 7.5 + MATERIAL= CD2 +% MATERIAL= Ti-t + X= 0 + Y= 0 + Z= 0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Helios + MField= 2. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 +HeliosDummyShape + X1_Y1= 11.5 11.5 150. + X128_Y1= 11.5 -11.5 150. + X1_Y128= 11.5 11.5 206. + X128_Y128= 11.5 -11.5 206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2 +HeliosDummyShape + X1_Y1= -11.5 11.5 150. + X128_Y1= 11.5 11.5 150. + X1_Y128= -11.5 11.5 206. + X128_Y128= 11.5 11.5 206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3 +HeliosDummyShape + X1_Y1= -11.5 -11.5 150. + X128_Y1= -11.5 11.5 150. + X1_Y128= -11.5 -11.5 206. + X128_Y128= -11.5 11.5 206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%4 +HeliosDummyShape + X1_Y1= 11.5 -11.5 150. + X128_Y1= -11.5 -11.5 150. + X1_Y128= 11.5 -11.5 206. + X128_Y128= -11.5 -11.5 206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 +HeliosDummyShape + X1_Y1= 11.5 11.5 208.4 + X128_Y1= 11.5 -11.5 208.4 + X1_Y128= 11.5 11.5 264.4 + X128_Y128= 11.5 -11.5 264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%6 +HeliosDummyShape + X1_Y1= -11.5 11.5 208.4 + X128_Y1= 11.5 11.5 208.4 + X1_Y128= -11.5 11.5 264.4 + X128_Y128= 11.5 11.5 264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%7 +HeliosDummyShape + X1_Y1= -11.5 -11.5 208.4 + X128_Y1= -11.5 11.5 208.4 + X1_Y128= -11.5 -11.5 264.4 + X128_Y128= -11.5 11.5 264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%8 +HeliosDummyShape + X1_Y1= 11.5 -11.5 208.4 + X128_Y1= -11.5 -11.5 208.4 + X1_Y128= 11.5 -11.5 264.4 + X128_Y128= -11.5 -11.5 264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9 +HeliosDummyShape + X1_Y1= 11.5 11.5 266.8 + X128_Y1= 11.5 -11.5 266.8 + X1_Y128= 11.5 11.5 322.8 + X128_Y128= 11.5 -11.5 322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%10 +HeliosDummyShape + X1_Y1= -11.5 11.5 266.8 + X128_Y1= 11.5 11.5 266.8 + X1_Y128= -11.5 11.5 322.8 + X128_Y128= 11.5 11.5 322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%11 +HeliosDummyShape + X1_Y1= -11.5 -11.5 266.8 + X128_Y1= -11.5 11.5 266.8 + X1_Y128= -11.5 -11.5 322.8 + X128_Y128= -11.5 11.5 322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%12 +HeliosDummyShape + X1_Y1= 11.5 -11.5 266.8 + X128_Y1= -11.5 -11.5 266.8 + X1_Y128= 11.5 -11.5 322.8 + X128_Y128= -11.5 -11.5 322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%13 +HeliosDummyShape + X1_Y1= 11.5 11.5 325.2 + X128_Y1= 11.5 -11.5 325.2 + X1_Y128= 11.5 11.5 381.2 + X128_Y128= 11.5 -11.5 381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%14 +HeliosDummyShape + X1_Y1= -11.5 11.5 325.2 + X128_Y1= 11.5 11.5 325.2 + X1_Y128= -11.5 11.5 381.2 + X128_Y128= 11.5 11.5 381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%15 +HeliosDummyShape + X1_Y1= -11.5 -11.5 325.2 + X128_Y1= -11.5 11.5 325.2 + X1_Y128= -11.5 -11.5 381.2 + X128_Y128= -11.5 11.5 381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%16 +HeliosDummyShape + X1_Y1= 11.5 -11.5 325.2 + X128_Y1= -11.5 -11.5 325.2 + X1_Y128= 11.5 -11.5 381.2 + X128_Y128= -11.5 -11.5 381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%17 +HeliosDummyShape + X1_Y1= 11.5 11.5 383.6 + X128_Y1= 11.5 -11.5 383.6 + X1_Y128= 11.5 11.5 439.6 + X128_Y128= 11.5 -11.5 439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%18 +HeliosDummyShape + X1_Y1= -11.5 11.5 383.6 + X128_Y1= 11.5 11.5 383.6 + X1_Y128= -11.5 11.5 439.6 + X128_Y128= 11.5 11.5 439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%19 +HeliosDummyShape + X1_Y1= -11.5 -11.5 383.6 + X128_Y1= -11.5 11.5 383.6 + X1_Y128= -11.5 -11.5 439.6 + X128_Y128= -11.5 11.5 439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%20 +HeliosDummyShape + X1_Y1= 11.5 -11.5 383.6 + X128_Y1= -11.5 -11.5 383.6 + X1_Y128= 11.5 -11.5 439.6 + X128_Y128= -11.5 -11.5 439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%21 +HeliosDummyShape + X1_Y1= 11.5 11.5 442. + X128_Y1= 11.5 -11.5 442. + X1_Y128= 11.5 11.5 498. + X128_Y128= 11.5 -11.5 498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%22 +HeliosDummyShape + X1_Y1= -11.5 11.5 442. + X128_Y1= 11.5 11.5 442. + X1_Y128= -11.5 11.5 498. + X128_Y128= 11.5 11.5 498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%23 +HeliosDummyShape + X1_Y1= -11.5 -11.5 442. + X128_Y1= -11.5 11.5 442. + X1_Y128= -11.5 -11.5 498. + X128_Y128= -11.5 11.5 498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%24 +HeliosDummyShape + X1_Y1= 11.5 -11.5 442. + X128_Y1= -11.5 -11.5 442. + X1_Y128= 11.5 -11.5 498 + X128_Y128= -11.5 -11.5 498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%25 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% AT backward angles: +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%25 +HeliosDummyShape + X1_Y1= 11.5 11.5 -150. + X128_Y1= 11.5 -11.5 -150. + X1_Y128= 11.5 11.5 -206. + X128_Y128= 11.5 -11.5 -206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%26 +HeliosDummyShape + X1_Y1= -11.5 11.5 -150. + X128_Y1= 11.5 11.5 -150. + X1_Y128= -11.5 11.5 -206. + X128_Y128= 11.5 11.5 -206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%27 +HeliosDummyShape + X1_Y1= -11.5 -11.5 -150. + X128_Y1= -11.5 11.5 -150. + X1_Y128= -11.5 -11.5 -206. + X128_Y128= -11.5 11.5 -206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%28 +HeliosDummyShape + X1_Y1= 11.5 -11.5 -150. + X128_Y1= -11.5 -11.5 -150. + X1_Y128= 11.5 -11.5 -206. + X128_Y128= -11.5 -11.5 -206. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%29 +HeliosDummyShape + X1_Y1= 11.5 11.5 -208.4 + X128_Y1= 11.5 -11.5 -208.4 + X1_Y128= 11.5 11.5 -264.4 + X128_Y128= 11.5 -11.5 -264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%30 +HeliosDummyShape + X1_Y1= -11.5 11.5 -208.4 + X128_Y1= 11.5 11.5 -208.4 + X1_Y128= -11.5 11.5 -264.4 + X128_Y128= 11.5 11.5 -264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%31 +HeliosDummyShape + X1_Y1= -11.5 -11.5 -208.4 + X128_Y1= -11.5 11.5 -208.4 + X1_Y128= -11.5 -11.5 -264.4 + X128_Y128= -11.5 11.5 -264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%32 +HeliosDummyShape + X1_Y1= 11.5 -11.5 -208.4 + X128_Y1= -11.5 -11.5 -208.4 + X1_Y128= 11.5 -11.5 -264.4 + X128_Y128= -11.5 -11.5 -264.4 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%33 +HeliosDummyShape + X1_Y1= 11.5 11.5 -266.8 + X128_Y1= 11.5 -11.5 -266.8 + X1_Y128= 11.5 11.5 -322.8 + X128_Y128= 11.5 -11.5 -322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%34 +HeliosDummyShape + X1_Y1= -11.5 11.5 -266.8 + X128_Y1= 11.5 11.5 -266.8 + X1_Y128= -11.5 11.5 -322.8 + X128_Y128= 11.5 11.5 -322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%35 +HeliosDummyShape + X1_Y1= -11.5 -11.5 -266.8 + X128_Y1= -11.5 11.5 -266.8 + X1_Y128= -11.5 -11.5 -322.8 + X128_Y128= -11.5 11.5 -322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%36 +HeliosDummyShape + X1_Y1= 11.5 -11.5 -266.8 + X128_Y1= -11.5 -11.5 -266.8 + X1_Y128= 11.5 -11.5 -322.8 + X128_Y128= -11.5 -11.5 -322.8 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%37 +HeliosDummyShape + X1_Y1= 11.5 11.5 -325.2 + X128_Y1= 11.5 -11.5 -325.2 + X1_Y128= 11.5 11.5 -381.2 + X128_Y128= 11.5 -11.5 -381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%38 +HeliosDummyShape + X1_Y1= -11.5 11.5 -325.2 + X128_Y1= 11.5 11.5 -325.2 + X1_Y128= -11.5 11.5 -381.2 + X128_Y128= 11.5 11.5 -381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%39 +HeliosDummyShape + X1_Y1= -11.5 -11.5 -325.2 + X128_Y1= -11.5 11.5 -325.2 + X1_Y128= -11.5 -11.5 -381.2 + X128_Y128= -11.5 11.5 -381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%40 +HeliosDummyShape + X1_Y1= 11.5 -11.5 -325.2 + X128_Y1= -11.5 -11.5 -325.2 + X1_Y128= 11.5 -11.5 -381.2 + X128_Y128= -11.5 -11.5 -381.2 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%41 +HeliosDummyShape + X1_Y1= 11.5 11.5 -383.6 + X128_Y1= 11.5 -11.5 -383.6 + X1_Y128= 11.5 11.5 -439.6 + X128_Y128= 11.5 -11.5 -439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%42 +HeliosDummyShape + X1_Y1= -11.5 11.5 -383.6 + X128_Y1= 11.5 11.5 -383.6 + X1_Y128= -11.5 11.5 -439.6 + X128_Y128= 11.5 11.5 -439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%43 +HeliosDummyShape + X1_Y1= -11.5 -11.5 -383.6 + X128_Y1= -11.5 11.5 -383.6 + X1_Y128= -11.5 -11.5 -439.6 + X128_Y128= -11.5 11.5 -439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%44 +HeliosDummyShape + X1_Y1= 11.5 -11.5 -383.6 + X128_Y1= -11.5 -11.5 -383.6 + X1_Y128= 11.5 -11.5 -439.6 + X128_Y128= -11.5 -11.5 -439.6 + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%45 +HeliosDummyShape + X1_Y1= 11.5 11.5 -442. + X128_Y1= 11.5 -11.5 -442. + X1_Y128= 11.5 11.5 -498. + X128_Y128= 11.5 -11.5 -498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%46 +HeliosDummyShape + X1_Y1= -11.5 11.5 -442. + X128_Y1= 11.5 11.5 -442. + X1_Y128= -11.5 11.5 -498. + X128_Y128= 11.5 11.5 -498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%47 +HeliosDummyShape + X1_Y1= -11.5 -11.5 -442. + X128_Y1= -11.5 11.5 -442. + X1_Y128= -11.5 -11.5 -498. + X128_Y128= -11.5 11.5 -498. + FIRSTSTAGE= 1 + VIS= all +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%48 +HeliosDummyShape + X1_Y1= 11.5 -11.5 -442. + X128_Y1= -11.5 -11.5 -442. + X1_Y128= 11.5 -11.5 -498 + X128_Y128= -11.5 -11.5 -498. + FIRSTSTAGE= 1 + VIS= all diff --git a/NPLib/Helios/Makefile b/NPLib/Helios/Makefile index a9645448b..1cfbe7db9 100644 --- a/NPLib/Helios/Makefile +++ b/NPLib/Helios/Makefile @@ -6,7 +6,7 @@ all: $(SHARELIB) #------------------------------------------------------------------------------ ############### Detector ############## -## MUST2 ## +## Helios ## libHelios.so: THeliosData.o THeliosDataDict.o Helios.o THeliosPhysics.o THeliosPhysicsDict.o $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ diff --git a/NPSimulation/include/Helios.hh b/NPSimulation/include/Helios.hh new file mode 100644 index 000000000..d0f0f1a7f --- /dev/null +++ b/NPSimulation/include/Helios.hh @@ -0,0 +1,73 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class manages different shapes of module for the Helios * + * detector. It allows to have Helios geometries with an * + * heterogeneous set of modules * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef Helios_h +#define Helios_h 1 + +// C++ headers +#include <vector> + +// NPTool header +#include "VDetector.hh" +#include "HeliosModule.hh" + +using namespace std; + + + +class Helios : public VDetector +{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// +public: + Helios(); + virtual ~Helios(); + + //////////////////////////////////////////////////// + ///////// 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(); + + // Read sensitive part and fill the Root tree. + // Called at in the EventAction::EndOfEventAvtion + void ReadSensitive(const G4Event* event); + +public: + // Initialize all scorers necessary for each detector + void InitializeScorers(); + +private: + vector<HeliosModule*> m_Modules; +}; +#endif diff --git a/NPSimulation/include/HeliosDetDummyShape.hh b/NPSimulation/include/HeliosDetDummyShape.hh new file mode 100644 index 000000000..3f2fa2332 --- /dev/null +++ b/NPSimulation/include/HeliosDetDummyShape.hh @@ -0,0 +1,170 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Define a dummy module for the Helios detector * + * The goal of this class is to be a starting point to create a * + * new shape to be added to the Helios detector. * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef HeliosDetDummyShape_h +#define HeliosDetDummyShape_h 1 + +// C++ headers +#include <vector> + +// NPTool header +#include "HeliosModule.hh" +#include "TInteractionCoordinates.h" + +using namespace std; + + + +class HeliosDetDummyShape : public HeliosModule +{ + //////////////////////////////////////////////////// + /////// Default Constructor and Destructor ///////// + //////////////////////////////////////////////////// +public: + HeliosDetDummyShape(); + virtual ~HeliosDetDummyShape(); + + //////////////////////////////////////////////////// + //////// Specific Function of this Class /////////// + //////////////////////////////////////////////////// +public: + // By Position Method + void AddModule(G4ThreeVector TL , + G4ThreeVector BL , + G4ThreeVector BR , + G4ThreeVector CT , + bool wFirstStage ); + + // By Angle Method + void AddModule(G4double R , + G4double Theta , + G4double Phi , + G4double beta_u , + G4double beta_v , + G4double beta_w , + bool wFirstStage ); + + // Effectively construct Volume + // Avoid to have two time same code for Angle and Point definition + void VolumeMaker(G4int TelescopeNumber , + G4ThreeVector MMpos , + G4RotationMatrix* MMrot , + bool wFirstStage , + G4LogicalVolume* world); + + + //////////////////////////////////////////////////// + //// Inherite from HeliosModule 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); + + // Give the static TInteractionCoordinates from VDetector to the classes + // deriving from HeliosModule + // This is mandatory since the Helios*** does not derive from VDetector + void SetInterCoordPointer(TInteractionCoordinates* interCoord); + TInteractionCoordinates* GetInterCoordPointer() {return ms_InterCoord;}; + + + //////////////////////////////////////////////////// + ///////////////Private intern Data////////////////// + //////////////////////////////////////////////////// +private: + // Interaction Coordinates coming from VDetector through the + // SetInteractionCoordinatesPointer method + TInteractionCoordinates* ms_InterCoord; + + // 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_Y128 ; // Bottom Left Corner Position Vector + vector<G4ThreeVector> m_X128_Y1 ; // Bottom Right Corner Position Vector + vector<G4ThreeVector> m_X128_Y128 ; // 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 ; // | +}; + + + +namespace HELIOSDUMMYSHAPE +{ + // Resolution + + const G4double ResoFirstStage = 0.0085; // = 20 keV of Resolution // Unit is MeV/2.35 + const G4double ResoTimeGpd = 0.4255; // = 1ns // 0.212765957;// = 500ps // Unit is ns/2.35 + const G4double ResoPosZ = 0.4255;// = 1mm for Helios // Unit is mm/2.35 + + // Geometry for the mother volume containing the different layers of your dummy shape module + const G4double FaceFrontWidth = 1.2*cm; + const G4double FaceFrontLength = 5.6*cm; + const G4double FaceBackWidth = 1.2*cm; + const G4double FaceBackLength = 5.6*cm; + const G4double Thickness = 0.1*cm; // ie: thickness 1 mm + const G4double InterStageDistance = 5*mm; + // for testing the excitation energy reconstruction +// const G4double Length = 4*cm; +// const G4double InterStageDistance = 15*mm; + + // First stage + const G4double FirstStageFaceWidth = 0.9*cm; + const G4double FirstStageFaceLength = 5.05*cm; + const G4double FirstStageThickness = 700*micrometer; +// const G4double FirstStageThickness = 2*mm; +// for testing the excitation energy reconstruction +// const G4double FirstStageThickness = 1.3*cm; + const G4int NumberOfStrips = 1; // PSD resistive strip + //const G4int NumberOfStrips = 500; // 100 um strip pitch + + + // Starting at the front of the first stage and going to the third stage + //const G4double FirstStage_PosZ = Thickness* -0.5 + 0.5*FirstStageThickness; + const G4double FirstStage_PosZ = 0.; +} + +#endif diff --git a/NPSimulation/include/HeliosModule.hh b/NPSimulation/include/HeliosModule.hh new file mode 100644 index 000000000..a8771ddb8 --- /dev/null +++ b/NPSimulation/include/HeliosModule.hh @@ -0,0 +1,95 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class is an Abstract Base Class (ABC) from which should * + * derive all different modules from the Helios detector. * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef HeliosModule_h +#define HeliosModule_h 1 + +// C++ headers +#include <string> +#include <vector> + +// G4 headers +#include "G4LogicalVolume.hh" +#include "G4Event.hh" +#include "G4MultiFunctionalDetector.hh" + +// NPTool - ROOT headers +#include "TInteractionCoordinates.h" +#include "THeliosData.h" + +using namespace std; + + + +class HeliosModule +{ +public: + HeliosModule(); + virtual ~HeliosModule(); + +public: + // Read stream at Configfile to pick-up parameters of detector (Position,...) + virtual void ReadConfiguration(string Path) = 0; + + // Construct detector and inialise sensitive part. + virtual void ConstructDetector(G4LogicalVolume* world) = 0; + + // Read sensitive part of a the HeliosModule detector and fill the Root tree. + virtual void ReadSensitive(const G4Event* event) = 0; + + // Add Detector branch to the ROOT output tree + virtual void InitializeRootOutput(); + + // Initialize all scorers necessary for each detector + virtual void InitializeScorers() = 0; + + // Give the static TInteractionCoordinates from VDetector to the classes + // deriving from HeliosModule + // This is mandatory since the Helios*** does not derive from VDetector + virtual void SetInterCoordPointer(TInteractionCoordinates* interCoord) = 0; + virtual TInteractionCoordinates* GetInterCoordPointer() = 0; + + // Initialize the Index map for the different modules of Gaspard + void InitializeIndex(); + +public: + THeliosData* GetEventPointer() {return ms_Event;}; + +protected: + // Class to store the result of simulation + static THeliosData* ms_Event; + + // Set to true if you want this stage on you telescope + vector<bool> m_wFirstStage; + + // Set to true if you want to see Telescope Frame in your visualisation + bool m_non_sensitive_part_visiualisation; + +protected: + // First stage Associate Scorer + G4MultiFunctionalDetector* m_FirstStageScorer; + + +protected: + map<string, int> m_index; +}; + +#endif diff --git a/NPSimulation/include/HeliosScorers.hh b/NPSimulation/include/HeliosScorers.hh new file mode 100644 index 000000000..0f635a54e --- /dev/null +++ b/NPSimulation/include/HeliosScorers.hh @@ -0,0 +1,249 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class holds all the scorers needed by the * + * Helios*** objects. * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef HELIOSScorer_h +#define HELIOSScorer_h 1 + +#include "G4VPrimitiveScorer.hh" +#include "G4THitsMap.hh" + +namespace HELIOSSCORERS +{ + +class HeliosScorerFirstStageEnergy : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageEnergy(G4String name, G4String volumeName, G4int depth = 0); + virtual ~HeliosScorerFirstStageEnergy(); + +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: + G4String m_VolumeName; + G4int HCID; + G4THitsMap<G4double>* EvtMap; +}; + + +class HeliosScorerFirstStageFrontStripDummyShape : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageFrontStripDummyShape(G4String name, G4int depth = 0, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageFrontStripDummyShape(); + +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 HeliosScorerFirstStageBackStripDummyShape : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageBackStripDummyShape(G4String name, G4int depth = 0, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageBackStripDummyShape(); + +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 HeliosScorerFirstStageFrontStripSquare : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageFrontStripSquare(G4String name, G4int depth = 0, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageFrontStripSquare(); + +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 HeliosScorerFirstStageBackStripSquare : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageBackStripSquare(G4String name, G4int depth = 0, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageBackStripSquare(); + +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 HeliosScorerFirstStageFrontStripTrapezoid : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageFrontStripTrapezoid(G4String name, G4int depth = 0, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageFrontStripTrapezoid(); + +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 HeliosScorerFirstStageBackStripTrapezoid : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageBackStripTrapezoid(G4String name, G4int depth = 0, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageBackStripTrapezoid(); + +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 HeliosScorerFirstStageFrontStripAnnular : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageFrontStripAnnular(G4String name, G4int depth = 0, G4double StripPlaneSize = 98, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageFrontStripAnnular(); + +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: + G4double m_StripPlaneSize; + G4int m_NumberOfStrip ; + G4int HCID; + G4THitsMap<G4double>* EvtMap; +}; + + + +class HeliosScorerFirstStageBackStripAnnular : public G4VPrimitiveScorer +{ +public: // with description + HeliosScorerFirstStageBackStripAnnular(G4String name, G4int depth = 0, G4double StripPlaneSize = 98, G4int NumberOfStrip = 128); + virtual ~HeliosScorerFirstStageBackStripAnnular(); + +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: + G4double m_StripPlaneSize; + G4int m_NumberOfStrip ; + G4int HCID; + G4THitsMap<G4double>* EvtMap; +}; + +} + +using namespace HELIOSSCORERS; +#endif diff --git a/NPSimulation/include/MyMagneticField.hh b/NPSimulation/include/MyMagneticField.hh new file mode 100644 index 000000000..252d61951 --- /dev/null +++ b/NPSimulation/include/MyMagneticField.hh @@ -0,0 +1,50 @@ + +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class defines the Magnetic field for Helios * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#ifndef MyMagneticField_H +#define MyMagneticField_H 1 + +#include "globals.hh" +#include "G4MagneticField.hh" +#include "G4ThreeVector.hh" + +class MyMagneticField : public G4MagneticField +{ + public: + MyMagneticField(const G4ThreeVector& FieldVector); + ~MyMagneticField(); + + void GetFieldValue( const G4double Point[3], + G4double *Bfield ) const; + + void SetFieldValue( const G4ThreeVector& newFieldValue ); + + + private: + G4double Bz; + G4double fFieldComponents[3]; + G4double rmax; + G4double zmax; + G4double zmin; +}; + +#endif diff --git a/NPSimulation/src/DetectorConstruction.cc b/NPSimulation/src/DetectorConstruction.cc index 3993aed59..d963b9995 100644 --- a/NPSimulation/src/DetectorConstruction.cc +++ b/NPSimulation/src/DetectorConstruction.cc @@ -39,6 +39,12 @@ #include "G4ios.hh" #include "G4String.hh" #include "G4RotationMatrix.hh" +#include "MyMagneticField.hh" +#include "G4FieldManager.hh" +#include "G4TransportationManager.hh" +#include "G4ChordFinder.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4SubtractionSolid.hh" // Detector class #include "../../NPLib/DetectorList.inc" @@ -48,6 +54,7 @@ #include "DummyDetector.hh" #include "Eurogam.hh" #include "GaspardTracker.hh" +#include "Helios.hh" #include "HydeTracker.hh" #include "MUST2Array.hh" #include "Paris.hh" @@ -166,6 +173,9 @@ void DetectorConstruction::ReadConfigurationFile(string Path){ string LineBuffer; string DataBuffer; + // needed for Magnetic field + double Bz=0.; + bool cAddThinSi = false; bool cComptonTelescope = false; bool cDummy = false; @@ -181,6 +191,8 @@ void DetectorConstruction::ReadConfigurationFile(string Path){ bool cSharc = false; bool cShield = false; bool cW1 = false; + bool cHelios = false; + bool check_MField = false; int VerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel(); @@ -475,7 +487,121 @@ void DetectorConstruction::ReadConfigurationFile(string Path){ AddDetector(myDetector) ; #endif } - + //////////////////////////////////////////// + //////////// Search for Helios //////////// + //////////////////////////////////////////// + else if (LineBuffer.compare(0, 6, "Helios") == 0 && cHelios == false) { +#ifdef INC_HELIOS + cHelios = true ; + G4cout << "//////// Helios detector ////////" << G4endl ; + + ConfigFile >> DataBuffer ; + if (DataBuffer.compare(0, 7, "MField=") == 0){ + check_MField = true; + ConfigFile >> DataBuffer ; + Bz = atof(DataBuffer.c_str()) ; + G4cout << "//////// Magentic Field set at Bz= " << Bz << " ////////" << G4endl ; + } + + // Instantiate the new array as a VDetector Object + VDetector* myDetector = new Helios() ; + + // Read Position of Telescope + ConfigFile.close() ; + myDetector->ReadConfiguration(Path) ; + ConfigFile.open(Path.c_str()) ; + + // Add array to the VDetector Vector + AddDetector(myDetector) ; + + //------------------------------world volume + // + + // Aluminium material + G4double a= 26.98 * g / mole; + G4double density = 2.7 * g / cm3; + G4double z = 13.; + G4Material* Aluminium = new G4Material("Aluminium", z, a, density); + + + + // Add the Aluminium rod + G4double Al_rod_x = 1. * cm; + G4double Al_rod_y = 1. * cm; + //G4double Al_rod_z = 20.0 * cm; + G4double Al_rod_z = 40.0 * cm; + G4Box* Al_rod_box + = new G4Box("Al_rod_box", Al_rod_x, Al_rod_y, Al_rod_z); + + G4Tubs* Al_rod_tub + = new G4Tubs("Al_rod_tub", 0, 0.5*cm, Al_rod_z+1.*mm, 0.*deg, 360*deg); + + G4SubtractionSolid* Al_rod=new G4SubtractionSolid("Rod",Al_rod_box, Al_rod_tub, 0, G4ThreeVector(0.,0.,0.)); + + G4LogicalVolume* Al_rod_log = new G4LogicalVolume(Al_rod, Aluminium, "Al_rod", 0, 0, 0); + + + new G4PVPlacement(0, G4ThreeVector(0.,0., Al_rod_z + 12.5*cm), Al_rod_log, "Al_rod", world_log, false, 0); + new G4PVPlacement(0, G4ThreeVector(0.,0., -(Al_rod_z + 12.5*cm)), Al_rod_log, "Al_rod", world_log, false, 1); + + + + // Add the Aluminium chamber + G4double Al_chamber_rmin = 50. * cm; + G4double Al_chamber_rmax = 55. * cm; + G4double Al_chamber_z = 100.0 * cm; + + //G4Tubs* Al_chamber_tub + // = new G4Tubs("Al_chamber_tub", Al_chamber_rmin, Al_chamber_rmax, Al_chamber_z, 0.*deg, 180*deg); + G4Tubs* Al_chamber_tub + = new G4Tubs("Al_chamber_tub", Al_chamber_rmin, Al_chamber_rmax, Al_chamber_z, 0.*deg, 360*deg); + + G4LogicalVolume* Al_chamber_log = new G4LogicalVolume(Al_chamber_tub, Aluminium, "Al_chamber", 0, 0, 0); + + G4RotationMatrix* RotZ = new G4RotationMatrix(); + RotZ->rotateZ(-90*deg); + + new G4PVPlacement(RotZ, G4ThreeVector(0.,0.,0.), Al_chamber_log, "Al_chamber", world_log, false, 0); + + + G4VisAttributes* VisAtt1 = new G4VisAttributes(G4Colour(0.2, 0.5, 0.8)); + Al_rod_log->SetVisAttributes(VisAtt1); + G4VisAttributes* VisAtt2 = new G4VisAttributes(G4Colour(0., 0.5, 0.3)); + Al_chamber_log->SetVisAttributes(VisAtt2); + + + //------------------------------------------------------------------------- + // add also My Magnetic field + //------------------------------------------------------------------------- + + + static G4bool fieldIsInitialized = false; + + if(!fieldIsInitialized) + { + MyMagneticField* myField = new MyMagneticField(G4ThreeVector(0.,0.,Bz)); + + G4FieldManager* fieldMgr + = G4TransportationManager::GetTransportationManager() + ->GetFieldManager(); + fieldMgr->SetDetectorField(myField); + + + /* + G4MagIntegratorStepper *pItsStepper; + G4ChordFinder* pChordFinder= new G4ChordFinder(myField, + 1.0e-2*mm, // stepper size + pItsStepper=0); + fieldMgr->SetChordFinder(pChordFinder); + */ + + fieldMgr->CreateChordFinder(myField); + + fieldIsInitialized = true; + } + +#endif + } //////////////////////////////////////////// //////////// Search for Target ///////////// //////////////////////////////////////////// diff --git a/NPSimulation/src/Helios.cc b/NPSimulation/src/Helios.cc new file mode 100644 index 000000000..63ed2f6c7 --- /dev/null +++ b/NPSimulation/src/Helios.cc @@ -0,0 +1,194 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class manages different shapes of module for the Helios * + * It allows to have Helios geometries with an * + * heterogeneous set of modules * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <fstream> +#include <sstream> +#include <string> +#include <cmath> + +// NPTool headers +#include "Helios.hh" +//#include "HeliosDetSquare.hh" +//#include "HeliosDetTrapezoid.hh" +//#include "HeliosDetAnnular.hh" +#include "HeliosDetDummyShape.hh" + +using namespace std; + + +Helios::Helios() +{ +} + + + +Helios::~Helios() +{ +} + + + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void Helios::ReadConfiguration(string Path) +{ + // open configuration file + ifstream ConfigFile; + ConfigFile.open(Path.c_str()); + + // bool HeliosSquare = false; + // bool HeliosTrapezoid = false; + // bool HeliosAnnular = false; + bool HeliosDummyShape = false; + + string LineBuffer; + while (!ConfigFile.eof()) { + getline(ConfigFile, LineBuffer); + /* + if (LineBuffer.compare(0, 12, "HeliosSquare") == 0 && HeliosSquare == false) { + HeliosSquare = true; + + // instantiate a new "detector" corresponding to the Square elements + HeliosModule* myDetector = new HeliosSquare(); + + // read part of the configuration file corresponding to square elements + ConfigFile.close(); + myDetector->ReadConfiguration(Path); + ConfigFile.open(Path.c_str()); + + // ms_InterCoord comes from VDetector + myDetector->SetInterCoordPointer(ms_InterCoord); + + // store HeliosSquare "detector" + m_Modules.push_back(myDetector); + } + else if (LineBuffer.compare(0, 15, "HeliosTrapezoid") == 0 && HeliosTrapezoid == false) { + HeliosTrapezoid = true; + + // instantiate a new "detector" corresponding to the Trapezoid elements + HeliosModule* myDetector = new HeliosTrapezoid(); + + // read part of the configuration file corresponding to trapezoid elements + ConfigFile.close(); + myDetector->ReadConfiguration(Path); + ConfigFile.open(Path.c_str()); + + // ms_InterCoord comes from VDetector + myDetector->SetInterCoordPointer(ms_InterCoord); + + // store HeliosTrapezoid "detector" + m_Modules.push_back(myDetector); + } + else if (LineBuffer.compare(0, 13, "HeliosAnnular") == 0 && HeliosAnnular == false) { + HeliosAnnular = true; + + // instantiate a new "detector" corresponding to the Trapezoid elements + HeliosModule* myDetector = new HeliosAnnular(); + + // read part of the configuration file corresponding to trapezoid elements + ConfigFile.close(); + myDetector->ReadConfiguration(Path); + ConfigFile.open(Path.c_str()); + + // ms_InterCoord comes from VDetector + myDetector->SetInterCoordPointer(ms_InterCoord); + + // store HeliosTrapezoid "detector" + m_Modules.push_back(myDetector); + } + else + */ + if (LineBuffer.compare(0, 16, "HeliosDummyShape") == 0 && HeliosDummyShape == false) { + HeliosDummyShape = true; + + // instantiate a new "detector" corresponding to the Shape elements + // The HeliosSquare class should be replaced by the + // HeliosShape class you need to define + HeliosModule* myDetector = new HeliosDetDummyShape(); + + // read part of the configuration file corresponding to shape elements + ConfigFile.close(); + myDetector->ReadConfiguration(Path); + ConfigFile.open(Path.c_str()); + + // ms_InterCoord comes from VDetector + myDetector->SetInterCoordPointer(ms_InterCoord); + + // store HeliosShape "detector" + m_Modules.push_back(myDetector); + } + } +} + + + +// Construct detector and initialize sensitive part. +// Called After DetecorConstruction::AddDetector Method +void Helios::ConstructDetector(G4LogicalVolume* world) +{ + // loop on sub-detectors belonging to Helios + int nbDetectors = m_Modules.size(); + for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ConstructDetector(world); +} + + + +// Connect the GaspardTrackingData class to the output TTree +// of the simulation +void Helios::InitializeRootOutput() +{ + // loop on sub-detectors belonging to Helios + int nbDetectors = m_Modules.size(); + for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeRootOutput(); +} + + + +// Initialize all scorers necessary for each detector +void Helios::InitializeScorers() +{ + // loop on sub-detectors belonging to Helios + int nbDetectors = m_Modules.size(); + for (int i = 0; i < nbDetectors; i++) m_Modules[i]->InitializeScorers(); +} + + + +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAction +void Helios::ReadSensitive(const G4Event* event) +{ + // Before looping on each sub-detector, clear the static variable + // ms_InterCoord + // This is done on the first element of the m_Modules vector. + // This should be done here since this variable (of type TIneractionCoordinates) + // deals with multiplicity of events > 1. + m_Modules[0]->GetInterCoordPointer()->Clear(); + + // We do the same for the static variable ms_Event + m_Modules[0]->GetEventPointer()->Clear(); + + // loop on sub-detectors belonging to Helios + int nbDetectors = m_Modules.size(); + for (int i = 0; i < nbDetectors; i++) m_Modules[i]->ReadSensitive(event); +} diff --git a/NPSimulation/src/HeliosDetDummyShape.cc b/NPSimulation/src/HeliosDetDummyShape.cc new file mode 100644 index 000000000..3d907a977 --- /dev/null +++ b/NPSimulation/src/HeliosDetDummyShape.cc @@ -0,0 +1,886 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: Define a dummy module for the Helios detector * + * The goal of this class is to be a starting point to create a * + * new shape to be added to the Helios detector. * + * * + *---------------------------------------------------------------------------* + * Comment: Here tyhe shape is a rectangular box (w~1cm, l~5cm, d~700um) * + * * + * * + *****************************************************************************/ + +// C++ headers +#include <sstream> +#include <string> +#include <cmath> + +// G4 Geometry headers +#include "G4Trd.hh" +#include "G4Box.hh" +#include "G4Trap.hh" +#include "G4SubtractionSolid.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 "HeliosDetDummyShape.hh" +#include "GeneralScorers.hh" +#include "HeliosScorers.hh" +#include "RootOutput.h" + +// CLHEP +#include "CLHEP/Random/RandGauss.h" + +using namespace std; +using namespace CLHEP; +using namespace HELIOSDUMMYSHAPE; + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +HeliosDetDummyShape::HeliosDetDummyShape() +{ + ms_InterCoord = 0; +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +HeliosDetDummyShape::~HeliosDetDummyShape() +{ +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void HeliosDetDummyShape::AddModule(G4ThreeVector X1_Y1 , + G4ThreeVector X128_Y1 , + G4ThreeVector X1_Y128 , + G4ThreeVector X128_Y128 , + bool wFirstStage ) +{ + m_DefinitionType.push_back(true) ; + + m_X1_Y1.push_back(X1_Y1) ; + m_X128_Y1.push_back(X128_Y1) ; + m_X1_Y128.push_back(X1_Y128) ; + m_X128_Y128.push_back(X128_Y128) ; + m_wFirstStage.push_back(wFirstStage) ; + + 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) ; + + +// m_wNumberStrip.push_back(wNumberStrip); +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void HeliosDetDummyShape::AddModule(G4double R , + G4double Theta , + G4double Phi , + G4double beta_u , + G4double beta_v , + G4double beta_w , + bool wFirstStage ) +{ + G4ThreeVector empty = G4ThreeVector(0, 0, 0); + + m_DefinitionType.push_back(false); + + m_X1_Y1.push_back(empty) ; + m_X128_Y1.push_back(empty) ; + m_X1_Y128.push_back(empty) ; + m_X128_Y128.push_back(empty) ; + + 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_wFirstStage.push_back(wFirstStage) ; +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +void HeliosDetDummyShape::VolumeMaker(G4int TelescopeNumber, + G4ThreeVector MMpos, + G4RotationMatrix* MMrot, + bool wFirstStage, + G4LogicalVolume* world) +{ + G4double NbrTelescopes = TelescopeNumber ; + G4String DetectorNumber ; + ostringstream Number ; + Number << NbrTelescopes ; + DetectorNumber = Number.str() ; + + ///////////////////////////////////////////////////////////////// + /////////////////Element Definition /////////////////////////// + //////////////////////////////////////////////////////////////// + G4String symbol; + G4double density = 0. , a = 0, z = 0; + G4int ncomponents = 0; + + //////////////////////////////////////////////////////////////// + /////////////////Material Definition /////////////////////////// + //////////////////////////////////////////////////////////////// + + // Si + a = 28.0855 * g / mole; + density = 2.321 * g / cm3; + G4Material* Silicon = new G4Material("Si", z = 14., a, density); + + // Copper + a = 63.546 * g / mole; + density = 8.96 * g / cm3; + G4Material* Copper; + Copper = new G4Material("Cu", z = 29., a, density); + + // Vacuum + 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); + + density = 0.000000001 * mg / cm3; + G4Material* Vacuum = new G4Material("Vacuum", density, ncomponents = 2); + Vacuum->AddElement(N, .7); + Vacuum->AddElement(O, .3); + + // PCB + G4Element* Si = new G4Element("Silicon" , symbol = "Si" , z = 14 , a = 28.0855 * g / mole); + G4Element* C = new G4Element("Carbon" , symbol = "C" , z = 6 , a = 12.011 * g / mole); + G4Element* H = new G4Element("Hydrogen" , symbol = "H" , z = 1 , a = 1.0079 * g / mole); + G4Element* Br = new G4Element("Bromine" , symbol = "Br" , z = 35 , a = 79.904 * g / mole); + + density = 1.7 * g / cm3; + G4Material* PCB = new G4Material("PCB", density, ncomponents = 5); + PCB->AddElement(Si, .181); + PCB->AddElement(O, .406); + PCB->AddElement(C, .278); + PCB->AddElement(H, .068); + PCB->AddElement(Br, .067); + + + //////////////////////////////////////////////////////////////// + ////////////// 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 = "HeliosDummyShape" + DetectorNumber ; + + G4Box* solidHeliosDummyShape = new G4Box(Name, 0.5*FaceFrontWidth, 0.5*FaceFrontLength, 0.5*Thickness); + G4LogicalVolume* logicHeliosDummyShape = new G4LogicalVolume(solidHeliosDummyShape, Vacuum, Name, 0, 0, 0); + + PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot, MMpos) , + logicHeliosDummyShape , + Name , + world , + false , + 0); + + //logicHeliosDummyShape->SetVisAttributes(G4VisAttributes::Invisible); + logicHeliosDummyShape->SetVisAttributes(G4VisAttributes(G4Colour(1, 0., 0.0))); // red + //if (m_non_sensitive_part_visiualisation) logicHeliosDummyShape->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90))); + + //Place two marker to identify the u and v axis on silicon face: + //marker are placed a bit before the silicon itself so they don't perturbate simulation + //Uncomment to help debugging or if you want to understand the way the code work. + //I should recommand to Comment it during simulation to avoid perturbation of simulation + //Remember G4 is limitationg step on geometry constraints. + /* + G4ThreeVector positionMarkerU = CT*0.98 + MMu*SiliconFace/4; + G4Box* solidMarkerU = new G4Box( "solidMarkerU" , SiliconFace/4 , 1*mm , 1*mm ) ; + G4LogicalVolume* logicMarkerU = new G4LogicalVolume( solidMarkerU , Vacuum , "logicMarkerU",0,0,0) ; + PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot,positionMarkerU),logicMarkerU,"MarkerU",world,false,0) ; + + G4VisAttributes* MarkerUVisAtt= new G4VisAttributes(G4Colour(0.,0.,0.5));//blue + logicMarkerU->SetVisAttributes(MarkerUVisAtt); + + G4ThreeVector positionMarkerV = CT*0.98 + MMv*SiliconFace/4; + G4Box* solidMarkerV = new G4Box( "solidMarkerU" , 1*mm , SiliconFace/4 , 1*mm ) ; + G4LogicalVolume* logicMarkerV = new G4LogicalVolume( solidMarkerV , Vacuum , "logicMarkerV",0,0,0) ; + PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot,positionMarkerV),logicMarkerV,"MarkerV",world,false,0) ; + + G4VisAttributes* MarkerVVisAtt= new G4VisAttributes(G4Colour(0.,0.5,0.5));//green + logicMarkerV->SetVisAttributes(MarkerVVisAtt); + */ + + //////////////////////////////////////////////////////////////// + ///////////////// First Stage Construction ///////////////////// + //////////////////////////////////////////////////////////////// + if (wFirstStage) { + // Silicon detector itself + G4ThreeVector positionFirstStage = G4ThreeVector(0, 0, FirstStage_PosZ); + + G4Box* solidFirstStage = new G4Box("solidFirstStage", 0.5*FirstStageFaceWidth, 0.5*FirstStageFaceLength, 0.5*FirstStageThickness); + G4LogicalVolume* logicFirstStage = new G4LogicalVolume(solidFirstStage, Silicon, "logicFirstStage", 0, 0, 0); + + PVPBuffer = new G4PVPlacement(0, + positionFirstStage, + logicFirstStage, +// "G" + DetectorNumber + "FirstStage", + Name + "_FirstStage", + logicHeliosDummyShape, + false, + 0); + + // Set First Stage sensible + logicFirstStage->SetSensitiveDetector(m_FirstStageScorer); + + ///Visualisation of FirstStage Strip + G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(1., 0., 0.)); // red + logicFirstStage->SetVisAttributes(FirstStageVisAtt); + + + /* + G4Box* solidFirstStageframeIn = new G4Box("solidFirstStageframeIn", 0.5*FaceFront, 0.5*FaceFront, 0.5*1*mm); + G4Box* solidFirstStageframeOut = new G4Box("solidFirstStageframeOut", 0.5*FirstStageFace, 0.5*FirstStageFace, 0.5*1.1*mm); + G4SubtractionSolid* solidFirstStagePCB = new G4SubtractionSolid("solidFirstStagePCB",solidFirstStageframeIn, solidFirstStageframeOut, 0, G4ThreeVector(0.,0.,0.)); + G4LogicalVolume* logicFirstStagePCB = new G4LogicalVolume(solidFirstStagePCB, PCB, "logicFirstStagePCB", 0, 0, 0); + + G4ThreeVector positionFirstStagePCB(0.,0.,-0.5*Length+0.5*1*mm); + PVPBuffer = new G4PVPlacement(0, + positionFirstStagePCB, + logicFirstStagePCB, +// "G" + DetectorNumber + "FirstStage", + Name + "_FirstStagePCB", + logicHeliosDummyShape, + false, + 0); + + + ///Visualisation of FirstStage Strip + G4VisAttributes* FirstStagePCBVisAtt = new G4VisAttributes(G4Colour(0., 0.8, 0.5)); // green + logicFirstStagePCB->SetVisAttributes(FirstStagePCBVisAtt); + */ + + } + + + +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// Virtual Method of VDetector class + +// Read stream at Configfile to pick-up parameters of detector (Position,...) +// Called in DetecorConstruction::ReadDetextorConfiguration Method +void HeliosDetDummyShape::ReadConfiguration(string Path) +{ + ifstream ConfigFile; + ConfigFile.open(Path.c_str()); + string LineBuffer, DataBuffer; + + // A:X1_Y1 --> X:1 Y:1 + // B:X128_Y1 --> X:128 Y:1 + // C:X1_Y128 --> X:1 Y:128 + // D:X128_Y128 --> X:128 Y:128 + + 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 ; + int FIRSTSTAGE = 0 , SECONDSTAGE = 0 , THIRDSTAGE = 0 ; + int NSTRIP = 128; + + bool ReadingStatus = false; + + bool check_A = false; + bool check_C = false; + bool check_B = false; + bool check_D = false; + + bool check_Theta = false; + bool check_Phi = false; + bool check_R = false; + bool check_beta = false; + + bool check_FirstStage = false; + bool check_SecondStage = false; + bool check_ThirdStage = false; + bool check_NStrip = false; + bool checkVis = false; + + while (!ConfigFile.eof()) { + getline(ConfigFile, LineBuffer); + if (LineBuffer.compare(0, 16, "HeliosDummyShape") == 0) { + G4cout << "///" << G4endl ; + G4cout << "DummyShape element found: " << G4endl ; + ReadingStatus = true ; + } + + while (ReadingStatus) { + ConfigFile >> DataBuffer; + // Comment Line + if (DataBuffer.compare(0, 1, "%") == 0) {/*do nothing */;} + + // 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, 8, "X128_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 << "X128 Y1 corner position : " << B << endl; + } + else if (DataBuffer.compare(0, 8, "X1_Y128=") == 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 Y128 corner position : " << C << endl; + } + else if (DataBuffer.compare(0, 10, "X128_Y128=") == 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 << "X128 Y128 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, 11, "FIRSTSTAGE=") == 0) { + check_FirstStage = true ; + ConfigFile >> DataBuffer; + FIRSTSTAGE = atoi(DataBuffer.c_str()) ; + } + else if (DataBuffer.compare(0, 12, "SECONDSTAGE=") == 0) { + check_SecondStage = true ; + ConfigFile >> DataBuffer; + SECONDSTAGE = atoi(DataBuffer.c_str()) ; + } + else if (DataBuffer.compare(0, 11, "THIRDSTAGE=") == 0) { + check_ThirdStage = true ; + ConfigFile >> DataBuffer; + THIRDSTAGE = atoi(DataBuffer.c_str()) ; + } + + else if (DataBuffer.compare(0, 7, "NSTRIP=") == 0) { + check_NStrip = true ; + ConfigFile >> DataBuffer; + NSTRIP = atoi(DataBuffer.c_str()) ; + } + + 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 G4cout << "WARNING: Wrong Token, HeliosDummyShape: DummyShape Element not added" << G4endl; + + // Add The previously define telescope + // With position method + if ((check_A && check_B && check_C && check_D && check_FirstStage && checkVis) && + !(check_Theta && check_Phi && check_R)) { + ReadingStatus = false; + check_A = false; + check_C = false; + check_B = false; + check_D = false; + check_FirstStage = false; + checkVis = false; + + AddModule(A, B, C, D, FIRSTSTAGE == 1); + } + + // With angle method + if ((check_Theta && check_Phi && check_R && check_FirstStage && checkVis) && + !(check_A && check_B && check_C && check_D)) { + ReadingStatus = false; + check_Theta = false; + check_Phi = false; + check_R = false; + check_beta = false; + check_FirstStage = false; + checkVis = false; + + AddModule(R, Theta, Phi, beta_u, beta_v, beta_w, FIRSTSTAGE == 1); + } + } + } +} + + + +// Construct detector and inialise sensitive part. +// Called After DetecorConstruction::AddDetector Method +void HeliosDetDummyShape::ConstructDetector(G4LogicalVolume* world) +{ + G4RotationMatrix* MMrot = NULL ; + G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ; + bool FirstStage = true ; + + G4int NumberOfTelescope = m_DefinitionType.size() ; + + for (G4int i = 0; i < NumberOfTelescope; 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 + MMu = m_X128_Y1[i] - m_X1_Y1[i]; + MMu = MMu.unit(); + + MMv = m_X1_Y128[i] - m_X1_Y1[i]; + MMv = MMv.unit(); + + // G4double MMscal = MMu.dot(MMv); + + MMw = MMu.cross(MMv); +// if (MMw.z() > 0) MMw = MMv.cross(MMu) ; + MMw = MMw.unit(); + + MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4; + + // Passage Matrix from Lab Referential to Telescope Referential + // MUST2 + MMrot = new G4RotationMatrix(MMu, MMv, MMw); + // translation to place Telescope + MMpos = MMw * Thickness * 0.5 + MMCenter; + } + + // 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); + MMw = G4ThreeVector(wX, wY, wZ); + + // vector corresponding to the center of the module + G4ThreeVector CT = MMw; + + // 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); + + MMw = MMw.unit(); + MMu = MMw.cross(Y); + MMv = MMw.cross(MMu); + MMv = MMv.unit(); + MMu = MMu.unit(); + + // Passage Matrix from Lab Referential to Telescope Referential + // MUST2 + MMrot = new G4RotationMatrix(MMu, MMv, MMw); + // Telescope is rotate of Beta angle around MMv axis. + MMrot->rotate(m_beta_u[i], MMu); + MMrot->rotate(m_beta_v[i], MMv); + MMrot->rotate(m_beta_w[i], MMw); + // translation to place Telescope + MMpos = MMw * Thickness * 0.5 + CT ; + } + + FirstStage = m_wFirstStage[i] ; + + VolumeMaker(i + 1, MMpos, MMrot, FirstStage, world); + } + + delete MMrot ; +} + + + +// Connect the GaspardTrackingData class to the output TTree +// of the simulation +void HeliosDetDummyShape::InitializeRootOutput() +{ +} + + + +// Set the TinteractionCoordinates object from VDetector to the present class +void HeliosDetDummyShape::SetInterCoordPointer(TInteractionCoordinates* interCoord) +{ + ms_InterCoord = interCoord; +} + + + +// Read sensitive part and fill the Root tree. +// Called at in the EventAction::EndOfEventAvtion +void HeliosDetDummyShape::ReadSensitive(const G4Event* event) +{ + ////////////////////////////////////////////////////////////////////////////////////// + //////////////////////// 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 the Scorer associate to the Silicon Strip + //Detector Number + G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/DetectorNumber") ; + DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ; + DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ; + + //Energy + G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/StripEnergy") ; + EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ; + Energy_itr = EnergyHitMap->GetMap()->begin() ; + + //Time of Flight + G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/StripTime") ; + TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ; + Time_itr = TimeHitMap->GetMap()->begin() ; + + //Strip Number X + G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/StripIDFront") ; + XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ; + X_itr = XHitMap->GetMap()->begin() ; + + //Strip Number Y + G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/StripIDBack"); + YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ; + Y_itr = YHitMap->GetMap()->begin() ; + + //Interaction Coordinate X + G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/InterCoordX") ; + PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ; + Pos_X_itr = PosXHitMap->GetMap()->begin() ; + + //Interaction Coordinate Y + G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/InterCoordY") ; + PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ; + Pos_Y_itr = PosYHitMap->GetMap()->begin() ; + + //Interaction Coordinate Z + G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/InterCoordZ") ; + PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ; + Pos_Z_itr = PosXHitMap->GetMap()->begin() ; + + //Interaction Coordinate Angle Theta + G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/InterCoordAngTheta") ; + AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ; + Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ; + + //Interaction Coordinate Angle Phi + G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("FirstStageScorerHeliosDummyShape/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(); + + + //G4cout << "sizeN:" << sizeN << G4endl; + + + if (sizeE != sizeT || sizeT != sizeX || sizeX != sizeY) { + //G4cout << "sizeE:" << sizeE << G4endl; + //G4cout << "sizeT:" << sizeT << G4endl; + //G4cout << "sizeX:" << sizeX << G4endl; + //G4cout << "sizeY:" << sizeY << G4endl; + //G4cout << "No match size Si Event Map: sE:" + // << sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << endl ; + return; + } + + // Loop on FirstStage number + for (G4int l = 0; l < sizeN; l++) { + G4int N = *(DetectorNumber_itr->second); + G4int NTrackID = DetectorNumber_itr->first - N; + + //G4cout <<"N:" <<N << G4endl; + //G4cout <<"DetectorNumber_itr->first:" << DetectorNumber_itr->first << G4endl; + //G4cout <<"NTrackID:" <<NTrackID << G4endl; + + + if (N > 0) { + // Fill detector number + ms_Event->SetHeliosFirstStageEDetectorNbr(m_index["DummyShape"] + N); + ms_Event->SetHeliosFirstStageTDetectorNbr(m_index["DummyShape"] + N); + + // Energy + for (G4int l = 0 ; l < sizeE ; l++) { + G4int ETrackID = Energy_itr->first - N; + G4double E = *(Energy_itr->second); + if (ETrackID == NTrackID) { + ms_Event->SetHeliosFirstStageEEnergy(RandGauss::shoot(E, ResoFirstStage)); + } + 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) { + T = RandGauss::shoot(T, ResoTimeGpd); + ms_Event->SetHeliosFirstStageTTime(RandGauss::shoot(T, ResoTimeGpd)); + } + Time_itr++; + } + + // Strip X + 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) { + ms_Event->SetHeliosFirstStageEStripNbr(int(X)); + ms_Event->SetHeliosFirstStageTStripNbr(int(X)); + } + X_itr++; + } + + + // Pos X + Pos_X_itr = PosXHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeX ; 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 < sizeX ; 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 < sizeX ; h++) { + G4int PosZTrackID = Pos_Z_itr->first - N ; + G4double PosZ = *(Pos_Z_itr->second) ; + if (PosZTrackID == NTrackID) { + ms_InterCoord->SetDetectedPositionZ(RandGauss::shoot(PosZ, ResoFirstStage)) ; // for helios !!!! + } + Pos_Z_itr++; + } + + // Angle Theta + Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); + for (G4int h = 0 ; h < sizeX ; 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 < sizeX ; h++) { + G4int AngPhiTrackID = Ang_Phi_itr->first - N ; + G4double AngPhi = *(Ang_Phi_itr->second) ; + if (AngPhiTrackID == NTrackID) { + ms_InterCoord->SetDetectedAnglePhi(AngPhi) ; + } + Ang_Phi_itr++; + } + + + DetectorNumber_itr++; + } + + // 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 HeliosDetDummyShape::InitializeScorers() +{ + // First stage Associate Scorer + m_FirstStageScorer = new G4MultiFunctionalDetector("FirstStageScorerHeliosDummyShape"); + G4VPrimitiveScorer* DetNbr = new GENERALSCORERS::PSDetectorNumber("DetectorNumber", "HeliosDummyShape", 0); + G4VPrimitiveScorer* TOF = new GENERALSCORERS::PSTOF("StripTime","HeliosDummyShape", 0); + G4VPrimitiveScorer* InteractionCoordinatesX = new GENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","HeliosDummyShape", 0); + G4VPrimitiveScorer* InteractionCoordinatesY = new GENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","HeliosDummyShape", 0); + G4VPrimitiveScorer* InteractionCoordinatesZ = new GENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","HeliosDummyShape", 0); + G4VPrimitiveScorer* InteractionCoordinatesAngleTheta = new GENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","HeliosDummyShape", 0); + G4VPrimitiveScorer* InteractionCoordinatesAnglePhi = new GENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","HeliosDummyShape", 0); + G4VPrimitiveScorer* Energy = new HeliosScorerFirstStageEnergy("StripEnergy", "HeliosDummyShape", 0); + G4VPrimitiveScorer* StripPositionX = new HeliosScorerFirstStageFrontStripDummyShape("StripIDFront", 0, NumberOfStrips); + G4VPrimitiveScorer* StripPositionY = new HeliosScorerFirstStageBackStripDummyShape("StripIDBack", 0, NumberOfStrips); + + //and register it to the multifunctionnal detector + m_FirstStageScorer->RegisterPrimitive(DetNbr); + m_FirstStageScorer->RegisterPrimitive(Energy); + m_FirstStageScorer->RegisterPrimitive(TOF); + m_FirstStageScorer->RegisterPrimitive(StripPositionX); + m_FirstStageScorer->RegisterPrimitive(StripPositionY); + m_FirstStageScorer->RegisterPrimitive(InteractionCoordinatesX); + m_FirstStageScorer->RegisterPrimitive(InteractionCoordinatesY); + m_FirstStageScorer->RegisterPrimitive(InteractionCoordinatesZ); + m_FirstStageScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta); + m_FirstStageScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi); + + + + + // Add All Scorer to the Global Scorer Manager + G4SDManager::GetSDMpointer()->AddNewDetector(m_FirstStageScorer); +} diff --git a/NPSimulation/src/HeliosModule.cc b/NPSimulation/src/HeliosModule.cc new file mode 100644 index 000000000..a8d2ef1b7 --- /dev/null +++ b/NPSimulation/src/HeliosModule.cc @@ -0,0 +1,63 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class is an Abstract Base Class (ABC) from which should * + * derive all different modules from the Helios detector. * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "HeliosModule.hh" +#include "RootOutput.h" + + +THeliosData *HeliosModule::ms_Event = 0; + + + +HeliosModule::HeliosModule() +{ + if (ms_Event == 0) ms_Event = new THeliosData(); + + InitializeRootOutput(); + InitializeIndex(); +} + + + +HeliosModule::~HeliosModule() +{ +} + + + +void HeliosModule::InitializeRootOutput() +{ + RootOutput *pAnalysis = RootOutput::getInstance(); + TTree *pTree = pAnalysis->GetTree(); + // if the branch does not exist yet, create it + if (!pTree->GetBranch("HELIOS")) + pTree->Branch("HELIOS", "THeliosData", &ms_Event); +} + + + +void HeliosModule::InitializeIndex() +{ + // m_index["Square"] = 0; + // m_index["Trapezoid"] = 100; + // m_index["Annular"] = 200; + m_index["DummyShape"] = 1000; +} diff --git a/NPSimulation/src/HeliosScorers.cc b/NPSimulation/src/HeliosScorers.cc new file mode 100644 index 000000000..5a64309ae --- /dev/null +++ b/NPSimulation/src/HeliosScorers.cc @@ -0,0 +1,642 @@ +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class holds all the scorers needed by the * + * Helios*** objects. * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +// G4 headers +#include "G4UnitsTable.hh" + +// NPTool headers +#include "GeneralScorers.hh" +#include "HeliosScorers.hh" + +#include "HeliosDetDummyShape.hh" +//#include "HeliosDetSquare.hh" +//#include "HeliosDetTrapezoid.hh" +//#include "HeliosDetAnnular.hh" +//using namespace HELIOSSQUARE; +//using namespace HELIOSTRAP; +//using namespace HELIOSANNULAR; +using namespace HELIOSDUMMYSHAPE; + +//....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...... +// FirstStage Energy Scorer (deal with multiple particle hit) +HeliosScorerFirstStageEnergy::HeliosScorerFirstStageEnergy(G4String name, G4String volumeName, G4int depth) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_VolumeName = volumeName; +} + +HeliosScorerFirstStageEnergy::~HeliosScorerFirstStageEnergy() +{ +} + +G4bool HeliosScorerFirstStageEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, m_VolumeName); + + // get energy + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->add(DetNbr + index, edep); + return TRUE; +} + +void HeliosScorerFirstStageEnergy::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageEnergy::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageEnergy::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageEnergy::DrawAll() +{ +} + +void HeliosScorerFirstStageEnergy::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl; + G4cout << " PrimitiveScorer " << GetName() << G4endl; + G4cout << " Number of entries " << EvtMap->entries() << G4endl; + std::map<G4int, G4double*>::iterator itr = EvtMap->GetMap()->begin(); + for (; itr != EvtMap->GetMap()->end(); itr++) { + G4cout << " copy no.: " << itr->first + << " energy deposit: " << G4BestUnit(*(itr->second), "Energy") + << G4endl; + } +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Front Strip position Scorer for DummyShape geometry +HeliosScorerFirstStageFrontStripDummyShape::HeliosScorerFirstStageFrontStripDummyShape(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip ; +} + +HeliosScorerFirstStageFrontStripDummyShape::~HeliosScorerFirstStageFrontStripDummyShape() +{ +} + +G4bool HeliosScorerFirstStageFrontStripDummyShape::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosDummyShape"); + + // get front strip number + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = HELIOSDUMMYSHAPE::FirstStageFaceWidth / m_NumberOfStrip; // longitudinal strip + + G4double temp = (POS(0) + HELIOSDUMMYSHAPE::FirstStageFaceWidth / 2.) / StripPitch ; + G4double X = int(temp) + 1 ; + + //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; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void HeliosScorerFirstStageFrontStripDummyShape::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageFrontStripDummyShape::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageFrontStripDummyShape::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageFrontStripDummyShape::DrawAll() +{ +} + +void HeliosScorerFirstStageFrontStripDummyShape::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Back Strip position Scorer for DummyShape geometry +HeliosScorerFirstStageBackStripDummyShape::HeliosScorerFirstStageBackStripDummyShape(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip ; +} + +HeliosScorerFirstStageBackStripDummyShape::~HeliosScorerFirstStageBackStripDummyShape() +{ +} + +G4bool HeliosScorerFirstStageBackStripDummyShape::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosDummyShape"); + + // get back strip number + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = HELIOSDUMMYSHAPE::FirstStageFaceLength / m_NumberOfStrip; // transversal strip + + G4double temp = (POS(1) + HELIOSDUMMYSHAPE::FirstStageFaceLength / 2.) / StripPitch ; + G4double X = int(temp) + 1 ; + //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; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void HeliosScorerFirstStageBackStripDummyShape::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageBackStripDummyShape::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageBackStripDummyShape::Clear() + +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageBackStripDummyShape::DrawAll() +{ +} + +void HeliosScorerFirstStageBackStripDummyShape::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + +/* +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Front Strip position Scorer for Square geometry +HeliosScorerFirstStageFrontStripSquare::HeliosScorerFirstStageFrontStripSquare(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip ; +} + +HeliosScorerFirstStageFrontStripSquare::~HeliosScorerFirstStageFrontStripSquare() +{ +} + +G4bool HeliosScorerFirstStageFrontStripSquare::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosSquare"); + + // get front strip + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = HeliosSQUARE::SiliconFace / m_NumberOfStrip; + + G4double temp = (POS(0) + HeliosSQUARE::SiliconFace / 2.) / StripPitch ; + G4double X = int(temp) + 1 ; + //Rare case where particle is close to edge of silicon plan + if (X == 129) X = 128; + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void HeliosScorerFirstStageFrontStripSquare::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageFrontStripSquare::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageFrontStripSquare::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageFrontStripSquare::DrawAll() +{ +} + +void HeliosScorerFirstStageFrontStripSquare::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Back Strip position Scorer for Square geometry +HeliosScorerFirstStageBackStripSquare::HeliosScorerFirstStageBackStripSquare(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip ; +} + +HeliosScorerFirstStageBackStripSquare::~HeliosScorerFirstStageBackStripSquare() +{ +} + +G4bool HeliosScorerFirstStageBackStripSquare::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosSquare"); + + // get back strip + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = HeliosSQUARE::SiliconFace / m_NumberOfStrip; + + G4double temp = (POS(1) + HeliosSQUARE::SiliconFace / 2.) / StripPitch ; + G4int temp2 = temp ; + G4double Y = temp2 + 1 ; + //Rare case where particle is close to edge of silicon plan + if (Y == 129) Y = 128; + + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, Y); + return TRUE; +} + +void HeliosScorerFirstStageBackStripSquare::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageBackStripSquare::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageBackStripSquare::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageBackStripSquare::DrawAll() +{ +} + +void HeliosScorerFirstStageBackStripSquare::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Front Strip position Scorer for Trapezoid geometry +HeliosScorerFirstStageFrontStripTrapezoid::HeliosScorerFirstStageFrontStripTrapezoid(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip; +} + +HeliosScorerFirstStageFrontStripTrapezoid::~HeliosScorerFirstStageFrontStripTrapezoid() +{ +} + +G4bool HeliosScorerFirstStageFrontStripTrapezoid::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosTrapezoid"); + + // get front strip + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = HeliosTRAP::FirstStageBaseLarge / m_NumberOfStrip; + + G4double temp = (POS(0) + HeliosTRAP::FirstStageBaseLarge / 2.) / StripPitch ; + G4double X = int(temp) + 1 ; + //Rare case where particle is close to edge of silicon plan + if (X == 129) X = 128; + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void HeliosScorerFirstStageFrontStripTrapezoid::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageFrontStripTrapezoid::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageFrontStripTrapezoid::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageFrontStripTrapezoid::DrawAll() +{ +} + +void HeliosScorerFirstStageFrontStripTrapezoid::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Back Strip position Scorer for Trapezoid geometry +HeliosScorerFirstStageBackStripTrapezoid::HeliosScorerFirstStageBackStripTrapezoid(G4String name, G4int depth, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_NumberOfStrip = NumberOfStrip; +} + +HeliosScorerFirstStageBackStripTrapezoid::~HeliosScorerFirstStageBackStripTrapezoid() +{ +} + +G4bool HeliosScorerFirstStageBackStripTrapezoid::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosTrapezoid"); + + // get back strip + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = HeliosTRAP::FirstStageHeight / m_NumberOfStrip; + + G4double temp = (POS(1) + HeliosTRAP::FirstStageHeight / 2.) / StripPitch ; + G4int temp2 = temp ; + G4double Y = temp2 + 1 ; + //Rare case where particle is close to edge of silicon plan + if (Y == 129) Y = 128; + + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, Y); + return TRUE; +} + +void HeliosScorerFirstStageBackStripTrapezoid::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageBackStripTrapezoid::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageBackStripTrapezoid::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageBackStripTrapezoid::DrawAll() +{ +} + +void HeliosScorerFirstStageBackStripTrapezoid::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Front Strip position Scorer for Annular geometry +HeliosScorerFirstStageFrontStripAnnular::HeliosScorerFirstStageFrontStripAnnular(G4String name, G4int depth, G4double StripPlaneSize, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_StripPlaneSize = StripPlaneSize ; + m_NumberOfStrip = NumberOfStrip ; +} + +HeliosScorerFirstStageFrontStripAnnular::~HeliosScorerFirstStageFrontStripAnnular() +{ +} + +G4bool HeliosScorerFirstStageFrontStripAnnular::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosAnnular"); + + // get front strip + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = m_StripPlaneSize / m_NumberOfStrip; + + G4double temp = (POS(0) + m_StripPlaneSize / 2.) / StripPitch ; + G4double X = int(temp) + 1 ; + //Rare case where particle is close to edge of silicon plan + if (X == 129) X = 128; + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, X); + return TRUE; +} + +void HeliosScorerFirstStageFrontStripAnnular::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageFrontStripAnnular::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageFrontStripAnnular::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageFrontStripAnnular::DrawAll() +{ +} + +void HeliosScorerFirstStageFrontStripAnnular::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// FirstStage Back Strip position Scorer for Annular geometry +HeliosScorerFirstStageBackStripAnnular::HeliosScorerFirstStageBackStripAnnular(G4String name, G4int depth, G4double StripPlaneSize, G4int NumberOfStrip) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + m_StripPlaneSize = StripPlaneSize ; + m_NumberOfStrip = NumberOfStrip ; +} + +HeliosScorerFirstStageBackStripAnnular::~HeliosScorerFirstStageBackStripAnnular() +{ +} + +G4bool HeliosScorerFirstStageBackStripAnnular::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + // get detector number + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "HeliosAnnular"); + + // get back strip + G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); + POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); + + G4double StripPitch = m_StripPlaneSize / m_NumberOfStrip; + + G4double temp = (POS(1) + m_StripPlaneSize / 2.) / StripPitch ; + G4int temp2 = temp ; + G4double Y = temp2 + 1 ; + //Rare case where particle is close to edge of silicon plan + if (Y == 129) Y = 128; + + G4double edep = aStep->GetTotalEnergyDeposit(); + //if (edep < 100*keV) return FALSE; + if (edep < 100*eV) return FALSE; // modified by mMarc + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(DetNbr + index, Y); + return TRUE; +} + +void HeliosScorerFirstStageBackStripAnnular::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void HeliosScorerFirstStageBackStripAnnular::EndOfEvent(G4HCofThisEvent*) +{ +} + +void HeliosScorerFirstStageBackStripAnnular::Clear() +{ + EvtMap->clear(); +} + +void HeliosScorerFirstStageBackStripAnnular::DrawAll() +{ +} + +void HelioScorerFirstStageBackStripAnnular::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} +*/ diff --git a/NPSimulation/src/MyMagneticField.cc b/NPSimulation/src/MyMagneticField.cc new file mode 100644 index 000000000..c26e7adcf --- /dev/null +++ b/NPSimulation/src/MyMagneticField.cc @@ -0,0 +1,71 @@ + +/***************************************************************************** + * Copyright (C) 2009 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: Marc Labiche contact address: marc.labiche@stfc.ac.uk * + * * + * Creation Date : 31/01/12 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class defines the Magnetic field for Helios * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "MyMagneticField.hh" + +MyMagneticField::MyMagneticField(const G4ThreeVector& FieldVector) +{ + + fFieldComponents[0]= FieldVector.x()*tesla; + fFieldComponents[1]= FieldVector.y()*tesla; + fFieldComponents[2]= FieldVector.z()*tesla; + + //Bz = 2.0*tesla; + //rmax = std::sqrt(625.)*cm; + rmax = 50.*cm; + zmin = 5.*cm; + zmax = 200.*cm; + + //G4cout << "rmax=" << rmax << G4endl; + +} + +MyMagneticField::~MyMagneticField() +{;} + +void MyMagneticField::GetFieldValue(const double Point[3],double *Bfield) const +{ + + // Uniform magnetic field along z axis confine in a cylinder + Bfield[0] = fFieldComponents[0]; + Bfield[1] = fFieldComponents[1]; + + //if(std::abs(Point[2])<zmax && (sqr(Point[0])+sqr(Point[1]))<rmax_sq) + //if(std::abs(Point[2])<zmax && std::abs(Point[2])>zmin && (std::sqrt(Point[0]*Point[0] + Point[1]*Point[1]))<rmax_sq) + if(std::abs(Point[2])<zmax && (std::sqrt(Point[0]*Point[0] + Point[1]*Point[1]))<rmax) + { Bfield[2] = fFieldComponents[2]; } + else + { Bfield[2] = 0.; } + +} + + + +void MyMagneticField::SetFieldValue(const G4ThreeVector &NewBfieldValue) +{ + + // Uniform magnetic field along z axis confine in a cylinder + fFieldComponents[0] = NewBfieldValue.x(); + fFieldComponents[1] = NewBfieldValue.y(); + fFieldComponents[2] = NewBfieldValue.z(); + +} -- GitLab