ZDC/AliGenZDC.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 /*
  17 $Log$
  18 Revision 1.12  2002/11/21 20:05:22  alibrary
  19 Removing AliMC and AliMCProcess
  20
  21 Revision 1.11  2001/11/11 17:27:50  hristov
  22 Minor corrections
  23
  24 Revision 1.10  2001/09/26 16:00:47  coppedis
  25 Minor change
  26
  27 Revision 1.8  2001/04/20 10:10:25  coppedis
  28 Minor changes
  29
  30 Revision 1.7  2001/03/15 16:13:28  coppedis
  31 Code review
  32
  33 Revision 1.6  2000/11/30 17:16:14  coppedis
  34 Changes suggested by fca
  35
  36 Revision 1.5  2000/11/22 11:30:12  coppedis
  37 Major code revision
  38
  39 Revision 1.4  2000/10/05 08:02:47  fca
  40 Correction of the generator direction
  41
  42 Revision 1.3  2000/10/02 21:28:20  fca
  43 Removal of useless dependecies via forward declarations
  44
  45 Revision 1.2  2000/07/11 11:12:34  fca
  46 Some syntax corrections for non standard HP aCC
  47
  48 Revision 1.1  2000/07/10 13:58:01  fca
  49 New version of ZDC from E.Scomparin & C.Oppedisano
  50
  51 Revision 1.7  2000/01/19 17:17:40  fca
  52
  53 Revision 1.6  1999/09/29 09:24:35  fca
  54 Introduction of the Copyright and cvs Log
  55
  56 */
  57 #include <assert.h>
  58
  59 #include <TDatabasePDG.h>
  60 #include <TLorentzVector.h>
  61 #include <TMCProcess.h>
  62 #include <TPDGCode.h>
  63 #include <TRandom.h>
  64 #include <TVector3.h>
  65
  66 #include "AliConst.h"
  67 #include "AliGenZDC.h"
  68 #include "AliRun.h"
  69
  70 ClassImp(AliGenZDC)
  71
  72 //_____________________________________________________________________________
  73 AliGenZDC::AliGenZDC()
  74    :AliGenerator()
  75 {
  76   //
  77   // Default constructor
  78   //
  79   fIpart = 0;
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliGenZDC::AliGenZDC(Int_t npart)
  84    :AliGenerator(npart)
  85 {
  86   //
  87   // Standard constructor
  88   //
  89   fName = "AliGenZDC";
  90   fTitle = "Generation of Test Particles for ZDCs";
  91   fIpart = kNeutron;
  92   fCosx  = 0.;
  93   fCosy  = 0.;
  94   fCosz  = 1.;
  95   fPseudoRapidity = 0.;
  96
  97   fFermiflag = 1;
  98   // LHC values for beam divergence and crossing angle
  99   fBeamDiv = 0.000032;
 100   fBeamCrossAngle = 0.0001;
 101   fBeamCrossPlane = 2;
 102
 103   Int_t i, j;
 104   for(i=0; i<201; i++){
 105      fProbintp[i] = 0;
 106      fProbintn[i] = 0;
 107   }
 108   for(j=0; j<3; j++){
 109      fPp[i] = 0;
 110   }
 111   fDebugOpt = 0;
 112 }
 113
 114 //_____________________________________________________________________________
 115 void AliGenZDC::Init()
 116 {
 117   printf("\n\n          AliGenZDC initialized with:\n");
 118   printf("   Fermi flag = %d, Beam Divergence = %f, Crossing Angle "
 119          "= %f, Crossing Plane = %d\n\n", fFermiflag, fBeamDiv, fBeamCrossAngle,
 120          fBeamCrossPlane);
 121
 122   //Initialize Fermi momentum distributions for Pb-Pb
 123   FermiTwoGaussian(207.,82.,fPp,fProbintp,fProbintn);
 124 }
 125
 126 //_____________________________________________________________________________
 127 void AliGenZDC::Generate()
 128 {
 129   //
 130   // Generate one trigger (n or p)
 131   //
 132   Int_t i;
 133
 134   Double_t Mass, pLab[3], fP0, fP[3], fBoostP[3], ddp[3], dddp0, dddp[3];
 135   Float_t  fPTrack[3], ptot = fPMin;
 136   Int_t nt;
 137
 138   if(fPseudoRapidity==0.){
 139     pLab[0] = ptot*fCosx;
 140     pLab[1] = ptot*fCosy;
 141     pLab[2] = ptot*fCosz;
 142   }
 143   else{
 144     Float_t scang = 2*TMath::ATan(TMath::Exp(-(fPseudoRapidity)));
 145     pLab[0] = -ptot*TMath::Sin(scang);
 146     pLab[1] = 0.;
 147     pLab[2] = ptot*TMath::Cos(scang);
 148   }
 149   for(i=0; i<=2; i++){
 150      fP[i] = pLab[i];
 151   }
 152
 153
 154   // Beam divergence and crossing angle
 155   if(fBeamCrossAngle!=0.) {
 156     BeamDivCross(1,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);
 157     for(i=0; i<=2; i++){
 158        fP[i] = pLab[i];
 159     }
 160   }
 161   if(fBeamDiv!=0.) {
 162     BeamDivCross(0,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);
 163     for(i=0; i<=2; i++){
 164        fP[i] = pLab[i];
 165     }
 166   }
 167
 168   // If required apply the Fermi momentum
 169   if(fFermiflag==1){
 170     if((fIpart==kProton) || (fIpart==kNeutron)){
 171       ExtractFermi(fIpart,fPp,fProbintp,fProbintn,ddp);
 172     }
 173     Mass=gAlice->PDGDB()->GetParticle(fIpart)->Mass();
 174     fP0 = TMath::Sqrt(fP[0]*fP[0]+fP[1]*fP[1]+fP[2]*fP[2]+Mass*Mass);
 175     for(i=0; i<=2; i++){
 176        dddp[i] = ddp[i];
 177     }
 178     dddp0 = TMath::Sqrt(dddp[0]*dddp[0]+dddp[1]*dddp[1]+dddp[2]*dddp[2]+Mass*Mass);
 179
 180     TVector3 b(fP[0]/fP0, fP[1]/fP0, fP[2]/fP0);
 181     TLorentzVector pFermi(dddp[0], dddp[1], dddp[2], dddp0);
 182
 183     pFermi.Boost(b);
 184     for(i=0; i<=2; i++){
 185        fBoostP[i] = pFermi[i];
 186        fP[i] = pFermi[i];
 187     }
 188
 189   }
 190
 191   for(i=0; i<=2; i++){
 192      fPTrack[i] = fP[i];
 193   }
 194
 195   Float_t polar[3] = {0,0,0};
 196   gAlice->SetTrack(fTrackIt,-1,fIpart,fPTrack,fOrigin.GetArray(),polar,0,
 197                    kPPrimary,nt);
 198   if(fDebugOpt == 1){
 199     printf("\n\n                Track momentum:\n");
 200     printf("\n   fPTrack = %f, %f, %f \n",fPTrack[0],fPTrack[1],fPTrack[2]);
 201   }
 202 }
 203
 204 //_____________________________________________________________________________
 205 void AliGenZDC::FermiTwoGaussian(Float_t A, Float_t Z, Double_t *fPp,
 206                 Double_t *fProbintp, Double_t *fProbintn)
 207 {
 208 //
 209 // Momenta distributions according to the "double-gaussian"
 210 // distribution (Ilinov) - equal for protons and neutrons
 211 //
 212
 213    fProbintp[0] = 0;
 214    fProbintn[0] = 0;
 215    Double_t sig1 = 0.113;
 216    Double_t sig2 = 0.250;
 217    Double_t alfa = 0.18*(TMath::Power((A/12.),(Float_t)1/3));
 218    Double_t xk = (2*k2PI)/((1.+alfa)*(TMath::Power(k2PI,1.5)));
 219
 220    for(Int_t i=1; i<=200; i++){
 221       Double_t p = i*0.005;
 222       fPp[i] = p;
 223       Double_t e1 = (p*p)/(2.*sig1*sig1);
 224       Double_t e2 = (p*p)/(2.*sig2*sig2);
 225       Double_t f1 = TMath::Exp(-(e1));
 226       Double_t f2 = TMath::Exp(-(e2));
 227       Double_t probp = xk*p*p*(f1/(TMath::Power(sig1,3.))+
 228                       alfa*f2/(TMath::Power(sig2,3.)))*0.005;
 229       fProbintp[i] = fProbintp[i-1] + probp;
 230       fProbintn[i] = fProbintp[i];
 231    }
 232    if(fDebugOpt == 1){
 233      printf("\n\n               Initialization of Fermi momenta distribution \n");
 234      for(Int_t i=0; i<=200; i++){
 235         printf(" fProbintp[%d] = %f, fProbintn[%d] = %f\n",i,fProbintp[i],i,fProbintn[i]);
 236      }
 237    }
 238 }
 239 //_____________________________________________________________________________
 240 void AliGenZDC::ExtractFermi(Int_t id, Double_t *fPp, Double_t *fProbintp,
 241                 Double_t *fProbintn, Double_t *ddp)
 242 {
 243 //
 244 // Compute Fermi momentum for spectator nucleons
 245 //
 246
 247   Int_t i;
 248   Float_t xx = gRandom->Rndm();
 249   assert ( id==kProton || id==kNeutron );
 250   if(id==kProton){
 251     for(i=1; i<=200; i++){
 252        if((xx>=fProbintp[i-1]) && (xx<fProbintp[i])) break;
 253        }
 254   }
 255   else {
 256     for(i=0; i<=200; i++){
 257        if((xx>=fProbintn[i-1]) && (xx<fProbintn[i])) break;
 258        }
 259    }
 260          Float_t pext = fPp[i]+0.001;
 261          Float_t phi = k2PI*(gRandom->Rndm());
 262          Float_t cost = (1.-2.*(gRandom->Rndm()));
 263          Float_t tet = TMath::ACos(cost);
 264          ddp[0] = pext*TMath::Sin(tet)*TMath::Cos(phi);
 265          ddp[1] = pext*TMath::Sin(tet)*TMath::Sin(phi);
 266          ddp[2] = pext*cost;
 267
 268   if(fDebugOpt == 1){
 269     printf("\n\n                Extraction of Fermi momentum\n");
 270     printf("\n  pxFermi = %f  pyFermi = %f  pzFermi = %f \n",ddp[0],ddp[1],ddp[2]);
 271   }
 272 }
 273
 274 //_____________________________________________________________________________
 275 void AliGenZDC::BeamDivCross(Int_t icross, Float_t fBeamDiv, Float_t fBeamCrossAngle,
 276                 Int_t fBeamCrossPlane, Double_t *pLab)
 277 {
 278   Double_t tetpart, fipart, tetdiv=0, fidiv=0, angleSum[2], tetsum, fisum;
 279   Double_t rvec;
 280
 281   Int_t i;
 282   Double_t pmq = 0.;
 283   for(i=0; i<=2; i++){
 284      pmq = pmq+pLab[i]*pLab[i];
 285   }
 286   Double_t pmod = TMath::Sqrt(pmq);
 287
 288   if(icross==0){
 289     rvec = gRandom->Gaus(0.0,1.0);
 290     tetdiv = fBeamDiv * TMath::Abs(rvec);
 291     fidiv = (gRandom->Rndm())*k2PI;
 292   }
 293   else if(icross==1){
 294     if(fBeamCrossPlane==0.){
 295       tetdiv = 0.;
 296       fidiv = 0.;
 297     }
 298     else if(fBeamCrossPlane==1.){
 299       tetdiv = fBeamCrossAngle;
 300       fidiv = 0.;
 301     }
 302     else if(fBeamCrossPlane==2.){
 303       tetdiv = fBeamCrossAngle;
 304       fidiv = k2PI/4.;
 305     }
 306   }
 307
 308   tetpart = TMath::ATan(TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1])/pLab[2]);
 309   if(pLab[1]!=0. || pLab[0]!=0.){
 310     fipart = TMath::ATan2(pLab[1],pLab[0]);
 311   }
 312   else{
 313     fipart = 0.;
 314   }
 315   if(fipart<0.) {fipart = fipart+k2PI;}
 316   tetdiv = tetdiv*kRaddeg;
 317   fidiv = fidiv*kRaddeg;
 318   tetpart = tetpart*kRaddeg;
 319   fipart = fipart*kRaddeg;
 320   AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum);
 321   tetsum = angleSum[0];
 322   fisum  = angleSum[1];
 323   tetsum = tetsum*kDegrad;
 324   fisum = fisum*kDegrad;
 325   pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum);
 326   pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum);
 327   pLab[2] = pmod*TMath::Cos(tetsum);
 328   if(fDebugOpt == 1){
 329     printf("\n\n                Beam divergence and crossing angle\n");
 330     for(i=0; i<=2; i++){
 331        printf("         pLab[%d] = %f\n",i,pLab[i]);
 332     }
 333   }
 334 }
 335
 336 //_____________________________________________________________________________
 337 void  AliGenZDC::AddAngle(Double_t theta1, Double_t phi1, Double_t theta2,
 338                Double_t phi2, Double_t *angleSum)
 339 {
 340   Double_t temp, conv, cx, cy, cz, ct1, st1, ct2, st2, cp1, sp1, cp2, sp2;
 341   Double_t rtetsum, tetsum, fisum;
 342
 343   temp = -1.;
 344   conv = 180./TMath::ACos(temp);
 345
 346   ct1 = TMath::Cos(theta1/conv);
 347   st1 = TMath::Sin(theta1/conv);
 348   cp1 = TMath::Cos(phi1/conv);
 349   sp1 = TMath::Sin(phi1/conv);
 350   ct2 = TMath::Cos(theta2/conv);
 351   st2 = TMath::Sin(theta2/conv);
 352   cp2 = TMath::Cos(phi2/conv);
 353   sp2 = TMath::Sin(phi2/conv);
 354   cx = ct1*cp1*st2*cp2+st1*cp1*ct2-sp1*st2*sp2;
 355   cy = ct1*sp1*st2*cp2+st1*sp1*ct2+cp1*st2*sp2;
 356   cz = ct1*ct2-st1*st2*cp2;
 357
 358   rtetsum = TMath::ACos(cz);
 359   tetsum = conv*rtetsum;
 360   if(tetsum==0. || tetsum==180.){
 361     fisum = 0.;
 362     return;
 363   }
 364   temp = cx/TMath::Sin(rtetsum);
 365   if(temp>1.) temp=1.;
 366   if(temp<-1.) temp=-1.;
 367   fisum = conv*TMath::ACos(temp);
 368   if(cy<0) {fisum = 360.-fisum;}
 369   angleSum[0] = tetsum;
 370   angleSum[1] = fisum;
 371 }
 372