///////////////////////////////////////////////////////////////////////////////
//
//
// File:     Mctest.cpp
//
// Purpose:  To test MC executables and fill some physics histograms
//
// Created:  15th Oct 2002 Jyothsna Rani
//
// History:  
//
//
///////////////////////////////////////////////////////////////////////////////

// Dependencies (#includes)
#include "mctest/Mctest.hpp"
#include "mctest/conejet.hpp"

//  includes for Event
#include "edm/IDSelector.hpp"
#include "edm/THandle.hpp"

//  includes for MC data
#include "mcpp/MCevt.hpp"   // include all headers for mcpp event containers 
#include "mcpp/MCparticleList.hpp" //  a useful class to manipulate list
#include "mcpp/MCparticleFilter.hpp"
#include "mcpp/ChargedParticleFilter.hpp" // filter to select charged particles
#include "mcpp/StableParticleFilter.hpp" // filter to select stable particles

#include "d0_mcpp/MCKineChunk.hpp"  // mcpp data
#include "kinem_util/AnglesUtil.hpp"

#include "identifiers/CollisionID.hpp"
#include "pileup_evt/MinBiasChunk.hpp"
#include "lumi_chunk/LumiInfoChunk.hpp"

#include <iostream>
#include <iomanip>

// Framework registry macro
#include "framework/Registry.hpp"
FWK_REGISTRY_IMPL(Mctest,"$Name:  $")

using namespace edm;
using namespace std;
using fwk::Context;
using fwk::Package;
using fwk::JobSummary;
using fwk::Result;

///////////////////////////////////////////////////////////////////////////////

