: EDProducer{p}

, fIsSimple(p.get<bool>("IsSimple", true))

, fFlashType(p.get<std::string>("FlashType", "simpleflash_pmt"))

, fUseOldMetrics(p.get<bool>("UseOldMetrics", false))

, fPandoraProducer(p.get<std::string>("PandoraProducer"))

, fSpacePointProducer(p.get<std::string>("SpacePointProducer"))

, fOpHitProducer(p.get<std::string>("OpHitProducer"))

, fOpHitARAProducer(p.get<std::string>("OpHitARAProducer", ""))

, fOpFlashProducer(p.get<std::vector<std::string>>("OpFlashProducer", {}))

, fOpFlashHitProducer(p.get<std::vector<std::string>>("OpFlashHitProducer", {}))

, fVerbose(p.get<bool>("Verbose", false))

// , fCaloProducer(p.get<std::string>("CaloProducer"))

// , fTrackProducer(p.get<std::string>("TrackProducer"))

, fBeamSpillTimeStart(p.get<double>("BeamSpillTimeStart")) //us

, fBeamSpillTimeEnd(p.get<double>("BeamSpillTimeEnd"))// us

, fFlashFindingTimeStart(p.get<double>("FlashFindingTimeStart")) //us

, fFlashFindingTimeEnd(p.get<double>("FlashFindingTimeEnd"))// us

, fFlashStart(p.get<double>("FlashStart")) // in us w.r.t. flash time

, fFlashEnd(p.get<double>("FlashEnd")) // in us w.r.t flash time

, fTimeBins(timeBins())

, fSelectNeutrino(p.get<bool>("SelectNeutrino", true)) // only attempt to match potential neutrino slices

, fForceConcurrence(p.get<bool>("ForceConcurrence", false)) // require light and charge to coincide, different requirements for SBND and ICARUS

, fUse3DMetrics(p.get<bool>("Use3DMetrics", false)) // use metrics that depend on (X,Y,Z)

, fUseOpCoords(p.get<bool>("UseOpCoords", true)) // Use precalculated OpFlash coordinates

, fCorrectDriftDistance(p.get<bool>("CorrectDriftDistance", false)) // require light and charge to coincide, different requirements for SBND and ICARUS

, fStoreMCInfo(p.get<bool>("StoreMCInfo", false))

// , fUseCalo(p.get<bool>("UseCalo", false))

, fRM(loadMetrics(p.get<std::string>("InputFileName")))

, fNoAvailableMetrics(p.get<bool>("NoAvailableMetrics", false))

, fMakeTree(p.get<bool>("MakeTree", false))

, fChargeToNPhotonsShower(p.get<double>("ChargeToNPhotonsShower", 1.0)) // ~40000/1600

, fChargeToNPhotonsTrack(p.get<double>("ChargeToNPhotonsTrack", 1.0)) // ~40000/1600

, fMinHitQ(p.get<double>("MinHitQ", 0.0))

, fMinSpacePointQ(p.get<double>("MinSpacePointQ", 0.0))

, fMinParticleQ(p.get<double>("MinParticleQ", 0.0))

, fMinSliceQ(p.get<double>("MinSliceQ", 0.0))

, fOpHitTime(p.get<std::string>("OpHitTime", "RiseTime"))

, fUseOpHitRiseTime("RiseTime" == fOpHitTime)

, fUseOpHitPeakTime("PeakTime" == fOpHitTime)

, fUseOpHitStartTime("StartTime" == fOpHitTime)

, fMinInTimeFlashes(p.get<unsigned>("MinInTimeFlashes", 1))

, fMaxFlashes(p.get<unsigned>("MaxFlashes", fMinInTimeFlashes))

, fMinOpHPE(p.get<double>("MinOpHPE", 0.0))

, fMinFlashPE(p.get<double>("MinFlashPE", 0.0))

, fFlashPEFraction(p.get<double>("FlashPEFraction", 0.8))

, fWireReadoutGeom{&art::ServiceHandle<geo::WireReadout>()->Get()}

, fDetector(detectorName(fGeometry->DetectorName()))

, fSBND((fDetector == "SBND") ? true : false )

, fICARUS((fDetector == "ICARUS") ? true : false )

, fPlaneList(p.get<std::vector<int>>("PlaneList", (fSBND) ? kSBNDPlanes : kICARUSPlanes)) // Use 0, 1, 2 for SBND and 0, 1, 3 for ICARUS

, fAllPlanes((fSBND && fPlaneList==kSBNDPlanes) || (fICARUS && fPlaneList==kICARUSPlanes))

, fPDMapAlgPtr(art::make_tool<opdet::PDMapAlg>(p.get<fhicl::ParameterSet>("PDMapAlg")))

, fNTPC(fGeometry->NTPC())

, fTPCPerDriftVolume(fNTPC/2) // 2 drift volumes: kRght and kLeft

, fCryostat(p.get<int>("Cryostat", 0)) //set =0 or =1 for ICARUS to match reco chain selection

, fGeoCryo(std::make_unique<geo::CryostatGeo>(fGeometry->Cryostat(geo::CryostatID(fCryostat))))

, fWiresX_gl(wiresXGl())

, fDriftDistance(driftDistance())

, fXBins(p.get<double>("XBins"))

, fXBinWidth(fDriftDistance/fXBins)// cm

// , fRR_TF1_fit(p.get<std::string>("rr_TF1_fit", "pol3")) // LEGACY

// , fRatio_TF1_fit(p.get<std::string>("ratio_TF1_fit", "pol3")) // LEGACY

, fYBins(p.get<unsigned>("YBins", 0.)) // roughly match the rows of opdets

, fZBins(p.get<unsigned>("ZBins", 0.)) // roughly match the columns of opdets

, fYLow(p.get<double>("YLow", 0.)) // lowest opdet position in cm

, fYHigh(p.get<double>("YHigh", 0.)) // highest opdet position in cm

, fZLow(p.get<double>("ZLow", 0.)) // most upstream opdet position in cm

, fZHigh(p.get<double>("ZHigh", 0.)) // most downstream opdet position in cm

