Index: trunk/src/SmearUtil.cc
===================================================================
--- trunk/src/SmearUtil.cc	(revision 369)
+++ trunk/src/SmearUtil.cc	(revision 374)
@@ -62,4 +62,7 @@
   ELG_Nfwd          = 0.0;               // N term for FCAL
   ELG_Cfwd          = 0.107;             // C term for FCAL
+  ELG_Szdc          = 0.70;              // S term for ZDC
+  ELG_Nzdc          = 0.0;               // N term for ZDC
+  ELG_Czdc          = 0.08;              // C term for ZDC 
 
   // Energy resolution for hadrons in ecal/hcal/hf
@@ -71,7 +74,15 @@
   HAD_Nhf           = 0.;                // N term for FCAL
   HAD_Chf           = 0.13;              // C term for FCAL
+  HAD_Szdc          = 1.38;              // S term for ZDC
+  HAD_Nzdc          = 0.;                // N term for ZDC 
+  HAD_Czdc          = 0.13;              // C term for ZDC 
 
   // Muon smearing
   MU_SmearPt        = 0.01;
+
+  // time resolution
+  ZDC_T_resolution   = 0;	// resolution for time measurement [s]
+  RP220_T_resolution = 0;
+  RP420_T_resolution = 0;
 
   // Tracking efficiencies
@@ -95,5 +106,5 @@
 
 
-  // Thresholds for reconstructed objetcs
+  // Thresholds for reconstructed objetcs (GeV)
   PTCUT_elec      = 10.0;
   PTCUT_muon      = 10.0;
@@ -101,4 +112,7 @@
   PTCUT_gamma     = 10.0;
   PTCUT_taujet    = 10.0;
+
+  ZDC_gamma_E     = 20; // GeV
+  ZDC_n_E	  = 50; // GeV
 
   // Isolation
@@ -199,6 +213,9 @@
   ELG_Sfwd          = DET.ELG_Sfwd;
   ELG_Nfwd          = DET.ELG_Nfwd;
-
-  // Energy resolution for hadrons in ecal/hcal/hf
+  ELG_Czdc          = DET.ELG_Czdc;
+  ELG_Szdc          = DET.ELG_Szdc;
+  ELG_Nzdc          = DET.ELG_Nzdc;
+
+  // Energy resolution for hadrons in ecal/hcal/hf/zdc
   HAD_Shcal         = DET.HAD_Shcal;
   HAD_Nhcal         = DET.HAD_Nhcal;
@@ -207,4 +224,12 @@
   HAD_Nhf           = DET.HAD_Nhf;
   HAD_Chf           = DET.HAD_Chf;
+  HAD_Szdc          = DET.HAD_Szdc;
+  HAD_Nzdc          = DET.HAD_Nzdc;
+  HAD_Czdc          = DET.HAD_Czdc;
+
+  // time resolution
+  ZDC_T_resolution   = DET.ZDC_T_resolution;       // resolution for time measurement [s]
+  RP220_T_resolution = DET.RP220_T_resolution;
+  RP420_T_resolution = DET.RP420_T_resolution;
 
   // Muon smearing
@@ -229,4 +254,7 @@
   PTCUT_gamma     = DET.PTCUT_gamma;
   PTCUT_taujet    = DET.PTCUT_taujet;
+
+  ZDC_gamma_E     = DET.ZDC_gamma_E;
+  ZDC_n_E         = DET.ZDC_n_E;
 
   // Isolation
@@ -321,4 +349,7 @@
   ELG_Sfwd          = DET.ELG_Sfwd;
   ELG_Nfwd          = DET.ELG_Nfwd;
+  ELG_Czdc          = DET.ELG_Czdc;
+  ELG_Szdc          = DET.ELG_Szdc;
+  ELG_Nzdc          = DET.ELG_Nzdc;
 
   // Energy resolution for hadrons in ecal/hcal/hf
@@ -329,4 +360,12 @@
   HAD_Nhf           = DET.HAD_Nhf;
   HAD_Chf           = DET.HAD_Chf;
+  HAD_Szdc          = DET.HAD_Szdc;
+  HAD_Nzdc          = DET.HAD_Nzdc;
+  HAD_Czdc          = DET.HAD_Czdc;
+
+  // time resolution
+  ZDC_T_resolution   = DET.ZDC_T_resolution;       // resolution for time measurement [s]
+  RP220_T_resolution = DET.RP220_T_resolution;
+  RP420_T_resolution = DET.RP420_T_resolution;
 
   // Muon smearing
