SFrameTools/src/Utils.cxx

#include "../include/Utils.h"


// global function to define a tagged jet

bool IsTagged(Jet & jet, E_BtagType type)
{
    if(type==e_CSVL && jet.btag_combinedSecondaryVertex()>0.244) return true;
    if(type==e_CSVM && jet.btag_combinedSecondaryVertex()>0.679) return true;
    if(type==e_CSVT && jet.btag_combinedSecondaryVertex()>0.898) return true;
    if(type==e_JPL && jet.btag_jetProbability()>0.275) return true;
    if(type==e_JPM && jet.btag_jetProbability()>0.545) return true;
    if(type==e_JPT && jet.btag_jetProbability()>0.790) return true;       
    
    return false;
}

//variable HEP Tagger from Rebekka

bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3)
{
    double mjet;
    double ptjet;
    int nsubjets;

    double topmass=172.3;
    double wmass=80.4;

    nsubjets=topjet.numberOfDaughters();

    LorentzVector allsubjets(0,0,0,0);

    for(int j=0; j<topjet.numberOfDaughters(); ++j) {
        allsubjets += topjet.subjets()[j].v4();
    }
    if(!allsubjets.isTimelike()) {
        mjet=0;
        return false;
    }

    mjet = allsubjets.M();
    ptjet= allsubjets.Pt();

    double m12, m13, m23;

    //The subjetcs have to be three
    if(nsubjets==3) {

        std::vector<Particle> subjets = topjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        m12 = 0;
        if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
            m12=(subjets[0].v4()+subjets[1].v4()).M();
        m13 = 0;
        if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
            m13=(subjets[0].v4()+subjets[2].v4()).M();
        m23 = 0;
        if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
            m23 = (subjets[1].v4()+subjets[2].v4()).M();

    } else {
        return false;
    }

    double rmin=massWindowLower*wmass/topmass;
    double rmax=massWindowUpper*wmass/topmass;

    int keep=0;

    //Conditions on the subjects: at least one has to be true
    //1 condition
    if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;

    //2 condition
    double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
    double cond2cent=1-pow(m23/mjet,2);
    double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));

    if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;

    //3 condition
    double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
    double cond3cent=1-pow(m23/mjet,2);
    double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));

    if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;

    //Final requirement: at least one of the three subjets conditions and total pt
    if(keep==1 && ptjet>ptJetMin) {
        return true;
    } else {
        return false;
    }

}


//HEP Tagger from Ivan

bool HepTopTag(TopJet topjet)
{

    double mjet;
    double ptjet;
    int nsubjets;

    double topmass=172.3;
    double wmass=80.4;

    nsubjets=topjet.numberOfDaughters();

    LorentzVector allsubjets(0,0,0,0);

    for(int j=0; j<topjet.numberOfDaughters(); ++j) {
        allsubjets += topjet.subjets()[j].v4();
    }
    if(!allsubjets.isTimelike()) {
        mjet=0;
        return false;
    }

    mjet = allsubjets.M();
    ptjet= allsubjets.Pt();

    double m12, m13, m23;

    //The subjetcs have to be three
    if(nsubjets==3) {

        std::vector<Particle> subjets = topjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        m12 = 0;
        if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
            m12=(subjets[0].v4()+subjets[1].v4()).M();
        m13 = 0;
        if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
            m13=(subjets[0].v4()+subjets[2].v4()).M();
        m23 = 0;
        if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
            m23 = (subjets[1].v4()+subjets[2].v4()).M();

    } else {
        return false;
    }

    double rmin=0.85*wmass/topmass;
    double rmax=1.15*wmass/topmass;

    int keep=0;

    //Conditions on the subjects: at least one has to be true
    //1 condition
    if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;

    //2 condition
    double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
    double cond2cent=1-pow(m23/mjet,2);
    double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));

    if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;

    //3 condition
    double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
    double cond3cent=1-pow(m23/mjet,2);
    double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));

    if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;

    //Final requirement: at least one of the three subjets conditions and total pt
    if(keep==1 && ptjet>200) {
        return true;
    } else {
        return false;
    }

}


