Changeset 558 in svn for trunk/src/SmearUtil.cc

Timestamp:

Feb 24, 2010, 6:26:47 PM (15 years ago)

Author:

Xavier Rouby

Message:

new b-tag

File:

• trunk/src/SmearUtil.cc (modified) (12 diffs)

trunk/src/SmearUtil.cc

@@ -38 +38 @@
 #include "SmearUtil.h"
 #include "TStopwatch.h"
+#include "TF2.h"
 
 #include <iostream>
@@ -142 +143 @@
 
   // Tagging definition
-  BTAG_b           = 40.;
-  BTAG_mistag_c    = 10.;
-  BTAG_mistag_l    = 1.;
+  BTAG_func_b      = "0.4  + 0*x + 0*y";        // efficiency functions (x=Pt, y=eta)  in ROOT::TF2 format
+  BTAG_func_c      = "0.1  + 0*x + 0*y";
+  BTAG_func_l      = "0.01 + 0*x + 0*y";
+  BTAG_func_g      = "0.01 + 0*x + 0*y";
+
 
   // FLAGS
@@ -306 +309 @@
 
   // Tagging definition
-  BTAG_b           = DET.BTAG_b;
-  BTAG_mistag_c    = DET.BTAG_mistag_c;
-  BTAG_mistag_l    = DET.BTAG_mistag_l;
+  BTAG_func_b      = DET.BTAG_func_b;
+  BTAG_func_c      = DET.BTAG_func_c;
+  BTAG_func_l      = DET.BTAG_func_l;
+  BTAG_func_g      = DET.BTAG_func_g;
+
 
   // FLAGS
@@ -462 +467 @@
 
   // Tagging definition
-  BTAG_b           = DET.BTAG_b;
-  BTAG_mistag_c    = DET.BTAG_mistag_c;
-  BTAG_mistag_l    = DET.BTAG_mistag_l;
+  BTAG_func_b      = DET.BTAG_func_b;
+  BTAG_func_c      = DET.BTAG_func_c;
+  BTAG_func_l      = DET.BTAG_func_l;
+  BTAG_func_g      = DET.BTAG_func_g;
 
   // FLAGS
@@ -669 +675 @@
     else if(strstr(temp_string.c_str(),"JET_Eflow"))        {curstring >> varname >> ivalue; JET_Eflow        = ivalue;}
  
-    else if(strstr(temp_string.c_str(),"BTAG_b"))           {curstring >> varname >> value;BTAG_b            = value;}
-    else if(strstr(temp_string.c_str(),"BTAG_mistag_c"))    {curstring >> varname >> value;BTAG_mistag_c     = value;}
-    else if(strstr(temp_string.c_str(),"BTAG_mistag_l"))    {curstring >> varname >> value;BTAG_mistag_l     = value;}
+    else if(strstr(temp_string.c_str(),"BTAG_func_b"))      {curstring >> varname >> svalue; BTAG_func_b     = svalue;}
+    else if(strstr(temp_string.c_str(),"BTAG_func_c"))      {curstring >> varname >> svalue; BTAG_func_c     = svalue;}
+    else if(strstr(temp_string.c_str(),"BTAG_func_l"))      {curstring >> varname >> svalue; BTAG_func_l     = svalue;}
+    else if(strstr(temp_string.c_str(),"BTAG_func_g"))      {curstring >> varname >> svalue; BTAG_func_g     = svalue;}
 
     else if(strstr(temp_string.c_str(),"FLAG_vfd"))         {curstring >> varname >> ivalue; FLAG_vfd         = ivalue;}
@@ -1185 +1192 @@
         << left << setw(5) <<TAU_track_pt                                <<""<< right << setw(20)<<"*"<<"\n";
   f_out<<"*                                                                    *"<<"\n";
-  f_out<<"#***************************                                         *"<<"\n";
-  f_out<<"# B-tagging efficiencies [%]                                         *"<<"\n";
-  f_out<<"#***************************                                         *"<<"\n";
-  f_out<<"*                                                                    *"<<"\n";
-  f_out << left << setw(50) <<"* Efficiency to tag a \"b\" as a b-jet: "<<""
-        << left << setw(10) <<BTAG_b               <<""<< right << setw(10)<<"*"<<"\n";
-  f_out << left << setw(50) <<"* Efficiency to mistag a c-jet as a b-jet: "<<""
-        << left << setw(10) <<BTAG_mistag_c            <<""<< right << setw(10)<<"*"<<"\n";
-  f_out << left << setw(50) <<"* Efficiency to mistag a light jet as a b-jet: "<<""
-        << left << setw(10) <<BTAG_mistag_l            <<""<< right << setw(10)<<"*"<<"\n";
+  f_out<<"#*******************************                                     *"<<"\n";
+  f_out<<"# B-tagging efficiency functions                                     *"<<"\n";
+  f_out<<"#*******************************                                     *"<<"\n";
+  f_out<<"*                                                                    *"<<"\n";
+  f_out << left << setw(50) <<"* Tagging a \"b\" as a b-jet: "<<""
+        << left << setw(18) <<BTAG_func_b            <<""<< right << setw(2)<<"*"<<"\n";
+  f_out << left << setw(50) <<"* Mistagging a c-jet as a b-jet: "<<""
+        << left << setw(18) <<BTAG_func_c            <<""<< right << setw(2)<<"*"<<"\n";
+  f_out << left << setw(50) <<"* Mistagging a gluon-jet as a b-jet: "<<""
+        << left << setw(18) <<BTAG_func_g            <<""<< right << setw(2)<<"*"<<"\n";
+  f_out << left << setw(50) <<"* Mistagging a light jet as a b-jet: "<<""
+        << left << setw(18) <<BTAG_func_l            <<""<< right << setw(2)<<"*"<<"\n";
+
+
   f_out<<"*                                                                    *"<<"\n";
   f_out<<"*                                                                    *"<<"\n";
