/*  Implementation: 
0. addapted from Ferenc Sickler LowPtAnalyzer
1. hit-by-hit track association
2. in case of muons: separate tracks+sta part matching
3. output: 3 TNtuples for genReco and recoGenTrigger, and recoGenAssociated (in case of dimuon-h) or 2 (for h-h)

*/


// framework includes
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/Exception.h"

// data formats
#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/Common/interface/RefToBase.h"
#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/MuonReco/interface/MuonTrackLinks.h"
#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"


// pixel
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"

// strip
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"

//geometry
#include "Geometry/CommonDetUnit/interface/GeomDetType.h"
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"

// sim data formats
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"

#include "SimMuon/MCTruth/interface/MuonAssociatorByHits.h"
#include "SimTracker/Records/interface/TrackAssociatorRecord.h"
#include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"



//tracking tools
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"

#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
#include "RecoVertex/KalmanVertexFit/interface/SingleTrackVertexConstraint.h"

// ROOT
#include "TROOT.h"
#include "TNtuple.h"
#include "TLorentzVector.h"

// miscellaneous  
#include <fstream>

using namespace std;
using namespace reco;
const double m_mu = .105658;

//__________________________________________________________________________
class McMatchAnalyzer : public edm::EDAnalyzer
{

 public:
  explicit McMatchAnalyzer(const edm::ParameterSet& pset);
  ~McMatchAnalyzer();
  
  virtual void analyze(const edm::Event& ev, const edm::EventSetup& es);
  virtual void beginJob(const edm::EventSetup& es);
  virtual void endJob();
  
private:
  void fillRecoDimuonTuple(const edm::RefToBase<Track>& trkRef,
			   reco::RecoToSimCollection& reco2Sim,
			   vector<float>& result);
  void fillRecoTrackTuple(const edm::RefToBase<Track>& trkRef,
			  vector<float>& result);
  void fillMatchedTrackTuple(const edm::RefToBase<Track>& trkRef,
			     Int_t& nmatches,vector<float>& result);

  void fillTrackTuple(const TrackingParticleRef& trkRef,
		      vector<float>& result);

  edm::RefToBase<Track> findRecoTrackMatch(const TrackingParticleRef& trk,
					   reco::SimToRecoCollection& reco2Sim,
					   Int_t& nmatches);
  TrackingParticleRef   findSimTrackMatch(const edm::RefToBase<Track>& trkRef,
					  reco::RecoToSimCollection& reco2Sim,
					  Int_t& nmatches);
 
  Int_t  getDetLayerId(const PSimHit& simHit);
  Int_t  getNumberOfPixelHits(const TrackingParticle& simTrack);
  Int_t  getNumberOfSimHits(const TrackingParticle& simTrack);
  Int_t  getSimParentId(const TrackingParticleRef& trk);
 
  
  
  Int_t  layerFromDetid(const DetId& detId); 
 
  void matchRecoMuons(edm::Handle<reco::MuonTrackLinksCollection>& muCollection,
		      reco::RecoToSimCollection& trk,reco::RecoToSimCollection& mu);
  void matchRecoDimuons(edm::Handle<reco::MuonTrackLinksCollection>& muCollection,
			reco::RecoToSimCollection& trk,	reco::RecoToSimCollection& mu);
  void matchSimDimuons(edm::Handle<TrackingParticleCollection>& simCollection,	
		       edm::Handle<reco::MuonTrackLinksCollection>& muCollection,
		       edm::Handle<edm::HepMCProduct>& hepEvt,	 
		       reco::SimToRecoCollection& trks, reco::SimToRecoCollection& mus); 
  void matchRecoTracks(edm::Handle<edm::View<Track> >& trackCollection,
		       reco::RecoToSimCollection& p);
  void matchSimTracks(edm::Handle<TrackingParticleCollection>& simCollection,
		      reco::SimToRecoCollection& q); 
  
  float refitWithVertex(const Track& recTrack);
  void dummVectorEntry(vector<float>& result, Int_t entries);

  // ----- member data -----
  bool                         doreco2sim;
  bool                         dosim2reco;
  bool                         matchmuons;
  bool                         matchtrackertracks;
  double                       matchhitfraction;
  double                       minpttrackertrkcut;
  Int_t                        proc, ntrk,nvtx;
  TNtuple                     *pnRecoSimDimuons;
  TNtuple                     *pnRecoSimMuons;
  TNtuple                     *pnRecoSimTracks;
  TNtuple                     *pnSimRecoDimuons;
  TNtuple                     *pnSimRecoMuons;
  TNtuple                     *pnSimRecoTracks;
  TNtuple                     *pnEventInfo;

  const TrackAssociatorByHits *theAssociatorByHits;
  const TrackerGeometry       *theTracker;
  const TransientTrackBuilder *theTTBuilder;
  const reco::BeamSpot        *theBeamSpot;
  const reco::VertexCollection *vertices;

  edm::InputTag                mutracktag;
  edm::InputTag                muonmaptag;
  edm::InputTag                simtracktag;
  edm::InputTag                statracktag; 
  edm::InputTag                trktracktag; 
  edm::ParameterSet            paramset;
  edm::Service<TFileService>   fs;
};


//___________________________________________________________________________
McMatchAnalyzer::McMatchAnalyzer(const edm::ParameterSet& pset):
  doreco2sim(pset.getUntrackedParameter<bool>("doReco2Sim") ),
  dosim2reco(pset.getUntrackedParameter<bool>("doSim2Reco") ),
  matchmuons(pset.getUntrackedParameter<bool>("matchMuons") ),
  matchtrackertracks(pset.getUntrackedParameter<bool>("matchTrackerTracks") ),
  matchhitfraction(pset.getParameter<double>("matchHitFraction") ),
  minpttrackertrkcut(pset.getParameter<double>("minPtTrackerTrkCut") ),
  mutracktag(pset.getUntrackedParameter<edm::InputTag>("mutracks") ),
  muonmaptag(pset.getParameter< edm::InputTag >("muonMapTag") ),
  simtracktag(pset.getUntrackedParameter<edm::InputTag>("simtracks") ),
  trktracktag(pset.getUntrackedParameter<edm::InputTag>("trktracks") ),
  paramset(pset)
{
  // constructor
 
}


//_____________________________________________________________________________
McMatchAnalyzer::~McMatchAnalyzer()
{
  // destructor

  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
}


