389 lines
16 KiB
C++
389 lines
16 KiB
C++
// Geant4 simulation for MuSR
|
|
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
|
|
// DATE : 2008-05
|
|
//
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
//cksdel #include "lem4RunAction.hh"
|
|
#include "lem4ScintSD.hh"
|
|
#include "G4HCofThisEvent.hh"
|
|
#include "G4Step.hh"
|
|
#include "G4ThreeVector.hh"
|
|
#include "G4SDManager.hh"
|
|
#include "G4ios.hh"
|
|
//#include "G4MagneticField.hh"
|
|
//#include "G4FieldManager.hh"
|
|
//#include "G4TransportationManager.hh"
|
|
#include <algorithm> // needed for the sort() function
|
|
#include "G4VProcess.hh" // needed for the degugging message of the process name
|
|
#include "G4RunManager.hh"
|
|
#include "G4Run.hh"
|
|
#include "lem4Parameters.hh"
|
|
#include "lem4ErrorMessage.hh"
|
|
#include <vector>
|
|
|
|
//bool myREMOVEfunction (int i,int j) { return (i<j); }
|
|
//bool timeOrdering (lem4ScintHit hit1, lem4ScintHit hit2) {
|
|
// return (hit1.GetGlobalTime()<hit2.GetGlobalTime());
|
|
//}
|
|
//
|
|
//bool timeOrdering2 (std::map<int,double>::iterator i1, std::map<int,double>::iterator m2) {
|
|
// return ( (*i1).first()<(*i2).second() );
|
|
//}
|
|
//
|
|
//bool timeOrdering2 (std::pair<int,double> p1, std::pair<int,double> p2) {
|
|
// return ( p1.first()<p2.second() );
|
|
//}
|
|
//
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
lem4ScintSD::lem4ScintSD(G4String name)
|
|
:G4VSensitiveDetector(name)
|
|
{
|
|
G4String HCname;
|
|
collectionName.insert(HCname="scintCollection");
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
lem4ScintSD::~lem4ScintSD(){ }
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void lem4ScintSD::Initialize(G4HCofThisEvent* HCE) {
|
|
if (verboseLevel>1) G4cout<<"VERBOSE 2: lem4ScintSD::Initialize\n";
|
|
// Positron_momentum_already_stored=0;
|
|
// comment the following lines to disable output of initial variables
|
|
//lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance();
|
|
//myRootOutput->ClearAllRootVariables();
|
|
|
|
scintCollection = new lem4ScintHitsCollection
|
|
(SensitiveDetectorName,collectionName[0]);
|
|
static G4int HCID = -1;
|
|
if(HCID<0) {
|
|
HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
|
|
if (verboseLevel>1) G4cout<<"VERBOSE 2: lem4ScintSD::HCID was <0\n, now HCID="<<HCID<<"\n";
|
|
}
|
|
HCE->AddHitsCollection( HCID, scintCollection );
|
|
myStoreOnlyEventsWithHits = lem4Parameters::storeOnlyEventsWithHits;
|
|
mySignalSeparationTime = lem4Parameters::signalSeparationTime;
|
|
myStoreOnlyTheFirstTimeHit= lem4Parameters::storeOnlyTheFirstTimeHit;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4bool lem4ScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
|
|
{
|
|
if (verboseLevel>1) G4cout<<"VERBOSE 2: lem4ScintSD::ProcessHits\n";
|
|
G4double edep = aStep->GetTotalEnergyDeposit();
|
|
if(edep==0.) {
|
|
return false;
|
|
}
|
|
|
|
G4Track* aTrack = aStep->GetTrack();
|
|
G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName();
|
|
// G4cout <<"actualVolume=="<<actualVolume<<G4endl;
|
|
///if (actualVolume=="Target") { // Perhaps not necessary! TS.
|
|
///return false;
|
|
///}
|
|
// else if ((strncmp(actualVolume.c_str(),"Shield",6)==0)||(strncmp(actualVolume.c_str(),"shield",6)==0)) {
|
|
// delete track if the particle is far away from the detector (i.e. in the "shield" volume)
|
|
// aTrack->SetTrackStatus(fSuspend);
|
|
// }
|
|
|
|
// If requested, store only the hit that happened first (usefull for some special studies, not for a serious simulation)
|
|
if (myStoreOnlyTheFirstTimeHit) {
|
|
G4int NbHits = scintCollection->entries();
|
|
if (NbHits>0) {
|
|
aTrack->SetTrackStatus(fStopAndKill);
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
//cks Ignore hits from the particles that had been created at the time when the muon
|
|
// was stopping in the target - not a correct thing to do, however no better way
|
|
// to suppres wrong timing information of the hit found yet.
|
|
// lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance();
|
|
// G4double muonArrivedInTarget = myRootOutput->GetTimeInTarget();
|
|
// if (muonArrivedInTarget>0.) {
|
|
// G4double timeOfThisTrack = aTrack->GetGlobalTime();
|
|
// if (fabs(muonArrivedInTarget-timeOfThisTrack)<(2*ns)) {
|
|
// return false;
|
|
// }
|
|
// }
|
|
// else {
|
|
// // G4cout<<"Kamil: in the hit before muon has arrived into the target"<<G4endl;
|
|
// }
|
|
// if ( (aTrack->GetGlobalTime()) < (0.03*microsecond) ) return false;
|
|
lem4ScintHit* newHit = new lem4ScintHit();
|
|
// G4RunManager* fRunManager = G4RunManager::GetRunManager();
|
|
// newHit->SetEventID(fRunManager->GetCurrentEvent()->GetEventID());
|
|
// newHit->SetRunID(fRunManager->GetCurrentRun()->GetRunID());
|
|
newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
|
|
newHit->SetParticleID (aTrack->GetDefinition()->GetPDGEncoding());
|
|
// if (aTrack->GetDefinition()->GetPDGEncoding()==-13) {
|
|
// if ((aTrack->GetGlobalTime())<(0.015*microsecond)) {
|
|
// G4cout<<"PROCESS: detected particle="<<aTrack->GetDefinition()->GetParticleName()
|
|
// <<" TIME="<<aTrack->GetGlobalTime()<<G4endl;
|
|
// const G4VProcess* interestingProcess = aStep->GetPostStepPoint()->GetProcessDefinedStep();
|
|
// // G4cout<<"GetProcessName()="<<interestingProcess->GetProcessName()<<G4endl;
|
|
// const G4VProcess* creatorProcess =aTrack->GetCreatorProcess();
|
|
// G4cout<<"Process that limitted this step: "; interestingProcess->DumpInfo();
|
|
// if (creatorProcess!=NULL) {G4cout<<"Process that created this particle: "; creatorProcess->DumpInfo();}
|
|
// }
|
|
newHit->SetTrackID (aTrack->GetTrackID());
|
|
newHit->SetEdep (edep);
|
|
newHit->SetPrePos (aStep->GetPreStepPoint()->GetPosition());
|
|
newHit->SetPostPos (aStep->GetPostStepPoint()->GetPosition());
|
|
newHit->SetPol (aTrack->GetPolarization());
|
|
newHit->SetLogVolName(aTrack->GetVolume()->GetLogicalVolume()->GetName());
|
|
newHit->SetGlobTime (aTrack->GetGlobalTime());
|
|
|
|
/// Printout to check proper assignment of Global Time!
|
|
//G4cout<<"StepTime="<<aStep->GetDeltaTime()<<"\t";
|
|
//G4cout<<"GlobalTime="<<aTrack->GetGlobalTime()<<" LocalTime="<<aTrack->GetLocalTime()<<" ProperTime="<<aTrack->GetProperTime();
|
|
//G4cout<<"\t trackID="<<aTrack->GetTrackID() <<"\t actualVolume="<<actualVolume<<G4endl;
|
|
|
|
// Warning - aStep->IsFirstStepInVolume() only available in Geant version >= 4.8.2 !
|
|
// newHit->SetFirstStepInVolumeFlag(aStep->IsFirstStepInVolume());
|
|
// newHit->SetLastStepInVolumeFlag(aStep->IsLastStepInVolume());
|
|
newHit->SetKineticEnergy(aTrack->GetVertexKineticEnergy());
|
|
// newHit->SetKineticEnergy(aTrack->GetKineticEnergy()+edep);
|
|
newHit->SetStepLength (aStep->GetStepLength());
|
|
// G4double BFieldAtOrigin[6] = {0.,0.,0.,0.,0.,0.};
|
|
// G4double Origin[4] ={0.,0.,0.,0.};
|
|
// G4FieldManager *fMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
|
|
// if (fMgr!=NULL) {
|
|
// if(!fMgr->DoesFieldChangeEnergy()) { //then we have a magnetic field
|
|
// fMgr->GetDetectorField()->GetFieldValue(Origin,BFieldAtOrigin);
|
|
// }
|
|
// else{
|
|
// }
|
|
// // G4cout<<"Kamil: pole="<<BFieldAtOrigin[0]/tesla<<" "<<BFieldAtOrigin[1]/tesla<<" "<<BFieldAtOrigin[2]/tesla
|
|
// // <<" "<<BFieldAtOrigin[3]/tesla<<" "<<BFieldAtOrigin[4]/tesla<<" "<<BFieldAtOrigin[5]/tesla<<G4endl;
|
|
// }
|
|
// newHit->SetBField(BFieldAtOrigin);
|
|
scintCollection->insert( newHit );
|
|
// newHit->Print();
|
|
newHit->Draw();
|
|
return true;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void lem4ScintSD::EndOfEvent(G4HCofThisEvent*) {
|
|
if (verboseLevel>1) {
|
|
G4cout<<"VERBOSE 2: lem4ScintSD::EndOfEvent"<<G4endl;
|
|
G4int NbHits = scintCollection->entries();
|
|
G4cout << "\n-------->Hits Collection: in this event they are " << NbHits
|
|
<< " hits in the scint chambers: " << G4endl;
|
|
// for (G4int i=0;i<NbHits;i++) (*scintCollection)[i]->Print();
|
|
}
|
|
|
|
// Positron_momentum_already_stored=0;
|
|
lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance();
|
|
|
|
G4RunManager* fRunManager = G4RunManager::GetRunManager();
|
|
myRootOutput->SetRunID(fRunManager->GetCurrentRun()->GetRunID());
|
|
myRootOutput->SetEventID(fRunManager->GetCurrentEvent()->GetEventID());
|
|
|
|
G4int NbHits = scintCollection->entries();
|
|
if ((NbHits<=0)&&(myStoreOnlyEventsWithHits)) {
|
|
myRootOutput->ClearAllRootVariables();
|
|
return;
|
|
}
|
|
|
|
// Sort out hits and fill them into root
|
|
if (NbHits>0) {
|
|
lem4ScintHit* firstHit = (*scintCollection)[0];
|
|
myRootOutput->SetGlobTime(firstHit->GetGlobalTime());
|
|
|
|
const G4int det_IDmax = lem4RootOutput::det_nMax;
|
|
G4double det_edep[det_IDmax];
|
|
G4int det_nsteps[det_IDmax];
|
|
G4double det_length[det_IDmax];
|
|
G4int det_ID[det_IDmax];
|
|
G4double det_edep_el[det_IDmax];
|
|
G4double det_edep_pos[det_IDmax];
|
|
G4double det_edep_gam[det_IDmax];
|
|
G4double det_edep_mup[det_IDmax];
|
|
G4double det_time_start[det_IDmax];
|
|
G4double det_time_end[det_IDmax];
|
|
G4double det_x[det_IDmax];
|
|
G4double det_y[det_IDmax];
|
|
G4double det_z[det_IDmax];
|
|
// for (G4int ii=0; ii<det_IDmax; ii++) {
|
|
// det_edep[ii]=0.;
|
|
// det_nsteps[ii]=0;
|
|
// det_length[ii]=0.;
|
|
// det_ID[ii]=-1000;
|
|
// det_edep_el[ii]=0.;
|
|
// det_edep_pos[ii]=0.;
|
|
// det_edep_gam[ii]=0.;
|
|
// det_edep_mup[ii]=0.;
|
|
// }
|
|
|
|
// const G4int subtrack_max=lem4RootOutput::maxNSubTracks;
|
|
// G4int particleID[subtrack_max];
|
|
// G4int trackID[subtrack_max];
|
|
// G4int logicalVolumeID[subtrack_max];
|
|
// G4double edep[subtrack_max];
|
|
// G4double x[subtrack_max];
|
|
// G4double y[subtrack_max];
|
|
// G4double z[subtrack_max];
|
|
// G4double post_x[subtrack_max];
|
|
// G4double post_y[subtrack_max];
|
|
// G4double post_z[subtrack_max];
|
|
// G4double kineticEnergy[subtrack_max];
|
|
// G4double stepLength[subtrack_max];
|
|
// G4int nrOfSubsteps[subtrack_max];
|
|
// for (G4int ii=0; ii<subtrack_max; ii++) {
|
|
// particleID[ii]=-1000;
|
|
// trackID[ii]=-1000;
|
|
// logicalVolumeID[ii]=-1000;
|
|
// edep[ii]=-1000;
|
|
// x[ii]=-1000;
|
|
// y[ii]=-1000;
|
|
// z[ii]=-1000;
|
|
// post_x[ii]=-1000;
|
|
// post_y[ii]=-1000;
|
|
// post_z[ii]=-1000;
|
|
// kineticEnergy[ii]=-1000;
|
|
// stepLength[ii]=-1000;
|
|
// nrOfSubsteps[ii]=-1000;
|
|
// }
|
|
|
|
// G4String aLogicalVolume_old="Unset";
|
|
// G4String aParticleName_old="Unset";
|
|
// G4int det_index=-1;
|
|
// G4int subtrack_index=-1;
|
|
// G4cout<<G4endl;
|
|
|
|
// int myints[] = {32,71,12,45,26,80,53,33};
|
|
// std::vector<int> myvector (myints, myints+8); // 32 71 12 45 26 80 53 33
|
|
// std::vector<int>::iterator it;
|
|
// sort (myvector.begin(), myvector.end(), myREMOVEfunction);
|
|
// // print out content:
|
|
// G4cout << "myvector contains:";
|
|
// for (it=myvector.begin(); it!=myvector.end(); ++it)
|
|
// G4cout << " " << *it;
|
|
// G4cout << G4endl;
|
|
|
|
|
|
// Sort hits according to the time. Using std::map is convenient, because it sorts
|
|
// its entries according to the key (the first variable of the pair).
|
|
std::map<G4double,G4int> myHitTimeMapping;
|
|
std::map<G4double,G4int>::iterator it;
|
|
for (G4int i=0; i<NbHits; i++) {
|
|
lem4ScintHit* aHit = (*scintCollection)[i];
|
|
G4double tmptime=aHit->GetGlobalTime();
|
|
// G4cout<<"Hit nr "<<i<<" at time="<<tmptime<<" with edep="<<aHit->GetEdep()/MeV
|
|
// <<" detID="<<myRootOutput->ConvertVolumeToID(aHit->GetLogVolName())<< G4endl;
|
|
myHitTimeMapping.insert ( std::pair<G4double,G4int>(tmptime,i) );
|
|
}
|
|
|
|
// Loop over all hits (which are sorted according to their time):
|
|
G4int nSignals=0;
|
|
for (it=myHitTimeMapping.begin(); it!=myHitTimeMapping.end(); it++) {
|
|
// G4cout << "Key:" << it->first;
|
|
// G4cout << " Value:" << it->second << "\n";
|
|
G4int ii = it->second; // ii is the index of the hits, which is sorted according to time
|
|
lem4ScintHit* aHit = (*scintCollection)[ii];
|
|
G4String aHitVolumeName = aHit->GetLogVolName();
|
|
G4int aHitVolumeID = myRootOutput->ConvertVolumeToID(aHitVolumeName);
|
|
G4double aHitTime = aHit->GetGlobalTime();
|
|
|
|
// Loop over all already defined signals and check whether the hit falls into any of them
|
|
G4bool signalAssigned=false;
|
|
for (G4int j=0; j<nSignals; j++) {
|
|
if ( (aHitVolumeID==det_ID[j]) && ((aHitTime-det_time_end[j])<mySignalSeparationTime) ) {
|
|
signalAssigned=true;
|
|
det_edep[j] += aHit->GetEdep();
|
|
det_nsteps[j]++;
|
|
det_length[j] += aHit->GetStepLength();
|
|
det_time_end[j] = aHitTime;
|
|
G4String aParticleName = aHit->GetParticleName();
|
|
if (aParticleName=="e-") {
|
|
det_edep_el[j] += aHit->GetEdep();
|
|
} else if (aParticleName=="e+") {
|
|
det_edep_pos[j] += aHit->GetEdep();
|
|
} else if (aParticleName=="gamma") {
|
|
det_edep_gam[j] += aHit->GetEdep();
|
|
} else if (aParticleName=="mu+") {
|
|
det_edep_mup[j] += aHit->GetEdep();
|
|
} else {
|
|
char message[200];
|
|
sprintf(message,"lem4ScintSD.cc::EndOfEvent(): untreated particle \"%s\" deposited energy.",aParticleName.c_str());
|
|
lem4ErrorMessage::GetInstance()->lem4Error(WARNING,message,true);
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
if (!signalAssigned) { // The hit does not belong to any existing signal --> create a new signal.
|
|
// Check, whether the maximum number of signals was not exceeded:
|
|
if ( nSignals >= (det_IDmax-1) ) {
|
|
char message[200];
|
|
sprintf(message,"lem4ScintSD.cc::EndOfEvent(): number of signals exceeds maximal allowed value.");
|
|
lem4ErrorMessage::GetInstance()->lem4Error(WARNING,message,true);
|
|
}
|
|
else {
|
|
det_edep[nSignals] = aHit->GetEdep();
|
|
det_nsteps[nSignals] = 1;
|
|
det_length[nSignals] = aHit->GetStepLength();
|
|
det_ID[nSignals] = aHitVolumeID;
|
|
det_time_start[nSignals] = aHitTime;
|
|
det_time_end[nSignals] = aHitTime;
|
|
det_edep_el[nSignals] = 0;
|
|
det_edep_pos[nSignals] = 0;
|
|
det_edep_gam[nSignals] = 0;
|
|
det_edep_mup[nSignals] = 0;
|
|
G4String aParticleName = aHit->GetParticleName();
|
|
if (aParticleName=="e-") {
|
|
det_edep_el[nSignals] += aHit->GetEdep();
|
|
} else if (aParticleName=="e+") {
|
|
det_edep_pos[nSignals] += aHit->GetEdep();
|
|
} else if (aParticleName=="gamma") {
|
|
det_edep_gam[nSignals] += aHit->GetEdep();
|
|
} else if (aParticleName=="mu+") {
|
|
det_edep_mup[nSignals] += aHit->GetEdep();
|
|
} else {
|
|
char message[200];
|
|
sprintf(message,"lem4ScintSD.cc::EndOfEvent(): UNTREATED PARTICLE \"%s\" deposited energy.",aParticleName.c_str());
|
|
lem4ErrorMessage::GetInstance()->lem4Error(WARNING,message,true);
|
|
}
|
|
G4ThreeVector prePos = aHit->GetPrePos();
|
|
det_x[nSignals]=prePos.x();
|
|
det_y[nSignals]=prePos.y();
|
|
det_z[nSignals]=prePos.z();
|
|
nSignals++;
|
|
}
|
|
} // end of "if (!signalAssigned)"
|
|
} // end of the for loop over the hits
|
|
|
|
// Sort the signals according to the energy (in decreasing order)
|
|
std::map<G4double,G4int> mySignalMapping;
|
|
std::map<G4double,G4int>::iterator itt;
|
|
for (G4int i=0; i<nSignals; i++) {
|
|
mySignalMapping.insert ( std::pair<G4double,G4int>(-det_edep[i],i) );
|
|
}
|
|
|
|
// Write out the signals (sorted according to energy) to the lem4RootOutput class:
|
|
G4int jj=-1;
|
|
for (itt=mySignalMapping.begin(); itt!=mySignalMapping.end(); itt++) {
|
|
jj++;
|
|
G4int ii = itt->second;
|
|
myRootOutput->SetDetectorInfo(jj,det_ID[ii],det_edep[ii],det_edep_el[ii],det_edep_pos[ii],
|
|
det_edep_gam[ii],det_edep_mup[ii],det_nsteps[ii],det_length[ii],
|
|
det_time_start[ii],det_time_end[ii],det_x[ii],det_y[ii],det_z[ii]);
|
|
}
|
|
|
|
} //end "if (NbHits>0)"
|
|
|
|
myRootOutput->FillEvent();
|
|
myRootOutput->ClearAllRootVariables();
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|