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

21.1.2011 Kamil Sedlak

Browse Source 
This version contains many changes!
1) Optical photon simulation is now possible - some work still may need to be done
   (e.g. the manual is not updated yet), but it should basically work already now.
2) Changes in the musrSimAna - correction of some bugs (mainly in the coincidence
    of coincidence and veto detectors) and some other improvements
This commit is contained in:
sedlak 2011-01-21 15:20:22 +00:00
parent dcc8c6119d
commit 4bfcc6aa29
17 changed files with 1266 additions and 169 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
doc/musrSim.tex
View File

@ -845,7 +845,7 @@ The list of variables that can be stored in the Root tree:
\item{\bf muIniMomX, muIniMomY, muIniMomZ} (Double\_t) -- initial momentum of the muon when it was generated (in MeV/c).
\item{\bf muIniPolX, muIniPolY, muIniPolZ} (Double\_t) -- initial polarisation of the muon when it was generated.
\item{\bf muDecayDetID} (Int\_t) -- ID number of the detector in which the muon stopped and decayed.
\item{\bf muDecayTime} (Double\_t) -- the time at which the muon stopped and decayed (in $\mu$s).
\item{\bf muDecayTime} (Double\_t) -- the time at which the muon decayed (in $\mu$s).
\item{\bf muDecayPosX, muDecayPosY, muDecayPosZ} (Double\_t) -- the position where the muon stopped and decayed (in mm).
\item{\bf muDecayPolX, muDecayPolY, muDecayPolZ} (Double\_t) -- polarisation of the muon when it stopped and decayed.
\item{\bf muTargetTime} (Double\_t) -- time at which the muon entered the volume whose name starts by ``target'' -- usually the sample (in $\mu$s).

8
include/musrDetectorConstruction.hh
View File

@ -29,6 +29,7 @@
#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4FieldManager.hh"
#include "G4OpticalSurface.hh"
#include <map>
class G4Tubs;
@ -47,7 +48,7 @@ class musrDetectorConstruction : public G4VUserDetectorConstruction
{
public:
musrDetectorConstruction();
musrDetectorConstruction(G4String steeringFileName);
~musrDetectorConstruction();
public:
@ -61,6 +62,7 @@ public:
void ReportProblemInStearingFile(char* myString);
G4Material* CharToMaterial(char myString[100]);
G4LogicalVolume* FindLogicalVolume(G4String LogicalVolumeName);
G4VPhysicalVolume* FindPhysicalVolume(G4String PhysicalVolumeName);
void SetColourOfLogicalVolume(G4LogicalVolume* pLogVol,char* colour);
private:
@ -74,6 +76,10 @@ private:
std::map<std::string,G4RotationMatrix*> pointerToRotationMatrix;
std::map<std::string,G4FieldManager*> pointerToField;
std::map<std::string,G4MaterialPropertiesTable*> materialPropertiesTableMap;
std::map<std::string,G4MaterialPropertiesTable*>::iterator itMPT;
private:
void DefineMaterials();
};

28
include/musrRootOutput.hh
View File

