htofhitf2.cc

Go to the documentation of this file.

#include "htofhitf2.h"
#include "hades.h"
#include "htofraw.h"
#include "htofrawsim.h"
#include "htofhit.h"
#include "htofhitsim.h"
#include "htofcalpar.h"
#include "htofwalkpar.h"
#include "hruntimedb.h"
#include "hcategory.h"
#include "hiterator.h"
#include "hdebug.h"
#include "tofdef.h"
#include "hevent.h"
#include "heventheader.h"
#include "hspectrometer.h"
#include "htofdetector.h"
#include "hgeomvolume.h"
#include "hgeomcompositevolume.h"
#include "hgeomtransform.h"
#include "htofgeompar.h"
#include "htofdigipar.h"
#include "hdetgeompar.h"
#include "hgeomvector.h"
#include "hspecgeompar.h"
#include "hstart2hit.h"
#include "hstartdef.h"



#include "TMath.h"

using namespace std;


#include <fstream>
#include <iomanip>

ClassImp(HTofHitF2)

//*-- Author :
//*-- Modified : 4/2013 by O. Svoboda
//*-- Modified : 14/11/2004 by P. Tlusty -
//*   tof hit lab position now related to [0,0,0]
//*-- Modified : 23/10/2002 by M. Sanchez
//*-- Modified : 09/05/2002 by D. Zovinec
//*-- Modified : 14/03/2002 by D. Zovinec
//*-- Modified : 21/09/2001 by D. Zovinec
//*-- Modified : 24/04/2001 by M. Sanchez
//*-- Modified : 08/03/2001 by M. Sanchez
//*-- Modified : 27/03/2000 by I. Koenig
//*-- Modified : 17/03/2000 by R. Holzmann
//*-- Modified : 30/11/99 by M. Sanchez
//*-- Modified : 02/11/99 by D. Vasiliev
//*-- Modified : 15/06/98 by M. Sanchez


//_HADES_CLASS_DESCRIPTION 
/////////////////////////////////////////////////////////////////////////
// HTofHitF2 Tof hit finder
//
// adapted from /u/halo/packages/abase/new by Manuel Sanchez (17/06/98)
//
// HTofHitF2 is the reconstructor class that iterates over RAW data
// and finds the TOF hits.
// The HIT is only created and filled if both left and right TDC time
// is available.
//
// Left (right) TDC time is used to find xpos and tof of hit (see HTofHit).
// Then HIT coordinates in the LAB, phi, theta and distance from the
// target are calculated by using HTofGeomPar.
//
// Left (right) ADC amplitude is used to find xposAdc and energy
// deposited in the rod . Additional information amplL and amplR
// is calculated (see HTofHit).
// The flagAdc have been added to distinguish levels of ADC info
// available:  no ADC         - flagAdc = 0
//             left ADC only  - flagAdc = -1
//             right ADC only - flagAdc = 1
//             both ADC's     - flagAdc = 2
// Important!!!: Note that if only one ADC is available (case flagAdc = -1 or 1)
// the xposAdc cannot be calculated. Also amplL (amplR) cannot be found
// when left (right) ADC has not provided the signal.
//
// Empirical dependence of tof on position along the rod included
// (see talk of P.Tlusty in Coll.Meeting XI at Smolenice, TOF detector performance.)
//
// Time walk correction is included in the tof reconstruction.
//
/////////////////////////////////////////////////////////////////////////

Bool_t  HTofHitF2::bforceKeepGeant=kFALSE;

HTofHitF2::HTofHitF2(void)
{
    fCalPar     = NULL;
    fCalParSim  = NULL;
    fTofWalkPar = NULL;
    fRawCat = fHitCat = NULL;
    fRawCatTmp = fHitCatTmp = NULL;
    fStartHitCat   = NULL;
    fTofSimulation = kFALSE;
    fLoc.set(3,0,0,0);
    iter     = NULL;
    iterTmp  = NULL;
    fDigiPar = NULL;

}

