CoincidenceTight.cc

Go to the documentation of this file.

//*****************************
//
// AD Coincidence Tag/Data
// chao@bnl.gov  2011-09-26
//
//*****************************

#include "CoincidenceTight.h"

#include "Context/Context.h"
#include "Context/TimeStamp.h"

#include "Event/HeaderObject.h"
#include "Event/RecHeader.h"
#include "Event/UserDataHeader.h"
#include "DataUtilities/DybArchiveList.h"
#include "DataSvc/IDaqDetailsSvc.h"

#include "Context/ServiceMode.h"
#include "GaudiKernel/SystemOfUnits.h"

#include <iomanip>
using namespace std;
using namespace Gaudi;
using namespace DayaBay;

CoincidenceTight::CoincidenceTight(const string& name, ISvcLocator* svcloc)
    : DybBaseAlg(name, svcloc),
    m_nTagged(0), m_nTaggedTriggers(0)
{
    declareProperty("RecLocation", m_recHeaderLoc="/Event/Rec/AdSimple",
            "the TES path of target event header location");
    declareProperty("Location", m_location="/Event/Data/Physics/CoincidenceTight",
            "the TES path of AD Coincidence Data location");
    declareProperty("TagLocation", m_tagLocation="/Event/Tag/Physics/CoincidenceTight",
            "the TES path of AD Coincidence Tag location");         
    declareProperty("TimeWindow", m_timeWindow=400*Gaudi::Units::microsecond,
        "AD coincidence time window");
    declareProperty("MaxGapWindow", m_maxGapWindow=0.2*Gaudi::Units::second,
        "Maximum Gap between two triggers in one AD");
}


CoincidenceTight::~CoincidenceTight()
{
}


StatusCode CoincidenceTight::initialize()
{
    info() << "initialize: coincidence time window: " << m_timeWindow/Gaudi::Units::microsecond << " microseconds" << endreq;
    m_daqDetailsSvc = svc<IDaqDetailsSvc>("DetailsFromRecordSvc", true);
    if (!m_daqDetailsSvc) {
        error() << "Failed to load IDaqDetailsSvc" << endreq;
        return StatusCode::FAILURE;
    }
         
    StatusCode sc = service("EventDataArchiveSvc", m_archiveSvc);
    if (sc == StatusCode::FAILURE) {
        error() << "Failed to retrieve EventDataArchiveSvc" << endreq;
        return StatusCode::FAILURE;
    }
    
    // initialize
    DetectorId::DetectorId_t detectorId[4] = {
        DetectorId::kAD1, DetectorId::kAD2, DetectorId::kAD3, DetectorId::kAD4
    };
    for (int i=0; i<4; i++) {
        DetectorId::DetectorId_t detector = detectorId[i];
        m_inputHeaders[detector] = vector<const IHeader*>();
        m_adScaledHeaders[detector] = vector<RecHeader*>();
        m_integratedLiveTime_sec[detector] = 0;
        m_integratedGapTime_sec[detector] = 0;
        m_extendTimeMultiplet_sec[detector] = 0;
    }
    m_detName[DetectorId::kAD1] = "AD1"; m_detName[DetectorId::kAD2] = "AD2";
    m_detName[DetectorId::kAD3] = "AD3"; m_detName[DetectorId::kAD4] = "AD4";
    
    
    m_integratedNGaps = 0;
    
    m_execNum = 0;
    return StatusCode::SUCCESS;
}



