/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/DQM/Physics/src/EwkMuLumiMonitorDQM.cc

Go to the documentation of this file.

#include "DQM/Physics/src/EwkMuLumiMonitorDQM.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/Common/interface/Handle.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"



#include "DataFormats/Common/interface/View.h"
#include "DataFormats/Common/interface/Handle.h"

#include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"

#include "FWCore/Common/interface/TriggerNames.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/HLTReco/interface/TriggerObject.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/Math/interface/deltaR.h"

#include "DataFormats/METReco/interface/MET.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"

#include "DataFormats/Math/interface/LorentzVector.h"

  
using namespace edm;
using namespace std;
using namespace reco;
using namespace isodeposit;

EwkMuLumiMonitorDQM::EwkMuLumiMonitorDQM( const ParameterSet & cfg ) :
  // Input collections
   trigTag_(cfg.getUntrackedParameter<edm::InputTag> ("TrigTag", edm::InputTag("TriggerResults::HLT"))),
  trigEv_(cfg.getUntrackedParameter<edm::InputTag> ("triggerEvent")),
  muonTag_(cfg.getUntrackedParameter<edm::InputTag>("muons")),
  trackTag_(cfg.getUntrackedParameter<edm::InputTag>("tracks")),
  caloTowerTag_(cfg.getUntrackedParameter<edm::InputTag>("calotower")),
  metTag_(cfg.getUntrackedParameter<edm::InputTag> ("metTag")),
  metIncludesMuons_(cfg.getUntrackedParameter<bool> ("METIncludesMuons")),
 // Main cuts 
  // massMin_(cfg.getUntrackedParameter<double>("MtMin", 20.)),
  //   massMax_(cfg.getUntrackedParameter<double>("MtMax", 2000.)) 
   //  hltPath_(cfg.getUntrackedParameter<std::string> ("hltPath")) ,
   //  L3FilterName_(cfg.getUntrackedParameter<std::string> ("L3FilterName")),           
  ptMuCut_(cfg.getUntrackedParameter<double>("ptMuCut")),
  etaMuCut_(cfg.getUntrackedParameter<double>("etaMuCut")),
  isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso")),
  isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso")),
  isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03")),
  //  deltaRTrk_(cfg.getUntrackedParameter<double>("deltaRTrk")),
  ptThreshold_(cfg.getUntrackedParameter<double>("ptThreshold")),
  //deltaRVetoTrk_(cfg.getUntrackedParameter<double>("deltaRVetoTrk")), 
  maxDPtRel_(cfg.getUntrackedParameter<double>("maxDPtRel")),
  maxDeltaR_(cfg.getUntrackedParameter<double>("maxDeltaR")),
  mtMin_(cfg.getUntrackedParameter<double>("mtMin")),
  mtMax_(cfg.getUntrackedParameter<double>("mtMax")),
  acopCut_(cfg.getUntrackedParameter<double>("acopCut")), 
  dxyCut_(cfg.getUntrackedParameter<double>("DxyCut")) {
 // just to initialize
  isValidHltConfig_ = false;

}



void EwkMuLumiMonitorDQM::beginRun(const Run& r, const EventSetup&iSetup) {
      nall =  0;
      nEvWithHighPtMu=0; 
      nInKinRange=  0;
      nsel = 0;
      niso = 0;
      nhlt =0;
      n1hlt =0;
      n2hlt =0;
      nNotIso =0;
      nGlbSta =0;
      nGlbTrk =0;
      nTMass =0;
      nW=0;

  // passed as parameter to HLTConfigProvider::init(), not yet used
  bool isConfigChanged = false;

  // isValidHltConfig_ used to short-circuit analyze() in case of problems
  isValidHltConfig_ = hltConfigProvider_.init( r, iSetup, trigTag_.process(), isConfigChanged );
  //std::cout << "hlt config trigger is valid??" << isValidHltConfig_ << std::endl; 

  }


void EwkMuLumiMonitorDQM::beginJob() {
      theDbe = Service<DQMStore>().operator->();
      theDbe->setCurrentFolder("Physics/EwkMuLumiMonitorDQM");

      init_histograms();
}

