flyopunta/src/src_reaz/reaprod.cpp

/***************************************************************************
                          greaprod.cpp  -  description
                             -------------------
    begin                : Sun Nov 26 2000
    copyright            : (C) 2000 by Giuseppe m Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************/
//--------------------------------------------------------------*/
//   Routine   di   B e a m  P r o d u c t i o n                             */
//--------------------------------------------------------------*/
//   Fbeam : beam production           GMP 21/7/92/23.04.99/

#include "flyoh.h"
using namespace GmpConst;
using namespace std;

//---------------------------------------------------
//        M a x S p e t t r o
//------------------------------------------------------
void Reazione::
MaxSpettro()
{
  Particella *pr=avo;
  Reazione* prea=up;
  double dW;
  Gout<<"\n\n   ----------> P r o d u z i o n e <------------  ";
  Gout<<"\n "<<pr->Get_Nome() <<" Beam pruduction probability "
  <<" for the reaction "<<header<<endl;



  if ( inc_pinc<=0.0 )
    { Gout << " \n Incident pinc non defined.... exit "; exit ( 0 );}


  Gout<<"\n Production parameters: \n inc_pinc="
  <<setprecision ( 2 ) <<fixed
  <<setw ( 8 ) <<inc_pinc
  <<" teta="<<setw ( 8 ) <<inc_teta;

  Gout<<"\n Virtual collimator position and dimension...";

  gvet dcoll  = Ccoll-Ctarg; dcoll.Norma();
  Ctarg.print ( "Targ" );
  Ccoll.print ( "Ccol" );Gout<< "  r="<<rcoll;
  dcoll.print ( "dcoll" );
  dW=PiGreco*rcoll*rcoll/ ( dcoll.z*dcoll.z ); // angolo solido coperto dal collimatore virtuale
  Gout<<"\nSolide angle   dW = "<< dW;


  //----------------------------------------------------------------
  while ( prea!=0 )
  {


    // verifico se lo stesso di una  reazione precedente....
    if ( avo->Get_Idm() ==prea->avo->Get_Idm() && Pinc.norma==prea->Pinc.norma )
    {
      maxspettro=prea->maxspettro;
      Dnsudp=prea->Dnsudp;
      prodrate=rate*Dnsudp;
      pincmx=prea->pincmx;

      Gout<<"\n -------";
      Gout<<"\n Maxspettro has found previuos reaction with the same beam!"
      << "   got the same parametrs...";


      Gout<<"\n\n "
      <<setprecision ( 7 ) <<fixed
      <<"  EV/PPP "
      <<setw ( 13 ) << prodrate<<" Max="
      <<setw ( 13 ) <<maxspettro  <<"  at "
      <<setprecision ( 2 ) <<fixed
      <<setw ( 6 ) <<pincmx<<" GeV  for "
      <<setw ( 6 ) <<avo->Get_Nome()
      <<" with "<< rate <<" of B.R.";


      Gout<<"\n ----------------------------------------\n";
      if ( next!=0 ) next->MaxSpettro();
      return;
    }
    prea=prea->up;
  }
  //========================== end  verifica  ========================---

  //loop
// Non !!! Integro sulla sup del collimatore virtuale



  //==========================================
  // integration loop
  //  calcolo il  momentum spectrum...nell'intervallo Pinc +- Wpinc
  //  for particles produced versus the collimator center direction...

  double stepp= wpinc*0.05;   //40 step in pinc
  ppb=Pinc.norma-wpinc;
  double mxpp= Pinc.norma+wpinc+stepp*0.5;
  Dnsudp=0.0;
  Gout<<"\n\n       T A B E L L A ";
  Gout<<"\n     teta         p       dN/dpdw    ";

  // fix parameters for beam particle  production
  // attenzione se l'angolo di produzione è diverso andrebbe corretto...
  vxb=vyb=0.0,vzb=1.0;
  double tett=0.0, Pmedio=0.0;

  while ( ppb<= mxpp )
  {
    /************************* n o t e *******************
    calcolo di wexp!
    la formula di Doble/Gatignon
    dn/do alla targhetta

    */
    if(ubeam==0) prob=Produci(); //Creare una routine utente.....
    else       prob=1000.0;   //prob è la sezione d'urto differenziale!!!!
   
    if ( maxspettro<prob ) { maxspettro=prob;pincmx=ppb;}
    Dnsudp+=prob; //Integro in funzione di p
    Pmedio+=prob*ppb;
    //   qui teta sempre nulla   tra la direzione del p e del  beam generato

    Gout<<"\n "
    <<setprecision ( 4 ) <<fixed
    <<setw ( 9 ) <<tett<<"  "
    <<setw ( 9 ) <<ppb<<"  "
    <<setw ( 9 ) <<prob;

    ppb+=stepp;
  }

  Pmedio= Pmedio/Dnsudp;

  Gout<<"\n\n Doble production max spettrum"<<setw ( 13 ) <<maxspettro
  <<"\n for "<<setw ( 6 ) <<avo->Get_Nome()
  <<"  at "<<setprecision ( 2 ) <<fixed<<setw ( 6 )
  <<pincmx<<" GeV "
  <<"\n with B.R  "<< rate <<" and  P_med "<<Pmedio
  <<setprecision ( 8 ) <<"  in Dw = "<< dW;

  //  Domega by the defining collimator.divided by 16 step of domega.. vedi loop up.
  //

  Dnsudp=Dnsudp*stepp*dW;
  prodrate=rate*Dnsudp;



  Gout<<"\n\n Ev/PPP "
  <<setprecision ( 7 ) <<fixed<<setw ( 13 ) << prodrate;


  Gout<<"\n ----------------------------------------\n";

  if ( next!=0 ) next->MaxSpettro();

}


