SFrameTools/src/Utils.cxx

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

#include <fastjet/JetDefinition.hh>
#include <fastjet/PseudoJet.hh>
#include <fastjet/ClusterSequence.hh>
#include <fastjet/ClusterSequenceArea.hh>
#include <fastjet/GhostedAreaSpec.hh>
namespace external {
#include "../include/HEPTopTagger.h"
}


//subjet b-tagger, returns number of b-tagged subjets

int subJetBTag(TopJet topjet, E_BtagType type){
  int nBTagsSub = 0;
  float discriminator_cut;
  if(type==e_CSVL) discriminator_cut = 0.244;
  if(type==e_CSVM) discriminator_cut = 0.679;
  if(type==e_CSVT) discriminator_cut = 0.898;

  //Create a vector of subjets
  std::vector<Particle> subjets_top;
  //Create a float vector of the subjets discriminators
  std::vector<float> btagsub_combinedSecondaryVertex_top;

  //Fill the vector of subjets with the subjets of a topjet
  subjets_top=topjet.subjets();
  //Fill the vector of discriminators with the discriminators of the subjets of a certain topjet
  btagsub_combinedSecondaryVertex_top=topjet.btagsub_combinedSecondaryVertex();

  //Looping over subjets and checking if they are b-tagged
  for(unsigned int i=0; i < btagsub_combinedSecondaryVertex_top.size(); ++i){
    float test=btagsub_combinedSecondaryVertex_top[i];
    if(test>discriminator_cut){
      nBTagsSub += 1;
      //This means it is b-tagged
    }
  }
  return nBTagsSub;
}


bool HiggsTag(TopJet topjet, E_BtagType type1, E_BtagType type2){
  int nBTagsSub1 = 0;
  int nBTagsSub2 = 0;
  float discriminator_cut1;
  float discriminator_cut2;
  if(type1==e_CSVL) discriminator_cut1 = 0.244;
  if(type1==e_CSVM) discriminator_cut1 = 0.679;
  if(type1==e_CSVT) discriminator_cut1 = 0.898;
  if(type2==e_CSVL) discriminator_cut2 = 0.244;
  if(type2==e_CSVM) discriminator_cut2 = 0.679;
  if(type2==e_CSVT) discriminator_cut2 = 0.898;
  // std::cout << "discriminator_cut1: " <<  discriminator_cut1 << " discriminator_cut2: "<< discriminator_cut1 << std::endl;
  //Create a vector of subjets
  std::vector<Particle> subjets_top;
  //Create a float vector of the subjets discriminators
  std::vector<float> btagsub_combinedSecondaryVertex_top;

  //Fill the vector of subjets with the subjets of a topjet
  subjets_top=topjet.subjets();
  //Fill the vector of discriminators with the discriminators of the subjets of a certain topjet
  btagsub_combinedSecondaryVertex_top=topjet.btagsub_combinedSecondaryVertex();

  //Looping over subjets and checking if they are b-tagged
  for(unsigned int i=0; i < btagsub_combinedSecondaryVertex_top.size(); ++i){
    float test=btagsub_combinedSecondaryVertex_top[i];
    if (nBTagsSub1 != 0 && test>discriminator_cut2) nBTagsSub2 =+ 1;
    if(test>discriminator_cut1){
      if(test>discriminator_cut2 && nBTagsSub2==0) nBTagsSub2+=1;
      else nBTagsSub1 += 1;      //This means it is b-tagged  
 
    }
  }
  if (nBTagsSub1!=0 && nBTagsSub2!=0) return true;
  else return false;
}


bool HepTopTagFull(TopJet topjet){

  //Transform the SFrame TopJet object in a fastjet::PseudoJet
 
  std::vector<fastjet::PseudoJet> jetpart;
  std::vector<Particle> pfconstituents_jet;

  fastjet::ClusterSequence* JetFinder;
  fastjet::JetDefinition* JetDef ;

  pfconstituents_jet=topjet.pfconstituents();
   
  for(unsigned int ic=0; ic<pfconstituents_jet.size(); ++ic){ 
   
    jetpart.push_back( fastjet::PseudoJet(
pfconstituents_jet[ic].pt()*cos(pfconstituents_jet[ic].phi()),
pfconstituents_jet[ic].pt()*sin(pfconstituents_jet[ic].phi()),
pfconstituents_jet[ic].pt()*sinh(pfconstituents_jet[ic].eta()),
pfconstituents_jet[ic].energy() ) );

  }
 
  //Clustering definition
  double conesize=3;
  JetDef = new
fastjet::JetDefinition(fastjet::cambridge_algorithm,conesize); 

  JetFinder = new fastjet::ClusterSequence(jetpart, *JetDef);

  std::vector<fastjet::PseudoJet> tops = JetFinder->inclusive_jets(10.);

  if (tops.size() != 1){
    std::cout << "Problem! Doesn't give exactly one jet!!!!!!!!!!!!!!Gives " << tops.size() << " jets" << std::endl;
    delete JetFinder;
    delete JetDef;
    return false;
  }

  std::vector<fastjet::PseudoJet> SortedJets = sorted_by_pt(tops);


  //Run the HEPTopTagger
  external::HEPTopTagger tagger(*JetFinder, SortedJets[0]);

  //Mass window to be applied in a second step
  tagger.set_top_range(0.0, 10000.0);
  tagger.set_mass_drop_threshold(0.8);
  tagger.set_max_subjet_mass(30);

  tagger.run_tagger();

  delete JetFinder;
  delete JetDef;
 
  if (tagger.is_masscut_passed()) return true;
  else return false;

  return true;
 
}


// 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)
{
  //call variable tagger with default parameters
  return variableHepTopTag(topjet);

}

//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)
{
  //call variable tagger with default parameters
  return variableTopTag(topjet, mjet, nsubjets, mmin);

}
 
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) {
        Jet ji = jets->at(i);
        if (fabs(p->pt() - ji.pt())<1e-8) continue; // skip identical particle
        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);
    if (!nextjet) return ptrel;

    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)
{
    Jet* j = nextJet(p,jets);
    double dr = 999.;
    if (j) dr = j->deltaR(*p);
    return dr;
}

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