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 nucluons in pA interactions.
  20   Source is modelled by a relativistic Maxwell distributions.
  21   Original code by  Ferenc Sikler  <sikler@rmki.kfki.hu>
  22  */
  23
  24 #include <TDatabasePDG.h>
  25 #include <TPDGCode.h>
  26 #include <TH2F.h>
  27 #include <TCanvas.h>
  28
  29 #include "AliCollisionGeometry.h"
  30 #include "AliGenSlowNucleons.h"
  31 #include "AliSlowNucleonModel.h"
  32
  33  ClassImp(AliGenSlowNucleons)
  34
  35  AliGenSlowNucleons::AliGenSlowNucleons():AliGenerator(-1)
  36 {
  37 // Default constructor
  38     fSlowNucleonModel = 0;
  39     fCollisionGeometry = 0;
  40 }
  41
  42 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
  43     :AliGenerator(npart)
  44 {
  45 // Constructor
  46     fName  = "SlowNucleons";
  47     fTitle = "Generator for gray particles in pA collisions";
  48     SetPmax();
  49     SetTarget();
  50     SetNominalCmsEnergy();
  51     SetCharge();
  52     SetTemperature();
  53     SetBetaSource();
  54     fSlowNucleonModel = new AliSlowNucleonModel();
  55     fCollisionGeometry = 0;
  56     fDebug = 0;
  57 }
  58
  59 //____________________________________________________________
  60 AliGenSlowNucleons::~AliGenSlowNucleons()
  61 {
  62 // Destructor
  63     delete  fSlowNucleonModel;
  64 }
  65
  66
  67 void AliGenSlowNucleons::Init()
  68 {
  69   //
  70   // Initialization
  71   //
  72     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
  73     fMomentum = fCMS/2. * fZTarget / fATarget;
  74     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
  75     if (fDebug) {
  76         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
  77         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
  78     }
  79 }
  80
  81 void AliGenSlowNucleons::FinishRun()
  82 {
  83     if (fDebug) {
  84         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
  85         c->Divide(2,1);
  86         c->cd(1);
  87         fDebugHist1->Draw();
  88         c->cd(2);
  89         fDebugHist2->Draw();
  90     }
  91 }
  92
  93
  94 void AliGenSlowNucleons::Generate()
  95 {
  96   //
  97   // Generate one event
  98   //
  99   //
 100   // Communication with Gray Particle Model
 101   //
 102     if (fCollisionGeometry) {
 103         Float_t b   = fCollisionGeometry->ImpactParameter();
 104         Int_t  nn   = fCollisionGeometry->NN();
 105         Int_t  nwn  = fCollisionGeometry->NwN();
 106         Int_t  nnw  = fCollisionGeometry->NNw();
 107         Int_t  nwnw = fCollisionGeometry->NwNw();
 108
 109         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 110         if (fDebug) {
 111             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
 112             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
 113             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 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         PushTrack(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         PushTrack(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         PushTrack(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         PushTrack(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