[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//
//  Generator for slow nucleons in pA interactions. 
//  Source is modelled by a relativistic Maxwell distributions.
//  This class cooparates with AliCollisionGeometry if used inside AliGenCocktail.
//  In this case the number of slow nucleons is determined from the number of wounded nuclei
//  using a realisation of AliSlowNucleonModel.
//  Original code by  Ferenc Sikler  <sikler@rmki.kfki.hu>
// 

#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TParticle.h>

#include "AliConst.h"
#include "AliCollisionGeometry.h"
#include "AliStack.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliGenSlowNucleons.h"
#include "AliSlowNucleonModel.h"

ClassImp(AliGenSlowNucleons)

    
AliGenSlowNucleons::AliGenSlowNucleons()
    :AliGenerator(-1),
     fCMS(0.),
     fMomentum(0.),
     fBeta(0.),
     fPmax (0.),
     fCharge(0),
     fProtonDirection(1.),
     fTemperatureG(0.), 
     fBetaSourceG(0.),
     fTemperatureB(0.),
     fBetaSourceB(0.),
     fNgp(0),
     fNgn(0),
     fNbp(0),
     fNbn(0),
     fDebug(0),
     fDebugHist1(0),
     fDebugHist2(0),
     fThetaDistribution(),
     fCosThetaGrayHist(),
     fCosTheta(),
     fBeamCrossingAngle(0.),
     fBeamDivergence(0.),
     fBeamDivEvent(0.),
     fSmearMode(2),
     fSlowNucleonModel(0)
{
// Default constructor
    fCollisionGeometry = 0;
}

AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
    :AliGenerator(npart),
     fCMS(14000.),
     fMomentum(0.),
     fBeta(0.),
     fPmax (10.),
     fCharge(1),
     fProtonDirection(1.),
     fTemperatureG(0.05), 
     fBetaSourceG(0.05),
     fTemperatureB(0.005),
     fBetaSourceB(0.),
     fNgp(0),
     fNgn(0),
     fNbp(0),
     fNbn(0),
     fDebug(0),
     fDebugHist1(0),
     fDebugHist2(0),
     fThetaDistribution(),
     fCosThetaGrayHist(),
     fCosTheta(),
     fBeamCrossingAngle(0.),
     fBeamDivergence(0.),
     fBeamDivEvent(0.),
     fSmearMode(2),
     fSlowNucleonModel(new AliSlowNucleonModel())

{
// Constructor
    fName  = "SlowNucleons";
    fTitle = "Generator for gray particles in pA collisions";
    fCollisionGeometry = 0;
}

//____________________________________________________________
AliGenSlowNucleons::~AliGenSlowNucleons()
{
// Destructor
    delete  fSlowNucleonModel;
}

void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
// Set direction of the proton to change between pA (1) and Ap (-1)
  fProtonDirection = dir / TMath::Abs(dir);
}

