Belle II Software  release-06-00-14
AnalysisPhase1StudyModule.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <beast/analysis/modules/AnalysisPhase1StudyModule.h>
10 #include <mdst/dataobjects/MCParticle.h>
11 #include <generators/SAD/dataobjects/SADMetaHit.h>
12 #include <framework/datastore/StoreArray.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
15 #include <cmath>
16 
17 #include <fstream>
18 #include <string>
19 
20 // ROOT
21 #include <TVector3.h>
22 #include <TH1.h>
23 #include <TH2.h>
24 
25 using namespace std;
26 
27 using namespace Belle2;
28 
29 //using namespace analysis;
30 
31 //-----------------------------------------------------------------
32 // Register the Module
33 //-----------------------------------------------------------------
34 REG_MODULE(AnalysisPhase1Study)
35 
36 //-----------------------------------------------------------------
37 // Implementation
38 //-----------------------------------------------------------------
39 
40 AnalysisPhase1StudyModule::AnalysisPhase1StudyModule() : HistoModule()
41 {
42  // Set module properties
43  setDescription("Study module for BEAST");
44 
45 }
46 
47 AnalysisPhase1StudyModule::~AnalysisPhase1StudyModule()
48 {
49 }
50 
51 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
52 void AnalysisPhase1StudyModule::defineHisto()
53 {
54  h_count = new TH1F("h_count", "", 10, 0., 10.);
55  for (int i = 0; i < 10; i++) {
56  h_prodvtx[i] = new TH2F(TString::Format("h_prodvtx_%d", i), "", 200, -400., 400., 200, 0., 400.);
57  h_decavtx[i] = new TH2F(TString::Format("h_decavtx_%d", i), "", 200, -400., 400., 200, 0., 400.);
58  h_kineticvz[i] = new TH2F(TString::Format("h_kineticvz_%d", i), "", 200, -400., 400., 1000, 0., 10.);
59  h_kineticvz1[i] = new TH2F(TString::Format("h_kineticvz1_%d", i), "", 200, -400., 400., 1000, 0., 10.);
60  h_kineticvz2[i] = new TH2F(TString::Format("h_kineticvz2_%d", i), "", 200, -400., 400., 1000, 0., 10.);
61  h_kineticvz_zoom[i] = new TH2F(TString::Format("h_kineticvz_zoom_%d", i), "", 200, -400., 400., 1000, 0., 0.01);
62  h_Wkineticvz[i] = new TH2F(TString::Format("h_Wkineticvz_%d", i), "", 200, -400., 400., 1000, 0., 10.);
63  h_Wkineticvz_zoom[i] = new TH2F(TString::Format("h_Wkineticvz_zoom_%d", i), "", 200, -400., 400., 1000, 0., 0.01);
64  h_phive[i] = new TH2F(TString::Format("h_phive_%d", i), "", 360, -180., 180., 1000, 0., 10.);
65  h_phivz[i] = new TH2F(TString::Format("h_phivz_%d", i), "", 200, -400., 400., 360, -180., 180.);
66  h_rve[i] = new TH2F(TString::Format("h_rve_%d", i), "", 200, 0., 200., 1000, 0., 10.);
67  h_thetavz[i] = new TH2F(TString::Format("h_thetavz_%d", i), "", 200, -400., 400., 180, 0., 180.);
68  h_kineticvz[i]->Sumw2();
69  h_kineticvz1[i]->Sumw2();
70  h_kineticvz2[i]->Sumw2();
71  h_kineticvz_zoom[i]->Sumw2();
72  h_Wkineticvz[i]->Sumw2();
73  h_Wkineticvz_zoom[i]->Sumw2();
74  h_phive[i]->Sumw2();
75  h_phivz[i]->Sumw2();
76  h_rve[i]->Sumw2();
77  }
78  h_g4_xy = new TH2F("h_g4_xy", "", 100, -5.99, 5.99, 100, -5.99, 5.99);
79  for (int i = 0; i < 2; i++) {
80  h_sad_xy[i] = new TH2F(TString::Format("h_sad_xy_%d", i), "", 100, -5.99, 5.99, 100, -5.99, 5.99);
81  h_sad_sir[i] = new TH1F(TString::Format("h_sad_sir_%d", i), "", 100, -3.99, 3.99);
82  h_z[i] = new TH1F(TString::Format("h_z_%d", i), "", 100, -3.99, 3.99);
83  h_sad_sall[i] = new TH1F(TString::Format("h_sad_sall_%d", i), "", 100, -1499.99, 1499.99);
84  h_sad_sE[i] = new TH2F(TString::Format("h_sad_sE_%d", i), "", 100, -3.99, 3.99, 1000, 0., 10.);
85  h_sad_sraw[i] = new TH1F(TString::Format("h_sad_sraw_%d", i), "", 100, -1499.99, 1499.99);
86  h_sad_E[i] = new TH1F(TString::Format("h_sad_E_%d", i), "", 8000, 0, 7.99);
87  h_sad_s[i] = new TH1F(TString::Format("h_sad_s_%d", i), "", 400, -3.99, 3.99);
88  }
89  h_dpx = new TH1F("h_dpx", "", 1000, -1., 1.);
90  h_dpy = new TH1F("h_dpy", "", 1000, -1., 1.);
91  h_dE = new TH1F("h_dE", "", 1000, -1., 1.);
92  h_px = new TH1F("h_px", "", 10000, -10., 10.);
93  h_py = new TH1F("h_py", "", 10000, -10., 10.);
94  h_pz = new TH1F("h_pz", "", 10000, -10., 10.);
95  h_dx = new TH1F("h_dx", "", 10000, -400., 400.);
96  h_dy = new TH1F("h_dy", "", 10000, -400., 400.);
97  h_dz = new TH1F("h_dz", "", 10000, -400., 400.);
98  h_E = new TH1F("h_E", "", 8000, 0., 7.99);
99  h_P = new TH1F("h_P", "", 8000, 0., 7.99);
100  h_E->Sumw2();
101  h_P->Sumw2();
102 }
103 
104 void AnalysisPhase1StudyModule::initialize()
105 {
106  B2INFO("AnalysisPhase1StudyModule: Initialize");
107 
108  REG_HISTOGRAM
109 
110 }
111 
112 void AnalysisPhase1StudyModule::beginRun()
113 {
114 }
115 
116 void AnalysisPhase1StudyModule::event()
117 {
118  double px = 0;
119  double py = 0;
120  double x = 0;
121  double y = 0;
122  double s = 0;
123  //double xraw = 0;
124  //double yraw = 0;
125  //double sraw = 0;
126  double E = 0;
127  double rate = 0;
128  StoreArray<SADMetaHit> sadMetaHits;
129  for (const auto& sadMetaHit : sadMetaHits) {
130  px = sadMetaHit.getpx();
131  py = sadMetaHit.getpy();
132  x = sadMetaHit.getx();
133  y = sadMetaHit.gety();
134  s = sadMetaHit.gets();
135  //xraw = sadMetaHit.getxraw();
136  //yraw = sadMetaHit.getyraw();
137  double sraw = sadMetaHit.getsraw();
138  E = sadMetaHit.getE();
139  rate = sadMetaHit.getrate();
140  h_sad_sir[0]->Fill(s / 100., rate);
141  h_sad_sall[0]->Fill(s / 100., rate);
142  h_sad_sE[0]->Fill(s / 100., E);
143  h_sad_s[0]->Fill(s / 100.);
144  h_sad_E[0]->Fill(E);
145  h_sad_sraw[0]->Fill(sraw / 100., rate);
146  if (-400. < s && s < 400.) {
147  h_sad_xy[1]->Fill(x, y);
148  h_sad_sir[1]->Fill(s / 100., rate);
149  h_sad_sall[1]->Fill(s / 100., rate);
150  h_sad_sraw[1]->Fill(sraw / 100., rate);
151  h_sad_sE[1]->Fill(s / 100., E);
152  h_sad_s[1]->Fill(s / 100.);
153  h_sad_E[1]->Fill(E);
154  }
155  }
156  int counter = 0;
157  StoreArray<MCParticle> mcParticles;
158  for (const auto& mcParticle : mcParticles) { // start loop over all MC particles
159  int PDG = mcParticle.getPDG();
160  //int* Mother = mcParticle.getMother();
161  float Mass = mcParticle.getMass();
162  float Energy = mcParticle.getEnergy();
163  float Kinetic = Energy - Mass;
164  float prodvtx[3];
165  prodvtx[0] = mcParticle.getProductionVertex().X();
166  prodvtx[1] = mcParticle.getProductionVertex().Y();
167  prodvtx[2] = mcParticle.getProductionVertex().Z();
168  float prodr = sqrt(prodvtx[0] * prodvtx[0] + prodvtx[1] * prodvtx[1]);
169  float decavtx[3];
170  decavtx[0] = mcParticle.getDecayVertex().X();
171  decavtx[1] = mcParticle.getDecayVertex().Y();
172  decavtx[2] = mcParticle.getDecayVertex().Z();
173  float decar = sqrt(decavtx[0] * decavtx[0] + decavtx[1] * decavtx[1]);
174  float mom[3];
175  mom[0] = mcParticle.getMomentum().X();
176  mom[1] = mcParticle.getMomentum().Y();
177  mom[2] = mcParticle.getMomentum().Z();
178  float theta = mcParticle.getMomentum().Theta() * TMath::RadToDeg();
179  float phi = mcParticle.getMomentum().Phi() * TMath::RadToDeg();
180  int partID[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
181 
182  if (PDG == Const::electron.getPDGCode()) partID[0] = 1; //positron
183  else if (PDG == -Const::electron.getPDGCode()) partID[1] = 1; //electron
184  else if (PDG == Const::photon.getPDGCode()) partID[2] = 1; //photon
185  else if (PDG == Const::neutron.getPDGCode()) partID[3] = 1; //neutron
186  else if (PDG == Const::proton.getPDGCode()) partID[4] = 1; //proton
187  else if (PDG == 1000080160) partID[5] = 1; // O
188  else if (PDG == 1000060120) partID[6] = 1; // C
189  else if (PDG == 1000020040) partID[7] = 1; // He
190  else partID[8] = 1;
191  for (int i = 0; i < 8; i++) {
192  if (partID[i] == 1) {
193  h_count->Fill(i);
194  h_prodvtx[i]->Fill(prodvtx[2], prodr);
195  h_decavtx[i]->Fill(decavtx[2], decar);
196  h_kineticvz[i]->Fill(prodvtx[2], Kinetic);
197  h_Wkineticvz[i]->Fill(prodvtx[2], Kinetic, rate);
198  h_kineticvz_zoom[i]->Fill(prodvtx[2], Kinetic);
199  h_Wkineticvz_zoom[i]->Fill(prodvtx[2], Kinetic, rate);
200  h_thetavz[i]->Fill(prodvtx[2], theta);
201  h_phivz[i]->Fill(prodvtx[2], phi);
202  h_phive[i]->Fill(phi, Kinetic);
203  h_rve[i]->Fill(prodr, Kinetic);
204  if (theta == 0) {
205  h_kineticvz1[i]->Fill(prodvtx[2], Kinetic);
206  }
207  if (phi == 0) {
208  h_kineticvz2[i]->Fill(prodvtx[2], Kinetic);
209  }
210  }
211  }
212 
213  //h_count->Fill(9);
214  h_prodvtx[8]->Fill(prodvtx[2], prodr);
215  h_decavtx[8]->Fill(decavtx[2], decar);
216  h_kineticvz[8]->Fill(prodvtx[2], Kinetic);
217  h_kineticvz_zoom[8]->Fill(prodvtx[2], Kinetic);
218  h_Wkineticvz[8]->Fill(prodvtx[2], Kinetic, rate);
219  h_Wkineticvz_zoom[8]->Fill(prodvtx[2], Kinetic, rate);
220  h_thetavz[8]->Fill(prodvtx[2], theta);
221  h_phivz[8]->Fill(prodvtx[2], phi);
222  h_phive[8]->Fill(phi, Kinetic);
223  h_rve[8]->Fill(prodr, Kinetic);
224  if (theta == 0) {
225  h_kineticvz1[8]->Fill(prodvtx[2], Kinetic);
226  }
227  if (phi == 0) {
228  h_kineticvz2[8]->Fill(prodvtx[2], Kinetic);
229  }
230  if (counter == 0) {
231  h_count->Fill(9);
232  h_prodvtx[9]->Fill(prodvtx[2], prodr);
233  h_decavtx[9]->Fill(decavtx[2], decar);
234  h_kineticvz[9]->Fill(prodvtx[2], Kinetic);
235  h_kineticvz_zoom[9]->Fill(prodvtx[2], Kinetic);
236  h_Wkineticvz[9]->Fill(prodvtx[2], Kinetic, rate);
237  h_Wkineticvz_zoom[9]->Fill(prodvtx[2], Kinetic, rate);
238  h_thetavz[9]->Fill(prodvtx[2], theta);
239  h_phivz[9]->Fill(prodvtx[2], phi);
240 
241  h_dx->Fill(x - prodvtx[0]);
242  h_dy->Fill(-y - prodvtx[1]);
243  h_dz->Fill(-s - prodvtx[2]);
244 
245  h_z[0]->Fill(prodvtx[2] / 100.);
246  h_z[1]->Fill(prodvtx[2] / 100., rate);
247 
248  h_dpx->Fill(px - mom[0]);
249  h_dpy->Fill(-py - mom[1]);
250  h_dE->Fill(E - Energy);
251 
252  h_g4_xy->Fill(prodvtx[0], prodvtx[1]);
253 
254  h_px->Fill(mom[0]);
255  h_py->Fill(mom[1]);
256  h_pz->Fill(mom[2]);
257  h_E->Fill(Kinetic);
258  h_P->Fill(sqrt(mom[0] * mom[0] + mom[1] * mom[1] + mom[2] * mom[2]));
259  }
260  counter++;
261  }
262 
263 }
264 
265 void AnalysisPhase1StudyModule::endRun()
266 {
267 
268 
269 }
270 
271 void AnalysisPhase1StudyModule::terminate()
272 {
273 }
274 
275 
Belle2::HistoModule
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17