Index: trunk/Resolutions.cpp
===================================================================
--- trunk/Resolutions.cpp	(revision 25)
+++ trunk/Resolutions.cpp	(revision 26)
@@ -52,5 +52,7 @@
 //------------------------------------------------------------------------------
 
-void PairingJet(TLorentzVector &JETSm, const TLorentzVector& JET,  vector<fastjet::PseudoJet> sorted_jetsS)
+
+
+void PairingJet(TLorentzVector &JETSm, TLorentzVector JET,  vector<fastjet::PseudoJet> sorted_jetsS)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
@@ -70,10 +72,9 @@
 }
 
-//------------------------------------------------------------------------------
 
 int main(int argc, char *argv[])
 {
   int appargc = 2;
-  char appName[100]; sprintf(appName,argv[0]); 
+  char *appName = "Smearing";
   char *appargv[] = {appName, "-b"};
   TApplication app(appName, &appargc, appargv);
@@ -97,11 +98,11 @@
   TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
-
+  
   string line;
   ifstream infile(inputFileList.c_str());
   infile >> line; // the first line determines the type of input files
- 
+  
   DataConverter *converter=0;
-   
+  
   if(strstr(line.c_str(),".hep"))
     {                           
@@ -129,5 +130,5 @@
     }
   else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};
-
+  
   TChain chain("GEN");
   chain.Add(outputfilename.c_str());
@@ -143,8 +144,8 @@
   ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
   ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
-
+  
   TRootGenParticle *particle;
   TRootETmis *etmisc;
-
+  
   RESOLELEC *elementElec;
   RESOLMUON *elementMuon;
@@ -152,6 +153,6 @@
   TAUHAD *elementTaujet;
   ETMIS *elementEtmis;
-
- 
+  
+  
   //read the datacard input file
   string DetDatacard("");
@@ -163,36 +164,36 @@
   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<PhysicsTower> towers;
-
+  
   fastjet::JetDefinition jet_def;
   fastjet::JetDefinition jet_defS;
   fastjet::JetDefinition::Plugin * plugins;
   fastjet::JetDefinition::Plugin * pluginsS;
-
+  
   // set up a CDF midpoint jet definition
-  #ifdef ENABLE_PLUGIN_CDFCONES
+#ifdef ENABLE_PLUGIN_CDFCONES
   plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
   jet_def = fastjet::JetDefinition(plugins);
-  #else
+#else
   plugins = NULL;
-  #endif
-
+#endif
+  
   // set up a CDF midpoint jet definition
-  #ifdef ENABLE_PLUGIN_CDFCONES
-  pluginsS = new fastjet::CDFJetCluPlugin(2,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
+#ifdef ENABLE_PLUGIN_CDFCONES
+  pluginsS = new fastjet::CDFJetCluPlugin(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
   jet_defS = fastjet::JetDefinition(pluginsS);
-  #else
+#else
   pluginsS = NULL;
-  #endif
-
- 
+#endif
+  
+  
   // Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
@@ -218,107 +219,88 @@
       sorted_jetsS.clear();
       towers.clear();
-
+      
       // Loop over all particles in event
       while( (particle = (TRootGenParticle*) itGen.Next()) )
 	{
           genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
-
+	  
 	  int pid = abs(particle->PID);
 	  float eta = fabs(particle->Eta);
-
+	  
           if(particle->Status == 1)
-             {
-               input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
-             }
-
+	    {
+	      input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+	    }
+	  
 	  // 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 != pMU)PTmis = PTmis + genMomentum;//ptmis
-		switch(pid) {
-		    
-		    case pE: // all electrons with eta < DET->MAX_CALO_FWD
-		      DET->SmearElectron(genMomentum);
-		      break; // case pE
-
-		    case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-		      DET->SmearElectron(genMomentum);
-		      break; // case pGAMMA
-
-		    case pMU: // all muons with eta < DET->MAX_MU
-		      DET->SmearMu(genMomentum);
-		      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);
-		      break; // case hadron
-		      
-		    default:   // all other final particles with eta < DET->MAX_CALO_FWD
-		      DET->SmearHadron(genMomentum, 1.0);
-		      break;
-		} // switch (pid)
-
-		// all final particles but muons and neutrinos	
-		// for calorimetric towers and mission PT
-                //if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis
-
-		if(pid != pMU)
-                  {
-                     PTmisS = PTmisS + genMomentum;
-                     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()));
-                  }
-		
-                // 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
+	       (pid == pMU  && eta < DET->MAX_MU) ||
+	       (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) )
+	      ) 
+	    {
+	      PTmis = PTmis + genMomentum;//ptmis
+	      switch(pid) {
+
+	      case pE: // all electrons with eta < DET->MAX_CALO_FWD
+		DET->SmearElectron(genMomentum);
+		break; // case pE
+	      case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
+		DET->SmearElectron(genMomentum);
+		break; // case pGAMMA
+	      case pMU: // all muons with eta < DET->MAX_MU
+		DET->SmearMu(genMomentum);
+		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);
+		break; // case hadron
+	      default:   // all other final particles with eta < DET->MAX_CALO_FWD
+		DET->SmearHadron(genMomentum, 1.0);
+		break;
+	      } // switch (pid)
+	      
+	      if(pid != pMU)
+		{
+//		  PTmisS = PTmisS + genMomentum;
+		  towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
+		  if(genMomentum.Et()>0.5)input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.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);
+		TrackCentral.push_back(genMomentum);
 	    } // switch
 	} // while
-
-  TLorentzVector toWerS(0,0,0,0);
-  //calcul de ETmis au depart des tours calorimetriques.. 
-/*  for(unsigned int i=0; i < towers.size(); i++) {
-      if(towers[i].fourVector.pt() < 0.5) continue;
-          toWerS.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
-          PTmisS = PTmisS + toWerS;
-          }
-*/
-//PTmis = (0,0,0,0)-PTmis;
-//PTmisS = (0,0,0,0)-PTmisS;
-
-//float ET=PTmis.Et()<<endl;
-//float ETS=PTmisS.Et()<<endl;
-
-cout<<"Ptmis "<<PTmis.Et()<<" PTmis smeare "<<PTmisS.Et()<<endl;
-cout<<"Ptmis "<<(-PTmis).Pt()<<" PTmis smeare "<<(-PTmisS).Pt()<<endl;
-
-  elementEtmis= (ETMIS*) branchetmis->NewEntry();
-  elementEtmis->Et = (PTmis).Et();
-  elementEtmis->EtSmeare = (PTmisS).Et();
-
-
-
+     
+     TLorentzVector TowerTLV(0.,0.,0.,0.);
+     for(unsigned int i=0; i < towers.size(); i++) 
+       {
+         TowerTLV.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
+         if(TowerTLV.Et()>0.5)PTmisS = PTmisS + TowerTLV;
+      }
+
+ 
+      elementEtmis= (ETMIS*) branchetmis->NewEntry();
+      elementEtmis->Et = (PTmis).Et();
+      elementEtmis->EtSmeare = (PTmisS).Et()-(PTmis).Et();
+      
       //*****************************
-
-     double ptmin=1;
-     if(input_particles.size()!=0)
+      
+      double ptmin=1;
+      if(input_particles.size()!=0)
         {
           fastjet::ClusterSequence clust_seq(input_particles, jet_def);
@@ -326,6 +308,6 @@
           sorted_jets = sorted_by_pt(inclusive_jets);
         }
-
-     if(input_particlesS.size()!=0)
+      
+      if(input_particlesS.size()!=0)
         {
           fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
@@ -333,7 +315,7 @@
           sorted_jetsS = sorted_by_pt(inclusive_jetsS);
         }
-
-    TLorentzVector JETSm(0,0,0,0);
-    for (unsigned int i = 0; i < sorted_jets.size(); i++) {
+      
+      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());
@@ -341,15 +323,12 @@
         if(JETSm.Pt()>1)
           {
-            float NonSmeareEt=(JET.E()*sin(JET.Theta()));
-            float SmeareEt=(JETSm.E()*sin(JETSm.Theta()));
-
             elementJet= (RESOLJET*) branchjet->NewEntry();
             elementJet->NonSmearePT = JET.Et();
             elementJet->SmearePT = JETSm.Et()/JET.Et();
-          
+	    
           }
-
-	} // for itJet : loop on all jets
-
+	
+      } // for itJet : loop on all jets
+      
       treeWriter->Fill();
     } // Loop over all events
