BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
BmnClusteringQa.cxx
Go to the documentation of this file.
1/**
2 * \file BmnClusteringQa.cxx
3 * \brief FairTask for clustering performance calculation.
4 * \author Andrey Lebedev <andrey.lebedev@gsi.de>
5 * \date 2011
6 */
7
8#include "BmnClusteringQa.h"
9
10#include "BmnClusteringQaReport.h"
11#include "BmnGemStripStation.h"
12#include "BmnGemStripStationSet_RunSpring2017.h"
13#include "BmnHistManager.h"
14#include "BmnHit.h"
15#include "BmnMatch.h"
16#include "BmnMath.h"
17#include "BmnSiliconStationSet.h"
18#include "CbmMCTrack.h"
19#include "FairMCEventHeader.h"
20#include "FairMCPoint.h"
21#include "FairRootManager.h"
22// #include "CbmStsAddress.h"
23
24#include "TClonesArray.h"
25#include "TF1.h"
26#include "TH1.h"
27#include "TH1D.h"
28#include "TH2.h"
29#include "TProfile.h"
30#include "TProfile2D.h"
31#include "TSystem.h"
32
33#include <boost/assign/list_of.hpp>
34#include <cmath>
35#include <fstream>
36#include <sstream>
37using boost::assign::list_of;
38using std::binary_search;
39using std::cout;
40using std::map;
41using std::vector;
42using namespace TMath;
43
44BmnClusteringQa::BmnClusteringQa()
45 : fHM(NULL)
46 , fOutputDir("./")
47 , fMCTracks(NULL)
48 , fGemPoints(NULL)
49 , fGemDigits(NULL)
50 ,
51 // fGemClusters(NULL),
52 fGemHits(NULL)
53 , fGemDigitMatches(NULL)
54 ,
55 // fGemClusterMatches(NULL),
56 fGemHitMatches(NULL)
57 , fGemDetector(NULL)
58 ,
59 // fTof1Points(NULL),
60 // fTof1Hits(NULL),
61 // fTof1HitMatches(NULL),
62 // fTof2Points(NULL),
63 // fTof2Hits(NULL),
64 // fTof2HitMatches(NULL),
65 fDch1Points(NULL)
66 , fDch1Hits(NULL)
67 , fDch1HitMatches(NULL)
68 ,
69 // fDch2Points(NULL),
70 // fDch2Hits(NULL),
71 // fDch2HitMatches(NULL)
72 fPrimes(kFALSE)
73{}
74
75BmnClusteringQa::~BmnClusteringQa()
76{
77 if (fHM)
78 delete fHM;
79}
80
81InitStatus BmnClusteringQa::Init()
82{
83 // Create histogram manager which is used throughout the program
84 fHM = new BmnHistManager();
85 TString gPathConfig = gSystem->Getenv("VMCWORKDIR");
86 TString gPathSiConfig = gPathConfig + "/silicon/XMLConfigs/";
87 TString confSi = "SiliconRunSpring2018.xml";
88 fSilDetector = new BmnSiliconStationSet(gPathSiConfig + confSi);
89 TString gPathGemConfig = gPathConfig + "/gem/XMLConfigs/";
90 TString confGem = "GemRunSpring2018.xml";
91 fGemDetector = new BmnGemStripStationSet(gPathGemConfig + confGem);
92
93 ReadDataBranches();
94
95 CreateHistograms();
96 return kSUCCESS;
97}
98
99void BmnClusteringQa::Exec(Option_t* opt)
100{
101 // FillEventCounterHistograms();
102 fHM->H1("hen_EventNo_ClusteringQa")->Fill(0.5);
103 ReadEventHeader();
104
105 ProcessPoints(fGemPoints, "Gem", kGEM);
106 // ProcessPoints(fTof1Points, "Tof1", kTOF1);
107 // ProcessPoints(fDch1Points, "Dch1", kDCH1);
108 // ProcessPoints(fDch2Points, "Dch2", kDCH2);
109 // ProcessPoints(fTof2Points, "Tof2", kTOF);
110
111 // ProcessDigis(fGemDigis, fGemDigiMatches, "Gem", kGEM);
112 // ProcessClusters(fGemClusters, fGemClusterMatches, "Gem", kGEM);
113
114 ProcessHits(fGemHits, fGemHitMatches, "Gem", kGEM);
115 // ProcessHits(fTof1Hits, fTof1HitMatches, "Tof1", kTOF1);
116 // ProcessHits(fDch1Hits, fDch1HitMatches, "Dch1", kDCH1);
117 // ProcessHits(fDch2Hits, fDch2HitMatches, "Dch2", kDCH2);
118 // ProcessHits(fTof2Hits, fTof2HitMatches, "Tof2", kTOF);
119
120 FillResidualAndPullHistograms(fGemPoints, fGemHits, fGemHitMatches, "Gem", kGEM);
121 // FillResidualAndPullHistograms(fTof1Points, fTof1Hits, fTof1HitMatches, "Tof1", kTOF1);
122 // FillResidualAndPullHistograms(fDch1Points, fDch1Hits, fDch1HitMatches, "Dch1", kDCH1);
123 // FillResidualAndPullHistograms(fDch2Points, fDch2Hits, fDch2HitMatches, "Dch2", kDCH2);
124 // FillResidualAndPullHistograms(fTof2Points, fTof2Hits, fTof2HitMatches, "Tof2", kTOF);
125
126 FillHitEfficiencyHistograms(fGemPoints, fGemHits, fGemHitMatches, "Gem", kGEM);
127 // FillHitEfficiencyHistograms(fTof1Points, fTof1Hits, fTof1HitMatches, "Tof1", kTOF1);
128 // FillHitEfficiencyHistograms(fDch1Points, fDch1Hits, fDch1HitMatches, "Dch1", kDCH1);
129 // FillHitEfficiencyHistograms(fDch2Points, fDch2Hits, fDch2HitMatches, "Dch2", kDCH2);
130 // FillHitEfficiencyHistograms(fTof2Points, fTof2Hits, fTof2HitMatches, "Tof2", kTOF);
131
132 // std::cout << "BmnClusteringQa::Exec: event=" << fHM->H1("hen_EventNo_ClusteringQa")->GetEntries() <<
133 // std::endl;
134}
135
136void BmnClusteringQa::Finish()
137{
138 fHM->WriteToFile();
139 // BmnSimulationReport* report = new BmnClusteringQaReport(fGemDetector->GetNStations(), 1);
140 BmnSimulationReport* report = new BmnClusteringQaReport(fGemDetector->GetNStations(), fSilDetector->GetNStations());
141 report->Create(fHM, fOutputDir);
142 delete report;
143}
144
145void BmnClusteringQa::ReadEventHeader()
146{
147 FairMCEventHeader* evHead = (FairMCEventHeader*)FairRootManager::Instance()->GetObject("MCEventHeader.");
148 fHM->H1("Impact parameter")->Fill(evHead->GetB());
149 // fHM->H1("Impact parameter")->Fill(evHead->GetRotZ());
150 fHM->H1("Multiplicity")->Fill(evHead->GetNPrim());
151 fHM->H2("Impact_Mult")->Fill(evHead->GetB(), evHead->GetNPrim());
152}
153
154void BmnClusteringQa::ReadDataBranches()
155{
156 FairRootManager* ioman = FairRootManager::Instance();
157 assert(ioman != NULL);
158
159 fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack");
160
161 fGemPoints = (TClonesArray*)ioman->GetObject("StsPoint");
162 fSilPoints = (TClonesArray*)ioman->GetObject("SiliconPoint");
163 // fTof1Points = (TClonesArray*) ioman->GetObject("TOF400Point");
164 // fTof2Points = (TClonesArray*) ioman->GetObject("TOF700Point");
165 fDch1Points = (TClonesArray*)ioman->GetObject("DCHPoint");
166 // fDch2Points = (TClonesArray*) ioman->GetObject("DCH2Point");
167
168 fGemHits = (TClonesArray*)ioman->GetObject("BmnGemStripHit");
169 fSilHits = (TClonesArray*)ioman->GetObject("BmnSiliconHit");
170 // fTof1Hits = (TClonesArray*) ioman->GetObject("TOF1Hit");
171 // fTof2Hits = (TClonesArray*) ioman->GetObject("BmnTof2Hit");
172 // fDch1Hits = (TClonesArray*) ioman->GetObject("BmnDchHit");
173 // fDch2Hits = (TClonesArray*) ioman->GetObject("BmnDch2Hit0");
174
175 fGemHitMatches = (TClonesArray*)ioman->GetObject("BmnGemStripHitMatch");
176 fSilHitMatches = (TClonesArray*)ioman->GetObject("BmnSiliconHitMatch");
177 // fTof1HitMatches = (TClonesArray*) ioman->GetObject("Tof1HitMatch");
178 // fTof2HitMatches = (TClonesArray*) ioman->GetObject("Tof2HitMatch");
179 // fDch1HitMatches = (TClonesArray*) ioman->GetObject("Dch1HitMatch");
180 // fDch2HitMatches = (TClonesArray*) ioman->GetObject("Dch2HitMatch");
181}
182
183Int_t BmnClusteringQa::GetStationId(Int_t address, DetectorId detId)
184{
185 // if (detId == kGEM) return CbmStsAddress::GetElementId(address, kGEMStation); //FIXME!!!!!!!!!!!!!
186 return 0;
187}
188
189void BmnClusteringQa::ProcessPoints(const TClonesArray* points, const string& detName, DetectorId detId)
190{
191 string histName = "_hno_NofObjects_" + detName + "Points_Station";
192 if (NULL == points || !fHM->Exists(histName))
193 return;
194 for (Int_t i = 0; i < points->GetEntriesFast(); i++) {
195 // const FairMCPoint* point = (const FairMCPoint*) (points->At(i));
196 // Int_t station = stationNumber("Gem", point->GetZ());
197 // fHM->H1(histName)->Fill(station); //FIXME! Get Station ID from point!
198 }
199}
200
201void BmnClusteringQa::ProcessDigis(const TClonesArray* digis,
202 const TClonesArray* digiMatches,
203 const string& detName,
204 DetectorId detId)
205{
206 // if (NULL == digis || !fHM->Exists("hno_NofObjects_" + detName + "Digis_Station")) return;
207 // for (Int_t i = 0; i < digis->GetEntriesFast(); i++) {
208 // const CbmDigi* digi = static_cast<const CbmDigi*> (digis->At(i));
209 // const BmnMatch* digiMatch = static_cast<const BmnMatch*> (digiMatches->At(i));
210 // Int_t stationId = GetStationId(digi->GetAddress(), detId);
211 // fHM->H1("hno_NofObjects_" + detName + "Digis_Station")->Fill(stationId);
212 // fHM->H1("hpa_" + detName + "Digi_NofPointsInDigi_H1")->Fill(digiMatch->GetNofLinks());
213 // fHM->H1("hpa_" + detName + "Digi_NofPointsInDigi_H2")->Fill(stationId, digiMatch->GetNofLinks());
214 // }
215}
216
217void BmnClusteringQa::ProcessClusters(const TClonesArray* clusters,
218 const TClonesArray* clusterMatches,
219 const string& detName,
220 DetectorId detId)
221{
222 // if (NULL != clusters && fHM->Exists("hno_NofObjects_" + detName + "Clusters_Station")) {
223 // for (Int_t i = 0; i < clusters->GetEntriesFast(); i++) {
224 // const CbmCluster* cluster = static_cast<const CbmCluster*> (clusters->At(i));
225 // const BmnMatch* clusterMatch = static_cast<const BmnMatch*> (clusterMatches->At(i));
226 // Int_t stationId = GetStationId(cluster->GetAddress(), detId);
227 // fHM->H1("hno_NofObjects_" + detName + "Clusters_Station")->Fill(stationId);
228 // fHM->H1("hpa_" + detName + "Cluster_NofDigisInCluster_H1")->Fill(cluster->GetNofDigis());
229 // fHM->H1("hpa_" + detName + "Cluster_NofDigisInCluster_H2")->Fill(stationId, cluster->GetNofDigis());
230 // fHM->H1("hpa_" + detName + "Cluster_NofPointsInCluster_H1")->Fill(clusterMatch->GetNofLinks());
231 // fHM->H1("hpa_" + detName + "Cluster_NofPointsInCluster_H2")->Fill(stationId,
232 // clusterMatch->GetNofLinks());
233 // }
234 // }
235}
236
237void BmnClusteringQa::ProcessHits(const TClonesArray* hits,
238 const TClonesArray* hitMatches,
239 const string& detName,
240 DetectorId detId)
241{
242 // if (NULL != hits && fHM->Exists("_hno_NofObjects_" + detName + "Hits_Station")) {
243 // for (Int_t i = 0; i < hits->GetEntriesFast(); i++) {
244 // const BmnMatch* hitMatch = (const BmnMatch*) (hitMatches->At(i));
245 // const BmnHit* hit = (const BmnHit*) (hits->At(i));
246 // Int_t stationId = hit->GetStation();
247 // fHM->H1("_hno_NofObjects_" + detName + "Hits_Station")->Fill(stationId);
248 // fHM->H1("_hpa_" + detName + "Hit_SigmaX_H1")->Fill(hit->GetDx());
249 // fHM->H1("_hpa_" + detName + "Hit_SigmaX_H2")->Fill(stationId, hit->GetDx());
250 // fHM->H1("_hpa_" + detName + "Hit_SigmaY_H1")->Fill(hit->GetDy());
251 // fHM->H1("_hpa_" + detName + "Hit_SigmaY_H2")->Fill(stationId, hit->GetDy());
252 // fHM->H1("_hpa_" + detName + "Hit_NofPointsInHit_H1")->Fill(hitMatch->GetNofLinks());
253 // fHM->H1("_hpa_" + detName + "Hit_NofPointsInHit_H2")->Fill(stationId, hitMatch->GetNofLinks());
254 // }
255 // }
256
257 for (Int_t iSt = 0; iSt < fGemDetector->GetNStations(); ++iSt) {
258 TString resXname = Form("ResX_%dst_gem", iSt);
259 TString resYname = Form("ResY_%dst_gem", iSt);
260 TString pntXhitXname = Form("PntX_vs_HitX_%dst_gem", iSt);
261 TString pntYhitYname = Form("PntY_vs_HitY_%dst_gem", iSt);
262 TString occupname = Form("Occupancy_%dst_gem", iSt);
263 TString pullXname = Form("PullX_%dst_gem", iSt);
264 TString pullYname = Form("PullY_%dst_gem", iSt);
265
266 /*for (Int_t i = 0; i < fGemPoints->GetEntriesFast(); ++i) {
267 const FairMCPoint* pnt = (const FairMCPoint*) (fGemPoints->At(i));
268 //const FairMCPoint* pnt = (const FairMCPoint*) (fGemPoints->At(hitMatch->GetMatchedLink().GetIndex()));
269// Int_t station = stationNumber("Gem", pnt->GetZ());
270// if (station != iSt) continue;
271 fHM->H2(occupname.Data())->Fill(pnt->GetX(), pnt->GetY());
272 }*/
273
274 for (Int_t i = 0; i < fGemHits->GetEntriesFast(); ++i) {
275 const BmnHit* hit = (const BmnHit*)(fGemHits->At(i));
276 if (hit->GetStation() != iSt)
277 continue;
278 const BmnMatch* hitMatch = (const BmnMatch*)(fGemHitMatches->At(i));
279 const FairMCPoint* pnt = (const FairMCPoint*)(fGemPoints->At(hitMatch->GetMatchedLink().GetIndex()));
280 const Float_t resX = pnt->GetX() - hit->GetX();
281 const Float_t resY = pnt->GetY() - hit->GetY();
282 fHM->H1(resXname.Data())->Fill(resX);
283 fHM->H1(resYname.Data())->Fill(resY);
284 fHM->H2(pntXhitXname.Data())->Fill(pnt->GetX(), hit->GetX());
285 fHM->H2(pntYhitYname.Data())->Fill(pnt->GetY(), hit->GetY());
286
287 fHM->H1(pullXname)->Fill(resX / hit->GetDx());
288 fHM->H1(pullYname)->Fill(resY / hit->GetDy());
289
290 fHM->H2(occupname.Data())->Fill(pnt->GetX(), pnt->GetY());
291 }
292 }
293
294 for (Int_t iSt = 0; iSt < fSilDetector->GetNStations(); ++iSt) {
295 // for (Int_t iSt = 0; iSt < 1; ++iSt){
296 TString occupname = Form("Occupancy_%dst_sil", iSt);
297 TString resXname = Form("ResX_%dst_sil", iSt);
298 TString resYname = Form("ResY_%dst_sil", iSt);
299 TString pntXhitXname = Form("PntX_vs_HitX_%dst_sil", iSt);
300 TString pntYhitYname = Form("PntY_vs_HitY_%dst_sil", iSt);
301 TString pullXname = Form("PullX_%dst_sil", iSt);
302 TString pullYname = Form("PullY_%dst_sil", iSt);
303
304 for (Int_t i = 0; i < fSilHits->GetEntriesFast(); ++i) {
305 const BmnHit* hit = (const BmnHit*)(fSilHits->At(i));
306 if (hit->GetStation() != iSt)
307 continue;
308 const BmnMatch* hitMatch = (const BmnMatch*)(fSilHitMatches->At(i));
309 const FairMCPoint* pnt = (const FairMCPoint*)(fSilPoints->At(hitMatch->GetMatchedLink().GetIndex()));
310
311 const Float_t resX = pnt->GetX() - hit->GetX();
312 const Float_t resY = pnt->GetY() - hit->GetY();
313
314 fHM->H1(resXname.Data())->Fill(resX);
315 fHM->H1(resYname.Data())->Fill(resY);
316 fHM->H2(pntXhitXname.Data())->Fill(pnt->GetX(), hit->GetX());
317 fHM->H2(pntYhitYname.Data())->Fill(pnt->GetY(), hit->GetY());
318
319 fHM->H1(pullXname)->Fill(resX / hit->GetDx());
320 fHM->H1(pullYname)->Fill(resY / hit->GetDy());
321
322 fHM->H2(occupname.Data())->Fill(pnt->GetX(), pnt->GetY());
323 }
324 }
325
326 // for (Int_t iSt = 0; iSt < 3; ++iSt){
327 // TString occupname = Form("Occupancy_%dst_dch1", iSt);
328 //
329 // for (Int_t i = 0; i < fDch1Hits->GetEntriesFast(); ++i) {
330 // printf("%d Number of Dch points \n",fDch1Hits->GetEntriesFast());
331 // const BmnHit* hit = (const BmnHit*) (fDch1Hits->At(i));
332 // if (hit->GetStation() != iSt) continue;
333 // const BmnMatch* hitMatch = (const BmnMatch*) (fDch1HitMatches->At(i));
334 // const FairMCPoint* pnt = (const FairMCPoint*)
335 // (fDch1Points->At(hitMatch->GetMatchedLink().GetIndex()));
336 //
337 // fHM->H2(occupname.Data())->Fill(pnt->GetX(), pnt->GetY());
338 // }
339 // }
340}
341
342void BmnClusteringQa::FillEventCounterHistograms()
343{
344 // if (NULL != fMvdPoints && fHM->Exists("hno_NofObjects_MvdPoints_Event"))
345 // fHM->H1("hno_NofObjects_MvdPoints_Event")->Fill(fMvdPoints->GetEntriesFast()); if (NULL != fMvdDigis &&
346 // fHM->Exists("hno_NofObjects_MvdDigis_Event"))
347 // fHM->H1("hno_NofObjects_MvdDigis_Event")->Fill(fMvdDigis->GetEntriesFast()); if (NULL != fMvdClusters &&
348 // fHM->Exists("hno_NofObjects_MvdClusters_Event"))
349 // fHM->H1("hno_NofObjects_MvdClusters_Event")->Fill(fMvdClusters->GetEntriesFast()); if (NULL != fMvdHits &&
350 // fHM->Exists("hno_NofObjects_MvdHits_Event"))
351 // fHM->H1("hno_NofObjects_MvdHits_Event")->Fill(fMvdHits->GetEntriesFast());
352
353 // if (NULL != fGemPoints && fHM->Exists("_hno_NofObjects_GemPoints_Event"))
354 // fHM->H1("_hno_NofObjects_GemPoints_Event")->Fill(fGemPoints->GetEntriesFast()); if (NULL != fGemDigits &&
355 // fHM->Exists("hno_NofObjects_StsDigis_Event"))
356 // fHM->H1("hno_NofObjects_StsDigis_Event")->Fill(fGemDigits->GetEntriesFast()); if (NULL != fGemClusters &&
357 // fHM->Exists("hno_NofObjects_StsClusters_Event"))
358 // fHM->H1("hno_NofObjects_StsClusters_Event")->Fill(fGemClusters->GetEntriesFast());
359 // if (NULL != fGemHits && fHM->Exists("_hno_NofObjects_GemHits_Event"))
360 // fHM->H1("_hno_NofObjects_GemHits_Event")->Fill(fGemHits->GetEntriesFast());
361
362 // if (NULL != fTof1Points && fHM->Exists("_hno_NofObjects_Tof1Points_Event"))
363 // fHM->H1("_hno_NofObjects_Tof1Points_Event")->Fill(fTof1Points->GetEntriesFast()); if (NULL != fTof2Points &&
364 // fHM->Exists("_hno_NofObjects_Tof2Points_Event"))
365 // fHM->H1("_hno_NofObjects_Tof2Points_Event")->Fill(fTof2Points->GetEntriesFast());
366 // if (NULL != fDch1Points && fHM->Exists("_hno_NofObjects_Dch1Points_Event"))
367 // fHM->H1("_hno_NofObjects_Dch1Points_Event")->Fill(fDch1Points->GetEntriesFast()); if (NULL != fDch2Points &&
368 // fHM->Exists("_hno_NofObjects_Dch2Points_Event"))
369 // fHM->H1("_hno_NofObjects_Dch2Points_Event")->Fill(fDch2Points->GetEntriesFast());
370 //
371 // if (NULL != fTof1Hits && fHM->Exists("_hno_NofObjects_Tof1Hits_Event"))
372 // fHM->H1("_hno_NofObjects_Tof1Hits_Event")->Fill(fTof1Hits->GetEntriesFast()); if (NULL != fTof2Hits &&
373 // fHM->Exists("_hno_NofObjects_Tof2Hits_Event"))
374 // fHM->H1("_hno_NofObjects_Tof2Hits_Event")->Fill(fTof2Hits->GetEntriesFast()); if (NULL != fDch1Hits &&
375 // fHM->Exists("_hno_NofObjects_Dch1Hits_Event"))
376 // fHM->H1("_hno_NofObjects_Dch1Hits_Event")->Fill(fDch1Hits->GetEntriesFast()); if (NULL != fDch2Hits &&
377 // fHM->Exists("_hno_NofObjects_Dch2Hits_Event"))
378 // fHM->H1("_hno_NofObjects_Dch2Hits_Event")->Fill(fDch2Hits->GetEntriesFast());
379}
380
381void BmnClusteringQa::FillResidualAndPullHistograms(const TClonesArray* points,
382 const TClonesArray* hits,
383 const TClonesArray* hitMatches,
384 const string& detName,
385 DetectorId detId)
386{
387 if (NULL == points || NULL == hits || NULL == hitMatches)
388 return;
389 string nameResidualX = "_hrp_" + detName + "_ResidualX_H2";
390 string nameResidualY = "_hrp_" + detName + "_ResidualY_H2";
391 string namePullX = "_hrp_" + detName + "_PullX_H2";
392 string namePullY = "_hrp_" + detName + "_PullY_H2";
393 if (!fHM->Exists(nameResidualX) || !fHM->Exists(nameResidualY) || !fHM->Exists(namePullX)
394 || !fHM->Exists(namePullY))
395 return;
396
397 for (Int_t iHit = 0; iHit < hits->GetEntriesFast(); iHit++) {
398 const BmnMatch* match = (const BmnMatch*)(hitMatches->At(iHit));
399 const FairMCPoint* point = (const FairMCPoint*)(points->At(match->GetMatchedLink().GetIndex()));
400 if (point == NULL)
401 continue;
402 const BmnHit* hit = (const BmnHit*)(hits->At(iHit));
403 Float_t residualX = point->GetX() - hit->GetX();
404 Float_t residualY = point->GetY() - hit->GetY();
405 Int_t stationId = hit->GetStation();
406 fHM->H2(nameResidualX)->Fill(stationId, residualX);
407 fHM->H2(nameResidualY)->Fill(stationId, residualY);
408 fHM->H2(namePullX)->Fill(stationId, residualX / hit->GetDx());
409 fHM->H2(namePullY)->Fill(stationId, residualY / hit->GetDy());
410 }
411}
412
413void BmnClusteringQa::FillHitEfficiencyHistograms(const TClonesArray* points,
414 const TClonesArray* hits,
415 const TClonesArray* hitMatches,
416 const string& detName,
417 DetectorId detId)
418{
419 if (NULL == points || NULL == hits || NULL == hitMatches)
420 return;
421 string accName = "_hhe_" + detName + "_All_Acc_Station";
422 if (NULL == points || !fHM->Exists(accName))
423 return;
424 for (Int_t i = 0; i < points->GetEntriesFast(); i++) {
425 const FairMCPoint* point = (const FairMCPoint*)(points->At(i));
426 fHM->H1(accName)->Fill(GetStationId(point->GetDetectorID(), detId));
427 }
428
429 string recName = "_hhe_" + detName + "_All_Rec_Station";
430 string cloneName = "_hhe_" + detName + "_All_Clone_Station";
431 set<Int_t> mcPointSet; // IDs of MC points
432 Int_t nofHits = hits->GetEntriesFast();
433 for (Int_t iHit = 0; iHit < nofHits; iHit++) {
434 const BmnMatch* match = (const BmnMatch*)(hitMatches->At(iHit));
435 const BmnHit* hit = (const BmnHit*)(hits->At(iHit));
436 if (mcPointSet.find(match->GetMatchedLink().GetIndex()) == mcPointSet.end()) {
437 fHM->H1(recName)->Fill(hit->GetStation());
438 mcPointSet.insert(match->GetMatchedLink().GetIndex());
439 } else {
440 fHM->H1(cloneName)->Fill(hit->GetStation());
441 }
442 }
443}
444
445void BmnClusteringQa::CreateHistograms()
446{
447
448 Int_t nStationsGEM = fGemDetector->GetNStations();
449 Int_t nStationsSil = fSilDetector->GetNStations();
450 Int_t nStationsDCH1 = 3;
451
452 for (Int_t iSt = 0; iSt < nStationsGEM; ++iSt) {
453 TString occupname = Form("Occupancy_%dst_gem", iSt);
454 TString resXname = Form("ResX_%dst_gem", iSt);
455 TString resYname = Form("ResY_%dst_gem", iSt);
456 TString pntXhitXname = Form("PntX_vs_HitX_%dst_gem", iSt);
457 TString pntYhitYname = Form("PntY_vs_HitY_%dst_gem", iSt);
458 TString pullXname = Form("PullX_%dst_gem", iSt);
459 TString pullYname = Form("PullY_%dst_gem", iSt);
460
461 CreateH1(resXname.Data(), "ResX, cm", "Counter", 200, -5, 5);
462 CreateH1(resYname.Data(), "ResY, cm", "Counter", 200, -5, 5);
463 CreateH2(pntXhitXname.Data(), "X_{MC}, cm", "X_{reco}, cm", "", 400, -70, 70, 400, -70, 70);
464 CreateH2(pntYhitYname.Data(), "Y_{MC}, cm", "Y_{reco}, cm", "", 400, -70, 70, 400, -70, 70);
465
466 CreateH2(occupname.Data(), "X_{MC}, cm", "Y_{MC}, cm", "Occupancy, #frac{part}{event * cm^{2}}", 50, -100, 100,
467 50, -30, 30);
468
469 CreateH1(pullXname.Data(), "PullX", "Counter", 200, -10, 10);
470 CreateH1(pullYname.Data(), "PullY", "Counter", 200, -10, 10);
471 }
472
473 for (Int_t iSt = 0; iSt < nStationsSil; ++iSt) {
474 TString occupname = Form("Occupancy_%dst_sil", iSt);
475 TString resXname = Form("ResX_%dst_sil", iSt);
476 TString resYname = Form("ResY_%dst_sil", iSt);
477 TString pntXhitXname = Form("PntX_vs_HitX_%dst_sil", iSt);
478 TString pntYhitYname = Form("PntY_vs_HitY_%dst_sil", iSt);
479 TString pullXname = Form("PullX_%dst_sil", iSt);
480 TString pullYname = Form("PullY_%dst_sil", iSt);
481
482 CreateH1(resXname.Data(), "ResX, cm", "Counter", 200, -1, 1);
483 CreateH1(resYname.Data(), "ResY, cm", "Counter", 200, -1, 1);
484 CreateH2(pntXhitXname.Data(), "X_{MC}, cm", "X_{reco}, cm", "", 400, -20, 20, 400, -20, 20);
485 CreateH2(pntYhitYname.Data(), "Y_{MC}, cm", "Y_{reco}, cm", "", 400, -20, 20, 400, -20, 20);
486
487 CreateH2(occupname.Data(), "X_{MC}, cm", "Y_{MC}, cm", "Occupancy, #frac{part}{event * cm^{2}}", 50, -20, 20,
488 50, -20, 20);
489
490 CreateH1(pullXname.Data(), "PullX", "Counter", 200, -10, 10);
491 CreateH1(pullYname.Data(), "PullY", "Counter", 200, -10, 10);
492 }
493
494 for (Int_t iSt = 0; iSt < nStationsDCH1; ++iSt) {
495
496 TString occupname = Form("Occupancy_%dst_dch1", iSt);
497
498 CreateH2(occupname.Data(), "X_{MC}, cm", "Y_{MC}, cm", "Occupancy, #frac{part}{event * cm^{2}}", 50, -20, 20,
499 50, -20, 20);
500 }
501
502 // Histogram stores number of events
503 CreateH1("hen_EventNo_ClusteringQa", "", "", 1, 0, 1.);
504 CreateH1("Impact parameter", "b, fm", "Counter", 50, 0.0, 0.0);
505 CreateH1("Multiplicity", "N_{prim}", "Counter", 60, 0.0, 0.0);
506 CreateH2("Impact_Mult", "b, fm", "N_{prim}", "", 400, 0.0, 0.0, 400, 0.0, 0.0);
507}
508
509void BmnClusteringQa::CreateH1(const string& name,
510 const string& xTitle,
511 const string& yTitle,
512 Int_t nofBins,
513 Double_t minBin,
514 Double_t maxBin)
515{
516 TH1F* h = new TH1F(name.c_str(), string(name + ";" + xTitle + ";" + yTitle).c_str(), nofBins, minBin, maxBin);
517 fHM->Add(name, h);
518}
519
520void BmnClusteringQa::CreateH2(const string& name,
521 const string& xTitle,
522 const string& yTitle,
523 const string& zTitle,
524 Int_t nofBinsX,
525 Double_t minBinX,
526 Double_t maxBinX,
527 Int_t nofBinsY,
528 Double_t minBinY,
529 Double_t maxBinY)
530{
531 TH2F* h = new TH2F(name.c_str(), (name + ";" + xTitle + ";" + yTitle + ";" + zTitle).c_str(), nofBinsX, minBinX,
532 maxBinX, nofBinsY, minBinY, maxBinY);
533 fHM->Add(name, h);
534}
535
536/*void BmnClusteringQa::CreateNofObjectsHistograms(DetectorId detId, const string& detName) {
537 if (!fDet.GetDet(detId)) return;
538 assert(detId == kTOF || detId == kGEM || detId == kTOF1 || detId == kDCH);
539 Int_t nofBins = 5000;
540 Double_t minX = 0;
541 Double_t maxX = 5000;
542 string name = "_hno_NofObjects_" + detName;
543 fHM->Create1<TH1F > (name + "Points_Event", name + "Points_Event;Points per event;Counter", nofBins, minX, maxX);
544 // fHM->Create1<TH1F > (name + "Digis_Event", name + "Digis_Event;Digis per event;Counter", nofBins, minX,
545maxX);
546 // fHM->Create1<TH1F > (name + "Clusters_Event", name + "Clusters_Event;Clusters per event;Counter", nofBins,
547minX, maxX); fHM->Create1<TH1F > (name + "Hits_Event", name + "Hits_Event;Hits per event;Counter", nofBins, minX, maxX);
548}*/
549
550/*void BmnClusteringQa::CreateNofObjectsHistograms(DetectorId detId, const string& detName, const string& parameter,
551const string& xTitle) { if (!fDet.GetDet(detId)) return; Int_t nofBins = 15; Double_t minX = -0.5; Double_t maxX = 14.5;
552 string name = "_hno_NofObjects_" + detName;
553 fHM->Create1<TH1F > (name + "Points_" + parameter, name + "Points_" + parameter + ";" + xTitle + ";Points per
554event", nofBins, minX, maxX); fHM->Create1<TH1F > (name + "Digis_" + parameter, name + "Digis_" + parameter + ";" +
555xTitle + ";Digis per event", nofBins, minX, maxX);
556 // fHM->Create1<TH1F > (name + "Clusters_" + parameter, name + "Clusters_" + parameter + ";" + xTitle +
557";Clusters per event", nofBins, minX, maxX); fHM->Create1<TH1F > (name + "Hits_" + parameter, name + "Hits_" + parameter
558+ ";" + xTitle + ";Hits per event", nofBins, minX, maxX);
559}*/
560
561/*void BmnClusteringQa::CreateClusterParametersHistograms(DetectorId detId, const string& detName) {
562 if (!fDet.GetDet(detId)) return;
563 Int_t nofBinsStation = 17;
564 Double_t minStation = 0;
565 Double_t maxStation = 17;
566 Int_t nofBins = 100;
567 Double_t min = -0.5;
568 Double_t max = 99.5;
569 Int_t nofBinsSigma = 100;
570 Double_t minSigma = 0.;
571 Double_t maxSigma = 10;
572 Int_t nofBinsResidual = 100;
573 Double_t minResidual = -5.0;
574 Double_t maxResidual = 5.0;
575 Int_t nofBinsPull = 100;
576 Double_t minPull = -2.0;
577 Double_t maxPull = 2.0;
578
579 string nameH1 = "_hpa_" + detName + "Cluster_NofDigisInCluster_H1";
580 // fHM->Create1<TH1F > (nameH1, nameH1 + ";Number of digis;Yield", nofBins, min, max);
581 string nameH2 = "_hpa_" + detName + "Cluster_NofDigisInCluster_H2";
582 // fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Number of digis;Yield", nofBinsStation, minStation, max,
583nofBins, min, max);
584 // nameH1 = "hpa_" + detName + "Cluster_NofPointsInCluster_H1";
585 // fHM->Create1<TH1F > (nameH1, nameH1 + ";Number of points;Yield", nofBins, min, max);
586 // nameH2 = "hpa_" + detName + "Cluster_NofPointsInCluster_H2";
587 // fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Number of points;Yield", nofBinsStation, minStation, max,
588nofBins, min, max);
589 // nameH1 = "hpa_" + detName + "Digi_NofPointsInDigi_H1";
590 // fHM->Create1<TH1F > (nameH1, nameH1 + ";Number of points;Yield", nofBins, min, max);
591 // nameH2 = "hpa_" + detName + "Digi_NofPointsInDigi_H2";
592 // fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Number of points;Yield", nofBinsStation, minStation,
593maxStation, nofBins, min, max); nameH1 = "_hpa_" + detName + "Hit_NofPointsInHit_H1"; fHM->Create1<TH1F > (nameH1,
594nameH1 + ";Number of points;Yield", nofBins, min, max); nameH2 = "_hpa_" + detName + "Hit_NofPointsInHit_H2";
595 fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Number of points;Yield", nofBinsStation, minStation, max, nofBins,
596min, max); nameH1 = "_hpa_" + detName + "Hit_SigmaX_H1"; fHM->Create1<TH1F > (nameH1, nameH1 + ";#sigma_{X} [cm];Yield",
597nofBinsSigma, minSigma, maxSigma); nameH2 = "_hpa_" + detName + "Hit_SigmaX_H2"; fHM->Create2<TH2F > (nameH2, nameH2 +
598";Station;#sigma_{X} [cm];Yield", nofBinsStation, minStation, maxStation, nofBinsSigma, minSigma, maxSigma); nameH1 =
599"_hpa_" + detName + "Hit_SigmaY_H1"; fHM->Create1<TH1F > (nameH1, nameH1 + ";#sigma_{Y} [cm];Yield", nofBinsSigma,
600minSigma, maxSigma); nameH2 = "_hpa_" + detName + "Hit_SigmaY_H2"; fHM->Create2<TH2F > (nameH2, nameH2 +
601";Station;#sigma_{Y} [cm];Yield", nofBinsStation, minStation, maxStation, nofBinsSigma, minSigma, maxSigma);
602
603 // Residual and pull histograms
604 nameH2 = "_hrp_" + detName + "_ResidualX_H2";
605 fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Residual X [cm];Yield", nofBinsStation, minStation, maxStation,
606nofBinsResidual, minResidual, maxResidual); nameH2 = "_hrp_" + detName + "_ResidualY_H2"; fHM->Create2<TH2F > (nameH2,
607nameH2 + ";Station;Residual Y [cm];Yield", nofBinsStation, minStation, maxStation, nofBinsResidual, minResidual,
608maxResidual); nameH2 = "_hrp_" + detName + "_PullX_H2"; fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Pull X;Yield",
609nofBinsStation, minStation, maxStation, nofBinsPull, minPull, maxPull); nameH2 = "_hrp_" + detName + "_PullY_H2";
610 fHM->Create2<TH2F > (nameH2, nameH2 + ";Station;Pull Y;Yield", nofBinsStation, minStation, maxStation, nofBinsPull,
611minPull, maxPull);
612}*/
613
614/*void BmnClusteringQa::CreateHitEfficiencyHistograms(DetectorId detId, const string& detName, const string& parameter,
615const string& xTitle, Int_t nofBins, Double_t minBin, Double_t maxBin) { if (!fDet.GetDet(detId)) return; vector<string>
616types = list_of("Acc")("Rec")("Eff")("Clone")("CloneProb"); vector<string> cat = list_of("All"); for (Int_t iCat = 0;
617iCat < cat.size(); iCat++) { for (Int_t iType = 0; iType < types.size(); iType++) { string yTitle = (types[iType] ==
618"Eff") ? "Efficiency [%]" : (types[iType] == "CloneProb") ? "Probability [%]" : "Counter"; string histName = "_hhe_" +
619detName + "_" + cat[iCat] + "_" + types[iType] + "_" + parameter; string histTitle = histName + ";" + xTitle + ";" +
620yTitle; fHM->Add(histName, new TH1F(histName.c_str(), histTitle.c_str(), nofBins, minBin, maxBin));
621 }
622 }
623}*/
i
int i
Definition P4_F32vec4.h:22
BmnClusteringQaReport.h
Simulation report for clustering QA.
BmnClusteringQa.h
FairTask for clustering performance calculation.
Int_t
int Int_t
Float_t
float Float_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
Option_t
const char Option_t
Form
char * Form(const char *fmt,...)
gSystem
R__EXTERN TSystem * gSystem
kGEM
@ kGEM
Definition BmnDetectorList.h:15
BmnClusteringQaReport
Simulation report for clustering QA.
Definition BmnClusteringQaReport.h:25
BmnClusteringQa::fGemDetector
BmnGemStripStationSet * fGemDetector
Definition BmnClusteringQa.h:166
BmnClusteringQa::ReadDataBranches
void ReadDataBranches()
Read data branches.
Definition BmnClusteringQa.cxx:154
BmnClusteringQa::fDch1Points
TClonesArray * fDch1Points
Definition BmnClusteringQa.h:173
BmnClusteringQa::ProcessClusters
void ProcessClusters(const TClonesArray *clusters, const TClonesArray *clusterMatches, const string &detName, DetectorId detId)
Definition BmnClusteringQa.cxx:217
BmnClusteringQa::fGemHitMatches
TClonesArray * fGemHitMatches
Definition BmnClusteringQa.h:164
BmnClusteringQa::fDch1Hits
TClonesArray * fDch1Hits
Definition BmnClusteringQa.h:174
BmnClusteringQa::fGemPoints
TClonesArray * fGemPoints
Definition BmnClusteringQa.h:158
BmnClusteringQa::FillEventCounterHistograms
void FillEventCounterHistograms()
Definition BmnClusteringQa.cxx:342
BmnClusteringQa::fSilHitMatches
TClonesArray * fSilHitMatches
Definition BmnClusteringQa.h:171
BmnClusteringQa::CreateH2
void CreateH2(const string &name, const string &xTitle, const string &yTitle, const string &zTitle, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY)
Definition BmnClusteringQa.cxx:520
BmnClusteringQa::FillResidualAndPullHistograms
void FillResidualAndPullHistograms(const TClonesArray *points, const TClonesArray *hits, const TClonesArray *hitMatches, const string &detName, DetectorId detId)
Definition BmnClusteringQa.cxx:381
BmnClusteringQa::fSilDetector
BmnSiliconStationSet * fSilDetector
Definition BmnClusteringQa.h:167
BmnClusteringQa::ProcessDigis
void ProcessDigis(const TClonesArray *digis, const TClonesArray *digiMatches, const string &detName, DetectorId detId)
Definition BmnClusteringQa.cxx:201
BmnClusteringQa::ProcessHits
void ProcessHits(const TClonesArray *hits, const TClonesArray *hitMatches, const string &detName, DetectorId detId)
Definition BmnClusteringQa.cxx:237
BmnClusteringQa::BmnClusteringQa
BmnClusteringQa()
Constructor.
Definition BmnClusteringQa.cxx:44
BmnClusteringQa::~BmnClusteringQa
virtual ~BmnClusteringQa()
Destructor.
Definition BmnClusteringQa.cxx:75
BmnClusteringQa::fMCTracks
TClonesArray * fMCTracks
Definition BmnClusteringQa.h:151
BmnClusteringQa::fSilHits
TClonesArray * fSilHits
Definition BmnClusteringQa.h:169
BmnClusteringQa::Exec
virtual void Exec(Option_t *opt)
Derived from FairTask.
Definition BmnClusteringQa.cxx:99
BmnClusteringQa::ProcessPoints
void ProcessPoints(const TClonesArray *points, const string &detName, DetectorId detId)
Definition BmnClusteringQa.cxx:189
BmnClusteringQa::fDch1HitMatches
TClonesArray * fDch1HitMatches
Definition BmnClusteringQa.h:175
BmnClusteringQa::Finish
virtual void Finish()
Derived from FairTask.
Definition BmnClusteringQa.cxx:136
BmnClusteringQa::FillHitEfficiencyHistograms
void FillHitEfficiencyHistograms(const TClonesArray *points, const TClonesArray *hits, const TClonesArray *hitMatches, const string &detName, DetectorId detId)
Definition BmnClusteringQa.cxx:413
BmnClusteringQa::fGemDigits
TClonesArray * fGemDigits
Definition BmnClusteringQa.h:159
BmnClusteringQa::GetStationId
Int_t GetStationId(Int_t address, DetectorId detId)
Definition BmnClusteringQa.cxx:183
BmnClusteringQa::fSilPoints
TClonesArray * fSilPoints
Definition BmnClusteringQa.h:170
BmnClusteringQa::CreateH1
void CreateH1(const string &name, const string &xTitle, const string &yTitle, Int_t nofBins, Double_t minBin, Double_t maxBin)
Definition BmnClusteringQa.cxx:509
BmnClusteringQa::fGemHits
TClonesArray * fGemHits
Definition BmnClusteringQa.h:161
BmnClusteringQa::fHM
BmnHistManager * fHM
Definition BmnClusteringQa.h:147
BmnClusteringQa::Init
virtual InitStatus Init()
Derived from FairTask.
Definition BmnClusteringQa.cxx:81
BmnClusteringQa::fGemDigitMatches
TClonesArray * fGemDigitMatches
Definition BmnClusteringQa.h:162
BmnHistManager
Histogram manager.
Definition BmnHistManager.h:42
BmnHit::GetStation
Short_t GetStation() const
Definition BmnHit.h:35
BmnMatch::GetMatchedLink
const BmnLink & GetMatchedLink() const
Definition BmnMatch.h:40
BmnSimulationReport
Base class for simulation reports.
Definition BmnSimulationReport.h:34
BmnSimulationReport::Create
void Create(BmnHistManager *histManager, const string &outputDir)
Main function which creates report data.
Definition BmnSimulationReport.cxx:35
TObjArray::GetEntriesFast
Int_t GetEntriesFast() const
TObjArray::At
TObject * At(Int_t idx) const override
TString::Data
const char * Data() const