BmnRoot
Loading...
Searching...
No Matches
CbmStsDigitizeQa.cxx
Go to the documentation of this file.
1#include "CbmStsDigitizeQa.h"
2#include "CbmStsDigitize.h"
3#include "CbmHistManager.h"
4#include "CbmMatch.h"
5#include "CbmStsDigi.h"
6#include "CbmStsAddress.h"
7#include "CbmStsModule.h"
8#include "CbmStsElement.h"
9#include "CbmStsSetup.h"
10#include "CbmMCDataManager.h"
11#include "CbmSimulationReport.h"
12#include "CbmStsDigitizeQaReport.h"
13#include "CbmStsDigitizeParameters.h"
14
15#include "FairRootManager.h"
16#include "FairLogger.h"
17#include "FairRun.h"
18#include "FairRuntimeDb.h"
19
20#include "TGeoPhysicalNode.h"
21#include "TGeoMatrix.h"
22
23#include "TClonesArray.h"
24#include "TH1.h"
25#include "TF1.h"
26#include "TH1D.h"
27#include "TH2.h"
28#include "TProfile.h"
29#include "TProfile2D.h"
30
31
32using std::cout;
33using std::endl;
34using std::vector;
35using std::set;
36using std::map;
37using std::pair;
38using std::string;
39
40CbmStsDigitizeQa::CbmStsDigitizeQa(CbmStsDigitize * /*digitizer*/):
41 FairTask("CbmStsDigitizeQa")
42 , fDigiPar(nullptr)
43 , fHM(NULL)
44 , fOutputDir(" ")
45 , fStsDigis(NULL)
46 , fStsPoints(NULL)
47 , fSetup(NULL)
48 , fNofStation(8)
49 , fMaxScale(0)
50 , fOutFile(NULL)
51 , fnOfDigisChip()
52{
53}
54
55CbmStsDigitizeQa::~CbmStsDigitizeQa(){
56 if ( fHM ) delete fHM;
57}
58
59void CbmStsDigitizeQa::SetParContainers()
60{
61 fDigiPar = static_cast<CbmStsDigitizeParameters*>(FairRun::Instance()->GetRuntimeDb()->getContainer("CbmStsDigitizeParameters"));
62}
63
64InitStatus CbmStsDigitizeQa::Init(){
65 fSetup = CbmStsSetup::Instance();
66 fNofStation = fSetup -> GetNofStations();
67 fHM = new CbmHistManager();
68 fnOfDigisChip.resize(fNofStation);
69 for (Int_t iStation = 0; iStation < fNofStation; iStation ++){
70 CbmStsElement * stat = fSetup -> GetDaughter(iStation);
71 fnOfDigisChip[iStation].resize(stat -> GetNofDaughters());
72 for (Int_t iLad = 0; iLad < stat -> GetNofDaughters(); iLad++) {
73 CbmStsElement* ladd = stat -> GetDaughter(iLad);
74 fnOfDigisChip[iStation][iLad].resize(ladd -> GetNofDaughters());
75 for (Int_t iHla = 0; iHla < ladd -> GetNofDaughters(); iHla++) {
76 CbmStsElement* hlad = ladd -> GetDaughter(iHla);
77 fnOfDigisChip[iStation][iLad][iHla].resize(hlad -> GetNofDaughters());
78 for (Int_t iMod = 0; iMod < hlad -> GetNofDaughters(); iMod++) {
79 CbmStsModule* modu = static_cast<CbmStsModule*>(hlad -> GetDaughter(iMod));
80 Int_t nOfChips = Int_t(modu -> GetNofChannels() / 128.);
81 fnOfDigisChip[iStation][iLad][iHla][iMod].resize(nOfChips);
82 for (Int_t iChip = 0; iChip < nOfChips; iChip++){
83 fnOfDigisChip[iStation][iLad][iHla][iMod][iChip] = 0;
84 }
85 }
86 }
87 }
88 }
89
90 ReadDataBranches();
91 CreateHistograms();
92 return kSUCCESS;
93}
94
95void CbmStsDigitizeQa::Exec(Option_t* /*opt*/){
96 ProcessDigisAndPoints(fStsDigis, fStsPoints);
97 fHM -> H1("h_EventNo_DigitizeQa") -> Fill(0.5);
98}
99
100void CbmStsDigitizeQa::Finish(){
101 ProcessAngles();
102 Int_t nofEvents = fHM -> H1("h_EventNo_DigitizeQa") -> GetEntries();
103 TString fileName = fOutputDir + "/digiRateChip";
104 fileName += nofEvents;
105 fileName += ".dat";
106 TString rmFile = "rm " + fileName;
107 gSystem -> Exec(rmFile);
108 fOutFile.open(Form("%s",fileName.Data()), std::ofstream::app);
109 for (Int_t iStation = 0; iStation < fNofStation; iStation ++){
110 CbmStsElement * stat = fSetup -> GetDaughter(iStation);
111 for (Int_t iLad = 0; iLad < stat -> GetNofDaughters(); iLad++) {
112 CbmStsElement* ladd = stat -> GetDaughter(iLad);
113 for (Int_t iHla = 0; iHla < ladd -> GetNofDaughters(); iHla++) {
114 CbmStsElement* hlad = ladd -> GetDaughter(iHla);
115 for (Int_t iMod = 0; iMod < hlad -> GetNofDaughters(); iMod++) {
116 CbmStsModule* modu = static_cast<CbmStsModule*>(hlad -> GetDaughter(iMod));
117 if(modu -> GetNofChannels() != 2048) cout << "nofChannels = " << modu -> GetNofChannels() << endl;
118 Int_t nOfChips = Int_t(modu -> GetNofChannels() / 128.);
119 for (Int_t iChip = 0; iChip < nOfChips; iChip++)
120 fOutFile << iStation << "\t" << iLad << "\t" << iHla << "\t" << iMod << "\t" << iChip
121 << "\t" << fnOfDigisChip[iStation][iLad][iHla][iMod][iChip] << endl;
122 }
123 }
124 }
125 }
126 gDirectory -> mkdir("STSDigitizeQA");
127 gDirectory -> cd("STSDigitizeQA");
128 fHM -> WriteToFile();
129 gDirectory -> cd("../");
130 CbmSimulationReport* report = new CbmStsDigitizeQaReport(fSetup, fDigiPar);
131 report -> Create(fHM, fOutputDir);
132 delete report;
133
134 /* Double_t matchedHits = 100. * (Double_t) fHM->H1("hno_NofObjects_MatchedHits_Station_" + type)->Integral() /
135 (Double_t) fHM->H1("hno_NofObjects_Hits_Station_" + type)->Integral();
136 Double_t efficiency = 100 * (Double_t) fHM->H1("hno_NofObjects_MatchedHits_Station_" + type)->Integral() /
137 (Double_t) fHM->H1("hno_NofObjects_Points_Station_" + type)->Integral();
138 Double_t ghost = 100 * ((Double_t) fHM->H1("hno_NofObjects_Hits_Station_" + type)->Integral() -
139 (Double_t) fHM->H1("hno_NofObjects_MatchedHits_Station_" + type)->Integral()) /
140 (Double_t) fHM->H1("hno_NofObjects_Points_Station_" + type)->Integral();
141
142 std::cout<<" -I- CbmStsTimeBasedQa: Hits: "<<fHM->H1("hno_NofObjects_Hits_Station_" + type)->Integral()
143 <<"\n -I- CbmStsTimeBasedQa: MatchedHits: "<<fHM->H1("hno_NofObjects_MatchedHits_Station_" + type)->Integral()
144 <<"\n -I- CbmStsTimeBasedQa: MatchedHits: "<<matchedHits<<" %"
145 <<"\n -I- CbmStsTimeBasedQa: Efficiency : "<<efficiency<<" %"
146 <<"\n -I- CbmStsTimeBasedQa: Ghost : "<<ghost<<" %";*/
147}
148
149void CbmStsDigitizeQa::ReadDataBranches(){
150 FairRootManager* ioman = FairRootManager::Instance();
151 if ( NULL == ioman )
152 LOG(fatal) << GetName() << ": No FairRootManager!";
153
154 fStsPoints = (TClonesArray*)ioman -> GetObject("StsPoint");
155 if ( NULL == fStsPoints )
156 //LOG(error) << GetName() << ": No StsPoint array!";
157
158 fStsDigis = (TClonesArray*) ioman -> GetObject("StsDigi");
159 if ( NULL == fStsDigis )
160 LOG(error) << GetName() << ": No StsDigi array!";
161}
162
163void CbmStsDigitizeQa::CreateHistograms(){
164 CreateNofObjectsHistograms();
165 CreateDigiHistograms();
166 fHM -> Create1<TH1F>("h_EventNo_DigitizeQa", "h_EventNo_DigitizeQa", 1, 0, 1.);
167}
168
169void CbmStsDigitizeQa::CreateNofObjectsHistograms(){
170 Int_t nofBins = 100;
171 Double_t minX = -0.5;
172 Double_t maxX = 49999.5;
173 string name = "h_NofObjects_";
174 fHM -> Create1<TH1F>(name + "Points", name + "Points;Objects per event;Entries", nofBins, minX, maxX);
175 fHM -> Create1<TH1F>(name + "Digis", name + "Digis;Objects per event;Entries", nofBins, minX, maxX);
176
177 nofBins = 8;
178 minX = -0.5;
179 maxX = 7.5;
180 fHM -> Create1<TH1F>(name + "Points_Station", name + "Points_Station;Station number;Objects per event", nofBins, minX, maxX);
181 fHM -> Create1<TH1F>(name + "Digis_Station", name + "Digis_Station;Station number;Oblects per enent", nofBins, minX, maxX);
182}
183
184void CbmStsDigitizeQa::CreateDigiHistograms(){
185 Int_t nofBins = 25;
186 Double_t minX = 0.5;
187 Double_t maxX = minX + nofBins;
188 fHM -> Create1<TH1F>("h_PointsInDigi", "PointsInDigi;Number of Points;Entries", nofBins, minX, maxX);
189 fHM -> Create1<TH1F>("h_PointsInDigiLog", "PointsInDigi;Number of Points;Entries", nofBins, minX, maxX);
190 fHM -> Create1<TH1F>("h_DigisByPoint", "DigisByPoint;Number of Digis;Entries" , nofBins, minX, maxX);
191 fHM -> Create1<TH1F>("h_DigisByPointLog", "DigisByPoint;Number of Digis;Entries" , nofBins, minX, maxX);
192 nofBins = fDigiPar->GetNofAdc();
193 fHM -> Create1<TH1F>("h_DigiCharge", "DigiCharge;Digi Charge, ADC;Entries", nofBins, 0., Double_t(nofBins));
194 for (Int_t stationId = 0; stationId < fNofStation; stationId++){
195 fHM -> Create2<TH2F>(Form("h_DigisPerChip_Station%i",stationId),
196 Form("Digis per Chip, Station %i;x, cm;y, cm", stationId), 400, -50, 50, 200, -50, 50);
197 fHM -> Create2<TH2F>(Form("h_PointsMap_Station%i",stationId),
198 Form("Points Map, Station %i;x, cm;y, cm", stationId), 100, -50, 50, 100, -50, 50);
199 fHM -> Create2<TH2F>(Form("h_MeanAngleMap_Station%i",stationId),
200 Form("Mean Angle Map, Station %i;x, cm;y, cm", stationId), 50, -50, 50, 50, -50, 50);
201 fHM -> Create2<TH2F>(Form("h_RMSAngleMap_Station%i",stationId),
202 Form("RMS Angle Map, Station %i;x, cm;y, cm", stationId), 50, -50, 50, 50, -50, 50);
203 }
204 Double_t local[3] = {0.,0.,0.};
205 Double_t global[3];
206 for (Int_t moduId = 0; moduId < fSetup -> GetNofModules(); moduId++){
207 CbmStsModule * modu = static_cast<CbmStsModule*>(fSetup -> GetModule(moduId));
208 TGeoPhysicalNode * node = modu -> CbmStsElement::GetDaughter(0) -> CbmStsElement::GetPnode();
209 if ( node ){
210 TGeoMatrix * matrix = node -> GetMatrix();
211 matrix -> LocalToMaster(local, global);
212 }
213 fHM -> Create1<TH1F>(Form("h_ParticleAngles_%s", modu -> GetName()),
214 Form("Particle Angles (%.0f cm, %.0f cm);Angle, deg;Entries", global[0], global[1]), 90, 0., 90.);
215 }
216}
217
218void CbmStsDigitizeQa::ProcessDigisAndPoints(const TClonesArray* digis, const TClonesArray * points){
219 if ( NULL != digis && fHM -> Exists("h_NofObjects_Digis") ) fHM -> H1("h_NofObjects_Digis") -> Fill(digis -> GetEntriesFast());
220 std::set<Double_t> pointIndexes;
221 std::map<Double_t, Int_t> stations;
222 std::map<Double_t, Int_t> digisByPoint;
223 std::map<Double_t, Int_t>::iterator map_it;
224 pointIndexes.clear();
225 Double_t local[3] = {0.,0.,0.};
226 Double_t global[3];
227 for(Int_t iDigi = 0; iDigi < digis -> GetEntriesFast(); iDigi++) {
228 const CbmStsDigi* stsDigi = static_cast<const CbmStsDigi*>(digis -> At(iDigi));
229 const CbmMatch* digiMatch = static_cast<const CbmMatch*>(stsDigi -> GetMatch());
230 Int_t stationId = fSetup->GetStationNumber(stsDigi->GetAddress());
231 Int_t iLad = CbmStsAddress::GetElementId(stsDigi -> GetAddress(), kStsLadder);
232 Int_t iHla = CbmStsAddress::GetElementId(stsDigi -> GetAddress(), kStsHalfLadder);
233 Int_t iMod = CbmStsAddress::GetElementId(stsDigi -> GetAddress(), kStsModule);
234 CbmStsModule * modu = static_cast<CbmStsModule*>(fSetup -> GetElement(stsDigi -> GetAddress(), kStsModule));
235 Int_t nOfChannelsM = modu -> GetNofChannels();
236 TGeoPhysicalNode * node = modu -> CbmStsElement::GetDaughter(0) -> CbmStsElement::GetPnode();
237 if ( node ){
238 TGeoMatrix * matrix = node -> GetMatrix();
239 matrix -> LocalToMaster(local, global);
240 }
241
242 Int_t iChan = stsDigi->GetChannel();
243 Int_t iChip = iChan / 128 ;
244 fnOfDigisChip[stationId][iLad][iHla][iMod][iChip]++;
245 fHM -> H2(Form("h_DigisPerChip_Station%i",stationId)) -> Fill(global[0] + 50. / 400. * ((iChip-8.) * 2. - 1.), global[1]);
246
247 for(Int_t iLink = 0; iLink < digiMatch -> GetNofLinks(); iLink++) {
248 const CbmLink link = digiMatch -> GetLink(iLink);
249 Double_t index = (1000 * link.GetIndex()) + (link.GetFile()) + (0.0001 * link.GetEntry());
250 pointIndexes.insert(index);
251 stations.insert(std::pair<Double_t, Int_t>(index, stationId));
252 Int_t channel = stsDigi -> GetChannel();
253
254 Int_t side = channel < Int_t(nOfChannelsM / 2.) ? 0 : 1;
255 map_it = digisByPoint.find(index + (side * 0.00001));
256 if ( map_it != digisByPoint.end() ) {
257 map_it -> second++;
258 } else {
259 digisByPoint.insert(std::pair<Double_t, Int_t>(index + (side * 0.00001), 1));
260 }
261 }
262 fHM -> H1("h_NofObjects_Digis_Station") -> Fill(stationId);
263 fHM -> H1("h_PointsInDigi") -> Fill(digiMatch -> GetNofLinks());
264 fHM -> H1("h_PointsInDigiLog") -> Fill(digiMatch -> GetNofLinks());
265 fHM -> H1("h_DigiCharge") -> Fill(stsDigi -> GetCharge());
266 }
267 fHM -> H1("h_NofObjects_Points") -> Fill(pointIndexes.size());
268 std::set<Double_t>::iterator set_it;
269 for(set_it = pointIndexes.begin(); set_it != pointIndexes.end(); ++set_it) {
270 fHM -> H1("h_NofObjects_Points_Station") -> Fill(stations[*set_it]);
271 fHM -> H1("h_DigisByPoint") -> Fill(digisByPoint[*set_it]);
272 fHM -> H1("h_DigisByPoint") -> Fill(digisByPoint[*set_it + 0.00001]);
273 fHM -> H1("h_DigisByPointLog") -> Fill(digisByPoint[*set_it]);
274 fHM -> H1("h_DigisByPointLog") -> Fill(digisByPoint[*set_it + 0.00001]);
275 }
276 if ( pointIndexes.size() > static_cast<size_t>(fMaxScale) ) fMaxScale = pointIndexes.size();
277
278 Double_t pointX, pointY; //, pointZ;
279 Double_t pointPX, pointPZ;
280 for (Int_t iPoint = 0; iPoint < points -> GetEntriesFast(); iPoint ++){
281 const FairMCPoint * stsPoint = static_cast<const FairMCPoint*>(points -> At(iPoint));
282 CbmStsModule * modu = static_cast<CbmStsModule*>(fSetup -> GetElement(stsPoint -> GetDetectorID(), kStsModule));
283 TGeoPhysicalNode * node = modu -> CbmStsElement::GetDaughter(0) -> CbmStsElement::GetPnode();
284 if ( node ){
285 TGeoMatrix * matrix = node -> GetMatrix();
286 matrix -> LocalToMaster(local, global);
287 }
288 pointX = stsPoint -> GetX();
289 pointY = stsPoint -> GetY();
290// pointZ = stsPoint -> GetZ();
291 pointPX = stsPoint -> GetPx();
292 pointPZ = stsPoint -> GetPz();
293 Int_t stationId = fSetup->GetStationNumber(stsPoint->GetDetectorID());
294 fHM -> H2(Form("h_PointsMap_Station%i", stationId)) -> Fill(pointX, pointY);
295 fHM -> H1(Form("h_ParticleAngles_%s", modu -> GetName())) -> Fill(TMath::Abs(TMath::ATan(pointPX / pointPZ)) * 180. / 3.1416);
296 }
297}
298
299
300void CbmStsDigitizeQa::ProcessAngles(){
301 Double_t local[3] = {0.,0.,0.};
302 Double_t global[3];
303 for (Int_t iStation = 0; iStation < fNofStation; iStation++){
304 CbmStsElement * stat = fSetup -> GetDaughter(iStation);
305 for (Int_t iLad = 0; iLad < stat -> GetNofDaughters(); iLad++) {
306 CbmStsElement* ladd = stat -> GetDaughter(iLad);
307 for (Int_t iHla = 0; iHla < ladd -> GetNofDaughters(); iHla++) {
308 CbmStsElement* hlad = ladd -> GetDaughter(iHla);
309 for (Int_t iMod = 0; iMod < hlad -> GetNofDaughters(); iMod++) {
310 CbmStsElement* modu = hlad -> GetDaughter(iMod);
311 Double_t mean = fHM -> H1(Form("h_ParticleAngles_%s",modu -> GetName())) -> GetMean();
312 Double_t rms = fHM -> H1(Form("h_ParticleAngles_%s",modu -> GetName())) -> GetRMS();
313 TGeoPhysicalNode * node = modu -> CbmStsElement::GetDaughter(0) -> CbmStsElement::GetPnode();
314 if ( node ){
315 TGeoMatrix * matrix = node -> GetMatrix();
316 matrix -> LocalToMaster(local, global);
317 }
318 fHM -> H2(Form("h_MeanAngleMap_Station%i", iStation)) -> Fill(global[0],global[1],mean);
319 fHM -> H2(Form("h_RMSAngleMap_Station%i", iStation)) -> Fill(global[0],global[1],rms);
320 }
321 }
322 }
323 }
324}
fDigiPar
BmnSsdDigitizeParameters * fDigiPar
Digitisation parameters.
Definition BmnSsdDigitize.cxx:65
GetNofModules
Int_t GetNofModules(TGeoNode *station)
Definition CbmBmnStsDigitize.cxx:489
kStsModule
@ kStsModule
Module.
Definition CbmStsAddress.h:24
kStsLadder
@ kStsLadder
Ladder.
Definition CbmStsAddress.h:22
kStsHalfLadder
@ kStsHalfLadder
Halfladder.
Definition CbmStsAddress.h:23
CbmStsAddress::GetElementId
static Int_t GetElementId(UInt_t address, Int_t level)
Definition CbmStsAddress.h:103
CbmStsDigitizeQa::Exec
virtual void Exec(Option_t *opt)
Definition CbmStsDigitizeQa.cxx:95
CbmStsDigitizeQa::ProcessDigisAndPoints
void ProcessDigisAndPoints(const TClonesArray *digis, const TClonesArray *points)
Definition CbmStsDigitizeQa.cxx:218
CbmStsDigitizeQa::CbmStsDigitizeQa
CbmStsDigitizeQa(CbmStsDigitize *digitizer=NULL)
Definition CbmStsDigitizeQa.cxx:40
CbmStsDigitizeQa::SetParContainers
virtual void SetParContainers()
Definition CbmStsDigitizeQa.cxx:59
CbmStsDigitizeQa::Finish
virtual void Finish()
Definition CbmStsDigitizeQa.cxx:100
CbmStsDigitizeQa::Init
virtual InitStatus Init()
Definition CbmStsDigitizeQa.cxx:64
CbmStsElement
Class representing an element of the STS setup.
Definition CbmStsElement.h:32
CbmStsElement::GetDaughter
CbmStsElement * GetDaughter(Int_t index) const
Definition CbmStsElement.cxx:62
CbmStsModule
Class representing an instance of a readout unit in the CBM-STS.
Definition CbmStsModule.h:41
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition CbmStsSetup.cxx:90