AltModuleSliceQualityCheck.cxx

Go to the documentation of this file.

#include <cassert>
#include <TFile.h>
#include <TCanvas.h>

#include "AltModuleSliceQualityCheck.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandSliceList.h"
#include "RecoBase/CandSliceHandle.h"
#include "DataUtil/TruthHelper.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "RecoBase/CandStripHandle.h"

ClassImp(AltModuleSliceQualityCheck)

//____________________________________________________________________________
CVSID("$Id: AltModuleSliceQualityCheck.cxx,v 1.4 2004/07/26 12:22:52 west Exp $");
//____________________________________________________________________________
JOBMODULE(AltModuleSliceQualityCheck, "AltModuleSliceQualityCheck", "slc QC");
//____________________________________________________________________________
AltModuleSliceQualityCheck::AltModuleSliceQualityCheck() :
fSlicePurity               ("fSlicePurity", "Slice Purity",
                                                                50, 0.0, 1.0),
fSliceCompleteness         ("fSliceCompleteness", "Slice Completeness",
                                                                50, 0.0, 1.0),
fSlicePurityVsCompleteness ("fSlicePurityVsCompleteness", "Pur. Vs. Compl.",
                                                  50, 0.0, 1.0, 50, 0.0, 1.0),
fSlicePurityVsEnu          ("fSlicePurityVsEnu", "Purity vs Enu",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSliceCompletenessVsEnu    ("fSliceCompletenessVsEnu", "Completeness vs Enu",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSlicePurityVsEvis         ("fSlicePurityVsEvis", "Purity vs Evis",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSliceCompletenessVsEvis   ("fSliceCompletenessVsEvis", "Completeness vs Evis",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSlicePurityVsVtxX         ("fSlicePurityVsVtxX", "Purity vs x-vtx",
                                                50, 0.0, 1.0, 200, -5.0, 5.0),
fSliceCompletenessVsVtxX   ("fSliceCompletenessVsVtxX", "Completeness vs x-vtx",
                                                50, 0.0, 1.0, 200, -5.0, 5.0),
fSlicePurityVsVtxY         ("fSlicePurityVsVtxY", "Purity vs y-vtx",
                                                50, 0.0, 1.0, 200, -5.0, 5.0),
fSliceCompletenessVsVtxY   ("fSliceCompletenessVsVtxY", "Completeness vs y-vtx",
                                                50, 0.0, 1.0, 200, -5.0, 5.0),
fSlicePurityVsVtxZ         ("fSlicePurityVsVtxZ", "Purity vs z-vtx",
                                                50, 0.0, 1.0, 200, 0.0, 10.0),
fSliceCompletenessVsVtxZ   ("fSliceCompletenessVsVtxZ", "Completeness vs z-vtx",
                                                50, 0.0, 1.0, 200, 0.0, 10.0),
fSlicePurityVsKinX         ("fSlicePurityVsKinX", "Purity vs x",
                                                  50, 0.0, 1.0, 50, 0.0, 1.0),
fSliceCompletenessVsKinX   ("fSliceCompletenessVsKinX", "Completeness vs x",
                                                  50, 0.0, 1.0, 50, 0.0, 1.0),
fSlicePurityVsKinY         ("fSlicePurityVsKinY", "Purity vs y",
                                                  50, 0.0, 1.0, 50, 0.0, 1.0),
fSliceCompletenessVsKinY   ("fSliceCompletenessVsKinY", "Completeness vs y",
                                                  50, 0.0, 1.0, 50, 0.0, 1.0),
fSlicePurityVsKinQ2        ("fSlicePurityVsKinQ2", "Purity vs Q2",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSliceCompletenessVsKinQ2  ("fSliceCompletenessVsKinQ2","Completeness vs Q2",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSlicePurityVsKinW         ("fSlicePurityVsKinW", "Purity vs Q2",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fSliceCompletenessVsKinW   ("fSliceCompletenessVsKinW", "Completeness vs Q2",
                                                 50, 0.0, 1.0, 30, 0.0, 30.0),
fWriteHistos       (true),
fDrawHistos        (false),
fRootFileName      ("slice_quality.root")
{
  MSG("AltModule", Msg::kVerbose) << "AltModuleSliceQualityCheck::Ctor\n";
}
//____________________________________________________________________________
AltModuleSliceQualityCheck::~AltModuleSliceQualityCheck()
{
  MSG("AltModule", Msg::kVerbose) << "AltModuleSliceQualityCheck::Dtor\n";
}
//____________________________________________________________________________
void AltModuleSliceQualityCheck::BeginJob()
{
  MSG("AltModule", Msg::kVerbose) << "AltModuleSliceQualityCheck::BeginJob\n";
}
//____________________________________________________________________________
const Registry & AltModuleSliceQualityCheck::DefaultConfig() const
{
  MSG("AltModule", Msg::kDebug)
                       << "AltModuleSliceQualityCheck::DefaultConfig" << endl;
  static Registry r;
  
  r.UnLockValues();
  r.Set("WriteHistos",  true);
  r.Set("DrawHistos",   false);
  r.Set("RootFileName", "slice_quality.root");
  r.LockValues();

  return r;
}
//____________________________________________________________________________
void AltModuleSliceQualityCheck::Config(const Registry & r)
{
  int          tmpi = 0;
  const char * tmpc = 0;
  //double       tmpd = 0;

  if (r.Get("fWriteHistos",  tmpi))  fWriteHistos  = tmpi;
  if (r.Get("fDrawHistos",   tmpi))  fDrawHistos   = tmpi;
  if (r.Get("fRootFileName", tmpc))  fRootFileName = string(tmpc);
}
//____________________________________________________________________________
JobCResult AltModuleSliceQualityCheck::Ana(const MomNavigator * mom)
{
  MSG("AltModule", Msg::kVerbose) << "AltModuleSliceQualityCheck::Ana\n";

  JobCResult result(JobCResult::kPassed);

  //-- Find PrimaryCandidateRecord fragment in MOM & check it is not null.

  CandRecord *cdrec = dynamic_cast<CandRecord *>
                    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (cdrec == 0) {
     MSG("AltModule", Msg::kWarning)
                              << "NULL PrimaryCandidateRecord pointer" << endl;
     result.SetWarning().SetFailed();
     return result;
  }

  //-- Get a handle to CandSliceList (previous step reconstruction objects).

  MSG("AltModule", Msg::kDebug)
          << "Asking PrimaryCandidateRecord for CaldSliceSRListHandle" << endl;

  CandSliceListHandle * cslh = dynamic_cast<CandSliceListHandle *>
                              (cdrec->FindCandHandle("CandSliceSRListHandle"));

  //-- Require CandSliceList to exist & number of CandSlices to be non-zero

  if (!cslh) {
     MSG("AltModule", Msg::kDebug)  << "NULL CandSliceListHandle ptr" << endl;
     result.SetFailed();
     return result;

  } else {
      if(cslh->GetNDaughters() < 1) {
          MSG("AltModule", Msg::kDebug) << "Empty CandSlice list." << endl;
          result.SetFailed();
          return result;
      }
  }

  //-- Get a Truth Helper -- obviously only use this module with Monte Carlo
  
  TruthHelper truth_helper(mom);
  
  //-- Get an iterator over slices, loop over them and fill in histograms
  
  CandSliceHandleItr slice_iter ( cslh->GetDaughterIterator() );

  int islice = 0;
  
  while( CandSliceHandle * slice = slice_iter.Ptr() ) {

     if( !IsRubbishSlice(slice) ) {

        int   neu_id       = truth_helper.GetBestSliceNeuMatch(*slice);

        float purity       = truth_helper.SlicePurity(*slice);
        float completeness = truth_helper.SliceCompleteness(*slice);
        float E            = (truth_helper.GetP4Neu(neu_id))[3];
        float Evis         = 0;
        float vtx_x        = 0;
        float vtx_y        = 0;
        float vtx_z        = 0;   // get these
        float x            = 0;
        float y            = 0;
        float Q2           = 0;
        float W            = 0;
          
        MSG("AltModule", Msg::kVerbose)
              << "- slice[" << islice++ << "] purity = " << purity
                            << ", completeness = " << completeness << endl;

        fSlicePurity.              Fill ( purity       );
        fSliceCompleteness.        Fill ( completeness );
        fSlicePurityVsCompleteness.Fill ( purity, completeness );
        fSlicePurityVsEnu.         Fill ( purity       , E     );
        fSliceCompletenessVsEnu.   Fill ( completeness , E     );
        fSlicePurityVsEvis.        Fill ( purity       , Evis  );
        fSliceCompletenessVsEvis.  Fill ( completeness , Evis  );
        fSlicePurityVsVtxX.        Fill ( purity       , vtx_x );
        fSliceCompletenessVsVtxX.  Fill ( completeness , vtx_x );
        fSlicePurityVsVtxY.        Fill ( purity       , vtx_y );
        fSliceCompletenessVsVtxY.  Fill ( completeness , vtx_y );
        fSlicePurityVsVtxZ.        Fill ( purity       , vtx_z );
        fSliceCompletenessVsVtxZ.  Fill ( completeness , vtx_z );
        fSlicePurityVsKinX.        Fill ( purity       , x     );
        fSliceCompletenessVsKinX.  Fill ( completeness , x     );
        fSlicePurityVsKinY.        Fill ( purity       , y     );
        fSliceCompletenessVsKinY.  Fill ( completeness , y     );
        fSlicePurityVsKinQ2.       Fill ( purity       , Q2    );
        fSliceCompletenessVsKinQ2. Fill ( completeness , Q2    );
        fSlicePurityVsKinW.        Fill ( purity       , W     );
        fSliceCompletenessVsKinW.  Fill ( completeness , W     );
     }
     slice_iter.Next();
  }

  if(fWriteHistos) WriteHistos();
  if(fDrawHistos)  DrawHistos ();
  
  return result;
}
//____________________________________________________________________________
void AltModuleSliceQualityCheck::HandleCommand(JobCommand * cmd)
{
  // Process configuration commands
  MSG("AltModule", Msg::kDebug)
                       << "AltModuleSliceQualityCheck::HandleCommand" << endl;

  if (cmd->HaveCmd() ) {
  MSG("AltModule", Msg::kWarning)
                          << "Use Config(const Registry & r) instead" << endl;
  }
}
//____________________________________________________________________________
void AltModuleSliceQualityCheck::WriteHistos(void)
{
  TFile root_file(fRootFileName.c_str(), "RECREATE");
  
  fSlicePurity.              Write ("slice_purity");
  fSliceCompleteness.        Write ("slice_completeness");
  fSlicePurityVsCompleteness.Write ("slice_purity_vs_completeness");
  fSlicePurityVsEnu.         Write ("slice_purity_Enu");
  fSliceCompletenessVsEnu.   Write ("slice_completeness_Enu");
  fSlicePurityVsEvis.        Write ("slice_purity_Evis");
  fSliceCompletenessVsEvis.  Write ("slice_completeness_Evis");
  fSlicePurityVsVtxX.        Write ("slice_purity_vtxx");
  fSliceCompletenessVsVtxX.  Write ("slice_completeness_vtxx");
  fSlicePurityVsVtxY.        Write ("slice_purity_vtxy");
  fSliceCompletenessVsVtxY.  Write ("slice_completeness_vtxy");
  fSlicePurityVsVtxZ.        Write ("slice_purity_vtxz");
  fSliceCompletenessVsVtxZ.  Write ("slice_completeness_vtxz");
  fSlicePurityVsKinX.        Write ("slice_purity_x");
  fSliceCompletenessVsKinX.  Write ("slice_completeness_x");
  fSlicePurityVsKinY.        Write ("slice_purity_y");
  fSliceCompletenessVsKinY.  Write ("slice_completeness_y");
  fSlicePurityVsKinQ2.       Write ("slice_purity_Q2");
  fSliceCompletenessVsKinQ2. Write ("slice_completeness_Q2");
  fSlicePurityVsKinW.        Write ("slice_purity_W");
  fSliceCompletenessVsKinW.  Write ("slice_completeness_W");
  
  root_file.Close();
}
//____________________________________________________________________________
void AltModuleSliceQualityCheck::DrawHistos(void) 
{
  TCanvas * canvas = new TCanvas("canvas","purity/completeness",20,20,600,600);

  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);

  canvas->Divide(1,2);

  canvas->cd(1);
  fSlicePurity.Draw();
  fSlicePurity.GetXaxis()->SetTitle("slice purity");
  
  canvas->cd(2);
  fSliceCompleteness.Draw();
  fSlicePurity.GetXaxis()->SetTitle("slice completeness");

  canvas->Update();
  
  delete canvas;
}
//____________________________________________________________________________
bool AltModuleSliceQualityCheck::IsRubbishSlice(CandSliceHandle * slice)
{
// Returns true if all hit strips in this slice are in the muon spectrometer.
// This can only happen for the 'rubbish slice' I am using to pass unmatched
// muon spectrometer hits to subsequent reconstructon stages.
// The user can switch off this 'feature' and make proper slices out of these
// hits... but since they are of no significance and they are not reconstructed,
// anyway, at later reco. stages, the default action is to pack them into the
// 'rubbish slice'. For the quality check purposes this slice is excluded.

  bool is_rubbish = true;
  
  CandStripHandleItr strip_iter ( slice->GetDaughterIterator() );

  while( CandStripHandle * strip = strip_iter.Ptr() ) {

    if(strip->GetPlane() < 120) return false;
    
    strip_iter.Next();
  }

  return is_rubbish;
}
//____________________________________________________________________________

---