#include "TMath.h" bool PrimaryEnergy_EffArea_CosTheta_MuTauNeutrino_AtmosMuon(){ //options int nstring = 80; bool tau = false;//this is true only with 80 strings but can be false with 80 strings if(nstring==9&&tau)tau=false; //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 basecut = (minimum_dom&&minimum_npe); //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.2 TCut energy_bin_norm_slope1 = "(generatorWeight_.maxEnergyLog_-generatorWeight_.minEnergyLog_)/0.2"; 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.2*log(10.0))"; TCut entry_norm = "10/(Entries$)"; TCut entry = "Entries$";//E1Tree->Draw(entry); //normalization const for geometry TCut area_norm = "TMath::Pi()*generatorWeight_.injectionSurfaceR_*generatorWeight_.injectionSurfaceR_/1.0e6";// [km^2] //TCut angle_norm = "4.0*TMath::Pi()";// [sr] TCut angle_norm = "";// [] 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 signal_npe = "(DetectorResponseEvent_.totalBestEstimatedNPE_ > 5.0e5)"; TCut signal_angle = "(generatorPrimary_.cosTheta_<0.1 && DetectorResponseEvent_.totalBestEstimatedNPE_ > 1.0e5)"; TCut angle_bin1 = "generatorTrack_.cosTheta_>-1.0 && generatorTrack_.cosTheta_ <= -0.8"; TCut angle_bin2 = "generatorTrack_.cosTheta_>-0.8 && generatorTrack_.cosTheta_ <= -0.6"; TCut angle_bin3 = "generatorTrack_.cosTheta_>-0.6 && generatorTrack_.cosTheta_ <= -0.4"; TCut angle_bin4 = "generatorTrack_.cosTheta_>-0.4 && generatorTrack_.cosTheta_ <= -0.2"; TCut angle_bin5 = "generatorTrack_.cosTheta_>-0.2 && generatorTrack_.cosTheta_ <= 0.0"; TCut angle_bin6 = "generatorTrack_.cosTheta_> 0.0 && generatorTrack_.cosTheta_ <= 0.2"; TCut angle_bin7 = "generatorTrack_.cosTheta_> 0.2 && generatorTrack_.cosTheta_ <= 0.4"; TCut angle_bin8 = "generatorTrack_.cosTheta_> 0.4 && generatorTrack_.cosTheta_ <= 0.6"; TCut angle_bin9 = "generatorTrack_.cosTheta_> 0.6 && generatorTrack_.cosTheta_ <= 0.8"; TCut angle_bin10 = "generatorTrack_.cosTheta_> 0.8 && generatorTrack_.cosTheta_ <= 1.0"; //normalized neutrino and atmospheric muonw weights TCut rate_neutrino_slope1 = prop_gzkNeutrino*energy_bin_norm_slope1*entry_norm*area_norm*angle_norm*sec_year; TCut rate_neutrino_slope2 = prop_gzkNeutrino*energy_bin_norm_slope2*entry_norm*area_norm*angle_norm*sec_year; TCut rate_atmMuon_slope1 = prop_atmMuon*energy_bin_norm_slope1*entry_norm*area_norm*angle_norm*sec_year; TCut rate_atmMuon_slope2 = prop_atmMuon*energy_bin_norm_slope2*entry_norm*area_norm*angle_norm*sec_year; TCut signalcut = signal_npe||signal_angle; TCut effArea_slope1 = energy_bin_norm_slope1*entry_norm*area_norm*angle_norm; TCut effArea_slope2 = energy_bin_norm_slope2*entry_norm*area_norm*angle_norm; TCut rate_neutrino_slope1 = effArea_slope1*prop_gzkNeutrino*sec_year; TCut rate_neutrino_slope2 = effArea_slope2*prop_gzkNeutrino*sec_year; TCut rate_atmMuon_slope1 = effArea_slope1*prop_atmMuon*sec_year; TCut rate_atmMuon_slope2 = effArea_slope2*prop_atmMuon*sec_year; //Define atmospherlic muon distribution //y-axis double area_min = 0; double area_max = 3; //x-axis double energy_min_log = 4.; double energy_max_log = 12.; double energy_bin_size = 0.2; int energy_nbin = (int)(energy_max_log - energy_min_log)/energy_bin_size; TF1 *highnpe = new TF1("highnpe","log10(5*pow(10,5))",0,100);highnpe->SetLineWidth(1);highnpe->SetLineStyle(1); TF1 *horizon = new TF1("horizon","5",0,100);horizon->SetLineWidth(1);horizon->SetLineStyle(2); TH1D* hist_E1_EffArea1 = new TH1D("hist_E1_EffArea1","Effective Area from E^{-1}, func. of Primary Energy", energy_nbin, energy_min_log, energy_max_log); TH1D* hist_E2_EffArea1 = new TH1D("hist_E2_EffArea1","Effective Area from E^{-2}, func. of Primary Muon Energy", energy_nbin, energy_min_log, energy_max_log); if(tau) TH1D* hist_E1Tau_EffArea1 = new TH1D("hist_E1Tau_EffArea1","Effective Area from E^{-1}, func. of Primary Tau Energy", energy_nbin, energy_min_log, energy_max_log); TAxis *xaxis, *yaxis; xaxis = hist_E1_EffArea1->GetXaxis(); yaxis = hist_E1_EffArea1->GetYaxis(); yaxis->SetRangeUser(area_min, area_max); xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); yaxis->SetTitle("Effective Area [km^{2}]"); xaxis = hist_E2_EffArea1->GetXaxis(); yaxis = hist_E2_EffArea1->GetYaxis(); yaxis->SetRangeUser(area_min, area_max); xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); yaxis->SetTitle("Effective Area [km^{2}]"); if(tau){ xaxis = hist_E1Tau_EffArea1->GetXaxis(); yaxis = hist_E1Tau_EffArea1->GetYaxis(); yaxis->SetRangeUser(area_min, area_max); xaxis->SetTitle("log_{10} (Primary Energy [GeV]) "); yaxis->SetTitle("Effective Area [km^{2}]"); } //Start Drawing TCanvas *canvas_EffArea = new TCanvas("EffArea","EffArea", 500,500); canvas_EffArea->Divide(1, 0.001,0.001); canvas_EffArea->SetHighLightColor(10); //(basecut&&contained)*neutrino_slope1 ///////////////// canvas_EffArea->cd(1); SetgPad(); ///////////////// E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea1", (basecut&&signalcut&&angle_bin1)*effArea_slope1, ""); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea2", (basecut&&signalcut&&angle_bin2)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea3", (basecut&&signalcut&&angle_bin3)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea4", (basecut&&signalcut&&angle_bin4)*effArea_slope1, "same"); E1Tree->SetLineColor(2); E1Tree->Draw("generatorPrimary_.energyLog_", (basecut&&signalcut&&angle_bin3)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea6", (basecut&&signalcut&&angle_bin6)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea7", (basecut&&signalcut&&angle_bin7)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea8", (basecut&&signalcut&&angle_bin8)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea9", (basecut&&signalcut&&angle_bin9)*effArea_slope1, "same"); E1Tree->SetLineColor(4); E1Tree->Draw("generatorPrimary_.energyLog_", (basecut&&signalcut&&angle_bin5)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea11", (basecut&&signalcut&&angle_bin11)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea12", (basecut&&signalcut&&angle_bin12)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea13", (basecut&&signalcut&&angle_bin13)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea14", (basecut&&signalcut&&angle_bin14)*effArea_slope1, "same"); E1Tree->SetLineColor(6); E1Tree->Draw("generatorPrimary_.energyLog_", (basecut&&signalcut&&angle_bin8)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea16", (basecut&&signalcut&&angle_bin16)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea17", (basecut&&signalcut&&angle_bin17)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea18", (basecut&&signalcut&&angle_bin18)*effArea_slope1, "same"); //E1Tree->Draw("generatorPrimary_.energyLog_ >> hist_E1_EffArea19", (basecut&&signalcut&&angle_bin19)*effArea_slope1, "same"); E1Tree->SetLineColor(8); E1Tree->Draw("generatorPrimary_.energyLog_", (basecut&&signalcut&&angle_bin10)*effArea_slope1, "same"); ///////////////// /* //E2Tree->Draw("generatorPrimary_.energyLog_ >> hist_E2_EffArea1", (basecut&&signalcut)*effArea_slope2, ""); if(nstring==80&&tau){ canvas_EffArea->cd(3); SetgPad(); hist_E1Tau_EffArea1->SetLineColor(kBlue); E1TauTree->Draw("generatorPrimary_.energyLog_ >> hist_E1Tau_EffArea1", (basecut&&signalcut)*effArea_slope1, ""); } //hist_E1_EffArea1->DrawCopy(""); if(tau) hist_E1Tau_EffArea1->DrawCopy("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(1); 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); }