LCDG4EventAction.cc

Go to the documentation of this file.

#include "LCDG4EventAction.hh"
#include "LCDG4TrackerHit.hh"
#include "LCDG4CalHit.hh"

#include "G4Event.hh"
#include "G4HCofThisEvent.hh"
#include "G4VVisManager.hh"
#include "G4SDManager.hh"
#include "G4UImanager.hh"
#include "G4ParticleTable.hh"

#include "Randomize.hh"
#include "CLHEP/Random/Random.h"
#include "CLHEP/Random/Randomize.h"
#include "Util_String.hh"

#include "SIO_stream.h"
#include "SIO_record.h"
#include "SIO_recordManager.h"
#include "SIO_blockManager.h"
#include "SIO_iodev.h"

#include "EVENT/LCIO.h"
#include "IMPL/LCCollectionVec.h"
#include "IMPL/SimCalorimeterHitImpl.h"
#include "IMPL/SimTrackerHitImpl.h"
#include "IMPL/LCFlagImpl.h"
using EVENT::LCIO;
using IMPL::LCCollectionVec;
using IMPL::SimCalorimeterHitImpl;
using IMPL::SimTrackerHitImpl;
using IMPL::LCFlagImpl;

#include <list>
#include <vector>
#include <iostream>
#include <fstream>
using std::vector;
using std::endl;
using std::ifstream;
using std::ofstream;
using std::ios;
using std::string;

// Norman's stuff
#include <vector>
#include "TkHit.h"
#include "MCParticleSIO.h"
#include "MCPrintSIO.h"
#include "VxdSIO.h"
#include "TrackerSIO.h"
#include "CalorimeterSIO.h"
#include "MuonHit.h"
#include "MuonSIO.h"
#include "LCDtowerID.h"
#include "LCDtowerIDNonProj.hh"

LCDG4EventAction::LCDG4EventAction(SIO_stream *stream_out, LCDG4McPartManager *mcmgr)
  : fixedSeed(0)
    , m_mcmgr(mcmgr)
    , _lcioHelper(NULL)
#ifdef TXT_MODE
    , outf_txt(NULL)
#endif
{
  m_vertexCollID = -9;
  m_trackerCollID = -9;
  m_ecalBCollID = -9;
  m_ecalECCollID = -9;
  m_hcalBCollID = -9;
  m_hcalECCollID = -9;
  m_muonCollID = -9;
  m_lumCollID = -9;
  m_event = 1;
  stream = stream_out;
  // create a new dump file
  ofstream dumpfile("debug-helper");
  dumpfile.close();
}

// This sets up LCIO output format
// File name is one of the arguments in command line
// Calling this (optional) method is all it takes to setup text output
void LCDG4EventAction::SetupOutputLcio(const char* filename)
{
  G4cout << "Setting up LCIO output file " << filename << G4endl;
  if (filename == "") {
    G4cerr<<"Argument for LCIO output file has not been set. Aborting.";
    exit(1);
  }
  _lcioHelper = LCDG4LcioHelper::getInstance();
  _lcioHelper->openLcioFile(string(filename));
}

#ifdef TXT_MODE
// This sets up ASCII output format
// File name is taken from TXTFILEPATH environment variable
// Calling this (optional) method is all it takes to setup text output
void LCDG4EventAction::SetupOutputTextFile()
{
  // ASCII output
  const char* fileenv = getenv("TXTFILEPATH");
  G4cout << "ASCII file from env: " << fileenv << G4endl;
  if (fileenv == "") {
    G4cerr<<"Env variable for ASCII output file has not been set. Aborting.";
    exit(1);
  }

  outf_txt = new ofstream(fileenv);
  if (outf_txt->fail() ) {
    G4cerr << "Failed to open output ASCII file. Aborting."<< G4endl;
    exit(1);
  }
}
#endif

LCDG4EventAction::~LCDG4EventAction() {
  if(_lcioHelper) delete _lcioHelper;
#ifdef TXT_MODE
  if(outf_txt) {
    outf_txt->close();
    delete outf_txt;
  }
#endif
}

