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