MJitterByDelay.cc

Go to the documentation of this file.

#include "MJitterByDelay.hh"
#include "QBaseType.hh"
#include "QVector.hh"
#include "QError.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QEventList.hh"
#include "QOFData.hh"
#include "QDetChannelCollectionHandle.hh"
#include "QJitter.hh"
#include "QJitterHandle.hh"
#include <vector>
#include <math.h>
#include <algorithm>
#include <TH1D.h>
 
REGISTER_MODULE(MJitterByDelay)
 
        using std::vector;
        using namespace Diana;
 
void MJitterByDelay::Init(QEvent& ev)
{
        fOutput = GetString("Output", "jitter.txt");
        fRefCh_IsSet = false;
        fRefCh = GetInt("ReferenceChannel", Q_INT_DEFAULT, false);
        fTower = GetInt("Tower", Q_INT_DEFAULT, false);
 
        Info("Tower is: %d", fTower);
 
        if (fRefCh!=Q_INT_DEFAULT) {
                Info("Reference Channel set to: %d", fRefCh);
                fRefCh_IsSet = true;
        }
 
        fOFLabel = GetString("OFLabel","COF"); // to maintain backwards compatibility
 
        ev.Require<QPulseInfo>("DAQ","PulseInfo");
        ev.Require<QHeader>("DAQ","Header");
        ev.Require<QOFData>(fOFLabel,"OFData");
 
        // clear quantities
        fJitterSum = 0.;
        fJitterErrorSum = 0.;
        fDelays.clear();
}
 
//void MJitterByDelay::Do(QEvent& ev, const QEventList& neighbours)
void MJitterByDelay::Do(QEvent& ev)
{
        const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo");
        const int& channel = pulseInfo.GetChannelId();
 
        fGoodChannels.insert(channel);
 
        const QOFData& ofd = ev.Get<QOFData>(fOFLabel.c_str(),"OFData");
        const double& delay = ofd.GetDelay();
 
        const QHeader& hdr = ev.Get<QHeader>("DAQ","Header");
        fRun = hdr.GetRun();
 
        // creates QVector of delays for each channel, if already set, return
        if(fDelays.find(channel) != fDelays.end()) 
        {     
                (fDelays.find(channel)->second).Append(delay);
                return;
        }
 
        QVector vecDelay;
        vecDelay.Append(delay);
        fDelays[channel] = vecDelay;
        vecDelay.Clear();
}
 
void MJitterByDelay::Done()
{
        // get dataset
        GlobalHandle<QInt> dHandle("Dataset");
        GlobalData().Get("",&dHandle,"");
        fDataset = dHandle.Get();
 
        QDetChannelCollectionHandle dccHandle;
        dccHandle.SetRun(fRun);
        GlobalData().Get("",&dccHandle,"");     
        const QDetChannelCollection& DCC = dccHandle.Get();
 
        QVector Delays;
        QVector Delays_red;
 
        double Delay = 0.;
        double DelayRMS = 0.;
 
        std::map<int,Diana::QVector > ::const_iterator mapIter;
        for (mapIter=fDelays.begin(); mapIter!=fDelays.end(); ++mapIter)
        {
                int ch = mapIter->first;
                Delays = mapIter->second;
 
                // calculates the mean value of vector fDelays for each channel
                // mean value is calculated after removing the outliers (cut on median +- 3*MAD)
 
                double medianDelays = Delays.GetMedian();
                double MADDelays = Delays.GetMedianAbsoluteDeviation();
 
                for (size_t j=0; j<Delays.Size(); j++)
                {
                        if(Delays[j]>(medianDelays-3.*MADDelays) && Delays[j]<(medianDelays+3.*MADDelays))
                        {
                                Delays_red.Append(Delays[j]);
                        }
                }
                if(Delays_red.Size()>0) {
 
                        Delay = (Delays_red.Sum(Delays_red.Size(),0))/((double) Delays_red.Size());
                        DelayRMS = (Delays_red.GetRMS(Delays_red.Size()))/(sqrt((double) Delays_red.Size()));         
                        fDelay[ch]=Delay;
                        fDelayRMS[ch]=DelayRMS; 
                }
 
                Delays_red.Clear();
                Delays.Clear();
                std::cout << ch << " " << fDelay[ch] << std::endl;
        }
 
        // loop over detector geometry
 
        const std::vector<int>& towers = DCC.GetTowers();
        fLastFloor = false;
 
        // if reference channel is not specified, identify reference channel as the channel in position 1, 
        // first floor, first tower
        //     |
        //   3 | 2
        // ----|-----
        //   4 | 1
        //     |
 
        // if reference channel is set find position of reference channel
 
 
 
        //FIXME temporary workaround, must check if there are cases where it doesn't work. See debugMJitterByDelay.txt
        int position = 0;  
 
        if (!fRefCh_IsSet) {
                for(position = 1; position <= 4; position++)    
                {
                        fRefCh = FindChanWithPosInTowerFloor(position,fTower,1);
                        Info("Reference Channel set to: %d, Position: %d", fRefCh,position);
                        if(fRefCh!=Q_INT_DEFAULT) continue;
                }
        }
 
        else if (fRefCh_IsSet) {
                position = DCC.Get(fRefCh).fTd.fTdPos;
        }
 
        for (size_t j=0; j<towers.size(); j++)
        {
                int tower = towers[j];
                const std::vector<int>& floors = DCC.GetFloorsInTower(tower);
 
                // check if tower has floors
                if (floors.size()>0) Info("Processing tower: %d",tower);
                if (tower!=fTower) {Info("Skip tower %d",tower); continue;} // process by tower
                else if (floors.size()==0) {Info("Tower %d has no floors",tower); continue;} //problem if size =0
 
                int last_floor = floors[floors.size()-1];
 
                // starts loop over floors
                for (size_t i=0; i<floors.size(); i++)
                {
                        int floor = floors[i];
                        if (floor == last_floor) fLastFloor = true;
 
                        Info("Processing floor: %d, position: %d",floor,position);
                        position = ProcessFloor(position,tower,floor);
                }
        }
}
 