, fSkewLimitY(p.get<double>("SkewLimitY", 10.)) //

, fSkewLimitZ(p.get<double>("SkewLimitZ", 10.)) //

, fOpDetNormalizer((fSBND) ? 4 : 1)

, fTermThreshold(p.get<double>("ThresholdTerm", 30.))

{

produces< std::vector<sbn::SimpleFlashMatch> >();

produces< art::Assns <recob::PFParticle, sbn::SimpleFlashMatch> >();

// TODO: check that

// fBeamSpillTime(Start/End) and fFlashFindingTime(Start/End) are sane

if(fMinInTimeFlashes > fMaxFlashes){

throw cet::exception("FlashPredict")

<< "Minimum number of flashes fMinInTimeFlashes: " << fMinInTimeFlashes << ",\n"

<< "has to be at least equal to the maximum number of flashes fMaxFlashes: " << fMaxFlashes;

}

// TODO: check that all params are sane

if(fFlashStart > 0. || fFlashEnd < 0.){

throw cet::exception("FlashPredict")

<< "fFlashStart has to be non-positive, "

<< "and fFlashEnd has to be non-negative.";

}

if(std::abs(p.get<double>("DriftDistance")-driftDistance()) > 0.001){

mf::LogError("FlashPredict")

<< "Provided driftDistance: " << p.get<double>("DriftDistance")

<< " is different than expected: " << driftDistance();

}

if(fSBND && !fICARUS) {

if(fCryostat != 0) {

throw cet::exception("FlashPredict")

<< "SBND has only one cryostat. \n"

<< "Check Detector and Cryostat parameter." << std::endl;

}

}

else if(fICARUS && !fSBND) {

if(fCryostat > 1) {

throw cet::exception("FlashPredict")

<< "ICARUS has only two cryostats. \n"

<< "Check Detector and Cryostat parameter." << std::endl;

}

}

else {

throw cet::exception("FlashPredict")

<< "Detector: " << fDetector

<< ", not supported. Stopping.\n";

}

if(fIsSimple && (fOpFlashProducer.empty() || fOpFlashHitProducer.empty())) {

throw cet::exception("FlashPrediction")

<< "Require OpFlashProducer and OpFlashHitProducer to score OpFlashes" << std::endl;

}

if (((fOpHitTime != "RiseTime" && fOpHitTime != "PeakTime" &&

fOpHitTime != "StartTime")) ||

(fUseOpHitRiseTime + fUseOpHitPeakTime + fUseOpHitStartTime > 1) ||

(fUseOpHitRiseTime + fUseOpHitPeakTime + fUseOpHitStartTime <= 0))

throw cet::exception("FlashPredict")

<< "Supported OpHit times are 'RiseTime', 'PeakTime' or 'StartTime', "

<< "selected OpHitTime parameter is: " << fOpHitTime << std::endl;

if (!fMakeTree && fStoreMCInfo)

throw cet::exception("FlashPredict")

<< "Conflicting options:\n"

<< "MakeTree: " << std::boolalpha << fMakeTree << ", and StoreMCInfo: "

<< std::boolalpha << fStoreMCInfo << "\n"

<< "Theres no point to store MC info if not making the tree" << std::endl;

if (fMakeTree) initTree();

consumes<std::vector<recob::PFParticle>>(fPandoraProducer);

consumes<std::vector<recob::Slice>>(fPandoraProducer);

consumes<art::Assns<recob::SpacePoint, recob::PFParticle>>(fPandoraProducer);

consumes<std::vector<recob::SpacePoint>>(fSpacePointProducer);

consumes<art::Assns<recob::Hit, recob::SpacePoint>>(fSpacePointProducer);

consumes<std::vector<recob::OpHit>>(fOpHitProducer);

{

// sFM is an alias for sbn::SimpleFlashMatch

std::unique_ptr< std::vector<sFM> >

sFM_v(new std::vector<sFM>);

std::unique_ptr< art::Assns <recob::PFParticle, sFM> >

pfp_sFM_assn_v(new art::Assns<recob::PFParticle, sFM>);

// reset TTree variables

if(fMakeTree){

_evt = evt.event();

_sub = evt.subRun();

_run = evt.run();

_slices = 0;

_is_nu = -9999;

_flash_time = -9999.;

_flash_pe = -9999.;

_flash_unpe = -9999.;

_flash_rr = -9999.;

_flash_ratio = -9999.;

_hypo_x = -9999.;

_hypo_x_err = -9999.;

_hypo_x_rr = -9999.;

_hypo_x_ratio = -9999.;

_score = -9999.;

_mcT0 = -9999.;

}

bk.events++;

// LEGACY

// if(fMakeTree && fStoreTrueNus){

// // _true_nus = trueNus(evt);

// _true_nus = false;

// }

// grab PFParticles in event

const auto pfps_h =

evt.getValidHandle<std::vector<recob::PFParticle>>(fPandoraProducer);

if (pfps_h->size() == 0) {

mf::LogWarning("FlashPredict")

<< "No recob:PFParticle on event. Skipping...";

bk.noslice++;

updateBookKeeping();

for(size_t pId=0; pId<pfps_h->size(); pId++) {

if(!pfps_h->at(pId).IsPrimary()) continue;

const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNoPFPInEvt)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

}

evt.put(std::move(sFM_v));

evt.put(std::move(pfp_sFM_assn_v));

return;

}

if (fSelectNeutrino &&

!pfpNeutrinoOnEvent(pfps_h)) {

mf::LogInfo("FlashPredict")

<< "No pfp neutrino on event. Skipping...";

bk.nopfpneutrino++;

updateBookKeeping();

for(size_t pId=0; pId<pfps_h->size(); pId++) {

if(!pfps_h->at(pId).IsPrimary()) continue;

const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNoPFPInEvt)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

}

evt.put(std::move(sFM_v));

evt.put(std::move(pfp_sFM_assn_v));

return;

}

