LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

23.11.2011 Kamil Sedlak

Browse Source 
Three major changes in the photon simulation:
1) Added timeC1 ... timeC20
2) Added odet_nPhotPrim[odet_n]  (number of primary photons)
3) Added all individual photons:  phot_time[nOptPhotDet]
Documentation has not been updated yet.
This commit is contained in:
sedlak 2011-11-23 15:35:49 +00:00
parent 4d069f990c
commit 455786e3fd
5 changed files with 326 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
include/musrRootOutput.hh
View File

@ -79,12 +79,14 @@ class musrRootOutput {
G4double ekVertex, G4double xVertex, G4double yVertex, G4double zVertex,
G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex, G4int particleID) ;
void SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeSecond,
void SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4int nPhotPrim, G4double timeFirst, G4double timeSecond,
G4double timeThird, G4double timeA, G4double timeB, G4double timeC, G4double timeD,
G4double timeMean, G4double timeLast, G4double timeCFD, G4double amplCFD);
void SetCFDSpecialInfo (G4int n, G4double time);
void SetTimeC1SpecialInfo (G4double* time);
void SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke, G4double x, G4double y, G4double z, G4double time,
G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) ;
@ -114,6 +116,7 @@ class musrRootOutput {
void SetTimeInM2(G4double time) {muM2Time = time/microsecond;}
void SetInitialPositronMomentum(G4ThreeVector mom) {posIniMomx=mom.x(); posIniMomy=mom.y(); posIniMomz=mom.z();}
void SetNOptPhot(G4int value) {nOptPhot=value;}
void SetPhotDetTime(G4double time);
void SetDecayTime(G4double time) {muDecayTime=time/microsecond;}
void SetNrFieldNomVal(G4int n) {nFieldNomVal = n;}
void SetFieldNomVal(G4int i, G4double value);
@ -179,6 +182,8 @@ class musrRootOutput {
static G4bool store_posIniMomY;
static G4bool store_posIniMomZ;
static G4bool store_nOptPhot;
static G4bool store_nOptPhotDet;
static G4bool store_phot_time;
static G4bool store_det_ID;
static G4bool store_det_edep;
static G4bool store_det_edep_el;
@ -218,6 +223,7 @@ class musrRootOutput {
static G4bool store_fieldIntegralBz3;
static G4bool store_odet_ID;
static G4bool store_odet_nPhot;
static G4bool store_odet_nPhotPrim;
static G4bool store_odet_timeFirst;
static G4bool store_odet_timeSecond;
static G4bool store_odet_timeThird;
@ -230,6 +236,7 @@ class musrRootOutput {
static G4bool store_odet_timeCFD;
static G4bool store_odet_amplCFD;
static G4bool store_odet_timeCFDarray;
static G4bool store_odet_timeC1;
static G4int oldEventNumberInG4EqEMFieldWithSpinFunction;
@ -270,7 +277,9 @@ class musrRootOutput {
Double_t muDecayPosX, muDecayPosY, muDecayPosZ;
Double_t muDecayTime;
Double_t posIniMomx, posIniMomy, posIniMomz;
Int_t nOptPhot;
Int_t nOptPhot, nOptPhotDet;
static const Int_t maxNOptPhotDet=10000;
Double_t phot_time[maxNOptPhotDet];
public:
static const Int_t maxNFieldnNominalValues=30;
@ -329,6 +338,7 @@ class musrRootOutput {
G4int odet_n;
G4int odet_ID[odet_nMax];
G4int odet_nPhot[odet_nMax];
G4int odet_nPhotPrim[odet_nMax];
G4double odet_timeFirst[odet_nMax];
G4double odet_timeSecond[odet_nMax];
G4double odet_timeThird[odet_nMax];
@ -405,6 +415,26 @@ class musrRootOutput {
G4double odet_timeCFD510[odet_nMax];
G4double odet_timeCFD511[odet_nMax];
G4double odet_timeCFD512[odet_nMax];
G4double odet_timeC1[odet_nMax];
G4double odet_timeC2[odet_nMax];
G4double odet_timeC3[odet_nMax];
G4double odet_timeC4[odet_nMax];
G4double odet_timeC5[odet_nMax];
G4double odet_timeC6[odet_nMax];
G4double odet_timeC7[odet_nMax];
G4double odet_timeC8[odet_nMax];
G4double odet_timeC9[odet_nMax];
G4double odet_timeC10[odet_nMax];
G4double odet_timeC11[odet_nMax];
G4double odet_timeC12[odet_nMax];
G4double odet_timeC13[odet_nMax];
G4double odet_timeC14[odet_nMax];
G4double odet_timeC15[odet_nMax];
G4double odet_timeC16[odet_nMax];
G4double odet_timeC17[odet_nMax];
G4double odet_timeC18[odet_nMax];
G4double odet_timeC19[odet_nMax];
G4double odet_timeC20[odet_nMax];
public:
static const Int_t save_nMax=1000;

65
include/musrScintSD.hh
View File

@ -31,32 +31,55 @@
class G4Step;
class G4HCofThisEvent;
class APDidClass {
public:
APDidClass(G4int xAPDcellID, G4int xAPDcellNphot) {
APDcellID = xAPDcellID;
APDcellNphot = xAPDcellNphot;
}
~APDidClass() {}
G4int GetAPDcellID() const {return APDcellID;}
G4int GetAPDcellNphot() const {return APDcellNphot;}
void IncreaseAPDcellNphot() {APDcellNphot++;}
private:
G4int APDcellID;
G4int APDcellNphot;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class signalInfo {
public:
signalInfo(G4int id, G4int nP, G4double tFirst, G4double tSecond, G4double tThird,
signalInfo(G4int id, G4int nP, G4int nPprim, G4double tFirst, G4double tSecond, G4double tThird,
G4double tA, G4double tB, G4double tC, G4double tD, G4double tE, G4double tLast,
G4double tCFD, G4double aCFD, G4double tCFDarray[1000])
G4double tCFD, G4double aCFD, G4double tCFDarray[1000], G4double tOPSAC1array[21])
{
detID=id; nPhot=nP; timeFirst=tFirst; timeSecond=tSecond; timeThird=tThird;
detID=id; nPhot=nP; nPhotPrim=nPprim; timeFirst=tFirst; timeSecond=tSecond; timeThird=tThird;
timeA=tA; timeB=tB; timeC=tC; timeD=tD; timeE=tE; timeLast=tLast;
timeCFD=tCFD; amplCFD=aCFD;
if (musrRootOutput::store_odet_timeCFDarray) {for(int i=0;i<1000;i++) {timeCFDarray[i]=tCFDarray[i];}}
if (musrRootOutput::store_odet_timeC1) {for(int i=1;i<21;i++) {timeC1array[i] =tOPSAC1array[i];}}
}
~signalInfo() {}
void transferDataToRoot(musrRootOutput* myRootOut, G4int nn) {
myRootOut->SetOPSAinfo(nn,detID,nPhot,timeFirst,timeSecond,timeThird,
myRootOut->SetOPSAinfo(nn,detID,nPhot,nPhotPrim,timeFirst,timeSecond,timeThird,
timeA,timeB,timeC,timeD,timeE,timeLast,timeCFD,amplCFD);
for (Int_t kk=0; kk<13; kk++) {
for (Int_t ll=0; ll<5; ll++) {
int index = (ll+1)*100+kk;
myRootOut -> SetCFDSpecialInfo(index,timeCFDarray[index]);
if (musrRootOutput::store_odet_timeCFDarray) {
for (Int_t kk=0; kk<13; kk++) {
for (Int_t ll=0; ll<5; ll++) {
int index = (ll+1)*100+kk;
myRootOut -> SetCFDSpecialInfo(index,timeCFDarray[index]);
}
}
}
if (musrRootOutput::store_odet_timeC1) myRootOut -> SetTimeC1SpecialInfo(timeC1array);
}
private:
G4int detID;
G4int nPhot;
G4int nPhotPrim;
G4int nPhot_abs;
G4int nPhot_refl;
G4int nPhot_other;
@ -72,6 +95,7 @@ class signalInfo {
G4double timeCFD;
G4double amplCFD;
G4double timeCFDarray[1000];
G4double timeC1array[1000];
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -135,7 +159,8 @@ class musrScintSD : public G4VSensitiveDetector
musrRootOutput* myRootOutput;
// for optical photon counting
typedef std::multimap<G4double,Int_t> optHitDetectorMapType;
// typedef std::multimap<G4double,Int_t> optHitDetectorMapType;
typedef std::multimap<G4double,APDidClass> optHitDetectorMapType;
typedef std::map<Int_t,optHitDetectorMapType*> optHitMapType;
optHitMapType optHitMap;
// for optical photon signal analysis (OPSA)
@ -175,8 +200,28 @@ class musrScintSD : public G4VSensitiveDetector
G4bool APDcrossTalkRequested;
G4double APDcrossTalkProb;
void FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time);
void FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time, G4int nTruePhe);
static const G4double OPSA_C1_threshold;
static const G4double OPSA_C2_threshold;
static const G4double OPSA_C3_threshold;
static const G4double OPSA_C4_threshold;
static const G4double OPSA_C5_threshold;
static const G4double OPSA_C6_threshold;
static const G4double OPSA_C7_threshold;
static const G4double OPSA_C8_threshold;
static const G4double OPSA_C9_threshold;
static const G4double OPSA_C10_threshold;
static const G4double OPSA_C11_threshold;
static const G4double OPSA_C12_threshold;
static const G4double OPSA_C13_threshold;
static const G4double OPSA_C14_threshold;
static const G4double OPSA_C15_threshold;
static const G4double OPSA_C16_threshold;
static const G4double OPSA_C17_threshold;
static const G4double OPSA_C18_threshold;
static const G4double OPSA_C19_threshold;
static const G4double OPSA_C20_threshold;
};

5
src/musrDetectorConstruction.cc
View File

@ -1365,6 +1365,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"posIniMomY")==0) {musrRootOutput::store_posIniMomY = false;}
if (strcmp(tmpString2,"posIniMomZ")==0) {musrRootOutput::store_posIniMomZ = false;}
if (strcmp(tmpString2,"nOptPhot")==0) {musrRootOutput::store_nOptPhot = false;}
if (strcmp(tmpString2,"nOptPhotDet")==0) {musrRootOutput::store_nOptPhotDet = false;}
if (strcmp(tmpString2,"fieldNomVal")==0) {musrRootOutput::store_fieldNomVal = false;}
if (strcmp(tmpString2,"det_ID")==0) {musrRootOutput::store_det_ID = false;}
if (strcmp(tmpString2,"det_edep")==0) {musrRootOutput::store_det_edep = false;}
@ -1398,6 +1399,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"det_VvvParticleID")==0){musrRootOutput::store_det_VvvParticleID = false;}
if (strcmp(tmpString2,"odet_ID")==0) {musrRootOutput::store_odet_ID = false;}
if (strcmp(tmpString2,"odet_nPhot")==0) {musrRootOutput::store_odet_nPhot = false;}
if (strcmp(tmpString2,"odet_nPhotPrim")==0) {musrRootOutput::store_odet_nPhotPrim = false;}
if (strcmp(tmpString2,"odet_timeFirst")==0) {musrRootOutput::store_odet_timeFirst = false;}
if (strcmp(tmpString2,"odet_timeSecond")==0) {musrRootOutput::store_odet_timeSecond = false;}
if (strcmp(tmpString2,"odet_timeThird")==0) {musrRootOutput::store_odet_timeThird = false;}
@ -1418,7 +1420,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"fieldIntegralBz2")==0){musrRootOutput::store_fieldIntegralBz2 = true;}
if (strcmp(tmpString2,"fieldIntegralBz3")==0){musrRootOutput::store_fieldIntegralBz3 = true;}
if (strcmp(tmpString2,"odet_timeCFDarray")==0){musrRootOutput::store_odet_timeCFDarray = true;}
if (strcmp(tmpString2,"odet_timeC1")==0) {musrRootOutput::store_odet_timeC1 = true;}
if (strcmp(tmpString2,"phot_time")==0) {musrRootOutput::store_phot_time = true;}
if ((musrRootOutput::store_fieldIntegralBx)||(musrRootOutput::store_fieldIntegralBy)||
(musrRootOutput::store_fieldIntegralBz)||(musrRootOutput::store_fieldIntegralBz1)||
(musrRootOutput::store_fieldIntegralBz2)||(musrRootOutput::store_fieldIntegralBz3) ) {

78
src/musrRootOutput.cc
View File

@ -114,6 +114,8 @@ G4bool musrRootOutput::store_posIniMomX = true;
G4bool musrRootOutput::store_posIniMomY = true;
G4bool musrRootOutput::store_posIniMomZ = true;
G4bool musrRootOutput::store_nOptPhot = true;
G4bool musrRootOutput::store_nOptPhotDet= true;
G4bool musrRootOutput::store_phot_time = false;
G4bool musrRootOutput::store_det_ID = true;
G4bool musrRootOutput::store_det_edep = true;
G4bool musrRootOutput::store_det_edep_el = true;
@ -153,6 +155,7 @@ G4bool musrRootOutput::store_fieldIntegralBz2 = false;
G4bool musrRootOutput::store_fieldIntegralBz3 = false;
G4bool musrRootOutput::store_odet_ID = true;
G4bool musrRootOutput::store_odet_nPhot = true;
G4bool musrRootOutput::store_odet_nPhotPrim = true;
G4bool musrRootOutput::store_odet_timeFirst = true;
G4bool musrRootOutput::store_odet_timeSecond = true;
G4bool musrRootOutput::store_odet_timeThird = true;
@ -165,6 +168,7 @@ G4bool musrRootOutput::store_odet_timeLast = true;
G4bool musrRootOutput::store_odet_timeCFD = true;
G4bool musrRootOutput::store_odet_amplCFD = true;
G4bool musrRootOutput::store_odet_timeCFDarray = false;
G4bool musrRootOutput::store_odet_timeC1 = false;
G4int musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction=-1;
@ -229,6 +233,8 @@ void musrRootOutput::BeginOfRunAction() {
if (store_posIniMomY) {rootTree->Branch("posIniMomY",&posIniMomy,"posIniMomY/D");}
if (store_posIniMomZ) {rootTree->Branch("posIniMomZ",&posIniMomz,"posIniMomZ/D");}
if (store_nOptPhot) {rootTree->Branch("nOptPhot",&nOptPhot,"nOptPhot/I");}
if (store_nOptPhotDet) {rootTree->Branch("nOptPhotDet",&nOptPhotDet,"nOptPhotDet/I");}
if (store_phot_time) {rootTree->Branch("phot_time",&phot_time,"phot_time[nOptPhotDet]/D");}
// if (store_globalTime) {rootTree->Branch("globalTime",&globalTime,"globalTime/D");}
// if (store_fieldValue) {rootTree->Branch("fieldValue",&fieldValue,"fieldValue/D");}
if (store_fieldNomVal) {
@ -291,11 +297,13 @@ void musrRootOutput::BeginOfRunAction() {
rootTree->Branch("save_polz",&save_polz,"save_polz[save_n]/D");
}
if (store_odet_ID || store_odet_nPhot || store_odet_timeFirst || store_odet_timeSecond || store_odet_timeThird || store_odet_timeA || store_odet_timeB ||
store_odet_timeC || store_odet_timeD || store_odet_timeMean || store_odet_timeLast || store_odet_timeCFD || store_odet_amplCFD)
if (store_odet_ID || store_odet_nPhot || store_odet_nPhotPrim || store_odet_timeFirst || store_odet_timeSecond ||
store_odet_timeThird || store_odet_timeA || store_odet_timeB || store_odet_timeC || store_odet_timeD ||
store_odet_timeMean || store_odet_timeLast || store_odet_timeCFD || store_odet_amplCFD || store_odet_timeC1)
{rootTree->Branch("odet_n",&odet_n,"odet_n/I");}
if (store_odet_ID) {rootTree->Branch("odet_ID",&odet_ID,"odet_ID[odet_n]/I");}
if (store_odet_nPhot) {rootTree->Branch("odet_nPhot",&odet_nPhot,"odet_nPhot[odet_n]/I");}
if (store_odet_nPhotPrim) {rootTree->Branch("odet_nPhotPrim",&odet_nPhotPrim,"odet_nPhotPrim[odet_n]/I");}
if (store_odet_timeFirst) {rootTree->Branch("odet_timeFirst",&odet_timeFirst,"odet_timeFirst[odet_n]/D");}
if (store_odet_timeSecond) {rootTree->Branch("odet_timeSecond",&odet_timeSecond,"odet_timeSecond[odet_n]/D");}
if (store_odet_timeThird) {rootTree->Branch("odet_timeThird",&odet_timeThird,"odet_timeThird[odet_n]/D");}
@ -374,6 +382,28 @@ void musrRootOutput::BeginOfRunAction() {
rootTree->Branch("odet_timeCFD511",&odet_timeCFD511,"odet_timeCFD511[odet_n]/D");
rootTree->Branch("odet_timeCFD512",&odet_timeCFD512,"odet_timeCFD512[odet_n]/D");
}
if (store_odet_timeC1) {
rootTree->Branch("odet_timeC1",&odet_timeC1,"odet_timeC1[odet_n]/D");
rootTree->Branch("odet_timeC2",&odet_timeC2,"odet_timeC2[odet_n]/D");
rootTree->Branch("odet_timeC3",&odet_timeC3,"odet_timeC3[odet_n]/D");
rootTree->Branch("odet_timeC4",&odet_timeC4,"odet_timeC4[odet_n]/D");
rootTree->Branch("odet_timeC5",&odet_timeC5,"odet_timeC5[odet_n]/D");
rootTree->Branch("odet_timeC6",&odet_timeC6,"odet_timeC6[odet_n]/D");
rootTree->Branch("odet_timeC7",&odet_timeC7,"odet_timeC7[odet_n]/D");
rootTree->Branch("odet_timeC8",&odet_timeC8,"odet_timeC8[odet_n]/D");
rootTree->Branch("odet_timeC9",&odet_timeC9,"odet_timeC9[odet_n]/D");
rootTree->Branch("odet_timeC10",&odet_timeC10,"odet_timeC10[odet_n]/D");
rootTree->Branch("odet_timeC11",&odet_timeC11,"odet_timeC11[odet_n]/D");
rootTree->Branch("odet_timeC12",&odet_timeC12,"odet_timeC12[odet_n]/D");
rootTree->Branch("odet_timeC13",&odet_timeC13,"odet_timeC13[odet_n]/D");
rootTree->Branch("odet_timeC14",&odet_timeC14,"odet_timeC14[odet_n]/D");
rootTree->Branch("odet_timeC15",&odet_timeC15,"odet_timeC15[odet_n]/D");
rootTree->Branch("odet_timeC16",&odet_timeC16,"odet_timeC16[odet_n]/D");
rootTree->Branch("odet_timeC17",&odet_timeC17,"odet_timeC17[odet_n]/D");
rootTree->Branch("odet_timeC18",&odet_timeC18,"odet_timeC18[odet_n]/D");
rootTree->Branch("odet_timeC19",&odet_timeC19,"odet_timeC19[odet_n]/D");
rootTree->Branch("odet_timeC20",&odet_timeC20,"odet_timeC20[odet_n]/D");
}
// htest1 = new TH1F("htest1","The debugging histogram 1",50,-4.,4.);
// htest2 = new TH1F("htest2","The debugging histogram 2",50,0.,3.142);
@ -455,6 +485,7 @@ void musrRootOutput::ClearAllRootVariables() {
muDecayTime=-1000;
posIniMomx=-1000;posIniMomy=-1000;posIniMomz=-1000;
nOptPhot=0;
nOptPhotDet=0;
BxIntegral = -1000; ByIntegral = -1000; BzIntegral = -1000;
BzIntegral1 = -1000; BzIntegral2 = -1000; BzIntegral3 = -1000;
det_n=0;
@ -610,7 +641,7 @@ void musrRootOutput::SetDetectorInfoVvv (G4int nDetectors,
}
void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst,
void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4int nPhotPrim, G4double timeFirst,
G4double timeSecond, G4double timeThird, G4double timeA, G4double timeB, G4double timeC,
G4double timeD, G4double timeMean, G4double timeLast, G4double timeCFD, G4double amplCFD)
{
@ -624,6 +655,7 @@ void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4
odet_n=nDetectors+1;
odet_ID[nDetectors]=ID;
odet_nPhot[nDetectors]=nPhot;
odet_nPhotPrim[nDetectors]=nPhotPrim;
odet_timeFirst[nDetectors]=timeFirst/microsecond;
odet_timeSecond[nDetectors]=timeSecond/microsecond;
odet_timeThird[nDetectors]=timeThird/microsecond;
@ -706,6 +738,32 @@ void musrRootOutput::SetCFDSpecialInfo (G4int n, G4double time) {
else if (n==512) {odet_timeCFD512[odet_n-1] = time/microsecond;}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::SetTimeC1SpecialInfo (G4double* time) {
odet_timeC1[odet_n-1] = time[1]/microsecond;
odet_timeC2[odet_n-1] = time[2]/microsecond;
odet_timeC3[odet_n-1] = time[3]/microsecond;
odet_timeC4[odet_n-1] = time[4]/microsecond;
odet_timeC5[odet_n-1] = time[5]/microsecond;
odet_timeC6[odet_n-1] = time[6]/microsecond;
odet_timeC7[odet_n-1] = time[7]/microsecond;
odet_timeC8[odet_n-1] = time[8]/microsecond;
odet_timeC9[odet_n-1] = time[9]/microsecond;
odet_timeC10[odet_n-1] = time[10]/microsecond;
odet_timeC11[odet_n-1] = time[11]/microsecond;
odet_timeC12[odet_n-1] = time[12]/microsecond;
odet_timeC13[odet_n-1] = time[13]/microsecond;
odet_timeC14[odet_n-1] = time[14]/microsecond;
odet_timeC15[odet_n-1] = time[15]/microsecond;
odet_timeC16[odet_n-1] = time[16]/microsecond;
odet_timeC17[odet_n-1] = time[17]/microsecond;
odet_timeC18[odet_n-1] = time[18]/microsecond;
odet_timeC19[odet_n-1] = time[19]/microsecond;
odet_timeC20[odet_n-1] = time[20]/microsecond;
// std::cout<<"odet_timeC1["<<odet_n-1<<"]="<<odet_timeC1[odet_n-1]<<std::endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::setRootOutputDirectoryName(char dirName[1000]) {
strcpy(rootOutputDirectoryName,dirName);
@ -713,3 +771,17 @@ void musrRootOutput::setRootOutputDirectoryName(char dirName[1000]) {
sprintf(message,"musrRootOutput.cc::setRootOutputDirectoryName: Root output file will be stored in directory %s",dirName);
musrErrorMessage::GetInstance()->musrError(INFO,message,false);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::SetPhotDetTime(G4double time) {
if (nOptPhotDet<maxNOptPhotDet) {
phot_time[nOptPhotDet] = time;
nOptPhotDet++;
}
else {
char message[200];
sprintf(message,"musrRootOutput.cc::SetPhotDetTime: number of individual photons larger than maxNOptPhotDet (=%d)",maxNOptPhotDet);
musrErrorMessage::GetInstance()->musrError(WARNING,message,false);
}
}

222
src/musrScintSD.cc
View File

@ -60,6 +60,27 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
const G4double musrScintSD::OPSA_C1_threshold (0.5);
const G4double musrScintSD::OPSA_C2_threshold (1.);
const G4double musrScintSD::OPSA_C3_threshold (2.);
const G4double musrScintSD::OPSA_C4_threshold (3.);
const G4double musrScintSD::OPSA_C5_threshold (4.);
const G4double musrScintSD::OPSA_C6_threshold (5.);
const G4double musrScintSD::OPSA_C7_threshold (6.);
const G4double musrScintSD::OPSA_C8_threshold (7.);
const G4double musrScintSD::OPSA_C9_threshold (8.);
const G4double musrScintSD::OPSA_C10_threshold (9.);
const G4double musrScintSD::OPSA_C11_threshold (10.);
const G4double musrScintSD::OPSA_C12_threshold (11.);
const G4double musrScintSD::OPSA_C13_threshold (12.);
const G4double musrScintSD::OPSA_C14_threshold (13.);
const G4double musrScintSD::OPSA_C15_threshold (14.);
const G4double musrScintSD::OPSA_C16_threshold (15.);
const G4double musrScintSD::OPSA_C17_threshold (16.);
const G4double musrScintSD::OPSA_C18_threshold (17.);
const G4double musrScintSD::OPSA_C19_threshold (18.);
const G4double musrScintSD::OPSA_C20_threshold (19.);
musrScintSD::musrScintSD(G4String name)
:G4VSensitiveDetector(name)
{
@ -261,15 +282,6 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
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
// Do APD cell variation time if requested. Even if not requested, generate the random numbers in order to
// have reproducible simulation.
G4double APDcellsTimeVariation = G4RandGauss::shoot(0,APDcellsTimeVariationSigma);
@ -278,49 +290,7 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
// If the simulation of the effect of the finite number of cells in APD is requested, find out
// whether a given cell already fired - if so, do not store the photon.
G4int APDcellID = (APDcellsEffectRequested) ? FindAPDcellID(aStep) : 0;
FireAPDcell(optHitDetectorMap,APDcellID,tmpTime);
// if (APDcellsEffectRequested) {
// APDcellID = FindAPDcellID(aStep);
// // G4cout<<"Cell ID="<<APDcellID<<G4endl;
// for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
// if ((it2->second)==APDcellID) {
// // G4cout<<"Already fired cell ="<<APDcellID<<G4endl;
// APDcellAlreadyFired = true;
// // Now lets check whether the new photon came earlier than the already saved one.
// // If so, delete the previous photon and store the time of the new one instead:
// if (tmpTime<(it2->first)) {
// optHitDetectorMap->erase(it2);
// optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
// }
// break;
// }
// }
// }
// if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(tmpTime,APDcellID));
//
// G4cout<<" tmpTime="<<tmpTime<<G4endl;
// 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;
// }
FireAPDcell(optHitDetectorMap,APDcellID,tmpTime,1);
}
}
@ -615,14 +585,14 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
G4double rand = G4UniformRand();
if (APDcrossTalkRequested) {
if (rand<APDcrossTalkProb) {
G4int APDcellID = it2->second;
G4int APDcellID = (it2->second).GetAPDcellID();
G4double tmpTime = it2->first;
if (!APDcellsEffectRequested) FireAPDcell(optHitDetectorMap,APDcellID,tmpTime);
if (!APDcellsEffectRequested) FireAPDcell(optHitDetectorMap,APDcellID,tmpTime,0);
else {
if (rand<(0.5*APDcrossTalkProb)) APDcellID--;
else APDcellID++;
if ((APDcellID<0)||(APDcellID>=APDcell_nx*APDcell_ny*APDcell_nz)) continue; // the cross-talk cell is outside range
FireAPDcell(optHitDetectorMap,APDcellID,tmpTime);
FireAPDcell(optHitDetectorMap,APDcellID,tmpTime,0);
}
}
}
@ -633,6 +603,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
optHitDetectorMapType::const_iterator it2_LAST = optHitDetectorMap->end(); it2_LAST--;
Double_t time = -1000, lastTime = -1000;
G4int OPSA_nPhot = 0;
G4int OPSA_nPhotPrim = 0;
G4double OPSA_timeFirst = -1000000;
G4double OPSA_timeSecond = -1000000;
G4double OPSA_timeThird = -1000000;
@ -710,6 +681,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
for (optHitDetectorMapType::const_iterator it3 = optHitDetectorSUBmap.begin(); it3 != optHitDetectorSUBmap.end(); it3++) {
nP++;
G4double timePhot = it3->first;
OPSA_nPhotPrim += (it3->second).GetAPDcellNphot();
// if (APDcellsTimeVariationRequested) timePhot += G4RandGauss::shoot(0,APDcellsTimeVariationSigma); // Shifted above
if (nP==1) OPSA_timeFirst = timePhot;
if (nP==2) OPSA_timeSecond= timePhot;
@ -743,6 +715,7 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
// if required, convert the histogram with photon times to histogram of electronic pulse shape
G4double timeCFDarray[1000];
G4double OPSA_timeC1[21];
if (bool_pulseShapeExists) {
// for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of photon histogram
// Double_t photonCount = OPSAhisto ->GetBinContent(iBin);
@ -813,6 +786,120 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
if (OPSA_C_time_found && OPSA_D_time_found) break; // times C and D found, no need to continue looping.
}
if (musrRootOutput::store_odet_timeC1) {
Double_t yValue=-1000;
for (Int_t i=1; i<=20; i++) {OPSA_timeC1[i] = -1000000.;}
G4bool OPSA_C1_time_found = false; G4bool OPSA_C2_time_found = false;
G4bool OPSA_C3_time_found = false; G4bool OPSA_C4_time_found = false;
G4bool OPSA_C5_time_found = false; G4bool OPSA_C6_time_found = false;
G4bool OPSA_C7_time_found = false; G4bool OPSA_C8_time_found = false;
G4bool OPSA_C9_time_found = false; G4bool OPSA_C10_time_found = false;
G4bool OPSA_C11_time_found = false; G4bool OPSA_C12_time_found = false;
G4bool OPSA_C13_time_found = false; G4bool OPSA_C14_time_found = false;
G4bool OPSA_C15_time_found = false; G4bool OPSA_C16_time_found = false;
G4bool OPSA_C17_time_found = false; G4bool OPSA_C18_time_found = false;
G4bool OPSA_C19_time_found = false; G4bool OPSA_C20_time_found = false;
Double_t OPSAshape_histoMax = OPSAshape->GetMaximum();
for (Int_t iBin=1; iBin<=OPSAhistoNbin; iBin++) { // loop over all bins of CFD histogram
Double_t xValue = (iBin-0.5)*OPSAhistoBinWidth;
Double_t oldYvalue = yValue;
yValue = OPSAshape->GetBinContent(iBin);
if ( (yValue>=OPSA_C1_threshold) && (!OPSA_C1_time_found) ) { // signal just moved above the threshold
OPSA_C1_time_found = true;
OPSA_timeC1[1] = OPSA_timeFirst + xValue - (yValue-OPSA_C1_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C2_threshold) && (!OPSA_C2_time_found) ) { // signal just moved above the threshold
OPSA_C2_time_found = true;
OPSA_timeC1[2] = OPSA_timeFirst + xValue - (yValue-OPSA_C2_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C3_threshold) && (!OPSA_C3_time_found) ) { // signal just moved above the threshold
OPSA_C3_time_found = true;
OPSA_timeC1[3] = OPSA_timeFirst + xValue - (yValue-OPSA_C3_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C4_threshold) && (!OPSA_C4_time_found) ) { // signal just moved above the threshold
OPSA_C4_time_found = true;
OPSA_timeC1[4] = OPSA_timeFirst + xValue - (yValue-OPSA_C4_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C5_threshold) && (!OPSA_C5_time_found) ) { // signal just moved above the threshold
OPSA_C5_time_found = true;
OPSA_timeC1[5] = OPSA_timeFirst + xValue - (yValue-OPSA_C5_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C6_threshold) && (!OPSA_C6_time_found) ) { // signal just moved above the threshold
OPSA_C6_time_found = true;
OPSA_timeC1[6] = OPSA_timeFirst + xValue - (yValue-OPSA_C6_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C7_threshold) && (!OPSA_C7_time_found) ) { // signal just moved above the threshold
OPSA_C7_time_found = true;
OPSA_timeC1[7] = OPSA_timeFirst + xValue - (yValue-OPSA_C7_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C8_threshold) && (!OPSA_C8_time_found) ) { // signal just moved above the threshold
OPSA_C8_time_found = true;
OPSA_timeC1[8] = OPSA_timeFirst + xValue - (yValue-OPSA_C8_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C9_threshold) && (!OPSA_C9_time_found) ) { // signal just moved above the threshold
OPSA_C9_time_found = true;
OPSA_timeC1[9] = OPSA_timeFirst + xValue - (yValue-OPSA_C9_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C10_threshold) && (!OPSA_C10_time_found) ) { // signal just moved above the threshold
OPSA_C10_time_found = true;
OPSA_timeC1[10] = OPSA_timeFirst + xValue - (yValue-OPSA_C10_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C11_threshold) && (!OPSA_C11_time_found) ) { // signal just moved above the threshold
OPSA_C11_time_found = true;
OPSA_timeC1[11] = OPSA_timeFirst + xValue - (yValue-OPSA_C11_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C12_threshold) && (!OPSA_C12_time_found) ) { // signal just moved above the threshold
OPSA_C12_time_found = true;
OPSA_timeC1[12] = OPSA_timeFirst + xValue - (yValue-OPSA_C12_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C13_threshold) && (!OPSA_C13_time_found) ) { // signal just moved above the threshold
OPSA_C13_time_found = true;
OPSA_timeC1[13] = OPSA_timeFirst + xValue - (yValue-OPSA_C13_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C14_threshold) && (!OPSA_C14_time_found) ) { // signal just moved above the threshold
OPSA_C14_time_found = true;
OPSA_timeC1[14] = OPSA_timeFirst + xValue - (yValue-OPSA_C14_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C15_threshold) && (!OPSA_C15_time_found) ) { // signal just moved above the threshold
OPSA_C15_time_found = true;
OPSA_timeC1[15] = OPSA_timeFirst + xValue - (yValue-OPSA_C15_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C16_threshold) && (!OPSA_C16_time_found) ) { // signal just moved above the threshold
OPSA_C16_time_found = true;
OPSA_timeC1[16] = OPSA_timeFirst + xValue - (yValue-OPSA_C16_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C17_threshold) && (!OPSA_C17_time_found) ) { // signal just moved above the threshold
OPSA_C17_time_found = true;
OPSA_timeC1[17] = OPSA_timeFirst + xValue - (yValue-OPSA_C17_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C18_threshold) && (!OPSA_C18_time_found) ) { // signal just moved above the threshold
OPSA_C18_time_found = true;
OPSA_timeC1[18] = OPSA_timeFirst + xValue - (yValue-OPSA_C18_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C19_threshold) && (!OPSA_C19_time_found) ) { // signal just moved above the threshold
OPSA_C19_time_found = true;
OPSA_timeC1[19] = OPSA_timeFirst + xValue - (yValue-OPSA_C19_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if ( (yValue>=OPSA_C20_threshold) && (!OPSA_C20_time_found) ) { // signal just moved above the threshold
OPSA_C20_time_found = true;
OPSA_timeC1[20] = OPSA_timeFirst + xValue - (yValue-OPSA_C20_threshold)/(yValue-oldYvalue)*OPSAhistoBinWidth; //linear interpolation
}
if (OPSA_C20_time_found) {
if (OPSA_C1_time_found && OPSA_C2_time_found && OPSA_C3_time_found && OPSA_C4_time_found &&
OPSA_C5_time_found && OPSA_C6_time_found && OPSA_C7_time_found && OPSA_C8_time_found &&
OPSA_C9_time_found && OPSA_C10_time_found && OPSA_C11_time_found && OPSA_C12_time_found &&
OPSA_C13_time_found && OPSA_C14_time_found && OPSA_C15_time_found && OPSA_C16_time_found &&
OPSA_C17_time_found && OPSA_C18_time_found && OPSA_C19_time_found) break; // all times found
}
if (yValue==OPSAshape_histoMax) break; // no need for further searching, because maximum was achieved.
}
// myRootOutput -> SetTimeC1SpecialInfo(OPSA_timeC1);
}
// G4cout<<" OPSA_CFD_time = "<< OPSA_CFD_time<<", OPSA_CFD_ampl="<<OPSA_CFD_ampl
// <<", OPSA_timeFirst="<<OPSA_timeFirst<<", OPSA_timeC="<<OPSA_timeC <<", OPSA_timeD="<<OPSA_timeD<<G4endl;
} // end if (pulseShapeExists)
@ -826,12 +913,13 @@ void musrScintSD::EndOfEvent_OptiacalPhotons() {
}
}
signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_timeFirst,OPSA_timeSecond,OPSA_timeThird,
signalInfo* mySignalInfo = new signalInfo(OPSA_detID,OPSA_nPhot,OPSA_nPhotPrim,OPSA_timeFirst,OPSA_timeSecond,OPSA_timeThird,
OPSA_timeA,OPSA_timeB,OPSA_timeC,OPSA_timeD,OPSA_timeMean,OPSA_timeLast,
OPSA_CFD_time,OPSA_CFD_ampl,timeCFDarray);
OPSA_CFD_time,OPSA_CFD_ampl,timeCFDarray,OPSA_timeC1);
OPSA_signal_Map.insert(std::pair<G4int,signalInfo*>(OPSA_nPhot,mySignalInfo) );
}
OPSA_nPhot = 0;
OPSA_nPhotPrim = 0;
OPSA_timeFirst = -1000000;
OPSA_timeSecond = -1000000;
OPSA_timeThird = -1000000;
@ -946,23 +1034,33 @@ G4int musrScintSD::FindAPDcellID(G4Step* aStep) {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time) {
void musrScintSD::FireAPDcell(optHitDetectorMapType* optHitDetectorMap, G4int APDcellID, G4double time, G4int nTruePhe) {
if ( (musrRootOutput::store_phot_time) && (nTruePhe>0) ) myRootOutput->SetPhotDetTime(time);
if (!APDcellsEffectRequested) {
optHitDetectorMap->insert(std::pair<G4double,G4int>(time,0));
APDidClass tmpAPDid(0,nTruePhe);
// optHitDetectorMap->insert(std::pair<G4double,G4int>(time,0));
optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid));
}
else {
G4bool APDcellAlreadyFired = false;
for (optHitDetectorMapType::iterator it2 = optHitDetectorMap->begin(); it2 != optHitDetectorMap->end(); it2++ ) {
if ((it2->second)==APDcellID) { // this cell already fired before ==> check times
if ((it2->second).GetAPDcellID()==APDcellID) { // this cell already fired before ==> check times
APDcellAlreadyFired = true;
if (time<(it2->first)) { // the new cell fired before the already saved cell ==> replace it
G4int tmpAPDcellNphot = (it2->second).GetAPDcellNphot() + nTruePhe;
APDidClass tmpAPDid(APDcellID,tmpAPDcellNphot);
optHitDetectorMap->erase(it2);
optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID));
// optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID));
optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid));
}
break;
}
}
if (!APDcellAlreadyFired) optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID));
if (!APDcellAlreadyFired) {
// optHitDetectorMap->insert(std::pair<G4double,G4int>(time,APDcellID));
APDidClass tmpAPDid(APDcellID,nTruePhe);
optHitDetectorMap->insert(std::pair<G4double,APDidClass>(time,tmpAPDid));
}
}
}