LDFTest.cc

Go to the documentation of this file.

/*
     StationTrigger ErrorCode
     0  ok
     2  no-data, but station ok

     53 strange
     56 strange
*/



#define DEBUG 1






#include "LDFTest.h"

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/StationRecData.h>
#include <sevt/StationTriggerData.h>
#include <sevt/EventTrigger.h>

#include <utl/AxialVector.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Reader.h>
#include <utl/RandomEngine.h>
#include <utl/UTMPoint.h>

#include <det/Detector.h>
#include <sdet/SDetector.h>

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

#include <CLHEP/Random/RandPoisson.h>
#include <CLHEP/Random/Randomize.h>
#include <CLHEP/Random/RandGauss.h>

#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TStyle.h>
#include <TProfile.h>
#include <TCanvas.h>

#include <iostream>

using CLHEP::RandFlat;
using CLHEP::RandPoisson;



using namespace std;
using namespace LDFTestKG;
using namespace evt;
using namespace fwk;
using namespace utl;
using namespace atm;


const double LDFTest::fNerlingA1a =  6.42522;
const double LDFTest::fNerlingA1b = -1.53183;
const double LDFTest::fNerlingA2a =  168.168;
const double LDFTest::fNerlingA2b = -42.1368;

//
static double kMeterPerVEM = 1.2*meter;      // [m/VEM]
static double kTrackPerGEV = 5.25*meter/GeV; //(see E. Zas, Malargue Nov 04)
//static double MeanEnergy = 25.*MeV;         // e+- / gammas
// 240.MeV v. muon dE (1.2m) 0.240GeV/1.2m -> 5m/GeV



LDFTest::LDFTest() {
  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics);
  fEvent = 0;
}

LDFTest::~LDFTest() {
}

VModule::ResultFlag LDFTest::Init() {

  CentralConfig* cc = CentralConfig::GetInstance();

  Branch topB =
     cc->GetTopBranch ("SdSimpleSim");
  if (topB == 0) {
    ERROR ("Could not find branch SdSimpleSim");
    return eFailure;
  }


  Branch n_sampleB =  topB.GetChild ("n_sample");
  if (n_sampleB == 0) {
    ERROR ("Could not find branch n_sample");
    return eFailure;
  }
  n_sampleB.GetData (fNSample);






  Branch MuonRScaleB =  topB.GetChild ("MuonRScale");
  if (MuonRScaleB == 0) {
    ERROR ("Could not find branch MuonRScale");
    return eFailure;
  }
  MuonRScaleB.GetData (fMuonRScale);




  Branch em_ageB =  topB.GetChild ("em_age_fact");
  if (em_ageB == 0) {
    ERROR ("Could not find branch em_age");
    return eFailure;
  }
  em_ageB.GetData (fEMAgeFactor);




  return eSuccess;
}



