| Attachment 2: |
veff_aeff2.C
//////////////////////////////////////
//To calculate various veff and aeff
//Run as ./veff_aeff2 $RADIUS $DEPTH $FILE
//For example ./veff_aeff2 8000 3000 /data/user/ypan/bin/AraSim/trunk/outputs/AraOut.setup_single_station_2_energy_20.run0.root
//
/////////////////////////////////////
#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include <time.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 <unistd.h>
#include "TVector3.h"
#include "TRotation.h"
#include "TSpline.h"
//#include "TObject.h"
#include "Tools.h"
#include "Constants.h"
#include "Vector.h"
#include "Position.h"
#include "EarthModel.h"
#include "IceModel.h"
#include "Efficiencies.h"
#include "Spectra.h"
#include "Event.h"
#include "Trigger.h"
#include "Detector.h"
#include "Settings.h"
#include "counting.hh"
#include "Primaries.h"
#include "signal.hh"
#include "secondaries.hh"
#include "Ray.h"
#include "RaySolver.h"
#include "Report.h"
using namespace std;
class EarthModel;
int main(int argc, char **argv){
//string readfile;
//readfile = string(argv[1]);
//readfile = "/data/user/ypan/bin/AraSim/branches/new_geom/outputs/AraOut20.0.root";
int ifile = 0;
double totweightsq;
double totweight;
int totnthrown;
int typetotnthrown[12];
double tottrigeff;
double sigma[12];
double typetotweight[12];
double typetotweightsq[12];
double totsigmaweight;
double totsigmaweightsq;
double volradius;
double voldepth;
const int nstrings = 9;
const int nantennas = 1;
double veff1 = 0.0;
double veff2 = 0.0;
double veff3 = 0.0;
double weight1 = 0.0;
double weight2 = 0.0;
double weight3 = 0.0;
double veffT[6], vefferrT[6], aeffT[6], aefferrT[6];
double veffF[3], vefferrF[3], aeffF[3], aefferrF[3];
double veffNu[2], vefferrNu[2], aeffNu[2], aefferrNu[2];
double veff, vefferr, aeff, aefferr, aeff2;
double pnu;
Detector *detector = 0;
//Settings *settings = 0;
//IceModel *icemodel = 0;
Event *event = 0;
Report *report = 0;
cout<<"construct detector"<<endl;
//TFile *AraFile=new TFile(readfile.c_str());
//TFile *AraFile=new TFile((outputdir+"/AraOut.root").c_str());
TChain *AraTree = new TChain("AraTree");
TChain *AraTree2 = new TChain("AraTree2");
TChain *eventTree = new TChain("eventTree");
//AraTree->SetBranchAddress("detector",&detector);
//AraTree->SetBranchAddress("settings",&settings);
//AraTree->SetBranchAddress("icemodel",&icemodel);
cout << "trees set" << endl;
for(ifile = 3; ifile < (argc - 1); ifile++){
AraTree->Add(string(argv[ifile + 1]).c_str());
AraTree2->Add(string(argv[ifile + 1]).c_str());
eventTree->Add(string(argv[ifile + 1]).c_str());
}
AraTree2->SetBranchAddress("event",&event);
AraTree2->SetBranchAddress("report",&report);
cout<<"branch detector"<<endl;
for(int i=0; i<12; i++) {
typetotweight[i] = 0.0;
typetotweightsq[i] = 0.0;
typetotnthrown[i] = 0;
}
totweightsq = 0.0;
totweight = 0.0;
totsigmaweight = 0.0;
totsigmaweightsq = 0.0;
totnthrown = AraTree2->GetEntries();
cout << "Total number of events: " << totnthrown << endl;
//totnthrown = settings->NNU;
//volradius = settings->POSNU_RADIUS;
volradius = atof(argv[1]);
voldepth = atof(argv[2]);
AraTree2->GetEntry(0);
pnu = event->pnu;
cout << "Energy " << pnu << endl;
for(int iEvt2=0; iEvt2<totnthrown; iEvt2++) {
AraTree2->GetEntry(iEvt2);
double sigm = event->Nu_Interaction[0].sigma;
int iflavor = (event->nuflavorint)-1;
int inu = event->nu_nubar;
int icurr = event->Nu_Interaction[0].currentint;
sigma[inu+2*icurr+4*iflavor] = sigm;
typetotnthrown[inu+2*icurr+4*iflavor]++;
if( (iEvt2 % 10000 ) == 0 ) cout << "*";
if(report->stations[0].Global_Pass<=0) continue;
double weight = event->Nu_Interaction[0].weight;
if(weight > 1.0){
cout << weight << "; " << iEvt2 << endl;
continue;
}
// cout << weight << endl;
totweightsq += pow(weight,2);
totweight += weight;
typetotweight[inu+2*icurr+4*iflavor] += weight;
typetotweightsq[inu+2*icurr+4*iflavor] += pow(weight,2);
totsigmaweight += weight*sigm;
totsigmaweightsq += pow(weight*sigm,2);
int trig1 = 0;
int trig2 = 0;
int trig3 = 0;
for (int i = 0; i < nstrings; i++){
if (i == 0 && report->stations[0].strings[i].antennas[0].Trig_Pass > 0) trig1++;
if (i > 0 && i < 5 && report->stations[0].strings[i].antennas[0].Trig_Pass > 0) trig2++;
if (i > 4 && i < 9 && report->stations[0].strings[i].antennas[0].Trig_Pass > 0) trig3++;
}
if ( trig1 > 0)//phase array
weight1 += event->Nu_Interaction[0].weight;
if ( trig2 > 1)//lpda
weight2 += event->Nu_Interaction[0].weight;
if (trig3 > 3)//bicone
weight3 += event->Nu_Interaction[0].weight;
}
tottrigeff = totweight / double(totnthrown);
double nnucleon = 5.54e29;
double vtot = PI * double(volradius) * double(volradius) * double(voldepth) / 1e9;
veff = vtot * totweight / double(totnthrown) * 4.0 * PI;
//vefferr = sqrt(SQ(sqrt(double(totnthrown))/double(totnthrown))+SQ(sqrt(totweightsq)/totweight));
vefferr = sqrt(totweightsq) / totweight * veff;
aeff = vtot * (1e3) * nnucleon * totsigmaweight / double(totnthrown);
//aefferr = sqrt(SQ(sqrt(double(totnthrown))/double(totnthrown))+SQ(sqrt(totsigmaweightsq)/totsigmaweight));
//aefferr = sqrt(SQ(sqrt(double(totnthrown))/double(totnthrown))+SQ(sqrt(totweightsq)/totweight));
aefferr = sqrt(totweightsq) / totweight * aeff;
double sigmaave = 0.0;
for(int iflavor=0; iflavor<3; iflavor++) {
double flavorweight = 0.0;
double flavorweightsq = 0.0;
double flavorsigmaave = 0.0;
int flavortotthrown = 0;
double temptotweightnu[2] = {0};
double tempsignu[2] = {0};
double temptotweight = 0.0;
for(int inu=0; inu<2; inu++) {
double tempsig = 0.0;
double tempweight = 0.0;
for(int icurr=0; icurr<2; icurr++) {
tempsig += sigma[inu+2*icurr+4*iflavor];
tempsignu[inu] += sigma[inu+2*icurr+4*iflavor];
tempweight += typetotweight[inu+2*icurr+4*iflavor];
flavorweight += typetotweight[inu+2*icurr+4*iflavor];
flavorweightsq += typetotweightsq[inu+2*icurr+4*iflavor];
temptotweight += typetotweight[inu+2*icurr+4*iflavor];
temptotweightnu[inu] += typetotweight[inu+2*icurr+4*iflavor];
flavortotthrown += typetotnthrown[inu+2*icurr+4*iflavor];
}
//printf("Temp Sigma: "); cout << tempsig << "\n";
sigmaave += tempsig*(tempweight/totweight);
}
flavorsigmaave += tempsignu[0]*(temptotweightnu[0]/temptotweight)+tempsignu[1]*(temptotweightnu[1]/temptotweight);
veffF[iflavor] = vtot*flavorweight/double(flavortotthrown);
vefferrF[iflavor] = sqrt(flavorweightsq)/flavorweight;
//printf("Volume: %.9f*%.9f/%.9f \n",vtot,flavorweight,double(totnthrown));
aeffF[iflavor] = veffF[iflavor]*(1e3)*nnucleon*flavorsigmaave;
aefferrF[iflavor] = sqrt(flavorweightsq)/flavorweight;
}
for(int inu=0; inu<2; inu++) {
double tempsig = 0.0;
double tempweight = 0.0;
double tempweightsq = 0.0;
int nutotthrown = 0;
for(int iflavor=0; iflavor<3; iflavor++) {
for(int icurr=0; icurr<2; icurr++) {
tempweight += typetotweight[inu+2*icurr+4*iflavor];
tempweightsq += typetotweightsq[inu+2*icurr+4*iflavor];
nutotthrown += typetotnthrown[inu+2*icurr+4*iflavor];
}
}
tempsig += sigma[inu+2*0+4*0];
tempsig += sigma[inu+2*1+4*0];
veffNu[inu] = vtot*tempweight/double(nutotthrown);
vefferrNu[inu] = sqrt(tempweightsq)/tempweight;
aeffNu[inu] = veffNu[inu]*(1e3)*nnucleon*tempsig;
aefferrNu[inu] = sqrt(tempweightsq)/tempweight;
}
double totalveff = 0.0;
double totalaeff = 0.0;
for(int inu=0; inu<2; inu++) {
for(int iflavor=0; iflavor<3; iflavor++) {
int typetotthrown = 0;
for(int icurr=0; icurr<2; icurr++) {
typetotthrown += typetotnthrown[inu+2*icurr+4*iflavor];
}
totalveff += veffT[iflavor+3*inu]*(double(typetotthrown)/double(totnthrown));
totalaeff += aeffT[iflavor+3*inu]*(double(typetotthrown)/double(totnthrown));
}
}
aeff2 = veff*(1e3)*nnucleon*sigmaave;
aeff = aeff2;
veff1 = weight1 / totnthrown * vtot * 4.0 * PI;
veff2 = weight2 / totnthrown * vtot * 4.0 * PI;
veff3 = weight3 / totnthrown * vtot * 4.0 * PI;
printf("\nvolthrown: %.6f; totweight: %.6f; Veff: %.6f +- %.6f\n", vtot, totweight, veff, vefferr);
printf("veff1: %.3f; veff2: %.3f; veff3: %.3f\n", veff1, veff2, veff3);
//string des = string(argv[4]);
char buf[100];
std::ostringstream stringStream;
stringStream << string(argv[3]);
std::string copyOfStr = stringStream.str();
snprintf(buf, sizeof(buf), "Veff_des%s.txt", copyOfStr.c_str());
FILE *fout = fopen(buf, "a+");
fprintf(fout, "%e, %.6f, %.6f, %.3f, %.3f, %.3f \n", pnu, veff, vefferr, veff1, veff2, veff3);
fclose(fout);
return 0;
}
|
| Attachment 4: |
make_plot.py
# -*- coding: utf-8 -*-
import numpy as np
import sys
import matplotlib.pyplot as plt
from pylab import setp
from matplotlib.pyplot import rcParams
import csv
import pandas as pd
rcParams['mathtext.default'] = 'regular'
def read_file(finame):
fi = open(finame, 'r')
rdr = csv.reader(fi, delimiter=',', skipinitialspace=True)
table = []
for row in rdr:
# print(row)
energy = float(row[0])
veff = float(row[1])
veff_err = float(row[2])
veff1 = float(row[3])
veff2 = float(row[4])
veff3 = float(row[5])
row = {'energy':energy, 'veff':veff, 'veff_err':veff_err, 'veff1':veff1, 'veff2':veff2, 'veff3':veff3}
table.append(row)
df=pd.DataFrame(table)
df_ordered=df.sort_values('energy',ascending=True)
# print(df_ordered)
return df_ordered
def beautify_veff(this_ax):
sizer=20
xlow = 1.e16 #the lower x limit
xup = 2.e20 #the uppper x limit
ylow =1e-3 #the lower x limit
yup = 6.e1 #the uppper x limit
this_ax.set_xlabel('Energy [eV]',size=sizer) #give it a title
this_ax.set_ylabel('[V$\Omega]_{eff}$ [km$^3$sr]',size=sizer)
this_ax.set_yscale('log')
this_ax.set_xscale('log')
this_ax.tick_params(labelsize=sizer)
this_ax.set_xlim([xlow,xup]) #set the x limits of the plot
this_ax.set_ylim([ylow,yup]) #set the y limits of the plot
this_ax.grid()
this_legend = this_ax.legend(loc='upper left')
setp(this_legend.get_texts(), fontsize=17)
setp(this_legend.get_title(), fontsize=17)
def main():
"""
arasim_energies = np.array([3.16e+16, 1e+17, 3.16e+17, 1e+18, 3.16e+18, 1e+19, 3.16e+19, 1e+20])
arasim_energies2 = np.array([3.16e+16, 1e+17, 3.16e+17, 1e+18, 1e+19, 3.16e+19, 1e+20])
#arasim_desA_veff = np.array([0.080894,0.290695,0.943223,2.388708,4.070498,6.824112,10.506490,13.969418])
arasim_desA_veff_s = np.array([0.067384,0.289591,0.996509,2.464464,4.945600,8.735506,13.357300,18.751915])
# arasim_desA_veff_ice = np.array([0.066291,0.303620,0.927647,2.427554,4.962093,8.465895,13.425852,18.706528])
arasim_desA_error = np.array([0.008367,0.017401,0.032222,0.084223,0.119240,0.221266,0.272460,0.320456])
# arasim_desA_error_ice = np.array([0.008345,0.017748,0.031295,0.083747,0.119370,0.217717,0.273075,0.319872])
# arasim_desB_veff = np.array([0.124111,0.417102,1.310555,3.648494,4.070E+0,11.675189,18.393961,13.688909])
arasim_desB_veff_s = np.array([0.064937,0.355747,1.289947,3.821705,8.002805,15.352981,25.391282,18.009545])
# arasim_desB_veff_ice = np.array([0.080070,0.340927,1.369081,3.938550,8.211407,15.190858,25.066541,18.021033])
arasim_desB_error = np.array([0.008268,0.019317,0.036926,0.105445,0.152488,0.296617,0.377997,0.314314])
# arasim_desB_error_ice = np.array([0.009220,0.019144,0.037966,0.107189,0.154430,0.293728,0.375827,0.314244])
arasim_desA_veff1_sm = np.array([0.054,0.231,0.868,2.239,4.586,8.339,12.910,18.301])
arasim_desA_veff2_sm = np.array([0.009,0.036,0.098,0.323,0.705,1.167,2.169,3.130])
arasim_desA_veff3_sm = np.array([0.011,0.074,0.193,0.664,1.208,2.347,3.689,4.814])
arasim_desA_veff1 = np.array([0.053,0.282,1.108,3.462,7.486,14.613,24.682,17.557])
arasim_desA_veff2 = np.array([0.006,0.040,0.121,0.322,0.636,1.308,1.961,2.708])
arasim_desA_veff3 = np.array([0.006,0.065,0.188,0.638,1.187,2.270,3.301,4.322])
"""
veff_A=read_file("Veff_desA.txt")
veff_B=read_file("Veff_desB.txt")
print("desA is \n", veff_A)
print("desB is \n", veff_B)
dfA = veff_A[['veff','veff_err']]
dfB = veff_B[['veff','veff_err']]
dfA.to_csv('veffA.csv', sep='\t',index=False)
dfB.to_csv('veffB.csv', sep='\t',index=False)
fig = plt.figure(figsize=(11,8.5))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(veff_A['energy'], veff_A['veff'],'bs-',label='Strawman',markersize=8,linewidth=2)
ax1.plot(veff_B['energy'], veff_B['veff'],'gs-',label='Punch @100 m',markersize=8,linewidth=2)
ax1.fill_between(veff_A['energy'], veff_A['veff']-veff_A['veff_err'], veff_A['veff']+veff_A['veff_err'], alpha=0.2, color='red')
ax1.fill_between(veff_B['energy'], veff_B['veff']-veff_B['veff_err'], veff_B['veff']+veff_B['veff_err'], alpha=0.2, color='red')
beautify_veff(ax1)
ax1.set_title("Punch vs Strawman, noise + signal",fontsize=20)
fig.savefig("desAB.png",edgecolor='none',bbox_inches="tight") #save the figure
fig2 = plt.figure(figsize=(11,8.5))
ax2 = fig2.add_subplot(1,1,1)
#Triggers plot
ax2.plot(veff_B['energy'], veff_B['veff1'],'g^-',label='Phased array (Punch)',markersize=8,linewidth=2)
ax2.plot(veff_A['energy'], veff_A['veff1'],'gs-',label='Phased array (Strawman)',markersize=8,linewidth=2)
ax2.plot(veff_B['energy'], veff_B['veff2'],'b^-',label='LPDAs (Punch)',markersize=8,linewidth=2)
ax2.plot(veff_A['energy'], veff_A['veff2'],'bs-',label='LPDAs (Strawman)',markersize=8,linewidth=2)
ax2.plot(veff_B['energy'], veff_B['veff3'],'y^-',label='Bicones (Punch)',markersize=8,linewidth=2)
ax2.plot(veff_A['energy'], veff_A['veff3'],'ys-',label='Bicones (Strawman)',markersize=8,linewidth=2)
ax2.set_title("Triggers contribution, noise + signal",fontsize=20)
beautify_veff(ax2)
fig2.savefig("desAB_triggers.png",edgecolor='none',bbox_inches="tight") #save the figure
main()
|