@@ -1386 +1397 @@
 
 
-//******************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
+//******************* Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
 
 bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootC::GenParticle> &subarray) {
-  float chance =grandom->Uniform()*100.;
-
-  /* old code -- for bookkeeping 
-  bool res1 = false;
-  if(      chance < BTAG_b        && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) { res1= true;} // b-tag of b-jets is 40%
-  else if( chance < BTAG_mistag_c && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) { res1= true;} // b-tag of c-jets is 10%
-  else if( chance < BTAG_mistag_l) {
-    int pid=Bjets(subarray,JET.Eta(),JET.Phi());
-    if(pid!=pB && pid!=pC && pid!=0) { res1=true; }// b-tag of light jets is 1%
-  }
-  else res1 = false;
-  //return res1;
-  */ 
-
-
-  bool res2 = false;   int jet;
-  if( chance > max(BTAG_b, max(BTAG_mistag_c, BTAG_mistag_l)) ) { // useless to try further
-        res2 = false;
-  } 
-  else {
-        jet = Bjets(subarray,JET.Eta(),JET.Phi()); // once for all
-        if(      chance < BTAG_b        && jet==pB ) { res2= true;} // b-tag of b-jets is 40%
-        else if( chance < BTAG_mistag_c && jet==pC ) { res2= true;} // b-tag of c-jets is 10%
-        else if( chance < BTAG_mistag_l && jet!=pB && jet!=pC && jet!=0) { res2=true; }// b-tag of light jets is 1%
-        else res2= false;
-  }
-  //if (res1 != res2) cout << "BUUUUUUUUUUUUUUGGGGGGG" << endl;
-
-  return res2;
+
+  const float ptmax = 9E6; // GeV/c
+  bool tag = false;
+  string efficiency_formula="";
+
+  // selects the correct formula
+  int jet_id = Bjets(subarray,JET.Eta(),JET.Phi()); // gets the particle id of the most energetic parton in the jet
+  switch (jet_id) {
+    case pB: efficiency_formula = BTAG_func_b; break;
+    case pC: efficiency_formula = BTAG_func_c; break;
+    case pGLUON: efficiency_formula = BTAG_func_g; break;
+    default: efficiency_formula = BTAG_func_l; break;
+  }
+
+  TF2 efficiency("efficiency",efficiency_formula.c_str(),0,ptmax,-CEN_max_tracker,CEN_max_tracker); // efficiency function wrt Pt
+  float x = grandom->Uniform(); // randomisation
+  tag = (x < efficiency.Eval(JET.Pt(),JET.Eta())) ? true : false;
+  //cout << "formula = " << efficiency_formula << "\t Pt = " << JET.Pt() 
+  //     << "\t value = " << efficiency.Eval(JET.Pt())  
+  //     << "\t x = " << x << "\t tag= " << tag << endl;
+
+  return tag;
 }
 
@@ -1450 +1454 @@
    float etiso = 0; // sum of all track pt in isolation cone
    for (unsigned int i=0; i< towers.size(); i++) {
-//       cout<<" eta part "<<iEta<<" eta tour "<<towers[i].getEta()<<endl;
-//       cout<<" phi part "<<iPhi<<" phi tour "<<towers[i].getPhi()<<endl;
-
-        if(DeltaR(iPhi,iEta,towers[i].getPhi(),towers[i].getEta()) < ISOL_calo_Cone) {
+        if( (DeltaR(iPhi,iEta,towers[i].getPhi(),towers[i].getEta()) < ISOL_calo_Cone) && (towers[i].getET() > ISOL_calo_ET) ) {
            etiso += towers[i].getET(); // dR cut on tracks
           }
   }
-
    return etiso;
 }
@@ -1463 +1463 @@
 
 
-// ******* Calorimetric isolation
+// ******* Calorimetric isolation -- for LHCO only
 float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower) {
   // iPhi and iEta are the coordinates of the center of the tower crossed by the particle
@@ -1469 +1469 @@
   // iTower is the index of the tower, in [0, n_tower]. iTower points only towers in positive range
 
-   if(ISOL_calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
    float et_sum=0;
 
@@ -1504 +1503 @@
             //     << "calMuon.getEta= " << calMuon.getEta() << "\tcalMuon.getPhi()= " << calMuon.getPhi() <<"\t";
 
-            if(calMuon.getEta() != UNDEFINED && calMuon.getET() > ISOL_calo_ET) {
+            //if(calMuon.getEta() != UNDEFINED && calMuon.getET() > ISOL_calo_ET) {
+            if(calMuon.getEta() != UNDEFINED && calMuon.getET() > 0.0) {
                et_sum += calMuon.getET();
                //cout << calMuon.getET() << " GeV";
@@ -1513 +1513 @@
    }  // undefined index
    else {
-     cout << "out of range" << endl;
      if (CEN_max_mu < CEN_max_calo_fwd) 
         cout << "**    ERROR in RESOLution::CaloIsolation: 'muon'-tower not found!  **" << endl; 
    }      // should never happen ! this would be a bug
-  
   return (part.Pt()==0)? UNDEFINED: et_sum/part.Pt();
 }