diff --git a/Projects/e793s/macro/CS2.h b/Projects/e793s/macro/CS2.h
deleted file mode 100644
index 3949ca6d5db02abc52ccbad7bf1f95351c325bea..0000000000000000000000000000000000000000
--- a/Projects/e793s/macro/CS2.h
+++ /dev/null
@@ -1,941 +0,0 @@
-/* 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);
-TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment);
-TList* peakFitList = new TList();
-
-/* Global variables */
-vector<Double_t> anglecentres, anglewidth;
-TGraph* currentThry;
-TGraphErrors* staticExp;
-int indexE;
-double globalS, globalSerr;
-int numAngleBins = 20;
-vector<int> removePoint;
-
-/* Output volume toggle */
-bool loud = 0;//1;
-
-/* Scale method toggle */
-bool scaleTogether = 1;
-
-/* Strings for image */
-string orbitalname;
-string orbital;
-
-/* Strings for SolidAngle input file */
-string statename;
-string inputdate;
-
-
-////////////////////////////////////////////////////////////////////////////////
-void canclone(TCanvas* major, int padNum, string name){
- string minorName = "./CS2_Figures/" + name + ".root";
- cout << minorName << endl;
- TFile* minorFile = new TFile(minorName.c_str(),"READ");
- TCanvas* minor = (TCanvas*) minorFile->FindObjectAny("c_peakFits");
-
- major->cd(padNum);
- TPad *pad = (TPad*)gPad;
- minor->cd();
- TObject *obj, *clone;
- pad->Range(minor->GetX1(),minor->GetY1(),minor->GetX2(),minor->GetY2());
- pad->SetTickx(minor->GetTickx());
- pad->SetTicky(minor->GetTicky());
- pad->SetGridx(minor->GetGridx());
- pad->SetGridy(minor->GetGridy());
- pad->SetLogx(minor->GetLogx());
- pad->SetLogy(minor->GetLogy());
- pad->SetLogz(minor->GetLogz());
- pad->SetBorderSize(minor->GetBorderSize());
- pad->SetBorderMode(minor->GetBorderMode());
- minor->TAttLine::Copy((TAttLine&)*pad);
- minor->TAttFill::Copy((TAttFill&)*pad);
- minor->TAttPad::Copy((TAttPad&)*pad);
-
- TIter next(minor->GetListOfPrimitives());
- gROOT->SetSelectedPad(pad);
- while ((obj=next())) {
- clone = obj->Clone();
- obj->Draw("SAME");
- pad->GetListOfPrimitives()->Add(clone,obj->GetDrawOption());
- }
- pad->Modified();
- pad->Update();
- major->cd(padNum);
- pad->Draw();
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void CS_Diagnosis(){
-
- auto majorCanv = new TCanvas("CompareCanv","CompareCanv",1500,1500);
- majorCanv->Divide(3,3);
- canclone(majorCanv, 1, "c_peakFits_110_112");
- canclone(majorCanv, 2, "c_peakFits_115_118");
- canclone(majorCanv, 3, "c_peakFits_120_122");
- canclone(majorCanv, 4, "c_peakFits_125_128");
- canclone(majorCanv, 5, "c_peakFits_130_132");
- canclone(majorCanv, 6, "c_peakFits_135_138");
- canclone(majorCanv, 7, "c_peakFits_140_142");
- canclone(majorCanv, 8, "c_peakFits_145_148");
- canclone(majorCanv, 9, "c_peakFits_150_152");
-
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void CS(){
-/* Overload function */
- cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
- cout << "---- 0.143, p3/2 = CS(0.143, 2, 1, 1.5) "<< endl;
- cout << "---- 0.279, p3/2 = CS(0.279, 2, 1, 1.5) "<< endl;
- cout << "---- 0.728, f7/2 = CS(0.728, 3, 3, 3.5) "<< endl;
- cout << "---- 0.968, p1/2 = CS(0.968, 0, 1, 0.5) "<< endl;
- cout << "---- 1.410, p3/2 = CS(1.410, 1, 1, 1.5) "<< endl;
- cout << "---- 1.981, p3/2 = CS(1.981, 1, 1, 0.5) "<< endl;
- cout << "---- 2.410, p3/2 = CS(2.410, 0, 1, 0.5) "<< endl;
- cout << "---- 3.2 , f7/2 = CS(3.2 , 3, 3, 3.5) "<< endl;
- cout << "---- 3.6 , f5/2 = CS(3.6 , 3, 3, 2.5) "<< endl;
- cout << "---- 3.8 , f5/2 = CS(3.8 , 3, 3, 2.5) "<< endl;
- cout << "---- 4.1 , f5/2 = CS(4.1 , 3, 3, 2.5) "<< endl;
- cout << "---- 4.4 , f5/2 = CS(4.4 , 3, 3, 2.5) "<< endl;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void CS(double Energy, double Spin, double spdf, double angmom){
- // p3/2 -> spdf = 1, angmom = 1.5
- // J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
- double J0 = 0.5;
- //double ElasticNorm = 5.8, ElasticNormErr = 1.3; // DeuteronNorm in elastics, 5.8 +- 1.3
- //double ElasticNorm = 2.70, ElasticNormErr = 0.46; // For thicker
- //double ElasticNorm = 2.8, ElasticNormErr = 0.7; // DeuteronNorm in elastics, 2.8 +- 0.7
- //double ElasticNorm = 2.363, ElasticNormErr = 0.000; // DeuteronNorm in elastics, reconstructed from determined thickness
- //double ElasticNorm = 3.7, ElasticNormErr = 0.000; // Estimated 'goal' normalization
-
- //double ElasticNorm = 0.000220, ElasticNormErr = 0.000; //14Oct22
- double ElasticNorm = 0.000234, ElasticNormErr = 0.000; //18Oct22
- inputdate = "18Oct22";
-
- // Extract spin orbit name
- orbitalname.clear();
- orbital.clear();
- if(spdf==1){
- if(angmom==0.5){
- orbitalname="p_{1/2}";
- orbital="p12";
- } else if (angmom==1.5){
- orbitalname="p_{3/2}";
- orbital="p32";
- } else {
- orbitalname;
- orbital="???";
- }
- } else if (spdf==3){
- if(angmom==2.5){
- orbitalname="f_{5/2}";
- orbital="f52";
- } else if (angmom==3.5) {
- orbitalname="f_{7/2}";
- orbital="f72";
- } else {
- orbitalname="?????";
- orbital="???";
- }
- } else {
- orbitalname="?????";
- orbital="???";
- }
-
- // Number of nodes
- double nodes;
- if(spdf==1){
- nodes=1;
- }
- else if(spdf==3){
- nodes=0;
- }
- else{
- cout << " INPUT NODES::" << endl;
- cin >> nodes;
- }
-
- /* Clean global variables, in case of multiple runs */
- indexE = 100;
- anglecentres.clear();
- anglewidth.clear();
- globalS=0.;
- globalSerr=0.;
- peakFitList->Clear();
-
- /* Retrieve array index of the entered peak energy */
- /* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
- bool found = 0;
- for(int i=0;i<numPeaks;i++){
- if(abs(Energy-means[i])<0.01){
- cout << "========================================================" << endl;
- cout << "Identified as state #" << i << ", E = " << means[i] << endl;
- indexE = i;
- found = 1;
- stringstream ss;
- ss << setfill('0') << setw(4) << (int)(means[i]*1000.);
- statename = ss.str();
- }
- }
- if(!found){
- cout << "========================================================" << endl;
- cout << "NO STATE AT THAT ENERGY INDENTIFIED!! CHECK KNOWN PEAKS!!" << endl;
- return;
- }
-
- /* Solid Angle (from simulation) */
- //auto file = new TFile("SolidAngle_HistFiles/SolidAngle_HistFile_10Aug22_TrueStripRemoval.root");
- string backupFileName = "SolidAngle_HistFiles/SolidAngle_HistFile_10Aug22_TrueStripRemoval.root";
- string saFileName = "SolidAngle_HistFiles/SAHF_";
- saFileName.append(inputdate);
- saFileName.append("_47Kdp_");
- saFileName.append(statename);
- saFileName.append(".root");
-
- TFile* file;
- ifstream infile(saFileName);
- ifstream backup(backupFileName);
- if(infile.good()){
- cout << "Opening file " << saFileName << endl;
- file = TFile::Open(saFileName.c_str());
- } else if (backup.good()){
- cout << BOLDRED << "FAILED TO OPEN " << saFileName << endl;
- cout << "Open standard backup file " << backupFileName << endl;
- cout << RESET;
- file = TFile::Open(backupFileName.c_str());
- } else {
- cout << BOLDRED << "FAILED TO OPEN MAIN OR BACKUP SOLID ANGLE FILE" << endl;
- cout << RED << "Check SolidAngle file exists..." << RESET << endl;
- }
-
- //auto file = new TFile("../SolidAngle_HistFile_New.root");
- TH1F* SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_Lab_MG");
- TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
- SolidAngle->Draw("HIST");
- SolidAngle->GetXaxis()->SetRangeUser(100.,160.);
- /* (canvas deleted after Area/SA calculation) */
-
- /* Area of experimental peaks */
- TCanvas* c_PeakArea = new TCanvas("c_PeakArea","c_PeakArea",1000,1000);
- vector<vector<double>> areaArray = GetExpDiffCross(means[indexE]);
- delete c_PeakArea;
-
- // Array: peakenergy, peakarea, areaerror, anglemin, anglemax
- if(loud){
- for(int i=0; i<areaArray.size();i++){
- cout << i << " "
- << areaArray[i][0] << " "
- << areaArray[i][1] << " "
- << areaArray[i][2] << " "
- << areaArray[i][3] << " "
- << areaArray[i][4] << endl;
- }
- }
-
- /* AoSA = Area/Solid Angle [#/msr] */
- /* dSdO = Experimental Diff. Cross Sect. (Area/Solid Angle *Norm) [mb/msr] */
- vector<Double_t> AoSA, AoSAerr;
- vector<Double_t> expDCS, expDCSerr;
- for(int i=0; i<areaArray.size();i++){
- int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001);
- int binmax = SolidAngle->FindBin(areaArray[i][4]-0.0001);
-
-//cout << "binmin " << binmin << " to binmax " << binmax << endl;
-
- 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
- SAerr = SAerr*1000.; //sr->msr
-
- /* Area over Solid angle ONLY */
- AoSA.push_back(areaArray[i][1]/SA);
- AoSAerr.push_back((areaArray[i][1]/SA)
- * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
- + pow(SAerr/SA,2)));
-
- /* Differential Cross Section */
- /* NOTE: DON'T INCLUDE NORM ERROR IN ERR PROPAGATION AS IT IS SYSTEMATIC! */
- double tempvalue = (areaArray[i][1]/SA)*ElasticNorm;
- expDCS.push_back(tempvalue);
- double temperror = tempvalue
- * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
- + pow(SAerr/SA,2)
- //+ pow(ElasticNormErr/ElasticNorm,2)
- );
- expDCSerr.push_back(temperror);
-
- if(loud){
- cout << "Angle " << areaArray[i][3] << " to " << areaArray[i][4]
- << ", centre " << anglecentres[i]
- << ": Area = " << areaArray[i][1] << " +- " << areaArray[i][2] << " cnts"
- << " SA = " << SA << " +- " << SAerr << " msr"
- << " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"
- << setprecision(5)
- << " Norm = " << ElasticNorm << " +- " << ElasticNormErr
- << " (don't include norm err in propagation)"
- << " dSdO = " << tempvalue << " +- " << temperror
- << setprecision(3)
- << endl;
- }
-
- // Exclude bad angles
- if(AoSA[i]<1E-6 ||AoSAerr[i]>1e2){
- removePoint.push_back(i);
- }
-
- }
- //delete c_SolidAngle;
-
- /* Graph of Area/Solid Angle*/
- TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
- c_AoSA->SetLogy();
- TGraphErrors* gAoSA = new TGraphErrors(
- anglecentres.size(), //n
- &(anglecentres[0]), &(AoSA[0]), //x, y
- //&(anglewidth[0]), &(AoSAerr[0]) );
- 0, &(AoSAerr[0]) ); //errX, errY
- gAoSA->SetTitle("Area/SolidAngle");
- gAoSA->GetXaxis()->SetTitle("ThetaLab [deg]");
- gAoSA->GetYaxis()->SetTitle("Counts/#Omega [counts/msr]");
- gAoSA->Draw();
-
- /* Graph of experimental diff. cross sect (dSigma/dOmega) */
- TCanvas* c_dSdO = new TCanvas("c_dSdO","c_dSdO",1000,1000);
- c_dSdO->SetLogy();
- TGraphErrors* gdSdO = new TGraphErrors(
- anglecentres.size(),
- &(anglecentres[0]), &(expDCS[0]),
- //&(anglewidth[0]), &(AoSAerr[0]) );
- 0, &(expDCSerr[0]) );
- gdSdO->SetTitle("Differential Cross Section");
- gdSdO->GetXaxis()->SetTitle("#theta_{lab} [deg]");
- gdSdO->GetXaxis()->SetTitleOffset(1.2);
- gdSdO->GetXaxis()->SetTitleSize(0.04);
- gdSdO->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
- gdSdO->GetYaxis()->SetTitleOffset(1.2);
- gdSdO->GetYaxis()->SetTitleSize(0.04);
- gdSdO->Draw();
- c_dSdO->Update();
-
-// for(int i=0; i<numAngleBins; i++){
-// if(binary_search(removePoint.begin(), removePoint.end(), i)){
-// cout << BLUE << "removing point " << i << RESET << endl;
-// gAoSA->RemovePoint(i);
-// gdSdO->RemovePoint(i);
-// }
-// }
-
-
- /* 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#Omega [mb/msr]"); //msr set in func above
- TheoryDiffCross->GetXaxis()->SetTitle("ThetaLab [deg]");
- TheoryDiffCross->Draw();
-
- /** TEMP **/
- cout << "UNSCALED THEORY DIFF CROSS SECTION EVALUATED AT DATA POINTS:::" << endl;
- cout << setprecision(6);
- for(int i=0; i<20; i++){
- cout << anglecentres.at(i) << "\t" << TheoryDiffCross->Eval(anglecentres.at(i)) << endl;
-
- }
- cout << setprecision(3);
- cout << "......................" << endl;
- /** **** **/
-
- /* Convert AoSA into Differential Cross Section */
- //cout << " SCALING BY NORMALIZATION = " << ElasticNorm << endl;
- //gAoSA->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
-
-
- /* Scaled and compared on same plot */
-/*
- cout << "USING BASIC SCALING METHOD..." << endl;
- TGraph* ScaledTWOFNR = TheoryDiffCross;
- TCanvas* c_Compare = new TCanvas("c_Compare","c_Compare",1000,1000);
- Scale(ScaledTWOFNR,gAoSA);
- c_Compare->SetLogy();
- gAoSA->SetLineColor(kRed);
- gAoSA->SetMarkerColor(kRed);
- gAoSA->SetMarkerStyle(21);
- gAoSA->Draw("AP");
- ScaledTWOFNR->Draw("SAME");
-*/
-
- /* Using Chi2 minimizaiton */
- cout << "USING CHI2 MINIMIZAITON..." << endl;
- TCanvas* c_Chi2Min = new TCanvas("c_Chi2Min","c_Chi2Min",1000,1000);
- gStyle->SetPadLeftMargin(0.12);
- gStyle->SetPadRightMargin(0.03);
- //c_Chi2Min->SetLogy();
-
- TPad *pad1 = new TPad("pad1","pad1",0,0.25,1,1);
- TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.25);
- pad1->SetTopMargin(0.1);
- pad1->SetBottomMargin(0.00001);
- pad1->SetBorderMode(0);
- pad1->SetLogy();
- pad1->SetTickx();
- pad1->SetTicky();
- pad2->SetTopMargin(0.00001);
- pad2->SetBottomMargin(0.3);
- pad2->SetBorderMode(0);
- pad2->SetTickx();
- pad2->SetTicky();
- pad1->Draw();
- pad2->Draw();
- pad1->cd();
-
- TGraph* Final = FindNormalisation(TheoryDiffCross,gdSdO);
- gdSdO->SetLineColor(kRed);
- gdSdO->SetMarkerColor(kRed);
- gdSdO->SetMarkerStyle(21);
- /* Construct file name string */
- /**/ ostringstream tempstream;
- /**/ if(means[indexE]<1.0){tempstream << 0;}
- /**/ tempstream << (int) (means[indexE]*1000);
- /**/ tempstream << "_" << orbital;
- /**/ tempstream << "_spin" << Spin;
- /**/ string tempstr = tempstream.str();
- /* Construct hist title string */
- /**/ ostringstream textstream;
- /**/ textstream << std::fixed << setprecision(3);
- /**/ textstream << " " << means[indexE];
- /**/ textstream << " MeV, ";
- /**/ textstream << orbitalname;
- /**/ textstream << ", spin " << (int)Spin;
- /**/ textstream << " --> S = " << globalS
- /**/ << " +- " << globalSerr;
- /**/ string textstring = textstream.str();
-
- gdSdO->SetTitle(textstring.c_str());
- gdSdO->GetYaxis()->SetTitleOffset(1.3);
- gdSdO->GetYaxis()->SetTitleSize(0.042);
- gdSdO->GetXaxis()->SetRangeUser(103.,157.);
- //gdSdO->GetXaxis()->SetRangeUser(105.,155.);
- gdSdO->GetXaxis()->SetNdivisions(512,kTRUE);
- gdSdO->Draw("AP");
- Final->Draw("SAME");
-
-
- pad2->cd();
- TGraphErrors* gResid = new TGraphErrors(*gdSdO);
- for(int n=0; n < gResid->GetN(); n++){
- double x = gdSdO->GetPointX(n);
- double residual = gdSdO->GetPointY(n) - Final->Eval(x);
- gResid->SetPoint(n,x,residual);
- gResid->SetPointError(n,0,gdSdO->GetErrorY(n));
- }
- TLine* markzero = new TLine(103.,0.,157.,0.);
- gResid->SetTitle("");
- gResid->GetXaxis()->SetRangeUser(103.,157.);
- //gResid->GetXaxis()->SetRangeUser(105.,155.);
- gResid->GetXaxis()->SetNdivisions(512,kTRUE);
- gResid->GetYaxis()->SetTitle("Residuals");
- gResid->GetYaxis()->SetTitleSize(0.15);
- gResid->GetYaxis()->SetTitleOffset(0.36);
- gResid->GetYaxis()->SetLabelSize(0.08);
- gResid->GetYaxis()->SetNdivisions(305);
- gResid->GetXaxis()->SetTitleSize(0.15);
- gResid->GetXaxis()->SetTitleOffset(0.8);
- gResid->GetXaxis()->SetLabelSize(0.1);
- gResid->GetXaxis()->SetTickLength(0.1);
- gResid->Draw();
- markzero->SetLineStyle(2);
- markzero->Draw("same");
-
- //pad1->cd();
- //TText* orb = new TText(0.5,0.2,"TESTING!!!!!!!!!!!");//orbitalname.c_str());
- //orb->Draw("SAME");
-
- string savestring1 = "./CS2_Figures/"+tempstr+".root";
- string savestring2 = "./CS2_Figures/"+tempstr+".pdf";
- c_Chi2Min->SaveAs(savestring1.c_str());
- c_Chi2Min->SaveAs(savestring2.c_str());
-
- cout << YELLOW
- << " ----- C2S = (2J+1)S = (2*" << (int)Spin << "+1)S = ("
- << (int)((2.*Spin)+1.) << ")S = "
- << ((2.*Spin)+1.) * globalS << " ----- "
- << RESET << endl;
-
-
- //delete c_AoSA;
- //delete c_dSdO;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-vector<vector<double>> GetExpDiffCross(double Energy){
- cout << "========================================================" << endl;
- vector<vector<double>> AllPeaks_OneGate;
- vector<vector<double>> OnePeak_AllGates;
- /****CHANGE ANGLE GATING****/
- double widthAngleBins = 2.5;
- double firstAngle = 105.;
- /***************************/
- double x[numAngleBins], y[numAngleBins];
- //TList* list = new TList();
-
- /* Determine scaling factor for PhaseSpace */
- TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
- double trackScale = 0.0;
- if(scaleTogether){
- TH1F* baseEx = PullThetaLabHist(0,firstAngle,widthAngleBins);
- TH1F* basePS = PullPhaseSpaceHist(0,firstAngle,widthAngleBins);
- for(int i=1; i<numAngleBins;i++){
- TH1F* addEx = PullThetaLabHist(i,firstAngle,widthAngleBins); baseEx->Add(addEx,1.);
- TH1F* addPS = PullPhaseSpaceHist(i,firstAngle,widthAngleBins); basePS->Add(addPS,1.);
- }
-
- /* Subtract flat background equal to smallest bin in range */
- baseEx->GetXaxis()->SetRange(baseEx->FindBin(-1.),baseEx->FindBin(1.));
- double minValueInRange = baseEx->GetBinContent(baseEx->GetMinimumBin());
- baseEx->GetXaxis()->UnZoom();
- cout << "Subtracting background of " << minValueInRange << endl;
- for(int b=1; b<baseEx->GetNbinsX() ; b++){
- baseEx->SetBinContent(b,baseEx->GetBinContent(b)-minValueInRange);
- }
-
- /* Begin scaling within range, track changes */
- basePS->Scale(0.1);
- trackScale = 0.1;
- int numAngleBinsScale = baseEx->GetNbinsX();
- int nbinlow = basePS->FindBin(4.); int nbinhigh = basePS->FindBin(8.0);
- for(int b=nbinlow; b<nbinhigh; b++){
- if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
- while(basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
- basePS->Scale(0.99999);
- trackScale *= 0.99999;
- }
- }
- }
- baseEx->Add(basePS,-1.);
- baseEx->SetName("ExSubPSpace");
- baseEx->SetTitle("ExSubPSpace");
- baseEx->Draw();
- cout << "PhaseSpace -> ExpData scaling = " << trackScale << endl;
- }
-
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- // TEMPORARY!!! REMOVE LAST THREE BINS ON HIGH ENERGY STATES!!!
- //if(means[indexE] > 3.0){numAngleBins-=3;}
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
-
-
- for(int i=0; i<numAngleBins;i++){
- double min = firstAngle + (i*widthAngleBins);
- double max = min + widthAngleBins;
- cout << "===================================" << endl;
- cout << "min: " << min << " max: " << max << endl;
-
- stringstream tmp; tmp << fixed << setprecision(0);
- tmp << "c_peakFits_" << min << "_" << max;
- string tmp2 = tmp.str();
- TCanvas* c_peakFits = new TCanvas("c_peakFits",tmp2.c_str(),1000,1000);
-
- /* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
- /* To change angle gates, run GateThetaLab_MultiWrite() */
- TH1F* gate = PullThetaLabHist(i,firstAngle,widthAngleBins);
- TH1F* pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins);
-
- /* Scale the Phase Space at this angle... */
- /* ... for all angles together */
- if(scaleTogether){
- gate->Add(pspace,-trackScale);
- }
- /* ... or seperately for each angular bin */
- /* NOTE THAT THIS DOES NOT ACCOUNT FOR FLAT BACKGROUND */
- else {
- if(pspace->Integral() > 50.){ // Non-garbage histogram
- pspace->Scale(0.01);
- trackScale=0.01;
- int numAngleBins = gate->GetNbinsX();
- for(int b=0; b<numAngleBins; b++){
- if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
- " > " << gate->GetBinContent(b);
- }
- while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
- pspace->Scale(0.9999);
- trackScale*=0.9999;
- }
- if(loud){cout << " TO " << pspace->GetBinContent(b) <<
- " > " << gate->GetBinContent(b) << endl;
- }
- }
- cout << " !!! SCALE FOR THIS ANGLE = " << trackScale << endl;
- gate->Add(pspace,-1);
- }
- }
-
- /* Subtract flat background equal to smallest bin in range */
- /* ????? */
- /*
- gate->GetXaxis()->SetRange(gate->FindBin(-1.),gate->FindBin(1.));
- double minValueInRange = gate->GetBinContent(gate->GetMinimumBin());
- gate->GetXaxis()->UnZoom();
- cout << "Subtracting background of " << minValueInRange << endl;
- for(int b=1; b<gate->GetNbinsX() ; b++){
- gate->SetBinContent(b,gate->GetBinContent(b)-minValueInRange);
- }
- */
-
- /* 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, 0.0); cout << "!!!!!!! NO SLIDING SHIFT!!!!!" << endl;
- //AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0); cout << "!!!!!!! WITH A VARIABLE SLIDING SHIFT!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
- //double slideshift = 0.00218*(((max-min)/2.)+min) - 0.29645; AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, slideshift); cout << "!!!!!!! WITH A FIXED SLIDING SHIFT!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
-
- /* Write PS-subtracted spectrum to list */
- //list->Add(gate);
- //list->Add(c_peakFits);
- gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
- string filename = "./CS2_Figures/";
- filename += tmp2 + ".root";
- c_peakFits->SaveAs(filename.c_str());
- //auto tempfile = new TFile(filename.c_str(),"UPDATE"); //Reopen newly made file
- //auto templist = new TList();
- //templist->Add();
-
-
- /* Check correct OneGate vector is selected */
- cout << "area of " << means[indexE] << " = "
- << AllPeaks_OneGate[indexE][1]
- << " +- " << AllPeaks_OneGate[indexE][2]
- << endl;
-
- /* Add min and max angle to end of relevant OneGate vector */
- AllPeaks_OneGate[indexE].push_back(min);
- AllPeaks_OneGate[indexE].push_back(max);
-
- /* Push relevant OneGate vector to end of AllGates vector */
- OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
- delete c_peakFits;
- }
-
- /* 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_11Jul22.root","READ");
- //TFile* file = new TFile("GateThetaLabHistograms_10Aug22_TrueStripRemoval.root","READ");
- //TFile* file = new TFile("GateThetaLabHistograms_29Aug22_TrueStripRemoval_0p05.root","READ");
- //TFile* file = new TFile("GateThetaLabHistograms_22Sep22_NoRun51-52.root","READ");
- TFile* file = new TFile("GateThetaLabHistograms_47Kdp_18Oct22_bin0p2.root","READ");
-
- string histname = "cThetaLabGate_"
- + 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");
- //TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_2p5degAngles_26May22v2.root","READ");
- //TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_11Jul22.root","READ");
- TFile* file = new TFile("GatePhaseSpaceThetaLabHistograms_29Aug22_TrueStripRemoval_0p05.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;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-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);
- }
- }
-
- //scaleTWOFNR = mean;
- cout << "SCALED THEORY BY " << mean << endl;
- cout << " therefore S = " << 1/mean << " ??" << endl;
-
- 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, const char* model){
-TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
- /* This function mved between directories in order to run TWOFNR in proper *
- * location. This is, weirdly, the least tempremental way of doing this. */
-
- cout << "========================================================" << endl;
- int johnson, tandyval;
- cout << "Using Johnson-Soper ..."; johnson=5; tandyval=0;
- //cout << "Using Johnson-Tandy 1 ..."; johnson=6; tandyval=1;
- //cout << "Using Johnson-Tandy 2 ..."; johnson=6; tandyval=2;
- //cout << "Using Johnson-Tandy 3 ..."; johnson=6; tandyval=3;
- //cout << "Using Johnson-Tandy 4 ..."; johnson=6; tandyval=4;
-
- int modelA,modelB;
-// switch (model):{
-// case 'K': case 'k':{
-// cout << " ... Koning-Delaroche." << endl; modelA=6; modelB=4;
-// }
-// case 'C': case 'c':{
- cout << " ... and Chapel-Hill." << endl; modelA=2; modelB=2;
-// }
-// case 'B': case 'b':{
-// cout << " ... Bechetti-Greenlees." << endl; modelA=1; modelB=1;
-// }
-// }
-
-
- char origDirchar[200];
- getcwd(origDirchar,200);
- string origDir{origDirchar};
- string twofnrDir = "/home/charlottepaxman/Programs/TWOFNR";
- cout << "Current directory " << origDir << endl;
- cout << "Moving to directory " << twofnrDir << endl;
- chdir(twofnrDir.c_str());
- //Check
- system("pwd");
- cout << "===================================" << endl;
-
- /* Delete existing tran.jjj & 24.jjj files */
- remove("tran.jjj");
- remove("24.jjj");
-
- double BeamEnergy = 7.7;
- double QValue = 2.274 - E;
-
- std::ofstream Front_Input("in.front");
- Front_Input << "jjj" << std::endl;
- Front_Input << "pipo" << std::endl;
- Front_Input << 2 << std::endl;
- Front_Input << 0 << std::endl;
- Front_Input << 0 << std::endl;
- Front_Input << BeamEnergy << std::endl;
- Front_Input << 47 << " " << 19 << 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;
- Front_Input << QValue << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << J0 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << johnson << std::endl;
- if(johnson==6){//JTandy selected, give version
- Front_Input << tandyval << std::endl;
- }
- Front_Input << 1 << std::endl;
- Front_Input << J << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << modelA << std::endl;
- Front_Input << modelB << std::endl;
- Front_Input << 1 << 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();
-
- cout << "Filled front20 input file." << endl;
- cout << "Executing front20..." << endl;
- system("(exec ~/Programs/TWOFNR/front20 < in.front > /dev/null)");
- ifstream checkfront("tran.jjj");
- if(!checkfront){
- //front20 fail!
- cout << " !! ERROR !! \n front20 failed to complete" << endl;
- return 0;
- } else {
- cout << "-> tran.jjj generated (success not guaranteed)" << endl;
- checkfront.close();
- }
-
- /* in.twofnr instructs twofnr20 to evaluate tran.jjj */
- cout << "Executing twofnr20..." << endl;
- system("(exec ~/Programs/TWOFNR/twofnr20 < in.twofnr > /dev/null)");
- ifstream checktwofnr("24.jjj");
- if(!checktwofnr){
- //twofnr20 fail!
- cout << " !! ERROR !! \n twofnr20 failed to complete" << endl;
- terminate();
-// return 0;
- } else {
- cout << "-> twofnr20 complete!" << endl;
- checktwofnr.close();
- }
-
- TGraph* CS = new TGraph("24.jjj");
-
- //mb/sr->mb/msr is x1/1000
- for(int i=0; i<CS->GetN(); i++){
- double x, newy;
- CS->GetPoint(i,x,newy);
- newy=newy/1000;
- CS->SetPoint(i,x,newy);
- }
-
-
- cout << "===================================" << endl;
- cout << "Current directory " << twofnrDir << endl;
- cout << "Moving to directory " << origDir << endl;
- chdir(origDir.c_str());
- system("pwd");
- cout << "========================================================" << endl;
-
- return CS;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double Chi2(TGraph* theory, TGraphErrors* exper){
- double Chi2 = 0;
- double chi = 0;
-
- //cout << setprecision(8);
- //for(int i = 1 ; i < exper->GetN() ; i++){
- for(int i = 0 ; i < exper->GetN() ; i++){
- if(exper->GetPointY(i)>1.0e-8){ //0){
- //chi=(exper->Eval(anglecentres[i])-theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
- chi=(exper->GetPointY(i) - theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
- //cout << "COMPARE::::: " << exper->Eval(anglecentres[i]) << " TO " << exper->GetPointY(i) << endl;
- Chi2 +=chi*chi;
- }
- }
- if(loud){cout << "Chi2 = " << Chi2 << endl;}
- return Chi2;
- //cout << setprecision(3);
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double ToMininize(const double* parameter){
- static int f = 0 ;
- TGraph* g = new TGraph();
- double* X = currentThry->GetX(); // gets valies from global g1 = tgraph passed to find norm
- double* Y = currentThry->GetY();
- for(int i = 0 ; i < currentThry->GetN() ; i++){
- g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i]); // scales this tgraph by parameter
- }
- double chi2 = Chi2(g,staticExp); //compares theory tgraph to experimental tgrapherrors by chi2
- return chi2;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){
- /* (dSdO)meas = S * (dSdO)calc */
-
- // Set global variable
- currentThry = theory;
- staticExp = experiment;
-
- // Construct minimiser
- const char* minName ="Minuit";const char* algoName="Migrad";
- ROOT::Math::Minimizer* min =
- ROOT::Math::Factory::CreateMinimizer(minName, algoName);
- min->SetValidError(true);
-
- // Number of parameters (should only be 1 for me)
- int mysize = 1;
-
- // Create funciton wrapper for minimizer
- // a IMultiGenFunction type
- ROOT::Math::Functor f(&ToMininize,mysize);
- min->SetFunction(f);
-
- // Set range of parameter(??)
- double* parameter = new double[mysize];
- for(unsigned int i = 0 ; i < mysize ; i++){
- parameter[i] = 0.8;
- char name[4];
- sprintf(name,"S%d",i+1);
- min->SetLimitedVariable(i,name,parameter[i],0.01,0,10);
- }
-
-
- ///// TO IMPROVE: FIND WAY OF OBTAINING NDF AND PRINT CHI2/NDF /////
-
- // Minimise
- min->Minimize();
- const double *xs = min->X();
- const double *err = min->Errors();
-
- // Write out
- for(int i = 0 ; i < mysize ; i++){
- cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
- }
-
- /* Store S value in global variable, to access for drawing on plots */
- globalS = xs[0];
- globalSerr = err[0];
-
- // Return the Fitted CS
- TGraph* g = new TGraph();
- double* X = theory->GetX();
- double* Y = theory->GetY();
- if(loud){
- cout << setprecision(8);
- cout << "Start: X[0] = " << theory->GetPointX(4) << " Y[0] = " << theory->GetPointY(4) << endl;
- cout << "multip by " << xs[0] << endl;
- }
-
- //for(int i=0; i<theory->GetN(); i++){ g->SetPoint(g->GetN(),X[i],xs[0]*Y[i]); }
- for(int i=0; i<theory->GetN(); i++){ g->SetPoint(i,X[i],xs[0]*Y[i]); }
-
- if(loud){
- cout << "End: X[0] = " << g->GetPointX(4) << " Y[0] = " << g->GetPointY(4) << endl;
- cout << setprecision(3);
- }
-
- return g;
-}
diff --git a/Projects/e793s/macro/CS2_dt.h b/Projects/e793s/macro/CS2_dt.h
deleted file mode 100644
index 87614b9b695c87058152dcfd1b1090e8c67736ac..0000000000000000000000000000000000000000
--- a/Projects/e793s/macro/CS2_dt.h
+++ /dev/null
@@ -1,924 +0,0 @@
-/* Predefine functions */
-vector<vector<double>> GetExpDiffCross(double Energy);
-TH1F* PullThetaCMHist(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);
-TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment);
-TList* peakFitList = new TList();
-
-/* Global variables */
-vector<Double_t> anglecentres, anglewidth;
-TGraph* currentThry;
-TGraphErrors* staticExp;
-int indexE;
-double globalS, globalSerr;
-int numAngleBins = 18;//10;
-int numAngleBinsBase = 18;//10;
-double widthAngleBins = 1.;//2.;
-double firstAngle = 2.;//3.;
-
-vector<int> removePoint;
-
-/* Output volume toggle */
-bool loud = 1;
-
-/* Scale method toggle */
-bool scaleTogether = 1;
-
-/* String for image */
-string orbitalname;
-string orbital;
-
-/* Strings for SolidAngle input file */
-string statename;
-string inputdate;
-
-////////////////////////////////////////////////////////////////////////////////
-void canclone(TCanvas* major, int padNum, string name){
- string minorName = "./CS2_Figures/" + name + ".root";
- cout << minorName << endl;
- TFile* minorFile = new TFile(minorName.c_str(),"READ");
- TCanvas* minor = (TCanvas*) minorFile->FindObjectAny("c_peakFits");
-
- major->cd(padNum);
- TPad *pad = (TPad*)gPad;
- minor->cd();
- TObject *obj, *clone;
- pad->Range(minor->GetX1(),minor->GetY1(),minor->GetX2(),minor->GetY2());
- pad->SetTickx(minor->GetTickx());
- pad->SetTicky(minor->GetTicky());
- pad->SetGridx(minor->GetGridx());
- pad->SetGridy(minor->GetGridy());
- pad->SetLogx(minor->GetLogx());
- pad->SetLogy(minor->GetLogy());
- pad->SetLogz(minor->GetLogz());
- pad->SetBorderSize(minor->GetBorderSize());
- pad->SetBorderMode(minor->GetBorderMode());
- minor->TAttLine::Copy((TAttLine&)*pad);
- minor->TAttFill::Copy((TAttFill&)*pad);
- minor->TAttPad::Copy((TAttPad&)*pad);
-
- TIter next(minor->GetListOfPrimitives());
- gROOT->SetSelectedPad(pad);
- while ((obj=next())) {
- clone = obj->Clone();
- obj->Draw("SAME");
- pad->GetListOfPrimitives()->Add(clone,obj->GetDrawOption());
- }
- pad->Modified();
- pad->Update();
- major->cd(padNum);
- pad->Draw();
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void CS_Diagnosis(){
-
- auto majorCanv = new TCanvas("CompareCanv","CompareCanv",1500,1500);
- majorCanv->Divide(3,3);
- canclone(majorCanv, 1, "c_peakFits_110_112");
- canclone(majorCanv, 2, "c_peakFits_115_118");
- canclone(majorCanv, 3, "c_peakFits_120_122");
- canclone(majorCanv, 4, "c_peakFits_125_128");
- canclone(majorCanv, 5, "c_peakFits_130_132");
- canclone(majorCanv, 6, "c_peakFits_135_138");
- canclone(majorCanv, 7, "c_peakFits_140_142");
- canclone(majorCanv, 8, "c_peakFits_145_148");
- canclone(majorCanv, 9, "c_peakFits_150_152");
-
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void CS(){
-/* Overload function */
- cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
- cout << "---- 1.945, s1/2 = CS(1.945, 1, 0, 0.5) "<< endl;
- cout << "---- 1.945, d3/2 = CS(1.945, 1, 2, 1.5) "<< endl;
- cout << "---- 0.691, f7/2 = CS(0.691, 4, 3, 3.5) "<< endl;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-void CS(double Energy, double Spin, double spdf, double angmom){
- // p3/2 -> spdf = 1, angmom = 1.5
- // J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
- double J0 = 0.5;
-
- //double ElasticNorm = 0.000220, ElasticNormErr = 0.000; //14Oct22
- double ElasticNorm = 0.000234, ElasticNormErr = 0.000016; //18Oct22
- inputdate = "18Oct22";
-
- // Extract spin orbit name
- orbitalname.clear();
- orbital.clear();
- if(spdf==0){
- if(angmom==0.5){
- orbitalname="s_{1/2}";
- orbital="s12";
- } else {
- orbitalname="?????";
- orbital="???";
- }
- } else if (spdf==2){
- if(angmom==1.5){
- orbitalname="d_{3/2}";
- orbital="d32";
- } else if (angmom==2.5) {
- orbitalname="d_{5/2}";
- orbital="d52";
- } else {
- orbitalname="?????";
- orbital="???";
- }
- } else if (spdf==3){
- if(angmom==3.5){
- orbitalname="f_{7/2}";
- orbital="f72";
- } else {
- orbitalname="?????";
- orbital="???";
- }
- } else {
- orbitalname="?????";
- orbital="???";
- }
-
- // Number of nodes
- double nodes;
- if(spdf==0){
- nodes=1;
- }
- else if(spdf==3 || spdf == 2){
- nodes=0;
- }
- else{
- cout << " INPUT NODES::" << endl;
- cin >> nodes;
- }
-
- /* Clean global variables, in case of multiple runs */
- indexE = 100;
- anglecentres.clear();
- anglewidth.clear();
- globalS=0.;
- globalSerr=0.;
- peakFitList->Clear();
- numAngleBins=numAngleBinsBase;
-
- /* Retrieve array index of the entered peak energy */
- /* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
- bool found = 0;
- for(int i=0;i<numPeaks;i++){
- if(abs(Energy-means_dt[i])<0.01){
- cout << "========================================================" << endl;
- cout << "Identified as state #" << i << ", E = " << means_dt[i] << endl;
- indexE = i;
- found = 1;
- stringstream ss;
- ss << setfill('0') << setw(4) << (int)(means_dt[i]*1000.);
- statename = ss.str();
- }
- }
- if(!found){
- cout << "========================================================" << endl;
- cout << "NO STATE AT THAT ENERGY INDENTIFIED!! CHECK KNOWN PEAKS!!" << endl;
- return;
- }
-
- /* Solid Angle (from simulation) */
- //auto file = new TFile("../SolidAngle_HistFile_18Oct22_47Kdt.root");
- string backupFileName = "SolidAngle_HistFiles/SolidAngle_HistFile_18Oct22_47Kdt.root";
- string saFileName = "SolidAngle_HistFiles/SAHF_";
- saFileName.append(inputdate);
- saFileName.append("_47Kdt_");
- saFileName.append(statename);
- saFileName.append(".root");
-
- TFile* file;
- ifstream infile(saFileName);
- ifstream backup(backupFileName);
- if(infile.good()){
- cout << "Opening file " << saFileName << endl;
- file = TFile::Open(saFileName.c_str());
- } else if (backup.good()){
- cout << BOLDRED << "FAILED TO OPEN " << saFileName << endl;
- cout << "Open standard backup file " << backupFileName << endl;
- cout << RESET;
- file = TFile::Open(backupFileName.c_str());
- } else {
- cout << BOLDRED << "FAILED TO OPEN MAIN OR BACKUP SOLID ANGLE FILE" << endl;
- cout << RED << "Check SolidAngle file exists..." << RESET << endl;
- }
-
- TH1F* SolidAngle = (TH1F*) file->FindObjectAny("SolidAngle_CM_MM");
- TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
- SolidAngle->Draw("HIST");
- SolidAngle->GetXaxis()->SetRangeUser(00.,30.);
- /* (canvas deleted after Area/SA calculation) */
-
- /* Area of experimental peaks */
- TCanvas* c_PeakArea = new TCanvas("c_PeakArea","c_PeakArea",1000,1000);
- vector<vector<double>> areaArray = GetExpDiffCross(means_dt[indexE]);
- delete c_PeakArea;
-
- // Array: peakenergy, peakarea, areaerror, anglemin, anglemax
- if(loud){
- for(int i=0; i<areaArray.size();i++){
- cout << i << " "
- << areaArray[i][0] << " "
- << areaArray[i][1] << " "
- << areaArray[i][2] << " "
- << areaArray[i][3] << " "
- << areaArray[i][4] << endl;
- }
- }
-
- /* AoSA = Area/Solid Angle [#/msr] */
- /* dSdO = Experimental Diff. Cross Sect. (Area/Solid Angle *Norm) [mb/msr] */
- vector<Double_t> AoSA, AoSAerr;
- vector<Double_t> expDCS, expDCSerr;
- for(int i=0; i<areaArray.size();i++){
- int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001);
- int binmax = SolidAngle->FindBin(areaArray[i][4]-0.0001);
-
-//cout << "binmin " << binmin << " to binmax " << binmax << endl;
-
- 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
- SAerr = SAerr*1000.; //sr->msr
-
- /* Area over Solid angle (exclude errors) */
- if(SA<=0.){
- SA = -20.; SAerr = -20.;
- AoSA.push_back(-20.); AoSAerr.push_back(-20.);
- } else {
- AoSA.push_back(areaArray[i][1]/SA);
- AoSAerr.push_back((areaArray[i][1]/SA)
- * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
- + pow(SAerr/SA,2)));
- }
-
- /* Differential Cross Section */
- /* NOTE: DON'T INCLUDE NORM ERROR IN ERR PROPAGATION AS IT IS SYSTEMATIC! */
- double tempvalue = (areaArray[i][1]/SA)*ElasticNorm;
- expDCS.push_back(tempvalue);
- double temperror = tempvalue
- * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
- + pow(SAerr/SA,2)
- //+ pow(ElasticNormErr/ElasticNorm,2)
- );
- expDCSerr.push_back(temperror);
-
- if(loud){
- cout << "Angle " << areaArray[i][3] << " to " << areaArray[i][4]
- << ", centre " << anglecentres[i]
- << ": Area = " << areaArray[i][1] << " +- " << areaArray[i][2] << " cnts"
- << " SA = " << SA << " +- " << SAerr << " msr"
- << " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"
- << setprecision(5)
- << " Norm = " << ElasticNorm << " +- " << ElasticNormErr
- << " (don't include norm err in propagation)"
- << " dSdO = " << tempvalue << " +- " << temperror
- << setprecision(3)
- << endl;
- }
-
- // Exclude bad angles
- if(AoSA[i]<1E-6 ||AoSAerr[i]>1e2){
- removePoint.push_back(i);
- }
-
- }
- //delete c_SolidAngle;
-
- /* Graph of Area/Solid Angle*/
- TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
- c_AoSA->SetLogy();
- TGraphErrors* gAoSA = new TGraphErrors(
- anglecentres.size(), //n
- &(anglecentres[0]), &(AoSA[0]), //x, y
- //&(anglewidth[0]), &(AoSAerr[0]) );
- 0, &(AoSAerr[0]) ); //errX, errY
- gAoSA->SetTitle("Area/SolidAngle");
- gAoSA->GetXaxis()->SetTitle("ThetaCM [deg]");
- gAoSA->GetYaxis()->SetTitle("Counts/#Omega [counts/msr]");
- gAoSA->Draw();
-
- /* Graph of experimental diff. cross sect (dSigma/dOmega) */
- TCanvas* c_dSdO = new TCanvas("c_dSdO","c_dSdO",1000,1000);
- c_dSdO->SetLogy();
- TGraphErrors* gdSdO = new TGraphErrors(
- anglecentres.size(),
- &(anglecentres[0]), &(expDCS[0]),
- //&(anglewidth[0]), &(AoSAerr[0]) );
- 0, &(expDCSerr[0]) );
- gdSdO->SetTitle("Differential Cross Section");
- gdSdO->GetXaxis()->SetTitle("#theta_{lab} [deg]");
- gdSdO->GetXaxis()->SetTitleOffset(1.2);
- gdSdO->GetXaxis()->SetTitleSize(0.04);
- gdSdO->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
- gdSdO->GetYaxis()->SetTitleOffset(1.2);
- gdSdO->GetYaxis()->SetTitleSize(0.04);
- gdSdO->Draw();
- c_dSdO->Update();
-
-// for(int i=0; i<numAngleBins; i++){
-// if(binary_search(removePoint.begin(), removePoint.end(), i)){
-// cout << BLUE << "removing point " << i << RESET << endl;
-// gAoSA->RemovePoint(i);
-// gdSdO->RemovePoint(i);
-// }
-// }
-
-
-
- /* 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_dt[indexE], J0, Spin, nodes, spdf, angmom);
- TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); //msr set in func above
- TheoryDiffCross->GetXaxis()->SetTitle("ThetaCM [deg]");
- TheoryDiffCross->Draw();
-
- /** TEMP **/
- cout << "UNSCALED THEORY DIFF CROSS SECTION EVALUATED AT DATA POINTS:::" << endl;
- cout << setprecision(6);
- for(int i=0; i<numAngleBins; i++){
- cout << anglecentres.at(i) << "\t" << TheoryDiffCross->Eval(anglecentres.at(i)) << endl;
-
- }
- cout << setprecision(3);
- cout << "......................" << endl;
- /** **** **/
-
- /* Convert AoSA into Differential Cross Section */
- //cout << " SCALING BY NORMALIZATION = " << ElasticNorm << endl;
- //gAoSA->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
-
-
- /* Scaled and compared on same plot */
-/*
- cout << "USING BASIC SCALING METHOD..." << endl;
- TGraph* ScaledTWOFNR = TheoryDiffCross;
- TCanvas* c_Compare = new TCanvas("c_Compare","c_Compare",1000,1000);
- Scale(ScaledTWOFNR,gAoSA);
- c_Compare->SetLogy();
- gAoSA->SetLineColor(kRed);
- gAoSA->SetMarkerColor(kRed);
- gAoSA->SetMarkerStyle(21);
- gAoSA->Draw("AP");
- ScaledTWOFNR->Draw("SAME");
-*/
-
- /* Using Chi2 minimizaiton */
- cout << "USING CHI2 MINIMIZAITON..." << endl;
- TCanvas* c_Chi2Min = new TCanvas("c_Chi2Min","c_Chi2Min",1000,1000);
- gStyle->SetPadLeftMargin(0.12);
- gStyle->SetPadRightMargin(0.03);
- //c_Chi2Min->SetLogy();
-
- TPad *pad1 = new TPad("pad1","pad1",0,0.25,1,1);
- TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.25);
- pad1->SetTopMargin(0.1);
- pad1->SetBottomMargin(0.00001);
- pad1->SetBorderMode(0);
- pad1->SetLogy();
- pad1->SetTickx();
- pad1->SetTicky();
- pad2->SetTopMargin(0.00001);
- pad2->SetBottomMargin(0.3);
- pad2->SetBorderMode(0);
- pad2->SetTickx();
- pad2->SetTicky();
- pad1->Draw();
- pad2->Draw();
- pad1->cd();
-
-
- double graphMin = (double) floor(gdSdO->GetPointX(0));
- double graphMax = (double) ceil(gdSdO->GetPointX(numAngleBins-1));
- cout << RED << graphMin << " " << graphMax << RESET << endl;
-
- TGraph* Final = FindNormalisation(TheoryDiffCross,gdSdO);
- gdSdO->SetLineColor(kRed);
- gdSdO->SetMarkerColor(kRed);
- gdSdO->SetMarkerStyle(21);
- /* Construct file name string */
- /**/ ostringstream tempstream;
- /**/ if(means_dt[indexE]<1.0){tempstream << 0;}
- /**/ tempstream << (int) (means_dt[indexE]*1000);
- /**/ tempstream << "_" << orbital;
- /**/ tempstream << "_spin" << Spin;
- /**/ string tempstr = tempstream.str();
- /* Construct hist title string */
- /**/ ostringstream textstream;
- /**/ textstream << std::fixed << setprecision(3);
- /**/ textstream << " " << means_dt[indexE];
- /**/ textstream << " MeV, ";
- /**/ textstream << orbitalname;
- /**/ textstream << ", spin " << (int)Spin;
- /**/ textstream << " --> S = " << globalS
- /**/ << " +- " << globalSerr;
- /**/ string textstring = textstream.str();
-
- gdSdO->SetTitle(textstring.c_str());
- gdSdO->GetYaxis()->SetTitleOffset(1.3);
- gdSdO->GetYaxis()->SetTitleSize(0.042);
- gdSdO->GetXaxis()->SetRangeUser(graphMin,graphMax);
- gdSdO->Draw("AP");
- Final->Draw("SAME");
-
-
- pad2->cd();
- TGraphErrors* gResid = new TGraphErrors(*gdSdO);
- for(int n=0; n < gResid->GetN(); n++){
- double x = gdSdO->GetPointX(n);
- double residual = gdSdO->GetPointY(n) - Final->Eval(x);
- gResid->SetPoint(n,x,residual);
- gResid->SetPointError(n,0,gdSdO->GetErrorY(n));
- }
- TLine* markzero = new TLine(graphMin,0.,graphMax,0.);
- gResid->SetTitle("");
- gResid->GetXaxis()->SetRangeUser(graphMin,graphMax);
- gResid->GetYaxis()->SetTitle("Residuals");
- gResid->GetYaxis()->SetTitleSize(0.15);
- gResid->GetYaxis()->SetTitleOffset(0.36);
- gResid->GetYaxis()->SetLabelSize(0.08);
- gResid->GetYaxis()->SetNdivisions(305);
- gResid->GetXaxis()->SetTitleSize(0.15);
- gResid->GetXaxis()->SetTitleOffset(0.8);
- gResid->GetXaxis()->SetLabelSize(0.1);
- gResid->GetXaxis()->SetTickLength(0.1);
- gResid->Draw();
- markzero->SetLineStyle(2);
- markzero->Draw("same");
-
- //pad1->cd();
- //TText* orb = new TText(0.5,0.2,"TESTING!!!!!!!!!!!");//orbitalname.c_str());
- //orb->Draw("SAME");
-
- string savestring1 = "./CS2_Figures/"+tempstr+".root";
- string savestring2 = "./CS2_Figures/"+tempstr+".pdf";
- c_Chi2Min->SaveAs(savestring1.c_str());
- c_Chi2Min->SaveAs(savestring2.c_str());
-
- cout << YELLOW
- << " ----- C2S = (2J+1)S = (2*" << (int)Spin << "+1)S = ("
- << (int)((2.*Spin)+1.) << ")S = "
- << ((2.*Spin)+1.) * globalS << " ----- "
- << RESET << endl;
-
-
- //delete c_AoSA;
- //delete c_dSdO;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-vector<vector<double>> GetExpDiffCross(double Energy){
- cout << "========================================================" << endl;
- vector<vector<double>> AllPeaks_OneGate;
- vector<vector<double>> OnePeak_AllGates;
- /***************************/
- double x[numAngleBins], y[numAngleBins];
- //TList* list = new TList();
-
- /* Determine scaling factor for PhaseSpace */
- TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
- double trackScale = 0.0;
- //if(scaleTogether){
- TH1F* baseEx = PullThetaCMHist(0,firstAngle,widthAngleBins);
- // TH1F* basePS = PullPhaseSpaceHist(0,firstAngle,widthAngleBins);
- for(int i=1; i<numAngleBins;i++){
- TH1F* addEx = PullThetaCMHist(i,firstAngle,widthAngleBins); baseEx->Add(addEx,1.);
- // TH1F* addPS = PullPhaseSpaceHist(i,firstAngle,widthAngleBins); basePS->Add(addPS,1.);
- }
-
- /* Subtract flat background equal to smallest bin in range */
- baseEx->GetXaxis()->SetRange(baseEx->FindBin(-1.),baseEx->FindBin(1.));
- double minValueInRange = baseEx->GetBinContent(baseEx->GetMinimumBin());
- baseEx->GetXaxis()->UnZoom();
- cout << "Subtracting background of " << minValueInRange << endl;
- for(int b=1; b<baseEx->GetNbinsX() ; b++){
- baseEx->SetBinContent(b,baseEx->GetBinContent(b)-minValueInRange);
- }
-
- /* Begin scaling within range, track changes */
- //basePS->Scale(0.1);
- //trackScale = 0.1;
- //int numAngleBinsScale = baseEx->GetNbinsX();
- //int nbinlow = basePS->FindBin(4.); int nbinhigh = basePS->FindBin(8.0);
- //for(int b=nbinlow; b<nbinhigh; b++){
- // if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
- // while(basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
- // basePS->Scale(0.99999);
- // trackScale *= 0.99999;
- // }
- // }
- //}
- //baseEx->Add(basePS,-1.);
- baseEx->SetName("Ex");//SubPSpace");
- baseEx->SetTitle("Ex");//SubPSpace");
- baseEx->Draw();
- cout << "PhaseSpace -> ExpData scaling = ZERO! NO PHASE SPACE!" << endl;//" << trackScale << endl;
- //}
-
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- /* TEMPORARY FIX FOR EX:THETACM PROBLEM!! */
- double thetaMax = 0.23*pow(means_dt[indexE],2)
- + 0.76*means_dt[indexE]
- + 9.01;
- thetaMax = floor(thetaMax);
-
- double thetaMin = 0.10*pow(means_dt[indexE],2)
- + 0.06*means_dt[indexE]
- + 1.58;
- thetaMin = ceil(thetaMin);
-
- numAngleBins = (int)thetaMax - (int)thetaMin;
- firstAngle = (int)thetaMin;
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
-
-
- for(int i=0; i<numAngleBins;i++){
- double min = firstAngle + (i*widthAngleBins);
- double max = min + widthAngleBins;
- cout << "===================================" << endl;
- cout << "min: " << min << " max: " << max << endl;
-
- stringstream tmp; tmp << fixed << setprecision(0);
- tmp << "c_peakFits_dt_" << min << "_" << max;
- string tmp2 = tmp.str();
- TCanvas* c_peakFits = new TCanvas("c_peakFits",tmp2.c_str(),1000,1000);
-
- /* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
- /* To change angle gates, run GateThetaLab_MultiWrite() */
- TH1F* gate = PullThetaCMHist(i,firstAngle,widthAngleBins);
- //TH1F* pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins);
-
- /* Scale the Phase Space at this angle... */
- /* ... for all angles together */
- //if(scaleTogether){
- // gate->Add(pspace,-trackScale);
- //}
- /* ... or seperately for each angular bin */
- /* NOTE THAT THIS DOES NOT ACCOUNT FOR FLAT BACKGROUND */
- //else {
- // if(pspace->Integral() > 50.){ // Non-garbage histogram
- // pspace->Scale(0.01);
- // trackScale=0.01;
- // int numAngleBins = gate->GetNbinsX();
- // for(int b=0; b<numAngleBins; b++){
- // if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
- // " > " << gate->GetBinContent(b);
- // }
- // while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
- // pspace->Scale(0.9999);
- // trackScale*=0.9999;
- // }
- // if(loud){cout << " TO " << pspace->GetBinContent(b) <<
- // " > " << gate->GetBinContent(b) << endl;
- // }
- // }
- // cout << " !!! SCALE FOR THIS ANGLE = " << trackScale << endl;
- // gate->Add(pspace,-1);
- // }
- //}
-
- /* Subtract flat background equal to smallest bin in range */
- /* ????? */
- /*
- gate->GetXaxis()->SetRange(gate->FindBin(-1.),gate->FindBin(1.));
- double minValueInRange = gate->GetBinContent(gate->GetMinimumBin());
- gate->GetXaxis()->UnZoom();
- cout << "Subtracting background of " << minValueInRange << endl;
- for(int b=1; b<gate->GetNbinsX() ; b++){
- gate->SetBinContent(b,gate->GetBinContent(b)-minValueInRange);
- }
- */
-
- /* Retrieve array containing all fits, for one angle gate. *
- * Specific peak of interest selected from the vector by *
- * global variable indexE */
-
- AllPeaks_OneGate = FitKnownPeaks_dt_RtrnArry(gate, 0.0); cout << "!!!!!!! NO SLIDING SHIFT!!!!!" << endl;
- //AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0); cout << "!!!!!!! WITH A VARIABLE SLIDING SHIFT!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
- //double slideshift = 0.00218*(((max-min)/2.)+min) - 0.29645; AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, slideshift); cout << "!!!!!!! WITH A FIXED SLIDING SHIFT!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
-
- /* Write PS-subtracted spectrum to list */
- //list->Add(gate);
- //list->Add(c_peakFits);
- gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
- string filename = "./CS2_Figures/";
- filename += tmp2 + ".root";
- c_peakFits->SaveAs(filename.c_str());
- //auto tempfile = new TFile(filename.c_str(),"UPDATE"); //Reopen newly made file
- //auto templist = new TList();
- //templist->Add();
-
-
- /* Check correct OneGate vector is selected */
- cout << "area of " << means_dt[indexE] << " = "
- << AllPeaks_OneGate[indexE][1]
- << " +- " << AllPeaks_OneGate[indexE][2]
- << endl;
-
- /* Add min and max angle to end of relevant OneGate vector */
- AllPeaks_OneGate[indexE].push_back(min);
- AllPeaks_OneGate[indexE].push_back(max);
-
- /* Push relevant OneGate vector to end of AllGates vector */
- OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
- delete c_peakFits;
- }
-
- return OnePeak_AllGates;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-TH1F* PullThetaCMHist(int i, double minTheta, double gatesize){
- //TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p2.root","READ");
- //TFile* file = new TFile("GateThetaCMHistograms_21Oct22_47Kdt.root","READ");
- //TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p1_10angles_0to20.root","READ");
- TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p1.root","READ");
-
- string histname = "cThetaCMGate_"
- + to_string((int) (minTheta+(i*gatesize))) + "-"
- + to_string((int) (minTheta+((i+1)*gatesize)));
- cout << "Loading " << histname << endl;
- TList *list = (TList*)file->Get("GateThetaCMHistograms");
- TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
-// file->Close();
- return hist;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-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);
- }
- }
-
- //scaleTWOFNR = mean;
- cout << "SCALED THEORY BY " << mean << endl;
- cout << " therefore S = " << 1/mean << " ??" << endl;
-
- 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, const char* model){
-TGraph* TWOFNR(double E, double J0, double J, double n, double l, double j){
- /* This function mved between directories in order to run TWOFNR in proper *
- * location. This is, weirdly, the least tempremental way of doing this. */
-
-// cout << "========================================================" << endl;
-// int johnson, tandyval;
-// cout << "Using Johnson-Soper ..."; johnson=5; tandyval=0;
- //cout << "Using Johnson-Tandy 1 ..."; johnson=6; tandyval=1;
- //cout << "Using Johnson-Tandy 2 ..."; johnson=6; tandyval=2;
- //cout << "Using Johnson-Tandy 3 ..."; johnson=6; tandyval=3;
- //cout << "Using Johnson-Tandy 4 ..."; johnson=6; tandyval=4;
-
-// int modelA,modelB;
-//// switch (model):{
-//// case 'K': case 'k':{
-//// cout << " ... Koning-Delaroche." << endl; modelA=6; modelB=4;
-//// }
-// case 'C': case 'c':{
-//// cout << " ... and Chapel-Hill." << endl; modelA=2; modelB=2;
-//// }
-//// case 'B': case 'b':{
-//// cout << " ... Bechetti-Greenlees." << endl; modelA=1; modelB=1;
-//// }
-//// }
-
-
- char origDirchar[200];
- getcwd(origDirchar,200);
- string origDir{origDirchar};
- string twofnrDir = "/home/charlottepaxman/Programs/TWOFNR";
- cout << "Current directory " << origDir << endl;
- cout << "Moving to directory " << twofnrDir << endl;
- chdir(twofnrDir.c_str());
- //Check
- system("pwd");
- cout << "===================================" << endl;
-
- /* Delete existing tran.jjj & 24.jjj files */
- remove("tran.jjj");
- remove("24.jjj");
-
- double BeamEnergy = 7.7;
- double QValue = -2.112 - E;
-
- std::ofstream Front_Input("in.front");
- Front_Input << "jjj" << std::endl;
- Front_Input << "pipo" << std::endl;
- Front_Input << 5 << std::endl;
- Front_Input << 0 << std::endl;
- Front_Input << 0 << std::endl;
- Front_Input << BeamEnergy << std::endl;
- Front_Input << 47 << " " << 19 << 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;
- Front_Input << QValue << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << J0 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << 1 << std::endl;
- Front_Input << J << std::endl;
- Front_Input << 1 << 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();
-
- cout << "Filled front20 input file." << endl;
- cout << "Executing front20..." << endl;
- system("(exec ~/Programs/TWOFNR/front20 < in.front > /dev/null)");
- ifstream checkfront("tran.jjj");
- if(!checkfront){
- //front20 fail!
- cout << " !! ERROR !! \n front20 failed to complete" << endl;
- return 0;
- } else {
- cout << "-> tran.jjj generated (success not guaranteed)" << endl;
- checkfront.close();
- }
-
- /* in.twofnr instructs twofnr20 to evaluate tran.jjj */
- cout << "Executing twofnr20..." << endl;
- system("(exec ~/Programs/TWOFNR/twofnr20 < in.twofnr > /dev/null)");
- ifstream checktwofnr("24.jjj");
- if(!checktwofnr){
- //twofnr20 fail!
- cout << " !! ERROR !! \n twofnr20 failed to complete" << endl;
- terminate();
-// return 0;
- } else {
- cout << "-> twofnr20 complete!" << endl;
- checktwofnr.close();
- }
-
- /* In file read, %lg means read this column, %*lg means ignore column */
- //TGraph* CS = new TGraph("24.jjj");
- TGraph* CS = new TGraph("21.jjj","%lg %lg %*lg");
-
- //mb/sr->mb/msr is x1/1000
- for(int i=0; i<CS->GetN(); i++){
- double x, newy;
- CS->GetPoint(i,x,newy);
- newy=newy/1000;
- CS->SetPoint(i,x,newy);
- }
-
-
- cout << "===================================" << endl;
- cout << "Current directory " << twofnrDir << endl;
- cout << "Moving to directory " << origDir << endl;
- chdir(origDir.c_str());
- system("pwd");
- cout << "========================================================" << endl;
-
- return CS;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double Chi2(TGraph* theory, TGraphErrors* exper){
- double Chi2 = 0;
- double chi = 0;
-
- //cout << setprecision(8);
- //for(int i = 1 ; i < exper->GetN() ; i++){
- for(int i = 0 ; i < exper->GetN() ; i++){
- if(exper->GetPointY(i)>1.0e-8){ //0){
- //chi=(exper->Eval(anglecentres[i])-theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
- chi=(exper->GetPointY(i) - theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
- cout << "COMPARE::::: Thr " << theory->Eval(anglecentres[i]) << " TO Exp " << exper->GetPointY(i) << endl;
- Chi2 +=chi*chi;
- }
- }
- if(loud){cout << "Chi2 = " << Chi2 << endl;}
- return Chi2;
- //cout << setprecision(3);
-}
-
-////////////////////////////////////////////////////////////////////////////////
-double ToMininize(const double* parameter){
- static int f = 0 ;
- TGraph* g = new TGraph();
- double* X = currentThry->GetX(); // gets valies from global g1 = tgraph passed to find norm
- double* Y = currentThry->GetY();
- for(int i = 0 ; i < currentThry->GetN() ; i++){
- g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i]); // scales this tgraph by parameter
- }
- double chi2 = Chi2(g,staticExp); //compares theory tgraph to experimental tgrapherrors by chi2
- return chi2;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){
- /* (dSdO)meas = S * (dSdO)calc */
-
- // Set global variable
- currentThry = theory;
- staticExp = experiment;
-
- // Construct minimiser
- const char* minName ="Minuit";const char* algoName="Migrad";
- ROOT::Math::Minimizer* min =
- ROOT::Math::Factory::CreateMinimizer(minName, algoName);
- min->SetValidError(true);
-
- // Number of parameters (should only be 1 for me)
- int mysize = 1;
-
- // Create funciton wrapper for minimizer
- // a IMultiGenFunction type
- ROOT::Math::Functor f(&ToMininize,mysize);
- min->SetFunction(f);
-
- // Set range of parameter(??)
- double* parameter = new double[mysize];
- for(unsigned int i = 0 ; i < mysize ; i++){
- parameter[i] = 0.8;
- char name[4];
- sprintf(name,"S%d",i+1);
- min->SetLimitedVariable(i,name,parameter[i],0.01,0,20);
- }
-
- ///// TO IMPROVE: FIND WAY OF OBTAINING NDF AND PRINT CHI2/NDF /////
-
- // Minimise
- min->Minimize();
- const double *xs = min->X();
- const double *err = min->Errors();
-
- // Write out
- for(int i = 0 ; i < mysize ; i++){
- cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
- }
-
- /* Store S value in global variable, to access for drawing on plots */
- globalS = xs[0];
- globalSerr = err[0];
-
- // Return the Fitted CS
- TGraph* g = new TGraph();
- double* X = theory->GetX();
- double* Y = theory->GetY();
- if(loud){
- cout << setprecision(8);
- cout << "Start: X[0] = " << theory->GetPointX(4) << " Y[0] = " << theory->GetPointY(4) << endl;
- cout << "multip by " << xs[0] << endl;
- }
-
- //for(int i=0; i<theory->GetN(); i++){ g->SetPoint(g->GetN(),X[i],xs[0]*Y[i]); }
- for(int i=0; i<theory->GetN(); i++){ g->SetPoint(i,X[i],xs[0]*Y[i]); }
-
- if(loud){
- cout << "End: X[0] = " << g->GetPointX(4) << " Y[0] = " << g->GetPointY(4) << endl;
- cout << setprecision(3);
- }
-
- return g;
-}
diff --git a/Projects/e793s/macro/CS2_master.h b/Projects/e793s/macro/CS2_master.h
new file mode 100644
index 0000000000000000000000000000000000000000..e7ea9dc0d67e97d01f9d051af1f1aa97916a891d
--- /dev/null
+++ b/Projects/e793s/macro/CS2_master.h
@@ -0,0 +1,2151 @@
+/* Predefine functions */
+void CS();
+void CS_Diagnosis();
+vector<vector<double>> GetExpDiffCross(double Energy);
+TH1F* PullThetaLabHist(int i, double minTheta, double gatesize);
+TH1F* PullThetaCMHist(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);
+TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment);
+TGraph* FindNormalisation(TGraph* theory, TGraph* theory2, TGraphErrors* experiment);
+TList* peakFitList = new TList();
+
+/* Global variables */
+vector<Double_t> anglecentres, anglewidth;
+TGraph* currentThry;
+TGraph* currentThry2;
+TGraphErrors* staticExp;
+int indexE;
+double globalS, globalSerr;
+double globalS2, globalSerr2;
+int numAngleBins;
+int numAngleBinsBase;
+int numPeaks_CS;
+double widthAngleBins;
+double firstAngle;
+double CSangleL, CSangleH;
+vector<double> means_CS;
+bool GammaGate = false;
+double globGammaGate;
+bool doMix = false;
+bool doDoublet = false;
+double globalBinning;
+
+/* Output volume toggle */
+bool loud = 1;
+
+/* Scale method toggle */
+bool scaleTogether = 1;
+
+/* Strings for image */
+string orbitalname;
+string orbital;
+
+/* Strings for SolidAngle input file */
+string statename;
+string inputdate;
+
+
+////////////////////////////////////////////////////////////////////////////////
+void WriteToFile(TGraphErrors* gr, string filename){
+ int npoints = gr->GetN();
+ ofstream file;
+ file.open (filename.c_str());
+
+ cout << " ============================ " << endl;
+ cout << " Writing to " << filename << " ..." << endl;
+
+ file << " i \t x \t y \t xe \t ye" << endl;
+
+ for(int i = 0; i < npoints; i++){
+ double x = gr->GetPointX(i);
+ double y = gr->GetPointY(i);
+ double xe = gr->GetErrorX(i);
+ double ye = gr->GetErrorY(i);
+
+ file << i << "\t"
+ << x << "\t"
+ << y << "\t"
+ << xe << "\t"
+ << ye << "\t"
+ << endl;
+
+ }
+
+ cout << " Complete! " << endl;
+ cout << " ============================ " << endl;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+TGraphErrors* GetExpDCS_ForFig(string filename){
+ vector<double> x, dx, y, dy;
+ ifstream f(filename.c_str());
+ string line;
+ while (getline(f, line)){
+ stringstream ss(line);
+ // i x y dx dy
+ double tA, tB, tC, tD, tE;
+ if (f >> tA >> tB >> tC >> tD >> tE){
+ //cout << tA << " " << tB << " " << tC << " " << tD << " " << tE << endl;
+ x.push_back(tB);
+ y.push_back(tC);
+ dx.push_back(0.0);
+ dy.push_back(tE);
+ }
+ }
+ f.close();
+
+ TGraphErrors* g = new TGraphErrors(x.size(), &(x[0]), &(y[0]), &(dx[0]), &(dy[0]));
+ return g;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* GetSimDCS_ForFig(string filename){
+ vector<double> x, dx, y, dy;
+ ifstream f(filename.c_str());
+ string line;
+ while (getline(f, line)){
+ stringstream ss(line);
+ // i x y dx dy
+ double tA, tB, tC, tD, tE;
+ if (f >> tA >> tB >> tC >> tD >> tE){
+ //cout << tA << " " << tB << " " << tC << " " << tD << " " << tE << endl;
+ x.push_back(tB);
+ y.push_back(tC);
+ dx.push_back(tD);
+ dy.push_back(tE);
+ }
+ }
+ f.close();
+
+ TGraph* g = new TGraph(x.size(), &(x[0]), &(y[0]));
+ return g;
+}
+
+void CanvasPartition(TCanvas *C,const Int_t Nx = 2,const Int_t Ny = 2,
+ Float_t lMargin = 0.15, Float_t rMargin = 0.05,
+ Float_t bMargin = 0.15, Float_t tMargin = 0.05);
+double XtoPad(double x);
+double YtoPad(double x);
+
+void CanvasPartition(TCanvas *C,const Int_t Nx,const Int_t Ny, Float_t lMargin, Float_t rMargin, Float_t bMargin, Float_t tMargin){
+ if (!C) return;
+
+ // Setup Pad layout:
+ Float_t vSpacing = 0.0;
+ Float_t vStep = (1.- bMargin - tMargin - (Ny-1) * vSpacing) / Ny;
+
+ Float_t hSpacing = 0.0;
+ Float_t hStep = (1.- lMargin - rMargin - (Nx-1) * hSpacing) / Nx;
+
+ Float_t vposd,vposu,vmard,vmaru,vfactor;
+ Float_t hposl,hposr,hmarl,hmarr,hfactor;
+
+ for (Int_t i=0;i<Nx;i++) {
+
+ if (i==0) {
+ hposl = 0.0;
+ hposr = lMargin + hStep;
+ hfactor = hposr-hposl;
+ hmarl = lMargin / hfactor;
+ hmarr = 0.0;
+ } else if (i == Nx-1) {
+ hposl = hposr + hSpacing;
+ hposr = hposl + hStep + rMargin;
+ hfactor = hposr-hposl;
+ hmarl = 0.0;
+ hmarr = rMargin / (hposr-hposl);
+ } else {
+ hposl = hposr + hSpacing;
+ hposr = hposl + hStep;
+ hfactor = hposr-hposl;
+ hmarl = 0.0;
+ hmarr = 0.0;
+ }
+
+ for (Int_t j=0;j<Ny;j++) {
+
+ if (j==0) {
+ vposd = 0.0;
+ vposu = bMargin + vStep;
+ vfactor = vposu-vposd;
+ vmard = bMargin / vfactor;
+ vmaru = 0.0;
+ } else if (j == Ny-1) {
+ vposd = vposu + vSpacing;
+ vposu = vposd + vStep + tMargin;
+ vfactor = vposu-vposd;
+ vmard = 0.0;
+ vmaru = tMargin / (vposu-vposd);
+ } else {
+ vposd = vposu + vSpacing;
+ vposu = vposd + vStep;
+ vfactor = vposu-vposd;
+ vmard = 0.0;
+ vmaru = 0.0;
+ }
+
+ C->cd(0);
+
+ auto name = TString::Format("pad_%d_%d",i,j);
+ auto pad = (TPad*) C->FindObject(name.Data());
+ if (pad) delete pad;
+ pad = new TPad(name.Data(),"",hposl,vposd,hposr,vposu);
+ pad->SetLeftMargin(hmarl);
+ pad->SetRightMargin(hmarr);
+ pad->SetBottomMargin(vmard);
+ pad->SetTopMargin(vmaru);
+
+ pad->SetFrameBorderMode(0);
+ pad->SetBorderMode(0);
+ pad->SetBorderSize(0);
+
+ pad->Draw();
+ }
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+void FigureTest(){
+ gStyle->SetOptStat(0);
+
+ auto C = (TCanvas*) gROOT->FindObject("C");
+ if (C) delete C;
+ C = new TCanvas("C","canvas",700,2000);
+ C->SetFillStyle(4000);
+
+ // Number of PADS
+ const Int_t Nx = 2;
+ //const Int_t Ny = 5;
+ const Int_t Ny = 3;
+
+ // Margins
+ Float_t lMargin = 0.15;
+ Float_t rMargin = 0.01;
+ Float_t bMargin = 0.15;
+ Float_t tMargin = 0.01;
+
+ // Canvas setup
+ CanvasPartition(C,Nx,Ny,lMargin,rMargin,bMargin,tMargin);
+
+ // Dummy histogram.
+ string baseEx = "./CS2_TextFiles/02Apr_Exp_";
+ string baseTp = "./CS2_TextFiles/02Apr_TheoryP_";
+ string baseTf = "./CS2_TextFiles/02Apr_TheoryF_";
+ string baseTm = "./CS2_TextFiles/02Apr_TheoryMixed_";
+ //string baseTp = "./CS2_TextFiles/24Mar22_TheoryP_";
+ //string baseTf = "./CS2_TextFiles/24Mar22_TheoryF_";
+ //string baseTm = "./CS2_TextFiles/24Mar22_TheoryMixed_";
+
+ vector<string> states =
+ //{"3601", "3830", "2407", "2909", "1409", "1978", "0143", "0967"};
+ //{"0143", "3601"};
+ {"2909", "2909", "2909"};
+ vector<int> l =
+ //{ 3, 3, 1, 2, 1, 1, 1, 1};
+ //{1, 3};
+ {1, 3, 2};
+ vector<string> S =
+ //{ "0.24(1)", "0.34(1)", "0.33(2)", "", "0.24(1)", "0.50(1)", "0.48(1)", "0.30(2)"};
+ //{ "0.48(1)","0.24(1)"};
+ { "0.48(1)","0.24(1)","both"};
+
+
+
+
+
+
+ TPad *pad[Nx][Ny];
+
+ for (Int_t i = 0; i < Nx; i++) {
+ //for (Int_t j = 0; j < Ny; j++) {
+ for (Int_t j = 0; j < Ny-1; j++) {
+ C->cd(0);
+
+ // Get the pads previously created.
+ pad[i][j] = (TPad*) C->FindObject(TString::Format("pad_%d_%d",i,j).Data());
+ pad[i][j]->Draw();
+ pad[i][j]->SetFillStyle(4000);
+ pad[i][j]->SetFrameFillStyle(4000);
+ pad[i][j]->SetLogy();
+ pad[i][j]->cd();
+
+ // Size factors
+ Float_t xFactor = pad[0][0]->GetAbsWNDC()/pad[i][j]->GetAbsWNDC();
+ Float_t yFactor = pad[0][0]->GetAbsHNDC()/pad[i][j]->GetAbsHNDC();
+
+
+
+// int i,j;
+// if(p%2==0){ i=0; }
+// else { i=1; }
+// j=(p-i)/2;
+
+ int p = (2*j)+i;
+
+ cout << p << " " << i << " " << j << endl;
+
+ vector<double> xTp, xTf, yTp, yTf;
+
+ string sEx = baseEx + states[p] + ".txt";
+ string sTp = baseTp + states[p] + ".txt";
+ string sTf = baseTf + states[p] + ".txt";
+ string sTm = baseTm + states[p] + ".txt";
+
+
+ cout << sEx << endl;
+ TGraphErrors* gEx = GetExpDCS_ForFig(sEx);
+ gEx->SetMarkerStyle(7);
+ gEx->SetTitle("");
+
+ TGraphErrors* gTp = GetExpDCS_ForFig(sTp);
+ gTp->SetLineColor(kRed);
+ gTp->SetTitle("");
+// gTp->GetXaxis()->SetRangeUser(102.,158.);
+// gTp->GetYaxis()->SetRangeUser(0.000008,0.015);
+// gTp->SetTitle(" ");
+// gTp->GetXaxis()->SetTitleSize(0.08);
+// gTp->GetXaxis()->SetLabelSize(0.05);
+// gTp->GetYaxis()->SetTitleSize(0.08);
+// gTp->GetYaxis()->SetLabelSize(0.05);
+// gTp->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+// gTp->GetXaxis()->CenterTitle(true);
+// gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+// gTp->GetYaxis()->CenterTitle(true);
+
+ TGraphErrors* gTf = GetExpDCS_ForFig(sTf);
+ gTf->SetLineColor(kBlue);
+ gTf->SetTitle("");
+
+ TGraphErrors* gTm;
+
+ gTp->SetLineWidth(2);
+ gTf->SetLineWidth(2);
+ if(l[p] == 1){
+ gTf->SetLineStyle(7);
+ } else if (l[p] == 3){
+ gTp->SetLineStyle(7);
+ } else if (l[p] == 2){
+ gTp->SetLineStyle(7);
+ gTf->SetLineStyle(7);
+
+ gTm = GetExpDCS_ForFig(sTm);
+ gTm->SetLineColor(kViolet);
+ gTm->SetLineWidth(2);
+ gTm->SetTitle("");
+
+ }
+
+
+
+
+ // y axis range
+ //gTp->SetMinimum(0.0001); // do not show 0
+ //gTp->SetMaximum(1.2*h->GetMaximum());
+
+ // Format for y axis
+ gTp->GetYaxis()->SetLabelFont(43);
+ gTp->GetYaxis()->SetLabelSize(16);
+ gTp->GetYaxis()->SetLabelOffset(0.02);
+ gTp->GetYaxis()->SetTitle("");
+ gTp->GetYaxis()->SetTitleFont(43);
+ gTp->GetYaxis()->SetTitleSize(25);
+ gTp->GetYaxis()->SetTitleOffset(6);
+ gTp->GetYaxis()->CenterTitle();
+ //gTp->GetYaxis()->SetRangeUser(0.000006,0.015);
+ //gTp->GetYaxis()->SetRangeUser(0.00009,0.006);
+ gTp->GetYaxis()->SetRangeUser(0.00002,0.0011);
+ gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+
+ // TICKS Y Axis
+ //gTp->GetYaxis()->SetTickLength(xFactor*0.04/yFactor);
+
+ // Format for x axis
+ gTp->GetXaxis()->SetLabelFont(43);
+ gTp->GetXaxis()->SetLabelSize(16);
+ gTp->GetXaxis()->SetLabelOffset(0.02);
+ gTp->GetXaxis()->SetTitleFont(43);
+ gTp->GetXaxis()->SetTitleSize(25);
+ gTp->GetXaxis()->SetTitleOffset(5);
+ gTp->GetXaxis()->CenterTitle();
+ gTp->GetXaxis()->SetRangeUser(102.,158.);
+ gTp->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+
+ // TICKS X Axis
+ //gTp->GetXaxis()->SetTickLength(yFactor*0.06/xFactor);
+
+ // Draw cloned histogram with individual settings
+
+
+
+
+
+
+
+
+
+
+
+
+
+ gTp->Draw("ac");
+ gTf->Draw("same c");
+ if(l[p]==2){gTm->Draw("same c");}
+ gEx->Draw("same P");
+
+
+ double e = (double)stoi(states[p])/1000.;
+
+
+ cout << e << endl;
+
+
+// string titletext = to_string(e) + " MeV, S = "+-"
+// + to_string(dS[p]);
+
+ std::ostringstream out;
+
+ TText text;
+ TText text2;
+ if(e==3.83){
+ text.SetTextAlign(21);
+ text.SetTextFont(43);
+ text.SetTextSize(20);
+ text2.SetTextAlign(21);
+ text2.SetTextFont(43);
+ text2.SetTextSize(20);
+
+ text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "3.792 MeV, S = 0.18(1)");
+ text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "3.868 MeV, S = 0.16(1)");
+
+ } else if (e==2.909){
+
+ text.SetTextAlign(21);
+ text.SetTextFont(43);
+ text.SetTextSize(20);
+ text2.SetTextAlign(21);
+ text2.SetTextFont(43);
+ text2.SetTextSize(20);
+
+ text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "2.908 MeV");
+ text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "L=1 S=0.02(1), L=3 S=0.07(3)");
+
+ } else {
+ out.precision(3);
+ out << std::fixed << e << " MeV, S = " + S[p];
+// out.precision(2);
+// out << S[p] << " +- "
+// << dS[p];
+
+ text.SetTextAlign(21);
+ text.SetTextFont(43);
+ text.SetTextSize(20);
+ text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), (out.str()).c_str());
+
+ }
+
+//
+//
+// TH1F *hFrame = (TH1F*) h->Clone(TString::Format("h_%d_%d",i,j).Data());
+//
+// // y axis range
+// hFrame->SetMinimum(0.0001); // do not show 0
+// hFrame->SetMaximum(1.2*h->GetMaximum());
+//
+// // Format for y axis
+// hFrame->GetYaxis()->SetLabelFont(43);
+// hFrame->GetYaxis()->SetLabelSize(16);
+// hFrame->GetYaxis()->SetLabelOffset(0.02);
+// hFrame->GetYaxis()->SetTitleFont(43);
+// hFrame->GetYaxis()->SetTitleSize(16);
+// hFrame->GetYaxis()->SetTitleOffset(2);
+//
+// hFrame->GetYaxis()->CenterTitle();
+// hFrame->GetYaxis()->SetNdivisions(505);
+//
+// // TICKS Y Axis
+// hFrame->GetYaxis()->SetTickLength(xFactor*0.04/yFactor);
+//
+// // Format for x axis
+// hFrame->GetXaxis()->SetLabelFont(43);
+// hFrame->GetXaxis()->SetLabelSize(16);
+// hFrame->GetXaxis()->SetLabelOffset(0.02);
+// hFrame->GetXaxis()->SetTitleFont(43);
+// hFrame->GetXaxis()->SetTitleSize(16);
+// hFrame->GetXaxis()->SetTitleOffset(1);
+// hFrame->GetXaxis()->CenterTitle();
+// hFrame->GetXaxis()->SetNdivisions(505);
+//
+// // TICKS X Axis
+// hFrame->GetXaxis()->SetTickLength(yFactor*0.06/xFactor);
+//
+// // Draw cloned histogram with individual settings
+// hFrame->Draw();
+//
+// TText text;
+// text.SetTextAlign(31);
+// text.SetTextFont(43);
+// text.SetTextSize(10);
+// text.DrawTextNDC(XtoPad(0.9), YtoPad(0.8), gPad->GetName());
+ }
+ }
+ C->cd();
+
+ C->SaveAs("FigureCS2s.pdf");
+ C->SaveAs("FigureCS2s.svg");
+}
+
+void FigureTest_dt(){
+ gStyle->SetOptStat(0);
+
+ auto C = (TCanvas*) gROOT->FindObject("C");
+ if (C) delete C;
+ C = new TCanvas("C","canvas",700,2000);
+ C->SetFillStyle(4000);
+
+ // Number of PADS
+ const Int_t Nx = 2;
+ //const Int_t Ny = 5;
+ const Int_t Ny = 3;
+
+ // Margins
+ Float_t lMargin = 0.15;
+ Float_t rMargin = 0.01;
+ Float_t bMargin = 0.15;
+ Float_t tMargin = 0.01;
+
+ // Canvas setup
+ CanvasPartition(C,Nx,Ny,lMargin,rMargin,bMargin,tMargin);
+
+ // Dummy histogram.
+ string baseEx = "./Exp_dt";
+ string baseTp = "./ThryS_dt";
+ string baseTf = "./ThryD_dt";
+ string baseTm = "./ThryF_dt";
+ //string baseTp = "./CS2_TextFiles/24Mar22_TheoryP_";
+ //string baseTf = "./CS2_TextFiles/24Mar22_TheoryF_";
+ //string baseTm = "./CS2_TextFiles/24Mar22_TheoryMixed_";
+
+ vector<string> states =
+ {"3344", "3344", "3344"};
+ vector<int> l =
+ //{ 3, 3, 1, 2, 1, 1, 1, 1};
+ //{1, 3};
+ {2, 3, 1};
+ vector<string> S =
+ //{ "0.24(1)", "0.34(1)", "0.33(2)", "", "0.24(1)", "0.50(1)", "0.48(1)", "0.30(2)"};
+ //{ "0.48(1)","0.24(1)"};
+ { " "," "," "};
+
+
+
+
+
+
+ TPad *pad[Nx][Ny];
+
+ for (Int_t i = 0; i < Nx; i++) {
+ //for (Int_t j = 0; j < Ny; j++) {
+ for (Int_t j = 0; j < Ny-1; j++) {
+ C->cd(0);
+
+ // Get the pads previously created.
+ pad[i][j] = (TPad*) C->FindObject(TString::Format("pad_%d_%d",i,j).Data());
+ pad[i][j]->Draw();
+ pad[i][j]->SetFillStyle(4000);
+ pad[i][j]->SetFrameFillStyle(4000);
+ pad[i][j]->SetLogy();
+ pad[i][j]->cd();
+
+ // Size factors
+ Float_t xFactor = pad[0][0]->GetAbsWNDC()/pad[i][j]->GetAbsWNDC();
+ Float_t yFactor = pad[0][0]->GetAbsHNDC()/pad[i][j]->GetAbsHNDC();
+
+
+
+// int i,j;
+// if(p%2==0){ i=0; }
+// else { i=1; }
+// j=(p-i)/2;
+
+ int p = (2*j)+i;
+
+ cout << p << " " << i << " " << j << endl;
+
+ vector<double> xTp, xTf, yTp, yTf;
+
+ string sEx = baseEx + states[p] + ".txt";
+ string sTp = baseTp + states[p] + ".txt";
+ string sTf = baseTf + states[p] + ".txt";
+ string sTm = baseTm + states[p] + ".txt";
+
+
+ cout << sEx << endl;
+ TGraphErrors* gEx = GetExpDCS_ForFig(sEx);
+ gEx->SetMarkerStyle(7);
+ gEx->SetTitle("");
+
+ cout << sTp << endl;
+ TGraphErrors* gTp = GetExpDCS_ForFig(sTp);
+ gTp->SetLineColor(kRed);
+ gTp->SetTitle("");
+// gTp->GetXaxis()->SetRangeUser(102.,158.);
+// gTp->GetYaxis()->SetRangeUser(0.000008,0.015);
+// gTp->SetTitle(" ");
+// gTp->GetXaxis()->SetTitleSize(0.08);
+// gTp->GetXaxis()->SetLabelSize(0.05);
+// gTp->GetYaxis()->SetTitleSize(0.08);
+// gTp->GetYaxis()->SetLabelSize(0.05);
+// gTp->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+// gTp->GetXaxis()->CenterTitle(true);
+// gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+// gTp->GetYaxis()->CenterTitle(true);
+
+ cout << sTf << endl;
+ TGraphErrors* gTf = GetExpDCS_ForFig(sTf);
+ gTf->SetLineColor(kBlue);
+ gTf->SetTitle("");
+
+ TGraphErrors* gTm;
+
+ gTp->SetLineWidth(2);
+ gTf->SetLineWidth(2);
+ if(l[p] == 1){
+// gTf->SetLineStyle(7);
+ } else if (l[p] == 3){
+// gTp->SetLineStyle(7);
+ } else if (l[p] == 2){
+// gTp->SetLineStyle(7);
+// gTf->SetLineStyle(7);
+
+ cout << sTm << endl;
+ gTm = GetExpDCS_ForFig(sTm);
+ gTm->SetLineColor(kViolet);
+ gTm->SetLineWidth(2);
+ gTm->SetTitle("");
+
+ }
+
+
+
+
+ // y axis range
+ //gTp->SetMinimum(0.0001); // do not show 0
+ //gTp->SetMaximum(1.2*h->GetMaximum());
+
+ // Format for y axis
+ gTp->GetYaxis()->SetLabelFont(43);
+ gTp->GetYaxis()->SetLabelSize(16);
+ gTp->GetYaxis()->SetLabelOffset(0.02);
+ gTp->GetYaxis()->SetTitle("");
+ gTp->GetYaxis()->SetTitleFont(43);
+ gTp->GetYaxis()->SetTitleSize(25);
+ gTp->GetYaxis()->SetTitleOffset(6);
+ gTp->GetYaxis()->CenterTitle();
+ //gTp->GetYaxis()->SetRangeUser(0.000006,0.015);
+ //gTp->GetYaxis()->SetRangeUser(0.00009,0.006);
+ //gTp->GetYaxis()->SetRangeUser(0.00002,0.0011);
+ gTp->GetYaxis()->SetRangeUser(0.0002,0.02);
+ gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+
+ // TICKS Y Axis
+ //gTp->GetYaxis()->SetTickLength(xFactor*0.04/yFactor);
+
+ // Format for x axis
+ gTp->GetXaxis()->SetLabelFont(43);
+ gTp->GetXaxis()->SetLabelSize(16);
+ gTp->GetXaxis()->SetLabelOffset(0.02);
+ gTp->GetXaxis()->SetTitleFont(43);
+ gTp->GetXaxis()->SetTitleSize(25);
+ gTp->GetXaxis()->SetTitleOffset(5);
+ gTp->GetXaxis()->CenterTitle();
+ //gTp->GetXaxis()->SetRangeUser(102.,158.);
+ gTp->GetXaxis()->SetRangeUser(3.,15.);
+ gTp->GetXaxis()->SetTitle("#theta_{CM} [deg]");
+
+ // TICKS X Axis
+ //gTp->GetXaxis()->SetTickLength(yFactor*0.06/xFactor);
+
+ // Draw cloned histogram with individual settings
+
+
+
+
+
+
+
+
+
+
+
+
+
+ gTp->Draw("ac");
+ gTf->Draw("same c");
+ if(l[p]==2){gTm->Draw("same c");}
+ gEx->Draw("same P");
+
+
+ double e = (double)stoi(states[p])/1000.;
+
+
+ cout << e << endl;
+
+
+// string titletext = to_string(e) + " MeV, S = "+-"
+// + to_string(dS[p]);
+
+ std::ostringstream out;
+
+ TText text;
+ TText text2;
+ if(e==3.83){
+ text.SetTextAlign(21);
+ text.SetTextFont(43);
+ text.SetTextSize(20);
+ text2.SetTextAlign(21);
+ text2.SetTextFont(43);
+ text2.SetTextSize(20);
+
+ text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "3.792 MeV, S = 0.18(1)");
+ text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "3.868 MeV, S = 0.16(1)");
+
+ } else if (e==2.909){
+
+ text.SetTextAlign(21);
+ text.SetTextFont(43);
+ text.SetTextSize(20);
+ text2.SetTextAlign(21);
+ text2.SetTextFont(43);
+ text2.SetTextSize(20);
+
+ text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), "2.908 MeV");
+ text2.DrawTextNDC(XtoPad(0.5), YtoPad(0.78), "L=1 S=0.02(1), L=3 S=0.07(3)");
+
+ } else {
+ out.precision(3);
+ out << std::fixed << e << " MeV, S = " + S[p];
+// out.precision(2);
+// out << S[p] << " +- "
+// << dS[p];
+
+ text.SetTextAlign(21);
+ text.SetTextFont(43);
+ text.SetTextSize(20);
+ text.DrawTextNDC(XtoPad(0.5), YtoPad(0.88), (out.str()).c_str());
+
+ }
+
+//
+//
+// TH1F *hFrame = (TH1F*) h->Clone(TString::Format("h_%d_%d",i,j).Data());
+//
+// // y axis range
+// hFrame->SetMinimum(0.0001); // do not show 0
+// hFrame->SetMaximum(1.2*h->GetMaximum());
+//
+// // Format for y axis
+// hFrame->GetYaxis()->SetLabelFont(43);
+// hFrame->GetYaxis()->SetLabelSize(16);
+// hFrame->GetYaxis()->SetLabelOffset(0.02);
+// hFrame->GetYaxis()->SetTitleFont(43);
+// hFrame->GetYaxis()->SetTitleSize(16);
+// hFrame->GetYaxis()->SetTitleOffset(2);
+//
+// hFrame->GetYaxis()->CenterTitle();
+// hFrame->GetYaxis()->SetNdivisions(505);
+//
+// // TICKS Y Axis
+// hFrame->GetYaxis()->SetTickLength(xFactor*0.04/yFactor);
+//
+// // Format for x axis
+// hFrame->GetXaxis()->SetLabelFont(43);
+// hFrame->GetXaxis()->SetLabelSize(16);
+// hFrame->GetXaxis()->SetLabelOffset(0.02);
+// hFrame->GetXaxis()->SetTitleFont(43);
+// hFrame->GetXaxis()->SetTitleSize(16);
+// hFrame->GetXaxis()->SetTitleOffset(1);
+// hFrame->GetXaxis()->CenterTitle();
+// hFrame->GetXaxis()->SetNdivisions(505);
+//
+// // TICKS X Axis
+// hFrame->GetXaxis()->SetTickLength(yFactor*0.06/xFactor);
+//
+// // Draw cloned histogram with individual settings
+// hFrame->Draw();
+//
+// TText text;
+// text.SetTextAlign(31);
+// text.SetTextFont(43);
+// text.SetTextSize(10);
+// text.DrawTextNDC(XtoPad(0.9), YtoPad(0.8), gPad->GetName());
+ }
+ }
+ C->cd();
+
+// C->SaveAs("FigureCS2s.pdf");
+// C->SaveAs("FigureCS2s.svg");
+}
+
+double XtoPad(double x){
+ double xl,yl,xu,yu;
+ gPad->GetPadPar(xl,yl,xu,yu);
+ double pw = xu-xl;
+ double lm = gPad->GetLeftMargin();
+ double rm = gPad->GetRightMargin();
+ double fw = pw-pw*lm-pw*rm;
+ return (x*fw+pw*lm)/pw;
+}
+
+double YtoPad(double y){
+ double xl,yl,xu,yu;
+ gPad->GetPadPar(xl,yl,xu,yu);
+ double ph = yu-yl;
+ double tm = gPad->GetTopMargin();
+ double bm = gPad->GetBottomMargin();
+ double fh = ph-ph*bm-ph*tm;
+ return (y*fh+bm*ph)/ph;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Figure_AllCSFigures(){
+ string baseEx = "./CS2_TextFiles/24Mar22_Exp_";
+ string baseTp = "./CS2_TextFiles/24Mar22_TheoryP_";
+ string baseTf = "./CS2_TextFiles/24Mar22_TheoryF_";
+ string baseTm = "./CS2_TextFiles/24Mar22_TheoryMixed_";
+
+ vector<string> states =
+ {"0143", "0967", "1409", "1978", "2407", "2909", "3601", "3830"};
+ vector<int> l =
+ { 1, 1, 1, 1, 1, 2, 3, 3};
+
+ TCanvas *c = new TCanvas("c","c",700,1000);
+
+ CanvasPartition(c, 2, 4, 0.15,0.10,0.15,0.10);
+ TPad *pad[2][4];
+
+// c->SetBottomMargin(0.2);
+// c->Divide(2,4,0,0);
+// c->cd();
+//
+// vector<TPad*> pads;
+// TPad * p1 = new TPad("p1","p1", 0.1, 0.7, 0.5, 1.0);
+// TPad * p2 = new TPad("p2","p2", 0.5, 0.7, 1.0, 1.0);
+// TPad * p3 = new TPad("p3","p3", 0.1, 0.4, 0.5, 0.7);
+// TPad * p4 = new TPad("p4","p4", 0.5, 0.4, 1.0, 0.7);
+// TPad * p5 = new TPad("p5","p5", 0.1, 0.1, 0.5, 0.4);
+// TPad * p6 = new TPad("p6","p6", 0.5, 0.1, 1.0, 0.4);
+// pads.push_back(p1);
+// pads.push_back(p2);
+// pads.push_back(p3);
+// pads.push_back(p4);
+// pads.push_back(p5);
+// pads.push_back(p6);
+//
+// for (int i = 0; i < states.size(); i++){
+// }
+
+
+ for (int p = 0; p < states.size(); p++){
+ //for (int ii = 0; i < 1; i++){
+
+ int i,j;
+ if(p%2==0){ i=0; }
+ else { i=1; }
+ j=(p-i)/2;
+
+ cout << p << " " << i << " " << j << endl;
+
+ pad[i][j] = (TPad*) c->FindObject(TString::Format("pad_%d_%d",i,j).Data());
+ pad[i][j]->Draw();
+ pad[i][j]->cd();
+
+ vector<double> xTp, xTf, yTp, yTf;
+
+ string sEx = baseEx + states[p] + ".txt";
+ string sTp = baseTp + states[p] + ".txt";
+ string sTf = baseTf + states[p] + ".txt";
+ string sTm = baseTm + states[p] + ".txt";
+
+
+ cout << sEx << endl;
+ TGraphErrors* gEx = GetExpDCS_ForFig(sEx);
+ gEx->SetMarkerStyle(7);
+ gEx->GetYaxis()->SetTitleSize(0.08);
+ gEx->GetYaxis()->SetLabelSize(0.05);
+
+ TGraphErrors* gTp = GetExpDCS_ForFig(sTp);
+ gTp->SetLineColor(kRed);
+ gTp->GetXaxis()->SetRangeUser(102.,158.);
+ gTp->GetYaxis()->SetRangeUser(0.000008,0.015);
+ gTp->SetTitle(" ");
+ gTp->GetXaxis()->SetTitleSize(0.08);
+ gTp->GetXaxis()->SetLabelSize(0.05);
+ gTp->GetYaxis()->SetTitleSize(0.08);
+ gTp->GetYaxis()->SetLabelSize(0.05);
+ gTp->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+ gTp->GetXaxis()->CenterTitle(true);
+ gTp->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+ gTp->GetYaxis()->CenterTitle(true);
+
+ TGraphErrors* gTf = GetExpDCS_ForFig(sTf);
+ gTf->SetLineColor(kBlue);
+
+ TGraphErrors* gTm;
+
+ gTp->SetLineWidth(2);
+ gTf->SetLineWidth(2);
+ if(l[p] == 1){
+ gTf->SetLineStyle(7);
+ } else if (l[p] == 3){
+ gTp->SetLineStyle(7);
+ } else if (l[p] == 2){
+ gTp->SetLineStyle(7);
+ gTf->SetLineStyle(7);
+
+ gTm = GetExpDCS_ForFig(sTm);
+ gTm->SetLineColor(kMagenta);
+ gTm->SetLineWidth(2);
+
+ }
+
+
+
+
+ c->cd(p+1);
+ c->cd(p+1)->SetLogy();
+ c->cd(p+1)->SetFixedAspectRatio();
+// pads[i]->Draw();
+// pads[i]->SetLogy();
+
+ if(p%2==0){
+ c->cd(p+1)->SetLeftMargin( 0.2);
+ c->cd(p+1)->SetRightMargin( 0.0);
+ } else {
+ c->cd(p+1)->SetLeftMargin( 0.0);
+ c->cd(p+1)->SetRightMargin( 0.1);
+ }
+
+// switch(i){
+// case 0:
+// p1->cd();
+// break;
+// case 1:
+// p2->cd();
+// break;
+// case 2:
+// p3->cd();
+// break;
+// case 3:
+// p4->cd();
+// break;
+// }
+
+ gTp->Draw("ac");
+ gTf->Draw("same c");
+ if(l[p]==2){gTm->Draw("same c");}
+ gEx->Draw("same P");
+ }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void WriteToFile(TGraph* gr, string filename){
+ int npoints = gr->GetN();
+ ofstream file;
+ file.open (filename.c_str());
+
+ cout << " ============================ " << endl;
+ cout << " Writing to " << filename << " ..." << endl;
+
+ file << " i \t x \t y \t xe \t ye" << endl;
+
+ for(int i = 0; i < npoints; i++){
+ double x = gr->GetPointX(i);
+ double y = gr->GetPointY(i);
+ double xe = gr->GetErrorX(i);
+ double ye = gr->GetErrorY(i);
+
+ file << i << "\t"
+ << x << "\t"
+ << y << "\t"
+ << xe << "\t"
+ << ye << "\t"
+ << endl;
+
+ }
+
+ cout << " Complete! " << endl;
+ cout << " ============================ " << endl;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void canclone(TCanvas* major, int padNum, string name){
+ string minorName = "./CS2_Figures/" + name + ".root";
+ TFile* minorFile = new TFile(minorName.c_str(),"READ");
+ TCanvas* minor = (TCanvas*) minorFile->FindObjectAny("c_peakFits");
+
+ major->cd(padNum);
+ TPad *pad = (TPad*)gPad;
+ minor->cd();
+ TObject *obj, *clone;
+ pad->Range(minor->GetX1(),minor->GetY1(),minor->GetX2(),minor->GetY2());
+ pad->SetTickx(minor->GetTickx());
+ pad->SetTicky(minor->GetTicky());
+ pad->SetGridx(minor->GetGridx());
+ pad->SetGridy(minor->GetGridy());
+ pad->SetLogx(minor->GetLogx());
+ pad->SetLogy(minor->GetLogy());
+ pad->SetLogz(minor->GetLogz());
+ pad->SetBorderSize(minor->GetBorderSize());
+ pad->SetBorderMode(minor->GetBorderMode());
+ minor->TAttLine::Copy((TAttLine&)*pad);
+ minor->TAttFill::Copy((TAttFill&)*pad);
+ minor->TAttPad::Copy((TAttPad&)*pad);
+
+ TIter next(minor->GetListOfPrimitives());
+ gROOT->SetSelectedPad(pad);
+ while ((obj=next())) {
+ clone = obj->Clone();
+ obj->Draw("SAME");
+ pad->GetListOfPrimitives()->Add(clone,obj->GetDrawOption());
+ }
+ pad->Modified();
+ pad->Update();
+ major->cd(padNum);
+ pad->Draw();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double GetNodes(double spdf){
+ double nodes;
+ if(spdf==1 || spdf==0){ nodes=1; }
+ else if(spdf==3 || spdf==2){ nodes=0; }
+ else{
+ cout << " INPUT NODES::" << endl;
+ cin >> nodes;
+ }
+ return nodes;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void GenerateResidual(TGraphErrors* gdSdO, TGraph* Final){
+
+ TGraphErrors* gResid = new TGraphErrors(*gdSdO);
+ for(int n=0; n < gResid->GetN(); n++){
+ double x = gdSdO->GetPointX(n);
+ double residual = gdSdO->GetPointY(n) - Final->Eval(x);
+ gResid->SetPoint(n,x,residual);
+ gResid->SetPointError(n,0,gdSdO->GetErrorY(n));
+ }
+ TLine* markzero = new TLine(CSangleL,0.,CSangleH,0.);
+ gResid->SetTitle("");
+ gResid->GetXaxis()->SetRangeUser(CSangleL,CSangleH);
+ gResid->GetXaxis()->SetNdivisions(512,kTRUE);
+ gResid->GetYaxis()->SetTitle("Residuals");
+ gResid->GetYaxis()->SetTitleSize(0.15);
+ gResid->GetYaxis()->SetTitleOffset(0.36);
+ gResid->GetYaxis()->SetLabelSize(0.08);
+ gResid->GetYaxis()->SetNdivisions(305);
+ gResid->GetXaxis()->SetTitleSize(0.15);
+ gResid->GetXaxis()->SetTitleOffset(0.8);
+ gResid->GetXaxis()->SetLabelSize(0.1);
+ gResid->GetXaxis()->SetTickLength(0.1);
+ gResid->Draw();
+ markzero->SetLineStyle(2);
+ markzero->Draw("same");
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* AddTGraphs(TGraph* a, TGraph* b){
+ if(a->GetN() != b->GetN()){
+ cout << "CANNOT ADD GRAPHS, UNEQUAL NUMBER OF POINTS" << endl;
+ cout << a->GetN() << " != " << b->GetN() << endl;
+ return 0;
+ }
+
+ TGraph* sum = (TGraph*) a->Clone();
+
+ int maxN = a->GetN();
+ for(int i = 0 ; i < maxN ; i++){
+ double xa = a->GetPointX(i);
+ double xb = b->GetPointX(i);
+ if(xa!=xb){cout << "FAIL! Unequal X!" << endl; return 0;}
+
+ double ya = a->GetPointY(i);
+ double yb = b->GetPointY(i);
+ sum->SetPoint(i,xa,ya+yb);
+ }
+
+ return sum;
+
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void CS(double Energy, double Spin, double spdf, double angmom, double binning, string option){
+ /* Clean global variables, in case of multiple runs */
+ indexE = 100;
+ anglecentres.clear();
+ anglewidth.clear();
+ globalS=0.;
+ globalSerr=0.;
+ globalS2=0.;
+ globalSerr2=0.;
+ peakFitList->Clear();
+ numAngleBins=numAngleBinsBase;
+ globalBinning = binning;
+
+ //gROOT->SetBatch(1);
+
+ double spdf2, angmom2 = 0;
+ doMix = false;
+ doDoublet = false;
+
+ if(option=="mixed"){
+ doMix = true;
+ cout << BOLDBLUE
+ << " || MIXING OPTION SELECTED || \n"
+ << " || GIVE DETAILS OF 2ND OPT || \n"
+ << " || L = ... || \n"
+ << endl;
+ cin >> spdf2;
+ cout << " || J = ... || \n"
+ << endl;
+ cin >> angmom2;
+ cout << RESET;
+ } else if (option=="doublet") {
+ doDoublet = true;
+ }
+
+ if(binning < 0.01){
+ cout << BOLDRED
+ << " BINNING IS BELOW 0.01 MeV! Exiting..."
+ << RESET
+ << endl;
+ return;
+ }
+
+ /* Assign local variables */
+ // p3/2 -> spdf = 1, angmom = 1.5
+ // J0 is incident spin, which is 47K g.s. therefore J0 = 1/2
+ double J0 = 0.5;
+ double ElasticNorm = 0.000234, ElasticNormErr = 0.000; //18Oct22
+ inputdate = "18Oct22";
+
+ string backupFileName = "SolidAngle_HistFiles/";
+ string saFileName = "SolidAngle_HistFiles/SAHF_";
+ saFileName.append(inputdate);
+
+ if(reactionName=="47K(d,p)"){
+ // Which angular bin set to pull from
+ //numAngleBins=20; numAngleBinsBase=20; widthAngleBins=2.5; firstAngle=105.0;
+ numAngleBins=26; numAngleBinsBase=26; widthAngleBins=2.0; firstAngle=104.0;
+
+ // Selecting (d,p) global values
+ numPeaks_CS = numPeaks;
+ means_CS=means;
+
+ // Additional...
+ backupFileName.append("SolidAngle_HistFile_10Aug22_TrueStripRemoval.root");
+ saFileName.append("_47Kdp_");
+ CSangleL = 100.;
+ CSangleH = 160.;
+ } else if(reactionName=="47K(d,t)"){
+ // Which angular bin set to pull from
+ numAngleBins=18; numAngleBinsBase=18; widthAngleBins=1.0; firstAngle=2.0;
+
+ // Selecting (d,p) global values
+ numPeaks_CS = numPeaks_dt;
+ means_CS=means_dt;
+
+ // Additional...
+ backupFileName.append("SolidAngle_HistFile_18Oct22_47Kdt.root");
+ saFileName.append("_47Kdt_");
+ CSangleL = 00.;
+ CSangleH = 30.;
+ }
+
+ // Extract spin orbit name
+ orbitalname.clear();
+ orbital.clear();
+
+ switch((int)spdf){
+ case 0:
+ orbitalname.append("s_{");
+ orbital.append("s");
+ break;
+ case 1:
+ orbitalname.append("p_{");
+ orbital.append("p");
+ break;
+ case 2:
+ orbitalname.append("d_{");
+ orbital.append("d");
+ break;
+ case 3:
+ orbitalname.append("f_{");
+ orbital.append("f");
+ break;
+ default:
+ cout << RED
+ <<"ERROR: SPDF SELECTION OUT OF RANGE"
+ << endl;
+ }
+ orbitalname.append(to_string((int)(2*angmom)));
+ orbitalname.append("/2}");
+ orbital.append(to_string((int)(2*angmom)));
+ orbital.append("2");
+
+ // Number of nodes
+ double nodes = GetNodes(spdf);
+ double nodes2;
+ if(doMix){
+ nodes2 = GetNodes(spdf2);
+ }
+
+ /* Retrieve array index of the entered peak energy */
+ /* numpeaks and Energy[] defined globally in KnownPeakFitter.h */
+ bool found = 0;
+ for(int i=0;i<numPeaks_CS;i++){
+ if(abs(Energy-means_CS.at(i))<0.01){
+ cout << "========================================================" << endl;
+ cout << "Identified as state #" << i << ", E = " << means_CS.at(i) << endl;
+ indexE = i;
+ found = 1;
+ stringstream ss;
+ ss << setfill('0') << setw(4) << (int)(means_CS.at(i)*1000.);
+ statename = ss.str();
+ }
+ }
+ if(!found){
+ cout << "========================================================" << endl;
+ cout << "NO STATE AT THAT ENERGY INDENTIFIED!! CHECK KNOWN PEAKS!!" << endl;
+ return;
+ }
+ saFileName.append(statename);
+ saFileName.append(".root");
+
+ /* Solid Angle (from simulation) */
+ TFile* file;
+ ifstream infile(saFileName);
+ ifstream backup(backupFileName);
+ if(infile.good()){
+ cout << GREEN << "Opening file " << saFileName << RESET << endl;
+ file = TFile::Open(saFileName.c_str());
+ } else if (backup.good()){
+ cout << BOLDRED << "FAILED TO OPEN " << saFileName << endl;
+ cout << "Open standard backup file " << backupFileName << endl;
+ cout << RESET;
+ file = TFile::Open(backupFileName.c_str());
+ } else {
+ cout << BOLDRED << "FAILED TO OPEN MAIN OR BACKUP SOLID ANGLE FILE" << endl;
+ cout << RED << "Check SolidAngle file exists..." << RESET << endl;
+ }
+
+ string objName;
+ if(reactionName=="47K(d,p)"){
+ objName.append("SolidAngle_Lab_MG");
+ }
+ else if(reactionName=="47K(d,t)"){
+ objName.append("SolidAngle_CM_MM");
+ }
+
+ TH1F* SolidAngle = (TH1F*) file->FindObjectAny(objName.c_str());
+ TCanvas* c_SolidAngle = new TCanvas("c_SolidAngle","c_SolidAngle",1000,1000);
+ SolidAngle->Draw("HIST");
+ SolidAngle->GetXaxis()->SetRangeUser(CSangleL,CSangleH);
+ /* (canvas deleted after Area/SA calculation) */
+
+ /* Area of experimental peaks */
+ TCanvas* c_PeakArea = new TCanvas("c_PeakArea","c_PeakArea",1000,1000);
+ vector<vector<double>> areaArray = GetExpDiffCross(means_CS.at(indexE));
+// delete c_PeakArea;
+
+ // Array: peakenergy, peakarea, areaerror, anglemin, anglemax
+ if(loud){
+ for(int i=0; i<areaArray.size();i++){
+ cout << i << " "
+ << areaArray[i][0] << " "
+ << areaArray[i][1] << " "
+ << areaArray[i][2] << " "
+ << areaArray[i][3] << " "
+ << areaArray[i][4] << endl;
+ }
+ }
+
+ /* AoSA = Area/Solid Angle [#/msr] */
+ /* dSdO = Experimental Diff. Cross Sect. (Area/Solid Angle *Norm) [mb/msr] */
+ vector<Double_t> AoSA, AoSAerr;
+ vector<Double_t> expDCS, expDCSerr;
+ for(int i=0; i<areaArray.size();i++){
+ int binmin = SolidAngle->FindBin(areaArray[i][3]+0.0001);
+ int binmax = SolidAngle->FindBin(areaArray[i][4]-0.0001);
+
+ 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
+ SAerr = SAerr*1000.; //sr->msr
+
+ /* Area over Solid angle ONLY */
+ AoSA.push_back(areaArray[i][1]/SA);
+ AoSAerr.push_back((areaArray[i][1]/SA)
+ * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
+ + pow(SAerr/SA,2)));
+
+ /* Differential Cross Section */
+ /* NOTE: DON'T INCLUDE NORM ERROR IN ERR PROPAGATION AS IT IS SYSTEMATIC! */
+ double tempvalue = (areaArray[i][1]/SA)*ElasticNorm;
+ double temperror = tempvalue
+ * sqrt( pow(areaArray[i][2]/areaArray[i][1],2)
+ + pow(SAerr/SA,2)
+ );
+ //if(temperror<tempvalue){ // exclude very large error bars
+ expDCS.push_back(tempvalue);
+ expDCSerr.push_back(temperror);
+ //}
+
+ if(loud){
+ cout << "Angle " << areaArray[i][3] << " to " << areaArray[i][4]
+ << ", centre " << anglecentres[i]
+ << ": Area = " << areaArray[i][1] << " +- " << areaArray[i][2] << " cnts"
+ << " SA = " << SA << " +- " << SAerr << " msr"
+ << " Area/SA = " << AoSA[i] << " +- " << AoSAerr[i] << " counts/msr"
+ << setprecision(5)
+ << " Norm = " << ElasticNorm << " +- " << ElasticNormErr
+ << " (don't include norm err in propagation)"
+ << " dSdO = " << tempvalue << " +- " << temperror
+ << setprecision(3)
+ << endl;
+ }
+
+ }
+ //delete c_SolidAngle;
+
+ //gROOT->SetBatch(0);
+
+ /* Graph of Area/Solid Angle*/
+ TCanvas* c_AoSA = new TCanvas("c_AoSA","c_AoSA",1000,1000);
+ c_AoSA->SetLogy();
+ TGraphErrors* gAoSA = new TGraphErrors(
+ anglecentres.size(), //n
+ &(anglecentres[0]), &(AoSA[0]), //x, y
+ 0, &(AoSAerr[0]) ); //errX, errY
+ gAoSA->SetTitle("Area/SolidAngle");
+ gAoSA->GetXaxis()->SetTitle("ThetaLab [deg]");
+ gAoSA->GetYaxis()->SetTitle("Counts/#Omega [counts/msr]");
+ gAoSA->Draw();
+
+ /* Graph of experimental diff. cross sect (dSigma/dOmega) */
+ TCanvas* c_dSdO = new TCanvas("c_dSdO","c_dSdO",1000,1000);
+ c_dSdO->SetLogy();
+ TGraphErrors* gdSdO = new TGraphErrors(
+ anglecentres.size(),
+ &(anglecentres[0]), &(expDCS[0]),
+ 0, &(expDCSerr[0]) );
+ gdSdO->SetTitle("Differential Cross Section");
+ if(reactionName=="47K(d,p)"){
+ gdSdO->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+ } else if (reactionName=="47K(d,t)"){
+ gdSdO->GetXaxis()->SetTitle("#theta_{CM} [deg]");
+ }
+ gdSdO->GetXaxis()->SetTitleOffset(1.2);
+ gdSdO->GetXaxis()->SetTitleSize(0.04);
+ gdSdO->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]");
+ gdSdO->GetYaxis()->SetTitleOffset(1.2);
+ gdSdO->GetYaxis()->SetTitleSize(0.04);
+ gdSdO->Draw();
+ c_dSdO->Update();
+
+ /* 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_CS.at(indexE), J0, Spin, nodes, spdf, angmom);
+ TheoryDiffCross->GetYaxis()->SetTitle("d#sigma/d#Omega [mb/msr]"); //msr set in func above
+ TheoryDiffCross->GetXaxis()->SetTitle("ThetaLab [deg]");
+ TheoryDiffCross->Draw();
+
+
+ TGraph* TheoryDiffCross2;
+ if(doMix){
+ TheoryDiffCross2 = TWOFNR(means_CS.at(indexE), J0, Spin, nodes2, spdf2, angmom2);
+ TheoryDiffCross2->SetLineColor(kBlue);
+ TheoryDiffCross2->Draw("same");
+ }
+
+ /** TEMP **/
+ if(loud){
+ cout << "UNSCALED THEORY DIFF CROSS SECTION EVALUATED AT DATA POINTS:::" << endl;
+ cout << setprecision(6);
+ for(int i=0; i<numAngleBins; i++){
+ cout << anglecentres.at(i) << "\t" << TheoryDiffCross->Eval(anglecentres.at(i)) << endl;
+ }
+ cout << setprecision(3);
+ cout << "......................" << endl;
+ }
+ /** **** **/
+
+ /* Using Chi2 minimizaiton */
+ if(loud){cout << "USING CHI2 MINIMIZAITON..." << endl;}
+ TCanvas* c_Chi2Min = new TCanvas("c_Chi2Min","c_Chi2Min",1000,1000);
+ gStyle->SetPadLeftMargin(0.12);
+ gStyle->SetPadRightMargin(0.03);
+
+ TPad *pad1 = new TPad("pad1","pad1",0,0.25,1,1);
+ TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.25);
+ pad1->SetTopMargin(0.1);
+ pad1->SetBottomMargin(0.00001);
+ pad1->SetBorderMode(0);
+ pad1->SetLogy();
+ pad1->SetTickx();
+ pad1->SetTicky();
+ pad2->SetTopMargin(0.00001);
+ pad2->SetBottomMargin(0.3);
+ pad2->SetBorderMode(0);
+ pad2->SetTickx();
+ pad2->SetTicky();
+ pad1->Draw();
+ pad2->Draw();
+ pad1->cd();
+
+ TGraphErrors* gdSdO2 = new TGraphErrors(*gdSdO);
+ gdSdO2->SetName("gdSdO2");
+ gdSdO2->SetTitle("gdSdO2");
+
+
+ TGraph* Final;
+
+ if(doMix){
+ Final = FindNormalisation(TheoryDiffCross,TheoryDiffCross2,gdSdO);
+ } else {
+ Final = FindNormalisation(TheoryDiffCross,gdSdO);
+ }
+
+ Final->SetName("Final");
+ Final->SetTitle("Final");
+
+ gdSdO->SetLineColor(kRed);
+ gdSdO->SetMarkerColor(kRed);
+ gdSdO->SetMarkerStyle(21);
+
+
+ /* Construct file name string */
+ /**/ ostringstream tempstream;
+ /**/ if(means_CS.at(indexE)<1.0){tempstream << 0;}
+ /**/ tempstream << (int) (means_CS.at(indexE)*1000);
+ /**/ tempstream << "_" << orbital;
+ /**/ tempstream << "_spin" << Spin;
+ /**/ string tempstr = tempstream.str();
+ /* Construct hist title string */
+ /**/ ostringstream textstream;
+ /**/ textstream << std::fixed << setprecision(3);
+ /**/ textstream << " " << means_CS.at(indexE);
+ /**/ textstream << " MeV, ";
+ /**/ textstream << orbitalname;
+ /**/ textstream << ", spin " << (int)Spin;
+ /**/ textstream << " --> S = " << globalS
+ /**/ << " +- " << globalSerr;
+ /**/ if(doMix){
+ /**/ textstream << " & " << globalS2
+ /**/ << " +- " << globalSerr2;
+ /**/ }
+ /**/ string textstring = textstream.str();
+
+ gdSdO->SetTitle(textstring.c_str());
+ gdSdO->GetYaxis()->SetTitleOffset(1.3);
+ gdSdO->GetYaxis()->SetTitleSize(0.042);
+ gdSdO->GetXaxis()->SetRangeUser(CSangleL,CSangleH);
+ //gdSdO->GetYaxis()->SetRangeUser(5e-6,5e-3); //Same vertical range for all? WNC suggestion
+ gdSdO->GetYaxis()->SetRangeUser(2e-6,5e-3); //Same vertical range for all? WNC suggestion
+ gdSdO->GetXaxis()->SetNdivisions(512,kTRUE);
+ gdSdO->Draw("AP");
+
+ Final->Draw("SAME");
+
+ pad2->cd();
+
+cout << "!!! !!! !!! !!! here1" << endl;
+
+ GenerateResidual(gdSdO,Final);
+
+cout << "!!! !!! !!! !!! here2" << endl;
+
+ string savestring1 = "./CS2_Figures/"+tempstr+".root";
+ string savestring2 = "./CS2_Figures/"+tempstr+".pdf";
+ c_Chi2Min->SaveAs(savestring1.c_str());
+ c_Chi2Min->SaveAs(savestring2.c_str());
+
+ cout << YELLOW
+ << " ----- C2S = (2J+1)S = (2*" << (int)Spin << "+1)S = ("
+ << (int)((2.*Spin)+1.) << ")S = "
+ << ((2.*Spin)+1.) * globalS << " ----- "
+ << RESET << endl;
+
+
+ //delete c_AoSA;
+ //delete c_dSdO;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+vector<vector<double>> GetExpDiffCross(double Energy){
+ cout << "========================================================" << endl;
+ vector<vector<double>> AllPeaks_OneGate;
+ vector<vector<double>> OnePeak_AllGates;
+ /****CHANGE ANGLE GATING****/
+ //double widthAngleBins = 2.5;
+ //double firstAngle = 105.;
+ /***************************/
+ double x[numAngleBins], y[numAngleBins];
+ //TList* list = new TList();
+
+ double trackScale = 0.0;
+ if(reactionName=="47K(d,p)"){
+ /* Determine scaling factor for PhaseSpace */
+ TCanvas* c_ExSubPSpace = new TCanvas("c_ExSubPSpace","c_ExSubPSpace",1000,1000);
+ if(scaleTogether){
+ TH1F* baseEx = PullThetaLabHist(0,firstAngle,widthAngleBins);
+ TH1F* basePS = PullPhaseSpaceHist(0,firstAngle,widthAngleBins);
+ for(int i=1; i<numAngleBins;i++){
+ TH1F* addEx = PullThetaLabHist(i,firstAngle,widthAngleBins); baseEx->Add(addEx,1.);
+ TH1F* addPS = PullPhaseSpaceHist(i,firstAngle,widthAngleBins); basePS->Add(addPS,1.);
+ }
+
+ /* Subtract flat background equal to smallest bin in range */
+ baseEx->GetXaxis()->SetRange(baseEx->FindBin(-1.),baseEx->FindBin(1.));
+ double minValueInRange = baseEx->GetBinContent(baseEx->GetMinimumBin());
+ baseEx->GetXaxis()->UnZoom();
+ cout << "Subtracting background of " << minValueInRange << endl;
+ for(int b=1; b<baseEx->GetNbinsX() ; b++){
+ baseEx->SetBinContent(b,baseEx->GetBinContent(b)-minValueInRange);
+ }
+
+ // Fix error where basePS = 600, baseEx=150
+ if(loud){
+ cout << "baseEx bins: " << baseEx->GetNbinsX() << endl;
+ cout << "basePS bins: " << basePS->GetNbinsX() << endl;
+ }
+ if(basePS->GetNbinsX()!=baseEx->GetNbinsX()){
+ cout << basePS->GetNbinsX() / baseEx->GetNbinsX() << endl;
+ int ratio = basePS->GetNbinsX() / baseEx->GetNbinsX();
+ basePS->Rebin(ratio);
+ }
+ if(loud){
+ cout << "baseEx bins: " << baseEx->GetNbinsX() << endl;
+ cout << "basePS bins: " << basePS->GetNbinsX() << endl;
+ }
+
+ /* Begin scaling within range, track changes */
+ basePS->Scale(0.1);
+ trackScale = 0.1;
+ int numAngleBinsScale = baseEx->GetNbinsX();
+ int nbinlow = basePS->FindBin(4.); int nbinhigh = basePS->FindBin(8.0);
+ for(int b=nbinlow; b<nbinhigh; b++){
+ if(baseEx->GetBinContent(b) > 0.0 && basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
+ while(basePS->GetBinContent(b) > baseEx->GetBinContent(b)){
+ basePS->Scale(0.99999);
+ trackScale *= 0.99999;
+ }
+ }
+ }
+
+ baseEx->Add(basePS,-1.);
+ baseEx->SetName("ExSubPSpace");
+ baseEx->SetTitle("ExSubPSpace");
+ baseEx->Draw();
+ if(loud){cout << "PhaseSpace -> ExpData scaling = " << trackScale << endl;}
+ }
+ }
+
+ /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
+ if(reactionName=="47K(d,p)"){
+ if(means_CS.at(indexE)>3.5){
+ //Find max angle, then round down to nearest integer
+ double thetaMax = -14.36*means_CS.at(indexE) + 205.50;
+ thetaMax = floor(thetaMax-5.0);//TEMP REMOVE MORE!!!!
+ if(loud){cout << RED << thetaMax << " rounds by " << widthAngleBins << " to ";}
+
+ //Round down to nearest angular bin (i.e. 2.5deg sections)
+ thetaMax = floor((double)thetaMax/widthAngleBins)*widthAngleBins;
+ if(loud){cout << thetaMax << RESET << endl;}
+
+ numAngleBins = (int)(thetaMax - firstAngle)/2.5;
+ }
+ }
+ else if(reactionName=="47K(d,t)"){
+ double thetaMax = 0.23*pow(means_dt[indexE],2)
+ + 0.76*means_dt[indexE]
+ + 9.01;
+ thetaMax = floor(thetaMax);
+
+ double thetaMin = 0.10*pow(means_dt[indexE],2)
+ + 0.06*means_dt[indexE]
+ + 1.58;
+ thetaMin = ceil(thetaMin);
+
+ numAngleBins = (int)thetaMax - (int)thetaMin;
+ firstAngle = (int)thetaMin;
+
+ cout << RED << "ANGLES:: " << firstAngle << " -> " << thetaMax << RESET << endl;
+
+ }
+
+ /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
+
+
+ //for(int i=0; i<1; i++){
+ for(int i=0; i<numAngleBins;i++){
+ double min = firstAngle + (i*widthAngleBins);
+ double max = min + widthAngleBins;
+ cout << "===================================" << endl;
+ cout << "min: " << min << " max: " << max << endl;
+
+ stringstream tmp; tmp << fixed << setprecision(0);
+ tmp << "c_peakFits_" << min << "_" << max;
+ string tmp2 = tmp.str();
+ TCanvas* c_peakFits = new TCanvas("c_peakFits",tmp2.c_str(),1000,1000);
+
+ /* Retrieve theta-gated Ex TH1F from file GateThetaLabHistograms.root */
+ /* To change angle gates, run GateThetaLab_MultiWrite() */
+ TH1F *gate, *pspace;
+ if(reactionName=="47K(d,p)"){
+ gate = PullThetaLabHist(i,firstAngle,widthAngleBins);
+ pspace = PullPhaseSpaceHist(i,firstAngle,widthAngleBins);
+
+ /* Scale the Phase Space at this angle... */
+ /* ... for all angles together */
+ if(scaleTogether){
+ if(loud){ cout << "gate bins: " << gate->GetNbinsX() << endl;
+ cout << "pspc bins: " << pspace->GetNbinsX() << endl; }
+ // Fix error where basePS = 600, baseEx=150
+ if(pspace->GetNbinsX()!=gate->GetNbinsX()){
+ if(loud){cout << pspace->GetNbinsX() / gate->GetNbinsX() << endl;}
+ int ratio = pspace->GetNbinsX() / gate->GetNbinsX();
+ pspace->Rebin(ratio);
+ }
+ if(loud){ cout << "gate bins: " << gate->GetNbinsX() << endl;
+ cout << "pspc bins: " << pspace->GetNbinsX() << endl; }
+ gate->Add(pspace,-trackScale);
+ }
+ /* ... or seperately for each angular bin */
+ /* NOTE THAT THIS DOES NOT ACCOUNT FOR FLAT BACKGROUND */
+ else {
+ if(pspace->Integral() > 50.){ // Non-garbage histogram
+ pspace->Scale(0.01);
+ trackScale=0.01;
+ int numAngleBins = gate->GetNbinsX();
+ for(int b=0; b<numAngleBins; b++){
+ if(loud){cout << " FROM " << pspace->GetBinContent(b) <<
+ " > " << gate->GetBinContent(b);
+ }
+ while(pspace->GetBinContent(b) > gate->GetBinContent(b)){
+ pspace->Scale(0.9999);
+ trackScale*=0.9999;
+ }
+ if(loud){cout << " TO " << pspace->GetBinContent(b) <<
+ " > " << gate->GetBinContent(b) << endl;
+ }
+ }
+ if(loud){cout << " !!! SCALE FOR THIS ANGLE = " << trackScale << 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, 0.0);
+ }
+ else if (reactionName=="47K(d,t)"){
+ gate = PullThetaCMHist(i,firstAngle,widthAngleBins);
+ //AllPeaks_OneGate = FitKnownPeaks_dt_RtrnArry(gate, 0.0);
+ AllPeaks_OneGate = FitKnownPeaks_RtrnArry(gate, 0.0);
+ }
+
+ /* Write PS-subtracted spectrum to list */
+ //list->Add(gate);
+ //list->Add(c_peakFits);
+ gate->GetXaxis()->SetRangeUser(-1.0,+8.0);
+ string filename = "./CS2_Figures/";
+ filename += tmp2 + ".root";
+ c_peakFits->SaveAs(filename.c_str());
+ //auto tempfile = new TFile(filename.c_str(),"UPDATE"); //Reopen newly made file
+ //auto templist = new TList();
+ //templist->Add();
+
+
+ /* Check correct OneGate vector is selected */
+ if(loud){
+ cout << "area of " << means_CS.at(indexE) << " = "
+ << AllPeaks_OneGate[indexE][1]
+ << " +- " << AllPeaks_OneGate[indexE][2]
+ << endl;
+ }
+
+ /* Add min and max angle to end of relevant OneGate vector */
+ AllPeaks_OneGate[indexE].push_back(min);
+ AllPeaks_OneGate[indexE].push_back(max);
+
+ /* Push relevant OneGate vector to end of AllGates vector */
+ OnePeak_AllGates.push_back(AllPeaks_OneGate[indexE]);
+ //delete c_peakFits;
+ }
+
+ /* 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){
+ //string name = "GateThetaLabHistograms_18Oct22_47Kdp.root";
+ string name = "GateThetaLabHistograms_18Oct22_47Kdp_26angles.root";
+ TFile* file = new TFile(name.c_str(),"READ");
+ string histname = "cThetaLabGate_"
+ + 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();
+
+ double pulledBin = hist->GetBinWidth(10);
+ int ratio = (int)(globalBinning / pulledBin);
+ hist->Rebin(ratio);
+
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TH1F* PullThetaCMHist(int i, double minTheta, double gatesize){
+ //TFile* file = new TFile("GateThetaCMHistograms_47Kdt_18Oct22_bin0p1.root","READ");
+ string name = "GateThetaCMHistograms_47Kdt_18Oct22_";
+ if(!GammaGate){
+ name = name + "bin0p1.root";
+ } else {
+ name = name + to_string((int)(globGammaGate*1000))
+ + "GammaGate.root";
+ }
+ TFile* file = new TFile(name.c_str(),"READ");
+
+ string histname = "cThetaCMGate_"
+ + to_string((int)(minTheta+(i*gatesize))) + "-"
+ + to_string((int)(minTheta+((i+1)*gatesize)));
+ cout << "Loading " << histname << endl;
+ TList *list = (TList*)file->Get("GateThetaCMHistograms");
+ TH1F* hist = (TH1F*)list->FindObject(histname.c_str());
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TH1F* PullPhaseSpaceHist(int i, double minTheta, double gatesize){
+ //string name = "GatePhaseSpaceThetaLabHistograms_18Oct22.root";
+ string name = "GatePhaseSpaceThetaLabHistograms_18Oct22_26angles.root";
+ TFile* file = new TFile(name.c_str(),"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();
+
+ double pulledBin = hist->GetBinWidth(10);
+ int ratio = (int)(globalBinning / pulledBin);
+ hist->Rebin(ratio);
+
+ return hist;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+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);
+ }
+ }
+
+ //scaleTWOFNR = mean;
+ if(loud){
+ cout << "SCALED THEORY BY " << mean << endl;
+ cout << " therefore S = " << 1/mean << " ??" << endl;
+ }
+
+ 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){
+ /* This function moves between directories in order to run TWOFNR in proper *
+ * location. This is, weirdly, the least tempremental way of doing this. */
+
+ cout << "========================================================" << endl;
+
+ double QValue, BeamEnergy = 7.7;
+ int johnson, tandyval, modelA, modelB;
+ string njjj;
+ if(reactionName=="47K(d,p)"){
+ cout << "Using Johnson-Soper ..."; johnson=5; tandyval=0;
+ cout << " ... and Chapel-Hill." << endl; modelA=2; modelB=2;
+ //QValue = 2.274 - E;
+ QValue = 2.419 - E;
+ njjj.append("24.jjj");
+ } else if (reactionName=="47K(d,t)"){
+ QValue = -2.112 - E;
+ johnson=1; //just makess it work below with the dp/dt matching
+ njjj.append("21.jjj");
+ }
+
+ char origDirchar[200];
+ getcwd(origDirchar,200);
+ string origDir{origDirchar};
+ string twofnrDir = "/home/charlottepaxman/Programs/TWOFNR";
+ cout << "Current directory " << origDir << endl;
+ cout << "Moving to directory " << twofnrDir << endl;
+ chdir(twofnrDir.c_str());
+ //Check
+ system("pwd");
+ cout << "===================================" << endl;
+
+ /* Delete existing tran.jjj & 24.jjj files */
+ remove("tran.jjj");
+ remove(njjj.c_str());
+
+ std::ofstream Front_Input("in.front");
+ Front_Input << "jjj" << std::endl;
+ Front_Input << "pipo" << std::endl;
+ if(reactionName=="47K(d,p)"){
+ Front_Input << 2 << std::endl;
+ } else if (reactionName=="47K(d,t)"){
+ Front_Input << 5 << std::endl;
+ }
+ Front_Input << 0 << std::endl;
+ Front_Input << 0 << std::endl;
+ Front_Input << BeamEnergy << std::endl;
+ Front_Input << 47 << " " << 19 << 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;
+ Front_Input << QValue << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << J0 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << johnson << std::endl;
+ if(johnson==6){//JTandy selected, give version
+ Front_Input << tandyval << std::endl;
+ }
+ Front_Input << 1 << std::endl;
+ Front_Input << J << std::endl;
+ if(reactionName=="47K(d,p)"){
+ Front_Input << 1 << std::endl;
+ Front_Input << modelA << std::endl;
+ Front_Input << modelB << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ Front_Input << 1 << std::endl;
+ } else if (reactionName=="47K(d,t)"){
+ Front_Input << 1 << 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();
+
+ cout << "Filled front20 input file." << endl;
+ cout << "Executing front20..." << endl;
+ system("(exec ~/Programs/TWOFNR/front20 < in.front > /dev/null)");
+ ifstream checkfront("tran.jjj");
+ if(!checkfront){
+ //front20 fail!
+ cout << " !! ERROR !! \n front20 failed to complete" << endl;
+ return 0;
+ } else {
+ cout << "-> tran.jjj generated (success not guaranteed)" << endl;
+ checkfront.close();
+ }
+
+ /* in.twofnr instructs twofnr20 to evaluate tran.jjj */
+ cout << "Executing twofnr20..." << endl;
+ system("(exec ~/Programs/TWOFNR/twofnr20 < in.twofnr > /dev/null)");
+ ifstream checktwofnr(njjj.c_str());
+ if(!checktwofnr){
+ //twofnr20 fail!
+ cout << " !! ERROR !! \n twofnr20 failed to complete" << endl;
+ terminate();
+// return 0;
+ } else {
+ cout << "-> twofnr20 complete!" << endl;
+ checktwofnr.close();
+ }
+
+ TGraph* CS = new TGraph(njjj.c_str());
+
+ //mb/sr->mb/msr is x1/1000
+ for(int i=0; i<CS->GetN(); i++){
+ double x, newy;
+ CS->GetPoint(i,x,newy);
+ newy=newy/1000;
+ CS->SetPoint(i,x,newy);
+ }
+
+
+ cout << "===================================" << endl;
+ cout << "Current directory " << twofnrDir << endl;
+ cout << "Moving to directory " << origDir << endl;
+ chdir(origDir.c_str());
+ system("pwd");
+ cout << "========================================================" << endl;
+
+ return CS;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double Chi2(TGraph* theory, TGraphErrors* exper){
+ double Chi2 = 0;
+ double chi = 0;
+
+ //cout << setprecision(8);
+ //for(int i = 1 ; i < exper->GetN() ; i++){
+ for(int i = 0 ; i < exper->GetN() ; i++){
+ if(exper->GetPointY(i)>1.0e-10){ //0){
+ ////chi=(exper->Eval(anglecentres[i])-theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
+ //chi=(exper->GetPointY(i) - theory->Eval(anglecentres[i]) ) / (exper->GetErrorY(i));
+ ////cout << "COMPARE::::: " << exper->Eval(anglecentres[i]) << " TO " << exper->GetPointY(i) << endl;
+ //Chi2 +=chi*chi;
+
+ Chi2 += pow((exper->GetPointY(i) - theory->Eval(anglecentres[i])), 2 ) / theory->Eval(anglecentres[i]);
+
+ }
+ }
+ if(loud){cout << "Chi2 = " << Chi2 << endl;}
+ return Chi2;
+ //cout << setprecision(3);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+double ToMininize(const double* parameter){
+ static int f = 0 ;
+ TGraph* g = new TGraph();
+ if(doMix){
+ double* X = currentThry->GetX(); // gets valies from global g1 = tgraph passed to find norm
+ double* Y1 = currentThry->GetY();
+ double* Y2 = currentThry2->GetY();
+ for(int i = 0 ; i < currentThry->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y1[i] + parameter[1]*Y2[i]); // scales this tgraph by parameter
+ }
+ } else {
+ double* X = currentThry->GetX(); // gets valies from global g1 = tgraph passed to find norm
+ double* Y = currentThry->GetY();
+ for(int i = 0 ; i < currentThry->GetN() ; i++){
+ g->SetPoint(g->GetN(),X[i],parameter[0]*Y[i]); // scales this tgraph by parameter
+ }
+ }
+
+ double chi2 = Chi2(g,staticExp); //compares theory tgraph to experimental tgrapherrors by chi2
+ return chi2;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* FindNormalisation(TGraph* theory, TGraphErrors* experiment){
+ /* (dSdO)meas = S * (dSdO)calc */
+
+
+cout << "!! !! !! TEST HERE1" << endl;
+
+
+ // Set global variable
+ currentThry = theory;
+ staticExp = experiment;
+
+ // Construct minimiser
+ const char* minName ="Minuit";const char* algoName="Migrad";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+cout << "!! !! !! TEST HERE2" << endl;
+
+ // Number of parameters (should only be 1 for me)
+ int mysize = 1;
+
+ // Create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+cout << "!! !! !! TEST HERE3" << endl;
+
+ // Set range of parameter(??)
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ parameter[i] = 0.8;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ min->SetLimitedVariable(i,name,parameter[i],0.01,0,10);
+ }
+
+ min->SetPrintLevel(3);
+
+ ///// TO IMPROVE: FIND WAY OF OBTAINING NDF AND PRINT CHI2/NDF /////
+
+ // Minimise
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ cout << "#################### " << min->ErrorDef() << endl;
+
+ // Write out
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
+ }
+
+ /* Store S value in global variable, to access for drawing on plots */
+ globalS = xs[0];
+ globalSerr = err[0];
+
+cout << "!! !! !! TEST HERE4" << endl;
+
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = theory->GetX();
+ double* Y = theory->GetY();
+ if(loud){
+ cout << setprecision(8);
+ cout << "Start: X[0] = " << theory->GetPointX(4) << " Y[0] = " << theory->GetPointY(4) << endl;
+ cout << "multip by " << xs[0] << endl;
+ }
+
+ //for(int i=0; i<theory->GetN(); i++){ g->SetPoint(g->GetN(),X[i],xs[0]*Y[i]); }
+ for(int i=0; i<theory->GetN(); i++){ g->SetPoint(i,X[i],xs[0]*Y[i]); }
+
+cout << "!! !! !! TEST HERE5" << endl;
+
+ if(loud){
+ cout << "End: X[0] = " << g->GetPointX(4) << " Y[0] = " << g->GetPointY(4) << endl;
+ cout << setprecision(3);
+ }
+
+ return g;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+TGraph* FindNormalisation(TGraph* theory, TGraph* theory2, TGraphErrors* experiment){
+ /* (dSdO)meas = S * (dSdO)calc */
+
+ // Set global variable
+ currentThry = theory;
+ currentThry2 = theory2;
+ staticExp = experiment;
+
+ // Construct minimiser
+ const char* minName ="Minuit";const char* algoName="Migrad";
+ ROOT::Math::Minimizer* min =
+ ROOT::Math::Factory::CreateMinimizer(minName, algoName);
+ min->SetValidError(true);
+
+ // Number of parameters (should only be 1 for me)
+ int mysize = 1;
+ if(doMix){mysize=2;}
+
+ // Create funciton wrapper for minimizer
+ // a IMultiGenFunction type
+ ROOT::Math::Functor f(&ToMininize,mysize);
+ min->SetFunction(f);
+
+ // Set range of parameter(??)
+ double* parameter = new double[mysize];
+ for(unsigned int i = 0 ; i < mysize ; i++){
+ parameter[i] = 0.8;
+ char name[4];
+ sprintf(name,"S%d",i+1);
+ min->SetLimitedVariable(i,name,parameter[i],0.01,0,10);
+ }
+
+ // Minimise
+ min->Minimize();
+ const double *xs = min->X();
+ const double *err = min->Errors();
+
+ // Write out
+ for(int i = 0 ; i < mysize ; i++){
+ cout << Form("S%d=",i+1) << xs[i] << "(" << err[i] << ")" << endl;
+ }
+
+ /* Store S value in global variable, to access for drawing on plots */
+ globalS = xs[0];
+ globalSerr = err[0];
+ globalS2 = xs[1];
+ globalSerr2 = err[1];
+
+ // Return the Fitted CS
+ TGraph* g = new TGraph();
+ double* X = theory->GetX();
+ double* Y1 = theory ->GetY();
+ double* Y2 = theory2->GetY();
+ if(loud){
+ cout << setprecision(8);
+ cout << "Start: X[0] = " << theory->GetPointX(4) << " Y[0] = " << theory->GetPointY(4) << endl;
+ cout << "multip by " << xs[0] << endl;
+ cout << "Start: X[1] = " << theory2->GetPointX(4) << " Y[0] = " << theory2->GetPointY(4) << endl;
+ cout << "multip by " << xs[1] << endl;
+ }
+
+ for(int i=0; i<theory->GetN(); i++){
+ g->SetPoint(i, X[i], xs[0]*Y1[i] + xs[1]*Y2[i]);
+ }
+
+ if(loud){
+ cout << "End: X = " << g->GetPointX(4) << " Y = " << g->GetPointY(4) << endl;
+ cout << setprecision(3);
+ }
+
+ return g;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void CS(double Energy, double Spin, double spdf, double angmom, double binning, double GammaGateEnergy){
+ GammaGate=true;
+ globGammaGate = GammaGateEnergy;
+
+ CS(Energy, Spin, spdf, angmom, binning, "");
+}
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 44eac574080ab87a6e67a30d3edd047ada187414..13f86830d24b27f72e5943a0f336fbdf7bd7ab2f 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -21,6 +21,9 @@ NPL::Reaction Tidt("47Ti(d,t)46Ti@355");
NPL::Reaction Tidd("47Ti(d,d)47Ti@355");
NPL::Reaction Ti12C12C("47Ti(12C,12C)47Ti@355");
+EColor myColours[5] = {(EColor)(kRed+1), (EColor)(kGreen+1), (EColor)(kBlue), (EColor)(kMagenta+1)};
+
+
void KnownLines_Ex(bool isVertical, double rangemin, double rangemax, Style_t lType, Color_t lColour);
void AddGammaLines(TH1F* hist, double particle, double ymax);
void AddPlacedGammas(TH1F* hist, double ymax);
@@ -28,6 +31,11 @@ void AddPlacedGammas(TH1F* hist, double ymax);
double tCentre;
double tRange;
+/* DRAWING FUNCTIONS */
+string timegate; /* defined by choice of dp or dt */
+string det_gate; /* defined by choice of dp or dt */
+string exclBmDcy = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003";
+
/* BASE FUNCTIONS */
TF1* f_efficAGATA(){
TF1 *f_E = new TF1("fit_1","TMath::Exp([0]+[1]*TMath::Log(x)+[2]*pow(TMath::Log(x),2.0)+[3]*pow(TMath::Log(x),3.0)+[4]*pow(TMath::Log(x),4.0))",10,6000);
@@ -71,8 +79,18 @@ void LoadChain47Kdp(){
/*************************/
- files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_PartI.root");
- files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_PartII.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_PartII.root");
+
+ //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_RetestForFailure_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdp_18Oct22_RetestForFailure_PartII.root");
+
+ //files.push_back("../../../Outputs/Analysis/20Feb23_NewGammaMatching_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/20Feb23_NewGammaMatching_PartII.root");
+
+ files.push_back("../../../Outputs/Analysis/23Feb23_AfterNPLGRITchanges_PartI.root");
+ files.push_back("../../../Outputs/Analysis/23Feb23_AfterNPLGRITchanges_PartII.root");
+
chain = Chain("PhysicsTree",files,true);
}
@@ -83,8 +101,11 @@ void LoadChain47Kdt(){
/* Offset MM1 timing by -5 */
/* Push MM1-4 +200mm in Z */
- files.push_back("../../../Outputs/Analysis/47Kdt_18Oct22_PartI.root");
- files.push_back("../../../Outputs/Analysis/47Kdt_18Oct22_PartII.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_18Oct22_PartI.root");
+ //files.push_back("../../../Outputs/Analysis/47Kdt_18Oct22_PartII.root");
+
+ files.push_back("../../../Outputs/Analysis/47Kdt_30Mar23_AfterNPLGRITChanges_PartI.root");
+ files.push_back("../../../Outputs/Analysis/47Kdt_30Mar23_AfterNPLGRITChanges_PartII.root");
chain = Chain("PhysicsTree",files,true);
}
@@ -161,29 +182,24 @@ void DrawParticleStates(TCanvas* canvas){
Sn->Draw();
-/* if(reactionName=="47K(d,p)"){
+ if(reactionName=="47K(d,p)"){
TLine** lines = new TLine*[numPeaks];
for(int j = 0; j < numPeaks; j++){
lines[j] = new TLine(means[j],0.0,means[j],max);
lines[j]->SetLineStyle(8);
- lines[j]->SetLineColor(kGray);
+ lines[j]->SetLineColor(kViolet);
+ lines[j]->Draw();
+ }
+ }
+ else if(reactionName=="47K(d,t)"){
+ TLine** lines = new TLine*[numPeaks_dt];
+ for(int j = 0; j < numPeaks_dt; j++){
+ lines[j] = new TLine(means_dt[j],0.0,means_dt[j],max);
+ lines[j]->SetLineStyle(8);
+ lines[j]->SetLineColor(kViolet);
lines[j]->Draw();
}
}
-*/
-// if(reactionName=="47K(d,t)"){
- TLine *l1945 = new TLine(1.945, 0.0, 1.945, max);
- l1945->SetLineStyle(kDashed);
- l1945->Draw();
-
- TLine *l2233 = new TLine(2.233, 0.0, 2.233, max);
- l2233->SetLineStyle(kDashed);
- l2233->Draw();
-
- TLine *l3344 = new TLine(3.344, 0.0, 3.344, max);
- l3344->SetLineStyle(kDashed);
- l3344->Draw();
-// }
}
void plot_kine(NPL::Reaction r, double Ex,Color_t c,int w, int s){
@@ -212,13 +228,6 @@ void AddTiStates(double E){
g->Draw("c");
}
-/* DRAWING FUNCTIONS */
-
-string timegate; /* defined by choice of dp or dt */
-string det_gate; /* defined by choice of dp or dt */
-string reactionName; /* defined by choice of dp or dt */
-string exclBmDcy = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003";
-
void Draw_1DGamma(){
string gate = timegate;
TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
@@ -246,7 +255,7 @@ void Draw_1DGamma_DetGate(){
TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
gStyle->SetOptStat(0);
- chain->Draw("AddBack_EDC>>Eg(5000,0,5)",gate.c_str(),"");
+ chain->Draw("AddBack_EDC>>Eg(9000,0,9)",gate.c_str(),"");
TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
Eg->GetYaxis()->SetTitle("Counts / 0.001 MeV");
@@ -259,15 +268,15 @@ void Draw_1DGamma_DetGate_x10x100(){
TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
gStyle->SetOptStat(0);
- chain->Draw("AddBack_EDC>>Eg(1250,0,5)",gate.c_str(),"");
+ chain->Draw("AddBack_EDC>>Eg(1000,0,5)",gate.c_str(),"");
TH1F* Eg = (TH1F*) gDirectory->Get("Eg");
Eg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
- Eg->GetYaxis()->SetTitle("Counts / 0.004 MeV");
+ Eg->GetYaxis()->SetTitle("Counts / 0.005 MeV");
- for(int b = 126; b<501; b++){
+ for(int b = Eg->FindBin(0.506); b< Eg->FindBin(2.0); b++){
Eg->SetBinContent(b,(Eg->GetBinContent(b) * 10.));
}
- for(int b = 501; b<1250; b++){
+ for(int b = Eg->FindBin(2.006); b< Eg->FindBin(5.0); b++){
Eg->SetBinContent(b,(Eg->GetBinContent(b) * 100.));
}
@@ -319,7 +328,7 @@ void Draw_1DGamma_DetGate_SplitCanv(){
void Load_1DGamma_MG(){
TCanvas *cEg = new TCanvas("cEg","cEg",1000,1000);
TH1F *hEgMG = new TH1F("hEg","Loaded 1D Gamma Spectrum, MG gated",600,-15,15);
- TFile *file = new TFile("LoadHistograms/Load_1DGamma_MG.root","READ");
+ TFile *file = new TFile("LoadHistograms_48K/Load_1DGamma_DetGate.root","READ");
hEgMG = (TH1F*)file->Get("Eg");
hEgMG->Draw();
}
@@ -342,8 +351,6 @@ void Load_1DParticle_SubPhaseSpace(){
}
void Draw_1DParticle(){
- string gate2;
- string gate3;
string gate = timegate
+ " && " + det_gate
+ " && Ex@.size()==1";
@@ -362,13 +369,13 @@ void Draw_1DParticle(){
Ep->GetYaxis()->SetTitle("Counts / 0.05 MeV");
}
- if(reactionName=="47K(d,p)"){DrawParticleStates(cEx);}
+// if(reactionName=="47K(d,p)"){DrawParticleStates(cEx);}
}
void Load_2DParticleGamma(){
TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000);
TH2F *hExEg = new TH2F("hExEg","Loaded 2D Particle-Gamma",600,-15,15,2500,0,5);
- TFile *file = new TFile("LoadHistograms/Load_2DParticleGamma.root","READ");
+ TFile *file = new TFile("LoadHistograms_15Feb23/Load_2DParticleGamma.root","READ");
hExEg = (TH2F*)file->Get("ExEg");
hExEg->Draw("colz");
}
@@ -384,7 +391,7 @@ void Draw_2DParticleGamma(){
TCanvas *cExEg = new TCanvas("cExEg","cExEg",1000,1000);
gStyle->SetOptStat(0);
//chain->Draw("AddBack_EDC:Ex>>ExEg(600,-15,15,2500,0,5)", gate.c_str(), "colz");
- chain->Draw("Ex:AddBack_EDC>>ExEg(5000,0,5,3000,-15,15)", gate.c_str(), "");
+ chain->Draw("Ex:AddBack_EDC>>ExEg(10000,0,5,6000,-15,15)", gate.c_str(), "");
TH1F* ExEg = (TH1F*) gDirectory->Get("ExEg");
ExEg->SetTitle("");
ExEg->GetYaxis()->SetTitle("Ex [MeV]");
@@ -422,22 +429,33 @@ void Load_2DGammaGamma(){
hEgEg->Draw("colz");
}
-void Draw_2DGammaGamma_ExcludeBeam(){
+void Draw_2DGammaGamma_New_ExcludeBeam(double binWidth){
string gate = timegate
+ " && " + exclBmDcy;
+ int bins = (int) (5./binWidth);
+
+ string draw1 = "AddBack_EDC_Event1:AddBack_EDC_Event2>>EgEg1("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+ string draw2 = "AddBack_EDC_Event2:AddBack_EDC_Event1>>EgEg2("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+
+ string titles = "Eg [Counts / " + to_string(binWidth) + " MeV]";
+
TCanvas *cEgEg = new TCanvas("cEgEg","cEgEg",1000,1000);
- chain->Draw("AddBack_EDC:AddBack_EDC2>>EgEg(2500,0,5,2500,0,5)",gate.c_str(),"colz");
- TH1F* EgEg = (TH1F*) gDirectory->Get("EgEg");
- chain->Draw("AddBack_EDC2:AddBack_EDC>>EgEg2(2500,0,5,2500,0,5)",gate.c_str(),"colz");
+ chain->Draw(draw1.c_str(),gate.c_str(),"colz");
+ TH1F* EgEg1 = (TH1F*) gDirectory->Get("EgEg1");
+ chain->Draw(draw2.c_str(),gate.c_str(),"colz");
TH1F* EgEg2 = (TH1F*) gDirectory->Get("EgEg2");
- EgEg->Add(EgEg2,1);
- EgEg->SetTitle("Egamma-Egamma");
- EgEg->GetXaxis()->SetTitle("Eg [Counts / 0.002 MeV]");
- EgEg->GetYaxis()->SetTitle("Eg [Counts / 0.002 MeV]");
- EgEg->GetXaxis()->SetRangeUser(0.005,5.0);
- EgEg->GetYaxis()->SetRangeUser(0.005,5.0);
- EgEg->Draw("colz");
+ EgEg1->Add(EgEg2,1);
+ EgEg1->SetTitle("Egamma-Egamma");
+ EgEg1->GetXaxis()->SetTitle(titles.c_str());
+ EgEg1->GetYaxis()->SetTitle(titles.c_str());
+// EgEg->GetXaxis()->SetRangeUser(0.005,5.0);
+// EgEg->GetYaxis()->SetRangeUser(0.005,5.0);
+ EgEg1->Draw("colz");
TLine *XeqY = new TLine(0,0,5,5);
XeqY->SetLineColor(kRed);
XeqY->SetLineStyle(kDashed);
@@ -658,35 +676,132 @@ void GateParticle_SeeGamma_WithBG(double particle, double width, double bg, doub
EgBG->SetFillStyle(3345);
}
-/*
-void GateGamma_SeeGamma_ExcludeBeamDecay(double gamma, double width){
- string gating
- = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003 && abs(AddBack_EDC2-"
- + to_string(gamma)
- + ")<"
- + to_string(width);
- string gating2
- = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003 && abs(AddBack_EDC-"
- + to_string(gamma)
- + ")<"
- + to_string(width);
+void GateGamma_SeeGamma_New(double gamma, double width, double binWidth){
+// string gating
+// = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003 && abs(AddBack_EDC2-"
+// + to_string(gamma)
+// + ")<"
+// + to_string(width);
+// string gating2
+// = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.003 && abs(AddBack_EDC-"
+// + to_string(gamma)
+// + ")<"
+// + to_string(width);
+//
+// string title = to_string(gamma-width) + " < Eg < " + to_string(gamma+width);
+// TCanvas *cEx_Gate = new TCanvas("cggGate","cggGate",1000,1000);
+//
+// chain->Draw("AddBack_EDC>>ggGate(999,0.005,5)",gating.c_str(),"");
+// TH1F* ggGate = (TH1F*) gDirectory->Get("ggGate");
+// ggGate->GetXaxis()->SetTitle("Eg [MeV]");
+// ggGate->GetYaxis()->SetTitle("Counts / 0.005 MeV");
+// ggGate->SetTitle(title.c_str());
+//
+// chain->Draw("AddBack_EDC2>>ggGate2(999,0.005,5)",gating2.c_str(),"");
+// TH1F* ggGate2 = (TH1F*) gDirectory->Get("ggGate2");
+// ggGate->Add(ggGate2,1);
+// ggGate->Draw();
- string title = to_string(gamma-width) + " < Eg < " + to_string(gamma+width);
- TCanvas *cEx_Gate = new TCanvas("cggGate","cggGate",1000,1000);
- chain->Draw("AddBack_EDC>>ggGate(999,0.005,5)",gating.c_str(),"");
+
+ string gate = timegate
+ + " && " + exclBmDcy
+ + " && abs(AddBack_EDC_Event2 - " + to_string(gamma)
+ + " )<" + to_string(width);
+
+ int bins = (int) (5./binWidth);
+
+ string draw1 = "AddBack_EDC_Event1>>ggGate("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+
+ string titles = "Eg [Counts / " + to_string(binWidth) + " MeV]";
+
+ TCanvas *cggGate = new TCanvas("cggGate","cggGate",1000,1000);
+ chain->Draw(draw1.c_str(),gate.c_str(),"");
TH1F* ggGate = (TH1F*) gDirectory->Get("ggGate");
- ggGate->GetXaxis()->SetTitle("Eg [MeV]");
- ggGate->GetYaxis()->SetTitle("Counts / 0.005 MeV");
- ggGate->SetTitle(title.c_str());
+ ggGate->GetXaxis()->SetTitle(titles.c_str());
+ ggGate->Draw("");
+
+
+
+
+
- chain->Draw("AddBack_EDC2>>ggGate2(999,0.005,5)",gating2.c_str(),"");
- TH1F* ggGate2 = (TH1F*) gDirectory->Get("ggGate2");
- ggGate->Add(ggGate2,1);
- ggGate->Draw();
}
-*/
+void GateGamma_SeeGamma_New_WithBG(double gamma, double width, double bg, double bgwidth, double binWidth){
+
+ string gate1 = timegate
+ + " && " + exclBmDcy
+ + " && abs(AddBack_EDC_Event2 - " + to_string(gamma)
+ + " )<" + to_string(width);
+ string gate2 = timegate
+ + " && " + exclBmDcy
+ + " && abs(AddBack_EDC_Event1 - " + to_string(gamma)
+ + " )<" + to_string(width);
+
+ string bggate1 = timegate
+ + " && " + exclBmDcy
+ + " && abs(AddBack_EDC_Event2 - " + to_string(bg)
+ + " )<" + to_string(bgwidth);
+
+ string bggate2 = timegate
+ + " && " + exclBmDcy
+ + " && abs(AddBack_EDC_Event1 - " + to_string(bg)
+ + " )<" + to_string(bgwidth);
+
+
+ int bins = (int) (5./binWidth);
+ double binScale = width/bgwidth;
+
+ string draw1 = "AddBack_EDC_Event1>>gg1Gate("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+ string draw2 = "AddBack_EDC_Event2>>gg2Gate("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+
+ string bgdraw1 = "AddBack_EDC_Event1>>bg1Gate("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+ string bgdraw2 = "AddBack_EDC_Event2>>bg2Gate("
+ + to_string(bins) + ",0,5,"
+ + to_string(bins) + ",0,5)";
+
+ string titles = "Eg [Counts / " + to_string(binWidth) + " MeV]";
+
+ TCanvas *cggGate = new TCanvas("cggGate","cggGate",1000,1000);
+ chain->Draw(draw1.c_str(),gate1.c_str(),"");
+ TH1F* gg1Gate = (TH1F*) gDirectory->Get("gg1Gate");
+ chain->Draw(draw2.c_str(),gate2.c_str(),"");
+ TH1F* gg2Gate = (TH1F*) gDirectory->Get("gg2Gate");
+ chain->Draw(bgdraw1.c_str(),bggate1.c_str(),"");
+ TH1F* bg1Gate = (TH1F*) gDirectory->Get("bg1Gate");
+ chain->Draw(bgdraw2.c_str(),bggate2.c_str(),"");
+ TH1F* bg2Gate = (TH1F*) gDirectory->Get("bg2Gate");
+
+ gg1Gate->Add(gg2Gate,1.);
+ bg1Gate->Add(bg2Gate,1.);
+
+ gg1Gate->GetXaxis()->SetTitle(titles.c_str());
+
+ gg1Gate->SetLineColor(kGreen);
+ bg1Gate->SetLineColor(kRed);
+
+ if(width!=bgwidth){
+ bg1Gate->Scale(binScale);
+ }
+
+ gg1Gate->Draw("");
+ bg1Gate->Draw("same");
+
+
+
+
+}
+
+/*
void GateGamma_SeeGamma(double gamma, double width){
string gating = "abs(AGATA_GammaE-2.013)>0.004 && abs(AGATA_GammaE-0.511)>0.003 && abs(AGATA_GammaE-0.564)>0.004 && abs(AGATA_GammaE-0.586)>0.00 && abs(AddBack_EDC2-"
+ to_string(gamma)
@@ -782,7 +897,7 @@ void GateGamma_SeeGamma_WithBG(double gamma, double width, double bg, double wid
ggGate->Draw();
ggBG->Draw("same");
}
-
+*/
void CompareExsAt4MeV(){
TCanvas *cExCompare = new TCanvas("cExCompare","cExCompare",1000,1000);
chain->Draw("AddBack_EDC>>gate3p0(1000,0,10)",
@@ -1231,7 +1346,6 @@ void ExThetaCM(){
}
}
-
void ExThetaLab(double gamma, double width){
TCanvas *diagnoseTheta = new TCanvas("diagnoseTheta","diagnoseTheta",1000,1000);
string gate = timegate
@@ -1292,6 +1406,8 @@ void ExThetaLab(double gamma, double width){
void ELabThetaLab(){
TCanvas *cELabTLab = new TCanvas("cELabTLab","cELabTLab",1000,1000);
gStyle->SetOptStat(0);
+ gStyle->SetPalette(kViridis);
+ gPad->SetFixedAspectRatio();
string gate = timegate
+ " && " + det_gate;
@@ -1300,7 +1416,7 @@ void ELabThetaLab(){
}
- chain->Draw("ELab:ThetaLab>>hKine(360,0,180,500,0,10)",gate.c_str(),"col");
+ chain->Draw("ELab:ThetaLab>>hKine(360,0,180,1000,0,10)",gate.c_str(),"col");
TH2F* hKine = (TH2F*) gDirectory->Get("hKine");
hKine->SetTitle("");
hKine->GetXaxis()->SetTitle("#theta_{lab} [deg]");
@@ -1335,8 +1451,30 @@ void ELabThetaLab(){
plot_kine(Kdp, 4.644, kBlack, 2, 1);
plot_kine(Cadp, 0.000, kRed, 2, 1);
plot_kine(Tidp, 0.000, kBlue, 2, 1);
- plot_kine(Scdp, 0.000, kGreen, 2, 1);
- plot_kine(K46dp, 0.000, kViolet, 2, 1);
+ plot_kine(Scdp, 0.000, kViolet, 2, 1);
+ //plot_kine(Scdp, 0.000, kGreen, 2, 1);
+ //plot_kine(K46dp, 0.000, kViolet, 2, 1);
+
+
+ TLatex *Ca = new TLatex(.5,.5,"^{47}Ca(d,p)");
+ Ca->SetTextColor(kRed);
+ Ca->SetTextSize(0.05);
+ Ca->SetX(130.);
+ Ca->SetY(5.);
+ Ca->Draw("same");
+ TLatex *Ti = new TLatex(.5,.5,"^{47}Ti(d,p)");
+ Ti->SetTextColor(kBlue);
+ Ti->SetTextSize(0.05);
+ Ti->SetX(140.);
+ Ti->SetY(4.);
+ Ti->Draw("same");
+ TLatex *Sc = new TLatex(.5,.5,"^{47}Sc(d,p)");
+ Sc->SetTextColor(kViolet);
+ Sc->SetTextSize(0.05);
+ Sc->SetX(150.);
+ Sc->SetY(3.);
+ Sc->Draw("same");
+
}
}
@@ -2091,8 +2229,6 @@ void GateThetaCM_MultiWrite(double startTheta, double finishTheta, int numGates,
file->ls();
}
-
-
void GateThetaLab_MultiWrite(double startTheta, double finishTheta, int numGates, double binsize, int MGX){
string core = timegate
+ " && " + det_gate
@@ -2169,13 +2305,14 @@ void CompareThetaLabGatesFromFile(){
//TFile* oldF = new TFile("GateThetaLabHistograms_11Apr22_20angles.root","READ");
//TFile* oldF = new TFile("GateThetaLabHistograms_11Jul22.root","READ");
//TFile* oldF = new TFile("GateThetaLabHistograms_29Aug22_TrueStripRemoval_0p05.root","READ");
- TFile* oldF = new TFile("GateThetaLabHistograms_14Oct22_bin0p05.root","READ");
+ TFile* oldF = new TFile("GateThetaLabHistograms_Store/GateThetaLabHistograms_14Oct22_bin0p05.root","READ");
TList* oldL = (TList*) oldF->FindObjectAny("GateThetaLabHistograms");
//TFile* newF = new TFile("GateThetaLabHistograms_11Jul22.root","READ");
//TFile* newF = new TFile("GateThetaLabHistograms_Test2um.root","READ");
//TFile* newF = new TFile("GateThetaLabHistograms_30Aug22_2p06um.root","READ");
- TFile* newF = new TFile("GateThetaLabHistograms_14Oct22_2_bin0p05.root","READ");
+ TFile* newF = new TFile("GateThetaLabHistograms_Store/GateThetaLabHistograms_14Oct22_2_bin0p05.root","READ");
+ //TFile* newF = new TFile("GateThetaLabHistograms_47Kdp_18Oct22_bin0p2.root","READ");
TList* newL = (TList*) newF->FindObjectAny("GateThetaLabHistograms");
double minTheta=105.;
@@ -2338,10 +2475,12 @@ void GatePhaseSpaceByThetaLab_MultiWrite(double startTheta, double finishTheta,
TList* list = new TList();
//TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_02Mar_47Kdp_PhaseSpace.root","READ");
- TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_PhaseSpace_11Jul22.root","READ");
+ //TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_PhaseSpace_11Jul22.root","READ");
+ TFile* psfile = new TFile("../../../Outputs/Analysis/Sim_18Oct22_PhaseSpace.root","READ");
TTree* PSTree = (TTree*) psfile->FindObjectAny("PhysicsTree");
for (int i=0; i<numGates; i++){
+ cout << GREEN << "Writing gate " << i+1 << "/" << numGates << RESET << endl;
double minTheta = startTheta + (i * gatesize);
string title = to_string((int) minTheta)+" < ThetaLab < "+to_string((int) (minTheta+gatesize));
string gating = core
@@ -2728,6 +2867,8 @@ void Figure_GateGamma_SeeParticle(double gamma, double width, double bg, double
+ to_string(bg)
+ ")<"
+ to_string(widthbg);
+
+ string gateBase = timegate + "&&" + det_gate + " && Ex@.size()==1";
string gammagate = timegate
+ " && " + det_gate
@@ -2739,7 +2880,8 @@ void Figure_GateGamma_SeeParticle(double gamma, double width, double bg, double
TCanvas *cFig_GGSP = new TCanvas("cFig_GGSP","cFig_GGSP",1000,1000);
gStyle->SetOptStat(0);
- cFig_GGSP->Divide(1,2,0.005,0.005,0);
+ //cFig_GGSP->Divide(1,2,0.005,0.005,0);
+ cFig_GGSP->Divide(2,1,0.005,0.005,0);
cFig_GGSP->cd(1);
gStyle->SetPadLeftMargin(0.10);
gStyle->SetPadRightMargin(0.01);
@@ -2752,8 +2894,14 @@ void Figure_GateGamma_SeeParticle(double gamma, double width, double bg, double
gPad->SetTicky();
cFig_GGSP->cd(2);
+ string drawBase = "Ex>>hBase(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(drawBase.c_str(),gateBase.c_str(),"");
+ TH1F* ExBase = (TH1F*) gDirectory->Get("hBase");
+ ExBase->Scale(0.05);
+
string draw1 = "Ex>>hEx(" + to_string(30./particleBinWidth) + ",-15,15)";
- chain->Draw(draw1.c_str(),gating.c_str(),"");
+ //chain->Draw(draw1.c_str(),gating.c_str(),"");
+ chain->Draw(draw1.c_str(),gating.c_str(),"same hist");
TH1F* ExGate = (TH1F*) gDirectory->Get("hEx");
ExGate->GetXaxis()->SetTitle("Ex [MeV]");
string nameExYaxis = "Counts / " + to_string(particleBinWidth) + " MeV";
@@ -2826,6 +2974,175 @@ void Figure_GateGamma_SeeParticle(double gamma, double width, double bg, double
}
+
+void Figure_GateGamma_SeeParticle(
+ double gamma1, double width1,
+ double gamma2, double width2,
+ double gamma3, double width3,
+ double gammaBinWidth, double particleBinWidth){
+
+ string gate1 = timegate + "&&" + det_gate + " && Ex@.size()==1 && abs(AddBack_EDC-"
+ + to_string(gamma1)
+ + ")<"
+ + to_string(width1);
+ string gate2 = timegate + "&&" + det_gate + " && Ex@.size()==1 && abs(AddBack_EDC-"
+ + to_string(gamma2)
+ + ")<"
+ + to_string(width2);
+ string gate3 = timegate + "&&" + det_gate + " && Ex@.size()==1 && abs(AddBack_EDC-"
+ + to_string(gamma3)
+ + ")<"
+ + to_string(width3);
+
+ string gateBase = timegate + "&&" + det_gate + " && Ex@.size()==1";
+
+ string gammagate = timegate
+ + " && " + det_gate
+ + " && " + exclBmDcy;
+
+ //double ratio = width/widthbg;
+
+ ////////////////////////////////////////////////
+
+ TF1* func = f_efficAGATA();
+
+
+ cout << " !!! " << 1./(func->Eval(gamma1*1000.)/100.) << endl;
+ cout << " !!! " << 1./(func->Eval(gamma2*1000.)/100.) << endl;
+ cout << " !!! " << 1./(func->Eval(gamma3*1000.)/100.) << endl;
+
+ TCanvas *cFig_GGSP = new TCanvas("cFig_GGSP","cFig_GGSP",1000,1000);
+ gStyle->SetOptStat(0);
+ //cFig_GGSP->Divide(1,2,0.005,0.005,0);
+ cFig_GGSP->Divide(2,1,0.005,0.005,0);
+ cFig_GGSP->cd(1);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.01);
+ gPad->SetTickx();
+ gPad->SetTicky();
+ cFig_GGSP->cd(2);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.01);
+ gPad->SetTickx();
+ gPad->SetTicky();
+
+ cFig_GGSP->cd(2);
+ string drawBase = "Ex>>hExBase(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(drawBase.c_str(),gateBase.c_str(),"");
+ TH1F* ExBase = (TH1F*) gDirectory->Get("hExBase");
+ ExBase->GetXaxis()->SetTitle("Ex [MeV]");
+ string nameExYaxis = "Counts / " + to_string(particleBinWidth) + " MeV";
+ ExBase->GetYaxis()->SetTitle(nameExYaxis.c_str());
+ ExBase->SetLineWidth(0);
+ ExBase->SetFillColor(kBlack);
+ ExBase->SetFillStyle(3003);
+ ExBase->SetTitle("");
+ ExBase->GetXaxis()->SetTitleSize(0.05);
+ ExBase->GetXaxis()->SetLabelSize(0.05);
+ ExBase->GetYaxis()->SetTitleSize(0.05);
+ ExBase->GetYaxis()->SetLabelSize(0.05);
+ ExBase->GetYaxis()->SetTitleOffset(0.5);
+ ExBase->GetYaxis()->SetNdivisions(520);
+ ExBase->GetXaxis()->SetRangeUser(-0.5,7.);
+ ExBase->GetXaxis()->SetRangeUser(0.0,150.);
+ //ExBase->Scale(0.07);
+ //ExBase->Scale(0.08);
+ //ExBase->Scale(0.03);
+
+ string draw1 = "Ex>>hEx1(" + to_string(30./particleBinWidth) + ",-15,15)";
+ //chain->Draw(draw1.c_str(),gating.c_str(),"");
+ chain->Draw(draw1.c_str(),gate1.c_str(),"same hist");
+ TH1F* Ex1 = (TH1F*) gDirectory->Get("hEx1");
+ Ex1->GetXaxis()->SetTitle("Ex [MeV]");
+ Ex1->GetYaxis()->SetTitle(nameExYaxis.c_str());
+ Ex1->SetLineColor(kGreen);
+ Ex1->SetFillColor(kGreen);
+ Ex1->SetFillStyle(3002);
+ Ex1->SetTitle("");
+ Ex1->GetXaxis()->SetTitleSize(0.05);
+ Ex1->GetXaxis()->SetLabelSize(0.05);
+ Ex1->GetYaxis()->SetTitleSize(0.05);
+ Ex1->GetYaxis()->SetLabelSize(0.05);
+ Ex1->GetYaxis()->SetTitleOffset(0.5);
+ Ex1->GetYaxis()->SetNdivisions(520);
+ Ex1->GetXaxis()->SetRangeUser(-1.,7.);
+ Ex1->Scale(1./(func->Eval(gamma1*1000.)/100.));
+// Ex1->Scale((1.0/funci));
+
+ string draw2 = "Ex>>hEx2(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(draw2.c_str(),gate2.c_str(),"same hist");
+ TH1F* Ex2 = (TH1F*) gDirectory->Get("hEx2");
+ Ex2->SetLineColor(kRed);
+ Ex2->SetFillColor(kRed);
+ Ex2->SetFillStyle(3002);
+ Ex2->SetTitle("");
+ Ex2->GetXaxis()->SetRangeUser(-1.,7.);
+ Ex2->Scale(1./(func->Eval(gamma2*1000.)/100.));
+ Ex2->Draw("same hist");
+
+ string draw3 = "Ex>>hEx3(" + to_string(30./particleBinWidth) + ",-15,15)";
+ chain->Draw(draw3.c_str(),gate3.c_str(),"same hist");
+ TH1F* Ex3 = (TH1F*) gDirectory->Get("hEx3");
+ Ex3->SetLineColor(kBlue);
+ Ex3->SetFillColor(kBlue);
+ Ex3->SetFillStyle(3002);
+ Ex3->SetTitle("");
+ Ex3->GetXaxis()->SetRangeUser(-1.,7.);
+ Ex3->Scale(1./(func->Eval(gamma3*1000.)/100.));
+ Ex3->Draw("same hist");
+
+ cFig_GGSP->Update();
+ double maxEx = cFig_GGSP->cd(2)->GetUymax();
+ TLine *Sn = new TLine(4.644, 0.0, 4.644, maxEx);
+ Sn->SetLineColor(kBlack); Sn->SetLineStyle(7);
+ Sn->Draw();
+ TLine *gs = new TLine(0.000, 0.0, 0.000, maxEx);
+ gs->SetLineColor(kBlack); gs->SetLineStyle(1);
+ gs->Draw();
+
+ cFig_GGSP->cd(1);
+ string drawg = "AddBack_EDC>>hEg(" + to_string(5./gammaBinWidth) + ",0,5)";
+ chain->Draw(drawg.c_str(),gammagate.c_str(),"");
+ TH1F* hEg = (TH1F*) gDirectory->Get("hEg");
+ hEg->SetLineColor(kBlack);
+ string nameEgYaxis = "Counts / " + to_string(gammaBinWidth) + " MeV";
+ hEg->GetXaxis()->SetTitle("E_{#gamma} [MeV]");
+ hEg->GetYaxis()->SetTitle(nameEgYaxis.c_str());
+ double zoomMin = min((double) floor((gamma1-width1)*10)/10., (double) floor((gamma2-width2)*10)/10.);
+ double zoomMax = max((double) ceil((gamma1+width1)*10)/10., (double) ceil((gamma2+width2)*10)/10.);
+ hEg->GetXaxis()->SetRangeUser(zoomMin,zoomMax);
+ hEg->SetTitle("");
+ hEg->GetXaxis()->SetTitleSize(0.05);
+ hEg->GetXaxis()->SetLabelSize(0.05);
+ hEg->GetYaxis()->SetTitleSize(0.05);
+ hEg->GetYaxis()->SetLabelSize(0.05);
+ hEg->GetYaxis()->SetTitleOffset(0.5);
+ hEg->Draw();
+
+ cFig_GGSP->Update();
+ double max = cFig_GGSP->cd(1)->GetUymax();
+ TLine* gate1L = new TLine(gamma1-width1,0,gamma1-width1,max);
+ gate1L->SetLineColor(kGreen); gate1L->SetLineStyle(2);
+ TLine* gate1H = new TLine(gamma1+width1,0,gamma1+width1,max);
+ gate1H->SetLineColor(kGreen); gate1H->SetLineStyle(2);
+ TLine* gate2L = new TLine(gamma2-width2,0,gamma2-width2,max);
+ gate2L->SetLineColor(kRed); gate2L->SetLineStyle(2);
+ TLine* gate2H = new TLine(gamma2+width2,0,gamma2+width2,max);
+ gate2H->SetLineColor(kRed); gate2H->SetLineStyle(2);
+ TLine* gate3L = new TLine(gamma3-width3,0,gamma3-width3,max);
+ gate3L->SetLineColor(kBlue); gate3L->SetLineStyle(2);
+ TLine* gate3H = new TLine(gamma3+width3,0,gamma3+width3,max);
+ gate3H->SetLineColor(kBlue); gate3H->SetLineStyle(2);
+
+ gate1L->Draw("SAME");
+ gate1H->Draw("SAME");
+ gate2L->Draw("SAME");
+ gate2H->Draw("SAME");
+ gate3L->Draw("SAME");
+ gate3H->Draw("SAME");
+}
+
+
void Figure_TopGamma_BottomParticle(double gammaBinWidth, double particleBinWidth){
string gating = timegate + "&&" + det_gate + " && Ex@.size()==1";
string gammagate = timegate
@@ -2897,7 +3214,6 @@ void Figure_TopGamma_BottomParticle(double gammaBinWidth, double particleBinWidt
cFig_GGSP->Update();
}
-
void Figure_SolidAngle_dt(){
//string histname = "SolidAngle_CM_MM";
string histname = "SolidAngle_Lab_MM";
@@ -3014,3 +3330,642 @@ void Figure_SolidAngle_dt(){
n6->Draw("same");
}
+
+TH2F* Func_LoadIn2DPartGamma(){
+ string loadname = "./";
+ if(reactionName=="47K(d,p)"){
+ loadname.append("LoadHistograms_48K/LoadMe_2DParticleGamma.root");
+ } else if (reactionName=="47K(d,t)"){
+ loadname.append("LoadHistograms_46K/LoadMe_2DParticleGamma.root");
+ }
+
+ TFile *file = new TFile(loadname.c_str(),"READ");
+// TH2F* h2D = new TH2F("h2D","Loaded 2D Ex-EGamma",600,-15,15);
+ TH2F* h2D = (TH2F*)file->Get("ExEg");
+
+ return h2D;
+}
+
+TH2F* Func_Rebin2DPartGamma(TH2F* h2D, double gammaBin, double particleBin){
+ double gammaBinMin = h2D->GetXaxis()->GetBinWidth(2);
+ double particleBinMin = h2D->GetYaxis()->GetBinWidth(2);
+
+ if(gammaBin<gammaBinMin || particleBin<particleBinMin){
+ cout << RED << endl;
+ cout << "BINNING TOO SMALL!" << endl;
+ cout << "Minimum binning values are..." << endl;
+ cout << " Eg: " << gammaBinMin << endl;
+ cout << " Ex: " << particleBinMin << endl;
+ cout << "Try again..." << endl;
+ cout << RESET << endl;
+ return 0;
+ }
+
+ double rebinEg = (gammaBin / gammaBinMin);
+ double rebinEx = (particleBin / particleBinMin);
+
+ cout << rebinEg << " " << rebinEx << endl;
+
+ h2D->RebinY(rebinEx);
+ h2D->RebinX(rebinEg);
+
+ return h2D;
+ //COULD MAKE FASTER USING POINTER TO POINTER, BUT CBA
+}
+
+
+TH2F* Func_Rebin2DGammaGamma(TH2F* h2D, double gammaBin1, double gammaBin2){
+ double gammaBinMin = h2D->GetXaxis()->GetBinWidth(2);
+
+ if(gammaBin1<gammaBinMin || gammaBin2<gammaBinMin){
+ cout << RED << endl;
+ cout << "BINNING TOO SMALL!" << endl;
+ cout << "Minimum binning values are..." << endl;
+ cout << " Eg: " << gammaBinMin << endl;
+ cout << "Try again..." << endl;
+ cout << RESET << endl;
+ return 0;
+ }
+
+ double rebinEg1 = (gammaBin1 / gammaBinMin);
+ double rebinEg2 = (gammaBin2 / gammaBinMin);
+
+ cout << rebinEg1 << " " << rebinEg2 << endl;
+
+ h2D->RebinY(rebinEg2);
+ h2D->RebinX(rebinEg1);
+
+ return h2D;
+ //COULD MAKE FASTER USING POINTER TO POINTER, BUT CBA
+}
+
+
+TH1F* Func_Gater(TH2F* h2D, const char* axis, double gateMean, double gateWidth, Color_t gateColour){
+
+ TH1F* output;
+
+ if(strncmp(axis,"X",1)==0){
+ output = (TH1F*)h2D->ProjectionX("output", h2D->GetYaxis()->FindBin(gateMean-gateWidth), h2D->GetYaxis()->FindBin(gateMean+gateWidth));
+ output->GetXaxis()->SetRangeUser(0.0, 5.0);
+ } else if(strncmp(axis,"Y",1)==0){
+ output = (TH1F*)h2D->ProjectionY("output", h2D->GetXaxis()->FindBin(gateMean-gateWidth), h2D->GetXaxis()->FindBin(gateMean+gateWidth));
+ output->GetXaxis()->SetRangeUser(-1.0, 7.0);
+ } else {
+ cout << " ERROR: AXIS OPTION INCORRECT" << endl;
+ return 0;
+ }
+ output->SetLineColor(gateColour);
+
+ return output;
+}
+
+
+TH2F* Func_LoadIn2DGammaGamma(){
+ string loadname = "./";
+ if(reactionName=="47K(d,p)"){
+ loadname.append("LoadHistograms_48K/LoadMe_2DGammaGamma.root");
+ }
+
+ TFile *file = new TFile(loadname.c_str(),"READ");
+ TH2F* h2D = (TH2F*)file->Get("EgEg");
+
+ return h2D;
+}
+
+void Func_AddGatingLines(double centre, double width, double height, EColor col){
+ TLine *lLow = new TLine(centre-width, 0.0, centre-width, height);
+ lLow->SetLineColor(col);
+ lLow->SetLineStyle(7);
+ lLow->Draw("SAME");
+ TLine *lHig = new TLine(centre+width, 0.0, centre+width, height);
+ lHig->SetLineColor(col);
+ lHig->SetLineStyle(7);
+ lHig->Draw("SAME");
+}
+
+void GateGamma_SeeParticle_LoadFromFile(double gamma, double width, double gammaBin, double particleBin){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,1000);
+ cGGSP->Divide(1,2);
+ double zoom = 0.2;
+
+ cGGSP->cd(1);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",0,10000);
+ hEg->GetXaxis()->SetRangeUser(gamma-zoom, gamma+zoom);
+ hEg->SetLineColor(kBlack);
+ hEg->Draw();
+ double height = hEg->GetMaximum();
+ Func_AddGatingLines(gamma, width, height, kRed);
+
+
+
+ cGGSP->cd(2);
+
+ TH1F* hEx1 = Func_Gater(h2D, "Y", gamma, width, kBlack);
+
+ //TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",
+ // h2D->GetXaxis()->FindBin(gamma-width),
+ // h2D->GetXaxis()->FindBin(gamma+width)
+ // );
+ //hEx->GetXaxis()->SetRangeUser(-1.0, 7.0);
+ //hEx->SetLineColor(kBlack);
+ hEx1->Draw();
+
+
+}
+
+void GateGamma_SeeParticle_LoadFromFile(double gamma1, double width1, double gamma2, double width2, double gammaBin, double particleBin){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,1000);
+ cGGSP->Divide(1,2);
+ double zoom = 0.2;
+
+ cGGSP->cd(1);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",0,10000);
+ hEg->GetXaxis()->SetRangeUser(gamma1-zoom, gamma1+zoom);
+ hEg->SetLineColor(kBlack);
+ hEg->Draw();
+ double height = hEg->GetMaximum();
+ Func_AddGatingLines(gamma1, width1, height, kRed);
+ Func_AddGatingLines(gamma2, width2, height, kBlue);
+
+ cGGSP->cd(2);
+
+ //// WHY DOES THIS NOT WORK??? //
+ //TH1F* hEx1 = Func_Gater(h2D, "Y", gamma1, width1, kRed);
+ //TH1F* hEx2 = Func_Gater(h2D, "Y", gamma2, width2, kBlue);
+ //hEx1->Draw();
+ //hEx2->Draw("same");
+ TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",
+ h2D->GetXaxis()->FindBin(gamma1-width1),
+ h2D->GetXaxis()->FindBin(gamma1+width1)
+ );
+ hEx->GetXaxis()->SetRangeUser(-1,7);
+ hEx->SetLineColor(kRed);
+ hEx->Draw();
+
+ TH1F* hEx2 = (TH1F*)h2D->ProjectionY("hEx2",
+ h2D->GetXaxis()->FindBin(gamma2-width2),
+ h2D->GetXaxis()->FindBin(gamma2+width2)
+ );
+ hEx2->GetXaxis()->SetRangeUser(-1,7);
+ hEx2->SetLineColor(kBlue);
+ hEx2->Draw("same");
+
+}
+
+void GateGamma_SeeParticle_LoadFromFile(double gamma, double width, double gamma2, double width2, double gamma3, double width3, double gammaBin, double particleBin){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+ TCanvas *cGGSP = new TCanvas("cGGSP","cGGSP",1000,1000);
+ cGGSP->Divide(1,2);
+ double zoom = 0.2;
+
+ cGGSP->cd(1);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",0,10000);
+ hEg->GetXaxis()->SetRangeUser(gamma-zoom, gamma+zoom);
+ hEg->SetLineColor(kBlack);
+ hEg->Draw();
+ double height = hEg->GetMaximum();
+ Func_AddGatingLines(gamma, width, height, kRed);
+ Func_AddGatingLines(gamma2, width2, height, kBlue);
+ Func_AddGatingLines(gamma3, width3, height, kGreen);
+
+ cGGSP->cd(2);
+ TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",
+ h2D->GetXaxis()->FindBin(gamma-width),
+ h2D->GetXaxis()->FindBin(gamma+width)
+ );
+ hEx->SetLineColor(kRed);
+ hEx->Draw();
+
+ TH1F* hEx2 = (TH1F*)h2D->ProjectionY("hEx2",
+ h2D->GetXaxis()->FindBin(gamma2-width2),
+ h2D->GetXaxis()->FindBin(gamma2+width2)
+ );
+ hEx2->SetLineColor(kBlue);
+ hEx2->Draw("same");
+
+ TH1F* hEx3 = (TH1F*)h2D->ProjectionY("hEx3",
+ h2D->GetXaxis()->FindBin(gamma3-width3),
+ h2D->GetXaxis()->FindBin(gamma3+width3)
+ );
+ hEx3->SetLineColor(kGreen);
+ hEx3->Draw("same");
+
+}
+
+void GateParticle_SeeGamma_LoadFromFile(double gammaBin, double particleBin){
+
+ //vector<double> edges = { 3.0, 3.4, 3.8, 4.2, 4.6, 5.0 };
+ //vector<double> edges = { 3.0, 3.2, 3.4, 3.6, 3.8, 4.0 };
+ //vector<double> edges = { 4.0, 4.1, 4.2, 4.3 };
+ vector<double> edges = { 4.065, 4.195, 4.285, 4.375 };
+
+ int numSlices = edges.size()-1;
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000);
+ cGPSG->Divide(1,2);
+
+ cGPSG->cd(1);
+ TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",0,6000);
+ hEx->SetLineColor(kBlack);
+ hEx->Draw();
+ double height = hEx->GetMaximum();
+
+ for(int i = 0; i < numSlices; i++){
+ double width = (edges[i+1]-edges[i])/2.0;
+ Func_AddGatingLines(edges[i]+width, width, height, (EColor)(i+2));
+ }
+
+ cGPSG->cd(2);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",
+ h2D->GetYaxis()->FindBin(edges.front()),
+ h2D->GetYaxis()->FindBin(edges.back())
+ );
+ hEg->SetLineColor(kBlack);
+
+ vector<TH1F*> hEgs;
+ for(int i = 0; i < numSlices; i++){
+ string name = "hEgs"+to_string(i);
+
+ TH1F* hEgTemp = (TH1F*)h2D->ProjectionX(name.c_str(),
+ h2D->GetYaxis()->FindBin(edges[i]),
+ h2D->GetYaxis()->FindBin(edges[i+1])
+ );
+ hEgTemp->SetLineColor(i+2);
+ hEgs.push_back(hEgTemp);
+ }
+
+ hEg->Draw();
+
+ for(int i = 0; i < numSlices; i++){
+ hEgs[i]->Draw("same");
+ }
+
+}
+
+void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double gammaBin, double particleBin){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000);
+ cGPSG->Divide(1,2);
+ double zoom = 0.2;
+
+ cGPSG->cd(1);
+ TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",0,6000);
+ hEx->SetLineColor(kBlack);
+ hEx->Draw();
+ double height = hEx->GetMaximum();
+ Func_AddGatingLines(particle, width, height, kRed);
+
+
+ cGPSG->cd(2);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",
+ h2D->GetYaxis()->FindBin(particle-width),
+ h2D->GetYaxis()->FindBin(particle+width)
+ );
+ hEg->SetLineColor(kBlack);
+ hEg->Draw();
+ double gamheight = hEg->GetMaximum();
+
+ AddGammaLines(hEg, particle, gamheight);
+
+}
+
+void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double particle2, double width2, double gammaBin, double particleBin){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000);
+ cGPSG->Divide(1,2);
+
+ cGPSG->cd(1);
+ TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",0,6000);
+ hEx->SetLineColor(kBlack);
+ hEx->Draw();
+ double height = hEx->GetMaximum();
+ Func_AddGatingLines(particle, width, height, kRed);
+ Func_AddGatingLines(particle2, width2, height, kBlue);
+
+
+ cGPSG->cd(2);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",
+ h2D->GetYaxis()->FindBin(particle-width),
+ h2D->GetYaxis()->FindBin(particle+width)
+ );
+ hEg->SetLineColor(kRed);
+
+ TH1F* hEg2 = (TH1F*)h2D->ProjectionX("hEg2",
+ h2D->GetYaxis()->FindBin(particle2-width2),
+ h2D->GetYaxis()->FindBin(particle2+width2)
+ );
+ hEg2->SetLineColor(kBlue);
+
+
+ hEg->Draw();
+ hEg2->Draw("same");
+}
+
+void GateParticle_SeeGamma_LoadFromFile(double particle, double width, double particle2, double width2, double particle3, double width3, double gammaBin, double particleBin){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DPartGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DPartGamma(h2D, gammaBin, particleBin);
+
+ TCanvas *cGPSG = new TCanvas("cGPSG","cGPSG",1000,1000);
+ cGPSG->Divide(1,2);
+
+ cGPSG->cd(1);
+ TH1F* hEx = (TH1F*)h2D->ProjectionY("hEx",0,6000);
+ hEx->SetLineColor(kBlack);
+ hEx->Draw();
+ double height = hEx->GetMaximum();
+ Func_AddGatingLines(particle, width, height, kRed);
+ Func_AddGatingLines(particle2, width2, height, kBlue);
+ Func_AddGatingLines(particle3, width3, height, kGreen);
+
+ cGPSG->cd(2);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",
+ h2D->GetYaxis()->FindBin(particle-width),
+ h2D->GetYaxis()->FindBin(particle+width)
+ );
+ hEg->SetLineColor(kRed);
+
+ TH1F* hEg2 = (TH1F*)h2D->ProjectionX("hEg2",
+ h2D->GetYaxis()->FindBin(particle2-width2),
+ h2D->GetYaxis()->FindBin(particle2+width2)
+ );
+ hEg2->SetLineColor(kBlue);
+
+ TH1F* hEg3 = (TH1F*)h2D->ProjectionX("hEg3",
+ h2D->GetYaxis()->FindBin(particle3-width3),
+ h2D->GetYaxis()->FindBin(particle3+width3)
+ );
+ hEg3->SetLineColor(kGreen);
+
+ hEg->Draw();
+ hEg2->Draw("same");
+ hEg3->Draw("same");
+}
+
+void GateGamma_SeeGamma_LoadFromFile(double gamma, double width, double gammaBin1, double gammaBin2){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DGammaGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DGammaGamma(h2D, gammaBin1, gammaBin2);
+
+ TCanvas *cGGSG = new TCanvas("cGGSG","cGGSG",1000,1000);
+ cGGSG->Divide(1,2);
+ double zoom = 0.2;
+
+ cGGSG->cd(1);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",0,10000);
+ hEg->GetXaxis()->SetRangeUser(gamma-zoom, gamma+zoom);
+ hEg->SetLineColor(kBlack);
+ hEg->Draw();
+ double height = hEg->GetMaximum();
+ Func_AddGatingLines(gamma, width, height, kRed);
+
+ cGGSG->cd(2);
+ TH1F* hEgg = (TH1F*)h2D->ProjectionY("hEgg",
+ h2D->GetXaxis()->FindBin(gamma-width),
+ h2D->GetXaxis()->FindBin(gamma+width)
+ );
+ hEgg->SetLineColor(kBlack);
+ hEgg->Draw();
+
+
+}
+
+void GateGamma_SeeGamma_LoadFromFile(double gamma, double width, double gamma2, double width2, double gammaBin1, double gammaBin2){
+
+ TH2F* h2D = (TH2F*)Func_LoadIn2DGammaGamma();
+ gStyle->SetOptStat(0);
+ h2D = Func_Rebin2DGammaGamma(h2D, gammaBin1, gammaBin2);
+
+ TCanvas *cGGSG = new TCanvas("cGGSG","cGGSG",1000,1000);
+ cGGSG->Divide(1,2);
+ double zoom = 0.2;
+
+ cGGSG->cd(1);
+ TH1F* hEg = (TH1F*)h2D->ProjectionX("hEg",0,10000);
+ hEg->GetXaxis()->SetRangeUser(gamma-zoom, gamma+zoom);
+ hEg->SetLineColor(kBlack);
+ hEg->Draw();
+ double height = hEg->GetMaximum();
+ Func_AddGatingLines(gamma, width, height, kRed);
+ Func_AddGatingLines(gamma2, width2, height, kBlue);
+
+ cGGSG->cd(2);
+ TH1F* hEgg = (TH1F*)h2D->ProjectionY("hEgg",
+ h2D->GetXaxis()->FindBin(gamma-width),
+ h2D->GetXaxis()->FindBin(gamma+width)
+ );
+ hEgg->SetLineColor(kRed);
+
+ TH1F* hEgg2 = (TH1F*)h2D->ProjectionY("hEgg2",
+ h2D->GetXaxis()->FindBin(gamma2-width2),
+ h2D->GetXaxis()->FindBin(gamma2+width2)
+ );
+ hEgg2->SetLineColor(kBlue);
+
+ hEgg->Draw();
+ hEgg2->Draw("same");
+}
+
+/*
+void ExRangesCompareGammas(){
+ TCanvas* cERCG = new TCanvas("cERCG","cERCG",1000,1000);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.03);
+
+ cERCG->Divide(2,3);
+
+ cERCG->cd(1);
+ chain->Draw("AddBack_EDC>>g01(1000,0,5)",
+ "abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && Mugast.TelescopeNumber<8 && abs(Ex-0.5)<0.5",
+ "");
+ TH2F* g01 = (TH2F*) gDirectory->Get("g01");
+
+ cERCG->cd(2);
+ chain->Draw("AddBack_EDC>>g12(1000,0,5)",
+ "abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && Mugast.TelescopeNumber<8 && abs(Ex-1.5)<0.5",
+ "");
+ TH2F* g12 = (TH2F*) gDirectory->Get("g12");
+
+ cERCG->cd(3);
+ chain->Draw("AddBack_EDC>>g23(1000,0,5)",
+ "abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && Mugast.TelescopeNumber<8 && abs(Ex-2.5)<0.5",
+ "");
+ TH2F* g23 = (TH2F*) gDirectory->Get("g23");
+
+ cERCG->cd(4);
+ chain->Draw("AddBack_EDC>>g34(1000,0,5)",
+ "abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && Mugast.TelescopeNumber<8 && abs(Ex-3.5)<0.5",
+ "");
+ TH2F* g34 = (TH2F*) gDirectory->Get("g34");
+
+ cERCG->cd(5);
+ chain->Draw("AddBack_EDC>>g45(1000,0,5)",
+ "abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && Mugast.TelescopeNumber<8 && abs(Ex-4.5)<0.5",
+ "");
+ TH2F* g45 = (TH2F*) gDirectory->Get("g45");
+
+ cERCG->cd(6);
+ chain->Draw("Ex>>gEx(800,-1,7)",
+ "abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && Mugast.TelescopeNumber<8",
+ "");
+ TH2F* gEx = (TH2F*) gDirectory->Get("gEx");
+
+}
+*/
+
+void ExRangesCompareGammas(){
+ TCanvas* cERCG = new TCanvas("cERCG","cERCG",1000,1000);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.03);
+
+ int s1 = 80, s2 = 60, s3 = 40, s4 = 20;
+
+ string gateBase = timegate
+ + " && " + det_gate
+ + " && Ex@.size()==1";
+ if(reactionName=="47K(d,t)"){
+ gateBase = gateBase + " && cutTritons && cutTime";
+ }
+
+ string gate01 = gateBase + "&& abs(Ex-0.5)<0.5";
+ chain->Draw("AddBack_EDC>>g01(2000,0,5)",gate01.c_str(),"");
+ TH1F* g01 = (TH1F*) gDirectory->Get("g01");
+
+ string gate12 = gateBase + "&& abs(Ex-1.5)<0.5";
+ chain->Draw("AddBack_EDC>>g12(2000,0,5)",gate12.c_str(),"");
+ TH1F* g12 = (TH1F*) gDirectory->Get("g12");
+
+ string gate23 = gateBase + "&& abs(Ex-2.5)<0.5";
+ chain->Draw("AddBack_EDC>>g23(2000,0,5)",gate23.c_str(),"");
+ TH1F* g23 = (TH1F*) gDirectory->Get("g23");
+
+ string gate34 = gateBase + "&& abs(Ex-3.5)<0.5";
+ chain->Draw("AddBack_EDC>>g34(2000,0,5)",gate34.c_str(),"");
+ TH1F* g34 = (TH1F*) gDirectory->Get("g34");
+
+ string gate45 = gateBase + "&& abs(Ex-4.5)<0.5";
+ chain->Draw("AddBack_EDC>>g45(2000,0,5)",gate45.c_str(),"");
+ TH1F* g45 = (TH1F*) gDirectory->Get("g45");
+
+ TH1F *t01 = (TH1F*)g01->Clone("t01");
+ TH1F *t12 = (TH1F*)g01->Clone("t12");
+ TH1F *t23 = (TH1F*)g01->Clone("t23");
+ TH1F *t34 = (TH1F*)g01->Clone("t34");
+
+ for(int i=0; i<g01->GetNbinsX()+1; i++){
+ t01->SetBinContent(i,s1);
+ t12->SetBinContent(i,s2);
+ t23->SetBinContent(i,s3);
+ t34->SetBinContent(i,s4);
+ }
+
+ g01->Add(t01,1);
+ g01->SetLineColor(kBlack);
+ g01->Draw();
+ TLine *l01 = new TLine(0.0, s1, 5.0, s1);
+ l01->SetLineColor(kBlack);
+ l01->SetLineStyle(kDotted);
+ l01->Draw();
+
+
+ g12->Add(t12,1);
+ g12->SetLineColor(kRed);
+ g12->Draw("SAME");
+ TLine *l12 = new TLine(0.0, s2, 5.0, s2);
+ l12->SetLineColor(kRed);
+ l12->SetLineStyle(kDotted);
+ l12->Draw();
+
+ g23->Add(t23,1);
+ g23->SetLineColor(kBlue);
+ g23->Draw("SAME");
+ TLine *l23 = new TLine(0.0, s3, 5.0, s3);
+ l23->SetLineColor(kBlue);
+ l23->SetLineStyle(kDotted);
+ l23->Draw();
+
+ g34->Add(t34,1);
+ g34->SetLineColor(kGreen);
+ g34->Draw("SAME");
+ TLine *l34 = new TLine(0.0, s4, 5.0, s4);
+ l34->SetLineColor(kGreen);
+ l34->SetLineStyle(kDotted);
+ l34->Draw();
+
+ g45->SetLineColor(kMagenta);
+ g45->Draw("SAME");
+ TLine *l45 = new TLine(0.0, 0.0, 5.0, 0.0);
+ l45->SetLineColor(kMagenta);
+ l45->SetLineStyle(kDotted);
+ l45->Draw();
+
+ g01->GetYaxis()->SetRangeUser(0,100);
+
+ auto legend = new TLegend(0.2, 0.2, .8, .8);
+ legend->AddEntry(g01, "0<Ex<1", "l");
+ legend->AddEntry(g12, "1<Ex<2", "l");
+ legend->AddEntry(g23, "2<Ex<3", "l");
+ legend->AddEntry(g34, "3<Ex<4", "l");
+ legend->AddEntry(g45, "4<Ex<5", "l");
+ legend->Draw();
+}
+
+void RandomStateSubtraction(double highlightMe){
+ bool wantHighlight = true;
+ int count = 0;
+ vector<double> temp;
+ if(reactionName=="47K(d,p)"){
+ temp = means;
+ } else if (reactionName=="47K(d,t)"){
+ temp = means_dt;
+ }
+
+ if(highlightMe<0.){wantHighlight=false;}
+
+ for(int i=temp.size(); i>=0; i--){
+ cout << "~~~ " << i << " ~~~" << endl;
+ for(int j=i; j>0; j--){
+ double val = temp.at(i-1) - temp.at(j-1);
+ if(abs(highlightMe-val)<0.010 && wantHighlight){cout << BOLDGREEN; count=+1;}
+ cout << temp.at(i-1) << " - " << temp.at(j-1) << " = " << val << endl;
+ if(abs(highlightMe-val)<0.010 && wantHighlight){cout << RESET;}
+ }
+ }
+
+ if(wantHighlight){cout << "FOUND " << count << " MATCH(ES)!" << endl;}
+}
+
+void RandomStateSubtraction(){
+ RandomStateSubtraction(-500.);
+}
+
+
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index f4f66e89c990396dd31411fc02525c27dca15d35..ed655dd4a4675d118352fe8c0f1db7e8957af650 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -3,40 +3,52 @@
#include <cmath>
#include "stdlib.h"
-const int numPeaks = 17;
-array<double,numPeaks> means = { 0.000,
+vector<double> means = { 0.000,
0.143,
0.279,
0.728,
- 0.968,
- 1.410,
- 1.981,//1.952,//1.981,
- 2.412,
- 2.910,
- 3.253,
- 3.605,
- 3.795,//Split in two?
- 3.870,//Split in two?
- 4.045,//4.1,
- 4.393,
- //4.51//,
- 5.15,
- 5.82
+ 0.967,
+ 1.409,
+ 1.978,//1.952,//1.981,
+ 2.407,
+ 2.909,
+ 3.254,
+ //3.455,//Possible???? TEMP!!!
+ 3.601,
+ //3.794,//Split in two?
+ //3.868,//Split in two?
+ 3.830,
+ //4.0,
+ 4.18,
+ //4.2,
+ 4.3,
+ 4.4
+ //4.16,
+ //4.23,
+ ////4.316,
+ //4.387
+ ////4.317,
+ ////4.42
};
+const int numPeaks = means.size();
-const int numPeaks_dt = 8;//9;
-array<double,numPeaks_dt> means_dt = { 0.000,
- //0.587, //from lit
+
+//const int numPeaks_dt = 8;//9;
+//array<double,numPeaks_dt> means_dt = { 0.000,
+vector<double> means_dt = { 0.000,
+ 0.587, //from lit
0.691, //from lit
//0.886, //from lit
- 1.738, //from lit
- 1.945,
- 2.26 , //from lit
- 2.76 , //from lit
+ //1.738, //from lit
+ 1.944,
+ 2.233,
+ 2.732,
3.344,
- 4.2
+ 3.410,
+ 4.297
};
+const int numPeaks_dt = means_dt.size();
array<double,27> knowngammas = { 0.143,
@@ -88,17 +100,22 @@ Double_t f_peak(Double_t *x, Double_t *par){
}
*/
+string f_full_string(){
+ string result = "gaus(0)";
+ for(int pk=1; pk<numPeaks; pk++){
+ result += "+gaus(" + to_string(3*pk) + ")";
+ }
+ return result;
+}
+
Double_t f_full(Double_t *x, Double_t *par) {
float xx = x[0];
double result, norm;
- // Flat background
-// result = par[0];
result = 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));
- * exp(-0.5*pow((xx-par[2+(pk*3)]-par[0])/par[1+(pk*3)],2)); //added par 0 as shift in energy
+ result += (par[1+(pk*3)]*exp(-0.5*pow(((xx-par[2+(pk*3)]-par[0])/par[3+(pk*3)]),2))); //REGULAR GAUS
}
return result;
}
@@ -121,259 +138,292 @@ Double_t f_full_dt(Double_t *x, Double_t *par) {
vector<vector<double>> FitKnownPeaks_RtrnArry(TH1F* hist, double slideshift){
double minFit=-1.0, maxFit=8.0;
double binWidth = hist->GetXaxis()->GetBinWidth(3);
- double sigma = 0.14;
-
- 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(), minFit, maxFit);
- //allPeaks[i] = new TF1(nameHere.c_str(), f_peak, -1, 5);
- allPeaks[i]->SetLineColor(kBlack);
- allPeaks[i]->SetLineStyle(7);
- allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth");
- }
-
- //Subtract flat background equal to smallest bin in range
-
- hist->GetXaxis()->SetRange(hist->FindBin(-0.9), hist->FindBin(-0.4));
- double bgmin = hist->GetBinContent(hist->GetMinimumBin());
- hist->GetXaxis()->UnZoom();
- cout << "Subtracting background of " << bgmin << endl;
- for(int b=1; b<hist->GetNbinsX() ; b++){
- hist->SetBinContent(b,hist->GetBinContent(b)-bgmin);
- }
-
-
- //Build background function
- TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
- bg->SetLineColor(kGreen);
- bg->SetLineStyle(9);
- bg->SetParNames("Background");
-
- //Build IMPROVED total function
- 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((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->SetParLimits((i*3)+3,10.0,1e5);
- }
- //full->SetParLimits(0,0.,40.); /* FOR TOTAL SPECTRUM FITTING */
- //full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
- //full->FixParameter(0,0.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
- full->FixParameter(9,0.); /* FIX 0.279 AREA TO ZERO */
- //full->SetParLimits(0,-0.5,+0.5); /* FOR WHEN PAR[0] IS VARIABLE ENERGY CENTRIOD SLIDER */
- //full->FixParameter(0,slideshift); /* FOR WHEN PAR[0] IS FIXED ENERGY CENTRIOD SLIDER */
- full->FixParameter(0,0.0); /* FOR WHEN FIXING ENERGY SLIDER TO ZERO */
- full->SetParameter(52,0.39); /* SET 5.8MeV SIGMA */
- full->SetParLimits(52,0.34,0.44); /* SET 5.8MeV SIGMA */
-
- allPeaks[14]->SetLineColor(kOrange);
-
- //Fit full function to histogram
- hist->Fit(full, "RWQB", "", minFit, maxFit);
- hist->Draw();
+ //double sigma = 0.14;
+
+ double sigma = 0.0;// = 0.14;
+ int numPeaks_temp = 0;
+ vector<double> means_temp;
+
+ //Gradient on mean value
+ double gradient;
+ if(reactionName=="47K(d,p)"){
+ sigma = 0.14;
+ gradient = 0.989;
+ numPeaks_temp = numPeaks;
+ means_temp = means;
+
+ } else if (reactionName=="47K(d,t)"){
+ sigma = 0.18;
+ gradient = 1.000;
+ numPeaks_temp = numPeaks_dt;
+ means_temp = means_dt;
- //Extract fitted variables, assign them to individual fits, and draw them
- 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[3+(i*3)]);
- }
- bg->SetParameter(0,finalPar[0]);
- bg->Draw("SAME");
- full->Draw("SAME");
-
- for (int i=0; i<numPeaks; i++){
- allPeaks[i]->Draw("SAME");
}
- /* Error propogation:
- * (Abin) +- deltaAbin, B+-0 (no uncertainty)
- * A = Abin/B
- * deltaA/A = deltaAbin/Abin
- * deltaA = A x deltaAbin/Abin
- */
-
- //Write to screen
- cout << "===========================" << endl;
- cout << "== PEAK =========== AREA ==" << endl;
-
- vector<vector<double>> allpeaks;
- for(int i=0; i<numPeaks; i++){
- double A = finalPar[(i*3)+3]/binWidth;
- double deltaA = A * (finalErr[(i*3)+3]/finalPar[(i*3)+3]);
-
- cout << fixed << setprecision(3)
- << " #" << i << " "
- << 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[(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;
-}
-
-vector<vector<double>> FitKnownPeaks_dt_RtrnArry(TH1F* hist, double slideshift){
- double minFit=-1.0, maxFit=8.0;
- double binWidth = hist->GetXaxis()->GetBinWidth(3);
- double sigma = 0.18;
+ cout << BOLDGREEN
+ << "FITTING WITH A REAL:OBSERVED Ex GRADIENT GRADIENT OF "
+ << gradient
+ << RESET
+ << endl;
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);
- **/
+ /** REMOVED FOR NOW **/
//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_dt];
- for(int i=0; i<numPeaks_dt; i++) {
+ string function = "gaus";//"[0]*exp(-0.5*pow(((x-[1])/[2]),2))";//REGULAR GAUS
+ TF1 **allPeaks = new TF1*[numPeaks_temp];
+ for(int i=0; i<numPeaks_temp; i++) {
string nameHere = nameBase;
nameHere +=to_string(i);
-
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);
- allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth");
+ allPeaks[i]->SetNpx(500);
}
//Subtract flat background equal to smallest bin in range
+ /** REMOVED FOR NOW **/
- hist->GetXaxis()->SetRange(hist->FindBin(-0.9), hist->FindBin(-0.4));
- double bgmin = hist->GetBinContent(hist->GetMinimumBin());
- hist->GetXaxis()->UnZoom();
- cout << "Subtracting background of " << bgmin << endl;
- for(int b=1; b<hist->GetNbinsX() ; b++){
- hist->SetBinContent(b,hist->GetBinContent(b)-bgmin);
- }
-
+ ////Gradient on mean value
+ //double gradient;
+ //gradient = 0.989;
+ //cout << BOLDGREEN
+ // << "FITTING WITH A REAL:OBSERVED Ex GRADIENT GRADIENT OF "
+ // << gradient
+ // << RESET
+ // << endl;
//Build background function
- TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
- bg->SetLineColor(kGreen);
- bg->SetLineStyle(9);
- bg->SetParNames("Background");
+ /** REMOVED FOR NOW **/
- //Build IMPROVED total function
- TF1 *full = new TF1("fitAllPeaks", f_full_dt, minFit, maxFit, (int) 1+(3*numPeaks_dt));
+ //Build total function
+ string s_full = f_full_string();
+ TF1 *full = new TF1("fitAllPeaks", s_full.c_str(), minFit, maxFit);//, (int) (3*numPeaks));
full->SetLineColor(kRed);
- const int numParams = (numPeaks_dt*3)+1;
- for(int i=0; i<numPeaks_dt; i++) {
- full->FixParameter((i*3)+1,sigma);
- full->FixParameter((i*3)+2,means_dt.at(i));
- full->SetParameter((i*3)+3,1e1);
- full->SetParLimits((i*3)+3,0.0,1e5);
- //full->SetParLimits((i*3)+3,10.0,1e5);
+ const int numParams = (numPeaks_temp*3);
+ for(int i=0; i<numPeaks_temp; i++) {
+ full->FixParameter((i*3)+1,means_temp.at(i)*gradient); //GRADIENT REMOVED FOR NOW
+ full->FixParameter((i*3)+2,sigma);
+ full->SetParameter((i*3)+0,1e1);
+ full->SetParLimits((i*3)+0,0.0,1e5);
}
- //full->SetParLimits(0,0.,40.); /* FOR TOTAL SPECTRUM FITTING */
- //full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
- //full->FixParameter(0,0.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
-
+ full->SetNpx(500);
+
+ //Set and fix specific parameters
+ /** REMOVED FOR NOW **/
+
//Fit full function to histogram
- hist->Fit(full, "RWQB", "", minFit, maxFit);
+ //hist->Fit(full, "RWQB", "", minFit, maxFit);
+ hist->Fit(full, "RWQBL", "", minFit, maxFit); cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl;
hist->Draw();
//Extract fitted variables, assign them to individual fits, and draw them
const Double_t* finalPar = full->GetParameters();
const Double_t* finalErr = full->GetParErrors();
- for (int i=0; i<numPeaks_dt; i++){
- allPeaks[i]->SetParameters(sigma, means_dt.at(i), finalPar[3+(i*3)]);
+ for (int i=0; i<numPeaks_temp; i++){
+ allPeaks[i]->SetParameters(finalPar[0+(i*3)], finalPar[1+(i*3)], finalPar[2+(i*3)]);
}
- bg->SetParameter(0,finalPar[0]);
- bg->Draw("SAME");
full->Draw("SAME");
- for (int i=0; i<numPeaks_dt; i++){
+ //Make irresolvable doublet clear
+ /** REMOVED FOR NOW **/
+
+ //Draw all peaks
+ for (int i=0; i<numPeaks_temp; i++){
allPeaks[i]->Draw("SAME");
}
- /* Error propogation:
- * (Abin) +- deltaAbin, B+-0 (no uncertainty)
- * A = Abin/B
- * deltaA/A = deltaAbin/Abin
- * deltaA = A x deltaAbin/Abin
- */
-
//Write to screen
cout << "===========================" << endl;
cout << "== PEAK =========== AREA ==" << endl;
+ //Loop over every peak
vector<vector<double>> allpeaks;
- for(int i=0; i<numPeaks_dt; i++){
- double A = finalPar[(i*3)+3]/binWidth;
- double deltaA = A * (finalErr[(i*3)+3]/finalPar[(i*3)+3]);
+ for(int i=0; i<numPeaks_temp; i++){
+ //For AREA = HEIGHT * SIGMA * SQRT(2*PI)
+ double A = (finalPar[0+(i*3)] * finalPar[2+(i*3)] * sqrt(2*pi)) / binWidth;
+// cout << "AREA = " << finalPar[0+(i*3)] << " * "
+// << finalPar[2+(i*3)] << " * "
+// << sqrt(2*pi) << " = "
+// << A << endl;
+
+ double deltaA = A * sqrt( pow( finalErr[0+(i*3)]/finalPar[0+(i*3)] ,2) + pow( finalErr[2+(i*3)]/finalPar[2+(i*3)] ,2) );
+ cout << "DELTAAREA = "
+ << finalErr[0+(i*3)] << " / "
+ << finalPar[0+(i*3)] << " and "
+ << finalErr[2+(i*3)] << " / "
+ << finalPar[2+(i*3)] << " = "
+ << deltaA << endl;
cout << fixed << setprecision(3)
<< " #" << i << " "
- << finalPar[(i*3)+2] << "\t" << setprecision(0)
+ << finalPar[(i*3)+1] << "\t" << setprecision(0)
<< A << "\t+- "
- << deltaA << setprecision(3)
- << endl;
+ << deltaA << setprecision(3);
+ cout << " SQRT: " << sqrt(A) << endl;
vector<double> onepeak; //energy, area and error for one peak
- onepeak.push_back(finalPar[(i*3)+2]);
+ onepeak.push_back(finalPar[(i*3)+1]);
onepeak.push_back(A);
onepeak.push_back(deltaA);
allpeaks.push_back(onepeak);
+ cout << "-------------" << endl;
}
- cout << " BG " << full->GetParameter(0)
- << " +- " << full->GetParError(0) << endl;
return allpeaks;
}
+////vector<vector<double>> FitKnownPeaks_dt_RtrnArry(TH1F* hist, double slideshift){
+//// double minFit=-1.0, maxFit=8.0;
+//// double binWidth = hist->GetXaxis()->GetBinWidth(3);
+//// double sigma = 0.18;
+////
+//// 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_dt];
+//// for(int i=0; i<numPeaks_dt; i++) {
+//// string nameHere = nameBase;
+//// nameHere +=to_string(i);
+////
+//// 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);
+//// allPeaks[i]->SetParNames("Sigma", "Mean", "Area*BinWidth");
+//// }
+////
+//// //Subtract flat background equal to smallest bin in range
+//// hist->GetXaxis()->SetRange(hist->FindBin(-0.8), hist->FindBin(-0.0));
+//// double bgmin = hist->GetBinContent(hist->GetMinimumBin());
+//// hist->GetXaxis()->UnZoom();
+//// cout << "Subtracting background of " << bgmin << endl;
+//// for(int b=1; b<hist->GetNbinsX() ; b++){
+//// hist->SetBinContent(b,hist->GetBinContent(b)-bgmin);
+//// }
+////
+//// ////Subtract flat background equal to average bin in range
+//// //double bgavg = 0;
+//// //int bgcount = 0;
+//// //for(int i = hist->FindBin(-0.8); i<hist->FindBin(0.1); i++){
+//// // bgavg += hist->GetBinContent(i);
+//// // bgcount++;
+//// //}
+//// //bgavg = bgavg/(double)bgcount;
+//// //hist->GetXaxis()->UnZoom();
+//// //cout << "Subtracting background of " << bgavg << endl;
+//// //for(int b=1; b<hist->GetNbinsX() ; b++){
+//// // hist->SetBinContent(b,hist->GetBinContent(b)-bgavg);
+//// //}
+////
+////
+//// //Build background function
+//// TF1 *bg = new TF1 ("bg","[0]",minFit, maxFit);
+//// bg->SetLineColor(kGreen);
+//// bg->SetLineStyle(9);
+//// bg->SetParNames("Background");
+////
+//// //Build IMPROVED total function
+//// TF1 *full = new TF1("fitAllPeaks", f_full_dt, minFit, maxFit, (int) 1+(3*numPeaks_dt));
+//// full->SetLineColor(kRed);
+//// const int numParams = (numPeaks_dt*3)+1;
+//// for(int i=0; i<numPeaks_dt; i++) {
+//// //full->FixParameter((i*3)+1,sigma);
+//// //full->FixParameter((i*3)+2,means_dt.at(i));
+//// //full->SetParameter((i*3)+3,1e1);
+//// //full->SetParLimits((i*3)+3,0.0,1e5);
+//// ////full->SetParLimits((i*3)+3,10.0,1e5);
+//// //
+//// //
+//// //
+//// //
+//// full->FixParameter((i*3)+1,means.at(i)); //GRADIENT REMOVED FOR NOW
+//// full->FixParameter((i*3)+2,sigma);
+//// full->SetParameter((i*3)+0,1e1);
+//// full->SetParLimits((i*3)+0,0.0,1e5);
+////
+////
+//// }
+//// //full->SetParLimits(0,0.,40.); /* FOR TOTAL SPECTRUM FITTING */
+//// //full->SetParLimits(0,0.,10.); /* FOR ANGLE GATED FITTING */
+//// //full->FixParameter(0,0.); /* FOR ANGLE GATED FITTING WITH BG SUBTRACTED */
+////
+//// //Fit full function to histogram
+//// //hist->Fit(full, "RWQB", "", minFit, maxFit);
+//// hist->Fit(full, "RWQBL", "", minFit, maxFit); cout << GREEN << "FITTING WITH MAX LIKELYHOOD METHOD" << RESET << endl;
+//// hist->Draw();
+////
+//// //Extract fitted variables, assign them to individual fits, and draw them
+//// const Double_t* finalPar = full->GetParameters();
+//// const Double_t* finalErr = full->GetParErrors();
+//// for (int i=0; i<numPeaks_dt; i++){
+//// allPeaks[i]->SetParameters(sigma, means_dt.at(i), finalPar[3+(i*3)]);
+//// }
+//// bg->SetParameter(0,finalPar[0]);
+//// bg->Draw("SAME");
+//// full->Draw("SAME");
+////
+//// for (int i=0; i<numPeaks_dt; i++){
+//// allPeaks[i]->Draw("SAME");
+//// }
+////
+//// /* Error propogation:
+//// * (Abin) +- deltaAbin, B+-0 (no uncertainty)
+//// * A = Abin/B
+//// * deltaA/A = deltaAbin/Abin
+//// * deltaA = A x deltaAbin/Abin
+//// */
+////
+//// //Write to screen
+//// cout << "===========================" << endl;
+//// cout << "== PEAK =========== AREA ==" << endl;
+////
+//// vector<vector<double>> allpeaks;
+//// for(int i=0; i<numPeaks_dt; i++){
+//// double A = finalPar[(i*3)+3]/binWidth;
+//// double deltaA = A * (finalErr[(i*3)+3]/finalPar[(i*3)+3]);
+////
+//// cout << fixed << setprecision(3)
+//// << " #" << i << " "
+//// << 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[(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;
+////}
+
void FitKnownPeaks(TH1F* hist){
//Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0);
@@ -381,6 +431,7 @@ void FitKnownPeaks(TH1F* hist){
void FitKnownPeaks_dt(TH1F* hist){
//Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
- vector<vector<double>> shell = FitKnownPeaks_dt_RtrnArry(hist,0.0);
+ //vector<vector<double>> shell = FitKnownPeaks_dt_RtrnArry(hist,0.0);
+ vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0);
}
diff --git a/Projects/e793s/macro/Plots_47Kdd.C b/Projects/e793s/macro/Plots_47Kdd.C
index 1e7e13c2ef312582629c522002942f20962db092..c656ffe7a60c2d9b5c72b2a26f4043f9f4c6553d 100644
--- a/Projects/e793s/macro/Plots_47Kdd.C
+++ b/Projects/e793s/macro/Plots_47Kdd.C
@@ -1,3 +1,4 @@
+string reactionName; /* defined by choice of dp or dt */
#include "DefineColours.h"
#include "GausFit.h"
#include "KnownPeakFitter.h"
diff --git a/Projects/e793s/macro/Plots_47Kdp.C b/Projects/e793s/macro/Plots_47Kdp.C
index ac9f2d4dc0a0afcce02ada0ffe8c5b173fdfed6f..db46a8bf1646e5add0463b12078340e1b4198ed7 100644
--- a/Projects/e793s/macro/Plots_47Kdp.C
+++ b/Projects/e793s/macro/Plots_47Kdp.C
@@ -1,19 +1,22 @@
+string reactionName; /* defined by choice of dp or dt */
#include "DefineColours.h"
#include "GausFit.h"
#include "KnownPeakFitter.h"
#include "DrawPlots.h"
-#include "CS2.h"
+#include "CS2_master.h"
+//#include "CS2.h"
//#include "CS2_MGX.h"
#include "ThreeBodyBreakup.h"
#include "ThreeBodyBreakup_FitPhaseSpace.h"
-#include "20Oct22_CompareYield.h"
+//#include "20Oct22_CompareYield.h"
void AddGammaLines(TH1F* hist, double particle, double ymax){
+
string base = "sub ";
- for(int i=1; i<means.size();i++){
+ for(int i=0; i<means.size();i++){
string name = base + to_string(means.at(i));
TLine *line = new TLine(particle-means.at(i), 0.0, particle-means.at(i), ymax);
line->SetLineColor(kBlack); line->SetLineStyle(kDotted);
@@ -97,18 +100,127 @@ void CompareCountsInThetaLab(){
}
+void CS_Diagnosis(){
+ auto majorCanv = new TCanvas("CompareCanv","CompareCanv",1500,1500);
+ majorCanv->Divide(3,3);
+ canclone(majorCanv, 1, "c_peakFits_110_112");
+ canclone(majorCanv, 2, "c_peakFits_115_118");
+ canclone(majorCanv, 3, "c_peakFits_120_122");
+ canclone(majorCanv, 4, "c_peakFits_125_128");
+ canclone(majorCanv, 5, "c_peakFits_130_132");
+ canclone(majorCanv, 6, "c_peakFits_135_138");
+ canclone(majorCanv, 7, "c_peakFits_140_142");
+ canclone(majorCanv, 8, "c_peakFits_145_148");
+ canclone(majorCanv, 9, "c_peakFits_150_152");
+}
+
+void CS(){
+// Overload function //
+ cout << "- CS(stateE, stateSp, orb_l, orb_j, binEx, options) " << endl;
+ cout << "|-----------------| GOOD |-----------------|" << endl;
+ cout << "---- 0.143, p3/2 = CS(0.143, 2, 1, 1.5, 0.05, \"\") " << endl;
+ cout << "---- 0.279, NOT POPULATED " << endl;
+ cout << "---- 0.728, WEAKLY POPULATED " << endl;
+ cout << "---- 0.967, p1/2 = CS(0.967, 0, 1, 0.5, 0.05, \"\") " << endl;
+ cout << "---- 1.409, p3/2 = CS(1.409, 1, 1, 1.5, 0.05, \"\") " << endl;
+ cout << "---- 1.978, p3/2 = CS(1.978, 1, 1, 0.5, 0.05, \"\") " << endl;
+ cout << "---- 2.407, p3/2 = CS(2.407, 0, 1, 0.5, 0.05, \"\") " << endl;
+ cout << "---- 2.908, p3/2 = CS(2.908, 2, 1, 1.5, 0.05, \"mixed\") " << endl;
+ cout << "---- 2.908, f5/2 = CS(2.908, 2, 3, 2.5, 0.05, \"mixed\") " << endl;
+ cout << "---- 3.254, f5/2 = CS(3.254, 3, 3, 2.5, 0.05, \"\") " << endl;
+ cout << "---- 3.601, f5/2 = CS(3.601, 3, 3, 2.5, 0.05, \"\") " << endl;
+ cout << endl;
+ cout << "---- 3.792 (x-?) " << endl;
+ cout << " FIT TOGETHER = CS(3.830, 2, 3, 2.5, 0.05, \"\")" << endl;
+ cout << "---- 3.868 (2-?) " << endl;
+ cout << endl;
+ cout << "|----------------| UNSURE |----------------|" << endl;
+ cout << "---- 4.061, f5/2 = CS(4.061, 3, 3, 2.5, 0.05, \"\") " << endl;
+ cout << "---- 4.387, f5/2 = CS(4.387, 3, 3, 2.5, 0.05, \"\") " << endl;
+}
+
+void CompareGammas_ParticleRegions_48K(double binwidth, double minEx, double maxEx, double stepEx){
+
+ vector<TH1F*> Eg;
+ vector<string> Names;
+ int numHists = (int) ((maxEx-minEx)/stepEx);
+ cout << numHists << endl;
+
+ string draw = "AddBack_EDC>>EgTemp(" + to_string((int) (5./binwidth)) + ",0,5)";
+ cout << draw << endl;
+
+ auto cTEMP = new TCanvas("cTEMP","cTEMP",500,500);
+ for(int i=0; i<numHists; i++){
+ cout << "====================================" << endl;
+ cout << i << endl;
+ double centre = minEx + (stepEx*(double)i) + (0.5*stepEx);
+
+ string gate = "Mugast.TelescopeNumber>0 && abs(T_MUGAST_VAMOS-2700)<400 && Ex@.size()==1 && abs(Ex-"
+ + to_string(centre)
+ + ")<"
+ + to_string(stepEx/2.0);
+
+ string nameTemp = "Ex = "
+ + to_string(centre-(0.5*stepEx))
+ + " to "
+ + to_string(centre+(0.5*stepEx));
+
+ string histTemp = "Ex"
+ + to_string(i);
+// + to_string((int)centre-(0.5*stepEx))
+// + "to"
+// + to_string(centre+(0.5*stepEx));
+
+ cout << gate << endl;
+
+ chain->Draw(draw.c_str(),gate.c_str(),"");
+ auto EgTemp = (TH1F*) gDirectory->Get("EgTemp");
+ EgTemp->SetLineColor(i+1);
+ EgTemp->SetNameTitle(histTemp.c_str(),histTemp.c_str());
+ Eg.push_back(EgTemp);
+ Names.push_back(nameTemp);
+
+ }
+ delete cTEMP;
+
+ cout << "out" << endl;
+
+
+ auto canv = new TCanvas("cCompareGammas","cCompareGammas",1000,1000);
+ Eg[0]->Draw();
+ for (int i=1; i<numHists; i++){
+ Eg[i]->Draw("same");
+ }
+
+ cout << "out2" << endl;
+
+ auto legend = new TLegend(0.1,0.7,0.48,0.9);
+ for (int i=0; i<numHists; i++){
+ cout << "----------" << endl;
+ cout << i << endl;
+ legend->AddEntry(Eg[i],Names[i].c_str(),"f");
+ }
+ legend->Draw("same");
+
+}
+
+
/* MAIN FUNCTION */
void Plots_47Kdp(){
+cout << "test" << endl;
+
LoadChain47Kdp();
gStyle->SetOptStat("nemMrRi");
+cout << "test" << endl;
tCentre = 2700; tRange = 400;
timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange);
det_gate = "Mugast.TelescopeNumber>0 && Mugast.TelescopeNumber<8";
reactionName = "47K(d,p)";
+cout << "test" << endl;
cout << "==============================================" << endl;
cout << "=============== (d,p) reaction ===============" << endl;
cout << "==============================================" << endl;
diff --git a/Projects/e793s/macro/Plots_47Kdt.C b/Projects/e793s/macro/Plots_47Kdt.C
index 74e7d97e9a489f4f3f769289951789ff45bef11f..6efbad443cea0a2b9dc03bc4b6fa6a7333222814 100644
--- a/Projects/e793s/macro/Plots_47Kdt.C
+++ b/Projects/e793s/macro/Plots_47Kdt.C
@@ -1,9 +1,11 @@
+string reactionName; /* defined by choice of dp or dt */
#include "DefineColours.h"
#include "GausFit.h"
#include "KnownPeakFitter.h"
#include "DrawPlots.h"
-#include "CS2_dt.h"
+//#include "CS2_dt.h"
+#include "CS2_master.h"
//#include "ThreeBodyBreakup.h"
//#include "ThreeBodyBreakup_FitPhaseSpace.h"
@@ -32,6 +34,53 @@ void AddPlacedGammas(TH1F* hist, double ymax){
// }
}
+void Figure_ELabThetaLabAll(){
+
+ TCanvas* cELabTLab = new TCanvas("cELabTLab","cELabTLab",1000,1000);
+
+ chain->Draw("ELab:ThetaLab>>kEl(360,0,180,500,0,10)","abs(T_MUGAST_VAMOS-2750)<200 && MUST2.TelescopeNumber==5","colz");
+ TH2F* kEl = (TH2F*) gDirectory->Get("kEl");
+ kEl->SetTitle("");
+ kEl->GetXaxis()->SetTitle("#theta_{lab} [deg]");
+ kEl->GetYaxis()->SetTitle("E_{lab} [MeV]");
+
+
+ chain->Draw("ELab:ThetaLab>>kdp(360,0,180,500,0,10)","abs(T_MUGAST_VAMOS-2700)<400 && Mugast.TelescopeNumber>0","colz same");
+ TH2F* kdp = (TH2F*) gDirectory->Get("kdp");
+
+ chain->Draw("ELab:ThetaLab>>kdt(360,0,180,500,0,10)","abs(T_MUGAST_VAMOS-2750)<350 && MUST2.TelescopeNumber<5 && cutTime && cutTritons","colz same");
+ TH2F* kdt = (TH2F*) gDirectory->Get("kdt");
+
+ kEl->Draw("");
+ kdp->Draw("same");
+ kdt->Draw("same");
+
+}
+
+void CS_Diagnosis(){
+ auto majorCanv = new TCanvas("CompareCanv","CompareCanv",1500,1500);
+ majorCanv->Divide(3,3);
+ canclone(majorCanv, 1, "c_peakFits_110_112");
+ canclone(majorCanv, 2, "c_peakFits_115_118");
+ canclone(majorCanv, 3, "c_peakFits_120_122");
+ canclone(majorCanv, 4, "c_peakFits_125_128");
+ canclone(majorCanv, 5, "c_peakFits_130_132");
+ canclone(majorCanv, 6, "c_peakFits_135_138");
+ canclone(majorCanv, 7, "c_peakFits_140_142");
+ canclone(majorCanv, 8, "c_peakFits_145_148");
+ canclone(majorCanv, 9, "c_peakFits_150_152");
+}
+
+void CS(){
+/* Overload function */
+ cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
+ cout << "---- 0.691, f7/2 = CS(0.691, 4, 3, 3.5) "<< endl;
+ cout << "---- '' s1/2 = CS(0.691, 1, 0, 0.5) "<< endl;
+ cout << "---- '' d3/2 = CS(0.691, 1, 2, 1.5) "<< endl;
+ cout << "---- 1.945, s1/2 = CS(1.945, 1, 0, 0.5) "<< endl;
+ cout << "---- 3.344, s1/2 = CS(3.344, 1, 0, 0.5) "<< endl;
+}
+
/* MAIN FUNCTION */
void MM_Timing_Comparison(){
diff --git a/Projects/e793s/macro/Plots_47Kpp.C b/Projects/e793s/macro/Plots_47Kpp.C
new file mode 100644
index 0000000000000000000000000000000000000000..1f11fa6395547f22d0f5be44ab73fa41aae71746
--- /dev/null
+++ b/Projects/e793s/macro/Plots_47Kpp.C
@@ -0,0 +1,53 @@
+string reactionName; /* defined by choice of dp or dt */
+#include "DefineColours.h"
+#include "GausFit.h"
+#include "KnownPeakFitter.h"
+#include "DrawPlots.h"
+#include "ElasticsFitELabGates.h"
+
+//#include "CS2.h"
+//#include "ThreeBodyBreakup.h"
+//#include "ThreeBodyBreakup_FitPhaseSpace.h"
+
+void AddGammaLines(TH1F* hist, double particle, double ymax){
+// string base = "sub ";
+//
+// for(int i=1; i<means.size();i++){
+// string name = base + to_string(means.at(i));
+// TLine *line = new TLine(particle-means.at(i), 0.0, particle-means.at(i), ymax);
+// line->SetLineColor(kBlack); line->SetLineStyle(kDotted);
+// line->Draw();
+// TText *text = new TText((1.-(means.at(i)/particle))*particle,0.8*ymax,name.c_str());
+// text->SetTextAngle(90);
+// //text->SetTextSize(40);
+// text->Draw();
+// }
+}
+
+void AddPlacedGammas(TH1F* hist, double ymax){
+// hist->Draw();
+// for(int i=0; i<knowngammas.size();i++){
+// TLine *line = new TLine(knowngammas.at(i), 0.0, knowngammas.at(i), ymax);
+// line->SetLineColor(kBlack); line->SetLineStyle(kDotted);
+// line->Draw();
+// }
+}
+
+
+/* MAIN FUNCTION */
+
+void Plots_47Kpp(){
+
+ LoadChain47Kpp();
+ gStyle->SetOptStat("nemMrRi");
+
+ tCentre = 2750; tRange = 200;
+ timegate = "abs(T_MUGAST_VAMOS-" + to_string(tCentre) + ")<" + to_string(tRange);
+ det_gate = "MUST2.TelescopeNumber==5";
+ reactionName = "47K(p,p)";
+
+ cout << "==============================================" << endl;
+ cout << "=============== (p,p) reaction ===============" << endl;
+ cout << "==============================================" << endl;
+
+}