Moving required CMake version from 2.8.4 to 2.8.8
[u/mrichter/AliRoot.git] / EVGEN / AliGenMC.cxx
index a6f8963..c7f1b95 100644 (file)
 #include <TMath.h>
 #include <TPDGCode.h>
 #include <TParticle.h>
+#include <TLorentzVector.h>
+#include <TVector3.h>
 
 #include "AliGenMC.h"
+#include "AliGenEventHeader.h"
+#include "AliRun.h"
 #include "AliGeometry.h"
+#include "AliDecayer.h"
 
 ClassImp(AliGenMC)
 
 AliGenMC::AliGenMC()
     :AliGenerator(),
-     fParticles(0),
+     fParticles("TParticle", 1000),
      fParentSelect(8),
      fChildSelect(8),
      fCutOnChild(0),
@@ -51,23 +56,24 @@ AliGenMC::AliGenMC()
      fXingAngleY(0.),
      fForceDecay(kAll),
      fMaxLifeTime(1.e-15),
-     fAProjectile(1),
-     fZProjectile(1),
-     fATarget(1),
-     fZTarget(1),
-     fProjectile("P"),
-     fTarget("P"),    
      fDyBoost(0.),
      fGeometryAcceptance(0),
      fPdgCodeParticleforAcceptanceCut(0),
-     fNumberOfAcceptedParticles(2)
+     fNumberOfAcceptedParticles(0),
+     fNprimaries(0)
 {
 // Default Constructor
+  fAProjectile = 1;
+  fZProjectile = 1;
+  fATarget = 1;
+  fZTarget = 1;
+  fProjectile = "P";
+  fTarget = "P";
 }
 
 AliGenMC::AliGenMC(Int_t npart)
     :AliGenerator(npart),
-     fParticles(0),
+     fParticles("TParticle", 1000),
      fParentSelect(8),
      fChildSelect(8),
      fCutOnChild(0),
@@ -85,20 +91,21 @@ AliGenMC::AliGenMC(Int_t npart)
      fXingAngleY(0.),
      fForceDecay(kAll),
      fMaxLifeTime(1.e-15),
-     fAProjectile(1),
-     fZProjectile(1),
-     fATarget(1),
-     fZTarget(1),
-     fProjectile("P"),
-     fTarget("P"),    
      fDyBoost(0.),
      fGeometryAcceptance(0),
      fPdgCodeParticleforAcceptanceCut(0),
-     fNumberOfAcceptedParticles(2)
+     fNumberOfAcceptedParticles(0),
+     fNprimaries(0)
 {
 //  Constructor
 // 
-    for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0;
+  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()
@@ -116,12 +123,16 @@ void AliGenMC::Init()
     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:
@@ -129,6 +140,7 @@ void AliGenMC::Init()
     case kWToMuon:
     case kWToCharmToMuon:
     case kZDiMuon:
+    case kZDiElectron:
        fChildSelect[0]=kMuonMinus;
        break;
     case kWToCharm:
@@ -142,14 +154,36 @@ void AliGenMC::Init()
        fChildSelect[0]=kKPlus;
        break;
     case kBJpsi:
-       fChildSelect[0]=443;
+    case kBJpsiUndecayed:
+       fChildSelect[0]= 443;
        break;
-    case kOmega:       
+   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;
     }
 
@@ -181,17 +215,21 @@ 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 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;
-    Float_t phi   = particle->Phi();
+    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);
+
     
     Double_t y, y0;
 
@@ -348,16 +386,16 @@ void AliGenMC::Boost()
 //
 
     Double_t beta  = TMath::TanH(fDyBoost);
-    Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
+    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", 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();
@@ -370,3 +408,54 @@ void AliGenMC::Boost()
        iparticle->SetMomentum(px, py, pzb, eb);
     }
 }
+
+void AliGenMC::BeamCrossAngle()
+{
+  // Applies a boost in the y-direction in order to take into account the 
+  // beam crossing angle
+
+  Double_t thetaPr0, phiPr0, pyPr2, pzPr2;
+  TVector3 beta;
+  
+  thetaPr0 = fXingAngleY / 2.;
+  phiPr0 = 0;
+
+  // 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);      
+    }
+}