HTofHitF2::HTofHitF2(const Text_t *name,const Text_t *title) : HReconstructor (name,title)
{
    fCalPar     = NULL;
    fCalParSim  = NULL;
    fTofWalkPar = NULL;
    fDigiPar    = NULL;
    fRawCat = fHitCat = NULL;
    fRawCatTmp = fHitCatTmp = NULL;
    fStartHitCat   = NULL;
    fTofSimulation = kFALSE;
    fLoc.set(3,0,0,0);
    iter    = NULL;
    iterTmp = NULL;

}

HTofHitF2::~HTofHitF2(void)
{
    if(iter)    delete iter;
    if(iterTmp) delete iterTmp;
}

void HTofHitF2::initParContainer(HSpectrometer *spec, HRuntimeDb *rtdb)
{
    // Adds to the runtime database "rtdb" the containers needed by the Hit Finder
    //spec is used to get information of the spectrometer setup.

    fCalPar       = (HTofCalPar *)  rtdb->getContainer("TofCalPar");
    fTofGeometry  = (HTofGeomPar *) rtdb->getContainer("TofGeomPar");
    fSpecGeometry = (HSpecGeomPar *)rtdb->getContainer("SpecGeomPar");
    fDigiPar      = (HTofDigiPar *) rtdb->getContainer("TofDigiPar");


    if(!fTofSimulation || ( fTofSimulation && gHades->getEmbeddingMode()>0) ){
        // walk correction not needed for simulation
        fTofWalkPar   = (HTofWalkPar *) rtdb->getContainer("TofWalkPar");
        fTofWalkPar->getContainers(); // need other container for calculation
    }
    return;

}

Bool_t HTofHitF2::init(void)
{
    // Initializes data levels and container parameters for operation of the
    //hit finder
    Bool_t r=kTRUE; // return value
    HSpectrometer *spec = gHades->getSetup();
    HRuntimeDb    *rtdb = gHades->getRuntimeDb();
    HEvent *ev          = gHades->getCurrentEvent(); // Event structure

    if (spec && rtdb && ev) {
        HTofDetector *tof =(HTofDetector*) spec->getDetector("Tof");

        if(ev->getCategory(catGeantKine)) fTofSimulation=kTRUE;

        if (tof) {
            // The real work starts here
            initParContainer(spec,rtdb);

            //-------------------------------------------------
            // input categories
            fRawCat = ev->getCategory(catTofRaw);
            if (!fRawCat) {
                if (!fRawCat) return kFALSE;
            }
            iter=(HIterator*)fRawCat->MakeIterator("native");

            // Get Start Hit category. If StartHit is not defined
            // a Warning is displayed and fStartHitCat is set to 0
            fStartHitCat = ev->getCategory(catStart2Hit);
            if (!fStartHitCat) Warning("init","Start hit level not defined; setting start time to 0");


            //-------------------------------------------------
            // output categories
            if(!fTofSimulation){ // real data
                fHitCat = ev->getCategory(catTofHit);
                if (!fHitCat) {
                    fHitCat = tof->buildCategory(catTofHit);
                    if (!fHitCat) return kFALSE;
                    else ev->addCategory(catTofHit,fHitCat,"Tof");
                }
            } else { // simulation or embedding
                fHitCat=ev->getCategory(catTofHit);
                if (!fHitCat) {
                    fHitCat=tof->buildMatrixCategory("HTofHitSim",0.5F);
                    if (!fHitCat) return kFALSE;
                    else ev->addCategory(catTofHit,fHitCat,"Tof");
                } else {
                    if (fHitCat->getClass()!=HTofHitSim::Class()) {
                        Error("HTofHitfSim::init()","Misconfigured output category");
                        return kFALSE;
                    }
                }
            }

            if(fTofSimulation && gHades->getEmbeddingMode()>0)
            {  // embedding mode
                fRawCatTmp =ev->getCategory(catTofRawTmp);
                if(!fRawCatTmp){
                    Error("init()","No catTofRawTmp in input! You are in embedding mode!");
                    return kFALSE;
                }
                iterTmp = (HIterator*)fRawCatTmp->MakeIterator();  // cannot be "native", because
                // iter->getLocation() is used

                fHitCatTmp=ev->getCategory(catTofHitTmp);
                if (!fHitCatTmp) {
                    fHitCatTmp=tof->buildMatrixCategory("HTofHitSimTmp",0.5F);
                    if (!fHitCatTmp) {
                        Error("init()","Could no build catTofHitTmp !");
                        return kFALSE;
                    }
                    else ev->addCategory(catTofHitTmp,fHitCatTmp,"Tof");
                    fHitCatTmp->setPersistency(kFALSE);
                }

                fCalParSim = new HTofCalPar("TofCalPar_Sim_Embedding","","TofCalProductionSimEmbedding");
                if(fCalParSim)rtdb->addContainer(fCalParSim);
                else {
                    Error("init()","Could no create HTofCalPar for sim params, needed for embedding !");
                }

            }
            //-------------------------------------------------



        } else {
            Error("init","TOF setup is not defined");
            r = kFALSE; // Notify error to task manager
        }

    } else {
        Error("init","Setup, RuntimeDb or event structure not found");
        r = kFALSE; // Notify error to task manager
    }

    return r;

}