//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper)
{

    nsubjets=topjet.numberOfDaughters();

    LorentzVector allsubjets(0,0,0,0);

    for(int j=0; j<topjet.numberOfDaughters(); ++j) {
        allsubjets += topjet.subjets()[j].v4();
    }
    if(!allsubjets.isTimelike()) {
        mjet=0;
        mmin=0;
//  mminLower=50;
//   mjetLower=140;
//   mjetUpper=250;
        return false;
    }

    mjet = allsubjets.M();

    if(nsubjets>=3) {

        std::vector<Particle> subjets = topjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        double m01 = 0;
        if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
            m01=(subjets[0].v4()+subjets[1].v4()).M();
        double m02 = 0;
        if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
            m02=(subjets[0].v4()+subjets[2].v4()).M();
        double m12 = 0;
        if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
            m12 = (subjets[1].v4()+subjets[2].v4()).M();

        //minimum pairwise mass
        mmin = std::min(m01,std::min(m02,m12));
    }

    //at least 3 sub-jets
    if(nsubjets<3) return false;
    //minimum pairwise mass > 50 GeV/c^2
    if(mmin<mminLower) return false;
    //jet mass between 140 and 250 GeV/c^2
    if(mjet<mjetLower || mjet>mjetUpper) return false;

    return true;
}


bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin)
{

    nsubjets=topjet.numberOfDaughters();

    LorentzVector allsubjets(0,0,0,0);

    for(int j=0; j<topjet.numberOfDaughters(); ++j) {
        allsubjets += topjet.subjets()[j].v4();
    }
    if(!allsubjets.isTimelike()) {
        mjet=0;
        mmin=0;
        return false;
    }

    mjet = allsubjets.M();

    if(nsubjets>=3) {

        std::vector<Particle> subjets = topjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        double m01 = 0;
        if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
            m01=(subjets[0].v4()+subjets[1].v4()).M();
        double m02 = 0;
        if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
            m02=(subjets[0].v4()+subjets[2].v4()).M();
        double m12 = 0;
        if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
            m12 = (subjets[1].v4()+subjets[2].v4()).M();

        //minimum pairwise mass
        mmin = std::min(m01,std::min(m02,m12));
    }

    //at least 3 sub-jets
    if(nsubjets<3) return false;
    //minimum pairwise mass > 50 GeV/c^2
    if(mmin<50) return false;
    //jet mass between 140 and 250 GeV/c^2
    if(mjet<140 || mjet>250) return false;

    return true;
}

Jet* nextJet(const Particle *p, std::vector<Jet> *jets)
{

    double deltarmin = double_infinity();
    Jet* nextjet=0;
    for(unsigned int i=0; i<jets->size(); ++i) {
        if(jets->at(i).deltaR(*p) < deltarmin) {
            deltarmin = jets->at(i).deltaR(*p);
            nextjet = &jets->at(i);
        }
    }

    return nextjet;
}

bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop)
{

    nsubjets=prunedjet.numberOfDaughters();

    mjet = 0;
    if(prunedjet.v4().isTimelike())
        mjet = prunedjet.v4().M();

    //calculate mass drop for first sub-jet ordered in pt
    massdrop = 0;
    if(nsubjets>=1 && mjet>0) {

        std::vector< Particle > subjets = prunedjet.subjets();
        sort(subjets.begin(), subjets.end(), HigherPt());

        double m1 = 0;
        if(subjets[0].v4().isTimelike())
            m1 = subjets[0].v4().M();

        massdrop = m1/mjet;
    }

    //at least 2 sub-jets
    if(nsubjets<2) return false;
    //60 GeV < pruned jet mass < 100 GeV
    if(mjet <= 60 || mjet >= 100) return false;
    //mass drop < 0.4
    if(massdrop>=0.4) return false;

    return true;

}

