flyopunta/src/src_dev/devmsct.cpp

/***************************************************************************
                          gdevmsct.cpp  -  description
                             -------------------
    begin                : Sun Nov 26 2000
    copyright            : (C) 2000 by Giuseppe m Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************
 *                                                                         *
 *   NA48  simulation program.                                             *
 *                                                                         *
 ***************************************************************************/


/***   888 andrebbe messa con device.cpp   !!! ***/


#define ESTERNO  1
#define INTERNO  0
#define INBUCO   -1

#include "parm.h"
#include "apparato.h"
#include "particella.h"
#include "evento.h"

using namespace std;

//-------------------------------------------------
//              M S C A T T E R
//-------------------------------------------------

// =======  Multiple Scattering Routine =====
// at t e n z i o n e ... aggiustare per magnete e sfera!!!

int Device::
Mscatter (Particella * pr, double tratok)

{
    if(Scatter==0) return 0;     // no multiple scattering required
    lrad=nowpos?lradb:lradi;

    if (lrad <= 0.0             ||   // no multiscat.
            pr->charg == 0.0           ||   // no charge
            tratok < 0.001          )       // no crossed space
        return 0;

// qui lrad >0

// Cut on  Energia minore di qualche mev....se troppo basso potrebbe dare problemi
    double ee=pr->P.e;
    double Te=ee-pr->massa;
    double zz= nowpos?Zbuc:Zint;
    double Ec= .8/(1.2 +zz);   // Ec= 800 MeV/(Z+1.2)

    if (ee <Ec)
    {
//        pr->P.print("P1");
//     Gout<< " " <<zz<<" Te "<<Te<<" Ec "<< Ec;
//
        pr->P.setqn(0.,0.,0.,pr->massa);        // aggiorno il quadrivettore P
       Hits[nhit].e_rivela+=Te;  // aggiorno l'energia persa...
        pr->Ifato=-17;
        return -17;
    }

// calcolo dle multiplescatterig
    nrad= tratok/lrad;
    double px,py,pz;
    px=pr->P.x;
    py=pr->P.y;
    pz=pr->P.z;
    double p = pr->P.norma;

// Gout<<"\n Ms: ev  "<< setprecision(3)<<Eventi_Fatti << "  "<<nome <<"   " <<pr->nome
//       << "  E "<< eq<<" tratok "<<tratok
//      << "  rad. "<< lrad<< "  Nrad. "<< nrad ;

// get the p direction angles
    double pt = sqrt (px*px+py*py);
    sint = pt / p;
    cost = pz /p;
    if (pt > 0.)
    {
        sinf = py / pt;
        cosf = px / pt;
    }
    else
    {
        sinf = 0.;
        cosf = 1.;
    }





// get rotation matrix (versors)  to the p direction
    gvet Rx,Ry,Rz;
    Rx.setv(cosf * cost,  -sinf   ,   cosf * sint);
    Ry.setv(sinf * cost,   cosf   ,   sinf * sint);
    Rz.setv( -sint     ,    0.    ,          cost);


// get the teta scattering rms value (charge set to 1 by default)

    double z=pr->charg;
    tet0 = 0.0136 * z*z*sqrt (nrad) * (1. + 0.038 * log(nrad));
    if (pr->massa<0.001) tet0*=1./ee; // electrons....
    else tet0*=ee/(p*p);




    gvet S,DPX;
// generate scattering displacements and angles in x and y plane
    z1= Pgauss(1.,0.);
    z2= Pgauss(1.,0.);
    xplane=tratok*tet0*(z1*0.28867 + z2*0.5);
    txplane=p*z2*tet0;
    z1= Pgauss(1.,0.);
    z2= Pgauss(1.,0.);
    yplane=tratok*tet0*(z1*0.28867 + z2*0.5);
    typlane=p*z2*tet0;
// teta nello spazio ==>teta=sqrt(txplane*txplane+typlane*typlane);
// generate the scattered momentum vector taking as initial
// direction the  z axe ....... and...
    S.setv(txplane,typlane,0.0 );
// then Rotate to the final p direction
    DPX.setv(Rx%S,Ry%S,Rz%S);
    pr->P.x+=DPX.x;
    pr->P.y+=DPX.y;
    pr->P.z+=DPX.z;
    pr->Vers= pr->P.NVerso();         // memorizzo la nuova direzione.
    pr->P.setvn(p*pr->Vers, pr->massa);  //correzione (infinitesima!!) impulso al valore originario: la perdita è dopo.
// Coordinate corrections
// rotate to the initial P direction and add to
// the initial coordinates
    S.setv(xplane,yplane,0.0);
    DPX.setv(Rx%S,Ry%S,Rz%S); // multiple scatter disturb..
    pr->X+=DPX;

//debug
    if (Debugon==1)
    {
        Gout<<"\nSctr "
        <<setprecision(3)<<fixed
        <<setw(5)<<Eventi_Fatti<<" "
        <<setw(5)<< nome<<" "
        <<setw(4)<< pr->nome
        <<" t "<<setw(6)<<tratok
        <<" dx "<<setw(6)<<xplane
        <<" dy "<<setw(6)<<yplane;

    }
//-------------------------------------------------
//              D e g r a d o

// =======  Multiple Scattering energy loss =====

    int crepa=0;
// do nothing if ....
    if (Perdita ==0)                       return 0;  // no action
    if (  (Dedx=nowpos?Dedxbuc:Dedxint)==0) return 0; // no loss

    double persa= Dedx*tratok; //in GeV

//  double ee=pr->P.e;    Te=ee-pr->massa;              vedi più su
    if (persa>=Te)
    {
        persa=Te;
        pr->P.setvn(0.,0.,0.,pr->massa);
        crepa=-16;
    }
    else
    {
        ee-=persa;
        p=sqrt(ee*ee-pr->massa*pr->massa);
        pr->P.setvn(  pr->Vers*p ,pr->massa);


    }
    Hits[nhit].e_rivela+=persa; // aggiorno l'energia rivelata

//     pr->P.print("P2");
//     pr->X.print("X2");
//     Gout<<"\n\n "<<Eventi_Fatti<<" " <<pr->ido<<" " <<pr->P.e<<" persa " <<persa<<" ee  "<<ee<<" p "<<p<<"  cr " <<crepa<<"\n" ;
    return crepa;
}




double Device::
Msperdita (double tratto)
{
    if (  (Dedx=nowpos?Dedxbuc:Dedxint)==0) return 0; // no loss
// si suppone sempre perdita al minimo
    return Dedx*tratto;
}