source: trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc@ 1349

/* ======================================================================== *\
!
! *
! * 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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!   Author(s): Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

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

#include <TF1.h>
#include <TMath.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, MParticle *phot, TList *list)
{
    //
    //  gamma: primary particle from source
    //  phot:  infrared photon from background. (angle = interaction angle)
    //
    const Double_t E0     = 511e-6;              // [GeV]

    const Double_t theta  = phot->GetAngle();    // [2pi]
    const Double_t Ep     = phot->GetEnergy();   // [GeV]
    const Double_t Eg     = gamma->GetEnergy();  // [GeV]

    const Double_t ctheta = cos(theta);

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

    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;

    MElectron *p[2];
    p[0] = new MElectron(E+dE, gamma->GetZ());
    p[1] = new MElectron(E-dE, gamma->GetZ());

    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);

    p[0]->SetBeta(beta1);
    p[1]->SetBeta(beta2);

    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));

    p[0]->SetAngle(atan(tan1));
    p[1]->SetAngle(atan(tan2));

    list->Add(p[0]);
    list->Add(p[1]);

    return kTRUE;
}