//
//
// -*- C++ -*-
//
// Package:    Muons
// Class:      Muons
// 
/**\class Muons Muons.cc Muons/Muons/src/Muons.cc

 Description: <one line class summary>

 Implementation:
     <Notes on implementation>
*/
//
// Original Author:  Jorge Robles
//         Created:  Wed Aug 20 22:50:37 CEST 2008
// $Id$
//
//


// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

// system include files
#include <iostream>
#include <cmath>
#include <fstream>
#include <stdio.h>
#include <math.h>
#include <vector>


//Analysis Files
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include <DataFormats/TrackReco/interface/Track.h>
#include <DataFormats/HepMCCandidate/interface/GenParticleCandidate.h>
//Root
#include <TFile.h>
#include <TFolder.h>
#include <TH1.h>
#include <TH2.h>
#include <TDirectory.h>
#include <TMath.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include <TF1.h>
#include <TH2F.h>
#include <TString.h>


//Particle.h includes
#include "DataFormats/Math/interface/Point3D.h"
#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/Common/interface/BoolCache.h"
#include "DataFormats/Candidate/interface/Particle.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include <DataFormats/MuonReco/interface/Muon.h>
#include <DataFormats/Candidate/interface/LeafCandidate.h>
#include <DataFormats/RecoCandidate/interface/RecoCandidate.h>
#include <DataFormats/MuonReco/interface/MuonEnergy.h>

//GenParticle includes



//
// class decleration
//

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




class Muons : public edm::EDAnalyzer {
   public:
      explicit Muons(const edm::ParameterSet&);
      ~Muons();
     
     TH1 *gMuonEta_h;
     TH1 *gMuonPt_h;
     TH1 *gZ0_InvMass_h;
     TH1 *gZ0_Eta_h;
     TH1 *gZ0_Rap_h;
     TH1 *gZ0_Pt_h;
     TH2 *gZ0_RapVSPt_h;
     TH2 *gZ0_EtaVSPhi_h;

     TH1 *rMuonEta_h;
     TH1 *rMuonPt_h;
     TH1 *rZ0_InvMass_h;
     TH1 *rZ0_Eta_h;
     TH1 *rZ0_Rap_h;
     TH1 *rZ0_Pt_h;
     TH2 *rZ0_RapVSPt_h;
     TH2 *rZ0_EtaVSPhi_h;

     TH1 *EffMuonEta_h;
     TH1 *EffMuonPt_h;
     TH1 *EffZ0_InvMass_h;
     TH1 *EffZ0_Eta_h;
     TH1 *EffZ0_Rap_h;
     TH1 *EffZ0_Pt_h;

     TH1 *NChambers_h;
     TH1 *NSegments_h;
     TH1 *CalEnergy_h;
    // TH1 *MuEnergy_h;
 


