flyopunta/src/src_reaz/reagas.cpp

/***************************************************************************
 *   Copyright (C) 2009 by giuseppe Pierazzini                                      *
 *   giuseppe@pierazzini.it                                                               *
 *                                                                         *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *
 *                                                                         *
 *  Dipartimento di Fisica  E.Fermi / INFN Pisa Italy                      *
 ***************************************************************************/
#include "flyoh.h"

Reagas::Reagas()
{

  Gascls= new Frammenti[1];
  gletti=0;
  rec=0;

  Gout<<"\n  Reagas Done...";
  npp_gas=np_gas=idnow=0;

  if ( BloccoGas>evento_.Mxev ) BloccoGas=evento_.Mxev+1;

  cout<<"\n Gas: Reading block of "<<BloccoGas<< " Gas_ev"<<endl;
  Gout<<"\n Gas: Reading block of "<<BloccoGas<< " Gas_ev"<<endl;
}


Reagas::~Reagas()
{
}
//===============================================
void Reagas::Read()

{

  char line[128];
  string label,simbolo;
  int  npp=0,ntp=0;
  int  nev=0, p_gas=0;
  int Gasskip=0;


//  if ( aperta==0 )

  if ( !Card_gas.is_open() )
    {

      Card_gas.open ( "../../Flukadati/KsuO.txt", ios::in );
      if ( !Card_gas )
        {
          Gout<<"\n Input file out.txt NOT opened ";
          Gout<<"\n no input files ..program stop !!! \n"<<endl;
          exit ( -1 );
        }
      else
        {
          cout<<"\n REAGAS:  Reading data gas  opened.."<<endl;
          Gout<<"\n REAGAS:  Reading data gas  opened.."<<endl;
        }

    }

// qui definisco lo spazio per i dati letti
  delete[] Gascls;

  if ( BloccoGas> ( evento_.Mxev-Eventi_Fatti ) ) BloccoGas= ( evento_.Mxev-Eventi_Fatti )+2;

  Gascls= new Frammenti[BloccoGas];

  p_gas=-1;


  Card_gas.getline ( line, 120 );
  while ( !Card_gas.eof() )
    {

      // Decodifico la linea
      stringstream scheda;
      scheda<<line;
      scheda>>simbolo;

      //  cout<<"\n Gas: "<<p_gas<<"  "<<ntp<<"  "<< line<< " == labl: "<<simbolo<<flush;

      if ( simbolo.size() <1 )
        {
          Card_gas.getline ( line, 120 );
          continue;
        }
      if ( simbolo.size() ==0 );
      else  if ( simbolo.size() != 0 )
        {

          if ( simbolo=="Evento" )
            {
              gletti++;
              if ( gletti>=GasBegin )
                {

                  p_gas++;
//                  if ( p_gas>BloccoGas ) break;

                  scheda>> nev;
                  scheda>> npp;
                  npp-=2;
                  Gascls[p_gas].Np= npp;     // ricavo il numero delle particelle del cluster
                  Gascls[p_gas].Pgas =  new NataInGas[npp]; // Creo il buffer per npp particelle
                  ntp=0;
                }
              else Gasskip++;
            }
          // memorizzo le particelle
// cerco il puntatore standard della particella e memorizzo il momento
          else  if ( ntp<npp )
            {
             
              double px,py,pz,mas;
              scheda>>px>>py>>pz>>mas;
              Gascls[p_gas].Pgas[ntp].mp= Matter->Pntmatter ( simbolo.c_str() );

              Gascls[p_gas].Pgas[ntp].px= px;
              Gascls[p_gas].Pgas[ntp].py= py;
              Gascls[p_gas].Pgas[ntp].pz= pz;
              Gascls[p_gas].Pgas[ntp].massa= abs ( mas );
              ntp++;
            }
        }
      if ( p_gas> ( BloccoGas-2 ) ) break;
      rec++;
      Card_gas.getline ( line, 120 );
    }

  Gout<<"\n REAGAS: Ev "<<Eventi_Fatti<<"  Gas event read  "<<nev<<" new data  blocco "<<p_gas+1<< " skipped "<<Gasskip<<endl;
  cout<<"\n REAGAS: Ev "<<Eventi_Fatti<<"  Gas event read  "<<nev<<" new data  blocco "<<p_gas+1<< " skipped "<<Gasskip<<endl;

  if ( p_gas<BloccoGas-2)
    {
      int evegas_disponibili =Eventi_Fatti+p_gas;
      if ( evento_.Mxev>evegas_disponibili )
        {
          evento_.Mxev=evegas_disponibili; // to avoid crash...
          
          Gout<<"\n REAGAS:  Event gas at disposal  " << evegas_disponibili
          << " Num. maxev  redifined to: "<< evento_.Mxev<<flush;
          cout<<"\n REAGAS:  Event gas at disposal  " << evegas_disponibili
          << " Num. maxev  redifined to: "<< evento_.Mxev<<flush; 
        }
    }
}
//=============================================