if(fMakeTree){

auto const& slice_h = evt.getValidHandle<std::vector<recob::Slice>>(fPandoraProducer);

_slices = slice_h.product()->size();

if (_slices == 0) {

mf::LogWarning("FlashPredict")

<< "No recob:Slice on event. Skipping...";

bk.noslice++;

updateBookKeeping();

for(size_t pId=0; pId<pfps_h->size(); pId++) {

if(!pfps_h->at(pId).IsPrimary()) continue;

const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNoSlcInEvt)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

}

evt.put(std::move(sFM_v));

evt.put(std::move(pfp_sFM_assn_v));

return;

}

}

// load OpHits previously created

std::vector<FlashMetrics> flashMetrics;

if(fIsSimple) { // Metrics for SimpleFlashes

std::vector<recob::OpHit> opHitsRght, opHitsLeft, opHits;

art::Handle<std::vector<recob::OpHit>> ophits_h;

evt.getByLabel(fOpHitProducer, ophits_h);

if(!ophits_h.isValid()) {

mf::LogError("FlashPredict")

<< "No optical hits from producer module "

<< fOpHitProducer;

bk.nonvalidophit++;

updateBookKeeping();

for(size_t pId=0; pId<pfps_h->size(); pId++) {

if(!pfps_h->at(pId).IsPrimary()) continue;

const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNoOpHInEvt)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

}

evt.put(std::move(sFM_v));

evt.put(std::move(pfp_sFM_assn_v));

return;

}

opHits.resize(ophits_h->size());

copyOpHitsInFlashFindingWindow(opHits, ophits_h);

mf::LogInfo("FlashPredict")

<< "OpHits found: " << opHits.size() << std::endl;

const std::vector<SimpleFlash> simpleFlashes = (fSBND) ?

makeSimpleFlashes(opHits, opHitsRght, opHitsLeft) : makeSimpleFlashes(opHits);

auto is_flash_in_time = [this](const SimpleFlash& f) -> bool

{ return (fBeamSpillTimeStart<=f.maxpeak_time &&

f.maxpeak_time<=fBeamSpillTimeEnd); };

auto flash_in_time = std::find_if(simpleFlashes.begin(), simpleFlashes.end(),

is_flash_in_time);

if(simpleFlashes.empty() ||

flash_in_time == simpleFlashes.end()){

mf::LogWarning("FlashPredict")

<< "No SimpleFlashes in beam window [" << fBeamSpillTimeStart << ", " << fBeamSpillTimeEnd << "], "

<< "\nor the sum of PE is less or equal to " << fMinFlashPE << " or 0."

<< "\nSkipping...";

bk.nullophittime++;

updateBookKeeping();

for(size_t pId=0; pId<pfps_h->size(); pId++) {

if(!pfps_h->at(pId).IsPrimary()) continue;

const art::Ptr<recob::PFParticle> pfp_ptr(pfps_h, pId);

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNoOpHInEvt)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

}

evt.put(std::move(sFM_v));

evt.put(std::move(pfp_sFM_assn_v));

return;

} else {

for(auto& sf : simpleFlashes) flashMetrics.push_back(getFlashMetrics(sf));

}

}

else

{ // Get Metrics from OpFlashes

for(unsigned i_tpc=0;i_tpc < fOpFlashProducer.size(); i_tpc++) {

art::Handle<std::vector<recob::OpFlash>> opflashes_h;

evt.getByLabel(fOpFlashProducer[i_tpc], opflashes_h);

if(opflashes_h->empty()) continue;

mf::LogInfo("FlashPredict")

<< "OpFlashes: " << opflashes_h->size() << std::endl;

// Check that there are OpFlashes

art::FindManyP<recob::OpHit> OpFlashToOpHitAssns(opflashes_h, evt, fOpFlashHitProducer[i_tpc]);

for(unsigned opf=0;opf<opflashes_h->size();opf++) {

auto& opflash = (*opflashes_h)[opf];

double optime = -9999.;

if(fSBND) optime = opflash.AbsTime();

else if(fICARUS) optime = opflash.Time();

if(optime<fFlashFindingTimeStart || optime>fFlashFindingTimeEnd) continue;

std::vector<art::Ptr<recob::OpHit>> ophit_v = OpFlashToOpHitAssns.at(opf);

flashMetrics.push_back(getFlashMetrics(opflash, ophit_v, opf));

}

}

}

mf::LogInfo("FlashPredict")

<< "FlashMetrics found: " << flashMetrics.size() << std::endl;

ChargeDigestMap chargeDigestMap = makeChargeDigest(evt, pfps_h);

for(auto& chargeDigestPair : chargeDigestMap) {

const auto& chargeDigest = chargeDigestPair.second;

const auto& pfp_ptr = chargeDigest.pfp_ptr;

const unsigned hitsInVolume = chargeDigest.hitsInVolume;

bk.pfp_to_score++;

if(chargeDigestPair.first < 0.){

mf::LogDebug("FlashPredict") << "Not a nu candidate slice. Skipping...";

bk.no_nu_candidate++;

mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNotANuScr)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

continue;

}

ChargeMetrics charge = computeChargeMetrics(chargeDigest);

if(!charge.metric_ok){

mf::LogWarning("FlashPredict")

<< "Clusters with No Charge.\n"

<< "the charge computed in the digest: " << chargeDigestPair.first

<< "\nSkipping...";

bk.no_charge++;

mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kNoChrgScr)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

continue;

}

FlashMetrics flash = {};

Score score = {std::numeric_limits<double>::max()};

bool hits_ophits_concurrence = false;

for(auto& origFlash : flashMetrics) {

unsigned ophsInVolume = origFlash.activity;

if(!isConcurrent(ophsInVolume, hitsInVolume)) continue;

hits_ophits_concurrence = true;

Score score_tmp = (fUse3DMetrics) ? computeScore3D(charge, origFlash) :

computeScore(charge, origFlash);

if(0. <= score_tmp.total && score_tmp.total < score.total

&& origFlash.metric_ok){

score = score_tmp;

flash = origFlash;

// // TODO: create charge.xb and/or charge.xb_gl

// if (fCorrectDriftDistance){

// charge.x = driftCorrection(charge.xb, flash.time);

// charge.x_gl = xGlCorrection(charge.x_gl, charge.xb, flash.time);

// }

}

} // for simpleFlashes

