CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

HcalTB06Analysis Class Reference

#include <SimG4CMS/HcalTestBeam/interface/HcalTB06Analysis.h>

Inheritance diagram for HcalTB06Analysis:
SimProducer Observer< const BeginOfRun * > Observer< const BeginOfEvent * > Observer< const EndOfEvent * > Observer< const G4Step * > SimWatcher

List of all members.

Public Member Functions

 HcalTB06Analysis (const edm::ParameterSet &p)
virtual void produce (edm::Event &, const edm::EventSetup &)
virtual ~HcalTB06Analysis ()

Private Member Functions

void clear ()
void fillBuffer (const EndOfEvent *evt)
void fillEvent (PHcalTB06Info &)
void finalAnalysis ()
 HcalTB06Analysis (const HcalTB06Analysis &)
void init ()
const HcalTB06Analysisoperator= (const HcalTB06Analysis &)
void update (const BeginOfEvent *evt)
 This routine will be called when the appropriate signal arrives.
void update (const BeginOfRun *run)
 This routine will be called when the appropriate signal arrives.
void update (const EndOfEvent *evt)
 This routine will be called when the appropriate signal arrives.
void update (const G4Step *step)
 This routine will be called when the appropriate signal arrives.

Private Attributes

G4RotationMatrix * beamline_RM
double beamOffset
int count
std::vector< CaloHitecalHitCache
double eecals
double ehcals
double etaInit
double etots
int evNum
std::vector< CaloHithcalHitCache
std::vector< CaloHithcalHitLayer
HcalTB06Histohisto
int iceta
int icphi
std::vector< std::string > names
int nPrimary
int particleType
double phiInit
double pInit
bool pvFound
G4ThreeVector pvMomentum
G4ThreeVector pvPosition
int pvType
G4ThreeVector pvUVW
std::vector< double > secEkin
std::vector< G4ThreeVector > secMomentum
std::vector< int > secPartID
std::vector< int > secTrackID
std::vector< int > shortLivedSecondaries

Detailed Description

Description: Analysis of 2004 Hcal Test beam simulation

Usage: A Simwatcher class and can be activated from Oscarproducer module

Definition at line 44 of file HcalTB06Analysis.h.


Constructor & Destructor Documentation

HcalTB06Analysis::HcalTB06Analysis ( const edm::ParameterSet p)

Definition at line 47 of file HcalTB06Analysis.cc.

References beamline_RM, beamOffset, funct::exp(), edm::ParameterSet::getParameter(), histo, iceta, icphi, init(), and names.

                                                          : histo(0) {

  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB06Analysis");
  names          = m_Anal.getParameter<std::vector<std::string> >("Names");
  beamOffset     =-m_Anal.getParameter<double>("BeamPosition")*cm;
  double fMinEta = m_Anal.getParameter<double>("MinEta");
  double fMaxEta = m_Anal.getParameter<double>("MaxEta");
  double fMinPhi = m_Anal.getParameter<double>("MinPhi");
  double fMaxPhi = m_Anal.getParameter<double>("MaxPhi");
  double beamEta = (fMaxEta+fMinEta)/2.;
  double beamPhi = (fMaxPhi+fMinPhi)/2.;
  double beamThet= 2*atan(exp(-beamEta));
  if (beamPhi < 0) beamPhi += twopi;
  iceta          = (int)(beamEta/0.087) + 1;
  icphi          = (int)(fabs(beamPhi)/0.087) + 5;
  if (icphi > 72) icphi -= 73;

  produces<PHcalTB06Info>();

  beamline_RM = new G4RotationMatrix;
  beamline_RM->rotateZ(-beamPhi);
  beamline_RM->rotateY(-beamThet);
 
  edm::LogInfo("HcalTBSim") << "HcalTB06:: Initialised as observer of BeginOf"
                            << "Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
                            << " with Parameter values:\n \tbeamOffset = " 
                            << beamOffset << "\ticeta = " << iceta 
                            << "\ticphi = " << icphi << "\n\tbeamline_RM = "
                            << *beamline_RM;

  init();

  histo  = new HcalTB06Histo(m_Anal);
} 
HcalTB06Analysis::~HcalTB06Analysis ( ) [virtual]

Definition at line 82 of file HcalTB06Analysis.cc.

References count, and histo.

                                    {

  edm::LogInfo("HcalTBSim") << "\n -------->  Total number of selected entries"
                            << " : " << count << "\nPointers:: Histo " <<histo;
  if (histo)   {
    delete histo;
    histo  = 0;
  }
}
HcalTB06Analysis::HcalTB06Analysis ( const HcalTB06Analysis ) [private]

Member Function Documentation

void HcalTB06Analysis::clear ( void  ) [private]

Definition at line 519 of file HcalTB06Analysis.cc.

References ecalHitCache, etaInit, hcalHitCache, nPrimary, particleType, phiInit, pInit, pvFound, pvMomentum, pvPosition, pvType, pvUVW, secEkin, secMomentum, secPartID, secTrackID, and shortLivedSecondaries.

Referenced by init(), and update().

                             {

  pvFound = false;
  pvType  =-2;
  pvPosition = G4ThreeVector();
  pvMomentum = G4ThreeVector();
  pvUVW      = G4ThreeVector();
  secTrackID.clear();
  secPartID.clear();
  secMomentum.clear();
  secEkin.clear();
  shortLivedSecondaries.clear();

  ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end()); 
  hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end()); 
  nPrimary = particleType = 0;
  pInit = etaInit = phiInit = 0;
}
void HcalTB06Analysis::fillBuffer ( const EndOfEvent evt) [private]

Definition at line 264 of file HcalTB06Analysis.cc.

References ecalHitCache, eta(), etaInit, evNum, CaloG4Hit::getEnergyDeposit(), CaloG4Hit::getEntry(), CaloG4Hit::getTimeSlice(), CaloG4Hit::getTrackID(), CaloG4Hit::getUnitID(), hcalHitCache, i, j, funct::log(), LogDebug, max(), min, names, npart, nPrimary, AlCaHLTBitMon_ParallelJobs::p, particleType, phi, phiInit, pInit, pos, funct::pow(), python::multivaluedict::sort(), mathSSE::sqrt(), funct::tan(), theta(), and cond::rpcobgas::time.

Referenced by update().

                                                        {

  std::vector<CaloHit> hhits;
  int                  idHC, j;
  CaloG4HitCollection* theHC;
  std::map<int,float,std::less<int> > primaries;
  double               etot1=0, etot2=0;

  // Look for the Hit Collection of HCal
  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
  std::string sdName = names[0];
  idHC  = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
  theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hit Collection for " << sdName
                        << " of ID " << idHC << " is obtained at " << theHC;

  if (idHC >= 0 && theHC > 0) {
    hhits.reserve(theHC->entries());
    for (j = 0; j < theHC->entries(); j++) {
      CaloG4Hit* aHit = (*theHC)[j]; 
      double e        = aHit->getEnergyDeposit()/GeV;
      double time     = aHit->getTimeSlice();
      math::XYZPoint pos  = aHit->getEntry();
      unsigned int id = aHit->getUnitID();
      double theta    = pos.theta();
      double eta      = -log(tan(theta * 0.5));
      double phi      = pos.phi();
      CaloHit hit(2,1,e,eta,phi,time,id);
      hhits.push_back(hit);
      primaries[aHit->getTrackID()]+= e;
      etot1 += e;
#ifdef ddebug
      LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hcal Hit i/p " << j 
                            << "  ID 0x" << std::hex << id << std::dec 
                            << " time " << std::setw(6) << time << " theta "
                            << std::setw(8) << theta << " eta " << std::setw(8)
                            << eta << " phi " << std::setw(8) << phi << " e " 
                            << std::setw(8) << e;
#endif
    }
  }

  // Add hits in the same channel within same time slice
  std::vector<CaloHit>::iterator itr;
  int nHit = hhits.size();
  std::vector<CaloHit*> hits(nHit);
  for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
    hits[j] = &hhits[j];
  }
  sort(hits.begin(),hits.end(),CaloHitIdMore());
  std::vector<CaloHit*>::iterator k1, k2;
  int nhit = 0;
  for (k1 = hits.begin(); k1 != hits.end(); k1++) {
    int      det    = (**k1).det();
    int      layer  = (**k1).layer();
    double   ehit   = (**k1).e();
    double   eta    = (**k1).eta();
    double   phi    = (**k1).phi();
    double   jitter = (**k1).t();
    uint32_t unitID = (**k1).id();
    int      jump  = 0;
    for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
           unitID==(**k2).id(); k2++) {
      ehit += (**k2).e();
      jump++;
    }
    nhit++;
    CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
    hcalHitCache.push_back(hit);
    etot2 += ehit;
    k1    += jump;
