CbmPlutoGenerator.cxx

Go to the documentation of this file.

// -------------------------------------------------------------------------
// -----                CbmPlutoGenerator source file                  -----
// -----          Created 13/07/04  by V. Friese / D.Bertini           -----
// -------------------------------------------------------------------------
#include "CbmPlutoGenerator.h"
 
#include "FairLogger.h"
#include "FairPrimaryGenerator.h"  // for FairPrimaryGenerator
 
#include "PDataBase.h"    // for PDataBase
#include "PParticle.h"    // for PParticle
#include "PStaticData.h"  // for PStaticData
 
#include "TChain.h"          // for TChain
#include "TClonesArray.h"    // for TClonesArray
#include "TDatabasePDG.h"    // for TDatabasePDG
#include "TFile.h"           // for TFile
#include "TLorentzVector.h"  // for TLorentzVector
#include "TTree.h"           // for TTree
#include "TVector3.h"        // for TVector3
#include <iosfwd>            // for ostream
 
//#include <stddef.h>                     // for NULL
#include <iostream>  // for operator<<, basic_ostream, etc
#include <sys/stat.h>
 
// -----   Default constructor   ------------------------------------------
CbmPlutoGenerator::CbmPlutoGenerator()
  : FairGenerator()
  , fdata(makeStaticData())
  , fbase(makeDataBase())
  , iEvent(0)
  , fFileName(nullptr)
  , fInputChain(nullptr)
  , fParticles(nullptr)
  , fPDGmanual(0) {}
// ------------------------------------------------------------------------
 
// -----   Standard constructor   -----------------------------------------
CbmPlutoGenerator::CbmPlutoGenerator(const Char_t* fileName)
  : FairGenerator()
  , fdata(makeStaticData())
  , fbase(makeDataBase())
  , iEvent(0)
  , fFileName(fileName)
  , fInputChain(nullptr)
  , fParticles(new TClonesArray("PParticle", 100))
  , fPDGmanual(0) {
  fInputChain = new TChain("data");
  CheckFileExist(fileName);
 
  fInputChain->Add(fileName);
  fInputChain->SetBranchAddress("Particles", &fParticles);
}
// ------------------------------------------------------------------------
 
// -----  Constructor with file list   -----------------------------------------
CbmPlutoGenerator::CbmPlutoGenerator(std::vector<std::string> fileNames)
  : FairGenerator()
  , fdata(makeStaticData())
  , fbase(makeDataBase())
  , iEvent(0)
  , fFileName()
  , fInputChain(nullptr)
  , fParticles(new TClonesArray("PParticle", 100))
  , fPDGmanual(0) {
  fInputChain = new TChain("data");
  for (auto& name : fileNames) {
    CheckFileExist(name);
    fInputChain->Add(name.c_str());
  }
  fInputChain->SetBranchAddress("Particles", &fParticles);
}
 
 
// -----   Destructor   ---------------------------------------------------
CbmPlutoGenerator::~CbmPlutoGenerator() {
  // remove Pluto database
  delete fdata;
  delete fbase;
  CloseInput();
}
// ------------------------------------------------------------------------
 
 
// -----   Public method ReadEvent   --------------------------------------
Bool_t CbmPlutoGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
 
  // Check for input file
  if (!fInputChain) {
    LOG(error) << "CbmPlutoGenerator: Input file not open!";
    return kFALSE;
  }
 
  // Check for number of events in input file
  if (iEvent > fInputChain->GetEntries()) {
    LOG(error) << "CbmPlutoGenerator: No more events in input file!";
    CloseInput();
    return kFALSE;
  }
  fInputChain->GetEntry(iEvent++);
 
  // Get PDG database
  TDatabasePDG* dataBase = TDatabasePDG::Instance();
 
  // Get number of particle in TClonesrray
  Int_t nParts = fParticles->GetEntriesFast();
 
  // Loop over particles in TClonesArray
  for (Int_t iPart = 0; iPart < nParts; iPart++) {
    PParticle* part = (PParticle*) fParticles->At(iPart);
 
    Int_t* pdgType = 0x0;
    Bool_t found = fbase->GetParamInt("pid", part->ID(), "pythiakf", &pdgType);
    // TODO: replace by fdata->GetParticleKF(part->ID()); as soon as FairSoft uses pluto version 5.43 or higher and remove fbase
 
    // Check if particle type is known to database
    if (!found) {
      LOG(warn) << "CbmPlutoGenerator: Unknown type " << part->ID()
                << ", skipping particle.";
      continue;
    }
    LOG(info) << iPart << " Particle (geant " << part->ID() << " PDG "
              << *pdgType << " -> "
              << dataBase->GetParticle(*pdgType)->GetName();
 
    // set PDG by hand for pluto dilepton pairs and other not defined codes in pluto
    Int_t dielectron = 99009911;
    Int_t dimuon     = 99009913;
    if (fPDGmanual && *pdgType == 0) {
      pdgType = &fPDGmanual;
      LOG(warn) << "PDG code changed by user defintion to " << *pdgType;
    } else if (part->ID() == 51)
      pdgType = &dielectron;
    else if (part->ID() == 50)
      pdgType = &dimuon;
 
    // get the mother
    Int_t parIdx = part->GetParentIndex();
    // get daughter
    Int_t idx = part->GetDaughterIndex();
 
    TLorentzVector mom = part->Vect4();
    Double_t px        = mom.Px();
    Double_t py        = mom.Py();
    Double_t pz        = mom.Pz();
    Double_t ee        = mom.E();
 
    TVector3 vertex = part->getVertex();
    Double_t vx     = vertex.x();
    Double_t vy     = vertex.y();
    Double_t vz     = vertex.z();
 
    Bool_t wanttracking = kTRUE;
    if (idx > -1) wanttracking = kFALSE;  // only tracking for decay products
    Int_t parent = parIdx;
    LOG(info) << "Add particle with parent at index " << parIdx
              << " and do tracking " << wanttracking;
 
    // Give track to PrimaryGenerator
    primGen->AddTrack(
      *pdgType, px, py, pz, vx, vy, vz, parent, wanttracking, ee);
 
  }  //  Loop over particle in event
 
  return kTRUE;
}
// ------------------------------------------------------------------------
 
 
// -----   Private method CloseInput   ------------------------------------
void CbmPlutoGenerator::CloseInput() {
  if (fInputChain) {
    LOG(info) << "CbmPlutoGenerator: Closing input file " << fFileName;
    delete fInputChain;
  }
  fInputChain = nullptr;
}
// ------------------------------------------------------------------------
 
void CbmPlutoGenerator::CheckFileExist(std::string filename) {
  struct stat buffer;
  if (stat(filename.c_str(), &buffer) == 0) {
    LOG(info) << "CbmPlutoGenerator: Add file " << filename
              << " to input chain";
  } else {
    LOG(fatal) << "Input File " << filename << " not found";
  }
}
 
ClassImp(CbmPlutoGenerator)