10 #include <svd/modules/svdSpacePointCreator/SpacePointHelperFunctions.h>
11 #include <svd/modules/svdSpacePointCreator/SVDSpacePointQICalibrationModule.h>
26 setDescription(
"Generate PDFs used in assigning quality to SpacePoints.");
27 setPropertyFlags(c_ParallelProcessingCertified);
30 addParam(
"SVDClusters", m_svdClustersName,
31 "SVDCluster collection name",
string(
""));
33 addParam(
"RecoTracks", m_recoTracksName,
34 "RecoTracks collection name",
string(
""));
37 addParam(
"NameOfInstance", m_nameOfInstance,
38 "allows the user to set an identifier for this module. Usefull if one wants to use several instances of that module",
string(
""));
39 addParam(
"binSize", m_binSize,
"Number of bins in charge distribution.",
42 addParam(
"maxClusterSize", m_maxClusterSize,
"Max number of strips the PDF are separated into.",
45 addParam(
"useLegacyNaming", m_useLegacyNaming,
46 "Use legacy pdf naming",
bool(
true));
48 addParam(
"outputFileName", m_outputFileName,
49 "Name of output file containing pdfs", std::string(
"spacePointQICalibration.root"));
54 void SVDSpacePointQICalibrationModule::initialize()
57 m_svdClusters.isRequired(m_svdClustersName);
58 m_recoTracks.isRequired(m_recoTracksName);
63 for (
auto& layers : geo.
getLayers(VXD::SensorInfoBase::SVD)) {
65 for (
auto& sensors : geo.
getSensors(ladders)) {
67 for (
int uSize = 1; uSize <= m_maxClusterSize; uSize++) {
68 for (
int vSize = 1; vSize <= m_maxClusterSize; vSize++) {
69 std::string sensorName =
"";
70 std::string errorStringName =
"";
73 spPDFName(sensors, uSize, vSize, m_maxClusterSize, sensorName, errorStringName, m_useLegacyNaming);
75 std::string signalHistName = sensorName +
"signal";
76 std::string backgroundHistName = sensorName +
"background";
78 if (signalHistMap.count(sensorName) == 0) {
79 TH2F* sigHist =
new TH2F(signalHistName.c_str(),
"", m_binSize, 0, 200000, m_binSize, 0, 200000);
80 TH2F* bkgHist =
new TH2F(backgroundHistName.c_str(),
"", m_binSize, 0, 200000, m_binSize, 0, 200000);
82 signalHistMap[sensorName] = sigHist;
83 backgroundHistMap[sensorName] = bkgHist;
90 TH2F* timeSignal =
new TH2F(
"timeSignal",
"", 40, -100, 100, 40, -100, 100);
91 TH2F* timeBackground =
new TH2F(
"timeBackground",
"", 40, -100, 100, 40, -100, 100);
92 TH2F* sizeSignal =
new TH2F(
"sizeSignal",
"", 20, 0, 20, 20, 0, 20);
93 TH2F* sizeBackground =
new TH2F(
"sizeBackground",
"", 20, 0, 20, 20, 0, 20);
95 signalHistMap[
"timeSignal"] = timeSignal;
96 signalHistMap[
"sizeSignal"] = sizeSignal;
97 backgroundHistMap[
"timeBackground"] = timeBackground;
98 backgroundHistMap[
"sizeBackground"] = sizeBackground;
104 void SVDSpacePointQICalibrationModule::event()
106 for (
auto& track : m_recoTracks) {
107 if (track.wasFitSuccessful() == 1) {
108 if (track.hasSVDHits() == 1) {
109 for (
auto& svdHit : track.getSVDHitList()) {
110 if (svdHit->isUCluster()) {
111 for (
auto& svdHit2 : track.getSVDHitList()) {
112 if (svdHit2->isUCluster() == 0) {
113 if (svdHit->getSensorID() == svdHit2->getSensorID()) {
115 std::string errorStringName;
116 spPDFName(svdHit->getSensorID(), svdHit->getSize(), svdHit2->getSize(), m_maxClusterSize, pdfName, errorStringName,
118 auto sigHist = signalHistMap.at(pdfName);
119 sigHist->Fill(svdHit->getCharge(), svdHit2->getCharge());
120 auto timeSigHist = signalHistMap.at(
"timeSignal");
121 timeSigHist->Fill(svdHit->getClsTime(), svdHit2->getClsTime());
122 auto sizeSigHist = signalHistMap.at(
"sizeSignal");
123 sizeSigHist->Fill(svdHit->getSize(), svdHit2->getSize());
134 for (
auto& uCluster : m_svdClusters) {
135 if (uCluster.isUCluster() == 1) {
136 for (
auto& vCluster : m_svdClusters) {
137 if (vCluster.isUCluster() == 0) {
138 if (uCluster.getSensorID() == vCluster.getSensorID()) {
140 std::string errorStringName;
141 spPDFName(uCluster.getSensorID(), uCluster.getSize(), vCluster.getSize(), m_maxClusterSize, pdfName, errorStringName,
143 auto bkgHist = backgroundHistMap.at(pdfName);
144 bkgHist->Fill(uCluster.getCharge(), vCluster.getCharge());
145 auto timeBkgHist = backgroundHistMap.at(
"timeBackground");
146 timeBkgHist->Fill(uCluster.getClsTime(), vCluster.getClsTime());
147 auto sizeBkgHist = backgroundHistMap.at(
"sizeBackground");
148 sizeBkgHist->Fill(uCluster.getSize(), vCluster.getSize());
158 void SVDSpacePointQICalibrationModule::terminate()
161 TFile* f =
new TFile(m_outputFileName.c_str(),
"RECREATE");
163 std::vector<std::string> usedSensors;
165 for (
auto& layers : geo.
getLayers(VXD::SensorInfoBase::SVD)) {
166 for (
auto& ladders : geo.
getLadders(layers)) {
167 for (
auto& sensors : geo.
getSensors(ladders)) {
168 for (
int uSize = 1; uSize <= m_maxClusterSize; uSize++) {
169 for (
int vSize = 1; vSize <= m_maxClusterSize; vSize++) {
170 std::string sensorName;
171 std::string errorName;
174 spPDFName(sensors, uSize, vSize, m_maxClusterSize, sensorName, errorName, m_useLegacyNaming);
176 if (std::find(usedSensors.begin(), usedSensors.end(), sensorName.c_str()) == usedSensors.end()) {
177 usedSensors.push_back(sensorName);
178 TH2F* probHist =
new TH2F(sensorName.c_str(),
"", m_binSize, 0, 200000, m_binSize, 0, 200000);
179 TH2F* errorHist =
new TH2F(errorName.c_str(),
"", m_binSize, 0, 200000, m_binSize, 0, 200000);
180 calculateProb(signalHistMap[sensorName], backgroundHistMap[sensorName], probHist);
181 calculateError(signalHistMap[sensorName], backgroundHistMap[sensorName], errorHist);
191 TH2F* timeProb =
new TH2F(
"timeProb",
"", 40, -100, 100, 40, -100, 100);
192 TH2F* timeError =
new TH2F(
"timeError",
"", 40, -100, 100, 40, -100, 100);
193 TH2F* sizeProb =
new TH2F(
"sizeProb",
"", 20, 0, 20, 20, 0, 20);
194 TH2F* sizeError =
new TH2F(
"sizeError",
"", 20, 0, 20, 20, 0, 20);
195 calculateProb(signalHistMap[
"timeSignal"], backgroundHistMap[
"timeBackground"], timeProb);
196 calculateError(signalHistMap[
"timeSignal"], backgroundHistMap[
"timeBackground"], timeError);
197 calculateProb(signalHistMap[
"sizeSignal"], backgroundHistMap[
"sizeBackground"], sizeProb);
198 calculateError(signalHistMap[
"sizeSignal"], backgroundHistMap[
"sizeBackground"], sizeError);
208 void SVDSpacePointQICalibrationModule::calculateProb(TH2F* signal, TH2F* background, TH2F* probability)
210 signal->Divide(background);
211 probability->Add(signal);
214 void SVDSpacePointQICalibrationModule::calculateError(TH2F* signal, TH2F* background, TH2F* error)
216 int imax = signal->GetXaxis()->GetNbins();
217 int jmax = signal->GetYaxis()->GetNbins();
218 for (
int i = 1; i <= imax; i++) {
219 for (
int j = 1; j <= jmax; j++) {
220 int bkg = background->GetBinContent(i, j);
221 int sig = signal->GetBinContent(i, j);
222 double var = ((sig + 1) * (sig + 2)) / ((bkg + 2) * (bkg + 3)) -
223 ((sig + 1) * (sig + 1)) / ((bkg + 2) * (bkg + 2));
224 double err = sqrt(var);
225 error->SetBinContent(i, j, err);