// $Id: RecoShowerTrack.cc,v 1.22 2006/03/28 14:26:23 naumov Exp $
// Author: Dmitry.Naumov   2005/10/11

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: RecoShowerTrack                                                      *
 *  Package: Energy                                                          *
 *  Coordinator: Dmitry.Naumov                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// RecoShowerTrack
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "RecoShowerTrack.hh"
#include "EConst.hh"
#include "SlastShowerGenerator.hh"
#include "ShowerTrack.hh"
#include "FluoCalculator.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "LightSourceFactory.hh"
#include "Ground.hh"
#include "O1_ClearSkyPropagator.hh"
#include "DetectorGeometry.hh"
#include "RadiativeFactory.hh"
#include "EnergyModule.hh"
#include "RecoEvent.hh"
#include "EEvent.hh"
#include "ETruth.hh"
#include "Atmosphere.hh"

#include <TMath.h>
#include <TCanvas.h>
#include <TProfile.h>
#include <TList.h>
#include <TPolyMarker3D.h>

using namespace sou;
using namespace TMath;
using namespace EConst;

ClassImp(RecoShowerTrack)

// constructors, destructor
//_____________________________________________________________________________
 RecoShowerTrack::RecoShowerTrack(): EsafMsgSource() {
    //
    // Constructor
    //
    Constructor();
}

//_____________________________________________________________________________
 RecoShowerTrack::RecoShowerTrack(Float_t t, Float_t p, Float_t h, Float_t tfov, Float_t pfov): EsafMsgSource() {
    //
    // Constructor
    //
    Constructor();
    SetThetaPhi(t,p);
    SetThetaPhiFoV(tfov,pfov);
    SetHmax(h);
    
}
//_____________________________________________________________________________
 void RecoShowerTrack::Constructor() {
    //
    // This is default ctor
    fHmax = -100;
    fShowerDir.SetXYZ(0,0,0);
    fEusoDir.SetXYZ(0,0,0);
    string fluoname         = Config::Get()->GetCF("Reco","EnergyModule")->GetStr("EnergyModule.fFluoCalculator");
    fFluoCalculator         = LightSourceFactory::Get()->GetFluoCalculator( fluoname );    
    fEusoAlt                = Config::Get()->GetCF("General","Euso")->GetNum("Euso.fAltitude")*km;
    fEnergyDistributionName = Config::Get()->GetCF("Reco","EnergyModule")->GetStr("EnergyModule.fEnergyDistributionName");
    fEsafSpectrum    = new EsafSpectrum();
    fShowerGenerator = new SlastShowerGenerator(kTRUE);
    fShowerGenerator->SetQuiet();
    fDetectorGeometry = new DetectorGeometry();
    Double_t radius   = Config::Get()->GetCF("Optics","ConicBaffle")->GetNum("ConicBaffle.fTopRadius")*mm;
    fDetectorGeometry->SetRadius(radius);
    fDetectorGeometry->SetPos(TVector3(0,0,fEusoAlt));
    fGround             = RadiativeFactory::Get()->GetGround();
    fClearSkyPropagator = new O1_ClearSkyPropagator(fGround);
    fClearSkyPropagator->CopyDetectorGeometry(fDetectorGeometry,false);
    fEnergyModule = NULL;
    fEnergy = 1.e21;
    fDebugList = new TList;
}

//_____________________________________________________________________________
 RecoShowerTrack::~RecoShowerTrack() {
    //
    // Destructor
    //
    SafeDelete(fShowerGenerator);
    SafeDelete(fDetectorGeometry);
    SafeDelete(fFluoCalculator);
    SafeDelete(fClearSkyPropagator);
    SafeDelete(fDebugList);
}
//_____________________________________________________________________________
 void RecoShowerTrack::Clear() {
    //
    //
    //
    for (UInt_t i(0); i < fSteps.size(); i++) {
        RecoShowerStep& rs = fSteps[i];    
        rs.fSpectrum->Delete(); 
    }
    fSteps.clear();
    for (Int_t i(0); i<fDebugList->GetEntries(); i++) {
         TObject* h = (TObject*)fDebugList->At(i);
         SafeDelete(h);
    }
    fDebugList->Clear();
    fDebugList = new TList;
    
}

//______________________________________________________________________________
 const RecoShowerStep& RecoShowerTrack::GetStepThetaPhi(Float_t Theta, Float_t Phi) const {
    // 
    // Return the step of the shower corresponding to the given time

    Int_t index = 0;
    Double_t mmin = -1.;
    TVector3 nFoV;

    nFoV.SetMagThetaPhi(1.,Theta,Phi);

    for (UInt_t i(0); i < fSteps.size(); i++) {
      const RecoShowerStep& rs = fSteps[i];
      EarthVector dir  = (rs.fDir).Unit();
      Double_t mag = dir*nFoV;
      if (mag > mmin) {mmin = mag; index = i;}
    }

    return (fSteps[index]);
}