VModule::ResultFlag LDFTest::Run (evt::Event& event) {

    fEvent++;

    INFO(" LDFTest -> Run ");

    // check event structure
    if (!event.HasSimShower())
        return eContinueLoop;

    evt::ShowerSimData& SimShower = event.GetSimShower();

    if (!SimShower.HasDirection())
        return eContinueLoop;

    if (!SimShower.HasPosition())
        return eContinueLoop;

    if (!SimShower.HasLongitudinalProfile())
        return eContinueLoop;

    if (!SimShower.HasGHParameters())
        return eContinueLoop;

    bool HasMuonProfile =
        SimShower.HasLongitudinalProfile (ShowerSimData::eMuon);

    bool HasGammaProfile =
        SimShower.HasLongitudinalProfile (ShowerSimData::ePhoton);

    bool HasElectronProfile =
        SimShower.HasLongitudinalProfile (ShowerSimData::eElectron);


    if (!HasMuonProfile ||
        !HasGammaProfile ||
        !HasElectronProfile) {
        ERROR ("Input file not suitable for LDFTest");
        return eContinueLoop;
    }


    // make new time
    const utl::TimeStamp& SimTime = SimShower.GetTimeStamp();
    utl::TimeStamp TriggerTime = SimTime;

    // get the detector definitions
    det::Detector& Det = det::Detector::GetInstance();
    const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();


    const atm::ProfileResult depthProfile
        = Det.GetAtmosphere().EvaluateDepthVsHeight();

    const atm::ProfileResult temperatureProfile
        = Det.GetAtmosphere().EvaluateTemperatureVsHeight();

    const atm::ProfileResult heightProfile
        = Det.GetAtmosphere().EvaluateHeightVsDepth();

/*
    const ProfileResult& slantDepthVsHeight =
        Det.GetAtmosphere().EvaluateSlantDepthVsHeight();


    const ProfileResult& heightVsSlantDepth =
        Det.GetAtmosphere().EvaluateHeightVsSlantDepth();

    const ProfileResult& distanceVsHeight =
        Det.GetAtmosphere().EvaluateDistanceVsHeight();
*/

    // get the simulated shower info
    double Xmax = SimShower.GetGHParameters().GetXMax();
    float simEnergy = SimShower.GetEnergy();
        const utl::Point& SimPos = SimShower.GetPosition();
    const utl::TabulatedFunction& SimProfile
        = SimShower.GetLongitudinalProfile();
    double Zenith = SimShower.GetZenith();
    double cosTheta = cos (Zenith);



    const utl::TabulatedFunction& SimMuonProfile =
        SimShower.GetLongitudinalProfile (ShowerSimData::eMuon);

    const utl::TabulatedFunction& SimGammaProfile =
        SimShower.GetLongitudinalProfile (ShowerSimData::ePhoton);

    const utl::TabulatedFunction SimElectronProfile =
        SimShower.GetLongitudinalProfile (ShowerSimData::eElectron);


    // Maximum profile slant-depth
    double ProfileMaxDepth = SimProfile [SimProfile. GetNPoints() - 1]. X();


    // calculate
    double HeightObs = wgs84.PointToLatitudeLongitudeHeight(SimPos).get<2>();


    // atmospheric parameters
    double XobsVert = depthProfile.Y (HeightObs);
    double Temperature = temperatureProfile.Y (HeightObs);
    static double g_earth = 9.81 * m/(s*s);
    double Pressure = XobsVert * g_earth;



    cout << " start " << endl;


    double rmin = 200*m;
    double rmax = 4100*m;
    int nBins = 50;
    double dr = (rmax-rmin) / (nBins);


    TGraph ldfEl (nBins);
    TGraph ldfGa (nBins);
    TGraph ldfMu (nBins);
    TGraph ldfSum (nBins);

    TGraph rhoEl (nBins);
    TGraph rhoGa (nBins);
    TGraph rhoMu (nBins);
    TGraph rhoSum (nBins);


    int iPoint = 0;
    for (double r=rmin; r<rmax; r+=dr) {

        // calcualte S1000 ////////////////////

        double Xobs = XobsVert / cosTheta;

        // if the profile range is exceeded, make a try with the last available profile bin
        // and if station did not triggered, we are fine - otherwise complain !
        bool ShowerProfileRangeExceeded = false;
        if (Xobs>ProfileMaxDepth) {

            ShowerProfileRangeExceeded = true;
            Xobs = ProfileMaxDepth;
        }


        double Nch=0;
        double Ne=0;
        double Ng=0;
        double Nmu=0;

        Nch = SimProfile.Y (Xobs); // charged particles
        if (HasMuonProfile) {
            Nmu = SimMuonProfile.Y (Xobs);
        } /*else {
            Nmu = Ne * .5;                // stupid muon number assumption
            }*/

        if (HasGammaProfile) {
            Ng = SimGammaProfile.Y (Xobs);
        } else {
            Ng = Ne * 10.;                // stupid gamma number assumption
        }

        if (HasElectronProfile) {
            Ne = SimElectronProfile.Y (Xobs);
        } else {
            Ne = Nch - Nmu;             // electron number assumption
            if (Ne<0) Ne = 0;
        }

        double Age = utl::ShowerAge (Xobs, Xmax);
        double Rm = utl::MoliereRadius (Temperature, Pressure, cosTheta);
        /*
        cout << " AGe " << Age << " rm " << Rm/m << " " << Ne << " " << Ng
             << " " << Nmu << " " << Xobs/g*cm*cm << endl;
        */
        double S = CalculateTankSignal (r, 1.8*m, 1.2*m, Zenith,
                                        Age, Rm,
                                        Ne, Ng, Nmu);


        if (fSmu) ldfMu.SetPoint (iPoint, (r), (fSmu));
        if (fSel) ldfEl.SetPoint (iPoint, (r), (fSel));
        if (fSga) ldfGa.SetPoint (iPoint, (r), (fSga));
        if (fSsum) ldfSum.SetPoint (iPoint, (r), (fSsum));

        rhoMu.SetPoint (iPoint, (r), (fRhomu));
        rhoEl.SetPoint (iPoint, (r), (fRhoel));
        rhoGa.SetPoint (iPoint, (r), (fRhoga));
        rhoSum.SetPoint (iPoint, (r), (fRhosum));

        iPoint++;


#ifdef DEBUG
        cout << " S/vem " << S << " Energy/EeV " << simEnergy/EeV << endl;
#endif

    }

    ldfSum.SetLineColor (1);
    ldfSum.SetLineWidth (2);
    ldfSum.SetMarkerColor (1);
    ldfMu.SetLineColor (2);
    ldfMu.SetLineWidth (2);
    ldfMu.SetMarkerColor (2);
    ldfEl.SetLineColor (3);
    ldfEl.SetLineWidth (2);
    ldfEl.SetMarkerColor (3);
    ldfGa.SetLineColor (4);
    ldfGa.SetLineWidth (2);
    ldfGa.SetMarkerColor (4);

    rhoSum.SetLineColor (1);
    rhoSum.SetLineWidth (2);
    rhoSum.SetMarkerColor (1);
    rhoMu.SetLineColor (2);
    rhoMu.SetLineWidth (2);
    rhoMu.SetMarkerColor (2);
    rhoEl.SetLineColor (3);
    rhoEl.SetLineWidth (2);
    rhoEl.SetMarkerColor (3);
    rhoGa.SetLineColor (4);
    rhoGa.SetLineWidth (2);
    rhoGa.SetMarkerColor (4);

    TLegend *l = new TLegend (.5,.6,.95,.9);

    //l->AddEntry (&ldfSum, "sum");
    l->AddEntry (&ldfMu, "muons");
    l->AddEntry (&ldfEl, "electrons");
    l->AddEntry (&ldfGa, "gammas");


    ostringstream cName;
    cName << "c_" << fEvent;
    TCanvas c (cName.str().c_str());
    c.Divide (2);

    TH2F axis1 ("axis1","axis",2,200,4100,2,0.001,200);
    TH2F axis2 ("axis2","axis",2,200,4100,2,0.001,200);

    gStyle->SetOptStat (0);
    gStyle->SetOptTitle (0);

    axis1.SetXTitle ("r/m");
    axis1.SetYTitle ("S/VEM");

    c.cd (1);
    gPad->SetLogy (1);
    gPad->SetGridx (1);
    gPad->SetGridy (1);
    axis1.Draw();
    ldfSum.Draw ("pl");
    ldfMu.Draw ("pl");
    ldfEl.Draw ("pl");
    ldfGa.Draw ("pl");
    l->Draw();

    c.cd (2);
    axis2.SetXTitle ("r/m");
    axis2.SetYTitle ("1/m^2");
    gPad->SetLogy (1);
    gPad->SetGridx (1);
    gPad->SetGridy (1);
    axis2.Draw();
    //rhoSum.Draw ("apl");
    rhoGa.Draw ("pl");
    rhoMu.Draw ("pl");
    rhoEl.Draw ("pl");


    TFile f ("ldf.root", "update");
    c.Write ("",TFile::kOverwrite);
    f.Close();

    return eSuccess;
}




