---

reconRootWriterAlg.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"

#include "GlastEvent/TopLevel/Event.h"

#include "GlastSvc/GlastDetSvc/IGlastDetSvc.h"

#include "reconRootWriterAlg.h"

#include "reconRootData/Recon.h"

#include "reconRootData/TkrSiCluster.h"
#include "reconRootData/TkrTrack.h"
#include "reconRootData/TkrRecon.h"
#include "reconRootData/ReconHeader.h"

#include "reconRootData/TkrLocator.h"
#include "reconRootData/TkrGamma.h"
#include "reconRootData/TkrHit.h"

#include "reconRootData/CalLogEne.h"
#include "reconRootData/CalCluster.h"
#include "reconRootData/CalRecon.h"

#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"
#include "TDirectory.h"

#include "TkrRecon/SiRecObjs.h"
#include "TkrRecon/SiClusters.h"
#include "CalRecon/CalRecLogs.h"
#include "CalRecon/CsIClusters.h"

static const AlgFactory<reconRootWriterAlg>  Factory;
const IAlgFactory& reconRootWriterAlgFactory = Factory;

static const int bufsize = 64000;
static const int split = 1;          // split objects into sub branches


reconRootWriterAlg::reconRootWriterAlg(const std::string& name, ISvcLocator* pSvcLocator) : 
Algorithm(name, pSvcLocator)
{
    declareProperty("reconRootFile",m_nameFile);
    declareProperty("inputSource",m_inputSource="SIM");
}



StatusCode reconRootWriterAlg::initialize()
{
    
    
    StatusCode sc = StatusCode::SUCCESS;
    MsgStream log(msgSvc(), name());
    
    // Use the Job options service to set the Algorithm's parameters
    setProperties();
    
    if(0 == gROOT)
    {
        static TROOT root("root","Root I/O");
        log << MSG::INFO << "Initializing ROOT" << endreq;
    } else {
        log << MSG::INFO << "Root already initialized(not an error)" <<endreq;
    }
    
    TDirectory *saveDir = gDirectory;
    
    recFile = new TFile(m_nameFile.c_str(), "RECREATE");
    if (!recFile->IsOpen()) sc = StatusCode::FAILURE;
    recFile->cd();
    recTree = new TTree("T", "TreeBase");
    m_rec = new Recon();
    m_rec->Create();
    recTree->Branch("Recon","Recon",&m_rec, bufsize, split);
    
    saveDir->cd();
    return sc;
    
}

StatusCode reconRootWriterAlg::execute()
{
    StatusCode sc = StatusCode::SUCCESS;
    
    sc = retrieve();
    
    saveHeader();

    save(m_SiClusters);
    save(m_SiRecObjs);
    save(m_CalRecLogs);
    save(m_CsIClusterList);
    
    writeEvent();
    return sc;
}

StatusCode reconRootWriterAlg::saveHeader() {
    MsgStream log(msgSvc(), name());
    // the location of the header info in the TDS depends upon
    // where the digi data is coming from...since when reading
    // from the Root files via the digiRootReaderAlg, we do not
    // have a proper Event object as the Root of the TDS.
    std::string headerPathStr;
    if (strcmp(m_inputSource.c_str(),"ROOT")==0) 
        headerPathStr = "/Event/Header";
    else 
        headerPathStr = "/Event";

    SmartDataPtr<Event> event(eventSvc(),headerPathStr.c_str());
    if( 0==event) { log << MSG::ERROR << "could not find the event header" << endreq;
       return StatusCode::FAILURE;
    }

    m_rec->getReconFlags()->setEventId(event->event());
    m_rec->getReconFlags()->setRunId(event->run());
    return StatusCode::SUCCESS;
}


void reconRootWriterAlg::writeEvent() 
{
    TDirectory *saveDir = gDirectory;
    recFile->cd();
    recTree->Fill();
    
    m_rec->Clean();
    m_rec->Create();
    
    saveDir->cd();
}

void reconRootWriterAlg::close() 
{
    TDirectory *saveDir = gDirectory;
    recFile->cd();
    recFile->Write(0, TObject::kOverwrite);
    recFile->Close();
    saveDir->cd();
}

StatusCode reconRootWriterAlg::finalize()
{
    close();
    
    StatusCode sc = StatusCode::SUCCESS;
    return sc;
}