void LCDG4EventAction::BeginOfEventAction(const G4Event* evt) {

  if(evt==0) return;

  // For debugging a specific event ==================
  // create file debevt.dat with just an event number inside it (start from 1).
  // This number will be read into ievtdgb.
  // debugging flag in McPartMgr will be set just for event chosen
  static int ievt=0;
  static int ievtdbg;
  if(ievt==0) {
    ifstream debevtfile("debevt.dat");
    debevtfile>> ievtdbg >> fixedSeed;
    debevtfile.close();
  }

  // increment event counter: 1st event is 1, not zero
  ievt++;
  if(ievt==ievtdbg) {
    m_mcmgr->SetDebugFlag(0x1f0);
    G4UImanager* UI = G4UImanager::GetUIpointer();
    UI->ApplyCommand("/tracking/verbose 1");
  }

  // At first event only
  if(ievt==1) {
    // Initialize the random number generator
//   HepRandom::createInstance();
//   HepRandomEngine* anEngine = new RanluxEngine((unsigned)time(NULL),3);
//   HepRandom::setTheEngine(anEngine);

    // Set the random seed for best randomness
    G4int seed = (unsigned)time(NULL);
    if(fixedSeed<0) {
      seed=-fixedSeed;
      G4cout<<"***** Using user provided seed at init only: seed="
            << seed << G4endl;
    }
    else {
      G4cout<<"***** Using time for best randomness: seed="  << seed << G4endl;
    }
    SetRandSeed((G4int)seed);
  }

  if(fixedSeed>0) {
    // Use a fixed random seed for all events: for easy reproducibility
    G4cout<< "***** Using fixed random seed = " << fixedSeed << G4endl;
    SetRandSeed(fixedSeed);
  }

  G4SDManager * SDman = G4SDManager::GetSDMpointer();
  G4String colNam;
  if(m_vertexCollID < -1) {
    m_vertexCollID = SDman->GetCollectionID(colNam="vxdCollection");
  }
  if(m_trackerCollID < -1) {
    m_trackerCollID = SDman->GetCollectionID(colNam="trackerCollection");
  }
  if(m_ecalBCollID < -1) {
    m_ecalBCollID = SDman->GetCollectionID(colNam="ecalBCollection");    
  }
  if(m_ecalECCollID < -1) {
    m_ecalECCollID = SDman->GetCollectionID(colNam="ecalECCollection");
  }
  if(m_hcalBCollID < -1) {
    m_hcalBCollID = SDman->GetCollectionID(colNam="hcalBCollection");
  }
  if(m_hcalECCollID < -1) {
    m_hcalECCollID = SDman->GetCollectionID(colNam="hcalECCollection");
  }
  if(m_muonCollID < -1) {
    m_muonCollID = SDman->GetCollectionID(colNam="muonCollection");
  }
  if(m_lumCollID < -1) {
    m_lumCollID = SDman->GetCollectionID(colNam="lumCollection");
  }

  G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
  if(pVVisManager) {
    G4UImanager::GetUIpointer()->ApplyCommand("/vis/drawTree/current");
  }
}

void LCDG4EventAction::EndOfEventAction(const G4Event* evt) {

  ofstream dumpfile("debug-helper",ios::app);

  G4HCofThisEvent*  HCE = evt->GetHCofThisEvent();
  LCDG4TrackerHitsCollection*     VHC = NULL;
  LCDG4TrackerHitsCollection* THC = NULL;
  LCDG4CalHitsCollection*     EBHC = NULL;
  LCDG4CalHitsCollection*     EECHC = NULL;
  LCDG4CalHitsCollection*     HBHC = NULL;
  LCDG4CalHitsCollection*     HECHC = NULL;
  LCDG4CalHitsCollection*    MHC = NULL;
  LCDG4CalHitsCollection*     LHC = NULL;

  if(HCE) {
    if (m_vertexCollID >= 0) {
      VHC = (LCDG4TrackerHitsCollection*)(HCE->GetHC(m_vertexCollID));
    }
    if (m_trackerCollID >= 0) {
      THC = (LCDG4TrackerHitsCollection*)(HCE->GetHC(m_trackerCollID));
    }
    if(m_ecalBCollID >= 0) {
      EBHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_ecalBCollID));
    }
    if(m_ecalECCollID >= 0) {
      EECHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_ecalECCollID));
    }
    if(m_hcalBCollID >= 0) {
      HBHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_hcalBCollID));
    }
    if(m_hcalECCollID >= 0) {
      HECHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_hcalECCollID));
    }
    if (m_muonCollID >= 0) {
      MHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_muonCollID));
    }
    if (m_lumCollID >= 0) {
      LHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_lumCollID));
    }
  }

#ifdef TXT_MODE
  if(outf_txt) {
    *outf_txt << "EVENT " << evt->GetEventID() << G4endl;
  }
