TimeModelTest.cc

Go to the documentation of this file.

#include "TimeModelTest.h"
#include "MuonTimeModel.h"

#include <det/Detector.h>
using namespace det;

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
using namespace atm;

#include <utl/TabulatedFunction.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
using namespace utl;

#include <fwk/CentralConfig.h>
using namespace fwk;

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
using namespace evt;

#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TGraphErrors.h>
#include <TimeModel.h>

#include <iostream>
#include <sstream>
using namespace std;

using namespace TimeModelTestKG;

TimeModelTest::TimeModelTest() {}

TimeModelTest::~TimeModelTest() {}

VModule::ResultFlag TimeModelTest::Init() {
    
    TCanvas c;
    c.Print ("test.ps[");
    return eSuccess;
}

VModule::ResultFlag TimeModelTest::Finish() {

    TCanvas c;
    c.Print ("test.ps]");    
    return eSuccess;
}

VModule::ResultFlag TimeModelTest::Run (Event& event) {


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


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

    Detector& theDetector = Detector::GetInstance();

    Point  showerCore = simShower.GetPosition();
    Vector showerDir  = simShower.GetDirection();
    
    theDetector.GetAtmosphere().InitSlantProfileModel (showerCore,
                                                        showerDir,
                                                        5.*g/cm/cm);

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

    bool HasMuonProductionProfile = 
        simShower.HasLongitudinalProfile (ShowerSimData::eMuonProduction);



    utl::TabulatedFunction SimMuonProfile;
    if (HasMuonProfile) 
        SimMuonProfile = 
            simShower.GetLongitudinalProfile (ShowerSimData::eMuon);
    else {
        //ERROR ("No muon profile. Using: Nmu = .5*Ne !!!!!!!!!");
        ERROR ("No muon profile !!!!!!!!!");
    }

    utl::TabulatedFunction SimMuonProductionProfile;
    if (HasMuonProductionProfile)
        SimMuonProductionProfile = 
            simShower.GetLongitudinalProfile (ShowerSimData::eMuonProduction);
    else {
        ERROR ("No muon production profile  !!!!!!!!!");
    }


    int simPrimary = simShower.GetPrimaryParticle();
    double zenith = simShower.GetZenith();
    float simEnergy = simShower.GetEnergy();
    double cosTheta = cos (zenith);
    double logE = log10 (simEnergy/eV);

    int lorenzoPrimary;
    switch (simPrimary) {
        case utl::Particle::ePhoton:
            lorenzoPrimary = 0;
            break;
            
        case utl::Particle::eProton:
            lorenzoPrimary = 1;
            break;
            
        case utl::Particle::eIron:
            lorenzoPrimary = 2;
            break;
            
        default:
            ERROR ("Unknown LorenzoPrimary!");
            return eFailure;
            break;
    }           


    cout << " " << zenith << " " << lorenzoPrimary << endl;

    int nPoints = 20;
    double rMin =   100.*m;
    double rMax =  2000.*m;
    double dR = (rMax-rMin)/nPoints;


    /*
    TFile *file = new TFile (fname.c_str());
    if (!file->IsOpen()) {
        
        cerr << " file " << fname << " corrupt " << endl;
        return;
    }
    TTree *tree = (TTree*) file->Get ("Shower");
    TTree *header = (TTree*) file->Get ("Header");
    if (not tree || not header) {

        cerr << " tree in file " << fname << " corrupt " << endl;
        return;
    }
    float dMu [2000];
    float X [2000];
    int nX;
    tree->Print();
    header->Print();
    tree->SetBranchAddress ("dMu", dMu);
    header->SetBranchAddress ("X", X);
    header->SetBranchAddress ("nX", &nX);
    tree->GetEntry (0);
    header->GetEntry (0);
    
    TabulatedFunction dMu_tab;
    double BinWidth = 0;
    for (int i=0; i<nX; i++) {
        
        double SlantDepthMin = X [i] * g/cm/cm;
        if (i<nX-1) {
            double SlantDepthMax = X [i+1] * g/cm/cm;
            BinWidth = SlantDepthMax-SlantDepthMin;
        }
        
        cout << i << " " << (SlantDepthMin+BinWidth/2)/g*cm*cm
             << " " << dMu [i] << endl;
        dMu_tab.PushBack (SlantDepthMin+BinWidth/2, dMu [i]);
    }
    */

    int iPoint = 0;
    double r = 0;
    double xi = 0;

    cout << " model 1 " << endl;
    cout << " zenith " << zenith/deg 
         << " logE " << logE
         << " primary " << lorenzoPrimary 
         << endl;


    TGraphErrors *g1 = new TGraphErrors (nPoints);
    MuonTimeModel *mod1 = new MuonTimeModel (zenith, logE, lorenzoPrimary);
    mod1->Info();

    for (r=rMin; r<=rMax; r+=dR) {
        
        int n = 1000;
        
        mod1->SetCoordinates (r, xi);
        double time, timeErr = 0;
        mod1->GetMeanAndRMSOfFirstTime (time, timeErr, n);
        //time = mod1->GetFirstTime (n);
        //double time= mod1->GetMeanTime (n);
        time += tan(zenith)*r / kSpeedOfLight / ns;
                
        g1->SetPoint (iPoint, r/1000, time);
        g1->SetPointError (iPoint, 0, timeErr);
        iPoint ++;
    }

    delete mod1;



    cout << " model 2 " << endl;

    double XobsVert = 850.*g/cm/cm;
    /*
    MuonTimeModel *mod2 = new MuonTimeModel (zenith, dMu_tab, Xobs, 
                                     ShowerSimData::eMuonProduction);
    */

    TabulatedFunction logZdist = GetLogZdist (SimMuonProductionProfile, 
                                              cosTheta,
                                              XobsVert);

    /*
    cout << "--------------"<<endl;
    for (unsigned int iBin=0; iBin<logZdist.GetNPoints(); iBin++) {
        
        cout << " ->> " << logZdist[iBin].X() << " " << logZdist[iBin].Y() 
             << endl;
    }
    cout << "--------------"<<endl;
    */
    MuonTimeModel *mod2 = new MuonTimeModel (zenith, &logZdist,
                                             true, true);
    mod2->Info();

    TGraphErrors *g2 = new TGraphErrors (nPoints);

    iPoint = 0;
    for (r=rMin; r<=rMax; r+=dR) {
        
        int n = 1000;
        
        mod2->SetCoordinates (r, xi);
        double time, timeErr = 0;
        mod2->GetMeanAndRMSOfFirstTime (time, timeErr, n);
        //time= mod2->GetFirstTime (n);
        //double time= mod2->GetMeanTime (n);
        time += tan(zenith)*r / kSpeedOfLight;
        g2->SetPoint (iPoint, r/1000, time);
        g2->SetPointError (iPoint, 0, timeErr);
        iPoint ++;
    }
    
    delete mod2;
    
    
    TLegend *l = new TLegend (0.15, 0.7, 0.4, 0.99);
    l->AddEntry (g1, "using model", "l");
    l->AddEntry (g2, "using MuonProduction profile", "l");
    
    TCanvas c2 ("sdlf");
    g2->SetLineColor (2);
    g2->Draw ("al");
    g1->SetLineColor (4);
    g1->Draw ("l");
    l->Draw();
    c2.Print ("test.ps");
    
    static int iEvent = 1;
    ostringstream g1Name, g2Name;
    g1Name << "model_" << int(zenith/deg) << "deg_" << iEvent;
    g2Name << "dndz_" << int(zenith/deg) << "deg_" << iEvent;
    
    TFile file ("test.root", "UPDATE");
    g1->Write (g1Name.str().c_str(), TFile::kOverwrite);
    g2->Write (g2Name.str().c_str(), TFile::kOverwrite);
    file.Close();
    
    iEvent++;

    return eSuccess;
}