StatusCode reconRootWriterAlg::retrieve()
{
    StatusCode sc = StatusCode::SUCCESS;

    m_SiClusters = SmartDataPtr<SiClusters>(eventSvc(),"/Event/TkrRecon/SiClusters"); 
    m_SiRecObjs  = SmartDataPtr<SiRecObjs>(eventSvc(),"/Event/TkrRecon/SiRecObjs");
    m_CalRecLogs  = SmartDataPtr<CalRecLogs>(eventSvc(),"/Event/CalRecon/CalRecLogs");
    m_CsIClusterList = SmartDataPtr<CsIClusterList>(eventSvc(),"/Event/CalRecon/CsIClusterList"); 

    return sc;
}

void reconRootWriterAlg::save(SiClusters* sicl)
{
    if (sicl)
    {
        int nhits = sicl->nHits();
        unsigned int ihit = 0;
        // pointer to the SiCluster TObjArray
        TObjArray *clusterArr = m_rec->getTkrRecon()->getSiClusters();
        // loop over all hits 
        for (ihit = 0; ihit < nhits; ihit ++) {
            SiCluster*    scl = sicl->getHit(ihit);
            TkrSiCluster* rcl = converter(scl);
            clusterArr->Add(rcl);
        }
    }
}

TkrSiCluster* reconRootWriterAlg::converter(const SiCluster* scl)
{
    TkrSiCluster* rcl = new TkrSiCluster();
    rcl->setCenterStrip(scl->strip());
    rcl->setId(scl->id());
    rcl->setLayer(scl->plane());
    rcl->setNumStrips(scl->size());
    Point p = scl->position();
    Float_t x = (scl->v() == 0? p.x():p.y());
    rcl->setPosition(x);
    TkrSiCluster::TKRAxes axis;
    axis = (scl->v()==0? TkrSiCluster::X : TkrSiCluster::Y);
    rcl->setXY(axis);
    rcl->setZPosition(p.z());
    return rcl;
}
//---------------- Tracks --------------------
void reconRootWriterAlg::save(SiRecObjs* sirec)
{       
    if (sirec)
    {
        importGamma(sirec);
        importTracks(sirec);
    }
}

void reconRootWriterAlg::importGamma(SiRecObjs* sirec)
{
    
    int ngamma = sirec->numGammas();
    if (ngamma !=1 ) return;
    
    // pointer to TrkGamma
        TkrGamma* rgamma = new TkrGamma(1);
        m_rec->getTkrRecon()->setGamma(rgamma);
        
        GFgamma* gamma = sirec->Gamma(0);
        TkrLocator* rloc = converter(gamma);
        ((TObjArray*) rgamma->getLocator())->Add(rloc);
        TObjArray* rgammaTracks = rgamma->getTracks();
        
        TObjArray* rtracks = m_rec->getTkrRecon()->getTracks();
        
        m_idtrack = 0;
        
        const GFtrack* xbest = gamma->getBest(SiCluster::X);
        TkrTrack* rbxtrack = converter(xbest,0,true);
        rbxtrack->setId((UShort_t) m_idtrack++);
        rtracks->Add(rbxtrack);
        rgammaTracks->Add(rbxtrack);
        
        const GFtrack* ybest = gamma->getBest(SiCluster::Y);
        TkrTrack* rbytrack = converter(ybest,0,true);
        rbytrack->setId( (UShort_t) m_idtrack++);
        rtracks->Add(rbytrack);
        rgammaTracks->Add(rbytrack);
        
        const GFtrack* xpair = gamma->getPair(SiCluster::X);
        if (!xpair->empty()) {
            TkrTrack* rpxtrack = converter(xpair,1,true);
            rpxtrack->setId( (UShort_t) m_idtrack++);
            rtracks->Add(rpxtrack);
            rgammaTracks->Add(rpxtrack);
        }
        
        const GFtrack* ypair = gamma->getPair(SiCluster::Y);
        if (!ypair->empty()) {
            TkrTrack* rpytrack = converter(ypair,1,true);
            rpytrack->setId( (UShort_t) m_idtrack++);
            rtracks->Add(rpytrack);
            rgammaTracks->Add(rpytrack);
        }
        
        if (rtracks->GetEntries() <=0) {
            int badbad = 1;
        }
}


void reconRootWriterAlg::importTracks(SiRecObjs* sirec)
{
    int nparticles = sirec->numParticles();
    if (nparticles <=0) return;
    
    TObjArray* rTracks = m_rec->getTkrRecon()->getTracks();
    for (int iparticle = 0; iparticle < nparticles; iparticle++) {
        GFparticle* particle = sirec->Particle(iparticle);
        const GFtrack* xtrack = particle->getXGFtrack();
        TkrTrack* rxtrack = converter(xtrack,iparticle);
        rxtrack->setId( (UShort_t) m_idtrack++);
        rTracks->Add(rxtrack);
        const GFtrack* ytrack = particle->getYGFtrack();
        TkrTrack* rytrack = converter(ytrack,iparticle);
        rytrack->setId( (UShort_t) m_idtrack++);
        rTracks->Add(rytrack);
    }
}


void reconRootWriterAlg::importHits(const GFtrack* track, TkrTrack* rtrack)
{
    int nhits = track->nhits();
    if (nhits <=2) return;
    
    SiCluster::view a = track->getAxis();
    TObjArray* rHits = rtrack->getHits();
    for (int ihit = 0; ihit < nhits; ihit++) {
        const KalPlane& kplane = track->kplanelist[ihit];
        TkrHit* rhit = converter(kplane,a);
        rHits->Add(rhit);
    }
    if (rHits->GetEntries() <=0 ) {
        int badbad = 1;
    }
}
//------------------ converters ---------------
TkrLocator* reconRootWriterAlg::converter(const GFdata* gf)
{
    
    TkrLocator* rloc = new TkrLocator(1);
    Float_t x = gf->vertex().x();
    Float_t y = gf->vertex().y();
    Float_t z = gf->vertex().z();
    rloc->setPosition(new TVector3(x,y,z));
    
    Float_t tx = gf->direction().x();
    Float_t ty = gf->direction().y();
    Float_t tz = gf->direction().z();
    Float_t slopex  = (tz != 0. ? slopex = tx/tz : 0);
    Float_t slopey  = (tz != 0. ? slopey = ty/tz : 0);
    rloc->setDirection(new TVector2(slopex,slopey));
    return rloc;
}

TkrTrack* reconRootWriterAlg::converter(const GFtrack* track, int id, bool mother)
{
    TkrTrack* rtrack = new TkrTrack(1);
    
    // set the locator at the vertex
    //TkrLocator* rloc = new TkrLocator(1);
    TkrLocator *rloc = converter(track);
    rtrack->getLocator()->Add(rloc);
    
    // set the other information
    Float_t chi2 = track->chiSquare();
    rtrack->setChi2(chi2);
    
    Float_t eneDet = track->RCenergy();
    rtrack->setEnergyDet(eneDet);
    
    Float_t eneIn  = track->inputEnergy();
    rtrack->setEnergyInput(eneIn);
    
    UInt_t fstLayer = track->firstLayer();
    rtrack->setFirstLayer(fstLayer);
    
    rtrack->setId(id);
    
    if (mother) rtrack->setMother(TkrTrack::GAMMA);
    else rtrack->setMother(TkrTrack::NONE);
    
    Float_t quality = track->Q();
    rtrack->setQuality(quality);
    
    Float_t residual = 0.;
    rtrack->setResidual(residual);
    
    rtrack->setTrackType(TkrTrack::ELECTRON);
    
    importHits(track, rtrack);
    
    return rtrack;
}

TkrHit* reconRootWriterAlg::converter(const KalPlane& kp, SiCluster::view a)
{
    TkrHit* rHit = new TkrHit();
    
    Float_t chi2 = kp.getDeltaChiSq(KalHit::SMOOTH);
    rHit->setChi2(chi2);
    
    Float_t xmea = kp.getHit(KalHit::MEAS).getPar().getPosition();
    Float_t xfit = kp.getHit(KalHit::SMOOTH).getPar().getPosition();
    Float_t res = xfit-xmea;
    rHit->setResidual(res);
    
    TkrLocator* rloc = converter(kp,a,1);
    rHit->setLocator(rloc);
    
    int ihit = kp.getIDHit();
    TkrSiCluster* rcl = m_rec->getTkrRecon()->getSiCluster(ihit);
    rHit->setSiCluster(rcl);
    
    return rHit;
}


