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

// Description: <one line class summary>

// Implementation:
//     <Notes on implementation>
//
//
// Original Author:  Jorge Robles
//         Created:  Thu Jul 10 21:11:46 CEST 2008
// $Id$
//
//
//////////////////////////////////////////////////////////
//							//
// 		ATTENTION!!!!!!!!!!!!!!			//
//							//
//							//
//	THIS ONE WORKS DONT MESS WITH IT		//
//							//	
//							//
//////////////////////////////////////////////////////////



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

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"
//Analysis Files
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include <DataFormats/TrackReco/interface/Track.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>

//Default Includes
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
//Partile.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"

//why not?
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"

#include <vector>
#include <iostream>

#include "PhysicsTools/CandUtils/interface/CandCombiner.h"
#include "PhysicsTools/CandUtils/interface/CandCombiner.h"
#include "DataFormats/Candidate/interface/CompositeRefCandidate.h"
#include "DataFormats/Candidate/interface/CompositeCandidate.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutRecord.h"
#include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
#include "L1Trigger/L1ExtraFromDigis/interface/L1ExtraParticlesProd.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/Candidate/interface/Particle.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Math/interface/Point3D.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"

#include <TH1.h>
#include <TH2F.h>
#include <TFile.h>
#include <TString.h>
#include <stdio.h>
#include <TLorentzVector.h>




//
// class declaration
//

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




class Zmumu : public edm::EDAnalyzer {
   public:
  TH1 *massID;
  TH1 *trackPt;
  TH1 *muonET;
  TH1 *invmass;
  TH1 *bkg_signal;
  TH2F *map_h;
  TH1 *bkg_pp;
  TH1 *bkg_pp_cut;
  TH1 *bkg_mm;
  TH1 *bkg;
  TH1 *clean;
  TH1 *invmass_cut;
  TH2F *accept;
  TH2F *accept2;
  TH1 *pz;
  TH1 *energy;
  TH1 *rapid;
  TH1 *MuonEta;
  TH1 *MuonInvmass;
  TH1 *MuonPt;   
  TH1 *MuonEta_PM;
  TH1 *MuonInvmass_PM;
  TH1 *MuonPt_PM;
  TH1 *MCMuonRap_PM;
TF1 *myfit;
      explicit Zmumu(const edm::ParameterSet&);
      ~Zmumu();


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