Particella * Reagas::Get_GasInt ( int fiot, Particella *pup )
{

  Particella *pr=0;
//    Partname= Dum;


  if ( np_gas==0 ) // first call in the cluster
    {
      npp_gas=  Gascls[fiot].Np;

    }
  else if ( np_gas>=npp_gas )
    {
      np_gas=0;
      return 0;
    }

  Part_db *pmt;
  pmt=Gascls[fiot].Pgas[np_gas].mp;

  if ( pmt==0 )
    {
      np_gas=0;
      return 0;
    }

// trasformo i K0 in Kl/Ks al 50%

  Partname=pmt->Get_Nome();
  if ( Partname==Kap0||Partname==A_Kap0 )
    {
      if ( Pran() >0.5 ) Partname=Kaps;
      else Partname=Kapl;
    }

  pr= SelRea->NuovaPart ( Partname );
  pr->Set_Fato ( Nata,1 );
  pup->cls=pr;
  pr->up=pup;
  pr->cls=pr->dw=pr->lf=pr->rg=0;
  pr->Gx=pr->X=pup->X;       // all born in the same place...
  double gpx=Gascls[fiot].Pgas[np_gas].px;
  double gpy=Gascls[fiot].Pgas[np_gas].py;
  double gpz=Gascls[fiot].Pgas[np_gas].pz;
  double mas=Gascls[fiot].Pgas[np_gas].massa;  // aggiustare la massa ...
  pr->Gp.setvn ( gpx,gpy,gpz,mas );
  pr->P=pr->Gp;
  pr->Vers=pr->P.NVerso();
  pr->Pbar=pr->P;
  np_gas++;
  //pr->Stampa_questa_particella();





  return pr;
}

//==================================================
void Reagas::Print ( int prima,int fine )
{
  Part_db *pmt=0;

  Gout<<"\n  B U F F E R  FLUKA  blocco debug   da "
  << prima <<" a "<<fine <<" dei "<<gletti<<" creati"<<endl;

  for ( int curr=prima; curr<fine;curr++ )
    {
      if ( curr>BloccoGas ) break;
      int npp=Gascls[curr].Np;

      Gout<<"\n  Cls_ev "<< curr <<"  Particelle "<<npp;
      for ( int n=0; n<npp;n++ )
        {
          pmt=Gascls[curr].Pgas[n].mp;
          if ( pmt>0 )
            {
              Gout<<"\n  "<< setw ( 4 ) <<n
              <<setw ( 7 ) <<pmt->Get_Nome()
              <<" "<< setw ( 8 ) <<Gascls[curr].Pgas[n].px
              <<" "<< setw ( 8 ) <<Gascls[curr].Pgas[n].py
              <<" "<< setw ( 8 ) <<Gascls[curr].Pgas[n].pz;
              if ( pmt==ptr_boo ) Gout<<"  particella sconosciuta! ";
            }

        }
      Gout<<"\n fine "<<endl;
    }

}

//=====================================================================
//                    M k C l u s t e r
//====================================================================

void Particella::
MkCluster ()
{
  Particella *pr,*pup;
// qui si leggono i dati prodotti da fluka.
  if ( ClsG%BloccoGas==0 )
    {
      if ( Gasdati==0 ) Gasdati = new Reagas;
      Gasdati->Read();
      cout<<"\n Gasdati letti..."<<std::endl<<std::endl;
//       Gasdati->Print ( 30,40 );
      ClsG=0;
    }
// Fluka();
// qui chiamare fluka???????????????????? o un generatore di particelle
//define the particle cluster...according with Fluka  generation;
//  Aggiorna_link (Particella * pup, Particella * prg, Particella * pdw, Particella * plf)

  pup = this;
  dw=0;
// Generate the read list in particles (level zero)
  while ( pup )
    {
      pr= Gasdati->Get_GasInt ( ClsG, pup );  //Ricavo una particella dalla interazione gas n. ClsG.
      if ( pr==0 ) break;   // fine delle particelle disponibili nel cluster selezionato
//    Gout<<"\n Evento "<<Eventi_Fatti<<" Cls "<<ClsG<<" sele "<<pr->nome;
      pup=pr;
    }
  ClsG++;
  SelRea->Gasint_Done++;


// indice dell'ultima particella letta dal file
  GasIdo=SelRea->Lastpart->Get_Ido();
//ultima particella letta dal file per questa interazione
  Particella *GasUltima=SelRea->Lastpart;
  
// Generate the decay of the  read particles if one ( level one)

  pr=SelRea->LastpartSave->next;
  while ( pr!=0&&pr!=GasUltima )
    {
      pr->Make_decay_chain ();
      pr=pr->next;
    }

}
//=============
void Particella::Make_decay_chain ()
{


  // costruisce la reazione iterativamente generando i decadimenti di
  // ciascuna particella coinvolta nella reazione.

  if ( dw!=0 ) return;  // no decay particles

  int canali = pointer_db->ncan; //punta alla lista dei decadimenti della particella padre [this]
  if ( canali==0 ) return; //no decay!

  Particella *pa,*pb=0,*pc=0;
  double Prob= Pran();
  double Tprob=0.0;

  int icn=0;
  for ( int i=0;i<canali;i++ )
    {
      Tprob+=pointer_db->Canale[i].rate;
      icn=i;
      if ( Prob<Tprob ) break;  // attenzione meglio se Tprob_max>=1. !!!
    }

  if ( pointer_db->Canale[icn].frammenti==2 )
    {

      pa=Decay_in_2 ( pointer_db->Canale[icn].nomea,pointer_db->Canale[icn].nomeb );
      pb=pa->rg;
      pc=0;

    }
  else if ( pointer_db->Canale[icn].frammenti==3 )
    {
      pa = Decay_in_3 ( pointer_db->Canale[icn].nomea, pointer_db->Canale[icn].nomeb,pointer_db->Canale[icn].nomec );
      pb=pa->rg->dw;
      pc=pb->rg;
    }
  else
    {
      Gout<<"\n Make_decay_chain: Error in decay data base partcle  list  !!! "<<nome<<endl;
      exit ( 0 );
    }

  if ( pa!=0 ) pa->Make_decay_chain ();
  if ( pb!=0 ) pb->Make_decay_chain ();
  if ( pc!=0 ) pc->Make_decay_chain ();

}