CastorShowerLibraryMaker.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     ShowerLibraryProducer
// Class  :     CastorShowerLibraryMaker
//
// Implementation:
//     <Notes on implementation>
//
// Original Author: P. Katsas
//         Created: 02/2007 
//
// Adapted by W. Carvalho , 02/2009
//
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "SimG4Core/Notification/interface/SimG4Exception.h"
#include "SimG4CMS/Calo/interface/CaloG4Hit.h"
#include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
#include "DataFormats/Math/interface/Point3D.h"

#include "SimG4CMS/ShowerLibraryProducer/interface/CastorShowerLibraryMaker.h"
#include "SimG4CMS/Forward/interface/CastorNumberingScheme.h"

#include "SimDataFormats/CaloHit/interface/CastorShowerEvent.h"

#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/GlobalPhysicalConstants.h"

#include "TFile.h"
#include <cmath>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <stdlib.h>

CastorShowerLibraryMaker::CastorShowerLibraryMaker(const edm::ParameterSet &p) : 
                             NPGParticle(0),DoHadSL(false),DoEmSL(false),DeActivatePhysicsProcess(false),
                             emShower(NULL) , hadShower(NULL) {

  MapOfSecondaries.clear();
  hadInfo = NULL;
  emInfo  = NULL;
  edm::ParameterSet p_SLM   = p.getParameter<edm::ParameterSet>("CastorShowerLibraryMaker");
  verbosity                 = p_SLM.getParameter<int>("Verbosity");
  eventNtFileName           = p_SLM.getParameter<std::string>("EventNtupleFileName");
  hadSLHolder.nEvtPerBinPhi = p_SLM.getParameter<int>("nhadEvents");
  emSLHolder.nEvtPerBinPhi  = p_SLM.getParameter<int>("nemEvents");
  hadSLHolder.SLEnergyBins  = p_SLM.getParameter<std::vector<double> >("SLhadEnergyBins");
  hadSLHolder.SLEtaBins     = p_SLM.getParameter<std::vector<double> >("SLhadEtaBins");
  hadSLHolder.SLPhiBins     = p_SLM.getParameter<std::vector<double> >("SLhadPhiBins");
  emSLHolder.SLEnergyBins   = p_SLM.getParameter<std::vector<double> >("SLemEnergyBins");
  emSLHolder.SLEtaBins      = p_SLM.getParameter<std::vector<double> >("SLemEtaBins");
  emSLHolder.SLPhiBins      = p_SLM.getParameter<std::vector<double> >("SLemPhiBins");
  PGParticleIDs             = p_SLM.getParameter<std::vector<int> >("PartID");
  DeActivatePhysicsProcess  = p_SLM.getParameter<bool>("DeActivatePhysicsProcess");
  MaxPhi                    = p_SLM.getParameter<double>("SLMaxPhi");
  MaxEta                    = p_SLM.getParameter<double>("SLMaxEta");
//
  NPGParticle               = PGParticleIDs.size(); 
  for(unsigned int i=0;i<PGParticleIDs.size();i++) {
     switch (int(fabs(PGParticleIDs.at(i)))) {
        case 11:
        case 22:
               DoEmSL = true;
               break;
        default:
               DoHadSL = true;
     }
  }
  hadSLHolder.nEvtPerBinEta = (hadSLHolder.nEvtPerBinPhi)*(hadSLHolder.SLPhiBins.size());
  hadSLHolder.nEvtPerBinE   = (hadSLHolder.nEvtPerBinEta)*(hadSLHolder.SLEtaBins.size());
  emSLHolder.nEvtPerBinEta = (emSLHolder.nEvtPerBinPhi)*(emSLHolder.SLPhiBins.size());
  emSLHolder.nEvtPerBinE   = (emSLHolder.nEvtPerBinEta)*(emSLHolder.SLEtaBins.size());

  std::cout << "============================================================================"<<std::endl;
  std::cout << "CastorShowerLibraryMaker:: Initialized as observer" << std::endl;
  std::cout << " Event Ntuple will be created" << std::endl;
  std::cout << " Event Ntuple file: " << eventNtFileName << std::endl;
  std::cout << " Number of Hadronic events in E   bins: " << hadSLHolder.nEvtPerBinE << std::endl;
  std::cout << " Number of Hadronic events in Eta bins: " << hadSLHolder.nEvtPerBinEta << std::endl;
  std::cout << " Number of Hadronic events in Phi bins: " << hadSLHolder.nEvtPerBinPhi << std::endl;
  std::cout << " Number of Electromag. events in E   bins: " << emSLHolder.nEvtPerBinE << std::endl;
  std::cout << " Number of Electromag. events in Eta bins: " << emSLHolder.nEvtPerBinEta << std::endl;
  std::cout << " Number of Electromag. events in Phi bins: " << emSLHolder.nEvtPerBinPhi << std::endl;
  std::cout << "============================================================================"<<std::endl;
  std::cout << std::endl;

  // Initializing the SL collections
  InitSLHolder(hadSLHolder);
  InitSLHolder(emSLHolder);
}
void CastorShowerLibraryMaker::InitSLHolder(ShowerLib& showerholder)
{
  int nBinsE,nBinsEta,nBinsPhi,nEvtPerBinPhi;
  nBinsE      = showerholder.SLEnergyBins.size();
  nBinsEta    = showerholder.SLEtaBins.size();
  nBinsPhi    = showerholder.SLPhiBins.size();
  nEvtPerBinPhi=showerholder.nEvtPerBinPhi;
//
// Info
//
  showerholder.SLInfo.Energy.setNEvts(nEvtPerBinPhi*nBinsPhi*nBinsEta*nBinsE);
  showerholder.SLInfo.Energy.setNEvtPerBin(nEvtPerBinPhi*nBinsPhi*nBinsEta);
  showerholder.SLInfo.Energy.setNBins(nBinsE);
  showerholder.SLInfo.Energy.setBin(showerholder.SLEnergyBins);
//
  showerholder.SLInfo.Eta.setNEvts(nEvtPerBinPhi*nBinsPhi*nBinsEta);
  showerholder.SLInfo.Eta.setNEvtPerBin(nEvtPerBinPhi*nBinsPhi);
  showerholder.SLInfo.Eta.setNBins(nBinsEta);
  showerholder.SLInfo.Eta.setBin(showerholder.SLEtaBins);
//
  showerholder.SLInfo.Phi.setNEvts(nEvtPerBinPhi*nBinsPhi);
  showerholder.SLInfo.Phi.setNEvtPerBin(nEvtPerBinPhi);
  showerholder.SLInfo.Phi.setNBins(nBinsPhi);
  showerholder.SLInfo.Phi.setBin(showerholder.SLPhiBins);
//
// Shower
  showerholder.SLCollection.assign(nBinsE,std::vector<std::vector<std::vector<CastorShowerEvent> > >());
  showerholder.nEvtInBinE.assign(nBinsE,0);
  showerholder.nEvtInBinEta.assign(nBinsE,std::vector<int>(0));
  showerholder.nEvtInBinPhi.assign(nBinsE,std::vector<std::vector<int> >());
  for(int i=0;i<nBinsE;i++) {
     showerholder.SLCollection.at(i).assign(nBinsEta,std::vector<std::vector<CastorShowerEvent> >());
     showerholder.nEvtInBinEta.at(i).assign(nBinsEta,0);
     showerholder.nEvtInBinPhi.at(i).assign(nBinsEta,std::vector<int>(0));
     for(int j=0;j<nBinsEta;j++) {
        showerholder.SLCollection.at(i).at(j).assign(nBinsPhi,std::vector<CastorShowerEvent>());
        showerholder.nEvtInBinPhi.at(i).at(j).assign(nBinsPhi,0);
        for(int k=0;k<nBinsPhi;k++) 
           showerholder.SLCollection.at(i).at(j).at(k).assign(nEvtPerBinPhi,CastorShowerEvent());
     }
  }
}

