TEPEMGEN/AliGenQEDBg.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2002, 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  * Copyright(c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken          *
  17  * See $ALICE_ROOT/EpEmGen/diffcross.f for full Copyright notice          *
  18  *                                                                        *
  19  *                                                                        *
  20  * Copyright(c) 2002 Kai Hencken, Yuri Kharlov, Serguei Sadovsky          *
  21  * See $ALICE_ROOT/EpEmGen/epemgen.f for full Copyright notice            *
  22  *                                                                        *
  23  **************************************************************************/
  24
  25 /* $Id$ */
  26
  27 // Event generator for background from e+e- pair production in ultraperipheral PbPb collisions
  28 // at 5.5 TeV/nucleon, integrated over specific readout cycle
  29 // Derived from AliGenEpEmv1
  30 // Author: ruben.shahoyan@cern.ch
  31 //%
  32 // References:
  33 // [1] "Multiple electromagnetic electron positron pair production in
  34 //      relativistic heavy ion collisions".
  35 //      Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur,
  36 //      Phys. Rev. A55 (1997) 396.
  37 // [2] K.Hencken, Yu.Kharlov, S.Sadovsky, Internal ALICE Note 2002-27.
  38 //%
  39 // Usage:
  40 // Initialization:
  41 //    AliGenQEDBg *gener = new AliGenQEDBg();
  42 //    gener->SetXXXRange(); // Set kinematics range
  43 //    gener->SetLumiIntTime(double lumi, double sec); // luminosity and intergration time in seconds
  44 //    gener->Init();
  45 // Event generation:
  46 //    gener->Generate(); // Produce poissonian number of e+e- pair with average number
  47 //    corresponding to requested integration time with given beam luminosity
  48 //    Each pair has its own vertex, and time randomly distributed at originT
  49 //    and originT+ integration time
  50 //
  51 //    For details of pair generation see AliGenEpEmv1.cxx by Yuri.Kharlov@cern.ch
  52
  53 /*
  54   The most useful way of using it is in the overlay with other generators (using the cocktail):
  55   In the Config.C
  56   // load the library:
  57   gSystem->Load("libTEPEMGEN");
  58   //
  59   // add to AliGenCocktail as:
  60   AliGenCocktail *cocktail = new AliGenCocktail();
  61   // ... setup other stuff
  62   //
  63   // QED background
  64   AliGenQEDBg*  genBg = new AliGenQEDBg();
  65   genBg->SetEnergyCMS(5500);
  66   genBg->SetProjectile("A", 208, 82);
  67   genBg->SetTarget    ("A", 208, 82);
  68   genBg->SetYRange(-6.,3);
  69   genBg->SetPtRange(1.e-3,1.0);      // Set pt limits (GeV) for e+-: 1MeV corresponds to max R=13.3mm at 5kGaus
  70   genBg->SetLumiIntTime(6.e27,20e-6); // luminosity and integration time
  71   //
  72   cocktail->AddGenerator(genBg,"QEDep",1);
  73   // We need to generate independent vertex for every event, this must be set after adding to cocktail
  74   // Note that the IO origin and sigma's is transferred from AliGenCocktail to daughter generators
  75   genBg->SetVertexSource(kInternal);
  76 */
  77
  78 #include "AliLog.h"
  79 #include "AliGenQEDBg.h"
  80 #include <TParticle.h>
  81 #include <TParticlePDG.h>
  82 #include <TDatabasePDG.h>
  83 #include <TEpEmGen.h>
  84
  85 ClassImp(AliGenQEDBg)
  86
  87 //------------------------------------------------------------
  88
  89 AliGenQEDBg::AliGenQEDBg()
  90 :  fLumi(0)
  91   ,fXSection(0)
  92   ,fXSectionEps(1e-2)
  93   ,fIntTime(0)
  94   ,fPairsInt(-1)
  95   ,fMinXSTest(1e3)
  96   ,fMaxXSTest(1e7)
  97 {
  98 }
  99
 100 //____________________________________________________________
 101 AliGenQEDBg::~AliGenQEDBg()
 102 {
 103   // Destructor
 104 }
 105
 106 //____________________________________________________________
 107 void AliGenQEDBg::Init()
 108 {
 109   // Initialisation:
 110   AliInfo(Form("Will estimate QED bg. for L=%e cm^-2*s^-1 and Integration Time of %e s.",fLumi,fIntTime));
 111   if (fLumi<=0 || fIntTime<=0) {
 112     AliWarning("One of parameters is not set properly, no pairs will be generated");
 113     return;
 114   }
 115   //
 116   // initialize the generator of e+e- pair production
 117   AliGenEpEmv1::Init();
 118   //
 119   fPairsInt = 0;
 120   int ngen = 0;
 121   AliInfo(Form("Estimating x-section with min.relative precision of %f and min/max test: %d/%d",
 122                fXSectionEps,int(fMinXSTest),int(fMaxXSTest)));
 123   //
 124   double yElectron,yPositron,xElectron,xPositron,phi12,weight,err=0;
 125   fXSection = -1;
 126   do {
 127     fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
 128     if (++ngen>fMinXSTest) { // ensure min number of tests
 129       fXSection = fEpEmGen->GetXsection();
 130       err = fEpEmGen->GetDsection();
 131     }
 132   } while(!((fXSection>0 && err/fXSection<fXSectionEps) || ngen>fMaxXSTest));
 133   //
 134   if (fXSection<=0) {
 135     AliError(Form("X-section = %e after %d trials, cannot generate",fXSection,ngen));
 136     return;
 137   }
 138   fPairsInt = fXSection*1e-21*fLumi*fIntTime; // xsestion is in kbarn!
 139   AliInfo(Form("Estimated x-secion: %e+-%ekb in %d tests, <Npairs>=%e per %e time interval",
 140                fXSection,err,ngen,fPairsInt,fIntTime));
 141   //
 142 }
 143
 144 //____________________________________________________________
 145 void AliGenQEDBg::Generate()
 146 {
 147   //
 148   // Generate poissian <fPairsInt> e+e- pairs, each one with its vertex
 149   //
 150   Float_t polar[3]= {0,0,0};
 151   Float_t origin[3];
 152   Float_t time = 0.;
 153   Float_t p[3];
 154   //
 155   int npairs,nt,id;;
 156   if (fPairsInt<=0 || (npairs=gRandom->Poisson(fPairsInt))<1) return;
 157   //
 158   Double_t ptElectron,ptPositron, phiElectron,phiPositron, mt, ms2 = fMass*fMass;
 159   Double_t phi12=0,xElectron=0,xPositron=0,yElectron=0,yPositron=0,weight=0;
 160   //
 161   for (int i=0;i<npairs;i++) {
 162     // each pair has its own vertex and time
 163     Vertex();
 164     for (int j=0;j<3;j++) origin[j] = fVertex[j];
 165     time = fTimeOrigin+gRandom->Rndm()*fIntTime;
 166     //
 167     fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
 168     ptElectron  = TMath::Power(10,xElectron) * 1.e-03;;
 169     ptPositron  = TMath::Power(10,xPositron) * 1.e-03;;
 170     phiElectron = fPhiMin + gRandom->Rndm() * (fPhiMax-fPhiMin);
 171     phiPositron = phiElectron + phi12;
 172     // Produce electron
 173     mt = TMath::Sqrt(ptElectron*ptElectron + ms2);
 174     p[0] = ptElectron*TMath::Cos(phiElectron);
 175     p[1] = ptElectron*TMath::Sin(phiElectron);
 176     p[2] = mt*TMath::SinH(yElectron);
 177     id =  11;
 178     PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,1);
 179     //
 180     // Produce positron
 181     mt = TMath::Sqrt(ptPositron*ptPositron + ms2);
 182     p[0] = ptPositron*TMath::Cos(phiPositron);
 183     p[1] = ptPositron*TMath::Sin(phiPositron);
 184     p[2] = mt*TMath::SinH(yPositron);
 185     id = -11;
 186     PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,1);
 187     //
 188   }
 189   fEvent++;
 190   //
 191   fHeader.SetNProduced(2*npairs);
 192   fHeader.SetEventWeight(1);
 193   fHeader.SetInteractionTime(fTimeOrigin);
 194   AddHeader(&fHeader);
 195 }
 196
 197 //__________________________________________________________
 198 void AliGenQEDBg::SetLumiIntTime(double lumi, double intTime)
 199 {
 200   // assign luminosity and integration time
 201   if (lumi<=0)    AliFatal(Form("Luminosity must be positive, in cm^-2*s^-1, %e asked",lumi));
 202   if (intTime<=0) AliFatal(Form("Integration time must be positive, in seconnds, %e asked",intTime));
 203   fLumi = lumi;
 204   fIntTime = intTime;
 205   //
 206 }
 207
 208 //__________________________________________________________
 209 void AliGenQEDBg::SetMinMaxXSTest(double mn,double mx)
 210 {
 211   // set min,max number of generator calls for xsection estimates
 212   fMinXSTest = mn>100 ? mn : 100.;
 213   fMaxXSTest = mx>mx ? mx : mx+100.;
 214 }