@@ -351,4 +390,7 @@
   PTCUT_gamma     = DET.PTCUT_gamma;
   PTCUT_taujet    = DET.PTCUT_taujet;
+
+  ZDC_gamma_E     = DET.ZDC_gamma_E;
+  ZDC_n_E         = DET.ZDC_n_E;
 
   // Isolation
@@ -475,4 +517,8 @@
     else if(strstr(temp_string.c_str(),"ELG_Cfwd"))         {curstring >> varname >> value; ELG_Cfwd          = value;}
     else if(strstr(temp_string.c_str(),"ELG_Nfwd"))         {curstring >> varname >> value; ELG_Nfwd          = value;}
+    else if(strstr(temp_string.c_str(),"ELG_Szdc"))         {curstring >> varname >> value; ELG_Szdc          = value;}
+    else if(strstr(temp_string.c_str(),"ELG_Czdc"))         {curstring >> varname >> value; ELG_Czdc          = value;}
+    else if(strstr(temp_string.c_str(),"ELG_Nzdc"))         {curstring >> varname >> value; ELG_Nzdc          = value;}
+
     else if(strstr(temp_string.c_str(),"HAD_Shcal"))        {curstring >> varname >> value; HAD_Shcal         = value;}
     else if(strstr(temp_string.c_str(),"HAD_Nhcal"))        {curstring >> varname >> value; HAD_Nhcal         = value;}
@@ -481,4 +527,10 @@
     else if(strstr(temp_string.c_str(),"HAD_Nhf"))          {curstring >> varname >> value; HAD_Nhf           = value;}
     else if(strstr(temp_string.c_str(),"HAD_Chf"))          {curstring >> varname >> value; HAD_Chf           = value;}
+    else if(strstr(temp_string.c_str(),"HAD_Szdc"))         {curstring >> varname >> value; HAD_Szdc          = value;}
+    else if(strstr(temp_string.c_str(),"HAD_Nzdc"))         {curstring >> varname >> value; HAD_Nzdc          = value;}
+    else if(strstr(temp_string.c_str(),"HAD_Czdc"))         {curstring >> varname >> value; HAD_Czdc          = value;}
+    else if(strstr(temp_string.c_str(),"ZDC_T_resolution")) {curstring >> varname >> value; ZDC_T_resolution  = value;}
+    else if(strstr(temp_string.c_str(),"RP220_T_resolution")) {curstring >> varname >> value; RP220_T_resolution  = value;}
+    else if(strstr(temp_string.c_str(),"RP420_T_resolution")) {curstring >> varname >> value; RP420_T_resolution  = value;}
     else if(strstr(temp_string.c_str(),"MU_SmearPt"))       {curstring >> varname >> value; MU_SmearPt        = value;}
     
@@ -503,4 +555,6 @@
     else if(strstr(temp_string.c_str(),"PTCUT_gamma"))      {curstring >> varname >> value; PTCUT_gamma       = value;}
     else if(strstr(temp_string.c_str(),"PTCUT_taujet"))     {curstring >> varname >> value; PTCUT_taujet      = value;}
+    else if(strstr(temp_string.c_str(),"ZDC_gamma_E"))      {curstring >> varname >> value; ZDC_gamma_E      = value;}
+    else if(strstr(temp_string.c_str(),"ZDC_n_E"))          {curstring >> varname >> value; ZDC_n_E      = value;}
 
     else if(strstr(temp_string.c_str(),"ISOL_PT"))          {curstring >> varname >> value; ISOL_PT           = value;}
@@ -694,4 +748,11 @@
   f_out << left << setw(30) <<"* C term for FCAL: "<<""
 	<< left << setw(30) <<ELG_Cfwd       <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(30) <<"* S term for ZDC: "<<""
+        << left << setw(30) <<ELG_Szdc       <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(30) <<"* N term for ZDC: "<<""
+        << left << setw(30) <<ELG_Nzdc       <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(30) <<"* C term for ZDC: "<<""
+        << left << setw(30) <<ELG_Czdc       <<""<< right << setw(10)<<"*"<<"\n";
+
   f_out<<"*                                                                    *"<<"\n";
   f_out<<"#*****************************                                       *"<<"\n";
@@ -711,4 +772,24 @@
   f_out << left << setw(30) <<"* C term for FCAL: "<<""
 	<< left << setw(30) <<HAD_Chf         <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(30) <<"* S term for ZDC: "<<""