StatusCode CoincidenceTight::execute()
{   
    // m_execNum++; 
    // info() << "exec#: " << m_execNum << endreq;

    RecHeader* recHeader = get<RecHeader>(m_recHeaderLoc);
    if ( !recHeader ) {
        warning() << "Cannot find header at " << m_recHeaderLoc << endreq;
        return StatusCode::FAILURE;
    } 

    RecHeader* adScaledHeader = get<RecHeader>("/Event/Rec/AdScaled");
    
    TimeStamp ts = recHeader->context().GetTimeStamp();
    DetectorId::DetectorId_t detectorId = recHeader->context().GetDetId();

    TimeStamp lastTS = ts;
    
    if(m_lastTimeStamp.find(detectorId) != m_lastTimeStamp.end()){
        lastTS = m_lastTimeStamp[detectorId];
    }       
    m_lastTimeStamp[detectorId] = ts;
        
    // only look for AD Coincidence
    if (   detectorId == DetectorId::kIWS 
        || detectorId == DetectorId::kOWS 
        || detectorId == DetectorId::kRPC) {
        return StatusCode::SUCCESS; 
    }

    TimeStamp dt = TimeStamp(ts);
    dt.Subtract(lastTS);
    double dt_sec = dt.GetSeconds();
    // info() << "AD: " << detectorId << ", DtSec: " << dt_sec << endreq;

    if (dt_sec*Gaudi::Units::second > m_maxGapWindow) {
        // long gap in DAQ, considered as deadtime in coincidence tree.
        m_integratedGapTime_sec[detectorId] += dt_sec;
        m_integratedNGaps++;
    }
    else {
        m_integratedLiveTime_sec[detectorId] += dt_sec;
    }

    if (dt_sec*Gaudi::Units::second > m_timeWindow) {
        int multiplicity = m_inputHeaders[detectorId].size();
        if (multiplicity < 2) {
            m_inputHeaders[detectorId].clear();
            m_adScaledHeaders[detectorId].clear();
        }
        else {
            // found a gap
            
            if (dt_sec*Gaudi::Units::second > m_maxGapWindow) {
                // protection: long gap in DAQ, the previous multiplets may not in memory any more
                warning() << "Event gap in AD " << detectorId << ": " << dt_sec << " sec" << endreq;
                m_integratedLiveTime_sec[detectorId] -= m_extendTimeMultiplet_sec[detectorId];
            } 
            else {
                // info() << "start saving: dtLast: " << dt_sec << endreq;
                // found a conincident multiplet
                vector<const IHeader*>& inputHeaders = m_inputHeaders[detectorId];
                
                // info() << "making user tag: size of input headers: " << inputHeaders.size() << endreq;
                // save a (redundant) tag header object into the tag location        
                HeaderObject *tag = MakeHeader<HeaderObject>(inputHeaders);                
                put(tag, m_tagLocation);
                                
                // make the UserData                
                UserDataHeader* data = MakeHeader<UserDataHeader>(inputHeaders);
                
                // store user-defined data
                SaveRecData(inputHeaders, data);
                SaveAdScaledData(m_adScaledHeaders[detectorId], data);
                SaveCalibStatsData(inputHeaders, data);
                SaveMuonData(data);
                SaveOtherData(inputHeaders, data);
                
                data->set("integratedLiveTime_sec", m_integratedLiveTime_sec[detectorId]);
                data->set("integratedGapTime_sec", m_integratedGapTime_sec[detectorId]);
                data->set("integratedNGaps", m_integratedNGaps);
                data->set("extendTimeMultiplet_sec", m_extendTimeMultiplet_sec[detectorId]);
                
                if (PreSelection(data)) {
                    // save UserData to TES
                    put(data, m_location);
                    m_nTagged++;
                    m_nTaggedTriggers += multiplicity;
                    // info() << "found one: multiplicty " << multiplicity 
                    //        << " execNumber: " << execNumber << endreq;
                }
                else {
                    delete data;
                }
            }
            
            m_inputHeaders[detectorId].clear();
            m_adScaledHeaders[detectorId].clear();
            m_extendTimeMultiplet_sec[detectorId] = 0;
        }
    }
    else {
        m_extendTimeMultiplet_sec[detectorId] += dt_sec;
        // info() << "extented time AD " << detectorId << " : " << m_extendTimeMultiplet_sec[detectorId] << endreq;
    }
    
    m_inputHeaders[detectorId].push_back(recHeader);
    m_adScaledHeaders[detectorId].push_back(adScaledHeader);
     
    return StatusCode::SUCCESS;
}


StatusCode CoincidenceTight::finalize()
{
    info() << m_nTagged << " multiplets (" << m_nTaggedTriggers 
           << " triggers) tagged at " << m_location  << endreq;
    return StatusCode::SUCCESS;
}