//_____________________________________________________________________________
 Bool_t RecoShowerTrack::MakeTrack() {
    //
    //
    //
    Int_t progress = 0;
    Clear();
    if (fHmax == -100 || !fShowerDir.Mag() || !fEusoDir.Mag()) {
        if (fHmax == -100)     Msg(EsafMsg::Warning) << "Hmax is not setupn" << MsgDispatch;
        if (!fShowerDir.Mag()) Msg(EsafMsg::Warning) << "Shower Direction is not setupn" << MsgDispatch;
        if (!fEusoDir.Mag())   Msg(EsafMsg::Warning) << "Shower Direction is not setupn" << MsgDispatch;
                               Msg(EsafMsg::Panic)   << "Can not do a showern" << MsgDispatch;
    }

    // find 3D point of shower maximum

    Float_t AB = EarthRadius() + fEusoAlt;
    Float_t AC = EarthRadius() + fHmax;
    Float_t cs = Cos(fEusoDir.Theta());
    Float_t sn = Sin(fEusoDir.Theta());
    Float_t z = fEusoAlt - AB*cs*cs + cs*Sqrt(AC*AC - AB*AB*sn*sn);
    Float_t rho = (fEusoAlt-z)*Tan(fEusoDir.Theta());
    Float_t x   = -rho*Cos(fEusoDir.Phi());
    Float_t y   = -rho*Sin(fEusoDir.Phi());

    Float_t tMin = 0, tMax = 0;
    Float_t age, TotalYield, Attenuation;
    EarthVector MaximumPosition(x,y,z);
    fDepthMax = Atmosphere::Get()->Grammage(MaximumPosition,fShowerDir,"dir");
    EarthVector ISS(0,0,fEusoAlt);
    EarthVector fInitPos;
    // estimate for the current energy and assuming GIL parametrization the depth at the shower maximum (should be around 850 g/cm2)
    // and jump from the point of the shower maximum along the track direction for that in the atmosphere which is far by the shower
    // maximum depth
    Double_t depthToMaximum = fShowerGenerator->GetXmax();
    Int_t status = Atmosphere::Get()->InvertGrammage(MaximumPosition,fShowerDir,depthToMaximum,fInitPos,101*km);
    Msg(EsafMsg::Debug) << "<GetFirstPoint> InvertGrammage() status is " << status << MsgDispatch;

    fShowerGenerator->SetShower(fEnergy,fShowerDir.Theta(),fShowerDir.Phi(),fInitPos);
    fShowerGenerator->SetQuiet(kTRUE);
    fShowerGenerator->SetDepthStep(0.2);
    ShowerTrack *fShowerTrack = (ShowerTrack*)fShowerGenerator->Get();

    if (!fShowerTrack->Size()) return kFALSE; // do not process further if shower track object is empty

    // Define tMin and tMax of shower evolution

    {
      const ShowerStep& s = fShowerTrack->GetStep(0);
      Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
      EarthVector dir  = (ISS - 0.5*(s.GetXYZi() + s.GetXYZf()));
      tMin = (time + dir.Mag()/Clight())*1.e-3; // in microsecond
    }

    {
     const ShowerStep& s = fShowerTrack->GetStep(fShowerTrack->Size()-1);
      Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
      EarthVector dir  = (ISS - 0.5*(s.GetXYZi() + s.GetXYZf()));
      tMax = (time + dir.Mag()/Clight())*1.e-3; // in microsecond
    }

    Float_t GtuStep = 2.5;
    if (fEnergyModule) 
      GtuStep=fEnergyModule->GetRecoEvent()->GetHeader().GetGtuLength()*1.e-3;

    Int_t Nbins = (Int_t)Floor((tMax - tMin)/GtuStep);

    TH1F* h[10];
    for (Int_t i(0); i < 10; i++) {
         TString title = Form("EnergyModule_Step%d",i);
         TH1F* hh = NULL;
         if ((hh = (TH1F*)gDirectory->FindObject(title))) hh->Delete();
         h[i] = new TH1F(title,"",Nbins,0,tMax - tMin);
    }

    // fill the track, find new Hmax
    Double_t newHmax(-1), fluxMax(0);
    for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
      const ShowerStep& s = fShowerTrack->GetStep(i);

      // get position and compute for it: fluo yield + attenuation + time at EUSO
      EarthVector pos  = 0.5*(s.GetXYZi() + s.GetXYZf());
      Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
      Float_t     X    = 0.5*(s.GetXi() + s.GetXf());
      EarthVector dir  = ISS - pos;
      Float_t     EusoTime = (time + dir.Mag()/Clight())*1.e-3 - tMin; // in microsecond

      if (!(i % 50)) {
	age = 0.5*(s.GetAgei() + s.GetAgef());
	TF2* EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
	TF12 EnergyDistributionProj("EDS_name",EnergyDistribution,age,"X");
        Double_t temp = Atmosphere::Get()->Temperature(pos.Zv());
        TotalYield  = 0;
        Attenuation = 0;
        if (temp<=0) continue;
	TotalYield = fFluoCalculator->GetFluoYield(pos.Zv(),&EnergyDistributionProj,fEsafSpectrum);
	Attenuation=fClearSkyPropagator->GoToDetector(*fEsafSpectrum,pos);
      }

      Float_t  Length = (s.GetXYZi() - s.GetXYZf()).Mag();
      Double_t SolidAngle = 1./(4.*Pi()*dir.Mag2());
      Double_t flux  = SolidAngle*s.GetNelectrons()*Attenuation*TotalYield*Length;

      if (flux > fluxMax) {
          newHmax = pos.Zv();
          fluxMax = flux;
      } 

      h[0]->Fill(EusoTime,s.GetNelectrons());
      h[1]->Fill(EusoTime,flux);
      h[2]->Fill(EusoTime,pos.X()*flux);
      h[3]->Fill(EusoTime,pos.Y()*flux);
      h[4]->Fill(EusoTime,pos.Z()*flux);
      h[5]->Fill(EusoTime,EusoTime*flux);
      h[6]->Fill(EusoTime,TotalYield*flux);
      h[7]->Fill(EusoTime,Attenuation*flux);
      h[8]->Fill(EusoTime,age*flux);
      h[9]->Fill(EusoTime,X*flux);
      
      if (100.*(i+1)/fShowerTrack->Size() >= progress) {
	Msg(EsafMsg::Info).SetProgress(progress);
	Msg(EsafMsg::Info) << "MakeTrack():" << MsgCount;
	progress += 10;
      }
    }
    
    if (progress <= 100) Msg(EsafMsg::Info) << MsgDispatch;

    for (Int_t i = 1; i <= Nbins; i++) {
         Double_t flux   = h[1]->GetBinContent(i);
         if (!flux) continue;
         Double_t ne     = h[0]->GetBinContent(i);
	 Double_t x      = h[2]->GetBinContent(i)/flux;
         Double_t y      = h[3]->GetBinContent(i)/flux;
         Double_t z      = h[4]->GetBinContent(i)/flux;
         Double_t time   = h[5]->GetBinContent(i)/flux;
         Double_t yield  = h[6]->GetBinContent(i)/flux;
         Double_t atten  = h[7]->GetBinContent(i)/flux;
         Double_t age    = h[8]->GetBinContent(i)/flux;
	 Double_t X      = h[9]->GetBinContent(i)/flux;
         EarthVector pos(x,y,z);

	 EarthVector dir = ISS - pos;
	 TVector3 nFoV;
	 nFoV.SetMagThetaPhi(1.,dir.Theta(),dir.Phi());
	 Float_t length = 1.*Clight()*GtuStep*1.e+3/(1.+fShowerDir*nFoV);

	 RecoShowerStep rs;
         rs.fEusoTime    = time;
         rs.fPos         = pos;
         rs.fDir         = dir;
         rs.fFluoYield   = yield;
         rs.fAttenuation = atten;
         rs.fNe          = ne;
         rs.fFlux        = flux;
         rs.fLength      = length;
	 rs.fX           = X*cm2/g;

         TF2 *EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
         TF12 EnergyDistributionProj("EDS_name",EnergyDistribution,age,"X");
         Double_t temp = Atmosphere::Get()->Temperature(pos.Zv());
         if (temp<=0) continue;
         fFluoCalculator->GetFluoYield(pos.Zv(),&EnergyDistributionProj,
				       fEsafSpectrum);
         fClearSkyPropagator->GoToDetector(*fEsafSpectrum,pos.Zv());

	 rs.fSpectrum = new TH1F(Form("RecoShowerTrackSpectrum%d",i),"",
				 100,300*nm,400*nm);

	 for (Int_t nt = 0; nt < 1000; nt++) 
	   rs.fSpectrum->Fill(fEsafSpectrum->GetLambda());

         Add(rs);
    }

    for (Int_t i(0); i < 10; i++) h[i]->Delete();

    fShowerTrack->Reset();
    fShowerGenerator->Reset();

    Msg(EsafMsg::Debug) << "Hmax set up " << fHmax/km << " Hmax found " << newHmax/km << MsgDispatch;
    fHmax = newHmax;

    return kTRUE;
}

