Belle II Software development
CDCDedxValidation.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <reconstruction/modules/CDCDedxValidation/CDCDedxValidation.h>
10
11#include <mdst/dataobjects/Track.h>
12#include <mdst/dataobjects/ECLCluster.h>
13
14#include <TH2D.h>
15
16using namespace Belle2;
17
18
19REG_MODULE(CDCDedxValidation);
20
21//-------------------------------------------------
24 fD0Window(1.0),
25 fZ0Window(1.0),
26 fnRunCounter(0),
27 fiRun(0),
28 fnBinsdedx(120),
29 fnBinsdedxLE(0.4),
30 fnBinsdedxUE(1.6)
31{
32 setDescription("Make data quality monitoring plots for CDC dE/dx");
33 addParam("outputFileName", fOutFileName, "Name for output file", std::string("CDCdEdxValidation.root"));
34 addParam("SampleType", fCollType, "Switch to hadron (false) vs bhabha files", std::string("temp"));
35 addParam("fnRuns", fnRuns, "Number of input runs");
36
37}
38
39//-----------------------------------------
41{
42
43 m_cdcDedxTracks.isRequired();
44
45 if (fCollType != "radbhabha" and fCollType != "bhabha" and fCollType != "hadron") {
46 printf("Wrong input file type (%s): choose bhabha or radbhabha or hadron \n", fCollType.data());
47 return;
48 } else {
49 fBasic = new TList(); fBasic->SetOwner(); fBasic->SetName("AllRunBasics");
50 fPRdEdx = new TList(); fPRdEdx->SetOwner(); fPRdEdx->SetName("PerRunDEdx");
51 fPRdEdxinP = new TList(); fPRdEdxinP->SetOwner(); fPRdEdxinP->SetName("PerRundEdxinP");
52 DefineHistograms("AR", 0);
53 }
54}
55
56
57//---------------------------------------
59{
60
61 StoreObjPtr<EventMetaData> eventMetaDataPtr;
62 fCurrentRunNum = eventMetaDataPtr->getRun();
63
65 fcRunGain = -1.0;
66
68
70}
71
72
73//------------------------------------
75{
76
77 //Loop over CDC Tracks
78 for (Int_t idedx = 0; idedx < m_cdcDedxTracks.getEntries(); idedx++) {
79
80 CDCDedxTrack* dedxTrack = m_cdcDedxTracks[idedx];
81 if (!dedxTrack) {
82 ((TH1D*)fBasic->FindObject("hTrkPerEvtStats"))->Fill(0.0);
83 continue;
84 }
85
86 const Track* track = dedxTrack->getRelatedFrom<Track>();
87 if (!track) {
88 ((TH1D*)fBasic->FindObject("hTrkPerEvtStats"))->Fill(1.0);
89 continue;
90 }
91
92 const TrackFitResult* mTrack = NULL;
93 if (fCollType == "bhabha" || fCollType == "radbhabha") {
94 mTrack = track->getTrackFitResultWithClosestMass(Const::electron);
95 } else {
96 mTrack = track->getTrackFitResultWithClosestMass(Const::pion);
97 }
98 if (!mTrack) {
99 ((TH1D*)fBasic->FindObject("hTrkPerEvtStats"))->Fill(2.0);
100 continue;
101 }
102
103 bool IsTrkSelected = IsSelectedTrack(mTrack);
104 if (!IsTrkSelected) {
105 ((TH1D*)fBasic->FindObject("hTrkPerEvtStats"))->Fill(3.0);
106 continue;
107 }
108
109 if (dedxTrack->getNLayerHits() <= 20) continue;
110
111 const ECLCluster* eclCluster = track->getRelated<ECLCluster>();
112 if (eclCluster and eclCluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
114 if (fCollType == "bhabha" || fCollType == "radbhabha") {
115 if (abs(fTrkEoverP - 1.0) >= 0.2)continue;
116 ((TH1D*)fBasic->FindObject(Form("hEOverP_AR")))->Fill(double(fTrkEoverP));
117 }
118 } else {
119 ((TH1D*)fBasic->FindObject("hTrkPerEvtStats"))->Fill(4.0);
120 continue;
121
122 }
123
124 ((TH1D*)fBasic->FindObject("hTrkPerEvtStats"))->Fill(5.0);
125 FillHistograms(dedxTrack, mTrack);
126 }
127}
128
129//----------------------------------------------------------------------------------------------------
131{
132
133 fcRunGain = dedxTrack->getRunGain();
134
135 Int_t TrkCharge = mTrack->getChargeSign();
136 Double_t TrkdEdxnosat = dedxTrack->getDedxNoSat();
137 Double_t TrkdEdx = dedxTrack->getDedx();
138 Double_t TrkCosTheta = dedxTrack->getCosTheta();
139 Double_t TrkMom = dedxTrack->getMomentum();
140
141 if (TrkMom >= 8.00)TrkMom = 7.999;
142 Double_t BWMom = 0.250; //in GeV
143 Int_t iMomBin = Int_t(TrkMom / BWMom);
144
145 if (fCollType == "radbhabha" or fCollType == "bhabha") {
146
147 if (TrkCharge > 0) {
148 ((TH1D*)fBasic->FindObject(Form("hP_Positron_AR")))->Fill(TrkMom);
149 ((TH1D*)fBasic->FindObject(Form("hdEdx_Positron_AR")))->Fill(TrkdEdxnosat);
150 ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", iMomBin)))->Fill(TrkdEdxnosat);
151 } else if (TrkCharge < 0) {
152 ((TH1D*)fBasic->FindObject(Form("hP_Electron_AR")))->Fill(TrkMom);
153 ((TH1D*)fBasic->FindObject(Form("hdEdx_Electron_AR")))->Fill(TrkdEdxnosat);
154 ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Elec_Pbin_AR%d", iMomBin)))->Fill(TrkdEdxnosat);
155 }
156
157 ((TH1D*)fBasic->FindObject(Form("hdEdx_AR")))->Fill(double(TrkdEdxnosat));
158 ((TH2D*)fBasic->FindObject(Form("hdEdxvsPhi_AR")))->Fill(double(mTrack->getPhi()), double(TrkdEdxnosat));
159 ((TH2D*)fBasic->FindObject(Form("hPvsdEdx_AR")))->Fill(TrkMom * TrkCharge, double(TrkdEdxnosat));
160 ((TH2D*)fBasic->FindObject(Form("hPvsCosth_AR")))->Fill(TrkMom * TrkCharge, double(TrkCosTheta));
161
162 hdEdx_PR[fiRun]->Fill(double(TrkdEdxnosat));
163
164 } else if (fCollType == "hadron") {
165
166 // double ChiE = dedxTrack->getChi(0);
167 // double ChiMu = dedxTrack->getChi(1);
168 double ChiPi = dedxTrack->getChi(2);
169 double ChiK = dedxTrack->getChi(3);
170 double ChiP = dedxTrack->getChi(4);
171 // double ChiD = dedxTrack->getChi(5);
172
173 ((TH2D*)fBasic->FindObject(Form("hPvsdEdx_hadAR")))->Fill(TrkMom, TrkdEdx);
174
175 if ((TrkMom < 0.40) && (fTrkEoverP < 0.4) && (TrkdEdx < (0.6 + 0.10 / (TrkMom * TrkMom)))
176 && (TrkdEdx > (0.4 + 0.012 / (TrkMom * TrkMom)))) {
177 ((TH1D*)fBasic->FindObject(Form("hPionChiallP")))->Fill(ChiPi);
178 if (TrkMom < 0.300)((TH1D*)fBasic->FindObject(Form("hPionChiLowP")))->Fill(ChiPi);
179 else ((TH1D*)fBasic->FindObject(Form("hPionChiHighP")))->Fill(ChiPi);
180 ((TH2D*)fBasic->FindObject(Form("hPvsdEdxPion_hadAR")))->Fill(TrkMom, TrkdEdx);
181 }
182
183 if ((TrkMom < 0.40) && (TrkdEdx > 1.35) && (TrkdEdx < (0.6 + 0.40 / (TrkMom * TrkMom)))
184 && (TrkdEdx > (0.6 + 0.10 / (TrkMom * TrkMom)))) {
185 ((TH1D*)fBasic->FindObject(Form("hKaonChiallP")))->Fill(ChiK);
186 if (TrkMom < 0.350)((TH1D*)fBasic->FindObject(Form("hKaonChiLowP")))->Fill(ChiK);
187 else ((TH1D*)fBasic->FindObject(Form("hKaonChiHighP")))->Fill(ChiK);
188 ((TH2D*)fBasic->FindObject(Form("hPvsdEdxKaon_hadAR")))->Fill(TrkMom, TrkdEdx);
189 }
190
191 if ((TrkMom < 0.80) && (TrkdEdx > 1.35) && (TrkdEdx < (0.6 + 1.20 / (TrkMom * TrkMom)))
192 && (TrkdEdx > (0.6 + 0.40 / (TrkMom * TrkMom)))) {
193 ((TH1D*)fBasic->FindObject(Form("hProtonChiallP")))->Fill(ChiP);
194 if (TrkMom < 0.600)((TH1D*)fBasic->FindObject(Form("hProtonChiLowP")))->Fill(ChiP);
195 else ((TH1D*)fBasic->FindObject(Form("hProtonChiHighP")))->Fill(ChiP);
196 ((TH2D*)fBasic->FindObject(Form("hPvsdEdxProton_hadAR")))->Fill(TrkMom, TrkdEdx);
197 }
198 }
199
200}
201
202
203
204//----------------------------------
206{
207 if (fCollType != "hadron")ExtractHistograms("PR");
208}
209
210
211//----------------------------------------------------------------------------------
212void CDCDedxValidationModule::DefineHistograms(TString level = "XR", Int_t iR = 123)
213{
214 if (level == "AR") {
215
216 TH1D* hTrkPerEvtStats = new TH1D("hTrkPerEvtStats", "Track selections", 6, -0.5, 5.5);
217 hTrkPerEvtStats->GetXaxis()->SetBinLabel(1, "no dedxobject");
218 hTrkPerEvtStats->GetXaxis()->SetBinLabel(2, "no assoc-track");
219 hTrkPerEvtStats->GetXaxis()->SetBinLabel(3, "no TrackFit");
220 hTrkPerEvtStats->GetXaxis()->SetBinLabel(4, "no pass cuts");
221 hTrkPerEvtStats->GetXaxis()->SetBinLabel(5, "no eclCluster");
222 hTrkPerEvtStats->GetXaxis()->SetBinLabel(6, "Selected");
223 hTrkPerEvtStats->SetFillColor(kRed);
224 hTrkPerEvtStats->SetFillStyle(3015);
225 hTrkPerEvtStats->SetMinimum(0);
226 fBasic->Add(hTrkPerEvtStats);
227
228 if (fCollType == "radbhabha" or fCollType == "bhabha") {
229
230 TH1D* hdEdx_AR = new TH1D("hdEdx_AR", "dE/dx (nohad sat)", fnBinsdedx, fnBinsdedxLE, fnBinsdedxUE);
231 hdEdx_AR->GetXaxis()->SetTitle(Form("dE/dx truncMean of %s tracks", fCollType.data()));
232 hdEdx_AR->GetYaxis()->SetTitle("Entries");
233 fBasic->Add(hdEdx_AR);
234
235 TH1D* hEOverP_AR = new TH1D("hEOverP_AR", "E/p distribution", 100, 0.5, 1.5);
236 hEOverP_AR->GetXaxis()->SetTitle("E/p distribution");
237 hEOverP_AR->GetYaxis()->SetTitle("Entries");
238 fBasic->Add(hEOverP_AR);
239
240 TH1D* hRunGainPR = new TH1D("hRunGainPR", "bla-bla", fnRuns, -0.5, fnRuns - 0.5);
241 hRunGainPR->SetTitle("Run gain variation vs. RunNumber;Run Numbers;dE/dx mean");
242 hRunGainPR->GetYaxis()->SetRangeUser(0.85, 1.15);
243 fBasic->Add(hRunGainPR);
244
245 TH1D* hP_Electron_AR = new TH1D("hP_Electron_AR", "bla-bla", 320, 0.0, 8.0);
246 hP_Electron_AR->SetTitle("Momentum distribution of e-; Momentum of (e-); Entries");
247 fBasic->Add(hP_Electron_AR);
248
249 TH1D* hP_Positron_AR = new TH1D("hP_Positron_AR", "bla-bla", 320, 0.0, 8.0);
250 hP_Positron_AR->SetTitle("Momentum distribution of e+;Momentum of (e+);Entries");
251 fBasic->Add(hP_Positron_AR);
252
253 TH1D* hdEdx_Electron_AR = new TH1D("hdEdx_Electron_AR", "bla-bla", fnBinsdedx, fnBinsdedxLE, fnBinsdedxUE);
254 hdEdx_Electron_AR->SetTitle("dE/dx (nohad sat) of e- ;dE/dx distribution (e-);Entries");
255 fBasic->Add(hdEdx_Electron_AR);
256
257 TH1D* hdEdx_Positron_AR = new TH1D("hdEdx_Positron_AR", "bla-bla", fnBinsdedx, fnBinsdedxLE, fnBinsdedxUE);
258 hdEdx_Positron_AR->SetTitle("dE/dx (nohad sat) of e+;dE/dx distribution (e+);Entries");
259 fBasic->Add(hdEdx_Positron_AR);
260
261 TH2D* hPvsdEdx_AR = new TH2D("hPvsdEdx_AR", "bla-bla", 320, -8.0, 8.0, 100, 0.0, 2.0);
262 hPvsdEdx_AR->SetTitle("dE/dx band plots for e+ and e-; Momentum of (e+(right)) and e-);dE/dx");
263 fBasic->Add(hPvsdEdx_AR);
264
265 TH2D* hdEdxvsPhi_AR = new TH2D("hdEdxvsPhi_AR", "dE/dx (no had sat) vs #phi", 64, -3.14, 3.14, 200, 0., 2.0);
266 hdEdxvsPhi_AR->SetTitle("dE/dx (no Had Sat) vs #phi;track #phi;dE/dx");
267 fBasic->Add(hdEdxvsPhi_AR);
268
269 TH2D* hPvsCosth_AR = new TH2D("hPvsCosth_AR", "cos(#theta) vs. p: all Runs", 2 * 48, -10., 10., 60, -1.2, 1.2);
270 hPvsCosth_AR->GetXaxis()->SetTitle(Form("Momentum of %s tracks", fCollType.data()));
271 hPvsCosth_AR->GetYaxis()->SetTitle("cos(#theta)");
272 fBasic->Add(hPvsCosth_AR);
273
274
275 TH1D* hdEdx_Posi_Pbin_AR[32], *hdEdx_Elec_Pbin_AR[32];
276 for (int ip = 0; ip < 32; ip++) {
277
278 hdEdx_Posi_Pbin_AR[ip] = new TH1D(Form("hdEdx_Posi_Pbin_AR%d", ip), Form("hdEdx_Posi_Pbin_AR%d", ip), fnBinsdedx, fnBinsdedxLE,
280 hdEdx_Posi_Pbin_AR[ip]->GetXaxis()->SetTitle("dE/dx distribution (e+)");
281 hdEdx_Posi_Pbin_AR[ip]->GetYaxis()->SetTitle("Entries");
282 hdEdx_Posi_Pbin_AR[ip]->SetTitle(Form("Momentum range %0.03f to %0.03f", ip * 0.250, (ip + 1) * 0.250));
283 fPRdEdxinP->Add(hdEdx_Posi_Pbin_AR[ip]);
284
285 hdEdx_Elec_Pbin_AR[ip] = new TH1D(Form("hdEdx_Elec_Pbin_AR%d", ip), Form("hdEdx_Elec_Pbin_AR%d", ip), fnBinsdedx, fnBinsdedxLE,
287 hdEdx_Elec_Pbin_AR[ip]->GetXaxis()->SetTitle("dE/dx distribution (e-)");
288 hdEdx_Elec_Pbin_AR[ip]->GetYaxis()->SetTitle("Entries");
289 hdEdx_Elec_Pbin_AR[ip]->SetTitle(Form("Momentum range %0.03f to %0.03f", ip * 0.250, (ip + 1) * 0.250));
290
291 }
292
293 for (int ip = 0; ip < 32; ip++) fPRdEdxinP->Add(hdEdx_Elec_Pbin_AR[ip]); //Adding later for simplicity
294
295 hdEdx_PR.reserve(fnRuns);
296 }
297
298 if (fCollType == "hadron") {
299
300 TH2D* hPvsdEdx_hadAR = new TH2D("hPvsdEdx_hadAR", "bla-bla", 500, 0.10, 15.0, 750, 0.05, 15);
301 hPvsdEdx_hadAR->SetTitle("dE/dx band plot; Momentum;dE/dx");
302 fBasic->Add(hPvsdEdx_hadAR);
303
304 TH2D* hPvsdEdxPion_hadAR = new TH2D("hPvsdEdxPion_hadAR", "bla-bla", 500, 0.10, 15.0, 750, 0.05, 15);
305 hPvsdEdxPion_hadAR->SetTitle("dE/dx band plot (Pion); Momentum;dE/dx");
306 hPvsdEdxPion_hadAR->SetMarkerColor(kRed);
307 fBasic->Add(hPvsdEdxPion_hadAR);
308
309 TH2D* hPvsdEdxKaon_hadAR = new TH2D("hPvsdEdxKaon_hadAR", "bla-bla", 500, 0.10, 15.0, 750, 0.05, 15);
310 hPvsdEdxKaon_hadAR->SetTitle("dE/dx band plot (Kaon); Momentum;dE/dx");
311 hPvsdEdxKaon_hadAR->SetMarkerColor(kGreen);
312 fBasic->Add(hPvsdEdxKaon_hadAR);
313
314 TH2D* hPvsdEdxProton_hadAR = new TH2D("hPvsdEdxProton_hadAR", "bla-bla", 500, 0.10, 15.0, 750, 0.05, 15);
315 hPvsdEdxProton_hadAR->SetTitle("dE/dx band plot (Proton); Momentum;dE/dx");
316 hPvsdEdxKaon_hadAR->SetMarkerColor(kBlue);
317 fBasic->Add(hPvsdEdxProton_hadAR);
318
319 //Pions chi values
320 TH1D* hPionChiallP = new TH1D("hPionChiallP", "bla-bla", 240, -6.0, 6.0);
321 hPionChiallP->SetTitle("Chi value (Pion);chi value; Entries");
322 fBasic->Add(hPionChiallP);
323
324 TH1D* hPionChiLowP = new TH1D("hPionChiLowP", "bla-bla", 240, -6.0, 6.0);
325 hPionChiLowP->SetTitle("Chi value (Pion), Momentum (0-300) MeV; chi value; Entries");
326 fBasic->Add(hPionChiLowP);
327
328 TH1D* hPionChiHighP = new TH1D("hPionChiHighP", "bla-bla", 240, -6.0, 6.0);
329 hPionChiHighP->SetTitle("Chi value (Pion), Momentum (300-400) MeV; chi value; Entries");
330 fBasic->Add(hPionChiHighP);
331
332 //Kaons chi values
333 TH1D* hKaonChiallP = new TH1D("hKaonChiallP", "bla-bla", 240, -6.0, 6.0);
334 hKaonChiallP->SetTitle("Chi value (Kaon);chi value; Entries");
335 fBasic->Add(hKaonChiallP);
336
337 TH1D* hKaonChiLowP = new TH1D("hKaonChiLowP", "bla-bla", 240, -6.0, 6.0);
338 hKaonChiLowP->SetTitle("Chi value (Kaon), Momentum (0-350) MeV; chi value; Entries");
339 fBasic->Add(hKaonChiLowP);
340
341 TH1D* hKaonChiHighP = new TH1D("hKaonChiHighP", "bla-bla", 240, -6.0, 6.0);
342 hKaonChiHighP->SetTitle("Chi value (Kaon), Momentum (350-800) MeV; chi value; Entries");
343 fBasic->Add(hKaonChiHighP);
344
345 //Protons chi values
346 TH1D* hProtonChiallP = new TH1D("hProtonChiallP", "bla-bla", 240, -6.0, 6.0);
347 hProtonChiallP->SetTitle("Chi value (Proton);chi value; Entries");
348 fBasic->Add(hProtonChiallP);
349
350 TH1D* hProtonChiLowP = new TH1D("hProtonChiLowP", "bla-bla", 240, -6.0, 6.0);
351 hProtonChiLowP->SetTitle("Chi value (Proton), Momentum (0-600) MeV; chi value; Entries");
352 fBasic->Add(hProtonChiLowP);
353
354 TH1D* hProtonChiHighP = new TH1D("hProtonChiHighP", "bla-bla", 240, -6.0, 6.0);
355 hProtonChiHighP->SetTitle("Chi value (Proton), Momentum (600-800) MeV; chi value; Entries");
356 fBasic->Add(hProtonChiHighP);
357 }
358 } else if (level == "PR") {
359 if (fCollType != "hadron") {
360 hdEdx_PR[iR] = new TH1D(Form("hdEdx_Run%d", fCurrentRunNum), Form("dE/dx (no had sat): Run # = %d", fCurrentRunNum), fnBinsdedx,
362 hdEdx_PR[iR]->GetXaxis()->SetTitle(Form("dE/dx trucMean of %s tracks", fCollType.data()));
363 hdEdx_PR[iR]->GetYaxis()->SetTitle("Entries");
364 fPRdEdx->Add(hdEdx_PR[iR]);
365 }
366 } else {
367 B2ERROR("Run Gain: Enter AR or PR mode only");
368 }
369}
370
371
372//-----------------------------------------------------------------------
374{
375
376 if (level == "PR") {
377
378 Double_t mean = 0., meanError = 0.;
379 Double_t sigma = 0., sigmaError = 0.;
380
381 if (hdEdx_PR[fiRun]->GetEntries() > 100) {
382 Int_t fitStatus = -1;
383 fitStatus = hdEdx_PR[fiRun]->Fit("gaus", "Q"); //Q = No printing
384 if (fitStatus == 0) {
385 TF1* fit = (TF1*)hdEdx_PR[fiRun]->GetFunction("gaus");
386 mean = fit->GetParameter(1);
387 meanError = fit->GetParError(1);
388 sigma = fit->GetParameter(2);
389 sigmaError = fit->GetParError(2);
390 hdEdx_PR[fiRun]->GetXaxis()->SetRangeUser(mean - 7 * sigma, mean + 7 * sigma);
391 fit->Clear();
392 }
393 }
394
395 TotRunN.push_back(fCurrentRunNum);
396 TotMean.push_back(mean);
397 TotMeanE.push_back(meanError);
398 TotSigma.push_back(sigma);
399 TotSigmaE.push_back(sigmaError);
400
401 if (fiRun % 10 == 0)((TH1D*)fBasic->FindObject("hRunGainPR"))->GetXaxis()->SetBinLabel(fiRun + 1, Form("%d", TotRunN.at(fiRun)));
402 ((TH1D*)fBasic->FindObject("hRunGainPR"))->SetBinContent(fiRun + 1, fcRunGain);
403 ((TH1D*)fBasic->FindObject("hRunGainPR"))->SetBinError(fiRun + 1, 0.001 * fcRunGain); // no meaning but histogramming only
404
405 } else if (level == "AR") {
406
407 const Int_t allNRun = TotMean.size();
408 ((TH1D*)fBasic->FindObject("hRunGainPR"))->GetXaxis()->SetRange(1, allNRun);
409
410 TH1D* hFitdEdxMeanPR = new TH1D("hFitdEdxMeanPR", "dE/dx(nohad-sat) #mu via fit vs. Runs", allNRun, 0, allNRun);
411 hFitdEdxMeanPR->GetYaxis()->SetRangeUser(0.90, 1.10);
412 hFitdEdxMeanPR->SetMarkerStyle(21);
413 hFitdEdxMeanPR->SetMarkerColor(kRed);
414 hFitdEdxMeanPR->SetMarkerSize(1);
415 hFitdEdxMeanPR->GetXaxis()->SetTitle("Run numbers");
416 hFitdEdxMeanPR->GetYaxis()->SetTitle("dEdx mean (fit)");
417 hFitdEdxMeanPR->GetXaxis()->LabelsOption("v");
418
419 TH1D* hFitdEdxSigmaPR = new TH1D("hFitdEdxSigmaPR", "dE/dx(nohad-sat) #sigma via fit vs. Runs", allNRun, 0, allNRun);
420 hFitdEdxSigmaPR->GetYaxis()->SetRangeUser(0, 0.30);
421 hFitdEdxSigmaPR->SetMarkerStyle(21);
422 hFitdEdxSigmaPR->SetMarkerColor(kRed);
423 hFitdEdxSigmaPR->SetMarkerSize(1);
424 hFitdEdxSigmaPR->GetXaxis()->SetTitle("Run numbers");
425 hFitdEdxSigmaPR->GetYaxis()->SetTitle("dEdx sigma (fit)");
426 hFitdEdxSigmaPR->GetXaxis()->LabelsOption("v");
427
428 for (Int_t i = 0; i < allNRun; i++) {
429
430 if (i % 10 == 0)hFitdEdxMeanPR->GetXaxis()->SetBinLabel(i + 1, Form("%d", TotRunN.at(i)));
431 hFitdEdxMeanPR->SetBinContent(i + 1, TotMean.at(i));
432 hFitdEdxMeanPR->SetBinError(i + 1, TotMeanE.at(i));
433
434 if (i % 10 == 0)hFitdEdxSigmaPR->GetXaxis()->SetBinLabel(i + 1, Form("%d", TotRunN.at(i)));
435 hFitdEdxSigmaPR->SetBinContent(i + 1, TotSigma.at(i));
436 hFitdEdxSigmaPR->SetBinError(i + 1, TotSigmaE.at(i));
437 }
438
439 fBasic->Add(hFitdEdxMeanPR);
440 fBasic->Add(hFitdEdxSigmaPR);
441
442 TH1D* hdEdxFit_allRun = (TH1D*)(fBasic->FindObject(Form("hdEdx_AR"))->Clone("hdEdxFit_allRun"));
443 if (hdEdxFit_allRun->GetEntries() > 100) {
444 hdEdxFit_allRun->Fit("gaus", "Q");
445 TF1* hGfit = (TF1*)hdEdxFit_allRun->GetFunction("gaus");
446 Double_t meanGFit = hGfit->GetParameter(1);
447 Double_t sigmaGFit = hGfit->GetParameter(2);
448 hdEdxFit_allRun->GetXaxis()->SetRangeUser(meanGFit - 7 * sigmaGFit, meanGFit + 7 * sigmaGFit);
449 hdEdxFit_allRun->SetFillColor(kYellow);
450 hdEdxFit_allRun->SetStats(kTRUE);
451 fBasic->Add(hdEdxFit_allRun);
452 }
453
454 TH1D* hdEdxMeanVsMomentum = new TH1D("hdEdxMeanVsMomentum", "dEdx-mean vs P bins (BW = 250MeV)", 64, -8.0, 8.0);
455 hdEdxMeanVsMomentum->GetXaxis()->SetTitle("Track Momentum");
456 hdEdxMeanVsMomentum->GetYaxis()->SetTitle("dEdx Mean");
457
458 TH1D* hdEdxSigmaVsMomentum = new TH1D("hdEdxSigmaVsMomentum", "dEdx-sigma vs P bins (BW = 250MeV)", 64, -8.0, 8.0);
459 hdEdxSigmaVsMomentum->GetXaxis()->SetTitle("Track Momentum");
460 hdEdxSigmaVsMomentum->GetYaxis()->SetTitle("dEdx Sigma");
461
462 for (int ip = 0; ip < 32; ip++) {
463 Int_t nTrack = ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", ip)))->GetEntries();
464 ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", ip)))->SetFillColor(kYellow);
465 Double_t iPMean = 1.0, iPSigma = 0.0;
466 if (nTrack > 100) {
467 ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", ip)))->Fit("gaus", "0");
468 iPMean = ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", ip)))->GetFunction("gaus")->GetParameter(1);
469 iPSigma = ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", ip)))->GetFunction("gaus")->GetParameter(2);
470 }
471 hdEdxMeanVsMomentum->SetBinContent(32 + ip + 1, iPMean);
472 hdEdxSigmaVsMomentum->SetBinContent(32 + ip + 1, iPSigma);
473 }
474
475 for (int ip = 0; ip < 32; ip++) {
476 Int_t nTrack = ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Posi_Pbin_AR%d", ip)))->GetEntries();
477 ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Elec_Pbin_AR%d", ip)))->SetFillColor(kYellow);
478 Double_t iPMean = 1.0, iPSigma = 0.0;
479
480 if (nTrack > 100) {
481 ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Elec_Pbin_AR%d", ip)))->Fit("gaus", "0");
482 iPMean = ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Elec_Pbin_AR%d", ip)))->GetFunction("gaus")->GetParameter(1);
483 iPSigma = ((TH1D*)fPRdEdxinP->FindObject(Form("hdEdx_Elec_Pbin_AR%d", ip)))->GetFunction("gaus")->GetParameter(2);
484 }
485 hdEdxMeanVsMomentum->SetBinContent(32 - ip, iPMean);
486 hdEdxSigmaVsMomentum->SetBinContent(32 - ip, iPSigma);
487 }
488
489 fBasic->Add(hdEdxSigmaVsMomentum);
490 fBasic->Add(hdEdxMeanVsMomentum);
491 } else {
492 B2ERROR("RunGain >> NO-REQUEST-FOUND for PR or AR level plots, exiting..");
493 }
494}
495
496//-----------------------------------------
498{
499
500 B2INFO("Terminating plots for all runs ");
501 fFileOutput = new TFile(Form("fvalidate%s", fOutFileName.data()), "RECREATE");
502 if (fCollType != "hadron")ExtractHistograms("AR");
503 fFileOutput->cd();
504 fBasic->Write("ARBasics", 1);
505 if (fCollType == "radbhabha" or fCollType == "bhabha") {
506 fPRdEdx->Write("PRdedx", 1);
507 fPRdEdxinP->Write("MeanSigmavsP", 1);
508 }
509 fFileOutput->Close();
510}
511
512//-----------------------------------------
514{
515
516 if (std::abs(mTrack->getD0()) >= fD0Window || std::abs(mTrack->getZ0()) >= fZ0Window)return kFALSE;
517 return kTRUE;
518
519}
Debug output for CDCDedxPID module.
Definition: CDCDedxTrack.h:25
double getDedx() const
Get dE/dx truncated mean for this track.
Definition: CDCDedxTrack.h:103
int getNLayerHits() const
Return the number of layer hits for this track.
Definition: CDCDedxTrack.h:174
double getCosTheta() const
Return cos(theta) for this track.
Definition: CDCDedxTrack.h:121
double getRunGain() const
Return the run gain for this track.
Definition: CDCDedxTrack.h:142
double getChi(int i) const
Return the PID (chi) value.
Definition: CDCDedxTrack.h:279
double getDedxNoSat() const
Get dE/dx truncated mean without the saturation correction for this track.
Definition: CDCDedxTrack.h:106
double getMomentum() const
Return the track momentum valid in the CDC.
Definition: CDCDedxTrack.h:118
std::vector< Double_t > TotMean
Mean of dedx by Fit.
Double_t fnBinsdedxUE
up edge of dedx
void FillHistograms(CDCDedxTrack *dedxTrack, const TrackFitResult *mTrack)
Filling histograms This will fill histogram defined histograms in above function.
TFile * fFileOutput
Write final objects to file for RG.
std::string fOutFileName
name of output ROOT file
Double_t fcRunGain
existing run gain
Int_t fCurrentRunNum
current run number
Int_t fnRunCounter
Total runs used counter.
TList * fBasic
List of basic histos.
virtual void initialize() override
Initialize This is inherited from base class.
std::vector< TH1D * > hdEdx_PR
histogram array per run
virtual void event() override
function to execute event (event by event) This is inherited from base class
std::vector< Double_t > TotSigma
Sigma of dedx by Fit.
CDCDedxValidationModule()
Default constructor contain list of members with initial values.
TList * fPRdEdxinP
list per run dedx in P histos
Double_t fTrkEoverP
E/p ratio for cut.
virtual void endRun() override
function is called after each event This is inherited from base class
std::vector< Int_t > TotRunN
veector array of runs processed
virtual void terminate() override
Terminate after all data processed This is inherited from base class.
std::vector< Double_t > TotSigmaE
Sigma Error of dedx by Fit.
std::vector< Double_t > TotMeanE
Mean Error of dedx by Fit.
virtual void beginRun() override
Function to execute each run This is inherited from base class.
StoreArray< CDCDedxTrack > m_cdcDedxTracks
Data members for objects, cuts and others.
Double_t fnBinsdedxLE
low edge of dedx
TList * fPRdEdx
List of per run dedx histos.
void ExtractHistograms(TString level)
Extracting histogram and some calculation Higher level histograms are filled after each run or full p...
void DefineHistograms(TString level, Int_t iR)
Definition of histograms This contain a list of histogram for validation.
Int_t fnRuns
Number of runs ref.
Int_t fiRun
Current run number.
Int_t fnBinsdedx
nbin of dedx range
std::string fCollType
collision type
Bool_t IsSelectedTrack(const TrackFitResult *mTrack)
Track selection A clean way to implement selections on tracks (so far few only)
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
static const ChargedStable electron
electron particle
Definition: Const.h:659
ECL cluster data.
Definition: ECLCluster.h:27
bool hasHypothesis(EHypothesisBit bitmask) const
Return if specific hypothesis bit is set.
Definition: ECLCluster.h:351
double getEnergy(EHypothesisBit hypothesis) const
Return Energy (GeV).
Definition: ECLCluster.cc:23
@ c_nPhotons
CR is split into n photons (N1)
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
FROM * getRelatedFrom(const std::string &name="", const std::string &namedRelation="") const
Get the object from which this object has a relation.
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
Values of the result of a track fit with a given particle hypothesis.
double getPhi() const
Getter for phi0 with CDF naming convention.
short getChargeSign() const
Return track charge (1 or -1).
double getD0() const
Getter for d0.
double getZ0() const
Getter for z0.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Class that bundles various TrackFitResults.
Definition: Track.h:25
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.