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