X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=EVGEN%2FAliGenMC.cxx;h=accdcb5d4a1456196b8fcdf84c061e52e2defda4;hb=1a468a3b24d18b41a25bc96975c834da2091b77b;hp=539d6c8cfc0c599b8249dc2e4ae6e87c9e325811;hpb=3b945a60dc1bdfa9caa06efbe8166a0d80fa79d3;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenMC.cxx b/EVGEN/AliGenMC.cxx index 539d6c8cfc0..accdcb5d4a1 100644 --- a/EVGEN/AliGenMC.cxx +++ b/EVGEN/AliGenMC.cxx @@ -13,55 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.14 2003/01/14 10:50:19 alibrary -Cleanup of STEER coding conventions - -Revision 1.13 2002/10/14 14:55:35 hristov -Merging the VirtualMC branch to the main development branch (HEAD) - -Revision 1.5.4.2 2002/07/24 08:56:28 alibrary -Updating EVGEN on TVirtulaMC - -Revision 1.12 2002/07/19 11:42:33 morsch -Use CalcMass() - -Revision 1.11 2002/06/06 15:26:24 morsch -Correct child-selection for kPhiKK - -Revision 1.10 2002/06/05 14:05:46 morsch -Decayer option kPhiKK for forced phi->K+K- decay added. - -Revision 1.9 2002/05/30 14:58:29 morsch -Add pointer to AliGeometry to handle geometrical acceptance. (G. MArtinez) - -Revision 1.8 2002/04/26 10:42:35 morsch -Case kNoDecayHeavy added. (N. Carrer) - -Revision 1.7 2002/04/17 10:32:32 morsch -Coding Rule violations corrected. - -Revision 1.6 2002/03/26 14:19:36 morsch -Saver calculation of rapdity. - -Revision 1.5 2002/03/12 17:02:20 morsch -Change in calculation of rapidity, include case in which numerically e == pz. - -Revision 1.4 2001/11/27 13:13:07 morsch -Maximum lifetime for long-lived particles to be put on the stack is parameter. -It can be set via SetMaximumLifetime(..). - -Revision 1.3 2001/10/16 08:48:56 morsch -Common vertex related code moved to base class AliGenerator. - -Revision 1.2 2001/10/15 08:15:51 morsch -Event vertex and vertex truncation setting moved into AliMC. - -Revision 1.1 2001/07/13 10:56:00 morsch -AliGenMC base class for AliGenParam and AliGenPythia commonalities. - -*/ +/* $Id$ */ // Base class for generators using external MC generators. // For example AliGenPythia using Pythia. @@ -69,57 +21,91 @@ AliGenMC base class for AliGenParam and AliGenPythia commonalities. // decay products and particle selection. // andreas.morsch@cern.ch +#include #include #include #include +#include +#include #include "AliGenMC.h" +#include "AliGenEventHeader.h" +#include "AliRun.h" +#include "AliGeometry.h" +#include "AliDecayer.h" - ClassImp(AliGenMC) +ClassImp(AliGenMC) AliGenMC::AliGenMC() - :AliGenerator() + :AliGenerator(), + fParticles("TParticle", 1000), + fParentSelect(8), + fChildSelect(8), + fCutOnChild(0), + fChildPtMin(0.), + fChildPtMax(1.e10), + fChildPMin(0.), + fChildPMax(1.e10), + fChildPhiMin(0.), + fChildPhiMax(2. * TMath::Pi()), + fChildThetaMin(0.), + fChildThetaMax(TMath::Pi()), + fChildYMin(-12.), + fChildYMax(12.), + fXingAngleX(0.), + fXingAngleY(0.), + fForceDecay(kAll), + fMaxLifeTime(1.e-15), + fDyBoost(0.), + fGeometryAcceptance(0), + fPdgCodeParticleforAcceptanceCut(0), + fNumberOfAcceptedParticles(0), + fNprimaries(0) { // Default Constructor - SetCutOnChild(); - SetChildMomentumRange(); - SetChildPtRange(); - SetChildPhiRange(); - SetChildThetaRange(); - SetChildYRange(); - SetMaximumLifetime(); - SetGeometryAcceptance(); - SetPdgCodeParticleforAcceptanceCut(); - SetNumberOfAcceptedParticles(); - SetTarget(); - SetProjectile(); + fAProjectile = 1; + fZProjectile = 1; + fATarget = 1; + fZTarget = 1; + fProjectile = "P"; + fTarget = "P"; } AliGenMC::AliGenMC(Int_t npart) - :AliGenerator(npart) + :AliGenerator(npart), + fParticles("TParticle", 1000), + fParentSelect(8), + fChildSelect(8), + fCutOnChild(0), + fChildPtMin(0.), + fChildPtMax(1.e10), + fChildPMin(0.), + fChildPMax(1.e10), + fChildPhiMin(0.), + fChildPhiMax(2. * TMath::Pi()), + fChildThetaMin(0.), + fChildThetaMax(TMath::Pi()), + fChildYMin(-12.), + fChildYMax(12.), + fXingAngleX(0.), + fXingAngleY(0.), + fForceDecay(kAll), + fMaxLifeTime(1.e-15), + fDyBoost(0.), + fGeometryAcceptance(0), + fPdgCodeParticleforAcceptanceCut(0), + fNumberOfAcceptedParticles(0), + fNprimaries(0) { // Constructor - SetCutOnChild(); - SetChildMomentumRange(); - SetChildPtRange(); - SetChildPhiRange(); - SetChildThetaRange(); - SetChildYRange(); // - fParentSelect.Set(8); - fChildSelect.Set(8); - for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0; - SetMaximumLifetime(); - SetGeometryAcceptance(); - SetPdgCodeParticleforAcceptanceCut(); - SetNumberOfAcceptedParticles(); - SetTarget(); - SetProjectile(); -} - -AliGenMC::AliGenMC(const AliGenMC & mc) -{ -// copy constructor + fAProjectile = 1; + fZProjectile = 1; + fATarget = 1; + fZTarget = 1; + fProjectile = "P"; + fTarget = "P"; + for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0; } AliGenMC::~AliGenMC() @@ -132,32 +118,82 @@ void AliGenMC::Init() // // Initialization switch (fForceDecay) { + case kBSemiElectronic: case kSemiElectronic: case kDiElectron: case kBJpsiDiElectron: case kBPsiPrimeDiElectron: + case kElectronEM: + case kDiElectronEM: fChildSelect[0] = kElectron; break; + case kHardMuons: + case kBSemiMuonic: + case kDSemiMuonic: case kSemiMuonic: case kDiMuon: + case kJpsiDiMuon: case kBJpsiDiMuon: case kBPsiPrimeDiMuon: case kPiToMu: case kKaToMu: + case kWToMuon: + case kWToCharmToMuon: + case kZDiMuon: + case kZDiElectron: fChildSelect[0]=kMuonMinus; break; + case kWToCharm: + break; case kHadronicD: + case kHadronicDWithout4Bodies: + case kHadronicDWithV0: + case kHadronicDWithout4BodiesWithV0: fChildSelect[0]=kPiPlus; fChildSelect[1]=kKPlus; break; case kPhiKK: fChildSelect[0]=kKPlus; - case kOmega: + break; + case kBJpsi: + case kBJpsiUndecayed: + fChildSelect[0]= 443; + break; + case kChiToJpsiGammaToMuonMuon: + fChildSelect[0]= 22; + fChildSelect[1]= 13; + break; + case kChiToJpsiGammaToElectronElectron: + fChildSelect[0]= 22; + fChildSelect[1]= 11; + break; + case kLambda: + fChildSelect[0]= kProton; + fChildSelect[1]= 211; + break; + case kPsiPrimeJpsiDiElectron: + fChildSelect[0]= 211; + fChildSelect[1]= 11; + break; + case kGammaEM: + fChildSelect[0] = kGamma; + break; + case kOmega: case kAll: + case kAllMuonic: case kNoDecay: case kNoDecayHeavy: + case kNoDecayBeauty: + case kNeutralPion: + case kBeautyUpgrade: break; } + + if (fZTarget != 0 && fAProjectile != 0) + { + fDyBoost = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / + (Double_t(fZTarget) * Double_t(fAProjectile))); + } } @@ -181,20 +217,22 @@ Bool_t AliGenMC::ChildSelected(Int_t ip) const return kFALSE; } -Bool_t AliGenMC::KinematicSelection(TParticle *particle, Int_t flag) const +Bool_t AliGenMC::KinematicSelection(const TParticle *particle, Int_t flag) const { // Perform kinematic selection - Float_t px = particle->Px(); - Float_t py = particle->Py(); - Float_t pz = particle->Pz(); - Float_t e = particle->Energy(); - Float_t pt = particle->Pt(); - Float_t p = particle->P(); - Float_t theta = particle->Theta(); - Float_t mass = particle->GetCalcMass(); - Float_t mt2 = pt * pt + mass * mass; + Double_t pz = particle->Pz(); + Double_t pt = particle->Pt(); + Double_t p = particle->P(); + Double_t theta = particle->Theta(); + Double_t mass = particle->GetCalcMass(); + Double_t mt2 = pt * pt + mass * mass; + Double_t phi = particle->Phi(); + Double_t e = particle->Energy(); + + if (e == 0.) + e = TMath::Sqrt(p * p + mass * mass); + - Float_t phi = Float_t(TMath::ATan2(Double_t(py),Double_t(px))); Double_t y, y0; if (TMath::Abs(pz) < e) { @@ -288,24 +326,26 @@ Bool_t AliGenMC::KinematicSelection(TParticle *particle, Int_t flag) const Bool_t AliGenMC::CheckAcceptanceGeometry(Int_t np, TClonesArray* particles) { - Bool_t Check ; // All fPdgCodeParticleforAcceptanceCut particles are in in the fGeometryAcceptance acceptance - Int_t NumberOfPdgCodeParticleforAcceptanceCut=0; - Int_t NumberOfAcceptedPdgCodeParticleforAcceptanceCut=0; +// Check the geometrical acceptance for particle. + + Bool_t check ; + Int_t numberOfPdgCodeParticleforAcceptanceCut = 0; + Int_t numberOfAcceptedPdgCodeParticleforAcceptanceCut = 0; TParticle * particle; Int_t i; - for (i=0; iAt(i); if( TMath::Abs( particle->GetPdgCode() ) == TMath::Abs( fPdgCodeParticleforAcceptanceCut ) ) { - NumberOfPdgCodeParticleforAcceptanceCut++; - if (fGeometryAcceptance->Impact(particle)) NumberOfAcceptedPdgCodeParticleforAcceptanceCut++; + numberOfPdgCodeParticleforAcceptanceCut++; + if (fGeometryAcceptance->Impact(particle)) numberOfAcceptedPdgCodeParticleforAcceptanceCut++; } } - if ( NumberOfAcceptedPdgCodeParticleforAcceptanceCut > (fNumberOfAcceptedParticles-1) ) - Check = kTRUE; + if ( numberOfAcceptedPdgCodeParticleforAcceptanceCut > (fNumberOfAcceptedParticles-1) ) + check = kTRUE; else - Check = kFALSE; + check = kFALSE; - return Check; + return check; } Int_t AliGenMC::CheckPDGCode(Int_t pdgcode) const @@ -341,23 +381,23 @@ Int_t AliGenMC::CheckPDGCode(Int_t pdgcode) const return pdgcode; } -void AliGenMC::Boost(Float_t dy) +void AliGenMC::Boost() { // // Boost cms into LHC lab frame // - Double_t beta = TMath::TanH(dy); - Double_t gamma = 1./TMath::Sqrt(1.-beta*beta); + Double_t beta = TMath::TanH(fDyBoost); + Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta)); Double_t gb = gamma * beta; - printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma); + // printf("\n Boosting particles to lab frame %f %f %f", fDyBoost, beta, gamma); Int_t i; - Int_t np = fParticles->GetEntriesFast(); + Int_t np = fParticles.GetEntriesFast(); for (i = 0; i < np; i++) { - TParticle* iparticle = (TParticle*) fParticles->At(i); + TParticle* iparticle = (TParticle*) fParticles.At(i); Double_t e = iparticle->Energy(); Double_t px = iparticle->Px(); @@ -371,11 +411,52 @@ void AliGenMC::Boost(Float_t dy) } } - - -AliGenMC& AliGenMC::operator=(const AliGenMC& rhs) +void AliGenMC::BeamCrossAngle() { -// Assignment operator - return *this; + // Applies a boost in the y-direction in order to take into account the + // beam crossing angle + + Double_t thetaPr0, pyPr2, pzPr2; + TVector3 beta; + + thetaPr0 = fXingAngleY / 2.; + + // Momentum of the CMS system + pyPr2 = TMath::Sqrt(fEnergyCMS * fEnergyCMS/ 4 - 0.938 * 0.938) * TMath::Sin(thetaPr0); + pzPr2 = TMath::Sqrt(fEnergyCMS * fEnergyCMS/ 4 - 0.938 * 0.938) * TMath::Cos(thetaPr0); + + TLorentzVector proj1Vect, proj2Vect, projVect; + proj1Vect.SetPxPyPzE(0., pyPr2, pzPr2, fEnergyCMS/2); + proj2Vect.SetPxPyPzE(0., pyPr2,-pzPr2, fEnergyCMS/2); + projVect = proj1Vect + proj2Vect; + beta=(1. / projVect.E()) * (projVect.Vect()); + + Int_t i; + Int_t np = fParticles.GetEntriesFast(); + for (i = 0; i < np; i++) + { + TParticle* iparticle = (TParticle*) fParticles.At(i); + + Double_t e = iparticle->Energy(); + Double_t px = iparticle->Px(); + Double_t py = iparticle->Py(); + Double_t pz = iparticle->Pz(); + + TLorentzVector partIn; + partIn.SetPxPyPzE(px,py,pz,e); + partIn.Boost(beta); + iparticle->SetMomentum(partIn.Px(),partIn.Py(),partIn.Pz(),partIn.E()); + } } + +void AliGenMC::AddHeader(AliGenEventHeader* header) +{ + // Passes header either to the container or to gAlice + if (fContainer) { + header->SetName(fName); + fContainer->AddHeader(header); + } else { + gAlice->SetGenEventHeader(header); + } +}