Miniclean annotated analysis script
Jump to navigation
Jump to search
Miniclean annotated analysis code by Nathan
annotated_analysis_script.C
//These #include statements load packages that will be used later on in the script//
#include <string>
#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH3D.h"
#include "TCanvas.h"
#include "RAT/DS/Root.hh"
#include "RAT/DS/MCParticle.hh"
#include "RAT/DS/MCTrack.hh"
#include "RAT/DS/MCTrackStep.hh"
//initializes the standard use for these commands and the string type//
using std::cout;
using std::endl;
using std::string;
//specifies rat namespace//
using namespace RAT::DS;
void makehistogram(const std::string &outFileName="weighted_hists.root")
{
//Making ambe spectrum for weighting (each number is the number of neutrons emitted at it's given energy level)://
Double_t ambe_mag[] = {0.000000 ,13.000000 ,14.000000 ,13.500000 ,12.000000 ,10.000000 ,8.800000 ,8.000000 ,7.800000 ,7.000000 ,6.500000 ,6.300000 ,6.000000 ,5.800000 ,5.200000 ,4.500000 ,4.300000 ,5.500000 ,6.300000 ,7.000000 ,7.500000 ,7.600000 ,7.300000 ,7.100000 ,7.800000 ,10.000000 ,11.000000 ,14.000000 ,16.000000 ,18.500000 ,18.000000 ,16.500000 ,15.500000 ,15.000000 ,14.300000 ,14.000000 ,13.800000 ,13.600000 ,13.300000 ,13.200000 ,13.200000 ,13.700000 ,14.000000 ,14.500000 ,15.000000 ,15.500000 ,15.300000 ,14.800000 ,14.000000 ,13.000000 ,12.500000 ,12.300000 ,11.800000 ,11.000000 ,10.000000 ,8.500000 ,7.800000 ,7.700000 ,8.000000 ,8.200000 ,7.900000 ,7.700000 ,6.600000 ,6.300000 ,6.100000 ,6.200000 ,6.500000 ,7.000000 ,7.500000 ,8.000000 ,8.300000 ,8.350000 ,8.200000 ,7.800000 ,7.300000 ,6.800000 ,5.800000 ,5.000000 ,4.200000 ,3.500000 ,2.800000 ,2.400000 ,1.800000 ,1.700000 ,1.800000 ,2.200000 ,2.700000 ,3.100000 ,3.500000 ,3.700000 ,3.500000 ,3.300000 ,3.100000 ,2.800000 ,2.700000 ,2.300000 ,1.800000 ,1.600000 ,1.100000 ,0.800000 ,0.400000 ,0.000000};
//Summing up array above to get the integral of the spectrum//
Double_t mag_tot = 0.0;
for(int i=0;i<102;i++){
mag_tot += ambe_mag[i];
}
//Dividing each spectrum element by integral to get the normalized weight for each energy//
for(int j=0;j<102;j++){
ambe_mag[j]*=1/mag_tot;
}
//Create new Chain to group all outputs together//
TChain* T = new TChain("T");
//Add all output files to chain in order to access them all//
T->Add("1400_pos_0.root");
T->Add("1400_pos_1.root");
T->Add("1400_pos_2.root");
T->Add("1400_pos_3.root");
T->Add("1400_pos_4.root");
T->Add("1400_pos_5.root");
T->Add("1400_pos_6.root");
T->Add("1400_pos_7.root");
T->Add("1400_pos_8.root");
T->Add("1400_pos_9.root");
T->Add("1400_pos_10.root");
T->Add("1400_pos_11.root");
T->Add("1400_pos_12.root");
T->Add("1400_pos_13.root");
T->Add("1400_pos_14.root");
T->Add("1400_pos_15.root");
T->Add("1400_pos_16.root");
T->Add("1400_pos_17.root");
T->Add("1400_pos_18.root");
T->Add("1400_pos_19.root");
//Specifiy path to the root tree//
RAT::DS::Root *rds = new RAT::DS::Root();
T->SetBranchAddress("ds", &rds);
//Initialize output file//
TFile *outFile=new TFile(outFileName.c_str(), "RECREATE");
//Initialize all desired histograms (some are commented out to be accessible for future use)//
TH1F *totalPE = new TH1F("totalPE","Total Photoelectrons", 100.0,0.0,550.0);
TH1F *totalPE_weight = new TH1F("totalPE_weight","Total Photoelectrons", 100.0,0.0,550.0);
TH1F *radialPos = new TH1F("radialPos","Event Distance Away From Origin",100.0,0.0,600.0);
TH1F *radialPos_weight = new TH1F("radialPos_weight","Event Distance Away From Origin",100.0,0.0,600.0);
TH1F *kE = new TH1F("kE","Kinetic Energy",100.0,0.0,12.0);
TH1F *kE_weight = new TH1F("kE_weight","Kinetic Energy",100.0,0.0,12.0);
//TH1F *lrecoil = new TH1F("lrecoil","lrecoil",100.0,0.0,1.0);
//TH1F *lrecoil_weight = new TH1F("lrecoil","lrecoil",100.0,0.0,1.0);
TH1F *fprompt = new TH1F("Fprompt","Fprompt",100.0,0.0,1.1);
TH1F *fprompt_weight = new TH1F("Fprompt_weight","Fprompt",100.0,0.0,1.1);
//TH1F *EndTime = new TH1F("Hit Time","Hit Time",100.0,0.0,1000);
//TH1F *EndTime_weight = new TH1F("Hit Time","Hit Time",100.0,0.0,1000);
//TH3D *neutronStartPos=new TH3D("neutronStartPos", "neutronStartPos", 20, -100.0, 1000.0, 30, -1000.0, 1000.0, 100.0, 0.0, 2600.0);
//TH3D *neutronEndPos=new TH3D("neutronEndPos", "neutronEndPos", 70.0, -3500.0, 3500.0, 90, -4500.0, 4500.0, 35.0, -1750.0, 1750.0);
TH3F *lre_fpr_pos = new TH3F("lrecoil_fprompt_vs_pos","lrecoil_fprompt_vs_pos", 100.0,0.0,1.0,100.0,0.0,1.1,200.0,0.0,600.0);
TH3F *lre_fpr_pos_weight = new TH3F("lrecoil_fprompt_vs_pos_weight","lrecoil_fprompt_vs_pos", 100.0,0.0,1.0,100.0,0.0,1.1,200.0,0.0,600.0);
TH3F *totPE_kE_pos =new TH3F("totPE_kE_pos","totPE_kE_pos",100.0,0.0,2000.0,50.0,0.0,0.25,200.0,0.0,600.0);
TH3F *totPE_kE_pos_weight =new TH3F("totPE_kE_pos_weight","totPE_kE_pos",100.0,0.0,2000.0,50.0,0.0,0.25,200.0,0.0,600.0);
TH3F *lre_fpr_kE = new TH3F("lrecoil_fprompt_vs_kE","lrecoil_fprompt_vs_kE",100.0,0.0,1.0,100.0,0.0,1.1,50.0,0.0,0.25);
TH3F *lre_fpr_kE_weight = new TH3F("lrecoil_fprompt_vs_kE_weight","lrecoil_fprompt_vs_kE",100.0,0.0,1.0,100.0,0.0,1.1,50.0,0.0,0.25);
TH3F *lre_fpr_pos_cut = new TH3F("lrecoil_fprompt_vs_pos_with_cuts","lrecoil_fprompt_vs_pos_with_cuts", 100.0,0.0,1.0,100.0,0.0,1.1,200.0,0.0,600.0);
TH3F *lre_fpr_pos_cut_weight = new TH3F("lrecoil_fprompt_vs_pos_with_cuts_weight","lrecoil_fprompt_vs_pos_with_cuts", 100.0,0.0,1.0,100.0,0.0,1.1,200.0,0.0,600.0);
TH3F *totPE_kE_pos_cut =new TH3F("totPE_kE_pos_with_cuts","totPE_kE_pos_with_cuts",100.0,0.0,2000.0,50.0,0.0,0.25,200.0,0.0,600.0);
TH3F *totPE_kE_pos_cut_weight =new TH3F("totPE_kE_pos_with_cuts_weight","totPE_kE_pos_with_cuts",100.0,0.0,2000.0,50.0,0.0,0.25,200.0,0.0,600.0);
TH3F *lre_fpr_kE_cut = new TH3F("lrecoil_fprompt_vs_kE_with_cuts","lrecoil_fprompt_vs_kE_with_cuts",100.0,0.0,1.0,100.0,0.0,1.1,50.0,0.0,0.25);
TH3F *lre_fpr_kE_cut_weight = new TH3F("lrecoil_fprompt_vs_kE_with_cuts_weight","lrecoil_fprompt_vs_kE_with_cuts",100.0,0.0,1.0,100.0,0.0,1.1,50.0,0.0,0.25);
/*
TH1F *x_mom = new TH1F("Initial x momentum","X_mom",50.0,-140.0,70.0);
TH1F *y_mom = new TH1F("Initial y momentum","Y_mom",50.0,-110.0,110.0);
TH1F *z_mom = new TH1F("Initial z momentum","Z_mom",50.0,-140.0,-10.0);
TH1F *x_mom_tot = new TH1F("Initial x momentum","X_mom_tot",50.0,-140.0,70.0);
TH1F *y_mom_tot = new TH1F("Initial y momentum","Y_mom_tot",50.0,-110.0,110.0);
TH1F *z_mom_tot = new TH1F("Initial z momentum","Z_mom_tot",50.0,-140.0,-10.0);
*/
TH1F *angle_hist = new TH1F("Triggered_Events","Initial angle",50.0,-0.25,1.0);
TH1F *angle_hist_weight = new TH1F("Triggered_Events_weight","Initial angle",50.0,-0.25,1.0);
TH1F *angle_hist_tot = new TH1F("All_Events","Initial angle",50.0,-0.25,1.0);
TH1F *angle_hist_tot_weight = new TH1F("All_Events_weight","Initial angle",50.0,-0.25,1.0);
//Get total number of events in chain//
int nEvents = T->GetEntries();
//Print out total number of entries into terminal//
cout << "number of events: " << nEvents << endl;
//Initializing variables to be used later to count how many events pass certain cuts//
int nlrecoilNuclear = 0;
int nlrecoilElec = 0;
int nradial = 0;
int ntotalPE = 0;
int nfprompt = 0;
int allcut = 0;
int rad_pe_cut = 0;
int rad_lr_cut = 0;
int rad_fp_cut = 0;
int pe_lr_cut = 0;
int pe_fp_cut = 0;
int lr_fp_cut = 0;
int pe_lr_fp_cut = 0;
int pe_lr_rad_cut = 0;
int pe_fp_rad_cut = 0;
int lr_fp_rad_cut = 0;
int N = 1;
int totalEvEvents = 0;
Double_t weight = 0.0;
//For loop interating over all events//
for(Long64_t iev=0; iev<nEvents; iev++){
//Get the ith event//
T->GetEntry(iev);
//Print out every 1000 events the number of event the for loop is on//
if(iev%1000==0) cout << "event " << iev << endl;
//Get the number of tracks for this event//
Long64_t ntracks=rds->GetMC()->GetMCTrackCount();
N++;
//set pathway to mc branch
RAT::DS::MC *rmc = rds->GetMC();
//Get the montecarlo (true) kinetic energy for this event//
Double_t kEnergy=rmc->GetMCParticle(0)->GetKE();
//then use that energy to get the wieght for this event//
weight=ambe_mag[(int)trunc(kEnergy/0.109091)];
//Fill two histograms, one with the weight one without//
kE->Fill(kEnergy);
kE_weight->Fill(kEnergy,weight);
//optional code to print out weight and energy//
//cout << weight << " " << kEnergy << endl;
//Get number of ev events (will be different than total number of events simulated above)//
int nEvEvents = rds->GetEVCount();
totalEvEvents += nEvEvents;
//iterating over ev events//
for(int ev_count=0; ev_count<nEvEvents; ev_count++){
//Get reconstructed (ev branch) variables of interest and fill various histograms (currently commented out but can be uncommented)//
Double_t kenergy = rds->GetEV(ev_count)->GetShellFit()->GetKE();
//kE->Fill(kenergy,weight);
Double_t FPrompt = rds->GetEV(ev_count)->GetFprompt();
//Counting number of events that pass this cut//
if(FPrompt>0.681){
nfprompt++;
//fprompt->Fill(FPrompt,weight);
}
//fprompt->Fill(FPrompt,weight);
Double_t LRecoil = rds->GetEV(ev_count)->GetLrecoil();
//Counting number of events that pass this cut//
if(LRecoil>0.373){
nlrecoilNuclear++;
//lrecoil->Fill(LRecoil,weight);
}
//lrecoil->Fill(LRecoil,weight);
Double_t totalpe = rds->GetEV(ev_count)->GetShellFit()->GetTotalPE();
//Counting number of events that pass this cut//
if(totalpe>20.0 && totalpe<500.0){
ntotalPE++;
//totalPE->Fill(totalpe,weight);
}
//totalPE->Fill(totalpe,weight);
//if(totalpe != 0){
//Double_t hit_time = rds->GetMC()->GetPMT(iev)->GetEndTime();
//EndTime->Fill(hit_time,weight);
//}
//calculating distance from event to center of detector
Double_t xpos = rds->GetEV(ev_count)->GetShellFit()->GetPosition().x();
Double_t ypos = rds->GetEV(ev_count)->GetShellFit()->GetPosition().y();
Double_t zpos = rds->GetEV(ev_count)->GetShellFit()->GetPosition().z();
Double_t pos = sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos));
/*
if(pos<400.0){
nradial++;
radialPos->Fill(sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)),weight);
}
if(kenergy<.25){
kE->Fill(kenergy);
//totalPE->Fill(totalpe,weight);
}
*/
//Filling 3D histograms with desired variables//
lre_fpr_pos->Fill(LRecoil,FPrompt,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)));
totPE_kE_pos->Fill(totalpe,kenergy,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)));
lre_fpr_kE->Fill(LRecoil,FPrompt,kenergy);
lre_fpr_pos_weight->Fill(LRecoil,FPrompt,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)),weight);
totPE_kE_pos_weight->Fill(totalpe,kenergy,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)),weight);
lre_fpr_kE_weight->Fill(LRecoil,FPrompt,kenergy,weight);
radialPos->Fill(sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)));
radialPos_weight->Fill(sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)),weight);
//Counting the rest of the cuts//
if(totalpe>20.0 && totalpe<500.0 && LRecoil>0.373 && FPrompt>0.681 && pos<400.0 && kenergy<.25){
allcut++;
}
if(totalpe>20.0 && totalpe<500.0 && LRecoil>0.373 && FPrompt>0.681 && kenergy<.25){
pe_lr_fp_cut++;
}
if(totalpe>20.0 && totalpe<500.0 && pos<400.0 && kenergy<.25){
totPE_kE_pos_cut->Fill(totalpe,kenergy,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)));
totPE_kE_pos_cut_weight->Fill(totalpe,kenergy,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)),weight);
rad_pe_cut++;
}
if(pos<400.0 && LRecoil>0.373 && kenergy<.25){
rad_lr_cut++;
}
if(pos<400.0 && FPrompt>0.681 && kenergy<.25){
rad_fp_cut++;
}
if(totalpe>20.0 && totalpe<500.0 && LRecoil>0.373 && kenergy<.25){
pe_lr_cut++;
}
if(totalpe>20.0 && totalpe<500.0 && FPrompt>0.681 && kenergy<.25){
pe_fp_cut++;
}
if(LRecoil>0.373 && FPrompt>0.681 && kenergy<.25){
lre_fpr_kE_cut->Fill(LRecoil,FPrompt,kenergy);
lre_fpr_kE_cut_weight->Fill(LRecoil,FPrompt,kenergy,weight);
lr_fp_cut++;
}
if(totalpe>20.0 && totalpe<500.0 && LRecoil>0.373 && FPrompt>0.681 && kenergy<.25){
//pe_lr_fp_cut++;
}
if(totalpe>20.0 && totalpe<500.0 && LRecoil>0.373 && pos<400.0 && kenergy<.25){
pe_lr_rad_cut++;
}
if(totalpe>20.0 && totalpe<500.0 && FPrompt>0.681 && pos<400.0 && kenergy<.25){
pe_fp_rad_cut++;
}
if(LRecoil>0.373 && FPrompt>0.681 && pos<400.0 && kenergy<.25){
lre_fpr_pos_cut->Fill(LRecoil,FPrompt,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)));
lre_fpr_pos_cut_weight->Fill(LRecoil,FPrompt,sqrt((xpos*xpos)+(ypos*ypos)+(zpos*zpos)),weight);
lr_fp_rad_cut++;
}
}
//The following section of code gets the volume the neutron gets captured in to understand where the non-triggering events end up. It does this by iterating through the tracks of each event.//
for(Int_t j=0; j<ntracks;j++){
//Specify path to the first track//
MCTrack* firstTrack=rds->GetMC()->GetMCTrack(j);
//Get particle name to only look at neutrons//
std::string trackname=firstTrack->GetParticleName();
if(trackname.compare("neutron")==0){
//Get number of steps in each track//
Int_t numSteps=firstTrack->GetMCTrackStepCount();
//std::string lastVol;
//for(Int_t iStep=0; iStep<numSteps; iStep++){
//Get initial position and volume name//
MCTrackStep* step = firstTrack->GetMCTrackStep(0);
std::string proc = step->GetProcess();
TVector3 curPos = step->GetEndpoint();
std::string firstVol = step->GetVolume();
if(j==0){
//cout << "First volume: " << firstVol << endl;
}
//If it is made by the initial process "start" then fill histogram with the start position//
if(proc.compare("start")==0){
//neutronStartPos->Fill(curPos.X(), curPos.Y(), curPos.Z(),weight);
}
//Now get the last step, or when the particle is captured, in this track//
MCTrackStep* step2 = firstTrack->GetMCTrackStep(numSteps-1);
std::string proc2 = step2->GetProcess();
TVector3 curPos2 = step2->GetEndpoint();
std::string lastVol2 = step2->GetVolume();
MCTrackStep* prevStep= firstTrack->GetMCTrackStep(numSteps-2);
//Since we only care about neutrons, only look at when the neutron is captured//
if(proc2.compare("nCapture")==0){
//cout << "First volume name: " << firstVol << " track: " << j << " step: " << numSteps-1 << " energy: " << prevStep->GetKE() << endl;
//neutronEndPos->Fill(curPos2.X(), curPos2.Y(), curPos2.Z(),weight);
//}
}
}
}
}
//Iterating over events in the mc branch(same as above)//
for(int i=0; i<nEvents;i++){
T->GetEntry(i);
RAT::DS::MC *rmc = rds->GetMC();
Double_t numberPE=rmc->GetNumPE();
//Get various variables//
Double_t xPos=rmc->GetMCParticle(0)->GetPosition().x();
Double_t yPos=rmc->GetMCParticle(0)->GetPosition().y();
Double_t zPos=rmc->GetMCParticle(0)->GetPosition().z();
Double_t kEnergy=rmc->GetMCParticle(0)->GetKE();
Double_t fpmpt=rmc->GetMCFprompt();
Double_t xmom = rmc->GetMCParticle(0)->GetMomentum().x();
Double_t ymom = rmc->GetMCParticle(0)->GetMomentum().y();
Double_t zmom = rmc->GetMCParticle(0)->GetMomentum().z();
Double_t lenSq1 = (xmom*xmom + ymom*ymom + zmom*zmom);
Double_t lenSq2 = (.580*.580 + 1.4*1.4);
Double_t dot = (-xmom*.580 - zmom*1.4);
Double_t angle = acos(dot/sqrt(lenSq1 * lenSq2));
//x_mom_tot->Fill(xmom,weight);
//y_mom_tot->Fill(ymom,weight);
//z_mom_tot->Fill(zmom,weight);
//Filling histograms with variables(with weights in some cases)//
angle_hist_tot->Fill(angle);
angle_hist_tot_weight->Fill(angle,weight);
//char particleName[] = rmc->GetMCParticle(0)->GetParticleName().c_str();
//cout<<particleName;
//This if statement makes it so only events that triggered a PMT response will be looked at//
if (numberPE != 0){
//x_mom->Fill(xmom,weight);
//y_mom->Fill(ymom,weight);
//z_mom->Fill(zmom,weight);
angle_hist->Fill(angle);
angle_hist_weight->Fill(angle,weight);
//radialPos->Fill(sqrt((xPos*xPos)+(yPos*yPos)+(zPos*zPos)));
//radialPos_weight->Fill(sqrt((xPos*xPos)+(yPos*yPos)+(zPos*zPos)),weight);
//numPE0->Fill(numberPE,weight);
//numPE1->Fill(numberPE,weight);
//kE->Fill(kEnergy,weight);
fprompt->Fill(fpmpt);
fprompt_weight->Fill(fpmpt,weight);
}
}
//This massive commented out section will display your histograms immediately and won't save them in a root file to be further manipulated//
cout << "test" << endl;
//TCanvas* c1 = new TCanvas("x_mom","Monte Carlo 1",600,600);
//x_mom->GetXaxis()->SetTitle("initial x momentum");
//x_mom->Draw();
//lrecoil->Write();
//gPad->SetLogy();
cout << "test2" << endl;
//TCanvas* c2 = new TCanvas("x_mom_tot","Monte Carlo 2",600,600);
//x_mom_tot->GetXaxis()->SetTitle("xmom");
//gPad->SetLogy();
//x_mom_tot->Draw();
//fprompt->Write();
//TCanvas* c3 = new TCanvas("c_3","Monte Carlo 3",600,600);
//radialPos->GetXaxis()->SetTitle("radial position (mm)");
//gPad->SetLogy();
//radialPos->Draw();
//radialPos->Write();
//TCanvas* c4 = new TCanvas("c_4","Monte Carlo 4",600,600);
//kE->GetXaxis()->SetTitle("kinetic energy (MeV)");
//gPad->SetLogy();
//kE->Draw();
//kE->Write();
//TCanvas* c5 = new TCanvas("c_5","Monte Carlo 5",600,600);
//totalPE->GetXaxis()->SetTitle("total number of photoelectrons");
//gPad->SetLogy();
//gStyle->SetOptStat(111111);
//totalPE->Draw();
//totalPE->Write();
//Double_t overflow = totalPE->GetBinContent(101);
//cout << "Total PE overflow: " << overflow << endl;
//TCanvas* c6 = new TCanvas("c_6","Monte Carlo 6",600,600);
//neutronEndPos->Draw();
//gPad->SetLogy();
//TCanvas* c7 = new TCanvas("c_7","Monte Carlo 7",600,600);
//neutronStartPos->Draw();
//Printing out the number of events that pass each combination of cuts//
cout << "Number of lrecoil: " << nlrecoilNuclear << endl << "Number of radial: " << nradial << endl << "Number of totalPE: " << ntotalPE << endl << "Number of fprompt: " << nfprompt << endl << "All cuts: " << allcut << endl << "pe_lr_fp: " << pe_lr_fp_cut << endl << "rad_pe: " << rad_pe_cut << endl << "rad_lr: " << rad_lr_cut << endl << "rad_fp: " << rad_fp_cut << endl << "pe_lr: " << pe_lr_cut << endl << "pe_fp: " << pe_fp_cut << endl << "lr_fp: " << lr_fp_cut << endl << "pe_lr_fp: " << pe_lr_fp_cut << endl << "pe_lr_rad: " << pe_lr_rad_cut << endl << "pe_fp_rad: " << pe_fp_rad_cut << endl << "lr_fp_rad: " << lr_fp_rad_cut << endl << "Total number of ev events: " << totalEvEvents << endl;
//Writing the histograms to be saved into the output root file//
totPE_kE_pos->Write();
lre_fpr_pos->Write();
lre_fpr_kE->Write();
totPE_kE_pos_cut->Write();
lre_fpr_pos_cut->Write();
lre_fpr_kE_cut->Write();
totPE_kE_pos_weight->Write();
lre_fpr_pos_weight->Write();
lre_fpr_kE_weight->Write();
totPE_kE_pos_cut_weight->Write();
lre_fpr_pos_cut_weight->Write();
lre_fpr_kE_cut_weight->Write();
totalPE->Write();
totalPE_weight->Write();
//lrecoil->Write();
//lrecoil_weight->Write();
fprompt->Write();
fprompt_weight->Write();
radialPos->Write();
radialPos_weight->Write();
//x_mom->Write();
//y_mom->Write();
//z_mom->Write();
//x_mom_tot->Write();
//y_mom_tot->Write();
//z_mom_tot->Write();
angle_hist_tot->Write();
angle_hist->Write();
kE->Write();
kE_weight->Write();
angle_hist_tot_weight->Write();
angle_hist_weight->Write();
//Close output file//
outFile->Close();
//Delete the chain to be complete and not have memory issues//
delete T;
}