@ -79,6 +79,9 @@ 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 timeA,
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast);
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) ;
@ -208,7 +211,15 @@ class musrRootOutput {
static G4bool store_fieldIntegralBz1;
static G4bool store_fieldIntegralBz2;
static G4bool store_fieldIntegralBz3;
static G4bool store_odet_ID;
static G4bool store_odet_nPhot;
static G4bool store_odet_timeFirst;
static G4bool store_odet_timeA;
static G4bool store_odet_timeB;
static G4bool store_odet_timeC;
static G4bool store_odet_timeD;
static G4bool store_odet_timeE;
static G4bool store_odet_timeLast;
static G4int oldEventNumberInG4EqEMFieldWithSpinFunction;
@ -299,6 +310,21 @@ class musrRootOutput {
G4int det_VvvTrackID[det_nMax];
G4int det_VvvParticleID[det_nMax];
public:
static const Int_t odet_nMax=det_nMax;
private:
G4int odet_n;
G4int odet_ID[odet_nMax];
G4int odet_nPhot[odet_nMax];
G4double odet_timeFirst[odet_nMax];
G4double odet_timeA[odet_nMax];
G4double odet_timeB[odet_nMax];
G4double odet_timeC[odet_nMax];
G4double odet_timeD[odet_nMax];
G4double odet_timeE[odet_nMax];
G4double odet_timeLast[odet_nMax];
public:
static const Int_t save_nMax=1000;

71
include/musrScintSD.hh
View File

@ -31,26 +31,87 @@
class G4Step;
class G4HCofThisEvent;
class signalInfo {
public:
signalInfo(G4int id, G4int nP, G4double tFirst,
G4double tA, G4double tB, G4double tC, G4double tD, G4double tE, G4double tLast) {
detID=id; nPhot=nP; timeFirst=tFirst;
timeA=tA; timeB=tB; timeC=tC; timeD=tD; timeE=tE; timeLast=tLast;}
~signalInfo() {}
void transferDataToRoot(musrRootOutput* myRootOut, G4int nn) {
myRootOut->SetOPSAinfo(nn,detID,nPhot,timeFirst,timeA,timeB,timeC,timeD,timeE,timeLast);
}
private:
G4int detID;
G4int nPhot;
G4int nPhot_abs;
G4int nPhot_refl;
G4int nPhot_other;
G4double timeFirst;
G4double timeA;
G4double timeB;
G4double timeC;
G4double timeD;
G4double timeE;
G4double timeLast;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class musrScintSD : public G4VSensitiveDetector
{
public:
musrScintSD(G4String);
~musrScintSD();
static musrScintSD* GetInstance();
musrScintSD(G4String);
~musrScintSD();
void Initialize(G4HCofThisEvent*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
void Initialize(G4HCofThisEvent*);
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
void EndOfEvent(G4HCofThisEvent*);
// Optical Photon Signal Analysis (OPSA)
void Set_OPSA_minNrOfDetectedPhotons(G4int val) {OPSA_minNrOfDetectedPhotons=val;}
void Set_OPSA_SignalSeparationTime(G4double val) {OPSA_signalSeparationTime=val;}
void Set_OPSA_frac(G4double a, G4double b, G4double c, G4double d, G4double e)
{OPSA_fracA=a; OPSA_fracB=b; OPSA_fracC=c; OPSA_fracD=d; OPSA_fracE=e;}
void AddEventIDToMultimapOfEventIDsForOPSAhistos (G4int ev_ID, G4int detector_ID)
{bool_multimapOfEventIDsForOPSAhistosEXISTS=true; multimapOfEventIDsForOPSAhistos.insert(std::pair<G4int,G4int>(ev_ID,detector_ID));}
void SetOPSAhistoBinning(Int_t nBins, Double_t min, Double_t max) {OPSAhistoNbin=nBins; OPSAhistoMin=min; OPSAhistoMax=max;}
void ProcessOpticalPhoton(G4Step*);
void EndOfEvent_OptiacalPhotons();
private:
static musrScintSD* pointer;
musrScintHitsCollection* scintCollection;
G4bool myStoreOnlyEventsWithHits;
G4int myStoreOnlyEventsWithHitInDetID;
G4double mySignalSeparationTime;
G4bool myStoreOnlyTheFirstTimeHit;
G4bool boolIsVvvInfoRequested;
musrRootOutput* myRootOutput;
// for optical photon counting
typedef std::multimap<G4double,Int_t> optHitDetectorMapType;
typedef std::map<Int_t,optHitDetectorMapType*> optHitMapType;
optHitMapType optHitMap;
// for optical photon signal analysis (OPSA)
G4int OPSA_minNrOfDetectedPhotons;
G4double OPSA_signalSeparationTime;
G4double OPSA_fracA;
G4double OPSA_fracB;
G4double OPSA_fracC;
G4double OPSA_fracD;
G4double OPSA_fracE;
typedef std::multimap<G4int,signalInfo*> OPSA_signal_MapType;
OPSA_signal_MapType OPSA_signal_Map;
G4bool bool_multimapOfEventIDsForOPSAhistosEXISTS;
typedef std::multimap<G4int,G4int> multimapOfEventIDsForOPSAhistos_Type;
multimapOfEventIDsForOPSAhistos_Type multimapOfEventIDsForOPSAhistos;
TH1D* OPSAhisto;
Int_t OPSAhistoNbin;
Double_t OPSAhistoMin, OPSAhistoMax;
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

4
include/musrSteppingAction.hh
View File

@ -23,7 +23,6 @@
#ifndef musrSteppingAction_h
#define musrSteppingAction_h 1
#include "G4UserSteppingAction.hh"
#include "G4ProcessManager.hh"
#include "globals.hh"
@ -32,6 +31,7 @@
class G4Timer;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class musrSteppingAction : public G4UserSteppingAction
@ -55,10 +55,10 @@ class musrSteppingAction : public G4UserSteppingAction
void SetCalculationOfFieldIntegralRequested(G4bool decision) {boolCalculateFieldIntegral = decision;}
private:
musrRootOutput* myRootOutput;
G4Timer* timer;
time_t realTimeWhenThisEventStarted;
static musrSteppingAction* pointer;
musrRootOutput* myRootOutput;
G4bool muAlreadyWasInTargetInThisEvent;
G4bool muAlreadyWasInM0InThisEvent;

15
musrSim.cc
View File

@ -3,6 +3,7 @@
#include "musrPrimaryGeneratorAction.hh"
#include "musrRunAction.hh"
#include "musrEventAction.hh"
//#include "StackingAction.hh"
#include "musrSteppingAction.hh"
#include "musrSteppingVerbose.hh"
@ -11,6 +12,10 @@
#include "G4UIterminal.hh"
#include "G4UItcsh.hh"
//#include <TApplication.h>
//#include <TSystem.h>
// The following two lines are needed to cope with the problem of
// "Error in <TPluginManager::FindHandler>: Cannot find plugin handler for TVirtualStreamerInfo!
// Does $ROOTSYS/etc/plugins/TVirtualStreamerInfo exist?"
@ -53,6 +58,8 @@ int main(int argc,char** argv) {
// Create class "musrErrorMessage"
musrErrorMessage* myErrorMessage = new musrErrorMessage();
// TApplication* myapp=new TApplication("myapp",0,0);
// Create Root class for storing the output of the Geant simulation
musrRootOutput* myRootOutput = new musrRootOutput();
@ -68,13 +75,14 @@ int main(int argc,char** argv) {
// UserInitialization classes (mandatory)
musrDetectorConstruction* musrdetector = new musrDetectorConstruction;
// musrDetectorConstruction* musrdetector = new musrDetectorConstruction;
if (argc>1) {
G4int myRunNr=atoi(argv[1]); // Get the run number from the name of the
// parameter file, if it starts with a number.
if (myRunNr>0) {runManager->SetRunIDCounter(myRunNr);}
musrdetector->SetInputParameterFileName(argv[1]);
// musrdetector->SetInputParameterFileName(argv[1]);
}
musrDetectorConstruction* musrdetector = new musrDetectorConstruction(steeringFileName);
runManager->SetUserInitialization(musrdetector);
runManager->SetUserInitialization(new musrPhysicsList);
@ -89,6 +97,7 @@ int main(int argc,char** argv) {
runManager->SetUserAction(new musrPrimaryGeneratorAction(musrdetector));
runManager->SetUserAction(new musrRunAction);
runManager->SetUserAction(new musrEventAction);
// runManager->SetUserAction(new StackingAction);
runManager->SetUserAction(new musrSteppingAction);
//Initialize G4 kernel
@ -133,6 +142,8 @@ int main(int argc,char** argv) {
}
}
// myapp->Run(kTRUE);
#ifdef G4VIS_USE
delete visManager;
#endif

295
musrSimAna/musrAnalysis.cxx
View File

@ -75,12 +75,12 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
exit(1);
}
std::cout << "Configuration file \"" << charV1190FileName << "\" was opened."<<std::endl;
char line[1001];
char line[10001];
// Write out the configuration file into the output file:
std::cout << "\n\n....oooOO0OOooo........oooOO0OOooo......Configuration file used for this run....oooOO0OOooo........oooOO0OOooo......"<<std::endl;
while (!feof(fSteeringFile)) {
fgets(line,1000,fSteeringFile);
fgets(line,10000,fSteeringFile);
// if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) // TS: Do not print comments
std::cout << line;
}
@ -120,7 +120,7 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
// }
// }
while (!feof(fSteeringFile)) {
fgets(line,1000,fSteeringFile);
fgets(line,10000,fSteeringFile);
if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')&&(line[0]!='$')&&(line[0]!='!')) {
char tmpString0[200]="Unset";
char tmpString1[200]="Unset";
@ -194,6 +194,18 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
sscanf(&line[strlen("MUSRMODE=")],"%s",&tmpString1);
musrMode = tmpString1[0];
}
else if (strcmp(tmpString0,"CLONECHANNEL")==0) {
int ichannel_orig_tmp, ichannel_new_tmp;
sscanf(&line[0],"%*s %d %d",&ichannel_orig_tmp,&ichannel_new_tmp);
bool_clonedChannelsMultimap_NotEmpty = true;
clonedChannelsMultimap.insert(std::pair<int,int>(ichannel_orig_tmp,ichannel_new_tmp));
}
else if (strcmp(tmpString0,"DEBUGEVENT")==0) {
int ieventToDebug_tmp, iLevelToDebug_tmp;
sscanf(&line[0],"%*s %d %d",&ieventToDebug_tmp,&iLevelToDebug_tmp);
bool_debugingRequired=true;
debugEventMap.insert(std::pair<int,int>(ieventToDebug_tmp,iLevelToDebug_tmp));
}
else if (strncmp(tmpString0,"musrTH",strlen("musrTH"))==0) {
// Definition of the histograms - either musrTH1D or musrTH2D
int beginningOfHistoTitle=0, endOfHistoTitle =0;
@ -355,6 +367,18 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
else if (strcmp(conditionNameTMP,"pileupEventCandidate")==0) conditionMap[iConditionTMP]=&pileupEventCandidate;
else if (strcmp(conditionNameTMP,"pileupEvent")==0) conditionMap[iConditionTMP]=&pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_det_AND_muonDecayedInSample_gen;
else if (strcmp(conditionNameTMP,"goodEvent_F_det")==0) conditionMap[iConditionTMP]=&goodEvent_F_det;
else if (strcmp(conditionNameTMP,"goodEvent_B_det")==0) conditionMap[iConditionTMP]=&goodEvent_B_det;
else if (strcmp(conditionNameTMP,"goodEvent_U_det")==0) conditionMap[iConditionTMP]=&goodEvent_U_det;
else if (strcmp(conditionNameTMP,"goodEvent_D_det")==0) conditionMap[iConditionTMP]=&goodEvent_D_det;
else if (strcmp(conditionNameTMP,"goodEvent_L_det")==0) conditionMap[iConditionTMP]=&goodEvent_L_det;
else if (strcmp(conditionNameTMP,"goodEvent_R_det")==0) conditionMap[iConditionTMP]=&goodEvent_R_det;
else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_D_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_D_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_L_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_L_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_R_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_R_det_AND_pileupEvent;
else {
std::cout<<" !!! ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx S T O P !!!"<<std::endl;
exit(1);
@ -404,6 +428,47 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT);
} while (nscan==2);
}
else if (strcmp(tmpString0,"counterGrouping")==0) {
int nscan; int N1; char NAME[100];
char *pch = line + strlen("counterGrouping");
char counterGroupName[100];
nscan = sscanf(pch,"%s",counterGroupName);
char* pch1 = strstr(pch,counterGroupName)+strlen(counterGroupName);
pch = pch1;
do {
nscan = sscanf(pch,"%d",&N1);
if (strcmp(counterGroupName,"F")==0) F_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"B")==0) B_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"U")==0) U_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"D")==0) D_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"L")==0) L_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"R")==0) R_posCounterList.push_back(N1);
else {
std::cout<<"\n\n UNKNOWN COUNTER GROUP REQUIRED !!! =====> S T O P F O R C E D"<<std::endl;
std::cout<<line<<std::endl;
exit(1);
}
std::cout<<"counterGroupName="<<counterGroupName<<" N1="<<N1<<std::endl;
nscan = sscanf(pch,"%s",NAME);
char* pch2 = strstr(pch ,NAME)+strlen(NAME);
pch=pch2;
nscan = sscanf(pch,"%d",&N1);
} while (nscan==1);
}
else if (strcmp(tmpString0,"artificiallyChangeMuDecayTime")==0) {
float min, max, mmmin, mmmax;
sscanf(&line[0],"%*s %g %g %g %g %g %g %g",&min,&max,&mmmin,&mmmax);
bool_muDecayTimeTransformation = true;
muDecayTime_t_min = min;
muDecayTime_t_max = max;
muDecayTime_Transformation_min = mmmin;
muDecayTime_Transformation_max = mmmax;
if ((muDecayTime_t_max <= muDecayTime_t_min) || (muDecayTime_Transformation_max <= muDecayTime_Transformation_min)) {
std::cout<<" ERROR! musrAnalysis: error when setting the \"artificiallyChangeMuDecayTime\" parameters! Min > Max !!!"<<std::endl;
std::cout<<" ==> S T O P "<<std::endl;
exit(1);
}
}
else if (strcmp(tmpString0,"fit")==0) {
char histoName[100]; char functionName[100]; float xMin; float xMax; float p0, p1, p2, p3, p4, p5, p6;
sscanf(&line[0],"%*s %s %s %g %g %g %g %g %g %g",histoName,functionName,&xMin,&xMax,&p0,&p1,&p2,&p3,&p4,&p5,&p6);
@ -453,7 +518,8 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
funct -> SetParameter(1,p1);
}
else if (strcmp(functionName,"rotFrameTime20")==0) {
funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) ");
// funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) ");
funct = new TF1("rotFrameTime20","[2]*cos(x*[0]+[1]) ");
funct -> SetParameter(0,p0);
funct -> SetParameter(1,p1);
funct -> SetParameter(2,p2);
@ -574,22 +640,8 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
pileupWindowBinMax = Long64_t(pileupWindowMax/tdcresolution);
dataWindowBinMin = Long64_t( dataWindowMin/tdcresolution);
dataWindowBinMax = Long64_t( dataWindowMax/tdcresolution);
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
char DetectorType = (it->second)->GetCounterType();
if (DetectorType=='M') {
(it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin);
(it->second)->SetCoincidenceTimeWindow(pileupWindowBinMin,pileupWindowBinMax);
(it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-mcoincwin,-mcoincwin,mcoincwin);
}
else if (DetectorType=='P') {
(it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin);
(it->second)->SetCoincidenceTimeWindow(dataWindowBinMin,dataWindowBinMax);
(it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-pcoincwin,-pcoincwin,pcoincwin);
}
else if (DetectorType=='V') {
(it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin);
(it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin);
}
int iChanNr = (it->second)->GetCounterNr();
for (std::multimap<int,int>::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) {
if ((*itCoinc).first == iChanNr) {
@ -606,6 +658,26 @@ void musrAnalysis::ReadInInputParameters(char* charV1190FileName) {
}
}
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
char DetectorType = (it->second)->GetCounterType();
if (DetectorType=='M') {
(it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin);
(it->second)->SetCoincidenceTimeWindow_M(pileupWindowBinMin,pileupWindowBinMax);
(it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin);
(it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin);
}
else if (DetectorType=='P') {
(it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin);
(it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax);
(it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin);
(it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin);
}
// else if (DetectorType=='V') {
// (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin);
// (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin);
// }
}
// for (int j=0; j<maxChannels; j++) {
// for (multimap<int,int>::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) {
// std::cout << " houby [" << (*itCoinc).first << ", " << (*itCoinc).second << "]" << std::endl;
@ -656,7 +728,6 @@ void musrAnalysis::CreateHistograms() {
// omega = 850.62*fieldValue; // value used in Geant ?
omega = 851.610*fieldValue; // value from the fits of the data
hInfo->Fill(1, fieldValue);
hInfo->Fill(6, runID);
hInfo->Fill(7, hGeantParameters->GetBinContent(7));
@ -731,12 +802,18 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
// Loop over several next event and preprocess them (i.e. fill
// them into the lists/maps of the class musrCounter).
if (bool_debugingRequired) {
if (debugEventMap[eventID]>0) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(before \"PreprocessEvent\"_________"<<std::endl;}
}
while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) {
Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
nextUnfilledEventTime+=deltaT;
};
fChain->GetEntry(iiiEntry); InitialiseEvent();
if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) PrintHitsInAllCounters();
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"PreprocessEvent\"_________"<<std::endl;}
}
// Loop over all interesting "moments", which are:
// 1) any hit in the muon counter
@ -751,7 +828,14 @@ void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
// for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
// (*it).second->myPrintThisCounter(eventID);
// }
if (bool_debugingRequired) {
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(before \"FillHistograms\"_________"<<std::endl;}
}
FillHistograms(iiiEntry);
if (bool_debugingRequired) {
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"FillHistograms\"_________"<<std::endl;}
}
currentTime+=timeToNextEvent*muonRateFactor;
@ -790,6 +874,14 @@ void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) {
//================================================================
void musrAnalysis::PrintHitsInAllCounters() {
std::cout<<"___________________\n";
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
(*it).second->myPrintThisCounter(eventID,0);
}
}
//================================================================
void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl;
@ -823,7 +915,10 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
// Check whether there was good hit in the Positron counter
// Long64_t dataBinMin = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMin : timeBinOfThePreviouslyProcessedHit-100000000;
// Long64_t dataBinMax = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMax : timeBinOfThePreviouslyProcessedHit+100000000;
pileup_eventID = -1001;
pileup_muDecayDetID_double = -1001;
pileup_muDecayPosZ = -1000000000;
pileup_muDecayPosR = -1000000000;
if (mCounterHitExistsForThisEventID) {
numberOfGoodMuons++;
Long64_t dataBinMin = timeBin0+dataWindowBinMin;
@ -833,7 +928,10 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
if (pCounterHitExistsForThisEventID && (kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) {
// This must be a pileup event (positron counter hit comes from the different event than the muon counter hit)
fChain->GetEntry(posEntry);
pileup_eventID = eventID;
pileup_muDecayDetID_double = muDecayDetID;
pileup_muDecayPosZ = muDecayPosZ;
pileup_muDecayPosR = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
// if (pileup_muDecayDetID_double==-1000) {
// std::cout<<"DEBUG: pileup_muDecayDetID_double==-1000, posEntry="<<posEntry<<", eventID="<<eventID<<", idetP_edep="<<idetP_edep<<", idetP="<<idetP<<", idetP_ID="<<idetP_ID<<std::endl;
// }
@ -875,7 +973,40 @@ void musrAnalysis::FillHistograms(Int_t iiiEntry) {
pileupEventCandidate = ((kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) ? true:false;
pileupEvent = pileupEventCandidate&&goodEvent_det;
goodEvent_det_AND_muonDecayedInSample_gen = goodEvent_det && muonDecayedInSample_gen;
// if (pileupEvent) std::cout<<"pileupEvent (eventID="<<eventID<<std::endl;
// posCounterList_Iterator = find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID);
// goodEvent_F_det = posCounterList_Iterator != F_posCounterList.end()
goodEvent_F_det = goodEvent_det && ( (find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID)) != F_posCounterList.end() );
goodEvent_B_det = goodEvent_det && ( (find(B_posCounterList.begin(), B_posCounterList.end(), idetP_ID)) != B_posCounterList.end() );
goodEvent_U_det = goodEvent_det && ( (find(U_posCounterList.begin(), U_posCounterList.end(), idetP_ID)) != U_posCounterList.end() );
goodEvent_D_det = goodEvent_det && ( (find(D_posCounterList.begin(), D_posCounterList.end(), idetP_ID)) != D_posCounterList.end() );
goodEvent_L_det = goodEvent_det && ( (find(L_posCounterList.begin(), L_posCounterList.end(), idetP_ID)) != L_posCounterList.end() );
goodEvent_R_det = goodEvent_det && ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() );
// std::cout<<"goodEvent_F_det="<<goodEvent_F_det<<std::endl;
if (pileupEvent&&goodEvent_F_det) {
std::cout<<" DEBUG: Pileup Event: eventID = "<<eventID<<" pileup_eventID = "<<pileup_eventID<<" det_time10 = "<<det_time10<<std::endl;
// debugEventMap.insert(std::pair<int,int>(eventID,10));
}
goodEvent_F_det_AND_pileupEvent = goodEvent_F_det && pileupEvent;
goodEvent_B_det_AND_pileupEvent = goodEvent_B_det && pileupEvent;
goodEvent_U_det_AND_pileupEvent = goodEvent_U_det && pileupEvent;
goodEvent_D_det_AND_pileupEvent = goodEvent_D_det && pileupEvent;
goodEvent_L_det_AND_pileupEvent = goodEvent_L_det && pileupEvent;
goodEvent_R_det_AND_pileupEvent = goodEvent_R_det && pileupEvent;
// if (bool_debugingRequired && muonTriggered_det) {
// std::cout<<"DEBUG: goodEvent_det: eventID="<<eventID<<std::endl;
// if (goodEvent_det) std::cout<<" ___DETECTED___"<<std::endl;
// MyPrintTree();
// MyPrintConditions();
// }
// if (pileupEvent) {
// std::cout<<"\n NEW: pileupEvent: eventID="<<eventID<<", kEntry="<<kEntry<<", posEntry="<<posEntry<< std::endl;
// std::cout<<"det_time10 = "<<det_time10<<std::endl;
// }
// Fill pileup-variables, but only if positron comes from different muon than the trigger signal
// if (mCounterHitExistsForThisEventID&&pCounterHitExistsForThisEventID)
@ -938,25 +1069,110 @@ void musrAnalysis::InitialiseEvent() {
Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
// std::cout<<"musrAnalysis::PreprocessEvent()"<<std::endl;
fChain->GetEntry(iEn); InitialiseEvent();
// Clone some channels into different one, if requested by user
// (This is usefull when e.g. user splits a signal from a veto
// and uses it in two different ways - e.g. once for vetoing
// muons, and second (with a different threshold) for validating
// a positron candidate. This is initiated by the
// keyword "CLONECHANNEL" in the *.v1190 file
if (bool_clonedChannelsMultimap_NotEmpty) {
// std::cout<<"det_n="<<det_n<<std::endl;
Int_t det_n_OLD=det_n;
for (Int_t i=0; i<det_n_OLD; i++) {
// std::cout<<" det_ID["<<i<<"]="<<det_ID[i]<<" edep="<<det_edep[i]<<std::endl;
clonedChannelsMultimapType::const_iterator it = clonedChannelsMultimap.find(det_ID[i]);
// std::cout<<" clonedChannelsMultimap[i]="<<clonedChannelsMultimap[i]<<std::endl;
if (it!=clonedChannelsMultimap.end()) {
int chNumTMP = it->first;
// std::cout<<"DEBUG: chNumTMP="<<chNumTMP<<" ";
std::pair<clonedChannelsMultimapType::const_iterator,clonedChannelsMultimapType::const_iterator> ret = clonedChannelsMultimap.equal_range(chNumTMP);
for (clonedChannelsMultimapType::const_iterator ittt=ret.first; ittt!=ret.second; ++ittt) {
// std::cout << " ittt->second=" << ittt->second;
int chNumNewTMP = ittt->second;
det_ID[det_n] = chNumNewTMP;
det_edep[det_n] = det_edep[i];
det_edep_el[det_n] = det_edep_el[i];
det_edep_pos[det_n] = det_edep_pos[i];
det_edep_gam[det_n] = det_edep_gam[i];
det_edep_mup[det_n] = det_edep_mup[i];
det_nsteps[det_n] = det_nsteps[i];
det_length[det_n] = det_length[i];
det_time_start[det_n] = det_time_start[i];
det_time_end[det_n] = det_time_end[i];
det_x[det_n] = det_x[i];
det_y[det_n] = det_y[i];
det_z[det_n] = det_z[i];
det_kine[det_n] = det_kine[i];
det_VrtxKine[det_n] = det_VrtxKine[i];
det_VrtxX[det_n] = det_VrtxX[i];
det_VrtxY[det_n] = det_VrtxY[i];
det_VrtxZ[det_n] = det_VrtxZ[i];
det_VrtxVolID[det_n] = det_VrtxVolID[i];
det_VrtxProcID[det_n] = det_VrtxProcID[i];
det_VrtxTrackID[det_n] = det_VrtxTrackID[i];
det_VrtxParticleID[det_n]= det_VrtxParticleID[i];
det_VvvKine[det_n] = det_VvvKine[i];
det_VvvX[det_n] = det_VvvX[i];
det_VvvY[det_n] = det_VvvY[i];
det_VvvZ[det_n] = det_VvvZ[i];
det_VvvVolID[det_n] = det_VvvVolID[i];
det_VvvProcID[det_n] = det_VvvProcID[i];
det_VvvTrackID[det_n] = det_VvvTrackID[i];
det_VvvParticleID[det_n] = det_VvvParticleID[i];
det_n++;
}
}
}
}
// std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<<eventID<<std::endl;
// MyPrintTree();
Double_t globTime = nextUnfilledEventTime;
for (Int_t i=0; i<det_n; i++) {
// // Int_t detID=det_ID[i];
std::map<int,musrCounter*>::iterator it;
it = allCounterMap.find(det_ID[i]);
if (it==allCounterMap.end()) {
// uncomment std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
// std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
}
else {
// Double_t omega = 851.372*fieldNomVal[0];
Double_t dPhaseShift = (omega==0) ? 0:phaseShiftMap[det_ID[i]]/omega;
//cDEL if (det_ID[i]<20) std::cout<<"phaseShift = "<<phaseShiftMap[det_ID[i]]<<" dPhaseShift="<<dPhaseShift<<" tdcresolution="<<tdcresolution<<std::endl;
Long64_t timeBin = Long64_t( (globTime+det_time_start[i] )/tdcresolution );
Long64_t timeBin2 = Long64_t( (globTime+det_time_start[i]+dPhaseShift)/tdcresolution );
//cDEL if (det_ID[i]<20) std::cout<<" timeBin = "<<timeBin<<" ("<<globTime+det_time_start[i]<<")"<<std::endl;
//cDEL if (det_ID[i]<20) std::cout<<" timeBin2 = "<<timeBin2<<" ("<<globTime+det_time_start[i]+dPhaseShift<<")"<<std::endl;
Double_t t1,t2;
if (bool_muDecayTimeTransformation) {
// THIS OPTION WAS SUPOSED TO ALLOW THE USER TO SOMEHOW (IN A TRICKY AND NOT 100% RELIABLE WAY)
// SHIFT THE POSITRON COUNTS TO A RESTRICTED TIME WINDOW, AS IF THE MUON DECAYED IN SOME RESTRICTED
// TIME INTERVAL. THIS, HOWEVER, DOES NOT WORK, IN THE WAY IT IS IMPLEMENTED HERE, BECAUSE THE
// MUON SPIN ROTATION AT REST IS IGNORED HERE - IT WOULD NEED SOME FURTHER WORK TO DO, AND
// BECOMES DIFFICULT TO INTERPRET LATER ==> WORK ON THIS WAS STOPPED.
Double_t diff = muDecayTime_t_max - muDecayTime_t_min;
Double_t ttt1 = det_time_start[i];
if ((ttt1 > muDecayTime_Transformation_min) && (ttt1 < muDecayTime_Transformation_max)) {
if (muDecayTime < muDecayTime_t_min) {
ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)+1) * diff + det_time_start[i];
}
else if (muDecayTime > muDecayTime_t_max) {
ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)) * diff + det_time_start[i];
}
}
t1 = (globTime+ttt1 )/tdcresolution;
t2 = (globTime+ttt1+dPhaseShift)/tdcresolution;
// std::cout<<"DEBUG: det_time_start[i]="<<det_time_start[i]<<" ttt1="<<ttt1<<" t1="<<t1<<" t2="<<t2<<std::endl;
}
else {
t1 = (globTime+det_time_start[i] )/tdcresolution;
t2 = (globTime+det_time_start[i]+dPhaseShift)/tdcresolution;
}
Long64_t timeBin = Long64_t(t1);
Long64_t timeBin2 = Long64_t(t2);
(*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i]);
}
}
@ -969,24 +1185,37 @@ Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
//================================================================
Bool_t musrAnalysis::PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Bool_t& doubleHit) {
if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID<<"\"PositronCounterHit\": pCounterMap.size()="<<pCounterMap.size()<<std::endl;}
}
if (pCounterMap.empty()) return false;
Bool_t positronHitFound = false;
Bool_t positronHitFound = false;
Bool_t goodPositronFound = false;
if (musrMode=='D') {
// Loop over all positron counters
for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
// std::cout<<" ===POSITRON==="<< pCounterMap.size()<<std::endl;
// Bool_t thereWasHit = (it->second)->GetNextGoodHitInThisEvent(evID,dataBinMin,dataBinMax,'P',tBin1,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
Bool_t thereWasHit = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
// Bool_t thereWasHit = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
Int_t positronQuality = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit);
// std::cout<<"000000000000 tBin1="<<tBin1<<" tBin2="<<tBin2<<std::endl;
if (doubleHit) {return false;} // There were two hits in the same positron counter
if (thereWasHit) {
if (positronHitFound) {doubleHit = true; return false;} // There were two hits in two different positron counters
if (positronQuality>0) {
if (positronHitFound) {
if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID<<"\"PositronCounterHit\": Coincidence with other positron candidate in other counter."<<std::endl;}
}
doubleHit = true;
return false;
} // There were two hits in two different positron counters
positronHitFound = true;
if (positronQuality>1) goodPositronFound = true;
}
}
}
return positronHitFound;
return goodPositronFound;
}
//================================================================