// helper methods
void CoincidenceTight::SaveRecData(vector<const IHeader*>& inputHeaders, UserDataHeader* dataHeader)
{
    // initialize data vectors
    vector<float> e, x, y, z, rawEvis;
    vector<int> triggerNumber, triggerType, t_s, t_ns, dt_ns, energyStatus, positionStatus;
    
    int t0_s = 0; 
    int t0_ns = 0; 
    for (size_t i=0; i<inputHeaders.size(); i++) {
        const RecHeader* recHeader = dynamic_cast<const RecHeader*>(inputHeaders[i]);
        RecTrigger result = recHeader->recTrigger();
        e.push_back(result.energy()/Units::MeV);
        x.push_back(result.position().x()/Units::mm);
        y.push_back(result.position().y()/Units::mm);
        z.push_back(result.position().z()/Units::mm);
        rawEvis.push_back(result.rawEvis()/Units::MeV);
        
        triggerNumber.push_back(result.triggerNumber());
        triggerType.push_back(result.triggerType());

        energyStatus.push_back(result.energyStatus());
        positionStatus.push_back(result.positionStatus());
        int thisT_s = recHeader->timeStamp().GetSec();
        int thisT_ns = recHeader->timeStamp().GetNanoSec();
        if (0 == i) {
            t0_s = thisT_s;
            t0_ns = thisT_ns;
        }
        t_s.push_back(thisT_s);
        t_ns.push_back(thisT_ns);
        dt_ns.push_back( int((thisT_s-t0_s)*1e9 + thisT_ns - t0_ns));
    }
    
    dataHeader->set("e", e);
    dataHeader->set("x", x);
    dataHeader->set("y", y);
    dataHeader->set("z", z);
    dataHeader->set("rawEvis", rawEvis);
    
    dataHeader->set("triggerNumber", triggerNumber);
    dataHeader->set("triggerType", triggerType);
    dataHeader->set("t_s", t_s);
    dataHeader->set("t_ns", t_ns);
    dataHeader->set("dt_ns", dt_ns);
    dataHeader->set("energyStatus", energyStatus);
    dataHeader->set("positionStatus", positionStatus);
}

void CoincidenceTight::SaveAdScaledData(vector<DayaBay::RecHeader*>& inputHeaders, UserDataHeader* dataHeader)
{
    // initialize data vectors
    vector<float> e, x, y, z, rawEvis;
    vector<int> energyStatus, positionStatus;
    
 
    for (size_t i=0; i<inputHeaders.size(); i++) {
        RecHeader* recHeader = inputHeaders[i];
        RecTrigger result = recHeader->recTrigger();
        e.push_back(result.energy()/Units::MeV);
        x.push_back(result.position().x()/Units::mm);
        y.push_back(result.position().y()/Units::mm);
        z.push_back(result.position().z()/Units::mm);
        rawEvis.push_back(result.rawEvis()/Units::MeV);
        energyStatus.push_back(result.energyStatus());
        positionStatus.push_back(result.positionStatus());
    }
    
    dataHeader->set("e_AdScaled", e);
    dataHeader->set("x_AdScaled", x);
    dataHeader->set("y_AdScaled", y);
    dataHeader->set("z_AdScaled", z);
    dataHeader->set("rawEvis_AdScaled", rawEvis);
    dataHeader->set("energyStatus_AdScaled", energyStatus);
    dataHeader->set("positionStatus_AdScaled", positionStatus);
}


void CoincidenceTight::SaveCalibStatsData(vector<const IHeader*>& inputHeaders, UserDataHeader* dataHeader)
{
    int nDataInt = 1;
    int nDataFloat = 23; // + dtLastAD/dtNextAD
    string dataStrInt[] = { "nHit" };
    string dataStrFloat[] = { 
        "nPESum", "nPERMS", "nPulseSum", "nPulseRMS", // PE distribution
        "tEarliest", "tLatest", "tMean", // time distribution
        "EarlyCharge", "NominalCharge", "LateCharge", "PeakCharge", // 4 diffrent charges
        "MaxQ", "Quadrant", "flasher_flag", "time_PSD", "time_PSD1", "MaxQ_2inchPMT", // for flashers
        "dtLastIWS_ms", "dtLastOWS_ms", "dtLastRPC_ms", "dtNextIWS_ms", "dtNextOWS_ms", "dtNextRPC_ms", // last/next trigger (deprecated)
    };
    vector< vector<int> > dataInt(nDataInt);
    vector< vector<float> > dataFloat(nDataFloat);
    vector<float> dtLastAD_ms;
    vector<float> dtNextAD_ms;
    
    DetectorId::DetectorId_t detectorId = dataHeader->context().GetDetId();
    string detName = m_detName[detectorId];
    string name_dtLastAD = string("dtLast")+detName+"_ms";
    string name_dtNextAD = string("dtNext")+detName+"_ms";

    for (size_t i=0; i<inputHeaders.size(); i++) {
        const RecHeader* recHeader = dynamic_cast<const RecHeader*>(inputHeaders[i]);
        vector<const IHeader*> inputHeaders = recHeader->findHeaders("/Event/Data/CalibStats");
        if ( 1 == inputHeaders.size() ) {
            // Fix:: UserDataHeader class need to make getBlahBlah() methods const. 
            UserDataHeader* calibStats =  const_cast<UserDataHeader*>(dynamic_cast<const UserDataHeader*>(inputHeaders[0]));
            for (int j=0; j<nDataInt; j++ ) {
                int value = calibStats->getInt(dataStrInt[j]);
                dataInt[j].push_back(value);
            }
            for (int j=0; j<nDataFloat; j++ ) {                
                float value = calibStats->getFloat(dataStrFloat[j]);
                dataFloat[j].push_back(value);                
            }
            dtLastAD_ms.push_back(calibStats->getFloat(name_dtLastAD));
            dtNextAD_ms.push_back(calibStats->getFloat(name_dtNextAD));
        }
        else {
            warning() << inputHeaders.size() << " Data at /Event/Data/CalibStats" << endreq;
        }
    }
    
    for (int j=0; j<nDataInt; j++ ) {
        dataHeader->set( string("calib_")+dataStrInt[j], dataInt[j] );
    }
    for (int j=0; j<nDataFloat; j++ ) {
        dataHeader->set( string("calib_")+dataStrFloat[j], dataFloat[j] );
    }            
    dataHeader->set("calib_dtLastAD_ms", dtLastAD_ms);
    dataHeader->set("calib_dtNextAD_ms", dtNextAD_ms);
}


