MFastCoincidence.cc

Go to the documentation of this file.

#include "MFastCoincidence.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QEventList.hh"
#include "QBaseType.hh"
#include "QCoincidenceData.hh"
#include "QOFData.hh"
#include "QRunDataHandle.hh"
#include "QDetChannelCollectionHandle.hh"
#include "QJitter.hh"
#include "QJitterHandle.hh"
 
#include <cmath>
#include <fstream>
#include <cstdlib>
#include <string>
#include <iostream>
#include <iomanip>
 
using namespace std;
using namespace Diana;
 
REGISTER_MODULE(MFastCoincidence);
 
inline bool FileExists(const std::string& name) { 
    if(FILE* file = fopen(name.c_str(), "r")) {
        fclose(file);
        return true;
    }
    else
        return false;
}
 
TimeNeighbour::TimeNeighbour(ULong64_t eventNumber, ULong64_t tag, int channel, double time, double energy) :
    fEventNumber(eventNumber), fTag(tag), fChannel(channel), fTime(time), fEnergy(energy) {
    }
 
bool timeSort(TimeNeighbour* a, TimeNeighbour* b) {
    return a->fTime < b->fTime;
}
 
void MFastCoincidence::Init(QEvent& ev) {
    
    fCoincidenceWindow = GetDouble("CoincidenceWindow", 0.100);
    // BEGIN space-based
    fCoincidenceRadius = GetDouble("CoincidenceRadius", 1200.); 
    Info("Loaded CoincidenceRadius = %lf [mm]", fCoincidenceRadius);
    if( fCoincidenceRadius < 0. )
      Panic("Negative CoincidenceRadius! No multiplet would be created!");
    fUseRunningRadius  = GetBool("UseRunningRadius", true); 
    Info("Set UseRunningRadius = %s", fUseRunningRadius ? "true" : "false");
    // END space-based
    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 {
        fTimeVariable = TV_RAW;
        Warn("Unrecognized time variable, setting it to Raw");
    }
 
    fDataset = -1;
 
    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";
            fSyncSource = GetString("SyncSource", "jitter.txt");
            fSyncTowerSource = GetString("SyncTowerSource", "jitter_tower.txt");
        }
    }
 
    ev.Add<QDouble>("TotalEnergy");
    ev.Add<QDouble>("Time");
 
    fUseRunningWindow = GetBool("UseRunningWindow", true);
    fStorePosition = GetBool("StorePosition", true);
    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");
    if( fTimeVariable == TV_OF )
        ev.Require<QOFData>(fOFLabel, "OFData");
 
    // BEGIN space-based
    // load position in xyz space of the center of the crystals for space-based coincidences
    fChannelPositionFile = GetString("ChannelPositionFile", "ChannelPosition.dat");
    fUseChannelPosition = GetBool("UseChannelPosition", false);
 
    if( fUseChannelPosition ){
      ev.Add<QDouble>("Radius"); // space-based
      if( FileExists(fChannelPositionFile) ) {
        ifstream fCPFile(fChannelPositionFile.c_str());
        string line;
        stringstream ss;
        int fChannel;
        while(getline(fCPFile, line)) {
            ss.str(line);
            ss >> fChannel;
            ss >> fChannelPositionX[fChannel] >> fChannelPositionY[fChannel] >> fChannelPositionZ[fChannel];
            ss.str("");
            ss.clear();
        }
        fCPFile.close();
      }else{
        Panic("UseChannelPosition = true but ASCII file %s does not exist. Modify ChannelPositionFile in cfg or set UseChannelPosition to false.", fChannelPositionFile.c_str() );
      }
    }else{
      Warn("UseChannelPosition = false. Will skip all space-based selections.");
    }// END space-based
    
}
 
