MCoincidenceMultiplicity.cc

Go to the documentation of this file.

#include "MCoincidenceMultiplicity.hh"
#include "QCalibData.hh"
#include "QDbDetector.hh"
#include "QError.hh"
#include "QEvent.hh"
#include "QFiltersData.hh"
#include "QMainPulse.hh"
#include "QRawEvent.hh"
#include "QTime.hh"
#include <set>
#include <vector>
 
REGISTER_MODULE(MCoincidenceMultiplicity)
 
using std::set;
using std::vector;
 
MCoincidenceMultiplicity::MCoincidenceMultiplicity(QSequence* s)
    : QModule("CoincidenceMultiplicity",s)
{
}
 
MCoincidenceMultiplicity::~MCoincidenceMultiplicity()
{
}
 
void MCoincidenceMultiplicity::Init()
{
    if ( GetIteration() == 1 ) {
        fCoincidenceWindow = GetDouble("CoincidenceWindow", 0.100);
                  
        fUseRunningWindow = GetBool("UseRunningWindow", true);
 
        fRelativeChannelsSource = GetString("RelativeChannelsSource",
                                            Q_STRING_DEFAULT);
        
        fThermistorRankingVariable = GetString("ThermistorRankingVariable",
                                               Q_STRING_DEFAULT, false);
        
        fPrefix = GetString("AuxDataPrefix", "", false);
    }
}
 
QEvent* MCoincidenceMultiplicity::Do(QEvent* ev)
{
    const QRawEvent& rawEvent = ev->GetRawEvent();
    
    const int& channel      = rawEvent.GetPulse().GetChannelId();
    const int& run          = rawEvent.GetRun();
    const int& event_number = rawEvent.GetEventNumber();
    const double& time_sec  = rawEvent.GetTime().GetFromStartRunNs() / 1e9;
    const double& energy    = ev->GetCalib().GetEnergy();
    
    // 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() ) {
        vector<int> allChannels( QDbDetector::ActiveChannelList(run) );
        
        set<int> channelsAndRelatives;
            // 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
        
        for ( vector<int>::const_iterator ch = allChannels.begin();
              ch != allChannels.end();
              ++ch ) {
            vector<int> relativeChannels;
            if (fRelativeChannelsSource != "NONE") {
                // by specifying RelativeChannelsSource = NONE in the cfg file,
                // one allows events on relative channels to be considered
                // multiple events
                try {
                    relativeChannels = QDbDetector::GetRelativeChannels(
                        run,
                        *ch,
                        fRelativeChannelsSource
                    );
                }
                catch (const QError& err) {
                    Warn( err.GetDescription().c_str() );
                }
            }
            
            if ( channelsAndRelatives.count(*ch) == 0 ) {
                fChannels.insert(*ch);
                channelsAndRelatives.insert(*ch);
                channelsAndRelatives.insert(
                    relativeChannels.begin(), relativeChannels.end()
                );
            }
        }
    }
    
    const bool useThermistorRankingVariable
        = (fThermistorRankingVariable != Q_STRING_DEFAULT);
    
    if ( GetIteration() == 1 && ev->GetFilters().GetPassed() ) {
        bool useEvent = false;
        
        if (useThermistorRankingVariable) {
            const int thermistorRanking
                = ev->AuxData().GetInt(fThermistorRankingVariable);
            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->AuxData().SetInt(fPrefix + "Multiplicity", multiplicity, "save");
        ev->AuxData().SetInt(fPrefix + "OrderInMultiple", orderInMultiple,
                             "save");
        ev->AuxData().SetDouble(fPrefix + "TotalEnergy", totalEnergy,
                                "save");
        ev->AuxData().SetInt(fPrefix + "FirstEventNumber", firstEventNumber,
                             "save");
    }
    return ev;
}
 
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);
    }
}