Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 49)
+++ trunk/Delphes.cpp	(revision 55)
@@ -19,8 +19,4 @@
 #include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
 
-#include "H_BeamParticle.h"
-#include "H_BeamLine.h"
-#include "H_RomanPot.h"
-
 #include "interface/DataConverter.h"
 #include "interface/HEPEVTConverter.h"
@@ -29,23 +25,11 @@
 
 #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>
-
-
+#include "interface/BFieldProp.h"
+#include "interface/TriggerUtil.h"
+#include "interface/VeryForward.h"
+#include "interface/JetUtils.h"
+
+#include <vector>
+#include <iostream>
 
 using namespace std;
@@ -68,5 +52,5 @@
 {
   int appargc = 2;
-  char *appName = "JetsSim";
+  char *appName = "Delphes";
   char *appargv[] = {appName, "-b"};
   TApplication app(appName, &appargc, appargv);
@@ -103,9 +87,24 @@
   string DetDatacard("");
   if(argc==4)  DetDatacard =argv[3];
+
+  //Smearing information
   RESOLution *DET = new RESOLution();
   DET->ReadDataCard(DetDatacard);
   DET->Logfile(LogName);
 
-  todo(LogName.c_str());
+  //Trigger information
+  Trigger *TRIG = new Trigger();
+  TRIG->TriggerReader("data/trigger.dat");
+
+  //Propagation of tracks in the B field
+  TrackPropagation *TRACP = new TrackPropagation(); 
+
+  //Jet information
+  JetsUtil *JETRUN = new JetsUtil();
+
+  //VFD information
+  VeryForward * VFD = new VeryForward();
+
+  //todo(LogName.c_str());
  
   DataConverter *converter=0;
@@ -163,19 +162,16 @@
   TRootMuon *elementMu;
   TRootPhoton *elementPhoton;
-  TRootJet *elementJet;
-  TRootTauJet *elementTauJet;
   TRootTracks *elementTracks;
   TRootCalo *elementCalo;
-  TRootZdcHits *elementZdc;
-  TRootRomanPotHits *elementRP220, *elementFP420;
   
   TLorentzVector genMomentum(0,0,0,0);
   TLorentzVector genMomentumCalo(0,0,0,0);
   LorentzVector jetMomentum;
+ 
+  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jets;
+
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
-  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
-  vector<fastjet::PseudoJet> inclusive_jets;
-  vector<fastjet::PseudoJet> sorted_jets;
   
   vector<TLorentzVector> electron;
@@ -185,85 +181,6 @@
   TSimpleArray<TRootGenParticle> NFCentralQ;
   
-  //Initialisation of Hector
-  extern bool relative_energy;
-  relative_energy = true; // should always be true
-  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
-  H_BeamLine* beamline1 = new H_BeamLine(1,500.);
-  beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
-  beamline1->offsetElements(120,-0.097);
-  H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",220,2000); // RP 220m, 2mm, beam 1
-  H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",420,4000); // RP 420m, 4mm, beam 1
-  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");
-  beamline2->offsetElements(120,+0.097);
-  H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",220,2000);// RP 220m, 2mm, beam 2
-  H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",420,4000);// RP 420m, 4mm, beam 2
-  beamline2->add(rp220_2);
-  beamline2->add(rp420_2);
-  
-  // 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: {
-    
-    // set up a CDF midpoint jet definition
-#ifdef ENABLE_PLUGIN_CDFCONES
-    plugins = new fastjet::CDFJetCluPlugin(DET->SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->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->SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXITERATIONS,DET->OVERLAPTHRESHOLD);
-    jet_def = fastjet::JetDefinition(plugins);
-#else
-    plugins = NULL;
-#endif
-  }
-    break;
-  case 3: {
-    // set up a siscone jet definition
-#ifdef ENABLE_PLUGIN_SISCONE
-    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->OVERLAPTHRESHOLD,DET->NPASS, DET->PROTOJET_PTMIN);
-    jet_def = fastjet::JetDefinition(plugins);
-#else
-    plugins = NULL;
-#endif
-  }
-    break;
-    
-  case 4: {
-    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
-  }
-    break;
-    
-  case 5: {
-    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,DET->CONERADIUS);
-  }
-    break;
-    
-  case 6: {
-    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,DET->CONERADIUS);
-  }
-    break;
-  }
-  
+
+
   // Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
@@ -286,6 +203,4 @@
       towers.clear();
       input_particles.clear();
-      inclusive_jets.clear();
-      sorted_jets.clear();
       
       // Loop over all particles in event
@@ -294,7 +209,7 @@
           genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
 	  
-	  int pid = abs(particle->PID);
-	  float eta = fabs(particle->Eta);
-	  
+          int pid = abs(particle->PID);
+          float eta = fabs(genMomentum.Eta());
+
           //subarray of the quarks (i.e. not final particles, i.e status not equal to 1) 
           // in the tracker (i.e. for those we know the tracks)
@@ -312,4 +227,6 @@
 	  // keeps only final particles, visible by the central detector, including the fiducial volume
 	  // the ordering of conditions have been optimised for speed : put first the STATUS condition
+	  //
+	  // 
           if( (particle->Status == 1)    &&
 	      ( 
@@ -318,5 +235,7 @@
 		)
 	      ) { 
-	    switch(pid) {
+          //  TRACP->Propagation(particle,genMomentum);
+	    eta = fabs(genMomentum.Eta());
+            switch(pid) {
 	      
 	    case pE: // all electrons with eta < DET->MAX_CALO_FWD
@@ -363,10 +282,9 @@
 	    }
 	    
-	    
 	    // all final charged particles
 	    if(
 	       ((rand()%100) < DET->TRACKING_EFF)  &&
 	       (genMomentum.E()!=0) &&
-	       (fabs(particle->Eta) < DET->MAX_TRACKER) && 
+	       (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) && 
 	       (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
 	       (pid != pGAMMA)    && 
@@ -386,79 +304,6 @@
 	  } // switch
 	  
-	  // Forward particles in 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);
-	    
-	  }
-	  
-	  // 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);
-		}
-	      }
-	      
-	      //		if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
-	      //			cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
-	    } // if forward proton
+          VFD->ZDC(treeWriter,branchZDC,particle);
+          VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
 	  
 	} // while
@@ -497,45 +342,19 @@
       //*****************************
       
-      // 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_jets = sorted_by_pt(inclusive_jets);
-	}
-      for (unsigned int i = 0; i < sorted_jets.size(); i++) {
-	TLorentzVector JET;
-	JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
-        // 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 ) &&
-	     ( JET.Pt() > DET->TAUJET_pt)
-	     ) {
-	    elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
-	    elementTauJet->Set(JET);
-	  } // if tau jet
-	} // if JET.eta < tracker - tau_cone : Tau jet identification
-	
-        if(JET.Pt() > DET->JET_pt)
-          {
-            elementJet = (TRootJet*) branchJet->NewEntry();
-	    elementJet->Set(JET);
-	    // b-jets
-	    bool btag=false;
-	    if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
-	    elementJet->Btag = btag;
-	  }
-      } // for itJet : loop on all jets
-      
+      sorted_jets=JETRUN->RunJets(input_particles);
+      JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
+      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
+
+
       treeWriter->Fill();
+
       // Add here the trigger
       // Should test all the trigger table on the event, based on reconstructed objects
+      
+ //     Trigger *TRIG = new Trigger();
+ //     TRIG->TriggerReader("data/trigger.dat");
+     
     } // Loop over all events
+
   treeWriter->Write();
   
