Belle II Software  release-08-01-10
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 
37 using namespace Belle2;
38 
39 TStyle* 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 
114 void 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 
123 void 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 
900 int 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;
913 badUsage:
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 informations 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
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:91