#endif

  // SIO
  SIO_record *record = SIO_recordManager::add("LCD_LCDG4_Event");
  SIO_iodev *event = new SIO_iodev(evt->GetEventID(),"EventMarker");
  SIO_blockManager::add(event);
  record->connect(event);

  VxdSIO *vxd = NULL;
  TrackerSIO *tracker = NULL;
  CalorimeterSIO *ecalsio = NULL;
  CalorimeterSIO *hcalsio = NULL;
  CalorimeterSIO *lum = NULL;
  MuonSIO *muon = NULL;
  MCPrintSIO* mcprint = new MCPrintSIO("MCPrint");

  //**************  MC Particle Information  *************
  vector<LCDMcPart*>& mclist = m_mcmgr->MCparticles;
  vector<LCDMcPart*>::const_iterator mciter = mclist.begin();

  int mapindex = 0;  // tags start at one, not zero

  // sio format requires that first element is parent of all other elements.
  // Its ParticleID is set to 99999999 (HePEVT id).
  // Sometimes the input stdhep file already has this HepEVT element, then
  // keep it.  Otherwise, create a new one for the .sio output
  //  if( (*mciter)->GetParticleID()!=99999999 ) { // no, it did not work...
  MCParticleSIO* mcparthead = new MCParticleSIO(99999999,0,0);
  mcparthead->setfsflag(300);
  mcpartmap[mapindex] = mcparthead;
  mapindex++;
  //  }

  dumpfile <<"  "<< evt->GetEventID();
  dumpfile <<"  "<< mclist.size();
  for(mciter = mclist.begin(); mciter != mclist.end(); mciter++, mapindex++) {
    LCDMcPart* pMcPart = *mciter;
    if(!pMcPart) continue;  // some McParts have been erased!

    MCParticleSIO* mcp;
    mcp = new MCParticleSIO(pMcPart->GetParticleID(),(int)pMcPart->GetCharge(),pMcPart->Get4MomentumPtr()->e());
    mcp->setMcPart(pMcPart);

    MCParticleSIO::Status Status = MCParticleSIO::ALIVE;

    switch(pMcPart->GetG4Status()) {
    case 0 :
      Status = MCParticleSIO::ALIVE;  break;
    case 1 :
      Status = MCParticleSIO::DECAYED;  break;
    case 2 :
      Status = MCParticleSIO::INTERACTED;  break;
    case 3 :
      Status = MCParticleSIO::LEFT;  break;
    case 4 :
      Status = MCParticleSIO::STOPPED;  break;
    case 5 :
      Status = MCParticleSIO::LOOPING;  break;
    case 6 :
      Status = MCParticleSIO::LOST;  break;
    case 7 :
      Status = MCParticleSIO::STUCK;  break;
    case 8 :
      Status = MCParticleSIO::PRIMARY;  break;
    case 9 :
      Status = MCParticleSIO::SHOWER;  break;
    case 10 :
      Status = MCParticleSIO::MAXSTEPS;  break;
    }
    mcp->setfsflag( 100*pMcPart->GetGenStatus() );


    mcp->setStatus(Status);
    // McPart position is set to where it dies
    // starting position is taken from where its parent dies
    // secondaries from middle of track lose info on starting point
    mcp->setPosition(pMcPart->GetTermPosition().x(),
                     pMcPart->GetTermPosition().y(),
                     pMcPart->GetTermPosition().z());
    mcp->setMomentum(pMcPart->Get4Momentum().px(),
                     pMcPart->Get4Momentum().py(),
                     pMcPart->Get4Momentum().pz());
    mcp->setCalPosition(pMcPart->GetCalorimeterEntrancePoint().x(),
                        pMcPart->GetCalorimeterEntrancePoint().y(),
                        pMcPart->GetCalorimeterEntrancePoint().z() );
    mcp->setCalMomentum(pMcPart->GetCalorimeterEntranceMomentum().x(),
                        pMcPart->GetCalorimeterEntranceMomentum().y(),
                        pMcPart->GetCalorimeterEntranceMomentum().z() );

//     if(mapindex==1) G4cout<< "\n **** Array of LCDMcPart's ****\n";
//     G4cout << "*** tag="<< mapindex
//         <<", id="<< pMcPart->GetParticleID()
//         <<", status="<< pMcPart->GetGenStatus()
//         <<' '<< pMcPart->GetG4Status()
//         <<", q="<< (int)pMcPart->GetCharge()
//         <<", E="<< pMcPart->Get4Momentum().e()
//         <<", pvec=("<< pMcPart->Get4Momentum().px()
//         <<"; "<< pMcPart->Get4Momentum().py()
//         <<"; "<< pMcPart->Get4Momentum().pz()
//         <<") par="<< pMcPart->GetParentTag()
//         <<", startPos:x,y,z=("<< pMcPart->GetPositionPtr()->x()/cm
//         <<' '<< pMcPart->GetPositionPtr()->y()/cm
//         <<' '<< pMcPart->GetPositionPtr()->z()/cm <<')'
//         <<", endPos:x,y,z=("<< pMcPart->GetTermPositionPtr()->x()/cm
//         <<' '<< pMcPart->GetTermPositionPtr()->y()/cm
//         <<' '<< pMcPart->GetTermPositionPtr()->z()/cm <<')'
//         <<", x,y,z=("<< pMcPart->GetCalorimeterEntrancePoint().x()
//         <<' '<< pMcPart->GetCalorimeterEntrancePoint().y()
//         <<' '<< pMcPart->GetCalorimeterEntrancePoint().z() <<')'
//         <<", mom=("<< pMcPart->GetCalorimeterEntranceMomentum().x()
//         <<' '<< pMcPart->GetCalorimeterEntranceMomentum().y()
//         <<' '<< pMcPart->GetCalorimeterEntranceMomentum().z() <<')'
//         << G4endl;

    mcpartmap[mapindex] = mcp;
  }

  // Setup hierarchy links (parents <--> children) based on McParts
  for(unsigned int k = 1; k < mcpartmap.size(); k++) {
    MCParticleSIO* pmcp = mcpartmap[k];
    if(!pmcp) continue;

    const LCDMcPart* pMcPart = pmcp->getMcPart();

    if( pMcPart->GetParentTag()==-1) {
      mcpartmap[0]->addChild(pmcp);
      pmcp->setParent(mcpartmap[0]);
    }
    else {
      int parentTag = pMcPart->GetParentTag();
      MCParticleSIO* parent = mcpartmap[parentTag];
      parent->addChild(pmcp);
      pmcp->setParent(parent);

      // set parent status to DECAYED
      if(parent->GetG4Status()==MCParticleSIO::ALIVE) {
        // for LCDMcPart
        if(parentTag>0) mclist[parentTag-1]->SetG4Status(1);
        // and for MCParticleSIO
        parent->SetG4Status(MCParticleSIO::DECAYED);
      }
    }
  }

  // set other status codes
  mcpartmap[0]->setStatus( MCParticleSIO::PRIMARY );
  for(unsigned int k = 1; k < mcpartmap.size(); k++) {
    MCParticleSIO* pmcp = mcpartmap[k];
    if(!pmcp) continue;

    LCDMcPart* pMcPart = mclist[k-1];

    // Documentation particles
    if( pmcp->getfsflag() == 300 ) {
      // HEPEVT status DOCUMENTATION
      pmcp->setStatus( MCParticleSIO::PRIMARY );
      pMcPart->SetG4Status(8);
//        continue;
    }
    else if( pmcp->getfsflag() == 200 ) {
      // HEPEVT status INTERMEDIATE
      if( pmcp->status() == MCParticleSIO::ALIVE ) {
        // several gluons and quarks show up here
        pmcp->setStatus( MCParticleSIO::PRIMARY );
        pMcPart->SetG4Status(8);
      }
    }

//     LCDMcPart* lcdmcpart = const_cast<LCDMcPart*>(pmcp->getMcPart());
//     Hep3Vector& hep3v = lcdmcpart->GetPosition();
//     if(k==1) G4cout<< "\n **** Array of MCParticleSIO's ****\n";
//     G4cout << pmcp
//         <<", tag="<< k <<", id="<< pmcp->id()
//         <<", stat(Gen,Sim)="<< pmcp->GetGenStatus() <<','<< pmcp->GetG4Status()
//         <<", q="<< pmcp->charge()
//         <<", E="<< pmcp->energy()
//                 <<", xyzProd=("<< hep3v.x()/mm
//                 <<' '<< hep3v.y()/mm
//                 <<' '<< hep3v.z()/mm
//         <<", xyzEnd=("<< pmcp->x()/mm
//         <<' '<< pmcp->y()/mm
//         <<' '<< pmcp->z()/mm <<')'
//         <<", pvec=("<< pmcp->px()/GeV
//         <<' '<< pmcp->py()/GeV
//         <<' '<< pmcp->pz()/GeV <<')'
//         <<", Calxyz=("<< pmcp->calX()/mm
//         <<' '<< pmcp->calY()/mm
//         <<' '<< pmcp->calZ()/mm <<')'
//         <<", Calmom=("<< pmcp->calPX()/GeV
//         <<' '<< pmcp->calPY()/GeV
//         <<' '<< pmcp->calPZ()/GeV <<')'
//         <<", parID="<< pmcp->parent()->id()
//         <<", parE="<< pmcp->parent()->energy()/GeV
//         <<' '<< pmcp->parent()
//         << G4endl << G4endl;

//     G4cout << "daulist ("<< pmcp->numChildren() <<" kids)";
//     for(int ii=0; ii<pmcp->numChildren(); ii++) {
//       const MCParticleSIO* kid = pmcp->child(ii);
//       G4cout <<" : "<< kid->id() <<' '<< kid;
//     }
//     G4cout << G4endl;

      }

  // Setup MCParticle table for output to SIO file
  mcprint->set(mcpartmap[0]);

  // Build the particle table for LCIO
  if(_lcioHelper) _lcioHelper->buildMCParticleColl(mclist);


  //**************  VXD hits  *************
  if(VHC) {
    int n_hit = VHC->entries();
    dumpfile << "  "<< n_hit;
    vector<TkHit> vxdhits;

    LCCollectionVec* lcioVxdColl = new LCCollectionVec(LCIO::SIMTRACKERHIT);

    for(int ihit = 0; ihit < n_hit; ihit++) {
      LCDG4TrackerHit* g4vxd=(*VHC)[ihit];

      // LCIO hit handling
      SimTrackerHitImpl* lcioVxdHit=_lcioHelper->getTrkHitLcio(*g4vxd);
      lcioVxdColl->push_back(lcioVxdHit);

      double edep = g4vxd->GetEdep();
      float tdep = g4vxd->GetTdep();
      int   trkid = g4vxd->GetTrackID();
      assert(trkid>0);
      G4int mctag = m_mcmgr->trkidToMctag(trkid);
      unsigned int layn = g4vxd->GetLayer();
      unsigned int sysn = g4vxd->GetSystem();
      unsigned int beflg = g4vxd->GetBarEnd();
      unsigned int nsflg = g4vxd->GetNorthSouth();
      G4ThreeVector xpos = g4vxd->GetPos();

      double pos[3];
      pos[0] = xpos.x();
      pos[1] = xpos.y();
      pos[2] = xpos.z();

      // pack hit info into tag
      unsigned int tag = layn;
      tag |= (sysn  << 28);
      tag |= (nsflg << 30);
      tag |= (beflg << 31);

      // Get MCParticle associated to this hit
      // make sure this is not a backscattering from calorimeter
      MCParticleSIO *mcp = NULL;
      int showerFlag = m_mcmgr->GetShowerFlag(trkid);
      if(showerFlag<=1) {
        mcp = mcpartmap[mctag];
      }

      // save hit into hits array
      TkHit tk(tag,pos[0],pos[1],pos[2],edep,mcp,tdep);
      vxdhits.push_back(tk);
    } // for ihit

    vxd = new VxdSIO("VXD",vxdhits);
    SIO_blockManager::add(vxd);

    // LCIO format
    // set flag for long format (including position)
    LCFlagImpl chFlag(0) ;
    lcioVxdColl->setFlag( chFlag.getFlag() );
    _lcioHelper->saveCollection("VXDtCollection",lcioVxdColl);
  } // if VHC


  //**************  Tracker hits  *************
  if(THC) {
    int n_hit = THC->entries();
    dumpfile << "  "<< n_hit;
    int ihit = 0;
    vector<TkHit> trkhits;

    LCCollectionVec* lcioTrkColl = new LCCollectionVec(LCIO::SIMTRACKERHIT);

    for(ihit = 0; ihit < n_hit; ihit++) {
      LCDG4TrackerHit* g4trk=(*THC)[ihit];

      // LCIO hit handling
      SimTrackerHitImpl* lcioTrkHit = _lcioHelper->getTrkHitLcio(*g4trk);
      lcioTrkColl->push_back(lcioTrkHit);

      double edep = g4trk->GetEdep();
      double tdep = g4trk->GetTdep();
      int trkid = g4trk->GetTrackID();
      assert(trkid>0);
      G4int mctag = m_mcmgr->trkidToMctag(trkid);
      int layn = g4trk->GetLayer();
      unsigned int sysn = g4trk->GetSystem();
      unsigned int beflg = g4trk->GetBarEnd();
      unsigned int nsflg = g4trk->GetNorthSouth();
      G4ThreeVector xpos = g4trk->GetPos();

      // pack hit info into tag
      unsigned int tag = layn;
      tag |= (sysn  << 28);
      tag |= (nsflg << 30);
      tag |= (beflg << 31);

      // Get MCParticle associated to this hit
      // make sure this is not a backscattering from calorimeter
      MCParticleSIO *mcp = NULL;
      int showerFlag = m_mcmgr->GetShowerFlag(trkid);
      if(showerFlag<=1) {
        mcp = mcpartmap[mctag];
      }
      else {
        // this is a backscattering hit!
//      G4cout<< "backscattering hit: mcp="<< mcp << ", mctag="<< mctag
//            << ", trkid="<< trkid <<", pos="<< xpos << G4endl;
      }

      double pos[3];
      pos[0] = xpos.x();
      pos[1] = xpos.y();
      pos[2] = xpos.z();

      TkHit tkhit(tag,pos[0],pos[1],pos[2],edep,mcp,tdep);
      trkhits.push_back(tkhit);
    } // for
    tracker = new TrackerSIO("Tracker",trkhits);
    SIO_blockManager::add(tracker);

    // LCIO format
    // set flag for long format (including position)
    LCFlagImpl chFlag(0);
    lcioTrkColl->setFlag( chFlag.getFlag() );
    _lcioHelper->saveCollection("CENtCollection",lcioTrkColl);

  } // if THC

  // ************* EMCAL  **********************
  // EM Calorimeter Hits
