diff --git a/NPSimulation/Detectors/TACTIC/TACTICScorer.cc b/NPSimulation/Detectors/TACTIC/TACTICScorer.cc index 1c64e10e04bcb679a0989778255f0acff4a28722..c874685b18c6516f034c2da3a9d1b5c867a90e72 100644 --- a/NPSimulation/Detectors/TACTIC/TACTICScorer.cc +++ b/NPSimulation/Detectors/TACTIC/TACTICScorer.cc @@ -1,8 +1,8 @@ //#define USE_Garfield //only use if compiled with Garfield #include "TACTICScorer.hh" -#include "G4UnitsTable.hh" #include "G4RunManager.hh" +#include "G4UnitsTable.hh" #include "TACTIC.hh" #ifdef USE_Garfield @@ -13,13 +13,14 @@ double excess; 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? -:G4VPrimitiveScorer(name, depth),HCID(-1){ +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? + : G4VPrimitiveScorer(name, depth), HCID(-1) { m_ScorerLength = ScorerLength; m_NumberOfSegments = NumberOfSegments; m_SegmentLength = m_ScorerLength / m_NumberOfSegments; m_Level = Level; - m_Position = G4ThreeVector(-1000,-1000,-1000); + m_Position = G4ThreeVector(-1000, -1000, -1000); m_SegmentNumber = -1; m_Index = -1; m_p0 = p0; @@ -29,122 +30,139 @@ Gas_Scorer::Gas_Scorer(G4String name,G4int Level,G4double ScorerLength,G4int Num 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]; - m_Position = aStep->GetPreStepPoint()->GetPosition(); - + m_Position = aStep->GetPreStepPoint()->GetPosition(); + Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID(); Infos[1] = aStep->GetTrack()->GetTrackID(); - Infos[2] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();; - + Infos[2] = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber(); + ; + Infos[3] = aStep->GetPreStepPoint()->GetGlobalTime(); Infos[4] = aStep->GetPreStepPoint()->GetKineticEnergy(); 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[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 == "Long_Chamber") Infos[8] = m_Position.y(); + if (m_Shape == "Cylindrical") + 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[10] = aStep->GetTrack()->GetVertexKineticEnergy(); 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(); int sector = -1; - if(Infos[8]/cm > 1.2) { - 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())) sector = 1; + if (Infos[8] / cm > 1.2) { + 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())) + sector = 1; } - 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())) sector = 3; + 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())) + sector = 3; } - 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())) sector = 5; + 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())) + sector = 5; } - 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())) sector = 7; + 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())) + sector = 7; } } - + 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(); m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level); m_Index = m_DetectorNumber * 1e3 + m_SegmentNumber * 1e6; - - if(isnan(Infos[10])) { + + if (std::isnan(Infos[10])) { aStep->GetTrack()->SetTrackStatus(fStopAndKill); return 0; } - if(aStep->IsFirstStepInVolume() == true) excess = 0.; - + if (aStep->IsFirstStepInVolume() == true) + excess = 0.; + map<G4int, G4double**>::iterator it; - it= EvtMap->GetMap()->find(m_Index); - if(it!=EvtMap->GetMap()->end()){ + it = EvtMap->GetMap()->find(m_Index); + if (it != EvtMap->GetMap()->end()) { 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; } - + #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 << 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.close(); + 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.close(); */ - excess = Infos[15]/eV; - + excess = Infos[15] / eV; + #endif EvtMap->set(m_Index, Infos); return TRUE; - } -void Gas_Scorer::Initialize(G4HCofThisEvent* HCE){ - EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName()); - if (HCID < 0) { - HCID = GetCollectionID(0); - } - HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +void Gas_Scorer::Initialize(G4HCofThisEvent* HCE) { + EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); } -void Gas_Scorer::EndOfEvent(G4HCofThisEvent*){} +void Gas_Scorer::EndOfEvent(G4HCofThisEvent*) {} -void Gas_Scorer::clear(){ - std::map<G4int, G4double**>::iterator MapIterator; - for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){ - delete *(MapIterator->second); - } +void Gas_Scorer::clear() { + std::map<G4int, G4double**>::iterator MapIterator; + for (MapIterator = EvtMap->GetMap()->begin(); MapIterator != EvtMap->GetMap()->end(); MapIterator++) { + 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(){ - G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; - G4cout << " PrimitiveScorer " << GetName() << G4endl ; - G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +void Gas_Scorer::PrintAll() { + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl; + G4cout << " PrimitiveScorer " << GetName() << G4endl; + G4cout << " Number of entries " << EvtMap->entries() << G4endl; }