10 #include "TClonesArray.h"
26 #include "CbmTofHit.h"
29 #include "FairDetector.h"
30 #include "FairLogger.h"
31 #include "FairMCPoint.h"
56 std::vector<std::pair<std::string, std::array<int, 4>>> decNames,
59 : FairTask(
"CbmRecoQa")
60 , pullresfile(nullptr)
61 , verbosity(verbose_l)
64 std::vector<std::vector<TH1F*>>(detectors.size(), std::vector<TH1F*>(6)))
80 fManager = FairRootManager::Instance();
85 LOG(error) <<
"Could not initialise MCDataManager" << std::endl;
89 LOG(warn) <<
"Could not initialise FairRootManager" << std::endl;
97 std::string px =
"Px_" + decName;
98 std::string py =
"Py_" + decName;
99 std::string pt =
"Pt_" + decName;
100 std::string rx =
"x_" + decName;
101 std::string ry =
"y_" + decName;
102 std::string rt =
"t_" + decName;
103 std::string pullx = decName +
" Pull x";
104 std::string pully = decName +
" Pull y";
105 std::string pullt = decName +
" Pull t";
106 std::string resx = decName +
" Residual x";
107 std::string resy = decName +
" Residual y";
108 std::string rest = decName +
" Residual t";
110 TH1F* pullxHist =
new TH1F(px.c_str(),
115 TH1F* pullyHist =
new TH1F(py.c_str(),
120 TH1F* pulltHist =
new TH1F(pt.c_str(),
126 TH1F* residualxHist =
new TH1F(rx.c_str(),
131 TH1F* residualyHist =
new TH1F(ry.c_str(),
136 TH1F* residualtHist =
new TH1F(rt.c_str(),
149 pullxHist->GetXaxis()->SetTitle(
"Pull x");
150 pullyHist->GetXaxis()->SetTitle(
"Pull y");
151 pulltHist->GetXaxis()->SetTitle(
"Pull t");
152 residualxHist->GetXaxis()->SetTitle(
"Residual x");
153 residualyHist->GetXaxis()->SetTitle(
"Residual y");
154 residualtHist->GetXaxis()->SetTitle(
"Residual t");
159 hists[
i][3] = residualxHist;
160 hists[
i][4] = residualyHist;
161 hists[
i][5] = residualtHist;
164 LOG(info) <<
"CbmRecoQa Success Initiliasing Histograms for: " << decName;
174 static int event = 0;
175 LOG(info) <<
"CbmRecoQa for Event " <<
event++;
176 for (
unsigned int k = 0; k <
detectors.size(); k++) {
185 std::string filename =
outname +
".qa.hists.root";
186 pullresfile =
new TFile(filename.c_str(),
"Recreate");
188 for (
unsigned int i = 0;
i <
hists.size();
i++) {
191 for (
unsigned int k = 0; k <
hists[
i].size(); k++) {
193 LOG(info) <<
"CbmRecoQa Histogramm " <<
hists[
i][k]->GetName()
194 <<
" Entries " <<
hists[
i][k]->GetEntries();
195 if (
hists[
i][k]->GetEntries() > 0)
hists[
i][k]->Write();
197 gDirectory->cd(
"..");
204 if (
verbosity > 1) LOG(info) <<
"CbmRecoQa Record called for: " << decName;
208 TClonesArray* listHits =
nullptr;
209 TClonesArray* listHitMatches =
nullptr;
213 TClonesArray* listClusterMatches =
nullptr;
215 if (decName ==
"sts") {
216 listHits = (TClonesArray*) (
fManager->GetObject(
"StsHit"));
217 listHitMatches = (TClonesArray*) (
fManager->GetObject(
"StsHitMatch"));
220 (TClonesArray*) (
fManager->GetObject(
"StsClusterMatch"));
222 if (listPoints ==
nullptr) { LOG(warn) <<
"No StsPoint data!"; }
225 if (listHits && listHitMatches) {
227 TVector3 hitPos(.0, .0, .0);
228 TVector3 mcPos(.0, .0, .0);
229 TVector3 hitErr(.0, .0, .0);
231 int nEnt = listHits->GetEntries();
233 LOG(info) <<
"CbmRecoQa for " << decName <<
" found " << nEnt
235 for (
int j = 0; j < nEnt; j++) {
237 float bestWeight = 0.f;
238 float totalWeight = 0.f;
245 if (listClusterMatches) {
253 LOG(info) <<
"Front Match: " << front_match->
ToString()
254 <<
" Back Match: " << back_match->
ToString();
256 for (
int frontlink_c = 0; frontlink_c < front_match->
GetNofLinks();
261 for (
int backlink_c = 0; backlink_c < back_match->
GetNofLinks();
266 LOG(info) <<
"FrontLink: " << frontLink.
ToString()
267 <<
" BackLink: " << backLink.
ToString();
268 if (frontLink == backLink) {
277 for (
int iLink = 0; iLink < curr_match.
GetNofLinks(); iLink++) {
279 totalWeight += tmpweight;
281 if (tmpweight > bestWeight) {
282 bestWeight = tmpweight;
284 link = curr_match.
GetLink(iLink);
285 if (
verbosity > 3) LOG(info) <<
"Found Link for current HIT";
292 if (curr_point == 0)
continue;
294 curr_point->GetTime()
299 if (
verbosity > 3) LOG(info) <<
"Calculated Position Error";
301 mcPos.SetX(curr_point->
GetX(hitPos.Z()));
302 mcPos.SetY(curr_point->
GetY(hitPos.Z()));
303 mcPos.SetZ(hitPos.Z());
304 if (
verbosity > 3) LOG(info) <<
"Calculated MCPos";
308 hists[decNum][0]->Fill((hitPos.X() - mcPos.X())
309 / curr_hit->
GetDx());
311 hists[decNum][1]->Fill((hitPos.Y() - mcPos.Y())
312 / curr_hit->
GetDy());
316 hists[decNum][3]->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
317 hists[decNum][4]->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
320 LOG(warn) <<
"CBMRECOQA WARNING :-- No Sts Cluster Matches found!"
325 LOG(warn) <<
"CBMRECOQA WARNING :-- NO Data for Reco QA found! "
326 <<
"Detector: " << decName << std::endl;
328 }
else if (decName ==
"mvd") {
329 listHits = (TClonesArray*) (
fManager->GetObject(
"MvdHit"));
330 listHitMatches = (TClonesArray*) (
fManager->GetObject(
"MvdHitMatch"));
333 TClonesArray* listMvdPoints =
334 (TClonesArray*) (
fManager->GetObject(
"MvdPoint"));
336 if (listHits && listHitMatches) {
338 TVector3 hitPos(.0, .0, .0);
339 TVector3 mcPos(.0, .0, .0);
340 TVector3 hitErr(.0, .0, .0);
342 int nEnt = listHits->GetEntries();
344 LOG(info) <<
"CbmRecoQa for " << decName <<
" found " << nEnt
346 for (
int j = 0; j < nEnt; j++) {
358 if (iMC < 0)
continue;
360 curr_point = (
CbmMvdPoint*) (listMvdPoints->At(iMC));
362 if (curr_point == 0)
continue;
366 if (
verbosity > 3) LOG(info) <<
"Calculated Position Error";
368 mcPos.SetX((curr_point->GetX() + curr_point->
GetXOut()) / 2.);
369 mcPos.SetY((curr_point->GetY() + curr_point->
GetYOut()) / 2.);
370 mcPos.SetZ(hitPos.Z());
371 if (
verbosity > 3) LOG(info) <<
"Calculated MCPos";
375 hists[decNum][0]->Fill((hitPos.X() - mcPos.X()) / curr_hit->
GetDx());
377 hists[decNum][1]->Fill((hitPos.Y() - mcPos.Y()) / curr_hit->
GetDy());
379 hists[decNum][3]->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
380 hists[decNum][4]->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
383 LOG(warn) <<
"CBMRECOQA WARNING :-- NO Data for Reco QA found! "
384 <<
"Detector: " << decName << std::endl;
386 }
else if (decName ==
"much") {
388 listHits = (TClonesArray*) (
fManager->GetObject(
"MuchPixelHit"));
389 listHitMatches = (TClonesArray*) (
fManager->GetObject(
"MuchPixelHitMatch"));
391 if (listPoints ==
nullptr) { LOG(warn) <<
"No MuchPoint data!"; }
393 if (listHits && listHitMatches) {
395 TVector3 hitPos(.0, .0, .0);
396 TVector3 mcPos(.0, .0, .0);
397 TVector3 hitErr(.0, .0, .0);
399 int nEnt = listHits->GetEntries();
401 LOG(info) <<
"CbmRecoQa for " << decName <<
" found " << nEnt
403 for (
int j = 0; j < nEnt; j++) {
405 float bestWeight = 0.f;
406 float totalWeight = 0.f;
414 for (
int iLink = 0; iLink < curr_match->
GetNofLinks(); iLink++) {
416 totalWeight += tmpweight;
418 if (tmpweight > bestWeight) {
419 bestWeight = tmpweight;
421 link = curr_match->
GetLink(iLink);
422 if (
verbosity > 3) LOG(info) <<
"Found Link for current HIT";
426 if (iMC < 0)
continue;
431 curr_point->GetTime()
435 if (curr_point == 0)
continue;
439 if (
verbosity > 3) LOG(info) <<
"Calculated Position Error";
441 mcPos.SetX((curr_point->
GetXIn() + curr_point->
GetXOut()) / 2.);
442 mcPos.SetY((curr_point->
GetYIn() + curr_point->
GetYOut()) / 2.);
443 mcPos.SetZ(hitPos.Z());
444 if (
verbosity > 3) LOG(info) <<
"Calculated MCPos";
449 hists[decNum][0]->Fill((hitPos.X() - mcPos.X()) / curr_hit->
GetDx());
451 hists[decNum][1]->Fill((hitPos.Y() - mcPos.Y()) / curr_hit->
GetDy());
453 hists[decNum][2]->Fill(
456 hists[decNum][3]->Fill((hitPos.X() - mcPos.X()));
457 hists[decNum][4]->Fill((hitPos.Y() - mcPos.Y()));
461 LOG(warn) <<
"CBMRECOQA WARNING :-- NO Data for Reco QA found! "
462 <<
"Detector: " << decName << std::endl;
466 LOG(warn) <<
"CBMRECOQA WARNING :-- NO matching Detector found ! "