TkrLocator* reconRootWriterAlg::converter(const KalPlane& kp, SiCluster::view a, int idummy)
{
    TkrLocator* rloc = new TkrLocator(1);
    Float_t x = 0;
    Float_t y = 0;
    Float_t xp = kp.getHit(KalHit::SMOOTH).getPar().getPosition();
    if (a == SiCluster::X) x = xp;
    else y = xp;
    Float_t z = kp.getZPlane();
    rloc->setPosition(new TVector3(x,y,z));
    
    Float_t slopex = 0.;
    Float_t slopey = 0.;
    Float_t slopep = kp.getHit(KalHit::SMOOTH).getPar().getSlope();
    if (a == SiCluster::X) slopex = slopep;
    else slopey = slopep;
    rloc->setDirection(new TVector2(slopex,slopey));
    
    Float_t ex = sqrt(kp.getHit(KalHit::SMOOTH).getCov().getsiga());
    Float_t es = sqrt(kp.getHit(KalHit::SMOOTH).getCov().getsigb());
    Float_t ec = kp.getHit(KalHit::SMOOTH).getCov().getcovab();
    
    if (a == SiCluster::X) rloc->setXmatrix(ex,es,ec);
    else rloc->setYmatrix(ex,es,ec);
    
    return rloc;
}



//--------------------------------------------------------------


void reconRootWriterAlg::save(CsIClusterList *ccl)
{
    if (ccl)
    {
        int nc = ccl->num();
        TObjArray* clusArr = m_rec->getCalRecon()->getCalClusters();
        for (int ic = 0; ic<nc; ic++){
            ICsICluster* csiclus = ccl->Cluster(ic);
            CalCluster* calclus = converter(csiclus);
            calclus->setId(ic);
            clusArr->Add(calclus);
        }
    }
}

CalCluster* reconRootWriterAlg::converter(const ICsICluster *cc)
{
    CalCluster* calclus = new CalCluster();     
    calclus->setEsum(cc->energySum());
    calclus->setEcorr(cc->energyCorrected());
    calclus->setEleak(cc->energyLeak());
    Point pos = cc->position();
    calclus->setPosition(pos.x(),pos.y(),pos.z());
    Vector dir = cc->direction();
    calclus->setDirection(dir.x(), dir.y(),dir.z());
    calclus->setCsiAlpha(cc->getCsiAlpha());
    calclus->setCsiLambda(cc->getCsiLambda());
    calclus->setCsiStart(cc->getCsiStart());
    calclus->setFitEnergy(cc->getFitEnergy());
    calclus->setProfChisq(cc->getProfChisq());
    for(int i=0;i<8;i++){
        calclus->setEneLayer(i,cc->getEneLayer(i));
        const Vector& posl = cc->getPosLayer(i);
        calclus->setXLayer(i,posl.x());
        calclus->setYLayer(i,posl.y());
    }
    // for the moment we have only one cluster and I include
    // all hited logs into its log list
    // this should be modified when real clusterisation
    // will be added
    TObjArray* logArr = m_rec->getCalRecon()->getCalLogs();
    TObjArray* cluslogs = calclus->getCalLog();
    TObjArrayIter it(logArr);
    TObject* log;
    while(log = it.Next())cluslogs->Add(log);
    
    return calclus;
}

void reconRootWriterAlg::save(CalRecLogs *crls)
{
    if (crls)
    {
        int nlogs = crls->num();
        TObjArray* logArr = m_rec->getCalRecon()->getCalLogs();
        
        for (int il=0; il<nlogs; il++){
            CalRecLog* crl = crls->Log(il);
            if(crl->posEnergy()>0 || crl->negEnergy()>0){
                CalLogEne* cle = converter(crl);
                logArr->Add(cle);
            }
        }
    }
}


CalLogEne* reconRootWriterAlg::converter(const CalRecLog* crl)
{
    CalLogEne* cle = new CalLogEne();
    cle->setePlus(crl->posEnergy());
    cle->seteMinus(crl->negEnergy());
    int tb_layer = crl->layer();
    int tb_col = crl->column();
    int tb_view = crl->view();
    int root_layer = 7-(2*tb_layer+tb_view);
    int root_col = tb_col;
    LogId* id = new LogId();
    id->setColumn(root_col);
    id->setLayer(root_layer);
    cle->setId(id);
    return cle; 
}

---