void MFastCoincidence::Do(QEvent& ev) {
 
    const QHeader& hdr = ev.Get<QHeader>("DAQ", "Header");
    const QPulseInfo& pi = ev.Get<QPulseInfo>("DAQ", "PulseInfo");
    const int run = hdr.GetRun();
    const int chan = pi.GetChannelId();
    QDouble energy = ev.GetByLabel<QDouble>(fEnergyLabel);
 
    // Dataset
    if( fDataset < 0 ) {
        GlobalHandle<QInt> dHandle("Dataset");
        GlobalData().Get("", &dHandle, "");
        fDataset = dHandle.Get();
    }
 
    if( fThisRun != run ) {
        fThisRun = run;
        QRunDataHandle rHandle(fThisRun);
        GlobalData().Get("DAQ", &rHandle, "");
        fRunData = rHandle.Get();
    }
 
    // Tower map
    if( ChannelTowerMap.find(chan) == ChannelTowerMap.end() ) {
        QDetChannelCollectionHandle dccHandle;
        dccHandle.SetRun(run);
        GlobalData().Get("",&dccHandle,"");
 
        int tower = dccHandle.Get().Get(chan).fTd.fTdTower;
        int floor = dccHandle.Get().Get(chan).fTd.fTdFloor;
        int position = dccHandle.Get().Get(chan).fTd.fTdPos;
        ChannelTowerMap[chan] = tower;
        ChannelFloorMap[chan] = floor;
        ChannelPositionMap[chan] = position;
    }
 
    double time = GetTime(ev);
    if( GetIteration() == 1 ) {
        // First iteration: fill neighbour map
        int tower = ChannelTowerMap[chan];
        ULong64_t evNumber = pi.GetMasterSample().GetEventNumber();
        ULong64_t tag = evNumber*100 + tower;
 
        timeNeighbours.push_back(new TimeNeighbour(evNumber, tag, chan, time, energy));
    } else {
        // Second iteration: coincidences
        // Initialize QCoincidenceData
        QCoincidenceData& coincData = ev.Get<QCoincidenceData>("CoincidenceData");
        coincData.fOrderInMultiple = 1;
        coincData.fMultiplicity = 1;
        coincData.fCoincidentChannels.clear();
        double LargestDeltaT = 0;
        double SmallestDeltaT = Q_DOUBLE_DEFAULT;
        double tot_energy = energy;
 
        ev.Get<QDouble>("Time") = time;
 
        // Find the current QEvent in the neighbours vector
        int tower = ChannelTowerMap[chan];
        ULong64_t thisTag = pi.GetMasterSample().GetEventNumber()*100 + tower;
        //std::vector<TimeNeighbour*>::iterator thisEv = std::find_if(timeNeighbours.begin(), timeNeighbours.end(), FindTimeNeighbour(thisTag));
 
        int p = 0; 
        int u = timeNeighbours.size()-1;
        int m;
        while(p <= u) {
            m = (p+u)/2;
            if( timeNeighbours[m]->fTag == thisTag )
                break;
            if( timeNeighbours[m]->fTime < time )
                p = m + 1;
            else
                u = m - 1;
        }
 
        TimeNeighbour* tn;
 
        // Main event time
        double eTime = time;
        std::set<int> tnSelected; // store here selected neighbors (their key in the timeNeighbours vector)
        std::set<int>::iterator tnSelIT;
        //      std::map<int> tnSelectedDeltaT; // store here DeltaT with respect to closest event wrt main (backwards compatibility)
 
        // Loop backwards on the timeNeighbours starting from the main (excluded)
        for( int s = m-1; s>= 0; s-- )
          {
            tn = timeNeighbours[s];
            double oTime  = tn->fTime;
            double deltaT = oTime - eTime;
            if( fabs(deltaT) < fCoincidenceWindow ){// time-coincident
              if( fUseRunningWindow )
                eTime = oTime;
              tnSelected.insert( s );
            }else{ // gone out of the coincidence window
              break;
            }
          }
        eTime = time;
        // Loop forward on the timeNeighbours starting from the main (excluded)
        for( size_t s = m+1; s < timeNeighbours.size(); s++ ) {
          tn = timeNeighbours[s];
          double oTime  = tn->fTime;
          double deltaT = oTime - eTime;
          if( fabs(deltaT) < fCoincidenceWindow){// time-coincident
            if( fUseRunningWindow )
              eTime = oTime;
            tnSelected.insert( s );
          }else{ // gone out of the coincidence window
            break;
          }
        }
        // END finding time neighbors.
 
        // BEGIN space-based selection 
        if( fUseChannelPosition ){// enable space-based selection
          std::set<int> tnSpaceBasedAccepted;
          for( tnSelIT = tnSelected.begin(); tnSelIT != tnSelected.end(); ++tnSelIT )
            {
              tn = timeNeighbours[ *tnSelIT ];
              int    oChan     = tn->fChannel;
              double oDistance = GetDistanceMmChannel(chan,oChan);
              if( oDistance <= fCoincidenceRadius )//add to the accepted event list
                tnSpaceBasedAccepted.insert( *tnSelIT );
            }
          int  EmergencyStop = -988; // DEBUG
          if( fUseRunningRadius ){// enabled running radius (tracking)
            bool AddedNewEvent;
            if( EmergencyStop >= 0 )
              Panic("Infinite loop detected! Due to bug in SpaceBased with RunningRadius");
            std::set<int>::iterator tnSBA_IT; // iterator on Space Based Accepted events
            do{
              AddedNewEvent = false;
              for( tnSBA_IT = tnSpaceBasedAccepted.begin(); tnSBA_IT != tnSpaceBasedAccepted.end(); ++tnSBA_IT )
                {// for each (space-based) accepted event ...
                  tn = timeNeighbours[ *tnSBA_IT ]; 
                  int sbChan = tn->fChannel;
                  for( tnSelIT = tnSelected.begin(); tnSelIT != tnSelected.end(); ++tnSelIT )
                    {// ... loop over all the time-selected events (except the main)
                      int oChan        = timeNeighbours[ *tnSelIT ]->fChannel;
                      double oDistance = GetDistanceMmChannel( sbChan, oChan );
                      if( oDistance < fCoincidenceRadius && tnSpaceBasedAccepted.end() == tnSpaceBasedAccepted.find( *tnSelIT ) )
                        {// if you find one that is WITHIN the radius buy NOT (yet) accepted ...
                          tnSpaceBasedAccepted.insert( *tnSelIT ); // ... accept it ...
                          AddedNewEvent = true; // ... and update this flag
                        }
                    }
                  if(AddedNewEvent){// DEBUG
                    //Warn("BREAK: Added new event because of RunningRadius to tag %d",tn->fTag);
                    EmergencyStop++;
                    break;
                  }
                }
            }while( AddedNewEvent );
          } // ------------------------------ END RUNNING RADIUS 
          // overwrite time selection with time && space selection
          tnSelected = tnSpaceBasedAccepted;
          // compute maximum distance in the multiplet
          double MaximumDistance = 0.;
          std::set<int>::iterator tnSelIT_i;
          std::set<int>::iterator tnSelIT_j;
          for( tnSelIT_i = tnSelected.begin(); tnSelIT_i != tnSelected.end(); ++tnSelIT_i)
            {
              int iChan = timeNeighbours[ *tnSelIT_i ]->fChannel;
              double iDistance = GetDistanceMmChannel( chan, iChan ); // distance from MAIN
              if( iDistance > MaximumDistance )
                MaximumDistance = iDistance;
              for( tnSelIT_j = tnSelected.begin(); tnSelIT_j != tnSelected.end(); ++tnSelIT_j)
                {
                  // make it faster: D_ij = D_ji and D_ii = 0 => compute only D_ij where i<j
                  if( *tnSelIT_i >= *tnSelIT_j )
                    continue;
                  int jChan = timeNeighbours[ *tnSelIT_j ]->fChannel;
                  double ijDistance = GetDistanceMmChannel( iChan, jChan );
                  if( ijDistance > MaximumDistance )
                    MaximumDistance = ijDistance;
                }
            }
          // fill "radius" with the maximum distance between coincident channels in this multiplet
          ev.Get<QDouble>("Radius") = MaximumDistance;
          // this is just a cross-check: self-coincident events happen when there is a trigger on baseline fluctuation and a pulse 1s after
          if( tnSelected.size()  > 0 && MaximumDistance < 1e-3){//multiplicity > 1 && Radius = 0
            Debug("Detected event with multiplicity > 1 and Radius < 0.001 mm. Tag %d. This is a normal behaviour.",thisTag);
            /*
            Debug("Main %d [tag=%d]\t NsTime/1e9 = %lf",chan, thisTag, double(hdr.GetTime().GetFromStartRunNs())/1e09);
            Debug("Main %d [tag=%d]\t jitter-corrected-time = %lf",chan, thisTag, time);
            for( tnSelIT_i = tnSelected.begin(); tnSelIT_i != tnSelected.end(); ++tnSelIT_i)
              Debug("Main %d [tag=%d]\t Coincident %d [tag=%d]",chan, thisTag, timeNeighbours[ *tnSelIT_i ]->fChannel, timeNeighbours[ *tnSelIT_i ]->fTag);
            */
          }
 
        }
        // END space-based
        
        // Write data on the QEvent
        coincData.fMultiplicity  = tnSelected.size() + 1; // # of coincident events + 1 main event
        for( tnSelIT = tnSelected.begin(); tnSelIT != tnSelected.end(); ++tnSelIT )
          { // looping over selected events
            tn = timeNeighbours[ *tnSelIT ];
            double oTime   = tn->fTime;
            double oEnergy = tn->fEnergy;
            int    oChan   = tn->fChannel;
            double deltaT  = oTime - time; // this is ALWAYS the time difference wrt the main (can be > fCoincidenceWindow)
            if( *tnSelIT < m )// raise OrderInMultiple just for events that happen BEFORE the main
              coincData.fOrderInMultiple++;
            if( fabs(deltaT) > fabs(LargestDeltaT) )
              LargestDeltaT = deltaT;
            if( fabs(deltaT) < fabs(SmallestDeltaT) )
              SmallestDeltaT = deltaT;
            tot_energy  += oEnergy; 
            QCoincidentChannel cc;
            cc.fChannelId = oChan;
            cc.fDeltaT    = deltaT;
            cc.fEnergy    = oEnergy;
            cc.fDeltaTower    = ChannelTowerMap[oChan] - ChannelTowerMap[chan];
            cc.fDeltaFloor    = ChannelFloorMap[oChan] - ChannelFloorMap[chan];
            cc.fDeltaPosition = ChannelPositionMap[oChan] - ChannelPositionMap[chan];
            cc.fEventNumber   = tn->fEventNumber;
            coincData.fCoincidentChannels.push_back(cc);
          }
        
        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;
        }
        ev.Get<QDouble>("TotalEnergy") = tot_energy;
    }
}
 
