From 8a7fbc48befd04e73eace3734fbdf093e2c58548 Mon Sep 17 00:00:00 2001
From: deserevi <deserevi@nptool>
Date: Thu, 7 Jan 2010 09:52:03 +0000
Subject: [PATCH] * Add 208Pb in target material
* Fix bug when writing G4 energy loss tables in EventGeneratorTransfert
* Update Gaspard Analysis and controls macro
---
.../gaspardTestSpheric.detector | 2 +-
Inputs/EventGenerator/132Sndp.reaction | 3 ++-
Inputs/EventGenerator/134Snpt.reaction | 4 +--
NPAnalysis/Gaspard/include/ObjectManager.hh | 12 ++++-----
.../Gaspard/macros/DisplayInputCrossSection.C | 25 ++++++++++++++-----
NPAnalysis/Gaspard/src/Analysis.cc | 25 +++++++++++++------
NPAnalysis/macros/ControlSimu.C | 14 +++++++++++
NPSimulation/src/EventGeneratorTransfert.cc | 3 ---
NPSimulation/src/Target.cc | 8 ++++++
9 files changed, 70 insertions(+), 26 deletions(-)
diff --git a/Inputs/DetectorConfiguration/gaspardTestSpheric.detector b/Inputs/DetectorConfiguration/gaspardTestSpheric.detector
index 8b10e2ead..b09aa7a8e 100644
--- a/Inputs/DetectorConfiguration/gaspardTestSpheric.detector
+++ b/Inputs/DetectorConfiguration/gaspardTestSpheric.detector
@@ -11,7 +11,7 @@ Target
THICKNESS= 0.00001
ANGLE= 0
RADIUS= 12
- MATERIAL= CD2
+ MATERIAL= Pb208
NBLAYERS= 50
X= 0
Y= 0
diff --git a/Inputs/EventGenerator/132Sndp.reaction b/Inputs/EventGenerator/132Sndp.reaction
index 9614773cb..6507dc83e 100644
--- a/Inputs/EventGenerator/132Sndp.reaction
+++ b/Inputs/EventGenerator/132Sndp.reaction
@@ -14,7 +14,8 @@ Transfert
SigmaY= 0
SigmaThetaX= 0
SigmaPhiY= 0
- CrossSectionPath= 132Sndp_10A_MeV_3p3_ZR_FRC.lis
+ CrossSectionPath= 132Sndp_10A_MeV_2f7_ZR_FRC.lis
+% CrossSectionPath= flat.txt
ShootLight= 1
ShootHeavy= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/134Snpt.reaction b/Inputs/EventGenerator/134Snpt.reaction
index 6ff89b7fe..d9b48e32f 100644
--- a/Inputs/EventGenerator/134Snpt.reaction
+++ b/Inputs/EventGenerator/134Snpt.reaction
@@ -10,8 +10,8 @@ Transfert
ExcitationEnergy= 0.0
BeamEnergy= 1340
BeamEnergySpread= 0
- SigmaX= 0
- SigmaY= 0
+ SigmaX= 0.851
+ SigmaY= 0.851
SigmaThetaX= 0
SigmaPhiY= 0
CrossSectionPath= CS_Ep10MeV_sn134pt_gs_1h9demi.dat
diff --git a/NPAnalysis/Gaspard/include/ObjectManager.hh b/NPAnalysis/Gaspard/include/ObjectManager.hh
index a4ec95b93..d6c973089 100644
--- a/NPAnalysis/Gaspard/include/ObjectManager.hh
+++ b/NPAnalysis/Gaspard/include/ObjectManager.hh
@@ -23,7 +23,7 @@
#include <TFile.h>
#include <TLeaf.h>
#include <TVector3.h>
-#include <TRandom.h>
+#include <TRandom3.h>
// NPL
#include "TMust2Data.h"
@@ -107,11 +107,11 @@ namespace ENERGYLOSS
// Declare your Energy loss here
// EnergyLoss LightTargetCD2 = EnergyLoss("proton_cd2.txt", 100, 1, 1); // LISE++
// For 132Sn(d,p)
- EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", "G4Table", 1000); // G4
- EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4
- // For 132Sn(d,p)
-// EnergyLoss LightTarget = EnergyLoss("triton_CH2.G4table", "G4Table", 1000); // G4
-// EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_CH2.G4table", "G4Table", 1000); // G4
+// EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", "G4Table", 1000); // G4
+// EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4
+ // For 134Sn(p,t)
+ EnergyLoss LightTarget = EnergyLoss("triton_CH2.G4table", "G4Table", 1000); // G4
+ EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_CH2.G4table", "G4Table", 1000); // G4
}
using namespace ENERGYLOSS ;
diff --git a/NPAnalysis/Gaspard/macros/DisplayInputCrossSection.C b/NPAnalysis/Gaspard/macros/DisplayInputCrossSection.C
index cd475ef12..7d584e715 100644
--- a/NPAnalysis/Gaspard/macros/DisplayInputCrossSection.C
+++ b/NPAnalysis/Gaspard/macros/DisplayInputCrossSection.C
@@ -23,29 +23,41 @@ void DisplayInputCrossSection()
TGraph *grpt2 = ReadCrossSectionPT("CS_Ep15MeV_sn134pt_gs_1h9demi.dat");
TGraph *grpt3 = ReadCrossSectionPT("CS_Ep20MeV_sn134pt_gs_1h9demi.dat");
+ // Read cross-section 132Sn(d,d)
+ // Angel
+ TGraph *grdd = new TGraph(path + "132Sndd_10A_MeV_ruth.dat");
+
// Draw cross-sections
TCanvas *can = new TCanvas("can");
can->SetLogy();
can->Draw();
// TH2F *hframe = new TH2F("hframe", "^{2}H(^{132}Sn,p)^{133}Sn", 180, 0, 180, 100, 1e-2, 100);
- TH2F *hframe = new TH2F("hframe", "^{1}H(^{134}Sn,t)^{132}Sn_{g.s.}", 180, 0, 180, 100, 1e-8, 1e-5);
+// TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-2, 100);
+// TH2F *hframe = new TH2F("hframe", "^{1}H(^{134}Sn,t)^{132}Sn_{g.s.}", 180, 0, 180, 100, 1e-8, 1e-5);
+// TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-8, 1e-5);
+ TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-3, 10);
hframe->Draw();
hframe->SetXTitle("#Theta_{c.m.} [deg]");
- hframe->SetYTitle("d#sigma/d#Omega [mb/sr]");
- grpt1->SetLineColor(kRed); grpt1->Draw("l");
+// hframe->SetYTitle("d#sigma/d#Omega [mb/sr]");
+ hframe->SetYTitle("d#sigma/d#Omega / (d#sigma/d#Omega)_{R}");
+// hframe->SetYTitle("#propto d#sigma/d#Omega");
+/* grpt1->SetLineColor(kRed); grpt1->Draw("l");
grpt2->SetLineColor(kMagenta); grpt2->Draw("l");
- grpt3->SetLineColor(kBlue); grpt3->Draw("l");
+ grpt3->SetLineColor(kBlue); grpt3->Draw("l");*/
/* gr1->SetLineColor(kRed); gr1->SetLineStyle(2); gr1->Draw("l");
gr2->SetLineColor(kRed); gr2->Draw("l");
gr3->SetLineColor(kBlue); gr3->SetLineStyle(2); gr3->Draw("l");
gr4->SetLineColor(kBlue); gr4->Draw("l");*/
// gr5->Draw("l");
+ grdd->Draw("l");
// TLegend
TLegend *leg = new TLegend(0.50, 0.64, 0.82, 0.84);
+ leg->AddEntry(grdd, "10 MeV/u", "l");
+/* TLegend *leg = new TLegend(0.50, 0.64, 0.82, 0.84);
leg->AddEntry(grpt1, "1h9/2 10 MeV/u", "l");
leg->AddEntry(grpt2, "1h9/2 15 MeV/u", "l");
- leg->AddEntry(grpt3, "1h9/2 20 MeV/u", "l");
+ leg->AddEntry(grpt3, "1h9/2 20 MeV/u", "l");*/
/* TLegend *leg = new TLegend(0.16, 0.17, 0.48, 0.37);
leg->AddEntry(gr1, "3p3/2 5 MeV/u", "l");
leg->AddEntry(gr2, "2f7/2 5 MeV/u", "l");
@@ -114,7 +126,8 @@ TGraph* ReadCrossSectionPT(const char* fname)
Double_t angle, sigma, dum;
Int_t nlines = 0;
TGraph *gr = new TGraph();
- while (fich >> angle >> sigma >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum) {
+// while (fich >> angle >> sigma >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum) {
+ while (fich >> angle >> sigma) {
gr->SetPoint(nlines++, angle, sigma);
}
diff --git a/NPAnalysis/Gaspard/src/Analysis.cc b/NPAnalysis/Gaspard/src/Analysis.cc
index d88e08b9c..5d9445567 100644
--- a/NPAnalysis/Gaspard/src/Analysis.cc
+++ b/NPAnalysis/Gaspard/src/Analysis.cc
@@ -36,6 +36,7 @@ int main(int argc,char** argv)
cout << BeamEnergyNominal << endl;
// Slow beam at target middle
Double_t BeamEnergy = BeamEnergyNominal - BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
+// Double_t BeamEnergy = 1293.56 * MeV;
cout << BeamEnergy << endl;
// Set energy beam at target middle
myReaction->SetBeamEnergy(BeamEnergy);
@@ -78,6 +79,9 @@ int main(int argc,char** argv)
double BeamTheta = 0;
double BeamPhi = 0;
+ // random generator
+ TRandom3 *gene = new TRandom3();
+
// Loop on all events
for (int i = 0; i < nentries; i ++) {
if (i%10000 == 0 && i!=0) cout << "\r" << i << " analyzed events" << flush;
@@ -91,8 +95,9 @@ int main(int argc,char** argv)
double E = GPDTrack->GetEnergyDeposit();
// if there is a hit in the detector array, treat it.
- double Theta, ThetaStrip;
+ double Theta, ThetaStrip, angle;
double DetecX, DetecY, DetecZ;
+ double r;
TVector3 A;
if (E > -1000) {
// Get exact scattering angle from TInteractionCoordinates object
@@ -123,17 +128,23 @@ int main(int argc,char** argv)
// Calculate scattering angle w.r.t. optical beam axis (do not take into account beam position on target)
// ThetaStrip = ThetaCalculation(A, TVector3(0,0,1));
// Theta = ThetaCalculation(Detec, TVector3(0, 0, 1));
- // Calculate scattering angle w.r.t. beam
- ThetaStrip = ThetaCalculation(HitDirection, BeamDirection);
- Theta = ThetaCalculation(Detec - TVector3(XTarget, YTarget, 0), BeamDirection);
+ // Calculate scattering angle w.r.t. beam (ideal case)
+// 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);
// Correct for energy loss in the target
E = LightTarget.EvaluateInitialEnergy(E, myDetector->GetTargetThickness()/2 * micrometer, ThetaStrip);
// Calculate excitation energy
-// if (Theta/deg > 35 && Theta/deg < 85) {
-// if (Theta/deg < 85) {
- if (Theta/deg > 90) {
+// if (Theta/deg > 150 && Theta/deg < 180) {
+ if (Theta/deg < 45) {
+// if (Theta/deg > 90) {
ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad);
Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
}
diff --git a/NPAnalysis/macros/ControlSimu.C b/NPAnalysis/macros/ControlSimu.C
index 7854bd567..eec454239 100644
--- a/NPAnalysis/macros/ControlSimu.C
+++ b/NPAnalysis/macros/ControlSimu.C
@@ -140,4 +140,18 @@ void ControlSimu(const char * fname = "mySimul")
hEmittedPhiWF->Draw();
// hControlPhi->SetLineColor(kRed);
// hControlPhi->Draw("same");
+
+ // Display histograms
+ TCanvas *canvasxi32 = new TCanvas("canvas3", "Emitted particle properties", 300, 600);
+ canvas3->Divide(1,3);
+ canvas3->cd(1);
+ hEmittedThetaCM->SetXTitle("#Theta_{c.m.}");
+ hEmittedThetaCM->Draw();
+ canvas3->cd(2);
+ hEmittedThetaWF->Draw();
+ hEmittedThetaWF->SetXTitle("#Theta_{lab}");
+ canvas3->cd(3);
+ hEmittedETheta->SetXTitle("#Theta_{lab}");
+ hEmittedETheta->SetYTitle("E [MeV]");
+ hEmittedETheta->Draw();
}
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index d4bb00cc4..092961997 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -69,9 +69,6 @@ void EventGeneratorTransfert::SetTarget(Target* Target)
{
if (Target != 0) {
m_Target = Target;
- G4int LightZ = m_Reaction->GetNucleus3()->GetZ();
- G4int LightA = m_Reaction->GetNucleus3()->GetA();
- m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZ,LightA, 0.) ,0, m_BeamEnergy);
}
}
diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc
index 4537ef5d8..bce572dc1 100644
--- a/NPSimulation/src/Target.cc
+++ b/NPSimulation/src/Target.cc
@@ -181,6 +181,14 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe
return myMaterial;
}
+ else if (MaterialName == "Pb208") {
+ G4Element* Pb = new G4Element("Lead" , "Pb" , 82. , 207.2*g / mole);
+
+ G4Material* myMaterial = new G4Material("Pb208", 11.342*g / cm3, 1);
+ myMaterial->AddElement(Pb , 1);
+ return myMaterial;
+ }
+
else {
G4cout << "No Matching Material in the Target Library Default is Vacuum" << G4endl;
G4Element* N = new G4Element("Nitrogen", "N", 7., 14.01*g / mole);
--
GitLab