   private:
      virtual void beginJob(const edm::EventSetup&) ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      // ----------member data ---------------------------
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
Muons::Muons(const edm::ParameterSet& iConfig)

{
   //now do what ever initialization is needed

}


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

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
Muons::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;
   using namespace reco;

//     HepMC
//     ->GenParticle
edm::Handle<edm::HepMCProduct>gMuons;
iEvent.getByLabel("source","",gMuons);

//for documentation of all this schwagg see
//http://projects.hepforge.org/rivet/code/hepmc/namespaceHepMC.html

int counter=0;
for (HepMC::GenEvent::particle_const_iterator gMu1 = gMuons->GetEvent()->particles_begin(); gMu1 !=gMuons->GetEvent()->particles_end();++gMu1)
//for (reco::GenParticle::const_iterator gMu1 = gMuons->GetEvent()->particles_begin(); gMu1 !=gMuons->GetEvent()->particles_end();++gMu1)
	{	
	HepMC::GenParticle* ith_gMu1 = *gMu1;
	int id1 = ith_gMu1->pdg_id();
	//Check that the generated particle is a Muon
	if(abs( id1 )==13)
		{
		counter++;
		HepMC::GenVertex* prodVtx1 = ith_gMu1->production_vertex();
		//Loop over te incoming particles to this vertex and look for a Z0
		for(HepMC::GenVertex::particles_in_const_iterator inPart1 = prodVtx1->particles_in_const_begin(); inPart1 != prodVtx1->particles_in_const_end();inPart1++)
			{
			HepMC::GenParticle* Parent1 = *inPart1;
			int Parent_id1 = Parent1->pdg_id();
			//cout<<"Mu's parent:  "<<Parent_id1<<endl;
			if(Parent_id1 == 23)
				{
				double massMu = 0.105658369;	
				double gMuonPx1  =ith_gMu1->momentum().px();
				double gMuonPy1  =ith_gMu1->momentum().py();
				double gMuonPz1  =ith_gMu1->momentum().pz();
				double gMuonE1   =ith_gMu1->momentum().e();
				double gMuonPt1  =sqrt(gMuonPx1*gMuonPx1+gMuonPy1*gMuonPy1);
				gMuonPt_h->Fill(gMuonPt1);
				gMuonEta_h->Fill(ith_gMu1->momentum().eta());
				//double gMuonEnergy1 = sqrt((gMuonMag1*gMuonMag1)+(massMu*massMu));
				//cout<<"The energy .e() is "<<gMuonE1 <<endl;
				//cout<<"Thecalculated energy is "<<gMuonEnergy1 <<endl;
				//both of the above are the same, just checking...
				TLorentzVector gMuonP1(gMuonPx1,gMuonPy1,gMuonPz1,gMuonE1);
				//Once you found one Muon loop over the other particles
				int flag=0; 
				for(HepMC::GenEvent::particle_const_iterator gMu2 = gMu1;gMu2 != gMuons->GetEvent()->particles_end();)
					{
					gMu2++;
					//This is because the regular loop was giving some problems, this seems to work
					if(gMu2==gMuons->GetEvent()->particles_end()) continue;
					HepMC::GenParticle* ith_gMu2 = *gMu2;
					int id2 = ith_gMu2->pdg_id();
					//Make sure the second particle is also a muon
					if(abs(id2)==13)
						{
						HepMC::GenVertex* prodVtx2 = ith_gMu2->production_vertex();
						//Loop over the incoming particles to the second vertex
						for(HepMC::GenVertex::particles_in_const_iterator inPart2 = prodVtx2->particles_in_const_begin();inPart2 != prodVtx2->particles_in_const_end(); inPart2++)
							{
							HepMC::GenParticle* Parent2=*inPart2;
							int Parent_id2 =Parent2->pdg_id();
							//Check that a Z0 went into the vertex
							if(Parent_id2==23)
								{
								//Make sure that the charges of the muons are opposite
								if(id1*id2<0 && flag==0)
									{
									flag=1;
									double gMuonPx2  =ith_gMu2->momentum().px();
									double gMuonPy2  =ith_gMu2->momentum().py();
									double gMuonPz2  =ith_gMu2->momentum().pz();
									double gMuonE2   =ith_gMu2->momentum().e();
									TLorentzVector gMuonP2(gMuonPx2,gMuonPy2,gMuonPz2,gMuonE2);
						

									TLorentzVector gZ0_P = gMuonP1 +gMuonP2;
									double gZ0_Rap= (0.5)*log((gZ0_P.E()+gZ0_P.Pz())/(gZ0_P.E()-gZ0_P.Pz()));
									gZ0_InvMass_h->Fill(gZ0_P.Mag());
									gZ0_Rap_h->Fill(gZ0_Rap);
									gZ0_Eta_h->Fill(gZ0_P.Eta());
									gZ0_Pt_h->Fill(gZ0_P.Pt());
									gZ0_RapVSPt_h->Fill(gZ0_Rap,gZ0_P.Pt());
									gZ0_EtaVSPhi_h->Fill(gZ0_P.Eta(),gZ0_P.Phi());

									}
								}
							}
						}
					}
				}
			}
		}
	}//Loop over particles in the event
cout<<"number of muons is : "<<counter<<endl;









//     RECO
//     ->RecoMuon
edm::Handle<std::vector<reco::Muon> >rMuons;
iEvent.getByLabel("paramMuons","ParamGlobalMuons",rMuons);
//iEvent.getByLabel("paramMuons","ParamL3Muons",rMuons);
//These two product labels can be found in data generated using famosWithMuons in the path of the .cfg file

for (std::vector<reco::Muon>::const_iterator rMu1 = rMuons->begin(); rMu1 !=rMuons->end();++rMu1)
	{
	
	//reco::RecoCandidate RCand = *rMu1;
	reco::LeafCandidate ith_candMu1 = *rMu1;
	reco::Particle ith_rMu1 = *rMu1;
	reco::Muon ith_recoMu = *rMu1;
	//reco::Track ith_rMu1_track =*rMu1;
        //reco::GenParticle ith_rMu1_GenPart =*rMu1;

	//reco::MuonEnergy recoE = ith_recoMu.getCalEnergy();
	//CalEnergy_h->Fill(recoE->em());
	NChambers_h->Fill(rMu1->numberOfChambers());
	//NSegments_h->Fill(rMu1->numberOfSegments());
	rMuonEta_h->Fill(ith_rMu1.eta());
	rMuonPt_h->Fill(ith_rMu1.pt());
	TLorentzVector rMuon_P1(ith_rMu1.px(),ith_rMu1.py(),ith_rMu1.pz(),ith_rMu1.energy());
	for (std::vector<reco::Muon>::const_iterator rMu2 = rMu1+1; rMu2 != rMuons->end(); ++rMu2)
		{
		reco::Particle ith_rMu2 = *rMu2;
		TLorentzVector rMuon_P2(ith_rMu2.px(),ith_rMu2.py(),ith_rMu2.pz(),ith_rMu2.energy());

		TLorentzVector rZ0_P = rMuon_P1 + rMuon_P2;
		rZ0_InvMass_h->Fill(rZ0_P.Mag());
		rZ0_Eta_h->Fill(rZ0_P.Eta());
		rZ0_Pt_h->Fill(rZ0_P.Pt());
		rZ0_Rap_h->Fill(rZ0_P.Rapidity());
		rZ0_RapVSPt_h->Fill(rZ0_P.Rapidity(),rZ0_P.Pt());
		rZ0_EtaVSPhi_h->Fill(rZ0_P.Eta(),rZ0_P.Phi());		

		}

	}


}


