Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 178)
+++ trunk/Delphes.cpp	(revision 184)
@@ -172,6 +172,7 @@
   TRootCalo *elementCalo;
   
-  TLorentzVector genMomentum(0,0,0,0);
-  TLorentzVector genMomentumCalo(0,0,0,0);
+  TLorentzVector genMomentum(0,0,0,0);        // four-momentum at the vertex
+  TLorentzVector genMomentumBfield(0,0,0,0);  // four-momentum at the exit of the tracks
+  TLorentzVector momentumCaloSegmentation(0,0,0,0);    // four-momentum in the calo, after applying the calo segmentation
   LorentzVector jetMomentum;
   
@@ -235,10 +236,10 @@
             { 
 	      genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+	      genMomentumBfield = genMomentum;
 
 	      // 2.1.2.1 ********************* central detector: magnetic field
-	      // genMomentum is then changed with respect to the magnetic field
-	      //if(DET->FLAG_bfield==1)  genMomentum = TRACP->Propagation(particle);
-	      if(DET->FLAG_bfield==1)  TRACP->Propagation(particle,genMomentum);
-	      float eta=fabs(genMomentum.Eta());
+	      // genMomentumBfield is then changed with respect to the magnetic field
+	      if(DET->FLAG_bfield==1)  TRACP->Propagation(particle,genMomentumBfield);
+	      float eta=fabs(genMomentumBfield.Eta());
 	     
  
@@ -247,27 +248,27 @@
 		
 	      case pE: // all electrons with eta < DET->MAX_CALO_FWD
-		DET->SmearElectron(genMomentum);
-		if(genMomentum.E()!=0 && eta < DET->CEN_max_tracker && genMomentum.Pt() > DET->PTCUT_elec){
-		  electron.push_back(ParticleUtil(genMomentum,particle->PID));
+		DET->SmearElectron(genMomentumBfield);
+		if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_elec){
+		  electron.push_back(ParticleUtil(genMomentumBfield,particle->PID));
 		}
 		break; // case pE
 	      case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-		DET->SmearElectron(genMomentum);
-		if(genMomentum.E()!=0 && eta < DET->CEN_max_tracker && genMomentum.Pt() > DET->PTCUT_gamma) {
-		  gamma.push_back(ParticleUtil(genMomentum,particle->PID));
+		DET->SmearElectron(genMomentumBfield);
+		if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_tracker && genMomentumBfield.Pt() > DET->PTCUT_gamma) {
+		  gamma.push_back(ParticleUtil(genMomentumBfield,particle->PID));
 		}
 		break; // case pGAMMA
 	      case pMU: // all muons with eta < DET->MAX_MU
-		DET->SmearMu(genMomentum);
-		if(genMomentum.E()!=0 && eta < DET->CEN_max_mu &&  genMomentum.Pt() > DET->PTCUT_muon){
-		  muon.push_back(ParticleUtil(genMomentum,particle->PID));
+		DET->SmearMu(genMomentumBfield);
+		if(genMomentumBfield.E()!=0 && eta < DET->CEN_max_mu &&  genMomentumBfield.Pt() > DET->PTCUT_muon){
+		  muon.push_back(ParticleUtil(genMomentumBfield,particle->PID));
 		}
 		break; // case pMU
 	      case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
 	      case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
-		DET->SmearHadron(genMomentum, 0.7);
+		DET->SmearHadron(genMomentumBfield, 0.7);
 		break; // case hadron
 	      default:   // all other final particles with eta < DET->MAX_CALO_FWD
-		DET->SmearHadron(genMomentum, 1.0);
+		DET->SmearHadron(genMomentumBfield, 1.0);
 		break;
 	      } // 2.1.2.2  switch (pid)
@@ -278,15 +279,15 @@
 	      // for calorimetric towers and missing PT
 	      int charge=Charge(pid); 
-	      if(genMomentum.E() !=0 && pid != pMU) {
+	      if(genMomentumBfield.E() !=0 && pid != pMU) {
 		// in case the Bfield is not simulated, checks that charged particles have enough pt to reach the calos
-		if ( !DET->FLAG_bfield && charge!=0 && genMomentum.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
+		if ( !DET->FLAG_bfield && charge!=0 && genMomentumBfield.Pt() <= DET->TRACK_ptmin ) { /* particules do not reach calos */ }
 		else { // particles reach calos	
 	  	  // applies the calo segmentation and returns iEta & iPhi
-		  DET->BinEtaPhi(genMomentum.Phi(), genMomentum.Eta(), iPhi, iEta); 
+		  DET->BinEtaPhi(genMomentumBfield.Phi(), genMomentumBfield.Eta(), iPhi, iEta); 
                   if(iEta != -100 && iPhi != -100) { 
-                      genMomentumCalo.SetPtEtaPhiE(genMomentum.Pt(),iEta,iPhi,genMomentum.E());
+                      momentumCaloSegmentation.SetPtEtaPhiE(genMomentumBfield.Pt(),iEta,iPhi,genMomentumBfield.E());
                       elementCalo = (TRootCalo*) branchCalo->NewEntry();
-                      elementCalo->Set(genMomentumCalo);
-		      PhysicsTower Tower(LorentzVector(genMomentumCalo.Px(),genMomentumCalo.Py(),genMomentumCalo.Pz(),genMomentumCalo.E()));
+                      elementCalo->Set(momentumCaloSegmentation);
+		      PhysicsTower Tower(LorentzVector(momentumCaloSegmentation.Px(),momentumCaloSegmentation.Py(),momentumCaloSegmentation.Pz(),momentumCaloSegmentation.E()));
                       towers.push_back(Tower);
                   } // if iEta != -100
@@ -298,7 +299,7 @@
 	      // all final charged particles
 	      if(
-		 (genMomentum.E()!=0) &&
-		 (fabs(genMomentum.Eta()) < DET->CEN_max_tracker) && 
-		 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentum.Pt() > DET->TRACK_ptmin )) &&     
+		 (genMomentumBfield.E()!=0) &&
+		 (fabs(genMomentumBfield.Eta()) < DET->CEN_max_tracker) && 
+		 (DET->FLAG_bfield || ( !DET->FLAG_bfield && genMomentumBfield.Pt() > DET->TRACK_ptmin )) &&     
 			// if bfield not simulated, pt should be high enough to be taken into account
 		 ((rand()%100) < DET->TRACK_eff)  &&
@@ -307,6 +308,12 @@
 		{
 		  elementTracks = (TRootTracks*) branchTracks->NewEntry();
-		  elementTracks->Set(genMomentum);
-		  TrackCentral.push_back(genMomentum);
+		  elementTracks->Set(genMomentum); // fills px,py,pz,pt,e,eta,phi at vertex
+	 	  elementTracks->Etaout = genMomentumBfield.Eta();
+	 	  elementTracks->Phiout = genMomentumBfield.Phi();
+		  // TODO!!! apply a smearing on the position of the origin of the track
+		  // elementTracks->SetPosition(particle->X,particle->Y,particle->Z);
+		  // uses the output of the bfield computation : Xout=... Yout=... Zout...
+		  // TODO!!! elementTrakcs->SetPositionOut(Xout,Yout,Zout);
+		  TrackCentral.push_back(genMomentum); // tracks at vertex!
 		} // 2.1.2.4 tracks
 	      
@@ -331,5 +338,5 @@
         elementElec->Set(electron[i].Px(),electron[i].Py(),electron[i].Pz(),electron[i].E());
         elementElec->Charge = sign(electron[i].PID());
-        elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);
+        elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);//isolation based on tracks
       }
       DET->SortedVector(muon);
