---

XANAMcInfo.cpp

Go to the documentation of this file.

#include "Utilities/Configuration/interface/Architecture.h"

#include <Utilities/Notification/interface/Observer.h>
#include "Utilities/Notification/interface/PackageInitializer.h"
#include <Utilities/UI/interface/SimpleConfigurable.h> 

#include "CARF/SimEvent/interface/SimEvent.h"
#include "GeneratorInterface/HepEvent/interface/RawHepEvent.h"
#include "GeneratorInterface/EventProxyInterface/interface/HepEventG3EventProxyReader.h"
#include "GeneratorInterface/Particle/interface/RawParticle.h"
#include "GeneratorInterface/Particle/interface/RawHepEventParticle.h"
#include "GeneratorInterface/HepPDT/interface/HepPDTable.h"
#include "GeneratorInterface/HepPDT/interface/HepParticleData.h"

#include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantVertex.h>
#include <XANADOO/XANAMcInfo/interface/XANAGeantTrack.h>
#include <XANADOO/XANAUtilities/interface/XANAToolkit.h>


class XANAMcInfoTester : private Observer<G3EventProxy *> {

 public:
 
  SimpleConfigurable<bool> readPileup_;
  SimpleConfigurable<int> minBunch_;
  SimpleConfigurable<int> maxBunch_;
  SimpleConfigurable<bool> readVertices_;
  
  XANAMcInfoTester()
   {    
    xana::out(1) << "[update] Constructing an XANAMcInfoTester" << endl ;
    Observer<G3EventProxy *>::init() ; 
    readPileup_=SimpleConfigurable<bool>(false,"HepG3EventProxyReader:ReadPileup"); 
    minBunch_=SimpleConfigurable<int>(0,"HepG3EventProxyReader:MinBunch"); 
    maxBunch_=SimpleConfigurable<int>(0,"HepG3EventProxyReader:MaxBunch"); 
    readVertices_=SimpleConfigurable<bool>(false,"HepG3EventProxyReader:ReadVertices");    
   }
   
  ~XANAMcInfoTester()
   {
   }