void MFastCoincidence::Done() {
  if( GetIteration() == 1 ) {
    SetRunAgain(true);
    std::sort(timeNeighbours.begin(), timeNeighbours.end(), timeSort);
  } else {
    SetRunAgain(false);
    for(size_t s = 0; s < timeNeighbours.size(); s++)
      delete timeNeighbours[s];
    
    timeNeighbours.clear();
    ChannelTowerMap.clear();
    ChannelFloorMap.clear();
    ChannelPositionMap.clear();
    fChannelPositionX.clear();
    fChannelPositionY.clear();
    fChannelPositionZ.clear();
    fSyncTimeDifference.clear();
  }
}
 
 
double MFastCoincidence::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; 
  
  switch(fTimeVariable) {
  case TV_RAW:
    {       
      break;
    }
  case TV_OF:
    {
      const QOFData& of = ev.Get<QOFData>(fOFLabel.c_str(),"OFData");
      double samplingFrequency = fRunData.GetChannelRunData(pi.GetChannelId()).fSamplingFrequency;
      time += of.GetDelay()/1000 - double(pi.GetMasterSample().GetSampleIndex())/samplingFrequency;
      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.;
      
    }
  }
  
  // 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.;
    }
  }
  return time;
}
 
// space-based 
double MFastCoincidence::GetDistanceMmChannel(const int ch1, const int ch2){
  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) );     
}