Adding sensitivity scripts

background_spectra/PlotBackground.cxx

@@ -0,0 +1,156 @@
+void Plot(TString BG);
+void PlotBackground();
+
+void PlotBackground() {
+
+  vector<TString> theBG = {"_all", "_P", "_TC", "_UC", "_C"};
+  for (int i=0; i<theBG.size(); i++) {
+    Plot( theBG[i] );
+  }
+
+}
+void Plot(TString BG){
+
+
+  TFile *f = new TFile(TString::Format("HadronicBackground.spectrum%s.root", BG.Data()).Data());
+  TCanvas * C = (TCanvas*)f->Get("c1");
+  TH1D * h_had = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+  f = new TFile(TString::Format("LeptonicBackground.spectrum%s.root",  BG.Data()).Data());
+  C = (TCanvas*)f->Get("c1");
+  TH1D * h_lep = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+  f = new TFile(TString::Format("PhotonicBackground.spectrum%s.root",  BG.Data()).Data());
+  C = (TCanvas*)f->Get("c1");
+  TH1D * h_pho = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+  f = new TFile(TString::Format("TrappedHadronicBackground.spectrum%s.root",  BG.Data()).Data());
+  C = (TCanvas*)f->Get("c1");
+  TH1D * h_thad = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+  f = new TFile(TString::Format("TrappedLeptonicBackground.spectrum%s.root",  BG.Data()).Data());
+  C = (TCanvas*)f->Get("c1");
+  TH1D * h_tlep = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+  f = new TFile(TString::Format("HadronicBackgroundDecay.spectrum%s.root",  BG.Data()).Data());
+  C = (TCanvas*)f->Get("c1");
+  TH1D * h_had_decay = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+  f = new TFile(TString::Format("TrappedHadronicBackgroundDecay.spectrum%s.root",  BG.Data()).Data());
+  C = (TCanvas*)f->Get("c1");
+  TH1D * h_thad_decay = (TH1D*)C->GetPrimitive("EnergySpectrum");
+
+
+
+	TCanvas * C_tot = new TCanvas("C_tot","C_tot", 1181, 802);
+	C_tot->SetLeftMargin(0.15);
+	C_tot->SetBottomMargin(0.15);
+
+	h_had->SetFillStyle(0);
+	h_had->SetLineColor(kRed);
+	h_had->SetLineWidth(2);
+	h_had->SetTitle("");
+	h_had->GetYaxis()->SetTitle("Flux (#gamma/keV/s)");
+	h_had->GetXaxis()->SetRangeUser(10, 1e6);
+	h_had->GetYaxis()->SetTitleSize(0.06);
+	h_had->GetYaxis()->SetTitleOffset(1.0);
+	h_had->GetYaxis()->SetLabelSize(0.05);
+	h_had->GetXaxis()->SetTitleSize(0.06);
+	h_had->GetXaxis()->SetLabelSize(0.05);
+	h_had->GetXaxis()->SetTitleOffset(1.1);
+
+	//Convert the x-axis to MeV instead of keV
+	//for (int i = 0; i < h_had->GetNBins(); i++) {
+	//	h_had->Get
+
+
+	double T = 3600; //time in simulated seconds
+	double scale = 1.0/3600.; //conversion to event rate (per second)
+
+	//Simulations were for 3600 seconds, so divide by time
+	h_had->Scale(scale);
+	//Multiply by 1000 to conver to /MeV
+	//h_had->Scale(1000.);
+	h_had->GetYaxis()->SetRangeUser(2e-6, 2);
+	h_had->Draw("hist");
+
+	h_lep->SetFillStyle(0);
+	h_lep->SetLineColor(kGreen+1);
+	h_lep->SetLineWidth(2);
+	h_lep->Scale(scale);
+	//h_lep->Scale(1000.);
+	h_lep->Draw("hist same");
+
+	h_pho->SetFillStyle(0);
+	h_pho->SetLineColor(kMagenta);
+	h_pho->SetLineWidth(2);
+	h_pho->Scale(scale);
+	//h_pho->Scale(1000.);
+	h_pho->Draw("hist same");
+
+	h_thad->SetFillStyle(0);
+	h_thad->SetLineColor(kBlue);
+	h_thad->SetLineWidth(2);
+	h_thad->Scale(scale);
+	//h_thad->Scale(1000.);
+	h_thad->Draw("hist same");
+
+	h_tlep->SetFillStyle(0);
+	h_tlep->SetLineColor(kOrange);
+	h_tlep->SetLineWidth(2);
+	h_tlep->Scale(scale);
+	//h_tlep->Scale(1000.);
+	h_tlep->Draw("hist same");
+
+        h_had_decay->SetFillStyle(0);
+        h_had_decay->SetLineColor(kRed+1);
+	h_had_decay->SetLineStyle(2);
+        h_had_decay->SetLineWidth(2);
+        h_had_decay->Scale(scale);
+        //h_pho->Scale(1000.);
+        h_had_decay->Draw("hist same");
+
+	h_thad_decay->SetFillStyle(0);
+	h_thad_decay->SetLineColor(kBlue+1);
+	h_thad_decay->SetLineStyle(2);
+	h_thad_decay->SetLineWidth(2);
+	h_thad_decay->Scale(scale);
+	//h_thad_decay->Scale(1000.);
+	h_thad_decay->Draw("hist same");
+
+
+	TH1D * h_tot = (TH1D*)h_had->Clone("h_tot");
+	h_tot->Add(h_lep);
+	h_tot->Add(h_pho);
+	h_tot->Add(h_thad);
+	h_tot->Add(h_tlep);
+	h_tot->Add(h_had_decay);
+	h_tot->Add(h_thad_decay);
+	h_tot->SetLineColor(kBlack);
+	h_tot->SetLineWidth(2);
+	h_tot->Draw("hist same");
+
+	C_tot->SetLogx();
+	C_tot->SetLogy();
+	gPad->RedrawAxis();
+
+	TLegend * leg = new TLegend(0.65, 0.59, 0.95, 0.95);
+	leg->AddEntry(h_had, "Hadrons", "l");
+	leg->AddEntry(h_had_decay, "Hadron Activation", "l");
+	leg->AddEntry(h_thad, "Trapped Hadrons", "l");
+	leg->AddEntry(h_thad_decay, "Trapped Hadron Activation", "l");
+	leg->AddEntry(h_lep, "Leptons", "l");
+	leg->AddEntry(h_tlep, "Trapped Leptons", "l");
+	leg->AddEntry(h_pho, "Photons", "l");
+	leg->AddEntry(h_tot, "Total", "l");
+	leg->Draw();
+
+	C_tot->SaveAs(TString::Format("bg_pixel%s.png", BG.Data()).Data() );
+	C_tot->SaveAs(TString::Format("bg_pixel%s.pdf", BG.Data()).Data() );
+	C_tot->SaveAs(TString::Format("bg_pixel%s.root", BG.Data()).Data() );
+
+	h_tot->SetName(TString::Format("totalBG%s", BG.Data() ).Data() );
+	h_tot->SaveAs(TString::Format("totalBG_pixel%s.root", BG.Data()).Data() );
+
+};
+