  void upDate( G3EventProxy * ev )
   {
    
    if (!ev) return ; 
    
    // first generator info
    vector<XANAGeneratorParticle> particles;
   
    const RawHepEventParticle * aHepParticle;
    HepEventG3EventProxyReader EventReader(ev);
    RawHepEvent MyHepEvent;
    MyHepEvent.read(&EventReader);
    
    xana::out(4) << std::endl;
    xana::out(4) << "[update] --- RawHepEvent ---" << std::endl;
    MyHepEvent.print();
    xana::out(4) << std::endl << "[update] --- total event tree ---" << std::endl;
    MyHepEvent.printTree();
    xana::out(4) << std::endl;
   
    XANAGeneratorEvent *aMcEvent = new XANAGeneratorEvent;
    XANAGeneratorParticle *aMcParticle;
    // !! in RawHepEvent, i runs from 1 to n() !!
    for(unsigned int i=1; i<=MyHepEvent.n(); i++) {
      aHepParticle = MyHepEvent.getParticle(i);
      HepLorentzVector momentum(aHepParticle->px(), aHepParticle->py(),
       aHepParticle->pz(), aHepParticle->e());
      aMcParticle = new XANAGeneratorParticle(aHepParticle->line(), aHepParticle->status(),
       aHepParticle->pid(), aHepParticle->mother1(), aHepParticle->mother2(),
       aHepParticle->daughter1(), aHepParticle->daughter2(), momentum,
       aHepParticle->mass(), aHepParticle->vertex());
      // assuming signal particles are the first ones
      if (i<=MyHepEvent.nSignal()) {
        aMcEvent->addSignalParticle(aMcParticle);
      } else {
        aMcEvent->addPileupParticle(aMcParticle);
      }
    }

    xana::out(4) << std::endl << "[update] --- XANAGeneratorEvent ---" << std::endl;
    aMcEvent->print();
    xana::out(4) << std::endl;
    
    // then Geant info

    xana::out(4) << "[update] --- GEANT vertices -----------------------------------" << std::endl;
    SimEvent::vertex_container GeantVertices; // this is an STL vector

    GeantVertices = ev->simSignal()->proxy()->vertices(); // get the vertex container
    xana::out(4) << "[update] This event has " << GeantVertices.size() << " GEANT vertices." <<std::endl;    
    xana::out(4).setf(std::ios::fixed, std::ios::floatfield);
    xana::out(2) << "[update] The annihilation process took place at " 
         << "x=" << std::setw(7) << std::setprecision(4) << GeantVertices[0].position().x()
         << ", y=" << std::setw(7) << std::setprecision(4) << GeantVertices[0].position().y()
         << ", z=" << std::setw(7) << std::setprecision(4) << GeantVertices[0].position().z() 
         << std::endl;

    xana::out(4) << "[update] --- GEANT tracks -----------------------------------" << std::endl;

    SimEvent::track_container GeantTracks; // this is an STL vector

    GeantTracks = ev->simSignal()->proxy()->tracks();
    xana::out(4) << "[update] This event has " << GeantTracks.size() << " GEANT tracks." <<std::endl;
    xana::out(2) << "[update] The first track is:" << std::endl;
    xana::out(2) << GeantTracks[0] << std::endl;

    SimEvent::genpart_container GenParts; 
    GenParts = ev->simSignal()->proxy()->genparts();
    if (! GeantTracks[0].noGenpart()) {
      int genpart = GeantTracks[0].genpartIndex();
      xana::out(2) << "[update] This corresponds to generator particle ";
      xana::out(2) << genpart << std::endl;
      GenParts[genpart].print();
    }
    if (! GeantTracks[0].noVertex()) {
      int vertex =  GeantTracks[0].vertIndex();
      xana::out(2) << "[update] This track comes from the Geant Vertex ";
      xana::out(2) << vertex << " at " << std::endl;
      xana::out(2) << "x=" << std::setw(7) << std::setprecision(4) << GeantVertices[vertex].position().x()
           << ", y=" << std::setw(7) << std::setprecision(4) << GeantVertices[vertex].position().y()
           << ", z=" << std::setw(7) << std::setprecision(4) << GeantVertices[vertex].position().z() 
           << std::endl;
    }

    XANAGeantEvent aGeantEvent;
    aGeantEvent.addGeneratorEvent(aMcEvent);
    XANAGeantVertex *vtx=0;
    XANAGeantTrack *trk=0;
    int vertex=0;
    
    for(unsigned int it=0; it<GeantTracks.size(); it++) {
      xana::out(2) << "[update] Track:" << it << " is " << std::endl;
      xana::out(2) << GeantTracks[it] << std::endl;
      // decode mother-daughter relationship
      xana::out(2) << "[update] track momentum " << GeantTracks[it].momentum() << std::endl;
      HepLorentzVector mom(GeantTracks[it].momentum());
      int mothindex = 0 ;
      if (mom.e() < 0.) {
        mothindex =  (int)(-GeantTracks[it].momentum().e()) ;
        xana::out(2) << "[update] track has daughter " << mothindex << std::endl;
        double energy2 = mom.perp2();
        double mass = HepPDT::theTable().getParticleData(GeantTracks[it].type())->mass();
        energy2 += mass*mass ;
        double energy = sqrt(energy2) ;
        mom.setE(energy);
      }
      xana::out(2) << "[update] recomputed geant track momentum " << mom << std::endl;
      // now create XANA track
      trk = new XANAGeantTrack(it+1,GeantTracks[it].type(), mom, mothindex);
      if (! GeantTracks[it].noVertex()) {
        vertex =  GeantTracks[it].vertIndex();
        xana::out(2) << "[update] This track comes from the Geant Vertex ";
        xana::out(2) << vertex << " at " << std::endl;
        xana::out(2) << "x=" << std::setw(7) << std::setprecision(4) << GeantVertices[vertex].position().x()
             << ", y=" << std::setw(7) << std::setprecision(4) << GeantVertices[vertex].position().y()
             << ", z=" << std::setw(7) << std::setprecision(4) << GeantVertices[vertex].position().z() 
             << endl;   
        vtx = new XANAGeantVertex(GeantVertices[vertex].position());
        aGeantEvent.addGeantVertex(vtx);
        aGeantEvent.addGeantTrack(trk, vtx);    
      } else
        aGeantEvent.addGeantTrack(trk);              
      
      xana::out(2) << "[update] Now setting link to generator particle"<< std::endl;
      XANAGeneratorParticle *parent;
      if (! GeantTracks[it].noGenpart()) {
        int genpart = GeantTracks[it].genpartIndex();
        xana::out(2) << "[update] This corresponds to generator particle ";
        xana::out(2) << genpart << std::endl;
        GenParts[genpart].print();
         // create generator particle with available info
        parent = new XANAGeneratorParticle(genpart+1, GenParts[genpart].status(),
         GenParts[genpart].pid(), GenParts[genpart].mother1(), GenParts[genpart].mother2(),
         GenParts[genpart].daughter1(), GenParts[genpart].daughter2(), 
         GenParts[genpart].fourmomentum(), 
         0., GeantVertices[vertex].position());
        aGeantEvent.setGeneratorParticle(trk, parent);              
      }

    }
    xana::out(4) << std::endl;
    xana::out(4) << "[update] Printing XANAGeantEvent "<< std::endl;
    aGeantEvent.print();
    xana::out(4) << std::endl;
    xana::out(4) << "[update] Related generator event is "<< std::endl;
    aGeantEvent.getGeneratorEvent()->print();
    xana::out(4) << std::endl;     
 
  }
   
};

static PKBuilder<XANAMcInfoTester> mybuilder("XANAMcInfoTester");

---