48
musrSimAna/musrAnalysis.hh
View File

@ -22,6 +22,13 @@
class musrTH;
//#include "musrSimGlobal.hh"
typedef std::map<int,int> debugEventMapType;
extern debugEventMapType debugEventMap;
extern Bool_t bool_debugingRequired;
class musrAnalysis {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
@ -199,7 +206,10 @@ public :
Double_t gen_time10;
Double_t det_time10_MINUS_gen_time10;
Double_t det_time20;
Double_t pileup_eventID;
Double_t pileup_muDecayDetID_double;
Double_t pileup_muDecayPosZ;
Double_t pileup_muDecayPosR;
typedef std::map<int, Bool_t*> conditionMapType;
@ -216,6 +226,18 @@ public :
Bool_t pileupEventCandidate;
Bool_t pileupEvent;
Bool_t goodEvent_det_AND_muonDecayedInSample_gen;
Bool_t goodEvent_F_det;
Bool_t goodEvent_B_det;
Bool_t goodEvent_U_det;
Bool_t goodEvent_D_det;
Bool_t goodEvent_L_det;
Bool_t goodEvent_R_det;
Bool_t goodEvent_F_det_AND_pileupEvent;
Bool_t goodEvent_B_det_AND_pileupEvent;
Bool_t goodEvent_U_det_AND_pileupEvent;
Bool_t goodEvent_D_det_AND_pileupEvent;
Bool_t goodEvent_L_det_AND_pileupEvent;
Bool_t goodEvent_R_det_AND_pileupEvent;
musrAnalysis(TTree *tree=0);
virtual ~musrAnalysis();
@ -234,6 +256,7 @@ public :
virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit);
// virtual void RewindAllTimeInfo(Double_t timeToRewind);
virtual void RewindAllTimeInfo(Long64_t timeBinsToRewind);
virtual void PrintHitsInAllCounters();
virtual void InitialiseEvent();
virtual Double_t PreprocessEvent(Long64_t iEn);
virtual Bool_t PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Bool_t& doubleHit);
@ -247,6 +270,10 @@ public :
TH1D* hGeantParameters;
TH1D* hInfo;
// typedef std::map<int,int> debugEventMapType;
// debugEventMapType debugEventMap;
// Bool_t bool_debugingRequired;
static const Double_t pi=3.14159265358979324;
private:
@ -290,6 +317,8 @@ public :
counterMapType allCounterMap;
Int_t testIVar1;
Double_t omega;
Bool_t bool_muDecayTimeTransformation;
Double_t muDecayTime_Transformation_min, muDecayTime_Transformation_max, muDecayTime_t_min, muDecayTime_t_max;
static const Double_t microsecond=1.;
static const Double_t nanosecond=0.001;
@ -327,6 +356,19 @@ private:
typedef std::multimap<Int_t,Int_t> humanDecayMultimapType;
humanDecayMapType humanDecayMap;
humanDecayMultimapType humanDecayMultimap;
typedef std::multimap<int,int> clonedChannelsMultimapType;
clonedChannelsMultimapType clonedChannelsMultimap;
Bool_t bool_clonedChannelsMultimap_NotEmpty;
// List of group of detectors: F,B,U,D,L,R:
std::list <Int_t> F_posCounterList;
std::list <Int_t> B_posCounterList;
std::list <Int_t> U_posCounterList;
std::list <Int_t> D_posCounterList;
std::list <Int_t> L_posCounterList;
std::list <Int_t> R_posCounterList;
// std::list <Int_t>::iterator posCounterList_Iterator;
};
#endif
@ -399,7 +441,10 @@ musrAnalysis::musrAnalysis(TTree *tree)
variableMap["det_time10"]=&det_time10;
variableMap["gen_time10"]=&gen_time10;
variableMap["det_time10_MINUS_gen_time10"]=&det_time10_MINUS_gen_time10;
variableMap["pileup_eventID"]=&pileup_eventID;
variableMap["pileup_muDecayDetID"]=&pileup_muDecayDetID_double;
variableMap["pileup_muDecayPosZ"]=&pileup_muDecayPosZ;
variableMap["pileup_muDecayPosR"]=&pileup_muDecayPosR;
variableMap["det_time20"]=&det_time20;
testIVar1=0;
@ -407,6 +452,9 @@ musrAnalysis::musrAnalysis(TTree *tree)
motherOfHumanDecayHistograms=NULL;
humanDecayPileupHistograms=NULL;
motherOfHumanDecayPileupHistograms=NULL;
bool_muDecayTimeTransformation = false;
bool_clonedChannelsMultimap_NotEmpty = false;
bool_debugingRequired = false;
// if parameter tree is not specified (or zero), connect the file
// used to generate this class and read the Tree.

