CDCDedxElectronCollectorModule.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#include <cdc/modules/CDCDedxElectronCollector/CDCDedxElectronCollectorModule.h>
 
#include <mdst/dataobjects/ECLCluster.h>
#include <TTree.h>
#include <TH1D.h>
#include <TH1I.h>
#include <TMath.h>
 
#include <cmath>
 
using namespace Belle2;
//-----------------------------------------------------------------
//                 Register the Module
//-----------------------------------------------------------------
REG_MODULE(CDCDedxElectronCollector);
 
//-----------------------------------------------------------------
//                 Implementation
//-----------------------------------------------------------------
CDCDedxElectronCollectorModule::CDCDedxElectronCollectorModule() : CalibrationCollectorModule()
{
  // Set module properties
  setDescription("A collector module for CDC dE/dx electron calibrations");
 
  // Parameter definitions
 
  addParam("cleanupCuts", m_cuts, "boolean to apply cleanup cuts", true);
  addParam("maxHits", m_maxHits, "maximum number of hits per track ", int(100));
  addParam("setEoP", m_setEoP, "Set E over p Cut values. ", double(0.25));
  addParam("isCosth", m_isCosth, "true for adding costh tree branch. ", false);
  addParam("isMom", m_isMom, "true for adding momentum tree branch. ", false);
  addParam("isPt", m_isPt, "true for adding transverse momentum tree branch. ", false);
  addParam("isCharge", m_isCharge, "true for charge dedx tree branch. ", false);
  addParam("isRun", m_isRun, "true for adding run number tree branch. ", false);
  addParam("isWire", m_isWire, "true for adding wires tree branch. ", false);
  addParam("isLayer", m_isLayer, "true for adding layers tree branch. ", false);
  addParam("isDoca", m_isDoca, "true for adding doca tree branch. ", false);
  addParam("isEnta", m_isEnta, "true for adding enta tree branch. ", false);
  addParam("isInjTime", m_isInjTime, "true for adding time var tree branch. ", false);
  addParam("isDocaRS", m_isDocaRS, "true for adding doca tree branch. ", false);
  addParam("isEntaRS", m_isEntaRS, "true for adding enta tree branch. ", false);
  addParam("isDedxhit", m_isDedxhit, "true for adding dedxhit tree branch. ", false);
  addParam("isADCcorr", m_isADCcorr, "true for adding adc tree branch. ", false);
  addParam("islLayer", m_islLayer, "true for adding layers for ldedx tree branch. ", false);
  addParam("islDedx", m_islDedx, "true for adding layer dedx tree branch. ", false);
  addParam("isBhabha", m_isBhabha, "true for bhabha events", true);
  addParam("isRadee", m_isRadee, "true for radee events", false);
  addParam("isTrgSel", m_isTrgSel, "true to enable trigger sel inside module", false);
}
 
//-----------------------------------------------------------------
//                 Create ROOT objects
//-----------------------------------------------------------------
 
void CDCDedxElectronCollectorModule::prepare()
{
  m_trgResult.isOptional();
  m_dedxTracks.isRequired();
  m_tracks.isRequired();
  m_trackFitResults.isRequired();
 
  // Data object creation
  auto means = new TH1D("means", "CDC dE/dx truncated means", 250, 0, 2.5);
  auto ttree = new TTree("tree", "Tree with dE/dx information");
 
  auto hestats = new TH1I("hestats", "Event Stats", 6, -0.5, 5.5);
  hestats->SetFillColor(kYellow);
  hestats->GetXaxis()->SetBinLabel(1, "all");
  hestats->GetXaxis()->SetBinLabel(2, "notrig");
  hestats->GetXaxis()->SetBinLabel(3, "noskim");
  hestats->GetXaxis()->SetBinLabel(4, "wrongskim");
  hestats->GetXaxis()->SetBinLabel(5, "unclean");
  hestats->GetXaxis()->SetBinLabel(6, "selected");
 
  auto htstats = new TH1I("htstats", "track Stats", 6, -0.5, 5.5);
  htstats->SetFillColor(kYellow);
  htstats->GetXaxis()->SetBinLabel(1, "alltrk");
  htstats->GetXaxis()->SetBinLabel(2, "vtx");
  htstats->GetXaxis()->SetBinLabel(3, "inCDC");
  htstats->GetXaxis()->SetBinLabel(4, "whits");
  htstats->GetXaxis()->SetBinLabel(5, "weop");
  htstats->GetXaxis()->SetBinLabel(6, "selected");
 
  if (m_isInjTime) {
    ttree->Branch<double>("injtime", &m_injTime);
    ttree->Branch<double>("injring", &m_injRing);
    ttree->Branch<double>("costh", &m_costh);
    ttree->Branch<int>("nhits", &m_nhits);
    ttree->Branch<double>("p", &m_p);
  }
 
  ttree->Branch<double>("dedx", &m_dedx);
  if (m_isCosth)ttree->Branch<double>("costh", &m_costh);
  if (m_isMom)ttree->Branch<double>("p", &m_p);
  if (m_isPt)ttree->Branch<double>("pt", &m_pt);
  if (m_isCharge)ttree->Branch<int>("charge", &m_charge);
  if (m_isRun)ttree->Branch<int>("run", &m_run);
 
  if (m_isWire)ttree->Branch("wire", &m_wire);
  if (m_isLayer)ttree->Branch("layer", &m_layer);
  if (m_isDoca)ttree->Branch("doca", &m_doca);
  if (m_isEnta)ttree->Branch("enta", &m_enta);
  if (m_isDocaRS)ttree->Branch("docaRS", &m_docaRS);
  if (m_isEntaRS)ttree->Branch("entaRS", &m_entaRS);
  if (m_isDedxhit)ttree->Branch("dedxhit", &m_dedxhit);
  if (m_isADCcorr)ttree->Branch("adccorr", &m_adccorr);
 
  if (m_islLayer)ttree->Branch("lLayer", &m_lLayer);
  if (m_islDedx)ttree->Branch("ldedx", &m_ldedx);
 
  // Collector object registration
  registerObject<TH1D>("means", means);
  registerObject<TTree>("tree", ttree);
  registerObject<TH1I>("hestats", hestats);
  registerObject<TH1I>("htstats", htstats);
}
 
