//
// This is an example reconstruction and simple analysis
// program for LCD events.
//
import java.util.*;

import hep.lcd.util.driver.Driver;
import hep.lcd.util.driver.AbstractProcessor;

import hep.lcd.event.LCDEvent;

import hep.lcd.mc.smear.SmearDriver;
import hep.lcd.recon.tracking.TrackReco;
import hep.lcd.event.Track;
import hep.lcd.event.TrackList;

import hep.lcd.recon.cluster.simple.SimpleClusterBuilder;
import hep.lcd.event.Cluster;
import hep.lcd.event.ClusterList;

import hep.physics.BasicHepLorentzVector;
import hep.physics.jets.AbstractJetFinder;
import hep.physics.jets.JadeEJetFinder;
import hep.physics.jets.DurhamJetFinder;

import hep.lcd.geometry.Detector;

final public class ExampleReconstruction extends Driver
{
	public ExampleReconstruction()
	{
 
		// Smear Tracker hits with default resolutions
		add(new SmearDriver());
		// Find tracks in combined inner tracker
		add(new TrackReco()); 
		// Find Clusters using simple clustering algorithm
		add(new SimpleClusterBuilder());
		// Find jets using tracks
		add(new TrackJetFinder());
		// Find jets using clusters
		add(new ClusterJetFinder());
 
	}
}

//Find Jets using reconstructed Tracks
class TrackJetFinder extends AbstractProcessor
{
	// The jet finder
	AbstractJetFinder _jetFinder;

	// Default Constructor
	public TrackJetFinder()
	{
		_jetFinder = new JadeEJetFinder(0.02);
	}

	// Fully qualified Constructor
	public TrackJetFinder( AbstractJetFinder jf)
	{
		_jetFinder = jf;
	}   

	public void process(final LCDEvent event)
	{
		// We will use the charged tracks in this example to find jets
		// Get the list of tracks from the event
		TrackList tracks = event.getTrackList();
		//feed them to the jet finder
		_jetFinder.setEvent(tracks.getTracks());
 
		// Fill some diagnostic histograms 
		histogram("Tracks Njets").fill(_jetFinder.njets());
 
		// Loop over the found jets
		double totalJetEnergy = 0;
		double totalParticleEnergy = 0;
		for (int ijet =0; ijet<_jetFinder.njets(); ijet++)
		{
			totalJetEnergy += _jetFinder.jet(ijet).t();
			histogram("Tracks Jet Energy").fill(_jetFinder.jet(ijet).t());
			histogram("Number of Tracks in Jet").fill(_jetFinder.nParticlesPerJet(ijet));
			// loop over all the tracks in the jet, and add up their energy
			// (should be the same as the jet energy)
			double etot = 0;
			Enumeration e = _jetFinder.particlesInJet(ijet);
			while (e.hasMoreElements()) 
				etot += ((Track) e.nextElement()).mag();
			histogram("Tracks Ptot").fill(etot);
			totalParticleEnergy += etot;
		}
		histogram("Total Track Jet Energy").fill(totalJetEnergy);
		histogram("Total Track Particle Momentum").fill(totalParticleEnergy);
	}
}

// Find Jets using Calorimeter Clusters
class ClusterJetFinder extends AbstractProcessor
{
	// the jet finder
	AbstractJetFinder _jetFinder;

	// Default Constructor
	public ClusterJetFinder()
	{
		_jetFinder = new DurhamJetFinder(0.02);
	}

	// Fully qualified Constructor
	public ClusterJetFinder( AbstractJetFinder jf)
	{
		_jetFinder = jf;
	}