#ifdef TXT_MODE
  if(outf_txt) *outf_txt << "EMCAL ";
#endif
  if(EBHC || EECHC) {
    // number of hits, Barrel + EndCap
    int n_hit =0;
    if(EBHC) {
      n_hit += EBHC->entries();
    }
    if(EECHC) {
      n_hit += EECHC->entries();
    }
    // write on dumpfile
    dumpfile << "  "<< n_hit;
#ifdef TXT_MODE
    // Write out number of hits for this event (ASCII format)
    if(outf_txt) *outf_txt << n_hit << endl;
#endif

    // collect barrel and endcap hits into a single collection
    vector<CalorimeterHit> calhits;
    LCCollectionVec* lcioEMColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
    if(EBHC) SaveCalorimeterHits( EBHC, calhits, lcioEMColl );
    if(EECHC) SaveCalorimeterHits( EECHC, calhits, lcioEMColl );

    // save hits into sio format
    ecalsio = new CalorimeterSIO("EMCalorimeter",calhits);
    SIO_blockManager::add(ecalsio);

    // LCIO format
    // set flag for long format (including position)
    // and PDG and cellid1
    LCFlagImpl chFlag(0) ;
    chFlag.setBit( LCIO::CHBIT_LONG );
    lcioEMColl->setFlag( chFlag.getFlag() );
    // save hits into lcio format
    _lcioHelper->saveCollection("EMcalCollection",lcioEMColl);

  } // if
