[u/mrichter/AliRoot.git] / TEPEMGEN / AliGenQEDBg.cxx

/**************************************************************************
 * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 *                                                                        *
 *                                                                        *
 * Copyright(c) 1997, 1998, 2002, Adrian Alscher and Kai Hencken          *
 * See $ALICE_ROOT/EpEmGen/diffcross.f for full Copyright notice          *
 *                                                                        *
 *                                                                        *
 * Copyright(c) 2002 Kai Hencken, Yuri Kharlov, Serguei Sadovsky          *
 * See $ALICE_ROOT/EpEmGen/epemgen.f for full Copyright notice            *
 *                                                                        *
 **************************************************************************/

/* $Id$ */

// Event generator for background from e+e- pair production in ultraperipheral PbPb collisions
// at 5.5 TeV/nucleon, integrated over specific readout cycle
// Derived from AliGenEpEmv1
// Author: ruben.shahoyan@cern.ch
//%
// References:
// [1] "Multiple electromagnetic electron positron pair production in
//      relativistic heavy ion collisions".
//      Adrian Alscher, Kai Hencken, Dirk Trautmann, and Gerhard Baur,
//      Phys. Rev. A55 (1997) 396.
// [2] K.Hencken, Yu.Kharlov, S.Sadovsky, Internal ALICE Note 2002-27.
//%
// Usage:
// Initialization:
//    AliGenQEDBg *gener = new AliGenQEDBg();
//    gener->SetXXXRange(); // Set kinematics range
//    gener->SetLumiIntTime(double lumi, double sec); // luminosity and intergration time in seconds
//    gener->Init();
// Event generation:
//    gener->Generate(); // Produce poissonian number of e+e- pair with average number
//    corresponding to requested integration time with given beam luminosity
//    Each pair has its own vertex, and time randomly distributed at originT 
//    and originT+ integration time                    
//
//    For details of pair generation see AliGenEpEmv1.cxx by Yuri.Kharlov@cern.ch

/*
  The most useful way of using it is in the overlay with other generators (using the cocktail):
  In the Config.C 
  // load the library:
  gSystem->Load("libTEPEMGEN");
  //
  // add to AliGenCocktail as:
  AliGenCocktail *cocktail = new AliGenCocktail();
  // ... setup other stuff
  // 
  // QED background
  AliGenQEDBg*  genBg = new AliGenQEDBg();
  genBg->SetEnergyCMS(5500);
  genBg->SetProjectile("A", 208, 82);
  genBg->SetTarget    ("A", 208, 82);
  genBg->SetYRange(-6.,3);
  genBg->SetPtRange(1.e-3,1.0);      // Set pt limits (GeV) for e+-: 1MeV corresponds to max R=13.3mm at 5kGaus
  genBg->SetLumiIntTime(6.e27,20e-6); // luminosity and integration time
  //
  cocktail->AddGenerator(genBg,"QEDep",1);
  // We need to generate independent vertex for every event, this must be set after adding to cocktail
  // Note that the IO origin and sigma's is transferred from AliGenCocktail to daughter generators
  genBg->SetVertexSource(kInternal);
*/

#include "AliLog.h"
#include "AliGenQEDBg.h"
#include <TParticle.h>
#include <TParticlePDG.h>
#include <TDatabasePDG.h>
#include <TEpEmGen.h>

ClassImp(AliGenQEDBg)

//------------------------------------------------------------

AliGenQEDBg::AliGenQEDBg()
:  fLumi(0)
  ,fXSection(0)
  ,fXSectionEps(1e-2)
  ,fIntTime(0)
  ,fPairsInt(-1)
  ,fMinXSTest(1e3)
  ,fMaxXSTest(1e7)
{
}

//____________________________________________________________
AliGenQEDBg::~AliGenQEDBg()
{
  // Destructor
}

