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 /*
  17 $Log$
  18 Revision 1.2  2003/03/25 09:55:24  morsch
  19 Numbers of slow nucleons either from model or user set.
  20
  21 Revision 1.1  2003/03/24 16:49:23  morsch
  22 Slow nucleon generator and model.
  23
  24 */
  25
  26 /*
  27   Generator for slow nucluons in pA interactions.
  28   Source is modelled by a relativistic Maxwell distributions.
  29   Original code by  Ferenc Sikler  <sikler@rmki.kfki.hu>
  30  */
  31
  32 #include <TDatabasePDG.h>
  33 #include <TPDGCode.h>
  34 #include <TH2F.h>
  35
  36 #include "AliCollisionGeometry.h"
  37 #include "AliGenSlowNucleons.h"
  38 #include "AliSlowNucleonModel.h"
  39
  40  ClassImp(AliGenSlowNucleons)
  41
  42  AliGenSlowNucleons::AliGenSlowNucleons():AliGenerator(-1)
  43 {
  44 // Default constructor
  45     fSlowNucleonModel = 0;
  46     fCollisionGeometry = 0;
  47 }
  48
  49 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
  50     :AliGenerator(npart)
  51 {
  52 // Constructor
  53     fName  = "SlowNucleons";
  54     fTitle = "Generator for gray particles in pA collisions";
  55     SetPmax();
  56     SetTarget();
  57     SetNominalCmsEnergy();
  58     SetCharge();
  59     SetTemperature();
  60     SetBetaSource();
  61     fSlowNucleonModel = new AliSlowNucleonModel();
  62     fCollisionGeometry = 0;
  63     fDebug = 0;
  64 }
  65
  66 //____________________________________________________________
  67 AliGenSlowNucleons::~AliGenSlowNucleons()
  68 {
  69 // Destructor
  70     delete  fSlowNucleonModel;
  71 }
  72
  73
  74 void AliGenSlowNucleons::Init()
  75 {
  76   //
  77   // Initialization
  78   //
  79     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
  80     fMomentum = fCMS/2. * fZTarget / fATarget;
  81     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
  82     if (fDebug) {
  83         fDebugHist = new TH2F("DebugHist", "N_heavy vs nu", 50, 0., 50., 50, 0., 50.);
  84     }
  85 }
  86
  87 void AliGenSlowNucleons::FinishRun()
  88 {
  89     if (fDebug) {
  90         fDebugHist->Draw();
  91     }
  92 }
  93
  94
  95 void AliGenSlowNucleons::Generate()
  96 {
  97   //
  98   // Generate one event
  99   //
 100   //
 101   // Communication with Gray Particle Model
 102   //
 103     if (fCollisionGeometry) {
 104         Float_t b   = fCollisionGeometry->ImpactParameter();
 105         Int_t  nn   = fCollisionGeometry->NN();
 106         Int_t  nwn  = fCollisionGeometry->NwN();
 107         Int_t  nnw  = fCollisionGeometry->NNw();
 108         Int_t  nwnw = fCollisionGeometry->NwNw();
 109
 110         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 111         if (fDebug) {
 112             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
 113             fDebugHist->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
 114             printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
 115                    b, nn, nwn, nnw, nwnw);
 116         }
 117     }
 118
 119    //
 120     Float_t p[3];
 121     Float_t origin[3] = {0., 0., 0.};
 122     Float_t polar [3] = {0., 0., 0.};
 123     Int_t nt, i, j;
 124     Int_t kf;
 125
 126     if(fVertexSmear == kPerEvent) {
 127         Vertex();
 128         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 129     } // if kPerEvent
 130 //
 131 //  Gray protons
 132 //
 133     fCharge = 1;
 134     kf = kProton;
 135     for(i = 0; i < fNgp; i++) {
 136         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
 137         SetTrack(fTrackIt, -1, kf, p, origin, polar,
 138                  0., kPNoProcess, nt, 1.);
 139         KeepTrack(nt);
 140     }
 141 //
 142 //  Gray neutrons
 143 //
 144     fCharge = 0;
 145     kf = kNeutron;
 146     for(i = 0; i < fNgn; i++) {
 147         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
 148         SetTrack(fTrackIt, -1, kf, p, origin, polar,
 149                  0., kPNoProcess, nt, 1.);
 150         KeepTrack(nt);
 151     }
 152 //
 153 //  Black protons
 154 //
 155     fCharge = 1;
 156     kf = kProton;
 157     for(i = 0; i < fNbp; i++) {
 158         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
 159         SetTrack(fTrackIt, -1, kf, p, origin, polar,
 160                  0., kPNoProcess, nt, 1.);
 161         KeepTrack(nt);
 162     }
 163 //
 164 //  Black neutrons
 165 //
 166     fCharge = 0;
 167     kf = kNeutron;
 168     for(i = 0; i < fNbn; i++) {
 169         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
 170         SetTrack(fTrackIt, -1, kf, p, origin, polar,
 171                  0., kPNoProcess, nt, 1.);
 172         KeepTrack(nt);
 173     }
 174 }
 175
 176
 177
 178
 179 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, Double_t beta, Float_t* q)
 180 /*
 181    Emit a slow nucleon with "temperature" T [GeV],
 182    from a source moving with velocity beta
 183    Three-momentum [GeV/c] is given back in q[3]
 184 */
 185
 186 {
 187  Double_t m, pmax, p, f, theta, phi;
 188
 189  TDatabasePDG * pdg = TDatabasePDG::Instance();
 190  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 191  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 192
 193  /* Select nucleon type */
 194  if(charge == 0) m = kMassNeutron;
 195  else m = kMassProton;
 196
 197  /* Momentum at maximum of Maxwell-distribution */
 198
 199  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
 200
 201  /* Try until proper momentum                                  */
 202  /* for lack of primitive function of the Maxwell-distribution */
 203  /* a brute force trial-accept loop, normalized at pmax        */
 204
 205  do
 206  {
 207      p = Rndm() * fPmax;
 208      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 209  }
 210  while(f < Rndm());
 211
 212  /* Spherical symmetric emission */
 213  theta = TMath::ACos(2. * Rndm() - 1.);
 214  phi   = 2. * TMath::Pi() * Rndm();
 215
 216  /* Determine momentum components in system of the moving source */
 217  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 218  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 219  q[2] = p * TMath::Cos(theta);
 220
 221  /* Transform to system of the target nucleus                             */
 222  /* beta is passed as negative, because the gray nucleons are slowed down */
 223  Lorentz(m, -beta, q);
 224
 225  /* Transform to laboratory system */
 226  Lorentz(m, fBeta, q);
 227 }
 228
 229 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 230 {
 231 /* Relativistic Maxwell-distribution */
 232     Double_t ekin;
 233     ekin = TMath::Sqrt(p*p+m*m)-m;
 234     return (p*p * exp(-ekin/T));
 235 }
 236
 237
 238 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 239 {
 240 /* Lorentz transform in the direction of q[2] */
 241
 242     Double_t gamma  = 1/sqrt(1-beta*beta);
 243     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 244     q[2] = gamma * (q[2] + beta*energy);
 245 }
 246
 247
 248
 249
 250
 251
 252