BmnRoot
Loading...
Searching...
No Matches
CbmStsMCQa.cxx
Go to the documentation of this file.
1#include "CbmStsMCQa.h"
2
3#include "CbmHistManager.h"
4#include "CbmStsAddress.h"
5#include "CbmStsModule.h"
6#include "CbmStsElement.h"
7#include "CbmStsSetup.h"
8#include "CbmStsPoint.h"
9#include "CbmSimulationReport.h"
10//#include "CbmStsMCQaReport.h"
11
12#include "FairRootManager.h"
13#include "FairLogger.h"
14
15#include "TClonesArray.h"
16#include "TH1.h"
17#include "TF1.h"
18#include "TH1D.h"
19#include "TH2.h"
20#include "TProfile.h"
21#include "TProfile2D.h"
22
23using std::vector;
24using std::map;
25
26CbmStsMCQa::CbmStsMCQa():
27 FairTask()
28 , fHM(new CbmHistManager())
29 , fStsPoints(NULL)
30 , fMCTracks(NULL)
31 , fSetup(NULL)
32 , fNofStation(8)
33{
34}
35
36CbmStsMCQa::~CbmStsMCQa()
37{
38 if ( fHM ) delete fHM;
39}
40
41InitStatus CbmStsMCQa::Init()
42{
43 fSetup = CbmStsSetup::Instance();
44 fNofStation = fSetup -> GetNofStations();
45 LOG(info) << "Sts Setup consist of " << fNofStation << " stations.";
46
47 ReadDataBranches();
48 CreateHistograms();
49 return kSUCCESS;
50}
51
52void CbmStsMCQa::ReadDataBranches()
53{
54 FairRootManager* ioman = FairRootManager::Instance();
55 if ( NULL == ioman )
56 LOG(fatal) << "No FairRootManager!";
57
58 fStsPoints = dynamic_cast<TClonesArray*>(ioman -> GetObject("StsPoint"));
59 if ( NULL == fStsPoints )
60 //LOG(error) << "No StsPoint array!";
61
62 fMCTracks = dynamic_cast<TClonesArray*>(ioman -> GetObject("MCTrack"));
63 if ( NULL == fMCTracks )
64 //LOG(error) << "No MCTrack array!";
65}
66
67void CbmStsMCQa::CreateHistograms()
68{
69 CreateNofObjectsHistograms();
70 CreatePointHistograms();
71 fHM -> Create1<TH1F>("h_sts_EventNo_MCQa", "h_stsEventNo_MCQa", 1, 0, 1.);
72}
73
74void CbmStsMCQa::CreateNofObjectsHistograms()
75{
76
77 Int_t nofBins = 100;
78 Double_t minX = -0.5;
79 Double_t maxX = 99.5;
80 string name = "h_sts_NofObjects_";
81 fHM -> Create1<TH1F>(name + "Points", name + "Points;Objects per event;Entries", nofBins, minX, maxX);
82
83 nofBins = fNofStation;
84 minX = -0.5;
85 maxX = static_cast<Float_t>(fNofStation)-0.5;
86 fHM -> Create1<TH1F>(name + "Points_Station", name + "Points_Station;Station number;Objects per event", nofBins, minX, maxX);
87}
88
89void CbmStsMCQa::CreatePointHistograms()
90{
91 for (Int_t stationId = 0; stationId < fNofStation; stationId++){
92 Int_t nofBins = 100;
93 Double_t minX = -0.5;
94 Double_t maxX = 99.5;
95 fHM -> Create1<TH1F>(Form("h_sts_MultPoints_Station%i",stationId),
96 Form("Mult, Station %i;Objects per event;Entries", stationId),
97 nofBins, minX, maxX);
98
99 fHM -> Create2<TH2F>(Form("h_sts_PointsMap_Station%i",stationId),
100 Form("StsPoint, Station %i;x, cm;y, cm", stationId),
101 200, -10., 10., 200, -10., 10.);
102 fHM -> Create2<TH2F>(Form("h_sts_PointsMap_NoOverlap_Station%i",stationId),
103 Form("StsPoint, Station %i;x, cm;y, cm", stationId),
104 200, -10., 10., 200, -10., 10.);
105
106 fHM -> Create2<TH2F>(Form("h_sts_PointsMapEvent_Station%i",stationId),
107 Form("StsPoint/cm^{2}, Station %i;x, cm;y, cm", stationId),
108 200, -10., 10., 200, -10., 10.);
109 fHM -> Create2<TH2F>(Form("h_sts_PointsMapEvent_NoOverlap_Station%i",stationId),
110 Form("StsPoint/cm^{2}, Station %i;x, cm;y, cm", stationId),
111 200, -10., 10., 200, -10., 10.);
112
113 fHM -> Create2<TH2F>(Form("h_sts_PointsMapRate_Station%i",stationId),
114 Form("StsPoint/cm^{2}/s, Station %i;x, cm;y, cm", stationId),
115 200, -10., 10., 200, -10., 10.);
116 fHM -> Create2<TH2F>(Form("h_sts_PointsMapRate_NoOverlap_Station%i",stationId),
117 Form("StsPoint/cm^{2}/s, Station %i;x, cm;y, cm", stationId),
118 200, -10., 10., 200, -10., 10.);
119 fHM -> Create1<TH1F>(Form("h_sts_XPos_Station%i",stationId),"X position;x, cm; Entries", 200, -10., 10.);
120 fHM -> Create1<TH1F>(Form("h_sts_YPos_Station%i",stationId),"Y position;y, cm; Entries", 200, -10., 10.);
121
122
123 }
124 fHM -> Create1<TH1F>("h_sts_XMom","momentum px; GeV/c; Entries", 100, -5., 5.);
125 fHM -> Create1<TH1F>("h_sts_YMom","momentum py; GeV/c; Entries", 100, -5., 5.);
126 fHM -> Create1<TH1F>("h_sts_ZMom","momentum pz; GeV/c; Entries", 500, -10., 40.);
127 fHM -> Create1<TH1F>("h_sts_ELoss","energy loss; ; Entries", 100, 0., 0.02);
128
129 fHM -> Create1<TH1F>("h_sts_XPos","X position;x, cm; Entries", 200, -10., 10.);
130 fHM -> Create1<TH1F>("h_sts_YPos","Y position;y, cm; Entries", 200, -10., 10.);
131}
132
133void CbmStsMCQa::Exec(Option_t*){
134 ProcessPoints(fStsPoints);
135 fHM -> H1("h_sts_EventNo_MCQa") -> Fill(0.5);
136}
137
138void CbmStsMCQa::ProcessPoints(const TClonesArray * points)
139{
140
141 fHM -> H1("h_sts_NofObjects_Points") -> Fill(points -> GetEntriesFast());
142
143 Double_t pointX=0.;
144 Double_t pointY=0.;
145
146 Double_t pX=0.;
147 Double_t pY=0.;
148 Double_t pZ=0.;
149
150 std::map<Int_t, vector<Int_t>> used_map = { {0, {}}, {1, {}} };
151
152 for(Int_t iPoint = 0; iPoint < points -> GetEntriesFast(); iPoint++) {
153 const CbmStsPoint* stsPoint = static_cast<const CbmStsPoint*>(points -> At(iPoint));
154 Int_t stationId = fSetup->GetStationNumber(stsPoint->GetDetectorID());
155 fHM -> H1("h_sts_NofObjects_Points_Station") -> Fill(stationId);
156
157 pointX = stsPoint -> GetXIn();
158 pointY = stsPoint -> GetYIn();
159
160 pX = stsPoint -> GetPx();
161 pY = stsPoint -> GetPy();
162 pZ = stsPoint -> GetPz();
163
164 fHM -> H1(Form("h_sts_XPos_Station%i", stationId)) -> Fill(pointX);
165 fHM -> H1(Form("h_sts_YPos_Station%i", stationId)) -> Fill(pointY);
166
167 fHM -> H2(Form("h_sts_PointsMap_Station%i", stationId)) -> Fill(pointX, pointY);
168 fHM -> H2(Form("h_sts_PointsMapEvent_Station%i", stationId)) -> Fill(pointX, pointY);
169 fHM -> H2(Form("h_sts_PointsMapRate_Station%i", stationId)) -> Fill(pointX, pointY);
170 fHM -> H1("h_sts_XPos") -> Fill(pointX);
171 fHM -> H1("h_sts_YPos") -> Fill(pointY);
172 fHM -> H1("h_sts_XMom") -> Fill(pX);
173 fHM -> H1("h_sts_YMom") -> Fill(pY);
174 fHM -> H1("h_sts_ZMom") -> Fill(pZ);
175 fHM -> H1("h_sts_ELoss") -> Fill(stsPoint -> GetEnergyLoss());
176
177 Int_t mcTrackID = stsPoint -> GetTrackID();
178
179 if (std::find(used_map[stationId].begin(), used_map[stationId].end(), mcTrackID) == used_map[stationId].end()) {
180 used_map[stationId].push_back(mcTrackID);
181 fHM -> H2(Form("h_sts_PointsMap_NoOverlap_Station%i", stationId)) -> Fill(pointX, pointY);
182 fHM -> H2(Form("h_sts_PointsMapEvent_NoOverlap_Station%i", stationId)) -> Fill(pointX, pointY);
183 fHM -> H2(Form("h_sts_PointsMapRate_NoOverlap_Station%i", stationId)) -> Fill(pointX, pointY);
184 }
185 }
186 fHM -> H1(Form("h_sts_MultPoints_Station%i",0)) -> Fill(used_map[0].size());
187 fHM -> H1(Form("h_sts_MultPoints_Station%i",1)) -> Fill(used_map[1].size());
188
189}
190
191void CbmStsMCQa::Finish(){
192
193 Int_t nofEvents = fHM -> H1("h_sts_EventNo_MCQa") -> GetEntries();
194
195 // Do here some scaling of the histograms to have MCPoint per cm^2
196
197 Int_t xbins = (fHM->H2("h_sts_PointsMap_Station0"))->GetXaxis()->GetNbins();
198 Float_t xmax = fHM -> H2("h_sts_PointsMapEvent_Station0") -> GetXaxis()->GetXmax();
199 Float_t xmin = fHM -> H2("h_sts_PointsMapEvent_Station0") -> GetXaxis()->GetXmin();
200 Float_t scaleX = static_cast<Float_t>(xbins)/(xmax - xmin);
201
202 LOG(info) << "scaleX: " << scaleX;
203
204 Int_t ybins = fHM -> H2("h_sts_PointsMapEvent_Station0") -> GetYaxis()->GetNbins();
205 Int_t ymax = fHM -> H2("h_sts_PointsMapEvent_Station0") -> GetYaxis()->GetXmax();
206 Int_t ymin = fHM -> H2("h_sts_PointsMapEvent_Station0") -> GetYaxis()->GetXmin();
207 Float_t scaleY = static_cast<Float_t>(ybins)/(ymax - ymin);
208
209 LOG(info) << "scaleY: " << scaleY;
210
211 Float_t scale = scaleX * scaleY;
212 scale=1.;
213
214 LOG(info) << "Scale factor to cm^2: " << scale;
215
216 for(Int_t i=0; i<fNofStation; ++i) {
217 fHM -> Scale(Form("h_sts_PointsMapEvent_Station%i", i),scale/nofEvents);
218 fHM -> Scale(Form("h_sts_PointsMapEvent_NoOverlap_Station%i", i),scale/nofEvents);
219 fHM -> Scale(Form("h_sts_PointsMapRate_Station%i", i),10000000.*scale/nofEvents);
220 fHM -> Scale(Form("h_sts_PointsMapRate_NoOverlap_Station%i", i),10000000*scale/nofEvents);
221 }
222
223 gDirectory -> mkdir("QA/StsMCQa");
224 gDirectory -> cd("QA/StsMCQa");
225 fHM -> WriteToFile();
226 gDirectory -> cd("../..");
227 // CbmSimulationReport* report = new CbmStsMCQaReport(fSetup, fDigitizer);
228 // report -> Create(fHM, fOutputDir);
229 // delete report;
230
231 // Compare results with defined benchmark results and raise an ERROR if there are differences
232 // Either get default histogramms from a benchmark qa file or as parameters from the parameter container
233}
i
int i
Definition P4_F32vec4.h:22
CbmStsMCQa::CreatePointHistograms
void CreatePointHistograms()
Definition CbmStsMCQa.cxx:89
CbmStsMCQa::ProcessPoints
void ProcessPoints(const TClonesArray *)
Definition CbmStsMCQa.cxx:138
CbmStsMCQa::Finish
virtual void Finish()
Definition CbmStsMCQa.cxx:191
CbmStsMCQa::CreateHistograms
void CreateHistograms()
Definition CbmStsMCQa.cxx:67
CbmStsMCQa::CreateNofObjectsHistograms
void CreateNofObjectsHistograms()
Definition CbmStsMCQa.cxx:74
CbmStsMCQa::~CbmStsMCQa
virtual ~CbmStsMCQa()
Definition CbmStsMCQa.cxx:36
CbmStsMCQa::Exec
virtual void Exec(Option_t *)
Definition CbmStsMCQa.cxx:133
CbmStsMCQa::Init
virtual InitStatus Init()
Definition CbmStsMCQa.cxx:41
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition CbmStsSetup.cxx:90