//--------------------------------------------------------------
//   Routine   di   B e a m  P r o d u c t i o n
//--------------------------------------------------------------
//--------------------------------------------------------------
double Reazione::Produci ()  {  return 1.0;}     //base procedure
//===============================================================

double ReaKld::Produci()


// calcola la probabilita' della produzione Kl di impulso p
// secondo la distribuzione per protone incidente
// Nota:  Kl=Ks=(K0+ K0a)/2 = (K+ + 3K-)/4

{

  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;



  //  old Niels code...end
  //  double A=0.16, B=8.5, C=3.0;
  /*
  double ppbs=pow(ppb,2.);
  double A=0.16, B=10.5, C=3.0;
  double c1= ( B/inc_pinc ) *exp ( -B*ppb/inc_pinc );
  double c2= ( C*ppbs/PiGreco ) *exp ( -C*ppb*ppb*tetq );
  double Kpfluss= A*c1*c2 ;

  // A=0.1, B=13., C=3.5;
  A=0.1, B=13., C=3.5;Inizio

  // spettro*=1. + 8.8*exp ( -.05*ppb );
  // spettro*=exp ( -5.*(ppb-60.)/inc_pinc  );

  c1= ( B/inc_pinc ) *exp ( -B*ppb/inc_pinc );
  c2= ( C*ppbs/PiGreco ) *exp ( -C*ppb*ppb*tetq );
  double Kmfluss= A*c1*c2 ;
  double spettro=(Kpfluss+3.*Kmfluss)*.25;

  */


  //  double C=3.5;


  // italo..

  //double slope =0.032840722 == 1./30.45;   org power ===  1.671;
  double slope =1./30.45;
  double power=1.671;          //org power ===  1.671;
  double spettro =  pow ( ppb,power ) *exp ( -slope*ppb );

  //correzione per cammino....
  double Inizio,Fine;
  Inizio=  Beg_fid - Ctarg.z;
  Fine=    End_fid - Ctarg.z;
  double dzdec=End_fid-Beg_fid;
  double camvol = ctau_avo_su_mass * ppb ;
  double  pathcor= exp ( -Inizio/camvol ) * ( 1.-exp ( -dzdec/camvol ) ) ;
  spettro*=pathcor;

  //  spettro*= ( 0.6325-0.2324e-2*ppb+0.2244e-4*ppb*ppb );  //standard...
  //  spettro*=(0.6325-0.2324e-2*ppb+0.2144e-4*ppb*ppb);


  // spettro*=(0.6325-0.1900e-2*ppb+0.2450e-4*ppb*ppb);  // quasi perfetto
  spettro*= ( 0.6325-0.1900e-2*ppb+0.27e-4*ppb*ppb );  // forse più perfetto
  return spettro;
}
double ReaKsd::Produci()

// ok for Ksdec. Vedi commento sopra.

{
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double A=0.12, B=16.2, C=3.2;
  double c1= ( B/inc_pinc ) *exp ( -B*ppb/inc_pinc );
  double c2= ( C*ppb*ppb/PiGreco ) *exp ( -C*ppb*ppb*tetq );
  double Kpfluss= A*c1*c2 ;
  double spettro=Kpfluss;

  return spettro;
}

