Tool_FeatureExtractor.cxx 74.7 KB
Newer Older
1
2
3
4
/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

5

6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
///////////////////////////////////////////////////////////////////
//   Implementation file for class Tool_FeatureExtractor
///////////////////////////////////////////////////////////////////
// (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
// Tool to extract jet features from seed
///////////////////////////////////////////////////////////////////
// sebastian.fleischmann@cern.ch
///////////////////////////////////////////////////////////////////

//! Helper classes
#include "xAODTau/TauJet.h"
#include "xAODTracking/Vertex.h"
#include "xAODTracking/TrackParticle.h"

21
22
//#include "TrkParameters/TrackParameters.h"
//#include "TrkVertexFitterInterfaces/ITrackToVertexIPEstimator.h"
23
24
25
26
27
28
29
30
31
32
33
34
35
36

//! ROOT includes
#include "TMath.h"
#include "TLorentzVector.h"
#include "TVector3.h"

//! C++ includes
#include <vector>
#include <map>
#include <math.h>

//! PanTau includes
#include "PanTauAlgs/Tool_FeatureExtractor.h"
#include "PanTauAlgs/Tool_InformationStore.h"
37
38
39
#include "PanTauAlgs/TauConstituent.h"
#include "PanTauAlgs/PanTauSeed.h"
#include "PanTauAlgs/TauFeature.h"
40
#include "PanTauAlgs/HelperFunctions.h"
41
42
43
44




45
bool        sortTauConstituentMVA(const PanTau::TauConstituent2* u, const PanTau::TauConstituent2* v) {
46
47
48
49
50
51
52
    double uBDT = u->getBDTValue();
    double vBDT = v->getBDTValue();
    return uBDT > vBDT;
}



53
54
55
bool        sortTauConstituentEt(const PanTau::TauConstituent2* u, const PanTau::TauConstituent2* v) {
    double uEt = u->p4().Et();
    double vEt = v->p4().Et();
56
57
58
59
60
61
62
63
    return uEt > vEt;
}





PanTau::Tool_FeatureExtractor::Tool_FeatureExtractor(
64
65
66
67
68
69
    const std::string& name ) :
        asg::AsgTool(name),
        m_Tool_InformationStore("PanTau::Tool_InformationStore/Tool_InformationStore"){
	//m_trackToVertexTool("Reco::TrackToVertex") {
    
    //declareProperty("TrackToVertexTool", m_trackToVertexTool);
70
    declareProperty("Tool_InformationStore",            m_Tool_InformationStore,            "Tool handle to the information store tool");
71
    declareProperty("Tool_InformationStoreName",        m_Tool_InformationStoreName,            "Tool handle to the information store tool");
72
73
74
75
76
77
78
79
    
}



StatusCode PanTau::Tool_FeatureExtractor::initialize() {

    ATH_MSG_INFO(" initialize()");
80
    m_init=true;
81
    
82
    //ATH_CHECK( m_trackToVertexTool.retrieve() );
83
    ATH_CHECK( HelperFunctions::bindToolHandle( m_Tool_InformationStore, m_Tool_InformationStoreName ) );
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
    ATH_CHECK( m_Tool_InformationStore.retrieve() );
    
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Sum",          m_varTypeName_Sum) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Ratio",        m_varTypeName_Ratio) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_EtInRing",     m_varTypeName_EtInRing) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Isolation",    m_varTypeName_Isolation) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Num",          m_varTypeName_Num) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Mean",         m_varTypeName_Mean) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_StdDev",       m_varTypeName_StdDev) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_HLV",          m_varTypeName_HLV) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Angle",        m_varTypeName_Angle) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_DeltaR",       m_varTypeName_DeltaR) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_JetMoment",    m_varTypeName_JetMoment) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Combined",     m_varTypeName_Combined) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_JetShape",     m_varTypeName_JetShape) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_ImpactParams", m_varTypeName_ImpactParams) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Basic",        m_varTypeName_Basic) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_PID",          m_varTypeName_PID) );
    
    ATH_CHECK( m_Tool_InformationStore->getInfo_Int("FeatureExtractor_UseEmptySeeds",                           m_Config_UseEmptySeeds) );
    
    ATH_CHECK( m_Tool_InformationStore->getInfo_VecDouble("CellBased_BinEdges_Eta",               m_Config_CellBased_BinEdges_Eta) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_VecDouble("CellBased_EtaBinned_Pi0MVACut_1prong", m_Config_CellBased_EtaBinned_Pi0MVACut_1prong) );
    ATH_CHECK( m_Tool_InformationStore->getInfo_VecDouble("CellBased_EtaBinned_Pi0MVACut_3prong", m_Config_CellBased_EtaBinned_Pi0MVACut_3prong) );
    
    return StatusCode::SUCCESS;
110
111
112
113
}



114
void PanTau::Tool_FeatureExtractor::fillVariantsSeedEt(std::vector<PanTau::TauConstituent2*> tauConstituents) {
115
116
117
118
119
120
121
122
123
124

    //use different approaches to calculate total energy of seed:
    m_Variants_SeedEt["EtAllConsts"] = 0.0;
    m_Variants_SeedEt["EtNeutLowA"]  = 0.0;
    m_Variants_SeedEt["EtNeutLowB"]  = 0.0;
    
    //loop over all constituents in seed
    for(unsigned int iConst = 0; iConst < tauConstituents.size(); iConst++) {
        
        //get current constituents
125
126
        PanTau::TauConstituent2* curConstituent  = tauConstituents.at(iConst);
        double                  curEt           = curConstituent->p4().Et();
127
128
        
        //update the different Et definitions
129
        if(curConstituent->isOfType(PanTau::TauConstituent2::t_Charged) == true) {
130
131
132
133
            m_Variants_SeedEt["EtAllConsts"]    += curEt;
            m_Variants_SeedEt["EtNeutLowA"]     += curEt;
            m_Variants_SeedEt["EtNeutLowB"]     += curEt;
        }
134
        if(curConstituent->isOfType(PanTau::TauConstituent2::t_Neutral) == true) {
135
136
137
            m_Variants_SeedEt["EtAllConsts"]    += curEt;
        }
        
138
        if(curConstituent->isOfType(PanTau::TauConstituent2::t_NeutLowA) == true) {
139
140
            m_Variants_SeedEt["EtNeutLowA"]     += curEt;
        }
141
        if(curConstituent->isOfType(PanTau::TauConstituent2::t_NeutLowB) == true) {
142
143
144
145
146
147
148
149
150
151
            m_Variants_SeedEt["EtNeutLowB"]     += curEt;
        }
        
    }//end loop over constituents in seed
    
    return;
}



152
void    PanTau::Tool_FeatureExtractor::addFeatureWrtSeedEnergy( PanTau::TauFeature2* targetMap,
153
154
155
156
157
158
159
160
161
162
163
164
165
166
                                                                std::string featName,
                                                                double numerator,
                                                                std::map<std::string, double>* denominatorMap) const {
    std::map<std::string, double>::iterator it = denominatorMap->begin();
    for(; it!=denominatorMap->end(); it++) {
        std::string FullName = featName + it->first;
        float       value    = (float)it->second;
        if(value <= 0. || isnan(value) || isinf(value) ) continue; 
        targetMap->addFeature(FullName, numerator / it->second);
    }
}



