Index: trunk/Resolutions.cpp
===================================================================
--- trunk/Resolutions.cpp	(revision 87)
+++ trunk/Resolutions.cpp	(revision 88)
@@ -29,18 +29,9 @@
 
 #include "interface/SmearUtil.h"
+#include "interface/JetUtils.h"
+#include "interface/BFieldProp.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>
@@ -109,5 +100,4 @@
 }
 
-
 int main(int argc, char *argv[])
 {
@@ -133,4 +123,5 @@
     return -1;
   }
+
   TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
@@ -196,54 +187,23 @@
   RESOLution *DET = new RESOLution();   
   DET->ReadDataCard(DetDatacard);
-  
-  TLorentzVector genMomentum(0,0,0,0);
+
+  //Jet information
+  JetsUtil *JETRUN = new JetsUtil();
+
+  //Propagation of tracks in the B field
+  TrackPropagation *TRACP = new TrackPropagation();
+   
+  TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
+  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing Reco level information
   LorentzVector jetMomentum;
   vector<TLorentzVector> TrackCentral;
   
-  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
-  vector<fastjet::PseudoJet> inclusive_jets;
-  vector<fastjet::PseudoJet> sorted_jets;
-  
-  vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
-  vector<fastjet::PseudoJet> inclusive_jetsS;
-  vector<fastjet::PseudoJet> sorted_jetsS;
-
-  vector<fastjet::PseudoJet> input_particlesT;//for FastJet algorithm
-  vector<fastjet::PseudoJet> inclusive_jetsT;
-  vector<fastjet::PseudoJet> sorted_jetsT;
+  vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jetsGEN;
+  
+  vector<fastjet::PseudoJet> input_particlesReco;//for FastJet algorithm
+  vector<fastjet::PseudoJet> sorted_jetsReco;
 
   vector<PhysicsTower> towers;
-  
-  fastjet::JetDefinition jet_def;
-  fastjet::JetDefinition jet_defS;
-  fastjet::JetDefinition jet_defT;
-  fastjet::JetDefinition::Plugin * plugins;
-  fastjet::JetDefinition::Plugin * pluginsS;
-  fastjet::JetDefinition::Plugin * pluginsT;
-  
-  // set up a CDF midpoint jet definition
-#ifdef ENABLE_PLUGIN_CDFCONES
-  plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
-  jet_def = fastjet::JetDefinition(plugins);
-#else
-  plugins = NULL;
-#endif
-  
-  // set up a CDF midpoint jet definition
-#ifdef ENABLE_PLUGIN_CDFCONES
-  pluginsS = new fastjet::CDFJetCluPlugin(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
-  jet_defS = fastjet::JetDefinition(pluginsS);
-#else
-  pluginsS = NULL;
-#endif
-  
- // set up a CDF midpoint jet definition
- #ifdef ENABLE_PLUGIN_CDFCONES
- pluginsT = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->OVERLAPTHRESHOLD);
- jet_defT = fastjet::JetDefinition(pluginsT);
- #else
- pluginsT = NULL;
- #endif
- 
   
   // Loop over all events
@@ -252,6 +212,6 @@
   for(entry = 0; entry < allEntries; ++entry)
     {
-      TLorentzVector PTmisS(0,0,0,0);
-      TLorentzVector PTmis(0,0,0,0);
+      TLorentzVector PTmisReco(0,0,0,0);
+      TLorentzVector PTmisGEN(0,0,0,0);
       treeReader->ReadEntry(entry);
       treeWriter->Clear();
@@ -263,15 +223,6 @@
       TSimpleArray<TRootGenParticle> NFCentralQ;
    
-      sorted_jetsS.clear();
-      input_particlesS.clear();
-      inclusive_jetsS.clear();
-
-      sorted_jetsT.clear();
-      input_particlesT.clear();
-      inclusive_jetsT.clear();
-
-      sorted_jets.clear();
-      input_particles.clear();
-      inclusive_jets.clear();
+      input_particlesReco.clear();
+      input_particlesGEN.clear();
       towers.clear();
       
@@ -283,118 +234,107 @@
 	  int pid = abs(particle->PID);
 	  float eta = fabs(particle->Eta);
-	  
+
+          //input generator level particle for jet algorithm	  
           if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
 	    {
-	      input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+	      input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
 	    }
-	  
-          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
-	  // 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)    &&
-	      ( 
-	       (pid == pMU  && eta < DET->MAX_MU) ||
-	       (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) )
-	      ) 
-	    {
-              //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
-
-              Float_t nonS=0;
+
+	  //Calculate ETMIS from generated particles
+          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
+
+	  if( (particle->Status == 1)    &&
+              ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
+              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
+	      )
+            {
+              recoMomentum = genMomentum;          
+              //use of the magnetic field propagation
+              TRACP->Propagation(particle,recoMomentum);
+              float eta=fabs(recoMomentum.Eta());
+
 	      switch(pid) {
+
 	      case pE: // all electrons with eta < DET->MAX_CALO_FWD
-                nonS = genMomentum.E();
-                DET->SmearElectron(genMomentum);
+                DET->SmearElectron(recoMomentum);
                 if(eta < DET->MAX_TRACKER){
                 elementElec=(RESOLELEC*) branchelec->NewEntry();
-                elementElec->NonSmeareE = nonS;
-                elementElec->SmeareE = genMomentum.E();}
+                elementElec->E = genMomentum.E();
+                elementElec->SmearedE = recoMomentum.E();}
 		break; // case pE
 	      case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-		DET->SmearElectron(genMomentum);
+		DET->SmearElectron(recoMomentum);
 		break; // case pGAMMA
 	      case pMU: // all muons with eta < DET->MAX_MU
+		DET->SmearMu(recoMomentum);
                 elementMuon = (RESOLMUON*) branchmuon->NewEntry();
-                elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
-		DET->SmearMu(genMomentum);
-                elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
+                elementMuon->OverPT = (1/genMomentum.Pt());
+                elementMuon->OverSmearedPT = (1/recoMomentum.Pt());
 		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(recoMomentum, 0.7);
 		break; // case hadron
 	      default:   // all other final particles with eta < DET->MAX_CALO_FWD
-		DET->SmearHadron(genMomentum, 1.0);
+		DET->SmearHadron(recoMomentum, 1.0);
 		break;
 	      } // switch (pid)
-	      
-	      if(pid != pMU && genMomentum.Et() !=0)
-		{
-		  towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
-		  input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+	     
+            //information to reconstruct jets from reconstructed towers
+            int charge=Charge(pid);
+            if(recoMomentum.E() !=0 && pid != pMU) {
+              if(charge == 0 || (charge !=0 && recoMomentum.Pt() >= DET->PT_TRACKS_MIN)){
+		  towers.push_back(PhysicsTower(LorentzVector(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E())));
+		  input_particlesReco.push_back(fastjet::PseudoJet(recoMomentum.Px(),recoMomentum.Py(),recoMomentum.Pz(), recoMomentum.E()));
 		}
+            }
 	      
 	      // all final charged particles
-	      if(
-		 ((rand()%100) < DET->TRACKING_EFF)  &&
-		 (genMomentum.E()!=0) &&
-		 (fabs(particle->Eta) < DET->MAX_TRACKER) && 
-		 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
-		 (pid != pGAMMA)    && 
-		 (pid != pPI0)      &&
-		 (pid != pK0L)      &&
-		 (pid != pN)        &&
-		 (pid != pSIGMA0)   &&
-		 (pid != pDELTA0)   &&
-		 (pid != pK0S)   // not charged particles : invisible by tracker
-                 )
-		TrackCentral.push_back(genMomentum);
+	         if(
+	           (genMomentum.E()!=0) &&
+	           (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) &&
+	           (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
+	           ((rand()%100) < DET->TRACKING_EFF)  &&
+	           (charge!=0)
+	           )
+	           {TrackCentral.push_back(genMomentum);}
+
 	    } // switch
 	} // while
-    
+
+     //compute missing transverse energy from calo towers    
      TLorentzVector Att(0.,0.,0.,0.);
      for(unsigned int i=0; i < towers.size(); i++) 
        {
          Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
-         PTmisS = PTmisS + Att;
+         PTmisReco = PTmisReco + Att;
        }
 
       elementEtmis= (ETMIS*) branchetmis->NewEntry();
-      elementEtmis->Et = (PTmis).Pt();
-      elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
+      elementEtmis->Et = (PTmisGEN).Pt();
+      elementEtmis->SmearedEt = (PTmisReco).Pt()-(PTmisGEN).Pt();
       
       //*****************************
       
-      double ptmin=1;
-      if(input_particles.size()!=0)
-        {
-          fastjet::ClusterSequence clust_seq(input_particles, jet_def);
-          inclusive_jets = clust_seq.inclusive_jets(ptmin);
-          sorted_jets = sorted_by_pt(inclusive_jets);
-        }
-      
-      if(input_particlesS.size()!=0)
-        {
-          fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
-          inclusive_jetsS = clust_seqS.inclusive_jets(ptmin);
-          sorted_jetsS = sorted_by_pt(inclusive_jetsS);
-        }
-      
+      sorted_jetsGEN=JETRUN->RunJets(input_particlesGEN);
+      sorted_jetsReco=JETRUN->RunJets(input_particlesReco);
+
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
 
-      TLorentzVector JETSm(0,0,0,0);
-      for (unsigned int i = 0; i < sorted_jets.size(); i++) {
-        TLorentzVector JET(0,0,0,0);
-        JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
-        PairingJet(JETSm,JET,sorted_jetsS);
-        if(JETSm.Pt()>1)
+      TLorentzVector JETreco(0,0,0,0);
+      for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
+        TLorentzVector JETgen(0,0,0,0);
+        JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
+        PairingJet(JETreco,JETgen,sorted_jetsReco);
+        if(JETreco.Pt()>1)
           {
             elementJet= (RESOLJET*) branchjet->NewEntry();
-            elementJet->NonSmearePT = JET.Et();
-            elementJet->SmearePT = JETSm.Et()/JET.Et();
+            elementJet->PT = JETgen.Et();
+            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
           }
       }
-      for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
+      for (unsigned int i = 0; i < sorted_jetsReco.size(); i++) {
         TLorentzVector JETT(0,0,0,0);
-        JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
+        JETT.SetPxPyPzE(sorted_jetsReco[i].px(),sorted_jetsReco[i].py(),sorted_jetsReco[i].pz(),sorted_jetsReco[i].E());
         if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) 
            {
Index: trunk/interface/TreeClasses.h
===================================================================
--- trunk/interface/TreeClasses.h	(revision 87)
+++ trunk/interface/TreeClasses.h	(revision 88)
@@ -13,6 +13,6 @@
 {
 public:
-Float_t SmearePT;
-Float_t NonSmearePT;
+Float_t SmearedPT;
+Float_t PT;
 
   static TCompare *fgCompare; //!
@@ -30,6 +30,6 @@
 {
 public:
-Float_t SmeareE;
-Float_t NonSmeareE;
+Float_t E;
+Float_t SmearedE;
 
   static TCompare *fgCompare; //!
@@ -46,6 +46,6 @@
 {
 public:
-Float_t OneoverSmearePT;
-Float_t OneoverNonSmearePT;
+Float_t OverSmearedPT;
+Float_t OverPT;
 
   static TCompare *fgCompare; //!
@@ -77,5 +77,5 @@
 public:
 Float_t Et;
-Float_t EtSmeare;
+Float_t SmearedEt;
 
  static TCompare *fgCompare; //!
