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;
134 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};
135 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.};
138 double MIN_LOGWEIGHT=-3;
139 double MAX_LOGWEIGHT=-1;
141 double eventsfound_binned[NBINS];
142 double eventsfound_binned_e[NBINS];
143 double eventsfound_binned_mu[NBINS];
144 double eventsfound_binned_tau[NBINS];
150 double error_e_plus=0;
151 double error_mu_plus=0;
152 double error_tau_plus=0;
153 double error_minus=0;
154 double error_e_minus=0;
155 double error_mu_minus=0;
156 double error_tau_minus=0;
158 double gain_dipole=2.15;
159 double changle_deg=0;
163 double RANDOMISEPOL=0.;
167 std::vector<int> NNU_per_source;
168 std::vector<double> eventsfound_per_source;
169 std::vector<double> eventsfound_prob_per_source;
171 double volume_thishorizon;
172 double realtime_this;
173 double longitude_this;
174 double latitude_this;
175 double altitude_this;
176 double heading_this=0.;
181 double pnu=pow(10., 20);
186 double SIGNALRADIUS=2.;
189 const double FREQ_LOW_DISCONES=120.E6;
190 const double FREQ_HIGH_DISCONES=1000.E6;
192 double bwslice_vnoise_thislayer[4];
193 int passes_thisevent=0;
194 int unmasked_thisevent=0;
196 int discones_passing;
198 double heff_discone=0;
199 double polarfactor_discone=0.;
201 double volts_discone=0.;
202 double vnoise_discone=0.;
205 double BW_DISCONES=300.E6-120.E6;
214 double t_coeff_pokey, t_coeff_slappy;
217 double e_comp_max1=0;
218 double h_comp_max1=0;
219 double e_comp_max2=0;
220 double h_comp_max2=0;
221 double e_comp_max3=0;
222 double h_comp_max3=0;
225 double diff_3tries=0;
228 double costheta_inc=0;
229 double costheta_exit=0;
230 double theta_rf_atbn;
231 double theta_rf_atbn_measured;
232 double theta_nnu_atbn;
233 double theta_nnu_rf_diff_atbn;
237 double phi_nnu_rf_diff_atbn;
245 double costhetanu=-1000;
250 double nearthlayers=0;
251 double weight_prob=0.;
256 double pieceofkm2sr=0;
259 int count_inthisloop1=0;
260 int count_inthisloop2=0;
261 int count_inthisloop3=0;
262 double averaging_thetas1=0;
263 double averaging_thetas2=0;
264 double averaging_thetas3=0;
269 int count_asktrigger=0;
270 int count_asktrigger_nfb=0;
273 double count_passestrigger_w=0;
276 int allcuts[2]={0, 0};
278 double allcuts_weighted[2]={0, 0};
279 double allcuts_weighted_polarization[3]={0, 0, 0};
283 int count_chanceofsurviving=0;
285 int count_chanceinhell0=0;
289 double count_chanceinhell2_w=0;
293 int count_chordgoodlength=0;
294 int count_d2goodlength=0;
298 double sum_frac_db[3];
300 const int NBINS_DISTANCE=28;
301 double eventsfound_binned_distance[NBINS_DISTANCE] = {0.};
302 int index_distance=0;
303 double km3sr_distance[NBINS_DISTANCE] = {0.};
304 double error_distance_plus[NBINS_DISTANCE] = {0.};
305 double error_distance_minus[NBINS_DISTANCE] = {0.};
306 int eventsfound_binned_distance_forerror[NBINS_DISTANCE][NBINS] = {{0}};
312 int count_passestrigger_nfb=0;
313 double percent_increase_db=0;
314 double percent_increase_nfb=0;
315 double percent_increase_total=0;
317 double error_km3sr_nfb=0;
318 double error_percent_increase_nfb=0;
323 int count_dbexitsice=0;
326 double eventsfound_nfb_binned[NBINS];
329 double heff_max=0.62639;
332 double scalefactor_distance=0;
333 double scalefactor_attenuation=0;
334 double MAX_ATTENLENGTH=1671;
337 double dviewangle_deg=0;
339 double forseckel[NVIEWANGLE][Anita::NFREQ];
340 double viewangles[NVIEWANGLE];
341 double GetAirDistance(
double altitude_bn,
double beta);
349 int max_antenna0 = -1;
350 int max_antenna1 = -1;
351 int max_antenna2 = -1;
352 double max_antenna_volts0 = 0;
353 double max_antenna_volts0_em = 0;
355 double max_antenna_volts1 = 0;
356 double max_antenna_volts2 = 0;
358 double rx0_signal_eachband[2][5];
359 double rx0_threshold_eachband[2][5];
360 double rx0_noise_eachband[2][5];
361 int rx0_passes_eachband[2][5];
368 double vmmhz1m_visible = 0;
369 int freq_bins = Anita::NFREQ;
370 double total_kgm2 = 0;
372 double r_in_array[3];
373 double nsurf_rfexit_array[3];
374 double nsurf_rfexit_db_array[3];
375 double r_bn_array[3];
376 double n_bn_array[3];
377 double posnu_array[3];
378 double nrf_iceside_array[5][3];
379 double nrf_iceside_db_array[5][3];
380 double ant_max_normal0_array[3];
381 double ant_max_normal1_array[3];
382 double ant_max_normal2_array[3];
383 double n_pol_array[3];
384 double n_exit2bn_array[5][3];
385 double r_enterice_array[3];
386 double n_exit2bn_db_array[5][3];
387 double rfexit_array[5][3];
388 double rfexit_db_array[5][3];
390 int times_crust_entered_det=0;
391 int times_mantle_entered_det=0;
392 int times_core_entered_det=0;
394 int mantle_entered=0;
396 int neutrinos_passing_all_cuts=0;
397 double sum_weights=0;
400 int xsecParam_nutype = 0;
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 while ((clswitch = getopt(argc, argv,
"t:i:o:r:n:e:x:")) != EOF) {
533 nnu_tmp=atoi(optarg);
534 cout <<
"Changed number of simulated neutrinos to " << nnu_tmp << endl;
537 trig_thresh=atof(optarg);
541 cout <<
"Changed input file to: " << input << endl;
545 cout <<
"Changed output directory to: " << string(outputdir.Data()) << endl;
546 stemp=
"mkdir -p " + string(outputdir.Data());
547 system(stemp.c_str());
550 exp_tmp=atof(optarg);
551 cout <<
"Changed neutrino energy exponent to " << exp_tmp << endl;
554 startNu=atoi(optarg);
555 cout <<
"Running icemc for just 1 event with eventNumber : " << startNu << endl;
559 stringstream convert(run_num);
566 settings1->SEED=settings1->SEED +run_no;
567 cout <<
"seed is " << settings1->SEED << endl;
570 stemp=string(outputdir.Data())+
"/nu_position"+run_num+
".txt";
571 ofstream nu_out(stemp.c_str(), ios::app);
573 stemp=string(outputdir.Data())+
"/veff"+run_num+
".txt";
574 ofstream veff_out(stemp.c_str(), ios::app);
576 stemp=string(outputdir.Data())+
"/distance"+run_num+
".txt";
577 ofstream distanceout(stemp.c_str());
579 stemp=string(outputdir.Data())+
"/debug"+run_num+
".txt";
580 fstream outfile(stemp.c_str(), ios::out);
582 stemp=string(outputdir.Data())+
"/forbrian"+run_num+
".txt";
583 fstream forbrian(stemp.c_str(), ios::out);
585 stemp=string(outputdir.Data())+
"/al_voltages_direct"+run_num+
".dat";
586 fstream al_voltages_direct(stemp.c_str(), ios::out);
588 stemp=string(outputdir.Data())+
"/events"+run_num+
".txt";
589 ofstream eventsthatpassfile(stemp.c_str());
591 stemp=string(outputdir.Data())+
"/numbers"+run_num+
".txt";
592 ofstream fnumbers(stemp.c_str());
594 stemp=string(outputdir.Data())+
"/output"+run_num+
".txt";
595 ofstream foutput(stemp.c_str(), ios::app);
597 stemp=string(outputdir.Data())+
"/slacviewangles"+run_num+
".dat";
598 ofstream fslac_viewangles(stemp.c_str());
600 stemp=string(outputdir.Data())+
"/slac_hitangles"+run_num+
".dat";
601 ofstream fslac_hitangles(stemp.c_str());
613 settings1->ReadInputs(input.c_str(), foutput, NNU, RANDOMISEPOL);
614 settings1->ApplyInputs(anita1, sec1, sig1, bn1, ray1);
620 settings1->SEED=settings1->SEED + run_no;
621 setSeed(settings1->SEED);
623 bn1->InitializeBalloon();
624 anita1->
Initialize(settings1, foutput, inu, outputdir);
632 anita1->powerthreshold[4]=trig_thresh;
636 settings1->EXPONENT=exp_tmp;
639 if (!strcasecmp(settings1->SOURCE.c_str(),
"CUSTOM"))
643 src_model->
addSource(
new Source(
"Custom Source", (settings1->CUSTOM_RA)/15, settings1->CUSTOM_DEC,
647 else if(strcasecmp(settings1->SOURCE.c_str(),
"NONE"))
651 SourceModel::Restriction::fromString(settings1->WHICH_START_TIME.c_str()),
652 SourceModel::Restriction::fromString(settings1->WHICH_END_TIME.c_str())));
655 if (settings1->SOURCE_USE_EXPONENT && settings1->EXPONENT < 22)
657 settings1->SOURCE_MIN_E=settings1->EXPONENT;
658 settings1->SOURCE_MAX_E=settings1->EXPONENT;
662 double src_min = TMath::Power(10,settings1->SOURCE_MIN_E-9);
663 double src_max = TMath::Power(10,settings1->SOURCE_MAX_E-9);
670 printf(
"Using Source Model %s\n", src_model->getName());
671 printf(
"Setting up weights: %g %g %g %g\n", bn1->min_time, bn1->max_time, src_min, src_max);
672 src_model->
setUpWeights(bn1->min_time, bn1->max_time , src_min, src_max);
678 rough1->SetRoughScale(settings1->ROUGHSIZE);
682 if(!src_model && spectra1->IsSpectrum())
684 cout <<
"Using an energy spectrum" << endl;
689 cout<<
"Source Model Bounds are 10^" << settings1->SOURCE_MIN_E <<
"eV to 10^" << settings1->SOURCE_MAX_E <<
"eV" << endl;
692 else if (!src_model && spectra1->IsMonoenergetic())
694 cout <<
"Using a monoenergetic exponent of: " << settings1->EXPONENT << endl;
697 std::cout <<
"----------------------" << std::endl;
704 time_t raw_start_time = time(NULL);
705 struct tm * start_time = localtime(&raw_start_time);
707 cout <<
"Date and time at start of run are: " << asctime (start_time) <<
"\n";
709 if (settings1->FORSECKEL)
710 SetupViewangles(sig1);
719 Tools::Zero(sum_frac, 3);
720 Tools::Zero(sum_frac_db, 3);
731 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;
736 double viewangle_temp=0;
737 double viewangle_deg=0;
738 double cosviewangle=0;
740 double nsigma_offaxis=0;
741 double theta_threshold=0;
742 double theta_threshold_deg=0;
744 double nsigma_em_threshold=0;
745 double nsigma_had_threshold=0;
746 double slopeyangle=0;
748 int beyondhorizon = 0;
751 double vmmhz1m_max=0;
752 double vmmhz_lowfreq=0.;
753 double bestcase_atten=0;
754 double vmmhz1m_fresneledonce=0;
755 double vmmhz1m_fresneledtwice=0;
756 double vmmhz[Anita::NFREQ];
759 double vmmhz_em[Anita::NFREQ];
760 double vmmhz_temp=1.;
761 double vmmhz_min_thatpasses=1000;
764 double deltheta_em[Anita::NFREQ], deltheta_had[Anita::NFREQ];
765 double deltheta_em_max, deltheta_had_max;
766 double deltheta_em_mid2, deltheta_had_mid2;
769 double emfrac, hadfrac, sumfrac;
770 int n_interactions=1;
771 double emfrac_db, hadfrac_db;
777 double d12=interaction1->
d1;
778 double d22=interaction1->
d2;
780 double logchord2=interaction1->
logchord;
785 double r_exit2bn2=interaction1->
r_exit2bn;
793 double volts_rx_max=0;
794 double volts_rx_ave=0;
795 double volts_rx_sum=0;
797 double volts_rx_max_highband;
798 double volts_rx_max_lowband;
799 double volts_rx_rfcm_lab_e_all[48][512];
800 double volts_rx_rfcm_lab_h_all[48][512];
803 static double e_component[Anita::NANTENNAS_MAX]={0};
804 static double h_component[Anita::NANTENNAS_MAX]={0};
805 static double n_component[Anita::NANTENNAS_MAX]={0};
807 static double e_component_kvector[Anita::NANTENNAS_MAX]={0};
808 static double h_component_kvector[Anita::NANTENNAS_MAX]={0};
809 static double n_component_kvector[Anita::NANTENNAS_MAX]={0};
811 Vector n_eplane = const_z;
812 Vector n_hplane = -const_y;
813 Vector n_normal = const_x;
818 double hitangle_e, hitangle_h;
819 double hitangle_e_all[Anita::NANTENNAS_MAX];
820 double hitangle_h_all[Anita::NANTENNAS_MAX];
823 double len_int_kgm2=0;
824 double eventsfound_db=0;
825 double eventsfound_nfb=0;
831 double dist_int_bn_2d=0;
832 double dist_int_bn_2d_chord=0;
840 double nuexitlength=0;
842 double nuentrancelength=0;
844 double icethickness=0;
845 double theta_pol_measured;
850 double nutauweightsum=0;
851 double tauweightsum=0;
852 double nutauweight=0;
854 int tauweighttrigger=0;
855 double sample_x=0, sample_y = 0;
867 Vector n_nutraject_ontheground;
875 int l2trig[
Anita::NPOL][Anita::NTRIGGERLAYERS_MAX];
877 int l1trig[
Anita::NPOL][Anita::NTRIGGERLAYERS_MAX];
884 int nchannels_triggered = 0;
885 int nchannels_perrx_triggered[48];
891 Tools::Zero(count1->npass, 2);
896 Tools::Zero(eventsfound_binned, NBINS);
897 Tools::Zero(eventsfound_binned_e, NBINS);
898 Tools::Zero(eventsfound_binned_mu, NBINS);
899 Tools::Zero(eventsfound_binned_tau, NBINS);
900 Tools::Zero(eventsfound_nfb_binned, NBINS);
905 TH2F *ref_int_coord=
new TH2F(
"ref_int_coord",
"", 600, -3000, 3000, 500, -2500, 2500);
906 ref_int_coord->SetMarkerSize(0.7);;
907 ref_int_coord->SetMarkerStyle(30);
908 ref_int_coord->SetMarkerColor(kBlue);
910 TH2F *dir_int_coord=
new TH2F(
"dir_int_coord",
"", 600, -3000, 3000, 500, -2500, 2500);
911 dir_int_coord->SetMarkerSize(0.7);
912 dir_int_coord->SetMarkerStyle(30);
914 TH1F *prob_eachphi_bn=
new TH1F(
"prob_eachphi_bn",
"prob_eachphi_bn", 100, 0., 6.3);
915 TH1F *prob_eachilon_bn=
new TH1F(
"prob_eachilon_bn",
"prob_eachilon_bn", 180, 0., 180.);
916 TH2F *h6=
new TH2F(
"theta_vs_hitangle_h",
"theta_vs_hitangle_h", 100, -3.14, 3.14, 100, -1.1, 1.1);
917 TH1F *h10=
new TH1F(
"hitangle_e",
"hitangle_e", 20, -1.6, 1.6);
918 TH1F *hy=
new TH1F(
"hy",
"hy", 100, 0., 1.);
919 TH1F *fraction_sec_muons =
new TH1F(
"fraction_sec_muons",
"fraction_sec_muons", 100, 0., 1.);
920 TH1F *fraction_sec_taus=
new TH1F(
"fraction_sec_taus",
"fraction_sec_taus", 100, 0., 1.);
921 TH1F *n_sec_muons=
new TH1F(
"n_sec_muons",
"n_sec_muons", 100, 0., 10.);
922 TH1F *n_sec_taus=
new TH1F(
"n_sec_taus",
"n_sec_taus", 100, 0., 10.);
923 TH1F *sampleweights=
new TH1F(
"sampleweights",
"sampleweights", 100, -5., 0.);
927 stemp=string(outputdir.Data())+
"/icefinal"+run_num+
".root";
928 TFile *hfile =
new TFile(stemp.c_str(),
"RECREATE",
"ice");
929 TTree *tree2 =
new TTree(
"h2000",
"h2000");
931 tree2->Branch(
"inu", &inu,
"inu/I");
932 tree2->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
933 tree2->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
934 tree2->Branch(
"scalefactor_distance", &scalefactor_distance,
"scalefactor_distance/D");
935 tree2->Branch(
"scalefactor_attenuation", &scalefactor_attenuation,
"scalefactor_attenuation/D");
937 TTree *tree3 =
new TTree(
"h3000",
"h3000");
938 tree3->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
939 tree3->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
940 tree3->Branch(
"theta_threshold_deg", &theta_threshold_deg,
"theta_threshold_deg/D");
941 tree3->Branch(
"nsigma_em_threshold", &nsigma_em_threshold,
"nsigma_em_threshold/D");
942 tree3->Branch(
"nsigma_had_threshold", &nsigma_had_threshold,
"nsigma_had_threshold/D");
943 tree3->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
944 tree3->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
945 tree3->Branch(
"vmmhz_max", &vmmhz_max,
"vmmhz_max/D");
946 tree3->Branch(
"vmmhz_min", &vmmhz_min,
"vmmhz_min/D");
947 tree3->Branch(
"dviewangle_deg", &dviewangle_deg,
"dviewangle_deg/D");
948 tree3->Branch(
"viewangle_deg", &viewangle_deg,
"viewangle_deg/D");
949 tree3->Branch(
"changle_deg", &changle_deg,
"changle_deg/D");
950 tree3->Branch(
"cosviewangle", &cosviewangle,
"cosviewangle/D");
951 tree3->Branch(
"emfrac", &emfrac,
"emfrac/D");
952 tree3->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
954 TTree *tree5 =
new TTree(
"h5000",
"h5000");
955 tree5->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
956 tree5->Branch(
"inu", &inu,
"inu/I");
957 tree5->Branch(
"nuexitlength", &nuexitlength,
"nuexitlength/D");
958 tree5->Branch(
"nuexitice", &nuexitice,
"nuexitice");
959 tree5->Branch(
"vmmhz_max", &vmmhz_max,
"vmmhz_max");
960 tree5->Branch(
"maxtaper", &maxtaper,
"maxtaper");
961 tree5->Branch(
"inu", &inu,
"inu/I");
962 tree5->Branch(
"whichray", &whichray,
"whichray/I");
963 tree5->Branch(
"pnu", &pnu,
"pnu/D");
964 tree5->Branch(
"costhetanu", &costhetanu,
"costhetanu/D");
965 tree5->Branch(
"viewangle", &viewangle,
"viewangle/D");
966 tree5->Branch(
"offaxis", &offaxis,
"offaxis/D");
967 tree5->Branch(
"nsigma_offaxis", &nsigma_offaxis,
"nsigma_offaxis/D");
968 tree5->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
969 tree5->Branch(
"emfrac", &emfrac,
"emfrac/D");
970 tree5->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
971 tree5->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
972 tree5->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
973 tree5->Branch(
"weight1", &weight1,
"weight1/D");
974 tree5->Branch(
"nearthlayers", &nearthlayers,
"nearthlayers/D");
975 tree5->Branch(
"logchord", &logchord2,
"interaction1->logchord/D");
976 tree5->Branch(
"diff_3tries", &diff_3tries,
"diff_3tries/D");
977 tree5->Branch(
"fresnel2", &fresnel2,
"fresnel2/D");
978 tree5->Branch(
"costheta_inc", &costheta_inc,
"costheta_inc/D");
979 tree5->Branch(
"costheta_exit", &costheta_exit,
"costheta_exit/D");
980 tree5->Branch(
"deltheta_em", &deltheta_em[0],
"deltheta_em/D");
981 tree5->Branch(
"deltheta_had", &deltheta_had[0],
"deltheta_had/F");
982 tree5->Branch(
"r_fromballoon", &interaction1->
r_fromballoon[0],
"r_fromballoon/F");
983 tree5->Branch(
"theta_in", &theta_in,
"theta_in/D");
984 tree5->Branch(
"lat_in", &lat_in,
"lat_in/D");
987 TTree *tree6 =
new TTree(
"h6000",
"h6000");
988 tree6->Branch(
"volts_rx_0", &volts_rx_0,
"volts_rx_0/D");
989 tree6->Branch(
"volts_rx_1", &volts_rx_1,
"volts_rx_1/D");
990 tree6->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
991 tree6->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
992 tree6->Branch(
"theta_in", &theta_in,
"theta_in/D");
993 tree6->Branch(
"chord_kgm2_bestcase", &chord_kgm2_bestcase2,
"chord_kgm2_bestcase/D");
994 tree6->Branch(
"chord_kgm2_ice", &interaction1->
chord_kgm2_ice,
"chord_kgm2_ice/D");
995 tree6->Branch(
"costheta_nutraject", &interaction1->
costheta_nutraject,
"costheta_nutraject/D");
996 tree6->Branch(
"weight1", &weight1,
"weight1/D");
997 tree6->Branch(
"weight_bestcase", &weight_bestcase2,
"weight_bestcase/D");
998 tree6->Branch(
"whichray", &whichray,
"whichray/I");
999 tree6->Branch(
"mybeta", &mybeta,
"mybeta/D");
1000 tree6->Branch(
"longitude", &longitude_this,
"longitude/D");
1003 TTree *tree6b =
new TTree(
"h6001",
"h6001");
1004 tree6b->Branch(
"bwslice_vnoise", bwslice_vnoise_thislayer,
"bwslice_vnoise_thislayer[4]/D");
1006 TTree *tree7 =
new TTree(
"h7000",
"h7000");
1007 tree7->Branch(
"emfrac", &emfrac,
"emfrac/D");
1008 tree7->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1009 tree7->Branch(
"current", &interaction1->
currentint,
"currentint/I");
1010 tree7->Branch(
"nuflavor", &interaction1->
nuflavorint,
"nuflavorint/I");
1011 tree7->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
1012 tree7->Branch(
"slopeyangle", &slopeyangle,
"slopeyangle/D");
1014 TTree *jaimetree=
new TTree(
"jaimetree",
"jaimetree");
1015 jaimetree->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
1016 jaimetree->Branch(
"emfrac", &emfrac,
"emfrac/D");
1017 jaimetree->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1018 jaimetree->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
1019 jaimetree->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
1020 jaimetree->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
1021 jaimetree->Branch(
"vmmhz1m_visible", &vmmhz1m_visible,
"vmmhz1m_visible/D");
1023 TTree *viewangletree=
new TTree(
"viewangletree",
"viewangletree");
1024 viewangletree->Branch(
"dviewangle_deg", &dviewangle_deg,
"dviewangle_deg/D");
1025 viewangletree->Branch(
"emfrac", &emfrac,
"emfrac/D");
1026 viewangletree->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1027 viewangletree->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
1028 viewangletree->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
1029 viewangletree->Branch(
"theta_threshold_deg", &theta_threshold_deg,
"theta_threshold_deg/D");
1030 viewangletree->Branch(
"dnutries", &interaction1->
dnutries,
"dnutries/D");
1031 viewangletree->Branch(
"viewangle", &viewangle,
"viewangle/D");
1032 viewangletree->Branch(
"chord", &interaction1->
chord,
"dnutries/D");
1034 TTree *neutrino_positiontree=
new TTree(
"neutrino_positiontree",
"neutrino_positiontree");
1035 neutrino_positiontree->Branch(
"nnu", &interaction1->
nnu,
"nnu[3]/D");
1036 neutrino_positiontree->Branch(
"dtryingdirection", &interaction1->
dtryingdirection,
"dtryingdirection/D");
1037 neutrino_positiontree->Branch(
"bn1->dtryingposition", &bn1->dtryingposition,
"bn1->dtryingposition/D");
1040 TTree *nupathtree=
new TTree(
"nupathtree",
"nupathtree");
1041 nupathtree->Branch(
"total_kgm2", &total_kgm2,
"total_kgm2/D");
1042 nupathtree->Branch(
"chord", &interaction1->
chord,
"chord/D");
1043 nupathtree->Branch(
"crust_entered", &crust_entered,
"crust_entered/I");
1044 nupathtree->Branch(
"mantle_entered", &mantle_entered,
"mantle_entered/I");
1045 nupathtree->Branch(
"core_entered", &core_entered,
"core_entered/I");
1046 nupathtree->Branch(
"mybeta", &mybeta,
"mybeta/D");
1047 nupathtree->Branch(
"costheta_nutraject", &interaction1->
costheta_nutraject,
"costheta_nutraject/D");
1049 TTree *finaltree =
new TTree(
"passing_events",
"passing_events");
1050 finaltree->Branch(
"inu", &inu,
"inu/I");
1051 finaltree->Branch(
"sample_x", &sample_x,
"sample_x/D");
1052 finaltree->Branch(
"sample_y", &sample_y,
"sample_y/D");
1053 finaltree->Branch(
"vmmhz_min", &vmmhz_min,
"vmmhz_min/D");
1054 finaltree->Branch(
"vmmhz_max", &vmmhz_max,
"vmmhz_max/D");
1055 finaltree->Branch(
"thresholdsAnt", &thresholdsAnt,
"thresholdsAnt[48][2][5]/D");
1056 finaltree->Branch(
"thresholdsAntPass", &thresholdsAntPass,
"thresholdsAntPass[48][2][5]/D");
1057 finaltree->Branch(
"deadTime", &anita1->
deadTime,
"deadTime/D");
1058 finaltree->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
1059 finaltree->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
1060 finaltree->Branch(
"horizcoord_bn", &bn1->horizcoord_bn,
"horizcoord_bn/D");
1061 finaltree->Branch(
"vertcoord_bn", &bn1->vertcoord_bn,
"vertcoord_bn/D");
1062 finaltree->Branch(
"r_bn", &r_bn_array,
"r_bn_array[3]/D");
1063 finaltree->Branch(
"n_bn", &n_bn_array,
"n_bn_array[3]/D");
1064 finaltree->Branch(
"longitude_bn", &longitude_this,
"longitude_bn/D");
1065 finaltree->Branch(
"heading_bn", &heading_this,
"heading_bn/D");
1066 finaltree->Branch(
"gps_offset", &gps_offset,
"gps_offset/D");
1068 finaltree->Branch(
"weight1", &weight1,
"weight1/D");
1070 finaltree->Branch(
"weight", &weight,
"weight/D");
1071 finaltree->Branch(
"logweight", &logweight,
"logweight/D");
1072 finaltree->Branch(
"posnu", &posnu_array,
"posnu_array[3]/D");
1073 finaltree->Branch(
"costheta_nutraject", &costheta_nutraject2,
"costheta_nutraject/D");
1074 finaltree->Branch(
"chord_kgm2_ice", &chord_kgm2_ice2,
"chord_kgm2_ice/D");
1075 finaltree->Branch(
"phi_nutraject", &phi_nutraject2,
"phi_nutraject/D");
1076 finaltree->Branch(
"altitude_int", &altitude_int2,
"altitude_int/D");
1077 finaltree->Branch(
"nnu", &nnu_array,
"nnu_array[3]/D");
1078 finaltree->Branch(
"n_exit2bn", &n_exit2bn_array,
"n_exit2bn_array[5][3]/D");
1079 finaltree->Branch(
"n_exit_phi", &n_exit_phi,
"n_exit_phi/D");
1080 finaltree->Branch(
"rfexit", &rfexit_array,
"rfexit_array[5][3]/D");
1082 finaltree->Branch(
"pnu", &pnu,
"pnu/D");
1083 finaltree->Branch(
"elast_y", &elast_y,
"elast_y/D");
1084 finaltree->Branch(
"emfrac", &emfrac,
"emfrac/D");
1085 finaltree->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1086 finaltree->Branch(
"sumfrac", &sumfrac,
"sumfrac/D");
1087 finaltree->Branch(
"nuflavor", &nuflavorint2,
"nuflavorint/I");
1088 finaltree->Branch(
"current", ¤tint2,
"currentint/I");
1089 finaltree->Branch(
"logchord", &logchord2,
"logchord/D");
1090 finaltree->Branch(
"nuexitice", &nuexitice,
"nuexitice/D");
1091 finaltree->Branch(
"weight_bestcase", &weight_bestcase2,
"weight_bestcase/D");
1092 finaltree->Branch(
"chord_kgm2_bestcase", &chord_kgm2_bestcase2,
"chord_kgm2_bestcase/D");
1093 finaltree->Branch(
"dtryingdirection", &dtryingdirection2,
"dtryingdirection/D");
1094 finaltree->Branch(
"l3trig", &l3trig,
"l3trig[2]/I");
1095 finaltree->Branch(
"l2trig", &l2trig,
"l2trig[2][3]/I");
1096 finaltree->Branch(
"l1trig", &l1trig,
"l1trig[2][3]/I");
1097 finaltree->Branch(
"phiTrigMask", &anita1->phiTrigMask,
"phiTrigMask/s");
1098 finaltree->Branch(
"phiTrigMaskH", &anita1->phiTrigMaskH,
"phiTrigMaskH/s");
1099 finaltree->Branch(
"l1TrigMask", &anita1->l1TrigMask,
"l1TrigMask/s");
1100 finaltree->Branch(
"l1TrigMaskH", &anita1->l1TrigMaskH,
"l1TrigMaskH/s");
1101 finaltree->Branch(
"max_antenna0", &max_antenna0,
"max_antenna0/I");
1102 finaltree->Branch(
"max_antenna1", &max_antenna1,
"max_antenna1/I");
1103 finaltree->Branch(
"max_antenna2", &max_antenna2,
"max_antenna2/I");
1105 finaltree->Branch(
"viewangle", &viewangle,
"viewangle/D");
1106 finaltree->Branch(
"offaxis", &offaxis,
"offaxis/D");
1107 finaltree->Branch(
"rx0_signal_eachband", &rx0_signal_eachband,
"rx0_signal_eachband[2][5]/D");
1108 finaltree->Branch(
"rx0_threshold_eachband", &rx0_threshold_eachband,
"rx0_threshold_eachband[2][5]/D");
1109 finaltree->Branch(
"rx0_noise_eachband", &rx0_noise_eachband,
"rx0_noise_eachband[2][5]/D");
1110 finaltree->Branch(
"rx0_passes_eachband", &rx0_passes_eachband,
"rx0_passes_eachband[2][5]/I");
1111 finaltree->Branch(
"e_component", &e_component,
"e_component[48]/D");
1112 finaltree->Branch(
"h_component", &h_component,
"h_component[48]/D");
1113 finaltree->Branch(
"dist_int_bn_2d", &dist_int_bn_2d,
"dist_int_bn_2d/D");
1114 finaltree->Branch(
"d1", &d12,
"d1/D");
1116 finaltree->Branch(
"cosalpha", &cosalpha,
"cosalpha/D");
1117 finaltree->Branch(
"mytheta", &mytheta,
"mytheta/D");
1118 finaltree->Branch(
"cosbeta0", &cosbeta0,
"cosbeta0/D");
1119 finaltree->Branch(
"mybeta", &mybeta,
"mybeta/D");
1120 finaltree->Branch(
"d1", &d12,
"d1/D");
1121 finaltree->Branch(
"d2", &d22,
"d2/D");
1124 finaltree->Branch(
"fresnel1", &fresnel1,
"fresnel1/D");
1125 finaltree->Branch(
"fresnel2", &fresnel2,
"fresnel2/D");
1126 finaltree->Branch(
"mag1", &mag1,
"mag1/D");
1127 finaltree->Branch(
"mag2", &mag2,
"mag2/D");
1128 finaltree->Branch(
"t_coeff_pokey", &t_coeff_pokey,
"t_coeff_pokey/D");
1129 finaltree->Branch(
"t_coeff_slappy", &t_coeff_slappy,
"t_coeff_slappy/D");
1130 finaltree->Branch(
"exponent", &settings1->EXPONENT,
"EXPONENT/D");
1132 finaltree->Branch(
"hitangle_e_all", &hitangle_e_all,
"hitangle_e_all[48]/D");
1133 finaltree->Branch(
"hitangle_h_all", &hitangle_h_all,
"hitangle_h_all[48]/D");
1135 finaltree->Branch(
"e_comp_max1", &e_comp_max1,
"e_comp_max1/D");
1136 finaltree->Branch(
"h_comp_max1", &h_comp_max1,
"h_comp_max1/D");
1137 finaltree->Branch(
"e_comp_max2", &e_comp_max2,
"e_comp_max2/D");
1138 finaltree->Branch(
"h_comp_max2", &h_comp_max2,
"h_comp_max2/D");
1139 finaltree->Branch(
"e_comp_max3", &e_comp_max3,
"e_comp_max3/D");
1140 finaltree->Branch(
"h_comp_max3", &h_comp_max3,
"h_comp_max3/D");
1141 finaltree->Branch(
"max_antenna_volts0", &max_antenna_volts0,
"max_antenna_volts0/D");
1142 finaltree->Branch(
"max_antenna_volts0_em", &max_antenna_volts0_em,
"max_antenna_volts0_em/D");
1143 finaltree->Branch(
"max_antenna_volts1", &max_antenna_volts1,
"max_antenna_volts1/D");
1144 finaltree->Branch(
"max_antenna_volts2", &max_antenna_volts2,
"max_antenna_volts2/D");
1145 finaltree->Branch(
"triggers", &nchannels_perrx_triggered,
"nchannels_perrx_triggered[48]/I");
1146 finaltree->Branch(
"nchannels_triggered", &nchannels_triggered,
"nchannels_triggered/I");
1147 finaltree->Branch(
"volts_rx_max", &volts_rx_max,
"volts_rx_max/D");
1148 finaltree->Branch(
"volts_rx_ave", &volts_rx_ave,
"volts_rx_ave/D");
1149 finaltree->Branch(
"volts_rx_sum", &volts_rx_sum,
"volts_rx_sum/D");
1151 finaltree->Branch(
"volts_rx_max_highband", &volts_rx_max_highband,
"volts_rx_max_highband/D");
1152 finaltree->Branch(
"volts_rx_max_lowband", &volts_rx_max_lowband,
"volts_rx_max_lowband/D");
1153 finaltree->Branch(
"theta_pol_measured", &theta_pol_measured,
"theta_pol_measured/D");
1154 finaltree->Branch(
"theta_rf_atbn", &theta_rf_atbn,
"theta_rf_atbn/D");
1155 finaltree->Branch(
"theta_rf_atbn_measured", &theta_rf_atbn_measured,
"theta_rf_atbn_measured/D");
1156 finaltree->Branch(
"theta_nnu_atbn",&theta_nnu_atbn,
"theta_nnu_atbn/D");
1157 finaltree->Branch(
"theta_nnu_rf_diff_atbn",&theta_nnu_rf_diff_atbn,
"theta_nnu_rf_diff_atbn/D");
1158 finaltree->Branch(
"phi_nnu_atbn",&phi_nnu_atbn,
"phi_nnu_atbn/D");
1159 finaltree->Branch(
"phi_rf_atbn",&phi_rf_atbn,
"phi_rf_atbn/D");
1160 finaltree->Branch(
"phi_nnu_rf_diff_atbn",&phi_nnu_rf_diff_atbn,
"phi_nnu_rf_diff_atbn/D");
1161 finaltree->Branch(
"voltage", &voltagearray,
"voltagearray[48]/D");
1162 finaltree->Branch(
"nlayers", &settings1->NLAYERS,
"settings1->NLAYERS/I");
1164 finaltree->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
1165 finaltree->Branch(
"vmmhz_lowfreq", &vmmhz_lowfreq,
"vmmhz_lowfreq/D");
1167 finaltree->Branch(
"deltheta_em_max", &deltheta_em_max,
"deltheta_em_max/D");
1168 finaltree->Branch(
"deltheta_had_max", &deltheta_had_max,
"deltheta_had_max/D");
1169 finaltree->Branch(
"r_enterice", &r_enterice_array,
"r_enterice_array[3]/D");
1170 finaltree->Branch(
"n_exit2bn_db", &n_exit2bn_db_array,
"n_exit2bn_db_array[5][3]/D");
1172 finaltree->Branch(
"rfexit_db", &rfexit_db_array,
"rfexit_db_array[5][3]/D");
1173 finaltree->Branch(
"r_in", &r_in_array,
"r_in_array[3]/D");
1174 finaltree->Branch(
"nsurf_rfexit", &nsurf_rfexit_array,
"nsurf_rfexit_array[3]/D");
1175 finaltree->Branch(
"nsurf_rfexit_db", &nsurf_rfexit_db_array,
"nsurf_rfexit_db_array[3]/D");
1176 finaltree->Branch(
"r_fromballoon", &r_fromballoon2,
"r_fromballoon/D");
1177 finaltree->Branch(
"r_fromballoon_db", &interaction1->
r_fromballoon_db,
"r_fromballoon_db/D");
1179 finaltree->Branch(
"nuexitlength", &nuexitlength,
"nuexitlength/D");
1180 finaltree->Branch(
"nuentrancelength", &nuentrancelength,
"nuentrancelength/D");
1181 finaltree->Branch(
"taulength", &taulength,
"taulength/D");
1182 finaltree->Branch(
"icethickness", &icethickness,
"icethickness/D");
1183 finaltree->Branch(
"nrf_iceside", &nrf_iceside_array,
"nrf_iceside_array[5][3]/D");
1184 finaltree->Branch(
"nrf_iceside_db", &nrf_iceside_db_array,
"nrf_iceside_db_array[5][3]/D");
1185 finaltree->Branch(
"ant_normal0", &ant_max_normal0_array,
"ant_max_normal0_array[3]/D");
1186 finaltree->Branch(
"ant_normal1", &ant_max_normal1_array,
"ant_max_normal1_array[3]/D");
1187 finaltree->Branch(
"ant_normal2", &ant_max_normal2_array,
"ant_max_normal2_array[3]/D");
1188 finaltree->Branch(
"vmmhz1m_visible", &vmmhz1m_visible,
"vmmhz1m_visible/D");
1189 finaltree->Branch(
"freq_bins", &freq_bins,
"freq_bins/I");
1190 finaltree->Branch(
"vmmhz", &vmmhz,
"vmmhz[freq_bins]/D");
1192 finaltree->Branch(
"dist_int_bn_2d_chord", &dist_int_bn_2d_chord,
"dist_int_bn_2d_chord/D");
1194 finaltree->Branch(
"dviewangle_deg", &dviewangle_deg,
"dviewangle_deg/D");
1195 finaltree->Branch(
"theta_threshold_deg", &theta_threshold_deg,
"theta_threshold_deg/D");
1196 finaltree->Branch(
"total_kgm2", &total_kgm2,
"total_kgm2/D");
1197 finaltree->Branch(
"chord", &interaction1->
chord,
"chord/D");
1198 finaltree->Branch(
"crust_entered", &crust_entered,
"crust_entered/I");
1199 finaltree->Branch(
"mantle_entered", &mantle_entered,
"mantle_entered/I");
1200 finaltree->Branch(
"core_entered", &core_entered,
"core_entered/I");
1201 finaltree->Branch(
"n_pol", &n_pol_array,
"n_pol_array[3]/D");
1202 finaltree->Branch(
"vmmhz_min_thatpasses", &vmmhz_min_thatpasses,
"vmmhz_min_thatpasses/D");
1204 finaltree->Branch(
"pieceofkm2sr", &pieceofkm2sr,
"pieceofkm2sr/D");
1206 finaltree->Branch(
"r_exit2bn", &r_exit2bn2,
"r_exit2bn/D");
1207 finaltree->Branch(
"r_exit2bn_measured", &r_exit2bn_measured2,
"r_exit2bn_measured/D");
1208 finaltree->Branch(
"scalefactor_attenuation", &scalefactor_attenuation,
"scalefactor_attenuation/D");
1209 finaltree->Branch(
"anita1->PHI_OFFSET", &anita1->
PHI_OFFSET,
"anita1->PHI_OFFSET/D");
1210 finaltree->Branch(
"igps", &bn1->igps,
"igyps/I");
1211 finaltree->Branch(
"volts_rx_rfcm_lab_e_all", &volts_rx_rfcm_lab_e_all,
"volts_rx_rfcm_lab_e_all[48][512]/D");
1212 finaltree->Branch(
"volts_rx_rfcm_lab_h_all", &volts_rx_rfcm_lab_h_all,
"volts_rx_rfcm_lab_h_all[48][512]/D");
1213 finaltree->Branch(
"ptaui", &ptaui,
"ptaui/D");
1214 finaltree->Branch(
"ptauf", &ptauf,
"ptauf/D");
1215 finaltree->Branch(
"sourceLon", &sourceLon,
"sourceLon/D");
1216 finaltree->Branch(
"sourceLat", &sourceLat,
"sourceLat/D");
1217 finaltree->Branch(
"sourceAlt", &sourceAlt,
"sourceAlt/D");
1218 finaltree->Branch(
"sourceMag", &sourceMag,
"sourceMag/D");
1222 TTree *mytaus_tree =
new TTree(
"mytaus",
"mytaus");
1223 mytaus_tree->Branch(
"taus", &TauPtr);
1232 double avgfreq_rfcm[Anita::NFREQ];
1233 double avgfreq_rfcm_lab[Anita::NFREQ];
1234 double freq[Anita::NFREQ];
1236 TTree *summarytree =
new TTree(
"summarytree",
"summarytree");
1237 summarytree->Branch(
"NNU", &NNU,
"NNU/I");
1238 summarytree->Branch(
"EXPONENT", &settings1->EXPONENT,
"EXPONENT/D");
1239 summarytree->Branch(
"eventsfound_beforetrigger", &eventsfound_beforetrigger,
"eventsfound_beforetrigger/D");
1240 summarytree->Branch(
"rms_rfcm_e", &rms_rfcm_e,
"rms_rfcm_e/D");
1241 summarytree->Branch(
"rms_rfcm_h", &rms_rfcm_h,
"rms_rfcm_h/D");
1242 summarytree->Branch(
"rms_lab_e", &rms_lab_e,
"rms_lab_e/D");
1243 summarytree->Branch(
"rms_lab_h", &rms_lab_h,
"rms_lab_h/D");
1244 summarytree->Branch(
"avgfreq_rfcm", &avgfreq_rfcm,
"avgfreq_rfcm[128]/D");
1245 summarytree->Branch(
"avgfreq_rfcm_lab", &avgfreq_rfcm_lab,
"avgfreq_rfcm_lab[128]/D");
1246 summarytree->Branch(
"freq", &freq,
"freq[128]/D");
1248 TTree *banana_tree =
new TTree(
"banana_tree",
"banana_tree");
1249 banana_tree->Branch(
"r_bn", &bn1->r_bn,
"r_bn[3]/D");
1251 TTree *ytree =
new TTree(
"ytree",
"ytree");
1252 ytree->Branch(
"elast_y", &elast_y,
"elast_y/D");
1264 TTree *icetree =
new TTree(
"icetree",
"icetree");
1265 icetree->Branch(
"icethck", &icethck,
"icethck/D");
1266 icetree->Branch(
"lon_ice", &lon_ice,
"lon_ice/D");
1267 icetree->Branch(
"lat_ice", &lat_ice,
"lat_ice/D");
1268 icetree->Branch(
"lon_water", &lon_water,
"lon_water/D");
1269 icetree->Branch(
"lat_water", &lat_water,
"lat_water/D");
1270 icetree->Branch(
"h20_depth", &h20_depth,
"h20_depth/D");
1272 TTree *groundtree =
new TTree(
"groundtree",
"groundtree");
1273 groundtree->Branch(
"elev", &elev,
"elev/D");
1274 groundtree->Branch(
"lon_ground", &lon_ground,
"lon_ground/D");
1275 groundtree->Branch(
"lat_ground", &lat_ground,
"lat_ground/D");
1279 TTree *tree11 =
new TTree(
"h11000",
"h11000");
1280 tree11->Branch(
"loctrig00", &loctrig[0][0],
"loctrig0/D");
1281 tree11->Branch(
"loctrig10", &loctrig[1][0],
"loctrig0/D");
1282 tree11->Branch(
"loctrig20", &loctrig[2][0],
"loctrig0/D");
1283 tree11->Branch(
"loctrig_nadironly0", &loctrig_nadironly[0],
"loctrig_nadironly0/D");
1284 tree11->Branch(
"loctrig01", &loctrig[0][1],
"loctrig1/D");
1285 tree11->Branch(
"loctrig11", &loctrig[1][1],
"loctrig1/D");
1286 tree11->Branch(
"loctrig21", &loctrig[2][1],
"loctrig1/D");
1287 tree11->Branch(
"loctrig_nadironly1", &loctrig_nadironly[1],
"loctrig0/D");
1289 TTree *tree16 =
new TTree(
"h16000",
"h16000");
1290 tree16->Branch(
"pnu", &pnu,
"pnu/D");
1291 tree16->Branch(
"ptau", &ptau,
"ptau/D");
1292 tree16->Branch(
"taulength", &taulength,
"taulength/D");
1293 tree16->Branch(
"weight1", &weight1,
"weight1/D");
1294 tree16->Branch(
"emfrac", &emfrac,
"emfrac/D");
1295 tree16->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1296 tree16->Branch(
"nuentrancelength", &nuentrancelength,
"nuentrancelength/D");
1300 TTree *tree18 =
new TTree(
"h18000",
"h18000");
1301 tree18->Branch(
"emfrac", &emfrac,
"emfrac/D");
1302 tree18->Branch(
"hadfrac", &hadfrac,
"hadfrac/D");
1303 tree18->Branch(
"pdgcode", &pdgcode,
"pdgcode/I");
1306 TH1D *h1mybeta =
new TH1D(
"betaforall",
"betaforall(deg)", 180, -15, 15);
1307 TH1D *h1mytheta=
new TH1D(
"mytheta",
"mytheta(deg)", 180, -90, 90);
1308 TH1F *hundogaintoheight_e=
new TH1F(
"undogaintoheight_e",
"undogaintoheight_e", 100, 0., 1.);
1309 TH1F *hundogaintoheight_h=
new TH1F(
"undogaintoheight_h",
"undogaintoheight_h", 100, 0., 1.);
1310 TH1F *rec_diff=
new TH1F(
"rec_diff",
"rec_diff", 100, -1., 1.);
1311 TH1F *recsum_diff=
new TH1F(
"recsum_diff",
"recsum_diff", 100, -1., 1.);
1312 TH1F *rec_diff0=
new TH1F(
"rec_diff0",
"rec_diff0", 100, -1., 1.);
1313 TH1F *rec_diff1=
new TH1F(
"rec_diff1",
"rec_diff1", 100, -1., 1.);
1314 TH1F *rec_diff2=
new TH1F(
"rec_diff2",
"rec_diff2", 100, -1., 1.);
1315 TH1F *rec_diff3=
new TH1F(
"rec_diff3",
"rec_diff3", 100, -1., 1.);
1317 TTree *vmmhz_tree =
new TTree(
"vmmhz_tree",
"vmmhz_tree");
1318 vmmhz_tree->Branch(
"freq_bins", &freq_bins,
"freq_bins/I");
1319 vmmhz_tree->Branch(
"vmmhz", &vmmhz,
"vmmhz[freq_bins]/D");
1321 TTree *tree1 =
new TTree(
"h1000",
"h1000");
1322 tree1->Branch(
"inu", &inu,
"inu/I");
1323 tree1->Branch(
"diffexit", &diffexit,
"diffexit/D");
1324 tree1->Branch(
"diffrefr", &diffrefr,
"diffrefr/D");
1325 tree1->Branch(
"horizcoord", &horizcoord,
"horizcoord/D");
1326 tree1->Branch(
"vertcoord", &vertcoord,
"vertcoord/D");
1327 tree1->Branch(
"costhetanu", &costhetanu,
"costhetanu/D");
1328 tree1->Branch(
"vmmhz1m_max", &vmmhz1m_max,
"vmmhz1m_max/D");
1329 tree1->Branch(
"volume_thishorizon", &volume_thishorizon,
"volume_thishorizon/D");
1330 tree1->Branch(
"realtime", &realtime_this,
"realtime/D");
1331 tree1->Branch(
"longitude", &longitude_this,
"longitude/D");
1332 tree1->Branch(
"latitude", &latitude_this,
"latitude/D");
1333 tree1->Branch(
"MAXHORIZON", &bn1->MAXHORIZON,
"MAXHORIZON/D");
1334 tree1->Branch(
"igps", &bn1->igps,
"igps/I");
1335 tree1->Branch(
"passes_thisevent", &passes_thisevent,
"passes_thisevent/I");
1336 tree1->Branch(
"igps", &bn1->igps,
"igps/I");
1337 tree1->Branch(
"weight", &weight,
"weight/D");
1338 tree1->Branch(
"r_exit2bn", &interaction1->
r_exit2bn,
"r_exit2bn/D");
1339 tree1->Branch(
"bn1->igps", &bn1->igps,
"bn1->igps/I");
1343 TTree *balloontree =
new TTree(
"balloon",
"balloon");
1344 balloontree->Branch(
"heading", &bn1->heading,
"heading/D");
1345 balloontree->Branch(
"pitch", &bn1->pitch,
"pitch/D");
1346 balloontree->Branch(
"roll", &bn1->roll,
"roll/D");
1347 balloontree->Branch(
"realTime_flightdata", &bn1->realTime_flightdata,
"realTime_flightdata/I");
1348 balloontree->Branch(
"latitude", &bn1->latitude,
"latitude/D");
1349 balloontree->Branch(
"longitude", &bn1->longitude,
"longitude/D");
1350 balloontree->Branch(
"altitude", &bn1->altitude,
"altitude/D");
1351 balloontree->Branch(
"horizcoord_bn", &bn1->horizcoord_bn,
"horizcoord_bn/D");
1352 balloontree->Branch(
"vertcoord_bn", &bn1->vertcoord_bn,
"vertcoord_bn/D");
1355 double undogaintoheight_e=0;
1356 double undogaintoheight_h=0;
1358 double undogaintoheight_e_array[4];
1359 double undogaintoheight_h_array[4];
1361 double nbins_array[4];
1363 double rec_efield=0;
1364 double true_efield=0;
1366 double rec_efield_array[4];
1367 double true_efield_array[4];
1373 IceModel *antarctica =
new IceModel(settings1->ICE_MODEL + settings1->NOFZ*10, settings1->CONSTANTICETHICKNESS * 1000 + settings1->CONSTANTCRUST * 100 + settings1->FIXEDELEVATION * 10 + 0, settings1->WEIGHTABSORPTION);
1374 cout <<
"area of the earth's surface covered by antarctic ice is " << antarctica->ice_area <<
"\n";
1376 ULong_t eventNumber;
1377 #ifdef ANITA_UTIL_EXISTS 1379 string outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaEventFile"+run_num+
".root";
1380 TFile *anitafileEvent =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1382 TTree *eventTree =
new TTree(
"eventTree",
"eventTree");
1383 eventTree->Branch(
"event", &realEvPtr );
1384 eventTree->Branch(
"run", &run_no,
"run/I" );
1385 eventTree->Branch(
"weight", &weight,
"weight/D");
1387 outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaHeadFile"+run_num+
".root";
1388 TFile *anitafileHead =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1390 TTree *headTree =
new TTree(
"headTree",
"headTree");
1391 headTree->Branch(
"header", &rawHeaderPtr );
1392 headTree->Branch(
"weight", &weight,
"weight/D");
1394 outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaGpsFile"+run_num+
".root";
1395 TFile *anitafileGps =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1397 TTree *adu5PatTree =
new TTree(
"adu5PatTree",
"adu5PatTree");
1398 adu5PatTree->Branch(
"pat", &Adu5PatPtr );
1399 adu5PatTree->Branch(
"eventNumber", &eventNumber,
"eventNumber/L");
1400 adu5PatTree->Branch(
"weight", &weight,
"weight/D" );
1402 #ifdef ANITA3_EVENTREADER 1405 AnitaVersion::set(settings1->ANITAVERSION);
1407 outputAnitaFile =string(outputdir.Data())+
"/SimulatedAnitaTruthFile"+run_num+
".root";
1408 TFile *anitafileTruth =
new TFile(outputAnitaFile.c_str(),
"RECREATE");
1411 TString icemcgitversion( EnvironmentVariable::ICEMC_VERSION(outputdir));
1412 printf(
"ICEMC GIT Repository Version: %s\n", icemcgitversion.Data());
1413 unsigned int timenow = time(NULL);
1415 TTree *configAnitaTree =
new TTree(
"configIcemcTree",
"Config file and settings information");
1416 configAnitaTree->Branch(
"gitversion", &icemcgitversion );
1417 configAnitaTree->Branch(
"nnu", &NNU );
1418 configAnitaTree->Branch(
"startTime", &timenow );
1420 configAnitaTree->Fill();
1422 TTree *triggerSettingsTree =
new TTree(
"triggerSettingsTree",
"Trigger settings");
1423 triggerSettingsTree->Branch(
"dioderms", anita1->bwslice_dioderms_fullband_allchan,
"dioderms[2][48][7]/D" );
1424 triggerSettingsTree->Branch(
"diodemean", anita1->bwslice_diodemean_fullband_allchan,
"diodemean[2][48][7]/D");
1425 triggerSettingsTree->Fill();
1428 TTree *summaryAnitaTree =
new TTree(
"summaryAnitaTree",
"summaryAnitaTree");
1429 summaryAnitaTree->Branch(
"EXPONENT", &settings1->EXPONENT,
"EXPONENT/D" );
1430 summaryAnitaTree->Branch(
"SELECTION_MODE", &settings1->UNBIASED_SELECTION,
"SELECTION_MODE/I" );
1431 summaryAnitaTree->Branch(
"totalnu_per_source", &NNU_per_source);
1432 summaryAnitaTree->Branch(
"weighted_nu_per_source", &eventsfound_per_source);
1433 summaryAnitaTree->Branch(
"weighted_prob_nu_per_source", &eventsfound_prob_per_source);
1434 summaryAnitaTree->Branch(
"total_nu", &NNU,
"total_nu/I" );
1435 summaryAnitaTree->Branch(
"total_nue", &count1->nnu_e,
"total_nue/I" );
1436 summaryAnitaTree->Branch(
"total_numu", &count1->nnu_mu,
"total_numu/I" );
1437 summaryAnitaTree->Branch(
"total_nutau", &count1->nnu_tau,
"total_nutau/I" );
1438 summaryAnitaTree->Branch(
"pass_nu", &count1->npass[0],
"pass_nu/I" );
1441 summaryAnitaTree->Branch(
"weighted_nu", &eventsfound,
"weighted_nu/D" );
1442 summaryAnitaTree->Branch(
"weighted_nue", &sum[0],
"weighted_nue/D" );
1443 summaryAnitaTree->Branch(
"weighted_numu", &sum[1],
"weighted_numu/D" );
1444 summaryAnitaTree->Branch(
"weighted_nutau", &sum[2],
"weighted_nutau/D" );
1446 summaryAnitaTree->Branch(
"weighted_prob_nu", &eventsfound_prob,
"weighted_prob_nu/D" );
1447 summaryAnitaTree->Branch(
"weighted_prob_nue", &sum_prob[0],
"weighted_prob_nue/D" );
1448 summaryAnitaTree->Branch(
"weighted_prob_numu", &sum_prob[1],
"weighted_prob_numu/D" );
1449 summaryAnitaTree->Branch(
"weighted_prob_nutau", &sum_prob[2],
"weighted_prob_nutau/D" );
1452 summaryAnitaTree->Branch(
"int_length", &len_int_kgm2,
"int_length/D");
1453 summaryAnitaTree->Branch(
"total_volume", &antarctica->volume,
"total_volume/D");
1454 summaryAnitaTree->Branch(
"rho_medium", &sig1->RHOMEDIUM,
"rho_medium/D");
1455 summaryAnitaTree->Branch(
"selection_box_half_length", &settings1->UNBIASED_PS_MAX_DISTANCE_KM,
"selection_box_half_length/D");
1458 summaryAnitaTree->Branch(
"effVol_nu", &km3sr,
"effVol_nu/D" );
1461 TTree *truthAnitaTree =
new TTree(
"truthAnitaTree",
"Truth Anita Tree");
1462 truthAnitaTree->Branch(
"truth", &truthEvPtr );
1465 TTree* truthAnitaNuTree = settings1->SAVE_TRUTH_NU_TREE ?
new TTree(
"truthAnitaNuTree",
"Truth ANITA Neutrino Tree (all nus)") : 0;
1466 if (settings1->SAVE_TRUTH_NU_TREE)
1468 truthAnitaNuTree->SetAutoFlush(10000);
1469 truthAnitaNuTree->Branch(
"truth_neutrino", &truthNuPtr );
1480 #ifdef ANITA3_EVENTCORRELATOR 1481 const int skyMapLat = 90;
1482 TMarker *astroObject =
new TMarker();
1486 if (bn1->WHICHPATH==3) {
1487 settings1->CONSTANTCRUST = 1;
1488 settings1->CONSTANTICETHICKNESS = 1;
1489 settings1->FIXEDELEVATION = 1;
1493 for (
int i=0;i<antarctica->nRows_ice;i++) {
1494 for (
int j=0;j<antarctica->nCols_ice;j++) {
1495 antarctica->IceENtoLonLat(j, i, lon_ice, lat_ice);
1496 icethck=antarctica->IceThickness(lon_ice, lat_ice);
1499 h20_depth=antarctica->h_water_depth.GetBinContent(j+1,i+1);
1500 if (settings1->HIST) icetree->Fill();
1503 for (
int i=0;i<antarctica->nRows_ground;i++) {
1504 for (
int j=0;j<antarctica->nCols_ground;j++) {
1505 antarctica->GroundENtoLonLat(j, i, lon_ground, lat_ground);
1506 elev=antarctica->SurfaceAboveGeoid(lon_ground, lat_ground);
1507 if (settings1->HIST) groundtree->Fill();
1512 anita1->GetBeamWidths(settings1);
1514 anita1->ReadGains();
1515 anita1->Set_gain_angle(settings1, sig1->NMEDIUM_RECEIVER);
1516 if(settings1->WHICH == 2 || settings1->WHICH == 6) anita1->SetDiffraction();
1518 TCanvas *cgains=
new TCanvas(
"cgains",
"cgains", 880, 800);
1519 TGraph *ggains=
new TGraph(anita1->NPOINTS_GAIN, anita1->frequency_forgain_measured, anita1->vvGaintoHeight);
1521 stemp=string(outputdir.Data())+
"/gains.eps";
1522 cgains->Print((TString)stemp);
1528 bn1->SetDefaultBalloonPosition(antarctica);
1530 anita1->SetNoise(settings1, bn1, antarctica);
1534 antarctica->GetMAXHORIZON(bn1);
1537 antarctica->CreateHorizons(settings1, bn1, bn1->theta_bn, bn1->phi_bn, bn1->altitude_bn, foutput);
1538 cout <<
"Done with CreateHorizons.\n";
1541 if (! src_model && spectra1->IsMonoenergetic() ){
1542 pnu=pow(10., settings1->EXPONENT);
1543 primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, xsecParam_nutype, xsecParam_nuint);
1544 cout <<
"pnu, sigma and len_int_kgm2 (for nu CC) are " << pnu <<
" " << sigma <<
" " << len_int_kgm2 <<
"\n";
1547 if (settings1->WRITEPOSFILE==1) {
1548 nu_out <<
"Neutrinos with energy " << pnu <<
"\n\n";
1549 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";
1552 if (bn1->WHICHPATH==3) {
1553 NNU = settings1->horizontal_banana_points*settings1->vertical_banana_points;
1556 settings1->CONSTANTCRUST=1;
1557 settings1->CONSTANTICETHICKNESS=1;
1558 settings1->FIXEDELEVATION=1;
1559 pnu=Interaction::pnu_banana;
1560 int_banana->surface_over_banana_nu=antarctica->Surface(int_banana->
nu_banana);
1564 cout <<
"reminder that I took out ChangeCoord.\n";
1567 anita1->GetPayload(settings1, bn1);
1569 if (settings1->TRIGGERSCHEME == 3 || settings1->TRIGGERSCHEME == 4 || settings1->TRIGGERSCHEME==5) {
1570 Vector plusz(0., 0., 1.);
1571 bn1->PickBalloonPosition(plusz, antarctica, settings1, anita1);
1572 anita1->calculate_all_offsets();
1573 double angle_theta=16.;
1574 double angle_phi=0.;
1575 Vector x =
Vector(cos(angle_theta * RADDEG) * cos((angle_phi+11.25) * RADDEG),
1576 cos(angle_theta * RADDEG) * sin((angle_phi+11.25) * RADDEG),
1577 sin(angle_theta * RADDEG));
1578 anita1->GetArrivalTimes(x,bn1,settings1);
1579 cout <<
"end of getarrivaltimes\n";
1583 if (settings1->SLAC)
1584 bn1->GetSlacPositions(anita1);
1586 switch (settings1->WHICHRAYS){
1588 settings1->MINRAY=0;
1589 settings1->MAXRAY=0;
1592 settings1->MINRAY=0;
1593 settings1->MAXRAY=1;
1596 settings1->MINRAY=1;
1597 settings1->MAXRAY=1;
1602 time_t raw_loop_start_time = time(NULL);
1603 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;
1606 TCanvas *ctest1 =
new TCanvas(
"ctest1",
"", 880, 800);
1609 if ( spectra1->IsSpectrum() ){
1610 spectra1->GetGEdNdEdAdt()->Draw(
"al");
1611 spectra1->GetGEdNdEdAdt()->GetHistogram()->SetMinimum(-0.2*spectra1->Getmaxflux());
1612 spectra1->GetSEdNdEdAdt()->SetLineColor(2);
1613 spectra1->GetSEdNdEdAdt()->Draw(
"l same");
1615 stemp=string(outputdir.Data())+
"/GetG_test1.pdf";
1616 ctest1->Print((TString)stemp);
1621 double average_altitude=0.;
1622 double average_rbn=0.;
1625 if ( spectra1->IsSpectrum() ){
1626 TCanvas *ctemp =
new TCanvas(
"ctemp");
1628 TFile *out =
new TFile(
"Temp.root",
"recreate");
1629 TGraph *g1 = spectra1->GetGEdNdEdAdt();
1631 ctemp->Print(
"Temp1.png");
1633 double *x = g1->GetX();
1635 for (
int i=0;i<n;i++) x2[i] = TMath::Power(10., x[i]);
1636 TGraph *g2 =
new TGraph(n, x2, g1->GetY());
1640 cout << g2->Integral() <<
" " << g2->Integral(1, n) << endl;
1641 ctemp->Print(
"Temp2.png");
1648 signal(SIGINT, interrupt_signal_handler);
1650 TString objName =
"noObject";
1653 if (settings1->WHICH==7){
1654 gps_offset=atan2(-0.7042,0.71)*DEGRAD;
1655 }
else if(settings1->WHICH==8){
1656 gps_offset=atan2(-0.7085,0.7056)*DEGRAD;
1657 }
else if (settings1->WHICH==9 || settings1->WHICH==10){
1659 }
else gps_offset=0;
1664 for (inu = startNu; inu < NNU; inu++) {
1667 if (inu % (NNU / 100) == 0)
1668 cout << inu <<
" neutrinos. " << (double(inu)/double(NNU)) * 100 <<
"% complete.\n";
1671 cout << inu <<
" neutrinos. " << (double(inu) / double(NNU)) * 100 <<
"% complete.\n";
1673 eventNumber=(ULong_t)(run_no)*NNU+inu;
1678 setSeed(eventNumber);
1684 std::string nunum = Form(
"%d",inu);
1686 double time_weight = 1;
1687 double src_time_weight = 1;
1689 for (whichray = settings1->MINRAY; whichray <= settings1->MAXRAY; whichray++) {
1690 anita1->passglobtrig[0]=0;
1691 anita1->passglobtrig[1]=0;
1693 unmasked_thisevent=1;
1694 vmmhz_min_thatpasses=1000;
1700 int which_source = -1;
1701 if (!src_model && spectra1->IsSpectrum() ){
1703 if(settings1->USEDARTBOARD) pnu=spectra1->GetNuEnergy();
1704 else pnu=spectra1->GetCDFEnergy();
1708 ierr=primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, xsecParam_nutype, interaction1->
currentint);
1711 len_int=1.0/(sigma*sig1->RHOH20*(1./M_NUCL)*1000);
1717 sec1->secondbang=
false;
1721 bn1->dtryingposition=0;
1722 for (
int i=0; i<Anita::NFREQ;i++) {
1727 int got_a_good_position = 0;
1732 bn1->PickBalloonPosition(antarctica, settings1, inu, anita1, getRNG(RNG_BALLOON_POSITION)->Rndm());
1737 average_altitude+=bn1->altitude_bn/(double)NNU;
1738 average_rbn+=bn1->r_bn.Mag()/(double)NNU;
1740 realtime_this=bn1->realTime_flightdata;
1741 longitude_this=bn1->longitude;
1742 latitude_this=bn1->latitude;
1743 altitude_this=bn1->altitude;
1744 heading_this=bn1->heading;
1748 time_weight = src_model->getTimeWeight(realtime_this);
1750 which_source = time_weight > 0 ? src_model->
getDirectionAndEnergy(&force_dir, realtime_this, pnu, src_min, src_max) : -1;
1751 if(which_source >= 0) {
1752 got_a_good_position = 1;
1753 src_time_weight = src_model->getPerSourceTimeWeight(realtime_this, which_source);
1755 if (NNU_per_source.size() < which_source+1)
1757 NNU_per_source.resize(which_source+1);
1758 eventsfound_per_source.resize(which_source+1);
1759 eventsfound_prob_per_source.resize(which_source+1);
1762 NNU_per_source[which_source]++;
1764 const Source * src = src_model->getSource(which_source);
1766 dec= src->getDec()* TMath::RadToDeg();
1767 objName = src->getName();
1769 ierr=primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, xsecParam_nutype, interaction1->
currentint);
1770 len_int=1.0/(sigma*sig1->RHOH20*(1./M_NUCL)*1000);
1781 got_a_good_position = 1;
1783 if (settings1->HORIZON_OFFSET > -999)
1785 double horizon_angle = acos(bn1->r_bn_shadow.Mag()/bn1->r_bn.Mag());
1786 double alt = -horizon_angle + TMath::DegToRad() * settings1->HORIZON_OFFSET;
1788 double lat = bn1->latitude*TMath::DegToRad();
1789 double sin_lat = sin(lat);
1790 double cos_lat = cos(lat);
1791 double sindec = sin_lat *sin(alt) - cos_lat*cos(alt) *cos(Az);
1792 double cosdec = sqrt(1-sindec*sindec);
1793 double h = atan2(sin(Az),-sin_lat*cos(Az) + cos_lat*tan(alt));
1794 force_dir =
Vector(cos(h)*cosdec,sin(h)*cosdec,sindec);
1798 }
while(!got_a_good_position && settings1->SOURCE_SKIP_WHEN_NONE) ;
1800 if (settings1->HIST && !settings1->ONLYFINAL
1801 && prob_eachphi_bn->GetEntries() < settings1->HIST_MAX_ENTRIES) {
1802 prob_eachphi_bn->Fill(bn1->phi_bn);
1803 prob_eachilon_bn->Fill(bn1->r_bn.Lon());
1804 balloontree->Fill();
1807 if (bn1->WHICHPATH==3) {
1809 bn1->setObservationLocation(int_banana, inu, antarctica, settings1);
1813 if (!got_a_good_position)
1815 #ifdef ANITA3_EVENTREADER 1817 if (truthNuPtr) truthNuPtr->setNoNu(bn1->r_bn.GetX(), bn1->r_bn.GetY(), bn1->r_bn.GetZ(), realtime_this);
1832 interaction1 =
new (interaction1)
Interaction(
"nu", primary1, settings1, whichray, count1);
1838 int taumodes = settings1->taumodes;
1844 if (taumodes==1 && interaction1->
nuflavor==
"nutau" && interaction1->
current==
"cc")
1849 bn1->PickDownwardInteractionPoint(interaction1, anita1, settings1, antarctica, ray1, beyondhorizon, len_int_kgm2, (src_model || settings1->HORIZON_OFFSET > -999) ? &force_dir : 0);
1851 bool havent_set_frac =
true;
1852 bool havent_set_weights =
true;
1853 bool havent_set_dir =
true;
1855 #ifdef ANITA3_EVENTREADER 1857 if (settings1->UNBIASED_SELECTION > 0 || src_model || settings1->HORIZON_OFFSET)
1859 havent_set_dir =
false;
1860 if (truthNuPtr) truthNuPtr->setDir(interaction1->
nnu.GetX(), interaction1->
nnu.GetY(), interaction1->
nnu.GetZ());
1866 truthNuPtr->setPos(interaction1->posnu.GetX(), interaction1->posnu.GetY(), interaction1->posnu.GetZ(),
1867 bn1->r_bn.GetX(), bn1->r_bn.GetY(), bn1->r_bn.GetZ(), realtime_this);
1868 truthNuPtr->setNu(pnu,pdgcode);
1872 #ifdef ANITA3_EVENTREADER 1873 #define DO_SKIP { if (truthNuPtr) {truthNuPtr->setSkipped(true,havent_set_frac,havent_set_weights,havent_set_dir); truthAnitaNuTree->Fill();} continue; } 1875 #define DO_SKIP continue ; 1880 if (interaction1->noway)
1884 count1->noway[whichray]++;
1886 if (interaction1->wheredoesitleave_err)
1890 count1->wheredoesitleave_err[whichray]++;
1892 if (interaction1->neverseesice)
1896 count1->neverseesice[whichray]++;
1898 if (interaction1->wheredoesitenterice_err)
1902 count1->wheredoesitenterice_err[whichray]++;
1904 if (interaction1->toohigh)
1908 count1->toohigh[whichray]++;
1910 if (interaction1->toolow)
1914 count1->toolow[whichray]++;
1916 if (bn1->WHICHPATH==3)
1917 interaction1=int_banana;
1923 count1->iceinteraction[whichray]++;
1925 if (beyondhorizon) {
1928 count1->inhorizon[whichray]++;
1931 if (settings1->FIRN) {
1932 sig1->SetNDepth(antarctica->GetN(interaction1->
altitude_int));
1934 changle_deg=sig1->changle*DEGRAD;
1938 if (settings1->FORSECKEL==1)
1939 sig1->SetChangle(acos(1/sig1->NICE));
1942 horizcoord=interaction1->posnu[0]/1000;
1943 vertcoord=interaction1->posnu[1]/1000;
1945 ray1->GetSurfaceNormal(settings1, antarctica, interaction1->posnu, slopeyangle, 0);
1951 costheta_inc=ray1->n_exit2bn[0]*ray1->nsurf_rfexit;
1954 costheta_exit=cos(ray1->rfexit[0].Theta());
1957 if (!ray1->TraceRay(settings1, anita1, 1, sig1->N_DEPTH)) {
1965 ray1->GetRFExit(settings1, anita1, whichray, interaction1->posnu, interaction1->posnu_down, bn1->r_bn, bn1->r_boresights, 1, antarctica);
1967 ray1->GetSurfaceNormal(settings1, antarctica, interaction1->posnu, slopeyangle, 1);
1969 if (!ray1->TraceRay(settings1, anita1, 2, sig1->N_DEPTH)) {;
1976 ray1->GetRFExit(settings1, anita1, whichray, interaction1->posnu, interaction1->posnu_down, bn1->r_bn, bn1->r_boresights, 2, antarctica);
1978 ray1->GetSurfaceNormal(settings1, antarctica, interaction1->posnu, slopeyangle, 2);
1980 if (bn1->WHICHPATH==4)
1981 ray1->PrintAnglesofIncidence();
1984 count1->nraypointsup1[whichray]++;
1986 sec1->GetTauDecay(interaction1->
nuflavor, interaction1->
current, taudecay, emfrac_db, hadfrac_db);
1989 elast_y=primary1->
pickY(settings1, pnu, xsecParam_nutype, interaction1->
currentint);
1990 if (settings1->CONSTANTY==1) {
1995 else if (settings1->CONSTANTY==2)
1997 if (settings1->SHOWERTYPE == 0)
2010 if (bn1->WHICHPATH==3)
2012 if (bn1->WHICHPATH==4)
2014 if (settings1->FORSECKEL==1) {
2015 if (settings1->SHOWERTYPE==0)
2017 if (settings1->SHOWERTYPE==1)
2021 if (ytree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2026 if ( tautrigger == 1 ) {
2027 if ( ( !settings1->UNBIASED_SELECTION ) && ( !settings1->SLAC ) && !src_model ) {
2028 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);
2031 else if (src_model || settings1->HORIZON_OFFSET > -999)
2033 interaction1->
nnu = force_dir;
2035 costhetanu = cos(force_dir.Theta());
2036 theta_threshold = 1;
2039 else if ( settings1->SLAC ) {
2040 Vector xaxis(1., 0., 0.);
2042 interaction1->
nnu = xaxis.RotateY(bn1->theta_bn-settings1->SLAC_HORIZDIST/EarthModel::R_EARTH);
2043 interaction1->
nnu = interaction1->
nnu.RotateZ(bn1->phi_bn);
2044 costhetanu=cos(interaction1->
nnu.Theta());
2046 if (settings1->BORESIGHTS) {
2047 fslac_viewangles << bn1->sslacpositions[bn1->islacposition] <<
"\n";
2048 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2049 for(
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
2050 viewangle_eachboresight[ilayer][ifold]=acos(interaction1->
nnu.Dot(ray1->nrf_iceside_eachboresight[4][ilayer][ifold]));
2051 fslac_viewangles << ilayer <<
"\t" << ifold <<
"\t" << (viewangle_eachboresight[ilayer][ifold]-sig1->changle)*DEGRAD <<
"\n";
2063 interaction1->
r_in = antarctica->WhereDoesItEnter(interaction1->posnu, interaction1->
nnu);
2065 taus1->
GetTauWeight(primary1, settings1, antarctica, interaction1, pnu, 1, ptauf, crust_entered);
2067 antarctica->Getchord(settings1, 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, nearthlayers, myair, total_kgm2, crust_entered, mantle_entered, core_entered);
2072 nutauweightsum +=nutauweight;
2073 tauweightsum +=tauweight;
2074 double xrndm=getRNG(RNG_XRNDM)->Rndm();
2085 TauPtr->
ptauf = ptauf;
2089 if (mytaus_tree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2090 mytaus_tree->Fill();
2098 sec1->GetEMFrac(settings1, interaction1->
nuflavor, interaction1->
current, taudecay, elast_y, hy, pnu, inu,emfrac, hadfrac, n_interactions, tauweighttrigger);
2099 havent_set_frac =
false;
2101 #ifdef ANITA3_EVENTREADER 2102 if (truthNuPtr) truthNuPtr->setFrac(hadfrac,emfrac);
2104 if (emfrac+hadfrac>1.000001) {
2105 cout <<
"Warning: " << inu <<
" " << emfrac+hadfrac <<
"\n";
2110 sumfrac=emfrac+hadfrac;
2112 if (tree7->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2115 vmmhz1m_visible = (emfrac+hadfrac)*vmmhz1m_max;
2118 if (interaction1->
nuflavor==
"numu" && bn1->WHICHPATH != 3 && !settings1->ONLYFINAL && settings1->HIST==1 && fraction_sec_muons->GetEntries()<settings1->HIST_MAX_ENTRIES) {
2119 fraction_sec_muons->Fill(emfrac+hadfrac, weight);
2120 n_sec_muons->Fill((
double)n_interactions);
2123 if (interaction1->
nuflavor==
"nutau" && bn1->WHICHPATH != 3 && !settings1->ONLYFINAL && settings1->HIST==1 && fraction_sec_taus->GetEntries()<settings1->HIST_MAX_ENTRIES) {
2124 fraction_sec_taus->Fill(emfrac+hadfrac, weight);
2125 n_sec_taus->Fill((
double)n_interactions);
2129 if(sec1->secondbang && sec1->interestedintaus) {
2130 ptau=(1-elast_y)*pnu;
2141 if(sec1->secondbang && sec1->interestedintaus) {
2142 vmmhz1m_max=sig1->GetVmMHz1m(ptau, anita1->
FREQ_HIGH);
2143 sig1->GetSpread(ptau, emfrac, hadfrac, anita1->
FREQ_LOW, deltheta_em_max, deltheta_had_max);
2146 vmmhz1m_max=sig1->GetVmMHz1m(pnu, anita1->
FREQ_HIGH);
2147 sig1->GetSpread(pnu, emfrac, hadfrac, anita1->
FREQ_LOW,
2148 deltheta_em_max, deltheta_had_max);
2151 if (jaimetree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2160 bestcase_atten=exp(interaction1->
altitude_int/MAX_ATTENLENGTH);
2165 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) {
2174 count1->nnottoosmall[whichray]++;
2180 chengji=ray1->nrf_iceside[4]*ray1->nrf_iceside[0];
2184 ray1->nrf_iceside[4] = ray1->nrf_iceside[4] - 2*chengji*ray1->nrf_iceside[0];
2187 if ( tautrigger == 0 ) {
2188 if ( ! src_model && ( !settings1->UNBIASED_SELECTION ) && ( !settings1->SLAC ) ) {
2189 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);
2192 else if (src_model || settings1->HORIZON_OFFSET > -999)
2194 interaction1->
nnu = force_dir;
2196 costhetanu = cos(force_dir.Theta());
2197 theta_threshold = 1;
2200 else if (settings1->SLAC) {
2201 Vector xaxis(1., 0., 0.);
2202 interaction1->
nnu = xaxis.RotateY(bn1->theta_bn-settings1->SLAC_HORIZDIST/EarthModel::R_EARTH);
2203 interaction1->
nnu = interaction1->
nnu.RotateZ(bn1->phi_bn);
2204 costhetanu=cos(interaction1->
nnu.Theta());
2206 if (settings1->BORESIGHTS) {
2207 fslac_viewangles << bn1->sslacpositions[bn1->islacposition] <<
"\n";
2208 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2209 for(
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
2210 viewangle_eachboresight[ilayer][ifold]=acos(interaction1->
nnu.Dot(ray1->nrf_iceside_eachboresight[4][ilayer][ifold]));
2211 fslac_viewangles << ilayer <<
"\t" << ifold <<
"\t" << (viewangle_eachboresight[ilayer][ifold]-sig1->changle)*DEGRAD <<
"\n";
2220 #ifdef ANITA3_EVENTREADER 2223 if (truthNuPtr) truthNuPtr->setDir(interaction1->
nnu.GetX(), interaction1->
nnu.GetY(), interaction1->
nnu.GetZ());
2224 havent_set_dir =
false;
2229 viewangle = GetViewAngle(ray1->nrf_iceside[4], interaction1->
nnu);
2230 if(viewangle>1.57 && !settings1->SKIPCUTS) {
2233 count1->nviewangle_lt_90[whichray]++;
2235 if (!Ray::WhereDoesItLeave(interaction1->posnu, interaction1->
nnu, antarctica, interaction1->
nuexit))
2240 GetBalloonLocation(interaction1, ray1, bn1, antarctica);
2246 theta_threshold_deg=theta_threshold*DEGRAD;
2249 n_nutraject_ontheground =
Vector(bn1->n_east*interaction1->
nnu, bn1->n_north*interaction1->
nnu, bn1->n_bn*interaction1->
nnu);
2251 cosviewangle=cos(viewangle);
2252 viewangle_deg=viewangle*DEGRAD;
2253 dviewangle_deg=(sig1->changle-viewangle)*DEGRAD;
2255 if (viewangletree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2256 viewangletree->Fill();
2258 if (neutrino_positiontree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2259 neutrino_positiontree->Fill();
2263 ray1->nrf_iceside[4] = ray1->nrf_iceside[4] + 2*chengji*ray1->nrf_iceside[0];
2267 count1->nbadfracs[whichray]++;
2268 cout<<
"err==0, so leaving.\n";
2271 count1->ngoodfracs[whichray]++;
2274 nsigma_em_threshold=theta_threshold/deltheta_em_max;
2275 nsigma_had_threshold=theta_threshold/deltheta_had_max;
2285 index_distance=(int)(bn1->r_bn.SurfaceDistance(interaction1->posnu, bn1->surface_under_balloon) / (700000./(double)NBINS_DISTANCE));
2289 interaction1->
r_in = antarctica->WhereDoesItEnter(interaction1->posnu, interaction1->
nnu);
2290 antarctica->Getchord(settings1, len_int_kgm2, interaction1->
r_in, interaction1->pathlength_inice, settings1->UNBIASED_SELECTION ==0, interaction1->posnu, inu, interaction1->
chord, interaction1->
weight_nu_prob,
2291 interaction1->
weight_nu, nearthlayers, myair, total_kgm2, crust_entered, mantle_entered, core_entered);
2296 double chord_kgm2_test=interaction1->posnu.
Distance(interaction1->
r_in)*sig1->RHOMEDIUM;
2298 double weight_test=0;
2301 IsAbsorbed(chord_kgm2_test, len_int_kgm2, weight_test);
2309 if ( bn1->WHICHPATH != 4 && settings1->FORSECKEL != 1 && !settings1->SKIPCUTS) {
2310 if (weight_test < settings1->CUTONWEIGHTS) {
2314 count_chanceofsurviving++;
2318 theta_in=interaction1->
r_in.Theta();
2319 lat_in=-90+theta_in*DEGRAD;
2325 myair=GetThisAirColumn(settings1, interaction1->
r_in, interaction1->
nnu, interaction1->posnu, col1, cosalpha, mytheta, cosbeta0, mybeta);
2330 if (!settings1->FORSECKEL && !settings1->UNBIASED_SELECTION) {
2331 if (!antarctica->WhereDoesItEnterIce(interaction1->posnu, interaction1->
nnu, len_int_kgm2/sig1->RHOMEDIUM/10., interaction1->
r_enterice)) {
2333 if (antarctica->OutsideAntarctica(interaction1->
r_enterice)) {
2334 cout<<
"Warning! Neutrino enters beyond continent, program is rejecting neutrino! inu = "<<inu<<endl;
2335 if (bn1->WHICHPATH==3)
2336 cout<<
"Warning! Neutrino enters beyond continent, program is rejecting neutrino!"<<endl;
2343 count1->nentersice[whichray]++;
2354 if(sec1->secondbang && sec1->interestedintaus) {
2355 cout <<
"Need to bring back GetFirstBang before you can simulate taus.\n";
2356 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";
2371 if (tree6->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
2387 if (bn1->WHICHPATH!=4 && interaction1->
weight_bestcase<settings1->CUTONWEIGHTS && !settings1->SKIPCUTS && !settings1->FORSECKEL) {
2388 if (bn1->WHICHPATH==3)
2389 cout<<
"Neutrino is getting absorbed and thrown out!"<<endl;
2394 count1->nabsorbed[whichray]++;
2397 count1->nraywithincontinent1[whichray]++;
2404 bestcase_atten=exp(-1*ray1->rfexit[1].
Distance(interaction1->posnu)/MAX_ATTENLENGTH);
2407 bestcase_atten=exp(-1*ray1->rfexit[1].
Distance(interaction1->posnu_down)/MAX_ATTENLENGTH);
2409 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) {
2410 if (bn1->WHICHPATH==3)
2411 cout<<
"Event rejected. Check."<<endl;
2415 count_chanceinhell0++;
2418 count1->nraypointsup2[whichray]++;
2420 double nbelowsurface;
2422 if (settings1->FIRN)
2423 nbelowsurface=NFIRN;
2425 nbelowsurface=sig1->NICE;
2428 if (!settings1->ROUGHNESS && TIR(ray1->nsurf_rfexit, ray1->nrf_iceside[3], nbelowsurface, sig1->N_AIR)) {
2431 count1->nnottir[whichray]++;
2435 ray1->GetRFExit(settings1, anita1, whichray, interaction1->posnu, interaction1->posnu_down, bn1->r_bn, bn1->r_boresights, 2, antarctica);
2437 count1->nraywithincontinent2[whichray]++;
2442 theta_nnu_atbn = interaction1->
nnu.Angle(bn1->r_bn);
2444 theta_rf_atbn = ray1->n_exit2bn[2].Angle(bn1->r_bn);
2446 theta_rf_atbn_measured = theta_rf_atbn+getRNG(RNG_THETA_RF_RESOLUTION)->Gaus()*anita1->SIGMA_THETA;
2447 interaction1->
r_exit2bn=bn1->r_bn.Distance(ray1->rfexit[2]);
2450 theta_nnu_rf_diff_atbn = theta_nnu_atbn - theta_rf_atbn;
2456 nnu_xy.SetX((interaction1->
nnu[0]) - (bn1->r_bn[0])*((interaction1->
nnu*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2457 nnu_xy.SetY((interaction1->
nnu[1]) - (bn1->r_bn[1])*((interaction1->
nnu*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2458 nnu_xy.SetZ((interaction1->
nnu[2]) - (bn1->r_bn[2])*((interaction1->
nnu*bn1->r_bn)/(bn1->r_bn*bn1->r_bn)));
2460 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)));
2461 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)));
2462 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)));
2464 r_bn_xy.SetX(bn1->r_bn[0]);
2465 r_bn_xy.SetY(bn1->r_bn[1]);
2467 phi_nnu_atbn = nnu_xy.Angle(r_bn_xy);
2468 phi_rf_atbn = rf_xy.Angle(r_bn_xy);
2478 phi_nnu_rf_diff_atbn = phi_nnu_atbn - phi_rf_atbn;
2485 if((settings1->WHICH == 2 || settings1->WHICH == 6) && theta_rf_atbn < 0.3790091) {
2489 if (!antarctica->AcceptableRfexit(ray1->nsurf_rfexit, ray1->rfexit[2], ray1->n_exit2bn[2])){
2490 if (bn1->WHICHPATH==3)
2491 cout<<
"Should look at this. Not expecting to be here."<<endl;
2496 count1->nacceptablerf[whichray]++;
2499 diff_3tries=ray1->rfexit[1].
Distance(ray1->rfexit[2]);
2501 if (tree5->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1 && bn1->WHICHPATH != 3)
2506 if (diff_3tries>10) {
2509 count1->nconverges[whichray]++;
2511 n_pol = GetPolarization(interaction1->
nnu, ray1->nrf_iceside[4]);
2514 if (settings1->BORESIGHTS) {
2515 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2516 for(
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
2517 n_pol_eachboresight[ilayer][ifold]=GetPolarization(interaction1->
nnu, ray1->nrf_iceside_eachboresight[4][ilayer][ifold]);
2523 if (settings1->FIRN){
2526 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);
2527 if (bn1->WHICHPATH==4)
2528 cout <<
"Lentenin factor is " << 1./mag1 <<
"\n";
2533 vmmhz1m_fresneledonce = vmmhz1m_max/mag1;
2536 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);
2539 vmmhz1m_fresneledtwice=vmmhz1m_fresneledonce*fresnel2*mag2;
2541 if (settings1->BORESIGHTS) {
2542 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2543 for(
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
2544 GetFresnel(rough1, settings1->ROUGHNESS, ray1->nsurf_rfexit, ray1->n_exit2bn_eachboresight[2][ilayer][ifold],
2545 n_pol_eachboresight[ilayer][ifold], ray1->nrf_iceside_eachboresight[3][ilayer][ifold],
2546 vmmhz1m_max, emfrac, hadfrac, deltheta_em_max, deltheta_had_max, t_coeff_pokey, t_coeff_slappy,
2547 fresnel1_eachboresight[ilayer][ifold], mag1_eachboresight[ilayer][ifold]);
2554 if (bn1->WHICHPATH==4)
2555 cout<<
"firn-air interface: fresnel2, mag2 are "<<fresnel2<<
" "<< mag2 <<
"\n";
2559 sig1->GetSpread(pnu, emfrac, hadfrac, (anita1->bwslice_min[2]+anita1->bwslice_max[2])/2., deltheta_em_mid2, deltheta_had_mid2);
2561 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);
2563 vmmhz1m_fresneledtwice = vmmhz1m_max*fresnel1*mag1;
2565 if (settings1->BORESIGHTS) {
2566 for(
int ilayer=0;ilayer<settings1->NLAYERS;ilayer++) {
2567 for(
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
2568 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]);
2575 if(settings1->ROUGHNESS){
2583 int num_validscreenpoints = 0;
2585 Vector vec_pos_current_to_balloon;
2589 Vector vec_nnu_to_impactPoint;
2592 double pol_perp_inc, pol_parl_inc;
2593 Vector vec_local_grnd_perp;
2594 Vector vec_local_grnd_parl;
2595 double pol_perp_trans, pol_parl_trans;
2599 double time_reference_specular, time_reference_local;
2600 double pathlength_local;
2601 double viewangle_local;
2602 double azimuth_local;
2604 double theta_0_local;
2605 double tcoeff_perp_polperp, tcoeff_parl_polperp;
2606 double tcoeff_perp_polparl, tcoeff_parl_polparl;
2607 double power_perp_polperp, power_parl_polperp;
2608 double power_perp_polparl, power_parl_polparl;
2609 power_perp_polperp = power_parl_polperp = power_perp_polparl = power_parl_polparl = 0.;
2610 double fresnel_r, mag_r;
2612 Vector npol_local_inc, npol_local_trans;
2620 if (settings1->FIRN)
2621 time_reference_specular = (interaction1->posnu.
Distance(ray1->rfexit[2])*NFIRN / CLIGHT) + (ray1->rfexit[2].
Distance(bn1->r_bn)/CLIGHT);
2623 time_reference_specular = (interaction1->posnu.
Distance(ray1->rfexit[2])*NICE / CLIGHT) + (ray1->rfexit[2].
Distance(bn1->r_bn)/CLIGHT);
2625 double slopeyx, slopeyy, slopeyz, rtemp;
2631 double basescreenedgelength = settings1->SCREENEDGELENGTH;
2632 double grd_stepsize = settings1->SCREENSTEPSIZE;
2634 if(settings1->ROUGHSIZE>0)
2635 grd_nsteps = int(basescreenedgelength/2. / grd_stepsize);
2650 vec_localnormal = antarctica->GetSurfaceNormal(ray1->rfexit[2]).Unit();
2654 panel1->
SetUnitY( (vec_localnormal.Cross(ray1->n_exit2bn[2])).Unit() );
2663 Emag_local = vmmhz1m_max;
2664 taperfactor = fresnel_r = mag_r = 1.;
2665 tcoeff_perp_polparl = tcoeff_parl_polparl = 0.;
2666 tcoeff_perp_polperp = tcoeff_parl_polperp = 0.;
2667 pos_projectedImpactPoint = panel1->
GetPosition(ii, jj);
2668 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] );
2671 vec_localnormal = antarctica->GetSurfaceNormal(pos_projectedImpactPoint).Unit();
2672 if (settings1->SLOPEY) {
2673 slopeyx=ray1->slopeyx;
2674 slopeyy=ray1->slopeyy;
2675 slopeyz=ray1->slopeyz;
2676 ntemp2 = vec_localnormal + slopeyx*xaxis + slopeyy*yaxis + slopeyz*zaxis;
2677 ntemp2 = ntemp2 / ntemp2.Mag();
2678 rtemp= ntemp2 * vec_localnormal;
2680 vec_localnormal = ntemp2;
2684 vec_nnu_to_impactPoint =
Vector( pos_projectedImpactPoint[0]-interaction1->posnu[0], pos_projectedImpactPoint[1]-interaction1->posnu[1], pos_projectedImpactPoint[2]-interaction1->posnu[2] ).Unit();
2686 vec_grndcomp2IP = (vec_nnu_to_impactPoint - (vec_nnu_to_impactPoint.Dot(vec_localnormal)*vec_localnormal)).Unit();
2687 vec_grndcomp2bln = (vec_pos_current_to_balloon - (vec_pos_current_to_balloon.Dot(vec_localnormal)*vec_localnormal)).Unit();
2688 temp_a = vec_localnormal.Cross(vec_pos_current_to_balloon).Unit();
2689 azimuth_local = vec_grndcomp2IP.Angle(vec_grndcomp2bln);
2690 if( temp_a.Dot(vec_nnu_to_impactPoint) < 0 )
2691 azimuth_local *= -1.;
2696 theta_local = vec_localnormal.Angle( (
const Vector)vec_pos_current_to_balloon );
2697 theta_0_local = vec_localnormal.Angle(vec_nnu_to_impactPoint);
2700 if( isnan(theta_local) | isnan(theta_0_local) | isnan(azimuth_local) ){
2703 viewangle_local = GetViewAngle(vec_nnu_to_impactPoint, interaction1->
nnu);
2706 deltheta_em[0]=deltheta_em_max*anita1->
FREQ_LOW/anita1->freq[0];
2707 deltheta_had[0]=deltheta_had_max*anita1->
FREQ_LOW/anita1->freq[0];
2708 sig1->TaperVmMHz(viewangle_local, deltheta_em[0], deltheta_had[0], emfrac, hadfrac, taperfactor, vmmhz_em[0]);
2717 if (settings1->FIRN)
2718 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);
2720 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);
2723 if( (power_perp_polperp==0.)&(power_parl_polperp==0.)&(power_perp_polparl==0.)&(power_parl_polparl==0.) ){
2727 if (settings1->FIRN){
2728 tcoeff_perp_polparl = sqrt(power_perp_polparl);
2729 tcoeff_parl_polparl = sqrt(power_parl_polparl);
2730 tcoeff_perp_polperp = sqrt(power_perp_polperp);
2731 tcoeff_parl_polperp = sqrt(power_parl_polperp);
2734 tcoeff_perp_polparl = sqrt(power_perp_polparl);
2735 tcoeff_parl_polparl = sqrt(power_parl_polparl);
2736 tcoeff_perp_polperp = sqrt(power_perp_polperp);
2737 tcoeff_parl_polperp = sqrt(power_parl_polperp);
2745 pathlength_local = interaction1->posnu.
Distance(pos_projectedImpactPoint) + pos_projectedImpactPoint.
Distance(bn1->r_bn);
2747 Emag_local /= pathlength_local ;
2749 Attenuate(antarctica, settings1, Emag_local, interaction1->posnu.
Distance(pos_projectedImpactPoint), interaction1->posnu);
2754 npol_local_inc = GetPolarization(interaction1->
nnu, vec_nnu_to_impactPoint).Unit();
2755 vec_inc_perp = (vec_localnormal.Cross(vec_nnu_to_impactPoint)).Unit();
2756 vec_inc_parl = (vec_nnu_to_impactPoint.Cross(vec_inc_perp)).Unit();
2757 pol_perp_inc = npol_local_inc * vec_inc_perp;
2758 pol_parl_inc = npol_local_inc * vec_inc_parl;
2760 pol_perp_trans = pol_perp_inc * tcoeff_perp_polperp + pol_perp_inc * tcoeff_perp_polparl;
2761 pol_parl_trans = pol_parl_inc * tcoeff_parl_polparl + pol_parl_inc * tcoeff_parl_polperp;
2763 vec_local_grnd_perp = (vec_localnormal.Cross(vec_pos_current_to_balloon)).Unit();
2764 vec_local_grnd_parl = (vec_pos_current_to_balloon.Cross(vec_local_grnd_perp)).Unit();
2767 npol_local_trans= (pol_perp_trans*vec_local_grnd_perp + pol_parl_trans*vec_local_grnd_parl).Unit();
2770 if( isnan(npol_local_trans[0]) ){
2775 fresnel_r = sqrt( pow(vmmhz1m_max*pol_perp_trans,2) + pow(vmmhz1m_max*pol_parl_trans,2) ) / vmmhz1m_max;
2776 mag_r = sqrt( tan(theta_0_local) / tan(theta_local) );
2777 Emag_local *= fresnel_r * mag_r;
2779 if (settings1->FIRN)
2780 time_reference_local = (interaction1->posnu.
Distance(pos_projectedImpactPoint)*NFIRN / CLIGHT) + (pos_projectedImpactPoint.
Distance(bn1->r_bn)/CLIGHT);
2782 time_reference_local = (interaction1->posnu.
Distance(pos_projectedImpactPoint)*NICE / CLIGHT) + (pos_projectedImpactPoint.
Distance(bn1->r_bn)/CLIGHT);
2784 num_validscreenpoints++;
2788 panel1->
AddVec2bln(vec_pos_current_to_balloon);
2789 panel1->
AddPol(npol_local_trans);
2790 panel1->
AddDelay( time_reference_specular - time_reference_local );
2814 double validScreenSummedArea = 0.;
2815 double vmmhz_local_array[Anita::NFREQ];
2818 sig1->GetVmMHz(panel1->
GetVmmhz0(jj), vmmhz1m_max, pnu, anita1->freq, anita1->
NOTCH_MIN, anita1->NOTCH_MAX, vmmhz_local_array, Anita::NFREQ);
2820 for (
int k=0;k<Anita::NFREQ;k++) {
2821 deltheta_em[k]=deltheta_em_max*anita1->
FREQ_LOW/anita1->freq[k];
2822 deltheta_had[k]=deltheta_had_max*anita1->
FREQ_LOW/anita1->freq[k];
2823 sig1->TaperVmMHz(panel1->
GetViewangle(jj), deltheta_em[k], deltheta_had[k], emfrac, hadfrac, vmmhz_local_array[k], vmmhz_em[k]);
2827 validScreenSummedArea += panel1->
GetWeight(jj);
2832 vmmhz_max = Tools::dMax(vmmhz_max, panel1->
GetVmmhz_freq(jj*Anita::NFREQ));
2845 if(settings1->CHANCEINHELL_FACTOR*vmmhz1m_fresneledtwice*heff_max*0.5*(anita1->bwmin/1.E6)<anita1->maxthreshold*anita1->VNOISE[0]/10.&& !settings1->SKIPCUTS) {
2846 if (bn1->WHICHPATH==3)
2847 cout<<
"Event is undetectable. Leaving loop."<<endl;
2851 count1->nchanceinhell_fresnel[whichray]++;
2856 diffexit=ray1->rfexit[0].
Distance(ray1->rfexit[1]);
2857 diffnorm=acos(ray1->nsurf_rfexit[0]*ray1->nsurf_rfexit[1]);
2858 diffrefr=acos(ray1->nrf_iceside[4]*ray1->nrf_iceside[0]);
2862 if (!settings1->ROUGHNESS) {
2866 vmmhz_max=ScaleVmMHz(vmmhz1m_fresneledtwice, interaction1->posnu, bn1->r_bn, ray1->rfexit[2]);
2868 vmmhz_max=ScaleVmMHz(vmmhz1m_fresneledtwice, interaction1->posnu_down, bn1->r_bn, ray1->rfexit[2]);
2877 if (!settings1->ROUGHNESS){
2878 if (settings1->CHANCEINHELL_FACTOR*vmmhz_max*heff_max*0.5*(anita1->bwmin/1.E6)<anita1->maxthreshold*anita1->VNOISE[0]/10. && !settings1->SKIPCUTS) {
2879 if (bn1->WHICHPATH==3)
2880 cout<<
"Event is undetectable. Leaving loop."<<endl;
2885 count1->nchanceinhell_1overr[whichray]++;
2888 if (!settings1->ROUGHNESS) {
2890 rflength=interaction1->posnu.
Distance(ray1->rfexit[2]);
2893 rflength=interaction1->posnu_down.
Distance(ray1->rfexit[2]);
2897 if (bn1->WHICHPATH==4)
2898 cout <<
"rflength is " << rflength <<
"\n";
2901 if (!settings1->ROUGHNESS) {
2903 Attenuate(antarctica, settings1, vmmhz_max, rflength, interaction1->posnu);
2905 Attenuate_down(antarctica, settings1, vmmhz_max, ray1->rfexit[2], interaction1->posnu, interaction1->posnu_down);
2912 if (tree2->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1 && bn1->WHICHPATH != 3)
2919 if (!settings1->ROUGHNESS){
2920 if (settings1->CHANCEINHELL_FACTOR*vmmhz_max*heff_max*0.5*(anita1->bwmin/1.E6)<anita1->maxthreshold*anita1->VNOISE[0]/10. && !settings1->SKIPCUTS) {
2921 if (bn1->WHICHPATH==3)
2922 cout<<
"Event is undetectable. Leaving loop."<<endl;
2928 count1->nchanceinhell[whichray]++;
2931 if (tree3->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1 && bn1->WHICHPATH != 3)
2944 if (!settings1->ROUGHNESS){
2945 if (settings1->FORSECKEL==1)
2946 sig1->SetNDepth(sig1->NICE);
2948 sig1->GetVmMHz(vmmhz_max, vmmhz1m_max, pnu, anita1->freq, anita1->
NOTCH_MIN, anita1->NOTCH_MAX, vmmhz, Anita::NFREQ);
2958 undogaintoheight_e=0;
2959 undogaintoheight_h=0;
2961 for (
int k=0;k<4;k++) {
2962 undogaintoheight_e_array[k]=0.;
2963 undogaintoheight_h_array[k]=0.;
2965 true_efield_array[k]=0.;
2966 rec_efield_array[k]=0.;
2973 if (!settings1->ROUGHNESS){
2975 double rtemp=Tools::dMin((viewangle-sig1->changle)/(deltheta_em_max), (viewangle-sig1->changle)/(deltheta_had_max));
2976 if (rtemp>Signal::VIEWANGLE_CUT && !settings1->SKIPCUTS) {
2980 count1->nviewanglecut[whichray]++;
2985 for (
int k=0;k<Anita::NFREQ;k++) {
2986 deltheta_em[k]=deltheta_em_max*anita1->
FREQ_LOW/anita1->freq[k];
2987 deltheta_had[k]=deltheta_had_max*anita1->
FREQ_LOW/anita1->freq[k];
2989 if (settings1->FORSECKEL==1) {
2990 for (
int iviewangle=0;iviewangle<NVIEWANGLE;iviewangle++) {
2992 vmmhz_temp=vmmhz[k]*vmmhz1m_max/vmmhz_max;
2994 viewangle_temp=viewangles[iviewangle];
2997 sig1->TaperVmMHz(viewangle_temp, deltheta_em[k], deltheta_had[k], emfrac, hadfrac, vmmhz_temp, djunk);
2998 forseckel[iviewangle][k]=vmmhz_temp;
3007 sig1->TaperVmMHz(viewangle, deltheta_em[k], deltheta_had[k], emfrac, hadfrac, vmmhz[k], vmmhz_em[k]);
3016 vmmhz_lowfreq=vmmhz[0];
3022 if (sig1->logscalefactor_taper>maxtaper)
3023 maxtaper=sig1->logscalefactor_taper;
3025 if (interaction1->
nuflavor==
"nue") pdgcode = 12;
3026 else if (interaction1->
nuflavor==
"numu") pdgcode = 14;
3027 else if (interaction1->
nuflavor==
"nutau") pdgcode = 16;
3029 if (settings1->HIST==1 && !settings1->ONLYFINAL && bn1->WHICHPATH != 3 && k==Anita::NFREQ/2 && tree18->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3034 if (bn1->WHICHPATH == 3)
3035 interaction1->
banana_volts += vmmhz[k]*(settings1->BW/(double)Anita::NFREQ/1.E6);
3043 if (bn1->WHICHPATH==3 && interaction1->
banana_volts != 0 && settings1->HIST && banana_tree->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3044 banana_tree->Fill();
3047 else if (bn1->WHICHPATH==3 && interaction1->
banana_volts == 0) {
3059 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) {
3068 if(!settings1->ROUGHNESS){
3069 vmmhz_max=Tools::dMax(vmmhz, Anita::NFREQ);
3070 vmmhz_min=Tools::dMin(vmmhz, Anita::NFREQ);
3073 count1->nchanceinhell2[whichray]++;
3080 if (settings1->USEDEADTIME){
3081 if ( (getRNG(RNG_DEADTIME)->Uniform(1)<anita1->
deadTime) ){
3083 if (settings1->MINBIAS!=1) {
3089 count1->ndeadtime[whichray]++;
3091 Tools::Zero(sumsignal_aftertaper, 5);
3102 if(!settings1->ROUGHNESS){
3104 for (
int k=0;k<Anita::NFREQ;k++) {
3114 if(settings1->FIRN){
3130 for (
int k=0;k<Anita::NFREQ;k++) {
3131 if (bn1->WHICHPATH==4)
3132 IntegrateBands(anita1, k, panel1, anita1->freq, vmmhz1m_max/(vmmhz_max*1.E6), sumsignal_aftertaper);
3141 if (!settings1->TRIGGEREFFSCAN){
3142 if(settings1->BORESIGHTS)
3143 anita1->GetArrivalTimesBoresights(ray1->n_exit2bn_eachboresight[2]);
3145 anita1->GetArrivalTimes(ray1->n_exit2bn[2],bn1,settings1);
3147 anita1->rx_minarrivaltime=Tools::WhichIsMin(anita1->arrival_times[0], settings1->NANTENNAS);
3150 for (
int i=0;i<settings1->NANTENNAS;i++) {
3158 max_antenna_volts0 =0;
3159 max_antenna_volts1 =0;
3160 max_antenna_volts2 =0;
3173 if (settings1->SLAC)
3174 fslac_hitangles << bn1->sslacpositions[bn1->islacposition] <<
"\n";
3177 double rotateangle=getRNG(RNG_RANDOMISE_POL)->Gaus(RANDOMISEPOL*RADDEG);
3178 n_pol=n_pol.Rotate(rotateangle, ray1->n_exit2bn[2]);
3181 if (bn1->WHICHPATH==4) {
3182 Tools::Zero(sumsignal, 5);
3183 for (
int k=0;k<Anita::NFREQ;k++)
3184 IntegrateBands(anita1, k, panel1, anita1->freq, bn1->r_bn.Distance(interaction1->posnu)/1.E6, sumsignal);
3188 bn1->CenterPayload(hitangle_e);
3191 if (settings1->MAKEVERTICAL) {
3195 Vector rotationaxis=ray1->n_exit2bn[2].Cross(bn1->n_bn);
3196 double rotateangle=PI/2.-ray1->n_exit2bn[2].Dot(bn1->n_bn);
3197 ray1->n_exit2bn[2]=ray1->n_exit2bn[2].Rotate(rotateangle, rotationaxis);
3199 for (
int ilayer=0; ilayer < settings1->NLAYERS; ilayer++) {
3201 for (
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
3202 Vector rotationaxis2=ray1->n_exit2bn_eachboresight[2][ilayer][ifold].Cross(n_pol_eachboresight[ilayer][ifold]);
3203 double rotateangle2=PI/2.-ray1->n_exit2bn_eachboresight[2][ilayer][ifold].Dot(n_pol_eachboresight[ilayer][ifold]);
3204 ray1->n_exit2bn_eachboresight[2][ilayer][ifold].Rotate(rotateangle2, rotationaxis2);
3209 globaltrig1->volts_rx_rfcm_trigger.assign(16, vector <vector <double> >(3, vector <double>(0)));
3212 if (!settings1->BORESIGHTS) {
3213 bn1->GetEcompHcompkvector(n_eplane, n_hplane, n_normal, ray1->n_exit2bn[2], e_component_kvector[0], h_component_kvector[0], n_component_kvector[0]);
3214 bn1->GetEcompHcompEvector(settings1, n_eplane, n_hplane, n_pol, e_component[0], h_component[0], n_component[0]);
3217 for (
int ilayer=0; ilayer < settings1->NLAYERS; ilayer++) {
3218 for (
int ifold=0;ifold<anita1->
NRX_PHI[ilayer];ifold++) {
3223 bn1->GetAntennaOrientation(settings1, anita1, ilayer, ifold, n_eplane, n_hplane, n_normal);
3225 if (settings1->BORESIGHTS){
3226 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]);
3227 bn1->GetEcompHcompEvector(settings1, n_eplane, n_hplane, n_pol_eachboresight[ilayer][ifold], e_component[count_rx], h_component[count_rx], n_component[count_rx]);
3228 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";
3230 else if (count_rx > 0)
3232 e_component_kvector[count_rx] = e_component_kvector[0];
3233 h_component_kvector[count_rx] = h_component_kvector[0];
3234 n_component_kvector[count_rx] = n_component_kvector[0];
3235 e_component[count_rx] = e_component[0];
3236 h_component[count_rx] = h_component[0];
3237 n_component[count_rx] = n_component[0];
3240 bn1->GetHitAngles(e_component_kvector[count_rx], h_component_kvector[count_rx], n_component_kvector[count_rx], hitangle_e, hitangle_h);
3242 hitangle_h_all[count_rx]=hitangle_h;
3243 hitangle_e_all[count_rx]=hitangle_e;
3245 if (h6->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3246 h6->Fill(hitangle_h, ray1->n_exit2bn[2]*bn1->n_bn);
3250 chantrig1->
ApplyAntennaGain(settings1, anita1, bn1, panel1, antNum, n_eplane, n_hplane, n_normal);
3252 chantrig1->
TriggerPath(settings1, anita1, antNum, bn1);
3280 chantrig1->
TimeShiftAndSignalFluct(settings1, anita1, ilayer, ifold, volts_rx_rfcm_lab_e_all, volts_rx_rfcm_lab_h_all);
3282 chantrig1->
saveTriggerWaveforms(anita1, justSignal_trig[0][antNum], justSignal_trig[1][antNum], justNoise_trig[0][antNum], justNoise_trig[1][antNum]);
3283 chantrig1->
saveDigitizerWaveforms(anita1, justSignal_dig[0][antNum], justSignal_dig[1][antNum], justNoise_dig[0][antNum], justNoise_dig[1][antNum]);
3285 Tools::Zero(sumsignal, 5);
3287 if (bn1->WHICHPATH==4 && ilayer==anita1->GetLayer(anita1->rx_minarrivaltime) && ifold==anita1->GetIfold(anita1->rx_minarrivaltime)) {
3288 for (
int ibw=0;ibw<5;ibw++) {
3289 cout <<
"Just after Taper, sumsignal is " << sumsignal_aftertaper[ibw] <<
"\n";
3290 cout <<
"Just after antennagain, sumsignal is " << sumsignal[ibw] <<
"\n";
3295 if (count_rx==anita1->rx_minarrivaltime) {
3296 undogaintoheight_e/=(double)Anita::NFREQ;
3297 undogaintoheight_h/=(double)Anita::NFREQ;
3298 for (
int k=0;k<4;k++) {
3299 undogaintoheight_e_array[k]/=(double)nbins_array[k];
3300 undogaintoheight_h_array[k]/=(double)nbins_array[k];
3303 if (settings1->SCALEDOWNLCPRX1)
3304 globaltrig1->volts[0][ilayer][0]=globaltrig1->volts[0][ilayer][0]/sqrt(2.);
3306 if (settings1->RCPRX2ZERO)
3307 globaltrig1->volts[1][ilayer][1]=0.;
3309 if (settings1->LCPRX2ZERO)
3310 globaltrig1->volts[0][ilayer][1]=0.;
3312 if (settings1->SIGNAL_FLUCT) {
3313 if (settings1->WHICH==0) {
3314 globaltrig1->volts[ilayer][ifold][0]+=getRNG(RNG_SIGNAL_FLUCT)->Gaus(0., anita1->VNOISE_ANITALITE[ifold]);
3315 globaltrig1->volts[ilayer][ifold][1]+=getRNG(RNG_SIGNAL_FLUCT)->Gaus(0., anita1->VNOISE_ANITALITE[ifold]);
3318 if (count_rx==anita1->rx_minarrivaltime) {
3319 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));
3320 for (
int ibw=0;ibw<4;ibw++) {
3321 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));
3322 bwslice_vnoise_thislayer[ibw]=anita1->bwslice_vnoise[ilayer][ibw];
3325 if (tree6b->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3331 chantrig1->
WhichBandsPass(settings1, anita1, globaltrig1, bn1, ilayer, ifold, viewangle-sig1->changle, emfrac, hadfrac, thresholdsAnt[antNum]);
3334 if (Anita::GetAntennaNumber(ilayer, ifold)==anita1->rx_minarrivaltime) {
3335 for (
int iband=0;iband<5;iband++) {
3336 for (
int ipol=0;ipol<2;ipol++) {
3337 rx0_signal_eachband[ipol][iband]=chantrig1->
signal_eachband[ipol][iband];
3339 rx0_noise_eachband[ipol][iband]=chantrig1->
noise_eachband[ipol][iband];
3340 rx0_passes_eachband[ipol][iband]=chantrig1->
passes_eachband[ipol][iband];
3346 if (ilayer == 0 && globaltrig1->volts[0][ilayer][ifold] > max_antenna_volts0) {
3347 max_antenna0 = count_rx;
3348 max_antenna_volts0 = globaltrig1->volts[0][ilayer][ifold];
3349 max_antenna_volts0_em=globaltrig1->volts_em[0][ilayer][ifold];
3350 ant_max_normal0 = ant_normal;
3351 e_comp_max1 = e_component[count_rx];
3352 h_comp_max1 = h_component[count_rx];
3354 else if (ilayer == 0 && globaltrig1->volts[0][ilayer][ifold] == max_antenna_volts0 && globaltrig1->volts[0][ilayer][ifold] != 0){
3355 cout<<
"Equal voltage on two antennas! Event : "<<inu<<endl;
3357 else if (ilayer == 1 && globaltrig1->volts[0][ilayer][ifold] > max_antenna_volts1) {
3358 max_antenna1 = count_rx;
3359 max_antenna_volts1 = globaltrig1->volts[0][ilayer][ifold];
3360 ant_max_normal1 = ant_normal;
3361 e_comp_max2 = e_component[count_rx];
3362 h_comp_max2 = h_component[count_rx];
3364 else if (ilayer == 1 && globaltrig1->volts[0][ilayer][ifold] == max_antenna_volts1 && globaltrig1->volts[0][ilayer][ifold] != 0){
3365 cout<<
"Equal voltage on two antennas! Event : "<<inu<<endl;
3367 else if (ilayer == 2 && globaltrig1->volts[0][ilayer][ifold] > max_antenna_volts2) {
3368 max_antenna2 = count_rx;
3369 max_antenna_volts2 = globaltrig1->volts[0][ilayer][ifold];
3370 ant_max_normal2 = ant_normal;
3371 e_comp_max3 = e_component[count_rx];
3372 h_comp_max3 = h_component[count_rx];
3374 else if (ilayer == 2 && globaltrig1->volts[0][ilayer][ifold] == max_antenna_volts2 && globaltrig1->volts[0][ilayer][ifold] != 0){
3375 cout<<
"Equal voltage on two antennas! Event : "<<inu<<endl;
3377 voltagearray[count_rx] = globaltrig1->volts[0][ilayer][ifold];
3383 if (settings1->TRIGTYPE==0 && ifold==1 && count_pass>=settings1->NFOLD) {
3384 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;
3393 if(!settings1->ROUGHNESS){
3394 if (settings1->DISCONES==1) {
3396 for (
int idiscone=0;NDISCONES;idiscone++) {
3399 polarfactor_discone=n_pol.Dot(bn1->n_bn);
3400 for (
int k=0;k<Anita::NFREQ;k++) {
3401 if (anita1->freq[k]>=FREQ_LOW_DISCONES && anita1->freq[k]<=FREQ_HIGH_DISCONES) {
3402 thislambda=CLIGHT/sig1->N_AIR/anita1->freq[k];
3403 heff_discone= thislambda*sqrt(2*Zr*gain_dipole/Z0/4/PI*sig1->N_AIR);
3405 volts_discone+=panel1->
GetVmmhz_freq(k)*0.5*heff_discone*((settings1->BW/1E6)/(
double)Anita::NFREQ)*polarfactor_discone;
3409 vnoise_discone=anita1->VNOISE[0]*sqrt(BW_DISCONES/settings1->BW_SEAVEYS);
3411 if (settings1->SIGNAL_FLUCT) {
3412 volts_discone+=getRNG(RNG_SIGNAL_FLUCT)->Gaus(0., vnoise_discone);
3415 if (fabs(volts_discone)/vnoise_discone>anita1->maxthreshold)
3422 for (
int irx=0;irx<settings1->NANTENNAS;irx++) {
3423 nchannels_perrx_triggered[irx]=globaltrig1->nchannels_perrx_triggered[irx];
3426 nchannels_triggered=Tools::iSum(globaltrig1->nchannels_perrx_triggered, settings1->NANTENNAS);
3427 volts_rx_ave=GetAverageVoltageFromAntennasHit(settings1, globaltrig1->nchannels_perrx_triggered, voltagearray, volts_rx_sum);
3432 if(sec1->secondbang && sec1->interestedintaus)
3433 count_asktrigger_nfb++;
3436 dist_int_bn_2d_chord = ray1->rfexit[0].
Distance(bn1->r_bn_shadow)/1000;
3438 dist_int_bn_2d = ray1->rfexit[0].
SurfaceDistance(bn1->r_bn_shadow, bn1->surface_under_balloon) / 1000;
3444 if (!antarctica->Getchord(settings1, len_int_kgm2, interaction1->
r_in, interaction1->pathlength_inice, settings1->UNBIASED_SELECTION ==0,
3445 interaction1->posnu, inu, interaction1->
chord, interaction1->
weight_nu_prob, interaction1->
weight_nu, nearthlayers, myair, total_kgm2, crust_entered, mantle_entered, core_entered)){
3449 if(tauweighttrigger==1) {
3461 weight = weight1 / interaction1->
dnutries * settings1->SIGMA_FACTOR * time_weight;
3462 weight_prob = weight_prob / interaction1->
dnutries * settings1->SIGMA_FACTOR * time_weight;
3470 #ifdef ANITA3_EVENTREADER 3471 if (truthNuPtr) truthNuPtr->setWeights(weight, 1./interaction1->
dnutries, time_weight);
3473 havent_set_weights =
false;
3476 if (weight < settings1->CUTONWEIGHTS || weight_prob < settings1->CUTONWEIGHTPROBS) {
3481 eventsfound_beforetrigger+=weight;
3493 globaltrig1->
PassesTrigger(settings1, anita1, discones_passing, 2, l3trignoise, l2trig, l1trig, settings1->antennaclump, loctrig, loctrig_nadironly, inu,
3494 thispassesnoise,
true);
3496 globaltrig1->
PassesTrigger(settings1, anita1, discones_passing, 2, l3trig, l2trig, l1trig, settings1->antennaclump, loctrig, loctrig_nadironly, inu,
3499 if ( (l3trignoise[0]>0 ) || (l3trignoise[1]>0 ) ){
3500 cout <<
"A thermal noise fluctuation generated this trigger!" << l3trignoise[0] <<
" " << l3trig[0] <<
" " << l3trignoise[1] <<
" " << l3trig[1] << endl;
3506 for (
int i=0;i<2;i++) {
3507 for (
int j=0;j<16;j++) {
3508 for (
int k=0;k<anita1->HALFNFOUR;k++) {
3509 count1->nl1triggers[i][whichray]+=anita1->l1trig_anita3and4_inanita[i][j][k];
3523 if ( (thispasses[0]==1 && anita1->pol_allowed[0]==1)
3524 || (thispasses[1]==1 && anita1->pol_allowed[1]==1)
3525 || (settings1->TRIGTYPE==0 && count_pass>=settings1->NFOLD)
3526 || (settings1->MINBIAS==1)){
3528 if (bn1->WHICHPATH==4)
3529 cout <<
"This event passes.\n";
3531 anita1->passglobtrig[0]=thispasses[0];
3532 anita1->passglobtrig[1]=thispasses[1];
3535 n_exit_phi = Tools::AbbyPhiCalc(ray1->n_exit2bn[2][0], ray1->n_exit2bn[2][1]);
3539 count1->npassestrigger[whichray]++;
3544 if (tree11->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3548 if(sec1->secondbang && sec1->interestedintaus)
3549 count_passestrigger_nfb++;
3559 if (nupathtree->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3563 count_chordgoodlength++;
3565 pieceofkm2sr=weight*antarctica->volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr/(double)NNU/len_int;
3566 if (h10->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST)
3567 h10->Fill(hitangle_e_all[0], weight);
3578 logweight=log10(weight);
3582 if (interaction1->
d2>1) {
3584 count_d2goodlength++;
3588 if(sec1->secondbang && sec1->interestedintaus) {
3589 eventsfound_nfb+=weight;
3590 index_weights=(int)(((logweight-MIN_LOGWEIGHT)/(MAX_LOGWEIGHT-MIN_LOGWEIGHT))*(double)NBINS);
3591 eventsfound_nfb_binned[index_weights]++;
3592 if (tree16->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3596 allcuts[whichray]++;
3597 allcuts_weighted[whichray]+=weight;
3598 if (thispasses[0] && thispasses[1]) {
3599 allcuts_weighted_polarization[2]+=weight;
3600 }
else if (thispasses[0]){
3601 allcuts_weighted_polarization[0]+=weight;
3602 }
else if (thispasses[1]){
3603 allcuts_weighted_polarization[1]+=weight;
3605 anita1->weight_inanita=weight;
3607 if (h1mybeta->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1)
3608 h1mybeta -> Fill(mybeta, weight);
3610 eventsfound+=weight;
3611 eventsfound_prob+=weight_prob;
3613 eventsfound_per_source[which_source]+=weight*src_time_weight/time_weight;
3614 eventsfound_prob_per_source[which_source]+=weight_prob*src_time_weight/time_weight;
3618 eventsfound_belowhorizon+=weight;
3620 count1->npass[whichray]++;
3624 if (logweight<MIN_LOGWEIGHT)
3626 else if (logweight>MAX_LOGWEIGHT)
3627 index_weights=NBINS-1;
3629 index_weights=(int)(((logweight-MIN_LOGWEIGHT)/(MAX_LOGWEIGHT-MIN_LOGWEIGHT))*(double)NBINS);
3632 if (index_weights<NBINS)
3633 eventsfound_binned[index_weights]++;
3636 if (index_distance<NBINS_DISTANCE)
3637 eventsfound_binned_distance[index_distance]+= weight;
3640 if (index_distance<NBINS_DISTANCE && index_weights<NBINS)
3641 eventsfound_binned_distance_forerror[index_distance][index_weights]++;
3644 eventsfound_weightgt01+=weight;
3647 if (nearthlayers==1)
3648 eventsfound_crust+=weight;
3650 if (h1mybeta->GetEntries()<settings1->HIST_MAX_ENTRIES && !settings1->ONLYFINAL && settings1->HIST==1) {
3651 h1mybeta -> Fill(mybeta, weight);
3652 h1mytheta -> Fill(mytheta, weight);
3658 double int_lon, int_lat;
3661 int_lon = interaction1->posnu.
Lon();
3662 int_lat = interaction1->posnu.
Lat();
3663 antarctica->LonLattoEN(int_lon, int_lat, E, N);
3665 dir_int_coord->Fill(E, N);
3667 ref_int_coord->Fill(E, N);
3671 offaxis=(double)fabs(viewangle-sig1->changle);
3672 nsigma_offaxis=offaxis/deltheta_had_max;
3674 hundogaintoheight_e->Fill(undogaintoheight_e, weight);
3675 hundogaintoheight_h->Fill(undogaintoheight_h, weight);
3676 rec_diff->Fill((rec_efield-true_efield)/true_efield, weight);
3677 rec_diff0->Fill((rec_efield_array[0]-true_efield_array[0])/true_efield_array[0], weight);
3678 rec_diff1->Fill((rec_efield_array[1]-true_efield_array[1])/true_efield_array[1], weight);
3679 rec_diff2->Fill((rec_efield_array[2]-true_efield_array[2])/true_efield_array[2], weight);
3680 rec_diff3->Fill((rec_efield_array[3]-true_efield_array[3])/true_efield_array[3], weight);
3681 recsum_diff->Fill((rec_efield_array[0]+rec_efield_array[1]+rec_efield_array[2]+rec_efield_array[3]-true_efield)/true_efield, weight);
3684 sourceLon = ray1->rfexit[2].
Lon() - 180;
3685 sourceLat = ray1->rfexit[2].
Lat() - 90;
3686 sourceAlt = antarctica->SurfaceAboveGeoid(sourceLon+180, sourceLat+90);
3689 if (settings1->HIST && finaltree->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3690 for (
int i=0;i<3;i++) {
3691 nnu_array[i] = interaction1->
nnu[i];
3692 r_in_array[i] = interaction1->
r_in[i];
3693 r_bn_array[i] = bn1->r_bn[i];
3694 n_bn_array[i] = bn1->n_bn[i];
3695 posnu_array[i] = interaction1->posnu[i];
3696 ant_max_normal0_array[i] = ant_max_normal0[i];
3697 ant_max_normal1_array[i] = ant_max_normal1[i];
3698 ant_max_normal2_array[i] = ant_max_normal2[i];
3699 n_pol_array[i] = n_pol[i];
3700 r_enterice_array[i] = interaction1->
r_enterice[i];
3701 nsurf_rfexit_array[i] = ray1->nsurf_rfexit[i];
3702 nsurf_rfexit_db_array[i] = ray1->nsurf_rfexit_db[i];
3704 for (
int j=0;j<5;j++) {
3705 for (
int i=0;i<3;i++) {
3706 nrf_iceside_array[j][i] = ray1->nrf_iceside[j][i];
3707 nrf_iceside_db_array[j][i] = nrf_iceside_db[j][i];
3708 n_exit2bn_array[j][i] = ray1->n_exit2bn[j][i];
3709 n_exit2bn_db_array[j][i] = n_exit2bn_db[j][i];
3710 rfexit_array[j][i] = ray1->rfexit[j][i];
3711 rfexit_db_array[j][i] = ray1->rfexit_db[j][i];
3714 if (settings1->HIST && vmmhz_tree->GetEntries()<20) {
3723 d12=interaction1->
d1;
3724 d22=interaction1->
d2;
3734 sourceMag = ray1->rfexit[2].Mag();
3735 sample_x = antarctica->getSampleX();
3736 sample_y = antarctica->getSampleY();
3739 count1->IncrementWeights_r_in(interaction1->
r_in, weight);
3742 #ifdef ANITA_UTIL_EXISTS 3746 Adu5PatPtr->
latitude= bn1->latitude;
3748 Adu5PatPtr->
altitude=bn1->altitude;
3749 Adu5PatPtr->
realTime=bn1->realTime_flightdata;
3750 Adu5PatPtr->
heading = bn1->heading;
3751 Adu5PatPtr->
pitch = bn1->pitch;
3752 Adu5PatPtr->roll = bn1->roll;
3753 Adu5PatPtr->run = run_no;
3756 memset(realEvPtr->
fVolts, 0,
sizeof(realEvPtr->
fVolts) );
3757 memset(realEvPtr->
fTimes, 0,
sizeof(realEvPtr->
fTimes) );
3759 int fNumPoints = 260;
3760 for (
int ichan=0; ichan<108; ichan++){
3763 for (
int j = 0; j < fNumPoints; j++) {
3765 realEvPtr->
fTimes[ichan][j] = j * anita1->TIMESTEP * 1.0E9;
3769 for (
int iant = 0; iant < settings1->NANTENNAS; iant++){
3771 int IceMCAnt = GetIceMCAntfromUsefulEventAnt(settings1, iant);
3776 realEvPtr->
chanId[UsefulChanIndexV] = UsefulChanIndexV;
3777 realEvPtr->
chanId[UsefulChanIndexH] = UsefulChanIndexH;
3779 for (
int j = 0; j < fNumPoints; j++) {
3784 realEvPtr->
fVolts[UsefulChanIndexH][j] = volts_rx_rfcm_lab_h_all[IceMCAnt][j+128]*1000;
3786 realEvPtr->
fVolts[UsefulChanIndexV][j] = volts_rx_rfcm_lab_e_all[IceMCAnt][j+128]*1000;
3797 if (settings1->MINBIAS==1)
3803 rawHeaderPtr->
run = run_no;
3813 if (settings1->WHICH<9){
3814 rawHeaderPtr->
phiTrigMask = (short) anita1->phiTrigMask;
3819 rawHeaderPtr->
realTime = bn1->realTime_flightdata;
3820 rawHeaderPtr->
triggerTime = bn1->realTime_flightdata;
3821 Adu5PatPtr->
latitude= bn1->latitude;
3823 Adu5PatPtr->
altitude=bn1->altitude;
3824 Adu5PatPtr->
realTime=bn1->realTime_flightdata;
3825 Adu5PatPtr->
heading = bn1->heading;
3826 Adu5PatPtr->
pitch = bn1->pitch;
3827 Adu5PatPtr->roll = bn1->roll;
3828 Adu5PatPtr->run = run_no;
3831 cout <<
"Neutrino (evNum = " << eventNumber <<
", weight = " << weight <<
" weight_prob = " << weight_prob <<
") passed" << endl;
3833 #ifdef ANITA3_EVENTREADER 3834 if (settings1->WHICH==9 || settings1->WHICH==10) {
3837 rawHeaderPtr->setMask( (
short) anita1->l1TrigMask, (
short) anita1->phiTrigMask,
AnitaPol::kVertical);
3838 rawHeaderPtr->setMask( (
short) anita1->l1TrigMaskH, (
short) anita1->phiTrigMaskH,
AnitaPol::kHorizontal);
3843 truthEvPtr->
realTime = bn1->realTime_flightdata;
3844 truthEvPtr->
run = run_no;
3845 truthEvPtr->
nuMom = pnu;
3846 truthEvPtr->
showerE = pnu * sumfrac;
3847 truthEvPtr->
nu_pdg = pdgcode;
3858 if(src_model){cout <<
"It originated from " << objName <<
" (" << RA*15 <<
"°, " << dec <<
"°)" << endl;}
3860 truthEvPtr->
weight = weight;
3861 truthEvPtr->
weight1 = weight1;
3866 truthEvPtr->
tuffIndex = (short)anita1->tuffIndex;
3870 #ifdef ANITA3_EVENTCORRELATOR
3871 if(settings1->ALL_SKY_MAP)
3873 astroObject->SetX(RA);
3874 astroObject->SetY(dec);
3875 astroObject->SetMarkerStyle(20);
3876 astroObject->SetMarkerSize(1.5);
3877 astroObject->SetMarkerColor(9);
3878 skyMapOut->addMarker(astroObject);
3884 truthEvPtr->RA = RA;
3885 truthEvPtr->
dec = dec;
3890 for (
int i=0;i<3;i++){
3893 truthEvPtr->
nuPos[i] = interaction1->posnu[i];
3894 truthEvPtr->
nuDir[i] = interaction1->
nnu[i];
3896 for (
int i=0;i<5;i++){
3897 for (
int j=0;j<3;j++){
3898 truthEvPtr->
rfExitNor[i][j] = ray1->n_exit2bn[i][j];
3899 truthEvPtr->
rfExitPos[i][j] = ray1->rfexit[i][j];
3902 for (
int i=0;i<48;i++){
3903 truthEvPtr->
hitangle_e[i] = hitangle_e_all[i];
3904 truthEvPtr->
hitangle_h[i] = hitangle_h_all[i];
3906 if(!settings1->ROUGHNESS){
3907 for (
int i=0;i<Anita::NFREQ;i++)
3924 for (
int iant = 0; iant < settings1->NANTENNAS; iant++){
3928 truthEvPtr->
SNRAtTrigger[UsefulChanIndexV] = Tools::calculateSNR(justSignal_trig[0][iant], justNoise_trig[0][iant]);
3929 truthEvPtr->
SNRAtTrigger[UsefulChanIndexH] = Tools::calculateSNR(justSignal_trig[1][iant], justNoise_trig[1][iant]);
3934 truthEvPtr->
SNRAtDigitizer[UsefulChanIndexV] = Tools::calculateSNR(justSignal_dig[0][iant], justNoise_dig[0][iant]);
3935 truthEvPtr->
SNRAtDigitizer[UsefulChanIndexH] = Tools::calculateSNR(justSignal_dig[1][iant], justNoise_dig[1][iant]);
3941 truthEvPtr->
thresholds[UsefulChanIndexV] = thresholdsAnt[iant][0][4];
3942 truthEvPtr->
thresholds[UsefulChanIndexH] = thresholdsAnt[iant][1][4];
3945 if (iant%2) irx = iant/2;
3946 else irx = iant/2 + 1;
3949 for (
int j = 0; j < fNumPoints; j++) {
3950 truthEvPtr->
fTimes[UsefulChanIndexV][j] = j * anita1->TIMESTEP * 1.0E9;
3951 truthEvPtr->
fTimes[UsefulChanIndexH][j] = j * anita1->TIMESTEP * 1.0E9;
3953 truthEvPtr->
fSignalAtTrigger[UsefulChanIndexV][j] = justSignal_trig[0][iant][j+128]*1000;
3954 truthEvPtr->
fSignalAtTrigger[UsefulChanIndexH][j] = justSignal_trig[1][iant][j+128]*1000;
3955 truthEvPtr->
fNoiseAtTrigger[UsefulChanIndexV][j] = justNoise_trig[0][iant][j+128]*1000;
3956 truthEvPtr->
fNoiseAtTrigger[UsefulChanIndexH][j] = justNoise_trig[1][iant][j+128]*1000;
3957 truthEvPtr->
fSignalAtDigitizer[UsefulChanIndexV][j] = justSignal_dig[0][iant][j+128]*1000;
3958 truthEvPtr->
fSignalAtDigitizer[UsefulChanIndexH][j] = justSignal_dig[1][iant][j+128]*1000;
3959 truthEvPtr->
fNoiseAtDigitizer[UsefulChanIndexV][j] = justNoise_dig[0][iant][j+128]*1000;
3960 truthEvPtr->
fNoiseAtDigitizer[UsefulChanIndexH][j] = justNoise_dig[1][iant][j+128]*1000;
3968 if (truthNuPtr) truthNuPtr->setSkipped(
false);
3972 if (truthAnitaNuTree) truthAnitaNuTree->Fill();
3973 truthAnitaTree->Fill();
3979 adu5PatTree->Fill();
3982 delete rawHeaderPtr;
3986 sum_weights+=weight;
3987 neutrinos_passing_all_cuts++;
3988 times_crust_entered_det+=crust_entered;
3989 times_mantle_entered_det+=mantle_entered;
3990 times_core_entered_det+=core_entered;
3992 if (settings1->WRITEPOSFILE==1)
3993 WriteNeutrinoInfo(interaction1->posnu, interaction1->
nnu, bn1->r_bn, interaction1->
altitude_int, interaction1->
nuflavor, interaction1->
current, elast_y, nu_out);
3996 if (settings1->HIST && !settings1->ONLYFINAL && sampleweights->GetEntries()<settings1->HIST_MAX_ENTRIES) {
3998 sampleweights->Fill(log10(weight));
4000 sampleweights->Fill(-6.);
4003 if (sampleweights->GetEntries()==1000) {
4004 double sum_sampleintegral=0.;
4005 double sum_sample=0.;
4007 for (
int k=sampleweights->GetNbinsX();k>=1;k--) {
4008 sum_sampleintegral+=sampleweights->GetBinContent(k)*pow(10., sampleweights->GetBinLowEdge(k));
4011 sum_sampleintegral+=sampleweights->GetBinContent(0)*pow(10., sampleweights->GetBinLowEdge(1));
4013 for (
int k=sampleweights->GetNbinsX();k>=1;k--) {
4014 sum_sample+=sampleweights->GetBinContent(k)*pow(10., sampleweights->GetBinLowEdge(k));
4015 if (sum_sample>0.99*sum_sampleintegral) {
4028 forbrian << interaction1->
costheta_nutraject <<
" " << n_nutraject_ontheground.Phi() <<
" " << bn1->phi_bn <<
" " << logweight <<
"\n";
4031 if (interaction1->
nuflavor==
"nue") {
4033 sum_prob[0]+=weight_prob;
4034 eventsfound_binned_e[index_weights]++;
4036 if (interaction1->
nuflavor==
"numu") {
4038 sum_prob[1]+=weight_prob;
4039 eventsfound_binned_mu[index_weights]++;
4041 if(!sec1->secondbang || !sec1->interestedintaus) {
4042 if (interaction1->
nuflavor==
"nutau") {
4044 sum_prob[2]+=weight_prob;
4045 eventsfound_binned_tau[index_weights]++;
4053 cout <<
"Chord is less than 1m.\n";
4056 if (settings1->HIST==1 && !settings1->ONLYFINAL && anita1->tglob->GetEntries()<settings1->HIST_MAX_ENTRIES) {
4058 anita1->tglob->Fill();
4059 anita1->
tdata->Fill();
4073 if (bn1->WHICHPATH==4)
4074 cout <<
"This event does not pass.\n";
4111 count_chanceinhell2_w += weight;
4114 count_passestrigger_w += weight;
4116 volume_thishorizon=antarctica->volume_inhorizon[bn1->Getibnposition()]/1.E9;
4118 if (settings1->HIST==1
4119 && !settings1->ONLYFINAL
4120 && tree1->GetEntries()<settings1->HIST_MAX_ENTRIES
4121 && bn1->WHICHPATH != 3){
4128 std::cout <<
"\n***********************************************************";
4129 std::cout <<
"\n* SIGINT received, aborting loop over events early.";
4130 std::cout <<
"\n* Stopped after event " << inu <<
" instead of " << NNU;
4131 std::cout <<
"\n* Any output which relied on NNU should be corrected for.";
4132 std::cout <<
"\n***********************************************************\n";
4133 foutput <<
"\n***********************************************************";
4134 foutput <<
"\n* SIGINT received, aborting loop over events early.";
4135 foutput <<
"\n* Stopped after event " << inu <<
" instead of " << NNU;
4136 foutput <<
"\n* Any output which relied on NNU should be corrected for.";
4137 foutput <<
"\n***********************************************************\n";
4145 cout <<
"about to close tsignals tree.\n";
4146 anita1->fsignals=anita1->tsignals->GetCurrentFile();
4147 anita1->fdata=anita1->
tdata->GetCurrentFile();
4149 anita1->fsignals->Write();
4150 anita1->fsignals->Close();
4152 anita1->fdata=anita1->tglob->GetCurrentFile();
4153 anita1->fdata->Write();
4154 anita1->fdata->Close();
4156 if (settings1->EVENTSMAP){
4158 TH2F *lat80deg=
new TH2F(
"lat80deg",
"", 600, -3000, 3000, 500, -2500, 2500);
4159 lat80deg->SetMarkerColor(kRed);
4161 for(
double lon=0;lon<360.;lon+=0.5){
4163 antarctica->LonLattoEN(lon, lat, E, N);
4164 lat80deg->Fill(E, N);
4168 if (bn1->WHICHPATH==4) {
4169 cout <<
"Earth radius at South Pole: " << antarctica->Geoid(0.) <<
"\n";
4171 cout <<
"Surface of ice at the South Pole: " << antarctica->Surface(posnu_temp) <<
"\n";
4172 cout <<
"Average balloon altitude is " << average_altitude <<
"\n";
4173 cout <<
"Average distance from earth center is " << average_rbn <<
"\n";
4174 cout <<
"Average height of balloon above ice surface is " << average_rbn-antarctica->Surface(bn1->r_bn) <<
"\n";
4176 cout <<
"Index of refraction at this depth is " << sig1->N_DEPTH <<
"\n";
4177 cout <<
"Cerenkov angle is " << sig1->changle*DEGRAD <<
"\n";
4178 cout <<
"Nadir angle to surface exit point is " << DEGRAD*bn1->r_bn.Angle(ray1->n_exit2bn[2]) <<
"\n";
4179 cout <<
"Distance from rfexit to balloon is " << (bn1->r_bn+ -1*ray1->rfexit[2]).Mag() <<
"\n";
4180 cout <<
"Payload zenith angle at event source is " << DEGRAD*ray1->rfexit[2].Angle(ray1->n_exit2bn[2]) <<
"\n";
4181 cout <<
"Angle of incidence just below surface is " << DEGRAD*ray1->rfexit[2].Angle(ray1->nrf_iceside[4]) <<
"\n";
4182 cout <<
"Angle between neutrino and surface normal is " << DEGRAD*ray1->rfexit[0].Angle(interaction1->
nnu)-90. <<
"\n";
4183 cout <<
"Angle of incidence below firn surface is " << DEGRAD*ray1->rfexit[2].Angle(ray1->nrf_iceside[3]) <<
"\n";
4185 cout <<
"about to Summarize.\n";
4187 anita1->
rms_rfcm[0] = sqrt(anita1->
rms_rfcm[0] / (
double)anita1->count_getnoisewaveforms)*1000.;
4188 anita1->
rms_rfcm[1] = sqrt(anita1->
rms_rfcm[1] / (
double)anita1->count_getnoisewaveforms)*1000.;
4189 anita1->
rms_lab[0] = sqrt(anita1->
rms_lab[0] / (
double)anita1->count_getnoisewaveforms)*1000.;
4190 anita1->
rms_lab[1] = sqrt(anita1->
rms_lab[1] / (
double)anita1->count_getnoisewaveforms)*1000.;
4192 cout <<
"RMS noise in rfcm e-pol is " << anita1->
rms_rfcm[0] <<
" mV.\n";
4193 cout <<
"RMS noise in rfcm h-pol is " << anita1->
rms_rfcm[1] <<
" mV.\n";
4194 cout <<
"RMS noise in lab e-pol is " << anita1->
rms_lab[0] <<
"mV.\n";
4195 cout <<
"RMS noise in lab h-pol is " << anita1->
rms_lab[1] <<
"mV.\n";
4196 for (
int i=0;i<Anita::NFREQ;i++) {
4197 anita1->avgfreq_rfcm[i]/=(double)anita1->count_getnoisewaveforms;
4198 anita1->avgfreq_rfcm_lab[i]/=(
double)anita1->count_getnoisewaveforms;
4205 for (
int i=0;i<Anita::NFREQ;i++) {
4206 avgfreq_rfcm[i]=anita1->avgfreq_rfcm[i];
4207 avgfreq_rfcm_lab[i]=anita1->avgfreq_rfcm_lab[i];
4208 freq[i]=anita1->freq[i];
4210 cout <<
"Filling summarytree. rms_rfcm_e is " << rms_rfcm_e <<
"\n";
4211 summarytree->Fill();
4215 primary1->
GetSigma(pnu, sigma, len_int_kgm2, settings1, xsecParam_nutype, xsecParam_nuint);
4216 len_int=1.0/(sigma*sig1->RHOH20*(1./M_NUCL)*1000);
4219 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);
4221 std::string spec_string = src_model ? settings1->SOURCE : std::to_string( settings1->EXPONENT);
4222 veff_out << spec_string <<
"\t" << km3sr <<
"\t" << km3sr_e <<
"\t" << km3sr_mu <<
"\t" << km3sr_tau <<
"\t" << settings1->SIGMA_FACTOR << endl;
4225 sum_frac[0]=sum[0]/eventsfound;
4226 sum_frac[1]=sum[1]/eventsfound;
4227 sum_frac[2]=sum[2]/eventsfound;
4230 sum_frac_db[0]=sum[0]/(eventsfound+eventsfound_db+eventsfound_nfb);
4231 sum_frac_db[1]=sum[1]/(eventsfound+eventsfound_db+eventsfound_nfb);
4232 sum_frac_db[2]=(sum[2]+eventsfound_db+eventsfound_nfb)/(eventsfound+eventsfound_db+eventsfound_nfb);
4239 #ifdef ANITA_UTIL_EXISTS 4241 anitafileEvent->cd();
4242 eventTree->Write(
"eventTree");
4243 anitafileEvent->Close();
4244 delete anitafileEvent;
4246 anitafileHead->cd();
4247 headTree->Write(
"headTree");
4248 anitafileHead->Close();
4249 delete anitafileHead;
4252 adu5PatTree->Write(
"adu5PatTree");
4253 anitafileGps->Close();
4254 delete anitafileGps;
4256 #ifdef ANITA3_EVENTREADER 4257 anitafileTruth->cd();
4258 configAnitaTree->Write(
"configAnitaTree");
4259 truthAnitaTree->Write(
"truthAnitaTree");
4260 if (truthAnitaNuTree) truthAnitaNuTree->Write();
4261 triggerSettingsTree->Write(
"triggerSettingsTree");
4262 summaryAnitaTree->Fill();
4263 summaryAnitaTree->Write(
"summaryAnitaTree");
4264 anitafileTruth->Close();
4265 delete anitafileTruth;
4271 #ifdef ANITA3_EVENTCORRELATOR 4272 if(settings1->ALL_SKY_MAP)
4274 TCanvas *skyMapC =
new TCanvas(
"skyMapC",
"skyMapC", 1000, 500);
4277 string outputSkyMap =string(outputdir.Data())+
"/allSkyMap"+run_num+
".png";
4278 skyMapC->SaveAs(outputSkyMap.c_str());
4282 cout <<
"closing file.\n";
4285 time_t raw_end_time = time(NULL);
4286 struct tm * end_time = localtime(&raw_end_time);
4287 cout <<
"Date and time at end of run are: " << asctime (end_time) <<
"\n";
4288 cout<<
"Total time elapsed is "<<(int)((raw_end_time - raw_start_time)/60)<<
":"<< ((raw_end_time - raw_start_time)%60)<<endl;
4290 foutput <<
"\nTotal time elapsed in run is " <<(int)((raw_end_time - raw_start_time)/60)<<
":"<< ((raw_end_time - raw_start_time)%60)<<endl;
4316 void IntegrateBands(
Anita *anita1,
int k,
Screen *panel1,
double *freq,
double scalefactor,
double *sumsignal) {
4317 for (
int j=0;j<5;j++) {
4320 if (anita1->bwslice_min[j]<=freq[k] && anita1->bwslice_max[j]>freq[k])
4321 sumsignal[j]+=panel1->
GetVmmhz_freq(jpt*Anita::NFREQ + k)*(freq[k+1]-freq[k])*scalefactor;
4328 void Integrate(
Anita *anita1,
int j,
int k,
double *vmmhz,
double *freq,
double scalefactor,
double sumsignal) {
4330 if (anita1->bwslice_min[j]<=freq[k] && anita1->bwslice_max[j]>freq[k])
4331 sumsignal+=vmmhz[k]*(freq[k+1]-freq[k])*scalefactor;
4336 void WriteNeutrinoInfo(
Position &posnu,
Vector &nnu,
Position &r_bn,
double altitude_int,
string nuflavor,
string current,
double elast_y, ofstream &nu_out) {
4337 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";
4342 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) {
4343 double rate_v_thresh[NTHRESHOLDS];
4344 double errorup_v_thresh[NTHRESHOLDS];
4345 double errordown_v_thresh[NTHRESHOLDS];
4346 double rate_h_thresh[NTHRESHOLDS];
4347 double errorup_h_thresh[NTHRESHOLDS];
4348 double errordown_h_thresh[NTHRESHOLDS];
4349 double zeroes[NTHRESHOLDS];
4350 Tools::Zero(zeroes, NTHRESHOLDS);
4353 for (
int i=0;i<NTHRESHOLDS;i++) {
4354 rate_v_thresh[i]=npass_v_thresh[i]/denom_v_thresh[i];
4356 if (npass_v_thresh[i]<=20) {
4357 errorup_v_thresh[i]=poissonerror_plus[(int)npass_v_thresh[i]]/denom_v_thresh[i];
4358 errordown_v_thresh[i]=poissonerror_minus[(int)npass_v_thresh[i]]/denom_v_thresh[i];
4361 if (settings1->WHICH==9 || settings1->WHICH==10){
4362 rate_h_thresh[i]=npass_h_thresh[i]/denom_h_thresh[i];
4363 if (npass_h_thresh[i]<=20) {
4364 errorup_h_thresh[i]=poissonerror_plus[(int)npass_h_thresh[i]]/denom_h_thresh[i];
4365 errordown_h_thresh[i]=poissonerror_minus[(int)npass_h_thresh[i]]/denom_h_thresh[i];
4371 string stemp=string(outputdir.Data())+
"/thresholds.root";
4372 TFile *fthresholds=
new TFile(stemp.c_str(),
"RECREATE");
4373 TCanvas *cthresh=
new TCanvas(
"cthresh",
"cthresh", 880, 800);
4376 TGraph *gnpass=
new TGraph(NTHRESHOLDS, thresholds, npass_v_thresh);
4377 gnpass->SetName(
"npass");
4378 TGraph *gdenom=
new TGraph(NTHRESHOLDS, thresholds, denom_v_thresh);
4379 gdenom->SetName(
"denom");
4381 TGraphAsymmErrors *g=
new TGraphAsymmErrors(NTHRESHOLDS, thresholds, rate_v_thresh, zeroes, zeroes, errorup_v_thresh, errordown_v_thresh);
4385 g->SetMarkerStyle(21);
4388 stemp = string(outputdir.Data())+
"/thresholds.eps";
4389 cthresh->Print((TString)stemp);
4395 if (settings1->WHICH==9 || settings1->WHICH==10){
4396 TGraph *gnpassH=
new TGraph(NTHRESHOLDS, thresholds, npass_h_thresh);
4397 gnpassH->SetName(
"npassH");
4398 TGraph *gdenomH=
new TGraph(NTHRESHOLDS, thresholds, denom_h_thresh);
4399 gdenomH->SetName(
"denomH");
4400 TGraphAsymmErrors *gH=
new TGraphAsymmErrors(NTHRESHOLDS, thresholds, rate_h_thresh, zeroes, zeroes, errorup_h_thresh, errordown_h_thresh);
4401 gH->SetName(
"rate");
4402 gH->SetLineWidth(2);
4403 gH->SetMarkerStyle(21);
4405 stemp = string(outputdir.Data())+
"/thresholds_HPOL.eps";
4406 cthresh->Print((TString)stemp);
4413 fthresholds->Write();
4414 fthresholds->Close();
4422 foutput <<
"Generated " << NNU <<
" neutrinos with energy " << pnu <<
"\n";
4423 foutput <<
"Number of unweighted direct, reflected that pass is: " << count1->npass[0] <<
"\t" << count1->npass[1] <<
"\n";
4424 foutput <<
"Number of (weighted) neutrinos that pass is: " << eventsfound <<
"\n";
4425 foutput <<
"Number of (weighted) neutrinos that pass, multiplied by prob. of interacting in the ice, is: " << eventsfound_prob <<
"\n";
4426 foutput <<
"Number of weighted direct, reflected that pass is: " << allcuts_weighted[0] <<
"\t" << allcuts_weighted[1] <<
"\n";
4427 foutput <<
"Number of (weighted) neutrinos that pass (with weight>0.001) is: " << eventsfound_weightgt01 <<
"\n";
4428 foutput <<
"Number of (weighted) neutrinos that only traverse the crust is " << eventsfound_crust <<
" -> " << eventsfound_crust/eventsfound*100 <<
"%\n\n";
4429 foutput <<
"Number of (weighted) neutrinos that pass only VPOL trigger is: " << allcuts_weighted_polarization[0] <<
"\n";
4430 foutput <<
"Number of (weighted) neutrinos that pass only HPOL trigger is: " << allcuts_weighted_polarization[1] <<
"\n";
4431 foutput <<
"Number of (weighted) neutrinos that pass both pol triggers is: " << allcuts_weighted_polarization[2] <<
"\n\n";
4433 cout <<
"Number of (weighted) neutrinos that pass (with weight>0.001) is: " << eventsfound_weightgt01 <<
"\n";
4434 cout <<
"Number of (weighted) neutrinos that pass, multiplied by prob. of interacting in the ice, is: " << eventsfound_prob <<
"\n";
4435 cout <<
"Number of (weighted) neutrinos that only traverse the crust is " << eventsfound_crust <<
" -> " << eventsfound_crust/eventsfound*100 <<
"%\n\n";
4436 cout <<
"Number of (weighted) neutrinos that pass only VPOL trigger is: " << allcuts_weighted_polarization[0] <<
"\n";
4437 cout <<
"Number of (weighted) neutrinos that pass only HPOL trigger is: " << allcuts_weighted_polarization[1] <<
"\n";
4438 cout <<
"Number of (weighted) neutrinos that pass both pol triggers is: " << allcuts_weighted_polarization[2] <<
"\n\n";
4440 foutput <<
"Total number of thrown electron neutrinos is : " << (double)count1->nnu_e <<
"\n";
4441 foutput <<
"Total number of thrown muon neutrinos is : " << (
double)count1->nnu_mu <<
"\n";
4442 foutput <<
"Total number of thrown tau neutrinos is : " << (double)count1->nnu_tau <<
"\n";
4443 foutput <<
"Number of (weighted) electron neutrinos that pass trigger is : " << sum[0] <<
"\n";
4444 foutput <<
"Number of (weighted) muon neutrinos that pass trigger is : " << sum[1] <<
"\n";
4445 foutput <<
"Number of (weighted) tau neutrinos that pass trigger is : " << sum[2] <<
"\n\n";
4447 cout <<
"Total number of thrown electron neutrinos is : " << (
double)count1->nnu_e <<
"\n";
4448 cout <<
"Total number of thrown muon neutrinos is : " << (double)count1->nnu_mu <<
"\n";
4449 cout <<
"Total number of thrown tau neutrinos is : " << (
double)count1->nnu_tau <<
"\n";
4450 cout <<
"Number of (weighted) electron neutrinos that pass trigger is : " << sum[0] <<
"\n";
4451 cout <<
"Number of (weighted) muon neutrinos that pass trigger is : " << sum[1] <<
"\n";
4452 cout <<
"Number of (weighted) tau neutrinos that pass trigger is : " << sum[2] <<
"\n\n";
4454 foutput <<
"Volume of ice is " << volume <<
"\n";
4455 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";
4471 double nevents_db=0;
4472 double nevents_nfb=0;
4475 nevents+=eventsfound;
4476 nevents_db+=eventsfound_db;
4477 nevents_nfb+=eventsfound_nfb;
4489 error_percent_increase_nfb=0;
4493 for (
int i=0;i<NBINS;i++) {
4494 if (eventsfound_binned[i]<=20) {
4495 error_plus+=pow(poissonerror_plus[(
int)eventsfound_binned[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4496 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);
4497 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);
4498 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);
4499 error_minus+=pow(poissonerror_minus[(
int)eventsfound_binned[i]]*pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4500 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);
4501 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);
4502 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);
4505 error_plus+=eventsfound_binned[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4506 error_e_plus+=eventsfound_binned_e[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4507 error_mu_plus+=eventsfound_binned_mu[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4508 error_tau_plus+=eventsfound_binned_tau[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*(MAX_LOGWEIGHT-MIN_LOGWEIGHT)+MIN_LOGWEIGHT), 2);
4509 error_minus=error_plus;
4510 error_e_minus=error_e_plus;
4511 error_mu_minus=error_mu_plus;
4512 error_tau_minus=error_tau_plus;
4515 error_nfb+=eventsfound_nfb_binned[i]*pow(pow(10., ((
double)i+0.5)/(
double)NBINS*4.+-5.), 2);
4516 for (
int j=0;j<NBINS_DISTANCE;j++) {
4517 if (eventsfound_binned_distance_forerror[j][i]<=20) {
4518 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);
4519 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);
4525 error_plus=sqrt(error_plus);
4526 error_e_plus=sqrt(error_e_plus);
4527 error_mu_plus=sqrt(error_mu_plus);
4528 error_tau_plus=sqrt(error_tau_plus);
4530 error_minus=sqrt(error_minus);
4531 error_e_minus=sqrt(error_e_minus);
4532 error_mu_minus=sqrt(error_mu_minus);
4533 error_tau_minus=sqrt(error_tau_minus);
4536 error_nfb=sqrt(error_nfb);
4537 for (
int j=0;j<NBINS_DISTANCE;j++) {
4538 error_distance_plus[j]=sqrt(error_distance_plus[j]);
4539 error_distance_minus[j]=sqrt(error_distance_minus[j]);
4544 if (NNU != 0 && nevents!=0) {
4545 km3sr=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*nevents/(double)NNU/settings1->SIGMA_FACTOR;
4547 cout << nevents <<
" events passed out of " << NNU <<
"\n";
4549 error_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_plus/(double)NNU/settings1->SIGMA_FACTOR;
4550 error_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_minus/(double)NNU/settings1->SIGMA_FACTOR;
4552 percent_increase_db=km3sr_db/km3sr*100;
4553 percent_increase_nfb=km3sr_nfb/km3sr*100;
4555 error_percent_increase_nfb=sqrt(pow(error_km3sr_nfb/km3sr_nfb, 2)+pow(error_plus/km3sr, 2))*percent_increase_nfb;
4557 percent_increase_total=percent_increase_db+percent_increase_nfb;
4559 km3sr_e = (pow(1.e-3, 3))*volume*sr*(sum[0]/(
double)count1->nnu_e)*sig1->RHOMEDIUM/sig1->RHOH20/settings1->SIGMA_FACTOR;
4561 cout << sum[0]/(
double)nevents*100. <<
"% are electron neutrinos\n";
4563 error_e_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_e_plus/(double)count1->nnu_e/settings1->SIGMA_FACTOR;
4564 error_e_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_e_minus/(double)count1->nnu_e/settings1->SIGMA_FACTOR;
4566 km3sr_mu = (pow(1.e-3, 3))*volume*sr*(sum[1]/(
double)count1->nnu_mu)*sig1->RHOMEDIUM/sig1->RHOH20/settings1->SIGMA_FACTOR;
4568 cout << sum[1]/(
double)nevents*100. <<
"% are muon neutrinos\n";
4570 error_mu_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_mu_plus/(double)count1->nnu_mu/settings1->SIGMA_FACTOR;
4571 error_mu_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_mu_minus/(double)count1->nnu_mu/settings1->SIGMA_FACTOR;
4573 km3sr_tau = (pow(1.e-3, 3))*volume*sr*(sum[2]/(
double)count1->nnu_tau)*sig1->RHOMEDIUM/sig1->RHOH20/settings1->SIGMA_FACTOR;
4575 cout << sum[2]/(
double)nevents*100. <<
"% are tau neutrinos\n";
4577 error_tau_plus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_tau_plus/(double)count1->nnu_tau/settings1->SIGMA_FACTOR;
4578 error_tau_minus=volume*pow(1.E-3, 3)*sig1->RHOMEDIUM/sig1->RHOH20*sr*error_tau_minus/(double)count1->nnu_tau/settings1->SIGMA_FACTOR;
4581 double sum_km3sr_error_plus=0;
4582 double sum_km3sr_error_minus=0;
4583 for (
int i=0;i<NBINS_DISTANCE;i++) {
4584 sum_km3sr+= (pow(1.e-3, 3))*volume*sr*eventsfound_binned_distance[i]*sig1->RHOMEDIUM/sig1->RHOH20/(
double)NNU/settings1->SIGMA_FACTOR;
4585 km3sr_distance[i]=sum_km3sr;
4586 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;
4587 error_distance_plus[i]=sqrt(sum_km3sr_error_plus);
4588 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;
4589 error_distance_minus[i]=sqrt(sum_km3sr_error_minus);
4590 distanceout << 700000./(double)NBINS_DISTANCE*(
double)i <<
"\t" << km3sr_distance[i] <<
"\t" << error_distance_plus[i] <<
"\t" << error_distance_minus[i] <<
"\n";
4593 foutput <<
"Total volume * solid angle is \t\t\t\t" << km3sr <<
" + " << error_plus <<
" - " << error_minus <<
" km^3 str\n";
4594 foutput <<
"Total volume * solid angle for electron neutrinos is \t" << km3sr_e <<
" + " << error_e_plus <<
" - " << error_e_minus <<
" km^3 str\n";
4595 foutput <<
"Total volume * solid angle for muon neutrinos is \t" << km3sr_mu <<
" + " << error_mu_plus <<
" - " << error_mu_minus <<
" km^3 str\n";
4596 foutput <<
"Total volume * solid angle for tau neutrinos is \t" << km3sr_tau <<
" + " << error_tau_plus <<
" - " << error_tau_minus <<
" km^3 str\n";
4597 cout <<
"Total volume * solid angle is \t\t\t\t" << km3sr <<
" + " << error_plus <<
" - " << error_minus <<
" km^3 str\n";
4598 cout <<
"Total volume * solid angle for electron neutrinos is \t" << km3sr_e <<
" + " << error_e_plus <<
" - " << error_e_minus <<
" km^3 str\n";
4599 cout <<
"Total volume * solid angle for muon neutrinos is \t" << km3sr_mu <<
" + " << error_mu_plus <<
" - " << error_mu_minus <<
" km^3 str\n";
4600 cout <<
"Total volume * solid angle for tau neutrinos is \t" << km3sr_tau <<
" + " << error_tau_plus <<
" - " << error_tau_minus <<
" km^3 str\n";
4602 if ( spectra1->IsMonoenergetic() ) {
4603 cout <<
"Cross section is " << sigma <<
"m^2\n";
4604 double len_int=1.0/(sigma*sig1->RHOH20*(1./M_NUCL)*1000);
4605 cout <<
"Interaction length is " << len_int <<
" km\n\n";
4606 foutput <<
"Interaction length is " << len_int <<
"\n";
4607 km2sr=km3sr/len_int;
4608 foutput <<
"Total area x steradians using km3sr/len_int is \t\t\t\t" << km2sr <<
" km^2 str\n\n";
4609 cout <<
"Total area x steradians is \t\t" << km2sr <<
" km^2 str\n\n";
4610 km2sr=ice_area/(1.E6)*PI*eventsfound_prob/(
double)NNU;
4611 foutput <<
"Total area x steradians using 4*PI*R_EARTH^2*eff. is \t" << km2sr <<
" km^2 str\n\n";
4612 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";
4619 foutput <<
"\t\t\t\t\t\t\tprobability of passing \t# events passing\n";
4621 foutput.precision(4);
4622 foutput <<
"No way this neutrino will see any ice \t\t\t" << (double)count1->noway[0]/(
double)count_total <<
"\t" <<
4623 (double)count1->noway[1]/(
double)count_total <<
"\t\t" <<
4624 count1->noway[0] <<
"\t" << count1->noway[1] <<
"\n";
4626 foutput.precision(4);
4627 foutput <<
"Wheredoesitleave in PickUnbiased gives an error \t" << (double)count1->wheredoesitleave_err[0]/(
double)count1->noway[0] <<
"\t" <<
4628 (double)count1->wheredoesitleave_err[1]/(
double)count1->noway[1] <<
"\t\t" <<
4629 count1->wheredoesitleave_err[0] <<
"\t" << count1->wheredoesitleave_err[1] <<
"\n";
4631 foutput.precision(4);
4632 foutput <<
"This neutrino direction never sees ice \t\t\t" << (double)count1->neverseesice[0]/(
double)count1->wheredoesitleave_err[0] <<
"\t" <<
4633 (double)count1->neverseesice[1]/(
double)count1->wheredoesitleave_err[1] <<
"\t\t" <<
4634 count1->neverseesice[0] <<
"\t" << count1->neverseesice[1] <<
"\n";
4637 foutput.precision(4);
4638 foutput <<
"WhereDoesItEnterIce in PickUnbiased gives an error \t\t\t" << (double)count1->wheredoesitenterice_err[0]/(
double)count1->neverseesice[0] <<
"\t" <<
4639 (double)count1->wheredoesitenterice_err[1]/(
double)count1->neverseesice[1] <<
"\t\t" <<
4640 count1->wheredoesitenterice_err[0] <<
"\t" << count1->wheredoesitenterice_err[1] <<
"\n";
4642 foutput.precision(4);
4643 foutput <<
"Interaction point too high \t\t\t\t" << (double)count1->toohigh[0]/(
double)count1->wheredoesitenterice_err[0] <<
"\t" <<
4644 (double)count1->toohigh[1]/(
double)count1->wheredoesitenterice_err[1] <<
"\t\t" <<
4645 count1->toohigh[0] <<
"\t" << count1->toohigh[1] <<
"\n";
4647 foutput.precision(4);
4648 foutput <<
"Interaction point too low \t\t\t\t" << (double)count1->toolow[0]/(
double)count1->toohigh[0] <<
"\t" <<
4649 (double)count1->toolow[1]/(
double)count1->toohigh[1] <<
"\t\t" <<
4650 count1->toolow[0] <<
"\t" << count1->toolow[1] <<
"\n";
4655 foutput.precision(4);
4656 foutput <<
"There is an interaction in ice \t\t\t\t" << (double)count1->iceinteraction[0]/(
double)count1->toolow[0] <<
"\t" <<
4657 (double)count1->iceinteraction[1]/(
double)count1->toolow[1] <<
"\t\t" <<
4658 count1->iceinteraction[0] <<
"\t" << count1->iceinteraction[1] <<
"\n";
4660 foutput.precision(4);
4661 foutput <<
"In horizon \t\t\t\t\t\t" << (double)count1->inhorizon[0]/(
double)count1->iceinteraction[0] <<
"\t" <<
4662 (double)count1->inhorizon[1]/(
double)count1->iceinteraction[1] <<
"\t\t" <<
4663 count1->inhorizon[0] <<
"\t" << count1->inhorizon[1] <<
"\n";
4667 foutput.precision(4);
4668 foutput <<
"From surface to balloon, ray not intersected by earth \t" << (double)count1->nraypointsup1[0]/(
double)count1->inhorizon[0] <<
"\t" <<
4669 (double)count1->nraypointsup1[1]/(
double)count1->inhorizon[1] <<
"\t\t" <<
4670 count1->nraypointsup1[0] <<
"\t" << count1->nraypointsup1[1] <<
"\n";
4672 foutput.precision(4);
4673 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" <<
4674 (double)count1->nnottoosmall[1]/(
double)count1->nraypointsup1[1] <<
"\t\t" 4675 << count1->nnottoosmall[0] <<
"\t" << count1->nnottoosmall[1] <<
"\n";
4678 foutput.precision(4);
4679 foutput <<
"Viewing angle lt 90 degrees\t\t\t" << (double)count1->nviewangle_lt_90[0]/(
double)count1->nnottoosmall[0] <<
"\t" <<
4680 (double)count1->nviewangle_lt_90[1]/(
double)count1->nnottoosmall[1] <<
"\t\t" 4681 << count1->nviewangle_lt_90[0] <<
"\t" << count1->nviewangle_lt_90[1] <<
"\n";
4684 foutput.precision(4);
4685 foutput <<
"Reality check: EM and hadronic fractions both nonzero\t" << (double)count1->ngoodfracs[0]/(
double)count1->nviewangle_lt_90[0] <<
"\t" <<
4686 (double)count1->ngoodfracs[1]/(
double)count1->nviewangle_lt_90[1] <<
"\t\t" 4687 << count1->ngoodfracs[0] <<
"\t" << count1->ngoodfracs[1] <<
"\n";
4688 foutput.precision(4);
4689 foutput <<
"\tBoth EM and hadronic fractions are zero\t\t" << (double)count1->nbadfracs[0]/(
double)count1->nviewangle_lt_90[0] <<
"\t" <<
4690 (double)count1->nbadfracs[1]/(
double)count1->nviewangle_lt_90[1] <<
"\t\t" <<
4691 count1->nbadfracs[0] <<
"\t" << count1->nbadfracs[1] <<
"\n";
4693 foutput.precision(4);
4694 foutput <<
"After finding neutrino direction, \n\tchance of making through Earth\t\t\t" << (double)count_chanceofsurviving/(
double)count1->ngoodfracs[0] <<
"\t\t\t";
4695 foutput.precision(10);
4696 foutput << count_chanceofsurviving <<
"\n";
4697 foutput.precision(4);
4698 foutput <<
"Neutrino enters ice south of 60deg S latitude\t\t" << (double)count1->nentersice[0]/(
double)count_chanceofsurviving <<
"\t" <<
4699 (double)count1->nentersice[1]/(
double)count_chanceofsurviving <<
4701 count1->nentersice[0] <<
"\t" << count1->nentersice[1] <<
"\n";
4703 foutput.precision(4);
4704 foutput <<
"Neutrino reasonably likely to survive trip through Earth " << (double)count1->nabsorbed[0]/(
double)count1->nentersice[0] <<
"\t" <<
4705 (double)count1->nabsorbed[1]/(
double)count1->nentersice[1] <<
"\t\t" 4706 << count1->nabsorbed[0] <<
"\t" << count1->nabsorbed[1] <<
"\n";
4707 foutput.precision(4);
4708 foutput <<
"Ray leaves the ice south of 60deg S latitude\t\t" << (double)count1->nraywithincontinent1[0]/(
double)count1->nabsorbed[0] <<
"\t" <<
4709 (double)count1->nraywithincontinent1[1]/(
double)count1->nabsorbed[1] <<
"\t" <<
4710 count1->nraywithincontinent1[0] <<
"\t" <<
4711 count1->nraywithincontinent1[1] <<
"\n";
4713 foutput.precision(4);
4714 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";
4715 foutput.precision(10);
4716 foutput <<count_chanceinhell0 <<
"\n";
4718 foutput.precision(4);
4719 foutput <<
"Ray is not totally internally reflected\t\t\t" << (double)count1->nnottir[0]/(
double)count_chanceinhell0 <<
"\t" <<
4720 (double)count1->nnottir[1]/(
double)count_chanceinhell0 <<
"\t\t" <<
4721 count1->nnottir[0] <<
"\t" << count1->nnottir[1] <<
"\n";
4723 foutput.precision(4);
4724 foutput <<
"From surface to balloon, ray not intersected by earth \t" << (double)count1->nraypointsup2[0]/(
double)count1->nnottir[0] <<
"\t" <<
4725 (double)count1->nraypointsup2[1]/(
double)count1->nnottir[1] <<
4727 << count1->nraypointsup2[0] <<
"\t" << count1->nraypointsup2[1] <<
"\n";
4728 foutput.precision(4);
4730 foutput <<
"Ray leaves the ice south of 60deg S latitude\t\t" << (double)count1->nraywithincontinent2[0]/(
double)count1->nraypointsup2[0] <<
"\t" <<
4731 (double)count1->nraywithincontinent2[0]/(
double)count1->nraypointsup2[1] <<
4732 "\t\t" << count1->nraywithincontinent2[0] <<
"\t" <<
4733 count1->nraywithincontinent2[1] <<
"\n";
4734 foutput.precision(4);
4736 foutput <<
"Ray leaves where there is ice\t\t\t\t" << (double)count1->nacceptablerf[0]/(
double)count1->nraywithincontinent2[0] <<
"\t" <<
4737 (double)count1->nacceptablerf[1]/(
double)count1->nraywithincontinent2[1] <<
"\t\t" 4738 << count1->nacceptablerf[0] <<
"\t" << count1->nacceptablerf[1] <<
"\n";
4739 foutput.precision(4);
4741 foutput <<
"Ray tracing converges to within 10 m\t\t\t" << (double)count1->nconverges[0]/(
double)count1->nacceptablerf[0] <<
"\t" <<
4742 (double)count1->nconverges[1]/(
double)count1->nacceptablerf[1] <<
4743 "\t\t" << count1->nconverges[0] <<
"\t" << count1->nconverges[1] <<
"\n";
4744 foutput.precision(4);
4746 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] <<
4747 "\t\t" <<count1->nchanceinhell_fresnel[0] <<
"\t" << count1->nchanceinhell_fresnel[1] <<
"\n";
4748 foutput.precision(4);
4749 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";
4750 foutput.precision(4);
4752 foutput <<
"After ice attenuation, \n\tmaximum signal is detectable\t\t\t" << (double)count1->nchanceinhell[0]/(
double)count1->nchanceinhell_1overr[0] <<
"\t" <<
4753 (double)count1->nchanceinhell[1]/(
double)count1->nchanceinhell_1overr[1] <<
"\t" <<count1->nchanceinhell[0] <<
"\t" << count1->nchanceinhell[1] <<
"\n";
4754 foutput.precision(4);
4756 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";
4758 foutput.precision(4);
4759 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";
4761 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";
4763 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";
4764 foutput <<
"Number of l1 triggers\t\t\t\t\t\t" << (
double)count1->nl1triggers[0][0] <<
"\t" << (double)count1->nl1triggers[1][0] <<
"\n";
4766 foutput <<
"Chord is good length\t\t\t\t\t" << (
double)count_chordgoodlength/(double)count1->npassestrigger[0] <<
"\t\t\t";
4767 foutput.precision(10);
4768 foutput <<count_chordgoodlength <<
"\n";
4769 foutput.precision(4);
4770 foutput <<
"Neutrino's path in ice more than 1m \t\t\t" << (double)count_d2goodlength/(
double)count_chordgoodlength <<
"\t\t\t";
4771 foutput.precision(10);
4772 foutput << count_d2goodlength <<
"\n";
4773 foutput.precision(4);
4774 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";
4775 foutput.precision(10);
4776 foutput <<count1->npass[0] <<
"\t" << count1->npass[1] <<
"\n";
4778 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";
4779 cout <<count1->npass[0] <<
"\t" << count1->npass[1] <<
"\n";
4785 TString str; str.Form(
"%s/source_times.txt", outputdir.Data());
4786 FILE * f = fopen(str.Data(),
"w");
4787 std::vector<double> times;
4789 for (
unsigned i = 0; i < times.size(); i++)
4791 fprintf(f,
"%f\n", times[i]);
4797 if ( src_model || spectra1->IsSpectrum() ) {
4798 double sum_events=0.;
4799 double thisenergy=0.;
4800 double thislen_int_kgm2=0.;
4808 double min_energy = src_model ? settings1->SOURCE_MIN_E : spectra1->Getenergy()[0];
4809 double max_energy = src_model ? settings1->SOURCE_MAX_E : spectra1->Getenergy()[spectra1->GetE_bin()-1];
4810 even_E = ( max_energy - min_energy ) / ( (
double) N_even_E );
4815 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);
4821 for (
int i=0;i<N_even_E;i++) {
4822 thisenergy=pow(10., (min_energy+((
double)i)*even_E));
4823 primary1->
GetSigma(thisenergy, sigma, thislen_int_kgm2, settings1, xsecParam_nutype, xsecParam_nuint);
4826 double EdNdEdAdt = src_model ? 1e9*src_flux->Interpolate(thisenergy*1e-9) : spectra1->GetEdNdEdAdt(log10(thisenergy));
4834 sum_events+=even_E*log(10.)*( EdNdEdAdt*1e4 )/(thislen_int_kgm2/sig1->RHOH20);
4835 integral+=even_E*log(10.)*(EdNdEdAdt);
4836 cout <<
"thisenergy, EdNdEdAdt is " << thisenergy <<
" " << EdNdEdAdt <<
"\n";
4852 cout <<
"SUM EVENTS IS " << sum_events << endl;
4853 cout <<
"INTEGRAL : " << integral << endl;
4854 sum_events*=volume*anita1->LIVETIME*sig1->RHOMEDIUM/sig1->RHOH20*nevents/(double)NNU*sr;
4856 foutput <<
"volume, LIVETIME, sig1->RHOMEDIUM, RHOH20, nevents, NNU, sr are " << volume <<
" " << anita1->LIVETIME <<
" " << sig1->RHOMEDIUM <<
" " << sig1->RHOH20 <<
" " << nevents <<
" " << NNU <<
" " << sr <<
"\n";
4857 foutput <<
"Total events observed is " << sum_events <<
"\n";
4864 double GetAirDistance(
double altitude_bn,
double beta) {
4865 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));
4870 double GetAverageVoltageFromAntennasHit(
Settings *settings1,
int *nchannels_perrx_triggered,
double *voltagearray,
double& volts_rx_sum) {
4872 int count_hitantennas=0;
4873 for (
int i=0;i<settings1->NANTENNAS;i++) {
4874 if (nchannels_perrx_triggered[i]>=3) {
4875 sum+=voltagearray[i];
4876 count_hitantennas++;
4880 sum = sum/(double)count_hitantennas;
4891 Vector n_bfield = nnu.Cross(nrf2_iceside);
4893 Vector n_pol = n_bfield.Cross(nrf2_iceside);
4894 n_pol = n_pol.Unit();
4896 if (nnu*nrf2_iceside>0 && n_pol*nnu>0){
4897 cout <<
"error in GetPolarization. Event is " << inu <<
"\n";
4904 void CloseTFile(TFile *hfile) {
4912 double IsItDoubleBang(
double exitlength,
double plepton) {
4913 double gamma=plepton/MTAU;
4914 return 1-exp(-1*exitlength/(TAUDECAY_TIME*CLIGHT*gamma));
4922 double gamma=pnu/MTAU;
4924 if (exp(-1*nuexitlength/(TAUDECAY_TIME*CLIGHT*gamma))>0.999){
4925 rnd1=getRNG(RNG_SECOND_BANG)->Rndm()*nuexitlength;
4928 while (rnd2>1-exp(-1*rnd1/(TAUDECAY_TIME*CLIGHT*gamma))) {
4929 rnd1=getRNG(RNG_SECOND_BANG)->Rndm()*nuexitlength;
4930 rnd2=getRNG(RNG_SECOND_BANG)->Rndm();
4933 posnu2 = posnu + rnd1*nnu;
4934 rfexit_db = antarctica1->Surface(posnu2)*posnu2.Unit();
4937 n_exit2bn_db = (r_bn - rfexit_db) / r_bn.
Distance(rfexit_db);
4939 double cosangle=(n_exit2bn_db * posnu2) / posnu2.Mag();
4950 double ATTENLENGTH=700;
4951 if(!settings1->VARIABLE_ATTEN){
4952 ATTENLENGTH=antarctica1->EffectiveAttenuationLength(settings1, posnu, 1);
4955 int position_in_iceshelves=antarctica1->IceOnWater(posnu);
4956 int position_in_rossexcept=antarctica1->RossExcept(posnu);
4959 double dtemp=posnu_down.
Distance(rfexit2)/ATTENLENGTH;
4963 if(position_in_iceshelves && (!position_in_rossexcept)){
4966 scalefactor_attenuation=0.71*exp(-dtemp);
4967 vmmhz_max=vmmhz_max*0.71*exp(-dtemp);
4969 else if(position_in_rossexcept){
4970 scalefactor_attenuation=0.1*exp(-dtemp);
4971 vmmhz_max=0.1*vmmhz_max*exp(-dtemp);
4974 scalefactor_attenuation=sqrt(0.001)*exp(-dtemp);
4975 vmmhz_max=sqrt(0.001)*vmmhz_max*exp(-dtemp);
4979 scalefactor_attenuation=0;
4986 void Attenuate(
IceModel *antarctica1,
Settings *settings1,
double& vmmhz_max,
double rflength,
const Position &posnu) {
4987 double ATTENLENGTH=700;
4988 if (!settings1->VARIABLE_ATTEN){
4989 ATTENLENGTH = antarctica1->EffectiveAttenuationLength(settings1, posnu, 0);
4992 double dtemp=(rflength/ATTENLENGTH);
4993 if(!settings1->ROUGHNESS){
4995 scalefactor_attenuation=exp(-dtemp);
4996 vmmhz_max=vmmhz_max*exp(-dtemp);
4999 scalefactor_attenuation=0;
5005 scalefactor_attenuation=exp(-dtemp);
5006 vmmhz_max=vmmhz_max*exp(-dtemp);
5009 scalefactor_attenuation=0;
5017 void IsAbsorbed(
double chord_kgm2,
double len_int_kgm2,
double &weight1) {
5023 rtemp=chord_kgm2/len_int_kgm2;
5025 weight1=exp(-rtemp);
5032 void GetSmearedIncidentAngle(
Vector &specular,
Vector &nrf_iceside,
Vector &n_exit2bn,
double SMEARINCIDENTANGLE){
5035 specular+=nrf_iceside;
5036 Vector parallel_to_surface;
5037 parallel_to_surface+=n_exit2bn;
5038 parallel_to_surface.Cross(specular);
5039 nrf_iceside.Rotate(SMEARINCIDENTANGLE*(2*getRNG(RNG_SMEARED_INCIDENT_ANGLE)->Rndm()-1.), parallel_to_surface);
5045 int GetRayIceSide(
const Vector &n_exit2rx,
const Vector &nsurf_rfexit,
double nexit,
double nenter,
Vector &nrf2_iceside) {
5048 double NRATIO=nexit/nenter;
5049 costh=(n_exit2rx*nsurf_rfexit)/(n_exit2rx.Mag() * nsurf_rfexit.Mag());
5055 double sinth=sqrt(1 - costh*costh);
5056 double factor=NRATIO*costh-sqrt(1-(NRATIO*sinth*NRATIO*sinth));
5057 nrf2_iceside = -factor*nsurf_rfexit + NRATIO*n_exit2rx;
5058 nrf2_iceside = nrf2_iceside.Unit();
5064 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) {
5072 double theta_test=0;
5073 double vmmhz1m_test=0;
5074 double costhetanu1 = 0;
5075 double costhetanu2 = 0;
5079 theta_threshold = 0;
5084 if (settings1->SKIPCUTS || !settings1->USEDIRECTIONWEIGHTS) {
5090 if (emfrac<=1.E-10 && deltheta_had >1.E-10) {
5091 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(hadfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle)>1)
5096 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);
5097 averaging_thetas1+=theta_threshold;
5099 count_inthisloop1++;
5102 if (emfrac>1.E-10 && deltheta_had <=1.E-10) {
5104 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(emfrac*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)*sin(sig1->changle)>1)
5109 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);
5110 averaging_thetas2+=theta_threshold;
5112 count_inthisloop2++;
5117 if (emfrac>1.E-10 && deltheta_had>1.E-10) {
5120 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/((hadfrac+emfrac)*vmmhz1m_max/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.) {
5122 theta_threshold=-1.;
5125 theta_test=deltheta_em;
5126 vmmhz1m_test=vmmhz1m_max;
5128 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5130 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.) {
5131 theta_threshold=theta_test;
5134 theta_test=1.5*deltheta_em;
5135 vmmhz1m_test=vmmhz1m_max;
5136 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5138 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.) {
5140 theta_threshold=theta_test;
5143 theta_test=2*deltheta_em;
5144 vmmhz1m_test=vmmhz1m_max;
5145 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5147 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.)
5148 theta_threshold=theta_test;
5150 theta_test=3*deltheta_em;
5151 vmmhz1m_test=vmmhz1m_max;
5152 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5155 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.)
5156 theta_threshold=theta_test;
5158 theta_test=deltheta_had;
5159 vmmhz1m_test=vmmhz1m_max;
5160 sig1->TaperVmMHz(sig1->changle+theta_test, deltheta_em, deltheta_had, emfrac, hadfrac, vmmhz1m_test, djunk);
5162 if (anita1->VNOISE[0]/10.*anita1->maxthreshold/(vmmhz1m_test/r_fromballoon*heff_max*anita1->bwmin/1.E6)>1.)
5164 theta_threshold=theta_test;
5166 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);
5176 count_inthisloop3++;
5177 averaging_thetas3+=theta_threshold;
5182 theta_threshold*=settings1->THETA_TH_FACTOR;
5183 if (theta_threshold>0) {
5184 if (sig1->changle-theta_threshold<0 && sig1->changle+theta_threshold> PI) {
5188 else if (sig1->changle-theta_threshold>0 && sig1->changle+theta_threshold> PI) {
5189 costhetanu2=cos(sig1->changle-theta_threshold);
5192 else if (sig1->changle-theta_threshold<0 && sig1->changle+theta_threshold< PI) {
5194 costhetanu1=cos(sig1->changle+theta_threshold);
5196 else if (sig1->changle-theta_threshold>0 && sig1->changle+theta_threshold< PI) {
5197 costhetanu2=cos(sig1->changle-theta_threshold);
5198 costhetanu1=cos(sig1->changle+theta_threshold);
5205 if (theta_threshold>0) {
5207 costhetanu=costhetanu1+getRNG(RNG_DIRECTION)->Rndm()*(costhetanu2-costhetanu1);
5209 double phinu=TWOPI*getRNG(RNG_DIRECTION)->Rndm();
5210 double sinthetanu=sqrt(1-costhetanu*costhetanu);
5212 nnu =
Vector(sinthetanu*cos(phinu), sinthetanu*sin(phinu), costhetanu);
5213 nnu = nnu.ChangeCoord(refr);
5221 double angle=(PI/2.-sig1->changle)-ray1->rfexit[0].Angle(ray1->nrf_iceside[4]);
5222 double thetaposnu=posnu.Theta();
5224 costhetanu=cos(PI/2+(thetaposnu-angle));
5225 sinthetanu=sqrt(1-costhetanu*costhetanu);
5229 nnu =
Vector(-1.*sinthetanu*cos(phinu), -1.*sinthetanu*sin(phinu), -1.*costhetanu);
5233 else if (theta_threshold==-1.) {
5234 cout <<
"theta_threshold is " << theta_threshold <<
"\n";
5237 else if (emfrac<=1.E-10 && deltheta_had <= 1.E-10) {
5238 cout <<
"Error: emfrac, hadfrac are (1st place)" << emfrac <<
" " << hadfrac <<
" " <<
"\n";
5277 double dtemp1 = r_bn.
Distance(rfexit);
5279 double dtemp2 = rfexit.
Distance(posnu1);
5283 double dtemp = dtemp1 + dtemp2;
5285 vmmhz1m_max= vmmhz1m_max/dtemp;
5287 scalefactor_distance=1/dtemp;
5298 void SetupViewangles(
Signal *sig1) {
5299 double viewangle_max=90.*RADDEG;
5300 double viewangle_min=30.*RADDEG;
5301 for (
int i=0;i<NVIEWANGLE-2;i++) {
5302 viewangles[i]=viewangle_max-(viewangle_max-viewangle_min)/(
double)(NVIEWANGLE-2)*(
double)i;
5304 viewangles[NVIEWANGLE-2]=acos(1/sig1->N_DEPTH);
5305 viewangles[NVIEWANGLE-1]=90.*RADDEG;
5310 double GetThisAirColumn(
Settings* settings1,
Position r_in,
Vector nnu,
Position posnu,
double *col1,
double& cosalpha,
double& mytheta,
double& cosbeta0,
double& mybeta) {
5313 cosalpha=(r_in * nnu) / r_in.Mag();
5314 mytheta=(double)(acos(cosalpha)*DEGRAD)-90.;
5317 if (settings1->ATMOSPHERE) {
5318 int index11=int(mytheta*10.);
5319 int index12=index11+1;
5321 myair=(col1[index11]+(col1[index12]-col1[index11])*(mytheta*10.-
double(index11)))*10.;
5328 cosbeta0= (posnu * nnu) / posnu.Mag();
5329 mybeta=(double)(acos(cosbeta0)*DEGRAD)-90.;
5335 void GetAir(
double *col1) {
5337 ifstream air1(ICEMC_SRC_DIR+
"/data/atmosphere.dat");
5340 for(
int iii=0;iii<900;iii++) {
5341 air1>>nothing>>col1[iii];
5347 int TIR(
const Vector &n_surf,
const Vector &nrf2_iceside,
double N_IN,
double N_OUT) {
5348 double test=sin(nrf2_iceside.Angle(n_surf))*N_IN/N_OUT;
5359 double GetViewAngle(
const Vector &nrf2_iceside,
const Vector &nnu) {
5362 double dtemp = nrf2_iceside*nnu;
5363 if (dtemp>=1 && dtemp<1.02)
5365 if (dtemp<=-1 && dtemp>-1.02)
5375 TStyle* RootStyle() {
5376 TStyle *RootStyle =
new TStyle(
"Root-Style",
"The Perfect Style for Plots ;-)");
5386 RootStyle->SetPaperSize(TStyle::kUSLetter);
5389 RootStyle->SetCanvasColor (0);
5390 RootStyle->SetCanvasBorderSize(10);
5391 RootStyle->SetCanvasBorderMode(0);
5392 RootStyle->SetCanvasDefH (600);
5393 RootStyle->SetCanvasDefW (600);
5394 RootStyle->SetCanvasDefX (10);
5395 RootStyle->SetCanvasDefY (10);
5398 RootStyle->SetPadColor (0);
5399 RootStyle->SetPadBorderSize (10);
5400 RootStyle->SetPadBorderMode (0);
5402 RootStyle->SetPadBottomMargin(0.16);
5403 RootStyle->SetPadTopMargin (0.08);
5404 RootStyle->SetPadLeftMargin (0.18);
5405 RootStyle->SetPadRightMargin (0.05);
5406 RootStyle->SetPadGridX (0);
5407 RootStyle->SetPadGridY (0);
5408 RootStyle->SetPadTickX (1);
5409 RootStyle->SetPadTickY (1);
5412 RootStyle->SetFrameFillStyle ( 0);
5413 RootStyle->SetFrameFillColor ( 0);
5414 RootStyle->SetFrameLineColor ( 1);
5415 RootStyle->SetFrameLineStyle ( 0);
5416 RootStyle->SetFrameLineWidth ( 2);
5417 RootStyle->SetFrameBorderSize(10);
5418 RootStyle->SetFrameBorderMode( 0);
5422 RootStyle->SetHistFillColor(0);
5423 RootStyle->SetHistFillStyle(1);
5424 RootStyle->SetHistLineColor(1);
5425 RootStyle->SetHistLineStyle(0);
5426 RootStyle->SetHistLineWidth(2);
5429 RootStyle->SetFuncColor(1);
5430 RootStyle->SetFuncStyle(0);
5431 RootStyle->SetFuncWidth(1);
5434 RootStyle->SetStatBorderSize(2);
5435 RootStyle->SetStatFont (42);
5437 RootStyle->SetOptStat (0);
5438 RootStyle->SetStatColor (0);
5439 RootStyle->SetStatX (0.93);
5440 RootStyle->SetStatY (0.90);
5441 RootStyle->SetStatFontSize (0.07);
5446 RootStyle->SetTickLength ( 0.015,
"X");
5447 RootStyle->SetTitleSize ( 0.10,
"X");
5448 RootStyle->SetTitleOffset( 1.20,
"X");
5449 RootStyle->SetTitleBorderSize(0);
5451 RootStyle->SetLabelOffset( 0.015,
"X");
5452 RootStyle->SetLabelSize ( 0.050,
"X");
5453 RootStyle->SetLabelFont ( 42 ,
"X");
5454 RootStyle->SetTickLength ( 0.015,
"Y");
5455 RootStyle->SetTitleSize ( 0.10,
"Y");
5456 RootStyle->SetTitleOffset( 0.600,
"Y");
5457 RootStyle->SetLabelOffset( 0.015,
"Y");
5458 RootStyle->SetLabelSize ( 0.050,
"Y");
5459 RootStyle->SetLabelFont ( 42 ,
"Y");
5460 RootStyle->SetTitleFont (42,
"XY");
5461 RootStyle->SetTitleColor (1);
5464 RootStyle->SetOptFit (0);
5465 RootStyle->SetMarkerStyle(20);
5466 RootStyle->SetMarkerSize(0.4);
5475 void GetFresnel(
Roughness *rough1,
int ROUGHNESS_SETTING,
const Vector &surface_normal,
5480 double emfrac,
double hadfrac,
double deltheta_em_max,
double deltheta_had_max,
5481 double &t_coeff_pokey,
double &t_coeff_slappy,
5486 double incident_angle = surface_normal.Angle(firn_rf);
5487 double transmitted_angle = surface_normal.Angle(air_rf);
5492 Vector perp = air_rf.Cross(surface_normal).Unit();
5494 Vector air_parallel = perp.Cross(air_rf).Unit();
5496 Vector firn_parallel = perp.Cross(firn_rf).Unit();
5499 double pol_perp_firn = pol*perp;
5500 double pol_parallel_firn = pol*firn_parallel;
5501 double pol_perp_air=0, pol_parallel_air=0;
5503 double r_coeff_pokey = tan(incident_angle - transmitted_angle) / tan(incident_angle + transmitted_angle);
5504 t_coeff_pokey = sqrt((1. - r_coeff_pokey*r_coeff_pokey));
5505 pol_parallel_air = t_coeff_pokey * pol_parallel_firn;
5507 double r_coeff_slappy = sin(incident_angle - transmitted_angle) / sin(incident_angle + transmitted_angle);
5508 t_coeff_slappy = sqrt((1. - r_coeff_slappy*r_coeff_slappy));
5509 pol_perp_air = t_coeff_slappy * pol_perp_firn;
5511 mag=sqrt( tan(incident_angle) / tan(transmitted_angle) );
5513 fresnel = sqrt( pow(efield * pol_perp_air, 2) + pow(efield * pol_parallel_air, 2)) / efield;
5515 pol = (pol_perp_air * perp + pol_parallel_air * air_parallel).Unit();
5539 const Vector nuvector = interaction1->
nnu;
5542 Vector zcoordvector = ray1->nsurf_rfexit;
5543 zcoordvector=zcoordvector.Unit();
5548 Vector xcoordvector = nuvector-(zcoordvector.Dot(nuvector))*zcoordvector;
5549 xcoordvector = xcoordvector.Unit();
5551 const Vector ycoordvector = zcoordvector.Cross(xcoordvector);
5555 if (interaction1->
nnu.Dot(zcoordvector)>0)
5556 origin_brian_tmp=interaction1->
nuexit;
5559 nnu_flipped=nnu_flipped-2.*nnu_flipped.Dot(zcoordvector)*zcoordvector;
5562 if (Ray::WhereDoesItLeave(interaction1->posnu,nnu_flipped,antarctica,nuexit_flipped))
5563 origin_brian_tmp=nuexit_flipped;
5567 r_bn_tmp=r_bn_tmp.ChangeCoord(xcoordvector,ycoordvector);
5570 double balloonphi = r_bn_tmp.Phi();
5572 balloonphi=balloonphi-2*PI;
5574 double balloontheta = r_bn_tmp.Theta();
5577 balloontheta = PI-balloontheta;
5588 void interrupt_signal_handler(
int sig){
5589 signal(sig, SIG_IGN);
5596 #ifdef ANITA_UTIL_EXISTS 5598 int GetIceMCAntfromUsefulEventAnt(
Settings *settings1,
int UsefulEventAnt){
5602 int IceMCAnt = UsefulEventAnt;
5603 if ((settings1->WHICH==9 || settings1->WHICH==10) && UsefulEventAnt<16) {
5604 IceMCAnt = (UsefulEventAnt%2==0)*UsefulEventAnt/2 + (UsefulEventAnt%2==1)*(UsefulEventAnt/2+8);
void AddTparallel_polPerpendicular(double A)
Appends a parallel transmission coefficient value to the fTcoeff_parl array.
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.
void SetUnitX(Vector a)
Sets an orientation vector of the screen.
Double_t rfExitPos[5][3]
Position where the RF exits the ice- 5 iterations, 3 dimensions each.
Vector GetUnitY()
Gets another orientation vector.
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.
void AddTperpendicular_polParallel(double A)
Appends a perpendicular transmission coefficient value to the fTcoeff_perp array. ...
void AddTparallel_polParallel(double A)
Appends a parallel transmission coefficient value to the fTcoeff_parl array.
void SetNvalidPoints(int i)
Sets the total number of points on the screen.
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)
void SetNormal(Vector a)
Sets the screen normal.
int currentint
Ditto - Stephen.
void AddViewangle(double A)
Appends a viewangle value to the fViewangle array.
void AddPol(Vector v)
Appends a vector to the fPols array.
Double_t showerE
Shower energy.
double GetViewangle(int i)
Get the viewangle value stored at the specified index.
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 SetCosineProjectionFactor(double a)
Sets the projection factor of the screen relative to the specular RF exit point.
double GetVmmhz0(int i)
Gets the Vmmhz0 value stored at the specified index.
void InitializeEachBand(Anita *anita1)
Initialize trigger bands.
void AddVec2bln(Vector v)
Appends a vector to the fVec2blns array.
Adu5Pat – The ADU5 Position and Attitude Data.
void TriggerPath(Settings *settings1, Anita *anita1, int ant, Balloon *bn1)
Apply trigger path.
double rms_rfcm_e_single_event
This is in Volts, not mV!
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.
void SetCentralPoint(Position a)
Sets the position of the central point of the screen.
Vector nnu_banana
Forced neutrino direction.
Position nu_banana_surface
The location of the surface above the forced neutrino interaction point.
double LAYER_VPOSITION[Anita::NLAYERS_MAX]
position of layers in z relative to vertical center of the payload
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.
int GetRxTriggerNumbering(int ilayer, int ifold)
get antenna number based on which layer and position it is
double THETA_ZENITH[NLAYERS_MAX]
how the antenna is tilted in theta (in radians with 0=up)
double chord
chord in m from earth entrance to rock-ice boundary
Int_t source_index
Name of the source.
double fVolts[12 *9][260]
Array of unwrapped (unless kNoCalib) voltages for each channel.
void SetUnitY(Vector a)
Sets another orientation vector of the screen.
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_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)
Double_t fSignalAtDigitizer[12 *9][260]
Array of signal at digitizer.
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.
void SetWeightNorm(double a)
Sets the normalization factor for the weights (so they sum to 1)
Double_t sourceTimeWeight
Relative Time weight for the given source assigned by icemc.
TTree * tdata
writing data out for the analysers
double GetVmmhz_freq(int i)
Get the Vmmhz value stored at the specified index.
int fNumPoints[12 *9]
Number of poins per channel.
UInt_t realTime
Time from the GPS unit.
const char * ICEMC_SRC_DIR()
void computeFluxTimeChanges(std::vector< double > *changes) const
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 GetNvalidPoints()
Gets the total number of points.
void Initialize(Settings *settings1, ofstream &foutput, int inu, TString outputdir)
initialize a bunch of stuff
double costheta_nutraject
theta of nnu with earth center to balloon as z axis
void AddImpactPt(Position p)
Appends a vector to the fImpactPt array.
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 SurfaceDistance(const Position &second, double local_surface) const
Returns "surface distance" between two positions.
void SetNsamples(int i)
Sets number of grid divisions for the base screen.
Class that handles the channel trigger.
void AddFacetLength(double A)
Appends a facet edge length value to the fFacetLength array.
Reads in and stores input settings for the run.
Double_t fNoiseAtDigitizer[12 *9][260]
Array of noise at digitizer.
double Lat() const
Returns latitude, where the +z direction is at 0 latitude.
double PHI_EACHLAYER[NLAYERS_MAX][NPHI_MAX]
phi of the center of each antenna on each layer
double PHI_OFFSET[NLAYERS_MAX]
antenna offset in phi for each layer (radians)
Functions you need to generate a primary interaction including cross sections and picking charged cur...
TString objName
Declination of source.
Double_t fNoiseAtTrigger[12 *9][260]
Array of noise at trigger.
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
void AddVmmhz0(double A)
Appends a Vmmhz value (for the lowest Anita frequency) to the fVmmhz0 array.
Double_t h_component[48]
H comp along polarization.
TTree * tgaryanderic
writing data out for the analysers
double rms_lab[2]
rms noise at lab chip
UInt_t realTime
unixTime of readout
Stores everything about a particular neutrino interaction. Interaction.
int fCapacitorNum[12 *9][260]
Array of capacitor numbers.
int GetNsamples()
Gets number of grid divisions for the base screen.
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 RRX[Anita::NLAYERS_MAX]
radius that the antenna sits from the axis of the payload (feedpoint)
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...
Float_t altitude
In metres.
Double_t hitangle_h[48]
Hit angles rel. to h plane stored for each antenna.
Position GetCentralPoint()
Gets the position of the screen's central point.
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
void DigitizerPath(Settings *settings1, Anita *anita1, int ant, Balloon *bn1)
Apply digitizer path.
Contains everything about positions within payload and signals it sees for each event, in both the trigger and signal paths.
int nuflavorint
Added by Stephen for output purposes.
void AddWeight(double a)
Appends a weight value to the fWeight array.
double GetEdgeLength()
Gets the screen length.
double rms_rfcm[2]
rms noise just after rfcm's
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)
double altitude_int
depth of interaction
void addSource(Source *source)
double FREQ_HIGH
highest frequency
void AddVmmhz_freq(double A)
Appends a Vmmhz value to the fVmmhz_freq array.
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 FREQ_LOW
lowest frequency
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.
double timedomain_output_allantennas[2][48][HALFNFOUR]
this is across all antennas, just the full band
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.
string nuflavor
neutrino flavor
Position GetPosition(int i, int j)
Calculates the physical position of the screen corresponding to the specified counter value...
Double_t deadTime
fractional deadTime
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
void SetEdgeLength(double a)
Sets the physical length of a side of the screen.
void AddTperpendicular_polPerpendicular(double A)
Appends a perpendicular transmission coefficient value to the fTcoeff_perp array. ...
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.
void ResetParameters()
Resets the following screen parameters (fNvalidpoints,fVmmhz_freq,fVmmhz0,fViewangle,fDelays,fVec2blns,fPols,fImpactPt,fWeight,fWeightNorm)
Double_t SNRAtTrigger[12 *9]
Array of SNR at trigger.
void AddIncidenceAngle(double A)
Appends an incidence angle value to the fIncAngles array.
Double_t SNRAtDigitizer[12 *9]
Array of SNR at digitizer.
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_t thresholds[12 *9]
Channel thresholds used in icemc.
double d2
ice-rock boundary to interaction point in m
void AddDelay(double A)
Appends a delay value to the fDelays array.
double Distance(const Position &second) const
Returns chord distance (direct distance between two vectors)
Double_t fSignalAtTrigger[12 *9][260]
Array of signal at trigger.
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 GetWeight(int i)
Gets the weight value at the specified index.
double NOTCH_MIN
low edge of notch filter. This is set in the input file
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 fTimes[12 *9][260]
Array of unwrapped (unless kNoCalib) times for each channel.
double d1
same as chord in m (earth entrance to rock-ice boundary)
Double_t fTimes[12 *9][260]
Array of unwrapped (unless kNoCalib) times for each channel.
int NRX_PHI[NLAYERS_MAX]
number of antennas around in each layer. (radians)
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 fDiodeOutput[12 *9][260]
Array of tunnel diode output.
Double_t balloonDir[3]
Balloon direction.
Ice thicknesses and water depth.
double Lon() const
Returns longitude.
void AddTransmissionAngle(double A)
Appends a transmission angle value to the fTransAngles array.
Vector nnu
direction of neutrino (+z in south pole direction)