	public void process(final LCDEvent event)
	{	
		// We will use the calorimeter clusters in this example to find jets	   
		Vector finals = new Vector();
		// Get the EM clusters from the event...
		ClusterList em = event.getEMClusters();
		Enumeration eem = em.getClusters();
		// Get the hadronic clusters from the event...
		ClusterList had = event.getHADClusters();
		Enumeration ehad = had.getClusters();

		// which detector are we using? Needed to apply sampling fractions.
		// This will be automatic in the future.
		Detector det = event.getCurrentDetector();
		// defaults...
		double emscale = 1/.0697; //cng
		double hadscale = 1/.0189; //cng		
		if(det.toString().equals("ldmar01"))
		{
			emscale = 1/.0263; 
			hadscale = 1/.028;
		}
		else if(det.toString().equals("sdmar01"))
		{
			emscale = 1/.0222;  
			hadscale = 1/.0555;  	
		}
		else if(det.toString().equals("pdmar01"))
		{
			emscale = 1/.0697;  
			hadscale = 1/.0189;  	
		}

		int emcount = 0;
		int hadcount = 0;
		int emcutcount = 0;
		int hadcutcount = 0;
		double emcutenergy = 0;
		double cEn = 0;
		double pt = 0;
		double px = 0;
		double py = 0;
		double pz = 0;
		double emEn = 0;
		double hadEn = 0;
		double eTot = 0;
		double xTot =0;
		double yTot = 0;
		double zTot = 0;

		// loop over the EM clusters
		while (eem.hasMoreElements())
		{
			emcount = emcount + 1;
			Cluster c = (Cluster) eem.nextElement();
			cEn = emscale*c.getEnergy();
			histogram("Total Cluster Energy").fill(cEn);
			histogram("EM Cluster Energy").fill(cEn);
			histogram("EM Theta").fill(c.getEnergyTheta());
			histogram("EM CosTheta").fill(Math.cos(c.getEnergyTheta()));
			pt = cEn;
			px = cEn*Math.sin(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
			py = cEn*Math.cos(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
			pz = cEn*Math.cos(c.getEnergyTheta());
			emEn = emEn + pt;
			eTot = eTot + pt;
			xTot = xTot + px;
			yTot = yTot + py;
			zTot = zTot + pz;
			// cut out low energy clusters...
			if ( cEn > 0.5 )
			{
				emcutcount++;
				emcutenergy = emcutenergy + cEn;
				finals.addElement(new BasicHepLorentzVector(pt,px,py,pz));	
			}
		}
		// Loop over the hadronic clusters
		double hadcutenergy = 0;
		while (ehad.hasMoreElements())
		{
			hadcount = hadcount + 1;
			Cluster c = (Cluster) ehad.nextElement();
			// scale the energy
			cEn = hadscale*c.getEnergy();
			histogram("HAD Cluster Energy").fill(cEn);
			histogram("Total Cluster Energy").fill(cEn);
			histogram("HAD Theta").fill(c.getEnergyTheta());
			histogram("HAD CosTheta").fill(Math.cos(c.getEnergyTheta()));
			pt = cEn;
			px = cEn*Math.sin(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
			py = cEn*Math.cos(c.getEnergyPhi())*Math.sin(c.getEnergyTheta());
			pz = cEn*Math.cos(c.getEnergyTheta());
			hadEn = hadEn + pt;
			eTot = eTot + pt;
			xTot = xTot + px;
			yTot = yTot + py;
			zTot = zTot + pz;
			if ( cEn > 0.5 )
			{
				hadcutcount++;
				hadcutenergy = hadcutenergy + cEn;
				finals.addElement(new BasicHepLorentzVector(pt + 0.137,px,py,pz));	
			}
		}  
 
		// find the jets... 
		_jetFinder.setEvent(finals.elements());
		int njets = _jetFinder.njets();
		histogram("Durham clusters Njets").fill(njets);
		double totalJetEnergy = 0;
		double totalParticleEnergy = 0;
		for (int ijet =0; ijet<_jetFinder.njets(); ijet++)
		{
			totalJetEnergy += _jetFinder.jet(ijet).t();
			histogram("Durham Cluster Jet Energy").fill(_jetFinder.jet(ijet).t());
			histogram("Durham Number of Clusters in Jet").fill(_jetFinder.nParticlesPerJet(ijet));
/*
			// loop over all the clusters in the jet, and add up their energy
			// (should be the same as the jet energy)
			double etot = 0;
			Enumeration e = _jetFinder.particlesInJet(ijet);
			while (e.hasMoreElements()) 
				etot += ((BasicHepLorentzVector) e.nextElement()).lorMag();
			histogram("Durham Ptot").fill(etot);
			totalParticleEnergy += etot;
			*/
		}
		histogram("Durham Total Jet Energy").fill(totalJetEnergy);
//		histogram("Durham Total Cluster Momentum").fill(totalParticleEnergy);		
	}
}