//===============================================================================================

CastorShowerLibraryMaker::~CastorShowerLibraryMaker() {

  Finish();

  std::cout << "CastorShowerLibraryMaker: End of process" << std::endl;

}

//=================================================================== per EVENT
void CastorShowerLibraryMaker::update(const BeginOfJob * job) {

  std::cout << " CastorShowerLibraryMaker::Starting new job " << std::endl;
}

//==================================================================== per RUN
void CastorShowerLibraryMaker::update(const BeginOfRun * run) {

  std::cout << std::endl << "CastorShowerLibraryMaker: Starting Run"<< std::endl; 

  std::cout << "CastorShowerLibraryMaker: output event root file created" << std::endl;

  TString eventfilename = eventNtFileName;
  theFile = new TFile(eventfilename,"RECREATE");
  theTree = new TTree("CastorCherenkovPhotons", "Cherenkov Photons");

  Int_t split = 1;
  Int_t bsize = 64000;
  emInfo      = new CastorShowerLibraryInfo();
  emShower    = new CastorShowerEvent();
  hadInfo     = new CastorShowerLibraryInfo();
  hadShower   = new CastorShowerEvent();
  // Create Branchs
  theTree->Branch("emShowerLibInfo.", "CastorShowerLibraryInfo", &emInfo, bsize, split);
  theTree->Branch("emParticles.", "CastorShowerEvent", &emShower, bsize, split);
  theTree->Branch("hadShowerLibInfo.", "CastorShowerLibraryInfo", &hadInfo, bsize, split);
  theTree->Branch("hadParticles.", "CastorShowerEvent", &hadShower, bsize, split);

// set the Info for electromagnetic shower
// set the energy bins info
  emInfo->Energy.setNEvts(emSLHolder.nEvtPerBinE*emSLHolder.SLEnergyBins.size());
  emInfo->Energy.setNBins(emSLHolder.SLEnergyBins.size());
  emInfo->Energy.setNEvtPerBin(emSLHolder.nEvtPerBinE);
  emInfo->Energy.setBin(emSLHolder.SLEnergyBins);
// set the eta bins info
  emInfo->Eta.setNEvts(emSLHolder.nEvtPerBinEta*emSLHolder.SLEtaBins.size());
  emInfo->Eta.setNBins(emSLHolder.SLEtaBins.size());
  emInfo->Eta.setNEvtPerBin(emSLHolder.nEvtPerBinEta);
  emInfo->Eta.setBin(emSLHolder.SLEtaBins);
// set the eta bins info
  emInfo->Phi.setNEvts(emSLHolder.nEvtPerBinPhi*emSLHolder.SLPhiBins.size());
  emInfo->Phi.setNBins(emSLHolder.SLPhiBins.size());
  emInfo->Phi.setNEvtPerBin(emSLHolder.nEvtPerBinPhi);
  emInfo->Phi.setBin(emSLHolder.SLPhiBins);
// The same for the hadronic shower
// set the energy bins info
  hadInfo->Energy.setNEvts(hadSLHolder.nEvtPerBinE*hadSLHolder.SLEnergyBins.size());
  hadInfo->Energy.setNBins(hadSLHolder.SLEnergyBins.size());
  hadInfo->Energy.setNEvtPerBin(hadSLHolder.nEvtPerBinE);
  hadInfo->Energy.setBin(hadSLHolder.SLEnergyBins);
// set the eta bins info
  hadInfo->Eta.setNEvts(hadSLHolder.nEvtPerBinEta*hadSLHolder.SLEtaBins.size());
  hadInfo->Eta.setNBins(hadSLHolder.SLEtaBins.size());
  hadInfo->Eta.setNEvtPerBin(hadSLHolder.nEvtPerBinEta);
  hadInfo->Eta.setBin(hadSLHolder.SLEtaBins);
// set the eta bins info
  hadInfo->Phi.setNEvts(hadSLHolder.nEvtPerBinPhi*hadSLHolder.SLPhiBins.size());
  hadInfo->Phi.setNBins(hadSLHolder.SLPhiBins.size());
  hadInfo->Phi.setNEvtPerBin(hadSLHolder.nEvtPerBinPhi);
  hadInfo->Phi.setBin(hadSLHolder.SLPhiBins);
  // int flag = theTree->GetBranch("CastorShowerLibInfo")->Fill();
  // Loop on all leaves of this branch to fill Basket buffer.
  // The function returns the number of bytes committed to the memory basket.
  // If a write error occurs, the number of bytes returned is -1.
  // If no data are written, because e.g. the branch is disabled,
  // the number of bytes returned is 0.
  // if(flag==-1) {
  //    edm::LogInfo("CastorAnalyzer") << " WARNING: Error writing to Branch \"CastorShowerLibInfo\" \n" ; 
  // } else 
  // if(flag==0) {
  //    edm::LogInfo("CastorAnalyzer") << " WARNING: No data written to Branch \"CastorShowerLibInfo\" \n" ; 
  // }

  // Initialize "accounting" variables

  eventIndex = 0;

}

