CMS 3D CMS Logo

HcalTB04Analysis Class Reference

Description: Analysis of 2004 Hcal Test beam simulation. More...

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

Inheritance diagram for HcalTB04Analysis:

SimProducer Observer< const BeginOfRun * > Observer< const BeginOfEvent * > Observer< const EndOfEvent * > Observer< const G4Step * > SimWatcher

List of all members.

Public Member Functions

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

Private Member Functions

void clear ()
void fillBuffer (const EndOfEvent *evt)
void fillEvent (PHcalTB04Info &)
void finalAnalysis ()
 HcalTB04Analysis (const HcalTB04Analysis &)
void init ()
const HcalTB04Analysis & operator= (const HcalTB04Analysis &)
void qieAnalysis ()
double scale (int det, int layer)
double timeOfFlight (int det, int layer, double eta)
int unitID (uint32_t id)
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.
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 xtalAnalysis ()

Private Attributes

G4RotationMatrix * beamline_RM
double beamOffset
int count
std::vector< CaloHit > ecalHitCache
double ecalNoise
double eecalq
double eecals
double ehcalq
double ehcals
std::vector< double > enois
std::vector< double > eqeta
std::vector< double > eqie
std::vector< double > eqlay
std::vector< double > eqphi
std::vector< double > eseta
std::vector< double > esime
std::vector< double > esimh
std::vector< double > eslay
std::vector< double > esphi
double etaInit
double etotq
double etots
int evNum
std::vector< CaloHit > hcalHitCache
std::vector< CaloHit > hcalHitLayer
bool hcalOnly
HcalTB04Histo * histo
int iceta
int icphi
std::vector< uint32_t > idEcal
std::vector< int > idHcal
std::vector< uint32_t > idTower
std::vector< int > idXtal
int mode
HcalQie * myQie
std::vector< std::string > names
int nCrystal
int nPrimary
int nTower
int particleType
double phiInit
double pInit
bool pvFound
G4ThreeVector pvMomentum
G4ThreeVector pvPosition
int pvType
G4ThreeVector pvUVW
double scaleHB0
double scaleHB16
double scaleHE0
double scaleHO
std::vector< double > secEkin
std::vector< G4ThreeVector > secMomentum
std::vector< int > secPartID
std::vector< int > secTrackID
std::vector< int > shortLivedSecondaries
int type

Detailed Description

Description: Analysis of 2004 Hcal Test beam simulation.

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

Definition at line 48 of file HcalTB04Analysis.h.

Constructor & Destructor Documentation

HcalTB04Analysis::HcalTB04Analysis ( const edm::ParameterSet &  p  ) 

Definition at line 57 of file HcalTB04Analysis.cc.

References beamline_RM, beamOffset, ecalNoise, funct::exp(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hcalOnly, histo, iceta, icphi, init(), int, mode, myQie, names, scaleHB0, scaleHB16, scaleHE0, and scaleHO.

00057                                                           : myQie(0),
00058                                                                 histo(0) {
00059 
00060   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("HcalTB04Analysis");
00061   hcalOnly       = m_Anal.getParameter<bool>("HcalOnly");
00062   mode           = m_Anal.getParameter<int>("Mode");
00063   type           = m_Anal.getParameter<int>("Type");
00064   ecalNoise      = m_Anal.getParameter<double>("EcalNoise");
00065   scaleHB0       = m_Anal.getParameter<double>("ScaleHB0");
00066   scaleHB16      = m_Anal.getParameter<double>("ScaleHB16");
00067   scaleHO        = m_Anal.getParameter<double>("ScaleHO");
00068   scaleHE0       = m_Anal.getParameter<double>("ScaleHE0");
00069   names          = m_Anal.getParameter<std::vector<std::string> >("Names");
00070   beamOffset     =-m_Anal.getUntrackedParameter<double>("BeamPosition",0.0)*cm;
00071   double fMinEta = m_Anal.getUntrackedParameter<double>("MinEta",-5.5);
00072   double fMaxEta = m_Anal.getUntrackedParameter<double>("MaxEta",5.5);
00073   double fMinPhi = m_Anal.getUntrackedParameter<double>("MinPhi",-3.14159265358979323846);
00074   double fMaxPhi = m_Anal.getUntrackedParameter<double>("MaxPhi", 3.14159265358979323846);
00075   double beamEta = (fMaxEta+fMinEta)/2.;
00076   double beamPhi = (fMaxPhi+fMinPhi)/2.;
00077   double beamThet= 2*atan(exp(-beamEta));
00078   if (beamPhi < 0) beamPhi += twopi;
00079   iceta          = (int)(beamEta/0.087) + 1;
00080   icphi          = (int)(fabs(beamPhi)/0.087) + 5;
00081   if (icphi > 72) icphi -= 73;
00082 
00083   produces<PHcalTB04Info>();
00084 
00085   beamline_RM = new G4RotationMatrix;
00086   beamline_RM->rotateZ(-beamPhi);
00087   beamline_RM->rotateY(-beamThet);
00088  
00089   edm::LogInfo("HcalTBSim") << "HcalTB04:: Initialised as observer of BeginOf"
00090                             << "Job/BeginOfRun/BeginOfEvent/G4Step/EndOfEvent"
00091                             << " with Parameter values:\n \thcalOnly = " 
00092                             << hcalOnly << "\tecalNoise = " << ecalNoise
00093                             << "\n\tMode = " << mode << " (0: HB2 Standard; "
00094                             << "1:HB2 Segmented)" << "\tType = " << type 
00095                             << " (0: HB; 1 HE; 2 HB+HE)\n\tbeamOffset = " 
00096                             << beamOffset << "\ticeta = " << iceta 
00097                             << "\ticphi = " << icphi << "\n\tbeamline_RM = "
00098                             << *beamline_RM;
00099 
00100   init();
00101 
00102   myQie  = new HcalQie(p);
00103   histo  = new HcalTB04Histo(m_Anal);
00104 }

HcalTB04Analysis::~HcalTB04Analysis (  )  [virtual]

Definition at line 106 of file HcalTB04Analysis.cc.

References count, histo, and myQie.

00106                                     {
00107 
00108   edm::LogInfo("HcalTBSim") << "\n -------->  Total number of selected entries"
00109                             << " : " << count << "\nPointers:: QIE " << myQie
00110                             << " Histo " << histo;
00111   if (myQie)   {
00112     delete myQie;
00113     myQie  = 0;
00114   }
00115   if (histo)   {
00116     delete histo;
00117     histo  = 0;
00118   }
00119 }

HcalTB04Analysis::HcalTB04Analysis ( const HcalTB04Analysis &   )  [private]

Member Function Documentation

void HcalTB04Analysis::clear ( void   )  [private]

Definition at line 951 of file HcalTB04Analysis.cc.

References ecalHitCache, enois, eqie, esime, esimh, etaInit, hcalHitCache, hcalHitLayer, i, nCrystal, nPrimary, nTower, particleType, phiInit, pInit, pvFound, pvMomentum, pvPosition, pvType, pvUVW, secEkin, secMomentum, secPartID, secTrackID, and shortLivedSecondaries.

Referenced by init(), and update().

00951                             {
00952   pvFound = false;
00953   pvType  =-2;
00954   pvPosition = G4ThreeVector();
00955   pvMomentum = G4ThreeVector();
00956   pvUVW      = G4ThreeVector();
00957   secTrackID.clear();
00958   secPartID.clear();
00959   secMomentum.clear();
00960   secEkin.clear();
00961   shortLivedSecondaries.clear();
00962 
00963   ecalHitCache.erase(ecalHitCache.begin(), ecalHitCache.end()); 
00964   hcalHitCache.erase(hcalHitCache.begin(), hcalHitCache.end()); 
00965   hcalHitLayer.erase(hcalHitLayer.begin(), hcalHitLayer.end());
00966   nPrimary = particleType = 0;
00967   pInit = etaInit = phiInit = 0;
00968 
00969   esimh.clear();
00970   eqie.clear();
00971   esimh.reserve(nTower);
00972   eqie.reserve(nTower);
00973   for (int i=0; i<nTower; i++) {
00974     esimh.push_back(0.);
00975     eqie.push_back(0.);
00976   }
00977   esime.clear();
00978   enois.clear();
00979   esime.reserve(nCrystal);
00980   enois.reserve(nCrystal);
00981   for (int i=0; i<nCrystal; i++) {
00982     esime.push_back(0.);
00983     enois.push_back(0.);
00984   }
00985 }

void HcalTB04Analysis::fillBuffer ( const EndOfEvent *  evt  )  [private]

Definition at line 398 of file HcalTB04Analysis.cc.

References e, ecalHitCache, eta, etaInit, evNum, CaloG4Hit::getEnergyDeposit(), CaloG4Hit::getEntry(), CaloG4Hit::getTimeSlice(), CaloG4Hit::getTrackID(), HcalTBNumberingScheme::getUnitID(), CaloG4Hit::getUnitID(), group, hcalHitCache, hcalHitLayer, i, int, j, funct::log(), LogDebug, max, min, mode, names, npart, nPrimary, p, particleType, phi, phiInit, pInit, funct::pow(), scale(), python::multivaluedict::sort(), funct::sqrt(), funct::tan(), theta, timeOfFlight(), unitID(), HcalTestNumbering::unpackHcalIndex(), and z.

Referenced by update().

00398                                                         {
00399 
00400   std::vector<CaloHit> hhits, hhitl;
00401   int                  idHC, j;
00402   CaloG4HitCollection* theHC;
00403   std::map<int,float,std::less<int> > primaries;
00404   double               etot1=0, etot2=0;
00405 
00406   // Look for the Hit Collection of HCal
00407   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
00408   std::string sdName = names[0];
00409   idHC  = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
00410   theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
00411   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName
00412                         << " of ID " << idHC << " is obtained at " << theHC;
00413 
00414   if (idHC >= 0 && theHC > 0) {
00415     hhits.reserve(theHC->entries());
00416     hhitl.reserve(theHC->entries());
00417     for (j = 0; j < theHC->entries(); j++) {
00418       CaloG4Hit* aHit = (*theHC)[j]; 
00419       double e        = aHit->getEnergyDeposit()/GeV;
00420       double time     = aHit->getTimeSlice();
00421       math::XYZPoint pos  = aHit->getEntry();
00422       unsigned int id = aHit->getUnitID();
00423       double theta    = pos.theta();
00424       double eta      = -log(tan(theta * 0.5));
00425       double phi      = pos.phi();
00426       int det, z, group, ieta, iphi, layer;
00427       HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
00428       double jitter   = time-timeOfFlight(det,layer,eta);
00429       if (jitter<0) jitter = 0;
00430       if (e < 0 || e > 1.) e = 0;
00431       double escl     = e * scale(det,layer);
00432       unsigned int idx= HcalTBNumberingScheme::getUnitID(id,mode);
00433       CaloHit hit(det,layer,escl,eta,phi,jitter,idx);
00434       hhits.push_back(hit);
00435       CaloHit hitl(det,layer,escl,eta,phi,jitter,id);
00436       hhitl.push_back(hitl);
00437       primaries[aHit->getTrackID()]+= e;
00438       etot1 += escl;
00439 #ifdef ddebug
00440       LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit i/p " << j 
00441                             << "  ID 0x" << std::hex << id << " 0x" << idx  
00442                             << std::dec << " time " << std::setw(6) << time 
00443                             << " "  << std::setw(6) << jitter << " theta " 
00444                             << std::setw(8) << theta << " eta " << std::setw(8)
00445                             << eta << " phi " << std::setw(8) << phi << " e " 
00446                             << std::setw(8) << e << " " << std::setw(8) <<escl;
00447 #endif
00448     }
00449   }
00450 
00451   // Add hits in the same channel within same time slice
00452   std::vector<CaloHit>::iterator itr;
00453   int nHit = hhits.size();
00454   std::vector<CaloHit*> hits(nHit);
00455   for (j = 0, itr = hhits.begin(); itr != hhits.end(); j++, itr++) {
00456     hits[j] = &hhits[j];
00457   }
00458   sort(hits.begin(),hits.end(),CaloHitIdMore());
00459   std::vector<CaloHit*>::iterator k1, k2;
00460   int nhit = 0;
00461   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
00462     int      det    = (**k1).det();
00463     int      layer  = (**k1).layer();
00464     double   ehit   = (**k1).e();
00465     double   eta    = (**k1).eta();
00466     double   phi    = (**k1).phi();
00467     double   jitter = (**k1).t();
00468     uint32_t unitID = (**k1).id();
00469     int      jump  = 0;
00470     for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
00471            unitID==(**k2).id(); k2++) {
00472       ehit += (**k2).e();
00473       jump++;
00474     }
00475     nhit++;
00476     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
00477     hcalHitCache.push_back(hit);
00478     etot2 += ehit;
00479     k1    += jump;
00480 #ifdef ddebug
00481     LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhit 
00482                           << "  ID 0x" << std::hex  << unitID  << std::dec 
00483                           << " time " << std::setw(6) << jitter << " eta "
00484                           << std::setw(8) << eta << " phi " << std::setw(8) 
00485                           << phi  << " e " << std::setw(8) << ehit;
00486 #endif
00487   }
00488   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " HCal hits"
00489                         << " from " << nHit << " input hits E(Hcal) " << etot1 
00490                         << " " << etot2;
00491 
00492   //Repeat for Hit in each layer (hhits and hhitl sizes are the same)
00493   for (j = 0, itr = hhitl.begin(); itr != hhitl.end(); j++, itr++) {
00494     hits[j] = &hhitl[j];
00495   }
00496   sort(hits.begin(),hits.end(),CaloHitIdMore());
00497   int    nhitl = 0;
00498   double etotl = 0;
00499   for (k1 = hits.begin(); k1 != hits.end(); k1++) {
00500     int      det    = (**k1).det();
00501     int      layer  = (**k1).layer();
00502     double   ehit   = (**k1).e();
00503     double   eta    = (**k1).eta();
00504     double   phi    = (**k1).phi();
00505     double   jitter = (**k1).t();
00506     uint32_t unitID = (**k1).id();
00507     int      jump  = 0;
00508     for (k2 = k1+1; k2 != hits.end() && fabs(jitter-(**k2).t())<1 &&
00509            unitID==(**k2).id(); k2++) {
00510       ehit += (**k2).e();
00511       jump++;
00512     }
00513     nhitl++;
00514     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
00515     hcalHitLayer.push_back(hit);
00516     etotl += ehit;
00517     k1    += jump;
00518 #ifdef ddebug
00519     LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hcal Hit store " << nhitl 
00520                           << "  ID 0x" << std::hex << unitID  << std::dec 
00521                           << " time " << std::setw(6) << jitter << " eta "
00522                           << std::setw(8) << eta << " phi " << std::setw(8) 
00523                           << phi << " e " << std::setw(8) << ehit;
00524 #endif
00525   }
00526   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhitl << " HCal "
00527                         << "hits from " << nHit << " input hits E(Hcal) " 
00528                         << etot1 << " " << etotl;
00529   
00530   // Look for the Hit Collection of ECal
00531   std::vector<CaloHit> ehits;
00532   sdName= names[1];
00533   idHC  = G4SDManager::GetSDMpointer()->GetCollectionID(sdName);
00534   theHC = (CaloG4HitCollection*) allHC->GetHC(idHC);
00535   etot1 = etot2 = 0;
00536   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Hit Collection for " << sdName
00537                         << " of ID " << idHC << " is obtained at " << theHC;
00538   if (idHC >= 0 && theHC > 0) {
00539     ehits.reserve(theHC->entries());
00540     for (j = 0; j < theHC->entries(); j++) {
00541       CaloG4Hit* aHit = (*theHC)[j]; 
00542       double e        = aHit->getEnergyDeposit()/GeV;
00543       double time     = aHit->getTimeSlice();
00544       math::XYZPoint pos  = aHit->getEntry();
00545       unsigned int id = aHit->getUnitID();
00546       double theta    = pos.theta();
00547       double eta      = -log(tan(theta * 0.5));
00548       double phi      = pos.phi();
00549       if (e < 0 || e > 100000.) e = 0;
00550       int det, z, group, ieta, iphi, layer;
00551       HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
00552       CaloHit hit(det,0,e,eta,phi,time,id);
00553       ehits.push_back(hit);
00554       primaries[aHit->getTrackID()]+= e;
00555       etot1 += e;
00556 #ifdef ddebug
00557       LogDebug("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit i/p " << j 
00558                             << "  ID 0x" << std::hex << id  << std::dec 
00559                             << " time " << std::setw(6) << time << " theta " 
00560                             << std::setw(8) << theta  << " eta " <<std::setw(8)
00561                             << eta  << " phi " << std::setw(8) << phi << " e "
00562                             << std::setw(8) << e;
00563 #endif
00564     }
00565   }
00566 
00567   // Add hits in the same channel within same time slice
00568   nHit = ehits.size();
00569   std::vector<CaloHit*> hite(nHit);
00570   for (j = 0, itr = ehits.begin(); itr != ehits.end(); j++, itr++) {
00571     hite[j] = &ehits[j];
00572   }
00573   sort(hite.begin(),hite.end(),CaloHitIdMore());
00574   nhit = 0;
00575   for (k1 = hite.begin(); k1 != hite.end(); k1++) {
00576     int      det    = (**k1).det();
00577     int      layer  = (**k1).layer();
00578     double   ehit   = (**k1).e();
00579     double   eta    = (**k1).eta();
00580     double   phi    = (**k1).phi();
00581     double   jitter = (**k1).t();
00582     uint32_t unitID = (**k1).id();
00583     int      jump  = 0;
00584     for (k2 = k1+1; k2 != hite.end() && fabs(jitter-(**k2).t())<1 &&
00585            unitID==(**k2).id(); k2++) {
00586       ehit += (**k2).e();
00587       jump++;
00588     }
00589     nhit++;
00590     CaloHit hit(det, layer, ehit, eta, phi, jitter, unitID);
00591     ecalHitCache.push_back(hit);
00592     etot2 += ehit;
00593     k1    += jump;
00594 #ifdef ddebug
00595     LogDebug("HcalTBSim") << "HcalTB04Analysis:: Ecal Hit store " << nhit
00596                           << "  ID 0x" << std::hex << unitID  << std::dec 
00597                           << " time " << std::setw(6) << jitter << " eta "
00598                           << std::setw(8) << eta << " phi " << std::setw(8)
00599                           << phi << " e " << std::setw(8) << ehit;
00600 #endif
00601   }
00602   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Stores " << nhit << " ECal hits"
00603                         << " from " << nHit << " input hits E(Ecal) " << etot1 
00604                         << " " << etot2;
00605 
00606   // Find Primary info:
00607   nPrimary    = (int)(primaries.size());
00608   int trackID = 0;
00609   G4PrimaryParticle* thePrim=0;
00610   int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
00611   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Event has " << nvertex 
00612                         << " verteices";
00613   if (nvertex<=0)
00614     edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERROR: no "
00615                               << "vertex found for event " << evNum;
00616 
00617   for (int i = 0 ; i<nvertex; i++) {
00618     G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
00619     if (avertex == 0) {
00620       edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: pointer "
00621                                 << "to vertex = 0 for event " << evNum;
00622     } else {
00623       LogDebug("HcalTBSim") << "HcalTB04Analysis::Vertex number :" << i << " "
00624                             << avertex->GetPosition();
00625       int npart = avertex->GetNumberOfParticle();
00626       if (npart == 0)
00627         edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::End Of Event ERR: "
00628                                      << "no primary!";
00629       if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
00630     }
00631   }
00632     
00633   if (thePrim != 0) {
00634     double px = thePrim->GetPx();
00635     double py = thePrim->GetPy();
00636     double pz = thePrim->GetPz();
00637     double p  = std::sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
00638     pInit     = p/GeV;
00639     if (p==0) 
00640       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis:: EndOfEvent ERR: "
00641                                    << "primary has p=0 ";
00642     else {
00643       double costheta = pz/p;
00644       double theta = acos(std::min(std::max(costheta,-1.),1.));
00645       etaInit = -log(tan(theta/2));
00646       if (px != 0 || py != 0) phiInit = atan2(py,px);  
00647     }
00648     particleType = thePrim->GetPDGcode();
00649   } else 
00650     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::EndOfEvent ERR: could "
00651                                  << "not find primary";
00652 
00653 }

void HcalTB04Analysis::fillEvent ( PHcalTB04Info &  product  )  [private]

Definition at line 877 of file HcalTB04Analysis.cc.

References ecalHitCache, eecalq, eecals, ehcalq, ehcals, enois, eqeta, eqie, eqlay, eqphi, eseta, esime, esimh, eslay, esphi, etaInit, etotq, etots, evNum, group, hcalHitCache, i, idHcal, idXtal, LogDebug, nPrimary, particleType, phiInit, pInit, pvMomentum, pvPosition, pvType, pvUVW, PHcalTB04Info::saveHit(), secEkin, secMomentum, secPartID, secTrackID, PHcalTB04Info::setEdep(), PHcalTB04Info::setEdepHcal(), PHcalTB04Info::setIDs(), PHcalTB04Info::setLongProf(), PHcalTB04Info::setPrimary(), PHcalTB04Info::setTrnsProf(), PHcalTB04Info::setVtxPrim(), PHcalTB04Info::setVtxSec(), HcalTestNumbering::unpackHcalIndex(), x, y, and z.

Referenced by produce().

00877                                                         {
00878 
00879   //Setup the ID's
00880   product.setIDs(idHcal, idXtal);
00881 
00882   //Beam Information
00883   product.setPrimary(nPrimary, particleType, pInit, etaInit, phiInit);
00884 
00885   //Energy deposits in the crystals and towers
00886   product.setEdepHcal(esimh, eqie);
00887   product.setEdepHcal(esime, enois);
00888 
00889   // Total Energy
00890   product.setEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
00891 
00892   // Lateral Profile
00893   product.setTrnsProf(eseta,eqeta,esphi,eqphi);
00894 
00895   // Longitudianl profile
00896   product.setLongProf(eslay, eqlay);
00897 
00898   //Save Hits
00899   int i, nhit=0;
00900   std::vector<CaloHit>::iterator itr;
00901   for (i=0, itr=ecalHitCache.begin(); itr!=ecalHitCache.end(); i++,itr++) {
00902     uint32_t id = itr->id();
00903     int det, z, group, ieta, iphi, lay;
00904     HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,lay);
00905     product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
00906     nhit++;
00907 #ifdef debug
00908     LogDebug("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) 
00909                           << i+1 << " ID 0x" << std::hex << group << std::dec 
00910                           << " "  << std::setw(2) << det << " " << std::setw(2)
00911                           << lay  << " " << std::setw(1) << z << " " 
00912                           << std::setw(3) << ieta << " " << std::setw(3) <<iphi
00913                           << " T " << std::setw(6) << itr->t() << " E " 
00914                           << std::setw(6) << itr->e();
00915 #endif
00916   }
00917   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit 
00918                         << " hits from Crystals";
00919   int hit = nhit;
00920   nhit = 0;
00921 
00922   for (i=hit, itr=hcalHitCache.begin(); itr!=hcalHitCache.end(); i++,itr++) {
00923     uint32_t id = itr->id();
00924     int det, z, group, ieta, iphi, lay;
00925     HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,lay);
00926     product.saveHit(det, lay, ieta, iphi, itr->e(), itr->t());
00927     nhit++;
00928 #ifdef debug
00929     LogDebug("HcalTBSim") << "HcalTB04Analysis:: Save Hit " << std::setw(3) 
00930                           << i+1 << " ID 0x" << std::hex << group << std::dec 
00931                           << " "  << std::setw(2) << det << " " << std::setw(2)
00932                           << lay  << " " << std::setw(1) << z << " " 
00933                           << std::setw(3) << ieta << " " << std::setw(3) <<iphi
00934                           << " T " << std::setw(6) << itr->t() << " E " 
00935                           << std::setw(6) << itr->e();
00936 #endif
00937   }
00938   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Saves " << nhit 
00939                         << " hits from HCal";
00940 
00941   //Vertex associated quantities
00942   product.setVtxPrim(evNum, pvType, pvPosition.x(), pvPosition.y(), 
00943                      pvPosition.z(), pvUVW.x(), pvUVW.y(), pvUVW.z(),
00944                      pvMomentum.x(), pvMomentum.y(), pvMomentum.z());
00945   for (unsigned int i = 0; i < secTrackID.size(); i++) {
00946     product.setVtxSec(secTrackID[i], secPartID[i], secMomentum[i].x(),
00947                       secMomentum[i].y(), secMomentum[i].z(), secEkin[i]);
00948   }
00949 }

void HcalTB04Analysis::finalAnalysis (  )  [private]

Definition at line 775 of file HcalTB04Analysis.cc.

References e1, e2, eecalq, eecals, ehcalq, ehcals, enois, eqeta, eqie, eqlay, eqphi, eseta, esime, esimh, eslay, esphi, etaInit, etotq, etots, HcalTB04Histo::fillEdep(), HcalTB04Histo::fillLongProf(), HcalTB04Histo::fillPrimary(), HcalTB04Histo::fillTrnsProf(), group, histo, i, iceta, icphi, id, idTower, LogDebug, nCrystal, nTower, phiInit, pInit, HcalTestNumbering::unpackHcalIndex(), and z.

Referenced by update().

00775                                      {
00776 
00777   //Beam Information
00778   histo->fillPrimary(pInit, etaInit, phiInit);
00779 
00780   // Total Energy
00781   eecals = ehcals = eecalq = ehcalq = 0.;
00782   for (int i=0; i<nTower; i++) {
00783     ehcals += esimh[i];
00784     ehcalq += eqie[i];
00785   }
00786   for (int i=0; i<nCrystal; i++) {
00787     eecals += esime[i];
00788     eecalq += enois[i];
00789   }
00790   etots = eecals + ehcals;
00791   etotq = eecalq + ehcalq;
00792   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Energy deposit at Sim Level "
00793                         << "(Total) " << etots << " (ECal) " << eecals 
00794                         << " (HCal) " << ehcals << "\nHcalTB04Analysis:: "
00795                         << "Energy deposit at Qie Level (Total) " << etotq
00796                         << " (ECal) " << eecalq << " (HCal) " << ehcalq;
00797   histo->fillEdep(etots, eecals, ehcals, etotq, eecalq, ehcalq);
00798 
00799   // Lateral Profile
00800   for (int i=0; i<5; i++) {
00801     eseta[i] = 0.;
00802     eqeta[i] = 0.;
00803   }
00804   for (int i=0; i<3; i++) {
00805     esphi[i] = 0.;
00806     eqphi[i] = 0.;
00807   }
00808   double e1=0, e2=0;
00809   unsigned int id;
00810   for (int i=0; i<nTower; i++) {
00811     int det, z, group, ieta, iphi, layer;
00812     id = idTower[i];
00813     HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
00814     iphi -= (icphi - 1);
00815     if (icphi > 4) {
00816       if (ieta == 0) ieta = 2;
00817       else           ieta =-1;
00818     } else {
00819       ieta = ieta - iceta + 2;
00820     }
00821     if (iphi >= 0 && iphi < 3 && ieta >= 0 && ieta < 5) {
00822       eseta[ieta] += esimh[i];
00823       esphi[iphi] += esimh[i];
00824       e1          += esimh[i];
00825       eqeta[ieta] += eqie[i];
00826       eqphi[iphi] += eqie[i];
00827       e2          += eqie[i];
00828     }
00829   }
00830   for (int i=0; i<3; i++) {
00831     if (e1>0) esphi[i] /= e1;
00832     if (e2>0) eqphi[i] /= e2;
00833   }
00834   for (int i=0; i<5; i++) {
00835     if (e1>0) eseta[i] /= e1;
00836     if (e2>0) eqeta[i] /= e2;
00837   }
00838   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Eta and" 
00839                         << " Phi (Sim/Qie)";
00840   for (int i=0; i<5; i++) 
00841     LogDebug("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Eta Sim = "
00842                           << eseta[i] << " Qie = " << eqeta[i] << " Phi Sim = "
00843                           << esphi[i] << " Qie = " << eqphi[i];
00844   histo->fillTrnsProf(eseta,eqeta,esphi,eqphi);
00845 
00846   // Longitudianl profile
00847   for (int i=0; i<20; i++) {
00848     eslay[i] = 0.;
00849     eqlay[i] = 0.;
00850   }
00851   e1=0; e2=0;
00852   for (int i=0; i<nTower; i++) {
00853     int det, z, group, ieta, iphi, layer;
00854     id = idTower[i];
00855     HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,layer);
00856     iphi  -= (icphi - 1);
00857     layer -= 1;
00858     if (iphi >= 0 && iphi < 3 && layer >= 0 && layer < 20) {
00859       eslay[layer] += esimh[i];
00860       e1           += esimh[i];
00861       eqlay[layer] += eqie[i];
00862       e2           += eqie[i];
00863     }
00864   }
00865   for (int i=0; i<20; i++) {
00866     if (e1>0) eslay[i] /= e1;
00867     if (e2>0) eqlay[i] /= e2;
00868   }
00869   LogDebug("HcalTBSim") << "HcalTB04Analysis:: Energy fraction along Layer";
00870   for (int i=0; i<20; i++)
00871     LogDebug("HcalTBSim") << "HcalTB04Analysis:: [" << i << "] Sim = " 
00872                           << eslay[i] << " Qie = " << eqlay[i];
00873   histo->fillLongProf(eslay, eqlay);
00874 }

void HcalTB04Analysis::init ( void   )  [private]

Definition at line 132 of file HcalTB04Analysis.cc.

References clear(), count, eqeta, eqlay, eqphi, eseta, eslay, esphi, evNum, HcalTBNumberingScheme::getUnitIDs(), hcalOnly, i, idEcal, idHcal, idTower, idXtal, LogDebug, mode, nCrystal, nTower, HcalTestNumbering::packHcalIndex(), and unitID().

Referenced by HcalTB04Analysis().

00132                             {
00133 
00134   idTower = HcalTBNumberingScheme::getUnitIDs(type, mode);
00135   nTower  = idTower.size();
00136   edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Save information from " 
00137                             << nTower << " HCal towers";
00138   idHcal.reserve(nTower);
00139   for (int i=0; i<nTower; i++) {
00140     int id = unitID(idTower[i]);
00141     idHcal.push_back(id);
00142     LogDebug("HcalTBSim") << "\tTower[" << i << "] Original " << std::hex
00143                           << idTower[i] << " Stored " << idHcal[i] << std::dec;
00144   }
00145 
00146   if (!hcalOnly) {
00147     int  det = 10;
00148     uint32_t id1;
00149     nCrystal = 0;
00150     for (int lay=1; lay<8; lay++) {
00151       for (int icr=1; icr<8; icr++) {
00152         id1    = HcalTestNumbering::packHcalIndex(det,0,1,icr,lay,1);
00153         int id = unitID(id1);
00154         idEcal.push_back(id1);
00155         idXtal.push_back(id);
00156         nCrystal++;
00157       }
00158     }
00159     edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Save information from " 
00160                               << nCrystal << " ECal Crystals";
00161     for (int i=0; i<nCrystal; i++) {
00162       LogDebug("HcalTBSim") << "\tCrystal[" << i << "] Original " << std::hex
00163                             << idEcal[i] << " Stored " << idXtal[i] <<std::dec;
00164     }
00165   }
00166   // Profile vectors
00167   eseta.reserve(5);
00168   eqeta.reserve(5); 
00169   esphi.reserve(3);
00170   eqphi.reserve(3);
00171   eslay.reserve(20);
00172   eqlay.reserve(20);
00173   for (int i=0; i<5; i++) {
00174     eseta.push_back(0.);
00175     eqeta.push_back(0.);
00176   }
00177   for (int i=0; i<3; i++) {
00178     esphi.push_back(0.);
00179     eqphi.push_back(0.);
00180   }
00181   for (int i=0; i<20; i++) {
00182     eslay.push_back(0.);
00183     eqlay.push_back(0.);
00184   }
00185 
00186   // counter 
00187   count = 0;
00188   evNum = 0;
00189   clear();
00190 }

const HcalTB04Analysis& HcalTB04Analysis::operator= ( const HcalTB04Analysis &   )  [private]

void HcalTB04Analysis::produce ( edm::Event &  e,
const edm::EventSetup &   
) [virtual]

Implements SimProducer.

Definition at line 125 of file HcalTB04Analysis.cc.

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

00125                                                                 {
00126 
00127   std::auto_ptr<PHcalTB04Info> product(new PHcalTB04Info);
00128   fillEvent(*product);
00129   e.put(product);
00130 }

void HcalTB04Analysis::qieAnalysis (  )  [private]

Definition at line 655 of file HcalTB04Analysis.cc.

References CaloHit::e(), edm::eq(), eqie, esimh, HcalQie::getCode(), HcalQie::getEnergy(), hcalHitCache, CaloHit::id(), idTower, LogDebug, myQie, and nTower.

Referenced by update().

