flyopunta/src/src_dev62/devmagn.cpp

/***************************************************************************
                          devmagn.cpp  -  description
                             -------------------
    begin                : Sat May 26 2001
    copyright            : (C) 2001 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#define ESTERNO  1
#define INTERNO  0
#define INBUCO   -1
#define EPSR     1.0E-3
#define MAXITER   500

#include "parm.h"
#include "part_db.h"
#include "evento.h"
#include "particella.h"
#include "reazione.h"
#include "devmagn.h"

using namespace std;

DevMagn::DevMagn()
{
  devtype=TypDevMagn;
  devclass="NA62 Experimental Magnet";                   
  fun=Magn;
Gout<<"\n\n < "<<nome <<" > Dev Id "<<idev<<" Device type <"<<devclass<<"> typ "<<devtype; 
  //  campo costante  campo=0 default;
  //  campo variabile campo>0; Definire.... ricorda!
  campo=0;
}


DevMagn::~DevMagn()
{}

// in Gvect.h e' definita l'algebra per i  vettori (gvet)) 3d....
// vettori di giuseppe
//===============================================================
//=================DevMagn==============================
//----------------- --------
void DevMagn::Prgeom()
{
  Gout<<"==================  R i v e l a t o r e =================";
  Gout<<"\n --> < "<<nome <<" >  [ "<<fun
  <<" ]  Magnete rettang. con buco rettang.";
  Device::Prgeom();
  Gout<<std::endl;
}
//======================================================
// =======  Magnet tracking

// vettori ==> gvet VW,PL,PT,X,XIN,PIN,P,KIK,RHO,CRHO;

int DevMagn::Traccia(Particella *pr)
{

  // 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);
}
//==================

double  DevMagn::Get_CamMag()
{
  //  Calcolo il cammino fino alla prima superficie di uscita.
  // mi riferisco al centro del buco interno
  // il punto deve essere interno al buco!!!!
  // P_Cin=X-Centrin;    gia' definito prima di chiamare


  Dind=P_Cin-Lin;
  Dins=-(Lin+P_Cin);
  V_dev=Pvers;


 
  return(CamBuco() );
}