+        << left << setw(30) <<HAD_Szdc        <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(30) <<"* N term for ZDC: "<<""
+        << left << setw(30) <<HAD_Nzdc        <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(30) <<"* C term for ZDC: "<<""
+        << left << setw(30) <<HAD_Czdc        <<""<< right << setw(10)<<"*"<<"\n";
+
+  f_out<<"*                                                                    *"<<"\n";
+  f_out<<"#*************************                                           *"<<"\n";
+  f_out<<"# Time smearing parameters                                           *"<<"\n";
+  f_out<<"#*************************                                           *"<<"\n";
+  f_out<<"*                                                                    *"<<"\n";
+  f_out << left << setw(55) <<"* Time resolution for ZDC        :              "<<""
+        << left << setw(5) <<ZDC_T_resolution       <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(55) <<"* Time resolution for RP220      :              "<<""
+        << left << setw(5) <<RP220_T_resolution     <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(55) <<"* Time resolution for RP420      :              "<<""
+        << left << setw(5) <<RP420_T_resolution     <<""<< right << setw(10)<<"*"<<"\n";
+  f_out<<"*                                                                    *"<<"\n";
+
   f_out<<"*                                                                    *"<<"\n";
   f_out<<"#*************************                                           *"<<"\n";
@@ -786,4 +867,9 @@
   f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
         << left << setw(20) <<PTCUT_gamma        <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(40) <<"* Minimum E for photons in ZDC: "<<""
+        << left << setw(20) <<ZDC_gamma_E        <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(40) <<"* Minimum E for neutrons in ZDC: "<<""
+        << left << setw(20) <<ZDC_n_E        <<""<< right << setw(10)<<"*"<<"\n";
+
   f_out<<"*                                                                    *"<<"\n";
   f_out<<"#*******************                                                 *"<<"\n";
Index: trunk/src/VeryForward.cc
===================================================================
--- trunk/src/VeryForward.cc	(revision 369)
+++ trunk/src/VeryForward.cc	(revision 374)
@@ -4,5 +4,5 @@
 **         |  Delphes, a framework for the fast simulation  |         **
 **         |       of a  generic collider experiment        |         **
-**          \----------------------------------------------/          **
+**          \-------------  arXiv:0903.2225v1  ------------/          **
 **                                                                    **
 **                                                                    **
@@ -39,44 +39,39 @@
 
 //------------------------------------------------------------------------------
-
-VeryForward::VeryForward() {
-   DET = new RESOLution();
-   beamline1 = new H_BeamLine(1,500.);
-   beamline2 = new H_BeamLine(1,500.);
-   init();
-  //Initialisation of Hector
-  relative_energy = true; // should always be true
-  kickers_on = 1;         // should always be 1
-
-}
-
-VeryForward::VeryForward(const string& DetDatacard) {
-   DET = new RESOLution();
+VeryForward::VeryForward() : 
+   DET(new RESOLution()), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)), 
+   beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
+   relative_energy(true), // should always be true
+   kickers_on(1)         // should always be 1
+   {
+   init(); //Initialisation of Hector
+}
+
+VeryForward::VeryForward(const string& DetDatacard) :
+	   DET(new RESOLution())
+   {
    DET->ReadDataCard(DetDatacard);
-   beamline1 = new H_BeamLine(1,500.);
-   beamline2 = new H_BeamLine(1,500.);
-   init();
-  //Initialisation of Hector
-  relative_energy = true; // should always be true
-  kickers_on = 1;         // should always be 1
-
-}
-
-VeryForward::VeryForward(const RESOLution * DetDatacard) {
-   DET = new RESOLution(*DetDatacard);
-   beamline2 = new H_BeamLine(1,500.);
-   beamline1 = new H_BeamLine(1,500.);
-
-   init();
-  //Initialisation of Hector
-  relative_energy = true; // should always be true
-  kickers_on = 1;         // should always be 1
-
-}
-
-VeryForward::VeryForward(const VeryForward& vf) {
-   DET = new RESOLution(*(vf.DET));
-   beamline1 = new H_BeamLine(*(vf.beamline1));
-   beamline2 = new H_BeamLine(*(vf.beamline2));
+   const float d_max = 1.+std::max(DET->RP_420_s,DET->RP_220_s);
+   beamline1 = new H_BeamLine(1,d_max);
+   beamline2 = new H_BeamLine(1,d_max);
+   init(); //Initialisation of Hector
+   relative_energy = true; // should always be true
+   kickers_on = 1;         // should always be 1
+}
+
+VeryForward::VeryForward(const RESOLution * DetDatacard) :
+	DET(new RESOLution(*DetDatacard)), d_max(1.+std::max(DET->RP_420_s,DET->RP_220_s)),
+   	beamline1(new H_BeamLine(1,d_max)), beamline2(new H_BeamLine(1,d_max)),
+        relative_energy(true), // should always be true
+        kickers_on(1)          // should always be 1
+   {
+   init();  //Initialisation of Hector
+}
+
+VeryForward::VeryForward(const VeryForward& vf) :
+   DET(new RESOLution(*(vf.DET))),   d_max(vf.d_max),
+   beamline1(new H_BeamLine(*(vf.beamline1))), beamline2(new H_BeamLine(*(vf.beamline2))),
+   relative_energy(vf.relative_energy),
+   kickers_on(vf.kickers_on) {
 }
 