VModule::ResultFlag LDFTest::Finish() {


  return eSuccess;
}


// ldf function from ldffinderOG
/*
inline double LDFTest::LDFFunction (const double r,
const double beta,
const double gamma) {

    if (LDFFinder::fgLDFType == LDFFinder::ePL) {
        const double nearCore = kNearRadius / k1000m;
        const double relR = r / k1000m;
        return pow(relR, beta + gamma*log(relR > nearCore ? relR : nearCore));
    } else
        return pow(r/k1000m, beta)*pow((700*meter + r)/(1700*meter), beta+gamma);
}
*/


double LDFTest::NKG (double N, double Rmoliere, double R, double s) {

    // normalize to Rmoliere
    R /= Rmoliere;

    double C = TMath::Gamma (4.5-s) /
        (TMath::Gamma (s) * TMath::Gamma (4.5-2*s));

    double rho = C * N/(2.*kPi*Rmoliere*Rmoliere);

    rho *= pow (R, s-2.);
    rho *= pow (1.+R, s-4.5);


    return rho;
}





double LDFTest::TankIntersection (double r,
                                  double phi,
                                  double z,
                                  double theta,
                                  double TankRadius) {

    double cosTheta = cos (theta);
    double sinTheta = sin (theta);

    double UntilBottom = 0;
    if (cosTheta!=0)
        UntilBottom = z / cos (theta);

    if (sinTheta==0)
        return UntilBottom;

    double UntilSide = 0;
    if (sinTheta!=0)
        UntilSide = (TankRadius + r * cos (phi) -
                     (TankRadius -
                      sqrt (TankRadius*TankRadius -
                            r * r * sin (phi) * sin (phi) )))
            / sin (theta);

    if (cosTheta==0)
        return UntilSide;


    return min (UntilSide, UntilBottom);
}