Int_t HTofHitF2::execute(void)
{
    // Find the hits.
    //See also HReconstructor::execute
#if DEBUG_LEVEL>2
    gDebuger->enterFunc("HTofHitF2::execute");
#endif

    //---------------------------------------------
    // read catTofRaw and create hits in catTofHit
    Bool_t sim=kFALSE;
    if      ( fTofSimulation) sim=kTRUE;
    else if (!fTofSimulation) sim=kFALSE;

    if( !( gHades->getEmbeddingMode()>0 && gHades->getEmbeddingDebug() == 1) ) fillHitCat(sim,kFALSE); // in sim or real or embeding allways fill here, only embedding+debug skipped
    //---------------------------------------------

    //---------------------------------------------
    // In embedding mode one has to
    // read catTofRawTmp and create hits in catTofHitTmp
    // in addition. Embedded Hits will be merged on hitlevel
    // later
    if(fTofSimulation&&gHades->getEmbeddingMode()>0){
        fillHitCat(kTRUE,kTRUE);
        mergeHitCats(sim,kTRUE);
    }
    //---------------------------------------------

#if DEBUG_LEVEL>2
    gDebuger->leaveFunc("HTofHitF2::execute");
#endif
    return 0;
}

void HTofHitF2::fillGeometry(HTofHit *hit) {
    // Fills in the LAB position parameters for the given HTofHit.
    //
    // This is done by transforming the point (x,0,0) in the ROD coordinate
    // system to the LAB coordinate system. Where x is the reconstructed
    // x position inside the hit.

    HGeomVector rLocal,r;
    Float_t d,phi,theta;
    Float_t rad2deg = 180./TMath::Pi();

    HModGeomPar *module=fTofGeometry->getModule(hit->getSector(),hit->getModule());
    HGeomTransform &trans = module->getLabTransform();
    HGeomVolume *rodVol=module->getRefVolume()->getComponent(hit->getCell());
    HGeomTransform &rodTrans=rodVol->getTransform();

    // Fill r with the hit coordinates in the ROD coordinate system
    // Since we do not have information about y,z coordinates of impact
    // y=0, z=0 is used. Note that (0,0,0) corresponds to the rod center.
    r.setX(hit->getXpos());
    r.setY(0.);
    r.setZ(0.);
    rLocal=rodTrans.transFrom(r);  // transform to module coordinate system
    r=trans.transFrom(rLocal);     // transform from module to LAB system
    HGeomVolume *tv=fSpecGeometry->getTarget(0);
    if (tv) r -= tv->getTransform().getTransVector();   // correct for target position

    // Fill output
    d = r.length();
    theta = (d>0.) ? (rad2deg * TMath::ACos(r.getZ() / d)) : 0.;
    phi = rad2deg * TMath::ATan2( r.getY(), r.getX());
    if (phi < 0.) phi += 360.;
    if (tv) r += tv->getTransform().getTransVector();   // correct for target position

    hit->setXYZLab(r.getX(), r.getY(), r.getZ());
    hit->setDistance( d );
    hit->setTheta(theta);
    hit->setPhi(phi);
}

void HTofHitF2::fillHitCat(Bool_t sim,Bool_t embed)
{


    HTofRaw *raw=NULL;
    HTofHit *hit=NULL;
    HStart2Hit *sH=NULL;

    Float_t atof;
    Float_t axpos;
    Float_t startTime=0.0;
    Float_t startTimeSmearing=0.0;
    Float_t subCl=0.0;      // raw left charge subtracted by left ADC pedestal
    Float_t subCr=0.0;      // raw right charge subtracted by right ADC pedestal
    Float_t slopeAdcL;   // left ADC slope
    Float_t slopeAdcR;  // right ADC slope
    Float_t xposAdc;
    Float_t depE;
    Float_t leftA;
    Float_t rightA;
    Float_t twalk,twalk_pos;
    Int_t flagadc;

    //--------------------------------------------------
    // start time extraction. For embedding some workd has to
    // be done
    startTime = 0.0;
    startTimeSmearing = 0.0;  // for simulation/embedding of sim tracks  pStartH->getResolution()
    // startTime  = startTimeSmearing  for simulation
    // startTime != startTimeSmearing  for embedding
    // startTimeSmearing == -1000      for real data
    if (fStartHitCat && fStartHitCat->getEntries()>0) {
        if((sH = (HStart2Hit *) fStartHitCat->getObject(0)) != NULL){
            startTime = sH->getTime();
            if(sH->getResolution()!=-1000) startTimeSmearing = sH->getResolution();
        }
    }
    //--------------------------------------------------

    HIterator*  iterLocal;
    HCategory*  catHitLocal;
    HTofCalPar* calparLocal;

    //--------------------------------------------------
    //  Are we running on real or sim data?
    //
    //  sim  && !embed  && gHades->getEmbeddingMode()==0  sim   data  1. call
    //  sim  && !embed  && gHades->getEmbeddingMode() >0  real  data  1. call embedding
    //  sim  &&  embed      sim   data  2. call embedding
    // !sim  && !embed      real  data  1. call real data

    Bool_t isRealData = kFALSE;
    if((sim && !embed && gHades->getEmbeddingMode() >0) ||
       (!sim && !embed)) isRealData = kTRUE;
    //--------------------------------------------------

    //--------------------------------------------------
    // switch data pointer depending if running for
    // sim or real data!
    if(sim&&embed)  iterLocal=iterTmp;
    else            iterLocal=iter;

    if(sim&&embed)  catHitLocal=fHitCatTmp;
    else            catHitLocal=fHitCat;

    if(sim&&embed)  calparLocal=fCalParSim;
    else            calparLocal=fCalPar;
    //--------------------------------------------------

    Float_t startTimeLocal = startTime;        // real works normal
    if(!isRealData) startTimeLocal = startTimeSmearing;// embedded hist have to use the


    // smearing since startTime is from real data
    iterLocal->Reset();
    while ( (raw=(HTofRaw *)iterLocal->Next())!=NULL) {
        fLoc[0]=raw->getSector();
        fLoc[1]=raw->getModule();
        fLoc[2]=raw->getCell();

        // Hit level is only filled if both left and right time are set
        if(raw->getLeftTime() && raw->getRightTime()){
            hit = (HTofHit *)catHitLocal->getSlot(fLoc);
            if (hit) {

                if(fTofSimulation) hit=new (hit) HTofHitSim;
                else               hit=new (hit) HTofHit;

                HTofCalParCell& cell  =(*calparLocal)[fLoc[0]][fLoc[1]][fLoc[2]];

                // New calibration with 20 bins and two side walk correction
                // The starting point starts at X=0 -> center of the rod.
                // After 6 iteration value converged (usually already at the third
                // iteration the minimum is reached.
                if (isRealData) { // walk correction has to be perform
                    fTofWalkPar->getTofPos(raw, atof, axpos, 0., startTimeLocal);
                    Float_t postemp = axpos;
                    fTofWalkPar->getTofPos(raw, atof, axpos, postemp, startTimeLocal);
                    postemp = axpos;
                    fTofWalkPar->getTofPos(raw, atof, axpos, postemp, startTimeLocal);
                    postemp = axpos;
                    fTofWalkPar->getTofPos(raw, atof, axpos, postemp, startTimeLocal);
                    postemp = axpos;
                    fTofWalkPar->getTofPos(raw, atof, axpos, postemp, startTimeLocal);
                    postemp = axpos;
                    fTofWalkPar->getTofPos(raw, atof, axpos, postemp, startTimeLocal);
                } else {  // for sim hits no walk correction
                    // The TDC's data is elaborated here.
                    atof = (raw->getLeftTime() * cell.getLeftK() +
                            raw->getRightTime()*cell.getRightK())/2.0 - cell.getTimK();

                    //Substract start time. Zero if no start hit level is filled
                    atof -= startTimeLocal;

                    axpos = cell.getVGroup()*(raw->getRightTime() * cell.getRightK() -
                                              raw->getLeftTime()*cell.getLeftK())/2.0 +cell.getPosK();
                }

                // The ADC's data is elaborated here.
                xposAdc=0.0;
                depE=0.0;
                leftA=0.0;
                rightA=0.0;
                flagadc=0;
                //twalk=twoffT;
                twalk=0.;
                twalk_pos=0.;

                // If at least one ADC signal then preliminary calculations.
                if (!isRealData) {  // sim data
                    if (raw->getLeftCharge())  subCl = (raw->getLeftCharge() - cell.getPedestalL());    //original version
                    if (raw->getRightCharge())  subCr = (raw->getRightCharge() - cell.getPedestalR());   //original version
                }

                if (isRealData) { // real data
                    subCl = tot2amp(raw->getLeftCharge() - cell.getPedestalL());    //modified O.S.
                    subCr = tot2amp(raw->getRightCharge() - cell.getPedestalR());   //modified O.S.
                }
                slopeAdcL = (cell.getEdepK())*(TMath::Exp((cell.getGainAsym())/(cell.getAttLen())));
                slopeAdcR = (cell.getEdepK())*(TMath::Exp(-(cell.getGainAsym())/(cell.getAttLen())));

                leftA=subCl*slopeAdcL;       //changed 25.4.2013
                rightA=subCr*slopeAdcR;
                // Individual cases.
                if(subCl>0){
                    flagadc=-1;
                    depE=(subCl*cell.getEdepK())*(TMath::Exp((cell.getGainAsym()-axpos)/(cell.getAttLen())));
                    // leftA=subCl*slopeAdcL;
                    twalk=-(cell.getTimeWalkC1()/(TMath::Sqrt(leftA))); // time walk correction (left ADC)
                    twalk_pos=(cell.getTimeWalkC1()/(TMath::Sqrt(leftA))); // time walk correction for position (left ADC)
                }
                else {
                    twalk=-(cell.getTimeWalkC1()/(TMath::Sqrt(0.35)));
                    twalk_pos=(cell.getTimeWalkC1()/(TMath::Sqrt(0.35)));
                }
                if(subCr>0){
                    flagadc=1;
                    depE=(subCr*cell.getEdepK())*(TMath::Exp((axpos-cell.getGainAsym())/(cell.getAttLen())));
                    //rightA=subCr*slopeAdcR;
                    twalk=twalk-(cell.getTimeWalkC2()/(TMath::Sqrt(rightA))); // time walk correction (right ADC)
                    twalk_pos=twalk_pos-(cell.getTimeWalkC2()/(TMath::Sqrt(rightA))); // time walk correction for position (right ADC)
                }
                else {
                    twalk=twalk-(cell.getTimeWalkC2()/(TMath::Sqrt(0.35)));
                    twalk_pos=twalk_pos-(cell.getTimeWalkC2()/(TMath::Sqrt(0.35)));
                }
                if(subCl>0 && subCr>0) {
                    flagadc=2;
                    xposAdc=(cell.getAttLen()/2.0)*(TMath::Log(subCl/subCr)) + cell.getGainAsym();
                    depE=(cell.getEdepK())*(TMath::Sqrt(subCl*subCr));
                }


                hit->setSector((Char_t) fLoc[0]);
                hit->setModule((Char_t) fLoc[1]);
                hit->setCell((Char_t) fLoc[2]);
                hit->setTof(atof);
                hit->setXpos(axpos);
                hit->setXposAdc(xposAdc);
                hit->setEdep(depE);
                hit->setLeftAmp(leftA);
                hit->setRightAmp(rightA);
                hit->setAdcFlag(flagadc);

                if(fTofSimulation) {
                    HTofHitSim *hs=(HTofHitSim *)hit;
                    HTofRawSim *rs=(HTofRawSim *)raw;
                    hs->setNTrack1(rs->getNTrack1());
                    hs->setNTrack2(rs->getNTrack2());
                }

                fillGeometry(hit);
            }
        }
    }
}
void HTofHitF2::mergeHitCats(Bool_t sim,Bool_t embed)
{
    //--------------------------------------------------
    HTofHitSim* hit;
    HTofHitSim* hittmp;

    HTofCalPar* calparLocal;
    if(sim&&embed)  calparLocal=fCalParSim;
    else            calparLocal=fCalPar;

    TIterator* hititer=fHitCatTmp->MakeIterator("native");  // sim data in embedding
    hititer->Reset();

    while ( (hittmp=(HTofHitSim *)hititer->Next())!=NULL)
    {
        fLoc[0]=hittmp->getSector();
        fLoc[1]=hittmp->getModule();
        fLoc[2]=hittmp->getCell();
        hit= (HTofHitSim*) fHitCat->getObject(fLoc);

        if(hit)
        {   // real hit is existing already
            // merging should be done here. At the moment
            // sim hits will be always transported


            if(gHades->getEmbeddingMode()==1 && !bforceKeepGeant)
            {
                HTofCalParCell&  cell    =(*calparLocal)[fLoc[0]][fLoc[1]][fLoc[2]];
                HTofDigiParCell& celldigi=(*fDigiPar)   [fLoc[0]][fLoc[1]][fLoc[2]];

                Float_t rTime1 = (celldigi.getHalfLen() + hit->getXpos())/cell.getVGroup() + hit->getTof();
                Float_t lTime1 = (celldigi.getHalfLen() - hit->getXpos())/cell.getVGroup() + hit->getTof();

                Float_t rTime2 = (celldigi.getHalfLen() + hittmp->getXpos())/cell.getVGroup() + hittmp->getTof();
                Float_t lTime2 = (celldigi.getHalfLen() - hittmp->getXpos())/cell.getVGroup() + hittmp->getTof();

                Float_t  rTime = rTime1;
                Float_t  lTime = lTime1;

                if(rTime>rTime2) rTime = rTime2;
                if(lTime>lTime2) lTime = lTime2;

                Float_t tof  = (lTime + rTime)/2. - celldigi.getHalfLen()/cell.getVGroup();
                Float_t xpos = cell.getVGroup() * (rTime - lTime)/2.;


                hit->setTof(tof);
                hit->setXpos(xpos);

                fillGeometry(hit);  // fill lab pos again

                hit->setEdep    (hit->getEdep()    + hittmp->getEdep());
                hit->setLeftAmp (hit->getLeftAmp() + hittmp->getLeftAmp());
                hit->setRightAmp(hit->getRightAmp()+ hittmp->getRightAmp());

                if(hittmp->getNTrack1()==hittmp->getNTrack2()){ // first track id is real track , second geant id
                    hit->setNTrack2(hittmp->getNTrack1());
                } else { // if we have 2 different numbers we want to keep it: drop real track id
                    hit->setNTrack1(hittmp->getNTrack1());
                    hit->setNTrack2(hittmp->getNTrack2());
                }

            }
            else if(gHades->getEmbeddingMode() == 2 || (gHades->getEmbeddingMode()==1 && bforceKeepGeant) ) {  // let GEANT particle survive
                new (hit) HTofHitSim(*hittmp);
            } else {
                Error("mergeHitCats()","Unknow embedding mode = %i",gHades->getEmbeddingMode()) ;
            }
        }
        else
        {   // cell was not fired by real hit before
            hit= (HTofHitSim*) fHitCat->getSlot(fLoc);
            if(hit)
            {
                new (hit)HTofHitSim(*hittmp);
            }
            else
            {
                Error("mergeHitCats()","Could not retrieve slot in catTofHit!");
            }
        }

    }

    delete hititer;
}


Float_t HTofHitF2::tot2amp(Float_t tot)
{
    Float_t amp = -1.;
    if(tot<150) amp =  3.0 * tot;
    else amp =  128. + 2.14 * tot;
    return amp;
}