#ifdef TXT_MODE
  else {
    if(outf_txt) *outf_txt << "0" << G4endl;
  }
#endif

  //**************  HDCal hits  *************
  // HAD Calorimeter Hits
#ifdef TXT_MODE
  if(outf_txt) *outf_txt << "HCAL ";
#endif
  if(HBHC || HECHC) {
    // number of hits, Barrel + EndCap
    int n_hit =0;
    if(HBHC) n_hit += HBHC->entries(); 
    if(HECHC) n_hit += HECHC->entries();

    // write on dumpfile
    dumpfile << "  "<< n_hit;
#ifdef TXT_MODE
    // Write out number of hits for this event (ASCII format)
    if(outf_txt) *outf_txt << n_hit << endl;
#endif

    // collect barrel and endcap hits into a single collection
    vector<CalorimeterHit> calhits;
    LCCollectionVec* lcioHADColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
    if(HBHC) SaveCalorimeterHits( HBHC, calhits, lcioHADColl );
    if(HECHC) SaveCalorimeterHits( HECHC, calhits, lcioHADColl );

    // save hits into sio format
    hcalsio = new CalorimeterSIO("HADCalorimeter",calhits);
    SIO_blockManager::add(hcalsio);

    // LCIO format
    // set flag for long format (including position)
    // and PDG and cellid1
    LCFlagImpl chFlag(0) ;
    chFlag.setBit( LCIO::CHBIT_LONG );
    lcioHADColl->setFlag( chFlag.getFlag() );
    _lcioHelper->saveCollection("HADcalCollection",lcioHADColl);

  } // if(HBHC || HECHC)
