Fork me on GitHub

Changeset 526 in svn


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/interface/SmearUtil.h

    r494 r526  
    4242#include <string>
    4343#include "TLorentzVector.h"
     44#include "TRandom3.h"
    4445
    4546#include "D_Constants.h"
     
    9091  RESOLution(const RESOLution & DET);
    9192  RESOLution& operator=(const RESOLution& DET);
    92   ~RESOLution() { delete [] TOWER_eta_edges; delete [] TOWER_dphi;};
     93  ~RESOLution() { delete [] TOWER_eta_edges; delete [] TOWER_dphi; delete grandom;};
    9394 
    9495  // Detector coverage for the central detector
     
    215216 
    216217 //tagging definition
    217  int BTAG_b;         
    218  int BTAG_mistag_c;
    219  int BTAG_mistag_l;
     218 float BTAG_b;         
     219 float BTAG_mistag_c;
     220 float BTAG_mistag_l;
    220221 
    221222
     
    236237 PdgTable PDGtable;
    237238 string inputfilelist, detectorcard, triggercard;
     239 TRandom3 * grandom;
    238240
    239241  // to sort a vector
  • trunk/src/SmearUtil.cc

    r494 r526  
    3737
    3838#include "SmearUtil.h"
    39 #include "TRandom.h"
    4039#include "TStopwatch.h"
    4140
     
    4645#include <map>
    4746#include <vector>
     47#include <cmath>
    4848using namespace std;
    4949
     
    140140
    141141  // Tagging definition
    142   BTAG_b           = 40;
    143   BTAG_mistag_c    = 10;
    144   BTAG_mistag_l    = 1;
     142  BTAG_b           = 40.;
     143  BTAG_mistag_c    = 10.;
     144  BTAG_mistag_l    = 1.;
    145145
    146146  // FLAGS
     
    218218  detectorcard      = "";
    219219  triggercard       = "";
     220
     221  grandom = new TRandom3(0);  // a new seed is set everytime Delphes is run
    220222}
    221223
     
    372374  detectorcard     = DET.detectorcard;
    373375  triggercard      = DET.triggercard;
     376
     377  grandom          = new TRandom3(*(DET.grandom));
    374378}
    375379
     
    528532  detectorcard     = DET.detectorcard;
    529533  triggercard      = DET.triggercard;
     534
     535  grandom          = new TRandom3(*(DET.grandom));
    530536
    531537  return *this;
     
    658664    else if(strstr(temp_string.c_str(),"JET_Eflow"))        {curstring >> varname >> ivalue; JET_Eflow        = ivalue;}
    659665 
    660     else if(strstr(temp_string.c_str(),"BTAG_b"))           {curstring >> varname >> ivalue;BTAG_b            = ivalue;}
    661     else if(strstr(temp_string.c_str(),"BTAG_mistag_c"))    {curstring >> varname >> ivalue;BTAG_mistag_c     = ivalue;}
    662     else if(strstr(temp_string.c_str(),"BTAG_mistag_l"))    {curstring >> varname >> ivalue;BTAG_mistag_l     = ivalue;}
     666    else if(strstr(temp_string.c_str(),"BTAG_b"))           {curstring >> varname >> value;BTAG_b            = value;}
     667    else if(strstr(temp_string.c_str(),"BTAG_mistag_c"))    {curstring >> varname >> value;BTAG_mistag_c     = value;}
     668    else if(strstr(temp_string.c_str(),"BTAG_mistag_l"))    {curstring >> varname >> value;BTAG_mistag_l     = value;}
    663669
    664670    else if(strstr(temp_string.c_str(),"FLAG_vfd"))         {curstring >> varname >> ivalue; FLAG_vfd         = ivalue;}
    665     else if(strstr(temp_string.c_str(),"FLAG_RP"))          {curstring >> varname >> ivalue; FLAG_RP         = ivalue;}
     671    else if(strstr(temp_string.c_str(),"FLAG_RP"))          {curstring >> varname >> ivalue; FLAG_RP          = ivalue;}
    666672    else if(strstr(temp_string.c_str(),"FLAG_trigger"))     {curstring >> varname >> ivalue; FLAG_trigger     = ivalue;}
    667673    else if(strstr(temp_string.c_str(),"FLAG_frog"))        {curstring >> varname >> ivalue; FLAG_frog        = ivalue;}
     
    11981204 
    11991205  if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
    1200     energyS = gRandom->Gaus(energy, sqrt(
     1206    energyS = grandom->Gaus(energy, sqrt(
    12011207                                         pow(ELG_Ncen,2) +
    12021208                                         pow(ELG_Ccen*energy,2) +
     
    12041210  }
    12051211  else if(fabs(electron.Eta()) >= CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_ec){
    1206     energyS = gRandom->Gaus(energy, sqrt(
     1212    energyS = grandom->Gaus(energy, sqrt(
    12071213                                         pow(ELG_Nec,2) +
    12081214                                         pow(ELG_Cec*energy,2) +
     
    12101216  }
    12111217  else if(fabs(electron.Eta()) >= CEN_max_calo_ec && fabs(electron.Eta()) < CEN_max_calo_fwd){
    1212     energyS = gRandom->Gaus(energy, sqrt(
     1218    energyS = grandom->Gaus(energy, sqrt(
    12131219                                         pow(ELG_Nfwd,2) +
    12141220                                         pow(ELG_Cfwd*energy,2) +
     
    12291235  if(fabs(muon.Eta()) < CEN_max_mu )
    12301236   {
    1231      ptS = gRandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
     1237     ptS = grandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
    12321238   }
    12331239     muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
     
    12521258  float energyS1,energyS2;
    12531259  if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
    1254    energyS1 = gRandom->Gaus(energy_hcal, sqrt(
     1260   energyS1 = grandom->Gaus(energy_hcal, sqrt(
    12551261                                              pow(HAD_Ncen,2) +
    12561262                                              pow(HAD_Ccen*energy_hcal,2) +
    12571263                                              pow(HAD_Scen*sqrt(energy_hcal),2) )) ;
    12581264     
    1259    energyS2 =  gRandom->Gaus(energy_ecal, sqrt(
     1265   energyS2 =  grandom->Gaus(energy_ecal, sqrt(
    12601266                                      pow(ELG_Ncen,2) +
    12611267                                      pow(ELG_Ccen*energy_ecal,2) +
     
    12651271  }
    12661272  else if(fabs(hadron.Eta()) >= CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_ec) {
    1267    energyS1 = gRandom->Gaus(energy_hcal, sqrt(
     1273   energyS1 = grandom->Gaus(energy_hcal, sqrt(
    12681274                                              pow(HAD_Nec,2) +
    12691275                                              pow(HAD_Cec*energy_hcal,2) +
    12701276                                              pow(HAD_Sec*sqrt(energy_hcal),2) )) ;
    12711277
    1272    energyS2 =  gRandom->Gaus(energy_ecal, sqrt(
     1278   energyS2 =  grandom->Gaus(energy_ecal, sqrt(
    12731279                                      pow(ELG_Nec,2) +
    12741280                                      pow(ELG_Cec*energy_ecal,2) +
     
    12781284  }
    12791285  else if(fabs(hadron.Eta()) >= CEN_max_calo_ec && fabs(hadron.Eta()) < CEN_max_calo_fwd){
    1280     energyS = gRandom->Gaus(energy, sqrt(
     1286    energyS = grandom->Gaus(energy, sqrt(
    12811287                            pow(HAD_Nfwd,2) +
    12821288                            pow(HAD_Cfwd*energy,2) +
     
    13731379
    13741380//******************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets****************
     1381
    13751382bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootC::GenParticle> &subarray) {
    1376   if(      rand()%100 < (BTAG_b+1)    && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) return true; // b-tag of b-jets is 40%
    1377   else if( rand()%100 < (BTAG_mistag_c+1) && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) return true; // b-tag of c-jets is 10%
    1378   else if( rand()%100 < (BTAG_mistag_l+1) && Bjets(subarray,JET.Eta(),JET.Phi())!=0)   return true; // b-tag of light jets is 1%
    1379   return false;
    1380 }
     1383  float chance =grandom->Uniform()*100.;
     1384
     1385  /* old code -- for bookkeeping
     1386  bool res1 = false;
     1387  if(      chance < BTAG_b        && Bjets(subarray,JET.Eta(),JET.Phi())==pB ) { res1= true;} // b-tag of b-jets is 40%
     1388  else if( chance < BTAG_mistag_c && Bjets(subarray,JET.Eta(),JET.Phi())==pC ) { res1= true;} // b-tag of c-jets is 10%
     1389  else if( chance < BTAG_mistag_l) {
     1390    int pid=Bjets(subarray,JET.Eta(),JET.Phi());
     1391    if(pid!=pB && pid!=pC && pid!=0) { res1=true; }// b-tag of light jets is 1%
     1392  }
     1393  else res1 = false;
     1394  //return res1;
     1395  */
     1396
     1397
     1398  bool res2 = false;   int jet;
     1399  if( chance > max(BTAG_b, max(BTAG_mistag_c, BTAG_mistag_l)) ) { // useless to try further
     1400        res2 = false;
     1401  }
     1402  else {
     1403        jet = Bjets(subarray,JET.Eta(),JET.Phi()); // once for all
     1404        if(      chance < BTAG_b        && jet==pB ) { res2= true;} // b-tag of b-jets is 40%
     1405        else if( chance < BTAG_mistag_c && jet==pC ) { res2= true;} // b-tag of c-jets is 10%
     1406        else if( chance < BTAG_mistag_l && jet!=pB && jet!=pC && jet!=0) { res2=true; }// b-tag of light jets is 1%
     1407        else res2= false;
     1408  }
     1409  //if (res1 != res2) cout << "BUUUUUUUUUUUUUUGGGGGGG" << endl;
     1410
     1411  return res2;
     1412}
     1413
    13811414
    13821415//***********************Isolation criteria***********************
Note: See TracChangeset for help on using the changeset viewer.