void AliGenSlowNucleons::Init()
{
  //
  // Initialization
  //
    Double_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
    fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
    fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
    //printf("  fMomentum %f    fBeta %1.10f\n",fMomentum, fBeta);
    if (fDebug) {
	fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
	fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
    	fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
    }
    //
    // non-uniform cos(theta) distribution
    //
    if(fThetaDistribution != 0) {
	fCosTheta = new TF1("fCosTheta",
			    "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
			    -1., 1.);
    }
   
    printf("\n  AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
}

void AliGenSlowNucleons::FinishRun()
{
// End of run action
// Show histogram for debugging if requested.
    if (fDebug) {
	TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
	c->Divide(2,1);
	c->cd(1);
	fDebugHist1->Draw("colz");
	c->cd(2);
	fDebugHist2->Draw();
	c->cd(3);
	fCosThetaGrayHist->Draw();
    }
}


void AliGenSlowNucleons::Generate()
{
  //
  // Generate one event
  //
  //
  // Communication with Gray Particle Model 
  // 
    if (fCollisionGeometry) {
	Float_t b   = fCollisionGeometry->ImpactParameter();
	//	Int_t  nn   = fCollisionGeometry->NN();
	//      Int_t  nwn  = fCollisionGeometry->NwN();
	//      Int_t  nnw  = fCollisionGeometry->NNw();
	//      Int_t  nwnw = fCollisionGeometry->NwNw();
	
	// (1) Sikler' model 
	if(fSmearMode==0) fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	// (2) Model inspired on exp. data at lower energy (Gallio-Oppedisano)
	// --- smearing the Ncoll fron generator used as input 
	else if(fSmearMode==1) fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	// --- smearing directly Nslow 
	else if(fSmearMode==2) fSlowNucleonModel->GetNumberOfSlowNucleons2s(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	if (fDebug) {
	    //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
	    printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
	    fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
	    fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);

	}
    }     

   //
    Float_t p[3] = {0., 0., 0.}, theta=0;
    Float_t origin[3] = {0., 0., 0.};
    Float_t time = 0.;
    Float_t polar [3] = {0., 0., 0.};    
    Int_t nt, i, j;
    Int_t kf;
     
    // Extracting 1 value per event for the divergence angle
    Double_t rvec = gRandom->Gaus(0.0, 1.0);
    fBeamDivEvent = fBeamDivergence * TMath::Abs(rvec);
    printf("\n  AliGenSlowNucleons: applying beam divergence %f mrad to slow nucleons\n",fBeamDivEvent*1000.);

    if(fVertexSmear == kPerEvent) {
	Vertex();
	for (j=0; j < 3; j++) origin[j] = fVertex[j];
	time = fTime;
    } // if kPerEvent
//
//  Gray protons
//
    fCharge = 1;
    kf = kProton;    
    for(i = 0; i < fNgp; i++) {
	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-2);
	KeepTrack(nt);
	SetProcessID(nt,kGrayProcess);
    }
//
//  Gray neutrons
//
    fCharge = 0;
    kf = kNeutron;    
    for(i = 0; i < fNgn; i++) {
	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-2);
	KeepTrack(nt);
	SetProcessID(nt,kGrayProcess);
    }
//
//  Black protons
//
    fCharge = 1;
    kf = kProton;    
    for(i = 0; i < fNbp; i++) {
	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-1);
	KeepTrack(nt);
	SetProcessID(nt,kBlackProcess);
    }
//
//  Black neutrons
//
    fCharge = 0;
    kf = kNeutron;    
    for(i = 0; i < fNbn; i++) {
	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
	PushTrack(fTrackIt, -1, kf, p, origin, polar,
		 time, kPNoProcess, nt, 1.,-1);
	KeepTrack(nt);
	SetProcessID(nt,kBlackProcess);
    }
}


void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, 
	Double_t beta, Float_t* q, Float_t &theta)

{
/* 
   Emit a slow nucleon with "temperature" T [GeV], 
   from a source moving with velocity beta         
   Three-momentum [GeV/c] is given back in q[3]    
*/

 //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
 
 Double_t m, pmax, p, f, phi;
 TDatabasePDG * pdg = TDatabasePDG::Instance();
 const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 
 /* Select nucleon type */
 if(charge == 0) m = kMassNeutron;
 else m = kMassProton;

 /* Momentum at maximum of Maxwell-distribution */

 pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));

 /* Try until proper momentum                                  */
 /* for lack of primitive function of the Maxwell-distribution */
 /* a brute force trial-accept loop, normalized at pmax        */

 do
 {
     p = Rndm() * fPmax;
     f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 }
 while(f < Rndm());

 /* Spherical symmetric emission for black particles (beta=0)*/
 if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 else if(fThetaDistribution!=0){
   Double_t costheta = fCosTheta->GetRandom();
   theta = TMath::ACos(costheta);
 }
 //
 phi   = 2. * TMath::Pi() * Rndm();


 /* Determine momentum components in system of the moving source */
 q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 q[2] = p * TMath::Cos(theta);
 //if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]);


 /* Transform to system of the target nucleus                             */
 /* beta is passed as negative, because the gray nucleons are slowed down */
 Lorentz(m, -beta, q);
 //if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]);

 /* Transform to laboratory system */
 Lorentz(m, fBeta, q);
 q[2] *= fProtonDirection; 
 if(fDebug==1)printf("\n Momentum after LHC boost: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
 
 if(fBeamCrossingAngle>0.) BeamCrossDivergence(1, q); // applying crossing angle
 if(fBeamDivergence>0.) BeamCrossDivergence(2, q);    // applying divergence
 
}

Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
{
/* Relativistic Maxwell-distribution */
    Double_t ekin;
    ekin = TMath::Sqrt(p*p+m*m)-m;
    return (p*p * exp(-ekin/T));
}