#ifdef ddebug
    LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hcal Hit store " << nhit 
                          << "  ID 0x" << std::hex  << unitID  << std::dec 
                          << " time " << std::setw(6) << jitter << " eta "
                          << std::setw(8) << eta << " phi " << std::setw(8) 
                          << phi  << " e " << std::setw(8) << ehit;
#endif
  }
  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Stores " << nhit << " HCal hits"
                        << " from " << nHit << " input hits E(Hcal) " << etot1 
                        << " " << etot2;
  
  // Look for the Hit Collection of ECal
  std::vector<CaloHit> ehits;
  sdName= names[1];
  idHC  = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
  theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
  etot1 = etot2 = 0;
  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Hit Collection for " << sdName
                        << " of ID " << idHC << " is obtained at " << theHC;
  if (idHC >= 0 && theHC > 0) {
    ehits.reserve(theHC->entries());
    for (j = 0; j < theHC->entries(); j++) {
      CaloG4Hit* aHit = (*theHC)[j]; 
      double e        = aHit->getEnergyDeposit()/GeV;
      double time     = aHit->getTimeSlice();
      math::XYZPoint pos  = aHit->getEntry();
      unsigned int id = aHit->getUnitID();
      double theta    = pos.theta();
      double eta      = -log(tan(theta * 0.5));
      double phi      = pos.phi();
      if (e < 0 || e > 100000.) e = 0;
      CaloHit hit(1,0,e,eta,phi,time,id);
      ehits.push_back(hit);
      primaries[aHit->getTrackID()]+= e;
      etot1 += e;
#ifdef ddebug
      LogDebug("HcalTBSim") << "HcalTB06Analysis:: Ecal Hit i/p " << j 
                            << "  ID 0x" << std::hex << id  << std::dec 
                            << " time " << std::setw(6) << time << " theta " 
                            << std::setw(8) << theta  << " eta " <<std::setw(8)
                            << eta  << " phi " << std::setw(8) << phi << " e "
                            << std::setw(8) << e;
#endif
    }
  }

  // Add hits in the same channel within same time slice
  nHit = ehits.size();
  std::vector<CaloHit*> hite(nHit);
  for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
    hite[j] = &ehits[j];
  }
  sort(hite.begin(),hite.end(),CaloHitIdMore());
  nhit = 0;
  for (k1 = hite.begin(); k1 != hite.end(); k1++) {
    int      det    = (**k1).det();
    int      layer  = (**k1).layer();
    double   ehit   = (**k1).e();
    double   eta    = (**k1).eta();
    double   phi    = (**k1).phi();
    double   jitter = (**k1).t();
    uint32_t unitID = (**k1).id();
    int      jump  = 0;
    for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
           unitID==(**k2).id(); k2++) {
      ehit += (**k2).e();
      jump++;
    }
    nhit++;
    CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
    ecalHitCache.push_back(hit);
    etot2 += ehit;
    k1    += jump;
#ifdef ddebug
    LogDebug("HcalTBSim") << "HcalTB06Analysis:: Ecal Hit store " << nhit
                          << "  ID 0x" << std::hex << unitID  << std::dec 
                          << " time " << std::setw(6) << jitter << " eta "
                          << std::setw(8) << eta << " phi " << std::setw(8)
                          << phi << " e " << std::setw(8) << ehit;
