MetaDataTree.cxx 38.1 KB
Newer Older
1
2
3
#include <XAMPPbase/MetaDataTree.h>
#include <iostream>

Johannes Junggeburth's avatar
Johannes Junggeburth committed
4
5
6
#include <XAMPPbase/IAnalysisHelper.h>

#include <XAMPPbase/AnalysisUtils.h>
Johannes Junggeburth's avatar
Johannes Junggeburth committed
7
#include <XAMPPbase/EventInfo.h>
Johannes Junggeburth's avatar
Johannes Junggeburth committed
8
9
#include <xAODCutFlow/CutBookkeeper.h>
#include <xAODCutFlow/CutBookkeeperContainer.h>
10

Johannes Junggeburth's avatar
Johannes Junggeburth committed
11
#include <xAODMetaData/FileMetaData.h>
12
#include <xAODTruth/TruthMetaDataContainer.h>
Johannes Junggeburth's avatar
Johannes Junggeburth committed
13
14

//#include <PMGTools/PMGTruthWeightTool.h>
Johannes Junggeburth's avatar
Johannes Junggeburth committed
15
16
17
#include <xAODLuminosity/LumiBlockRange.h>
#include <xAODLuminosity/LumiBlockRangeContainer.h>

Johannes Junggeburth's avatar
Johannes Junggeburth committed
18
#ifndef XAOD_STANDALONE
Johannes Junggeburth's avatar
Johannes Junggeburth committed
19
20
21
#include <AthAnalysisBaseComps/AthAnalysisHelper.h>
#include <EventInfo/EventStreamInfo.h>
#include <GaudiKernel/ITHistSvc.h>
Johannes Junggeburth's avatar
Johannes Junggeburth committed
22
23
#endif
namespace XAMPP {
24
25
    // List of prcoess IDs
    // https://twiki.cern.ch/twiki/bin/view/AtlasProtected/SUSYSignalUncertainties#Subprocess_IDs
xampp's avatar
xampp committed
26
27
28
29
30
31
32
33
34
35
36
37
38
    std::vector<unsigned int> SUSYprocessIDs() {
        static std::vector<unsigned int> proc_ids;
        if (proc_ids.empty()) {
            for (int i = 1; i <= 220; ++i) {
                // Skip all ranges not occupied by any SUSYprocess
                if ((i > 4 && i < 51) || (i > 52 && i < 61) || (i > 62 && i < 71) || (i > 78 && i < 81) || (i > 89 && i < 111)) continue;
                if ((i > 118 && i < 122) || (i > 128 && i < 133) || (i > 138 && i < 144) || (i > 158 && i < 167)) continue;
                if ((i > 168 && i < 201) || i == 215) continue;
                proc_ids.push_back(i);
            }
        }
        return proc_ids;
    }
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    std::string susyProcName(int ID) {
        if (ID == 1)
            return "gluino-squark";
        else if (ID == 2)
            return "gluino-gluino";
        else if (ID == 3)
            return "squark-squark";
        else if (ID == 4)
            return "squark-antisquark";

        else if (ID == 51)
            return "sbottom1-antisbottom1";
        else if (ID == 51)
            return "sbottom2-antisbottom2";

        else if (ID == 61)
            return "stop1-antistop1";
        else if (ID == 62)
            return "stop2-antistop2";

        else if (ID == 71)
            return "gluino-nino1";
        else if (ID == 72)
            return "gluino-nino2";
        else if (ID == 73)
            return "gluino-nino3";
        else if (ID == 74)
            return "gluino-nino4";

        else if (ID == 75)
            return "gluino-chinoplus1";
        else if (ID == 76)
            return "gluino-chinoplus2";
        else if (ID == 77)
            return "gluino-chinominus1";
        else if (ID == 78)
            return "gluino-chinominus2";

        else if (ID == 81)
            return "squark-nino1";
        else if (ID == 82)
            return "squark-nino2";
        else if (ID == 83)
            return "squark-nino3";
        else if (ID == 84)
            return "squark-nino4";

        else if (ID == 85)
            return "squark-chinoplus1";
        else if (ID == 86)
            return "squark-chinoplus2";
        else if (ID == 87)
            return "squark-chinominus1";
        else if (ID == 88)
            return "squark-chinominus2";

        else if (ID == 111)
            return "nino1-nino1";
        else if (ID == 112)
            return "nino1-nino2";
        else if (ID == 113)
            return "nino1-nino3";
        else if (ID == 114)
            return "nino1-nino4";
        else if (ID == 115)
            return "nino1-chinoplus1";
        else if (ID == 116)
            return "nino1-chinoplus2";
        else if (ID == 117)
            return "nino1-chinminus1";
        else if (ID == 118)
            return "nino1-chinminus2";

        else if (ID == 122)
            return "nino2-nino2";
        else if (ID == 123)
            return "nino2-nino3";
        else if (ID == 124)
            return "nino2-nino4";
        else if (ID == 125)
            return "nino2-chinoplus1";
        else if (ID == 126)
            return "nino2-chinoplus2";
        else if (ID == 127)
            return "nino2-chinminus1";
        else if (ID == 128)
            return "nino2-chinminus2";

        else if (ID == 133)
            return "nino3-nino3";
        else if (ID == 134)
            return "nino3-nino4";
        else if (ID == 135)
            return "nino3-chinoplus1";
        else if (ID == 136)
            return "nino3-chinoplus2";
        else if (ID == 137)
            return "nino3-chinminus1";
        else if (ID == 138)
            return "nino3-chinminus2";

        else if (ID == 144)
            return "nino4-nino4";
        else if (ID == 145)
            return "nino4-chinoplus1";
        else if (ID == 146)
            return "nino4-chinoplus2";
        else if (ID == 147)
            return "nino4-chinminus1";
        else if (ID == 148)
            return "nino4-chinminus2";

        else if (ID == 157)
            return "chiplus1-chinminus1";
        else if (ID == 158)
            return "chiplus1-chinminus2";

        else if (ID == 167)
            return "chiplus2-chinminus1";
        else if (ID == 168)
            return "chiplus2-chinminus2";

        else if (ID == 201)
            return "selectronL-selectronL";
        else if (ID == 202)
            return "selectronR-selectronR";
        else if (ID == 203)
            return "sneutrinoEl-sneutrinoEl";

        else if (ID == 204)
            return "sneutrinoEl-selectronPlus";
        else if (ID == 205)
            return "sneutrinoEl-selectronMinus";

        else if (ID == 206)
            return "stau1-stau1";
        else if (ID == 207)
            return "stau2-stau2";
        else if (ID == 208)
            return "stau1-stau2";
        else if (ID == 209)
            return "sneutrinoTau-sneutrinoTau";

        else if (ID == 210)
            return "stau1Plus-sneutrinoTau";
        else if (ID == 211)
            return "stau1Minus-sneutrinoTau";
        else if (ID == 212)
            return "stau2Plus-sneutrinoTau";
        else if (ID == 213)
            return "stau2Minus-sneutrinoTau";

        else if (ID == 216)
            return "smuonL-smuonL";
        else if (ID == 217)
            return "smuonR-smuonR";
        else if (ID == 218)
            return "sneutrinoMu-sneutrinoMu";

        else if (ID == 219)
            return "sneutrinoMu-smuonPlus";
        else if (ID == 220)
            return "sneutrinoMu-smuonMinus";
202
203
        return "UNKOWN";
    }
xampp's avatar
xampp committed
204
205
206
    //#################################################################
    //                          MetaDataTree
    //#################################################################
207
    MetaDataTree::MetaDataTree(const std::string& myname) :
xampp's avatar
xampp committed
208
        AsgMetadataTool(myname),
209
        m_TreeName("MetaDataTree"),
xampp's avatar
xampp committed
210
211
        m_UseFileMetaData(true),
        m_fillLHEWeights(false),
212
213
        m_shiftMetaDSID(false),
        m_tree(nullptr),
xampp's avatar
xampp committed
214
215
216
217
218
219
220
        m_histSvc("THistSvc", myname),
        m_isData(false),
        m_init(false),
        m_analysis_helper("AnalysisHelper"),
        m_XAMPPInfo("EventInfoHandler"),
        m_MetaDB(),
        m_ActDB(m_MetaDB.end()),
Philipp Gadow's avatar
Philipp Gadow committed
221
        m_DummyEntry(nullptr) {
xampp's avatar
xampp committed
222
223
224
225
226
        declareProperty("isData", m_isData);
        declareProperty("TreeName", m_TreeName);
        declareProperty("AnalysisHelper", m_analysis_helper);
        declareProperty("useFileMetaData", m_UseFileMetaData);
        declareProperty("fillLHEWeights", m_fillLHEWeights);
227
        declareProperty("SwitchOnDSIDshift", m_shiftMetaDSID);
xampp's avatar
xampp committed
228
229
        m_XAMPPInfo.declarePropertyFor(this, "EventInfoHandler", "The XAMPPInfo event Handler");
    }
230
    MetaDataTree::~MetaDataTree() {}
xampp's avatar
xampp committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
    StatusCode MetaDataTree::initialize() {
        ATH_MSG_INFO("Initialising...");
        ATH_CHECK(m_analysis_helper.retrieve());
        ATH_CHECK(m_XAMPPInfo.retrieve());
        if (m_init) {
            ATH_MSG_WARNING("The metadata tree is already intialized");
            return StatusCode::SUCCESS;
        }
        m_tree = new TTree(m_TreeName.c_str(), "MetaData Tree for Small Analysis Ntuples");
        m_tree->Branch("isData", &m_isData);
        ATH_CHECK(m_histSvc->regTree("/XAMPP/" + m_TreeName, m_tree));
        m_init = true;
        if (m_isData) m_fillLHEWeights = false;
        return StatusCode::SUCCESS;
245
    }
xampp's avatar
xampp committed
246
    void MetaDataTree::LoadMCMetaData(unsigned int mcChannel, unsigned int periodNumber) {
247
        MetaID Ident(mcChannel, periodNumber);
248
        if (m_ActDB == m_MetaDB.end() || m_ActDB->first != Ident) {
249
            m_ActDB = m_MetaDB.find(Ident);
250
            if (m_ActDB == m_MetaDB.end()) {
251
                ATH_MSG_INFO("Found new mcChannel " << mcChannel << ".  Create new MCMetaData entry.");
xampp's avatar
xampp committed
252
                m_MetaDB.insert(std::pair<MetaID, std::shared_ptr<MetaDataElement>>(
253
                    Ident, std::make_shared<MetaDataMC>(mcChannel, periodNumber, m_analysis_helper, m_XAMPPInfo.getHandle(), this)));
254
                m_ActDB = m_MetaDB.find(Ident);
255
                if (m_fillLHEWeights) { m_ActDB->second->setLHEWeightNames(getLHEWeightNames()); }
256
            }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
257
258
        }
    }
259
260
261
    void MetaDataTree::LoadRunMetaData(unsigned int run) {
        MetaID Ident(run, run);
        if (m_ActDB == m_MetaDB.end() || m_ActDB->first != Ident) {
262
            m_ActDB = m_MetaDB.find(Ident);
263
            if (m_ActDB == m_MetaDB.end()) {
xampp's avatar
xampp committed
264
265
                ATH_MSG_INFO("Found new runNumber " << run << ".  Create new runMetaData entry.");
                m_MetaDB.insert(std::pair<MetaID, std::shared_ptr<MetaDataElement>>(
266
                    Ident, std::make_shared<runMetaData>(run, m_XAMPPInfo.getHandle())));
267
268
                m_ActDB = m_MetaDB.find(Ident);
            }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
269
        }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
270
    }
271
272
273
    StatusCode MetaDataTree::beginEvent() {
        ATH_CHECK(m_analysis_helper->LoadContainers());
        if (m_XAMPPInfo->isMC() == m_isData) {
274
            ATH_MSG_FATAL("The analysis is configured to run over data while the input file is MC...");
275
276
277
278
            return StatusCode::FAILURE;
        }
        // The current event is MC
        if (m_XAMPPInfo->isMC()) {
xampp's avatar
xampp committed
279
            LoadMCMetaData(m_XAMPPInfo->mcChannelNumber(), m_XAMPPInfo->runNumber());
280
281
282
283
        }
        // The event is data
        else
            LoadRunMetaData(m_XAMPPInfo->runNumber());
Johannes Junggeburth's avatar
Johannes Junggeburth committed
284

285
286
287
288
289
        if (m_DummyEntry) {
            ATH_CHECK(m_ActDB->second->CopyStore(m_DummyEntry.get()));
            m_DummyEntry = std::shared_ptr<MetaDataElement>();
        }
        ATH_CHECK(m_ActDB->second->newEvent());
290
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
291
    }
292
    StatusCode MetaDataTree::finalize() {
293
        if (!m_tree) { return StatusCode::SUCCESS; }
xampp's avatar
xampp committed
294
        for (const auto& Meta : m_MetaDB) ATH_CHECK(Meta.second->finalize(m_tree));
295
        m_tree = nullptr;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
296
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
297
    }
xampp's avatar
xampp committed
298
299
300
    StatusCode MetaDataTree::CheckLumiBlockContainer(const std::string& Container, bool& HasCont) {
        if (!inputMetaStore()->contains<xAOD::LumiBlockRangeContainer>(Container)) {
            ATH_MSG_DEBUG("Lumi block range container " << Container << " not present");
301
302
303
304
305
306
            return StatusCode::SUCCESS;
        }
        const xAOD::LumiBlockRangeContainer* LumiBlocks = nullptr;
        ATH_CHECK(inputMetaStore()->retrieve(LumiBlocks, Container));
        for (const auto& lumi : *LumiBlocks) {
            if (lumi->startRunNumber() != lumi->stopRunNumber())
xampp's avatar
xampp committed
307
308
309
                ATH_MSG_WARNING("Found different start (" << lumi->startRunNumber() << " and end (" << lumi->stopRunNumber()
                                                          << ") runNumbers for lumiblock range " << lumi->startRunNumber() << "-"
                                                          << lumi->stopLumiBlockNumber());
310
311
312
313
            LoadRunMetaData(lumi->startRunNumber());
            ATH_CHECK(m_ActDB->second->newFile(lumi));
        }
        HasCont = true;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
314
315
        return StatusCode::SUCCESS;
    }
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
    StatusCode MetaDataTree::fillLHEMetaData(unsigned int Idx) {
        if (Idx == 0) {
            ATH_MSG_ERROR(
                "This method is meant to store the variational weights in the "
                "metadata. Not nominal");
            return StatusCode::FAILURE;
        }
        if (m_isData) {
            ATH_MSG_ERROR("Data does not have any weight variations");
            return StatusCode::FAILURE;
        }
        if (!m_fillLHEWeights) {
            ATH_MSG_DEBUG("No lhe meta-data should be stored");
            return StatusCode::SUCCESS;
        }
        const xAOD::EventInfo* info = m_XAMPPInfo->GetOrigInfo();
332

333
334
335
336
        if (Idx >= info->mcEventWeights().size()) {
            ATH_MSG_ERROR("Index " << Idx << " is out of range");
            return StatusCode::FAILURE;
        }
xampp's avatar
xampp committed
337
        LoadMCMetaData(m_XAMPPInfo->mcChannelNumber(), m_XAMPPInfo->runNumber());
338
339
        double W = m_XAMPPInfo->GetGenWeight(Idx);
        unsigned int finalState = m_analysis_helper->finalState();
xampp's avatar
xampp committed
340
        if (finalState > 0 && !m_ActDB->second->fillVariation((Idx + 1000) * 1.e4 + finalState, W).isSuccess()) {
341
342
343
            return StatusCode::FAILURE;
        }
        return m_ActDB->second->fillVariation(Idx + 1000, W);
344
    }
345
    std::vector<std::string> MetaDataTree::getLHEWeightNames() const {
346
347
        std::vector<std::string> ret;
        if (m_isData) return ret;
348
349
        if (inputMetaStore()->contains<xAOD::TruthMetaDataContainer>("TruthMetaData")) {
            const xAOD::TruthMetaDataContainer* truthmetadata = nullptr;
350
351
352
            if (inputMetaStore()->retrieve(truthmetadata, "TruthMetaData").isFailure()) {
                ATH_MSG_WARNING("Failed to pick up truth metadata, can not write LHE by name");
                return ret;
353
354
            };
            for (auto meta : *truthmetadata) {
355
                if (!meta->weightNames().empty()) ret.insert(ret.end(), meta->weightNames().begin(), meta->weightNames().end());
356
357
            }
        }
358
        return ret;
359
360
    }

361
362
363
364
365
366
367
368
    StatusCode MetaDataTree::beginInputFile() {
        if (!m_UseFileMetaData) {
            ATH_MSG_INFO("File meta data is disabled");
            return StatusCode::SUCCESS;
        }
        bool HasLumiCont = false;
        if (m_isData) {
            ATH_CHECK(CheckLumiBlockContainer("LumiBlocks", HasLumiCont));
xampp's avatar
xampp committed
369
370
            ATH_CHECK(CheckLumiBlockContainer("IncompleteLumiBlocks", HasLumiCont));
            for (auto& Meta : m_MetaDB) Meta.second->CutBookKeeperAvailable(HasLumiCont);
371
        }
xampp's avatar
xampp committed
372
        if (inputMetaStore()->contains<xAOD::CutBookkeeperContainer>("CutBookkeepers")) {
373
374
375
376
377
378
            const xAOD::CutBookkeeperContainer* bks = nullptr;
            ATH_CHECK(inputMetaStore()->retrieve(bks, "CutBookkeepers"));
            if (!m_isData) {
                std::shared_ptr<MetaDataElement> MC;
                const EventStreamInfo* esi = nullptr;
                ATH_CHECK(inputMetaStore()->retrieve(esi));
xampp's avatar
xampp committed
379
                if (esi->getEventTypes().size() > 1) ATH_MSG_WARNING("There seem to be more event types than one");
380
381
382
383
384
                if (m_shiftMetaDSID && !m_XAMPPInfo->applyDSIDShift() && m_XAMPPInfo->dsidOffSet() != 0) {
                    ATH_MSG_WARNING(
                        "The meta-data information will be shifted by "
                        << m_XAMPPInfo->dsidOffSet()
                        << ". Please make sure that the information in the meta-data is consistent with the information given in the tree");
385
                }
386
387
                LoadMCMetaData(esi->getEventTypes().begin()->mc_channel_number() +
                                   (m_shiftMetaDSID || m_XAMPPInfo->applyDSIDShift() ? m_XAMPPInfo->dsidOffSet() : 0),
388
                               (*esi->getRunNumbers().begin()));
389
                MC = m_ActDB->second;
390

391
                if (m_fillLHEWeights) {
xampp's avatar
xampp committed
392
393
394
                    if (inputMetaStore()->contains<xAOD::TruthMetaDataContainer>("TruthMetaData")) {
                        const xAOD::TruthMetaDataContainer* truthmetadata = nullptr;
                        ATH_CHECK(inputMetaStore()->retrieve(truthmetadata, "TruthMetaData"));
395
                        for (auto meta : *truthmetadata) {
xampp's avatar
xampp committed
396
                            if (!meta->weightNames().empty()) MC->setLHEWeightNames(meta->weightNames());
397
                        }
398
                    }
399
400
                }
                ATH_CHECK(MC->newFile(bks));
Johannes Junggeburth's avatar
Johannes Junggeburth committed
401
            }
402
403
            // Data instance and no lumi container is present
            else if (!HasLumiCont) {
xampp's avatar
xampp committed
404
                if (m_DummyEntry) ATH_MSG_WARNING("There is still an instance of the dummy meta");
405
                m_DummyEntry = std::make_shared<runMetaData>(-1, m_XAMPPInfo.getHandle());
406
407
408
409
                ATH_CHECK(m_DummyEntry->newFile(bks));
            }
        } else {
            ATH_MSG_DEBUG("No CutBookKeeper has been found.");
xampp's avatar
xampp committed
410
            for (auto& Meta : m_MetaDB) Meta.second->CutBookKeeperAvailable(false);
411
412
413
        }
        return StatusCode::SUCCESS;
    }
414
    void MetaDataTree::subtractEventFromMetaData(unsigned int index) {
415
        if (m_isData) return;
xampp's avatar
xampp committed
416
        ATH_MSG_WARNING("Subtract event " << m_XAMPPInfo->eventNumber() << " in sample " << m_XAMPPInfo->mcChannelNumber()
417
                                          << " from the Sum of weights at index " << index);
xampp's avatar
xampp committed
418
        LoadMCMetaData(m_XAMPPInfo->mcChannelNumber(), m_XAMPPInfo->runNumber());
419
        unsigned int finalState = m_analysis_helper->finalState();
420
421
422
423
424
425
426
427
428
        unsigned int LHE_Idx = (index == 0 ? finalState : (finalState > 0 ? (index + 1000) * 1.e4 + finalState : index + 1000));
        m_ActDB->second->SubtractEvent(LHE_Idx, m_XAMPPInfo->GetRawGenWeight(index));
    }
    void MetaDataTree::subtractEventFromMetaData() {
        if (m_isData) return;
        ATH_MSG_WARNING("Subtract event " << m_XAMPPInfo->eventNumber() << " in sample " << m_XAMPPInfo->mcChannelNumber()
                                          << " from the Sum of weights.");
        const xAOD::EventInfo* info = m_XAMPPInfo->GetOrigInfo();

429
        size_t nWeights = m_fillLHEWeights ? info->mcEventWeights().size() : 1;
430
        for (size_t Idx = 0; Idx < nWeights; ++Idx) { subtractEventFromMetaData(Idx); }
431
    }
432
433
434
435
436
    std::string MetaDataTree::getLHEWeightNameForTree(const std::string& w_name) const {
        std::string n_name = EraseWhiteSpaces(w_name);
        n_name = ReplaceExpInString(n_name, " ", "_");
        n_name = ReplaceExpInString(n_name, "=", "_eq_");
        n_name = ReplaceExpInString(n_name, "/", "_by_");
437
438
        n_name = ReplaceExpInString(n_name, ".", "p");
        n_name = ReplaceExpInString(n_name, ",", "");
Johannes Junggeburth's avatar
Johannes Junggeburth committed
439
        n_name = ReplaceExpInString(n_name, ":", "");
440

441
442
        while (n_name.size() && n_name.rfind("_") == n_name.size() - 1) n_name = n_name.substr(0, n_name.size() - 1);
        return n_name;
443
    }
444
445
446
    //########################################################################################################
    //                                 MetaDataElement
    //########################################################################################################
xampp's avatar
xampp committed
447
448
449
450
451
    MetaDataElement::MetaDataElement(ToolHandle<XAMPP::IEventInfo> XAMPPInfo) : m_XAMPPInfo(XAMPPInfo), m_LHE_WeightNames() {}
    void MetaDataElement::setLHEWeightNames(const std::vector<std::string>& weights) { m_LHE_WeightNames = weights; }
    size_t MetaDataElement::numOfWeights() const { return m_LHE_WeightNames.size(); }
    const xAOD::CutBookkeeper* MetaDataElement::FindCutBookKeeper(const xAOD::CutBookkeeperContainer* container,
                                                                  MetaDataElement::BookKeeperType Type, unsigned int procID) {
452
453
454
455
        static bool lhe3MessageShown = false;
        const xAOD::CutBookkeeper* all = nullptr;
        int maxCycle = -1;  // need to find the max cycle where input stream is
                            // StreamAOD and the name is AllExecutedEvents
456

457
        std::vector<std::string> cb_names;
458
        if (Type == MetaDataElement::BookKeeperType::SUSY) {
459
            cb_names.push_back(Form("SUSYWeight_ID_%d", procID));
Johannes Junggeburth's avatar
Johannes Junggeburth committed
460
        } else if (Type == MetaDataElement::BookKeeperType::LHE3) {
461
            if (procID >= m_LHE_WeightNames.size()) {
xampp's avatar
xampp committed
462
                Error("MetaDataElement::FindCutBookKeeper()", "Invalid variation ID: %d", procID);
463
                return nullptr;
464
            }
465
            cb_names.push_back("LHE3Weight_" + RemoveAllExpInStr(m_LHE_WeightNames.at(procID), "."));
466
            cb_names.push_back(RemoveAllExpInStr(cb_names.back(), " "));
467
            cb_names.push_back(Form("AllExecutedEvents_NonNominalMCWeight_%d", procID));
468
            cb_names.push_back(m_LHE_WeightNames.at(procID));
469
470
471
        } else
            cb_names.push_back("AllExecutedEvents");

472
        for (auto cbk : *container) {
473
474
475
            bool Stream = (cbk->inputStream() == "StreamAOD" || cbk->inputStream().find("DAOD_TRUTH") != std::string::npos) &&
                          std::find(cb_names.begin(), cb_names.end(), cbk->name()) != cb_names.end();
            if (Stream && cbk->cycle() > maxCycle) {
476
477
478
                maxCycle = cbk->cycle();
                all = cbk;
            }
479
        }
480
        if (Type == MetaDataElement::BookKeeperType::LHE3 && !lhe3MessageShown) {
481
            Info("MetaDataElement::FindCutBookKeeper()", "Found cutbookkeeper matching LHE3 variation %s which will be associated to ID %u",
482
                 all->name().c_str(), 1000 + procID);
Johannes Junggeburth's avatar
Johannes Junggeburth committed
483
        }
xampp's avatar
xampp committed
484
        if (Type == MetaDataElement::BookKeeperType::LHE3) lhe3MessageShown = true;
485
        return all;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
486
    }
487

488
489
490
    //########################################################################################################
    //                                 MetaDataMC
    //########################################################################################################
491

xampp's avatar
xampp committed
492
    MetaDataMC::MetaDataMC(unsigned int mcChannel, unsigned int periodNumber, ToolHandle<XAMPP::IAnalysisHelper>& helper,
493
                           ToolHandle<XAMPP::IEventInfo> info, const IMetaDataTree* tree_class) :
494
        MetaDataElement(info),
495
496
497
        m_MC(mcChannel),
        m_periodNumber(periodNumber),
        m_helper(helper),
498
        m_tree_inst(tree_class),
499
500
501
502
        m_Data(),
        m_ActMeta(m_Data.end()),
        m_Inclusive(0),
        m_init(false) {
xampp's avatar
xampp committed
503
        m_init = m_helper.retrieve().isSuccess() && m_XAMPPInfo.retrieve().isSuccess();
504
505
506
507
508
        if (m_init) LoadMetaData(0);
    }
    StatusCode MetaDataMC::newFile(const xAOD::LumiBlockRange*) {
        Warning("MetaDataMC::newFile()", "Lumi blocks not supported.");
        return StatusCode::SUCCESS;
509
    }
510
511
512
513
514
515
    void MetaDataMC::LoadMetaData(unsigned int ID) {
        // See whether the process ID is already Loaded;
        if (m_ActMeta == m_Data.end() || m_ActMeta->first != ID) {
            m_ActMeta = m_Data.find(ID);
            // No Meta data has been found
            if (m_ActMeta == m_Data.end()) {
516
                std::string procName = "";
517
                if (ID > 1000 && (ID - 1000) < m_LHE_WeightNames.size()) {
518
                    procName = m_tree_inst->getLHEWeightNameForTree(m_LHE_WeightNames.at(ID - 1000));
xampp's avatar
xampp committed
519
520
                    Info("MetaDataMC::LoadMetaData()", Form("New LHE variation %u found for mcChannelNumber %u "
                                                            "corresponding to %s",
521
                                                            ID, m_MC, procName.c_str()));
522
523

                } else {
xampp's avatar
xampp committed
524
                    Info("MetaDataMC::LoadMetaData()", Form("New process Id %u found for mcChannelNumber %u", ID, m_MC));
525
                    Info("MetaDataMC::LoadMetaData()", Form("have %lu names ", m_LHE_WeightNames.size()));
526
527
528
                    if (ID == 0) {
                        procName = "Nominal_Inclusive";
                    } else {
529
                        procName = susyProcName(ID);
530
531
                    }
                }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
532

533
534
                m_Data.insert(
                    std::pair<unsigned int, std::shared_ptr<MetaDataMC::MetaData>>(ID, std::make_shared<MetaDataMC::MetaData>((ID))));
535
                m_ActMeta = m_Data.find(ID);
536
                m_ActMeta->second->procName = m_tree_inst->getLHEWeightNameForTree(procName);
537
538
539
                // If the SUSY process id = 0 -> inclusive state
                if (ID == 0) m_Inclusive = m_ActMeta->second;
            }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
540
        }
541
    }
542

xampp's avatar
xampp committed
543
    StatusCode MetaDataMC::newFile(const xAOD::CutBookkeeperContainer* container) {
544
        if (!m_init) {
xampp's avatar
xampp committed
545
            Error("MetaDataMC::newFile()", "The component has not been intialized");
546
547
548
            return StatusCode::FAILURE;
        }
        if (!container) {
xampp's avatar
xampp committed
549
            Error("MetaDataMC::newFile()", "No CutBookkeeperContainer was given");
550
551
552
553
554
            return StatusCode::FAILURE;
        }
        CutBookKeeperAvailable(false);
        const xAOD::CutBookkeeper* all = FindCutBookKeeper(container);
        if (!all) {
xampp's avatar
xampp committed
555
            Error("MetaDataMC::newFile()", "Could not read out the CutBookKeeper");
556
557
            return StatusCode::FAILURE;
        }
558
        std::function<double(const xAOD::CutBookkeeper*)> sumW = [this](const xAOD::CutBookkeeper* cbk) {
559
560
561
562
            return m_XAMPPInfo->getOutlierWeightStrategy() != IEventInfo::resetWeight
                       ? cbk->sumOfEventWeights()
                       : cbk->sumOfEventWeights();  // MG: The truncated version is broken as of 19-03-18
                                                    //   : cbk->sumOfTruncatedEventWeights();
563
564
        };
        std::function<double(const xAOD::CutBookkeeper*)> sumW2 = [this](const xAOD::CutBookkeeper* cbk) {
565
566
567
568
            return m_XAMPPInfo->getOutlierWeightStrategy() != IEventInfo::resetWeight
                       ? cbk->sumOfEventWeightsSquared()
                       : cbk->sumOfEventWeightsSquared();  // MG: The truncated version is broken as of 19-03-18
                                                           //   : cbk->sumOfTruncatedEventWeightsSquared();
569
570
571
        };

        AddFileInformation(m_Inclusive, all->nAcceptedEvents(), sumW(all), sumW2(all));
572
573
        // Load the SUSY Cut BookKeepers
        for (auto& ID : SUSYprocessIDs()) {
xampp's avatar
xampp committed
574
            const xAOD::CutBookkeeper* proc = FindCutBookKeeper(container, MetaDataElement::BookKeeperType::SUSY, ID);
575
576
            if (proc == nullptr || !proc->sumOfEventWeights()) continue;
            LoadMetaData(ID);
577
            AddFileInformation(m_ActMeta->second, proc->nAcceptedEvents(), sumW(proc), sumW2(proc));
578
579
580
581
        }
        // Load the LHE3 CutBookKeepers
        if (numOfWeights() > 0) {
            for (size_t lhe = numOfWeights() - 1; lhe > 0; --lhe) {
xampp's avatar
xampp committed
582
583
                const xAOD::CutBookkeeper* lhe_keeper = FindCutBookKeeper(container, MetaDataElement::BookKeeperType::LHE3, lhe);
                if (lhe_keeper == nullptr || !lhe_keeper->sumOfEventWeights()) continue;
584
                LoadMetaData(1000 + lhe);
585
                AddFileInformation(m_ActMeta->second, lhe_keeper->nAcceptedEvents(), sumW(lhe_keeper), sumW2(lhe_keeper));
586
587
            }
        }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
588

589
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
590
    }
591
592
    StatusCode MetaDataMC::newEvent() {
        if (!m_init) {
xampp's avatar
xampp committed
593
            Error("MetaDataMC::newFile()", "The component has not been inialized");
594
595
596
597
            return StatusCode::FAILURE;
        }
        unsigned int ProcID = m_helper->finalState();
        double Weight = m_XAMPPInfo->GetGenWeight();
Johannes Junggeburth's avatar
Johannes Junggeburth committed
598

599
600
601
602
603
604
        if (ProcID > 0) {
            LoadMetaData(ProcID);
            AddEventInformation(m_ActMeta->second, Weight);
        }
        AddEventInformation(m_Inclusive, Weight);
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
605
    }
606
607
608
609
610
611
    void MetaDataMC::SubtractEvent(unsigned int ProcID, double W) {
        if (ProcID > 0) {
            LoadMetaData(ProcID);
            SubtractEvent(m_ActMeta->second, W);
        }
        if (ProcID <= 1000) SubtractEvent(m_Inclusive, W);
Johannes Junggeburth's avatar
Johannes Junggeburth committed
612
    }
613
614
    StatusCode MetaDataMC::fillVariation(unsigned int Id, double W) {
        if (Id <= 1000) {
xampp's avatar
xampp committed
615
            Error("MetaDataMC::fillVariation", Form("Process id %u is protected.", Id));
616
617
618
619
620
            return StatusCode::FAILURE;
        }
        LoadMetaData(Id);
        AddEventInformation(m_ActMeta->second, W);
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
621
    }
622

xampp's avatar
xampp committed
623
    void MetaDataMC::AddEventInformation(std::shared_ptr<MetaDataMC::MetaData> Meta, double GenWeight) {
624
625
626
        ++Meta->NumProcessedEvents;
        if (!Meta->MetaInit) {
            // Retrieve the processId
xampp's avatar
xampp committed
627
            unsigned int SUSYId = Meta->ProcID != 0 ? m_helper->finalState() : 0;
628
            Meta->xSec = m_helper->GetMCXsec(m_MC, SUSYId);
629
630
631
632
633
634
635
636
            double rel_var_xs_down = -1;
            double rel_var_xs_up = -1;
            m_helper->GetMCXsecErrors(Meta->has_xSec_err, rel_var_xs_down, rel_var_xs_up, m_MC, SUSYId);
            // wait with setting the meta since we promise to ignore the two double refs if there is no well defined XS err
            if (Meta->has_xSec_err) {
                Meta->xSec_err_down = rel_var_xs_down;
                Meta->xSec_err_up = rel_var_xs_up;
            }
637
638
            Meta->FilterEff = m_helper->GetMCFilterEff(m_MC, SUSYId);
            Meta->kFaktor = m_helper->GetMCkFactor(m_MC, SUSYId);
639
            Meta->luminosity = m_XAMPPInfo->GetPileUpLuminosity();
640

641
642
            Meta->MetaInit = true;
            double xSecTimes = Meta->xSec * Meta->kFaktor * Meta->FilterEff;
xampp's avatar
xampp committed
643
644
            Info("MetaDataMC::AddEventInformation()",
                 Form("Set xSection to %f for process %u with mcChannelNumber %u", xSecTimes, Meta->ProcID, m_MC));
645
646
647
648
649
650
        }
        if (Meta->KeeperAvailable) return;
        Meta->SumW = Meta->SumW + GenWeight;
        Meta->SumW2 = Meta->SumW2 + GenWeight * GenWeight;
        ++Meta->NumTotalEvents;
    }
xampp's avatar
xampp committed
651
    void MetaDataMC::SubtractEvent(std::shared_ptr<MetaDataMC::MetaData> Meta, double GenWeight) {
652
653
654
655
656
657
658
659
        Meta->SumW = Meta->SumW - GenWeight;
        Meta->SumW2 = Meta->SumW2 - GenWeight * GenWeight;
    }
    void MetaDataMC::CutBookKeeperAvailable(bool B) {
        for (auto& MC : m_Data) MC.second->KeeperAvailable = B;
    }
    StatusCode MetaDataMC::finalize(TTree* MetaDataTree) {
        if (!m_init) {
xampp's avatar
xampp committed
660
            Error("MetaDataMC::finalize()", "The component has not been inialized");
661
662
663
            return StatusCode::FAILURE;
        }
        if (!SetBranchAddress(MetaDataTree, "mcChannelNumber", m_MC)) {
xampp's avatar
xampp committed
664
            Error("MetaDataMC::finalize()", Form("Cannot finalize DSID %u", m_MC));
665
666
667
            return StatusCode::FAILURE;
        }
        if (!SetBranchAddress(MetaDataTree, "runNumber", m_periodNumber)) {
xampp's avatar
xampp committed
668
            Error("MetaDataMC::finalize()", Form("Cannot finalize DSID %u", m_MC));
Johannes Junggeburth's avatar
Johannes Junggeburth committed
669
670
            return StatusCode::FAILURE;
        }
671
672
        for (auto& MC : m_Data) {
            if (!SaveMetaDataInTree(MC.second, MetaDataTree)) {
xampp's avatar
xampp committed
673
                Error("MetaDataMC::finalize()", Form("Cannot finalize DSID %u", m_MC));
674
675
676
677
678
                return StatusCode::FAILURE;
            }
        }
        m_init = false;
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
679
    }
xampp's avatar
xampp committed
680
    bool MetaDataMC::SaveMetaDataInTree(std::shared_ptr<MetaDataMC::MetaData> Meta, TTree* tree) {
681
        if (!Meta) {
xampp's avatar
xampp committed
682
            Error("MetaDataMC::SaveMetaDataInTree()", "No element has been given");
683
684
            return false;
        }
685
686
687
        /// Somehow ROOT has decided not to remap the Basket anymore
        /// thus lets use some static variables
        static std::string my_stonjek;
688
689
        static TTree* not_much_sense = nullptr;
        if (not_much_sense != tree) {
690
691
692
693
            if (!SetBranchAddress(tree, "ProcessName", &my_stonjek)) return false;
            not_much_sense = tree;
        }
        my_stonjek = Meta->procName;
694

695
696
        if (!SetBranchAddress(tree, "ProcessID", Meta->ProcID)) return false;
        if (!SetBranchAddress(tree, "xSection", Meta->xSec)) return false;
697
698
699
        if (!SetBranchAddress(tree, "xSection_RelErrorDown", Meta->xSec_err_down)) return false;
        if (!SetBranchAddress(tree, "xSection_RelErrorUp", Meta->xSec_err_up)) return false;
        if (!SetBranchAddress(tree, "xSection_has_Uncertainties", Meta->has_xSec_err)) return false;
700
        if (!SetBranchAddress(tree, "kFactor", Meta->kFaktor)) return false;
xampp's avatar
xampp committed
701
702
        if (!SetBranchAddress(tree, "FilterEfficiency", Meta->FilterEff)) return false;
        if (!SetBranchAddress(tree, "TotalEvents", Meta->NumTotalEvents)) return false;
703
704
        if (!SetBranchAddress(tree, "TotalSumW", Meta->SumW)) return false;
        if (!SetBranchAddress(tree, "TotalSumW2", Meta->SumW2)) return false;
xampp's avatar
xampp committed
705
        if (!SetBranchAddress(tree, "ProcessedEvents", Meta->NumProcessedEvents)) return false;
706
        // Luminosity of the data to which the dataset is going to be assigned
707
        if (m_XAMPPInfo->ApplyPileUp()) {
708
            if (!SetBranchAddress(tree, "prwLuminosity", Meta->luminosity)) return false;
709
        }
710
        // Check that the weights in the meta-data tree are finite
711
712
        std::function<bool(const double&)> is_finite = [](const double& V) { return !std::isnan(V) && !std::isinf(V); };
        if (!is_finite(Meta->NumTotalEvents)) {
713
714
715
            Error("MetaDataMC()", "Total events is not a number");
            return false;
        }
716
        if (!is_finite(Meta->SumW)) {
717
718
719
            Error("MetaDataMC()", "The SumW is not a number.");
            return false;
        }
720
        if (!is_finite(Meta->SumW2)) {
721
722
            Error("MetaDataMC()", "The SumW2 is not a number.");
            return false;
723
        }
724
725
        tree->Fill();
        Info("MetaDataMC::SaveMetaDataInTree()",
726
727
728
             Form("Write metadata in file with DSID: %u, ProcessID : %u, process name: %s SumW: %f, SumW2: %f, TotalEvents: %llu, "
                  "ProcessedEvents: %llu",
                  m_MC, Meta->ProcID, Meta->procName.c_str(), Meta->SumW, Meta->SumW2, Meta->NumTotalEvents, Meta->NumProcessedEvents));
729
730
        return true;
    }
xampp's avatar
xampp committed
731
    void MetaDataMC::AddFileInformation(std::shared_ptr<MetaDataMC::MetaData> Meta, Long64_t TotEv, double SumW, double SumW2) {
732
733
734
735
        Meta->NumTotalEvents = Meta->NumTotalEvents + TotEv;
        Meta->SumW = Meta->SumW + SumW;
        Meta->SumW2 = Meta->SumW2 + SumW2;
        Meta->KeeperAvailable = true;
736
        if (Meta->ProcID > 1000 && Meta->ProcID < 1000 + m_LHE_WeightNames.size()) {
737
            Meta->procName = m_LHE_WeightNames.at(Meta->ProcID - 1000);
738
        } else if (Meta->ProcID == 0) {
739
            Meta->procName = "Nominal_Inclusive";
740
        } else {
741
            Meta->procName = susyProcName(Meta->ProcID);
742
        }
743
744
745
746
747
748
749
750
751
752
753
754
755
    }
    StatusCode MetaDataMC::CopyStore(const MetaDataElement* Store) {
        const MetaDataMC* Other = dynamic_cast<const MetaDataMC*>(Store);
        if (!Store) {
            Error("MetaDataMC::CopyStore()", "Could not copy the input store");
            return StatusCode::FAILURE;
        }
        if (Other->m_MC != m_MC && Other->m_MC != (unsigned int)-1) {
            Error("MetaDataMC::CopyStore()", "Wrong MC channel");
            return StatusCode::FAILURE;
        }
        for (auto& toCopy : Other->m_Data) {
            LoadMetaData(toCopy.second->ProcID);
xampp's avatar
xampp committed
756
            AddFileInformation(m_ActMeta->second, toCopy.second->NumTotalEvents, toCopy.second->SumW, toCopy.second->SumW2);
757
            m_ActMeta->second->procName = toCopy.second->procName;
758
759
        }
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
760
    }
761
762
763
764

