Belle II Software  release-06-02-00
FANGSStudyModule.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/fangs/modules/FANGSStudyModule.h>
10 #include <beast/fangs/dataobjects/FANGSSimHit.h>
11 #include <beast/fangs/dataobjects/FANGSHit.h>
12 #include <framework/datastore/StoreArray.h>
13 #include <framework/gearbox/GearDir.h>
14 #include <framework/logging/Logger.h>
15 #include <framework/gearbox/Const.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 int eventNum = 0;
27 
28 using namespace std;
29 
30 using namespace Belle2;
31 using namespace fangs;
32 
33 //-----------------------------------------------------------------
34 // Register the Module
35 //-----------------------------------------------------------------
36 REG_MODULE(FANGSStudy)
37 
38 //-----------------------------------------------------------------
39 // Implementation
40 //-----------------------------------------------------------------
41 
42 FANGSStudyModule::FANGSStudyModule() : HistoModule()
43 {
44  // Set module properties
45  setDescription("Study module for Fangs (BEAST)");
46 
47 }
48 
49 FANGSStudyModule::~FANGSStudyModule()
50 {
51 }
52 
53 //This module is a histomodule. Any histogram created here will be saved by the HistoManager module
54 void FANGSStudyModule::defineHisto()
55 {
56  h_time = new TH2F("h_time", "Detector # vs. time", 20, 0., 20., 1000, 0., 750.);
57  h_time->Sumw2();
58  h_timeWeighted = new TH2F("h_timeWeigthed", "Detector # vs. time weighted by the energy deposited", 20, 0., 20., 1000, 0., 750.);
59  h_timeWeighted->Sumw2();
60  h_timeThres = new TH2F("h_timeThres", "Detector # vs. time", 20, 0., 20., 750, 0., 750.);
61  h_timeThres->Sumw2();
62  h_timeWeightedThres = new TH2F("h_timeWeigthedThres", "Detector # vs. time weighted by the energy deposited", 20, 0., 20., 750, 0.,
63  750.);
64  h_timeWeightedThres->Sumw2();
65  h_edep = new TH2F("h_edep", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
66  h_edep->Sumw2();
67  h_edep1 = new TH2F("h_edep1", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
68  h_edep1->Sumw2();
69  h_edep2 = new TH2F("h_edep2", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
70  h_edep2->Sumw2();
71  h_edep3 = new TH2F("h_edep3", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
72  h_edep3->Sumw2();
73 
74  h_edepThres = new TH2F("h_edepThres", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
75  h_edepThres->Sumw2();
76  h_edepThres1 = new TH2F("h_edepThres1", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
77  h_edepThres1->Sumw2();
78  h_edepThres2 = new TH2F("h_edepThres2", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
79  h_edepThres2->Sumw2();
80  h_edepThres3 = new TH2F("h_edepThres3", "Time bin # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
81  h_edepThres3->Sumw2();
82  for (int i = 0; i < 3; i++) {
83  h_zvedep[i] = new TH1F(TString::Format("h_zvedep_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25.);
84  h_zvedep[i]->Sumw2();
85 
86  h_xvzvedep[i] = new TH2F(TString::Format("h_xvzvedep_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
87  h_xvzvedep[i]->Sumw2();
88 
89  h_yvzvedep[i] = new TH2F(TString::Format("h_yvzvedep_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
90  h_yvzvedep[i]->Sumw2();
91 
92  h_rvzvedep[i] = new TH2F(TString::Format("h_rvzvedep_%d", i) , "edep [MeV] vs. z [cm]", 2000, 0., 25., 2000, -25., 25.);
93  h_rvzvedep[i]->Sumw2();
94 
95  h_xvyvedep[i] = new TH2F(TString::Format("h_xvyvedep_%d", i) , "edep [MeV] vs. y [cm]", 2000, -25., 25., 2000, -25., 25.);
96  h_xvyvedep[i]->Sumw2();
97 
98  h_zvedepW[i] = new TH1F(TString::Format("h_zvedepW_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25.);
99  h_zvedepW[i]->Sumw2();
100 
101  h_xvzvedepW[i] = new TH2F(TString::Format("h_xvzvedepW_%d", i) , "edep [MeV] vs. x vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
102  h_xvzvedepW[i]->Sumw2();
103 
104  h_yvzvedepW[i] = new TH2F(TString::Format("h_yvzvedepW_%d", i) , "edep [MeV] vs. y vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
105  h_yvzvedepW[i]->Sumw2();
106 
107  h_xvyvedepW[i] = new TH2F(TString::Format("h_xvyvedepW_%d", i) , "edep [MeV] vs. x vs. y [cm]", 2000, -25., 25., 2000, -25., 25.);
108  h_xvyvedepW[i]->Sumw2();
109 
110  h_rvzvedepW[i] = new TH2F(TString::Format("h_rvzvedepW_%d", i) , "edep [MeV] vs. z [cm]", 2000, 0., 25., 2000, -25., 25.);
111  h_rvzvedepW[i]->Sumw2();
112 
113 
114  h_zvedepT[i] = new TH1F(TString::Format("h_zvedepT_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25.);
115  h_zvedepT[i]->Sumw2();
116 
117  h_xvzvedepT[i] = new TH2F(TString::Format("h_xvzvedepT_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
118  h_xvzvedepT[i]->Sumw2();
119 
120  h_yvzvedepT[i] = new TH2F(TString::Format("h_yvzvedepT_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
121  h_yvzvedepT[i]->Sumw2();
122 
123  h_rvzvedepT[i] = new TH2F(TString::Format("h_rvzvedepT_%d", i) , "edep [MeV] vs. z [cm]", 2000, 0., 25., 2000, -25., 25.);
124  h_rvzvedepT[i]->Sumw2();
125 
126  h_xvyvedepT[i] = new TH2F(TString::Format("h_xvyvedepT_%d", i) , "edep [MeV] vs. y [cm]", 2000, -25., 25., 2000, -25., 25.);
127  h_xvyvedepT[i]->Sumw2();
128 
129  h_zvedepWT[i] = new TH1F(TString::Format("h_zvedepWT_%d", i) , "edep [MeV] vs. z [cm]", 2000, -25., 25.);
130  h_zvedepWT[i]->Sumw2();
131 
132  h_xvzvedepWT[i] = new TH2F(TString::Format("h_xvzvedepWT_%d", i) , "edep [MeV] vs. x vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
133  h_xvzvedepWT[i]->Sumw2();
134 
135  h_yvzvedepWT[i] = new TH2F(TString::Format("h_yvzvedepWT_%d", i) , "edep [MeV] vs. y vs. z [cm]", 2000, -25., 25., 2000, -25., 25.);
136  h_yvzvedepWT[i]->Sumw2();
137 
138  h_xvyvedepWT[i] = new TH2F(TString::Format("h_xvyvedepWT_%d", i) , "edep [MeV] vs. x vs. y [cm]", 2000, -25., 25., 2000, -25., 25.);
139  h_xvyvedepWT[i]->Sumw2();
140 
141  h_rvzvedepWT[i] = new TH2F(TString::Format("h_rvzvedepWT_%d", i) , "edep [MeV] vs. z [cm]", 2000, 0., 25., 2000, -25., 25.);
142  h_rvzvedepWT[i]->Sumw2();
143  }
144  h_Edep = new TH2F("h_Edep", "det # # vs. energy deposited", 20, 0., 20., 1000, 0., 10.);
145  h_pxNb = new TH2F("h_pxNb", "det # # vs. nb pixel", 20, 0., 20., 1000, 0., 1000.);
146  for (int i = 0; i < 15; i++) {
147  h_cvr[i] = new TH2F(TString::Format("cvr_%d", i) , " col v. row", 80, 0., 80., 336, 0., 336.);
148  }
149 }
150 
151 
152 void FANGSStudyModule::initialize()
153 {
154  B2INFO("FANGSStudyModule: Initialize");
155 
156  REG_HISTOGRAM
157 
158  //convert sample time into rate in Hz
159  //rateCorrection = m_sampletime / 1e6;
160  //get FANGS paramters
161  getXMLData();
162 
163  fctQ_Calib1 = new TF1("fctQ_Calib1", "[0]*([1]*x-[2])/([3]-x)", 0., 15.);
164  fctQ_Calib1->SetParameters(m_TOTQ1, m_TOTC1, m_TOTA1 * m_TOTB1, m_TOTA1);
165 
166  fctQ_Calib2 = new TF1("fctQ_Calib2", "[0]*([1]*x-[2])/([3]-x)", 0., 15.);
167  fctQ_Calib2->SetParameters(m_TOTQ2, m_TOTC2, m_TOTA2 * m_TOTB2, m_TOTA2);
168 }
169 
170 void FANGSStudyModule::beginRun()
171 {
172 }
173 
174 void FANGSStudyModule::event()
175 {
176  //Here comes the actual event processing
177 
178  StoreArray<FANGSSimHit> SimHits;
179  StoreArray<FANGSHit> Hits;
180 
181  int olddetNb = -1;
182  int ipix = 0;
183  float esum = 0;
184  //number of entries in Hits
185  for (const auto& FANGSHit : Hits) {
186  int detNb = FANGSHit.getdetNb();
187  //int pdg = FANGSHit.getPDG();
188  //int trkID = FANGSHit.gettrkID();
189  int col = FANGSHit.getcolumn();
190  int row = FANGSHit.getrow();
191  int tot = FANGSHit.getTOT();
192  int bcid = FANGSHit.getBCID();
193 
194  if (olddetNb != detNb) {
195  if (esum > 0) {
196  h_Edep->Fill(detNb, esum);
197  h_pxNb->Fill(detNb, ipix);
198  }
199  ipix = 0;
200  esum = 0;
201  for (int j = 0; j < maxSIZE; j++) {
202  x[j] = 0;
203  y[j] = 0;
204  z[j] = 0;
205  e[j] = 0;
206  }
207  olddetNb = detNb;
208  }
209  x[ipix] = col * (2. * m_ChipColumnX / (float)m_ChipColumnNb) - m_ChipColumnX;
210  y[ipix] = row * (2. * m_ChipRowY / (float)m_ChipRowNb) - m_ChipRowY;
211  z[ipix] = (m_PixelTimeBin / 2. + m_PixelTimeBin * bcid) * m_v_sensor;
212  if (tot < 3) e[ipix] = fctQ_Calib1->Eval(tot) * m_Workfct * 1e-3;
213  else e[ipix] = fctQ_Calib2->Eval(tot) * m_Workfct * 1e-3;
214  esum += e[ipix];
215  h_cvr[detNb]->Fill(col, row);
216  ipix ++;
217  }
218  //number of entries in SimHits
219  int nSimHits = SimHits.getEntries();
220  //cout << nSimHits << endl;
221 
222  //loop over all SimHit entries
223  for (int i = 0; i < nSimHits; i++) {
224  FANGSSimHit* aHit = SimHits[i];
225  int lad = aHit->getLadder();
226  int sen = aHit->getSensor();
227  double adep = aHit->getEnergyDep();
228  double timeBin = aHit->getTime();
229  int pdg = aHit->getPDG();
230 
231  TVector3 position = aHit->getPosEntry();
232  double r = sqrt(position.X() * position.X() + position.Y() * position.Y());
233  int detNB = (lad - 1) * 5 + sen - 1;
234  //cout <<" lad " << lad << " sen " << sen << " detNB " << detNB << " time " << timeBin << " edep " << adep*1e3 << endl;
235  //cout <<" x " << position.X() << " y " << position.Y() << " z " << position.Z() << endl;
236  h_time->Fill(detNB, timeBin);
237  h_edep->Fill(detNB, adep * 1e3);
238  if (fabs(pdg) == Const::electron.getPDGCode())h_edep1->Fill(detNB, adep * 1e3);
239  if (pdg == Const::photon.getPDGCode())h_edep2->Fill(detNB, adep * 1e3);
240  if (pdg != Const::photon.getPDGCode() && fabs(pdg) != Const::electron.getPDGCode())h_edep3->Fill(detNB, adep * 1e3);
241  if (adep > 50.*1e-6) {
242  h_timeThres->Fill(detNB, timeBin);
243  h_edepThres->Fill(detNB, adep * 1e3);
244  if (fabs(pdg) == Const::electron.getPDGCode())h_edepThres1->Fill(detNB, adep * 1e3);
245  if (pdg == Const::photon.getPDGCode())h_edepThres2->Fill(detNB, adep * 1e3);
246  if (pdg != Const::photon.getPDGCode() && fabs(pdg) != Const::electron.getPDGCode())h_edepThres3->Fill(detNB, adep * 1e3);
247  }
248  h_zvedep[lad - 1]->Fill(position.Z());
249  h_xvzvedep[lad - 1]->Fill(position.X(), position.Z());
250  h_yvzvedep[lad - 1]->Fill(position.Y(), position.Z());
251  h_xvyvedep[lad - 1]->Fill(position.X(), position.Y());
252  h_rvzvedep[lad - 1]->Fill(r, position.Z());
253  h_zvedepW[lad - 1]->Fill(position.Z(), adep * 1e3);
254  h_xvzvedepW[lad - 1]->Fill(position.X(), position.Z(), adep * 1e3);
255  h_yvzvedepW[lad - 1]->Fill(position.Y(), position.Z(), adep * 1e3);
256  h_xvyvedepW[lad - 1]->Fill(position.X(), position.Y(), adep * 1e3);
257  h_rvzvedepW[lad - 1]->Fill(r, position.Z(), adep * 1e3);
258  if (adep > 50.*1e-6) {
259  h_zvedepT[lad - 1]->Fill(position.Z());
260  h_xvzvedepT[lad - 1]->Fill(position.X(), position.Z());
261  h_yvzvedepT[lad - 1]->Fill(position.Y(), position.Z());
262  h_xvyvedepT[lad - 1]->Fill(position.X(), position.Y());
263  h_rvzvedepT[lad - 1]->Fill(r, position.Z());
264  h_zvedepWT[lad - 1]->Fill(position.Z(), adep * 1e3);
265  h_xvzvedepWT[lad - 1]->Fill(position.X(), position.Z(), adep * 1e3);
266  h_yvzvedepWT[lad - 1]->Fill(position.Y(), position.Z(), adep * 1e3);
267  h_xvyvedepWT[lad - 1]->Fill(position.X(), position.Y(), adep * 1e3);
268  h_rvzvedepWT[lad - 1]->Fill(r, position.Z(), adep * 1e3);
269  }
270  }
271 
272  eventNum++;
273 }
274 
275 void FANGSStudyModule::getXMLData()
276 {
277  GearDir content = GearDir("/Detector/DetectorComponent[@name=\"FANGS\"]/Content/");
278 
279  m_PixelThreshold = content.getInt("PixelThreshold");
280  m_PixelThresholdRMS = content.getInt("PixelThresholdRMS");
281  m_PixelTimeBinNb = content.getInt("PixelTimeBinNb");
282  m_PixelTimeBin = content.getDouble("PixelTimeBin");
283  m_ChipColumnNb = content.getInt("ChipColumnNb");
284  m_ChipRowNb = content.getInt("ChipRowNb");
285  m_ChipColumnX = content.getDouble("ChipColumnX");
286  m_ChipRowY = content.getDouble("ChipRowY");
287  m_TOTA1 = content.getDouble("TOTA1");
288  m_TOTB1 = content.getDouble("TOTB1");
289  m_TOTC1 = content.getDouble("TOTC1");
290  m_TOTQ1 = content.getDouble("TOTQ1");
291  m_TOTA2 = content.getDouble("TOTA2");
292  m_TOTB2 = content.getDouble("TOTB2");
293  m_TOTC2 = content.getDouble("TOTC2");
294  m_TOTQ2 = content.getDouble("TOTQ2");
295  m_v_sensor = content.getDouble("v_sensor");
296  m_sensor_width = content.getDouble("sensor_width");
297  m_Workfct = content.getDouble("Workfct");
298 
299 }
300 
301 void FANGSStudyModule::endRun()
302 {
303 
304 
305 
306 }
307 
308 void FANGSStudyModule::terminate()
309 {
310 }
311 
312 
Belle2::FANGSHit
ClassFANGSHit - digitization simulated hit for the FANGS detector.
Definition: FANGSHit.h:26
Belle2::FANGSHit::getBCID
int getBCID() const
Return the BCID.
Definition: FANGSHit.h:54
Belle2::FANGSHit::getrow
int getrow() const
Return the row.
Definition: FANGSHit.h:52
Belle2::FANGSHit::getTOT
int getTOT() const
Return the TOT.
Definition: FANGSHit.h:56
Belle2::FANGSHit::getdetNb
int getdetNb() const
Return the TPC number.
Definition: FANGSHit.h:58
Belle2::FANGSHit::getcolumn
int getcolumn() const
Return the column.
Definition: FANGSHit.h:50
Belle2::FANGSSimHit
Class FANGSSimHit - Geant4 simulated hit for the FANGS detector.
Definition: FANGSSimHit.h:28
Belle2::FANGSSimHit::getTime
float getTime() const
Return the global time.
Definition: FANGSSimHit.h:74
Belle2::FANGSSimHit::getEnergyDep
float getEnergyDep() const
Return the energy deposition in electrons.
Definition: FANGSSimHit.h:76
Belle2::FANGSSimHit::getLadder
int getLadder() const
Return the Ladder number (starting at 1, increasing with phi)
Definition: FANGSSimHit.h:68
Belle2::FANGSSimHit::getSensor
int getSensor() const
Return the Sensor number (starting at 1, increasing with decreasing z)
Definition: FANGSSimHit.h:70
Belle2::FANGSSimHit::getPDG
int getPDG() const
Return the PDG number of the track.
Definition: FANGSSimHit.h:72
Belle2::FANGSSimHit::getPosEntry
TVector3 getPosEntry() const
Return the entry track position.
Definition: FANGSSimHit.h:78
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:31
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
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
Belle2::fangs::FANGSStudyModule
Study module for Fangs (BEAST)
Definition: FANGSStudyModule.h:33
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