Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TrackDirectionModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TrackDirectionModule.cc,v 1.22 2005/10/21 13:30:08 moreggia Exp $
// R. Pesce created May, 11 2004

#include "RecoPixelData.hh"
#include "RecoEvent.hh"
#include "TrackDirectionModule.hh"
#include "EConst.hh"
#include "MedianFit.hh"
#include "LeastSquaresFit.hh"
#include "RecoRootEvent.hh"
#include "Config.hh"

#include <TVector3.h>
#include <TObject.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TGraph2D.h>
#include <TDirectory.h>
#include <TKey.h>
#include <TMath.h>
#include <TMinuit.h>
#include <TTree.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TGraph2D.h>

using namespace sou;
using namespace EConst;

ClassImp(TrackData)
ClassImp(TrackDirectionModule)

// ChiSquare function for beta fitting
void ChiSquareTrack(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag);
void MyChiSquareTrack(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag);
Double_t TimeAtEuso(Double_t *x, Double_t *par);

//_____________________________________________________________________________
 TrackDirectionModule::TrackDirectionModule() : RecoModule("TrackDirection"), 
    fTimeAtEuso(0) {
    // ctor
    fAll            = NULL;

}

//_____________________________________________________________________________
 TrackDirectionModule::~TrackDirectionModule() {
    // dtor
    delete fAll;
}