      // ----------member data ---------------------------
edm::InputTag genParticles_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//

Zmumu::Zmumu(const edm::ParameterSet& iConfig):
genParticles_(iConfig.getUntrackedParameter<edm::InputTag>("src"))
{
   //now do what ever initialization is needed

}


Zmumu::~Zmumu()
{
 
   // 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
Zmumu::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  cout<<"Im here"<< endl;
  //  using namespace edm;
  ///****************
  //HepMC
  double massMu = 0.105658369;
  double counter =0; 
  edm::Handle<edm::HepMCProduct> genEvent;
  iEvent.getByLabel("source","",genEvent);
  
  for (HepMC::GenEvent::particle_const_iterator particle1 = genEvent->GetEvent()->particles_begin(); particle1 !=genEvent->GetEvent()->particles_end();++particle1)
    {   
      counter=counter+1;
    
      HepMC::GenParticle* genPart1=*particle1;
      
      double mass = genPart1->generated_mass();
      int id1 = genPart1->pdg_id();
      massID->Fill(mass);
      
    
       if(abs(id1)==13)
	{
	  
	  //for documentation in all this schwagg see
	  //http://projects.hepforge.org/rivet/code/hepmc/namespaceHepMC.html
	  HepMC::GenVertex* prodVtx=genPart1->production_vertex();
	  for(HepMC::GenVertex::particles_in_const_iterator inPart1 = prodVtx->particles_in_const_begin();inPart1 !=prodVtx->particles_in_const_end();inPart1++ )
	    {
	      HepMC::GenParticle* Parent1=*inPart1;
	      int Parent_id1 = Parent1->pdg_id();
	      //cout<<"Who's your daddy? :  "<<Parent_id<<endl;
	       if(Parent_id1==23)
	      	{ 
		  
		
		  cout<<"counter = "<<counter<<"*****************"<<endl;
		  //These are from the class HEpMC::FourVector which does not translate into TLorentzVector
		  double Muon_E1         = (*particle1)->momentum().e();
		  double Muon_pz1        = (*particle1)->momentum().z();
		  double Muon_py1        = (*particle1)->momentum().y();
		  double Muon_px1        = (*particle1)->momentum().x();
		  double Muon_pmag1      = (*particle1)->momentum().mag();
		  TVector3 mom1(Muon_px1,Muon_py1,Muon_pz1); 
		  //cout<<"The x_momentum is :  "<<Muon_px1<< " from TVector3: " << mom1.x() <<endl;
		  //cout<<"The y_momentum is :  "<<Muon_py1<< " from TVector3: " << mom1.y() <<endl;
		  //cout<<"The z_momentum is :  "<<Muon_pz1<< " from TVector3: " << mom1.z() <<endl;
		  //cout<<"The e_momentum is :  "<<Muon_E1<<endl;
		  //cout<<"Magnitude by def:  "<<Muon_pmag1<<"   Magnitude by root is :  "<<mom1.Mag()<<endl;
		  
		  double MuonEnergy1     = sqrt((Muon_pmag1*Muon_pmag1)+(massMu*massMu));
		  TLorentzVector Muon_p1(Muon_px1,Muon_py1,Muon_pz1,MuonEnergy1); 
		  
		  //cout<<"The culculated energy : "<<MuonEnergy<<endl;
		  //cout<<"The magnitude of momentum is :  "<<Muon_mom_mag <<endl;
		  
		  muonET->Fill(Muon_E1);
		  double Muon_pt1=sqrt((Muon_px1*Muon_px1)+(Muon_py1*Muon_py1));
		  //cout<<"Muon pt :  "<<Muon_pt1<<endl;
		  //cout<<"energy :  "<<Muon_E1<<endl;
		  double muon_eta1 =(0.5) *log( (Muon_pmag1+Muon_pz1)/(Muon_pmag1-Muon_pz1)  ) ;
		  //cout<<"The eta of muon is :  "<<muon_eta<<endl;
		  MuonEta->Fill(muon_eta1);
		  MuonPt->Fill(Muon_pt1);	   
		
		  //cout<<"Point five"<<endl;//For debugging purposes only
		  int flag=0;
		  for (HepMC::GenEvent::particle_const_iterator particle2 = particle1; particle2 !=genEvent->GetEvent()->particles_end();)
		    {
                      //cout<<"Point six"<<endl;//ditto
		      particle2++;
		      if (particle2==genEvent->GetEvent()->particles_end()) continue;
		      
		      HepMC::GenParticle* genPart2=*particle2;
		      int id2=genPart2->pdg_id();
		      if(abs(id2)==13)
			{
			  //cout<<"Point seven"<<endl;//ditto
			  //Loop over all the particle that come into the vertex that produced the particle
			  //check for thei id an meke sure they are the Z0(pdg_id=23)
			  HepMC::GenVertex* prodVtx2=genPart2->production_vertex();
			  for(HepMC::GenVertex::particles_in_const_iterator inPart2 = prodVtx2->particles_in_const_begin();inPart2 !=prodVtx2->particles_in_const_end();inPart2++ )
			    {
			      //cout<<"Point eigth"<<endl;//ditto
			      HepMC::GenParticle* Parent2=*inPart2;
			      int Parent_id2 = Parent2->pdg_id();
			      	      if (Parent_id2==23)
			      	{
				  double Muon_E2    	  = (*particle2)->momentum().e();
				  double Muon_pz2   	  = (*particle2)->momentum().z();
				  double Muon_py2   	  = (*particle2)->momentum().y();
				  double Muon_px2         = (*particle2)->momentum().x();
				  double Muon_pmag2 	  = (*particle2)->momentum().mag();
				  
				  TVector3 mom2(Muon_px2,Muon_py2,Muon_pz2);
				  double MuonEnergy2     = sqrt((Muon_pmag2*Muon_pmag2)+(massMu*massMu));
				  TLorentzVector Muon_p2(Muon_px2,Muon_py2,Muon_pz2,MuonEnergy2);
				  TLorentzVector Muon_p3 = Muon_p1 + Muon_p2;
				  double Muon_inv = Muon_p3.Mag();
				  
				  //cout<<"The invariant mass is :  "<<Muon_inv<<endl;
				  MuonInvmass->Fill(Muon_inv);
				  
				  
				  
	//Get +/- pairs
	//The flag is to make sure that the histograms only get filled only once with each muon from the first loop
				  if (id1*id2<0 && flag==0)
				    {
				      flag=1;
				      double MCMuon_rap3 = (0.5)*log((Muon_p3.E()+Muon_p3.Pz())/(Muon_p3.E()-Muon_p3.Pz()));
				      MCMuonRap_PM->Fill(MCMuon_rap3);
				      MuonInvmass_PM->Fill(Muon_inv);
				      MuonEta_PM->Fill(muon_eta1);
				      MuonPt_PM->Fill(Muon_pt1);
				    }
				  }//Checking parent id 2
			    }
			}//Checking id particle 2
		    }
		  }//Cheking Parent id 1
	    }
	}//Checking id particle 1
    }//This is from the 1st foor loop
  //*************/
  //RECO
  ///*****************************
  using namespace reco;
  //Reco Tracks
  
  //NOTE:In the Handle Template used to have <reco::TrackCollection>, this worked for 
  //CMSSW_1_6_8 but in CMSSW_1_7_5 they didn't bother to write the necesary lines of code.
  //reco::Particle* candidate=*recoTracks; 
  //if (candidate->pdgId==13 | candidate->pdgId==-13){
  
  Handle<std::vector<Track> >recoTracks;
  iEvent.getByLabel("gsWithMaterialTracks","",recoTracks);
  cout<<"Reco tarck size=  "<<recoTracks->size()<<endl;
  for(std::vector<Track>::const_iterator track1 = recoTracks->begin(); track1 != recoTracks->end(); ++track1)
    {
      //reco::Particle* candidate=*recoTracks; 
      //if (candidate->pdgId==13 | candidate->pdgId==-13){
      //cout<<"Im also here"<<endl;
      double massMu = 0.105658369;
      trackPt->Fill(track1->pt());
      double eta_j=(track1->outerEta());
      double phi_j=(track1->outerPhi());
      double pt_j=(track1->outerPt());
      map_h->Fill(eta_j,phi_j);
      accept->Fill(eta_j,pt_j);
      double px_j=(track1->outerPx());
      double py_j=(track1->outerPy());
      double pz_j=(track1->outerPz());
      double ch1=(track1->charge());
      TVector3 mom1(px_j,py_j,pz_j);
      double E1=sqrt(mom1.Mag2()+(massMu*massMu));
      TLorentzVector p1(mom1,E1);
      pz->Fill(pz_j);
      energy->Fill(E1);
      //cout<<"E1:"<<E1<<endl;
      //cout<<"Pz:"<<pz_j<<endl;
      double rapidity = (1/2)*log((E1+pz_j)/(E1-pz_j));
      rapid->Fill(rapidity);
      accept2->Fill(rapidity,pt_j);
      //if (reco::Particle->pdgId==13 |reco::Particle->pdgId==-13){	  
      for(std::vector<Track>::const_iterator track2 = track1+1; track2 != recoTracks->end(); ++track2)
	{
	  double px2_j=(track2->outerPx());
	  double py2_j=(track2->outerPy());
	  double pz2_j=(track2->outerPz());
	  double ch2=(track2->charge());
	  TVector3 mom2(px2_j,py2_j,pz2_j);
	  double E2=sqrt(mom2.Mag2()+(massMu*massMu));
	  TLorentzVector p2(mom2,E2);
	  
	  TLorentzVector p3=p1+p2;
	  
	  bkg_signal->Fill(p3.Mag());
	  if (ch1*ch2<0)
	    {
	      invmass->Fill(p3.Mag());
	      if (p1.Mag()>10 | p2.Mag()>10)
		{
		  continue;
		}
	      else
		{
		  invmass_cut->Fill(p3.Mag());
		  //  printf("***********\n**********\n******\n*******\n*********\n");
		}
	    }
	  if (ch1>0 && ch2>0){
	    bkg_pp->Fill(p3.Mag());
	    if (p1.Mag()>10 | p2.Mag()>10 )
	      {
		bkg_pp_cut->Fill(p3.Mag());
	      }
	  }
	  if(ch1<0 && ch2<0){
	    bkg_mm->Fill(p3.Mag());
	  }
	}	    
    }
  //	}
  invmass->SetTitle("Invariant mass;GeV/c^{2};Counts");
  accept2->SetTitle(";Rapidity;P_{t}[GeV]");
  accept->SetTitle("Acceptance;#eta;P_{t}[GeV]");
  //*******************************/
  // invmass_cut->Add(invmass_cut,bkg_pp_cut,1,-1.41);
  //}
  //Now look at the reco:: prticle information
  using namespace std;
  //Handle<std::vector<Particle> >particles;
  //Handle<std::vector<Track> >recoTracks;
  //iEvent.getByLabel("source","",particles);
  
  //genParticles taking
  
  /***** 21
	 int numberofMuons=0;
	 
	 cout<<"Point three"<<endl;
	 
	 edm::OwnVector<reco::CandidateCollection,edm::ClonePolicy<reco::CandidateCollection> >genParticles;
	 cout<<"Point four"<<endl;
	 iEvent.getByLabel("genParticles_", genParticles);
	 cout<<"Point five"<< endl;
	 
	 
	 //cout<<"Number of  particles=  "<< genParticles->size()<<endl;
	 //iEvent.getByLabel("source,"",genParticles");
	 //for(std::vector<Track>::const_iterator track1 = recoTracks->begin(); track1 != recoTracks->end(); ++track1)
	 
	 

	 unsigned int NOParticles =genParticles->size();
	 cout<<"im here too"<<endl;
	 cout<<"the number of particles is :  "<<NOParticles<<endl;
	 //static P * clone( const P & p );
	 for(size_t k=0; k<genParticles->size(); ++k)
	 {
	 const Candidate& p = (*genParticles)[k];
	 int id = p.pdgId();
	 if(abs(id)==13)
	 {
	 cout<<"i found a fucking muon  *************"<<endl;
	 numberofMuons++;
	 }
	 }
	 *****21**/
}

/****
     Handle<reco::Candidate>candidate;
     iEvent.getByLabel("candidate_",candidate);
     
     for(reco::Candidate::const_iterator particle1 = candidate->begin();particle1 != candidate->end(); ++particle1)
     {	
     int id = reco::Particle->pdgId();
     if (abs(id)==13)
     {
     cout<<"i found a muon"<<endl;
     }
     }
****/
//cout<<"Number of muons is:"<<numberofMuons<<endl;
//This one is to loop over all the tracks-----Don't mess with it


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

massID           = new TH1F("massID","Mass generated distribution",100,0,2); 
trackPt          = new TH1F("trackPt","Transverse momentum of tracks",100,0,15);
muonET           = new TH1F("muonE","Energy of muons",100,0,100);
map_h            = new TH2F("map_h","#eta #phi distribution",200,-4,4,100,0,4);
invmass          = new TH1F("invmass","Invariant mass//",40,80,100);
bkg_signal       = new TH1F("bkg_signal","Background and signal",100,0,100);
bkg_pp           = new TH1F("bkg_pp","Background for muons",100,0,100);
bkg_pp_cut       = new TH1F("bkg_pp_cut","Background for muons with cut",100,0,100);
bkg_mm           = new TH1F("bkg_mm","Background for muons",100,0,100);
bkg              = new TH1F("bkg","Background",100,0,100);
clean            = new TH1F("clean","clean signal",100,0,100);
invmass_cut      = new TH1F("invmass_cut","Cuts>20,<75 GeV",100,0,100);
accept           = new TH2F("accept","#eta vs P_{t} ",100,-5,5,100,0,50);
accept2          = new TH2F("accept2","acceptance",100,-1,1,100,0,50);
pz               = new TH1F("pz","pz",100,0,100);
energy           = new TH1F("energy","energy",100,0,100);
rapid            = new TH1F("rapid","rapidity",200,-100,100);
MuonEta          = new TH1F("MuonEta","Muons #eta MC",100,-5,5);
MuonInvmass      = new TH1F("MuonInvmass","Muons invariant mass",80,60,140);
MuonPt           = new TH1F("MuonPt","Pt of Muons",100,0,100);
MuonEta_PM       = new TH1F("MuonEta_PM","Muons #eta +- MC",100,-5,5);
MuonInvmass_PM   = new TH1F("MuonInvmass_PM","Muons invariant mass +-",80,60,140);
MuonPt_PM        = new TH1F("MuonPt_PM","Pt of Muons +-",100,0,100);
MCMuonRap_PM     = new TH1F("MCMuonRap_PM","Z0 rapidity",120,-6,6);
myfit            = new TF1("myfit","gaus(0)",80,100);
}


// ------------ method called once each job just after ending the event loop  ------------
void 
Zmumu::endJob() {
  cout<<"Writing file"<<endl;
 //Create the output file
  TFile * outputFile = new TFile("Zmumu_500_167_out.root","RECREATE");
  massID->Write();
  trackPt->Write();
  muonET->Write();
  map_h->Write();
  invmass->Write();
  bkg_signal->Write();
  bkg->Write();
  bkg_pp->Write();
  bkg_pp_cut->Write();
  bkg_mm->Write();
  clean->Write();
  invmass_cut->Write();
  accept->Write();
  accept2->Write();
  pz->Write();
  energy->Write();
  rapid->Write();
  MuonEta->Write();
  MuonInvmass->Write();
  MuonPt->Write();
  MuonEta_PM->Write();
  MuonInvmass_PM->Write();
  MuonPt_PM->Write();
  MCMuonRap_PM->Write();
  outputFile->Close();

}

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