167
StatusCode PanTau::Tool_FeatureExtractor::execute(PanTau::PanTauSeed2* inSeed) {
168
169
170
171
    
    ATH_MSG_DEBUG("Calculating features...");
    
    
172
173
174
    bool noAnyConstituents           = inSeed->isOfTechnicalQuality(PanTau::PanTauSeed2::t_NoConstituentsAtAll);
    bool noSelConstituents           = inSeed->isOfTechnicalQuality(PanTau::PanTauSeed2::t_NoSelectedConstituents);
    bool noValidInputTau             = inSeed->isOfTechnicalQuality(PanTau::PanTauSeed2::t_NoValidInputTau);
175
176
177
178
179
180
181
182
183
184
185
186
    bool isBadSeed                   = (noAnyConstituents || noSelConstituents || noValidInputTau);
    if(m_Config_UseEmptySeeds == true) isBadSeed = noValidInputTau;
    
    if(isBadSeed == true) {
        ATH_MSG_DEBUG("Seed is not valid for feature extraction (no constituents or no valid input tau) - just fill isPanTauCandidate feature");
        inSeed->getFeatures()->addFeature(inSeed->getNameInputAlgorithm() + "_" + m_varTypeName_Basic + "_isPanTauCandidate", 0);
        return StatusCode::SUCCESS;
    }
    inSeed->getFeatures()->addFeature(inSeed->getNameInputAlgorithm() + "_" + m_varTypeName_Basic + "_isPanTauCandidate", 1);
    
    
    ATH_MSG_DEBUG("Basic features");
187
    ATH_CHECK( calculateBasicFeatures(inSeed) );
188
189
190
    
    
    ATH_MSG_DEBUG("RawConstituent 4 vectors");
191
    ATH_CHECK( addConstituentMomenta(inSeed) );
192
193
194
195
196
197
198
    
    //first, calculate the Et variants for the seed
    fillVariantsSeedEt(inSeed->getConstituentsAsList_All());
    
    //loop through all types of Constituents in tau and calculate type features for them
    ATH_MSG_DEBUG("type specific features");
    //baseline
199
200
201
202
    ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent2::t_NoType) );  //=> all constituents
    ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent2::t_Charged) ); //=> charged ones in core
    ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent2::t_Neutral) ); //=> neutral ones in core
    ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent2::t_Pi0Neut) ); //=> pi0 tagged ones in core
203
    //for testing
204
205
    ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent2::t_NeutLowA) ); //=> same as neutral but with lower Et
    ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent2::t_NeutLowB) ); //=> same as neutral but with even lower et
206
207
208
209
    
    
    //fill the combined features
    ATH_MSG_DEBUG("combined features");
210
    ATH_CHECK( addCombinedFeatures(inSeed) );
211
212
213
    
    //fill the generic jet features
    ATH_MSG_DEBUG("generic jet features");
214
    ATH_CHECK( addGenericJetFeatures(inSeed) );
215
216
217
    
    //fill the impact paramter features
    ATH_MSG_DEBUG("impact parameter features");
218
    ATH_CHECK( addImpactParameterFeatures(inSeed) );
219
220
221
222
223
224
225
    
    ATH_MSG_DEBUG("Finished feature extraction");
    return StatusCode::SUCCESS;
}



226
StatusCode PanTau::Tool_FeatureExtractor::calculateBasicFeatures(PanTau::PanTauSeed2* inSeed) {
227
228
    
    ATH_MSG_DEBUG("calculating basic features");
229
    PanTau::TauFeature2* featureMap = inSeed->getFeatures();
230
231
232
233
234
235
236
237
238
239
240
241
    
    std::string featureAlg    = inSeed->getNameInputAlgorithm();
    std::string featurePrefix = m_varTypeName_Basic;
    
    
    
    double SumCharge = 0;
    double AbsCharge = 0;
    
    //! Loop over types to fill
    //!     - multiplicity of that type
    //!     - sum charge and abs charge
242
    for(int iType=0; iType<PanTau::TauConstituent2::t_nTypes; iType++) {
243
244
        
        bool foundIt = false;
245
        std::vector<TauConstituent2*> curList = inSeed->getConstituentsOfType(iType, foundIt);
246
247
248
        if(foundIt == false) continue;
        
        //store multiplicity of current type
249
        std::string     typeName        = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)iType);
250
251
252
253
254
        unsigned int    nConstituents   = curList.size();
        featureMap->addFeature(featureAlg + "_" + featurePrefix + "_N" + typeName + "Consts", nConstituents);
        
        
        //count charge, i.e. skip if not charged
255
        if(iType != (int)PanTau::TauConstituent2::t_Charged) continue;
256
        for(unsigned int iConst=0; iConst<nConstituents; iConst++) {
257
            PanTau::TauConstituent2* curConstituent = curList[iConst];
258
259
260
261
262
263
264
265
266
267
268
            SumCharge += curConstituent->getCharge();
            AbsCharge += fabs((double)curConstituent->getCharge());
        }
    }
    
    //store charge
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_SumCharge", SumCharge);
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_AbsCharge", AbsCharge);
    
    //! Fill multiplicity for any constituents
    //all constituents
269
    std::string                     typeNameAll         = PanTau::TauConstituent2::AllConstituentsName();
270
271
272
273
274
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_N" + typeNameAll + "Consts", inSeed->getConstituentsAsList_Core().size() + inSeed->getConstituentsAsList_Wide().size());
    
    //! Fill the proto vector (i.e. sum momentum of constituents)
    //proto 4-vector (just the sum of all constituents)
    // will have better four momentum after mode ID
275
276
277
278
279
    TLorentzVector tlv_ProtoMomentumCore = inSeed->getProtoMomentumCore();
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumCore_pt", tlv_ProtoMomentumCore.Perp());
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumCore_eta", tlv_ProtoMomentumCore.Eta());
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumCore_phi", tlv_ProtoMomentumCore.Phi());
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumCore_m", tlv_ProtoMomentumCore.M());
280
    
281
282
283
284
285
    TLorentzVector tlv_ProtoMomentumWide = inSeed->getProtoMomentumWide();
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumWide_pt", tlv_ProtoMomentumWide.Perp());
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumWide_eta", tlv_ProtoMomentumWide.Eta());
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumWide_phi", tlv_ProtoMomentumWide.Phi());
    featureMap->addFeature(featureAlg + "_" + featurePrefix + "_ProtoMomentumWide_m", tlv_ProtoMomentumWide.M());
286
287
288
289
290
291
292
293
    
    
    return StatusCode::SUCCESS;
}




294
StatusCode PanTau::Tool_FeatureExtractor::addConstituentMomenta(PanTau::PanTauSeed2* inSeed) {
295
    std::string inputAlgName  = inSeed->getNameInputAlgorithm();
296
    TauFeature2* tauFeatureMap = inSeed->getFeatures();
297
    std::string prefixVARType = m_varTypeName_HLV;
298
    for(int iType=0; iType<(int)PanTau::TauConstituent2::t_nTypes; iType++) {
299
        bool isOK;
300
        std::vector<PanTau::TauConstituent2*>    list_TypeConstituents   = inSeed->getConstituentsOfType(iType, isOK); // list of constituents of current type
301
        unsigned int                            n_Constituents_Type     = list_TypeConstituents.size();          // number of objects of current type
302
303
        TLorentzVector                          tlv_TypeConstituents    = inSeed->getSubsystemHLV(iType, isOK); // summed hlv of objects of current type
        std::string                             curTypeName             = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)iType);
304
        
305
        std::vector<PanTau::TauConstituent2*>    list_TypeConstituents_SortBDT = inSeed->getConstituentsOfType(iType, isOK);
306
307
308
        std::sort(list_TypeConstituents_SortBDT.begin(), list_TypeConstituents_SortBDT.end(), sortTauConstituentMVA);
        
        if(list_TypeConstituents.size() > 0) {
309
310
311
312
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumPt",  tlv_TypeConstituents.Perp());
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumEta", tlv_TypeConstituents.Eta());
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumPhi", tlv_TypeConstituents.Phi());
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumM",   tlv_TypeConstituents.M());
313
314
315
316
317
318
319
320
        }
        
        //store 4-vectors of current type (et sort);
        std::vector<double> curConsts_pt    = std::vector<double>(0);
        std::vector<double> curConsts_eta   = std::vector<double>(0);
        std::vector<double> curConsts_phi   = std::vector<double>(0);
        std::vector<double> curConsts_m     = std::vector<double>(0);
        for(unsigned int iConst=0; iConst<n_Constituents_Type; iConst++) {
321
322
323
324
325
            TLorentzVector tlv_curConst = list_TypeConstituents[iConst]->p4();
            curConsts_pt.push_back(tlv_curConst.Perp());
            curConsts_eta.push_back(tlv_curConst.Eta());
            curConsts_phi.push_back(tlv_curConst.Phi());
            curConsts_m.push_back(tlv_curConst.M());
326
327
328
329
330
331
        }
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtSort_Constituents_pt", curConsts_pt);
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtSort_Constituents_eta", curConsts_eta);
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtSort_Constituents_phi", curConsts_phi);
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtSort_Constituents_m", curConsts_m);
        
332

