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