//____________________________________________________________
void AliGenQEDBg::Init()
{
  // Initialisation:
  AliInfo(Form("Will estimate QED bg. for L=%e cm^-2*s^-1 and Integration Time of %e s.",fLumi,fIntTime));
  if (fLumi<=0 || fIntTime<=0) {
    AliWarning("One of parameters is not set properly, no pairs will be generated");
    return;
  }
  //
  // initialize the generator of e+e- pair production
  AliGenEpEmv1::Init();
  //
  fPairsInt = 0;
  int ngen = 0;
  AliInfo(Form("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 {
    fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
    if (++ngen>fMinXSTest) { // ensure min number of tests
      fXSection = fEpEmGen->GetXsection();
      err = fEpEmGen->GetDsection();
    }
  } while(!((fXSection>0 && err/fXSection<fXSectionEps) || ngen>fMaxXSTest));
  //
  if (fXSection<=0) {
    AliError(Form("X-section = %e after %d trials, cannot generate",fXSection,ngen));
    return;
  }
  fPairsInt = fXSection*1e-21*fLumi*fIntTime; // xsestion is in kbarn!
  AliInfo(Form("Estimated x-secion: %e+-%ekb in %d tests, <Npairs>=%e per %e time interval",
	       fXSection,err,ngen,fPairsInt,fIntTime));
  //
}

//____________________________________________________________
void AliGenQEDBg::Generate()
{
  //
  // Generate poissian <fPairsInt> e+e- pairs, each one with its vertex
  //
  Float_t polar[3]= {0,0,0};
  Float_t origin[3];
  Float_t time = 0.;
  Float_t p[3];
  //
  int npairs,nt,id;;
  if (fPairsInt<=0 || (npairs=gRandom->Poisson(fPairsInt))<1) return;
  //
  Double_t ptElectron,ptPositron, phiElectron,phiPositron, mt, ms2 = fMass*fMass;
  Double_t phi12=0,xElectron=0,xPositron=0,yElectron=0,yPositron=0,weight=0;
  //
  for (int i=0;i<npairs;i++) {
    // each pair has its own vertex and time
    Vertex();
    for (int j=0;j<3;j++) origin[j] = fVertex[j];
    time = fTimeOrigin+gRandom->Rndm()*fIntTime;
    //
    fEpEmGen->GenerateEvent(fYMin,fYMax,fPtMin,fPtMax,yElectron,yPositron,xElectron,xPositron,phi12,weight);
    ptElectron  = TMath::Power(10,xElectron) * 1.e-03;;
    ptPositron  = TMath::Power(10,xPositron) * 1.e-03;;
    phiElectron = fPhiMin + gRandom->Rndm() * (fPhiMax-fPhiMin);
    phiPositron = phiElectron + phi12;
    // Produce electron
    mt = TMath::Sqrt(ptElectron*ptElectron + ms2);
    p[0] = ptElectron*TMath::Cos(phiElectron);
    p[1] = ptElectron*TMath::Sin(phiElectron);
    p[2] = mt*TMath::SinH(yElectron);
    id =  11;
    PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,1);    
    //
    // Produce positron
    mt = TMath::Sqrt(ptPositron*ptPositron + ms2);
    p[0] = ptPositron*TMath::Cos(phiPositron);
    p[1] = ptPositron*TMath::Sin(phiPositron);
    p[2] = mt*TMath::SinH(yPositron);
    id = -11;
    PushTrack(fTrackIt,-1, id,p,origin,polar,time,kPPrimary,nt,1);
    //
  }
  fEvent++;
  //
  fHeader.SetNProduced(2*npairs);
  fHeader.SetEventWeight(1);
  fHeader.SetInteractionTime(fTimeOrigin);
  AddHeader(&fHeader);
}

//__________________________________________________________
void AliGenQEDBg::SetLumiIntTime(double lumi, double intTime)
{
  // assign luminosity and integration time
  if (lumi<=0)    AliFatal(Form("Luminosity must be positive, in cm^-2*s^-1, %e asked",lumi));
  if (intTime<=0) AliFatal(Form("Integration time must be positive, in seconnds, %e asked",intTime));
  fLumi = lumi;
  fIntTime = intTime;
  //
}

//__________________________________________________________
void AliGenQEDBg::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.;
}
Commit	Line	Data
0ebd3869	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	//
15a3d47f	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:
4070f709	57	gSystem->Load("libTEPEMGEN");
15a3d47f	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	*/
0ebd3869	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	//
bdc5420c	124	double yElectron,yPositron,xElectron,xPositron,phi12,weight,err=0;
0ebd3869	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 = fXSection1e-21fLumi*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	}