[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 "AliConst.h"
#include "AliCollisionGeometry.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.),
     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.),
     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();
	
	//fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
	fSlowNucleonModel->GetNumberOfSlowNucleons2(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);
    }
//
//  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);
    }
//
//  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);
    }
//
//  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);
    }
}


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