void EwkMuLumiMonitorDQM::init_histograms() {

      char chtitle[256] = "";
      for (int i=0; i<2; ++i) {
            snprintf(chtitle, 255, "Z mass [GeV/c^{2}]");
            mass2HLT_ = theDbe->book1D("Z_2HLT_MASS",chtitle,200,0.,200.);
            mass1HLT_ = theDbe->book1D("Z_1HLT_MASS",chtitle,200,0.,200.);
            massNotIso_ = theDbe->book1D("Z_NOTISO_MASS",chtitle,200,0.,200.);
            massGlbSta_ = theDbe->book1D("Z_GLBSTA_MASS",chtitle,200,0.,200.);
            massGlbTrk_ = theDbe->book1D("Z_GLBTRK_MASS",chtitle,200,0.,200.);
            massIsBothGlbTrkThanW_ = theDbe->book1D("Z_ISBOTHGLBTRKTHANW_MASS",chtitle,200,0.,200.);
          

            snprintf(chtitle, 255, "Z high mass [GeV/c^{2}]");
            highMass2HLT_ = theDbe->book1D("Z_2HLT_HIGHMASS",chtitle,2000,0.,2000.);
            highMass1HLT_ = theDbe->book1D("Z_1HLT_HIGHMASS",chtitle,2000,0.,2000.);              
            highMassNotIso_ = theDbe->book1D("Z_NOTISO_HIGHMASS",chtitle,2000,0.,2000.);              
            highMassGlbSta_ = theDbe->book1D("Z_GLBSTA_HIGHMASS",chtitle,2000,0.,2000.);              
            highMassGlbTrk_ = theDbe->book1D("Z_GLBTRK_HIGHMASS",chtitle,2000,0.,2000.);              
            highMassIsBothGlbTrkThanW_ = theDbe->book1D("Z_ISBOTHGLBTRKTHANW_HIGHMASS",chtitle,2000,0.,2000.);  
 
            /*      snprintf(chtitle, 255, "Highest muon p_{t}[GeV/c]");
            highest_mupt2HLT_ = theDbe->book1D("HIGHEST_MU_PT_2HLT",chtitle,200,0.,200.);
            highest_mupt1HLT_ = theDbe->book1D("HIGHEST_MU_PT_1HLT",chtitle,200,0.,200.);
            highest_muptNotIso_ = theDbe->book1D("HIGHEST_MU_PT_NOTISO",chtitle,200,0.,200.);
            highest_muptGlbSta_ = theDbe->book1D("HIGHEST_MU_PT_GLBSTA",chtitle,200,0.,200.);
            highest_muptGlbTrk_ = theDbe->book1D("HIGHEST_MU_PT_GLBTRK",chtitle,200,0.,200.);

            snprintf(chtitle, 255, "Lowest muon p_{t} [GeV/c]");
            lowest_mupt2HLT_ = theDbe->book1D("LOWEST_MU_PT_2HLT",chtitle,200,0.,200.);
            lowest_mupt1HLT_ = theDbe->book1D("LOWEST_MU_PT_1HLT",chtitle,200,0.,200.);
            lowest_muptNotIso_ = theDbe->book1D("LOWEST_MU_PT_NOTISO",chtitle,200,0.,200.);
            lowest_muptGlbSta_ = theDbe->book1D("LOWEST_MU_PT_GLBSTA",chtitle,200,0.,200.);
            lowest_muptGlbTrk_ = theDbe->book1D("LOWEST_MU_PT_GLBTRK",chtitle,200,0.,200.);
            */

            snprintf(chtitle, 255, "Transverse mass [GeV]" );
            TMass_ = theDbe->book1D("TMASS",chtitle,300,0.,300.);
            
      }
}



double EwkMuLumiMonitorDQM::muIso( const reco::Muon & mu )  {        
    double isovar = mu.isolationR03().sumPt;
    if (isCombinedIso_) {
      isovar += mu.isolationR03().emEt;
      isovar += mu.isolationR03().hadEt;
    }
    if (isRelativeIso_) isovar /= mu.pt();
    return isovar;
}