if(!hits_ophits_concurrence) {

std::string extra_message = (!fForceConcurrence) ? "" :

"\nConsider setting ForceConcurrence to false to lower requirements";

mf::LogInfo("FlashPredict")

<< "No OpHits where there's charge. Skipping..." << extra_message;

bk.no_oph_hits++;

mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(kQNoOpHScr)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

continue;

}

else if(!flash.metric_ok){

printMetrics("ERROR", charge, flash, 0, mf::LogError("FlashPredict"));

bk.no_flash_pe++;

mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";

sFM_v->emplace_back(sFM(kNoScr, kNoScrTime, Charge(kNoScrQ),

Flash(kNoScrPE), Score(k0VUVPEScr)));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

continue;

}

if(0. <= score.total &&

score.total < std::numeric_limits<double>::max()){

if(fMakeTree) {

_mcT0 = chargeDigest.mcT0;

_is_nu = chargeDigest.isNu;

_petoq = PEToQ(flash.pe, charge.q);

updateChargeMetrics(charge);

updateFlashMetrics(flash);

updateScore(score);

_flashmatch_nuslice_tree->Fill();

}

bk.scored_pfp++;

mf::LogDebug("FlashPredict") << "Creating sFM and PFP-sFM association";

Charge c{charge.q, TVector3(charge.x_gl, charge.y, charge.z),

TVector3(charge.x_glw, charge.yw, charge.zw)};

Flash f{flash.pe, TVector3(flash.x_gl, flash.y, flash.z),

TVector3(flash.xw, flash.yw, flash.zw)};

sFM_v->emplace_back(sFM(true, flash.time, c, f, score));

util::CreateAssn(*this, evt, *sFM_v, pfp_ptr, *pfp_sFM_assn_v);

}

else{

mf::LogError("FlashPredict")

<< "ERROR: score < 0. Dumping info.\n"

<< "score: " << score.total << "\n"

<< "score.y: " << score.y << "\t"

<< "score.z: " << score.z << "\n"

<< "score.rr: " << score.rr << "\t"

<< "score.ratio: " << score.ratio << "\n"

<< "score.slope: " << score.slope << "\t"

<< "score.petoq: " << score.petoq << "\n";

printMetrics("ERROR", charge, flash, 0, mf::LogError("FlashPredict"));

}

} // chargeDigestMap: PFparticles that pass criteria

bk.events_processed++;

updateBookKeeping();

evt.put(std::move(sFM_v));

evt.put(std::move(pfp_sFM_assn_v));

{

art::ServiceHandle<art::TFileService> tfs;

_flashmatch_nuslice_tree = tfs->make<TTree>("nuslicetree", "nu FlashPredict tree");

_flashmatch_nuslice_tree->Branch("evt", &_evt, "evt/I");

_flashmatch_nuslice_tree->Branch("run", &_run, "run/I");

_flashmatch_nuslice_tree->Branch("sub", &_sub, "sub/I");

_flashmatch_nuslice_tree->Branch("slices", &_slices, "slices/I");

_flashmatch_nuslice_tree->Branch("flash_id", &_flash_id, "flash_id/I");

_flashmatch_nuslice_tree->Branch("flash_activity", &_flash_activity, "flash_activity/I");

_flashmatch_nuslice_tree->Branch("flash_x", &_flash_x, "flash_x/D");

_flashmatch_nuslice_tree->Branch("flash_yb", &_flash_yb, "flash_yb/D");

_flashmatch_nuslice_tree->Branch("flash_zb", &_flash_zb, "flash_zb/D");

_flashmatch_nuslice_tree->Branch("flash_x_gl", &_flash_x_gl, "flash_x_gl/D");

_flashmatch_nuslice_tree->Branch("flash_y", &_flash_y, "flash_y/D");

_flashmatch_nuslice_tree->Branch("flash_z", &_flash_z, "flash_z/D");

_flashmatch_nuslice_tree->Branch("flash_xw", &_flash_xw, "flash_xw/D");

_flashmatch_nuslice_tree->Branch("flash_yw", &_flash_yw, "flash_yw/D");

_flashmatch_nuslice_tree->Branch("flash_zw", &_flash_zw, "flash_zw/D");

_flashmatch_nuslice_tree->Branch("flash_rr", &_flash_rr, "flash_rr/D");

_flashmatch_nuslice_tree->Branch("flash_ratio", &_flash_ratio, "flash_ratio/D");

_flashmatch_nuslice_tree->Branch("flash_slope", &_flash_slope, "flash_slope/D");

_flashmatch_nuslice_tree->Branch("flash_pe", &_flash_pe, "flash_pe/D");

_flashmatch_nuslice_tree->Branch("flash_unpe", &_flash_unpe, "flash_unpe/D");

_flashmatch_nuslice_tree->Branch("flash_time", &_flash_time, "flash_time/D");

_flashmatch_nuslice_tree->Branch("hypo_x", &_hypo_x, "hypo_x/D");

_flashmatch_nuslice_tree->Branch("hypo_x_err", &_hypo_x_err, "hypo_x_err/D");

_flashmatch_nuslice_tree->Branch("hypo_x_rr", &_hypo_x_rr, "hypo_x_rr/D");

_flashmatch_nuslice_tree->Branch("hypo_x_ratio", &_hypo_x_ratio, "hypo_x_ratio/D");

_flashmatch_nuslice_tree->Branch("y_skew", &_y_skew, "y_skew/D");

_flashmatch_nuslice_tree->Branch("z_skew", &_z_skew, "z_skew/D");

_flashmatch_nuslice_tree->Branch("y_kurt", &_y_kurt, "y_kurt/D");

_flashmatch_nuslice_tree->Branch("z_kurt", &_z_kurt, "z_kurt/D");

_flashmatch_nuslice_tree->Branch("charge_id", &_charge_id, "charge_id/I");

_flashmatch_nuslice_tree->Branch("charge_activity", &_charge_activity, "charge_activity/I");

_flashmatch_nuslice_tree->Branch("charge_x", &_charge_x, "charge_x/D");

_flashmatch_nuslice_tree->Branch("charge_x_gl", &_charge_x_gl, "charge_x_gl/D");

_flashmatch_nuslice_tree->Branch("charge_y", &_charge_y, "charge_y/D");

_flashmatch_nuslice_tree->Branch("charge_z", &_charge_z, "charge_z/D");

_flashmatch_nuslice_tree->Branch("charge_x_glw", &_charge_x_glw, "charge_xglw/D");

_flashmatch_nuslice_tree->Branch("charge_yw", &_charge_yw, "charge_yw/D");

_flashmatch_nuslice_tree->Branch("charge_zw", &_charge_zw, "charge_zw/D");

_flashmatch_nuslice_tree->Branch("charge_slope", &_charge_slope, "charge_slope/D");

_flashmatch_nuslice_tree->Branch("charge_q", &_charge_q, "charge_q/D");

_flashmatch_nuslice_tree->Branch("petoq", &_petoq, "petoq/D");

_flashmatch_nuslice_tree->Branch("score", &_score, "score/D");

_flashmatch_nuslice_tree->Branch("scr_y", &_scr_y, "scr_y/D");

_flashmatch_nuslice_tree->Branch("scr_z", &_scr_z, "scr_z/D");

_flashmatch_nuslice_tree->Branch("scr_rr", &_scr_rr, "scr_rr/D");

_flashmatch_nuslice_tree->Branch("scr_ratio", &_scr_ratio, "scr_ratio/D");

_flashmatch_nuslice_tree->Branch("scr_slope", &_scr_slope, "scr_slope/D");

_flashmatch_nuslice_tree->Branch("scr_petoq", &_scr_petoq, "scr_petoq/D");

_flashmatch_nuslice_tree->Branch("is_nu", &_is_nu, "is_nu/I");

_flashmatch_nuslice_tree->Branch("mcT0", &_mcT0, "mcT0/D");

const std::string inputFilename) const

