LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
musrsim/src/musrScintSD.cc
Kamil Sedlak d11e19a7d8 31.1.2011 Kamil Sedlak - changes for Geant 4.9.4 
Changes in musrPhysicsList.cc that were needed for new version of 
Geant = Geant4.9.4.  The changes remove G4MultipleScattering and
MultipleAndCoulombScattering,  because multiple scattering has to
be treated differently for different particles in Geant4.9.4 and
"G4MultipleScattering" is not supported anymore.  The documentation
has to be updated (i.e. comments about MultipleAndCoulombScattering
has to be removed from doc/musrSim.tex).
2011-01-31 10:40:44 +00:00

638 lines
28 KiB
C++
Raw Blame History

/***************************************************************************
* musrSim - the program for the simulation of (mainly) muSR instruments. *
* More info on http://lmu.web.psi.ch/simulation/index.html . *
* musrSim is based od Geant4 (http://geant4.web.cern.ch/geant4/) *
* *
* Copyright (C) 2009 by Paul Scherrer Institut, 5232 Villigen PSI, *
* Switzerland *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the Free Software *
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *
***************************************************************************/
#include "musrScintSD.hh"
#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
#include "G4ThreeVector.hh"
#include "G4SDManager.hh"
#include "G4ios.hh"
#include <algorithm> // needed for the sort() function
#include "G4VProcess.hh" // needed for the degugging message of the process name
#include "G4OpBoundaryProcess.hh" // Optical photon process
#include "G4RunManager.hh"
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
#include "musrSteppingAction.hh"
//#include "TCanvas.h"
#include "TMath.h"
#include "TF1.h"
#include <vector>
//bool myREMOVEfunction (int i,int j) { return (i<j); }
//bool timeOrdering (musrScintHit hit1, musrScintHit 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......
Double_t poissonf(Double_t* x, Double_t* par) {
return par[0]*TMath::Poisson(x[0],par[1]);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintSD::musrScintSD(G4String name)
:G4VSensitiveDetector(name)
{
pointer=this;
G4String HCname;
collectionName.insert(HCname="scintCollection");
OPSA_minNrOfDetectedPhotons = 1;
OPSA_signalSeparationTime = 10.;
OPSA_fracA = 0.01;
OPSA_fracB = 0.05;
OPSA_fracC = 0.10;
OPSA_fracD = 0.20;
OPSA_fracE = 0.5;
bool_multimapOfEventIDsForOPSAhistosEXISTS = false;
OPSAhistoNbin = 100;
OPSAhistoMin =0;
OPSAhistoMax = 10.;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintSD::~musrScintSD(){ }
musrScintSD* musrScintSD::pointer=NULL;
musrScintSD* musrScintSD::GetInstance() {return pointer;}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::Initialize\n";
scintCollection = new musrScintHitsCollection
(SensitiveDetectorName,collectionName[0]);
static G4int HCID = -1;
if(HCID<0) {
HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::HCID was <0\n, now HCID="<<HCID<<"\n";
}
HCE->AddHitsCollection( HCID, scintCollection );
myStoreOnlyEventsWithHits = musrParameters::storeOnlyEventsWithHits;
mySignalSeparationTime = musrParameters::signalSeparationTime;
myStoreOnlyTheFirstTimeHit= musrParameters::storeOnlyTheFirstTimeHit;
myStoreOnlyEventsWithHitInDetID = musrParameters::storeOnlyEventsWithHitInDetID;
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
boolIsVvvInfoRequested = myMusrSteppingAction->IsVvvInfoRequested();
myRootOutput = musrRootOutput::GetRootInstance();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
{
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::ProcessHits\n";
G4double edep = aStep->GetTotalEnergyDeposit();
if(edep==0.) {
return false;
}
G4Track* aTrack = aStep->GetTrack();
// G4LogicalVolume* hitLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
// G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName();
G4String hitLogicalVolumeName = aTrack->GetVolume()->GetLogicalVolume()->GetName();
G4String particleName=aTrack->GetDefinition()->GetParticleName();
// if (particleName=="opticalphoton") {G4cout<<"UFON JE TU: edep="<<edep<<G4endl; return false;}
if (particleName=="opticalphoton") {
ProcessOpticalPhoton(aStep);
return false;
}
// 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;
}
}
musrScintHit* newHit = new musrScintHit();
// newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
newHit->SetParticleName(particleName);
G4int particleID = aTrack->GetDefinition()->GetPDGEncoding();
newHit->SetParticleID (particleID);
newHit->SetEdep (edep);
newHit->SetPrePos (aStep->GetPreStepPoint()->GetPosition());
newHit->SetPostPos (aStep->GetPostStepPoint()->GetPosition());
newHit->SetPol (aTrack->GetPolarization());
// G4LogicalVolume* hitLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
newHit->SetLogVolName(hitLogicalVolumeName);
// newHit->SetLogVolName(hitLogicalVolume->GetName());
newHit->SetGlobTime(aTrack->GetGlobalTime());
// Warning - aStep->IsFirstStepInVolume() only available in Geant version >= 4.8.2 !
// newHit->SetFirstStepInVolumeFlag(aStep->IsFirstStepInVolume());
// newHit->SetLastStepInVolumeFlag(aStep->IsLastStepInVolume());
newHit->SetKineticEnergy(aTrack->GetKineticEnergy());
// newHit->SetKineticEnergy(aTrack->GetKineticEnergy()+edep);
G4double vertexKineticEnergy = aTrack->GetVertexKineticEnergy();
newHit->SetVertexKineticEnergy(vertexKineticEnergy);
G4ThreeVector vertexPosition = aTrack->GetVertexPosition();
newHit->SetVertexPosition(vertexPosition);
const G4LogicalVolume* vertexLogicalVolume = aTrack->GetLogicalVolumeAtVertex();
G4String vertexLogicalVolumeName = vertexLogicalVolume->GetName();
newHit->SetLogicalVolumeAtVertex(vertexLogicalVolumeName);
G4String processName;
if ((aTrack->GetCreatorProcess())!=0) { processName=aTrack->GetCreatorProcess()->GetProcessName(); }
else {processName="initialParticle";} //if no process found, the track comes from the generated particle
newHit->SetCreatorProcessName(processName);
G4int trackID = aTrack->GetTrackID();
newHit->SetTrackID (trackID);
newHit->SetStepLength (aStep->GetStepLength());
scintCollection->insert( newHit );
// newHit->Print();
newHit->Draw();
return true;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
// //Was the photon absorbed by the absorption process ?
// const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
// G4String processName = (process) ? process->GetProcessName() : "Unknown";
// G4cout<<"processName="<<processName<<G4endl;
// if (processName=="OpAbsorption") {
// G4cout<<"\n ABSORPTION\n";
// }
G4OpBoundaryProcessStatus boundaryStatus=Undefined;
static G4OpBoundaryProcess* boundaryProc=NULL;
//find the boundary process only once
if(!boundaryProc){
G4ProcessManager* pm = aStep->GetTrack()->GetDefinition()->GetProcessManager();
G4int nprocesses = pm->GetProcessListLength();
G4ProcessVector* pv = pm->GetProcessList();
for(G4int i=0; i<nprocesses; i++){
if((*pv)[i]->GetProcessName()=="OpBoundary"){
boundaryProc = (G4OpBoundaryProcess*)(*pv)[i];
break;
}
}
}
boundaryStatus=boundaryProc->GetStatus();
//Check to see if the partcile was actually at a boundary
//Otherwise the boundary status may not be valid
//Prior to Geant4.6.0-p1 this would not have been enough to check
if(aStep->GetPostStepPoint()->GetStepStatus()==fGeomBoundary){
// G4cout<<" boundaryStatus="<<boundaryStatus<<" ";
// G4String actualVolume = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
G4String actualVolume = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
Int_t detID = myRootOutput->ConvertVolumeToID(actualVolume);
// G4cout<<" detID ="<<detID<<" actualVolume="<<actualVolume<<G4endl;
optHitMapType::iterator iter = optHitMap.find(detID);
if (iter==optHitMap.end()) { // optHitDetectorMapType does not exist ==> create it
optHitDetectorMapType* optHitDetectorMapTmp = new optHitDetectorMapType;
optHitMap.insert(std::pair<Int_t,optHitDetectorMapType*>(detID,optHitDetectorMapTmp));
iter = optHitMap.find(detID);
}
optHitDetectorMapType* optHitDetectorMap = (*iter).second;
// optHitDetectorMapType optHitDetectorMap = optHitMap[detID];
G4double tmpTime = aStep->GetPreStepPoint()->GetGlobalTime();
// G4cout<<" tmpTime="<<tmpTime;
// optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
if (boundaryStatus!=Detection) {
char message[200];
sprintf(message,"musrScintSD.cc::ProcessOpticalPhoton(): Optical photon boundary status is not Detection but %i",boundaryStatus);
// musrErrorMessage::GetInstance()->musrError(FATAL,message,false);
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
// switch(boundaryStatus){
// case Absorption:
// G4cout<<" AAAAAAAAAAAAAAAAAA Absorption"<<G4endl;
// break;
// case Detection: //Note, this assumes that the volume causing detection
//is the photocathode because it is the only one with
//non-zero efficiency
optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,boundaryStatus));
// G4cout<<" Detection"<<G4endl;
// break;
// case FresnelReflection:
// case LambertianReflection:
// case LobeReflection:
// case TotalInternalReflection:
// case SpikeReflection:
// G4cout<<" Reflection"<<G4endl;
// break;
// case StepTooSmall:
// G4cout<<" StepTooSmall"<<G4endl;
// break;
// case NoRINDEX:
// G4cout<<" NoRINDEX"<<G4endl;
// break;
// // case :
// // G4cout<<" "<<G4endl;
// // break;
// default:
// G4cout<<" SomethingElse"<<G4endl;
// break;
// }
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
if (verboseLevel>1) {
G4cout<<"VERBOSE 2: musrScintSD::EndOfEvent"<<G4endl;
G4int NbHits = scintCollection->entries();
G4cout << "\n-------->Hits Collection: in this event they are " << NbHits
<< " hits in the scint chambers: " << G4endl;
}
G4int NbHits = scintCollection->entries();
if (myStoreOnlyEventsWithHits) {
if (NbHits<=0) {
return;
}
else if (myStoreOnlyEventsWithHitInDetID!=0) {
for (G4int i=0; i<NbHits; i++) {
musrScintHit* aHit = (*scintCollection)[i];
G4String aHitVolumeName = aHit->GetLogVolName();
G4int aHitVolumeID = myRootOutput->ConvertVolumeToID(aHitVolumeName);
if (aHitVolumeID==myStoreOnlyEventsWithHitInDetID) break; // hit in the requested detector was identified
if (i==(NbHits-1)) return; // no hit identified in the requested detector
}
}
}
// Sort out hits and fill them into root
if (NbHits>0) {
const G4int det_IDmax = musrRootOutput::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];
G4double det_kine[det_IDmax];
G4double det_VrtxKine[det_IDmax];
G4double det_VrtxX[det_IDmax];
G4double det_VrtxY[det_IDmax];
G4double det_VrtxZ[det_IDmax];
G4int det_VrtxVolID[det_IDmax];
G4int det_VrtxProcID[det_IDmax];
G4int det_VrtxTrackID[det_IDmax];
G4int det_VrtxParticleID[det_IDmax];
G4int det_VvvTrackSign[det_IDmax];
// 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::multimap<G4double,G4int> myHitTimeMapping; // "map" replaced by "multimap" (needed for radioactive decay,
// in which case times are huge and due to the limited rounding
// precision become equal --> map ignores the same "keys",
// multimap does not.
std::map<G4double,G4int>::iterator it;
for (G4int i=0; i<NbHits; i++) {
musrScintHit* 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
musrScintHit* aHit = (*scintCollection)[ii];
G4String aHitVolumeName = aHit->GetLogVolName();
G4int aHitVolumeID = myRootOutput->ConvertVolumeToID(aHitVolumeName);
G4double aHitTime = aHit->GetGlobalTime();
G4int aHitTrackID = aHit->GetTrackID();
// 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+")||(aParticleName=="mu-")) {
det_edep_mup[j] += aHit->GetEdep();
} else {
char message[200];
sprintf(message,"musrScintSD.cc::EndOfEvent(): untreated particle \"%s\" deposited energy.",aParticleName.c_str());
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
// Check whether the signals consits of more then just one hit, in which case make the track ID negative:
if (abs(det_VrtxTrackID[j])!=aHitTrackID) {
det_VrtxTrackID[j]=-1*abs(det_VrtxTrackID[j]);
if (boolIsVvvInfoRequested) {
if (det_VvvTrackSign[j]==1) {
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
if (!(myMusrSteppingAction->AreTracksCommingFromSameParent(aHitTrackID,abs(det_VrtxTrackID[j]),aHitVolumeName))) {det_VvvTrackSign[j]=-1;}
}
}
}
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,"musrScintSD.cc::EndOfEvent(): number of signals exceeds maximal allowed value.");
musrErrorMessage::GetInstance()->musrError(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+")||(aParticleName=="mu-")) {
det_edep_mup[nSignals] += aHit->GetEdep();
} else {
char message[200];
sprintf(message,"musrScintSD.cc::EndOfEvent(): UNTREATED PARTICLE \"%s\" deposited energy.",aParticleName.c_str());
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
G4ThreeVector prePos = aHit->GetPrePos();
det_x[nSignals]=prePos.x();
det_y[nSignals]=prePos.y();
det_z[nSignals]=prePos.z();
det_kine[nSignals] = aHit->GetKineticEnergy();
det_VrtxKine[nSignals] = aHit->GetVertexKineticEnergy();
G4ThreeVector VrtxPos = aHit->GetVertexPosition();
det_VrtxX[nSignals] = VrtxPos.x();
det_VrtxY[nSignals] = VrtxPos.y();
det_VrtxZ[nSignals] = VrtxPos.z();
G4String logicalVolumeAtVertex = aHit->GetLogicalVolumeAtVertex();
det_VrtxVolID[nSignals] = myRootOutput->ConvertVolumeToID(logicalVolumeAtVertex);
G4String creatorProcessName = aHit->GetCreatorProcessName();
det_VrtxProcID[nSignals] = myRootOutput->ConvertProcessToID(creatorProcessName);
det_VrtxTrackID[nSignals] = aHit->GetTrackID();
det_VrtxParticleID[nSignals] = aHit->GetParticleID();
det_VvvTrackSign[nSignals] = 1;
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 musrRootOutput class:
G4int jj=-1;
for (itt=mySignalMapping.begin(); itt!=mySignalMapping.end(); itt++) {
jj++;
G4int ii = itt->second;
myRootOutput->SetDetectorInfo(jj,det_ID[ii],det_VrtxParticleID[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],
det_kine[ii],det_VrtxKine[ii],det_VrtxX[ii],det_VrtxY[ii],det_VrtxZ[ii],
det_VrtxVolID[ii],det_VrtxProcID[ii],det_VrtxTrackID[ii] );
if (boolIsVvvInfoRequested) {
G4int oldTrackID = abs(det_VrtxTrackID[ii]);
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
G4int vvvParentTrackID= -1; G4int vvvPparticleID; G4double vvvKine; G4ThreeVector vvvPosition; G4String vvvLogVol; G4String vvvProcess;
G4int vvvLogVolID=-999;
G4bool oldTrackRetrievedOK;
do {
if (vvvParentTrackID>-1) {oldTrackID=vvvParentTrackID;}
oldTrackRetrievedOK = myMusrSteppingAction->GetInfoAboutOldTrack(oldTrackID,
vvvParentTrackID, vvvPparticleID, vvvKine, vvvPosition, vvvLogVol, vvvProcess);
if (oldTrackRetrievedOK) {
// G4cout<<"musrScintSD: Old Track seems to be achieved fine -> Lets save it to Root tree"<<G4endl;
vvvLogVolID=myRootOutput->ConvertVolumeToID(vvvLogVol);
}
else {
G4cout<<" oldTrackRetrievedOK is false"<<G4endl;
oldTrackID = -999;
vvvParentTrackID = -999;
vvvPparticleID = -999;
vvvKine = -999;
vvvPosition = G4ThreeVector(-999,-999,-999);
vvvLogVol = "undefined";
vvvProcess = "undefined";
}
} while (oldTrackRetrievedOK && (vvvLogVolID==det_ID[ii]));
if (oldTrackRetrievedOK) {
G4int vvvProcessID=myRootOutput->ConvertProcessToID(vvvProcess);
myRootOutput->SetDetectorInfoVvv(jj,vvvKine,vvvPosition.x(),vvvPosition.y(),vvvPosition.z(),
vvvLogVolID,vvvProcessID,oldTrackID*det_VvvTrackSign[ii],vvvPparticleID);
}
} //end "boolIsVvvInfoRequested"
}
} //end "if (NbHits>0)"
// Analyse optical photons if they were produced
if (musrParameters::boolG4OpticalPhotons) EndOfEvent_OptiacalPhotons();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::EndOfEvent_OptiacalPhotons() {
// G4double kCarTolerance = G4GeometryTolerance::GetInstance() ->GetSurfaceTolerance();
// G4cout<<" DEBUG 10: kCarTolerance="<<kCarTolerance<<G4endl;
if (optHitMap.empty()) return;
G4RunManager* fRunManager = G4RunManager::GetRunManager();
G4int eeeventID = fRunManager->GetCurrentEvent()->GetEventID();
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
G4bool boolStoreThisOPSAhist = false;
G4int OPSA_detID= it->first;
optHitDetectorMapType* optHitDetectorMap = it->second;
// Check whether OPSA histograming of times of optical photon detection is required for this eventID.
if (bool_multimapOfEventIDsForOPSAhistosEXISTS) {
if (multimapOfEventIDsForOPSAhistos.find(eeeventID)!=multimapOfEventIDsForOPSAhistos.end()) {
// Now check whether the histogramming is required for the currently analysed detector.
std::pair<multimapOfEventIDsForOPSAhistos_Type::iterator,multimapOfEventIDsForOPSAhistos_Type::iterator> retOPSAhist;
multimapOfEventIDsForOPSAhistos_Type::iterator itOPSAhist;
retOPSAhist = multimapOfEventIDsForOPSAhistos.equal_range(eeeventID);
for (itOPSAhist = retOPSAhist.first; itOPSAhist!=retOPSAhist.second; itOPSAhist++) {
// Store histograms if the second parameters of eventsForOPSAhistos (i.e. detector ID)
// is set to 0 or corresponds to the currently analysed sensitive detector.
if ( (itOPSAhist->second == 0) || (itOPSAhist->second == OPSA_detID) ) boolStoreThisOPSAhist = true;
}
}
}
if (optHitDetectorMap->empty()) continue;
optHitDetectorMapType::const_iterator it2_START = optHitDetectorMap->begin();
optHitDetectorMapType::const_iterator it2_STOP = it2_START;
optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end(); it2_LAST--;
Double_t time = -1000, lastTime = -1000;
G4int OPSA_nPhot = 0;
G4double OPSA_timeFirst = -1000;
G4double OPSA_timeA = -1000;
G4double OPSA_timeB = -1000;
G4double OPSA_timeC = -1000;
G4double OPSA_timeD = -1000;
G4double OPSA_timeE = -1000;
G4double OPSA_timeLast = -1000;
G4int iHistNr = -1;
TF1* poiss = NULL;
for (optHitDetectorMapType::const_iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
OPSA_nPhot++;
lastTime = time;
time = it2->first;
if (OPSA_nPhot==1) lastTime=time; // First photon does not have a proper "lastTime" defined
if ( (it2==it2_LAST) || ((time-lastTime) > OPSA_signalSeparationTime)) {
// The iterator it2 reached last element of the map optHitDetectorMap or
// the time difference between two subsequently detected photons is too big ==> divide the signal into more signals
it2_STOP = it2;
if (it2==it2_LAST) it2_STOP++; // if we are at the end of optHitDetectorMap,
// the it2_STOP should point just behind the last map element
else OPSA_nPhot--; // if we split the optHitDetectorMap, however, the latest optical photon
// should have not be counted, because it belong to the next signal.
optHitDetectorMapType optHitDetectorSUBmap(it2_START, it2_STOP);
it2_START = it2_STOP;
if (OPSA_nPhot >= OPSA_minNrOfDetectedPhotons) { // ignore hits with too low number of detected photons
G4double OPSA_f_nPhot = OPSA_nPhot;
G4int NA = int (OPSA_fracA * OPSA_f_nPhot + 0.5);
G4int NB = int (OPSA_fracB * OPSA_f_nPhot + 0.5);
G4int NC = int (OPSA_fracC * OPSA_f_nPhot + 0.5);
G4int ND = int (OPSA_fracD * OPSA_f_nPhot + 0.5);
G4int NE = int (OPSA_fracE * OPSA_f_nPhot + 0.5);
Int_t nP=0;
// Define OPSA histograms if required for this event
if (boolStoreThisOPSAhist) {
iHistNr++;
char nameHist[100]; sprintf(nameHist,"OPSAhist_%d_%d_%d",eeeventID,OPSA_detID,iHistNr);
char nameHistTitle[100]; sprintf(nameHistTitle,"OPSAhist_%d_%d_%d;time (ns);Nr of photons",eeeventID,OPSA_detID,iHistNr);
OPSAhisto = new TH1D(nameHist, nameHistTitle, OPSAhistoNbin, OPSAhistoMin, OPSAhistoMax);
poiss = new TF1("poiss",poissonf,0.,.5,2); // x in [0;300], 2
poiss->SetParameter(0,1);
poiss->SetParameter(1,1);
}
for (optHitDetectorMapType::const_iterator it3 = optHitDetectorSUBmap.begin(); it3 != optHitDetectorSUBmap.end(); it3++) {
nP++;
if (nP==1) OPSA_timeFirst = it3->first;
if (nP==NA) OPSA_timeA = it3->first;
if (nP==NB) OPSA_timeB = it3->first;
if (nP==NC) OPSA_timeC = it3->first;
if (nP==ND) OPSA_timeD = it3->first;
if (nP==NE) OPSA_timeE = it3->first;
if (nP==OPSA_nPhot) OPSA_timeLast = it3->first;
if (boolStoreThisOPSAhist) OPSAhisto->Fill((it3->first)-OPSA_timeFirst+0.00000000001);
}
signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_timeFirst,OPSA_timeA,OPSA_timeB,OPSA_timeC,
OPSA_timeD,OPSA_timeE,OPSA_timeLast);
OPSA_signal_Map.insert(std::pair<G4int,signalInfo*>(OPSA_nPhot,mySignalInfo) );
if (boolStoreThisOPSAhist) {
OPSAhisto->Fit(poiss,"Q");
OPSAhisto->Write();
// TCanvas* cc = new TCanvas();
// OPSAhisto->Draw();
// cc->Update();
}
}
OPSA_nPhot = 0;
OPSA_timeFirst = -1000;
OPSA_timeA = -1000;
OPSA_timeB = -1000;
OPSA_timeC = -1000;
OPSA_timeD = -1000;
OPSA_timeE = -1000;
OPSA_timeLast = -1000;
}
}
}
// Now delete all optHitDetectorMap*
for (optHitMapType::const_iterator it=optHitMap.begin() ; it != optHitMap.end(); it++ ) {
delete (it->second);
}
optHitMap.clear();
// Now store the information of all mySignalInfo* to rootOutputFile
G4int nn=0;
for (OPSA_signal_MapType::reverse_iterator ritOPSA = OPSA_signal_Map.rbegin(); ritOPSA != OPSA_signal_Map.rend(); ritOPSA++) {
(ritOPSA->second) -> transferDataToRoot(myRootOutput,nn);
nn++;
}
// Now delete all mySignalInfo* from OPSA_signal_Map
for (OPSA_signal_MapType::const_iterator itOPSA = OPSA_signal_Map.begin(); itOPSA != OPSA_signal_Map.end(); itOPSA++) {
delete (itOPSA->second);
}
OPSA_signal_Map.clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Reference in New Issue View Git Blame Copy Permalink