#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <time.h>
#include <cstdio>
#include <cstdlib>
#include <unistd.h>
#include <cstring>
#include <unistd.h>
#include "TTreeIndex.h"
#include "TChain.h"
#include "TH1.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TRandom.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TTree.h"
#include "TLegend.h"
#include "TLine.h"
#include "TROOT.h"
#include "TPostScript.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TText.h"
#include "TProfile.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TMath.h"
#include "TVector3.h"
#include "TRotation.h"
#include "TSpline.h"
#include "Math/InterpolationTypes.h"
#include "Math/Interpolator.h"
#include "Math/Integrator.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TLatex.h"
#include "Constants.h"
#include "Settings.h"
#include "Position.h"
#include "EarthModel.h"
#include "Tools.h"
#include "Vector.h"
#include "IceModel.h"
#include "Trigger.h"
#include "Spectra.h"
#include "signal.hh"
#include "secondaries.hh"
#include "Ray.h"
#include "counting.hh"
#include "Primaries.h"
#include "Efficiencies.h"
#include "Event.h"
#include "Detector.h"
#include "Report.h"
using namespace std;
/////////////////Plotting Script for nuflavorint for AraSIMQC
/////////////////Created by Kaeli Hughes, modified by Brian Clark
/////////////////Prepared on April 28 2016 as an "introductor script" to plotting with AraSim; for new users of queenbee
class EarthModel;
int main(int argc, char **argv)
{
gStyle->SetOptStat(111111); //this is a selection of statistics settings; you should do some googling and figure out exactly what this particular comibation does
gStyle->SetOptDate(1); //this tells root to put a date and timestamp on whatever plot we output
if(argc<2){ //if you don't have at least 1 file to run over, then you haven't given it a file to analyze; this checks this
cout<<"Not enough arguments! Abort run."<<endl;
}
//Create the histogram
TCanvas *c2 = new TCanvas("c2", "nuflavorint", 1100,850); //make a canvas on which to plot the data
TH1F *nuflavorint_hist = new TH1F("nuflavorint_hist", "nuflavorint histogram", 3, 0.5, 3.5); //create a histogram with three bins
nuflavorint_hist->GetXaxis()->SetNdivisions(3); //set the number of divisions
int total_thrown=0; // a variable to hold the grant total number of events thrown for all input files
for(int i=1; i<argc;i++){ //loop over the input files
string readfile; //create a variable called readfile that will hold the title of the simulation file
readfile = string(argv[i]); //set the readfile variable equal to the filename
Event *event = 0; //create a Event class pointer called event; note that it is set equal to zero to avoid creating a bald pointer
Report *report=0; //create a Event class pointer called event; note that it is set equal to zero to avoid creating a bald pointer
Settings *settings = 0; //create a Event class pointer called event; note that it is set equal to zero to avoid creating a bald pointer
TFile *AraFile=new TFile(( readfile ).c_str()); //make a new file called "AraFile" that will be simulation file we are reading in
if(!(AraFile->IsOpen())) return 0; //checking to see if the file we're trying to read in opened correctly; if not, bounce out of the program
int num_pass;//number of passing events
TTree *AraTree=(TTree*)AraFile->Get("AraTree"); //get the AraTree
TTree *AraTree2=(TTree*)AraFile->Get("AraTree2"); //get the AraTree2
AraTree2->SetBranchAddress("event",&event); //get the event branch
AraTree2->SetBranchAddress("report",&report); //get the report branch
AraTree->SetBranchAddress("settings",&settings); //get the settings branch
num_pass=AraTree2->GetEntries(); //get the number of passed events in the data file
AraTree->GetEvent(0); //get the first event; sometimes the tree does not instantiate properly if you'd explicitly "active" the first event
total_thrown+=(settings->NNU); //add the number of events from this file (The NNU variable) to the grand total; NNU is the number of THROWN neutrinos
for (int k=0; k<num_pass; k++){ //going to fill the histograms for as many events as were in this input file
AraTree2->GetEvent(k); //get the even from the tree
int nuflavorint; //make the container variable
double weight; //the weight of the event
int trigger; //the global trigger value for the event
nuflavorint=event->nuflavorint; //draw the event out; one of the objects in the event class is the nuflavorint, and this is the syntax for accessing it
weight=event->Nu_Interaction[0].weight; //draw out the weight of the event
trigger=report->stations[0].Global_Pass; //find out if the event was a triggered event or not
/*
if(trigger!=0){ //check if the event triggered
nuflavorint_hist->Fill(nuflavorint,weight); //fill the event into the histogram; the first argument of the fill (which is mandatory) is the value you're putting into the histogram; the second value is optional, and is the weight of the event in the histogram
//in this particular version of the code then, we are only plotting the TRIGGERED events; if you wanted to plot all of the events, you could instead remove this "if" condition and just Fill everything
}
*/
nuflavorint_hist->Fill(nuflavorint,weight); //fill the event into the histogram; the first argument of the fill (which is mandatory)
}
}
//After looping over all the files, make the plots and save them
//do some stuff to get the "total events thrown" text box ready; this is useful because on the plot itself you can then see how many events you THREW into the simulation; this is especially useful if you're only plotting passed events, but want to know what fraction of your total thrown that is
char *buffer= (char*) malloc (250); //declare a buffer pointer, and allocate it some chunk of memory
int a = snprintf(buffer, 250,"Total Events Thrown: %d",total_thrown); //print the words "Total Events Thrown: %d" to the variable "buffer" and tell the system how long that phrase is; the %d sign tells C++ to replace that "%d" with the next argument, or in this case, the number "total_thrown"
if(a>=250){ //if the phrase is longer than the pre-allocated space, increase the size of the buffer until it's big enough
buffer=(char*) realloc(buffer, a+1);
snprintf(buffer, 250,"Total Events Thrown: %d",total_thrown);
}
TLatex *u = new TLatex(.3,.01,buffer); //create a latex tex object which we can draw on the cavas
u->SetNDC(kTRUE); //changes the coordinate system for the tex object plotting
u->SetIndiceSize(.1); //set the size of the latex index
u->SetTextSize(.025); //set the size of the latex text
nuflavorint_hist->Draw(); //draw the histogram
nuflavorint_hist->GetXaxis()->SetTitle("Neutrino Flavor"); //set the x-axis label
nuflavorint_hist->GetYaxis()->SetTitle("Number of Events (weighted)"); //set the y-axis label
nuflavorint_hist->GetYaxis()->SetTitleOffset(1.5); //set the separation between the y-axis and it's label; root natively makes this smaller than is ideal
nuflavorint_hist->SetLineColor(kBlack); //set the color of the histogram to black, instead of root's default navy blue
u->Draw(); //draw the statistics box information
c2->SaveAs("outputs/plotting_example.png"); //save the canvas as a JPG file for viewing
c2->SaveAs("outputs/plotting_example.pdf"); //save the canvas as a PDF file for viewing
c2->SaveAs("outputs/plotting_example.root"); //save the canvas as a ROOT file for viewing or editing later
} //end main; this is the end of the script
|
NFOUR=1024
EXPONENT=21
NNU=300 // number of neutrino events
NNU_PASSED=10 // number of neutrino events that are allowed to pass the trigger
ONLY_PASSED_EVENTS=0 // 0 (default): AraSim throws NNU events whether or not they pass; 1: AraSim throws events until the number of events that pass the trigger is equal to NNU_PASSED (WARNING: may cause long run times if reasonable values are not chosen)
NOISE_WAVEFORM_GENERATE_MODE=0 // generate new noise waveforms for each events
NOISE_EVENTS=16 // number of pure noise waveforms
TRIG_ANALYSIS_MODE=0 // 0 = signal + noise, 1 = signal only, 2 = noise only
DETECTOR=1 // ARA stations 1 to 7
NOFZ=1
core_x=10000
core_y=10000
TIMESTEP=5.E-10 // value for 2GHz actual station value
TRIG_WINDOW=1.E-7 // 100ns which is actual testbed trig window
POWERTHRESHOLD=-6.06 // 100Hz global trig rate for 3 out of 16 ARA stations
POSNU_RADIUS=3000
V_MIMIC_MODE=0 // 0 : global trig is located center of readout windows
DATA_SAVE_MODE=0 // 2 : don't save any waveform informations at all
DATA_LIKE_OUTPUT=0 // 0 : don't save any waveform information to eventTree
BORE_HOLE_ANTENNA_LAYOUT=0
SECONDARIES=0
TRIG_ONLY_BH_ON=0
CALPULSER_ON=0
USE_MANUAL_GAINOFFSET=0
USE_TESTBED_RFCM_ON=0
NOISE_TEMP_MODE=0
TRIG_THRES_MODE=0
READGEOM=0 // reads geometry information from the sqlite file or not (0 : don't read)
TRIG_MODE=0 // use vpol, hpol separated trigger mode. by default N_TRIG_V=3, N_TRIG_H=3. You can change this values
number_of_stations=1
core_x=10000
core_y=10000
DETECTOR=1
DETECTOR_STATION=2
DATA_LIKE_OUTPUT=0
NOISE_WAVEFORM_GENERATE_MODE=0 // generates new waveforms for every event
NOISE=0 //flat thermal noise
NOISE_CHANNEL_MODE=0 //using different noise temperature for each channel
NOISE_EVENTS=16 // number of noise events which will be store in the trigger class for later use
ANTENNA_MODE=1
APPLY_NOISE_FIGURE=0
|