double pTrel(const Particle *p, std::vector<Jet> *jets)
{

    double ptrel=0;

    Jet* nextjet =  nextJet(p,jets);

    TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
    TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());

    if(p3.Mag()!=0 && jet3.Mag()!=0) {
        double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
        ptrel = p3.Mag()*sin_alpha;
    } else {
        std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
    }

    return ptrel;
}

double deltaRmin(const Particle *p, std::vector<Jet> *jets)
{
    return nextJet(p,jets)->deltaR(*p);
}

TVector3 toVector(LorentzVector v4)
{

    TVector3 v3(0,0,0);
    v3.SetX(v4.X());
    v3.SetY(v4.Y());
    v3.SetZ(v4.Z());
    return v3;
}

TVector3 toVector(LorentzVectorXYZE v4)
{

    TVector3 v3(0,0,0);
    v3.SetX(v4.X());
    v3.SetY(v4.Y());
    v3.SetZ(v4.Z());
    return v3;
}

LorentzVectorXYZE toXYZ(LorentzVector v4)
{

    LorentzVectorXYZE v4_new(0,0,0,0);
    v4_new.SetPx(v4.Px());
    v4_new.SetPy(v4.Py());
    v4_new.SetPz(v4.Pz());
    v4_new.SetE(v4.E());
    return v4_new;
}

LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
{

    LorentzVector v4_new(0,0,0,0);
    v4_new.SetPt(v4.Pt());
    v4_new.SetEta(v4.Eta());
    v4_new.SetPhi(v4.Phi());
    v4_new.SetE(v4.E());
    return v4_new;
}

double deltaR(LorentzVector v1, LorentzVector v2)
{

    Particle p1;
    p1.set_v4(v1);
    Particle p2;
    p2.set_v4(v2);
    return p1.deltaR(p2);
}

double double_infinity()
{
    return std::numeric_limits<double>::infinity() ;
}

int int_infinity()
{
    return std::numeric_limits<int>::max() ;
}

int myPow(int x, unsigned int p)
{
    int i = 1;
    for (unsigned int j = 1; j <= p; j++)  i *= x;
    return i;
}