#endif
  }
  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Stores " << nhit << " ECal hits"
                        << " from " << nHit << " input hits E(Ecal) " << etot1 
                        << " " << etot2;

  // Find Primary info:
  nPrimary    = (int)(primaries.size());
  int trackID = 0;
  G4PrimaryParticle* thePrim=0;
  int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Event has " << nvertex 
                        << " verteices";
  if (nvertex<=0)
    edm::LogInfo("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERROR: no "
                              << "vertex found for event " << evNum;

  for (int i = 0 ; i<nvertex; i++) {
    G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
    if (avertex == 0) {
      edm::LogInfo("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERR: pointer "
                                << "to vertex = 0 for event " << evNum;
    } else {
      LogDebug("HcalTBSim") << "HcalTB06Analysis::Vertex number :" << i << " "
                            << avertex->GetPosition();
      int npart = avertex->GetNumberOfParticle();
      if (npart == 0)
        edm::LogWarning("HcalTBSim") << "HcalTB06Analysis::End Of Event ERR: "
                                     << "no primary!";
      if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
    }
  }
    
  if (thePrim != 0) {
    double px = thePrim->GetPx();
    double py = thePrim->GetPy();
    double pz = thePrim->GetPz();
    double p  = std::sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
    pInit     = p/GeV;
    if (p==0) 
      edm::LogWarning("HcalTBSim") << "HcalTB06Analysis:: EndOfEvent ERR: "
                                   << "primary has p=0 ";
    else {
      double costheta = pz/p;
      double theta = acos(std::min(std::max(costheta,-1.),1.));
      etaInit = -log(tan(theta/2));
      if (px != 0 || py != 0) phiInit = atan2(py,px);  
    }
    particleType = thePrim->GetPDGcode();
  } else 
    edm::LogWarning("HcalTBSim") << "HcalTB06Analysis::EndOfEvent ERR: could "
                                 << "not find primary";

}
void HcalTB06Analysis::fillEvent ( PHcalTB06Info product) [private]

Definition at line 491 of file HcalTB06Analysis.cc.

References ecalHitCache, eecals, ehcals, etaInit, etots, evNum, hcalHitCache, i, nPrimary, particleType, phiInit, pInit, pvMomentum, pvPosition, pvType, pvUVW, PHcalTB06Info::saveHit(), secEkin, secMomentum, secPartID, secTrackID, PHcalTB06Info::setEdep(), PHcalTB06Info::setPrimary(), PHcalTB06Info::setVtxPrim(), PHcalTB06Info::setVtxSec(), x, detailsBasic3DVector::y, and z.

Referenced by produce().

                                                        {

  //Beam Information
  product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);

  // Total Energy
  product.setEdep(etots, eecals, ehcals);

  //Energy deposits in the crystals and towers
  for (unsigned int i=0; i<hcalHitCache.size(); i++) 
    product.saveHit(hcalHitCache[i].id(), hcalHitCache[i].eta(),
                    hcalHitCache[i].phi(), hcalHitCache[i].e(),
                    hcalHitCache[i].t());
  for (unsigned int i=0; i<ecalHitCache.size(); i++) 
    product.saveHit(ecalHitCache[i].id(), ecalHitCache[i].eta(),
                    ecalHitCache[i].phi(), ecalHitCache[i].e(),
                    ecalHitCache[i].t());

  //Vertex associated quantities
  product.setVtxPrim(evNum, pvType, pvPosition.x(), pvPosition.y(), 
                     pvPosition.z(), pvUVW.x(), pvUVW.y(), pvUVW.z(),
                     pvMomentum.x(), pvMomentum.y(), pvMomentum.z());
  for (unsigned int i = 0; i < secTrackID.size(); i++) {
    product.setVtxSec(secTrackID[i], secPartID[i], secMomentum[i].x(),
                      secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
  }
}
void HcalTB06Analysis::finalAnalysis ( ) [private]

Definition at line 470 of file HcalTB06Analysis.cc.

References ecalHitCache, eecals, ehcals, etaInit, etots, HcalTB06Histo::fillEdep(), HcalTB06Histo::fillPrimary(), hcalHitCache, histo, i, LogDebug, phiInit, and pInit.

Referenced by update().

                                     {

  //Beam Information
  histo->fillPrimary(pInit, etaInit, phiInit);

  // Total Energy
  eecals = ehcals = 0.;
  for (unsigned int i=0; i<hcalHitCache.size(); i++) {
    ehcals += hcalHitCache[i].e();
  }
  for (unsigned int i=0; i<ecalHitCache.size(); i++) {
    eecals += ecalHitCache[i].e();
  }
  etots = eecals + ehcals;
  LogDebug("HcalTBSim") << "HcalTB06Analysis:: Energy deposit at Sim Level "
                        << "(Total) " << etots << " (ECal) " << eecals 
                        << " (HCal) " << ehcals;
  histo->fillEdep(etots, eecals, ehcals);
}
void HcalTB06Analysis::init ( void  ) [private]

Definition at line 103 of file HcalTB06Analysis.cc.

References clear(), count, and evNum.

Referenced by HcalTB06Analysis().

                            {

  // counter 
  count = 0;
  evNum = 0;
  clear();
}
const HcalTB06Analysis& HcalTB06Analysis::operator= ( const HcalTB06Analysis ) [private]
void HcalTB06Analysis::produce ( edm::Event e,
const edm::EventSetup  
) [virtual]

Implements SimProducer.

Definition at line 96 of file HcalTB06Analysis.cc.

References fillEvent(), and edm::Event::put().

                                                                {

  std::auto_ptr<PHcalTB06Info> product(new PHcalTB06Info);
  fillEvent(*product);
  e.put(product);
}
void HcalTB06Analysis::update ( const BeginOfEvent ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 118 of file HcalTB06Analysis.cc.

References clear(), and evNum.

                                                      {
 
  evNum = (*evt) ()->GetEventID ();
  clear();
  edm::LogInfo("HcalTBSim") << "HcalTB06Analysis: =====> Begin of event = "
                            << evNum;
}
void HcalTB06Analysis::update ( const BeginOfRun ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfRun * >.

Definition at line 111 of file HcalTB06Analysis.cc.

                                                    {

  int irun = (*run)()->GetRunID();
  edm::LogInfo("HcalTBSim") <<" =====> Begin of Run = " << irun;
 
}
void HcalTB06Analysis::update ( const EndOfEvent ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 240 of file HcalTB06Analysis.cc.

References count, fillBuffer(), finalAnalysis(), and LogDebug.

                                                    {

  count++;

  //fill the buffer
  LogDebug("HcalTBSim") << "HcalTB06Analysis::Fill event " 
                        << (*evt)()->GetEventID();
  fillBuffer (evt);
  
  //Final Analysis
  LogDebug("HcalTBSim") << "HcalTB06Analysis::Final analysis";  
  finalAnalysis();

  int iEvt = (*evt)()->GetEventID();
  if (iEvt < 10) 
    edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
  else if ((iEvt < 100) && (iEvt%10 == 0)) 
    edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
  else if ((iEvt < 1000) && (iEvt%100 == 0)) 
    edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
  else if ((iEvt < 10000) && (iEvt%1000 == 0)) 
    edm::LogInfo("HcalTBSim") << "HcalTB06Analysis:: Event " << iEvt;
}
void HcalTB06Analysis::update ( const G4Step *  ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 126 of file HcalTB06Analysis.cc.

References LogDebug, NULL, evf::utils::pid, pos, position, pvFound, pvMomentum, pvPosition, pvType, pvUVW, secEkin, secMomentum, secPartID, secTrackID, and shortLivedSecondaries.

                                                  {

  if (aStep != NULL) {
    //Get Step properties
    G4ThreeVector thePreStepPoint  = aStep->GetPreStepPoint()->GetPosition();
    G4ThreeVector thePostStepPoint;

    // Get Tracks properties
    G4Track*      aTrack   = aStep->GetTrack();
    int           trackID  = aTrack->GetTrackID();
    int           parentID = aTrack->GetParentID();
    G4ThreeVector position = aTrack->GetPosition();
    G4ThreeVector momentum = aTrack->GetMomentum();
    G4String      partType = aTrack->GetDefinition()->GetParticleType();
    G4String      partSubType = aTrack->GetDefinition()->GetParticleSubType();
    int    partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
#ifdef ddebug
    bool   isPDGStable = aTrack->GetDefinition()->GetPDGStable();
#endif
    double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
    double gammaFactor = aStep->GetPreStepPoint()->GetGamma();

    if (!pvFound) { //search for v1
      double stepDeltaEnergy = aStep->GetDeltaEnergy ();
      double kinEnergy = aTrack->GetKineticEnergy ();
      
      // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
      if (trackID == 1 && parentID == 0 && 
          ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
        int pvType = -1;
        if (kinEnergy == 0.) {
          pvType = 0;
        } else {
          if (fabs (stepDeltaEnergy / kinEnergy) > 0.1) pvType = 1;
        }
        pvFound    = true;
        pvPosition = position;
        pvMomentum = momentum;
        // Rotated coord.system:
        pvUVW      = (*beamline_RM)*(pvPosition);

        //Volume name requires some checks:
        G4String thePostPVname = "NoName";
        G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
        if (thePostPoint) {
          thePostStepPoint = thePostPoint->GetPosition();
          G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
          if (thePostPV) thePostPVname = thePostPV->GetName ();
        }
#ifdef ddebug
        LogDebug("HcalTBSim") << "HcalTB06Analysis:: V1 found at: " 
                              << thePostStepPoint << " G4VPhysicalVolume: " 
                              << thePostPVname;
#endif      
        LogDebug("HcalTBSim") << "HcalTB06Analysis::fill_v1Pos: Primary Track "
                              << "momentum: " << pvMomentum << " psoition " 
                              << pvPosition << " u/v/w " << pvUVW;
      }
    } else { 
      // watch for secondaries originating @v1, including the surviving primary
      if ((trackID != 1 && parentID == 1 &&
           (aTrack->GetCurrentStepNumber () == 1) && 
           (thePreStepPoint == pvPosition)) || 
          (trackID == 1 && thePreStepPoint == pvPosition)) {
#ifdef ddebug
        LogDebug("HcalTBSim") << "HcalTB06Analysis::A secondary...  PDG:" 
                              << partPDGEncoding << " TrackID:" << trackID
                              << " ParentID:" << parentID << " stable: "  
                              << isPDGStable << " Tau: " << pDGlifetime 
                              << " cTauGamma=" 
                              << c_light*pDGlifetime*gammaFactor*1000.
                              << "um" << " GammaFactor: " << gammaFactor;
#endif      
        secTrackID.push_back(trackID);
        secPartID.push_back(partPDGEncoding);
        secMomentum.push_back(momentum);
        secEkin.push_back(aTrack->GetKineticEnergy());

        // Check for short-lived secondaries: cTauGamma<100um
        double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
        if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {//short-lived secondary
          shortLivedSecondaries.push_back(trackID);
      } else {//normal secondary - enter into the V1-calorimetric tree
        //          histos->fill_v1cSec (aTrack);
      }
      }
      // Also watch for tertiary particles coming from 
      // short-lived secondaries from V1
      if (aTrack->GetCurrentStepNumber() == 1) {
        if (shortLivedSecondaries.size() > 0) {
          int pid = parentID;
          std::vector<int>::iterator pos1= shortLivedSecondaries.begin();
          std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
          std::vector<int>::iterator pos;
          for (pos = pos1; pos != pos2; pos++) {
            if (*pos == pid) {//ParentID is on the list of short-lived 
              // secondary 
#ifdef ddebug
              LogDebug("HcalTBSim") << "HcalTB06Analysis:: A tertiary...  PDG:"
                                    << partPDGEncoding << " TrackID:" <<trackID
                                    << " ParentID:" << parentID << " stable: "
                                    << isPDGStable << " Tau: " << pDGlifetime
                                    << " cTauGamma=" 
                                    << c_light*pDGlifetime*gammaFactor*1000. 
                                    << "um GammaFactor: " << gammaFactor;
#endif
            }
          }
        }
      }
    }
  }
}

Member Data Documentation

G4RotationMatrix* HcalTB06Analysis::beamline_RM [private]

Definition at line 85 of file HcalTB06Analysis.h.

Referenced by HcalTB06Analysis().

double HcalTB06Analysis::beamOffset [private]

Definition at line 82 of file HcalTB06Analysis.h.

Referenced by HcalTB06Analysis().

int HcalTB06Analysis::count [private]

Definition at line 88 of file HcalTB06Analysis.h.

Referenced by init(), update(), and ~HcalTB06Analysis().

std::vector<CaloHit> HcalTB06Analysis::ecalHitCache [private]

Definition at line 93 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), fillEvent(), and finalAnalysis().

double HcalTB06Analysis::eecals [private]

Definition at line 95 of file HcalTB06Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

double HcalTB06Analysis::ehcals [private]

Definition at line 95 of file HcalTB06Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

double HcalTB06Analysis::etaInit [private]

Definition at line 92 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), fillEvent(), and finalAnalysis().

double HcalTB06Analysis::etots [private]

Definition at line 95 of file HcalTB06Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

int HcalTB06Analysis::evNum [private]

Definition at line 98 of file HcalTB06Analysis.h.

Referenced by fillBuffer(), fillEvent(), init(), and update().

std::vector<CaloHit> HcalTB06Analysis::hcalHitCache [private]

Definition at line 94 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), fillEvent(), and finalAnalysis().

std::vector<CaloHit> HcalTB06Analysis::hcalHitLayer [private]

Definition at line 94 of file HcalTB06Analysis.h.

Definition at line 79 of file HcalTB06Analysis.h.

Referenced by finalAnalysis(), HcalTB06Analysis(), and ~HcalTB06Analysis().

int HcalTB06Analysis::iceta [private]

Definition at line 83 of file HcalTB06Analysis.h.

Referenced by HcalTB06Analysis().

int HcalTB06Analysis::icphi [private]

Definition at line 83 of file HcalTB06Analysis.h.

Referenced by HcalTB06Analysis().

std::vector<std::string> HcalTB06Analysis::names [private]

Definition at line 84 of file HcalTB06Analysis.h.

Referenced by fillBuffer(), and HcalTB06Analysis().

Definition at line 91 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), and fillEvent().

Definition at line 91 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), and fillEvent().

