[ViewVC] Diff of: cvsroot/UserCode/UHHAnalysis/SFrameTools/src/Utils.cxx

Comparing UserCode/UHHAnalysis/SFrameTools/src/Utils.cxx (file contents):
Revision 1.5 by peiffer, Mon Jul 9 08:37:32 2012 UTC vs.
Revision 1.7 by peiffer, Fri Dec 7 10:56:32 2012 UTC

+#include "../include/Utils.h"
-+
+//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();
+  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;
-+
-+
+}

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing UserCode/UHHAnalysis/SFrameTools/src/Utils.cxx (file contents): Revision 1.5 by peiffer, Mon Jul 9 08:37:32 2012 UTC vs. Revision 1.7 by peiffer, Fri Dec 7 10:56:32 2012 UTC

Diff Legend

Comparing UserCode/UHHAnalysis/SFrameTools/src/Utils.cxx (file contents):
Revision 1.5 by peiffer, Mon Jul 9 08:37:32 2012 UTC vs.
Revision 1.7 by peiffer, Fri Dec 7 10:56:32 2012 UTC