//_____________________________________________________________________________
void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
{
/* Lorentz transform in the direction of q[2] */
 
    Double_t gamma  = 1./TMath::Sqrt((1.-beta)*(1.+beta)); 
    Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
    q[2] = gamma * (q[2] + beta*energy);
    //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
}

//_____________________________________________________________________________
void AliGenSlowNucleons::BeamCrossDivergence(Int_t iwhat, Float_t *pLab)
{
  // Applying beam divergence and crossing angle
  //
  Double_t pmod = TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]+pLab[2]*pLab[2]);

  Double_t tetdiv = 0.;
  Double_t fidiv = 0.;
  if(iwhat==1){
    tetdiv = fBeamCrossingAngle;
    fidiv = k2PI/4.;
  }
  else if(iwhat==2){
    tetdiv = fBeamDivEvent;
    fidiv = (gRandom->Rndm())*k2PI;
  }

  Double_t tetpart = TMath::ATan2(TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]), pLab[2]);
  Double_t fipart=0.;
  if(TMath::Abs(pLab[1])>0. || TMath::Abs(pLab[0])>0.) fipart = TMath::ATan2(pLab[1], pLab[0]);
  if(fipart<0.) {fipart = fipart+k2PI;}
  tetdiv = tetdiv*kRaddeg;
  fidiv = fidiv*kRaddeg;
  tetpart = tetpart*kRaddeg;
  fipart = fipart*kRaddeg;
  
  Double_t angleSum[2]={0., 0.};
  AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum);
  
  Double_t tetsum = angleSum[0];
  Double_t fisum  = angleSum[1];
  //printf("tetpart %f fipart %f tetdiv %f fidiv %f angleSum %f %f\n",tetpart,fipart,tetdiv,fidiv,angleSum[0],angleSum[1]);
  tetsum = tetsum*kDegrad;
  fisum = fisum*kDegrad;
  
  pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum);
  pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum);
  pLab[2] = pmod*TMath::Cos(tetsum);
  if(fDebug==1){
    if(iwhat==1) printf(" Beam crossing angle %f mrad ", fBeamCrossingAngle*1000.);
    else if(iwhat==2) printf(" Beam divergence %f mrad ", fBeamDivEvent*1000.);
    printf("  p = (%f, %f, %f)\n",pLab[0],pLab[1],pLab[2]);
  }
}
  
//_____________________________________________________________________________
void  AliGenSlowNucleons::AddAngle(Double_t theta1, Double_t phi1, 
	Double_t theta2, Double_t phi2, Double_t *angleSum)
{ 
  // Calculating the sum of 2 angles  
  Double_t temp = -1.;
  Double_t conv = 180./TMath::ACos(temp);
  
  Double_t ct1 = TMath::Cos(theta1/conv);
  Double_t st1 = TMath::Sin(theta1/conv);
  Double_t cp1 = TMath::Cos(phi1/conv);
  Double_t sp1 = TMath::Sin(phi1/conv);
  Double_t ct2 = TMath::Cos(theta2/conv);
  Double_t st2 = TMath::Sin(theta2/conv);
  Double_t cp2 = TMath::Cos(phi2/conv);
  Double_t sp2 = TMath::Sin(phi2/conv);
  Double_t cx = ct1*cp1*st2*cp2+st1*cp1*ct2-sp1*st2*sp2;
  Double_t cy = ct1*sp1*st2*cp2+st1*sp1*ct2+cp1*st2*sp2;
  Double_t cz = ct1*ct2-st1*st2*cp2;
  
  Double_t rtetsum = TMath::ACos(cz);
  Double_t tetsum = conv*rtetsum;
  if(TMath::Abs(tetsum)<1e-4 || tetsum==180.) return;

  temp = cx/TMath::Sin(rtetsum);
  if(temp>1.) temp=1.;
  if(temp<-1.) temp=-1.;
  Double_t fisum = conv*TMath::ACos(temp);
  if(cy<0) {fisum = 360.-fisum;}
  
  angleSum[0] = tetsum;
  angleSum[1] = fisum;
}  