@@ -84,6 +79,9 @@
    if (this==&vf) return *this;
    DET = new RESOLution(*(vf.DET));
+   d_max = vf.d_max;
    beamline1 = new H_BeamLine(*(vf.beamline1));
    beamline2 = new H_BeamLine(*(vf.beamline2));
+   relative_energy =vf.relative_energy;
+   kickers_on = vf.kickers_on;
    return *this;
 }
@@ -94,6 +92,4 @@
   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->fill(DET->RP_beam1Card,1,DET->RP_IP_name);                                
   beamline1->offsetElements(DET->RP_offsetEl_s,-DET->RP_offsetEl_x);
@@ -114,30 +110,53 @@
 
 
-void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC,TRootGenParticle *particle)
+void VeryForward::ZDC(ExRootTreeWriter *treeWriter, ExRootTreeBranch *branchZDC, TRootGenParticle *particle)
 {
-  int pid=abs(particle->PID);
   TRootZdcHits *elementZdc;
-  TLorentzVector genMomentum; 
+  float energy = particle->E;
+  //TLorentzVector genMomentum; 
+
   // Zero degree calorimeter, for forward neutrons and photons
-  if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && fabs(particle->Eta) > DET->VFD_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 !
+  if (particle->Status ==1 && ( (particle->PID==pN     && energy>DET->ZDC_n_E) || 
+                                (particle->PID==pGAMMA && energy>DET->ZDC_gamma_E) ) 
+			   && fabs(particle->Eta) > DET->VFD_min_zdc ) {
     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->VFD_s_zdc - particle->Z*(1E-6);
-    // assumes also that the emission angle is so small that 1/(cos theta) = 1
-    elementZdc->T = particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
-//cout << "ZDC: T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
+  
+    // 1) energy smearing
+    float energyS = -1.; 
+    if (particle->PID == pGAMMA) 
+        energyS = gRandom->Gaus(particle->E, sqrt( pow(DET->ELG_Nzdc,2) + 
+                                                   pow(DET->ELG_Czdc*particle->E,2) + 
+                                                   pow(DET->ELG_Szdc*sqrt(particle->E),2) ));
+    else // smearing with hadronic resolution
+        energyS =  gRandom->Gaus(particle->E, sqrt( pow(DET->HAD_Nzdc,2) + 
+                                                    pow(DET->HAD_Czdc*particle->E,2) + 
+                                                    pow(DET->HAD_Szdc*sqrt(particle->E),2) ));
+    elementZdc->E = energyS;
+ 
+
+    // 2) time of flight t is t = T + d/[ cos(theta) v ]
+    float cos_theta = 1;      //very good approximation, if eta_zdc >3
+    if (DET->VFD_min_zdc<3) { // if smaller eta -> make the complete calculation
+       double tx = atan(particle->Px/particle->Pz);
+       double ty = atan(particle->Py/particle->Pz);
+       double theta = sqrt( pow(tx,2) + pow(ty,2) );
+       //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
+       // NB: in practice, eta= 8 <-> theta 0.038Â° <-> 7x10^-4 rad <-> cos(theta) ~1
+       // eta = 2.6 <-> cos(theta) = 0.99
+       // eta = 3.0 <-> cos(theta) = 0.995
+       cos_theta = cos(theta);
+    }
+    // units from StdHEP : Z [mm] T[mm/c] 
+    // units from Delphes : VFD_s_zdc [m] speed_of_light [m/s]
+    double flight_distance = (DET->VFD_s_zdc - particle->Z*(1E-3))/cos_theta ;
+    double flight_time = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
+    double timeS = gRandom->Gaus(flight_time,DET->ZDC_T_resolution);
+    elementZdc->T = timeS;
+
+    // 3) side: which ZDC has been hit? 
     elementZdc->side = sign(particle->Eta);
-    
-    
+
+    // 4) object nature : e.m. (photon) or had (neutron) ?
+    elementZdc->hadronic_hit = (bool) (particle->PID==pN);
   }
   
@@ -148,5 +167,5 @@
   
   TRootRomanPotHits* elementRP220;
-  TRootRomanPotHits* elementFP420;
+  TRootForwardTaggerHits* elementFP420;
   
   TLorentzVector genMomentum;
@@ -177,15 +196,31 @@
 	  elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
 	  elementRP220->S = p1.getS();          // [m]
-	  // in first approximation only ! this number is always lower than the real distance-of-flight
-          double flight_distance = p1.getS() - particle->Z*(1E-6); 
-    	  elementRP220->T = particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
-//cout << "T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
-	  elementRP220->E = p1.getE();  	// not yet implemented
-	  elementRP220->q2 = -1;                // not yet implemented
-	  elementRP220->side = sign(particle->Eta);
+
+	    /* time of flight t is t = T + d/[ cos(theta) v ]
+	    // nb: here we assume a straight path to the detector, which is not the case!
+	    // this time estimate is always underestimated (while exact for the ZDC case)
+	    float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
+	    if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
+	       double tx = atan(particle->Px/particle->Pz);
+	       double ty = atan(particle->Py/particle->Pz);
+	       double theta = sqrt( pow(tx,2) + pow(ty,2) );
+	       //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
+	       // NB: in practice, eta= 8 <-> theta 0.038Â° <-> 7x10^-4 rad <-> cos(theta) ~1
+	       // eta = 2.6 <-> cos(theta) = 0.99
+	       // eta = 3.0 <-> cos(theta) = 0.995
+	       cos_theta = cos(theta);
+	    }
+	    // units from StdHEP : Z [mm] T[mm/c]
+	    // units from Delphes : p1.getS [m]   speed_of_light [m/s]
+	    //double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
+	    //elementRP220->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
+	    */
+	    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->RP_420_s);
-	  elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
+	  elementFP420 = (TRootForwardTaggerHits*) branchFP420->NewEntry();
 	  elementFP420->X  = (1E-6)*p1.getX();  // [m]
 	  elementFP420->Y  = (1E-6)*p1.getY();  // [m]