double EwkMuLumiMonitorDQM::tkIso( const reco::Track tk,  Handle< TrackCollection > tracks , Handle<CaloTowerCollection> calotower  )  {        
      double ptSum = 0;
      for (size_t i=0; i< tracks->size(); ++i){ 
        const reco::Track & elem = tracks->at(i);
        double elemPt = elem.pt();
        // same parameter used for muIsolation: dR [0.01, IsoCut03_], |dZ|<0.2, |d_r(xy)|<0.1
       double elemVx = elem.vx();
       double elemVy = elem.vy();
       double elemD0 = sqrt( elemVx * elemVx + elemVy * elemVy );
       if ( elemD0 > 0.2 ) continue;
       double dz = fabs( elem.vz() - tk.vz()); 
       if ( dz > 0.1) continue;
        // evaluate only for tracks with pt>ptTreshold
        if ( elemPt <  ptThreshold_ ) continue;
        double dR = deltaR( elem.eta(), elem.phi(), tk.eta(), tk.phi() );
        // isolation in a cone with dR=0.3, and vetoing the track itself 
        if ( (dR < 0.01) || (dR > 0.3)  )  continue;
        ptSum += elemPt;
      }
      if (isCombinedIso_) {
        // loop on clusters....
        for (CaloTowerCollection::const_iterator it=calotower->begin();  it!=calotower->end();++it ){
          double dR = deltaR( it->eta(), it->phi(),tk.outerEta(), tk.outerPhi() );
          // veto value is 0.1 for towers....
          if ( (dR < 0.1) || (dR > 0.3)  )  continue;
          ptSum += it->emEnergy();
          ptSum += it->hadEnergy();
        }
      }
    if (isRelativeIso_) ptSum /= tk.pt();
    return ptSum;
}



bool EwkMuLumiMonitorDQM::IsMuMatchedToHLTMu ( const reco::Muon & mu, std::vector<reco::Particle> HLTMu , double DR, double DPtRel ) {
  size_t dim =  HLTMu.size();
  size_t nPass=0;
  if (dim==0) return false;
  for (size_t k =0; k< dim; k++ ) {
    if (  (deltaR(HLTMu[k], mu) < DR)   && (fabs(HLTMu[k].pt() - mu.pt())/ HLTMu[k].pt()<DPtRel)){     
      nPass++ ;
    }
  }
  return (nPass>0);
}



void EwkMuLumiMonitorDQM::endJob() {
}

void EwkMuLumiMonitorDQM::endRun(const Run& r, const EventSetup&) {
  
  LogVerbatim("") << "\n>>>>>> Z/W SELECTION SUMMARY BEGIN >>>>>>>>>>>>>>>";
  LogVerbatim("") << "Total numer of events analyzed: " << nall << " [events]";
  LogVerbatim("") << "Total numer of events selected: " << nsel << " [events]";
  
  LogVerbatim("") << "Passing HLT criteria:           " << nhlt << " [events] " ; 
  LogVerbatim("") << "Passing 2 HLT match criteria:           " << n2hlt << " [events] " ;
  LogVerbatim("") << "Passing 1 HLT match criteria:           " << n1hlt << " [events] " ;
  LogVerbatim("") << "Passing Not Iso criteria:           " << nNotIso << " [events] " ;
  LogVerbatim("") << "Passing GlbSta criteria:           " << nGlbSta << " [events] " ;
  LogVerbatim("") << "Passing W criteria:           " << nTMass << " [events] " ;
  LogVerbatim("") << ">>>>>> Z/W SELECTION SUMMARY END   >>>>>>>>>>>>>>>\n";
}