double LDFTest::LTP (double r,
                     double theta,
                     double energy) {

    double Rltp = 1;
    double Wltp = 1;

    return 1. / (1. + exp ((r-Rltp)/Wltp));
}



double LDFTest::SampleEnergy (double Emin,
                              double Emax,
                              double Age) {

    Emin /= MeV;
    Emax /= MeV;

    //double a0 = 1.0;
    double a1 = 6.42522 - 1.53183 * Age;
    double a2 = 168.168 - 42.1368 * Age;

    double F2max = NerlingF2 (Emin, a2, Age);

    double E;
    double accept;
    double r1, r2;
#ifdef DEBUG
    int N = 0;
#endif
    do {

        // draw random number from Nerling F1
        r1 = RandFlat::shoot(&fRandomEngine->GetEngine(), 0., 1.);
        r2 = RandFlat::shoot(&fRandomEngine->GetEngine(), 0., 1.);

        E = exp (r1 * log (Emax+a1) +
                 (1-r1) * log (Emin+a1)) - a1;

        // check with Nerling F2
        accept = NerlingF2 (E, a2, Age);
#ifdef DEBUG
        N++;
#endif

    } while (accept/F2max < r2);

    E *= MeV;

#ifdef DEBUG
    //cout << " >E: " << E/MeV << " " << Age << " " << N << endl;
#endif

    return E;
}




/*

      The nerling dN/dE = F1*F2

*/

double LDFTest::NerlingF1 (double E,
                           double a1) {

    return 1. / (E+a1);
}


double LDFTest::NerlingF2 (double E,
                           double a2,
                           double s) {

    return 1. / pow ((E+a2), s);
}



