Port Greg/RD+'s J/psiK reweighting code for RDX run 2

yipengsun commented 2 years ago

Greg has shared his J/psi K reweighting code at:

/afs/cern.ch/user/g/gciezare/public/forRDRun2/JpsiK

We need to port this to RDX run 2.

Current status

[x] Produce J/psi K data sample ntuples locally
[x] Produce J/psi K MC sample ntuples locally
[x] Update Greg's selection
[x] Run Greg's fit
[x] Submit J/psi K data production to the GRID
[x] Submit J/psi K MC (12143001) production to the GRID
[x] Add standard weights to MC
[x] Compute efficiency ratios in the following sequence:
- b_OWNPV_NDOF and b_nTracks
- b_P and b_ETA
[x] Integrate this weight to our workflow

Possible improvements

Update Greg's selection
- Currently the selections are mostly aligned w/ run 1 analysis. We may need to update the selections once we finish our selection optimization
Apply a DiMuon weight. See here for more details
Apply an efficiency correction based on DiMuon / All Trigger efficiency, binned in P,PT of the B meson
Use an Ipatia function to describe signal
Apply a nSPDhits < 450 cut for data, a weight for MC. This can be derived s.t. nSPDhits = nSPDhits(nTracks) by looking at the L0DiMuon trigger line
- RD+ did this
- Not sure if we need it, because our L0HadronTOS effectively does this

Validations

Stripping: See line here. There's actually a PIDK > 0 in the K requirement. It is not ideal, but we have a tighter PIDK > 4 in the offline cut, and hopefully the generated weights will make this baked-in cut non-effective
Fit quality: See this comment
sWeight test: See this post

References

Slides on sWeights: https://indico.cern.ch/event/940874/contributions/3953530/attachments/2127199/3581560/sweights_ms.pdf
Paper: https://arxiv.org/abs/physics/0402083

yipengsun commented 2 years ago

I don't even see a stripping line name in Greg's code. Maybe this is intended for J/psi? Anyway I'm going back to RD+'s code for some cross reference.

yipengsun commented 2 years ago

Don't see anything useful by searching J/psi in RD+'s code.

yipengsun commented 2 years ago

Now Greg's script runs. But it leaves the output tree empty, suggesting NO event is being selected. I think the most likely case is that the J/psi K requires a specific input data.

yipengsun commented 2 years ago

There's a CHARM.MDST that looks promising. I think I'd like to ask Greg the following questions:

Which data files should I be looking at?
Confirm that indeed we don't need to filtering on a particular stripping line.

yipengsun commented 2 years ago

From RD+'s ANA, Sec 7.2.1, we see that the stripping line is BetaSBu2JpsiKDetachedLine, from a DIMUON sample. The MC sample used is 12143001.

The LFN on DIRAC is:

sim+std://LHCb/Collision16/Beam6500GeV-VeloClosed-MagDown/Real Data/Reco16/Stripping28r2p1/90000000/DIMUON.DST

yipengsun commented 2 years ago

I tried to use the DIMUON with the correct stripping line. The input .DST files have ~82k events, yet the very next selection has 0 event. I think this is due to some TES location changes.

yipengsun commented 2 years ago

Hmm, turns out that when I'm using latest stripping for 2016: v28r2p1, even for RDX the output ntuple was empty. I think we either stick with v28r2 or update our DaVinci to v46 (which should be the final version for run 1 & 2).

@manuelfs What do you think?

yipengsun commented 2 years ago

Personally I'm inclined to update DaVinci. But then the problem is: What should we do for the already produced tracker-only ntuples (with DaVinci v45r6)?

yipengsun commented 2 years ago

Or, it could be just there's a TES location change, and our DaVinci v45r6 would still work.

yipengsun commented 2 years ago

To summarize:

Our local DaVinci v45r6 can't produce ntuples with latest stripping (v28r2p1, done in winter 2021, so very very recent)
I have 3 hypotheses:
1. The latest stripping changed TES location (doubtful, as it's an incremental stripping, and I don't see a reason to change TES other than forward compatibility)
  - If this is the case, we just update TES location and everything's good
