162 G4double thisStepEnergyDeposit =
step->GetTotalEnergyDeposit() *
step->GetTrack()->GetWeight();
166 ATH_MSG_DEBUG(
"****** LArBarrelCalculator: Step energy,length "
167 << thisStepEnergyDeposit <<
" " << thisStepLength);
170 if(thisStepEnergyDeposit <= 0. || thisStepLength <= 0.)
179 G4StepPoint *thisStepPoint =
step->GetPreStepPoint();
180 G4StepPoint *thisStepBackPoint =
step->GetPostStepPoint();
181 G4ThreeVector startPoint = thisStepPoint->GetPosition();
182 G4ThreeVector endPoint = thisStepBackPoint->GetPosition();
186 << startPoint.x() <<
" " << startPoint.y() <<
" " << startPoint.z());
192 const G4NavigationHistory* g4navigation = thisStepPoint->GetTouchable()->GetHistory();
193 G4int ndep = g4navigation->GetDepth();
194 for (G4int ii=0;ii<=ndep;ii++) {
195 G4VPhysicalVolume* v1 = g4navigation->GetVolume(ii);
196 G4String vname = v1->GetName();
197 if ( vname.find(
"NegPhysical") != std::string::npos) zSide=-1;
204 const G4AffineTransform transformation =
205 thisStepPoint->GetTouchable()->GetHistory()->GetTopTransform();
207 G4ThreeVector startPointinLocal =
208 transformation.TransformPoint(startPoint);
209 G4ThreeVector endPointinLocal =
210 transformation.TransformPoint (endPoint);
214 << startPointinLocal.x() <<
" " << startPointinLocal.y() <<
" "
215 << startPointinLocal.z());
217 << endPointinLocal.x() <<
" " << endPointinLocal.y() <<
" "
218 << endPointinLocal.z());
221 G4double
energy =
step->GetTotalEnergyDeposit() *
step->GetTrack()->GetWeight();
224 const double EField = 10.;
225 const double wholeStepLengthCm =
step->GetStepLength() /
CLHEP::cm;
226 energy = (*m_birksLaw)(
energy, wholeStepLengthCm, EField);
237 G4double delta=1./((
double) nsub_step);
239 ATH_MSG_DEBUG(
" nsub_step,delta " << nsub_step <<
" " << delta);
246 for (G4int
i=0;
i<nsub_step;
i++) {
250 G4double fraction=(((G4double)
i)+0.5)*delta;
251 G4ThreeVector subinLocal=startPointinLocal*(1-fraction) + endPointinLocal*fraction;
252 G4double xloc = subinLocal.x();
253 G4double yloc = subinLocal.y();
254 G4double zloc = subinLocal.z();
255 G4double etaloc = subinLocal.pseudoRapidity();
256 G4double philoc = subinLocal.phi();
257 if(philoc<0.) philoc = philoc + 2.*
M_PI;
258 G4double radloc = sqrt( xloc*xloc + yloc*yloc );
261 << xloc <<
" " << yloc <<
" " << zloc);
274 if (currentCellData.
cellID == 0)
277 ATH_MSG_DEBUG(
"LArBarrelCalculator: Invalid hit CELLID == 0 ");
278 ATH_MSG_DEBUG(
"x,y,z local " << xloc <<
" " << yloc <<
" " << zloc);
279 ATH_MSG_DEBUG(
"radius " << radloc <<
" etaloc "<< etaloc <<
" philoc " << philoc);
283 G4int region = currentCellData.
region;
286 G4int sampling = currentCellData.
sampling;
291 if( sampling == 1 && region ==0 )
296 if( sampling == 1 && region ==1 )
308 ATH_MSG_DEBUG(
" region,side,sampling,eta,phi " << region <<
" " << zSide <<
" "
312 ATH_MSG_DEBUG(
" local coordinates " << currentCellData.
x0 <<
" " << currentCellData.
y0);
318 ATH_MSG_DEBUG(
" hit in absorber or electrode radius:" << radloc);
352 G4double Current_xt1,Current_xt2;
362 G4int nfold = currentCellData.
nfold;
363 G4double
x0=currentCellData.
x0;
364 G4double
y0=currentCellData.
y0;
365 if (x0<1500 || x0>1960 || y0>30 || y0<-30) {
369 G4double
rr = sqrt(x0*x0+y0*y0);
372 const CurrMap* mapOfCurrent(
nullptr);
376 x0 < m_accmap->GetXmax(nfold) &&
378 y0 < m_accmap->GetYmax(nfold) &&
379 (nfold != 13 || x0 < 1957.5) &&
380 (nfold>0 || x0 < 1504.6) ) {
384 ATH_MSG_DEBUG(
" Map for fold xmap,ymap " << nfold <<
" " << xmap <<
" " << ymap);
387 G4int sampMap = currentCellData.
sampMap;
393 ATH_MSG_WARNING(
" region,sampling,eta,fold " << region <<
" " << sampMap <<
" "
394 <<
etaMap <<
" " << nfold);
400 G4int nstraight = currentCellData.
nstraight;
401 if (nstraight%2==0)
n=22;
409 xmap = currentCellData.
xl;
412 if (nstraight==0) xmap = 0.5*(xmap+1.);
414 ATH_MSG_DEBUG(
" Map for straight xl,delec " << xmap <<
" " << ymap);
418 double current0,current1,current2;
420 mapOfCurrent->GetAll(xmap,ymap,&
gap,¤t0,¤t1,¤t2);
421 G4double gap2=std::fabs(currentCellData.
distElec)+std::fabs(currentCellData.
distAbs)
425 int ipm,ielec,ieta,iside;
428 ielec=currentCellData.
phiGap;
431 if(ielec < 0 ) ielec += 1024;
433 ieta=((
int) (etaloc*(1./0.2)));
436 if ((currentCellData.
distElec>0 && zSide==1)
437 || (currentCellData.
distElec<0 && zSide==-1) ) iside=1;
441 double hv=
m_hv[ipm][ielec][ieta][iside];
442 if (hv>1995.)
current=current0;
453 if (
gap>1
e-3 && gap2 >1
e-3) dgap=
gap/gap2-1;
464 if (!UseFold && std::fabs(currentCellData.
distElec)>2.1 && current0 < 0.1) {
467 <<
" str number " << currentCellData.
nstraight
468 <<
" sampling,eta " << sampling <<
" " <<
etaBin <<
" "
469 <<
" current/E " << current0);
480 if (fabs(etaloc-etaTrans) < 0.025) {
481 double x=std::fabs(zloc-radloc*sinh(etaTrans))/cosh(etaTrans);
482 double y=std::fabs(currentCellData.
distElec);
483 if ( x < m_etamap3->Xmax() ) {
486 Current = Current*resp;
497 double resp,xt0,xt1,xt2,deta;
501 &resp,&xt0,&xt1,&xt2);
503 ATH_MSG_DEBUG(
"hit in strip etaloc,etaBin,deta,delec,resp,xt0,xt1,xt2 "
504 << etaloc <<
" " <<
etaBin <<
" " << deta*1000 <<
" " << currentCellData.
distElec
505 <<
" " << resp <<
" " << xt0 <<
" " << xt1 <<
" " << xt2);
507 if (!InTrans) Current = Current*resp;
511 identifier_xt1.
clear();
512 identifier_xt1 << 4 << 1 << zSide << sampling << region << (
etaBin+1) <<
phiBin;
513 Current_xt1 = Current*xt1;
514 identifier_xt2.
clear();
515 identifier_xt2 << 4 << 1 << zSide << sampling << region << (
etaBin-1) <<
phiBin;
516 Current_xt2 = Current*xt2;
519 identifier_xt1.
clear();
520 identifier_xt1 << 4 << 1 << zSide << sampling << region << (
etaBin-1) <<
phiBin;
521 Current_xt1 = Current*xt1;
522 identifier_xt2.
clear();
523 identifier_xt2 << 4 << 1 << zSide << sampling << region << (
etaBin+1) <<
phiBin;
524 Current_xt2 = Current*xt2;
526 Current = Current*xt0;
529 else if (sampling==2 && !InTrans) {
534 ATH_MSG_DEBUG(
"hit in middle etaloc,etaBin,deta,delec,resp,xt0,xt1,xt2 "
535 << etaloc <<
" " <<
etaBin <<
" " << deta*1000 <<
" " << currentCellData.
distElec
538 Current = Current*resp;
549 for (
unsigned int i=0;
i<hdata.size();
i++) {
550 if (hdata[
i].
id==identifier2) {
551 hdata[
i].energy += Current;
552 hdata[
i].time +=
time*Current;
559 hdata.push_back(newdata);
563 for (
unsigned int i=0;
i<hdata.size();
i++) {
564 if (hdata[
i].
id==identifier_xt1) {
565 hdata[
i].energy += Current_xt1;
566 hdata[
i].time +=
time*Current_xt1;
572 LArHitData newdata = {identifier_xt1,
time*Current_xt1, Current_xt1};
573 hdata.push_back(newdata);
576 for (
unsigned int i=0;
i<hdata.size();
i++) {
577 if (hdata[
i].
id==identifier_xt2) {
578 hdata[
i].energy += Current_xt2;
579 hdata[
i].time +=
time*Current_xt2;
585 LArHitData newdata = {identifier_xt2,
time*Current_xt2, Current_xt2};
586 hdata.push_back(newdata);
594 ATH_MSG_DEBUG(
"Number of hits for this step " << m_nhits <<
" "
599 for (
unsigned int i=0;
i<hdata.size();
i++) {
600 if (std::fabs(hdata[
i].
energy)>1
e-6) hdata[
i].time=hdata[
i].time/hdata[
i].energy;
601 else hdata[
i].time=0.;
608 if (hdata.size()>0)
return true;