@@ -193,8 +228,23 @@
 	  elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
 	  elementFP420->S = p1.getS();          // [m]
-	  // in first approximation only ! this number is always lower than the real distance-of-flight
-          double flight_distance = p1.getS() - particle->Z*(1E-6);
-//cout << "T = " << particle->T << " ; " << flight_distance/speed_of_light << endl;
-          elementFP420->T = particle->T + flight_distance/speed_of_light;
+
+            // time of flight t is t = T + d/[ cos(theta) v ]
+            // nb: here we assume a straight path to the detector, which is not the case!
+            // this time estimate is always underestimated (while exact for the ZDC case)
+            float cos_theta = 1;      //very good approximation, if CEN_max_calo_fwd >3
+            if (DET->CEN_max_calo_fwd<3) { // if smaller eta -> make the complete calculation
+               double tx = atan(particle->Px/particle->Pz);
+               double ty = atan(particle->Py/particle->Pz);
+               double theta = sqrt( pow(tx,2) + pow(ty,2) );
+               //cout << "tx = " << tx << " ty = " << ty << " theta  = " << theta << " cos(theta) = " << cos(theta) << endl;
+               // NB: in practice, eta= 8 <-> theta 0.038Â° <-> 7x10^-4 rad <-> cos(theta) ~1
+               // eta = 2.6 <-> cos(theta) = 0.99
+               // eta = 3.0 <-> cos(theta) = 0.995
+               cos_theta = cos(theta);
+            }
+            // units from StdHEP : Z [mm] T[mm/c]
+            // units from Delphes : p1.getS [m]   speed_of_light [m/s]
+            double flight_distance = (p1.getS() - particle->Z*(1E-3))/cos_theta ;
+            elementFP420->T = (flight_distance + 1E-3 * particle->T )/speed_of_light; // assumes highly relativistic particles, [s]
 	  elementFP420->E = p1.getE();          // not yet implemented
 	  elementFP420->q2 = -1;                // not yet implemented
