Merge lp:~ma5/madanalysis5/v1.1.7beta into lp:madanalysis5

[Benjamin Fuks]

- LHCO output
 - Compatibility with pythia-8 output
 - Improvement of the jet clustering / object identification
 - Few bug fixes
 - Photons added in the reco mode

 - LHCO output
 - Compatibility with pythia-8 output
 - Improvement of the jet clustering / object identification
 - Few bug fixes
 - Photons added in the reco mode

[Eric Conte (eric-conte)]

release validated !

review: Approve

=== modified file 'README'
--- README	2013-05-17 17:10:55 +0000
+++ README	2013-06-13 12:22:26 +0000
@@ -15,10 +15,10 @@
 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 #  GNU General Public License for more details.
-#  
+#
 #  You should have received a copy of the GNU General Public License
 #  along with MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
-#  
+#
 ################################################################################
 
 
 
=== modified file 'bin/ma5' (properties changed: -x to +x)
--- bin/ma5	2013-05-17 16:41:08 +0000
+++ bin/ma5	2013-06-13 12:22:26 +0000
@@ -113,8 +113,8 @@
 
 # Release version
 # Do not touch it !!!!!  
-main.version = "1.1.6"
-main.date = "2013/05/17"
+main.version = "1.1.25"
+main.date = "2013/06/13"
 
 # Checking arguments
 import getopt
 
=== modified file 'madanalysis/IOinterface/job_writer.py'
--- madanalysis/IOinterface/job_writer.py	2013-05-17 07:04:14 +0000
+++ madanalysis/IOinterface/job_writer.py	2013-06-13 12:22:26 +0000
@@ -253,7 +253,10 @@
         if self.output!="":
             file.write('  //Getting pointer to the writer\n')
             file.write('  WriterBase* writer1 = \n')
-            file.write('      manager.InitializeWriter("lhe","'+self.output+'");\n')
+            if self.output.lower().endswith('lhe') or self.output.lower().endswith('lhe.gz'):
+                file.write('      manager.InitializeWriter("lhe","'+self.output+'");\n')
+            elif self.output.lower().endswith('lhco') or self.output.lower().endswith('lhco.gz'):
+                file.write('      manager.InitializeWriter("lhco","'+self.output+'");\n')
             file.write('  if (writer1==0) return 1;\n\n')
         if self.clustering.algorithm!="none":
             file.write('  //Getting pointer to the clusterer\n')
@@ -958,11 +961,15 @@
         if not dataset.weighted_events:
             weighted_events=" --no_event_weight"
 
+        # Release
+        release = ' --ma5_version="' + self.main.version + ';' + self.main.date + '"'
+
         # Running SampleAnalyzer 
         res=os.system('cd '\
                       +self.path+'/Output/'+name+';'\
                       +' ../../Build/'\
                       +'MadAnalysis5Job '+weighted_events +\
+                      release +\
                       ' ../../Input/'+name+'.list')
 
         return True
 
=== modified file 'madanalysis/IOinterface/multiparticle_reader.py'
--- madanalysis/IOinterface/multiparticle_reader.py	2013-05-17 06:19:22 +0000
+++ madanalysis/IOinterface/multiparticle_reader.py	2013-06-13 12:22:26 +0000
@@ -29,12 +29,13 @@
 
 class MultiparticleReader():
 
-    def __init__(self,path,cmd_define,level=MA5RunningType.PARTON):
+    def __init__(self,path,cmd_define,level=MA5RunningType.PARTON,forced=False):
         self.cmd_define = cmd_define
         self.npart      = 0
         self.path       = path
         self.isopen     = False
         self.level      = level
+        self.forced     = forced
 
     def Load(self):
         if self.level==MA5RunningType.PARTON:
@@ -54,6 +55,8 @@
                 return False
             self.Read()
             self.Close()
+            self.AddSpecialMultiparticles()
+            self.Close()
          
     def OpenPartonLevel(self):        
 
@@ -160,7 +163,7 @@
             logging.debug("Extracting a multiparticle labelled by ["+split[0]+"] with a PDG-id : "+str(split[2:]))
             
             #feed multiparticle
-            self.cmd_define.fill(split[0],split[2:])
+            self.cmd_define.fill(split[0],split[2:],self.forced)
             self.npart += 1
 
 
@@ -172,6 +175,7 @@
                 text += item + " "
             logging.error(text)    
 
+
     def AddSpecialMultiparticles(self):
         # Creating invisible and hadronic multiparticles (mandatory in parton and hadron level)
         if self.level==MA5RunningType.RECO:
 
=== modified file 'madanalysis/IOinterface/particle_reader.py'
--- madanalysis/IOinterface/particle_reader.py	2013-05-17 06:19:22 +0000
+++ madanalysis/IOinterface/particle_reader.py	2013-06-13 12:22:26 +0000
@@ -29,12 +29,13 @@
 
 class ParticleReader():
 
-    def __init__(self,path,cmd_define,level=MA5RunningType.PARTON):
+    def __init__(self,path,cmd_define,level=MA5RunningType.PARTON,forced=False):
         self.cmd_define = cmd_define
-        self.npart  = 0
-        self.path   = path
-        self.level  = level
-        self.isopen = False
+        self.npart      = 0
+        self.path       = path
+        self.level      = level
+        self.isopen     = False
+        self.forced     = forced
 
     def Load(self):
         if self.level==MA5RunningType.PARTON:
@@ -162,7 +163,7 @@
             logging.debug("Extracting a particle labelled by ["+split[1]+"] with a PDG-id : "+split[0])
             
             #feed particle
-            self.cmd_define.fill(split[1],[split[0]])
+            self.cmd_define.fill(split[1],[split[0]],self.forced)
             self.npart += 1
 
     @staticmethod        
 
=== modified file 'madanalysis/UpdateNotes.txt'
--- madanalysis/UpdateNotes.txt	2013-05-17 16:41:08 +0000
+++ madanalysis/UpdateNotes.txt	2013-06-13 12:22:26 +0000
@@ -1,5 +1,54 @@
 Update notes for MadAnalysis 5 (in reverse time order)
 
+76   1.1.25  (2013/06/13) econte:  fixing one bug related to progress bar occuring on some plateforms 
+
+75   1.1.24  (2013/06/12) econte:  adding compilation service + remove warnings 
+
+74   1.1.23  (2013/06/12) econte:  adding MA5/SampleAnalyzer version and date to LHE and LHCO output 
+
+73   1.1.22  (2013/06/12) econte:  change operation order in ma5commit 
+
+72   1.1.21  (2013/06/12) econte:  fixing bugs with label definition in RECO mode + adding version 
+                                   number to SampleAnalyzer 
+
+71   1.1.20  (2013/06/12) econte:  fix bug with usage of zlib from tools folder + release/date stamp 
+                                   inside SampleAnalyzer 
+
+70   1.1.19  (2013/06/12) econte:  small fix related to the progress bar display 
+
+69   1.1.18  (2013/06/12) econte:  fixing bugs (zlib + numpy) + display progress bar 
+
+70   1.1.19  (2013/06/12) econte:  fixing bugs (numpy, zlib) + progress bar display 
+
+69   1.1.18  (2013/06/12) econte:  fixing bugs with numpy + fixing bugs with gzip + fixing bugs with 
+                                   Progress bar [to validate] 
+
+68   1.1.17  (2013/06/10) bfuks:   Bug fixes: progress bar on macosx + empty histogram display 
+
+67   1.1.16  (2013/06/01) econte:  2 files missed in the previous commit 
+
+66   1.1.15  (2013/05/31) econte:  finalizing jet-clustering with exclusiveId option and full 
+                                   compatibility with Pythia8 
+
+65   1.1.14  (2013/05/26) econte:  validation LHCO format 
+
+64   1.1.13  (2013/05/24) econte:  improving jetclustering + adding photons in NPID plot + nb of read 
+                                   events + surprise 
+
+63   1.1.12  (2013/05/21) econte:  first implementation of the exclusive id 
+
+62   1.1.11  (2013/05/21) econte:  adding changes to LHCOWriter + fixing warning appearing in Adam's 
+                                   computer 
+
+61   1.1.10  (2013/05/19) econte:  fixing bug with LHCO reader + implementing LHCO output (only 
+                                   kinematics part is validated) + beautiful headers for LHCO + LHE 
+
+60   1.1.9   (2013/05/18) econte:  test 
+
+59   1.1.8   (2013/05/17) bfuks:   test 
+
+58   1.1.7   (2013/05/17) bfuks:   initializating 1.1.7 
+
 54   1.1.6   (2013/05/17) ma5team:  - Many bug fixes (special thanks to Adam Alloul and Jose Ruiz).
                                     - Simplified LHE output added, LHCO writer at the beta level.
                                     - Clustering and b, c and tau (mis)ident. improved.
 
=== modified file 'madanalysis/configuration/beauty_identification.py'
--- madanalysis/configuration/beauty_identification.py	2013-05-17 06:19:22 +0000
+++ madanalysis/configuration/beauty_identification.py	2013-06-13 12:22:26 +0000
@@ -26,7 +26,7 @@
 class BeautyIdentification():
 
     default_matching_dr = 0.5
-    default_exclusive   = False
+    default_exclusive   = True
     default_efficiency  = 1.
     default_misid_cjet  = 0.
     default_misid_ljet  = 0.
 
=== modified file 'madanalysis/configuration/clustering_configuration.py'
--- madanalysis/configuration/clustering_configuration.py	2013-05-17 06:19:22 +0000
+++ madanalysis/configuration/clustering_configuration.py	2013-06-13 12:22:26 +0000
@@ -50,7 +50,7 @@
         self.algorithm    = "none"
         self.beauty       = BeautyIdentification()
         self.tau          = TauIdentification()
-        self.exclusive_id = False
+        self.exclusive_id = True
 
         
     def Display(self):
@@ -189,7 +189,6 @@
         if parameter=="exclusive_id":
             if value=="true":
                 self.exclusive_id=True
-                logging.warning("This functionnality is not available yet.")
             elif value=="false":
                 self.exclusive_id=False
             else:
 
=== modified file 'madanalysis/configuration/tau_identification.py'
--- madanalysis/configuration/tau_identification.py	2013-05-17 06:19:22 +0000
+++ madanalysis/configuration/tau_identification.py	2013-06-13 12:22:26 +0000
@@ -26,41 +26,28 @@
 class TauIdentification():
 
     default_matching_dr = 0.5
-    default_exclusive   = False
+    default_exclusive   = True
     default_efficiency  = 1.
     default_misid_ljet  = 0.
-    
 
-    userVariables = { "tau_id.matching_dr" : [str(default_matching_dr)],\
-                      "tau_id.exclusive"   : [str(default_exclusive)],\
-                      "tau_id.efficiency"  : [str(default_efficiency)],\
+    userVariables = { "tau_id.efficiency"  : [str(default_efficiency)],\
                       "tau_id.misid_ljet"  : [str(default_misid_ljet)]\
                     }
 
+
     def __init__(self):
-        self.matching_dr = TauIdentification.default_matching_dr
-        self.exclusive   = TauIdentification.default_exclusive
         self.efficiency  = TauIdentification.default_efficiency
         self.misid_ljet  = TauIdentification.default_misid_ljet
 
         
     def Display(self):
         logging.info("  + hadronic-tau identification:")
-        self.user_DisplayParameter("tau_id.matching_dr")
-        self.user_DisplayParameter("tau_id.exclusive")
         self.user_DisplayParameter("tau_id.efficiency")
         self.user_DisplayParameter("tau_id.misid_ljet")
 
 
     def user_DisplayParameter(self,parameter):
-        if parameter=="tau_id.matching_dr":
-            logging.info("    + DeltaR matching = "+str(self.matching_dr))
-        elif parameter=="tau_id.exclusive":
-            msg="false"
-            if self.exclusive:
-                msg="true"
-            logging.info("    + exclusive algo = "+msg)
-        elif parameter=="tau_id.efficiency":
+        if parameter=="tau_id.efficiency":
             logging.info("    + id efficiency = "+str(self.efficiency))
         elif parameter=="tau_id.misid_ljet":
             logging.info("    + mis-id efficiency (light quarks) = "+str(self.misid_ljet))
@@ -70,12 +57,7 @@
 
     def SampleAnalyzerConfigString(self):
         mydict = {}
-        mydict['tau_id.matching_dr'] = str(self.matching_dr)
         mydict['tau_id.efficiency']  = str(self.efficiency)
-        if self.exclusive:
-            mydict['tau_id.exclusive'] = '1'
-        else:
-            mydict['tau_id.exclusive'] = '0'
         mydict['tau_id.misid_ljet']  = str(self.misid_ljet)
         return mydict
 
@@ -92,30 +74,8 @@
 
 
     def user_SetParameter(self,parameter,value):
-        # matching deltar
-        if parameter=="tau_id.matching_dr":
-            try:
-                number = float(value)
-            except:
-                logging.error("the 'matching deltaR' must be a float value.")
-                return False
-            if number<=0:
-                logging.error("the 'matching deltaR' cannot be negative or null.")
-                return False
-            self.matching_dr=number
-
-        # exclusive
-        elif parameter=="tau_id.exclusive":
-            if value == "true":
-                self.exclusive=True
-            elif value == "false":
-                self.exclusive=False
-            else:
-                logging.error("'exclusive' possible values are : 'true', 'false'")
-                return False
-
         # efficiency
-        elif parameter=="tau_id.efficiency":
+        if parameter=="tau_id.efficiency":
             try:
                 number = float(value)
             except:
 
=== modified file 'madanalysis/core/config_checker.py'
--- madanalysis/core/config_checker.py	2013-05-17 06:19:22 +0000
+++ madanalysis/core/config_checker.py	2013-06-13 12:22:26 +0000
@@ -314,6 +314,23 @@
                 os.environ['LIBRARY_PATH'] = os.environ['LIBRARY_PATH'] + \
                     		        ":" + item
 	        break
+
+        # If not, test if it is there locally 
+        if not find: 
+            item = self.ma5dir+'/tools/zlib/lib/'
+            files=glob.glob(item+"/libz.so")
+            files.extend(glob.glob(item+"/libz.a"))
+            files.extend(glob.glob(item+"/libz.dylib"))
+            if len(files)!=0:
+	        self.configLinux.libraries['ZLib']=files[0]+":"+str(os.stat(files[0]).st_mtime)
+   	        find=True
+                os.environ['LD_LIBRARY_PATH'] = os.environ['LD_LIBRARY_PATH'] + \
+                    		        ":" + item
+                os.environ['DYLD_LIBRARY_PATH'] = os.environ['DYLD_LIBRARY_PATH'] + \
+                    		        ":" + item
+                os.environ['LIBRARY_PATH'] = os.environ['LIBRARY_PATH'] + \
+                    		        ":" + item
+        # Not fail -> Warning
         if not find:
             self.PrintFAIL(warning=True)
 	    logging.warning("Library called 'libz' not found. Gzip format will be disabled.")
@@ -327,9 +344,16 @@
             if len(files)!=0:
   	        self.configLinux.headers['ZLib']=files[0]+":"+str(os.stat(files[0]).st_mtime)
 	        find=True
-                os.environ['CPLUS_INCLUDE_PATH'] = os.environ['CPLUS_INCLUDE_PATH'] + \
-                    		        ":" + item
+                os.environ['CPLUS_INCLUDE_PATH'] = os.environ['CPLUS_INCLUDE_PATH'] + ":" + item
    	        break
+
+        # If not, test if it is there locally 
+        if not find: 
+            find=os.path.isfile(self.ma5dir+'/tools/zlib/include/zlib.h')
+            if find:
+                os.environ['CPLUS_INCLUDE_PATH'] = os.environ['CPLUS_INCLUDE_PATH'] + ":" + \
+                                                   self.ma5dir+'/tools/zlib/include/'
+
         if not find:
             self.PrintFAIL(warning=True)
 	    logging.warning("Header file called 'zlib.h' not found. Gzip format will be disabled.")
@@ -369,6 +393,7 @@
             if(len(files))!=0:
                find = True
                self.configLinux.fastjet_version = commands.getstatusoutput('fastjet-config --version')[1]
+
         # If not, test if it is there locally 
         if not find: 
             find=os.path.isfile(self.ma5dir+'/tools/fastjet/bin/fastjet-config')
 
=== modified file 'madanalysis/core/expert_mode.py'
--- madanalysis/core/expert_mode.py	2013-05-17 06:19:22 +0000
+++ madanalysis/core/expert_mode.py	2013-06-13 12:22:26 +0000
@@ -35,7 +35,6 @@
     def __init__(self,main):
         self.main=main
         self.path=""
-        self.dir=""
         self.forbiddenpaths=[]
         self.forbiddenpaths.append(os.path.normpath(self.main.ma5dir+'/lib'))
         self.forbiddenpaths.append(os.path.normpath(self.main.ma5dir+'/bin'))
@@ -45,8 +44,10 @@
         logging.info("\nWelcome to the expert mode of MadAnalysis")
         logging.info("Please enter a folder for creating an empty SampleAnalyzer job")
         answer=raw_input("Answer: ")
-        self.dir = answer
-        self.path = os.path.normpath(self.main.currentdir+'/'+answer)
+        self.path = os.path.expanduser(answer)
+        if not self.path.startswith('/'):
+            self.path = self.main.currentdir+'/'+self.path
+        self.path = os.path.normpath(self.path)
 
         # Checking folder
         if self.path in self.forbiddenpaths:
 
=== modified file 'madanalysis/core/main.py'
--- madanalysis/core/main.py	2013-05-17 06:19:22 +0000
+++ madanalysis/core/main.py	2013-06-13 12:22:26 +0000
@@ -50,10 +50,13 @@
     userVariables = { "currentdir"      : [], \
                       "normalize"       : ["none","lumi","lumi_weight"], \
                       "lumi"            : [], \
-                      "SBratio"         : ['"S/B"','"B/S"','"S/(S+B)"','"B/(B+S)"','"S/sqrt(S+B)"','"B/sqrt(B+S)"'], \
+                      "SBratio"         : ['"S/B"','"B/S"',\
+                                           '"S/(S+B)"','"B/(B+S)"',\
+                                           '"S/sqrt(S+B)"','"B/sqrt(B+S)"'], \
                       "SBerror"         : [], \
                       "stacking_method" : ["stack","superimpose","normalize2one"], \
-                      "outputfile"      : ['"output.lhe.gz"'] }
+                      "outputfile"      : ['"output.lhe.gz"','"output.lhco.gz"'] \
+                      }
 
     SBformula = { 'S/B'         : '1./(B**2)*sqrt(B**2*ES**2+S**2*EB**2)', \
                   'S/(S+B)'     : '1./(S+B)**2*sqrt(B**2*ES**2+S**2*EB**2)', \
@@ -291,11 +294,41 @@
                 quoteTag=True
             if quoteTag:
                 value=value[1:-1]
-            if not value.lower().endswith(".lhe") and not value.lower().endswith(".lhe.gz"):
-                logging.error("LHE format (*.lhe) and compressed LHE format (*.lhe.gz) are the only possible formats for output.")
-                return False
-            self.output = value
-            return
+            valuemin = value.lower()
+
+            # Compressed file
+            if valuemin.endswith(".gz") and not self.libZIP:
+                logging.error("Compressed formats (*.gz) are not available. "\
+                              + "Please install zlib with the command line:")
+                logging.error(" install zlib")
+                return False
+
+            # LHE
+            if valuemin.endswith(".lhe") or valuemin.endswith(".lhe.gz"):
+                self.output = value
+                return
+
+            # LHCO
+            elif valuemin.endswith(".lhco") or valuemin.endswith(".lhco.gz"):
+                if self.mode == MA5RunningType.RECO:
+                    self.output = value
+                    return
+                elif self.mode == MA5RunningType.PARTON:
+                    logging.error("LHCO format is not available in PARTON mode.")
+                    return False
+                elif self.mode == MA5RunningType.HADRON:
+                    if self.clustering.algorithm == "none":
+                        logging.error("Please select a jet-clustering algorithm before requesting a LHCO file output.")
+                        logging.error("Command: set main.merging.algorithm = ")
+                        return False
+                    else:
+                        self.output = value
+                        return
+
+            else:
+                logging.error("Output format is not available. Extension allowed: " +\
+                              ".lhe .lhe.gz .lhco .lhco.gz")
+                return False
 
         # other    
         else:
 
=== modified file 'madanalysis/input/reco_default.txt'
--- madanalysis/input/reco_default.txt	2013-04-14 20:56:15 +0000
+++ madanalysis/input/reco_default.txt	2013-06-13 12:22:26 +0000
@@ -7,7 +7,8 @@
  15 ta-
 -15 ta+
  21 j
+ 22 a
   5 b
   1 nb
 100 met
-101 mht
 \ No newline at end of file
+101 mht
 
=== modified file 'madanalysis/interpreter/cmd_import.py'
--- madanalysis/interpreter/cmd_import.py	2013-05-17 06:19:22 +0000
+++ madanalysis/interpreter/cmd_import.py	2013-06-13 12:22:26 +0000
@@ -163,7 +163,7 @@
         # Loading particles
         input = ParticleReader(self.main.ma5dir,cmd_define,self.main.mode)
         input.Load()
-        input = MultiparticleReader(self.main.ma5dir,cmd_define,self.main.mode)
+        input = MultiparticleReader(self.main.ma5dir,cmd_define,self.main.mode,self.main.forced)
         input.Load()
 
         # Reset history
 
=== modified file 'madanalysis/interpreter/cmd_reset.py'
--- madanalysis/interpreter/cmd_reset.py	2013-05-17 06:19:22 +0000
+++ madanalysis/interpreter/cmd_reset.py	2013-06-13 12:22:26 +0000
@@ -85,7 +85,7 @@
         # Loading particles
         input = ParticleReader(self.main.ma5dir,cmd_define,self.main.mode)
         input.Load()
-        input = MultiparticleReader(self.main.ma5dir,cmd_define,self.main.mode)
+        input = MultiparticleReader(self.main.ma5dir,cmd_define,self.main.mode,self.main.forced)
         input.Load()
         
 
 
=== modified file 'madanalysis/interpreter/interpreter.py'
--- madanalysis/interpreter/interpreter.py	2013-05-17 06:19:22 +0000
+++ madanalysis/interpreter/interpreter.py	2013-06-13 12:22:26 +0000
@@ -282,7 +282,7 @@
         input.Load()
 
     def InitializeMultiparticle(self):
-        input = MultiparticleReader(self.main.ma5dir,self.cmd_define,self.main.mode)
+        input = MultiparticleReader(self.main.ma5dir,self.cmd_define,self.main.mode,self.main.forced)
         input.Load()
 
     # PreLoop
@@ -474,5 +474,5 @@
         if line[-1] in [' ','\t']:
             return '', line, len(line),len(enidx)
         return text, line, begidx, endidx
-        
+
         
 
=== modified file 'madanalysis/interpreter/interpreter_base.py'
--- madanalysis/interpreter/interpreter_base.py	2013-05-17 06:19:22 +0000
+++ madanalysis/interpreter/interpreter_base.py	2013-06-13 12:22:26 +0000
@@ -103,6 +103,12 @@
         # --> removing additionnal whitespace characters
         line = line.lstrip()
 
+        # pattern design
+        if len(line)==4 and \
+           line[0]=='m' and line[1]=='u' and line[2]=='f' and line[3]=='!':
+            self.pattern_design()
+            return ''
+
         # Add the line to the history
         # except for useless commands (empty history, help calls and comments)
         if line != "history" and \
@@ -334,3 +340,44 @@
             else:
                 out.append(data)
         return out
+
+
+    def pattern_design(self):
+        pattern=[]
+        pattern.append('32-32-32-32-32-32-32-32-95-95-95-95-95-32')
+        pattern.append('32-32-32-32-32-95-45-126-126-32-32-32-32-32-126-126-45-95-47-47-32')
+        pattern.append('32-32-32-47-126-32-32-32-32-32-32-32-32-32-32-32-32-32-126-92-32')
+        pattern.append('32-32-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-95-32-32-124-95-32')
+        pattern.append('32-124-32-32-32-32-32-32-32-32-32-95-45-45-126-126-126-32-41-126-126-32-41-95-95-95-32')
+        pattern.append('92-124-32-32-32-32-32-32-32-32-47-32-32-32-95-95-95-32-32-32-95-45-126-32-32-32-126-45-39-95-32')
+        pattern.append('92-32-32-32-32-32-32-32-32-32-32-95-45-126-32-32-32-126-45-95-32-32-32-32-32-32-32-32-32-92-32')
+        pattern.append('124-32-32-32-32-32-32-32-32-32-47-32-32-32-32-32-32-32-32-32-92-32-32-32-32-32-32-32-32-32-124-32')
+        pattern.append('124-32-32-32-32-32-32-32-32-124-32-32-32-32-32-32-32-32-32-32-32-124-32-32-32-32-32-40-79-32-32-124-32')
+        pattern.append('32-124-32-32-32-32-32-32-124-32-32-32-32-32-32-32-32-32-32-32-32-32-124-32-32-32-32-32-32-32-32-124-32')
+        pattern.append('32-32-124-32-32-32-32-32-32-124-32-32-32-79-41-32-32-32-32-32-32-32-32-124-32-32-32-32-32-32-32-124-32')
+        pattern.append('32-32-47-124-32-32-32-32-32-32-124-32-32-32-32-32-32-32-32-32-32-32-124-32-32-32-32-32-32-32-47-32')
+        pattern.append('32-32-47-32-92-32-95-45-45-95-32-92-32-32-32-32-32-32-32-32-32-47-45-95-32-32-32-95-45-126-41-32')
+        pattern.append('32-32-32-32-47-126-32-32-32-32-92-32-126-45-95-32-32-32-95-45-126-32-32-32-126-126-126-95-95-47-32')
+        pattern.append('32-32-32-124-32-32-32-124-92-32-32-126-45-95-32-126-126-126-32-95-45-126-126-45-45-45-126-32-32-92-32')
+        pattern.append('32-32-32-124-32-32-32-124-32-124-32-32-32-32-126-45-45-126-126-32-32-47-32-92-32-32-32-32-32-32-126-45-39-95-32')
+        pattern.append('32-32-32-124-32-32-32-92-32-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-126-45-39-95-32')
+        pattern.append('32-32-32-32-92-32-32-32-126-45-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-126-126-45-45-95-95-32-95-45-126-126-45-44-32')
+        pattern.append('32-32-32-32-32-126-45-95-32-32-32-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-47-32-32-32-32-32-124-32')
+        pattern.append('32-32-32-32-32-32-32-32-126-126-45-45-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-47-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-47-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-32-124-32-32-32-32-32-32-32-32-32-32-32-32-32-32-95-32-32-32-32-32-32-32-32-32-32-32-32-95-45-126-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-47-126-126-45-45-95-32-32-32-95-95-45-45-45-126-126-32-32-32-32-32-32-32-32-32-32-95-45-126-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-92-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-32-95-95-45-45-126-126-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-124-126-126-45-45-95-95-32-32-32-32-32-95-95-95-45-45-45-126-126-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-124-32-32-32-32-32-32-126-126-126-126-126-32')
+        pattern.append('32-32-32-32-32-32-32-32-32-32-124-32-32-124-32')
+
+        for word in pattern:
+            msg=""
+            words = word.split('-')
+            for i in range(len(words)):
+                msg+=chr((int(words[i])))
+            logging.info("\x1b[1m"+"\x1b[32m"+msg+"\x1b[0m")
+        logging.info("")
+
+
 
=== modified file 'madanalysis/job/job_execute.py'
--- madanalysis/job/job_execute.py	2013-05-17 06:19:22 +0000
+++ madanalysis/job/job_execute.py	2013-06-13 12:22:26 +0000
@@ -162,6 +162,24 @@
                    container+'.push_back(&(event.rec()->electrons()[i]));\n')
 
 
+def WriteFillWithPhotonContainer(part,file,rank,status):
+    
+    # If PTrank, no fill
+    if part.PTrank!=0:
+        return
+
+    # Skipping if already defined
+    if InstanceName.Find('P_'+part.name+rank+status):
+        return
+
+    # Getting container name
+    container=InstanceName.Get('P_'+part.name+rank+status)
+
+    # Put photon
+    if part.particle.Find(22):
+        file.write('      '+container+'.push_back(&(event.rec()->photons()[i]));\n')
+
+
 def WriteFillWithMuonContainer(part,file,rank,status):
     
     # If PTrank, no fill
@@ -299,6 +317,14 @@
         file.write('    }\n')
         InstanceName.Clear()
 
+        # Filling with photons
+        file.write('    for (UInt_t i=0;i<event.rec()->photons().size();i++)\n')
+        file.write('    {\n')
+        for item in part_list:
+            WriteFillWithPhotonContainer(item[0],file,item[1],item[2])
+        file.write('    }\n')
+        InstanceName.Clear()
+
         # Filling with electrons
         file.write('    for (UInt_t i=0;i<event.rec()->electrons().size();i++)\n')
         file.write('    {\n')
 
=== modified file 'madanalysis/job/job_initialize.py'
--- madanalysis/job/job_initialize.py	2013-05-17 06:19:22 +0000
+++ madanalysis/job/job_initialize.py	2013-06-13 12:22:26 +0000
@@ -45,13 +45,390 @@
     file.write('{\n')
 
     # mcConfig initialization
-    #if main.mode!=MA5RunningType.RECO: FOR ALL MODES !!!!
-    file.write('  // Initializing PhysicsService for MC\n') 
-    file.write('  PHYSICS->mcConfig().Reset();\n\n')
-    WriteHadronicList(file,main)
-    file.write('\n')
-    WriteInvisibleList(file,main)
-    file.write('\n')
+    if main.mode!=MA5RunningType.RECO:
+        file.write('  // Initializing PhysicsService for MC\n') 
+        file.write('  PHYSICS->mcConfig().Reset();\n\n')
+        WriteHadronicList(file,main)
+        file.write('\n')
+        WriteInvisibleList(file,main)
+        file.write('\n')
+    else:
+        file.write('  // Initializing PhysicsService for MC\n') 
+        file.write('  PHYSICS->mcConfig().Reset();\n\n')
+        file.write('\n')
+        file.write('  // definition of the multiparticle "hadronic"\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20543);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20533);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20523);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20513);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20433);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20423);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20413);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20323);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20313);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-20213);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10543);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10541);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10533);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10531);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10523);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10521);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10513);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10511);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10433);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10431);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10423);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10421);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10413);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10411);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10323);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10321);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10313);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10311);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10213);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-10211);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5554);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5544);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5542);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5534);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5532);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5524);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5522);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5514);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5512);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5503);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5444);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5442);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5434);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5432);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5424);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5422);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5414);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5412);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5403);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5401);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5342);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5334);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5332);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5324);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5322);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5314);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5312);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5303);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5301);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5242);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5232);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5222);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5201);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5142);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5132);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5122);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-5101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4444);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4434);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4432);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4424);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4422);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4414);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4412);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4403);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4334);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4332);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4324);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4322);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4314);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4312);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4303);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4301);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4232);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4222);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4201);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4132);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4122);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-4101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3334);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3324);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3322);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3314);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3312);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3303);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3222);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3201);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3122);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-3101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-2101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-1114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-1103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-545);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-543);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-541);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-535);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-533);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-531);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-525);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-523);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-521);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-515);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-513);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-511);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-435);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-433);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-431);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-425);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-423);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-421);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-415);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-413);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-411);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-325);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-323);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-321);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-315);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-313);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-311);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-215);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-213);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(-211);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(111);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(113);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(115);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(130);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(211);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(213);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(215);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(221);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(223);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(225);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(310);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(311);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(313);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(315);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(321);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(323);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(325);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(331);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(333);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(335);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(411);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(413);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(415);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(421);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(423);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(425);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(431);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(433);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(435);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(441);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(443);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(445);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(511);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(513);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(515);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(521);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(523);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(525);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(531);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(533);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(535);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(541);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(543);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(545);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(551);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(553);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(555);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(1103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(1114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(2224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3122);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3201);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3222);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3303);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3312);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3314);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3322);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3324);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(3334);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4122);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4132);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4201);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4222);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4232);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4301);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4303);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4312);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4314);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4322);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4324);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4332);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4334);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4403);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4412);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4414);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4422);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4424);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4432);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4434);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(4444);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5101);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5103);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5112);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5114);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5122);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5132);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5142);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5201);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5203);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5212);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5214);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5222);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5224);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5232);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5242);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5301);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5303);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5312);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5314);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5322);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5324);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5332);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5334);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5342);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5401);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5403);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5412);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5414);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5422);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5424);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5432);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5434);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5442);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5444);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5503);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5512);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5514);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5522);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5524);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5532);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5534);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5542);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5544);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(5554);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10111);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10113);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10211);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10213);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10221);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10223);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10311);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10313);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10321);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10323);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10331);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10333);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10411);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10413);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10421);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10423);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10431);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10433);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10441);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10443);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10511);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10513);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10521);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10523);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10531);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10533);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10541);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10543);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10551);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(10553);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20113);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20213);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20223);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20313);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20323);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20333);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20413);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20423);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20433);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20443);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20513);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20523);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20533);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20543);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(20553);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(100443);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(100553);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9900440);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9900441);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9900443);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9900551);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9900553);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9910441);\n')
+        file.write('  PHYSICS->mcConfig().AddHadronicId(9910551);\n')
+        file.write('\n')
+        file.write('  // definition of the multiparticle "invisible"\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(-16);\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(-14);\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(-12);\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(12);\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(14);\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(16);\n')
+        file.write('  PHYSICS->mcConfig().AddInvisibleId(1000022);\n')
+        file.write('\n')
+
 
     # recConfig initialization
     if main.mode==MA5RunningType.RECO:
 
=== modified file 'madanalysis/job/job_plot.py'
--- madanalysis/job/job_plot.py	2013-05-17 06:19:22 +0000
+++ madanalysis/job/job_plot.py	2013-06-13 12:22:26 +0000
@@ -88,6 +88,12 @@
     # RECO mode
     else:
 
+        # photons
+        file.write('  for (unsigned int i=0;i<event.rec()->photons().size();i++)\n')
+        file.write('  {\n')
+        file.write('    H'+str(ihisto)+'_->Fill(22,__event_weight__);\n')
+        file.write('  }\n')
+
         # electrons
         file.write('  for (unsigned int i=0;i<event.rec()->electrons().size();i++)\n')
         file.write('  {\n')
 
=== modified file 'madanalysis/layout/histogram.py'
--- madanalysis/layout/histogram.py	2013-05-17 06:19:22 +0000
+++ madanalysis/layout/histogram.py	2013-06-13 12:22:26 +0000
@@ -24,7 +24,6 @@
 
 from madanalysis.layout.histogram_core import HistogramCore
 import logging
-import numpy
 
 
 class Histogram:
@@ -49,6 +48,8 @@
 
     def FinalizeReading(self,main,dataset):
 
+        import numpy
+
         # Statistics
         self.summary.nevents   = self.positive.nevents   + self.negative.nevents
         self.summary.nentries   = self.positive.nentries   + self.negative.nentries
 
=== modified file 'madanalysis/layout/histogram_frequency_core.py'
--- madanalysis/layout/histogram_frequency_core.py	2013-05-17 06:19:22 +0000
+++ madanalysis/layout/histogram_frequency_core.py	2013-06-13 12:22:26 +0000
@@ -23,11 +23,11 @@
 
 
 import logging
-import numpy
 
 class HistogramFrequencyCore:
 
     def __init__(self):
+        import numpy
         self.integral  = 0
         self.nevents   = 0
         self.entries   = 0
 
=== modified file 'madanalysis/layout/plotflow.py'
--- madanalysis/layout/plotflow.py	2013-05-17 06:19:22 +0000
+++ madanalysis/layout/plotflow.py	2013-06-13 12:22:26 +0000
@@ -164,7 +164,6 @@
 
         # Loop over datasets and histos
         for ind in range(0,len(histos)):
-
             # Scaling 
             histos[ind].Scale(scales[ind])
             
@@ -396,7 +395,7 @@
         stack = THStack("mystack","")
         ntot = 0
         for item in histos:
-            ntot+=item.GetEntries()
+            ntot+=item.Integral()
             stack.Add(item)
 
         # Drawing
 
=== modified file 'tools/SampleAnalyzer/Analyzer/AnalyzerBase.h'
--- tools/SampleAnalyzer/Analyzer/AnalyzerBase.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Analyzer/AnalyzerBase.h	2013-06-13 12:22:26 +0000
@@ -94,10 +94,10 @@
 
   /// Initialize (common part to all analyses)
   bool PreInitialize(const std::string& outputName, 
-                     bool weighted_events=false)
+                     const Configuration* cfg)
   {
-    weighted_events_ = weighted_events;
-    out_.Initialize(outputName.c_str());
+    weighted_events_ = cfg->IsNoEventWeight();
+    out_.Initialize(cfg,outputName.c_str());
     return true;
   }
 
 
=== modified file 'tools/SampleAnalyzer/Analyzer/MergingPlots.cpp'
--- tools/SampleAnalyzer/Analyzer/MergingPlots.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Analyzer/MergingPlots.cpp	2013-06-13 12:22:26 +0000
@@ -27,6 +27,7 @@
 #include "SampleAnalyzer/Analyzer/MergingPlots.h"
 #include "SampleAnalyzer/Core/Configuration.h"
 #include "SampleAnalyzer/Service/Physics.h"
+#include "SampleAnalyzer/Service/CompilationService.h"
 
 //STL headers
 #include <sstream>
@@ -201,6 +202,7 @@
 void MergingPlots::ExtractDJRwithFortran(const std::vector<fastjet::PseudoJet>& inputs,std::vector<Double_t>& DJRvalues)
 {
   double PP[512][4];
+  UNUSED(PP);
   for (unsigned int i=0;i<inputs.size();i++)
   {
     PP[i][0]=inputs[i].px();
 
=== modified file 'tools/SampleAnalyzer/Core/Configuration.cpp'
--- tools/SampleAnalyzer/Core/Configuration.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Core/Configuration.cpp	2013-06-13 12:22:26 +0000
@@ -36,6 +36,14 @@
 
 
 // -----------------------------------------------------------------------------
+// Initializing static data members
+// -----------------------------------------------------------------------------
+// DO NOT TOUCH THESE LINES
+const std::string Configuration::sampleanalyzer_version_ = "1.1.25";
+const std::string Configuration::sampleanalyzer_date_    = "2013/06/13";
+// DO NOT TOUCH THESE LINES
+
+// -----------------------------------------------------------------------------
 // PrintSyntax
 // -----------------------------------------------------------------------------
 void Configuration::PrintSyntax()
@@ -94,6 +102,22 @@
     // weighted event
     else if (option=="--no_event_weight") no_event_weight_ = true;
 
+    // version
+    else if (option.find("--ma5_version=")==0)
+    {
+      std::string stamp = option.substr(14,std::string::npos);
+      std::size_t result = stamp.find(";");
+      if (result==std::string::npos)
+      {
+        WARNING << "MA5 version '" << stamp << "' is not valid." << std::endl;
+      }
+      else
+      {
+        pythoninterface_version_ = stamp.substr(0,result);
+        pythoninterface_date_ = stamp.substr(result+1,std::string::npos);
+      }
+    }
+
     // unknown option
     else
     {
@@ -116,6 +140,13 @@
 // -----------------------------------------------------------------------------
 void Configuration::Display()
 {
+  INFO << "      - version: " << sampleanalyzer_version_ << " (" << sampleanalyzer_date_ << ") ";
+  if ((sampleanalyzer_version_!=pythoninterface_version_ && pythoninterface_version_!="") || 
+      (sampleanalyzer_date_!=pythoninterface_date_ && pythoninterface_version_!=""))
+    INFO << "[ python interface version: " << pythoninterface_version_ 
+         << " (" << pythoninterface_date_ << ") ]";
+  INFO << endmsg;
+ 
   INFO << "      - general: ";
 
   // Is there option ?
 
=== modified file 'tools/SampleAnalyzer/Core/Configuration.h'
--- tools/SampleAnalyzer/Core/Configuration.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Core/Configuration.h	2013-06-13 12:22:26 +0000
@@ -44,6 +44,14 @@
   // -------------------------------------------------------------
   private:
 
+    /// SampleAnalyzer version
+    static const std::string sampleanalyzer_version_;
+    static const std::string sampleanalyzer_date_;  
+
+    /// Python interface version
+    std::string pythoninterface_version_;
+    std::string pythoninterface_date_;
+
     /// option : check event mode
     Bool_t check_event_;
 
@@ -76,6 +84,22 @@
     /// Help message
     void PrintSyntax();
 
+    /// Accessor to the sampleanalyzer date
+    const std::string& GetSampleAnalyzerDate() const
+    {return sampleanalyzer_date_;}
+
+    /// Accessor to the sampleanalyzer version
+    const std::string& GetSampleAnalyzerVersion() const
+    {return sampleanalyzer_version_;}
+
+    /// Accessor to the python interface date
+    const std::string& GetPythonInterfaceDate() const
+    {return sampleanalyzer_date_;}
+
+    /// Accessor to the python interface version
+    const std::string& GetPythonInterfaceVersion() const
+    {return sampleanalyzer_version_;}
+
     /// Accessor to the input name
     const std::string& GetInputFileName() const
     {return input_list_name_;}
 
=== modified file 'tools/SampleAnalyzer/Core/ProgressBar.cpp'
--- tools/SampleAnalyzer/Core/ProgressBar.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Core/ProgressBar.cpp	2013-06-13 12:22:26 +0000
@@ -23,6 +23,8 @@
 
 
 // STL headers
+#include <cstdlib>
+
 // SampleAnalyzer headers
 #include "SampleAnalyzer/Core/ProgressBar.h"
 #include "SampleAnalyzer/Service/LogStream.h"
@@ -30,7 +32,7 @@
 using namespace MA5;
 
 
-const std::string ProgressBar::header("        => progress ");
+const std::string ProgressBar::header("        => progress: ");
 
 
 // -----------------------------------------------------------------------------
@@ -63,14 +65,22 @@
   //  std::cout << "eric Indicator_=" << Indicator_ << std::endl;
 
   /// Saving previous stream buffer related to std::cout
-  oldstreambuf_ = std::cout.rdbuf();
+  oldstreambuf_cout_ = std::cout.rdbuf();
   std::cout.flush();
+  oldstreambuf_cerr_ = std::cerr.rdbuf();
+  std::cerr.flush();
+  oldstreambuf_clog_ = std::clog.rdbuf();
+  std::clog.flush();
 
   /// Assigning a new stream buffer with spy
-  newstreambuf_ = new SpyStreamBuffer(oldstreambuf_);
-  std::cout.rdbuf(newstreambuf_);
-
-  //  std::cout << "muf Indicator_=" << Indicator_ << std::endl;
+  newstreambuf_cout_ = new SpyStreamBuffer(oldstreambuf_cout_);
+  std::cout.rdbuf(newstreambuf_cout_);
+
+  newstreambuf_cerr_ = new SpyStreamBuffer(oldstreambuf_cerr_);
+  std::cout.rdbuf(newstreambuf_cerr_);
+
+  newstreambuf_clog_ = new SpyStreamBuffer(oldstreambuf_clog_);
+  std::cout.rdbuf(newstreambuf_clog_);
 
 }
 
@@ -120,22 +130,30 @@
   todisplay[ind+1]='>';
   todisplay[Nstep_+1]=']'; // overwrite the character '>' if ind=Nstep_
 
+  bool newline=false;
   // Go back to the line ? 
   if (ind!=0 && !FirstTime_)
   {
-    std::vector<char> last_chars = newstreambuf_->get_last_chars();
-    if (last_chars[0] != '\n') std::cout << std::endl;
-    else if (last_chars[5] == 93 && last_chars[4] == 27 && 
-             last_chars[3] == 91 && last_chars[2] == 48 && 
-             last_chars[1] == 109) std::cout << "\b\r";
+    if (!newstreambuf_cout_->GetProgressBarMode()) newline=true;
+    if (!newstreambuf_cerr_->GetProgressBarMode()) newline=true;
+    if (!newstreambuf_clog_->GetProgressBarMode()) newline=true;
   }
-  FirstTime_=false;
 
   // Adding header
+  unsigned int toremove = (header+todisplay).size();
   todisplay = header + "\x1b[34m"+ todisplay + "\x1b[0m";
 
   // Displaying
-  std::cout << todisplay << std::endl;
+  newstreambuf_cout_->SetProgressBarMode(false);
+  newstreambuf_cerr_->SetProgressBarMode(false);
+  newstreambuf_clog_->SetProgressBarMode(false);
+  if (newline) std::cout << std::endl;
+  else if (!FirstTime_) std::cout << std::string(toremove,'\b');
+  std::cout << todisplay;
+  newstreambuf_cout_->SetProgressBarMode(true);
+  newstreambuf_cerr_->SetProgressBarMode(true);
+  newstreambuf_clog_->SetProgressBarMode(true);
+  FirstTime_=false;
 }
 
 
@@ -149,11 +167,21 @@
 
   // Display the indicator with a status of 100% ?
   if (!MuteEnd_) Display(Nstep_);
+  newstreambuf_cout_->SetProgressBarMode(false);
+  newstreambuf_cerr_->SetProgressBarMode(false);
+  newstreambuf_clog_->SetProgressBarMode(false);
+  std::cout << std::endl;
 
   // Restoring the initial stream buffer
-  std::cout.rdbuf(oldstreambuf_);
-  oldstreambuf_=0;
-  if (newstreambuf_!=0) delete newstreambuf_;
+  std::cout.rdbuf(oldstreambuf_cout_);
+  std::cerr.rdbuf(oldstreambuf_cerr_);
+  std::clog.rdbuf(oldstreambuf_clog_);
+  oldstreambuf_cout_=0;
+  oldstreambuf_cerr_=0;
+  oldstreambuf_clog_=0;
+  if (newstreambuf_cout_!=0) delete newstreambuf_cout_;
+  if (newstreambuf_cerr_!=0) delete newstreambuf_cerr_;
+  if (newstreambuf_clog_!=0) delete newstreambuf_clog_;
 
   // Reset the progress bar
   Reset();
 
=== modified file 'tools/SampleAnalyzer/Core/ProgressBar.h'
--- tools/SampleAnalyzer/Core/ProgressBar.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Core/ProgressBar.h	2013-06-13 12:22:26 +0000
@@ -56,30 +56,36 @@
     public:
 
       /// Constructor
-      SpyStreamBuffer(std::streambuf* buf) : buf(buf)
+    SpyStreamBuffer(std::streambuf* buf) : buf_(buf)
       {
-        last_chars.resize(7,traits_type::eof());
+        add_endl_=false;
         // no buffering, overflow on every char
         setp(0, 0);
       }
 
-      /// Getting last character
-      const std::vector<char>& get_last_chars() const 
-      { return last_chars; }
+      /// Set ProgressBar mode
+      void SetProgressBarMode(bool status=true)
+      { add_endl_=status; }
+
+      /// Accessor to ProgressBar mode status
+      bool GetProgressBarMode() const
+      { return add_endl_; }
 
       /// Overflow method
       virtual int_type overflow(int_type c)
       {
-        buf->sputc(c);
-        for (unsigned int i=(last_chars.size()-1);i>0;i--) 
-           last_chars[i]=last_chars[i-1];
-        last_chars[0] = c;
+        if (add_endl_) 
+        {
+          buf_->sputc('\n');
+          add_endl_=false;
+        }
+        buf_->sputc(c);
         return c;
       }
     private:
 
-      std::streambuf* buf;
-      std::vector<char> last_chars;
+      std::streambuf* buf_;
+      bool add_endl_;
     };
 
   // -------------------------------------------------------------
@@ -112,10 +118,14 @@
     std::vector<Long64_t> Thresholds_;
 
     /// Pointer to the new stream buffer
-    SpyStreamBuffer* newstreambuf_;
+    SpyStreamBuffer* newstreambuf_cout_;
+    SpyStreamBuffer* newstreambuf_cerr_;
+    SpyStreamBuffer* newstreambuf_clog_;
 
     /// Pointer to the old stream buffer
-    std::streambuf* oldstreambuf_;
+    std::streambuf* oldstreambuf_cout_;
+    std::streambuf* oldstreambuf_cerr_;
+    std::streambuf* oldstreambuf_clog_;
 
     static const std::string header;
 
@@ -127,14 +137,22 @@
     /// Constructor without argument
     ProgressBar()
     {
-      newstreambuf_=0;
-      oldstreambuf_=0;
+      newstreambuf_cout_=0;
+      newstreambuf_cerr_=0;
+      newstreambuf_clog_=0;
+      oldstreambuf_cout_=0;
+      oldstreambuf_cerr_=0;
+      oldstreambuf_clog_=0;
       Reset(); 
     }
 
     /// Destructor 
     ~ProgressBar() 
-    { if (newstreambuf_!=0) delete newstreambuf_; }
+    { 
+      if (newstreambuf_cout_!=0) delete newstreambuf_cout_;  
+      if (newstreambuf_cerr_!=0) delete newstreambuf_cerr_; 
+      if (newstreambuf_clog_!=0) delete newstreambuf_clog_; 
+   }
 
     /// Reset
     void Reset()
 
=== modified file 'tools/SampleAnalyzer/Core/SampleAnalyzer.cpp'
--- tools/SampleAnalyzer/Core/SampleAnalyzer.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Core/SampleAnalyzer.cpp	2013-06-13 12:22:26 +0000
@@ -27,15 +27,14 @@
 #include "SampleAnalyzer/Writer/SAFWriter.h"
 #include "SampleAnalyzer/Service/ExceptionService.h"
 #include "SampleAnalyzer/Service/TimeService.h"
-#include "SampleAnalyzer/Service/Terminate.h"
 #include "SampleAnalyzer/Service/PDGService.h"
 #include "SampleAnalyzer/Service/Terminate.h"
+#include "SampleAnalyzer/Service/CompilationService.h"
 
 
 using namespace MA5;
 
-
-/// Constructor withtout arguments
+/// Constructor without arguments
 SampleAnalyzer::SampleAnalyzer() 
 {
   // Initializing service
@@ -49,7 +48,7 @@
   progressBar_=0;
 
   // Header
-  INFO << "    * SampleAnalyzer 2.0 for MadAnalysis 5 - Welcome.";
+  INFO << "    * SampleAnalyzer for MadAnalysis 5 - Welcome.";
   INFO << endmsg;
  
 }
@@ -74,7 +73,7 @@
   std::string filename = cfg_.GetInputListName();
 
   // Checks if a file has been provided
-  INFO << "      - extracting event samples..." << endmsg;
+  INFO << "      - extracting the list of event samples..." << endmsg;
   std::ifstream input(filename.c_str());
   if (!input)
   {
@@ -151,7 +150,7 @@
 
   // Initialize (common part to all analyses)
   if (!myAnalysis->PreInitialize(outputname,
-                                 cfg_.IsNoEventWeight()))
+                                 &cfg_))
   {
     ERROR << "problem during the pre-initialization of the analysis called '" 
           << name << "'" << endmsg;
@@ -195,7 +194,7 @@
 
   // Initialize (common part to all filters)
   if (!myFilter->PreInitialize(outputname,
-                               cfg_.IsNoEventWeight()))
+                               &cfg_))
   {
     ERROR << "problem during the pre-initialization of the filter called '" 
           << name << "'" << endmsg;
@@ -237,7 +236,7 @@
   writers_.push_back(myWriter);
 
   // Initializing
-  if (!myWriter->Initialize(outputname))
+  if (!myWriter->Initialize(&cfg_,outputname))
   {
     ERROR << "problem during the initialization of the writer called '" 
           << name << "'" << endmsg;
@@ -302,6 +301,10 @@
   if (progressBar_!=0)
   {
     progressBar_->Finalize();
+    INFO << "        => total number of events: " << counter_read_[file_index_-1] 
+         << " ( analyzed: " << counter_passed_[file_index_-1] 
+         << " ; skipped: " << counter_read_[file_index_-1] - counter_passed_[file_index_-1]
+         << " ) " << endmsg;
   }
 
   // Next file
 
=== modified file 'tools/SampleAnalyzer/Core/xdr_istream.cpp'
--- tools/SampleAnalyzer/Core/xdr_istream.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Core/xdr_istream.cpp	2013-06-13 12:22:26 +0000
@@ -115,7 +115,7 @@
 {
   if (eof()) return (*this);
 
-	UInt_t n;
+	UInt_t n=0;
 	(*this)>>n;
 	Float_t* vp = reinterpret_cast<Float_t*>(&n);
 	v=*vp;
@@ -130,7 +130,7 @@
 {
   if (eof()) return (*this);
 
-	ULong64_t n;
+	ULong64_t n=0;
 	(*this)>>n;
 	Double_t* vp = reinterpret_cast<Double_t*>(&n);
 	v=*vp;
 
=== modified file 'tools/SampleAnalyzer/DataFormat/MCParticleFormat.h'
--- tools/SampleAnalyzer/DataFormat/MCParticleFormat.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/DataFormat/MCParticleFormat.h	2013-06-13 12:22:26 +0000
@@ -70,7 +70,8 @@
                               /// 2 intermediate state, 1 final state)
   UInt_t 	        mothup1_;   /// first mother index
   UInt_t 	        mothup2_;   /// second mother index
-  Int_t           extra_;
+  Int_t           extra1_;
+  Int_t           extra2_;
 
   std::vector<MCParticleFormat*> Daughters_;
 
@@ -96,7 +97,8 @@
   {
     momentum_.SetPxPyPzE(0.,0.,0.,0.);
     ctau_=0.; spin_=0.; pdgid_=0; 
-    statuscode_=0; mothup1_=0; mothup2_=0; mother1_=0; mother2_=0; extra_=0;
+    statuscode_=0; mothup1_=0; mothup2_=0; mother1_=0; mother2_=0; 
+    extra1_=0; extra2_=0;
   }
 
   /// Print particle informations
 
=== modified file 'tools/SampleAnalyzer/DataFormat/RecEventFormat.h'
--- tools/SampleAnalyzer/DataFormat/RecEventFormat.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/DataFormat/RecEventFormat.h	2013-06-13 12:22:26 +0000
@@ -48,6 +48,7 @@
 class ROOTReader;
 class TauTagger;
 class bTagger;
+class JetClusteringFastJet;
 
 class RecEventFormat
 {
@@ -56,6 +57,7 @@
   friend class ROOTReader;
   friend class TauTagger;
   friend class bTagger;
+  friend class JetClusteringFastJet;
 
   // -------------------------------------------------------------
   //                        data members
@@ -90,19 +92,19 @@
   Double_t THT_;
 
   /// Monte Carlo taus decaying hadronically
-  std::vector<MCParticleFormat*> MCHadronicTaus_;
+  std::vector<const MCParticleFormat*> MCHadronicTaus_;
 
   /// Monte Carlo taus decaying into muon
-  std::vector<MCParticleFormat*> MCMuonicTaus_;
+  std::vector<const MCParticleFormat*> MCMuonicTaus_;
 
   /// Monte Carlo taus decaying into electron
-  std::vector<MCParticleFormat*> MCElectronicTaus_;
+  std::vector<const MCParticleFormat*> MCElectronicTaus_;
 
   /// Monte Carlo b-quarks
-  std::vector<MCParticleFormat*> MCBquarks_;
+  std::vector<const MCParticleFormat*> MCBquarks_;
 
   /// Monte Carlo c-quarks
-  std::vector<MCParticleFormat*> MCCquarks_;
+  std::vector<const MCParticleFormat*> MCCquarks_;
 
 
   // -------------------------------------------------------------
@@ -146,23 +148,23 @@
   const Double_t& THT() const {return THT_;}
 
   /// Accessor to the Monte Carlo taus decaying hadronically
-  const std::vector<MCParticleFormat*>& MCHadronicTaus() const
+  const std::vector<const MCParticleFormat*>& MCHadronicTaus() const
   {return MCHadronicTaus_;}
 
   /// Accessor to Monte Carlo taus decaying into muon
-  const std::vector<MCParticleFormat*>& MCMuonicTaus() const
+  const std::vector<const MCParticleFormat*>& MCMuonicTaus() const
   {return MCMuonicTaus_;}
 
   /// Accessor to Monte Carlo taus decaying into electron
-  const std::vector<MCParticleFormat*>& MCElectronicTaus() const
+  const std::vector<const MCParticleFormat*>& MCElectronicTaus() const
   {return MCElectronicTaus_;}
 
   /// Accessor to Monte Carlo b-quarks
-  const std::vector<MCParticleFormat*>& MCBquarks() const
+  const std::vector<const MCParticleFormat*>& MCBquarks() const
   {return MCBquarks_;}
 
   /// Accessor to Monte Carlo c-quarks
-  const std::vector<MCParticleFormat*>& MCCquarks() const
+  const std::vector<const MCParticleFormat*>& MCCquarks() const
   {return MCCquarks_;}
 
   /// Accessor to the electron collection
@@ -193,23 +195,23 @@
   Double_t& THT() {return THT_;}
 
   /// Accessor to the Monte Carlo taus decaying hadronically
-  std::vector<MCParticleFormat*>& MCHadronicTaus()
+  std::vector<const MCParticleFormat*>& MCHadronicTaus()
   {return MCHadronicTaus_;}
 
   /// Accessor to Monte Carlo taus decaying into muon
-  std::vector<MCParticleFormat*>& MCMuonicTaus()
+  std::vector<const MCParticleFormat*>& MCMuonicTaus()
   {return MCMuonicTaus_;}
 
   /// Accessor to Monte Carlo taus decaying into electron
-  std::vector<MCParticleFormat*>& MCElectronicTaus()
+  std::vector<const MCParticleFormat*>& MCElectronicTaus()
   {return MCElectronicTaus_;}
 
   /// Accessor to Monte Carlo b-quarks
-  std::vector<MCParticleFormat*>& MCBquarks()
+  std::vector<const MCParticleFormat*>& MCBquarks()
   {return MCBquarks_;}
 
   /// Accessor to Monte Carlo c-quarks
-  std::vector<MCParticleFormat*>& MCCquarks()
+  std::vector<const MCParticleFormat*>& MCCquarks()
   {return MCCquarks_;}
 
   /// Clearing all information
 
=== modified file 'tools/SampleAnalyzer/DataFormat/RecParticleFormat.h'
--- tools/SampleAnalyzer/DataFormat/RecParticleFormat.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/DataFormat/RecParticleFormat.h	2013-06-13 12:22:26 +0000
@@ -53,7 +53,7 @@
  protected:
    
   Float_t 	        HEoverEE_; /// hadronic energy over electromagnetic energy
-  MCParticleFormat* mc_ ;      /// mother generated particle
+  const MCParticleFormat* mc_ ;      /// mother generated particle
 
   // -------------------------------------------------------------
   //                      method members
@@ -94,7 +94,7 @@
   const MCParticleFormat* mc() const {return mc_;}
 
   /// Mutator relatied to matched Monte Carlo particle
-  void setMc(MCParticleFormat* mc) {mc_=mc;}
+  void setMc(const MCParticleFormat* mc) {mc_=mc;}
 
   /// Accessor to hadronic energy / electromagnetic energy ratio
   const Float_t& HEoverEE() const {return HEoverEE_;}
 
=== modified file 'tools/SampleAnalyzer/DataFormat/RecTauFormat.h'
--- tools/SampleAnalyzer/DataFormat/RecTauFormat.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/DataFormat/RecTauFormat.h	2013-06-13 12:22:26 +0000
@@ -98,13 +98,25 @@
   const Float_t  charge() const  
   { if (charge_) return +1.; else return -1.; }
 
+  /// Mutator to the electric charge
+  void setCharge(Float_t charge )
+  { if (charge>0) charge_=true; else charge_=false; }
+
   /// Accessor to the number of tracks
   const UShort_t ntracks() const 
   { return ntracks_; }
 
+  /// Mutator to the number of tracks
+  void setNtracks(UShort_t ntracks)
+  { ntracks_=ntracks; }
+
   /// Accessor to the decay mode
   const Int_t DecayMode() const
   { return DecayMode_; }
+
+  /// Mutator to the decay mode
+  void setDecayMode(Int_t mode)
+  { DecayMode_=mode; }
 };
 
 }
 
=== modified file 'tools/SampleAnalyzer/Filter/FilterBase.h'
--- tools/SampleAnalyzer/Filter/FilterBase.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Filter/FilterBase.h	2013-06-13 12:22:26 +0000
@@ -93,10 +93,10 @@
 
   /// Initialize (common part to all analyses)
   bool PreInitialize(const std::string& outputName, 
-                     bool weighted_events=false)
+                     const Configuration* cfg)
   { 
     outputName_ = outputName;
-    weighted_events_ = weighted_events;
+    weighted_events_ = cfg->IsNoEventWeight();
     return true;
   }
 
 
=== modified file 'tools/SampleAnalyzer/JetClustering/JetClustererBase.h'
--- tools/SampleAnalyzer/JetClustering/JetClustererBase.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/JetClustererBase.h	2013-06-13 12:22:26 +0000
@@ -55,7 +55,7 @@
     /// Is the jet clustering exclusive ?
     Bool_t Exclusive_;
 
-    /// Exclusive id for tau-elec-jet
+    /// Exclusive id for tau-elec-photon-jet
     Bool_t ExclusiveId_;
 
     /// Tagger
@@ -63,6 +63,12 @@
     cTagger*    myCtagger_;
     TauTagger*  myTAUtagger_;
 
+    UInt_t muon;
+    UInt_t electron;
+    UInt_t tauH;
+    UInt_t tauM;
+    UInt_t tauE;
+    UInt_t photon;
 
     //--------------------------------------------------------------------------
     //                              method members
@@ -81,6 +87,12 @@
       myBtagger_   = 0;
       myCtagger_   = 0;
       myTAUtagger_ = 0;
+     muon=0;
+     electron=0;
+     tauH=0;
+     tauM=0;
+     tauE=0;
+     photon=0;
     }
 
     /// Destructor
 
=== modified file 'tools/SampleAnalyzer/JetClustering/JetClusteringFastJet.cpp'
--- tools/SampleAnalyzer/JetClustering/JetClusteringFastJet.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/JetClusteringFastJet.cpp	2013-06-13 12:22:26 +0000
@@ -27,13 +27,230 @@
 
 using namespace MA5;
 
+
+void JetClusteringFastJet::GetFinalState(const MCParticleFormat* part, std::set<const MCParticleFormat*>& finalstates)
+{
+  for (unsigned int i=0; i<part->Daughters().size(); i++)
+  {
+    if (PHYSICS->IsFinalState(part->Daughters()[i])) finalstates.insert(part->Daughters()[i]);
+    else return GetFinalState(part->Daughters()[i],finalstates);
+  }
+}
+
+Bool_t JetClusteringFastJet::IsLast(const MCParticleFormat* part, EventFormat& myEvent)
+{
+  for (unsigned int i=0; i<part->Daughters().size(); i++)
+  {
+    if (part->Daughters()[i]->pdgid()==part->pdgid()) return false;
+  }
+  return true;
+}
+
+
+Bool_t JetClusteringFastJet::IrrelevantPhoton(const MCParticleFormat* part)
+{
+  if (part->mother1()==0) return false;
+
+  // Checking mother
+  UInt_t absid = std::abs(part->mother1()->pdgid());
+
+  if (absid==15) return true;
+  else return IrrelevantPhoton(part->mother1());
+}
+
+Bool_t JetClusteringFastJet::ComingFromHadronDecay(const MCParticleFormat* part)
+{
+  // Weird case ? Safety: removing this case
+  if (part->mother1()==0) return true;
+  //  std::cout << "part " << part->pdgid() << " <- " << part->mother1()->pdgid() << std::endl;
+  // Checking mother
+  Bool_t had = PHYSICS->IsHadronic(part->mother1()->pdgid()) && part->mother1()->pdgid()!=21;
+  if (had && part->mother1()->mother1()==0) return false;
+  else if (had) return true;
+  else return ComingFromHadronDecay(part->mother1());
+}
+
+
 bool JetClusteringFastJet::Execute(SampleFormat& mySample, EventFormat& myEvent)
 {
-  if (mySample.mc()==0 ||  myEvent.mc()==0) return false;
+  if (mySample.mc()==0 || myEvent.mc()==0) return false;
   if (mySample.rec()==0) mySample.InitializeRec();
   if (myEvent.rec() ==0) myEvent.InitializeRec();
+
+  // Reseting the reconstructed event
   myEvent.rec()->Reset();
 
+  // Veto
+  std::vector<bool> vetos(myEvent.mc()->particles().size(),false);
+  std::set<const MCParticleFormat*> vetos2;
+
+  // Filling the dataformat with electron/muon
+  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
+  {
+    const MCParticleFormat& part = myEvent.mc()->particles()[i];
+    UInt_t absid = std::abs(part.pdgid());
+
+    // Rejecting particle with a null pt (initial state ?)
+    if (part.pt()<1e-10) continue;
+
+    // Treating intermediate particles
+    if (PHYSICS->IsInterState(part))
+    {
+      // rejecting not interesting particles
+      if (absid!=5 && absid!=4 && absid!=15) continue;
+
+      // keeping the last particle with the same id in the decay chain
+      if (!IsLast(&part, myEvent)) continue;
+
+      // looking for b quarks
+      if (absid==5)
+      {
+        bool found=false;
+        for (unsigned int j=0;j<myEvent.rec()->MCBquarks_.size();j++)
+        {
+          if (myEvent.rec()->MCBquarks_[j]==&(part)) 
+          {found=true; break;}
+        }
+        if (!found) myEvent.rec()->MCBquarks_.push_back(&(part));
+      }
+
+      // looking for c quarks
+      else if (absid==4)
+      {
+        bool found=false;
+        for (unsigned int j=0;j<myEvent.rec()->MCCquarks_.size();j++)
+        {
+          if (myEvent.rec()->MCCquarks_[j]==&(part)) 
+          {found=true; break;}
+        }
+        if (!found) myEvent.rec()->MCCquarks_.push_back(&(part));
+      }
+
+      // looking for taus
+      else if (absid==15)
+      {
+        // rejecting particle if coming from hadronization
+        if (ComingFromHadronDecay(&part)) continue;
+
+        // Looking taus daughters id
+        bool leptonic   = true;
+        bool muonic     = false;
+        bool electronic = false;
+        for (unsigned int j=0;j<part.Daughters().size();j++)
+        {
+          UInt_t pdgid = std::abs(part.Daughters()[j]->pdgid());
+          if      (pdgid==13) muonic=true;
+          else if (pdgid==11) electronic=true;
+          else if (pdgid!=22 /*photons*/ &&
+                   !(pdgid>=11 && pdgid<=16) /*neutrinos*/) 
+                  leptonic=false;
+        }
+        if (!leptonic) {muonic=false; electronic=false;}
+
+        // Saving taus decaying into muons (only one copy)
+        if (muonic)
+        {
+          bool found=false;
+          for (unsigned int j=0;j<myEvent.rec()->MCMuonicTaus_.size();j++)
+          {
+            if (myEvent.rec()->MCMuonicTaus_[j]==&(part)) 
+            {found=true; break;}
+          }
+          if (!found) myEvent.rec()->MCMuonicTaus_.push_back(&(part));
+        }
+
+        // Saving taus decaying into electrons (only one copy)
+        else if (electronic)
+        {
+          bool found=false;
+          for (unsigned int j=0;j<myEvent.rec()->MCElectronicTaus_.size();j++)
+          {
+            if (myEvent.rec()->MCElectronicTaus_[j]==&(part)) 
+            {found=true; break;}
+          }
+          if (!found) myEvent.rec()->MCElectronicTaus_.push_back(&(part));
+        }
+
+        // Saving taus decaying into hadrons (only copy)
+        else
+        {
+          bool found=false;
+          for (unsigned int j=0;j<myEvent.rec()->MCHadronicTaus_.size();j++)
+          {
+            if (myEvent.rec()->MCHadronicTaus_[j]==&(part)) 
+            {found=true; break;}
+          }
+          if (!found) 
+          {
+            // Saving the hadrons in MC container
+            myEvent.rec()->MCHadronicTaus_.push_back(&(part));
+
+            // Applying efficiency
+            if (!myTAUtagger_->IsIdentified()) continue;
+
+            // Creating reco hadronic taus
+            RecTauFormat* myTau = myEvent.rec()->GetNewTau();
+            if (part.pdgid()>0) myTau->setCharge(-1.);
+            else myTau->setCharge(+1.);
+            myTau->setMomentum(part.momentum());
+            myTau->setMc(&part);
+            myTau->setDecayMode(PHYSICS->GetTauDecayMode(myTau->mc()));
+            if (myTau->DecayMode()<=0) myTau->setNtracks(0); // ERROR case
+            else if (myTau->DecayMode()==7 || 
+                     myTau->DecayMode()==9) myTau->setNtracks(3); // 3-Prong
+            else myTau->setNtracks(1); // 1-Prong
+
+            // Searching final state
+            GetFinalState(&part,vetos2);
+          }
+        }
+      }
+    }
+
+    // Keeping only final states
+    else if (PHYSICS->IsFinalState(part))
+    {
+      // rejecting not interesting particles
+      if (absid!=22 && absid!=11 && absid!=13) continue;
+
+      // rejecting particle if coming from hadronization
+      if (ExclusiveId_ && ComingFromHadronDecay(&part)) continue;
+
+      // Muons
+      if (absid==13)
+      {
+        vetos[i]=true;
+        RecLeptonFormat * muon = myEvent.rec()->GetNewMuon();
+        muon->setMomentum(part.momentum());
+        muon->setMc(&(part));
+        if (part.pdgid()==13) muon->SetCharge(-1);
+        else muon->SetCharge(+1);
+      }
+
+      // Electrons
+      else if (absid==11)
+      {
+        vetos[i]=true;
+        RecLeptonFormat * elec = myEvent.rec()->GetNewElectron();
+        elec->setMomentum(part.momentum());
+        elec->setMc(&(part));
+        if (part.pdgid()==11) elec->SetCharge(-1);
+        else elec->SetCharge(+1);
+      }
+
+      // Photons
+      else if (absid==22)
+      {
+        if (IrrelevantPhoton(&part)) continue;
+        vetos[i]=true;
+        RecPhotonFormat * photon = myEvent.rec()->GetNewPhoton();
+        photon->setMomentum(part.momentum());
+        photon->setMc(&(part));
+      }
+    }
+
+  }
+
   double & TET = myEvent.rec()->TET();
   double & THT = myEvent.rec()->THT();
 
@@ -41,10 +258,21 @@
   std::vector<fastjet::PseudoJet> inputs;
   for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
   {
+    const MCParticleFormat& part = myEvent.mc()->particles()[i];
+
     // Selecting input for jet clustering
     if (myEvent.mc()->particles()[i].statuscode()!=1)       continue;
     if (PHYSICS->IsInvisible(myEvent.mc()->particles()[i])) continue;
-    if (fabs(myEvent.mc()->particles()[i].pdgid())==13)     continue;
+
+    // ExclusiveId mode
+    if (ExclusiveId_)
+    {
+      if (vetos[i]) continue;
+      if (vetos2.find(&part)!=vetos2.end()) continue;
+    }
+
+    // NonExclusive Id mode
+    else if (std::abs(myEvent.mc()->particles()[i].pdgid())==13) continue;
 
     // Filling good particle for clustering
     inputs.push_back(
@@ -90,40 +318,35 @@
     for (unsigned int j=0;j<constituents.size();j++)
     {
       jet->AddConstituent(constituents[j].user_index());
-      //      if (fabs(myEvent.mc()->particles()[constituents[j].user_index()].pdgid())==11) continue;
+      //      if (std::abs(myEvent.mc()->particles()[constituents[j].user_index()].pdgid())==11) continue;
       if (PDG->IsCharged(myEvent.mc()->particles()[constituents[j].user_index()].pdgid())) tracks++;
     }
     jet->ntracks_ = tracks;
   }
 
-  // Filling the dataformat with electron/muon
-  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
+
+  if (ExclusiveId_)
   {
-    if (myEvent.mc()->particles()[i].statuscode()!=1) continue;
-    if (fabs(myEvent.mc()->particles()[i].pdgid())==13)
-    {
-      RecLeptonFormat * muon = myEvent.rec()->GetNewMuon();
-      muon->setMomentum(myEvent.mc()->particles()[i].momentum());
-      muon->setMc(&(myEvent.mc()->particles()[i]));
-      if (myEvent.mc()->particles()[i].pdgid()==13) muon->SetCharge(-1);
-      else muon->SetCharge(+1);
-    }
-    else if (fabs(myEvent.mc()->particles()[i].pdgid())==11)
-    {
-      RecLeptonFormat * elec = myEvent.rec()->GetNewElectron();
-      elec->setMomentum(myEvent.mc()->particles()[i].momentum());
-      elec->setMc(&(myEvent.mc()->particles()[i]));
-      if (myEvent.mc()->particles()[i].pdgid()==11) elec->SetCharge(-1);
-      else elec->SetCharge(+1);
+    for (unsigned int i=0;i<myEvent.rec()->electrons().size();i++)
+    {
+      (*MET) -= myEvent.rec()->electrons()[i].momentum();
+      TET += myEvent.rec()->electrons()[i].pt();
+    }
+    for (unsigned int i=0;i<myEvent.rec()->photons().size();i++)
+    {
+      (*MET) -= myEvent.rec()->photons()[i].momentum();
+      TET += myEvent.rec()->photons()[i].pt();
+    }
+    for (unsigned int i=0;i<myEvent.rec()->taus().size();i++)
+    {
+      (*MET) -= myEvent.rec()->taus()[i].momentum();
+      TET += myEvent.rec()->taus()[i].pt();
     }
   }
 
-
   for (unsigned int i=0;i<myEvent.rec()->muons().size();i++)
   {
     (*MET) -= myEvent.rec()->muons()[i].momentum();
-    //    MET->momentum().SetPx(MET->momentum().Px() - myEvent.rec()->muons()[i].px());
-    //    MET->momentum().SetPy(MET->momentum().Py() - myEvent.rec()->muons()[i].py());
     TET += myEvent.rec()->muons()[i].pt();
   }
 
 
=== modified file 'tools/SampleAnalyzer/JetClustering/JetClusteringFastJet.h'
--- tools/SampleAnalyzer/JetClustering/JetClusteringFastJet.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/JetClusteringFastJet.h	2013-06-13 12:22:26 +0000
@@ -55,6 +55,7 @@
     /// Jet definition
     fastjet::JetDefinition JetDefinition_;
 
+
 //---------------------------------------------------------------------------------
 //                                method members
 //---------------------------------------------------------------------------------
@@ -74,6 +75,13 @@
     /// Initialization
     virtual bool Initialize(const std::map<std::string,std::string>& options)=0;
 
+ private:
+    Bool_t IsLast(const MCParticleFormat* part, EventFormat& myEvent);
+    Bool_t ComingFromHadronDecay(const MCParticleFormat* part);
+    Bool_t IrrelevantPhoton(const MCParticleFormat* part);
+    void GetFinalState(const MCParticleFormat* part, std::set<const MCParticleFormat*>& finalstates);
+ 
+
 };
 
 }
 
=== modified file 'tools/SampleAnalyzer/JetClustering/TaggerBase.cpp'
--- tools/SampleAnalyzer/JetClustering/TaggerBase.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/TaggerBase.cpp	2013-06-13 12:22:26 +0000
@@ -35,6 +35,8 @@
 }
 
 
+
+
 void TaggerBase::SetParameter(const std::string& key, 
                               const std::string& value,
                               std::string header)
 
=== modified file 'tools/SampleAnalyzer/JetClustering/TaggerBase.h'
--- tools/SampleAnalyzer/JetClustering/TaggerBase.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/TaggerBase.h	2013-06-13 12:22:26 +0000
@@ -79,7 +79,7 @@
     }
 
     /// Destructor
-    ~TaggerBase()
+    virtual ~TaggerBase()
     {}
 
     /// Matching using dr
 
=== modified file 'tools/SampleAnalyzer/JetClustering/TauTagger.cpp'
--- tools/SampleAnalyzer/JetClustering/TauTagger.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/TauTagger.cpp	2013-06-13 12:22:26 +0000
@@ -28,116 +28,10 @@
 
 void TauTagger::Method1 (SampleFormat& mySample, EventFormat& myEvent)
 {
-  // loop on the particles searching for tau
-  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
-  {
-    // Keeping only taus
-    if (fabs(myEvent.mc()->particles()[i].pdgid())!=15) continue;
-
-    // Removing initial states
-    if (PHYSICS->IsInitialState(myEvent.mc()->particles()[i])) continue;
-
-    // Removing final states
-    if (PHYSICS->IsFinalState(myEvent.mc()->particles()[i])) continue;
-
-    // Keeping the last taus in the decay chain
-    if (!IsLast(&myEvent.mc()->particles()[i], myEvent)) continue;
-
-    // Removing taus decaying to leptons
-    bool leptonic=true;
-    bool muonic=false;
-    bool electronic=false;
-    for (unsigned int j=0;j<myEvent.mc()->particles()[i].Daughters().size();j++)
-    {
-      unsigned int pdgid = 
-        std::abs(myEvent.mc()->particles()[i].Daughters()[j]->pdgid());
-      if (pdgid==13) muonic=true;
-      else if (pdgid==11) electronic=true;
-      if (pdgid!=22 && !(pdgid>=11 && pdgid<=18)) {leptonic=false;}
-    }
-    if (leptonic && muonic)
-    {
-      bool found=false;
-      for (unsigned int j=0;j<myEvent.rec()->MCMuonicTaus_.size();j++)
-      {
-        if (myEvent.rec()->MCMuonicTaus_[j]==&(myEvent.mc()->particles()[i])) 
-        {found=true; break;}
-      }
-      if (!found) 
-        myEvent.rec()->MCMuonicTaus_.push_back(&(myEvent.mc()->particles()[i]));
-    }
-    else if (leptonic && electronic)
-    {
-      bool found=false;
-      for (unsigned int j=0;j<myEvent.rec()->MCElectronicTaus_.size();j++)
-      {
-        if (myEvent.rec()->MCElectronicTaus_[j]==&(myEvent.mc()->particles()[i])) 
-        {found=true; break;}
-      }
-      if (!found) 
-        myEvent.rec()->MCElectronicTaus_.push_back(&(myEvent.mc()->particles()[i]));
-    }
-    else if (!leptonic)
-    {
-      bool found=false;
-      for (unsigned int j=0;j<myEvent.rec()->MCHadronicTaus_.size();j++)
-      {
-        if (myEvent.rec()->MCHadronicTaus_[j]==&(myEvent.mc()->particles()[i])) 
-        {found=true; break;}
-      }
-      if (!found) 
-        myEvent.rec()->MCHadronicTaus_.push_back(&(myEvent.mc()->particles()[i]));
-
-    }
-  }
-
-  // Matching MCtaus and RECtaus
-  std::vector<std::pair<RecJetFormat*,MCParticleFormat*> > Candidates;
-  for (unsigned int i=0;i<myEvent.rec()->MCHadronicTaus_.size();i++)
-  {
-    RecJetFormat* tag = 0;
-    Double_t DeltaRmax = DeltaRmax_;
-
-    // loop on the jets
-    for (unsigned int j=0;j<myEvent.rec()->jets().size();j++)
-    {
-      //      if (myEvent.rec()->jets()[j].ntracks()!=1 && 
-      //    myEvent.rec()->jets()[j].ntracks()!=3) continue;
-
-      // Calculating Delta R
-      Float_t DeltaR = 
-         myEvent.rec()->MCHadronicTaus_[i]->dr(myEvent.rec()->jets()[j]);
-
-      if (DeltaR <= DeltaRmax)
-      {
-        if(Exclusive_)
-        {
-          tag = &(myEvent.rec()->jets()[j]);
-          DeltaRmax = DeltaR;
-        }
-        else Candidates.push_back( 
-             std::make_pair(& myEvent.rec()->jets()[j],
-                            myEvent.rec()->MCHadronicTaus_[i]) );
-      }
-    }
-    if (Exclusive_ && tag!=0) 
-         Candidates.push_back(std::make_pair(tag,
-                                             myEvent.rec()->MCHadronicTaus_[i]));
-  }
-
-  // Performing efficiency
-  for (unsigned int i=0;i<Candidates.size();i++)
-  {
-    if (!IsIdentified()) continue;
-    Candidates[i].first->mc_ = Candidates[i].second;
-    RecTauFormat* myTau = myEvent.rec()->GetNewTau();
-    Jet2Tau(Candidates[i].first, myTau, myEvent);
-    //      myEvent.rec()->jets().erase((std::vector<RecJetFormat>::iterator) Candidates[j-1]);
-  }
-
   // Performing mis-id
   if (doMisefficiency_)
   {
+    std::vector<unsigned int> toRemove;
     for (unsigned int i=0;i<myEvent.rec()->jets().size();i++)
     {
       // keeping only light jets
@@ -149,8 +43,14 @@
       {
         RecTauFormat* myTau = myEvent.rec()->GetNewTau();
         Jet2Tau(&myEvent.rec()->jets()[i], myTau, myEvent);
+        toRemove.push_back(i);
       }
     }
+    for (unsigned int i =0;i<toRemove.size();i++)
+    {
+      myEvent.rec()->jets().erase(myEvent.rec()->jets().begin()
+                                  + toRemove[i]);
+    }
   }
 
 }
@@ -338,8 +238,8 @@
     charge += PDG->GetCharge(myEvent.mc()->particles()[myJet->Constituents_[i]].pdgid());
   }
 
-  if (charge == 3) myTau->charge_ = true;
-  else if (charge == -3) myTau->charge_ = false;
+  if (charge>0) myTau->charge_ = true;
+  else if (charge<=0) myTau->charge_ = false;
 }
 
 
@@ -378,8 +278,6 @@
 std::string TauTagger::GetParameters()
 {
   std::stringstream str;
-  str << "dR=" << DeltaRmax_ << " ; ";
-  if (Exclusive_) str << "Exclusive ; "; else str << "Inclusive ; ";
   str << "IDeff=" << Efficiency_;
   str << " ; MisID(q)=" << misid_ljet_;
   return str.str();
 
=== modified file 'tools/SampleAnalyzer/JetClustering/TauTagger.h'
--- tools/SampleAnalyzer/JetClustering/TauTagger.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/TauTagger.h	2013-06-13 12:22:26 +0000
@@ -27,7 +27,7 @@
 
 //SampleAnalyser headers
 #include "SampleAnalyzer/JetClustering/TaggerBase.h"
-
+#include <TH1F.h>
 
 namespace MA5
 {
@@ -44,7 +44,7 @@
 
   /// Apply misefficiency
   Bool_t doMisefficiency_;
-
+  //  TH1F* htest;
 
 //---------------------------------------------------------------------------------
 //                                method members
@@ -56,10 +56,12 @@
     {
       misid_ljet_=0.0;
       doMisefficiency_=false;
+      //      htest = new TH1F("eric","eric",10,0,10);
     }
 
     /// Destructor
-    ~TauTagger () { }
+    virtual ~TauTagger () 
+    { /*htest->SaveAs("Alloul.root");*/ }
 
     /// Matching using dr
     virtual void Method1 (SampleFormat& mySample, EventFormat& myEvent);
 
=== modified file 'tools/SampleAnalyzer/JetClustering/bTagger.cpp'
--- tools/SampleAnalyzer/JetClustering/bTagger.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/bTagger.cpp	2013-06-13 12:22:26 +0000
@@ -29,49 +29,6 @@
 {
   std::vector<RecJetFormat*> Candidates;
 
-  // loop on the particles searching for b
-  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
-  {
-    // Keeping only b-quark
-    if (std::abs(myEvent.mc()->particles()[i].pdgid())!=5 &&
-        std::abs(myEvent.mc()->particles()[i].pdgid())!=4) continue;
-
-    // Removing initial states
-    if (PHYSICS->IsInitialState(myEvent.mc()->particles()[i])) continue;
-
-    // Removing final states
-    if (PHYSICS->IsFinalState(myEvent.mc()->particles()[i])) continue;
-
-    // Keeping the last taus in the decay chain
-    if (!IsLast(&myEvent.mc()->particles()[i], myEvent)) continue;
-
-    // b-quark already defined ?
-    if (std::abs(myEvent.mc()->particles()[i].pdgid())==5)
-    {
-      bool found=false;
-      for (unsigned int j=0;j<myEvent.rec()->MCBquarks_.size();j++)
-      {
-        if (myEvent.rec()->MCBquarks_[j]==&(myEvent.mc()->particles()[i])) 
-        {found=true; break;}
-      }
-      if (!found) 
-        myEvent.rec()->MCBquarks_.push_back(&(myEvent.mc()->particles()[i]));
-    }
-
-    // c-quark already defined ?
-    else if (std::abs(myEvent.mc()->particles()[i].pdgid())==4)
-    {
-      bool found=false;
-      for (unsigned int j=0;j<myEvent.rec()->MCCquarks_.size();j++)
-      {
-        if (myEvent.rec()->MCCquarks_[j]==&(myEvent.mc()->particles()[i])) 
-        {found=true; break;}
-      }
-      if (!found) 
-        myEvent.rec()->MCCquarks_.push_back(&(myEvent.mc()->particles()[i]));
-    }
-  }
-
   // Matching b-quarks to jets
   for (unsigned int i=0;i<myEvent.rec()->MCBquarks_.size();i++)
   {
@@ -81,6 +38,7 @@
     // loop on the jets
     for (unsigned int j=0;j<myEvent.rec()->jets().size();j++)
     {
+      if (myEvent.rec()->jets()[j].pt()<1e-10) continue;
       Float_t DeltaR = myEvent.rec()->MCBquarks_[i]->dr(myEvent.rec()->jets()[j]);
 
       if (DeltaR <= DeltaRmax) 
@@ -112,6 +70,8 @@
     // loop on the jets
     for (unsigned int j=0;j<myEvent.rec()->jets().size();j++)
     {
+      if (myEvent.rec()->jets()[j].pt()<1e-10) continue;
+
       Float_t DeltaR = 
           myEvent.rec()->MCCquarks_[i]->dr(myEvent.rec()->jets()[j]);
 
 
=== modified file 'tools/SampleAnalyzer/JetClustering/bTagger.h'
--- tools/SampleAnalyzer/JetClustering/bTagger.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/bTagger.h	2013-06-13 12:22:26 +0000
@@ -62,7 +62,7 @@
     }
 
     /// Destructor
-    ~bTagger() {}
+    virtual ~bTagger() {}
 
     /// Matching using dr
     virtual void Method1 (SampleFormat& mySample, EventFormat& myEvent);
 
=== modified file 'tools/SampleAnalyzer/JetClustering/cTagger.h'
--- tools/SampleAnalyzer/JetClustering/cTagger.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/JetClustering/cTagger.h	2013-06-13 12:22:26 +0000
@@ -45,7 +45,7 @@
     cTagger() {}
 
     /// Destructor
-    ~cTagger () {}
+    virtual ~cTagger () {}
 
     /// Matching using dr
     virtual void Method1 (SampleFormat& mySample, EventFormat& myEvent);
 
=== removed file 'tools/SampleAnalyzer/Makefile'
--- tools/SampleAnalyzer/Makefile	2013-04-19 17:21:55 +0000
+++ tools/SampleAnalyzer/Makefile	1970-01-01 00:00:00 +0000
@@ -1,89 +0,0 @@
-################################################################################
-#                 MAKEFILE DEVOTED TO SAMPLE ANALYZER LIBRARY                  #
-################################################################################
-
-# Compilators
-GCC = g++
-
-# Options for compilation
-CXXFASTJET = $(shell fastjet-config --cxxflags --plugins)
-CXXFLAGS = -Wall -O3 $(shell root-config --cflags) -I./../ -DZIP_USE -DFASTJET_USE $(CXXFASTJET)
-
-# Files
-SRCS = $(wildcard */*.cpp)
-HDRS = $(wildcard */*.h)
-OBJS = $(SRCS:.cpp=.o)
-
-# Name of the library
-PROGRAM = SampleAnalyzer
-
-# Defining colours
-GREEN  = "\\033[1;32m"
-RED    = "\\033[1;31m"
-PINK   = "\\033[1;35m"
-BLUE   = "\\033[1;34m"
-YELLOW = "\\033[1;33m"
-CYAN   = "\\033[1;36m"
-NORMAL = "\\033[0;39m"
-
-# All target
-all: header compile_header compile link_header link
-
-# Header target
-header:
-	@echo -e $(YELLOW)"--------------------------------------------------"
-	@echo -e "         Building SampleAnalyzer library          "
-	@echo -e "--------------------------------------------------"$(NORMAL)
-
-# Compile_header target
-compile_header:
-	@echo -e $(YELLOW)"--------------------------------------------------"
-	@echo -e "                   Compilation                    "
-	@echo -e "--------------------------------------------------"$(NORMAL)
-
-# Link_header target
-link_header:
-	@echo -e $(YELLOW)"--------------------------------------------------"
-	@echo -e "                  Final linking                   "
-	@echo -e "--------------------------------------------------"$(NORMAL)
-
-# clean_header target
-clean_header:
-	@echo -e $(YELLOW)"--------------------------------------------------"
-	@echo -e "    Removing intermediate files from building     "
-	@echo -e "--------------------------------------------------"$(NORMAL)
-
-# mrproper_header target
-mrproper_header:
-	@echo -e $(YELLOW)"--------------------------------------------------"
-	@echo -e "             Cleaning all the project             "
-	@echo -e "--------------------------------------------------"$(NORMAL)
-
-# Precompile target
-precompile:
-
-# Compile target
-compile: precompile $(OBJS)
-
-# Link target
-link: $(OBJS)
-	ar -ruc Lib/lib$(PROGRAM).a $(OBJS)
-
-# Phony target
-.PHONY: do_clean header compile_header link_header
-
-# Clean target
-clean: clean_header do_clean
-
-# Do clean target
-do_clean: 
-	@rm -f $(OBJS) Reader/FACdict.*
-
-# Mr Proper target 
-mrproper: mrproper_header do_mrproper
-
-# Do Mr Proper target 
-do_mrproper: do_clean
-	@rm -f Lib/lib$(PROGRAM).a
-	@rm -f compilation.log linking.log cleanup.log *~ */*~
-
 
=== modified file 'tools/SampleAnalyzer/Reader/HEPMCReader.cpp'
--- tools/SampleAnalyzer/Reader/HEPMCReader.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Reader/HEPMCReader.cpp	2013-06-13 12:22:26 +0000
@@ -133,6 +133,7 @@
 }
 
 
+
 // -----------------------------------------------------------------------------
 // FinalizeEvent
 // -----------------------------------------------------------------------------
@@ -150,11 +151,39 @@
   mySample.mc()->set_xsection(value);
   mySample.mc()->set_xsection_error(err_value);
 
-  // Mother pointer assignment
+  // Computing met, mht, ...
   for (unsigned int i=0; i<myEvent.mc()->particles_.size();i++)
   {
     MCParticleFormat& part = myEvent.mc()->particles_[i];
 
+    // Setting mother
+    if (part.extra1_!=part.extra2_)
+    {
+      unsigned int nmother=0;
+      for (unsigned int j=0; j < myEvent.mc()->particles_.size();j++)
+      {
+        if (i==j) continue;
+        if (part.extra1_ == myEvent.mc()->particles_[j].extra2_)
+        {
+          // set daughter
+          myEvent.mc()->particles_[j].Daughters_.push_back(&part);
+
+          // set mother
+          nmother++;
+          if      (nmother==1) part.mother1_=&(myEvent.mc()->particles()[j]);
+          else if (nmother==2) part.mother2_=&(myEvent.mc()->particles()[j]);
+          else 
+          { 
+            if (warnmother_) 
+            {
+              WARNING << "Number of mothers greather than 2 : " << nmother << endmsg; 
+              warnmother_=false; 
+            }
+          }
+        }
+      }
+    }
+
     // MET, MHT, TET, THT
     if (part.statuscode()==1 && !PHYSICS->IsInvisible(part))
     {
@@ -410,12 +439,12 @@
 // FillEventParticleLine
 // -----------------------------------------------------------------------------
 void HEPMCReader::FillEventParticleLine(const std::string& line,
-                                      EventFormat& myEvent)
+                                        EventFormat& myEvent)
 {
   std::stringstream str;
   str << line;
 
-  double   	tmp;    // temporary variable to fill in LorentzVector
+  double tmp;    // temporary variable to fill in LorentzVector
 
   // Get a new particle
   MCParticleFormat * part = myEvent.mc()->GetNewParticle();
@@ -432,9 +461,8 @@
   str >> part->statuscode_;
   str >> tmp; 
   str >> tmp; 
-  str >> part->extra_;
-
-  SetMother(part, myEvent);
+  str >> part->extra2_;
+  part->extra1_=current_vertex_.barcode_;
 }
 
 // -----------------------------------------------------------------------------
@@ -460,19 +488,45 @@
 //--------------------------------------------------------------------------
 // SetMother
 //--------------------------------------------------------------------------
-void HEPMCReader::SetMother(MCParticleFormat* part, EventFormat& myEvent)
+void HEPMCReader::SetMother(MCParticleFormat* const part, EventFormat& myEvent)
 {
+  /*  std::cout << current_vertex_.barcode_ << std::endl;
+
+  // No history
   if (myEvent.mc()->particles().size()==0) return;
-  
+
+  //ERIC orphan special treatment 
+  if (part->extra_==current_vertex_.barcode_) return;
+
+  std::cout << "---------------------------------------------" << std::endl;
+  for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
+  {
+    std::cout << "i=" << i << "\t" << myEvent.mc()->particles()[i].pdgid() << " from ";
+    if (myEvent.mc()->particles()[i].mother1()!=0) std::cout <<  myEvent.mc()->particles()[i].mother1()->pdgid();
+    std::cout << " extra=" << myEvent.mc()->particles()[i].extra_ << " current=" << current_vertex_.barcode_ << std::endl;
+    std::cout << std::endl;
+  }
+
+
   unsigned int nmother=0;
-  for (unsigned int i =0; i<(myEvent.mc()->particles().size()-1);i++)
+  for (unsigned int i=0;i<(myEvent.mc()->particles().size()-1);i++)
   {
     if(myEvent.mc()->particles()[i].extra_==current_vertex_.barcode_)
     {
       nmother++;
-      if(nmother==1)      part->mother1_=&(myEvent.mc()->particles()[i]);
-      else if(nmother==2) part->mother2_=&(myEvent.mc()->particles()[i]);
-      else { if (warnmother_) {WARNING << "Number of mothers greather than 2 : " << nmother << endmsg; warnmother_=false;} }
+      if      (nmother==1) part->mother1_=&(myEvent.mc()->particles()[i]);
+      else if (nmother==2) part->mother2_=&(myEvent.mc()->particles()[i]);
+      else 
+      { 
+        if (warnmother_) 
+        {
+          WARNING << "Number of mothers greather than 2 : " << nmother << endmsg; 
+          warnmother_=false; 
+        }
+      }
     }
   }
+
+    if (part->pdgid()==-24) {std::cout << "MMMMHHH" << std::endl; ComingFromHadronDecay(part); }
+  */
 }
 
=== modified file 'tools/SampleAnalyzer/Reader/HEPMCReader.h'
--- tools/SampleAnalyzer/Reader/HEPMCReader.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Reader/HEPMCReader.h	2013-06-13 12:22:26 +0000
@@ -96,7 +96,7 @@
   void FillEventPDFInfo(const std::string& line, SampleFormat& mySample, EventFormat& myEvent);
   void FillEventParticleLine(const std::string& line, EventFormat& myEvent);
   void FillEventVertexLine(const std::string& line, EventFormat& myEvent);
-  void SetMother(MCParticleFormat* part, EventFormat& myEvent);
+  void SetMother(MCParticleFormat* const part, EventFormat& myEvent);
   Bool_t FillWeightNames(const std::string& line);
   Bool_t FillHeavyIons(const std::string& line);
 
 
=== modified file 'tools/SampleAnalyzer/Reader/LHCOReader.cpp'
--- tools/SampleAnalyzer/Reader/LHCOReader.cpp	2013-05-17 07:13:44 +0000
+++ tools/SampleAnalyzer/Reader/LHCOReader.cpp	2013-06-13 12:22:26 +0000
@@ -65,49 +65,49 @@
 
   bool EndOfLoop = false;
   
-// Declarging a new string for line
+  // Declarging a new string for line
   std::string line;
   std::string tmp;
      
   do
-    {
-      std::stringstream str;
-      str.str("");
-      
-      
-      if (!ReadLine(line))
-	{
-	  EndOfFile_=true;
-	  return StatusCode::KEEP;
-	}
-      str << line;
-      str >> tmp;
-
-      if (saved_)
-	{
-	  FillEventInitLine(savedline_,myEvent);
-	  saved_=false;
-	}
-
-      if(tmp=="0")
-	{
-	  if(firstevent_ )
+  {
+    std::stringstream str;
+    str.str("");
+      
+    if (!ReadLine(line))
+    {
+      EndOfFile_=true;
+      return StatusCode::KEEP;
+    }
+    str << line;
+    str >> tmp;
+
+    if (saved_)
+    {
+      FillEventInitLine(savedline_,myEvent);
+      saved_=false;
+    }
+
+    if(tmp=="0")
+    {
+      if(firstevent_ )
 	    {
 	      FillEventInitLine(line,myEvent);
+        firstevent_=false;
+        continue;
 	    }
-	  else 
+      else 
 	    {
 	      EndOfLoop = true;
 	      savedline_=line;
 	      saved_=true;
 	      continue;
 	    }
-	}
+    }
 
-      FillEventParticleLine(line,myEvent);
+    FillEventParticleLine(line,myEvent);
       
-      firstevent_=false;
-    }
+  }
   while(!EndOfLoop);
   
   // Normal end
@@ -280,6 +280,7 @@
     // 9th column : isolation
     str >> tmp;
     muon->sumPT_isol_=std::floor(tmp);
+
     Float_t ET_PT=tmp-muon->sumPT_isol_;
     Bool_t test=false;
     for (unsigned int j=0;j<5;j++)
@@ -291,7 +292,7 @@
         break;
       }
     }
-    if (test) muon->sumET_isol_=std::floor(ET_PT)*muon->sumPT_isol_;
+    if (test) muon->sumET_isol_=std::floor(ET_PT)*muon->pt();
     else muon->sumET_isol_=0;
   }
 
 
=== modified file 'tools/SampleAnalyzer/Reader/ReaderTextBase.cpp'
--- tools/SampleAnalyzer/Reader/ReaderTextBase.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Reader/ReaderTextBase.cpp	2013-06-13 12:22:26 +0000
@@ -233,6 +233,8 @@
 Long64_t ReaderTextBase::GetPosition()
 {
   if (input_==0) return 0;
+#ifdef ZIP_USE
   if (compress_) return gzinput_->tellg();
+#endif
   else  return input_->tellg();
 }
 
=== added file 'tools/SampleAnalyzer/Service/CompilationService.h'
--- tools/SampleAnalyzer/Service/CompilationService.h	1970-01-01 00:00:00 +0000
+++ tools/SampleAnalyzer/Service/CompilationService.h	2013-06-13 12:22:26 +0000
@@ -0,0 +1,30 @@
+////////////////////////////////////////////////////////////////////////////////
+//  
+//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
+//  The MadAnalysis development team, email: <[e-mail]>
+//  
+//  This file is part of MadAnalysis 5.
+//  Official website: <https://launchpad.net/madanalysis5>
+//  
+//  MadAnalysis 5 is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//  
+//  MadAnalysis 5 is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//  GNU General Public License for more details.
+//  
+//  You should have received a copy of the GNU General Public License
+//  along with MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
+//  
+////////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef TIMER_SERVICE_H
+#define TIMER_SERVICE_H
+
+#define UNUSED(expr) do {(void)(expr);} while(0);
+
+#endif
 
=== modified file 'tools/SampleAnalyzer/Service/Physics.h'
--- tools/SampleAnalyzer/Service/Physics.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Service/Physics.h	2013-06-13 12:22:26 +0000
@@ -132,7 +132,8 @@
   /// Is Initial State
   Bool_t IsInitialState(const MCParticleFormat& part) const
   {
-    return (part.statuscode()==initialstate_);
+    //    return (part.statuscode()==initialstate_);
+    return (part.statuscode()==-1 || (part.statuscode()>=11 && part.statuscode()<=19));
   }
 
   /// Is Final State
@@ -144,7 +145,7 @@
   /// Is Inter State
   Bool_t IsInterState(const MCParticleFormat& part) const
   {
-    return (part.statuscode()!=finalstate_ && part.statuscode()!=initialstate_);
+    return (!IsInitialState(part) && !IsFinalState(part));
   }
 
   /// Is Initial State
@@ -225,6 +226,13 @@
     if (found==mcConfig_.hadronic_ids_.end()) return false; else return true;
   }
 
+  /// Is hadronic ?
+  inline bool IsHadronic(Int_t pdgid) const
+  {
+    std::set<Int_t>::iterator found = mcConfig_.hadronic_ids_.find(pdgid);
+    if (found==mcConfig_.hadronic_ids_.end()) return false; else return true;
+  }
+
   /// Is invisible ?
   inline bool IsInvisible(const MCParticleFormat& part) const
   {
 
=== added file 'tools/SampleAnalyzer/Writer/LHCOParticleFormat.cpp'
--- tools/SampleAnalyzer/Writer/LHCOParticleFormat.cpp	1970-01-01 00:00:00 +0000
+++ tools/SampleAnalyzer/Writer/LHCOParticleFormat.cpp	2013-06-13 12:22:26 +0000
@@ -0,0 +1,98 @@
+////////////////////////////////////////////////////////////////////////////////
+//  
+//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
+//  The MadAnalysis development team, email: <[e-mail]>
+//  
+//  This file is part of MadAnalysis 5.
+//  Official website: <https://launchpad.net/madanalysis5>
+//  
+//  MadAnalysis 5 is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//  
+//  MadAnalysis 5 is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//  GNU General Public License for more details.
+//  
+//  You should have received a copy of the GNU General Public License
+//  along with MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
+//  
+////////////////////////////////////////////////////////////////////////////////
+
+
+// STL headers
+#include <sstream>
+#include <iomanip>
+
+// SampleAnalyzer headers
+#include "SampleAnalyzer/Writer/LHCOParticleFormat.h"
+
+using namespace MA5;
+
+// -----------------------------------------------------------------------------
+// Header
+// -----------------------------------------------------------------------------
+const std::string LHCOParticleFormat::header = "  #  typ    eta     phi        pt    jmas    ntrk    btag  had/em   dum1   dum2";
+
+
+// -----------------------------------------------------------------------------
+// Print a particle line
+// -----------------------------------------------------------------------------
+void LHCOParticleFormat::Print(UInt_t num, std::ostream* out)
+{
+  *out << std::setw(3) << std::right << num;
+  *out << "  ";
+  *out << std::setw(2) << std::right << id;
+  *out << "  ";
+  *out << std::fixed;
+  unsigned int pres = out->precision();
+  out->precision(3);
+  double myeta = eta;
+  if (myeta>9) myeta=9.999;
+  else if (myeta<-9) myeta=-9.999;
+  *out << std::setw(6) << std::right << myeta;  // -X.YYY
+  *out << "  ";
+  *out << std::setw(6) << std::right << phi;    // -X.YYY
+  *out << "  ";
+  *out << std::setw(8) << std::right << pt;     // XXXX.YYY 
+  *out << "  ";
+  *out << std::setw(6) << std::right << jmass;  // -X.YYY
+  *out << "  ";
+  double myntrk = ntrk;
+  if (myntrk>9) myntrk=9.999;
+  else if (myntrk<-9) myntrk=-9.999;
+  *out << std::setw(6) << std::right << myntrk;   // -X.YYY
+  *out << "   ";
+  *out << std::setw(5) << std::right << btag;   // X.YYY
+  *out << "   ";
+  out->precision(2);
+  *out << std::setw(5) << std::right << hadem;  // XX.YY
+  *out << "  ";
+  *out << std::setw(5) << std::right << 0.;     // -X.YY
+  *out << "  ";
+  *out << std::setw(5) << std::right << 0.;     // -X.YY
+  out->precision(pres);
+  *out << std::endl;
+}
+
+
+// -----------------------------------------------------------------------------
+// Print the event header
+// -----------------------------------------------------------------------------
+void LHCOParticleFormat::WriteEventHeader(unsigned int numEvent,
+                                          std::ostream* out)
+{
+  // Particle number
+  *out << std::setw(3) << std::right << 0;
+  *out << "  ";
+
+  // Event number
+  *out << std::setw(10) << std::right << numEvent+1;
+  *out << "  ";
+
+  // Trigger word
+  *out << std::setw(6) << std::right << 0;
+  *out << std::endl;
+}
 
=== added file 'tools/SampleAnalyzer/Writer/LHCOParticleFormat.h'
--- tools/SampleAnalyzer/Writer/LHCOParticleFormat.h	1970-01-01 00:00:00 +0000
+++ tools/SampleAnalyzer/Writer/LHCOParticleFormat.h	2013-06-13 12:22:26 +0000
@@ -0,0 +1,64 @@
+////////////////////////////////////////////////////////////////////////////////
+//  
+//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
+//  The MadAnalysis development team, email: <[e-mail]>
+//  
+//  This file is part of MadAnalysis 5.
+//  Official website: <https://launchpad.net/madanalysis5>
+//  
+//  MadAnalysis 5 is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//  
+//  MadAnalysis 5 is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//  GNU General Public License for more details.
+//  
+//  You should have received a copy of the GNU General Public License
+//  along with MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
+//  
+////////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef LHCO_PARTICLE_FORMAT_h
+#define LHCO_PARTICLE_FORMAT_h
+
+// STL headers
+#include <iostream>
+#include <string>
+
+// ROOT headers
+#include <Rtypes.h>
+
+namespace MA5
+{
+
+class LHCOParticleFormat
+{
+
+  // -------------------------------------------------------------
+  //                        data members
+  // -------------------------------------------------------------
+ public:
+
+  UInt_t id;
+  Float_t eta;
+  Float_t phi;
+  Float_t pt;
+  Float_t jmass;
+  Float_t ntrk;
+  Float_t btag;
+  Float_t hadem;
+
+  static const std::string header;
+
+  void Print(UInt_t num, std::ostream* out);
+  static void WriteEventHeader(unsigned int numEvent,std::ostream* out);
+};
+
+
+}
+
+#endif
 
=== modified file 'tools/SampleAnalyzer/Writer/LHCOWriter.cpp'
--- tools/SampleAnalyzer/Writer/LHCOWriter.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/LHCOWriter.cpp	2013-06-13 12:22:26 +0000
@@ -31,105 +31,93 @@
 
 using namespace MA5;
 
-std::string LHCOWriter::FortranFormat_SimplePrecision(Float_t value,UInt_t precision)
-{
-  std::stringstream str;
-  str.precision(precision);
-  std::string word;
-
-  Bool_t negative=false;
-  if (value<0) {negative=true; value*=-1;}
-
-  Int_t exponent = 0;
-  if (value!=0)
-  {
-    for (; value > 1.0; exponent++) value/=10.;
-    for (; value < 0.0; exponent--) value*=10.;
-  }
-
-  str << std::uppercase << std::fixed << value << "E";
-  if (exponent>=0) str << "+"; else str << "-";
-  if (abs(exponent)<10) str << "0";
-  str << abs(exponent);
-  str >> word;
-  if (!negative) return word;
-  else return "-"+word;
-}
-
-
-std::string LHCOWriter::FortranFormat_DoublePrecision(Double_t value,UInt_t precision)
-{
-  std::stringstream str;
-  str.precision(precision);
-  std::string word;
-
-  Bool_t negative=false;
-  if (value<0) {negative=true; value*=-1;}
-
-  Int_t exponent = 0;
-  if (value!=0)
-  {
-    for (; value > 1.0; exponent++) value/=10.;
-    for (; value < 0.0; exponent--) value*=10.;
-  }
-
-  str << std::uppercase << std::fixed << value << "E";
-  if (exponent>=0) str << "+"; else str << "-";
-  if (abs(exponent)<10) str << "0";
-  str << abs(exponent);
-  str >> word;
-  if (!negative) return word;
-  else return "-"+word;
-}
-
 
 /// Read the sample
 bool LHCOWriter::WriteHeader(const SampleFormat& mySample)
 {
-  // Opening tag
-  *output_ << "<LesHouchesEvents version=""1.0"">" << std::endl;
-
-  // Header block
-  *output_ << "<header>" << std::endl;
-  *output_ << "<!--" << std::endl;
-  *output_ << "#*********************************************************************" << std::endl;
-  *output_ << "#                                                                    *" << std::endl;
-  *output_ << "#           This file has been produced by MadAnalysis 5             *" << std::endl;
-  *output_ << "#                                                                    *" << std::endl;
-  *output_ << "#....................................................................*" << std::endl;
-  *output_ << "" << std::endl;
+  // MA5 logo
+  WriteMA5header();
+
+  // LHCO format tag
+  *output_ << "#<MA5Format> LHCO format </MA5Format>" << std::endl;
+
+  // Python interface version
+  *output_ << "#<MadAnalysis5Version> " << cfg_->GetPythonInterfaceVersion() 
+           << " " << cfg_->GetPythonInterfaceDate() 
+           << "</MadAnalysis5Version>" << std::endl;
+
+  // SampleAnalyzer version
+  *output_ << "#<SampleAnalyzerVersion> "<< cfg_->GetSampleAnalyzerVersion()
+           << " " << cfg_->GetSampleAnalyzerVersion()
+           << " </SampleAnalyzerVersion>" << std::endl;
+
+  // Explanation about the LHCO
+  *output_ << "#<FormatDescription>" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
+  *output_ << "# Each event begins with a header row. Meaning of the different columns:        #" << std::endl;
+  *output_ << "# the following:                                                                #" << std::endl;
+  *output_ << "# - column 1: label '0' indicating the beginning a new event.                   #" << std::endl;
+  *output_ << "# - column 2: event number.                                                     #" << std::endl;
+  *output_ << "# - column 3: triggering information. For files produced by MadAnalysis, this   #" << std::endl;
+  *output_ << "#             value is set to 0.                                                #" << std::endl;
+  *output_ << "#                                                                               #" << std::endl;
+  *output_ << "# The rest of the rows are the physics objects in the events. Meaning of the    #" << std::endl;
+  *output_ << "# different columns:                                                            #" << std::endl;
+  *output_ << "# - column 1: row number                                                        #" << std::endl;
+  *output_ << "# - column 2: type of physics object. The possible values are:                  #" << std::endl;
+  *output_ << "#               0 = photon, 1 = electron, 2 = muon,                             #" << std::endl;
+  *output_ << "#               3 = hadronically-decaying tau, 4 = jet, 6 = MET                 #" << std::endl;
+  *output_ << "# - column 3: pseudorapidity (eta) of the physics object                        #" << std::endl;
+  *output_ << "# - column 4: azimuthal angle (phi) of the physics object                       #" << std::endl;
+  *output_ << "# - column 5: transverse momentum (pt) of the physics object                    #" << std::endl;
+  *output_ << "# - column 6: mass of the physics object. Required information to               #" << std::endl;
+  *output_ << "#             reconstruct the four-vector momentum of jets.                     #" << std::endl;
+  *output_ << "# - column 7:                                                                   #" << std::endl;
+  *output_ << "#   + case of a jet: number of tracks associated with the object.               #" << std::endl;
+  *output_ << "#   + case of an electron/muon: electric charge.                                #" << std::endl;
+  *output_ << "#   + case of an hadronically-decaying tau: electric charge times the number of #" << std::endl;
+  *output_ << "#                                           tracks associated with the object.  #" << std::endl;
+  *output_ << "#                                           Usual values: +/-1 , +/-3.          #" << std::endl;
+  *output_ << "#   + other cases: 0.                                                           #" << std::endl;
+  *output_ << "# - column 8:                                                                   #" << std::endl;
+  *output_ << "#   + case of a jet: b-jet tag. Allowed values are:                             #" << std::endl;
+  *output_ << "#                      0 = no identified b-jet,                                 #" << std::endl;
+  *output_ << "#                      1 = b-jet identified by a loose tagger algorithm,        #" << std::endl;
+  *output_ << "#                      2 = b-jet identified by a tight tagger algorithm         #" << std::endl;
+  *output_ << "#                          (the value 2 is not used by MadAnalysis 5).          #" << std::endl;
+  *output_ << "#   + case of a muon: the value indicates the row number of the closest jet     #" << std::endl;
+  *output_ << "#                     in DeltaR.                                                #" << std::endl;
+  *output_ << "#   + other cases: 0.                                                           #" << std::endl;
+  *output_ << "# - column 9:                                                                   #" << std::endl;
+  *output_ << "#   + case of jet/electron/photon: ratio of the hadronic versus                 #" << std::endl;
+  *output_ << "#                                  electromagnetic energy deposited in the      #" << std::endl;
+  *output_ << "#                                  calorimeters cells associated to the object. #" << std::endl;
+  *output_ << "#                                  It is typically > 1 for a jet and < 1 for an #" << std::endl;
+  *output_ << "#                                  electron or a photon.                        #" << std::endl;
+  *output_ << "#   + case of muon: value related to muon isolation. The format of this value   #" << std::endl;
+  *output_ << "#                   is XXX.YYY. The integer part (XXX) is the summed transverse #" << std::endl;
+  *output_ << "#                   momentum of tracks in a DeltaR-cone (typically DeltaR=0.4), #" << std::endl;
+  *output_ << "#                   excluding the muon track. The decimal part (YYY) is the     #" << std::endl;
+  *output_ << "#                   ratio of the transverse energy in the same DeltaR-cone to   #" << std::endl;
+  *output_ << "#                   the transverse momentum of the muon.                        #" << std::endl;
+  *output_ << "#   + in other cases: 0.                                                        #" << std::endl;
+  *output_ << "# - column 10: empty column, added for possible use.                            #" << std::endl;
+  *output_ << "# - column 11: empty column, added for possible use.                            #" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
+  *output_ << "#</FormatDescription>" << std::endl;
 
   if (mySample.mc()!=0)
   {
-    *output_ << "Original header:" << std::endl;
+    *output_ << "#Original header:" << std::endl;
     *output_ << "" << std::endl;
 
     for (unsigned int i=0;i<mySample.mc()->header().size();i++)
-      *output_ << mySample.mc()->header()[i] << std::endl;
+      *output_ << "#" << mySample.mc()->header()[i] << std::endl;
   }
 
+  // One line to remind the meaning of each line
   *output_ << "  #" << std::endl;
-  // *output_ << " ## More info on LHCO files: "
-  //        << "http://cp3wks05.fynu.ucl.ac.be/Manual/lhco.html" << std::endl;
-  *output_ << "  #        eventNum triggerWord" << std::endl;
-  *output_ << "  # typ ";
-  *output_ << std::setw(18) << std::right << "eta";
-  *output_ << " ";
-  *output_ << std::setw(18) << std::right << "phi";
-  *output_ << " ";
-  *output_ << std::setw(18) << std::right << "pt";
-  *output_ << " ";
-  *output_ << std::setw(18) << std::right << "jmass";
-  *output_ << " ";
-  *output_ << std::setw(5) << std::right << "ntracks";
-  *output_ << " ";
-  *output_ << std::setw(4) << std::right << "btag";
-  *output_ << " ";
-  *output_ << std::setw(7) << std::right << "had/em";
-  *output_ << " ";
-  *output_ << std::setw(4) << std::right << "dum1";
-  *output_ << " ";
-  *output_ << std::setw(4) << std::right << "dum2" << std::endl;
+  *output_ << LHCOParticleFormat::header << std::endl;
 
   return true;
 }
@@ -138,67 +126,78 @@
 bool LHCOWriter::WriteEvent(const EventFormat& myEvent, 
                            const SampleFormat& mySample)
 {
-  // Event header
-  *output_ << "0" << std::endl;
+  // FirstEvent
+  if (counter_==0)
+  {
+    WriteHeader(mySample);
+  }
 
   // Skipping event if no reconstructed object
   if (myEvent.rec()==0) return true;
 
-  // Counter of particles
-  unsigned int npart=0;
-
-  // Writing electrons
+  // Event header
+  LHCOParticleFormat::WriteEventHeader(counter_,output_);
+
+  // Particle container in LHCO format
+  std::vector<LHCOParticleFormat> PartTable;
+  PartTable.reserve( myEvent.rec()->photons().size() +
+                     myEvent.rec()->electrons().size() +
+                     myEvent.rec()->muons().size() + 
+                     myEvent.rec()->taus().size() + 
+                     myEvent.rec()->jets().size() + 1 /*MET*/);
+
+  // Writing photons (=0)
+  for (unsigned int i=0;i<myEvent.rec()->photons().size();i++)
+  {
+    PartTable.push_back(LHCOParticleFormat());
+    WritePhoton(myEvent.rec()->photons()[i],&PartTable.back());
+  } 
+
+  // Writing electrons (=1)
   for (unsigned int i=0;i<myEvent.rec()->electrons().size();i++)
   {
-    npart++;
-    WriteElectron(myEvent.rec()->electrons()[i],npart);
+    PartTable.push_back(LHCOParticleFormat());
+    WriteElectron(myEvent.rec()->electrons()[i],&PartTable.back());
   } 
 
-  // Writing muons
+  // Writing muons (=2)
   for (unsigned int i=0;i<myEvent.rec()->muons().size();i++)
   {
-    npart++;
-    WriteMuon(myEvent.rec()->muons()[i],npart);
+    PartTable.push_back(LHCOParticleFormat());
+    WriteMuon(myEvent.rec()->muons()[i],&PartTable.back(),myEvent.rec(),
+              myEvent.rec()->photons().size()+
+              myEvent.rec()->electrons().size()+
+              myEvent.rec()->muons().size());
   }
 
-  // Writing taus
+  // Writing taus (=3)
   for (unsigned int i=0;i<myEvent.rec()->taus().size();i++)
   {
-    npart++;
-    WriteTau(myEvent.rec()->taus()[i],npart);
+    PartTable.push_back(LHCOParticleFormat());
+    WriteTau(myEvent.rec()->taus()[i],&PartTable.back());
   }
     
-  // Writing jets
+  // Writing jets (=4)
   for (unsigned int i=0;i<myEvent.rec()->jets().size();i++)
   {
-    npart++;
-    WriteJet(myEvent.rec()->jets()[i],npart);
+    PartTable.push_back(LHCOParticleFormat());
+    WriteJet(myEvent.rec()->jets()[i],&PartTable.back());
   }
 
-  // Writing MET
-  npart++;
-  WriteMET(myEvent.rec()->MET(),npart);
-
-  return true;
-}
-
-
-/// Finalize the event
-bool LHCOWriter::WriteEventHeader(const SampleFormat& mySample,unsigned int numEvent)
-{
-  // Particle number
-  *output_ << std::setw(3) << std::right << 0;
-  *output_ << "     ";
-
-  // Event number
-  *output_ << std::setw(18) << std::right << numEvent;
-  *output_ << " ";
-
-  // Trigger word
-  *output_ << std::setw(18) << std::right << 0;
-
-  return true;
-}
+  // Writing MET (=6)
+  PartTable.push_back(LHCOParticleFormat());
+  WriteMET(myEvent.rec()->MET(),&PartTable.back());
+
+  // Printing the table
+  for (UInt_t i=0;i<PartTable.size();i++) PartTable[i].Print(i+1,output_);
+
+  // Incremeting event counter
+  counter_++;
+
+  return true;
+}
+
+
 
 
 /// Finalize the event
@@ -208,93 +207,48 @@
 }
 
 
-void LHCOWriter::WriteJet(const RecJetFormat& jet,unsigned int partNum)
+void LHCOWriter::WriteJet(const RecJetFormat& jet,
+                          LHCOParticleFormat* lhco)
 {
-  // Particle number
-  *output_ << std::setw(3) << std::right << partNum;
-  *output_ << " ";
-
-  // Particle type
-  *output_ << std::setw(3) << std::right << 4;
-  *output_ << " ";
-
-  // eta
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(jet.momentum().Eta()) 
-           << " ";
-
-  // phi
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(jet.momentum().Phi()) 
-           << " ";
-
-  // pt
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(jet.momentum().Pt())
-           << " ";
-
-  // mass
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(jet.momentum().M())
-           << " ";
-
-  // ntracks
-  *output_ << std::setw(5) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(jet.ntracks()) 
-           << " ";
-
-  // b tagging
-  if (jet.btag()) *output_ << " 1.0 ";
-  else *output_ << " 0.0 ";
-
-  // HE/EE
-  *output_ << std::setw(7) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(jet.HEoverEE()) 
-           << " ";
-
-  // dummy1 & dummy2
-  *output_ << "0.00 0.00";
-
-  // end
-  *output_ << std::endl;
+  lhco->id    = 4;
+  lhco->eta   = jet.momentum().Eta();
+  lhco->phi   = jet.momentum().Phi();
+  lhco->pt    = jet.momentum().Pt();
+  lhco->jmass = jet.momentum().M();
+  lhco->ntrk  = jet.ntracks();
+  if (jet.btag()) lhco->btag = 1.0; else lhco->btag = 0.0;
+  lhco->hadem = jet.HEoverEE();
 }
 
-void LHCOWriter::WriteMuon(const RecLeptonFormat& muon,unsigned int partNum)
+void LHCOWriter::WriteMuon(const RecLeptonFormat& muon,
+                           LHCOParticleFormat* lhco, 
+                           const RecEventFormat* myEvent,
+                           unsigned int npart)
 {
-  // Particle number
-  *output_ << std::setw(3) << std::right << partNum;
-  *output_ << " ";
-
-  // Particle type
-  *output_ << std::setw(3) << std::right << 2;
-  *output_ << " ";
-
-  // eta
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(muon.momentum().Eta())
-           << " ";
-
-  // phi
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(muon.momentum().Phi())
-           << " ";
-
-  // pt
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(muon.momentum().Pt())
-           << " ";
-
-  // mass
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(muon.momentum().M())
-           << " ";
-
-  // electric charge
-  if (muon.charge()>0) *output_ << " +1.0 ";
-  else *output_ << " -1.0 ";
-
-  // null
-  *output_ << "0.0 ";
+  lhco->id    = 2;
+  lhco->eta   = muon.momentum().Eta();
+  lhco->phi   = muon.momentum().Phi();
+  lhco->pt    = muon.momentum().Pt();
+  lhco->jmass = muon.momentum().M();
+  if (muon.charge()>0) lhco->ntrk=+1.; else lhco->ntrk=-1.;
+
+  //------------- the closest jet ---------------
+  unsigned int theClosestJet=0;
+  Double_t minDeltaR=-1;
+  for (unsigned int i=0;i<myEvent->jets().size();i++)
+  {
+    if (myEvent->jets()[i].pt()==0) continue;
+    Double_t DeltaR=muon.dr(myEvent->jets()[i]);
+    if (i==0 || DeltaR<minDeltaR)
+    {
+      theClosestJet=i;
+      minDeltaR=DeltaR;
+    }
+  }  
+  if (minDeltaR<0) lhco->btag = 0.;
+  else lhco->btag = theClosestJet+npart+1;
+
+  //---------------- isolation ------------------
 
   // isolation : sumPT_isol
   double isolation = 0;
@@ -302,147 +256,59 @@
 
   // isolation : sumET_isol
   double ET_PT = 0;
-  if (muon.sumPT_isol()!=0) ET_PT=muon.sumET_isol()/muon.sumPT_isol();
-  ET_PT=std::floor(ET_PT);
+  if (muon.pt()!=0) ET_PT=muon.sumET_isol()/muon.pt();
 
   // gathering isolation variables
-  bool test=false;
-  for (unsigned int j=0;j<5;j++)
-  {
-    ET_PT/=10;
-    if (ET_PT<1.)
-    {
-      test=true;
-      break;
-    }
-  }
-  if (!test) ET_PT=0;
+  if (ET_PT>100) ET_PT=0.99; else ET_PT=ET_PT/100.;
+
   isolation+=ET_PT;
-  *output_ << std::setw(7) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(isolation) 
-           << " ";
-
-  // dummy1 & dummy2
-  *output_ << "0.00 0.00";
-
-  // end
-  *output_ << std::endl;
-}
-
-void LHCOWriter::WriteElectron(const RecLeptonFormat& electron, unsigned int partNum)
-{
-  // Particle number
-  *output_ << std::setw(3) << std::right << partNum;
-  *output_ << " ";
-
-  // Particle type
-  *output_ << std::setw(3) << std::right << 1;
-  *output_ << " ";
-
-  // eta
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(electron.momentum().Eta()) << " ";
-
-  // phi
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(electron.momentum().Phi()) << " ";
-
-  // pt
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(electron.momentum().Pt()) << " ";
-
-  // mass
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(electron.momentum().M())  << " ";
-
-  // electric charge
-  if (electron.charge()>0) *output_ << " +1.0 ";
-  else *output_ << " -1.0 ";
-
-  // null
-  *output_ << "0.0 ";
-
-  // HE/EE
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(electron.HEoverEE()) 
-           << " ";
-
-  // dummy1 & dummy2
-  *output_ << "0.00 0.00";
-
-  // end
-  *output_ << std::endl;
-}
-
-void LHCOWriter::WriteTau(const RecTauFormat& tau, unsigned int partNum)
-{
-  // Particle number
-  *output_ << std::setw(3) << std::right << partNum;
-  *output_ << " ";
-
-  // Particle type
-  *output_ << std::setw(3) << std::right << 3;
-  *output_ << " ";
-
-  // eta
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(tau.momentum().Eta()) << " ";
-
-  // phi
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(tau.momentum().Phi()) << " ";
-
-  // pt
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(tau.momentum().Pt()) << " ";
-
-  // mass
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(tau.momentum().M())  << " ";
-
-  // electric charge
-  if (tau.charge()>0) *output_ << " +1.0 ";
-  else *output_ << " -1.0 ";
-
-  // HE/EE
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(tau.HEoverEE()) << " ";
-
-  // dummy1 & dummy2
-  *output_ << "0.00 0.00";
-
-  // end
-  *output_ << std::endl;
-}
-
-void LHCOWriter::WriteMET(const ParticleBaseFormat& met, unsigned int partNum)
-{
-  // Particle number
-  *output_ << std::setw(3) << std::right << partNum;
-  *output_ << " ";
-
-  // Particle type
-  *output_ << std::setw(3) << std::right << 6;
-  *output_ << " ";
-
-  // no eta for met
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(0) 
-           << " ";
-
-  // phi
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(met.phi())
-           << " ";
-
-  // et
-  *output_ << std::setw(18) << std::right 
-           << LHCOWriter::FortranFormat_DoublePrecision(met.et())
-           << " ";
-
-  // dummy1 & dummy2
-  *output_ << "0.00 0.00";
-
-  // end
-  *output_ << std::endl;
+  lhco->hadem = isolation;
+}
+
+void LHCOWriter::WriteElectron(const RecLeptonFormat& electron, 
+                               LHCOParticleFormat* lhco)
+{
+  lhco->id    = 1;
+  lhco->eta   = electron.momentum().Eta();
+  lhco->phi   = electron.momentum().Phi();
+  lhco->pt    = electron.momentum().Pt();
+  lhco->jmass = electron.momentum().M();
+  if (electron.charge()>0) lhco->ntrk=+1.; else lhco->ntrk=-1.;
+  lhco->btag = 0.;
+  lhco->hadem = electron.HEoverEE();
+}
+
+void LHCOWriter::WritePhoton(const RecPhotonFormat& photon, 
+                             LHCOParticleFormat* lhco)
+{
+  lhco->id    = 0;
+  lhco->eta   = photon.momentum().Eta();
+  lhco->phi   = photon.momentum().Phi();
+  lhco->pt    = photon.momentum().Pt();
+  lhco->jmass = photon.momentum().M();
+  lhco->ntrk  = 0.;
+  lhco->btag  = 0.;
+  lhco->hadem = photon.HEoverEE();
+}
+
+void LHCOWriter::WriteTau(const RecTauFormat& tau, 
+                          LHCOParticleFormat* lhco)
+{
+  lhco->id    = 3;
+  lhco->eta   = tau.momentum().Eta();
+  lhco->phi   = tau.momentum().Phi();
+  lhco->pt    = tau.momentum().Pt();
+  lhco->jmass = tau.momentum().M();
+  if (tau.charge()>0) lhco->ntrk=tau.ntracks(); else lhco->ntrk=-tau.ntracks();
+  lhco->btag = 0.;
+  lhco->hadem = tau.HEoverEE();
+}
+
+void LHCOWriter::WriteMET(const ParticleBaseFormat& met, 
+                          LHCOParticleFormat* lhco)
+{
+  lhco->id    = 6;
+  lhco->eta   = 0.;
+  lhco->phi   = met.phi();
+  lhco->pt    = met.pt();
 }
 
=== modified file 'tools/SampleAnalyzer/Writer/LHCOWriter.h'
--- tools/SampleAnalyzer/Writer/LHCOWriter.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/LHCOWriter.h	2013-06-13 12:22:26 +0000
@@ -32,6 +32,7 @@
 
 // SampleAnalyzer headers
 #include "SampleAnalyzer/Writer/WriterTextBase.h"
+#include "SampleAnalyzer/Writer/LHCOParticleFormat.h"
 
 namespace MA5
 {
@@ -44,7 +45,7 @@
   // -------------------------------------------------------------
  protected:
 
-
+  UInt_t counter_;
 
   // -------------------------------------------------------------
   //                       method members
@@ -53,7 +54,7 @@
 
   /// Constructor without argument
   LHCOWriter() : WriterTextBase()
-  {}
+  { counter_=0; }
 
 	/// Destructor
   virtual ~LHCOWriter()
@@ -71,17 +72,15 @@
  
  private:
 
-  static std::string FortranFormat_SimplePrecision(Float_t value,UInt_t precision=7); 
-  static std::string FortranFormat_DoublePrecision(Double_t value,UInt_t precision=11); 
 
   // Writing a reconstructed jet
 
-  bool WriteEventHeader(const SampleFormat& mySample,unsigned int);
-  void WriteJet(const RecJetFormat& jet,unsigned int);
-  void WriteMuon(const RecLeptonFormat& muon,unsigned int);
-  void WriteElectron(const RecLeptonFormat& electron,unsigned int);
-  void WriteTau(const RecTauFormat& tau,unsigned int);
-  void WriteMET(const ParticleBaseFormat& met,unsigned int);
+  void WriteJet(const RecJetFormat& jet, LHCOParticleFormat* lhco);
+  void WriteMuon(const RecLeptonFormat& muon, LHCOParticleFormat* lhco, const RecEventFormat* myEvent, unsigned int npart);
+  void WriteElectron(const RecLeptonFormat& electron, LHCOParticleFormat* lhco);
+  void WritePhoton(const RecPhotonFormat& photon, LHCOParticleFormat* lhco);
+  void WriteTau(const RecTauFormat& tau, LHCOParticleFormat* lhco);
+  void WriteMET(const ParticleBaseFormat& met, LHCOParticleFormat* lhco);
 
 
 };
 
=== added file 'tools/SampleAnalyzer/Writer/LHEParticleFormat.cpp'
--- tools/SampleAnalyzer/Writer/LHEParticleFormat.cpp	1970-01-01 00:00:00 +0000
+++ tools/SampleAnalyzer/Writer/LHEParticleFormat.cpp	2013-06-13 12:22:26 +0000
@@ -0,0 +1,115 @@
+////////////////////////////////////////////////////////////////////////////////
+//  
+//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
+//  The MadAnalysis development team, email: <[e-mail]>
+//  
+//  This file is part of MadAnalysis 5.
+//  Official website: <https://launchpad.net/madanalysis5>
+//  
+//  MadAnalysis 5 is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//  
+//  MadAnalysis 5 is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//  GNU General Public License for more details.
+//  
+//  You should have received a copy of the GNU General Public License
+//  along with MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
+//  
+////////////////////////////////////////////////////////////////////////////////
+
+
+// STL headers
+#include <sstream>
+#include <iomanip>
+#include <cstdlib>
+
+// SampleAnalyzer headers
+#include "SampleAnalyzer/Writer/LHEParticleFormat.h"
+
+using namespace MA5;
+
+
+// -----------------------------------------------------------------------------
+// FortranFormat_SimplePrecision
+// -----------------------------------------------------------------------------
+std::string LHEParticleFormat::FortranFormat_SimplePrecision(Float_t value,UInt_t precision)
+{
+  std::stringstream str;
+  str.precision(precision);
+  std::string word;
+
+  Bool_t negative=false;
+  if (value<0) {negative=true; value*=-1;}
+
+  Int_t exponent = 0;
+  if (value!=0)
+  {
+    for (; value > 1.0; exponent++) value/=10.;
+    for (; value < 0.0; exponent--) value*=10.;
+  }
+
+  str << std::uppercase << std::fixed << value << "E";
+  if (exponent>=0) str << "+"; else str << "-";
+  if (std::abs(exponent)<10) str << "0";
+  str << std::abs(exponent);
+  str >> word;
+  if (!negative) return word;
+  else return "-"+word;
+}
+
+
+// -----------------------------------------------------------------------------
+// FortranFormat_DoublePrecision
+// -----------------------------------------------------------------------------
+std::string LHEParticleFormat::FortranFormat_DoublePrecision(Double_t value,UInt_t precision)
+{
+  std::stringstream str;
+  str.precision(precision);
+  std::string word;
+
+  Bool_t negative=false;
+  if (value<0) {negative=true; value*=-1;}
+
+  Int_t exponent = 0;
+  if (value!=0)
+  {
+    for (; value > 1.0; exponent++) value/=10.;
+    for (; value < 0.0; exponent--) value*=10.;
+  }
+
+  str << std::uppercase << std::fixed << value << "E";
+  if (exponent>=0) str << "+"; else str << "-";
+  if (std::abs(exponent)<10) str << "0";
+  str << std::abs(exponent);
+  str >> word;
+  if (!negative) return word;
+  else return "-"+word;
+}
+
+
+// -----------------------------------------------------------------------------
+// Print a particle line
+// -----------------------------------------------------------------------------
+void LHEParticleFormat::Print(UInt_t num, std::ostream* out)
+{
+  *out << std::setw(9)  << std::right << id       << " ";
+  *out << std::setw(4)  << std::right << status   << " ";
+  *out << std::setw(4)  << std::right << mother1  << " ";
+  *out << std::setw(4)  << std::right << mother2  << " ";
+  *out << std::setw(4)  << std::right << color1   << " ";
+  *out << std::setw(4)  << std::right << color2   << " ";
+  *out << std::setw(18) << std::right << FortranFormat_DoublePrecision(px) << " ";
+  *out << std::setw(18) << std::right << FortranFormat_DoublePrecision(py) << " ";
+  *out << std::setw(18) << std::right << FortranFormat_DoublePrecision(pz) << " ";
+  *out << std::setw(18) << std::right << FortranFormat_DoublePrecision(e)  << " ";
+  *out << std::setw(18) << std::right << FortranFormat_DoublePrecision(m)  << " ";
+  *out << std::setw(2)  << std::right << std::showpoint << ctau << " ";
+  *out << std::setw(3)  << std::right << std::showpoint << spin;
+  *out << std::endl;
+}
+
+
 
=== added file 'tools/SampleAnalyzer/Writer/LHEParticleFormat.h'
--- tools/SampleAnalyzer/Writer/LHEParticleFormat.h	1970-01-01 00:00:00 +0000
+++ tools/SampleAnalyzer/Writer/LHEParticleFormat.h	2013-06-13 12:22:26 +0000
@@ -0,0 +1,71 @@
+////////////////////////////////////////////////////////////////////////////////
+//  
+//  Copyright (C) 2012-2013 Eric Conte, Benjamin Fuks
+//  The MadAnalysis development team, email: <[e-mail]>
+//  
+//  This file is part of MadAnalysis 5.
+//  Official website: <https://launchpad.net/madanalysis5>
+//  
+//  MadAnalysis 5 is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//  
+//  MadAnalysis 5 is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//  GNU General Public License for more details.
+//  
+//  You should have received a copy of the GNU General Public License
+//  along with MadAnalysis 5. If not, see <http://www.gnu.org/licenses/>
+//  
+////////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef LHE_PARTICLE_FORMAT_h
+#define LHE_PARTICLE_FORMAT_h
+
+// STL headers
+#include <iostream>
+#include <string>
+
+// ROOT headers
+#include <Rtypes.h>
+
+namespace MA5
+{
+
+class LHEParticleFormat
+{
+
+  // -------------------------------------------------------------
+  //                        data members
+  // -------------------------------------------------------------
+ public:
+
+  Int_t   id;
+  Int_t   status;
+  UInt_t  mother1;
+  UInt_t  mother2;
+  UInt_t  color1;
+  UInt_t  color2;
+  Float_t px;
+  Float_t py;
+  Float_t pz;
+  Float_t e;
+  Float_t m;
+  Float_t ctau;
+  Float_t spin;
+
+  static std::string FortranFormat_SimplePrecision(Float_t value,
+                                                   UInt_t precision=7); 
+  static std::string FortranFormat_DoublePrecision(Double_t value,
+                                                   UInt_t precision=11); 
+
+  void Print(UInt_t num, std::ostream* out);
+};
+
+
+}
+
+#endif
 
=== modified file 'tools/SampleAnalyzer/Writer/LHEWriter.cpp'
--- tools/SampleAnalyzer/Writer/LHEWriter.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/LHEWriter.cpp	2013-06-13 12:22:26 +0000
@@ -88,16 +88,14 @@
   // Opening tag
   *output_ << "<LesHouchesEvents version=""1.0"">" << std::endl;
 
-  // Header block
+  // Header tag
   *output_ << "<header>" << std::endl;
   *output_ << "<!--" << std::endl;
-  *output_ << "#*********************************************************************" << std::endl;
-  *output_ << "#                                                                    *" << std::endl;
-  *output_ << "#           This file has been produced by MadAnalysis 5             *" << std::endl;
-  *output_ << "#                                                                    *" << std::endl;
-  *output_ << "#....................................................................*" << std::endl;
-  *output_ << "" << std::endl;
-  // Explanation about the Simplified LHE
+
+  // MA5 logo
+  WriteMA5header();
+
+  // LHE format
   if (mySample.rec()!=0)
   {
     *output_ << "<MA5Format> Simplified LHE format </MA5Format>" << std::endl;
@@ -106,79 +104,106 @@
   {
     *output_ << "<MA5Format> LHE format </MA5Format>" << std::endl;
   }
+
+  // Python interface version
+  *output_ << "<MadAnalysis5Version> " << cfg_->GetPythonInterfaceVersion() 
+           << " " << cfg_->GetPythonInterfaceDate() 
+           << "</MadAnalysis5Version>" << std::endl;
+
+  // SampleAnalyzer version
+  *output_ << "<SampleAnalyzerVersion> "<< cfg_->GetSampleAnalyzerVersion()
+           << " " << cfg_->GetSampleAnalyzerVersion()
+           << " </SampleAnalyzerVersion>" << std::endl;
+
+  // Explanation about the LHE
   *output_ << "<FormatDescription>" << std::endl;
-  *output_ << "The original Les Houches Event (LHE) format is defined in hep-ph/0609017" << std::endl;
-  *output_ << "" << std::endl;
-  *output_ << "The <init> ... </init> block contains global information about the samples given as a single line:" << std::endl;
-  *output_ << " IDBMUP1 IDBMUP2 EBMUP1 EBMUP2 PDFGUP1 PDFGUP2 PDFSUP1 PDFSUP1 PDFSUP2 IDWTUP NPRUP" << std::endl;
-  *output_ << "with:" << std::endl;
-  *output_ << " - IDBMUP1: PDG code of the first beam" << std::endl;
-  *output_ << " - IDBMUP2: PDG code of the second beam" << std::endl;
-  *output_ << " - EBMUP1:  energy of the first beam" << std::endl;
-  *output_ << " - EBMUP2:  energy of the second beam" << std::endl;
-  *output_ << " - PDFGUP1: author group of the PDF employed for the first beam" << std::endl;
-  *output_ << " - PDFGUP2: author group of the PDF employed for the second beam" << std::endl;
-  *output_ << " - PDFSUP1: id of the PDF set employed for the first beam" << std::endl;
-  *output_ << " - PDFSUP2: id of the PDF set employed for the second beam" << std::endl;
-  *output_ << " - IDWTUP:  weighting strategy" << std::endl;
-  *output_ << " - NPRUP:   number of physics processes involved during the generation of the sample" << std::endl;
-  *output_ << "; the following lines give detailed process information (one line for each process):" << std::endl;
-  *output_ << " XSECUP XERRUP XMAXUP LPRUP" << std::endl;
-  *output_<<  "with:" << std::endl;
-  *output_ << " -  XSECUP: cross section" << std::endl;
-  *output_ << " -  XERRUP: cross section error" << std::endl;
-  *output_ << " -  XMAXUP: maximum event weight" << std::endl;
-  *output_ << " -  LPRUP:  process id" << std::endl;
-  *output_ << "" << std::endl;
-  *output_ << "Each event is described by an <event> ... </event> block. This block always starts by a single line containing general information on the event:" << std::endl;
-  *output_ << " NUP IDPRUP XWGTUP SCALUP AQEDUP AQCDUP" << std::endl;
-  *output_ << "with:" << std::endl;
-  *output_ << " - NUP:    number of particles" << std::endl;
-  *output_ << " - IDPRUP: process id" << std::endl;
-  *output_ << " - XWGTUP: event weight" << std::endl;
-  *output_ << " - SCALUP: scale" << std::endl;
-  *output_ << " - AQEDUP: alpha QED" << std::endl;
-  *output_ << " - AQCDUP: alpha QCD" << std::endl;
-  *output_ << "This line is then followed by one line for each particle in the event:" << std::endl;
-  *output_ << " IDUP ISTUP MOTHUP1 MOTHUP2 ICOLUP1 ICOLUP2 PUP1 PUP2 PUP3 PUP4 PUP5 VTIMUP SPINUP" << std::endl;
-  *output_ << "with:" << std::endl;
-  *output_ << " - IDUP:    PDG code" << std::endl;
-  *output_ << " - ISTUP:   status code" << std::endl;
-  *output_ << " - MOTHUP1: index of the first mother particle" << std::endl;
-  *output_ << " - MOTHUP2: index of the second mother particle" << std::endl;
-  *output_ << " - ICOLUP1: first color tag" << std::endl;
-  *output_ << " - ICOLUP2: second color tag" << std::endl;
-  *output_ << " - PUP1:    PX [GeV/c]" << std::endl;
-  *output_ << " - PUP2:    PY [GeV/c]" << std::endl;
-  *output_ << " - PUP3:    PZ [GeV/c]" << std::endl;
-  *output_ << " - PUP4:    E  [GeV]" << std::endl;
-  *output_ << " - PUP5:    M  [GeV/c^2] (a space-like virtuality is denoted by a negative mass)" << std::endl;
-  *output_ << " - VTIMUP:  c tau" << std::endl;
-  *output_ << " - SPINUP:  cosine of the angle between the spin vector of the particle and its three-momentum" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
+  *output_ << "# The original Les Houches Event (LHE) format is defined in hep-ph/0609017      #" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
+  *output_ << "# The <init> ... </init> block contains global information about the samples    #" << std::endl;
+  *output_ << "# given as a single line:                                                       #" << std::endl;
+  *output_ << "#    IDBM1 IDBM2 EBM1 EBM2 PDFG1 PDFG2 PDFS1 PDFS1 PDFS2 IDWT NPR               #" << std::endl;
+  *output_ << "# with:                                                                         #" << std::endl;
+  *output_ << "#   - IDBM1: PDG code of the first beam.                                        #" << std::endl;
+  *output_ << "#   - IDBM2: PDG code of the second beam.                                       #" << std::endl;
+  *output_ << "#   - EBM1:  energy of the first beam.                                          #" << std::endl;
+  *output_ << "#   - EBM2:  energy of the second beam.                                         #" << std::endl;
+  *output_ << "#   - PDFG1: author group of the PDF employed for the first beam.               #" << std::endl;
+  *output_ << "#   - PDFG2: author group of the PDF employed for the second beam.              #" << std::endl;
+  *output_ << "#   - PDFS1: id of the PDF set employed for the first beam.                     #" << std::endl;
+  *output_ << "#   - PDFS2: id of the PDF set employed for the second beam.                    #" << std::endl;
+  *output_ << "#   - IDWT:  weighting strategy.                                                #" << std::endl;
+  *output_ << "#   - NPR:   number of physics processes involved during the generation of      #" << std::endl;
+  *output_ << "#            the sample.                                                        #" << std::endl;
+  *output_ << "# The following lines give detailed process information (one line for each      #" << std::endl;
+  *output_ << "# process):                                                                     #" << std::endl;
+  *output_ << "#    XSEC XERR XMAX LPR                                                         #" << std::endl;
+  *output_ << "# with:                                                                         #" << std::endl;
+  *output_ << "#   -  XSEC: cross section                                                      #" << std::endl;
+  *output_ << "#   -  XERR: cross section error                                                #" << std::endl;
+  *output_ << "#   -  XMAX: maximum event weight                                               #" << std::endl;
+  *output_ << "#   -  LPR:  process id                                                         #" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
+  *output_ << "# Each event is described by an <event> ... </event> block. This block always   #" << std::endl;
+  *output_ << "# starts by a single line containing general information on the event:          #" << std::endl;
+  *output_ << "#    N IDPR XWGT SCAL AQED AQCD                                                 #" << std::endl;
+  *output_ << "# with:                                                                         #" << std::endl;
+  *output_ << "#   - N:    number of particles                                                 #" << std::endl;
+  *output_ << "#   - IDPR: process id                                                          #" << std::endl;
+  *output_ << "#   - XWGT: event weight                                                        #" << std::endl;
+  *output_ << "#   - SCAL: scale                                                               #" << std::endl;
+  *output_ << "#   - AQED: alpha QED                                                           #" << std::endl;
+  *output_ << "#   - AQCD: alpha QCD                                                           #" << std::endl;
+  *output_ << "# This line is then followed by one line for each particle in the event:        #" << std::endl;
+  *output_ << "#    ID IST MOTH1 MOTH2 ICOL1 ICOL2 P1 P2 P3 P4 P5 VTIM SPIN                    #" << std::endl;
+  *output_ << "# with:                                                                         #" << std::endl;
+  *output_ << "#   - ID:    PDG code                                                           #" << std::endl;
+  *output_ << "#   - IST:   status code                                                        #" << std::endl;
+  *output_ << "#   - MOTH1: row number corresponding to the first mother particle              #" << std::endl;
+  *output_ << "#   - MOTH2: row number corresponding to the second mother particle             #" << std::endl;
+  *output_ << "#   - ICOL1: first color tag                                                    #" << std::endl;
+  *output_ << "#   - ICOL2: second color tag                                                   #" << std::endl;
+  *output_ << "#   - P1:    PX                                                                 #" << std::endl;
+  *output_ << "#   - P2:    PY                                                                 #" << std::endl;
+  *output_ << "#   - P3:    PZ                                                                 #" << std::endl;
+  *output_ << "#   - P4:    E                                                                  #" << std::endl;
+  *output_ << "#   - P5:    M (a space-like virtuality is denoted by a negative mass)          #" << std::endl;
+  *output_ << "#   - VTIM:  c tau                                                              #" << std::endl;
+  *output_ << "#   - SPIN:  cosine of the angle between the spin vector of the particle and    #" << std::endl;
+  *output_ << "#            its three-momentum                                                 #" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
 
   // Explanation about the Simplified LHE
   if (mySample.rec()!=0)
   {
-    *output_ << "" << std::endl;
-    *output_ << "In the 'simplified LHE' format, there are three types of objects classified according to their statuscode:" << std::endl;
-    *output_ << " - Objects with StatusCode = -1: initial interacting partons" << std::endl;
-    *output_ << " - Objects with StatusCode = +3: particles produced during the hard process" << std::endl;
-    *output_ << " - Objects with StatusCode = +1: reconstructed objects (after applying a fast simulation of a (possibly perfect) detector and a jet-clustering algorithm" << std::endl;
-    *output_ << "When MadAnalysis is in charge of the reconstruction (i.e., applying the jet-clustering algorithm), the particle codes follow the conventions:" << std::endl;
-    *output_ << " - particle with a PDG code = +11 or -11: electrons and positrons." <<std::endl;
-    *output_ << "   They can be isolated or not as well aspossibly issued from the hadronization process." << std::endl;
-    *output_ << " - particle with a PDG code = +13 or -13: muons and antimuons." << std::endl;
-    *output_ << "   They can be isolated or not as well aspossibly issued from the hadronization process." << std::endl;
-    *output_ << " - particle with a PDG code = +15 or -15: hadronically decaying (anti)taus." << std::endl;
-    *output_ << "   These consist of jets matching a hadronically decaying tau when inspecting the Monte Carlo history." << std::endl;
-    *output_ << "   (Mis)Identification efficiency can be possibly included." << std::endl;
-    *output_ << " - particle with a PDG code = 5: b-jets." << std::endl;
-    *output_ << "   These consist of jets matching a b-quark when inspecting the Monte Carlo history." << std::endl;
-    *output_ << "   (Mis)Identification efficiency can be possibly included." << std::endl;
-    *output_ << " - particle with a PDG code = 21: jets which are not b-tagged and taus which are not tau-tagged." << std::endl;
-    *output_ << "   the jet collection includes also electrons collection and hadronic taus collection" << std::endl;
-    *output_ << " - particle with a PDG code = 12: the missing transverse energy." << std::endl;
-    *output_ << "   The missing transverse energy is computed as opposite to the sum of the four-momenta of all jets, electrons, muons and hadronic taus." << std::endl;
+    *output_ << "# In the 'simplified LHE' format, there are three types of objects classified   #" << std::endl;
+    *output_ << "# according to their statuscode:                                                #" << std::endl;
+    *output_ << "#   - objects with StatusCode = -1: initial interacting partons.                #" << std::endl;
+    *output_ << "#   - objects with StatusCode = +3: particles produced during the hard process. #" << std::endl;
+    *output_ << "#   - objects with StatusCode = +1: physics objects reconstructed by a fast     #" << std::endl;
+    *output_ << "#                                   detector simulation (or perfect detector).  #" << std::endl;
+    *output_ << "# When MadAnalysis is in charge of the reconstruction (i.e., applying the       #" << std::endl;
+    *output_ << "# jet-clustering algorithm), the particle codes follow the conventions:         #" << std::endl;
+    *output_ << "#   - particle with a PDG code = +11 or -11: electrons and positrons.           #" << std::endl;
+    *output_ << "#     They can be isolated or not as well as possibly issued from the           #" << std::endl;
+    *output_ << "#     hadronization process.                                                    #" << std::endl;
+    *output_ << "#   - particle with a PDG code = +13 or -13: muons and antimuons.               #" << std::endl;
+    *output_ << "#     They can be isolated or not as well as possibly issued from the           #" << std::endl;
+    *output_ << "#     hadronization process.                                                    #" << std::endl;
+    *output_ << "#   - particle with a PDG code = +15 or -15: hadronically decaying (anti)taus.  #" << std::endl;
+    *output_ << "#     These consist of jets matching a hadronically-decaying tau when           #" << std::endl;
+    *output_ << "#     inspecting the Monte Carlo history. (Mis)Identification efficiency can be #" << std::endl;
+    *output_ << "#     possibly included.                                                        #" << std::endl;
+    *output_ << "#   - particle with a PDG code = 5: b-jets.                                     #" << std::endl;
+    *output_ << "#     These consist of jets matching a b-quark when inspecting the Monte Carlo  #" << std::endl;
+    *output_ << "#     history. (Mis)Identification efficiency can be possibly included.         #" << std::endl;
+    *output_ << "#   - particle with a PDG code = 21: jets which are not b-tagged and taus which #" << std::endl;
+    *output_ << "#     are not tau-tagged. The jet collection includes also electrons collection #" << std::endl;
+    *output_ << "#     and hadronic taus collection                                              #" << std::endl;
+    *output_ << "#   - particle with a PDG code = 12: the missing transverse energy.             #" << std::endl;
+    *output_ << "#     The missing transverse energy is computed as opposite to the sum of the   #" << std::endl;
+    *output_ << "#     four-momenta of all jets, electrons, muons and hadronic taus.             #" << std::endl;
+    *output_ << "#################################################################################" << std::endl;
   }
   *output_ << "</FormatDescription>" << std::endl;
   if (mySample.mc()!=0)
@@ -279,59 +304,43 @@
   // Event header
   *output_ << "<event>" << std::endl;
 
-  unsigned int counter = 0;
+  // Container for particles
+  std::vector<LHEParticleFormat> particles;
+  UInt_t counter=0;
 
   // Writing MC particles : only MC info case
   // -> hypothesis : input = LHE
   if (myEvent.mc()!=0 && myEvent.rec()==0)
   {
-    for (unsigned int i=0;i<myEvent.mc()->particles().size();i++) counter ++;
+    counter += myEvent.mc()->particles().size();
   }
 
   // Writing MC particles : MC+REC info case
   // -> hypothesis : input = HEP
-  if (myEvent.mc()!=0 && myEvent.rec()!=0)
+  else if (myEvent.mc()!=0 && myEvent.rec()!=0)
   {
-    for (unsigned int i=4;i<myEvent.mc()->particles().size();i++)
+    for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
+    {
       if (myEvent.mc()->particles()[i].statuscode()==3 || 
           ( myEvent.mc()->particles()[i].statuscode()>=21 &&
             myEvent.mc()->particles()[i].statuscode()<=29)
           ) counter++;
+    }
+
   }
 
   // Writing REC particles
   if (myEvent.rec()!=0)
   {
-    if (myEvent.mc()==0)
-    {
-      for (unsigned int i=0;i<myEvent.rec()->muons().size();i++)
-        counter++;
-      for (unsigned int i=0;i<myEvent.rec()->electrons().size();i++)
-        counter++;
-      for (unsigned int i=0;i<myEvent.rec()->taus().size();i++)
-        counter++;
-    }
-    else
-    {
-      for (unsigned int i=6;i<myEvent.mc()->particles().size();i++)
-      {
-        if ( (myEvent.mc()->particles()[i].statuscode()==3 || 
-              ( myEvent.mc()->particles()[i].statuscode()>=21 &&
-                myEvent.mc()->particles()[i].statuscode()<=29))&& (
-            fabs(myEvent.mc()->particles()[i].pdgid())==11 ||
-            fabs(myEvent.mc()->particles()[i].pdgid())==13 ||
-            fabs(myEvent.mc()->particles()[i].pdgid())==15 ))
-          {
-            counter++;
-          }
-      }
-    }
-    for (unsigned int i=0;i<myEvent.rec()->jets().size();i++)
-      counter++;
-    counter++;
+    counter += myEvent.rec()->muons().size() + 
+               myEvent.rec()->electrons().size() +
+               myEvent.rec()->taus().size() + 
+               myEvent.rec()->photons().size() +
+               myEvent.rec()->jets().size() + 1 /*MET*/;
   }
 
   // Writing event global information
+  particles.reserve(counter);
   WriteEventHeader(myEvent,counter);
 
   // Writing MC particles : only MC info case
@@ -339,9 +348,14 @@
   if (myEvent.mc()!=0 && myEvent.rec()==0)
   {
     for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
+    {
+      particles.push_back(LHEParticleFormat());
       WriteParticle(myEvent.mc()->particles()[i],
                     myEvent.mc()->particles()[i].mothup1_,
-                    myEvent.mc()->particles()[i].mothup2_);
+                    myEvent.mc()->particles()[i].mothup2_,
+                    0,
+                    particles.back());
+    }
   }
 
   // Writing MC particles : MC+REC info case
@@ -350,20 +364,22 @@
   {
     for (unsigned int i=0;i<myEvent.mc()->particles().size();i++)
     {
-      if ( myEvent.mc()->particles()[i].statuscode()>=11 &&
-           myEvent.mc()->particles()[i].statuscode()<=19 )
-         WriteParticle(myEvent.mc()->particles()[i],0,0,-1);
-
-      else if (i>3 && i<6 && myEvent.mc()->particles()[i].statuscode()==3)
-        WriteParticle(myEvent.mc()->particles()[i],0,0,-1);
+      if ( myEvent.mc()->particles()[i].statuscode()==-1 ||
+           ( myEvent.mc()->particles()[i].statuscode()>=11 &&
+             myEvent.mc()->particles()[i].statuscode()<=19 ))
+      {
+        particles.push_back(LHEParticleFormat());
+        WriteParticle(myEvent.mc()->particles()[i],0,0,-1, particles.back());
+      }
 
       else if (myEvent.mc()->particles()[i].statuscode()==3 || 
            ( myEvent.mc()->particles()[i].statuscode()>=21 &&
              myEvent.mc()->particles()[i].statuscode()<=29))
       {
+          particles.push_back(LHEParticleFormat());
           Int_t moth1 = GetMotherIndex(myEvent.mc()->particles()[i].mothup1_);
           Int_t moth2 = GetMotherIndex(myEvent.mc()->particles()[i].mothup2_);
-          WriteParticle(myEvent.mc()->particles()[i],moth1,moth2,3);
+          WriteParticle(myEvent.mc()->particles()[i],moth1,moth2,3, particles.back());
       }
     } 
   }
@@ -371,40 +387,39 @@
   // Writing REC particles
   if (myEvent.rec()!=0)
   {
-    //    if (myEvent.mc()==0)
-    {
-      for (unsigned int i=0;i<myEvent.rec()->muons().size();i++)
-        WriteMuon(myEvent.rec()->muons()[i]);
-      for (unsigned int i=0;i<myEvent.rec()->electrons().size();i++)
-        WriteElectron(myEvent.rec()->electrons()[i]);
-      for (unsigned int i=0;i<myEvent.rec()->taus().size();i++)
-        WriteTau(myEvent.rec()->taus()[i]);
-    }
-    /*
-    else
-    {
-      for (unsigned int i=6;i<myEvent.mc()->particles().size();i++)
-      {
-        if ((myEvent.mc()->particles()[i].statuscode()==3 || 
-            ( myEvent.mc()->particles()[i].statuscode()>=21 &&
-              myEvent.mc()->particles()[i].statuscode()<=29)) && (
-            fabs(myEvent.mc()->particles()[i].pdgid())==11 ||
-            fabs(myEvent.mc()->particles()[i].pdgid())==13 ||
-            fabs(myEvent.mc()->particles()[i].pdgid())==15 ))
-          {
-            WriteParticle(myEvent.mc()->particles()[i],0,0,1);
-          }
-      }
-    }*/
-
+    for (unsigned int i=0;i<myEvent.rec()->muons().size();i++)
+    {
+      particles.push_back(LHEParticleFormat());
+      WriteMuon(myEvent.rec()->muons()[i],particles.back());
+    }
+    for (unsigned int i=0;i<myEvent.rec()->electrons().size();i++)
+    {
+      particles.push_back(LHEParticleFormat());
+      WriteElectron(myEvent.rec()->electrons()[i],particles.back());
+    }
+    for (unsigned int i=0;i<myEvent.rec()->taus().size();i++)
+    {
+      particles.push_back(LHEParticleFormat());
+      WriteTau(myEvent.rec()->taus()[i],particles.back());
+    }
     for (unsigned int i=0;i<myEvent.rec()->jets().size();i++)
-      WriteJet(myEvent.rec()->jets()[i]);
-    WriteMET(myEvent.rec()->MET());
+    {
+      particles.push_back(LHEParticleFormat());
+      WriteJet(myEvent.rec()->jets()[i],particles.back());
+    }
+    for (unsigned int i=0;i<myEvent.rec()->photons().size();i++)
+    {
+      particles.push_back(LHEParticleFormat());
+      WritePhoton(myEvent.rec()->photons()[i],particles.back());
+    }
+    particles.push_back(LHEParticleFormat());
+    WriteMET(myEvent.rec()->MET(),particles.back());
   }
 
   // Event foot
+  for (unsigned int i=0;i<particles.size();i++) particles[i].Print(i+1, output_);
   *output_ << "</event>" << std::endl;
-  return true;
+  return true; 
 }
 
 
@@ -447,117 +462,129 @@
 
 
 /// Writing a particle
-bool LHEWriter::WriteParticle(const MCParticleFormat& myPart, Int_t mother1, Int_t mother2, Int_t statuscode )
-{
-  *output_ << std::setw(9)  << std::right << myPart.pdgid_ << " ";
-  if (statuscode!=0)   *output_ << std::setw(4)  << std::right << statuscode << " ";
-  else *output_ << std::setw(4)  << std::right << myPart.statuscode_ << " ";
-  *output_ << std::setw(4)  << std::right << mother1 /*myPart.mothup1_*/ << " ";
-  *output_ << std::setw(4)  << std::right << mother2 /*myPart.mothup2_*/ << " ";
-  *output_ << std::setw(4)  << std::right << 0 << " ";
-  *output_ << std::setw(4)  << std::right << 0 << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(myPart.momentum_.Px()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(myPart.momentum_.Py()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(myPart.momentum_.Pz()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(myPart.momentum_.E())  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(myPart.momentum_.M())  << " "; 
-  *output_ << std::setw(2)  << std::right << std::showpoint << myPart.ctau_ << " ";
-  *output_ << std::setw(3)  << std::right << std::showpoint << myPart.spin_;
-  *output_ << std::endl;
-  return true;
-}
-
-
-void LHEWriter::WriteJet(const RecJetFormat& jet)
-{
-  if (jet.btag()) *output_ << std::setw(9)  << std::right << 5 << " ";
-  else *output_ << std::setw(9)  << std::right << 21 << " ";
-  *output_ << std::setw(4)  << std::right << 1  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(jet.momentum().Px()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(jet.momentum().Py()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(jet.momentum().Pz()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(jet.momentum().E())  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(jet.momentum().M())  << " "; 
-  *output_ << std::setw(2)  << std::right << std::showpoint << 0. << " ";
-  *output_ << std::setw(3)  << std::right << std::showpoint << 0.;
-  *output_ << std::endl;
-}
-
-void LHEWriter::WriteMuon(const RecLeptonFormat& muon)
-{
-  if (muon.charge()>0) *output_ << std::setw(9)  << std::right << -13 << " ";
-  else *output_ << std::setw(9)  << std::right << +13 << " ";
-  *output_ << std::setw(4)  << std::right << 1  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(muon.momentum().Px()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(muon.momentum().Py()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(muon.momentum().Pz()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(muon.momentum().E())  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(muon.momentum().M())  << " "; 
-  *output_ << std::setw(2)  << std::right << std::showpoint << 0. << " ";
-  *output_ << std::setw(3)  << std::right << std::showpoint << 0.;
-  *output_ << std::endl;
-}
-
-void LHEWriter::WriteElectron(const RecLeptonFormat& electron)
-{
-  if (electron.charge()>0) *output_ << std::setw(9)  << std::right << -11 << " ";
-  else *output_ << std::setw(9)  << std::right << +11 << " ";
-  *output_ << std::setw(4)  << std::right << 1  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(electron.momentum().Px()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(electron.momentum().Py()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(electron.momentum().Pz()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(electron.momentum().E())  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(electron.momentum().M())  << " "; 
-  *output_ << std::setw(2)  << std::right << std::showpoint << 0. << " ";
-  *output_ << std::setw(3)  << std::right << std::showpoint << 0.;
-  *output_ << std::endl;
-}
-
-void LHEWriter::WriteTau(const RecTauFormat& tau)
-{
-  if (tau.charge()>0) *output_ << std::setw(9)  << std::right << -15 << " ";
-  else *output_ << std::setw(9)  << std::right << +15 << " ";
-  *output_ << std::setw(4)  << std::right << 1  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(tau.momentum().Px()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(tau.momentum().Py()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(tau.momentum().Pz()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(tau.momentum().E())  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(tau.momentum().M())  << " "; 
-  *output_ << std::setw(2)  << std::right << std::showpoint << 0. << " ";
-  *output_ << std::setw(3)  << std::right << std::showpoint << 0.;
-  *output_ << std::endl;
-}
-
-void LHEWriter::WriteMET(const ParticleBaseFormat& met)
-{
-  *output_ << std::setw(9)  << std::right << 12 << " ";
-  *output_ << std::setw(4)  << std::right << 1  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(4)  << std::right << 0  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(met.px()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(met.py()) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(0.) << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(met.pt())  << " ";
-  *output_ << std::setw(18) << std::right << LHEWriter::FortranFormat_DoublePrecision(0.)  << " "; 
-  *output_ << std::setw(2)  << std::right << std::showpoint << 0. << " ";
-  *output_ << std::setw(3)  << std::right << std::showpoint << 0.;
-  *output_ << std::endl;
+void LHEWriter::WriteParticle(const MCParticleFormat& myPart, 
+                              Int_t mother1, Int_t mother2, 
+                              Int_t statuscode, LHEParticleFormat& lhe)
+{
+  if (statuscode!=0) lhe.status = statuscode;
+  else lhe.status = myPart.statuscode_;
+  lhe.id      = myPart.pdgid_;
+  lhe.mother1 = mother1;
+  lhe.mother2 = mother2;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = myPart.momentum().Px();
+  lhe.py      = myPart.momentum().Py();
+  lhe.pz      = myPart.momentum().Pz();
+  lhe.e       = myPart.momentum().E();
+  lhe.m       = myPart.momentum().M();
+  lhe.ctau    = myPart.ctau_;
+  lhe.spin    = myPart.spin_;
+}
+
+
+void LHEWriter::WriteJet(const RecJetFormat& jet, LHEParticleFormat& lhe)
+{
+  if (jet.btag()) lhe.id = 5; else lhe.id = 21;
+  lhe.status  = 1;
+  lhe.mother1 = 0;
+  lhe.mother2 = 0;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = jet.momentum().Px();
+  lhe.py      = jet.momentum().Py();
+  lhe.pz      = jet.momentum().Pz();
+  lhe.e       = jet.momentum().E();
+  lhe.m       = jet.momentum().M();
+  lhe.ctau    = 0.;
+  lhe.spin    = 0.;
+}
+
+
+void LHEWriter::WriteMuon(const RecLeptonFormat& muon, LHEParticleFormat& lhe)
+{
+  if (muon.charge()>0) lhe.id = -13; else lhe.id = +13;
+  lhe.status  = 1;
+  lhe.mother1 = 0;
+  lhe.mother2 = 0;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = muon.momentum().Px();
+  lhe.py      = muon.momentum().Py();
+  lhe.pz      = muon.momentum().Pz();
+  lhe.e       = muon.momentum().E();
+  lhe.m       = muon.momentum().M();
+  lhe.ctau    = 0.;
+  lhe.spin    = 0.;
+}
+
+void LHEWriter::WriteElectron(const RecLeptonFormat& electron, LHEParticleFormat& lhe)
+{
+  if (electron.charge()>0) lhe.id = -11; else lhe.id = +11;
+  lhe.status  = 1;
+  lhe.mother1 = 0;
+  lhe.mother2 = 0;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = electron.momentum().Px();
+  lhe.py      = electron.momentum().Py();
+  lhe.pz      = electron.momentum().Pz();
+  lhe.e       = electron.momentum().E();
+  lhe.m       = electron.momentum().M();
+  lhe.ctau    = 0.;
+  lhe.spin    = 0.;
+}
+
+
+void LHEWriter::WritePhoton(const RecPhotonFormat& photon, LHEParticleFormat& lhe)
+{
+  lhe.id      = 22;
+  lhe.status  = 1;
+  lhe.mother1 = 0;
+  lhe.mother2 = 0;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = photon.momentum().Px();
+  lhe.py      = photon.momentum().Py();
+  lhe.pz      = photon.momentum().Pz();
+  lhe.e       = photon.momentum().E();
+  lhe.m       = photon.momentum().M();
+  lhe.ctau    = 0.;
+  lhe.spin    = 0.;
+}
+
+
+void LHEWriter::WriteTau(const RecTauFormat& tau, LHEParticleFormat& lhe)
+{
+  if (tau.charge()>0) lhe.id = -15; else lhe.id = +15;
+  lhe.status  = 1;
+  lhe.mother1 = 0;
+  lhe.mother2 = 0;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = tau.momentum().Px();
+  lhe.py      = tau.momentum().Py();
+  lhe.pz      = tau.momentum().Pz();
+  lhe.e       = tau.momentum().E();
+  lhe.m       = tau.momentum().M();
+  lhe.ctau    = 0.;
+  lhe.spin    = 0.;
+}
+
+
+void LHEWriter::WriteMET(const ParticleBaseFormat& met, LHEParticleFormat& lhe)
+{
+  lhe.id      = 12;
+  lhe.status  = 1;
+  lhe.mother1 = 0;
+  lhe.mother2 = 0;
+  lhe.color1  = 0;
+  lhe.color2  = 0;
+  lhe.px      = met.px();
+  lhe.py      = met.py();
+  lhe.pz      = 0.;
+  lhe.e       = met.pt();
+  lhe.m       = 0.;
+  lhe.ctau    = 0.;
+  lhe.spin    = 0.;
 }
 
=== modified file 'tools/SampleAnalyzer/Writer/LHEWriter.h'
--- tools/SampleAnalyzer/Writer/LHEWriter.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/LHEWriter.h	2013-06-13 12:22:26 +0000
@@ -32,6 +32,7 @@
 
 // SampleAnalyzer headers
 #include "SampleAnalyzer/Writer/WriterTextBase.h"
+#include "SampleAnalyzer/Writer/LHEParticleFormat.h"
 
 namespace MA5
 {
@@ -76,17 +77,19 @@
                         unsigned int nevents);
 
   /// Writing a particle
-  bool WriteParticle(const MCParticleFormat& myPart, Int_t mother1, Int_t mother2, Int_t statuscode=0);
+  void WriteParticle(const MCParticleFormat& myPart, Int_t mother1, Int_t mother2, 
+                     Int_t statuscode, LHEParticleFormat& lhe);
 
   static std::string FortranFormat_SimplePrecision(Float_t value,UInt_t precision=7); 
   static std::string FortranFormat_DoublePrecision(Double_t value,UInt_t precision=11); 
 
   // Writing a reconstructed jet
-  void WriteJet(const RecJetFormat& jet);
-  void WriteMuon(const RecLeptonFormat& muon);
-  void WriteElectron(const RecLeptonFormat& electron);
-  void WriteTau(const RecTauFormat& tau);
-  void WriteMET(const ParticleBaseFormat& met);
+  void WriteJet(const RecJetFormat& jet, LHEParticleFormat& lhe);
+  void WriteMuon(const RecLeptonFormat& muon, LHEParticleFormat& lhe);
+  void WriteElectron(const RecLeptonFormat& electron, LHEParticleFormat& lhe);
+  void WritePhoton(const RecPhotonFormat& photon, LHEParticleFormat& lhe);
+  void WriteTau(const RecTauFormat& tau, LHEParticleFormat& lhe);
+  void WriteMET(const ParticleBaseFormat& met, LHEParticleFormat& lhe);
 
 
 };
 
=== modified file 'tools/SampleAnalyzer/Writer/WriterBase.h'
--- tools/SampleAnalyzer/Writer/WriterBase.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/WriterBase.h	2013-06-13 12:22:26 +0000
@@ -36,6 +36,7 @@
 #include "SampleAnalyzer/DataFormat/EventFormat.h"
 #include "SampleAnalyzer/DataFormat/SampleFormat.h"
 #include "SampleAnalyzer/Service/Physics.h"
+#include "SampleAnalyzer/Core/Configuration.h"
 
 // ROOT headers
 #include <TVector.h>
@@ -77,7 +78,8 @@
   }
 
   /// Initialize (virtual pure)
-  virtual bool Initialize(const std::string& filename) = 0;
+  virtual bool Initialize(const Configuration* cfg,
+                          const std::string& filename) = 0;
 
   /// Read the sample (virtual pure)
   virtual bool WriteHeader(const SampleFormat& mySample) = 0;
 
=== modified file 'tools/SampleAnalyzer/Writer/WriterTextBase.cpp'
--- tools/SampleAnalyzer/Writer/WriterTextBase.cpp	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/WriterTextBase.cpp	2013-06-13 12:22:26 +0000
@@ -37,8 +37,12 @@
 // -----------------------------------------------------------------------------
 // Initialize
 // -----------------------------------------------------------------------------
-bool WriterTextBase::Initialize(const std::string& rawfilename)
+bool WriterTextBase::Initialize(const Configuration* cfg,
+                                const std::string& rawfilename)
 {
+  // Saving configuration file
+  cfg_ = cfg;
+
   // Is the file stored in Rfio
   rfio_ = ReaderBase::IsRfioMode(rawfilename);
 
@@ -147,3 +151,26 @@
   return true;
 }
 
+
+// -----------------------------------------------------------------------------
+// Header
+// -----------------------------------------------------------------------------
+void WriterTextBase::WriteMA5header()
+{
+  *output_ << "#################################################################################" << std::endl;
+  *output_ << "#                 THIS FILE HAS BEEN PRODUCED BY MADANALYSIS 5                  #" << std::endl;
+  *output_ << "#                                   ______  ______                              #" << std::endl;
+  *output_ << "#                           /'\\_/`\\/\\  __ \\/\\  ___\\                             #" << std::endl;
+  *output_ << "#                          /\\      \\ \\ \\_\\ \\ \\ \\__/                             #" << std::endl;
+  *output_ << "#                          \\ \\ \\__\\ \\ \\  __ \\ \\___``\\                           #" << std::endl;
+  *output_ << "#                           \\ \\ \\_/\\ \\ \\ \\/\\ \\/\\ \\_\\ \\                          #" << std::endl;
+  *output_ << "#                            \\ \\_\\\\ \\_\\ \\_\\ \\_\\ \\____/                          #" << std::endl;
+  *output_ << "#                             \\/_/ \\/_/\\/_/\\/_/\\/___/                           #" << std::endl;
+  *output_ << "#                                                                               #" << std::endl;
+  //  *output_ << "#   MA5 release : " + "%-24s" % main.version + "%+15s" % main.date  + "         #" << std::endl;
+  //  *output_ << "#                                                                               #" << std::endl;
+  *output_ << "#                    Comput. Phys. Commun. 184 (2013) 222-256                   #" << std::endl;
+  *output_ << "#             The MadAnalysis Development Team - Please visit us at             #" << std::endl;
+  *output_ << "#                      https://launchpad.net/madanalysis5                       #" << std::endl;
+  *output_ << "#################################################################################" << std::endl;
+}
 
=== modified file 'tools/SampleAnalyzer/Writer/WriterTextBase.h'
--- tools/SampleAnalyzer/Writer/WriterTextBase.h	2013-05-17 06:19:22 +0000
+++ tools/SampleAnalyzer/Writer/WriterTextBase.h	2013-06-13 12:22:26 +0000
@@ -52,6 +52,8 @@
   /// Streaming for reading input
   std::ostream* output_;
 
+  /// Configuration
+  const Configuration* cfg_;
 
   // -------------------------------------------------------------
   //                       method members
@@ -71,7 +73,8 @@
   }
 
   /// Initialize
-  virtual bool Initialize(const std::string& filename);
+  virtual bool Initialize(const Configuration* cfg,
+                          const std::string& filename);
 
   /// Read the sample (virtual pure)
   virtual bool WriteHeader(const SampleFormat& mySample) = 0;
@@ -85,6 +88,9 @@
 
   /// Finalize
   virtual bool Finalize();
+
+  /// MA5 logo
+  void WriteMA5header();
  
 };