From 2697dbffa1caaa494fc41b7d4828dd7cde818b0f Mon Sep 17 00:00:00 2001 From: deserevi <deserevi@nptool> Date: Fri, 1 Oct 2010 08:45:23 +0000 Subject: [PATCH] * Fix bug in read sensitive part for TInteractionCoordinates + The loop was not done properly (wrong size) + Changes have been made for Gaspard + *** Similar changes should be done for other detectors *** * Gaspard and Paris simulations behaves as before merging --- .../gaspardTestSpheric.detector | 4 +-- Inputs/EventGenerator/gamma.source | 6 ++-- NPAnalysis/Gaspard/include/ObjectManager.hh | 8 ++--- NPAnalysis/Gaspard/src/Analysis.cc | 32 +++++++++---------- NPSimulation/src/GaspardTrackerAnnular.cc | 12 ++++--- NPSimulation/src/GaspardTrackerDummyShape.cc | 19 +++++------ NPSimulation/src/GaspardTrackerSquare.cc | 12 ++++--- NPSimulation/src/GaspardTrackerTrapezoid.cc | 12 ++++--- NPSimulation/src/PhysicsList.cc | 6 ++-- NPSimulation/vis.mac | 7 ++-- 10 files changed, 60 insertions(+), 58 deletions(-) diff --git a/Inputs/DetectorConfiguration/gaspardTestSpheric.detector b/Inputs/DetectorConfiguration/gaspardTestSpheric.detector index fb396b487..fb2cb947a 100644 --- a/Inputs/DetectorConfiguration/gaspardTestSpheric.detector +++ b/Inputs/DetectorConfiguration/gaspardTestSpheric.detector @@ -8,10 +8,10 @@ GeneralTarget %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Target - THICKNESS= 1.0 + THICKNESS= 10.0 ANGLE= 0 RADIUS= 12 - MATERIAL= CH2 + MATERIAL= CD2 NBLAYERS= 50 X= 0 Y= 0 diff --git a/Inputs/EventGenerator/gamma.source b/Inputs/EventGenerator/gamma.source index 7ec9ce537..877c06c32 100644 --- a/Inputs/EventGenerator/gamma.source +++ b/Inputs/EventGenerator/gamma.source @@ -5,9 +5,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic EnergyLow= .1 - EnergyHigh= .1 - HalfOpenAngleMin= 90 - HalfOpenAngleMax= 90 + EnergyHigh= 5. + HalfOpenAngleMin= 0 + HalfOpenAngleMax= 180 x0= 0 y0= 0 z0= 0 diff --git a/NPAnalysis/Gaspard/include/ObjectManager.hh b/NPAnalysis/Gaspard/include/ObjectManager.hh index 467f5e45e..f9e4038e8 100644 --- a/NPAnalysis/Gaspard/include/ObjectManager.hh +++ b/NPAnalysis/Gaspard/include/ObjectManager.hh @@ -108,8 +108,8 @@ namespace ENERGYLOSS // EnergyLoss LightTargetCD2 = EnergyLoss("proton_cd2.txt", 100, 1, 1); // LISE++ // For 132Sn(d,p) // CD2 -// EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", "G4Table", 1000); // G4 -// EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4 + EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", "G4Table", 1000); // G4 + EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4 // solid D2 // EnergyLoss LightTarget = EnergyLoss("proton_D2_solid.G4table", "G4Table", 1000); // G4 // EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_D2_solid.G4table", "G4Table", 1000); // G4 @@ -118,8 +118,8 @@ namespace ENERGYLOSS // EnergyLoss LightTarget = EnergyLoss("triton_CH2.G4table", "G4Table", 1000); // G4 // EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_CH2.G4table", "G4Table", 1000); // G4 // solid H2 - EnergyLoss LightTarget = EnergyLoss("triton_H2_solid.G4table", "G4Table", 1000); // G4 - EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_H2_solid.G4table", "G4Table", 1000); // G4 +// EnergyLoss LightTarget = EnergyLoss("triton_H2_solid.G4table", "G4Table", 1000); // G4 +// EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_H2_solid.G4table", "G4Table", 1000); // G4 // For 132Sn(d,t) // EnergyLoss LightTarget = EnergyLoss("triton_CD2.G4table", "G4Table", 1000); // G4 // EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4 diff --git a/NPAnalysis/Gaspard/src/Analysis.cc b/NPAnalysis/Gaspard/src/Analysis.cc index 6d46e9a27..b8dcba8a1 100644 --- a/NPAnalysis/Gaspard/src/Analysis.cc +++ b/NPAnalysis/Gaspard/src/Analysis.cc @@ -41,14 +41,8 @@ int main(int argc,char** argv) // Set energy beam at target middle myReaction->SetBeamEnergy(BeamEnergy); - // nominal beam energy - Double_t BeamEnergyNominal = myReaction->GetBeamEnergy() * MeV; - cout << BeamEnergyNominal << endl; - // slow beam at target middle - Double_t BeamEnergy = BeamEnergyNominal - BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0); - cout << BeamEnergy << endl; - // set energy beam at target middle - myReaction->SetBeamEnergy(BeamEnergy); + // Print target thickness + cout << myDetector->GetTargetThickness() << endl; // Attach more branch to the output double Ex = 0 ; double ExNoStrips = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ; @@ -77,6 +71,7 @@ int main(int argc,char** argv) // Analysis is here! int nentries = chain->GetEntries(); +// nentries = 106; cout << "Number of entries to be analysed: " << nentries << endl; // Default initialization @@ -99,6 +94,7 @@ int main(int argc,char** argv) // Get total energy double E = GPDTrack->GetEnergyDeposit(); +// cout << i << " " << E << endl; // if there is a hit in the detector array, treat it. double Theta, ThetaStrip, angle, ThetaCM; @@ -111,9 +107,13 @@ int main(int argc,char** argv) // Get exact scattering angle from TInteractionCoordinates object // Theta = interCoord->GetDetectedAngleTheta(0) * deg; +// cout << interCoord << endl; +// interCoord->Dump(); +// cout << i << " mult: " << interCoord->GetDetectedMultiplicity() << endl; DetecX = interCoord->GetDetectedPositionX(0); DetecY = interCoord->GetDetectedPositionY(0); DetecZ = interCoord->GetDetectedPositionZ(0); +// cout << DetecX << " " << DetecY << " " << DetecZ << endl; TVector3 Detec(DetecX, DetecY, DetecZ); // Get interaction position in detector @@ -141,12 +141,12 @@ int main(int argc,char** argv) // ThetaStrip = ThetaCalculation(HitDirection, BeamDirection); // Theta = ThetaCalculation(Detec - TVector3(XTarget, YTarget, 0), BeamDirection); // Calculate scattering angle w.r.t. beam (finite spatial resolution) -/* double resol = 800; // in micrometer - angle = gene->Rndm() * 2*3.14; - r = fabs(gene->Gaus(0, resol)) * micrometer; - ThetaStrip = ThetaCalculation(A - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection); - Theta = ThetaCalculation(Detec - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection); -*/ +// double resol = 800; // in micrometer +// angle = gene->Rndm() * 2*3.14; +// r = fabs(gene->Gaus(0, resol)) * micrometer; +// ThetaStrip = ThetaCalculation(A - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection); +// Theta = ThetaCalculation(Detec - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection); +// // Correct for energy loss in the target E = LightTarget.EvaluateInitialEnergy(E, myDetector->GetTargetThickness()/2 * micrometer, ThetaStrip); @@ -155,8 +155,8 @@ int main(int argc,char** argv) // if (Theta/deg < 60 && ThetaCM/deg < 90) { // if (Theta/deg > 35 && Theta/deg < 45 && E/MeV < 17) { // if (Theta/deg < 45) { - if (E/MeV < 38) { // for (p,t) reaction -// if (Theta/deg > 90) { // for (d,p) reaction +// if (E/MeV < 38) { // for (p,t) reaction + if (Theta/deg > 90) { // for (d,p) reaction ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad); Ex = myReaction->ReconstructRelativistic(E, ThetaStrip); } diff --git a/NPSimulation/src/GaspardTrackerAnnular.cc b/NPSimulation/src/GaspardTrackerAnnular.cc index b6014d4d9..a9e6efb95 100644 --- a/NPSimulation/src/GaspardTrackerAnnular.cc +++ b/NPSimulation/src/GaspardTrackerAnnular.cc @@ -16,6 +16,8 @@ *---------------------------------------------------------------------------* * Comment: * * + 12/10/09: Change scorer scheme (N. de Sereville) * + * + 01/10/10: Fix bug with TInteractionCoordinate map size in Read * + * Sensitive (N. de Sereville) * * * * * *****************************************************************************/ @@ -669,7 +671,7 @@ void GaspardTrackerAnnular::ReadSensitive(const G4Event* event) // Pos X Pos_X_itr = PosXHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosXHitMap->entries(); h++) { G4int PosXTrackID = Pos_X_itr->first - N ; G4double PosX = *(Pos_X_itr->second) ; if (PosXTrackID == NTrackID) { @@ -680,7 +682,7 @@ void GaspardTrackerAnnular::ReadSensitive(const G4Event* event) // Pos Y Pos_Y_itr = PosYHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosYHitMap->entries(); h++) { G4int PosYTrackID = Pos_Y_itr->first - N ; G4double PosY = *(Pos_Y_itr->second) ; if (PosYTrackID == NTrackID) { @@ -691,7 +693,7 @@ void GaspardTrackerAnnular::ReadSensitive(const G4Event* event) // Pos Z Pos_Z_itr = PosZHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosZHitMap->entries(); h++) { G4int PosZTrackID = Pos_Z_itr->first - N ; G4double PosZ = *(Pos_Z_itr->second) ; if (PosZTrackID == NTrackID) { @@ -702,7 +704,7 @@ void GaspardTrackerAnnular::ReadSensitive(const G4Event* event) // Angle Theta Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngThetaHitMap->entries(); h++) { G4int AngThetaTrackID = Ang_Theta_itr->first - N ; G4double AngTheta = *(Ang_Theta_itr->second) ; if (AngThetaTrackID == NTrackID) { @@ -713,7 +715,7 @@ void GaspardTrackerAnnular::ReadSensitive(const G4Event* event) // Angle Phi Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngPhiHitMap->entries(); h++) { G4int AngPhiTrackID = Ang_Phi_itr->first - N; G4double AngPhi = *(Ang_Phi_itr->second); if (AngPhiTrackID == NTrackID) { diff --git a/NPSimulation/src/GaspardTrackerDummyShape.cc b/NPSimulation/src/GaspardTrackerDummyShape.cc index 5c83ef469..1867e685e 100644 --- a/NPSimulation/src/GaspardTrackerDummyShape.cc +++ b/NPSimulation/src/GaspardTrackerDummyShape.cc @@ -20,6 +20,8 @@ * + 07/09/09: Fix bug for placing module with (r,theta,phi) method. * * (N. de Sereville) * * + 12/10/09: Change scorer scheme (N. de Sereville) * + * + 01/10/10: Fix bug with TInteractionCoordinate map size in Read * + * Sensitive (N. de Sereville) * * * * * *****************************************************************************/ @@ -758,7 +760,6 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) //G4cout << "sizeN:" << sizeN << G4endl; - if (sizeE != sizeT || sizeT != sizeX || sizeX != sizeY) { G4cout << "No match size Si Event Map: sE:" << sizeE << " sT:" << sizeT << " sX:" << sizeX << " sY:" << sizeY << endl ; @@ -770,11 +771,6 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) G4double 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->SetGPDTrkFirstStageFrontEDetectorNbr(m_index["DummyShape"] + N); @@ -832,7 +828,7 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) // Pos X Pos_X_itr = PosXHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosXHitMap->entries(); h++) { G4int PosXTrackID = Pos_X_itr->first - N ; G4double PosX = *(Pos_X_itr->second) ; if (PosXTrackID == NTrackID) { @@ -843,7 +839,7 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) // Pos Y Pos_Y_itr = PosYHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0 ; h < PosYHitMap->entries(); h++) { G4int PosYTrackID = Pos_Y_itr->first - N ; G4double PosY = *(Pos_Y_itr->second) ; if (PosYTrackID == NTrackID) { @@ -854,7 +850,7 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) // Pos Z Pos_Z_itr = PosZHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0 ; h < PosZHitMap->entries(); h++) { G4int PosZTrackID = Pos_Z_itr->first - N ; G4double PosZ = *(Pos_Z_itr->second) ; if (PosZTrackID == NTrackID) { @@ -865,7 +861,7 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) // Angle Theta Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngThetaHitMap->entries(); h++) { G4int AngThetaTrackID = Ang_Theta_itr->first - N ; G4double AngTheta = *(Ang_Theta_itr->second) ; if (AngThetaTrackID == NTrackID) { @@ -876,7 +872,7 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) // Angle Phi Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngPhiHitMap->entries(); h++) { G4int AngPhiTrackID = Ang_Phi_itr->first - N ; G4double AngPhi = *(Ang_Phi_itr->second) ; if (AngPhiTrackID == NTrackID) { @@ -922,6 +918,7 @@ void GaspardTrackerDummyShape::ReadSensitive(const G4Event* event) DetectorNumber_itr++; } + // clear map for next event DetectorNumberHitMap -> clear(); EnergyHitMap -> clear(); diff --git a/NPSimulation/src/GaspardTrackerSquare.cc b/NPSimulation/src/GaspardTrackerSquare.cc index 5fbeb423a..71e021a8e 100644 --- a/NPSimulation/src/GaspardTrackerSquare.cc +++ b/NPSimulation/src/GaspardTrackerSquare.cc @@ -18,6 +18,8 @@ * + 07/09/09: Fix bug for placing module with (r,theta,phi) method. * * (N. de Sereville) * * + 12/10/09: Change scorer scheme (N. de Sereville) * + * + 01/10/10: Fix bug with TInteractionCoordinate map size in Read * + * Sensitive (N. de Sereville) * * * * * *****************************************************************************/ @@ -1007,7 +1009,7 @@ void GaspardTrackerSquare::ReadSensitive(const G4Event* event) // Pos X Pos_X_itr = PosXHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosXHitMap->entries(); h++) { G4int PosXTrackID = Pos_X_itr->first - N ; G4double PosX = *(Pos_X_itr->second) ; if (PosXTrackID == NTrackID) { @@ -1018,7 +1020,7 @@ void GaspardTrackerSquare::ReadSensitive(const G4Event* event) // Pos Y Pos_Y_itr = PosYHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosYHitMap->entries(); h++) { G4int PosYTrackID = Pos_Y_itr->first - N ; G4double PosY = *(Pos_Y_itr->second) ; if (PosYTrackID == NTrackID) { @@ -1029,7 +1031,7 @@ void GaspardTrackerSquare::ReadSensitive(const G4Event* event) // Pos Z Pos_Z_itr = PosZHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosZHitMap->entries(); h++) { G4int PosZTrackID = Pos_Z_itr->first - N ; G4double PosZ = *(Pos_Z_itr->second) ; if (PosZTrackID == NTrackID) { @@ -1040,7 +1042,7 @@ void GaspardTrackerSquare::ReadSensitive(const G4Event* event) // Angle Theta Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngThetaHitMap->entries(); h++) { G4int AngThetaTrackID = Ang_Theta_itr->first - N ; G4double AngTheta = *(Ang_Theta_itr->second) ; if (AngThetaTrackID == NTrackID) { @@ -1051,7 +1053,7 @@ void GaspardTrackerSquare::ReadSensitive(const G4Event* event) // Angle Phi Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngPhiHitMap->entries(); h++) { G4int AngPhiTrackID = Ang_Phi_itr->first - N ; G4double AngPhi = *(Ang_Phi_itr->second) ; if (AngPhiTrackID == NTrackID) { diff --git a/NPSimulation/src/GaspardTrackerTrapezoid.cc b/NPSimulation/src/GaspardTrackerTrapezoid.cc index a0b6abe10..1a501934a 100644 --- a/NPSimulation/src/GaspardTrackerTrapezoid.cc +++ b/NPSimulation/src/GaspardTrackerTrapezoid.cc @@ -16,6 +16,8 @@ *---------------------------------------------------------------------------* * Comment: * * + 12/10/09: Change scorer scheme (N. de Sereville) * + * + 01/10/10: Fix bug with TInteractionCoordinate map size in Read * + * Sensitive (N. de Sereville) * * * * * *****************************************************************************/ @@ -891,7 +893,7 @@ void GaspardTrackerTrapezoid::ReadSensitive(const G4Event* event) // Pos X Pos_X_itr = PosXHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosXHitMap->entries(); h++) { G4int PosXTrackID = Pos_X_itr->first - N ; G4double PosX = *(Pos_X_itr->second) ; if (PosXTrackID == NTrackID) { @@ -902,7 +904,7 @@ void GaspardTrackerTrapezoid::ReadSensitive(const G4Event* event) // Pos Y Pos_Y_itr = PosYHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosYHitMap->entries(); h++) { G4int PosYTrackID = Pos_Y_itr->first - N ; G4double PosY = *(Pos_Y_itr->second) ; if (PosYTrackID == NTrackID) { @@ -913,7 +915,7 @@ void GaspardTrackerTrapezoid::ReadSensitive(const G4Event* event) // Pos Z Pos_Z_itr = PosZHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < PosZHitMap->entries(); h++) { G4int PosZTrackID = Pos_Z_itr->first - N ; G4double PosZ = *(Pos_Z_itr->second) ; if (PosZTrackID == NTrackID) { @@ -924,7 +926,7 @@ void GaspardTrackerTrapezoid::ReadSensitive(const G4Event* event) // Angle Theta Ang_Theta_itr = AngThetaHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngThetaHitMap->entries(); h++) { G4int AngThetaTrackID = Ang_Theta_itr->first - N ; G4double AngTheta = *(Ang_Theta_itr->second) ; if (AngThetaTrackID == NTrackID) { @@ -935,7 +937,7 @@ void GaspardTrackerTrapezoid::ReadSensitive(const G4Event* event) // Angle Phi Ang_Phi_itr = AngPhiHitMap->GetMap()->begin(); - for (G4int h = 0 ; h < sizeX ; h++) { + for (G4int h = 0; h < AngPhiHitMap->entries(); h++) { G4int AngPhiTrackID = Ang_Phi_itr->first - N ; G4double AngPhi = *(Ang_Phi_itr->second) ; if (AngPhiTrackID == NTrackID) { diff --git a/NPSimulation/src/PhysicsList.cc b/NPSimulation/src/PhysicsList.cc index a31200545..464c52a63 100644 --- a/NPSimulation/src/PhysicsList.cc +++ b/NPSimulation/src/PhysicsList.cc @@ -310,9 +310,9 @@ void PhysicsList::SetCuts() SetCutsWithDefault(); // for gamma-rays - SetCutValue(0.01*mm, "gamma"); - SetCutValue(0.01*mm, "e-"); - SetCutValue(0.01*mm, "e+"); +// SetCutValue(0.01*mm, "gamma"); +// SetCutValue(0.01*mm, "e-"); +// SetCutValue(0.01*mm, "e+"); // Retrieve verbose level SetVerboseLevel(temp); diff --git a/NPSimulation/vis.mac b/NPSimulation/vis.mac index de2be93c3..a5d707fb4 100644 --- a/NPSimulation/vis.mac +++ b/NPSimulation/vis.mac @@ -1,10 +1,10 @@ #verbose level -/control/verbose 1 -/material/verbose 1 +/control/verbose 0 +/material/verbose 0 /tracking/verbose 0 /process/verbose 0 /event/verbose 0 -/run/verbose 1 +/run/verbose 0 # choose a graphic system ##/vis/open OGLIX @@ -25,4 +25,3 @@ # run event #/run/beamOn 0 #/run/beamOn 10000 -/run/beamOn 5000 -- GitLab