// ------------ method called once each job just before starting event loop  ------------
void 
Muons::beginJob(const edm::EventSetup&)
{

gMuonEta_h          = new TH1F("gMuonEta_h","MuonMC #eta ",14,-7,7);
gMuonPt_h           = new TH1F("gMuonPt_h","MuonMC Pt ",20,0,100);
gZ0_InvMass_h       = new TH1F("gZ0_InvMass_h","MC Z^{0} ->#mu + #mu",80,60,140);
gZ0_Eta_h           = new TH1F("gZ0_Eta_h","MC Z^{0} #eta",14,-7,7);
gZ0_Rap_h           = new TH1F("gZ0_Rap_h","MC Z^{0} Rapidity",10,-5,5);
gZ0_Pt_h            = new TH1F("gZ0_Pt_h","MC Z^{0} Pt",20,0,100);
gZ0_RapVSPt_h       = new TH2F("gZ0_RapVSPt_h","MC Z^{0}   Rap VS Pt",50,-5,5,100,0,100);
gZ0_EtaVSPhi_h      = new TH2F("gZ0_EtaVSPhi_h","MC Z^{0}   #eta VS #phi",50,-5,5,50,-5,5);

rMuonEta_h          = new TH1F("rMuonEta_h","MuonReco #eta ",14,-7,7);
rMuonPt_h           = new TH1F("rMuonPt_h","MuonReco Pt ",20,0,100);
rZ0_InvMass_h       = new TH1F("rZ0_InvMass_h","Reco Z^{0} ->#mu + #mu",80,60,140);
rZ0_Eta_h           = new TH1F("rZ0_Eta_h","Reco Z^{0} #eta",14,-7,7);
rZ0_Rap_h           = new TH1F("rZ0_Rap_h","Reco Z^{0} Rapidity",10,-5,5);
rZ0_Pt_h            = new TH1F("rZ0_Pt_h","Reco Z^{0} Pt",20,0,100);
rZ0_RapVSPt_h       = new TH2F("rZ0_RapVSPt_h","Reco Z^{0}   Rap VS Pt",50,-5,5,100,0,100);
rZ0_EtaVSPhi_h      = new TH2F("rZ0_EtaVSPhi_h","Reco Z^{0}   #eta VS #phi",50,-5,5,50,-5,5);

EffMuonEta_h        = new TH1F("EffMuonEta_h","MuonEfficiency #eta ",14,-7,7);
EffMuonPt_h         = new TH1F("EffMuonPt_h","MuonEfficiency Pt ",20,0,100);
EffZ0_InvMass_h     = new TH1F("EffZ0_InvMass_h","Efficiency Z^{0} ->#mu + #mu",80,60,140);
EffZ0_Eta_h         = new TH1F("EffZ0_Eta_h","Efficiency Z^{0} #eta",14,-7,7);
EffZ0_Rap_h         = new TH1F("EffZ0_Rap_h","Efficiency Z^{0} Rapidity",10,-5,5);
EffZ0_Pt_h          = new TH1F("EffZ0_Pt_h","Efficiency Z^{0} Pt",20,0,100);

NChambers_h         = new TH1F("NChambers_h","Number of chambers of reco muon",10,0,10);
CalEnergy_h         = new TH1F("CalEnergy_h","Energy deposited in Calorimeters",1000,0,1000);
//NSegments_h         = new TH1F("NSegments","Numbers of segments reco muon",100,0,100);
}