//=================================================================== per EVENT
void CastorShowerLibraryMaker::update(const BeginOfEvent * evt) {

  eventIndex++;
  stepIndex = 0;
  InsideCastor = false;
  PrimaryMomentum.clear();
  PrimaryPosition.clear();
  int NAccepted =0;
// reset the pointers to the shower objects
  SLShowerptr = NULL;
  MapOfSecondaries.clear();
  thePrims.clear();
  G4EventManager *e_mgr = G4EventManager::GetEventManager();
  if (IsSLReady()) {
      printSLstatus(-1,-1,-1);
      update((EndOfRun*)NULL);
      return;
  }

  thePrims = GetPrimary((*evt)());
  for(unsigned int i=0;i<thePrims.size();i++) {
     G4PrimaryParticle* thePrim = thePrims.at(i);
     int particleType = thePrim->GetPDGcode();

     std::string SLType("");
     if (particleType==11) {
        SLShowerptr = &emSLHolder;
        SLType = "Electromagnetic";
     }
     else {
        SLShowerptr = &hadSLHolder;
        SLType = "Hadronic";
     }
     double px=0., py=0., pz=0., pInit = 0., eta = 0., phi = 0.;
     GetKinematics(thePrim,px,py,pz,pInit,eta,phi);
     int ebin = FindEnergyBin(pInit);
     int etabin= FindEtaBin(eta);
     int phibin = FindPhiBin(phi);
     if (verbosity)
        std::cout << "\n========================================================================"
                  << "\nBeginOfEvent: E   : " << pInit <<"\t  bin : " << ebin
                  << "\n              Eta : " << eta   <<"\t  bin : " << etabin
                  << "\n              Phi : " << phi   <<"\t  bin : " << phibin
                  << "\n========================================================================"
                  << std::endl;

     if (ebin<0||etabin<0||phibin<0) continue;
     bool accept=false;
     if (!(SLacceptEvent(ebin,etabin,phibin))) {
/*
// To increase the chance of a particle arriving at CASTOR inside a not full bin,
// check if there is available phase space in the neighboring bins
        unsigned int ebin_min = std::max(0,ebin-3);
        unsigned int eta_bin_min = std::max(0,etabin-2);
        unsigned int eta_bin_max = std::min(etabin,etabin+2);
        unsigned int phi_bin_min = std::max(0,phibin-2);
        unsigned int phi_bin_max = std::min(phibin,phibin+2);
        for(unsigned int i_ebin=ebin_min;i_ebin<=(unsigned int)ebin;i_ebin++) {
           for (unsigned int i_etabin=eta_bin_min;i_etabin<=eta_bin_max;i_etabin++) {
               for (unsigned int i_phibin=phi_bin_min;i_phibin<=phi_bin_max;i_phibin++) {
                   if (SLacceptEvent((int)i_ebin,(int)i_etabin,(int)i_phibin)) {accept=true;break;}
               }
               if (accept) break;
           }
           if (accept) break;
        }
*/
        if (!accept) edm::LogInfo("CastorShowerLibraryMaker") << "Event not accepted for ebin="
             << ebin<<",etabin="<<etabin<<",phibin="<<phibin<<std::endl;
     }
     else {
        accept=true;
     }
     if (accept) NAccepted++;
  }
  
  if (NAccepted==0) {
     const_cast<G4Event*>((*evt)())->SetEventAborted();
     const_cast<G4Event*>((*evt)())->KeepTheEvent((G4bool)false);
     e_mgr->AbortCurrentEvent();
  }
  SLShowerptr=NULL;
//
  std::cout << "CastorShowerLibraryMaker: Processing Event Number: " << eventIndex << std::endl;
}

