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     Float_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     if (fDebug) {
 120         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
 121         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
 122         fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
 123     }
 124     //
 125     // non-uniform cos(theta) distribution
 126     //
 127     if(fThetaDistribution != 0) {
 128         fCosTheta = new TF1("fCosTheta",
 129                             "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
 130                             -1., 1.);
 131     }
 132 }
 133
 134 void AliGenSlowNucleons::FinishRun()
 135 {
 136 // End of run action
 137 // Show histogram for debugging if requested.
 138     if (fDebug) {
 139         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
 140         c->Divide(2,1);
 141         c->cd(1);
 142         fDebugHist1->Draw();
 143         c->cd(2);
 144         fDebugHist2->Draw();
 145         c->cd(3);
 146         fCosThetaGrayHist->Draw();
 147     }
 148 }
 149
 150
 151 void AliGenSlowNucleons::Generate()
 152 {
 153   //
 154   // Generate one event
 155   //
 156   //
 157   // Communication with Gray Particle Model
 158   //
 159     if (fCollisionGeometry) {
 160         Float_t b   = fCollisionGeometry->ImpactParameter();
 161         //      Int_t  nn   = fCollisionGeometry->NN();
 162         //      Int_t  nwn  = fCollisionGeometry->NwN();
 163         //      Int_t  nnw  = fCollisionGeometry->NNw();
 164         //      Int_t  nwnw = fCollisionGeometry->NwNw();
 165
 166         fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 167         if (fDebug) {
 168             printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
 169             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
 170             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
 171             //printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
 172         }
 173     }
 174
 175    //
 176     Float_t p[3], theta=0;
 177     Float_t origin[3] = {0., 0., 0.};
 178     Float_t time = 0.;
 179     Float_t polar [3] = {0., 0., 0.};
 180     Int_t nt, i, j;
 181     Int_t kf;
 182
 183     if(fVertexSmear == kPerEvent) {
 184         Vertex();
 185         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 186         time = fTime;
 187     } // if kPerEvent
 188 //
 189 //  Gray protons
 190 //
 191     fCharge = 1;
 192     kf = kProton;
 193     for(i = 0; i < fNgp; i++) {
 194         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 195         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 196         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 197                  time, kPNoProcess, nt, 1.);
 198         KeepTrack(nt);
 199     }
 200 //
 201 //  Gray neutrons
 202 //
 203     fCharge = 0;
 204     kf = kNeutron;
 205     for(i = 0; i < fNgn; i++) {
 206         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 207         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 208         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 209                  time, kPNoProcess, nt, 1.);
 210         KeepTrack(nt);
 211     }
 212 //
 213 //  Black protons
 214 //
 215     fCharge = 1;
 216     kf = kProton;
 217     for(i = 0; i < fNbp; i++) {
 218         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 219         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 220                  time, kPNoProcess, nt, 1.);
 221         KeepTrack(nt);
 222     }
 223 //
 224 //  Black neutrons
 225 //
 226     fCharge = 0;
 227     kf = kNeutron;
 228     for(i = 0; i < fNbn; i++) {
 229         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 230         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 231                  time, kPNoProcess, nt, 1.);
 232         KeepTrack(nt);
 233     }
 234 }
 235
 236
 237
 238
 239 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 240         Double_t beta, Float_t* q, Float_t &theta)
 241
 242 {
 243 /*
 244    Emit a slow nucleon with "temperature" T [GeV],
 245    from a source moving with velocity beta
 246    Three-momentum [GeV/c] is given back in q[3]
 247 */
 248
 249  //printf("Generating slow nuc. with: charge %d. temp. %1.3f, beta %1.2f\n",charge,T,beta);
 250
 251  Double_t m, pmax, p, f, phi;
 252  TDatabasePDG * pdg = TDatabasePDG::Instance();
 253  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 254  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 255
 256  /* Select nucleon type */
 257  if(charge == 0) m = kMassNeutron;
 258  else m = kMassProton;
 259
 260  /* Momentum at maximum of Maxwell-distribution */
 261
 262  pmax = TMath::Sqrt(2*T*(T+sqrt(T*T+m*m)));
 263
 264  /* Try until proper momentum                                  */
 265  /* for lack of primitive function of the Maxwell-distribution */
 266  /* a brute force trial-accept loop, normalized at pmax        */
 267
 268  do
 269  {
 270      p = Rndm() * fPmax;
 271      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 272  }
 273  while(f < Rndm());
 274
 275  /* Spherical symmetric emission for black particles (beta=0)*/
 276  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 277  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 278  else if(fThetaDistribution!=0){
 279    Double_t costheta = fCosTheta->GetRandom();
 280    theta = TMath::ACos(costheta);
 281  }
 282  //
 283  phi   = 2. * TMath::Pi() * Rndm();
 284
 285  /* Determine momentum components in system of the moving source */
 286  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 287  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 288  q[2] = p * TMath::Cos(theta);
 289
 290  /* Transform to system of the target nucleus                             */
 291  /* beta is passed as negative, because the gray nucleons are slowed down */
 292  Lorentz(m, -beta, q);
 293
 294  /* Transform to laboratory system */
 295  Lorentz(m, fBeta, q);
 296  q[2] *= fProtonDirection;
 297 }
 298
 299 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 300 {
 301 /* Relativistic Maxwell-distribution */
 302     Double_t ekin;
 303     ekin = TMath::Sqrt(p*p+m*m)-m;
 304     return (p*p * exp(-ekin/T));
 305 }
 306
 307
 308 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 309 {
 310 /* Lorentz transform in the direction of q[2] */
 311
 312     Double_t gamma  = 1./sqrt((1.-beta)*(1.+beta));
 313     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 314     q[2] = gamma * (q[2] + beta*energy);
 315 }
 316