background_spectra/PlotTotalBG.C

@@ -0,0 +1,101 @@
+{
+    TCanvas * c = new TCanvas("c", "c", 1000, 600);
+
+    TFile * theFile = new TFile( "pixel_rightOrbit/totalBG_pixel.root");
+    //TFile * theFile = new TFile( "base_wrongOrbit/totalBG_DSSD.root");
+    theFile->cd();
+
+    TH1D* totalBG_all = (TH1D*)theFile->Get("totalBG_all");
+    TH1D* totalBG_P = (TH1D*)theFile->Get("totalBG_P");
+    TH1D* totalBG_UC = (TH1D*)theFile->Get("totalBG_UC");
+    TH1D* totalBG_TC = (TH1D*)theFile->Get("totalBG_TC");
+
+    totalBG_all->Rebin(10);
+    totalBG_all->Scale(0.1);
+    totalBG_P->Rebin(10);
+    totalBG_P->Scale(0.1);
+    totalBG_TC->Rebin(10);
+    totalBG_TC->Scale(0.1);
+    totalBG_UC->Rebin(10);
+    totalBG_UC->Scale(0.1);
+
+    totalBG_all->GetXaxis()->SetRangeUser(50, 1e6);
+
+    totalBG_all->SetLineColor(kGray);
+    totalBG_all->SetLineWidth(3);
+    totalBG_all->DrawCopy("hist l");
+    totalBG_P->SetLineColor(kRed);
+    totalBG_P->SetLineWidth(3);
+    totalBG_P->DrawCopy("hist l same ");
+    totalBG_TC->SetLineWidth(3);
+    totalBG_UC->SetLineWidth(3);
+    totalBG_UC->SetLineColor(kBlue);
+    totalBG_TC->SetLineColor(kGreen);
+    totalBG_TC->DrawCopy("hist l same ");
+    totalBG_UC->DrawCopy("hist l same ");
+    c->SetLogy();
+    c->SetLogx();
+
+    c->SetLeftMargin(0.15);
+    c->SetRightMargin(0.05);
+    c->SetBottomMargin(0.15);
+    c->SetTopMargin(0.05);
+
+    TLegend * legend = new TLegend(0.65,0.7,0.95,0.95);
+    //legend->SetHeader("The Legend Title","C"); // option "C" allows to center the header
+    legend->AddEntry(totalBG_all ,"Total BG","l");
+    legend->AddEntry(totalBG_P ,"Fake pair events","l");
+    legend->AddEntry(totalBG_TC ,"Fake tracked Compton","l");
+    legend->AddEntry(totalBG_UC ,"Fake untracked Compton","l");
+    legend->Draw();
+
+
+    //TFile * theFile = new TFile( "pixel_rightOrbit/totalBG_pixel.root");
+    TFile * theFile2 = new TFile( "pixel_EHCut/totalBG_pixel_EHCut.root");
+    theFile2->cd();
+
+    theFile2->ls();
+
+    TH1D* totalBG_all_CUT = (TH1D*)theFile2->Get("totalBG_all");
+    TH1D* totalBG_P_CUT = (TH1D*)theFile2->Get("totalBG_P");
+    TH1D* totalBG_UC_CUT = (TH1D*)theFile2->Get("totalBG_UC");
+    TH1D* totalBG_TC_CUT = (TH1D*)theFile2->Get("totalBG_TC");
+
+    totalBG_all_CUT->Rebin(10);
+    totalBG_all_CUT->Scale(0.1);
+    totalBG_P_CUT->Rebin(10);
+    totalBG_P_CUT->Scale(0.1);
+    totalBG_TC_CUT->Rebin(10);
+    totalBG_TC_CUT->Scale(0.1);
+    totalBG_UC_CUT->Rebin(10);
+    totalBG_UC_CUT->Scale(0.1);
+
+
+    totalBG_all_CUT->SetLineStyle(7);
+    totalBG_P_CUT->SetLineStyle(7);
+    totalBG_UC_CUT->SetLineStyle(7);
+    totalBG_TC_CUT->SetLineStyle(7);
+
+    totalBG_all_CUT->SetLineColor(kGray);
+    totalBG_all_CUT->SetLineWidth(3);
+    totalBG_all_CUT->DrawCopy("hist l same");
+    totalBG_P_CUT->SetLineColor(kRed);
+    totalBG_P_CUT->SetLineWidth(3);
+    totalBG_P_CUT->DrawCopy("hist l same ");
+    totalBG_TC_CUT->SetLineWidth(3);
+    totalBG_UC_CUT->SetLineWidth(3);
+    totalBG_UC_CUT->SetLineColor(kBlue);
+    totalBG_TC_CUT->SetLineColor(kGreen);
+    totalBG_TC_CUT->DrawCopy("hist l same ");
+    totalBG_UC_CUT->DrawCopy("hist l same ");
+
+
+    legend->AddEntry(totalBG_all_CUT ,"after horizon cut","l");
+//    legend->AddEntry(totalBG_P_CUT ,"Fake pair events","l");
+//    legend->AddEntry(totalBG_TC_CUT ,"Fake tracked Compton","l");
+//    legend->AddEntry(totalBG_UC_CUT ,"Fake untracked Compton","l");
+//    legend->Draw();
+
+    c->SaveAs("AMEGO-X-Pixel-BG_EHCut.png");
+
+}

