flyopunta/src/src_par/partlepton.cpp

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

using namespace std;
using namespace GmpConst;

PartLepton::PartLepton()

{
    type=Lepstar;

    Ecut=0.250;
    // attenzione questo taglio determina l'energia minima di soglia per la produzion del Brems.
    // Nella distribuzione della energia di ionizzazione rilasciata dagli elettrini  nel radiatore appare un picco
    // alla stessa energia di soglia su scelta.
    // Questa soglia determina anche il numero delle particelle (gamma. elettroni e positroni)
    // generati nella simulazione per ciascun evento.

    double mrap= ptr_Elec->Get_Massa()/massa;
    brems_factor=mrap*mrap;  //brems_factor=0.01;
    if (Abilita_Brems) Brems_Production=1;
    else Brems_Production=0;

    if (evento_.Gen.Event<1)
    {
        Gout<<"\n Classe PartLepton  : "<<nome;
        if (Abilita_Brems) Gout<<"*  done!  with Brems_factor: "<<setprecision(6)<<brems_factor;
        else Gout<<"  done! ";
    }

}


PartLepton::~PartLepton()
{
}

int   PartLepton::
Get_Brems(Device *dev,double &path,double &rad)
{
// questa procedura e' chiamata da Device::traccia()
// se siamo qui....c' e' possibilita' di bremsstrhalung
// nell'interno o nel buco del dev..
//-----------------------------------------------
// vedi l'articolo di Y. S. Tsai, Pair Production and Bremsstrahlung of Charged Leptons,
// Rev. Moder. Phys. 46,1974

// procedura:
// 1- si accettano solo particelle con energia superiore a Ecut
// 2- in ogni intervallo dx=step=frazione di lrad si
//    verifica se viene emesso un gamma
// 3- se si, se ne genera l'energia
// 4- si aggiorna l'energia della particella madre e si
//    passa all step successivo.
// 5- l'energia di tutti i gamma emessi si sommano insieme
// 6- flyo poi genera un gamma di brems unico a partire
//    da un punto spaziale medio
// see the blue book or the article of Tsai del 1974 Rev.Moder. Phys. 46,4
//-------------------------------------------------
    double step,e,Emax, rapmx,rapmxq,rap,rapq,pmax;
    double eb=0.,prob,prand;
    double path_fatto,path_free,path_resta;
    double rad_length=rad;




//  rad deve essere > 0 !!

    if (path<.1) return 0;   //ritorna subito per cammini insignificanti


    if (massa>0.0005) rad_length/=brems_factor;
    if (path/rad_length<0.000001) return 0;

    step=rad_length*0.1;
    e=P.e;
    if (e<Ecut) return 0;  // attenzione a questo cut... vedi il costruttore Partlepton()
    ebrem=mbrem=0.0;
    path_fatto=0.0;
    path_free=0.0;

//  start loop
    while (path_fatto<path)
    {
        path_resta=path-path_fatto;
        if (step>path_resta) step=path_resta;


        e-=dev->Msperdita(step);       // stima della energia persa per msc.

        Emax=e-massa;
        rapmx= (Emax-Ecut)/e;
        rapmxq=rapmx*(Emax+Ecut)/e;
        pmax=(1.33333*(log(Emax/Ecut)-rapmx) + .5*rapmxq);

//test to make a brems

        if (Pran()<pmax*step/rad_length) // makes a gam...
        {
//  get gam energy  
//  con eb =(Ecut+alfa(Emax-Ecut)
//  rap = (eb-Ecut)/e = (Ecut+alfa(Emax-Ecut)-Ecut)/e = alfa*rapmx
            prand=Pran()*pmax;
            for (double alfa=0.0;alfa<1.;alfa+=0.001) // in step del permille
            {
                rap=alfa*rapmx;
                rapq=alfa*alfa*rapmxq;
                eb =Ecut+e*rap;
                prob=(1.33333*(log(eb/Ecut)-rap) + .5*rapq);
                if (prand<prob) break;
            }
            ebrem+=eb;
            e-=eb;
            path_free+=(path_fatto+Pran()*step)*eb;
        }

        if (e<=Ecut) break;   // vedere meglio????
        path_fatto+=step;
    }

    if (ebrem>0.0)      // bremsstrhalung ok...
    {
        Brems_Production=2;  // brems fatto....
        fato=Decay;// si fa per dire....si fa decare?!
        Ifato=-13;
        path_free/=ebrem;     //cammino libero prima di decadere: punto piu' probabile..
        SelRea->Brem_done++;
        dev->Brem_done++;
       
        if (Debugon==1)
        {
            P.print( nome);
            X.print(nome);
            Gout<<"\n Brlm "<<evento_.Gen.Event<<" "
            <<fato<< " dw "<<dw
            << " dev " << dev->Get_Nome()<<" Brem_done "<<dev->Brem_done
            << "  path " <<path<<"  free "<<path_free;
        }
        path=path_free;

// generate the new partcles: un nuovo leptone ed un gamma
// max defined into device cards...
        Particella *pb,*pc;
// electron o muon      test se  possono fare ancora brems...
        Partname=nome;
        if (Abilita_Brems>0&&dev->Brem_done<=dev->maxBrems) {
            pb = new PartLepton();
        }
        else pb= new Particella();  //no brems....

// gamma di brems...test se  può fare pair production ...
        Partname=Gam;
        if (Abilita_Pair>0&& dev->Pair_done<=dev->maxPair)   {
            pc = new PartGam();  // no pair production
        }
        else pc= new Particella();
// links: attention here the link is broken to the older particle... but see Restore_link
// 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();

        return -13;  //brems ok

    }
    return 0;  // no brems
}
//-------------------------------------------------
//     B r e m s s t r h a l u n g
//     produzione di gamma da elettroni
//--------------------------------------------------------------
// Get the brems mass  (attenzione al decadimento)