TabulatedFunction TimeModelTest::GetLogZdist (const TabulatedFunction& prof,
                                              double cosTheta,
                                              double XobsVert) {
    

    // convert X (slant depth) to z / m
    // taking care of earth curvature
    det::Detector& Det = det::Detector::GetInstance();
    
    const atm::ProfileResult heightProfile 
        = Det.GetAtmosphere().EvaluateHeightVsDepth();
    
    const atm::ProfileResult depthProfile 
        = Det.GetAtmosphere().EvaluateDepthVsHeight();
    
    const ProfileResult& heightVsSlantDepth = 
        Det.GetAtmosphere().EvaluateHeightVsSlantDepth();

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


    double atmDepthMin = heightVsSlantDepth.MinX();
    double atmDepthMax = heightVsSlantDepth.MaxX();

    //double HeightXobs = heightProfile.Y (XobsVert);

    vector <double> ValueZ;
    vector <double> ValuedNdZ;

    //double AtmDepthMin = depthProfile.Y (depthProfile.MaxX());
    //double AtmDepthMax = depthProfile.Y (depthProfile.MinX());


    // one bin in depth:
    double dXslant = abs (prof [1].X() - prof [0].X());
    
    for (unsigned int iBin=0; iBin<prof.GetNPoints(); iBin++) {
        
        double nDmu = prof [iBin].Y();
        double xSlant = prof [iBin].X();

        double xSlantLow  = xSlant-dXslant/2;
        double xSlantHigh = xSlant+dXslant/2;
        
        //cout << " slant " << xSlant/g*cm*cm;

        if (xSlantLow<atmDepthMin || xSlantHigh>atmDepthMax) {
/*
            cout << " continue "
                 << " " << dXslant/g*cm*cm
                 << " atmMin " << atmDepthMin/g*cm*cm
                 << " atmMax " << atmDepthMax/g*cm*cm
                 << endl;
*/
            continue;
        }
        
        double heightLow = heightVsSlantDepth.Y (xSlantHigh);
        double heightHigh = heightVsSlantDepth.Y (xSlantLow);
        double height = heightVsSlantDepth.Y (xSlant);
        double z = -distanceVsHeight.Y (height);
        double zLow = -distanceVsHeight.Y (heightLow);
        double zHigh = -distanceVsHeight.Y (heightHigh);
        
        
        //cout << " profile: height " << height << " ndmu " << nDmu << flush;
        
        
        
        double dz = abs (zHigh-zLow);
        
        //cout << " dz " <<  dz << flush;
        
        // dont't go below ground leel
        if (z<0) {
            //cout << endl;
            continue; 
        }

        // only look at shower not before
        if (nDmu==0) {
            //cout << endl;
            continue; 
        }
        
        // normalisation is not required !!! ????
        double logz = log10 (z);
        double dLogz = log10 (dz);
        double dMuLogz = (nDmu*dXslant/dLogz) / g*cm*cm * m;

        //cout << " -> log10(z) " << logz  << " " << dMuLogz << endl;

        ValueZ.insert (ValueZ.begin(), logz);
        ValuedNdZ.insert (ValuedNdZ.begin(), dMuLogz);
        
    }

/*
    for (int i=0; i< ValueZ.size(); i++ ) {

        cout << ValueZ [i] << " " << ValuedNdZ [i] << endl;
    }
*/

/*
    TabulatedFunction result (ValueZ, ValuedNdZ);
    
    cout << " ---------- result -------- " << endl;
    for (int i=0; i<result.GetNPoints(); i++) {
        
        cout << " ->> " << result[i].X() << " " << result[i].Y() << endl;
    }
    cout << " ---------- result -------- " << endl;
    return result;
*/

    return TabulatedFunction (ValueZ, ValuedNdZ);
    
}

