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