Changeset 526 in svn

Timestamp:

Aug 14, 2009, 2:20:19 PM (15 years ago)

Author:

Xavier Rouby

Message:

moving to TRandom3 for random numbers; BTAG_* are now 'float' and not 'int'

Location:

trunk

Files:

• interface/SmearUtil.h (modified) (4 diffs)
• src/SmearUtil.cc (modified) (15 diffs)

trunk/interface/SmearUtil.h

@@ -42 +42 @@
 #include <string>
 #include "TLorentzVector.h"
+#include "TRandom3.h"
 
 #include "D_Constants.h"
@@ -90 +91 @@
   RESOLution(const RESOLution & DET);
   RESOLution& operator=(const RESOLution& DET);
-  ~RESOLution() { delete [] TOWER_eta_edges; delete [] TOWER_dphi;};
+  ~RESOLution() { delete [] TOWER_eta_edges; delete [] TOWER_dphi; delete grandom;};
  
   // Detector coverage for the central detector
@@ -215 +216 @@
  
  //tagging definition
- int BTAG_b;         
- int BTAG_mistag_c;
- int BTAG_mistag_l;
+ float BTAG_b;         
+ float BTAG_mistag_c;
+ float BTAG_mistag_l;
  
 
@@ -236 +237 @@
  PdgTable PDGtable;
  string inputfilelist, detectorcard, triggercard;
+ TRandom3 * grandom;
 
   // to sort a vector

trunk/src/SmearUtil.cc

@@ -37 +37 @@
 
 #include "SmearUtil.h"
-#include "TRandom.h"
 #include "TStopwatch.h"
 
@@ -46 +45 @@
 #include <map>
 #include <vector>
+#include <cmath>
 using namespace std;
 
@@ -140 +140 @@
 
   // Tagging definition
-  BTAG_b           = 40;
-  BTAG_mistag_c    = 10;
-  BTAG_mistag_l    = 1;
+  BTAG_b           = 40.;
+  BTAG_mistag_c    = 10.;
+  BTAG_mistag_l    = 1.;
 
   // FLAGS
@@ -218 +218 @@
   detectorcard      = "";
   triggercard       = "";
+
+  grandom = new TRandom3(0);  // a new seed is set everytime Delphes is run
 }
 
@@ -372 +374 @@
   detectorcard     = DET.detectorcard;
   triggercard      = DET.triggercard;
+
+  grandom          = new TRandom3(*(DET.grandom));
 }
 
@@ -528 +532 @@
   detectorcard     = DET.detectorcard;
   triggercard      = DET.triggercard;
+
+  grandom          = new TRandom3(*(DET.grandom));
 
   return *this;
@@ -658 +664 @@
     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 >> ivalue;BTAG_b            = ivalue;}
-    else if(strstr(temp_string.c_str(),"BTAG_mistag_c"))    {curstring >> varname >> ivalue;BTAG_mistag_c     = ivalue;}
-    else if(strstr(temp_string.c_str(),"BTAG_mistag_l"))    {curstring >> varname >> ivalue;BTAG_mistag_l     = 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(),"FLAG_vfd"))         {curstring >> varname >> ivalue; FLAG_vfd         = ivalue;}
-    else if(strstr(temp_string.c_str(),"FLAG_RP"))          {curstring >> varname >> ivalue; FLAG_RP         = ivalue;}
+    else if(strstr(temp_string.c_str(),"FLAG_RP"))          {curstring >> varname >> ivalue; FLAG_RP          = ivalue;}
     else if(strstr(temp_string.c_str(),"FLAG_trigger"))     {curstring >> varname >> ivalue; FLAG_trigger     = ivalue;}
     else if(strstr(temp_string.c_str(),"FLAG_frog"))        {curstring >> varname >> ivalue; FLAG_frog        = ivalue;}
@@ -1198 +1204 @@
   
   if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
-    energyS = gRandom->Gaus(energy, sqrt( 
+    energyS = grandom->Gaus(energy, sqrt( 
                                          pow(ELG_Ncen,2) + 
                                          pow(ELG_Ccen*energy,2) + 
@@ -1204 +1210 @@
   } 
   else if(fabs(electron.Eta()) >= CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_ec){
-    energyS = gRandom->Gaus(energy, sqrt(
+    energyS = grandom->Gaus(energy, sqrt(
                                          pow(ELG_Nec,2) +
                                          pow(ELG_Cec*energy,2) +
@@ -1210 +1216 @@
   }
   else if(fabs(electron.Eta()) >= CEN_max_calo_ec && fabs(electron.Eta()) < CEN_max_calo_fwd){
-    energyS = gRandom->Gaus(energy, sqrt(
+    energyS = grandom->Gaus(energy, sqrt(
                                          pow(ELG_Nfwd,2) +
                                          pow(ELG_Cfwd*energy,2) +
@@ -1229 +1235 @@
   if(fabs(muon.Eta()) < CEN_max_mu )
    {
-     ptS = gRandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
+     ptS = grandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
    }
      muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
@@ -1252 +1258 @@
   float energyS1,energyS2;
   if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
-   energyS1 = gRandom->Gaus(energy_hcal, sqrt(
+   energyS1 = grandom->Gaus(energy_hcal, sqrt(
                                               pow(HAD_Ncen,2) +
                                               pow(HAD_Ccen*energy_hcal,2) + 
                                               pow(HAD_Scen*sqrt(energy_hcal),2) )) ;
       
-   energyS2 =  gRandom->Gaus(energy_ecal, sqrt(
+   energyS2 =  grandom->Gaus(energy_ecal, sqrt(
                                       pow(ELG_Ncen,2) +
                                       pow(ELG_Ccen*energy_ecal,2) + 
@@ -1265 +1271 @@
   } 
   else if(fabs(hadron.Eta()) >= CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_ec) {
-   energyS1 = gRandom->Gaus(energy_hcal, sqrt(
+   energyS1 = grandom->Gaus(energy_hcal, sqrt(
                                               pow(HAD_Nec,2) +
                                               pow(HAD_Cec*energy_hcal,2) +
                                               pow(HAD_Sec*sqrt(energy_hcal),2) )) ;
 
-   energyS2 =  gRandom->Gaus(energy_ecal, sqrt(
+   energyS2 =  grandom->Gaus(energy_ecal, sqrt(
                                       pow(ELG_Nec,2) +
                                       pow(ELG_Cec*energy_ecal,2) +
@@ -1278 +1284 @@
   }
   else if(fabs(hadron.Eta()) >= CEN_max_calo_ec && fabs(hadron.Eta()) < CEN_max_calo_fwd){ 
-    energyS = gRandom->Gaus(energy, sqrt(
+    energyS = grandom->Gaus(energy, sqrt(
                             pow(HAD_Nfwd,2) +
                             pow(HAD_Cfwd*energy,2) +
@@ -1373 +1379 @@
 
 //******************** 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) {
-  if(      rand()%100 < (BTAG_b+1)    && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
-  else if( rand()%100 < (BTAG_mistag_c+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
-  else if( rand()%100 < (BTAG_mistag_l+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0)   return true; // b-tag of light jets is 1%
-  return false;
-}
+  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;
+}
+
 
 //***********************Isolation criteria***********************