flyopunta/src/src_reaz/reabeam.cpp

/***************************************************************************
 ***
 *                                                                         *
 *   NA48  simulation program.                                             *
 *                                                                         *
 ***************************************************************************/

//
// beam routine   gmp Tirrenia 29-10-99  da...
// Tirrenia  28.08.00

#include "flyoh.h"
using namespace GmpConst;


//--------------------------------------------------------------*/
//          G e t b e a m                                       */
//--------------------------------------------------------------*/
// get the beam ..... versione base   10.02.00   gmp
// per le classe Reazione,Reapmd..............;
// va bene per i decadimenti neutri e carichi... attenzione i carichi
//

int Reazione::
Get_Beam ()
{
  //--------   genration loop
  //metodo
  // loop per l'estrazione di un impulso che soddisfi
  // all'ottica del fascio....  si verifica con
  // la distribuzione di Double.

  int iter=0;
  double prb,fi,ro;
  double dwpb = 2.*wpinc;
  double volo=0.0;
  
  double pathd=0.0;

  // punto targhetta...
  ro=rtarg * sqrt ( Pran() );
  fi= Pran() * DuePiGreco;
  Xtarg.setv ( ro*cos ( fi ), ro*sin ( fi ),dztarg* ( Pran() - 0.5 ) );
  Xtarg+=Ctarg;
        


  // calcolo la direzione di volo ricordando che deve passare per il
  // defining collimetor
  // ro= rcoll *pow(Pran(),.5);
  ro= rcoll *sqrt ( Pran() );
  fi= Pran() *DuePiGreco;
  //  User beam parameters;
  //  ubeam=0:si sceglie il mom nell'intervallo fissato con distribuzione uniformewmente
  //  ubeam=1 si sceglie  con distribuzione gaussiana.....e forma ellittica
  if ( ubeam==0 )
  {
    Vbeam.setv ( ro*cos ( fi ), ro*sin ( fi ), 0.0 );  // coordinate al collimatore virtuale di def.

  }
  else
  {
    // coordinate al coll virt. con distr. ellitica  (quando??)
 //   Vbeam.setv ( ro*cos ( fi ), 0.8*ro*sin ( fi ), 0.0 );
    Vbeam.setv ( ro*cos ( fi ), ro*sin ( fi ), 0.0 );

  }

  Vbeam-=Xtarg-Ccoll;
  Vbeam.Norma();       // calcola la norma
  Vbeam=Vbeam.Verso(); // versore del beam..
  Vbeam.putv( vxb,vyb,vzb );

// selezione del momento

 if ( ubeam==0 )
  {
  while ( 1 )  // se la reazione dipende da teta spostare questo while in alto...
  {
    iter++;
      ppb = Pinc.norma + dwpb* ( Pran()-0.5 ); // dwpb intervallo in impulso...
      prb=Produci() /maxspettro;
      double rrr=Pran();
      if ( rrr<prb ) break;

      if ( iter>1000 )
      {
        Gout<<"\nEv "<<evento_.Gen.Event<<" "<<avo->Get_Nome() <<" in rea_Id "<<id<<" Beam_error iter="
        <<iter<<" Prb="<< prb<<" Pmx="<<maxspettro;
        break;  // esco per forza
      }
    }
}
    else ppb = Pgauss ( wpinc, Pinc.norma );  // beam a forma gaussiana(gvet_0)

 

  pathd=Zona;

  if ( ctau_avo_su_mass>0.0 )
  {
    // la particella ha una massa ed una vita media
    volo = ctau_avo_su_mass * ppb;

                while(1)
                {
                        evento_.Gen.tempo=TBurst->Get_ProdTempo();
                        
                        pathd= - volo*log ( Pran() )- Xtarg.z;
                        if(Beg_fid<=pathd&&pathd<=End_fid)break;
                
                
                }
                                        
  /*      
      // gli eventi sono generati solo nella regione di fiducia
      // per rendere il programma piu' efficiente, ma attenzione ai fondi....e normalizzazioni varie
      // conviene allora prendere una regione non limitata superiormente...
   */
    
    
    
  }
  else evento_.Gen.tempo=TBurst->Get_ProdTempo(); // se non decade.....

  evento_.Gen.tdif=TBurst->Get_Tdif();  // differenza tra due eventi accettati.....             
  evento_.Gen.burst=TBurst->Get_Nburst();  // differenza tra due eventi accettati......

  trasmessi++;

  //-----------------------------------------------------------------
  avo->X = Xtarg;
  Pbeam= Vbeam*ppb;
  avo->P.setvn ( Pbeam,avo->Get_Massa() ); // P quadrivettore impulso energia
  avo->Vers=Vbeam;
  avo->path_max=pathd;
  avo->X2Gx();
  avo->P2Gp();
  avo->Set_Fato ( Nata,-1 );
  avo->Set_Nata();  //nata=1;

  return ( 1 );
}