12 #include <klm/bklm/modules/bklmAna/BKLMAnaModule.h>
15 #include <framework/dataobjects/EventMetaData.h>
16 #include <framework/datastore/StoreObjPtr.h>
17 #include <mdst/dataobjects/Track.h>
18 #include <mdst/dataobjects/TrackFitResult.h>
19 #include <klm/dataobjects/KLMMuidLikelihood.h>
22 #include <CLHEP/Units/SystemOfUnits.h>
26 using namespace CLHEP;
53 for (
int i = 0; i < 200; ++i) {
58 for (
int i = 0; i < 15; ++i) {
59 m_totalThephi[i] =
nullptr;
60 m_passThephi[i] =
nullptr;
61 m_effiThephi[i] =
nullptr;
62 m_totalTrkThephi[i] =
nullptr;
63 m_passTrkThephi[i] =
nullptr;
64 m_effiTrkThephi[i] =
nullptr;
66 setDescription(
"analyze bklm efficiency associated to CDC, check performance of bklm et al.");
67 addParam(
"filename", m_filename,
"Output root filename",
string(
"bklmana.root"));
70 BKLMAnaModule::~BKLMAnaModule()
74 void BKLMAnaModule::initialize()
81 m_file =
new TFile(m_filename.c_str(),
"RECREATE");
83 m_extTree =
new TTree(
"exthit",
"ext hit");
84 m_extTree->Branch(
"run", &m_run,
"m_run/I");
85 m_extTree->Branch(
"nExtHit", &m_nExtHit,
"m_nExtHit/I");
86 m_extTree->Branch(
"x", &m_extx,
"m_extx[m_nExtHit]/F");
87 m_extTree->Branch(
"y", &m_exty,
"m_exty[m_nExtHit]/F");
88 m_extTree->Branch(
"z", &m_extz,
"m_extz[m_nExtHit]/F");
94 float thetaMin = 35.0;
95 float thetaMax = 130.0;
97 for (
int iL = 0; iL < 15 ; iL ++) {
99 sprintf(hname,
"denominator_Layer%i", iL + 1);
100 m_totalThephi[iL] =
new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
101 m_totalThephi[iL]->GetXaxis()->SetTitle(
"phi");
102 m_totalThephi[iL]->GetYaxis()->SetTitle(
"theta");
103 sprintf(hname,
"numerator_Layer%i", iL + 1);
104 m_passThephi[iL] =
new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
105 sprintf(hname,
"effi_Layer%i", iL + 1);
106 m_passThephi[iL]->GetXaxis()->SetTitle(
"phi");
107 m_passThephi[iL]->GetYaxis()->SetTitle(
"theta");
108 m_effiThephi[iL] =
new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
109 m_effiThephi[iL]->GetXaxis()->SetTitle(
"phi");
110 m_effiThephi[iL]->GetYaxis()->SetTitle(
"theta");
112 sprintf(hname,
"Denominator_Layer%i", iL + 1);
113 m_totalTrkThephi[iL] =
new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
114 m_totalTrkThephi[iL]->GetXaxis()->SetTitle(
"trk.phi");
115 m_totalTrkThephi[iL]->GetYaxis()->SetTitle(
"trk.theta");
116 sprintf(hname,
"Numerator_Layer%i", iL + 1);
117 m_passTrkThephi[iL] =
new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
118 sprintf(hname,
"Effi_Layer%i", iL + 1);
119 m_passTrkThephi[iL]->GetXaxis()->SetTitle(
"trk.phi");
120 m_passTrkThephi[iL]->GetYaxis()->SetTitle(
"trk.theta");
121 m_effiTrkThephi[iL] =
new TH2F(hname, hname, phiBins, phiMin, phiMax, thetaBins, thetaMin, thetaMax);
122 m_effiTrkThephi[iL]->GetXaxis()->SetTitle(
"trk.phi");
123 m_effiTrkThephi[iL]->GetYaxis()->SetTitle(
"trk.theta");
132 m_hdistance =
new TH1F(
"m_hdistance",
" distance between mathced extHit and bklmHit2d ", 100, 0, 30);
134 m_totalMom =
new TH1F(
"m_totalMom",
" denominator vs. p", mNbin, mommin, mommax);
135 m_passMom =
new TH1F(
"m_passMom",
" numerator vs. p", mNbin, mommin, mommax);
136 m_effiMom =
new TH1F(
"m_effiMom",
" effi. vs. p", mNbin, mommin, mommax);
137 m_effiMom->GetXaxis()->SetTitle(
"p (GeV)");
138 m_effiMom->GetYaxis()->SetTitle(
"Efficiency (GeV)");
140 m_totalYX =
new TH2F(
"m_totalYX",
" denominator Y vs. X", gNbin, gmin, gmax, gNbin, gmin, gmax);
141 m_passYX =
new TH2F(
"m_passYX",
" numerator Y vs. X", gNbin, gmin, gmax, gNbin, gmin, gmax);
142 m_totalYZ =
new TH2F(
"m_totalYZ",
" denominator Y vs. Z", gNbin, gmin, gmax, gNbin, gmin, gmax);
143 m_passYZ =
new TH2F(
"m_passYZ",
" numerator Y vs. Z", gNbin, gmin, gmax, gNbin, gmin, gmax);
144 m_effiYX =
new TH2F(
"m_effiYX",
" effi. Y vs. X", gNbin, gmin, gmax, gNbin, gmin, gmax);
145 m_effiYZ =
new TH2F(
"m_effiYZ",
" effi. Y vs. Z", gNbin, gmin, gmax, gNbin, gmin, gmax);
146 m_effiYX->GetXaxis()->SetTitle(
"x (cm)");
147 m_effiYX->GetYaxis()->SetTitle(
"y (cm)");
148 m_effiYZ->GetXaxis()->SetTitle(
"z (cm)");
149 m_effiYZ->GetYaxis()->SetTitle(
"y (cm)");
153 void BKLMAnaModule::beginRun()
157 void BKLMAnaModule::event()
161 unsigned long runNumber = eventMetaData->getRun();
176 for (
int t = 0; t < extHits.getEntries(); t++) {
177 ExtHit* exthit = extHits[t];
178 if (exthit->
getDetectorID() != Const::EDetector::BKLM)
continue;
183 if (nExtHit > 199)
break;
189 for (
int k = 0; k < tracks.getEntries(); k++) {
190 Track* track = tracks[k];
197 double trkphi = p4.Vect().Phi() * 180.0 / CLHEP::pi;
198 double trktheta = p4.Vect().Theta() * 180.0 / CLHEP::pi;
199 if (trkphi < 0) trkphi = trkphi + 360.0;
202 for (
unsigned int t = 0; t < relatedExtHit.size(); t++) {
203 ExtHit* exthit = relatedExtHit[t];
204 if (exthit->
getDetectorID() != Const::EDetector::BKLM)
continue;
206 int section = BKLMElementNumbers::getSectionByModule(module);
207 int sector = BKLMElementNumbers::getSectorByModule(module);
208 int layer = BKLMElementNumbers::getLayerByModule(module);
209 bool crossed =
false;
212 int extPattern = muid->getExtLayerPattern();
213 if ((extPattern & (1 << (layer - 1))) != 0) crossed =
true;
214 if (!crossed)
continue;
218 bool matched =
false;
219 m_totalYX->Fill(extVec[0], extVec[1]);
220 m_totalYZ->Fill(extVec[2], extVec[1]);
221 float phi = extVec.Phi() * 180.0 / CLHEP::pi;
222 float theta = extVec.Theta() * 180.0 / CLHEP::pi;
223 if (phi < 0) phi = phi + 360.0;
224 m_totalThephi[layer - 1]->Fill(phi, theta);
225 m_totalTrkThephi[layer - 1]->Fill(trkphi, trktheta);
226 m_totalMom->Fill(mom);
230 for (
int mHit = 0; mHit < hits2D.getEntries(); mHit++) {
233 if (hit->getSection() != section)
continue;
234 if (hit->getSector() != sector)
continue;
235 if (hit->getLayer() != layer)
continue;
236 TVector3 position = hit->getGlobalPosition();
237 TVector3 distance = extVec - position;
240 m_hdistance->Fill(distance.Mag());
241 if (distance.Mag() < 20) matched =
true;
246 m_passYX->Fill(extVec[0], extVec[1]);
247 m_passYZ->Fill(extVec[2], extVec[1]);
248 m_passTrkThephi[layer - 1]->Fill(trkphi, trktheta);
249 m_passThephi[layer - 1]->Fill(phi, theta);
250 m_passMom->Fill(mom);
257 void BKLMAnaModule::endRun()
261 void BKLMAnaModule::terminate()
263 for (
int iL = 0; iL < 15; iL ++) {
264 for (
int i = 0; i < m_totalThephi[iL]->GetNbinsX(); i++) {
265 for (
int j = 0; j < m_totalThephi[iL]->GetNbinsY(); j++) {
266 float num = m_passThephi[iL]->GetBinContent(i + 1, j + 1);
267 float denom = m_totalThephi[iL]->GetBinContent(i + 1, j + 1);
269 m_effiThephi[iL]->SetBinContent(i + 1, j + 1, num / denom);
270 m_effiThephi[iL]->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
272 m_effiThephi[iL]->SetBinContent(i + 1, j + 1, 0);
273 m_effiThephi[iL]->SetBinError(i + 1, j + 1, 0);
279 for (
int iL = 0; iL < 15; iL ++) {
280 for (
int i = 0; i < m_totalTrkThephi[iL]->GetNbinsX(); i++) {
281 for (
int j = 0; j < m_totalTrkThephi[iL]->GetNbinsY(); j++) {
282 float num = m_passTrkThephi[iL]->GetBinContent(i + 1, j + 1);
283 float denom = m_totalTrkThephi[iL]->GetBinContent(i + 1, j + 1);
285 m_effiTrkThephi[iL]->SetBinContent(i + 1, j + 1, num / denom);
286 m_effiTrkThephi[iL]->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
288 m_effiTrkThephi[iL]->SetBinContent(i + 1, j + 1, 0);
289 m_effiTrkThephi[iL]->SetBinError(i + 1, j + 1, 0);
296 for (
int i = 0; i < m_totalYX->GetNbinsX(); i++) {
297 for (
int j = 0; j < m_totalYX->GetNbinsY(); j++) {
298 float num = m_passYX->GetBinContent(i + 1, j + 1);
299 float denom = m_totalYX->GetBinContent(i + 1, j + 1);
301 m_effiYX->SetBinContent(i + 1, j + 1, num / denom);
302 m_effiYX->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
304 m_effiYX->SetBinContent(i + 1, j + 1, 0);
305 m_effiYX->SetBinError(i + 1, j + 1, 0);
308 num = m_passYZ->GetBinContent(i + 1, j + 1);
309 denom = m_totalYZ->GetBinContent(i + 1, j + 1);
311 m_effiYZ->SetBinContent(i + 1, j + 1, num / denom);
312 m_effiYZ->SetBinError(i + 1, j + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
314 m_effiYZ->SetBinContent(i + 1, j + 1, 0);
315 m_effiYZ->SetBinError(i + 1, j + 1, 0);
320 for (
int i = 0; i < m_totalMom->GetNbinsX(); i++) {
321 float num = m_passMom->GetBinContent(i + 1);
322 float denom = m_totalMom->GetBinContent(i + 1);
324 m_effiMom->SetBinContent(i + 1, num / denom);
325 m_effiMom->SetBinError(i + 1, sqrt(num * (denom - num) / (denom * denom * denom)));
327 m_effiMom->SetBinContent(i + 1, 0);
328 m_effiMom->SetBinError(i + 1, 0);
333 m_hdistance->Write();
343 for (
int i = 0; i < 15; i++) {
344 m_totalThephi[i]->Write();
345 m_passThephi[i]->Write();
346 m_effiThephi[i]->Write();
347 m_totalTrkThephi[i]->Write();
348 m_passTrkThephi[i]->Write();
349 m_effiTrkThephi[i]->Write();