MAlphaDelayedCoincidence.cc

Go to the documentation of this file.

#include "MAlphaDelayedCoincidence.hh"
#include "QEvent.hh"
#include "QRawEvent.hh"
#include "QCoincidenceData.hh"
#include "QBaseType.hh"
#include "QOFData.hh"
#include "QMatrix.hh"
#include "QDetChannelCollectionHandle.hh"
#include <cmath>
 
 
using namespace Diana;
using namespace std;
 
REGISTER_MODULE(MAlphaDelayedCoincidence)
 
 
void MAlphaDelayedCoincidence::Init(QEvent& ev)
{
        // called before event loop
        fEnergyLabel = GetString("EnergyLabel","ApplyCalibration@Energy");  
 
        vector<string> defaults;
        fPSALabels = GetVectorString("PSALabels",defaults);
        fCoincidenceDataOwner = GetString("CoincidenceDataOwner","Coincidence",false);
        fMotherSelectionLabel = GetString("MotherSelectionLabel",Q_STRING_DEFAULT);
        fDaughterSelectionLabel = GetString("DaughterSelectionLabel",Q_STRING_DEFAULT);
        // events flagged as mother nuclides
        ev.Add<QCoincidenceData>("CoincidenceData");
        ev.RequireByLabel<QBool>(fMotherSelectionLabel);
        ev.RequireByLabel<QBool>(fDaughterSelectionLabel);
        if(!fPSALabels.empty())  {
                for(size_t p = 0; p < fPSALabels.size(); p++) {
                        ev.RequireByLabel<QDouble>(fPSALabels[p]);
                }
        }
 
        GetBool("SkipEvents",false,false);
 
}
 
void MAlphaDelayedCoincidence::Do(QEvent& ev)
{
        const bool isMother = ev.GetByLabel<QBool>(fMotherSelectionLabel);
        const bool isDaughter = ev.GetByLabel<QBool>(fDaughterSelectionLabel);
        if(!isMother && !isDaughter) return;
 
        // get info on this event
        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();
        QDetChannelCollectionHandle dccHandle;
        dccHandle.SetRun(run);
        GlobalData().Get("",&dccHandle,"");
        int floor = dccHandle.Get().Get(chan).fTd.fTdFloor;
        const QCoincidenceData& coincData = ev.Get<QCoincidenceData>(fCoincidenceDataOwner.c_str(),"CoincidenceData");
        const double time = ev.Get<QDouble>(fCoincidenceDataOwner.c_str(),"Time");
 
        // Act on this event as daughter
 
        if(isDaughter){
                // reset data related to mothers
                QCoincidenceData& backData = ev.Get<QCoincidenceData>("CoincidenceData");
                backData.fOrderInMultiple = -1; // meaningless
                backData.fMultiplicity = 1;
                backData.fCoincidentChannels.clear();
                backData.fNearestCoincidentIndex = 0;
                backData.fFarthestCoincidentIndex = 0;
 
                // look for mothers of this multiplet closest, with time precedence
                std::map<double,int>::reverse_iterator timeOrderedFloorIter = fTimeOrderedFloors.rbegin();
                while(timeOrderedFloorIter != fTimeOrderedFloors.rend()) {
                        const int o_floor = timeOrderedFloorIter->second;
                        const double deltaT = timeOrderedFloorIter->first-time;
                        // veto 2 seconds
                        if(fabs(deltaT) < 2.) {
                                timeOrderedFloorIter++;
                                continue;
                        }
                        int deltaFloor = o_floor-floor;
 
                        // in case this event is a multiplet, check whether deltaFloor 
                        // of an event in the multiplet is smaller than 3
                        const std::vector<QCoincidentChannel>& multChannels = coincData.fCoincidentChannels;
                        int numMultChannels = multChannels.size();
                        int countMultChannels = numMultChannels;
                        while(abs(deltaFloor) > 2 && countMultChannels>0) {
                                const QCoincidentChannel& multChannel = multChannels[numMultChannels-countMultChannels];
                                deltaFloor = o_floor-(floor+multChannel.fDeltaFloor);
                                countMultChannels--;
                        }
 
                        if(abs(deltaFloor) <=2) {
                                // store in this event the info of the back alpha(s)
                                QCoincidentChannel backCoinc = fFloors[o_floor];
                                backCoinc.fDeltaT -= time;
                                backCoinc.fDeltaFloor = deltaFloor;
                                backData.fMultiplicity++;
                                backData.fCoincidentChannels.push_back(backCoinc);
                                backData.fFarthestCoincidentIndex =  backData.fMultiplicity-2;
                        }
                        timeOrderedFloorIter++;
                }
        }
 
        // Act on this event as mother
 
        if(isMother && coincData.fMultiplicity <= 2 ) {
                const double total_energy = ev.Get<QDouble>(fCoincidenceDataOwner.c_str(),"TotalEnergy");
                QCoincidentChannel backChannel;
                backChannel.fChannelId = chan;
                backChannel.fDeltaT = time;
                backChannel.fEnergy = total_energy;
                for(size_t p = 0; p < fPSALabels.size(); p++) {
                        backChannel.fPSAVariables.push_back(ev.GetByLabel<QDouble>(fPSALabels[p]));
                }
                backChannel.fDeltaFloor = -1; 
                backChannel.fDeltaPosition = -1;
                if(coincData.fMultiplicity == 1) {
                        fFloors[floor] = backChannel;
                } else if(coincData.fMultiplicity == 2) {
                        const double energy = ev.GetByLabel<QDouble>(fEnergyLabel); 
                        const double energy2 = coincData.fCoincidentChannels[0].fEnergy;
                        // the crystal with less energy contains the contamination so we save only this among the two
                        if(energy < energy2) {
                                fFloors[floor]  = backChannel;
                        }
                }
                // update time ordering of floors
                fTimeOrderedFloors.clear();
                std::map<int,QCoincidentChannel>::const_iterator floorIter = fFloors.begin();
                while(floorIter != fFloors.end()) {
                        const int backFloor = floorIter->first;
                        const double backTime = floorIter->second.fDeltaT;
                        fTimeOrderedFloors[backTime] = backFloor;
                        floorIter++;
                }
        }
        return;
}
 
void MAlphaDelayedCoincidence::Done()
{
        // called at the end of the event loop
}