MJitterByCoincidence.cc

Go to the documentation of this file.

#include "MJitterByCoincidence.hh"
#include "QBaseType.hh"
#include "QVector.hh"
#include "QError.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QEventList.hh"
#include "QCoincidenceData.hh"
#include "QDetChannelCollectionHandle.hh"
#include "QJitter.hh"
#include "QJitterHandle.hh"
#include <vector>
#include <math.h>
#include <algorithm>
#include <TH1D.h>
 
using namespace std;
 
REGISTER_MODULE(MJitterByCoincidence)
 
using std::vector;
using namespace Diana;
 
void MJitterByCoincidence::Init(QEvent& ev)
{
  fCoincidenceDataOwner = GetString("CoincidenceDataOwner", "Coincidence");
  fOutput = GetString("Output", "jitter.txt");
  fRefCh_IsSet = false;
  fRefCh = GetInt("ReferenceChannel", Q_INT_DEFAULT, false);
  fTower = GetInt("Tower", Q_INT_DEFAULT, false);
  fJitterByDelaySource = GetString("JitterByDelaySource",Q_STRING_DEFAULT, false);
 
  Info("Tower is: %d", fTower);
 
  if (fRefCh!=Q_INT_DEFAULT) {
    Info("Reference Channel set to: %d", fRefCh);
    fRefCh_IsSet = true;
  }
 
  fEnergyLabel = GetString("EnergyLabel","ApplyCalibration@Energy");
  ev.RequireByLabel<QDouble>(fEnergyLabel);
 
  ev.Require<QPulseInfo>("DAQ","PulseInfo");
  ev.Require<QHeader>("DAQ","Header");
  ev.Require<QCoincidenceData>(fCoincidenceDataOwner,"CoincidenceData");
  ev.Require<QDouble>(fCoincidenceDataOwner,"TotalEnergy");
 
  // clear quantities
  fJitterSum = 0.;
  fJitterErrorSum = 0.;
  fDeltaT.clear();
}
 
//void MJitterByCoincidence::Do(QEvent& ev, const QEventList& neighbours)
void MJitterByCoincidence::Do(QEvent& ev)
{
    const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const int& channel = pulseInfo.GetChannelId();
 
    fGoodChannels.insert(channel);
 
    const QCoincidenceData& cd = ev.Get<QCoincidenceData>(fCoincidenceDataOwner.c_str(),"CoincidenceData");
    
    const QHeader& hdr = ev.Get<QHeader>("DAQ","Header");
    fRun = hdr.GetRun();
    
    QDetChannelCollectionHandle dccHandle;
    dccHandle.SetRun(fRun);
    GlobalData().Get("",&dccHandle,"");
    const QDetChannelCollection& DCC = dccHandle.Get();
    const std::vector<int>& boloRelChan = DCC.GetRelativeChannels(channel);
   
    if(cd.fMultiplicity==2) {
    // select multiplicity 2 event and order the coincident channels according to channel number
    const int o_channel = cd.fCoincidentChannels[0].fChannelId;
 
        if(channel>o_channel) {        
            const int deltaTower = cd.fCoincidentChannels[0].fDeltaTower;
            const int deltaFloor = cd.fCoincidentChannels[0].fDeltaFloor;
            const int deltaPos = cd.fCoincidentChannels[0].fDeltaPosition;
            
            // if coincident channel is not on same bolometer (relative channel)
            if (find(boloRelChan.begin(),boloRelChan.end(),o_channel)==boloRelChan.end()
                // if coincident channel is in same tower 
                && deltaTower==0
                // if coincident channel is in same floor OR in same column (+1 floor and same position)
                && ((abs(deltaFloor)==1 && deltaPos==0)
                || deltaFloor==0))
                {
                    const double deltaT = cd.fCoincidentChannels[0].fDeltaT;
                    // fills the map fDeltaT with a vector of time differences for each couple of channels
                    fDeltaT[channel][o_channel].push_back(deltaT);
                }
        }
    }
}
 
