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 #include <TParticle.h>
  34
  35 #include "AliConst.h"
  36 #include "AliCollisionGeometry.h"
  37 #include "AliStack.h"
  38 #include "AliRun.h"
  39 #include "AliMC.h"
  40 #include "AliGenSlowNucleons.h"
  41 #include "AliSlowNucleonModel.h"
  42
  43 ClassImp(AliGenSlowNucleons)
  44
  45
  46 AliGenSlowNucleons::AliGenSlowNucleons()
  47     :AliGenerator(-1),
  48      fCMS(0.),
  49      fMomentum(0.),
  50      fBeta(0.),
  51      fPmax (0.),
  52      fCharge(0),
  53      fProtonDirection(1.),
  54      fTemperatureG(0.),
  55      fBetaSourceG(0.),
  56      fTemperatureB(0.),
  57      fBetaSourceB(0.),
  58      fNgp(0),
  59      fNgn(0),
  60      fNbp(0),
  61      fNbn(0),
  62      fDebug(0),
  63      fDebugHist1(0),
  64      fDebugHist2(0),
  65      fThetaDistribution(),
  66      fCosThetaGrayHist(),
  67      fCosTheta(),
  68      fBeamCrossingAngle(0.),
  69      fBeamDivergence(0.),
  70      fBeamDivEvent(0.),
  71      fSlowNucleonModel(0)
  72 {
  73 // Default constructor
  74     fCollisionGeometry = 0;
  75 }
  76
  77 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
  78     :AliGenerator(npart),
  79      fCMS(14000.),
  80      fMomentum(0.),
  81      fBeta(0.),
  82      fPmax (10.),
  83      fCharge(1),
  84      fProtonDirection(1.),
  85      fTemperatureG(0.05),
  86      fBetaSourceG(0.05),
  87      fTemperatureB(0.005),
  88      fBetaSourceB(0.),
  89      fNgp(0),
  90      fNgn(0),
  91      fNbp(0),
  92      fNbn(0),
  93      fDebug(0),
  94      fDebugHist1(0),
  95      fDebugHist2(0),
  96      fThetaDistribution(),
  97      fCosThetaGrayHist(),
  98      fCosTheta(),
  99      fBeamCrossingAngle(0.),
 100      fBeamDivergence(0.),
 101      fBeamDivEvent(0.),
 102      fSlowNucleonModel(new AliSlowNucleonModel())
 103
 104 {
 105 // Constructor
 106     fName  = "SlowNucleons";
 107     fTitle = "Generator for gray particles in pA collisions";
 108     fCollisionGeometry = 0;
 109 }
 110
 111 //____________________________________________________________
 112 AliGenSlowNucleons::~AliGenSlowNucleons()
 113 {
 114 // Destructor
 115     delete  fSlowNucleonModel;
 116 }
 117
 118 void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
 119 // Set direction of the proton to change between pA (1) and Ap (-1)
 120   fProtonDirection = dir / TMath::Abs(dir);
 121 }
 122
 123 void AliGenSlowNucleons::Init()
 124 {
 125   //
 126   // Initialization
 127   //
 128     Double_t kMass  = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
 129     fMomentum = fCMS/2. * Float_t(fZTarget) / Float_t(fATarget);
 130     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
 131     //printf("  fMomentum %f    fBeta %1.10f\n",fMomentum, fBeta);
 132     if (fDebug) {
 133         fDebugHist1 = new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.);
 134         fDebugHist2 = new TH2F("DebugHist2", "b  vs N_slow", 100, 0., 100., 15, 0., 15.);
 135         fCosThetaGrayHist = new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.);
 136     }
 137     //
 138     // non-uniform cos(theta) distribution
 139     //
 140     if(fThetaDistribution != 0) {
 141         fCosTheta = new TF1("fCosTheta",
 142                             "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./3.14159265358979312))*exp(2.*x/3.14159265358979312)",
 143                             -1., 1.);
 144     }
 145
 146     printf("\n  AliGenSlowNucleons: applying crossing angle %f mrad to slow nucleons\n",fBeamCrossingAngle*1000.);
 147 }
 148
 149 void AliGenSlowNucleons::FinishRun()
 150 {
 151 // End of run action
 152 // Show histogram for debugging if requested.
 153     if (fDebug) {
 154         TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
 155         c->Divide(2,1);
 156         c->cd(1);
 157         fDebugHist1->Draw("colz");
 158         c->cd(2);
 159         fDebugHist2->Draw();
 160         c->cd(3);
 161         fCosThetaGrayHist->Draw();
 162     }
 163 }
 164
 165
 166 void AliGenSlowNucleons::Generate()
 167 {
 168   //
 169   // Generate one event
 170   //
 171   //
 172   // Communication with Gray Particle Model
 173   //
 174     if (fCollisionGeometry) {
 175         Float_t b   = fCollisionGeometry->ImpactParameter();
 176         //      Int_t  nn   = fCollisionGeometry->NN();
 177         //      Int_t  nwn  = fCollisionGeometry->NwN();
 178         //      Int_t  nnw  = fCollisionGeometry->NNw();
 179         //      Int_t  nwnw = fCollisionGeometry->NwNw();
 180
 181         //fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 182         fSlowNucleonModel->GetNumberOfSlowNucleons2(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
 183         if (fDebug) {
 184             //printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw);
 185             printf("Slow nucleons: %d grayp  %d grayn  %d blackp  %d blackn \n", fNgp, fNgn, fNbp, fNbn);
 186             fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NNw(), 1.);
 187             fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
 188
 189         }
 190     }
 191
 192    //
 193     Float_t p[3] = {0., 0., 0.}, theta=0;
 194     Float_t origin[3] = {0., 0., 0.};
 195     Float_t time = 0.;
 196     Float_t polar [3] = {0., 0., 0.};
 197     Int_t nt, i, j;
 198     Int_t kf;
 199
 200     // Extracting 1 value per event for the divergence angle
 201     Double_t rvec = gRandom->Gaus(0.0, 1.0);
 202     fBeamDivEvent = fBeamDivergence * TMath::Abs(rvec);
 203     printf("\n  AliGenSlowNucleons: applying beam divergence %f mrad to slow nucleons\n",fBeamDivEvent*1000.);
 204
 205     if(fVertexSmear == kPerEvent) {
 206         Vertex();
 207         for (j=0; j < 3; j++) origin[j] = fVertex[j];
 208         time = fTime;
 209     } // if kPerEvent
 210 //
 211 //  Gray protons
 212 //
 213     fCharge = 1;
 214     kf = kProton;
 215     for(i = 0; i < fNgp; i++) {
 216         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 217         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 218         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 219                  time, kPNoProcess, nt, 1.,-2);
 220         KeepTrack(nt);
 221         SetProcessID(nt,kGrayProcess);
 222     }
 223 //
 224 //  Gray neutrons
 225 //
 226     fCharge = 0;
 227     kf = kNeutron;
 228     for(i = 0; i < fNgn; i++) {
 229         GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
 230         if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
 231         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 232                  time, kPNoProcess, nt, 1.,-2);
 233         KeepTrack(nt);
 234         SetProcessID(nt,kGrayProcess);
 235     }
 236 //
 237 //  Black protons
 238 //
 239     fCharge = 1;
 240     kf = kProton;
 241     for(i = 0; i < fNbp; i++) {
 242         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 243         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 244                  time, kPNoProcess, nt, 1.,-1);
 245         KeepTrack(nt);
 246         SetProcessID(nt,kBlackProcess);
 247     }
 248 //
 249 //  Black neutrons
 250 //
 251     fCharge = 0;
 252     kf = kNeutron;
 253     for(i = 0; i < fNbn; i++) {
 254         GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
 255         PushTrack(fTrackIt, -1, kf, p, origin, polar,
 256                  time, kPNoProcess, nt, 1.,-1);
 257         KeepTrack(nt);
 258         SetProcessID(nt,kBlackProcess);
 259     }
 260 }
 261
 262
 263
 264
 265 void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 266         Double_t beta, Float_t* q, Float_t &theta)
 267
 268 {
 269 /*
 270    Emit a slow nucleon with "temperature" T [GeV],
 271    from a source moving with velocity beta
 272    Three-momentum [GeV/c] is given back in q[3]
 273 */
 274
 275  //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
 276
 277  Double_t m, pmax, p, f, phi;
 278  TDatabasePDG * pdg = TDatabasePDG::Instance();
 279  const Double_t kMassProton  = pdg->GetParticle(kProton) ->Mass();
 280  const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 281
 282  /* Select nucleon type */
 283  if(charge == 0) m = kMassNeutron;
 284  else m = kMassProton;
 285
 286  /* Momentum at maximum of Maxwell-distribution */
 287
 288  pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
 289
 290  /* Try until proper momentum                                  */
 291  /* for lack of primitive function of the Maxwell-distribution */
 292  /* a brute force trial-accept loop, normalized at pmax        */
 293
 294  do
 295  {
 296      p = Rndm() * fPmax;
 297      f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 298  }
 299  while(f < Rndm());
 300
 301  /* Spherical symmetric emission for black particles (beta=0)*/
 302  if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * Rndm() - 1.);
 303  /* cos theta distributed according to experimental results for gray particles (beta=0.05)*/
 304  else if(fThetaDistribution!=0){
 305    Double_t costheta = fCosTheta->GetRandom();
 306    theta = TMath::ACos(costheta);
 307  }
 308  //
 309  phi   = 2. * TMath::Pi() * Rndm();
 310
 311
 312  /* Determine momentum components in system of the moving source */
 313  q[0] = p * TMath::Sin(theta) * TMath::Cos(phi);
 314  q[1] = p * TMath::Sin(theta) * TMath::Sin(phi);
 315  q[2] = p * TMath::Cos(theta);
 316  //if(fDebug==1) printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
 317
 318
 319  /* Transform to system of the target nucleus                             */
 320  /* beta is passed as negative, because the gray nucleons are slowed down */
 321  Lorentz(m, -beta, q);
 322  //if(fDebug==1) printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
 323
 324  /* Transform to laboratory system */
 325  Lorentz(m, fBeta, q);
 326  q[2] *= fProtonDirection;
 327  if(fDebug==1)printf("\n Momentum after LHC boost: p = (%f, %f, %f)\n",q[0],q[1],q[2]);
 328
 329  if(fBeamCrossingAngle>0.) BeamCrossDivergence(1, q); // applying crossing angle
 330  if(fBeamDivergence>0.) BeamCrossDivergence(2, q);    // applying divergence
 331
 332 }
 333
 334 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
 335 {
 336 /* Relativistic Maxwell-distribution */
 337     Double_t ekin;
 338     ekin = TMath::Sqrt(p*p+m*m)-m;
 339     return (p*p * exp(-ekin/T));
 340 }
 341
 342
 343 //_____________________________________________________________________________
 344 void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
 345 {
 346 /* Lorentz transform in the direction of q[2] */
 347
 348     Double_t gamma  = 1./TMath::Sqrt((1.-beta)*(1.+beta));
 349     Double_t energy = TMath::Sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
 350     q[2] = gamma * (q[2] + beta*energy);
 351     //printf(" \t beta %1.10f gamma %f energy %f -> p_z = %f\n",beta, gamma, energy,q[2]);
 352 }
 353
 354 //_____________________________________________________________________________
 355 void AliGenSlowNucleons::BeamCrossDivergence(Int_t iwhat, Float_t *pLab)
 356 {
 357   // Applying beam divergence and crossing angle
 358   //
 359   Double_t pmod = TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]+pLab[2]*pLab[2]);
 360
 361   Double_t tetdiv = 0.;
 362   Double_t fidiv = 0.;
 363   if(iwhat==1){
 364     tetdiv = fBeamCrossingAngle;
 365     fidiv = k2PI/4.;
 366   }
 367   else if(iwhat==2){
 368     tetdiv = fBeamDivEvent;
 369     fidiv = (gRandom->Rndm())*k2PI;
 370   }
 371
 372   Double_t tetpart = TMath::ATan2(TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]), pLab[2]);
 373   Double_t fipart=0.;
 374   if(TMath::Abs(pLab[1])>0. || TMath::Abs(pLab[0])>0.) fipart = TMath::ATan2(pLab[1], pLab[0]);
 375   if(fipart<0.) {fipart = fipart+k2PI;}
 376   tetdiv = tetdiv*kRaddeg;
 377   fidiv = fidiv*kRaddeg;
 378   tetpart = tetpart*kRaddeg;
 379   fipart = fipart*kRaddeg;
 380
 381   Double_t angleSum[2]={0., 0.};
 382   AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum);
 383
 384   Double_t tetsum = angleSum[0];
 385   Double_t fisum  = angleSum[1];
 386   //printf("tetpart %f fipart %f tetdiv %f fidiv %f angleSum %f %f\n",tetpart,fipart,tetdiv,fidiv,angleSum[0],angleSum[1]);
 387   tetsum = tetsum*kDegrad;
 388   fisum = fisum*kDegrad;
 389
 390   pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum);
 391   pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum);
 392   pLab[2] = pmod*TMath::Cos(tetsum);
 393   if(fDebug==1){
 394     if(iwhat==1) printf(" Beam crossing angle %f mrad ", fBeamCrossingAngle*1000.);
 395     else if(iwhat==2) printf(" Beam divergence %f mrad ", fBeamDivEvent*1000.);
 396     printf("  p = (%f, %f, %f)\n",pLab[0],pLab[1],pLab[2]);
 397   }
 398 }
 399
 400 //_____________________________________________________________________________
 401 void  AliGenSlowNucleons::AddAngle(Double_t theta1, Double_t phi1,
 402         Double_t theta2, Double_t phi2, Double_t *angleSum)
 403 {
 404   // Calculating the sum of 2 angles
 405   Double_t temp = -1.;
 406   Double_t conv = 180./TMath::ACos(temp);
 407
 408   Double_t ct1 = TMath::Cos(theta1/conv);
 409   Double_t st1 = TMath::Sin(theta1/conv);
 410   Double_t cp1 = TMath::Cos(phi1/conv);
 411   Double_t sp1 = TMath::Sin(phi1/conv);
 412   Double_t ct2 = TMath::Cos(theta2/conv);
 413   Double_t st2 = TMath::Sin(theta2/conv);
 414   Double_t cp2 = TMath::Cos(phi2/conv);
 415   Double_t sp2 = TMath::Sin(phi2/conv);
 416   Double_t cx = ct1*cp1*st2*cp2+st1*cp1*ct2-sp1*st2*sp2;
 417   Double_t cy = ct1*sp1*st2*cp2+st1*sp1*ct2+cp1*st2*sp2;
 418   Double_t cz = ct1*ct2-st1*st2*cp2;
 419
 420   Double_t rtetsum = TMath::ACos(cz);
 421   Double_t tetsum = conv*rtetsum;
 422   if(TMath::Abs(tetsum)<1e-4 || tetsum==180.) return;
 423
 424   temp = cx/TMath::Sin(rtetsum);
 425   if(temp>1.) temp=1.;
 426   if(temp<-1.) temp=-1.;
 427   Double_t fisum = conv*TMath::ACos(temp);
 428   if(cy<0) {fisum = 360.-fisum;}
 429
 430   angleSum[0] = tetsum;
 431   angleSum[1] = fisum;
 432 }
 433
 434 //_____________________________________________________________________________
 435 void AliGenSlowNucleons::SetProcessID(Int_t nt, UInt_t process)
 436 {
 437   // Tag the particle as
 438   // gray or black
 439   if (fStack)
 440     fStack->Particle(nt)->SetUniqueID(process);
 441   else
 442     gAlice->GetMCApp()->Particle(nt)->SetUniqueID(process);
 443 }