double HcalTB06Analysis::phiInit [private]

Definition at line 92 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), fillEvent(), and finalAnalysis().

double HcalTB06Analysis::pInit [private]

Definition at line 92 of file HcalTB06Analysis.h.

Referenced by clear(), fillBuffer(), fillEvent(), and finalAnalysis().

bool HcalTB06Analysis::pvFound [private]

Definition at line 97 of file HcalTB06Analysis.h.

Referenced by clear(), and update().

G4ThreeVector HcalTB06Analysis::pvMomentum [private]

Definition at line 99 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

G4ThreeVector HcalTB06Analysis::pvPosition [private]

Definition at line 99 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

int HcalTB06Analysis::pvType [private]

Definition at line 98 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

G4ThreeVector HcalTB06Analysis::pvUVW [private]

Definition at line 99 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

std::vector<double> HcalTB06Analysis::secEkin [private]

Definition at line 102 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

std::vector<G4ThreeVector> HcalTB06Analysis::secMomentum [private]

Definition at line 101 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

std::vector<int> HcalTB06Analysis::secPartID [private]

Definition at line 100 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

std::vector<int> HcalTB06Analysis::secTrackID [private]

Definition at line 100 of file HcalTB06Analysis.h.

Referenced by clear(), fillEvent(), and update().

std::vector<int> HcalTB06Analysis::shortLivedSecondaries [private]

Definition at line 103 of file HcalTB06Analysis.h.

Referenced by clear(), and update().