{

// Many changes needed everywhere

ReferenceMetrics rm;

// read histograms and fill vectors for match score calculation

std::string fname;

cet::search_path sp("FW_SEARCH_PATH");

if(!sp.find_file(inputFilename, fname)) {

mf::LogError("FlashPredict")

<< "Could not find the light-charge match root file '"

<< inputFilename << "' on FW_SEARCH_PATH\n";

throw cet::exception("FlashPredict")

<< "Could not find the light-charge match root file '"

<< inputFilename << "' on FW_SEARCH_PATH\n";

}

TFile *topfile = new TFile(fname.c_str(), "READ");

auto dirsinfile = topfile->GetListOfKeys();

if(!dirsinfile->Contains(fFlashType.c_str()) && !fUseOldMetrics) {

throw cet::exception("FlashPredict")

<< "Metrics file " << inputFilename

<< " doesn't contain TDirectoryFile " << fFlashType << "\n"

<< "Set UseOldMetrics to true or update FW_SEARCH_PATH \n";

}

else if(dirsinfile->Contains(fFlashType.c_str()) && fUseOldMetrics) {

throw cet::exception("FlashPredict")

<< "Metrics file " << inputFilename

<< " contains TDirectoryFile " << fFlashType << " but UseOldMetrics is set to true \n"

<< "Set UseOldMetrics to false or update FW_SEARCH_PATH \n";

}

TDirectory *infile = (fUseOldMetrics) ? topfile : (TDirectory*)topfile->Get(fFlashType.c_str());

auto metricsInFile = infile->GetListOfKeys();

if(!metricsInFile->Contains("dy_h1") ||

!metricsInFile->Contains("dz_h1") ||

!metricsInFile->Contains("rr_h1") ||

!metricsInFile->Contains("ratio_h1") ||

!metricsInFile->Contains("slope_h1") ||

!metricsInFile->Contains("petoq_h1") ||

// LEGACY

// !metricsInFile->Contains("rr_fit_l") ||

// !metricsInFile->Contains("rr_fit_m") ||

// !metricsInFile->Contains("rr_fit_h") ||

// !metricsInFile->Contains("ratio_fit_l") ||

// !metricsInFile->Contains("ratio_fit_m") ||

// !metricsInFile->Contains("ratio_fit_h") ||

// LEGACY

!metricsInFile->Contains("pol_coeffs_y") ||

!metricsInFile->Contains("pol_coeffs_z"))

{

mf::LogError("FlashPredict")

<< "The metrics file '" << fname << "' lacks at least one metric.";

throw cet::exception("FlashPredict")

<< "The metrics file '" << fname << "'lacks at least one metric.";

}

//

mf::LogInfo("FlashPredict")

<< "Checked all metrics are present" << std::endl;

TH1 *temphisto = (TH1*)infile->Get("dy_h1");

mf::LogInfo("FlashPredict")

<< "Loaded dy" << std::endl;

int bins = temphisto->GetNbinsX();

for (int ib = 1; ib <= bins; ++ib) {

rm.dYMeans.push_back(temphisto->GetBinContent(ib));

rm.dYSpreads.push_back(temphisto->GetBinError(ib));

}

//

mf::LogInfo("FlashPredict")

<< "Read dy hist" << std::endl;

temphisto = (TH1*)infile->Get("dz_h1");

bins = temphisto->GetNbinsX();

for (int ib = 1; ib <= bins; ++ib) {

rm.dZMeans.push_back(temphisto->GetBinContent(ib));

rm.dZSpreads.push_back(temphisto->GetBinError(ib));

}

//

temphisto = (TH1*)infile->Get("rr_h1");

bins = temphisto->GetNbinsX();

for (int ib = 1; ib <= bins; ++ib) {

rm.RRMeans.push_back(temphisto->GetBinContent(ib));

rm.RRSpreads.push_back(temphisto->GetBinError(ib));

}

// LEGACY, now broken

// unsigned s = 0;

// for(auto& rrF : rm.RRFits){

// std::string nold = "rr_fit_" + kSuffixes[s];

// std::string nnew = "rrFit_" + kSuffixes[s];

// TF1* tempF1 = (TF1*)infile->Get(nold.c_str());

// auto params = tempF1->GetParameters();

// rrF.f = std::make_unique<TF1>(nnew.c_str(), fRR_TF1_fit.c_str(),

// 0., fDriftDistance);

// rrF.f->SetParameters(params);

// rrF.min = rrF.f->GetMinimum(0., fDriftDistance, kEps);

// rrF.max = rrF.f->GetMaximum(0., fDriftDistance, kEps);

// s++;

// }

//

temphisto = (TH1*)infile->Get("ratio_h1");

bins = temphisto->GetNbinsX();

for (int ib = 1; ib <= bins; ++ib) {

rm.RatioMeans.push_back(temphisto->GetBinContent(ib));

rm.RatioSpreads.push_back(temphisto->GetBinError(ib));

}

// LEGACY, now broken

// s = 0;

// for(auto& ratioF : rm.RatioFits){

// std::string nold = "ratio_fit_" + kSuffixes[s];

// std::string nnew = "ratioFit_" + kSuffixes[s];

// TF1* tempF1 = (TF1*)infile->Get(nold.c_str());

// auto params = tempF1->GetParameters();

// ratioF.f = std::make_unique<TF1>(nnew.c_str(), fRatio_TF1_fit.c_str(),

// 0., fDriftDistance);

// ratioF.f->SetParameters(params);

// ratioF.min = ratioF.f->GetMinimum(0., fDriftDistance, kEps);

// ratioF.max = ratioF.f->GetMaximum(0., fDriftDistance, kEps);

// s++;

// }

//

temphisto = (TH1*)infile->Get("slope_h1");

bins = temphisto->GetNbinsX();

for (int ib = 1; ib <= bins; ++ib) {

rm.SlopeMeans.push_back(temphisto->GetBinContent(ib));

rm.SlopeSpreads.push_back(temphisto->GetBinError(ib));

}

//

temphisto = (TH1*)infile->Get("petoq_h1");

bins = temphisto->GetNbinsX();

for (int ib = 1; ib <= bins; ++ib) {

rm.PEToQMeans.push_back(temphisto->GetBinContent(ib));

rm.PEToQSpreads.push_back(temphisto->GetBinError(ib));

}

TH2* tmp1_h2 = (TH2*)infile->Get("rr_h2");

rm.RRH2 = (TH2D*)tmp1_h2->Clone("RRH2");

TH2* tmp2_h2 = (TH2*)infile->Get("ratio_h2");

rm.RatioH2 = (TH2D*)tmp2_h2->Clone("RatioH2");

std::vector<double>* polCoeffsY_p;

infile->GetObject("pol_coeffs_y", polCoeffsY_p);

rm.PolCoeffsY = std::vector(*polCoeffsY_p);

std::string tty = "Polynomial correction coefficients for Y: ( ";

for(double c : rm.PolCoeffsY) tty += std::to_string(c) + " ";

tty += ")";

mf::LogInfo("FlashPredict") << tty;

std::vector<double>* polCoeffsZ_p;

infile->GetObject("pol_coeffs_z", polCoeffsZ_p);

rm.PolCoeffsZ = std::vector(*polCoeffsZ_p);

std::string ttz = "Polynomial correction coefficients for Z: ( ";

for(double c : rm.PolCoeffsZ) ttz += std::to_string(c) + " ";

ttz += ")";

mf::LogInfo("FlashPredict") << ttz;

// BIG TODO: Metrics should depend on X,Y,Z.

// TODO: Test!

// TODO: store 3D-arrays of means and spreads, instead of the

// TProfile3D

TProfile3D* tmp1_prof3 = (TProfile3D*)infile->Get("dy_prof3");

rm.dYP3 = (TProfile3D*)tmp1_prof3->Clone("dYP3");

TProfile3D* tmp2_prof3 = (TProfile3D*)infile->Get("dz_prof3");

rm.dZP3 = (TProfile3D*)tmp2_prof3->Clone("dZP3");

TProfile3D* tmp3_prof3 = (TProfile3D*)infile->Get("rr_prof3");

rm.RRP3 = (TProfile3D*)tmp3_prof3->Clone("RRP3");

TProfile3D* tmp4_prof3 = (TProfile3D*)infile->Get("ratio_prof3");

rm.RatioP3 = (TProfile3D*)tmp4_prof3->Clone("RatioP3");

TProfile3D* tmp5_prof3 = (TProfile3D*)infile->Get("slope_prof3");

rm.SlopeP3 = (TProfile3D*)tmp5_prof3->Clone("SlopeP3");

TProfile3D* tmp6_prof3 = (TProfile3D*)infile->Get("petoq_prof3");

rm.PEToQP3 = (TProfile3D*)tmp6_prof3->Clone("PEToQP3");

infile->Close();

delete infile;

mf::LogInfo("FlashPredict") << "Finish loading metrics";

return rm;

const std::vector<art::Ptr<recob::Hit>>& hits,

const std::vector<art::Ptr<simb::MCParticle>>& mcParticles) const