333
334
335
336
337
338
        //store 4-vectors of current type (bdt sort)
        std::vector<double> curConstsBDT_pt    = std::vector<double>(0);
        std::vector<double> curConstsBDT_eta   = std::vector<double>(0);
        std::vector<double> curConstsBDT_phi   = std::vector<double>(0);
        std::vector<double> curConstsBDT_m     = std::vector<double>(0);
        for(unsigned int iConst=0; iConst<n_Constituents_Type; iConst++) {
339
340
341
342
343
            TLorentzVector tlv_curConstBDT = list_TypeConstituents_SortBDT[iConst]->p4();
            curConstsBDT_pt.push_back(tlv_curConstBDT.Perp());
            curConstsBDT_eta.push_back(tlv_curConstBDT.Eta());
            curConstsBDT_phi.push_back(tlv_curConstBDT.Phi());
            curConstsBDT_m.push_back(tlv_curConstBDT.M());
344
345
346
347
348
349
350
351
352
353
354
355
356
        }
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTSort_Constituents_pt", curConstsBDT_pt);
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTSort_Constituents_eta", curConstsBDT_eta);
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTSort_Constituents_phi", curConstsBDT_phi);
        tauFeatureMap->addVecFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTSort_Constituents_m", curConstsBDT_m);
    
    } //end loop over constituent types
    
    return StatusCode::SUCCESS;
}



357
StatusCode PanTau::Tool_FeatureExtractor::calculateFeatures(PanTau::PanTauSeed2* inSeed,
358
359
360
361
362
                                                            int                 tauConstituentType) {
    //
    
    
    
363
364
365
    std::string                             curTypeName                     = PanTau::TauConstituent2::getTypeName( (PanTau::TauConstituent2::Type)tauConstituentType );
    std::string                             curTypeName_All                 = PanTau::TauConstituent2::AllConstituentsName();
    PanTau::TauFeature2*                     tauFeatureMap                   = inSeed->getFeatures();
366
    std::string                             inputAlgName                    = inSeed->getNameInputAlgorithm();
367
    TLorentzVector                          tlv_Reference                   = inSeed->getProtoMomentumCore();
368
    
369
    std::vector<PanTau::TauConstituent2*>    list_AllConstituents            = inSeed->getConstituentsAsList_Core();
370
371
    
    bool                                    foundIt                         = false;
372
373
374
    std::vector<PanTau::TauConstituent2*>    list_TypeConstituents;
    if(tauConstituentType != PanTau::TauConstituent2::t_NoType) list_TypeConstituents = inSeed->getConstituentsOfType(tauConstituentType, foundIt);
    if(tauConstituentType == PanTau::TauConstituent2::t_NoType) list_TypeConstituents = list_AllConstituents;
375
376
377
378
379
380
381
382
383
384
385
    if(foundIt == false) return StatusCode::SUCCESS;

    unsigned int                            n_Constituents_All              = list_AllConstituents.size();
    unsigned int                            n_Constituents_Type             = list_TypeConstituents.size();
    
    ATH_MSG_DEBUG("Calculating features for " << list_TypeConstituents.size() << " constituents of type " << tauConstituentType << "(" << curTypeName << ")");
    
    //sort the lists by Et
    std::sort(list_AllConstituents.begin(),     list_AllConstituents.end(),     sortTauConstituentEt);
    std::sort(list_TypeConstituents.begin(),    list_TypeConstituents.end(),    sortTauConstituentEt);
    
386
387
388
    TLorentzVector  tlv_1st_Et;
    TLorentzVector  tlv_2nd_Et;
    TLorentzVector  tlv_3rd_Et;
389
    
390
391
392
    if(list_TypeConstituents.size() > 0) tlv_1st_Et = list_TypeConstituents[0]->p4();
    if(list_TypeConstituents.size() > 1) tlv_2nd_Et = list_TypeConstituents[1]->p4();
    if(list_TypeConstituents.size() > 2) tlv_3rd_Et = list_TypeConstituents[2]->p4();
393
394
    
    
395
396
    TLorentzVector tlv_Last_Et;
    if(list_TypeConstituents.size() > 0) tlv_Last_Et = list_TypeConstituents.back()->p4();
397
398
    
    //make an additional list of constituents, but now ordered by BDT value
399
    std::vector<PanTau::TauConstituent2*>    list_TypeConstituents_SortBDT = list_TypeConstituents;
400
401
    std::sort(list_TypeConstituents_SortBDT.begin(), list_TypeConstituents_SortBDT.end(), sortTauConstituentMVA);
    
402
403
404
    TLorentzVector  tlv_1st_BDT;
    TLorentzVector  tlv_2nd_BDT;
    TLorentzVector  tlv_3rd_BDT;
405
    
406
407
408
    if(list_TypeConstituents_SortBDT.size() > 0) tlv_1st_BDT = list_TypeConstituents_SortBDT[0]->p4();
    if(list_TypeConstituents_SortBDT.size() > 1) tlv_2nd_BDT = list_TypeConstituents_SortBDT[1]->p4();
    if(list_TypeConstituents_SortBDT.size() > 2) tlv_3rd_BDT = list_TypeConstituents_SortBDT[2]->p4();
409
410
411
412
413
    
    
    
    
    
414
    //! //////////////////////////////////////////                  
415
416
417
418
    //! Prepare variables for information from eflow Objects
    //! //////////////////////////////////////////
    
    // ===> hlv for the leading EFOs and the summed HLV
419
    TLorentzVector              tlv_TypeConstituents;
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
    // ===> Sum of DeltaR to jet axis
    double                      sum_DRToReference             = 0;
    double                      sum_DR2ToReference            = 0;
    double                      sum_DRToLeading             = 0;
    double                      sum_DR2ToLeading            = 0;
    // ===> Sum of Et, Et^2, E and E^2
    double                      sum_Et                      = 0;
    double                      sum_Et2                     = 0;
    double                      sum_E                       = 0;
    double                      sum_E2                      = 0;
    // ===> Sum of Et (and E) times DeltaR, DeltaR', Angle 
    double                      sum_EtxDR                   = 0;
    double                      sum_EtxDR2                  = 0;
    double                      sum_EtxDRprime              = 0;
    double                      sum_EtxAngle                = 0;
    double                      sum_ExDR                    = 0;
    double                      sum_ExDR2                   = 0;
    double                      sum_ExDRprime               = 0;
    double                      sum_ExAngle                 = 0;
    // ===> Isolation rings
    double                      sum_EtInRing00To01          = 0;
    double                      sum_EtInRing01To02          = 0;
    double                      sum_EtInRing02To03          = 0;
    double                      sum_EtInRing03To04          = 0;
    double                      sum_EtInRing04To05          = 0;
    // ===> Multiplicities
    unsigned int                num_EFOs                    = 0;
    unsigned int                num_ConstsIn00To01             = 0;
    unsigned int                num_ConstsIn01To02             = 0;
    unsigned int                num_ConstsIn02To03             = 0;
    unsigned int                num_ConstsIn03To04             = 0;
    unsigned int                num_ConstsIn04To05             = 0;
    // ===> Maximal values
    double                      max_DeltaR                  = 0;
    
    
    //! //////////////////////////////////////////
    //! Loop over selected constituents and collect information
    //! //////////////////////////////////////////
    for(unsigned int iTypeConst=0; iTypeConst<list_TypeConstituents.size(); iTypeConst++) {
        
        //get hep lorentz vector
462
        TLorentzVector tlv_curConst = list_TypeConstituents.at(iTypeConst)->p4();
463
464
        
        //final check (nan & inf)
465
        if (isnan(tlv_curConst.Pt()) || isinf(tlv_curConst.Pt())) continue;
466
467
468
469
        
        //ready to calc stuff
        //basically update each of the prepared sum_* and num_* variables above,
        // the sum HLV and the pointers to 1st, 2nd, 3rd leading constituents of current type
470
        tlv_TypeConstituents += tlv_curConst;
471
472
        
        //helpers to reduce function calls
473
        double hlp_Et               = tlv_curConst.Et();
474
        double hlp_Et2              = hlp_Et * hlp_Et;
475
        double hlp_E                = tlv_curConst.E();
476
        double hlp_E2               = hlp_E * hlp_E;
477
        double hlp_DeltaR           = tlv_Reference.DeltaR(tlv_curConst);
478
        double hlp_DeltaR2          = hlp_DeltaR * hlp_DeltaR;
479
        double hlp_DeltaRLeading    = (tlv_1st_Et.Pt() == 0 ? 0 : tlv_1st_Et.DeltaR(tlv_curConst));
480
        double hlp_DeltaR2Leading   = hlp_DeltaRLeading * hlp_DeltaRLeading;
481
482
        double hlp_DeltaRprime      = m_HelperFunctions.deltaRprime(tlv_Reference.Vect(), tlv_curConst.Vect());
        double hlp_Angle            = tlv_Reference.Angle(tlv_curConst.Vect());
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
        
        // update sum of DeltaR to jet axis
        sum_DRToReference           += hlp_DeltaR;
        sum_DR2ToReference          += hlp_DeltaR2;
        sum_DRToLeading             += hlp_DeltaRLeading;
        sum_DR2ToLeading            += hlp_DeltaR2Leading;
        // update Sum of Et, Et^2, E and E^2
        sum_Et                      += hlp_Et;
        sum_Et2                     += hlp_Et2;
        sum_E                       += hlp_E;
        sum_E2                      += hlp_E2;
        // update Sum of Et (and E) times DeltaR, DeltaR', Angle 
        sum_EtxDR                   += hlp_Et * hlp_DeltaR;
        sum_EtxDR2                  += hlp_Et * hlp_DeltaR2;
        sum_EtxDRprime              += hlp_Et * hlp_DeltaRprime;
        sum_EtxAngle                += hlp_Et * hlp_Angle;
        sum_ExDR                    += hlp_E  * hlp_DeltaR;
        sum_ExDR2                   += hlp_E  * hlp_DeltaR2;
        sum_ExDRprime               += hlp_E  * hlp_DeltaRprime;
        sum_ExAngle                 += hlp_E  * hlp_Angle;
        // update Isolation rings
        if(hlp_DeltaR >= 0.0 && hlp_DeltaR < 0.1) sum_EtInRing00To01 += hlp_Et;
        if(hlp_DeltaR >= 0.1 && hlp_DeltaR < 0.2) sum_EtInRing01To02 += hlp_Et;
        if(hlp_DeltaR >= 0.2 && hlp_DeltaR < 0.3) sum_EtInRing02To03 += hlp_Et;
        if(hlp_DeltaR >= 0.3 && hlp_DeltaR < 0.4) sum_EtInRing03To04 += hlp_Et;
        if(hlp_DeltaR >= 0.4 && hlp_DeltaR < 0.5) sum_EtInRing04To05 += hlp_Et;
        // update Multiplicities
        num_EFOs++;
        if(hlp_DeltaR >= 0.0 && hlp_DeltaR < 0.1) num_ConstsIn00To01++;
        if(hlp_DeltaR >= 0.1 && hlp_DeltaR < 0.2) num_ConstsIn01To02++;
        if(hlp_DeltaR >= 0.2 && hlp_DeltaR < 0.3) num_ConstsIn02To03++;
        if(hlp_DeltaR >= 0.3 && hlp_DeltaR < 0.4) num_ConstsIn03To04++;
        if(hlp_DeltaR >= 0.4 && hlp_DeltaR < 0.5) num_ConstsIn04To05++;
        // update Max values
        if(hlp_DeltaR > max_DeltaR) max_DeltaR = hlp_DeltaR;
    }//end loop over selected EFOs
    
    
    
    //! //////////////////////////////////////////
    //! Calculate & fill features
    //! //////////////////////////////////////////
    
    // naming string
    std::string prefixVARType = "";
    
    //! Multiplicities ///////////////////////////////////////////
    //if there are no EFOs of this type in the seed, all other variables cannot be calculated to have any reasonable value
    if( num_EFOs == 0 ) {
        return StatusCode::SUCCESS;
    }
    
    prefixVARType = m_varTypeName_Num;
    tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ConstsIn00To01", num_ConstsIn00To01);
    tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ConstsIn01To02", num_ConstsIn01To02);
    tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ConstsIn02To03", num_ConstsIn02To03);
    tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ConstsIn03To04", num_ConstsIn03To04);
    
    
    
    //! Substructure particle ID features ///////////////////////////////////////////
    prefixVARType = m_varTypeName_PID;
    
    // Sorted by highest BDT score
    double value_sumBDT_BDTSort = 0;
    for(unsigned int iTypeConst=0; iTypeConst<n_Constituents_Type; iTypeConst++) {
        
        double value_BDT = list_TypeConstituents_SortBDT[iTypeConst]->getBDTValue();
        if( isnan(value_BDT) || isinf(value_BDT) ) continue;
        
        //correct BDT value based on BDT cut
554
        if(tauConstituentType != PanTau::TauConstituent2::t_Charged) {
555
            double mvaCorrection = 0.0;
556
557
            double  etaCurConst = list_TypeConstituents_SortBDT[iTypeConst]->p4().Eta();
            int     etaBinIndex = m_HelperFunctions.getBinIndex(m_Config_CellBased_BinEdges_Eta, fabs(etaCurConst));
558
            bool    isOK;
559
            int     numTrack    = inSeed->getConstituentsOfType(PanTau::TauConstituent2::t_Charged, isOK).size();
560
561
562
563
564
565
            if(numTrack == 1) { mvaCorrection = m_Config_CellBased_EtaBinned_Pi0MVACut_1prong.at(etaBinIndex); }
            else              { mvaCorrection = m_Config_CellBased_EtaBinned_Pi0MVACut_3prong.at(etaBinIndex); }
            value_BDT = value_BDT - mvaCorrection;
        }
        
        value_sumBDT_BDTSort += value_BDT;
566
        std::string iConst = m_HelperFunctions.convertNumberToString((double)(iTypeConst+1));
567
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTValues_BDTSort_" + iConst, value_BDT);
568
	//ATH_MSG_DEBUG("\t\tAdded variable " << inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_nPhotons_BDTSort_" + iConstStr << " with value " << totalPhotonsInNeutral);
569
570
571
572
573
574
575
576
577
578
579
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTValuesSum_BDTSort_" + iConst, value_sumBDT_BDTSort);
    }
    
    //sorted by highest Et
    double value_sumBDT_EtSort = 0;
    for(unsigned int iTypeConst=0; iTypeConst<n_Constituents_Type; iTypeConst++) {
        
        double value_BDT = list_TypeConstituents[iTypeConst]->getBDTValue();
        if( isnan(value_BDT) || isinf(value_BDT) ) continue;
        
        value_sumBDT_EtSort += value_BDT;
580
        std::string iConst = m_HelperFunctions.convertNumberToString((double)(iTypeConst+1));
581
582
583
584
585
586
587
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTValues_EtSort_" + iConst, value_BDT);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BDTValuesSum_EtSort_" + iConst, value_sumBDT_EtSort);
    }
    
    
    //! Shot information ///////////////////////////////////////////
    prefixVARType = PanTau::Tool_FeatureExtractor::varTypeName_Shots();
