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      fCharge(0),
  48      fProtonDirection(1.),
  49      fTemperatureG(0.),
  50      fBetaSourceG(0.),
  51      fTemperatureB(0.),
  52      fBetaSourceB(0.),
  53      fNgp(0),
  54      fNgn(0),
  55      fNbp(0),
  56      fNbn(0),
  57      fDebug(0),
  58      fDebugHist1(0),
  59      fDebugHist2(0),
  60      fThetaDistribution(),
  61      fCosThetaGrayHist(),
  62      fCosTheta(),
  63      fSlowNucleonModel(0)
  64 {
  65 // Default constructor
  66     fCollisionGeometry = 0;
  67 }
  68
  69 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
  70     :AliGenerator(npart),
  71      fCMS(14000.),
  72      fMomentum(0.),
  73      fBeta(0.),
  74      fPmax (10.),
  75      fCharge(1),
  76      fProtonDirection(1.),
  77      fTemperatureG(0.05),
  78      fBetaSourceG(0.05),
  79      fTemperatureB(0.005),
  80      fBetaSourceB(0.),
  81      fNgp(0),
  82      fNgn(0),
  83      fNbp(0),
  84      fNbn(0),
  85      fDebug(0),
  86      fDebugHist1(0),
  87      fDebugHist2(0),
  88      fThetaDistribution(),
  89      fCosThetaGrayHist(),
  90      fCosTheta(),
  91      fSlowNucleonModel(new AliSlowNucleonModel())
  92 {
  93 // Constructor
  94     fName  = "SlowNucleons";
  95     fTitle = "Generator for gray particles in pA collisions";
  96     fCollisionGeometry = 0;
  97 }
  98
  99 //____________________________________________________________
 100 AliGenSlowNucleons::~AliGenSlowNucleons()
 101 {
 102 // Destructor
 103     delete  fSlowNucleonModel;
 104 }
 105
 106 void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
 107 // Set direction of the proton to change between pA (1) and Ap (-1)
 108   fProtonDirection = dir / TMath::Abs(dir);
 109 }
 110
 111 void AliGenSlowNucleons::Init()
 112 {
 113   //
 114   // Initialization
 115   //
 116     Double_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
 117     fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
 118     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
 119     //printf("  fMomentum %f    fBeta %1.10f\n",fMomentum, fBeta);
 120     if (fDebug) {
 121         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
 122         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
 123         fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
 124     }
 125     //
 126     // non-uniform cos(theta) distribution
 127     //
 128     if(fThetaDistribution != 0) {
 129         fCosTheta = new TF1("fCosTheta",
 130                             "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
 131                             -1., 1.);
 132     }
 133 }
 134
 135 void AliGenSlowNucleons::FinishRun()
 136 {
 137 // End of run action
 138 // Show histogram for debugging if requested.
 139     if (fDebug) {
 140         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
 141         c->Divide(2,1);
 142         c->cd(1);
 143         fDebugHist1->Draw();
 144         c->cd(2);
 145         fDebugHist2->Draw();
 146         c->cd(3);
 147         fCosThetaGrayHist->Draw();
 148     }
 149 }
 150
 151
 152 void AliGenSlowNucleons::Generate()
 153 {
 154   //
 155   // Generate one event
 156   //
 157   //
 158   // Communication with Gray Particle Model
 159   //
 160     if (fCollisionGeometry) {
 161         Float_t b   = fCollisionGeometry->ImpactParameter();
 162         //      Int_t  nn   = fCollisionGeometry->NN();
 163         //      Int_t  nwn  = fCollisionGeometry->NwN();
 164         //      Int_t  nnw  = fCollisionGeometry->NNw();
 165         //      Int_t  nwnw = fCollisionGeometry->NwNw();
 166
 167         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 168         if (fDebug) {
 169             //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
 170             printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
 171             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
 172             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
 173         }
 174     }
 175
 176    //
 177     Float_t p[3] = {0., 0., 0.}, theta=0;
 178     Float_t origin[3] = {0., 0., 0.};
 179     Float_t time = 0.;
 180     Float_t polar [3] = {0., 0., 0.};
 181     Int_t nt, i, j;
 182     Int_t kf;
 183
 184     if(fVertexSmear == kPerEvent) {
 185         Vertex();
 186         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 187         time = fTime;
 188     } // if kPerEvent
 189 //
 190 //  Gray protons
 191 //
 192     fCharge = 1;
 193     kf = kProton;
 194     for(i = 0; i < fNgp; i++) {
 195         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 196         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 197         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 198                  time, kPNoProcess, nt, 1.,-1);
 199         KeepTrack(nt);
 200     }
 201 //
 202 //  Gray neutrons
 203 //
 204     fCharge = 0;
 205     kf = kNeutron;
 206     for(i = 0; i < fNgn; i++) {
 207         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 208         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 209         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 210                  time, kPNoProcess, nt, 1.,-1);
 211         KeepTrack(nt);
 212     }
 213 //
 214 //  Black protons
 215 //
 216     fCharge = 1;
 217     kf = kProton;
 218     for(i = 0; i < fNbp; i++) {
 219         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 220         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 221                  time, kPNoProcess, nt, 1.,-1);
 222         KeepTrack(nt);
 223     }
 224 //
 225 //  Black neutrons
 226 //
 227     fCharge = 0;
 228     kf = kNeutron;
 229     for(i = 0; i < fNbn; i++) {
 230         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 231         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 232                  time, kPNoProcess, nt, 1.,-1);
 233         KeepTrack(nt);
 234     }
 235 }
 236
 237
 238
 239
 240 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 241         Double_t beta, Float_t* q, Float_t &theta)
 242
 243 {
 244 /*
 245    Emit a slow nucleon with "temperature" T [GeV],
 246    from a source moving with velocity beta
 247    Three-momentum [GeV/c] is given back in q[3]
 248 */
 249
 250  //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
 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+TMath::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
 287  /* Determine momentum components in system of the moving source */
 288  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 289  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 290  q[2] = p * TMath::Cos(theta);
 291
 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
 303 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 304 {
 305 /* Relativistic Maxwell-distribution */
 306     Double_t ekin;
 307     ekin = TMath::Sqrt(p*p+m*m)-m;
 308     return (p*p * exp(-ekin/T));
 309 }
 310
 311
 312 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 313 {
 314 /* Lorentz transform in the direction of q[2] */
 315
 316     Double_t gamma  = 1./TMath::Sqrt((1.-beta)*(1.+beta));
 317     Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 318     q[2] = gamma * (q[2] + beta*energy);
 319     //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
 320 }
 321