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