//=================================================================== per STEP
void CastorShowerLibraryMaker::update(const G4Step * aStep) {
   static thread_local int CurrentPrimary = 0;
   G4Track *trk = aStep->GetTrack();
   if (trk->GetCurrentStepNumber()==1) {
      if (trk->GetParentID()==0) {
         CurrentPrimary = (int)trk->GetDynamicParticle()->GetPDGcode();
         if (CurrentPrimary==0) 
             SimG4Exception("CastorShowerLibraryMaker::update(G4Step) -> Primary particle undefined");
         InsideCastor=false;
// Deactivate the physics process
         if (DeActivatePhysicsProcess) {
            G4ProcessManager *p_mgr = trk->GetDefinition()->GetProcessManager();
            G4ProcessVector *pvec = p_mgr->GetProcessList();
            for (int i = 0;i<pvec->size();i++) {
                G4VProcess *proc = (*pvec)(i);
                if (proc->GetProcessName()!="Transportation"&&
                    proc->GetProcessName()!="Decay") {
                    std::cout << "DeActivating process: " << proc->GetProcessName() << std::endl;
                    p_mgr->SetProcessActivation(proc,false);
                }
            } 
         }
// move track to z of CASTOR
         G4ThreeVector pos;
         pos.setZ(-14390);
         double t = abs((pos.z()-trk->GetPosition().z()))/trk->GetVelocity();
         double r = (pos.z()-trk->GetPosition().z())/trk->GetMomentum().cosTheta();
         pos.setX(r*sin(trk->GetMomentum().theta())*cos(trk->GetMomentum().phi())+trk->GetPosition().x());
         pos.setY(r*sin(trk->GetMomentum().theta())*sin(trk->GetMomentum().phi())+trk->GetPosition().y());
         trk->SetPosition(pos);
         trk->SetGlobalTime(trk->GetGlobalTime()+t);
         trk->AddTrackLength(r);
      }
      else if (!InsideCastor) {
         std::cout<<"CastorShowerLibraryMaker::update(G4Step) -> Killing spurious track" << std::endl;
         trk->SetTrackStatus(fKillTrackAndSecondaries);
         return;
      }
      MapOfSecondaries[CurrentPrimary].insert((int)trk->GetTrackID());
   }
// Checks if primary already inside CASTOR
   std::string CurVolume = trk->GetVolume()->GetName();
   if (!InsideCastor&&(
       //CurVolume=="C3EF"||CurVolume=="C4EF"||CurVolume=="CAEL"||
       //CurVolume=="CAHL"||CurVolume=="C3HF"||CurVolume=="C4HF")) {
       //CurVolume=="CastorB"||
       CurVolume=="CAST")) {
       //CurVolume=="CAIR")) {
       InsideCastor=true;
// Activate the physics process
       if (trk->GetParentID()==0&&DeActivatePhysicsProcess) {
          G4ProcessManager *p_mgr = trk->GetDefinition()->GetProcessManager();
          G4ProcessVector *pvec = p_mgr->GetProcessList();
          for (int i = 0;i<pvec->size();i++) {
              G4VProcess *proc = (*pvec)(i);
              if (proc->GetProcessName()!="Transportation"&&
                 proc->GetProcessName()!="Decay") {
                 std::cout << "  Activating process: " << proc->GetProcessName() << std::endl;
                 p_mgr->SetProcessActivation(proc,true);
              }
          }
       }
       //PrimaryMomentum[CurrentPrimary]=aStep->GetPreStepPoint()->GetMomentum();
// check fiducial eta and phi
       if (trk->GetMomentum().phi()>MaxPhi||trk->GetMomentum().eta()>MaxEta) {
          trk->SetTrackStatus(fKillTrackAndSecondaries);
          InsideCastor=false;
          return;
       }
       PrimaryMomentum[CurrentPrimary]=trk->GetMomentum();
       PrimaryPosition[CurrentPrimary]=trk->GetPosition();
       KillSecondaries(aStep);
       return;
   } 
// Kill the secondaries if they have been produced before entering castor
   if (CurrentPrimary!=0&&trk->GetParentID()==0&&!InsideCastor) {
      KillSecondaries(aStep);
      if (verbosity) {
         double pre_phi = aStep->GetPreStepPoint()->GetMomentum().phi();
         double cur_phi = trk->GetMomentum().phi();
         if (pre_phi!=cur_phi) {
            std::cout << "Primary track phi :  " << pre_phi << " changed in current step: "
                      << cur_phi << " by processes:" << std::endl;
            const G4VProcess *proc = aStep->GetPreStepPoint()->GetProcessDefinedStep();
            std::cout << "           " << proc->GetProcessName()
                      << "  In volume " << CurVolume << std::endl;
         }
      }
   }

//==============================================
/*
*/
/*
  if(aStep->IsFirstStepInVolume()) { 
    edm::LogInfo("CastorShowerLibraryMaker") << "CastorShowerLibraryMaker::update(const G4Step * aStep):"
                                             << "\n IsFirstStepInVolume , " 
                                             << "time = " << aStep->GetTrack()->GetGlobalTime() ; 
  }
  stepIndex++;
*/
}

//================= End of EVENT ===============
void CastorShowerLibraryMaker::update(const EndOfEvent * evt) {

// check if the job is done!
  if ((*evt)()->IsAborted()) {
     std::cout << "\n========================================================================"
               << "\nEndOfEvent: EVENT ABORTED"
               << "\n========================================================================"
               << std::endl;
     return;
  }
  //DynamicRangeFlatRandomEGunProducer* pgun = edm::DynamicRangeFlatRandomEGunKernel::get_instance();
  //std::cout << pgun->EGunMaxE() << std::endl;
/*
  std::cout << "Minimum energy in Particle Gun : " << pgun->EGunMinE() << "\n"
            << "Maximum energy in Particle Gun : " << pgun->EGunMaxE() << std::endl;
*/
  if (verbosity) 
     std::cout << "CastorShowerLibraryMaker: End of Event: " << eventIndex << std::endl;
// Get the pointer to the primary particle
  if (thePrims.size() == 0) {
     edm::LogInfo("CastorShowerLibraryMaker") << "No valid primary particle found. Skipping event" << std::endl;
     return;
  }
  // access to the G4 hit collections 
  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
  int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
  CaloG4HitCollection* theCAFI = (CaloG4HitCollection*) allHC->GetHC(CAFIid);
  if (verbosity) edm::LogInfo("CastorShowerLibraryMaker") << " update(*evt) --> accessed all HC ";
  edm::LogInfo("CastorShowerLibraryMaker") << "Found "<<theCAFI->entries() << " hits in G4HitCollection";
  if (theCAFI->entries()== 0) {
     edm::LogInfo("CastorShowerLibraryMaker") << "\n Empty G4HitCollection";
     return;
  }

// Loop over primaries
  int NEvtAccepted = 0;
  int NHitInEvent  = 0;
  for(unsigned int i=0;i<thePrims.size();i++) {
     G4PrimaryParticle* thePrim = thePrims.at(i);
     if (!thePrim) {
        edm::LogInfo("CastorShowerLibraryMaker") << "NULL Pointer to the primary" << std::endl;
        continue;
     }
// Check primary particle type
     int particleType = thePrim->GetPDGcode();

// set the pointer to the shower collection
     std::string SLType("");
     if (particleType==11) {
        SLShowerptr = &emSLHolder;
        SLType = "Electromagnetic";
     }
     else {
        SLShowerptr = &hadSLHolder;
        SLType = "Hadronic";
     }
     edm::LogInfo("CastorShowerLibraryMaker") << "\n Primary (thePrim) trackID is " << thePrim->GetTrackID() << "\n" ;

// Obtain primary particle's initial momentum (pInit)
     double px=0., py=0., pz=0., pInit = 0., eta = 0., phi = 0.;
     GetKinematics(particleType,px,py,pz,pInit,eta,phi);
// Check if current event falls into any bin
// first: energy
     if (pInit==0) {
        edm::LogInfo("CastorShowerLibraryMaker") << "Primary did not hit CASTOR" << std::endl;
        continue;
     }
     int ebin = FindEnergyBin(pInit);
     int etabin= FindEtaBin(eta);
     int phibin = FindPhiBin(phi);
     std::cout << SLType << std::endl;
     printSLstatus(ebin,etabin,phibin);
     if (!SLacceptEvent(ebin,etabin,phibin)) {
        edm::LogInfo("CastorShowerLibraryMaker") << "Event not accepted for ebin="
             << ebin<<",etabin="<<etabin<<",phibin="<<phibin
             <<"("<< pInit <<","<<eta<<","<<phi<<")" <<std::endl;
        continue;
     }
//
// event passed. Fill the vector accordingly
//  
// Look for the Hit Collection 
     edm::LogInfo("CastorShowerLibraryMaker")
        << "\n CastorShowerLibraryMaker::update(EndOfEvent * evt) - event #" 
        << (*evt)()->GetEventID() ;

/*
     std::cout << "Number of collections : " << allHC->GetNumberOfCollections() << std::endl;
     for(int ii = 0;ii<allHC->GetNumberOfCollections();ii++) 
        std::cout << "Name of collection " << ii << " : " << allHC->GetHC(ii)->GetName() << std::endl;
*/

     CastorShowerEvent* shower=NULL;
     int cur_evt_idx = SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin);
     shower = &(SLShowerptr->SLCollection.at(ebin).at(etabin).at(phibin).at(cur_evt_idx));

// Get Hit information
     if (FillShowerEvent(theCAFI,shower,particleType)) { 
//  Primary particle information
/*
        edm::LogInfo("CastorShowerLibraryMaker") << "New SL event: Primary = " << particleType
             << "; Energy = " << pInit << "; Eta = " << eta << "; Phi = " << phi
             << "; Nhits = " << shower->getNhit() << std::endl;
*/
        shower->setPrimE(pInit);
        shower->setPrimEta(eta);
        shower->setPrimPhi(phi);
        shower->setPrimX(PrimaryPosition[particleType].x());
        shower->setPrimY(PrimaryPosition[particleType].y());
        shower->setPrimZ(PrimaryPosition[particleType].z());
        SLnEvtInBinE(ebin)++;
        SLnEvtInBinEta(ebin,etabin)++;
        SLnEvtInBinPhi(ebin,etabin,phibin)++;
        NHitInEvent+=shower->getNhit();
        NEvtAccepted++;
      }
      else {shower->Clear();}
  }
