/ElectroWeakAnalysis/ZMuMu/plugins/ZMuMu_radiative_analysis.cc
C++ | 558 lines | 476 code | 59 blank | 23 comment | 79 complexity | cd7c36f24c1efd326ae42582bcc49626 MD5 | raw file
1/* \class ZMuMu_Radiative_analyzer 2 * 3 * author: Pasquale Noli 4 * 5 * ZMuMu Radiative analyzer 6 * 7 * 8 */ 9#include "DataFormats/HepMCCandidate/interface/GenParticle.h" 10#include "DataFormats/RecoCandidate/interface/RecoCandidate.h" 11#include "DataFormats/MuonReco/interface/Muon.h" 12#include "DataFormats/Common/interface/AssociationVector.h" 13#include "DataFormats/Candidate/interface/CandMatchMap.h" 14#include "DataFormats/Candidate/interface/Particle.h" 15#include "DataFormats/Candidate/interface/Candidate.h" 16#include "DataFormats/Candidate/interface/CandidateFwd.h" 17#include "DataFormats/Candidate/interface/OverlapChecker.h" 18#include "DataFormats/PatCandidates/interface/Muon.h" 19#include "DataFormats/PatCandidates/interface/Lepton.h" 20#include "DataFormats/PatCandidates/interface/GenericParticle.h" 21#include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h" 22#include "DataFormats/PatCandidates/interface/PATObject.h" 23#include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h" 24#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h" 25#include "DataFormats/RecoCandidate/interface/IsoDeposit.h" 26#include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h" 27#include "DataFormats/PatCandidates/interface/Isolation.h" 28#include "DataFormats/Common/interface/Handle.h" 29#include "DataFormats/Common/interface/ValueMap.h" 30#include "DataFormats/Candidate/interface/CandAssociation.h" 31#include "FWCore/Framework/interface/EDAnalyzer.h" 32#include "FWCore/ParameterSet/interface/ParameterSet.h" 33#include "FWCore/Framework/interface/Event.h" 34#include "FWCore/Framework/interface/EventSetup.h" 35#include "FWCore/Utilities/interface/InputTag.h" 36#include "FWCore/ServiceRegistry/interface/Service.h" 37#include "CommonTools/UtilAlgos/interface/TFileService.h" 38#include <iostream> 39#include <iterator> 40#include <cmath> 41#include <vector> 42#include "TH1.h" 43#include "TH2.h" 44#include "TH3.h" 45 46 47using namespace edm; 48using namespace std; 49using namespace reco; 50using namespace isodeposit; 51 52class ZMuMu_Radiative_analyzer : public edm::EDAnalyzer { 53public: 54 ZMuMu_Radiative_analyzer(const edm::ParameterSet& pset); 55private: 56 virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override; 57 virtual void endJob() override; 58 59 edm::InputTag zMuMu_, zMuMuMatchMap_, zMuTk_, zMuTkMatchMap_,zMuSa_, zMuSaMatchMap_; 60 double dRVeto_, dRTrk_, ptThreshold_; 61 //histograms 62 TH1D *h_zmass_FSR,*h_zmass_no_FSR; 63 TH1D *h_zMuSamass_FSR,*h_zMuSamass_no_FSR; 64 TH1D *h_zMuTkmass_FSR,*h_zMuTkmass_no_FSR; 65 TH1D *h_Iso_,*h_Iso_FSR_ ; 66 TH3D *h_Iso_3D_, *h_Iso_FSR_3D_; 67 TH2D *h_staProbe_pt_eta_no_FSR_, *h_staProbe_pt_eta_FSR_, *h_ProbeOk_pt_eta_no_FSR_, *h_ProbeOk_pt_eta_FSR_; 68 TH1D *h_trackProbe_eta_no_FSR, *h_trackProbe_pt_no_FSR, *h_staProbe_eta_no_FSR, *h_staProbe_pt_no_FSR, *h_ProbeOk_eta_no_FSR, *h_ProbeOk_pt_no_FSR; 69 TH1D *h_trackProbe_eta_FSR, *h_trackProbe_pt_FSR, *h_staProbe_eta_FSR, *h_staProbe_pt_FSR, *h_ProbeOk_eta_FSR, *h_ProbeOk_pt_FSR; 70 //boolean 71 bool FSR_mu, FSR_tk, FSR_mu0, FSR_mu1; 72 bool trig0found, trig1found; 73 //counter 74 int zmmcounter , zmscounter, zmtcounter, evntcounter; 75 }; 76 77 78typedef edm::ValueMap<float> IsolationCollection; 79 80ZMuMu_Radiative_analyzer::ZMuMu_Radiative_analyzer(const ParameterSet& pset) : 81 zMuMu_(pset.getParameter<InputTag>("zMuMu")), 82 zMuMuMatchMap_(pset.getParameter<InputTag>("zMuMuMatchMap")), 83 zMuTk_(pset.getParameter<InputTag>("zMuTk")), 84 zMuTkMatchMap_(pset.getParameter<InputTag>("zMuTkMatchMap")), 85 zMuSa_(pset.getParameter<InputTag>("zMuSa")), 86 zMuSaMatchMap_(pset.getParameter<InputTag>("zMuSaMatchMap")), 87 dRVeto_(pset.getUntrackedParameter<double>("veto")), 88 dRTrk_(pset.getUntrackedParameter<double>("deltaRTrk")), 89 ptThreshold_(pset.getUntrackedParameter<double>("ptThreshold")){ 90 zmmcounter=0; 91 zmscounter=0; 92 zmtcounter=0; 93 evntcounter=0; 94 Service<TFileService> fs; 95 96 // general histograms 97 h_zmass_FSR= fs->make<TH1D>("h_zmass_FRS","Invariant Z mass distribution",200,0,200); 98 h_zmass_no_FSR= fs->make<TH1D>("h_zmass_no_FSR","Invariant Z mass distribution",200,0,200); 99 h_zMuSamass_FSR= fs->make<TH1D>("h_zMuSamass_FRS","Invariant Z mass distribution",200,0,200); 100 h_zMuSamass_no_FSR= fs->make<TH1D>("h_zMuSamass_no_FSR","Invariant Z mass distribution",200,0,200); 101 h_zMuTkmass_FSR= fs->make<TH1D>("h_zMuTkmass_FRS","Invariant Z mass distribution",200,0,200); 102 h_zMuTkmass_no_FSR= fs->make<TH1D>("h_zMuTkmass_no_FSR","Invariant Z mass distribution",200,0,200); 103 h_Iso_= fs->make<TH1D>("h_iso","Isolation distribution of muons without FSR",100,0,20); 104 h_Iso_FSR_= fs->make<TH1D>("h_iso_FSR","Isolation distribution of muons with FSR ",100,0,20); 105 h_Iso_3D_= fs->make<TH3D>("h_iso_3D","Isolation distribution of muons without FSR",100,20,100,100,-2.0,2.0,100,0,20); 106 h_Iso_FSR_3D_= fs->make<TH3D>("h_iso_FSR_3D","Isolation distribution of muons with FSR ",100,20,100,100,-2.0,2.0,100,0,20); 107 h_staProbe_pt_eta_no_FSR_= fs->make<TH2D>("h_staProbe_pt_eta_no_FSR","Pt vs Eta StandAlone without FSR ",100,20,100,100,-2.0,2.0); 108 h_staProbe_pt_eta_FSR_= fs->make<TH2D>("h_staProbe_pt_eta_FSR","Pt vs Eta StandAlone with FSR ",100,20,100,100,-2.0,2.0); 109 h_ProbeOk_pt_eta_no_FSR_= fs->make<TH2D>("h_ProbeOk_pt_eta_no_FSR","Pt vs Eta probeOk without FSR ",100,20,100,100,-2.0,2.0); 110 h_ProbeOk_pt_eta_FSR_= fs->make<TH2D>("h_ProbeOk_pt_eta_FSR","Pt vs Eta probeOk with FSR ",100,20,100,100,-2.0,2.0); 111 112 h_trackProbe_eta_no_FSR = fs->make<TH1D>("trackProbeEta_no_FSR","Eta of tracks",100,-2.0,2.0); 113 h_trackProbe_pt_no_FSR = fs->make<TH1D>("trackProbePt_no_FSR","Pt of tracks",100,0,100); 114 h_staProbe_eta_no_FSR = fs->make<TH1D>("standAloneProbeEta_no_FSR","Eta of standAlone",100,-2.0,2.0); 115 h_staProbe_pt_no_FSR = fs->make<TH1D>("standAloneProbePt_no_FSR","Pt of standAlone",100,0,100); 116 h_ProbeOk_eta_no_FSR = fs->make<TH1D>("probeOkEta_no_FSR","Eta of probe Ok",100,-2.0,2.0); 117 h_ProbeOk_pt_no_FSR = fs->make<TH1D>("probeOkPt_no_FSR","Pt of probe ok",100,0,100); 118 119 h_trackProbe_eta_FSR = fs->make<TH1D>("trackProbeEta_FSR","Eta of tracks",100,-2.0,2.0); 120 h_trackProbe_pt_FSR = fs->make<TH1D>("trackProbePt_FSR","Pt of tracks",100,0,100); 121 h_staProbe_eta_FSR = fs->make<TH1D>("standAloneProbeEta_FSR","Eta of standAlone",100,-2.0,2.0); 122 h_staProbe_pt_FSR = fs->make<TH1D>("standAloneProbePt_FSR","Pt of standAlone",100,0,100); 123 h_ProbeOk_eta_FSR = fs->make<TH1D>("probeOkEta_FSR","Eta of probe Ok",100,-2.0,2.0); 124 h_ProbeOk_pt_FSR = fs->make<TH1D>("probeOkPt_FSR","Pt of probe ok",100,0,100); 125 126 127 128} 129 130void ZMuMu_Radiative_analyzer::analyze(const Event& event, const EventSetup& setup) { 131 evntcounter++; 132 Handle<CandidateView> zMuMu; //Collection of Z made by Mu global + Mu global 133 Handle<GenParticleMatch> zMuMuMatchMap; //Map of Z made by Mu global + Mu global with MC 134 event.getByLabel(zMuMu_, zMuMu); 135 Handle<CandidateView> zMuTk; //Collection of Z made by Mu global + Track 136 Handle<GenParticleMatch> zMuTkMatchMap; 137 event.getByLabel(zMuTk_, zMuTk); 138 Handle<CandidateView> zMuSa; //Collection of Z made by Mu global + Sa 139 Handle<GenParticleMatch> zMuSaMatchMap; 140 event.getByLabel(zMuSa_, zMuSa); 141 cout << "********** New Event ***********"<<endl; 142 // ZMuMu 143 if (zMuMu->size() > 0 ) { 144 event.getByLabel(zMuMuMatchMap_, zMuMuMatchMap); 145 for(unsigned int i = 0; i < zMuMu->size(); ++i) { //loop on candidates 146 147 const Candidate & zMuMuCand = (*zMuMu)[i]; //the candidate 148 CandidateBaseRef zMuMuCandRef = zMuMu->refAt(i); 149 150 151 CandidateBaseRef dau0 = zMuMuCand.daughter(0)->masterClone(); 152 CandidateBaseRef dau1 = zMuMuCand.daughter(1)->masterClone(); 153 const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon 154 const pat::Muon& mu1 = dynamic_cast<const pat::Muon&>(*dau1); 155 156 double zmass= zMuMuCand.mass(); 157 double pt0 = mu0.pt(); 158 double pt1 = mu1.pt(); 159 double eta0 = mu0.eta(); 160 double eta1 = mu1.eta(); 161 if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){ 162 GenParticleRef zMuMuMatch = (*zMuMuMatchMap)[zMuMuCandRef]; 163 if(zMuMuMatch.isNonnull()) { // ZMuMu matched 164 zmmcounter++; 165 cout<<" Zmumu cuts && matched" <<endl; 166 FSR_mu0 = false; 167 FSR_mu1 = false; 168 169 //Isodeposit 170 const pat::IsoDeposit * mu0TrackIso =mu0.isoDeposit(pat::TrackIso); 171 const pat::IsoDeposit * mu1TrackIso =mu1.isoDeposit(pat::TrackIso); 172 Direction mu0Dir = Direction(mu0.eta(),mu0.phi()); 173 Direction mu1Dir = Direction(mu1.eta(),mu1.phi()); 174 175 reco::IsoDeposit::AbsVetos vetos_mu0; 176 vetos_mu0.push_back(new ConeVeto( mu0Dir, dRVeto_ )); 177 vetos_mu0.push_back(new ThresholdVeto( ptThreshold_ )); 178 179 reco::IsoDeposit::AbsVetos vetos_mu1; 180 vetos_mu1.push_back(new ConeVeto( mu1Dir, dRVeto_ )); 181 vetos_mu1.push_back(new ThresholdVeto( ptThreshold_ )); 182 183 double Tracker_isovalue_mu0 = mu0TrackIso->sumWithin(dRTrk_,vetos_mu0); 184 double Tracker_isovalue_mu1 = mu1TrackIso->sumWithin(dRTrk_,vetos_mu1); 185 186 //trigger study 187 trig0found = false; 188 trig1found = false; 189 const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 190 mu0.triggerObjectMatchesByPath( "HLT_Mu9" ); 191 const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 192 mu1.triggerObjectMatchesByPath( "HLT_Mu9" ); 193 if( mu0HLTMatches.size()>0 ) 194 trig0found = true; 195 if( mu1HLTMatches.size()>0 ) 196 trig1found = true; 197 198 //MonteCarlo Study 199 const reco::GenParticle * muMc0 = mu0.genLepton(); 200 const reco::GenParticle * muMc1 = mu1.genLepton(); 201 const Candidate * motherMu0 = muMc0->mother(); 202 int num_dau_muon0 = motherMu0->numberOfDaughters(); 203 const Candidate * motherMu1 = muMc1->mother(); 204 int num_dau_muon1 = motherMu1->numberOfDaughters(); 205 cout<<" muone0"<<endl; 206 cout<<" num di daughters = "<< num_dau_muon0 <<endl; 207 if( num_dau_muon0 > 1 ){ 208 for(int j = 0; j < num_dau_muon0; ++j){ 209 int id =motherMu0 ->daughter(j)->pdgId(); 210 cout<<" dauther["<<j<<"] pdgId = "<<id<<endl; 211 if(id == 22) FSR_mu0=true; 212 } 213 }//end check of gamma 214 215 cout<<" muone1"<<endl; 216 cout<<" num di daughters = "<< num_dau_muon1 <<endl; 217 if( num_dau_muon1 > 1 ){ 218 for(int j = 0; j < num_dau_muon1; ++j){ 219 int id = motherMu1->daughter(j)->pdgId(); 220 cout<<" dauther["<<j<<"] pdgId = "<<id<<endl; 221 if(id == 22) FSR_mu1=true; 222 } 223 }//end check of gamma 224 225 if(FSR_mu0 || FSR_mu1 )h_zmass_FSR->Fill(zmass); 226 else h_zmass_no_FSR->Fill(zmass); 227 228 if (trig1found) { // check efficiency of muon0 not imposing the trigger on it 229 cout<<"muon 1 is triggered "<<endl; 230 if(FSR_mu0){ 231 cout<< "and muon 0 does FSR"<<endl; 232 h_trackProbe_eta_FSR->Fill(eta0); 233 h_trackProbe_pt_FSR->Fill(pt0); 234 h_staProbe_eta_FSR->Fill(eta0); 235 h_staProbe_pt_FSR->Fill(pt0); 236 h_staProbe_pt_eta_FSR_->Fill(pt0,eta0); 237 h_ProbeOk_eta_FSR->Fill(eta0); 238 h_ProbeOk_pt_FSR->Fill(pt0); 239 h_ProbeOk_pt_eta_FSR_->Fill(pt0,eta0); 240 }else{ 241 cout<<"and muon 0 doesn't FSR"<<endl; 242 h_trackProbe_eta_no_FSR->Fill(eta0); 243 h_trackProbe_pt_no_FSR->Fill(pt0); 244 h_staProbe_eta_no_FSR->Fill(eta0); 245 h_staProbe_pt_no_FSR->Fill(pt0); 246 h_staProbe_pt_eta_no_FSR_->Fill(pt0,eta0); 247 h_ProbeOk_eta_no_FSR->Fill(eta0); 248 h_ProbeOk_pt_no_FSR->Fill(pt0); 249 h_ProbeOk_pt_eta_no_FSR_->Fill(pt0,eta0); 250 } 251 if(FSR_mu0){ 252 h_Iso_FSR_->Fill(Tracker_isovalue_mu0); 253 h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0); 254 } 255 else{ 256 h_Iso_->Fill(Tracker_isovalue_mu0); 257 h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0); 258 } 259 } 260 if (trig0found) { // check efficiency of muon1 not imposing the trigger on it 261 cout<<"muon 0 is triggered"<<endl; 262 if(FSR_mu1){ 263 cout<<"and muon 1 does FSR"<<endl; 264 h_trackProbe_eta_FSR->Fill(eta1); 265 h_staProbe_eta_FSR->Fill(eta1); 266 h_trackProbe_pt_FSR->Fill(pt1); 267 h_staProbe_pt_FSR->Fill(pt1); 268 h_ProbeOk_eta_FSR->Fill(eta1); 269 h_ProbeOk_pt_FSR->Fill(pt1); 270 h_staProbe_pt_eta_FSR_->Fill(pt1,eta1); 271 h_ProbeOk_pt_eta_FSR_->Fill(pt1,eta1); 272 273 }else{ 274 cout<<"and muon 1 doesn't FSR"<<endl; 275 h_trackProbe_eta_no_FSR->Fill(eta1); 276 h_staProbe_eta_no_FSR->Fill(eta1); 277 h_trackProbe_pt_no_FSR->Fill(pt1); 278 h_staProbe_pt_no_FSR->Fill(pt1); 279 h_ProbeOk_eta_no_FSR->Fill(eta1); 280 h_ProbeOk_pt_no_FSR->Fill(pt1); 281 h_staProbe_pt_eta_no_FSR_->Fill(pt1,eta1); 282 h_ProbeOk_pt_eta_no_FSR_->Fill(pt1,eta1); 283 284 285 } 286 if(FSR_mu1){ 287 h_Iso_FSR_->Fill(Tracker_isovalue_mu1); 288 h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1); 289 }else{ 290 h_Iso_->Fill(Tracker_isovalue_mu1); 291 h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1); 292 } 293 } 294 }// end MC match 295 }//end of cuts 296 }// end loop on ZMuMu cand 297 }// end if ZMuMu size > 0 298 299 // ZMuSa 300 if (zMuSa->size() > 0 ) { 301 event.getByLabel(zMuSaMatchMap_, zMuSaMatchMap); 302 for(unsigned int i = 0; i < zMuSa->size(); ++i) { //loop on candidates 303 304 const Candidate & zMuSaCand = (*zMuSa)[i]; //the candidate 305 CandidateBaseRef zMuSaCandRef = zMuSa->refAt(i); 306 const Candidate * lep0 =zMuSaCand.daughter(0); 307 const Candidate * lep1 =zMuSaCand.daughter(1); 308 CandidateBaseRef dau0 = lep0->masterClone(); 309 CandidateBaseRef dau1 = lep1->masterClone(); 310 const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon 311 const pat::Muon& mu1 = dynamic_cast<const pat::Muon&>(*dau1); 312 313 double zmass= zMuSaCand.mass(); 314 double pt0 = mu0.pt(); 315 double pt1 = mu1.pt(); 316 double eta0 = mu0.eta(); 317 double eta1 = mu1.eta(); 318 if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){ 319 GenParticleRef zMuSaMatch = (*zMuSaMatchMap)[zMuSaCandRef]; 320 if(zMuSaMatch.isNonnull()) { // ZMuSa matched 321 cout<<" Zmusa cuts && matched" <<endl; 322 FSR_mu0 = false; 323 FSR_mu1 = false; 324 zmscounter++; 325 //Isodeposit 326 const pat::IsoDeposit * mu0TrackIso =mu0.isoDeposit(pat::TrackIso); 327 const pat::IsoDeposit * mu1TrackIso =mu1.isoDeposit(pat::TrackIso); 328 Direction mu0Dir = Direction(mu0.eta(),mu0.phi()); 329 Direction mu1Dir = Direction(mu1.eta(),mu1.phi()); 330 331 reco::IsoDeposit::AbsVetos vetos_mu0; 332 vetos_mu0.push_back(new ConeVeto( mu0Dir, dRVeto_ )); 333 vetos_mu0.push_back(new ThresholdVeto( ptThreshold_ )); 334 335 reco::IsoDeposit::AbsVetos vetos_mu1; 336 vetos_mu1.push_back(new ConeVeto( mu1Dir, dRVeto_ )); 337 vetos_mu1.push_back(new ThresholdVeto( ptThreshold_ )); 338 339 double Tracker_isovalue_mu0 = mu0TrackIso->sumWithin(dRTrk_,vetos_mu0); 340 double Tracker_isovalue_mu1 = mu1TrackIso->sumWithin(dRTrk_,vetos_mu1); 341 342 // HLT match (check just dau0 the global) 343 const pat::TriggerObjectStandAloneCollection mu0HLTMatches = 344 mu0.triggerObjectMatchesByPath( "HLT_Mu9" ); 345 const pat::TriggerObjectStandAloneCollection mu1HLTMatches = 346 mu1.triggerObjectMatchesByPath( "HLT_Mu9" ); 347 trig0found = false; 348 trig1found = false; 349 if( mu0HLTMatches.size()>0 ) 350 trig0found = true; 351 if( mu1HLTMatches.size()>0 ) 352 trig1found = true; 353 354 //MonteCarlo Study 355 const reco::GenParticle * muMc0 = mu0.genLepton(); 356 const reco::GenParticle * muMc1 = mu1.genLepton(); 357 const Candidate * motherMu0 = muMc0->mother(); 358 const Candidate * motherMu1 = muMc1->mother(); 359 int num_dau_muon0 = motherMu0->numberOfDaughters(); 360 int num_dau_muon1 = motherMu1->numberOfDaughters(); 361 cout<<" muone0"<<endl; 362 cout<<" num di daughters = "<< num_dau_muon0 <<endl; 363 if( num_dau_muon0 > 1 ){ 364 for(int j = 0; j < num_dau_muon0; ++j){ 365 int id =motherMu0 ->daughter(j)->pdgId(); 366 cout<<" dauther["<<j<<"] pdgId = "<<id<<endl; 367 if(id == 22) FSR_mu0=true; 368 } 369 }//end check of gamma 370 371 cout<<" muone1"<<endl; 372 cout<<" num di daughters = "<< num_dau_muon1 <<endl; 373 if( num_dau_muon1 > 1 ){ 374 for(int j = 0; j < num_dau_muon1; ++j){ 375 int id = motherMu1->daughter(j)->pdgId(); 376 cout<<" dauther["<<j<<"] pdgId = "<<id<<endl; 377 if(id == 22) FSR_mu1=true; 378 } 379 }//end check of gamma 380 if(FSR_mu0 || FSR_mu1 )h_zMuSamass_FSR->Fill(zmass); 381 else h_zMuSamass_no_FSR->Fill(zmass); 382 if(lep0->isGlobalMuon() && trig0found){ 383 if(FSR_mu1){ 384 h_staProbe_eta_FSR->Fill(eta1); 385 h_staProbe_pt_FSR->Fill(pt1); 386 h_staProbe_pt_eta_FSR_->Fill(pt1,eta1); 387 388 }else{ 389 h_staProbe_eta_no_FSR->Fill(eta1); 390 h_staProbe_pt_no_FSR->Fill(pt1); 391 h_staProbe_pt_eta_no_FSR_->Fill(pt1,eta1); 392 393 } 394 if(FSR_mu1){ 395 h_Iso_FSR_->Fill(Tracker_isovalue_mu1); 396 h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1); 397 } 398 else{ 399 h_Iso_->Fill(Tracker_isovalue_mu1); 400 h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_mu1); 401 } 402 } 403 if(lep1->isGlobalMuon() && trig1found){ 404 if(FSR_mu0){ 405 h_staProbe_eta_FSR->Fill(eta0); 406 h_staProbe_pt_FSR->Fill(pt0); 407 h_staProbe_pt_eta_FSR_->Fill(pt0,eta0); 408 409 }else{ 410 h_staProbe_eta_no_FSR->Fill(eta0); 411 h_staProbe_pt_no_FSR->Fill(pt0); 412 h_staProbe_pt_eta_FSR_->Fill(pt0,eta0); 413 414 } 415 if(FSR_mu0){ 416 h_Iso_FSR_->Fill(Tracker_isovalue_mu0); 417 h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0); 418 } 419 else{ 420 h_Iso_->Fill(Tracker_isovalue_mu0); 421 h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu0); 422 } 423 } 424 }// end MC match 425 }//end of cuts 426 }// end loop on ZMuSa cand 427 }// end if ZMuSa size > 0 428 429 //ZMuTk 430 if (zMuTk->size() > 0 ) { 431 event.getByLabel(zMuTkMatchMap_, zMuTkMatchMap); 432 for(unsigned int i = 0; i < zMuTk->size(); ++i) { //loop on candidates 433 const Candidate & zMuTkCand = (*zMuTk)[i]; //the candidate 434 CandidateBaseRef zMuTkCandRef = zMuTk->refAt(i); 435 CandidateBaseRef dau0 = zMuTkCand.daughter(0)->masterClone(); 436 CandidateBaseRef dau1 = zMuTkCand.daughter(1)->masterClone(); 437 const pat::Muon& mu0 = dynamic_cast<const pat::Muon&>(*dau0);//cast in patMuon 438 const pat::GenericParticle& mu1 = dynamic_cast<const pat::GenericParticle &>(*dau1); 439 440 441 double zmass= zMuTkCand.mass(); 442 double pt0 = mu0.pt(); 443 double pt1 = mu1.pt(); 444 double eta0 = mu0.eta(); 445 double eta1 = mu1.eta(); 446 if(pt0>20 && pt1 > 20 && abs(eta0)<2 && abs(eta1)<2 && zmass > 20 && zmass < 200){//kinematical cuts 447 GenParticleRef zMuTkMatch = (*zMuTkMatchMap)[zMuTkCandRef]; 448 if(zMuTkMatch.isNonnull()) { // ZMuTk matched 449 FSR_mu = false; 450 FSR_tk = false; 451 cout<<" ZmuTk cuts && matched"<<endl; 452 zmtcounter++; 453 //Isodeposit 454 const pat::IsoDeposit * muTrackIso =mu0.isoDeposit(pat::TrackIso); 455 const pat::IsoDeposit * tkTrackIso =mu1.isoDeposit(pat::TrackIso); 456 Direction muDir = Direction(mu0.eta(),mu0.phi()); 457 Direction tkDir = Direction(mu1.eta(),mu1.phi()); 458 459 IsoDeposit::AbsVetos vetos_mu; 460 vetos_mu.push_back(new ConeVeto( muDir, dRVeto_ )); 461 vetos_mu.push_back(new ThresholdVeto( ptThreshold_ )); 462 463 reco::IsoDeposit::AbsVetos vetos_tk; 464 vetos_tk.push_back(new ConeVeto( tkDir, dRVeto_ )); 465 vetos_tk.push_back(new ThresholdVeto( ptThreshold_ )); 466 467 double Tracker_isovalue_mu = muTrackIso->sumWithin(dRTrk_,vetos_mu); 468 double Tracker_isovalue_tk = tkTrackIso->sumWithin(dRTrk_,vetos_tk); 469 470 //MonteCarlo Study 471 const reco::GenParticle * muMc0 = mu0.genLepton(); 472 const reco::GenParticle * muMc1 = mu1.genParticle() ; 473 const Candidate * motherMu0 = muMc0->mother(); 474 const Candidate * motherMu1 = muMc1->mother(); 475 int num_dau_muon0 = motherMu0->numberOfDaughters(); 476 int num_dau_muon1 = motherMu1->numberOfDaughters(); 477 cout<<"numero di figli muone0 = " << num_dau_muon0 <<endl; 478 cout<<"numero di figli muone1 = " << num_dau_muon1 <<endl; 479 480 cout<<" muon"<<endl; 481 cout<<" num di daughters = "<< num_dau_muon0 <<endl; 482 if( num_dau_muon0 > 1 ){ 483 for(int j = 0; j < num_dau_muon0; ++j){ 484 int id = motherMu0->daughter(j)->pdgId(); 485 cout<<" dau["<<j<<"] pdg ID = "<<id<<endl; 486 if(id == 22) { 487 FSR_mu=true; 488 } 489 } 490 }//end check of gamma 491 else cout<<" dau[0] pdg ID = "<<motherMu0->daughter(0)->pdgId()<<endl; 492 cout<<" traccia"<<endl; 493 cout<<" num di daughters = "<< num_dau_muon1 <<endl; 494 if( num_dau_muon1 > 1 ){ 495 for(int j = 0; j < num_dau_muon1; ++j){ 496 int id = motherMu1->daughter(j)->pdgId(); 497 cout<<" dau["<<j<<"] pdg ID = "<<id<<endl; 498 if(id == 22) { 499 FSR_tk=true; 500 } 501 } 502 }//end check of gamma 503 else cout<<" dau[0] pdg ID = "<<motherMu1->daughter(0)->pdgId()<<endl; 504 cout<<"Mu Isolation = "<< Tracker_isovalue_mu <<endl; 505 cout<<"Track Isolation = "<< Tracker_isovalue_tk <<endl; 506 if(FSR_mu){ 507 h_Iso_FSR_->Fill(Tracker_isovalue_mu); 508 h_Iso_FSR_3D_->Fill(pt0,eta0,Tracker_isovalue_mu); 509 } 510 else{ 511 h_Iso_->Fill( Tracker_isovalue_mu); 512 h_Iso_3D_->Fill(pt0,eta0,Tracker_isovalue_mu); 513 514 } 515 if(FSR_tk){ 516 h_Iso_FSR_->Fill(Tracker_isovalue_tk); 517 h_Iso_FSR_3D_->Fill(pt1,eta1,Tracker_isovalue_tk); 518 h_trackProbe_eta_FSR->Fill(eta1); 519 h_trackProbe_pt_FSR->Fill(pt1); 520 } 521 else{ 522 h_Iso_->Fill( Tracker_isovalue_tk); 523 h_Iso_3D_->Fill(pt1,eta1,Tracker_isovalue_tk); 524 h_trackProbe_eta_no_FSR->Fill(eta1); 525 h_trackProbe_pt_no_FSR->Fill(pt1); 526 } 527 }// end MC match 528 }//end Kine-cuts 529 }// end loop on ZMuTk cand 530 }// end if ZMuTk size > 0 531}// end analyze 532 533 534 535void ZMuMu_Radiative_analyzer::endJob() { 536 cout<<" ============= Summary =========="<<endl; 537 cout <<" Numero di eventi "<< evntcounter << endl; 538 cout <<" 1)Numero di ZMuMu matched dopo i tagli cinematici = "<< zmmcounter << endl; 539 cout <<" 2)Numero di ZMuSa matched dopo i tagli cinematici = "<< zmscounter << endl; 540 cout <<" 3)Numero di ZMuTk matched dopo i tagli cinematici = "<< zmtcounter << endl; 541 double n1= h_Iso_FSR_->Integral(); 542 double icut1= h_Iso_FSR_->Integral(0,15); 543 double eff_iso_FSR = (double)icut1/(double)n1; 544 double err_iso_FSR = sqrt(eff_iso_FSR*(1-eff_iso_FSR)/n1); 545 double n2= h_Iso_->Integral(); 546 double icut2= h_Iso_->Integral(0,15); 547 double eff_iso= (double)icut2/(double)n2; 548 double err_iso = sqrt(eff_iso*(1-eff_iso)/n2); 549 cout<<" ============= Isolation Efficiecy =========="<<endl; 550 cout<<"Isolation Efficiency = "<< eff_iso <<" +/- "<< err_iso <<endl; 551 cout<<"Isolation Efficiency with FSR = "<< eff_iso_FSR <<" +/- "<< err_iso_FSR <<endl; 552 553 } 554 555#include "FWCore/Framework/interface/MakerMacros.h" 556 557DEFINE_FWK_MODULE(ZMuMu_Radiative_analyzer); 558