diff --git a/NPLib/Physics/NPFissionDecay.cxx b/NPLib/Physics/NPFissionDecay.cxx
index ff47f954de0c2970f9ee2ff1255f766c18530853..9083a63ac3f95d90fee31e8387e50a1b3c572da7 100644
--- a/NPLib/Physics/NPFissionDecay.cxx
+++ b/NPLib/Physics/NPFissionDecay.cxx
@@ -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;
diff --git a/NPSimulation/Process/FissionDecay.cc b/NPSimulation/Process/FissionDecay.cc
index a3e8a1509369e2d7e669f7b87cb6d0135dd9b8c0..45cd4a22d6bdf39fe9fa6a58a9acb1a187d62f10 100644
--- a/NPSimulation/Process/FissionDecay.cc
+++ b/NPSimulation/Process/FissionDecay.cc
@@ -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);
diff --git a/Projects/Sofia/sofia.detector b/Projects/Sofia/sofia.detector
index 80ae67ba7cf207f40388763e72de87e3f21fb44e..a10a046b1bbefd7fc7b675e1bcf62c7ea557ccbb 100644
--- a/Projects/Sofia/sofia.detector
+++ b/Projects/Sofia/sofia.detector
@@ -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
diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
index 6a429bb2052dcd3373e1a7f1e01396cfa8880f72..89cc9dffb3f025d9c57ea7cc3d40708517ac0422 100644
--- a/Projects/e793s/macro/CS2.h
+++ b/Projects/e793s/macro/CS2.h
@@ -1,6 +1,7 @@
/* 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 ;
diff --git a/Projects/e793s/macro/DrawPlots.C b/Projects/e793s/macro/DrawPlots.C
index 290f6a48b998a9f389e95aea98ecff8d2392b94a..a7d5258fd3971e624b963a07a451ed37e3ea7d93 100755
--- a/Projects/e793s/macro/DrawPlots.C
+++ b/Projects/e793s/macro/DrawPlots.C
@@ -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;
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 55122547f19667b1ee35a2117aa8cd710b569ea2..3e7ee21499cb88287dc4e92b08d8007981bc4543 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -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();
+}
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index c50f71d2cb5ff7300ff804477ab392d33dd0a9b3..c70c4d6d06a68935caa978f520bb7eca77588790 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -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;
}
diff --git a/Projects/e793s/macro/ThreeBodyBreakup.h b/Projects/e793s/macro/ThreeBodyBreakup.h
index 5bdbe366e5519a29bc968ee0aecbc199d02dda8c..3df02dc9cc8c2d2b1211a0a78ac0c01cafa63e61 100644
--- a/Projects/e793s/macro/ThreeBodyBreakup.h
+++ b/Projects/e793s/macro/ThreeBodyBreakup.h
@@ -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);