Merge pull request alisw#162 from hcab14/master

added reader for Les Houches event record format

EVGEN/AliGenReaderLHE.cxx

@@ -0,0 +1,102 @@
+#include <TVirtualMC.h>
+#include <TDatabasePDG.h>
+#include <TLorentzVector.h>
+#include <TParticle.h>
+
+#include "AliLog.h"
+#include "AliGenReaderLHE.h"
+#include "AliRun.h"
+#include "AliStack.h"
+
+ClassImp(AliGenReaderLHE);
+
+AliGenReaderLHE::~AliGenReaderLHE() {
+  if (fXMLDoc) {
+    fXMLEngine.FreeDoc(fXMLDoc);
+    fXMLDoc = NULL;
+  }
+}
+
+void AliGenReaderLHE::Init()
+{
+  if (fXMLDoc)
+    fXMLEngine.FreeDoc(fXMLDoc);
+
+  fXMLDoc = fXMLEngine.ParseFile(fFileName);
+  if (!fXMLDoc)
+    AliFatalF("fXMLEngine.ParseFile('%s') failed", fFileName);
+
+  XMLNodePointer_t xmlMainNode = fXMLEngine.DocGetRootElement(fXMLDoc);
+
+  if (TString(fXMLEngine.GetNodeName(xmlMainNode)) != "LesHouchesEvents")
+    AliFatal("main node != 'LesHouchesEvents'");
+
+  XMLAttrPointer_t attr = fXMLEngine.GetFirstAttr(xmlMainNode);
+  if (!attr)
+    AliFatal("cannot determine version of Les Houches file format");
+
+  if (TString(fXMLEngine.GetAttrName(attr))  != "version" ||
+      TString(fXMLEngine.GetAttrValue(attr)) != "1.0")
+    AliFatal("version != 1.0");
+
+  Bool_t foundEvent = kFALSE;
+  fXMLChild = fXMLEngine.GetChild(xmlMainNode);
+  while (fXMLChild && !foundEvent) {
+    foundEvent = (TString(fXMLEngine.GetNodeName(fXMLChild)) == "event");
+    fXMLChild  = fXMLEngine.GetNext(fXMLChild);
+  }
+  if (!fXMLChild || !foundEvent)
+    AliFatal("no events found");
+}
+
+Int_t AliGenReaderLHE::NextEvent()
+{
+  const char* content = fXMLEngine.GetNodeContent(fXMLChild);
+  if (!content) {
+    AliError("empty event");
+    return 0;
+  }
+  fStrStream.str(content);
+
+  Int_t nup=0, idprup=0;
+  Double_t xwgtup=0, scalup=0, aqedup=0, aqcdup=0;
+  fStrStream >> nup >> idprup >> xwgtup >> scalup >> aqedup >> aqcdup;
+  if (!fStrStream)
+    AliFatal("malformed event record (1st line)");
+
+  // todo: xwgtup, scalup, aqedup, aqcdup -> event header
+  fPosTracksBegin = fStrStream.tellg();
+  return nup;
+}
+
+TParticle* AliGenReaderLHE::NextParticle()
+{
+  if (!fStrStream) {
+    AliError("missing particle");
+    return NULL;
+  }
+  Int_t idup=0, istup=0, mothup[2]={0,0}, icolup[2]={0,0};
+  Double_t p[4] = { 0,0,0,0 }, mass=0, vtimup=0, spinup=0;
+  fStrStream >> idup >> istup >> mothup[0] >> mothup[1] >> icolup[0] >> icolup[1]
+	     >> p[0] >> p[1] >> p[2] >> p[3] >> mass >> vtimup >> spinup;
+  if (!fStrStream) {
+    AliError("malformed particle record");
+    return NULL;
+  }
+
+  const TLorentzVector vp(p);
+  fParticle.SetPdgCode(idup);
+  fParticle.SetStatusCode(istup);
+  fParticle.SetFirstMother(mothup[0]);
+  fParticle.SetLastMother (mothup[1]);
+  fParticle.SetMomentum(vp);
+  fParticle.SetProductionVertex(mass ? vp*(vtimup/mass) : TLorentzVector(0,0,0,0)); // TODO: check unit
+  fParticle.SetBit(kTransportBit, istup == 1);
+
+  return &fParticle;
+}
+
+void AliGenReaderLHE::RewindEvent()
+{
+  fStrStream.seekg(fPosTracksBegin);
+}

EVGEN/AliGenReaderLHE.h

@@ -0,0 +1,51 @@
+// -*- C++ -*-
+#ifndef ALIGENREADER_LHE_H
+#define ALIGENREADER_LHE_H
+
+// Les Houches record format reader
+// see, e.g.,
+//  * https://arxiv.org/pdf/hep-ph/0609017.pdf
+//  * https://arxiv.org/pdf/hep-ph/0109068.pdf
+
+#include <sstream>
+
+#include <TXMLEngine.h>
+#include <TParticle.h>
+
+#include "AliGenReader.h"
+
+class AliGenReaderLHE : public AliGenReader
+{
+public:
+  AliGenReaderLHE()
+    : AliGenReader()
+    , fXMLEngine()
+    , fXMLDoc(NULL)
+    , fXMLChild(NULL)
+    , fStrStream()
+    , fPosTracksBegin()
+    , fParticle() {}
+
+  virtual ~AliGenReaderLHE();
+
+  virtual void       Init();
+  virtual Int_t      NextEvent();
+  virtual TParticle* NextParticle();
+  virtual void       RewindEvent();
+
+protected:
+private:
+  AliGenReaderLHE(const AliGenReaderLHE&);
+  AliGenReaderLHE& operator=(const AliGenReaderLHE&);
+
+  TXMLEngine        fXMLEngine;      //!
+  XMLDocPointer_t   fXMLDoc;         //!
+  XMLNodePointer_t  fXMLChild;       //!
+  std::stringstream fStrStream;      //! contains event record
+  std::streampos    fPosTracksBegin; //! used in RewindEvent() method
+  TParticle         fParticle;       //!
+
+  ClassDef(AliGenReaderLHE, 1);
+} ;
+
+#endif // ALIGENREADER_LHE_H

EVGEN/CMakeLists.txt

@@ -99,6 +99,7 @@ set(SRCS
     AliGenPerformance.cxx
     AliGenReadersEMD.cxx
     AliGenToy.cxx
+    AliGenReaderLHE.cxx
 )
 
 # Headers from sources

EVGEN/EVGENLinkDef.h

@@ -77,5 +77,6 @@
 #pragma link C++ class  AliGenPerformance+;
 #pragma link C++ class  AliGenReadersEMD+;
 #pragma link C++ class  AliGenToy+;
+#pragma link C++ class  AliGenReaderLHE+;
 
 #endif