    //########################################################################################################
    //                                 runMetaData
    //########################################################################################################
xampp's avatar
xampp committed
765
    runMetaData::runMetaData(unsigned int runNumber, ToolHandle<XAMPP::IEventInfo> info) :
766
        MetaDataElement(info),
767
768
769
770
771
772
        m_runNumber(runNumber),
        m_NumTotalEvents(0),
        m_NumProcessedEvents(0),
        m_ProcessedBlocks(),
        m_TotalBlocks(),
        m_KeeperAvailable(false) {}
xampp's avatar
xampp committed
773
    StatusCode runMetaData::newFile(const xAOD::CutBookkeeperContainer* container) {
774
775
        const xAOD::CutBookkeeper* all = FindCutBookKeeper(container);
        if (!all) {
xampp's avatar
xampp committed
776
            Error("MetaDataMC::newFile()", "Could not read out the CutBookKeeper");
777
778
779
780
781
            return StatusCode::FAILURE;
        }
        CutBookKeeperAvailable(true);
        m_NumTotalEvents = m_NumTotalEvents + all->nAcceptedEvents();
        return StatusCode::SUCCESS;
782
    }
783
784
785
786
787
788
789
790
    StatusCode runMetaData::newEvent() {
        ++m_NumProcessedEvents;
        m_ProcessedBlocks.insert(m_XAMPPInfo->GetOrigInfo()->lumiBlock());
        if (!m_KeeperAvailable) {
            ++m_NumTotalEvents;
            m_TotalBlocks.insert(m_XAMPPInfo->GetOrigInfo()->lumiBlock());
        }
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
791
    }
792
793
    void runMetaData::CutBookKeeperAvailable(bool B) { m_KeeperAvailable = B; }
    StatusCode runMetaData::finalize(TTree* MetaDataTree) {
xampp's avatar
xampp committed
794
795
796
        if (!SetBranchAddress(MetaDataTree, "runNumber", m_runNumber)) return StatusCode::FAILURE;
        if (!SetBranchAddress(MetaDataTree, "TotalEvents", m_NumTotalEvents)) return StatusCode::FAILURE;
        if (!SetBranchAddress(MetaDataTree, "ProcessedEvents", m_NumProcessedEvents)) return StatusCode::FAILURE;
797
        std::set<unsigned int>* ProcLumiPtr = &m_ProcessedBlocks;
xampp's avatar
xampp committed
798
        if (!SetBranchAddress(MetaDataTree, "ProcessedLumiBlocks", ProcLumiPtr)) return StatusCode::FAILURE;
799
        std::set<unsigned int>* TotalLumiPtr = &m_TotalBlocks;
xampp's avatar
xampp committed
800
        if (!SetBranchAddress(MetaDataTree, "TotalLumiBlocks", TotalLumiPtr)) return StatusCode::FAILURE;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
801

xampp's avatar
xampp committed
802
803
804
        Info("runMetaData::finalize()", Form("Fill metadata tree with runNumber: %u, TotalEvents: %llu, "
                                             "ProcessedEvents: %llu",
                                             m_runNumber, m_NumTotalEvents, m_NumProcessedEvents));
805
806
        MetaDataTree->Fill();
        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
807
    }
808
809
810
811
812
813
    StatusCode runMetaData::CopyStore(const MetaDataElement* Store) {
        const runMetaData* toCopy = dynamic_cast<const runMetaData*>(Store);
        if (!toCopy) {
            Error("runMetaData::CopyStore()", "Failed to copy information");
            return StatusCode::FAILURE;
        }
xampp's avatar
xampp committed
814
        if (toCopy->m_runNumber != m_runNumber && toCopy->m_runNumber != (unsigned int)-1) {
815
816
817
818
819
            Error("runMetaData::CopyStore()", "Wrong run to copy");
            return StatusCode::FAILURE;
        }
        m_NumProcessedEvents += toCopy->m_NumProcessedEvents;
        m_NumTotalEvents += toCopy->m_NumTotalEvents;
xampp's avatar
xampp committed
820
        for (auto& BCID : toCopy->m_ProcessedBlocks) m_ProcessedBlocks.insert(BCID);
821
822
        for (auto& BCID : toCopy->m_TotalBlocks) m_TotalBlocks.insert(BCID);
        return StatusCode::SUCCESS;
823
    }
824
825
826
827
828
829
    StatusCode runMetaData::newFile(const xAOD::LumiBlockRange* LumiBlock) {
        if (!LumiBlock) {
            Error("runMetaData::newFile()", "No LumiBlock element given");
            return StatusCode::FAILURE;
        }
        m_NumTotalEvents = m_NumTotalEvents + LumiBlock->eventsSeen();
xampp's avatar
xampp committed
830
        for (unsigned int L = LumiBlock->startLumiBlockNumber(); L <= LumiBlock->stopLumiBlockNumber(); ++L) m_TotalBlocks.insert(L);
831
832

        return StatusCode::SUCCESS;
Johannes Junggeburth's avatar
Johannes Junggeburth committed
833
    }
xampp's avatar
xampp committed
834
    StatusCode runMetaData::fillVariation(unsigned int, double) { return StatusCode::FAILURE; }
Johannes Junggeburth's avatar
Johannes Junggeburth committed
835
836

}  // namespace XAMPP