// ------------ method called once each job just after ending the event loop  ------------
void 
Muons::endJob() {
cout<<"Writing file....."<<endl;
//Create the output file
TFile * outputFile = new TFile("MuonReco_10_175_OUT.root","RECREATE");

EffMuonEta_h->Divide(rMuonEta_h,gMuonEta_h,1,1);
EffMuonPt_h->Divide(rMuonPt_h,gMuonPt_h,1,1);
EffZ0_InvMass_h->Divide(rZ0_InvMass_h,gZ0_InvMass_h,1,1);
EffZ0_Eta_h->Divide(rZ0_Eta_h,gZ0_Eta_h,1,1);
EffZ0_Rap_h->Divide(rZ0_Rap_h,gZ0_Rap_h,1,1);
EffZ0_Pt_h->Divide(rZ0_Pt_h,gZ0_Pt_h,1,1);


gMuonEta_h->Write();
gMuonPt_h->Write();
gZ0_InvMass_h->Write();
gZ0_Eta_h->Write();
gZ0_Rap_h->Write();
gZ0_Pt_h->Write();
gZ0_RapVSPt_h->Write();
gZ0_EtaVSPhi_h->Write();

rMuonEta_h->Write();
rMuonPt_h->Write();
rZ0_InvMass_h->Write();
rZ0_Eta_h->Write();
rZ0_Rap_h->Write();
rZ0_Pt_h->Write();
rZ0_RapVSPt_h->Write();
rZ0_EtaVSPhi_h->Write();

EffMuonEta_h->Write();
EffMuonPt_h->Write();
EffZ0_InvMass_h->Write();
EffZ0_Eta_h->Write();
EffZ0_Rap_h->Write();
EffZ0_Pt_h->Write();

NChambers_h->Write();
CalEnergy_h->Write();
//NSegments_h->Write();


outputFile->Close();
}

//define this as a plug-in
DEFINE_FWK_MODULE(Muons);