//_____________________________________________________________________
void McMatchAnalyzer::analyze(const edm::Event& ev, const edm::EventSetup& es)
{
  // method called each event  
 
  edm::LogInfo("MCMatchAnalyzer")<<"Start analyzing each event ...";

  // Get generated event
  edm::Handle<edm::HepMCProduct> hepEv;
  ev.getByLabel("source",hepEv);
  proc = hepEv->GetEvent()->signal_process_id();
  edm::LogInfo("McMatchAnalyzer::analyze()")<<"Process ID= " << hepEv->GetEvent()->signal_process_id();

  // Get simulated particles
  edm::Handle<TrackingParticleCollection> simCollection;
  ev.getByLabel(simtracktag,simCollection);
  edm::LogInfo("McMatchAnalyzer::analyze()")<< "Size of simCollection = " << simCollection.product()->size();  

  // Get reconstructed tracks
  edm::Handle<edm::View<Track> >  trackCollection;
  ev.getByLabel(trktracktag, trackCollection); 
  edm::LogInfo("McMatchAnalyzer::analyze()")<< "Size of trackCollection = " << trackCollection.product()->size();
 
  //get beamSpot
  edm::Handle<reco::BeamSpot>      beamSpotHandle;
  ev.getByLabel("offlineBeamSpot", beamSpotHandle);
  theBeamSpot = beamSpotHandle.product();

   // Get vertices
  edm::Handle<reco::VertexCollection> vertexCollection;
  ev.getByLabel("pixelVertices",vertexCollection);
  vertices = vertexCollection.product();
  edm::LogInfo("McMatchAnalyzer::analyze()")<<"Number of vertices in the event = "<< vertexCollection.product()->size();

  // fill the event info TNtuple
  ntrk = trackCollection.product()->size();
  nvtx = vertexCollection.product()->size();
  {
    vector<float> result;
    result.push_back(proc); // proc
    result.push_back(ntrk); // ntrkr
    result.push_back(nvtx); // nvtxr
    
    pnEventInfo->Fill(&result[0]);
  }  

  // Associators
  //tracker trakcs
  reco::SimToRecoCollection trkSim2Reco;
  reco::RecoToSimCollection trkReco2Sim;
  if(dosim2reco)
    {
      trkSim2Reco = theAssociatorByHits->associateSimToReco(trackCollection,simCollection,&ev);
      edm::LogInfo("McMatchAnalyzer::analyze()")<<"Got track sim2reco association maps for tracks! ";
    }
  if(doreco2sim)
    {
      trkReco2Sim = theAssociatorByHits->associateRecoToSim(trackCollection,simCollection,&ev);
      edm::LogInfo("McMatchAnalyzer::analyze()")<<"Got track reco2sim association maps for tracks! ";
    }
  // Analyze
  if(matchtrackertracks)
    {
      if(dosim2reco) matchSimTracks(simCollection, trkSim2Reco);
      if(doreco2sim) matchRecoTracks(trackCollection, trkReco2Sim);
    }

  if(matchmuons)// dimu-h correlation
    {
      edm::LogInfo("McMatchAnalyzer::analyze()")<<"Using muons ...";

      // Get muon tracks only
      edm::Handle<reco::MuonTrackLinksCollection> muTrkLksCollection;
      ev.getByLabel(mutracktag,muTrkLksCollection);
      const reco::MuonTrackLinksCollection mutracks = *(muTrkLksCollection.product());
      edm::LogInfo("McMatchAnalyzer::analyze()") << "##### Size of MuonTrakLinksCollection = "<< mutracks.size()<<endl;

      // Get muon asscociation maps
      // reco2sim
      edm::Handle<reco::RecoToSimCollection> muRecoSimCollection;
      ev.getByLabel(muonmaptag,muRecoSimCollection);
      if (muRecoSimCollection.isValid() && doreco2sim) 
	{
	  edm::LogInfo("McMatchAnalyzer::analyze()") <<"Doing muRecoSim analysis!";
	  reco::RecoToSimCollection muReco2Sim = *(muRecoSimCollection.product());
	  edm::LogInfo("McMatchAnalyzer::analyze()") <<"##### Size of muReco2Sim collection = "<<muReco2Sim.size()<<endl;
	  
	  	  matchRecoMuons(muTrkLksCollection,trkReco2Sim,muReco2Sim);
		 
	  if( mutracks.size() >= 2 )
	    matchRecoDimuons(muTrkLksCollection,trkReco2Sim,muReco2Sim);
	  else 
	    LogDebug("McMatchAnalyzer::analyze()") <<"##### Size of muReco2Sim <2! No dimuons whatsoever";
	} 
      else 
	  LogDebug("McMatchAnalyzer::analyze()") << "##### NO RecoToSimCollection found for muons!";
        
      //sim2reco
	  edm::Handle<reco::SimToRecoCollection> muSimRecoCollection;
	  ev.getByLabel(muonmaptag,muSimRecoCollection);

	  if (muSimRecoCollection.isValid() && dosim2reco) 
	    {
	      edm::LogInfo("McMatchAnalyzer::analyze()") <<"Doing muSimReco analysis!";
	      reco::SimToRecoCollection muSim2Reco = *(muSimRecoCollection.product());
	      LogDebug("McMatchAnalyzer::analyze()") << "##### Size of muSim2Reco collection = "<<muSim2Reco.size();
	      matchSimDimuons(simCollection, muTrkLksCollection,hepEv,muSim2Reco, trkSim2Reco);
	    } 
	  else 
	    LogDebug("McMatchAnalyzer::analyze()") << "##### NO SimToRecoCollection found for muons!";
	  
	 // if dimu-h correlation 

  edm::LogInfo("McMatchAnalyzer::analyze()") << "[McMatchAnalyzer] done with event# " << ev.id();
  edm::LogInfo("McMatchAnalyzer::analyze()") << "----------------------------------------------" ;
    }
}