//_____________________________________________________________________________
void AliGenSlowNucleons::SetProcessID(Int_t nt, UInt_t process)
{
  // Tag the particle as
  // gray or black
  if (fStack)
    fStack->Particle(nt)->SetUniqueID(process);
  else 
    gAlice->GetMCApp()->Particle(nt)->SetUniqueID(process);
}
Commit	Line	Data
c9a8628a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
88cb7938	16	/* $Id$ */
c9a8628a	17
f687abfe	18	//
	19	// Generator for slow nucleons in pA interactions.
	20	// Source is modelled by a relativistic Maxwell distributions.
	21	// This class cooparates with AliCollisionGeometry if used inside AliGenCocktail.
	22	// In this case the number of slow nucleons is determined from the number of wounded nuclei
	23	// using a realisation of AliSlowNucleonModel.
	24	// Original code by Ferenc Sikler <sikler@rmki.kfki.hu>
	25	//
c9a8628a	26
	27	#include <TDatabasePDG.h>
	28	#include <TPDGCode.h>
	29	#include <TH2F.h>
5a9b7f8e	30	#include <TH1F.h>
5a9b7f8e	31	#include <TF1.h>
c8f7f6f9	32	#include <TCanvas.h>
6218004b	33	#include <TParticle.h>
c9a8628a	34
74cfba91	35	#include "AliConst.h"
c9a8628a	36	#include "AliCollisionGeometry.h"
6218004b	37	#include "AliStack.h"
87f95355	38	#include "AliRun.h"
87f95355	39	#include "AliMC.h"
c9a8628a	40	#include "AliGenSlowNucleons.h"
	41	#include "AliSlowNucleonModel.h"
	42
706938e6	43	ClassImp(AliGenSlowNucleons)
1c56e311	44
c9a8628a	45
1c56e311	46	AliGenSlowNucleons::AliGenSlowNucleons()
	47	:AliGenerator(-1),
	48	fCMS(0.),
	49	fMomentum(0.),
	50	fBeta(0.),
	51	fPmax (0.),
1c56e311	52	fCharge(0),
61cf877f	53	fProtonDirection(1.),
1c56e311	54	fTemperatureG(0.),
	55	fBetaSourceG(0.),
	56	fTemperatureB(0.),
	57	fBetaSourceB(0.),
	58	fNgp(0),
	59	fNgn(0),
	60	fNbp(0),
	61	fNbn(0),
	62	fDebug(0),
	63	fDebugHist1(0),
	64	fDebugHist2(0),
	65	fThetaDistribution(),
	66	fCosThetaGrayHist(),
	67	fCosTheta(),
74cfba91	68	fBeamCrossingAngle(0.),
	69	fBeamDivergence(0.),
	70	fBeamDivEvent(0.),
c3e58d2c	71	fSmearMode(2),
1c56e311	72	fSlowNucleonModel(0)
c9a8628a	73	{
c9a8628a	74	// Default constructor
c9a8628a	75	fCollisionGeometry = 0;
	76	}
	77
	78	AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
1c56e311	79	:AliGenerator(npart),
	80	fCMS(14000.),
	81	fMomentum(0.),
	82	fBeta(0.),
	83	fPmax (10.),
