MpdPlutoGenerator.cxx

Go to the documentation of this file.

// -------------------------------------------------------------------------
// -----                MpdPlutoGenerator source file                  -----
// -----          Created 13/07/04  by V. Friese / D.Bertini           -----
// -----      Modified from FairPlutoGenerator for MPD by V. Zhezher   -----
// -------------------------------------------------------------------------
#include "MpdPlutoGenerator.h"
 
#include "FairPrimaryGenerator.h"
#include "PParticle.h"
#include "TClonesArray.h"
#include "TDatabasePDG.h"
#include "TFile.h"
#include "TLorentzVector.h"
#include "TTree.h"
#include "TVector3.h"
 
#include <iostream>
 
const Double_t kProtonMass = 0.938271998;
const Double_t kElectronMass = 0.00051098892;
 
// -----   Default constructor   ------------------------------------------
MpdPlutoGenerator::MpdPlutoGenerator()
    : FairGenerator()
    , iEvent(0)
    , fEkin(0)
    , fFileName(NULL)
    , fInputFile(NULL)
    , fInputTree(NULL)
    , fParticles(NULL)
{}
// ------------------------------------------------------------------------
 
// -----   Standard constructor   -----------------------------------------
MpdPlutoGenerator::MpdPlutoGenerator(const Char_t* fileName, Double_t ekin = 25.0)
    : FairGenerator()
    , iEvent(0)
    , fFileName(fileName)
    , fInputFile(new TFile(fileName))
    , fInputTree(NULL)
    , fParticles(new TClonesArray("PParticle", 100))
{
    fEkin = ekin;
    fInputTree = (TTree*)fInputFile->Get("data");
    fInputTree->SetBranchAddress("Particles", &fParticles);
}
// ------------------------------------------------------------------------
 
// -----   Destructor   ---------------------------------------------------
MpdPlutoGenerator::~MpdPlutoGenerator()
{
    CloseInput();
}
// ------------------------------------------------------------------------
 
// -----   Public method ReadEvent   --------------------------------------
Bool_t MpdPlutoGenerator::ReadEvent(FairPrimaryGenerator* primGen)
{
 
    // Check for input file
    if (!fInputFile) {
        cout << "-E MpdPlutoGenerator: Input file nor open!" << endl;
        return kFALSE;
    }
 
    // Check for number of events in input file
    if (iEvent > fInputTree->GetEntries()) {
        cout << "-E MpdPlutoGenerator: No more events in input file!" << endl;
        CloseInput();
        return kFALSE;
    }
    TFile* g = gFile;
    fInputFile->cd();
    fInputTree->GetEntry(iEvent++);
    g->cd();
 
    // 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 = dataBase->ConvertGeant3ToPdg(part->ID());
 
        // Check if particle type is known to database
        if (!pdgType) {
            cout << "-W MpdPlutoGenerator: Unknown type " << part->ID() << ", skipping particle." << endl;
            continue;
        }
 
        TLorentzVector mom = part->Vect4();
        Double_t px = mom.Px();
        Double_t py = mom.Py();
        Double_t pz = mom.Pz();
 
        TVector3 vertex = part->getVertex();
        Double_t vx = vertex.x();
        Double_t vy = vertex.y();
        Double_t vz = vertex.z();
 
        // ---> Calculate beta and gamma for Lorentztransformation
        Double_t eBeam = fEkin + kProtonMass;
        Double_t pBeam = TMath::Sqrt(eBeam * eBeam - kProtonMass * kProtonMass);
        Double_t betaCM = pBeam / (eBeam + kProtonMass);
        Double_t gammaCM = TMath::Sqrt(1. / (1. - betaCM * betaCM));
 
        // AZ Double_t mass = kElectronMass;
        Double_t mass = dataBase->GetParticle(pdgType)->Mass();
        Double_t e = sqrt(mass * mass + px * px + py * py + pz * pz);
        pz = gammaCM * (pz - betaCM * e);
 
        // Give track to PrimaryGenerator
        // primGen->AddTrack(pdgType, px, py, pz, vx, vy, vz);
        // AZ - keep track of parents
        Int_t parentIndx = part->GetParentIndex();
        primGen->AddTrack(pdgType, px, py, pz, vx, vy, vz, -parentIndx - 10);
        // AZ
 
    }   //  Loop over particle in event
 
    return kTRUE;
}
// ------------------------------------------------------------------------
 
// -----   Public method SkipEvents   --------------------------------------
 
Bool_t MpdPlutoGenerator::SkipEvents(Int_t count)
{
    if (count <= 0)
        return kTRUE;
 
    iEvent = count;
 
    std::cout << "-I MpdPlutoGenerator: Skipped " << count << " Events!" << std::endl;
 
    return kTRUE;
}
 
// -----   Private method CloseInput   ------------------------------------
void MpdPlutoGenerator::CloseInput()
{
    if (fInputFile) {
        cout << "-I MpdPlutoGenerator: Closing input file " << fFileName << endl;
        fInputFile->Close();
        delete fInputFile;
    }
    fInputFile = NULL;
}
// ------------------------------------------------------------------------