Skip to content
Snippets Groups Projects
Commit 924ec09d authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

* fixing ambiguous isnan call

parent 9b5f6cc0
No related branches found
No related tags found
No related merge requests found
Pipeline #197776 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
NPSimulation/Detectors/TACTIC/TACTICScorer.cc
+ 83
− 65
View file @ 924ec09d
//#define USE_Garfield //only use if compiled with Garfield //#define USE_Garfield //only use if compiled with Garfield
#include "TACTICScorer.hh" #include "TACTICScorer.hh"
#include "G4UnitsTable.hh"
#include "G4RunManager.hh" #include "G4RunManager.hh"
#include "G4UnitsTable.hh"
#include "TACTIC.hh" #include "TACTIC.hh"
#ifdef USE_Garfield #ifdef USE_Garfield
...@@ -13,13 +13,14 @@ double excess; ...@@ -13,13 +13,14 @@ double excess;
using namespace TACTICScorer; using namespace TACTICScorer;
Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int NumberOfSegments, G4int depth, G4double p0, G4double p1, G4double p2, G4double p3,string Shape) //what do level and depth do? Gas_Scorer::Gas_Scorer(G4String name, G4int Level, G4double ScorerLength, G4int NumberOfSegments, G4int depth,
:G4VPrimitiveScorer(name, depth),HCID(-1){ G4double p0, G4double p1, G4double p2, G4double p3, string Shape) // what do level and depth do?
: G4VPrimitiveScorer(name, depth), HCID(-1) {
m_ScorerLength = ScorerLength; m_ScorerLength = ScorerLength;
m_NumberOfSegments = NumberOfSegments; m_NumberOfSegments = NumberOfSegments;
m_SegmentLength = m_ScorerLength / m_NumberOfSegments; m_SegmentLength = m_ScorerLength / m_NumberOfSegments;
m_Level = Level; m_Level = Level;
m_Position = G4ThreeVector(-1000,-1000,-1000); m_Position = G4ThreeVector(-1000, -1000, -1000);
m_SegmentNumber = -1; m_SegmentNumber = -1;
m_Index = -1; m_Index = -1;
m_p0 = p0; m_p0 = p0;
...@@ -29,122 +30,139 @@ Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int Num ...@@ -29,122 +30,139 @@ Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int Num
m_Shape = Shape; m_Shape = Shape;
} }
Gas_Scorer::~Gas_Scorer(){} Gas_Scorer::~Gas_Scorer() {}
G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
G4double* Infos = new G4double[16]; G4double* Infos = new G4double[16];
m_Position = aStep->GetPreStepPoint()->GetPosition(); m_Position = aStep->GetPreStepPoint()->GetPosition();
Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID(); Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
Infos[1] = aStep->GetTrack()->GetTrackID(); Infos[1] = aStep->GetTrack()->GetTrackID();
Infos[2] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();; Infos[2] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();
;
Infos[3] = aStep->GetPreStepPoint()->GetGlobalTime(); Infos[3] = aStep->GetPreStepPoint()->GetGlobalTime();
Infos[4] = aStep->GetPreStepPoint()->GetKineticEnergy(); Infos[4] = aStep->GetPreStepPoint()->GetKineticEnergy();
Infos[5] = aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit(); Infos[5] = aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit();
m_SegmentNumber = (int)((m_Position.z() + m_ScorerLength / 2.) / m_SegmentLength ) + 1; //Pad number m_SegmentNumber = (int)((m_Position.z() + m_ScorerLength / 2.) / m_SegmentLength) + 1; // Pad number
Infos[6] = m_SegmentNumber; Infos[6] = m_SegmentNumber;
Infos[7] = m_Position.z(); Infos[7] = m_Position.z();
if(m_Shape == "Cylindrical") Infos[8] = pow(pow(m_Position.x(),2) + pow(m_Position.y(),2),0.5); //R if (m_Shape == "Cylindrical")
if(m_Shape == "Long_Chamber") Infos[8] = m_Position.y(); Infos[8] = pow(pow(m_Position.x(), 2) + pow(m_Position.y(), 2), 0.5); // R
if (m_Shape == "Long_Chamber")
Infos[8] = m_Position.y();
Infos[9] = aStep->GetTrack()->GetVertexPosition()[2]; Infos[9] = aStep->GetTrack()->GetVertexPosition()[2];
Infos[10] = aStep->GetTrack()->GetVertexKineticEnergy(); Infos[10] = aStep->GetTrack()->GetVertexKineticEnergy();
G4ThreeVector p_vec = aStep->GetTrack()->GetVertexMomentumDirection(); G4ThreeVector p_vec = aStep->GetTrack()->GetVertexMomentumDirection();
Infos[11] = acos(p_vec[2]/pow(pow(p_vec[0],2)+pow(p_vec[1],2)+pow(p_vec[2],2),0.5))/deg; //angle relative to z axis (theta); Infos[11] = acos(p_vec[2] / pow(pow(p_vec[0], 2) + pow(p_vec[1], 2) + pow(p_vec[2], 2), 0.5)) /
deg; // angle relative to z axis (theta);
Infos[12] = aStep->GetTrack()->GetTrackLength(); Infos[12] = aStep->GetTrack()->GetTrackLength();
int sector = -1; int sector = -1;
if(Infos[8]/cm > 1.2) { if (Infos[8] / cm > 1.2) {
if(m_Position.x() > 0 and m_Position.y() > 0) { if (m_Position.x() > 0 and m_Position.y() > 0) {
if(abs(m_Position.x()) < abs(m_Position.y())) sector = 0; if (abs(m_Position.x()) < abs(m_Position.y()))
if(abs(m_Position.x()) > abs(m_Position.y())) sector = 1; sector = 0;
if (abs(m_Position.x()) > abs(m_Position.y()))
sector = 1;
} }
if(m_Position.x() > 0 and m_Position.y() < 0) { if (m_Position.x() > 0 and m_Position.y() < 0) {
if(abs(m_Position.x()) > abs(m_Position.y())) sector = 2; if (abs(m_Position.x()) > abs(m_Position.y()))
if(abs(m_Position.x()) < abs(m_Position.y())) sector = 3; sector = 2;
if (abs(m_Position.x()) < abs(m_Position.y()))
sector = 3;
} }
if(m_Position.x() < 0 and m_Position.y() < 0) { if (m_Position.x() < 0 and m_Position.y() < 0) {
if(abs(m_Position.x()) < abs(m_Position.y())) sector = 4; if (abs(m_Position.x()) < abs(m_Position.y()))
if(abs(m_Position.x()) > abs(m_Position.y())) sector = 5; sector = 4;
if (abs(m_Position.x()) > abs(m_Position.y()))
sector = 5;
} }
if(m_Position.x() < 0 and m_Position.y() > 0) { if (m_Position.x() < 0 and m_Position.y() > 0) {
if(abs(m_Position.x()) > abs(m_Position.y())) sector = 6; if (abs(m_Position.x()) > abs(m_Position.y()))
if(abs(m_Position.x()) < abs(m_Position.y())) sector = 7; sector = 6;
if (abs(m_Position.x()) < abs(m_Position.y()))
sector = 7;
} }
} }
Infos[13] = sector; Infos[13] = sector;
Infos[14] = m_p0 + m_p1*Infos[8] + m_p2*Infos[8]*Infos[8] + m_p3*Infos[8]*Infos[8]*Infos[8]; Infos[14] = m_p0 + m_p1 * Infos[8] + m_p2 * Infos[8] * Infos[8] + m_p3 * Infos[8] * Infos[8] * Infos[8];
G4ThreeVector delta_Position = aStep->GetDeltaPosition(); G4ThreeVector delta_Position = aStep->GetDeltaPosition();
m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6; m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6;
if(isnan(Infos[10])) { if (std::isnan(Infos[10])) {
aStep->GetTrack()->SetTrackStatus(fStopAndKill); aStep->GetTrack()->SetTrackStatus(fStopAndKill);
return 0; return 0;
} }
if(aStep->IsFirstStepInVolume() == true) excess = 0.; if (aStep->IsFirstStepInVolume() == true)
excess = 0.;
map<G4int, G4double**>::iterator it; map<G4int, G4double**>::iterator it;
it= EvtMap->GetMap()->find(m_Index); it = EvtMap->GetMap()->find(m_Index);
if(it!=EvtMap->GetMap()->end()){ if (it != EvtMap->GetMap()->end()) {
G4double* dummy = *(it->second); G4double* dummy = *(it->second);
if(Infos[1]==dummy[1]) Infos[5]+=dummy[5]; //accumulate ionisation energy deposit to get total accross pad if (Infos[1] == dummy[1])
Infos[5] += dummy[5]; // accumulate ionisation energy deposit to get total accross pad
delete dummy; delete dummy;
} }
#ifdef USE_Garfield #ifdef USE_Garfield
Infos[15] = GARFDRIFT(((aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV+excess), Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0], (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV, m_Shape)*eV; Infos[15] = GARFDRIFT(((aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / eV + excess),
/* Infos[3], m_Position / cm, delta_Position / cm, Infos[8] / cm, Infos[6], Infos[2],
m_ScorerLength / cm, m_SegmentLength / cm, Infos[0],
(aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / eV, m_Shape) *
eV;
/*
file.open("excess_test.dat",std::ios::app); file.open("excess_test.dat",std::ios::app);
file << Infos[6] << "\t" << "\t" << aStep->IsFirstStepInVolume() << "\t" << excess << "\t" << Infos[15]/eV << "\t" << (int)((((aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV+excess) / 41.1)*0.01) << "\t" << Infos[8] << endl; file << Infos[6] << "\t" << "\t" << aStep->IsFirstStepInVolume() << "\t" << excess << "\t" << Infos[15]/eV << "\t"
file.close(); << (int)((((aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit())/eV+excess) / 41.1)*0.01) << "\t" <<
Infos[8] << endl; file.close();
*/ */
excess = Infos[15]/eV; excess = Infos[15] / eV;
#endif #endif
EvtMap->set(m_Index, Infos); EvtMap->set(m_Index, Infos);
return TRUE; return TRUE;
} }
void Gas_Scorer::Initialize(G4HCofThisEvent* HCE){ void Gas_Scorer::Initialize(G4HCofThisEvent* HCE) {
EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName()); EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
if (HCID < 0) { if (HCID < 0) {
HCID = GetCollectionID(0); HCID = GetCollectionID(0);
} }
HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
} }
void Gas_Scorer::EndOfEvent(G4HCofThisEvent*){} void Gas_Scorer::EndOfEvent(G4HCofThisEvent*) {}
void Gas_Scorer::clear(){ void Gas_Scorer::clear() {
std::map<G4int, G4double**>::iterator MapIterator; std::map<G4int, G4double**>::iterator MapIterator;
for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){ for (MapIterator = EvtMap->GetMap()->begin(); MapIterator != EvtMap->GetMap()->end(); MapIterator++) {
delete *(MapIterator->second); delete *(MapIterator->second);
} }
EvtMap->clear(); EvtMap->clear();
} }
void Gas_Scorer::DrawAll() {}
void Gas_Scorer::DrawAll(){} //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Gas_Scorer::PrintAll(){ void Gas_Scorer::PrintAll() {
G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
G4cout << " PrimitiveScorer " << GetName() << G4endl ; G4cout << " PrimitiveScorer " << GetName() << G4endl;
G4cout << " Number of entries " << EvtMap->entries() << G4endl ; G4cout << " Number of entries " << EvtMap->entries() << G4endl;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment