Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 58)
+++ trunk/Delphes.cpp	(revision 59)
@@ -207,22 +207,14 @@
       while( (particle = (TRootGenParticle*) itGen.Next()) )
 	{
-          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
-	  
           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)
-          // with enough p_T
-          //// This subarray is needed for the B-jet algorithm
+         //// This subarray is needed for the B-jet algorithm
 	  // optimization for speed : put first PID condition, then ETA condition, then either pt or status
 	  if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
-              eta < DET->MAX_TRACKER &&
+              fabs(particle->Eta) < DET->MAX_TRACKER &&
               particle->Status != 1 &&
               particle->PT > DET->PT_QUARKS_MIN ) {
             NFCentralQ.Add(particle);
           }
-	  
-	  
+	 	  
 	  // 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
@@ -230,11 +222,11 @@
 	  // 
           if( (particle->Status == 1)    &&
-	      ( 
-	       (pid == pMU  && eta < DET->MAX_MU) ||
-	       (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
-		)
-	      ) { 
-          //  TRACP->Propagation(particle,genMomentum);
-	    eta = fabs(genMomentum.Eta());
+	      ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
+	    ) 
+            { 
+            genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+            TRACP->Propagation(particle,genMomentum);
+            float eta=fabs(genMomentum.Eta());
             switch(pid) {
 	      
@@ -264,5 +256,5 @@
 	      break;
 	    } // switch (pid)
-	    
+
 	    // all final particles but muons and neutrinos	
 	    // for calorimetric towers and mission PT
@@ -275,5 +267,5 @@
 		// back of the input_particles vector
 		input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
-                //genMomentumCalo.SetPtEtaPhiE(CaloTower.Et(),CaloTower.iEta(),CaloTower.iPhi(),CaloTower.fourVector.E);
+                
                 genMomentumCalo.SetPxPyPzE(CaloTower.fourVector.px,CaloTower.fourVector.py,CaloTower.fourVector.pz,CaloTower.fourVector.E);
 		elementCalo = (TRootCalo*) branchCalo->NewEntry();
@@ -301,5 +293,5 @@
 		TrackCentral.push_back(genMomentum);
 	      }
-	    
+
 	  } // switch
 	  
@@ -340,5 +332,7 @@
       elementEtmis->Px = (-PTmis).Px();
       elementEtmis->Py = (-PTmis).Py();
+
       //*****************************
+      treeWriter->Fill();
       
       sorted_jets=JETRUN->RunJets(input_particles);
@@ -346,12 +340,6 @@
       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
Index: trunk/src/BFieldProp.cc
===================================================================
--- trunk/src/BFieldProp.cc	(revision 58)
+++ trunk/src/BFieldProp.cc	(revision 59)
@@ -45,5 +45,5 @@
   
   //out of trackibg coverage?
-  if(sqrt(pow(Xvertex1,2)+pow(Yvertex1,2)) > TRACKING_RADIUS){return;}
+  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > TRACKING_RADIUS){return;}
   if(fabs(Zvertex1) > TRACKING_LENGTH){return;}
   
@@ -55,79 +55,69 @@
   double py = Py / 0.003;
   double pz = Pz / 0.003;
-  double pt = sqrt(pow(px,2)+pow(py,2));
+  double pt = sqrt(px*px+py*py);
+  double p  = sqrt(px*px+py*py+pz*pz);
+ 
+  if(q!=0){
+     double e  = Part->E  / 0.003;
+//     double M  = sqrt(e*e - (px*px + py*py + pz*pz) );
+     double M  = Part->M;
+/*if(fabs(Part->PID)==11)
+{
+     cout<<"genMomentum.M() "<<genMomentum.M()<<"  "<<Part->M<<endl;
+     cout<<"e*e - (px*px + py*py + pz*pz) "<<e*e - (px*px + py*py + pz*pz)<<endl;
+}*/
+     double vx = px/M;
+     double vy = py/M;
+     double vz = pz/M;
 
-  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;
-	}	
+     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);
+/*
+cout<<"M "<<M<<" et le pid "<<Part->PID<<endl;
+cout<<"ax "<<ax<<" ay "<<ay<<" az "<<az<<endl;
+*/
+     double dt = 1/p;
+     if(pt<266 && vz < 0.0012)       dt = fabs(0.001/vz);
+ 
+     double xold=Xvertex1;         double x=xold;
+     double yold=Yvertex1;         double y=yold;
+     double zold=Zvertex1;         double z=zold;
+
+     double VTold = sqrt(vx*vx+vy*vy);
+ 
+     int k = 0;
+     
+     while(k < MAXITERATION){
+        k++;
+
+        vx += ax*dt;
+        vy += ay*dt;
+        vz += az*dt;
+
+        double VTratio = VTold/sqrt(vx*vx+vy*vy);
+        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;
+
+       xold = x;
+       yold = y;
+       zold = z;
+     }
+ 
       if(x!=0 && y!=0 && z!=0)
 	{
@@ -139,5 +129,5 @@
 	  genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
 	}
-      cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
+     // cout<<"genMomentum final "<<genMomentum.Pt()<<endl;
     }
 }
