Skip to content
Snippets Groups Projects
Commit 5e45de1a authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

Merge branch 'NPTool.2.dev' of https://gitlab.in2p3.fr/np/nptool into NPTool.2.dev

parents 7f9d5300 fca2f86f
No related branches found
No related tags found
No related merge requests found
Pipeline #165004 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Showing
with 283 additions and 198 deletions
NPLib/Physics/NPFissionDecay.cxx
+ 24
− 14
View file @ 5e45de1a
......@@ -125,6 +125,7 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
FFl.SetKineticEnergy(KEl);
FFh.SetKineticEnergy(KEh);
}
FissionFragments.push_back(FFl);
FissionFragments.push_back(FFh);
Ex.push_back(0);
......@@ -145,6 +146,11 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
double Phil = m_FissionModel->GetPhffl();
double Phih = m_FissionModel->GetPhffh();
if(Thetal != Thetal)
worked=false;
if(Thetah != Thetah)
worked=false;
TVector3 uxy = TVector3(cos(TMath::Pi()/2-PhiCN), sin(TMath::Pi()/2-PhiCN), 0);
TVector3 Momentuml = Pl * TVector3(sin(Thetal)*cos(Phil),
sin(Thetal)*sin(Phil),
......@@ -155,7 +161,7 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
sin(Thetah)*sin(Phih),
cos(Thetah));
Momentumh.Rotate(-ThetaCN, uxy);
DPx.push_back(Momentuml.X());
DPx.push_back(Momentumh.X());
DPy.push_back(Momentuml.Y());
......@@ -168,19 +174,23 @@ bool NPL::FissionDecay::GenerateEvent(string CompoundName, double MEx,double MEK
KE2 = m_FissionModel->GetKE2();
// Neutron and gamma emission
float* En1;
float* Eg1;
En1 = m_FissionModel->GetNeutronEnergyFrag1();
//cout << "----- Neutron energy: " << endl;
//for(int i=0; i<51; i++) {
// cout << "En= " << En1[i] << endl;
//}
Eg1 = m_FissionModel->GetGammaEnergyFrag1();
//cout << "----- Gamma energy: " << endl;
//for(int i=0; i<101; i++){
// cout << "Eg= " << Eg1[i] << endl;
//}
if(m_shoot_neutron){
float* En1;
En1 = m_FissionModel->GetNeutronEnergyFrag1();
//cout << "----- Neutron energy: " << endl;
//for(int i=0; i<51; i++) {
// cout << "En= " << En1[i] << endl;
//}
}
if(m_shoot_gamma){
float* Eg1;
Eg1 = m_FissionModel->GetGammaEnergyFrag1();
//cout << "----- Gamma energy: " << endl;
//for(int i=0; i<101; i++){
// cout << "Eg= " << Eg1[i] << endl;
//}
}
}
}
return worked;
......
NPSimulation/Process/FissionDecay.cc
+ 5
− 8
View file @ 5e45de1a
......@@ -107,13 +107,13 @@ G4bool FissionDecay::ModelTrigger(const G4FastTrack& fastTrack) {
const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
int Parent_ID = PrimaryTrack->GetParentID();
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
m_FissionConditions->Clear();
if(Parent_ID>=0){
Trigger = true;
m_PreviousEnergy=fastTrack.GetPrimaryTrack()->GetKineticEnergy();
}
return Trigger;
}
......@@ -150,11 +150,11 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
std::vector<double> DPz;
double TKE, KE1, KE2;
G4bool IsFissionDecay = m_FissionDecay.GenerateEvent(NPL::ChangeNameFromG4Standard(m_CurrentName),m_ExcitationEnergy,energy,
pdirection.x(),pdirection.y(),pdirection.z(),
FissionFragment, Ex,DEK,DPx,DPy,DPz,
TKE, KE1, KE2);
if(IsFissionDecay){
/////////////////////////////////////////////////
// Fillion the attached Fission condition Tree //
......@@ -187,7 +187,7 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
int FFZ = FissionFragment[i].GetZ();
int FFA = FissionFragment[i].GetA();
FissionFragmentDef=NULL;
// Set the momentum direction
G4ThreeVector Momentum (DPx[i],DPy[i],DPz[i]);
Momentum=Momentum.unit();
......@@ -197,7 +197,6 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
m_FissionConditions->SetFragmentZ(FFZ);
m_FissionConditions->SetFragmentA(FFA);
//m_FissionConditions->SetFragmentKineticEnergy(DEK[i]);
m_FissionConditions->SetFragmentKineticEnergy(KineticEnergy);
m_FissionConditions->SetFragmentBrho(Brho);
m_FissionConditions->SetFragmentTheta(Momentum.theta()/deg);
......@@ -222,13 +221,11 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
// the rest
else
FissionFragmentDef=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(FFZ, FFA, Ex[i]);
G4DynamicParticle DynamicFissionFragment(FissionFragmentDef,Momentum,DEK[i]);
fastStep.CreateSecondaryTrack(DynamicFissionFragment, localPosition, time);
}
if(size){
if(size>0){
// Set the end of the step conditions
fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0,pdirection);
fastStep.SetPrimaryTrackFinalPosition(worldPosition);
......
Projects/Sofia/sofia.detector
+ 1
− 1
View file @ 5e45de1a
......@@ -15,7 +15,7 @@ SofTofW
Build_GLAD= 0
GLAD_TiltAngle= 14 deg
GLAD_DistanceFromTarget= 3.05 m
Build_MagneticField= 0
Build_MagneticField= 1
GLAD_MagField= 1.70 T
Build_VacuumPipe= 1
VacuumPipeX= 0 cm
......
Projects/e793s/macro/CS2.h
+ 95
− 9
View file @ 5e45de1a
/* Predefine functions */
vector<vector<double>> GetExpDiffCross(double Energy);
TH1F* PullThetaLabHist(int i, double minTheta, double gatesize);
TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize);
void Scale(TGraph* g , TGraphErrors* ex);
TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j);
double ToMininize(const double* parameter);
......@@ -13,7 +14,10 @@ TGraphErrors* staticExp;
int indexE;
/* Output volume toggle */
bool loud = 1;
bool loud = 0;
/* Scale method toggle */
bool scaleTogether = 1;
////////////////////////////////////////////////////////////////////////////////
void CS(){
......@@ -113,6 +117,14 @@ void CS(double Energy, double Spin, double spdf, double angmom){
anglecentres.push_back(((areaArray[i][4]-areaArray[i][3])/2.)+areaArray[i][3]);
anglewidth.push_back(areaArray[i][4]-areaArray[i][3]);
double tempsum=0, tempsumerr=0;
for(int x=binmin;x<binmax+1;x++){
tempsum += SolidAngle->GetBinContent(x);
tempsumerr += SolidAngle->GetBinError(x);
if(loud){cout << x << endl;}
}
if(loud){cout << "TEST CHECK " << tempsum << " +- " << tempsumerr << endl;}
double SAerr;
double SA = SolidAngle->IntegralAndError(binmin,binmax,SAerr);
SA = SA*1000.; //sr->msr
......@@ -162,20 +174,20 @@ void CS(double Energy, double Spin, double spdf, double angmom){
0, &(expDCSerr[0]) );
gdSdO->SetTitle("Differential Cross Section");
gdSdO->GetXaxis()->SetTitle("ThetaLab [deg]");
gdSdO->GetYaxis()->SetTitle("d#sigma/d#Theta [mb/msr]");
gdSdO->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
gdSdO->Draw();
/* TWOFNR diff. cross section, in mb/msr */
TCanvas* c_TWOFNR = new TCanvas("c_TWOFNR","c_TWOFNR",1000,1000);
c_TWOFNR->SetLogy();
TGraph* TheoryDiffCross = TWOFNR(means[indexE], J0, Spin, nodes, spdf, angmom);
TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Theta [mb/msr]"); //msr set in func above
TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); //msr set in func above
TheoryDiffCross->GetXaxis()->SetTitle("ThetaLab [deg]");
TheoryDiffCross->Draw();
/* Convert AoSA into Differential Cross Section */
//cout << " SCALING BY NORMALIZATION = " << ElasticNorm << endl;
//gAoSA->GetYaxis()->SetTitle("d#sigma/d#Theta [mb/msr]");
//gAoSA->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
/* Scaled and compared on same plot */
......@@ -211,6 +223,32 @@ vector<vector<double>> GetExpDiffCross(double Energy){
vector<vector<double>> OnePeak_AllGates;
int numbins = 10;
double x[numbins], y[numbins];
TList* list = new TList();
/* Determine scaling factor for PhaseSpace */
double trackScale = 0.0;
if(scaleTogether){
TH1F* baseEx = PullThetaLabHist(0,105.,5.);
TH1F* basePS = PullPhaseSpaceHist(0,105.,5.);
for(int i=1; i<numbins;i++){
TH1F* addEx = PullThetaLabHist(i,105.,5.); baseEx->Add(addEx,1.);
TH1F* addPS = PullPhaseSpaceHist(i,105.,5.); basePS->Add(addPS,1.);
}
basePS->Scale(0.01);
trackScale = 0.01;
int numbinsScale = baseEx->GetNbinsX();
for(int b=0; b<numbinsScale; b++){
while(basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
basePS->Scale(0.99999);
trackScale *= 0.99999;
}
}
baseEx->Add(basePS,-1.);
baseEx->SetName("AllAngles");
list->Add(baseEx);
cout << " !!!!!!!!!!!!!!!FINAL SCALING = " << trackScale << endl;
}
for(int i=0; i<numbins;i++){
double bin = 5.;
double min = 105. + (i*bin);
......@@ -221,12 +259,41 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
/* To change angle gates, run GateThetaLab_MultiWrite() */
TH1F* gate = PullThetaLabHist(i,105.,5.);
TH1F* pspace = PullPhaseSpaceHist(i,105.,5.);
/* Scale Phase Space... */
/* ... for all angles together */
if(scaleTogether){
gate->Add(pspace,-trackScale);
}
/* ... or seperately for each angular bin */
else {
if(pspace->Integral() > 50.){ // Non-garbage histogram
pspace->Scale(0.01);
int numbins = gate->GetNbinsX();
for(int b=0; b<numbins; b++){
if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
" > " << gate->GetBinContent(b);
}
while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
pspace->Scale(0.9999);
}
if(loud){cout << " TO " << pspace->GetBinContent(b) <<
" > " << gate->GetBinContent(b) << endl;
}
}
gate->Add(pspace,-1);
}
}
/* Retrieve array containing all fits, for one angle gate. *
* Specific peak of interest selected from the vector by *
* global variable indexE */
AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate);
/* Write PS-subtracted spectrum to list */
list->Add(gate);
/* Check correct OneGate vector is selected */
cout << "area of " << means[indexE] << " = "
<< AllPeaks_OneGate[indexE][1]
......@@ -240,17 +307,37 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* Push relevant OneGate vector to end of AllGates vector */
OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
}
/* Write PS-subtracted spectrum to file */
TFile* outfile = new TFile("GateThetaLab_ExMinusGatePhaseSpace.root","RECREATE");
list->Write("GateThetaLab_ExMinusPhaseSpace",TObject::kSingleKey);
return OnePeak_AllGates;
}
////////////////////////////////////////////////////////////////////////////////
TH1F* PullThetaLabHist(int i, double minTheta, double gatesize){
TFile* file = new TFile("GateThetaLabHistograms.root","READ");
TFile* file = new TFile("GateThetaLabHistograms_ReadMe.root","READ");
string histname = "cThetaLabGate_"
+ to_string(minTheta+(i*gatesize)) + "-"
+ to_string(minTheta+((i+1)*gatesize));
+ to_string((int) (minTheta+(i*gatesize))) + "-"
+ to_string((int) (minTheta+((i+1)*gatesize)));
cout << "Loading " << histname << endl;
TList *list = (TList*)file->Get("GateThetaLabHistograms");
TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
// file->Close();
return hist;
}
////////////////////////////////////////////////////////////////////////////////
TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize){
TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_ReadMe.root","READ");
string histname = "cPSpaceThetaLabGate_"
+ to_string((int) (minTheta+(i*gatesize))) + "-"
+ to_string((int) (minTheta+((i+1)*gatesize)));
cout << "Loading " << histname << endl;
TList *list = (TList*)file->Get("GatePhaseSpaceThetaLabHistograms");
TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
file->Close();
return hist;
}
......@@ -389,7 +476,7 @@ TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
}
////////////////////////////////////////////////////////////////////////////////
double Chi2(TGraph* theory , TGraphErrors* exper){
double Chi2(TGraph* theory, TGraphErrors* exper){
double Chi2 = 0 ;
double chi;
//for(int i = 1 ; i < exper->GetN() ; i++){
......@@ -403,7 +490,6 @@ double Chi2(TGraph* theory , TGraphErrors* exper){
return Chi2;
}
////////////////////////////////////////////////////////////////////////////////
double ToMininize(const double* parameter){
static int f = 0 ;
......
Projects/e793s/macro/DrawPlots.C
+ 3
− 0
View file @ 5e45de1a
......@@ -4,6 +4,7 @@
#include "CS2.h"
#include "ThreeBodyBreakup.h"
#include "ThreeBodyBreakup_FitPhaseSpace.h"
/* USE THIS SPACE TO TEST NEW FEATURES */
void thickness(){
......@@ -523,10 +524,12 @@ void DrawPlots(){
cout << "========== AVAILABLE FUNCTIONS ===========" << endl;
cout << " 2D Matrices " << endl;
cout << "\t- Draw_2DParticleGamma() "<< endl;
cout << "\t- Load_2DParticleGamma() "<< endl;
cout << "\t- Draw_2DGammaGamma() "<< endl;
cout << ""<< endl;
cout << " Ungated histograms " << endl;
cout << "\t- Draw_1DParticle() "<< endl;
cout << "\t- Load_1DParticle() "<< endl;
cout << "\t- Draw_1DParticle_MUST2() "<< endl;
cout << "\t- Draw_1DGamma() "<< endl;
cout << "\t- Draw_1DGamma_MG() "<< endl;
......
Projects/e793s/macro/DrawPlots.h
+ 72
− 30
View file @ 5e45de1a
......@@ -40,22 +40,13 @@ TChain* Chain(std::string TreeName, std::vector<std::string>& file, bool EventLi
void LoadChainNP(){
vector<string> files;
//files.push_back("../../../Outputs/Analysis/OriginalValues_ptI.root");
//files.push_back("../../../Outputs/Analysis/OriginalValues_ptII.root");
//files.push_back("../../../Outputs/Analysis/OriginalValues_ptIII.root");
//files.push_back("../../../Outputs/Analysis/47K_Full_07Sep_PartI.root");
//files.push_back("../../../Outputs/Analysis/47K_Full_07Sep_PartII.root");
//files.push_back("../../../Outputs/Analysis/47Kdp_11Oct_PartI.root");
//files.push_back("../../../Outputs/Analysis/47Kdp_11Oct_PartII.root");
//files.push_back("../../../Outputs/Analysis/47Kdp_13Oct_wildTry_PartI.root");
//files.push_back("../../../Outputs/Analysis/47Kdp_13Oct_wildTry_PartII.root");
files.push_back("../../../Outputs/Analysis/47Kdp_08Nov_PartI.root");
files.push_back("../../../Outputs/Analysis/47Kdp_08Nov_PartII.root");
// files.push_back("../../../Outputs/Analysis/47Kdp_17Feb_AGATA_RotateBYpi.root");
// files.push_back("../../../Outputs/Analysis/47Kdp_19Feb_AGATA_RotateBXYZ.root");
chain = Chain("PhysicsTree",files,true);
}
......@@ -183,9 +174,16 @@ void Draw_1DGamma_MM(){
Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
}
void Load_1DParticle(){
TH1F *hEx = new TH1F("hEx","Loaded 1D Particle Spectrum",200,-1,9);
TFile *file = new TFile("LoadHistograms/Load_1DParticle.root","READ");
hEx = (TH1F*)file->Get("Ep");
hEx->Draw();
}
void Draw_1DParticle(){
TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000);
chain->Draw("Ex>>Ep(160,-1,7)",
chain->Draw("Ex>>Ep(200,-1,9)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"");
TH1F* Ep = (TH1F*) gDirectory->Get("Ep");
......@@ -198,7 +196,7 @@ void Draw_1DParticle(){
void Draw_1DParticleMUST2(){
TCanvas *cEx = new TCanvas("cEx","cEx",1000,1000);
chain->Draw("Ex>>Ep(160,-1,7)",
chain->Draw("Ex>>Ep(200,-1,9)",
"abs(T_MUGAST_VAMOS-2777)<600 && MUST2.TelescopeNumber>0 && MUST2.TelescopeNumber<4",
"");
TH1F* Ep = (TH1F*) gDirectory->Get("Ep");
......@@ -207,9 +205,16 @@ void Draw_1DParticleMUST2(){
Ep->GetYaxis()->SetTitle("Counts / 0.05 MeV");
}
void Load_2DParticleGamma(){
TH2F *hExEg = new TH2F("hExEg","Loaded 2D Particle-Gamma",200,-1,9,2500,0,5);
TFile *file = new TFile("LoadHistograms/Load_2DParticleGamma.root","READ");
hExEg = (TH2F*)file->Get("ExEg");
hExEg->Draw();
}
void Draw_2DParticleGamma(){
TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000);
chain->Draw("AddBack_EDC:Ex>>ExEg(140,-1,7,2500,0,5)",
chain->Draw("AddBack_EDC:Ex>>ExEg(200,-1,9,2500,0,5)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"colz");
TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg");
......@@ -287,7 +292,7 @@ cout << " NO MG3 IN THIS ONE!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
string title = to_string(gamma-width)+" < Eg < "+to_string(gamma+width);
string ytitle = "Counts / " + to_string(binsize) + " MeV";
string draw = "Ex>>ExGate(" + to_string(8.0/binsize) + ",-1,7)";
string draw = "Ex>>ExGate(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000);
//chain->Draw("Ex>>ExGate(60,-1,5)",gating.c_str(),"colz");
......@@ -314,7 +319,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){
+ " BG: "+to_string(bg-width)+" to "+to_string(bg+width)+".";
TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000);
chain->Draw("Ex>>ExGate(80,-1,7)",gating.c_str(),"");
chain->Draw("Ex>>ExGate(100,-1,9)",gating.c_str(),"");
//chain->Draw("Ex>>ExGate(120,-1,5)",gating.c_str(),"");
TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate");
ExGate->GetXaxis()->SetTitle("Ex [MeV]");
......@@ -325,7 +330,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg){
ExGate->SetFillStyle(3154);
ExGate->SetTitle(title.c_str());
chain->Draw("Ex>>ExBG(80,-1,7)",bggate.c_str(),"same");
chain->Draw("Ex>>ExBG(100,-1,9)",bggate.c_str(),"same");
//chain->Draw("Ex>>ExBG(120,-1,5)",bggate.c_str(),"same");
TH1F* ExBG = (TH1F*) gDirectory->Get("ExBG");
ExBG->SetLineColor(kRed);
......@@ -351,7 +356,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg, double
+ " BG: "+to_string(bg-width)+" to "+to_string(bg+width)+".";
TCanvas *cEx_Gate = new TCanvas("cEx_Gate","cEx_Gate",1000,1000);
chain->Draw("Ex>>ExGate(80,-1,7)",gating.c_str(),"");
chain->Draw("Ex>>ExGate(100,-1,9)",gating.c_str(),"");
//chain->Draw("Ex>>ExGate(120,-1,5)",gating.c_str(),"");
TH1F* ExGate = (TH1F*) gDirectory->Get("ExGate");
ExGate->GetXaxis()->SetTitle("Ex [MeV]");
......@@ -362,7 +367,7 @@ void GateGamma_SeeParticle_WithBG(double gamma, double width, double bg, double
ExGate->SetFillStyle(3154);
ExGate->SetTitle(title.c_str());
chain->Draw("Ex>>ExBG(80,-1,7)",bggate.c_str(),"same");
chain->Draw("Ex>>ExBG(100,-1,9)",bggate.c_str(),"same");
//chain->Draw("Ex>>ExBG(120,-1,5)",bggate.c_str(),"same");
TH1F* ExBG = (TH1F*) gDirectory->Get("ExBG");
ExBG->Scale(ratio);
......@@ -621,7 +626,7 @@ void CompareSimExp(){
TCanvas *cSimExp = new TCanvas("cSimExp","cSimExp",1000,1000);
gStyle->SetOptStat(0);
chain->Draw("Ex>>hexp(70,-1,6)","abs(T_MUGAST_VAMOS-2777)<600","");
chain->Draw("Ex>>hexp(100,-1,9)","abs(T_MUGAST_VAMOS-2777)<600","");
TH1F* hexp = (TH1F*) gDirectory->Get("hexp");
hexp->SetTitle("Comparing Simulation to Experiment");
hexp->GetXaxis()->SetTitle("Ex [MeV]");
......@@ -632,7 +637,7 @@ void CompareSimExp(){
//TFile* simfile = new TFile("../../../Outputs/Analysis/SimTest_Jun22_TWOFNR_p32.root",
"READ");
TTree* simtree = (TTree*) simfile->FindObjectAny("PhysicsTree");
simtree->Draw("Ex>>hsimMGp32(70,-1,6)",
simtree->Draw("Ex>>hsimMGp32(100,-1,9)",
"Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8","same");
TH1F* hsimMGp32 = (TH1F*) gDirectory->Get("hsimMGp32");
hsimMGp32->SetLineColor(kBlue);
......@@ -890,7 +895,7 @@ void ExPhiLab_ForPoster(){
void ExThetaLab(){
TCanvas *diagnoseTheta = new TCanvas("diagnoseTheta","diagnoseTheta",1000,1000);
chain->Draw(
"Ex:ThetaLab>>thetaHist(60,100,160,120,-1,5)",
"Ex:ThetaLab>>thetaHist(60,100,160,100,-1,9)",
"abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8",
"colz");
TH1F* thetaHist = (TH1F*) gDirectory->Get("thetaHist");
......@@ -939,7 +944,7 @@ void ExThetaLab(double gamma, double width){
string gating = "abs(T_MUGAST_VAMOS-2777)<600 && Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && abs(AddBack_EDC-"
+ to_string(gamma) + ") < " + to_string(width);
chain->Draw("Ex:ThetaLab>>thetaHist(60,100,160,60,-1,5)", gating.c_str(), "colz");
chain->Draw("Ex:ThetaLab>>thetaHist(60,100,160,100,-1,9)", gating.c_str(), "colz");
TH1F* thetaHist = (TH1F*) gDirectory->Get("thetaHist");
thetaHist->GetXaxis()->SetTitle("Theta (degrees)");
thetaHist->GetYaxis()->SetTitle("Ex [MeV]");
......@@ -1181,7 +1186,7 @@ void GateThetaCM(double minTheta, double maxTheta, double binsize){
string title = to_string(minTheta)+" < ThetaCM < "+to_string(maxTheta);
string ytitle = "Counts / " + to_string(binsize) + " MeV";
string draw = "Ex>>Ex_ThetaCMGate(" + to_string(8.0/binsize) + ",-1,7)";
string draw = "Ex>>Ex_ThetaCMGate(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaCMGate = new TCanvas("cEx_ThetaCMGate","cEx_ThetaCMGate",1000,1000);
chain->Draw(draw.c_str(),gating.c_str(),"colz");
......@@ -1201,7 +1206,7 @@ void GateThetaLab(double minTheta, double maxTheta, double binsize){
string title = to_string(minTheta)+" < ThetaLab < "+to_string(maxTheta);
string ytitle = "Counts / " + to_string(binsize) + " MeV";
string draw = "Ex>>Ex_ThetaLabGate(" + to_string(8.0/binsize) + ",-1,7)";
string draw = "Ex>>Ex_ThetaLabGate(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaLabGate = new TCanvas("cEx_ThetaLabGate","cEx_ThetaLabGate",1000,1000);
chain->Draw(draw.c_str(),gating.c_str(),"colz");
......@@ -1221,19 +1226,20 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
for (int i=0; i<numGates; i++){
double minTheta = startTheta + (i * gatesize);
string title = to_string(minTheta)+" < ThetaLab < "+to_string(minTheta+gatesize);
string title = to_string((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize));
string gating = core
+ to_string(minTheta)
+ " && ThetaLab < "
+ to_string(minTheta+gatesize);
string histname = "cThetaLabGate_" + to_string(minTheta) + "-" + to_string(minTheta+gatesize);
string draw = "Ex>>" + histname + "(" + to_string(8.0/binsize) + ",-1,7)";
string histname = "cThetaLabGate_" + to_string((int) minTheta) + "-" + to_string((int) (minTheta+gatesize));
string draw = "Ex>>" + histname + "(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cEx_ThetaLabGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000);
chain->Draw(draw.c_str(),gating.c_str(),"colz");
TH1F* Ex_ThetaLabGate = (TH1F*) gDirectory->Get(histname.c_str());
Ex_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
Ex_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
Ex_ThetaLabGate->Sumw2();
Ex_ThetaLabGate->SetTitle(title.c_str());
list->Add(Ex_ThetaLabGate);
delete cEx_ThetaLabGate;
......@@ -1243,3 +1249,39 @@ void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates
list->Write("GateThetaLabHistograms",TObject::kSingleKey);
file->ls();
}
void GatePhaseSpaceByThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize){
string core = "EventWeight*(Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8 && ThetaLab > ";
string ytitle = "Counts / " + to_string(binsize) + " MeV";
double gatesize = (finishTheta-startTheta)/numGates;
TList* list = new TList();
TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_02Mar_47Kdp_PhaseSpace.root","READ");
TTree* PSTree = (TTree*) psfile->FindObjectAny("PhysicsTree");
for (int i=0; i<numGates; i++){
double minTheta = startTheta + (i * gatesize);
string title = to_string((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize));
string gating = core
+ to_string(minTheta)
+ " && ThetaLab < "
+ to_string(minTheta+gatesize)
+ ")";
string histname = "cPSpaceThetaLabGate_" + to_string((int) minTheta) + "-" + to_string((int) (minTheta+gatesize));
string draw = "Ex>>" + histname + "(" + to_string(10.0/binsize) + ",-1,9)";
TCanvas *cPSpace_ThetaLabGate = new TCanvas(histname.c_str(),histname.c_str(),1000,1000);
PSTree->Draw(draw.c_str(),gating.c_str(),"colz");
TH1F* PSpace_ThetaLabGate = (TH1F*) gDirectory->Get(histname.c_str());
PSpace_ThetaLabGate->GetXaxis()->SetTitle("Ex [MeV]");
PSpace_ThetaLabGate->GetYaxis()->SetTitle(ytitle.c_str());
PSpace_ThetaLabGate->Sumw2();
PSpace_ThetaLabGate->SetTitle(title.c_str());
list->Add(PSpace_ThetaLabGate);
delete cPSpace_ThetaLabGate;
}
TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms.root","RECREATE");
list->Write("GatePhaseSpaceThetaLabHistograms",TObject::kSingleKey);
file->ls();
}
Projects/e793s/macro/KnownPeakFitter.h
+ 76
− 134
View file @ 5e45de1a
......@@ -3,23 +3,24 @@
#include <cmath>
#include "stdlib.h"
const int numPeaks = 14;
array<double,14> means = { 0.000,
0.143,
0.279,
0.728,
0.968,
1.410,
1.981,
2.410,
2.910,
3.2,
3.605,
3.792,
4.1,
4.4
const int numPeaks = 15;
array<double,numPeaks> means = { 0.000,
0.143,
0.279,
0.728,
0.968,
1.410,
1.981,
2.410,
2.910,
3.2,
3.605,
3.792,
4.1,
4.4,
5.24
};
/*
Double_t f_bg(Double_t *x, Double_t *par){
// Flat bg [0] + semicircle [1]*sqrt(6.183^2 - (x-10.829)^2)
Float_t xx = x[0];
......@@ -32,72 +33,58 @@ Double_t f_bg(Double_t *x, Double_t *par){
}
Double_t f_peak(Double_t *x, Double_t *par){
Float_t xx = x[0];
Double_t f = (par[2]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[1])/par[0],2));
float xx = x[0];
double f = (par[2]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[1])/par[0],2));
return f;
}
Double_t f_full(Double_t *x, Double_t *par){
Float_t xx = x[0];
Double_t f;
Double_t peaks = (par[2]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[1])/par[0],2))
+ (par[4]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[3])/par[0],2))
+ (par[6]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[5])/par[0],2))
+ (par[8]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[7])/par[0],2))
+ (par[10]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[9])/par[0],2))
+ (par[12]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[11])/par[0],2))
+ (par[14]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[13])/par[0],2))
+ (par[16]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[15])/par[0],2))
+ (par[18]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[17])/par[0],2))
+ (par[20]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[19])/par[0],2))
+ (par[22]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[21])/par[0],2))
+ (par[24]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[23])/par[0],2))
+ (par[26]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[25])/par[0],2))
+ (par[28]/(par[0]*sqrt(2*pi)))*exp(-0.5*pow((xx-par[27])/par[0],2));
Double_t bg;
Double_t a = TMath::Power(6.183,2);
Double_t b = TMath::Power(xx-10.829,2);
if(a > b){ bg = par[29] + (par[30]*TMath::Sqrt(a-b)); }
else{ bg = par[29]; }
f = peaks + bg;
return f;
*/
Double_t f_full(Double_t *x, Double_t *par) {
float xx = x[0];
double result, norm;
// Flat background
result = par[0];
// Add N peaks
for(int pk=0; pk<numPeaks; pk++){
result += (par[3+(pk*3)]/(par[1+(pk*3)]*sqrt(2*pi)))
* exp(-0.5*pow((xx-par[2+(pk*3)])/par[1+(pk*3)],2));
}
return result;
}
vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
double minFit=-1.0, maxFit=5.0;
double minFit=-1.0, maxFit=8.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
double sigma = 0.14;
string nameBase = "Peak ";
string function = "([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2))";
//Assign peak centroids
// const int numPeaks = 14;
// array<double,14> means = { 0.000,
// 0.143,
// 0.279,
// 0.728,
// 0.968,
// 1.410,
// 1.981,
// 2.410,
// 2.910,
// 3.2,
// 3.605,
// 3.792,
// 4.1,
// 4.4
// };
hist->Sumw2();
/* Construct flat BG to subtract */
/**
cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
cout << " REMOVING FLAT BG OF 36 COUNTS!!!!" << endl;
double ConstBG = 36.0; double ErrBG = 1.0;
int xbins = hist->GetXaxis()->GetNbins();
double xmin = hist->GetXaxis()->GetXmin();
double xmax = hist->GetXaxis()->GetXmax();
TH1F *FlatBG = new TH1F("FlatBG","FlatBG", xbins, xmin, xmax);
for(int i=0; i<xbins;i++){
FlatBG->SetBinContent(i,ConstBG);
FlatBG->SetBinError(i,ErrBG);
}
hist->Add(FlatBG,-1);
**/
//Build individual peak fit functions
string nameBase = "Peak ";
string function = "([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2))";
TF1 **allPeaks = new TF1*[numPeaks];
for(int i=0; i<numPeaks; i++) {
string nameHere = nameBase;
nameHere +=to_string(i);
allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), -1, 5);
allPeaks[i] = new TF1(nameHere.c_str(), function.c_str(), minFit, maxFit);
//allPeaks[i] = new TF1(nameHere.c_str(), f_peak, -1, 5);
allPeaks[i]->SetLineColor(kBlack);
allPeaks[i]->SetLineStyle(7);
......@@ -110,71 +97,24 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
bg->SetLineStyle(9);
bg->SetParNames("Background");
//Build IMPROVED background function
//TF1 *bg = new TF1 ("bg",f_bg, minFit, maxFit);
//bg->SetLineColor(kGreen);
//bg->SetLineStyle(9);
//bg->SetParNames("Background","BreakupScale");
//Build total function
TF1 *full = new TF1("fitAllPeaks",
" ([2]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[1])/[0],2)) "
"+ ([4]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[3])/[0],2)) "
"+ ([6]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[5])/[0],2)) "
"+ ([8]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[7])/[0],2)) "
"+ ([10]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[9])/[0],2)) "
"+ ([12]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[11])/[0],2)) "
"+ ([14]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[13])/[0],2)) "
"+ ([16]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[15])/[0],2)) "
"+ ([18]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[17])/[0],2)) "
"+ ([20]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[19])/[0],2)) "
"+ ([22]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[21])/[0],2)) "
"+ ([24]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[23])/[0],2)) "
"+ ([26]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[25])/[0],2)) "
"+ ([28]/([0]*sqrt(2*pi)))*exp(-0.5*pow((x-[27])/[0],2)) "
"+ [29]" , minFit, maxFit);
full->SetLineColor(kRed);
//Build IMPROVED total function
//TF1 *full = new TF1("fitAllPeaks", f_full, minFit, maxFit);
//full->SetLineColor(kRed);
//Annoyingly long parameter name assignment
//(SetParNames only works for up to 11 variables)
full->SetParName(0,"Sigma");
full->SetParName(1,"Mean 01"); full->SetParName(2,"Area*BinWidth 01");
full->SetParName(3,"Mean 02"); full->SetParName(4,"Area*BinWidth 02");
full->SetParName(5,"Mean 03"); full->SetParName(6,"Area*BinWidth 03");
full->SetParName(7,"Mean 04"); full->SetParName(8,"Area*BinWidth 04");
full->SetParName(9,"Mean 05"); full->SetParName(10,"Area*BinWidth 05");
full->SetParName(11,"Mean 06"); full->SetParName(12,"Area*BinWidth 06");
full->SetParName(13,"Mean 07"); full->SetParName(14,"Area*BinWidth 07");
full->SetParName(15,"Mean 08"); full->SetParName(16,"Area*BinWidth 08");
full->SetParName(17,"Mean 09"); full->SetParName(18,"Area*BinWidth 09");
full->SetParName(19,"Mean 10"); full->SetParName(20,"Area*BinWidth 10");
full->SetParName(21,"Mean 11"); full->SetParName(22,"Area*BinWidth 11");
full->SetParName(23,"Mean 12"); full->SetParName(24,"Area*BinWidth 12");
full->SetParName(25,"Mean 13"); full->SetParName(26,"Area*BinWidth 13");
full->SetParName(27,"Mean 14"); full->SetParName(28,"Area*BinWidth 14");
full->SetParName(29,"Background");
full->SetParName(30,"BreakupScale");
//Fix sigma & centroid, only allow area to vary
const int numParams = (numPeaks*2)+2;
full->FixParameter(0, sigma);
TF1 *full = new TF1("fitAllPeaks", f_full, minFit, maxFit, (int) 1+(3*numPeaks));
full->SetLineColor(kRed);
const int numParams = (numPeaks*3)+1;
for(int i=0; i<numPeaks; i++) {
full->FixParameter(2*i+1,means.at(i));
full->SetParameter(2*i+2,1.0);
full->SetParLimits(2*i+2,0.0,1e5);
full->FixParameter((i*3)+1,sigma);
full->FixParameter((i*3)+2,means.at(i));
full->SetParameter((i*3)+3,1e1);
full->SetParLimits((i*3)+3,0.0,1e5);
}
//full->SetParameter(0,30.);
//full->SetParLimits(0,0.,40.);
full->FixParameter(0,0.);
//TESTING!!!!!
//full->FixParameter(6,0.);
//full->FixParameter(numParams-1,0.0);
full->SetParameter(numParams-1,1.0);
full->SetParLimits(numParams-1,0.0,1e1);
// Specific limits
// Set max 279 counts to 100
full->SetParLimits(9,0.0,1e2); // Doesnt work???
allPeaks[14]->SetLineColor(kOrange);
//Fit full function to histogram
hist->Fit(full, "RQ", "", minFit, maxFit);
......@@ -184,9 +124,9 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
const Double_t* finalPar = full->GetParameters();
const Double_t* finalErr = full->GetParErrors();
for (int i=0; i<numPeaks; i++){
allPeaks[i]->SetParameters(sigma, means.at(i), finalPar[2*i+2]);
allPeaks[i]->SetParameters(sigma, means.at(i), finalPar[3+(i*3)]);
}
bg->SetParameter(0,finalPar[numParams-1]);
bg->SetParameter(0,finalPar[0]);
bg->Draw("SAME");
full->Draw("SAME");
......@@ -207,22 +147,24 @@ vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist){
vector<vector<double>> allpeaks;
for(int i=0; i<numPeaks; i++){
double A = finalPar[2*i+2]/binWidth;
double deltaA = A * (finalErr[2*i+2]/finalPar[2*i+2]);
double A = finalPar[(i*3)+3]/binWidth;
double deltaA = A * (finalErr[(i*3)+3]/finalPar[(i*3)+3]);
cout << fixed << setprecision(3)
<< " #" << i << " "
<< finalPar[2*i+1] << "\t" << setprecision(0)
<< finalPar[(i*3)+2] << "\t" << setprecision(0)
<< A << "\t+- "
<< deltaA << setprecision(3)
<< endl;
vector<double> onepeak; //energy, area and error for one peak
onepeak.push_back(finalPar[2*i+1]);
onepeak.push_back(finalPar[(i*3)+2]);
onepeak.push_back(A);
onepeak.push_back(deltaA);
allpeaks.push_back(onepeak);
}
cout << " BG " << full->GetParameter(0)
<< " +- " << full->GetParError(0) << endl;
return allpeaks;
}
......
Projects/e793s/macro/ThreeBodyBreakup.h
+ 7
− 2
View file @ 5e45de1a
......@@ -8,8 +8,12 @@ Double_t f_semi(Double_t *x, Double_t *par){
return f;
}
void ThreeBodyBreakup(){
TFile *fweight = new TFile("../../../Outputs/Simulation/EventWeight_25Feb.root");
TCanvas *cweight = (TCanvas*)fweight->Get("c1");
TH1F *hweight = (TH1F*)cweight->GetPrimitive("hEventWeight");
// hweight->Draw();
TGenPhaseSpace PS1p;
//47K(d,p)47K+p+n
......@@ -42,10 +46,11 @@ void ThreeBodyBreakup(){
double ex=R.ReconstructRelativistic(ELab,theta);
//cout << ELab << " " << ex << endl;
Ex1p->Fill(ex,weight);
}
}
TCanvas* c_Ex1p = new TCanvas("c_Ex1p","c_Ex1p",1000,500);
Ex1p->Draw();
Ex1p->Multiply(hweight);
// TCanvas* c_semi = new TCanvas("c_semi","c_semi",1000,500);
TF1 *f1 = new TF1("f_semi",f_semi,-10,40,3);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment