This routine will be called when the appropriate signal arrives.
90 CLHEP::HepRandomEngine* engine = G4Random::getTheEngine();
91 CLHEP::RandGaussQ randGauss(*engine);
94 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis::Fill event " 95 << (*evt)()->GetEventID();
98 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
103 int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
105 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis :: Hit Collection for " << sd
106 <<
" of ID " << HCHCid <<
" is obtained at " <<theHCHC;
109 nentries = theHCHC->entries();
110 if (nentries==0)
return;
119 XTALSid = G4SDManager::GetSDMpointer()->GetCollectionID(sd);
121 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis :: Hit Collection for " << sd
122 <<
" of ID " << XTALSid <<
" is obtained at " 124 xentries = theXTHC->entries();
125 if (xentries==0)
return;
128 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis :: There are " << nentries
129 <<
" HCal hits, and" << xentries <<
" xtal hits";
131 float ETot=0., xETot=0.;
132 int scintID=0, xtalID=0;
137 if (HCHCid >= 0 && theHCHC !=
nullptr) {
138 for ( ihit = 0; ihit < nentries; ihit++) {
143 float enEm = aHit->
getEM();
163 float TowerEne[8][18], TowerEneCF[8][18], LayerEne[19], EnRing[100];
164 for (
int iphi=0 ; iphi<8; iphi++) {
166 TowerEne[iphi][
jeta]=0.;
167 TowerEneCF[iphi][
jeta]=0.;
171 for (
int ilayer=0; ilayer<19; ilayer++) LayerEne[ilayer]=0.;
172 for (
int iring=0; iring<100; iring++) EnRing[iring]=0.;
181 if ( (layer!=17) && (layer!=18) ) {
187 TowerEne[(
int)phi][(
int)eta] += ETot;
189 TowerEneCF[(
int)phi][(
int)eta] += ETot*(1.+ 0.1*randGauss.fire() );
190 double dR=0.08727*
std::
sqrt( (eta-8.)*(eta-8.) + (phi-3.)*(phi-3.) );
191 EnRing[(
int)(dR/0.01)] += ETot;
194 LayerEne[layer] += ETot;
198 for (
int ilayer=0 ; ilayer<19 ; ilayer++) {
202 for (
int iring=0; iring<100; iring++)
205 for (
int iphi=0 ; iphi<8; iphi++) {
208 SEnergyN += TowerEneCF[iphi][
jeta] + 3.*randGauss.fire();
210 double Rand = 3.*randGauss.fire();
212 if ( (iphi>=0) && (iphi<7) ) {
218 if ( (iphi>=1) && (iphi<6) ) {
237 G4PrimaryParticle* thePrim=
nullptr;
238 G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
239 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis :: Event has " << nvertex
243 <<
"ERROR: no vertex";
245 for (
int i = 0 ;
i<nvertex;
i++) {
247 G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(
i);
248 if (avertex ==
nullptr) {
250 <<
"ERROR: pointer to vertex = 0";
252 int npart = avertex->GetNumberOfParticle();
253 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis::Vertex number :" <<
i 254 <<
" with " << npart <<
" particles";
255 if (thePrim==
nullptr) thePrim=avertex->GetPrimary(trackID);
259 double px=0.,py=0.,pz=0.;
261 if (thePrim !=
nullptr) {
262 px = thePrim->GetPx();
263 py = thePrim->GetPy();
264 pz = thePrim->GetPz();
268 <<
" ERROR: primary has p=0 ";
270 float costheta = pz/
pInit;
273 if (px != 0) phi = atan(py/px);
277 LogDebug(
"HcalTBSim") <<
"HcalTB02Analysis:: End Of Event ERROR: could" 278 <<
" not find primary ";
290 if (XTALSid >= 0 && theXTHC !=
nullptr) {
291 for (
int xihit = 0; xihit < xentries; xihit++) {
295 float xenEm = xaHit->
getEM();
296 float xenhad = xaHit->
getHadr();
302 float xCrysEne[7][7];
303 for (
int irow=0 ; irow<7; irow++) {
304 for (
int jcol=0 ; jcol<7; jcol++) {
305 xCrysEne[irow][jcol]=0.;
311 int xtalID = (*is).first;
312 xETot = (*is).second;
314 int irow = (
int)(xtalID/100.);
315 int jcol = (
int)(xtalID-100.*irow);
318 xCrysEne[irow][jcol] = xETot;
320 float dR=
std::sqrt( 0.01619*0.01619*(jcol-3)*(jcol-3) +
321 0.01606*0.01606*(irow-3)*(irow-3) );
324 if ( (irow>0) && (irow<6) ) {
325 if ( (jcol>0) && (jcol<6) ) {
328 xE5x5MatrixN += xCrysEne[irow][jcol] + 108.5*randGauss.fire();
330 if ( (irow>1) && (irow<5) ) {
331 if ( (jcol>1) && (jcol<5) ) {
333 xE3x3MatrixN += xCrysEne[irow][jcol] +108.5*randGauss.fire();
343 if ( theXTHC !=
nullptr ) {
353 int iEvt = (*evt)()->GetEventID();
355 std::cout <<
" Event " << iEvt << std::endl;
356 else if ((iEvt < 100) && (iEvt%10 == 0))
357 std::cout <<
" Event " << iEvt << std::endl;
358 else if ((iEvt < 1000) && (iEvt%100 == 0))
359 std::cout <<
" Event " << iEvt << std::endl;
360 else if ((iEvt < 10000) && (iEvt%1000 == 0))
361 std::cout <<
" Event " << iEvt << std::endl;
void fillTransProf(float u, float v)
Geom::Theta< T > theta() const
int getphiID(int sID) const
double getIncidentEnergy() const
std::map< int, float > energyInCrystals
int getlayerID(int sID) const
void fillAllTime(float v)
Tan< T >::type tan(const T &t)
std::map< int, float > primaries
std::map< int, float > energyInScints
int getTimeSliceID() const
std::vector< std::string > names
G4THitsCollection< CaloG4Hit > CaloG4HitCollection
void fillProfile(int ilayer, float value)
int getetaID(int sID) const
uint32_t getUnitID() const
Power< A, B >::type pow(const A &a, const B &b)