Index: trunk/src/BFieldProp.cc
===================================================================
--- trunk/src/BFieldProp.cc	(revision 53)
+++ trunk/src/BFieldProp.cc	(revision 53)
@@ -0,0 +1,143 @@
+/*
+ *   ---- Delphes ----
+ *     A Fast Simulator for general purpose LHC detector
+ *       S. Ovyn ~~~~ severine.ovyn@uclouvain.be
+ *
+ *     Center for Particle Physics and Phenomenology (CP3)
+ *     Universite Catholique de Louvain (UCL)
+ *     Louvain-la-Neuve, Belgium
+ *             */
+
+#include "interface/BFieldProp.h"
+#include "TRandom.h"
+
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <iomanip>
+
+#include<cmath>
+
+
+using namespace std;
+
+
+//------------------------------------------------------------------------------
+
+TrackPropagation::TrackPropagation() {
+
+ TRACKING_RADIUS = 129;
+ TRACKING_LENGTH = 300;
+ MAXITERATION = 20000;
+ MINSEGLENGTH = 70;
+
+}
+
+void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum)
+{
+
+  double q  = Charge(Part->PID);
+  if(q==0)return;
+  
+  float Xvertex1 = Part->X;
+  float Yvertex1 = Part->Y;
+  float Zvertex1 = Part->Z;
+  
+  //out of trackibg coverage?
+  if(sqrt(pow(Xvertex1,2)+pow(Yvertex1,2)) > TRACKING_RADIUS){return;}
+  if(fabs(Zvertex1) > TRACKING_LENGTH){return;}
+  
+  float Px = Part->Px;
+  float Py = Part->Py;
+  float Pz = Part->Pz;
+    
+  double px = Px / 0.003;
+  double py = Py / 0.003;
+  double pz = Pz / 0.003;
+  double pt = sqrt(pow(px,2)+pow(py,2));
+
+  double p  = sqrt(pow(px,2)+pow(py,2)+pow(pz,2));
+  
+  double m = sqrt(pow(Part->E,2) -  pow(Part->Px,2)+ pow(Part->Py,2)+ pow(Part->Pz,2));
+  
+  if(q!=0)
+    {
+      double e  = Part->E  / 0.003;
+      double M  = sqrt(e*e - (px*px + py*py + pz*pz) );
+      double vx = px/M;
+      double vy = py/M;
+      double vz = pz/M;
+      
+      /*double Bx = BField_[0];
+	double By = BField_[1];
+	double Bz = BField_[2];
+      */
+      
+      double Bx = 0;
+      double By = 0;
+      double Bz = 3.8;
+      
+      double ax =  (q/M)*(Bz*vy - By*vz);
+      double ay =  (q/M)*(Bx*vz - Bz*vx);
+      double az =  (q/M)*(By*vx - Bx*vy);
+      
+      double xold=Xvertex1;        
+      double yold=Yvertex1;        
+      double zold=Zvertex1;       
+      
+      double dt = 1/p;
+            
+      double x = xold;
+      double y = yold;        
+      double z = zold;  
+      
+      double VTold = sqrt(pow(vx,2)+pow(vy,2));
+      
+      int k = 0;
+      
+      while(k < MAXITERATION)
+	{
+	  k++;
+	  
+	  //
+	  vx += ax*dt;
+	  vy += ay*dt;
+	  vz += az*dt;
+	  
+	  double VTratio = VTold/sqrt(pow(vx,2)+pow(vy,2));
+	  vx *= VTratio;
+	  vy *= VTratio;
+	  
+	  ax = (q/M)*(Bz*vy - By*vz);
+	  ay = (q/M)*(Bx*vz - Bz*vx);
+	  az = (q/M)*(By*vx - Bx*vy);
+	  
+	  x  += vx*dt;
+	  y  += vy*dt;
+	  z  += vz*dt;
+	  
+	  if( (x*x+y*y) > (TRACKING_RADIUS*TRACKING_RADIUS) )
+	    {
+	      x /= (x*x+y*y)/ (TRACKING_RADIUS*TRACKING_RADIUS); 
+	      y /= (x*x+y*y)/ (TRACKING_RADIUS*TRACKING_RADIUS); 
+	      break;
+	    }
+	  if( fabs(z) > TRACKING_LENGTH)break;
+	  if( fabs((pow(x,2)+pow(y,2)+pow(z,2)) - (pow(xold,2)+pow(yold,2)+pow(zold,2))) < MINSEGLENGTH )continue;
+	  
+	  xold = x;
+	  yold = y;
+	  zold = z;
+	}	
+      if(x!=0 && y!=0 && z!=0)
+	{
+	  double theta = fabs(atan(y/z));
+	  double eta;
+	  if(z > 0)eta = -log(tan(theta/2));
+	  else eta = -(-log(tan(theta/2)));
+	  double phi = atan2(y,x);
+	  genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
+	}
+      cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
+    }
+}
Index: trunk/src/JetUtils.cc
===================================================================
--- trunk/src/JetUtils.cc	(revision 53)
+++ trunk/src/JetUtils.cc	(revision 53)
@@ -0,0 +1,136 @@
+/*
+ *   ---- Delphes ----
+ *     A Fast Simulator for general purpose LHC detector
+ *       S. Ovyn ~~~~ severine.ovyn@uclouvain.be
+ *
+ *         Center for Particle Physics and Phenomenology (CP3)
+ *           Universite Catholique de Louvain (UCL)
+ *             Louvain-la-Neuve, Belgium
+ *             */
+
+// \brief Trigger class, and some generic definitions
+
+#include "interface/JetUtils.h"
+
+using namespace std;
+
+JetsUtil::JetsUtil() {
+  
+  switch(JETALGO) {
+  default:
+  case 1: {
+    
+    // set up a CDF midpoint jet definition
+#ifdef ENABLE_PLUGIN_CDFCONES
+    plugins = new fastjet::CDFJetCluPlugin(SEEDTHRESHOLD,CONERADIUS,C_ADJACENCYCUT,C_MAXITERATIONS,C_IRATCH,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 (SEEDTHRESHOLD,CONERADIUS,M_CONEAREAFRACTION,M_MAXPAIRSIZE,M_MAXITERATIONS,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 (CONERADIUS,OVERLAPTHRESHOLD,NPASS, PROTOJET_PTMIN);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+    
+  case 4: {
+    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, CONERADIUS);
+  }
+    break;
+    
+  case 5: {
+    jet_def = fastjet::JetDefinition(fastjet::cambridge_algorithm,CONERADIUS);
+  }
+    break;
+    
+  case 6: {
+    jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm,CONERADIUS);
+  }
+    break;
+  }
+  
+}
+
+vector<fastjet::PseudoJet> JetsUtil::RunJets(const vector<fastjet::PseudoJet> & input_particles)
+{
+  inclusive_jets.clear();
+  sorted_jets.clear();
+  
+  // 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);
+    }
+  
+  return sorted_jets;                                                    
+}
+
+void JetsUtil::RunJetBtagging(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchJet,const vector<fastjet::PseudoJet> & sorted_jets,const TSimpleArray<TRootGenParticle>& NFCentralQ)
+{
+  TRootJet *elementJet;
+  TLorentzVector JET;
+  for (unsigned int i = 0; i < sorted_jets.size(); i++) {
+    JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
+    if(JET.Pt() > JET_pt)
+      {
+	elementJet = (TRootJet*) branchJet->NewEntry();
+	elementJet->Set(JET);
+	// b-jets
+	bool btag=false;
+	if((fabs(JET.Eta()) < MAX_TRACKER && Btaggedjet(JET, NFCentralQ)))btag=true;
+	elementJet->Btag = btag;
+      }
+  } // for itJet : loop on all jets
+ treeWriter->Fill();
+  
+}
+
+void JetsUtil::RunTauJets(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchTauJet,const vector<fastjet::PseudoJet> & sorted_jets,const vector<PhysicsTower> & towers, const vector<TLorentzVector> & TrackCentral)
+{
+  TRootTauJet *elementTauJet;
+  TLorentzVector JET;
+  for (unsigned int i = 0; i < sorted_jets.size(); i++) {
+    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()) < (MAX_TRACKER - TAU_CONE_TRACKS)) {
+      double Energie_tau_central = EnergySmallCone(towers,JET.Eta(),JET.Phi());
+      if(
+	 ( Energie_tau_central/JET.E() > TAU_EM_COLLIMATION ) &&
+	 ( NumTracks(TrackCentral,PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 ) &&
+	 ( JET.Pt() > TAUJET_pt)
+	 ) {
+	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();
+  
+}
+
Index: trunk/src/VeryForward.cc
===================================================================
--- trunk/src/VeryForward.cc	(revision 53)
+++ trunk/src/VeryForward.cc	(revision 53)
@@ -0,0 +1,155 @@
+ /*
+ *   ---- Delphes ----
+ *     A Fast Simulator for general purpose LHC detector
+ *       S. Ovyn ~~~~ severine.ovyn@uclouvain.be
+ *
+ *     Center for Particle Physics and Phenomenology (CP3)
+ *     Universite Catholique de Louvain (UCL)
+ *     Louvain-la-Neuve, Belgium
+ *             */
+
+#include "interface/VeryForward.h"
+#include "interface/SmearUtil.h"
+#include "TRandom.h"
+
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <iomanip>
+
+#include<cmath>
+
+
+using namespace std;
+
+
+//------------------------------------------------------------------------------
+
+VeryForward::VeryForward() {
+  
+  //Initialisation of Hector
+  relative_energy = true; // should always be true
+  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
+  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);
+  
+  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); 
+  
+  
+}
+
+void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC,TRootGenParticle *particle)
+{
+  int pid=abs(particle->PID);
+  float eta=fabs(particle->Eta);
+  
+  
+  TRootZdcHits *elementZdc;
+  TLorentzVector genMomentum; 
+  // Zero degree calorimeter, for forward neutrons and photons
+  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > MIN_ZDC ) {
+    genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+    // !!!!!!!!! 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 = (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(eta);
+    
+    treeWriter->Fill();
+    
+  }
+  
+}
+void VeryForward::RomanPots(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchRP220,ExRootTreeBranch *branchFP420,TRootGenParticle *particle) 
+{
+  int pid=abs(particle->PID);
+  float eta=fabs(particle->Eta);
+  
+  TRootRomanPotHits* elementRP220;
+  TRootRomanPotHits* elementFP420;
+  
+  TLorentzVector genMomentum;
+  genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+  // if forward proton
+  if( (pid == pP)  && (particle->Status == 1) &&  (fabs(genMomentum.Eta()) > 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(genMomentum.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(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(eta);
+	  
+	} else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
+	  p1.propagate(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(eta);
+	}
+
+      }
+     treeWriter->Fill();
+      //                if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
+      //                        cout << "Eloss ="  << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
+    } // if forward proton
+}
+
+    // Forward particles in CASTOR ?
+    //           /*      if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
+    //                             && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
+    //
+    //
+    //                                               } // CASTOR
+    //                                                         */
+    //
