Belle II Software development
b2klm-likelihood-parameters.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/* KLM headers. */
10#include <klm/dataobjects/KLMMuidLikelihood.h>
11#include <klm/muid/MuidElementNumbers.h>
12
13/* Basf2 headers. */
14#include <analysis/dataobjects/Particle.h>
15#include <analysis/variables/AcceptanceVariables.h>
16#include <mdst/dataobjects/MCParticle.h>
17
18/* ROOT headers. */
19#include <Math/Vector4D.h>
20#include <TCanvas.h>
21#include <TChain.h>
22#include <TChainElement.h>
23#include <TF1.h>
24#include <TFile.h>
25#include <TH1D.h>
26#include <TMath.h>
27#include <TROOT.h>
28#include <TStyle.h>
29
30/* C++ headers. */
31#include <cmath>
32#include <cstdio>
33#include <iostream>
34
35#define NOUTCOME 67
36
37using namespace Belle2;
38
39TStyle* Belle2Style()
40{
41 TStyle* belle2Style = new TStyle("BELLE2", "Belle2 style");
42
43 // use plain black on white colors
44 Int_t icol = 0; // WHITE
45 belle2Style->SetFrameBorderMode(icol);
46 belle2Style->SetFrameFillColor(icol);
47 belle2Style->SetCanvasBorderMode(icol);
48 belle2Style->SetCanvasColor(icol);
49 belle2Style->SetPadBorderMode(icol);
50 belle2Style->SetPadColor(icol);
51 belle2Style->SetStatColor(icol);
52 //belle2Style->SetFillColor(icol); // don't use: white fill color for *all* objects
53
54 // set the paper & margin sizes
55 belle2Style->SetPaperSize(20, 26);
56
57 // set margin sizes
58 belle2Style->SetPadTopMargin(0.05);
59 belle2Style->SetPadRightMargin(0.05);
60 belle2Style->SetPadBottomMargin(0.16);
61 belle2Style->SetPadLeftMargin(0.16);
62
63 // set title offsets (for axis label)
64 belle2Style->SetTitleXOffset(1.4);
65 belle2Style->SetTitleYOffset(1.4);
66
67 // use large fonts
68 //Int_t font=72; // Helvetica italics
69 Int_t font = 42; // Helvetica
70 Double_t tsize = 0.05;
71 belle2Style->SetTextFont(font);
72
73 belle2Style->SetTextSize(tsize);
74 belle2Style->SetLabelFont(font, "x");
75 belle2Style->SetTitleFont(font, "x");
76 belle2Style->SetLabelFont(font, "y");
77 belle2Style->SetTitleFont(font, "y");
78 belle2Style->SetLabelFont(font, "z");
79 belle2Style->SetTitleFont(font, "z");
80
81 belle2Style->SetLabelSize(tsize, "x");
82 belle2Style->SetTitleSize(tsize, "x");
83 belle2Style->SetLabelSize(tsize, "y");
84 belle2Style->SetTitleSize(tsize, "y");
85 belle2Style->SetLabelSize(tsize, "z");
86 belle2Style->SetTitleSize(tsize, "z");
87
88 // use bold lines and markers
89 belle2Style->SetMarkerStyle(20);
90 belle2Style->SetMarkerSize(1.2);
91 belle2Style->SetHistLineWidth(2.);
92 belle2Style->SetLineStyleString(2, "[12 12]"); // postscript dashes
93
94 // get rid of X error bars
95 //belle2Style->SetErrorX(0.001);
96 // get rid of error bar caps
97 belle2Style->SetEndErrorSize(0.);
98
99 // do not display any of the standard histogram decorations
100 belle2Style->SetOptTitle(0);
101 //belle2Style->SetOptStat(1111);
102 belle2Style->SetOptStat(0);
103 //belle2Style->SetOptFit(1111);
104 belle2Style->SetOptFit(0);
105
106 // put tick marks on top and RHS of plots
107 belle2Style->SetPadTickX(1);
108 belle2Style->SetPadTickY(1);
109
110 return belle2Style;
111
112}
113
114void SetBelle2Style()
115{
116 const static TStyle* belle2Style = 0;
117 std::cout << "\nApplying BELLE2 style settings...\n" << std::endl ;
118 if (belle2Style == 0) belle2Style = Belle2Style();
119 gROOT->SetStyle("BELLE2");
120 gROOT->ForceStyle();
121}
122
123void makeprob(TChain* chain, FILE* output, const char* xmllabel)
124{
125
126 //Check first 2 char of the strings if equal (==0) or not
127 bool isMuon = (strncmp(xmllabel, "Mu", 2) == 0);
128 bool isPion = (strncmp(xmllabel, "Pi", 2) == 0);
129 bool isKaon = (strncmp(xmllabel, "Ka", 2) == 0);
130 bool isElec = (strncmp(xmllabel, "El", 2) == 0) || (strncmp(xmllabel, "Po", 2) == 0);
131 bool isProt = (strncmp(xmllabel, "Pr", 2) == 0) || (strncmp(xmllabel, "An", 2) == 0);
132 bool isDeut = (strncmp(xmllabel, "De", 2) == 0);
133
134 char label[80], title[80];
135 const char detectorName[3][40] = { "BarrelAndEndcap", "BarrelOnly", "EndcapOnly" };
136 //List of all the possible outcome and relative explanations
137 const char outcomeName[NOUTCOME][40] = { "",
138 "Barrel stop",
139 "Fwd Endcap stop",
140 "Barrel exit",
141 "Fwd Endcap exit",
142 "Bwd Endcap stop",
143 "Bwd Endcap exit",
144 "Fwd Endcap stop B0+E",
145 "Fwd Endcap stop B1+E",
146 "Fwd Endcap stop B2+E",
147 "Fwd Endcap stop B3+E",
148 "Fwd Endcap stop B4+E",
149 "Fwd Endcap stop B5+E",
150 "Fwd Endcap stop B6+E",
151 "Fwd Endcap stop B7+E",
152 "Fwd Endcap stop B8+E",
153 "Fwd Endcap stop B9+E",
154 "Fwd Endcap stop B10+E",
155 "Fwd Endcap stop B11+E",
156 "Fwd Endcap stop B12+E",
157 "Fwd Endcap stop B13+E",
158 "Fwd Endcap stop B14+E",
159 "Bwd Endcap stop B0+E",
160 "Bwd Endcap stop B1+E",
161 "Bwd Endcap stop B2+E",
162 "Bwd Endcap stop B3+E",
163 "Bwd Endcap stop B4+E",
164 "Bwd Endcap stop B5+E",
165 "Bwd Endcap stop B6+E",
166 "Bwd Endcap stop B7+E",
167 "Bwd Endcap stop B8+E",
168 "Bwd Endcap stop B9+E",
169 "Bwd Endcap stop B10+E",
170 "Bwd Endcap stop B11+E",
171 "Bwd Endcap stop B12+E",
172 "Bwd Endcap stop B13+E",
173 "Bwd Endcap stop B14+E",
174 "Fwd Endcap exit B0+E",
175 "Fwd Endcap exit B1+E",
176 "Fwd Endcap exit B2+E",
177 "Fwd Endcap exit B3+E",
178 "Fwd Endcap exit B4+E",
179 "Fwd Endcap exit B5+E",
180 "Fwd Endcap exit B6+E",
181 "Fwd Endcap exit B7+E",
182 "Fwd Endcap exit B8+E",
183 "Fwd Endcap exit B9+E",
184 "Fwd Endcap exit B10+E",
185 "Fwd Endcap exit B11+E",
186 "Fwd Endcap exit B12+E",
187 "Fwd Endcap exit B13+E",
188 "Fwd Endcap exit B14+E",
189 "Bwd Endcap exit B0+E",
190 "Bwd Endcap exit B1+E",
191 "Bwd Endcap exit B2+E",
192 "Bwd Endcap exit B3+E",
193 "Bwd Endcap exit B4+E",
194 "Bwd Endcap exit B5+E",
195 "Bwd Endcap exit B6+E",
196 "Bwd Endcap exit B7+E",
197 "Bwd Endcap exit B8+E",
198 "Bwd Endcap exit B9+E",
199 "Bwd Endcap exit B10+E",
200 "Bwd Endcap exit B11+E",
201 "Bwd Endcap exit B12+E",
202 "Bwd Endcap exit B13+E",
203 "Bwd Endcap exit B14+E"
204 };
205 const char outcomeComment[NOUTCOME][100] = { "",
206 "soft tracks that range out (stop) within the barrel KLM", // outcome=1
207 "soft tracks that range out (stop) within the forward endcap KLM (no barrel hits)", // outcome=2
208 "hard tracks that escape from the barrel KLM", // outcome=3
209 "hard tracks that escape from the forward endcap KLM (no barrel hits)", // outcome=4
210 "soft tracks that range out (stop) within the backward endcap KLM (no barrel hits)", // outcome=5
211 "hard tracks that escape from the backward endcap KLM (no barrel hits)", // outcome=6
212 "soft tracks that range out (stop) within the forward endcap KLM with B0 hit(s)", // outcome=7
213 "soft tracks that range out (stop) within the forward endcap KLM with B1 hit(s)", // outcome=8
214 "soft tracks that range out (stop) within the forward endcap KLM with B2 hit(s)", // outcome=9
215 "soft tracks that range out (stop) within the forward endcap KLM with B3 hit(s)", // outcome=10
216 "soft tracks that range out (stop) within the forward endcap KLM with B4 hit(s)", // outcome=11
217 "soft tracks that range out (stop) within the forward endcap KLM with B5 hit(s)", // outcome=12
218 "soft tracks that range out (stop) within the forward endcap KLM with B6 hit(s)", // outcome=13
219 "soft tracks that range out (stop) within the forward endcap KLM with B7 hit(s)", // outcome=14
220 "soft tracks that range out (stop) within the forward endcap KLM with B8 hit(s)", // outcome=15
221 "soft tracks that range out (stop) within the forward endcap KLM with B9 hit(s)", // outcome=16
222 "soft tracks that range out (stop) within the forward endcap KLM with B10 hit(s)", // outcome=17
223 "soft tracks that range out (stop) within the forward endcap KLM with B11 hit(s)", // outcome=18
224 "soft tracks that range out (stop) within the forward endcap KLM with B12 hit(s)", // outcome=19
225 "soft tracks that range out (stop) within the forward endcap KLM with B13 hit(s)", // outcome=20
226 "soft tracks that range out (stop) within the forward endcap KLM with B14 hit(s)", // outcome=21
227 "soft tracks that range out (stop) within the backward endcap KLM with B0 hit(s)", // outcome=22
228 "soft tracks that range out (stop) within the backward endcap KLM with B1 hit(s)", // outcome=23
229 "soft tracks that range out (stop) within the backward endcap KLM with B2 hit(s)", // outcome=24
230 "soft tracks that range out (stop) within the backward endcap KLM with B3 hit(s)", // outcome=25
231 "soft tracks that range out (stop) within the backward endcap KLM with B4 hit(s)", // outcome=26
232 "soft tracks that range out (stop) within the backward endcap KLM with B5 hit(s)", // outcome=27
233 "soft tracks that range out (stop) within the backward endcap KLM with B6 hit(s)", // outcome=28
234 "soft tracks that range out (stop) within the backward endcap KLM with B7 hit(s)", // outcome=29
235 "soft tracks that range out (stop) within the backward endcap KLM with B8 hit(s)", // outcome=30
236 "soft tracks that range out (stop) within the backward endcap KLM with B9 hit(s)", // outcome=31
237 "soft tracks that range out (stop) within the backward endcap KLM with B10 hit(s)", // outcome=32
238 "soft tracks that range out (stop) within the backward endcap KLM with B11 hit(s)", // outcome=33
239 "soft tracks that range out (stop) within the backward endcap KLM with B12 hit(s)", // outcome=34
240 "soft tracks that range out (stop) within the backward endcap KLM with B13 hit(s)", // outcome=35
241 "soft tracks that range out (stop) within the backward endcap KLM with B14 hit(s)", // outcome=36
242 "hard tracks that escape from the forward endcap KLM with B0 hit(s)", // outcome=37
243 "hard tracks that escape from the forward endcap KLM with B1 hit(s)", // outcome=38
244 "hard tracks that escape from the forward endcap KLM with B2 hit(s)", // outcome=39
245 "hard tracks that escape from the forward endcap KLM with B3 hit(s)", // outcome=40
246 "hard tracks that escape from the forward endcap KLM with B4 hit(s)", // outcome=41
247 "hard tracks that escape from the forward endcap KLM with B5 hit(s)", // outcome=42
248 "hard tracks that escape from the forward endcap KLM with B6 hit(s)", // outcome=43
249 "hard tracks that escape from the forward endcap KLM with B7 hit(s)", // outcome=44
250 "hard tracks that escape from the forward endcap KLM with B8 hit(s)", // outcome=45
251 "hard tracks that escape from the forward endcap KLM with B9 hit(s)", // outcome=46
252 "hard tracks that escape from the forward endcap KLM with B10 hit(s)", // outcome=47
253 "hard tracks that escape from the forward endcap KLM with B11 hit(s)", // outcome=48
254 "hard tracks that escape from the forward endcap KLM with B12 hit(s)", // outcome=49
255 "hard tracks that escape from the forward endcap KLM with B13 hit(s)", // outcome=50
256 "hard tracks that escape from the forward endcap KLM with B14 hit(s)", // outcome=51
257 "hard tracks that escape from the backward endcap KLM with B0 hit(s)", // outcome=52
258 "hard tracks that escape from the backward endcap KLM with B1 hit(s)", // outcome=53
259 "hard tracks that escape from the backward endcap KLM with B2 hit(s)", // outcome=54
260 "hard tracks that escape from the backward endcap KLM with B3 hit(s)", // outcome=55
261 "hard tracks that escape from the backward endcap KLM with B4 hit(s)", // outcome=56
262 "hard tracks that escape from the backward endcap KLM with B5 hit(s)", // outcome=57
263 "hard tracks that escape from the backward endcap KLM with B6 hit(s)", // outcome=58
264 "hard tracks that escape from the backward endcap KLM with B7 hit(s)", // outcome=59
265 "hard tracks that escape from the backward endcap KLM with B8 hit(s)", // outcome=60
266 "hard tracks that escape from the backward endcap KLM with B9 hit(s)", // outcome=61
267 "hard tracks that escape from the backward endcap KLM with B10 hit(s)", // outcome=62
268 "hard tracks that escape from the backward endcap KLM with B11 hit(s)", // outcome=63
269 "hard tracks that escape from the backward endcap KLM with B12 hit(s)", // outcome=64
270 "hard tracks that escape from the backward endcap KLM with B13 hit(s)", // outcome=65
271 "hard tracks that escape from the backward endcap KLM with B14 hit(s)"
272 }; // outcome=66
273 snprintf(label, sizeof(label), "Muon ID PDFs for %s", xmllabel);
274
275 //Histograms definition
276 TH1D* layerHitU[15][NOUTCOME];
277 TH1D* layerHitV[15][NOUTCOME];
278 TH1D* layerHitW[15][NOUTCOME];
279 TH1D* rchisqA[21][3];
280 TH1D* rchisqB[21][3];
281 TH1D* discardA;
282 TH1D* discardB;
283
284 //prepare histograms: title, bins and intervals
285 //Sumw2 -> store
sum
of
squares
of
weights
286 discardA = new TH1D("discardA", "discard fraction", 70, -0.5, 69.5);
287 discardB = new TH1D("discardB", "discard denominator", 70, -0.5, 69.5);
288 for (int outcome = 1; outcome < NOUTCOME; ++outcome) {
289 for (int layer = 0; layer < 15; ++layer) {
290 snprintf(label, sizeof(label), "LayerHitU-layer%02d-outcome%d", layer, outcome);
291 snprintf(title, sizeof(title), "LayerHitU for %s: layer %02d (%s)", xmllabel, layer, outcomeName[outcome]);
292 layerHitU[layer][outcome] = new TH1D(label, title, 32, -0.5, 31.5);
293 layerHitU[layer][outcome]->Sumw2(true);
294 snprintf(label, sizeof(label), "LayerHitV-layer%02d-outcome%d", layer, outcome);
295 snprintf(title, sizeof(title), "LayerHitV for %s: layer %02d (%s)", xmllabel, layer, outcomeName[outcome]);
296 layerHitV[layer][outcome] = new TH1D(label, title, 32, -0.5, 31.5);
297 layerHitV[layer][outcome]->Sumw2(true);
298 snprintf(label, sizeof(label), "LayerHitW-layer%02d-outcome%d", layer, outcome);
299 snprintf(title, sizeof(title), "LayerHitW for %s: layer %02d (%s)", xmllabel, layer, outcomeName[outcome]);
300 layerHitW[layer][outcome] = new TH1D(label, title, 32, -0.5, 31.5);
301 layerHitW[layer][outcome]->Sumw2(true);
302 }
303 }
304
305 for (int detector = 0; detector <= 2; ++detector) {
306 for (int halfNdof = 0; halfNdof <= 20; ++halfNdof) {
307 snprintf(label, sizeof(label), "ReducedChiSquaredA-%02ddof-detector%d", 2 * halfNdof, detector);
308 snprintf(title, sizeof(title), "Reduced Chi-squared for %s: %02d dof (%s)", xmllabel, 2 * halfNdof, detectorName[detector]);
309 rchisqA[halfNdof][detector] = new TH1D(label, title, 100, 0.0, 10.0);
310 rchisqA[halfNdof][detector]->Sumw2(true);
311 snprintf(label, sizeof(label), "ReducedChiSquaredB-%02ddof-detector%d", 2 * halfNdof, detector);
312 snprintf(title, sizeof(title), "Reduced Chi-squared for %s: %02d dof (%s)", xmllabel, 2 * halfNdof, detectorName[detector]);
313 rchisqB[halfNdof][detector] = new TH1D(label, title, 1000, 0.0, 100.0); // 10x wider range for integration
314 rchisqB[halfNdof][detector]->Sumw2(true);
315 }
316 }
317
318 printf("Finished defining histograms\n");
319
320 TObjArray* fileElements = chain->GetListOfFiles();
321 TIter next(fileElements);
322 TChainElement* chEl = 0;
323 while ((chEl = static_cast<TChainElement*>(next()))) {
324 //printf("Processing file %s ...\n", chEl->GetTitle());
325 TFile file(chEl->GetTitle());
326 TTree* tree = static_cast<TTree*>(file.Get("tree"));
327 tree->SetBranchStatus("*", 0);
328 tree->SetBranchStatus("MCParticles*", 1);
329 tree->SetBranchStatus("KLMMuidLikelihoods*", 1);
330 TClonesArray* mcParticles = 0;
331 tree->SetBranchAddress("MCParticles", &mcParticles);
332 TClonesArray* muids = 0;
333 tree->SetBranchAddress("KLMMuidLikelihoods", &muids);
334
335 printf("Number of events: %lld in file %s\n", tree->GetEntriesFast(), chEl->GetTitle());
336 ROOT::Math::PxPyPzEVector decayMomentum(0.0, 0.0, 0.0, 0.0);
337 //ROOT::Math::PxPyPzEVector temp(0.0, 0.0, 0.0, 0.0);
338 for (int event = 0; event < tree->GetEntriesFast(); ++event) {
339 tree->GetEntry(event);
340 if (mcParticles->GetEntriesFast() == 0) continue;
341 Belle2::MCParticle* mcp = static_cast<Belle2::MCParticle*>((*mcParticles)[0]);
342 Particle* particle = new Particle(mcp);
343 // look for muon decay-in-flight or hard radiation/delta production by examining daughters
344 if (isMuon) {
345 if (mcp->getNDaughters() > 0) {
346 decayMomentum.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
347 int firstDaughter = mcp->getFirstDaughter();
348 int lastDaughter = mcp->getLastDaughter();
349 for (int d = firstDaughter - 1; d < lastDaughter; ++d) {
350 const Belle2::MCParticle* daughter = static_cast<Belle2::MCParticle*>((*mcParticles)[d]);
351 int dPDG = abs(daughter->getPDG());
352 if ((dPDG != 11) && (dPDG != 12) && (dPDG != 14) && (dPDG != 22)) continue; // electron or neutrino or gamma?
353 //temp.SetVectM(daughter->getMomentum()*1000.0, (dPDG == 11 ? Const::electron.getMass() : 0.0));
354 //decayMomentum += temp;
355 ROOT::Math::XYZVector p3 = daughter->getMomentum();
356 decayMomentum.SetPx(decayMomentum.Px() + p3.X());
357 decayMomentum.SetPy(decayMomentum.Py() + p3.Y());
358 decayMomentum.SetPz(decayMomentum.Pz() + p3.Z());
359 }
360 if (decayMomentum.P() > 0.05) {
361 //debug printf("Event %5d: p(mu) = %f\n", event, decayMomentum.P());
362 continue;
363 }
364 }
365 }
366 // end decay-in-flight check
367
368 //Check acceptances and outcomes
369 double mcMomentumZ = 1000.0 * mcp->getMomentum().Z(); // MeV/c
370 //double mcMomentumR = 1000.0 * mcp->getMomentum().Rho(); // MeV/c
371 double mcTheta = mcp->getMomentum().Theta() * 180.0 / TMath::Pi(); // degrees
372 bool isForward = (mcMomentumZ > 0.0);
373 // if ((mcTheta > 31.5) && (mcTheta < 33.1)) continue; // avoid forward ECL crack
374 // if ((mcTheta > 126.9) && (mcTheta < 128.4)) continue; // avoid backward ECL crack
375 if (Variable::thetaInECLAcceptance(particle) == 0) continue; // avoid forward & backward ECL crack
376 double mcPhi = mcp->getMomentum().Phi() * 180.0 / TMath::Pi(); // degrees
377 if ((mcTheta > 110.0) && (mcTheta < 125.0) && (mcPhi > 60.0) && (mcPhi < 120.0)) continue; // avoid chimney
378 for (int m = 0; m < muids->GetEntriesFast(); ++m) {
379 const Belle2::KLMMuidLikelihood* muid = (Belle2::KLMMuidLikelihood*)((*muids)[m]);
380 // if (fabs(mcMomentumZ - muid->getMomentum().Z()) > 60.0) continue;
381 // if (fabs(mcMomentumR - muid->getMomentum().Perp()) > 60.0) continue;
382 int outcome = muid->getOutcome();
383 if (outcome <= 0) continue;
384 unsigned int extLayerPattern = muid->getExtLayerPattern();
385 unsigned int hitLayerPattern = muid->getHitLayerPattern();
386 int blayer = muid->getBarrelExtLayer();
387 int elayer = muid->getEndcapExtLayer();
388 int layer = blayer;
389 if (outcome == MuidElementNumbers::c_StopInForwardEndcap) { // forward endcap stop (no barrel hits). outcome=2
390 layer = elayer;
391 if (blayer < 0) {
392 outcome = isForward ? MuidElementNumbers::c_StopInForwardEndcap :
393 MuidElementNumbers::c_StopInBackwardEndcap; // backward endcap stop (no barrel hits)
394 } else {
395 outcome = (isForward ? MuidElementNumbers::c_CrossBarrelStopInForwardMin : MuidElementNumbers::c_CrossBarrelStopInBackwardMin) +
396 blayer; // forward endcap stop (B+E)
397 }
398 } else if (outcome == MuidElementNumbers::c_ExitForwardEndcap) { // forward endcap exit (no barrel hits)
399 layer = elayer;
400 if (blayer < 0) {
401 outcome = isForward ? MuidElementNumbers::c_ExitForwardEndcap :
402 MuidElementNumbers::c_ExitBackwardEndcap; // backward endcap exit (no barrel hits)
403 } else {
404 outcome = (isForward ? MuidElementNumbers::c_CrossBarrelExitForwardMin : MuidElementNumbers::c_CrossBarrelExitBackwardMin) +
405 blayer; // forward endcap exit (B+E)
406 }
407 }
408 int ndof = muid->getDegreesOfFreedom();
409 double rchisq = (ndof > 0 ? muid->getChiSquared() / ndof : -1.0);
410 if (ndof >= 40) continue;
411 int detector = 0; // crossed barrel and endcap
412 if (elayer < 0) {
413 detector = 1; // crossed barrel only
414 } else if (blayer < 0) {
415 detector = 2; // crossed endcap only
416 }
417 rchisqA[ndof / 2][detector]->Fill(rchisq);
418 unsigned int testbit = 1;
419 unsigned int mask = 0;
420 for (int bit = 0; bit <= blayer; ++bit) {
421 mask |= testbit;
422 testbit <<= 1;
423 }
424 // testbit = 1 << 15;
425 testbit = 1 << (MuidElementNumbers::getMaximalBarrelLayer() + 1);
426 for (int bit = 0; bit <= elayer; ++bit) {
427 mask |= testbit;
428 testbit <<= 1;
429 }
430 discardB->Fill(outcome);
431 if (extLayerPattern == mask) {
432 testbit = 1;
433 for (int bit = 0; bit < 32; ++bit) {
434 if ((testbit & extLayerPattern) != 0) {
435 double xLyr = (double)(bit);
436 layerHitV[layer][outcome]->Fill(xLyr);
437 if ((testbit & hitLayerPattern) != 0) {
438 layerHitU[layer][outcome]->Fill(xLyr);
439 layerHitW[layer][outcome]->Fill(xLyr);
440 }
441 }
442 testbit <<= 1;
443 }
444 } else {
445 discardA->Fill(outcome);
446 }
447 }
448 delete particle;
449 }
450
451 } // end of loop over files in the chain
452
453 TF1* rchisqFunc = new TF1("rchisqFunc", "[0]*[2]*std::pow(x*[1]*[2],[2]-1)/TMath::Gamma([2])*std::exp(-x*[1]*[2])", 0.0, 10.0);
454 TF1* layerHitFunc = new TF1("layerHitFunc", "[0]*std::exp(-[1]*(x-[2]))", -0.5, 28.5);
455 fprintf(output, " <%s>\n", xmllabel);
456 for (int outcome = 1; outcome <= 6; ++outcome) {
457 for (int layer = 0; layer < MuidElementNumbers::getMaximalBarrelLayer() + 1; ++layer) {
458 if ((outcome == MuidElementNumbers::c_StopInBarrel) && (layer > MuidElementNumbers::getMaximalBarrelLayer() - 1)) continue;
459 if ((outcome == MuidElementNumbers::c_StopInForwardEndcap)
460 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue;
461 if ((outcome == MuidElementNumbers::c_ExitBarrel) && (layer > MuidElementNumbers::getMaximalBarrelLayer())) continue;
462 if ((outcome == MuidElementNumbers::c_ExitForwardEndcap) && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer())) continue;
463 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap)
464 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue;
465 if ((outcome == MuidElementNumbers::c_ExitBackwardEndcap)
466 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue;
467 int bin0 = ((outcome == MuidElementNumbers::c_StopInBarrel)
468 || (outcome == MuidElementNumbers::c_ExitBarrel)) ? 0 : MuidElementNumbers::getMaximalBarrelLayer() + 1;
469 if (layerHitV[layer][outcome]->GetEntries() > 0.0) {
470 layerHitU[layer][outcome]->Divide(layerHitV[layer][outcome]);
471 }
472 printf("Fitting spectrum for outcome %d layer %d bin0 %d ...\n", outcome, layer, bin0);
473 double p0 = 1.0;
474 double p0Min = 0.01;
475 double p0Max = 2.0;
476 double p1 = 0.0;
477 double p1Min = 0.1;
478 double p1Max = 1.0;
479 double xMin = bin0 + 0.5;
480 if ((outcome == MuidElementNumbers::c_StopInBarrel) || (outcome == MuidElementNumbers::c_ExitBarrel)) xMin = bin0 + 1.5;
481 double xMax = bin0 - 0.5 + layer;
482 if (isMuon) {
483 p0 = 0.97; p0Min = 0.9; p0Max = 1.0; p1 = 0.005; p1Min = 0.0; p1Max = 0.02;
484 } else if (isPion) {
485 p0 = 0.29; p0Min = 0.02; p0Max = 0.5; p1 = 0.30; p1Min = 0.2; p1Max = 0.5;
486 } else if (isKaon) {
487 p0 = 0.33; p0Min = 0.02; p0Max = 0.5; p1 = 0.23; p1Min = 0.17; p1Max = 0.6;
488 } else if (isProt) {
489 p0 = 0.05; p0Min = 0.001; p0Max = 0.4; p1 = 0.36; p1Min = 0.3; p1Max = 0.7;
490 if ((outcome == MuidElementNumbers::c_StopInBarrel) && (layer >= 12)) xMin = bin0 + 3.5;
491 if ((outcome == MuidElementNumbers::c_StopInForwardEndcap) && (layer >= 8)) xMin = bin0 + 1.5;
492 if ((outcome == MuidElementNumbers::c_StopInForwardEndcap) && (layer >= 10)) xMin = bin0 + 2.5;
493 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap) && (layer >= 8)) xMin = bin0 + 1.5;
494 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap) && (layer >= 10)) xMin = bin0 + 2.5;
495 } else if (isDeut) {
496 p0 = 0.20; p0Min = 0.0002; p0Max = 0.26; p1 = 0.45; p1Min = 0.40; p1Max = 0.50;
497 if ((outcome == MuidElementNumbers::c_StopInBarrel) && (layer <= 4)) { p0 = 0.00076; p1 = 0.745; }
498 if ((outcome == MuidElementNumbers::c_StopInForwardEndcap) && (layer <= 5)) { p0 = 0.00027; p1 = 0.202; }
499 if ((outcome == MuidElementNumbers::c_ExitForwardEndcap) && (layer <= 4)) { p0 = 0.03; p1 = 0.83; }
500 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap) && (layer <= 4)) { p0 = 0.00021; p1 = 0.26; }
501 } else if (isElec) {
502 p0 = 0.035; p0Min = 0.005; p0Max = 0.25; p1 = 1.65; p1Min = 1.5; p1Max = 2.5;
503 if ((outcome == MuidElementNumbers::c_ExitForwardEndcap) && (layer >= 11)) xMin = bin0 + 1.5;
504 if ((outcome == MuidElementNumbers::c_ExitForwardEndcap) && (layer <= 10)) { p0 = 0.012; p1 = 1.186; }
505 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap) && (layer <= 10)) { p0 = 0.010; p1 = 2.33; }
506 if ((outcome == MuidElementNumbers::c_ExitBackwardEndcap) && (layer <= 10)) { p0 = 0.008; p1 = 1.94; }
507 }
508 double e0 = 0.0;
509 double e1 = 0.0;
510 printf(",B,%d,%d,%g,%g,%g,%g\n", outcome, layer, p0, 0.0, p1, 0.0);
511 layerHitFunc->SetParameters(p0, p1, (double)bin0);
512 layerHitFunc->SetParLimits(0, p0Min, p0Max);
513 layerHitFunc->SetParLimits(1, p1Min, p1Max);
514 layerHitFunc->FixParameter(2, (double)bin0);
515 layerHitFunc->ReleaseParameter(1);
516 if (int(xMax - xMin + 0.5) < layerHitFunc->GetNumberFreeParameters()) layerHitFunc->FixParameter(1, p1);
517 if (int(xMax - xMin + 0.5) >= layerHitFunc->GetNumberFreeParameters()) {
518 layerHitU[layer][outcome]->Fit("layerHitFunc", "WL", "", xMin, xMax);
519 p0 = layerHitFunc->GetParameter(0);
520 p1 = layerHitFunc->GetParameter(1);
521 e0 = layerHitFunc->GetParError(0);
522 e1 = layerHitFunc->GetParError(1);
523 printf(",A,%d,%d,%g,%g,%g,%g\n", outcome, layer, p0, e0, p1, e1);
524 }
525 for (int bin = bin0 + 1; bin <= bin0 + 15; ++bin) {
526 double y = layerHitU[layer][outcome]->GetBinContent(bin) * layerHitV[layer][outcome]->GetBinContent(bin);
527 double yFit = p0 * std::exp(-p1 * (bin - 1 - bin0)) * layerHitV[layer][outcome]->GetBinContent(bin);
528 if ((TMath::Sq(y - yFit) < 9.0 * (yFit + TMath::Max(1.0, y))) || (y <= 2.0)) {
529 layerHitW[layer][outcome]->SetBinContent(bin, yFit);
530 layerHitW[layer][outcome]->SetBinError(bin, yFit * TMath::Sqrt(TMath::Sq(e0 / p1) + TMath::Sq((bin - 1 - bin0)*e1)));
531 }
532 }
533 }
534 }
535 for (int outcome = MuidElementNumbers::c_CrossBarrelStopInForwardMin; outcome < NOUTCOME; ++outcome) {
536 for (int layer = 0; layer < MuidElementNumbers::getMaximalBarrelLayer() + 1; ++layer) {
537 // the first condition is implied by the loop bound but kept for symmetry with the checks below
538 // cppcheck-suppress knownConditionTrueFalse
539 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInForwardMin)
540 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInForwardMax)
541 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue; // like outcome == 2
542 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInBackwardMin)
543 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInBackwardMax)
544 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue; // like outcome == 5
545 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitForwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitForwardMax)
546 && (layer > MuidElementNumbers::getMaximalBarrelLayer() - 1)) continue; // like outcome == 4
547 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitBackwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitBackwardMax)
548 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue; // like outcome == 6
549 // Smooth out layerHitW, avoid zeroes
550 int elayer = layer;
551 int blayer = outcome - 7;
552 int outcomeRef = 2;
553 if (outcome >= MuidElementNumbers::c_CrossBarrelStopInBackwardMin) { outcomeRef = MuidElementNumbers::c_StopInBackwardEndcap; blayer = outcome - MuidElementNumbers::c_CrossBarrelStopInBackwardMin; }
554 if (outcome >= MuidElementNumbers::c_CrossBarrelExitForwardMin) { outcomeRef = MuidElementNumbers::c_ExitForwardEndcap; blayer = outcome - MuidElementNumbers::c_CrossBarrelExitForwardMin; }
555 if (outcome >= MuidElementNumbers::c_CrossBarrelExitBackwardMin) { outcomeRef = MuidElementNumbers::c_ExitBackwardEndcap; blayer = outcome - MuidElementNumbers::c_CrossBarrelExitBackwardMin; }
556 double wSum = 0.0;
557 double vSum = 0.0;
558 double wSumRef = 0.0;
559 double vSumRef = 0.0;
560 for (int bin = 1; bin <= blayer + 1; ++bin) {
561 wSum += layerHitW[layer][outcome]->GetBinContent(bin);
562 vSum += layerHitV[layer][outcome]->GetBinContent(bin);
563 // blayer is reduced to [0, getMaximalBarrelLayer()] by the assignments above,
564 // but cppcheck evaluates the intermediate 'outcome - 7' value
565 // cppcheck-suppress arrayIndexOutOfBoundsCond
566 wSumRef += layerHitW[blayer][3]->GetBinContent(bin);
567 // cppcheck-suppress arrayIndexOutOfBoundsCond
568 vSumRef += layerHitV[blayer][3]->GetBinContent(bin);
569 }
570 if (wSum == 0.0) {
571 wSum = 0.2;
572 if (vSum == 0.0) vSum = wSum;
573 }
574 double factor = (wSum / vSum) / (wSumRef / vSumRef);
575 for (int bin = 1; bin <= blayer + 1; ++bin) {
576 double w = layerHitW[layer][outcome]->GetBinContent(bin);
577 double v = layerHitV[layer][outcome]->GetBinContent(bin);
578 double v1 = TMath::Max(v, 1.0);
579 double wRef = layerHitW[blayer][3]->GetBinContent(bin);
580 double vRef = layerHitV[blayer][3]->GetBinContent(bin);
581 double wRef1 = TMath::Max(wRef, 1.0);
582 double wNew = wRef * factor * (v1 / vRef);
583 wNew = TMath::Min((isMuon ? 0.95 : 0.6667) * v1, wNew);
584 double e = (w == 0.0) ? 1.0 : w;
585 double eNew = factor * factor * wRef1 * v1 * (v1 * (vRef - wRef) + wRef * vRef) / (vRef * vRef * vRef);
586 wNew = (w / e + wNew / eNew) / (1.0 / e + 1.0 / eNew);
587 if (!isMuon && (bin > 1)) {
588 wNew = TMath::Min(wNew, layerHitW[layer][outcome]->GetBinContent(bin - 1));
589 }
590 //if ((TMath::Sq(w - wNew) < 9.0 * (e + eNew)) || (w <= 2.0)) {
591 layerHitW[layer][outcome]->SetBinContent(bin, wNew);
592 //}
593 printf("%d,%d,%d,%d,%g,%g,%g,%g,%g,%g,%g\n", outcome, layer, blayer, bin, factor, w, v, wRef, vRef, wNew,
594 TMath::Sq(w - wNew) / (e + eNew));
595 }
596 wSum = 0.0;
597 vSum = 0.0;
598 wSumRef = 0.0;
599 vSumRef = 0.0;
600 for (int bin = 1; bin <= elayer + 1; ++bin) {
601 wSum += layerHitW[layer][outcome]->GetBinContent(bin + 15);
602 vSum += layerHitV[layer][outcome]->GetBinContent(bin + 15);
603 wSumRef += layerHitW[layer][outcomeRef]->GetBinContent(bin + 15);
604 vSumRef += layerHitV[layer][outcomeRef]->GetBinContent(bin + 15);
605 }
606 if (wSum <= 1.0) {
607 factor *= (layerHitW[blayer][3]->GetBinContent(blayer + 1) / layerHitV[blayer][3]->GetBinContent(blayer + 1))
608 / (layerHitW[layer][outcomeRef]->GetBinContent(1 + 15) / layerHitV[layer][outcomeRef]->GetBinContent(1 + 15));
609 } else {
610 factor = (wSum / vSum) / (wSumRef / vSumRef);
611 }
612 for (int bin = 1; bin <= elayer + 1; ++bin) {
613 double w = layerHitW[layer][outcome]->GetBinContent(bin + 15);
614 double v = layerHitV[layer][outcome]->GetBinContent(bin + 15);
615 double v1 = TMath::Max(v, 1.0);
616 double wRef = layerHitW[layer][outcomeRef]->GetBinContent(bin + 15);
617 double vRef = layerHitV[layer][outcomeRef]->GetBinContent(bin + 15);
618 double wRef1 = TMath::Max(wRef, 1.0);
619 double wNew = wRef * factor * (v1 / vRef);
620 wNew = TMath::Min((isMuon ? 0.95 : 0.6667) * v1, wNew);
621 double e = (w == 0.0) ? 1.0 : w;
622 double eNew = factor * factor * wRef1 * v1 * (v1 * (vRef - wRef) + wRef * vRef) / (vRef * vRef * vRef);
623 wNew = (w / e + wNew / eNew) / (1.0 / e + 1.0 / eNew);
624 if (!isMuon && (bin > 1)) {
625 wNew = TMath::Min(wNew, layerHitW[layer][outcome]->GetBinContent(bin + 15 - 1));
626 }
627 //if ((TMath::Sq(w - wNew) < 9.0 * (e + eNew)) || (w <= 2.0)) {
628 layerHitW[layer][outcome]->SetBinContent(bin + 15, wNew);
629 //}
630 printf("%d,%d,%d,%d,%g,%g,%g,%g,%g,%g,%g\n", outcome, layer, blayer, bin, factor, w, v, wRef, vRef, wNew,
631 TMath::Sq(w - wNew) / (e + eNew));
632 }
633 }
634 }
635 fprintf(output, " <LayerProfile>\n");
636 for (int outcome = 1; outcome < NOUTCOME; ++outcome) {
637 fprintf(output, " <!-- %s -->\n", outcomeComment[outcome]);
638 fprintf(output, " <Outcome outcome=\"%d\">\n", outcome);
639 for (int layer = 0; layer < MuidElementNumbers::getMaximalBarrelLayer() + 1; ++layer) {
640 if ((outcome == MuidElementNumbers::c_StopInBarrel) && (layer > MuidElementNumbers::getMaximalBarrelLayer() - 1)) continue;
641 if ((outcome == MuidElementNumbers::c_StopInForwardEndcap)
642 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue;
643 if ((outcome == MuidElementNumbers::c_ExitBarrel) && (layer > MuidElementNumbers::getMaximalBarrelLayer())) continue;
644 if ((outcome == MuidElementNumbers::c_ExitForwardEndcap) && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer())) continue;
645 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap)
646 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue;
647 if ((outcome == MuidElementNumbers::c_ExitBackwardEndcap)
648 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue;
649 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInForwardMin)
650 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInForwardMax)
651 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue; // like outcome == 2
652 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInBackwardMin)
653 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInBackwardMax)
654 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue; // like outcome == 5
655 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitForwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitForwardMax)
656 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer())) continue; // like outcome == 4
657 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitBackwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitBackwardMax)
658 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue; // like outcome == 6
659 int lastBin = layer + 1;
660 if ((outcome != MuidElementNumbers::c_StopInBarrel)
661 && (outcome != MuidElementNumbers::c_ExitBarrel)) lastBin += MuidElementNumbers::getMaximalBarrelLayer() + 1;
662 // the enclosing loop bound is 'outcome < NOUTCOME', so outcome stays in range
663 // cppcheck-suppress arrayIndexOutOfBoundsCond
664 if (layerHitV[layer][outcome]->GetEntries() > 0.0) {
665 if (outcome >= MuidElementNumbers::c_CrossBarrelStopInForwardMin) layerHitU[layer][outcome]->Divide(layerHitV[layer][outcome]);
666 layerHitW[layer][outcome]->Divide(layerHitV[layer][outcome]);
667 }
668 fprintf(output, " <LastLayer layer=\"%d\">\n ", layer);
669 for (int bin = 1; bin <= lastBin; ++bin) {
670 fprintf(output, "%16.8E", layerHitW[layer][outcome]->GetBinContent(bin));
671 if ((bin % 8 == 0) && (bin != lastBin)) fprintf(output, "\n ");
672 }
673 fprintf(output, "\n </LastLayer>\n");
674 }
675 fprintf(output, " </Outcome>\n");
676 }
677 fprintf(output, " </LayerProfile>\n");
678
679 // --------------------------------------------------------------------------
680
681 int nBins = rchisqA[0][1]->GetNbinsX();
682 double binWidth = rchisqA[0][1]->GetBinWidth(1);
683 double lowEdge = rchisqA[0][1]->GetBinLowEdge(1);
684 fprintf(output, " <TransversePDF>\n");
685 for (int detector = 0; detector <= 2; ++detector) {
686 fprintf(output, " <%s>\n", detectorName[detector]);
687 double p1 = 1.0; // successful best-fit value will be inherited for next iteration of the ndof loop
688 for (int halfNdof = 1; halfNdof <= 18; ++halfNdof) {
689 //double p0Err = 0.0;
690 //double p1Err = 0.0;
691 fprintf(output, " <DegreesOfFreedom ndof=\"%d\">\n", 2 * halfNdof);
692 double entries = rchisqA[halfNdof][detector]->GetEntries();
693 double p0 = entries * binWidth;
694 if (entries == 0) {
695 //int bin = (int)((1.0 - lowEdge) / binWidth + 0.5) + 1;
696 p0 = 1.0 * binWidth;
697 // DIVOT p1 = 1.0;
698 }
699 double p2 = (double)(halfNdof);
700 rchisqFunc->SetParameters(p0, p1, p2);
701 rchisqFunc->ReleaseParameter(0);
702 rchisqFunc->ReleaseParameter(1);
703 rchisqFunc->FixParameter(2, p2);
704 int peakBin = rchisqA[halfNdof][detector]->GetMaximumBin();
705 double xMin = 0.0;
706 double xMax = lowEdge + binWidth * (peakBin + 10);
707 double p0Left = p0;
708 double p1Left = p1;
709 if ((halfNdof > 1) && (entries > 0)) {
710 rchisqA[halfNdof][detector]->Fit("rchisqFunc", "WL", "", xMin, xMax);
711 p0Left = rchisqFunc->GetParameter(0);
712 p1Left = rchisqFunc->GetParameter(1);
713 }
714 if (halfNdof == 2) {
715 if (!isMuon) p0 = entries * (isElec ? 1.0 : 10.0);
716 p1 = 1.4;
717 }
718 xMax = 10.0;
719 if ((entries > 500.0) && (halfNdof <= 15)) {
720 double peak = rchisqA[halfNdof][detector]->GetBinContent(peakBin);
721 //double peakX = lowEdge + binWidth * peakBin;
722 int binLow = peakBin;
723 double threshold = ((isMuon || isPion || isKaon) ? 0.1 : 0.2) * peak;
724 if (halfNdof >= 10) { threshold *= 2.0; }
725 for (int bin = peakBin; bin <= nBins; ++bin) {
726 if (rchisqA[halfNdof][detector]->GetBinContent(bin) < threshold) break;
727 binLow = bin - 2;
728 if (binLow < 1) { binLow = 1; }
729 }
730 xMin = lowEdge + binWidth * binLow;
731 /*
732 if (isMuon) {
733 if (halfNdof == 1) { xMin = 6.4; }
734 } else {
735 if (halfNdof == 1) { xMin = 8.0; }
736 if (halfNdof == 2) { xMin = (isElec ? 3.8 : 5.4); }
737 if (halfNdof == 3) { xMin = (isElec ? 2.4 : 4.0); }
738 if (halfNdof == 4) { xMin = (isElec ? 0.0 : 3.4); }
739 }
740 */
741 for (int bin = rchisqA[halfNdof][detector]->GetNbinsX(); bin > 0; --bin) {
742 if (rchisqA[halfNdof][detector]->GetBinContent(bin) > 0.0) {
743 int bin2 = bin - (33 - 3 * halfNdof);
744 if (bin2 > bin - 2) { bin2 = bin - 2; }
745 if (bin2 < peakBin) {
746 xMin = 0.0;
747 } else {
748 double xMin2 = lowEdge + binWidth * bin2;
749 if (xMin2 < xMin) {
750 xMin = xMin2;
751 }
752 }
753 break;
754 }
755 }
756 if (isMuon) {
757 if (halfNdof == 1) { xMin = 6.4; }
758 } else {
759 if (halfNdof == 1) { xMin = 8.0; }
760 if (halfNdof == 2) { xMin = (isElec ? 3.8 : 4.6); }
761 if (halfNdof == 3) { xMin = (isElec ? 2.4 : 4.0); }
762 if (halfNdof == 4) { xMin = (isElec ? 0.0 : 3.4); }
763 }
764 }
765 /*
766 if (halfNdof == 1) {
767 rchisqFunc->SetParLimits(1, 0.25, 3.0);
768 } else {
769 rchisqFunc->SetParLimits(1, 0.25, 2.0);
770 }
771 */
772 if (entries < 50.0) {
773 rchisqFunc->FixParameter(1, p1);
774 }
775 if ((entries > 0) && (int((xMax - xMin) / binWidth + 0.5) >= rchisqFunc->GetNumberFreeParameters())) {
776 rchisqA[halfNdof][detector]->Fit("rchisqFunc", "WL", "", xMin, xMax);
777 p0 = rchisqFunc->GetParameter(0);
778 //p0Err = rchisqFunc->GetParError(0);
779 p1 = rchisqFunc->GetParameter(1);
780 //p1Err = rchisqFunc->GetParError(1);
781 }
782
783 double yNorm = 1.0;
784 if ((entries > 0.0) && (rchisqA[halfNdof][detector]->GetEntries() > 0.0)) {
785 yNorm = rchisqA[halfNdof][detector]->GetEntries() / entries;
786 }
787 double yScaleLeft = yNorm * p0Left * p2 / TMath::Gamma(p2);
788 double yScale = yNorm * p0 * p2 / TMath::Gamma(p2);
789 for (int bin = 1; bin <= nBins; ++bin) {
790 double x = lowEdge + (bin - 0.5) * binWidth;
791 double y = rchisqA[halfNdof][detector]->GetBinContent(bin) * yNorm;
792 if ((rchisqA[halfNdof][detector]->GetBinContent(bin) == 0.0) || (x > xMin)) {
793 if ((bin < peakBin) && (xMin > 0.0)) {
794 x *= p1Left * p2;
795 y = yScaleLeft * std::pow(x, p2 - 1.0) * std::exp(-x);
796 } else {
797 x *= p1 * p2;
798 y = yScale * std::pow(x, p2 - 1.0) * std::exp(-x);
799 }
800 }
801 rchisqB[halfNdof][detector]->SetBinContent(bin, y);
802 }
803 for (int bin = nBins + 1; bin <= 100; ++bin) {
804 double x = (lowEdge + (bin - 0.5) * binWidth) * p1 * p2;
805 double y = yScale * std::pow(x, p2 - 1.0) * std::exp(-x);
806 rchisqB[halfNdof][detector]->SetBinContent(bin, y);
807 }
808 double sum = 0.0;
809 for (int bin = 100; bin > 0; --bin) {
810 sum += rchisqB[halfNdof][detector]->GetBinContent(bin);
811 }
812 sum *= binWidth;
813 fprintf(output, " <Tail>\n");
814 fprintf(output, " <Threshold>%16.7E</Threshold>\n", xMin);
815 fprintf(output, " <ScaleY>%16.7E</ScaleY>\n", yScale / sum);
816 fprintf(output, " <ScaleX>%16.7E</ScaleX>\n", p1);
817 fprintf(output, " </Tail>\n");
818 fprintf(output, " <Histogram>\n");
819 for (int bin = 1; bin <= nBins; bin += 10) {
820 fprintf(output, " %16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E\n",
821 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin) / sum),
822 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 1) / sum),
823 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 2) / sum),
824 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 3) / sum),
825 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 4) / sum),
826 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 5) / sum),
827 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 6) / sum),
828 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 7) / sum),
829 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 8) / sum),
830 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 9) / sum));
831 }
832 fprintf(output, " </Histogram>\n");
833 fprintf(output, " </DegreesOfFreedom>\n");
834 }
835 fprintf(output, " </%s>\n", detectorName[detector]);
836 }
837 fprintf(output, " </TransversePDF>\n");
838 fprintf(output, " </%s>\n", xmllabel);
839
840 //---------------------------------------------------------------------
841
842 TCanvas* canvas = new TCanvas("tree", label, 900, 900);
843 //gStyle->SetPaperSize(1); // BAD BAD BAD!
844 gStyle->SetOptStat(10);
845 gStyle->SetOptFit(0);
846 canvas->GetPad(0)->SetLogy(0);
847 canvas->GetPad(0)->SetGrid(1, 1);
848 canvas->Clear();
849 canvas->Divide(1, 1);
850 canvas->GetPad(1)->SetLogy(1);
851 snprintf(label, sizeof(label), "Muid-%s.pdf[", xmllabel); // open file without drawing
852 canvas->SaveAs(label);
853 snprintf(label, sizeof(label), "Muid-%s.pdf", xmllabel); // save each canvas to same file
854 canvas->cd(1);
855 discardA->Divide(discardB);
856 discardA->SetLineColor(kRed);
857 discardA->Draw("");
858 canvas->SaveAs(label);
859 for (int outcome = 1; outcome < NOUTCOME; ++outcome) {
860 for (int layer = 0; layer < MuidElementNumbers::getMaximalBarrelLayer() + 1; ++layer) {
861 canvas->cd(1);
862 layerHitU[layer][outcome]->SetMaximum(2.0);
863 layerHitU[layer][outcome]->SetMinimum(1.0E-4);
864 if (isProt or isDeut) layerHitU[layer][outcome]->SetMinimum(1.0E-6);
865 if (isElec) layerHitU[layer][outcome]->SetMinimum(1.0E-10);
866 layerHitU[layer][outcome]->GetXaxis()->SetTitle("KLM Layer");
867 layerHitU[layer][outcome]->GetYaxis()->SetTitle("Frequency");
868 layerHitU[layer][outcome]->Draw("E1");
869 layerHitW[layer][outcome]->SetLineColor(kMagenta);
870 layerHitW[layer][outcome]->Draw("HIST,SAME");
871 layerHitU[layer][outcome]->Draw("E1,SAME");
872 canvas->SaveAs(label);
873 }
874 }
875
876 canvas->Clear();
877 canvas->Divide(1, 2);
878 canvas->GetPad(1)->SetLogy(1);
879 canvas->GetPad(2)->SetLogy(1);
880 gStyle->SetOptStat(110);
881 gStyle->SetOptFit(2);
882 for (int detector = 0; detector <= 2; ++detector) {
883 for (int halfNdof = 1; halfNdof <= 18; ++halfNdof) {
884 int padNo = (halfNdof - 1) % 2 + 1;
885 canvas->cd(padNo);
886 gPad->Clear();
887 rchisqA[halfNdof][detector]->SetMinimum(0.001); // for vertical log axis
888 if (rchisqA[halfNdof][detector]->GetEntries() == 0.0) { rchisqA[halfNdof][detector]->SetMaximum(2.0); }
889 rchisqA[halfNdof][detector]->GetXaxis()->SetTitle("#chi^{2}/ndf");
890 rchisqA[halfNdof][detector]->GetYaxis()->SetTitle("Frequency");
891 rchisqA[halfNdof][detector]->Draw("");
892 rchisqB[halfNdof][detector]->SetLineColor(kMagenta);
893 rchisqB[halfNdof][detector]->Draw("SAME");
894 rchisqA[halfNdof][detector]->Draw("SAME,E1");
895 if (padNo == 2) {
896 canvas->SaveAs(label);
897 }
898 }
899 }
900
901 snprintf(label, sizeof(label), "Muid-%s.pdf]", xmllabel); // close file without drawing
902 canvas->SaveAs(label);
903
904 printf("Done.\n");
905
906}
907
908int main(int argc, char** argv)
909{
910 TChain* chain;
911 FILE* output;
912 if (argc < 4)
913 goto badUsage;
914 chain = new TChain;
915 chain->Add(argv[1]);
916 output = std::fopen(argv[2], "a");
917 if (output == nullptr) {
918 std::printf("Cannot open output file %s\n", argv[2]);
919 delete chain;
920 return 1;
921 }
922 makeprob(chain, output, argv[3]);
923 std::fclose(output);
924 delete chain;
925 return 0;
926badUsage:
927 std::printf(
928 "Usage:\n"
929 "b2klm-likelihood-parameters input_file output_file xml_label\n");
930 return -1;
931}
Belle2::KLMMuidLikelihood
Class to store the likelihoods from KLM with additional information related to the extrapolation.
Definition KLMMuidLikelihood.h:28
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
Belle2::MCParticle::getNDaughters
int getNDaughters() const
Return number of daughter MCParticles.
Definition MCParticle.cc:73
Belle2::MCParticle::getLastDaughter
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
Definition MCParticle.h:249
Belle2::MCParticle::getMomentum
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition MCParticle.h:187
Belle2::MCParticle::getFirstDaughter
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
Definition MCParticle.h:241
Belle2::MuidElementNumbers::getMaximalEndcapForwardLayer
static constexpr int getMaximalEndcapForwardLayer()
Get maximal endcap-forward layer number (0-based).
Definition MuidElementNumbers.h:138
Belle2::MuidElementNumbers::getMaximalEndcapBackwardLayer
static constexpr int getMaximalEndcapBackwardLayer()
Get maximal endcap-forward layer number (0-based).
Definition MuidElementNumbers.h:146
Belle2::MuidElementNumbers::getMaximalBarrelLayer
static constexpr int getMaximalBarrelLayer()
Get maximal barrel layer number (0-based).
Definition MuidElementNumbers.h:130
Belle2::Particle
Class to store reconstructed particles.
Definition Particle.h:76
Belle2::ECLElementNumbers::isForward
bool isForward(int cellId)
Check whether the crystal is in forward ECL.
Definition ECLElementNumbers.h:41
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17