CbmEventGenerator.cxx

Go to the documentation of this file.

/* Copyright (C) 2019-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Volker Friese [committer] */

 
#include "CbmEventGenerator.h"
 
#include "CbmBeam.h"
 
#include <FairGenericStack.h>
#include <FairMCEventHeader.h>
#include <Logger.h>
 
#include <TDatabasePDG.h>
#include <TMath.h>
#include <TRandom.h>
#include <TVector3.h>
 
#include <cassert>
 
using std::cerr;
using std::cout;
using std::endl;
 
 
// -----   Constructor   ---------------------------------------------------
CbmEventGenerator::CbmEventGenerator()
  : FairPrimaryGenerator()
  , fBeamProfile()
  , fTarget()
  , fForceVertexInTarget(kTRUE)
  , fForceVertexAtZ(kFALSE)
  , fVertexZ(0.)
{
  // Mother class members
  fName      = "EventGenerator";
  fBeamAngle = kTRUE;
}
// -------------------------------------------------------------------------
 
 
// -----   Destructor   ----------------------------------------------------
CbmEventGenerator::~CbmEventGenerator() {}
// -------------------------------------------------------------------------
 
 
// -----   Add a track to the stack   --------------------------------------
// Code copied from fairroot/base/sim/FairPrimary Generator.cxx,
// with the implementation of the beam angle rotation being modified.
void CbmEventGenerator::AddTrack(Int_t pdgid, Double_t px_raw, Double_t py_raw, Double_t pz_raw, Double_t vx,
                                 Double_t vy, Double_t vz, Int_t parent, Bool_t wanttracking, Double_t e, Double_t tof,
                                 Double_t weight, TMCProcess proc)
{
  // ---> Add event vertex to track vertex
  vx += fVertex.X();
  vy += fVertex.Y();
  vz += fVertex.Z();
 
  TVector3 mom(px_raw, py_raw, pz_raw);
 
  // ---> Rotate the track momentum in the transverse (x-y) plane by the event plane angle.
  if (fEventPlane) mom.RotateZ(fPhi);
 
  // ---> Rotate the track momentum to account for the beam direction, which may not be
  // in the z axis. This implementation deviates from that in fairroot, where TVector3::RotateUz() is used.
  // We found this to behave wrongly and use here TVector3::Rotate() instead.
  if (fBeamAngle) {
    TVector3 axis(-fBeamDirection.Unit().Y(), fBeamDirection.Unit().X(), 0);
    auto angle = TMath::ACos(fBeamDirection.Unit().Z());
    mom.Rotate(angle, axis);
  }
 
  // ---> Convert K0 and AntiK0 into K0s and K0L (random decision with equal probability)
  if (pdgid == 311 || pdgid == -311) {
    Double_t test = gRandom->Uniform(0., 1.);
    if (test >= 0.5) {
      pdgid = 310;
    }  // K0S
    else {
      pdgid = 130;
    }  // K0L
  }
 
  // ---> Check whether particle type is in PDG Database
  TDatabasePDG* pdgBase = TDatabasePDG::Instance();
  if (!pdgBase) {
    Fatal("CbmEventGenerator", "No TDatabasePDG instantiated");
  }
  else {
    TParticlePDG* pdgPart = pdgBase->GetParticle(pdgid);
    if (!pdgPart) {
      if (e < 0) {
        cerr << "\033[5m\033[31m -E CbmEventGenerator: PDG code " << pdgid << " not found in database.\033[0m " << endl;
        cerr << "\033[5m\033[31m -E CbmEventGenerator: Discarding particle "
                "\033[0m "
             << endl;
        cerr << "\033[5m\033[31m -E CbmEventGenerator: now MC Index is "
                "corrupted \033[0m "
             << endl;
        return;
      }
      else {
        cout << "\033[5m\033[31m -W CbmEventGenerator: PDG code " << pdgid
             << " not found in database. This warning can be savely "
                "ignored.\033[0m "
             << endl;
      }
    }
  }
  // ---> Get mass and calculate energy of particle
  if (e < 0) {
    Double_t mass = pdgBase->GetParticle(pdgid)->Mass();
    e             = TMath::Sqrt(mom.Mag2() + mass * mass);
  }  // else, use the value of e given to the function
 
  // ---> Set all other parameters required by PushTrack
  Int_t doTracking = 0;  // Go to tracking
  if (fdoTracking && wanttracking) {
    doTracking = 1;
  }
  Int_t dummyparent = -1;  // Primary particle (now the value is -1 by default)
  Double_t polx     = 0.;  // Polarisation
  Double_t poly     = 0.;
  Double_t polz     = 0.;
  Int_t ntr         = 0;  // Track number; to be filled by the stack
  Int_t status      = 0;  // Generation status
 
  if (parent != -1) {
    parent += fMCIndexOffset;
  }  // correct for tracks which are in list before generator is called
  // Add track to stack
  fStack->PushTrack(doTracking, dummyparent, pdgid, mom.X(), mom.Y(), mom.Z(), e, vx, vy, vz, tof, polx, poly, polz,
                    proc, ntr, weight, status, parent);
  fNTracks++;
}
// -------------------------------------------------------------------------
 
 
// -----   Force vertex at a given z position   ----------------------------
void CbmEventGenerator::ForceVertexAtZ(Double_t zVertex)
{
  fForceVertexAtZ      = kTRUE;
  fForceVertexInTarget = kFALSE;
  fVertexZ             = zVertex;
}
// -------------------------------------------------------------------------
 
 
// -----   Generate the input event   --------------------------------------
Bool_t CbmEventGenerator::GenerateEvent(FairGenericStack* stack)
{
 
  std::cout << std::endl;
  LOG(info) << GetName() << ": Generate event " << ++fEventNr;
 
  // An MCEventHeader must be present
  assert(fEvent);
 
  // Set the stack
  assert(stack);
  fStack = stack;
 
  // Reset MCEventHeader and track counter
  fEvent->Reset();
  fNTracks = 0;
 
  // Generate the event vertex and the beam angle
  MakeVertex();
 
  // Generate the event plane angle
  if (fEventPlane) {
    MakeEventPlane();
    LOG(info) << GetName() << ": Rotate event by " << fPhi << " rad";
  }
 
  // Call the registered generator classes
  fListIter->Reset();
  TObject* item = nullptr;
  while ((item = fListIter->Next())) {
    FairGenerator* generator = dynamic_cast<FairGenerator*>(item);
    assert(generator);
    fMCIndexOffset = fNTracks;  // Number of tracks before call of generator
    if (!generator->ReadEvent(this)) {
      LOG(info) << GetName() << ": ReadEvent failed for generator " << generator->GetName();
      return kFALSE;
    }
  }
 
  // Update event header
  if (!fEvent->IsSet()) fEvent->SetEventID(fEventNr);
  fEvent->SetNPrim(fNTracks);
  fEvent->SetVertex(fVertex);
  fEvent->SetRotX(fBeamAngleX);
  fEvent->SetRotY(fBeamAngleY);
  Double_t phiGen = fEvent->GetRotZ();  // from concrete generators
  fEvent->SetRotZ(phiGen + fPhi);       // total event plane angle
  // Event info to screen
  LOG(info) << GetName() << ": Event ID " << fEvent->GetEventID() << ", tracks " << fNTracks << ", vertex ("
            << fVertex.X() << ", " << fVertex.Y() << ", " << fVertex.Z() << ") cm";
  LOG(info) << GetName() << ": Beam angle (" << fBeamAngleX << ", " << fBeamAngleY << ") rad, event plane angle "
            << fPhi << " rad (total " << fEvent->GetRotZ() << " rad) change";
 
  // Update global track counter (who knows what it's good for ...)
  fTotPrim += fNTracks;
 
  return kTRUE;
}
// -------------------------------------------------------------------------
 
 
// -----   Generate event vertex   -----------------------------------------
void CbmEventGenerator::MakeVertex()
{
  if (fForceVertexAtZ) {
    MakeVertexAtZ();
  }
  else if (fForceVertexInTarget && fTarget) {
    MakeVertexInTarget();
  }
  else {
    MakeVertexInFocalPlane();
  }
}
// -------------------------------------------------------------------------
 
 
// -----   Generate event vertex in the beam focal plane   -----------------
void CbmEventGenerator::MakeVertexAtZ()
{
 
  std::unique_ptr<CbmBeam> beam;
  beam = fBeamProfile.GenerateBeam();
  TVector3 point(0., 0., fVertexZ);  
  TVector3 norm(0, 0., 1.);          
  fVertex        = beam->ExtrapolateToPlane(point, norm);
  fBeamAngleX    = beam->GetThetaX();
  fBeamAngleY    = beam->GetThetaY();
  fBeamDirection = beam->GetDirection();
}
// -------------------------------------------------------------------------
 
 
// -----   Generate event vertex in the beam focal plane   -----------------
void CbmEventGenerator::MakeVertexInFocalPlane()
{
 
  std::unique_ptr<CbmBeam> beam;
  beam           = fBeamProfile.GenerateBeam();
  fVertex        = beam->GetPosition();
  fBeamAngleX    = beam->GetThetaX();
  fBeamAngleY    = beam->GetThetaY();
  fBeamDirection = beam->GetDirection();
}
// -------------------------------------------------------------------------
 
 
// -----   Generate event vertex   -----------------------------------------
void CbmEventGenerator::MakeVertexInTarget()
{
 
  assert(fTarget);
 
  // Target surface centres and normal vector
  TVector3 surf1 = fTarget->GetSurfaceCentreUp();
  TVector3 surf2 = fTarget->GetSurfaceCentreDown();
  TVector3 norm  = fTarget->GetNormal();
 
  // Declare some vectors
  TVector3 point1;  // Intersection of beam with first target surface
  TVector3 point2;  // Intersection of beam with second target surface
 
  // It is possible that a generated beam trajectory misses the target.
  // In that case, the beam will be re-sampled until it hits the target.
  Bool_t isInTarget = kFALSE;
  Int_t nSamples    = 0;
  std::unique_ptr<CbmBeam> beam;
  while (!isInTarget) {
 
    // Abort if too many beam samples
    if (nSamples > 1000.)
      LOG(fatal) << GetName() << ": Aborting after " << nSamples << " beam samplings. Adjust beam and target.";
 
    // Sample a beam trajectory
    beam = fBeamProfile.GenerateBeam();
    nSamples++;
 
    // Intersections of beam trajectory with target surfaces
    point1 = beam->ExtrapolateToPlane(surf1, norm);
    point2 = beam->ExtrapolateToPlane(surf2, norm);
 
    // Check whether both intersections are inside the target
    // (i.e., the beam crosses both target surfaces)
    Bool_t check1 = ((point1 - surf1).Mag() < 0.5 * fTarget->GetDiameter());
    Bool_t check2 = ((point2 - surf2).Mag() < 0.5 * fTarget->GetDiameter());
    isInTarget    = check1 && check2;
 
  }  //? Beam not in target
 
  LOG(debug) << beam->ToString() << ", generated after " << nSamples << (nSamples > 1 ? " samplings: " : " sampling");
 
  // Sample a point between entry and exit of beam in target
  Double_t scale = 0.5;
  if (fSmearVertexZ) scale = gRandom->Uniform();
  fVertex = point1 + scale * (point2 - point1);
 
  // Set the beam angles (member variables of base class)
  fBeamAngleX    = beam->GetThetaX();
  fBeamAngleY    = beam->GetThetaY();
  fBeamDirection = beam->GetDirection();
}
// -------------------------------------------------------------------------
 
 
// -----   Info   ----------------------------------------------------------
void CbmEventGenerator::Print(Option_t*) const
{
  LOG(info) << fBeamProfile.ToString();
  if (fTarget) LOG(info) << fTarget->ToString();
  LOG(info) << "Vertex smearing along beam " << (fSmearVertexZ ? "ON" : "OFF");
  if (fEventPlane) {
    LOG(info) << "Random event plane angle between " << fPhiMin << " and " << fPhiMax << " rad";
  }
  else {
    LOG(info) << "Fixed event plane angle = 0";
  }
  LOG(info) << "Number of generators " << fGenList->GetEntries();
  for (Int_t iGen = 0; iGen < fGenList->GetEntries(); iGen++) {
    fGenList->At(iGen)->Print();
  }
}
// -------------------------------------------------------------------------
 
 
// -----   Set beam angle   ------------------------------------------------
void CbmEventGenerator::SetBeamAngle(Double_t meanThetaX, Double_t meanThetaY, Double_t sigmaThetaX,
                                     Double_t sigmaThetaY)
{
  fBeamProfile.SetAngle(meanThetaX, meanThetaY, sigmaThetaX, sigmaThetaY);
}
// -------------------------------------------------------------------------
 
 
// -----   Set beam position   ---------------------------------------------
void CbmEventGenerator::SetBeamPosition(Double_t meanX, Double_t meanY, Double_t sigmaX, Double_t sigmaY, Double_t zF)
{
  fBeamProfile.SetPosition(meanX, meanY, sigmaX, sigmaY, zF);
}
// -------------------------------------------------------------------------
 
 
ClassImp(CbmEventGenerator)