Mctest::Mctest(Context* context): 
  Package(context), Process(context), Analyze(context), JobSummary(context),
  _processed_events(0){
  // Default constructor
  
  cout << endl
       << "+------------------------------------------+" << endl
       << "|      MCTEST - package initialization     |" << endl
       << "+------------------------------------------+" << endl
       << endl;
  
  // Access the RCP information
  RCP rcp = packageRCP();

  // Create histograms
  _do_histograms = rcp.getBool("FillHistograms");
  
  // Selection cuts for histograms -- read from rcp file
  _ElectronPTmin = rcp.getFloat("ElectronPTmin");
  _MuonPTmin = rcp.getFloat("MuonPTmin");
  _LightPartonPTmin = rcp.getFloat("LightPartonPTmin");

  // Initialize ROOT
  if ( _do_histograms ) {
    
    _histo_saved = false;
    
    // Open the file which will contain the histograms
    string histoFile = rcp.getString("HistogramFile");
    cout << " Histograms will be written to : " << histoFile << endl;
    _hepMgr = new HepRootFileManager(histoFile.c_str(), HepFileManager::HEP_RENAME);

    // Book the user histograms
    // No. of MC particles
    _h1nmc   = new TH1F("nmc","Number of MC particles",100,0,1000);
    // No. of Stable particles 
    _h1nstab = new TH1F("nstab","Number of stable particles",100,0,1000);

    // pt, eta and phi of Stable Particles
    _h1mcpt  = new TH1F("mcpt","P_{T} of stable particles",150,0.,150.);
    _h1mceta = new TH1F("mceta","#eta of stable particles",60,-6.,6.);
    _h1mcphi = new TH1F("mcphi","#phi of stable particles",65,-0.1,6.4);

    // pt, eta and phi of electrons, and parent ID
    _h1nel_all = new TH1F("nel_all","Number of electrons (no cut)",20,0,20);
    _h1elpt_all = new TH1F("elpt_all","P_{T} of electrons (no cut)",150,0.,150.);
    char title[200];
    sprintf(title, "Number of electrons (P_{T}>%d GeV)", (int)_ElectronPTmin);
    _h1nel   = new TH1F("nel",title,20,0,20);
    sprintf(title, "P_{T} of electrons (P_{T}>%d GeV)", (int)_ElectronPTmin);
    _h1elpt  = new TH1F("elpt",title,150,0.,150.);
    _h1eleta = new TH1F("eleta","#eta of electrons",60,-6.,6.);
    _h1elphi = new TH1F("elphi","#phi of electrons",65,-0.1,6.4);
    _h1elpid = new TH1F("elpid","abs(Parent ID) of electrons",100,0,100);
    _h1elpid2 = new TH1F("elpid2","abs(Parent ID) of electrons (overflows)",500,100,600);
    // muons
    _h1nmu_all = new TH1F("nmu_all","Number of muons (no cut)",20,0,20);
    _h1mupt_all = new TH1F("mupt_all","P_{T} of muons (no cut)",150,0.,150.);
    sprintf(title, "Number of muons (P_{T}>%d GeV)", (int)_MuonPTmin);
    _h1nmu   = new TH1F("nmu",title,20,0,20);
    sprintf(title, "P_{T} of muons (P_{T}>%d GeV)", (int)_MuonPTmin);
    _h1mupt  = new TH1F("mupt",title,150,0.,150.);
    _h1mueta = new TH1F("mueta","#eta of muons",60,-6.,6.);
    _h1muphi = new TH1F("muphi","#phi of muons",65,-0.1,6.4);
    _h1mupid = new TH1F("mupid","abs(Parent ID) of muons",100,0,100);
    _h1mupid2 = new TH1F("mupid2","abs(Parent ID) of muons (overflows)",500,100,600);

    // event information
    _h1xsec  = new TH1F("xsec","Cross section vs. run number",100,0,100);
    _h1pdfx1 = new TH1F("pdfx1","Parton x1",100,0,1);
    _h1pdfx2 = new TH1F("pdfx2","Parton x2",100,0,1);
    _h1nvtx  = new TH1F("nvtx","Number of vertices",100,0,1000);
    _h1pvx   = new TH1F("pvx","Primary vertex x",100,-0.1,0.1);
    _h1pvy   = new TH1F("pvy","Primary vertex y",100,-0.1,0.1);
    _h1pvz   = new TH1F("pvz","Primary vertex z",100,-100,100);
    _h1lumi  = new TH1F("lumi","Inst. lumi. of overlaid events",100,0,10);

    // generated W and Z
    _h1_w_n   = new TH1F("w_n","Number of generated W bosons",5,0,5);
    _h1_w_pt  = new TH1F("w_pt","P_{T} of generated W bosons",100,0.,50.);
    _h1_w_pt2 = new TH1F("w_pt2","P_{T} of generated W bosons (log)",100,0.,300.);
    _h1_w_eta = new TH1F("w_eta","Rapidity of generated W bosons",60,-6.,6.);
    _h1_w_phi = new TH1F("w_phi","#phi of generated W bosons",65,-0.1,6.4);
    _h1_w_mass = new TH1F("w_mass","Mass of generated W bosons",100,40,140);
    _h1_z_n   = new TH1F("z_n","Number of generated Z bosons",5,0,5);
    _h1_z_pt  = new TH1F("z_pt","P_{T} of generated Z bosons",100,0.,50.);
    _h1_z_pt2 = new TH1F("z_pt2","P_{T} of generated Z bosons (log)",100,0.,300.);
    _h1_z_eta = new TH1F("z_eta","Rapidity of generated Z bosons",60,-6.,6.);
    _h1_z_phi = new TH1F("z_phi","#phi of generated Z bosons",65,-0.1,6.4);
    _h1_z_mass = new TH1F("z_mass","Mass of generated Z bosons",100,0,600);

    // generated top quarks
    _h1_t_n   = new TH1F("t_n","Number of generated t quarks",5,0,5);
    _h1_t_pt  = new TH1F("t_pt","P_{T} of generated t quarks",100,0.,400.);
    _h1_t_eta = new TH1F("t_eta","Rapidity of generated t quarks",60,-6.,6.);
    _h1_t_phi = new TH1F("t_phi","#phi of generated t quarks",65,-0.1,6.4);
    _h1_t_mass = new TH1F("t_mass","Mass of generated t quarks",50,150,200);

    // generated b, c quarks
    _h1_b_n   = new TH1F("b_n","Number of generated b quarks",10,0,10);
    _h1_b_pt  = new TH1F("b_pt","P_{T} of generated b quarks",100,0.,100.);
    _h1_b_eta = new TH1F("b_eta","#eta of generated b quarks",60,-6.,6.);
    _h1_b_phi = new TH1F("b_phi","#phi of generated b quarks",65,-0.1,6.4);
    _h1_b1_pt = new TH1F("b1_pt","P_{T} of generated leading b",100,0.,100.);
    _h1_b2_pt = new TH1F("b2_pt","P_{T} of generated second b",100,0.,100.);
    _h1_b3_pt = new TH1F("b3_pt","P_{T} of generated third b",100,0.,100.);
    _h1_b4_pt = new TH1F("b4_pt","P_{T} of generated fourth b",100,0.,100.);
    _h1_b_drmin = new TH1F("b_drmin","#DeltaR_{min} of b quarks",50,0.,5.);
    _h1_c_n   = new TH1F("c_n","Number of generated c quarks",10,0,10);
    _h1_c_pt  = new TH1F("c_pt","P_{T} of generated c quarks",100,0.,100.);
    _h1_c_eta = new TH1F("c_eta","#eta of generated c quarks",60,-6.,6.);
    _h1_c_phi = new TH1F("c_phi","#phi of generated c quarks",65,-0.1,6.4);
    _h1_c1_pt = new TH1F("c1_pt","P_{T} of generated leading c",100,0.,100.);
    _h1_c2_pt = new TH1F("c2_pt","P_{T} of generated second c",100,0.,100.);
    _h1_c3_pt = new TH1F("c3_pt","P_{T} of generated third c",100,0.,100.);
    _h1_c4_pt = new TH1F("c4_pt","P_{T} of generated fourth c",100,0.,100.);
    _h1_c_drmin = new TH1F("c_drmin","#DeltaR_{min} of c quarks",50,0.,5.);

    // generated light-parton quarks
    _h1_lp_n_all = new TH1F("lp_n_all","Number of generated lp (no cut)",20,0,20);
    _h1_lp_pt_all = new TH1F("lp_pt_all","P_{T} of generated lp (no cut)",100,0.,100.);
    sprintf(title, "Number of generated lp (P_{T}>%d GeV)", (int)_LightPartonPTmin);
    _h1_lp_n   = new TH1F("lp_n",title,10,0,10);
    sprintf(title, "P_{T} of generated lp (P_{T}>%d GeV)", (int)_LightPartonPTmin);
    _h1_lp_pt  = new TH1F("lp_pt",title,100,0.,100.);
    _h1_lp_eta = new TH1F("lp_eta","#eta of generated light partons",60,-6.,6.);
    _h1_lp_phi = new TH1F("lp_phi","#phi of generated light partons",65,-0.1,6.4);
    _h1_lp1_pt = new TH1F("lp1_pt","P_{T} of generated leading lp",100,0.,100.);
    _h1_lp2_pt = new TH1F("lp2_pt","P_{T} of generated second lp",100,0.,100.);
    _h1_lp3_pt = new TH1F("lp3_pt","P_{T} of generated third lp",100,0.,100.);
    _h1_lp4_pt = new TH1F("lp4_pt","P_{T} of generated fourth lp",100,0.,100.);
    _h1_lp5_pt = new TH1F("lp5_pt","P_{T} of generated fifth lp",100,0.,100.);
    _h1_lp_drmin = new TH1F("lp_drmin","#DeltaR_{min} of light partons",50,0.,5.);

    // parton jets by cone algorithm
    _h1_jet_n   = new TH1F("jet_n","Number of parton jets",50,0,50);
    _h1_jet_pt  = new TH1F("jet_pt","P_{T} of parton jets",100,0.,200.);
    _h1_jet_pt2 = new TH1F("jet_pt2","P_{T} of parton jets (log)",100,0.,1000.);
    _h1_jet_n_no_e = new TH1F("jet_n_no_e","Number of parton jets (no electron)",50,0,50);
    _h1_jet_pt_no_e = new TH1F("jet_pt_no_e","P_{T} of parton jets (no electron)",100,0.,200.);
    _h1_jet_pt2_no_e = new TH1F("jet_pt2_no_e","P_{T} of parton jets (no electron, log)",100,0.,1000.);
    _h1_jet_eta = new TH1F("jet_eta","#eta of parton jets",60,-6.,6.);
    _h1_jet_phi = new TH1F("jet_phi","#phi of parton jets",65,-0.1,6.4);
    _h1_jet_mjj = new TH1F("jet_mjj","Mass of two leading jets",100,0.,300.);
    _h1_jet_dphi = new TH1F("jet_dphi","#Delta#phi of two leading jets",40,0.,4.);


  }// Initialize ROOT 

  // to average cross section
  _last_evt = 0;
  _last_xsec = 0;
  _runnum = 0;


  // Finish initialization
  cout << endl
       << "+------------------------------------------+" << endl
       << "|         Initialization finished          |" << endl
       << "+------------------------------------------+" << endl
       << endl;
  
}// Default constructor

//------------------------------------------------------------//
Result Mctest::processEvent(Event &event) {

  // Process one event
  // Increase the counter for the number of events processed
  ++_processed_events;
  return Result::success;

}// processEvent

//------------------------------------------------------------//
Result Mctest::analyzeEvent(const Event &event) {
  
  // ----------------------------------------------
  //            Access the MCKine chunk
  // ----------------------------------------------
  //  extract pointer to MC information from MCKineChunk
  std::list<edm::THandle<MCKineChunk> > mc_chunks=_mcKey.findAll(event);
  std::list<edm::THandle<MCKineChunk> >::const_iterator imc_chunk;

  int numMC   = 0;
  int numStab = 0;

  bool evt_done = false;
  // loop over MCKineChunks, grab info from 0th chunk (as long as it's valid) 
  // and end
  for (imc_chunk=mc_chunks.begin(); imc_chunk != mc_chunks.end(); imc_chunk++) {
    edm::THandle<MCKineChunk> mcChunk=*imc_chunk;     
    if ( mcChunk.isValid() ) {
      MCevt* ptevt = mcChunk->ptrMCevt();  
      if (!evt_done) {
	list<MCparticle*> listp;
	ptevt->getParticles(listp);   //  get complete list of particles

	// Number of particles
	numMC = listp.size();


	////////////////////////////////////////////////////////////
	// histograms for generated: W, Z, t, b, c, lp

	int num_w = 0;
	int num_z = 0;
	int num_t = 0;
	int num_b = 0;
	int num_c = 0;
	int num_lp = 0;
	int num_lp_all = 0;
	vector<pair<float, float> > ang_b;
	vector<pair<float, float> > ang_c;
	vector<pair<float, float> > ang_lp;

	for (list<MCparticle*>::iterator it=listp.begin(); it!=listp.end(); ++it) {
	  int pdgid = (*it)->pdgId();
	  L4vec vec4 = (*it)->vec4();
	  float pt  = vec4.pT();
	  float eta = vec4.eta();
	  float phi = vec4.phi();
	  float mass = vec4.mass();
	  float rapidity = (*it)->y();
	  MCvtx *vtx = (*it)->vtx();
	  MCvtx *vtxend = (*it)->vtxend();
	  MCparticle *parentp = vtx->aparent();
	  MCparticle *daughtp = (vtxend!=0) ? vtxend->adaughter() : 0;
	  int parentid = (parentp!=0) ? parentp->pdgId() : 0;
	  int daughtid = (daughtp!=0) ? daughtp->pdgId() : 0;

	  if (abs(pdgid)==24 && parentid==0 && daughtid!=0) { // W
	    ++num_w;
	    if (_do_histograms) {
	      _h1_w_pt->Fill(pt);
	      _h1_w_pt2->Fill(pt);
	      _h1_w_eta->Fill(rapidity);
	      _h1_w_phi->Fill(phi);
	      _h1_w_mass->Fill(mass);
	    }
	  }
	  if (abs(pdgid)==23 && parentid==0 && daughtid!=0) { // Z
	    ++num_z;
	    if (_do_histograms) {
	      _h1_z_pt->Fill(pt);
	      _h1_z_pt2->Fill(pt);
	      _h1_z_eta->Fill(rapidity);
	      _h1_z_phi->Fill(phi);
	      _h1_z_mass->Fill(mass);
	    }
	  }
	  if (abs(pdgid)==6 && parentid==0 && daughtid!=0) { // t
	    ++num_t;
	    if (_do_histograms) {
	      _h1_t_pt->Fill(pt);
	      _h1_t_eta->Fill(rapidity);
	      _h1_t_phi->Fill(phi);
	      _h1_t_mass->Fill(mass);
	    }
	  }
	  if (abs(pdgid)==5 && parentid==0 && daughtid!=0) { // b
	    ++num_b;
	    ang_b.push_back(make_pair(eta, phi));
	    if (_do_histograms) {
	      _h1_b_pt->Fill(pt);
	      _h1_b_eta->Fill(eta);
	      _h1_b_phi->Fill(phi);
	      if (num_b==1) _h1_b1_pt->Fill(pt);
	      if (num_b==2) _h1_b2_pt->Fill(pt);
	      if (num_b==3) _h1_b3_pt->Fill(pt);
	      if (num_b==4) _h1_b4_pt->Fill(pt);
	    }
	  }
	  if (abs(pdgid)==4 && parentid==0 && daughtid!=0 && mass>1.0) { // c
	    ++num_c;
	    ang_c.push_back(make_pair(eta, phi));
	    if (_do_histograms) {
	      _h1_c_pt->Fill(pt);
	      _h1_c_eta->Fill(eta);
	      _h1_c_phi->Fill(phi);
	      if (num_c==1) _h1_c1_pt->Fill(pt);
	      if (num_c==2) _h1_c2_pt->Fill(pt);
	      if (num_c==3) _h1_c3_pt->Fill(pt);
	      if (num_c==4) _h1_c4_pt->Fill(pt);
	    }
	  }
	  if ((abs(pdgid)>=1&&abs(pdgid)<=3||abs(pdgid)==9||abs(pdgid)==21
	       || abs(pdgid)==4 && mass<1.0)
	      && parentid==0 && daughtid!=0) { // lp
	    ++num_lp_all;
	    if (_do_histograms) _h1_lp_pt_all->Fill(pt);
	    if (pt>_LightPartonPTmin) {
	      ++num_lp;
	      ang_lp.push_back(make_pair(eta, phi));
	      if (_do_histograms) {
		_h1_lp_pt->Fill(pt);
		_h1_lp_eta->Fill(eta);
		_h1_lp_phi->Fill(phi);
		if (num_lp==1) _h1_lp1_pt->Fill(pt);
		if (num_lp==2) _h1_lp2_pt->Fill(pt);
		if (num_lp==3) _h1_lp3_pt->Fill(pt);
		if (num_lp==4) _h1_lp4_pt->Fill(pt);
		if (num_lp==5) _h1_lp5_pt->Fill(pt);
	      }
	    }
	  }

	}

	if (_do_histograms) {
	  _h1_w_n->Fill(num_w);
	  _h1_z_n->Fill(num_z);
	  _h1_t_n->Fill(num_t);
	  _h1_b_n->Fill(num_b);
	  _h1_c_n->Fill(num_c);
	  _h1_lp_n->Fill(num_lp);
	  _h1_lp_n_all->Fill(num_lp_all);
	}

	// Delta_R_min between (j, j), (Q, Q)
	float drmin_bb = -1; int nb = ang_b.size();
	float drmin_cc = -1; int nc = ang_c.size();
	float drmin_jj = -1; int nlp = ang_lp.size();
	for (int i=0; i<nb-1; ++i) {
	  for (int j=i+1; j<nb; ++j) {
	    float dr = kinem::delta_R(ang_b[i].first, ang_b[i].second, ang_b[j].first, ang_b[j].second);
	    if (dr<drmin_bb || drmin_bb<0) drmin_bb = dr;
	  }
	}
	for (int i=0; i<nc-1; ++i) {
	  for (int j=i+1; j<nc; ++j) {
	    float dr = kinem::delta_R(ang_c[i].first, ang_c[i].second, ang_c[j].first, ang_c[j].second);
	    if (dr<drmin_cc || drmin_cc<0) drmin_cc = dr;
	  }
	}
	for (int i=0; i<nlp-1; ++i) {
	  for (int j=i+1; j<nlp; ++j) {
	    float dr = kinem::delta_R(ang_lp[i].first, ang_lp[i].second, ang_lp[j].first, ang_lp[j].second);
	    if (dr<drmin_jj || drmin_jj<0) drmin_jj = dr;
	  }
	}
	if (_do_histograms) {
	  _h1_b_drmin->Fill(drmin_bb);
	  _h1_c_drmin->Fill(drmin_cc);
	  _h1_lp_drmin->Fill(drmin_jj);
	}


	// Find the stable Particles
	StableParticleFilter stabfilter; // selects stable particles
	// define lists
	list<MCparticle*> stab;  // list will contain stable particles

	// filter list to pick stable particles, store results in an MCparticleList
	MCparticleList stablist(&stabfilter, listp);  
	stablist.getParticles(stab);

	numStab = stab.size();

	list<MCparticle*>::iterator jj;

	Int_t nelectron = 0;
	Int_t nmuon = 0;
	Int_t nelectron_all = 0;
	Int_t nmuon_all = 0;

	// Loop over stable particles
	for ( jj=stab.begin(); jj!=stab.end(); ++jj ) {

	  L4vec vec4 = (*jj)->vec4();
	  float Pt  = vec4.pT();
	  float Eta = vec4.eta();
	  float Phi = vec4.phi();

	  Int_t pdgid = (*jj)->pdgId();

	  // find parent particle ID
	  MCvtx *vtx = (*jj)->vtx();
	  MCparticle *parentp = vtx->aparent();
	  Int_t parentid = 0;
	  if (parentp!=0) parentid = parentp->pdgId();

	  if ( _do_histograms ) {
	    _h1mcpt->Fill(Pt);
	    _h1mceta->Fill(Eta);
	    _h1mcphi->Fill(Phi);
	  }

	  if (abs(pdgid)==11) {
	    ++nelectron_all;
	    if (_do_histograms) _h1elpt_all->Fill(Pt);
	    if (Pt>_ElectronPTmin) {
	      ++nelectron;
              if (_do_histograms) {
		_h1elpt->Fill(Pt);
		_h1eleta->Fill(Eta);
		_h1elphi->Fill(Phi);
		_h1elpid->Fill(abs(parentid));
		_h1elpid2->Fill(abs(parentid));
	      }
	    }
	  }

	  if (abs(pdgid)==13) {
	    ++nmuon_all;
	    if (_do_histograms) _h1mupt_all->Fill(Pt);
	    if (Pt>_MuonPTmin) {
	      ++nmuon;
	      if (_do_histograms) {
		_h1mupt->Fill(Pt);
		_h1mueta->Fill(Eta);
		_h1muphi->Fill(Phi);
		_h1mupid->Fill(abs(parentid));
		_h1mupid2->Fill(abs(parentid));
	      }
	    }
	  }

	}// Loop

	if ( _do_histograms ) {
	  _h1nel_all->Fill(nelectron_all);
	  _h1nmu_all->Fill(nmuon_all);
	  _h1nel->Fill(nelectron);
	  _h1nmu->Fill(nmuon);
	}

	////////////////////////////////////////////////////////////
	// Extract more information from MCKineChunk

	// some event information
	Int_t run =0, evtid =0, nreac =0, flav1=0, flav2=0; 
	Float_t xsect=0, weight=0, qsq=0, shat=0, that=0, uhat=0, x1=0., x2=0.;
	Float_t xyzpv[3];

	MCevtInfo mcei=ptevt->evtInfo();
	weight  = mcei.weight();
	xsect   = mcei.crossSection();
	mcei.interaction(qsq, shat, that, uhat);
	mcei.ids(run, evtid, nreac);
	mcei.parton_x(x1,x2);
	mcei.parton_flav(flav1,flav2);
	ptevt->getPrimaryVtx(xyzpv);

	// vertex
	std::list<const MCvtx*> listv;
	ptevt->getVertices(listv);
	Int_t nvtx = listv.size();

	if ( _do_histograms ) {
	  _h1pdfx1->Fill(x1);
	  _h1pdfx2->Fill(x2);
	  _h1nvtx->Fill(nvtx);
	  _h1pvx->Fill(xyzpv[0]);
	  _h1pvy->Fill(xyzpv[1]);
	  _h1pvz->Fill(xyzpv[2]);
	}

	// to average cross section
	if (evtid < _last_evt) {
	  if (_do_histograms && _runnum<_h1xsec->GetNbinsX())
	    _h1xsec->Fill(_runnum, _last_xsec);
	  ++_runnum;
	}
	_last_evt = evtid;
	_last_xsec = xsect;


	////////////////////////////////////////////////////////////
	// Try parton jets
	list<vector<double> > jetparts;
	list<vector<double> > jetparts_no_e;
	for ( jj=stab.begin(); jj!=stab.end(); ++jj ) {
	  int pdgid = (*jj)->pdgId();
	  // remove neutrinos from list
	  if (abs(pdgid)==12 || abs(pdgid)==14 || abs(pdgid)==16) continue;
	  // remove muons which don't interact in calorimeter
	  if (abs(pdgid)==13) continue;
	  L4vec vec4 = (*jj)->vec4();
	  vector<double> momenta;
	  momenta.push_back(vec4.px());
	  momenta.push_back(vec4.py());
	  momenta.push_back(vec4.pz());
	  momenta.push_back(vec4.E());
	  jetparts.push_back(momenta);
	  if (abs(pdgid)!=11) jetparts_no_e.push_back(momenta);
	}
	ConeJet partonjet(jetparts, 0.5, 10, 0.5);
	list<vector<double> > jets;
	partonjet.GetConeJets(jets);

	ConeJet partonjet_no_e(jetparts_no_e, 0.5, 10, 0.5);
	list<vector<double> > jets_no_e;
	partonjet_no_e.GetConeJets(jets_no_e);

	if ( _do_histograms ) {
	  _h1_jet_n->Fill(jets.size());
	  for (list<vector<double> >::iterator it=jets.begin(); it!=jets.end(); ++it) {
	    _h1_jet_pt->Fill((*it)[6]);
	    _h1_jet_pt2->Fill((*it)[6]);
	    _h1_jet_eta->Fill((*it)[4]);
	    _h1_jet_phi->Fill((*it)[5]);
	  }

	  _h1_jet_n_no_e->Fill(jets_no_e.size());
	  for (list<vector<double> >::iterator it=jets_no_e.begin(); it!=jets_no_e.end(); ++it) {
	    _h1_jet_pt_no_e->Fill((*it)[6]);
	    _h1_jet_pt2_no_e->Fill((*it)[6]);
	  }

	  // Mjj and Delat_PHI_jj
	  if (jets.size()>1) {
	    list<vector<double> >::iterator it = jets.begin();
	    list<vector<double> >::iterator jt = jets.begin(); ++jt;
	    double mass = ((*it)[3]+(*jt)[3])*((*it)[3]+(*jt)[3])
	      - ((*it)[0]+(*jt)[0])*((*it)[0]+(*jt)[0])
	      - ((*it)[1]+(*jt)[1])*((*it)[1]+(*jt)[1])
	      - ((*it)[2]+(*jt)[2])*((*it)[2]+(*jt)[2]);
	    if (mass>0) mass = sqrt(mass);
	    double dphi = kinem::delta_phi((*it)[5], (*jt)[5]);
	    _h1_jet_mjj->Fill(mass);
	    _h1_jet_dphi->Fill(dphi);
	  }
	}


	evt_done=true;
      } //evt_done

    }// isValid()
  } // loop over chunks

  //cout << "MC Particles     : " << numMC << endl;
  //cout << "Stable Particles : " << numStab << endl;

  if ( _do_histograms ) {
    _h1nmc->Fill(float(numMC));
    _h1nstab->Fill(static_cast<float>(numStab));
  }

  ////////////////////////////////////////////////////////////
  // Extract data from MinBias Chunck (if available)
  Long_t overlayrun=0, overlayevtid=0;
  Int_t  overlaylumblk=0, overlaytick=0;
  Float_t inst_lum = -1.0;

  edm::TKey<MinBiasChunk> mb_key;
  edm::THandle<MinBiasChunk> handle = mb_key.find(event);
  if(handle.isValid()) {
    overlayrun    = handle->overlayEvent().runNumber();
    overlayevtid  = handle->overlayEvent().eventNumber();
    overlaylumblk = handle->overlayLumblk();
    overlaytick   = handle->overlayTick();

    // instantaneous luminosity
    inst_lum = -999.;
    edm::TKey<fwk::LumiInfoChunk> lumiKey;
    edm::THandle<fwk::LumiInfoChunk> lumiChunk = lumiKey.find(event);

    if(lumiChunk.isValid()) {
      inst_lum = lumiChunk->get_lumi_ticknum(overlaytick);
    } else{
      inst_lum = -1.;
    }
  }
  if ( _do_histograms ) _h1lumi->Fill(inst_lum);


  return Result::success;
}// analyzeEvent

//------------------------------------------------------------//
Result Mctest::jobSummary() {

  // Print a summary page and save the histograms.
  cout << endl
       << "+------------------------------------------+" << endl
       << "|            Summary of MCTEST             |" << endl
       << "+------------------------------------------+" << endl
       << "| Number of events processed " << setw(13) << _processed_events << " |" << endl
       << "+------------------------------------------+" << endl;

  if (_do_histograms) {
    // add the cross section of the last event
    if (_runnum<_h1xsec->GetNbinsX())
      _h1xsec->Fill(_runnum, _last_xsec);

    SaveHistograms();
  }

  cout << "|           Histograms saved               |" << endl
       << "+------------------------------------------+" << endl;

  return Result::success;

}// jobSummary

//------------------------------------------------------------//
Mctest::~Mctest() {

  // Default destructor, ensure that histograms have been saved.
  if (_do_histograms) {
    if ( !_histo_saved ) SaveHistograms();
    delete _hepMgr;
  }

}// Default destructor

//------------------------------------------------------------//
void Mctest::SaveHistograms() {

  // Save the histograms.
  if (_do_histograms) {
    _hepMgr->write();
    _histo_saved = true;
  } 

}// SaveHistograms