22 #include "G4ThreeVector.hh"
23 #include "G4StepPoint.hh"
26 #include "G4AffineTransform.hh"
27 #include "G4NavigationHistory.hh"
28 #include "G4VTouchable.hh"
29 #include "G4TouchableHistory.hh"
34 #include "GaudiKernel/ISvcLocator.h"
35 #include "GaudiKernel/Bootstrap.h"
54 ATH_MSG_DEBUG(
"LArBarrelCalculator: Beginning initialization ");
56 const double Birks_LAr_density = 1.396;
91 +
parameters->GetValue(
"LArEMBThickElecKapton"));
108 ATH_MSG_DEBUG(
" LArBarrelCalculator: start reading of current maps");
113 ATH_MSG_DEBUG(
" LArBarrelCalculator: end of reading current maps");
127 if(
m_IflCur)
ATH_MSG_DEBUG(
" LArBarrelCalculator: Deposited Energy dE/dX Corrected ==> CURRENT Option ON");
128 else ATH_MSG_DEBUG(
" LArBarrelCalculator: Crude Deposited Energy dE/dX NO CURRENT option");
130 else ATH_MSG_DEBUG(
" LArBarrelCalculator: Ignore effect of E field around eta=0.8 ");
141 return StatusCode::SUCCESS;
148 return StatusCode::SUCCESS;
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;
639 float defaultHvVal=2000.;
641 for (
int ipm=0;ipm<2;ipm++) {
642 for (
int ielec=0;ielec<1024;ielec++) {
643 for (
int ieta=0;ieta<7;ieta++) {
644 for (
int iside=0;iside<2;iside++) {
645 m_hv[ipm][ielec][ieta][iside] = defaultHvVal;
653 SmartIF<StoreGateSvc> pDetStore{Gaudi::svcLocator()->service(
"DetectorStore")};
655 ATH_MSG_WARNING(
"LArBarrelCalculator::InitHV() unable to get Detector Store! Use default HV values.");
661 if (pDetStore->retrieve(
manager)==StatusCode::SUCCESS) {
666 for (
unsigned int iSide=0;iSide<2;iSide++) {
668 for (
unsigned int iSector=0;iSector<2;iSector++) {
671 for (
unsigned int ielec=0;ielec<32;ielec++) {
673 unsigned jElec = ielec+32*iSector+64*
iPhi;
674 for (
unsigned int iGap=0;iGap<2;iGap++) {
675 double hv = hvdata.
voltage (electrode, iGap);
676 ATH_MSG_DEBUG(
" iSide,jElec,iEta,iGap,hv " << iSide <<
" " << jElec <<
" " <<
iEta <<
" " << iGap <<
" " << hv);
677 if (hv>-999.)
m_hv[iSide][jElec][
iEta][iGap] = hv;
705 if (std::fabs(curr1)>1
e-6) {
706 double x=curr0/curr1;
714 if (std::fabs(curr2)>1
e-6) {
715 double x=curr1/curr2;
723 if (std::fabs(curr2)>1
e-6) {