1+ #include < benchmark/benchmark.h>
2+
3+ #include < fcntl.h>
4+ #include < sys/stat.h>
5+ #include < sys/types.h>
6+ #include < unistd.h>
7+
8+ #include < algorithm>
9+ #include < cassert>
10+ #include < cstdio>
11+ #include < iostream>
12+ #include < future>
13+ #include < memory>
14+ #include < string>
15+ #include < vector>
16+
17+ #include < ROOT/RDataFrame.hxx>
18+ #include < ROOT/RNTuple.hxx>
19+ #include < ROOT/RNTupleDS.hxx>
20+ #include < ROOT/RNTupleOptions.hxx>
21+ #include < Compression.h>
22+ #include < TApplication.h>
23+ #include < TBranch.h>
24+ #include < TCanvas.h>
25+ #include < TClassTable.h>
26+ #include < TF1.h>
27+ #include < TFile.h>
28+ #include < TH1D.h>
29+ #include < TLatex.h>
30+ #include < TLegend.h>
31+ #include < TLorentzVector.h>
32+ #include < TMath.h>
33+ #include < TStyle.h>
34+ #include < TSystem.h>
35+ #include < TTree.h>
36+ #include < TTreeReader.h>
37+ #include < TTreePerfStats.h>
38+
39+ #include < Math/Vector4D.h>
40+
41+ #include < rootbench/RBConfig.h>
42+
43+
44+ static float ComputeInvariantMass (float pt0, float pt1, float eta0, float eta1, float phi0, float phi1, float e0 , float e1 )
45+ {
46+ ROOT::Math::PtEtaPhiEVector p1 (pt0, eta0, phi0, e0 );
47+ ROOT::Math::PtEtaPhiEVector p2 (pt1, eta1, phi1, e1 );
48+ return (p1 + p2).mass () / 1000.0 ;
49+ }
50+
51+
52+ static void ProcessNTuple (ROOT::Experimental::RNTupleReader *ntuple, TH1D *hMass, bool isMC)
53+ {
54+
55+ auto viewTrigP = ntuple->GetView <bool >(" trigP" );
56+ auto viewPhotonN = ntuple->GetView <std::uint32_t >(" photon_n" );
57+ auto viewPhotonIsTightId = ntuple->GetView <std::vector<bool >>(" photon_isTightID" );
58+ auto viewPhotonPt = ntuple->GetView <std::vector<float >>(" photon_pt" );
59+ auto viewPhotonEta = ntuple->GetView <std::vector<float >>(" photon_eta" );
60+ auto viewPhotonPhi = ntuple->GetView <std::vector<float >>(" photon_phi" );
61+ auto viewPhotonE = ntuple->GetView <std::vector<float >>(" photon_E" );
62+ auto viewPhotonPtCone30 = ntuple->GetView <std::vector<float >>(" photon_ptcone30" );
63+ auto viewPhotonEtCone20 = ntuple->GetView <std::vector<float >>(" photon_etcone20" );
64+
65+ auto viewScaleFactorPhoton = ntuple->GetView <float >(" scaleFactor_PHOTON" );
66+ auto viewScaleFactorPhotonTrigger = ntuple->GetView <float >(" scaleFactor_PhotonTRIGGER" );
67+ auto viewScaleFactorPileUp = ntuple->GetView <float >(" scaleFactor_PILEUP" );
68+ auto viewMcWeight = ntuple->GetView <float >(" mcWeight" );
69+
70+ unsigned nevents = 0 ;
71+ for (auto e : ntuple->GetEntryRange ()) {
72+ nevents++;
73+
74+ if (!viewTrigP (e)) continue ;
75+
76+ std::vector<size_t > idxGood;
77+ auto isTightId = viewPhotonIsTightId (e);
78+ auto pt = viewPhotonPt (e);
79+ auto eta = viewPhotonEta (e);
80+
81+ for (size_t i = 0 ; i < viewPhotonN (e); ++i) {
82+ if (!isTightId[i]) continue ;
83+ if (pt[i] <= 25000 .) continue ;
84+ if (abs (eta[i]) >= 2.37 ) continue ;
85+ if (abs (eta[i]) >= 1.37 && abs (eta[i]) <= 1.52 ) continue ;
86+ idxGood.push_back (i);
87+ }
88+ if (idxGood.size () != 2 ) continue ;
89+
90+ auto ptCone30 = viewPhotonPtCone30 (e);
91+ auto etCone20 = viewPhotonEtCone20 (e);
92+
93+ bool isIsolatedPhotons = true ;
94+ for (int i = 0 ; i < 2 ; ++i) {
95+ if ((ptCone30[idxGood[i]] / pt[idxGood[i]] >= 0.065 ) ||
96+ (etCone20[idxGood[i]] / pt[idxGood[i]] >= 0.065 ))
97+ {
98+ isIsolatedPhotons = false ;
99+ break ;
100+ }
101+ }
102+ if (!isIsolatedPhotons) continue ;
103+
104+ auto phi = viewPhotonPhi (e);
105+ auto E = viewPhotonE (e);
106+
107+ float myy = ComputeInvariantMass (
108+ pt[idxGood[0 ]], pt[idxGood[1 ]],
109+ eta[idxGood[0 ]], eta[idxGood[1 ]],
110+ phi[idxGood[0 ]], phi[idxGood[1 ]],
111+ E[idxGood[0 ]], E[idxGood[1 ]]);
112+
113+ if (pt[idxGood[0 ]] / 1000 . / myy <= 0.35 ) continue ;
114+ if (pt[idxGood[1 ]] / 1000 . / myy <= 0.25 ) continue ;
115+ if (myy <= 105 ) continue ;
116+ if (myy >= 160 ) continue ;
117+
118+ if (isMC) {
119+ auto weight = viewScaleFactorPhoton (e) * viewScaleFactorPhotonTrigger (e) *
120+ viewScaleFactorPileUp (e) * viewMcWeight (e);
121+ hMass->Fill (myy, weight);
122+ } else {
123+ hMass->Fill (myy);
124+ }
125+ }
126+ }
127+
128+
129+ static void BM_RNTuple_ATLAS (benchmark::State &state, const std::string &comprAlgorithm)
130+ {
131+ using RNTupleReader = ROOT::Experimental::RNTupleReader;
132+ std::string path = RB::GetDataDir () + " /atlas-" + comprAlgorithm + " .ntuple" ;
133+
134+ for (auto _ : state) {
135+ auto hData = new TH1D (" " , " Diphoton invariant mass; m_{#gamma#gamma} [GeV];Events" , 30 , 105 , 160 );
136+
137+ auto ntuple = RNTupleReader::Open (" mini" , path);
138+ ProcessNTuple (ntuple.get (), hData, false /* isMC */ );
139+
140+ delete hData;
141+ }
142+ }
143+ BENCHMARK_CAPTURE (BM_RNTuple_ATLAS, BM_RNTuple_ATLAS_LZ4, " lz4" )->Unit(benchmark::kMicrosecond )->Iterations(5 );
144+ BENCHMARK_CAPTURE (BM_RNTuple_ATLAS, BM_RNTuple_ATLAS_ZLIB, " zlib" )->Unit(benchmark::kMicrosecond )->Iterations(5 );
145+ BENCHMARK_CAPTURE (BM_RNTuple_ATLAS, BM_RNTuple_ATLAS_LZMA, " lzma" )->Unit(benchmark::kMicrosecond )->Iterations(5 );
146+ BENCHMARK_CAPTURE (BM_RNTuple_ATLAS, BM_RNTuple_ATLAS_ZSTD, " zstd" )->Unit(benchmark::kMicrosecond )->Iterations(5 );
147+ BENCHMARK_CAPTURE (BM_RNTuple_ATLAS, BM_RNTuple_ATLAS_None, " none" )->Unit(benchmark::kMicrosecond )->Iterations(5 );
148+
149+ BENCHMARK_MAIN ();
0 commit comments