4
4
#include " Artus/Utility/interface/DefaultValues.h"
5
5
#include " Artus/Utility/interface/Utility.h"
6
6
7
+ #include < boost/algorithm/string.hpp>
8
+ #include < boost/algorithm/string/trim.hpp>
9
+
7
10
8
11
std::string MELAProducer::GetProducerId () const
9
12
{
@@ -14,12 +17,11 @@ void MELAProducer::Init(setting_type const& settings, metadata_type& metadata)
14
17
{
15
18
ProducerBase<HttTypes>::Init (settings, metadata);
16
19
20
+ m_higgsProductionMode = HttEnumTypes::ToMELAHiggsProductionMode (boost::algorithm::to_lower_copy (boost::algorithm::trim_copy (settings.GetMELAHiggsProductionMode ())));
21
+
17
22
// https://github.com/cms-analysis/HiggsAnalysis-ZZMatrixElement/blob/v2.1.1/MELA/interface/Mela.h
18
23
// https://github.com/cms-analysis/HiggsAnalysis-ZZMatrixElement/blob/v2.1.1/MELA/interface/TVar.hh
19
- if (metadata.m_mela == nullptr )
20
- {
21
- metadata.m_mela = new Mela (13.0 , 125.0 , TVar::SILENT);
22
- }
24
+ m_mela = std::unique_ptr<Mela>(new Mela (13.0 , 125.0 , TVar::SILENT));
23
25
}
24
26
25
27
@@ -30,7 +32,7 @@ void MELAProducer::Produce(event_type const& event, product_type& product,
30
32
{
31
33
SimpleParticleCollection_t daughters; // Higgs boson or two tau leptons
32
34
SimpleParticleCollection_t associated; // additional reconstructed jets
33
- SimpleParticleCollection_t mothers; // incoming partons in case of gen level mode
35
+ // SimpleParticleCollection_t mothers; // incoming partons in case of gen level mode
34
36
35
37
TLorentzVector higgsLV = Utility::ConvertPtEtaPhiMLorentzVector<RMFLV, TLorentzVector>(*product.m_svfitResults .fittedHiggsLV );
36
38
daughters.emplace_back (DefaultValues::pdgIdH, higgsLV);
@@ -44,14 +46,60 @@ void MELAProducer::Produce(event_type const& event, product_type& product,
44
46
}
45
47
}
46
48
47
- metadata.m_mela ->setInputEvent (&daughters, &associated, &mothers, false );
49
+ m_mela->setInputEvent (&daughters, &associated, nullptr /* &mothers*/ , false );
50
+
51
+ // CP even
52
+ if (m_higgsProductionMode == HttEnumTypes::MELAHiggsProductionMode::GGH)
53
+ {
54
+ m_mela->setProcess (TVar::HSMHiggs, TVar::JHUGen, TVar::JJQCD);
55
+ }
56
+ else if (m_higgsProductionMode == HttEnumTypes::MELAHiggsProductionMode::VBF)
57
+ {
58
+ m_mela->setProcess (TVar::HSMHiggs, TVar::JHUGen, TVar::JJVBF);
59
+ }
60
+ float probCPEven = 0.0 ;
61
+ m_mela->computeProdP (probCPEven, false );
62
+
63
+ // CP odd
64
+ if (m_higgsProductionMode == HttEnumTypes::MELAHiggsProductionMode::GGH)
65
+ {
66
+ m_mela->setProcess (TVar::H0minus, TVar::JHUGen, TVar::JJQCD);
67
+ }
68
+ else if (m_higgsProductionMode == HttEnumTypes::MELAHiggsProductionMode::VBF)
69
+ {
70
+ m_mela->setProcess (TVar::H0minus, TVar::JHUGen, TVar::JJVBF);
71
+ }
72
+ float probCPOdd = 0.0 ;
73
+ m_mela->computeProdP (probCPOdd, false );
74
+
75
+ // CP mixing (maximum)
76
+ if (m_higgsProductionMode == HttEnumTypes::MELAHiggsProductionMode::GGH)
77
+ {
78
+ m_mela->setProcess (TVar::SelfDefine_spin0, TVar::JHUGen, TVar::JJQCD);
79
+ m_mela->selfDHggcoupl [0 ][gHIGGS_GG_2 ][0 ] = 1 ; // a1
80
+ m_mela->selfDHggcoupl [0 ][gHIGGS_GG_4 ][0 ] = 1 ; // a3
81
+ }
82
+ else if (m_higgsProductionMode == HttEnumTypes::MELAHiggsProductionMode::VBF)
83
+ {
84
+ m_mela->setProcess (TVar::SelfDefine_spin0, TVar::JHUGen, TVar::JJVBF);
85
+ m_mela->selfDHzzcoupl [0 ][gHIGGS_VV_1 ][0 ] = 1 ; // a1
86
+ m_mela->selfDHzzcoupl [0 ][gHIGGS_VV_4 ][0 ] = 0.297979 ; // a3
87
+ }
88
+ float probCPMix = 0.0 ;
89
+ m_mela->computeProdP (probCPMix, false );
90
+
91
+ LOG (WARNING) << " probabilities: " << probCPEven << " , " << probCPOdd << " , " << probCPMix;
48
92
49
- metadata.m_mela ->setProcess (TVar::HSMHiggs, TVar::JHUGen, TVar::ZZGG); // TODO
50
- float result = 0.0 ;
51
- metadata.m_mela ->computeProdP (result, true );
52
- LOG (INFO) << " production matrix element: " << result;
93
+ float discriminatorD0minus = DefaultValues::UndefinedFloat;
94
+ float discriminatorDCP = DefaultValues::UndefinedFloat;
95
+ if ((probCPEven + probCPOdd) != 0.0 )
96
+ {
97
+ discriminatorD0minus = probCPEven / (probCPEven + probCPOdd) if ;
98
+ discriminatorDCP = (probCPMix - probCPEven - probCPOdd) / (probCPEven + probCPOdd);
99
+ }
100
+ LOG (WARNING) << " discriminators: " << discriminatorD0minus << " , " << discriminatorDCP;
53
101
54
- metadata. m_mela ->resetInputEvent ();
102
+ m_mela->resetInputEvent ();
55
103
}
56
104
}
57
105
0 commit comments