void MJitterByCoincidence::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 vecDeltaT;
    QVector vecDeltaT_red;
 
    double DeltaT = 0.;
    double DeltaTRMS = 0.;
 
    std::map<int,std::map<int,std::vector<double> > >::const_iterator iter;
    std::map<int,std::vector<double> > ::const_iterator mapIter;
 
    // calculates the mean value of vector fDeltaT for each couple of channels
    // mean value is calculated after removing the outliers (cut on median +- 3*MAD)
 
    iter = fDeltaT.begin();
    while(iter != fDeltaT.end()) {
        int ch = iter->first;
        mapIter = (iter->second).begin();
        while(mapIter != (iter->second).end()) {
            int o_ch = mapIter->first;
            if(fDeltaT[ch][o_ch].size()>0) {
                vecDeltaT.Resize(fDeltaT[ch][o_ch].size());
                vecDeltaT_red.Resize(0);
                for (size_t i=0; i<vecDeltaT.Size(); i++) {
                    vecDeltaT[i] = (fDeltaT[ch][o_ch])[i];
                }
                double medianDeltaT = vecDeltaT.GetMedian();
                double MADDeltaT = vecDeltaT.GetMedianAbsoluteDeviation();
                for (size_t j=0; j<vecDeltaT.Size(); j++) {
                    if(vecDeltaT[j]>(medianDeltaT-3.*MADDeltaT) && vecDeltaT[j]<(medianDeltaT+3.*MADDeltaT)) {
                        vecDeltaT_red.Append(vecDeltaT[j]);
                    }
                }
                if(vecDeltaT_red.Size()>0) {
                    DeltaT = vecDeltaT_red.Sum(vecDeltaT_red.Size(),0)/vecDeltaT_red.Size();
                    DeltaTRMS = vecDeltaT_red.GetRMS(vecDeltaT_red.Size())/sqrt(vecDeltaT_red.Size());      
                    fJitter[ch][o_ch] = DeltaT;
                    fJitterRMS[ch][o_ch] = DeltaTRMS;
                    // debug printing
                    std::cout<< ch << " " << o_ch <<" "<<fJitter[ch][o_ch]<<" "<< fJitterRMS[ch][o_ch] << std::endl;
                }
            }
            mapIter++;
        }
        iter++;
    }
 
    // loop over detector geometry
 
    const std::vector<int>& towers = DCC.GetTowers();
    fLastFloor = false;
  
    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) {
            // process by tower
            Info("Skip tower %d",tower); 
            continue;
        } else if (floors.size()==0) {
            // problem if size = 0
            Info("Tower %d has no floors",tower); 
            continue;
        }
      
        int last_floor = floors[floors.size()-1];
 
        // identify reference channel as the first valid channel in the first floor of the current tower
        //     |
        //   3 | 2
        // ----|-----
        //   4 | 1
        //     |
      
        int position=1;     
        int flfl = 1;
        while (!fRefCh_IsSet) {
            fRefCh = FindChanWithPosInTowerFloor(position,tower,flfl);
            std::map<int,std::map<int,double > >::const_iterator riter = fJitter.find(fRefCh); 
            if (riter!=fJitter.end()) {
                if(fJitterByDelaySource!=Q_STRING_DEFAULT) {
                    // get jitter by delay handle for this channel
                    QJitterHandle jdHandle(fRefCh,fDataset);
                    GlobalData().Get("JitterByDelay",&jdHandle,fJitterByDelaySource);
      
                    if(jdHandle.IsValid()) {
                        fRefCh_IsSet=true;
                    } else {
                        Warn("Channel %d, has no jitter by delay: %s",fRefCh,jdHandle.GetError().GetDescription().c_str());
                        position = position+1;
                    } 
                } else {
                    fRefCh_IsSet=true;
                }
            } else {
                position = position+1;
            }
            if (position>4) {
                Warn("Tower %d has no valid reference channels in floor 1",tower);
                position = 1;
                flfl++;
            }
            if( flfl > last_floor )
                Panic("Tower %d has no valid reference channels anywhere", tower);
        }
        Info("Reference Channel set to: %d", fRefCh);
     
        position = DCC.Get(fRefCh).fTd.fTdPos;
 
        // starts loop over floors
        for (size_t i=0; i<floors.size(); i++) {
            int floor = floors[i];
            if (floor == last_floor) 
                fLastFloor = true;
 
            if( floor < flfl )
                continue;
 
            Info("Processing floor: %d, position: %d",floor,position);
            position = ProcessFloor(position,tower,floor);
        }
    }
}
 
double MJitterByCoincidence::FindJitter(int ch, int o_ch)
{
  double jitter = 0.;
  if (ch>o_ch) 
    {
      std::map<int,std::map<int,double > >::const_iterator miter = fJitter.find(ch); 
      if (miter!=fJitter.end()) jitter = fJitter[ch][o_ch];
    }
  else if (ch<o_ch) 
    {
      std::map<int,std::map<int,double > >::const_iterator miter = fJitter.find(o_ch); 
      if (miter != fJitter.end()) {
std::map<int,double>::const_iterator viter = (miter->second).find(ch);
if (viter != (miter->second).end()) jitter = -fJitter[o_ch][ch];
      }
    }
  return jitter;
}
 