void CoincidenceTight::SaveMuonData(UserDataHeader* dataHeader)
{
    vector<int> t_s = dataHeader->getIntArray("t_s");
    vector<int> t_ns = dataHeader->getIntArray("t_ns");
         
    int firstCoinc_t_s = t_s.at(0);
    int firstCoinc_t_ns = t_ns.at(0);
    int multiplicity = t_s.size();

    SmartDataPtr<DybArchiveList> muons(m_archiveSvc, "/Event/Data/Physics/MuonData");
    if (!muons) {
        debug() << "no Muons in archive?? " << endreq;
        dataHeader->set("isGap", 1);
        dataHeader->set("nMuon", 0);
        dataHeader->set("mu_t_s", vector<int>());
        dataHeader->set("mu_t_ns", vector<int>());
        dataHeader->set("mu_nHitIWS", vector<int>());
        dataHeader->set("mu_nHitOWS", vector<int>());
        dataHeader->set("mu_hitRPC", vector<int>());
        dataHeader->set("mu_nPESumAD", vector<float>());
        dataHeader->set("mu_triggerNumberIWS", vector<int>());
        dataHeader->set("mu_triggerNumberOWS", vector<int>());
        dataHeader->set("mu_triggerNumberRPC", vector<int>());
        dataHeader->set("mu_triggerNumberAD", vector<int>());
        dataHeader->set("dtLastWPMuon_ms", vector<float>(multiplicity, 1e7));
        dataHeader->set("dtLastADMuon_ms", vector<float>(multiplicity, 1e7));
        dataHeader->set("dtLastRPCMuon_ms", vector<float>(multiplicity, 1e7));
        dataHeader->set("dtLastShowerMuon_ms", vector<float>(multiplicity, 1e7));
        dataHeader->set("nPESumLastShowerMuon", vector<float>(multiplicity, 0));
        return;
    }
    
    
    DetectorId::DetectorId_t detectorId = dataHeader->context().GetDetId();
    
    string detName = m_detName[detectorId];
    string name_nPESumAD = string("calib_nPESum_")+detName;
    string name_triggerNumberAD = string("triggerNumber_")+detName;
    
    vector<int> mu_t_s, mu_t_ns, mu_nHit_IWS, mu_nHit_OWS, mu_hitRPC;
    vector<int> mu_triggerNumberIWS, mu_triggerNumberOWS, mu_triggerNumberRPC,  mu_triggerNumberAD;
    vector<float> mu_nPESum_AD;
     
    int nMuon=0;
    DybArchiveList::const_iterator iter = muons->begin();
    for (; iter != muons->end(); iter++) {
        UserDataHeader* muonData = dynamic_cast<UserDataHeader*> (*iter);
        int tMu_s = muonData->getInt("tMu_s");
        
        int tMu_ns = muonData->getInt("tMu_ns");
        int calib_nHit_IWS = muonData->getInt("calib_nHit_IWS");
        int calib_nHit_OWS = muonData->getInt("calib_nHit_OWS");
        int hitRPC = muonData->getInt("hitRPC");
        float calib_nPESum_AD = muonData->getFloat(name_nPESumAD);
        int triggerNumberIWS = muonData->getInt("triggerNumber_IWS");
        int triggerNumberOWS = muonData->getInt("triggerNumber_OWS");
        int triggerNumberRPC = muonData->getInt("triggerNumber_RPC");
        int triggerNumberAD  = muonData->getInt(name_triggerNumberAD);

        double dtMax_ms = (double(firstCoinc_t_s-tMu_s)*1e9 + (firstCoinc_t_ns-tMu_ns))/1e6;
        if (dtMax_ms > 2) {
            // store shower muons between (2ms, 1s)
            if (dtMax_ms>1000) break;
            if (calib_nPESum_AD > 2e5) {
                mu_t_s.push_back(tMu_s);
                mu_t_ns.push_back(tMu_ns);
                mu_nHit_IWS.push_back(calib_nHit_IWS);
                mu_nHit_OWS.push_back(calib_nHit_OWS);
                mu_hitRPC.push_back(hitRPC);
                mu_nPESum_AD.push_back(calib_nPESum_AD);
                mu_triggerNumberIWS.push_back(triggerNumberIWS);
                mu_triggerNumberOWS.push_back(triggerNumberOWS);
                mu_triggerNumberRPC.push_back(triggerNumberRPC);
                mu_triggerNumberAD.push_back(triggerNumberAD);
                nMuon++;                
            }
            continue;
        }
        else { // store non-shower muons between (~-0.4ms, 2ms)
            // cout << tMu_s << ", " << tMu_ns << endl;
            // cout << firstCoinc_t_s << ", " << firstCoinc_t_ns << endl;
            // cout << dtMax_ms << endl;
            // cout << "---" << endl;
            mu_t_s.push_back(tMu_s);
            mu_t_ns.push_back(tMu_ns);
            mu_nHit_IWS.push_back(calib_nHit_IWS);
            mu_nHit_OWS.push_back(calib_nHit_OWS);
            mu_hitRPC.push_back(hitRPC);
            mu_nPESum_AD.push_back(calib_nPESum_AD);
            mu_triggerNumberIWS.push_back(triggerNumberIWS);
            mu_triggerNumberOWS.push_back(triggerNumberOWS);
            mu_triggerNumberRPC.push_back(triggerNumberRPC);
            mu_triggerNumberAD.push_back(triggerNumberAD);
            nMuon++;
        }

    }
    // info() << "# Muon: " << nMuon << endreq; // 1ms before first to 200us after last
    
    // for convenience to use in TTree::Draw(), however if need redefintion of muons, please don't use those
    vector<float> dtLastWPMuon_ms(multiplicity); 
    vector<float> dtLastADMuon_ms(multiplicity); 
    vector<float> dtLastRPCMuon_ms(multiplicity); 
    vector<float> dtLastShowerMuon_ms(multiplicity); 
    vector<float> nPESumLastShowerMuon(multiplicity); 
    
    // cout << "multiplicity: " << multiplicity << endl;
    for (int i=0; i<multiplicity; i++) {
        dtLastWPMuon_ms[i] = 1e7; // any large number
        dtLastADMuon_ms[i] = 1e7; 
        dtLastRPCMuon_ms[i] = 1e7;
        dtLastShowerMuon_ms[i] = 1e7;
        nPESumLastShowerMuon[i] = 0;
        
        bool hasSetWPMuon = false;
        bool hasSetADMuon = false;
        bool hasSetRPCMuon = false;
        bool hasSetShowerMuon = false;
        for (int j=0; j<nMuon; j++) {
            double dtMuon_ms = (double(t_s[i]-mu_t_s[j])*1e9 + (t_ns[i]-mu_t_ns[j]))/1e6;
            if (dtMuon_ms < -2e-3) { continue; } // too far in the future, must not be correlated
            
            bool isWPMuon = mu_nHit_IWS[j]>12 || mu_nHit_OWS[j]>12;
            bool isADMuon = mu_nPESum_AD[j] > 3000;
            bool isRPCMuon = mu_hitRPC[j] == 1; 
            bool isShowerMuon = mu_nPESum_AD[j] > 5.5e5;
            
            if (isWPMuon && !hasSetWPMuon) { dtLastWPMuon_ms[i] = dtMuon_ms; hasSetWPMuon = true; }
            if (isADMuon && !hasSetADMuon) { dtLastADMuon_ms[i] = dtMuon_ms; hasSetADMuon = true; }
            if (isRPCMuon && !hasSetRPCMuon) { dtLastRPCMuon_ms[i] = dtMuon_ms; hasSetRPCMuon = true; }
            if (isShowerMuon && !hasSetShowerMuon) { 
                dtLastShowerMuon_ms[i] = dtMuon_ms; 
                nPESumLastShowerMuon[i] = mu_nPESum_AD[j]; 
                hasSetShowerMuon = true; 
            }
            
        }
        
    }
    dataHeader->set("isGap", 0);
    dataHeader->set("nMuon", nMuon);
    
    dataHeader->set("mu_t_s", mu_t_s);
    dataHeader->set("mu_t_ns", mu_t_ns);
    dataHeader->set("mu_nHitIWS", mu_nHit_IWS);
    dataHeader->set("mu_nHitOWS", mu_nHit_OWS);
    dataHeader->set("mu_hitRPC", mu_hitRPC);
    dataHeader->set("mu_nPESumAD", mu_nPESum_AD);
    dataHeader->set("mu_triggerNumberIWS", mu_triggerNumberIWS);
    dataHeader->set("mu_triggerNumberOWS", mu_triggerNumberOWS);
    dataHeader->set("mu_triggerNumberRPC", mu_triggerNumberRPC);
    dataHeader->set("mu_triggerNumberAD", mu_triggerNumberAD);
    
    dataHeader->set("dtLastWPMuon_ms", dtLastWPMuon_ms);
    dataHeader->set("dtLastADMuon_ms", dtLastADMuon_ms);
    dataHeader->set("dtLastRPCMuon_ms", dtLastRPCMuon_ms);
    dataHeader->set("dtLastShowerMuon_ms", dtLastShowerMuon_ms);
    dataHeader->set("nPESumLastShowerMuon", nPESumLastShowerMuon);
    
}


