Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 29)
+++ trunk/Delphes.cpp	(revision 30)
@@ -123,5 +123,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,6 +143,6 @@
   ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
   ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
-
- 
+  
+  
   TRootGenParticle *particle; 
   TRootETmis *elementEtmis;
@@ -163,13 +163,18 @@
   DET->ReadDataCard(DetDatacard);
   
-  
   TLorentzVector genMomentum(0,0,0,0);
   LorentzVector jetMomentum;
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
-
+  
   vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
   vector<fastjet::PseudoJet> inclusive_jets;
   vector<fastjet::PseudoJet> sorted_jets;
+  
+  vector<TLorentzVector> electron;
+  vector<int> elecPID;
+  vector<TLorentzVector> muon;
+  vector<int> muonPID;
+  TSimpleArray<TRootGenParticle> NFCentralQ;
   
   //Initialisation of Hector
@@ -196,57 +201,57 @@
   beamline2->add(rp220_2);
   beamline2->add(rp420_2);
-
+  
   // we will have four jet definitions, and the first three will be
   // plugins
   fastjet::JetDefinition jet_def;
   fastjet::JetDefinition::Plugin * plugins;
-
+  
   switch(DET->JETALGO) {
   default:
   case 1: {
-
-   // set up a CDF midpoint jet definition
-   #ifdef ENABLE_PLUGIN_CDFCONES
-   plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
-   jet_def = fastjet::JetDefinition(plugins);
-   #else
-   plugins = NULL;
-   #endif
-   }
-   break;
-
-   case 2: {
-
-   // set up a CDF midpoint jet definition
-   #ifdef ENABLE_PLUGIN_CDFCONES
-   plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
-   jet_def = fastjet::JetDefinition(plugins);
-   #else
-   plugins = NULL;
-   #endif
-   }
-   break;
-   case 3: {
-   // set up a siscone jet definition
-   #ifdef ENABLE_PLUGIN_SISCONE
-   int npass = 0;               // do infinite number of passes
-   double protojet_ptmin = 0.0; // use all protojets
-   plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
-   jet_def = fastjet::JetDefinition(plugins);
-   #else
-   plugins = NULL;
-   #endif
-   }
-   break;
-
+    
+    // set up a CDF midpoint jet definition
+#ifdef ENABLE_PLUGIN_CDFCONES
+    plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+    
+  case 2: {
+    
+    // set up a CDF midpoint jet definition
+#ifdef ENABLE_PLUGIN_CDFCONES
+    plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+  case 3: {
+    // set up a siscone jet definition
+#ifdef ENABLE_PLUGIN_SISCONE
+    int npass = 0;               // do infinite number of passes
+    double protojet_ptmin = 0.0; // use all protojets
+    plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
+    jet_def = fastjet::JetDefinition(plugins);
+#else
+    plugins = NULL;
+#endif
+  }
+    break;
+    
   case 4: {
-  jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
-  //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
-  //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
-  }
-  break;
-  }
- 
-   // Loop over all events
+    jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
+    //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
+    //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
+  }
+    break;
+  }
+  
+  // Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
   cout << "** Chain contains " << allEntries << " events" << endl;
@@ -259,22 +264,25 @@
       if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
       
-      TSimpleArray<TRootGenParticle> electron;
-      TSimpleArray<TRootGenParticle> muon;
+      electron.clear();
+      muon.clear();
+      elecPID.clear();
+      muonPID.clear();
+      NFCentralQ.Clear();
+      
       itGen.Reset();
       TrackCentral.clear();
       towers.clear();
       input_particles.clear();
-      TSimpleArray<TRootGenParticle> NFCentralQ;
       inclusive_jets.clear();
       sorted_jets.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);
-
+	  
           //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)
@@ -288,99 +296,78 @@
             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
           if( (particle->Status == 1)    &&
 	      ( 
-		(pid == pMU  && eta < DET->MAX_MU) ||
-                (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
-              )
-            ) { 
-		switch(pid) {
-		    
-		    case pE: // all electrons with eta < DET->MAX_CALO_FWD
-		      DET->SmearElectron(genMomentum);
-                      electron.Add(particle);
-		      break; // case pE
-		    case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-		      DET->SmearElectron(genMomentum);
-		      if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
-			  elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
-			  elementPhoton->Set(genMomentum);
-		      }
-		      break; // case pGAMMA
-		    case pMU: // all muons with eta < DET->MAX_MU
-		      DET->SmearMu(genMomentum);
-                      muon.Add(particle);
-		      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) {
-		    towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
-		    // create a fastjet::PseudoJet with these components and put it onto
-		    // back of the input_particles vector
-		    input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
-		    elementCalo = (TRootCalo*) branchCalo->NewEntry();
-		    elementCalo->Set(genMomentum);
-		  }
-		}
-		
-		
-		// 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
-		   )
-		  {
-		    elementTracks = (TRootTracks*) branchTracks->NewEntry();
-		    elementTracks->Set(genMomentum);
-		    TrackCentral.push_back(genMomentum);
-		  }
+	       (pid == pMU  && eta < DET->MAX_MU) ||
+	       (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) 
+		)
+	      ) { 
+	    switch(pid) {
+	      
+	    case pE: // all electrons with eta < DET->MAX_CALO_FWD
+              DET->SmearElectron(genMomentum);
+	      electron.push_back(genMomentum);
+              elecPID.push_back(particle->PID);
+	      break; // case pE
+	    case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
+	      DET->SmearElectron(genMomentum);
+	      if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
+		elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
+		elementPhoton->Set(genMomentum);
+	      }
+	      break; // case pGAMMA
+	    case pMU: // all muons with eta < DET->MAX_MU
+	      DET->SmearMu(genMomentum);
+              muonPID.push_back(particle->PID);
+	      muon.push_back(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) {
+		towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
+		// create a fastjet::PseudoJet with these components and put it onto
+		// back of the input_particles vector
+		input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
+		elementCalo = (TRootCalo*) branchCalo->NewEntry();
+		elementCalo->Set(genMomentum);
+	      }
+	    }
+	    
+	    
+	    // 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
+	       )
+	      {
+		elementTracks = (TRootTracks*) branchTracks->NewEntry();
+		elementTracks->Set(genMomentum);
+		TrackCentral.push_back(genMomentum);
+	      }
 	  } // switch
-
-          //branchElectron->Reset();
-          
-          if(electron.GetEntries()>0) {
-          for(int i=0; i < electron.GetEntries();i++) { 
-               genMomentum.SetPxPyPzE(electron[i]->Px,electron[i]->Py,electron[i]->Pz,electron[i]->E);
-               if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_TRACKER) {
-               elementElec = (TRootElectron*) branchElectron->NewEntry();
-               elementElec->Set(genMomentum);
-               elementElec->Charge = sign(electron[i]->PID);
-               if(fabs(genMomentum.Eta()))elementElec->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
-               }
-          }}
-          if(muon.GetEntries()>0) {
-          for(int i=0; i < muon.GetEntries();i++) {
-               genMomentum.SetPxPyPzE(muon[i]->Px,muon[i]->Py,muon[i]->Pz,muon[i]->E);
-               if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_MU) {
-                 elementMu = (TRootMuon*) branchMuon->NewEntry();
-                 elementMu->Charge = sign(particle->PID);
-                 elementMu->Set(genMomentum);
-                 elementMu->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
-               }
-          }}
-
+	  
 	  // Forward particles in CASTOR ?
 	  /*	  if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD) 
@@ -461,4 +448,24 @@
 	} // while
       
+      for(unsigned int i=0; i < electron.size(); i++) {
+	if(electron[i].E()!=0 && fabs(electron[i].Eta()) < DET->MAX_TRACKER)
+	  {
+	    elementElec = (TRootElectron*) branchElectron->NewEntry();
+	    elementElec->Set(electron[i]);
+	    elementElec->Charge = sign(elecPID[i]);
+	    elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,2.0);
+	  }
+      }
+      for(unsigned int i=0; i < muon.size(); i++) {
+	if(muon[i].E()!=0 && fabs(muon[i].Eta()) < DET->MAX_MU) 
+	  {
+	    elementMu = (TRootMuon*) branchMuon->NewEntry();
+	    elementMu->Charge = sign(muonPID[i]);
+	    elementMu->Set(muon[i]);
+	    elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,2.0);
+	  }
+      }
+      
+      
       // computes the Missing Transverse Momentum
       elementEtmis = (TRootETmis*) branchETmis->NewEntry();
@@ -469,5 +476,5 @@
       
       //*****************************
-
+      
       // run the jet clustering with the above jet definition
       if(input_particles.size()!=0)