1c56e311	84	fCharge(1),
61cf877f	85	fProtonDirection(1.),
230d85c6	86	fTemperatureG(0.05),
1c56e311	87	fBetaSourceG(0.05),
230d85c6	88	fTemperatureB(0.005),
1c56e311	89	fBetaSourceB(0.),
	90	fNgp(0),
	91	fNgn(0),
	92	fNbp(0),
	93	fNbn(0),
	94	fDebug(0),
	95	fDebugHist1(0),
	96	fDebugHist2(0),
	97	fThetaDistribution(),
	98	fCosThetaGrayHist(),
	99	fCosTheta(),
74cfba91	100	fBeamCrossingAngle(0.),
	101	fBeamDivergence(0.),
	102	fBeamDivEvent(0.),
c3e58d2c	103	fSmearMode(2),
1c56e311	104	fSlowNucleonModel(new AliSlowNucleonModel())
74cfba91	105
c9a8628a	106	{
	107	// Constructor
	108	fName = "SlowNucleons";
	109	fTitle = "Generator for gray particles in pA collisions";
c9a8628a	110	fCollisionGeometry = 0;
c9a8628a	111	}
	112
	113	//____________________________________________________________
	114	AliGenSlowNucleons::~AliGenSlowNucleons()
	115	{
	116	// Destructor
	117	delete fSlowNucleonModel;
	118	}
	119
83b3cd6c	120	void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
	121	// Set direction of the proton to change between pA (1) and Ap (-1)
	122	fProtonDirection = dir / TMath::Abs(dir);
	123	}
c9a8628a	124
	125	void AliGenSlowNucleons::Init()
	126	{
	127	//
	128	// Initialization
	129	//
88a9f852	130	Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
271e6aae	131	fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
c9a8628a	132	fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
88a9f852	133	//printf(" fMomentum %f fBeta %1.10f\n",fMomentum, fBeta);
c9a8628a	134	if (fDebug) {
c8f7f6f9	135	fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
c8f7f6f9	136	fDebugHist2 = new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.);
5a9b7f8e	137	fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
	138	}
	139	//
	140	// non-uniform cos(theta) distribution
	141	//
	142	if(fThetaDistribution != 0) {
	143	fCosTheta = new TF1("fCosTheta",
	144	"(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))exp(2.x/3.14159265358979312)",
	145	-1., 1.);
c9a8628a	146	}
74cfba91	147
74cfba91	148	printf("\n AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
c9a8628a	149	}
	150
	151	void AliGenSlowNucleons::FinishRun()
	152	{
f687abfe	153	// End of run action
f687abfe	154	// Show histogram for debugging if requested.
c9a8628a	155	if (fDebug) {
c8f7f6f9	156	TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
	157	c->Divide(2,1);
	158	c->cd(1);
58776b75	159	fDebugHist1->Draw("colz");
c8f7f6f9	160	c->cd(2);
c8f7f6f9	161	fDebugHist2->Draw();
5a9b7f8e	162	c->cd(3);
5a9b7f8e	163	fCosThetaGrayHist->Draw();
c9a8628a	164	}
	165	}
	166
	167
	168	void AliGenSlowNucleons::Generate()
	169	{
	170	//
	171	// Generate one event
	172	//
	173	//
	174	// Communication with Gray Particle Model
	175	//
1c22a81e	176	if (fCollisionGeometry) {
1c22a81e	177	Float_t b = fCollisionGeometry->ImpactParameter();
230d85c6	178	// Int_t nn = fCollisionGeometry->NN();
	179	// Int_t nwn = fCollisionGeometry->NwN();
	180	// Int_t nnw = fCollisionGeometry->NNw();
	181	// Int_t nwnw = fCollisionGeometry->NwNw();
88a9f852	182
c3e58d2c	183	// (1) Sikler' model
	184	if(fSmearMode==0) fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	185	// (2) Model inspired on exp. data at lower energy (Gallio-Oppedisano)
	186	// --- smearing the Ncoll fron generator used as input
	187	else if(fSmearMode==1) fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	188	// --- smearing directly Nslow
	189	else if(fSmearMode==2) fSlowNucleonModel->GetNumberOfSlowNucleons2s(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
1c22a81e	190	if (fDebug) {
88a9f852	191	//printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
230d85c6	192	printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, fNgn, fNbp, fNbn);
58776b75	193	fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
c8f7f6f9	194	fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
74cfba91	195
1c22a81e	196	}
1c22a81e	197	}
c9a8628a	198
c9a8628a	199	//
88a9f852	200	Float_t p[3] = {0., 0., 0.}, theta=0;
c9a8628a	201	Float_t origin[3] = {0., 0., 0.};
21391258	202	Float_t time = 0.;
c9a8628a	203	Float_t polar [3] = {0., 0., 0.};
99f5114f	204	Int_t nt, i, j;
c9a8628a	205	Int_t kf;
74cfba91	206
	207	// Extracting 1 value per event for the divergence angle
	208	Double_t rvec = gRandom->Gaus(0.0, 1.0);
	209	fBeamDivEvent = fBeamDivergence * TMath::Abs(rvec);
	210	printf("\n AliGenSlowNucleons: applying beam divergence %f mrad to slow nucleons\n",fBeamDivEvent*1000.);
	211
99f5114f	212	if(fVertexSmear == kPerEvent) {
	213	Vertex();
	214	for (j=0; j < 3; j++) origin[j] = fVertex[j];
21391258	215	time = fTime;
99f5114f	216	} // if kPerEvent
c9a8628a	217	//
	218	// Gray protons
	219	//
	220	fCharge = 1;
	221	kf = kProton;
1c22a81e	222	for(i = 0; i < fNgp; i++) {
5a9b7f8e	223	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
5a9b7f8e	224	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf	225	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	226	time, kPNoProcess, nt, 1.,-2);
c9a8628a	227	KeepTrack(nt);
87f95355	228	SetProcessID(nt,kGrayProcess);
c9a8628a	229	}
	230	//
	231	// Gray neutrons
	232	//
	233	fCharge = 0;
	234	kf = kNeutron;