#ifdef TXT_MODE
  else {
     if(outf_txt) *outf_txt << "0" << endl;
  }
#endif

  // ******************** LUM Cal Hits ********************
  if(LHC) {
    int n_hit = LHC->entries();
    dumpfile << "  "<< n_hit;
    vector<CalorimeterHit> lumhits;

    LCCollectionVec* lcioLumColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);

    int ihit = 0;
    for(ihit = 0; ihit < n_hit; ihit++) {
      LCDG4CalHit* g4lum=(*LHC)[ihit];
      double edeps = g4lum->GetEdep();
      double edepa = g4lum->GetEdepAbs();

      int layn = g4lum->GetLayer();
      int sysn = g4lum->GetSystem();
//        int beflg = g4lum->GetBarEnd();
      int ithe = g4lum->GetTheta();
      int iphi = g4lum->GetPhi();
      int nmc = g4lum->GetNMC();

      if(nmc != 0) {
        // LCIO hit handling
        SimCalorimeterHitImpl* lcioHit = _lcioHelper->getCalHitLcio(*g4lum);
        lcioLumColl->push_back(lcioHit);

        // SIO hit handling
        LCDtowerID *towerID = new LCDtowerID(0);
        towerID->SetTheta(ithe);
        towerID->SetPhi(iphi);
        towerID->SetLayer(layn);
        towerID->SetSystem(sysn);

        CalorimeterHit lumhit(towerID->GetTag());
        int imc = 0;
        for(imc = 0; imc < nmc; imc++) {
          int trkid = g4lum->GetTrkID(imc);
          double emc = g4lum->GetEMC(imc);
          G4int mctag = m_mcmgr->trkidToMctag(trkid);
          if(mctag != -1) {
            MCParticleSIO *mcp = mcpartmap[mctag];
            lumhit.addMcHit(mcp,emc);
          }
          else {
            G4cout <<" EndOfEvent LUM problem: mctag=-1 for imc,trkid="
                   << imc <<' '<< trkid << G4endl;
            G4cout <<" ihit,ir,ith,iphi,nmc,edeps,abs = "
                   << ihit <<' '<< layn <<' '<< ithe <<' '<< iphi <<' '
                   << nmc <<' '<< edeps/MeV <<' '<< edepa/MeV << endl;
          }
        } // for nmc
        lumhits.push_back(lumhit);
        delete towerID;
      } // if nmc
    } // for hit
    lum = new CalorimeterSIO("LUMCalorimeter",lumhits);
    SIO_blockManager::add(lum);

    // LCIO format
    // set flag for long format (including position)
    LCFlagImpl chFlag(0) ;
    chFlag.setBit( LCIO::CHBIT_LONG );
    lcioLumColl->setFlag( chFlag.getFlag() );
    _lcioHelper->saveCollection("LUMcalCollection",lcioLumColl);

  } // if LUM Cal Hits


  // ******************** Muon Cal Hits ********************
  if(MHC) {
    int n_hit = MHC->entries();
    dumpfile << "  "<< n_hit;
    vector<MuonHit> calhits;

    LCCollectionVec* lcioMuonColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
    int ihit = 0;
    for(ihit = 0; ihit < n_hit; ihit++) {
      LCDG4CalHit *g4muon = (*MHC)[ihit];
      double edeps = g4muon->GetEdep();
      double edepa = g4muon->GetEdepAbs();

      int layn = g4muon->GetLayer();
//        int sysn = g4muon->GetSystem();
//        int beflg = g4muon->GetBarEnd();
      int ithe = g4muon->GetTheta();
      int iphi = g4muon->GetPhi();
      int nmc = g4muon->GetNMC();
      if (nmc != 0) {
        // LCIO hit handling
        SimCalorimeterHitImpl* lcioMuonHit=_lcioHelper->getCalHitLcio(*g4muon);
        lcioMuonColl->push_back(lcioMuonHit);

        // SIO hit handling
        LCDtowerID *towerID = new LCDtowerID(0);
        towerID->SetTheta(g4muon->GetTheta());
        towerID->SetPhi(g4muon->GetPhi());
        towerID->SetLayer(g4muon->GetLayer());
        towerID->SetSystem(g4muon->GetSystem());

        MuonHit calhit(towerID->GetTag(),towerID->GetTag());

        for(int imc = 0; imc < nmc; imc++) {
          int trkid = g4muon->GetTrkID(imc);
          double emc = g4muon->GetEMC(imc);
          G4int mctag = m_mcmgr->trkidToMctag(trkid);
          if(mctag != -1) {
            MCParticleSIO *mcp = mcpartmap[mctag];
            calhit.addMcHit(mcp,emc);
          }
          else {
            G4cout <<" EndOfEvent MUON problem: mctag=-1 for imc,trkid="
                   << imc <<' '<< trkid << G4endl;
            G4cout <<" ihit,ir,ith,iphi,nmc,edeps,abs = "
                   << ihit <<' '<< layn <<' '<< ithe <<' '<< iphi <<' '
                   << nmc <<' '<< edeps/MeV <<' '<< edepa/MeV << endl;
          }
        } // for imc<nmc
        calhits.push_back(calhit);
        delete towerID;
      } // if nmc!=0
    } // for geant4 hits
    muon = new MuonSIO("MUCalorimeter",calhits);
    SIO_blockManager::add(muon);

    // LCIO format
    // set flag for long format (including position)
    LCFlagImpl chFlag(0) ;
    chFlag.setBit( LCIO::CHBIT_LONG );
    lcioMuonColl->setFlag( chFlag.getFlag() );
    _lcioHelper->saveCollection("MUONcalCollection",lcioMuonColl);
  } // if MHC

  dumpfile << endl;
  dumpfile.close();

  // SIO
  stream->write(record);

  delete vxd;
  delete tracker;
  delete ecalsio;
  delete hcalsio;
  delete lum;
  delete muon;
  delete mcprint;
  delete event;

  SIO_recordManager::remove(record);

  int evnum = evt->GetEventID();
  int necal=0, nhcal=0;
  // ECal
  if(EBHC||EECHC) {
    if(EBHC) necal += EBHC->entries();
    if(EECHC) necal += EECHC->entries();
  }
  // HCal
  if(HBHC||HECHC) {
    if(HBHC) nhcal += HBHC->entries();
    if(HECHC) nhcal += HECHC->entries();
  }

  if( evnum%10 == 0 ) {
    G4cout<<" Event  #Hits: VXD  Tracker    ECal    HCal  Muon   Lum"
          << G4endl;
  }
  G4cout << myFormat( 5, evnum )
         << myFormat(13, ( VHC ? VHC->entries() : 0 ) )
         << myFormat( 9, ( THC ? THC->entries() : 0 ) )
         << myFormat( 8, necal )
         << myFormat( 8, nhcal )
         << myFormat( 6, ( MHC ? MHC->entries() : 0 ) )
         << myFormat( 6, ( LHC ? LHC->entries() : 0 ) )
         << G4endl;

  // Write out event to LCIO output
  _lcioHelper->writeEvent(evt);

  // Clear LCDG4McPartManager arrays
  m_mcmgr->Clear();

  // Clear mcpartmap collection
  for(unsigned int i=0; i<mcpartmap.size(); ++i) {
    if(mcpartmap[i]) delete mcpartmap[i];
  }
  mcpartmap.clear();
}