{

auto clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();

int pidMaxEnergy =

TruthMatchUtils::TrueParticleIDFromTotalTrueEnergy(clockData, hits, true);

for(auto& mcp: mcParticles){

if(mcp->TrackId() == pidMaxEnergy){

double mcT0 = mcp->Position().T()/1000.;

art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;

art::Ptr<simb::MCTruth> truth = pi_serv->TrackIdToMCTruth_P(pidMaxEnergy);

bool isNu = (truth->Origin() == simb::kBeamNeutrino);

return {mcT0, isNu};

}

}

return {-9999.999, false};

const ChargeDigest& chargeDigest) const

{

const unsigned pId = chargeDigest.pId;

const int pfpPDGC = chargeDigest.pfpPDGC;

const unsigned hitsInVolume = chargeDigest.hitsInVolume;

const auto& qClusters = chargeDigest.qClusters;

using Q_t = const flashmatch::QPoint_t&;

double charge_q = std::accumulate(qClusters.begin(), qClusters.end(), 0.,

[](double x, Q_t q) {return x+q.q;});

if (charge_q <= 0.) return {};

ChargeMetrics charge;

charge.id = pId; charge.activity = hitsInVolume; charge.pdgc = pfpPDGC;

charge.q = charge_q; charge.metric_ok = true;

charge.x_gl =

(std::accumulate(qClusters.begin(), qClusters.end(), 0.,

[](double x, Q_t q) {return x+q.q*q.x;}))/charge_q;

charge.x = foldXGl(charge.x_gl);

charge.y =

(std::accumulate(qClusters.begin(), qClusters.end(), 0.,

[](double y, Q_t q) {return y+q.q*q.y;}))/charge_q;

charge.z =

(std::accumulate(qClusters.begin(), qClusters.end(), 0.,

[](double z, Q_t q) {return z+q.q*q.z;}))/charge_q;

unsigned charge_center_in_volume = fGeometry->PositionToTPCID(

{charge.x_gl, charge.y, charge.z}).TPC/fTPCPerDriftVolume;

if (charge_center_in_volume>2 && fSBND){

// TODO HACK

charge_center_in_volume = fGeometry->PositionToTPCID(

{charge.x_gl/2., charge.y, charge.z}).TPC/fTPCPerDriftVolume;

}

unsigned activity = 0;

if (charge_center_in_volume == kRght) activity = kActivityInRght;

else if (charge_center_in_volume == kLeft) activity = kActivityInLeft;

if (charge.activity < kActivityInBoth &&

(activity != charge.activity) && fVerbose){

mf::LogError("FlashPredict")

<< "INFO: The charge center: ("

<< charge.x_gl << ", " << charge.y << ", " << charge.z << "), "

<< "belonging to volume " << charge_center_in_volume << "\n"

<< "is not inside the volume of the chargeDigest activity "

<< charge.activity << ".\n"

<< "Will update the charge activity to have both volumes.";

charge.activity = kActivityInBoth;

}

// Compute the widths

// TODO: unharcode... but these numbers work for SBND and ICARUS

std::unique_ptr<TH1F> qX = std::make_unique<TH1F>("qX", "", 160, -400., 400.);

std::unique_ptr<TH1F> qY = std::make_unique<TH1F>("qY", "", 84, -210., 210.);

std::unique_ptr<TH1F> qZ = std::make_unique<TH1F>("qZ", "", 364, -910., 910.);

for (size_t i=0; i<qClusters.size(); ++i) {

// double q2 = qClusters[i].q * qClusters[i].q;

double q = qClusters[i].q;

qX->Fill(qClusters[i].x, q);

qY->Fill(qClusters[i].y, q);

qZ->Fill(qClusters[i].z, q);

}

charge.x_glw = qX->GetStdDev();

charge.yw = qY->GetStdDev();

charge.zw = qZ->GetStdDev();

// Now fractional widths

double csize = std::sqrt(charge.x_glw*charge.x_glw + charge.yw*charge.yw + charge.zw*charge.zw);

// double cx_fac = charge.x_glw / csize;

double cy_fac = charge.yw / csize;

double cz_fac = charge.zw / csize;

double min_axis = std::min(cy_fac, cz_fac);

double max_axis = std::max(cy_fac, cz_fac);

double dir = ((cz_fac>=cy_fac)) ? 1. : -1.;

charge.slope = dir * std::sqrt((max_axis*max_axis) - (min_axis*min_axis));

// Store fractional widths?

// charge.x_glw = cx_fac;

// charge.yw = cy_fac;

// charge.zw = cz_fac;

return charge;

const std::vector<recob::OpHit>& ophits) const