//_____________________________________________________________________________
 Bool_t TrackDirectionModule::Init() {
    fMinuitOutputLevel = (Int_t)Conf()->GetNum("TrackDirectionModule.fMinuitOutputLevel");
    if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3 )
        Msg(EsafMsg::Panic) << "Wrong config value TrackDirectionModule.fMinuitOutputLevel" << MsgDispatch;
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirectionModule::PreProcess() {
    Double_t kHUGE  = 1.e20;
    fAll            = NULL;
    fTheta          = -kHUGE;
    fPhi            = -kHUGE;
    fErrorTheta     = -kHUGE;
    fErrorPhi       = -kHUGE;
    fErrorDirection = -kHUGE;
    

    fNumPoints = 0;
    fNumHits = 0;
    fNum = -1;
    fBetaData.Clear();
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirectionModule::Process(RecoEvent *ev) {

    fEv = ev;
    const RecoModuleData *gcm = fEv->GetModuleData("GTUClustering");
    // check if GTUClusteringModule is not called and try to call HoughTransformModule
    if ( gcm == NULL ) {
         gcm = fEv->GetModuleData("HoughTransform");
    }
    if ( gcm == NULL ) {
        Msg(EsafMsg::Panic) << "TrackDirectionModule: not found GTUClusteringModule neither HoughTransformModule." << MsgDispatch;
    }
    else {
        if ( gcm->GetObj("CluPixels") == NULL ) {
            Msg(EsafMsg::Warning) << "TrackDirectionModule: no data in event" << MsgDispatch;
            return kFALSE;
        }
        fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
        if (fPointsId.size() == 0) {
            fNumPoints = 0;
            Msg(EsafMsg::Warning) << "TrackDirectionModule: no pixels selected: terminated." << MsgDispatch;
            return kFALSE;
        } else {
            fNumPoints = fPointsId.size();
            Msg(EsafMsg::Debug) << "TrackDirectionModule: processing " << fNumPoints << " points" << MsgDispatch;
        }
    }

    const RecoModuleData *tdpm = fEv->GetModuleData("TrackDetectorPlane");
    if ( gcm==NULL ) {
        Msg(EsafMsg::Panic) << "TrackDirectionModule: not found TrackDetectorModule." << MsgDispatch;
    } else {
        Double_t nx = tdpm->GetDouble("NormX");
        Double_t ny = tdpm->GetDouble("NormY");
        Double_t nz = tdpm->GetDouble("NormZ");
        Double_t wx = tdpm->GetDouble("WAxisX");
        Double_t wy = tdpm->GetDouble("WAxisY");
        Double_t wz = tdpm->GetDouble("WAxisZ");
        Double_t ux = tdpm->GetDouble("UAxisX");
        Double_t uy = tdpm->GetDouble("UAxisY");
        Double_t uz = tdpm->GetDouble("UAxisZ");
        Double_t xx = tdpm->GetDouble("XAxisX");
        Double_t xy = tdpm->GetDouble("XAxisY");
        Double_t xz = tdpm->GetDouble("XAxisZ");
        Double_t yx = tdpm->GetDouble("YAxisX");
        Double_t yy = tdpm->GetDouble("YAxisY");
        Double_t yz = tdpm->GetDouble("YAxisZ");
        fNormTDP.SetXYZ(nx,ny,nz);
        fWAxisTDP.SetXYZ(wx,wy,wz);
        fUAxisTDP.SetXYZ(ux,uy,uz);
        fXAxisTDP.SetXYZ(xx,xy,xz);
        fYAxisTDP.SetXYZ(yx,yy,yz);
        vector<TVector3> *spPoints = (vector<TVector3>*)tdpm->GetObj("SpPoints");
        fSpPos = *spPoints;
        fNumHits = tdpm->GetInt("NumHits");
    }
    fNum = fEv->GetHeader().GetNum();

    ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
    fHISS = pConfig->GetNum("Euso.fAltitude")*km;
    fEarthRadius = EarthRadius(); 

    FillBetaData();
    if (Conf()->GetStr("TrackDirectionModule.fMethod")=="approx") {
        UseApprox();
    }
    else if (Conf()->GetStr("TrackDirectionModule.fMethod")=="exact") {
        UseApprox();
        UseExact();
    } 
    else {
        Msg(EsafMsg::Panic) << "Wrong config value TrackDirectionModule.fMethod" << MsgDispatch;
    }
    PrintBeta();
    fTheta = TMath::Pi() - fTheta;   // DN FIXME! Roberto uses a definition of theta opposite to default
    MyData()->Add("Theta",fTheta);
    MyData()->Add("Phi",fPhi);
    
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirectionModule::PostProcess() {
    if ( fNumPoints == 0 ) return kTRUE;

    fCenterPos = fEv->GetHeader().GetTrueShowerMaxPos();
    Double_t dummyZ = fCenterPos.Z();
    fCenterPos.SetZ(-dummyZ+fHISS);
    fCenterPos = -fCenterPos;
    fCenterX = fCenterPos.X();
    fCenterY = fCenterPos.Y();
    fDist = TMath::Sqrt(TMath::Power(fCenterPos.X(),2)+TMath::Power(fCenterPos.Y(),2));
    fErrorPlane = TMath::ACos(fNormTDP*fTrueDir);
    fVisualAngle = TMath::ACos(fTrueDir*(fCenterPos.Unit()));

    fEASCenterX.push_back(fCenterPos.X());
    fEASCenterY.push_back(fCenterPos.Y());
    fDirectionErrors.push_back(fErrorDirection*TMath::RadToDeg());
    fBetaErrors.push_back(fErrorBeta*TMath::RadToDeg());
    fThetaErrors.push_back(fErrorTheta*TMath::RadToDeg());
    fPhiErrors.push_back(fErrorPhi*TMath::RadToDeg());
    fHMaxErrors.push_back(fErrorHMax/km);
    fNumFastOr.push_back(fEv->GetHeader().GetTrueNumFastOR());
    fNumPhEl.push_back(fNumHits);
    fPlaneErrors.push_back(fErrorPlane*TMath::RadToDeg());
    
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirectionModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
    fRecoRootEvent->GetRecoTrackDirection().SetTheta(fTheta);
    fRecoRootEvent->GetRecoTrackDirection().SetErrorTheta(fErrorTheta);
    fRecoRootEvent->GetRecoTrackDirection().SetPhi(fPhi);
    fRecoRootEvent->GetRecoTrackDirection().SetErrorPhi(fErrorPhi);
    fRecoRootEvent->GetRecoTrackDirection().SetErrorDirection(fErrorDirection);
    fRecoRootEvent->GetRecoTrackDirection().SetQuality(0); // Put here a quality variable.
    fRecoRootEvent->GetRecoTrackDirection().SetCheckInfo(fAll);
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t TrackDirectionModule::Done() {

    Int_t numEv = (Int_t)fDirectionErrors.size();
    if (numEv < 2) return kTRUE;

    Float_t *xc, *yc, *dErr, *tErr, *pErr, *hErr, *bErr;
    Float_t *fOr, *r, *plErr, *pEl;
    xc = new Float_t[numEv];
    yc = new Float_t[numEv];
    hErr = new Float_t[numEv];
    dErr = new Float_t[numEv];
    bErr = new Float_t[numEv];
    tErr = new Float_t[numEv];
    pErr = new Float_t[numEv];
    fOr = new Float_t[numEv];
    pEl = new Float_t[numEv];
    r = new Float_t[numEv];
    plErr = new Float_t[numEv];
    
    for( Int_t i=0; i<numEv; i++ ) {
        xc[i] = fEASCenterX[i];
        yc[i] = fEASCenterY[i];
        dErr[i] = fDirectionErrors[i];
        bErr[i] = fBetaErrors[i];
        tErr[i] = fThetaErrors[i];
        pErr[i] = fPhiErrors[i];
        hErr[i] = fHMaxErrors[i];
        fOr[i] = (Float_t)fNumFastOr[i];
        pEl[i] = (Float_t)fNumPhEl[i];
        r[i] = TMath::Sqrt(fEASCenterX[i]*fEASCenterX[i]+fEASCenterY[i]*fEASCenterY[i])/km;
        plErr[i] = fPlaneErrors[i];
    }

     // events statistics
     Int_t bestCounter(0);
     Int_t medCounter(0);
     Int_t worstCounter(0);
     for( Int_t i=0; i<numEv; i++ ) {
         if ( dErr[i] < 5.0 ) bestCounter++;
         else if ( dErr[i] >= 5.0 && dErr[i] <= 10.0 ) medCounter++;
         else 
             worstCounter++;
     }
     Msg(EsafMsg::Debug) << Form("STATISTICS on %d events",numEv) << MsgDispatch;
     Msg(EsafMsg::Debug) << Form("Events with reco error under 5 deg = %d (%.3f %)",
                  bestCounter,(Float_t)bestCounter/numEv*100.) << MsgDispatch;
     Msg(EsafMsg::Debug) << Form("Events with reco error between 5 and 10 deg = %d (%.3f %)",
                  medCounter,(Float_t)medCounter/numEv*100.) << MsgDispatch;
     Msg(EsafMsg::Debug) << Form("Events with reco error above 10 deg = %d (%.3f %)",
                  worstCounter,(Float_t)worstCounter/numEv*100.) << MsgDispatch;
 


    /*cout << "n-----Updating Error Histogram-----n" << endl;

    if ( !fGeneralDir ) {
        gDirectory->cd("/");
        cout << "GeneralDir undefined. Looking for the key..." << flush;
        TKey *kGeneralDir = gDirectory->FindKey("general");
        if ( !kGeneralDir ) {
            char pathname[] = "general";
            cout << " search failed.nCreating new "<< pathname << " directory..." << flush;
            fGeneralDir = new TDirectory(pathname,pathname);
        } else {
            cout << " succeded.nLoading object..." << flush;
            fGeneralDir = (TDirectory*)kGeneralDir->ReadObj();
        }

        cout << ( fGeneralDir ? " done" : " failed" ) << endl;
    }
    fGeneralDir->cd();

    
    cout << "Plotting " << numEv << " events..." << flush;
    
    TGraph2D *recoerr = new TGraph2D(numEv, xc, yc, dErr);
    recoerr->SetNameTitle("recoerr","Error in Reco Direction");
    recoerr->Write();
    delete recoerr;

    TGraph2D *betaerr = new TGraph2D(numEv, xc, yc, bErr);
    betaerr->SetNameTitle("betaerr","Error in Reco Beta");
    betaerr->Write();
    delete betaerr;

    TGraph2D *thetaerr = new TGraph2D(numEv, xc, yc, tErr);
    thetaerr->SetNameTitle("thetaerr","Error in Reco Theta");
    thetaerr->Write();
    delete thetaerr;
   
    TGraph2D *phierr = new TGraph2D(numEv, xc, yc, pErr);
    phierr->SetNameTitle("phierr","Error in Reco Phi");
    phierr->Write();
    delete phierr;
   
    TGraph2D *hmaxerr = new TGraph2D(numEv, xc, yc, hErr);
    hmaxerr->SetNameTitle("hmaxerr","Error in Reco HMax");
    hmaxerr->Write();
    delete hmaxerr;
   
    TGraph *fastor = new TGraph(numEv,fOr,dErr);
    fastor->SetNameTitle("fastor","RecoError vs NumFastOr");
    fastor->SetMarkerSize(.5);
    fastor->SetMarkerColor(kRed);
    fastor->SetMarkerStyle(23);
    fastor->Write();
    delete fastor;

    TGraph *phel = new TGraph(numEv,pEl,dErr);
    phel->SetNameTitle("phel","RecoError vs Num P.E.");
    phel->SetMarkerSize(.5);
    phel->SetMarkerColor(kRed);
    phel->SetMarkerStyle(23);
    phel->Write();
    delete phel;

    TGraph *dist = new TGraph(numEv,r,dErr);
    dist->SetNameTitle("dist","RecoError vs Dist from centre of FOV");
    dist->SetMarkerSize(.5);
    dist->SetMarkerColor(kRed);
    dist->SetMarkerStyle(23);
    dist->Write();
    delete dist;

    TGraph *plane = new TGraph(numEv,plErr,dErr);
    plane->SetNameTitle("plane","RecoError vs PlaneError");
    plane->SetMarkerSize(.5);
    plane->SetMarkerColor(kRed);
    plane->SetMarkerStyle(23);
    plane->Write();
    delete plane;

    TGraph *planefastor = new TGraph(numEv,fOr,plErr);
    planefastor->SetNameTitle("planefastor","PlaneError vs Num FastOr");
    planefastor->SetMarkerSize(.5);
    planefastor->SetMarkerColor(kRed);
    planefastor->SetMarkerStyle(23);
    planefastor->Write();
    delete plane;


    TGraph *planedist = new TGraph(numEv,r,plErr);
    planedist->SetNameTitle("planedist","PlaneError vs Dist from center");
    planedist->SetMarkerSize(.5);
    planedist->SetMarkerColor(kRed);
    planedist->SetMarkerStyle(23);
    planedist->Write();
    delete planedist;

    TGraph *planephel = new TGraph(numEv,pEl,plErr);
    planephel->SetNameTitle("planephel","PlaneError vs Num P.E.");
    planephel->SetMarkerSize(.5);
    planephel->SetMarkerColor(kRed);
    planephel->SetMarkerStyle(23);
    planephel->Write();
    delete planephel;

    TH1F *numErr = new TH1F("numErr", "Num Events vs RecoError", 100, 0, 180 );
    TH1F *plnumErr = new TH1F("plnumErr", "Num Events vs PlaneError", 100, 0, 180 );
    TH1F *nThetaErr = new TH1F("nThetaErr", "Num Events vs ThetaError", 100, -180, 180 );
    TH1F *nPhiErr = new TH1F("nPhiErr", "Num Events vs PhiError", 100, -180, 180 );
    for( Int_t i=0; i<numEv; i++ ) {
        numErr->Fill(dErr[i]);
        plnumErr->Fill(plErr[i]);
        nThetaErr->Fill(tErr[i]);
        nPhiErr->Fill(pErr[i]);
    }
    numErr->Write();
    plnumErr->Write();
    nThetaErr->Write();
    nPhiErr->Write();
    delete numErr;
    delete plnumErr;
    delete nThetaErr;
    delete nPhiErr;
*/
    delete [] xc;
    delete [] yc;
    delete [] dErr;
    delete [] tErr;
    delete [] pErr;
    delete [] hErr;
    delete [] bErr;
    delete [] fOr;
    delete [] r;
    delete [] pEl;
    delete [] plErr;

    return kTRUE;
}

//_____________________________________________________________________________
 void TrackDirectionModule::UserMemoryClean() {
    fAll = NULL;
}

//_____________________________________________________________________________
 void TrackDirectionModule::FillBetaData() {
    fBetaData.fNumPoints = fNumPoints;
    fBetaData.fNumHits = fNumHits;
    for( Int_t i=0; i<fNumPoints; i++ ) {
        RecoPixelData *pix = fEv->GetRecoPixelData(fPointsId[i]);
        fBetaData.fAlpha.push_back(TMath::ATan2(fSpPos[i]*fUAxisTDP,fSpPos[i]*fWAxisTDP));
        fBetaData.fTime.push_back((Double_t)2500.*pix->GetGtu()*ns);
        fBetaData.fHits.push_back(pix->GetCounts());
        fBetaData.fErrors.push_back(2500.*ns);
    }
    // search gtu min and max and with max nuber of hits
    Double_t tmin = fBetaData.fTime[0];
    Double_t tmax = fBetaData.fTime[0];
    Int_t maxhits = fBetaData.fHits[0];
    Int_t minid = 0;
    Int_t maxid = 0;
    Int_t mostpopid = 0;
    for( Int_t i=1; i<fNumPoints; i++ ) {
        if ( fBetaData.fTime[i] < tmin ) {
            tmin = fBetaData.fTime[i];
            minid = i;
        }
        if ( fBetaData.fTime[i] > tmax ) {
            tmax = fBetaData.fTime[i];
            maxid = i;
        }
        if ( fBetaData.fHits[i] > maxhits ) {
            maxhits = fBetaData.fHits[i];
            mostpopid = i;
        }
    }

    vector<Int_t> maxCounts;
    for( Int_t i=0; i<fNumPoints; i++ ) {
        if ( fBetaData.fHits[i] == maxhits ) maxCounts.push_back(i);
    }
    
    Double_t sum(0);
    for( size_t i=0; i<maxCounts.size(); i++ ) {
        sum += fBetaData.fAlpha[maxCounts[i]];
    }
    sum /= maxCounts.size();
    vector<Double_t> diffAlpha;
    for( size_t i=0; i<maxCounts.size(); i++ ) {
        diffAlpha.push_back(fBetaData.fAlpha[maxCounts[i]]-sum);
    }
    Int_t mpid(0);
    for( size_t i=1; i<maxCounts.size(); i++ ) {
        if ( diffAlpha[i] < diffAlpha[mpid] ) mpid=i;
    }
    mostpopid = maxCounts[mpid];
   
    Double_t alphamin = TMath::Abs(fBetaData.fAlpha[minid]);
    Double_t alphamax = TMath::Abs(fBetaData.fAlpha[maxid]);
    if ( alphamin < alphamax) {
        for( Int_t i=0; i<fNumPoints; i++ ) {
            if ( fBetaData.fTime[i] == tmin ) {
                if ( fBetaData.fAlpha[i] < alphamin ) {
                    alphamin = fBetaData.fAlpha[i];
                    minid = i;
                }
            }
            if ( fBetaData.fTime[i] == tmax ) {
                if ( fBetaData.fAlpha[i] > alphamax ) {
                    alphamax = fBetaData.fAlpha[i];
                    maxid = i;
                }
            }
        }
    } else if ( alphamax > alphamin ) {
        for( Int_t i=0; i<fNumPoints; i++ ) {
            if ( fBetaData.fTime[i] == tmin ) {
                if ( fBetaData.fAlpha[i] > alphamin ) {
                    alphamin = fBetaData.fAlpha[i];
                    minid = i;
                }
            }
            if ( fBetaData.fTime[i] == tmax ) {
                if ( fBetaData.fAlpha[i] < alphamax ) {
                    alphamax = fBetaData.fAlpha[i];
                    maxid = i;
                }
            }
        }
    }
    
    TVector3 pMax = fSpPos[mostpopid];
    fBetaData.fIdMostPop = mostpopid;
    fAlphaMax = fBetaData.fAlpha[fBetaData.fIdMostPop];
    fTMax = fBetaData.fTime[fBetaData.fIdMostPop];

    // MC truth
    fTrueEnergy = fEv->GetHeader().GetTrueEnergy();
    //Double_t zInit = fEv->GetHeader().GetTrueInitPos().Z();
    fTrueHMax = fEv->GetHeader().GetTrueShowerMaxPos().Z();
    fTrueThetaInc = fEv->GetHeader().GetTrueTheta();
    fTrueTheta = TMath::Pi()-fEv->GetHeader().GetTrueTheta();
    fTruePhi = fEv->GetHeader().GetTruePhi();
    fTrueDir.SetMagThetaPhi(1.,fTrueTheta,fTruePhi);
    fTrueBeta = TMath::ATan2(fTrueDir*fUAxisTDP, fTrueDir*fWAxisTDP);
    //TVector3 perp = fNormTDP.Cross(fTrueDir);
    //TVector3 initpos = fEv->GetHeader().GetTrueInitPos();
    
    //fRMax = CalculateRMax(fTrueHMax);
    fHMax = 5.*km;
    fRMax = CalculateRMax(fHMax);

    // create the function used to fit
    fTimeAtEuso = new TF1("timeteo", TimeAtEuso, -TMath::Pi(), TMath::Pi(), 5);
    // range's too wide. Get the correct interval
    if ( fBetaData.fAlpha[mostpopid] > 0)  
        fTimeAtEuso->SetRange(TMath::Pi()/4., 3./4.*TMath::Pi());
    else 
        fTimeAtEuso->SetRange(-3./4.*TMath::Pi(), -TMath::Pi()/4.);

    // set parameters
    fTimeAtEuso->SetParName(0,"h_max");
    fTimeAtEuso->SetParName(1,"beta");
    fTimeAtEuso->SetParName(2,"alpha_max");
    fTimeAtEuso->SetParName(3,"t_max");
    fTimeAtEuso->SetParName(4,"r_max");

    // fix const params

    fTimeAtEuso->SetParameter(0, 5.*km);
    fTimeAtEuso->SetParameter(1, fTrueBeta);
    fTimeAtEuso->FixParameter(2, fAlphaMax);
    fTimeAtEuso->FixParameter(3, fTMax);
    fTimeAtEuso->FixParameter(4, fRMax);


        
    Msg(EsafMsg::Debug) << Form("fTimeAtEuso range: alpha_max = %.3f, min = %.3f, max = %.3f",
                  fBetaData.fAlpha[mostpopid], fTimeAtEuso->GetXmin(), fTimeAtEuso->GetXmax()) << MsgDispatch;
    
    fBetaData.fIdMin = minid;
    fBetaData.fIdMax = maxid;
    fBetaData.fNorm = fNormTDP;
    fBetaData.fWAxis = fWAxisTDP;
    fBetaData.fUAxis = fUAxisTDP;
    fBetaData.fTimeAtEuso = fTimeAtEuso;

    fTimeAtEuso->Eval(fBetaData.fAlpha[mostpopid]);

    fBetaData.fCsi = TMath::ACos((pMax.Unit()).Z());

    //calculate omega mean
    fOmegaMean = 0;
    Int_t used = 0;
    for( Int_t i(0); i<fNumPoints; i++ ) {
        if (fBetaData.fTime[i]!=fTMax) {
            fOmegaMean += (fBetaData.fAlpha[i]-fAlphaMax)/(fBetaData.fTime[i]-fTMax);
            used++;
        }
    }
    fOmegaMean /= used;
    
    //calculate omega first and last
    fOmegaFirst = (fBetaData.fAlpha[minid]-fAlphaMax)/(fBetaData.fTime[minid]-fTMax); 
    fOmegaLast = (fBetaData.fAlpha[maxid]-fAlphaMax)/(fBetaData.fTime[maxid]-fTMax); 

}

//_____________________________________________________________________________
 void TrackDirectionModule::UseApprox() {
    // median fit in order to find an approximate value of beta to initialize fit
    LeastSquaresFit *lBeta = new LeastSquaresFit(fNumPoints, fBetaData.fAlpha, fBetaData.fTime, fBetaData.fErrors);
    MedianFit *mBeta = new MedianFit(fNumPoints, fBetaData.fTime, fBetaData.fAlpha);
    Double_t visual = 2*TMath::ATan(1/(mBeta->GetSlope()*fRMax/Clight()));
    fBeta = fBetaData.fAlpha[fBetaData.fIdMostPop] - visual - TMath::Pi();
    if ( fBeta > 1.e11 ) fBeta -= TMath::TwoPi() *1.e10;
    Int_t bPi = (Int_t) (fBeta/TMath::TwoPi());
    fBeta -= TMath::TwoPi()*bPi;
    fBetaInit = fBeta;
    fOmega = mBeta->GetSlope();
    fOmegaLeast = 1./lBeta->GetSlope();

    MsgForm(EsafMsg::Debug, "__Beta Pre-Fit_____________________________________________________");
    MsgForm(EsafMsg::Debug, "Slope = %.3e", mBeta->GetSlope()); 
    MsgForm(EsafMsg::Debug, "Offset = %.3e", mBeta->GetOffset());
    MsgForm(EsafMsg::Debug, "Absolute Deviation = %.3f", mBeta->GetAbsoluteDeviation());
    MsgForm(EsafMsg::Debug, "RMax init = %.3f km", fRMax/km);
    MsgForm(EsafMsg::Debug, "Visual Angle = %.3f deg", visual*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "Beta Init = %.3f deg", fBeta*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "____________________________________________________________________");
    
    fDir = TMath::Cos(fBeta) * fWAxisTDP + TMath::Sin(fBeta) * fUAxisTDP;
    fTheta = fDir.Theta();
    fPhi = fDir.Phi();
    if ( fPhi < - TMath::Pi() ) fPhi += TMath::Pi();
    if ( fPhi > TMath::Pi() ) fPhi -= TMath::Pi();

    fThetaInit = fTheta;
    fPhiInit = fPhi;

   
    // Error in reco direction
    fErrorDirection = TMath::ACos(TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
                      *TMath::Cos(fPhi)*TMath::Cos(fTruePhi) + 
                      TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
                      *TMath::Sin(fPhi)*TMath::Sin(fTruePhi) +
                      TMath::Cos(fTheta)*TMath::Cos(fTrueTheta));
    fErrorDirectionInit = fErrorDirection;
    fErrorBetaInit = fBetaInit - fTrueBeta;
    fErrorThetaInit = fThetaInit - fTrueTheta;
    fErrorPhiInit = fPhiInit - fTruePhi;

    MsgForm(EsafMsg::Debug, "__Angular Velocity__________________________________________________");
    MsgForm(EsafMsg::Debug, "Fitted (median) = %.3e rad/ns", fOmega);
    MsgForm(EsafMsg::Debug, "Fitted (least squares) = %.3e rad/ns", fOmegaLeast);
    MsgForm(EsafMsg::Debug, "Mean = %.3e rad/ns", fOmegaMean);
    MsgForm(EsafMsg::Debug, "First = %.3e rad/ns", fOmegaFirst);
    MsgForm(EsafMsg::Debug, "Last = %.3e rad/ns", fOmegaLast);
    MsgForm(EsafMsg::Debug, "____________________________________________________________________");

}

//_____________________________________________________________________________
 void TrackDirectionModule::UseExact() {
    
    // correction of HMax
    fHMax = 7.5*km * TMath::Log(1033./850./TMath::Abs(TMath::Cos(fTheta)));
    fRMax = CalculateRMax(fHMax);
    fTimeAtEuso->FixParameter(4, fRMax);

    
    // initialize first fit
    TMinuit *minbeta = new TMinuit(2);
    minbeta->SetPrintLevel(fMinuitOutputLevel);
    minbeta->SetObjectFit(&fBetaData);
    minbeta->SetFCN(MyChiSquareTrack);
    minbeta->SetErrorDef(1.);

    // initialize parameters
    Double_t par[2];
    
    // For debug purpose: initialize parameters with true values
    // par[0] = fTrueHMax;
    // par[1] = fTrueBeta;

    par[0] = fHMax;
    par[1] = fBeta;
    
    Double_t step[2] = {0.000, 0.001};
    Double_t min[2] = {0.0, 0.0};
    Double_t max[2] = {0.0, 0.0};
    string cpar[2] = {"H_max","Beta"};
    
    for( Int_t i=0; i<2; i++ ) {
        minbeta->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
    }
    
    //Double_t plist[1] = {2.};
    //Int_t ierflg;
    //minbeta->mnexcm("SET STRategy",plist,1,ierflg);
    //plist[0] = 1;
    //minbeta->mnexcm("FIX",plist,1,ierflg);
    minbeta->Migrad();
    
    // retrieve results
    Double_t outpar[2],err[2];
    for( Int_t i=0; i<2; i++ ) {
        minbeta->GetParameter(i,outpar[i],err[i]);
    }

    fBeta = outpar[1];
    fHMax = outpar[0];
    if ( fBeta > 1.e15 ) fBeta = fBetaInit;
    Int_t bPi = (Int_t) (fBeta/TMath::TwoPi());
    fBeta -= TMath::TwoPi()*bPi;

    MsgForm(EsafMsg::Debug, "n___First beta fit: hmax = %.3f km", fHMax/km);
    MsgForm(EsafMsg::Debug, "___First beta fit: beta = %.3f degn", fBeta*TMath::RadToDeg());
    
    fDir = TMath::Cos(fBeta) * fWAxisTDP + TMath::Sin(fBeta) * fUAxisTDP;
    fTheta = fDir.Theta();
    fPhi = fDir.Phi();
    if ( fPhi < - TMath::Pi() ) fPhi += TMath::Pi();
    if ( fPhi > TMath::Pi() ) fPhi -= TMath::Pi();
    
    // correction of HMax
    fHMax = 7.5*km * TMath::Log(1033./850./TMath::Abs(TMath::Cos(fTheta)));
    fRMax = CalculateRMax(fHMax);
    fTimeAtEuso->FixParameter(4, fRMax);

    // second fit
    par[0] = fHMax;
    par[1] = fBeta;
    
    for( Int_t i=0; i<2; i++ ) {
        minbeta->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
    }
    
    minbeta->Migrad();
    
    for( Int_t i=0; i<2; i++ ) {
        minbeta->GetParameter(i,outpar[i],err[i]);
    }

    fBeta = outpar[1];
    fHMax = outpar[0];
    if ( fBeta > 1.e11 ) fBeta = fBetaInit;
    bPi = (Int_t) (fBeta/TMath::TwoPi());
    fBeta -= TMath::TwoPi()*bPi;
    
    fDir = TMath::Cos(fBeta) * fWAxisTDP + TMath::Sin(fBeta) * fUAxisTDP;
    fTheta = fDir.Theta();
    fPhi = fDir.Phi();
    if ( fPhi < - TMath::Pi() ) fPhi += TMath::Pi();
    if ( fPhi > TMath::Pi() ) fPhi -= TMath::Pi();
    
    // Error in reco direction
    fErrorDirection = TMath::ACos(TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
                      *TMath::Cos(fPhi)*TMath::Cos(fTruePhi) + 
                      TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
                      *TMath::Sin(fPhi)*TMath::Sin(fTruePhi) +
                      TMath::Cos(fTheta)*TMath::Cos(fTrueTheta));                   
}

//______________________________________________________________________________
 Double_t TrackDirectionModule::CalculateRMax(Double_t hmax) {

    Double_t rmax;
    Double_t thmax = fEv->GetRecoPixelData(fPointsId[fBetaData.fIdMostPop])->GetTheta();
    rmax = (fEarthRadius+fHISS)*TMath::Cos(thmax)-
           TMath::Sqrt(TMath::Power(fEarthRadius+hmax,2)
              -TMath::Power((fEarthRadius+fHISS)*TMath::Sin(thmax),2));
    return rmax;
}

//______________________________________________________________________________
 void TrackDirectionModule::PrintBeta() {

    MsgForm(EsafMsg::Debug, "_EAS Direction Fit Results_____________________________________________________");
    MsgForm(EsafMsg::Debug, "    TrueEnergy = %.3e eV ", fTrueEnergy);
    MsgForm(EsafMsg::Debug, "    TrueDir = (%.3f, %.3f, %.3f)", fTrueDir[0],fTrueDir[1],fTrueDir[2]);
    MsgForm(EsafMsg::Debug, "    Dir = (%.3f, %.3f, %.3f)", fDir[0],fDir[1],fDir[2]);
    MsgForm(EsafMsg::Debug, "    Theta = %.3f deg   TrueTheta = %.3f deg", 
                 fTheta*TMath::RadToDeg(), fTrueTheta*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    Phi = %.3f deg   TruePhi = %.3f deg", 
                 fPhi*TMath::RadToDeg(), fTruePhi*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    Beta = %.3f deg   TrueBeta = %.3f deg",
                 fBeta*TMath::RadToDeg(), fTrueBeta*TMath::RadToDeg());

    MsgForm(EsafMsg::Debug, "    Error in Direction = %f deg", fErrorDirection*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    HMax = %.3f km   TrueHMax = %.3f km", fHMax/km, fTrueHMax/km);
    MsgForm(EsafMsg::Debug, "_______________________________________________________________________________");
    
    fThetaInc = TMath::Pi()-fTheta;
    fErrorTheta = fTheta - fTrueTheta;
    fErrorPhi = fPhi - fTruePhi;
    if ( fErrorPhi < -TMath::Pi() ) fErrorPhi += TMath::TwoPi();
    fErrorHMax = fHMax - fTrueHMax;
    fErrorBeta = fBeta - fTrueBeta;

    Double_t t_ref = fBetaData.fTime[fBetaData.fIdMostPop];
    Double_t t_prop, Dt, t_prop_tr, Dt_tr;

    TVector3 initpos = fEv->GetHeader().GetTrueInitPos();
    Double_t dummyZ = initpos.Z();
    initpos.SetZ(-dummyZ+fHISS);
    TVector3 trueMaxPos = fEv->GetHeader().GetTrueShowerMaxPos();
    dummyZ = trueMaxPos.Z();
    trueMaxPos.SetZ(-dummyZ+fHISS);
    trueMaxPos = -trueMaxPos;
    TVector3 perp = (fBetaData.fNorm.Cross(fDir)).Unit();
    TVector3 truenorm = (fTrueDir.Cross(initpos)).Unit();
    if ( truenorm.Z() < 0 ) truenorm *= -1;
    TVector3 trueperp = (truenorm.Cross(fTrueDir)).Unit();
    Double_t trueCsi = TMath::ACos((trueMaxPos.Unit()).Z());
    fTrueAlphaMax = TMath::ATan2(trueMaxPos*fUAxisTDP, trueMaxPos*fWAxisTDP);
    fTrueRMax = (-(fEarthRadius+fHISS)*TMath::Cos(trueCsi)-
                    TMath::Sqrt(TMath::Power(fEarthRadius+fTrueHMax,2)
                      -TMath::Power((fEarthRadius+fHISS)*TMath::Sin(trueCsi),2)));
    fTrueR0 = fTrueRMax*TMath::Sin(fTrueAlphaMax-fTrueBeta);
    fR0 = fRMax*TMath::Sin(fAlphaMax-fBeta);
    fAlpha0 = TMath::ATan2(perp*fBetaData.fUAxis, perp*fBetaData.fWAxis);
    fTrueAlpha0 = TMath::ATan2(trueperp*fBetaData.fUAxis, trueperp*fBetaData.fWAxis);
    Double_t dAlpha = fAlphaMax-fAlpha0;
    Double_t dAlpha_tr = fTrueAlphaMax-fTrueAlpha0;
    fTau0  = t_ref-TMath::Abs(fR0/Clight()*(1/TMath::Cos(dAlpha)))
                    +fR0/Clight()*TMath::Tan(dAlpha);
    fTrueTau0 = t_ref-TMath::Abs(fTrueR0/Clight()*(1/TMath::Cos(dAlpha_tr)))
                        +fTrueR0/Clight()*TMath::Tan(dAlpha_tr);
    //TVector3 testdir;
    //testdir.SetMagThetaPhi(1.,TMath::Pi()-fTrueDir.Theta(),fTrueDir.Phi());
    
    MsgForm(EsafMsg::Debug, "_______________________________________________________________________________");
    MsgForm(EsafMsg::Debug, "    Norm = (%.3f, %.3f, %.3f)", fNormTDP[0],fNormTDP[1],fNormTDP[2]);
    MsgForm(EsafMsg::Debug, "    TrueNorm = (%.3f, %.3f, %.3f)", truenorm[0],truenorm[1],truenorm[2]);
    MsgForm(EsafMsg::Debug, "    Perp = (%.3f, %.3f, %.3f)", perp[0],perp[1],perp[2]);
    MsgForm(EsafMsg::Debug, "    TruePerp = (%.3f, %.3f, %.3f)", trueperp[0],trueperp[1],trueperp[2]);
    MsgForm(EsafMsg::Debug, "    TrueMax = (%.3f, %.3f, %.3f)", trueMaxPos[0],trueMaxPos[1],trueMaxPos[2]);
    MsgForm(EsafMsg::Debug, "    RMax = %.3f km", fRMax/km);
    MsgForm(EsafMsg::Debug, "    True RMax = %.3f km", fTrueRMax/km);
    MsgForm(EsafMsg::Debug, "    R0 = %.3f km", fR0/km);
    MsgForm(EsafMsg::Debug, "    True R0 = %.3f km", fTrueR0/km);
    MsgForm(EsafMsg::Debug, "    test: t_max = %.3f microsecond", t_ref/microsecond);
    MsgForm(EsafMsg::Debug, "    Alpha_max = %.3f deg", fAlphaMax*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    True Alpha_max = %.3f deg", fTrueAlphaMax*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    Csi = %.3f deg", fBetaData.fCsi*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    True Csi = %.3f deg", trueCsi*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    Alpha_0 = %.3f deg", fAlpha0*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    True Alpha_0 = %.3f deg", fTrueAlpha0*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    dAlpha = %.3f deg", dAlpha*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    True dAlpha = %.3f deg", dAlpha_tr*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    tau0 = %.3f microsecond", fTau0/microsecond);
    MsgForm(EsafMsg::Debug, "    True tau0 = %.3f microsecond", fTrueTau0/microsecond);
    MsgForm(EsafMsg::Debug, "    Visual Angle = %.3f deg",(fAlphaMax-fBeta-TMath::Pi())*TMath::RadToDeg() );
    MsgForm(EsafMsg::Debug, "    True Visual Angle 1 = %.3f deg",(fTrueAlphaMax-fTrueBeta-TMath::Pi())*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "    True Visual Angle 2 = %.3f deg",(TMath::ACos((-trueMaxPos.Unit())*fTrueDir)-TMath::TwoPi())*TMath::RadToDeg());
    MsgForm(EsafMsg::Debug, "_______________________________________________________________________________");

    Double_t *alpha, *t_teo, *t_sper, *t_trueteo, *hits; 
    alpha = new Double_t[fNumPoints];
    t_teo = new Double_t[fNumPoints];
    t_sper = new Double_t[fNumPoints];
    t_trueteo = new Double_t[fNumPoints];
    hits = new Double_t[fNumPoints];

    TF1* trueTimeAtEuso = new TF1(*fTimeAtEuso);
    trueTimeAtEuso->SetParameter(0, fTrueHMax);
    trueTimeAtEuso->SetParameter(1, fTrueBeta);
    trueTimeAtEuso->SetParameter(2, fTrueAlphaMax);
    trueTimeAtEuso->SetParameter(3, t_ref);
    trueTimeAtEuso->SetParameter(4, fTrueRMax);

    fTimeAtEuso->SetParameter(0, fHMax);
    fTimeAtEuso->SetParameter(1, fBeta);

    for (Int_t i(0); i<fNumPoints ; i++) {
        dAlpha = fBetaData.fAlpha[i]-fAlpha0;
        dAlpha_tr = fBetaData.fAlpha[i]-fTrueAlpha0;
        t_prop = TMath::Abs(fR0/Clight()*1/TMath::Cos(dAlpha));
        Dt = fR0/Clight()*(-TMath::Sin(dAlpha))/TMath::Cos(dAlpha); 
        t_prop_tr = TMath::Abs(fTrueR0/Clight()*1/TMath::Cos(dAlpha_tr));
        Dt_tr = fTrueR0/Clight()*(-TMath::Sin(dAlpha_tr))/TMath::Cos(dAlpha_tr); 
        
        hits[i] = fBetaData.fHits[i];
        alpha[i] = fBetaData.fAlpha[i]*TMath::RadToDeg();
        t_sper[i] = fBetaData.fTime[i];
        
//        t_teo[i] = tau0 + Dt + t_prop;
//        t_trueteo[i] = tau0_tr + Dt_tr + t_prop_tr;
        t_teo[i] = fTimeAtEuso->Eval(fBetaData.fAlpha[i]);
        t_trueteo[i] = trueTimeAtEuso->Eval(fBetaData.fAlpha[i]);
    }


    // sort arrays
    Int_t j,k;
    Double_t a, tt, ts, ttt;
    
    for(k=1; k<fNumPoints; k++) {
    
        a = alpha[k];
        tt = t_teo[k];
        ts = t_sper[k];
        ttt = t_trueteo[k];
        j=k-1;
        while( j>=0 && alpha[j]>a ){
            alpha[j+1]=alpha[j];
            t_teo[j+1]=t_teo[j];
            t_sper[j+1]=t_sper[j];
            t_trueteo[j+1]=t_trueteo[j];
            j--;
        }
        alpha[j+1]=a;
        t_teo[j+1]=tt;
        t_sper[j+1]=ts;
        t_trueteo[j+1]=ttt;
    
    }
    
    TGraph *gTimeTeo = new TGraph(fNumPoints, alpha, t_teo);
    gTimeTeo->SetLineWidth(3);
    gTimeTeo->SetLineColor(kBlue);
//    gTimeTeo->SetMarkerStyle(20);

    TGraph *gTimeSper = new TGraph(fNumPoints, alpha, t_sper);
    gTimeSper->SetMarkerSize(.5);
    gTimeSper->SetMarkerColor(kRed);
    gTimeSper->SetMarkerStyle(23);

    TGraph *gTimeTrueTeo = new TGraph(fNumPoints, alpha, t_trueteo);
    gTimeTrueTeo->SetLineWidth(1);
    gTimeTrueTeo->SetLineColor(kGreen);

    if (!fAll)
        fAll = new TMultiGraph("BetaCheck","Comparison between t_sper and t_teo"); 
    else
        fAll->Clear();        
    fAll->Add(gTimeSper,"p");
    fAll->Add(gTimeTeo,"l");
    fAll->Add(gTimeTrueTeo,"l");
    
    
 /*   delete gTimeTeo;
    delete gTimeSper;
    delete gTimeTrueTeo;
 */   
    delete [] t_teo;
    delete [] t_trueteo;
    
//    TGraph2D *gHits = new TGraph2D(fNumPoints, alpha, t_sper, hits);
//    gHits->SetNameTitle("hits", "Hits of Pixels");
//    gHits->Write();
    delete [] alpha;
    delete [] t_sper;
    delete [] hits;
//    delete gHits;
}


//______________________________________________________________________________
/*void ChiSquareTrack(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
    f = 0;

    TrackData *tData = (TrackData*) gMinuit->GetObjectFit();
    
    Double_t t_teo;
    Double_t t_ref = tData->fTime[tData->fIdMostPop];
    Double_t alpha_ref = tData->fAlpha[tData->fIdMostPop];

    for(Int_t i(0); i<tData->fNumPoints; i++) {
        t_teo = t_ref - (tData->fHISS - par[0])/(TMath::Cos(tData->fCsi)*Clight()) *
                TMath::Sin(alpha_ref-par[1])/TMath::Tan(par[1]/2.) * (tData->fAlpha[i] - alpha_ref);
        f += tData->fHits[i]*TMath::Power(((tData->fTime[i]-t_teo)/tData->fErrors[i]),2);
    }
    
    f = (Double_t) f/(tData->fNumHits - 2.);
}*/

//______________________________________________________________________________
void ChiSquareTrack(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
    f = 0;
    TrackData *bData = (TrackData*) gMinuit->GetObjectFit();
    
    Double_t t_teo;
    //TVector3 dir = TMath::Cos(par[1])*bData->fWAxis+TMath::Sin(par[1])*bData->fUAxis;
    //TVector3 perp = bData->fNorm.Cross(dir);
    //Double_t alpha0 = TMath::ATan2(perp*bData.fUAxis, perp*bData.fWAxis);

    Double_t alpha_ref = bData->fAlpha[bData->fIdMostPop];
    Double_t alpha0 = par[1] + TMath::PiOver2();
    Double_t dAlpha = alpha_ref - alpha0;
    Double_t t_ref = bData->fTime[bData->fIdMostPop];
    Double_t R0 = (430.*km - par[0])/TMath::Cos(bData->fCsi)*TMath::Sin(alpha_ref-par[1]); 
    Double_t Dt = R0/Clight()*(-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);
    Double_t t_prop = TMath::Abs(R0/Clight()*(1)/TMath::Cos(dAlpha));
    Double_t tau0  = t_ref-(Dt+t_prop);
    //Double_t tau0 = t_ref-R0/Clight()*(1-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);

//    Double_t tau0  = t_ref-par[0]/Clight()*(1-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);

    for (Int_t i(0); i<bData->fNumPoints ; i++) {
        dAlpha = bData->fAlpha[i]-alpha0;
        t_prop = TMath::Abs(R0/Clight()*(1)/TMath::Cos(dAlpha));
        Dt = R0/Clight()*(-TMath::Sin(dAlpha))/TMath::Cos(dAlpha); 
        t_teo = tau0+Dt+t_prop;
        f+=bData->fHits[i]*TMath::Power(((bData->fTime[i]-t_teo)/bData->fErrors[i]),2);
    }
    f = (Double_t) f/(bData->fNumHits - 2.);
}

//______________________________________________________________________________
void MyChiSquareTrack(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
    f = 0;
    TrackData *bData = (TrackData*) gMinuit->GetObjectFit();
    
    Double_t t_teo;
    bData->fTimeAtEuso->SetParameter("h_max", par[0]);
    bData->fTimeAtEuso->SetParameter("beta", par[1]);
    
    for (Int_t i(0); i<bData->fNumPoints ; i++) {
        t_teo = bData->fTimeAtEuso->Eval(bData->fAlpha[i]) ;
        f+=bData->fHits[i]*TMath::Power(((bData->fTime[i]-t_teo)/bData->fErrors[i]),2);
    }
    f = (Double_t) f/(bData->fNumHits - 2.);
}

//_____________________________________________________________________________
void TrackData::Clear() {
    // TrackData clear
    fTime.clear();
    fAlpha.clear();
    fHits.clear();
    fErrors.clear();
    fNumPoints = 0;
    fNumHits = 0;
    fIdMin = 0;
    fIdMax = 0;
    fNorm.SetXYZ(0,0,0);
    fWAxis.SetXYZ(0,0,0);
    fUAxis.SetXYZ(0,0,0);
}

//_____________________________________________________________________________
Double_t TimeAtEuso(Double_t *x, Double_t *par){
    // par[0] = H_max
    // par[1] = Beta
    // par[2] = alpha_max (alpha_ref)
    // par[3] = t_max
    // par[4] = RMax
    // x[0] = alpha

    Double_t t_teo;

    Double_t alpha0 = par[1] + TMath::PiOver2();
    Double_t dAlpha = par[2]  - alpha0;
    Double_t R0 = par[4]*TMath::Sin(par[2] -par[1]); 
    Double_t t_prop = TMath::Abs(R0/Clight()*1/TMath::Cos(dAlpha));
    Double_t Dt = -R0/Clight()*TMath::Tan(dAlpha);
    Double_t tau0  = par[3]-(Dt+t_prop);

    /*
    cout << Form("_______________________________________________________________________________") << endl;
    cout << "dAlpha = " << dAlpha*TMath::RadToDeg()  << " deg"<< endl;
    cout << "R0 = " << R0/km << " km"<< endl;
    cout << "Dt = " << Dt/microsecond << " microsecond"<< endl;
    cout << "t_prop = " << t_prop/microsecond<< " microsecond" << endl;
    cout << "tau0  " <<  tau0/microsecond<< " microsecond" << endl;
    cout << "par[0] = "  << par[0]/km << " km" << endl;
    cout << "par[1] = " << par[1]*TMath::RadToDeg()  << " deg"  << endl;
    cout << "par[2] = " << par[2]*TMath::RadToDeg()  << " deg"  << endl;
    cout << "par[3] = " << par[3]/microsecond << " microsecond" << endl;
    cout << "par[4] = " << par[4]*TMath::RadToDeg()  << " deg"  << endl;
    cout << "par[5] = " << par[5]/km << " km"<< endl;
    cout << "x[0] = " << x[0] << endl;
    
    cout << Form("_______________________________________________________________________________") << endl;
     */
    dAlpha = x[0]-alpha0;
    t_prop = TMath::Abs(R0/Clight()*(1)/TMath::Cos(dAlpha));
    Dt = -R0/Clight()*TMath::Tan(dAlpha); 
    t_teo = tau0+Dt+t_prop;
    return t_teo;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005