/*
TabulatedFunction TimeModelTest::GetLogZdist (TabulatedFunction prof,
                                              double cosTheta,
                                              double XobsVert) {
    


    // convert X (slant depth) to z / m
    // taking care of earth curvature
    det::Detector& Det = det::Detector::GetInstance();
    
    const ProfileResult heightProfile 
        = Det.GetAtmosphere().EvaluateHeightVsDepth();
    const ProfileResult depthProfile 
        = Det.GetAtmosphere().EvaluateDepthVsHeight();

    double HeightXobs = heightProfile.Y (XobsVert);

    vector <double> ValueZ;
    vector <double> ValuedNdZ;

    double AtmDepthMin = depthProfile.Y (depthProfile.MaxX());
    double AtmDepthMax = depthProfile.Y (depthProfile.MinX());


    // one bin in depth:
    double dXslant = abs (prof [1].X() - prof [0].X());
    double dXvert = dXslant * cosTheta;

    for (unsigned int iBin=0; iBin<prof.GetNPoints(); iBin++) {
        
        double nDmu = prof[iBin].Y();
        double Xslant = prof[iBin].X();
        double Xvert = Xslant * cosTheta;


        if (Xvert-dXvert/2<AtmDepthMin || Xvert+dXvert/2>AtmDepthMax) {

            continue;
        }

        double Height = heightProfile.Y (Xvert);
        double z = (Height-HeightXobs)/cosTheta;

        double Height1 = heightProfile.Y (Xvert-dXvert/2);
        double z1 = (Height1-HeightXobs)/cosTheta;
        
        double Height2 = heightProfile.Y (Xvert+dXvert/2);
        double z2 = (Height2-HeightXobs)/cosTheta;
        
        double dz = abs (z1-z2);

        //cout << " dz " <<  dz << flush;

        // dont't go below ground leel
        if (z<0) {
            //cout << endl;
            continue; 
        }

        // only look at shower not before
        if (nDmu==0) {
            //cout << endl;
            continue; 
        }
        
        // normalisation is not required !!! ????
        double logz = log10 (z);
        double dLogz = log10 (dz);
        double dMuLogz = (nDmu*dXslant/dLogz) / g*cm*cm * m;

        //cout << " -> log10(z) " << logz  << " " << dMuLogz << endl;

        ValueZ.insert (ValueZ.begin(), logz);
        ValuedNdZ.insert (ValuedNdZ.begin(), dMuLogz);
        
    }

    return TabulatedFunction (ValueZ, ValuedNdZ);
    
}
*/


double TimeModelTest::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;
}