//____________________________________________________________________________
void McMatchAnalyzer::beginJob(const edm::EventSetup& es)
{
 // method called once each job just before starting event loop
  edm::LogInfo("MCMatchAnalyzer::beginJob()")<<"Begin job initialization ...";

  // Get tracker geometry
  edm::ESHandle<TrackerGeometry> tracker;
  es.get<TrackerDigiGeometryRecord>().get(tracker);
  theTracker = tracker.product();
 
  // Get track associator by hits
  edm::ESHandle<TrackAssociatorBase> theHitsAssociator;
  es.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",theHitsAssociator);
  theAssociatorByHits     = (const TrackAssociatorByHits*)theHitsAssociator.product();

  // Get transient track builder
  edm::ESHandle<TransientTrackBuilder> builder;
  es.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
  theTTBuilder = builder.product();
 
  // Bookkeepping
  // first global muon information, then sta muon, tracker muon info, and last simu info
  // the muon matching: tracker + sta hit-by-hit

  pnRecoSimMuons   = fs->make<TNtuple>("pnRecoSimMuons","pnRecoSimMuons",
				       "charge:chi2ndof:dxy:dxyerr:dz:dzerr:nvalidhits:eta:phi:pt:" // reco global muon
				       "chargesta:chi2ndofsta:dxy:dxyerr:dz:dzerr:nvalidhitssta:etasta:phista:ptsta:"// reco muon sta 
				       "nstamatch:" // #of reco sta matched to one sim sta
				       "idsimsta:idparentsimsta:chargesimsta:etasimsta:phisimsta:ptsimsta:"// matched sim sta
				       "chargetrk:chi2ndoftrk:dxy:dxyerr:dz:dzerr:nvalidhitstrk:etatrk:phitrk:pttrk:" //reco muon track
				       "ntrkmatch:" // #of reco trks amtched to one sim track
				       "idsimtrk:idparentsimtrk:chargesimtrk:etasimtrk:phisimtrk:ptsimtrk");// matched sim track
  pnRecoSimDimuons = fs->make<TNtuple>("pnRecoSimDimuons","pnRecoSimDimuons",
				       "eta:minv:phi:pt:" // dimuon 
				       "eta1:phi1:pt1:nstamatch1:idsimsta1:idparentsimsta1:ntrkmatch1:idsimtrk1:idparentsimtrk1:"//mu1
				       "eta2:phi2:pt2:nstamatch2:idsimsta2:idparentsimsta2:ntrkmatch2:idsimtrk2:idparentsimtrk2");//mu2
			  

  pnSimRecoDimuons = fs->make<TNtuple>("pnSimRecoDimuons","pnSimRecoDimuons",
				       "eta:minv:phi:pt:" // sim Z0
				       "charge1:eta1:phi1:pt1:" // sim muon 1
				       "nmatchsta1:chargestareco1:etastareco1:phistareco1:ptstareco1:" // reco sta  muon1
				       "nmatchtrk1:chargetrkreco1:etatrkreco1:phitrkreco1:pttrkreco1:" // reco trk  muon1  
				       "charge2:eta2:phi2:pt2:" // sim muon 2
				       "nmatchsta2:chargestareco2:etastareco2:phistareco2:ptstareco2:" // reco sta muon2
				       "nmatchtrk2:chargetrkreco2:etatrkreco2:phitrkreco2:pttrkreco2:" // reco trk muon2
				       "whomatch1:whomatch2:etareco1:phireco1:ptreco1:etareco2:phireco2:ptreco2:"
				       "etareco:minvreco:phireco:ptreco" ); // reconstructed Z0
				     

  pnRecoSimTracks  = fs->make<TNtuple>("pnRecoSimTracks","pnRecoSimTracks",
				       "charge:chi2ndof:dxy:dxyerr:dz:dzerr:nvalidhits:eta:phi:pt:" //reco track
				       "nmatch:"// #of reco trks matched to one sim track
				       "idsim:idparentsim:chargesim:etasim:phisim:ptsim"); // matched sim track

  pnSimRecoTracks  = fs->make<TNtuple>("pnSimRecoTracks","pnSimRecoTracks",
				       "id:idparent:charge:eta:phi:pt:npixelhits:"  //sim track
				       "nmatch:" //# of reco tracks matched to one sim track
				       "chargereco:etareco:phireco:ptreco"); // matched reco track
  pnEventInfo      = fs->make<TNtuple>("pnEventInfo","pnEventInfo","proc:ntrk:nvtx");
  
  edm::LogInfo("MCMatchAnalyzer::beginJob()")<<"Ended job initialization ...";

}


//_________________________________________________________________________
void McMatchAnalyzer::endJob()
{
  //method called once each job just after ending the event loop 

}


//______________________________________________
void McMatchAnalyzer::matchRecoDimuons(edm::Handle<reco::MuonTrackLinksCollection>& muTrkLksCollection,
				       reco::RecoToSimCollection& trkReco2Sim, 
				       reco::RecoToSimCollection& muReco2Sim )
{
  // fill the TNtuple pnRecoSimDimuons, 
  // with reco and corresponding sim tracks information

  edm::LogInfo("McMatchAnalyzer::matchRecoDimuons()")<<"Start reconstructing Dimuons!";

  
  for ( reco::MuonTrackLinksCollection::const_iterator links1 = muTrkLksCollection->begin(); 
	links1 != muTrkLksCollection->end()-1; links1++) 
    {
      edm::RefToBase<Track> glbRef1   = edm::RefToBase<Track>(links1->globalTrack() );
      if( glbRef1.isNull() ) continue;
   
      for ( reco::MuonTrackLinksCollection::const_iterator links2 = links1+1; links2 != muTrkLksCollection->end(); links2++) 
	{
	  vector<float> result;
	  
	  edm::RefToBase<Track> glbRef2   = edm::RefToBase<Track>(links2->globalTrack() );
	  if( glbRef2.isNull() ) continue;

	  if(glbRef1->charge()*glbRef2->charge()>0) continue; // just oppposite sign dimuons
	  
	  LogDebug("McMatchAnalyzer::matchRecoDimuons()")<<"Opposite sign dimuons!";
	  TLorentzVector child1, child2;
	  double en1 = sqrt(glbRef1->p()*glbRef1->p()+m_mu*m_mu);
	  double en2 = sqrt(glbRef2->p()*glbRef2->p()+m_mu*m_mu);
	  
	  child1.SetPxPyPzE(glbRef1->px(),glbRef1->py(),glbRef1->pz(),en1);
	  child2.SetPxPyPzE(glbRef2->px(),glbRef2->py(),glbRef2->pz(),en2);
	  
	  TLorentzVector dimuon;
	  dimuon = child1 + child2;
	 	  
	  result.push_back(dimuon.Eta());
	  result.push_back(dimuon.M());
	  result.push_back(dimuon.Phi());
	  result.push_back(dimuon.Pt());
	  
	  // ######### first muon  
	  result.push_back(glbRef1->eta());
	  result.push_back(glbRef1->phi());
	  result.push_back(glbRef1->pt());
	   // sta1  
	  edm::RefToBase<Track> staRef1   = edm::RefToBase<Track>(links1->standAloneTrack());  
	  if( !staRef1.isNull() )
	    fillRecoDimuonTuple(staRef1,muReco2Sim,result);
	  else dummVectorEntry(result,3);

	   // trk1
	  edm::RefToBase<Track> trkMuRef1 = edm::RefToBase<Track>(links1->trackerTrack());  
	  if( !trkMuRef1.isNull() )
	    fillRecoDimuonTuple(trkMuRef1,trkReco2Sim,result);
	  else dummVectorEntry(result,3);
	  LogDebug("McMatchAnalyzer::matchRecoDimuons()")<<"First sta,trk muon filled!";
	  //------------------------------------------------	 
	  // second muon
	  result.push_back(glbRef2->eta());
	  result.push_back(glbRef2->phi());
	  result.push_back(glbRef2->pt());
	  
	  // sta
	  edm::RefToBase<Track> staRef2   = edm::RefToBase<Track>(links2->standAloneTrack()); 
	  if( !staRef2.isNull() )
	    fillRecoDimuonTuple(staRef2,muReco2Sim,result);
	  else dummVectorEntry(result,3);
	  //trk2
	  edm::RefToBase<Track> trkMuRef2 = edm::RefToBase<Track>(links2->trackerTrack()); 
	  if( !trkMuRef2.isNull() )
	    fillRecoDimuonTuple(trkMuRef2,trkReco2Sim,result);
	  else dummVectorEntry(result,3);

	  LogDebug("McMatchAnalyzer::matchRecoDimuons()")<<"First sta,trk muon filled!";
	  pnRecoSimDimuons->Fill(&result[0]);

	  LogDebug("McMatchAnalyzer::matchRecoDimuons():result")<<"eta mu1 = "<<result[4]<<"\t eta mu2 = "<< result[13];
	  LogDebug("McMatchAnalyzer::matchRecoDimuons():direct")<<"eta mu1 = " << glbRef1->eta() <<"\t ta mu2 = " << glbRef2->eta();
	  result.clear();
       }// second muon loop
    }// first muon loop
  edm::LogInfo("McMatchAnalyzer::matchRecoDimuons()")<<"Done matchRecoDimuons()!";
}


//________________________________________________________________________
void McMatchAnalyzer::matchRecoMuons(edm::Handle<reco::MuonTrackLinksCollection>& muTrkLksCollection,
				     reco::RecoToSimCollection& trkReco2Sim, reco::RecoToSimCollection& muReco2Sim )
{
  // fill the TNtuple pnRecoSimMuons, 
  // with reco and corresponding sim tracks information
  edm::LogInfo("MCMatchAnalyzer::matchRecoMuons()")<<"Matching recoMuons ...";

  for ( reco::MuonTrackLinksCollection::const_iterator linkses = muTrkLksCollection->begin(); 
	linkses != muTrkLksCollection->end(); linkses++) 
    {    
      vector<float> result1;

      // global reco muon
      edm::RefToBase<Track> glbRef = edm::RefToBase<Track>(linkses->globalTrack() );
      if( glbRef.isNull() ) continue;
    
      LogDebug("McMatchAnalyzer::matchRecoMuons()")<<"##### Global track: charge "<<glbRef->charge()<<"\t pT= "<<glbRef->pt()<<endl;
      fillRecoTrackTuple(glbRef,result1);
 
      // stand alone (aka Sta) muon
      edm::RefToBase<Track> staRef = edm::RefToBase<Track>(linkses->standAloneTrack()); 
      if( !staRef.isNull() )
	{ 
	  LogDebug("McMatchAnalyzer::matchRecoMuons()")<< "##### STA track: charge = "<<staRef->charge()<<"\t pT= "<<staRef->pt()<<endl;
     
	  fillRecoTrackTuple(staRef,result1);

	  // ... matched sim
	  Int_t nStaSim = 0;
	  TrackingParticleRef matchedSimStaTrack = findSimTrackMatch(staRef,muReco2Sim,nStaSim);
	  
	  result1.push_back(nStaSim); // # of sim Sta matched to one reco sta  
	  Int_t ids      = -999;
	  Int_t parent = -999;
	  if( !matchedSimStaTrack.isNull() )
	    {
	      ids      = matchedSimStaTrack->pdgId();
	      parent = getSimParentId(matchedSimStaTrack);
	    }
	  result1.push_back(ids);
	  result1.push_back(parent);
	  fillTrackTuple(matchedSimStaTrack,result1);
	}
      else dummVectorEntry(result1,17);
      // tracker muon
      edm::RefToBase<Track> trkMuRef = edm::RefToBase<Track>(linkses->trackerTrack()); 
      if( !trkMuRef.isNull() ) 
	{      
	  LogDebug(" McMatchAnalyzer::matchRecoMuons()")<<"##### Track muon: charge = "<<trkMuRef->charge()<<"\t pT= "<<trkMuRef->pt()<<endl;
	  fillRecoTrackTuple(trkMuRef,result1);  
	  // ... matched sim 
	  Int_t nTrkSim = 0;
	  TrackingParticleRef matchedSimTrkTrack = findSimTrackMatch(trkMuRef,trkReco2Sim,nTrkSim);
	  result1.push_back(nTrkSim); // # of sim trks matched to one reco track
	  Int_t ids = -999;
	  Int_t parent = -999;
	  if( !matchedSimTrkTrack.isNull() )
	    {
	      ids      = matchedSimTrkTrack->pdgId();
	      parent   = getSimParentId(matchedSimTrkTrack);
	    }
	  result1.push_back(ids);
	  result1.push_back(parent);

	  fillTrackTuple(matchedSimTrkTrack,result1);
	}
      else dummVectorEntry(result1,17);

      // fill
      pnRecoSimMuons->Fill(&result1[0]);
      result1.clear();
    }// muon loop

  edm::LogInfo("MCMatchAnalyzer::matchRecoMuons()")<<"Done matchRecoMuons()!";
}


//________________________________________________________________________
void McMatchAnalyzer::matchRecoTracks(edm::Handle<edm::View<Track> >& trackCollection,
				      reco::RecoToSimCollection& p)
{
  // fill the TNtuple pnRecoSimTracks, 
  // with reco tracks and their corresponding sim tracks information
  
  edm::LogInfo("MCMatchAnalyzer::matchRecoTracks()")<<"Start matching reco tracks ...";

  for(edm::View<Track> ::size_type i=0; i < trackCollection.product()->size(); ++i)
  {
    edm::RefToBase<Track> recTrack(trackCollection, i);
   
    vector<float> result2;
    Int_t ids      = -999;
    Int_t parentId = -999;
    Int_t nSim = 0;
    // reco tracker
    if( !recTrack.isNull() )
      {
	fillRecoTrackTuple(recTrack,result2);
	TrackingParticleRef matchedSimTrack = findSimTrackMatch(recTrack,p,nSim);
	result2.push_back(nSim); // # matched sim track to one reco track
	if( !matchedSimTrack.isNull() )
	  {
	    //	result2.push_back(refitWithVertex(*recTrack));                 // ptv
	    // matched sim 
	    
	    ids      = matchedSimTrack->pdgId();
	    parentId = getSimParentId(matchedSimTrack);
	    //===========
	    //  if( !(matchedSimTrack->parentVertex().isNull()) )
// 	      {
// 		if(matchedSimTrack->parentVertex()->nSourceTracks() == 0)
// 		  parentId = 0;
// 		else
// 		  parentId = (*(matchedSimTrack->parentVertex()->sourceTracks_begin()))->pdgId();
// 	      }
// 	    else parentId = 999;
	     //==============
	  
	  }// there is a matched track
	result2.push_back(ids);
	result2.push_back(parentId);
	fillTrackTuple(matchedSimTrack,result2);
	
	pnRecoSimTracks->Fill(&result2[0]);
	result2.clear();
      }//valid tracker track 
  
  }//i trackCollection
  edm::LogInfo("MCMatchAnalyzer::matchRecoTracks()")<<"Done matchRecoTracks().";
}


