Reference Manual

Go to the documentation of this file.
00001 //  \class MuScleFitPlotter
00002 //  Plotter for simulated,generated and reco info of muons
00003 //
00004 //  $Date: 2010/10/22 17:47:44 $
00005 //  $Revision: 1.22 $
00006 //  \author  C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo, M.De Mattia - INFN Padova
00007 //
00008 // ----------------------------------------------------------------------------------
00009 
00010 #include "MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitPlotter.h"
00011 #include "MuonAnalysis/MomentumScaleCalibration/interface/Histograms.h"
00012 #include "MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitUtils.h"
00013 #include "FWCore/Framework/interface/MakerMacros.h"
00014 #include "FWCore/Framework/interface/Frameworkfwd.h"
00015 #include "FWCore/Framework/interface/Event.h"
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 //#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00018 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00019 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00020 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00021 #include "DataFormats/TrackReco/interface/Track.h"
00022 #include "DataFormats/MuonReco/interface/Muon.h"
00023 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00024 #include "DataFormats/Candidate/interface/LeafCandidate.h"
00025 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00026 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00027 #include <CLHEP/Vector/LorentzVector.h>
00028 
00029 #include "TFile.h"
00030 #include "TTree.h"
00031 #include "TMinuit.h"
00032 #include <vector>
00033 
00034 // Constructor
00035 // ----------
00036 MuScleFitPlotter::MuScleFitPlotter(std::string theGenInfoRootFileName){
00037   outputFile = new TFile(theGenInfoRootFileName.c_str(),"RECREATE");
00038   fillHistoMap();
00039 }
00040 
00041 MuScleFitPlotter::~MuScleFitPlotter(){
00042   outputFile->cd();
00043   writeHistoMap();
00044   outputFile->Close();
00045 }
00046 
00047 // Find and store in histograms the generated resonance and muons
00048 // --------------------------------------------------------------
00049 void MuScleFitPlotter::fillGen1(const reco::GenParticleCollection* genParticles, bool PATmuons)
00050 {
00051   //  bool prova = false;
00052   //Loop on generated particles
00053   std::pair<reco::Particle::LorentzVector,reco::Particle::LorentzVector> muFromRes;
00054   reco::Particle::LorentzVector genRes;
00055 
00056   int mothersFound[] = {0, 0, 0, 0, 0, 0};
00057 
00058   for( reco::GenParticleCollection::const_iterator mcIter=genParticles->begin(); mcIter!=genParticles->end(); ++mcIter ) {
00059     int status = mcIter->status();
00060     int pdgId = std::abs(mcIter->pdgId());
00061     //Check if it's a resonance
00062     if( status == 2 &&
00063         ( pdgId==23  || pdgId==443    || pdgId==100443 ||
00064           pdgId==553 || pdgId==100553 || pdgId==200553 ) ) {
00065       genRes = mcIter->p4();
00066       // std::cout << "mother's mother = " << mcIter->mother()->pdgId() << std::endl;
00067       if( pdgId == 23 ) mapHisto["hGenResZ"]->Fill(genRes);
00068       else if( pdgId == 443 || pdgId == 100443 ) mapHisto["hGenResJPsi"]->Fill(genRes);
00069       else if( pdgId == 553 || pdgId == 100553 || pdgId == 200553 ) mapHisto["hGenResUpsilon1S"]->Fill(genRes);
00070     }
00071     //Check if it's a muon from a resonance
00072     if( status==1 && pdgId==13 && !PATmuons) {
00073       int momPdgId = std::abs(mcIter->mother()->pdgId());
00074       if( momPdgId==23  || momPdgId==443    || momPdgId==100443 || 
00075           momPdgId==553 || momPdgId==100553 || momPdgId==200553 ) {
00076         if( momPdgId == 23 ) mothersFound[0] = 1;
00077         if( momPdgId == 443 || momPdgId == 100443 ) mothersFound[5] = 1;
00078         if( momPdgId == 553 || momPdgId == 100553 || momPdgId == 200553 ) mothersFound[3] = 1;
00079         mapHisto["hGenMu"]->Fill(mcIter->p4());
00080         std::cout<<"genmu "<<mcIter->p4()<<std::endl;
00081         if(mcIter->charge()>0){
00082           muFromRes.first = mcIter->p4();
00083           // prova = true;
00084         }
00085         else muFromRes.second = mcIter->p4();
00086       }
00087     }//if PATmuons you don't have the info of the mother !!! Here I assume is a JPsi
00088     if( status==1 && pdgId==13 && PATmuons) {
00089       mothersFound[5] = 1;
00090       mapHisto["hGenMu"]->Fill(mcIter->p4());
00091       std::cout<<"genmu "<<mcIter->p4()<<std::endl;
00092       if(mcIter->charge()>0){
00093         muFromRes.first = mcIter->p4();
00094         // prova = true;
00095       }
00096       else muFromRes.second = mcIter->p4();
00097     }
00098   }
00099   //   if(!prova)
00100   //     std::cout<<"hgenmumu not found"<<std::endl;
00101 
00102   if( mothersFound[0] == 1 ) {
00103     mapHisto["hGenMuMuZ"]->Fill(muFromRes.first+muFromRes.second);
00104     mapHisto["hGenResVSMuZ"]->Fill( muFromRes.first, genRes, 1 );
00105     mapHisto["hGenResVSMuZ"]->Fill( muFromRes.second,genRes, -1 );
00106   }
00107   if( mothersFound[3] == 1 ) {
00108     mapHisto["hGenMuMuUpsilon1S"]->Fill(muFromRes.first+muFromRes.second);
00109     mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.first, genRes, 1 );
00110     mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.second,genRes, -1 );
00111   }
00112   if( mothersFound[5] == 1 ) {
00113     mapHisto["hGenMuMuJPsi"]->Fill(muFromRes.first+muFromRes.second);
00114     mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.first, genRes, 1 );
00115     mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.second,genRes, -1 );
00116   }
00117 
00118   mapHisto["hGenResVsSelf"]->Fill( genRes, genRes, 1 );
00119 }
00120 
00121 // Find and store in histograms the generated resonance and muons
00122 // --------------------------------------------------------------
00123 void MuScleFitPlotter::fillGen2(const edm::HepMCProduct* evtMC, bool sherpaFlag_)
00124 {
00125   //Loop on generated particles
00126   const HepMC::GenEvent* Evt = evtMC->GetEvent();
00127   std::pair<reco::Particle::LorentzVector,reco::Particle::LorentzVector> muFromRes; 
00128   reco::Particle::LorentzVector genRes;
00129   
00130   int mothersFound[] = {0, 0, 0, 0, 0, 0};
00131   
00132   if( sherpaFlag_ ) {
00133     
00134     for (HepMC::GenEvent::particle_const_iterator part=Evt->particles_begin(); 
00135          part!=Evt->particles_end(); part++) {
00136       if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {//looks for muon in the final state
00137         bool fromRes = false;
00138         for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(HepMC::ancestors);//loops on the mother of the final state muons
00139              mother != (*part)->production_vertex()->particles_end(HepMC::ancestors); ++mother) {
00140           unsigned int motherPdgId = (*mother)->pdg_id();
00141           if( motherPdgId == 13 && (*mother)->status() == 3 ) fromRes = true;
00142         }
00143         if(fromRes){
00144           if((*part)->pdg_id()==13){
00145             muFromRes.first = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00146                                              (*part)->momentum().pz(),(*part)->momentum().e()));
00147           }
00148           else if((*part)->pdg_id()==-13){
00149             muFromRes.second = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00150                                               (*part)->momentum().pz(),(*part)->momentum().e()));           
00151           }
00152         }
00153       }
00154       
00155     }
00156     mapHisto["hGenResZ"]->Fill(muFromRes.first+muFromRes.second);   
00157   }
00158   else{
00159     for (HepMC::GenEvent::particle_const_iterator part=Evt->particles_begin(); 
00160          part!=Evt->particles_end(); part++) {
00161       int status = (*part)->status();
00162       int pdgId = std::abs((*part)->pdg_id());
00163       //std::cout<<"PDG ID "<< (*part)->pdg_id() <<"    status "<< (*part)->status()
00164       //<<"   pt "<<(*part)->momentum().perp()<< "     eta  "<<(*part)->momentum().eta()<<std::endl    ;
00165       //Check if it's a resonance       
00166       // For sherpa the resonance is not saved. The muons from the resonance can be identified
00167       // by having as mother a muon of status 3.
00168       
00169       if (pdgId==13 && status==1) {  
00170         if( status==2 && 
00171             ( pdgId==23  || pdgId==443    || pdgId==100443 ||
00172               pdgId==553 || pdgId==100553 || pdgId==200553 ) ) {
00173           genRes = reco::Particle::LorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00174                                                  (*part)->momentum().pz(),(*part)->momentum().e());
00175           
00176           if( pdgId == 23 ) mapHisto["hGenResZ"]->Fill(genRes);
00177           if( pdgId == 443 ) mapHisto["hGenResJPsi"]->Fill(genRes);
00178           if( pdgId == 553 ) {
00179             // std::cout << "genRes mass = " << CLHEP::HepLorentzVector(genRes.x(),genRes.y(),genRes.z(),genRes.t()).m() << std::endl;
00180             mapHisto["hGenResUpsilon1S"]->Fill(genRes);
00181           }
00182         }
00183       
00184         //Check if it's a muon from a resonance
00185         if (pdgId==13 && status==1) {      
00186           bool fromRes=false;
00187           for (HepMC::GenVertex::particle_iterator mother = 
00188                  (*part)->production_vertex()->particles_begin(HepMC::ancestors);
00189                mother != (*part)->production_vertex()->particles_end(HepMC::ancestors); ++mother) {
00190             int motherPdgId = (*mother)->pdg_id();
00191             if (motherPdgId==23  || motherPdgId==443    || motherPdgId==100443 || 
00192                 motherPdgId==553 || motherPdgId==100553 || motherPdgId==200553) {
00193               fromRes=true;
00194               if( motherPdgId == 23 ) mothersFound[0] = 1;
00195               if( motherPdgId == 443 ) mothersFound[3] = 1;
00196               if( motherPdgId == 553 ) mothersFound[5] = 1;
00197             }
00198           }
00199           
00200           if(fromRes) { 
00201             mapHisto["hGenMu"]->Fill(reco::Particle::LorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00202                                                                    (*part)->momentum().pz(),(*part)->momentum().e()));
00203             mapHisto["hGenMuVSEta"]->Fill(reco::Particle::LorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00204                                                                         (*part)->momentum().pz(),(*part)->momentum().e()));
00205             if((*part)->pdg_id()==-13)
00206               muFromRes.first = (reco::Particle::LorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00207                                                                (*part)->momentum().pz(),(*part)->momentum().e()));
00208             else
00209               muFromRes.second = (reco::Particle::LorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
00210                                                                 (*part)->momentum().pz(),(*part)->momentum().e()));
00211           }
00212         }
00213       }
00214     }
00215   }
00216   if( mothersFound[0] == 1 ) {
00217     mapHisto["hGenMuMuZ"]->Fill(muFromRes.first+muFromRes.second);
00218     mapHisto["hGenResVSMuZ"]->Fill( muFromRes.first, genRes, 1 );
00219     mapHisto["hGenResVSMuZ"]->Fill( muFromRes.second,genRes, -1 );
00220   }
00221   if( mothersFound[3] == 1 ) {
00222     mapHisto["hGenMuMuUpsilon1S"]->Fill(muFromRes.first+muFromRes.second);
00223     mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.first, genRes, 1 );
00224     mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.second,genRes, -1 );
00225   }
00226   if( mothersFound[5] == 1 ) {
00227     mapHisto["hGenMuMuJPsi"]->Fill(muFromRes.first+muFromRes.second);
00228     mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.first, genRes, 1 );
00229     mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.second,genRes, -1 );
00230   }
00231   mapHisto["hGenResVsSelf"]->Fill( genRes, genRes, 1 );
00232 }
00233 
00234 // Find and store in histograms the simulated resonance and muons
00235 // --------------------------------------------------------------
00236 void MuScleFitPlotter::fillSim(edm::Handle<edm::SimTrackContainer> simTracks)
00237 {
00238   std::vector<SimTrack> simMuons;
00239 
00240   //Loop on simulated tracks
00241   for( edm::SimTrackContainer::const_iterator simTrack=simTracks->begin(); simTrack!=simTracks->end(); ++simTrack ) {
00242     // Select the muons from all the simulated tracks
00243     if (fabs((*simTrack).type())==13) {
00244       simMuons.push_back(*simTrack);      
00245       mapHisto["hSimMu"]->Fill((*simTrack).momentum());
00246     }
00247   }
00248   mapHisto["hSimMu"]->Fill(simMuons.size());
00249 
00250   // Recombine all the possible Z from simulated muons
00251   if( simMuons.size() >= 2 ) {
00252     for( std::vector<SimTrack>::const_iterator  imu=simMuons.begin(); imu != simMuons.end(); ++imu ) {   
00253       for( std::vector<SimTrack>::const_iterator imu2=imu+1; imu2!=simMuons.end(); ++imu2 ) {
00254         if (imu==imu2) continue;
00255 
00256         // Try all the pairs with opposite charge
00257         if (((*imu).charge()*(*imu2).charge())<0) {
00258           reco::Particle::LorentzVector Z = (*imu).momentum()+(*imu2).momentum();
00259           mapHisto["hSimMuPMuM"]->Fill(Z); 
00260         }
00261       }
00262     }
00263 
00264     // Plots for the best possible simulated resonance
00265     std::pair<SimTrack,SimTrack> simMuFromBestRes = MuScleFitUtils::findBestSimuRes(simMuons);
00266     reco::Particle::LorentzVector bestSimZ = (simMuFromBestRes.first).momentum()+(simMuFromBestRes.second).momentum();
00267     mapHisto["hSimBestRes"]->Fill(bestSimZ);
00268     if (fabs(simMuFromBestRes.first.momentum().eta())<2.5 && fabs(simMuFromBestRes.second.momentum().eta())<2.5 &&
00269         simMuFromBestRes.first.momentum().pt()>2.5 && simMuFromBestRes.second.momentum().pt()>2.5) {
00270       mapHisto["hSimBestResVSMu"]->Fill (simMuFromBestRes.first.momentum(), bestSimZ, int(simMuFromBestRes.first.charge()));
00271       mapHisto["hSimBestResVSMu"]->Fill (simMuFromBestRes.second.momentum(),bestSimZ, int(simMuFromBestRes.second.charge()));
00272     }
00273   }
00274 }
00275 
00276 // Find and store in histograms the RIGHT simulated resonance and muons
00277 // --------------------------------------------------------------
00278 void MuScleFitPlotter::fillGenSim(edm::Handle<edm::HepMCProduct> evtMC, edm::Handle<edm::SimTrackContainer> simTracks)
00279 {
00280   std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> simMuFromRes = 
00281     MuScleFitUtils::findSimMuFromRes(evtMC,simTracks);
00282   //Fill resonance info
00283   reco::Particle::LorentzVector rightSimRes = (simMuFromRes.first)+(simMuFromRes.second);
00284   mapHisto["hSimRightRes"]->Fill(rightSimRes);
00285   /*if ((fabs(simMuFromRes.first.Eta())<2.5 && fabs(simMuFromRes.second.Eta())<2.5) 
00286     && simMuFromRes.first.Pt()>2.5 && simMuFromRes.second.Pt()>2.5) {
00287   }*/
00288 }
00289 
00290 // Find and store in histograms the reconstructed resonance and muons
00291 // --------------------------------------------------------------
00292 void MuScleFitPlotter::fillRec(std::vector<reco::LeafCandidate>& muons)
00293 {
00294   for(std::vector<reco::LeafCandidate>::const_iterator mu1 = muons.begin(); mu1!=muons.end(); mu1++){
00295     mapHisto["hRecMu"]->Fill(mu1->p4());
00296     mapHisto["hRecMuVSEta"]->Fill(mu1->p4());
00297     for(std::vector<reco::LeafCandidate>::const_iterator mu2 = muons.begin(); mu2!=muons.end(); mu2++){  
00298       if (mu1->charge()<0 || mu2->charge()>0)
00299         continue;
00300       reco::Particle::LorentzVector Res (mu1->p4()+mu2->p4());
00301       mapHisto["hRecMuPMuM"]->Fill(Res);          
00302     } 
00303   }
00304   mapHisto["hRecMu"]->Fill(muons.size());
00305 }
00306 
00308 void MuScleFitPlotter::fillRec( const std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> > & savedPairs )
00309 {
00310   std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >::const_iterator muonPair = savedPairs.begin();
00311   for( ; muonPair != savedPairs.end(); ++muonPair ) {
00312     mapHisto["hRecMu"]->Fill(muonPair->first);
00313     mapHisto["hRecMuVSEta"]->Fill(muonPair->first);
00314     mapHisto["hRecMu"]->Fill(muonPair->second);
00315     mapHisto["hRecMuVSEta"]->Fill(muonPair->second);
00316     reco::Particle::LorentzVector Res( muonPair->first+muonPair->second );
00317     mapHisto["hRecMuPMuM"]->Fill(Res);
00318     mapHisto["hRecMu"]->Fill(savedPairs.size());
00319   }
00320 }
00321 
00327 void MuScleFitPlotter::fillGen( const std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> > & genPairs )
00328 {
00329   std::vector<std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> >::const_iterator genPair = genPairs.begin();
00330   for( ; genPair != genPairs.end(); ++genPair ) {
00331     reco::Particle::LorentzVector genRes(genPair->first+genPair->second);
00332     mapHisto["hGenResZ"]->Fill(genRes);
00333     mapHisto["hGenMu"]->Fill(genPair->first);
00334     mapHisto["hGenMuVSEta"]->Fill(genPair->first);
00335     mapHisto["hGenMuMuZ"]->Fill(genRes);
00336     mapHisto["hGenResVSMuZ"]->Fill( genPair->first, genRes, 1 );
00337     mapHisto["hGenResVSMuZ"]->Fill( genPair->second, genRes, -1 );
00338     mapHisto["hGenMuMuUpsilon1S"]->Fill(genRes);
00339     mapHisto["hGenResVSMuUpsilon1S"]->Fill( genPair->first, genRes, 1 );
00340     mapHisto["hGenResVSMuUpsilon1S"]->Fill( genPair->second, genRes, -1 );
00341     mapHisto["hGenMuMuJPsi"]->Fill(genRes);
00342     mapHisto["hGenResVSMuJPsi"]->Fill( genPair->first, genRes, 1 );
00343     mapHisto["hGenResVSMuJPsi"]->Fill( genPair->second, genRes, -1 );
00344     mapHisto["hGenResVsSelf"]->Fill( genRes, genRes, 1 );
00345   }
00346 }
00347 
00348 // Histogram booking
00349 // -----------------
00350 void MuScleFitPlotter::fillHistoMap() {
00351 
00352   // Generated Z and muons
00353   // ---------------------
00354   mapHisto["hGenResJPsi"]      = new HParticle   ("hGenResJPsi", 3., 3.1);
00355   mapHisto["hGenResUpsilon1S"] = new HParticle   ("hGenResUpsilon1S", 9., 11.);
00356   mapHisto["hGenResZ"]         = new HParticle   ("hGenResZ", 60., 120.);
00357   mapHisto["hGenMu"]      = new HParticle  ("hGenMu");
00358   mapHisto["hGenMuVSEta"] = new HPartVSEta ("hGenMuVSEta");
00359 
00360   mapHisto["hGenMuMuJPsi"]      = new HParticle   ("hGenMuMuJPsi",3., 3.1 );
00361   mapHisto["hGenResVSMuJPsi"]   = new HMassVSPart ("hGenResVSMuJPsi",3., 3.1);
00362   mapHisto["hGenMuMuUpsilon1S"]      = new HParticle   ("hGenMuMuUpsilon1S", 9., 11.);
00363   mapHisto["hGenResVSMuUpsilon1S"]   = new HMassVSPart ("hGenResVSMuUpsilon1S", 9., 11.);
00364   mapHisto["hGenMuMuZ"]      = new HParticle   ("hGenMuMuZ", 60., 120.);
00365   mapHisto["hGenResVSMuZ"]   = new HMassVSPart ("hGenResVSMuZ", 60., 120.);
00366 
00367   mapHisto["hGenResVsSelf"] = new HMassVSPart ("hGenResVsSelf");
00368 
00369   // Simulated resonance and muons
00370   // -----------------------------
00371   mapHisto["hSimMu"]      = new HParticle ("hSimMu");
00372 
00373   mapHisto["hSimMuPMuM"]      = new HParticle ("hSimMuPMuM");      
00374                                                                  
00375   mapHisto["hSimBestMu"]      = new HParticle ("hSimBestMu");
00376   mapHisto["hSimBestRes"]         = new HParticle  ("hSimBestRes");
00377   mapHisto["hSimBestResVSMu"]  = new HMassVSPart ("hSimBestResVSMu");
00378     
00379   mapHisto["hSimRightRes"]         = new HParticle  ("hSimRightZ");
00380  
00381   // Reconstructed resonance and muons
00382   // -----------------------------  
00383   mapHisto["hRecMu"]      = new HParticle ("hRecMu");
00384   mapHisto["hRecMuVSEta"]      = new HPartVSEta ("hRecMuVSEta");
00385   mapHisto["hRecMuPMuM"]         = new HParticle  ("hRecMuPMuM");
00386 }  
00387 
00388 
00389 // Histogram saving
00390 // -----------------
00391 void MuScleFitPlotter::writeHistoMap() {
00392   outputFile->cd();
00393   for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto.begin(); 
00394        histo!=mapHisto.end(); histo++) {
00395     (*histo).second->Write();
00396   }
00397 }
00398
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitPlotter.cc
