flyopunta/src/src_par/partgam.cpp

/***************************************************************************
 *   Copyright (C) 2005 by Giuseppe Pierazzini          Pisa 20.02.05      *
 *   pierazzini@unipi.it                                                   *
 *                                     F l y o                             *
 *   Epsi/NA48                                                             *
 ***************************************************************************/

#define EPSR     1.0E-200

#include "flyoh.h"

using namespace GmpConst;
using namespace std;

//-------------------------------------------------------------
//            P a r t P a i r    per pair prod.
//--------------------------------------------------------------
// Get the pair quadri momentum

PartGam::PartGam()
{
    type=Gamstar;
    prnorm=log(0.13497/0.0010201);
    Pair_Production=1;
    if (evento_.Gen.Event<1) Gout<<"\n Classe PartGam: "<<nome<<"*  done!";
}


void  PartGam::
DecParm_Pair()
{
// Attenzione questa routine vale solo per la
// produzione di  pair nella interazione dei gamma con la materia
// si suppone una distribuzione che va come 1/m**2;
// limiti tra 2 volte la massa dell'elettrone e 100-200 mev.
// mass electron= 0.00051 ,max  mass .100
//
//   P.print( nome);Gout<<" "<<fato<<" Dc "<<evento_.Event;
//   Qpatre.print( nome);Gout<<"  Dc "<<evento_.Event<< " dw "<<dw;

    if (fato!=Decay || dw==0)return ;
//  prnorm=log(0.100/(2.*melec)) =4.585367
    Qpatre.m= masspair=0.0010201*exp(Pran()*5.27851);
    Qpatre.Energia();
//  Qpatre.print("Pout");Gout<<"  Dc "<<evento_.Event;
// ritorna un quadrivettore massivo di pari impulso...
    SelRea->mspair=Qpatre.m;

// calcolo la direzione di decadimento dei figli.
    double cost = 1.0;
    while (1)
    {
        cost = 1.0 -  2.0 * Pran();

        if (Pran()<(0.66666+0.33333*cost*cost)) break;

// double csq=cost*cost;
// if(dw->Pran()<((1.-csq)*exp(2.*csq)/1.36)) break;  // forma piu' ok??


    }
    /***********/

    double   fi = DuePiGreco * Pran();
    dw->cstar=cost;
    dw->fistar=fi;
    dw->rg->cstar=-cost;
    dw->rg->fistar=fi+PiGreco;

    return;
}
//-----------------------------------------------------
//    P a i r   p r o d u c t i o n
//--------------------------------------------------
void   PartGam::
Get_Pair(Device *dev,double &path, double &rad_length)
{


    if (path<0.01) return;
    
    
    
    
    double xrd = - log(Pran()+EPSR);
    double trt = rad_length * 1.285714 * xrd;  // Lunghezza di radiazione  reale*9/7 //rivedre

if (trt<path)
    {
        path=trt;
        Pair_Production=2;   //pair pruduction done....
        SelRea->Pair_done++;
        dev->Pair_done++;
        fato=Decay;
        Ifato=-12;  // dead per pair production 

        /*    Gout<<setprecision(4)<<"\n Pair: "<<evento_.Event<<"  "<<dev->Get_Nome()<<"  "<<nome<<"  path "
             <<setw(8)<<path<<"  trt "
             <<setw(8)<<trt<<" rad "
             <<setw(8)<<rad_length
             <<" Pair_done= " <<SelRea->Pair_done;
        */
//    P.print( nome);Gout<<" "<<fato<<" Pr2 "<<evento_.Event;
        Particella *pb,*pc;
        if (Abilita_Brems>0 && dev->Brem_done<dev->maxBrems)
        {
            Partname=Elec;
            pb = new PartLepton();
            Partname=Elep;
            pc = new PartLepton();
        }
        else
        {
            Partname=Elec;
            pb = new Particella(); 
            Partname=Elep;
            pc = new Particella();
        }
// attenzione non alterare la sequenza dei link....altrimenti cambiano gli Id
        dw=pb;
        pb->up=this;
        pb->rg=pc;
        pc->lf=pb;
        pc->up=this;
        pb->Set_Id();
        pc->Set_Id(); 

    }
       

}