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(const AliGenSlowNucleons & sn):
 105     AliGenerator(sn),
 106     fCMS(0.),
 107     fMomentum(0.),
 108     fBeta(0.),
 109     fPmax (0.),
 110     fATarget (0.),
 111     fZTarget (0.),
 112     fCharge(0),
 113     fProtonDirection(0.),
 114     fTemperatureG(0.),
 115     fBetaSourceG(0.),
 116     fTemperatureB(0.),
 117     fBetaSourceB(0.),
 118     fNgp(0),
 119     fNgn(0),
 120     fNbp(0),
 121     fNbn(0),
 122     fDebug(0),
 123     fDebugHist1(0),
 124     fDebugHist2(0),
 125     fThetaDistribution(),
 126     fCosThetaGrayHist(),
 127     fCosTheta(),
 128     fSlowNucleonModel(0)
 129 {
 130 // Copy constructor
 131     sn.Copy(*this);
 132 }
 133
 134 //____________________________________________________________
 135 AliGenSlowNucleons::~AliGenSlowNucleons()
 136 {
 137 // Destructor
 138     delete  fSlowNucleonModel;
 139 }
 140
 141 void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
 142 // Set direction of the proton to change between pA (1) and Ap (-1)
 143   fProtonDirection = dir / TMath::Abs(dir);
 144 }
 145
 146 void AliGenSlowNucleons::Init()
 147 {
 148   //
 149   // Initialization
 150   //
 151     Float_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
 152     fMomentum = fCMS/2. * fZTarget / fATarget;
 153     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
 154     if (fDebug) {
 155         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
 156         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
 157         fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
 158     }
 159     //
 160     // non-uniform cos(theta) distribution
 161     //
 162     if(fThetaDistribution != 0) {
 163         fCosTheta = new TF1("fCosTheta",
 164                             "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
 165                             -1., 1.);
 166     }
 167 }
 168
 169 void AliGenSlowNucleons::FinishRun()
 170 {
 171 // End of run action
 172 // Show histogram for debugging if requested.
 173     if (fDebug) {
 174         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
 175         c->Divide(2,1);
 176         c->cd(1);
 177         fDebugHist1->Draw();
 178         c->cd(2);
 179         fDebugHist2->Draw();
 180         c->cd(3);
 181         fCosThetaGrayHist->Draw();
 182     }
 183 }
 184
 185
 186 void AliGenSlowNucleons::Generate()
 187 {
 188   //
 189   // Generate one event
 190   //
 191   //
 192   // Communication with Gray Particle Model
 193   //
 194     if (fCollisionGeometry) {
 195         Float_t b   = fCollisionGeometry->ImpactParameter();
 196         Int_t  nn   = fCollisionGeometry->NN();
 197         Int_t  nwn  = fCollisionGeometry->NwN();
 198         Int_t  nnw  = fCollisionGeometry->NNw();
 199         Int_t  nwnw = fCollisionGeometry->NwNw();
 200
 201         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 202         if (fDebug) {
 203             printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
 204             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
 205             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
 206             printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n",
 207                    b, nn, nwn, nnw, nwnw);
 208         }
 209     }
 210
 211    //
 212     Float_t p[3], theta=0;
 213     Float_t origin[3] = {0., 0., 0.};
 214     Float_t polar [3] = {0., 0., 0.};
 215     Int_t nt, i, j;
 216     Int_t kf;
 217
 218     if(fVertexSmear == kPerEvent) {
 219         Vertex();
 220         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 221     } // if kPerEvent
 222 //
 223 //  Gray protons
 224 //
 225     fCharge = 1;
 226     kf = kProton;
 227     for(i = 0; i < fNgp; i++) {
 228         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 229         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 230         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 231                  0., kPNoProcess, nt, 1.);
 232         KeepTrack(nt);
 233     }
 234 //
 235 //  Gray neutrons
 236 //
 237     fCharge = 0;
 238     kf = kNeutron;
 239     for(i = 0; i < fNgn; i++) {
 240         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 241         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 242         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 243                  0., kPNoProcess, nt, 1.);
 244         KeepTrack(nt);
 245     }
 246 //
 247 //  Black protons
 248 //
 249     fCharge = 1;
 250     kf = kProton;
 251     for(i = 0; i < fNbp; i++) {
 252         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 253         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 254                  0., kPNoProcess, nt, 1.);
 255         KeepTrack(nt);
 256     }
 257 //
 258 //  Black neutrons
 259 //
 260     fCharge = 0;
 261     kf = kNeutron;
 262     for(i = 0; i < fNbn; i++) {
 263         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 264         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 265                  0., kPNoProcess, nt, 1.);
 266         KeepTrack(nt);
 267     }
 268 }
 269
 270
 271
 272
 273 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 274         Double_t beta, Float_t* q, Float_t &theta)
 275
 276 {
 277 /*
 278    Emit a slow nucleon with "temperature" T [GeV],
 279    from a source moving with velocity beta
 280    Three-momentum [GeV/c] is given back in q[3]
 281 */
 282
 283  Double_t m, pmax, p, f, phi;
 284  TDatabasePDG * pdg = TDatabasePDG::Instance();
 285  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 286  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 287
 288  /* Select nucleon type */
 289  if(charge == 0) m = kMassNeutron;
 290  else m = kMassProton;
 291
 292  /* Momentum at maximum of Maxwell-distribution */
 293
 294  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
 295
 296  /* Try until proper momentum                                  */
 297  /* for lack of primitive function of the Maxwell-distribution */
 298  /* a brute force trial-accept loop, normalized at pmax        */
 299
 300  do
 301  {
 302      p = Rndm() * fPmax;
 303      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 304  }
 305  while(f < Rndm());
 306
 307  /* Spherical symmetric emission for black particles (beta=0)*/
 308  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 309  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 310  else if(fThetaDistribution!=0){
 311    Double_t costheta = fCosTheta->GetRandom();
 312    theta = TMath::ACos(costheta);
 313  }
 314  //
 315  phi   = 2. * TMath::Pi() * Rndm();
 316
 317  /* Determine momentum components in system of the moving source */
 318  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 319  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 320  q[2] = p * TMath::Cos(theta);
 321
 322  /* Transform to system of the target nucleus                             */
 323  /* beta is passed as negative, because the gray nucleons are slowed down */
 324  Lorentz(m, -beta, q);
 325
 326  /* Transform to laboratory system */
 327  Lorentz(m, fBeta, q);
 328  q[2] *= fProtonDirection;
 329 }
 330
 331 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 332 {
 333 /* Relativistic Maxwell-distribution */
 334     Double_t ekin;
 335     ekin = TMath::Sqrt(p*p+m*m)-m;
 336     return (p*p * exp(-ekin/T));
 337 }
 338
 339
 340 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 341 {
 342 /* Lorentz transform in the direction of q[2] */
 343
 344     Double_t gamma  = 1/sqrt(1-beta*beta);
 345     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 346     q[2] = gamma * (q[2] + beta*energy);
 347 }
 348
 349
 350 AliGenSlowNucleons& AliGenSlowNucleons::operator=(const  AliGenSlowNucleons& rhs)
 351 {
 352 // Assignment operator
 353     rhs.Copy(*this);
 354     return *this;
 355 }
 356
 357 void AliGenSlowNucleons::Copy(TObject&) const
 358 {
 359     //
 360     // Copy
 361     //
 362     Fatal("Copy","Not implemented!\n");
 363 }
 364
 365
 366
 367
 368
 369
 370