void CoincidenceTight::SaveOtherData(vector<const IHeader*>& inputHeaders, UserDataHeader* dataHeader)
{
    // initialize data vectors
    vector<int> I, J;
    
    int multiplicity = inputHeaders.size();
    for (int i=0; i<multiplicity-1; i++) {
        for (int j=i+1; j<multiplicity; j++) {
            I.push_back(i);
            J.push_back(j);
        }
    }
    
    dataHeader->set("multiplicity", multiplicity);
    dataHeader->set("I", I);
    dataHeader->set("J", J);
    
    // get from CoincidenceLoose if exist (a hack)
    // SmartDataPtr<DybArchiveList> coinc(m_archiveSvc, "/Event/Data/Physics/CoincidenceLoose");
    // if (coinc) {
    //    DybArchiveList::const_iterator iter = coinc->begin();
    //    UserDataHeader* first = dynamic_cast<UserDataHeader*> (*iter);
    //    dataHeader->set("runNo", first->getInt("runNo"));
    //    dataHeader->set("fileNo", first->getInt("fileNo"));
    //    return;
    // }
    
    // get from file header (raw file)
    int runNo = 0;
    int fileNo = 0;
    ServiceMode svcMode = ServiceMode(dynamic_cast<const HeaderObject*>(inputHeaders.at(0))->context(), 0);
    const DaqRunDetails& runDetails = m_daqDetailsSvc->runDetails(svcMode);
    if (&runDetails) {
        runNo = runDetails.runNumber();
    }
    const DaqFileDetails& fileDetails = m_daqDetailsSvc->fileDetails(svcMode);
    if (&fileDetails) {
        fileNo = fileDetails.fileNumber();
    }
    dataHeader->set("runNo", runNo);
    dataHeader->set("fileNo", fileNo);
}

// returns false for obvious useless multiplets
// (some obvious non-IBD multipets are still useful for other studies)
bool CoincidenceTight::PreSelection(UserDataHeader* dataHeader)
{
    int multiplicity = dataHeader->getInt("multiplicity");
    
    // Pre Energy cut, only true if all events in multipets satisfy
    float cutElow = 0.5; float cutEHigh = 80; bool allOutsideCutE = true;
    vector<float> rawEvis = dataHeader->getFloatArray("rawEvis");
    for (int i=0; i<multiplicity; i++) {
        if (rawEvis.at(i) > cutElow && rawEvis.at(i) < cutEHigh) {
            allOutsideCutE = false;
            break;
        }
    }
    if (allOutsideCutE) return false;
    
    return true;
}