43 : base_class(
name, pSvcLocator)
72 for (G4int idat = 0; idat <
m_Nbrt1 ; idat++) {
107 return StatusCode::SUCCESS;
131 return StatusCode::SUCCESS;
153 const double &yhit,
const int &PhiCell,
int &Num_Straight,
154 const int &Num_Coude,
double &xl)
const
167 const G4double radfold = sqrt(Xc[0]*Xc[0]+Xc[1]*Xc[1]);
168 const G4double radhit = sqrt(xhit*xhit+yhit*yhit);
172 if (Num_Coude == Num_Straight && radhit <radfold) {
173 if (Num_Straight>0) { Num_Straight = Num_Straight-1; }
176 if (Num_Coude == (Num_Straight+1) && radhit > radfold) {
177 if (Num_Straight<12) { Num_Straight = Num_Straight+1; }
186 double dx = xhit - Xm[0];
187 double dy = yhit - Xm[1];
191 const double hit =
dx*
u[0] +
dy*
u[1];
198 if(std::fabs(hit) < Half_Elec) {
209 if (Num_Coude==Num_Straight) { xl=-1.; }
227 const double &yhit,
const int &PhiCell,
const int &Num_Straight,
228 const int &Num_Coude)
const
244 double dx = xhit - Xm[0];
double dy = yhit - Xm[1];
248 const double hit =
dx*
u[0] +
dy*
u[1];
294 G4int& iregion, G4int& isampling, G4int& ieta,
295 G4int& isamp2, G4int& ieta2)
const
299 G4double Eta_max,R_max_acc,Z_max_acc,R_min_acc,Z_max_lowr,dzdr;
300 G4double zmax1,zmax2,zmax3,zmax4,zmax5,zmax6,zmax7,rmax1,rmax2,rmax3,rmax4;
304 static const Geo
g = []() {
309 g.Eta_max =
parameters->GetValue(
"LArEMBMaxEtaAcceptance");
311 g.R_min_acc=
parameters->GetValue(
"LArEMBRadiusInnerAccordion");
313 g.R_max_acc =
parameters->GetValue(
"LArEMBFiducialRmax");
316 g.Z_max_acc =
parameters->GetValue(
"LArEMBfiducialMothZmax");
318 const G4double R_min_highz=1548.;
321 const G4double deltaz=7.;
322 g.zmax1=
g.Z_max_acc-deltaz;
323 g.zmax2=
g.Z_max_acc-2.*deltaz;
324 g.zmax3=
g.Z_max_acc-3.*deltaz;
325 g.zmax4=
g.Z_max_acc-4.*deltaz;
326 g.zmax5=
g.Z_max_acc-5.*deltaz;
327 g.zmax6=
g.Z_max_acc-6.*deltaz;
328 g.zmax7=
g.Z_max_acc-7.*deltaz;
329 g.rmax1=
g.R_max_acc-deltaz;
330 g.rmax2=
g.R_max_acc-2.*deltaz;
331 g.rmax3=
g.R_max_acc-3.*deltaz;
332 g.rmax4=
g.R_max_acc-4.*deltaz;
335 g.Z_max_lowr = sinh(
g.Eta_max)*
g.R_min_acc;
337 g.dzdr = (
g.Z_max_acc-
g.Z_max_lowr)/(R_min_highz-
g.R_min_acc);
348 const G4double Dr_s12=1.1;
349 const G4double Eta_max_s1=1.4;
350 const G4double Eta_max_s3=1.325;
351 const G4double deta=0.025;
357 const G4double aeta=std::fabs(eta);
368 if (aeta<Eta_max_s1) {
371 istrip=(
int) (aeta/deta*8.);
372 if (istrip<0 || istrip >=448) {
373 ATH_MSG_ERROR(
" Problem aeta,istrip " << aeta <<
" " << istrip);
379 imid = (
int) (aeta/deta);
380 if (imid <0 || imid >=56) {
392 if (ieta==0) iactive=0;
398 else if (
radius < (r12+Dr_s12)) {
409 if (aeta<Eta_max_s3) {
424 else if (
z>
g.zmax2) {
428 else if (
z>
g.zmax3) {
432 else if (
z>
g.zmax4) {
436 else if (aeta<1.3 && z>
g.zmax5) {
440 else if (aeta<1.3 && z>
g.zmax6) {
444 else if (
radius>
g.rmax4 && z<g.zmax5 && aeta>1.2) {
485 if (isampling==3 &&
z<
g.zmax4 && (
radius<
g.rmax4 || aeta<1.2) ) {
486 const double etastr = (imid%2==0) ? 0.025*imid : 0.025*(imid+1);
487 const double delta=
radius*(sinh(etastr)-sinh(aeta))/cosh(etastr);
489 if (aeta<0.475) { deltastr=1.5;}
490 else if (aeta<0.80) { deltastr=2.75;}
491 else if (aeta<0.85) { deltastr=1.5;}
492 else if (aeta<1.1) { deltastr=2.75;}
493 else { deltastr=3.25;}
495 if (std::fabs(delta)<deltastr) {
507 else if (
z>
g.zmax2 && aeta<1.375) {
518 if (aeta>=Eta_max_s1 && aeta<
g.Eta_max) {
520 r23=
g.Z_max_acc/sinh(aeta);
522 const double zmax =
g.Z_max_lowr +
g.dzdr*(
radius-
g.R_min_acc);
527 ieta=
int((aeta-Eta_max_s1)/deta);
528 if (
z>
zmax) { iactive=0; }
530 else if (
radius < (r12+Dr_s12)) {
532 ieta=
int((aeta-Eta_max_s1)/deta);
538 if (
z>
zmax) { iactive=0; }
550 if (aeta>
g.Eta_max) {
567 if (
z>
g.Z_max_acc ||
radius>
g.R_max_acc || radius<g.R_min_acc || aeta >
g.Eta_max) iactive=0;
580 const double &xPosition,
581 const double &yPosition,
582 const double &zPosition,
583 const double &aRadius,
586 const bool MapDetail)
const
589 currentCellData.cellID = 0;
599 currentCellData.nstraight=0;
601 if (
m_rc[
i] > aRadius) {
break; }
602 currentCellData.nstraight++;
604 if (currentCellData.nstraight <0 || currentCellData.nstraight >=
m_Nbrt) {
605 ATH_MSG_ERROR(
"Invalid straight number " << currentCellData.nstraight <<
" " << aRadius);
610 currentCellData.nfold=currentCellData.nstraight;
611 if (std::fabs(aRadius-
m_rc[currentCellData.nfold]) > std::fabs(aRadius-
m_rc[currentCellData.nfold+1]) ) {
612 currentCellData.nfold +=1;
614 if (currentCellData.nfold <0 || currentCellData.nfold >=
m_Nbrt1) {
615 ATH_MSG_ERROR(
"Invalid fold number " << currentCellData.nfold);
621 ATH_MSG_VERBOSE(
" BarrelGeometry: radius,eta,phi " << aRadius <<
" " << anEta <<
" ");
622 ATH_MSG_VERBOSE(
" Straight/Fold numbers " << currentCellData.nstraight <<
" " << currentCellData.nfold);
626 G4int ireg,isamp,ieta,isamp2,ieta2;
627 currentCellData.cellID = this->
SampSeg(anEta,aRadius,zPosition,ireg,isamp,ieta,isamp2,ieta2);
629 currentCellData.etaBin = ieta;
630 currentCellData.sampling = isamp;
631 currentCellData.region = ireg;
632 currentCellData.etaMap = ieta2;
633 currentCellData.sampMap = isamp2;
636 int phicell = this->
PhiGap(aRadius,xPosition,yPosition);
637 if (phicell<0) phicell=0;
643 currentCellData.cellID=0;
652 int sampling_phi_nGaps=4;
653 if (currentCellData.region==0 && currentCellData.sampling==1) { sampling_phi_nGaps=16; }
655 if (currentCellData.cellID==0) {
656 currentCellData.phiBin = (G4int) ( phicell/sampling_phi_nGaps );
657 currentCellData.distElec=9999.;
663 G4int nstr = currentCellData.nstraight;
664 const G4double distElec = this->
Distance_Ele(xPosition,yPosition,phicell,nstr,currentCellData.nfold,xl);
671 if (std::fabs(distElec) > 2.5) {
673 double dElecMin=distElec;
675 int phicellmin=phicell;
676 for (
int ii=-2;ii<3;ii++) {
677 if (ii==0) {
continue; }
678 int phicellnew = phicell+ii;
681 if (phicellnew < 0) { phicellnew +=
m_NCellTot; }
684 int nstr2=currentCellData.nstraight;
685 double dElec =
Distance_Ele(xPosition,yPosition,phicellnew,nstr2,currentCellData.nfold,xln);
686 if (std::fabs(dElec)<std::fabs(dElecMin)) {
687 phicellmin=phicellnew;
693 currentCellData.phiGap = phicellmin;
694 currentCellData.distElec = dElecMin;
695 currentCellData.xl = xlmin;
696 currentCellData.nstraight = nstr;
699 currentCellData.phiGap=phicell;
700 currentCellData.distElec=distElec;
701 currentCellData.xl=xl;
702 currentCellData.nstraight=nstr;
706 ATH_MSG_VERBOSE(
" final phiGap,distElec,xl " << currentCellData.phiGap <<
" " << currentCellData.distElec <<
" "
707 << currentCellData.xl);
713 if (currentCellData.distElec<0) nabs=currentCellData.phiGap;
714 else nabs=currentCellData.phiGap+1;
716 currentCellData.distAbs =
Distance_Abs(xPosition,yPosition,nabs,currentCellData.nstraight,currentCellData.nfold);
718 ATH_MSG_VERBOSE(
" nabs,distAbs " << nabs <<
" " << currentCellData.distAbs);
723 if ((currentCellData.distAbs>0. && currentCellData.distElec>0) ||
724 (currentCellData.distAbs<0. && currentCellData.distElec<0) ) {
727 if (std::fabs(currentCellData.distElec)>std::fabs(currentCellData.distAbs)) {
728 if (currentCellData.distAbs>0) { currentCellData.phiGap += 1; }
729 if (currentCellData.distAbs<0) { currentCellData.phiGap -= 1; }
731 if (currentCellData.phiGap <0) { currentCellData.phiGap=0; }
735 if (currentCellData.phiGap < 0) { currentCellData.phiGap +=
m_NCellTot; }
738 currentCellData.distElec =
Distance_Ele(xPosition,yPosition,currentCellData.phiGap,currentCellData.nstraight,currentCellData.nfold,currentCellData.xl);
743 currentCellData.phiBin = (G4int) ( currentCellData.phiGap/sampling_phi_nGaps );
753 currentCellData.x0 = dx1 +
m_xc[currentCellData.nfold];
754 currentCellData.y0 = dy1 +
m_yc[currentCellData.nfold];
755 if (
m_parity==1) { currentCellData.y0 = -1*currentCellData.y0; }
765 const G4double
dl=0.001;
766 const G4double inv_dl = 1. /
dl;
767 G4double cenx[15],ceny[15];
769 G4double sum1[5000],sumx[5000];
778 const G4double inv_rint = 1. / rint;
779 const G4double
dt=
dl * inv_rint;
780 const G4double inv_dt = 1. /
dt;
786 for (G4int
i=0;
i<15;
i++) {
791 for (G4int
i=0;
i<15;
i++) {
828 const G4int nstep=
int((phi1-
phi0)*inv_dt)+1;
829 for (
int ii=0;ii<nstep;ii++) {
831 const G4double
phi=
phi0+
dt*((G4double)ii);
832 const G4double
x=cenx[
i]+rint*
cos(phi);
833 const G4double
y=ceny[
i]+rint*
sin(phi);
836 const G4double phid=
atan(
y/
x);
846 const G4double
dx=cenx[
i+1]-cenx[
i];
847 const G4double
dy=ceny[
i+1]-ceny[
i];
848 const G4double along=std::sqrt(
dx*
dx+
dy*
dy-4.*rint*rint);
849 const G4double
x0=0.5*(cenx[
i+1]+cenx[
i]);
850 const G4double
y0=0.5*(ceny[
i+1]+ceny[
i]);
852 const G4double
x1=
x0-0.5*along*
cos(phi);
853 const G4double
y1=
y0-0.5*along*
sin(phi);
855 const int nstep=
int(along*inv_dl)+1;
856 for (
int ii=0;ii<nstep;ii++) {
858 const G4double
x=
x1+
dl*((G4double)ii)*
cos(phi);
859 const G4double
y=
y1+
dl*((G4double)ii)*
sin(phi);
862 const G4double phid=
atan(
y/
x);
905 const G4double phi_hit=atan2(yhit,xhit);
906 G4double dphi=phi_hit-
phi0;
908 if (dphi<0) dphi=dphi+2*
M_PI;
910 dphi=dphi/(2*
M_PI)*1024;
911 const G4int ngap=((
int) dphi);
928 const G4NavigationHistory* g4navigation = a_step->GetPreStepPoint()->GetTouchable()->GetHistory();
929 const G4int ndep = g4navigation->GetDepth();
932 G4int indECAM = -999;
935 for (G4int ii=0;ii<=ndep;ii++) {
936 const G4String& vname = g4navigation->GetVolume(ii)->GetName();
938 if ( indECAM<0 && vname ==
m_ecamName ) indECAM=ii;
939 if ( !inSTAC && vname.find(
"STAC") !=std::string::npos) inSTAC=
true;
940 if ( vname.find(
"NegPhysical") != std::string::npos) zside=-1;
945 ATH_MSG_ERROR(
"LArBarrel::Geometry::CalculateIdentifier ECAM volume not found in hierarchy");
1008 const G4StepPoint *thisStepPoint = a_step->GetPreStepPoint();
1009 const G4StepPoint *thisStepBackPoint = a_step->GetPostStepPoint();
1010 const G4ThreeVector startPoint = thisStepPoint->GetPosition();
1011 const G4ThreeVector endPoint = thisStepBackPoint->GetPosition();
1012 const G4ThreeVector
p = (thisStepPoint->GetPosition() + thisStepBackPoint->GetPosition()) * 0.5;
1015 ATH_MSG_VERBOSE(
"Position of the step in the ATLAS frame (x,y,z) --> " <<
p.x() <<
" " <<
p.y() <<
" " <<
p.z());
1016 ATH_MSG_VERBOSE(
"Eta and Phi in the ATLAS frame --> " <<
p.eta() <<
" " <<
p.phi());
1021 const G4NavigationHistory* g4navigation = thisStepPoint->GetTouchable()->GetHistory();
1022 const G4AffineTransform transformation = g4navigation->GetTransform(indECAM);
1023 const G4ThreeVector startPointinLocal = transformation.TransformPoint(startPoint);
1024 const G4ThreeVector endPointinLocal = transformation.TransformPoint (endPoint);
1025 const G4ThreeVector midinLocal = (startPointinLocal+endPointinLocal)*0.5;
1028 ATH_MSG_VERBOSE(
"Position of the step in the LOCAL frame (x,y,z) --> " << midinLocal.x() <<
" " << midinLocal.y() <<
" " << midinLocal.z());
1029 ATH_MSG_VERBOSE(
"Eta and Phi of the step in LOCAL frame --> " << midinLocal.eta() <<
" " << midinLocal.phi());
1034 const G4double xZpos = midinLocal.x();
1035 const G4double yZpos = midinLocal.y();
1036 const G4double zZpos = midinLocal.z();
1037 const G4double etaZpos = midinLocal.pseudoRapidity();
1038 const G4double phiZpos = (midinLocal.phi()<0.) ? midinLocal.phi() + 2.*
M_PI : midinLocal.phi();
1039 const G4double radius2Zpos = xZpos*xZpos + yZpos*yZpos;
1040 const G4double radiusZpos = sqrt(radius2Zpos);
1044 currentCellData.
zSide = 1;
1047 currentCellData.
zSide = zside;
1064 const bool MapDetail(
false);
1065 this->
findCell( currentCellData, xZpos, yZpos, zZpos, radiusZpos, etaZpos, phiZpos, MapDetail );
1069 if( currentCellData.
zSide == -1 )
1071 if( currentCellData.
sampling == 1 && currentCellData.
region ==0 )
1074 if(currentCellData.
phiBin < 0 ) currentCellData.
phiBin += 64;
1076 if( currentCellData.
sampling == 1 && currentCellData.
region ==1 )
1079 if(currentCellData.
phiBin < 0 ) currentCellData.
phiBin += 256;
1081 if( currentCellData.
sampling >= 2 )
1084 if(currentCellData.
phiBin < 0 ) currentCellData.
phiBin += 256;
1092 currentCellData.
etaBin=1;
1097 << currentCellData.
zSide
1099 << currentCellData.
region
1100 << currentCellData.
etaBin
1101 << currentCellData.
phiBin;
1111 ATH_MSG_VERBOSE(
"And also etafirst ----> " << thisStepPoint->GetPosition().pseudoRapidity());
1128 const G4int numDeadPhiBins = 64;
1129 double abs_eta = std::fabs(etaZpos);
1130 const double DM1EtaWidth = 0.1 ;
1131 const double DM1PhiWidth = 2.*
M_PI / numDeadPhiBins ;
1132 currentCellData.
etaBin = (G4int) ( abs_eta * (1./DM1EtaWidth) ) ;
1133 currentCellData.
phiBin = (G4int) (phiZpos/ DM1PhiWidth );
1136 if (currentCellData.
phiBin==numDeadPhiBins) currentCellData.
phiBin=currentCellData.
phiBin-1;
1140 if ( currentCellData.
zSide == -1 ) {
1142 if (currentCellData.
phiBin < 0 ) currentCellData.
phiBin +=64 ;
1149 if (currentCellData.
etaBin > 14) currentCellData.
etaBin=14;
1152 ATH_MSG_VERBOSE(
"This hit is in the ECAM volume in front of the accordion (DEAD MATERIAL) !!!!! ");
1158 if (abs_eta < 1.0 ) {
1161 ATH_MSG_VERBOSE(
"This hit is in the ECAM volume behind accordion (DEAD MATERIAL 0) !!!!! ");
1163 }
else if ( abs_eta >= 1.0 && abs_eta < 1.5) {
1167 ATH_MSG_VERBOSE(
"This hit is in the ECAM volume behind accordion (DEAD MATERIAL 2) !!!!! ");
1170 ATH_MSG_ERROR(
" LArBarrelGeometry: hit behind accordion at eta>1.5 !!! ");
1172 currentCellData.
etaBin = 4;
1178 const G4int phisave=currentCellData.
phiBin;
1179 const G4bool MapDetail(
false);
1180 this->
findCell( currentCellData, xZpos, yZpos, zZpos, radiusZpos, etaZpos, phiZpos, MapDetail );
1181 sampling = currentCellData.
sampling;
1182 currentCellData.
etaBin=0;
1183 currentCellData.
phiBin=phisave;
1186 if (abs_eta >1.0 && abs_eta < 1.5) {
1190 }
else if (abs_eta < 1.6) {
1193 currentCellData.
etaBin=0;
1198 currentCellData.
etaBin=0;
1202 const G4double thisStepEnergyDeposit = a_step->GetTotalEnergyDeposit() * a_step->GetTrack()->GetWeight();
1203 std::ostringstream dmLog;
1204 dmLog <<
"LArBarrelGeometry: cannot find region for DM hit..." << std::endl;
1205 dmLog <<
"LArBarrelGeometry: cannot find region for DM hit..." << std::endl;
1206 dmLog <<
"m_zMinBarrel: " <<
m_zMinBarrel << std::endl;
1207 dmLog <<
"m_zMaxBarrel: " <<
m_zMaxBarrel << std::endl;
1210 dmLog <<
"r,z,eta,phi " << radiusZpos <<
" " << zZpos <<
" " << etaZpos <<
" " << phiZpos << std::endl;
1211 dmLog <<
"x,y,z (Atlas) " <<
p.x() <<
" " <<
p.y() <<
" " <<
p.z() << std::endl;
1212 dmLog <<
" inSTAC " << inSTAC << std::endl;
1213 dmLog <<
" eDeposited " << thisStepEnergyDeposit << std::endl;
1214 const G4VPhysicalVolume* vol = thisStepPoint->GetPhysicalVolume();
1215 const G4String volName = vol->GetName();
1216 dmLog <<
" volName " << volName << std::endl;
1217 const G4int ndep = g4navigation->GetDepth();
1218 for (G4int ii=0;ii<=ndep;ii++) {
1219 const G4VPhysicalVolume* v1 = g4navigation->GetVolume(ii);
1220 const G4String vname = v1->GetName();
1221 dmLog <<
"vname " << vname << std::endl;
1233 G4bool MapDetail=
false;
1234 this->
findCell( currentCellData, xZpos, yZpos, zZpos, radiusZpos, etaZpos, phiZpos, MapDetail );
1240 currentCellData.
etaBin=1;
1246 << currentCellData.
zSide
1248 << currentCellData.
region
1249 << currentCellData.
etaBin
1250 << currentCellData.
phiBin;
1256 << currentCellData.
zSide * 4
1260 << currentCellData.
etaBin
1261 << currentCellData.
phiBin;
1264 ATH_MSG_VERBOSE(
"Here the identifier for the barrel DEAD materials ---->");
1289 if (type <1 || type > 2)
return false;
1291 if (sampling<1 || sampling>2)
return false;
1293 if (region!=3 && region !=4)
return false;
1294 if (phi<0 || phi>63)
return false;
1296 if (eta<0 || eta>14)
return false;
1299 if (eta !=0)
return false;
1303 if (region !=0 && region !=2)
return false;
1304 if (phi<0 || phi>63)
return false;
1306 if (eta<0 || eta>9)
return false;
1309 if (eta<0 || eta>4)
return false;
1314 if (sampling<1 || sampling >3)
return false;
1315 if (region !=0)
return false;
1316 if (eta!=0)
return false;
1317 if (phi<0 || phi>63)
return false;
1325 if (sampling<0 || sampling >3)
return false;
1327 if (region!=0)
return false;
1328 if (eta<0 || eta>60)
return false;
1329 if (phi<0 || phi>63)
return false;
1332 if (region<0 || region >1)
return false;
1334 if (eta<1 || eta>447)
return false;
1335 if (phi<0 || phi>63)
return false;
1338 if (eta<0 || eta>2)
return false;
1339 if (phi<0 || phi>255)
return false;
1343 if (region<0 || region >1)
return false;
1345 if (eta<0 || eta>55)
return false;
1346 if (phi<0 || phi>255)
return false;
1349 if (eta!=0)
return false;
1350 if (phi<0 || phi>255)
return false;
1354 if (region !=0)
return false;
1355 if (eta<0 || eta>26)
return false;
1356 if (phi<0 || phi>255)
return false;