int JetFlavor(Jet *jet)
{

    EventCalc* calc = EventCalc::Instance();

    std::vector< GenParticle >* genparticles = calc->GetGenParticles();
    if(genparticles) {

        //fill pdg IDs of all matched GenParticles
        std::vector<int> matched_genparticle_ids;
        for(unsigned int i=0; i<genparticles->size(); ++i) {

            GenParticle genp = genparticles->at(i);

            //only take status 3 particles into account
            if(genp.status()!=3) continue;

            if(jet->deltaR(genp)<0.5)
                matched_genparticle_ids.push_back(genp.pdgId());
        }

        //search for b quarks first
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
        }
        //no b quark -> search for c quarks
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
        }
        for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
            if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
        }

    }

    //no matched GenParticle -> return default value 0
    return 0;

}
Revision:	1.8
Committed:	Wed Dec 19 00:02:35 2012 UTC (12 years, 4 months ago) by bazterra
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Jan-17-2013-v2, Jan-17-2013-v1, Jan-16-2012-v1, Jan-09-2012-v2, Jan-09-2012-v1, Dec-26-2012-v1, Dec-20-2012-v1
Changes since 1.7:	+373 -348 lines
Log Message:	Adding common function to check if a jet is btagged.
#	User	Rev	Content
1	peiffer	1.1	#include "../include/Utils.h"
2
3
4	bazterra	1.8	// global function to define a tagged jet
5
6			bool IsTagged(Jet & jet, E_BtagType type)
7			{
8			if(type==e_CSVL && jet.btag_combinedSecondaryVertex()>0.244) return true;
9			if(type==e_CSVM && jet.btag_combinedSecondaryVertex()>0.679) return true;
10			if(type==e_CSVT && jet.btag_combinedSecondaryVertex()>0.898) return true;
11			if(type==e_JPL && jet.btag_jetProbability()>0.275) return true;
12			if(type==e_JPM && jet.btag_jetProbability()>0.545) return true;
13			if(type==e_JPT && jet.btag_jetProbability()>0.790) return true;
14
15			return false;
16			}
17
18	hoeing	1.6	//variable HEP Tagger from Rebekka
19
20	bazterra	1.8	bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3)
21			{
22			double mjet;
23			double ptjet;
24			int nsubjets;
25
26			double topmass=172.3;
27			double wmass=80.4;
28
29			nsubjets=topjet.numberOfDaughters();
30	hoeing	1.6
31	bazterra	1.8	LorentzVector allsubjets(0,0,0,0);
32	hoeing	1.6
33	bazterra	1.8	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
34			allsubjets += topjet.subjets()[j].v4();
35			}
36			if(!allsubjets.isTimelike()) {
37			mjet=0;
38			return false;
39			}
40	hoeing	1.6
41	bazterra	1.8	mjet = allsubjets.M();
42			ptjet= allsubjets.Pt();
43	hoeing	1.6
44	bazterra	1.8	double m12, m13, m23;
45	hoeing	1.6
46	bazterra	1.8	//The subjetcs have to be three
47			if(nsubjets==3) {
48	hoeing	1.6
49	bazterra	1.8	std::vector<Particle> subjets = topjet.subjets();
50			sort(subjets.begin(), subjets.end(), HigherPt());
51
52			m12 = 0;
53			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
54			m12=(subjets[0].v4()+subjets[1].v4()).M();
55			m13 = 0;
56			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
57			m13=(subjets[0].v4()+subjets[2].v4()).M();
58			m23 = 0;
59			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
60			m23 = (subjets[1].v4()+subjets[2].v4()).M();
61	hoeing	1.6
62	bazterra	1.8	} else {
63			return false;
64			}
65	hoeing	1.6
66	bazterra	1.8	double rmin=massWindowLower*wmass/topmass;
67			double rmax=massWindowUpper*wmass/topmass;
68	hoeing	1.6
69	bazterra	1.8	int keep=0;
70	hoeing	1.6
71	bazterra	1.8	//Conditions on the subjects: at least one has to be true
72			//1 condition
73			if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
74
75			//2 condition
76			double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
77			double cond2cent=1-pow(m23/mjet,2);
78			double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
79
80			if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;
81
82			//3 condition
83			double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
84			double cond3cent=1-pow(m23/mjet,2);
85			double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
86
87			if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;
88
89			//Final requirement: at least one of the three subjets conditions and total pt
90			if(keep==1 && ptjet>ptJetMin) {
91			return true;
92			} else {
93			return false;
94			}
95	hoeing	1.6
96			}
97
98
99
100
101
102
103			//HEP Tagger from Ivan
104
105	bazterra	1.8	bool HepTopTag(TopJet topjet)
106			{
107
108			double mjet;
109			double ptjet;
110			int nsubjets;
111
112			double topmass=172.3;
113			double wmass=80.4;
114
115			nsubjets=topjet.numberOfDaughters();
116	hoeing	1.6
117	bazterra	1.8	LorentzVector allsubjets(0,0,0,0);
118	hoeing	1.6
119	bazterra	1.8	for(int j=0; j<topjet.numberOfDaughters(); ++j) {
120			allsubjets += topjet.subjets()[j].v4();
121			}
122			if(!allsubjets.isTimelike()) {
123			mjet=0;
124			return false;
125			}
126	hoeing	1.6
127	bazterra	1.8	mjet = allsubjets.M();
128			ptjet= allsubjets.Pt();
129	hoeing	1.6
130	bazterra	1.8	double m12, m13, m23;
131	hoeing	1.6
132	bazterra	1.8	//The subjetcs have to be three
133			if(nsubjets==3) {
134	hoeing	1.6
135	bazterra	1.8	std::vector<Particle> subjets = topjet.subjets();
136			sort(subjets.begin(), subjets.end(), HigherPt());
137
138			m12 = 0;
139			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
140			m12=(subjets[0].v4()+subjets[1].v4()).M();
141			m13 = 0;
142			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
143			m13=(subjets[0].v4()+subjets[2].v4()).M();
144			m23 = 0;
145			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
146			m23 = (subjets[1].v4()+subjets[2].v4()).M();
147	hoeing	1.6
148	bazterra	1.8	} else {
149			return false;
150			}
151	hoeing	1.6
152	bazterra	1.8	double rmin=0.85*wmass/topmass;
153			double rmax=1.15*wmass/topmass;
154	hoeing	1.6
155	bazterra	1.8	int keep=0;
156	hoeing	1.6
157	bazterra	1.8	//Conditions on the subjects: at least one has to be true
158			//1 condition
159			if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
160
161			//2 condition
162			double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
163			double cond2cent=1-pow(m23/mjet,2);
164			double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
165
166			if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;
167
168			//3 condition
169			double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
170			double cond3cent=1-pow(m23/mjet,2);
171			double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
172
173			if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;
174
175			//Final requirement: at least one of the three subjets conditions and total pt
176			if(keep==1 && ptjet>200) {
177			return true;
178			} else {
179			return false;
180			}
181	hoeing	1.6
182			}
183
184
185			//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
186	bazterra	1.8	bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper)
187			{
188
189			nsubjets=topjet.numberOfDaughters();
190
191			LorentzVector allsubjets(0,0,0,0);
192
193			for(int j=0; j<topjet.numberOfDaughters(); ++j) {
194			allsubjets += topjet.subjets()[j].v4();
195			}
196			if(!allsubjets.isTimelike()) {
197			mjet=0;
198			mmin=0;
199			// mminLower=50;
200	hoeing	1.6	// mjetLower=140;
201			// mjetUpper=250;
202	bazterra	1.8	return false;
203			}
204
205			mjet = allsubjets.M();
206
207			if(nsubjets>=3) {
208
209			std::vector<Particle> subjets = topjet.subjets();
210			sort(subjets.begin(), subjets.end(), HigherPt());
211
212			double m01 = 0;
213			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
214			m01=(subjets[0].v4()+subjets[1].v4()).M();
215			double m02 = 0;
216			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
217			m02=(subjets[0].v4()+subjets[2].v4()).M();
218			double m12 = 0;
219			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
220			m12 = (subjets[1].v4()+subjets[2].v4()).M();
221
222			//minimum pairwise mass
223			mmin = std::min(m01,std::min(m02,m12));
224			}
225
226			//at least 3 sub-jets
227			if(nsubjets<3) return false;
228			//minimum pairwise mass > 50 GeV/c^2
229			if(mmin<mminLower) return false;
230			//jet mass between 140 and 250 GeV/c^2
231			if(mjet<mjetLower \|\| mjet>mjetUpper) return false;
232
233			return true;
234			}
235
236
237			bool TopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin)
238			{
239
240			nsubjets=topjet.numberOfDaughters();
241
242			LorentzVector allsubjets(0,0,0,0);
243
244			for(int j=0; j<topjet.numberOfDaughters(); ++j) {
245			allsubjets += topjet.subjets()[j].v4();
246			}
247			if(!allsubjets.isTimelike()) {
248			mjet=0;
249			mmin=0;
250			return false;
251			}
252
253			mjet = allsubjets.M();
254
255			if(nsubjets>=3) {
256
257			std::vector<Particle> subjets = topjet.subjets();
258			sort(subjets.begin(), subjets.end(), HigherPt());
259
260			double m01 = 0;
261			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
262			m01=(subjets[0].v4()+subjets[1].v4()).M();
263			double m02 = 0;
264			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
265			m02=(subjets[0].v4()+subjets[2].v4()).M();
266			double m12 = 0;
267			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
268			m12 = (subjets[1].v4()+subjets[2].v4()).M();
269
270			//minimum pairwise mass
271			mmin = std::min(m01,std::min(m02,m12));
272			}
273	peiffer	1.1
274	bazterra	1.8	//at least 3 sub-jets
275			if(nsubjets<3) return false;
276			//minimum pairwise mass > 50 GeV/c^2
277			if(mmin<50) return false;
278			//jet mass between 140 and 250 GeV/c^2
279			if(mjet<140 \|\| mjet>250) return false;
280	peiffer	1.1
281	bazterra	1.8	return true;
282	peiffer	1.1	}
283	peiffer	1.2
284	bazterra	1.8	Jet* nextJet(const Particle p, std::vector<Jet> jets)
285			{
286	peiffer	1.2
287	bazterra	1.8	double deltarmin = double_infinity();
288			Jet* nextjet=0;
289			for(unsigned int i=0; i<jets->size(); ++i) {
290			if(jets->at(i).deltaR(*p) < deltarmin) {
291			deltarmin = jets->at(i).deltaR(*p);
292			nextjet = &jets->at(i);
293			}
294	peiffer	1.2	}
295
296	bazterra	1.8	return nextjet;
297	peiffer	1.2	}
298
299	bazterra	1.8	bool WTag(TopJet prunedjet, double& mjet, int &nsubjets, double& massdrop)
300			{
301	peiffer	1.5
302	bazterra	1.8	nsubjets=prunedjet.numberOfDaughters();
303	peiffer	1.5
304	bazterra	1.8	mjet = 0;
305			if(prunedjet.v4().isTimelike())
306			mjet = prunedjet.v4().M();
307	peiffer	1.5
308	bazterra	1.8	//calculate mass drop for first sub-jet ordered in pt
309			massdrop = 0;
310			if(nsubjets>=1 && mjet>0) {
311	peiffer	1.5
312	bazterra	1.8	std::vector< Particle > subjets = prunedjet.subjets();
313			sort(subjets.begin(), subjets.end(), HigherPt());
314	peiffer	1.5
315	bazterra	1.8	double m1 = 0;
316			if(subjets[0].v4().isTimelike())
317			m1 = subjets[0].v4().M();
318	peiffer	1.5
319	bazterra	1.8	massdrop = m1/mjet;
320			}
321	peiffer	1.5
322	bazterra	1.8	//at least 2 sub-jets
323			if(nsubjets<2) return false;
324			//60 GeV < pruned jet mass < 100 GeV
325			if(mjet <= 60 \|\| mjet >= 100) return false;
326			//mass drop < 0.4
327			if(massdrop>=0.4) return false;
328
329			return true;
330	peiffer	1.5
331			}
332
333	bazterra	1.8	double pTrel(const Particle p, std::vector<Jet> jets)
334			{
335	peiffer	1.2
336	bazterra	1.8	double ptrel=0;
337	peiffer	1.2
338	bazterra	1.8	Jet* nextjet = nextJet(p,jets);
339	peiffer	1.2
340	bazterra	1.8	TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
341			TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
342	peiffer	1.2
343	bazterra	1.8	if(p3.Mag()!=0 && jet3.Mag()!=0) {
344			double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
345			ptrel = p3.Mag()*sin_alpha;
346			} else {
347			std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
348			}
349	peiffer	1.2
350	bazterra	1.8	return ptrel;
351	peiffer	1.2	}
352
353	bazterra	1.8	double deltaRmin(const Particle p, std::vector<Jet> jets)
354			{
355			return nextJet(p,jets)->deltaR(*p);
356	peiffer	1.2	}
357
358	bazterra	1.8	TVector3 toVector(LorentzVector v4)
359			{
360	peiffer	1.4
361	bazterra	1.8	TVector3 v3(0,0,0);
362			v3.SetX(v4.X());
363			v3.SetY(v4.Y());
364			v3.SetZ(v4.Z());
365			return v3;
366	peiffer	1.4	}
367
368	bazterra	1.8	TVector3 toVector(LorentzVectorXYZE v4)
369			{
370	peiffer	1.4
371	bazterra	1.8	TVector3 v3(0,0,0);
372			v3.SetX(v4.X());
373			v3.SetY(v4.Y());
374			v3.SetZ(v4.Z());
375			return v3;
376	peiffer	1.4	}
377
378	bazterra	1.8	LorentzVectorXYZE toXYZ(LorentzVector v4)
379			{
380	peiffer	1.4
381	bazterra	1.8	LorentzVectorXYZE v4_new(0,0,0,0);
382			v4_new.SetPx(v4.Px());
383			v4_new.SetPy(v4.Py());
384			v4_new.SetPz(v4.Pz());
385			v4_new.SetE(v4.E());
386			return v4_new;
387	peiffer	1.4	}
388
389	bazterra	1.8	LorentzVector toPtEtaPhi(LorentzVectorXYZE v4)
390			{
391	peiffer	1.4
392	bazterra	1.8	LorentzVector v4_new(0,0,0,0);
393			v4_new.SetPt(v4.Pt());
394			v4_new.SetEta(v4.Eta());
395			v4_new.SetPhi(v4.Phi());
396			v4_new.SetE(v4.E());
397			return v4_new;
398	peiffer	1.4	}
399
400	bazterra	1.8	double deltaR(LorentzVector v1, LorentzVector v2)
401			{
402	peiffer	1.4
403	bazterra	1.8	Particle p1;
404			p1.set_v4(v1);
405			Particle p2;
406			p2.set_v4(v2);
407			return p1.deltaR(p2);
408	peiffer	1.4	}
409	peiffer	1.2
410	bazterra	1.8	double double_infinity()
411			{
412			return std::numeric_limits<double>::infinity() ;
413	peiffer	1.2	}
414
415	bazterra	1.8	int int_infinity()
416			{
417			return std::numeric_limits<int>::max() ;
418	peiffer	1.2	}
419	peiffer	1.4
420	bazterra	1.8	int myPow(int x, unsigned int p)
421			{
422			int i = 1;
423			for (unsigned int j = 1; j <= p; j++) i *= x;
424			return i;
425	peiffer	1.4	}
426	peiffer	1.7
427	bazterra	1.8	int JetFlavor(Jet *jet)
428			{
429
430			EventCalc* calc = EventCalc::Instance();
431	peiffer	1.7
432	bazterra	1.8	std::vector< GenParticle >* genparticles = calc->GetGenParticles();
433			if(genparticles) {
434	peiffer	1.7
435	bazterra	1.8	//fill pdg IDs of all matched GenParticles
436			std::vector<int> matched_genparticle_ids;
437			for(unsigned int i=0; i<genparticles->size(); ++i) {
438	peiffer	1.7
439	bazterra	1.8	GenParticle genp = genparticles->at(i);
440
441			//only take status 3 particles into account
442			if(genp.status()!=3) continue;
443	peiffer	1.7
444	bazterra	1.8	if(jet->deltaR(genp)<0.5)
445			matched_genparticle_ids.push_back(genp.pdgId());
446			}
447	peiffer	1.7
448	bazterra	1.8	//search for b quarks first
449			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
450			if(abs(matched_genparticle_ids[i])==5) return matched_genparticle_ids[i];
451			}
452			//no b quark -> search for c quarks
453			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
454			if(abs(matched_genparticle_ids[i])==4) return matched_genparticle_ids[i];
455			}
456			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
457			if(abs(matched_genparticle_ids[i])==3) return matched_genparticle_ids[i];
458			}
459			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
460			if(abs(matched_genparticle_ids[i])==2) return matched_genparticle_ids[i];
461			}
462			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
463			if(abs(matched_genparticle_ids[i])==1) return matched_genparticle_ids[i];
464			}
465			for(unsigned int i=0; i<matched_genparticle_ids.size(); ++i) {
466			if(abs(matched_genparticle_ids[i])==21) return matched_genparticle_ids[i];
467			}
468	peiffer	1.7
469			}
470
471	bazterra	1.8	//no matched GenParticle -> return default value 0
472			return 0;
473	peiffer	1.7
474			}