588
589
    //only execute if the constituent type is neutral
    if(PanTau::TauConstituent2::isNeutralType(tauConstituentType) == true) {
590
591
        ATH_MSG_DEBUG("---> Dumping shot information from " << list_TypeConstituents_SortBDT.size() << " constituents of type " << curTypeName << " in tau");
        
592
        TLorentzVector          totalTLV_SumShots       = TLorentzVector(0., 0., 0., 0.);
593
594
595
596
597
598
599
        unsigned int            totalPhotonsInSeed      = 0;
        unsigned int            totalShotsInSeed        = 0;
        double                  maxDeltaRSumShotToConst = -999;
        double                  minDeltaRSumShotToConst = 999;
        double                  maxDeltaRSumShotToTau   = -999;
        double                  minDeltaRSumShotToTau   = 999;
        
600
        std::vector<TLorentzVector> allShotTLVs = std::vector<TLorentzVector>(0);
601
602
603
604
        
        for(unsigned int iConst=0; iConst<list_TypeConstituents_SortBDT.size(); iConst++) {
            ATH_MSG_DEBUG("\tConstituent " << iConst << " / " << list_TypeConstituents_SortBDT.size());
            
605
606
607
            PanTau::TauConstituent2*                 curConst            = list_TypeConstituents_SortBDT.at(iConst);
            TLorentzVector                          tlv_CurConst        = curConst->p4();
            std::vector<PanTau::TauConstituent2*>    shotConstituents    = curConst->getShots();
608
            unsigned int                            nShots              = shotConstituents.size();
609
            ATH_MSG_DEBUG("\t\tConstituent has Pt/Eta/Phi/M: " << tlv_CurConst.Pt() << " / " << tlv_CurConst.Eta() << " / " << tlv_CurConst.Phi() << " / " << tlv_CurConst.M());
610
611
612
            ATH_MSG_DEBUG("\t\tShots in this constituent: " << nShots);
            
            unsigned int totalPhotonsInNeutral = 0;
613
            TLorentzVector tlv_SumShots = TLorentzVector(0., 0., 0., 0.);
614
615
            
            for(unsigned int iShot=0; iShot<nShots; iShot++) {
616
                PanTau::TauConstituent2* curShot     = shotConstituents.at(iShot);
617
618
619
                unsigned int            curNPhotons = curShot->getNPhotonsInShot();
                ATH_MSG_DEBUG("\t\t\tPhotons in this shot: " << curNPhotons);
                totalPhotonsInNeutral += curShot->getNPhotonsInShot();
620
621
622
                ATH_MSG_DEBUG("\t\t\tPt/Eta/Phi/M of this shot: " << curShot->p4().Pt() << " / " << curShot->p4().Eta() << " / " << curShot->p4().Phi() << " / " << curShot->p4().M() );
                tlv_SumShots += curShot->p4();
                allShotTLVs.push_back(curShot->p4());
623
624
            }//end loop over shots
            totalShotsInSeed    += nShots;
625
            totalTLV_SumShots   += tlv_SumShots;
626
627
628
            totalPhotonsInSeed  += totalPhotonsInNeutral;
            
            ATH_MSG_DEBUG("\t\tTotal Photons: " << totalPhotonsInNeutral);
629
            ATH_MSG_DEBUG("\t\tPt/Eta/Phi/M of combined shots: " << tlv_SumShots.Pt() << " / " << tlv_SumShots.Eta() << " / " << tlv_SumShots.Phi() << " / " << tlv_SumShots.M() );
630
            
631
632
            std::string iConstStr = m_HelperFunctions.convertNumberToString((double)(iConst+1));
                       
633
634
635
636
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_nPhotons_BDTSort_" + iConstStr, totalPhotonsInNeutral);
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_nShots_BDTSort_" + iConstStr, nShots);
            
            //the et/eta/phi/m of the hlv of all shots combined for this neutral-type constituent
637
638
639
640
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumShots_Et_BDTSort_" + iConstStr, tlv_SumShots.Et());
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumShots_Eta_BDTSort_" + iConstStr, tlv_SumShots.Eta());
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumShots_Phi_BDTSort_" + iConstStr, tlv_SumShots.Phi());
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SumShots_M_BDTSort_" + iConstStr, tlv_SumShots.M());
641
642
            
            //energy ratio, deltaR of sumShots and constituent
643
            double deltaRSumShotToConst = tlv_CurConst.DeltaR(tlv_SumShots);
644
645
646
            if(deltaRSumShotToConst > maxDeltaRSumShotToConst) maxDeltaRSumShotToConst = deltaRSumShotToConst;
            if(deltaRSumShotToConst < minDeltaRSumShotToConst) minDeltaRSumShotToConst = deltaRSumShotToConst;
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ConstDeltaRToSumShots_BDTSort_" + iConstStr, deltaRSumShotToConst);
647
            if(tlv_CurConst.Et() > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtSumShotsOverConstEt_BDTSort_" + iConstStr, tlv_SumShots.Et() / tlv_CurConst.Et());
648
649
            
            //energy ratio, deltaR of shots and tauSeed
650
            double deltaRSumShotToTau = tlv_Reference.DeltaR(tlv_SumShots);
651
652
653
            if(deltaRSumShotToTau > maxDeltaRSumShotToTau) maxDeltaRSumShotToTau = deltaRSumShotToTau;
            if(deltaRSumShotToTau < minDeltaRSumShotToTau) minDeltaRSumShotToTau = deltaRSumShotToTau;
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_TauDeltaRToSumShots_BDTSort_" + iConstStr, deltaRSumShotToTau);
654
            if(tlv_Reference.Et() > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtSumShotsOverTauEt_BDTSort_" + iConstStr, tlv_SumShots.Et() / tlv_Reference.Et());
655
656
657
658
659
660
661
662
663
664
665
            
        }//end loop over constituents in tau
        
        ATH_MSG_DEBUG("\tLoop over constituents in tau finished!");
        ATH_MSG_DEBUG("\tTotal photons from shots: " << totalPhotonsInSeed);
        
        //delta R values
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MaxDeltaRSumShotToConst", maxDeltaRSumShotToConst);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MinDeltaRSumShotToConst", minDeltaRSumShotToConst);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MaxDeltaRSumShotToTau", maxDeltaRSumShotToTau);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MinDeltaRSumShotToTau", minDeltaRSumShotToTau);
666
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_DeltaRAllShotsToTau", tlv_Reference.DeltaR(totalTLV_SumShots));
667
668
        
        //et ratio
669
        if(tlv_Reference.Et() > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtAllShotsOverEtTau", totalTLV_SumShots.Et() / tlv_Reference.Et());
670
671
672
673
674
675
676
677
678
679
680
        
        //number of shots in seed
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_NShotsInSeed", totalShotsInSeed);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_NPhotonsInSeed", totalPhotonsInSeed);

        //build di-Shot mass
        ATH_MSG_DEBUG("\tBuild di-shot masses and check difference to pion mass");
        double maxDiShotMass    = -200;
        double minDiShotMass    = 99999;
        double bestDiShotMass   = -200;
        double bestPi0Diff      = 99999;
681
682
        for(unsigned int iShot=0; iShot<allShotTLVs.size(); iShot++) {
            TLorentzVector cur_iShot = allShotTLVs.at(iShot);
683
            
684
685
            for(unsigned int jShot=iShot+1; jShot<allShotTLVs.size(); jShot++) {
                TLorentzVector cur_jShot = allShotTLVs.at(jShot);
686
687
                
                ATH_MSG_DEBUG("\t\tBuilding di-shot mass of shots " << iShot << " & " << jShot);
688
689
                TLorentzVector          tlv_DiShot    = cur_iShot + cur_jShot;
                double                  curDiShotMass = tlv_DiShot.M();
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
                double                  curpi0Diff    = fabs(curDiShotMass - 134.98);
                ATH_MSG_DEBUG("\t\tit is: " << curDiShotMass);
                if(curpi0Diff < bestPi0Diff) bestDiShotMass = curDiShotMass;
                if(curDiShotMass > maxDiShotMass) maxDiShotMass = curDiShotMass;
                if(curDiShotMass < minDiShotMass) minDiShotMass = curDiShotMass;
            }
        }
        ATH_MSG_DEBUG("\tBest di-shot mass: " << bestDiShotMass);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_BestDiShotMass", bestDiShotMass);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MaxDiShotMass", maxDiShotMass);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MinDiShotMass", minDiShotMass);
        
    }//end if check for shot info dumping
    
    
    //! Ratios ///////////////////////////////////////////
    prefixVARType = m_varTypeName_Ratio;
    
    if(curTypeName != curTypeName_All) addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtOver", sum_Et, &m_Variants_SeedEt);
    