background_spectra/README.md

@@ -0,0 +1,20 @@
+Various scripts and config files to obtain simulated background spectra from `.tra` files.
+
+`RemoveBottomTrackerLeptonicInteractions.py`: This should be called once after the background sims are run. Removes background component of trapped leptonic events that come in via the gap below the ACD.
+
+`mimrec_bg_spectra.sh`: Runs mimrec over tra files to produce root files with background spectra. Input/output directory, config files and detector geometry are all hardcoded and should be adjusted. This produces one root file per component (TrappedLeptonic, Photonic etc) and event type (Tracked Compton, Pair etc). The spectra are not normalized by exposure time yet!
+
+`PlotBackground.cxx`: Reads in the root files from the previous step to normalize and plot the spectra and calculate the total background rates (sum over all components). Again, input and output files and simulated time are hardcoded. Output: plots (png, pdf) of the reconstructed background spectrum for each component, root files with said spectra, and root files with just the total background histogram for each event type. Use the `hadd` macro to combine those into one file if desired. 
+
+`PlotTotalBG.C`: Root script to plot total Untracked Compton/Tracked Compton/Pair background spectra (summed over all components), using the output from the previous script.
+
+`mimrec_AMEGOX_PlotBkgSpectrum_all.cfg`: Config file, selecting both Compton and pair events.
+
+`mimrec_AMEGOX_PlotBkgSpectrum_C.cfg`: Config file, Compton only.
+
+`mimrec_AMEGOX_PlotBkgSpectrum_P.cfg`: Config file, Pair only.
+
+`mimrec_AMEGOX_PlotBkgSpectrum_TC.cfg`: Config file, tracked Compton only.
+
+`mimrec_AMEGOX_PlotBkgSpectrum_UC.cfg`: Config file, untracked Compton only.
+