//-----------------------------------------------------------------
//                 Fill ROOT objects
//-----------------------------------------------------------------
void CDCDedxElectronCollectorModule::collect()
{
 
  auto hestats = getObjectPtr<TH1I>("hestats");
  hestats->Fill(0);
 
  if (m_isTrgSel) {
    if (!m_trgResult.isValid()) {
      B2WARNING("SoftwareTriggerResult required to select bhabha/radee event is not found");
      hestats->Fill(1);
      return;
    }
 
    //release05: bhabha_all is grand skim = bhabha+bhabhaecl+radee
    const std::map<std::string, int>& fresults = m_trgResult->getResults();
    if (fresults.find("software_trigger_cut&skim&accept_bhabha") == fresults.end() and
        (fresults.find("software_trigger_cut&skim&accept_radee") == fresults.end()
         || fresults.find("software_trigger_cut&skim&accept_bhabha_cdc") == fresults.end())) {
      B2WARNING("Can't find required bhabha/radee trigger identifiers");
      hestats->Fill(2);
      return;
    }
 
    const bool eBhabha = (m_trgResult->getResult("software_trigger_cut&skim&accept_bhabha") ==
                          SoftwareTriggerCutResult::c_accept);
 
    bool eRadBhabha = false;
 
    const std::string radee =
      "software_trigger_cut&skim&accept_radee";
 
    const std::string bhabhaCDC =
      "software_trigger_cut&skim&accept_bhabha_cdc";
 
    if (fresults.find(bhabhaCDC) != fresults.end()) {
      // If bhabha_cdc exists, use only bhabha_cdc
      eRadBhabha = (m_trgResult->getResult(bhabhaCDC) == SoftwareTriggerCutResult::c_accept);
    } else if (fresults.find(radee) != fresults.end()) {
      // Otherwise use radee
      eRadBhabha = (m_trgResult->getResult(radee) == SoftwareTriggerCutResult::c_accept);
    }
 
 
    if (!m_isBhabha && !m_isRadee) {
      B2WARNING("requested not-supported event type: going back");
      hestats->Fill(3);
      return;
    } else if (m_isBhabha && !m_isRadee && !eBhabha) {
      B2WARNING("requested bhabha only but event not found: going back");
      hestats->Fill(3);
      return;
    } else  if (m_isRadee && !m_isBhabha && !eRadBhabha) {
      B2WARNING("requested radee only but event not found: going back");
      hestats->Fill(3);
      return;
    }
  } else {
    hestats->GetXaxis()->SetBinLabel(2, "inact1");
    hestats->GetXaxis()->SetBinLabel(3, "inact2");
    hestats->GetXaxis()->SetBinLabel(4, "inact3");
  }
 
  StoreObjPtr<EventMetaData> eventMetaDataPtr;
  int run = eventMetaDataPtr->getRun();
  if (m_isRun)m_run = run;
  int nTracks = m_dedxTracks.getEntries();
  if (nTracks >= 4) {
    B2WARNING("too many tracks: unclean bhabha or radee event: " << nTracks);
    hestats->Fill(4);
    return;
  }
 
  hestats->Fill(5);
 
  //Collector object access
  auto tree = getObjectPtr<TTree>("tree");
  auto htstats = getObjectPtr<TH1I>("htstats");
  auto hmeans = getObjectPtr<TH1D>("means");
 
  for (int idedx = 0; idedx < nTracks; idedx++) {
 
    CDCDedxTrack* dedxTrack = m_dedxTracks[idedx];
    if (!dedxTrack) {
      B2WARNING("No dedx track: Going back: " << idedx);
      continue;
    }
 
    const Track* track = dedxTrack->getRelatedFrom<Track>();
    if (!track) {
      B2WARNING("No track: Going back: " << idedx);
      continue;
    }
 
    const TrackFitResult* fitResult = track->getTrackFitResultWithClosestMass(Const::pion);
    if (!fitResult) {
      B2WARNING("No related fit for this track...");
      continue;
    }
 
    m_dedx = dedxTrack->getDedxNoSat();
    m_p = dedxTrack->getMomentum();
    m_costh = dedxTrack->getCosTheta();
    m_charge = fitResult->getChargeSign();
    m_pt = fitResult->getTransverseMomentum();
    m_injTime = dedxTrack->getInjectionTime();
    m_injRing = dedxTrack->getInjectionRing();
    m_nhits = dedxTrack->size();
    m_nlhits = dedxTrack->getNLayerHits();
 
    htstats->Fill(0);
 
    if (m_cuts) {
      // apply cleanup cuts
      if (std::abs(fitResult->getD0()) >= 1.0)continue;
      if (std::abs(fitResult->getZ0()) >= 1.0) continue;
      htstats->Fill(1);
 
      //if outside CDC
      if (m_costh < TMath::Cos(150.0 * TMath::DegToRad()))continue; //-0.866
      if (m_costh > TMath::Cos(17.0 * TMath::DegToRad())) continue; //0.95
      htstats->Fill(2);
 
      if (m_nhits > m_maxHits) continue;
 
      //making some cuts based on acceptance
      if (m_costh > -0.55 && m_costh < 0.820) {
        if (dedxTrack->getNLayerHits() < 25)continue; //all CDC layer available here
      } else {
        if (m_costh <= -0.62 || m_costh >= 0.880) {
          if (dedxTrack->getNLayerHits() < 10)continue; //less layer available here
          if (m_costh > 0 && dedxTrack->getNLayerHits() < 13)continue;
        } else {
          if (dedxTrack->getNLayerHits() < 18)continue;
        }
      }
      htstats->Fill(3);
 
      const ECLCluster* eclCluster = track->getRelated<ECLCluster>();
      if (eclCluster and eclCluster->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons)) {
        double TrkEoverP = (eclCluster->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons)) / (fitResult->getMomentum().R());
        if (std::abs(TrkEoverP - 1.0) > m_setEoP)continue;
      }
      htstats->Fill(4);
    }
 
    // Make sure to remove all the data in vectors from the previous track
    if (m_isWire)m_wire.clear();
    if (m_isLayer)m_layer.clear();
    if (m_isDoca)m_doca.clear();
    if (m_isEnta)m_enta.clear();
    if (m_isDocaRS)m_docaRS.clear();
    if (m_isEntaRS)m_entaRS.clear();
    if (m_isDedxhit)m_dedxhit.clear();
    if (m_isADCcorr)m_adccorr.clear();
 
    // Simple numbers don't need to be cleared
    // make sure to use the truncated mean without the hadron saturation correction
 
    for (int i = 0; i < m_nhits; ++i) {
      // if (m_DBWireGains->getWireGain(dedxTrack->getWire(i)) == 0)continue;
      if (m_isWire)m_wire.push_back(dedxTrack->getWire(i));
      if (m_isLayer)m_layer.push_back(dedxTrack->getHitLayer(i));
      if (m_isDoca)m_doca.push_back(dedxTrack->getDoca(i));
      if (m_isEnta)m_enta.push_back(dedxTrack->getEnta(i));
      if (m_isDocaRS)m_docaRS.push_back(dedxTrack->getDocaRS(i) / dedxTrack->getCellHalfWidth(i));
      if (m_isEntaRS)m_entaRS.push_back(dedxTrack->getEntaRS(i));
      if (m_isDedxhit)m_dedxhit.push_back(dedxTrack->getDedx(i));
      if (m_isADCcorr)m_adccorr.push_back(dedxTrack->getADCCount(i));
    }
 
    // Make sure to remove all the data in vectors from the previous track
    if (m_islLayer)m_lLayer.clear();
    if (m_islDedx)m_ldedx.clear();
 
    for (int i = 0; i < m_nlhits; ++i) {
      if (m_islLayer)m_lLayer.push_back(dedxTrack->getLayer(i));
      if (m_islDedx)m_ldedx.push_back(dedxTrack->getLayerDedx(i));
    }
    // Track and/or hit information filled as per config
    htstats->Fill(5);
    hmeans->Fill(m_dedx);
    tree->Fill();
  }
}