double LDFTest::CalculateTankSignal (double CoreDistance,
                                     double TankRadius, double TankHeight, double Zenith,
                                     double Age, double Rm,
                                     double Ne, double Ng, double Nmu) {

        double cosTheta = cos (Zenith);
        double sinTheta = sin (Zenith);

        double TankProjTop = kPi * TankRadius * TankRadius * cosTheta;
        double TankProjSide =  TankHeight * 2.* TankRadius * sinTheta;
        double TankProjArea = TankProjTop + TankProjSide;


        double RhoElectrons = NKG (Ne, Rm, CoreDistance, Age/fEMAgeFactor);
        double RhoGammas = NKG (Ng, Rm, CoreDistance, Age/fEMAgeFactor);
        double RhoMuons = NKG (Nmu, fMuonRScale, CoreDistance, Age);

/*
        double RhoElectrons = NKG (Ne, Rm/2.3, CoreDistance, Age);
        double RhoGammas = NKG (Ng, Rm/2.3, CoreDistance, Age);
        double RhoMuons = NKG (Nmu, fMuonRScale*4.5, CoreDistance, Age);
*/

        fRhoga = RhoGammas;
        fRhoel = RhoElectrons;
        fRhomu = RhoMuons;
        fRhosum = 0;


        /*
        int nElectrons = int (RandPoisson::shoot (fRandomEngine->GetEngine(),
                                                  RhoElectrons)
                              * TankProjArea);

        int nGammas = int (RandPoisson::shoot (fRandomEngine->GetEngine(),
                                               RhoGammas)
                           * TankProjArea);

        nMuons = int (RandPoisson::shoot (fRandomEngine->GetEngine(),
                                              RhoMuons)
                          * TankProjArea);
                          */

        double nElectrons = TankProjArea * RhoElectrons;
        double nGammas = TankProjArea * RhoGammas;
        double nMuons = TankProjArea * RhoMuons;

        cout << " Zenith " << Zenith/deg << endl;
        cout << " A " << TankProjArea/m/m << endl;
        cout << " rho   " << RhoElectrons << " " << RhoGammas << " " << RhoMuons << endl;
        cout << " nPart " << nElectrons << " " << nGammas << " " << nMuons << endl;


        fSga = 0;
        fSel = 0;
        fSmu = 0;
        fSsum = 0;

        double TankSignal = 0;
        int nSample;

        //fSmu += nMuons * kPi*TankRadius*TankRadius*TankHeight/TankProjArea / kMeterPerVEM;
        //TankSignal = fSmu;

        nSample = 100;//min (fNSample, nMuons);
        if (nSample) {
            double Weight = nMuons/nSample;
            for (int iSample=0; iSample<nSample; iSample++) {

              double p = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 1);

                double r, phi, z;
                if (p<TankProjTop/TankProjArea) { // hit tank-top

                    // sample on top
                    r = sqrt (RandFlat::shoot(&fRandomEngine->GetEngine(),
                                               0., TankRadius));

                    phi = RandFlat::shoot(&fRandomEngine->GetEngine(),
                                           0., 2.*kPi);

                    z = TankHeight;


                } else {                           // hit tank-side

                    // sample on side
                    r = TankRadius;

                    phi = RandFlat::shoot (&fRandomEngine->GetEngine(),
                                           -kPi/2., +kPi/2.);

                    z = RandFlat::shoot (&fRandomEngine->GetEngine(),
                                         0., TankHeight);

                }

                double TrackLength = TankIntersection (r,
                                                       phi,
                                                       z,
                                                       Zenith,
                                                       TankRadius);

                /*
                cout << " muon track/m: "
                     << TrackLength/m
                     << endl;
                */

                fSmu += Weight * TrackLength / kMeterPerVEM;
                TankSignal = fSmu;
            }
        } // if nMuon>0


        double nEM = nElectrons;// + nGammas;
        nSample = 100;//min (fNSample, nEM);
        if (nSample) {
            double Weight = nEM/nSample;
            for (int iSample=0; iSample<nSample; iSample++) {

                //double Eem = MeanEnergy;
                static double Emin = 1.e-2*MeV;
                static double Emax = 1.e+5*MeV;
                double Eem = SampleEnergy (Emin, Emax, Age);

                fSel += Weight * Eem  * kTrackPerGEV / kMeterPerVEM;
                TankSignal = fSel;

            }
        } // if nMuon>0

        nEM = nGammas;
        nSample = 100;//min (fNSample, nEM);
        if (nSample) {
            double Weight = nEM/nSample;
            for (int iSample=0; iSample<nSample; iSample++) {

                //double Eem = MeanEnergy;
                static double Emin = 1.e-2*MeV;
                static double Emax = 1.e+5*MeV;
                double Eem = SampleEnergy (Emin, Emax, Age);

                fSga += Weight * Eem  * kTrackPerGEV / kMeterPerVEM;
                TankSignal = fSga;

            }
        } // if nMuon>0


        fSsum = TankSignal;

        fSsum /= TankProjArea;
        fSga /= TankProjArea;
        fSel /= TankProjArea;
        fSmu /= TankProjArea;

        //cout << " signl " << fSga << " " << fSel << " " << fSmu << " " << fSsum << endl;

        return TankSignal;
}



