TPHIC/AliGenTPHIC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, 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 /* $Id$ */
  18
  19 // Event generator of two-photon processes
  20 // in ultra-peripheral ion collisions.
  21 // 5 two-photon process are implemented, see comments to SetProcess().
  22 //%
  23 // The example of the generator initialization for the process
  24 // gamma gamma -> X in CaCa collisions at 7000 TeV/nucleon is the following:
  25 //%
  26 //    AliGenTPHIC *gener = new AliGenTPHIC();
  27 //    gener->SetProcess(1);
  28 //    gener->SetBeamEnergy(3500.);
  29 //    gener->SetBeamZ(20);
  30 //    gener->SetBeamA(40);
  31 //    gener->SetMggRange(70.,200.);
  32 //    gener->SetYggRange(-5.,5.);
  33 //    gener->SetLumFunName("lum_ca_70_200.dat");
  34 //    gener->SetLumFunFlag(-1);
  35 //    gener->Init();
  36 //%
  37 // The two-photon luminosity function needed for the process cross section
  38 // calculation is time-consuming, therefore it can be calculated once for a
  39 // selected two-photon energy and rapidity range and selected ion species,
  40 // saved into a file and then can be reused in further calculation.
  41 //%
  42 // The integral cross section of the process is calculated in each event
  43 // by the Monte Carlo integration procedure, the differential cross section
  44 // of each thrown event is assigned as a weight to each track of the event
  45 // which can be then used during the event analysis by retreiving this weight
  46 // from any of the event tracks.
  47 //
  48 // The manual of the fortran generator is available in
  49 // $ALICE_ROOT/TPHIC/TPHIC_doc.ps.gz
  50 // For the two-photon physics in heavy ion collisions see
  51 // G.Baur et al, Phys.Rep. 364 (2002), 359.
  52 //%
  53 // Author: Yuri.Kharlov@cern.ch
  54 // 15 April 2003
  55
  56 #include <TParticle.h>
  57 #include <TParticlePDG.h>
  58 #include <TDatabasePDG.h>
  59 #include <TClonesArray.h>
  60
  61 #include "AliPythia.h"
  62 #include "AliRun.h"
  63 #include <AliGenTPHIC.h>
  64 #include <TPHICgen.h>
  65 #include <TPHICcommon.h>
  66
  67 ClassImp(AliGenTPHIC)
  68
  69 //------------------------------------------------------------
  70
  71 AliGenTPHIC::AliGenTPHIC() :
  72   AliGenMC(),
  73   fTPHICgen(0x0),
  74   fPythia(0x0),
  75   fParticles(0x0),
  76   fEvent(-1),
  77   fDebug(0),
  78   fDebugEventFirst(-1),
  79   fDebugEventLast(-1)
  80 {
  81   // Constructor: create generator instance,
  82   // create particle array
  83   // Set TPHIC parameters to default values:
  84   // eta_b production in Ca-Ca collisions at 7 A*TeV
  85
  86   SetMC(new TPHICgen());
  87   fTPHICgen = (TPHICgen*) fMCEvGen;
  88   fPythia = AliPythia::Instance();
  89   fParticles = new TClonesArray("TParticle",100);
  90
  91   SetProcess   ();
  92   SetBeamEnergy();
  93   SetBeamZ     ();
  94   SetBeamA     ();
  95   SetYggRange  ();
  96   SetMggRange  ();
  97   SetNgridY    ();
  98   SetNgridM    ();
  99   SetLumFunName();
 100   SetLumFunFlag();
 101   SetKfOnium   ();
 102 }
 103
 104 //____________________________________________________________
 105 AliGenTPHIC::~AliGenTPHIC()
 106 {
 107   // Destroys the object, deletes and disposes all TParticles currently on list.
 108   if (fParticles) {
 109     fParticles->Delete();
 110     delete fParticles;
 111     fParticles = 0;
 112   }
 113 }
 114
 115 //____________________________________________________________
 116 void AliGenTPHIC::Init()
 117 {
 118   // Initialize the generator TPHIC
 119
 120   fTPHICgen->Initialize();
 121   fEvent = 0;
 122 }
 123
 124 //____________________________________________________________
 125 void AliGenTPHIC::Generate()
 126 {
 127   // Generate one event of two-photon process.
 128   // Gaussian smearing on the vertex is done if selected.
 129   // All particles from the TPHIC/PYTHIA event listing
 130   // are stored in TreeK, and only final-state particles are tracked.
 131   // The event differectial cross section is assigned as a weight
 132   // to each track of the event.
 133
 134   Float_t polar[3]   = {0,0,0};
 135   Float_t origin0[3] = {0,0,0};
 136   Float_t time0 = 0.;
 137   Float_t origin[3]  = {0,0,0};
 138   Float_t p[3], tof;
 139   Double_t weight;
 140
 141   Int_t    ks,kf,iparent,nt, trackIt;
 142   Float_t  random[6];
 143   const Float_t kconv=0.001/2.999792458e8;
 144
 145   fTPHICgen->GenerateEvent();
 146   weight = GetXSectionCurrent();
 147   if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 148       gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
 149   fPythia->ImportParticles(fParticles,"All");
 150
 151   if (fDebug == 1)
 152     Info("Generate()","one event is produced");
 153
 154   Int_t j;
 155   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 156   time0 = fTimeOrigin;
 157   if(fVertexSmear==kPerEvent) {
 158     Rndm(random,6);
 159     for (j=0;j<3;j++) {
 160       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 161         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 162     }
 163     Rndm(random,2);
 164     time0 += fOsigma[2]/TMath::Ccgs()*
 165       TMath::Cos(2*random[0]*TMath::Pi())*
 166       TMath::Sqrt(-2*TMath::Log(random[1]));
 167   }
 168
 169   Int_t ip;
 170   Int_t np = fParticles->GetEntriesFast();
 171   TParticle *iparticle;
 172 //   Int_t* pParent   = new Int_t[np];
 173   for (ip=0; ip<np; ip++) {
 174     iparticle = (TParticle *) fParticles->At(ip);
 175     ks = iparticle->GetStatusCode();
 176 //     // No initial state partons
 177 //     if (ks==21) continue;
 178     p[0] = iparticle->Px();
 179     p[1] = iparticle->Py();
 180     p[2] = iparticle->Pz();
 181     origin[0] = origin0[0]+iparticle->Vx()/10.;
 182     origin[1] = origin0[1]+iparticle->Vy()/10.;
 183     origin[2] = origin0[2]+iparticle->Vz()/10.;
 184     kf = CheckPDGCode(iparticle->GetPdgCode());
 185     iparent = -1;
 186     tof     = time0 + kconv*iparticle->T();
 187     if (ks == 1) trackIt = 1;
 188     else         trackIt = 0;
 189     PushTrack(fTrackIt*trackIt,iparent,kf,p,origin,polar,tof,kPPrimary,nt,weight,ks);
 190     KeepTrack(nt);
 191
 192     if (fDebug == 2)
 193       printf("id=%+4d, parent=%3d, ks=%d, p = (%+11.4e,%+11.4e,%+11.4e) GeV\n",
 194              kf,iparent,fTrackIt*trackIt,p[0],p[1],p[2]);
 195   }
 196   fEvent++;
 197   fTPHICgen->SetNEVENT(fEvent);
 198   if (fDebug == 1 && fEvent%100 == 0) {
 199     Info("Generate","Event %d\n",fEvent);
 200     fTPHICgen->Finish();
 201   }
 202 }
 203
 204 //____________________________________________________________
 205 void AliGenTPHIC::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 206 {
 207   // Set a range of event numbers, for which a table
 208   // of generated particle will be printed
 209   fDebugEventFirst = eventFirst;
 210   fDebugEventLast  = eventLast;
 211   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 212 }
 213
 214 //____________________________________________________________
 215 void AliGenTPHIC::SetProcess       (Int_t   proc  )
 216 {
 217   // Set process number:
 218   // proc=1 - gamma gamma -> X
 219   // proc=2 - gamma gamma -> quarkonium
 220   // proc=3 - gamma gamma -> fermion+ fermion-
 221   // proc=4 - gamma gamma -> W+ W-
 222   // proc=5 - not implemented
 223   // proc=6 - gamma gamma -> V1 V2 (vector meson pair)
 224
 225   fTPHICgen->SetIPROC(proc);
 226 }
 227 //____________________________________________________________
 228   void AliGenTPHIC::SetBeamEnergy  (Float_t energy)
 229 {
 230   // Set energy of the beam ion per nucleon in GeV
 231   fTPHICgen->SetEBMN(energy);
 232 }
 233 //____________________________________________________________
 234   void AliGenTPHIC::SetBeamZ       (Int_t   z     )
 235 {
 236   // Set beam ion charge
 237   fTPHICgen->SetIZ(z);
 238 }
 239 //____________________________________________________________
 240   void AliGenTPHIC::SetBeamA       (Int_t   a     )
 241 {
 242   // Set beam ion atomic number
 243   fTPHICgen->SetIA(a);
 244 }
 245 //____________________________________________________________
 246   void AliGenTPHIC::SetYggRange    (Float_t ymin, Float_t ymax)
 247 {
 248   // Set rapidity range of 2-photon system for the
 249   // luminosity function calculation
 250   fTPHICgen->SetYMIN(ymin);
 251   fTPHICgen->SetYMAX(ymax);
 252 }
 253 //____________________________________________________________
 254   void AliGenTPHIC::SetMggRange    (Float_t mmin, Float_t mmax)
 255 {
 256   // Set invariant mass range of 2-photon system for the
 257   // luminosity function calculation
 258   fTPHICgen->SetAMIN(mmin);
 259   fTPHICgen->SetAMAX(mmax);
 260 }
 261 //____________________________________________________________
 262   void AliGenTPHIC::SetNgridY      (Int_t   ny    )
 263 {
 264   // Set number of nodes on the grid along the rapidity axis
 265   // to calculate the 2-photon luminosity function
 266   fTPHICgen->SetNY(ny);
 267 }
 268 //____________________________________________________________
 269   void AliGenTPHIC::SetNgridM      (Int_t   nm    )
 270 {
 271   // Set number of nodes on the grid along the mass axis
 272   // to calculate the 2-photon luminosity function
 273   fTPHICgen->SetNMAS(nm);
 274 }
 275 //____________________________________________________________
 276   void AliGenTPHIC::SetLumFunName  (TString name  )
 277 {
 278   // Set filename to store the 2-photon luminosity
 279   // function calculated on the grid
 280   fTPHICgen->SetLUMFIL(name);
 281 }
 282 //____________________________________________________________
 283   void AliGenTPHIC::SetLumFunFlag  (Int_t   flag  )
 284 {
 285   // Set flag to calculate the 2-photon luminosity function:
 286   // flag=-1 if a new lumimosity function to be calculated
 287   //         and stored in a file
 288   // flag=+1 if a previously calculated function to be read
 289   //         from a file
 290   fTPHICgen->SetILUMF(flag);
 291 }
 292 //____________________________________________________________
 293   void AliGenTPHIC::SetKfFermion   (Int_t   kf    )
 294 {
 295   // Set a PDG flavour code of a fermion for the process 3,
 296   // gamma gamma -> fermion+ fermion-
 297   fTPHICgen->SetKFERM(kf);
 298 }
 299 //____________________________________________________________
 300   void AliGenTPHIC::SetKfOnium    (Int_t   kf    )
 301 {
 302   // Set a PDG flavour code of a quarkonium for the process 2,
 303   // gamma gamma -> quarkonium
 304   fTPHICgen->SetKFONIUM(kf);
 305 }
 306 //____________________________________________________________
 307   void AliGenTPHIC::SetMassOnium   (Float_t mass  )
 308 {
 309   // Set a quarkonium mass [GeV] for the process 2 if it
 310   // differes from the Pythia's particle table.
 311   // For the well-known quarkonia no need to set this mass
 312   // because it will be taken from the Pythia table
 313   fTPHICgen->SetXMRES(mass);
 314 }
 315 //____________________________________________________________
 316   void AliGenTPHIC::SetGGwidthOnium(Float_t width )
 317 {
 318   // Set 2-photon partial width [GeV] of the quarkonium for
 319   // the process 2 if it differes fromthe Pythia's particle table.
 320   // For the well-known quarkonia no need to set this width
 321   // because it will be taken from the Pythia table
 322   fTPHICgen->SetXGGRES(width);
 323 }
 324 //____________________________________________________________
 325   void AliGenTPHIC::SetKfVmesons   (Int_t kf1, Int_t kf2)
 326 {
 327   // Set PDG flavour codes of the two vector vesons
 328   // for the process 2:  gamma gamma -> V1 V2
 329   // So far this processes is implemented for the following
 330   // mesons and their combinations only:
 331   // pho0 (113), omega (223), phi (333), J/psi (443)
 332   fTPHICgen->SetKV1(kf1);
 333   fTPHICgen->SetKV2(kf2);
 334 }
 335 //____________________________________________________________
 336   Float_t AliGenTPHIC::GetGGmass    ()
 337 {
 338   // Get invariant mass of generated 2-photon system
 339   return fTPHICgen->GetWSQ();
 340 }
 341 //____________________________________________________________
 342   Float_t AliGenTPHIC::GetGGrapidity()
 343 {
 344   // Get rapidity of generated 2-photon system
 345   return fTPHICgen->GetYGG();
 346 }
 347 //____________________________________________________________
 348   Float_t AliGenTPHIC::GetG1mass    ()
 349 {
 350   // Get a mass of the first virtual photon of
 351   // the 2-photon process (-sqrt(q1^2)).
 352   return fTPHICgen->GetXMG1();
 353 }
 354 //____________________________________________________________
 355   Float_t AliGenTPHIC::GetG2mass    ()
 356 {
 357   // Get a mass of the second virtual photon of
 358   // the 2-photon process (-sqrt(q2^2)).
 359   return fTPHICgen->GetXMG2();
 360 }
 361 //____________________________________________________________
 362   TClonesArray*  AliGenTPHIC::GetParticleList ()
 363 {
 364   // Get array of particles of the event listing
 365   return fParticles;
 366 }
 367 //____________________________________________________________
 368   TLorentzVector AliGenTPHIC::MomentumRecNucl1()
 369 {
 370   // Get 4-momentum of the first recoil nucleus after
 371   // the 2-photon process.
 372   return TLorentzVector(fTPHICgen->GetPTAG1(0),
 373                         fTPHICgen->GetPTAG1(1),
 374                         fTPHICgen->GetPTAG1(2),
 375                         fTPHICgen->GetPTAG1(3));
 376 }
 377 //____________________________________________________________
 378   TLorentzVector AliGenTPHIC::MomentumRecNucl2()
 379 {
 380   // Get 4-momentum of the first recoil nucleus after
 381   // the 2-photon process.
 382   return TLorentzVector(fTPHICgen->GetPTAG2(0),
 383                         fTPHICgen->GetPTAG2(1),
 384                         fTPHICgen->GetPTAG2(2),
 385                         fTPHICgen->GetPTAG2(3));
 386 }
 387 //____________________________________________________________
 388   Float_t AliGenTPHIC::GetXSectionCurrent()
 389 {
 390   // Get the cross section of the produced event for the
 391   // Monte Carlo integral cross section calculation
 392   return fTPHICgen->GetXSCUR();
 393 }
 394 //____________________________________________________________
 395   Float_t AliGenTPHIC::GetXSection       ()
 396 {
 397   // Get the integral cross section of the process
 398   // calculated so far
 399   return fTPHICgen->GetXSTOT();
 400 }
 401 //____________________________________________________________
 402   Float_t AliGenTPHIC::GetXSectionError  ()
 403 {
 404   // Get the error of the integral cross section of the process
 405   // calculated so far
 406   return fTPHICgen->GetXSTOTE();
 407 }