1c22a81e	235	for(i = 0; i < fNgn; i++) {
5a9b7f8e	236	GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
5a9b7f8e	237	if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
642f15cf	238	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	239	time, kPNoProcess, nt, 1.,-2);
c9a8628a	240	KeepTrack(nt);
87f95355	241	SetProcessID(nt,kGrayProcess);
c9a8628a	242	}
	243	//
	244	// Black protons
	245	//
	246	fCharge = 1;
	247	kf = kProton;
1c22a81e	248	for(i = 0; i < fNbp; i++) {
5a9b7f8e	249	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf	250	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	251	time, kPNoProcess, nt, 1.,-1);
c9a8628a	252	KeepTrack(nt);
87f95355	253	SetProcessID(nt,kBlackProcess);
c9a8628a	254	}
	255	//
	256	// Black neutrons
	257	//
	258	fCharge = 0;
	259	kf = kNeutron;
1c22a81e	260	for(i = 0; i < fNbn; i++) {
5a9b7f8e	261	GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
642f15cf	262	PushTrack(fTrackIt, -1, kf, p, origin, polar,
58776b75	263	time, kPNoProcess, nt, 1.,-1);
c9a8628a	264	KeepTrack(nt);
87f95355	265	SetProcessID(nt,kBlackProcess);
c9a8628a	266	}
	267	}
	268
	269
	270
	271