00655                                    {
00656 
00657   int hittot = hcalHitCache.size();
00658   if (hittot<=0) hittot = 1;
00659   std::vector<CaloHit> hits(hittot);
00660   std::vector<int>     todo(nTower,0);
00661 
00662   LogDebug("HcalTBSim") << "HcalTB04Analysis::qieAnalysis: Size " 
00663                         << hits.size() << " " << todo.size() << " " 
00664                         << idTower.size() << " " << esimh.size() << " " 
00665                         << eqie.size();
00666   // Loop over all HCal hits
00667   for (unsigned int k1 = 0; k1 < hcalHitCache.size(); k1++) {
00668     CaloHit hit = hcalHitCache[k1];
00669     uint32_t id = hit.id();
00670     int    nhit = 0;
00671     double esim = hit.e();
00672     hits[nhit]  = hit;
00673     for (unsigned int k2 = k1+1; k2 < hcalHitCache.size(); k2++) {
00674       hit = hcalHitCache[k2];
00675       if (hit.id() == id) {
00676         nhit++;
00677         hits[nhit] = hit;
00678         esim += hit.e();
00679       }
00680     }
00681     k1 += nhit;
00682     nhit++;
00683     std::vector<int> cd = myQie->getCode(nhit,hits);
00684     double eq = myQie->getEnergy(cd);
00685     LogDebug("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << id 
00686                           << std::dec << " registers " << esim << " energy "
00687                           << "from " << nhit << " hits starting with hit # " 
00688                           << k1 << " energy with noise " << eq;
00689     for (int k2 = 0; k2 < nTower; k2++) {
00690       if (id == idTower[k2]) {
00691         todo[k2]  = 1;
00692         esimh[k2] = esim;
00693         eqie[k2]  = eq;
00694       }
00695     }
00696   }
00697   
00698   // Towers with no hit
00699   for (int k2 = 0; k2 < nTower; k2++) {
00700     if (todo[k2] == 0) {
00701       std::vector<int> cd = myQie->getCode(0,hits);
00702       double eq = myQie->getEnergy(cd);
00703       esimh[k2] = 0;
00704       eqie[k2]  = eq;
00705 #ifdef ddebug
00706       LogDebug("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex 
00707                             << idTower[k2] << std::dec << " registers " 
00708                             << esimh[k2] << " energy from hits and energy "
00709                             << "after QIE analysis " << eqie[k2];
00710 #endif
00711     }
00712   }
00713 }

double HcalTB04Analysis::scale ( int  det,
int  layer 
) [private]

Definition at line 999 of file HcalTB04Analysis.cc.

References HcalBarrel, scaleHB0, scaleHB16, scaleHE0, scaleHO, and tmp.

Referenced by fillBuffer().

00999                                                  {
01000 
01001   double tmp = 1.;
01002   if (det == static_cast<int>(HcalBarrel)) {
01003     if (layer == 1)       tmp = scaleHB0;
01004     else if (layer == 17) tmp = scaleHB16;
01005     else if (layer > 17)  tmp = scaleHO;
01006   } else {
01007     if (layer <= 2)       tmp = scaleHE0;
01008   }
01009   return tmp;
01010 }

double HcalTB04Analysis::timeOfFlight ( int  det,
int  layer,
double  eta 
) [private]

Definition at line 1012 of file HcalTB04Analysis.cc.

References beamOffset, funct::cos(), dist(), funct::exp(), HcalBarrel, LogDebug, funct::sin(), theta, and tmp.

Referenced by fillBuffer().

01012                                                                     {
01013 
01014   double theta = 2.0*atan(exp(-eta));
01015   double dist  = beamOffset;
01016   if (det == static_cast<int>(HcalBarrel)) {
01017     const double rLay[19] = {
01018       1836.0, 1902.0, 1962.0, 2022.0, 2082.0, 2142.0, 2202.0, 2262.0, 2322.0, 
01019       2382.0, 2448.0, 2514.0, 2580.0, 2646.0, 2712.0, 2776.0, 2862.5, 3847.0,
01020       4052.0};
01021     if (layer>0 && layer<=19) dist += rLay[layer-1]*mm/sin(theta);
01022   } else {
01023     const double zLay[19] = {
01024       4034.0, 4032.0, 4123.0, 4210.0, 4297.0, 4384.0, 4471.0, 4558.0, 4645.0, 
01025       4732.0, 4819.0, 4906.0, 4993.0, 5080.0, 5167.0, 5254.0, 5341.0, 5428.0,
01026       5515.0};
01027     if (layer>0 && layer<=19) dist += zLay[layer-1]*mm/cos(theta);
01028   }
01029 
01030   double tmp = dist/c_light/ns;
01031 #ifdef ddebug
01032   LogDebug("HcalTBSim") << "HcalTB04Analysis::timeOfFlight " << tmp 
01033                         << " for det/lay " << det  << " " << layer 
01034                         << " eta/theta " << eta << " " << theta/deg 
01035                         << " dist " << dist;
01036 #endif
01037   return tmp;
01038 }

int HcalTB04Analysis::unitID ( uint32_t  id  )  [private]

Definition at line 987 of file HcalTB04Analysis.cc.

References group, HcalTestNumbering::unpackHcalIndex(), and z.

Referenced by fillBuffer(), and init().

00987                                         {
00988 
00989   int det, z, group, ieta, iphi, lay;
00990   HcalTestNumbering::unpackHcalIndex(id,det,z,group,ieta,iphi,lay);
00991   group  = (det&15)<<20;
00992   group += ((lay-1)&31)<<15;
00993   group += (z&1)<<14;
00994   group += (ieta&127)<<7;
00995   group += (iphi&127);
00996   return group;
00997 }

void HcalTB04Analysis::update ( const EndOfEvent *   )  [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 362 of file HcalTB04Analysis.cc.

References count, ecalHitCache, fillBuffer(), finalAnalysis(), hcalHitCache, hcalOnly, LogDebug, qieAnalysis(), and xtalAnalysis().

00362                                                     {
00363 
00364   count++;
00365 
00366   //fill the buffer
00367   LogDebug("HcalTBSim") << "HcalTB04Analysis::Fill event " 
00368                         << (*evt)()->GetEventID();
00369   fillBuffer (evt);
00370 
00371   //QIE analysis
00372   LogDebug("HcalTBSim") << "HcalTB04Analysis::Do QIE analysis with " 
00373                         << hcalHitCache.size() << " hits";
00374   qieAnalysis();
00375 
00376   //Energy in Crystal Matrix
00377   if (!hcalOnly) {
00378     LogDebug("HcalTBSim") << "HcalTB04Analysis::Do Xtal analysis with " 
00379                           << ecalHitCache.size() << " hits";
00380     xtalAnalysis();
00381   }
00382   
00383   //Final Analysis
00384   LogDebug("HcalTBSim") << "HcalTB04Analysis::Final analysis";  
00385   finalAnalysis();
00386 
00387   int iEvt = (*evt)()->GetEventID();
00388   if (iEvt < 10) 
00389     edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
00390   else if ((iEvt < 100) && (iEvt%10 == 0)) 
00391     edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
00392   else if ((iEvt < 1000) && (iEvt%100 == 0)) 
00393     edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
00394   else if ((iEvt < 10000) && (iEvt%1000 == 0)) 
00395     edm::LogInfo("HcalTBSim") << "HcalTB04Analysis:: Event " << iEvt;
00396 }

void HcalTB04Analysis::update ( const G4Step *   )  [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 248 of file HcalTB04Analysis.cc.

References LogDebug, NULL, pvFound, pvMomentum, pvPosition, pvType, pvUVW, secEkin, secMomentum, secPartID, secTrackID, and shortLivedSecondaries.

00248                                                   {
00249 
00250   if (aStep != NULL) {
00251     //Get Step properties
00252     G4ThreeVector thePreStepPoint  = aStep->GetPreStepPoint()->GetPosition();
00253     G4ThreeVector thePostStepPoint;
00254 
00255     // Get Tracks properties
00256     G4Track*      aTrack   = aStep->GetTrack();
00257     int           trackID  = aTrack->GetTrackID();
00258     int           parentID = aTrack->GetParentID();
00259     G4ThreeVector position = aTrack->GetPosition();
00260     G4ThreeVector momentum = aTrack->GetMomentum();
00261     G4String      partType = aTrack->GetDefinition()->GetParticleType();
00262     G4String      partSubType = aTrack->GetDefinition()->GetParticleSubType();
00263     int    partPDGEncoding = aTrack->GetDefinition()->GetPDGEncoding();
00264 #ifdef ddebug
00265     bool   isPDGStable = aTrack->GetDefinition()->GetPDGStable();
00266 #endif
00267     double pDGlifetime = aTrack->GetDefinition()->GetPDGLifeTime();
00268     double gammaFactor = aStep->GetPreStepPoint()->GetGamma();
00269 
00270     if (!pvFound) { //search for v1
00271       double stepDeltaEnergy = aStep->GetDeltaEnergy ();
00272       double kinEnergy = aTrack->GetKineticEnergy ();
00273       
00274       // look for DeltaE > 10% kinEnergy of particle, or particle death - Ek=0
00275       if (trackID == 1 && parentID == 0 && 
00276           ((kinEnergy == 0.) || (fabs (stepDeltaEnergy / kinEnergy) > 0.1))) {
00277         int pvType = -1;
00278         if (kinEnergy == 0.) {
00279           pvType = 0;
00280         } else {
00281           if (fabs (stepDeltaEnergy / kinEnergy) > 0.1) pvType = 1;
00282         }
00283         pvFound    = true;
00284         pvPosition = position;
00285         pvMomentum = momentum;
00286         // Rotated coord.system:
00287         pvUVW      = (*beamline_RM)*(pvPosition);
00288 
00289         //Volume name requires some checks:
00290         G4String thePostPVname = "NoName";
00291         G4StepPoint * thePostPoint = aStep->GetPostStepPoint ();
00292         if (thePostPoint) {
00293           thePostStepPoint = thePostPoint->GetPosition();
00294           G4VPhysicalVolume * thePostPV = thePostPoint->GetPhysicalVolume ();
00295           if (thePostPV) thePostPVname = thePostPV->GetName ();
00296         }
00297 #ifdef ddebug
00298         LogDebug("HcalTBSim") << "HcalTB04Analysis:: V1 found at: " 
00299                               << thePostStepPoint << " G4VPhysicalVolume: " 
00300                               << thePostPVname;
00301 #endif      
00302         LogDebug("HcalTBSim") << "HcalTB04Analysis::fill_v1Pos: Primary Track "
00303                               << "momentum: " << pvMomentum << " psoition " 
00304                               << pvPosition << " u/v/w " << pvUVW;
00305       }
00306     } else { 
00307       // watch for secondaries originating @v1, including the surviving primary
00308       if ((trackID != 1 && parentID == 1 &&
00309            (aTrack->GetCurrentStepNumber () == 1) && 
00310            (thePreStepPoint == pvPosition)) || 
00311           (trackID == 1 && thePreStepPoint == pvPosition)) {
00312 #ifdef ddebug
00313         LogDebug("HcalTBSim") << "HcalTB04Analysis::A secondary...  PDG:" 
00314                               << partPDGEncoding << " TrackID:" << trackID
00315                               << " ParentID:" << parentID << " stable: "  
00316                               << isPDGStable << " Tau: " << pDGlifetime 
00317                               << " cTauGamma=" 
00318                               << c_light*pDGlifetime*gammaFactor*1000.
00319                               << "um" << " GammaFactor: " << gammaFactor;
00320 #endif      
00321         secTrackID.push_back(trackID);
00322         secPartID.push_back(partPDGEncoding);
00323         secMomentum.push_back(momentum);
00324         secEkin.push_back(aTrack->GetKineticEnergy());
00325 
00326         // Check for short-lived secondaries: cTauGamma<100um
00327         double ctaugamma_um = c_light * pDGlifetime * gammaFactor * 1000.;
00328         if ((ctaugamma_um>0.) && (ctaugamma_um<100.)) {//short-lived secondary
00329           shortLivedSecondaries.push_back(trackID);
00330       } else {//normal secondary - enter into the V1-calorimetric tree
00331         //          histos->fill_v1cSec (aTrack);
00332       }
00333       }
00334       // Also watch for tertiary particles coming from 
00335       // short-lived secondaries from V1
00336       if (aTrack->GetCurrentStepNumber() == 1) {
00337         if (shortLivedSecondaries.size() > 0) {
00338           int pid = parentID;
00339           std::vector<int>::iterator pos1= shortLivedSecondaries.begin();
00340           std::vector<int>::iterator pos2 = shortLivedSecondaries.end();
00341           std::vector<int>::iterator pos;
00342           for (pos = pos1; pos != pos2; pos++) {
00343             if (*pos == pid) {//ParentID is on the list of short-lived 
00344               // secondary 
00345 #ifdef ddebug
00346               LogDebug("HcalTBSim") << "HcalTB04Analysis:: A tertiary...  PDG:"
00347                                     << partPDGEncoding << " TrackID:" <<trackID
00348                                     << " ParentID:" << parentID << " stable: "
00349                                     << isPDGStable << " Tau: " << pDGlifetime
00350                                     << " cTauGamma=" 
00351                                     << c_light*pDGlifetime*gammaFactor*1000. 
00352                                     << "um GammaFactor: " << gammaFactor;
00353 #endif
00354             }
00355           }
00356         }
00357       }
00358     }
00359   }
00360 }

void HcalTB04Analysis::update ( const BeginOfEvent *   )  [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 240 of file HcalTB04Analysis.cc.

References clear(), and evNum.

00240                                                       {
00241  
00242   evNum = (*evt) ()->GetEventID ();
00243   clear();
00244   edm::LogInfo("HcalTBSim") << "HcalTB04Analysis: =====> Begin of event = "
00245                             << evNum;
00246 }

void HcalTB04Analysis::update ( const BeginOfRun *   )  [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfRun * >.

Definition at line 192 of file HcalTB04Analysis.cc.

References hcalOnly, names, ECalSD::setNumberingScheme(), and HCalSD::setNumberingScheme().

00192                                                     {
00193 
00194   int irun = (*run)()->GetRunID();
00195   edm::LogInfo("HcalTBSim") <<" =====> Begin of Run = " << irun;
00196  
00197   G4SDManager* sd     = G4SDManager::GetSDMpointerIfExist();
00198   if (sd != 0) {
00199     std::string  sdname = names[0];
00200     G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
00201     if (aSD==0) {
00202       edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
00203                                    << " with name " << sdname << " in this "
00204                                    << "Setup";
00205     } else {
00206       HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
00207       edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD "
00208                                 << "with name " << theCaloSD->GetName() 
00209                                 << " in this Setup";
00210       HcalNumberingScheme* org = new HcalTestNumberingScheme(false);
00211       theCaloSD->setNumberingScheme(org);
00212       edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a "
00213                                 << "new numbering scheme";
00214     }
00215     if (!hcalOnly) {
00216       sdname = names[1];
00217       aSD = sd->FindSensitiveDetector(sdname);
00218       if (aSD==0) {
00219         edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: No SD"
00220                                      << " with name " << sdname << " in this "
00221                                      << "Setup";
00222       } else {
00223         ECalSD* theCaloSD = dynamic_cast<ECalSD*>(aSD);
00224         edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Finds SD "
00225                                   << "with name " << theCaloSD->GetName() 
00226                                   << " in this Setup";
00227         EcalNumberingScheme* org = new HcalTB04XtalNumberingScheme();
00228         theCaloSD->setNumberingScheme(org);
00229         edm::LogInfo("HcalTBSim") << "HcalTB04Analysis::beginOfRun: set a "
00230                                   << "new numbering scheme";
00231       }
00232     }
00233   } else {
00234     edm::LogWarning("HcalTBSim") << "HcalTB04Analysis::beginOfRun: Could "
00235                                  << "not get SD Manager!";
00236   }
00237 
00238 }

void HcalTB04Analysis::xtalAnalysis (  )  [private]

Definition at line 715 of file HcalTB04Analysis.cc.

References ecalHitCache, ecalNoise, enois, edm::eq(), esime, Exception, idEcal, edm::Service< T >::isAvailable(), LogDebug, and nCrystal.

Referenced by update().

00715                                     {
00716 
00717   edm::Service<edm::RandomNumberGenerator> rng;
00718   if ( ! rng.isAvailable()) {
00719     throw cms::Exception("Configuration")
00720       << "HcalTB04Analysis requires the RandomNumberGeneratorService\n"
00721       << "which is not present in the configuration file. "
00722       << "You must add the service\n in the configuration file or "
00723       << "remove the modules that require it.";
00724   }
00725   CLHEP::RandGaussQ  randGauss(rng->getEngine());
00726 
00727   // Crystal Data
00728   std::vector<int> iok(nCrystal,0);
00729   LogDebug("HcalTBSim") << "HcalTB04Analysis::xtalAnalysis: Size " <<iok.size()
00730                         << " " << idEcal.size() << " " << esime.size() << " " 
00731                         << enois.size();
00732   for (unsigned int k1 = 0; k1 < ecalHitCache.size(); k1++) {
00733     uint32_t id   = ecalHitCache[k1].id();
00734     int      nhit = 0;
00735     double   esim = ecalHitCache[k1].e();
00736     for (unsigned int k2 = k1+1; k2 < ecalHitCache.size(); k2++) {
00737       if (ecalHitCache[k2].id() == id) {
00738         nhit++;
00739         esim += ecalHitCache[k2].e();
00740       }
00741     }
00742     k1 += nhit;
00743     nhit++;
00744     double eq = esim + randGauss.fire(0., ecalNoise);
00745 #ifdef ddebug
00746     LogDebug("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex << id 
00747                           << std::dec << " registers " << esim << " energy "
00748                           << "from " << nhit << " hits starting with hit # "
00749                           << k1 << " energy with noise " << eq;
00750 #endif
00751     for (int k2 = 0; k2 < nCrystal; k2++) {
00752       if (id == idEcal[k2]) {
00753         iok[k2]   = 1;
00754         esime[k2] = esim;
00755         enois[k2] = eq;
00756       }
00757     }
00758   }
00759     
00760   // Crystals with no hit
00761   for (int k2 = 0; k2 < nCrystal; k2++) {
00762     if (iok[k2] == 0) {
00763       esime[k2] = 0;
00764       enois[k2] = randGauss.fire(0., ecalNoise);
00765 #ifdef ddebug
00766       LogDebug("HcalTBSim") << "HcalTB04Analysis::  ID 0x" << std::hex 
00767                             << idEcal[k2] << std::dec << " registers " 
00768                             << esime[k2] << " energy from hits and energy from"
00769                             << " noise " << enois[k2];
00770 #endif
00771     }
00772   }
00773 }

Member Data Documentation

G4RotationMatrix* HcalTB04Analysis::beamline_RM [private]

Definition at line 98 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis().

double HcalTB04Analysis::beamOffset [private]

Definition at line 94 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), and timeOfFlight().

int HcalTB04Analysis::count [private]

Definition at line 101 of file HcalTB04Analysis.h.

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

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

Definition at line 109 of file HcalTB04Analysis.h.

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

double HcalTB04Analysis::ecalNoise [private]

Definition at line 94 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), and xtalAnalysis().

double HcalTB04Analysis::eecalq [private]

Definition at line 113 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

double HcalTB04Analysis::eecals [private]

Definition at line 113 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

double HcalTB04Analysis::ehcalq [private]

Definition at line 113 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

double HcalTB04Analysis::ehcals [private]

Definition at line 113 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

std::vector<double> HcalTB04Analysis::enois [private]

Definition at line 111 of file HcalTB04Analysis.h.

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

std::vector<double> HcalTB04Analysis::eqeta [private]

Definition at line 112 of file HcalTB04Analysis.h.

Referenced by fillEvent(), finalAnalysis(), and init().

std::vector<double> HcalTB04Analysis::eqie [private]

Definition at line 111 of file HcalTB04Analysis.h.

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

std::vector<double> HcalTB04Analysis::eqlay [private]

Definition at line 112 of file HcalTB04Analysis.h.

Referenced by fillEvent(), finalAnalysis(), and init().

std::vector<double> HcalTB04Analysis::eqphi [private]

Definition at line 112 of file HcalTB04Analysis.h.

Referenced by fillEvent(), finalAnalysis(), and init().

std::vector<double> HcalTB04Analysis::eseta [private]

Definition at line 112 of file HcalTB04Analysis.h.

Referenced by fillEvent(), finalAnalysis(), and init().

std::vector<double> HcalTB04Analysis::esime [private]

Definition at line 111 of file HcalTB04Analysis.h.

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

std::vector<double> HcalTB04Analysis::esimh [private]

Definition at line 111 of file HcalTB04Analysis.h.

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

std::vector<double> HcalTB04Analysis::eslay [private]

Definition at line 112 of file HcalTB04Analysis.h.

Referenced by fillEvent(), finalAnalysis(), and init().

std::vector<double> HcalTB04Analysis::esphi [private]

Definition at line 112 of file HcalTB04Analysis.h.

Referenced by fillEvent(), finalAnalysis(), and init().

double HcalTB04Analysis::etaInit [private]

Definition at line 108 of file HcalTB04Analysis.h.

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

double HcalTB04Analysis::etotq [private]

Definition at line 113 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

double HcalTB04Analysis::etots [private]

Definition at line 113 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and finalAnalysis().

int HcalTB04Analysis::evNum [private]

Definition at line 116 of file HcalTB04Analysis.h.

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

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

Definition at line 110 of file HcalTB04Analysis.h.

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

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

Definition at line 110 of file HcalTB04Analysis.h.

Referenced by clear(), and fillBuffer().

bool HcalTB04Analysis::hcalOnly [private]

Definition at line 92 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), init(), and update().

HcalTB04Histo* HcalTB04Analysis::histo [private]

Definition at line 89 of file HcalTB04Analysis.h.

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

int HcalTB04Analysis::iceta [private]

Definition at line 95 of file HcalTB04Analysis.h.

Referenced by finalAnalysis(), and HcalTB04Analysis().

int HcalTB04Analysis::icphi [private]

Definition at line 95 of file HcalTB04Analysis.h.

Referenced by finalAnalysis(), and HcalTB04Analysis().

std::vector<uint32_t> HcalTB04Analysis::idEcal [private]

Definition at line 104 of file HcalTB04Analysis.h.

Referenced by init(), and xtalAnalysis().

std::vector<int> HcalTB04Analysis::idHcal [private]

Definition at line 103 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and init().

std::vector<uint32_t> HcalTB04Analysis::idTower [private]

Definition at line 104 of file HcalTB04Analysis.h.

Referenced by finalAnalysis(), init(), and qieAnalysis().

std::vector<int> HcalTB04Analysis::idXtal [private]

Definition at line 103 of file HcalTB04Analysis.h.

Referenced by fillEvent(), and init().

int HcalTB04Analysis::mode [private]

Definition at line 93 of file HcalTB04Analysis.h.

Referenced by fillBuffer(), HcalTB04Analysis(), and init().

HcalQie* HcalTB04Analysis::myQie [private]

Definition at line 88 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), qieAnalysis(), and ~HcalTB04Analysis().

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

Definition at line 97 of file HcalTB04Analysis.h.

Referenced by fillBuffer(), HcalTB04Analysis(), and update().

int HcalTB04Analysis::nCrystal [private]

Definition at line 102 of file HcalTB04Analysis.h.

Referenced by clear(), finalAnalysis(), init(), and xtalAnalysis().

int HcalTB04Analysis::nPrimary [private]

Definition at line 107 of file HcalTB04Analysis.h.

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

int HcalTB04Analysis::nTower [private]

Definition at line 102 of file HcalTB04Analysis.h.

Referenced by clear(), finalAnalysis(), init(), and qieAnalysis().

int HcalTB04Analysis::particleType [private]

Definition at line 107 of file HcalTB04Analysis.h.

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

double HcalTB04Analysis::phiInit [private]

Definition at line 108 of file HcalTB04Analysis.h.

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

double HcalTB04Analysis::pInit [private]

Definition at line 108 of file HcalTB04Analysis.h.

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

bool HcalTB04Analysis::pvFound [private]

Definition at line 115 of file HcalTB04Analysis.h.

Referenced by clear(), and update().

G4ThreeVector HcalTB04Analysis::pvMomentum [private]

Definition at line 117 of file HcalTB04Analysis.h.

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

G4ThreeVector HcalTB04Analysis::pvPosition [private]

Definition at line 117 of file HcalTB04Analysis.h.

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

int HcalTB04Analysis::pvType [private]

Definition at line 116 of file HcalTB04Analysis.h.

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

G4ThreeVector HcalTB04Analysis::pvUVW [private]

Definition at line 117 of file HcalTB04Analysis.h.

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

double HcalTB04Analysis::scaleHB0 [private]

Definition at line 96 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), and scale().

double HcalTB04Analysis::scaleHB16 [private]

Definition at line 96 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), and scale().

double HcalTB04Analysis::scaleHE0 [private]

Definition at line 96 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), and scale().

double HcalTB04Analysis::scaleHO [private]

Definition at line 96 of file HcalTB04Analysis.h.

Referenced by HcalTB04Analysis(), and scale().

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

Definition at line 120 of file HcalTB04Analysis.h.

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

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

Definition at line 119 of file HcalTB04Analysis.h.

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

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

Definition at line 118 of file HcalTB04Analysis.h.

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

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

Definition at line 118 of file HcalTB04Analysis.h.

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

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

Definition at line 121 of file HcalTB04Analysis.h.

Referenced by clear(), and update().

int HcalTB04Analysis::type [private]

Definition at line 93 of file HcalTB04Analysis.h.

The documentation for this class was generated from the following files:
Generated on Tue Jun 9 18:23:58 2009 for CMSSW by  doxygen 1.5.4