{

std::unique_ptr<TH1F> ophY = std::make_unique<TH1F>("ophY", "", fYBins, fYLow, fYHigh);

std::unique_ptr<TH1F> ophZ = std::make_unique<TH1F>("ophZ", "", fZBins, fZLow, fZHigh);

std::unique_ptr<TH1F> oph2Y = std::make_unique<TH1F>("oph2Y", "", fYBins, fYLow, fYHigh);

std::unique_ptr<TH1F> oph2Z = std::make_unique<TH1F>("oph2Z", "", fZBins, fZLow, fZHigh);

double peSumMax_wallX = wallXWithMaxPE(ophits);

double sum = 0.;

double sum_PE = 0.;

double sum_PE2 = 0.;

double sum_unPE = 0.;

double sum_PE2Y = 0.; double sum_PE2Z = 0.;

double sum_PE2Y2 = 0.; double sum_PE2Z2 = 0.;

for(auto& oph : ophits) {

int opChannel = oph.OpChannel();

auto opDetXYZ = fWireReadoutGeom->OpDetGeoFromOpChannel(opChannel).GetCenter();

bool is_pmt_vis = false, is_ara_vis = false;

if(fSBND){// because VIS light

auto op_type = fPDMapAlgPtr->pdType(opChannel);

if(op_type == "pmt_uncoated") {

is_pmt_vis = true, is_ara_vis = false;

}

else if(op_type == "xarapuca_vis" || op_type == "arapuca_vis") {

is_pmt_vis = false, is_ara_vis = true;

}

}

double ophPE = oph.PE();

double ophPE2 = ophPE * ophPE;

sum += 1.0;

sum_PE += ophPE;

sum_PE2 += ophPE2;

sum_PE2Y += ophPE2 * opDetXYZ.Y();

sum_PE2Z += ophPE2 * opDetXYZ.Z();

sum_PE2Y2 += ophPE2 * opDetXYZ.Y() * opDetXYZ.Y();

sum_PE2Z2 += ophPE2 * opDetXYZ.Z() * opDetXYZ.Z();

ophY->Fill(opDetXYZ.Y(), ophPE);

ophZ->Fill(opDetXYZ.Z(), ophPE);

oph2Y->Fill(opDetXYZ.Y(), ophPE2);

oph2Z->Fill(opDetXYZ.Z(), ophPE2);

if(fICARUS){

if(std::abs(peSumMax_wallX-opDetXYZ.X()) > 5.) sum_unPE += ophPE;

}

else {// fSBND

if(is_pmt_vis || is_ara_vis) {

sum_unPE += ophPE;

}

}

} // for opHits

if (sum_PE > 0.) {

FlashMetrics flash;

flash.metric_ok = true;

flash.pe = sum_PE;

flash.unpe = sum_unPE;

flash.y_skew = ophY->GetSkewness();

flash.z_skew = ophZ->GetSkewness();

flash.y_kurt = ophY->GetKurtosis();

flash.z_kurt = ophZ->GetKurtosis();

// Flash widths

// flash.xw = fractTimeWithFractionOfLight(orig_flash, flash.pe, fFlashPEFraction);

flash.xw = fractTimeWithFractionOfLight(ophits, sum_PE2, fFlashPEFraction, true);

// flash.xw = fractTimeWithFractionOfLight(orig_flash, flash.unpe, fFlashPEFraction, false, true); // TODO:

// TODO: low values of flash.xw <0.5 are indicative of good

// mcT0-flash_time matching, so akin to matching to prompt light

// Note that around the middle of the detector (~(+-100, 0, 250)

// cm for SBND) the values of flash.xw are slightly larger, this

// is natural and has to do with the way light disperses

flash.yw = oph2Y->GetStdDev();

flash.zw = oph2Z->GetStdDev();

flash.ratio = fOpDetNormalizer * flash.unpe / flash.pe;

if(fSBND && (fFlashType == "simpleflash_ara" || fFlashType == "opflash_ara")) {

flash.ratio = flash.unpe / (flash.pe + flash.unpe);

}

flash.yb = sum_PE2Y / sum_PE2;

flash.zb = sum_PE2Z / sum_PE2;

flash.rr = std::sqrt(

std::abs(sum_PE2Y2 + sum_PE2Z2 + sum_PE2 * (flash.yb * flash.yb + flash.zb * flash.zb)

- 2.0 * (flash.yb * sum_PE2Y + flash.zb * sum_PE2Z) ) / sum_PE2);

std::tie(flash.h_x, flash.h_xerr, flash.h_xrr, flash.h_xratio) =

hypoFlashX_H2(flash.rr, flash.ratio);

// TODO: using _hypo_x make further corrections to _flash_time to

// account for light transport time and/or rising edge

// double flash_time = timeCorrections(orig_flash.maxpeak_time, hypo_x);

flash.x = peSumMax_wallX;

flash.x_gl = flashXGl(flash.h_x, flash.x);

// Fractional widths

double fsize = std::sqrt(flash.yw*flash.yw + flash.zw*flash.zw);

double fy_fac = flash.yw / fsize;

double fz_fac = flash.zw / fsize;

double min_axis = std::min(fy_fac, fz_fac);

double max_axis = std::max(fy_fac, fz_fac);

double dir = (fz_fac>=fy_fac) ? 1. : -1.;

flash.slope = dir * std::sqrt((max_axis*max_axis) - (min_axis*min_axis));

// store fractional widths?

// flash.yw = fy_fac;

// flash.zw = fz_fac;

flash.y = flash.yb - polynomialCorrection(flash.y_skew, flash.h_x,

fRM.PolCoeffsY, fSkewLimitY);

flash.z = flash.zb - polynomialCorrection(flash.z_skew, flash.h_x,

fRM.PolCoeffsZ, fSkewLimitZ);

return flash;

}

else {

std::string channels;

for(auto oph : ophits) channels += std::to_string(oph.OpChannel()) + ' ';

mf::LogError("FlashPredict")

<< "Really odd that I landed here, this shouldn't had happen.\n"

<< "sum: \t" << sum << "\n"

<< "sum_PE: \t" << sum_PE << "\n"

<< "sum_unPE: \t" << sum_unPE << "\n"

<< "opHits size: \t" << ophits.size() << "\n"

<< "channels: \t" << channels << std::endl;

return {};

}

const FlashPredict::SimpleFlash& sf) const

