Primaries.h

//class Primaries:

#ifndef PRIMARIES_H_
#define PRIMARIES_H_

#include "TRandom3.h" 
#include <iostream>

#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TF3.h"
#include "TH2D.h"
#include "vector.hh"
#include "position.hh"

class Vector;
class Position;
class Interaction;
class Primaries;
class IceModel;
class Counting;
class Settings;

//using namespace std;
using std::string;
//
using std::cout;




class Y {

  

private:
    TF1* ffrac; 
    

    TF1* fC1_high[2][2]; 
    
    TF1 *fC1_low; 
    
    TF1 *fC2; 

    
    TF3* fy0_low;   
    TF2* fy0_high;     
    TRandom3 Rand3;

    double pickYConnollyetal2011(int NU,int CURRENT,double e); 
    // NU=0: nubar, NU=1: nu
    // CURRENT=0: CC, CURRENT-1: NC

    static constexpr double ymin_low = 0.00002;
    constexpr static double ymax_low=0.001;
    constexpr static double ymin_high=0.001;
    constexpr static double ymax_high=1.;
    constexpr static double dy_low=0.00002;
    constexpr static double dy_high=0.001;

    double pickYGandhietal(); 

    
    constexpr static double R1=0.36787944;  
    constexpr static double R2=0.63212056;  

    
public:
    Y();
    double pickY(Settings *settings1,double pnu,int nu_nubar,int currentint);
    double Getyweight(double pnu,double y,int nu_nubar,int currentint);
    
};//Y
class Primaries {
    
private:
    TRandom3 Rand3;
    
    TH2D *m_hsigma; 
    TCanvas *m_csigma;
    Y *m_myY; 
    int run_old_code;
    
public:
    double pickY(Settings *settings1,double pnu,int nu_nubar,int currentint);
    double Getyweight(double pnu,double y,int nu_nubar,int currentint);


    double A_low[4];
    double A0_high[2][2];
    double A1_high[2][2];
    double A2_high[2][2];
    double A3_high[2][2];
    double b0; 
    double b1; 
    
    TF1* m_fy[2][2];
    TF1* m_fsigma[2][2];
    
    double c0[2][2];
    double c1[2][2];
    double c2[2][2];
    double c3[2][2];
    double c4[2][2];
    
    static constexpr int NSIGMAS=2;
    double mine[NSIGMAS];
    double maxe[NSIGMAS]; 
    
    Primaries();
    ~Primaries();
    
    int GetSigma(double pnu,double& sigma,double &len_int_kgm2,Settings *settings1,int nu_nubar,int currentint);



    // string GetCurrent();
    string GetNuFlavor();
protected:
};

class Interaction  {
    
private:
    
    
    
    Vector tmp_banana; 
    
    
    static constexpr double nu_banana_theta_angle=-0.413 * 3.14159/180.;
    
    
    static constexpr double altitude_nu_banana=-400.;
    
    
    static constexpr double lat_nu_banana=0.; 
    static constexpr double lon_nu_banana=0.;
    
    
    static constexpr double banana_slopey=0.;
    static constexpr double nu_banana_phi_angle=0. * 3.14159/180.; 
    
    
    
public:
    
    static constexpr double phi_nu_banana=3.14159/4; 
    
    static constexpr double banana_observation_distance=600000.;
    static constexpr double theta_nu_banana=170.*3.14159/180.;

    static constexpr int kcc=0;
    static constexpr int knc=1;

    double banana_phi_obs;
    Vector banana_obs; 
    Interaction(string inttype,Primaries *primary1,Settings *settings1,int whichray,Counting *count1);

    void PickAnyDirection();
    
    int noway;
    int wheredoesitleave_err;
    int neverseesice;
    int wheredoesitenterice_err;
    int toohigh;
    int toolow;
    
    double pathlength_inice;
    
    Vector nnu;  
    double costheta_nutraject; 
    double phi_nutraject; 

    double weight_nu; 
    double weight_nu_prob; 
    
    Position r_in; 
    Position r_enterice; 
    Position nuexit; 
    Position nuexitice; 
    double chord;  
    double logchord; 
    double weight_bestcase; 
    double chord_kgm2_bestcase; 
    double chord_kgm2_ice; 
    double d1;  
    double d2;  
    
    
    static constexpr double pnu_banana=2.00E19;
    static constexpr double banana_y=0.2;
    double banana_weight;
    double banana_theta_obs;
    double banana_volts;
    static constexpr double banana_signal_fluct=0.;
    static constexpr double banana_sigma=0.;
    
    
    void  setNuFlavor(Primaries *primary1,Settings *settings1,int whichray,Counting *count1);
    string GetCurrent();
    void setCurrent();
    Position posnu;
    Position posnu_down;
    string  nuflavor;                   
    string  current;                    
    int nuflavorint;                
    int currentint;                 
    
    
    double surface_over_banana_nu;
    
    string banana_flavor; 
    string banana_current;  
    
    Vector nnu_banana; 
    
    Position nu_banana;  
    Position nu_banana_surface; 
    
    double dtryingdirection; 
    double dnutries; 

    double d_effective_area; 
    
    double altitude_int;
    double altitude_int_mirror;
    
    double r_fromballoon[2]; 
    
    double r_fromballoon_db; 
    double r_exit2bn; 
    double r_exit2bn_measured; 
    int iceinteraction;
    
    
protected:
};
#endif