65 if(hitsCollection ==
nullptr)
67 G4ExceptionDescription msg;
69 G4Exception(
"Par03EventAction::GetHitsCollection()",
"MyCode0001",
73 auto analysisManager = G4AnalysisManager::Instance();
83 auto primaryVertex = G4EventManager::GetEventManager()
84 ->GetConstCurrentEvent()
86 auto primaryParticle = primaryVertex->GetPrimary(0);
87 G4double primaryEnergy = primaryParticle->GetTotalEnergy();
90 auto primaryDirection = primaryParticle->GetMomentumDirection();
91 auto primaryEntrance = primaryVertex->GetPosition() -
92 primaryVertex->GetPosition().z() * primaryDirection;
93 G4double cosDirection = std::cos(primaryDirection.theta());
94 G4double sinDirection = std::sin(primaryDirection.theta());
99 G4double totalEnergy = 0;
103 G4double tDistance = 0., rDistance = 0.;
104 G4double tFirstMoment = 0., tSecondMoment = 0.;
105 G4double rFirstMoment = 0., rSecondMoment = 0.;
106 for(
size_t iHit = 0; iHit < hitsCollection->entries(); iHit++)
108 hit =
static_cast<Par03Hit*
>(hitsCollection->GetHit(iHit));
115 totalEnergy += hitEn;
118 (hitRho *
fCellSizeRho - primaryEntrance.perp()) * sinDirection;
121 (hitRho *
fCellSizeRho - primaryEntrance.perp()) * cosDirection;
122 tFirstMoment += hitEn * tDistance;
123 rFirstMoment += hitEn * rDistance;
124 analysisManager->FillH1(4, tDistance, hitEn);
125 analysisManager->FillH1(5, rDistance, hitEn);
126 analysisManager->FillH1(10, hitType);
129 tFirstMoment /= totalEnergy;
130 rFirstMoment /= totalEnergy;
131 analysisManager->FillH1(0, primaryEnergy / GeV);
132 analysisManager->FillH1(1, totalEnergy / GeV);
133 analysisManager->FillH1(2, totalEnergy / primaryEnergy);
134 analysisManager->FillH1(3,
fTimer.GetRealElapsed());
135 analysisManager->FillH1(6, tFirstMoment);
136 analysisManager->FillH1(7, rFirstMoment);
139 for(
size_t iHit = 0; iHit < hitsCollection->entries(); iHit++)
141 hit =
static_cast<Par03Hit*
>(hitsCollection->GetHit(iHit));
147 tDistance = hitZ *
fCellSizeZ * cosDirection +
148 (hitRho *
fCellSizeRho - primaryEntrance.r()) * sinDirection;
149 rDistance = hitZ *
fCellSizeZ * (-sinDirection) +
150 (hitRho *
fCellSizeRho - primaryEntrance.r()) * cosDirection;
151 tSecondMoment += hitEn * std::pow(tDistance - tFirstMoment, 2);
152 rSecondMoment += hitEn * std::pow(rDistance - rFirstMoment, 2);
155 tSecondMoment /= totalEnergy;
156 rSecondMoment /= totalEnergy;
157 analysisManager->FillH1(8, tSecondMoment);
158 analysisManager->FillH1(9, rSecondMoment);