11 #include "TTreeIndex.h" 19 #include "icemc_random.h" 24 #include "TPostScript.h" 30 #include "TGraphErrors.h" 31 #include "TGraphAsymmErrors.h" 37 #include "TRotation.h" 39 #include "Math/InterpolationTypes.h" 40 #include "Math/Interpolator.h" 46 #include "Constants.h" 48 #include "position.hh" 50 #include "earthmodel.hh" 53 #include "roughness.hh" 56 #include "icemodel.hh" 60 #include "secondaries.hh" 62 #include "counting.hh" 63 #include "Primaries.h" 64 #include "Taumodel.hh" 66 #include "GlobalTrigger.h" 67 #include "ChanTrigger.h" 68 #include "SimulatedSignal.h" 69 #include "EnvironmentVariable.h" 75 #if __cplusplus > 199711L 76 #define isnan std::isnan 77 #include <type_traits> 84 #ifdef ANITA_UTIL_EXISTS 85 #include "UsefulAnitaEvent.h" 86 #ifdef ANITA3_EVENTCORRELATOR 89 #include "AnitaGeomTool.h" 90 #include "AnitaConventions.h" 91 #include "RawAnitaHeader.h" 94 #include "UsefulAdu5Pat.h" 98 #ifdef ANITA3_EVENTREADER 99 #include "TruthAnitaEvent.h" 121 const int NVIEWANGLE=100;
125 double eventsfound_beforetrigger=0.;
126 double eventsfound_crust=0;
127 double eventsfound_weightgt01=0;
128 double eventsfound_belowhorizon=0;
129 double eventsfound=0;
130 double eventsfound_prob=0;
131 double eventsfound_divided_by_lint = 0;
135 double poissonerror_minus[21] = {0.-0.00, 1.-0.37, 2.-0.74, 3.-1.10, 4.-2.34, 5.-2.75, 6.-3.82, 7.-4.25, 8.-5.30, 9.-6.33, 10.-6.78, 11.-7.81, 12.-8.83, 13.-9.28, 14.-10.30, 15.-11.32, 16.-12.33, 17.-12.79, 18.-13.81, 19.-14.82, 20.-15.83};
136 double poissonerror_plus[21] = {1.29-0., 2.75-1., 4.25-2., 5.30-3., 6.78-4., 7.81-5., 9.28-6., 10.30-7., 11.32-8., 12.79-9., 13.81-10., 14.82-11., 16.29-12., 17.30-13., 18.32-14., 19.32-15., 20.80-16., 21.81-17., 22.82-18., 23.82-19., 25.30-20.};
139 double MIN_LOGWEIGHT=-3;
140 double MAX_LOGWEIGHT=-1;
142 double eventsfound_binned[NBINS];
143 double eventsfound_binned_e[NBINS];
144 double eventsfound_binned_mu[NBINS];
145 double eventsfound_binned_tau[NBINS];
152 double error_e_plus=0;
153 double error_mu_plus=0;
154 double error_tau_plus=0;
155 double error_minus=0;
156 double error_e_minus=0;
157 double error_mu_minus=0;
158 double error_tau_minus=0;
160 double gain_dipole=2.15;
161 double changle_deg=0;
165 double RANDOMISEPOL=0.;
169 std::vector<int> NNU_per_source;
170 std::vector<std::string> source_names;
171 std::vector<double> eventsfound_per_source;
172 std::vector<double> eventsfound_prob_per_source;
174 double volume_thishorizon;
175 double realtime_this;
176 double longitude_this;
177 double latitude_this;
178 double altitude_this;
179 double heading_this=0.;
184 double pnu=pow(10., 20);
189 double SIGNALRADIUS=2.;
192 const double FREQ_LOW_DISCONES=120.E6;
193 const double FREQ_HIGH_DISCONES=1000.E6;
195 double bwslice_vnoise_thislayer[4];
196 int passes_thisevent=0;
197 int unmasked_thisevent=0;
199 int discones_passing;
201 double heff_discone=0;
202 double polarfactor_discone=0.;
204 double volts_discone=0.;
205 double vnoise_discone=0.;
208 double BW_DISCONES=300.E6-120.E6;
217 double t_coeff_pokey, t_coeff_slappy;
220 double e_comp_max1=0;
221 double h_comp_max1=0;
222 double e_comp_max2=0;
223 double h_comp_max2=0;
224 double e_comp_max3=0;
225 double h_comp_max3=0;
228 double diff_3tries=0;
231 double costheta_inc=0;
232 double costheta_exit=0;
233 double theta_rf_atbn;
234 double theta_rf_atbn_measured;
235 double theta_nnu_atbn;
236 double theta_nnu_rf_diff_atbn;
240 double phi_nnu_rf_diff_atbn;
248 double costhetanu=-1000;
253 double nearthlayers=0;
254 double weight_prob=0.;
259 double pieceofkm2sr=0;
262 int count_inthisloop1=0;
263 int count_inthisloop2=0;
264 int count_inthisloop3=0;
265 double averaging_thetas1=0;
266 double averaging_thetas2=0;
267 double averaging_thetas3=0;
272 int count_asktrigger=0;
273 int count_asktrigger_nfb=0;
276 double count_passestrigger_w=0;
279 int allcuts[2]={0, 0};
281 double allcuts_weighted[2]={0, 0};
282 double allcuts_weighted_polarization[3]={0, 0, 0};
286 int count_chanceofsurviving=0;
288 int count_chanceinhell0=0;
292 double count_chanceinhell2_w=0;
296 int count_chordgoodlength=0;
297 int count_d2goodlength=0;
301 double sum_frac_db[3];
303 const int NBINS_DISTANCE=28;
304 double eventsfound_binned_distance[NBINS_DISTANCE] = {0.};
305 int index_distance=0;
306 double km3sr_distance[NBINS_DISTANCE] = {0.};
307 double error_distance_plus[NBINS_DISTANCE] = {0.};
308 double error_distance_minus[NBINS_DISTANCE] = {0.};
309 int eventsfound_binned_distance_forerror[NBINS_DISTANCE][NBINS] = {{0}};
315 int count_passestrigger_nfb=0;
316 double percent_increase_db=0;
317 double percent_increase_nfb=0;
318 double percent_increase_total=0;
320 double error_km3sr_nfb=0;
321 double error_percent_increase_nfb=0;
326 int count_dbexitsice=0;
329 double eventsfound_nfb_binned[NBINS];
332 double heff_max=0.62639;
335 double scalefactor_distance=0;
336 double scalefactor_attenuation=0;
337 double MAX_ATTENLENGTH=1671;
340 double dviewangle_deg=0;
342 double forseckel[NVIEWANGLE][Anita::NFREQ];
343 double viewangles[NVIEWANGLE];
344 double GetAirDistance(
double altitude_bn,
double beta);
352 int max_antenna0 = -1;
353 int max_antenna1 = -1;
354 int max_antenna2 = -1;
355 double max_antenna_volts0 = 0;
356 double max_antenna_volts0_em = 0;
358 double max_antenna_volts1 = 0;
359 double max_antenna_volts2 = 0;
361 double rx0_signal_eachband[2][5];
362 double rx0_threshold_eachband[2][5];
363 double rx0_noise_eachband[2][5];
364 int rx0_passes_eachband[2][5];
371 double vmmhz1m_visible = 0;
372 int freq_bins = Anita::NFREQ;
374 double r_in_array[3];
375 double nsurf_rfexit_array[3];
376 double nsurf_rfexit_db_array[3];
377 double r_bn_array[3];
378 double n_bn_array[3];
379 double posnu_array[3];
380 double nrf_iceside_array[5][3];
381 double nrf_iceside_db_array[5][3];
382 double ant_max_normal0_array[3];
383 double ant_max_normal1_array[3];
384 double ant_max_normal2_array[3];
385 double n_pol_array[3];
386 double n_exit2bn_array[5][3];
387 double r_enterice_array[3];
388 double n_exit2bn_db_array[5][3];
389 double rfexit_array[5][3];
390 double rfexit_db_array[5][3];
392 int times_crust_entered_det=0;
393 int times_mantle_entered_det=0;
394 int times_core_entered_det=0;
395 int neutrinos_passing_all_cuts=0;
396 double sum_weights=0;
399 int xsecParam_nutype = 0;
400 int thisNuType = xsecParam_nutype;
401 int xsecParam_nuint = 1;
404 double justNoise_trig[2][48][512];
405 double justSignal_trig[2][48][512];
406 double justNoise_dig[2][48][512];
407 double justSignal_dig[2][48][512];
412 void SetupViewangles(
Signal *sig1);
414 void GetAir(
double *col1);
415 double GetThisAirColumn(
Settings*,
Position r_in,
Vector nnu,
Position posnu,
double *col1,
double& cosalpha,
double& mytheta,
double& cosbeta0,
double& mybeta);
420 double IsItDoubleBang(
double exitlength,
double plepton);
423 void GetCurrent(
string& current);
426 double GetAverageVoltageFromAntennasHit(
Settings *settings1,
int *nchannels_perrx_triggered,
double *voltagearray,
double& volts_rx_sum);
431 void Attenuate(
IceModel *antartica1,
Settings *settings1,
double& vmmhz_max,
double rflength,
const Position &posnu);
435 void IsAbsorbed(
double chord_kgm2,
double len_int_kgm2,
double& weight);
439 int GetRayIceSide(
const Vector &n_exit2rx,
const Vector &nsurf_rfexit,
double nexit,
double nenter,
Vector &nrf2_iceside);
441 int GetDirection(
Settings *settings1,
Interaction *interaction1,
const Vector &refr,
double deltheta_em,
double deltheta_had,
double emfrac,
double hadfrac,
double vmmhz1m_max,
double r_fromballoon,
Ray *ray1,
Signal *sig1,
Position posnu, Anita *anita1,
Balloon *bn1,
Vector &nnu,
double& costhetanu,
double& theta_threshold);
443 void GetFresnel(
Roughness *rough1,
int ROUGHNESS_SETTING,
const Vector &nsurf_rfexit,
const Vector &n_exit2rx,
Vector &n_pol,
const Vector &nrf2_iceside,
double efield,
double emfrac,
double hadfrac,
double deltheta_em,
double deltheta_had,
double &t_coeff_pokey,
double &t_coeff_slappy,
double &fresnel,
double &mag);
445 double GetViewAngle(
const Vector &nrf2_iceside,
const Vector &nnu);
447 int TIR(
const Vector &n_surf,
const Vector &nrf2_iceside,
double N_IN,
double N_OUT);
449 void IntegrateBands(Anita *anita1,
int k, Screen *panel1,
double *freq,
double scalefactor,
double *sumsignal);
451 void Integrate(Anita *anita1,
int j,
int k,
double *vmmhz,
double *freq,
double scalefactor,
double sumsignal);
453 void interrupt_signal_handler(
int);
455 bool ABORT_EARLY =
false;
457 void Summarize(
Settings *settings1, Anita* anita1,
Counting *count1,
Spectra *spectra1,
Signal *sig1,
Primaries *primary1,
double,
double eventsfound,
double,
double,
double,
double*,
double,
double,
double&,
double&,
double&,
double&, ofstream&, ofstream&, TString,
Balloon*,
SourceModel * src_model);
461 void CloseTFile(TFile *hfile);
463 int Getmine(
double*,
double*,
double*,
double*);
465 void Getearth(
double*,
double*,
double*,
double*);
468 #ifdef ANITA_UTIL_EXISTS 469 int GetIceMCAntfromUsefulEventAnt(
Settings *settings1,
int UsefulEventAnt);
476 double thresholdsAnt[48][2][5];
477 double thresholdsAntPass[48][2][5];
481 double threshold_start=-1.;
482 double threshold_end=-6.;
483 const int NTHRESHOLDS=20;
484 double threshold_step=(threshold_end-threshold_start)/(
double)NTHRESHOLDS;
486 double npass_v_thresh[NTHRESHOLDS]={0.};
487 double denom_v_thresh[NTHRESHOLDS]={0.};
488 double npass_h_thresh[NTHRESHOLDS]={0.};
489 double denom_h_thresh[NTHRESHOLDS]={0.};
490 double thresholds[NTHRESHOLDS];
492 int main(
int argc,
char **argv) {
503 #ifdef ICEMC_FEEXCEPT 504 feenableexcept(FE_INVALID | FE_DIVBYZERO);
508 double sumsignal[5]={0.};
509 double sumsignal_aftertaper[5]={0.};
515 string input= ICEMC_SRC_DIR +
"/inputs.conf";
519 if( (argc%3!=1)&&(argc%2!=1)) {
520 cout <<
"Syntax for options: -i inputfile -o outputdir -r run_number\n";
525 double trig_thresh=0.;
530 bool useEventList=
false;
531 ifstream eventListFile;
532 string eventListName;
535 while ((clswitch = getopt(argc, argv,
"t:i:o:r:n:e:x:l:")) != EOF) {
538 nnu_tmp=atoi(optarg);
539 cout <<
"Changed number of simulated neutrinos to " << nnu_tmp << endl;
542 trig_thresh=atof(optarg);
546 cout <<
"Changed input file to: " << input << endl;
550 cout <<
"Changed output directory to: " << string(outputdir.Data()) << endl;
551 stemp=
"mkdir -p " + string(outputdir.Data());
552 system(stemp.c_str());
555 exp_tmp=atof(optarg);
556 cout <<
"Changed neutrino energy exponent to " << exp_tmp << endl;
559 startNu=atoi(optarg);
560 cout <<
"Running icemc for just 1 event with eventNumber : " << startNu << endl;
564 eventListName=optarg;
565 cout <<
"Using event list from " << eventListName << endl;
566 eventListFile.open(eventListName.c_str());
567 if (!eventListFile.is_open()) {
568 cout <<
"This is not a valid filename " << endl;
574 stringstream convert(run_num);
581 settings1->SEED=settings1->SEED +run_no;
582 cout <<
"seed is " << settings1->SEED << endl;
585 stemp=string(outputdir.Data())+
"/nu_position"+run_num+
".txt";
586 ofstream nu_out(stemp.c_str(), ios::app);
588 stemp=string(outputdir.Data())+
"/veff"+run_num+
".txt";
589 ofstream veff_out(stemp.c_str(), ios::app);
591 stemp=string(outputdir.Data())+
"/distance"+run_num+
".txt";
592 ofstream distanceout(stemp.c_str());
594 stemp=string(outputdir.Data())+
"/debug"+run_num+
".txt";
595 fstream outfile(stemp.c_str(), ios::out);
597 stemp=string(outputdir.Data())+
"/forbrian"+run_num+
".txt";
598 fstream forbrian(stemp.c_str(), ios::out);
600 stemp=string(outputdir.Data())+
"/al_voltages_direct"+run_num+
".dat";
601 fstream al_voltages_direct(stemp.c_str(), ios::out);
603 stemp=string(outputdir.Data())+
"/events"+run_num+
".txt";
604 ofstream eventsthatpassfile(stemp.c_str());
606 stemp=string(outputdir.Data())+
"/numbers"+run_num+
".txt";
607 ofstream fnumbers(stemp.c_str());
609 stemp=string(outputdir.Data())+
"/output"+run_num+
".txt";
610 ofstream foutput(stemp.c_str(), ios::app);
612 stemp=string(outputdir.Data())+
"/slacviewangles"+run_num+
".dat";
613 ofstream fslac_viewangles(stemp.c_str());
615 stemp=string(outputdir.Data())+
"/slac_hitangles"+run_num+
".dat";
616 ofstream fslac_hitangles(stemp.c_str());
619 Anita *anita1=
new Anita();
628 settings1->ReadInputs(input.c_str(), foutput, NNU, RANDOMISEPOL);
629 settings1->ApplyInputs(anita1, sec1, sig1, bn1, ray1);
635 settings1->SEED=settings1->SEED + run_no;
636 setSeed(settings1->SEED);
638 bn1->InitializeBalloon();
639 anita1->Initialize(settings1, foutput, inu, outputdir);
647 anita1->powerthreshold[4]=trig_thresh;
651 settings1->EXPONENT=exp_tmp;
654 if (!strcasecmp(settings1->SOURCE.c_str(),
"CUSTOM"))
658 src_model->
addSource(
new Source(
"Custom Source", (settings1->CUSTOM_RA)/15, settings1->CUSTOM_DEC,
662 else if(strcasecmp(settings1->SOURCE.c_str(),
"NONE"))
666 SourceModel::Restriction::fromString(settings1->WHICH_START_TIME.c_str()),
667 SourceModel::Restriction::fromString(settings1->WHICH_END_TIME.c_str())));
670 if (settings1->SOURCE_USE_EXPONENT && settings1->EXPONENT < 22)
672 settings1->SOURCE_MIN_E=settings1->EXPONENT;
673 settings1->SOURCE_MAX_E=settings1->EXPONENT;
677 double src_min = TMath::Power(10,settings1->SOURCE_MIN_E-9);
678 double src_max = TMath::Power(10,settings1->SOURCE_MAX_E-9);
685 printf(
"Using Source Model %s\n", src_model->getName());
686 printf(
"Setting up weights: %g %g %g %g\n", bn1->min_time, bn1->max_time, src_min, src_max);
687 src_model->
setUpWeights(bn1->min_time, bn1->max_time , src_min, src_max);
693 rough1->SetRoughScale(settings1->ROUGHSIZE);
695 Screen *panel1 =
new Screen(0);
697 if(!src_model && spectra1->IsSpectrum())
699 cout <<
"Using an energy spectrum" << endl;
704 cout<<
"Source Model Bounds are 10^" << settings1->SOURCE_MIN_E <<
"eV to 10^" << settings1->SOURCE_MAX_E <<
"eV" << endl;
707 else if (!src_model && spectra1->IsMonoenergetic())
709 cout <<
"Using a monoenergetic exponent of: " << settings1->EXPONENT << endl;
712 std::cout <<
"----------------------" << std::endl;
719 time_t raw_start_time = time(NULL);
720 struct tm * start_time = localtime(&raw_start_time);
722 cout <<
"Date and time at start of run are: " << asctime (start_time) <<
"\n";
724 if (settings1->FORSECKEL)
725 SetupViewangles(sig1);
734 Tools::Zero(sum_frac, 3);
735 Tools::Zero(sum_frac_db, 3);
740 Tools::Zero(anita1->PHI_OFFSET, Anita::NLAYERS_MAX);
741 Tools::Zero(anita1->THETA_ZENITH, Anita::NLAYERS_MAX);
742 Tools::Zero(anita1->LAYER_VPOSITION, Anita::NLAYERS_MAX);
743 Tools::Zero(anita1->RRX, Anita::NLAYERS_MAX);
746 al_voltages_direct<<
"antenna #"<<
" "<<
"volts chan 1"<<
" "<<
"volts chan 2"<<
" "<<
"volts chan 3"<<
" "<<
"volts chan 4"<<
" "<<
"noise chan 1"<<
" "<<
"noise chan 2"<<
" "<<
"noise chan 3"<<
" "<<
"noise chan 4"<<
" "<<
"weight"<<endl;
751 double viewangle_temp=0;
752 double viewangle_deg=0;
753 double cosviewangle=0;
755 double nsigma_offaxis=0;
756 double theta_threshold=0;
757 double theta_threshold_deg=0;
759 double nsigma_em_threshold=0;
760 double nsigma_had_threshold=0;
761 double slopeyangle=0;
763 int beyondhorizon = 0;
766 double vmmhz1m_max=0;
767 double vmmhz_lowfreq=0.;
768 double bestcase_atten=0;
769 double vmmhz1m_fresneledonce=0;
770 double vmmhz1m_fresneledtwice=0;
771 double vmmhz[Anita::NFREQ];
774 double vmmhz_em[Anita::NFREQ];
775 double vmmhz_temp=1.;
776 double vmmhz_min_thatpasses=1000;
779 double deltheta_em[Anita::NFREQ], deltheta_had[Anita::NFREQ];
780 double deltheta_em_max, deltheta_had_max;
781 double deltheta_em_mid2, deltheta_had_mid2;
784 double emfrac, hadfrac, sumfrac;
785 int n_interactions=1;
786 double emfrac_db, hadfrac_db;
792 double d12=interaction1->
d1;
793 double d22=interaction1->
d2;
795 double logchord2=interaction1->
logchord;
800 double r_exit2bn2=interaction1->
r_exit2bn;
808 double volts_rx_max=0;
809 double volts_rx_ave=0;
810 double volts_rx_sum=0;
812 double volts_rx_max_highband;
813 double volts_rx_max_lowband;
814 double volts_rx_rfcm_lab_e_all[48][512];
815 double volts_rx_rfcm_lab_h_all[48][512];
818 static double e_component[Anita::NANTENNAS_MAX]={0};
819 static double h_component[Anita::NANTENNAS_MAX]={0};
820 static double n_component[Anita::NANTENNAS_MAX]={0};
822 static double e_component_kvector[Anita::NANTENNAS_MAX]={0};
823 static double h_component_kvector[Anita::NANTENNAS_MAX]={0};
824 static double n_component_kvector[Anita::NANTENNAS_MAX]={0};
826 Vector n_eplane = const_z;
827 Vector n_hplane = -const_y;
828 Vector n_normal = const_x;
833 double hitangle_e, hitangle_h;
834 double hitangle_e_all[Anita::NANTENNAS_MAX];
835 double hitangle_h_all[Anita::NANTENNAS_MAX];
838 double len_int_kgm2=0;
839 double eventsfound_db=0;
840 double eventsfound_nfb=0;
846 double dist_int_bn_2d=0;
847 double dist_int_bn_2d_chord=0;
855 double nuexitlength=0;
857 double nuentrancelength=0;
859 double icethickness=0;
860 double theta_pol_measured;
865 double nutauweightsum=0;
866 double tauweightsum=0;
867 double nutauweight=0;
869 int tauweighttrigger=0;
870 double sample_x=0, sample_y = 0;
882 Vector n_nutraject_ontheground;
890 int l2trig[
Anita::NPOL][Anita::NTRIGGERLAYERS_MAX];
892 int l1trig[
Anita::NPOL][Anita::NTRIGGERLAYERS_MAX];
899 int nchannels_triggered = 0;
900 int nchannels_perrx_triggered[48];
906 Tools::Zero(count1->npass, 2);
911 Tools::Zero(eventsfound_binned, NBINS);
912 Tools::Zero(eventsfound_binned_e, NBINS);
913 Tools::Zero(eventsfound_binned_mu, NBINS);
914 Tools::Zero(eventsfound_binned_tau, NBINS);
915 Tools::Zero(eventsfound_nfb_binned, NBINS);
920 TH2F *ref_int_coord=
new TH2F(
"ref_int_coord",
"", 600, -3000, 3000, 500, -2500, 2500);
921 ref_int_coord->SetMarkerSize(0.7);;
922 ref_int_coord->SetMarkerStyle(30);
923 ref_int_coord->SetMarkerColor(kBlue);
925 TH2F *dir_int_coord=
new TH2F(
"dir_int_coord",
"", 600, -3000, 3000, 500, -2500, 2500);
926 dir_int_coord->SetMarkerSize(0.7);
927 dir_int_coord->SetMarkerStyle(30);
929 TH1F *prob_eachphi_bn=
new TH1F(
"prob_eachphi_bn",
"prob_eachphi_bn", 100, 0., 6.3);
930 TH1F *prob_eachilon_bn=
new TH1F(
"prob_eachilon_bn",
"prob_eachilon_bn", 180, 0., 180.);
931 TH2F *h6=
new TH2F(
"theta_vs_hitangle_h",
"theta_vs_hitangle_h", 100, -3.14, 3.14, 100, -1.1, 1.1);
932 TH1F *h10=
new TH1F(
"hitangle_e",
"hitangle_e", 20, -1.6, 1.6);
933 TH1F *hy=
new TH1F(
"hy",
"hy", 100, 0., 1.);
934 TH1F *fraction_sec_muons =
new TH1F(
"fraction_sec_muons",
"fraction_sec_muons", 100, 0., 1.);
935 TH1F *fraction_sec_taus=
new TH1F(
"fraction_sec_taus",
"fraction_sec_taus", 100, 0., 1.);
936 TH1F *n_sec_muons=
new TH1F(
"n_sec_muons",
"n_sec_muons", 100, 0., 10.);
937 TH1F *n_sec_taus=
new TH1F(
"n_sec_taus",
"n_sec_taus", 100, 0., 10.);
938 TH1F *sampleweights=
new TH1F(
"sampleweights",
"sampleweights", 100, -5., 0.);
942 stemp=string(outputdir.Data())+
"/icefinal"+run_num+
".root";
943 TFile *hfile =
new TFile(stemp.c_str(),
"RECREATE",
"ice");
944 TTree *tree2 =
new TTree(
"h2000",
"h2000");
946 tree2->Branch(
"inu", &inu,
"inu/I");
947 tree2->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
948 tree2->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
949 tree2->Branch(
"scalefactor_distance", &scalefactor_distance,
"scalefactor_distance/D");
950 tree2->Branch(
"scalefactor_attenuation", &scalefactor_attenuation,
"scalefactor_attenuation/D");
952 TTree *tree3 =
new TTree(
"h3000",
"h3000");
953 tree3->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
954 tree3->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
955 tree3->Branch(
"theta_threshold_deg", &theta_threshold_deg,
"theta_threshold_deg/D");
956 tree3->Branch(
"nsigma_em_threshold", &nsigma_em_threshold,
"nsigma_em_threshold/D");
957 tree3->Branch(
"nsigma_had_threshold", &nsigma_had_threshold,
"nsigma_had_threshold/D");
958 tree3->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
959 tree3->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
960 tree3->Branch(
"vmmhz_max", &vmmhz_max,
"vmmhz_max/D");
961 tree3->Branch(
"vmmhz_min", &vmmhz_min,
"vmmhz_min/D");
962 tree3->Branch(
"dviewangle_deg", &dviewangle_deg,
"dviewangle_deg/D");
963 tree3->Branch(
"viewangle_deg", &viewangle_deg,
"viewangle_deg/D");
964 tree3->Branch(
"changle_deg", &changle_deg,
"changle_deg/D");
965 tree3->Branch(
"cosviewangle", &cosviewangle,
"cosviewangle/D");
966 tree3->Branch(
"emfrac", &emfrac,
"emfrac/D");
967 tree3->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
969 TTree *tree5 =
new TTree(
"h5000",
"h5000");
970 tree5->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
971 tree5->Branch(
"inu", &inu,
"inu/I");
972 tree5->Branch(
"nuexitlength", &nuexitlength,
"nuexitlength/D");
973 tree5->Branch(
"nuexitice", &nuexitice,
"nuexitice");
974 tree5->Branch(
"vmmhz_max", &vmmhz_max,
"vmmhz_max");
975 tree5->Branch(
"maxtaper", &maxtaper,
"maxtaper");
976 tree5->Branch(
"inu", &inu,
"inu/I");
977 tree5->Branch(
"whichray", &whichray,
"whichray/I");
978 tree5->Branch(
"pnu", &pnu,
"pnu/D");
979 tree5->Branch(
"costhetanu", &costhetanu,
"costhetanu/D");
980 tree5->Branch(
"viewangle", &viewangle,
"viewangle/D");
981 tree5->Branch(
"offaxis", &offaxis,
"offaxis/D");
982 tree5->Branch(
"nsigma_offaxis", &nsigma_offaxis,
"nsigma_offaxis/D");
983 tree5->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
984 tree5->Branch(
"emfrac", &emfrac,
"emfrac/D");
985 tree5->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
986 tree5->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
987 tree5->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
988 tree5->Branch(
"weight1", &weight1,
"weight1/D");
989 tree5->Branch(
"nearthlayers", &nearthlayers,
"nearthlayers/D");
990 tree5->Branch(
"logchord", &logchord2,
"interaction1->logchord/D");
991 tree5->Branch(
"diff_3tries", &diff_3tries,
"diff_3tries/D");
992 tree5->Branch(
"fresnel2", &fresnel2,
"fresnel2/D");
993 tree5->Branch(
"costheta_inc", &costheta_inc,
"costheta_inc/D");
994 tree5->Branch(
"costheta_exit", &costheta_exit,
"costheta_exit/D");
995 tree5->Branch(
"deltheta_em", &deltheta_em[0],
"deltheta_em/D");
996 tree5->Branch(
"deltheta_had", &deltheta_had[0],
"deltheta_had/F");
997 tree5->Branch(
"r_fromballoon", &interaction1->
r_fromballoon[0],
"r_fromballoon/F");
998 tree5->Branch(
"theta_in", &theta_in,
"theta_in/D");
999 tree5->Branch(
"lat_in", &lat_in,
"lat_in/D");
1002 TTree *tree6 =
new TTree(
"h6000",
"h6000");
1003 tree6->Branch(
"volts_rx_0", &volts_rx_0,
"volts_rx_0/D");
1004 tree6->Branch(
"volts_rx_1", &volts_rx_1,
"volts_rx_1/D");
1005 tree6->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
1006 tree6->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
1007 tree6->Branch(
"theta_in", &theta_in,
"theta_in/D");
1008 tree6->Branch(
"chord_kgm2_bestcase", &chord_kgm2_bestcase2,
"chord_kgm2_bestcase/D");
1009 tree6->Branch(
"chord_kgm2_ice", &interaction1->
chord_kgm2_ice,
"chord_kgm2_ice/D");
1010 tree6->Branch(
"costheta_nutraject", &interaction1->
costheta_nutraject,
"costheta_nutraject/D");
1011 tree6->Branch(
"weight1", &weight1,
"weight1/D");
1012 tree6->Branch(
"weight_bestcase", &weight_bestcase2,
"weight_bestcase/D");
1013 tree6->Branch(
"whichray", &whichray,
"whichray/I");
1014 tree6->Branch(
"mybeta", &mybeta,
"mybeta/D");
1015 tree6->Branch(
"longitude", &longitude_this,
"longitude/D");
1018 TTree *tree6b =
new TTree(
"h6001",
"h6001");
1019 tree6b->Branch(
"bwslice_vnoise", bwslice_vnoise_thislayer,
"bwslice_vnoise_thislayer[4]/D");
1021 TTree *tree7 =
new TTree(
"h7000",
"h7000");
1022 tree7->Branch(
"emfrac", &emfrac,
"emfrac/D");
1023 tree7->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1024 tree7->Branch(
"current", &interaction1->
currentint,
"currentint/I");
1025 tree7->Branch(
"nuflavor", &interaction1->
nuflavorint,
"nuflavorint/I");
1026 tree7->Branch(
"nunubar", &thisNuType,
"nunubar/I");
1027 tree7->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
1028 tree7->Branch(
"slopeyangle", &slopeyangle,
"slopeyangle/D");
1030 TTree *jaimetree=
new TTree(
"jaimetree",
"jaimetree");
1031 jaimetree->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
1032 jaimetree->Branch(
"emfrac", &emfrac,
"emfrac/D");
1033 jaimetree->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1034 jaimetree->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
1035 jaimetree->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
1036 jaimetree->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
1037 jaimetree->Branch(
"vmmhz1m_visible", &vmmhz1m_visible,
"vmmhz1m_visible/D");
1039 TTree *viewangletree=
new TTree(
"viewangletree",
"viewangletree");
1040 viewangletree->Branch(
"dviewangle_deg", &dviewangle_deg,
"dviewangle_deg/D");
1041 viewangletree->Branch(
"emfrac", &emfrac,
"emfrac/D");
1042 viewangletree->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1043 viewangletree->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
1044 viewangletree->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
1045 viewangletree->Branch(
"theta_threshold_deg", &theta_threshold_deg,
"theta_threshold_deg/D");
1046 viewangletree->Branch(
"dnutries", &interaction1->
dnutries,
"dnutries/D");
1047 viewangletree->Branch(
"viewangle", &viewangle,
"viewangle/D");
1048 viewangletree->Branch(
"chord", &interaction1->
chord,
"dnutries/D");
1050 TTree *neutrino_positiontree=
new TTree(
"neutrino_positiontree",
"neutrino_positiontree");
1051 neutrino_positiontree->Branch(
"nnu", &interaction1->
nnu,
"nnu[3]/D");
1052 neutrino_positiontree->Branch(
"dtryingdirection", &interaction1->
dtryingdirection,
"dtryingdirection/D");
1053 neutrino_positiontree->Branch(
"bn1->dtryingposition", &bn1->dtryingposition,
"bn1->dtryingposition/D");
1056 TTree *nupathtree=
new TTree(
"nupathtree",
"nupathtree");
1057 nupathtree->Branch(
"total_kgm2", &interaction1->
total_kgm2,
"total_kgm2/D");
1058 nupathtree->Branch(
"chord", &interaction1->
chord,
"chord/D");
1059 nupathtree->Branch(
"crust_entered", &interaction1->crust_entered,
"crust_entered/I");
1060 nupathtree->Branch(
"mantle_entered", &interaction1->mantle_entered,
"mantle_entered/I");
1061 nupathtree->Branch(
"core_entered", &interaction1->core_entered,
"core_entered/I");
1062 nupathtree->Branch(
"mybeta", &mybeta,
"mybeta/D");
1063 nupathtree->Branch(
"costheta_nutraject", &interaction1->
costheta_nutraject,
"costheta_nutraject/D");
1065 TTree *finaltree =
new TTree(
"passing_events",
"passing_events");
1066 finaltree->Branch(
"inu", &inu,
"inu/I");
1067 finaltree->Branch(
"sample_x", &sample_x,
"sample_x/D");
1068 finaltree->Branch(
"sample_y", &sample_y,
"sample_y/D");
1069 finaltree->Branch(
"vmmhz_min", &vmmhz_min,
"vmmhz_min/D");
1070 finaltree->Branch(
"vmmhz_max", &vmmhz_max,
"vmmhz_max/D");
1071 finaltree->Branch(
"thresholdsAnt", &thresholdsAnt,
"thresholdsAnt[48][2][5]/D");
1072 finaltree->Branch(
"thresholdsAntPass", &thresholdsAntPass,
"thresholdsAntPass[48][2][5]/D");
1073 finaltree->Branch(
"deadTime", &anita1->deadTime,
"deadTime/D");
1074 finaltree->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
1075 finaltree->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
1076 finaltree->Branch(
"horizcoord_bn", &bn1->horizcoord_bn,
"horizcoord_bn/D");
1077 finaltree->Branch(
"vertcoord_bn", &bn1->vertcoord_bn,
"vertcoord_bn/D");
1078 finaltree->Branch(
"r_bn", &r_bn_array,
"r_bn_array[3]/D");
1079 finaltree->Branch(
"n_bn", &n_bn_array,
"n_bn_array[3]/D");
1080 finaltree->Branch(
"longitude_bn", &longitude_this,
"longitude_bn/D");
1081 finaltree->Branch(
"heading_bn", &heading_this,
"heading_bn/D");
1082 finaltree->Branch(
"gps_offset", &gps_offset,
"gps_offset/D");
1084 finaltree->Branch(
"weight1", &weight1,
"weight1/D");
1086 finaltree->Branch(
"weight", &weight,
"weight/D");
1087 finaltree->Branch(
"logweight", &logweight,
"logweight/D");
1088 finaltree->Branch(
"posnu", &posnu_array,
"posnu_array[3]/D");
1089 finaltree->Branch(
"costheta_nutraject", &costheta_nutraject2,
"costheta_nutraject/D");
1090 finaltree->Branch(
"chord_kgm2_ice", &chord_kgm2_ice2,
"chord_kgm2_ice/D");
1091 finaltree->Branch(
"phi_nutraject", &phi_nutraject2,
"phi_nutraject/D");
1092 finaltree->Branch(
"altitude_int", &altitude_int2,
"altitude_int/D");
1093 finaltree->Branch(
"nnu", &nnu_array,
"nnu_array[3]/D");
1094 finaltree->Branch(
"n_exit2bn", &n_exit2bn_array,
"n_exit2bn_array[5][3]/D");
1095 finaltree->Branch(
"n_exit_phi", &n_exit_phi,
"n_exit_phi/D");
1096 finaltree->Branch(
"rfexit", &rfexit_array,
"rfexit_array[5][3]/D");
1098 finaltree->Branch(
"pnu", &pnu,
"pnu/D");
1099 finaltree->Branch(
"elast_y", &elast_y,
"elast_y/D");
1100 finaltree->Branch(
"emfrac", &emfrac,
"emfrac/D");
1101 finaltree->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1102 finaltree->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
1103 finaltree->Branch(
"nuflavor", &nuflavorint2,
"nuflavorint/I");
1104 finaltree->Branch(
"current", ¤tint2,
"currentint/I");
1105 finaltree->Branch(
"logchord", &logchord2,
"logchord/D");
1106 finaltree->Branch(
"nuexitice", &nuexitice,
"nuexitice/D");
1107 finaltree->Branch(
"weight_bestcase", &weight_bestcase2,
"weight_bestcase/D");
1108 finaltree->Branch(
"chord_kgm2_bestcase", &chord_kgm2_bestcase2,
"chord_kgm2_bestcase/D");
1109 finaltree->Branch(
"dtryingdirection", &dtryingdirection2,
"dtryingdirection/D");
1110 finaltree->Branch(
"l3trig", &l3trig,
"l3trig[2]/I");
1111 finaltree->Branch(
"l2trig", &l2trig,
"l2trig[2][3]/I");
1112 finaltree->Branch(
"l1trig", &l1trig,
"l1trig[2][3]/I");
1113 finaltree->Branch(
"phiTrigMask", &anita1->phiTrigMask,
"phiTrigMask/s");
1114 finaltree->Branch(
"phiTrigMaskH", &anita1->phiTrigMaskH,
"phiTrigMaskH/s");
1115 finaltree->Branch(
"l1TrigMask", &anita1->l1TrigMask,
"l1TrigMask/s");
1116 finaltree->Branch(
"l1TrigMaskH", &anita1->l1TrigMaskH,
"l1TrigMaskH/s");
1117 finaltree->Branch(
"max_antenna0", &max_antenna0,
"max_antenna0/I");
1118 finaltree->Branch(
"max_antenna1", &max_antenna1,
"max_antenna1/I");
1119 finaltree->Branch(
"max_antenna2", &max_antenna2,
"max_antenna2/I");
1121 finaltree->Branch(
"viewangle", &viewangle,
"viewangle/D");
1122 finaltree->Branch(
"offaxis", &offaxis,
"offaxis/D");
1123 finaltree->Branch(
"rx0_signal_eachband", &rx0_signal_eachband,
"rx0_signal_eachband[2][5]/D");
1124 finaltree->Branch(
"rx0_threshold_eachband", &rx0_threshold_eachband,
"rx0_threshold_eachband[2][5]/D");
1125 finaltree->Branch(
"rx0_noise_eachband", &rx0_noise_eachband,
"rx0_noise_eachband[2][5]/D");
1126 finaltree->Branch(
"rx0_passes_eachband", &rx0_passes_eachband,
"rx0_passes_eachband[2][5]/I");
1127 finaltree->Branch(
"e_component", &e_component,
"e_component[48]/D");
1128 finaltree->Branch(
"h_component", &h_component,
"h_component[48]/D");
1129 finaltree->Branch(
"dist_int_bn_2d", &dist_int_bn_2d,
"dist_int_bn_2d/D");
1130 finaltree->Branch(
"d1", &d12,
"d1/D");
1132 finaltree->Branch(
"cosalpha", &cosalpha,
"cosalpha/D");
1133 finaltree->Branch(
"mytheta", &mytheta,
"mytheta/D");
1134 finaltree->Branch(
"cosbeta0", &cosbeta0,
"cosbeta0/D");
1135 finaltree->Branch(
"mybeta", &mybeta,
"mybeta/D");
1136 finaltree->Branch(
"d1", &d12,
"d1/D");
1137 finaltree->Branch(
"d2", &d22,
"d2/D");
1140 finaltree->Branch(
"fresnel1", &fresnel1,
"fresnel1/D");
1141 finaltree->Branch(
"fresnel2", &fresnel2,
"fresnel2/D");
1142 finaltree->Branch(
"mag1", &mag1,
"mag1/D");
1143 finaltree->Branch(
"mag2", &mag2,
"mag2/D");
1144 finaltree->Branch(
"t_coeff_pokey", &t_coeff_pokey,
"t_coeff_pokey/D");
1145 finaltree->Branch(
"t_coeff_slappy", &t_coeff_slappy,
"t_coeff_slappy/D");
1146 finaltree->Branch(
"exponent", &settings1->EXPONENT,
"EXPONENT/D");
1148 finaltree->Branch(
"hitangle_e_all", &hitangle_e_all,
"hitangle_e_all[48]/D");
1149 finaltree->Branch(
"hitangle_h_all", &hitangle_h_all,
"hitangle_h_all[48]/D");
1151 finaltree->Branch(
"e_comp_max1", &e_comp_max1,
"e_comp_max1/D");
1152 finaltree->Branch(
"h_comp_max1", &h_comp_max1,
"h_comp_max1/D");
1153 finaltree->Branch(
"e_comp_max2", &e_comp_max2,
"e_comp_max2/D");
1154 finaltree->Branch(
"h_comp_max2", &h_comp_max2,
"h_comp_max2/D");
1155 finaltree->Branch(
"e_comp_max3", &e_comp_max3,
"e_comp_max3/D");
1156 finaltree->Branch(
"h_comp_max3", &h_comp_max3,
"h_comp_max3/D");
1157 finaltree->Branch(
"max_antenna_volts0", &max_antenna_volts0,
"max_antenna_volts0/D");
1158 finaltree->Branch(
"max_antenna_volts0_em", &max_antenna_volts0_em,
"max_antenna_volts0_em/D");
1159 finaltree->Branch(
"max_antenna_volts1", &max_antenna_volts1,
"max_antenna_volts1/D");
1160 finaltree->Branch(
"max_antenna_volts2", &max_antenna_volts2,
"max_antenna_volts2/D");
1161 finaltree->Branch(
"triggers", &nchannels_perrx_triggered,
"nchannels_perrx_triggered[48]/I");
1162 finaltree->Branch(
"nchannels_triggered", &nchannels_triggered,
"nchannels_triggered/I");
1163 finaltree->Branch(
"volts_rx_max", &volts_rx_max,
"volts_rx_max/D");
1164 finaltree->Branch(
"volts_rx_ave", &volts_rx_ave,
"volts_rx_ave/D");
1165 finaltree->Branch(
"volts_rx_sum", &volts_rx_sum,
"volts_rx_sum/D");
1167 finaltree->Branch(
"volts_rx_max_highband", &volts_rx_max_highband,
"volts_rx_max_highband/D");
1168 finaltree->Branch(
"volts_rx_max_lowband", &volts_rx_max_lowband,
"volts_rx_max_lowband/D");
1169 finaltree->Branch(
"theta_pol_measured", &theta_pol_measured,
"theta_pol_measured/D");
1170 finaltree->Branch(
"theta_rf_atbn", &theta_rf_atbn,
"theta_rf_atbn/D");
1171 finaltree->Branch(
"theta_rf_atbn_measured", &theta_rf_atbn_measured,
"theta_rf_atbn_measured/D");
1172 finaltree->Branch(
"theta_nnu_atbn",&theta_nnu_atbn,
"theta_nnu_atbn/D");
1173 finaltree->Branch(
"theta_nnu_rf_diff_atbn",&theta_nnu_rf_diff_atbn,
"theta_nnu_rf_diff_atbn/D");
1174 finaltree->Branch(
"phi_nnu_atbn",&phi_nnu_atbn,
"phi_nnu_atbn/D");
1175 finaltree->Branch(
"phi_rf_atbn",&phi_rf_atbn,
"phi_rf_atbn/D");
1176 finaltree->Branch(
"phi_nnu_rf_diff_atbn",&phi_nnu_rf_diff_atbn,
"phi_nnu_rf_diff_atbn/D");
1177 finaltree->Branch(
"voltage", &voltagearray,
"voltagearray[48]/D");
1178 finaltree->Branch(
"nlayers", &settings1->NLAYERS,
"settings1->NLAYERS/I");
1180 finaltree->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
1181 finaltree->Branch(
"vmmhz_lowfreq", &vmmhz_lowfreq,
"vmmhz_lowfreq/D");
1183 finaltree->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
1184 finaltree->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
1185 finaltree->Branch(
"r_enterice", &r_enterice_array,
"r_enterice_array[3]/D");
1186 finaltree->Branch(
"n_exit2bn_db", &n_exit2bn_db_array,
"n_exit2bn_db_array[5][3]/D");
1188 finaltree->Branch(
"rfexit_db", &rfexit_db_array,
"rfexit_db_array[5][3]/D");
1189 finaltree->Branch(
"r_in", &r_in_array,
"r_in_array[3]/D");
1190 finaltree->Branch(
"nsurf_rfexit", &nsurf_rfexit_array,
"nsurf_rfexit_array[3]/D");
1191 finaltree->Branch(
"nsurf_rfexit_db", &nsurf_rfexit_db_array,
"nsurf_rfexit_db_array[3]/D");
1192 finaltree->Branch(
"r_fromballoon", &r_fromballoon2,
"r_fromballoon/D");
1193 finaltree->Branch(
"r_fromballoon_db", &interaction1->
r_fromballoon_db,
"r_fromballoon_db/D");
1195 finaltree->Branch(
"nuexitlength", &nuexitlength,
"nuexitlength/D");
1196 finaltree->Branch(
"nuentrancelength", &nuentrancelength,
"nuentrancelength/D");
1197 finaltree->Branch(
"taulength", &taulength,
"taulength/D");
1198 finaltree->Branch(
"icethickness", &icethickness,
"icethickness/D");
1199 finaltree->Branch(
"nrf_iceside", &nrf_iceside_array,
"nrf_iceside_array[5][3]/D");
1200 finaltree->Branch(
"nrf_iceside_db", &nrf_iceside_db_array,
"nrf_iceside_db_array[5][3]/D");
1201 finaltree->Branch(
"ant_normal0", &ant_max_normal0_array,
"ant_max_normal0_array[3]/D");
1202 finaltree->Branch(
"ant_normal1", &ant_max_normal1_array,
"ant_max_normal1_array[3]/D");
1203 finaltree->Branch(
"ant_normal2", &ant_max_normal2_array,
"ant_max_normal2_array[3]/D");
1204 finaltree->Branch(
"vmmhz1m_visible", &vmmhz1m_visible,
"vmmhz1m_visible/D");
1205 finaltree->Branch(
"freq_bins", &freq_bins,
"freq_bins/I");
1206 finaltree->Branch(
"vmmhz", &vmmhz,
"vmmhz[freq_bins]/D");
1208 finaltree->Branch(
"dist_int_bn_2d_chord", &dist_int_bn_2d_chord,
"dist_int_bn_2d_chord/D");
1210 finaltree->Branch(
"dviewangle_deg", &dviewangle_deg,
"dviewangle_deg/D");
1211 finaltree->Branch(
"theta_threshold_deg", &theta_threshold_deg,
"theta_threshold_deg/D");
1212 finaltree->Branch(
"total_kgm2", &interaction1->
total_kgm2,
"total_kgm2/D");
1213 finaltree->Branch(
"chord", &interaction1->
chord,
"chord/D");
1214 finaltree->Branch(
"crust_entered", &interaction1->crust_entered,
"crust_entered/I");
1215 finaltree->Branch(
"mantle_entered", &interaction1->mantle_entered,
"mantle_entered/I");
1216 finaltree->Branch(
"core_entered", &interaction1->core_entered,
"core_entered/I");
1217 finaltree->Branch(
"n_pol", &n_pol_array,
"n_pol_array[3]/D");
1218 finaltree->Branch(
"vmmhz_min_thatpasses", &vmmhz_min_thatpasses,
"vmmhz_min_thatpasses/D");
1220 finaltree->Branch(
"pieceofkm2sr", &pieceofkm2sr,
"pieceofkm2sr/D");
1222 finaltree->Branch(
"r_exit2bn", &r_exit2bn2,
"r_exit2bn/D");
1223 finaltree->Branch(
"r_exit2bn_measured", &r_exit2bn_measured2,
"r_exit2bn_measured/D");
1224 finaltree->Branch(
"scalefactor_attenuation", &scalefactor_attenuation,
"scalefactor_attenuation/D");
1225 finaltree->Branch(
"anita1->PHI_OFFSET", &anita1->PHI_OFFSET,
"anita1->PHI_OFFSET/D");
1226 finaltree->Branch(
"igps", &bn1->igps,
"igyps/I");
1227 finaltree->Branch(
"volts_rx_rfcm_lab_e_all", &volts_rx_rfcm_lab_e_all,
"volts_rx_rfcm_lab_e_all[48][512]/D");
1228 finaltree->Branch(
"volts_rx_rfcm_lab_h_all", &volts_rx_rfcm_lab_h_all,
"volts_rx_rfcm_lab_h_all[48][512]/D");
1229 finaltree->Branch(
"ptaui", &ptaui,
"ptaui/D");
1230 finaltree->Branch(
"ptauf", &ptauf,
"ptauf/D");
1231 finaltree->Branch(
"sourceLon", &sourceLon,
"sourceLon/D");
1232 finaltree->Branch(
"sourceLat", &sourceLat,
"sourceLat/D");
1233 finaltree->Branch(
"sourceAlt", &sourceAlt,
"sourceAlt/D");
1234 finaltree->Branch(
"sourceMag", &sourceMag,
"sourceMag/D");
1238 TTree *mytaus_tree =
new TTree(
"mytaus",
"mytaus");
1239 mytaus_tree->Branch(
"taus", &TauPtr);
1248 double avgfreq_rfcm[Anita::NFREQ];
1249 double avgfreq_rfcm_lab[Anita::NFREQ];
1250 double freq[Anita::NFREQ];
1252 TTree *summarytree =
new TTree(
"summarytree",
"summarytree");
1253 summarytree->Branch(
"NNU", &NNU,
"NNU/I");
1254 summarytree->Branch(
"EXPONENT", &settings1->EXPONENT,
"EXPONENT/D");
1255 summarytree->Branch(
"eventsfound_beforetrigger", &eventsfound_beforetrigger,
"eventsfound_beforetrigger/D");
1256 summarytree->Branch(
"rms_rfcm_e", &rms_rfcm_e,
"rms_rfcm_e/D");
1257 summarytree->Branch(
"rms_rfcm_h", &rms_rfcm_h,
"rms_rfcm_h/D");
1258 summarytree->Branch(
"rms_lab_e", &rms_lab_e,
"rms_lab_e/D");
1259 summarytree->Branch(
"rms_lab_h", &rms_lab_h,
"rms_lab_h/D");
1260 summarytree->Branch(
"avgfreq_rfcm", &avgfreq_rfcm,
"avgfreq_rfcm[128]/D");
1261 summarytree->Branch(
"avgfreq_rfcm_lab", &avgfreq_rfcm_lab,
"avgfreq_rfcm_lab[128]/D");
1262 summarytree->Branch(
"freq", &freq,
"freq[128]/D");
1264 TTree *banana_tree =
new TTree(
"banana_tree",
"banana_tree");
1265 banana_tree->Branch(
"r_bn", &bn1->r_bn,
"r_bn[3]/D");
1267 TTree *ytree =
new TTree(
"ytree",
"ytree");
1268 ytree->Branch(
"elast_y", &elast_y,
"elast_y/D");
1280 TTree *icetree =
new TTree(
"icetree",
"icetree");
1281 icetree->Branch(
"icethck", &icethck,
"icethck/D");
1282 icetree->Branch(
"lon_ice", &lon_ice,
"lon_ice/D");
1283 icetree->Branch(
"lat_ice", &lat_ice,
"lat_ice/D");
1284 icetree->Branch(
"lon_water", &lon_water,
"lon_water/D");
1285 icetree->Branch(
"lat_water", &lat_water,
"lat_water/D");
1286 icetree->Branch(
"h20_depth", &h20_depth,
"h20_depth/D");
1288 TTree *groundtree =
new TTree(
"groundtree",
"groundtree");
1289 groundtree->Branch(
"elev", &elev,
"elev/D");
1290 groundtree->Branch(
"lon_ground", &lon_ground,
"lon_ground/D");
1291 groundtree->Branch(
"lat_ground", &lat_ground,
"lat_ground/D");
1295 TTree *tree11 =
new TTree(
"h11000",
"h11000");
1296 tree11->Branch(
"loctrig00", &loctrig[0][0],
"loctrig0/D");
1297 tree11->Branch(
"loctrig10", &loctrig[1][0],
"loctrig0/D");
1298 tree11->Branch(
"loctrig20", &loctrig[2][0],
"loctrig0/D");
1299 tree11->Branch(
"loctrig_nadironly0", &loctrig_nadironly[0],
"loctrig_nadironly0/D");
1300 tree11->Branch(
"loctrig01", &loctrig[0][1],
"loctrig1/D");
1301 tree11->Branch(
"loctrig11", &loctrig[1][1],
"loctrig1/D");
1302 tree11->Branch(
"loctrig21", &loctrig[2][1],
"loctrig1/D");
1303 tree11->Branch(
"loctrig_nadironly1", &loctrig_nadironly[1],
"loctrig0/D");
1305 TTree *tree16 =
new TTree(
"h16000",
"h16000");
1306 tree16->Branch(
"pnu", &pnu,
"pnu/D");
1307 tree16->Branch(
"ptau", &ptau,
"ptau/D");
1308 tree16->Branch(
"taulength", &taulength,
"taulength/D");
1309 tree16->Branch(
"weight1", &weight1,
"weight1/D");
1310 tree16->Branch(
"emfrac", &emfrac,
"emfrac/D");
1311 tree16->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1312 tree16->Branch(
"nuentrancelength", &nuentrancelength,
"nuentrancelength/D");
1316 TTree *tree18 =
new TTree(
"h18000",
"h18000");
1317 tree18->Branch(
"emfrac", &emfrac,
"emfrac/D");
1318 tree18->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1319 tree18->Branch(
"pdgcode", &pdgcode,
"pdgcode/I");
1322 TH1D *h1mybeta =
new TH1D(
"betaforall",
"betaforall(deg)", 180, -15, 15);
1323 TH1D *h1mytheta=
new TH1D(
"mytheta",
"mytheta(deg)", 180, -90, 90);
1324 TH1F *hundogaintoheight_e=
new TH1F(
"undogaintoheight_e",
"undogaintoheight_e", 100, 0., 1.);
1325 TH1F *hundogaintoheight_h=
new TH1F(
"undogaintoheight_h",
"undogaintoheight_h", 100, 0., 1.);
1326 TH1F *rec_diff=
new TH1F(
"rec_diff",
"rec_diff", 100, -1., 1.);
1327 TH1F *recsum_diff=
new TH1F(
"recsum_diff",
"recsum_diff", 100, -1., 1.);
1328 TH1F *rec_diff0=
new TH1F(
"rec_diff0",
"rec_diff0", 100, -1., 1.);
1329 TH1F *rec_diff1=
new TH1F(
"rec_diff1",
"rec_diff1", 100, -1., 1.);
1330 TH1F *rec_diff2=
new TH1F(
"rec_diff2",
"rec_diff2", 100, -1., 1.);
1331 TH1F *rec_diff3=
new TH1F(
"rec_diff3",
"rec_diff3", 100, -1., 1.);
1333 TTree *vmmhz_tree =
new TTree(
"vmmhz_tree",
"vmmhz_tree");
1334 vmmhz_tree->Branch(
"freq_bins", &freq_bins,
"freq_bins/I");
1335 vmmhz_tree->Branch(
"vmmhz", &vmmhz,
"vmmhz[freq_bins]/D");
1337 TTree *tree1 =
new TTree(
"h1000",
"h1000");
1338 tree1->Branch(
"inu", &inu,
"inu/I");
1339 tree1->Branch(
"diffexit", &diffexit,
"diffexit/D");
1340 tree1->Branch(
"diffrefr", &diffrefr,
"diffrefr/D");
1341 tree1->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
1342 tree1->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
1343 tree1->Branch(
"costhetanu", &costhetanu,
"costhetanu/D");
1344 tree1->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
1345 tree1->Branch(
"volume_thishorizon", &volume_thishorizon,
"volume_thishorizon/D");
1346 tree1->Branch(
"realtime", &realtime_this,
"realtime/D");
1347 tree1->Branch(
"longitude", &longitude_this,
"longitude/D");
1348 tree1->Branch(
"latitude", &latitude_this,
"latitude/D");
1349 tree1->Branch(
"MAXHORIZON", &bn1->MAXHORIZON,
"MAXHORIZON/D");
1350 tree1->Branch(
"igps", &bn1->igps,
"igps/I");
1351 tree1->Branch(
"passes_thisevent", &passes_thisevent,
"passes_thisevent/I");
1352 tree1->Branch(
"igps", &bn1->igps,
"igps/I");
1353 tree1->Branch(
"weight", &weight,
"weight/D");
1354 tree1->Branch(
"r_exit2bn", &interaction1->
r_exit2bn,
"r_exit2bn/D");
1355 tree1->Branch(
"bn1->igps", &bn1->igps,
"bn1->igps/I");
1359 TTree *balloontree =
new TTree(
"balloon",
"balloon");
1360 balloontree->Branch(
"heading", &bn1->heading,
"heading/D");
1361 balloontree->Branch(
"pitch", &bn1->pitch,
"pitch/D");
1362 balloontree->Branch(
"roll", &bn1->roll,
"roll/D");
1363 balloontree->Branch(
"realTime_flightdata", &bn1->realTime_flightdata,
"realTime_flightdata/I");
1364 balloontree->Branch(
"latitude", &bn1->latitude,
"latitude/D");
1365 balloontree->Branch(
"longitude", &bn1->longitude,
"longitude/D");
1366 balloontree->Branch(
"altitude", &bn1->altitude,
"altitude/D");
1367 balloontree->Branch(
"horizcoord_bn", &bn1->horizcoord_bn,
"horizcoord_bn/D");
1368 balloontree->Branch(
"vertcoord_bn", &bn1->vertcoord_bn,
"vertcoord_bn/D");
1371 double undogaintoheight_e=0;
1372 double undogaintoheight_h=0;
1374 double undogaintoheight_e_array[4];
1375 double undogaintoheight_h_array[4];
1377 double nbins_array[4];
1379 double rec_efield=0;
1380 double true_efield=0;
1382 double rec_efield_array[4];
1383 double true_efield_array[4];
1389 IceModel *antarctica =
new IceModel(settings1->ICE_MODEL + settings1->NOFZ*10, settings1->CONSTANTICETHICKNESS * 1000 + settings1->CONSTANTCRUST * 100 + settings1->FIXEDELEVATION * 10 + 0, settings1->WEIGHTABSORPTION);
1390 cout <<
"area of the earth's surface covered by antarctic ice is " << antarctica->ice_area <<
"\n";
1392 ULong_t eventNumber;
1393 #ifdef ANITA_UTIL_EXISTS 1395 string outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaEventFile"+run_num+
".root";
1396 TFile *anitafileEvent =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1398 TTree *eventTree =
new TTree(
"eventTree",
"eventTree");
1399 eventTree->Branch(
"event", &realEvPtr );
1400 eventTree->Branch(
"run", &run_no,
"run/I" );
1401 eventTree->Branch(
"weight", &weight,
"weight/D");
1403 outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaHeadFile"+run_num+
".root";
1404 TFile *anitafileHead =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1406 TTree *headTree =
new TTree(
"headTree",
"headTree");
1407 headTree->Branch(
"header", &rawHeaderPtr );
1408 headTree->Branch(
"weight", &weight,
"weight/D");
1410 outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaGpsFile"+run_num+
".root";
1411 TFile *anitafileGps =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1413 TTree *adu5PatTree =
new TTree(
"adu5PatTree",
"adu5PatTree");
1414 adu5PatTree->Branch(
"pat", &Adu5PatPtr );
1415 adu5PatTree->Branch(
"eventNumber", &eventNumber,
"eventNumber/L");
1416 adu5PatTree->Branch(
"weight", &weight,
"weight/D" );
1418 #ifdef ANITA3_EVENTREADER 1421 AnitaVersion::set(settings1->ANITAVERSION);
1423 outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaTruthFile"+run_num+
".root";
1424 TFile *anitafileTruth =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1427 TString icemcgitversion( EnvironmentVariable::ICEMC_VERSION(outputdir));
1428 printf(
"ICEMC GIT Repository Version: %s\n", icemcgitversion.Data());
1429 unsigned int timenow = time(NULL);
1431 TTree *configAnitaTree =
new TTree(
"configIcemcTree",
"Config file and settings information");
1432 configAnitaTree->Branch(
"gitversion", &icemcgitversion );
1433 configAnitaTree->Branch(
"nnu", &NNU );
1434 configAnitaTree->Branch(
"startTime", &timenow );
1436 configAnitaTree->Fill();
1438 TTree *triggerSettingsTree =
new TTree(
"triggerSettingsTree",
"Trigger settings");
1439 triggerSettingsTree->Branch(
"dioderms", anita1->bwslice_dioderms_fullband_allchan,
"dioderms[2][48][7]/D" );
1440 triggerSettingsTree->Branch(
"diodemean", anita1->bwslice_diodemean_fullband_allchan,
"diodemean[2][48][7]/D");
1441 triggerSettingsTree->Fill();
1444 TTree *summaryAnitaTree =
new TTree(
"summaryAnitaTree",
"summaryAnitaTree");
1445 summaryAnitaTree->Branch(
"EXPONENT", &settings1->EXPONENT,
"EXPONENT/D" );
1446 summaryAnitaTree->Branch(
"SELECTION_MODE", &settings1->UNBIASED_SELECTION,
"SELECTION_MODE/I" );
1447 summaryAnitaTree->Branch(
"totalnu_per_source", &NNU_per_source);
1448 summaryAnitaTree->Branch(
"source_names", &source_names);
1449 summaryAnitaTree->Branch(
"weighted_nu_per_source", &eventsfound_per_source);
1450 summaryAnitaTree->Branch(
"weighted_prob_nu_per_source", &eventsfound_prob_per_source);
1451 summaryAnitaTree->Branch(
"total_nu", &NNU,
"total_nu/I" );
1452 summaryAnitaTree->Branch(
"total_nue", &count1->nnu_e,
"total_nue/I" );
1453 summaryAnitaTree->Branch(
"total_numu", &count1->nnu_mu,
"total_numu/I" );
1454 summaryAnitaTree->Branch(
"total_nutau", &count1->nnu_tau,
"total_nutau/I" );
1455 summaryAnitaTree->Branch(
"pass_nu", &count1->npass[0],
"pass_nu/I" );
1458 summaryAnitaTree->Branch(
"weighted_nu", &eventsfound,
"weighted_nu/D" );
1459 summaryAnitaTree->Branch(
"weighted_nue", &sum[0],
"weighted_nue/D" );
1460 summaryAnitaTree->Branch(
"weighted_numu", &sum[1],
"weighted_numu/D" );
1461 summaryAnitaTree->Branch(
"weighted_nutau", &sum[2],
"weighted_nutau/D" );
1463 summaryAnitaTree->Branch(
"weighted_prob_nu", &eventsfound_prob,
"weighted_prob_nu/D" );
1464 summaryAnitaTree->Branch(
"weighted_prob_nue", &sum_prob[0],
"weighted_prob_nue/D" );
1465 summaryAnitaTree->Branch(
"weighted_prob_numu", &sum_prob[1],
"weighted_prob_numu/D" );
1466 summaryAnitaTree->Branch(
"weighted_prob_nutau", &sum_prob[2],
"weighted_prob_nutau/D" );
1469 summaryAnitaTree->Branch(
"int_length", &len_int_kgm2,
"int_length/D");
1470 summaryAnitaTree->Branch(
"total_volume", &antarctica->volume,
"total_volume/D");
1471 summaryAnitaTree->Branch(
"rho_medium", &sig1->RHOMEDIUM,
"rho_medium/D");
1472 summaryAnitaTree->Branch(
"selection_box_half_length", &settings1->UNBIASED_PS_MAX_DISTANCE_KM,
"selection_box_half_length/D");
1475 summaryAnitaTree->Branch(
"effVol_nu", &km3sr,
"effVol_nu/D" );
1476 summaryAnitaTree->Branch(
"effArea_nu", &km2sr,
"effArea_nu/D");
1478 TTree *truthAnitaTree =
new TTree(
"truthAnitaTree",
"Truth Anita Tree");
1479 truthAnitaTree->Branch(
"truth", &truthEvPtr );
1482 TTree* truthAnitaNuTree = settings1->SAVE_TRUTH_NU_TREE ?
new TTree(
"truthAnitaNuTree",
"Truth ANITA Neutrino Tree (all nus)") : 0;
1483 if (settings1->SAVE_TRUTH_NU_TREE)
1485 truthAnitaNuTree->SetAutoFlush(10000);
1486 truthAnitaNuTree->Branch(
"truth_neutrino", &truthNuPtr );
1497 #ifdef ANITA3_EVENTCORRELATOR 1498 const int skyMapLat = 90;
1499 TMarker *astroObject =
new TMarker();
1503 if (bn1->WHICHPATH==3) {
1504 settings1->CONSTANTCRUST = 1;
1505 settings1->CONSTANTICETHICKNESS = 1;
1506 settings1->FIXEDELEVATION = 1;
1510 for (
int i=0;i<antarctica->nRows_ice;i++) {
1511 for (
int j=0;j<antarctica->nCols_ice;j++) {
1512 antarctica->IceENtoLonLat(j, i, lon_ice, lat_ice);
1513 icethck=antarctica->IceThickness(lon_ice, lat_ice);
1516 h20_depth=antarctica->h_water_depth.GetBinContent(j+1,i+1);
1517 if (settings1->HIST) icetree->Fill();
1520 for (
int i=0;i<antarctica->nRows_ground;i++) {
1521 for (
int j=0;j<antarctica->nCols_ground;j++) {
1522 antarctica->GroundENtoLonLat(j, i, lon_ground, lat_ground);
1523 elev=antarctica->SurfaceAboveGeoid(lon_ground, lat_ground);
1524 if (settings1->HIST) groundtree->Fill();
1529 anita1->GetBeamWidths(settings1);
1531 anita1->ReadGains();
1532 anita1->Set_gain_angle(settings1, sig1->NMEDIUM_RECEIVER);
1533 if(settings1->WHICH == 2 || settings1->WHICH == 6) anita1->SetDiffraction();
1535 TCanvas *cgains=
new TCanvas(
"cgains",
"cgains", 880, 800);
1536 TGraph *ggains=
new TGraph(anita1->NPOINTS_GAIN, anita1->frequency_forgain_measured, anita1->vvGaintoHeight);
1538 stemp=string(outputdir.Data())+
"/gains.eps";
1539 cgains->Print((TString)stemp);
1545 bn1->SetDefaultBalloonPosition(antarctica);
1547 anita1->SetNoise(settings1, bn1, antarctica);
1551 antarctica->GetMAXHORIZON(bn1);
1554 antarctica->CreateHorizons(settings1, bn1, bn1->theta_bn, bn1->phi_bn, bn1->altitude_bn, foutput);
1555 cout <<
"Done with CreateHorizons.\n";
1558 if (! src_model && spectra1->IsMonoenergetic() ){
1559 pnu=pow(10., settings1->EXPONENT);
1560 primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, xsecParam_nutype, Interaction::ktotal);
1561 cout <<
"pnu, sigma and len_int_kgm2 " << pnu <<
" " << sigma <<
" " << len_int_kgm2 <<
"\n";
1564 if (settings1->WRITEPOSFILE==1) {
1565 nu_out <<
"Neutrinos with energy " << pnu <<
"\n\n";
1566 nu_out <<
"nu #, position of nu interaction, depth of int., Direction of nu momentum, Position of balloon, nu flavour, current type, elasticity\n\n\n\n";
1569 if (bn1->WHICHPATH==3) {
1570 NNU = settings1->horizontal_banana_points*settings1->vertical_banana_points;
1573 settings1->CONSTANTCRUST=1;
1574 settings1->CONSTANTICETHICKNESS=1;
1575 settings1->FIXEDELEVATION=1;
1576 pnu=Interaction::pnu_banana;
1577 int_banana->surface_over_banana_nu=antarctica->Surface(int_banana->
nu_banana);
1581 cout <<
"reminder that I took out ChangeCoord.\n";
1584 anita1->GetPayload(settings1, bn1);
1586 if (settings1->TRIGGERSCHEME == 3 || settings1->TRIGGERSCHEME == 4 || settings1->TRIGGERSCHEME==5) {
1587 Vector plusz(0., 0., 1.);
1588 bn1->PickBalloonPosition(plusz, antarctica, settings1, anita1);
1589 anita1->calculate_all_offsets();
1590 double angle_theta=16.;
1591 double angle_phi=0.;
1592 Vector x =
Vector(cos(angle_theta * RADDEG) * cos((angle_phi+11.25) * RADDEG),
1593 cos(angle_theta * RADDEG) * sin((angle_phi+11.25) * RADDEG),
1594 sin(angle_theta * RADDEG));
1595 anita1->GetArrivalTimes(x,bn1,settings1);
1596 cout <<
"end of getarrivaltimes\n";
1600 if (settings1->SLAC)
1601 bn1->GetSlacPositions(anita1);
1603 switch (settings1->WHICHRAYS){
1605 settings1->MINRAY=0;
1606 settings1->MAXRAY=0;
1609 settings1->MINRAY=0;
1610 settings1->MAXRAY=1;
1613 settings1->MINRAY=1;
1614 settings1->MAXRAY=1;
1619 time_t raw_loop_start_time = time(NULL);
1620 cout<<
"Starting loop over events. Time required for setup is "<<(int)((raw_loop_start_time - raw_start_time)/60)<<
":"<< ((raw_loop_start_time - raw_start_time)%60)<<endl;
1623 TCanvas *ctest1 =
new TCanvas(
"ctest1",
"", 880, 800);
1626 if ( spectra1->IsSpectrum() ){
1627 spectra1->GetGEdNdEdAdt()->Draw(
"al");
1628 spectra1->GetGEdNdEdAdt()->GetHistogram()->SetMinimum(-0.2*spectra1->Getmaxflux());
1629 spectra1->GetSEdNdEdAdt()->SetLineColor(2);
1630 spectra1->GetSEdNdEdAdt()->Draw(
"l same");
1632 stemp=string(outputdir.Data())+
"/GetG_test1.pdf";
1633 ctest1->Print((TString)stemp);
1638 double average_altitude=0.;
1639 double average_rbn=0.;
1642 if ( spectra1->IsSpectrum() ){
1643 TCanvas *ctemp =
new TCanvas(
"ctemp");
1645 TFile *out =
new TFile(
"Temp.root",
"recreate");
1646 TGraph *g1 = spectra1->GetGEdNdEdAdt();
1648 ctemp->Print(
"Temp1.png");
1650 double *x = g1->GetX();
1652 for (
int i=0;i<n;i++) x2[i] = TMath::Power(10., x[i]);
1653 TGraph *g2 =
new TGraph(n, x2, g1->GetY());
1657 cout << g2->Integral() <<
" " << g2->Integral(1, n) << endl;
1658 ctemp->Print(
"Temp2.png");
1665 signal(SIGINT, interrupt_signal_handler);
1667 TString objName =
"noObject";
1670 if (settings1->WHICH==7){
1671 gps_offset=atan2(-0.7042,0.71)*DEGRAD;
1672 }
else if(settings1->WHICH==8){
1673 gps_offset=atan2(-0.7085,0.7056)*DEGRAD;
1674 }
else if (settings1->WHICH==9 || settings1->WHICH==10){
1676 }
else gps_offset=0;
1682 for (inu = startNu; inu < NNU; inu++) {
1685 if (inu % (NNU / 100) == 0)
1686 cout << inu <<
" neutrinos. " << (double(inu)/double(NNU)) * 100 <<
"% complete.\n";
1689 cout << inu <<
" neutrinos. " << (double(inu) / double(NNU)) * 100 <<
"% complete.\n";
1691 eventNumber=(ULong_t)(run_no)*NNU+inu;
1695 if (!eventListFile.eof()){
1696 eventListFile >> eventNumber;
1697 cout <<
"Looking at event number " << eventNumber << endl;
1703 eventNumber=(ULong_t)(run_no)*NNU+inu;
1706 setSeed(eventNumber);
1709 panel1->ResetParameters();
1712 std::string nunum = Form(
"%d",inu);
1714 double time_weight = 1;
1715 double src_time_weight = 1;
1717 for (whichray = settings1->MINRAY; whichray <= settings1->MAXRAY; whichray++) {
1718 anita1->passglobtrig[0]=0;
1719 anita1->passglobtrig[1]=0;
1721 unmasked_thisevent=1;
1722 vmmhz_min_thatpasses=1000;
1728 int which_source = -1;
1729 if (!src_model && spectra1->IsSpectrum() ){
1731 if(settings1->USEDARTBOARD) pnu=spectra1->GetNuEnergy();
1732 else pnu=spectra1->GetCDFEnergy();
1738 interaction1 =
new (interaction1)
Interaction(
"nu", primary1, settings1, whichray, count1);
1740 thisNuType = getRNG(RNG_INTERACTION)->Integer(2);
1742 ierr=primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, thisNuType, Interaction::ktotal);
1743 if (!src_model) interaction1->
currentint = primary1->GetCurrent(pnu,settings1,thisNuType);
1747 len_int=1.0/(sigma*sig1->RHOICE*(1./M_NUCL)*1000);
1753 sec1->secondbang=
false;
1757 bn1->dtryingposition=0;
1758 for (
int i=0; i<Anita::NFREQ;i++) {
1763 int got_a_good_position = 0;
1769 bn1->PickBalloonPosition(antarctica, settings1, inu, anita1, getRNG(RNG_BALLOON_POSITION)->Rndm());
1774 average_altitude+=bn1->altitude_bn/(double)NNU;
1775 average_rbn+=bn1->r_bn.Mag()/(double)NNU;
1777 realtime_this=bn1->realTime_flightdata;
1778 longitude_this=bn1->longitude;
1779 latitude_this=bn1->latitude;
1780 altitude_this=bn1->altitude;
1781 heading_this=bn1->heading;
1785 time_weight = src_model->getTimeWeight(realtime_this);
1787 which_source = time_weight > 0 ? src_model->
getDirectionAndEnergy(&force_dir, realtime_this, pnu, src_min, src_max) : -1;
1788 if(which_source >= 0) {
1789 got_a_good_position = 1;
1790 src_time_weight = src_model->getPerSourceTimeWeight(realtime_this, which_source);
1792 if (
int(NNU_per_source.size()) < which_source+1)
1794 NNU_per_source.resize(which_source+1);
1795 source_names.resize(which_source+1);
1796 eventsfound_per_source.resize(which_source+1);
1797 eventsfound_prob_per_source.resize(which_source+1);
1800 NNU_per_source[which_source]++;
1802 const Source * src = src_model->getSource(which_source);
1804 dec= src->getDec()* TMath::RadToDeg();
1805 objName = src->getName();
1806 if (!source_names[which_source].size())
1808 source_names[which_source] = objName.Data();
1812 ierr=primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, thisNuType, Interaction::ktotal);
1813 interaction1->
currentint = primary1->GetCurrent(pnu,settings1,thisNuType);
1814 len_int=1.0/(sigma*sig1->RHOICE*(1./M_NUCL)*1000);
1825 got_a_good_position = 1;
1827 if (settings1->HORIZON_OFFSET > -999)
1829 double horizon_angle = acos(bn1->r_bn_shadow.Mag()/bn1->r_bn.Mag());
1830 double alt = -horizon_angle + TMath::DegToRad() * settings1->HORIZON_OFFSET;
1832 double lat = bn1->latitude*TMath::DegToRad();
1833 double sin_lat = sin(lat);
1834 double cos_lat = cos(lat);
1835 double sindec = sin_lat *sin(alt) - cos_lat*cos(alt) *cos(Az);
1836 double cosdec = sqrt(1-sindec*sindec);
1837 double h = atan2(sin(Az),-sin_lat*cos(Az) + cos_lat*tan(alt));
1838 force_dir =
Vector(cos(h)*cosdec,sin(h)*cosdec,sindec);
1842 }
while(!got_a_good_position && settings1->SOURCE_SKIP_WHEN_NONE) ;
1844 if (settings1->HIST && !settings1->ONLYFINAL
1845 && prob_eachphi_bn->GetEntries() < settings1->HIST_MAX_ENTRIES) {
1846 prob_eachphi_bn->Fill(bn1->phi_bn);
1847 prob_eachilon_bn->Fill(bn1->r_bn.Lon());
1848 balloontree->Fill();
1851 if (bn1->WHICHPATH==3) {
1853 bn1->setObservationLocation(int_banana, inu, antarctica, settings1);
1857 if (!got_a_good_position)
1859 #ifdef ANITA3_EVENTREADER 1861 if (truthNuPtr) truthNuPtr->setNoNu(bn1->r_bn.GetX(), bn1->r_bn.GetY(), bn1->r_bn.GetZ(), realtime_this);
1879 int taumodes = settings1->taumodes;
1885 if (taumodes==1 && interaction1->nuflavor()==
"nutau" && interaction1->current()==
"cc")
1890 bn1->PickDownwardInteractionPoint(interaction1, anita1, settings1, antarctica, ray1, beyondhorizon, len_int_kgm2, (src_model || settings1->HORIZON_OFFSET > -999) ? &force_dir : 0);
1892 bool havent_set_frac =
true;
1893 bool havent_set_weights =
true;
1894 bool havent_set_dir =
true;
1896 #ifdef ANITA3_EVENTREADER 1898 if (settings1->UNBIASED_SELECTION > 0 || src_model || settings1->HORIZON_OFFSET)
1900 havent_set_dir =
false;
1901 if (truthNuPtr) truthNuPtr->setDir(interaction1->
nnu.GetX(), interaction1->
nnu.GetY(), interaction1->
nnu.GetZ());
1907 truthNuPtr->setPos(interaction1->posnu.GetX(), interaction1->posnu.GetY(), interaction1->posnu.GetZ(),
1908 bn1->r_bn.GetX(), bn1->r_bn.GetY(), bn1->r_bn.GetZ(), realtime_this);
1909 truthNuPtr->setNu(pnu,pdgcode);
1913 #ifdef ANITA3_EVENTREADER 1914 #define DO_SKIP { if (truthNuPtr) {truthNuPtr->setSkipped(true,havent_set_frac,havent_set_weights,havent_set_dir); truthAnitaNuTree->Fill();} continue; } 1916 #define DO_SKIP continue ; 1921 if (interaction1->noway)
1925 count1->noway[whichray]++;
1927 if (interaction1->wheredoesitleave_err)
1931 count1->wheredoesitleave_err[whichray]++;
1933 if (interaction1->neverseesice)
1937 count1->neverseesice[whichray]++;
1939 if (interaction1->wheredoesitenterice_err)
1943 count1->wheredoesitenterice_err[whichray]++;
1945 if (interaction1->toohigh)
1949 count1->toohigh[whichray]++;
1951 if (interaction1->toolow)
1955 count1->toolow[whichray]++;
1957 if (bn1->WHICHPATH==3)
1958 interaction1=int_banana;
1964 count1->iceinteraction[whichray]++;
1966 if (beyondhorizon) {
1969 count1->inhorizon[whichray]++;
1972 if (settings1->FIRN) {
1973 sig1->SetNDepth(antarctica->GetN(interaction1->
altitude_int));
1975 changle_deg=sig1->changle*DEGRAD;
1979 if (settings1->FORSECKEL==1)
1980 sig1->SetChangle(acos(1/sig1->NICE));
1983 horizcoord=interaction1->posnu[0]/1000;
1984 vertcoord=interaction1->posnu[1]/1000;
1986 ray1->GetSurfaceNormal(settings1, antarctica, interaction1->posnu, slopeyangle, 0);
1992 costheta_inc=ray1->n_exit2bn[0]*ray1->nsurf_rfexit;
1995 costheta_exit=cos(ray1->rfexit[0].Theta());
1998 if (!ray1->TraceRay(settings1, anita1, 1, sig1->N_DEPTH)) {
2006 ray1->GetRFExit(settings1, anita1, whichray, interaction1->posnu, interaction1->posnu_down, bn1->r_bn, bn1->r_boresights, 1, antarctica);
2008 ray1->GetSurfaceNormal(settings1, antarctica, interaction1->posnu, slopeyangle, 1);
2010 if (!ray1->TraceRay(settings1, anita1, 2, sig1->N_DEPTH)) {;
2017 ray1->GetRFExit(settings1, anita1, whichray, interaction1->posnu, interaction1->posnu_down, bn1->r_bn, bn1->r_boresights, 2, antarctica);
2019 ray1->GetSurfaceNormal(settings1, antarctica, interaction1->posnu, slopeyangle, 2);
2021 if (bn1->WHICHPATH==4)
2022 ray1->PrintAnglesofIncidence();
2025 count1->nraypointsup1[whichray]++;
2027 sec1->GetTauDecay(interaction1->nuflavor(), interaction1->current(), taudecay, emfrac_db, hadfrac_db);
2030 elast_y=primary1->
pickY(settings1, pnu, thisNuType, interaction1->
currentint);
2031 if (settings1->CONSTANTY==1) {
2036 else if (settings1->CONSTANTY==2)
2038 if (settings1->SHOWERTYPE == 0)
2051 if (bn1->WHICHPATH==3)
2053 if (bn1->WHICHPATH==4)
2055 if (settings1->FORSECKEL==1) {
2056 if (settings1->SHOWERTYPE==0)
2058 if (settings1->SHOWERTYPE==1)
2062 if (ytree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2067 if ( tautrigger == 1 ) {
2068 if ( ( !settings1->UNBIASED_SELECTION ) && ( !settings1->SLAC ) && !src_model ) {
2069 err = GetDirection(settings1, interaction1, ray1->nrf_iceside[4], deltheta_em_max, deltheta_had_max, emfrac, hadfrac, vmmhz1m_max*bestcase_atten, interaction1->
r_fromballoon[whichray], ray1, sig1, interaction1->posnu, anita1, bn1, interaction1->
nnu, costhetanu, theta_threshold);
2072 else if (src_model || settings1->HORIZON_OFFSET > -999)
2074 interaction1->
nnu = force_dir;
2076 costhetanu = cos(force_dir.Theta());
2077 theta_threshold = 1;
2080 else if ( settings1->SLAC ) {
2081 Vector xaxis(1., 0., 0.);
2083 interaction1->
nnu = xaxis.RotateY(bn1->theta_bn-settings1->SLAC_HORIZDIST/EarthModel::R_EARTH);
2084 interaction1->
nnu = interaction1->
nnu.RotateZ(bn1->phi_bn);
2085 costhetanu=cos(interaction1->
nnu.Theta());
2087 if (settings1->BORESIGHTS) {
2088 fslac_viewangles << bn1->sslacpositions[bn1->islacposition] <<
"\n";
2089 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2090 for(
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
2091 viewangle_eachboresight[ilayer][ifold]=acos(interaction1->
nnu.Dot(ray1->nrf_iceside_eachboresight[4][ilayer][ifold]));
2092 fslac_viewangles << ilayer <<
"\t" << ifold <<
"\t" << (viewangle_eachboresight[ilayer][ifold]-sig1->changle)*DEGRAD <<
"\n";
2104 if (!settings1->UNBIASED_SELECTION)
2106 interaction1->
r_in = antarctica->WhereDoesItEnter(interaction1->posnu, interaction1->
nnu);
2110 taus1->
GetTauWeight(primary1, settings1, antarctica, interaction1, pnu, 1, ptauf, interaction1->crust_entered);
2113 if (!settings1->UNBIASED_SELECTION)
2115 antarctica->Getchord(
false, len_int_kgm2, interaction1->
r_in, interaction1->pathlength_inice, settings1->UNBIASED_SELECTION == 0, interaction1->posnu, inu, interaction1->
chord, interaction1->
weight_nu_prob, interaction1->
weight_nu, interaction1->nearthlayers, myair, interaction1->
total_kgm2, interaction1->crust_entered, interaction1->mantle_entered, interaction1->core_entered);
2121 nutauweightsum +=nutauweight;
2122 tauweightsum +=tauweight;
2123 double xrndm=getRNG(RNG_XRNDM)->Rndm();
2134 TauPtr->
ptauf = ptauf;
2138 if (mytaus_tree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2139 mytaus_tree->Fill();
2147 sec1->GetEMFrac(settings1, interaction1->nuflavor(), interaction1->current(), taudecay, elast_y, hy, pnu, inu,emfrac, hadfrac, n_interactions, tauweighttrigger);
2148 havent_set_frac =
false;
2150 #ifdef ANITA3_EVENTREADER 2151 if (truthNuPtr) truthNuPtr->setFrac(hadfrac,emfrac);
2153 if (emfrac+hadfrac>1.000001) {
2154 cout <<
"Warning: " << inu <<
" " << emfrac+hadfrac <<
"\n";
2159 sumfrac=emfrac+hadfrac;
2161 if (tree7->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2164 vmmhz1m_visible = (emfrac+hadfrac)*vmmhz1m_max;
2167 if (interaction1->nuflavor()==
"numu" && bn1->WHICHPATH != 3 && !settings1->ONLYFINAL && settings1->HIST==1 && fraction_sec_muons->GetEntries()<settings1->HIST_MAX_ENTRIES) {
2168 fraction_sec_muons->Fill(emfrac+hadfrac, weight);
2169 n_sec_muons->Fill((
double)n_interactions);
2172 if (interaction1->nuflavor()==
"nutau" && bn1->WHICHPATH != 3 && !settings1->ONLYFINAL && settings1->HIST==1 && fraction_sec_taus->GetEntries()<settings1->HIST_MAX_ENTRIES) {
2173 fraction_sec_taus->Fill(emfrac+hadfrac, weight);
2174 n_sec_taus->Fill((
double)n_interactions);
2178 if(sec1->secondbang && sec1->interestedintaus) {
2179 ptau=(1-elast_y)*pnu;
2190 if(sec1->secondbang && sec1->interestedintaus) {
2191 vmmhz1m_max=sig1->GetVmMHz1m(ptau, anita1->FREQ_HIGH);
2192 sig1->GetSpread(ptau, emfrac, hadfrac, anita1->FREQ_LOW, deltheta_em_max, deltheta_had_max);
2195 vmmhz1m_max=sig1->GetVmMHz1m(pnu, anita1->FREQ_HIGH);
2196 sig1->GetSpread(pnu, emfrac, hadfrac, anita1->FREQ_LOW,
2197 deltheta_em_max, deltheta_had_max);
2200 if (jaimetree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2209 bestcase_atten=exp(interaction1->
altitude_int/MAX_ATTENLENGTH);
2214 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/((hadfrac+emfrac)*vmmhz1m_max*bestcase_atten/interaction1->
r_fromballoon[whichray]*heff_max*anita1->bwmin/1.E6)>settings1->CHANCEINHELL_FACTOR && !settings1->SKIPCUTS) {
2223 count1->nnottoosmall[whichray]++;
2229 chengji=ray1->nrf_iceside[4]*ray1->nrf_iceside[0];
2233 ray1->nrf_iceside[4] = ray1->nrf_iceside[4] - 2*chengji*ray1->nrf_iceside[0];
2236 if ( tautrigger == 0 ) {
2237 if ( ! src_model && ( !settings1->UNBIASED_SELECTION ) && ( !settings1->SLAC ) ) {
2238 err = GetDirection(settings1, interaction1, ray1->nrf_iceside[4], deltheta_em_max, deltheta_had_max, emfrac, hadfrac, vmmhz1m_max*bestcase_atten, interaction1->
r_fromballoon[whichray], ray1, sig1, interaction1->posnu, anita1, bn1, interaction1->
nnu, costhetanu, theta_threshold);
2241 else if (src_model || settings1->HORIZON_OFFSET > -999)
2243 interaction1->
nnu = force_dir;
2245 costhetanu = cos(force_dir.Theta());
2246 theta_threshold = 1;
2249 else if (settings1->SLAC) {
2250 Vector xaxis(1., 0., 0.);
2251 interaction1->
nnu = xaxis.RotateY(bn1->theta_bn-settings1->SLAC_HORIZDIST/EarthModel::R_EARTH);
2252 interaction1->
nnu = interaction1->
nnu.RotateZ(bn1->phi_bn);
2253 costhetanu=cos(interaction1->
nnu.Theta());
2255 if (settings1->BORESIGHTS) {
2256 fslac_viewangles << bn1->sslacpositions[bn1->islacposition] <<
"\n";
2257 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2258 for(
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
2259 viewangle_eachboresight[ilayer][ifold]=acos(interaction1->
nnu.Dot(ray1->nrf_iceside_eachboresight[4][ilayer][ifold]));
2260 fslac_viewangles << ilayer <<
"\t" << ifold <<
"\t" << (viewangle_eachboresight[ilayer][ifold]-sig1->changle)*DEGRAD <<
"\n";
2269 #ifdef ANITA3_EVENTREADER 2272 if (truthNuPtr) truthNuPtr->setDir(interaction1->
nnu.GetX(), interaction1->
nnu.GetY(), interaction1->
nnu.GetZ());
2273 havent_set_dir =
false;
2278 viewangle = GetViewAngle(ray1->nrf_iceside[4], interaction1->
nnu);
2279 if(viewangle>1.57 && !settings1->SKIPCUTS) {
2282 count1->nviewangle_lt_90[whichray]++;
2284 if (!Ray::WhereDoesItLeave(interaction1->posnu, interaction1->
nnu, antarctica, interaction1->
nuexit))
2289 GetBalloonLocation(interaction1, ray1, bn1, antarctica);
2295 theta_threshold_deg=theta_threshold*DEGRAD;
2298 n_nutraject_ontheground =
Vector(bn1->n_east*interaction1->
nnu, bn1->n_north*interaction1->
nnu, bn1->n_bn*interaction1->
nnu);
2300 cosviewangle=cos(viewangle);
2301 viewangle_deg=viewangle*DEGRAD;
2302 dviewangle_deg=(sig1->changle-viewangle)*DEGRAD;
2304 if (viewangletree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2305 viewangletree->Fill();
2307 if (neutrino_positiontree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2308 neutrino_positiontree->Fill();
2312 ray1->nrf_iceside[4] = ray1->nrf_iceside[4] + 2*chengji*ray1->nrf_iceside[0];
2316 count1->nbadfracs[whichray]++;
2317 cout<<
"err==0, so leaving.\n";
2320 count1->ngoodfracs[whichray]++;
2323 nsigma_em_threshold=theta_threshold/deltheta_em_max;
2324 nsigma_had_threshold=theta_threshold/deltheta_had_max;
2334 index_distance=(int)(bn1->r_bn.SurfaceDistance(interaction1->posnu, bn1->surface_under_balloon) / (700000./(double)NBINS_DISTANCE));
2337 if (tautrigger==0 && !settings1->UNBIASED_SELECTION){
2338 interaction1->
r_in = antarctica->WhereDoesItEnter(interaction1->posnu, interaction1->
nnu);
2339 antarctica->Getchord(
false, len_int_kgm2, interaction1->
r_in, interaction1->pathlength_inice, settings1->UNBIASED_SELECTION ==0, interaction1->posnu, inu, interaction1->
chord, interaction1->
weight_nu_prob,
2340 interaction1->
weight_nu, interaction1->nearthlayers, myair, interaction1->
total_kgm2, interaction1->crust_entered, interaction1->mantle_entered, interaction1->core_entered);
2345 double chord_kgm2_test=interaction1->posnu.
Distance(interaction1->
r_in)*sig1->RHOMEDIUM;
2347 double weight_test=0;
2350 IsAbsorbed(chord_kgm2_test, len_int_kgm2, weight_test);
2358 if ( bn1->WHICHPATH != 4 && settings1->FORSECKEL != 1 && !settings1->SKIPCUTS) {
2359 if (weight_test < settings1->CUTONWEIGHTS) {
2363 count_chanceofsurviving++;
2367 theta_in=interaction1->
r_in.Theta();
2368 lat_in=-90+theta_in*DEGRAD;
2372 myair=GetThisAirColumn(settings1, interaction1->
r_in, interaction1->
nnu, interaction1->posnu, col1, cosalpha, mytheta, cosbeta0, mybeta);
2377 if (!settings1->FORSECKEL && !settings1->UNBIASED_SELECTION) {
2378 if (!antarctica->WhereDoesItEnterIce(interaction1->posnu, interaction1->
nnu, len_int_kgm2/sig1->RHOMEDIUM/10., interaction1->
r_enterice)) {
2380 if (antarctica->OutsideAntarctica(interaction1->
r_enterice)) {
2381 cout<<
"Warning! Neutrino enters beyond continent, program is rejecting neutrino! inu = "<<inu<<endl;
2382 if (bn1->WHICHPATH==3)
2383 cout<<
"Warning! Neutrino enters beyond continent, program is rejecting neutrino!"<<endl;
2390 count1->nentersice[whichray]++;
2401 if(sec1->secondbang && sec1->interestedintaus) {
2402 cout <<
"Need to bring back GetFirstBang before you can simulate taus.\n";
2403 cout <<
"I removed it because it required EarthModel and I wanted Secondaries to be a stand-alone class to use in the embedded simulation.\n";
2418 if (tree6->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2434 if (bn1->WHICHPATH!=4 && interaction1->
weight_bestcase<settings1->CUTONWEIGHTS && !settings1->SKIPCUTS && !settings1->FORSECKEL) {
2435 if (bn1->WHICHPATH==3)
2436 cout<<
"Neutrino is getting absorbed and thrown out!"<<endl;
2441 count1->nabsorbed[whichray]++;
2444 count1->nraywithincontinent1[whichray]++;
2451 bestcase_atten=exp(-1*ray1->rfexit[1].
Distance(interaction1->posnu)/MAX_ATTENLENGTH);
2454 bestcase_atten=exp(-1*ray1->rfexit[1].
Distance(interaction1->posnu_down)/MAX_ATTENLENGTH);
2456 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/((hadfrac+emfrac)*vmmhz1m_max*bestcase_atten/interaction1->
r_fromballoon[whichray]*heff_max*anita1->bwmin/1.E6)>settings1->CHANCEINHELL_FACTOR && !settings1->SKIPCUTS && !settings1->FORSECKEL) {
2457 if (bn1->WHICHPATH==3)
2458 cout<<
"Event rejected. Check."<<endl;
2462 count_chanceinhell0++;
2465 count1->nraypointsup2[whichray]++;
2467 double nbelowsurface;
2469 if (settings1->FIRN)
2470 nbelowsurface=NFIRN;
2472 nbelowsurface=sig1->NICE;
2475 if (!settings1->ROUGHNESS && TIR(ray1->nsurf_rfexit, ray1->nrf_iceside[3], nbelowsurface, sig1->N_AIR)) {
2478 count1->nnottir[whichray]++;
2482 ray1->GetRFExit(settings1, anita1, whichray, interaction1->posnu, interaction1->posnu_down, bn1->r_bn, bn1->r_boresights, 2, antarctica);
2484 count1->nraywithincontinent2[whichray]++;
2489 theta_nnu_atbn = interaction1->
nnu.Angle(bn1->r_bn);
2491 theta_rf_atbn = ray1->n_exit2bn[2].Angle(bn1->r_bn);
2493 theta_rf_atbn_measured = theta_rf_atbn+getRNG(RNG_THETA_RF_RESOLUTION)->Gaus()*anita1->SIGMA_THETA;
2494 interaction1->
r_exit2bn=bn1->r_bn.Distance(ray1->rfexit[2]);
2497 theta_nnu_rf_diff_atbn = theta_nnu_atbn - theta_rf_atbn;
2503 nnu_xy.SetX((interaction1->
nnu[0]) - (bn1->r_bn[0])*((interaction1->
nnu*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2504 nnu_xy.SetY((interaction1->
nnu[1]) - (bn1->r_bn[1])*((interaction1->
nnu*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2505 nnu_xy.SetZ((interaction1->
nnu[2]) - (bn1->r_bn[2])*((interaction1->
nnu*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2507 rf_xy.SetX((ray1->n_exit2bn[2][0]) - (bn1->r_bn[0])*((ray1->n_exit2bn[2]*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2508 rf_xy.SetY((ray1->n_exit2bn[2][1]) - (bn1->r_bn[1])*((ray1->n_exit2bn[2]*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2509 rf_xy.SetZ((ray1->n_exit2bn[2][2]) - (bn1->r_bn[2])*((ray1->n_exit2bn[2]*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2511 r_bn_xy.SetX(bn1->r_bn[0]);
2512 r_bn_xy.SetY(bn1->r_bn[1]);
2514 phi_nnu_atbn = nnu_xy.Angle(r_bn_xy);
2515 phi_rf_atbn = rf_xy.Angle(r_bn_xy);
2525 phi_nnu_rf_diff_atbn = phi_nnu_atbn - phi_rf_atbn;
2532 if((settings1->WHICH == 2 || settings1->WHICH == 6) && theta_rf_atbn < 0.3790091) {
2536 if (!antarctica->AcceptableRfexit(ray1->nsurf_rfexit, ray1->rfexit[2], ray1->n_exit2bn[2])){
2537 if (bn1->WHICHPATH==3)
2538 cout<<
"Should look at this. Not expecting to be here."<<endl;
2543 count1->nacceptablerf[whichray]++;
2546 diff_3tries=ray1->rfexit[1].
Distance(ray1->rfexit[2]);
2548 if (tree5->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1 && bn1->WHICHPATH != 3)
2553 if (diff_3tries>10) {
2556 count1->nconverges[whichray]++;
2558 n_pol = GetPolarization(interaction1->
nnu, ray1->nrf_iceside[4]);
2561 if (settings1->BORESIGHTS) {
2562 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2563 for(
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
2564 n_pol_eachboresight[ilayer][ifold]=GetPolarization(interaction1->
nnu, ray1->nrf_iceside_eachboresight[4][ilayer][ifold]);
2570 if (settings1->FIRN){
2573 GetFresnel(rough1, settings1->ROUGHNESS, ray1->nsurf_rfexit, ray1->nrf_iceside[3], n_pol, ray1->nrf_iceside[4], vmmhz1m_max, emfrac, hadfrac, deltheta_em_max, deltheta_had_max, t_coeff_pokey, t_coeff_slappy, fresnel1, mag1);
2574 if (bn1->WHICHPATH==4)
2575 cout <<
"Lentenin factor is " << 1./mag1 <<
"\n";
2580 vmmhz1m_fresneledonce = vmmhz1m_max/mag1;
2583 GetFresnel(rough1, settings1->ROUGHNESS, ray1->nsurf_rfexit, ray1->n_exit2bn[2], n_pol, ray1->nrf_iceside[3], vmmhz1m_fresneledonce, emfrac, hadfrac, deltheta_em_max, deltheta_had_max, t_coeff_pokey, t_coeff_slappy, fresnel2, mag2);
2586 vmmhz1m_fresneledtwice=vmmhz1m_fresneledonce*fresnel2*mag2;
2588 if (settings1->BORESIGHTS) {
2589 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2590 for(
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
2591 GetFresnel(rough1, settings1->ROUGHNESS, ray1->nsurf_rfexit, ray1->n_exit2bn_eachboresight[2][ilayer][ifold],
2592 n_pol_eachboresight[ilayer][ifold], ray1->nrf_iceside_eachboresight[3][ilayer][ifold],
2593 vmmhz1m_max, emfrac, hadfrac, deltheta_em_max, deltheta_had_max, t_coeff_pokey, t_coeff_slappy,
2594 fresnel1_eachboresight[ilayer][ifold], mag1_eachboresight[ilayer][ifold]);
2601 if (bn1->WHICHPATH==4)
2602 cout<<
"firn-air interface: fresnel2, mag2 are "<<fresnel2<<
" "<< mag2 <<
"\n";
2606 sig1->GetSpread(pnu, emfrac, hadfrac, (anita1->bwslice_min[2]+anita1->bwslice_max[2])/2., deltheta_em_mid2, deltheta_had_mid2);
2608 GetFresnel(rough1, settings1->ROUGHNESS, ray1->nsurf_rfexit, ray1->n_exit2bn[2], n_pol, ray1->nrf_iceside[4], vmmhz1m_max, emfrac, hadfrac, deltheta_em_mid2, deltheta_had_mid2, t_coeff_pokey, t_coeff_slappy, fresnel1, mag1);
2610 vmmhz1m_fresneledtwice = vmmhz1m_max*fresnel1*mag1;
2612 if (settings1->BORESIGHTS) {
2613 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2614 for(
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
2615 GetFresnel(rough1, settings1->ROUGHNESS, ray1->nsurf_rfexit, ray1->n_exit2bn_eachboresight[2][ilayer][ifold], n_pol_eachboresight[ilayer][ifold], ray1->nrf_iceside_eachboresight[4][ilayer][ifold], vmmhz1m_max, emfrac, hadfrac, deltheta_em_max, deltheta_had_max, t_coeff_pokey, t_coeff_slappy, fresnel1_eachboresight[ilayer][ifold], mag1_eachboresight[ilayer][ifold]);
2622 if(settings1->ROUGHNESS){
2630 int num_validscreenpoints = 0;
2632 Vector vec_pos_current_to_balloon;
2636 Vector vec_nnu_to_impactPoint;
2639 double pol_perp_inc, pol_parl_inc;
2640 Vector vec_local_grnd_perp;
2641 Vector vec_local_grnd_parl;
2642 double pol_perp_trans, pol_parl_trans;
2646 double time_reference_specular, time_reference_local;
2647 double pathlength_local;
2648 double viewangle_local;
2649 double azimuth_local;
2651 double theta_0_local;
2652 double tcoeff_perp_polperp, tcoeff_parl_polperp;
2653 double tcoeff_perp_polparl, tcoeff_parl_polparl;
2654 double power_perp_polperp, power_parl_polperp;
2655 double power_perp_polparl, power_parl_polparl;
2656 power_perp_polperp = power_parl_polperp = power_perp_polparl = power_parl_polparl = 0.;
2657 double fresnel_r, mag_r;
2659 Vector npol_local_inc, npol_local_trans;
2667 if (settings1->FIRN)
2668 time_reference_specular = (interaction1->posnu.
Distance(ray1->rfexit[2])*NFIRN / CLIGHT) + (ray1->rfexit[2].
Distance(bn1->r_bn)/CLIGHT);
2670 time_reference_specular = (interaction1->posnu.
Distance(ray1->rfexit[2])*NICE / CLIGHT) + (ray1->rfexit[2].
Distance(bn1->r_bn)/CLIGHT);
2672 double slopeyx, slopeyy, slopeyz, rtemp;
2678 double basescreenedgelength = settings1->SCREENEDGELENGTH;
2679 double grd_stepsize = settings1->SCREENSTEPSIZE;
2681 if(settings1->ROUGHSIZE>0)
2682 grd_nsteps = int(basescreenedgelength/2. / grd_stepsize);
2691 panel1->ResetParameters();
2693 panel1->SetNsamples( grd_nsteps );
2694 panel1->SetEdgeLength( grd_stepsize );
2696 panel1->SetCentralPoint( ray1->rfexit[2] );
2697 vec_localnormal = antarctica->GetSurfaceNormal(ray1->rfexit[2]).Unit();
2698 panel1->SetNormal( vec_localnormal );
2699 panel1->SetCosineProjectionFactor( 1. );
2701 panel1->SetUnitY( (vec_localnormal.Cross(ray1->n_exit2bn[2])).Unit() );
2702 panel1->SetUnitX( (panel1->GetUnitY().Cross(vec_localnormal)).Unit() );
2705 for (
int ii= -2*panel1->GetNsamples(); ii< 2*panel1->GetNsamples()+1; ii++){
2706 for (
int jj= -2*panel1->GetNsamples(); jj< 2*panel1->GetNsamples()+1; jj++){
2710 Emag_local = vmmhz1m_max;
2711 taperfactor = fresnel_r = mag_r = 1.;
2712 tcoeff_perp_polparl = tcoeff_parl_polparl = 0.;
2713 tcoeff_perp_polperp = tcoeff_parl_polperp = 0.;
2714 pos_projectedImpactPoint = panel1->GetPosition(ii, jj);
2715 vec_pos_current_to_balloon =
Vector( bn1->r_bn[0] - pos_projectedImpactPoint[0], bn1->r_bn[1] - pos_projectedImpactPoint[1], bn1->r_bn[2] - pos_projectedImpactPoint[2] );
2718 vec_localnormal = antarctica->GetSurfaceNormal(pos_projectedImpactPoint).Unit();
2719 if (settings1->SLOPEY) {
2720 slopeyx=ray1->slopeyx;
2721 slopeyy=ray1->slopeyy;
2722 slopeyz=ray1->slopeyz;
2723 ntemp2 = vec_localnormal + slopeyx*xaxis + slopeyy*yaxis + slopeyz*zaxis;
2724 ntemp2 = ntemp2 / ntemp2.Mag();
2725 rtemp= ntemp2 * vec_localnormal;
2727 vec_localnormal = ntemp2;
2731 vec_nnu_to_impactPoint =
Vector( pos_projectedImpactPoint[0]-interaction1->posnu[0], pos_projectedImpactPoint[1]-interaction1->posnu[1], pos_projectedImpactPoint[2]-interaction1->posnu[2] ).Unit();
2733 vec_grndcomp2IP = (vec_nnu_to_impactPoint - (vec_nnu_to_impactPoint.Dot(vec_localnormal)*vec_localnormal)).Unit();
2734 vec_grndcomp2bln = (vec_pos_current_to_balloon - (vec_pos_current_to_balloon.Dot(vec_localnormal)*vec_localnormal)).Unit();
2735 temp_a = vec_localnormal.Cross(vec_pos_current_to_balloon).Unit();
2736 azimuth_local = vec_grndcomp2IP.Angle(vec_grndcomp2bln);
2737 if( temp_a.Dot(vec_nnu_to_impactPoint) < 0 )
2738 azimuth_local *= -1.;
2739 if( panel1->GetCentralPoint().Distance(pos_projectedImpactPoint)<0.75*grd_stepsize ){
2743 theta_local = vec_localnormal.Angle( (
const Vector)vec_pos_current_to_balloon );
2744 theta_0_local = vec_localnormal.Angle(vec_nnu_to_impactPoint);
2747 if( isnan(theta_local) | isnan(theta_0_local) | isnan(azimuth_local) ){
2750 viewangle_local = GetViewAngle(vec_nnu_to_impactPoint, interaction1->
nnu);
2753 deltheta_em[0]=deltheta_em_max*anita1->FREQ_LOW/anita1->freq[0];
2754 deltheta_had[0]=deltheta_had_max*anita1->FREQ_LOW/anita1->freq[0];
2755 sig1->TaperVmMHz(viewangle_local, deltheta_em[0], deltheta_had[0], emfrac, hadfrac, taperfactor, vmmhz_em[0]);
2764 if (settings1->FIRN)
2765 rough1->InterpolatePowerValue(power_perp_polperp, power_parl_polperp, power_perp_polparl, power_parl_polparl, theta_0_local*180./PI, theta_local*180./PI, azimuth_local *180./PI);
2767 rough1->InterpolatePowerValue(power_perp_polperp, power_parl_polperp, power_perp_polparl, power_parl_polparl, theta_0_local*180./PI, theta_local*180./PI, azimuth_local *180./PI);
2770 if( (power_perp_polperp==0.)&(power_parl_polperp==0.)&(power_perp_polparl==0.)&(power_parl_polparl==0.) ){
2774 if (settings1->FIRN){
2775 tcoeff_perp_polparl = sqrt(power_perp_polparl);
2776 tcoeff_parl_polparl = sqrt(power_parl_polparl);
2777 tcoeff_perp_polperp = sqrt(power_perp_polperp);
2778 tcoeff_parl_polperp = sqrt(power_parl_polperp);
2781 tcoeff_perp_polparl = sqrt(power_perp_polparl);
2782 tcoeff_parl_polparl = sqrt(power_parl_polparl);
2783 tcoeff_perp_polperp = sqrt(power_perp_polperp);
2784 tcoeff_parl_polperp = sqrt(power_parl_polperp);
2792 pathlength_local = interaction1->posnu.
Distance(pos_projectedImpactPoint) + pos_projectedImpactPoint.
Distance(bn1->r_bn);
2794 Emag_local /= pathlength_local ;
2796 Attenuate(antarctica, settings1, Emag_local, interaction1->posnu.
Distance(pos_projectedImpactPoint), interaction1->posnu);
2801 npol_local_inc = GetPolarization(interaction1->
nnu, vec_nnu_to_impactPoint).Unit();
2802 vec_inc_perp = (vec_localnormal.Cross(vec_nnu_to_impactPoint)).Unit();
2803 vec_inc_parl = (vec_nnu_to_impactPoint.Cross(vec_inc_perp)).Unit();
2804 pol_perp_inc = npol_local_inc * vec_inc_perp;
2805 pol_parl_inc = npol_local_inc * vec_inc_parl;
2807 pol_perp_trans = pol_perp_inc * tcoeff_perp_polperp + pol_perp_inc * tcoeff_perp_polparl;
2808 pol_parl_trans = pol_parl_inc * tcoeff_parl_polparl + pol_parl_inc * tcoeff_parl_polperp;
2810 vec_local_grnd_perp = (vec_localnormal.Cross(vec_pos_current_to_balloon)).Unit();
2811 vec_local_grnd_parl = (vec_pos_current_to_balloon.Cross(vec_local_grnd_perp)).Unit();
2814 npol_local_trans= (pol_perp_trans*vec_local_grnd_perp + pol_parl_trans*vec_local_grnd_parl).Unit();
2817 if( isnan(npol_local_trans[0]) ){
2822 fresnel_r = sqrt( pow(vmmhz1m_max*pol_perp_trans,2) + pow(vmmhz1m_max*pol_parl_trans,2) ) / vmmhz1m_max;
2823 mag_r = sqrt( tan(theta_0_local) / tan(theta_local) );
2824 Emag_local *= fresnel_r * mag_r;
2826 if (settings1->FIRN)
2827 time_reference_local = (interaction1->posnu.
Distance(pos_projectedImpactPoint)*NFIRN / CLIGHT) + (pos_projectedImpactPoint.
Distance(bn1->r_bn)/CLIGHT);
2829 time_reference_local = (interaction1->posnu.
Distance(pos_projectedImpactPoint)*NICE / CLIGHT) + (pos_projectedImpactPoint.
Distance(bn1->r_bn)/CLIGHT);
2831 num_validscreenpoints++;
2834 panel1->AddVmmhz0(Emag_local);
2835 panel1->AddVec2bln(vec_pos_current_to_balloon);
2836 panel1->AddPol(npol_local_trans);
2837 panel1->AddDelay( time_reference_specular - time_reference_local );
2838 panel1->AddImpactPt(pos_projectedImpactPoint);
2839 panel1->AddViewangle(viewangle_local);
2840 panel1->AddIncidenceAngle(theta_0_local);
2841 panel1->AddTransmissionAngle(theta_local);
2842 panel1->AddWeight( (panel1->GetEdgeLength() / panel1->GetNsamples()) * (panel1->GetEdgeLength() / panel1->GetNsamples()) );
2843 panel1->AddFacetLength(panel1->GetEdgeLength() / panel1->GetNsamples());
2844 panel1->AddTparallel_polParallel(tcoeff_parl_polparl);
2845 panel1->AddTperpendicular_polParallel(tcoeff_perp_polparl);
2846 panel1->AddTparallel_polPerpendicular(tcoeff_parl_polperp);
2847 panel1->AddTperpendicular_polPerpendicular(tcoeff_perp_polperp);
2858 panel1->SetNvalidPoints(num_validscreenpoints);
2861 double validScreenSummedArea = 0.;
2862 double vmmhz_local_array[Anita::NFREQ];
2863 for (
int jj=0; jj<panel1->GetNvalidPoints(); jj++){
2865 sig1->GetVmMHz(panel1->GetVmmhz0(jj), vmmhz1m_max, pnu, anita1->freq, anita1->NOTCH_MIN, anita1->NOTCH_MAX, vmmhz_local_array, Anita::NFREQ);
2867 for (
int k=0;k<Anita::NFREQ;k++) {
2868 deltheta_em[k]=deltheta_em_max*anita1->FREQ_LOW/anita1->freq[k];
2869 deltheta_had[k]=deltheta_had_max*anita1->FREQ_LOW/anita1->freq[k];
2870 sig1->TaperVmMHz(panel1->GetViewangle(jj), deltheta_em[k], deltheta_had[k], emfrac, hadfrac, vmmhz_local_array[k], vmmhz_em[k]);
2871 panel1->AddVmmhz_freq(vmmhz_local_array[k]);
2874 validScreenSummedArea += panel1->GetWeight(jj);
2876 panel1->SetWeightNorm(validScreenSummedArea);
2878 for(
int jj=0; jj<panel1->GetNvalidPoints(); jj++){
2879 vmmhz_max = Tools::dMax(vmmhz_max, panel1->GetVmmhz_freq(jj*Anita::NFREQ));
2885 if( settings1->ROUGHNESS && !panel1->GetNvalidPoints() ){
2892 if(settings1->CHANCEINHELL_FACTOR*vmmhz1m_fresneledtwice*heff_max*0.5*(anita1->bwmin/1.E6)<anita1->maxthreshold*anita1->VNOISE[0]/10.&& !settings1->SKIPCUTS) {
2893 if (bn1->WHICHPATH==3)
2894 cout<<
"Event is undetectable. Leaving loop."<<endl;
2898 count1->nchanceinhell_fresnel[whichray]++;
2903 diffexit=ray1->rfexit[0].
Distance(ray1->rfexit[1]);
2904 diffnorm=acos(ray1->nsurf_rfexit[0]*ray1->nsurf_rfexit[1]);
2905 diffrefr=acos(ray1->nrf_iceside[4]*ray1->nrf_iceside[0]);
2909 if (!settings1->ROUGHNESS) {
2913 vmmhz_max=ScaleVmMHz(vmmhz1m_fresneledtwice, interaction1->posnu, bn1->r_bn, ray1->rfexit[2]);
2915 vmmhz_max=ScaleVmMHz(vmmhz1m_fresneledtwice, interaction1->posnu_down, bn1->r_bn, ray1->rfexit[2]);
2924 if (!settings1->ROUGHNESS){
2925 if (settings1->CHANCEINHELL_FACTOR*vmmhz_max*heff_max*0.5*(anita1->bwmin/1.E6)<anita1->maxthreshold*anita1->VNOISE[0]/10. && !settings1->SKIPCUTS) {
2926 if (bn1->WHICHPATH==3)
2927 cout<<
"Event is undetectable. Leaving loop."<<endl;
2932 count1->nchanceinhell_1overr[whichray]++;
2935 if (!settings1->ROUGHNESS) {
2937 rflength=interaction1->posnu.
Distance(ray1->rfexit[2]);
2940 rflength=interaction1->posnu_down.
Distance(ray1->rfexit[2]);
2944 if (bn1->WHICHPATH==4)
2945 cout <<
"rflength is " << rflength <<
"\n";
2948 if (!settings1->ROUGHNESS) {
2950 Attenuate(antarctica, settings1, vmmhz_max, rflength, interaction1->posnu);
2952 Attenuate_down(antarctica, settings1, vmmhz_max, ray1->rfexit[2], interaction1->posnu, interaction1->posnu_down);
2959 if (tree2->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1 && bn1->WHICHPATH != 3)
2966 if (!settings1->ROUGHNESS){
2967 if (settings1->CHANCEINHELL_FACTOR*vmmhz_max*heff_max*0.5*(anita1->bwmin/1.E6)<anita1->maxthreshold*anita1->VNOISE[0]/10. && !settings1->SKIPCUTS) {
2968 if (bn1->WHICHPATH==3)
2969 cout<<
"Event is undetectable. Leaving loop."<<endl;
2975 count1->nchanceinhell[whichray]++;
2978 if (tree3->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1 && bn1->WHICHPATH != 3)
2991 if (!settings1->ROUGHNESS){
2992 if (settings1->FORSECKEL==1)
2993 sig1->SetNDepth(sig1->NICE);
2995 sig1->GetVmMHz(vmmhz_max, vmmhz1m_max, pnu, anita1->freq, anita1->NOTCH_MIN, anita1->NOTCH_MAX, vmmhz, Anita::NFREQ);
3005 undogaintoheight_e=0;
3006 undogaintoheight_h=0;
3008 for (
int k=0;k<4;k++) {
3009 undogaintoheight_e_array[k]=0.;
3010 undogaintoheight_h_array[k]=0.;
3012 true_efield_array[k]=0.;
3013 rec_efield_array[k]=0.;
3020 if (!settings1->ROUGHNESS){
3022 double rtemp=Tools::dMin((viewangle-sig1->changle)/(deltheta_em_max), (viewangle-sig1->changle)/(deltheta_had_max));
3023 if (rtemp>Signal::VIEWANGLE_CUT && !settings1->SKIPCUTS) {
3027 count1->nviewanglecut[whichray]++;
3032 for (
int k=0;k<Anita::NFREQ;k++) {
3033 deltheta_em[k]=deltheta_em_max*anita1->FREQ_LOW/anita1->freq[k];
3034 deltheta_had[k]=deltheta_had_max*anita1->FREQ_LOW/anita1->freq[k];
3036 if (settings1->FORSECKEL==1) {
3037 for (
int iviewangle=0;iviewangle<NVIEWANGLE;iviewangle++) {
3039 vmmhz_temp=vmmhz[k]*vmmhz1m_max/vmmhz_max;
3041 viewangle_temp=viewangles[iviewangle];
3044 sig1->TaperVmMHz(viewangle_temp, deltheta_em[k], deltheta_had[k], emfrac, hadfrac, vmmhz_temp, djunk);
3045 forseckel[iviewangle][k]=vmmhz_temp;
3054 sig1->TaperVmMHz(viewangle, deltheta_em[k], deltheta_had[k], emfrac, hadfrac, vmmhz[k], vmmhz_em[k]);
3063 vmmhz_lowfreq=vmmhz[0];
3069 if (sig1->logscalefactor_taper>maxtaper)
3070 maxtaper=sig1->logscalefactor_taper;
3072 if (interaction1->nuflavor()==
"nue") pdgcode = 12;
3073 else if (interaction1->nuflavor()==
"numu") pdgcode = 14;
3074 else if (interaction1->nuflavor()==
"nutau") pdgcode = 16;
3076 if (settings1->HIST==1 && !settings1->ONLYFINAL && bn1->WHICHPATH != 3 && k==Anita::NFREQ/2 && tree18->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3081 if (bn1->WHICHPATH == 3)
3082 interaction1->
banana_volts += vmmhz[k]*(settings1->BW/(double)Anita::NFREQ/1.E6);
3090 if (bn1->WHICHPATH==3 && interaction1->
banana_volts != 0 && settings1->HIST && banana_tree->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3091 banana_tree->Fill();
3094 else if (bn1->WHICHPATH==3 && interaction1->
banana_volts == 0) {
3106 if (settings1->CHANCEINHELL_FACTOR*Tools::dMax(vmmhz, Anita::NFREQ)*heff_max*0.5*(anita1->bwmin/1.E6) < anita1->maxthreshold*anita1->VNOISE[0]/10. && !settings1->SKIPCUTS) {
3115 if(!settings1->ROUGHNESS){
3116 vmmhz_max=Tools::dMax(vmmhz, Anita::NFREQ);
3117 vmmhz_min=Tools::dMin(vmmhz, Anita::NFREQ);
3120 count1->nchanceinhell2[whichray]++;
3127 if (settings1->USEDEADTIME){
3128 if ( (getRNG(RNG_DEADTIME)->Uniform(1)<anita1->deadTime) ){
3130 if (settings1->MINBIAS!=1) {
3136 count1->ndeadtime[whichray]++;
3138 Tools::Zero(sumsignal_aftertaper, 5);
3149 if(!settings1->ROUGHNESS){
3150 panel1->SetNvalidPoints(1);
3151 for (
int k=0;k<Anita::NFREQ;k++) {
3153 panel1->AddVmmhz_freq(vmmhz[k]);
3155 panel1->AddVmmhz0(vmmhz[0]);
3156 panel1->AddVec2bln(ray1->n_exit2bn[2]);
3157 panel1->AddPol(n_pol);
3158 panel1->AddDelay( 0. );
3159 panel1->AddImpactPt(ray1->rfexit[2]);
3160 panel1->AddViewangle(viewangle);
3161 if(settings1->FIRN){
3162 panel1->AddIncidenceAngle(ray1->nsurf_rfexit.Angle(ray1->nrf_iceside[3]));
3165 panel1->AddIncidenceAngle(ray1->nsurf_rfexit.Angle(ray1->nrf_iceside[4]));
3167 panel1->AddTransmissionAngle(ray1->nsurf_rfexit.Angle(ray1->n_exit2bn[2]));
3168 panel1->AddWeight( 1. );
3169 panel1->SetWeightNorm( 1. );
3170 panel1->AddFacetLength( 1. );
3171 panel1->AddTparallel_polParallel(t_coeff_pokey);
3172 panel1->AddTperpendicular_polPerpendicular(t_coeff_slappy);
3174 panel1->AddTparallel_polPerpendicular(0.);
3175 panel1->AddTperpendicular_polParallel(0.);
3177 for (
int k=0;k<Anita::NFREQ;k++) {
3178 if (bn1->WHICHPATH==4)
3179 IntegrateBands(anita1, k, panel1, anita1->freq, vmmhz1m_max/(vmmhz_max*1.E6), sumsignal_aftertaper);
3186 Tools::Zero(anita1->arrival_times[0], Anita::NLAYERS_MAX*
Anita::NPHI_MAX);
3187 Tools::Zero(anita1->arrival_times[1], Anita::NLAYERS_MAX*
Anita::NPHI_MAX);
3188 if (!settings1->TRIGGEREFFSCAN){
3189 if(settings1->BORESIGHTS)
3190 anita1->GetArrivalTimesBoresights(ray1->n_exit2bn_eachboresight[2]);
3192 anita1->GetArrivalTimes(ray1->n_exit2bn[2],bn1,settings1);
3194 anita1->rx_minarrivaltime=Tools::WhichIsMin(anita1->arrival_times[0], settings1->NANTENNAS);
3197 for (
int i=0;i<settings1->NANTENNAS;i++) {
3205 max_antenna_volts0 =0;
3206 max_antenna_volts1 =0;
3207 max_antenna_volts2 =0;
3220 if (settings1->SLAC)
3221 fslac_hitangles << bn1->sslacpositions[bn1->islacposition] <<
"\n";
3224 double rotateangle=getRNG(RNG_RANDOMISE_POL)->Gaus(RANDOMISEPOL*RADDEG);
3225 n_pol=n_pol.Rotate(rotateangle, ray1->n_exit2bn[2]);
3228 if (bn1->WHICHPATH==4) {
3229 Tools::Zero(sumsignal, 5);
3230 for (
int k=0;k<Anita::NFREQ;k++)
3231 IntegrateBands(anita1, k, panel1, anita1->freq, bn1->r_bn.Distance(interaction1->posnu)/1.E6, sumsignal);
3235 bn1->CenterPayload(hitangle_e);
3238 if (settings1->MAKEVERTICAL) {
3242 Vector rotationaxis=ray1->n_exit2bn[2].Cross(bn1->n_bn);
3243 double rotateangle=PI/2.-ray1->n_exit2bn[2].Dot(bn1->n_bn);
3244 ray1->n_exit2bn[2]=ray1->n_exit2bn[2].Rotate(rotateangle, rotationaxis);
3246 for (
int ilayer=0; ilayer < settings1->NLAYERS; ilayer++) {
3248 for (
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
3249 Vector rotationaxis2=ray1->n_exit2bn_eachboresight[2][ilayer][ifold].Cross(n_pol_eachboresight[ilayer][ifold]);
3250 double rotateangle2=PI/2.-ray1->n_exit2bn_eachboresight[2][ilayer][ifold].Dot(n_pol_eachboresight[ilayer][ifold]);
3251 ray1->n_exit2bn_eachboresight[2][ilayer][ifold].Rotate(rotateangle2, rotationaxis2);
3256 globaltrig1->volts_rx_rfcm_trigger.assign(16, vector <vector <double> >(3, vector <double>(0)));
3257 anita1->rms_rfcm_e_single_event = 0;
3259 if (!settings1->BORESIGHTS) {
3260 bn1->GetEcompHcompkvector(n_eplane, n_hplane, n_normal, ray1->n_exit2bn[2], e_component_kvector[0], h_component_kvector[0], n_component_kvector[0]);
3261 bn1->GetEcompHcompEvector(settings1, n_eplane, n_hplane, n_pol, e_component[0], h_component[0], n_component[0]);
3264 for (
int ilayer=0; ilayer < settings1->NLAYERS; ilayer++) {
3265 for (
int ifold=0;ifold<anita1->NRX_PHI[ilayer];ifold++) {
3270 bn1->GetAntennaOrientation(settings1, anita1, ilayer, ifold, n_eplane, n_hplane, n_normal);
3272 if (settings1->BORESIGHTS){
3273 bn1->GetEcompHcompkvector(n_eplane, n_hplane, n_normal, ray1->n_exit2bn_eachboresight[2][ilayer][ifold], e_component_kvector[count_rx], h_component_kvector[count_rx], n_component_kvector[count_rx]);
3274 bn1->GetEcompHcompEvector(settings1, n_eplane, n_hplane, n_pol_eachboresight[ilayer][ifold], e_component[count_rx], h_component[count_rx], n_component[count_rx]);
3275 if (settings1->SLAC) fslac_hitangles << ilayer <<
"\t" << ifold <<
"\t" << hitangle_e <<
"\t" << hitangle_h <<
"\t" << e_component_kvector <<
"\t" << h_component_kvector <<
"\t" << fresnel1_eachboresight[ilayer][ifold] <<
" " << mag1_eachboresight[ilayer][ifold] <<
"\n";
3277 else if (count_rx > 0)
3279 e_component_kvector[count_rx] = e_component_kvector[0];
3280 h_component_kvector[count_rx] = h_component_kvector[0];
3281 n_component_kvector[count_rx] = n_component_kvector[0];
3282 e_component[count_rx] = e_component[0];
3283 h_component[count_rx] = h_component[0];
3284 n_component[count_rx] = n_component[0];
3287 bn1->GetHitAngles(e_component_kvector[count_rx], h_component_kvector[count_rx], n_component_kvector[count_rx], hitangle_e, hitangle_h);
3289 hitangle_h_all[count_rx]=hitangle_h;
3290 hitangle_e_all[count_rx]=hitangle_e;
3292 if (h6->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3293 h6->Fill(hitangle_h, ray1->n_exit2bn[2]*bn1->n_bn);
3295 antNum = anita1->GetRxTriggerNumbering(ilayer, ifold);
3297 chantrig1->
ApplyAntennaGain(settings1, anita1, bn1, panel1, antNum, n_eplane, n_hplane, n_normal);
3299 chantrig1->
TriggerPath(settings1, anita1, antNum, bn1);
3327 chantrig1->
TimeShiftAndSignalFluct(settings1, anita1, ilayer, ifold, volts_rx_rfcm_lab_e_all, volts_rx_rfcm_lab_h_all);
3329 chantrig1->
saveTriggerWaveforms(anita1, justSignal_trig[0][antNum], justSignal_trig[1][antNum], justNoise_trig[0][antNum], justNoise_trig[1][antNum]);
3330 chantrig1->
saveDigitizerWaveforms(anita1, justSignal_dig[0][antNum], justSignal_dig[1][antNum], justNoise_dig[0][antNum], justNoise_dig[1][antNum]);
3332 Tools::Zero(sumsignal, 5);
3334 if (bn1->WHICHPATH==4 && ilayer==anita1->GetLayer(anita1->rx_minarrivaltime) && ifold==anita1->GetIfold(anita1->rx_minarrivaltime)) {
3335 for (
int ibw=0;ibw<5;ibw++) {
3336 cout <<
"Just after Taper, sumsignal is " << sumsignal_aftertaper[ibw] <<
"\n";
3337 cout <<
"Just after antennagain, sumsignal is " << sumsignal[ibw] <<
"\n";
3342 if (count_rx==anita1->rx_minarrivaltime) {
3343 undogaintoheight_e/=(double)Anita::NFREQ;
3344 undogaintoheight_h/=(double)Anita::NFREQ;
3345 for (
int k=0;k<4;k++) {
3346 undogaintoheight_e_array[k]/=(double)nbins_array[k];
3347 undogaintoheight_h_array[k]/=(double)nbins_array[k];
3350 if (settings1->SCALEDOWNLCPRX1)
3351 globaltrig1->volts[0][ilayer][0]=globaltrig1->volts[0][ilayer][0]/sqrt(2.);
3353 if (settings1->RCPRX2ZERO)
3354 globaltrig1->volts[1][ilayer][1]=0.;
3356 if (settings1->LCPRX2ZERO)
3357 globaltrig1->volts[0][ilayer][1]=0.;
3359 if (settings1->SIGNAL_FLUCT) {
3360 if (settings1->WHICH==0) {
3361 globaltrig1->volts[ilayer][ifold][0]+=getRNG(RNG_SIGNAL_FLUCT)->Gaus(0., anita1->VNOISE_ANITALITE[ifold]);
3362 globaltrig1->volts[ilayer][ifold][1]+=getRNG(RNG_SIGNAL_FLUCT)->Gaus(0., anita1->VNOISE_ANITALITE[ifold]);
3365 if (count_rx==anita1->rx_minarrivaltime) {
3366 rec_efield=sqrt(pow(globaltrig1->volts_original[0][ilayer][ifold]/(undogaintoheight_e*0.5), 2)+pow(globaltrig1->volts_original[1][ilayer][ifold]/(undogaintoheight_h*0.5), 2));
3367 for (
int ibw=0;ibw<4;ibw++) {
3368 rec_efield_array[ibw]=sqrt(pow(chantrig1->
bwslice_volts_pole[ibw]/(undogaintoheight_e_array[ibw]*0.5), 2)+pow(chantrig1->
bwslice_volts_polh[ibw]/(undogaintoheight_h_array[ibw]*0.5), 2));
3369 bwslice_vnoise_thislayer[ibw]=anita1->bwslice_vnoise[ilayer][ibw];
3372 if (tree6b->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3378 chantrig1->
WhichBandsPass(settings1, anita1, globaltrig1, bn1, ilayer, ifold, viewangle-sig1->changle, emfrac, hadfrac, thresholdsAnt[antNum]);
3381 if (Anita::GetAntennaNumber(ilayer, ifold)==anita1->rx_minarrivaltime) {
3382 for (
int iband=0;iband<5;iband++) {
3383 for (
int ipol=0;ipol<2;ipol++) {
3384 rx0_signal_eachband[ipol][iband]=chantrig1->
signal_eachband[ipol][iband];
3386 rx0_noise_eachband[ipol][iband]=chantrig1->
noise_eachband[ipol][iband];
3387 rx0_passes_eachband[ipol][iband]=chantrig1->
passes_eachband[ipol][iband];
3393 if (ilayer == 0 && globaltrig1->volts[0][ilayer][ifold] > max_antenna_volts0) {
3394 max_antenna0 = count_rx;
3395 max_antenna_volts0 = globaltrig1->volts[0][ilayer][ifold];
3396 max_antenna_volts0_em=globaltrig1->volts_em[0][ilayer][ifold];
3397 ant_max_normal0 = ant_normal;
3398 e_comp_max1 = e_component[count_rx];
3399 h_comp_max1 = h_component[count_rx];
3401 else if (ilayer == 0 && globaltrig1->volts[0][ilayer][ifold] == max_antenna_volts0 && globaltrig1->volts[0][ilayer][ifold] != 0){
3402 cout<<
"Equal voltage on two antennas! Event : "<<inu<<endl;
3404 else if (ilayer == 1 && globaltrig1->volts[0][ilayer][ifold] > max_antenna_volts1) {
3405 max_antenna1 = count_rx;
3406 max_antenna_volts1 = globaltrig1->volts[0][ilayer][ifold];
3407 ant_max_normal1 = ant_normal;
3408 e_comp_max2 = e_component[count_rx];
3409 h_comp_max2 = h_component[count_rx];
3411 else if (ilayer == 1 && globaltrig1->volts[0][ilayer][ifold] == max_antenna_volts1 && globaltrig1->volts[0][ilayer][ifold] != 0){
3412 cout<<
"Equal voltage on two antennas! Event : "<<inu<<endl;
3414 else if (ilayer == 2 && globaltrig1->volts[0][ilayer][ifold] > max_antenna_volts2) {
3415 max_antenna2 = count_rx;
3416 max_antenna_volts2 = globaltrig1->volts[0][ilayer][ifold];
3417 ant_max_normal2 = ant_normal;
3418 e_comp_max3 = e_component[count_rx];
3419 h_comp_max3 = h_component[count_rx];
3421 else if (ilayer == 2 && globaltrig1->volts[0][ilayer][ifold] == max_antenna_volts2 && globaltrig1->volts[0][ilayer][ifold] != 0){
3422 cout<<
"Equal voltage on two antennas! Event : "<<inu<<endl;
3424 voltagearray[count_rx] = globaltrig1->volts[0][ilayer][ifold];
3430 if (settings1->TRIGTYPE==0 && ifold==1 && count_pass>=settings1->NFOLD) {
3431 al_voltages_direct<<
"0 0 0"<<
" "<<
" "<<globaltrig1->volts_original[1][0][0]<<
" "<<(globaltrig1->volts_original[0][0][0]/sqrt(2.))<<
" "<<globaltrig1->volts_original[1][0][1]<<
" "<<globaltrig1->volts_original[0][0][1]<<
" "<<anita1->VNOISE[0]<<
" "<<anita1->VNOISE[0]<<
" "<<anita1->VNOISE[0]<<
" "<<anita1->VNOISE[0]<<
" "<<weight<<endl;
3438 anita1->rms_rfcm_e_single_event = sqrt(anita1->rms_rfcm_e_single_event / (anita1->HALFNFOUR * settings1->NANTENNAS));
3440 if(!settings1->ROUGHNESS){
3441 if (settings1->DISCONES==1) {
3443 for (
int idiscone=0;NDISCONES;idiscone++) {
3446 polarfactor_discone=n_pol.Dot(bn1->n_bn);
3447 for (
int k=0;k<Anita::NFREQ;k++) {
3448 if (anita1->freq[k]>=FREQ_LOW_DISCONES && anita1->freq[k]<=FREQ_HIGH_DISCONES) {
3449 thislambda=CLIGHT/sig1->N_AIR/anita1->freq[k];
3450 heff_discone= thislambda*sqrt(2*Zr*gain_dipole/Z0/4/PI*sig1->N_AIR);
3452 volts_discone+=panel1->GetVmmhz_freq(k)*0.5*heff_discone*((settings1->BW/1E6)/(
double)Anita::NFREQ)*polarfactor_discone;
3456 vnoise_discone=anita1->VNOISE[0]*sqrt(BW_DISCONES/settings1->BW_SEAVEYS);
3458 if (settings1->SIGNAL_FLUCT) {
3459 volts_discone+=getRNG(RNG_SIGNAL_FLUCT)->Gaus(0., vnoise_discone);
3462 if (fabs(volts_discone)/vnoise_discone>anita1->maxthreshold)
3469 for (
int irx=0;irx<settings1->NANTENNAS;irx++) {
3470 nchannels_perrx_triggered[irx]=globaltrig1->nchannels_perrx_triggered[irx];
3473 nchannels_triggered=Tools::iSum(globaltrig1->nchannels_perrx_triggered, settings1->NANTENNAS);
3474 volts_rx_ave=GetAverageVoltageFromAntennasHit(settings1, globaltrig1->nchannels_perrx_triggered, voltagearray, volts_rx_sum);
3479 if(sec1->secondbang && sec1->interestedintaus)
3480 count_asktrigger_nfb++;
3483 dist_int_bn_2d_chord = ray1->rfexit[0].
Distance(bn1->r_bn_shadow)/1000;
3485 dist_int_bn_2d = ray1->rfexit[0].
SurfaceDistance(bn1->r_bn_shadow, bn1->surface_under_balloon) / 1000;
3501 if(tauweighttrigger==1) {
3513 weight = weight1 / interaction1->
dnutries * settings1->SIGMA_FACTOR * time_weight;
3514 weight_prob = weight_prob / interaction1->
dnutries * settings1->SIGMA_FACTOR * time_weight;
3522 #ifdef ANITA3_EVENTREADER 3523 if (truthNuPtr) truthNuPtr->setWeights(weight, 1./interaction1->
dnutries, time_weight);
3525 havent_set_weights =
false;
3528 if (weight < settings1->CUTONWEIGHTS || weight_prob < settings1->CUTONWEIGHTPROBS) {
3533 eventsfound_beforetrigger+=weight;
3545 globaltrig1->
PassesTrigger(settings1, anita1, discones_passing, 2, l3trignoise, l2trig, l1trig, settings1->antennaclump, loctrig, loctrig_nadironly, inu,
3546 thispassesnoise,
true);
3548 globaltrig1->
PassesTrigger(settings1, anita1, discones_passing, 2, l3trig, l2trig, l1trig, settings1->antennaclump, loctrig, loctrig_nadironly, inu,
3551 if ( (l3trignoise[0]>0 ) || (l3trignoise[1]>0 ) ){
3552 cout <<
"A thermal noise fluctuation generated this trigger!" << l3trignoise[0] <<
" " << l3trig[0] <<
" " << l3trignoise[1] <<
" " << l3trig[1] << endl;
3558 for (
int i=0;i<2;i++) {
3559 for (
int j=0;j<16;j++) {
3560 for (
int k=0;k<anita1->HALFNFOUR;k++) {
3561 count1->nl1triggers[i][whichray]+=anita1->l1trig_anita3and4_inanita[i][j][k];
3575 if ( (thispasses[0]==1 && anita1->pol_allowed[0]==1)
3576 || (thispasses[1]==1 && anita1->pol_allowed[1]==1)
3577 || (settings1->TRIGTYPE==0 && count_pass>=settings1->NFOLD)
3578 || (settings1->MINBIAS==1)){
3580 if (bn1->WHICHPATH==4)
3581 cout <<
"This event passes.\n";
3583 anita1->passglobtrig[0]=thispasses[0];
3584 anita1->passglobtrig[1]=thispasses[1];
3587 n_exit_phi = Tools::AbbyPhiCalc(ray1->n_exit2bn[2][0], ray1->n_exit2bn[2][1]);
3591 count1->npassestrigger[whichray]++;
3596 if (tree11->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3600 if(sec1->secondbang && sec1->interestedintaus)
3601 count_passestrigger_nfb++;
3603 interaction1->crust_entered=0;
3604 interaction1->mantle_entered=0;
3605 interaction1->core_entered=0;
3611 if (nupathtree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3615 count_chordgoodlength++;
3617 pieceofkm2sr=weight*antarctica->volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr/(double)NNU/len_int;
3618 if (h10->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST)
3619 h10->Fill(hitangle_e_all[0], weight);
3630 logweight=log10(weight);
3634 if (interaction1->
d2>1) {
3636 count_d2goodlength++;
3640 if(sec1->secondbang && sec1->interestedintaus) {
3641 eventsfound_nfb+=weight;
3642 index_weights=(int)(((logweight-MIN_LOGWEIGHT)/(MAX_LOGWEIGHT-MIN_LOGWEIGHT))*(double)NBINS);
3643 eventsfound_nfb_binned[index_weights]++;
3644 if (tree16->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3648 allcuts[whichray]++;
3649 allcuts_weighted[whichray]+=weight;
3650 if (thispasses[0] && thispasses[1]) {
3651 allcuts_weighted_polarization[2]+=weight;
3652 }
else if (thispasses[0]){
3653 allcuts_weighted_polarization[0]+=weight;
3654 }
else if (thispasses[1]){
3655 allcuts_weighted_polarization[1]+=weight;
3657 anita1->weight_inanita=weight;
3659 if (h1mybeta->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3660 h1mybeta -> Fill(mybeta, weight);
3662 eventsfound+=weight;
3663 eventsfound_prob+=weight_prob;
3664 eventsfound_divided_by_lint+=weight/len_int;
3666 if (NNU_per_source.size())
3668 eventsfound_per_source[which_source]+=weight*src_time_weight/time_weight;
3669 eventsfound_prob_per_source[which_source]+=weight_prob*src_time_weight/time_weight;
3674 eventsfound_belowhorizon+=weight;
3676 count1->npass[whichray]++;
3680 if (logweight<MIN_LOGWEIGHT)
3682 else if (logweight>MAX_LOGWEIGHT)
3683 index_weights=NBINS-1;
3685 index_weights=(int)(((logweight-MIN_LOGWEIGHT)/(MAX_LOGWEIGHT-MIN_LOGWEIGHT))*(double)NBINS);
3688 if (index_weights<NBINS)
3689 eventsfound_binned[index_weights]++;
3692 if (index_distance<NBINS_DISTANCE)
3693 eventsfound_binned_distance[index_distance]+= weight;
3696 if (index_distance<NBINS_DISTANCE && index_weights<NBINS)
3697 eventsfound_binned_distance_forerror[index_distance][index_weights]++;
3700 eventsfound_weightgt01+=weight;
3703 if (nearthlayers==1)
3704 eventsfound_crust+=weight;
3706 if (h1mybeta->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1) {
3707 h1mybeta -> Fill(mybeta, weight);
3708 h1mytheta -> Fill(mytheta, weight);
3714 double int_lon, int_lat;
3717 int_lon = interaction1->posnu.
Lon();
3718 int_lat = interaction1->posnu.
Lat();
3719 antarctica->LonLattoEN(int_lon, int_lat, E, N);
3721 dir_int_coord->Fill(E, N);
3723 ref_int_coord->Fill(E, N);
3727 offaxis=(double)fabs(viewangle-sig1->changle);
3728 nsigma_offaxis=offaxis/deltheta_had_max;
3730 hundogaintoheight_e->Fill(undogaintoheight_e, weight);
3731 hundogaintoheight_h->Fill(undogaintoheight_h, weight);
3732 rec_diff->Fill((rec_efield-true_efield)/true_efield, weight);
3733 rec_diff0->Fill((rec_efield_array[0]-true_efield_array[0])/true_efield_array[0], weight);
3734 rec_diff1->Fill((rec_efield_array[1]-true_efield_array[1])/true_efield_array[1], weight);
3735 rec_diff2->Fill((rec_efield_array[2]-true_efield_array[2])/true_efield_array[2], weight);
3736 rec_diff3->Fill((rec_efield_array[3]-true_efield_array[3])/true_efield_array[3], weight);
3737 recsum_diff->Fill((rec_efield_array[0]+rec_efield_array[1]+rec_efield_array[2]+rec_efield_array[3]-true_efield)/true_efield, weight);
3740 sourceLon = ray1->rfexit[2].
Lon() - 180;
3741 sourceLat = ray1->rfexit[2].
Lat() - 90;
3742 sourceAlt = antarctica->SurfaceAboveGeoid(sourceLon+180, sourceLat+90);
3745 if (settings1->HIST && finaltree->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3746 for (
int i=0;i<3;i++) {
3747 nnu_array[i] = interaction1->
nnu[i];
3748 r_in_array[i] = interaction1->
r_in[i];
3749 r_bn_array[i] = bn1->r_bn[i];
3750 n_bn_array[i] = bn1->n_bn[i];
3751 posnu_array[i] = interaction1->posnu[i];
3752 ant_max_normal0_array[i] = ant_max_normal0[i];
3753 ant_max_normal1_array[i] = ant_max_normal1[i];
3754 ant_max_normal2_array[i] = ant_max_normal2[i];
3755 n_pol_array[i] = n_pol[i];
3756 r_enterice_array[i] = interaction1->
r_enterice[i];
3757 nsurf_rfexit_array[i] = ray1->nsurf_rfexit[i];
3758 nsurf_rfexit_db_array[i] = ray1->nsurf_rfexit_db[i];
3760 for (
int j=0;j<5;j++) {
3761 for (
int i=0;i<3;i++) {
3762 nrf_iceside_array[j][i] = ray1->nrf_iceside[j][i];
3763 nrf_iceside_db_array[j][i] = nrf_iceside_db[j][i];
3764 n_exit2bn_array[j][i] = ray1->n_exit2bn[j][i];
3765 n_exit2bn_db_array[j][i] = n_exit2bn_db[j][i];
3766 rfexit_array[j][i] = ray1->rfexit[j][i];
3767 rfexit_db_array[j][i] = ray1->rfexit_db[j][i];
3770 if (settings1->HIST && vmmhz_tree->GetEntries()<20) {
3779 d12=interaction1->
d1;
3780 d22=interaction1->
d2;
3790 sourceMag = ray1->rfexit[2].Mag();
3791 sample_x = antarctica->getSampleX();
3792 sample_y = antarctica->getSampleY();
3795 count1->IncrementWeights_r_in(interaction1->
r_in, weight);
3798 #ifdef ANITA_UTIL_EXISTS 3802 Adu5PatPtr->
latitude= bn1->latitude;
3804 Adu5PatPtr->
altitude=bn1->altitude;
3805 Adu5PatPtr->
realTime=bn1->realTime_flightdata;
3806 Adu5PatPtr->
heading = bn1->heading;
3807 Adu5PatPtr->
pitch = bn1->pitch;
3808 Adu5PatPtr->roll = bn1->roll;
3809 Adu5PatPtr->run = run_no;
3812 memset(realEvPtr->
fVolts, 0,
sizeof(realEvPtr->
fVolts) );
3813 memset(realEvPtr->
fTimes, 0,
sizeof(realEvPtr->
fTimes) );
3815 int fNumPoints = 260;
3816 for (
int ichan=0; ichan<108; ichan++){
3819 for (
int j = 0; j < fNumPoints; j++) {
3821 realEvPtr->
fTimes[ichan][j] = j * anita1->TIMESTEP * 1.0E9;
3825 for (
int iant = 0; iant < settings1->NANTENNAS; iant++){
3827 int IceMCAnt = GetIceMCAntfromUsefulEventAnt(settings1, iant);
3832 realEvPtr->
chanId[UsefulChanIndexV] = UsefulChanIndexV;
3833 realEvPtr->
chanId[UsefulChanIndexH] = UsefulChanIndexH;
3835 for (
int j = 0; j < fNumPoints; j++) {
3840 realEvPtr->
fVolts[UsefulChanIndexH][j] = volts_rx_rfcm_lab_h_all[IceMCAnt][j+128]*1000;
3842 realEvPtr->
fVolts[UsefulChanIndexV][j] = volts_rx_rfcm_lab_e_all[IceMCAnt][j+128]*1000;
3853 if (settings1->MINBIAS==1)
3859 rawHeaderPtr->
run = run_no;
3869 if (settings1->WHICH<9){
3870 rawHeaderPtr->
phiTrigMask = (short) anita1->phiTrigMask;
3875 rawHeaderPtr->
realTime = bn1->realTime_flightdata;
3876 rawHeaderPtr->
triggerTime = bn1->realTime_flightdata;
3877 Adu5PatPtr->
latitude= bn1->latitude;
3879 Adu5PatPtr->
altitude=bn1->altitude;
3880 Adu5PatPtr->
realTime=bn1->realTime_flightdata;
3881 Adu5PatPtr->
heading = bn1->heading;
3882 Adu5PatPtr->
pitch = bn1->pitch;
3883 Adu5PatPtr->roll = bn1->roll;
3884 Adu5PatPtr->run = run_no;
3887 cout <<
"Neutrino (evNum = " << eventNumber <<
", weight = " << weight <<
" weight_prob = " << weight_prob <<
") passed" << endl;
3889 #ifdef ANITA3_EVENTREADER 3890 if (settings1->WHICH==9 || settings1->WHICH==10) {
3893 rawHeaderPtr->setMask( (
short) anita1->l1TrigMask, (
short) anita1->phiTrigMask,
AnitaPol::kVertical);
3894 rawHeaderPtr->setMask( (
short) anita1->l1TrigMaskH, (
short) anita1->phiTrigMaskH,
AnitaPol::kHorizontal);
3899 truthEvPtr->
realTime = bn1->realTime_flightdata;
3900 truthEvPtr->
run = run_no;
3901 truthEvPtr->
nuMom = pnu;
3902 truthEvPtr->
showerE = pnu * sumfrac;
3903 truthEvPtr->
nu_pdg = pdgcode;
3914 if(src_model){cout <<
"It originated from " << objName <<
" (" << RA*15 <<
"°, " << dec <<
"°)" << endl;}
3916 truthEvPtr->
weight = weight;
3917 truthEvPtr->
weight1 = weight1;
3922 truthEvPtr->
l_int = len_int;
3923 truthEvPtr->
tuffIndex = (short)anita1->tuffIndex;
3927 #ifdef ANITA3_EVENTCORRELATOR
3928 if(settings1->ALL_SKY_MAP)
3930 astroObject->SetX(RA);
3931 astroObject->SetY(dec);
3932 astroObject->SetMarkerStyle(20);
3933 astroObject->SetMarkerSize(1.5);
3934 astroObject->SetMarkerColor(9);
3935 skyMapOut->addMarker(astroObject);
3941 truthEvPtr->RA = RA;
3942 truthEvPtr->
dec = dec;
3947 for (
int i=0;i<3;i++){
3950 truthEvPtr->
nuPos[i] = interaction1->posnu[i];
3951 truthEvPtr->
nuDir[i] = interaction1->
nnu[i];
3953 for (
int i=0;i<5;i++){
3954 for (
int j=0;j<3;j++){
3955 truthEvPtr->
rfExitNor[i][j] = ray1->n_exit2bn[i][j];
3956 truthEvPtr->
rfExitPos[i][j] = ray1->rfexit[i][j];
3959 for (
int i=0;i<48;i++){
3960 truthEvPtr->
hitangle_e[i] = hitangle_e_all[i];
3961 truthEvPtr->
hitangle_h[i] = hitangle_h_all[i];
3963 if(!settings1->ROUGHNESS){
3964 for (
int i=0;i<Anita::NFREQ;i++)
3965 truthEvPtr->
vmmhz[i] = panel1->GetVmmhz_freq(i);
3981 for (
int iant = 0; iant < settings1->NANTENNAS; iant++){
3985 truthEvPtr->
SNRAtTrigger[UsefulChanIndexV] = Tools::calculateSNR(justSignal_trig[0][iant], justNoise_trig[0][iant]);
3986 truthEvPtr->
SNRAtTrigger[UsefulChanIndexH] = Tools::calculateSNR(justSignal_trig[1][iant], justNoise_trig[1][iant]);
3991 truthEvPtr->
SNRAtDigitizer[UsefulChanIndexV] = Tools::calculateSNR(justSignal_dig[0][iant], justNoise_dig[0][iant]);
3992 truthEvPtr->
SNRAtDigitizer[UsefulChanIndexH] = Tools::calculateSNR(justSignal_dig[1][iant], justNoise_dig[1][iant]);
3998 truthEvPtr->
thresholds[UsefulChanIndexV] = thresholdsAnt[iant][0][4];
3999 truthEvPtr->
thresholds[UsefulChanIndexH] = thresholdsAnt[iant][1][4];
4002 if (iant%2) irx = iant/2;
4003 else irx = iant/2 + 1;
4006 if (settings1->WRITE_WAVEFORMS)
4008 for (
int j = 0; j < fNumPoints; j++) {
4009 truthEvPtr->
fTimes[UsefulChanIndexV][j] = j * anita1->TIMESTEP * 1.0E9;
4010 truthEvPtr->
fTimes[UsefulChanIndexH][j] = j * anita1->TIMESTEP * 1.0E9;
4012 truthEvPtr->
fSignalAtTrigger[UsefulChanIndexV][j] = justSignal_trig[0][iant][j+128]*1000;
4013 truthEvPtr->
fSignalAtTrigger[UsefulChanIndexH][j] = justSignal_trig[1][iant][j+128]*1000;
4014 truthEvPtr->
fNoiseAtTrigger[UsefulChanIndexV][j] = justNoise_trig[0][iant][j+128]*1000;
4015 truthEvPtr->
fNoiseAtTrigger[UsefulChanIndexH][j] = justNoise_trig[1][iant][j+128]*1000;
4016 truthEvPtr->
fSignalAtDigitizer[UsefulChanIndexV][j] = justSignal_dig[0][iant][j+128]*1000;
4017 truthEvPtr->
fSignalAtDigitizer[UsefulChanIndexH][j] = justSignal_dig[1][iant][j+128]*1000;
4018 truthEvPtr->
fNoiseAtDigitizer[UsefulChanIndexV][j] = justNoise_dig[0][iant][j+128]*1000;
4019 truthEvPtr->
fNoiseAtDigitizer[UsefulChanIndexH][j] = justNoise_dig[1][iant][j+128]*1000;
4021 truthEvPtr->
fDiodeOutput[UsefulChanIndexV][j] = anita1->timedomain_output_allantennas[0][irx][j];
4022 truthEvPtr->
fDiodeOutput[UsefulChanIndexH][j] = anita1->timedomain_output_allantennas[1][irx][j];
4028 if (truthNuPtr) truthNuPtr->setSkipped(
false);
4032 if (truthAnitaNuTree) truthAnitaNuTree->Fill();
4033 truthAnitaTree->Fill();
4039 adu5PatTree->Fill();
4042 delete rawHeaderPtr;
4046 sum_weights+=weight;
4047 neutrinos_passing_all_cuts++;
4048 times_crust_entered_det+=interaction1->crust_entered;
4049 times_mantle_entered_det+=interaction1->mantle_entered;
4050 times_core_entered_det+=interaction1->core_entered;
4052 if (settings1->WRITEPOSFILE==1)
4053 WriteNeutrinoInfo(interaction1->posnu, interaction1->
nnu, bn1->r_bn, interaction1->
altitude_int, interaction1->nuflavor(), interaction1->current(), elast_y, nu_out);
4056 if (settings1->HIST && !settings1->ONLYFINAL && sampleweights->GetEntries()<settings1->HIST_MAX_ENTRIES) {
4058 sampleweights->Fill(log10(weight));
4060 sampleweights->Fill(-6.);
4063 if (sampleweights->GetEntries()==1000) {
4064 double sum_sampleintegral=0.;
4065 double sum_sample=0.;
4067 for (
int k=sampleweights->GetNbinsX();k>=1;k--) {
4068 sum_sampleintegral+=sampleweights->GetBinContent(k)*pow(10., sampleweights->GetBinLowEdge(k));
4071 sum_sampleintegral+=sampleweights->GetBinContent(0)*pow(10., sampleweights->GetBinLowEdge(1));
4073 for (
int k=sampleweights->GetNbinsX();k>=1;k--) {
4074 sum_sample+=sampleweights->GetBinContent(k)*pow(10., sampleweights->GetBinLowEdge(k));
4075 if (sum_sample>0.99*sum_sampleintegral) {
4088 forbrian << interaction1->
costheta_nutraject <<
" " << n_nutraject_ontheground.Phi() <<
" " << bn1->phi_bn <<
" " << logweight <<
"\n";
4091 if (interaction1->nuflavor()==
"nue") {
4093 sum_prob[0]+=weight_prob;
4094 eventsfound_binned_e[index_weights]++;
4096 if (interaction1->nuflavor()==
"numu") {
4098 sum_prob[1]+=weight_prob;
4099 eventsfound_binned_mu[index_weights]++;
4101 if(!sec1->secondbang || !sec1->interestedintaus) {
4102 if (interaction1->nuflavor()==
"nutau") {
4104 sum_prob[2]+=weight_prob;
4105 eventsfound_binned_tau[index_weights]++;
4113 cout <<
"Chord is less than 1m.\n";
4116 if (settings1->HIST==1 && !settings1->ONLYFINAL && anita1->tglob->GetEntries()<settings1->HIST_MAX_ENTRIES) {
4118 anita1->tglob->Fill();
4119 anita1->tdata->Fill();
4120 anita1->tgaryanderic->Fill();
4133 if (bn1->WHICHPATH==4)
4134 cout <<
"This event does not pass.\n";
4171 count_chanceinhell2_w += weight;
4174 count_passestrigger_w += weight;
4176 volume_thishorizon=antarctica->volume_inhorizon[bn1->Getibnposition()]/1.E9;
4178 if (settings1->HIST==1
4179 && !settings1->ONLYFINAL
4180 && tree1->GetEntries()<settings1->HIST_MAX_ENTRIES
4181 && bn1->WHICHPATH != 3){
4188 std::cout <<
"\n***********************************************************";
4189 std::cout <<
"\n* SIGINT received, aborting loop over events early.";
4190 std::cout <<
"\n* Stopped after event " << inu <<
" instead of " << NNU;
4191 std::cout <<
"\n* Any output which relied on NNU should be corrected for.";
4192 std::cout <<
"\n***********************************************************\n";
4193 foutput <<
"\n***********************************************************";
4194 foutput <<
"\n* SIGINT received, aborting loop over events early.";
4195 foutput <<
"\n* Stopped after event " << inu <<
" instead of " << NNU;
4196 foutput <<
"\n* Any output which relied on NNU should be corrected for.";
4197 foutput <<
"\n***********************************************************\n";
4205 cout <<
"about to close tsignals tree.\n";
4206 anita1->fsignals=anita1->tsignals->GetCurrentFile();
4207 anita1->fdata=anita1->tdata->GetCurrentFile();
4208 anita1->fdata=anita1->tgaryanderic->GetCurrentFile();
4209 anita1->fsignals->Write();
4210 anita1->fsignals->Close();
4212 anita1->fdata=anita1->tglob->GetCurrentFile();
4213 anita1->fdata->Write();
4214 anita1->fdata->Close();
4216 if (settings1->EVENTSMAP){
4218 TH2F *lat80deg=
new TH2F(
"lat80deg",
"", 600, -3000, 3000, 500, -2500, 2500);
4219 lat80deg->SetMarkerColor(kRed);
4221 for(
double lon=0;lon<360.;lon+=0.5){
4223 antarctica->LonLattoEN(lon, lat, E, N);
4224 lat80deg->Fill(E, N);
4228 if (bn1->WHICHPATH==4) {
4229 cout <<
"Earth radius at South Pole: " << antarctica->Geoid(0.) <<
"\n";
4231 cout <<
"Surface of ice at the South Pole: " << antarctica->Surface(posnu_temp) <<
"\n";
4232 cout <<
"Average balloon altitude is " << average_altitude <<
"\n";
4233 cout <<
"Average distance from earth center is " << average_rbn <<
"\n";
4234 cout <<
"Average height of balloon above ice surface is " << average_rbn-antarctica->Surface(bn1->r_bn) <<
"\n";
4235 cout <<
"theta_zenith are " << anita1->THETA_ZENITH[0]*DEGRAD <<
" " << anita1->THETA_ZENITH[1]*DEGRAD <<
" " << anita1->THETA_ZENITH[2]*DEGRAD <<
"\n";
4236 cout <<
"Index of refraction at this depth is " << sig1->N_DEPTH <<
"\n";
4237 cout <<
"Cerenkov angle is " << sig1->changle*DEGRAD <<
"\n";
4238 cout <<
"Nadir angle to surface exit point is " << DEGRAD*bn1->r_bn.Angle(ray1->n_exit2bn[2]) <<
"\n";
4239 cout <<
"Distance from rfexit to balloon is " << (bn1->r_bn+ -1*ray1->rfexit[2]).Mag() <<
"\n";
4240 cout <<
"Payload zenith angle at event source is " << DEGRAD*ray1->rfexit[2].Angle(ray1->n_exit2bn[2]) <<
"\n";
4241 cout <<
"Angle of incidence just below surface is " << DEGRAD*ray1->rfexit[2].Angle(ray1->nrf_iceside[4]) <<
"\n";
4242 cout <<
"Angle between neutrino and surface normal is " << DEGRAD*ray1->rfexit[0].Angle(interaction1->
nnu)-90. <<
"\n";
4243 cout <<
"Angle of incidence below firn surface is " << DEGRAD*ray1->rfexit[2].Angle(ray1->nrf_iceside[3]) <<
"\n";
4245 cout <<
"about to Summarize.\n";
4247 anita1->rms_rfcm[0] = sqrt(anita1->rms_rfcm[0] / (
double)anita1->count_getnoisewaveforms)*1000.;
4248 anita1->rms_rfcm[1] = sqrt(anita1->rms_rfcm[1] / (
double)anita1->count_getnoisewaveforms)*1000.;
4249 anita1->rms_lab[0] = sqrt(anita1->rms_lab[0] / (
double)anita1->count_getnoisewaveforms)*1000.;
4250 anita1->rms_lab[1] = sqrt(anita1->rms_lab[1] / (
double)anita1->count_getnoisewaveforms)*1000.;
4252 cout <<
"RMS noise in rfcm e-pol is " << anita1->rms_rfcm[0] <<
" mV.\n";
4253 cout <<
"RMS noise in rfcm h-pol is " << anita1->rms_rfcm[1] <<
" mV.\n";
4254 cout <<
"RMS noise in lab e-pol is " << anita1->rms_lab[0] <<
"mV.\n";
4255 cout <<
"RMS noise in lab h-pol is " << anita1->rms_lab[1] <<
"mV.\n";
4256 for (
int i=0;i<Anita::NFREQ;i++) {
4257 anita1->avgfreq_rfcm[i]/=(double)anita1->count_getnoisewaveforms;
4258 anita1->avgfreq_rfcm_lab[i]/=(
double)anita1->count_getnoisewaveforms;
4261 rms_rfcm_e=anita1->rms_rfcm[0];
4262 rms_rfcm_h=anita1->rms_rfcm[1];
4263 rms_lab_e=anita1->rms_lab[0];
4264 rms_lab_h=anita1->rms_lab[1];
4265 for (
int i=0;i<Anita::NFREQ;i++) {
4266 avgfreq_rfcm[i]=anita1->avgfreq_rfcm[i];
4267 avgfreq_rfcm_lab[i]=anita1->avgfreq_rfcm_lab[i];
4268 freq[i]=anita1->freq[i];
4270 cout <<
"Filling summarytree. rms_rfcm_e is " << rms_rfcm_e <<
"\n";
4271 summarytree->Fill();
4275 primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, xsecParam_nutype, Interaction::ktotal);
4276 len_int=1.0/(sigma*sig1->RHOH20*(1./M_NUCL)*1000);
4279 Summarize(settings1, anita1, count1, spectra1, sig1, primary1, pnu, eventsfound, eventsfound_db, eventsfound_nfb, sigma, sum, antarctica->volume, antarctica->ice_area, km3sr, km3sr_e, km3sr_mu, km3sr_tau, foutput, distanceout, outputdir, bn1, src_model);
4281 std::string spec_string = src_model ? settings1->SOURCE : std::to_string( settings1->EXPONENT);
4282 veff_out << spec_string <<
"\t" << km3sr <<
"\t" << km3sr_e <<
"\t" << km3sr_mu <<
"\t" << km3sr_tau <<
"\t" << settings1->SIGMA_FACTOR << endl;
4285 sum_frac[0]=sum[0]/eventsfound;
4286 sum_frac[1]=sum[1]/eventsfound;
4287 sum_frac[2]=sum[2]/eventsfound;
4290 sum_frac_db[0]=sum[0]/(eventsfound+eventsfound_db+eventsfound_nfb);
4291 sum_frac_db[1]=sum[1]/(eventsfound+eventsfound_db+eventsfound_nfb);
4292 sum_frac_db[2]=(sum[2]+eventsfound_db+eventsfound_nfb)/(eventsfound+eventsfound_db+eventsfound_nfb);
4299 #ifdef ANITA_UTIL_EXISTS 4301 anitafileEvent->cd();
4302 eventTree->Write(
"eventTree");
4303 anitafileEvent->Close();
4304 delete anitafileEvent;
4306 anitafileHead->cd();
4307 headTree->Write(
"headTree");
4308 anitafileHead->Close();
4309 delete anitafileHead;
4312 adu5PatTree->Write(
"adu5PatTree");
4313 anitafileGps->Close();
4314 delete anitafileGps;
4316 #ifdef ANITA3_EVENTREADER 4317 anitafileTruth->cd();
4318 configAnitaTree->Write(
"configAnitaTree");
4319 truthAnitaTree->Write(
"truthAnitaTree");
4320 if (truthAnitaNuTree) truthAnitaNuTree->Write();
4321 triggerSettingsTree->Write(
"triggerSettingsTree");
4322 summaryAnitaTree->Fill();
4323 summaryAnitaTree->Write(
"summaryAnitaTree");
4324 anitafileTruth->Close();
4325 delete anitafileTruth;
4331 #ifdef ANITA3_EVENTCORRELATOR 4332 if(settings1->ALL_SKY_MAP)
4334 TCanvas *skyMapC =
new TCanvas(
"skyMapC",
"skyMapC", 1000, 500);
4337 string outputSkyMap =string(outputdir.Data())+
"/allSkyMap"+run_num+
".png";
4338 skyMapC->SaveAs(outputSkyMap.c_str());
4342 cout <<
"closing file.\n";
4345 time_t raw_end_time = time(NULL);
4346 struct tm * end_time = localtime(&raw_end_time);
4347 cout <<
"Date and time at end of run are: " << asctime (end_time) <<
"\n";
4348 cout<<
"Total time elapsed is "<<(int)((raw_end_time - raw_start_time)/60)<<
":"<< ((raw_end_time - raw_start_time)%60)<<endl;
4350 foutput <<
"\nTotal time elapsed in run is " <<(int)((raw_end_time - raw_start_time)/60)<<
":"<< ((raw_end_time - raw_start_time)%60)<<endl;
4376 void IntegrateBands(Anita *anita1,
int k, Screen *panel1,
double *freq,
double scalefactor,
double *sumsignal) {
4377 for (
int j=0;j<5;j++) {
4379 for (
int jpt=0; jpt<panel1->GetNvalidPoints(); jpt++){
4380 if (anita1->bwslice_min[j]<=freq[k] && anita1->bwslice_max[j]>freq[k])
4381 sumsignal[j]+=panel1->GetVmmhz_freq(jpt*Anita::NFREQ + k)*(freq[k+1]-freq[k])*scalefactor;
4388 void Integrate(Anita *anita1,
int j,
int k,
double *vmmhz,
double *freq,
double scalefactor,
double sumsignal) {
4390 if (anita1->bwslice_min[j]<=freq[k] && anita1->bwslice_max[j]>freq[k])
4391 sumsignal+=vmmhz[k]*(freq[k+1]-freq[k])*scalefactor;
4396 void WriteNeutrinoInfo(
Position &posnu,
Vector &nnu,
Position &r_bn,
double altitude_int,
string nuflavor,
string current,
double elast_y, ofstream &nu_out) {
4397 nu_out <<
"\n" << inu <<
"\t" << posnu[0] <<
" " << posnu[1] <<
" " << posnu[2] <<
"\t" << altitude_int <<
"\t" << nnu[0] <<
" " << nnu[1] <<
" " << nnu[2] <<
"\t" << r_bn[0] <<
" " << r_bn[1] <<
" " << r_bn[2] <<
"\t" << nuflavor <<
"\t" << current <<
"\t" << elast_y <<
"\n\n";
4402 void Summarize(
Settings *settings1, Anita* anita1,
Counting *count1,
Spectra *spectra1,
Signal *sig1,
Primaries *primary1,
double pnu,
double eventsfound,
double eventsfound_db,
double eventsfound_nfb,
double sigma,
double* sum,
double volume,
double ice_area,
double& km3sr,
double& km3sr_e,
double& km3sr_mu,
double& km3sr_tau, ofstream &foutput, ofstream &distanceout, TString outputdir,
Balloon * bn1,
SourceModel * src_model) {
4403 double rate_v_thresh[NTHRESHOLDS];
4404 double errorup_v_thresh[NTHRESHOLDS];
4405 double errordown_v_thresh[NTHRESHOLDS];
4406 double rate_h_thresh[NTHRESHOLDS];
4407 double errorup_h_thresh[NTHRESHOLDS];
4408 double errordown_h_thresh[NTHRESHOLDS];
4409 double zeroes[NTHRESHOLDS];
4410 Tools::Zero(zeroes, NTHRESHOLDS);
4413 for (
int i=0;i<NTHRESHOLDS;i++) {
4414 rate_v_thresh[i]=npass_v_thresh[i]/denom_v_thresh[i];
4416 if (npass_v_thresh[i]<=20) {
4417 errorup_v_thresh[i]=poissonerror_plus[(int)npass_v_thresh[i]]/denom_v_thresh[i];
4418 errordown_v_thresh[i]=poissonerror_minus[(int)npass_v_thresh[i]]/denom_v_thresh[i];
4421 if (settings1->WHICH==9 || settings1->WHICH==10){
4422 rate_h_thresh[i]=npass_h_thresh[i]/denom_h_thresh[i];
4423 if (npass_h_thresh[i]<=20) {
4424 errorup_h_thresh[i]=poissonerror_plus[(int)npass_h_thresh[i]]/denom_h_thresh[i];
4425 errordown_h_thresh[i]=poissonerror_minus[(int)npass_h_thresh[i]]/denom_h_thresh[i];
4431 string stemp=string(outputdir.Data())+
"/thresholds.root";
4432 TFile *fthresholds=
new TFile(stemp.c_str(),
"RECREATE");
4433 TCanvas *cthresh=
new TCanvas(
"cthresh",
"cthresh", 880, 800);
4436 TGraph *gnpass=
new TGraph(NTHRESHOLDS, thresholds, npass_v_thresh);
4437 gnpass->SetName(
"npass");
4438 TGraph *gdenom=
new TGraph(NTHRESHOLDS, thresholds, denom_v_thresh);
4439 gdenom->SetName(
"denom");
4441 TGraphAsymmErrors *g=
new TGraphAsymmErrors(NTHRESHOLDS, thresholds, rate_v_thresh, zeroes, zeroes, errorup_v_thresh, errordown_v_thresh);
4445 g->SetMarkerStyle(21);
4448 stemp = string(outputdir.Data())+
"/thresholds.eps";
4449 cthresh->Print((TString)stemp);
4455 if (settings1->WHICH==9 || settings1->WHICH==10){
4456 TGraph *gnpassH=
new TGraph(NTHRESHOLDS, thresholds, npass_h_thresh);
4457 gnpassH->SetName(
"npassH");
4458 TGraph *gdenomH=
new TGraph(NTHRESHOLDS, thresholds, denom_h_thresh);
4459 gdenomH->SetName(
"denomH");
4460 TGraphAsymmErrors *gH=
new TGraphAsymmErrors(NTHRESHOLDS, thresholds, rate_h_thresh, zeroes, zeroes, errorup_h_thresh, errordown_h_thresh);
4461 gH->SetName(
"rate");
4462 gH->SetLineWidth(2);
4463 gH->SetMarkerStyle(21);
4465 stemp = string(outputdir.Data())+
"/thresholds_HPOL.eps";
4466 cthresh->Print((TString)stemp);
4473 fthresholds->Write();
4474 fthresholds->Close();
4481 foutput <<
"Generated " << NNU <<
" neutrinos with energy " << pnu <<
"\n";
4482 foutput <<
"Number of unweighted direct, reflected that pass is: " << count1->npass[0] <<
"\t" << count1->npass[1] <<
"\n";
4483 foutput <<
"Number of (weighted) neutrinos that pass is: " << eventsfound <<
"\n";
4484 foutput <<
"Number of (weighted) neutrinos that pass, multiplied by prob. of interacting in the ice, is: " << eventsfound_prob <<
"\n";
4485 foutput <<
"Number of (weighted) neutrinos that pass, divided by lint is: " << eventsfound_divided_by_lint <<
"\n";
4486 foutput <<
"Number of weighted direct, reflected that pass is: " << allcuts_weighted[0] <<
"\t" << allcuts_weighted[1] <<
"\n";
4487 foutput <<
"Number of (weighted) neutrinos that pass (with weight>0.001) is: " << eventsfound_weightgt01 <<
"\n";
4488 foutput <<
"Number of (weighted) neutrinos that only traverse the crust is " << eventsfound_crust <<
" -> " << eventsfound_crust/eventsfound*100 <<
"%\n\n";
4489 foutput <<
"Number of (weighted) neutrinos that pass only VPOL trigger is: " << allcuts_weighted_polarization[0] <<
"\n";
4490 foutput <<
"Number of (weighted) neutrinos that pass only HPOL trigger is: " << allcuts_weighted_polarization[1] <<
"\n";
4491 foutput <<
"Number of (weighted) neutrinos that pass both pol triggers is: " << allcuts_weighted_polarization[2] <<
"\n\n";
4493 cout <<
"Number of (weighted) neutrinos that pass (with weight>0.001) is: " << eventsfound_weightgt01 <<
"\n";
4494 cout <<
"Number of (weighted) neutrinos that pass, multiplied by prob. of interacting in the ice, is: " << eventsfound_prob <<
"\n";
4495 cout <<
"Number of (weighted) neutrinos that pass, divided by lint is: " << eventsfound_divided_by_lint <<
"\n";
4496 cout <<
"Number of (weighted) neutrinos that only traverse the crust is " << eventsfound_crust <<
" -> " << eventsfound_crust/eventsfound*100 <<
"%\n\n";
4497 cout <<
"Number of (weighted) neutrinos that pass only VPOL trigger is: " << allcuts_weighted_polarization[0] <<
"\n";
4498 cout <<
"Number of (weighted) neutrinos that pass only HPOL trigger is: " << allcuts_weighted_polarization[1] <<
"\n";
4499 cout <<
"Number of (weighted) neutrinos that pass both pol triggers is: " << allcuts_weighted_polarization[2] <<
"\n\n";
4501 foutput <<
"Total number of thrown electron neutrinos is : " << (double)count1->nnu_e <<
"\n";
4502 foutput <<
"Total number of thrown muon neutrinos is : " << (
double)count1->nnu_mu <<
"\n";
4503 foutput <<
"Total number of thrown tau neutrinos is : " << (double)count1->nnu_tau <<
"\n";
4504 foutput <<
"Number of (weighted) electron neutrinos that pass trigger is : " << sum[0] <<
"\n";
4505 foutput <<
"Number of (weighted) muon neutrinos that pass trigger is : " << sum[1] <<
"\n";
4506 foutput <<
"Number of (weighted) tau neutrinos that pass trigger is : " << sum[2] <<
"\n\n";
4508 cout <<
"Total number of thrown electron neutrinos is : " << (
double)count1->nnu_e <<
"\n";
4509 cout <<
"Total number of thrown muon neutrinos is : " << (double)count1->nnu_mu <<
"\n";
4510 cout <<
"Total number of thrown tau neutrinos is : " << (
double)count1->nnu_tau <<
"\n";
4511 cout <<
"Number of (weighted) electron neutrinos that pass trigger is : " << sum[0] <<
"\n";
4512 cout <<
"Number of (weighted) muon neutrinos that pass trigger is : " << sum[1] <<
"\n";
4513 cout <<
"Number of (weighted) tau neutrinos that pass trigger is : " << sum[2] <<
"\n\n";
4515 foutput <<
"Volume of ice is " << volume <<
"\n";
4516 foutput <<
"Value of 4*pi*pi*r_earth*r_earth in km^2 " << 4*PI*PI*(EarthModel::R_EARTH*EarthModel::R_EARTH/1.E6) <<
"\n";
4532 double nevents_db=0;
4533 double nevents_nfb=0;
4536 nevents+=eventsfound;
4537 nevents_db+=eventsfound_db;
4538 nevents_nfb+=eventsfound_nfb;
4550 error_percent_increase_nfb=0;
4554 for (
int i=0;i<NBINS;i++) {
4555 if (eventsfound_binned[i]<=20) {
4556 error_plus+=pow(poissonerror_plus[(
int)eventsfound_binned[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4557 error_e_plus+=pow(poissonerror_plus[(
int)eventsfound_binned_e[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4558 error_mu_plus+=pow(poissonerror_plus[(
int)eventsfound_binned_mu[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4559 error_tau_plus+=pow(poissonerror_plus[(
int)eventsfound_binned_tau[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4560 error_minus+=pow(poissonerror_minus[(
int)eventsfound_binned[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4561 error_e_minus+=pow(poissonerror_minus[(
int)eventsfound_binned_e[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4562 error_mu_minus+=pow(poissonerror_minus[(
int)eventsfound_binned_mu[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4563 error_tau_minus+=pow(poissonerror_minus[(
int)eventsfound_binned_tau[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4566 error_plus+=eventsfound_binned[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4567 error_e_plus+=eventsfound_binned_e[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4568 error_mu_plus+=eventsfound_binned_mu[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4569 error_tau_plus+=eventsfound_binned_tau[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4570 error_minus=error_plus;
4571 error_e_minus=error_e_plus;
4572 error_mu_minus=error_mu_plus;
4573 error_tau_minus=error_tau_plus;
4576 error_nfb+=eventsfound_nfb_binned[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*4.+-5.), 2);
4577 for (
int j=0;j<NBINS_DISTANCE;j++) {
4578 if (eventsfound_binned_distance_forerror[j][i]<=20) {
4579 error_distance_plus[j]+=pow(poissonerror_plus[(
int)eventsfound_binned_distance_forerror[j][i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4580 error_distance_minus[j]+=pow(poissonerror_minus[(
int)eventsfound_binned_distance_forerror[j][i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4586 error_plus=sqrt(error_plus);
4587 error_e_plus=sqrt(error_e_plus);
4588 error_mu_plus=sqrt(error_mu_plus);
4589 error_tau_plus=sqrt(error_tau_plus);
4591 error_minus=sqrt(error_minus);
4592 error_e_minus=sqrt(error_e_minus);
4593 error_mu_minus=sqrt(error_mu_minus);
4594 error_tau_minus=sqrt(error_tau_minus);
4597 error_nfb=sqrt(error_nfb);
4598 for (
int j=0;j<NBINS_DISTANCE;j++) {
4599 error_distance_plus[j]=sqrt(error_distance_plus[j]);
4600 error_distance_minus[j]=sqrt(error_distance_minus[j]);
4605 if (NNU != 0 && nevents!=0) {
4606 km3sr=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*nevents/(double)NNU/settings1->SIGMA_FACTOR;
4608 km2sr = volume*pow(1.e-3,3)*eventsfound_divided_by_lint*sr/double(NNU)/settings1->SIGMA_FACTOR;
4610 cout << nevents <<
" events passed out of " << NNU <<
"\n";
4612 error_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_plus/(double)NNU/settings1->SIGMA_FACTOR;
4613 error_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_minus/(double)NNU/settings1->SIGMA_FACTOR;
4615 percent_increase_db=km3sr_db/km3sr*100;
4616 percent_increase_nfb=km3sr_nfb/km3sr*100;
4618 error_percent_increase_nfb=sqrt(pow(error_km3sr_nfb/km3sr_nfb, 2)+pow(error_plus/km3sr, 2))*percent_increase_nfb;
4620 percent_increase_total=percent_increase_db+percent_increase_nfb;
4622 km3sr_e = (pow(1.e-3, 3))*volume*sr*(sum[0]/(
double)count1->nnu_e)*sig1->RHOMEDIUM/sig1->RHOH20/settings1->SIGMA_FACTOR;
4624 cout << sum[0]/(
double)nevents*100. <<
"% are electron neutrinos\n";
4626 error_e_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_e_plus/(double)count1->nnu_e/settings1->SIGMA_FACTOR;
4627 error_e_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_e_minus/(double)count1->nnu_e/settings1->SIGMA_FACTOR;
4629 km3sr_mu = (pow(1.e-3, 3))*volume*sr*(sum[1]/(
double)count1->nnu_mu)*sig1->RHOMEDIUM/sig1->RHOH20/settings1->SIGMA_FACTOR;
4631 cout << sum[1]/(
double)nevents*100. <<
"% are muon neutrinos\n";
4633 error_mu_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_mu_plus/(double)count1->nnu_mu/settings1->SIGMA_FACTOR;
4634 error_mu_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_mu_minus/(double)count1->nnu_mu/settings1->SIGMA_FACTOR;
4636 km3sr_tau = (pow(1.e-3, 3))*volume*sr*(sum[2]/(
double)count1->nnu_tau)*sig1->RHOMEDIUM/sig1->RHOH20/settings1->SIGMA_FACTOR;
4638 cout << sum[2]/(
double)nevents*100. <<
"% are tau neutrinos\n";
4640 error_tau_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_tau_plus/(double)count1->nnu_tau/settings1->SIGMA_FACTOR;
4641 error_tau_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_tau_minus/(double)count1->nnu_tau/settings1->SIGMA_FACTOR;
4644 double sum_km3sr_error_plus=0;
4645 double sum_km3sr_error_minus=0;
4646 for (
int i=0;i<NBINS_DISTANCE;i++) {
4647 sum_km3sr+= (pow(1.e-3, 3))*volume*sr*eventsfound_binned_distance[i]*sig1->RHOMEDIUM/sig1->RHOH20/(
double)NNU/settings1->SIGMA_FACTOR;
4648 km3sr_distance[i]=sum_km3sr;
4649 sum_km3sr_error_plus+=pow(pow(1.E-3, 3)*volume*error_distance_plus[i]*sig1->RHOMEDIUM/sig1->RHOH20*sr/(
double)NNU, 2)/settings1->SIGMA_FACTOR;
4650 error_distance_plus[i]=sqrt(sum_km3sr_error_plus);
4651 sum_km3sr_error_minus+=pow(pow(1.E-3, 3)*volume*error_distance_minus[i]*sig1->RHOMEDIUM/sig1->RHOH20*sr/(
double)NNU, 2)/settings1->SIGMA_FACTOR;
4652 error_distance_minus[i]=sqrt(sum_km3sr_error_minus);
4653 distanceout << 700000./(double)NBINS_DISTANCE*(
double)i <<
"\t" << km3sr_distance[i] <<
"\t" << error_distance_plus[i] <<
"\t" << error_distance_minus[i] <<
"\n";
4656 foutput <<
"Total volume * solid angle is \t\t\t\t" << km3sr <<
" + " << error_plus <<
" - " << error_minus <<
" km^3 str\n";
4657 foutput <<
"Total volume * solid angle for electron neutrinos is \t" << km3sr_e <<
" + " << error_e_plus <<
" - " << error_e_minus <<
" km^3 str\n";
4658 foutput <<
"Total volume * solid angle for muon neutrinos is \t" << km3sr_mu <<
" + " << error_mu_plus <<
" - " << error_mu_minus <<
" km^3 str\n";
4659 foutput <<
"Total volume * solid angle for tau neutrinos is \t" << km3sr_tau <<
" + " << error_tau_plus <<
" - " << error_tau_minus <<
" km^3 str\n";
4660 cout <<
"Total volume * solid angle is \t\t\t\t" << km3sr <<
" + " << error_plus <<
" - " << error_minus <<
" km^3 str\n";
4661 cout <<
"Total volume * solid angle for electron neutrinos is \t" << km3sr_e <<
" + " << error_e_plus <<
" - " << error_e_minus <<
" km^3 str\n";
4662 cout <<
"Total volume * solid angle for muon neutrinos is \t" << km3sr_mu <<
" + " << error_mu_plus <<
" - " << error_mu_minus <<
" km^3 str\n";
4663 cout <<
"Total volume * solid angle for tau neutrinos is \t" << km3sr_tau <<
" + " << error_tau_plus <<
" - " << error_tau_minus <<
" km^3 str\n";
4665 if ( spectra1->IsMonoenergetic() ) {
4666 cout <<
"Cross section is " << sigma <<
"m^2\n";
4667 double len_int=1.0/(sigma*sig1->RHOH20*(1./M_NUCL)*1000);
4668 cout <<
"Interaction length is " << len_int <<
" km\n\n";
4669 foutput <<
"Interaction length is " << len_int <<
"\n";
4670 double mean_km2sr=km3sr/len_int;
4671 foutput <<
"Total area x steradians using km3sr/ CClen_int is \t\t\t\t" << mean_km2sr <<
" km^2 str\n\n";
4672 foutput <<
"Total area x steradians using sum(dkm3sr/lint) is \t\t\t\t" <<km2sr <<
" km^2 str\n\n";
4673 cout <<
"Total area x steradians is \t\t" << km2sr <<
" km^2 str\n\n";
4674 double alt_km2sr=ice_area/(1.E6)*PI*eventsfound_prob/(
double)NNU;
4675 foutput <<
"Total area x steradians using 4*PI*R_EARTH^2*eff. is \t" << alt_km2sr <<
" km^2 str\n\n";
4676 foutput <<
"These are not the same because we are not throwing all directions on all points of the surface. Believe the first one as an approximation, we are working on this for high cross sections.\n";
4683 foutput <<
"\t\t\t\t\t\t\tprobability of passing \t# events passing\n";
4685 foutput.precision(4);
4686 foutput <<
"No way this neutrino will see any ice \t\t\t" << (double)count1->noway[0]/(
double)count_total <<
"\t" <<
4687 (double)count1->noway[1]/(
double)count_total <<
"\t\t" <<
4688 count1->noway[0] <<
"\t" << count1->noway[1] <<
"\n";
4690 foutput.precision(4);
4691 foutput <<
"Wheredoesitleave in PickUnbiased gives an error \t" << (double)count1->wheredoesitleave_err[0]/(
double)count1->noway[0] <<
"\t" <<
4692 (double)count1->wheredoesitleave_err[1]/(
double)count1->noway[1] <<
"\t\t" <<
4693 count1->wheredoesitleave_err[0] <<
"\t" << count1->wheredoesitleave_err[1] <<
"\n";
4695 foutput.precision(4);
4696 foutput <<
"This neutrino direction never sees ice \t\t\t" << (double)count1->neverseesice[0]/(
double)count1->wheredoesitleave_err[0] <<
"\t" <<
4697 (double)count1->neverseesice[1]/(
double)count1->wheredoesitleave_err[1] <<
"\t\t" <<
4698 count1->neverseesice[0] <<
"\t" << count1->neverseesice[1] <<
"\n";
4701 foutput.precision(4);
4702 foutput <<
"WhereDoesItEnterIce in PickUnbiased gives an error \t\t\t" << (double)count1->wheredoesitenterice_err[0]/(
double)count1->neverseesice[0] <<
"\t" <<
4703 (double)count1->wheredoesitenterice_err[1]/(
double)count1->neverseesice[1] <<
"\t\t" <<
4704 count1->wheredoesitenterice_err[0] <<
"\t" << count1->wheredoesitenterice_err[1] <<
"\n";
4706 foutput.precision(4);
4707 foutput <<
"Interaction point too high \t\t\t\t" << (double)count1->toohigh[0]/(
double)count1->wheredoesitenterice_err[0] <<
"\t" <<
4708 (double)count1->toohigh[1]/(
double)count1->wheredoesitenterice_err[1] <<
"\t\t" <<
4709 count1->toohigh[0] <<
"\t" << count1->toohigh[1] <<
"\n";
4711 foutput.precision(4);
4712 foutput <<
"Interaction point too low \t\t\t\t" << (double)count1->toolow[0]/(
double)count1->toohigh[0] <<
"\t" <<
4713 (double)count1->toolow[1]/(
double)count1->toohigh[1] <<
"\t\t" <<
4714 count1->toolow[0] <<
"\t" << count1->toolow[1] <<
"\n";
4719 foutput.precision(4);
4720 foutput <<
"There is an interaction in ice \t\t\t\t" << (double)count1->iceinteraction[0]/(
double)count1->toolow[0] <<
"\t" <<
4721 (double)count1->iceinteraction[1]/(
double)count1->toolow[1] <<
"\t\t" <<
4722 count1->iceinteraction[0] <<
"\t" << count1->iceinteraction[1] <<
"\n";
4724 foutput.precision(4);
4725 foutput <<
"In horizon \t\t\t\t\t\t" << (double)count1->inhorizon[0]/(
double)count1->iceinteraction[0] <<
"\t" <<
4726 (double)count1->inhorizon[1]/(
double)count1->iceinteraction[1] <<
"\t\t" <<
4727 count1->inhorizon[0] <<
"\t" << count1->inhorizon[1] <<
"\n";
4731 foutput.precision(4);
4732 foutput <<
"From surface to balloon, ray not intersected by earth \t" << (double)count1->nraypointsup1[0]/(
double)count1->inhorizon[0] <<
"\t" <<
4733 (double)count1->nraypointsup1[1]/(
double)count1->inhorizon[1] <<
"\t\t" <<
4734 count1->nraypointsup1[0] <<
"\t" << count1->nraypointsup1[1] <<
"\n";
4736 foutput.precision(4);
4737 foutput <<
"After 1/r scaling and best case attenuation, \n\tMaximum signal is detectable\t\t\t" << (double)count1->nnottoosmall[0]/(
double)count1->nraypointsup1[0] <<
"\t" <<
4738 (double)count1->nnottoosmall[1]/(
double)count1->nraypointsup1[1] <<
"\t\t" 4739 << count1->nnottoosmall[0] <<
"\t" << count1->nnottoosmall[1] <<
"\n";
4742 foutput.precision(4);
4743 foutput <<
"Viewing angle lt 90 degrees\t\t\t" << (double)count1->nviewangle_lt_90[0]/(
double)count1->nnottoosmall[0] <<
"\t" <<
4744 (double)count1->nviewangle_lt_90[1]/(
double)count1->nnottoosmall[1] <<
"\t\t" 4745 << count1->nviewangle_lt_90[0] <<
"\t" << count1->nviewangle_lt_90[1] <<
"\n";
4748 foutput.precision(4);
4749 foutput <<
"Reality check: EM and hadronic fractions both nonzero\t" << (double)count1->ngoodfracs[0]/(
double)count1->nviewangle_lt_90[0] <<
"\t" <<
4750 (double)count1->ngoodfracs[1]/(
double)count1->nviewangle_lt_90[1] <<
"\t\t" 4751 << count1->ngoodfracs[0] <<
"\t" << count1->ngoodfracs[1] <<
"\n";
4752 foutput.precision(4);
4753 foutput <<
"\tBoth EM and hadronic fractions are zero\t\t" << (double)count1->nbadfracs[0]/(
double)count1->nviewangle_lt_90[0] <<
"\t" <<
4754 (double)count1->nbadfracs[1]/(
double)count1->nviewangle_lt_90[1] <<
"\t\t" <<
4755 count1->nbadfracs[0] <<
"\t" << count1->nbadfracs[1] <<
"\n";
4757 foutput.precision(4);
4758 foutput <<
"After finding neutrino direction, \n\tchance of making through Earth\t\t\t" << (double)count_chanceofsurviving/(
double)count1->ngoodfracs[0] << "\t\t\t";
4759 foutput.precision(10);
4760 foutput << count_chanceofsurviving << "\n";
4761 foutput.precision(4);
4762 foutput << "Neutrino enters ice south of 60deg S latitude\t\t" << (
double)count1->nentersice[0]/(
double)count_chanceofsurviving << "\t" <<
4763 (
double)count1->nentersice[1]/(
double)count_chanceofsurviving <<
4765 count1->nentersice[0] << "\t" << count1->nentersice[1] << "\n";
4767 foutput.precision(4);
4768 foutput << "Neutrino reasonably likely
to survive trip through Earth " << (
double)count1->nabsorbed[0]/(
double)count1->nentersice[0] << "\t" <<
4769 (
double)count1->nabsorbed[1]/(
double)count1->nentersice[1] << "\t\t"
4770 << count1->nabsorbed[0] << "\t" << count1->nabsorbed[1] << "\n";
4771 foutput.precision(4);
4772 foutput << "
Ray leaves the ice south of 60deg S latitude\t\t" << (
double)count1->nraywithincontinent1[0]/(
double)count1->nabsorbed[0] << "\t" <<
4773 (
double)count1->nraywithincontinent1[1]/(
double)count1->nabsorbed[1] << "\t" <<
4774 count1->nraywithincontinent1[0] << "\t" <<
4775 count1->nraywithincontinent1[1] << "\n";
4777 foutput.precision(4);
4778 foutput << "After 1/r, best guess ice attenuation, \n\tmaximum signal is detectable\t\t\t" << (
double)count_chanceinhell0/(
double)count1->nraywithincontinent1[0] << "\t\t\t";
4779 foutput.precision(10);
4780 foutput <<count_chanceinhell0 << "\n";
4782 foutput.precision(4);
4783 foutput << "
Ray is not totally internally reflected\t\t\t" << (
double)count1->nnottir[0]/(
double)count_chanceinhell0 << "\t" <<
4784 (
double)count1->nnottir[1]/(
double)count_chanceinhell0 << "\t\t" <<
4785 count1->nnottir[0] << "\t" << count1->nnottir[1] << "\n";
4787 foutput.precision(4);
4788 foutput << "From surface
to balloon, ray not intersected by earth \t" << (
double)count1->nraypointsup2[0]/(
double)count1->nnottir[0] << "\t" <<
4789 (
double)count1->nraypointsup2[1]/(
double)count1->nnottir[1] <<
4791 << count1->nraypointsup2[0] << "\t" << count1->nraypointsup2[1] << "\n";
4792 foutput.precision(4);
4794 foutput << "
Ray leaves the ice south of 60deg S latitude\t\t" << (
double)count1->nraywithincontinent2[0]/(
double)count1->nraypointsup2[0] <<"\t" <<
4795 (
double)count1->nraywithincontinent2[0]/(
double)count1->nraypointsup2[1] <<
4796 "\t\t" << count1->nraywithincontinent2[0] << "\t" <<
4797 count1->nraywithincontinent2[1] << "\n";
4798 foutput.precision(4);
4800 foutput << "
Ray leaves where there is ice\t\t\t\t" << (
double)count1->nacceptablerf[0]/(
double)count1->nraywithincontinent2[0] << "\t" <<
4801 (
double)count1->nacceptablerf[1]/(
double)count1->nraywithincontinent2[1] << "\t\t"
4802 << count1->nacceptablerf[0] << "\t" << count1->nacceptablerf[1] << "\n";
4803 foutput.precision(4);
4805 foutput << "
Ray tracing converges
to within 10 m\t\t\t" << (
double)count1->nconverges[0]/(
double)count1->nacceptablerf[0] << "\t" <<
4806 (
double)count1->nconverges[1]/(
double)count1->nacceptablerf[1] <<
4807 "\t\t" << count1->nconverges[0] << "\t" << count1->nconverges[1] << "\n";
4808 foutput.precision(4);
4810 foutput << "After fresnel coefficient, \n\tmaximum signal is detectable\t\t\t" << (
double)count1->nchanceinhell_fresnel[0]/(
double)count1->nconverges[0] << "\t" << (
double)count1->nchanceinhell_fresnel[1]/(
double)count1->nconverges[1] <<
4811 "\t\t" <<count1->nchanceinhell_fresnel[0] << "\t" << count1->nchanceinhell_fresnel[1] << "\n";
4812 foutput.precision(4);
4813 foutput << "After 1/r, \n\tmaximum signal is detectable\t\t\t" << (
double)count1->nchanceinhell_1overr[0]/(
double)count1->nchanceinhell_fresnel[0] << "\t" << (
double)count1->nchanceinhell_1overr[1]/(
double)count1->nchanceinhell_fresnel[1] << "\t\t" <<count1->nchanceinhell_1overr[0] << "\t" << count1->nchanceinhell_1overr[1] << "\n";
4814 foutput.precision(4);
4816 foutput << "After ice attenuation, \n\tmaximum signal is detectable\t\t\t" << (
double)count1->nchanceinhell[0]/(
double)count1->nchanceinhell_1overr[0] << "\t" <<
4817 (
double)count1->nchanceinhell[1]/(
double)count1->nchanceinhell_1overr[1] << "\t" <<count1->nchanceinhell[0] << "\t" << count1->nchanceinhell[1] << "\n";
4818 foutput.precision(4);
4820 foutput << "After viewing angle cut, \t\t\t\t" << (
double)count1->nviewanglecut[0]/(
double)count1->nchanceinhell[0] << "\t" << (
double)count1->nviewanglecut[1]/(
double)count1->nchanceinhell[1] << "\t\t" << count1->nviewanglecut[0] << " " << count1->nviewanglecut[1] << "\n";
4822 foutput.precision(4);
4823 foutput << "After factoring in off-Cerenkov cone tapering, \n\tMaximum signal is detectable\t\t\t" << (
double)count1->nchanceinhell2[0]/(
double)count1->nviewanglecut[0] << "\t" << (
double)count1->nchanceinhell2[1]/(
double)count1->nviewanglecut[1] << "\t\t" << count1->nchanceinhell2[0] << " " << count1->nchanceinhell2[1] << "\n";
4825 foutput << "Survive dead time \t\t\t\t\t" << (
double)count1->ndeadtime[0]/(
double)count1->nchanceinhell2[0] << "\t" << (
double)count1->ndeadtime[1]/(
double)count1->nchanceinhell2[1] << "\t\t" << (
double)count1->ndeadtime[0] << " " << count1->ndeadtime[1] << "\n";
4827 foutput << "Passes trigger\t\t\t\t\t\t" << (
double)count1->npassestrigger[0]/(
double)count1->ndeadtime[0] << "\t" << (
double)count1->npassestrigger[1]/(
double)count1->ndeadtime[1] << "\t\t" << count1->npassestrigger[0] << "\t" << count1->npassestrigger[1] << "\n";
4828 foutput << "Number of l1 triggers\t\t\t\t\t\t" << (
double)count1->nl1triggers[0][0] << "\t" << (
double)count1->nl1triggers[1][0] << "\n";
4830 foutput << "Chord is good length\t\t\t\t\t" << (
double)count_chordgoodlength/(
double)count1->npassestrigger[0] << "\t\t\t";
4831 foutput.precision(10);
4832 foutput <<count_chordgoodlength << "\n";
4833 foutput.precision(4);
4834 foutput << "Neutrino's path in ice more than 1m \t\t\t" << (
double)count_d2goodlength/(
double)count_chordgoodlength << "\t\t\t";
4835 foutput.precision(10);
4836 foutput << count_d2goodlength << "\n";
4837 foutput.precision(4);
4838 foutput << "Events that pass all cuts\t\t\t\t" << (
double)count1->npass[0]/(
double)count_d2goodlength << "\t" << (
double)count1->npass[1]/(
double)count_d2goodlength << "\t\t";
4839 foutput.precision(10);
4840 foutput <<count1->npass[0] << "\t" << count1->npass[1] << "\n";
4842 cout << "Events that pass all cuts\t\t\t\t" << (
double)count1->npass[0]/(
double)count_d2goodlength << "\t" << (
double)count1->npass[1]/(
double)count_d2goodlength << "\t\t";
4843 cout <<count1->npass[0] << "\t" << count1->npass[1] << "\n";
4849 TString str; str.Form(
"%s/source_times.txt", outputdir.Data());
4850 FILE * f = fopen(str.Data(),
"w");
4851 std::vector<double> times;
4852 src_model->computeFluxTimeChanges(×);
4853 for (
unsigned i = 0; i < times.size(); i++)
4855 fprintf(f,
"%f\n", times[i]);
4861 if ( src_model || spectra1->IsSpectrum() ) {
4862 double sum_events=0.;
4863 double thisenergy=0.;
4864 double thislen_int_kgm2=0.;
4872 double min_energy = src_model ? settings1->SOURCE_MIN_E : spectra1->Getenergy()[0];
4873 double max_energy = src_model ? settings1->SOURCE_MAX_E : spectra1->Getenergy()[spectra1->GetE_bin()-1];
4874 even_E = ( max_energy - min_energy ) / ( (
double) N_even_E );
4879 src_flux = src_model->estimateFlux(bn1->min_time,bn1->max_time, TMath::Power(10,min_energy-9), TMath::Power(10,max_energy-9), 2*N_even_E, 10000);
4885 for (
int i=0;i<N_even_E;i++) {
4886 thisenergy=pow(10., (min_energy+((
double)i)*even_E));
4887 primary1->
GetSigma(thisenergy, sigma, thislen_int_kgm2, settings1, xsecParam_nutype, Interaction::ktotal);
4890 double EdNdEdAdt = src_model ? 1e9*src_flux->Interpolate(thisenergy*1e-9) : spectra1->GetEdNdEdAdt(log10(thisenergy));
4898 sum_events+=even_E*log(10.)*( EdNdEdAdt*1e4 )/(thislen_int_kgm2/sig1->RHOH20);
4899 integral+=even_E*log(10.)*(EdNdEdAdt);
4900 cout <<
"thisenergy, EdNdEdAdt is " << thisenergy <<
" " << EdNdEdAdt <<
"\n";
4916 cout <<
"SUM EVENTS IS " << sum_events << endl;
4917 cout <<
"INTEGRAL : " << integral << endl;
4918 sum_events*=volume*anita1->LIVETIME*sig1->RHOMEDIUM/sig1->RHOH20*nevents/(double)NNU*sr;
4920 foutput <<
"volume, LIVETIME, sig1->RHOMEDIUM, RHOH20, nevents, NNU, sr are " << volume <<
" " << anita1->LIVETIME <<
" " << sig1->RHOMEDIUM <<
" " << sig1->RHOH20 <<
" " << nevents <<
" " << NNU <<
" " << sr <<
"\n";
4921 foutput <<
"Total events observed is " << sum_events <<
"\n";
4928 double GetAirDistance(
double altitude_bn,
double beta) {
4929 return EarthModel::R_EARTH*acos((altitude_bn+EarthModel::R_EARTH)/EarthModel::R_EARTH*(1-sin(beta)*sin(beta))+1/EarthModel::R_EARTH*sin(beta)*sqrt((altitude_bn+EarthModel::R_EARTH)*(altitude_bn+EarthModel::R_EARTH)*sin(beta)*sin(beta)-2*EarthModel::R_EARTH*altitude_bn-altitude_bn*altitude_bn));
4934 double GetAverageVoltageFromAntennasHit(
Settings *settings1,
int *nchannels_perrx_triggered,
double *voltagearray,
double& volts_rx_sum) {
4936 int count_hitantennas=0;
4937 for (
int i=0;i<settings1->NANTENNAS;i++) {
4938 if (nchannels_perrx_triggered[i]>=3) {
4939 sum+=voltagearray[i];
4940 count_hitantennas++;
4944 sum = sum/(double)count_hitantennas;
4955 Vector n_bfield = nnu.Cross(nrf2_iceside);
4957 Vector n_pol = n_bfield.Cross(nrf2_iceside);
4958 n_pol = n_pol.Unit();
4960 if (nnu*nrf2_iceside>0 && n_pol*nnu>0){
4961 cout <<
"error in GetPolarization. Event is " << inu <<
"\n";
4968 void CloseTFile(TFile *hfile) {
4976 double IsItDoubleBang(
double exitlength,
double plepton) {
4977 double gamma=plepton/MTAU;
4978 return 1-exp(-1*exitlength/(TAUDECAY_TIME*CLIGHT*gamma));
4986 double gamma=pnu/MTAU;
4988 if (exp(-1*nuexitlength/(TAUDECAY_TIME*CLIGHT*gamma))>0.999){
4989 rnd1=getRNG(RNG_SECOND_BANG)->Rndm()*nuexitlength;
4992 while (rnd2>1-exp(-1*rnd1/(TAUDECAY_TIME*CLIGHT*gamma))) {
4993 rnd1=getRNG(RNG_SECOND_BANG)->Rndm()*nuexitlength;
4994 rnd2=getRNG(RNG_SECOND_BANG)->Rndm();
4997 posnu2 = posnu + rnd1*nnu;
4998 rfexit_db = antarctica1->Surface(posnu2)*posnu2.Unit();
5001 n_exit2bn_db = (r_bn - rfexit_db) / r_bn.
Distance(rfexit_db);
5003 double cosangle=(n_exit2bn_db * posnu2) / posnu2.Mag();
5014 double ATTENLENGTH=700;
5015 if(!settings1->VARIABLE_ATTEN){
5016 ATTENLENGTH=antarctica1->EffectiveAttenuationLength(settings1, posnu, 1);
5019 int position_in_iceshelves=antarctica1->IceOnWater(posnu);
5020 int position_in_rossexcept=antarctica1->RossExcept(posnu);
5023 double dtemp=posnu_down.
Distance(rfexit2)/ATTENLENGTH;
5027 if(position_in_iceshelves && (!position_in_rossexcept)){
5030 scalefactor_attenuation=0.71*exp(-dtemp);
5031 vmmhz_max=vmmhz_max*0.71*exp(-dtemp);
5033 else if(position_in_rossexcept){
5034 scalefactor_attenuation=0.1*exp(-dtemp);
5035 vmmhz_max=0.1*vmmhz_max*exp(-dtemp);
5038 scalefactor_attenuation=sqrt(0.001)*exp(-dtemp);
5039 vmmhz_max=sqrt(0.001)*vmmhz_max*exp(-dtemp);
5043 scalefactor_attenuation=0;
5050 void Attenuate(
IceModel *antarctica1,
Settings *settings1,
double& vmmhz_max,
double rflength,
const Position &posnu) {
5051 double ATTENLENGTH=700;
5052 if (!settings1->VARIABLE_ATTEN){
5053 ATTENLENGTH = antarctica1->EffectiveAttenuationLength(settings1, posnu, 0);
5056 double dtemp=(rflength/ATTENLENGTH);
5057 if(!settings1->ROUGHNESS){
5059 scalefactor_attenuation=exp(-dtemp);
5060 vmmhz_max=vmmhz_max*exp(-dtemp);
5063 scalefactor_attenuation=0;
5069 scalefactor_attenuation=exp(-dtemp);
5070 vmmhz_max=vmmhz_max*exp(-dtemp);
5073 scalefactor_attenuation=0;
5081 void IsAbsorbed(
double chord_kgm2,
double len_int_kgm2,
double &weight1) {
5087 rtemp=chord_kgm2/len_int_kgm2;
5089 weight1=exp(-rtemp);
5096 void GetSmearedIncidentAngle(
Vector &specular,
Vector &nrf_iceside,
Vector &n_exit2bn,
double SMEARINCIDENTANGLE){
5099 specular+=nrf_iceside;
5100 Vector parallel_to_surface;
5101 parallel_to_surface+=n_exit2bn;
5102 parallel_to_surface.Cross(specular);
5103 nrf_iceside.Rotate(SMEARINCIDENTANGLE*(2*getRNG(RNG_SMEARED_INCIDENT_ANGLE)->Rndm()-1.), parallel_to_surface);
5109 int GetRayIceSide(
const Vector &n_exit2rx,
const Vector &nsurf_rfexit,
double nexit,
double nenter,
Vector &nrf2_iceside) {
5112 double NRATIO=nexit/nenter;
5113 costh=(n_exit2rx*nsurf_rfexit)/(n_exit2rx.Mag() * nsurf_rfexit.Mag());
5119 double sinth=sqrt(1 - costh*costh);
5120 double factor=NRATIO*costh-sqrt(1-(NRATIO*sinth*NRATIO*sinth));
5121 nrf2_iceside = -factor*nsurf_rfexit + NRATIO*n_exit2rx;
5122 nrf2_iceside = nrf2_iceside.Unit();
5128 int GetDirection(
Settings *settings1,
Interaction *interaction1,
const Vector &refr,
double deltheta_em,
double deltheta_had,
double emfrac,
double hadfrac,
double vmmhz1m_max,
double r_fromballoon,
Ray *ray1,
Signal *sig1,
Position posnu, Anita *anita1,
Balloon *bn1,
Vector &nnu,
double& costhetanu,
double& theta_threshold) {
5136 double theta_test=0;
5137 double vmmhz1m_test=0;
5138 double costhetanu1 = 0;
5139 double costhetanu2 = 0;
5143 theta_threshold = 0;
5148 if (settings1->SKIPCUTS || !settings1->USEDIRECTIONWEIGHTS) {
5154 if (emfrac<=1.E-10 && deltheta_had >1.E-10) {
5155 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(hadfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle)>1)
5160 theta_threshold=sqrt(-1*deltheta_had*deltheta_had*log(anita1->VNOISE[0]/10.*anita1->maxthreshold/(hadfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle))/ALOG2);
5161 averaging_thetas1+=theta_threshold;
5163 count_inthisloop1++;
5166 if (emfrac>1.E-10 && deltheta_had <=1.E-10) {
5168 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(emfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle)>1)
5173 theta_threshold=sqrt(-1*deltheta_em*deltheta_em*log(anita1->VNOISE[0]/10.*anita1->maxthreshold/(emfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle))/0.5);
5174 averaging_thetas2+=theta_threshold;
5176 count_inthisloop2++;
5181 if (emfrac>1.E-10 && deltheta_had>1.E-10) {
5184 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/((hadfrac+emfrac)*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.) {
5186 theta_threshold=-1.;
5189 theta_test=deltheta_em;
5190 vmmhz1m_test=vmmhz1m_max;
5192 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5194 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.) {
5195 theta_threshold=theta_test;
5198 theta_test=1.5*deltheta_em;
5199 vmmhz1m_test=vmmhz1m_max;
5200 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5202 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.) {
5204 theta_threshold=theta_test;
5207 theta_test=2*deltheta_em;
5208 vmmhz1m_test=vmmhz1m_max;
5209 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5211 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.)
5212 theta_threshold=theta_test;
5214 theta_test=3*deltheta_em;
5215 vmmhz1m_test=vmmhz1m_max;
5216 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5219 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.)
5220 theta_threshold=theta_test;
5222 theta_test=deltheta_had;
5223 vmmhz1m_test=vmmhz1m_max;
5224 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5226 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.)
5228 theta_threshold=theta_test;
5230 theta_threshold=sqrt(-1*deltheta_had*deltheta_had*log(anita1->VNOISE[0]/10.*anita1->maxthreshold/(hadfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle))/0.5);
5240 count_inthisloop3++;
5241 averaging_thetas3+=theta_threshold;
5246 theta_threshold*=settings1->THETA_TH_FACTOR;
5247 if (theta_threshold>0) {
5248 if (sig1->changle-theta_threshold<0 && sig1->changle+theta_threshold> PI) {
5252 else if (sig1->changle-theta_threshold>0 && sig1->changle+theta_threshold> PI) {
5253 costhetanu2=cos(sig1->changle-theta_threshold);
5256 else if (sig1->changle-theta_threshold<0 && sig1->changle+theta_threshold< PI) {
5258 costhetanu1=cos(sig1->changle+theta_threshold);
5260 else if (sig1->changle-theta_threshold>0 && sig1->changle+theta_threshold< PI) {
5261 costhetanu2=cos(sig1->changle-theta_threshold);
5262 costhetanu1=cos(sig1->changle+theta_threshold);
5269 if (theta_threshold>0) {
5271 costhetanu=costhetanu1+getRNG(RNG_DIRECTION)->Rndm()*(costhetanu2-costhetanu1);
5273 double phinu=TWOPI*getRNG(RNG_DIRECTION)->Rndm();
5274 double sinthetanu=sqrt(1-costhetanu*costhetanu);
5276 nnu =
Vector(sinthetanu*cos(phinu), sinthetanu*sin(phinu), costhetanu);
5277 nnu = nnu.ChangeCoord(refr);
5285 double angle=(PI/2.-sig1->changle)-ray1->rfexit[0].Angle(ray1->nrf_iceside[4]);
5286 double thetaposnu=posnu.Theta();
5288 costhetanu=cos(PI/2+(thetaposnu-angle));
5289 sinthetanu=sqrt(1-costhetanu*costhetanu);
5293 nnu =
Vector(-1.*sinthetanu*cos(phinu), -1.*sinthetanu*sin(phinu), -1.*costhetanu);
5297 else if (theta_threshold==-1.) {
5298 cout <<
"theta_threshold is " << theta_threshold <<
"\n";
5301 else if (emfrac<=1.E-10 && deltheta_had <= 1.E-10) {
5302 cout <<
"Error: emfrac, hadfrac are (1st place)" << emfrac <<
" " << hadfrac <<
" " <<
"\n";
5341 double dtemp1 = r_bn.
Distance(rfexit);
5343 double dtemp2 = rfexit.
Distance(posnu1);
5347 double dtemp = dtemp1 + dtemp2;
5349 vmmhz1m_max= vmmhz1m_max/dtemp;
5351 scalefactor_distance=1/dtemp;
5362 void SetupViewangles(
Signal *sig1) {
5363 double viewangle_max=90.*RADDEG;
5364 double viewangle_min=30.*RADDEG;
5365 for (
int i=0;i<NVIEWANGLE-2;i++) {
5366 viewangles[i]=viewangle_max-(viewangle_max-viewangle_min)/(
double)(NVIEWANGLE-2)*(
double)i;
5368 viewangles[NVIEWANGLE-2]=acos(1/sig1->N_DEPTH);
5369 viewangles[NVIEWANGLE-1]=90.*RADDEG;
5374 double GetThisAirColumn(
Settings* settings1,
Position r_in,
Vector nnu,
Position posnu,
double *col1,
double& cosalpha,
double& mytheta,
double& cosbeta0,
double& mybeta) {
5377 cosalpha=(r_in * nnu) / r_in.Mag();
5378 mytheta=(double)(acos(cosalpha)*DEGRAD)-90.;
5381 if (settings1->ATMOSPHERE) {
5382 int index11=int(mytheta*10.);
5383 int index12=index11+1;
5385 myair=(col1[index11]+(col1[index12]-col1[index11])*(mytheta*10.-
double(index11)))*10.;
5392 cosbeta0= (posnu * nnu) / posnu.Mag();
5393 mybeta=(double)(acos(cosbeta0)*DEGRAD)-90.;
5399 void GetAir(
double *col1) {
5401 ifstream air1(ICEMC_SRC_DIR+
"/data/atmosphere.dat");
5404 for(
int iii=0;iii<900;iii++) {
5405 air1>>nothing>>col1[iii];
5411 int TIR(
const Vector &n_surf,
const Vector &nrf2_iceside,
double N_IN,
double N_OUT) {
5412 double test=sin(nrf2_iceside.Angle(n_surf))*N_IN/N_OUT;
5423 double GetViewAngle(
const Vector &nrf2_iceside,
const Vector &nnu) {
5426 double dtemp = nrf2_iceside*nnu;
5427 if (dtemp>=1 && dtemp<1.02)
5429 if (dtemp<=-1 && dtemp>-1.02)
5439 TStyle* RootStyle() {
5440 TStyle *RootStyle =
new TStyle(
"Root-Style",
"The Perfect Style for Plots ;-)");
5450 RootStyle->SetPaperSize(TStyle::kUSLetter);
5453 RootStyle->SetCanvasColor (0);
5454 RootStyle->SetCanvasBorderSize(10);
5455 RootStyle->SetCanvasBorderMode(0);
5456 RootStyle->SetCanvasDefH (600);
5457 RootStyle->SetCanvasDefW (600);
5458 RootStyle->SetCanvasDefX (10);
5459 RootStyle->SetCanvasDefY (10);
5462 RootStyle->SetPadColor (0);
5463 RootStyle->SetPadBorderSize (10);
5464 RootStyle->SetPadBorderMode (0);
5466 RootStyle->SetPadBottomMargin(0.16);
5467 RootStyle->SetPadTopMargin (0.08);
5468 RootStyle->SetPadLeftMargin (0.18);
5469 RootStyle->SetPadRightMargin (0.05);
5470 RootStyle->SetPadGridX (0);
5471 RootStyle->SetPadGridY (0);
5472 RootStyle->SetPadTickX (1);
5473 RootStyle->SetPadTickY (1);
5476 RootStyle->SetFrameFillStyle ( 0);
5477 RootStyle->SetFrameFillColor ( 0);
5478 RootStyle->SetFrameLineColor ( 1);
5479 RootStyle->SetFrameLineStyle ( 0);
5480 RootStyle->SetFrameLineWidth ( 2);
5481 RootStyle->SetFrameBorderSize(10);
5482 RootStyle->SetFrameBorderMode( 0);
5486 RootStyle->SetHistFillColor(0);
5487 RootStyle->SetHistFillStyle(1);
5488 RootStyle->SetHistLineColor(1);
5489 RootStyle->SetHistLineStyle(0);
5490 RootStyle->SetHistLineWidth(2);
5493 RootStyle->SetFuncColor(1);
5494 RootStyle->SetFuncStyle(0);
5495 RootStyle->SetFuncWidth(1);
5498 RootStyle->SetStatBorderSize(2);
5499 RootStyle->SetStatFont (42);
5501 RootStyle->SetOptStat (0);
5502 RootStyle->SetStatColor (0);
5503 RootStyle->SetStatX (0.93);
5504 RootStyle->SetStatY (0.90);
5505 RootStyle->SetStatFontSize (0.07);
5510 RootStyle->SetTickLength ( 0.015,
"X");
5511 RootStyle->SetTitleSize ( 0.10,
"X");
5512 RootStyle->SetTitleOffset( 1.20,
"X");
5513 RootStyle->SetTitleBorderSize(0);
5515 RootStyle->SetLabelOffset( 0.015,
"X");
5516 RootStyle->SetLabelSize ( 0.050,
"X");
5517 RootStyle->SetLabelFont ( 42 ,
"X");
5518 RootStyle->SetTickLength ( 0.015,
"Y");
5519 RootStyle->SetTitleSize ( 0.10,
"Y");
5520 RootStyle->SetTitleOffset( 0.600,
"Y");
5521 RootStyle->SetLabelOffset( 0.015,
"Y");
5522 RootStyle->SetLabelSize ( 0.050,
"Y");
5523 RootStyle->SetLabelFont ( 42 ,
"Y");
5524 RootStyle->SetTitleFont (42,
"XY");
5525 RootStyle->SetTitleColor (1);
5528 RootStyle->SetOptFit (0);
5529 RootStyle->SetMarkerStyle(20);
5530 RootStyle->SetMarkerSize(0.4);
5539 void GetFresnel(
Roughness *rough1,
int ROUGHNESS_SETTING,
const Vector &surface_normal,
5544 double emfrac,
double hadfrac,
double deltheta_em_max,
double deltheta_had_max,
5545 double &t_coeff_pokey,
double &t_coeff_slappy,
5550 double incident_angle = surface_normal.Angle(firn_rf);
5551 double transmitted_angle = surface_normal.Angle(air_rf);
5556 Vector perp = air_rf.Cross(surface_normal).Unit();
5558 Vector air_parallel = perp.Cross(air_rf).Unit();
5560 Vector firn_parallel = perp.Cross(firn_rf).Unit();
5563 double pol_perp_firn = pol*perp;
5564 double pol_parallel_firn = pol*firn_parallel;
5565 double pol_perp_air=0, pol_parallel_air=0;
5567 double r_coeff_pokey = tan(incident_angle - transmitted_angle) / tan(incident_angle + transmitted_angle);
5568 t_coeff_pokey = sqrt((1. - r_coeff_pokey*r_coeff_pokey));
5569 pol_parallel_air = t_coeff_pokey * pol_parallel_firn;
5571 double r_coeff_slappy = sin(incident_angle - transmitted_angle) / sin(incident_angle + transmitted_angle);
5572 t_coeff_slappy = sqrt((1. - r_coeff_slappy*r_coeff_slappy));
5573 pol_perp_air = t_coeff_slappy * pol_perp_firn;
5575 mag=sqrt( tan(incident_angle) / tan(transmitted_angle) );
5577 fresnel = sqrt( pow(efield * pol_perp_air, 2) + pow(efield * pol_parallel_air, 2)) / efield;
5579 pol = (pol_perp_air * perp + pol_parallel_air * air_parallel).Unit();
5603 const Vector nuvector = interaction1->
nnu;
5606 Vector zcoordvector = ray1->nsurf_rfexit;
5607 zcoordvector=zcoordvector.Unit();
5612 Vector xcoordvector = nuvector-(zcoordvector.Dot(nuvector))*zcoordvector;
5613 xcoordvector = xcoordvector.Unit();
5615 const Vector ycoordvector = zcoordvector.Cross(xcoordvector);
5619 if (interaction1->
nnu.Dot(zcoordvector)>0)
5620 origin_brian_tmp=interaction1->
nuexit;
5623 nnu_flipped=nnu_flipped-2.*nnu_flipped.Dot(zcoordvector)*zcoordvector;
5626 if (Ray::WhereDoesItLeave(interaction1->posnu,nnu_flipped,antarctica,nuexit_flipped))
5627 origin_brian_tmp=nuexit_flipped;
5631 r_bn_tmp=r_bn_tmp.ChangeCoord(xcoordvector,ycoordvector);
5634 double balloonphi = r_bn_tmp.Phi();
5636 balloonphi=balloonphi-2*PI;
5638 double balloontheta = r_bn_tmp.Theta();
5641 balloontheta = PI-balloontheta;
5652 void interrupt_signal_handler(
int sig){
5653 signal(sig, SIG_IGN);
5660 #ifdef ANITA_UTIL_EXISTS 5662 int GetIceMCAntfromUsefulEventAnt(
Settings *settings1,
int UsefulEventAnt){
5666 int IceMCAnt = UsefulEventAnt;
5667 if ((settings1->WHICH==9 || settings1->WHICH==10) && UsefulEventAnt<16) {
5668 IceMCAnt = (UsefulEventAnt%2==0)*UsefulEventAnt/2 + (UsefulEventAnt%2==1)*(UsefulEventAnt/2+8);
UInt_t eventNumber
Event number from software.
double weight_bestcase
what weight1 would be if whole earth had density of crust - for quick and dirty calculation of best c...
int getDirectionAndEnergy(Vector *nudir, double t, double &nuE, double minE=1e9, double maxE=1e12)
Double_t h_component_k[48]
Component of the e-field along the rx h-plane.
Double_t rfExitPos[5][3]
Position where the RF exits the ice- 5 iterations, 3 dimensions each.
Double_t fNoiseAtDigitizer[12 *9][260]
Array of noise at digitizer.
int WHICHPATH
0=fixed balloon position,1=randomized,2=ANITA-lite GPS data,3=banana plot
void PassesTrigger(Settings *settings1, Anita *anita1, int discones_passing, int mode, int *l3trig, int l2trig[Anita::NPOL][Anita::NTRIGGERLAYERS_MAX], int l1trig[Anita::NPOL][Anita::NTRIGGERLAYERS_MAX], int antennaclump, int loctrig[Anita::NPOL][Anita::NLAYERS_MAX][Anita::NPHI_MAX], int loctrig_nadironly[Anita::NPOL][Anita::NPHI_MAX], int inu, int *thispasses, bool noiseOnly=false)
Evaluate the full trigger simulation for a given payload configuration.
Double_t weight
Weight assigned by icemc.
Int_t nu_pdg
Neutrino PDG code.
Position r_bn
position of balloon
Double_t sourceLat
RF position when leaving the ice: Latitude (using icemc model)
int currentint
Ditto - Stephen.
double fVolts[12 *9][260]
Array of unwrapped (unless kNoCalib) voltages for each channel.
Double_t showerE
Shower energy.
double r_fromballoon[2]
distance from interaction to balloon for each ray
Double_t timeWeight
Relative Time weight assigned by icemc.
double chord_kgm2_ice
from ice entrance to interaction point
void InitializeEachBand(Anita *anita1)
Initialize trigger bands.
Adu5Pat – The ADU5 Position and Attitude Data.
double Distance(const Position &second) const
Returns chord distance (direct distance between two vectors)
void TriggerPath(Settings *settings1, Anita *anita1, int ant, Balloon *bn1)
Apply trigger path.
static const int NPOL
number of polarizations
Double_t n_component_k[48]
Component of the e-field along the normal.
int GetSigma(double pnu, double &sigma, double &len_int_kgm2, Settings *settings1, int nu_nubar, int currentint)
Neutrino-nucleon cross-sections using model chosen.
Double_t fDiodeOutput[12 *9][260]
Array of tunnel diode output.
Vector nnu_banana
Forced neutrino direction.
Position nu_banana_surface
The location of the surface above the forced neutrino interaction point.
Float_t latitude
In degrees.
void ApplyAntennaGain(Settings *settings1, Anita *anita1, Balloon *bn1, Screen *panel1, int ant, Vector &n_eplane, Vector &n_hplane, Vector &n_normal)
Apply the antenna gain.
double chord
chord in m from earth entrance to rock-ice boundary
Int_t source_index
Name of the source.
double weight_tau_prob
Weight for tau neutrino to interact, create a tau, tau survives and decays in the ice...
double r_exit2bn
exit to balloon
Radiation from interaction.
double dnutries
product of dtryingdirection and dtryingposition
double fTimes[12 *9][260]
Array of unwrapped (unless kNoCalib) times for each channel.
Double_t sourceLon
RF position when leaving the ice: Longitude (using icemc model)
Double_t sourceAlt
RF position when leaving the ice: Altitude (using icemc model)
void TimeShiftAndSignalFluct(Settings *settings1, Anita *anita1, int ilayer, int ifold, double volts_rx_rfcm_lab_e_all[48][512], double volts_rx_rfcm_lab_h_all[48][512])
Time shift and fluctuate signal.
Float_t longitude
In degrees.
Double_t thresholds[12 *9]
Channel thresholds used in icemc.
Double_t sourceTimeWeight
Relative Time weight for the given source assigned by icemc.
UInt_t realTime
Time from the GPS unit.
const char * ICEMC_SRC_DIR()
Double_t n_component[48]
Normal comp along polarization.
Double_t balloonPos[3]
Balloon position.
double dtryingdirection
weighting factor: how many equivalent tries each neutrino counts for after having reduced angular pha...
void saveTriggerWaveforms(Anita *anita1, double sig0[48], double sig1[48], double noise0[48], double noise1[48])
Save signal and noise waveforms at trigger.
double costheta_nutraject
theta of nnu with earth center to balloon as z axis
double weight_nu
Weight for neutrino that survives to posnu.
static constexpr double banana_sigma
NSIGMA in the case of a banana plot.
double noise_eachband[2][Anita::NBANDS_MAX]
Noise in each band.
Double_t vmmhz[128]
V/m/MHz at balloon (128 frequency bins)
Double_t l_int
Interaction length.
Double_t SNRAtDigitizer[12 *9]
Array of SNR at digitizer.
double Lat() const
Returns latitude, where the +z direction is at 0 latitude.
Class that handles the channel trigger.
Reads in and stores input settings for the run.
Double_t fSignalAtTrigger[12 *9][260]
Array of signal at trigger.
Functions you need to generate a primary interaction including cross sections and picking charged cur...
TString objName
Declination of source.
Double_t maxSNRAtTriggerV
Max SNR at trigger V-POL.
int CENTER
whether or not to center one phi sector of the payload on the incoming signal (for making signal effi...
This class is a 3-vector that represents a position on the Earth's surface.
void saveDigitizerWaveforms(Anita *anita1, double sig0[48], double sig1[48], double noise0[48], double noise1[48])
Save signal and noise waveforms at digitizer.
Short_t tuffIndex
TUFF configuration index.
Double_t dec
Right ascension of source.
double logchord
log_10 of chord length earth entrance to where it enters ice
Double_t h_component[48]
H comp along polarization.
UInt_t realTime
unixTime of readout
Stores everything about a particular neutrino interaction. Interaction.
Double_t rfExitNor[5][3]
Normal vector in direction of exit point to balloon - 5 iterations.
UsefulAnitaEvent – The Calibrated Useful Anita Event object.
double pickY(Settings *settings1, double pnu, int nu_nubar, int currentint)
pick inelasticity y according to chosen model
double GetTauWeight(Primaries *primary1, Settings *settings1, IceModel *antarctica1, Interaction *interaction1, double pnu, int nu_nubar, double &ptauf, int &crust_entered)
GetTauWeight is the function that will calculate the probability that a tau neutrino will interact al...
double SurfaceDistance(const Position &second, double local_surface) const
Returns "surface distance" between two positions.
Float_t altitude
In metres.
Double_t hitangle_h[48]
Hit angles rel. to h plane stored for each antenna.
double total_kgm2
number of earth layers traversed
static SourceModel * getSourceModel(const char *key, Restriction r=Restriction())
double altitude_int_mirror
depth of the mirror point of interaction.
double r_exit2bn_measured
exit to balloon deduced from measured theta
Double_t fSignalAtDigitizer[12 *9][260]
Array of signal at digitizer.
void DigitizerPath(Settings *settings1, Anita *anita1, int ant, Balloon *bn1)
Apply digitizer path.
int nuflavorint
Added by Stephen for output purposes.
Position r_in
position where neutrino enters the earth
double phi_nutraject
phi of nnu with earth center to balloon as z axis
TruthAnitaEvent – The Truth ANITA Event.
Double_t e_component[48]
E comp along polarization.
double threshold_eachband[2][Anita::NBANDS_MAX]
Threshold in each band.
Shape of the earth, ice thicknesses, profiles of earth layers, densities, neutrino absorption...
Double_t weight1
Absorption weight assigned by icemc.
Double_t weight_prob
weight including probability of interacting!
static const int NLAYERS_MAX
max number of layers (in smex design, it's 4)
int fCapacitorNum[12 *9][260]
Array of capacitor numbers.
double altitude_int
depth of interaction
void addSource(Source *source)
double weight_nu_prob
Weight for neutrino that survives to posnu and interacts in the ice.
int passes_eachband[2][Anita::NBANDS_MAX]
Whether the signal passes or not each band.
double banana_volts
Total voltage measured at a spot on the sky.
Double_t fNoiseAtTrigger[12 *9][260]
Array of noise at trigger.
Double_t nuDir[3]
Neutrino direction.
double chord_kgm2_bestcase
the chord the neutrino would traverse if it all was crust density
static constexpr double banana_y
Elasticity. 0.2 is an average number.
int iceinteraction
whether or not there is an interaction in the ice
double bwslice_volts_polh[5]
Sum voltage for each slice in bandwidth for the h polarization.
void WhichBandsPass(Settings *settings1, Anita *anita1, GlobalTrigger *globaltrig1, Balloon *bn1, int ilayer, int ifold, double dangle, double emfrac, double hadfrac, double thresholds[2][5])
Which bands passes the trigger.
Double_t e_component_k[48]
Component of e-field along the rx e-plane.
Double_t maxSNRAtDigitizerV
Max SNR at digitizer V-POL.
Double_t nuMom
Neutrino momentum.
int fNumPoints[12 *9]
Number of poins per channel.
Int_t fRFSpike
Flag raised if the ADC value is too large or small in RF.
Handles everything related to balloon positions, payload orientation over the course of a flight...
Position nuexitice
place where neutrino would have left the ice
Position nuexit
place where neutrino would have left the earth
This class represents a three-vector. Operators are overloaded to provide for the familiar operations...
Position nu_banana
The forced interaction point of the neutrino for the banana plots.
static const int NPHI_MAX
max number of antennas around in phi (in smex, 16)
double bwslice_volts_pole[5]
Sum voltage for each slice in bandwidth for the e polarization.
double ptauf
Final energy of the tau.
Float_t heading
0 is facing north, 180 is facing south
Double_t maxSNRAtTriggerH
Max SNR at trigger H-POL.
Handles event counting as cuts are made.
double d2
ice-rock boundary to interaction point in m
UInt_t eventNumber
Software event number.
double signal_eachband[2][Anita::NBANDS_MAX]
Signal in each band.
double r_fromballoon_db
same, for double bangs
Double_t nuPos[3]
Neutrino position.
Double_t fTimes[12 *9][260]
Array of unwrapped (unless kNoCalib) times for each channel.
Double_t hitangle_e[48]
Hit angles rel. to e plane stored for each antenna.
Position r_enterice
position where neutrino enters the ice
Double_t SNRAtTrigger[12 *9]
Array of SNR at trigger.
double Lon() const
Returns longitude.
double d1
same as chord in m (earth entrance to rock-ice boundary)
Double_t maxSNRAtDigitizerH
Max SNR at digitizer H-POL.
static constexpr double banana_signal_fluct
Turn off noise for banana plots (settings1->SIGNAL_FLUCT) (shouldn't matter)
void setUpWeights(double t0, double t1, double minE=1e9, double maxE=1e12, int N=1e6)
Double_t phaseWeight
Phase weight assigned by icemc.
Double_t balloonDir[3]
Balloon direction.
Ice thicknesses and water depth.
Vector nnu
direction of neutrino (+z in south pole direction)