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(0.),
  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(0.),
  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 polar [3] = {0., 0., 0.};
 184     Int_t nt, i, j;
 185     Int_t kf;
 186
 187     if(fVertexSmear == kPerEvent) {
 188         Vertex();
 189         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 190     } // if kPerEvent
 191 //
 192 //  Gray protons
 193 //
 194     fCharge = 1;
 195     kf = kProton;
 196     for(i = 0; i < fNgp; i++) {
 197         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 198         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 199         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 200                  0., kPNoProcess, nt, 1.);
 201         KeepTrack(nt);
 202     }
 203 //
 204 //  Gray neutrons
 205 //
 206     fCharge = 0;
 207     kf = kNeutron;
 208     for(i = 0; i < fNgn; i++) {
 209         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 210         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 211         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 212                  0., kPNoProcess, nt, 1.);
 213         KeepTrack(nt);
 214     }
 215 //
 216 //  Black protons
 217 //
 218     fCharge = 1;
 219     kf = kProton;
 220     for(i = 0; i < fNbp; i++) {
 221         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 222         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 223                  0., kPNoProcess, nt, 1.);
 224         KeepTrack(nt);
 225     }
 226 //
 227 //  Black neutrons
 228 //
 229     fCharge = 0;
 230     kf = kNeutron;
 231     for(i = 0; i < fNbn; i++) {
 232         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 233         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 234                  0., kPNoProcess, nt, 1.);
 235         KeepTrack(nt);
 236     }
 237 }
 238
 239
 240
 241
 242 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 243         Double_t beta, Float_t* q, Float_t &theta)
 244
 245 {
 246 /*
 247    Emit a slow nucleon with "temperature" T [GeV],
 248    from a source moving with velocity beta
 249    Three-momentum [GeV/c] is given back in q[3]
 250 */
 251
 252  Double_t m, pmax, p, f, phi;
 253  TDatabasePDG * pdg = TDatabasePDG::Instance();
 254  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 255  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 256
 257  /* Select nucleon type */
 258  if(charge == 0) m = kMassNeutron;
 259  else m = kMassProton;
 260
 261  /* Momentum at maximum of Maxwell-distribution */
 262
 263  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
 264
 265  /* Try until proper momentum                                  */
 266  /* for lack of primitive function of the Maxwell-distribution */
 267  /* a brute force trial-accept loop, normalized at pmax        */
 268
 269  do
 270  {
 271      p = Rndm() * fPmax;
 272      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 273  }
 274  while(f < Rndm());
 275
 276  /* Spherical symmetric emission for black particles (beta=0)*/
 277  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 278  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 279  else if(fThetaDistribution!=0){
 280    Double_t costheta = fCosTheta->GetRandom();
 281    theta = TMath::ACos(costheta);
 282  }
 283  //
 284  phi   = 2. * TMath::Pi() * Rndm();
 285
 286  /* Determine momentum components in system of the moving source */
 287  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 288  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 289  q[2] = p * TMath::Cos(theta);
 290
 291  /* Transform to system of the target nucleus                             */
 292  /* beta is passed as negative, because the gray nucleons are slowed down */
 293  Lorentz(m, -beta, q);
 294
 295  /* Transform to laboratory system */
 296  Lorentz(m, fBeta, q);
 297  q[2] *= fProtonDirection;
 298 }
 299
 300 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 301 {
 302 /* Relativistic Maxwell-distribution */
 303     Double_t ekin;
 304     ekin = TMath::Sqrt(p*p+m*m)-m;
 305     return (p*p * exp(-ekin/T));
 306 }
 307
 308
 309 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 310 {
 311 /* Lorentz transform in the direction of q[2] */
 312
 313     Double_t gamma  = 1/sqrt((1.-beta)*(1.+beta));
 314     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 315     q[2] = gamma * (q[2] + beta*energy);
 316 }