5a9b7f8e	272	void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
5a9b7f8e	273	Double_t beta, Float_t* q, Float_t &theta)
f687abfe	274
f687abfe	275	{
c9a8628a	276	/*
	277	Emit a slow nucleon with "temperature" T [GeV],
	278	from a source moving with velocity beta
	279	Three-momentum [GeV/c] is given back in q[3]
	280	*/
	281
88a9f852	282	//printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
230d85c6	283
5a9b7f8e	284	Double_t m, pmax, p, f, phi;
c9a8628a	285	TDatabasePDG * pdg = TDatabasePDG::Instance();
	286	const Double_t kMassProton = pdg->GetParticle(kProton) ->Mass();
	287	const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
	288
	289	/* Select nucleon type */
	290	if(charge == 0) m = kMassNeutron;
	291	else m = kMassProton;
	292
	293	/* Momentum at maximum of Maxwell-distribution */
	294
88a9f852	295	pmax = TMath::Sqrt(2T(T+TMath::Sqrt(TT+mm)));
c9a8628a	296
	297	/* Try until proper momentum */
	298	/* for lack of primitive function of the Maxwell-distribution */
	299	/* a brute force trial-accept loop, normalized at pmax */
	300
	301	do
	302	{
	303	p = Rndm() * fPmax;
	304	f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
	305	}
	306	while(f < Rndm());
	307
5a9b7f8e	308	/* Spherical symmetric emission for black particles (beta=0)*/
	309	if(beta==0 \|\| fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
	310	/* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
	311	else if(fThetaDistribution!=0){
	312	Double_t costheta = fCosTheta->GetRandom();
	313	theta = TMath::ACos(costheta);
	314	}
	315	//
c9a8628a	316	phi = 2. * TMath::Pi() * Rndm();
c9a8628a	317
88a9f852	318
c9a8628a	319	/* Determine momentum components in system of the moving source */
	320	q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
	321	q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
	322	q[2] = p * TMath::Cos(theta);
74cfba91	323	//if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
c9a8628a	324
88a9f852	325
c9a8628a	326	/* Transform to system of the target nucleus */
	327	/* beta is passed as negative, because the gray nucleons are slowed down */
	328	Lorentz(m, -beta, q);
74cfba91	329	//if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
74cfba91	330
c9a8628a	331	/* Transform to laboratory system */
c9a8628a	332	Lorentz(m, fBeta, q);
83b3cd6c	333	q[2] *= fProtonDirection;
74cfba91	334	if(fDebug==1)printf("\n Momentum after LHC boost: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
	335
	336	if(fBeamCrossingAngle>0.) BeamCrossDivergence(1, q); // applying crossing angle
	337	if(fBeamDivergence>0.) BeamCrossDivergence(2, q); // applying divergence
	338
c9a8628a	339	}
	340
	341	Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
	342	{
	343	/* Relativistic Maxwell-distribution */
	344	Double_t ekin;
	345	ekin = TMath::Sqrt(pp+mm)-m;
	346	return (pp exp(-ekin/T));
	347	}
	348
	349
74cfba91	350	//_____________________________________________________________________________
c9a8628a	351	void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
	352	{
	353	/* Lorentz transform in the direction of q[2] */
	354
88a9f852	355	Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
88a9f852	356	Double_t energy = TMath::Sqrt(mm + q[0]q[0] + q[1]q[1] + q[2]q[2]);
c9a8628a	357	q[2] = gamma * (q[2] + beta*energy);
88a9f852	358	//printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
c9a8628a	359	}
230d85c6	360
74cfba91	361	//_____________________________________________________________________________
	362	void AliGenSlowNucleons::BeamCrossDivergence(Int_t iwhat, Float_t *pLab)
	363	{
	364	// Applying beam divergence and crossing angle
	365	//
	366	Double_t pmod = TMath::Sqrt(pLab[0]pLab[0]+pLab[1]pLab[1]+pLab[2]*pLab[2]);
	367
	368	Double_t tetdiv = 0.;
	369	Double_t fidiv = 0.;
	370	if(iwhat==1){
	371	tetdiv = fBeamCrossingAngle;
	372	fidiv = k2PI/4.;
	373	}
	374	else if(iwhat==2){
	375	tetdiv = fBeamDivEvent;
	376	fidiv = (gRandom->Rndm())*k2PI;
	377	}
	378
	379	Double_t tetpart = TMath::ATan2(TMath::Sqrt(pLab[0]pLab[0]+pLab[1]pLab[1]), pLab[2]);
	380	Double_t fipart=0.;
	381	if(TMath::Abs(pLab[1])>0. \|\| TMath::Abs(pLab[0])>0.) fipart = TMath::ATan2(pLab[1], pLab[0]);
	382	if(fipart<0.) {fipart = fipart+k2PI;}
	383	tetdiv = tetdiv*kRaddeg;
	384	fidiv = fidiv*kRaddeg;
	385	tetpart = tetpart*kRaddeg;
	386	fipart = fipart*kRaddeg;
	387
	388	Double_t angleSum[2]={0., 0.};
	389	AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum);
	390
	391	Double_t tetsum = angleSum[0];
	392	Double_t fisum = angleSum[1];
	393	//printf("tetpart %f fipart %f tetdiv %f fidiv %f angleSum %f %f\n",tetpart,fipart,tetdiv,fidiv,angleSum[0],angleSum[1]);
	394	tetsum = tetsum*kDegrad;
	395	fisum = fisum*kDegrad;
	396
	397	pLab[0] = pmodTMath::Sin(tetsum)TMath::Cos(fisum);
	398	pLab[1] = pmodTMath::Sin(tetsum)TMath::Sin(fisum);
	399	pLab[2] = pmod*TMath::Cos(tetsum);
	400	if(fDebug==1){
	401	if(iwhat==1) printf(" Beam crossing angle %f mrad ", fBeamCrossingAngle*1000.);
	402	else if(iwhat==2) printf(" Beam divergence %f mrad ", fBeamDivEvent*1000.);
	403	printf(" p = (%f, %f, %f)\n",pLab[0],pLab[1],pLab[2]);
	404	}
	405	}
	406
	407	//_____________________________________________________________________________
	408	void AliGenSlowNucleons::AddAngle(Double_t theta1, Double_t phi1,
	409	Double_t theta2, Double_t phi2, Double_t *angleSum)
	410	{
	411	// Calculating the sum of 2 angles
	412	Double_t temp = -1.;
	413	Double_t conv = 180./TMath::ACos(temp);
	414
	415	Double_t ct1 = TMath::Cos(theta1/conv);
	416	Double_t st1 = TMath::Sin(theta1/conv);
	417	Double_t cp1 = TMath::Cos(phi1/conv);
	418	Double_t sp1 = TMath::Sin(phi1/conv);
	419	Double_t ct2 = TMath::Cos(theta2/conv);
	420	Double_t st2 = TMath::Sin(theta2/conv);
	421	Double_t cp2 = TMath::Cos(phi2/conv);
	422	Double_t sp2 = TMath::Sin(phi2/conv);
	423	Double_t cx = ct1cp1st2cp2+st1cp1ct2-sp1st2*sp2;
	424	Double_t cy = ct1sp1st2cp2+st1sp1ct2+cp1st2*sp2;
425	Double_t cz = ct1ct2-st1st2*cp2;
426
427	Double_t rtetsum = TMath::ACos(cz);
428	Double_t tetsum = conv*rtetsum;
429	if(TMath::Abs(tetsum)<1e-4 \|\| tetsum==180.) return;
430
431	temp = cx/TMath::Sin(rtetsum);
432	if(temp>1.) temp=1.;
433	if(temp<-1.) temp=-1.;
434	Double_t fisum = conv*TMath::ACos(temp);
435	if(cy<0) {fisum = 360.-fisum;}
436
437	angleSum[0] = tetsum;
438	angleSum[1] = fisum;
439	}
440
87f95355	441	//_____________________________________________________________________________
	442	void AliGenSlowNucleons::SetProcessID(Int_t nt, UInt_t process)
	443	{
	444	// Tag the particle as
	445	// gray or black
	446	if (fStack)
	447	fStack->Particle(nt)->SetUniqueID(process);
	448	else
	449	gAlice->GetMCApp()->Particle(nt)->SetUniqueID(process);
	450	}