flyopunta/src/src_dev/devflt.cpp

/***************************************************************************
                          devflt  -  description
                             -------------------
    begin                : Sun Aug 15 2004
    copyright            : (C) 2004 by gmp
    email                : peppe@newpeppe
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/

#define ESTERNO  1
#define INTERNO  0
#define INBUCO   -1
#define EPSR     1.0E-3
#define MAXITER   500

#include "parm.h"
#include "evento.h"
#include "particella.h"
#include "reazione.h"
#include "devflt.h"

using namespace std;

DevFlt::DevFlt()
 : DevRtCl()
{
   devtype=TypDevFlt;
   devclass="Magn_Filter";
//  campo costante  campo=0 default;
//  campo variabile campo>0; Definire.... ricorda! 
   campo=0;
}


DevFlt::~DevFlt()
{
}

//----------------- --------
void DevFlt::Prgeom()
 {

  Gout<<"==================  R i v e l a t o r e =================";
  Gout<<"\n --> < "<<nome <<" >  [ "<<fun
      <<" ]  Mag in FE con buco cilindrico senza campo."<<std::endl;
  Device::Prgeom();
  Gout<<"\n";
}
//======================================================
// =======  Magnet tracking

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

int DevFlt::Traccia(Particella *pr)
 {
   
// see the DevMagn:: Traccia()
 // 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) )

  const double ec=2.99792458e-7;
//****

  
// qui il campo è differente da zero solo nel ferro... il campo nel buco è nullo!!!  
// si vuole filtrare via i muoni....
  
  if(nowpos==ESTERNO || nowpos==INBUCO)   // se esterno o nel buco traccia regolarmente...
      { Device::Traccia(pr);  return  0;}
// allora e' INTERNO: si traccia in campo magnetico.      
  if (pr->P.norma == 0.0)     return  0;    // no motion ....
  if (pr->Get_Charg() == 0.0)               // neutra traccio regolamente
      { Device::Traccia(pr);  return  0;}

 
//------------
  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
   P_Cin=X;
  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(rho<2.*freestep) return 1;    // la particella resta intrappolata....

// definizione dello step.
  step =rho/3000.;
  if(step<2.)step=2.;
  else if(step>25.)step=25.;
  if(step>freestep) step= freestep+0.04;
  dse = step*ec*pr->Get_Charg();
  path=0.0;
  iter_mag=0;
  
  
  
     
 if(Debugon==1)
  {// ok nel caso di campo costante ==> valore significativo
    double forkik=ec*Field.norma*cam_max*pr->Get_Charg();
    Gout<<"\n\n Ev="<<evento_.Gen.Event<<" "<<pr->Get_Nome()<<" in  "<<nome
        <<"  pathok="<<pr->pathok<<" Cammx="<< cam_max
        <<"  freestep "<<freestep        
        <<" Kik="<<forkik<<"  rho="<<rho;
    Field.print("Field");
    
  }
  
     
 
  
// 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

      if (fabs(X[0])> Lout[0] || fabs(X[1])> Lout[1] ||fabs(X[2])> Lout[2])   break;
      
//      if (fabs(P_Cin[0])> Loin[0] || fabs(P_Cin[1])> Lin[1] ||fabs(P_Cin[2])> Lin[2])   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.04;
         if(step<0.06) break;
         dse = step*ec*pr->Get_Charg();
        }
      }

      
   }
    
  

// tracing end ====================================
 
  if (iter_mag> 499)
     {
      Gout<<"\n\n Ev="<<evento_.Gen.Event<<" Trac. in "<<nome
          <<" iter "<<iter_mag<<" for "<<  pr->Get_Nome()<< " P norma "<< pr->P.norma
          <<" camx "<< cam_max<<"  rho="<<rho<<" step "<<step;    
      return (1);
     }

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

/*****************************/
 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()<<"  Patok="<<freestep<<" Path "
         <<path<<" Iter="<< iter_mag<<std::endl;
    }
/*****************************/


  return (0);
}
//==================

double  DevFlt::Get_CamMag()
{
//  Calcolo il cammino fino alla prima superficie di uscita.
// dall'interno del dev 
// (In particolare qui il punto del mu  deve essere interno nel ferro!!!!
// P_Cin=X-Centrin;    gia' definito prima di chiamare
 Dind=P_Cin-Lin;
 Dins=-(Lin+P_Cin);
 V_dev=Pvers;
 return(CamInter() );
}