MCoincidenceMultiplicity.cc

Go to the documentation of this file.

#include "MCoincidenceMultiplicity.hh"
#include "QBaseType.hh"
#include "QError.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QRunDataHandle.hh"
#include "QDetChannelCollectionHandle.hh"
#include "QTime.hh"
#include <set>
#include <vector>
 
REGISTER_MODULE(MCoincidenceMultiplicity)
 
using std::set;
using std::vector;
 
using namespace Diana;
 
 
void MCoincidenceMultiplicity::Init(QEvent& ev)
{
    if ( GetIteration() == 1 ) {
        fCoincidenceWindow = GetDouble("CoincidenceWindow", 0.100);
                  
        fUseRunningWindow = GetBool("UseRunningWindow", true);
 
 
        fEnergyLabel = GetString("EnergyLabel","ApplyCalibration@Energy");
        
        fUseThermistorRanking = false;
 
        std::string tmp = GetString("ThermistorRankingLabel",Q_STRING_DEFAULT);
        if(tmp != Q_STRING_DEFAULT) {
          fThermistorRankingLabel = tmp;
          fUseThermistorRanking = true;
        }
 
        ev.Add<QInt>("Multiplicity");
        ev.Add<QInt>("OrderInMultiple");
        ev.Add<QDouble>("TotalEnergy");
        ev.Add<QInt>("FirstEventNumber");
 
 
    }
    ev.Require<QPulseInfo>("DAQ","PulseInfo");
    ev.Require<QHeader>("DAQ","Header");
}
 
void MCoincidenceMultiplicity::Do(QEvent& ev)
{
 
    const QPulseInfo& pulseInfo = ev.Get<QPulseInfo>("DAQ","PulseInfo");
    const QHeader& header = ev.Get<QHeader>("DAQ","Header");
 
    const int& run     = header.GetRun();
    const int& channel = pulseInfo.GetChannelId();
    const double  time_sec = header.GetTime().GetFromStartRunNs() / 1e9;
 
    const double energy = ev.GetByLabel<QDouble>(fEnergyLabel);
    const int& event_number = header.GetEventNumber();
    
    // fChannels will contain a list of the first channel number on each crystal
    // If no ThermistorRankingVariable is specified in the cfg file, only
    // channels in fChannels will be used to build multiple events.
    if ( fChannels.empty() ) {
        QRunDataHandle rHandle(run);
        GlobalData().Get("",&rHandle,"");
        vector<int> allChannels = rHandle.Get().fBolometerChannels;
        
            // as the fChannels vector is built up with, those channels and
            // their relatives are put into this set to avoid putting any
            // relatives into fChannels
        set<int> channelsAndRelatives;
        QDetChannelCollectionHandle dccHandle;
        dccHandle.SetRun(run);
        GlobalData().Get("",&dccHandle,"");
 
        for ( vector<int>::const_iterator ch = allChannels.begin();
              ch != allChannels.end();
              ++ch ) {
           vector<int> relativeChannels = dccHandle.Get().GetRelativeChannels(channel);
 
            if ( channelsAndRelatives.count(*ch) == 0 ) {
                fChannels.insert(*ch);
                channelsAndRelatives.insert(*ch);
                channelsAndRelatives.insert(
                    relativeChannels.begin(), relativeChannels.end()
                );
            }
        }
    }
    
    const bool useThermistorRankingVariable = fUseThermistorRanking;
       
    if ( GetIteration() == 1) {
        bool useEvent = false;
        
        if (useThermistorRankingVariable) {
            const int thermistorRanking= ev.GetByLabel<QInt>(fThermistorRankingLabel);
 
            if (thermistorRanking <= 1) {
                useEvent = true;
            }
        }
        else if ( fChannels.count(channel) != 0 ) {
            useEvent = true;
        }
        
        if (useEvent) {
            EventInfo event;
            event.fChannel     = channel;
            event.fEnergy      = energy;
            event.fEventNumber = event_number;
            event.fTime        = time_sec;
            
            fEventInfos.push_back(event);
            fEventInfosIndex[event_number] = fEventInfos.size() - 1;
        }
    }
    else if ( GetIteration() == 2 ) {
        int multiplicity = -1;
        int orderInMultiple = -1;
        double totalEnergy = -1;
        int firstEventNumber = -1;
        if ( fEventInfosIndex.count(event_number) != 0 ) {
            int index = fEventInfosIndex[event_number];
            
            multiplicity     = fEventInfos[index].fMultiplicity;
            orderInMultiple  = fEventInfos[index].fOrderInMultiple;
            totalEnergy      = fEventInfos[index].fTotalEnergy;
            firstEventNumber = fEventInfos[index].fFirstEventNumber;
        }
 
        ev.Get<QInt>("Multiplicity") = multiplicity;
        ev.Get<QInt>("OrderInMultiple") = orderInMultiple;
        ev.Get<QInt>("FirstEventNumber") = firstEventNumber;
        ev.Get<QDouble>("TotalEnergy") = totalEnergy;
 
 
 
    }
 
}
 
void MCoincidenceMultiplicity::Done()
{
    if ( GetIteration() == 1 ) {
        int indexFirst = 0; // index of first event in mulitiple
        int indexLast  = 0; // index of last event in multiple
                  double DeltaT = 0;
        for (int index = 0; index < (int) fEventInfos.size(); ++index) {
                          
            if(fUseRunningWindow) {
                DeltaT = fEventInfos[index].fTime
                       - fEventInfos[indexLast].fTime;
            }
            else {
                DeltaT = fEventInfos[index].fTime
                       - fEventInfos[indexFirst].fTime;
            }
            
            if (DeltaT < fCoincidenceWindow) {
                indexLast = index;
            }
            else {
                double totalEnergy = 0;
                for (int i = indexFirst; i <= indexLast; ++i) {
                    totalEnergy += fEventInfos[i].fEnergy;
                }
 
                int order = 1;
                for (int i = indexFirst; i <= indexLast; ++i) {
                    fEventInfos[i].fMultiplicity = indexLast - indexFirst + 1;
                    fEventInfos[i].fOrderInMultiple = order;
                    ++order;
                    fEventInfos[i].fFirstEventNumber
                        = fEventInfos[indexFirst].fEventNumber;
                    fEventInfos[i].fTotalEnergy = totalEnergy;
                }
 
                indexFirst = index;
                indexLast  = index;
            }
        }
 
        // handle last event(s)
        double totalEnergy = 0;
        for (int i = indexFirst; i <= indexLast; ++i) {
            totalEnergy += fEventInfos[i].fEnergy;
        }
 
        int order = 1;
        for (int i = indexFirst; i <= indexLast; ++i) {
            fEventInfos[i].fMultiplicity = indexLast - indexFirst + 1;
            fEventInfos[i].fOrderInMultiple = order;
            ++order;
            fEventInfos[i].fFirstEventNumber
                = fEventInfos[indexFirst].fEventNumber;
            fEventInfos[i].fTotalEnergy = totalEnergy;
        }
        
        SetRunAgain(true);
    }
}