/*
  See GAP note 2005-059
*/

static double fgTriggerThetaBins [5] = {15, 25, 35, 45, 55};
static double fgTriggerS1000Bins [4] = {4.5, 13.5, 40.5, 121.5};

static double fgTriggerPropShalf [5] [4] = {
    {1.9, 2.2, 3.6, 5.1},
    {1.9, 3.0, 3.5, 4.8},
    {2.5, 3.3, 4.1, 5.1},
    {3.5, 3.8, 4.2, 5.4},
    {4.2, 4.8, 4.8, 5.0}};

/*
static double fgTriggerPropShalfError [5] [4] = {
    {.6, .9, .4, .2},
    {.9, .5, .5, .3},
    {.7, .6, .5, .2},
    {.4, .5, .5, .2},
    {.9, .4, .4, .2}};
*/

double LDFTest::T1TriggerProbability (double signal, double S1000, double theta) {

    cout << " tgheta " << theta << endl;

    static double n = 2.8; // or 3

    if (S1000>fgTriggerS1000Bins [3]) {
        S1000 = fgTriggerS1000Bins [3];
    }

    if (S1000<fgTriggerS1000Bins [0]) {
        S1000 = fgTriggerS1000Bins [0];
    }

    if (theta>fgTriggerThetaBins [4]) {
        theta = fgTriggerThetaBins [4];
    }

    if (theta<fgTriggerThetaBins [0]) {
        theta = fgTriggerThetaBins [0];
    }

#ifdef DEBUG
    cout << " pT1: " << theta << " " << S1000 << " " << signal << endl;
#endif

    int thetaBin;
    for (thetaBin=0; thetaBin<4; thetaBin++) {
        if (theta>=fgTriggerThetaBins [thetaBin] && theta<=fgTriggerThetaBins [thetaBin+1]) {
            break;
        }
    }

    int s1000Bin;
    for (s1000Bin=0; s1000Bin<3; s1000Bin++) {
        if (S1000>=fgTriggerS1000Bins [s1000Bin] && S1000<=fgTriggerS1000Bins [s1000Bin+1]) {
            break;
        }
    }

    double theta_s1000_11 = fgTriggerPropShalf [thetaBin] [s1000Bin];
    double theta_s1000_12 = fgTriggerPropShalf [thetaBin] [s1000Bin+1];
    double theta_s1000_21 = fgTriggerPropShalf [thetaBin+1] [s1000Bin];
    double theta_s1000_22 = fgTriggerPropShalf [thetaBin+1] [s1000Bin+1];

    // interpolate in theta direction

    double delta_s1000 = (S1000 - fgTriggerS1000Bins [s1000Bin]) / (fgTriggerS1000Bins [s1000Bin+1]-fgTriggerS1000Bins [s1000Bin]);
    double s1000_1 = theta_s1000_11 + (theta_s1000_12-theta_s1000_11) * delta_s1000;
    double s1000_2 = theta_s1000_21 + (theta_s1000_22-theta_s1000_21) * delta_s1000;

    double delta_theta = (theta - fgTriggerThetaBins [thetaBin]) / (fgTriggerThetaBins [thetaBin+1]-fgTriggerThetaBins [thetaBin]);
    double S_half = s1000_1 + (s1000_2-s1000_1) * delta_theta;

#ifdef DEBUG
    cout << " S_0.5 " << S_half << endl;
#endif

    double p = pow (signal,n);
    p /= (p + pow (S_half,n));
    return p;
}



// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C ..
// End: