Index: trunk/Delphes.cpp
===================================================================
--- trunk/Delphes.cpp	(revision 306)
+++ trunk/Delphes.cpp	(revision 307)
@@ -47,4 +47,5 @@
 #include "LHEFConverter.h"
 #include "STDHEPConverter.h"
+#include "LHCOConverter.h"
 
 #include "SmearUtil.h"
@@ -77,5 +78,5 @@
 int main(int argc, char *argv[])
 {
-
+  
   int appargc = 2;
   char *appName= new char[20];
@@ -87,5 +88,5 @@
   delete [] appName;
   delete [] appOpt;
-
+  
   if(argc != 3 && argc != 4 && argc != 5) {
     cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl;
@@ -96,68 +97,68 @@
     exit(1);
   }
-
-
-  
- 
-cout << endl << endl;
-
-cout <<"*********************************************************************"<< endl;
-cout <<"*********************************************************************"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**                            Welcome to                           **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**     .ddddddd-            lL             hH                      **"<< endl;
-cout <<"**     -Dd` `dD:           Ll             hH`                      **"<< endl;
-cout <<"**     dDd   dDd    eeee.  lL  .pp+pp    Hh+hhh`   -eeee- `sssss   **"<< endl;
-cout <<"**    -Dd   `DD   ee. ee  Ll   .Pp. PP   Hh. HH.  ee. ee  sSs      **"<< endl;
-cout <<"**    dD`   dDd  eEeee:   lL.  pP.  pP  hH   hH` eEeee:` -sSSSs.   **"<< endl;
-cout <<"**   .Dd  :dd    eE.     LlL   PpppPP   Hh  Hh   eE         sSS    **"<< endl;
-cout <<"**   dddddd:.     eee+:  lL.  pp.      hh.  hh    eee+  sssssS     **"<< endl;
-cout <<"**                            Pp                                   **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**           Delphes, a framework for the fast simulation          **"<< endl;
-cout <<"**                 of a  generic collider experiment               **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**                --- Version 1.4beta of Delphes ---               **"<< endl;
-cout <<"**               Last date of change: 9 February 2009              **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**     This package uses:                                          **"<< endl;
-cout <<"**     ------------------                                          **"<< endl;
-cout <<"**     FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
-cout <<"**     Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]  **"<< endl;
-cout <<"**     FROG:                               [hep-ex/0901.2718v1]    **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**-----------------------------------------------------------------**"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**   Main authors:                                                 **"<< endl;
-cout <<"**   -------------                                                 **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**              SÃ©verine Ovyn               Xavier Rouby           **"<< endl;
-cout <<"**       severine.ovyn@uclouvain.be      xavier.rouby@cern         **"<< endl;
-cout <<"**       Center for Particle Physics and Phenomenology (CP3)       **"<< endl;
-cout <<"**       Universite Catholique de Louvain (UCL)                    **"<< endl;
-cout <<"**       Louvain-la-Neuve, Belgium                                 **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**-----------------------------------------------------------------**"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**   Former Delphes versions and documentation can be found on :   **"<< endl;
-cout <<"**               http://www.fynu.ucl.ac.be/delphes.html            **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"**   Disclaimer: this program is a beta version of Delphes and     **"<< endl;
-cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
-cout <<"**        give us your feedbacks about potential bugs              **"<< endl;
-cout <<"**                                                                 **"<< endl;
-cout <<"*********************************************************************"<< endl;
-cout <<"*********************************************************************"<< endl;
- 
-// 1. ********** initialisation ***********
-
+  
+  
+  
+  
+  cout << endl << endl;
+  
+  cout <<"*********************************************************************"<< endl;
+  cout <<"*********************************************************************"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                            Welcome to                           **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**     .ddddddd-            lL             hH                      **"<< endl;
+  cout <<"**     -Dd` `dD:           Ll             hH`                      **"<< endl;
+  cout <<"**     dDd   dDd    eeee.  lL  .pp+pp    Hh+hhh`   -eeee- `sssss   **"<< endl;
+  cout <<"**    -Dd   `DD   ee. ee  Ll   .Pp. PP   Hh. HH.  ee. ee  sSs      **"<< endl;
+  cout <<"**    dD`   dDd  eEeee:   lL.  pP.  pP  hH   hH` eEeee:` -sSSSs.   **"<< endl;
+  cout <<"**   .Dd  :dd    eE.     LlL   PpppPP   Hh  Hh   eE         sSS    **"<< endl;
+  cout <<"**   dddddd:.     eee+:  lL.  pp.      hh.  hh    eee+  sssssS     **"<< endl;
+  cout <<"**                            Pp                                   **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**           Delphes, a framework for the fast simulation          **"<< endl;
+  cout <<"**                 of a  generic collider experiment               **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                --- Version 1.4beta of Delphes ---               **"<< endl;
+  cout <<"**               Last date of change: 9 February 2009              **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**     This package uses:                                          **"<< endl;
+  cout <<"**     ------------------                                          **"<< endl;
+  cout <<"**     FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **"<< endl;
+  cout <<"**     Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]  **"<< endl;
+  cout <<"**     FROG:                               [hep-ex/0901.2718v1]    **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**-----------------------------------------------------------------**"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**   Main authors:                                                 **"<< endl;
+  cout <<"**   -------------                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**              SÃ©verine Ovyn               Xavier Rouby           **"<< endl;
+  cout <<"**       severine.ovyn@uclouvain.be      xavier.rouby@cern         **"<< endl;
+  cout <<"**       Center for Particle Physics and Phenomenology (CP3)       **"<< endl;
+  cout <<"**       Universite Catholique de Louvain (UCL)                    **"<< endl;
+  cout <<"**       Louvain-la-Neuve, Belgium                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**-----------------------------------------------------------------**"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**   Former Delphes versions and documentation can be found on :   **"<< endl;
+  cout <<"**               http://www.fynu.ucl.ac.be/delphes.html            **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**   Disclaimer: this program is a beta version of Delphes and     **"<< endl;
+  cout <<"** therefore comes without guarantees. Beware of errors and please **"<< endl;
+  cout <<"**        give us your feedbacks about potential bugs              **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"*********************************************************************"<< endl;
+  cout <<"*********************************************************************"<< endl;
+  
+  // 1. ********** initialisation ***********
+  
   srand (time (NULL));         /* Initialisation du gÃ©nÃ©rateur */
   TStopwatch globalwatch, loopwatch, triggerwatch, frogwatch;
   globalwatch.Start();
-
+  
   
   //read the output TROOT file
@@ -182,5 +183,5 @@
   string DetDatacard("data/DetectorCard.dat"); //for detector smearing parameters
   string TrigDatacard("data/TriggerCard.dat"); //for trigger selection
- 
+  
   string lineCard1,lineCard2;
   bool detecCard=false,trigCard=false;
@@ -190,8 +191,8 @@
       infile1 >> lineCard1; // the first line determines the type of input files
       if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==true) 
-           cerr <<"**         ERROR: A DETECTOR card has already been loaded          **"<< endl;
+	cerr <<"**         ERROR: A DETECTOR card has already been loaded          **"<< endl;
       else if(strstr(lineCard1.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[3]; detecCard=true;}
       else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==true)
-           cerr <<"**         ERROR: A TRIGGER card has already been loaded           **"<< endl;
+	cerr <<"**         ERROR: A TRIGGER card has already been loaded           **"<< endl;
       else if(strstr(lineCard1.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[3]; trigCard=true;}
     }
@@ -201,8 +202,8 @@
       infile2 >> lineCard2; // the first line determines the type of input files
       if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==true)
-           cerr <<"**         ERROR: A DETECTOR card has already been loaded          **"<< endl;
+	cerr <<"**         ERROR: A DETECTOR card has already been loaded          **"<< endl;
       else if(strstr(lineCard2.c_str(),"DETECTOR") && detecCard==false){DetDatacard =argv[4]; detecCard=true;}
       else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==true)
-           cerr <<"**         ERROR: A TRIGGER card has already been loaded           **"<< endl;
+	cerr <<"**         ERROR: A TRIGGER card has already been loaded           **"<< endl;
       else if(strstr(lineCard2.c_str(),"TRIGGER") && trigCard==false){TrigDatacard =argv[4]; trigCard=true;}
     }
@@ -221,12 +222,12 @@
        << right << setw(2) <<"**"<<""<<endl;
   cout <<"**                                                                 **"<< endl;
-
+  
   //Trigger information
   cout <<"**        ########### Start reading TRIGGER card ##########        **"<< endl;
   if(trigCard==false)
-     {
-       cout <<"**        WARNING: Datacard  not found, use default card           **" << endl;
-       TrigDatacard="data/TriggerCard.dat";
-     }
+    {
+      cout <<"**        WARNING: Datacard  not found, use default card           **" << endl;
+      TrigDatacard="data/TriggerCard.dat";
+    }
   TriggerTable *TRIGT = new TriggerTable();
   TRIGT->TriggerCardReader(TrigDatacard.c_str());
@@ -248,10 +249,10 @@
   //VFD information
   VeryForward * VFD = new VeryForward(DET);
- 
+  
   // data converters
   DataConverter *converter=NULL;
   cout <<"**                                                                 **"<<endl;
   cout <<"**        ####### Start convertion to TRoot format ########        **"<< endl;
