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