// Check for unassociated energy
     
  if (NEvtAccepted==int(thePrims.size())&&theCAFI->entries()!=NHitInEvent){
     std::cout << "WARNING: Inconsistent Number of Hits -> Hits in collection: "<< theCAFI->entries()
          << "   Hits in the showers: " << NHitInEvent << std::endl;
     double miss_energy=0;
     double tot_energy=0;
     GetMissingEnergy(theCAFI,miss_energy,tot_energy);
     if (miss_energy>0) {
        std::cout << "Total missing energy: " << miss_energy
               << " for an incident energy: " << tot_energy
               << std::endl;
     }
  }
         
/*
  for (int i=emSLHolder.SLEnergyBins.size()-1;i>0;i--) {
      if (emSLHolder.nEvtInBinE.at(i)==(int)emSLHolder.nEvtPerBinE) {
         std::ostringstream out;
         out << emSLHolder.SLEnergyBins.at(i);
         cout << "Bin Limit: " << out.str() << std::endl;
         setenv("CASTOR_SL_PG_MAXE",out.str().c_str(),1);
       }
       break;
   }
*/
  //int iEvt = (*evt)()->GetEventID();
  //double xint;
/*
  if (modf(log10(iEvt),&xint)==0) 
    std::cout << " CastorShowerLibraryMaker Event " << iEvt << std::endl;
*/
  // std::cout << std::endl << "===>>> Done writing user histograms " << std::endl;
}

//========================= End of RUN ======================
void CastorShowerLibraryMaker::update(const EndOfRun * run)
{
// Fill the tree with the collected objects
  if (!IsSLReady()) SimG4Exception("\n\nShower Library  NOT READY.\n\n");

  unsigned int  ibine,ibineta,ibinphi,ievt; // indexes for em shower
  unsigned int  jbine,jbineta,jbinphi,jevt;// indexes for had shower

  ibine=ibineta=ibinphi=ievt=jbine=jbineta=jbinphi=jevt=0;

  int  nEvtInTree = 0;
  int  nEMevt  = emSLHolder.nEvtPerBinE*emSLHolder.SLEnergyBins.size();
  int  nHadevt = hadSLHolder.nEvtPerBinE*hadSLHolder.SLEnergyBins.size();
  int  maxEvtInTree=std::max(nEMevt,nHadevt);

  emInfo = &emSLHolder.SLInfo;
  hadInfo= &hadSLHolder.SLInfo;

  while(nEvtInTree<maxEvtInTree) {
    if (emShower) emShower->Clear();
    if (hadShower) hadShower->Clear();
    while(ibine<emSLHolder.SLEnergyBins.size()&&nEMevt>0){
      emShower = &(emSLHolder.SLCollection.at(ibine).at(ibineta).at(ibinphi).at(ievt));
      ievt++;
      if (ievt==emSLHolder.nEvtPerBinPhi) {ievt=0;ibinphi++;}
      if (ibinphi==emSLHolder.SLPhiBins.size()) {ibinphi=0;ibineta++;}
      if (ibineta==emSLHolder.SLEtaBins.size()) {ibineta=0;ibine++;}
      break;
    }
    while(jbine<hadSLHolder.SLEnergyBins.size()&&nHadevt>0){
      hadShower = &(hadSLHolder.SLCollection.at(jbine).at(jbineta).at(jbinphi).at(jevt));
      jevt++;
      if (jevt==hadSLHolder.nEvtPerBinPhi) {jevt=0;jbinphi++;}
      if (jbinphi==hadSLHolder.SLPhiBins.size()) {jbinphi=0;jbineta++;}
      if (jbineta==hadSLHolder.SLEtaBins.size()) {jbineta=0;jbine++;}
      break;
    }
    theTree->Fill();
    nEvtInTree++;
    if (nEvtInTree==1) {
       theTree->SetBranchStatus("emShowerLibInfo.",0);
       theTree->SetBranchStatus("hadShowerLibInfo.",0);
    }
  }
// check if run is NULL and exit
  if (run==NULL) throw SimG4Exception("\n\nNumber of needed trigger events reached in CastorShowerLibraryMaker\n\n");
}

