MCoincidence.cc

Go to the documentation of this file.

#include "MCoincidence.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QEventList.hh"
#include "QBaseType.hh"
#include "QCoincidenceData.hh"
#include "QRunDataHandle.hh"
#include "QJitter.hh"
#include "QJitterHandle.hh"
#include <cmath>
#include <cstdlib>
#include <string>
 
using std::vector;
using std::string;
using std::ifstream;
using namespace Diana;
 
REGISTER_MODULE(MCoincidence)
 
 
void MCoincidence::Init(QEvent& ev)
{
    // called before event loop
    fCoincidenceWindow = GetDouble("CoincidenceWindow", 0.100);
    fCoincidenceRadius = GetDouble("CoincidenceRadius", 1200.);
    ev.Add<QCoincidenceData>("CoincidenceData");
    fEnergyLabel = GetString("EnergyLabel","ApplyCalibration@Energy");
    vector<string> defaults;
    fPSALabels = GetVectorString("PSALabels",defaults);
    std::string timeVariableStr = GetString("TimeVariable","Raw");
    if(timeVariableStr == "Raw") {
        fTimeVariable = TV_RAW;
    } else if(timeVariableStr == "OF") {
        fTimeVariable = TV_OF;
    } else {
        Warn("Unrecognized time variable, setting it to Raw");
        fTimeVariable = TV_RAW;
    }
 
    // get dataset
    GlobalHandle<QInt> dHandle("Dataset");
    GlobalData().Get("",&dHandle,"");
    fDataset = dHandle.Get();
 
    fSyncTime = GetBool("SyncTime",false);
    if(fSyncTime) {
        fSyncType = GetString("SyncType","coincidence");
        if (fSyncType == "coincidence")
        {
            fSyncSource = GetString("SyncSource","jitter.txt");
            fSyncTowerSource = GetString("SyncTowerSource","jitter_tower.txt");
        }
        else if (fSyncType == "delay")
        {
            fSyncTowerSource = GetString("SyncTowerSource","jitter_tower.txt");
        }
        else {
            Warn("Unrecognized SyncType variable, setting it to default (coincidence)");
            fSyncType = "coincidence";
        }
    }
    fOkInvalidTotalEnergy = GetBool("OkInvalidTotalEnergy",false,false);
    ev.Add<QDouble>("TotalEnergy");
    ev.Add<QDouble>("Time");
    ev.Add<QDouble>("Radius"); // DEBUG
    fUseRunningWindow = GetBool("UseRunningWindow",true);
    fStorePosition = GetBool("StorePosition",true,false);
    fThisRun = -1;
 
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Require<QHeader>("DAQ","Header");
    if(!fPSALabels.empty())  {
        for(size_t p = 0; p < fPSALabels.size(); p++) {
            ev.RequireByLabel<QDouble>(fPSALabels[p]);
        }
    }
    ev.RequireByLabel<QDouble>(fEnergyLabel);
 
    fOFLabel = GetString("OFLabel","COF"); // to maintain backwards compatibility
    if(fTimeVariable == TV_OF) ev.Require<QDouble>(fOFLabel,"Jitter");
 
}
 
