Merge pull request #1 from shahor02/fix_qed

Move Xsection (barn) calculation to TGenEpEmv1

TEPEMGEN/TGenEpEmv1.cxx

@@ -89,7 +89,11 @@ TGenEpEmv1::TGenEpEmv1():
   fYMin(-100.), fYMax(100.),
   fPhiMin(0.), fPhiMax(360.),
   fPtMin(0.), fPtMax(1.e10),
-  fTimeOrigin(0.)
+  fTimeOrigin(0.),
+  fXSection(-1.),
+  fXSectionEps(1e-2),
+  fMinXSTest(1000),
+  fMaxXSTest(10000000)
 {
   // Default constructor
   for (Int_t i = 0; i < 3; ++i) {
@@ -115,8 +119,46 @@ void TGenEpEmv1::Init()
   if (fPtMin == 0) fPtMin = 1.E-04; // avoid zero pT
   Initialize(fYMin, fYMax, fPtMin, fPtMax);
   fEvent = 0;
+  //
+  // calculate XSection
+  double err = 0;
+  fXSection = CalcXSection(fXSectionEps,fMinXSTest,fMaxXSTest,err);
+  if (fXSection<=0 || err/fXSection>fXSectionEps) {
+    abort();
+  }
+  fXSectionEps = err/fXSection;
 }
 
+//____________________________________________________________
+double TGenEpEmv1::CalcXSection(double eps, int triMin, int triMax, double& err)
+{
+  if (eps<1e-4) {
+    eps = 1e-4;
+  }
+  int ngen = 0;
+  printf("Estimating x-section with min.relative precision of %f and min/max test: %d/%d\n",
+	 eps,triMin,triMax);
+  double yElectron,yPositron,xElectron,xPositron,phi12,weight;
+  double xSect = -1;
+  err = -1;
+  //
+  do {
+    TEpEmGen::GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
+    if (++ngen>triMin) { // ensure min number of tests
+      xSect = TEpEmGen::GetXsection()*1000;
+      err = TEpEmGen::GetDsection()*1000;
+    }
+  } while(!((xSect>0 && err/xSect<eps) || ngen>triMax));
+  //
+  if (xSect<=0) {
+    printf("Failed to estimate X-section after %d trials\n",ngen);
+    abort();
+  }
+  printf("X-section = %e with %e error after %d trials",xSect,err,ngen);
+  return xSect;
+}
+
+
 //____________________________________________________________
 void TGenEpEmv1::GenerateEvent()
 {

TEPEMGEN/TGenEpEmv1.h

@@ -32,14 +32,19 @@ class TGenEpEmv1 : public TEpEmGen {
   void SetOrigin(Float_t ox, Float_t oy, Float_t oz) {fOrigin[0]=ox; fOrigin[1]=oy; fOrigin[2]=oz;};
   void SetSigma(Float_t sx, Float_t sy, Float_t sz) {fOsigma[0]=sx; fOsigma[1]=sy; fOsigma[2]=sz;};
   void SetTimeOrigin(Float_t timeorig) {fTimeOrigin = timeorig;};
+  void SetXSectionEps(double eps=1e-2)  {fXSectionEps = eps>0 ? eps:1e-2;}
+  void SetMinMaxXSTest(int mn,int mx);
+  double CalcXSection(double eps, int triMin, int triMax, double& err);
+  double GetXSection()            const {return fXSection;}
+  double GetXSectionEps()         const {return fXSectionEps;}
 
  protected:
   TGenEpEmv1(const TGenEpEmv1 & gen);
   TGenEpEmv1 & operator=(const TGenEpEmv1 & gen);
   void GeneratePair(Double_t vx, Double_t vy, Double_t vz, Double_t vt);
   void Rndm(Float_t *array, Int_t n) {gRandom->RndmArray(n, array);};
   void Rndm(Double_t *array, Int_t n) {gRandom->RndmArray(n, array);};
-
+  
   Float_t    fMass;    // electron mass
   Int_t      fDebug;   // debug level
   Int_t      fEvent;   // internal event number
@@ -49,6 +54,12 @@ class TGenEpEmv1 : public TEpEmGen {
   Double_t fPtMin, fPtMax;
   Double_t fTimeOrigin, fOrigin[3], fOsigma[3];
 
+  Double_t   fXSection;     // estimated cross section in barns
+  Double_t   fXSectionEps;  // error with wich Xsection is calculated
+  int        fMinXSTest;    // min number of generator calls for Xsection estimate
+  int        fMaxXSTest;    // max number of generator calls for Xsection estimate
+
+
   ClassDef(TGenEpEmv1,1) // Generator of single e+e- pair production in PbPb ultra-peripheral collisions
 };
 #endif

TEPEMGEN/TGenQEDBg.cxx

@@ -89,12 +89,8 @@ ClassImp(TGenQEDBg);
 
 TGenQEDBg::TGenQEDBg()
 :  fLumi(0)
-  ,fXSection(0)
-  ,fXSectionEps(1e-2)
   ,fIntTime(0)
   ,fPairsInt(-1)
-  ,fMinXSTest(1e3)
-  ,fMaxXSTest(1e7)
 {
 }
 
@@ -118,27 +114,9 @@ void TGenQEDBg::Init()
   TGenEpEmv1::Init();
   //
   fPairsInt = 0;
-  int ngen = 0;
-  printf("Estimating x-section with min.relative precision of %f and min/max test: %d/%d",
-	 fXSectionEps,int(fMinXSTest),int(fMaxXSTest));
-  //
-  double yElectron,yPositron,xElectron,xPositron,phi12,weight,err=0;
-  fXSection = -1;
-  do {
-    TEpEmGen::GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
-    if (++ngen>fMinXSTest) { // ensure min number of tests
-      fXSection = GetXsection();
-      err = GetDsection();
-    }
-  } while(!((fXSection>0 && err/fXSection<fXSectionEps) || ngen>fMaxXSTest));
-  //
-  if (fXSection<=0) {
-    printf("X-section = %e after %d trials, cannot generate",fXSection,ngen);
-    abort();
-  }
-  fPairsInt = fXSection*1e-21*fLumi*fIntTime; // xsestion is in kbarn!
-  printf("Estimated x-secion: %e+-%ekb in %d tests, <Npairs>=%e per %e time interval",
-	 fXSection,err,ngen,fPairsInt,fIntTime);
+  fPairsInt = fXSection*1e-24*fLumi*fIntTime; // xsestion is in barn!
+  printf("Estimated x-secion: %e+-%eb, <Npairs>=%e per %e time interval",
+	 fXSection,fXSectionEps*fXSection,fPairsInt,fIntTime);
   //
 }
 
@@ -192,11 +170,3 @@ void TGenQEDBg::SetLumiIntTime(double lumi, double intTime)
   fIntTime = intTime;
   //
 }
-
-//__________________________________________________________
-void TGenQEDBg::SetMinMaxXSTest(double mn,double mx)
-{
-  // set min,max number of generator calls for xsection estimates
-  fMinXSTest = mn>100 ? mn : 100.;
-  fMaxXSTest = mx>mx ? mx : mx+100.;
-}

TEPEMGEN/TGenQEDBg.h

@@ -27,25 +27,17 @@ class TGenQEDBg : public TGenEpEmv1
   //
   Double_t  GetLuminosity()       const {return fLumi;}
   Double_t  GetIntegrationTime()  const {return fIntTime;}
-  Double_t  GetXSection()         const {return fXSection;}
-  Double_t  GetXSectionEps()      const {return fXSectionEps;}
   Double_t  GetMeanNPairs()       const {return fPairsInt;}
   //
   void      SetLumiIntTime(double lumi, double intTime);
-  void      SetXSectionEps(double eps=1e-2)  {fXSectionEps = eps>0 ? eps:1e-2;}
-  void      SetMinMaxXSTest(double mn,double mx);
   //
  protected:
   TGenQEDBg(const TGenQEDBg & gen);
   TGenQEDBg & operator=(const TGenQEDBg & gen);
   //
   Double_t   fLumi;         // beam luminsity
-  Double_t   fXSection;     // estimated cross section in k-barns
-  Double_t   fXSectionEps;  // error with wich Xsection is calculated
   Double_t   fIntTime;      // integration time in seconds
   Double_t   fPairsInt;     // estimated average number of pairs in IntTime
-  Double_t   fMinXSTest;    // min number of generator calls for Xsection estimate
-  Double_t   fMaxXSTest;    // max number of generator calls for Xsection estimate
 
   //
   ClassDef(TGenQEDBg,1) // Generator e+e- pair background from PbPb QED interactions