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