//_________________________________________________________________________
void McMatchAnalyzer::matchSimDimuons(edm::Handle<TrackingParticleCollection>& simCollection,
				      edm::Handle<reco::MuonTrackLinksCollection>& muTrkLksCollection,
				      edm::Handle<edm::HepMCProduct>& evt,
				      reco::SimToRecoCollection& staSim2Reco, reco::SimToRecoCollection& trkSim2Reco)
{
  // fill TNtuple pnSimRecoDimuons with gen info for the gen Z0 and muon daughters and reco info for the reco muons

  // status == 1 : stable particle at pythia level, but which can be decayed by geant
  // status == 3 : outgoing partons from hard scattering

  // The barcode is the particle's reference number, every vertex in the
  // event has a unique barcode. Particle barcodes are positive numbers,
  // vertex barcodes are negative numbers.

 // ("pnSimRecoDimuons","pnSimRecoDimuons",
//   "eta:minv:phi:pt:" //parent
//   "charge1:eta1:phi1:pt1:" // sim muon 1
//   "nmatchsta1:chargestareco1:etastareco1:phistareco1:ptstareco1:" // reco sta  muon1
//   "nmatchtrk1:chargetrkreco1:etatrkreco1:phitrkreco1:pttrkreco1:" // reco trk  muon1  
//   "charge2:eta2:phi2:pt2:" // sim muon 2
//   "nmatchsta2:chargestareco2:etastareco2:phistareco2:ptstareco2:" // reco sta muon2
//   "nmatchtrk2:chargetrkreco2:etatrkreco2:phitrkreco2:pttrkreco2"); // reco trk muon2
//   "whomatch1:whomatch2:ptreco1:phireco1:etareco1:ptreco2:phireco2:etareco2"); // reco muons
//   "etareco:minvreco:phireco:ptreco:" // reconstructed Z0


  edm::LogInfo("MCMatchAnalyzer::matchSimDimuons()")<<"Start matching dimuons ...";
  
  const HepMC::GenEvent * myGenEvent = evt->GetEvent();

  Bool_t foundmuons = false;
  vector<Int_t> ixmu;
  vector<float> result3;
  for ( HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p ) 
    { 
      if ( ( (*p)->pdg_id() == 23 ) && (*p)->status() == 3 ) 
	{ 
	  for( HepMC::GenVertex::particle_iterator aDaughter=(*p)->end_vertex()->particles_begin(HepMC::descendants); 
	       aDaughter !=(*p)->end_vertex()->particles_end(HepMC::descendants);
	       aDaughter++)
	    {
	      if ( abs((*aDaughter)->pdg_id())==13 && (*aDaughter)->status()==1 )
		    
		{
		  ixmu.push_back((*aDaughter)->barcode());
		  LogDebug("MCMatchAnalyzer::matchSimDimuons()") << "Stable muon from Z0" << "\tindex= "<< (*aDaughter)->barcode(); 
		  foundmuons = true;
		}//  z0 descendentes
	    }
	  if (foundmuons)
	    {
	      result3.push_back((*p)->momentum().eta());
	      result3.push_back((*p)->momentum().m());
	      result3.push_back((*p)->momentum().phi());
	      result3.push_back((*p)->momentum().perp());
	    }
	}// Z0
    }// loop over genParticles

  if( foundmuons )
    {
      vector<edm::RefToBase<Track> > stastub;
      vector<edm::RefToBase<Track> > trkstub;
      Int_t ndaughters  = 0;

      for(TrackingParticleCollection::size_type i=0; i < simCollection.product()->size(); i++)
	{
	  const TrackingParticleRef simTrk(simCollection, i);
	  if( simTrk.isNull() ) continue;
	  const SimTrack *gTrack = &(*simTrk->g4Track_begin());
	  if( gTrack == NULL ) continue;

	  if(abs(gTrack->type())==13 && (gTrack->genpartIndex()==ixmu[0] || gTrack->genpartIndex()==ixmu[1] ))
	    {
	    
	      // sim track
	      fillTrackTuple(simTrk,result3);
	     
	      // staMatched_muon
	      Int_t nSta = 0;      
	      edm::RefToBase<Track>  matchedRecSta = findRecoTrackMatch(simTrk,staSim2Reco,nSta);
	      result3.push_back(nSta); // # of reco tracks matched to one sim track

	      if( !matchedRecSta.isNull() )
		{
		  fillMatchedTrackTuple( matchedRecSta,nSta,result3);
		}
	      else 
		dummVectorEntry(result3,4);
		
	      // trackerMatched_muon		 
	      // rec
	      Int_t nTrk = 0;
	      edm::RefToBase<Track> matchedRecTrk = findRecoTrackMatch(simTrk,trkSim2Reco,nTrk);
	      result3.push_back(nTrk); // # of reco tracks matched to one sim track

	      if( !matchedRecTrk.isNull() )
		{
		  fillMatchedTrackTuple( matchedRecTrk,nTrk,result3);  
		}
	      else 
		dummVectorEntry(result3,4);

	      if( !matchedRecSta.isNull() || !matchedRecTrk.isNull() )
		{
		  stastub.push_back(matchedRecSta);
		  trkstub.push_back(matchedRecTrk);
		  ndaughters++;
		}
	    }// muon id
	  
	}// TrackCollection loop

      if(ndaughters == 1) // not found second daughter
	{
	  LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Only one sim daughter found to match a reco track"<<endl;
	  dummVectorEntry(result3,14); 

	}
      if(ndaughters == 0) // no daughter found
	{
	  //	  LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<
	  //	  cout<<"No sim daughter reconstructed"<<endl;
	  dummVectorEntry(result3,28);
	}

      // get the reco Z0
      if(ndaughters == 2)
	{
	  LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Both sim daughters found to match a reco track"<<endl;

	  Int_t nsplit = 0;
	  vector<edm::RefToBase<Track> > glbchild;
	  for ( reco::MuonTrackLinksCollection::const_iterator links = muTrkLksCollection->begin(); 
		links != muTrkLksCollection->end(); links++) 
	    {
	      for( Int_t istubs = 0; istubs < ndaughters; istubs++)
		{
		  bool stamatch = !stastub[istubs].isNull() && !links->standAloneTrack().isNull() &&
		    (links->standAloneTrack()->pt()  ==  stastub[istubs]->pt() &&
		     links->standAloneTrack()->phi() ==  stastub[istubs]->phi() &&
		     links->standAloneTrack()->eta() ==  stastub[istubs]->eta()
		     ); 
			    
		  bool trkmatch = !links->trackerTrack().isNull() && !trkstub[istubs].isNull()  &&
		    (links->trackerTrack()->pt()    == trkstub[istubs]->pt()
		     && links->trackerTrack()->phi() == trkstub[istubs]->phi() 
		     && links->trackerTrack()->eta() == trkstub[istubs]->eta() 
		     );
		
		  if(!stamatch && !trkmatch)
		    {
		      LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Sim matched to none of the links!!."<<endl;
		      continue;
		    }

		  if(stamatch && trkmatch)
		    { 
		      LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Sim matched to both sta and trk."<<endl;
		      result3.push_back(1.); // both ssta and trk match == whomatch
		    }
		  if(stamatch && !trkmatch)
		    {
			  
		      LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Sim matched only to sta."<<endl;
		      result3.push_back(2.); // only sta match
		    }
		  if(!stamatch && trkmatch)
		    {
		      LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Sim matched only to trk."<<endl;
		      result3.push_back(3.); // only sta match
		    }

		  edm::RefToBase<Track> glbmu = edm::RefToBase<Track>(links->globalTrack() );
		  glbchild.push_back(glbmu);
		  nsplit++;
		  if(nsplit==1) break;
		}// each daughter
	      
	    }//each muontracklink

	  if(glbchild.size()==2)
	    {
	  
	      edm::RefToBase<Track> ch1 = glbchild[0];
	      edm::RefToBase<Track> ch2 = glbchild[1];
	      
	      // fill the individual reco global muon for the whole reconstructed muons
	      result3.push_back(ch1->eta());
	      result3.push_back(ch1->phi());
	      result3.push_back(ch1->pt());
	      result3.push_back(ch2->eta());
	      result3.push_back(ch2->phi());
	      result3.push_back(ch2->pt());
	      
	      // reconstructed Z0
	      TLorentzVector child1, child2;
	      double en1 = sqrt(ch1->p()*ch1->p()+m_mu*m_mu); 
	      double en2 = sqrt(ch2->p()*ch2->p()+m_mu*m_mu);
	      
	      child1.SetPxPyPzE(ch1->px(),ch1->py(),ch1->pz(),en1);
	      child2.SetPxPyPzE(ch2->px(),ch2->py(),ch2->pz(),en2);
	      
	      
	      TLorentzVector recoz0;
	      recoz0 = child1 + child2;
	      result3.push_back(recoz0.Eta());
	      result3.push_back(recoz0.M());
	      result3.push_back(recoz0.Phi());
	      result3.push_back(recoz0.Pt());
	      
	      LogDebug("MCMatchAnalyzer::matchSimDimuons()")<<"Filled the matched reco Z0 info";
	      glbchild.clear();
	      stastub.clear();
	      trkstub.clear();
	    
	    }
	  else 
	    {
	      if(glbchild.size()==1) dummVectorEntry(result3,11);
	      else dummVectorEntry(result3,12);
	    }
	}
      else // no recoZ0
	{
	  dummVectorEntry(result3,12); 
	}
      // 2 (which part of the single muon was matched)+4(recoZ0)+6(individual reco muons)
  
      pnSimRecoDimuons->Fill(&result3[0]);
      result3.clear();

    }// there are sim muons from Z0
 
  edm::LogInfo("MCMatchAnalyzer::matchSimDimuons()")<<"Done matchSimDimuons()!.";
}


