#include "TMath.h" bool NPE_EventRate_MuTauNeutrino_AtmosMuon_BaseCut(){ //options bool tau = false;//this is true only with 80 strings but can be false with 80 strings int nstring = 9; if(nstring==9)tau=false; bool integralrate = false; bool use_firstguess=false;//if this is false, mc truth cos theta is used bool plot_firstguess=false;//if this is false, plotted as mc truth cos theta bool doWeight=true; //gSystem->Load("libweighting-module"); /*style*/ SetStyle(); gROOT->ForceStyle(); //minor option for now bool dowrite = false; string indir = "/home/aya/IceCube/ehe/March06/RootReducedLeaves/"; if(nstring==80){ string file_slope2 = "PhysicsWeightedWiscE-2.MergedReduced.root"; string file_slope1 = "PhysicsWeightedE-1.MergedReduced.root"; if(tau)string file_slope1_tau = "JulietTauGZKMergedAll.reduced.root"; } else if(nstring==9){ string file_slope1 = "JulietWeighedWiscE1-S9.MergedReduced.root"; string file_slope2 = "JulietWeighedWiscE2-S9.MergedReduced.root"; if(tau)string file_slope1_tau = ""; } /* Read in E^-1 file */ const char *filenameE1 = (const *char)(indir+file_slope1).c_str(); TFile *fileE1 = (TFile *) TFile::Open(filenameE1); TTree *E1Tree = (TTree *) fileE1->Get("JulietTree"); /* Read in E^-2 file */ const char *filenameE2 = (const *char)(indir+file_slope2).c_str(); TFile *fileE2 = (TFile *) TFile::Open(filenameE2); TTree *E2Tree = (TTree *) fileE2->Get("JulietTree"); if(tau){ /* Read in E^-1 Tau file */ const char *filenameE1Tau = (const *char)(indir+file_slope1_tau).c_str(); TFile *fileE1Tau = (TFile *) TFile::Open(filenameE1Tau); TTree *E1TauTree = (TTree *) fileE1Tau->Get("JulietTree"); } cout << "*************************************************" << endl; cout <<"Energy E-1 file " << filenameE1 << " is " << E1Tree->GetEntries() << endl; cout <<"Energy E-2 file " << filenameE2 << " is " << E2Tree->GetEntries() << endl; if(tau)cout <<"Energy E-1 Tau file " << filenameE1Tau << " is " << E1TauTree->GetEntries() << endl; cout << "Number of string is " << nstring << endl; cout << "*************************************************" << endl; //Define lowest level cuts --- Our Base Cut TCut minimum_dom = "DetectorResponseEvent_.numberOfDomWithRecoHit_>4"; TCut minimum_npe = "DetectorResponseEvent_.totalBestEstimatedNPE_>200"; TCut minimum_fgcostheta = "firstguessTrack_.cosTheta_>-5"; TCut thres_nch = "DetectorResponseEvent_.numberOfDomWithRecoHit_>=53"; TCut thres_npe = "DetectorResponseEvent_.totalBestEstimatedNPE_>10000"; TCut basecut = (minimum_dom&&minimum_npe); TCut basebasecut = "DetectorResponseEvent_.numberOfDomWithRecoHit_>0&&DetectorResponseEvent_.totalBestEstimatedNPE_>0"; //TCut basecut = (minimum_npe); //TCut basecut = (minimum_dom); //Propagation matrix TCut prop_gzkNeutrino = "generatorWeight_.julietPropagationMatrixNeutrinoFlux_"; TCut prop_atmMuon = "generatorWeight_.julietPropagationMatrixAtmMuonFlux_"; //normalization const. for energy we fixed energy bin size to log(Delta E) = 0.1 TCut energy_bin_norm_slope1 = "(generatorWeight_.maxEnergyLog_-generatorWeight_.minEnergyLog_)/0.1"; TCut energy_bin_norm_slope2 = "(pow(10.0,generatorWeight_.minEnergyLog_*(1.0-generatorWeight_.powerLawIndex_))-pow(10.0,generatorWeight_.maxEnergyLog_*(1.0-generatorWeight_.powerLawIndex_)))*pow(generatorPrimary_.energy_,(generatorWeight_.powerLawIndex_-1.0))/(0.1*log(10.0))"; TCut entry_norm = "1/Entries$"; TCut entry = "Entries$";//E1Tree->Draw(entry); //normalization const for geometry TCut area_norm = "TMath::Pi()*generatorWeight_.injectionSurfaceR_*generatorWeight_.injectionSurfaceR_";// [m^2] TCut angle_norm = "4.0*TMath::Pi()";// [sr] TCut sec_year = "365.0*24.0*60.0*60.0"; // [sec/year] TCut cm2_m2 = "1.0e4"; //[cm^2/m^2] TCut m2_km2 = "1.0e6"; //[m^2/km^2] TCut km2_m2 = "1.0e-6"; //[km^2/m^2] TCut e20 = "1.0e-20"; //[] TCut signal_npe = "(DetectorResponseEvent_.totalBestEstimatedNPE_ > 5.0e5)"; if(use_firstguess)TCut signal_angle = "(firstguessTrack_.cosTheta_<0.0 && DetectorResponseEvent_.totalBestEstimatedNPE_ > 1.0e4 )"; //else TCut signal_angle = "(generatorPrimary_.cosTheta_<0.0 && DetectorResponseEvent_.totalBestEstimatedNPE_ > 1.0e4 )"; else TCut signal_angle = "(generatorPrimary_.cosTheta_<0.1 && DetectorResponseEvent_.totalBestEstimatedNPE_ > 1.0e5)"; //normalized neutrino and atmospheric muonw weights TCut signalcut = signal_npe||signal_angle; TCut effArea_slope1 = energy_bin_norm_slope1*entry_norm*area_norm*km2_m2*angle_norm; TCut effArea_slope2 = energy_bin_norm_slope2*entry_norm*area_norm*km2_m2*angle_norm; TCut rate_neutrino_slope1 = effArea_slope1*prop_gzkNeutrino*m2_km2*cm2_m2*sec_year; TCut rate_neutrino_slope2 = effArea_slope2*prop_gzkNeutrino*m2_km2*cm2_m2*sec_year; TCut rate_atmMuon_slope1 = effArea_slope1*prop_atmMuon*m2_km2*cm2_m2*sec_year; TCut rate_atmMuon_slope2 = effArea_slope2*prop_atmMuon*m2_km2*cm2_m2*sec_year; TCut rate_slope1 = entry_norm*area_norm*km2_m2*angle_norm*e20*m2_km2*cm2_m2*sec_year; TCut rate_slope2 = entry_norm*area_norm*km2_m2*angle_norm*e20*m2_km2*cm2_m2*sec_year; //Define atmospherlic muon distribution //y-axis double rate_min = 0.0; double rate_max = 0.5; //x-axis double npe_min_log = 2; double npe_max_log = 8; double npe_bin_size = 0.05; int npe_nbin = (int)(npe_max_log - npe_min_log)/npe_bin_size; //int npe_nbin = 88; TH1D* hist_E1_EventRate = new TH1D("hist_E1_EventRate","Event Rate with signal domain cut", npe_nbin, npe_min_log, npe_max_log); TH1D* hist_E2_EventRate = new TH1D("hist_E2_EventRate","atmospheric Muon Event Rate from E^{-2}, func. of Primary Muon Energy", npe_nbin, npe_min_log, npe_max_log); hist_E1_EventRate->SetLineColor(kBlue); hist_E2_EventRate->SetLineColor(kRed); hist_E1_EventRate->SetLineWidth(1); hist_E2_EventRate->SetLineWidth(1); TAxis *xaxis, *yaxis; xaxis = hist_E1_EventRate->GetXaxis(); yaxis = hist_E1_EventRate->GetYaxis(); yaxis->SetRangeUser(rate_min, rate_max); yaxis->SetTitle("Event Rate [a.u.]"); xaxis->SetTitle("log_{10} NPE"); xaxis = hist_E2_EventRate->GetXaxis(); yaxis = hist_E2_EventRate->GetYaxis(); yaxis->SetRangeUser(rate_min, rate_max); xaxis->SetTitle("log_{10} NPE"); yaxis->SetTitle("Event Rate [a.u.]"); //Start Drawing TCanvas *canvas_EventRate = new TCanvas("EventRate","EventRate", 900,500); canvas_EventRate->Divide(2, 0.001,0.001); canvas_EventRate->SetHighLightColor(10); //(basecut&&contained)*neutrino_slope ///////////////// canvas_EventRate->cd(1); SetgPad(); ///////////////// E1Tree->SetLineColor(kBlue); E2Tree->SetLineColor(kRed); E1Tree->SetLineWidth(1); E2Tree->SetLineWidth(1); E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut)*rate_atmMuon_slope2, ""); TH1D *h1 = (TH1D*)gPad->GetPrimitive("htemp"); h1->SetTitle(""); TAxis *xaxis = htemp->GetXaxis(); TAxis *yaxis = htemp->GetYaxis(); yaxis->SetTitle("Event Rate [year^{-1}] with basecuts"); xaxis->SetTitle("log_{10} NPE"); //E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut)*rate_atmMuon_slope2, "same"); E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut)*rate_neutrino_slope1, "same"); E2Tree->SetLineStyle(2); E1Tree->SetLineStyle(2); E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut&&thres_nch)*rate_atmMuon_slope2, "same"); E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut&&thres_nch)*rate_neutrino_slope1, "same"); E2Tree->SetLineStyle(3); E1Tree->SetLineStyle(3); E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut&&thres_nch&&thres_npe)*rate_atmMuon_slope2, "same"); E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_)", (basecut&&thres_nch&&thres_npe)*rate_neutrino_slope1, "same"); ///////////////// canvas_EventRate->cd(2); SetgPad(); ///////////////// E2Tree->SetLineStyle(1); E1Tree->SetLineStyle(1); E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1_EventRate", (basecut)*rate_slope1, ""); E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E2_EventRate", (basecut)*rate_slope2, "same"); E2Tree->SetLineStyle(2); E1Tree->SetLineStyle(2); E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1_EventRate1", (basecut&&thres_nch)*rate_slope1, "same"); E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E2_EventRate1", (basecut&&thres_nch)*rate_slope2, "same"); E2Tree->SetLineStyle(3); E1Tree->SetLineStyle(3); E1Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E1_EventRate2", (basecut&&thres_nch&&thres_npe)*rate_slope1, "same"); E2Tree->Draw("log10(DetectorResponseEvent_.totalBestEstimatedNPE_) >> hist_E2_EventRate2", (basecut&&thres_nch&&thres_npe)*rate_slope2, "same"); ///////////////////////////////////////////////////////////////////////////////////////// //// output to root file if you want to /////// if(dowrite) { string outfile = "outhist.root"; TFile *hist_outfile = new TFile(outfile.c_str(),"recreate"); hist_outfile->Close(); } cout << "exit" << endl; } //////////////////////////////// void SetStyle(){ gStyle->SetOptStat(1111); gStyle->SetPadGridX(true); gStyle->SetPadGridY(true); gStyle->SetStatW(0.4); gStyle->SetStatFontSize(0.03); gStyle->SetTitleFontSize(18); gStyle->SetLabelSize(0.03,"x"); gStyle->SetLabelSize(0.03,"y"); gStyle->SetLabelSize(0.03,"z"); gStyle->SetLabelOffset(0.002,"x"); gStyle->SetNdivisions(1025,"x"); gStyle->SetNdivisions(1010,"y"); gStyle->SetTitleOffset(1.4,"y"); gStyle->SetTitleOffset(0.8,"x"); gStyle->SetOptStat(0000); gStyle->SetTitleBorderSize(0); gStyle->SetTitleW(0.6); gStyle->SetTitleH(0.06); gStyle->SetOptTitle(0); gStyle->SetPalette(1); } void SetgPad(){ double right = 0.12; //double right = 0.05; double left = 0.12; double top = 0.1; double bottom = 0.08; gPad->SetRightMargin(right); gPad->SetLeftMargin(left); gPad->SetTopMargin(top); gPad->SetBottomMargin(bottom); }