Reference Manual

Go to the documentation of this file.
00001 #include "TauAnalysis/MCEmbeddingTools/plugins/SelectReplacementCandidates.h"
00002 #include "DataFormats/Math/interface/deltaR.h"
00003 
00004 SelectReplacementCandidates::SelectReplacementCandidates(const edm::ParameterSet& iConfig)
00005 {
00006         using namespace edm;
00007         using namespace std;
00008         produces< vector<uint32_t> >();
00009         produces< vector<uint32_t> >("assocHitsWithHCAL");
00010         produces< std::vector<reco::Muon> >();
00011         
00012         targetParticleMass_=iConfig.getUntrackedParameter<double>("targetParticlesMass",1.77690);
00013         targetParticlePdgID_=iConfig.getUntrackedParameter<int>("targetParticlesPdgID",15);
00014         
00015         muonInputTag_ = iConfig.getParameter<edm::InputTag>("muonInputTag");
00016         edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00017         parameters_.loadParameters( parameters );
00018         
00019         trackAssociator_.useDefaultPropagator();
00020 }
00021 
00022 SelectReplacementCandidates::~SelectReplacementCandidates()
00023 {
00024 
00025 }
00026 
00027 void SelectReplacementCandidates::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00028 {
00029         using namespace std;
00030         using namespace edm;
00031         reco::Muon muon1, muon2;
00032 
00033         // determine the muons to be used
00034         // these can be reconstructed muons or muons from a HepMC::Event
00035         if (determineMuonsToUse(iEvent, iSetup, &muon1, &muon2)!=0)
00036                 return;
00037         
00038         vector<uint32_t> hits ; 
00039         getRawIDsAdvanced(iEvent, iSetup, &hits, &muon1, false);
00040         getRawIDsAdvanced(iEvent, iSetup, &hits, &muon2, false);
00041         std::auto_ptr< vector<uint32_t> > selectedHitsAutoPtr(new vector<uint32_t>(hits) );
00042         iEvent.put( selectedHitsAutoPtr );
00043 
00044         vector<uint32_t> assoc_hits_withHCAL;   
00045         getRawIDsAdvanced(iEvent, iSetup, &assoc_hits_withHCAL, &muon1, true);
00046         getRawIDsAdvanced(iEvent, iSetup, &assoc_hits_withHCAL, &muon2, true);
00047         std::cout << "found in total " << assoc_hits_withHCAL.size() << " cells with associated hits with hcal\n";
00048         std::auto_ptr< vector<uint32_t> > selectedAssocHitsWithHCALAutoPtr(new vector<uint32_t>(assoc_hits_withHCAL) );
00049         iEvent.put( selectedAssocHitsWithHCALAutoPtr, "assocHitsWithHCAL");
00050                                                         
00051         std::vector<reco::Muon> muons;
00052         transformMuMu2TauTau(&muon1, &muon2);
00053         muons.push_back(muon1);
00054         muons.push_back(muon2);
00055         std::auto_ptr< std::vector<reco::Muon> > selectedMuonsPtr(new std::vector<reco::Muon>(muons) );
00056         iEvent.put( selectedMuonsPtr );
00057 
00058 }
00059 
00060 void SelectReplacementCandidates::beginJob()
00061 {
00062 
00063 }
00064 
00065 void SelectReplacementCandidates::endJob()
00066 {
00067 
00068 }
00069 
00070 
00071 int SelectReplacementCandidates::determineMuonsToUse(const edm::Event& iEvent, const edm::EventSetup& iConfig, reco::Muon * muon1, reco::Muon * muon2)
00072 {
00073         using namespace edm;
00074         using namespace reco;
00075         using namespace std;
00076 
00077         Handle< edm::RefToBaseVector<reco::Candidate> > zCandidate_handle;
00078         if (!iEvent.getByLabel("dimuonsGlobal", zCandidate_handle))
00079         {
00080                 std::cout << "Could not find product: " << "dimuonsGlobal" << "\n";    
00081     std::vector< edm::Handle< edm::RefToBaseVector<reco::Candidate> >  > allHandles;
00082     iEvent.getManyByType(allHandles);
00083     std::vector< edm::Handle< edm::RefToBaseVector<reco::Candidate> > >::iterator it;
00084     for (it = allHandles.begin(); it != allHandles.end(); it++)
00085     {
00086       std::cout << "available product: " << (*it).provenance()->moduleLabel() << ", " << (*it).provenance()->productInstanceName() << ", " << (*it).provenance()->processName();
00087     }
00088 
00089           std::cout << "Objekt nicht gefunden: dimuonsGloal\n";
00090           return -1;
00091         }
00092 //      std::cout << zCandidate_handle->size() << " Kandidaten gefunden!\n";
00093 
00094         unsigned int nMuons = zCandidate_handle->size();
00095         if (nMuons==0)
00096                 return -1;
00097 
00098         for (edm::RefToBaseVector<reco::Candidate>::const_iterator z = zCandidate_handle->begin(); z!=zCandidate_handle->end(); ++z)
00099         {
00100                 reco::Particle::LorentzVector muon1p4 = z->get()->daughter(0)->p4();
00101                 reco::Particle::LorentzVector muon2p4 = z->get()->daughter(1)->p4();
00102 
00103                 edm::Handle<edm::View<reco::Muon> > trackCollection;
00104                 iEvent.getByLabel(muonInputTag_, trackCollection);
00105                 const edm::View<reco::Muon>& muons = * trackCollection;
00106 
00107                 bool found1=false, found2=false;
00108                 for (unsigned int i=0;i<muons.size() && !(found1 && found2);i++)
00109                 {
00110                         if (deltaR(muon1p4,muons[i].p4())<0.1)
00111                         {
00112                                 *muon1 = muons[i];
00113                                 found1=true;
00114                         }
00115                         if (deltaR(muon2p4,muons[i].p4())<0.1)
00116                         {
00117                                 *muon2 = muons[i];
00118                                 found2=true;
00119                         }
00120                 }
00121 
00122                 break;
00123         }
00124         return 0;
00125 }
00126 
00127 
00128 int SelectReplacementCandidates::determineMuonsToUse_old(const edm::Event& iEvent, const edm::EventSetup& iConfig, reco::Muon * muon1, reco::Muon * muon2)
00129 {
00130     using namespace edm;
00131     using namespace reco;
00132     using namespace std;
00133 
00134         Handle<MuonCollection> muonHandle;
00135         if (!iEvent.getByLabel(muonInputTag_,muonHandle))
00136                 ProductNotFound<reco::MuonCollection>(iEvent, muonInputTag_);
00137         
00138         edm::Handle<edm::View<reco::Muon> > trackCollection;
00139         iEvent.getByLabel(muonInputTag_, trackCollection);
00140         const edm::View<reco::Muon>& muons = * trackCollection;
00141 
00142         unsigned int nMuons = muons.size();
00143         if (nMuons<2)
00144                 return -1;
00145         
00146         *muon1 = muons[0];
00147         *muon2 = muons[1];
00148         std::cout << muons[0].p4() << "\n";
00149         std::cout << muons[1].p4() << "\n";
00150     
00151   return 0;
00152 }
00153 
00154 void SelectReplacementCandidates::getRawIDsAdvanced(const edm::Event& iEvent, const edm::EventSetup& iConfig, std::vector<uint32_t> * L, reco::Muon * muon, bool includeHCAL)
00155 {
00156   using namespace edm;
00157   using namespace reco;
00158   using namespace std;
00159 
00160         if (muon->bestTrackRef().isNonnull())
00161         {
00162                 TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iConfig, *(muon->bestTrackRef().get()), parameters_);
00163                 if (includeHCAL)
00164                 {
00165                         for(std::vector<const EcalRecHit*>::const_iterator hit = info.ecalRecHits.begin(); hit != info.ecalRecHits.end(); ++hit)
00166         L->push_back((*hit)->detid().rawId());
00167                                                                                                 
00168                         for(std::vector<const HBHERecHit*>::const_iterator hit = info.hcalRecHits.begin(); hit != info.hcalRecHits.end(); ++hit)
00169                                 L->push_back((*hit)->detid().rawId());
00170                                                                                                                                                               }
00171                         int recs = (muon->bestTrackRef().get())->recHitsSize();
00172                         for (int i=0; i<recs; i++)
00173                                 L->push_back(((muon->bestTrackRef().get())->recHit(i).get())->geographicalId().rawId());
00174         }
00175         else
00176                 std::cout << "ERROR: Muon has no bestTrackRef!!\n";
00177 
00178     std::cout << " with " <<  L->size() << " muon hits found\n";
00179 //     return L;
00180 }
00181 
00182 
00183 template < typename T > void SelectReplacementCandidates::ProductNotFound(const edm::Event& iEvent, edm::InputTag inputTag)
00184 {
00185     std::cout << "--- TauHitSelector ------------------------------------\n";
00186     std::cout << "Could not find product with:\n"<< inputTag << "\n";
00187     std::vector< edm::Handle< T > > allHandles; 
00188     iEvent.getManyByType(allHandles);
00189     typename std::vector< edm::Handle< T > >::iterator it;
00190     for (it = allHandles.begin(); it != allHandles.end(); it++)
00191     {
00192             std::cout << "module label:        " << (*it).provenance()->moduleLabel() << "\n";
00193             std::cout << "productInstanceName: " << (*it).provenance()->productInstanceName() << "\n";
00194             std::cout << "processName:         " << (*it).provenance()->processName() << "\n\n";
00195     }
00196     std::cout << "-------------------------------------------------------\n";           
00197                 return;
00198 }
00199 
00200 
00202 void SelectReplacementCandidates::transformMuMu2TauTau(reco::Muon * muon1, reco::Muon * muon2)
00203 {
00204         using namespace edm;
00205         using namespace reco;
00206         using namespace std;
00207         
00208         reco::Particle::LorentzVector muon1_momentum = muon1->p4();
00209         reco::Particle::LorentzVector muon2_momentum =  muon2->p4();
00210         reco::Particle::LorentzVector z_momentum = muon1_momentum + muon2_momentum;
00211 
00212         ROOT::Math::Boost booster(z_momentum.BoostToCM());
00213         ROOT::Math::Boost invbooster(booster.Inverse());
00214         
00215         reco::Particle::LorentzVector Zb = booster(z_momentum);
00216 
00217         reco::Particle::LorentzVector muon1b = booster(muon1_momentum);
00218         reco::Particle::LorentzVector muon2b = booster(muon2_momentum);
00219         
00220         double tau_mass2 = targetParticleMass_*targetParticleMass_;
00221 
00222         double muonxb_mom2 = muon1b.x()*muon1b.x() + muon1b.y()*muon1b.y() + muon1b.z() * muon1b.z();
00223         double tauxb_mom2 = 0.25 * Zb.t() * Zb.t() - tau_mass2;
00224 
00225         float scaling1 = sqrt(tauxb_mom2 / muonxb_mom2);
00226         float scaling2 = scaling1;
00227 
00228         float tauEnergy= Zb.t() / 2.;
00229 
00230         if (tauEnergy*tauEnergy<tau_mass2)
00231                 return;
00232         
00233         reco::Particle::LorentzVector tau1b_mom = reco::Particle::LorentzVector(scaling1*muon1b.x(),scaling1*muon1b.y(),scaling1*muon1b.z(),tauEnergy);
00234         reco::Particle::LorentzVector tau2b_mom = reco::Particle::LorentzVector(scaling2*muon2b.x(),scaling2*muon2b.y(),scaling2*muon2b.z(),tauEnergy);
00235 
00236         // some checks
00237         // the following test guarantees a deviation
00238         // of less than 0.1% for phi and theta for the
00239         // original muons and the placed taus
00240         // (in the centre-of-mass system of the z boson)
00241         assert((muon1b.phi()-tau1b_mom.phi())/muon1b.phi()<0.001);
00242         assert((muon2b.phi()-tau2b_mom.phi())/muon2b.phi()<0.001);
00243         assert((muon1b.theta()-tau1b_mom.theta())/muon1b.theta()<0.001);
00244         assert((muon2b.theta()-tau2b_mom.theta())/muon2b.theta()<0.001);        
00245 
00246         reco::Particle::LorentzVector tau1_mom = (invbooster(tau1b_mom));
00247         reco::Particle::LorentzVector tau2_mom = (invbooster(tau2b_mom));
00248         
00249         // some additional checks
00250         // the following tests guarantee a deviation of less
00251         // than 0.1% for the following values of the original
00252         // muons and the placed taus
00253         //      invariant mass
00254         //      transverse momentum
00255         assert(((muon1_momentum+muon1_momentum).mass()-(tau1_mom+tau2_mom).mass())/(muon1_momentum+muon1_momentum).mass()<0.001);
00256         assert(((muon1_momentum+muon2_momentum).pt()-(tau1_mom+tau2_mom).pt())/(muon1_momentum+muon1_momentum).pt()<0.001);
00257 
00258         muon1->setP4(tau1_mom);
00259         muon2->setP4(tau2_mom);
00260 
00261         muon1->setPdgId(targetParticlePdgID_*muon1->pdgId()/abs(muon1->pdgId()));
00262         muon2->setPdgId(targetParticlePdgID_*muon2->pdgId()/abs(muon2->pdgId()));
00263 
00264         muon1->setStatus(1);
00265         muon2->setStatus(1);
00266 
00267         return;
00268 }
00269 
00270 #include "FWCore/Framework/interface/MakerMacros.h"
00271 
00272 DEFINE_FWK_MODULE(SelectReplacementCandidates);
00273
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/TauAnalysis/MCEmbeddingTools/plugins/SelectReplacementCandidates.cc
