|
|
|
#include <TtHadLRSignalSelObservables.h>
Classes | |
| struct | CompareBdisc |
| struct | CompareDouble |
| struct | CompareET |
Public Member Functions | |
| void | operator() (TtHadEvtSolution &) |
| TtHadLRSignalSelObservables () | |
| ~TtHadLRSignalSelObservables () | |
Private Attributes | |
| CompareBdisc | BdiscComparator |
| CompareDouble | dComparator |
| CompareET | EtComparator |
| std::vector< std::pair < unsigned int, double > > | evtselectVarVal |
Definition at line 26 of file TtHadLRSignalSelObservables.h.
| TtHadLRSignalSelObservables::TtHadLRSignalSelObservables | ( | ) |
Definition at line 3 of file TtHadLRSignalSelObservables.cc.
{
}
| TtHadLRSignalSelObservables::~TtHadLRSignalSelObservables | ( | ) |
Definition at line 7 of file TtHadLRSignalSelObservables.cc.
{
}
| void TtHadLRSignalSelObservables::operator() | ( | TtHadEvtSolution & | TS | ) |
Definition at line 11 of file TtHadLRSignalSelObservables.cc.
References BdiscComparator, funct::C, funct::cos(), funct::D, dComparator, relval_parameters_module::energy, EtComparator, evtselectVarVal, TtHadEvtSolution::getHadb(), TtHadEvtSolution::getHadbbar(), TtHadEvtSolution::getHadj(), TtHadEvtSolution::getHadk(), TtHadEvtSolution::getHadp(), TtHadEvtSolution::getHadq(), i, edm::detail::isnan(), funct::log(), WElectronSkim_cff::MT, MultiGaussianStateTransform::N, PI, funct::pow(), ExpressReco_HICollisions_FallBack::pt, TtHadEvtSolution::setLRSignalEvtObservables(), funct::sin(), python::multivaluedict::sort(), and mathSSE::sqrt().
{
evtselectVarVal.clear();
std::vector<pat::Jet> topJets;
topJets.clear();
topJets.push_back(TS.getHadp());
topJets.push_back(TS.getHadq());
topJets.push_back(TS.getHadj());
topJets.push_back(TS.getHadk());
topJets.push_back(TS.getHadb());
topJets.push_back(TS.getHadbbar());
//Assume the highest Et jets are the b-jets, the others the jets from the W - but how to work out which two are from which W? FIXME!!! for now don't sort jets in Et before after this...
TLorentzVector *Hadp = new TLorentzVector();
TLorentzVector *Hadq = new TLorentzVector();
TLorentzVector *Hadj = new TLorentzVector();
TLorentzVector *Hadk = new TLorentzVector();
TLorentzVector *Hadb = new TLorentzVector();
TLorentzVector *Hadbbar = new TLorentzVector();
Hadp->SetPxPyPzE(topJets[0].px(),topJets[0].py(),topJets[0].pz(),topJets[3].energy());
Hadq->SetPxPyPzE(topJets[1].px(),topJets[1].py(),topJets[1].pz(),topJets[2].energy());
Hadj->SetPxPyPzE(topJets[2].px(),topJets[2].py(),topJets[2].pz(),topJets[2].energy());
Hadk->SetPxPyPzE(topJets[3].px(),topJets[3].py(),topJets[3].pz(),topJets[3].energy());
Hadb->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].energy());
Hadbbar->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].energy());
//sort the topJets in Et
std::sort(topJets.begin(),topJets.end(),EtComparator);
//Et-Sum of the lightest jets
double EtSum = topJets[5].et()+topJets[5].et();
double Obs1 = (EtSum>0 ? EtSum : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(1,Obs1));
//Log of the difference in Bdisc between the 2nd and the 3rd jets (ordered in Bdisc) - does that still
//make sense for the fully hadronic channel as well?FIXME!!!
//sort the jets in Bdiscriminant
std::sort(topJets.begin(),topJets.end(),BdiscComparator);
double BGap = topJets[1].bDiscriminator("trackCountingJetTags") - topJets[2].bDiscriminator("trackCountingJetTags");
double Obs2 = (BGap>0 ? log(BGap) : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(2,Obs2));
//Circularity of the event
double N=0,D=0,C_tmp=0,C=1000;
double nx,ny,nz;
// C = 2min(E(pt.n)^2/E(pt)^2) = 2*N/D but it is theorically preferable to use C'=PI/2*min(E|pt.n|/E|pt|), sum over all jets+lepton+MET (cf PhysRevD 48 R3953(Nov 1993))
for(unsigned int i=0;i<6;i++){
D += fabs(topJets[i].pt());
}
if((D>0)){
// Loop over all the unit vectors in the transverse plane in order to find the miminum :
for(unsigned int i=0; i<360; i++){
nx = cos((2*PI/360)*i);
ny = sin((2*PI/360)*i);
nz = 0;
N=0;
for(unsigned int i=0;i<4;i++){
N += fabs(topJets[i].px()*nx+topJets[i].py()*ny+topJets[i].pz()*nz);
}
C_tmp = 2*N/D;
if(C_tmp<C) C = C_tmp;
}
}
double Obs3 = ( C!=1000 ? C : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(3,Obs3));
//HT variable (Et-sum of the six jets)
double HT=0;
for(unsigned int i=0;i<6;i++){
HT += topJets[i].et();
}
double Obs4 = ( HT!=0 ? HT : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(4,Obs4));
//Transverse Mass of the system
math::XYZTLorentzVector pjets;
// for the six jets
for(unsigned int i=0;i<6;i++){
pjets += topJets[i].p4();
}
double MT = sqrt(std::pow(pjets.mass(),2));
double Obs5 = ( MT>0 ? MT : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(5,Obs5));
//CosTheta(Hadp,Hadq) and CosTheta(Hadj,Hadq)
//sort the lightJets in Et
std::sort(topJets.begin(),topJets.end(),EtComparator);
double px1 = topJets[2].px(); double px2 = topJets[3].px();
double py1 = topJets[2].py(); double py2 = topJets[3].py();
double pz1 = topJets[2].pz(); double pz2 = topJets[3].pz();
double E1 = topJets[2].energy(); double E2 = topJets[3].energy();
double px3 = topJets[4].px(); double px4 = topJets[5].px();
double py3 = topJets[4].py(); double py4 = topJets[5].py();
double pz3 = topJets[4].pz(); double pz4 = topJets[5].pz();
double E3 = topJets[4].energy(); double E4 = topJets[5].energy();
TLorentzVector *LightJet1 = new TLorentzVector();
TLorentzVector *LightJet2 = new TLorentzVector();
TLorentzVector *LightJet3 = new TLorentzVector();
TLorentzVector *LightJet4 = new TLorentzVector();
LightJet1->SetPxPyPzE(px1,py1,pz1,E1);
//LightJet1->Boost(BoostBackToCM);
LightJet2->SetPxPyPzE(px2,py2,pz2,E2);
//LightJet2->Boost(BoostBackToCM);
LightJet3->SetPxPyPzE(px3,py3,pz3,E3);
//LightJet1->Boost(BoostBackToCM);
LightJet4->SetPxPyPzE(px4,py4,pz4,E4);
//LightJet2->Boost(BoostBackToCM);
double CosTheta1 = cos(LightJet2->Angle(LightJet1->Vect()));
double CosTheta2 = cos(LightJet4->Angle(LightJet3->Vect()));
double Obs6 = ( -1<CosTheta1 ? CosTheta1 : -2);
evtselectVarVal.push_back(std::pair<unsigned int,double>(6,Obs6));
double Obs7 = ( -1<CosTheta2 ? CosTheta2 : -2);
evtselectVarVal.push_back(std::pair<unsigned int,double>(7,Obs7));
delete LightJet1;
delete LightJet2;
delete LightJet3;
delete LightJet4;
// try to find out more powerful observables related to the b-disc
//sort the jets in Bdiscriminant
std::sort(topJets.begin(),topJets.end(),BdiscComparator);
double BjetsBdiscSum = topJets[0].bDiscriminator("trackCountingJetTags") + topJets[1].bDiscriminator("trackCountingJetTags");
double Ljets1BdiscSum = topJets[2].bDiscriminator("trackCountingJetTags") + topJets[3].bDiscriminator("trackCountingJetTags");
double Ljets2BdiscSum = topJets[4].bDiscriminator("trackCountingJetTags") + topJets[5].bDiscriminator("trackCountingJetTags");
//std::cout<<"BjetsBdiscSum = "<<BjetsBdiscSum<<std::endl;
//std::cout<<"LjetsBdiscSum = "<<LjetsBdiscSum<<std::endl;
double Obs8 = (Ljets1BdiscSum !=0 ? log(BjetsBdiscSum/Ljets1BdiscSum) : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(8,Obs8));
double Obs9 = (Ljets2BdiscSum !=0 ? log(BjetsBdiscSum/Ljets2BdiscSum) : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(9,Obs9));
double Obs10 = (BGap>0 ? log(BjetsBdiscSum*BGap) : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(10,Obs10));
// Et-Ratio between the two highest in Et jets and four highest jets
double Obs11 = (HT!=0 ? (HT-EtSum)/HT : -1);
evtselectVarVal.push_back(std::pair<unsigned int,double>(11,Obs11));
//Sphericity and Aplanarity without boosting back the system to CM frame
TMatrixDSym Matrix(3);
double PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+
std::pow(Hadj->Px(),2)+std::pow(Hadk->Px(),2)+std::pow(Hadbbar->Px(),2);
double PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+
std::pow(Hadj->Py(),2)+std::pow(Hadk->Py(),2)+std::pow(Hadbbar->Py(),2);
double PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+
std::pow(Hadj->Pz(),2)+std::pow(Hadk->Pz(),2)+std::pow(Hadbbar->Pz(),2);
double P2 = PX2+PY2+PZ2;
double PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+
Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py();
double PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+
Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz();
double PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+
Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz();
Matrix(0,0) = PX2/P2; Matrix(0,1) = PXY/P2; Matrix(0,2) = PXZ/P2;
Matrix(1,0) = PXY/P2; Matrix(1,1) = PY2/P2; Matrix(1,2) = PYZ/P2;
Matrix(2,0) = PXZ/P2; Matrix(2,1) = PYZ/P2; Matrix(2,2) = PZ2/P2;
TMatrixDSymEigen pTensor(Matrix);
std::vector<double> EigValues;
EigValues.clear();
for(int i=0;i<3;i++){
EigValues.push_back(pTensor.GetEigenValues()[i]);
}
std::sort(EigValues.begin(),EigValues.end(),dComparator);
double Sphericity = 1.5*(EigValues[1]+EigValues[2]);
double Aplanarity = 1.5*EigValues[2];
double Obs12 = (isnan(Sphericity) ? -1 : Sphericity);
evtselectVarVal.push_back(std::pair<unsigned int,double>(12,Obs12));
double Obs13 = (isnan(Aplanarity) ? -1 : Aplanarity);
evtselectVarVal.push_back(std::pair<unsigned int,double>(13,Obs13));
//Sphericity and Aplanarity with boosting back the system to CM frame
TLorentzVector *TtbarSystem = new TLorentzVector();
TtbarSystem->SetPx(Hadp->Px()+Hadq->Px()+Hadb->Px()+Hadj->Px()+Hadk->Px()+Hadbbar->Px());
TtbarSystem->SetPy(Hadp->Py()+Hadq->Py()+Hadb->Py()+Hadj->Py()+Hadk->Py()+Hadbbar->Py());
TtbarSystem->SetPz(Hadp->Pz()+Hadq->Pz()+Hadb->Pz()+Hadj->Pz()+Hadk->Pz()+Hadbbar->Pz());
TtbarSystem->SetE(Hadp->E()+Hadq->E()+Hadb->E()+Hadj->E()+Hadk->E()+Hadbbar->E());
TVector3 BoostBackToCM = -(TtbarSystem->BoostVector());
Hadp->Boost(BoostBackToCM);
Hadq->Boost(BoostBackToCM);
Hadb->Boost(BoostBackToCM);
Hadj->Boost(BoostBackToCM);
Hadk->Boost(BoostBackToCM);
Hadbbar->Boost(BoostBackToCM);
double BOOST_PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+
std::pow(Hadj->Px(),2)+std::pow(Hadk->Px(),2)+std::pow(Hadbbar->Px(),2);
double BOOST_PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+
std::pow(Hadj->Py(),2)+std::pow(Hadk->Py(),2)+std::pow(Hadbbar->Py(),2);
double BOOST_PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+
std::pow(Hadj->Pz(),2)+std::pow(Hadk->Pz(),2)+std::pow(Hadbbar->Pz(),2);
double BOOST_P2 = BOOST_PX2+BOOST_PY2+BOOST_PZ2;
double BOOST_PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+
Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py();
double BOOST_PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+
Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz();
double BOOST_PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+
Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz();
TMatrixDSym BOOST_Matrix(3);
BOOST_Matrix(0,0) = BOOST_PX2/BOOST_P2; BOOST_Matrix(0,1) = BOOST_PXY/BOOST_P2; BOOST_Matrix(0,2) = BOOST_PXZ/BOOST_P2;
BOOST_Matrix(1,0) = BOOST_PXY/BOOST_P2; BOOST_Matrix(1,1) = BOOST_PY2/BOOST_P2; BOOST_Matrix(1,2) = BOOST_PYZ/BOOST_P2;
BOOST_Matrix(2,0) = BOOST_PXZ/BOOST_P2; BOOST_Matrix(2,1) = BOOST_PYZ/BOOST_P2; BOOST_Matrix(2,2) = BOOST_PZ2/BOOST_P2;
TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
std::vector<double> BOOST_EigValues;
BOOST_EigValues.clear();
for(int i=0;i<3;i++){
BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]);
}
std::sort(BOOST_EigValues.begin(),BOOST_EigValues.end(),dComparator);
double BOOST_Sphericity = 1.5*(BOOST_EigValues[1]+BOOST_EigValues[2]);
double BOOST_Aplanarity = 1.5*BOOST_EigValues[2];
double Obs14 = ( isnan(BOOST_Sphericity) ? -1 : BOOST_Sphericity );
evtselectVarVal.push_back(std::pair<unsigned int,double>(14,Obs14));
double Obs15 = ( isnan(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity );
evtselectVarVal.push_back(std::pair<unsigned int,double>(15,Obs15));
// Centrality of the six jets
double H=0;
for(unsigned int i=0;i<6;i++){
H += topJets[i].energy();
}
double Obs16 = ( H != 0 ? HT/H : -1 );
evtselectVarVal.push_back(std::pair<unsigned int,double>(16,Obs16));
// distance from the origin in the (BjetsBdiscSum, Ljets1BdiscSum) and (BjetsBdiscSum, Ljets2BdiscSum)
double Obs17 = ( BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets1BdiscSum) : -1 );
evtselectVarVal.push_back(std::pair<unsigned int,double>(17,Obs17));
double Obs18 = ( BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets2BdiscSum) : -1 );
evtselectVarVal.push_back(std::pair<unsigned int,double>(18,Obs18));
// Ratio of the Et Sum of the two lightest jets over HT=Sum of the Et of the six highest jets in Et
double Obs19 = ( HT != 0 ? EtSum/HT : -1 );
evtselectVarVal.push_back(std::pair<unsigned int,double>(19,Obs19));
// Put the vector in the TtHadEvtSolution
TS.setLRSignalEvtObservables(evtselectVarVal);
}
Definition at line 55 of file TtHadLRSignalSelObservables.h.
Referenced by operator()().
Definition at line 65 of file TtHadLRSignalSelObservables.h.
Referenced by operator()().
Definition at line 45 of file TtHadLRSignalSelObservables.h.
Referenced by operator()().
std::vector<std::pair<unsigned int,double> > TtHadLRSignalSelObservables::evtselectVarVal [private] |
Definition at line 67 of file TtHadLRSignalSelObservables.h.
Referenced by operator()().
1.7.3