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