void MCoincidence::Do(QEvent& ev, const QEventList& n_neighbours)
{
    const QHeader& hdr = ev.Get<QHeader>("DAQ","Header");
    const int run = hdr.GetRun();
   // const int chan = pi.GetChannelId();
 
    if(fThisRun != run) {
        fThisRun = run;
        QRunDataHandle rHandle(fThisRun);
        GlobalData().Get("DAQ",&rHandle,"");
        fRunData = rHandle.Get();
    }
 
 
    // Get the energy of the main event
    bool validTotalEnergy = true;
    double tot_ene;
    GetEventEnergy(ev,tot_ene,validTotalEnergy);
 
    QCoincidenceData& coincData = ev.Get<QCoincidenceData>("CoincidenceData");
    coincData.fOrderInMultiple = 1;
    coincData.fMultiplicity = 1;
    coincData.fCoincidentChannels.clear();
    double LargestDeltaT = 0;
    double SmallestDeltaT = Q_DOUBLE_DEFAULT;
    const double time = GetTime(ev);// depends on the cfg options what "time" is used
    ev.Get<QDouble>("Time") = time;
 
    // order neighbours according to our time variable;
    std::map<double,int> indices; // map of <time,index>
    for(size_t i = 0; i < n_neighbours.Size(); i++) {
        double o_time = GetTime(n_neighbours[i]);
        if(indices.find(o_time) != indices.end()) o_time+=(i+1)*1.e-09;
        indices[o_time] = i;
    }
    QEventList neighbours;
    double radius = 0.;
    std::map<double,int>::const_iterator index = indices.begin();
    while(index != indices.end()) {
        // implement here radial cut
        // will be possible to compute maximum radius of events after time selection is performed (i.e. leftIndex and rightIndex are defined)
        radius = GetDistanceMm(ev, &n_neighbours[index->second]);
        if( radius <= fCoincidenceRadius )neighbours.Push(&n_neighbours[index->second]);
        index++;
    }
 
    // create left and right index of the neighbours list == perform time selection
    int leftIndex = 1, rightIndex = -1;
    double l_time = time, r_time = time;
    for(int i = int(neighbours.Size())-1; i >= 0; i--) {// inverse loop on neigh.
        const double o_time = GetTime(neighbours[i]);
        const double deltaT = o_time - l_time;
        if(deltaT < 0) {
            if(fabs(deltaT) < fCoincidenceWindow) {
                if(fUseRunningWindow) 
                    l_time = GetTime(neighbours[i]);
                leftIndex = i;
            }
        } else if(fabs(deltaT) < fCoincidenceWindow) { 
            leftIndex = i;
        }
    }
    for(size_t i = leftIndex; i < neighbours.Size(); i++) {
        const double o_time = GetTime(neighbours[i]);
        const double deltaT = o_time - r_time;
        if(deltaT > 0) {
            if(fabs(deltaT) < fCoincidenceWindow) {
                if(fUseRunningWindow)
                    r_time = GetTime(neighbours[i]);
                rightIndex = i;
            }
        } else if(fabs(deltaT) < fCoincidenceWindow) { 
            rightIndex = i;
        }
    }
    // end of time selection
 
    int tower = -1;
    int floor = -1;
    int position = -1;
    if(fStorePosition) 
    {
//        tower = dccHandle.Get().Get(chan).fTd.fTdTower;
//        floor = dccHandle.Get().Get(chan).fTd.fTdFloor;
//        position = dccHandle.Get().Get(chan).fTd.fTdPos;
    }
 
 
    // fill informations
    double maxradius = 0.;
    for(int i = leftIndex; i <= rightIndex; i++) {
        const QEvent& other = neighbours[i];
        const QPulseInfo& o_pi = other.Get<QPulseInfo>("","PulseInfo");
        const double o_time = GetTime(neighbours[i]);
        const int o_channel = o_pi.GetChannelId();
        const double deltaT = o_time - time;
 
        double o_energy;
        bool validEnergy;
        GetEventEnergy(other,o_energy,validEnergy);
        validTotalEnergy&=validEnergy;
 
        int o_tower = -1; 
        int o_floor = -1; 
        int o_position = -1; 
        if(fStorePosition) 
        {
//            o_tower = dccHandle.Get().Get(o_channel).fTd.fTdTower;
//            o_floor = dccHandle.Get().Get(o_channel).fTd.fTdFloor;
//            o_position = dccHandle.Get().Get(o_channel).fTd.fTdPos;
        }
        if(deltaT < 0.) coincData.fOrderInMultiple++;
        // this is to get unique ordering when channels have same time.
        if(deltaT == 0. && ev.GetReadNumber() > other.GetReadNumber()) {
            coincData.fOrderInMultiple++;
        }
 
        coincData.fMultiplicity++;
        if(validEnergy)
            tot_ene += o_energy;
 
        // compute maximum distance btw main && coincident channels FIXME <-- maybe energy center vs {coincidents + main}?
        double o_radius = GetDistanceMm(ev,&other); // dist(Main <-> Coincident)
        if(o_radius > maxradius)maxradius = o_radius;
 
        if(fabs(deltaT) > fabs(LargestDeltaT))
            LargestDeltaT = deltaT;
        if(fabs(deltaT) < fabs(SmallestDeltaT))
            SmallestDeltaT = deltaT;
        QCoincidentChannel cc;
        cc.fChannelId = o_channel;
        cc.fDeltaT = deltaT;
        cc.fEnergy = o_energy;
        for(size_t p = 0; p < fPSALabels.size(); p++) {
            cc.fPSAVariables.push_back(other.GetByLabel<QDouble>(fPSALabels[p]));
        }
        cc.fDeltaTower = o_tower-tower;
        cc.fDeltaFloor = o_floor-floor;
        cc.fDeltaPosition = o_position-position;
        cc.fEventNumber = other.Get<QHeader>("DAQ","Header").GetEventNumber();
        coincData.fCoincidentChannels.push_back(cc);
    }
    ev.Get<QDouble>("Radius") = maxradius; // write radius into event
 
    coincData.fNearestCoincidentIndex = -1;
    coincData.fFarthestCoincidentIndex = -1;
    for(size_t c = 0; c < coincData.fCoincidentChannels.size(); c++) {
        if(coincData.fCoincidentChannels[c].fDeltaT == SmallestDeltaT) 
            coincData.fNearestCoincidentIndex = c;
        if(coincData.fCoincidentChannels[c].fDeltaT == LargestDeltaT) 
            coincData.fFarthestCoincidentIndex = c;
    }
    if(validTotalEnergy)
        ev.Get<QDouble>("TotalEnergy") = tot_ene;
}
 
