Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 2)
+++ trunk/Delphes.cpp	(revision 11)
@@ -19,9 +19,4 @@
 #include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
 
-#include "Utilities/CDFCones/interface/JetCluAlgorithm.h"
-#include "Utilities/CDFCones/interface/MidPointAlgorithm.h"
-#include "Utilities/CDFCones/interface/PhysicsTower.h"
-#include "Utilities/CDFCones/interface/Cluster.h"
-
 #include "H_BeamParticle.h"
 #include "H_BeamLine.h"
@@ -34,4 +29,23 @@
 
 #include "interface/SmearUtil.h"
+#include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
+#include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"
+
+// get info on how fastjet was configured
+#include "Utilities/Fastjet/include/fastjet/config.h"
+
+// make sure we have what is needed
+#ifdef ENABLE_PLUGIN_SISCONE
+#  include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
+#endif
+#ifdef ENABLE_PLUGIN_CDFCONES
+#  include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
+#  include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
+#endif
+
+#include<vector>
+#include<iostream>
+
+
 
 using namespace std;
@@ -154,6 +168,9 @@
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
-  vector<Cluster> jets; 
-
+
+  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
+  vector<fastjet::PseudoJet> inclusive_jets;
+  vector<fastjet::PseudoJet> sorted_jets;
+  
   //Initialisation of Hector
   extern bool relative_energy;
@@ -161,7 +178,7 @@
   extern int kickers_on;
   kickers_on = 1;	  // should always be 1
-
-  	// user should provide : (1) optics file for each beamline, and IPname, 
-	  // and offset data (s,x) for optical elements
+  
+  // user should provide : (1) optics file for each beamline, and IPname, 
+  // and offset data (s,x) for optical elements
   H_BeamLine* beamline1 = new H_BeamLine(1,500.);
   beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
@@ -171,5 +188,5 @@
   beamline1->add(rp220_1);
   beamline1->add(rp420_1);
-
+  
   H_BeamLine* beamline2 = new H_BeamLine(1,500.);
   beamline2->fill("data/LHCB1IR5_v6.500.tfs",-1,"IP5");
@@ -197,4 +214,5 @@
       TrackCentral.clear();
       towers.clear();
+      input_particles.clear();
       TSimpleArray<TRootGenParticle> NFCentralQ;
 
@@ -266,116 +284,119 @@
 		      break;
 		} // switch (pid)
-
+		
 		// all final particles but muons and neutrinos	
 		// for calorimetric towers and mission PT
 		if(genMomentum.E()!=0) {
-			PTmis = PTmis + genMomentum;//ptmis
-			if(pid !=pMU) {
-			      towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
-			      elementCalo = (TRootCalo*) branchCalo->NewEntry();
-			      elementCalo->Set(genMomentum);
-			}
+		  PTmis = PTmis + genMomentum;//ptmis
+		  if(pid !=pMU) {
+		    towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
+		    // create a fastjet::PseudoJet with these components and put it onto
+		    // back of the input_particles vector
+		    input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
+		    elementCalo = (TRootCalo*) branchCalo->NewEntry();
+		    elementCalo->Set(genMomentum);
+		  }
 		}
-
-
+		
+		
 		// all final charged particles
 		if(
-		     ((rand()%100) < DET->TRACKING_EFF)  &&
-		     (genMomentum.E()!=0) &&
-		     (fabs(particle->Eta) < DET->MAX_TRACKER) && 
-		     (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
-		     (pid != pGAMMA)    && 
-		     (pid != pPI0)      &&
-		     (pid != pK0L)      &&
-		     (pid != pN)        &&
-		     (pid != pSIGMA0)   &&
-		     (pid != pDELTA0)   &&
-		     (pid != pK0S)   // not charged particles : invisible by tracker
-                 )
-		    {
-		      elementTracks = (TRootTracks*) branchTracks->NewEntry();
-		      elementTracks->Set(genMomentum);
-		      TrackCentral.push_back(genMomentum);
-		    }
-	    } // switch
-
+		   ((rand()%100) < DET->TRACKING_EFF)  &&
+		   (genMomentum.E()!=0) &&
+		   (fabs(particle->Eta) < DET->MAX_TRACKER) && 
+		   (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
+		   (pid != pGAMMA)    && 
+		   (pid != pPI0)      &&
+		   (pid != pK0L)      &&
+		   (pid != pN)        &&
+		   (pid != pSIGMA0)   &&
+		   (pid != pDELTA0)   &&
+		   (pid != pK0S)   // not charged particles : invisible by tracker
+		   )
+		  {
+		    elementTracks = (TRootTracks*) branchTracks->NewEntry();
+		    elementTracks->Set(genMomentum);
+		    TrackCentral.push_back(genMomentum);
+		  }
+	  } // switch
+	  
 	  // Forward particles in CASTOR ?
