/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////
#include "MPairProduction.h"

#include <TF1.h>
#include <TMath.h>
#include <TRandom.h>

#include "TList.h"
#include "MElectron.h"

ClassImp(MPairProduction);

Double_t AngleDistrib(Double_t *x, Double_t *k)
{

    const Double_t c = x[0]; // cos(alpha)
    const Double_t b = k[0]; // sqrt(1-1/s)

    const Double_t b2 = b*b;
    const Double_t b4 = b2*b2;

    const Double_t c2 = c*c;
    const Double_t c4 = c2*c2;

    const Double_t u = 1 - b4*c4 +2*b2*(1-b2)*(1-c2);
    const Double_t d = 1-b2*c2;

    return u/(d*d);
}

// --------------------------------------------------------------------------
MPairProduction::MPairProduction()
{
    fAngle = new TF1("AngleDistrib", AngleDistrib, -1, 1, 1);
}

MPairProduction::~MPairProduction()
{
    delete fAngle;
}

#include <iostream.h>

// --------------------------------------------------------------------------
Bool_t MPairProduction::Process(MParticle *gamma, const Double_t Ep, const Double_t theta, TList *list)
{
    //
    //  gamma: primary particle from source
    //  phot:  infrared photon from background. (angle = interaction angle)
    //  Ep: Energy photon [GeV]
    //  theta: interaction angle [rad]
    //
    const Double_t E0     = 511e-6;              // [GeV]
    const Double_t Eg     = gamma->GetEnergy();  // [GeV]

    const Double_t ctheta = cos(theta);

    const Double_t s      = Ep*Eg*2*(1-ctheta); //[1]
    const Double_t S      = s/(E0*E0*4); //[1]
    if (S<1)
        return kFALSE;

    const Double_t stheta = sin(theta);

    const Double_t sqrbetah = s/((Eg+Ep)*(Eg+Ep)) + 1;
    const Double_t sqrbetae = 1.-1./S;

    const Double_t GammaH = (Eg+Ep)/sqrt(s);

    const Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));

    fAngle->SetParameter(0, sqrt(sqrbetae));

    const Double_t alpha = psi-acos(fAngle->GetRandom());

    const Double_t salpha = sin(alpha);
    const Double_t calpha = cos(alpha);

    const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi)

    const Double_t bb = sqrt(sqrbetah/sqrbetae);

    const Double_t s1 = calpha/GammaH;
    const Double_t s2 = tphi*s1 - salpha - bb;

    const Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
    const Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;

    const Double_t E = (Eg+Ep)/2;;
    const Double_t f = sqrt(sqrbetah*sqrbetae)*salpha;

    // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush;

    MElectron &p0 = *new MElectron(*gamma, E*(1.-f), kTRUE);
    MElectron &p1 = *new MElectron(*gamma, E*(1.+f), kFALSE);

    Double_t rnd = gRandom->Uniform(TMath::Pi() * 2);
    p0.SetNewDirection(atan(+tan1), rnd);
    p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi());

    list->Add(&p0);
    list->Add(&p1);

    /*
    const Double_t Epg    = Ep/Eg; // 1+Epg ~ 1
    const Double_t Egp    = Eg/Ep; // 1+Egp ~ Egp

    const Double_t phi    = atan(sin(theta)/(Egp+ctheta));

    const Double_t sphi   = sin(phi);
    const Double_t cphi   = cos(phi);

    const Double_t alpha  = theta-phi;
    const Double_t calpha = cos(alpha);
    const Double_t salpha = sin(alpha);

    const Double_t beta = (Eg*cphi+Ep*calpha)/(Ep+Eg);

    //
    // gamma1 = 1/gamma = sqrt(1-beta*beta)
    //
    const Double_t gamma1 = sqrt((sphi*phi+Epg*salpha*Epg*salpha+2*Epg*(1-cphi*calpha)));

    const Double_t Beta   = sqrt(1-1/s);                //[1]

    fAngle->SetParameter(0, Beta);

    const Double_t psi    = atan(gamma1*sphi/(cphi-beta));
    const Double_t delta  = acos(fAngle->GetRandom()) - psi;

    const Double_t Bcosd  = Beta*cos(delta);
    const Double_t Bsind  = Beta*sin(delta);

    const Double_t E = sqrt(s)*E0/gamma1;
    const Double_t dE = E*Bcosd;

    const Double_t E1 = E0/(E+dE);
    const Double_t E2 = E0/(E-dE);

    const Double_t beta1 = sqrt(1.-E1*E1);
    const Double_t beta2 = sqrt(1.-E2*E2);

    const Double_t Bscp = Bsind*cphi;
    const Double_t spg  = sphi/gamma1;
    const Double_t cpg  = cphi/gamma1;

    const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp));
    const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp));

    MElectron &p0 = *new MElectron(E+dE, gamma->GetZ());
    MElectron &p1 = *new MElectron(E-dE, gamma->GetZ());
    p0 = *gamma; // Set Position and direction
    p1 = *gamma; // Set Position and direction

    p0.SetBeta(beta1);
    p1.SetBeta(beta2);

    static TRandom rand(0);
    Double_t rnd = rand.Uniform(TMath::Pi() * 2);
    p0.SetNewDirection(atan(tan1), rnd);
    p1.SetNewDirection(atan(tan2), rnd+TMath::Pi());

    list->Add(&p0);
    list->Add(&p1);
    */
    return kTRUE;
}