-
+  
   if(strstr(line.c_str(),".hep"))
     {                           
@@ -273,18 +274,18 @@
     }
   else {
-      cerr << left  << setw(4) <<"**  "<<""
-           << left  << setw(63) << line.c_str() <<""
-           << right << setw(2) <<"**"<<endl;
-      cerr <<"**         ERROR: File format not identified --  Exiting...        **"<< endl;
-      cout <<"**                                                                 **"<< endl;
-      cout <<"*********************************************************************"<< endl;
-      return -1;};
+    cerr << left  << setw(4) <<"**  "<<""
+	 << left  << setw(63) << line.c_str() <<""
+	 << right << setw(2) <<"**"<<endl;
+    cerr <<"**         ERROR: File format not identified --  Exiting...        **"<< endl;
+    cout <<"**                                                                 **"<< endl;
+    cout <<"*********************************************************************"<< endl;
+    return -1;};
   cout <<"**                       Exiting conversion...                     **"<< endl;
-
+  
   TChain chain("GEN");
   chain.Add(outputfilename.c_str());
   ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
   const TClonesArray *branchGen = treeReader->UseBranch("Particle");
-
+  
   TIter itGen((TCollection*)branchGen);
   
@@ -322,31 +323,33 @@
   vector<D_Particle> muon;
   vector<D_Particle> gamma;
-
-bool FLAG_lhco = true;
-
+ 
+  vector<int> NTrackJet;
+
+  bool FLAG_lhco = true;
+  
   TSimpleArray<TRootGenParticle> NFCentralQ;
-
-D_CaloList list_of_calorimeters;
-D_CaloElement CentralCalo("centralcalo", 
-	-DET->CEN_max_calo_cen, DET->CEN_max_calo_cen,
-	 DET->ELG_Ccen,  DET->ELG_Ncen,  DET->ELG_Scen,
-         DET->HAD_Chcal, DET->HAD_Nhcal, DET->HAD_Shcal);
-D_CaloElement ForwardCalo("forwardcalo",
-	 DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd, 
-         DET->ELG_Cfwd,  DET->ELG_Nfwd,  DET->ELG_Sfwd,
-         DET->HAD_Chf,   DET->HAD_Nhf,   DET->HAD_Shf  );
-D_CaloElement BackwardCalo("backwardcalo",
-	-DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen,
-         DET->ELG_Cfwd,  DET->ELG_Nfwd,  DET->ELG_Sfwd,
-         DET->HAD_Chf,   DET->HAD_Nhf,   DET->HAD_Shf  );
-//D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0);
-list_of_calorimeters.addElement(CentralCalo);
-list_of_calorimeters.addElement(ForwardCalo);
-list_of_calorimeters.addElement(BackwardCalo);
-//list_of_calorimeters.addElement(CastorCalo);
-list_of_calorimeters.sortElements();
-
-
-// 2. ********** Loop over all events ***********
+  
+  D_CaloList list_of_calorimeters;
+  D_CaloElement CentralCalo("centralcalo", 
+			    -DET->CEN_max_calo_cen, DET->CEN_max_calo_cen,
+			    DET->ELG_Ccen,  DET->ELG_Ncen,  DET->ELG_Scen,
+			    DET->HAD_Chcal, DET->HAD_Nhcal, DET->HAD_Shcal);
+  D_CaloElement ForwardCalo("forwardcalo",
+			    DET->CEN_max_calo_cen, DET->CEN_max_calo_fwd, 
+			    DET->ELG_Cfwd,  DET->ELG_Nfwd,  DET->ELG_Sfwd,
+			    DET->HAD_Chf,   DET->HAD_Nhf,   DET->HAD_Shf  );
+  D_CaloElement BackwardCalo("backwardcalo",
+			     -DET->CEN_max_calo_fwd, -DET->CEN_max_calo_cen,
+			     DET->ELG_Cfwd,  DET->ELG_Nfwd,  DET->ELG_Sfwd,
+			     DET->HAD_Chf,   DET->HAD_Nhf,   DET->HAD_Shf  );
+  //D_CaloElement CastorCalo("castor",5.5,6.6,1,0,0,1,0,0);
+  list_of_calorimeters.addElement(CentralCalo);
+  list_of_calorimeters.addElement(ForwardCalo);
+  list_of_calorimeters.addElement(BackwardCalo);
+  //list_of_calorimeters.addElement(CastorCalo);
+  list_of_calorimeters.sortElements();
+  
+  
+  // 2. ********** Loop over all events ***********
   Long64_t entry, allEntries = treeReader->GetEntries();
   cout <<"**                                                                 **"<<endl;
@@ -355,9 +358,9 @@
        << left  << setw(15) << allEntries <<""
        << right << setw(2) <<"**"<<endl;
-
+  
   ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
-
+  
   loopwatch.Start();
-
+  
   // loop on all events
   for(entry = 0; entry < allEntries; ++entry)
@@ -376,20 +379,21 @@
       towers.clear();
       input_particles.clear();
-
+      NTrackJet.clear();     
+ 
       {
-      D_CaloTowerList list_of_active_towers; 
-      D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates
-
-      // 2.1a Loop over all particles in event, to fill the towers
-      itGen.Reset();
-      GenParticle *particleG; 
-      while( (particleG = (GenParticle*) itGen.Next()) ) {
-
-      TRootGenParticle *particle = new TRootGenParticle(particleG);
-      int pid = abs(particle->PID);
-      particle->Charge=ChargeVal(particle->PID);
-      particle->setFractions(); // init
-
-
+	D_CaloTowerList list_of_active_towers; 
+	D_CaloTowerList list_of_towers_with_photon; // to speed up the code: will only look in interesting towers for gamma candidates
+	
+	// 2.1a Loop over all particles in event, to fill the towers
+	itGen.Reset();
+	GenParticle *particleG; 
+	while( (particleG = (GenParticle*) itGen.Next()) ) {
+	  
+	  TRootGenParticle *particle = new TRootGenParticle(particleG);
+	  int pid = abs(particle->PID);
+	  particle->Charge=ChargeVal(particle->PID);
+	  particle->setFractions(); // init
+	  
+	  
 	  // 2.1a.1********************* preparation for the b-tagging
 	  //// This subarray is needed for the B-jet algorithm
@@ -401,8 +405,8 @@
             NFCentralQ.Add(particle);
           }
-
+	  
 	  // 2.1a.2********************* visible particles only
 	  if( (particle->Status == 1) && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) ){
-
+	    
             // 2.1a.2.1 Central solenoidal magnetic field 
 	    TRACP->bfield(particle); // fills in particle->EtaCalo et particle->PhiCalo
@@ -412,22 +416,22 @@
                 particle->Charge==0 ||
                 (!DET->FLAG_bfield && particle->Charge!=0 && particle->PT > DET->TRACK_ptmin)) 
-	    if(
-		(particle->EtaCalo > list_of_calorimeters.getEtamin() ) &&
-		(particle->EtaCalo < list_of_calorimeters.getEtamax() )
-              ) {
+	      if(
+		 (particle->EtaCalo > list_of_calorimeters.getEtamin() ) &&
+		 (particle->EtaCalo < list_of_calorimeters.getEtamax() )
+		 ) {
 		float iEta=UNDEFINED, iPhi=UNDEFINED;
 		DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
 		if (iEta != UNDEFINED && iPhi != UNDEFINED) 
-                {
-		  D_CaloTower tower(iEta,iPhi);	                               // new tower
-		  tower.Set_Eem_Ehad_E_ET(particle->E*particle->getFem() , particle->E*particle->getFhad() );
-		
-		  list_of_active_towers.addTower(tower);
-		  // this list may contain several times the same calotower, as several particles 
-		  // may leave some energy in the same calotower
-		  // After the loop on particles, identical cells in the list should be merged 
-		} // iEta and iPhi must be defined
-	    }
-
+		  {
+		    D_CaloTower tower(iEta,iPhi);	                               // new tower
+		    tower.Set_Eem_Ehad_E_ET(particle->E*particle->getFem() , particle->E*particle->getFhad() );
+		    
+		    list_of_active_towers.addTower(tower);
+		    // this list may contain several times the same calotower, as several particles 
+		    // may leave some energy in the same calotower
+		    // After the loop on particles, identical cells in the list should be merged 
+		  } // iEta and iPhi must be defined
+	      }
+	    
  	    // 2.1a.2.3 charged particles in tracker: energy flow
 	    // if bfield not simulated, pt should be high enough to be taken into account
@@ -435,71 +439,71 @@
 	    if( particle->Charge !=0 &&
 	        fabs(particle->EtaCalo)< DET->CEN_max_tracker && // stays in the tracker -> track available
-	       ( DET->FLAG_bfield || 
+		( DET->FLAG_bfield || 
                   (!DET->FLAG_bfield && particle->PT > DET->TRACK_ptmin)
-                )
-              ) {	
-		// 2.1a.2.3.1 Filling the particle properties + smearing
-		   // Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction
-		   // This is the EnergyFlow hypothesis
-		   particle->SetEtaPhi(particle->Eta,particle->Phi);
-		   float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks
-
-		   // 2.1a.2.3.2 Muons
-		   if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu) {
-                        TLorentzVector p;
-                        float sPT = gRandom->Gaus(particle->PT, DET->MU_SmearPt*particle->PT );
-			if (sPT > 0 && sPT > DET->PTCUT_muon) {
-                           p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta));
-                           muon.push_back(D_Particle(p,pMU,particle->EtaCalo,particle->PhiCalo));
-                        }
-			sET = (sPT >0)? sPT : 0;
-                   }  
-		   // 2.1a.2.3.3 Electrons
-		   else if (pid == pE) {
-		     // Finds in which calorimeter the particle has gone, to know its resolution
-
-		     D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
-		     if(currentCalo.getName() == dummyCalo.getName()) {
-			cout << "** Warning: the calo coverage behind the tracker is not complete! **" << endl; }
-
-		     // final smeared EM energy   // electromagnetic fraction F_em =1 for electrons;
-	             float sE  = currentCalo.getElectromagneticResolution().Smear(particle->E);
-		     if (sE>0) {
- 		       sET = sE/cosh(particle->Eta); 
-		       // NB: ET is found via the calorimetry and not via the track curvature
-
-		       TLorentzVector p;
-		       p.SetPtEtaPhiE(sET,particle->Eta,particle->Phi,sE);
-		       if (sET > DET->PTCUT_elec)
-		         electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
-                     } else { sET=0;} // if negative smeared energy -- needed for the tracks
-		   }  
-		   // 2.1a.2.3.4 Other charged particles : smear them for the tracks!
-		   else {  //other particles
-			D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
-			float sEem  = currentCalo.getElectromagneticResolution().Smear(particle->E * particle->getFem());
-			float sEhad = currentCalo.getHadronicResolution().Smear(particle->E * particle->getFhad());
-			float sE = ( (sEem>0)? sEem : 0 ) + ( (sEhad>0)? sEhad : 0 );
-			sET = sE/cosh(particle->EtaCalo);
-		   }
-
-		   // 2.1a.2.3.5 Tracks
-		   if( (rand()%100) < DET->TRACK_eff && sET!=0) {
-   		     elementTrack = (TRootTracks*) branchTrack->NewEntry();
-		     elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge);
-		     //TrackCentral.push_back(elementTrack->GetFourVector()); // tracks at vertex!
-		     TrackCentral.push_back(*elementTrack); // tracks at vertex!
-		     // TODO!!! apply a smearing on the position of the origin of the track
-		     // TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout);
-		   }
+		  )
+		) {	
+	      // 2.1a.2.3.1 Filling the particle properties + smearing
+	      // Hypothesis: the final eta/phi are the ones from the generator, thanks to the track reconstruction
+	      // This is the EnergyFlow hypothesis
+	      particle->SetEtaPhi(particle->Eta,particle->Phi);
+	      float sET=UNDEFINED; // smeared ET, computed from the smeared E -> needed for the tracks
+	      
+	      // 2.1a.2.3.2 Muons
+	      if (pid == pMU && fabs(particle->EtaCalo)< DET->CEN_max_mu) {
+		TLorentzVector p;
+		float sPT = gRandom->Gaus(particle->PT, DET->MU_SmearPt*particle->PT );
+		if (sPT > 0 && sPT > DET->PTCUT_muon) {
+		  p.SetPtEtaPhiE(sPT,particle->Eta,particle->Phi,sPT*cosh(particle->Eta));
+		  muon.push_back(D_Particle(p,pMU,particle->EtaCalo,particle->PhiCalo));
+		}
+		sET = (sPT >0)? sPT : 0;
+	      }  
+	      // 2.1a.2.3.3 Electrons
+	      else if (pid == pE) {
+		// Finds in which calorimeter the particle has gone, to know its resolution
+		
+		D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
+		if(currentCalo.getName() == dummyCalo.getName()) {
+		  cout << "** Warning: the calo coverage behind the tracker is not complete! **" << endl; }
+		
+		// final smeared EM energy   // electromagnetic fraction F_em =1 for electrons;
+		float sE  = currentCalo.getElectromagneticResolution().Smear(particle->E);
+		if (sE>0) {
+		  sET = sE/cosh(particle->Eta); 
+		  // NB: ET is found via the calorimetry and not via the track curvature
+		  
+		  TLorentzVector p;
+		  p.SetPtEtaPhiE(sET,particle->Eta,particle->Phi,sE);
+		  if (sET > DET->PTCUT_elec)
+		    electron.push_back(D_Particle(p,particle->PID,particle->EtaCalo,particle->PhiCalo));
+		} else { sET=0;} // if negative smeared energy -- needed for the tracks
+	      }  
+	      // 2.1a.2.3.4 Other charged particles : smear them for the tracks!
+	      else {  //other particles
+		D_CaloElement currentCalo = list_of_calorimeters.getElement(particle->EtaCalo);
+		float sEem  = currentCalo.getElectromagneticResolution().Smear(particle->E * particle->getFem());
+		float sEhad = currentCalo.getHadronicResolution().Smear(particle->E * particle->getFhad());
+		float sE = ( (sEem>0)? sEem : 0 ) + ( (sEhad>0)? sEhad : 0 );
+		sET = sE/cosh(particle->EtaCalo);
+	      }
+	      
+	      // 2.1a.2.3.5 Tracks
+	      if( (rand()%100) < DET->TRACK_eff && sET!=0) {
+		elementTrack = (TRootTracks*) branchTrack->NewEntry();
+		elementTrack->Set(particle->Eta, particle->Phi, particle->EtaCalo, particle->PhiCalo, sET, particle->Charge);
+		//TrackCentral.push_back(elementTrack->GetFourVector()); // tracks at vertex!
+		TrackCentral.push_back(*elementTrack); // tracks at vertex!
+		// TODO!!! apply a smearing on the position of the origin of the track
+		// TODO!!! elementTracks->SetPositionOut(Xout,Yout,Zout);
+	      }
             } // 2.1a.2.3 : if tracker/energy-flow
  	    // 2.1a.2.4 Photons
-               // stays in the tracker -> track available -> gamma ID
+	    // stays in the tracker -> track available -> gamma ID
 	    else if( (pid == pGAMMA) && fabs(particle->EtaCalo)< DET->CEN_max_tracker ) {
-		float iEta=UNDEFINED, iPhi=UNDEFINED;
-                DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
-                D_CaloTower tower(iEta,iPhi); 
-		// stores the list of towers where to apply the photon ID algorithm. Just a trick for a faster search
-		list_of_towers_with_photon.addTower(tower); 
+	      float iEta=UNDEFINED, iPhi=UNDEFINED;
+	      DET->BinEtaPhi(particle->PhiCalo,particle->EtaCalo,iPhi,iEta); // fills in iPhi and iEta
+	      D_CaloTower tower(iEta,iPhi); 
+	      // stores the list of towers where to apply the photon ID algorithm. Just a trick for a faster search
+	      list_of_towers_with_photon.addTower(tower); 
 	    }
 	    // 2.1a.2.5 : very forward detectors
@@ -507,49 +511,49 @@
               {
                 if (DET->FLAG_RP==1) {
-		   // for the moment, only protons are transported
-		   // BUT !!! could be a beam of other particles! (heavy ions?)
-		   // BUT ALSO !!! if very forward muons, or others!
-                   VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
+		  // for the moment, only protons are transported
+		  // BUT !!! could be a beam of other particles! (heavy ions?)
+		  // BUT ALSO !!! if very forward muons, or others!
+		  VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
 	        }
 	        // 2.1a.2.6: Zero degree calorimeter
 	        if(DET->FLAG_vfd==1) {
-                   VFD->ZDC(treeWriter,branchZDC,particle);
+		  VFD->ZDC(treeWriter,branchZDC,particle);
 	        }
               }
-
+	    
           } // 2.1a.2 : if visible particle
-        delete particle;
-      } // loop on all particles 2.1a
+	  delete particle;
+	} // loop on all particles 2.1a
 	
-
-      // 2.1b loop on all (activated) towers
-      // at this stage, list_of_active_towers may contain several times the same tower
-      // first step is to merge identical towers, by matching their (iEta,iPhi)  
-      list_of_active_towers.mergeDuplicates();
-      // Calotower smearing 
-      list_of_active_towers.smearTowers(list_of_calorimeters);
-
-      for(unsigned int i=0; i<list_of_active_towers.size(); i++) {
-	float iEta = list_of_active_towers[i].getEta();
-	float iPhi = list_of_active_towers[i].getPhi();
-	float e = list_of_active_towers[i].getE(); 
-	if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0) {
-           elementCalo = (TRootCalo*) branchCalo->NewEntry();
-	   elementCalo->set(list_of_active_towers[i]);
-           // not beautiful : should be improved!
-	   TLorentzVector p;
-	   p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e );
-           PhysicsTower Tower(LorentzVector(p.Px(),p.Py(),p.Pz(),p.E()));
-           towers.push_back(Tower);
-	}
-      } // loop on towers
-
-      // 2.1c photon ID
-      // list_of_towers_with_photon is the list of towers with photon candidates
-      // already smeared !
-      // sorts the vector and smears duplicates 
-      list_of_towers_with_photon.mergeDuplicates(); 
-
-      for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) {
+	
+	// 2.1b loop on all (activated) towers
+	// at this stage, list_of_active_towers may contain several times the same tower
+	// first step is to merge identical towers, by matching their (iEta,iPhi)  
+	list_of_active_towers.mergeDuplicates();
+	// Calotower smearing 
+	list_of_active_towers.smearTowers(list_of_calorimeters);
+	
+	for(unsigned int i=0; i<list_of_active_towers.size(); i++) {
+	  float iEta = list_of_active_towers[i].getEta();
+	  float iPhi = list_of_active_towers[i].getPhi();
+	  float e = list_of_active_towers[i].getE(); 
+	  if(iEta != UNDEFINED && iPhi != UNDEFINED && e!=0) {
+	    elementCalo = (TRootCalo*) branchCalo->NewEntry();
+	    elementCalo->set(list_of_active_towers[i]);
+	    // not beautiful : should be improved!
+	    TLorentzVector p;
+	    p.SetPtEtaPhiE(list_of_active_towers[i].getET(), iEta, iPhi, e );
+	    PhysicsTower Tower(LorentzVector(p.Px(),p.Py(),p.Pz(),p.E()));
+	    towers.push_back(Tower);
+	  }
+	} // loop on towers
+	
+	// 2.1c photon ID
+	// list_of_towers_with_photon is the list of towers with photon candidates
+	// already smeared !
+	// sorts the vector and smears duplicates 
+	list_of_towers_with_photon.mergeDuplicates(); 
+	
+	for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) {
 	  float eta = list_of_towers_with_photon[i].getEta(); 
 	  float phi = list_of_towers_with_photon[i].getPhi();
@@ -560,11 +564,11 @@
 	    if (cal.getET() > DET->PTCUT_gamma) { gamma.push_back(D_Particle(p,pGAMMA,p.Eta(),p.Phi())); }
           }
-      } // for -- list of photons
-
-}  // IF NEW ALGORITHM with energy flow
-
-
-
-
+	} // for -- list of photons
+	
+      }  // IF NEW ALGORITHM with energy flow
+      
+      
+      
+      
       // 2.2 ********** Output preparation & complex objects  ***********
       // 2.2.1 ********************* sorting collections by decreasing pt
@@ -599,10 +603,10 @@
 	  Att.SetPxPyPzE(towers[i].fourVector.px, towers[i].fourVector.py, towers[i].fourVector.pz, towers[i].fourVector.E);
 	  if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
-	  {
-            PTmis = PTmis + Att;
-            // create a fastjet::PseudoJet with these components and put it onto
-            // back of the input_particles vector
-            input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
-	  }
+	    {
+	      PTmis = PTmis + Att;
+	      // create a fastjet::PseudoJet with these components and put it onto
+	      // back of the input_particles vector
+	      input_particles.push_back(fastjet::PseudoJet(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E));
+	    }
         }
       elementEtmis = (TRootETmis*) branchETmis->NewEntry();
@@ -613,22 +617,23 @@
       
       // 2.2.3 ************* jets, B-tag, tau jets
-      sorted_jets=JETRUN->RunJets(input_particles);
-      JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ);
-      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral);
+      vector<int> NTrackJet;
+      sorted_jets=JETRUN->RunJets(input_particles, TrackCentral,NTrackJet);
+      JETRUN->RunJetBtagging(treeWriter, branchJet,sorted_jets,NFCentralQ,NTrackJet);
+      JETRUN->RunTauJets(treeWriter,branchTauJet,sorted_jets,towers, TrackCentral,NTrackJet);
       
       treeWriter->Fill();
     } // 2. Loop over all events ('for' loop)
-
-
+  
+  
   cout <<"**                                                                 **"<< endl;
   cout <<"**                 Exiting detector simulation...                  **"<< endl;
-
- 
+  
+  
   treeWriter->Write();
   delete treeWriter;
   loopwatch.Stop();
- 
-
- 
+  
+  
+  
   // 3. ********** Trigger  & Frog ***********
   // 3.1 ************ running the trigger in case the FLAG trigger is put to 1 in the datacard
@@ -636,7 +641,7 @@
   if(DET->FLAG_trigger == 1)
     {
-       cout <<"**                                                                 **"<<endl;
-       cout <<"**        ########### Start Trigger selection ###########          **"<< endl;
-
+      cout <<"**                                                                 **"<<endl;
+      cout <<"**        ########### Start Trigger selection ###########          **"<< endl;
+      
       // input
       TChain chainT("Analysis");
@@ -654,13 +659,13 @@
       ExRootTreeWriter *treeWriterT = new ExRootTreeWriter(outputfilename, "Trigger");
       ExRootTreeBranch *branchTrigger = treeWriterT->NewBranch("TrigResult", TRootTrigger::Class());
-     
- 
+      
+      
       Long64_t entryT, allEntriesT = treeReaderT->GetEntries();
       // loop on all entries
       for(entryT = 0; entryT < allEntriesT; ++entryT) {
-          treeWriterT->Clear();
-          treeReaderT->ReadEntry(entryT);
-          TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger);
-          treeWriterT->Fill();
+	treeWriterT->Clear();
+	treeReaderT->ReadEntry(entryT);
+	TRIGT->GetGlobalResult(branchElecTrig, branchMuonTrig,branchJetTrig, branchTauJetTrig,branchPhotonTrig, branchETmisTrig,branchTrigger);
+	treeWriterT->Fill();
       } // loop on all entries
       cout <<"**                 Exiting trigger simulation...                   **"<< endl;
@@ -670,31 +675,48 @@
       delete treeReaderT;
     } // trigger
-    triggerwatch.Stop();
- 
- 
+  triggerwatch.Stop();
+  
+  
   // 3.2 ************** FROG display
   frogwatch.Start();
   if(DET->FLAG_frog == 1) {
-      cout <<"**                                                                 **"<<endl;
-      cout <<"**        ################## Start FROG #################          **"<< endl;
-
-      FrogDisplay *FROG = new FrogDisplay(DET);
-      FROG->BuildEvents(outputfilename);
-      FROG->BuildGeom();
-      delete FROG;
+    cout <<"**                                                                 **"<<endl;
+    cout <<"**        ################## Start FROG #################          **"<< endl;
+    
+    FrogDisplay *FROG = new FrogDisplay(DET);
+    FROG->BuildEvents(outputfilename);
+    FROG->BuildGeom();
+    delete FROG;
   }
   frogwatch.Stop();
   
-
-
-
-// 4. ********** End & Exit ***********
-
+  // 3.3 *************** LHCO output
+  if(DET->FLAG_lhco == 1){
+    cout <<"**                                                                 **"<<endl;
+    cout <<"**        ############ Start LHCO convertion ############          **"<< endl;
+    /*
+    //create output log-file name
+    string forLogLHCO = outputfilename;
+    string LogNameLHCO = forLogLHCO.erase(forLogLHCO.find(".root"));
+    LogNameLHCO = LogNameLHCO+"_events.lhco";
+    */
+    
+    //LHCOConverter *LHCO = new LHCOConverter(outputfilename,LogNameLHCO);
+    LHCOConverter *LHCO = new LHCOConverter(outputfilename,"");
+    LHCO->CopyRunLogFile();
+    LHCO->ConvertExRootAnalysisToLHCO();
+    delete LHCO;
+  }
+  
+  
+  
+  // 4. ********** End & Exit ***********
+  
   globalwatch.Stop();
   cout <<"**                                                                 **"<< endl;
   cout <<"**        ################## Time report #################         **"<< endl;
   cout << left  << setw(32) <<"**              Time report for "<<""
-              << left  << setw(15) << allEntries <<""
-              << right << setw(22) <<"events         **"<<endl;
+       << left  << setw(15) << allEntries <<""
+       << right << setw(22) <<"events         **"<<endl;
   cout <<"**                                                                 **"<< endl;
   cout << left  << setw(10) <<"**"<<""
@@ -715,27 +737,27 @@
   if(DET->FLAG_trigger == 1)
     {
-     cout << left  << setw(10) <<"**"<<""
-          << left  << setw(15) <<" +  Trigger:"<<""
-          << right << setw(15) <<triggerwatch.CpuTime()<<""
-          << right << setw(15) <<triggerwatch.RealTime()<<""
-          << right << setw(14)  <<"**"<<endl;
+      cout << left  << setw(10) <<"**"<<""
+	   << left  << setw(15) <<" +  Trigger:"<<""
+	   << right << setw(15) <<triggerwatch.CpuTime()<<""
+	   << right << setw(15) <<triggerwatch.RealTime()<<""
+	   << right << setw(14)  <<"**"<<endl;
     }
   if(DET->FLAG_frog == 1)
     {
-     cout << left  << setw(10) <<"**"<<""
-          << left  << setw(15) <<" +  Frog:"<<""
-          << right << setw(15) <<frogwatch.CpuTime()<<""
-          << right << setw(15) <<frogwatch.RealTime()<<""
-          << right << setw(14)  <<"**"<<endl;
-
-    }
-  
-   cout <<"**                                                                 **"<< endl;
-   cout <<"**                        Exiting Delphes ...                      **"<< endl;
-   cout <<"**                                                                 **"<< endl;
-   cout <<"*********************************************************************"<< endl;
-   cout <<"*********************************************************************"<< endl;
-
-
+      cout << left  << setw(10) <<"**"<<""
+	   << left  << setw(15) <<" +  Frog:"<<""
+	   << right << setw(15) <<frogwatch.CpuTime()<<""
+	   << right << setw(15) <<frogwatch.RealTime()<<""
+	   << right << setw(14)  <<"**"<<endl;
+      
+    }
+  
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                        Exiting Delphes ...                      **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"*********************************************************************"<< endl;
+  cout <<"*********************************************************************"<< endl;
+  
+  
   delete treeReader;
   delete DET;
@@ -746,4 +768,4 @@
   delete converter;
   
-//  todo("TODO");
+  //  todo("TODO");
 }
