Belle II Software  release-05-01-25
PXDClusterChargeCollectorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Benjamin Schwenker *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <pxd/modules/pxdClusterChargeCollector/PXDClusterChargeCollectorModule.h>
12 
13 #include <framework/datastore/RelationArray.h>
14 
15 #include <vxd/geometry/GeoCache.h>
16 #include <pxd/geometry/SensorInfo.h>
17 #include <pxd/reconstruction/PXDGainCalibrator.h>
18 
19 #include <TTree.h>
20 #include <TH2I.h>
21 
22 #include <boost/format.hpp>
23 
24 using namespace std;
25 using boost::format;
26 using namespace Belle2;
27 
28 //-----------------------------------------------------------------
29 // Register the Module
30 //-----------------------------------------------------------------
31 REG_MODULE(PXDClusterChargeCollector)
32 
33 //-----------------------------------------------------------------
34 // Implementation
35 //-----------------------------------------------------------------
36 
38  , m_signal(0), m_run(0), m_exp(0)
39 {
40  // Set module properties
41  setDescription("Calibration collector module for cluster charge related calibrations.");
42  setPropertyFlags(c_ParallelProcessingCertified);
43 
44  addParam("clustersName", m_storeClustersName, "Name of the collection to use for PXDClusters", string(""));
45  addParam("mcParticlesName", m_storeMCParticlesName, "Name of the collection to use for MCParticles", string("MCParticles"));
46  addParam("minClusterCharge", m_minClusterCharge, "Minimum cluster charge cut", int(0));
47  addParam("minClusterSize", m_minClusterSize, "Minimum cluster size cut ", int(2));
48  addParam("maxClusterSize", m_maxClusterSize, "Maximum cluster size cut ", int(6));
49  addParam("nBinsU", m_nBinsU, "Number of gain corrections per sensor along u side", int(4));
50  addParam("nBinsV", m_nBinsV, "Number of gain corrections per sensor along v side", int(6));
51  addParam("chargePayloadName", m_chargeName, "Payload name for Cluster Charge to be read from DB", string(""));
52  addParam("gainPayloadName", m_gainName, "Payload name for Gain to be read from DB", string(""));
53  addParam("fillChargeHistogram", m_fillChargeHistogram, "Flag to fill Charge histograms", bool(false));
54  addParam("matchTrack", m_matchTrack,
55  "Flag to use track matched clusters (=1) and apply theta angle projection to cluster charge (=2)", int(0));
56  addParam("mcSamples", m_mcSamples, "Flag to deal with MC samples", bool(false));
57  addParam("relationCheck", m_relationCheck, "Flag to check relations between PXDClusters and PXDClustersFromTracks", bool(false));
58 
59 }
60 
61 void PXDClusterChargeCollectorModule::prepare() // Do your initialise() stuff here
62 {
63  m_pxdClusters.isRequired(m_storeClustersName);
64  if (m_relationCheck) {
65  std::string storeClustersName2;
66  storeClustersName2 = (m_storeClustersName == "PXDClusters") ?
67  "PXDClustersFromTracks" : "PXDClusters";
68  StoreArray<PXDCluster> pxdClusters2;
69  pxdClusters2.isRequired(storeClustersName2);
70  }
71  m_tracks.isOptional(); //m_storeTracksName);
72  m_recoTracks.isOptional(); //m_storeRecoTracksName);
73 
74  if (m_mcSamples && m_matchTrack > 0) {
75  m_mcParticles.isRequired(m_storeMCParticlesName);
76  } else {
77  m_mcParticles.isOptional();
78  }
79 
80  if (m_nBinsU == 0) {
81  B2WARNING("Number of bins along u side incremented from 0->1");
82  m_nBinsU = 1;
83  }
84 
85  if (m_nBinsV == 0) {
86  B2WARNING("Number of bins along v side incremented from 0->1");
87  m_nBinsV = 1;
88  }
89 
90  auto gTools = VXD::GeoCache::getInstance().getGeoTools();
91  int nPXDSensors = gTools->getNumberOfPXDSensors();
92 
93  //-------------------------------------------------------------------------------------
94  // PXDClusterCounter: Count the number of PXDClusters and store charge for each uBin/vBin pair
95  //-------------------------------------------------------------------------------------
96 
97  auto hPXDClusterCounter = new TH1I("hPXDClusterCounter", "Number of clusters found in data sample",
98  m_nBinsU * m_nBinsV * nPXDSensors, 0,
99  m_nBinsU * m_nBinsV * nPXDSensors);
100  hPXDClusterCounter->GetXaxis()->SetTitle("bin id");
101  hPXDClusterCounter->GetYaxis()->SetTitle("Number of clusters");
102  for (int iSensor = 0; iSensor < nPXDSensors; iSensor++) {
103  for (int uBin = 0; uBin < m_nBinsU; uBin++) {
104  for (int vBin = 0; vBin < m_nBinsV; vBin++) {
105  VxdID id = gTools->getSensorIDFromPXDIndex(iSensor);
106  string sensorDescr = id;
107  hPXDClusterCounter->GetXaxis()->SetBinLabel(iSensor * m_nBinsU * m_nBinsV + uBin * m_nBinsV + vBin + 1,
108  str(format("%1%_%2%_%3%") % sensorDescr % uBin % vBin).c_str());
109  }
110  }
111  }
112  registerObject<TH1I>("PXDClusterCounter", hPXDClusterCounter);
113  if (m_fillChargeHistogram) {
114  auto hPXDClusterCharge = new TH2I("hPXDClusterCharge", "Charge of clusters found in data sample",
115  m_nBinsU * m_nBinsV * nPXDSensors, 0,
116  m_nBinsU * m_nBinsV * nPXDSensors,
117  250, 0, 250);
118  hPXDClusterCharge->GetXaxis()->SetTitle("bin id");
119  hPXDClusterCharge->GetYaxis()->SetTitle("Cluster charge [ADU]");
120  registerObject<TH2I>("PXDClusterCharge", hPXDClusterCharge);
121  }
122 
123  //----------------------------------------------------------------------
124  // PXDTrees: One tree to store the calibration data for each grid bin
125  //----------------------------------------------------------------------
126 
127  for (int iSensor = 0; iSensor < nPXDSensors; iSensor++) {
128  for (int uBin = 0; uBin < m_nBinsU; uBin++) {
129  for (int vBin = 0; vBin < m_nBinsV; vBin++) {
130  VxdID id = gTools->getSensorIDFromPXDIndex(iSensor);
131  auto layerNumber = id.getLayerNumber();
132  auto ladderNumber = id.getLadderNumber();
133  auto sensorNumber = id.getSensorNumber();
134  string treename = str(format("tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
135  auto tree = new TTree(treename.c_str(), treename.c_str());
136  tree->Branch<int>("signal", &m_signal);
137  registerObject<TTree>(treename, tree);
138  }
139  }
140  }
141 
142  auto dbtree = new TTree("dbtree", "dbtree");
143  dbtree->Branch<int>("run", &m_run);
144  dbtree->Branch<int>("exp", &m_exp);
145  dbtree->Branch<PXDClusterChargeMapPar>("chargeMap", &m_chargeMap);
146  dbtree->Branch<PXDGainMapPar>("gainMap", &m_gainMap);
147  registerObject<TTree>("dbtree", dbtree);
148 }
149 
150 void PXDClusterChargeCollectorModule::startRun() // Do your beginRun() stuff here
151 {
152  m_run = m_evtMetaData->getRun();
153  m_exp = m_evtMetaData->getExperiment();
154  if (m_chargeName.length()) {
155  DBObjPtr<PXDClusterChargeMapPar> chargeMap(m_chargeName);
156  m_chargeMap = *chargeMap;
157  } else {
158  m_chargeMap = PXDClusterChargeMapPar();
159  }
160  if (m_gainName.length()) {
161  DBObjPtr<PXDGainMapPar> gainMap(m_gainName);
162  m_gainMap = *gainMap;
163  } else {
164  m_gainMap = PXDGainMapPar();
165  }
166  getObjectPtr<TTree>("dbtree")->Fill();
167 }
168 
169 void PXDClusterChargeCollectorModule::collect() // Do your event() stuff here
170 {
171  // If no input, nothing to do
172  if (!m_pxdClusters) return;
173 
174  auto gTools = VXD::GeoCache::getInstance().getGeoTools();
175 
176  if (!m_matchTrack || m_mcSamples) { // all clusters for data
177 
178  for (auto& cluster : m_pxdClusters) {
179  if (m_relationCheck) {
180  if (m_storeClustersName == "PXDClusters") {
181  RelationVector<PXDCluster> relPXDClusters = cluster.getRelationsTo<PXDCluster>("PXDClustersFromTracks");
182  if (!relPXDClusters.size()) continue;
183  } else if (m_storeClustersName == "PXDClustersFromTracks") {
184  RelationVector<PXDCluster> relPXDClusters = DataStore::getRelationsWithObj<PXDCluster>(&cluster, "PXDClusters");
185  if (!relPXDClusters.size()) continue;
186  }
187  }
188  double correction = 1.0;
189  if (m_matchTrack > 0) { // mc samples
190  RelationVector<MCParticle> mcParticlesPXD = DataStore::getRelationsWithObj<MCParticle>(&cluster, m_storeMCParticlesName);
191  if (!mcParticlesPXD.size()) continue;
192  if (m_matchTrack == 2) {
193  correction = sin(mcParticlesPXD[0]->getMomentum().Theta());
194  }
195  }
196  // Apply cluster selection cuts
197  if (cluster.getCharge() >= m_minClusterCharge && cluster.getSize() >= m_minClusterSize && cluster.getSize() <= m_maxClusterSize) {
198 
199  VxdID sensorID = cluster.getSensorID();
200  const PXD::SensorInfo& Info = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(sensorID));
201  auto uID = Info.getUCellID(cluster.getU());
202  auto vID = Info.getVCellID(cluster.getV());
203  auto iSensor = gTools->getPXDSensorIndex(sensorID);
204  auto layerNumber = sensorID.getLayerNumber();
205  auto ladderNumber = sensorID.getLadderNumber();
206  auto sensorNumber = sensorID.getSensorNumber();
207  auto uBin = PXD::PXDGainCalibrator::getInstance().getBinU(sensorID, uID, vID, m_nBinsU);
208  auto vBin = PXD::PXDGainCalibrator::getInstance().getBinV(sensorID, vID, m_nBinsV);
209  string treename = str(format("tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
210 
211  // Compute variables from cluster needed for gain estimation
212  m_signal = cluster.getCharge();
213  // Fill variabels into tree
214  getObjectPtr<TTree>(treename)->Fill();
215  // Increment the counter & store charge (optional)
216  getObjectPtr<TH1I>("PXDClusterCounter")->Fill(iSensor * m_nBinsU * m_nBinsV + uBin * m_nBinsV + vBin);
217  if (m_fillChargeHistogram)
218  getObjectPtr<TH2I>("PXDClusterCharge")->Fill(iSensor * m_nBinsU * m_nBinsV + uBin * m_nBinsV + vBin,
219  cluster.getCharge()*correction);
220  }
221 
222  }
223 
224  } else { // store only clusters matched with tracks
225 
226  // If no input, nothing to do
227  if (!m_tracks) return;
228 
229  for (const Track& track : m_tracks) {
230 
231  RelationVector<RecoTrack> recoTrack = track.getRelationsTo<RecoTrack>(m_storeRecoTracksName);
232  if (!recoTrack.size()) continue;
233  RelationVector<PXDCluster> pxdClustersTrack = DataStore::getRelationsWithObj<PXDCluster>(recoTrack[0], m_storeClustersName);
234 
235  const TrackFitResult* tfr = track.getTrackFitResultWithClosestMass(Const::pion);
236  double correction = 1.0;
237  if (tfr && m_matchTrack == 2) correction = sin(tfr->getMomentum().Theta());
238 
239  for (auto& cluster : pxdClustersTrack) {
240  if (m_relationCheck) {
241  if (m_storeClustersName == "PXDClustersFromTracks") {
242  RelationVector<PXDCluster> relPXDClusters = DataStore::getRelationsWithObj<PXDCluster>(&cluster, "PXDClusters");
243  if (!relPXDClusters.size()) continue;
244  } else {
245  B2WARNING("Relation check for data only works when clustersName is set to PXDClustersFromTracks");
246  }
247  }
248 
249  // Apply cluster selection cuts
250  if (cluster.getCharge() >= m_minClusterCharge && cluster.getSize() >= m_minClusterSize && cluster.getSize() <= m_maxClusterSize) {
251 
252  VxdID sensorID = cluster.getSensorID();
253  const PXD::SensorInfo& Info = dynamic_cast<const PXD::SensorInfo&>(VXD::GeoCache::get(sensorID));
254  auto uID = Info.getUCellID(cluster.getU());
255  auto vID = Info.getVCellID(cluster.getV());
256  auto iSensor = gTools->getPXDSensorIndex(sensorID);
257  auto layerNumber = sensorID.getLayerNumber();
258  auto ladderNumber = sensorID.getLadderNumber();
259  auto sensorNumber = sensorID.getSensorNumber();
260  auto uBin = PXD::PXDGainCalibrator::getInstance().getBinU(sensorID, uID, vID, m_nBinsU);
261  auto vBin = PXD::PXDGainCalibrator::getInstance().getBinV(sensorID, vID, m_nBinsV);
262  string treename = str(format("tree_%1%_%2%_%3%_%4%_%5%") % layerNumber % ladderNumber % sensorNumber % uBin % vBin);
263 
264  // Compute variables from cluster needed for gain estimation
265  m_signal = cluster.getCharge();
266  // Fill variabels into tree
267  getObjectPtr<TTree>(treename)->Fill();
268  // Increment the counter & store charge (optional)
269  getObjectPtr<TH1I>("PXDClusterCounter")->Fill(iSensor * m_nBinsU * m_nBinsV + uBin * m_nBinsV + vBin);
270  if (m_fillChargeHistogram)
271  getObjectPtr<TH2I>("PXDClusterCharge")->Fill(iSensor * m_nBinsU * m_nBinsV + uBin * m_nBinsV + vBin,
272  cluster.getCharge()*correction);
273 
274  }
275 
276  }
277 
278  }
279 
280  }
281 
282 }
Belle2::CalibrationCollectorModule
Calibration collector module base class.
Definition: CalibrationCollectorModule.h:44
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::VXD::SensorInfoBase::getUCellID
int getUCellID(double u, double v=0, bool clamp=false) const
Return the corresponding pixel/strip ID of a given u coordinate.
Definition: SensorInfoBase.h:203
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:116
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::PXDGainMapPar
The payload class for PXD gain corrections.
Definition: PXDGainMapPar.h:53
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::PXD::SensorInfo
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Definition: SensorInfo.h:34
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::VXD::SensorInfoBase::getVCellID
int getVCellID(double v, bool clamp=false) const
Return the corresponding pixel/strip ID of a given v coordinate.
Definition: SensorInfoBase.h:213
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::PXDClusterChargeMapPar
The payload class for PXD cluster charge calibrations.
Definition: PXDClusterChargeMapPar.h:40
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::StoreArray< PXDCluster >
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::PXDClusterChargeCollectorModule
Collector module for PXD gain calibration.
Definition: PXDClusterChargeCollectorModule.h:52