135
musrSimAna/musrCounter.cxx
View File

@ -2,6 +2,10 @@
#include "musrCounter.hh"
#include "TCanvas.h"
typedef std::map<int,int> debugEventMapType;
debugEventMapType debugEventMap;
Bool_t bool_debugingRequired;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT) {
@ -17,8 +21,12 @@ musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TY
TDC_t1=0;
TDC_t2=0;
TDC_histoNrAdd=0;
coincidenceTimeWindowMin=0;
coincidenceTimeWindowMax=0;
antiCoincidenceTimeWindowMin=0;
antiCoincidenceTimeWindowMax=0;
coincidenceTimeWindowMin_M=0;
coincidenceTimeWindowMax_M=0;
coincidenceTimeWindowMin_P=0;
coincidenceTimeWindowMax_P=0;
maxCoincidenceTimeWindow=0;
strcpy(TDC_histoNameAdd,"Unset");
doubleHitN=0;
@ -74,7 +82,7 @@ void musrCounter::RemoveHitsInCounter(Long64_t timeBinLimit) {
// myPrintThisCounter();
// if (counterNr==1) {std::cout<<"ooooo1 timeBinLimit="<<timeBinLimit<<std::endl; myPrintThisCounter();}
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
// std::cout<<"musrCounter::RemoveHitsInCounter: counterNr="<<counterNr<<" hitMap.size()="<<hitMap.size()<<" maxCoincidenceTimeWindow="<<maxCoincidenceTimeWindow<<" bins of tdcresolution."<<std::endl;
// std::cout<<" musrCounter::RemoveHitsInCounter: counterNr="<<counterNr<<" timeBinLimit="<<timeBinLimit<<" maxCoincidenceTimeWindow="<<maxCoincidenceTimeWindow<<" counterTimeShift="<<counterTimeShift<<std::endl;
if ((it->first)>(timeBinLimit+maxCoincidenceTimeWindow-counterTimeShift)) return; //note that maxCoincidenceTimeWindow is usually negative number
else {
// std::cout<<" Deleting hit from counter "<<counterNr<<", time bin = "<<(it->first)<<std::endl;
@ -104,19 +112,31 @@ void musrCounter::RewindHitsInCounter(Long64_t timeBinsToRewind) {
}
//================================================================
Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, Bool_t ignoreHitsAtBinZero, Long64_t timeBinMinimum, Long64_t timeBinMaximum){
Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, char motherCounter, Bool_t ignoreHitsAtBinZero, Long64_t timeBinMinimum, Long64_t timeBinMaximum){
// timeBin ... time bin, at which the coincidence is searched
// counterTimeShiftOfRequestingCounter ... time shift (in bin units) of the counter, for which the coincidence is searched
// ignoreHitsAtBinZero ... if "true", hits at timeBin will be ignored (needed for searching of coincidence of M counter
// with other M counters or P counters with other P counters)
// "false" should be used for coincidence detectors and vetos.
if (hitMap.empty()) return false;
// Long64_t timeBinToTest = timeBin + counterTimeShift - counterTimeShiftOfRequestingCounter;
Long64_t timeBinToTest = timeBin;
// Long64_t timeBinToTest = timeBin;
Long64_t timeBinMin;
Long64_t timeBinMax;
// If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors).
// Otherwise use timeBinMinimum and timeBinMaximum (e.g.coincidence of a positron counter with other positron counters).
Long64_t timeBinMin = (timeBinMinimum==-123456789) ? timeBinToTest + coincidenceTimeWindowMin : timeBinToTest + timeBinMinimum;
Long64_t timeBinMax = (timeBinMaximum==-123456789) ? timeBinToTest + coincidenceTimeWindowMax : timeBinToTest + timeBinMaximum;
if (timeBinMinimum!=-123456789) timeBinMin = timeBin + timeBinMinimum; // time window requested through "timeBinMinimum"
else if (counterType == 'V') timeBinMin = timeBin + antiCoincidenceTimeWindowMin; // this is veto detector
else if (motherCounter=='M') timeBinMin = timeBin + coincidenceTimeWindowMin_M; // this is coinc. detector (or M counter) connected to M
else if (motherCounter=='P') timeBinMin = timeBin + coincidenceTimeWindowMin_P; // this is coinc. detector connected to P
if (timeBinMaximum!=-123456789) timeBinMax = timeBin + timeBinMaximum; // time window requested through "timeBinMinimum"
else if (counterType == 'V') timeBinMax = timeBin + antiCoincidenceTimeWindowMax; // this is veto detector
else if (motherCounter=='M') timeBinMax = timeBin + coincidenceTimeWindowMax_M; // this is coinc. detector (or M counter) connected to M
else if (motherCounter=='P') timeBinMax = timeBin + coincidenceTimeWindowMax_P; // this is coinc. detector connected to P
// Long64_t timeBinMin = (timeBinMinimum==-123456789) ? timeBin + coincidenceTimeWindowMin : timeBin + timeBinMinimum;
// Long64_t timeBinMax = (timeBinMaximum==-123456789) ? timeBin + coincidenceTimeWindowMax : timeBin + timeBinMaximum;
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
Long64_t timeBinOfCount_tmp = it->first;
if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) {
@ -157,17 +177,14 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t&
// std::cout<<"*** "<<evtID<<" eventNumber="<<eventNumber<<" canditas="<<numberOfMuonCandidates<<std::endl;
// Hit candidate was found. Now check its coincidences and vetos
for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) {
// if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) )) goto endOfThisHit; // no coincidence found ==> skip hit
if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) )) goto endOfThisHit; // no coincidence found ==> skip hit
if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') )) goto endOfThisHit; // no coincidence found ==> skip hit
}
for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) {
// if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) ) goto endOfThisHit; // coincidence with veto found ==> skip hit
if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) ) goto endOfThisHit; // coincidence with veto found ==> skip hit
if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'M') ) goto endOfThisHit; // coincidence with veto found ==> skip hit
}
// Check coincidences with other hits in the M counter
// it is expected that there is only one M counter
// if (this->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift,true) ) { // coincidence with another M-counter hit ==> skip hit
if (this->IsInCoincidence(timeBinOfCount_tmp,true) ) { // coincidence with another M-counter hit ==> skip hit
if (this->IsInCoincidence(timeBinOfCount_tmp,'M',true) ) { // coincidence with another M-counter hit ==> skip hit
doubleHitFound = true;
// std::cout<<"UJGUR double hit found"<<std::endl;
goto endOfThisHit;
@ -187,35 +204,58 @@ Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t&
}
//================================================================
Bool_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound) {
Int_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound) {
// INPUT PARAMETERS: evtID, timeBinMin
// OUTPUT PARAMETERS: timeBinOfNextGoodHit
//
// positronQuality = 0 ... no positron candidate found
// = 1 ... positron candidate found, but vetoed or not in a required coincidence with other detector
// = 2 ... good positron found
// Loop over the hits in the counter
if (hitMap.empty()) return false;
Int_t positronQuality=0;
if (bool_debugingRequired) {if (debugEventMap[evtID]>4) myPrintThisCounter(evtID);}
if (hitMap.empty()) return 0;
if (counterType!='P') {std::cout<<"\n!!! FATAL ERROR !!! musrCounter::GetNextGoodPositron: not the positron counter! ==> S T O P !!!\n"; exit(1);}
doubleHitFound = false;
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
Int_t eventNumber = (it->second)->eventIDnumber;
Long64_t timeBinOfCount_tmp = it->first;
if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) continue; // This hit is out of the data interval ==> skip it
if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) {
if (bool_debugingRequired) {
if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Hit out of data interval"<<std::endl;}
}
continue; // This hit is out of the data interval ==> skip it
}
// Hit candidate was found. Now check its coincidences and vetos
positronQuality=1;
for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) {
// if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) )) goto endOfThisHit; // no coincidence found ==> skip hit
if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) )) goto endOfThisHit; // no coincidence found ==> skip hit
if (bool_debugingRequired) {
if (debugEventMap[evtID]>4) { (itCounter->second)->myPrintThisCounter(evtID); }
}
if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') )) {
if (bool_debugingRequired) {
if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Coincidence required but not found"<<std::endl;}
}
goto endOfThisHit; // no coincidence found ==> skip hit
}
}
for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) {
// if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) ) goto endOfThisHit; // coincidence with veto found ==> skip hit
if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) ) goto endOfThisHit; // coincidence with veto found ==> skip hit
if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,'P') ) {
if (bool_debugingRequired) {
if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Coincidence vith veto detector found"<<std::endl;}
}
goto endOfThisHit; // coincidence with veto found ==> skip hit
}
}
// Check coincidences with other P counters
// Coincidences with other P counters must be checked in musrAnalysis.cxx
// if (this->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift,true) ) { // coincidence with another P-counter hit ==> skip hit
// hovno
//
// ADD HERE THE CHECK THAT THE POSITRON IS FOUND WITHIN THE DATA INTERVAL, i.e. within <timeBinMin,timeBinMax>.
// AND ALSO CHECK THAT THIS IS TRUE FOR THE COINCIDENCE WITH OTHER POSITRON COUNTERS!!!
if (this->IsInCoincidence(timeBinOfCount_tmp,true,timeBinMin,timeBinMax) ) { // coincidence with another P-counter hit ==> skip hit
if (this->IsInCoincidence(timeBinOfCount_tmp,'P',true,timeBinMin,timeBinMax) ) { // coincidence with another P-counter hit ==> skip hit
if (bool_debugingRequired) {
if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Coincidence with other positron candidate in the same counter."<<std::endl;}
}
doubleHitFound = true;
// std::cout<<"tttttttttttttttttttt doubleHitN="<<doubleHitN++<<std::endl;
goto endOfThisHit;
@ -227,50 +267,57 @@ Bool_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64
idetEdep = (it->second)->det_edep;
timeBinOfNextGoodHit = timeBinOfCount_tmp;
timeBinOfNextGoodHit_phaseShifted = (it->second) -> timeBin2;
if (bool_debugingRequired) {
if (debugEventMap[evtID]>3) {std::cout<<"DEBUGEVENT:"<<evtID<<"\"GetNextGoodPositron\": Good positron candidate found in this counter."<<std::endl;}
}
//cDEL std::cout<<"timeBinOfNextGoodHit ="<<timeBinOfNextGoodHit<<" timeBinOfNextGoodHit_phaseShifted = "<<timeBinOfNextGoodHit_phaseShifted<<std::endl;
return true;
return 2;
endOfThisHit:
continue;
;
// continue;
}
return false;
return positronQuality;
}
//================================================================
void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max) {
void musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max) {
// std::cout<<"QQQQQQQQQQQQQQQQQQQQQ koincidenceCounterMap.size()="<<koincidenceCounterMap.size()<<std::endl;
for (counterMapType::const_iterator it = koincidenceCounterMap.begin(); it!=koincidenceCounterMap.end(); ++it) {
if ( ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=0) &&
( ((it->second)->GetMaxCoincidenceTimeWindow()) !=maxCoinc) ) {
std::cout<<" !!!! ERROR SetCoincidenceTimeWindowOfAllCoincidenceDetectors : coincidenceTimeWindow set multiple times! ==> S T O P !!!"<<std::endl;
Long64_t maxCoinc_AlreadySet = ((it->second)->GetMaxCoincidenceTimeWindow());
if (maxCoinc < maxCoinc_AlreadySet) (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
if (motherCounter=='M') (it->second)->SetCoincidenceTimeWindow_M(min,max);
else if (motherCounter=='P') (it->second)->SetCoincidenceTimeWindow_P(min,max);
else {
std::cout<<"musrCounter::SetCoincidenceTimeWindowOfAllCoincidenceDetectors ERROR: Strange motherCounter "
<<motherCounter<<"\n ==> S T O P "<<std::endl;
exit(1);
}
(it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
(it->second)->SetCoincidenceTimeWindow(min,max);
}
}
//================================================================
void musrCounter::SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max) {
for (counterMapType::const_iterator it = vetoCounterMap.begin(); it!=vetoCounterMap.end(); ++it) {
if ( ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=0) &&
( ((it->second)->GetMaxCoincidenceTimeWindow()) !=maxCoinc) ) {
std::cout<<" !!!! ERROR SetCoincidenceTimeWindowOfAllVetoDetectors : coincidenceTimeWindow set multiple times! ==> S T O P !!!"<<std::endl;
exit(1);
}
(it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
(it->second)->SetCoincidenceTimeWindow(min,max);
Long64_t maxCoinc_AlreadySet = ((it->second)->GetMaxCoincidenceTimeWindow());
if (maxCoinc < maxCoinc_AlreadySet) (it->second)->SetMaxCoincidenceTimeWindow(maxCoinc);
(it->second)->SetAntiCoincidenceTimeWindow(min,max);
}
}
//================================================================
void musrCounter::myPrintThisCounter(Int_t evtID) {
void musrCounter::myPrintThisCounter(Int_t evtID, Int_t detail) {
Bool_t eventMixing=false;
std::cout<<"musrCounter::myPrintThisCounter: counterNr="<<counterNr<<": ";
if ((hitMap.begin()==hitMap.end()) && (detail<=1) ) return;
if (detail>1) std::cout<<"musrCounter::myPrintThisCounter: counterNr = "<<counterNr<<": ";
else std::cout<<" counter = "<<counterNr<<": ";
for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) {
std::cout<<"\t"<<it->first<<" "<<(it->second)->eventIDnumber<<",";
if (evtID != (it->second)->eventIDnumber) {eventMixing=true;}
}
if (eventMixing) {std::cout<<" Potential event mmmmmmmmixing";}
if (eventMixing) {std::cout<<" Potential event mixing";}
std::cout<<std::endl;
}

21
musrSimAna/musrCounter.hh
View File

@ -18,6 +18,11 @@ public:
Int_t det_id;
Double_t det_edep;
Long64_t timeBin2;
// extern double GlobalKamil;
// typedef std::map<int,int> debugEventMapType;
// extern debugEventMapType debugEventMap;
// extern Bool_t bool_debugingRequired;
};
class musrCounter {
@ -39,19 +44,21 @@ class musrCounter {
}
void SetMaxCoincidenceTimeWindow(Long64_t val) {maxCoincidenceTimeWindow=val;}
Long64_t GetMaxCoincidenceTimeWindow() {return maxCoincidenceTimeWindow;}
void SetCoincidenceTimeWindowOfAllCoincidenceDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max);
void SetCoincidenceTimeWindowOfAllCoincidenceDetectors(char motherCounter, Long64_t maxCoinc, Long64_t min, Long64_t max);
void SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max);
void SetCoincidenceTimeWindow(Long64_t min, Long64_t max) {coincidenceTimeWindowMin=min; coincidenceTimeWindowMax=max;}
void SetCoincidenceTimeWindow_M(Long64_t min, Long64_t max) {coincidenceTimeWindowMin_M=min; coincidenceTimeWindowMax_M=max;}
void SetCoincidenceTimeWindow_P(Long64_t min, Long64_t max) {coincidenceTimeWindowMin_P=min; coincidenceTimeWindowMax_P=max;}
void SetAntiCoincidenceTimeWindow(Long64_t min, Long64_t max) {antiCoincidenceTimeWindowMin=min; antiCoincidenceTimeWindowMax=max;}
void SetTDChistogram(char hName[200],int t0,int t1,int t2,int hNr,char hNameAdd[200]);
void FillTDChistogram(Double_t variable, Double_t vaha);
void DrawTDChistogram();
void FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID);
void RemoveHitsInCounter(Long64_t timeBinLimit);
void RewindHitsInCounter(Long64_t timeBinsToRewind);
Bool_t IsInCoincidence(Long64_t timeBin, Bool_t ignoreHitsAtBinZero=false, Long64_t timeBinMinimum=-123456789, Long64_t timeBinMaximum=-123456789);
Bool_t IsInCoincidence(Long64_t timeBin, char motherCounter, Bool_t ignoreHitsAtBinZero=false, Long64_t timeBinMinimum=-123456789, Long64_t timeBinMaximum=-123456789);
Bool_t GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
Bool_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
void myPrintThisCounter(Int_t evtID);
Int_t GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound);
void myPrintThisCounter(Int_t evtID, Int_t detail=2);
Long64_t GetNumberOfMuonCandidates(){return numberOfMuonCandidates;}
@ -69,8 +76,8 @@ class musrCounter {
int TDC_t0, TDC_t1, TDC_t2, TDC_histoNrAdd;
char TDC_histoNameAdd[200];
TH1D* histTDC;
Long64_t coincidenceTimeWindowMin;
Long64_t coincidenceTimeWindowMax;
Long64_t antiCoincidenceTimeWindowMin, coincidenceTimeWindowMin_M, coincidenceTimeWindowMin_P;
Long64_t antiCoincidenceTimeWindowMax, coincidenceTimeWindowMax_M, coincidenceTimeWindowMax_P;
Long64_t maxCoincidenceTimeWindow;
// typedef std::map<Long64_t,Int_t> hitMap_TYPE; // Long64_t = timeBin, Int_t=eventID
typedef std::map<Long64_t,hitInfo*> hitMap_TYPE; // Long64_t = timeBin, hitInfo = eventID and det_i

5
musrSimAna/musrTH.cxx
View File

@ -47,7 +47,8 @@ void musrTH::FillTH1D(Double_t vaha, Bool_t* cond){
for (Int_t i=0; i<musrAnalysis::nrConditions; i++) {
if (bool_rotating_reference_frame) {
// Double_t var = *variableToBeFilled_X;
if (cond[i]) histArray1D[i]->Fill(*variableToBeFilled_X,vaha*cos(2*musrAnalysis::pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase));
Double_t waha = vaha*exp((*variableToBeFilled_X)/2.19703)*cos(2*musrAnalysis::pi*rot_ref_frequency*(*variableToBeFilled_X)+rot_ref_phase);
if (cond[i]) histArray1D[i]->Fill(*variableToBeFilled_X,waha);
}
else {
if (cond[i]) histArray1D[i]->Fill(*variableToBeFilled_X,vaha);
@ -223,7 +224,7 @@ void musrTH::FitHistogramsIfRequired(Double_t omega) {
for (Int_t i=0; i<musrAnalysis::nrConditions; i++) {
// std::cout<<"fitted histogram pointer="<<histArray1D[i]<<std::endl;
histArray1D[i]->Fit(funct,"","",funct_xMin,funct_xMax);
histArray1D[i]->Fit(funct,"WW","",funct_xMin,funct_xMax);
// histArray1D[i]->Fit(funct,"LL","",funct_xMin,funct_xMax);
}
// if (strcmp(funct->GetName(),"simpleExpoPLUSconst")==0) {

378
src/musrDetectorConstruction.cc
View File

@ -64,6 +64,9 @@
#include "G4RegionStore.hh"
#include "G4ProductionCuts.hh"
//#include "G4OpticalSurface.hh"
#include "G4LogicalBorderSurface.hh"
#include "musrRootOutput.hh" //cks for storing some info in the Root output file
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
@ -71,9 +74,10 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrDetectorConstruction::musrDetectorConstruction()
:parameterFileName("Unset"), checkOverlap(true), aScintSD(0)
musrDetectorConstruction::musrDetectorConstruction(G4String steeringFileName)
:checkOverlap(true), aScintSD(0)
{
parameterFileName = steeringFileName;
DefineMaterials();
detectorMessenger = new musrDetectorMessenger(this);
}
@ -102,7 +106,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
// G4double roundingErr=0.01*mm; // 0.01mm precision is OK for displaying subtracted volumes, while 0.001mm is not
//----------------------
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
musrSteppingAction* mySteppingAction = musrSteppingAction::GetInstance();
G4VPhysicalVolume* pointerToWorldVolume=NULL;
// Read detector configuration parameters from the steering file:
@ -483,7 +489,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strstr(name,"save")!=NULL) {
if (volumeID!=0) { // do not store hits in special "save" volume if ID of this volume is 0
// (due to difficulties to distinguish between ID=0 and no save volume when using std::map)
musrSteppingAction* mySteppingAction = musrSteppingAction::GetInstance();
mySteppingAction->SetLogicalVolumeAsSpecialSaveVolume(logicName,volumeID);
musrRootOutput::GetRootInstance()->SetSpecialSaveVolumeDefined();
}
@ -562,6 +567,164 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
}
// cks: Optical Boundary (needed only when simulating optical photons)
else if (strcmp(tmpString1,"opticalSurface")==0) {
if (musrParameters::boolG4OpticalPhotons) {
char optSurfaceName[100];
char physVolName1[100];
char physVolName2[100];
char type[100];
char finish[100];
char model[100];
char materialPropertiesTableName[100]="Undefined";
sscanf(&line[0],"%*s %*s %s %s %s %s %s %s %s",optSurfaceName,physVolName1,physVolName2,type,finish,model,materialPropertiesTableName);
G4VPhysicalVolume* pPhysVol1 = FindPhysicalVolume(physVolName1);
G4VPhysicalVolume* pPhysVol2 = FindPhysicalVolume(physVolName2);
if ((pPhysVol1==NULL)||(pPhysVol2==NULL)) {
G4cout << "ERROR! musrDetectorConstruction::Construct(): Physical Volume not found!"<<G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
ReportGeometryProblem(line);
exit(1);
}
G4OpticalSurface* optSurfTMP = new G4OpticalSurface(optSurfaceName);
new G4LogicalBorderSurface(optSurfaceName,pPhysVol1,pPhysVol2,optSurfTMP);
std::map<std::string,G4SurfaceType> OpticalTypeMap;
std::map<std::string,G4OpticalSurfaceModel> OpticalModelMap;
std::map<std::string,G4OpticalSurfaceFinish> OpticalFinishMap;
OpticalTypeMap["dielectric_metal"]=dielectric_metal; // dielectric-metal interface
OpticalTypeMap["dielectric_dielectric"]=dielectric_dielectric; // dielectric-dielectric interface
OpticalTypeMap["dielectric_LUT"]=dielectric_LUT; // dielectric-Look-Up-Table interface
OpticalTypeMap["firsov"]=firsov; // for Firsov Process
OpticalTypeMap["x_ray"]=x_ray; // for x-ray mirror process
OpticalModelMap["glisur"]=glisur; // original GEANT3 model
OpticalModelMap["unified"]=unified; // UNIFIED model
OpticalModelMap["LUT"]=LUT; // Look-Up-Table model
OpticalFinishMap["polished"]=polished; // smooth perfectly polished surface
OpticalFinishMap["polishedfrontpainted"]=polishedfrontpainted; // smooth top-layer (front) paint
OpticalFinishMap["polishedbackpainted"]=polishedbackpainted; // same is 'polished' but with a back-paint
OpticalFinishMap["ground"]=ground; // rough surface
OpticalFinishMap["groundfrontpainted"]=groundfrontpainted; // rough top-layer (front) paint
OpticalFinishMap["groundbackpainted"]=groundbackpainted; // same as 'ground' but with a back-paint
OpticalFinishMap["polishedlumirrorair"]=polishedlumirrorair; // mechanically polished surface, with lumirror
OpticalFinishMap["polishedlumirrorglue"]=polishedlumirrorglue; // mechanically polished surface, with lumirror & meltmount
OpticalFinishMap["polishedair"]=polishedair; // mechanically polished surface
OpticalFinishMap["polishedteflonair"]=polishedteflonair; // mechanically polished surface, with teflon
OpticalFinishMap["polishedtioair"]=polishedtioair; // mechanically polished surface, with tio paint
OpticalFinishMap["polishedtyvekair"]=polishedtyvekair; // mechanically polished surface, with tyvek
OpticalFinishMap["polishedvm2000air"]=polishedvm2000air; // mechanically polished surface, with esr film
OpticalFinishMap["polishedvm2000glue"]=polishedvm2000glue; // mechanically polished surface, with esr film & meltmount
OpticalFinishMap["etchedlumirrorair"]=etchedlumirrorair; // chemically etched surface, with lumirror
OpticalFinishMap["etchedlumirrorglue"]=etchedlumirrorglue; // chemically etched surface, with lumirror & meltmount
OpticalFinishMap["etchedair"]=etchedair; // chemically etched surface
OpticalFinishMap["etchedteflonair"]=etchedteflonair; // chemically etched surface, with teflon
OpticalFinishMap["etchedtioair"]=etchedtioair; // chemically etched surface, with tio paint
OpticalFinishMap["etchedtyvekair"]=etchedtyvekair; // chemically etched surface, with tyvek
OpticalFinishMap["etchedvm2000air"]=etchedvm2000air; // chemically etched surface, with esr film
OpticalFinishMap["etchedvm2000glue"]=etchedvm2000glue; // chemically etched surface, with esr film & meltmount
OpticalFinishMap["groundlumirrorair"]=groundlumirrorair; // rough-cut surface, with lumirror
OpticalFinishMap["groundlumirrorglue"]=groundlumirrorglue; // rough-cut surface, with lumirror & meltmount
OpticalFinishMap["groundair"]=groundair; // rough-cut surface
OpticalFinishMap["groundteflonair"]=groundteflonair; // rough-cut surface, with teflon
OpticalFinishMap["groundtioair"]=groundtioair; // rough-cut surface, with tio paint
OpticalFinishMap["groundtyvekair"]=groundtyvekair; // rough-cut surface, with tyvek
OpticalFinishMap["groundvm2000air"]=groundvm2000air; // rough-cut surface, with esr film
OpticalFinishMap["groundvm2000glue"]=groundvm2000glue; // rough-cut surface, with esr film & meltmount
G4SurfaceType OpticalType = OpticalTypeMap[type];
G4OpticalSurfaceModel OpticalModel = OpticalModelMap[model];
G4OpticalSurfaceFinish OpticalFinish = OpticalFinishMap[finish];
if ((OpticalType==0)&&(strcmp(type,"dielectric_metal")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical type \""<<type<<"\" not found!"<<G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
exit(1);
}
if ((OpticalModel==0)&&(strcmp(model,"glisur")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical surface model \""<<model<<"\" not found!"<<G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
exit(1);
}
if ((OpticalFinish==0)&&(strcmp(finish,"polished")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical surface finish \""<<finish<<"\" not found!"<<G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
exit(1);
}
optSurfTMP->SetType(OpticalType);
optSurfTMP->SetFinish(OpticalFinish);
optSurfTMP->SetModel(OpticalModel);
// Assign the "Material properties table" if required by the user:
G4cout<<"materialPropertiesTableName="<<materialPropertiesTableName<<G4endl;
if (strcmp(materialPropertiesTableName,"Undefined")!=0) {
G4MaterialPropertiesTable* MPT_tmp=NULL;
itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
if (itMPT==materialPropertiesTableMap.end()) { // G4MaterialPropertiesTable of this name does not exist --> problem
G4cout<<"\n\n musrDetectorConstruction(): Material Properties Table \""<<materialPropertiesTableName
<<"\" should be assigned to G4OpticalSurface \""<<optSurfaceName<<"\""<<G4endl;
G4cout<<" but the table was not defined yet (by command /musr/command materialPropertiesTable )"<<G4endl;
G4cout << "S T O P F O R C E D"<<G4endl;
G4cout << line << G4endl;
exit(1);
}
else {
MPT_tmp = itMPT->second;
}
optSurfTMP->SetMaterialPropertiesTable(MPT_tmp);
G4cout<<optSurfTMP<<G4endl;
optSurfTMP->GetMaterialPropertiesTable()->DumpTable();
}
G4cout<<"Optical surface \""<<optSurfaceName<<"\" created. OpticalType="<<OpticalType
<<" OpticalFinish="<<OpticalFinish<<" OpticalModel="<<OpticalModel<<G4endl;
}
}
else if (strcmp(tmpString1,"OPSA")==0){ // optical photon signal analysis
if (musrParameters::boolG4OpticalPhotons) {
musrScintSD* myMusrScintSD = musrScintSD::GetInstance();
if (myMusrScintSD==NULL) {
sprintf(eMessage,"musrDetectorConstruction.cc::Construct(): musrScintSD::GetInstance() is NULL - no musr/ScintSD set?");
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
char varName[100];
float fVarValue;
int iVarValue;
sscanf(&line[0],"%*s %*s %s",varName);
if (strcmp(varName,"minNrOfDetectedPhotons")==0) {
sscanf(&line[0],"%*s %*s %*s %d",&iVarValue);
myMusrScintSD -> Set_OPSA_minNrOfDetectedPhotons(iVarValue);
}
else if (strcmp(varName,"signalSeparationTime")==0) {
sscanf(&line[0],"%*s %*s %*s %g",&fVarValue);
myMusrScintSD -> Set_OPSA_SignalSeparationTime(fVarValue*nanosecond);
}
else if (strcmp(varName,"photonFractions")==0) {
float a, b, c, d, e;
sscanf(&line[0],"%*s %*s %*s %g %g %g %g %g",&a, &b, &c, &d, &e);
myMusrScintSD -> Set_OPSA_frac(a,b,c,d,e);
}
else if (strcmp(varName,"eventsForOPSAhistos")==0) {
int i_eventID, i_detectorID;
sscanf(&line[0],"%*s %*s %*s %d %d",&i_eventID,&i_detectorID);
myMusrScintSD -> AddEventIDToMultimapOfEventIDsForOPSAhistos(i_eventID,i_detectorID);
}
else if (strcmp(varName,"OPSAhist")==0) {
int nBins;
float min, max;
sscanf(&line[0],"%*s %*s %*s %d %g %g",&nBins,&min,&max);
myMusrScintSD -> SetOPSAhistoBinning(nBins,min*nanosecond,max*nanosecond);
}
}
}
// cks: Implementation of the Global Field (to allow overlapping fields) based
// on the Peter Gumplinger's implementation of G4BeamLine code into Geant4.
//
@ -823,7 +986,6 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
int eventWeight;
char tmpLogVolName[100];
sscanf(&line[0],"%*s %*s %*s %s %d",tmpLogVolName,&eventWeight);
musrSteppingAction* mySteppingAction = musrSteppingAction::GetInstance();
mySteppingAction -> SetVolumeForMuonEventReweighting(G4String(tmpLogVolName),eventWeight);
}
else {
@ -979,6 +1141,15 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"det_VvvProcID")==0){musrRootOutput::store_det_VvvProcID = false;}
if (strcmp(tmpString2,"det_VvvTrackID")==0){musrRootOutput::store_det_VvvTrackID = false;}
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_timeFirst")==0) {musrRootOutput::store_odet_timeFirst = false;}
if (strcmp(tmpString2,"odet_timeA")==0) {musrRootOutput::store_odet_timeA = false;}
if (strcmp(tmpString2,"odet_timeB")==0) {musrRootOutput::store_odet_timeB = false;}
if (strcmp(tmpString2,"odet_timeC")==0) {musrRootOutput::store_odet_timeC = false;}
if (strcmp(tmpString2,"odet_timeD")==0) {musrRootOutput::store_odet_timeD = false;}
if (strcmp(tmpString2,"odet_timeE")==0) {musrRootOutput::store_odet_timeE = false;}
if (strcmp(tmpString2,"odet_timeLast")==0) {musrRootOutput::store_odet_timeLast = false;}
}
else if(strcmp(storeIt,"on")==0) {
if (strcmp(tmpString2,"fieldIntegralBx")==0){musrRootOutput::store_fieldIntegralBx = true;}
@ -991,12 +1162,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if ((musrRootOutput::store_fieldIntegralBx)||(musrRootOutput::store_fieldIntegralBy)||
(musrRootOutput::store_fieldIntegralBz)||(musrRootOutput::store_fieldIntegralBz1)||
(musrRootOutput::store_fieldIntegralBz2)||(musrRootOutput::store_fieldIntegralBz3) ) {
musrSteppingAction::GetInstance()->SetCalculationOfFieldIntegralRequested(true);
mySteppingAction -> SetCalculationOfFieldIntegralRequested(true);
}
}
}
else if (strcmp(tmpString1,"process")==0) {
else if ((strcmp(tmpString1,"process")==0)||(strcmp(tmpString1,"G4OpticalPhotons")==0)
||(strcmp(tmpString1,"materialPropertiesTable")==0)||(strcmp(tmpString1,"setMaterialPropertiesTable")==0)) {
; // processes are interpreded later in musrPhysicsList.cc
}
@ -1133,49 +1305,140 @@ void musrDetectorConstruction::DefineMaterials()
new G4Material("ArgonGas", z= 18., a= 39.95*g/mole, density= 0.00000000001*mg/cm3);
if (musrParameters::boolG4OpticalPhotons) {
G4NistManager* man = G4NistManager::Instance();
G4Material* scintik = man->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
// G4Material* scintik = G4Material::GetMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
G4cout<<"scintik="<<scintik<<G4endl;
if (scintik!=NULL) {
const G4int nEntries = 14;
G4double PhotonEnergy[nEntries] =
{ 2.695*eV, 2.75489*eV, 2.8175*eV, 2.88302*eV, // 460, 450, 440, 430 nm
2.95167*eV, 3.02366*eV, 3.09925*eV, 3.17872*eV, 3.26237*eV, // 420, 410, 400, 390, 380 nm
3.30587*eV, 3.35054*eV, 3.44361*eV, 3.542*eV, 3.64618*eV }; // 375, 370, 360, 350, 340 nm
G4double RefractiveIndex[nEntries] =
{ 1.58, 1.58, 1.58, 1.58,
1.58, 1.58, 1.58, 1.58, 1.58,
1.58, 1.58, 1.58, 1.58, 1.58 };
G4double Absorption[nEntries] =
{ 8*cm, 8*cm, 8*cm, 8*cm,
8*cm, 8*cm, 8*cm, 8*cm, 8*cm,
8*cm, 8*cm, 8*cm, 8*cm, 8*cm };
G4double ScintilFast[nEntries] =
{ 0.01, 0.07, 0.15, 0.26,
0.375, 0.52, 0.65, 0.80, 0.95,
1, 0.88, 0.44, 0.08, 0.01 };
G4double ScintilSlow[nEntries] =
{ 0.01, 0.07, 0.15, 0.26,
0.375, 0.52, 0.65, 0.80, 0.95,
1, 0.88, 0.44, 0.08, 0.01 };
G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
myMPT1->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex, nEntries);
myMPT1->AddProperty("ABSLENGTH", PhotonEnergy, Absorption, nEntries);
myMPT1->AddProperty("FASTCOMPONENT", PhotonEnergy, ScintilFast, nEntries);
myMPT1->AddProperty("SLOWCOMPONENT", PhotonEnergy, ScintilSlow, nEntries);
myMPT1->AddConstProperty("SCINTILLATIONYIELD", 8400./MeV);
myMPT1->AddConstProperty("RESOLUTIONSCALE",1.0);
myMPT1->AddConstProperty("FASTTIMECONSTANT",1.6*ns);
myMPT1->AddConstProperty("SLOWTIMECONSTANT",1.6*ns);
myMPT1->AddConstProperty("YIELDRATIO",1.0);
scintik->SetMaterialPropertiesTable(myMPT1);
scintik->GetMaterialPropertiesTable()->DumpTable();
if (musrParameters::boolG4OpticalPhotons) {
FILE *fSteeringFile=fopen(parameterFileName.c_str(),"r");
if (fSteeringFile) {
char line[5001];
while (!feof(fSteeringFile)) {
fgets(line,5000,fSteeringFile);
if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) {
char tmpString0[100]="Unset";
sscanf(&line[0],"%s",tmpString0);
if ( (strcmp(tmpString0,"/musr/ignore")!=0) &&(strcmp(tmpString0,"/musr/command")!=0) ) continue;
char tmpString1[100]="Unset",tmpString2[100]="Unset";
sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2);
if (strcmp(tmpString1,"materialPropertiesTable")==0){
std::string materialPropertiesTableName=tmpString2;
G4MaterialPropertiesTable* MPT_tmp;
itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
if (itMPT==materialPropertiesTableMap.end()) { // G4MaterialPropertiesTable of this name does not exist yet --> create it
MPT_tmp = new G4MaterialPropertiesTable();
materialPropertiesTableMap.insert ( std::pair<std::string,G4MaterialPropertiesTable*>(materialPropertiesTableName,MPT_tmp) );
}
else {
MPT_tmp = itMPT->second;
}
char propertyName[100];
int nEntries;
sscanf(&line[0],"%*s %*s %*s %s %d",propertyName,&nEntries);
std::cout<<" Optical Material Def: MPT_tmp="<<MPT_tmp<<", materialPropertiesTableName="<<materialPropertiesTableName
<<", propertyName="<<propertyName<<", nEntries="<<nEntries<<std::endl;
if (nEntries==0) { // AddConstProperty
float value;
sscanf(&line[0],"%*s %*s %*s %*s %*d %g",&value);
MPT_tmp->AddConstProperty(propertyName,value);
}
else { // AddProperty
char* pch = strstr(line,propertyName)+strlen(propertyName);
float value;
G4double photonEnergyArray[100];
G4double valueArray[100];
char dummy[100];
sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
for (int i=0; i<nEntries; i++) {
sscanf(pch,"%g",&value);
// G4cout<<" DDDDD var1="<<value<<" &pch="<<&pch<<G4endl;
photonEnergyArray[i]=value;
sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
}
for (int i=0; i<nEntries; i++) {
sscanf(pch,"%g",&value);
// G4cout<<" DDDDD var2="<<value<<" &pch="<<&pch<<G4endl;
valueArray[i]=value;
sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
}
MPT_tmp->AddProperty(propertyName,photonEnergyArray,valueArray,nEntries);
}
}
else if (strcmp(tmpString1,"setMaterialPropertiesTable")==0){
char tmpString3[100]="Unset";
sscanf(&line[0],"%*s %*s %*s %s",tmpString3);
std::string materialPropertiesTableName=tmpString2;
std::string materialName = tmpString3;
G4NistManager* man = G4NistManager::Instance();
G4Material* myMaterial = man->FindOrBuildMaterial(materialName);
G4MaterialPropertiesTable* MPT_tmp=NULL;
itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
if (itMPT==materialPropertiesTableMap.end()) { // G4MaterialPropertiesTable of this name does not exist --> problem
G4cout<<"\n\n musrDetectorConstruction::DefineMaterials(): Material Properties Table \""<<materialPropertiesTableName
<<"\" should be assigned to material \""<<materialName<<"\""<<G4endl;
G4cout<<" but the table was not defined yet (by command /musr/command materialPropertiesTable )"<<G4endl;
G4cout <<" ===> S T O P F O R C E D"<<G4endl;
G4cout<<line<<G4endl;
exit(1);
}
else {
MPT_tmp = itMPT->second;
}
myMaterial->SetMaterialPropertiesTable(MPT_tmp);
G4cout<<myMaterial<<G4endl;
myMaterial->GetMaterialPropertiesTable()->DumpTable();
}
}
}
}
G4cout<<"OK konec"<<G4endl;
exit(0);
fclose(fSteeringFile);
// G4cout<<" 1eV = "<<eV<<G4endl;
// G4cout<<" 1MeV = "<<MeV<<G4endl;
// G4cout<<" 1ns = "<<ns<<G4endl;
//
// G4NistManager* man = G4NistManager::Instance();
// G4Material* scintik = man->FindOrBuildMaterial("G4_PLASTIC_SC_VINYLTOLUENE");
// G4cout<<"scintik="<<scintik<<G4endl;
// if (scintik!=NULL) {
// const G4int nEntries = 14;
// G4double PhotonEnergy[nEntries] =
// { 2.695*eV, 2.75489*eV, 2.8175*eV, 2.88302*eV, // 460, 450, 440, 430 nm
// 2.95167*eV, 3.02366*eV, 3.09925*eV, 3.17872*eV, 3.26237*eV, // 420, 410, 400, 390, 380 nm
// 3.30587*eV, 3.35054*eV, 3.44361*eV, 3.542*eV, 3.64618*eV }; // 375, 370, 360, 350, 340 nm
// G4double RefractiveIndex[nEntries] =
// { 1.58, 1.58, 1.58, 1.58,
// 1.58, 1.58, 1.58, 1.58, 1.58,
// 1.58, 1.58, 1.58, 1.58, 1.58 };
// G4double Absorption[nEntries] =
// { 8*cm, 8*cm, 8*cm, 8*cm,
// 8*cm, 8*cm, 8*cm, 8*cm, 8*cm,
// 8*cm, 8*cm, 8*cm, 8*cm, 8*cm };
// G4double ScintilFast[nEntries] =
// { 0.01, 0.07, 0.15, 0.26,
// 0.375, 0.52, 0.65, 0.80, 0.95,
// 1, 0.88, 0.44, 0.08, 0.01 };
// G4double ScintilSlow[nEntries] =
// { 0.01, 0.07, 0.15, 0.26,
// 0.375, 0.52, 0.65, 0.80, 0.95,
// 1, 0.88, 0.44, 0.08, 0.01 };
// G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
// myMPT1->AddProperty("RINDEX", PhotonEnergy, RefractiveIndex, nEntries);
// myMPT1->AddProperty("ABSLENGTH", PhotonEnergy, Absorption, nEntries);
// myMPT1->AddProperty("FASTCOMPONENT", PhotonEnergy, ScintilFast, nEntries);
// myMPT1->AddProperty("SLOWCOMPONENT", PhotonEnergy, ScintilSlow, nEntries);
// myMPT1->AddConstProperty("SCINTILLATIONYIELD", 8400./MeV);
// myMPT1->AddConstProperty("RESOLUTIONSCALE",1.0);
// myMPT1->AddConstProperty("FASTTIMECONSTANT",1.6*ns);
// myMPT1->AddConstProperty("SLOWTIMECONSTANT",1.6*ns);
// myMPT1->AddConstProperty("YIELDRATIO",1.0);
// scintik->SetMaterialPropertiesTable(myMPT1);
// scintik->GetMaterialPropertiesTable()->DumpTable();
// }
}
}
@ -1240,6 +1503,25 @@ G4LogicalVolume* musrDetectorConstruction::FindLogicalVolume(G4String LogicalVol
return NULL;
}
G4VPhysicalVolume* musrDetectorConstruction::FindPhysicalVolume(G4String PhysicalVolumeName) {
G4PhysicalVolumeStore* pPhysStore = G4PhysicalVolumeStore::GetInstance();
if (pPhysStore==NULL) {
G4cout<<"ERROR: musrDetectorConstruction.cc: G4PhysicalVolumeStore::GetInstance() not found!"<<G4endl;
}
else {
for (unsigned int i=0; i<pPhysStore->size(); i++) {
G4VPhysicalVolume* pPhysVol=pPhysStore->at(i);
G4String iPhysName=pPhysVol->GetName();
if (iPhysName==PhysicalVolumeName) {
return pPhysVol;
}
}
}
G4cout<<"\n musrDetectorConstruction.cc::FindPhysicalVolume: \n ERROR !!! Physical volume "
<<PhysicalVolumeName<<" not found !!!"<<G4endl;
return NULL;
}
void musrDetectorConstruction::SetColourOfLogicalVolume(G4LogicalVolume* pLogVol,char* colour) {
if (pLogVol!=NULL) {
if (strcmp(colour,"red" )==0) {pLogVol->SetVisAttributes(G4Colour(1,0,0));}

2
src/musrEventAction.cc
View File

@ -78,6 +78,8 @@ void musrEventAction::EndOfEventAction(const G4Event* evt) {
// cout << ":." << flush;
long thisEventNr = (long) evt->GetEventID();
// musrSteppingAction::GetInstance()->DoAtTheEndOfEvent();
// write out the root tree:
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
G4RunManager* fRunManager = G4RunManager::GetRunManager();

44
src/musrPhysicsList.cc
View File

@ -46,6 +46,8 @@
#include "G4StepLimiter.hh"
#include "G4UserSpecialCuts.hh"
#include "G4OpticalPhysics.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrPhysicsList::musrPhysicsList(): G4VUserPhysicsList()
@ -55,6 +57,26 @@ musrPhysicsList::musrPhysicsList(): G4VUserPhysicsList()
cutForElectron = 0.1*mm;
cutForMuon = 0.01*mm;
SetVerboseLevel(0);
if (musrParameters::boolG4OpticalPhotons) {
// // * Optical Physics
// G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
// RegisterPhysics( opticalPhysics );
//
// // adjust some parameters for the optical physics
// opticalPhysics->SetWLSTimeProfile("delta");
//
// opticalPhysics->SetScintillationYieldFactor(1.0);
// opticalPhysics->SetScintillationExcitationRatio(0.0);
//
// opticalPhysics->SetMaxNumPhotonsPerStep(100);
// opticalPhysics->SetMaxBetaChangePerStep(10.0);
//
// opticalPhysics->SetTrackSecondariesFirst(true);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -71,6 +93,13 @@ void musrPhysicsList::ConstructParticle()
ConstructLeptons();
ConstructMesons();
ConstructBaryons();
if (musrParameters::boolG4OpticalPhotons) {
// optical photon
G4OpticalPhoton::OpticalPhotonDefinition();
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -234,6 +263,9 @@ void musrPhysicsList::ConstructProcess()
// For a simple Muonium "scattering" when Mu hits solid materials
#include "musrMuScatter.hh"
// For optical photons
#include "G4Scintillation.hh"
// For models
#include "G4ProcessTable.hh"
#include "G4WentzelVIModel.hh"
@ -273,6 +305,11 @@ void musrPhysicsList::ConstructEM()
G4ParticleDefinition* particleDefinition = G4ParticleTable::GetParticleTable() -> FindParticle(stringParticleName);
// G4cout<<"particleDefinition of "<<stringParticleName<<" = "<<particleDefinition<<G4endl;
if (particleDefinition==NULL) {
if ((stringParticleName=="opticalphoton")&& (!(musrParameters::boolG4OpticalPhotons))) {
musrErrorMessage::GetInstance()->musrError(INFO,
"musrPhysicsList: Ignoring optical photon process definition because G4OpticalPhotons in *.mac file is set to false",false);
continue;
}
sprintf(eMessage,"musrPhysicsList: Partile \"%s\" not found in G4ParticleTable when trying to find or assign process \"%s\".",
charParticleName,charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
@ -294,6 +331,11 @@ void musrPhysicsList::ConstructEM()
else if (stringProcessName=="G4LowEnergyRayleigh") pManager->AddDiscreteProcess(new G4LowEnergyRayleigh);
else if (stringProcessName=="G4CoulombScattering") pManager->AddDiscreteProcess(new G4CoulombScattering);
else if (stringProcessName=="G4OpAbsorption") pManager->AddDiscreteProcess(new G4OpAbsorption);
else if (stringProcessName=="G4OpRayleigh") pManager->AddDiscreteProcess(new G4OpRayleigh);
else if (stringProcessName=="G4OpBoundaryProcess") pManager->AddDiscreteProcess(new G4OpBoundaryProcess);
else if (stringProcessName=="G4OpWLS") pManager->AddDiscreteProcess(new G4OpWLS);
else {
sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess. It can be easily added.",
charProcessName);
@ -325,6 +367,7 @@ void musrPhysicsList::ConstructEM()
// else if (stringProcessName=="G4hIonisation") pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
// else if (stringProcessName=="G4hLowEnergyIonisation") pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
else if (stringProcessName=="musrMuFormation") pManager->AddProcess(new musrMuFormation,nr1,nr2,nr3);
else if (stringProcessName=="G4Scintillation") pManager->AddProcess(new G4Scintillation,nr1,nr2,nr3);
// cks: musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments
// in his original "musrPhysicsList.cc about implementing musrMuScatter.
// else if (stringProcessName=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3);
@ -394,7 +437,6 @@ void musrPhysicsList::ConstructEM()
G4MuMultipleScattering* mmm = (G4MuMultipleScattering*) processTable->FindProcess("muMsc",particleDefinition);
mmm->AddEmModel(modelPriority, new G4UrbanMscModel93());
}
else {
sprintf(eMessage,"musrPhysicsList: Model \"%s\" is not implemented for \"%s\" in musrPhysicsList.cc for addModel. It can be easily added.",
charModelName,charProcessName);

51
src/musrRootOutput.cc
View File

@ -147,6 +147,16 @@ G4bool musrRootOutput::store_fieldIntegralBz = false;
G4bool musrRootOutput::store_fieldIntegralBz1 = false;
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_timeFirst = true;
G4bool musrRootOutput::store_odet_timeA = true;
G4bool musrRootOutput::store_odet_timeB = true;
G4bool musrRootOutput::store_odet_timeC = true;
G4bool musrRootOutput::store_odet_timeD = true;
G4bool musrRootOutput::store_odet_timeE = true;
G4bool musrRootOutput::store_odet_timeLast = true;
G4int musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction=-1;
@ -266,6 +276,19 @@ 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_timeA || store_odet_timeB ||
store_odet_timeC || store_odet_timeD || store_odet_timeE || store_odet_timeLast)
{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_timeFirst) {rootTree->Branch("odet_timeFirst",&odet_timeFirst,"odet_timeFirst[odet_n]/D");}
if (store_odet_timeA) {rootTree->Branch("odet_timeA",&odet_timeA,"odet_timeA[odet_n]/D");}
if (store_odet_timeB) {rootTree->Branch("odet_timeB",&odet_timeB,"odet_timeB[odet_n]/D");}
if (store_odet_timeC) {rootTree->Branch("odet_timeC",&odet_timeC,"odet_timeC[odet_n]/D");}
if (store_odet_timeD) {rootTree->Branch("odet_timeD",&odet_timeD,"odet_timeD[odet_n]/D");}
if (store_odet_timeE) {rootTree->Branch("odet_timeE",&odet_timeE,"odet_timeE[odet_n]/D");}
if (store_odet_timeLast) {rootTree->Branch("odet_timeLast",&odet_timeLast,"odet_timeLast[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);
htest1 = new TH2F("htest1","x, y",50,-200.,200.,50,-200.,200.);
@ -317,7 +340,7 @@ void musrRootOutput::FillEvent() {
htest5->Fill(atan2(posIniMomy,posIniMomx));
htest6->Fill(atan2(sqrt(posIniMomx*posIniMomx+posIniMomy*posIniMomy),posIniMomz));
if (weight>0.) {
if ( !((musrParameters::storeOnlyEventsWithHits)&&(det_n<=0)) ) {
if ( !((musrParameters::storeOnlyEventsWithHits)&&(det_n<=0)&&(odet_n<=0)) ) {
rootTree->Fill();
}
}
@ -349,7 +372,7 @@ void musrRootOutput::ClearAllRootVariables() {
BzIntegral1 = -1000; BzIntegral2 = -1000; BzIntegral3 = -1000;
det_n=0;
save_n=0;
odet_n=0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -498,3 +521,27 @@ void musrRootOutput::SetDetectorInfoVvv (G4int nDetectors,
det_VvvParticleID[nDetectors]=particleID;
}
}
void musrRootOutput::SetOPSAinfo (G4int nDetectors, G4int ID, G4int nPhot, G4double timeFirst, G4double timeA,
G4double timeB, G4double timeC, G4double timeD, G4double timeE, G4double timeLast)
{
if ((nDetectors<0)||(nDetectors>=(odet_nMax-1))) {
char message[200];
sprintf(message,"musrRootOutput.cc::SetOPSAInfo: nDetectors %i is larger than det_nMax = %i",nDetectors,det_nMax);
musrErrorMessage::GetInstance()->musrError(SERIOUS,message,false);
return;
}
else {
odet_n=nDetectors+1;
odet_ID[nDetectors]=ID;
odet_nPhot[nDetectors]=nPhot;
odet_timeFirst[nDetectors]=timeFirst/microsecond;
odet_timeA[nDetectors]=timeA/microsecond;
odet_timeB[nDetectors]=timeB/microsecond;
odet_timeC[nDetectors]=timeC/microsecond;
odet_timeD[nDetectors]=timeD/microsecond;
odet_timeE[nDetectors]=timeE/microsecond;
odet_timeLast[nDetectors]=timeLast/microsecond;
}
}

275
src/musrScintSD.cc
View File

@ -29,9 +29,14 @@
#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); }
@ -49,17 +54,38 @@
//
//....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) {
@ -78,6 +104,7 @@ void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
myStoreOnlyEventsWithHitInDetID = musrParameters::storeOnlyEventsWithHitInDetID;
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
boolIsVvvInfoRequested = myMusrSteppingAction->IsVvvInfoRequested();
myRootOutput = musrRootOutput::GetRootInstance();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -91,7 +118,15 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
}
G4Track* aTrack = aStep->GetTrack();
G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName();
// 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) {
@ -102,16 +137,20 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
}
}
musrScintHit* newHit = new musrScintHit();
newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
// 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(hitLogicalVolume->GetName());
// 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());
@ -141,6 +180,94 @@ G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
//....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();
Int_t detID = myRootOutput->ConvertVolumeToID(actualVolume);
// G4cout<<" detID ="<<detID<<" actualVolume="<<actualVolume;
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);
}
// 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;
@ -149,8 +276,6 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
<< " hits in the scint chambers: " << G4endl;
}
// Positron_momentum_already_stored=0;
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
G4int NbHits = scintCollection->entries();
@ -367,6 +492,144 @@ void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
}
} //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......

49
src/musrSteppingAction.cc
View File

@ -29,25 +29,35 @@
#include "G4MagneticField.hh" // needed for storing the magnetic field to the Root class
#include "G4FieldManager.hh" // ---------------||------------------
#include "G4TransportationManager.hh" // ---------------||------------------
//#include "G4OpBoundaryProcess.hh" // Optical photon process
#include "musrErrorMessage.hh"
#include "musrParameters.hh"
#include "F04GlobalField.hh"
//#include "TCanvas.h"
//#include "TMath.h"
//#include "TF1.h"
//#include "G4GeometryTolerance.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//Double_t poissonf(Double_t* x, Double_t* par) {
// return par[0]*TMath::Poisson(x[0],par[1]);
//}
musrSteppingAction::musrSteppingAction() {
pointer=this;
pointer=this;
boolIsAnySpecialSaveVolumeDefined = false;
boolIsVvvInfoRequested = false;
boolMuonEventReweighting = false;
boolCalculateFieldIntegral = false;
myRootOutput = musrRootOutput::GetRootInstance();
if (myRootOutput == NULL) {
musrErrorMessage::GetInstance()->musrError(FATAL,
"musrSteppingAction::musrSteppingAction(): pointer to the musrRootOutput class not found! ==> EXECUTION STOPPED",true);
}
lastActualVolume="Unset";
boolIsAnySpecialSaveVolumeDefined = false;
boolIsVvvInfoRequested = false;
boolMuonEventReweighting = false;
boolCalculateFieldIntegral = false;
myRootOutput = musrRootOutput::GetRootInstance();
if (myRootOutput == NULL) {
musrErrorMessage::GetInstance()->musrError(FATAL,
"musrSteppingAction::musrSteppingAction(): pointer to the musrRootOutput class not found! ==> EXECUTION STOPPED",true);
}
lastActualVolume="Unset";
}
musrSteppingAction::~musrSteppingAction() {
@ -140,7 +150,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
}
// Save info about the old tracks, if the user wish to have Vvv info in the output Root Tree.
// Save info about the old tracks, if the user wishes to have Vvv info in the output Root Tree.
if (boolIsVvvInfoRequested) {
G4VPhysicalVolume* nextVolume = aTrack->GetNextVolume();
if (nextVolume!=NULL) {
@ -375,6 +385,7 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ -472,3 +483,17 @@ G4bool musrSteppingAction::AreTracksCommingFromSameParent(G4int trackID1, G4int
// G4cout<<"\t\t\t\ttrack1="<<track1<<"\ttrack2="<<track2<<G4endl;
return false;
}
// //Double_t musrSteppingAction::poissonf(Double_t* x, Double_t* par)
// //{
//Double_t musrSteppingAction::poissonf(Double_t* x, Double_t* par) {
// // if (x<0)
// // return 0;
// // else if (x == 0.0)
// // return 1./Math::Exp(par);
// // else {
// // Double_t lnpoisson = x*log(par)-par-LnGamma(x+1.);
// // return Exp(lnpoisson);
// // }
// return par[0]*TMath::Poisson(x[0],par[1]);
//}