-
Notifications
You must be signed in to change notification settings - Fork 5
/
readPhotonHits.cc
153 lines (152 loc) · 6.77 KB
/
readPhotonHits.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//---------------------------------------------------------------------
//* |\___/| *
//* ) ( *
//* =\ /= *
//* )===( *
//* / \ CaTS: Calorimeter and Tracker Simulation *
//* | | is a flexible and extend-able framework *
//* / \ for the simulation of various detector *
//* \ / systems *
//* \__ _/ https://github.com/hanswenzel/CaTS *
//* ( ( *
//* ) ) *
//* (_( *
//* CaTS also serves as an example that demonstrates how to use *
//* opticks from within Geant4 for the creation and propagation of *
//* optical photons. *
//* see https://bitbucket.org/simoncblyth/opticks.git). *
//* Ascii Art by Joan Stark: https://www.asciiworld.com/-Cats-2-.html *
//---------------------------------------------------------------------
//
/// \file readPhotonHits.cc
/// \brief example how to read the CaTS::PhotonHits
//
// Root headers
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TSystem.h"
#include "TTree.h"
// Project headers
#include "Event.hh"
#include "PhotonHit.hh"
#include "lArTPCHit.hh"
int main(int argc, char** argv)
{
// initialize ROOT
TSystem ts;
gSystem->Load("libCaTSClassesDict");
if(argc < 4)
{
G4cout << "Program requires 3 arguments: name of input file, name of "
"output file, Volume that sensitive detector is attached to"
<< G4endl;
exit(1);
}
TFile* outfile = new TFile(argv[2], "RECREATE");
outfile->cd();
TH2F* pos2 = new TH2F("position", "position of Photon Hits", 400, -1000., 1000., 400, -500, 500);
TH1F* pos0 = new TH1F("position0", "position of Photon Hits", 400, -100., 100.);
TH1F* time = new TH1F("time", "timing of photon hits", 1000, 0., 2500);
TH1F* time0 = new TH1F("time0", "timing of photon hits", 1000, 0., 2500.);
TH1F* time1 = new TH1F("time1", "timing of photon hits", 1000, 0., 2500.);
TH1F* time2 = new TH1F("time2", "timing of photon hits", 1000, 0., 2500.);
TH1F* time3 = new TH1F("time3", "timing of photon hits", 1000, 0., 2500.);
TH1F* time4 = new TH1F("time4", "timing of photon hits", 1000, 0., 2500.);
TH1F* wl = new TH1F("wl", "wavelength of detected photons", 1000, 0., 1000.);
TH1F* wlce = new TH1F("wlce", "wavelength of detected Cerenkov photons", 1000, 0., 1000.);
TH1F* wlsc = new TH1F("wlsc", "wavelength of detected Scintillation photons", 1000, 0., 1000.);
TH1F* np = new TH1F("np", "number of detected photons", 100, 0., 50.);
TFile fo(argv[1]);
fo.GetListOfKeys()->Print();
Event* event = new Event();
TTree* Tevt = (TTree*) fo.Get("Events");
Tevt->SetBranchAddress("event.", &event);
TBranch* fevtbranch = Tevt->GetBranch("event.");
Int_t nevent = fevtbranch->GetEntries();
G4cout << "Nr. of Events: " << nevent << G4endl;
std::string CollectionName = argv[3];
CollectionName = CollectionName + "_Photondetector_HC";
for(Int_t i = 0; i < nevent; i++)
{
fevtbranch->GetEntry(i);
auto* hcmap = event->GetHCMap();
for(const auto& ele : *hcmap)
{
G4cout << ele.first << " " << CollectionName << G4endl;
auto hits = ele.second;
if(ele.first.compare(CollectionName) == 0)
{
auto hits = ele.second;
G4int NbHits = hits.size();
G4cout << "Event: " << i << " Number of Hits: " << NbHits << G4endl;
np->Fill(NbHits);
for(G4int ii = 0; ii < NbHits; ii++)
{
PhotonHit* photonHit = dynamic_cast<PhotonHit*>(hits.at(ii));
time->Fill(photonHit->GetTime());
wl->Fill(photonHit->GetWavelength());
std::cout << "sensor ID: " << photonHit->GetId() << std::endl;
switch(photonHit->GetId())
{
case 0:
time0->Fill(photonHit->GetTime());
break;
case 1:
time1->Fill(photonHit->GetTime());
break;
case 2:
time2->Fill(photonHit->GetTime());
break;
case 3:
time3->Fill(photonHit->GetTime());
break;
case 4:
time4->Fill(photonHit->GetTime());
break;
}
switch(photonHit->GetPid())
{
case 0:
wlce->Fill(photonHit->GetWavelength());
break;
case 1:
wlsc->Fill(photonHit->GetWavelength());
break;
default:
wl->Fill(photonHit->GetWavelength());
}
pos2->Fill(photonHit->GetPosition().getZ(), photonHit->GetPosition().getY());
pos0->Fill(photonHit->GetPosition().getZ());
}
}
}
}
outfile->cd();
outfile->Write();
}