source: trunk/examples/novice/gemc/src/MEventAction.cc@ 1323

// %%%%%%%%%%
// G4 headers
// %%%%%%%%%%
#include "G4Event.hh"

// %%%%%%%%%%%%%
// gemc headers
// %%%%%%%%%%%%%
#include "MEventAction.h"
#include "MHit.h"

#include <iostream>
using namespace std;

MEventAction::MEventAction(gemc_opts opts) {gemcOpt = opts;}
MEventAction::~MEventAction(){;}

void MEventAction::BeginOfEventAction(const G4Event* evt)
{
 string hd_msg  = gemcOpt.args["LOG_MSG"].args + " Event Action: >> ";
 int Modulo    = (int) gemcOpt.args["PRINT_EVENT"].arg ;

 if(evtN%Modulo == 0 )
    cout << hd_msg << " Begin of event " << evtN << endl;

}

void MEventAction::EndOfEventAction(const G4Event* evt)
{
 string hd_msg  = gemcOpt.args["LOG_MSG"].args + " Event Action: >> ";
 int Modulo    = (int) gemcOpt.args["PRINT_EVENT"].arg ;
 double    VERB = gemcOpt.args["OUT_VERBOSITY"].arg ;
 string catch_v = gemcOpt.args["CATCH"].args;


 if(evtN%Modulo == 0 )
    cout << hd_msg << "   End of event " << evtN << endl << endl;

 MHitCollection* MHC;
 MHit* aHit;
 int nhits;
 MPHBaseClass     *ProcessHitRoutine = NULL;
 MOutputBaseClass *ProcessOutput     = NULL;

 PH_output output;
 int bitcid;

 double Etot;

 // Making sure the output routine exists in the ProcessOutput map
 map<string, MOutput_Factory>::iterator ito = Out->find(MOut->outType);
 if(ito == Out->end())
 {
    if(MOut->outType != "no") 
       cout << hd_msg << " Warning: output type <" << MOut->outType << "> is not registered in the <MOutput_Factory>. " << endl
                      << "      This event will not be written out." << endl;
    return;
 }
 ProcessOutput     = GetMOutputClass((*Out), MOut->outType);

 // Header Bank contains event number
 // Need to change this to read DB header bank
 ProcessOutput->SetOutpHeader(evtN, MOut);

 // Getting Generated Particles info
 vector<MGeneratedParticle> MPrimaries;
 for(int pv=0; pv< evt->GetNumberOfPrimaryVertex(); pv++)
 {
    MGeneratedParticle Mparticle;
    G4PrimaryVertex* MPV = evt->GetPrimaryVertex(pv);
    Mparticle.vertex = MPV->GetPosition();
    for(int pp = 0; pp<MPV->GetNumberOfParticle(); pp++)
    {

      G4PrimaryParticle *PP = MPV->GetPrimary(pp);
      Mparticle.momentum    = PP->GetMomentum();
      Mparticle.PID         = PP->GetPDGcode();
    }
    MPrimaries.push_back(Mparticle)  ;
 }
 ProcessOutput-> WriteGenerated(MOut, MPrimaries);

 for(map<string, MSensitiveDetector*>::iterator it = SeDe_Map.begin(); it!= SeDe_Map.end(); it++)
 {
    string hitType;
    ProcessOutput->SetBankHeader(it->second->SDID.id, it->first, MOut);

    MHC = it->second->GetMHitCollection();
    nhits = MHC->GetSize();
    // The same ProcessHit Routine must apply to all the hits  in this HitCollection. 
    // Instantiating the ProcessHitRoutine only once for the first hit.
    if(nhits)
    {
       aHit = (*MHC)[0];
       string vname = aHit->GetDetector().name;
       hitType = it->second->GetDetectorHitType(vname);

       // Making sure the Routine exists in the MProcessHit_Map
       map<string, MPHB_Factory>::iterator itp = MProcessHit_Map->find(hitType);
       if(itp == MProcessHit_Map->end())
       {
          cout << hd_msg << " Warning: hit Type <" << hitType << "> is not registered in the <MProcessHit_Map>. " << endl
                         << "      This hit collection will not be processed." << endl;
         return;
       }
       ProcessHitRoutine = GetMPHClass((*MProcessHit_Map), hitType);

       for(int h=0; h<nhits; h++)
       {
          aHit = (*MHC)[h];
          output = ProcessHitRoutine->ProcessHit(aHit, gemcOpt);
          bitcid = it->second->get_bit_compressed_id(output.identity);
          ProcessOutput->ProcessOutput(bitcid, output, MOut, (*MBank_Map)[hitType]);

          string vname = aHit->GetId()[aHit->GetId().size()-1].name;
          if(VERB > 4 || vname.find(catch_v) != string::npos)
          {
             cout << hd_msg << " Hit " << h + 1 << " --  total number of steps this hit: " << aHit->GetPos().size() << endl;
             cout << aHit->GetId();
             Etot = 0;
             for(int e=0; e<aHit->GetPos().size(); e++) Etot = Etot + aHit->GetEdep()[e];
             cout << "   Total energy deposited: " << Etot/MeV << " MeV" << endl;
          }
       }
       //       cout << it->second->HCname << " " << evtN << " " << nhits << endl;
    }
    ProcessOutput->RecordAndClear(MOut, (*MBank_Map)[hitType]);

 }
 ProcessOutput->WriteEvent(MOut);

 delete ProcessOutput;
 delete ProcessHitRoutine;

 evtN++;    
}