double MJitterByCoincidence::FindJitterRMS(int ch, int o_ch)
{
  double jitterRMS = 0.;
  if (ch>o_ch) 
    {
      std::map<int,std::map<int,double > >::const_iterator miter = fJitterRMS.find(ch); 
      if (miter!=fJitterRMS.end()) jitterRMS = fJitterRMS[ch][o_ch];
    }
  else if (ch<o_ch) 
    {
      std::map<int,std::map<int,double > >::const_iterator miter = fJitterRMS.find(o_ch); 
      if (miter != fJitterRMS.end()) {
std::map<int,double>::const_iterator viter = (miter->second).find(ch);
if (viter != (miter->second).end()) jitterRMS = fJitterRMS[o_ch][ch];
      }
    }
  return jitterRMS;
}
 
int MJitterByCoincidence::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]; 
    }
 
  return chan;
}
 
void MJitterByCoincidence::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);
 
    GlobalHandle<QString> clabelh("CalibrationLabel");
    clabelh.SetDataset(fDataset);
    clabelh.SetChannel(ch);
    GlobalData().Get(fEnergyLabel.owner,&clabelh,"CurrentReader");
    if(clabelh.IsValid()) {
      jitterHandle.SetCalibLabel(clabelh.Get());
    } else {
      jitterHandle.SetCalibLabel("");
    }
 
    jitterHandle.Set(jitterObj);
    GlobalData().Set(&jitterHandle,fOutput);
}
 
