Belle II Software  release-05-02-19
ClawStudyModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - 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/claw/modules/ClawStudyModule.h>
12 #include <beast/claw/dataobjects/ClawSimHit.h>
13 #include <beast/claw/dataobjects/ClawHit.h>
14 #include <generators/SAD/dataobjects/SADMetaHit.h>
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/gearbox/GearDir.h>
17 #include <framework/logging/Logger.h>
18 #include <cmath>
19 
20 #include <fstream>
21 #include <string>
22 
23 // ROOT
24 #include <TH1.h>
25 #include <TH2.h>
26 
27 int eventNum = 0;
28 
29 using namespace std;
30 
31 using namespace Belle2;
32 using namespace claw;
33 
34 //-----------------------------------------------------------------
35 // Register the Module
36 //-----------------------------------------------------------------
37 REG_MODULE(ClawStudy)
38 
39 //-----------------------------------------------------------------
40 // Implementation
41 //-----------------------------------------------------------------
42 
44 {
45  // Set module properties
46  setDescription("Study module for Claws (BEAST)");
47 
48  addParam("Ethres", m_Ethres, "Energy threshold in MeV", 0.0);
49 }
50 
51 ClawStudyModule::~ClawStudyModule()
52 {
53 }
54 
55 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
56 void ClawStudyModule::defineHisto()
57 {
58  for (int i = 0; i < 8; i++) {
59  h_claws_Evtof1[i] = new TH2F(TString::Format("claws_Evtof1_%d", i), "Energy deposited [MeV] vs TOF [ns] - all", 5000, 0., 1000.,
60  1000, 0., 10.);
61  h_claws_Evtof2[i] = new TH2F(TString::Format("claws_Evtof2_%d", i), "Energy deposited [MeV] vs TOF [ns] - only photons", 5000, 0.,
62  1000., 1000, 0., 10.);
63  h_claws_Evtof3[i] = new TH2F(TString::Format("claws_Evtof3_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 5000, 0.,
64  1000., 1000, 0., 10.);
65  h_claws_Evtof4[i] = new TH2F(TString::Format("claws_Evtof4_%d", i), "Energy deposited [MeV] vs TOF [ns] - only e+/e-", 5000, 0.,
66  1000., 1000, 0., 10.);
67  h_claws_edep[i] = new TH1F(TString::Format("claws_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 10.);
68  h_Wclaws_edep[i] = new TH1F(TString::Format("Wclaws_edep_%d", i), "Energy deposited [MeV]", 5000, 0., 10.);
69  }
70 
71  h_claws_hitrate1 = new TH1F("claws_hitrate1", "Hit distributions", 8, 0., 8.);
72  h_claws_hitrate2 = new TH1F("claws_hitrate2", "Hit distributions", 8, 0., 8.);
73  h_claws_hitrate1W = new TH1F("claws_hitrate1W", "Hit distributions", 8, 0., 8.);
74  h_claws_hitrate2W = new TH1F("claws_hitrate2W", "Hit distributions", 8, 0., 8.);
75 
76  h_claws_hitrate1->Sumw2();
77  h_claws_hitrate1W->Sumw2();
78  h_claws_hitrate2->Sumw2();
79  h_claws_hitrate2W->Sumw2();
80 
81  h_claws_rs_hitrate1 = new TH2F("claws_rs_hitrate1", "Hit distributions vs rs", 8, 0., 8., 12, 0., 12.);
82  h_claws_rs_hitrate2 = new TH2F("claws_rs_hitrate2", "Hit distributions vs rs", 8, 0., 8., 12, 0., 12.);
83  h_claws_rs_hitrate1W = new TH2F("claws_rs_hitrate1W", "Hit distributions vs rs", 8, 0., 8., 12, 0., 12.);
84  h_claws_rs_hitrate2W = new TH2F("claws_rs_hitrate2W", "Hit distributions vs rs", 8, 0., 8., 12, 0., 12.);
85 
86  h_claws_rs_hitrate1->Sumw2();
87  h_claws_rs_hitrate1W->Sumw2();
88  h_claws_rs_hitrate2->Sumw2();
89  h_claws_rs_hitrate2W->Sumw2();
90 
91  for (int i = 0; i < 8; i++) {
92  h_claws_rate1[i] = new TH1F(TString::Format("claws_rate1_%d", i), "PE distributions", 5000, 0., 5000.);
93  h_claws_rate2[i] = new TH1F(TString::Format("claws_rate2_%d", i), "PE distributions", 5000, 0., 5000.);
94  h_claws_rate1W[i] = new TH1F(TString::Format("claws_rate1W_%d", i), "PE distributions", 5000, 0., 5000.);
95  h_claws_rate2W[i] = new TH1F(TString::Format("claws_rate2W_%d", i), "PE distributions", 5000, 0., 5000.);
96  h_claws_pe1[i] = new TH2F(TString::Format("claws_pe1_%d", i), "PE distributions", 5000, 0., 5000., 1000, 0., 1000.);
97  h_claws_pe2[i] = new TH2F(TString::Format("claws_pe2_%d", i), "PE distributions", 5000, 0., 5000., 1000, 0., 1000.);
98  h_claws_pe1W[i] = new TH2F(TString::Format("claws_pe1W_%d", i), "PE distributions", 5000, 0., 5000., 1000, 0., 1000.);
99  h_claws_pe2W[i] = new TH2F(TString::Format("claws_pe2W_%d", i), "PE distributions", 5000, 0., 5000., 1000, 0., 1000.);
100 
101  h_claws_rs_rate1[i] = new TH2F(TString::Format("claws_rs_rate1_%d", i), "PE distributions", 5000, 0., 5000., 12, 0., 12.);
102  h_claws_rs_rate2[i] = new TH2F(TString::Format("claws_rs_rate2_%d", i), "PE distributions", 5000, 0., 5000., 12, 0., 12.);
103  h_claws_rs_rate1W[i] = new TH2F(TString::Format("claws_rs_rate1W_%d", i), "PE distributions", 5000, 0., 5000., 12, 0., 12.);
104  h_claws_rs_rate2W[i] = new TH2F(TString::Format("claws_rs_rate2W_%d", i), "PE distributions", 5000, 0., 5000., 12, 0., 12.);
105 
106  h_claws_rate1[i]->Sumw2();
107  h_claws_rate2[i]->Sumw2();
108  h_claws_rate1W[i]->Sumw2();
109  h_claws_rate2W[i]->Sumw2();
110  h_claws_rs_rate1[i]->Sumw2();
111  h_claws_rs_rate2[i]->Sumw2();
112  h_claws_rs_rate1W[i]->Sumw2();
113  h_claws_rs_rate2W[i]->Sumw2();
114  h_claws_pe1[i]->Sumw2();
115  h_claws_pe2[i]->Sumw2();
116  h_claws_pe1W[i]->Sumw2();
117  h_claws_pe2W[i]->Sumw2();
118  }
119 }
120 
121 
122 void ClawStudyModule::initialize()
123 {
124  B2INFO("ClawStudyModule: Initialize");
125 
126  REG_HISTOGRAM
127 
128  //get CLAWS paramters ie energy threshold
129  getXMLData();
130 }
131 
132 void ClawStudyModule::beginRun()
133 {
134 }
135 
136 void ClawStudyModule::event()
137 {
138  //Here comes the actual event processing
139  StoreArray<ClawSimHit> SimHits;
140  StoreArray<ClawHit> Hits;
141  StoreArray<SADMetaHit> MetaHits;
142 
143  double rate = 0;
144  int ring_section = -1;
145  for (const auto& MetaHit : MetaHits) {
146  rate = MetaHit.getrate();
147  ring_section = MetaHit.getring_section() - 1;
148  }
149 
150  //number of entries in SimHits
151  //int nSimHits = SimHits.getEntries();
152 
153  //loop over all SimHit entries
154  for (int i = 0; i < (int) SimHits.getEntries(); i++) {
155  ClawSimHit* aHit = SimHits[i];
156  int detNB = aHit->getCellId();
157  if (detNB < 8) {
158  //int trkID = aHit->getTrackId();
159  int pdg = aHit->getPDGCode();
160  double Edep = aHit->getEnergyDep() * 1e3; //GeV -> MeV
161  double tof = aHit->getFlightTime(); //ns
162 
163  h_claws_Evtof1[detNB]->Fill(tof, Edep);
164  if (pdg == 22) h_claws_Evtof2[detNB]->Fill(tof, Edep);
165  else if (fabs(pdg) == 11) h_claws_Evtof3[detNB]->Fill(tof, Edep);
166  else h_claws_Evtof4[detNB]->Fill(tof, Edep);
167  if (Edep > m_Ethres) {
168  h_claws_edep[detNB]->Fill(Edep);
169  h_Wclaws_edep[detNB]->Fill(Edep, rate);
170  }
171  }
172  }
173 
174  for (const auto& Hit : Hits) {
175  const int detNb = Hit.getdetNb();
176  if (detNb < 8) {
177  const int timebin = Hit.gettime();
178  const float edep = Hit.getedep();
179  const float pe = Hit.getPE();
180  h_claws_hitrate1->Fill(detNb);
181  h_claws_hitrate1W->Fill(detNb, rate);
182  h_claws_rate1[detNb]->Fill(pe);
183  h_claws_rate1W[detNb]->Fill(pe, rate);
184  h_claws_rs_rate1[detNb]->Fill(pe, ring_section);
185  h_claws_rs_rate1W[detNb]->Fill(pe, ring_section, rate);
186  h_claws_rs_hitrate1->Fill(detNb, ring_section);
187  h_claws_rs_hitrate1W->Fill(detNb, ring_section, rate);
188  h_claws_pe1[detNb]->Fill(timebin, pe);
189  h_claws_pe1W[detNb]->Fill(timebin, pe, rate);
190  if (edep > m_Ethres) {
191  h_claws_hitrate2->Fill(detNb);
192  h_claws_hitrate2W->Fill(detNb, rate);
193  h_claws_rate2[detNb]->Fill(pe);
194  h_claws_rate2W[detNb]->Fill(pe, rate);
195  h_claws_rs_rate2[detNb]->Fill(pe, ring_section);
196  h_claws_rs_rate2W[detNb]->Fill(pe, ring_section, rate);
197  h_claws_rs_hitrate2->Fill(detNb, ring_section);
198  h_claws_rs_hitrate2W->Fill(detNb, ring_section, rate);
199  h_claws_pe2[detNb]->Fill(timebin, pe);
200  h_claws_pe2W[detNb]->Fill(timebin, pe, rate);
201  }
202  }
203  }
204 
205 }
206 
207 //read energy threshold from CLAW.xml
208 void ClawStudyModule::getXMLData()
209 {
210  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"CLAW\"]/Content/");
211  m_Ethres = content.getDouble("Ethres");
212 
213  B2INFO("ClawStudy");
214 }
215 
216 
217 void ClawStudyModule::endRun()
218 {
219 
220 
221 
222 }
223 
224 void ClawStudyModule::terminate()
225 {
226 }
227 
228 
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ClawSimHit::getCellId
int getCellId() const
Get Cell ID.
Definition: ClawSimHit.h:106
Belle2::ClawSimHit::getFlightTime
double getFlightTime() const
Get Flight time from IP.
Definition: ClawSimHit.h:121
Belle2::ClawSimHit::getEnergyDep
double getEnergyDep() const
Get Deposit energy.
Definition: ClawSimHit.h:126
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::ClawSimHit::getPDGCode
int getPDGCode() const
Get Particle PDG (can be one of secondaries)
Definition: ClawSimHit.h:116
Belle2::claw::ClawStudyModule
Study module for Claws (BEAST)
Definition: ClawStudyModule.h:39
Belle2::ClawSimHit
ClassClawSimHit - Geant4 simulated hit for the Claw crystal in beast.
Definition: ClawSimHit.h:41
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
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::Hit
Structure to hold some of the calpulse data.
Definition: TOPTimeBaseCalibratorModule.h:40