//_________________________________________________________________________
void McMatchAnalyzer::matchSimTracks(edm::Handle<TrackingParticleCollection>& simCollection,
				     reco::SimToRecoCollection& q)
{
  // fill TNtuple pnSimRecoTracks with sim track and corresponding reco tracks info
  edm::LogInfo("MCMatchAnalyzer::matchSimTracks()")<<"Start matching simTracks ...";
 
  for(TrackingParticleCollection::size_type i=0; i < simCollection.product()->size(); i++)
    {
      const TrackingParticleRef simTrack(simCollection,i);
      if( simTrack.isNull() ) continue;
      if( simTrack->charge() == 0 || simTrack->pt() < minpttrackertrkcut ) continue;
       
      Int_t nPixelLayers = getNumberOfPixelHits(*simTrack);
      if( nPixelLayers < 3 ) continue;

      vector<float> result4;
      
      // sim track
    
      result4.push_back(simTrack->pdgId());               // id
      Int_t parent = getSimParentId(simTrack);
      
      //============
     //  if(simTrack->parentVertex()->nSourceTracks() > 0)
// 	{
// 	  // track is not primary, has a parent
// 	  TrackingVertex::tp_iterator iv = simTrack->parentVertex()->sourceTracks_begin();
// 	  result4.push_back((*iv)->pdgId());
// 	}
//       else result4.push_back(0) ; // parent id
      //==============

      result4.push_back(parent) ;
      fillTrackTuple(simTrack,result4);
      result4.push_back(nPixelLayers); // nhits
      
      // recoTracker matched 
      Int_t nRec = 0;
      edm::RefToBase<Track> matchedRecTrack = findRecoTrackMatch(simTrack,q,nRec);
      result4.push_back(nRec); // # of reco tracks matched to one sim track
      if( !matchedRecTrack.isNull() )
	  fillMatchedTrackTuple(matchedRecTrack,nRec,result4);
      else 
	  dummVectorEntry(result4,4);
      
      // fill
      pnSimRecoTracks->Fill(&result4[0]);
      result4.clear();
    }

  edm::LogInfo("MCMatchAnalyzer::matchSimTracks()")<<"Done matchSimTracks()!";
}


//_______________________________________________________________________
void McMatchAnalyzer::fillRecoDimuonTuple(const edm::RefToBase<Track>& trkRef,
					  reco::RecoToSimCollection& reco2Sim,
					  vector<float>& result)
{
  // get the match, 
  // fill the result with the id of match, id of parent, and and number of matches
  edm::LogInfo("MCMatchAnalyzer::fillRecoDimuonTuple()")<<"Filling dimuon tuple ...";

  Int_t nTrkSim = 0;
  TrackingParticleRef matchedSim = findSimTrackMatch(trkRef,reco2Sim,nTrkSim); 
  Int_t ids = -999;
  Int_t parentID = -999;

  if( !matchedSim.isNull() )
    {
      ids = matchedSim->pdgId();
      if(nTrkSim > 0)
	parentID = getSimParentId(matchedSim);
     
    }
 
  result.push_back(nTrkSim);
  result.push_back(ids);
  result.push_back(parentID);
  
  edm::LogInfo("MCMatchAnalyzer::fillRecoDimuonTuple()")<<"Done fillRecoDimuonTuple()!";
 }


//_______________________________________________________________________
void McMatchAnalyzer::fillRecoTrackTuple(const edm::RefToBase<Track>& trkRef, vector<float>& result)
{
  // fill reco track info to be fed later to the tuple
  edm::LogInfo("MCMatchAnalyzer::fillRecoTrackTuple()")<<"Filling reco track info ...";
  // transverse and z DCA calculated with respect to the beamSpot plus the error

  if( !trkRef.isNull())
    {
      double dt      = trkRef->dxy(theBeamSpot->position());
      double sigmaDt = sqrt(trkRef->dxyError() * trkRef->dxyError() +
			    theBeamSpot->BeamWidth() * theBeamSpot->BeamWidth());
      
      double dz      = trkRef->dz(theBeamSpot->position());
      double sigmaDz = sqrt(trkRef->dzError() * trkRef->dzError() +
			    theBeamSpot->BeamWidth() * theBeamSpot->BeamWidth());
      
      result.push_back(trkRef->charge());                         // charge
      result.push_back(trkRef->normalizedChi2());                 // normalized Chi2
      result.push_back(dt);                                       // transverse DCA
      result.push_back(sigmaDt);                                  // Sigma_Dca_t
      result.push_back(dz);                                       // z DCA
      result.push_back(sigmaDz);                                  // Sigma_Dca_z
      result.push_back(trkRef->numberOfValidHits());              // number of hits found 
      result.push_back(trkRef->eta());                            // eta
      result.push_back(trkRef->phi());                            // phi
      result.push_back(trkRef->pt());                             // pt
    }
  else dummVectorEntry(result,10);
  edm::LogInfo("MCMatchAnalyzer::fillRecoTrackTuple()")<<"Done fillRecoTrackTuple()!";
}


//____________________________________________________________________
void McMatchAnalyzer::fillMatchedTrackTuple(const edm::RefToBase<Track>& recoTrack,
					Int_t& nMatches, vector<float>& result)
{
  // the info part of the matched sim track
  edm::LogInfo("MCMatchAnalyzer::fillMatchedTrackTuple()")<<"Filling matched reco track info ...";

  if(nMatches > 0 )
    {
      result.push_back(recoTrack->charge()); // chargereco
      result.push_back(recoTrack->eta());    // eta
      result.push_back(recoTrack->phi());    // phireco
      result.push_back(recoTrack->pt());     // ptreco
    }
  else dummVectorEntry(result,4);
  edm::LogInfo("MCMatchAnalyzer::fillMatchedTrackTuple()")<<"Done fillMatchedTrackTuple()!"; 
}


//____________________________________________________________________
void McMatchAnalyzer::fillTrackTuple(const TrackingParticleRef& recoTrack,
				     vector<float>& result)
{
  // the info part of the matched sim track
  edm::LogInfo("MCMatchAnalyzer::fillMatchedTrackTuple()")<<"Filling matched reco track info ...";

  if( !recoTrack.isNull() )
    {
      result.push_back(recoTrack->charge()); // chargereco
      result.push_back(recoTrack->eta());    // eta
      result.push_back(recoTrack->phi());    // phireco
      result.push_back(recoTrack->pt());     // ptreco
    }
  else dummVectorEntry(result,4);
  edm::LogInfo("MCMatchAnalyzer::fillMatchedTrackTuple()")<<"Done fillMatchedTrackTuple()!"; 
}


//________________________________________________________________________
edm::RefToBase<Track> McMatchAnalyzer::findRecoTrackMatch(const TrackingParticleRef& simTrack,
							  reco::SimToRecoCollection& sim2Reco, Int_t& nMatches)
{
  // return the matched reco track
  edm::LogInfo("MCMatchAnalyzer::findRecoTrackMatch()")<<"Finding reco track match ...";

  edm::RefToBase<Track> recoTrackMatch;

  if(sim2Reco.find(simTrack) != sim2Reco.end()) 	 
    {
      vector<pair<edm::RefToBase<Track>, double> > recTracks = sim2Reco[simTrack];
      for(vector<pair<edm::RefToBase<Track>,double> >::const_iterator it = recTracks.begin(); 
	  it != recTracks.end(); it++)
	{
	  edm::RefToBase<Track> recTrack = it->first;
	  float fraction = it->second;
	  
	  if(fraction > matchhitfraction )
	    { 
	      recoTrackMatch = recTrack; 
	      nMatches++; 
	    }
	}
    }
  edm::LogInfo("MCMatchAnalyzer::findRecoTrackMatch()")<<"Done findRecoTrackMatch().";
  return recoTrackMatch;
}


//________________________________________________________________________
TrackingParticleRef McMatchAnalyzer::findSimTrackMatch(const edm::RefToBase<Track>& recoTrack,
						       reco::RecoToSimCollection& reco2Sim, Int_t& nMatches)
{
  // return the matched sim track
  edm::LogInfo("MCMatchAnalyzer::findSimTrackMatch()")<<"Finding sim track match ...";

  TrackingParticleRef simTrackMatch;

  if(reco2Sim.find(recoTrack) != reco2Sim.end())  
    {
      LogDebug("McMatchAnalyzer::matchRecoMuons()")<<"Found Track match !!!!";
      vector<pair<TrackingParticleRef, double> > simTracks = reco2Sim[recoTrack];
      
      for(vector<pair<TrackingParticleRef, double> >::const_iterator it = simTracks.begin(); it != simTracks.end(); ++it)
	{
	  
	  TrackingParticleRef simTrack = it->first;
	  float fraction = it->second;
	  
	  if(fraction > matchhitfraction)
	    { 
	      simTrackMatch = simTrack; 
	      nMatches++; 
	    }
	}
    }

  edm::LogInfo("MCMatchAnalyzer::findSimTrackMatch()")<<"Done finding sim track match!";

  return simTrackMatch;
}