//_____________________________________________________________________________
 void RecoShowerTrack::DoDebugHistos() {
    //
    //
    //
    Float_t tMin, tMax;

    // Set tMin and tMax

    {
      RecoShowerStep& rs = fSteps[0];
      tMin = rs.fEusoTime;
    }

    {
      RecoShowerStep& rs = fSteps[fSteps.size()-1];
      tMax = rs.fEusoTime;
    }
    
    TH1* h[12];
    for (Int_t i(0); i < 11; i++) {
         TString title = Form("debug%d",i);
         TProfile* hh = NULL;
         if (( hh = (TProfile*)gDirectory->FindObject(title))) hh->Delete();
	 if (i >= 9)  h[i] = new TProfile(title,"",150,0,30);
         else         h[i] = new TProfile(title,"",150,tMin,tMax);
         fDebugList->Add(h[i]);
    }
    h[11] = new TH1F("debug11","",150,0,30);
    fDebugList->Add(h[11]);

    TPolyMarker3D *track = new TPolyMarker3D;
    fDebugList->Add(track);

    Double_t DetectorRadius = 1250*mm;
    Double_t Area = TMath::Pi()*DetectorRadius*DetectorRadius;
    EEvent  *SimuEvent  = fEnergyModule->GetRecoEvent()->GetSimuEvent();
    Double_t EnergyScale=1;
    if (SimuEvent)
        EnergyScale = SimuEvent->GetTruth()->GetTrueEnergy()/(fEnergy*1.e-6);

    Int_t points(0);
    fXmin=fYmin=fZmin=1.e10; 
    fXmax=fYmax=fZmax=-1.e10; 

    for (UInt_t i(0); i < fSteps.size(); i++) {
        RecoShowerStep& rs = fSteps[i];
        h[0]->Fill(rs.fEusoTime,rs.fNe);
        h[1]->Fill(rs.fEusoTime,rs.fPos.X()/km);
        h[2]->Fill(rs.fEusoTime,rs.fPos.Y()/km);
        h[3]->Fill(rs.fEusoTime,rs.fPos.Z()/km);
        h[4]->Fill(rs.fEusoTime,rs.fFluoYield*m);
        h[5]->Fill(rs.fEusoTime,rs.fAttenuation);
        h[6]->Fill(rs.fEusoTime,rs.fFlux);
        h[7]->Fill(rs.fEusoTime,rs.fLength/km);
        h[8]->Fill(rs.fEusoTime,rs.fX);

        h[9]->Fill(rs.fPos.Zv()/km,rs.fFluoYield*m);
        h[10]->Fill(rs.fPos.Zv()/km,rs.fAttenuation);
        Double_t fluoEP = rs.fFlux*Area*EnergyScale*Cos(rs.fDir.Theta());
        Int_t    nfluoEP = (Int_t)fluoEP;
        h[11]->Fill(rs.fPos.Zv()/km,fluoEP);
        if (nfluoEP) for (Int_t j(0); j<nfluoEP; j++) track->SetPoint(points++,rs.fPos.X()/km,rs.fPos.Y()/km,rs.fPos.Z()/km);
        if (fXmin > rs.fPos.X()/km) fXmin = rs.fPos.X()/km;
        if (fYmin > rs.fPos.Y()/km) fYmin = rs.fPos.Y()/km;
        if (fZmin > rs.fPos.Z()/km) fZmin = rs.fPos.Z()/km;
        if (fXmax < rs.fPos.X()/km) fXmax = rs.fPos.X()/km;
        if (fYmax < rs.fPos.Y()/km) fYmax = rs.fPos.Y()/km;
        if (fZmax < rs.fPos.Z()/km) fZmax = rs.fPos.Z()/km;

    }

}

//_____________________________________________________________________________
 void RecoShowerTrack::Debug(TCanvas *c, Int_t id, Option_t *opt) {
    //
    //
    //
    if (!c) {
        Warning("RecoShowerTrack::Debug","Please specify the pad where I have to draw");
        return;
    }
    
    if (id < 0 || id > 11) {
        Warning("RecoShowerTrack::Debug","Please specify id between 0 and 11");
        return;
    }
    
    TH1* h = (TH1*)fDebugList->At(id);
    h->Draw(opt);
}