From f01b4b471cf294d92729bd83d8770c2b0864e7ab Mon Sep 17 00:00:00 2001
From: adrien-matta <a.matta@surrey.ac.uk>
Date: Tue, 8 Mar 2016 17:59:44 +0000
Subject: [PATCH] * Adding SharcEfficiency project + update on S1554 analysis
---
Inputs/DetectorConfiguration/S1107.detector | 82 +----
NPSimulation/ressources/macro/vis.mac | 2 +-
Outputs/Simulation/.gitignore | 0
Projects/S1554/Analysis.cxx | 10 +-
Projects/S1554/Analysis.h | 1 +
Projects/S1554/EXD.txt | 6 +
Projects/S1554/Reaction/28Mgdp_0.reaction | 40 ++
Projects/S1554/Reaction/28Mgdp_1.reaction | 40 ++
Projects/S1554/Reaction/28Mgdp_2.reaction | 40 ++
Projects/S1554/Reaction/28Mgdp_3.reaction | 40 ++
Projects/S1554/Reaction/28Mgdp_4.reaction | 40 ++
Projects/S1554/Reaction/28Mgdp_5.reaction | 40 ++
Projects/S1554/Reaction/28Mgdp_EXD.reaction | 30 ++
Projects/S1554/Reaction/28Mgdp_gs.reaction | 30 ++
Projects/S1554/macro/CS.cxx | 202 ++++++++++
Projects/S1554/macro/CheckGamma.cxx | 23 ++
Projects/S1554/macro/GammaRes.cxx | 34 ++
Projects/S1554/macro/Mg.cxx | 122 +++++++
Projects/S1554/macro/MgCS.cxx | 344 ++++++++++++++++++
Projects/S1554/macro/MgCS2.cxx | 238 ++++++++++++
Projects/S1554/macro/MgCSLab.cxx | 238 ++++++++++++
Projects/S1554/macro/MgCheck.cxx | 35 ++
Projects/S1554/macro/Mgdd.cxx | 59 +++
Projects/S1554/macro/Resolution.cxx | 61 ++++
Projects/S1554/macro/SM.cxx | 143 ++++++++
Projects/S1554/macro/ShowTrifoil.cxx | 23 ++
Projects/S1554/macro/Si.cxx | 57 +++
Projects/S1554/macro/SiCS.cxx | 211 +++++++++++
Projects/S1554/macro/SiCheck.cxx | 20 +
Projects/S1554/macro/SiCheck2.cxx | 22 ++
Projects/S1554/macro/SiCheck3.cxx | 19 +
Projects/S1554/macro/SiCheck4.cxx | 16 +
Projects/S1554/macro/Sidd.cxx | 44 +++
Projects/S1554/macro/makeMgCS.cxx | 89 +++++
Projects/S1554/macro/makeMgCSLab.cxx | 89 +++++
Projects/S1554/macro/pipo.cxx | 11 +
Projects/S1554/macro/test.cxx | 25 ++
Projects/S1554/mg29olapc.lsf | 94 +++++
Projects/SharcEfficiency/Analysis.cxx | 308 ++++++++++++++++
Projects/SharcEfficiency/Analysis.h | 117 ++++++
Projects/SharcEfficiency/CMakeLists.txt | 2 +
.../SharcEfficiency/configs/ConfigSharc.dat | 99 +++++
42 files changed, 3071 insertions(+), 75 deletions(-)
delete mode 100644 Outputs/Simulation/.gitignore
create mode 100644 Projects/S1554/EXD.txt
create mode 100644 Projects/S1554/Reaction/28Mgdp_0.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_1.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_2.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_3.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_4.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_5.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_EXD.reaction
create mode 100644 Projects/S1554/Reaction/28Mgdp_gs.reaction
create mode 100644 Projects/S1554/macro/CS.cxx
create mode 100644 Projects/S1554/macro/CheckGamma.cxx
create mode 100644 Projects/S1554/macro/GammaRes.cxx
create mode 100644 Projects/S1554/macro/Mg.cxx
create mode 100644 Projects/S1554/macro/MgCS.cxx
create mode 100644 Projects/S1554/macro/MgCS2.cxx
create mode 100644 Projects/S1554/macro/MgCSLab.cxx
create mode 100644 Projects/S1554/macro/MgCheck.cxx
create mode 100644 Projects/S1554/macro/Mgdd.cxx
create mode 100644 Projects/S1554/macro/Resolution.cxx
create mode 100644 Projects/S1554/macro/SM.cxx
create mode 100644 Projects/S1554/macro/ShowTrifoil.cxx
create mode 100644 Projects/S1554/macro/Si.cxx
create mode 100644 Projects/S1554/macro/SiCS.cxx
create mode 100644 Projects/S1554/macro/SiCheck.cxx
create mode 100644 Projects/S1554/macro/SiCheck2.cxx
create mode 100644 Projects/S1554/macro/SiCheck3.cxx
create mode 100644 Projects/S1554/macro/SiCheck4.cxx
create mode 100644 Projects/S1554/macro/Sidd.cxx
create mode 100644 Projects/S1554/macro/makeMgCS.cxx
create mode 100644 Projects/S1554/macro/makeMgCSLab.cxx
create mode 100644 Projects/S1554/macro/pipo.cxx
create mode 100644 Projects/S1554/macro/test.cxx
create mode 100644 Projects/S1554/mg29olapc.lsf
create mode 100644 Projects/SharcEfficiency/Analysis.cxx
create mode 100644 Projects/SharcEfficiency/Analysis.h
create mode 100644 Projects/SharcEfficiency/CMakeLists.txt
create mode 100644 Projects/SharcEfficiency/configs/ConfigSharc.dat
diff --git a/Inputs/DetectorConfiguration/S1107.detector b/Inputs/DetectorConfiguration/S1107.detector
index 5d70af4f6..1cbe728d0 100644
--- a/Inputs/DetectorConfiguration/S1107.detector
+++ b/Inputs/DetectorConfiguration/S1107.detector
@@ -58,74 +58,14 @@ Sharc
ThicknessPAD2= 1000
ThicknessPAD3= 1000
ThicknessPAD4= 1000
-%%%%%%%%%%%%%%%%%%%%%
-%Tigress Array
-Tigress
-%%%%%%%%%%%%%%%%%%%%
-%Corona(90) Ring
-%%%%%%%%%%%%%%%%%%%%
-TigressClover
- CloverId= 5
- R= 145
- Phi= 22.5
- Theta= 90
-TigressClover
- CloverId= 6
- R= 145
- Phi= 67.5
- Theta= 90
-TigressClover
- CloverId= 7
- R= 145
- Phi= 112.5
- Theta= 90
-TigressClover
- CloverId= 8
- R= 145
- Phi= 157.5
- Theta= 90
-TigressClover
- CloverId= 9
- R= 145
- Phi= 202.5
- Theta= 90
-TigressClover
- CloverId= 10
- R= 145
- Phi= 247.5
- Theta= 90
-TigressClover
- CloverId= 11
- R= 145
- Phi= 292.5
- Theta= 90
-TigressClover
- CloverId= 12
- R= 145
- Phi= 337.5
- Theta= 90
-%%%%%%%%%%%%%%%%%%%%
-%Backwards(145) Ring
-%%%%%%%%%%%%%%%%%%%%
-TigressClover
- CloverId= 13
- R= 145
- Phi= 22.5
- Theta= 135
-TigressClover
- CloverId= 14
- R= 145
- Phi= 112.5
- Theta= 135
-TigressClover
- CloverId= 15
- R= 145
- Phi= 202.5
- Theta= 135
-TigressClover
- CloverId= 16
- R= 145
- Phi= 292.5
- Theta= 135
-%%%%%%%%%%%%%%%%%%%%
-Trifoil
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Plastic
+ X= 0
+ Y= 0
+ Z= 400
+ Thickness= 0.010
+ Shape= Square
+ Height= 40
+ Width= 40
+ Scintillator= BC400
+ LeadThickness= 0.009
diff --git a/NPSimulation/ressources/macro/vis.mac b/NPSimulation/ressources/macro/vis.mac
index 45df6a2cb..bb1c0b15b 100644
--- a/NPSimulation/ressources/macro/vis.mac
+++ b/NPSimulation/ressources/macro/vis.mac
@@ -47,7 +47,7 @@
/vis/viewer/set/autoRefresh true
/vis/verbose 0
-/vis/viewer/set/background white
+/vis/viewer/set/background black
# print Option
#/vis/ogl/set/printMode vectored
diff --git a/Outputs/Simulation/.gitignore b/Outputs/Simulation/.gitignore
deleted file mode 100644
index e69de29bb..000000000
diff --git a/Projects/S1554/Analysis.cxx b/Projects/S1554/Analysis.cxx
index 6570ce369..4641555cb 100644
--- a/Projects/S1554/Analysis.cxx
+++ b/Projects/S1554/Analysis.cxx
@@ -41,19 +41,19 @@ void Analysis::Init(){
myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
string ligth,heavy;
- if(myReaction->GetNucleus3()->GetZ()==1 && myReaction->GetNucleus3()->GetA()==1)
+ if(myReaction->GetNucleus3().GetZ()==1 && myReaction->GetNucleus3().GetA()==1)
ligth = "proton";
else
ligth = "deuteron";
- if(myReaction->GetNucleus4()->GetZ()==12)
+ if(myReaction->GetNucleus4().GetZ()==12)
heavy = "Mg28";
else
heavy = "Si28";
Sharc = (TSharcPhysics*) m_DetectorManager -> GetDetector("Sharc");
- Tigress = (TTigressPhysics*) m_DetectorManager -> GetDetector("Tigress");
+ //Tigress = (TTigressPhysics*) m_DetectorManager -> GetDetector("Tigress");
LightCD2 = EnergyLoss(ligth+"_CD2.G4table","G4Table",10);
LightAl = EnergyLoss(ligth+"_Al.G4table","G4Table",10);
@@ -176,6 +176,8 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("RunNumber",&RunNumber,"RunNumber/I");
RootOutput::getInstance()->GetTree()->Branch("RunNumberMinor",&RunNumberMinor,"RunNumberMinor/I");
+ RootOutput::getInstance()->GetTree()->Branch("Run",&Run,"Run/I");
+
}
@@ -183,6 +185,8 @@ void Analysis::InitOutputBranch(){
void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchAddress("RunNumber",&RunNumber);
RootInput::getInstance()->GetChain()->SetBranchAddress("RunNumberMinor",&RunNumberMinor);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("Run",&Run);
+
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
diff --git a/Projects/S1554/Analysis.h b/Projects/S1554/Analysis.h
index ebe1b6d41..78b7664c2 100644
--- a/Projects/S1554/Analysis.h
+++ b/Projects/S1554/Analysis.h
@@ -75,6 +75,7 @@ TInitialConditions* myInit ;
int DetectorNumber ;
int RunNumber;
int RunNumberMinor;
+ int Run;
double ThetaNormalTarget;
double Energy ;
diff --git a/Projects/S1554/EXD.txt b/Projects/S1554/EXD.txt
new file mode 100644
index 000000000..ea82e214b
--- /dev/null
+++ b/Projects/S1554/EXD.txt
@@ -0,0 +1,6 @@
+0 1
+2 1
+4 1
+6 1
+8 1
+10 1
diff --git a/Projects/S1554/Reaction/28Mgdp_0.reaction b/Projects/S1554/Reaction/28Mgdp_0.reaction
new file mode 100644
index 000000000..7977482d3
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_0.reaction
@@ -0,0 +1,40 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 1.0
+ CrossSectionPath= flat.txt CSR0
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GammaDecay 29Mg
+ Cascade
+ BranchingRatio= 70
+ Energies= 1
+ DifferentialCrossSection= flat.txt Gamma29Mg
+ Cascade
+ BranchingRatio= 30
+ Energies= 0.5 0.5
diff --git a/Projects/S1554/Reaction/28Mgdp_1.reaction b/Projects/S1554/Reaction/28Mgdp_1.reaction
new file mode 100644
index 000000000..0d933064a
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_1.reaction
@@ -0,0 +1,40 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 2.2
+ CrossSectionPath= CS_28Mgdp_d.txt CSR1
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+GammaDecay 29Mg
+ Cascade
+ BranchingRatio= 70
+ Energies= 2.2
+ DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma25Na
+ Cascade
+ BranchingRatio= 30
+ Energies= 1.0 1.2
diff --git a/Projects/S1554/Reaction/28Mgdp_2.reaction b/Projects/S1554/Reaction/28Mgdp_2.reaction
new file mode 100644
index 000000000..b953d35c6
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_2.reaction
@@ -0,0 +1,40 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 1.1
+ CrossSectionPath= CS_28Mgdp_p.txt CSR2
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%GammaDecay 25Na
+ Cascade
+ BranchingRatio= 70
+ Energies= 2.2
+ DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma25Na
+ Cascade
+ BranchingRatio= 30
+ Energies= 1.0 1.2
diff --git a/Projects/S1554/Reaction/28Mgdp_3.reaction b/Projects/S1554/Reaction/28Mgdp_3.reaction
new file mode 100644
index 000000000..79de41691
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_3.reaction
@@ -0,0 +1,40 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 1.4
+ CrossSectionPath= CS_28Mgdp_f.txt CSR3
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%GammaDecay 25Na
+ Cascade
+ BranchingRatio= 70
+ Energies= 2.2
+ DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma25Na
+ Cascade
+ BranchingRatio= 30
+ Energies= 1.0 1.2
diff --git a/Projects/S1554/Reaction/28Mgdp_4.reaction b/Projects/S1554/Reaction/28Mgdp_4.reaction
new file mode 100644
index 000000000..507121a94
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_4.reaction
@@ -0,0 +1,40 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 2.3
+ CrossSectionPath= CS_28Mgdp_f.txt CSR4
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%GammaDecay 25Na
+ Cascade
+ BranchingRatio= 70
+ Energies= 2.2
+ DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma25Na
+ Cascade
+ BranchingRatio= 30
+ Energies= 1.0 1.2
diff --git a/Projects/S1554/Reaction/28Mgdp_5.reaction b/Projects/S1554/Reaction/28Mgdp_5.reaction
new file mode 100644
index 000000000..75ce15731
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_5.reaction
@@ -0,0 +1,40 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 4.4
+ CrossSectionPath= CS_28Mgdp_f.txt CSR5
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%GammaDecay 25Na
+ Cascade
+ BranchingRatio= 70
+ Energies= 2.2
+ DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma25Na
+ Cascade
+ BranchingRatio= 30
+ Energies= 1.0 1.2
diff --git a/Projects/S1554/Reaction/28Mgdp_EXD.reaction b/Projects/S1554/Reaction/28Mgdp_EXD.reaction
new file mode 100644
index 000000000..f47723f6f
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_EXD.reaction
@@ -0,0 +1,30 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0.1635909
+ MeanY= 0.910980
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyDistribution= EXD.txt EX
+ CrossSectionPath= CS_28Mgdp_g.txt CSR0
+ ShootLight= 1
+ ShootHeavy= 1
diff --git a/Projects/S1554/Reaction/28Mgdp_gs.reaction b/Projects/S1554/Reaction/28Mgdp_gs.reaction
new file mode 100644
index 000000000..0628c58da
--- /dev/null
+++ b/Projects/S1554/Reaction/28Mgdp_gs.reaction
@@ -0,0 +1,30 @@
+%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 28Mg
+ ExcitationEnergy= 0
+ Energy= 224
+ SigmaEnergy= 0.448
+ SigmaThetaX= 0.01
+ SigmaPhiY= 0.01
+ SigmaX= 0.5
+ SigmaY= 0.5
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+ %EnergyProfilePath=
+ %XThetaXProfilePath=
+ %YPhiYProfilePath=
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+TwoBodyReaction
+ Beam= 28Mg
+ Target= 2H
+ Light= 1H
+ Heavy= 29Mg
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 0.0
+ CrossSectionPath= CS_28Mgdp_d.txt CSR0
+ ShootLight= 1
+ ShootHeavy= 1
diff --git a/Projects/S1554/macro/CS.cxx b/Projects/S1554/macro/CS.cxx
new file mode 100644
index 000000000..d9ad57428
--- /dev/null
+++ b/Projects/S1554/macro/CS.cxx
@@ -0,0 +1,202 @@
+TGraph* thdae;
+TGraphErrors* gdd;
+double minX = 1e30;
+double minC = 1e30;
+void NumericalMinimization(const char * minName ,const char *algoName);
+double Chi2(const double* parameter);
+////////////////////////////////////////////////////////////////////////////////
+void CS(){
+ TFile* file = new TFile("macro/Kine.root","READ");
+ TH2F* h2 = (TH2F*)file->FindObjectAny("hEXT2");
+
+ file = new TFile("SharcEfficiency_28Mgdd.root","READ");
+ TH1D* eff = (TH1D*) file->FindObjectAny("SolidAngleCM");
+
+ file = new TFile("Mg28dd_dae.root","READ");
+ thdae = (TGraph*) file->FindObjectAny("Graph");
+
+ file = new TFile("Mg28dd_pap.root","READ");
+ TGraph* thpap = (TGraph*) file->FindObjectAny("Graph");
+ thpap->SetLineColor(kGreen-3);
+
+ file = new TFile("Mg28dd_loh.root","READ");
+ TGraph* thloh = (TGraph*) file->FindObjectAny("Graph");
+ thloh->SetLineColor(kMagenta-3);
+
+ NPL::Reaction dd("28Mg(d,d)28Mg@222");
+ NPL::Reaction pp("28Mg(p,p)28Mg@222");
+
+ TGraph* kdd = dd.GetKinematicLine3();
+ TGraph* kpp = pp.GetKinematicLine3();
+ TGraph* tdd = dd.GetELabVersusThetaCM();
+ TGraph* tpp = pp.GetELabVersusThetaCM();
+
+ const double* xdd= kdd->GetX();
+ const double* ydd = kdd->GetY();
+ kdd= new TGraph(kdd->GetN(),ydd,xdd);
+
+ const double* xpp= kpp->GetX();
+ const double* ypp = kpp->GetY();
+ kpp= new TGraph(kpp->GetN(),ypp,xpp);
+
+ const double* xdd1= tdd->GetX();
+ const double* ydd1 = tdd->GetY();
+ tdd= new TGraph(tdd->GetN()/2,xdd1,ydd1);
+
+ const double* xpp2= tpp->GetX();
+ const double* ypp2 = tpp->GetY();
+ tpp= new TGraph(tpp->GetN()/2,xpp2,ypp2);
+
+ TH1D* hdd = new TH1D("hdd","hdd",300,0,180);
+ TH1D* hpp = new TH1D("hpp","hpp",300,0,180);
+
+ gdd = new TGraphErrors();
+ TGraphErrors* gpp = new TGraphErrors();
+
+ TF1* f = new TF1("f","gaus(0)+gaus(3)+pol4(6)",0,90);
+ f->SetNpx(1000);
+ f->SetParameter(0,20);
+ f->SetParameter(2,1);
+
+ f->SetParameter(3,20);
+ f->SetParameter(5,1);
+
+ f->SetParameter(6,3);
+ f->SetParameter(7,-0.01);
+
+ f->SetParLimits(0,0,50);
+ f->SetParLimits(3,0,50);
+
+ f->SetParLimits(2,1,2);
+ f->SetParLimits(5,1,2);
+
+ TF1* g = new TF1("g","gaus(0)",0,90);
+
+ unsigned int step = 100;
+ double eps = 30./step;
+ cout << "EPSILON = " << eps <<" MeV" << endl;
+ double E=0;
+ TH1D* p;
+ double* r;
+ for(unsigned int i = 0 ; i < step ; i++){
+ E = i*eps;
+
+ double lTCMd = tdd->Eval(E-eps*0.5);
+ double uTCMd = tdd->Eval(E+eps*0.5);
+
+ p = h2->ProjectionX(Form("p%i",i),h2->GetYaxis()->FindBin(E-eps*0.5),h2->GetYaxis()->FindBin(E+eps*0.5));
+ if(p->GetEntries()>300){
+ double Td = kdd->Eval(E);
+ double Tp = kpp->Eval(E);
+ double TCMd = tdd->Eval(E);
+ double TCMp = tpp->Eval(E);
+
+ p->Rebin(4);
+
+ p->Draw();
+
+ f->SetParameter(1,Td);
+ f->SetParameter(4,Tp);
+
+ f->SetParLimits(1,Td,Td);
+ f->SetParLimits(4,Tp,Tp);
+ double pi,po;
+ f->GetParLimits(3,pi,po);
+ p->Fit(f,"RBFQ");
+ gPad->Update();
+
+ r = f->GetParameters();
+ g->SetParameter(0,r[0]);
+ g->SetParameter(1,r[1]);
+ g->SetParameter(2,r[2]);
+
+ double dd = g->Integral(r[1]-3*r[2],r[1]+3*r[2]);
+ double ddo = dd ;
+ double dde = sqrt(dd);
+ double mynorm = 1;
+ double EFF = eff->Integral(eff->FindBin(uTCMd),eff->FindBin(lTCMd))/abs(uTCMd-lTCMd);
+ dd = mynorm*dd / eff->GetBinContent(eff->FindBin(TCMd));
+ dde = mynorm*dde / eff->GetBinContent(eff->FindBin(TCMd));
+ double ddx = 0.5*abs(uTCMd-lTCMd);
+ if(dd<1e50 && ddo>10 && abs(Td-Tp)>1){
+ hdd->Fill(tdd->Eval(E),dd);
+ gdd->SetPoint(gdd->GetN(),TCMd,dd);
+ gdd->SetPointError(gdd->GetN()-1,ddx,dde);
+ }
+
+ g->SetParameter(0,r[3]);
+ g->SetParameter(1,r[4]);
+ g->SetParameter(2,r[5]);
+ double pp = g->Integral(0,90);
+ double ppe = sqrt(pp);
+
+ hpp->Fill(tpp->Eval(E),pp);
+ gpp->SetPoint(gpp->GetN(),TCMp,pp);
+ gpp->SetPointError(gpp->GetN()-1,0,ppe);
+// gPad->WaitPrimitive();
+ }
+ }
+
+ new TCanvas();
+ gdd->Draw("ap");
+ gdd->GetXaxis()->SetTitle("#theta_{CM} (deg)");
+ gdd->GetYaxis()->SetTitle("d#sigma/d#Omega (mb/sr)");
+ gdd->GetYaxis()->SetRangeUser(1,1e5);
+ thdae->Draw("c");
+ thpap->Draw("c");
+ thloh->Draw("c");
+
+ gPad->SetLogy();
+ NumericalMinimization("Minuit2","Fumili2");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Chi2(const double* parameter){
+ double Chi2 = 0 ;
+ for (int i = 0; i < gdd->GetN() ; i++) {
+ if(gdd->GetX()[i]<33 && gdd->GetX()[i]>15 && gdd->GetY()[i]>1){
+ double chi = ( thdae->Eval( gdd->GetX()[i] ) - parameter[0]*gdd->GetY()[i]) / (parameter[0]*gdd->GetEY()[i] );
+ Chi2 += chi*chi ;
+ }
+ }
+ if(Chi2<minC){
+ minC= Chi2;
+ minX= parameter[0];
+ }
+ return Chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NumericalMinimization(const char * minName ,const char *algoName){
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+
+ // set tolerance , etc...
+ min->SetMaxFunctionCalls(100000); // for Minuit/Minuit2
+ min->SetPrecision(1e-6);
+ min->SetValidError(true);
+ min->SetPrintLevel(1);
+
+ // create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&Chi2,1);
+
+ min->SetFunction(f);
+ double* parameter = new double[1];
+ parameter[0] = 1e-4;
+ // Set the free variables to be minimized!
+ min->SetLimitedVariable(0,"Normalistation",parameter[0],1e-5,1e-5,1);
+
+ // do the minimization
+ min->Minimize();
+
+ const double *xs = min->X();
+ std::cout << "Minimum = " << xs[0] << " but MinX = " << minX << std::endl;
+
+ for(int i = 0; i < gdd->GetN() ; i++ ){
+ gdd->SetPoint(i,gdd->GetX()[i],minX*gdd->GetY()[i]);
+ gdd->SetPointError(i,gdd->GetEX()[i],minX*gdd->GetEY()[i]);
+
+ }
+ return 0;
+}
diff --git a/Projects/S1554/macro/CheckGamma.cxx b/Projects/S1554/macro/CheckGamma.cxx
new file mode 100644
index 000000000..72591a18b
--- /dev/null
+++ b/Projects/S1554/macro/CheckGamma.cxx
@@ -0,0 +1,23 @@
+
+
+void CheckGamma(){
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Si_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ tree->Draw("AddBack_DC/1000>>h6(250,0,8)","AddBack_Theta < 125 && AddBack_Theta > 120 ","");
+ tree->Draw("AddBack_DC/1000>>h1(250,0,8)","AddBack_Theta < 80 && AddBack_Theta > 70 ","same");
+ tree->Draw("AddBack_DC/1000>>h2(250,0,8)","AddBack_Theta < 85 && AddBack_Theta > 80 ","same");
+ tree->Draw("AddBack_DC/1000>>h3(250,0,8)","AddBack_Theta < 100 && AddBack_Theta > 90 ","same");
+ tree->Draw("AddBack_DC/1000>>h4(250,0,8)","AddBack_Theta < 105.5 && AddBack_Theta > 104 ","same");
+ tree->Draw("AddBack_DC/1000>>h5(250,0,8)","AddBack_Theta < 107 && AddBack_Theta > 106 ","same");
+ tree->Draw("AddBack_DC/1000>>h7(250,0,8)","AddBack_Theta < 150 && AddBack_Theta > 140 ","same");
+
+ TH1* h ;
+ for(unsigned int i = 0 ; i < 7 ; i++){
+ h = (TH1*) gDirectory->FindObjectAny(Form("h%d",i+1));
+ h->SetFillStyle(0);
+ h->SetLineColor(kAzure+i+1);
+ }
+
+
+
diff --git a/Projects/S1554/macro/GammaRes.cxx b/Projects/S1554/macro/GammaRes.cxx
new file mode 100644
index 000000000..0f839fef3
--- /dev/null
+++ b/Projects/S1554/macro/GammaRes.cxx
@@ -0,0 +1,34 @@
+
+void GammaRes(){
+ TGraph* Res = new TGraph();
+ TLorentzVector GammaLV;
+ TVector3 position;
+ TVector3 beta(0,0,0.10) ;
+ for(unsigned int i = 0 ; i < 1000 ; i++){
+
+ double E = i*10;
+ position =TVector3(0,110+45,0);
+ position.SetMag(1);
+ GammaLV.SetPx(E*position.X());
+ GammaLV.SetPy(E*position.Y());
+ GammaLV.SetPz(E*position.Z());
+ GammaLV.SetE(E);
+ GammaLV.Boost(-beta);
+ double E1 = GammaLV.Energy();
+ position =TVector3(0,110+45,33.4*2);
+ position.SetMag(1);
+ GammaLV.SetPx(E*position.X());
+ GammaLV.SetPy(E*position.Y());
+ GammaLV.SetPz(E*position.Z());
+ GammaLV.SetE(E);
+ GammaLV.Boost(-beta);
+ double E2 = GammaLV.Energy();
+
+ Res->SetPoint(Res->GetN(),E,abs(E1-E2)*0.5);
+ }
+
+ Res->Draw("al");
+ cout << Res->Eval(1000)<<endl ;
+ cout << Res->Eval(5000)<<endl ;
+
+}
diff --git a/Projects/S1554/macro/Mg.cxx b/Projects/S1554/macro/Mg.cxx
new file mode 100644
index 000000000..6a3c1ec97
--- /dev/null
+++ b/Projects/S1554/macro/Mg.cxx
@@ -0,0 +1,122 @@
+void Mg(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Mgdp_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ file = new TFile("cuts/Cut_pp.root");
+ TCutG* Cut_pp = (TCutG*) file->FindObjectAny("Cut_pp");
+
+ tree->SetAlias("GoodTrifoil","(Trifoil.Time>150 && Trifoil.Time<160 && RunNumber>35110) && ELab > 0 ") ;
+ tree->SetAlias("GoodProton","ThetaLab>90 && Sharc.Strip_E>0.8");
+ TCanvas* c = new TCanvas("Mg","Mg");
+ c->Divide(3,3);
+
+ // Setting Trifoil entry list for fast drawing
+ TFile* fileEL = new TFile("macro/GoodTrifoil.root");
+ TEntryList* EL = (TEntryList*) fileEL->FindObjectAny("TrifoilEL");
+ if(!EL){
+ cout << "Trifoil Entry List does not exist, generating : ";
+ tree->Draw(">>TrifoilEL","GoodTrifoil","entrylist");
+ EL = (TEntryList*) gDirectory->FindObjectAny("TrifoilEL");
+ cout << " done" << endl;
+ EL->SaveAs("macro/GoodTrifoil.root");
+ }
+
+ cout <<"Setting Trifoil entry list :" ;
+ tree->SetEntryList(EL);
+ cout << " Done" << endl;
+
+
+ c->cd(1);
+ tree->Draw("Ex:AddBack_DC/1000.>>h2(2000,0,8,200,-1,8)","ThetaLab>90 && GoodProton","colz");
+ TH1* h2 = (TH1*) gDirectory->FindObjectAny("h2");
+ h2->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
+ h2->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
+ gPad->SetLogz();
+ TLine* line = new TLine(0,0,8,8);
+ line -> Draw("c");
+
+ c->cd(2);
+ tree->Draw("Ex>>h(200,-1,8)","GoodTrifoil && GoodProton ","colz");
+ TH1* h = (TH1*) gDirectory->FindObjectAny("h");
+ h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
+ h->GetYaxis()->SetTitle("counts/40 keV");
+
+ int maximum = h->GetMaximum();
+ TLine* state1= new TLine(0,0,0,maximum);
+ state1->Draw("l");
+
+ TLine* state2= new TLine(0.054,0,0.054,maximum);
+ state2->Draw("l");
+
+ TLine* state3= new TLine(1.095,0,1.095,maximum);
+ state3->Draw("l");
+
+ TLine* state4= new TLine(1.431,0,1.431,maximum);
+ state4->Draw("l");
+
+ TLine* state5= new TLine(2.266,0,2.266,maximum);
+ state5->Draw("l");
+
+ TLine* state6= new TLine(2.5,0,2.5,maximum);
+ state6->Draw("l");
+
+ TLine* state7= new TLine(3.2,0,3.2,maximum);
+ state7->Draw("l");
+
+ TLine* state8= new TLine(4.43,0,4.43,maximum);
+ state8->Draw("l");
+
+ c->cd(3);
+ tree->Draw("ELab:ThetaLab>>hK(180,90,180,900,0,14)"," GoodProton","colz");
+ NPL::Reaction r("28Mg(d,p)29Mg@222");
+ r.GetKinematicLine3()->Draw("c");
+
+ NPL::Reaction r2("28Mg(d,d)28Mg@222");
+ r2.GetKinematicLine3()->Draw("c");
+
+ NPL::Reaction r3("28Mg(p,p)28Mg@222");
+ r3.GetKinematicLine3()->Draw("c");
+
+
+ TH2* hk = (TH2*) gDirectory->FindObjectAny("hK");
+ hk->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ hk->GetYaxis()->SetTitle("E_{Lab} (MeV)");
+
+ c->cd(4);
+ tree->Draw("AddBack_DC/1000>>hDC(300,0,3)","ELab>0&& ThetaLab>90 && GoodProton && Ex > 0.5 && Ex<2","");
+ gPad->SetLogz();
+ h = (TH1*) gDirectory->FindObjectAny("hDC");
+ h->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
+ h->GetYaxis()->SetTitle("counts / 10 keV");
+
+ c->cd(5);
+ tree->Draw("Ex>>h11(50,-1,8)","GoodProton && AddBack_DC/1000. > 1.0 && AddBack_DC/1000. < 1.06 ","");
+ h = (TH1*) gDirectory->FindObjectAny("h11");
+ h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
+ tree->Draw("Ex>>h12(50,-1,8)","GoodProton && AddBack_DC/1000. > 0.3 && AddBack_DC/1000. < 0.36 ","same");
+ h = (TH1*) gDirectory->FindObjectAny("h12");
+ h->SetFillColor(kBlack);
+ h->SetLineColor(kBlack);
+ h->SetFillStyle(3004);
+
+ c->cd(6);
+ tree->Draw("AddBack_DC/1000>>hDC2(300,0,3)","ELab>0&& ThetaLab>90 && GoodProton && Ex > 2 && Ex<3","");
+ gPad->SetLogz();
+ h = (TH1*) gDirectory->FindObjectAny("hDC2");
+ h->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
+ h->GetYaxis()->SetTitle("counts / 10 keV");
+
+
+ c->cd(7);
+ tree->Draw("Ex>>h13(50,-1+0.1,8+0.1)","GoodProton && AddBack_DC/1000. > 2.2 && AddBack_DC/1000. < 2.28 ","");
+ h = (TH1*) gDirectory->FindObjectAny("h11");
+ h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
+
+ tree->Draw("Ex>>h14(50,-1+0.1,8+0.1)","GoodProton && AddBack_DC/1000. > 2.35 && AddBack_DC/1000. < 2.65 ","same");
+ h = (TH1*) gDirectory->FindObjectAny("h14");
+ h->SetFillColor(kBlack);
+ h->SetLineColor(kBlack);
+ h->SetFillStyle(3004);
+}
+
diff --git a/Projects/S1554/macro/MgCS.cxx b/Projects/S1554/macro/MgCS.cxx
new file mode 100644
index 000000000..860ad5d12
--- /dev/null
+++ b/Projects/S1554/macro/MgCS.cxx
@@ -0,0 +1,344 @@
+void Scale(TGraph* g , TGraphErrors* ex);
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
+
+////////////////////////////////////////////////////////////////////////////////
+void MgCS(){
+
+
+ file = new TFile("SharcEfficiency_29Mg.root");
+ TH2* Eff = (TH2*) file->FindObjectAny("SolidAngleCM_2D");
+
+ //TCanvas* c0 = new TCanvas("Mg0","Mg0");
+
+ Goodfile = new TFile("GoodExThetaCM.root");
+ TH2D* goodEx = (TH2D*) Goodfile->FindObjectAny("hexcmG");
+
+ Badfile = new TFile("BadExThetaCM.root");
+ TH2D* badEx = (TH2D*) Badfile->FindObjectAny("hexcmB");
+
+ TH2D* h = goodEx;
+ badEx->Scale(7./(145-1+300-160));
+ // h->Sumw2();
+ h->SetMarkerSize(0.5);
+ h->SetMarkerColor(kAzure+7);
+ // badEx->Sumw2();
+// h->Draw("colz");
+ h->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("#theta_{CM} (deg)");
+ h->GetXaxis()->SetRangeUser(-1,50);
+
+ vector<double> E,W;
+ E.push_back(0.0);
+ E.push_back(0.054);
+// E.push_back(0.590);
+
+ E.push_back(1.094);
+ E.push_back(1.431);
+// E.push_back(1.638);
+
+ E.push_back(2.266);
+ E.push_back(2.500);
+// E.push_back(2.8);
+ E.push_back(3.223);
+ E.push_back(3.674);
+ E.push_back(3.985);
+ E.push_back(4.280);
+ E.push_back(4.43);
+
+
+ string formula="[0]+[1]*x+";
+
+ for (unsigned int i = 0 ; i < E.size() ; i++){
+ //formula+=Form("gaus(%i)",i*3);
+ formula+=Form("[%i]*exp( -(x-[5]-[%i])*(x-[5]-[%i])/(([2]+[3]*x+[4]*x*x)*([2]+[3]*x+[4]*x*x)))",i*2+6,i*2+6+1,i*2+6+1);
+ if(i!=E.size()-1)
+ formula+="+";
+ W.push_back(1.1*(0.144206+0.00644644*E[i]+0.0006305*E[i]*E[i]));
+ }
+
+ TGraph* g_shift = new TGraph();
+ vector<TGraph*> g_width;
+ for (unsigned int i = 0 ; i < E.size() ; i++){
+ g_width.push_back(new TGraph);
+ }
+ TF1* bck = new TF1("bck","[0]+[1]*x",-6,-1);
+ bck->SetLineColor(kBlack);
+ TF1* func = new TF1("func",formula.c_str(),-2,5);
+ func->FixParameter(0,0);
+ func->FixParameter(1,0);
+
+ func->FixParameter(2,0.226217);
+ func->FixParameter(3,0.00294652);
+ func->FixParameter(4,0);
+ func->FixParameter(5,-0.055);
+
+ // func->SetParLimits(2,0.1,0.3);
+ // func->SetParLimits(3,0,1);
+ // func->SetParLimits(4,0,1e-31);
+ // func->SetParLimits(5,-0.06,0.06);
+
+ func->SetNpx(500);
+
+ double Energy;
+ for (unsigned int i = 0 ; i < E.size() ; i++){
+ func->SetParameter(i*2+6,10);
+ func->SetParLimits(i*2+6,0,100);
+ func->FixParameter(i*2+6+1,E[i]);
+ }
+
+ //////////
+ vector<TGraphErrors*> CS;
+ CS.resize(E.size(),NULL);
+
+ for(unsigned int i = 0 ; i < CS.size() ; i++){
+ CS[i]=new TGraphErrors();
+ }
+
+ const double* par;
+ unsigned int nbin = 10;
+ double minCM = 2;
+ double maxCM = 42;
+ double step = (maxCM-minCM)/nbin;
+ Eff->RebinX(step);
+
+
+ TH1D* p; TH1D* bbb;
+ double ThetaCM;
+ int pad = 1;
+
+ TCanvas* c4 = new TCanvas("Fit","Fit",1800,1000);
+ c4->Divide(4,2,1e-4,1e-4);
+ c4->cd(1);
+ for(unsigned int i = 0 ; i < nbin ; i++){
+ ThetaCM = minCM+i*step;
+ p = h->ProjectionY(Form("p%i",i),h->GetXaxis()->FindBin(ThetaCM),h->GetXaxis()->FindBin(ThetaCM+step));
+ p->Rebin(2);
+ // Background
+ bbb = badEx->ProjectionY(Form("b%i",i),badEx->GetXaxis()->FindBin(ThetaCM),badEx->GetXaxis()->FindBin(ThetaCM+step));
+ bbb->Rebin(2);
+ // Substract the Background
+ p->Add(bbb,-1);
+ // Presentation
+ p->SetMarkerSize(0.5);
+ p->SetMarkerColor(kAzure+7);
+ bbb->SetMarkerSize(0.1);
+ bbb->SetLineColor(kOrange+7);
+ bbb->SetFillColor(kOrange+7);
+ bbb->SetFillStyle(3004);
+
+ p->GetYaxis()->SetRangeUser(0,30);
+ p->GetXaxis()->SetRangeUser(-1,6);
+
+ // Limit the heigh of the gaussian to the bin at the state position
+ for(unsigned int n = 0 ; n < CS.size() ; n++){
+ func->SetParameter(n*2+6,p->GetBinContent(p->FindBin(func->GetParameter(n*2+6+1))));
+ double integral = p->Integral(p->FindBin(E[n]-2*W[n]),p->FindBin(E[n]+2*W[n]));
+ if(integral >0)
+ func->SetParLimits(n*2+6,0, integral+1e-20);
+ else
+ func->SetParLimits(n*2+6,0, 10);
+ }
+ c4->cd(pad++);
+ p->Draw();
+ // p->Fit(bck,"RBFQ");
+ bbb->Draw("same");
+
+ // func->FixParameter(0,bck->GetParameter(0));
+ // func->FixParameter(1,bck->GetParameter(1));
+ p->Fit(func,"RBFQ+");
+ par = func->GetParameters();
+ TLatex* tex = new TLatex(2,25,Form("%.1f deg",ThetaCM+0.5*step));
+ tex->Draw();
+ // bck->Draw("same");
+
+ g_shift->SetPoint(g_shift->GetN(),ThetaCM,func->GetParameter(5));
+ if(p->Integral(p->FindBin(-1),p->FindBin(6))>-1000 ){
+ for(unsigned int n = 0 ; n < CS.size() ; n++){
+ double SolidAngle = Eff->ProjectionX("px",Eff->GetYaxis()->FindBin(E[n]-W[n]),Eff->GetYaxis()->FindBin(E[n]+W[n]))->Interpolate(ThetaCM);
+ double Width = func->GetParameter(2)+func->GetParameter(3)*E[n]+func->GetParameter(4)*E[n]*E[n];
+ g_width[n]->SetPoint(g_width[n]->GetN(),ThetaCM,Width);
+ CS[n]->SetPoint(CS[n]->GetN(),ThetaCM+0.5*step,par[n*2+6]/SolidAngle);
+ CS[n]->SetPointError(CS[n]->GetN()-1,0,func->GetParError(n*2+6)/SolidAngle);
+ }
+ }
+ else{
+ pad= pad-1;
+ }
+ }
+
+ // Draw Width of each state vs ThetaCM
+ c4->cd(pad++);
+ const UInt_t Number = 2;
+ Double_t Red[Number] = { 0, 0 };
+ Double_t Green[Number] = { 0, 0.8 };
+ Double_t Blue[Number] = { 0, 1.00 };
+
+ Double_t Length[Number] = { 0,1.00 };
+ Int_t nb=g_width.size() ;
+ Int_t color = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
+
+ for(unsigned int n = 0 ; n < g_width.size() ; n++){
+ g_width[n]->SetMarkerColor(color++) ;
+ g_width[n]->SetMarkerSize(0.4);
+ if(n==0){
+ g_width[n]->Draw("ap");
+ g_width[n]->GetYaxis()->SetRangeUser(0.3,1);
+ }
+ else
+ g_width[n]->Draw("p");
+ }
+
+ // Draw the shift of Ex spectrum for each ThetaCM
+ c4->cd(pad++);
+ g_shift->SetLineColor(kAzure+7);
+ g_shift->Draw("acp");
+
+
+
+ TCanvas* c = new TCanvas("Mg","Mg",2000,2000*2./5.);
+ c->Divide(5,2,0.00001,0.00001);
+
+
+ for(unsigned int i = 0 ; i < CS.size() ; i++){
+ c->cd(i+1);
+
+ CS[i]->Draw("ap");
+ gPad->SetLogy();
+
+ TGraph* gs = TWOFNR(E[i], 0,0.5,1,0, 0.5);
+ TGraph* gp = TWOFNR(E[i], 0,1.5,1,1, 1.5);
+ TGraph* gd = TWOFNR(E[i], 0,1.5,0,2, 1.5);
+ TGraph* gf = TWOFNR(E[i], 0,3.5,0,3, 3.5);
+ TGraph* gg = TWOFNR(E[i], 0,4.5,0,4, 4.5);
+
+ Scale(gs,CS[i]);
+ Scale(gp,CS[i]);
+ Scale(gd,CS[i]);
+ Scale(gf,CS[i]);
+ Scale(gg,CS[i]);
+
+ gs->SetLineColor(kGreen);
+ gp->SetLineColor(kBlue);
+ gd->SetLineColor(kRed);
+ gf->SetLineColor(kBlack);
+ gg->SetLineColor(kAzure+7);
+
+ if(i!=0){
+ gp->Draw("c");
+ gd->Draw("c");
+ gf->Draw("c");
+ gs->Draw("c");
+ gg->Draw("c");
+ TLegend* leg = new TLegend(0.7,0.6,1,0.9);
+ leg->SetBorderSize(0);
+ leg->SetFillStyle(0);
+ leg->AddEntry(gs,"L=0","l");
+ leg->AddEntry(gp,"L=1","l");
+ leg->AddEntry(gd,"L=2","l");
+ leg->AddEntry(gf,"L=3","l");
+ leg->AddEntry(gg,"L=4","l");
+ leg->Draw();
+ }
+
+ else{
+ gs->Draw("c");
+ gd->Draw("c");
+ TGraph* gm = new TGraph(*gs);
+ double* XX = gm->GetX();
+ double* YY = gm->GetY();
+ double* ZZ = gd->GetY();
+ for(int n = 0 ; n < gs->GetN() ; n++)
+ gm->SetPoint(n,XX[n],YY[n]+ZZ[n]);
+ Scale(gm,CS[i]);
+ gm->SetLineStyle(2);
+ gm->SetLineColor(kGreen-3);
+ gm->Draw("c");
+ TLegend* leg = new TLegend(0.7,0.6,1,0.9);
+ leg->SetBorderSize(0);
+ leg->SetFillStyle(0);
+ leg->AddEntry(gs,"L=0","l");
+ leg->AddEntry(gd,"L=2","l");
+ leg->AddEntry(gm,"L=2+0","l");
+ leg->Draw();
+ }
+
+ TLatex* tex = new TLatex(10,1.7e3,Form("%.2f MeV",E[i]));
+ tex->Draw();
+ CS[i]->GetYaxis()->SetRangeUser(1,1e4);
+ CS[i]->GetXaxis()->SetRangeUser(0,40);
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+void Scale(TGraph* g , TGraphErrors* ex){
+ double scale;
+ double mean = 0 ;
+ double* eX = ex->GetX();
+ double* eY = ex->GetY();
+ double totalW = 0;
+ double W = 0;
+ for(int i = 0 ; i < ex->GetN() ; i++){
+ if(eY[i]>1 && eY[i] <1e4){
+ // Incremental Error weighted average
+ W = 1./ex->GetErrorY(i);
+ scale = eY[i]/g->Eval(eX[i]);
+ totalW +=W;
+ mean = mean + W*(scale - mean)/(totalW);
+ }
+ }
+
+ double* x = g->GetX();
+ double* y = g->GetY();
+ for(unsigned int i = 0 ; i < g->GetN() ; i++)
+ g->SetPoint(i,x[i],y[i]*mean);
+}
+////////////////////////////////////////////////////////////////////////////////
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
+ double BeamEnergy = 8;
+
+ NPL::Reaction r("28Mg(d,p)29Mg@224");
+ r.SetExcitationHeavy(E);
+ double QValue = r.GetQValue();
+
+ std::ofstream Front_Input("in.front");
+ Front_Input << "jjj" << std::endl;
+ Front_Input << "pipo" << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << BeamEnergy << std::endl;
+ Front_Input << 28 << " " << 12 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << "0 0 0" << std::endl;
+ Front_Input << l << " " << j << std::endl;
+ Front_Input << n << std::endl;
+ Front_Input << 2 << std::endl;
+
+ // unbound case:
+ if( QValue < 0 )
+ Front_Input << 0.01 << std::endl;
+ else
+ Front_Input << QValue << std::endl;
+
+ Front_Input << 1 << std::endl;
+ Front_Input << J0 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 5 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << J << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1.25 << " " << 0.65 << std::endl;
+ Front_Input << 6.0 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << 0 << std::endl;
+
+ Front_Input.close() ;
+
+ system("(exec FRONT < in.front &> /dev/null)");
+ system("(exec echo tran.jjj | TWOFNR &> /dev/null)");
+
+ TGraph* CS = new TGraph("22.jjj");
+ return CS;
+}
diff --git a/Projects/S1554/macro/MgCS2.cxx b/Projects/S1554/macro/MgCS2.cxx
new file mode 100644
index 000000000..13b54c2fe
--- /dev/null
+++ b/Projects/S1554/macro/MgCS2.cxx
@@ -0,0 +1,238 @@
+void Scale(TGraph* g , TGraphErrors* ex);
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
+TH1* Region(double Emin, double Emax,TH2* h, TCanvas* c, int pad, TH2* Eff);
+double ToMininize(const double* parameter);
+TGraph* FindNormalisation(TGraph* cs1, TGraph* cs2, TH1* hexp);
+
+TGraph* g1;
+TGraph* g2;
+TH1* current;
+
+////////////////////////////////////////////////////////////////////////////////
+void MgCS2(){
+ file = new TFile("SharcEfficiency_29Mg.root");
+ TH2* Eff = (TH2*) file->FindObjectAny("SolidAngleCM_2D");
+
+ Goodfile = new TFile("GoodExThetaCM.root");
+ TH2D* goodEx = (TH2D*) Goodfile->FindObjectAny("hexcmG");
+
+ Badfile = new TFile("BadExThetaCM.root");
+ TH2D* badEx = (TH2D*) Badfile->FindObjectAny("hexcmB");
+
+ TH2D* h = goodEx;
+ badEx->Scale(7./(145-1+300-160));
+ h->Sumw2();
+ h->SetMarkerSize(0.5);
+ h->SetMarkerColor(kAzure+7);
+ badEx->Sumw2();
+ h->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("#theta_{CM} (deg)");
+ h->Add(badEx,-1);
+ Eff->Sumw2();
+ h->Scale(0.010048*10/7.37);
+ h->Divide(Eff);
+
+ TH1D* p;
+ int nn=0;
+ // Full Ex
+ p = h->ProjectionY(Form("p%i",nn++),h->GetXaxis()->FindBin(0.),h->GetXaxis()->FindBin(50.));
+ p->Draw();
+
+ new TCanvas();
+ h->Draw("colz");
+ // CS
+ TCanvas* c = new TCanvas("CS","CS",1000,1000);
+ c->Divide(2,2);
+
+ // Region
+ current = Region(-0.1,0.4,h,c,nn++,Eff);
+ g1 = TWOFNR(0.000, 0, 0.5, 1, 0, 0.5);
+ g2 = TWOFNR(0.054, 0, 1.5, 0, 2, 1.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+ current = Region(0.6,1.8,h,c,nn++,Eff);
+ g1 = TWOFNR(1.095, 0, 1.5, 1, 1, 1.5);
+ g2 = TWOFNR(1.431, 0, 3.5, 0, 3, 3.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+ current = Region(2,3,h,c,nn++,Eff);
+ g1 = TWOFNR(2.266, 0, 0.5, 1, 1, 0.5);
+ g2 = TWOFNR(2.500, 0, 1.5, 0, 2, 1.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+ current = Region(3,4.8,h,c,nn++,Eff);
+ g1 = TWOFNR(3.2, 0, 2.5, 0, 2, 2.5);
+ g2 = TWOFNR(3.2, 0, 3.5, 0, 3, 3.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+
+}
+////////////////////////////////////////////////////////////////////////////////
+TH1* Region(double Emin, double Emax,TH2* h, TCanvas* c, int pad , TH2* Eff){
+ c->cd(pad);
+ TH1D* p = h->ProjectionX(Form("p%i",pad),h->GetYaxis()->FindBin(Emin),h->GetYaxis()->FindBin(Emax));
+ p->Draw();
+ gPad->SetLogy();
+// p->GetYaxis()->SetRangeUser(1e1,1e5);
+ return p;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Chi2(TGraph* g , TH1* h){
+ double Chi2 = 0 ;
+ double chi;
+
+ for(int i = 1 ; i < h->GetNbinsX() ; i++){
+ if(h->GetBinContent(i)>0) {
+ chi = (h->GetBinContent(i) - g->Eval(h->GetBinCenter(i)) ) / (h->GetBinError(i));
+ Chi2 +=chi*chi;
+ }
+ }
+
+ return Chi2;
+}
+////////////////////////////////////////////////////////////////////////////////
+double ToMininize(const double* parameter){
+static int f = 0 ;
+ TGraph* g = new TGraph();
+ double* X = g1->GetX();
+ double* Y = g1->GetY();
+ for(int i = 0 ; i < g1->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i] + parameter[1]*g2->Eval(X[i]) );
+ }
+ double chi2 = Chi2(g,current);
+ return chi2;
+}
+////////////////////////////////////////////////////////////////////////////////
+void Scale(TGraph* g , TGraphErrors* ex){
+ double scale;
+ double mean = 0 ;
+ double* eX = ex->GetX();
+ double* eY = ex->GetY();
+ double totalW = 0;
+ double W = 0;
+ for(int i = 0 ; i < ex->GetN() ; i++){
+ if(eY[i]>1 && eY[i] <1e4){
+ // Incremental Error weighted average
+ W = 1./ex->GetErrorY(i);
+ scale = eY[i]/g->Eval(eX[i]);
+ totalW +=W;
+ mean = mean + W*(scale - mean)/(totalW);
+ }
+ }
+
+ double* x = g->GetX();
+ double* y = g->GetY();
+ for(unsigned int i = 0 ; i < g->GetN() ; i++)
+ g->SetPoint(i,x[i],y[i]*mean);
+}
+////////////////////////////////////////////////////////////////////////////////
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
+ double BeamEnergy = 8;
+
+ NPL::Reaction r("28Mg(d,p)29Mg@224");
+ r.SetExcitationHeavy(E);
+ double QValue = r.GetQValue();
+
+ std::ofstream Front_Input("in.front");
+ Front_Input << "jjj" << std::endl;
+ Front_Input << "pipo" << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << BeamEnergy << std::endl;
+ Front_Input << 28 << " " << 12 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << "0 0 0" << std::endl;
+ Front_Input << l << " " << j << std::endl;
+ Front_Input << n << std::endl;
+ Front_Input << 2 << std::endl;
+
+ // unbound case:
+// if( QValue < 0 )
+// Front_Input << 0.01 << std::endl;
+// else
+ Front_Input << QValue << std::endl;
+
+ Front_Input << 1 << std::endl;
+ Front_Input << J0 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 5 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << J << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1.25 << " " << 0.65 << std::endl;
+ Front_Input << 6.0 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << 0 << std::endl;
+
+ Front_Input.close() ;
+
+ system("(exec FRONT < in.front &> /dev/null)");
+ system("(exec echo tran.jjj | TWOFNR &> /dev/null)");
+ // system("exec FRONT < in.front");
+ // system("(exec echo tran.jjj | TWOFNR)");
+
+
+ TGraph* CS = new TGraph("22.jjj");
+ return CS;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* FindNormalisation(TGraph* cs1, TGraph* cs2, TH1* hexp){
+ // Set global variable
+ g1 = cs1;
+ g2 = cs2;
+ current = hexp;
+
+ const char* minName ="Minuit";const char* algoName="Fumili2";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+ // Number of parameter
+ mysize = 2;
+
+ // create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ parameter[i] = 0.5;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ min->SetLimitedVariable(i,name,parameter[i],0.1,0,1000);
+ }
+
+ // do the minimization
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] <<"("<<err[i] << ") ";
+ }
+ cout << endl;
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = cs1->GetX();
+ double* Y = cs1->GetY();
+ for(int i = 0 ; i < cs1->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],xs[0]*Y[i] + xs[1]*cs2->Eval(X[i]) );
+ }
+ return g;
+}
+
diff --git a/Projects/S1554/macro/MgCSLab.cxx b/Projects/S1554/macro/MgCSLab.cxx
new file mode 100644
index 000000000..dad050953
--- /dev/null
+++ b/Projects/S1554/macro/MgCSLab.cxx
@@ -0,0 +1,238 @@
+void Scale(TGraph* g , TGraphErrors* ex);
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
+TH1* Region(double Emin, double Emax,TH2* h, TCanvas* c, int pad, TH2* Eff);
+double ToMininize(const double* parameter);
+TGraph* FindNormalisation(TGraph* cs1, TGraph* cs2, TH1* hexp);
+
+TGraph* g1;
+TGraph* g2;
+TH1* current;
+
+////////////////////////////////////////////////////////////////////////////////
+void MgCSLab(){
+ file = new TFile("SharcEfficiency_29Mg.root");
+ TH2* Eff = (TH2*) file->FindObjectAny("SolidAngleLab_2D");
+
+ Goodfile = new TFile("GoodExThetaLab.root");
+ TH2D* goodEx = (TH2D*) Goodfile->FindObjectAny("hexcmG");
+
+ Badfile = new TFile("BadExThetaLab.root");
+ TH2D* badEx = (TH2D*) Badfile->FindObjectAny("hexcmB");
+
+ TH2D* h = goodEx;
+ badEx->Scale(7./(145-1+300-160));
+ h->Sumw2();
+ h->SetMarkerSize(0.5);
+ h->SetMarkerColor(kAzure+7);
+ badEx->Sumw2();
+ h->GetYaxis()->SetTitle("E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("#theta_{CM} (deg)");
+ h->Add(badEx,-1);
+ Eff->Sumw2();
+ h->Scale(0.010048*10/7.37);
+ h->Divide(Eff);
+
+ TH1D* p;
+ int nn=0;
+ // Full Ex
+ p = h->ProjectionY(Form("p%i",nn++),h->GetXaxis()->FindBin(0.),h->GetXaxis()->FindBin(50.));
+ p->Draw();
+
+ new TCanvas();
+ h->Draw("colz");
+ // CS
+ TCanvas* c = new TCanvas("CS","CS",1000,1000);
+ c->Divide(2,2);
+
+ // Region
+ current = Region(-0.1,0.4,h,c,nn++,Eff);
+ g1 = TWOFNR(0.000, 0, 0.5, 1, 0, 0.5);
+ g2 = TWOFNR(0.054, 0, 1.5, 0, 2, 1.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+ current = Region(0.6,1.8,h,c,nn++,Eff);
+ g1 = TWOFNR(1.095, 0, 1.5, 1, 1, 1.5);
+ g2 = TWOFNR(1.431, 0, 3.5, 0, 3, 3.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+ current = Region(2,3,h,c,nn++,Eff);
+ g1 = TWOFNR(2.266, 0, 0.5, 1, 1, 0.5);
+ g2 = TWOFNR(2.500, 0, 1.5, 0, 2, 1.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+ current = Region(3,4.8,h,c,nn++,Eff);
+ g1 = TWOFNR(3.2, 0, 2.5, 0, 2, 2.5);
+ g2 = TWOFNR(3.2, 0, 3.5, 0, 3, 3.5);
+ g1->SetLineStyle(2);g1->Draw("c");
+ g2->SetLineStyle(3);g2->Draw("c");
+ FindNormalisation(g1,g2,current)->Draw("c");
+
+
+}
+////////////////////////////////////////////////////////////////////////////////
+TH1* Region(double Emin, double Emax,TH2* h, TCanvas* c, int pad , TH2* Eff){
+ c->cd(pad);
+ TH1D* p = h->ProjectionX(Form("p%i",pad),h->GetYaxis()->FindBin(Emin),h->GetYaxis()->FindBin(Emax));
+ p->Draw();
+ gPad->SetLogy();
+// p->GetYaxis()->SetRangeUser(1e1,1e5);
+ return p;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Chi2(TGraph* g , TH1* h){
+ double Chi2 = 0 ;
+ double chi;
+
+ for(int i = 1 ; i < h->GetNbinsX() ; i++){
+ if(h->GetBinContent(i)>0) {
+ chi = (h->GetBinContent(i) - g->Eval(h->GetBinCenter(i)) ) / (h->GetBinError(i));
+ Chi2 +=chi*chi;
+ }
+ }
+
+ return Chi2;
+}
+////////////////////////////////////////////////////////////////////////////////
+double ToMininize(const double* parameter){
+static int f = 0 ;
+ TGraph* g = new TGraph();
+ double* X = g1->GetX();
+ double* Y = g1->GetY();
+ for(int i = 0 ; i < g1->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i] + parameter[1]*g2->Eval(X[i]) );
+ }
+ double chi2 = Chi2(g,current);
+ return chi2;
+}
+////////////////////////////////////////////////////////////////////////////////
+void Scale(TGraph* g , TGraphErrors* ex){
+ double scale;
+ double mean = 0 ;
+ double* eX = ex->GetX();
+ double* eY = ex->GetY();
+ double totalW = 0;
+ double W = 0;
+ for(int i = 0 ; i < ex->GetN() ; i++){
+ if(eY[i]>1 && eY[i] <1e4){
+ // Incremental Error weighted average
+ W = 1./ex->GetErrorY(i);
+ scale = eY[i]/g->Eval(eX[i]);
+ totalW +=W;
+ mean = mean + W*(scale - mean)/(totalW);
+ }
+ }
+
+ double* x = g->GetX();
+ double* y = g->GetY();
+ for(unsigned int i = 0 ; i < g->GetN() ; i++)
+ g->SetPoint(i,x[i],y[i]*mean);
+}
+////////////////////////////////////////////////////////////////////////////////
+TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
+ double BeamEnergy = 8;
+
+ NPL::Reaction r("28Mg(d,p)29Mg@224");
+ r.SetExcitationHeavy(E);
+ double QValue = r.GetQValue();
+
+ std::ofstream Front_Input("in.front");
+ Front_Input << "jjj" << std::endl;
+ Front_Input << "pipo" << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << BeamEnergy << std::endl;
+ Front_Input << 28 << " " << 12 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << "0 0 0" << std::endl;
+ Front_Input << l << " " << j << std::endl;
+ Front_Input << n << std::endl;
+ Front_Input << 2 << std::endl;
+
+ // unbound case:
+// if( QValue < 0 )
+// Front_Input << 0.01 << std::endl;
+// else
+ Front_Input << QValue << std::endl;
+
+ Front_Input << 1 << std::endl;
+ Front_Input << J0 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 5 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << J << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 2 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1.25 << " " << 0.65 << std::endl;
+ Front_Input << 6.0 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << 0 << std::endl;
+
+ Front_Input.close() ;
+
+ system("(exec FRONT < in.front &> /dev/null)");
+ system("(exec echo tran.jjj | TWOFNR &> /dev/null)");
+ // system("exec FRONT < in.front");
+ // system("(exec echo tran.jjj | TWOFNR)");
+
+
+ TGraph* CS = new TGraph("24.jjj");
+ return CS;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* FindNormalisation(TGraph* cs1, TGraph* cs2, TH1* hexp){
+ // Set global variable
+ g1 = cs1;
+ g2 = cs2;
+ current = hexp;
+
+ const char* minName ="Minuit";const char* algoName="Fumili2";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+ // Number of parameter
+ mysize = 2;
+
+ // create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ parameter[i] = 0.5;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ min->SetLimitedVariable(i,name,parameter[i],0.1,0,1000);
+ }
+
+ // do the minimization
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] <<"("<<err[i] << ") ";
+ }
+ cout << endl;
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = cs1->GetX();
+ double* Y = cs1->GetY();
+ for(int i = 0 ; i < cs1->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],xs[0]*Y[i] + xs[1]*cs2->Eval(X[i]) );
+ }
+ return g;
+}
+
diff --git a/Projects/S1554/macro/MgCheck.cxx b/Projects/S1554/macro/MgCheck.cxx
new file mode 100644
index 000000000..d1da673de
--- /dev/null
+++ b/Projects/S1554/macro/MgCheck.cxx
@@ -0,0 +1,35 @@
+void MgCheck(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Mg_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ tree->SetAlias("GoodTrifoil","((Trifoil.Time>150 && Trifoil.Time<160 && RunNumber>35110)||(Trifoil.Time>50 && Trifoil.Time<60 && RunNumber<35110)) &&ELab > 0 && Sharc.Strip_E>0.4") ;
+ TCanvas* c = new TCanvas("Mg","Mg",1500,1000);
+ tree->Draw("Ex>>h(200,-10,10)","GoodTrifoil && ThetaLab>90 ","");
+ TH1* h = (TH1*) gDirectory->FindObjectAny("h");
+ h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
+ int maximum = h->GetMaximum();
+ TLine* state1= new TLine(0,0,0,maximum);
+ state1->Draw("l");
+
+ TLine* state2= new TLine(0.054,0,0.054,maximum);
+ state2->Draw("l");
+
+ TLine* state3= new TLine(1.095,0,1.095,maximum);
+ state3->Draw("l");
+
+ TLine* state4= new TLine(1.431,0,1.431,maximum);
+ state4->Draw("l");
+
+ TLine* state5= new TLine(2.266,0,2.266,maximum);
+ state5->Draw("l");
+
+ TLine* state6= new TLine(2.5,0,2.5,maximum);
+ state6->Draw("l");
+
+ TLine* state7= new TLine(3.2,0,3.2,maximum);
+ state7->Draw("l");
+
+ TLine* state8= new TLine(4.43,0,4.43,maximum);
+ state8->Draw("l");
+}
diff --git a/Projects/S1554/macro/Mgdd.cxx b/Projects/S1554/macro/Mgdd.cxx
new file mode 100644
index 000000000..2cf94c352
--- /dev/null
+++ b/Projects/S1554/macro/Mgdd.cxx
@@ -0,0 +1,59 @@
+void Mgdd(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Mgdd_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+ tree->SetAlias("GoodTrifoil","Trifoil.Time>150 && Trifoil.Time<160 && ELab >0") ;
+ NPL::Reaction r("28Mg(d,d)28Mg@221.6");
+ NPL::Reaction s("28Si(d,d)28Si@221.6");
+
+ // Load Cuts
+ fileCut = new TFile("cuts/Cut_Kine_dd.root");
+ TCutG* Cut_Kine_dd = (TCutG*) fileCut->FindObjectAny("Cut_Kine_dd");
+
+ // Create Canvas
+ TCanvas* c = new TCanvas("Result","Result",1200,1000);
+ c->Divide(1,2);
+ TH1* h;
+ TH1F* h2;
+
+ // Setting Trifoil entry list for fast drawing
+ TFile* fileEL = new TFile("macro/GoodTrifoil.root");
+ TEntryList* EL = (TEntryList*) fileEL->FindObjectAny("TrifoilEL");
+ if(!EL){
+ cout << "Trifoil Entry List does not exist, generating : ";
+ tree->Draw(">>TrifoilEL","GoodTrifoil","entrylist");
+ EL = (TEntryList*) gDirectory->FindObjectAny("TrifoilEL");
+ cout << " done" << endl;
+ EL->SaveAs("macro/GoodTrifoil.root");
+ }
+ cout <<"Setting Trifoil entry list :" ;
+ tree->SetEntryList(EL);
+ cout << " Done" << endl;
+
+ c->cd(1);
+ tree->Draw("ELab:ThetaLab>>hEXT2(1000,0,180,1200,0,60)","ELab>0","colz");
+
+ r.GetKinematicLine3()->Draw("c");
+ NPL::Reaction p("28Mg(p,p)28Mg@221.6");
+ p.GetKinematicLine3()->Draw("c");
+ Cut_Kine_dd->Draw("same");
+
+ gPad->SetLogz();
+
+ h2 = (TH1F*) gDirectory->FindObjectAny("hEXT2");
+ h2->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ h2->GetYaxis()->SetTitle("E_{Lab} (MeV)");
+
+ TVirtualPad* pad = c->cd(2);
+ pad->Divide(3);
+
+ pad->cd(2);
+ tree->Draw("Y_Trifoil:X_Trifoil>>htr(200,-50,50,200,-50,50)","ELab>0 && Cut_Kine_dd","colz") ;
+ pad->cd(1);
+ tree->Draw("ELab:ThetaLab>>hEXTA(360,0,180,1200,0,60)","ELab>0 && Sharc.DetectorNumber==12","colz") ;
+ r.SetExcitationHeavy(0);
+ r.GetKinematicLine3()->Draw("c");
+ pad->cd(3);
+ tree->Draw("ELab:ThetaLab>>hEXTB(360,0,180,1200,0,60)","ELab>0 && Sharc.DetectorNumber==10","colz") ;
+ r.GetKinematicLine3()->Draw("c");
+ }
diff --git a/Projects/S1554/macro/Resolution.cxx b/Projects/S1554/macro/Resolution.cxx
new file mode 100644
index 000000000..2254035a7
--- /dev/null
+++ b/Projects/S1554/macro/Resolution.cxx
@@ -0,0 +1,61 @@
+void Resolution(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/PhysicsTree.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ TCanvas* c = new TCanvas("ExTheta","ExTheta");
+ c->Divide(3,3);
+ TH2* h ;
+ vector<TGraph*> g;
+ TGraph* gt = new TGraph();
+ double* pp;
+
+
+ const double* par;
+ unsigned int nbin = 10;
+ double minCM = 4;
+ double maxCM = 34;
+ double step = (maxCM-minCM)/nbin;
+ TH1D* p;
+ double ThetaCM;
+ g.resize(nbin,NULL);
+ int point = 0 ;
+
+ for (unsigned int i = 0 ; i < 7 ; i++){
+ c->cd(i+1);
+
+ int entries = tree->Draw(Form("Ex:ThetaCM>>h%i(180,0,180,1000,-1,10)",i),Form("Run==%i && Plastic.Energy@.size()==1 && ThetaLab>90&& Sharc.Strip_E > 0.8",i+1));
+ h = (TH2*) gDirectory->FindObjectAny(Form("h%i",i));
+
+ for(unsigned int b = 0 ; b < nbin ; b++){
+ if(g[b]==NULL)
+ g[b]=new TGraph();
+
+ ThetaCM = minCM+b*step;
+ p = h->ProjectionY(Form("p%i",b*10+i),h->GetXaxis()->FindBin(ThetaCM),h->GetXaxis()->FindBin(ThetaCM+step));
+ if(p->Integral()>0){
+ p->Fit("gaus","Q");
+ pp = p->GetFunction("gaus")->GetParameters();
+ if(pp[2]<0.18)
+ g[b]->SetPoint(g[b]->GetN(),pp[1],pp[2]);
+ gt->SetPoint(gt->GetN(),pp[1],pp[2]);
+
+ }
+ }
+ }
+
+
+ TCanvas* cg = new TCanvas("Res","Res");
+ cg->Divide(4,4);
+ for(unsigned int i = 0 ; i < nbin; i++){
+ cg->cd(i+1);
+ g[i]->Draw("ap");
+ g[i]->Fit("pol2","F");
+ }
+ cg->cd(nbin+1);
+ gt->Draw("ap");
+ gt->Fit("pol2","F");
+ cg->cd(nbin+2);
+ tree->Draw(Form("Ex:ThetaCM>>h%i(180,0,180,1000,-1,10)",1000),"Plastic.Energy@.size()==1 && ThetaLab>90 && Sharc.Strip_E > 0.8","colz");
+
+}
diff --git a/Projects/S1554/macro/SM.cxx b/Projects/S1554/macro/SM.cxx
new file mode 100644
index 000000000..935467787
--- /dev/null
+++ b/Projects/S1554/macro/SM.cxx
@@ -0,0 +1,143 @@
+#include"RSShellModelCollection.h"
+#include"RSShellModelState.h"
+#include"RSExperimentalCrossSection.h"
+void PlotLevel(string name,int offset,bool LR, RS::ShellModelCollection Collection, vector<int> matched);
+vector<TPad*> thePad;
+double Emin = -0.2;
+double Emax = 5.2;
+double ScaleMin = 0;
+double ScaleMax = 7;
+double LabelSize = 0.025;
+int nbr = 2;
+unsigned int pad =0 ;
+
+TCanvas* c = new TCanvas() ;
+TCanvas* c2 ;
+TPad* tPAD;
+int currentPad =0 ;
+
+vector<RS::ShellModelState> SingleFit;
+vector<RS::ShellModelState> DoubleFit;
+vector<double> dErr;
+vector<double> sErr;
+vector<int> dorder;
+vector<double> dSMs;
+vector<double> sSMs;
+vector<double> eSM;
+vector<string> ExpFileName;
+
+
+void SM(){
+ ifstream infile("mg29olapc.lsf");
+ RS::ShellModelCollection WBP("WBP");
+ WBP.LoadCollectionFromNushell("mg29olapc.lsf");
+ RS::ShellModelCollection Exp("Exp");
+ Exp.LoadCollectionFromSimpleFile("Level.txt");
+
+ vector<int> matched;
+ matched.resize(100,1);
+ PlotLevel("WBP",1,1,WBP,matched);
+ PlotLevel("Exp",0,0,Exp,matched);
+
+}
+////////////////////////////////////////////////////////////////////////////////
+void PlotLevel(string name,int offset,bool LR, RS::ShellModelCollection Collection, vector<int> matched){
+ c->cd(1);
+ if(name!=""){
+ TLatex* t = new TLatex((double) offset/nbr+0.25/nbr,0.01,name.c_str());
+ t->SetTextSize(LabelSize+0.01);
+ t->Draw();
+ }
+
+ unsigned int mysize = Collection.GetNumberOfState();
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ RS::ShellModelState state = Collection.GetState(i);
+ if(Collection.GetStatus(i)==1 && state.GetNumberOfOrbital()>0){
+
+ double start, finnish;
+ if(LR){
+ start = (double)(offset)/nbr+0.05/nbr;
+ finnish = (double)(offset+1)/nbr-0.25/nbr;
+ }
+
+ else{
+ start = (double)(offset+1)/nbr-0.05/nbr;
+ finnish = (double)(offset)/nbr+0.25/nbr;
+ }
+
+ double line_length = abs(finnish-start);
+ double last = start;
+
+ double position = (state.GetEnergy()-Emin)/(Emax-Emin);
+ double labelpos ;
+
+ string final_label;
+ if(state.GetOrder()!=0){
+ if(state.GetParity()>0)
+ final_label = Form("%d^{+}_{%d} %s %.2f" ,(int)state.GetJ() ,state.GetOrder(),state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ else
+ final_label = Form("%d^{-}_{%d} %s %.2f" ,(int)state.GetJ() ,state.GetOrder(),state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ }
+ else{
+ if(state.GetParity()>0)
+ final_label = Form("%d^{+} %s %.2f" ,(int)state.GetJ() ,state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ else
+ final_label = Form("%d^{-} %s %.2f" ,(int)state.GetJ() ,state.GetOrbitalString(state.GetMainOrbital(),"SM").c_str(),state.GetEnergy());
+ }
+
+
+ if(LR)
+ labelpos = start+0.8*line_length;
+ else
+ labelpos = finnish-0.13*line_length;
+
+ TLatex* l = new TLatex(labelpos,position-0.003,final_label.c_str());
+ l->SetTextSize(LabelSize);
+ if(matched.size() ==0 || matched[i]==1)
+ l->Draw();
+
+ unsigned int orbsize = state.GetNumberOfOrbital();
+ // Where the last one stopped
+ int thickness = 4;
+ for(unsigned int orb = 0 ; orb < orbsize ; orb++){
+ TLine* Level;
+ double SF = state.GetOrbitalS(orb);
+ int color = kBlack;
+
+ string orbital = state.GetOrbitalString(orb);
+ if(orbital.find("s")!=string::npos){
+ color = kRed;
+ }
+ else if(orbital.find("p")!=string::npos){
+ color = kBlue;
+ }if(orbital.find("d")!=string::npos){
+ color = kGreen+2;
+ }if(orbital.find("f")!=string::npos){
+ color = kOrange+7;
+ }if(orbital.find("g")!=string::npos){
+ color = kMagenta+2;
+ }
+
+
+ if(LR){
+ Level =
+ new TLine(last,position,last+(SF*line_length),position);
+ last = last+(SF*line_length);
+ }
+
+ else{
+ Level =
+ new TLine(last,position,last-(SF*line_length),position);
+ last = last-(SF*line_length);
+ }
+
+ Level->SetLineStyle(1);
+ Level->SetLineWidth(thickness--);
+ Level->SetLineColor(color);
+ Level->Draw();
+ }
+ }
+ }
+
+}
+
diff --git a/Projects/S1554/macro/ShowTrifoil.cxx b/Projects/S1554/macro/ShowTrifoil.cxx
new file mode 100644
index 000000000..0d85780bc
--- /dev/null
+++ b/Projects/S1554/macro/ShowTrifoil.cxx
@@ -0,0 +1,23 @@
+
+void ShowTrifoil(){
+
+ TFile* file = new TFile("~/Desktop/S1554/Root/good/data35216_001.root");
+ TTree* tree = (TTree*) file->FindObjectAny("DataS1554");
+ TTrifoilData* tr = new TTrifoilData;
+ tree->SetBranchAddress("Trifoil",&tr);
+
+ string user;
+ int i = 0 ;
+
+ while(cin >> user ){
+ if(user == "q")
+ break;
+
+ tree->GetEntry(i++);
+ if(tr->GetMultiplicity()>0){
+ tr->GetWaveform(0)->Draw();
+ gPad->Update();
+ }
+ }
+
+}
diff --git a/Projects/S1554/macro/Si.cxx b/Projects/S1554/macro/Si.cxx
new file mode 100644
index 000000000..e0cb5611a
--- /dev/null
+++ b/Projects/S1554/macro/Si.cxx
@@ -0,0 +1,57 @@
+void Si(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Sidp_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+ tree->SetAlias("GoodTiming","Gamma_Time-Sharc.Strip_T>-1000 && Gamma_Time-Sharc.Strip_T<1000") ;
+ NPL::Reaction r("28Si(d,p)29Si@221.623");
+ NPL::Reaction d("28Si(d,d)28Si@221.623");
+ NPL::Reaction p("28Si(p,p)28Si@221.623");
+
+ TCanvas* c = new TCanvas("Result","Result",1000,1000);
+ c->Divide(3,3);
+ TH1* h;
+
+// c->cd(1);
+ tree->Draw("Ex:AddBack_DC/1000.>>h2(2000,0,8,300,0,8)","ELab>0 && ThetaLab>90&&GoodTiming && AddBack_DC@.size()==1","colz");
+ h = (TH1*) gDirectory->FindObjectAny("h2");
+ h->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
+ h->GetYaxis()->SetTitle("E_{29Si} (MeV)");
+ gPad->SetLogz();
+ TLine* line = new TLine(0,0,8,8);
+ line -> Draw("c");
+return;
+ c->cd(2);
+ tree->Draw("Ex:ThetaLab>>hEXT(180,0,180,300,-5,8)","ELab>0","colz");
+ gPad->SetLogz();
+ h = (TH1*) gDirectory->FindObjectAny("hEXT");
+ h->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ h->GetYaxis()->SetTitle("E_{29Si} (MeV)");
+
+ c->cd(3);
+ tree->Draw("ELab:ThetaLab>>hEXT2(360,0,180,1200,0,20)","ELab>0","colz");
+ gPad->SetLogz();
+ r.GetKinematicLine3()->Draw("c");
+ r.SetExcitationHeavy(4.9336);
+ r.GetKinematicLine3()->Draw("c");
+ d.GetKinematicLine3()->Draw("c");
+ p.GetKinematicLine3()->Draw("c");
+
+ h = (TH1*) gDirectory->FindObjectAny("hEXT2");
+ h->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ h->GetYaxis()->SetTitle("E_{Lab} (MeV)");
+
+ c->cd(4);
+ tree->Draw("AddBack_DC/1000>>hDC(800,0,8)","ELab>0&& ThetaLab>90 && Ex>0.8 && Ex<1.6 && AddBack_Theta >110","");
+ gPad->SetLogz();
+ h = (TH1*) gDirectory->FindObjectAny("hDC");
+ h->GetXaxis()->SetTitle("E_{#gamma} (MeV)");
+ h->GetYaxis()->SetTitle("counts / 10 keV");
+
+ c->cd(5);
+ tree->Draw("AddBack_DC/1000>>hDC2(1600,0,8)","ELab>0&& ThetaLab>90 && Ex>3 && Ex<5.5 && AddBack_Theta >110 && AddBack_Theta < 130","");
+
+ gPad->SetLogz();
+ h = (TH1*) gDirectory->FindObjectAny("hDC2");
+ h->GetYaxis()->SetTitle("E_{#gamma} (MeV)");
+ h->GetXaxis()->SetTitle("#theta_{gamma}");
+}
diff --git a/Projects/S1554/macro/SiCS.cxx b/Projects/S1554/macro/SiCS.cxx
new file mode 100644
index 000000000..3e8eddaf5
--- /dev/null
+++ b/Projects/S1554/macro/SiCS.cxx
@@ -0,0 +1,211 @@
+TGraph* thdae;
+TGraphErrors* gdd;
+double minX = 1e30;
+double minC = 1e30;
+void NumericalMinimization(const char * minName ,const char *algoName);
+double Chi2(const double* parameter);
+////////////////////////////////////////////////////////////////////////////////
+void SiCS(){
+ TFile* file = new TFile("macro/SiKine.root","READ");
+ TH2F* h2 = (TH2F*)file->FindObjectAny("hEXT4");
+
+ file = new TFile("SharcEfficiency_28Sidd.root","READ");
+ TH1D* eff = (TH1D*) file->FindObjectAny("SolidAngleCM");
+
+ file = new TFile("Si28dd_dae.root","READ");
+ thdae = (TGraph*) file->FindObjectAny("Graph");
+
+ NPL::Reaction dd("28Si(d,d)28Si@222");
+ NPL::Reaction ddp("28Si(d,d)28Si@222");
+ ddp.SetExcitationHeavy(1.779);
+ NPL::Reaction pp("28Si(p,p)28Si@222");
+
+ TGraph* kdd = dd.GetKinematicLine3();
+ TGraph* kddp = ddp.GetKinematicLine3();
+
+ TGraph* kpp = pp.GetKinematicLine3();
+ TGraph* tdd = dd.GetELabVersusThetaCM();
+ TGraph* tddp = ddp.GetELabVersusThetaCM();
+
+ TGraph* tpp = pp.GetELabVersusThetaCM();
+
+ const double* xdd= kdd->GetX();
+ const double* ydd = kdd->GetY();
+ kdd= new TGraph(kdd->GetN(),ydd,xdd);
+
+ const double* xddp= kddp->GetX();
+ const double* yddp = kddp->GetY();
+ kddp= new TGraph(kddp->GetN(),yddp,xddp);
+
+
+ const double* xpp= kpp->GetX();
+ const double* ypp = kpp->GetY();
+ kpp= new TGraph(kpp->GetN(),ypp,xpp);
+
+ const double* xdd1= tdd->GetX();
+ const double* ydd1 = tdd->GetY();
+ tdd= new TGraph(tdd->GetN()/2,xdd1,ydd1);
+
+ const double* xddp1= tddp->GetX();
+ const double* yddp1 = tddp->GetY();
+ tddp= new TGraph(tddp->GetN()/2,xddp1,yddp1);
+
+
+ const double* xpp2= tpp->GetX();
+ const double* ypp2 = tpp->GetY();
+ tpp= new TGraph(tpp->GetN()/2,xpp2,ypp2);
+
+ gdd = new TGraphErrors();
+ TGraphErrors* gpp = new TGraphErrors();
+
+ TF1* f = new TF1("f","gaus(0)+gaus(3)+gaus(6)+pol4(9)",0,90);
+ f->SetNpx(1000);
+ f->SetParameter(0,1000);
+ f->SetParameter(2,1);
+
+ f->SetParameter(3,1000);
+ f->SetParameter(5,1);
+
+ f->SetParameter(6,1000);
+ f->SetParameter(8,1);
+
+ f->SetParLimits(0,0,1e4);
+ f->SetParLimits(3,0,1e4);
+ f->SetParLimits(6,0,1e4);
+
+ f->SetParLimits(2,1,2);
+ f->SetParLimits(5,1,2);
+ f->SetParLimits(8,1,2);
+
+ TF1* g = new TF1("g","gaus(0)",0,90);
+
+ unsigned int step = 100;
+ double eps = 30./step;
+ cout << "EPSILON = " << eps <<" MeV" << endl;
+ double E=0;
+ TH1D* p;
+ double* r;
+ for(unsigned int i = 0 ; i < step ; i++){
+ E = i*eps;
+
+
+ double lTCMd = tdd->Eval(E-eps*0.5);
+ double uTCMd = tdd->Eval(E+eps*0.5);
+ p = h2->ProjectionX(Form("p%i",i),h2->GetYaxis()->FindBin(E-eps*0.5),h2->GetYaxis()->FindBin(E+eps*0.5));
+ if(p->GetEntries()>100){
+ double Td = kdd->Eval(E);
+ double Tddp = kddp->Eval(E);
+
+ double Tp = kpp->Eval(E);
+
+
+ double TCMd = tdd->Eval(E);
+ double TCMddp = tddp->Eval(E);
+
+ double TCMp = tpp->Eval(E);
+
+ cout << E << " " << Td << endl;
+ p->Draw();
+ f->SetParameter(1,Td);
+ f->SetParameter(4,Tp);
+ f->SetParameter(7,Tddp);
+
+ f->SetParLimits(1,Td,Td);
+ f->SetParLimits(4,Tp,Tp);
+ f->SetParLimits(7,Tddp,Tddp);
+
+ double pi,po;
+ f->GetParLimits(3,pi,po);
+ p->Fit(f,"RBF");
+ gPad->Update();
+
+ r = f->GetParameters();
+ g->SetParameter(0,r[0]);
+ g->SetParameter(1,r[1]);
+ g->SetParameter(2,r[2]);
+
+ double dd = g->Integral(r[1]-3*r[2],r[1]+3*r[2]);
+ double ddo = dd ;
+ double dde = sqrt(dd);
+ double mynorm = 1;
+ double EFF = eff->Integral(eff->FindBin(uTCMd),eff->FindBin(lTCMd))/abs(uTCMd-lTCMd);
+ dd = mynorm*dd / eff->GetBinContent(eff->FindBin(TCMd));
+ dde = mynorm*dde / eff->GetBinContent(eff->FindBin(TCMd));
+ double ddx = 0.5*abs(uTCMd-lTCMd);
+ if(EFF>0.0005 && E>0.8 && dd<1e50 && ddo>10 && abs(Td-Tp)>3&& abs(Td-Tddp)>3 && Td>60 && Td<80){
+ gdd->SetPoint(gdd->GetN(),TCMd,dd);
+ gdd->SetPointError(gdd->GetN()-1,ddx,dde);
+ }
+
+ g->SetParameter(0,r[3]);
+ g->SetParameter(1,r[4]);
+ g->SetParameter(2,r[5]);
+ double pp = g->Integral(0,90);
+ double ppe = sqrt(pp);
+
+ gpp->SetPoint(gpp->GetN(),TCMp,pp);
+ gpp->SetPointError(gpp->GetN()-1,0,ppe);
+ }
+ }
+
+ new TCanvas();
+ gdd->Draw("ap");
+ gdd->GetXaxis()->SetTitle("#theta_{CM} (deg)");
+ gdd->GetYaxis()->SetTitle("d#sigma/d#Omega (mb/sr)");
+ gdd->GetYaxis()->SetRangeUser(1,1e5);
+ thdae->Draw("c");
+
+ gPad->SetLogy();
+ NumericalMinimization("Minuit2","Fumili2");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Chi2(const double* parameter){
+ double Chi2 = 0 ;
+ for (int i = 0; i < gdd->GetN() ; i++) {
+ if(gdd->GetX()[i]<33 && gdd->GetX()[i]>15 && gdd->GetY()[i]>1){
+ double chi = ( thdae->Eval( gdd->GetX()[i] ) - parameter[0]*gdd->GetY()[i]) / (parameter[0]*gdd->GetEY()[i] );
+ Chi2 += chi*chi ;
+ }
+ }
+ if(Chi2<minC){
+ minC= Chi2;
+ minX= parameter[0];
+ }
+ return Chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void NumericalMinimization(const char * minName ,const char *algoName){
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+
+ // set tolerance , etc...
+ min->SetMaxFunctionCalls(100000); // for Minuit/Minuit2
+ min->SetPrecision(1e-6);
+ min->SetValidError(true);
+ min->SetPrintLevel(1);
+
+ // create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&Chi2,1);
+
+ min->SetFunction(f);
+ double* parameter = new double[1];
+ parameter[0] = 1e-4;
+ // Set the free variables to be minimized!
+ min->SetLimitedVariable(0,"Normalistation",parameter[0],1e-5,1e-5,1);
+
+ // do the minimization
+ min->Minimize();
+
+ const double *xs = min->X();
+ std::cout << "Minimum = " << xs[0] << " but MinX = " << minX << std::endl;
+
+ for(int i = 0; i < gdd->GetN() ; i++ ){
+ gdd->SetPoint(i,gdd->GetX()[i],minX*gdd->GetY()[i]);
+ gdd->SetPointError(i,gdd->GetEX()[i],minX*gdd->GetEY()[i]);
+
+ }
+ return 0;
+}
diff --git a/Projects/S1554/macro/SiCheck.cxx b/Projects/S1554/macro/SiCheck.cxx
new file mode 100644
index 000000000..cd8ad17e6
--- /dev/null
+++ b/Projects/S1554/macro/SiCheck.cxx
@@ -0,0 +1,20 @@
+void SiCheck(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Si_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ int current=16;
+ TCanvas* c;
+ for(unsigned int i = 35090 ; i < 35300 ; i++){
+
+ if(current == 16){
+ c = new TCanvas(Form("Check%i",i),Form("Check%i",i),900,900);
+ c->Divide(4,4);
+ current = 1;
+ }
+ c->cd(current++);
+ tree->Draw(Form("ELab:ThetaLab>>hK%i(180,0,180,150,0,10)",i),Form("ELab>0&&RunNumber == %i",i),"colz");
+ gPad->SetLogz();
+ c->Update();
+ }
+}
diff --git a/Projects/S1554/macro/SiCheck2.cxx b/Projects/S1554/macro/SiCheck2.cxx
new file mode 100644
index 000000000..254a8dd57
--- /dev/null
+++ b/Projects/S1554/macro/SiCheck2.cxx
@@ -0,0 +1,22 @@
+void SiCheck2(){
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Si_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+ NPL::Reaction r("28Si(d,p)29Si@222");
+ r.SetExcitationHeavy(4.9336);
+ TGraph* line1 = r.GetKinematicLine3();
+ tree->SetAlias("GoodEvent","ELab>0 && AddBack_DC/1000. >4.6 && AddBack_DC/1000.<5.4 && Ex<6.5");
+ new TCanvas();
+ tree->Draw(Form("ELab:ThetaLab>>hEX%i(180,0,180,400,0,10)",5)," DetectorNumber == 5 && GoodEvent","colz");
+ line1->Draw("c");
+ new TCanvas();
+ tree->Draw(Form("ELab:ThetaLab>>hEX%i(180,0,180,400,0,10)",6),"DetectorNumber == 6 && GoodEvent","colz");
+ line1->Draw("c");
+ new TCanvas();
+ tree->Draw(Form("ELab:ThetaLab>>hEX%i(180,0,180,400,0,10)",7),"ELab>0 && DetectorNumber == 7 && GoodEvent","colz");
+ line1->Draw("c");
+ new TCanvas();
+ tree->Draw(Form("ELab:ThetaLab>>hEX%i(180,0,180,400,0,10)",8),"ELab>0 && DetectorNumber == 8 && GoodEvent","colz");
+ line1->Draw("c");
+}
+
+
diff --git a/Projects/S1554/macro/SiCheck3.cxx b/Projects/S1554/macro/SiCheck3.cxx
new file mode 100644
index 000000000..4525b61c8
--- /dev/null
+++ b/Projects/S1554/macro/SiCheck3.cxx
@@ -0,0 +1,19 @@
+void SiCheck3(){
+
+ TFile* file = new TFile("../../../Outputs/Analysis/S1554_28Si_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ TCanvas* c = new TCanvas("Result","Result",900,900);
+ c->Divide(4,4);
+
+ string hist,cond;
+ for(unsigned int i = 0 ; i < 16 ; i ++){
+ c->cd(i+1);
+ hist = Form("Ex:Gamma_Energy/1000.>>h%i(2000,0,10,300,-5,15)",i+1);
+ cond = Form("ELab>0 && ThetaLab>90 && Clover_Number==%i",i+1);
+ tree->Draw(hist.c_str(),cond.c_str(),"colz");
+ gPad->SetLogz();
+ }
+
+
+}
diff --git a/Projects/S1554/macro/SiCheck4.cxx b/Projects/S1554/macro/SiCheck4.cxx
new file mode 100644
index 000000000..cb75410b1
--- /dev/null
+++ b/Projects/S1554/macro/SiCheck4.cxx
@@ -0,0 +1,16 @@
+void SiCheck4(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Si_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ TCanvas* c = new TCanvas("Result","Result",900,900);
+ c->Divide(4,4);
+
+ for(unsigned int i = 0 ; i < 16 ; i++){
+ c->cd(i+1);
+ if(i+1>4)
+ tree->Draw(Form("Gamma_Energy/1000.>>h%i(10000,0,10)",i),Form("ELab>0 && ThetaLab>90 && Clover_Number==%i",i+1),"colz");
+ c->Update();
+ }
+
+}
diff --git a/Projects/S1554/macro/Sidd.cxx b/Projects/S1554/macro/Sidd.cxx
new file mode 100644
index 000000000..7a4547419
--- /dev/null
+++ b/Projects/S1554/macro/Sidd.cxx
@@ -0,0 +1,44 @@
+void Sidd(){
+
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Sidd_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ NPL::Reaction r("28Si(d,d)28Si@224");
+ NPL::Reaction p("28Si(p,p)28Si@224");
+
+ // Setting Trifoil entry list for fast drawing
+ // Load Cuts
+ TFile* fileCut = new TFile("cuts/Cut_dd.root");
+ TCutG* Cut_dd = (TCutG*) fileCut->FindObjectAny("Cut_dd");
+
+ TFile* fileCut1 = new TFile("cuts/cut_elgs_det10.root");
+ TCutG* Cut_dd10 = (TCutG*) fileCut1->FindObjectAny("cut_elgs_det10");
+
+ TFile* fileCut2 = new TFile("cuts/cut_elgs_det12.root");
+ TCutG* Cut_dd12 = (TCutG*) fileCut2->FindObjectAny("cut_elgs_det12");
+
+ // Create Canvas
+ TCanvas* c = new TCanvas("Result","Result",900,900);
+ c->Divide(2,2);
+ TH1* h;
+ TH1F* h2;
+
+ c->cd(1);
+ tree->Draw("ELab:ThetaLab>>hEXT4(360,0,180,1200,0,40)","ELab>0 ","colz");
+
+ gPad->SetLogz();
+ h2 = (TH1F*) gDirectory->FindObjectAny("hEXT4");
+ h2->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ h2->GetYaxis()->SetTitle("E_{Lab} (MeV)");
+ r.GetKinematicLine3()->Draw("c");
+ p.GetKinematicLine3()->Draw("c");
+
+ c->cd(2);
+ tree->Draw("ThetaCM>>hCM(180,0,180)","ELab>0 && (Cut_dd||Sharc.PAD_E<0) && Ex>-0.8 && Ex<0.8","");
+
+ c->cd(3);
+ tree->Draw("ELab:ThetaLab>>hEXTA(360,0,180,1200,0,40)","ELab>0 && (Cut_dd||Sharc.PAD_E<0) && Sharc.DetectorNumber==12","colz");
+ c->cd(4);
+ tree->Draw("ELab:ThetaLab>>hEXTB(360,0,180,1200,0,40)","ELab>0 && (Cut_dd||Sharc.PAD_E<0)&& Sharc.DetectorNumber ==10","colz");
+
+}
diff --git a/Projects/S1554/macro/makeMgCS.cxx b/Projects/S1554/macro/makeMgCS.cxx
new file mode 100644
index 000000000..064a39628
--- /dev/null
+++ b/Projects/S1554/macro/makeMgCS.cxx
@@ -0,0 +1,89 @@
+////////////////////////////////////////////////////////////////////////////////
+void makeMgCS(){
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Mgdp_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ tree->SetAlias("GoodEvent","(RunNumber>35110 && ELab > 0)");
+ tree->SetAlias("GoodTrifoil","(Trifoil.Time@.size()==1 && Trifoil.Time>145 && Trifoil.Time<160)") ;
+ tree->SetAlias("BadTrifoil","(Trifoil.Time@.size()==1 && (Trifoil.Time<145 || Trifoil.Time>160) && Trifoil.Time<300) ") ;
+ tree->SetAlias("GoodProton","(ThetaLab>90 &&( (Sharc.Strip_E>0.8&&ThetaCM<10)||(Sharc.Strip_E>1 &&ThetaCM>10) ))");
+
+ TCanvas* c0 = new TCanvas("Mg0","Mg0",800,800*2./5.);
+ c0->Divide(5,2);
+ c0->cd(1);
+ tree->Draw("Ex:ThetaCM>>hexcmG(72,0,180,400,-8,8)","GoodTrifoil && GoodProton && GoodEvent","colz");
+ TH2D* h = (TH2D*) gDirectory->FindObjectAny("hexcmG");
+ h->GetYaxis()->SetTitle("GOOD E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("GOOD #theta_{CM} (deg)");
+ h->GetXaxis()->SetRangeUser(-1,50);
+ h->SaveAs("GoodExThetaCM.root");
+
+ c0->cd(2);
+ tree->Draw("Ex:ThetaCM>>hexcmB(72,0,180,400,-8,8)","BadTrifoil && GoodProton && GoodEvent","colz");
+ h = (TH2D*) gDirectory->FindObjectAny("hexcmB");
+ h->GetYaxis()->SetTitle("BAD E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("BAD #theta_{CM} (deg)");
+ h->GetXaxis()->SetRangeUser(-1,50);
+ // h->Scale(7./(145-1+300-160));
+ h->SaveAs("BadExThetaCM.root");
+
+ c0->cd(3);
+ tree->Draw("Ex:ThetaCM>>hexcmB2(72,0,180,400,-8,8)","Trifoil.Time@.size()==0 && GoodProton && GoodEvent","colz");
+ h = (TH2D*) gDirectory->FindObjectAny("hexcmB2");
+ h->GetYaxis()->SetTitle("BAD E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("BAD #theta_{CM} (deg)");
+ h->GetXaxis()->SetRangeUser(-1,50);
+ // h->Scale(7./(145-1+300-160));
+ h->SaveAs("BadExThetaCM2.root");
+
+/*
+ c0->cd(4);
+ tree->Draw("Ex>>hexG(400,-8,8)","GoodTrifoil && GoodProton && GoodEvent","");
+ TH1D* hh = (TH1D*) gDirectory->FindObjectAny("hexG");
+ h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
+ tree->Draw("Ex>>hexB(400,-8,8)","BadTrifoil && GoodProton && GoodEvent","same");
+ hh = (TH1D*) gDirectory->FindObjectAny("hexB");
+ hh->SetFillColor(kOrange+7);
+ hh->SetLineColor(kOrange+7);
+ hh->SetFillStyle(3004);
+// hh->Scale(7./(145-1+300-160));
+
+ c0->cd(5);
+ tree->Draw("Trifoil.Time>>htt(300,0,300)","GoodProton && GoodEvent","");
+ tree->Draw("Trifoil.Time>>hgt(300,0,300)","GoodTrifoil && GoodProton && GoodEvent ","same");
+ tree->Draw("Trifoil.Time>>hbt(300,0,300)","BadTrifoil && GoodProton && GoodEvent ","same");
+
+ TH1D* hgt = (TH1D*) gDirectory->FindObjectAny("hgt");
+ hgt->SetFillColor(kOrange+7);
+ hgt->SetLineColor(kOrange+7);
+ hgt->SetFillStyle(3004);
+ TH1D* hbt = (TH1D*) gDirectory->FindObjectAny("hbt");
+ hbt->SetFillColor(kGreen-3);
+ hbt->SetLineColor(kGreen-3);
+ hbt->SetFillStyle(3004);
+
+
+
+ c0->cd(6);
+ tree->Draw("Trifoil.Time:ThetaLab>>hte(360,0,360,300,0,300)","GoodProton && GoodEvent","colz");
+ TH1D* hte = (TH1D*) gDirectory->FindObjectAny("hte");
+ hte->GetXaxis()->SetTitle("Sharc E (MeV)");
+ hte->GetYaxis()->SetTitle("Trifoil Time (a.u.)");
+
+ c0->cd(7);
+ tree->Draw("Ex>>hexSi(400,-8,8)","Trifoil.Time@.size()==0 && GoodProton && GoodEvent","");
+ hgt = (TH1D*) gDirectory->FindObjectAny("hexSi");
+ TH1D* hhh = new TH1D(*hh);
+ hhh->Scale(hgt->GetEntries()/hhh->GetEntries());
+ hhh->Draw("same");
+
+ c0->cd(8);
+ tree->Draw("Ex>>hexSi2(400,-8,8)","Trifoil.Time@.size()==0 && GoodProton && GoodEvent && AddBack_DC/1000.>4 && AddBack_DC/1000.<6","");
+
+ // new TCanvas();
+// tree->Draw("AddBack_DC>>hGAL(2000,0,8000)","GoodTrifoil && Ex<0 && Ex>-1000 && ThetaLab>90");
+new TCanvas();
+// tree->Draw("Ex:ThetaCM>>hexcmG(90,0,180,400,-8,8)","GoodTrifoil && GoodProton && GoodEvent","colz");
+ tree->Draw("Ex","GoodTrifoil && GoodProton && GoodEvent && AddBack_DC/1000.> 1 && AddBack_DC/1000.<1.2 && Ex>-1000","");
+ */
+}
diff --git a/Projects/S1554/macro/makeMgCSLab.cxx b/Projects/S1554/macro/makeMgCSLab.cxx
new file mode 100644
index 000000000..903abdcbd
--- /dev/null
+++ b/Projects/S1554/macro/makeMgCSLab.cxx
@@ -0,0 +1,89 @@
+////////////////////////////////////////////////////////////////////////////////
+void makeMgCSLab(){
+ TFile* file = new TFile("../../Outputs/Analysis/S1554_28Mgdp_Phy.root");
+ TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+
+ tree->SetAlias("GoodEvent","(RunNumber>35110 && ELab > 0)");
+ tree->SetAlias("GoodTrifoil","(Trifoil.Time@.size()==1 && Trifoil.Time>145 && Trifoil.Time<160)") ;
+ tree->SetAlias("BadTrifoil","(Trifoil.Time@.size()==1 && (Trifoil.Time<145 || Trifoil.Time>160) && Trifoil.Time<300) ") ;
+ tree->SetAlias("GoodProton","(ThetaLab>90 && Sharc.Strip_E>0.8)");
+
+ TCanvas* c0 = new TCanvas("Mg0","Mg0",800,800*2./5.);
+ c0->Divide(5,2);
+ c0->cd(1);
+ tree->Draw("Ex:ThetaLab>>hexcmG(72,0,180,400,-8,8)","GoodTrifoil && GoodProton && GoodEvent","colz");
+ TH2D* h = (TH2D*) gDirectory->FindObjectAny("hexcmG");
+ h->GetYaxis()->SetTitle("GOOD E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("GOOD #theta_{Lab} (deg)");
+ h->GetXaxis()->SetRangeUser(90,180);
+ h->SaveAs("GoodExThetaLab.root");
+
+ c0->cd(2);
+ tree->Draw("Ex:ThetaLab>>hexcmB(72,0,180,400,-8,8)","BadTrifoil && GoodProton && GoodEvent","colz");
+ h = (TH2D*) gDirectory->FindObjectAny("hexcmB");
+ h->GetYaxis()->SetTitle("BAD E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("BAD #theta_{Lab} (deg)");
+ h->GetXaxis()->SetRangeUser(90,180);
+ // h->Scale(7./(145-1+300-160));
+ h->SaveAs("BadExThetaLab.root");
+
+ c0->cd(3);
+ tree->Draw("Ex:ThetaLab>>hexcmB2(72,0,180,400,-8,8)","Trifoil.Time@.size()==0 && GoodProton && GoodEvent","colz");
+ h = (TH2D*) gDirectory->FindObjectAny("hexcmB2");
+ h->GetYaxis()->SetTitle("BAD E_{29Mg} (MeV)");
+ h->GetXaxis()->SetTitle("BAD #theta_{Lab} (deg)");
+ h->GetXaxis()->SetRangeUser(90,180);
+ // h->Scale(7./(145-1+300-160));
+ h->SaveAs("BadExThetaLab2.root");
+
+/*
+ c0->cd(4);
+ tree->Draw("Ex>>hexG(400,-8,8)","GoodTrifoil && GoodProton && GoodEvent","");
+ TH1D* hh = (TH1D*) gDirectory->FindObjectAny("hexG");
+ h->GetXaxis()->SetTitle("E_{29Mg} (MeV)");
+ tree->Draw("Ex>>hexB(400,-8,8)","BadTrifoil && GoodProton && GoodEvent","same");
+ hh = (TH1D*) gDirectory->FindObjectAny("hexB");
+ hh->SetFillColor(kOrange+7);
+ hh->SetLineColor(kOrange+7);
+ hh->SetFillStyle(3004);
+// hh->Scale(7./(145-1+300-160));
+
+ c0->cd(5);
+ tree->Draw("Trifoil.Time>>htt(300,0,300)","GoodProton && GoodEvent","");
+ tree->Draw("Trifoil.Time>>hgt(300,0,300)","GoodTrifoil && GoodProton && GoodEvent ","same");
+ tree->Draw("Trifoil.Time>>hbt(300,0,300)","BadTrifoil && GoodProton && GoodEvent ","same");
+
+ TH1D* hgt = (TH1D*) gDirectory->FindObjectAny("hgt");
+ hgt->SetFillColor(kOrange+7);
+ hgt->SetLineColor(kOrange+7);
+ hgt->SetFillStyle(3004);
+ TH1D* hbt = (TH1D*) gDirectory->FindObjectAny("hbt");
+ hbt->SetFillColor(kGreen-3);
+ hbt->SetLineColor(kGreen-3);
+ hbt->SetFillStyle(3004);
+
+
+
+ c0->cd(6);
+ tree->Draw("Trifoil.Time:ThetaLab>>hte(360,0,360,300,0,300)","GoodProton && GoodEvent","colz");
+ TH1D* hte = (TH1D*) gDirectory->FindObjectAny("hte");
+ hte->GetXaxis()->SetTitle("Sharc E (MeV)");
+ hte->GetYaxis()->SetTitle("Trifoil Time (a.u.)");
+
+ c0->cd(7);
+ tree->Draw("Ex>>hexSi(400,-8,8)","Trifoil.Time@.size()==0 && GoodProton && GoodEvent","");
+ hgt = (TH1D*) gDirectory->FindObjectAny("hexSi");
+ TH1D* hhh = new TH1D(*hh);
+ hhh->Scale(hgt->GetEntries()/hhh->GetEntries());
+ hhh->Draw("same");
+
+ c0->cd(8);
+ tree->Draw("Ex>>hexSi2(400,-8,8)","Trifoil.Time@.size()==0 && GoodProton && GoodEvent && AddBack_DC/1000.>4 && AddBack_DC/1000.<6","");
+
+ // new TCanvas();
+// tree->Draw("AddBack_DC>>hGAL(2000,0,8000)","GoodTrifoil && Ex<0 && Ex>-1000 && ThetaLab>90");
+new TCanvas();
+// tree->Draw("Ex:ThetaLab>>hexcmG(90,0,180,400,-8,8)","GoodTrifoil && GoodProton && GoodEvent","colz");
+ tree->Draw("Ex","GoodTrifoil && GoodProton && GoodEvent && AddBack_DC/1000.> 1 && AddBack_DC/1000.<1.2 && Ex>-1000","");
+ */
+}
diff --git a/Projects/S1554/macro/pipo.cxx b/Projects/S1554/macro/pipo.cxx
new file mode 100644
index 000000000..4fe13d211
--- /dev/null
+++ b/Projects/S1554/macro/pipo.cxx
@@ -0,0 +1,11 @@
+void pipo(){
+
+
+ double* x = new double[2];
+ double* y = new double[2];
+ x[0] = 0.027; y[0] = -0.029;
+ x[1] = 1.431; y[1] = 1.376;
+
+TGraph* g= new TGraph(2,x,y);
+g->Draw("ac");
+}
diff --git a/Projects/S1554/macro/test.cxx b/Projects/S1554/macro/test.cxx
new file mode 100644
index 000000000..595f63869
--- /dev/null
+++ b/Projects/S1554/macro/test.cxx
@@ -0,0 +1,25 @@
+{
+ c0 = new TCanvas("c1","multigraph L3",200,10,700,500);
+ // c0->SetFrameFillColor(30);
+ TMultiGraph *mg = new TMultiGraph();
+ TGraph *gr1 = new TGraph(); gr1->SetLineColor(kBlue);
+ TGraph *gr2 = new TGraph(); gr2->SetLineColor(kRed);
+ TGraph *gr3 = new TGraph(); gr3->SetLineColor(kGreen);
+ TGraph *gr4 = new TGraph(); gr4->SetLineColor(kOrange);
+ Double_t dx = 6.28/100;
+ Double_t x = -3.14;
+ for (int i=0; i<=100; i++) {
+ x = x+dx;
+ gr1->SetPoint(i,x,2.*TMath::Sin(x));
+ gr2->SetPoint(i,x,TMath::Cos(x));
+ gr3->SetPoint(i,x,TMath::Cos(x*x));
+ gr4->SetPoint(i,x,TMath::Cos(x*x*x));
+ }
+ mg->Add(gr4); gr4->SetTitle("Cos(x*x*x)"); gr4->SetLineWidth(3);
+ mg->Add(gr3); gr3->SetTitle("Cos(x*x)") ; gr3->SetLineWidth(3);
+ mg->Add(gr2); gr2->SetTitle("Cos(x)") ; gr2->SetLineWidth(3);
+ mg->Add(gr1); gr1->SetTitle("2*Sin(x)") ; gr1->SetLineWidth(3);
+
+ mg->Draw("a l3d");
+ return c0;
+}
diff --git a/Projects/S1554/mg29olapc.lsf b/Projects/S1554/mg29olapc.lsf
new file mode 100644
index 000000000..120a18e02
--- /dev/null
+++ b/Projects/S1554/mg29olapc.lsf
@@ -0,0 +1,94 @@
+! model space = spsdpf
+! interaction = spsdpfpc
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 1 0.3966 -265.061 -268.637 -1.418 0.000 0.000
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 2 0.0004 -265.061 -265.732 0.000 0.000 2.905
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 3 0.0001 -265.061 -263.462 0.000 0.000 5.175
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 4 0.0144 -265.061 -262.423 0.038 0.000 6.214
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 5 0.0015 -265.061 -261.565 0.005 0.000 7.072
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 6 0.0000 -265.061 -261.240 0.000 0.000 7.397
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 7 0.0033 -265.061 -261.196 0.013 0.000 7.441
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 8 0.0004 -265.061 -260.558 0.002 0.000 8.079
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 9 0.0000 -265.061 -260.194 0.000 0.000 8.443
+( 28 2.0) ( 29 2.5) ( n 2 0 1) 0.0+ 0.5+ 1 10 0.0040 -265.061 -260.049 0.020 0.000 8.588
+ sum 0.4207 -1.341
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 1 0.4013 -265.061 -268.547 -1.399 0.000 0.090
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 2 0.2505 -265.061 -266.368 -0.327 0.000 2.269
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 3 0.0043 -265.061 -265.018 0.000 0.000 3.619
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 4 0.0000 -265.061 -263.969 0.000 0.000 4.668
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 5 0.0009 -265.061 -263.018 0.002 0.000 5.619
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 6 0.0076 -265.061 -262.692 0.018 0.000 5.945
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 7 0.0033 -265.061 -262.567 0.008 0.000 6.070
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 8 0.0019 -265.061 -262.264 0.005 0.000 6.373
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 9 0.0001 -265.061 -262.004 0.000 0.000 6.633
+( 28 2.0) ( 29 2.5) ( n 1 2 3) 0.0+ 1.5+ 1 10 0.0009 -265.061 -261.676 0.003 0.000 6.961
+ sum 0.6708 -1.689
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 1 0.0023 -265.061 -267.026 -0.005 0.000 1.611
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 2 0.0544 -265.061 -265.490 -0.023 0.000 3.147
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 3 0.0000 -265.061 -265.009 0.000 0.000 3.628
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 4 0.0008 -265.061 -264.384 0.001 0.000 4.253
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 5 0.0003 -265.061 -263.477 0.000 0.000 5.160
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 6 0.0001 -265.061 -263.238 0.000 0.000 5.399
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 7 0.0003 -265.061 -263.035 0.001 0.000 5.602
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 8 0.0010 -265.061 -262.724 0.002 0.000 5.913
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 9 0.0007 -265.061 -262.411 0.002 0.000 6.226
+( 28 2.0) ( 29 2.5) ( n 1 2 5) 0.0+ 2.5+ 1 10 0.0000 -265.061 -262.054 0.000 0.000 6.583
+ sum 0.0599 -0.022
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 1 0.3037 -265.061 -266.216 -0.351 0.000 2.421
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 2 0.4153 -265.061 -264.991 0.029 0.000 3.646
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 3 0.0001 -265.061 -263.484 0.000 0.000 5.153
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 4 0.0424 -265.061 -263.179 0.080 0.000 5.458
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 5 0.0004 -265.061 -262.953 0.001 0.000 5.684
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 6 0.0178 -265.061 -262.364 0.048 0.000 6.273
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 7 0.0056 -265.061 -261.866 0.018 0.000 6.771
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 8 0.0004 -265.061 -261.743 0.001 0.000 6.894
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 9 0.0002 -265.061 -261.575 0.001 0.000 7.062
+( 28 2.0) ( 29 2.5) ( n 2 1 1) 0.0+ 0.5- 1 10 0.0011 -265.061 -261.161 0.004 0.000 7.476
+ sum 0.7870 -0.169
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 1 0.6262 -265.061 -267.287 -1.394 0.000 1.350
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 2 0.0121 -265.061 -265.157 -0.001 0.000 3.480
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 3 0.0848 -265.061 -264.664 0.034 0.000 3.973
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 4 0.0000 -265.061 -263.763 0.000 0.000 4.874
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 5 0.0497 -265.061 -263.166 0.094 0.000 5.471
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 6 0.0818 -265.061 -262.932 0.174 0.000 5.705
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 7 0.0003 -265.061 -262.625 0.001 0.000 6.012
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 8 0.0004 -265.061 -262.443 0.001 0.000 6.194
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 9 0.0085 -265.061 -262.363 0.023 0.000 6.274
+( 28 2.0) ( 29 2.5) ( n 2 1 3) 0.0+ 1.5- 1 10 0.0009 -265.061 -262.154 0.003 0.000 6.483
+ sum 0.8647 -1.066
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 1 0.0253 -265.061 -265.564 -0.013 0.000 3.073
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 2 0.0183 -265.061 -264.274 0.014 0.000 4.363
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 3 0.0325 -265.061 -263.520 0.050 0.000 5.117
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 4 0.0095 -265.061 -263.298 0.017 0.000 5.339
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 5 0.0106 -265.061 -263.070 0.021 0.000 5.567
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 6 0.0005 -265.061 -262.853 0.001 0.000 5.784
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 7 0.0248 -265.061 -262.718 0.058 0.000 5.919
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 8 0.0179 -265.061 -262.521 0.045 0.000 6.116
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 9 0.0000 -265.061 -262.204 0.000 0.000 6.433
+( 28 2.0) ( 29 2.5) ( n 1 3 5) 0.0+ 2.5- 1 10 0.0000 -265.061 -262.083 0.000 0.000 6.554
+ sum 0.1394 0.194
+ Exi Exf
+( Ai Tzi) ( Af Tzf) (type n,l,2j) Ji Jf ni nf C^2S Ei Ef C^2S*(Ef-Ei) Exi Exf
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 1 0.4246 -265.061 -266.770 -0.726 0.000 1.867
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 2 0.3390 -265.061 -264.480 0.197 0.000 4.157
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 3 0.0000 -265.061 -263.731 0.000 0.000 4.906
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 4 0.0093 -265.061 -263.190 0.017 0.000 5.447
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 5 0.0007 -265.061 -262.999 0.001 0.000 5.638
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 6 0.0096 -265.061 -262.840 0.021 0.000 5.797
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 7 0.0011 -265.061 -262.559 0.003 0.000 6.078
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 8 0.0011 -265.061 -262.263 0.003 0.000 6.374
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 9 0.0103 -265.061 -262.037 0.031 0.000 6.600
+( 28 2.0) ( 29 2.5) ( n 1 3 7) 0.0+ 3.5- 1 10 0.0007 -265.061 -261.834 0.002 0.000 6.803
+ sum 0.7964 -0.449
+ total sum 3.7389
diff --git a/Projects/SharcEfficiency/Analysis.cxx b/Projects/SharcEfficiency/Analysis.cxx
new file mode 100644
index 000000000..be0550555
--- /dev/null
+++ b/Projects/SharcEfficiency/Analysis.cxx
@@ -0,0 +1,308 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2014 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: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2015 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include<iostream>
+using namespace std;
+#include"Analysis.h"
+#include"NPAnalysisFactory.h"
+#include"NPDetectorManager.h"
+#include"NPOptionManager.h"
+#include "NPFunction.h"
+#include"TF1.h"
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis(){
+}
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init(){
+ InitOutputBranch();
+ InitInputBranch();
+ Sharc = (TSharcPhysics*) m_DetectorManager -> GetDetector("Sharc");
+ Trifoil = (TPlasticPhysics*) m_DetectorManager -> GetDetector("Plastic");
+
+ myReaction = new NPL::Reaction();
+ myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+
+ // target thickness
+ TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
+ string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+
+ // energy losses
+ string light=NPL::ChangeNameToG4Standard(myReaction->GetNucleus3().GetName());
+ string beam=NPL::ChangeNameToG4Standard(myReaction->GetNucleus1().GetName());
+
+ LightCD2 = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
+ LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100);
+ LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
+ BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
+
+ OriginalBeamEnergy = myReaction->GetBeamEnergy();
+ Rand = TRandom3();
+ DetectorNumber = 0 ;
+ ThetaNormalTarget = 0 ;
+ Energy = 0;
+
+ RunNumber = 0;
+ RunNumberMinor=0;
+
+ ThetaSharcSurface = 0;
+ X_Sharc = 0 ;
+ Y_Sharc = 0 ;
+ Z_Sharc = 0 ;
+ X_Trifoil = 0;
+ Y_Trifoil = 0 ;
+
+
+ Si_E_Sharc = 0 ;
+ E_Sharc = 0;
+ ThetaDetector = 0 ;
+ BeamDirection = TVector3(0,0,1);
+ TargetPosition = TVector3(0.1635909,0.910980,m_DetectorManager->GetTargetZ() );
+ double finalEnergy = BeamCD2.Slow(224,TargetThickness*0.5,0);
+ myReaction->SetBeamEnergy(finalEnergy);
+ cout << "Set Beam energy to: " << finalEnergy << " MeV" << endl;
+
+ myInit = new TInitialConditions();
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&myInit);
+
+ HistoFile = new TFile("SharcEfficiency.root","RECREATE");
+ ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",72,0,180);
+ ThetaCM_detected = new TH1F("ThetaCM_detected","ThetaCM_detected",72,0,180);
+ ThetaLab_emmitted = new TH1F("ThetaLab_emmitted","ThetaLab_emmitted",72,0,180);
+ ThetaLab_detected = new TH1F("ThetaLab_detected","ThetaLab_detected",72,0,180);
+
+ ThetaCM_emmitted_2D = new TH2F("ThetaCM_emmitted_2D","ThetaCM_emmitted_2D",72,0,180,400,-8,8);
+ ThetaCM_detected_2D = new TH2F("ThetaCM_detected_2D","ThetaCM_detected_2D",72,0,180,400,-8,8);
+ ThetaLab_emmitted_2D = new TH2F("ThetaLab_emmitted_2D","ThetaLab_emmitted_2D",72,0,180,400,-8,8);
+ ThetaLab_detected_2D = new TH2F("ThetaLab_detected_2D","ThetaLab_detected_2D",72,0,180,400,-8,8);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent(){
+ // Reinitiate calculated variable
+ ReInitValue();
+ ////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////// LOOP on Sharc//////////////////
+ // Fill Initial condition Histo
+ ThetaCM_emmitted->Fill(myInit->GetThetaCM(0));
+ ThetaLab_emmitted->Fill(myInit->GetThetaLab_WorldFrame(0));
+
+ double EXD =myReaction->ReconstructRelativistic(myInit->GetKineticEnergy(0),myInit->GetThetaLab_WorldFrame(0)*deg);
+ ThetaCM_emmitted_2D->Fill(myInit->GetThetaCM(0),EXD);
+ ThetaLab_emmitted_2D->Fill(myInit->GetThetaLab_WorldFrame(0),EXD);
+
+ if(Sharc->Strip_E.size()==1){
+ /************************************************/
+ // Part 1 : Impact Angle
+ TVector3 HitDirection = Sharc -> GetPositionOfInteraction(0,true)-TargetPosition;
+ X_Sharc = Sharc -> GetPositionOfInteraction(0,true).X();
+ Y_Sharc = Sharc -> GetPositionOfInteraction(0,true).Y();
+ Z_Sharc = Sharc -> GetPositionOfInteraction(0,true).Z();
+
+ ThetaLab = HitDirection.Angle( BeamDirection );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ ThetaDetector = HitDirection.Angle(-Sharc->GetDetectorNormal(0));
+ /************************************************/
+
+ /************************************************/
+ // Part 2 : Impact Energy
+ Energy = ELab = 0;
+ if(Sharc->PAD_E[0]>0 && Sharc->PAD_E[0] <100){
+ Energy = Sharc->PAD_E[0];
+ }
+ if(Sharc->Strip_E[0]>0)
+ Energy += Sharc->Strip_E[0];
+
+ Energy = LightAl.EvaluateInitialEnergy(Energy,Sharc->GetDeadLayer(0)*micrometer,0);
+
+ // Target Correction
+ ELab = Energy;
+ ELab = LightCD2.EvaluateInitialEnergy( Energy ,TargetThickness*0.5, ThetaNormalTarget);
+
+ /************************************************/
+
+ /************************************************/
+ // Part 3 : Excitation Energy Calculation
+ Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
+ /************************************************/
+
+ /************************************************/
+ // Part 4 : Theta CM Calculation
+ ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ ThetaLab=ThetaLab/deg;
+ /************************************************/
+ // Part 5: Compute the X and Y coordinate of the Heavy in the Trifoil plan
+ myReaction->SetThetaCM(ThetaCM*deg);
+ myReaction->SetExcitationHeavy(Ex);
+
+ double d1,d2,d3,t;
+ myReaction->KineRelativistic(d1,d2,t,d3);
+ double x = tan(t)*400;
+ TVector3 P;
+ P.SetXYZ(x,0,400);
+ P.RotateZ(HitDirection.Phi()+3.14159);
+ X_Trifoil = P.X();
+ Y_Trifoil = P.Y();
+
+ myReaction->SetThetaCM(0);
+ myReaction->SetExcitationHeavy(0);
+ if(Trifoil->Energy.size()>0){
+ if( abs(Ex-EXD)<0.5 && Trifoil->Energy[0]>0 && Sharc->Strip_E[0]>0.8) {
+ // if( abs(Ex-EXD)<0.5 && Trifoil->Energy.size()>0 && Trifoil->Energy[0]>0 /*&& ((Sharc->Strip_E[0]>2 && Sharc->DetectorNumber[0]!=1) ||(Sharc->Strip_E[0]>0.7 && Sharc->DetectorNumber[0]==1))*/ ) { // for S1107
+
+ ThetaCM_detected->Fill(myInit->GetThetaCM(0));
+ ThetaLab_detected->Fill(myInit->GetThetaLab_WorldFrame(0));
+ ThetaCM_detected_2D->Fill(myInit->GetThetaCM(0),Ex);
+ ThetaLab_detected_2D->Fill(myInit->GetThetaLab_WorldFrame(0),Ex);
+ }
+ }
+ /************************************************/
+ }//end loop Sharc
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End(){
+ Efficiency_CM = new TH1F(*ThetaCM_detected);
+ Efficiency_CM->SetName("EfficiencyCM");
+ Efficiency_CM->SetTitle("EfficiencyCM");
+ Efficiency_CM->Sumw2();
+
+ Efficiency_Lab = new TH1F(*ThetaLab_detected);
+ Efficiency_Lab->SetName("EfficiencyLab");
+ Efficiency_Lab->SetTitle("EfficiencyLab");
+ Efficiency_Lab->Sumw2();
+
+ Efficiency_CM_2D = new TH2F(*ThetaCM_detected_2D);
+ Efficiency_CM_2D->SetName("EfficiencyCM_2D");
+ Efficiency_CM_2D->SetTitle("EfficiencyCM_2D");
+ Efficiency_CM_2D->Sumw2();
+
+ Efficiency_Lab_2D = new TH2F(*ThetaLab_detected_2D);
+ Efficiency_Lab_2D->SetName("EfficiencyLab_2D");
+ Efficiency_Lab_2D->SetTitle("EfficiencyLab_2D");
+ Efficiency_Lab_2D->Sumw2();
+
+
+ ThetaCM_detected->Sumw2();
+ ThetaLab_detected->Sumw2();
+
+ SolidAngle_CM = new TH1F(*ThetaCM_detected);
+ SolidAngle_CM->SetName("SolidAngleCM");
+ SolidAngle_CM->SetTitle("SolidAngleCM");
+
+ SolidAngle_Lab = new TH1F(*ThetaLab_detected);
+ SolidAngle_Lab->SetName("SolidAngleLab");
+ SolidAngle_Lab->SetTitle("SolidAngleLab");
+
+ SolidAngle_CM_2D = new TH2F(*ThetaCM_detected_2D);
+ SolidAngle_CM_2D->SetName("SolidAngleCM_2D");
+ SolidAngle_CM_2D->SetTitle("SolidAngleCM_2D");
+ SolidAngle_CM_2D->Sumw2();
+
+ SolidAngle_Lab_2D = new TH2F(*ThetaLab_detected_2D);
+ SolidAngle_Lab_2D->SetName("SolidAngleLab_2D");
+ SolidAngle_Lab_2D->SetTitle("SolidAngleLab_2D");
+ SolidAngle_Lab_2D->Sumw2();
+
+ Efficiency_CM->Divide(ThetaCM_emmitted);
+ Efficiency_Lab->Divide(ThetaLab_emmitted);
+
+ Efficiency_CM_2D->Divide(ThetaCM_emmitted_2D);
+ Efficiency_Lab_2D->Divide(ThetaLab_emmitted_2D);
+
+ double dt = 180./Efficiency_Lab->GetNbinsX();
+ cout << "Angular infinitesimal = " << dt << "deg " << endl;
+ TF1* C = new TF1("C",Form("1./(2*%f*sin(x*%f/180.)*%f*%f/180.)",M_PI,M_PI,dt,M_PI),0,180);
+
+
+ SolidAngle_CM->Divide(ThetaCM_emmitted);
+ SolidAngle_CM->Divide(C,1);
+
+ SolidAngle_Lab->Divide(ThetaLab_emmitted);
+ SolidAngle_Lab->Divide(C,1);
+
+ SolidAngle_CM_2D->Divide(ThetaCM_emmitted_2D);
+ SolidAngle_CM_2D->Divide(C,1);
+
+ SolidAngle_Lab_2D->Divide(ThetaLab_emmitted_2D);
+ SolidAngle_Lab_2D->Divide(C,1);
+
+
+ HistoFile->Write();
+ HistoFile->Close();
+}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch(){
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaDetector",&ThetaDetector,"ThetaDetector/D");
+ RootOutput::getInstance()->GetTree()->Branch("X_Sharc",&X_Sharc,"X_Sharc/D");
+ RootOutput::getInstance()->GetTree()->Branch("Y_Sharc",&Y_Sharc,"Y_Sharc/D");
+ RootOutput::getInstance()->GetTree()->Branch("Z_Sharc",&Z_Sharc,"Z_Sharc/D");
+ RootOutput::getInstance()->GetTree()->Branch("X_Trifoil",&X_Trifoil,"X_Trifoil/D");
+ RootOutput::getInstance()->GetTree()->Branch("Y_Trifoil",&Y_Trifoil,"Y_Trifoil/D");
+
+
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ RootInput::getInstance()->GetChain()->SetBranchAddress("RunNumber",&RunNumber);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("RunNumberMinor",&RunNumberMinor);
+}
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+ Ex = -1000 ;
+ ELab = -1000;
+ ThetaLab = -1000;
+ ThetaCM = -1000;
+ ThetaDetector=-1000;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the AnalysisFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct(){
+ return (NPL::VAnalysis*) new Analysis();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+};
+
+proxy p;
+}
+
diff --git a/Projects/SharcEfficiency/Analysis.h b/Projects/SharcEfficiency/Analysis.h
new file mode 100644
index 000000000..c028e132c
--- /dev/null
+++ b/Projects/SharcEfficiency/Analysis.h
@@ -0,0 +1,117 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 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: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include"NPVAnalysis.h"
+#include"NPEnergyLoss.h"
+#include"NPReaction.h"
+#include"RootOutput.h"
+#include"RootInput.h"
+#include "TSharcPhysics.h"
+#include "TPlasticPhysics.h"
+#include "TInitialConditions.h"
+#include <TRandom3.h>
+#include <TVector3.h>
+#include <TMath.h>
+#include<TCutG.h>
+
+#include<fstream>
+#include<vector>
+class Analysis: public NPL::VAnalysis{
+ public:
+ Analysis();
+ ~Analysis();
+
+ public:
+ void Init();
+ void TreatEvent();
+ void End();
+
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+ private:
+ double Ex;
+ double ELab;
+ double ThetaLab;
+ double ThetaCM;
+ NPL::Reaction* myReaction;
+ TInitialConditions* myInit ;
+ // Energy loss table: the G4Table are generated by the simulation
+ EnergyLoss LightCD2;
+ EnergyLoss LightAl;
+ EnergyLoss LightSi;
+ EnergyLoss BeamCD2;
+ TVector3 BeamDirection;
+ TVector3 TargetPosition;
+
+
+ double TargetThickness ;
+ // Beam Energy
+ double OriginalBeamEnergy ; // AMEV
+ // intermediate variable
+ TRandom3 Rand ;
+ int DetectorNumber ;
+ int RunNumber;
+ int RunNumberMinor;
+ double ThetaNormalTarget;
+ double Energy ;
+
+ double ThetaSharcSurface ;
+ double X_Sharc ;
+ double Y_Sharc ;
+ double Z_Sharc ;
+ double X_Trifoil;
+ double Y_Trifoil;
+ double ThetaDetector;
+ double Si_E_Sharc ;
+ double E_Sharc ;
+ TSharcPhysics* Sharc;
+ TPlasticPhysics* Trifoil;
+
+ TFile* HistoFile;
+ TH1F* ThetaCM_emmitted;
+ TH1F* ThetaCM_detected;
+ TH1F* ThetaLab_emmitted;
+ TH1F* ThetaLab_detected;
+
+ TH1F* Efficiency_CM;
+ TH1F* Efficiency_Lab;
+
+ TH1F* SolidAngle_CM;
+ TH1F* SolidAngle_Lab;
+
+ TH2F* ThetaCM_emmitted_2D;
+ TH2F* ThetaCM_detected_2D;
+ TH2F* ThetaLab_emmitted_2D;
+ TH2F* ThetaLab_detected_2D;
+
+ TH2F* Efficiency_CM_2D;
+ TH2F* Efficiency_Lab_2D;
+
+ TH2F* SolidAngle_CM_2D;
+ TH2F* SolidAngle_Lab_2D;
+
+};
+#endif
diff --git a/Projects/SharcEfficiency/CMakeLists.txt b/Projects/SharcEfficiency/CMakeLists.txt
new file mode 100644
index 000000000..6806d778c
--- /dev/null
+++ b/Projects/SharcEfficiency/CMakeLists.txt
@@ -0,0 +1,2 @@
+cmake_minimum_required (VERSION 2.8)
+include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
diff --git a/Projects/SharcEfficiency/configs/ConfigSharc.dat b/Projects/SharcEfficiency/configs/ConfigSharc.dat
new file mode 100644
index 000000000..9c7ea5ce7
--- /dev/null
+++ b/Projects/SharcEfficiency/configs/ConfigSharc.dat
@@ -0,0 +1,99 @@
+ConfigSharc
+ STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA 20
+ STRIP_ENERGY_MATCHING_SIGMA 0.060
+% Detector 1
+ DISABLE_CHANNEL SH01STRF7
+ DISABLE_CHANNEL SH01STRF13
+% Detector 2
+ DISABLE_CHANNEL SH02STRF1
+ DISABLE_CHANNEL SH02STRB1
+% Detector 3
+ DISABLE_CHANNEL SH03STRB24
+% Detector 4
+ DISABLE_CHANNEL SH04STRF8
+ DISABLE_CHANNEL SH04STRF9
+ DISABLE_CHANNEL SH04STRF10
+ DISABLE_CHANNEL SH04STRF11
+ DISABLE_CHANNEL SH04STRB5
+ DISABLE_CHANNEL SH04STRB17
+ DISABLE_CHANNEL SH04STRB18
+ DISABLE_CHANNEL SH04STRB19
+ DISABLE_CHANNEL SH04STRB20
+ DISABLE_CHANNEL SH04STRB21
+% Detector 5
+ DISABLE_CHANNEL SH05STRF9
+ DISABLE_CHANNEL SH05STRF14
+ DISABLE_CHANNEL SH05STRF17
+ DISABLE_CHANNEL SH05STRB9
+ DISABLE_CHANNEL SH05STRB10
+ DISABLE_CHANNEL SH05STRB11
+ DISABLE_CHANNEL SH05STRB12
+ DISABLE_CHANNEL SH05STRB13
+ DISABLE_CHANNEL SH05STRB14
+ DISABLE_CHANNEL SH05STRB17
+ DISABLE_CHANNEL SH05STRB18
+ DISABLE_CHANNEL SH05STRB22
+ DISABLE_CHANNEL SH05STRB23
+ DISABLE_CHANNEL SH05STRB24
+ DISABLE_CHANNEL SH05STRB33
+% Detector 6
+ DISABLE_CHANNEL SH06STRF9
+ DISABLE_CHANNEL SH06STRF24
+ DISABLE_CHANNEL SH06STRB8
+ DISABLE_CHANNEL SH06STRB17
+ DISABLE_CHANNEL SH06STRB21
+ DISABLE_CHANNEL SH06STRB27
+% Detector 7
+ DISABLE_CHANNEL SH07STRB2
+ DISABLE_CHANNEL SH07STRB4
+ DISABLE_CHANNEL SH07STRB18
+ DISABLE_CHANNEL SH07STRB22
+ DISABLE_CHANNEL SH07STRB24
+ DISABLE_CHANNEL SH07STRB29
+ DISABLE_CHANNEL SH07STRB31
+ DISABLE_CHANNEL SH07STRB44
+ DISABLE_CHANNEL SH07STRB48
+% Detector 8
+ DISABLE_CHANNEL SH08STRF6
+ DISABLE_CHANNEL SH08STRF7
+ DISABLE_CHANNEL SH08STRF8
+ DISABLE_CHANNEL SH08STRF9
+ DISABLE_CHANNEL SH08STRF10
+ DISABLE_CHANNEL SH08STRF14
+ DISABLE_CHANNEL SH08STRF15
+ DISABLE_CHANNEL SH08STRB7
+ DISABLE_CHANNEL SH08STRB13
+ DISABLE_CHANNEL SH08STRB32
+ DISABLE_CHANNEL SH08STRB33
+ DISABLE_CHANNEL SH08STRB44
+ DISABLE_CHANNEL SH08STRB47
+ DISABLE_CHANNEL SH08STRB48
+% Detector 9
+ DISABLE_ALL SH09
+% Detector 10
+ DISABLE_CHANNEL SH10STRF7
+ DISABLE_CHANNEL SH10STRB1
+ DISABLE_CHANNEL SH10STRB2
+ DISABLE_CHANNEL SH10STRB3
+ DISABLE_CHANNEL SH10STRB4
+ DISABLE_CHANNEL SH10STRB8
+ DISABLE_CHANNEL SH10STRB17
+ DISABLE_CHANNEL SH10STRB36
+ % Detector 11
+ DISABLE_ALL SH11
+% Detector 12
+ DISABLE_CHANNEL SH12STRF4
+ DISABLE_CHANNEL SH12STRF7
+ DISABLE_CHANNEL SH12STRF24
+ DISABLE_CHANNEL SH12STRB1
+ DISABLE_CHANNEL SH12STRB17
+ DISABLE_CHANNEL SH12STRB18
+ DISABLE_CHANNEL SH12STRB28
+ DISABLE_CHANNEL SH12STRB29
+ DISABLE_CHANNEL SH12STRB36
+ DISABLE_CHANNEL SH12STRB45
+ DISABLE_CHANNEL SH12STRB46
+ DISABLE_CHANNEL SH12STRB47
+ DISABLE_CHANNEL SH12STRB48
+% PAD threshold
+ PAD_E_THRESHOLD 0.3
--
GitLab