int MJitterByCoincidence::ProcessFloor(int start_pos, int tower, int floor)
{
 
  int ndet = 4;
  int pos[ndet];
  int ch_pos[ndet];
  int ch_pos_next_floor[ndet];
  int best_pos = 1; 
  double best_jitter = 0.;
  double best_jitterRMS = 0.;
 
  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);
      ch_pos_next_floor[i] = FindChanWithPosInTowerFloor(pos[i],tower,floor+1);
 
      double temp_jitterRMS = FindJitterRMS(ch_pos_next_floor[i],ch_pos[i]);
      if (temp_jitterRMS>0.) best_jitterRMS = temp_jitterRMS;
    }
 
   // finds path to next floor with lower error and calculates jitter (only if floor != last floor)
  if (!fLastFloor) {
     for (int i=0; i<ndet; i++)
      {
double floor_jitterRMS = FindJitterRMS(ch_pos_next_floor[i],ch_pos[i]);
double floor_jitter = FindJitter(ch_pos_next_floor[i],ch_pos[i]);
 
if(floor_jitterRMS!=0 && floor_jitterRMS<=best_jitterRMS) 
  {
    best_jitterRMS = floor_jitterRMS;
    best_jitter = floor_jitter;
    best_pos=pos[i];
  }
      }
  }
 
  // calculates jitter with respect to ref channel on this floor
 
  double jitter_rel[ndet];
  double jitterRMS_rel[ndet];
 
  double jitter11,jitter21,jitter41,jitter32,jitter34,jitter31;
  double jitter11RMS,jitter21RMS,jitter41RMS,jitter32RMS,jitter34RMS,jitter31RMS;
  double error321, error341;  
 
  jitter11 = FindJitter(ch_pos[0],ch_pos[0]);
  jitter11RMS = FindJitterRMS(ch_pos[0],ch_pos[0]);
 
  jitter21 = FindJitter(ch_pos[1],ch_pos[0]);
  jitter21RMS = FindJitterRMS(ch_pos[1],ch_pos[0]);
  
  jitter41 = FindJitter(ch_pos[3],ch_pos[0]);  
  jitter41RMS = FindJitterRMS(ch_pos[3],ch_pos[0]);
  
  jitter32 = FindJitter(ch_pos[2],ch_pos[1]);
  jitter32RMS = FindJitterRMS(ch_pos[2],ch_pos[1]);
  
  jitter34 = FindJitter(ch_pos[2],ch_pos[3]);  
  jitter34RMS = FindJitterRMS(ch_pos[2],ch_pos[3]);
  
  jitter31 = FindJitter(ch_pos[2],ch_pos[0]);  
  jitter31RMS = FindJitterRMS(ch_pos[2],ch_pos[0]);
 
  if (jitter32RMS!=0. && jitter21RMS!=0.) {error321 = sqrt(pow(jitter32RMS,2)+pow(jitter21RMS,2));}
  else error321 = 9999.;
  if (jitter34RMS!=0. && jitter41RMS!=0.) {error341 = sqrt(pow(jitter34RMS,2)+pow(jitter41RMS,2));}
  else error341 = 9999.;
 
  // position 3  (finds path with lower error)
 
  if (error321<error341) {
    jitter_rel[2] = fJitterSum + jitter32 + jitter21;
    jitterRMS_rel[2] = fJitterErrorSum + pow(error321,2);
    //    std::cout<< ch_pos[2]<<" "<<ch_pos[0]<<" "<<jitter_rel[2]<<" "<< sqrt(jitterRMS_rel[2])<<std::endl;
    SaveJitter(ch_pos[2],jitter_rel[2],sqrt(jitterRMS_rel[2]));
  }
  else if (error341<error321) {
    jitter_rel[2] = fJitterSum + jitter34 + jitter41;
    jitterRMS_rel[2] = fJitterErrorSum + pow(error341,2);
    //   std::cout<< ch_pos[2]<<" "<<ch_pos[0]<<" "<<jitter_rel[2]<<" "<< sqrt(jitterRMS_rel[2])<<std::endl;
    SaveJitter(ch_pos[2],jitter_rel[2],sqrt(jitterRMS_rel[2]));
  }
  else {
    if (jitter31!=0.) {  
      jitter_rel[2] = fJitterSum + jitter31;
      jitterRMS_rel[2] = fJitterErrorSum + pow(jitter31RMS,2);
      //      std::cout<< ch_pos[2]<<" "<<ch_pos[0]<<" "<<jitter_rel[2]<<" "<< sqrt(jitterRMS_rel[2])<<std::endl;
      SaveJitter(ch_pos[2],jitter_rel[2],sqrt(jitterRMS_rel[2]));
    }
  }
 
  // position 2
  if (jitter21!=0.) {  
    jitter_rel[1] = fJitterSum + jitter21;
    jitterRMS_rel[1] = fJitterErrorSum + pow(jitter21RMS,2);
    //   std::cout<< ch_pos[1]<<" "<<ch_pos[0]<<" "<<jitter_rel[1]<<" "<< sqrt(jitterRMS_rel[1])<<std::endl;
    SaveJitter(ch_pos[1],jitter_rel[1],sqrt(jitterRMS_rel[1])); 
 }
  
  // position 4
  if (jitter41!=0.) {  
    jitter_rel[3] = fJitterSum + jitter41;
    jitterRMS_rel[3] = fJitterErrorSum + pow(jitter41RMS,2);
    //    std::cout<< ch_pos[3]<<" "<<ch_pos[0]<<" "<<jitter_rel[3]<<" "<< sqrt(jitterRMS_rel[3])<<std::endl;
    SaveJitter(ch_pos[3],jitter_rel[3],sqrt(jitterRMS_rel[3]));   
   }
 
  // position 1
  jitter_rel[0] = fJitterSum + jitter11;
  jitterRMS_rel[0] = fJitterErrorSum + pow(jitter11RMS,2);
  //  std::cout<< ch_pos[0]<<" "<<ch_pos[0]<<" "<<jitter_rel[0]<<" "<< sqrt(jitterRMS_rel[0])<<std::endl;
  SaveJitter(ch_pos[0],jitter_rel[0],sqrt(jitterRMS_rel[0]));
 
  // calculates jitter of channel with best position (to go to lower floor) with respect to reference channel;
  // increments fJitterSum and fJitterErrorSum
 
  if (!fLastFloor) {
 
    if (best_pos!=start_pos)
      {
int ch = 0;
if (best_pos>start_pos){
  ch = best_pos - start_pos;
}
if (best_pos<start_pos){
  ch = 4 - (start_pos - best_pos);
}
 
fJitterSum = jitter_rel[ch] + best_jitter;
fJitterErrorSum = fJitterErrorSum + pow(jitterRMS_rel[ch],2) + pow(best_jitterRMS,2); 
      }
    else
      {
fJitterSum = fJitterSum + best_jitter; 
fJitterErrorSum = fJitterErrorSum + pow(best_jitterRMS,2); 
      }
  }
 
  // return position of channel with lower error on jitter
  return best_pos;
}