double MJitterByDelay::FindJitter(int ch, int o_ch)
{
        double jitter = 0.;
        std::map<int,double>::const_iterator miter = fDelay.find(ch); 
        std::map<int,double>::const_iterator viter = fDelay.find(o_ch); 
        if (miter!=fDelay.end() && viter!=fDelay.end()) 
        {
                //     jitter = fDelay[ch]-fDelay[o_ch];
                jitter = fDelay[ch];
        }
        return jitter;
}
 
double MJitterByDelay::FindJitterRMS(int ch, int o_ch)
{
        double jitterRMS = 0.;
        std::map<int,double>::const_iterator miter = fDelayRMS.find(ch); 
        std::map<int,double>::const_iterator viter = fDelayRMS.find(o_ch); 
        if (miter!=fDelayRMS.end() && viter!=fDelayRMS.end()) 
        {
                //      jitterRMS = fDelayRMS[ch]-fDelayRMS[o_ch];
                jitterRMS = fDelayRMS[ch];
        }
        return jitterRMS;
}
 
int MJitterByDelay::FindChanWithPosInTowerFloor(int pos, int tower, int floor)
{
        int chan = Q_INT_DEFAULT;
        std::vector<int> chans;
 
        QDetChannelCollectionHandle dccHandle;
        dccHandle.SetRun(fRun);
        GlobalData().Get("",&dccHandle,"");     
        const QDetChannelCollection& DCC = dccHandle.Get();
 
        const std::vector<int>& chTF = DCC.GetChannelsInTowerFloor(tower,floor);
 
        for (size_t i=0; i<chTF.size(); i++)
        {
                if (DCC.Get(chTF[i]).fTd.fTdPos == pos) chans.push_back(chTF[i]); 
        }
 
        for (size_t j=0; j<chans.size(); j++)
        {
                if (fGoodChannels.count(chans[j])>0) 
                {
                        chan = chans[j];
                        std::cout << "good channels : " << chans[j] << std::endl; 
                }
        }
 
        return chan;
}
 
void MJitterByDelay::SaveJitter(int ch, double jitter, double jitter_error)
{
        QJitter jitterObj;
        jitterObj.Clear();  
 
        // set the object jitterObj
        jitterObj.fRefChannel = fRefCh; 
        jitterObj.fJitter = jitter; 
        jitterObj.fJitterError = jitter_error; 
 
        QJitterHandle jitterHandle(ch,fDataset);
        jitterHandle.SetCalibVersion(GetConfig().fVersionTag);
        jitterHandle.SetCalibLabel("");
 
        jitterHandle.Set(jitterObj);
        GlobalData().Set(&jitterHandle,fOutput);
}
 
int MJitterByDelay::ProcessFloor(int start_pos, int tower, int floor)
{
        int ndet = 4;
        int pos[ndet];
        int ch_pos[ndet];
        int best_pos = 1; 
 
        double jitter_rel;
        double jitterRMS_rel;
 
        for (int i=0; i<ndet; i++)
        {
                pos[i] = start_pos+i;
                if (pos[i]>4) pos[i]=pos[i]-4;
 
                ch_pos[i] = FindChanWithPosInTowerFloor(pos[i],tower,floor);
                jitter_rel = FindJitter(ch_pos[i],fRefCh);
                jitterRMS_rel = FindJitterRMS(ch_pos[i],fRefCh);
                if (jitter_rel!=0.) SaveJitter(ch_pos[i],jitter_rel,jitterRMS_rel);     
        }
        return best_pos;
}