void EwkMuLumiMonitorDQM::analyze (const Event & ev, const EventSetup &) {
      nall++;  
      bool iso_sel = false;
      bool hlt_sel = false;
      double iso1=-1;
      double iso2=-1;
      bool isMu1Iso = false;
      bool isMu2Iso = false;
      bool singleTrigFlag1 = false;
      bool singleTrigFlag2 = false; 
      isZGolden1HLT_=false;  
      isZGolden2HLT_=false; 
      isZGoldenNoIso_=false;
      isZGlbSta_=false;
      isZGlbTrk_=false;      
      isW_=false;
      // Trigger
      bool trigger_fired = false;
      
      Handle<TriggerResults> triggerResults;
      if (!ev.getByLabel(trigTag_, triggerResults)) {
        //LogWarning("") << ">>> TRIGGER collection does not exist !!!";
        return;
      }
      
      ev.getByLabel(trigTag_, triggerResults); 
      /*
        const edm::TriggerNames & trigNames = ev.triggerNames(*triggerResults);

        
      for (size_t i=0; i<triggerResults->size(); i++) {
      std::string trigName = trigNames.triggerName(i);
      //std::cout << " trigName == " << trigName << std::endl;
        if ( trigName == hltPath_ && triggerResults->accept(i)) {
        trigger_fired = true;
        hlt_sel=true;
        nhlt++;
        }
        } 
      */
      // see the trigger single muon which are present
      string lowestMuonUnprescaledTrig = "";
      bool lowestMuonUnprescaledTrigFound = false;
      const std::vector<std::string>& triggerNames = hltConfigProvider_.triggerNames();
      for (size_t ts = 0; ts< triggerNames.size() ; ts++){
        string trig = triggerNames[ts];
        size_t f = trig.find("HLT_Mu");
        if ( (f != std::string::npos)  ) {
          // std::cout << "single muon trigger present: " << trig << std::endl; 
          // See if the trigger is prescaled; 
          bool prescaled = false;
          const unsigned int prescaleSize = hltConfigProvider_.prescaleSize() ;
          for (unsigned int ps= 0; ps<  prescaleSize; ps++){
            const unsigned int prescaleValue = hltConfigProvider_.prescaleValue(ps, trig) ;
            if (prescaleValue != 1) prescaled =true;
            //  std::cout<< " prescaleValue[" << ps << "] =" << prescaleValue <<std::endl;
          }
          if (!prescaled){
            //looking now for the lowest hlt path not prescaled, with name of the form HLT_MuX or HLTMuX_vY
            for (unsigned int n=9; n<100 ; n++ ){
              string lowestTrig= "HLT_Mu";
              string lowestTrigv0 = "copy";
              std::stringstream out;
              out << n;
              std::string s = out.str(); 
              lowestTrig.append(s);
              lowestTrigv0 = lowestTrig;  
              for (unsigned int v = 1; v<10 ; v++ ){
                lowestTrig.append("_v");
                std::stringstream oout;
                oout << v;
                std::string ss = oout.str(); 
                lowestTrig.append(ss);
                if (trig==lowestTrig) lowestMuonUnprescaledTrig = trig ;
                if (trig==lowestTrig) lowestMuonUnprescaledTrigFound = true ;
                if (trig==lowestTrig) break ;
              } 
              if (lowestMuonUnprescaledTrigFound) break; 
             
              lowestTrig = lowestTrigv0;
              if (trig==lowestTrig) lowestMuonUnprescaledTrig = trig ;
              //      if (trig==lowestTrig) {std::cout << " before break, lowestTrig lowest single muon trigger present unprescaled: " << lowestTrig << std::endl; }
              if (trig==lowestTrig) lowestMuonUnprescaledTrigFound = true ;
              if (trig==lowestTrig) break ;
            }
            if (lowestMuonUnprescaledTrigFound) break; 
          }
        }
      }
      //  std::cout << "after break, lowest single muon trigger present unprescaled: " << lowestMuonUnprescaledTrig << std::endl; 
      unsigned int triggerIndex; // index of trigger path
      
      // See if event passed trigger paths
      std::string hltPath_ = lowestMuonUnprescaledTrig;
 
      triggerIndex = hltConfigProvider_.triggerIndex(hltPath_);
      if (triggerIndex < triggerResults->size()) trigger_fired = triggerResults->accept(triggerIndex);
      std::string L3FilterName_="";
      if (trigger_fired){
        const std::vector<std::string>& moduleLabs = hltConfigProvider_.moduleLabels(hltPath_); 
      /*for (size_t k =0; k < moduleLabs.size()-1 ; k++){
        std::cout << "moduleLabs[" << k << "] == " << moduleLabs[k] << std::endl;
      }
      */
      // the l3 filter name is just the last module.... 
        size_t moduleLabsSizeMinus2 = moduleLabs.size() - 2 ;
        //      std::cout<<"moduleLabs[" << moduleLabsSizeMinus2 << "]== "<< moduleLabs[moduleLabsSizeMinus2] << std::endl;        
         
        L3FilterName_ = moduleLabs[moduleLabsSizeMinus2];
       
      }

      
      edm::Handle< trigger::TriggerEvent > handleTriggerEvent;
      LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << hltPath_ << ")";
      if ( ! ev.getByLabel( trigEv_, handleTriggerEvent ))  {
        //LogWarning( "errorTriggerEventValid" ) << "trigger::TriggerEvent product with InputTag " << trigEv_.encode() << " not in event";
        return;
      }
      ev.getByLabel( trigEv_, handleTriggerEvent );
      const trigger::TriggerObjectCollection & toc(handleTriggerEvent->getObjects());
      size_t nMuHLT =0;
      std::vector<reco::Particle>  HLTMuMatched; 
      for ( size_t ia = 0; ia < handleTriggerEvent->sizeFilters(); ++ ia) {
        std::string fullname = handleTriggerEvent->filterTag(ia).encode();
        //std::cout<< "fullname::== " << fullname<< std::endl;
        std::string name;
        size_t p = fullname.find_first_of(':');
        if ( p != std::string::npos) {
          name = fullname.substr(0, p);
        }
        else {
          name = fullname;
        }
        if ( &toc !=0 ) {
          const trigger::Keys & k = handleTriggerEvent->filterKeys(ia);
          for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
            // looking at all the single muon l3 trigger present, for example hltSingleMu15L3Filtered15.....
            if (name == L3FilterName_  ) {
              // trigger_fired = true;
              hlt_sel=true;
              nhlt++; 
              HLTMuMatched.push_back(toc[*ki].particle());
              nMuHLT++;     
            }
          }    
        }
      }
      
      // Beam spot
      Handle<reco::BeamSpot> beamSpotHandle;
      if (!ev.getByLabel(InputTag("offlineBeamSpot"), beamSpotHandle)) {
        //LogWarning("") << ">>> No beam spot found !!!";
        return;
      }

      
        //  looping on muon....
      Handle<View<Muon> >   muons;     
      if (!ev.getByLabel(muonTag_, muons)) {           
        //LogError("") << ">>> muon collection does not exist !!!";     
        return;     
      }
      
      ev.getByLabel(muonTag_, muons);  
      //saving only muons with pt> ptMuCut and eta<etaMuCut, and dxy<dxyCut  
      std::vector<reco::Muon>  highPtGlbMuons; 
      std::vector<reco::Muon>  highPtStaMuons; 

      for (size_t i=0; i<muons->size(); i++ ){  
        const reco::Muon & mu = muons->at(i);
        double pt = mu.pt();
        double eta = mu.eta();
        if (pt> ptMuCut_ && fabs(eta)< etaMuCut_ ) {
          if (mu.isGlobalMuon()) { 
            //check the dxy....
            double dxy = mu.innerTrack()->dxy(beamSpotHandle->position()); 
            if (fabs(dxy)>dxyCut_) continue;       
            highPtGlbMuons.push_back(mu);
          }
          if (mu.isGlobalMuon()) continue;
          // if is not, look if it is a standalone.... 
          if(mu.isStandAloneMuon()) highPtStaMuons.push_back(mu); 
          nEvWithHighPtMu++;     
        }
      }
      size_t nHighPtGlbMu = highPtGlbMuons.size();
      size_t nHighPtStaMu = highPtStaMuons.size();
      if ( hlt_sel &&  (nHighPtGlbMu> 0) ) {
        // loop on high pt muons if there's at least two with opposite charge build a Z, more then one z candidate is foreseen.........
        // stop the loop after 10 cicles....  
        (nHighPtGlbMu> 10)?   nHighPtGlbMu=10 : 1; 
        // Z selection
        if (nHighPtGlbMu>1 ){
          for(unsigned int i =0 ; i < nHighPtGlbMu ; i++) {
            reco::Muon muon1 = highPtGlbMuons[i];
            math::XYZTLorentzVector mu1(muon1.p4());
            //double pt1= muon1.pt();
            for (unsigned int j =i+1; j <nHighPtGlbMu ; ++j ){
              reco::Muon muon2 = highPtGlbMuons[j];
              math::XYZTLorentzVector mu2(muon2.p4());
              //double pt2= muon2.pt();
              if (muon1.charge() == muon2.charge()) continue; 
              math::XYZTLorentzVector pair = mu1 + mu2;
              double mass = pair.M();
                // checking isolation and hlt maching
              iso1 = muIso ( muon1 );
              iso2 = muIso ( muon2 );
              isMu1Iso = (iso1 < isoCut03_);
              isMu2Iso = (iso2 < isoCut03_);
              singleTrigFlag1 = IsMuMatchedToHLTMu ( muon1,  HLTMuMatched ,maxDeltaR_, maxDPtRel_ );
              singleTrigFlag2 = IsMuMatchedToHLTMu ( muon2,  HLTMuMatched ,maxDeltaR_, maxDPtRel_ );
              if (singleTrigFlag1 && singleTrigFlag2)  isZGolden2HLT_ = true; 
              if ( (singleTrigFlag1 && !singleTrigFlag2) || (!singleTrigFlag1 && singleTrigFlag2) )  isZGolden1HLT_ = true;  
              // Z Golden passing all criteria, with 2 or 1 muon matched to an HLT muon. Using the two cathegories as a control sample for the HLT matching efficiency 
              if (isZGolden2HLT_ || isZGolden1HLT_){
                // a Z golden has been found, let's remove the two muons from the list, dome for avoiding resolution effect enter muons in the standalone and tracker collections.........
                highPtGlbMuons.erase(highPtGlbMuons.begin() + i );
                highPtGlbMuons.erase(highPtGlbMuons.begin() + j );
                // and updating the number of high pt muons....
                nHighPtGlbMu = highPtGlbMuons.size();

                if (  isMu1Iso &&  isMu2Iso  ){
                  niso++;
                  iso_sel=true;
                  if (isZGolden2HLT_) {
                    n2hlt++; 
                    mass2HLT_->Fill(mass);
                    highMass2HLT_->Fill(mass);
                    /*
                      if (pt1 > pt2) {
                      highest_mupt2HLT_ -> Fill (pt1);
                      lowest_mupt2HLT_ -> Fill (pt2);
                      } else {
                      highest_mupt2HLT_ -> Fill (pt2);
                      lowest_mupt2HLT_ -> Fill (pt1);
                      } 
                    */  
                  }
                  if (isZGolden1HLT_) {
                    n1hlt++; 
                    mass1HLT_->Fill(mass);
                    highMass1HLT_->Fill(mass);
                    /*
                      if (pt1 >pt2) {
                      highest_mupt1HLT_ -> Fill (pt1);
                      lowest_mupt1HLT_ -> Fill (pt2);
                    } else {
                      highest_mupt1HLT_ -> Fill (pt2);
                      lowest_mupt1HLT_ -> Fill (pt1);
                    }
                    */
                  } 
                } else {
                  // ZGlbGlb when at least one of the two muon is not isolated and at least one HLT matched, used as control sample for the isolation efficiency
                  // filling events with one muon not isolated and both hlt mathced  
                  if (((isMu1Iso &&  !isMu2Iso) || (!isMu1Iso &&  isMu2Iso )) && (isZGolden2HLT_ && isZGolden1HLT_)  ) {
                    isZGoldenNoIso_=true;
                    nNotIso++;
                    massNotIso_->Fill(mass);
                    highMassNotIso_->Fill(mass);
                  }
                  /*
                    if (pt1 > pt2) {
                    highest_muptNotIso_ -> Fill (pt1);
                    lowest_muptNotIso_ -> Fill (pt2);
                    } else {
                    highest_muptNotIso_ -> Fill (pt2);
                    lowest_muptNotIso_ -> Fill (pt1);
                    } 
                  */ 
                }
              }
            }
          }
        }
// W->MuNu selection (if at least one muon has been selected)
        // looking for a W if a Z is not found.... let's think if we prefer to exclude zMuMuNotIso or zMuSta....
        if ( !(isZGolden2HLT_ || isZGolden1HLT_ )){
          Handle<View<MET> > metCollection;     
          if (!ev.getByLabel(metTag_, metCollection)) {           
            //LogError("") << ">>> MET collection does not exist !!!";           
            return;     
          }     
          const MET& met = metCollection->at(0);   
          nW=0;
          for(unsigned int i =0 ; i < nHighPtGlbMu ; i++) {
            reco::Muon muon1 = highPtGlbMuons[i];
            /*      
                   quality cut not implemented
                Quality Cuts           double dxy = gm->dxy(beamSpotHandle->position()); 
                Cut in 0.2           double trackerHits = gm->hitPattern().numberOfValidTrackerHits(); 
                Cut in 11           bool quality = fabs(dxy)<dxyCut_  && muon::isGoodMuon(mu,muon::GlobalMuonPromptTight) && trackerHits>=trackerHitsCut_ && mu.isTrackerMuon() ;  
  if (!quality) continue;          
            */
            // isolation cut and hlt maching
            iso1 = muIso ( muon1 );
            isMu1Iso = (iso1 < isoCut03_);
            if (!isMu1Iso) continue;
            // checking if muon is matched to any HLT muon....
            singleTrigFlag1 = IsMuMatchedToHLTMu ( muon1,  HLTMuMatched ,maxDeltaR_, maxDPtRel_ );
            if (!singleTrigFlag1) continue;
            //std::cout << " is GlobMu macthecd to HLT: " << singleTrigFlag1 << std::endl; 
            // MT cuts        
            double met_px = met.px();
            double met_py = met.py();
            
            if (!metIncludesMuons_) {
              for(unsigned int i =0 ; i < nHighPtGlbMu ; i++) {
                reco::Muon muon1 = highPtGlbMuons[i];
                met_px -= muon1.px();
                met_py -= muon1.py();
              }
            }
            double met_et = met.pt() ; //sqrt(met_px*met_px+met_py*met_py);
            LogTrace("") << ">>> MET, MET_px, MET_py: " << met_et << ", " << met_px << ", " << met_py << " [GeV]";
            double w_et = met_et+ muon1.pt();           
            double w_px = met_px+ muon1.px();           
            double w_py = met_py+muon1.py();           
            double massT = w_et*w_et - w_px*w_px - w_py*w_py;           
            massT = (massT>0) ? sqrt(massT) : 0;
            // Acoplanarity cuts           
            Geom::Phi<double> deltaphi(muon1.phi()-atan2(met_py,met_px));           
            double acop = M_PI - fabs(deltaphi.value());  
            //std::cout << " is acp of W candidate less then cut: " << (acop< acopCut_) << std::endl; 
            if (acop > acopCut_) continue ; // Cut in 2.0    
            TMass_->Fill(massT);
            if (massT<mtMin_ || massT>mtMax_) continue;    // Cut in (50,200) GeV           
            // we give the number of W only in the tMass selected but we have a look at mass tails to check the QCD background 
            isW_=true;
            nW++;
            nTMass++;     
          }
        }
// if a zGlobal is not selected, look at the dimuonGlobalOneStandAlone and dimuonGlobalOneTrack...., used as a control sample for the tracking efficiency  and muon system efficiency
        if ( !(isZGolden2HLT_ || isZGolden1HLT_ || isZGoldenNoIso_ )){
          for(unsigned int i =0 ; i < nHighPtGlbMu ; ++i) {
            reco::Muon glbMuon = highPtGlbMuons[i];
            math::XYZTLorentzVector mu1(glbMuon.p4());
            //double pt1= glbMuon.pt();
            // checking that the global muon is hlt matched otherwise skip the event
            singleTrigFlag1 = IsMuMatchedToHLTMu ( glbMuon,  HLTMuMatched ,maxDeltaR_, maxDPtRel_ );
            if (!singleTrigFlag1) continue;
            // checking that the global muon is isolated matched otherwise skip the event
            iso1 = muIso ( glbMuon );
            isMu1Iso = (iso1 < isoCut03_);
            if (!isMu1Iso) continue;
            // look at the standalone muon ....
            // stop the loop after 10 cicles....  
            (nHighPtStaMu> 10)?   nHighPtStaMu=10 : 1; 
            for (unsigned int j =0; j <nHighPtStaMu ; ++j ){
              reco::Muon staMuon = highPtStaMuons[j];
              math::XYZTLorentzVector mu2(staMuon.p4());
              //double pt2= staMuon.pt();
              if (glbMuon.charge() == staMuon.charge()) continue; 
              math::XYZTLorentzVector pair = mu1 + mu2;
              double mass = pair.M();
              iso2 = muIso ( staMuon);
              isMu2Iso = (iso2 < isoCut03_);
              LogTrace("") << "\t... isolation value" << iso1 <<", isolated? " << isMu1Iso ;
              LogTrace("") << "\t... isolation value" << iso2 <<", isolated? " << isMu2Iso ;
              // requiring theat the global mu is mathed to the HLT  and both are isolated 
              if (isMu2Iso ) {    
                isZGlbSta_=true;
                nGlbSta++; 
                massGlbSta_->Fill(mass);
                highMassGlbSta_->Fill(mass);
                /*
                  if (pt1 > pt2) {
                  highest_muptGlbSta_ -> Fill (pt1);
                  lowest_muptGlbSta_ -> Fill (pt2);
                } else {
                  highest_muptGlbSta_ -> Fill (pt2);
                  lowest_muptGlbSta_ -> Fill (pt1);
                }
                */
              }
            }
            // look at the tracks....
            Handle< TrackCollection >   tracks;     
            if (!ev.getByLabel(trackTag_, tracks)) {           
              //LogError("") << ">>> track collection does not exist !!!";     
              return;     
            }
            ev.getByLabel(trackTag_, tracks);        
            Handle< CaloTowerCollection >   calotower;     
            if (!ev.getByLabel(caloTowerTag_, calotower)) {           
              //LogError("") << ">>> calotower collection does not exist !!!";     
              return;     
            } 
            ev.getByLabel(caloTowerTag_, calotower);
            // avoid to loop on more than 5000 trks
            size_t nTrk = tracks->size();
            (nTrk> 5000)?   nTrk=5000 : 1; 
            for (unsigned int j=0; j<nTrk; j++ ){       
              const reco::Track & tk = tracks->at(j);
              if (glbMuon.charge() == tk.charge()) continue; 
              double pt2 = tk.pt();
              double eta = tk.eta();
              double dxy = tk.dxy(beamSpotHandle->position());
              if (pt2< ptMuCut_ || fabs(eta) > etaMuCut_ || fabs(dxy)>dxyCut_) continue;
              //assuming that the track is a mu....
              math::XYZTLorentzVector mu2(tk.px(),tk.py(), tk.pz(), sqrt( (tk.p() * tk.p())  + ( 0.10566 * 0.10566))) ;
              math::XYZTLorentzVector pair = mu1 + mu2;
              double mass = pair.M();
              // now requiring that the tracks is isolated....... 
              iso2 = tkIso( tk, tracks, calotower);
              isMu2Iso = (iso2 < isoCut03_);
              //              std::cout << "found a track with rel/comb iso: " << iso2 << std::endl; 
              if (isMu2Iso) {    
                isZGlbTrk_=true;
                nGlbTrk++; 
                massGlbTrk_->Fill(mass);
                highMassGlbTrk_->Fill(mass);
                if (!isW_) continue;
                massIsBothGlbTrkThanW_->Fill(mass);
                highMassIsBothGlbTrkThanW_->Fill(mass);
                /*
                  if (pt1 > pt2) {
                  highest_muptGlbTrk_ -> Fill (pt1);
                  lowest_muptGlbTrk_ -> Fill (pt2);
                  } else {
                  highest_muptGlbTrk_ -> Fill (pt2);
                  lowest_muptGlbTrk_ -> Fill (pt1);
                } 
                */
              }
            }
          }
        }
      }
      if ( (hlt_sel || isZGolden2HLT_ || isZGolden1HLT_ || isZGoldenNoIso_ || isZGlbSta_ || isZGlbTrk_ || isW_) ) {   
        nsel++;  
        LogTrace("") << ">>>> Event ACCEPTED";
      } else {
        LogTrace("") << ">>>> Event REJECTED";
      }
      return;
}