{

std::vector<recob::OpHit> ophits(sf.opH_beg, sf.opH_end);

std::sort(ophits.begin(), ophits.end(),

[this] (const recob::OpHit& oph1, const recob::OpHit& oph2)

{ return (opHitTime(oph1) < opHitTime(oph2)); });

auto fm = computeFlashMetrics(ophits);

fm.id = sf.flashId;

fm.activity = sf.ophsInVolume;

fm.time = opHitTime(ophits[0]);

return fm;

const recob::OpFlash& opflash,

std::vector<art::Ptr<recob::OpHit>> ophit_v,

unsigned id) const

{

std::vector<recob::OpHit> ophits;

for(unsigned i=0; i< ophit_v.size(); i++) {

ophits.emplace_back(*(ophit_v)[i]);

}

std::sort(ophits.begin(), ophits.end(),

[this] (const recob::OpHit& oph1, const recob::OpHit& oph2)

{ return (opHitTime(oph1) < opHitTime(oph2)); });

// Search for ophits in either TPC. X-coordinate is flipped by convention

auto fm = computeFlashMetrics(ophits);

fm.id = id;

if(fSBND) fm.time = opflash.AbsTime();

else if(fICARUS) fm.time = opflash.Time();

bool in_left = false, in_right = false;

if(fSBND) {

for(auto const& oph : ophits) {

auto ch = oph.OpChannel();

auto opDetX = fWireReadoutGeom->OpDetGeoFromOpChannel(ch).GetCenter().X();

if(opDetX >= 0.) in_left = true;

else in_right = true;

}

}

else if(fICARUS) {

for(auto const& oph : ophits) {

auto ch = oph.OpChannel();

auto opDetXYZ = fWireReadoutGeom->OpDetGeoFromOpChannel(ch).GetCenter();

if(!fGeoCryo->ContainsPosition(opDetXYZ)) continue;

unsigned t = icarusPDinTPC(ch);

if(t/fTPCPerDriftVolume == kRght) in_right = true;

else if(t/fTPCPerDriftVolume == kLeft) in_left = true;

}

}

if(in_left) {

fm.activity = (in_right) ? kActivityInBoth : kActivityInLeft;

}

else fm.activity = kActivityInRght;

if(fUseOpCoords) {

if(opflash.hasXCenter()) {

fm.x = opflash.XCenter();

fm.xw = opflash.XWidth();

}

fm.y = opflash.YCenter();

fm.yw = opflash.YWidth();

fm.z = opflash.ZCenter();

fm.zw = opflash.ZWidth();

}