//_________________________________________________________________________
Int_t McMatchAnalyzer::getSimParentId(const TrackingParticleRef& match)
{
  // return the particle ID of the matched sim track
 // genParticle corresponding to the matched TrackingParticle
  // it is a final, stable particle (status=1): after final state radiaiton

  // same particle before final state radiation (status=3); from this one get to the real parent;
  // this is the 'parent' of the final state particle, p; same pdg_id

  // eg: look for the parent in the chain pion->mu+x
  // first parent, will be still the muon, before the final state radiation: pion->mu*+x->mu+x
  // from this parent get the parent that produced the muon, the pion in this case

  edm::LogInfo("McMatchAnalyzer::getSimParentId")<<"Getting the parent Id ...";
  
  Int_t parentId = -999; // primary vertex, default
  
  if( match.isNull() ) return parentId;
   
  for (TrackingParticle::genp_iterator b = match->genParticle_begin();
       b != match->genParticle_end(); ++b)
    {
      
      const HepMC::GenParticle *mother = (*b)->production_vertex() && 
	(*b)->production_vertex()->particles_in_const_begin() != 
	(*b)->production_vertex()->particles_in_const_end() ?
	*((*b)->production_vertex()->particles_in_const_begin()):0;
      
      if( mother!=0 && 
	  (abs(mother->pdg_id())<10 || (mother->pdg_id()<=100 && mother->pdg_id()>=80)) ) return 0; 
     
      if( mother!=0 && !(isnan(mother->momentum().mag())) && !isnan(abs(mother->pdg_id())) &&
	  (mother->momentum().mag() < 1e+5) && (mother->momentum().mag() > 1e-5) && (abs(mother->pdg_id())<1e+6)
	  )
	{
	
	  parentId = mother->pdg_id();
	  LogDebug("McMatchAnalyzer::getSimParentId")<<" 1 parentId = "<<parentId<<"\t parent_pt= "<<mother->momentum().perp();
	  LogDebug("McMatchAnalyzer::getSimParentId")<<" 1 kidId = "<<match->pdgId()<<"\t kid_pt= "<<match->pt(); 
	 
	  // real parent after final state radiation
	  const HepMC::GenParticle *motherMother = mother->production_vertex() &&
	    mother->production_vertex()->particles_in_const_begin() !=
	    mother->production_vertex()->particles_in_const_end() ?
	    *(mother->production_vertex()->particles_in_const_begin()) : 0 ;
	 
	  if( motherMother!=0 && 
	      (abs(motherMother->pdg_id())<10 || (motherMother->pdg_id()<=100 && motherMother->pdg_id()>=80)) ) return parentId;
	  
	  if( motherMother!=0 
	      && !(isnan(motherMother->momentum().perp()) ) && !(isnan(motherMother->pdg_id()) ) 
	      && (motherMother->momentum().perp() < 1e+5) && (motherMother->momentum().perp() > 1e-5) && (motherMother->pdg_id() < 1e+6)
	      ) 
	    {
	      parentId = motherMother->pdg_id();
	      LogDebug("McMatchAnalyzer::getSimParentId")<<" 2 parentId = "<<parentId<<"\t parent_pt= "<<motherMother->momentum().perp();
	    }
	  motherMother = 0;
	}// valid values for first particle in the vtx
      
      mother = 0;
    }
  // LogDebug("McMatchAnalyzer::getSimParentId ()")
  
  edm::LogInfo("McMatchAnalyzer::getSimParentId")<<"Done with getSimParentId!\t"<<parentId;
  return parentId;
}


//________________________________________________________________________
Int_t McMatchAnalyzer::getDetLayerId(const PSimHit& simHit)
{
  // return the detector layer ID
  LogDebug("MCMatchAnalyzer::getDetLayerId()")<<"Get detector layer ID ...";
 
  Int_t layerId;
  unsigned int id = simHit.detUnitId();

  if(theTracker->idToDet(id) && theTracker->idToDetUnit(id)->subDetector() == GeomDetEnumerators::PixelBarrel)
  {
    PXBDetId pid(id);
    layerId = pid.layer() - 1;
  }
  else
  {
    PXFDetId pid(id);
    layerId = 2 + pid.disk();
  }
 
  LogDebug("MCMatchAnalyzer::getDetLayerId()")<<"Done";
  return layerId;
}


//___________________________________________________________________________
Int_t McMatchAnalyzer::getNumberOfPixelHits(const TrackingParticle& simTrack)
{
  // return number of pixel hits
  // pixel traker: 3 barel layers + 2 endcap disks
  edm::LogInfo("MCMatchAnalyzer::getNumberOfPixelHits()")<<"Get number of pixel hits ...";
 
  const Int_t nLayers = 5;
  vector<bool> filled(nLayers,false);
  for(std::vector<PSimHit>::const_iterator simHit = simTrack.pSimHit_begin();
                                           simHit!= simTrack.pSimHit_end();
                                           simHit++)
    {    
      unsigned int id = simHit->detUnitId();
      //  cout<<"id = "<< id<<endl;
     
      if( (theTracker->idToDet(id) && theTracker->idToDet(id)->subDetector() == GeomDetEnumerators::PixelBarrel) ||
	  (theTracker->idToDet(id) && theTracker->idToDet(id)->subDetector() == GeomDetEnumerators::PixelEndcap) )
	{
	  filled[getDetLayerId(*simHit)] = true;
	}
    }
  
  // Count the number of filled pixel layers
  Int_t fLayers = 0;
  for(Int_t i=0; i<nLayers; i++)
    {
      if(filled[i] == true) fLayers++;
    }
 
  LogDebug("MCMatchAnalyzer::getNumberOfPixelHits()")<<"Done!";

  return fLayers;
}


//______________________________________________________________________________
float McMatchAnalyzer::refitWithVertex(const Track & recTrack)
{ 
  // returns the momentum of the track updated at th primary vertex
  edm::LogInfo("MCMatchAnalyzer::refitWithVertex()")<<"Refit with primary vertex ...";
 
  TransientTrack theTransientTrack = theTTBuilder->build(recTrack);
  
  // If there are vertices found
  if(vertices->size() > 0)
  { 
    float dzmin = -1.; 
    const reco::Vertex * closestVertex = 0;

    // Look for the closest vertex in z
    for(reco::VertexCollection::const_iterator
        vertex = vertices->begin(); vertex!= vertices->end(); vertex++)
    {
      float dz = fabs(recTrack.vertex().z() - vertex->position().z());
      if(vertex == vertices->begin() || dz < dzmin)
      { dzmin = dz ; closestVertex = &(*vertex); }
    }

    // Get vertex position and error matrix
    GlobalPoint vertexPosition(closestVertex->position().x(),
                               closestVertex->position().y(),
                               closestVertex->position().z());

    float beamSize = 15e-4; // 15 um
    GlobalError vertexError(beamSize*beamSize, 0,
                            beamSize*beamSize, 0,
                            0,closestVertex->covariance(2,2));

    // Refit track with vertex constraint
    SingleTrackVertexConstraint stvc;
    pair<TransientTrack, float> result =
      stvc.constrain(theTransientTrack, vertexPosition, vertexError);

    edm::LogInfo("MCMatchAnalyzer::refitWithVertex()")<<"Done";
    return result.first.impactPointTSCP().pt();
  }
  else
    {
      edm::LogInfo("MCMatchAnalyzer::refitWithVertex()")<<"Done";
      return recTrack.pt();
    }
}


//_______________________________________________________________________
void McMatchAnalyzer::dummVectorEntry(vector<float>& result, Int_t nEntries)
{
  // add to vector nentries of '-999' 
  for(Int_t i = 1; i<= nEntries; i++)
    result.push_back(-999);
}


//________________________________________________________________________

DEFINE_FWK_MODULE(McMatchAnalyzer);