-/*	  if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 
-                                    && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
-		
-
-	  } // CASTOR
-*/
+	  /*	  if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 
+		  && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
+		  
+		  
+		  } // CASTOR
+	  */
 	  // Zero degree calorimeter, for forward neutrons and photons
 	  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->MIN_ZDC ) {
-// !!!!!!!!! vÃ©rifier que particle->Z est bien en micromÃštres!!!
-// !!!!!!!!! vÃ©rifier que particle->T est bien en secondes!!!
-// !!!!!!!!! pas de smearing ! on garde trop d'info !
-		elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
-		elementZdc->Set(genMomentum);
-
-		// time of flight t is t = T + d/[ cos(theta) v ]
-		//double tx = acos(particle->Px/particle->Pz);
-		//double ty = acos(particle->Py/particle->Pz);
-		//double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
-		//double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
-		double flight_distance = (DET->ZDC_S - particle->Z*(1E-6)); 
-			// assumes also that the emission angle is so small that 1/(cos theta) = 1
-		elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles 
-		elementZdc->side = sign(particle->Eta);
-
+	    // !!!!!!!!! vÃ©rifier que particle->Z est bien en micromÃštres!!!
+	    // !!!!!!!!! vÃ©rifier que particle->T est bien en secondes!!!
+	    // !!!!!!!!! pas de smearing ! on garde trop d'info !
+	    elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
+	    elementZdc->Set(genMomentum);
+	    
+	    // time of flight t is t = T + d/[ cos(theta) v ]
+	    //double tx = acos(particle->Px/particle->Pz);
+	    //double ty = acos(particle->Py/particle->Pz);
+	    //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
+	    //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
+	    double flight_distance = (DET->ZDC_S - particle->Z*(1E-6)); 
+	    // assumes also that the emission angle is so small that 1/(cos theta) = 1
+	    elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles 
+	    elementZdc->side = sign(particle->Eta);
+	    
 	  }
-
+	  
 	  // if forward proton
           if( (pid == pP)  && (particle->Status == 1) &&  (fabs(particle->Eta) > DET->MAX_CALO_FWD) ) 
             {
-// !!!!!!!! put here particle->CHARGE and particle->MASS
-		H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
-		p1.smearAng();
-                p1.smearPos();
-                p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
-		p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
-
-		H_BeamLine *beamline;
-		if(particle->Eta >0) beamline = beamline1;
-		else beamline = beamline2;
-
-		p1.computePath(beamline,1);
-
-		if(p1.stopped(beamline)) {
- 		    if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
-			p1.propagate(DET->RP220_S);
-			elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
-	                elementRP220->X  = (1E-6)*p1.getX();  // [m]
-        	        elementRP220->Y  = (1E-6)*p1.getY();  // [m]
-                	elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
-	                elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
-        	        elementRP220->S = p1.getS();          // [m]
-                	elementRP220->T = -1; 		// not yet implemented
-                        elementRP220->E = p1.getE(); 	// not yet implemented
-                        elementRP220->q2 = -1; 		// not yet implemented
-                	elementRP220->side = sign(particle->Eta); 
-
-		    } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
-			p1.propagate(DET->FP420_S);
-                        elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
-                        elementFP420->X  = (1E-6)*p1.getX();  // [m]
-                        elementFP420->Y  = (1E-6)*p1.getY();  // [m]
-                        elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
-                        elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
-                        elementFP420->S = p1.getS();	      // [m]
-                        elementFP420->T = -1; 		      // not yet implemented
-                        elementFP420->E = p1.getE();          // not yet implemented
-                        elementFP420->q2 = -1;                // not yet implemented
-                        elementFP420->side = sign(particle->Eta);
-                    }
+	      // !!!!!!!! put here particle->CHARGE and particle->MASS
+	      H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
+	      p1.smearAng();
+	      p1.smearPos();
+	      p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
+	      p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
+	      
+	      H_BeamLine *beamline;
+	      if(particle->Eta >0) beamline = beamline1;
+	      else beamline = beamline2;
+	      
+	      p1.computePath(beamline,1);
+	      
+	      if(p1.stopped(beamline)) {
+		if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
+		  p1.propagate(DET->RP220_S);
+		  elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
+		  elementRP220->X  = (1E-6)*p1.getX();  // [m]
+		  elementRP220->Y  = (1E-6)*p1.getY();  // [m]
+		  elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
+		  elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
+		  elementRP220->S = p1.getS();          // [m]
+		  elementRP220->T = -1; 		// not yet implemented
+		  elementRP220->E = p1.getE(); 	// not yet implemented
+		  elementRP220->q2 = -1; 		// not yet implemented
+		  elementRP220->side = sign(particle->Eta); 
+		  
+		} else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
+		  p1.propagate(DET->FP420_S);
+		  elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
+		  elementFP420->X  = (1E-6)*p1.getX();  // [m]
+		  elementFP420->Y  = (1E-6)*p1.getY();  // [m]
+		  elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
+		  elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
+		  elementFP420->S = p1.getS();	      // [m]
+		  elementFP420->T = -1; 		      // not yet implemented
+		  elementFP420->E = p1.getE();          // not yet implemented
+		  elementFP420->q2 = -1;                // not yet implemented
+		  elementFP420->side = sign(particle->Eta);
 		}
-
-//		if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
-//			cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
+	      }
+	      
+	      //		if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
+	      //			cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
 	    } // if forward proton
-
+	  
 	} // while
-
+      
       // computes the Missing Transverse Momentum
       elementEtmis = (TRootETmis*) branchETmis->NewEntry();
@@ -384,56 +405,96 @@
       elementEtmis->Px = (-PTmis).Px();
       elementEtmis->Py = (-PTmis).Py();
-
+      
       //*****************************
-      jets.clear();
+
+      // we will have four jet definitions, and the first three will be
+      // plugins
+      fastjet::JetDefinition jet_def;
+      fastjet::JetDefinition::Plugin * plugins;
+      
       switch(DET->JETALGO) {
-	default:
-	case 1: {
-	    JetCluAlgorithm jetAlgoC(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
-	    jetAlgoC.run(towers, jets);
-	  } 
-	  break;
-
-        case 2: {
-	    MidPointAlgorithm jetAlgoM(DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXITERATIONS,DET->M_OVERLAPTHRESHOLD);
-	    jetAlgoM.run(towers, jets);
-	  }
-	  break;
-/*
-	case 3: {
-		FastJet
-  	  }
-	  break;
-*/
-       } // switch
-
-      // Loop on all jets
-      // Dealing with jets, tau-jets, b-jets
-      vector<Cluster>::iterator itJet;
-      for(itJet = jets.begin(); itJet != jets.end(); ++itJet) {
-	  elementJet = (TRootJet*) branchJet->NewEntry();
-          jetMomentum = itJet->fourVector;
-          TLorentzVector JET;
-          JET.SetPxPyPzE(jetMomentum.px,jetMomentum.py,jetMomentum.pz,jetMomentum.E);
-	  elementJet->Set(JET);
-         
-	  // b-jets 
-          bool btag=false;
-          if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
-          elementJet->Btag = btag;
-
-          // Tau jet identification : 1! track and electromagnetic collimation
-	  if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
-	      double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
-	      if( 
-		  ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
-		  ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 )
-                ) {
-		  elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
-		  elementTauJet->Set(JET);
-	       } // if tau jet
-	  } // if JET.eta < tracker - tau_cone : Tau jet identification
-	} // for itJet : loop on all jets
-
+      default:
+      case 1: {
+	
+	// set up a CDF midpoint jet definition
+        #ifdef ENABLE_PLUGIN_CDFCONES
+	plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
+	jet_def = fastjet::JetDefinition(plugins);
+        #else
+	plugins = NULL;
+        #endif
+	
+      } 
+	break;
+	
+      case 2: {
+	
+        // set up a CDF midpoint jet definition
+        #ifdef ENABLE_PLUGIN_CDFCONES
+        plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
+	jet_def = fastjet::JetDefinition(plugins);
+        #else
+	plugins = NULL;
+        #endif
+      }
+	break;
+	
+      case 3: {
+	// set up a siscone jet definition
+        #ifdef ENABLE_PLUGIN_SISCONE
+	int npass = 0;               // do infinite number of passes
+	double protojet_ptmin = 0.0; // use all protojets
+	plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
+	jet_def = fastjet::JetDefinition(plugins);
+        #else
+        plugins = NULL;
+        #endif
+      }
+        break; 
+	
+      case 4: {
+        jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
+        //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
+        //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
+      }
+        break;
+      }
+      // run the jet clustering with the above jet definition
+      if(input_particles.size()!=0)
+	{
+	  fastjet::ClusterSequence clust_seq(input_particles, jet_def);
+	  
+	  
+	  // extract the inclusive jets with pt > 5 GeV
+	  double ptmin = 5.0;
+	  inclusive_jets = clust_seq.inclusive_jets(ptmin);
+	  
+	  // sort jets into increasing pt
+	  sorted_by_pt(inclusive_jets);
+	}
+      
+      for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
+	elementJet = (TRootJet*) branchJet->NewEntry();
+	TLorentzVector JET;
+	JET.SetPxPyPzE(inclusive_jets[i].px(),inclusive_jets[i].py(),inclusive_jets[i].pz(),inclusive_jets[i].E());
+	elementJet->Set(JET);
+	// b-jets
+	bool btag=false;
+	if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
+	elementJet->Btag = btag;
+	
+	// Tau jet identification : 1! track and electromagnetic collimation
+	if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
+	  double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
+	  if(
+	     ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
+	     ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 )
+	     ) {
+	    elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
+	    elementTauJet->Set(JET);
+	  } // if tau jet
+	} // if JET.eta < tracker - tau_cone : Tau jet identification
+      } // for itJet : loop on all jets
+      
       treeWriter->Fill();
       // Add here the trigger