void LCDG4EventAction::SaveCalorimeterHits(
     const LCDG4CalHitsCollection* pHitColl,
     vector<CalorimeterHit>& calhits,
     LCCollectionVec* lcioColl) const
{
  // loop over hits
  for(int i=0; i < pHitColl->entries(); ++i) {
    LCDG4CalHit* g4hit=(*pHitColl)[i];

    double edeps = g4hit->GetEdep();
    double edepa = g4hit->GetEdepAbs();

    int layn = g4hit->GetLayer();
    int sysn = g4hit->GetSystem();
//        int beflg = g4hit->GetBarEnd();
    int ithe = g4hit->GetTheta();
    int iphi = g4hit->GetPhi();
    int nmc = g4hit->GetNMC();

#ifdef TXT_MODE
    // added emcal text output for NP JM030625
    G4ThreeVector pos = g4hit->GetPos();
    if(outf_txt) {
      *outf_txt << layn <<' '<< ithe <<' '<< iphi <<' '
                << pos.getX()/mm <<' '
                << pos.getY()/mm <<' '
                << pos.getZ()/mm <<' '
                << edeps/MeV <<' '<< edepa/MeV << endl; // edeps in MeV
      outf_txt->flush();
    }
#endif

    if(nmc != 0) { // remove hits with deposition only in absorber

      // LCIO hit handling
      if(lcioColl) {
        SimCalorimeterHitImpl* lcioHit = _lcioHelper->getCalHitLcio(*g4hit);
        lcioColl->push_back(lcioHit);
      }

      // SIO hit handling
      LCDtowerID *towerID = new LCDtowerID(0);
      towerID->SetTheta(ithe);
      towerID->SetPhi(iphi);
      towerID->SetLayer(layn);
      towerID->SetSystem(sysn);
//       G4cout<<"Cal hit: sysn="<< sysn <<", layn="<< layn
//        <<", ithe="<< ithe <<", iphi="<< iphi << G4endl;

      assert( (towerID->GetTag()&0x0fffffff) == (g4hit->GetCellID()&0x0fffffff) );

      CalorimeterHit calhit(towerID->GetTag());
      int imc = 0;
      for(imc = 0; imc < nmc; imc++) {
        int trkid = g4hit->GetTrkID(imc);
        double emc = g4hit->GetEMC(imc);
        G4int mctag = m_mcmgr->trkidToMctag(trkid);

        if(mctag != -1) {
          MCParticleSIO *mcp = mcpartmap.find(mctag)->second;
          calhit.addMcHit(mcp,emc);
//          G4cout << "ir,ith,iphi,PDGID = " << layn <<' '<< ithe <<' '<< iphi
//                 <<' '<< mcp->id()
//                 <<":"<< mctag
//                 << endl;
        }
        else {
          G4cout <<" EndOfEvent EM problem: mctag=-1 for imc,trkid="
                 << imc <<' '<< trkid << G4endl;
          G4cout <<" ihit,ir,ith,iphi,nmc,edeps,abs = "
                 << i <<' '<< layn <<' '<< ithe <<' '<< iphi <<' '
                 << nmc <<' '<< edeps/MeV <<' '<< edepa/MeV << G4endl;
        }
      } // for imc<nmc
      calhits.push_back(calhit);
      delete towerID;
    } // if(nmc != 0)
  } // for hit
}

---