Taumodel.cc

#include "vector.hh"
#include "Settings.h"
#include "vector.hh"
#include "position.hh"
#include "signal.hh"
#include "Primaries.h"
#include "secondaries.hh"
#include "icemodel.hh"
#include "Tools.h"
#include <sstream>
#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>

#include "TH1F.h"
#include "Constants.h"
#include "Settings.h"
#include "TTreeIndex.h"
#include "TChain.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TTree.h"
#include "TLegend.h"
#include "TLine.h"
#include "TROOT.h"
#include "TPostScript.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TText.h"
#include "TProfile.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TMath.h"
#include <unistd.h>
#include "TVector3.h"
#include "TRotation.h"
#include "TSpline.h"
#include "Taumodel.hh"
#include "icemc_random.h" 
using std::cout;
using std::stringstream;
using std::setprecision;
using std::accumulate;
using std::max_element;
using std::partial_sum;
using std::max;

Taumodel::Taumodel() {

  B0=1.2*pow(10.,-7); 
    B1=0.16*pow(10.,-7);
    E0=pow(10.,19);     
    mT=1.777E9;    
    cT=0.00008693;     
    Mn=1.672622E-24;   
    A=1.;              
        
    mydensityvector.clear();
    myavgdensityvector.clear();
    myenergyvector.clear();
    myPsurvvector.clear();
    etaufarray.clear();
    PDFarray.clear();
    
}//

//-------------------------------------------------------------

double Taumodel::GetTauWeight(Primaries *primary1, Settings *settings1,IceModel *antarctica1,Interaction *interaction1, double pnu, int nu_nubar, double& ptauf, int& crust_entered){ // 1 or 0 
                  // int& mantle_entered, // 1 or 0
                  // int& core_entered){//add secondaries?

  Vector chord3;
  Vector nchord;
  
  
  int currentint = interaction1->currentint;
  const Position earth_in = interaction1->r_in;
  double TauWeight = 0;
  const Position r_enterice = interaction1->r_enterice;
  const Position nuexitice = interaction1->nuexitice;
  int inu=4;
  //cout<<"inu is "<<inu<<"\n";


  chord3 = nuexitice - earth_in;
  double Distance=chord3.Mag();
  nchord = chord3 / Distance;
  
  
  double Etaui,Etau_final;
  double y=0;
  double yweight;
  
  double  zdistance=0; //m
  TauWeight=0;
  double TauWeight_tmp=0;
  double prob_at_z=0; 
  
  
  double Prob_Nu_Surv;
    
  double tau_surv;
  double sigma = 0;
  double len_int_kgm2 =0;
  
  primary1->GetSigma(pnu,sigma,len_int_kgm2,settings1,nu_nubar,currentint);
  
  double step=Tools::dMin(len_int_kgm2/densities[1]/10,25); 
  
  Position posnunow;
  double avgdensity=0;
  
  // double Etau_now;//=Etau_final;
  double Emin=1E15;
 
  int i=0;
 
  double taudensity=0;
  double dEtauidEtauf=0;
   
  int etauint =0;
  double startingz=0;
  
  double totaltaudistance=0;
  int totalnusteps=0;
  Vector nchord1;
  double enter_ice_mag = r_enterice.Distance(earth_in);
  
  mydensityvector.clear();
  myavgdensityvector.clear();
  this->GetDensityVectors(antarctica1,interaction1,nchord,step,Distance,totalnusteps, crust_entered);
  
  
  for(double logEtau_final=log10(Emin);logEtau_final<=log10(pnu);logEtau_final+=.01){//integral over energy (in log space?)
    Etau_final= pow(10,logEtau_final);
    i=0;
    int totalsteps=0;
    TauWeight_tmp=0;
    // Etau_now=Etau_final;
    double gamma = Etau_final/mT;
   
    //calculate the initial energy needed at the step so the tau will end at the correct final energy
    
    this->GetEnergyVector(Etau_final,step,totalnusteps,totalsteps,totaltaudistance, pnu);
    
    this->GetTauSurvVector(step,totalsteps);
    
    startingz = Distance-totaltaudistance;
    if(inu==7857){
       //cout<<"total nu steps is "<<totalnusteps<<"\n";
       //cout<<"density_Vector is "<<mydensityvector[mydensityvector.size()-1]<<"\n";
       //  cout<<"Distance is "<<Distance<<"\n";
       //cout<<"total nu steps is "<<totalnusteps<<"\n";
       //cout<<"vector size is "<<myenergyvector.size()<<"\n";
       //cout<<"totaltaudistance is "<<totaltaudistance<<"\n";
       //cout<<"crust, mantle, core are "<<crust_entered<<","<<mantle_entered<<","<<core_entered<<"\n";
       //cout<<"*********************************************************************************************** \n";
    }
    for(zdistance = startingz; zdistance<=Distance; zdistance +=step){
      int nustep = (int)(zdistance/step);
      int tauenergystep=totalsteps-i;
           
      nchord1 =  zdistance*nchord;
      posnunow = earth_in + nchord1; 
      
      avgdensity = myavgdensityvector[nustep];
      taudensity=mydensityvector[nustep];
      Etaui=myenergyvector[tauenergystep];
      tau_surv = myPsurvvector[i];
      y=1.-Etaui/pnu;
      
      if(zdistance >=enter_ice_mag){
        tau_surv=1;
      }
      
      if (Etaui>=pnu || Etaui<Etau_final || Etaui!=Etaui){//to catch anything that might make it through. Precaution
        prob_at_z =0;
      }
      else{
        yweight=primary1->Getyweight(pnu,y,nu_nubar,currentint);
        Prob_Nu_Surv = avgdensity/len_int_kgm2*exp(-zdistance*avgdensity*1./len_int_kgm2);
        dEtauidEtauf = exp(B1*taudensity*(Distance-zdistance))*Etaui/Etau_final;
        prob_at_z = Prob_Nu_Surv*yweight*1./pnu*dEtauidEtauf*step*tau_surv;
    /*  if(i==100 && Etaui > 1E18){
      cout<<"Prob_Nu_Surv is "<<Prob_Nu_Surv<<"\n";
      cout<<"avgdensity is "<<avgdensity<<"\n";
      cout<<"len_int_kgm2 is "<<len_int_kgm2<<"\n";
      cout<<"zdistance is "<<zdistance<<"\n";
      cout<<"yweight is "<<yweight<<"\n";
      cout<<"dE is "<<dEtauidEtauf<<"\n";
      cout<<"prob at z is "<<prob_at_z<<"\n";
      cout<<"TauWeight_tmp is "<<TauWeight_tmp<<"\n";
      
      }*/
        if(zdistance<=enter_ice_mag){
          prob_at_z*=1.-exp(-1.*(r_enterice.Distance(nuexitice)/(gamma*cT)));
        }
        else if(zdistance>enter_ice_mag){
          prob_at_z*=1.-exp(-1.*(posnunow.Distance(nuexitice)/(gamma*cT)));
        }
        TauWeight_tmp+=prob_at_z;
    //cout<<"prob at z is "<<prob_at_z<<"\n";

        i++;
      }//etaui<pnu
    }//zdist loop
     
     
    TauWeight_tmp*=log(10)*Etau_final;
    TauWeight_tmp*=.01; 
    
    TauWeight+=TauWeight_tmp;
    // cout<<"TauWeight_tmp is "<<TauWeight_tmp<<" at Etau_final is "<<Etau_final<<"\n";
    //prob_at_z_vector.push_back(row);
    etaufarray.push_back(Etau_final);
    PDFarray.push_back(TauWeight);
    etauint++;
    
  }//Etau_final loop
  
   // }//xloop
  double xrandom= TauWeight*getRNG(RNG_XRNDM)->Rndm();
 
  //cout<<"TauWeight is "<<TauWeight<<" xrandom is "<<xrandom<<"\n";
  for(size_t loopthrough =0; loopthrough<=PDFarray.size();loopthrough++){
    if(xrandom >PDFarray[loopthrough] && xrandom <PDFarray[loopthrough+1]){
      ptauf=etaufarray[loopthrough];
      break;
    }//if xrandom
    
  }//loopthrough
   //cout<<"ptauf is "<<ptauf<<"\n";
  weight_tau_prob = TauWeight;
  return 1;
} //GetTauWeight


