|
FlyoDoc_2011 Pisa 2011 by GmP --- 011
|
|
FlyoDoc_2011 Pisa 2011 by GmP --- 011
|
#include <devmagn.h>
Public Member Functions | |
| double | Get_CamMag () |
| virtual int | Traccia (Particella *) |
| ================ T R A C C I A ==================================== | |
| virtual void | Prgeom () |
Public Attributes | |
| int | iter_mag |
| double | step |
| double | dse |
| double | path |
| double | rho |
| double | freestep |
| double | vol_max |
| double | cam_max |
Protected Attributes | |
| gvet | B |
| gvet | X |
| gvet | XIN |
| gvet | PIN |
| gvet | VSIN |
| gvet | Puls |
| gvet | KIK |
| gvet | Pvers |
| gvet | Dvers |
| int DevMagn::Traccia | ( | Particella * | pr | ) | [virtual] |
================ T R A C C I A ====================================
Traccia le particelle nell'apparato per un tratto come definito in Get_camm(). Nel tracciare la particella si tiene conto della possibilita' della pair production e del bremsstrhalung.
====================================================================================
Reimplemented from Device.
Definition at line 60 of file devmagn.cpp.
{
// 18/03/03 update: general tracing.. the magnetic field now could be variable
// and without any restriction..
// *****attention! inside the magnet no brems, nor pairprod for the moment..
// the track tracing is iterative with step size small enough to keep the
// tracing errors to very low values
//
//*****
// La carica elementare nel sistema Ibrido e'
// ec=.299792458[ Gev/c * T**(-1) m **(-1) ]
// va qui moltiplicata per 10e-6 poiche'
// si usa il Gauss ed il Centimetro; segue...
// r(cm)=[p cos(lambda)(GeV)]/(2.9979258e-7 * H(Gauss) )
//**
// note 2007 the field is in the hole!
const double ec=2.99792458e-7;
if(nowpos==ESTERNO) return Device::Traccia(pr); // se esterno traccia regolarmente...
else if(nowpos==INTERNO)
{ // se interno (!no nel buco) traccio regolarmente per i mu.. le altre particelle sono perse...
if (pr->Get_Nome()==Mum||pr->Get_Nome()==Mup) return Device::Traccia(pr);
else return -1; // traccia lost in the frame
}
if (pr->P.norma == 0.0) return 0; // no motion ....
if (pr->Get_Charg() == 0.0) return Device::Traccia(pr); // neutra traccio regolamente dritto
//------------ tracing in the field....
X=Lab2cDev(pr->X); // Coordinata di impatto sul dev
PIN=Puls=Lab2Dev(pr->P); // impulso vettore in device....
VSIN=Pvers=Puls.Verso();
XIN=X; // device frame coordinates
vol_max=pr->pathok; // cammino potenziale prima di decadere!
P_Cin=X-Centrin; //Occorre in Get_CamMag(), specie per buchi non simmetrici
freestep=cam_max=Get_CamMag(); // cammino potenziale in linea retta per uscire dal magnete
if(freestep>vol_max) freestep=vol_max;
/* Campo in Field
campo=0 B field costante (va definito nella classe)
campo>0 B field variabile
*/
if(campo==0)
{
if (Field.norma == 0.0)
{
Gout<<"\n *** error*** Standard magnet "<<nome<<" null field???";
exit(0);
}
}
// se variabile ridefinisco il campo iniziale
else Get_Field(P_Cin);
// nota il sistema di riferimento dei vettori e' quello del dev.
// predico il raggio di curvatura.... massimo...
rho= 9999999. ;
if(Field.norma>0.0001)rho = Puls.Norma()/(ec*Field.norma);
// il calcolo corretto e' rho= (P&H).Norma/(H_norma**2*ec)
if(Debugon==1)
// if( (idev==9||idev==10)&&pr->Get_Idm()==14)
{// ok nel caso di campo costante ==> valore significativo
double forkik=ec*Field.norma*cam_max*pr->Get_Charg();
Gout<<"\n\n Traccia Ev="
<<setprecision(0)<<fixed
<<setw(6)<<Eventi_Fatti<<" "
<<setw(6)<<pr->Get_Nome()<<" in"
<<setw(6)<<nome<<" pth="
<<setprecision(1)<<fixed
<<setw(7)<<pr->pathok<<" Cmx="
<<setw(7)<<cam_max<<" rho="
<<setw(7)<<rho<<" kik="
<<setprecision(4)<<fixed
<<setw(6)<<forkik;
Field.print("Field");
X.print("X_in");
PIN.print("p_in");
}
/*************/
if(rho<2.*freestep) return 1; // la particella resta intrappolata....
// definizione dello step.
step =rho/3000.;
// if(step<0.04)step=0.04;
double stepmin=freestep/400.;
if(step<stepmin)step=stepmin;
else if(step>3.)step=3.;
if(step>freestep) step= freestep+0.01;
dse = step*ec*pr->Get_Charg();
path=0.0;
iter_mag=0;
// tracing main loop =======================
while (iter_mag++<MAXITER)
{
if(campo>0) Get_Field(P_Cin); //variabile
/*
evento_.Event,nome,pr->Get_Nome(),pr->pathok,cam_max,rho);
Pvers.print("pvers");
P_Cin.print(pr->Get_Nome());
Field.print("Field")
*/
Puls+= Pvers&Field * dse; // Aggiorno P dopo un cammino ds
Pvers=Puls.NVerso(); // calcola la nuova direzione di volo
P_Cin+= step*Pvers; // aggiorno la posizione
path+=step; // aggiorno il cammino
// check for end iteration...
if(step>freestep) break; // decaduta...o uscita dal mag...
freestep-=step; //aggiorno il cammino da fare
// test if still inside!!Rectangular shape ! for different shape not yet implemented
if (fabs(P_Cin[0])> Lin[0] || fabs(P_Cin[1])> Lin[1] ||fabs(P_Cin[2])> Lin[2]) break;
/*{
if(idev==36)
{
P_Cin.print("X1");Gout<<" Ev "<<evento_.Event<<" iter "<<iter_mag;
}
break;
}*/
if(step>freestep) // se supero con uno step il cammino restante allora...
{
// allora calcolo il cammino potenziale restante...
freestep=Get_CamMag();
if(freestep>vol_max-path)freestep=vol_max-path;
if(step>freestep)
{
step= freestep+0.01;
if(step<0.05)break;
dse = step*ec*pr->Get_Charg();
}
}
}
// tracing end ====================================
// ===> go back to the apparatus system
pr->P.setvn(Dev2Lab(Puls),pr->Get_Massa()); // nota calcola anche la norma
pr->Vers=pr->P.Verso();
X=P_Cin+Centrin;
pr->X=Devc2Lab(X)+ 0.01*pr->Vers; // non norma... + piccolo spostamento
pr->pathok=path;
pr->path_done+=path;
// si deve tener conto del multiploscattering ??
Mscatter(pr,path);
if (iter_mag> 400)
{
Gout<<"\n\n Even="<<evento_.Gen.Event<<" "<<nome
<<" iter "<<iter_mag<<" for "<< pr->Get_Nome()<< " P_norma "<< pr->P.norma
<<" Camx "<< cam_max<<" rho="<<rho<<" step "<<step;
KIK=Puls-PIN;
Devc2Lab(XIN).print("Xin");
pr->X.print("Xout");
PIN.print("Pin");
Puls.print("Pout");
KIK.Norma();
KIK.print("KIK");
// Field.print("Last Field");
Gout<<"\n "<<nome<<" for "<<pr->Get_Nome()<<" Zstp="
<<setprecision(2)<<fixed
<<setw(6)<<freestep<<" Path="
<<setprecision(1)<<fixed
<<setw(7)<<path<<" Iter="
<<setw(5)<<iter_mag<<" end tracing"<<std::endl;
return (1);
}
/*****************************/
if(Debugon==1)
{
KIK=Puls-PIN;
Devc2Lab(XIN).print("Xin");
pr->X.print("Xout");
PIN.print("Pin");
Puls.print("Pout");
KIK.Norma();
KIK.print("KIK");
// Field.print("Last Field");
Gout<<"\n "<<nome<<" for "<<pr->Get_Nome()<<" Zstp="
<<std::setprecision(2)<<fixed
<<setw(6)<<freestep<<" Path="
<<setprecision(1)<<fixed
<<setw(7)<<path<<" Iter="
<<setw(5)<<iter_mag<<" end tracing"<<std::endl;
}
/*****************************/
return (0);
}