void  PartLepton::
DecParm_Brems()
{
// Attenzione questa routine vale solo per la produzione
// di gamma da bremsstrahlung
// Procedura
// Si ridefiniscono i parametri della particella che fa bresstrhalung
// immaginandola una particella con massa virtuale che decade in lepton piu gamma
// di bremsstrhalung di energia ebrem definata in Get_Brems().
//-------------
// (1)
// calcolo la massa invariante del sistema e-gam con forma esponenziale..
// partendo da un valore limite minimo che tiene conto dell'energia del gamma emesso
// (2)
// ridefinizione del quadrimpulso della particella padre
// si aggiorna la massa (virtuale) e si ricalcola l'energia
// (3)
// nota l'energia del gamma nel lab, calcolo il coseno* di decadimento
// compatibile con l'energia (ebrem) del gamma nel lab.
    if (Debugon==1)
    {
        P.print( "Qpatre");
        Gout<<" Brem_ev "<<evento_.Gen.Event;
    }

    if (fato!=Decay || dw==0) return ;
    double Et,Ele,pe,masmin;
    Et= Qpatre.e;
    Ele=Et-ebrem;

// caso perdita energia per multiple scattering aggiusta...
    if (Ele<Ecut)
    {
//        Gout<<"\n *** ERRORE *** Ev "<<evento_.Gen.Event<<setprecision(4)<<" ebrem "<< ebrem<<"  Ele "<<Ele<<"Q.e  "<< Qpatre.e;
        Ele=Ecut;
        ebrem=Et-Ecut;
    }


    pe=sqrt(Ele*Ele-massa*massa);
    masmin = sqrt( massa*massa + 2.*(Ele -pe)*ebrem);

//   Gout<<"\nDecBrem:  Ev "<<evento_.Event<<"  Ele "<<Ele<<" ebrem="<<ebrem<<" ms="<<masmin;

// genero la massa a random esponenzialmente partendo dalla minima... ...piccola
//    mbrem= -0.001*log(Pran()+EPSI)+masmin;  //larghezza ... mev
    mbrem= -0.01*log(Pran()+EPSI)+masmin;  //larghezza ... mev
    mbremq= mbrem*mbrem;
    SelRea->mbrem=mbrem;  //ricorda l'energia in reazione...
    Qpatre.m=mbrem;
    Qpatre.mq=mbremq;
// piccola correzione (rinormalizzazione)  dell 'impulso... della particella padre
// a energia totale costante....

    Qpatre.Impulso();


// calcolo il coseno* di decadimento
    double  gamma = Qpatre.e/Qpatre.m;
    double  beta  = Qpatre.norma/Qpatre.e;

    double    pstar = 0.5*(mbremq-massa*massa)/mbrem;   // = Estar del gamma

//    double cost = (Ele-ebrem-gamma*massa*massa/mbrem)/(2.*beta*gamma*pstar);
    double cost = -(ebrem-gamma*pstar)/(beta*gamma*pstar);

//      Gout<<"\n ebrem "<< ebrem<<"  mbrem "<<mbrem<<"  pstar "<<pstar<<"  cos "<<cost;
    /*
    //if(evento_.Event>83500)
    {
       char sbuf[128];
       sprintf(sbuf,"\nDecBrem: Ev %6d id  %2d ET %10.6f Ele %10.6f pe %10.6f Min  %8.6f  mb %8.6f  eb %10.6f cs  %8.6f z %8.1f ",
                       evento_.Event,id,Et,Ele,pe,masmin,mbrem,ebrem,cost,X.z);
       Gout<<sbuf;
       P.print(evento_.Event,nome,id,"DecBrem");
       Qpatre.print(evento_.Event,nome,id,"qPatre");Gout<<endl;
    }
    */

// Ridefinisco Pbar....Attenzione se si dimentica si sbaglia lorentz....
    Pbar=Qpatre;
//----------------------------------------------------------------
// memorizzo la direzione di decadimento dei figli.
    double   fi = DuePiGreco * Pran();
    dw->cstar=cost;
    dw->fistar=fi ;             //fi  elettrone
    dw->rg->cstar=-cost;
    dw->rg->fistar=fi+PiGreco; //fi+pigrec  gam
}