//==========================
double ReaKpd::Produci()

// calcola la probabilita' della produzione Kp di impulso p
// secondo la distribuzione per protone incidente
// dNk/Np=(A/inc_pinc)*p^3*exp(-B*p/inc_pinc -C*(p*teta)^2 )[dp/p*dOmega]

{
  double A=0.16, B=8.5, C=3.0;

  //  vedi articolo
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double c1= ( B/inc_pinc ) *exp ( -B*ppb/inc_pinc );
  double c2= ( C*ppb*ppb/PiGreco ) *exp ( -C*ppb*ppb*tetq );
  return ( A*c1*c2 );



}
//------------------
double ReaKmd::Produci()
// calcola la probabilita' della produzione Km di impulso p
// secondo la distribuzione per protone incidente
// dNk/Np=(A/inc_pinc)*p^3*exp(-B*p/inc_pinc -C*(p*teta)^2 )[dp/p*dOmega]
// con  A=0.58/inc_pinc, B=13.0, C=3.5, inc_pinc=400. Gev e teta=-2 10-3 rad.

{
  double A=.1, B=13.0, C=3.5;
  //  vedi articolo
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double c1= ( B/inc_pinc ) *exp ( -B*ppb/inc_pinc );
  double c2= ( C*ppb*ppb/PiGreco ) *exp ( -C*ppb*ppb*tetq );
  return ( A*c1*c2 );
}
//------------------
double Reappd::Produci()
// calcola la probabilita' della produzione pip di impulso p

//............................
{
  double A=.8, B=11.5, C=5.0;
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double c1= ( B/inc_pinc ) *exp ( -B*ppb/inc_pinc );
  double c2= ( C*ppb*ppb/PiGreco ) *exp ( -C*ppb*ppb*tetq );
  return ( A*c1*c2 );

}
//------------------------------
double Reapmd::Produci()

// calcola la probabilita' della produzione pim di impulso p
// con  A=5.86/inc_pinc, B=11.5, C=5.0, inc_pinc=400. Gev e teta=2 10-3 mrad.
//............................
{
  double A=5.86, B=11.5, C=5.0;
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double a= - B*ppb/inc_pinc;
  double b= - C*ppb*ppb*tetq;
  return ( A* ppb*ppb*exp ( a+b ) /inc_pinc );
}
//------------------------------

double ReaEtad::Produci()

// calcola la probabilita' della produzione eta

//............................
{
  double  mmq = 0.5473* 0.5473;
  double ener = sqrt ( ppb*ppb+mmq );
  double teta=acos ( vzb );
  if ( teta>0.02||teta<0.00001 ) return ( 0.0 );
  double w1, w2;
  w1=-46.309+2.6560*ener-.36088e-1*ener*ener+.18706e-3*
     ener*ener*ener-.33651e-6*ener*ener*ener*ener;
  w2=teta*exp ( -36646.*teta*teta );

  return ( w1*w2 );

}
//------------------------------
double ReaPi0d::Produci()


// come per il pim ??????????
//............................
{
  double A=5.86, B=11.5, C=5.0;
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double a= - B*ppb/inc_pinc;
  double b= - C*ppb*ppb*tetq;
  return ( A* ppb*ppb*exp ( a+b ) /inc_pinc );
}
//---------------------------------------
double ReaEleb::Produci()

// calcola la probabilita' della produzione elec/elep
//......... spettro come vole iacopini...................
{
  double A=40.,B=50.,  C=5.;
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double a=2.- ( ppb-B ) *.05;if ( a<0. ) a=0.; // spettro lineare
  double b= - C*ppb*ppb*tetq;
  return ( A*a*exp ( b ) /inc_pinc );
  //    return 1. ;
}
//------------------------------
double ReaPb::Produci()
// calcola la probabilita' della produzione proton di impulso p
// con  A=.14, C=3.5, inc_pinc=400. Gev e teta=2 10-3 m
//............................
{
  double A=0.14, C=3.5;
  gvet vtet= Vbeam&Vinc_ppp; vtet.Norma();
  double  tetq= asin ( vtet.norma );tetq*=tetq;
  double a=0.4*log ( ppb/inc_pinc );
  double b= - C*ppb*ppb*tetq;
  return ( A* ppb*exp ( a+b ) );
}