710
711
    if(tlv_1st_Et.Pt() != 0) addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stEtOver",   tlv_1st_Et.Et(), &m_Variants_SeedEt);
    if(tlv_1st_BDT.Pt() != 0) addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stBDTEtOver",   tlv_1st_BDT.Et(), &m_Variants_SeedEt);
712
    
713
    if(tlv_Last_Et.Pt() != 0) addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_SoftestEtOver",   tlv_Last_Et.Et(), &m_Variants_SeedEt);
714
    
715
716
    if(tlv_1st_Et.Pt() != 0 && sum_Et > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stEtOverTypeEt",    tlv_1st_Et.Et() / sum_Et);
    if(tlv_1st_BDT.Pt() != 0 && sum_Et > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stBDTEtOverTypeEt",    tlv_1st_BDT.Et() / sum_Et);
717
718
719
720
    
    
    
    if(n_Constituents_All != 0 && curTypeName != curTypeName_All)   tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EFOsOverTotalEFOs", (double)(((double)num_EFOs) / ((double)n_Constituents_All)));
721
722
723
    if(tlv_1st_Et.Pt() != 0 && tlv_2nd_Et.Pt() != 0) {
        if(tlv_1st_Et.Et() > 0. && tlv_2nd_Et.Et() > 0. ) {
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Log1stEtOver2ndEt", TMath::Log10(tlv_1st_Et.Et() / tlv_2nd_Et.Et()));
724
725
        }
    }
726
727
728
    if(tlv_1st_Et.Pt() != 0 && tlv_3rd_Et.Pt() != 0) {
        if(tlv_1st_Et.Et() > 0. && tlv_3rd_Et.Et() > 0.) {
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Log1stEtOver3rdEt", TMath::Log10(tlv_1st_Et.Et() / tlv_3rd_Et.Et()));
729
730
        }
    }
731
732
733
    if(tlv_2nd_Et.Pt() != 0 && tlv_3rd_Et.Pt() != 0) {
        if(tlv_2nd_Et.Et() > 0. && tlv_3rd_Et.Et() > 0.) {
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Log2ndEtOver3rdEt", TMath::Log10(tlv_2nd_Et.Et() / tlv_3rd_Et.Et()));
734
735
736
737
        }
    }
    
    //and for the BDT score ordered EFOs
738
739
740
    if(tlv_1st_BDT.Pt() != 0 && tlv_2nd_BDT.Pt() != 0) {
        if(tlv_1st_BDT.Et() > 0. && tlv_2nd_BDT.Et() > 0. ) {
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Log1stEtOver2ndEt_BDTSort", TMath::Log10(tlv_1st_BDT.Et() / tlv_2nd_BDT.Et()));
741
742
        }
    }
743
744
745
    if(tlv_1st_BDT.Pt() != 0 && tlv_3rd_BDT.Pt() != 0) {
        if(tlv_1st_BDT.Et() > 0. && tlv_3rd_BDT.Et() > 0.) {
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Log1stEtOver3rdEt_BDTSort", TMath::Log10(tlv_1st_BDT.Et() / tlv_3rd_BDT.Et()));
746
747
        }
    }
748
749
750
    if(tlv_2nd_BDT.Pt() != 0 && tlv_3rd_BDT.Pt() != 0) {
        if(tlv_2nd_BDT.Et() > 0. && tlv_3rd_BDT.Et() > 0.) {
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Log2ndEtOver3rdEt_BDTSort", TMath::Log10(tlv_2nd_BDT.Et() / tlv_3rd_BDT.Et()));
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
        }
    }
    
    
    //! EtRings  ///////////////////////////////////////////
    if(curTypeName == curTypeName_All) {
        prefixVARType = m_varTypeName_EtInRing;
        
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_00To01", sum_EtInRing00To01);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_01To02", sum_EtInRing01To02);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_02To03", sum_EtInRing02To03);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_03To04", sum_EtInRing03To04);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_04To05", sum_EtInRing04To05);
    }
    
    
    //! Isolations ///////////////////////////////////////////
    if(curTypeName == curTypeName_All) {
        prefixVARType = m_varTypeName_Isolation;
        
        double iso_EtIn01 = sum_EtInRing00To01;
        double iso_EtIn02 = iso_EtIn01 + sum_EtInRing01To02;
        double iso_EtIn03 = iso_EtIn02 + sum_EtInRing02To03;
        double iso_EtIn04 = iso_EtIn03 + sum_EtInRing03To04;
        
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn01Over", iso_EtIn01, &m_Variants_SeedEt);
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn02Over", iso_EtIn02, &m_Variants_SeedEt);
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn03Over", iso_EtIn03, &m_Variants_SeedEt);
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn04Over", iso_EtIn04, &m_Variants_SeedEt);
        
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn00To02Over", (sum_EtInRing00To01 + sum_EtInRing01To02), &m_Variants_SeedEt);
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn02To04Over", (sum_EtInRing02To03 + sum_EtInRing03To04), &m_Variants_SeedEt);
        
        if(iso_EtIn02>0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn01OverEtIn02", iso_EtIn01 / iso_EtIn02);
        if(iso_EtIn04>0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtIn01OverEtIn04", iso_EtIn01 / iso_EtIn04);
    }
    
    
    //! Means ///////////////////////////////////////////
    prefixVARType = m_varTypeName_Mean;
    
    if(num_EFOs > 0) {
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Et_Wrt",           (sum_Et / num_EFOs), &m_Variants_SeedEt);
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_DRToJetAxis_Wrt",  (sum_DRToReference / num_EFOs), &m_Variants_SeedEt);
        addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_DRToLeading_Wrt",  (sum_DRToLeading / num_EFOs), &m_Variants_SeedEt);
    }
    
    
    //! Standard deviations ///////////////////////////////////////////
    prefixVARType =  m_varTypeName_StdDev;

802
803
804
805
    double stddev_E             = m_HelperFunctions.stddev(sum_E2, sum_E, num_EFOs);
    double stddev_Et            = m_HelperFunctions.stddev(sum_Et2, sum_Et, num_EFOs);
    double stddev_DRToJetAxis   = m_HelperFunctions.stddev(sum_DR2ToReference, sum_DRToReference, num_EFOs);
    double stddev_DRToLeading   = m_HelperFunctions.stddev(sum_DRToLeading, sum_DR2ToLeading, num_EFOs);
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
    
    if(stddev_E > 0.)           tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_E",              stddev_E);
    if(stddev_Et > 0.)          tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Et",             stddev_Et);
    if(stddev_Et > 0.)          addFeatureWrtSeedEnergy(tauFeatureMap, inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Et_Wrt", stddev_Et, &m_Variants_SeedEt);
    if(stddev_DRToJetAxis > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_DRToJetAxis",    stddev_DRToJetAxis);
    if(stddev_DRToLeading > 0.) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_DRToLeading",    stddev_DRToLeading);
        
    
    //! Angles ///////////////////////////////////////////
    prefixVARType = m_varTypeName_Angle;
    
    double angle_12 = 0;
    double angle_13 = 0;
    double angle_23 = 0;
    
821
822
823
824
825
826
827
    if(curTypeName != curTypeName_All) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ToJetAxis",   tlv_Reference.Angle(tlv_TypeConstituents.Vect()));
    if(tlv_1st_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stToJetAxis",   tlv_Reference.Angle(tlv_1st_Et.Vect()));
    if(tlv_2nd_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_2ndToJetAxis",   tlv_Reference.Angle(tlv_2nd_Et.Vect()));
    if(tlv_3rd_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_3rdToJetAxis",   tlv_Reference.Angle(tlv_3rd_Et.Vect()));
    if(tlv_1st_Et.Pt() != 0) {
        if(tlv_2nd_Et.Pt() != 0) {
            angle_12 = tlv_1st_Et.Angle(tlv_2nd_Et.Vect());
828
829
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stTo2nd", angle_12);
        }
830
831
        if(tlv_3rd_Et.Pt() != 0) {
            angle_13 = tlv_1st_Et.Angle(tlv_3rd_Et.Vect());
832
833
834
            tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stTo3rd", angle_13);
        }
    }
835
836
    if(tlv_2nd_Et.Pt() != 0 && tlv_3rd_Et.Pt() != 0) {
        angle_23 = tlv_2nd_Et.Angle(tlv_3rd_Et.Vect());
837
838
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_2ndTo3rd", angle_23);
    }
839
840
    if(num_EFOs > 2 && tlv_1st_Et.Pt() != 0 && tlv_2nd_Et.Pt() != 0 && tlv_3rd_Et.Pt() != 0) {
        double angle_Planes = ( tlv_1st_Et.Vect().Cross(tlv_2nd_Et.Vect()) ).Angle( tlv_1st_Et.Vect().Cross(tlv_3rd_Et.Vect()) );
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
        double angle_max    = 0;
        if(angle_12 > angle_13) {
            if(angle_12 > angle_23) angle_max = angle_12;
            else angle_max =angle_23;
        } else {
            if(angle_13 > angle_23) angle_max =angle_13;
            else angle_max =angle_23;
        }
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MaxToJetAxis", angle_max);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MeanValue123", (angle_12 + angle_13 + angle_23)/3.);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_Btw1213Planes", angle_Planes);
        
    }
    
    
    //! DeltaR ///////////////////////////////////////////
    prefixVARType = m_varTypeName_DeltaR;
    
859
    if(curTypeName != curTypeName_All) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_ToJetAxis",      tlv_Reference.DeltaR(tlv_TypeConstituents));
860
    tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_MaxToJetAxis_EtSort",   max_DeltaR);
861
862
863
864
865
866
    if(tlv_1st_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stToJetAxis_EtSort",   tlv_Reference.DeltaR(tlv_1st_Et));
    if(tlv_2nd_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_2ndToJetAxis_EtSort",   tlv_Reference.DeltaR(tlv_2nd_Et));
    if(tlv_3rd_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_3rdToJetAxis_EtSort",   tlv_Reference.DeltaR(tlv_3rd_Et));
    if(tlv_1st_Et.Pt() != 0) {
        if(tlv_2nd_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stTo2nd_EtSort",   tlv_1st_Et.DeltaR(tlv_2nd_Et));
        if(tlv_3rd_Et.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stTo3rd_EtSort",   tlv_1st_Et.DeltaR(tlv_3rd_Et));
867
    }
868
869
    if(tlv_2nd_Et.Pt() != 0 && tlv_3rd_Et.Pt() != 0) {
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_2ndTo3rd_EtSort",   tlv_2nd_Et.DeltaR(tlv_3rd_Et));
870
871
    }
    
872
873
874
875
876
877
    if(tlv_1st_BDT.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stToJetAxis_BDTSort",   tlv_Reference.DeltaR(tlv_1st_BDT));
    if(tlv_2nd_BDT.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_2ndToJetAxis_BDTSort",   tlv_Reference.DeltaR(tlv_2nd_BDT));
    if(tlv_3rd_BDT.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_3rdToJetAxis_BDTSort",   tlv_Reference.DeltaR(tlv_3rd_BDT));
    if(tlv_1st_BDT.Pt() != 0) {
        if(tlv_2nd_BDT.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stTo2nd_BDTSort",   tlv_1st_BDT.DeltaR(tlv_2nd_BDT));
        if(tlv_3rd_BDT.Pt() != 0) tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_1stTo3rd_BDTSort",   tlv_1st_BDT.DeltaR(tlv_3rd_BDT));
878
    }
879
880
    if(tlv_2nd_BDT.Pt() != 0 && tlv_3rd_BDT.Pt() != 0) {
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_2ndTo3rd_BDTSort",   tlv_2nd_BDT.DeltaR(tlv_3rd_BDT));
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
    }
    
    
    //! JetMoment ///////////////////////////////////////////
    //  Moment wrt X = Sum( Et * X ) / Sum(Et)
    //  ==> Named as Eflow_EFOType_JetMoment_EtX"   when using Et as weight
    //  ==> Named as Eflow_EFOType_JetMoment_EX"    when using E  as weight
    prefixVARType = m_varTypeName_JetMoment;
    
    if(sum_Et > 0.) {
        //Using transverse energy
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtDR",            sum_EtxDR / sum_Et);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtDRprime",       sum_EtxDRprime / sum_Et);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtDR2",           sum_EtxDR2 / sum_Et);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtAngle",         sum_EtxAngle / sum_Et);
        tauFeatureMap->addFeature(inputAlgName + "_" + curTypeName + "_" + prefixVARType + "_EtDRxTotalEt",    (sum_EtxDR / sum_Et) * m_Variants_SeedEt["EtAllConsts"]);
    }
    
    
    return StatusCode::SUCCESS;
}



905
StatusCode PanTau::Tool_FeatureExtractor::addCombinedFeatures(PanTau::PanTauSeed2* inSeed) {
906
907
908
909
    
    //! //////////////////////////////////////////
    //! Prepare some short notations for variables
    //! //////////////////////////////////////////
910
    PanTau::TauFeature2* tauFeatures     = inSeed->getFeatures();
911
912
913
914
    std::string         inputAlgName    = inSeed->getNameInputAlgorithm();
    
    
    //et: EFO Type
915
916
917
918
    int et_Charged      = PanTau::TauConstituent2::t_Charged;
    int et_Pi0Neut      = PanTau::TauConstituent2::t_Pi0Neut;
    int et_Neutral      = PanTau::TauConstituent2::t_Neutral;
    int et_All          = PanTau::TauConstituent2::t_NoType;
919
920
    
    bool foundIt;
921
    std::vector<PanTau::TauConstituent2*>    list_NeutralConstituents = inSeed->getConstituentsOfType(et_Neutral, foundIt);
922
923
924
925
926
927
    
    //! //////////////////////////////////////////
    //! Prepare the list of names for EFO Types & 
    //! the 4 momenta of the different sub systems 
    //! (ie. charged, neutral subsystem, etc...)
    //! //////////////////////////////////////////
928
929
930
931
932
    std::string                 name_EFOType[PanTau::TauConstituent2::t_nTypes];
    double                      num_EFOs[PanTau::TauConstituent2::t_nTypes];
    TLorentzVector     tlv_System[PanTau::TauConstituent2::t_nTypes];
    TLorentzVector     tlv_1stEFO[PanTau::TauConstituent2::t_nTypes];
    TLorentzVector     tlv_2ndEFO[PanTau::TauConstituent2::t_nTypes];
933
934
935
936
    
    //! //////////////////////////////////////////
    //! get input objects to calc combined features
    //! //////////////////////////////////////////
937
938
939
    bool tlv_Sys_OK[PanTau::TauConstituent2::t_nTypes];
    bool tlv_1st_OK[PanTau::TauConstituent2::t_nTypes];
    bool tlv_2nd_OK[PanTau::TauConstituent2::t_nTypes];
940
941
942
    
    
    //initialize arrays with default values
943
    for(unsigned int iType=0; iType<(unsigned int)PanTau::TauConstituent2::t_nTypes; iType++) {
944
945
        name_EFOType[iType] = "";
        num_EFOs[iType]     = 0.;
946
947
948
949
950
951
        tlv_System[iType]   = TLorentzVector();
        tlv_1stEFO[iType]   = TLorentzVector();
        tlv_2ndEFO[iType]   = TLorentzVector();
        tlv_Sys_OK[iType]   = false;
        tlv_1st_OK[iType]   = false;
        tlv_2nd_OK[iType]   = false;
952
953
    }
    
954
955
    for(int iType=0; iType<(int)PanTau::TauConstituent2::t_nTypes; iType++) {
        name_EFOType[iType] = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)iType);
956
        
957
958
        tlv_System[iType] = inSeed->getSubsystemHLV(iType, tlv_Sys_OK[iType]);
        if(tlv_Sys_OK[iType] == false) continue;
959
        
960
961
962
        std::vector<TauConstituent2*> typeConstituents = inSeed->getConstituentsOfType(iType, tlv_Sys_OK[iType]);
        if(typeConstituents.size() == 0) tlv_Sys_OK[iType] = false;
        if(tlv_Sys_OK[iType] == false) continue;
963
964
965
966
        
        num_EFOs[iType] = typeConstituents.size();
        
        if(typeConstituents.size() > 0) {
967
968
            tlv_1stEFO[iType] = typeConstituents.at(0)->p4();
            tlv_1st_OK[iType] = true;
969
        } else {
970
            tlv_1st_OK[iType] = false;
971
972
973
        }
        
        if(typeConstituents.size() > 1) {
974
975
            tlv_2ndEFO[iType] = typeConstituents.at(1)->p4();
            tlv_2nd_OK[iType] = true;
976
        } else {
977
            tlv_2nd_OK[iType] = false;
978
979
980
981
982
983
984
985
986
987
988
989
990
        }
        
    }
    
    
    //! //////////////////////////////////////////
    //! From the extracted input, calc combined features
    //! //////////////////////////////////////////
    std::string prefixVARType = m_varTypeName_Combined;
    
    
    //! Combined-Single Features ////////////////////////////////////////
    // Ratios of numbers (heavily spiked, just keep them for validation)
991
    if(tlv_Sys_OK[et_Charged] == true && tlv_Sys_OK[et_Neutral] && num_EFOs[et_Neutral] > 0.) {
992
993
        tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_NumChargedOverNumNeutral", num_EFOs[et_Charged] / num_EFOs[et_Neutral]);
    }
994
    if(tlv_Sys_OK[et_Charged] == true && tlv_Sys_OK[et_All] && num_EFOs[et_All] > 0.) {
995
996
997
998
        tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_NumChargedOverNumTotal",   num_EFOs[et_Charged] / num_EFOs[et_All]);
    }
    
    if(num_EFOs[et_Charged]>0. && num_EFOs[et_Neutral]>1.) {
999
1000
1001
1002
        if(tlv_1st_OK[et_Charged] && tlv_1st_OK[et_Neutral] && tlv_2nd_OK[et_Neutral]) {
            TVector3 axis_Plane_cn1 = (tlv_1stEFO[et_Charged].Vect()).Cross( tlv_1stEFO[et_Neutral].Vect() );
            TVector3 axis_Plane_cn2 = (tlv_1stEFO[et_Charged].Vect()).Cross( tlv_2ndEFO[et_Neutral].Vect() );
            double anglePlanes = axis_Plane_cn1.Angle(axis_Plane_cn2);
1003
1004
1005
1006
1007
            tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_AnglePlane1stCharged1st2ndNeutral", anglePlanes);
        }
    }
    
    
1008
    PanTau::TauConstituent2* tauConst_NeutralLargestAngle = m_HelperFunctions.getNeutralConstWithLargestAngle(tlv_System[et_Charged],
1009
1010
                                                                                                                   list_NeutralConstituents);
    if(tauConst_NeutralLargestAngle != 0) {
1011
        TLorentzVector tlv_NeutralLargestAngle = tauConst_NeutralLargestAngle->p4();
1012
        
1013
        tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_FarthestNeutral_AngleToCharged", tlv_System[et_Charged].Angle(tlv_NeutralLargestAngle.Vect()) );
1014
        tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_FarthestNeutral_BDTScore", tauConst_NeutralLargestAngle->getBDTValue());
1015
        if(tlv_System[et_Charged].Et() > 0) tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_FarthestNeutral_EtOverChargedEt", tlv_NeutralLargestAngle.Et() / tlv_System[et_Charged].Et());
1016
1017
1018
1019
1020
1021
1022
    }
    
    //! Combined Type-vs-Type Features ////////////////////////////////////////
    //Loop over all EFO types...
    //  for every type, loop over all other types (but neglect permutations, i.e.   A,B   is the same as   B,A)
    //     calculate ratios, angles etc...
    
1023
    for(int iType=0; iType<PanTau::TauConstituent2::t_nTypes; iType++) {
1024
        
1025
        if(iType == (int)PanTau::TauConstituent2::t_NoType) continue;
1026
1027
        int type_Denom    = iType;
        
1028
        for(int jType=0; jType<PanTau::TauConstituent2::t_nTypes; jType++) {
1029
            
1030
            if(jType == (int)PanTau::TauConstituent2::t_NoType) continue;
1031
1032
1033
1034
            int type_Nom   = jType;
            
            if(jType == iType) continue;
            
1035
1036
            std::string typeName_Nom    = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)type_Nom);
            std::string typeName_Denom  = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)type_Denom);
1037
1038
1039
            
            double sum_Et_Nom   = 0.0;
            double sum_Et_Denom = 0.0;
1040
1041
1042
            if(tlv_Sys_OK[type_Nom] && tlv_Sys_OK[type_Denom]) {
                sum_Et_Nom   = tlv_System[type_Nom].Et();
                sum_Et_Denom = tlv_System[type_Denom].Et();
1043
1044
1045
1046
            }
            
            //Fraction of leading EFO of system A wrt complete system B
            // this is not symmetric wrt system A and B, hence do this for all combinations
1047
1048
            if(tlv_1st_OK[type_Nom]) {
                double LeadEFO_Et_Nom = tlv_1stEFO[type_Nom].Et();
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
                if(LeadEFO_Et_Nom > 0. && sum_Et_Denom > 0.) {
                    tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_Log1st" + typeName_Nom + "EtOver" + typeName_Denom + "Et", TMath::Log10(LeadEFO_Et_Nom / sum_Et_Denom) );
                }
            }
            
            //following features symmetric in system A-B, hence skip multi calculation
            if(jType <= iType) continue;
            
            //Energy ratios between systems
            if(sum_Et_Denom > 0. && sum_Et_Nom > 0.) {
                tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_Log" + typeName_Nom + "EtOver" + typeName_Denom + "Et",   TMath::Log10(sum_Et_Nom / sum_Et_Denom) );
            }//end check for div by zero
            
            //Angles between systems
1063
1064
1065
1066
            if(tlv_Sys_OK[type_Nom] && tlv_Sys_OK[type_Denom]) {
	        double angle_system = tlv_System[type_Nom].Angle( tlv_System[type_Denom].Vect() );
                m_HelperFunctions.dumpFourMomentum(tlv_System[type_Nom]);
                m_HelperFunctions.dumpFourMomentum(tlv_System[type_Denom]);
1067
1068
1069
1070
                tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_Angle" + typeName_Nom + "To" + typeName_Denom, angle_system );
            }//end check for valid system pointer
            
            
1071
            if(tlv_1st_OK[type_Nom] && tlv_1st_OK[type_Denom]) {
1072
                //Delta R between 1st and 1st EFO
1073
	        double deltaR_1st1st = tlv_1stEFO[type_Nom].DeltaR( tlv_1stEFO[type_Denom] );
1074
1075
1076
                tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_DeltaR1st" + typeName_Nom + "To1st" + typeName_Denom, deltaR_1st1st );
                
                //Angles between 1st and 1st EFO
1077
                double angle_1st1st = tlv_1stEFO[type_Nom].Angle( tlv_1stEFO[type_Denom].Vect() );
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
                tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_Angle1st" + typeName_Nom + "To1st" + typeName_Denom, angle_1st1st );
            } //end check for valid leading efo
            
        }//end loop over system B
    }//end loop over system A
    
    
    //! Combined Selected-Type Features ////////////////////////////////////////

    //setup arrays for combination of selected charged and neutral combinations
    const int cTypes = 1;
    const int nTypes = 2;
    int index_charged[cTypes] = {et_Charged};
    int index_neutral[nTypes] = {et_Pi0Neut, et_Neutral};
    
    for(int cType=0; cType<cTypes; cType++) {
        for(int nType=0; nType<nTypes; nType++) {
            
            int et_c = index_charged[cType];
            int et_n = index_neutral[nType];
            
1099
1100
            std::string name_cType = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)et_c);
            std::string name_nType = PanTau::TauConstituent2::getTypeName((PanTau::TauConstituent2::Type)et_n);
1101
            
1102
            if(tlv_Sys_OK[et_c]==false || tlv_Sys_OK[et_n]==false) continue;
1103
1104
1105
            
            //mean Et fraction of charged+neutral system wrt total ET
            if(num_EFOs[et_c] + num_EFOs[et_n] > 0.) {
1106
                double mean_cTypenTypeEt = ( tlv_System[et_c].Et() + tlv_System[et_n].Et() ) / (num_EFOs[et_c] + num_EFOs[et_n]);
1107
1108
1109
1110
                addFeatureWrtSeedEnergy(tauFeatures, inputAlgName + "_" + prefixVARType + "_Mean" + name_cType + name_nType + "Et_Wrt", mean_cTypenTypeEt, &m_Variants_SeedEt);
            }
            
            //invariant masses
1111
            double mass_cTypenType = ( tlv_System[et_c]  +  tlv_System[et_n] ).M();
1112
1113
1114
            tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_InvMass"  + name_cType + name_nType,   mass_cTypenType );
            
            //angle 1st charged to second neutral
1115
            if(tlv_2nd_OK[et_n]) {
1116
                //Angles between 1st and 2nd EFO
1117
	        double angle_1st2nd = tlv_1stEFO[et_c].Angle( tlv_2ndEFO[et_n].Vect() );
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
                tauFeatures->addFeature(inputAlgName + "_" + prefixVARType + "_Angle1st2nd" + name_cType + name_nType, angle_1st2nd );
            } //end check for valid 2nd EFOs
            
        }//end loop neutral types
    }//end loop charged types
    
    
    return StatusCode::SUCCESS;
}



1130
StatusCode PanTau::Tool_FeatureExtractor::addGenericJetFeatures(PanTau::PanTauSeed2* inSeed) const {
1131
1132
    
    std::string                     inputAlgName    = inSeed->getNameInputAlgorithm();
1133
1134
    std::vector<TauConstituent2*>    allConstituents = inSeed->getConstituentsAsList_Core();
    PanTau::TauFeature2*             tauFeatures     = inSeed->getFeatures();
1135
1136
1137
1138
1139
    
    const std::string