From 1d8d4770cef95ef7dbee17c2f96dfb5a64d06f46 Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Fri, 18 Sep 2009 08:53:51 +0000
Subject: [PATCH] * Upgrade Gaspard analysis - Take into account energy loss
in the target for reconstruction of the excitation energy - change
makefile to recompile when header files are changed
* Clean some part of the code to remove warnings at compilation time
---
Inputs/EventGenerator/132Sndp.reaction | 4 +--
NPAnalysis/Gaspard/include/ObjectManager.hh | 14 ++++------
NPAnalysis/Gaspard/src/Analysis.cc | 10 ++++---
NPAnalysis/Gaspard/src/GNUmakefile | 3 ++-
NPLib/GASPARD/GaspardTracker.cxx | 6 +++++
NPLib/GASPARD/GaspardTracker.h | 2 +-
NPLib/GASPARD/TGaspardTrackerPhysics.cxx | 7 ++---
NPLib/VDetector/VDetector.h | 2 +-
.../include/GaspardTrackerDummyShape.hh | 26 ++++++++-----------
NPSimulation/include/Target.hh | 7 +++++
NPSimulation/src/AnnularS1.cc | 12 ++++-----
NPSimulation/src/GaspardTrackerAnnular.cc | 12 ++++-----
NPSimulation/src/GaspardTrackerSquare.cc | 6 ++---
NPSimulation/src/GaspardTrackerTrapezoid.cc | 12 ++++-----
NPSimulation/src/GeneralScorers.cc | 2 +-
NPSimulation/src/MUST2Array.cc | 1 -
NPSimulation/src/PlasticScorer.cc | 6 ++---
TODO | 3 +--
18 files changed, 72 insertions(+), 63 deletions(-)
diff --git a/Inputs/EventGenerator/132Sndp.reaction b/Inputs/EventGenerator/132Sndp.reaction
index 0618fe3f6..0ff5e1462 100644
--- a/Inputs/EventGenerator/132Sndp.reaction
+++ b/Inputs/EventGenerator/132Sndp.reaction
@@ -12,8 +12,8 @@ Transfert
BeamEnergySpread= 0
SigmaX= 0.851
SigmaY= 0.851
- SigmaThetaX= 1
- SigmaPhiY= 1
+ SigmaThetaX= 0
+ SigmaPhiY= 0
CrossSectionPath= sn132dp_gs_10AMeV.txt
ShootLight= 1
ShootHeavy= 0
diff --git a/NPAnalysis/Gaspard/include/ObjectManager.hh b/NPAnalysis/Gaspard/include/ObjectManager.hh
index 87249677a..49ee8e46c 100644
--- a/NPAnalysis/Gaspard/include/ObjectManager.hh
+++ b/NPAnalysis/Gaspard/include/ObjectManager.hh
@@ -102,15 +102,11 @@ using namespace CUT ;
#include "NPEnergyLoss.h"
using namespace NPL ;
namespace ENERGYLOSS
- {
-
- // Declare your Energy loss here :
- /* EnergyLoss ProtonTarget = EnergyLoss ( "CD2.txt" ,
- 100 ,
- 1 );
- */
- }
-
+{
+ // Declare your Energy loss here
+ EnergyLoss DeutonTargetCD2 = EnergyLoss("deuton_cd2.txt", 100, 1, 2);
+}
+
using namespace ENERGYLOSS ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/NPAnalysis/Gaspard/src/Analysis.cc b/NPAnalysis/Gaspard/src/Analysis.cc
index 8e4125f0c..f7c83fdf9 100644
--- a/NPAnalysis/Gaspard/src/Analysis.cc
+++ b/NPAnalysis/Gaspard/src/Analysis.cc
@@ -31,8 +31,9 @@ int main(int argc,char** argv)
myDetector->ReadConfigurationFile(detectorfileName);
// Attach more branch to the output
- double Ex = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
+ double Ex = 0 ; double ExNoStrips = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&Ex,"Ex/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergyNoStrips",&ExNoStrips,"ExNoStrips/D") ;
RootOutput::getInstance()->GetTree()->Branch("E",&EE,"EE/D") ;
RootOutput::getInstance()->GetTree()->Branch("A",&TT,"TT/D") ;
RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ;
@@ -76,10 +77,13 @@ int main(int argc,char** argv)
// Calculate scattering angle
ThetaStrip = ThetaCalculation (A ,TVector3(0,0,1));
+ // Correct for energy loss in the target
+ E = DeutonTargetCD2.EvaluateInitialEnergy(E, 4.85*micrometer, ThetaStrip);
+
// Calculate excitation energy
if (Theta/deg > 90) {
-// Ex = myReaction->ReconstructRelativistic(E, Theta / rad);
- Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
+ ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad);
+ Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
}
else Ex = -200;
}
diff --git a/NPAnalysis/Gaspard/src/GNUmakefile b/NPAnalysis/Gaspard/src/GNUmakefile
index c18b108aa..ee550bf80 100644
--- a/NPAnalysis/Gaspard/src/GNUmakefile
+++ b/NPAnalysis/Gaspard/src/GNUmakefile
@@ -24,12 +24,13 @@ LDFLAGS+= -L$(NPLIB)/lib -lVDetector -lIORoot -lReaction -lEnergyLoss \
LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB)
SRC= $(wildcard *.cc)
+INC= $(wildcard $(NPAINCLUDES)/*.hh)
OBJ=$(SRC:.cc=.o)
#all:$(EXEC)
# @$(CPP) -o $@ -c $< $(CXXFLAGS)
-Analysis:$(OBJ)
+Analysis:$(OBJ) $(INC)
@$(CPP) $(CXXFLAGS) -o $@ $^ $(LDFLAGS)
mv Analysis ../Analysis
diff --git a/NPLib/GASPARD/GaspardTracker.cxx b/NPLib/GASPARD/GaspardTracker.cxx
index 32874e9e3..ee430b6cd 100644
--- a/NPLib/GASPARD/GaspardTracker.cxx
+++ b/NPLib/GASPARD/GaspardTracker.cxx
@@ -536,6 +536,9 @@ void GaspardTracker::AddModuleSquare(TVector3 C_X1_Y1,
{
m_NumberOfModule++;
+ // remove warning using C_X128_Y128
+ C_X128_Y128.Unit();
+
// Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
TVector3 U = C_X128_Y1 - C_X1_Y1;
U = U.Unit();
@@ -692,6 +695,9 @@ void GaspardTracker::AddModuleDummyShape(TVector3 C_X1_Y1,
{
m_NumberOfModule++;
+ // remove warning using C_X128_Y128
+ C_X128_Y128.Unit();
+
// Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
TVector3 U = C_X128_Y1 - C_X1_Y1;
U = U.Unit();
diff --git a/NPLib/GASPARD/GaspardTracker.h b/NPLib/GASPARD/GaspardTracker.h
index e0d25a08b..5de55bd49 100644
--- a/NPLib/GASPARD/GaspardTracker.h
+++ b/NPLib/GASPARD/GaspardTracker.h
@@ -37,7 +37,7 @@ class GaspardTracker : public NPA::VDetector
{
public:
GaspardTracker();
- ~GaspardTracker();
+ virtual ~GaspardTracker();
public:
/////////////////////////////////////
diff --git a/NPLib/GASPARD/TGaspardTrackerPhysics.cxx b/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
index 921d30947..ae48b30d5 100644
--- a/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
+++ b/NPLib/GASPARD/TGaspardTrackerPhysics.cxx
@@ -65,11 +65,12 @@ void TGaspardTrackerPhysics::BuildPhysicalEvent(TGaspardTrackerData* Data)
bool Check_FirstStage = false ;bool Check_SecondStage = false ; bool Check_ThirdStage = false ;
// Thresholds
+/*
double FirstStage_Front_E_Threshold = 0; double FirstStage_Front_T_Threshold = 0;
double FirstStage_Back_E_Threshold = 0; double FirstStage_Back_T_Threshold = 0;
double SecondStage_E_Threshold = 0; double SecondStage_T_Threshold = 0;
double ThirdStage_E_Threshold = 0; double ThirdStage_T_Threshold = 0;
-
+*/
// calculate multipicity in the first stage
int multXE = Data->GetGPDTrkFirstStageFrontEMult();
int multYE = Data->GetGPDTrkFirstStageBackEMult();
@@ -133,8 +134,8 @@ void TGaspardTrackerPhysics::BuildPhysicalEvent(TGaspardTrackerData* Data)
// get time from strips and store it
double TimeStripFront = Data->GetGPDTrkFirstStageFrontEEnergy(0);
double TimeStripBack = Data->GetGPDTrkFirstStageBackEEnergy(0);
- double TimeStrip = 0.5 * (EnergyStripFront + EnergyStripBack);
-// double TimeStrip = EnergyStripFront;
+ double TimeStrip = 0.5 * (TimeStripFront + TimeStripBack);
+// double TimeStrip = TimeStripFront;
// if (TimeStripBack > TimeStrip) TimeStrip = TimeStripBack;
FirstStage_T.push_back(TimeStrip);
diff --git a/NPLib/VDetector/VDetector.h b/NPLib/VDetector/VDetector.h
index c10a4f39c..20fd8d3f9 100644
--- a/NPLib/VDetector/VDetector.h
+++ b/NPLib/VDetector/VDetector.h
@@ -38,7 +38,7 @@ namespace NPA
// Default Constructor and destructor
VDetector() ;
- ~VDetector() ;
+ virtual ~VDetector() ;
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
diff --git a/NPSimulation/include/GaspardTrackerDummyShape.hh b/NPSimulation/include/GaspardTrackerDummyShape.hh
index 8285639eb..49eb5e54c 100644
--- a/NPSimulation/include/GaspardTrackerDummyShape.hh
+++ b/NPSimulation/include/GaspardTrackerDummyShape.hh
@@ -142,30 +142,26 @@ private:
namespace GPDDUMMYSHAPE
{
// Resolution
-// const G4double ResoFirstStage = 0 ;// = 52keV of Resolution // Unit is MeV/2.35
- const G4double ResoFirstStage = 0.022 ;// = 52keV of Resolution // Unit is MeV/2.35
- const G4double ResoSecondStage = 0 ;// = 130 keV of resolution // Unit is MeV/2.35
-// const G4double ResoSecondStage = 0.055 ;// = 130 keV of resolution // Unit is MeV/2.35
- const G4double ResoThirdStage = 0 ;// = 100 keV of resolution // Unit is MeV/2.35
-// const G4double ResoThirdStage = 0.043 ;// = 100 kev of resolution // Unit is MeV/2.35
- const G4double ResoTimeGpd = 0.212765957 ;// = 500ps // Unit is ns/2.35
+ const G4double ResoFirstStage = 0.0213; // = 50 keV of Resolution // Unit is MeV/2.35
+ const G4double ResoSecondStage = 0.0213; // = 50 keV of resolution // Unit is MeV/2.35
+ const G4double ResoThirdStage = 0.0213; // = 50 keV of resolution // Unit is MeV/2.35
+ const G4double ResoTimeGpd = 0.212765957;// = 500ps // Unit is ns/2.35
// Geometry for the mother volume containing the different layers of your dummy shape module
const G4double FaceFront = 5.1*cm;
const G4double FaceBack = 5.1*cm;
-// const G4double Length = 1.5*cm;
-// const G4double InterStageDistance = 5*mm;
+ const G4double Length = 1.5*cm;
+ const G4double InterStageDistance = 5*mm;
// for testing the excitation energy reconstruction
- const G4double Length = 4*cm;
- const G4double InterStageDistance = 15*mm;
+// const G4double Length = 4*cm;
+// const G4double InterStageDistance = 15*mm;
// First stage
const G4double FirstStageFace = 5.0*cm;
-// const G4double FirstStageThickness = 300*micrometer ;
+ const G4double FirstStageThickness = 300*micrometer;
// for testing the excitation energy reconstruction
- const G4double FirstStageThickness = 1.3*cm;
-// const G4int NumberOfStrips = 128;
- const G4int NumberOfStrips = 25;
+// const G4double FirstStageThickness = 1.3*cm;
+ const G4int NumberOfStrips = 25; // 2mm strip pitch
// Second stage
const G4double SecondStageFace = FirstStageFace;
diff --git a/NPSimulation/include/Target.hh b/NPSimulation/include/Target.hh
index dea5c2a46..5fcee069b 100644
--- a/NPSimulation/include/Target.hh
+++ b/NPSimulation/include/Target.hh
@@ -67,6 +67,13 @@ public:
void ReadSensitive(const G4Event* event);
+public:
+ // method for debug purpose (still to be implemented)
+ // This method should check if the results of the beam interaction within the target
+ // (interaction coordinates) are well located inside the target volume
+ bool IsInsideTarget() {return false;};
+
+
public:
G4Material* GetMaterialFromLibrary(G4String MaterialName, G4double Temperature = 0, G4double Pressure = 0);
diff --git a/NPSimulation/src/AnnularS1.cc b/NPSimulation/src/AnnularS1.cc
index e24141ea2..e99912dd5 100644
--- a/NPSimulation/src/AnnularS1.cc
+++ b/NPSimulation/src/AnnularS1.cc
@@ -128,14 +128,14 @@ void AnnularS1::VolumeMaker(G4int TelescopeNumber ,
G4Material* Silicon = new G4Material("Si", z = 14., a, density);
// Al
- density = 2.702 * g / cm3;
- a = 26.98 * g / mole;
- G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
+// density = 2.702 * g / cm3;
+// a = 26.98 * g / mole;
+// G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
// Iron
- density = 7.874 * g / cm3;
- a = 55.847 * g / mole;
- G4Material* Iron = new G4Material("Iron", z = 26., a, density);
+// density = 7.874 * g / cm3;
+// a = 55.847 * g / mole;
+// G4Material* Iron = new G4Material("Iron", z = 26., a, density);
// CsI
density = 4.51 * g / cm3;
diff --git a/NPSimulation/src/GaspardTrackerAnnular.cc b/NPSimulation/src/GaspardTrackerAnnular.cc
index 9432cb7dc..ba566a4c2 100644
--- a/NPSimulation/src/GaspardTrackerAnnular.cc
+++ b/NPSimulation/src/GaspardTrackerAnnular.cc
@@ -135,14 +135,14 @@ void GaspardTrackerAnnular::VolumeMaker(G4int TelescopeNumber ,
G4Material* Silicon = new G4Material("Si", z = 14., a, density);
// Al
- density = 2.702 * g / cm3;
- a = 26.98 * g / mole;
- G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
+// density = 2.702 * g / cm3;
+// a = 26.98 * g / mole;
+// G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
// Iron
- density = 7.874 * g / cm3;
- a = 55.847 * g / mole;
- G4Material* Iron = new G4Material("Iron", z = 26., a, density);
+// density = 7.874 * g / cm3;
+// a = 55.847 * g / mole;
+// G4Material* Iron = new G4Material("Iron", z = 26., a, density);
// CsI
density = 4.51 * g / cm3;
diff --git a/NPSimulation/src/GaspardTrackerSquare.cc b/NPSimulation/src/GaspardTrackerSquare.cc
index 307c3ce24..39c2f6bde 100644
--- a/NPSimulation/src/GaspardTrackerSquare.cc
+++ b/NPSimulation/src/GaspardTrackerSquare.cc
@@ -188,9 +188,9 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber,
G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
// Iron
- density = 7.874 * g / cm3;
- a = 55.847 * g / mole;
- G4Material* Iron = new G4Material("Iron", z = 26., a, density);
+// density = 7.874 * g / cm3;
+// a = 55.847 * g / mole;
+// G4Material* Iron = new G4Material("Iron", z = 26., a, density);
// CsI
density = 4.51 * g / cm3;
diff --git a/NPSimulation/src/GaspardTrackerTrapezoid.cc b/NPSimulation/src/GaspardTrackerTrapezoid.cc
index 714a28c9d..b5f91d1f6 100644
--- a/NPSimulation/src/GaspardTrackerTrapezoid.cc
+++ b/NPSimulation/src/GaspardTrackerTrapezoid.cc
@@ -179,14 +179,14 @@ void GaspardTrackerTrapezoid::VolumeMaker(G4int TelescopeNumber ,
G4Material* Silicon = new G4Material("Si", z = 14., a, density);
// Al
- density = 2.702 * g / cm3;
- a = 26.98 * g / mole;
- G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
+// density = 2.702 * g / cm3;
+// a = 26.98 * g / mole;
+// G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density);
// Iron
- density = 7.874 * g / cm3;
- a = 55.847 * g / mole;
- G4Material* Iron = new G4Material("Iron", z = 26., a, density);
+// density = 7.874 * g / cm3;
+// a = 55.847 * g / mole;
+// G4Material* Iron = new G4Material("Iron", z = 26., a, density);
// CsI
density = 4.51 * g / cm3;
diff --git a/NPSimulation/src/GeneralScorers.cc b/NPSimulation/src/GeneralScorers.cc
index 9fe510c3b..f532c4495 100644
--- a/NPSimulation/src/GeneralScorers.cc
+++ b/NPSimulation/src/GeneralScorers.cc
@@ -429,7 +429,7 @@ G4bool PSDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
int numberOfCharacterInDetectorNumber = name.length() - (int)found ;
- for( int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
+ for(unsigned int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
nbr += name[i] ;
G4int DetNbr = atoi( nbr.c_str() ) ;
diff --git a/NPSimulation/src/MUST2Array.cc b/NPSimulation/src/MUST2Array.cc
index 3040dec99..5b854ff90 100644
--- a/NPSimulation/src/MUST2Array.cc
+++ b/NPSimulation/src/MUST2Array.cc
@@ -192,7 +192,6 @@ void MUST2Array::VolumeMaker(G4int TelescopeNumber ,
PVPBuffer = new G4PVPlacement(0,G4ThreeVector(0,0,0),logicDegrader,"Degrader",logicVacBox,false,0) ;
- /*
//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.
diff --git a/NPSimulation/src/PlasticScorer.cc b/NPSimulation/src/PlasticScorer.cc
index 4cf0904c2..ba2ca272b 100644
--- a/NPSimulation/src/PlasticScorer.cc
+++ b/NPSimulation/src/PlasticScorer.cc
@@ -45,7 +45,7 @@ G4bool PSEnergy::ProcessHits(G4Step* aStep, G4TouchableHistory*)
int numberOfCharacterInDetectorNumber = name.length() - (int)found ;
- for( int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
+ for(unsigned int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
nbr += name[i] ;
G4int DetNbr = atoi( nbr.c_str() ) ;
@@ -119,7 +119,7 @@ G4bool PSDetectorNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*)
int numberOfCharacterInDetectorNumber = name.length() - (int)found ;
- for( int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
+ for(unsigned int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
nbr += name[i] ;
G4int DetNbr = atoi( nbr.c_str() ) ;
@@ -193,7 +193,7 @@ G4bool PSTOF::ProcessHits(G4Step* aStep, G4TouchableHistory*)
int numberOfCharacterInDetectorNumber = name.length() - (int)found ;
- for( int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
+ for(unsigned int i = found ; i < found + numberOfCharacterInDetectorNumber ; i++ )
nbr += name[i] ;
G4int DetNbr = atoi( nbr.c_str() ) ;
diff --git a/TODO b/TODO
index f675b8725..3bf175015 100644
--- a/TODO
+++ b/TODO
@@ -31,7 +31,6 @@ TODO for NPTool:
(Adrien)
+ Possibility to tilt the target with a given angle and to deal with
the emittance
- + Shooting in the target within a dedicated method in VEventGenerator
+ IsInsideTarget method for debugging purposes
+ Add a dedicated class to deal with materials (see example from G4 tutorial)
(still under debate)
@@ -43,7 +42,6 @@ TODO for NPTool:
TODO for NPAnalysis:
--------------------
+ Add comment support for RunToTreat.txt (Adrien)
- + Check Makefiles dependency (mostly for Gaspard) (Nicolas)
+ Add calibration manager class from Julien (Adrien)
+ Use only one class for the analysis per detector instead of two currently
(ROOT feature in private member definition ( //!))
@@ -51,6 +49,7 @@ TODO for NPAnalysis:
TODO for AnnularS1 detector:
----------------------------
+ Update the way the scorers are dealt to the new way (no more loop in ReadSensitive)
+ (Nicolas)
TODO for Gaspard:
-----------------
--
GitLab