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 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, 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 char detectorName[3][40] = { "BarrelAndEndcap", "BarrelOnly", "EndcapOnly" };
136 //List of all the possible outcome and relative explanations
137 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 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 sprintf(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 sprintf(label, "LayerHitU-layer%02d-outcome%d", layer, outcome);
291 sprintf(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 sprintf(label, "LayerHitV-layer%02d-outcome%d", layer, outcome);
295 sprintf(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 sprintf(label, "LayerHitW-layer%02d-outcome%d", layer, outcome);
299 sprintf(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 sprintf(label, "ReducedChiSquaredA-%02ddof-detector%d", 2 * halfNdof, detector);
308 sprintf(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 sprintf(label, "ReducedChiSquaredB-%02ddof-detector%d", 2 * halfNdof, detector);
312 sprintf(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 = (TChainElement*)next())) {
324 //printf("Processing file %s ...\n", chEl->GetTitle());
325 TFile file(chEl->GetTitle());
326 TTree* tree = (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 = (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 Belle2::MCParticle* daughter = (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 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 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInForwardMin)
538 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInForwardMax)
539 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue; // like outcome == 2
540 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInBackwardMin)
541 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInBackwardMax)
542 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue; // like outcome == 5
543 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitForwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitForwardMax)
544 && (layer > MuidElementNumbers::getMaximalBarrelLayer() - 1)) continue; // like outcome == 4
545 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitBackwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitBackwardMax)
546 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue; // like outcome == 6
547 // Smooth out layerHitW, avoid zeroes
548 int elayer = layer;
549 int blayer = outcome - 7;
550 int outcomeRef = 2;
551 if (outcome >= MuidElementNumbers::c_CrossBarrelStopInBackwardMin) { outcomeRef = MuidElementNumbers::c_StopInBackwardEndcap; blayer = outcome - MuidElementNumbers::c_CrossBarrelStopInBackwardMin; }
552 if (outcome >= MuidElementNumbers::c_CrossBarrelExitForwardMin) { outcomeRef = MuidElementNumbers::c_ExitForwardEndcap; blayer = outcome - MuidElementNumbers::c_CrossBarrelExitForwardMin; }
553 if (outcome >= MuidElementNumbers::c_CrossBarrelExitBackwardMin) { outcomeRef = MuidElementNumbers::c_ExitBackwardEndcap; blayer = outcome - MuidElementNumbers::c_CrossBarrelExitBackwardMin; }
554 double wSum = 0.0;
555 double vSum = 0.0;
556 double wSumRef = 0.0;
557 double vSumRef = 0.0;
558 for (int bin = 1; bin <= blayer + 1; ++bin) {
559 wSum += layerHitW[layer][outcome]->GetBinContent(bin);
560 vSum += layerHitV[layer][outcome]->GetBinContent(bin);
561 wSumRef += layerHitW[blayer][3]->GetBinContent(bin);
562 vSumRef += layerHitV[blayer][3]->GetBinContent(bin);
563 }
564 if (wSum == 0.0) {
565 wSum = 0.2;
566 if (vSum == 0.0) vSum = wSum;
567 }
568 double factor = (wSum / vSum) / (wSumRef / vSumRef);
569 for (int bin = 1; bin <= blayer + 1; ++bin) {
570 double w = layerHitW[layer][outcome]->GetBinContent(bin);
571 double v = layerHitV[layer][outcome]->GetBinContent(bin);
572 double v1 = TMath::Max(v, 1.0);
573 double wRef = layerHitW[blayer][3]->GetBinContent(bin);
574 double vRef = layerHitV[blayer][3]->GetBinContent(bin);
575 double wRef1 = TMath::Max(wRef, 1.0);
576 double wNew = wRef * factor * (v1 / vRef);
577 wNew = TMath::Min((isMuon ? 0.95 : 0.6667) * v1, wNew);
578 double e = (w == 0.0) ? 1.0 : w;
579 double eNew = factor * factor * wRef1 * v1 * (v1 * (vRef - wRef) + wRef * vRef) / (vRef * vRef * vRef);
580 wNew = (w / e + wNew / eNew) / (1.0 / e + 1.0 / eNew);
581 if (!isMuon && (bin > 1)) {
582 wNew = TMath::Min(wNew, layerHitW[layer][outcome]->GetBinContent(bin - 1));
583 }
584 //if ((TMath::Sq(w - wNew) < 9.0 * (e + eNew)) || (w <= 2.0)) {
585 layerHitW[layer][outcome]->SetBinContent(bin, wNew);
586 //}
587 printf("%d,%d,%d,%d,%g,%g,%g,%g,%g,%g,%g\n", outcome, layer, blayer, bin, factor, w, v, wRef, vRef, wNew,
588 TMath::Sq(w - wNew) / (e + eNew));
589 }
590 wSum = 0.0;
591 vSum = 0.0;
592 wSumRef = 0.0;
593 vSumRef = 0.0;
594 for (int bin = 1; bin <= elayer + 1; ++bin) {
595 wSum += layerHitW[layer][outcome]->GetBinContent(bin + 15);
596 vSum += layerHitV[layer][outcome]->GetBinContent(bin + 15);
597 wSumRef += layerHitW[layer][outcomeRef]->GetBinContent(bin + 15);
598 vSumRef += layerHitV[layer][outcomeRef]->GetBinContent(bin + 15);
599 }
600 if (wSum <= 1.0) {
601 factor *= (layerHitW[blayer][3]->GetBinContent(blayer + 1) / layerHitV[blayer][3]->GetBinContent(blayer + 1))
602 / (layerHitW[layer][outcomeRef]->GetBinContent(1 + 15) / layerHitV[layer][outcomeRef]->GetBinContent(1 + 15));
603 } else {
604 factor = (wSum / vSum) / (wSumRef / vSumRef);
605 }
606 for (int bin = 1; bin <= elayer + 1; ++bin) {
607 double w = layerHitW[layer][outcome]->GetBinContent(bin + 15);
608 double v = layerHitV[layer][outcome]->GetBinContent(bin + 15);
609 double v1 = TMath::Max(v, 1.0);
610 double wRef = layerHitW[layer][outcomeRef]->GetBinContent(bin + 15);
611 double vRef = layerHitV[layer][outcomeRef]->GetBinContent(bin + 15);
612 double wRef1 = TMath::Max(wRef, 1.0);
613 double wNew = wRef * factor * (v1 / vRef);
614 wNew = TMath::Min((isMuon ? 0.95 : 0.6667) * v1, wNew);
615 double e = (w == 0.0) ? 1.0 : w;
616 double eNew = factor * factor * wRef1 * v1 * (v1 * (vRef - wRef) + wRef * vRef) / (vRef * vRef * vRef);
617 wNew = (w / e + wNew / eNew) / (1.0 / e + 1.0 / eNew);
618 if (!isMuon && (bin > 1)) {
619 wNew = TMath::Min(wNew, layerHitW[layer][outcome]->GetBinContent(bin + 15 - 1));
620 }
621 //if ((TMath::Sq(w - wNew) < 9.0 * (e + eNew)) || (w <= 2.0)) {
622 layerHitW[layer][outcome]->SetBinContent(bin + 15, wNew);
623 //}
624 printf("%d,%d,%d,%d,%g,%g,%g,%g,%g,%g,%g\n", outcome, layer, blayer, bin, factor, w, v, wRef, vRef, wNew,
625 TMath::Sq(w - wNew) / (e + eNew));
626 }
627 }
628 }
629 fprintf(output, " <LayerProfile>\n");
630 for (int outcome = 1; outcome < NOUTCOME; ++outcome) {
631 fprintf(output, " <!-- %s -->\n", outcomeComment[outcome]);
632 fprintf(output, " <Outcome outcome=\"%d\">\n", outcome);
633 for (int layer = 0; layer < MuidElementNumbers::getMaximalBarrelLayer() + 1; ++layer) {
634 if ((outcome == MuidElementNumbers::c_StopInBarrel) && (layer > MuidElementNumbers::getMaximalBarrelLayer() - 1)) continue;
635 if ((outcome == MuidElementNumbers::c_StopInForwardEndcap)
636 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue;
637 if ((outcome == MuidElementNumbers::c_ExitBarrel) && (layer > MuidElementNumbers::getMaximalBarrelLayer())) continue;
638 if ((outcome == MuidElementNumbers::c_ExitForwardEndcap) && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer())) continue;
639 if ((outcome == MuidElementNumbers::c_StopInBackwardEndcap)
640 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue;
641 if ((outcome == MuidElementNumbers::c_ExitBackwardEndcap)
642 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue;
643 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInForwardMin)
644 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInForwardMax)
645 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer() - 1)) continue; // like outcome == 2
646 if ((outcome >= MuidElementNumbers::c_CrossBarrelStopInBackwardMin)
647 && (outcome <= MuidElementNumbers::c_CrossBarrelStopInBackwardMax)
648 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer() - 1)) continue; // like outcome == 5
649 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitForwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitForwardMax)
650 && (layer > MuidElementNumbers::getMaximalEndcapForwardLayer())) continue; // like outcome == 4
651 if ((outcome >= MuidElementNumbers::c_CrossBarrelExitBackwardMin) && (outcome <= MuidElementNumbers::c_CrossBarrelExitBackwardMax)
652 && (layer > MuidElementNumbers::getMaximalEndcapBackwardLayer())) continue; // like outcome == 6
653 int lastBin = layer + 1;
654 if ((outcome != MuidElementNumbers::c_StopInBarrel)
655 && (outcome != MuidElementNumbers::c_ExitBarrel)) lastBin += MuidElementNumbers::getMaximalBarrelLayer() + 1;
656 if (layerHitV[layer][outcome]->GetEntries() > 0.0) {
657 if (outcome >= MuidElementNumbers::c_CrossBarrelStopInForwardMin) layerHitU[layer][outcome]->Divide(layerHitV[layer][outcome]);
658 layerHitW[layer][outcome]->Divide(layerHitV[layer][outcome]);
659 }
660 fprintf(output, " <LastLayer layer=\"%d\">\n ", layer);
661 for (int bin = 1; bin <= lastBin; ++bin) {
662 fprintf(output, "%16.8E", layerHitW[layer][outcome]->GetBinContent(bin));
663 if ((bin % 8 == 0) && (bin != lastBin)) fprintf(output, "\n ");
664 }
665 fprintf(output, "\n </LastLayer>\n");
666 }
667 fprintf(output, " </Outcome>\n");
668 }
669 fprintf(output, " </LayerProfile>\n");
670
671 // --------------------------------------------------------------------------
672
673 int nBins = rchisqA[0][1]->GetNbinsX();
674 double binWidth = rchisqA[0][1]->GetBinWidth(1);
675 double lowEdge = rchisqA[0][1]->GetBinLowEdge(1);
676 fprintf(output, " <TransversePDF>\n");
677 for (int detector = 0; detector <= 2; ++detector) {
678 fprintf(output, " <%s>\n", detectorName[detector]);
679 double p1 = 1.0; // successful best-fit value will be inherited for next iteration of the ndof loop
680 for (int halfNdof = 1; halfNdof <= 18; ++halfNdof) {
681 //double p0Err = 0.0;
682 //double p1Err = 0.0;
683 fprintf(output, " <DegreesOfFreedom ndof=\"%d\">\n", 2 * halfNdof);
684 double entries = rchisqA[halfNdof][detector]->GetEntries();
685 double p0 = entries * binWidth;
686 if (entries == 0) {
687 //int bin = (int)((1.0 - lowEdge) / binWidth + 0.5) + 1;
688 p0 = 1.0 * binWidth;
689 // DIVOT p1 = 1.0;
690 }
691 double p2 = (double)(halfNdof);
692 rchisqFunc->SetParameters(p0, p1, p2);
693 rchisqFunc->ReleaseParameter(0);
694 rchisqFunc->ReleaseParameter(1);
695 rchisqFunc->FixParameter(2, p2);
696 int peakBin = rchisqA[halfNdof][detector]->GetMaximumBin();
697 double xMin = 0.0;
698 double xMax = lowEdge + binWidth * (peakBin + 10);
699 double p0Left = p0;
700 double p1Left = p1;
701 if ((halfNdof > 1) && (entries > 0)) {
702 rchisqA[halfNdof][detector]->Fit("rchisqFunc", "WL", "", xMin, xMax);
703 p0Left = rchisqFunc->GetParameter(0);
704 p1Left = rchisqFunc->GetParameter(1);
705 }
706 if (halfNdof == 2) {
707 if (!isMuon) p0 = entries * (isElec ? 1.0 : 10.0);
708 p1 = 1.4;
709 }
710 xMax = 10.0;
711 if ((entries > 500.0) && (halfNdof <= 15)) {
712 double peak = rchisqA[halfNdof][detector]->GetBinContent(peakBin);
713 //double peakX = lowEdge + binWidth * peakBin;
714 int binLow = peakBin;
715 double threshold = ((isMuon || isPion || isKaon) ? 0.1 : 0.2) * peak;
716 if (halfNdof >= 10) { threshold *= 2.0; }
717 for (int bin = peakBin; bin <= nBins; ++bin) {
718 if (rchisqA[halfNdof][detector]->GetBinContent(bin) < threshold) break;
719 binLow = bin - 2;
720 if (binLow < 1) { binLow = 1; }
721 }
722 xMin = lowEdge + binWidth * binLow;
723 /*
724 if (isMuon) {
725 if (halfNdof == 1) { xMin = 6.4; }
726 } else {
727 if (halfNdof == 1) { xMin = 8.0; }
728 if (halfNdof == 2) { xMin = (isElec ? 3.8 : 5.4); }
729 if (halfNdof == 3) { xMin = (isElec ? 2.4 : 4.0); }
730 if (halfNdof == 4) { xMin = (isElec ? 0.0 : 3.4); }
731 }
732 */
733 for (int bin = rchisqA[halfNdof][detector]->GetNbinsX(); bin > 0; --bin) {
734 if (rchisqA[halfNdof][detector]->GetBinContent(bin) > 0.0) {
735 int bin2 = bin - (33 - 3 * halfNdof);
736 if (bin2 > bin - 2) { bin2 = bin - 2; }
737 if (bin2 < peakBin) {
738 xMin = 0.0;
739 } else {
740 double xMin2 = lowEdge + binWidth * bin2;
741 if (xMin2 < xMin) {
742 xMin = xMin2;
743 }
744 }
745 break;
746 }
747 }
748 if (isMuon) {
749 if (halfNdof == 1) { xMin = 6.4; }
750 } else {
751 if (halfNdof == 1) { xMin = 8.0; }
752 if (halfNdof == 2) { xMin = (isElec ? 3.8 : 4.6); }
753 if (halfNdof == 3) { xMin = (isElec ? 2.4 : 4.0); }
754 if (halfNdof == 4) { xMin = (isElec ? 0.0 : 3.4); }
755 }
756 }
757 /*
758 if (halfNdof == 1) {
759 rchisqFunc->SetParLimits(1, 0.25, 3.0);
760 } else {
761 rchisqFunc->SetParLimits(1, 0.25, 2.0);
762 }
763 */
764 if (entries < 50.0) {
765 rchisqFunc->FixParameter(1, p1);
766 }
767 if ((entries > 0) && (int((xMax - xMin) / binWidth + 0.5) >= rchisqFunc->GetNumberFreeParameters())) {
768 rchisqA[halfNdof][detector]->Fit("rchisqFunc", "WL", "", xMin, xMax);
769 p0 = rchisqFunc->GetParameter(0);
770 //p0Err = rchisqFunc->GetParError(0);
771 p1 = rchisqFunc->GetParameter(1);
772 //p1Err = rchisqFunc->GetParError(1);
773 }
774
775 double yNorm = 1.0;
776 if ((entries > 0.0) && (rchisqA[halfNdof][detector]->GetEntries() > 0.0)) {
777 yNorm = rchisqA[halfNdof][detector]->GetEntries() / entries;
778 }
779 double yScaleLeft = yNorm * p0Left * p2 / TMath::Gamma(p2);
780 double yScale = yNorm * p0 * p2 / TMath::Gamma(p2);
781 for (int bin = 1; bin <= nBins; ++bin) {
782 double x = lowEdge + (bin - 0.5) * binWidth;
783 double y = rchisqA[halfNdof][detector]->GetBinContent(bin) * yNorm;
784 if ((rchisqA[halfNdof][detector]->GetBinContent(bin) == 0.0) || (x > xMin)) {
785 if ((bin < peakBin) && (xMin > 0.0)) {
786 x *= p1Left * p2;
787 y = yScaleLeft * std::pow(x, p2 - 1.0) * std::exp(-x);
788 } else {
789 x *= p1 * p2;
790 y = yScale * std::pow(x, p2 - 1.0) * std::exp(-x);
791 }
792 }
793 rchisqB[halfNdof][detector]->SetBinContent(bin, y);
794 }
795 for (int bin = nBins + 1; bin <= 100; ++bin) {
796 double x = (lowEdge + (bin - 0.5) * binWidth) * p1 * p2;
797 double y = yScale * std::pow(x, p2 - 1.0) * std::exp(-x);
798 rchisqB[halfNdof][detector]->SetBinContent(bin, y);
799 }
800 double sum = 0.0;
801 for (int bin = 100; bin > 0; --bin) {
802 sum += rchisqB[halfNdof][detector]->GetBinContent(bin);
803 }
804 sum *= binWidth;
805 fprintf(output, " <Tail>\n");
806 fprintf(output, " <Threshold>%16.7E</Threshold>\n", xMin);
807 fprintf(output, " <ScaleY>%16.7E</ScaleY>\n", yScale / sum);
808 fprintf(output, " <ScaleX>%16.7E</ScaleX>\n", p1);
809 fprintf(output, " </Tail>\n");
810 fprintf(output, " <Histogram>\n");
811 for (int bin = 1; bin <= nBins; bin += 10) {
812 fprintf(output, " %16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E%16.7E\n",
813 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin) / sum),
814 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 1) / sum),
815 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 2) / sum),
816 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 3) / sum),
817 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 4) / sum),
818 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 5) / sum),
819 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 6) / sum),
820 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 7) / sum),
821 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 8) / sum),
822 std::log(rchisqB[halfNdof][detector]->GetBinContent(bin + 9) / sum));
823 }
824 fprintf(output, " </Histogram>\n");
825 fprintf(output, " </DegreesOfFreedom>\n");
826 }
827 fprintf(output, " </%s>\n", detectorName[detector]);
828 }
829 fprintf(output, " </TransversePDF>\n");
830 fprintf(output, " </%s>\n", xmllabel);
831
832 //---------------------------------------------------------------------
833
834 TCanvas* canvas = new TCanvas("tree", label, 900, 900);
835 //gStyle->SetPaperSize(1); // BAD BAD BAD!
836 gStyle->SetOptStat(10);
837 gStyle->SetOptFit(0);
838 canvas->GetPad(0)->SetLogy(0);
839 canvas->GetPad(0)->SetGrid(1, 1);
840 canvas->Clear();
841 canvas->Divide(1, 1);
842 canvas->GetPad(1)->SetLogy(1);
843 sprintf(label, "Muid-%s.pdf[", xmllabel); // open file without drawing
844 canvas->SaveAs(label);
845 sprintf(label, "Muid-%s.pdf", xmllabel); // save each canvas to same file
846 canvas->cd(1);
847 discardA->Divide(discardB);
848 discardA->SetLineColor(kRed);
849 discardA->Draw("");
850 canvas->SaveAs(label);
851 for (int outcome = 1; outcome < NOUTCOME; ++outcome) {
852 for (int layer = 0; layer < MuidElementNumbers::getMaximalBarrelLayer() + 1; ++layer) {
853 canvas->cd(1);
854 layerHitU[layer][outcome]->SetMaximum(2.0);
855 layerHitU[layer][outcome]->SetMinimum(1.0E-4);
856 if (isProt or isDeut) layerHitU[layer][outcome]->SetMinimum(1.0E-6);
857 if (isElec) layerHitU[layer][outcome]->SetMinimum(1.0E-10);
858 layerHitU[layer][outcome]->GetXaxis()->SetTitle("KLM Layer");
859 layerHitU[layer][outcome]->GetYaxis()->SetTitle("Frequency");
860 layerHitU[layer][outcome]->Draw("E1");
861 layerHitW[layer][outcome]->SetLineColor(kMagenta);
862 layerHitW[layer][outcome]->Draw("HIST,SAME");
863 layerHitU[layer][outcome]->Draw("E1,SAME");
864 canvas->SaveAs(label);
865 }
866 }
867
868 canvas->Clear();
869 canvas->Divide(1, 2);
870 canvas->GetPad(1)->SetLogy(1);
871 canvas->GetPad(2)->SetLogy(1);
872 gStyle->SetOptStat(110);
873 gStyle->SetOptFit(2);
874 for (int detector = 0; detector <= 2; ++detector) {
875 for (int halfNdof = 1; halfNdof <= 18; ++halfNdof) {
876 int padNo = (halfNdof - 1) % 2 + 1;
877 canvas->cd(padNo);
878 gPad->Clear();
879 rchisqA[halfNdof][detector]->SetMinimum(0.001); // for vertical log axis
880 if (rchisqA[halfNdof][detector]->GetEntries() == 0.0) { rchisqA[halfNdof][detector]->SetMaximum(2.0); }
881 rchisqA[halfNdof][detector]->GetXaxis()->SetTitle("#chi^{2}/ndf");
882 rchisqA[halfNdof][detector]->GetYaxis()->SetTitle("Frequency");
883 rchisqA[halfNdof][detector]->Draw("");
884 rchisqB[halfNdof][detector]->SetLineColor(kMagenta);
885 rchisqB[halfNdof][detector]->Draw("SAME");
886 rchisqA[halfNdof][detector]->Draw("SAME,E1");
887 if (padNo == 2) {
888 canvas->SaveAs(label);
889 }
890 }
891 }
892
893 sprintf(label, "Muid-%s.pdf]", xmllabel); // close file without drawing
894 canvas->SaveAs(label);
895
896 printf("Done.\n");
897
898}
899
900int main(int argc, char** argv)
901{
902 TChain* chain;
903 FILE* output;
904 if (argc < 4)
905 goto badUsage;
906 chain = new TChain;
907 chain->Add(argv[1]);
908 output = std::fopen(argv[2], "a");
909 makeprob(chain, output, argv[3]);
910 std::fclose(output);
911 delete chain;
912 return 0;
913badUsage:
914 std::printf(
915 "Usage:\n"
916 "b2klm-likelihood-parameters input_file output_file xml_label\n");
917 return -1;
918}
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:75
Belle2::MCParticle::getLastDaughter
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
Definition: MCParticle.h:259
Belle2::MCParticle::getMomentum
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition: MCParticle.h:198
Belle2::MCParticle::getFirstDaughter
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
Definition: MCParticle.h:251
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:75
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
VXDTFFilterTest::v1
const std::vector< double > v1
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:28