//============================================================
void CastorShowerLibraryMaker::Finish() {

  // if (doNTcastorevent) {

  theFile->cd();
  theTree->Write("",TObject::kOverwrite);
  std::cout << "CastorShowerLibraryMaker: Ntuple event written" << std::endl;   
  theFile->Close();
  std::cout << "CastorShowerLibraryMaker: Event file closed" << std::endl;

  // Delete pointers to objects, now that TTree has been written and TFile closed
//  delete      info; 
//  delete  emShower;
//  delete hadShower;
  // }
} 
int CastorShowerLibraryMaker::FindEnergyBin(double energy) {
  //
  // returns the integer index of the energy bin, taken from SLenergies vector
  // returns -1 if ouside valid range
  //
  if (!SLShowerptr) {
     edm::LogInfo("CastorShowerLibraryMaker") << "\n\nFindEnergyBin can be called only after BeginOfEvent\n\n";
     throw SimG4Exception("\n\nNULL Pointer to the shower library.\n\n");
  }
  const std::vector<double>& SLenergies = SLShowerptr->SLEnergyBins;
  if (energy >= SLenergies.back()) return SLenergies.size()-1;

  unsigned int i = 0;
  for(;i<SLenergies.size()-1;i++) 
    if (energy >= SLenergies.at(i) && energy < SLenergies.at(i+1)) return (int)i;

  // now i points to the last but 1 bin
  if (energy>=SLenergies.at(i)) return (int)i;
  // energy outside bin range
  return -1;
}
int CastorShowerLibraryMaker::FindEtaBin(double eta) {
  //
  // returns the integer index of the eta bin, taken from SLetas vector
  // returns -1 if ouside valid range
  //
  if (!SLShowerptr) {
     edm::LogInfo("CastorShowerLibraryMaker") << "\n\nFindEtaBin can be called only after BeginOfEvent\n\n";
     throw SimG4Exception("\n\nNULL Pointer to the shower library.\n\n");
  }
  const std::vector<double>& SLetas = SLShowerptr->SLEtaBins;
  if (eta>=SLetas.back()) return SLetas.size()-1;
  unsigned int i = 0;
  for(;i<SLetas.size()-1;i++)
     if (eta >= SLetas.at(i) && eta < SLetas.at(i+1)) return (int)i;
  // now i points to the last but 1 bin
  if (eta>=SLetas.at(i)) return (int)i;
  // eta outside bin range
  return -1;
}
int CastorShowerLibraryMaker::FindPhiBin(double phi) {
  //
  // returns the integer index of the phi bin, taken from SLphis vector
  // returns -1 if ouside valid range
  //
  // needs protection in case phi is outside range -pi,pi
  //
  if (!SLShowerptr) {
     edm::LogInfo("CastorShowerLibraryMaker") << "\n\nFindPhiBin can be called only after BeginOfEvent\n\n";
     throw SimG4Exception("\n\nNULL Pointer to the shower library.\n\n");
  }
  const std::vector<double>& SLphis = SLShowerptr->SLPhiBins;
  if (phi>=SLphis.back()) return SLphis.size()-1;
  unsigned int i = 0;
  for(;i<SLphis.size()-1;i++)
     if (phi >= SLphis.at(i) && phi < SLphis.at(i+1)) return (int)i;
  // now i points to the last but 1 bin
  if (phi>=SLphis.at(i)) return (int)i;
  // phi outside bin range
  return -1;
}
bool CastorShowerLibraryMaker::IsSLReady()
{
// at this point, the pointer to the shower library should be NULL
  if (SLShowerptr) {
     edm::LogInfo("CastorShowerLibraryMaker") << "\n\nIsSLReady must be called when a new event starts.\n\n";
     throw SimG4Exception("\n\nNOT NULL Pointer to the shower library.\n\n");
  }
// it is enough to check if all the energy bin is filled
  if (DoEmSL) {
     SLShowerptr = &emSLHolder;
     for(unsigned int i=0;i<SLShowerptr->SLEnergyBins.size();i++) {
        if (!SLisEBinFilled(i)) {
           SLShowerptr=NULL;
           return false;
        }
     }
  }
  if (DoHadSL) {
     SLShowerptr = &hadSLHolder;
     for(unsigned int i=0;i<SLShowerptr->SLEnergyBins.size();i++) {
        if (!SLisEBinFilled(i)) {
           SLShowerptr=NULL;
           return false;
        }
     }
  }
  SLShowerptr=NULL;
  return true;
}
void CastorShowerLibraryMaker::GetKinematics(int thePrim,double& px, double& py, double& pz, double& pInit, double& eta, double& phi)
{
    if (thePrim==0) return;
    if (PrimaryMomentum.find(thePrim)==PrimaryMomentum.end()) return;
    px = PrimaryMomentum[thePrim].x()/GeV;
    py = PrimaryMomentum[thePrim].y()/GeV;
    pz = PrimaryMomentum[thePrim].z()/GeV;
    pInit=PrimaryMomentum[thePrim].mag()/GeV;
    if (pInit==0) return;
    double costheta = pz/pInit;
    double theta = acos(std::min(std::max(costheta,double(-1.)),double(1.)));
    eta = -log(tan(theta/2.0));
    phi = (px==0 && py==0) ? 0 : atan2(py,px); // the recommended way of calculating phi
    phi = PrimaryMomentum[thePrim].phi();
}
void CastorShowerLibraryMaker::GetKinematics(G4PrimaryParticle* thePrim,double& px, double& py, double& pz, double& pInit, double& eta, double& phi)
{
    px=py=pz=phi=eta=0.0;
    if (thePrim==0) return;
    px = thePrim->GetMomentum().x()/GeV;
    py = thePrim->GetMomentum().y()/GeV;
    pz = thePrim->GetMomentum().z()/GeV;
    pInit = thePrim->GetMomentum().mag()/GeV;
    //pInit = sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
    if (pInit==0) return;
    double costheta = pz/pInit;
    double theta = acos(std::min(std::max(costheta,double(-1.)),double(1.)));
    eta = -log(tan(theta/2.0));
    phi = (px==0 && py==0) ? 0 : atan2(py,px); // the recommended way of calculating phi
    phi = thePrim->GetMomentum().phi();
    //if (px!=0) phi=atan(py/px);
}
std::vector<G4PrimaryParticle*> CastorShowerLibraryMaker::GetPrimary(const G4Event * evt)
{
  // Find Primary info:
  int trackID = 0;
  std::vector<G4PrimaryParticle*> thePrims;
  G4PrimaryParticle* thePrim = 0;
  G4int nvertex = evt->GetNumberOfPrimaryVertex();
  edm::LogInfo("CastorShowerLibraryMaker")  << "Event has " << nvertex << " vertex";   
  if (nvertex!=1) {
     edm::LogInfo("CastorShowerLibraryMaker") << "CastorShowerLibraryMaker::GetPrimary ERROR: no vertex";
     return thePrims;
  }

  for (int i = 0 ; i<nvertex; i++) {
    G4PrimaryVertex* avertex = evt->GetPrimaryVertex(i);
    if (avertex == 0) {
       edm::LogInfo("CastorShowerLibraryMaker")
          << "CastorShowerLibraryMaker::GetPrimary ERROR: pointer to vertex = 0";
       continue;
    }
    unsigned int npart = avertex->GetNumberOfParticle();
    if (npart!=NPGParticle) continue;
    for (unsigned int j=0;j<npart;j++) {
        unsigned int k = 0;
        //int test_pID = 0;
        trackID = j;
        thePrim=avertex->GetPrimary(trackID);
        while(k<NPGParticle&&PGParticleIDs.at(k++)!=thePrim->GetPDGcode()){;};
        if (k>NPGParticle) continue; // ID not found in the requested particles
        thePrims.push_back(thePrim);
    }
  }
  return thePrims;
}
void CastorShowerLibraryMaker::printSLstatus(int ebin,int etabin,int phibin)
{
  if (!SLShowerptr) {
     edm::LogInfo("CastorShowerLibraryInfo") << "NULL shower pointer. Printing both";
     std::cout << "Electromagnetic" << std::endl;
     SLShowerptr = &emSLHolder;
     this->printSLstatus(ebin,etabin,phibin);
     std::cout << "Hadronic" << std::endl;
     SLShowerptr = &hadSLHolder;
     this->printSLstatus(ebin,etabin,phibin);
     SLShowerptr = NULL;
     return;
  }
  int nBinsE  =SLShowerptr->SLEnergyBins.size();
  int nBinsEta=SLShowerptr->SLEtaBins.size();
  int nBinsPhi=SLShowerptr->SLPhiBins.size();
  std::vector<double> SLenergies = SLShowerptr->SLEnergyBins;
  for(int n=0;n<11+(nBinsEta*nBinsPhi);n++) std::cout << "=";
  std::cout << std::endl;
  for(int i=0;i<nBinsE;i++) {
     std::cout << "E bin " << std::setw(6) << SLenergies.at(i) << " : ";
     for(int j=0;j<nBinsEta;j++) {
        for(int k=0;k<nBinsPhi;k++) {
           (SLisPhiBinFilled(i,j,k))?std::cout << "1":std::cout << "-";
        }
        if (j<nBinsEta-1) std::cout << "|";
     }
     std::cout << " (" << SLnEvtInBinE(i) << " events)";
     std::cout << std::endl;
     if (ebin!=i) continue;
     std::cout << "               ";
     for(int j=0;j<nBinsEta;j++) {
        for(int k=0;k<nBinsPhi;k++) {
           (ebin==i&&etabin==j&&phibin==k)?std::cout <<  "^":std::cout << " ";
        }
        if (j<nBinsEta-1) std::cout << " ";
     }
     std::cout << std::endl;
  }
  for(int n=0;n<11+(nBinsEta*nBinsPhi);n++) std::cout << "=";
  std::cout << std::endl;
}
bool CastorShowerLibraryMaker::SLacceptEvent(int ebin, int etabin, int phibin)
{
     if (SLShowerptr==NULL) {
        edm::LogInfo("CastorShowerLibraryMaker::SLacceptEvent:") << "Error. NULL pointer to CastorShowerEvent";
        return false;
     }
     if (ebin<0||ebin>=int(SLShowerptr->SLEnergyBins.size())) return false;
     if (SLisEBinFilled(ebin)) return false;

     if (etabin<0||etabin>=int(SLShowerptr->SLEtaBins.size())) return false;
     if (SLisEtaBinFilled(ebin,etabin)) return false;

     if (phibin<0||phibin>=int(SLShowerptr->SLPhiBins.size())) return false;
     if (SLisPhiBinFilled(ebin,etabin,phibin)) return false;
     return true;
}
bool CastorShowerLibraryMaker::FillShowerEvent(CaloG4HitCollection* theCAFI, CastorShowerEvent* shower,int ipart)
{
     unsigned int volumeID=0;
     double en_in_fi = 0.;
     //double totalEnergy = 0;

     int nentries = theCAFI->entries();

// Compute Total Energy in CastorFI volume
/*
     for(int ihit = 0; ihit < nentries; ihit++) {
       CaloG4Hit* aHit = (*theCAFI)[ihit];
       totalEnergy += aHit->getEnergyDeposit();
     }
*/
     if (!shower) {
        edm::LogInfo("CastorShowerLibraryMaker") << "Error. NULL pointer to CastorShowerEvent";
        return false;
     }

     CastorNumberingScheme *theCastorNumScheme = new CastorNumberingScheme();
// Hit position
     math::XYZPoint entry;
     math::XYZPoint position;
     int nHits;
     nHits=0;
     for (int ihit = 0; ihit < nentries; ihit++) {
       CaloG4Hit* aHit  = (*theCAFI)[ihit];
       int hit_particleID = aHit->getTrackID();
       if (MapOfSecondaries[ipart].find(hit_particleID)==MapOfSecondaries[ipart].end()) {
          if (verbosity) edm::LogInfo("CastorShowerLibraryMaker")
              << "Skipping hit from trackID " << hit_particleID;
          continue;
       }
       volumeID         = aHit->getUnitID();
       double hitEnergy = aHit->getEnergyDeposit();
       en_in_fi        += aHit->getEnergyDeposit();
       float time       = aHit->getTimeSlice();
       int zside, sector, zmodule;
       theCastorNumScheme->unpackIndex(volumeID, zside, sector,zmodule);
       entry    = aHit->getEntry();
       position = aHit->getPosition();
       if (verbosity) edm::LogInfo("CastorShowerLibraryMaker")
          << "\n side , sector , module = " << zside << " , " 
          << sector << " , " << zmodule 
          << "\n nphotons = " << hitEnergy ; 

       if (verbosity) edm::LogInfo("CastorShowerLibraryMaker")
          << "\n packIndex = " 
          << theCastorNumScheme->packIndex(zside, sector,zmodule);

       if(verbosity&&time>100.) {
         edm::LogInfo("CastorShowerLibraryMaker")
            << "\n nentries = " << nentries 
            << "\n     time[" << ihit << "] = " << time
            << "\n  trackID[" << ihit << "] = " << aHit->getTrackID()
            << "\n volumeID[" << ihit << "] = " << volumeID 
            << "\n nphotons[" << ihit << "] = " << hitEnergy 
            << "\n side, sector, module  = " << zside <<", " << sector<<", " << zmodule
            << "\n packIndex " << theCastorNumScheme->packIndex(zside,sector,zmodule)
            << "\n X,Y,Z = " << entry.x() << ","<< entry.y() << "," << entry.z();
       }
       if (verbosity) edm::LogInfo("CastorShowerLibraryMaker") << "\n    Incident Energy = "  
                                                << aHit->getIncidentEnergy() << " \n" ;

//  CaloG4Hit information 
       shower->setDetID(volumeID);
       shower->setHitPosition(position);
       shower->setNphotons(hitEnergy);
       shower->setTime(time);
       nHits++;
     }
// Write number of hits to CastorShowerEvent instance
     if (nHits==0) {
        edm::LogInfo("CastorShowerLibraryMaker") << "No hits found for this track (trackID=" << ipart << ")." << std::endl;
        if (theCastorNumScheme) delete theCastorNumScheme;
        return false;
     }
     shower->setNhit(nHits);

// update the event counters
     if (theCastorNumScheme) delete theCastorNumScheme;
     return true;
}
int& CastorShowerLibraryMaker::SLnEvtInBinE(int ebin)
{
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinE can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.");
   }
   return SLShowerptr->nEvtInBinE.at(ebin);
}

