Pythia8Generator.cxx

Go to the documentation of this file.

/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *         GNU Lesser General Public Licence version 3 (LGPL) version 3,        *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
// -------------------------------------------------------------------------
// -----                  M. Al-Turany   June 2014                     -----
// -------------------------------------------------------------------------
 
#include <FairPrimaryGenerator.h>
 
#include <TROOT.h>
 
#include "Pythia.h"
#include <math.h>
//#include <FairGenerator.h>
 
#include "Pythia8Generator.h"
 
using namespace Pythia8;
 
// -----   Default constructor   -------------------------------------------
Pythia8Generator::Pythia8Generator()
{
   fUseRandom1 = kFALSE;
   fUseRandom3 = kTRUE;
   fId = 2212; // proton
   fMom = 400; // proton
   fHNL = 0;   // HNL  if set to !=0, for example 9900014, only track
}
// -------------------------------------------------------------------------
 
// -----   Default constructor   -------------------------------------------
Bool_t Pythia8Generator::Init()
{
   if (fUseRandom1)
      fRandomEngine = new PyTr1Rng();
   if (fUseRandom3)
      fRandomEngine = new PyTr3Rng();
 
   fPythia.setRndmEnginePtr(fRandomEngine);
 
   cout << "Beam Momentum " << fMom << endl;
   fPythia.init(fId, 2212, 0., 0., fMom, 0., 0., 0.);
   return kTRUE;
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
Pythia8Generator::~Pythia8Generator() {}
// -------------------------------------------------------------------------
 
// -----   Passing the event   ---------------------------------------------
Bool_t Pythia8Generator::ReadEvent(FairPrimaryGenerator *cpg)
{
   Int_t npart = 0;
   while (npart == 0) {
      fPythia.next();
      for (int i = 0; i < fPythia.event.size(); i++) {
         if (fPythia.event[i].isFinal()) {
            // only send HNL decay products to G4
            if (fHNL != 0) {
               Int_t im = fPythia.event[i].mother1();
               if (fPythia.event[im].id() == fHNL) {
                  // for the moment, hardcode 110m is maximum decay length
                  Double_t z = fPythia.event[i].zProd();
                  Double_t x = abs(fPythia.event[i].xProd());
                  Double_t y = abs(fPythia.event[i].yProd());
                  // cout<<"debug HNL decay pos "<<x<<" "<< y<<" "<< z <<endl;
                  if (z < 11000. && z > 7000. && x < 250. && y < 250.) {
                     npart++;
                  }
               }
            } else {
               npart++;
            }
         };
      };
      // happens if a charm particle being produced which does decay without producing a HNL. Try another event.
      //       if (npart == 0){ fPythia.event.list();}
   };
   // cout<<"debug p8 event 0 " << fPythia.event[0].id()<< " "<< fPythia.event[1].id()<< " "<< fPythia.event[2].id()<< "
   // "<< npart <<endl;
   for (Int_t ii = 0; ii < fPythia.event.size(); ii++) {
      if (fPythia.event[ii].isFinal()) {
         Bool_t wanttracking = true;
         if (fHNL != 0) {
            Int_t im = fPythia.event[ii].mother1();
            if (fPythia.event[im].id() != fHNL) {
               wanttracking = false;
            }
         }
         if (wanttracking) {
            Double_t z = fPythia.event[ii].zProd();
            Double_t x = fPythia.event[ii].xProd();
            Double_t y = fPythia.event[ii].yProd();
            Double_t pz = fPythia.event[ii].pz();
            Double_t px = fPythia.event[ii].px();
            Double_t py = fPythia.event[ii].py();
            cpg->AddTrack((Int_t)fPythia.event[ii].id(), px, py, pz, x, y, z, (Int_t)fPythia.event[ii].mother1(),
                          wanttracking);
            // cout<<"debug p8->geant4 "<< wanttracking << " "<< ii <<  " " << fPythia.event[ii].id()<< " "<<
            // fPythia.event[ii].mother1()<<" "<<x<<" "<< y<<" "<< z <<endl;
         }
         //    virtual void AddTrack(Int_t pdgid, Double_t px, Double_t py, Double_t pz,
         //                          Double_t vx, Double_t vy, Double_t vz, Int_t parent=-1,Bool_t
         //                          wanttracking=true,Double_t e=-9e9);
      };
      if (fHNL != 0 && fPythia.event[ii].id() == fHNL) {
         Int_t im = (Int_t)fPythia.event[ii].mother1();
         Double_t z = fPythia.event[ii].zProd();
         Double_t x = fPythia.event[ii].xProd();
         Double_t y = fPythia.event[ii].yProd();
         Double_t pz = fPythia.event[ii].pz();
         Double_t px = fPythia.event[ii].px();
         Double_t py = fPythia.event[ii].py();
         cpg->AddTrack((Int_t)fPythia.event[im].id(), px, py, pz, x, y, z, 0, false);
         cpg->AddTrack((Int_t)fPythia.event[ii].id(), px, py, pz, x, y, z, im, false);
         // cout<<"debug p8->geant4 "<< 0 << " "<< ii <<  " " << fake<< " "<< fPythia.event[ii].mother1()<<endl;
      };
   }
 
   // make separate container ??
   //    FairRootManager *ioman =FairRootManager::Instance();
 
   return kTRUE;
}
// -------------------------------------------------------------------------
void Pythia8Generator::SetParameters(char *par)
{
   // Set Parameters
   fPythia.readString(par);
   cout << "fPythia.readString(\"" << par << "\")" << endl;
}
 
// -------------------------------------------------------------------------
void Pythia8Generator::Print()
{
   fPythia.settings.listAll();
}
// -------------------------------------------------------------------------
void Pythia8Generator::GetPythiaInstance(int arg)
{
   fPythia.particleData.list(arg);
   cout << "canDecay " << fPythia.particleData.canDecay(arg) << " " << fPythia.particleData.mayDecay(arg) << endl;
}
// -------------------------------------------------------------------------
 
ClassImp(Pythia8Generator)