void Taumodel::GetDensityVectors(IceModel *antarctica1,Interaction *interaction1, Vector nchord, double step, double Distance, int &totalnusteps,int &crust_entered){
   
    Vector nchord1;
    double avgdensity =0;//initilize average density.
    double density_total=0;//initilize running sum
    double density_now=0;//density at this step
    Position posnunow;
    Position postaunow;
    double altitude_tau;
    double lat_tau;
    const Position earth_in = interaction1->r_in;
    ofstream myNewFile;
    
   
  for(double taudistance=0;taudistance<=Distance;taudistance+=step){
     nchord1=taudistance*nchord;
     postaunow=earth_in+nchord1;
     lat_tau=postaunow.Lat();
     altitude_tau = postaunow.Mag()-antarctica1->Geoid(lat_tau);
     density_now=antarctica1->GetDensity(altitude_tau,postaunow,crust_entered);
     mydensityvector.push_back(density_now);
     
     density_total+=density_now;
     avgdensity=density_total/(totalnusteps+1);
     myavgdensityvector.push_back(avgdensity);
     totalnusteps++;
     
  }
  
}//Get Density Vectors

void Taumodel::GetEnergyVector(double Etau_final, double step,int totalnusteps, int &totalsteps, double &totaltaudistance, double pnu){
  
  myenergyvector.clear();
  myenergyvector.push_back(Etau_final);
  double Etau_now=Etau_final;
  double density_now;
  ofstream myNewFile_1;
 
  // double pnu;//SET UP EVENT CLASS
  for(int energysteps=totalnusteps-1;energysteps>0;energysteps--){
    density_now=mydensityvector[energysteps];
    
    Etau_now=Etau_now + (B0+B1*log(Etau_now/E0))*density_now*Etau_now*step;
        
    if(Etau_now <=pnu){
      totaltaudistance=(totalnusteps-energysteps)*step;
      totalsteps++;
      myenergyvector.push_back(Etau_now);
    }//Etau_now<=pnu
    else if(Etau_now >pnu){
      break;
    }
   
  } //energy steps


}//energy vector

void Taumodel::GetTauSurvVector(double step, int totalsteps){
  myPsurvvector.clear();
  double Etau_now;
  double tau_surv;
  for(int k1=totalsteps;k1>=0;k1--){
    tau_surv=1;
    for(int k2=k1-1;k2>=0;k2--){
      Etau_now=myenergyvector[k2];
      if (Etau_now >0)
        tau_surv=tau_surv*(1-step*mT/(cT*Etau_now));
      else
        tau_surv = 0;
    }
    myPsurvvector.push_back(tau_surv);
  }//tau surv

}//Tau surv vector