2. The LOCAL DaVinci is problematic because it lacks latest data tag (DaVinci v45r6 was released end of 2020, the stripping was end of 2021, so we definitely don't have the tag locally)
  - Then we can just update our DaVinci docker image to v46 so that we can produce ntuples locally, while still use v45r6 on the GRID
  - I think this is the most likely scenario
3. The current DaVinci (v45r6) is fundamentally incompatible with latest stripping
  - Then we have to decide if we want to update DaVinci, or stick with older stripping
  - If we choose to update DaVinci, we need to think about what should we do for the already produced TO MC ntuples

FYI @manuelfs

yipengsun commented 2 years ago

This is the announcement for the v28r2p1 incremental stripping: mailing list archive.

I think the idea is that our lines (SL and J/psi K lines) don't get changed in this incremental stripping, therefore they don't exist in v28r2p1. I should try v28r2 and see if it works.

yipengsun commented 2 years ago

Tried the v28r2 DIMUON .DST files. Still doesn't work. The main problem is that DaVinci can't see any event:

TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
TimingAuditor.T...   INFO This machine has a speed about   3.85 times the speed of a 2.8 GHz Xeon.
TimingAuditor.T...   INFO Algorithm                              (millisec) |    <user> |   <clock> |      min       max    sigma | entries | total (s) |
TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
TimingAuditor.T...   INFO EVENT LOOP                                        |     0.042 |     0.042 |    0.019     329.5     1.83 |   82250 |     3.493 |
TimingAuditor.T...   INFO  DaVinciEventSeq                                  |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   DaVinciInitAlg                                  |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   FilteredEventSeq                                |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    DaVinciEventInitSeq                            |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     PhysInitSeq                                   |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     AnalysisInitSeq                               |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    DaVinciAnalysisSeq                             |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     DaVinciUserSequence                           |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO      MyProtoPMaker                                |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO      StdNoPIDsVeloPions                           |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO      TupleBminus                                  |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     MonitoringSequence                            |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   LumiSeq                                         |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    EventAccount                                   |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   IntegratorSeq                                   |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    IntegrateBeamCrossing                          |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    GetIntegratedLuminosity                        |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------

manuelfs commented 2 years ago

Will provided his J/Psi K script https://gitlab.cern.ch/LHCb-BnoC/lhcb-ana-2020-022-analysis-preservation/-/blob/master/TUPLES/ntuple_Bjpsikstrip.py

yipengsun commented 2 years ago

The only material difference between Will's script and Greg's is the stripping line: Will is using Bs2MuMuLinesBu2JPsiKLine.

I tried to swap stripping line to use Will's, and I still don't get any event:

TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
TimingAuditor.T...   INFO This machine has a speed about   4.17 times the speed of a 2.8 GHz Xeon.
TimingAuditor.T...   INFO Algorithm                              (millisec) |    <user> |   <clock> |      min       max    sigma | entries | total (s) |
TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
TimingAuditor.T...   INFO EVENT LOOP                                        |     0.041 |     0.042 |    0.019     335.0     1.81 |   82250 |     3.481 |
TimingAuditor.T...   INFO  DaVinciEventSeq                                  |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   DaVinciInitAlg                                  |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   FilteredEventSeq                                |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    DaVinciEventInitSeq                            |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     PhysInitSeq                                   |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     AnalysisInitSeq                               |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    DaVinciAnalysisSeq                             |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     DaVinciUserSequence                           |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO      MyProtoPMaker                                |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO      StdNoPIDsVeloPions                           |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO      TupleBminus                                  |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO     MonitoringSequence                            |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   LumiSeq                                         |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    EventAccount                                   |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO   IntegratorSeq                                   |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    IntegrateBeamCrossing                          |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO    GetIntegratedLuminosity                        |     0.000 |     0.000 |    0.000       0.0     0.00 |       0 |     0.000 |
TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
NTupleSvc            INFO NTuples saved successfully

Next I'll try to just download Will's DV script and do minimal changes to see if it works locally.

yipengsun commented 2 years ago

I think the most interesting lines from the logs are:

DaVinciInitAlg    SUCCESS ==================================================================
DaVinciInitAlg    SUCCESS 82250 events processed
DaVinciInitAlg    SUCCESS ==================================================================
StdNoPIDsVeloPions   INFO ********************************************************************************
StdNoPIDsVeloPions   INFO  created 'pi+' and 'pi- : 0 per 0 calls (0+-0)/event
GetIntegratedLu...   INFO The FSR were checked : 1
GetIntegratedLu...   INFO number of events seen: 0
GetIntegratedLu...   INFO Integrated Luminosity  : 0 +/-  0 [pb-1]

So DaVinci is NOT REALLY seeing any event.

yipengsun commented 2 years ago

I finally found the problem: The previous dst files don't contain OUR line!

This is unexpected, as they contain 82250 events and they are listed under the same DIMUON bookkeeping location. Now I just randomly select a few more in the same LFN location and do it by trial-and-error:

TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
TimingAuditor.T...   INFO This machine has a speed about   4.00 times the speed of a 2.8 GHz Xeon.
TimingAuditor.T...   INFO Algorithm                              (millisec) |    <user> |   <clock> |      min       max    sigma | entries | total (s) |
TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------
TimingAuditor.T...   INFO EVENT LOOP                                        |     4.897 |     4.898 |    0.020    1350.9    10.45 |  113254 |   554.797 |
TimingAuditor.T...   INFO  DaVinciEventSeq                                  |     7.215 |     7.211 |    0.712    1350.8    11.93 |   76051 |   548.444 |
TimingAuditor.T...   INFO   DaVinciInitAlg                                  |     0.029 |     0.033 |    0.025       0.6     0.01 |   76051 |     2.577 |
TimingAuditor.T...   INFO   FilteredEventSeq                                |     7.161 |     7.157 |    0.663    1350.8    11.93 |   76051 |   544.323 |
TimingAuditor.T...   INFO    DaVinciEventInitSeq                            |     0.004 |     0.005 |    0.004       0.1     0.00 |   76051 |     0.393 |
TimingAuditor.T...   INFO     PhysInitSeq                                   |     0.000 |     0.001 |    0.001       0.0     0.00 |   76051 |     0.109 |
TimingAuditor.T...   INFO     AnalysisInitSeq                               |     0.000 |     0.001 |    0.001       0.1     0.00 |   76051 |     0.103 |
TimingAuditor.T...   INFO    DaVinciAnalysisSeq                             |     7.154 |     7.149 |    0.656    1350.8    11.93 |   76051 |   543.730 |
TimingAuditor.T...   INFO     DaVinciUserSequence                           |     7.150 |     7.145 |    0.653    1350.8    11.93 |   76051 |   543.420 |
TimingAuditor.T...   INFO      MyProtoPMaker                                |     3.533 |     3.532 |    0.176      31.6     4.87 |   76051 |   268.628 |
TimingAuditor.T...   INFO      StdNoPIDsVeloPions                           |     0.029 |     0.033 |    0.001       0.7     0.02 |   76051 |     2.580 |
TimingAuditor.T...   INFO      TupleBminus                                  |     0.910 |     0.913 |    0.001    1347.6    10.34 |   76051 |    69.457 |
TimingAuditor.T...   INFO     MonitoringSequence                            |     0.000 |     0.001 |    0.001       0.0     0.00 |   76051 |     0.109 |
TimingAuditor.T...   INFO   LumiSeq                                         |     0.003 |     0.004 |    0.003       0.1     0.00 |   76051 |     0.315 |
TimingAuditor.T...   INFO    EventAccount                                   |     0.000 |     0.001 |    0.001       0.0     0.00 |   76051 |     0.109 |
TimingAuditor.T...   INFO   IntegratorSeq                                   |     0.011 |     0.010 |    0.008       0.1     0.00 |   76051 |     0.793 |
TimingAuditor.T...   INFO    IntegrateBeamCrossing                          |     0.001 |     0.001 |    0.001       0.0     0.00 |   76051 |     0.122 |
TimingAuditor.T...   INFO    GetIntegratedLuminosity                        |     0.000 |     0.001 |    0.001       0.0     0.00 |   76051 |     0.102 |
TimingAuditor.T...   INFO *      UnpackTrack                                |     3.259 |     3.257 |    0.042      31.3     4.87 |   76051 |   247.754 |
TimingAuditor.T...   INFO *     UnpackRecVertex                             |     0.628 |     0.620 |    0.022      19.6     2.05 |   75987 |    47.160 |
TimingAuditor.T...   INFO *      unpackFittedVeloTracks                     |     0.542 |     0.537 |    0.016      15.2     2.03 |   75953 |    40.819 |
TimingAuditor.T...   INFO *     Dimuon_PsAndVsUnpack                        |     0.050 |     0.051 |    0.001       3.4     0.04 |  608911 |    31.173 |
TimingAuditor.T...   INFO *      UnpackChargedProtosSeq                     |     6.789 |     7.236 |    0.343     195.7     7.23 |     947 |     6.853 |
TimingAuditor.T...   INFO *       UnpackChargedProtos                       |     5.417 |     5.740 |    0.075     194.3     7.20 |     947 |     5.436 |
TimingAuditor.T...   INFO *       ProtoParticleCombDLLs                     |     1.372 |     1.488 |    0.230       5.2     0.40 |     947 |     1.409 |
TimingAuditor.T...   INFO *        CheckChargedProtosExist                  |     0.000 |     0.003 |    0.002       0.0     0.00 |     947 |     0.003 |
TimingAuditor.T...   INFO *        ChargedProtoPAddMuon                     |     0.401 |     0.403 |    0.093       3.8     0.17 |     947 |     0.383 |
TimingAuditor.T...   INFO *        ChargedProtoPAddRich                     |     0.908 |     0.991 |    0.088       4.1     0.30 |     947 |     0.939 |
TimingAuditor.T...   INFO *        ChargedProtoPCombDLL                     |     0.052 |     0.080 |    0.018       0.3     0.03 |     947 |     0.076 |
TimingAuditor.T...   INFO *         UnpackMuonPIDSeq                        |     0.348 |     0.315 |    0.042       3.7     0.17 |     947 |     0.299 |
TimingAuditor.T...   INFO *          UnpackMuonPIDs                         |     0.327 |     0.305 |    0.034       3.6     0.16 |     947 |     0.290 |
TimingAuditor.T...   INFO *          CheckMuonPIDs                          |     0.021 |     0.003 |    0.003       0.0     0.00 |     947 |     0.003 |
TimingAuditor.T...   INFO *         UnpackRichPIDSeq                        |     0.802 |     0.896 |    0.044       4.0     0.29 |     947 |     0.849 |
TimingAuditor.T...   INFO *          UnpackRichPIDs                         |     0.791 |     0.887 |    0.037       4.0     0.29 |     947 |     0.840 |
TimingAuditor.T...   INFO *          CheckRichPIDs                          |     0.000 |     0.002 |    0.002       0.0     0.00 |     947 |     0.002 |
TimingAuditor.T...   INFO *      UnpackMuonTracks                           |     0.454 |     0.497 |    0.025      15.0     0.50 |     947 |     0.471 |
TimingAuditor.T...   INFO *      L0DecReportsMaker                          |     4.403 |     4.668 |    0.271      21.5     2.60 |     947 |     4.421 |
TimingAuditor.T...   INFO *       L0DUFromRaw                               |     4.382 |     4.630 |    0.246      21.5     2.58 |     947 |     4.385 |
TimingAuditor.T...   INFO *      L0SelReportsMaker                          |     1.045 |     0.979 |    0.274     607.1    19.72 |     947 |     0.927 |
TimingAuditor.T...   INFO *       L0MuonFromRaw                             |     0.190 |     0.173 |    0.137       0.8     0.03 |     947 |     0.164 |
TimingAuditor.T...   INFO *       L0CaloFromRaw                             |     0.073 |     0.053 |    0.041       0.7     0.03 |     947 |     0.051 |
TimingAuditor.T...   INFO *      Hlt1DecReportsDecoder                      |     0.031 |     0.023 |    0.020       0.2     0.01 |     947 |     0.022 |
TimingAuditor.T...   INFO *      Hlt1SelReportsDecoder                      |     0.084 |     0.066 |    0.022       0.4     0.03 |     947 |     0.063 |
TimingAuditor.T...   INFO *      Hlt2SelReportsDecoder                      |     0.623 |     0.735 |    0.043     415.0    13.48 |     947 |     0.696 |
TimingAuditor.T...   INFO *      Hlt2DecReportsDecoder                      |     0.168 |     0.145 |    0.134       0.8     0.03 |     947 |     0.137 |
TimingAuditor.T...   INFO *      StdAllNoPIDsPions                          |     0.095 |     0.058 |    0.013       0.2     0.03 |     947 |     0.056 |
TimingAuditor.T...   INFO *      StdNoPIDsUpPions                           |     0.000 |     0.022 |    0.003       0.1     0.01 |     947 |     0.021 |
TimingAuditor.T...   INFO *      UnpackMCParticles                          |     0.018 |     0.019 |    0.007       0.3     0.03 |    8276 |     0.158 |
TimingAuditor.T...   INFO *      createODIN                                 |     0.021 |     0.009 |    0.007       0.0     0.00 |     947 |     0.009 |
TimingAuditor.T...   INFO ----------------------------------------------------------------------------------------------------------------------

Now we successfully reco'ed 947 events out of 113254.

The lesson here is: For the latest stripping (v28r2), the stripping lines are NOT EVENLY distributed (s.t. multiple files may not contain the line we needed)!

yipengsun commented 2 years ago

Will's line: Bs2MuMuLinesBu2JPsiKLine
Greg/RD+'s line: BetaSBu2JpsiKDetachedLine

@manuelfs

manuelfs commented 2 years ago

Will's B -> J/Psi K ntuples are at /eos/lhcb/wg/BnoC/B2Kpi0_ANA-2020-22. The kaon goes down to 250 MeV, and the J/Psi is not clear. However, it looks like they don't have L0MuonDecision, so if the separation between L0 types is important, these ntuples can only be useful to play with them

yipengsun commented 2 years ago

I think the path forward would be:

Play around my local production to make sure all required trigger branches are there
Set up the fit and sWeight procedure locally with my local ntuples
Once both are working, submit a GRID job

yipengsun commented 2 years ago

Found a very interesting sWeight python script: https://gist.github.com/alexpearce/9230249003efbc8eae51

yipengsun commented 2 years ago

Some interesting RooFit sample code:

yipengsun commented 2 years ago

DiMuon weight

The weight depends on variable sqrt(PT_mup * PT_mum).

This looks like a momentum correction based on trigger categories, but I am VERY CONFUSED on Greg's > -0.5 cut:

  data->Draw("sqrt(muplus_PT*muminus_PT)>>allPT",cutstringData+"&&muminus_L0MuonDecision_TOS > -0.5");
  data->Draw("sqrt(muplus_PT*muminus_PT)>>TOSPT",cutstringData+"&&Jpsi_L0DiMuonDecision_TOS > 0.5");
  TH1F dataEff = TOSPT / allPT;

  mc->Draw("sqrt(muplus_PT*muminus_PT)>>allPTMC","("+cutstringMC+"&&muminus_L0MuonDecision_TOS > -0.5"+")*"+MCweight);
  mc->Draw("sqrt(muplus_PT*muminus_PT)>>TOSPTMC","("+cutstringMC+"&&Jpsi_L0DiMuonDecision_TOS > 0.5"+")*"+MCweight);
  TH1F mcEff = TOSPTMC / allPTMC;

I also don't see the motivation. A naive guess would be: Take a look at all track PT vs. PT with some L0 trigger cuts, and see that some of the efficiencies are not modeled well in MC.

yipengsun commented 2 years ago

Greg also reweight different trigger categories differently:

   categories.push_back("Bplus_L0MuonDecision_TIS > 0.5");
   categories.push_back("Bplus_L0Global_TIS > 0.5 && Bplus_L0MuonDecision_TIS < 0.5");
   categories.push_back("Bplus_L0Global_TIS < 0.5 ");

But for run 2, the only trigger on B is L0Global TIS, so I think we don't need to reweight by category.

yipengsun commented 2 years ago

I feel run 2 RD+'s reweighting strategy would be more useful to us in this case, because we are using the same trigger paths.

yipengsun commented 2 years ago

This is what fit pdfs look like before the fit:

22_02_16_03_56-fit-JpsiK-init

No wonder the fit doesn't converge.

I'm also confused by the tail component. I think the Crystal Ball powerlaw tail might be sufficient already.

yipengsun commented 2 years ago

Fit validation

I added another Gaussian for the signal so that now it's double Gaussian + Crystal Ball. I also used finer binning. The fit now looks like this:

fit_final_log_scale

The fit results are:

╒═════════╤═════════════╤══════════════════╤═════════╤═════════════╕
│ valid   │ converged   │ param at limit   │ edm     │ min value   │
╞═════════╪═════════════╪══════════════════╪═════════╪═════════════╡
│ True    │ True        │ False            │ 0.00011 │ -3454       │
╘═════════╧═════════════╧══════════════════╧═════════╧═════════════╛

Parameters
name               value         minuit_minos    at limit
-------------  ---------  -------------------  ----------
yld_bkg_ratio     0.2322  - 0.0021   +  0.002       False
yld_sig_ratio     0.7678  - 0.0021   + 0.0022       False
expc           -0.001743  -8.3e-05   +  8e-05       False
frac_sig_cb       0.4252  -  0.028   +  0.031       False
frac_sig_g          0.53  -  0.033   +   0.03       False
alpha              1.857  -  0.045   +  0.049       False
peak                5280  -  0.011   +  0.011       False
n_cb                 150  -1.1e+02   +  0.022       False
width_cb           11.66  -   0.28   +   0.28       False
width_g_sig        7.351  -    0.1   +  0.093       False
width_g2_sig       34.03  -    2.4   +    2.5       False

Old

@afernez Here's the updated fit plot:

fit_final_log_scale

PDF version: fit_final_log_scale.pdf

The signal is modeled by a Gaussian + Crystal Ball, background just an exponential.

The final fit result is:

╒═════════╤═════════════╤══════════════════╤═════════╤═════════════╕
│ valid   │ converged   │ param at limit   │ edm     │ min value   │
╞═════════╪═════════════╪══════════════════╪═════════╪═════════════╡
│ True    │ True        │ False            │ 5.1e-05 │ -1.454e+04  │
╘═════════╧═════════════╧══════════════════╧═════════╧═════════════╛

Parameters
name               value         minuit_minos    at limit
-------------  ---------  -------------------  ----------
yld_bkg_ratio     0.2429  -0.00081   +0.00079       False
yld_sig_ratio     0.7571  -  0.001   +  0.001       False
expc           -0.001164  -3.9e-05   +3.9e-05       False
frac_sig_cb_g     0.3181  -  0.012   +  0.013       False
alpha              1.711  -  0.049   +  0.048       False
peak                5280  -  0.011   +  0.011       False
n_cb               6.907  -    1.6   +    2.6       False
width_cb           13.74  -    0.2   +    0.2       False
width_g_sig        7.756  -  0.047   +  0.045       False

yipengsun commented 2 years ago

sWeights test

If we naively apply sWeights back to B mass on real data, we get the following plot:

splot_sig

In short, we cannot apply sWeights back to the mass, because we derived sWeights from the mass fit, and sWeights only works on (totally) uncorrelated variables (r.s.t. the variables deriving it).

For more info, see this paper, and this forum post.

yipengsun commented 2 years ago

Interesting that some of the data bins are negative after applying sWeights. Take a look at B_OWNPV_NDOF & nTracks:

index	data yld	MC yld	data/MC wt	expected MC yld rwt	MC yld rwt	diff
(10, 2)	-0.130365	0.80233	-0.0771535	-0.0619026	-801.528	801.466
(19, 5)	-0.244527	0	0	0	0	0

yipengsun commented 2 years ago

These bins have large relative differences:

index	data yld	MC yld	data/MC wt	expected MC yld rwt	MC yld rwt	rel diff
(1, 19)	102.906	9.06387	5.37239	48.6946	45.47	0.0709176
(2, 17)	203.276	22.8756	4.20488	96.189	97.251	-0.0109201
(5, 19)	335.221	24.2448	6.54264	158.625	154.528	0.0265115
(6, 19)	420.071	30.2384	6.5736	198.775	187.253	0.0615321
(7, 19)	459.643	29.1825	7.45311	217.501	212.331	0.0243488
(12, 19)	588.037	33.522	8.3007	278.256	256.993	0.082739
(13, 19)	573.458	29.9832	9.0503	271.357	264.044	0.0276978
(14, 19)	557.77	32.7843	8.05061	263.934	238.293	0.107602
(15, 19)	509.507	38.8194	6.2107	241.096	232.854	0.0353943
(16, 19)	445.676	24.383	8.64912	210.891	187.405	0.125322
(17, 19)	357.053	20.1945	8.3664	168.956	155.58	0.0859727
(19, 18)	349.68	29.9247	5.52943	165.467	167.753	-0.0136327
(19, 19)	254.492	15.1202	7.96446	120.424	114.577	0.0510339

yipengsun commented 2 years ago

One difference between ROOT's TH1 and numpy's histogram:

ROOT excludes the maximum edge: https://root.cern.ch/doc/master/classTH1.html#binning
numpy includes it: https://numpy.org/doc/stable/reference/generated/numpy.histogram.html

yipengsun commented 2 years ago

What I have checked so far:

numpy and ROOT can get weights consistently: https://github.com/umd-lhcb/lhcb-ntuples-gen/tree/master/studies/test-numpy_vs_root_histo
numpy and ROOT can build histogram consistently (provided that the open-close-ness of the last bin doesn't come into play): https://github.com/umd-lhcb/lhcb-ntuples-gen/tree/master/studies/test-JpsiK_weight_apply

So, perhaps my weight application script has a bug. I'll look into that.

yipengsun commented 2 years ago

Actually, in the step-2 babymaker template, I used tree->BuildIndex("runNumber", "eventNumber"). Given that the writings of output is sequential, so a re-indexing is not really needed, maybe I should just comment out that line.

Edit: I checked 2 auxiliary ntuples with numpy, and the branch orderings are consistent between the 2.

I also checked the duplicated UIDs of one of the ntuples:

Num of events: 2478, Num of IDs: 2375, Num of UIDs: 2280
Num of duplicated IDs: 95, Num of duplicated events: 103, duplicate rate: 4.16%

Given that runNumber-eventNumber is not unique, it will likely confuse the BuildIndex("runNumber", "eventNumber")

yipengsun commented 2 years ago

I disabled the index building for friend trees. Now I got something like this:

index	data yld	MC yld	data/MC wt	exp. MC yld rwt	MC yld rwt	rel diff
(0, 0)	476.835	2.33769	0.313152	0.732052	0.731825	0.000310958
(0, 1)	874.139	2.81657	0.51863	1.46076	1.46031	0.000310958
(0, 2)	612.826	1.5432	0.562428	0.86794	0.86767	0.000310958
(0, 3)	530.333	1.58764	0.626665	0.994917	0.994608	0.000310958

Noting that rel diff = (hMc[i, j] * hRatio[i, j] - hResult[i, j]) / hResult[i, j] = const

So, hMC[i, j] = (1 + const) / hRatio[i, j] * hResult[i, j]. I feel there's misalignment by 1 error somewhere.

yipengsun commented 2 years ago

I think I have a better idea now:

The inconsistencies are due to the fact that I'm NOT using the same histogram between result comparison and step-2 ntuple postprocessing! This is because I tried to be smart in the postprocessing ntuple and match files in a fuzzy way, which ends up picking up the old, archived ntuple containing outdated histograms!

yipengsun commented 2 years ago

After carefully handle all possible inconsistencies beween ROOT and numpy histograms, here's the result after reweighting:

`B_OWNPV_NDOF` & `nTracks`

b_ownpv_ndof ntracks

`B_ETA` & `B_PT`

b_eta b_pt

yipengsun commented 2 years ago

Consider this done, noting that for events outside the binning range, we probably want to apply weights of its nearest neighbor.

yipengsun commented 2 years ago

I'll check the distributions of related variables for both J/psi K and RDX data to see if we need extra bins.

yipengsun commented 2 years ago

@manuelfs I documented current J/psi K reweighting specifications at https://github.com/umd-lhcb/rdx-run2-analysis/blob/master/docs/reweight/JpsiK_reweight.md

In that doc, I listed the binning "efficiencies", that is, percent of events out of binning range.

I think we could make the PV NDOF and PT max bin larger, and keep the rest as-is, this way, we can minimize the event loss, and we DON'T use the nearest bins to reduce systematic uncertanties.

Note that some of the variables (e.g. PT of the B mesons) can have very large values (104.7 actual, binned at 25 GeV), but there's only a few lying outside (If I bin at 30 GeV, only 0.149% of events is outside the binning range), so putting all those outliers in a single overflow bin is perhaps not a good idea.

manuelfs commented 2 years ago

Thanks for the summary. Why do think the systematic uncertainties from applying the correction found in the last bin to the overflow would be large?

yipengsun commented 2 years ago

Because the overflow bin covers a very large kinematic region. For example, the PT bin is [0, 25 GeV], with 20 bins. The overflow bin (a single bin) would be [25 GeV, 100+ GeV].

If we apply the weights from the last bin, the assumption would be, the weights are more or less constant in this binning range, and the range in the last nominal bin, which is probably false, because we used a much finer binning for the nominal region.

Edit: Missed 1 crucial point.

manuelfs commented 2 years ago

That is certainly a fair concern: we can't really know what's happening in those ranges since it hasn't been measured. However, we do expect the corrections to vary reasonably smooth. For instance, if the correction is mostly small or tends to a constant value. Thus, in many cases we can make the case that the last bin is a good approximation within a reasonably small uncertainty.

In any case, what is your suggestion, to cut those events? If I understand this correctly we cannot do that as we do not know the B_PT or B_ETA in B -> D(*) ell nu data.

yipengsun commented 2 years ago

My suggestions was: enlarge the ranges for the more significantly-affected variables, then cut the rest off. But indeed we cannot meaningfully cut off events, because the B_PT, B_ETA are NOT measured very accurately in real data.

So, I think maybe we should:

Enlarge the ranges for B_PT and B_ETA
Apply a cut on B_PV_NDOF for real data
Use nearest bins for the rest of the variables.

manuelfs commented 2 years ago

For 1 I guess you mean B_PT and B_ETA and nearest bin for nTracks? Can you add the ratio plot to the plots above so that we can see the trend in the corrections?

yipengsun commented 2 years ago

I'd also look at the weight histograms directly to see the bin-by-bin trend (disregarding the error bars):

(The nTracks & PV NDOF is varying much more smoothly)

`nTracks & PV NDOF`

nTracks \ PV NDOF	1 (1.0,10.9)	2 (20.9)	3 (30.8)	4 (40.8)	5 (50.8)	6 (60.7)	7 (70.6)	8 (80.6)	9 (90.5)	10 (100.5)	11 (110.4)	12 (120.4)	13 (130.3)	14 (140.3)	15 (150.2)	16 (160.2)	17 (170.1)	18 (180.1)	19 (190.0)	20 (200.0)
1 (0.0,22.5)	0.31 ± 0.56	0.30 ± 0.55	0.26 ± 0.51	0.22 ± 0.47	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
2 (45.0)	0.52 ± 0.72	0.45 ± 0.67	0.40 ± 0.63	0.35 ± 0.60	0.32 ± 0.57	0.28 ± 0.53	0.16 ± 0.40	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
3 (67.5)	0.56 ± 0.75	0.65 ± 0.81	0.72 ± 0.85	0.68 ± 0.82	0.59 ± 0.77	0.50 ± 0.71	0.41 ± 0.64	0.31 ± 0.56	0.28 ± 0.53	0.32 ± 0.56	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
4 (90.0)	0.63 ± 0.79	0.66 ± 0.81	0.79 ± 0.89	0.93 ± 0.96	0.96 ± 0.98	0.89 ± 0.94	0.69 ± 0.83	0.55 ± 0.74	0.42 ± 0.65	0.35 ± 0.59	0.28 ± 0.53	0.20 ± 0.45	0.29 ± 0.54	0.22 ± 0.47	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
5 (112.5)	0.59 ± 0.76	0.64 ± 0.80	0.81 ± 0.90	0.96 ± 0.98	1.09 ± 1.04	1.10 ± 1.05	1.06 ± 1.03	0.91 ± 0.96	0.74 ± 0.86	0.59 ± 0.77	0.44 ± 0.66	0.35 ± 0.59	0.28 ± 0.53	0.23 ± 0.48	0.15 ± 0.38	0.17 ± 0.42	0.31 ± 0.56	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
6 (135.0)	0.67 ± 0.82	0.64 ± 0.80	0.80 ± 0.90	0.98 ± 0.99	1.11 ± 1.05	1.26 ± 1.12	1.26 ± 1.12	1.23 ± 1.11	1.12 ± 1.06	0.94 ± 0.97	0.77 ± 0.88	0.59 ± 0.77	0.46 ± 0.68	0.36 ± 0.60	0.28 ± 0.53	0.22 ± 0.47	0.18 ± 0.42	0.16 ± 0.40	0.16 ± 0.39	0.00 ± 0.00
7 (157.5)	0.73 ± 0.85	0.73 ± 0.85	0.81 ± 0.90	0.95 ± 0.97	1.13 ± 1.06	1.26 ± 1.12	1.33 ± 1.15	1.36 ± 1.17	1.33 ± 1.15	1.21 ± 1.10	1.16 ± 1.08	1.00 ± 1.00	0.79 ± 0.89	0.60 ± 0.78	0.49 ± 0.70	0.37 ± 0.61	0.29 ± 0.54	0.22 ± 0.47	0.17 ± 0.42	0.10 ± 0.31
8 (180.0)	0.81 ± 0.90	0.83 ± 0.91	0.91 ± 0.96	1.07 ± 1.03	1.16 ± 1.07	1.22 ± 1.11	1.33 ± 1.15	1.42 ± 1.19	1.45 ± 1.20	1.44 ± 1.20	1.36 ± 1.17	1.26 ± 1.12	1.21 ± 1.10	1.02 ± 1.01	0.85 ± 0.92	0.70 ± 0.84	0.53 ± 0.73	0.41 ± 0.64	0.32 ± 0.56	0.23 ± 0.48
9 (202.5)	1.19 ± 1.09	0.88 ± 0.94	0.99 ± 1.00	1.18 ± 1.09	1.30 ± 1.14	1.28 ± 1.13	1.40 ± 1.18	1.56 ± 1.25	1.53 ± 1.24	1.52 ± 1.23	1.51 ± 1.23	1.47 ± 1.21	1.44 ± 1.20	1.29 ± 1.14	1.24 ± 1.11	1.08 ± 1.04	0.94 ± 0.97	0.72 ± 0.85	0.59 ± 0.77	0.46 ± 0.68
10 (225.0)	1.36 ± 1.17	1.06 ± 1.03	1.23 ± 1.11	1.36 ± 1.17	1.37 ± 1.17	1.47 ± 1.21	1.46 ± 1.21	1.68 ± 1.30	1.66 ± 1.29	1.68 ± 1.30	1.66 ± 1.29	1.63 ± 1.28	1.53 ± 1.24	1.53 ± 1.23	1.42 ± 1.19	1.45 ± 1.20	1.33 ± 1.15	1.22 ± 1.10	0.98 ± 0.99	0.90 ± 0.95
11 (247.5)	1.58 ± 1.26	1.14 ± 1.07	1.29 ± 1.14	1.48 ± 1.22	1.61 ± 1.27	1.57 ± 1.25	1.71 ± 1.31	1.61 ± 1.27	1.65 ± 1.29	1.78 ± 1.33	1.75 ± 1.32	1.64 ± 1.28	1.65 ± 1.29	1.66 ± 1.29	1.61 ± 1.27	1.66 ± 1.29	1.50 ± 1.23	1.44 ± 1.20	1.38 ± 1.17	1.34 ± 1.16
12 (270.0)	1.62 ± 1.27	1.40 ± 1.18	1.50 ± 1.23	1.63 ± 1.27	1.66 ± 1.29	1.93 ± 1.39	1.91 ± 1.38	1.95 ± 1.40	1.85 ± 1.36	2.07 ± 1.44	2.00 ± 1.41	1.95 ± 1.40	1.94 ± 1.39	1.91 ± 1.38	1.84 ± 1.36	1.81 ± 1.35	1.83 ± 1.35	1.78 ± 1.33	1.64 ± 1.28	1.57 ± 1.25
13 (292.5)	3.31 ± 1.82	1.45 ± 1.20	1.84 ± 1.35	2.16 ± 1.47	2.11 ± 1.45	2.07 ± 1.44	2.63 ± 1.62	2.63 ± 1.62	2.38 ± 1.54	2.50 ± 1.58	2.17 ± 1.47	2.46 ± 1.57	2.22 ± 1.49	2.30 ± 1.52	2.21 ± 1.49	2.26 ± 1.50	2.18 ± 1.48	2.13 ± 1.46	2.31 ± 1.52	2.11 ± 1.45
14 (315.0)	2.75 ± 1.66	1.96 ± 1.40	2.50 ± 1.58	2.74 ± 1.66	3.01 ± 1.73	2.68 ± 1.64	3.08 ± 1.75	3.04 ± 1.74	2.74 ± 1.65	2.55 ± 1.60	3.17 ± 1.78	2.95 ± 1.72	2.72 ± 1.65	2.65 ± 1.63	2.83 ± 1.68	2.73 ± 1.65	2.50 ± 1.58	2.75 ± 1.66	2.58 ± 1.61	2.68 ± 1.64
15 (337.5)	2.48 ± 1.57	2.90 ± 1.70	2.77 ± 1.66	4.05 ± 2.01	3.95 ± 1.99	3.62 ± 1.90	4.38 ± 2.09	3.59 ± 1.90	3.72 ± 1.93	3.55 ± 1.89	3.84 ± 1.96	3.62 ± 1.90	3.64 ± 1.91	3.27 ± 1.81	3.26 ± 1.81	3.18 ± 1.78	3.23 ± 1.80	3.23 ± 1.80	4.15 ± 2.04	3.68 ± 1.92
16 (360.0)	4.50 ± 2.12	3.13 ± 1.77	2.75 ± 1.66	4.52 ± 2.12	3.76 ± 1.94	4.84 ± 2.20	4.04 ± 2.01	4.65 ± 2.16	4.44 ± 2.11	4.01 ± 2.00	4.59 ± 2.14	3.73 ± 1.93	4.32 ± 2.08	4.03 ± 2.01	4.30 ± 2.07	3.70 ± 1.92	3.95 ± 1.99	4.08 ± 2.02	4.45 ± 2.11	3.47 ± 1.86
17 (382.5)	3.81 ± 1.95	2.98 ± 1.73	4.32 ± 2.08	5.60 ± 2.37	5.55 ± 2.36	4.41 ± 2.10	4.95 ± 2.22	4.95 ± 2.23	4.80 ± 2.19	4.98 ± 2.23	4.95 ± 2.22	4.43 ± 2.10	4.87 ± 2.21	5.40 ± 2.32	4.26 ± 2.06	4.97 ± 2.23	5.54 ± 2.35	4.44 ± 2.11	4.63 ± 2.15	5.69 ± 2.38
18 (405.0)	2.74 ± 1.65	4.12 ± 2.03	4.12 ± 2.03	5.94 ± 2.44	7.20 ± 2.68	5.54 ± 2.35	5.15 ± 2.27	5.91 ± 2.43	5.16 ± 2.27	6.10 ± 2.47	5.43 ± 2.33	5.79 ± 2.41	5.22 ± 2.28	5.06 ± 2.25	6.16 ± 2.48	5.26 ± 2.29	4.93 ± 2.22	5.51 ± 2.35	4.65 ± 2.16	5.35 ± 2.31
19 (427.5)	4.15 ± 2.04	4.64 ± 2.15	4.62 ± 2.15	5.08 ± 2.25	7.30 ± 2.70	7.56 ± 2.75	7.08 ± 2.66	5.28 ± 2.30	6.56 ± 2.56	6.78 ± 2.60	7.18 ± 2.68	6.61 ± 2.57	7.03 ± 2.65	7.07 ± 2.66	7.66 ± 2.77	6.17 ± 2.48	6.12 ± 2.47	7.12 ± 2.67	8.00 ± 2.83	5.53 ± 2.35
20 (450.0)	11.12 ± 3.33	5.11 ± 2.26	6.54 ± 2.56	5.57 ± 2.36	8.79 ± 2.97	6.48 ± 2.55	6.63 ± 2.58	7.41 ± 2.72	7.77 ± 2.79	7.02 ± 2.65	7.97 ± 2.82	8.70 ± 2.95	8.54 ± 2.92	8.91 ± 2.99	8.61 ± 2.93	6.24 ± 2.50	9.46 ± 3.08	8.86 ± 2.98	6.92 ± 2.63	8.15 ± 2.85

`PT & ETA`

η \ pt	1 (0.0,1250.0)	2 (2500.0)	3 (3750.0)	4 (5000.0)	5 (6250.0)	6 (7500.0)	7 (8750.0)	8 (10000.0)	9 (11250.0)	10 (12500.0)	11 (13750.0)	12 (15000.0)	13 (16250.0)	14 (17500.0)	15 (18750.0)	16 (20000.0)	17 (21250.0)	18 (22500.0)	19 (23750.0)	20 (25000.0)
1 (2.0,2.3)	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.57 ± 0.75	0.73 ± 0.85	0.72 ± 0.85	0.64 ± 0.80	0.57 ± 0.76	0.66 ± 0.81	0.65 ± 0.80	0.64 ± 0.80	0.59 ± 0.77	0.66 ± 0.81	0.62 ± 0.79	0.59 ± 0.77	0.66 ± 0.81	0.69 ± 0.83	0.65 ± 0.80	0.65 ± 0.80	0.69 ± 0.83
2 (2.7)	0.00 ± 0.00	0.00 ± 0.00	0.85 ± 0.92	0.92 ± 0.96	0.86 ± 0.93	0.78 ± 0.88	0.74 ± 0.86	0.69 ± 0.83	0.67 ± 0.82	0.67 ± 0.82	0.66 ± 0.81	0.66 ± 0.81	0.61 ± 0.78	0.60 ± 0.77	0.60 ± 0.78	0.60 ± 0.77	0.58 ± 0.76	0.61 ± 0.78	0.68 ± 0.83	0.75 ± 0.87
3 (3.0)	0.00 ± 0.00	0.96 ± 0.98	1.02 ± 1.01	1.00 ± 1.00	0.89 ± 0.94	0.80 ± 0.90	0.74 ± 0.86	0.71 ± 0.84	0.69 ± 0.83	0.69 ± 0.83	0.68 ± 0.82	0.65 ± 0.81	0.64 ± 0.80	0.71 ± 0.84	0.81 ± 0.90	0.89 ± 0.94	1.13 ± 1.06	1.27 ± 1.13	1.37 ± 1.17	1.62 ± 1.27
4 (3.3)	0.00 ± 0.00	1.05 ± 1.02	1.10 ± 1.05	1.05 ± 1.02	0.91 ± 0.95	0.83 ± 0.91	0.76 ± 0.87	0.74 ± 0.86	0.73 ± 0.86	0.73 ± 0.86	0.82 ± 0.91	1.00 ± 1.00	1.33 ± 1.15	1.63 ± 1.28	2.44 ± 1.56	3.30 ± 1.82	6.71 ± 2.59	11.44 ± 3.38	32.32 ± 5.68	75.92 ± 8.71
5 (3.7)	0.81 ± 0.90	1.15 ± 1.07	1.14 ± 1.07	1.04 ± 1.02	0.94 ± 0.97	0.82 ± 0.90	0.82 ± 0.91	0.85 ± 0.92	1.07 ± 1.04	1.41 ± 1.19	2.42 ± 1.56	3.91 ± 1.98	7.84 ± 2.80	27.46 ± 5.24	112.23 ± 10.59	356.42 ± 18.88	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
6 (4.0)	1.01 ± 1.01	1.17 ± 1.08	1.19 ± 1.09	1.09 ± 1.04	1.01 ± 1.00	1.06 ± 1.03	1.39 ± 1.18	2.53 ± 1.59	5.18 ± 2.28	14.83 ± 3.85	68.01 ± 8.25	658.25 ± 25.66	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
7 (4.3)	1.09 ± 1.05	1.21 ± 1.10	1.21 ± 1.10	1.27 ± 1.13	1.80 ± 1.34	3.65 ± 1.91	11.65 ± 3.41	91.71 ± 9.58	3985.49 ± 63.13	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
8 (4.7)	1.12 ± 1.06	1.28 ± 1.13	1.42 ± 1.19	2.55 ± 1.60	10.81 ± 3.29	93.97 ± 9.69	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00
9 (5.0)	1.20 ± 1.10	1.37 ± 1.17	2.47 ± 1.57	15.12 ± 3.89	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00	0.00 ± 0.00

yipengsun commented 2 years ago

Here's the ratio plots:

b_eta b_ownpv_ndof b_pt ntracks

yipengsun commented 2 years ago

I think in this case the plots are misleading: Because it only shows the 1D ratio, and all 4 variables follow a smooth pattern. But if we look at the 2D histograms, some of the bins (edge bins) have much larger efficiencies.

yipengsun commented 2 years ago

So, the PID efficiencies are always non-0, BUT there's a significant portion of events that falls OUTSIDE the binning range.

I'm going to first plot the P range for pidcalib samples, and see if we can meaningfully increase the max bin edge. If not, I'll just round it to nearest bin.

yipengsun commented 2 years ago

`P` distribution for `K` PIDCalib sample

Brunel_P

`P` distribution for `Mu_nopt` PIDCalib sample

Brunel_P

yipengsun commented 2 years ago

I think it is OK to put max P at 300 GeV.

yipengsun commented 2 years ago

`K` PID efficiency of a `ETA` slice:

At P > 200 GeV, PID efficiencies are close to 0. I think in this case these are real 0's because efficiency has VERY small errors.

η, 5 (5.0)

p \ nTrk	1 (0.0,50.0)	2 (200.0)	3 (300.0)	4 (500.0)
1 (3000.0,9300.0)	nan	0.24804 ± 0.24452	0.36279 ± 0.25997	1.59017 ± 1.68936
2 (15600.0)	0.93889 ± 0.02419	0.69755 ± 0.00448	0.56763 ± 0.00744	0.46221 ± 0.01439
3 (19600.0)	0.95854 ± 0.01589	0.79297 ± 0.00320	0.67831 ± 0.00556	0.56692 ± 0.01071
4 (24400.0)	0.91049 ± 0.01034	0.83317 ± 0.00223	0.76369 ± 0.00413	0.66820 ± 0.00786
5 (29800.0)	0.87638 ± 0.00925	0.83054 ± 0.00173	0.77938 ± 0.00309	0.73869 ± 0.00609
6 (35200.0)	0.85049 ± 0.00899	0.81743 ± 0.00155	0.78395 ± 0.00274	0.73292 ± 0.00529
7 (40600.0)	0.82855 ± 0.00896	0.79931 ± 0.00149	0.76522 ± 0.00260	0.73517 ± 0.00504
8 (46000.0)	0.79330 ± 0.00912	0.76308 ± 0.00149	0.72917 ± 0.00259	0.70257 ± 0.00480
9 (51400.0)	0.77349 ± 0.00985	0.71951 ± 0.00153	0.68364 ± 0.00260	0.65527 ± 0.00478
10 (56800.0)	0.70825 ± 0.01028	0.67437 ± 0.00157	0.63824 ± 0.00264	0.60432 ± 0.00476
11 (62200.0)	0.68610 ± 0.01044	0.62988 ± 0.00161	0.59027 ± 0.00269	0.54390 ± 0.00480
12 (67600.0)	0.64492 ± 0.01104	0.57379 ± 0.00166	0.54026 ± 0.00273	0.49717 ± 0.00471
13 (73000.0)	0.60126 ± 0.01157	0.51343 ± 0.00168	0.47615 ± 0.00275	0.44328 ± 0.00474
14 (78400.0)	0.52176 ± 0.01225	0.45504 ± 0.00171	0.41520 ± 0.00275	0.37498 ± 0.00467
15 (83800.0)	0.47010 ± 0.01275	0.39283 ± 0.00171	0.35633 ± 0.00273	0.32095 ± 0.00460
16 (89200.0)	0.38623 ± 0.01295	0.33116 ± 0.00171	0.29523 ± 0.00266	0.26431 ± 0.00450
17 (94600.0)	0.31218 ± 0.01257	0.27187 ± 0.00168	0.23785 ± 0.00259	0.20321 ± 0.00429
18 (100000.0)	0.26707 ± 0.01322	0.21794 ± 0.00163	0.19127 ± 0.00250	0.16014 ± 0.00409
19 (200000.0)	0.07291 ± 0.00252	0.05506 ± 0.00031	0.04599 ± 0.00047	0.04071 ± 0.00077
20 (250000.0)	0.00000 ± 0.00000	0.00000 ± 0.00000	0.00000 ± 0.00000	0.00000 ± 0.00000
21 (300000.0)	0.00000 ± 0.00000	0.00000 ± 0.00000	0.00000 ± 0.00000	0.00000 ± 0.00000

`Mu_nopt` PID efficiency of a `ETA` slice

η, 5 (5.0)

p \ nTrk	1 (0.0,50.0)	2 (200.0)	3 (300.0)	4 (500.0)
1 (3000.0,9300.0)	0.78962 ± 0.10483	0.90174 ± 0.05020	0.80972 ± 0.11653	0.77745 ± 0.21505
2 (15600.0)	0.89096 ± 0.04431	0.86826 ± 0.01371	0.85746 ± 0.03258	0.79838 ± 0.06667
3 (19600.0)	0.86254 ± 0.04477	0.80078 ± 0.01053	0.71774 ± 0.02147	0.63754 ± 0.03480
4 (24400.0)	0.82087 ± 0.03401	0.75763 ± 0.00703	0.68452 ± 0.01388	0.64443 ± 0.02355
5 (29800.0)	0.79997 ± 0.02437	0.71427 ± 0.00509	0.63695 ± 0.00961	0.56879 ± 0.01529
6 (35200.0)	0.78125 ± 0.01963	0.73460 ± 0.00436	0.65307 ± 0.00802	0.61616 ± 0.01385
7 (40600.0)	0.80668 ± 0.01681	0.74170 ± 0.00382	0.69004 ± 0.00733	0.64196 ± 0.01214
8 (46000.0)	0.80924 ± 0.01516	0.74536 ± 0.00352	0.70788 ± 0.00682	0.64501 ± 0.01070
9 (51400.0)	0.79482 ± 0.01450	0.76589 ± 0.00351	0.73182 ± 0.00664	0.69814 ± 0.01092
10 (56800.0)	0.76873 ± 0.01420	0.74432 ± 0.00338	0.70234 ± 0.00632	0.69185 ± 0.01051
11 (62200.0)	0.78313 ± 0.01474	0.72608 ± 0.00340	0.69489 ± 0.00630	0.67549 ± 0.00991
12 (67600.0)	0.76683 ± 0.01472	0.72622 ± 0.00349	0.70756 ± 0.00655	0.67262 ± 0.00996
13 (73000.0)	0.71764 ± 0.01537	0.70003 ± 0.00358	0.67894 ± 0.00640	0.67122 ± 0.01003
14 (78400.0)	0.72514 ± 0.01634	0.69370 ± 0.00376	0.68278 ± 0.00678	0.65875 ± 0.01036
15 (83800.0)	0.69980 ± 0.01804	0.67078 ± 0.00384	0.66128 ± 0.00699	0.63626 ± 0.01052
16 (89200.0)	0.66907 ± 0.01857	0.66616 ± 0.00410	0.64624 ± 0.00718	0.63366 ± 0.01094
17 (94600.0)	0.66213 ± 0.02061	0.64277 ± 0.00419	0.63145 ± 0.00746	0.62133 ± 0.01103
18 (100000.0)	0.62902 ± 0.02181	0.63461 ± 0.00446	0.62616 ± 0.00797	0.60542 ± 0.01163
19 (200000.0)	0.55328 ± 0.00823	0.55779 ± 0.00158	0.55757 ± 0.00271	0.55464 ± 0.00391
20 (250000.0)	0.50822 ± 0.02927	0.49160 ± 0.00464	0.50479 ± 0.00784	0.50343 ± 0.01130
21 (300000.0)	0.49919 ± 0.04745	0.48649 ± 0.00711	0.50672 ± 0.01198	0.47060 ± 0.01720

umd-lhcb / lhcb-ntuples-gen