background_spectra/RemoveBottomTrackerLeptonicInteractions.py

@@ -0,0 +1,113 @@
+import gzip
+import re
+import os
+import glob
+
+#trafile = gzip.open('TraFiles_100420/R1_Files/TrappedLeptonicBackground/TrappedLeptonicBackground.p2.inc30.id1.tra.gz', 'rb')
+#newtrafile = gzip.open('TrappedLeptonicBackground_mod.p2.inc30.id1.tra.gz', 'wb')
+
+def RemoveEvents(filename, newtrafilename):
+
+	trafile = gzip.open(filename, 'rb')
+	print('Reading in file ', filename) 
+
+	newtrafile = gzip.open(newtrafilename, 'wb')
+
+
+	#Write the first 7 lines of the tra file into the new condensed file
+	for x in range(7):
+		line = trafile.readline()
+		newtrafile.write(line)
+
+	#Read the whole file, and append the event details to new_event. If the position is not consistent with the last tracker layer, save the event into the new trafile
+	while True:
+		line = trafile.readline()
+		if not line:
+			break
+		#Find the start of an event with the SE flag
+		if re.match('SE', line.decode('utf-8')):
+			new_event = []
+			#append SE to be the first element in the new_event list
+			new_event.append(line)
+
+			#Read the next line and see if it matches with a CO event (since these are the only ones in the TrappedLeptonic file other than unknowns)
+			line = trafile.readline()
+			if re.match('ET CO', line.decode('utf-8')):
+				new_event.append(line)
+
+
+	#			line = trafile.readline()
+	#			if re.match('ID 5349', line.decode('utf-8')):
+	#				print('found ID 5349')
+	#				new_event.append(line)
+
+				#If it's a compton event, append the next 9 lines of the event details to the new_event list to get to the position info
+				for x in range(7):
+					line = trafile.readline()
+					new_event.append(line)
+				#Some events have the TQ flag which makes them one line longer, so we have to check for that
+				if re.match('TQ', line.decode('utf-8')):
+					for x in range(3):
+						line = trafile.readline()
+						new_event.append(line)
+				else: 
+					for x in range(2):
+						line = trafile.readline()
+						new_event.append(line)
+
+
+
+				#Read the CH line and see if it's consistent with the first event being in the bottom layer of the tracker
+				line = trafile.readline()
+				if re.match('CH 0 -?\d+\.\d+ -?\d+\.\d+ 0.75', line.decode('utf-')):
+					#print('Got a bad one!')
+					pass
+
+				#If the event isn't consistent, then write it to the new tra file
+				else:	
+					new_event.append(line)
+					line = trafile.readline()
+					while re.match('CH \d ', line.decode('utf-8')):
+						new_event.append(line)
+						line = trafile.readline()
+					#print(new_event)
+					for i in new_event:
+						newtrafile.write(i)
+
+		#If you've reached the last event, copy over the final lines of the file to include the FT info
+		if re.match('EN', line.decode('utf-8')):
+			newtrafile.write(line)
+			while True:
+				line = trafile.readline()
+				newtrafile.write(line)
+				if not line:
+					break
+
+
+
+	trafile.close()
+	newtrafile.close()
+
+	return
+
+###########################################333
+
+
+
+
+#Change the path here!
+#path = 'TraFiles_040520/R1_Files/TrappedLeptonicBackground/'
+path = '/data/slag2/hfleisc1/amego_sims/simfiles/pixel/background/TrappedLeptonicBackground/'
+path = '/data/slag2/hfleisc1/amego_sims/simfiles/base/background/TrappedLeptonicBackground/'
+path = '/data/slag2/hfleisc1/amego_sims/simfiles/pixel/BG_above_0.01MeV/TrappedLeptonicBackground/'
+#for filename in glob.glob(os.path.join(path, 'TrappedLeptonicBackground.p*.inc*.tra.gz')):
+for filename in os.listdir(path):
+        if os.path.isfile(os.path.join(path, filename)):
+                if re.match('TrappedLeptonicBackground.p\d.inc\d+.id1.tra.gz', filename):
+                        newtrafilename = filename.replace('TrappedLeptonicBackground','TrappedLeptonicBackground_old')
+                        print (filename, newtrafilename)
+                        os.rename( os.path.join(path, filename), os.path.join(path, newtrafilename))
+                        RemoveEvents(os.path.join(path, newtrafilename), os.path.join(path,filename))
+
+
+