int& CastorShowerLibraryMaker::SLnEvtInBinEta(int ebin, int etabin)
{
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinEta can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.");
   }
   return SLShowerptr->nEvtInBinEta.at(ebin).at(etabin);
}

int& CastorShowerLibraryMaker::SLnEvtInBinPhi(int ebin, int etabin, int phibin)
{
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinPhi can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.");
   }
   return SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin);
}
bool CastorShowerLibraryMaker::SLisEBinFilled(int ebin)
{
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLisEBinFilled can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.");
   }
   if (SLShowerptr->nEvtInBinE.at(ebin)<(int)SLShowerptr->nEvtPerBinE) return false;
   return true;
}
bool CastorShowerLibraryMaker::SLisEtaBinFilled(int ebin,int etabin)
{
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLisEtaBinFilled can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.");
   }
   if (SLShowerptr->nEvtInBinEta.at(ebin).at(etabin)<(int)SLShowerptr->nEvtPerBinEta) return false;
   return true;
}
bool CastorShowerLibraryMaker::SLisPhiBinFilled(int ebin,int etabin,int phibin)
{
   if (!SLShowerptr) {
      edm::LogInfo("CastorShowerLibraryMaker") << "\n\nSLisPhiBinFilled can be called only after BeginOfEvent\n\n";
      throw SimG4Exception("\n\nNULL Pointer to the shower library.");
   }
   if (SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin)<(int)SLShowerptr->nEvtPerBinPhi) return false;
   return true;
}
void CastorShowerLibraryMaker::KillSecondaries(const G4Step * aStep) {
   G4TrackVector *p_sec = const_cast<G4TrackVector*>(aStep->GetSecondary());
   for(int i=0;i<int(p_sec->size());i++) {
      /*if (verbosity)*/ std::cout << "Killing track ID " << p_sec->at(i)->GetTrackID() << " and its secondaries"
          << " Produced by Process " << p_sec->at(i)->GetCreatorProcess()->GetProcessName()
          << " in the volume " << aStep->GetTrack()->GetVolume()->GetName()
          << std::endl;
      p_sec->at(i)->SetTrackStatus(fKillTrackAndSecondaries);
   }
}

void CastorShowerLibraryMaker::GetMissingEnergy(CaloG4HitCollection* theCAFI,double& miss_energy,double& tot_energy){
// Get the total deposited energy and the one from hit not associated to a primary
     miss_energy = 0;
     tot_energy = 0;
     for (int ihit = 0; ihit < theCAFI->entries(); ihit++) {
         int id = (*theCAFI)[ihit]->getTrackID();
         tot_energy+=(*theCAFI)[ihit]->getEnergyDeposit();
         int hit_prim = 0;
         for(unsigned int i=0;i<thePrims.size();i++) {
            int particleType = thePrims.at(i)->GetPDGcode();
            if (MapOfSecondaries[particleType].find(id)!=MapOfSecondaries[particleType].end())
               hit_prim = particleType;
         }
         if (hit_prim==0) {
            std::cout << "Track ID " << id << " produced a hit which is not associated with a primary."
            << std::endl;
            miss_energy+=(*theCAFI)[ihit]->getEnergyDeposit();
         }
     }
}