void MCoincidence::Done()
{
    // called at the end of the event loop
}
 
double MCoincidence::GetTime(const QEvent& ev)
{
    const QHeader& header = ev.Get<QHeader>("DAQ","Header");
    const QPulseInfo& pi = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    double time = double(header.GetTime().GetFromStartRunNs())/1e09; 
    double preTriggerTime = 0;
 
    switch(fTimeVariable) {
        case TV_RAW:
            {       
                break;
            }
        case TV_OF:
            {
                const QDouble& ofjitter = ev.Get<QDouble>(fOFLabel.c_str(),"Jitter");
                time += ofjitter/1000;
                break;
            }
 
    }
 
    // response sinchronisation using by default jitter by coincidence and jitter by delay only for the reference channels
    if(fSyncTime && fSyncType == "coincidence") {
        const int channel = pi.GetChannelId();
 
        // get jitter handle for this channel
        QJitterHandle jHandle(channel,fDataset);
        GlobalData().Get("JitterByCoincidence",&jHandle,fSyncSource);
 
        if(!jHandle.IsValid()) {
            if(fBlackList.count(channel)) return time;
            fBlackList.insert(channel);   
            Error("Channel %d, has no jitter. Cannot be synced: %s",channel,jHandle.GetError().GetDescription().c_str());
            return time;
        }
        else {
            const QJitter& jitters = jHandle.Get(); 
 
            // get jitter for the reference channel of the given channel
            QJitterHandle jdhandle(jitters.fRefChannel,fDataset);
            GlobalData().Get("JitterByDelay",&jdhandle,fSyncTowerSource);
 
            if(!jdhandle.IsValid()) {
                Error("Reference channel %d, has no jitter. The relative tower cannot be synced: %s",jitters.fRefChannel,jdhandle.GetError().GetDescription().c_str());
                return time;
            }
 
            const QJitter& jittersTower = jdhandle.Get(); 
            time = time + jitters.fJitter - jittersTower.fJitter/1000. + preTriggerTime;
        }
    }
 
    // response sinchronisation using only jitter by delay 
    if(fSyncTime && fSyncType == "delay") {
        const int channel = pi.GetChannelId();
 
        // get jitter handle for this channel
        QJitterHandle jHandle(channel,fDataset);
        GlobalData().Get("JitterByDelay",&jHandle,fSyncTowerSource);
 
        if(!jHandle.IsValid()) {
            if(fBlackList.count(channel)) return time;
            fBlackList.insert(channel);   
            Error("Channel %d, has no jitter. Cannot be synced: %s",channel,jHandle.GetError().GetDescription().c_str());
            return time;
        }
        else {
            const QJitter& jitters = jHandle.Get(); 
            time = time - jitters.fJitter/1000. + preTriggerTime;
        }
    }
 
    return time;
}
 
void MCoincidence::GetEventEnergy(const Diana::QEvent &ev,double &energy,bool &valid) const {
 
    Diana::QEventLabel branchToSelect = fEnergyLabel;
    Diana::ReadHandle<Diana::QDouble> evHandle(branchToSelect.GetName().c_str());
 
    ev.Get(branchToSelect.GetOwner().c_str(),evHandle);
    valid=true;
    energy = Q_DOUBLE_DEFAULT;
 
    // If its Ok to have invalid total energy, only access the event
    // energy if its valid.
    if(fOkInvalidTotalEnergy) {
        if(evHandle.IsValid()) {
            energy = evHandle.Get();
        }
        else
            valid = false;
    }
    else 
        // If its not ok to have invalid total energy, then this will
        // throw if the value is invalid
        energy = evHandle.Get();
}
 
double MCoincidence::GetDistanceMm(const Diana::QEvent &ev,const Diana::QEvent* coincident){
    // get channel ids
    const int ch1 = ev.Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
    const int ch2 = coincident->Get<QPulseInfo>("DAQ","PulseInfo").GetChannelId();
    /* compute distance
       double x1,y1,z1,x2,y2,z2;
       QChannelUtils::ChannelPosition(ch1,x1,y1,z1);
       QChannelUtils::ChannelPosition(ch2,x2,y2,z2);
       if(!check)Error("Could not load position of channel %d and %d.",ch1,ch2);
       return sqrt( pow(x1-x2,2) + pow(y1-y2,2) +  pow(z1-z2,2) );
       */
    double x1,y1,z1,x2,y2,z2;
    x1 = fChannelPositionX[ch1];
    y1 = fChannelPositionY[ch1];
    z1 = fChannelPositionZ[ch1];
    x2 = fChannelPositionX[ch2];
    y2 = fChannelPositionY[ch2];
    z2 = fChannelPositionZ[ch2];
    return sqrt( pow(x1-x2,2) + pow(y1-y2,2) +  pow(z1-z2,2) );     
}