EVGEN/AliGenSlowNucleons.cxx

   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
  16 /* $Id$ */
  17
  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 //
  26
  27 #include <TDatabasePDG.h>
  28 #include <TPDGCode.h>
  29 #include <TH2F.h>
  30 #include <TH1F.h>
  31 #include <TF1.h>
  32 #include <TCanvas.h>
  33
  34 #include "AliCollisionGeometry.h"
  35 #include "AliGenSlowNucleons.h"
  36 #include "AliSlowNucleonModel.h"
  37
  38 ClassImp(AliGenSlowNucleons)
  39
  40
  41 AliGenSlowNucleons::AliGenSlowNucleons()
  42     :AliGenerator(-1),
  43      fCMS(0.),
  44      fMomentum(0.),
  45      fBeta(0.),
  46      fPmax (0.),
  47      fATarget (0.),
  48      fZTarget (0.),
  49      fCharge(0),
  50      fProtonDirection(1.),
  51      fTemperatureG(0.),
  52      fBetaSourceG(0.),
  53      fTemperatureB(0.),
  54      fBetaSourceB(0.),
  55      fNgp(0),
  56      fNgn(0),
  57      fNbp(0),
  58      fNbn(0),
  59      fDebug(0),
  60      fDebugHist1(0),
  61      fDebugHist2(0),
  62      fThetaDistribution(),
  63      fCosThetaGrayHist(),
  64      fCosTheta(),
  65      fSlowNucleonModel(0)
  66 {
  67 // Default constructor
  68     fCollisionGeometry = 0;
  69 }
  70
  71 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
  72     :AliGenerator(npart),
  73      fCMS(14000.),
  74      fMomentum(0.),
  75      fBeta(0.),
  76      fPmax (10.),
  77      fATarget (208.),
  78      fZTarget (82.),
  79      fCharge(1),
  80      fProtonDirection(1.),
  81      fTemperatureG(0.04),
  82      fBetaSourceG(0.05),
  83      fTemperatureB(0.004),
  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(),
  95      fSlowNucleonModel(new AliSlowNucleonModel())
  96 {
  97 // Constructor
  98     fName  = "SlowNucleons";
  99     fTitle = "Generator for gray particles in pA collisions";
 100     fCollisionGeometry = 0;
 101 }
 102
 103 //____________________________________________________________
 104 AliGenSlowNucleons::~AliGenSlowNucleons()
 105 {
 106 // Destructor
 107     delete  fSlowNucleonModel;
 108 }
 109
 110 void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
 111 // Set direction of the proton to change between pA (1) and Ap (-1)
 112   fProtonDirection = dir / TMath::Abs(dir);
 113 }
 114
 115 void AliGenSlowNucleons::Init()
 116 {
 117   //
 118   // Initialization
 119   //
 120     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
 121     fMomentum = fCMS/2. * fZTarget / fATarget;
 122     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
 123     if (fDebug) {
 124         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
 125         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
 126         fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
 127     }
 128     //
 129     // non-uniform cos(theta) distribution
 130     //
 131     if(fThetaDistribution != 0) {
 132         fCosTheta = new TF1("fCosTheta",
 133                             "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
 134                             -1., 1.);
 135     }
 136 }
 137
 138 void AliGenSlowNucleons::FinishRun()
 139 {
 140 // End of run action
 141 // Show histogram for debugging if requested.
 142     if (fDebug) {
 143         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
 144         c->Divide(2,1);
 145         c->cd(1);
 146         fDebugHist1->Draw();
 147         c->cd(2);
 148         fDebugHist2->Draw();
 149         c->cd(3);
 150         fCosThetaGrayHist->Draw();
 151     }
 152 }
 153
 154
 155 void AliGenSlowNucleons::Generate()
 156 {
 157   //
 158   // Generate one event
 159   //
 160   //
 161   // Communication with Gray Particle Model
 162   //
 163     if (fCollisionGeometry) {
 164         Float_t b   = fCollisionGeometry->ImpactParameter();
 165         Int_t  nn   = fCollisionGeometry->NN();
 166         Int_t  nwn  = fCollisionGeometry->NwN();
 167         Int_t  nnw  = fCollisionGeometry->NNw();
 168         Int_t  nwnw = fCollisionGeometry->NwNw();
 169
 170         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 171         if (fDebug) {
 172             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
 173             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
 174             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
 175             printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
 176                    b, nn, nwn, nnw, nwnw);
 177         }
 178     }
 179
 180    //
 181     Float_t p[3], theta=0;
 182     Float_t origin[3] = {0., 0., 0.};
 183     Float_t time = 0.;
 184     Float_t polar [3] = {0., 0., 0.};
 185     Int_t nt, i, j;
 186     Int_t kf;
 187
 188     if(fVertexSmear == kPerEvent) {
 189         Vertex();
 190         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 191         time = fTime;
 192     } // if kPerEvent
 193 //
 194 //  Gray protons
 195 //
 196     fCharge = 1;
 197     kf = kProton;
 198     for(i = 0; i < fNgp; i++) {
 199         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 200         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 201         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 202                  time, kPNoProcess, nt, 1.);
 203         KeepTrack(nt);
 204     }
 205 //
 206 //  Gray neutrons
 207 //
 208     fCharge = 0;
 209     kf = kNeutron;
 210     for(i = 0; i < fNgn; i++) {
 211         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 212         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 213         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 214                  time, kPNoProcess, nt, 1.);
 215         KeepTrack(nt);
 216     }
 217 //
 218 //  Black protons
 219 //
 220     fCharge = 1;
 221     kf = kProton;
 222     for(i = 0; i < fNbp; i++) {
 223         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 224         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 225                  time, kPNoProcess, nt, 1.);
 226         KeepTrack(nt);
 227     }
 228 //
 229 //  Black neutrons
 230 //
 231     fCharge = 0;
 232     kf = kNeutron;
 233     for(i = 0; i < fNbn; i++) {
 234         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 235         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 236                  time, kPNoProcess, nt, 1.);
 237         KeepTrack(nt);
 238     }
 239 }
 240
 241
 242
 243
 244 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 245         Double_t beta, Float_t* q, Float_t &theta)
 246
 247 {
 248 /*
 249    Emit a slow nucleon with "temperature" T [GeV],
 250    from a source moving with velocity beta
 251    Three-momentum [GeV/c] is given back in q[3]
 252 */
 253
 254  Double_t m, pmax, p, f, phi;
 255  TDatabasePDG * pdg = TDatabasePDG::Instance();
 256  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 257  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 258
 259  /* Select nucleon type */
 260  if(charge == 0) m = kMassNeutron;
 261  else m = kMassProton;
 262
 263  /* Momentum at maximum of Maxwell-distribution */
 264
 265  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
 266
 267  /* Try until proper momentum                                  */
 268  /* for lack of primitive function of the Maxwell-distribution */
 269  /* a brute force trial-accept loop, normalized at pmax        */
 270
 271  do
 272  {
 273      p = Rndm() * fPmax;
 274      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 275  }
 276  while(f < Rndm());
 277
 278  /* Spherical symmetric emission for black particles (beta=0)*/
 279  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 280  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 281  else if(fThetaDistribution!=0){
 282    Double_t costheta = fCosTheta->GetRandom();
 283    theta = TMath::ACos(costheta);
 284  }
 285  //
 286  phi   = 2. * TMath::Pi() * Rndm();
 287
 288  /* Determine momentum components in system of the moving source */
 289  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 290  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 291  q[2] = p * TMath::Cos(theta);
 292
 293  /* Transform to system of the target nucleus                             */
 294  /* beta is passed as negative, because the gray nucleons are slowed down */
 295  Lorentz(m, -beta, q);
 296
 297  /* Transform to laboratory system */
 298  Lorentz(m, fBeta, q);
 299  q[2] *= fProtonDirection;
 300 }
 301
 302 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 303 {
 304 /* Relativistic Maxwell-distribution */
 305     Double_t ekin;
 306     ekin = TMath::Sqrt(p*p+m*m)-m;
 307     return (p*p * exp(-ekin/T));
 308 }
 309
 310
 311 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 312 {
 313 /* Lorentz transform in the direction of q[2] */
 314
 315     Double_t gamma  = 1./sqrt((1.-beta)*(1.+beta));
 316     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 317     q[2] = gamma * (q[2] + beta*energy);
 318 }