-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpstat_module.cc
176 lines (148 loc) · 5.47 KB
/
pstat_module.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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
////////////////////////////////////////////////////////////////////////
// Class: NeutronCaptureMC
// Plugin Type: analyzer (art v3_01_02)
// File: NeutronCaptureMC_module.cc
//
// Generated at Tue Mar 26 11:44:15 2019 by Jingbo Wang using cetskelgen
// from cetlib version v3_05_01.
////////////////////////////////////////////////////////////////////////
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "TTree.h"
#include "art_root_io/TFileService.h"
#include "TH1F.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "larcore/CoreUtils/ServiceUtil.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardataobj/Simulation/SimEnergyDeposit.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Slice.h"
#include "lardataobj/Simulation/SimPhotons.h"
#include "larsim/IonizationScintillation/ISCalcSeparate.h"
#include "larsim/PhotonPropagation/PhotonVisibilityService.h"
#include "larsim/Simulation/LArG4Parameters.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "larcore/Geometry/Geometry.h"
#include "TGeoMaterial.h"
#include "TGeoElement.h"
namespace test {
class NeutronCaptureMC;
}
class test::NeutronCaptureMC : public art::EDAnalyzer {
public:
explicit NeutronCaptureMC(fhicl::ParameterSet const& p);
// The compiler-generated destructor is fine for non-base
// classes without bare pointers or other resource use.
// Plugins should not be copied or assigned.
NeutronCaptureMC(NeutronCaptureMC const&) = delete;
NeutronCaptureMC(NeutronCaptureMC&&) = delete;
NeutronCaptureMC& operator=(NeutronCaptureMC const&) = delete;
NeutronCaptureMC& operator=(NeutronCaptureMC&&) = delete;
// Required functions.
void analyze(art::Event const& e) override;
void beginJob();
void endJob();
void reconfigure(fhicl::ParameterSet const &p);
private:
//analyzer;
// Declare member data here.
TTree *fTree; ///< My TTree
int fEvent;
std::string fTruthLabel;
// Reco variables
std::vector<float> en_n;
std::vector<float> en_e;
std::vector<float> en_mu;
std::vector<int> pdg;
std::vector<float> Emuv;
std::vector<float> nnv;
double fEkGen;
};
test::NeutronCaptureMC::NeutronCaptureMC(fhicl::ParameterSet const& p)
: EDAnalyzer(p)
// More initializers here.
{
this->reconfigure(p);
// Call appropriate consumes<>() for any products to be retrieved by this module.
}
void test::NeutronCaptureMC::beginJob() {
//create output tree
//art::ServiceHandle<geo::Geometry> geo_serv;
art::ServiceHandle<art::TFileService> tfs;
fTree = tfs->make<TTree>("mytree", "My Tree");
fTree->Branch("event", &fEvent, "event/I");
fTree->Branch("en_e",&en_e);
fTree->Branch("en_mu",&en_mu);
fTree->Branch("en_n",&en_n);
fTree->Branch("pdg",&pdg);
fTree->Branch("Emuv",&Emuv);
fTree->Branch("nnv",&nnv);
}
void test::NeutronCaptureMC::endJob() {
}
void test::NeutronCaptureMC::reconfigure(fhicl::ParameterSet const & p)
{
fTruthLabel = p.get<std::string>("TruthLabel");
}
void test::NeutronCaptureMC::analyze(art::Event const& e)
{
//art::ServiceHandle<geo::Geometry> geo_serv;
// Implementation of required member function here
en_n.clear();
en_e.clear();
en_mu.clear();
pdg.clear();
Emuv.clear();
nnv.clear();
fEvent = e.id().event();
double Emu=0;
double nn=0;
if(!e.isRealData()) {
// Get the list of MC particles from GEANT
auto mcParticles = e.getValidHandle<std::vector<simb::MCParticle>>(fTruthLabel);
if(mcParticles.isValid()) {
for(auto &trueParticle : *mcParticles) {
fEkGen = (std::sqrt(trueParticle.P()*trueParticle.P() + trueParticle.Mass()*trueParticle.Mass()) - trueParticle.Mass()) * 1000; // MeVs
pdg.push_back(trueParticle.PdgCode());
if(trueParticle.Process() == "primary") {
std::cout<<"PdgCode, Process, Total E(GeV), KineticE(MeV) = "<<trueParticle.PdgCode()<<", "<<trueParticle.Process()<<", "<<trueParticle.E()<<", "<<fEkGen<<std::endl;
}
if(trueParticle.Process() == "primary"&&trueParticle.PdgCode() ==2112) {
en_n.push_back(fEkGen);
}
if(trueParticle.Process() == "primary"&&trueParticle.PdgCode() ==11) {
en_e.push_back(fEkGen);
}
if(trueParticle.Process() == "primary"&&trueParticle.PdgCode() ==13) {
en_mu.push_back(fEkGen);
Emu=fEkGen;
}
if(trueParticle.Process() == "primary"&&trueParticle.PdgCode() ==-13) {
en_mu.push_back(fEkGen);
}
if(trueParticle.PdgCode() ==2112) {
en_n.push_back(fEkGen);
nn++;
}
}
Emuv.push_back(Emu);
nnv.push_back(nn);
std::cout<<"muon energy: "<<Emu<<" number of neutron: "<<nn<<std::endl;
}
}
// Analysis goes here...
fTree->Fill();
}
DEFINE_ART_MODULE(test::NeutronCaptureMC)