|
|
|
#include <TtSemiLRSignalSelObservables.h>
Classes | |
| struct | CompareBdisc |
| struct | CompareDouble |
| struct | CompareET |
Public Member Functions | |
| void | operator() (TtSemiEvtSolution &, const std::vector< pat::Jet > &) |
| TtSemiLRSignalSelObservables () | |
| ~TtSemiLRSignalSelObservables () | |
Private Attributes | |
| CompareBdisc | BdiscComparator |
| CompareDouble | dComparator |
| CompareET | EtComparator |
| std::vector< std::pair < unsigned int, double > > | evtselectVarVal |
Definition at line 29 of file TtSemiLRSignalSelObservables.h.
| TtSemiLRSignalSelObservables::TtSemiLRSignalSelObservables | ( | ) |
Definition at line 3 of file TtSemiLRSignalSelObservables.cc.
{
}
| TtSemiLRSignalSelObservables::~TtSemiLRSignalSelObservables | ( | ) |
Definition at line 7 of file TtSemiLRSignalSelObservables.cc.
{
}
| void TtSemiLRSignalSelObservables::operator() | ( | TtSemiEvtSolution & | TS, |
| const std::vector< pat::Jet > & | SelectedTopJets | ||
| ) |
Definition at line 11 of file TtSemiLRSignalSelObservables.cc.
References BdiscComparator, funct::C, MEzCalculator::Calculate(), funct::cos(), gather_cfg::cout, funct::D, dComparator, DEBUG, reco::LeafCandidate::energy(), relval_parameters_module::energy, reco::LeafCandidate::et(), ExpressReco_HICollisions_FallBack::et, EtComparator, evtselectVarVal, TtSemiEvtSolution::getDecay(), TtSemiEvtSolution::getElectron(), TtSemiEvtSolution::getMuon(), TtSemiEvtSolution::getNeutrino(), i, edm::detail::isnan(), j, gen::k, funct::log(), WElectronSkim_cff::MT, MultiGaussianStateTransform::N, n, p4, PI, funct::pow(), ExpressReco_HICollisions_FallBack::pt, reco::LeafCandidate::px(), reco::LeafCandidate::py(), reco::LeafCandidate::pz(), MEzCalculator::SetLepton(), TtSemiEvtSolution::setLRSignalEvtObservables(), MEzCalculator::SetMET(), funct::sin(), and python::multivaluedict::sort().
{
evtselectVarVal.clear();
// choice the B-tag algorithm
const char *BtagAlgo = "trackCountingJetTags";
// activate the "debug mode"
bool DEBUG = false;
if(DEBUG) std::cout<<"---------- Start calculating the LR observables ----------"<<std::endl;
std::vector<pat::Jet> TopJets;
TopJets.clear();
for(unsigned int i = 0 ; i<SelectedTopJets.size() ; i++)
{
TopJets.push_back(SelectedTopJets[i]);
}
//sort the TopJets in Et
std::sort(TopJets.begin(),TopJets.end(),EtComparator);
TLorentzVector * Hadp = new TLorentzVector();
Hadp->SetPxPyPzE(TopJets[3].px(),TopJets[3].py(),TopJets[3].pz(),TopJets[3].energy());
TLorentzVector * Hadq = new TLorentzVector();
Hadq->SetPxPyPzE(TopJets[2].px(),TopJets[2].py(),TopJets[2].pz(),TopJets[2].energy());
TLorentzVector * Hadb = new TLorentzVector();
Hadb->SetPxPyPzE(TopJets[1].px(),TopJets[1].py(),TopJets[1].pz(),TopJets[1].energy());
TLorentzVector * Lepb = new TLorentzVector();
Lepb->SetPxPyPzE(TopJets[0].px(),TopJets[0].py(),TopJets[0].pz(),TopJets[0].energy());
TLorentzVector * Lept = new TLorentzVector();
if(TS.getDecay() == "muon") Lept->SetPxPyPzE(TS.getMuon().px(),TS.getMuon().py(),TS.getMuon().pz(),TS.getMuon().energy());
else if(TS.getDecay() == "electron") Lept->SetPxPyPzE(TS.getElectron().px(),TS.getElectron().py(),TS.getElectron().pz(),TS.getElectron().energy());
TLorentzVector *Lepn = new TLorentzVector();
Lepn->SetPxPyPzE(TS.getNeutrino().px(),TS.getNeutrino().py(),TS.getNeutrino().pz(),TS.getNeutrino().energy());
// Calculation of the pz of the neutrino due to W-mass constraint
MEzCalculator *Mez = new MEzCalculator();
Mez->SetMET(TS.getNeutrino());
if(TS.getDecay() == "muon") Mez->SetLepton(TS.getMuon());
else Mez->SetLepton(TS.getElectron(), false);
double MEZ = Mez->Calculate();
Lepn->SetPz(MEZ);
// Pt of the lepton
double Obs1 = Lept->Pt();
evtselectVarVal.push_back(std::pair<unsigned int,double>(1,Obs1));
if(DEBUG) std::cout<<"------ LR observable 1 "<<Obs1<<" calculated ------"<<std::endl;
// Missing transverse energy
double Obs2 = TS.getNeutrino().et();
evtselectVarVal.push_back(std::pair<unsigned int,double>(2,Obs2));
if(DEBUG) std::cout<<"------ LR observable 2 "<<Obs2<<" calculated ------"<<std::endl;
//HT variable (Et-sum of the four jets)
double HT=0;
for(unsigned int i=0;i<4;i++)
{
HT += TopJets[i].et();
}
double Obs3 = HT;
evtselectVarVal.push_back(std::pair<unsigned int,double>(3,Obs3));
if(DEBUG) std::cout<<"------ LR observable 3 "<<Obs3<<" calculated ------"<<std::endl;
//Et-Sum of the lightest jets
double EtSum = TopJets[2].et()+TopJets[3].et();
double Obs4 = EtSum;
evtselectVarVal.push_back(std::pair<unsigned int,double>(4,Obs4));
if(DEBUG) std::cout<<"------ LR observable 4 "<<Obs4<<" calculated ------"<<std::endl;
// Et-Ratio between the two lowest jets in Et and four highest jets
double Obs5 = EtSum/HT;
evtselectVarVal.push_back(std::pair<unsigned int,double>(5,Obs5));
if(DEBUG) std::cout<<"------ LR observable 5 "<<Obs5<<" calculated ------"<<std::endl;
// Et-Ratio between the two highest jets in Et and four highest jets
double Obs6 = (HT-EtSum)/HT;
evtselectVarVal.push_back(std::pair<unsigned int,double>(6,Obs6));
if(DEBUG) std::cout<<"------ LR observable 6 "<<Obs6<<" calculated ------"<<std::endl;
// Transverse Mass of the system
TLorentzVector TtbarSystem = (*Hadp)+(*Hadq)+(*Hadb)+(*Lepb)+(*Lept)+(*Lepn);
double MT = TtbarSystem.Mt();
double Obs7 = MT;
evtselectVarVal.push_back(std::pair<unsigned int,double>(7,Obs7));
if(DEBUG) std::cout<<"------ LR observable 7 "<<Obs7<<" calculated ------"<<std::endl;
// Observables related to the b-disc (jets ordered in Bdisc)
//sort the TopJets in Bdiscriminant
std::sort(TopJets.begin(),TopJets.end(),BdiscComparator);
double Obs8;
double Obs9;
double Obs10;
double Obs11;
// Difference in bdisc between the 2nd and the 3rd jets
double BGap = TopJets[1].bDiscriminator(BtagAlgo) - TopJets[2].bDiscriminator(BtagAlgo);
// Sum of the bdisc of the two highest/lowest jets
double BjetsBdiscSum = TopJets[0].bDiscriminator(BtagAlgo) + TopJets[1].bDiscriminator(BtagAlgo);
double LjetsBdiscSum = TopJets[2].bDiscriminator(BtagAlgo) + TopJets[3].bDiscriminator(BtagAlgo);
Obs8 = (TopJets[2].bDiscriminator(BtagAlgo) > -10 ? log(BGap) : -5);
if(DEBUG) std::cout<<"------ LR observable 8 "<<Obs8<<" calculated ------"<<std::endl;
Obs9 = (BjetsBdiscSum*BGap);
if(DEBUG) std::cout<<"------ LR observable 9 "<<Obs9<<" calculated ------"<<std::endl;
Obs10 = (BjetsBdiscSum/LjetsBdiscSum);
if(DEBUG) std::cout<<"------ LR observable 10 "<<Obs10<<" calculated ------"<<std::endl;
// Distance from the origin in the (BjetsBdiscSum, LjetsBdiscSum) plane
Obs11 = 0.707*((BjetsBdiscSum+LjetsBdiscSum)/2 +10);
if(DEBUG) std::cout<<"------ LR observable 11 "<<Obs11<<" calculated ------"<<std::endl;
evtselectVarVal.push_back(std::pair<unsigned int,double>(8,Obs8));
evtselectVarVal.push_back(std::pair<unsigned int,double>(9,Obs9));
evtselectVarVal.push_back(std::pair<unsigned int,double>(10,Obs10));
evtselectVarVal.push_back(std::pair<unsigned int,double>(11,Obs11));
//sort the TopJets in Et
std::sort(TopJets.begin(),TopJets.end(),EtComparator);
// Circularity of the event
double N=0,D=0,C_tmp=0,C=10000;
double N_NoNu=0,D_NoNu=0,C_NoNu_tmp=0,C_NoNu=10000;
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<4;i++)
{
D += fabs(TopJets[i].pt());
}
// modified the 26th of September to calculate also the circularity without the neutrino contribution
D += fabs(Lept->Pt());
D_NoNu = D;
D += fabs(Lepn->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<720; i++)
{
nx = cos((2*PI/720)*i);
ny = sin((2*PI/720)*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);
}
// modified the 26th of September to calculate also the circularity without the neutrino contribution
N += fabs(Lept->Px()*nx+Lept->Py()*ny+Lept->Pz()*nz);
N_NoNu = N;
N += fabs(Lepn->Px()*nx+Lepn->Py()*ny+Lepn->Pz()*nz);
C_tmp = 2*N/D;
// modified the 26th of September to calculate also the circularity without the neutrino contribution
C_NoNu_tmp = 2*N_NoNu/D_NoNu;
if(C_tmp<C) C = C_tmp;
// modified the 26th of September to calculate also the circularity without the neutrino contribution
if(C_NoNu_tmp<C_NoNu) C_NoNu = C_NoNu_tmp;
}
}
double Obs12 = ( C!=10000 ? C : 0);
evtselectVarVal.push_back(std::pair<unsigned int,double>(12,Obs12));
if(DEBUG) std::cout<<"------ LR observable 12 "<<Obs12<<" calculated ------"<<std::endl;
// Circularity of the event without neutrino
double Obs13 = ( C_NoNu != 10000 ? C_NoNu : 0);
evtselectVarVal.push_back(std::pair<unsigned int,double>(13,Obs13));
if(DEBUG) std::cout<<"------ LR observable 13 "<<Obs13<<" calculated ------"<<std::endl;
// Centrality of the four highest jets
double H=0;
for(unsigned int i=0;i<4;i++)
{
H += TopJets[i].energy();
}
double Obs14 = HT/H;
evtselectVarVal.push_back(std::pair<unsigned int,double>(14,Obs14));
if(DEBUG) std::cout<<"------ LR observable 14 "<<Obs14<<" calculated ------"<<std::endl;
// 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(Lepb->Px(),2)+std::pow(Lept->Px(),2)+std::pow(Lepn->Px(),2);
double PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2)+std::pow(Lept->Py(),2)+std::pow(Lepn->Py(),2);
double PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2)+std::pow(Lept->Pz(),2)+std::pow(Lepn->Pz(),2);
double P2 = PX2+PY2+PZ2;
double PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py()+Lept->Px()*Lept->Py()+Lepn->Px()*Lepn->Py();
double PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz()+Lept->Px()*Lept->Pz()+Lepn->Px()*Lepn->Pz();
double PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz()+Lept->Py()*Lept->Pz()+Lepn->Py()*Lepn->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 Obs15 = (isnan(Sphericity) ? 0 : Sphericity);
evtselectVarVal.push_back(std::pair<unsigned int,double>(15,Obs15));
if(DEBUG) std::cout<<"------ LR observable 15 "<<Obs15<<" calculated ------"<<std::endl;
double Obs16 = (isnan(Aplanarity) ? 0 : Aplanarity);
evtselectVarVal.push_back(std::pair<unsigned int,double>(16,Obs16));
if(DEBUG) std::cout<<"------ LR observable 16 "<<Obs16<<" calculated ------"<<std::endl;
// Sphericity and Aplanarity with boosting back the system to CM frame
TVector3 BoostBackToCM = -(TtbarSystem.BoostVector());
Hadp->Boost(BoostBackToCM);
Hadq->Boost(BoostBackToCM);
Hadb->Boost(BoostBackToCM);
Lepb->Boost(BoostBackToCM);
Lept->Boost(BoostBackToCM);
Lepn->Boost(BoostBackToCM);
double BOOST_PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2)+std::pow(Lept->Px(),2)+std::pow(Lepn->Px(),2);
double BOOST_PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2)+std::pow(Lept->Py(),2)+std::pow(Lepn->Py(),2);
double BOOST_PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2)+std::pow(Lept->Pz(),2)+std::pow(Lepn->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()+Lepb->Px()*Lepb->Py()+Lept->Px()*Lept->Py()+Lepn->Px()*Lepn->Py();
double BOOST_PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz()+Lept->Px()*Lept->Pz()+Lepn->Px()*Lepn->Pz();
double BOOST_PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz()+Lept->Py()*Lept->Pz()+Lepn->Py()*Lepn->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 Obs17 = ( isnan(BOOST_Sphericity) ? 0 : BOOST_Sphericity );
evtselectVarVal.push_back(std::pair<unsigned int,double>(17,Obs17));
if(DEBUG) std::cout<<"------ LR observable 17 "<<Obs17<<" calculated ------"<<std::endl;
double Obs18 = ( isnan(BOOST_Aplanarity) ? 0 : BOOST_Aplanarity );
evtselectVarVal.push_back(std::pair<unsigned int,double>(18,Obs18));
if(DEBUG) std::cout<<"------ LR observable 18 "<<Obs18<<" calculated ------"<<std::endl;
//ratio between ET of the fifth jet and HT
double Obs19 = (TopJets.size() > 4) ? TopJets[4].et()/HT : 1.0;
evtselectVarVal.push_back(std::pair<unsigned int,double>(19,Obs19));
if(DEBUG) std::cout<<"------ LR observable 19 "<<Obs19<<" calculated ------"<<std::endl;
// HT variable calculated with all the jets in the event.
double HT_alljets = 0;
double H_alljets = 0;
for(unsigned int i=0;i<TopJets.size();i++)
{
HT_alljets += TopJets[i].et();
H_alljets += TopJets[i].energy();
}
double Obs20 = HT_alljets;
evtselectVarVal.push_back(std::pair<unsigned int,double>(20,Obs20));
if(DEBUG) std::cout<<"------ LR observable 20 "<<Obs20<<" calculated ------"<<std::endl;
// HT3 = HT calculated with all jets except the two leading jets
double HT3_alljets = 0;
for(unsigned int i=2;i<TopJets.size();i++)
{
HT3_alljets += TopJets[i].et();
}
double Obs21 = HT3_alljets;
evtselectVarVal.push_back(std::pair<unsigned int,double>(21,Obs21));
if(DEBUG) std::cout<<"------ LR observable 21 "<<Obs21<<" calculated ------"<<std::endl;
// ET0, ratio of the Et of the leading and HT_alljets
double ET0 = TopJets[0].et()/HT_alljets;
double Obs22 = ET0;
evtselectVarVal.push_back(std::pair<unsigned int,double>(22,Obs22));
if(DEBUG) std::cout<<"------ LR observable 22 "<<Obs22<<" calculated ------"<<std::endl;
// Centrality of the event computed with all jets
double Obs23 = ( (H_alljets>0) ? HT_alljets/H_alljets : 0);
evtselectVarVal.push_back(std::pair<unsigned int,double>(23,Obs23));
if(DEBUG) std::cout<<"------ LR observable 23 "<<Obs23<<" calculated ------"<<std::endl;
//Fox-Wolfram momenta (1st to 6th), modified for hadron collider and using a Legendre polynomials expansion
double FW_momentum_0=0, FW_momentum_1=0, FW_momentum_2=0, FW_momentum_3=0, FW_momentum_4=0, FW_momentum_5=0, FW_momentum_6=0;
for(unsigned int i=0;i<TopJets.size();i++)
{
for(unsigned int j=0;j<TopJets.size();j++)
{
double ET_ij_over_ETSum2= TopJets[i].et()*TopJets[j].et()/(std::pow(HT_alljets,2));
double cosTheta_ij = (TopJets[i].px()*TopJets[j].px()+
TopJets[i].py()*TopJets[j].py()+
TopJets[i].pz()*TopJets[j].pz())
/(TopJets[i].p4().R()*TopJets[j].p4().R());
FW_momentum_0 += ET_ij_over_ETSum2;
FW_momentum_1 += ET_ij_over_ETSum2 * cosTheta_ij;
FW_momentum_2 += ET_ij_over_ETSum2 * 0.5 * ( 3*std::pow(cosTheta_ij,2)- 1);
FW_momentum_3 += ET_ij_over_ETSum2 * 0.5 * ( 5*std::pow(cosTheta_ij,3)- 3*cosTheta_ij);
FW_momentum_4 += ET_ij_over_ETSum2 * 0.125 * ( 35*std::pow(cosTheta_ij,4)- 30*std::pow(cosTheta_ij,2)+3);
FW_momentum_5 += ET_ij_over_ETSum2 * 0.125 * ( 63*std::pow(cosTheta_ij,5)- 70*std::pow(cosTheta_ij,3)+15*cosTheta_ij);
FW_momentum_6 += ET_ij_over_ETSum2 * 0.0625* (231*std::pow(cosTheta_ij,6)-315*std::pow(cosTheta_ij,4)+105*std::pow(cosTheta_ij,2)-5);
}
}
double Obs24 = FW_momentum_0;
evtselectVarVal.push_back(std::pair<unsigned int,double>(24,Obs24));
if(DEBUG) std::cout<<"------ LR observable 24 "<<Obs24<<" calculated ------"<<std::endl;
double Obs25 = FW_momentum_1;
evtselectVarVal.push_back(std::pair<unsigned int,double>(25,Obs25));
if(DEBUG) std::cout<<"------ LR observable 25 "<<Obs25<<" calculated ------"<<std::endl;
double Obs26 = FW_momentum_2;
evtselectVarVal.push_back(std::pair<unsigned int,double>(26,Obs26));
if(DEBUG) std::cout<<"------ LR observable 26 "<<Obs26<<" calculated ------"<<std::endl;
double Obs27 = FW_momentum_3;
evtselectVarVal.push_back(std::pair<unsigned int,double>(27,Obs27));
if(DEBUG) std::cout<<"------ LR observable 27 "<<Obs27<<" calculated ------"<<std::endl;
double Obs28 = FW_momentum_4;
evtselectVarVal.push_back(std::pair<unsigned int,double>(28,Obs28));
if(DEBUG) std::cout<<"------ LR observable 28 "<<Obs28<<" calculated ------"<<std::endl;
double Obs29 = FW_momentum_5;
evtselectVarVal.push_back(std::pair<unsigned int,double>(29,Obs29));
if(DEBUG) std::cout<<"------ LR observable 29 "<<Obs29<<" calculated ------"<<std::endl;
double Obs30 = FW_momentum_6;
evtselectVarVal.push_back(std::pair<unsigned int,double>(30,Obs30));
if(DEBUG) std::cout<<"------ LR observable 30 "<<Obs30<<" calculated ------"<<std::endl;
// Thrust, thrust major and thrust minor
TVector3 n(0,0,0), n_tmp(0,0,0);
double Thrust_D = 0, Thrust_N = 0;
double Thrust = -1, Thrust_tmp =0;
double Thrust_Major = -1, Thrust_Major_tmp =0;
double Thrust_Minor = -1;
for(unsigned int i=0;i<TopJets.size();i++)
{
Thrust_D += TopJets[i].p();
}
Thrust_D += Lept->P();
if(Thrust_D>0)
{
// Calculation of the thrust :
for(unsigned int i=0;i<720;i++)
{
for(unsigned int j=0;j<360;j++)
{
n_tmp.SetX(cos((2*PI/720)*i)*sin((PI/360)*j));
n_tmp.SetY(sin((2*PI/720)*i)*sin((PI/360)*j));
n_tmp.SetZ(cos((PI/360)*j));
for(unsigned int k=0;k<TopJets.size();k++)
{
Thrust_N += fabs(TopJets[k].px()*(n_tmp.x())+TopJets[k].py()*(n_tmp.y())+TopJets[k].pz()*(n_tmp.z()));
}
Thrust_N += fabs(Lept->Px()*(n_tmp.x())+Lept->Py()*(n_tmp.y())+Lept->Pz()*(n_tmp.z()));
Thrust_tmp = Thrust_N/Thrust_D;
Thrust_N = 0;
if(Thrust_tmp > Thrust)
{
Thrust = Thrust_tmp;
n.SetXYZ(n_tmp.x(),n_tmp.y(),n_tmp.z());
}
}
}
// Calculation of the thrust major :
TVector3 nT = n.Orthogonal();
nT = nT.Unit();
for(unsigned int i=0;i<720;i++)
{
nT.Rotate((2*PI/720)*i,n);
for(unsigned int j=0;j<TopJets.size();j++)
{
Thrust_N += fabs(TopJets[j].px()*(nT.x())+TopJets[j].py()*(nT.y())+TopJets[j].pz()*(nT.z()));
}
Thrust_N += fabs(Lept->Px()*nT.x()+Lept->Py()*(nT.y())+Lept->Pz()*(nT.z()));
Thrust_Major_tmp = Thrust_N/Thrust_D;
Thrust_N = 0;
if(Thrust_Major_tmp > Thrust_Major)
{
Thrust_Major = Thrust_Major_tmp;
}
}
// Calculation of the thrust minor :
TVector3 nMinor = nT.Cross(n);
nMinor = nMinor.Unit();
for(unsigned int i=0;i<TopJets.size();i++)
{
Thrust_N += fabs(TopJets[i].px()*(nMinor.x())+TopJets[i].py()*(nMinor.y())+TopJets[i].pz()*(nMinor.z()));
}
Thrust_N += fabs(Lept->Px()*nMinor.x()+Lept->Py()*(nMinor.y())+Lept->Pz()*(nMinor.z()));
Thrust_Minor = Thrust_N/Thrust_D;
}
double Obs31 = Thrust;
evtselectVarVal.push_back(std::pair<unsigned int,double>(31,Obs31));
if(DEBUG) std::cout<<"------ LR observable 31 "<<Obs31<<" calculated ------"<<std::endl;
double Obs32 = Thrust_Major;
evtselectVarVal.push_back(std::pair<unsigned int,double>(32,Obs32));
if(DEBUG) std::cout<<"------ LR observable 32 "<<Obs32<<" calculated ------"<<std::endl;
double Obs33 = Thrust_Minor;
evtselectVarVal.push_back(std::pair<unsigned int,double>(33,Obs33));
if(DEBUG) std::cout<<"------ LR observable 33 "<<Obs33<<" calculated ------"<<std::endl;
// Oblateness
double Obs34 = Thrust_Major-Thrust_Minor;
evtselectVarVal.push_back(std::pair<unsigned int,double>(34,Obs34));
if(DEBUG) std::cout<<"------ LR observable 34 "<<Obs34<<" calculated ------"<<std::endl;
// Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino
TMatrixDSym Matrix_NoNu(3);
double PX2_NoNu = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2)+std::pow(Lept->Px(),2);
double PY2_NoNu = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2)+std::pow(Lept->Py(),2);
double PZ2_NoNu = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2)+std::pow(Lept->Pz(),2);
double P2_NoNu = PX2_NoNu+PY2_NoNu+PZ2_NoNu;
double PXY_NoNu = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py()+Lept->Px()*Lept->Py();
double PXZ_NoNu = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz()+Lept->Px()*Lept->Pz();
double PYZ_NoNu = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz()+Lept->Py()*Lept->Pz();
Matrix_NoNu(0,0) = PX2_NoNu/P2_NoNu; Matrix_NoNu(0,1) = PXY_NoNu/P2_NoNu; Matrix_NoNu(0,2) = PXZ_NoNu/P2_NoNu;
Matrix_NoNu(1,0) = PXY_NoNu/P2_NoNu; Matrix_NoNu(1,1) = PY2_NoNu/P2_NoNu; Matrix_NoNu(1,2) = PYZ_NoNu/P2_NoNu;
Matrix_NoNu(2,0) = PXZ_NoNu/P2_NoNu; Matrix_NoNu(2,1) = PYZ_NoNu/P2_NoNu; Matrix_NoNu(2,2) = PZ2_NoNu/P2_NoNu;
TMatrixDSymEigen pTensor_NoNu(Matrix_NoNu);
std::vector<double> EigValues_NoNu;
EigValues_NoNu.clear();
for(int i=0;i<3;i++)
{
EigValues_NoNu.push_back(pTensor_NoNu.GetEigenValues()[i]);
}
std::sort(EigValues_NoNu.begin(),EigValues_NoNu.end(),dComparator);
double Sphericity_NoNu = 1.5*(EigValues_NoNu[1]+EigValues_NoNu[2]);
double Aplanarity_NoNu = 1.5*EigValues_NoNu[2];
double Obs35 = (isnan(Sphericity_NoNu) ? 0 : Sphericity_NoNu);
evtselectVarVal.push_back(std::pair<unsigned int,double>(35,Obs35));
if(DEBUG) std::cout<<"------ LR observable 35 "<<Obs35<<" calculated ------"<<std::endl;
double Obs36 = (isnan(Aplanarity_NoNu) ? 0 : Aplanarity_NoNu);
evtselectVarVal.push_back(std::pair<unsigned int,double>(36,Obs36));
if(DEBUG) std::cout<<"------ LR observable 36 "<<Obs36<<" calculated ------"<<std::endl;
// Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino or lepton
TMatrixDSym Matrix_NoNuNoLep(3);
double PX2_NoNuNoLep = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+std::pow(Lepb->Px(),2);
double PY2_NoNuNoLep = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+std::pow(Lepb->Py(),2);
double PZ2_NoNuNoLep = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+std::pow(Lepb->Pz(),2);
double P2_NoNuNoLep = PX2_NoNuNoLep+PY2_NoNuNoLep+PZ2_NoNuNoLep;
double PXY_NoNuNoLep = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+Lepb->Px()*Lepb->Py();
double PXZ_NoNuNoLep = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+Lepb->Px()*Lepb->Pz();
double PYZ_NoNuNoLep = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+Lepb->Py()*Lepb->Pz();
Matrix_NoNuNoLep(0,0) = PX2_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(0,1) = PXY_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(0,2) = PXZ_NoNuNoLep/P2_NoNuNoLep;
Matrix_NoNuNoLep(1,0) = PXY_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(1,1) = PY2_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(1,2) = PYZ_NoNuNoLep/P2_NoNuNoLep;
Matrix_NoNuNoLep(2,0) = PXZ_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(2,1) = PYZ_NoNuNoLep/P2_NoNuNoLep; Matrix_NoNuNoLep(2,2) = PZ2_NoNuNoLep/P2_NoNuNoLep;
TMatrixDSymEigen pTensor_NoNuNoLep(Matrix_NoNuNoLep);
std::vector<double> EigValues_NoNuNoLep;
EigValues_NoNuNoLep.clear();
for(int i=0;i<3;i++)
{
EigValues_NoNuNoLep.push_back(pTensor_NoNuNoLep.GetEigenValues()[i]);
}
std::sort(EigValues_NoNuNoLep.begin(),EigValues_NoNuNoLep.end(),dComparator);
double Sphericity_NoNuNoLep = 1.5*(EigValues_NoNuNoLep[1]+EigValues_NoNuNoLep[2]);
double Aplanarity_NoNuNoLep = 1.5*EigValues_NoNuNoLep[2];
double Obs37 = (isnan(Sphericity_NoNuNoLep) ? 0 : Sphericity_NoNuNoLep);
evtselectVarVal.push_back(std::pair<unsigned int,double>(37,Obs37));
if(DEBUG) std::cout<<"------ LR observable 37 "<<Obs37<<" calculated ------"<<std::endl;
double Obs38 = (isnan(Aplanarity_NoNuNoLep) ? 0 : Aplanarity_NoNuNoLep);
evtselectVarVal.push_back(std::pair<unsigned int,double>(38,Obs38));
if(DEBUG) std::cout<<"------ LR observable 38 "<<Obs38<<" calculated ------"<<std::endl;
// Put the vector in the TtSemiEvtSolution
TS.setLRSignalEvtObservables(evtselectVarVal);
if(DEBUG) std::cout<<"------ Observable values stored in the event ------"<<std::endl;
delete Hadp;
if(DEBUG) std::cout<<"------ Pointer to Hadp deleted ------"<<std::endl;
delete Hadq;
if(DEBUG) std::cout<<"------ Pointer to Hadq deleted ------"<<std::endl;
delete Hadb;
if(DEBUG) std::cout<<"------ Pointer to Hadb deleted ------"<<std::endl;
delete Lepb;
if(DEBUG) std::cout<<"------ Pointer to Lepb deleted ------"<<std::endl;
delete Lept;
if(DEBUG) std::cout<<"------ Pointer to Lepn deleted ------"<<std::endl;
delete Lepn;
if(DEBUG) std::cout<<"------ Pointer to Mez deleted ------"<<std::endl;
delete Mez;
if(DEBUG) std::cout<<"------------ LR observables calculated -----------"<<std::endl;
}
Definition at line 60 of file TtSemiLRSignalSelObservables.h.
Referenced by operator()().
Definition at line 70 of file TtSemiLRSignalSelObservables.h.
Referenced by operator()().
Definition at line 50 of file TtSemiLRSignalSelObservables.h.
Referenced by operator()().
std::vector<std::pair<unsigned int,double> > TtSemiLRSignalSelObservables::evtselectVarVal [private] |
Definition at line 40 of file TtSemiLRSignalSelObservables.h.
Referenced by operator()().
1.7.3