PlotMuonClusterModule.cxx

Go to the documentation of this file.

// $Id: PlotMuonClusterModule.cxx,v 1.1 2002/06/24 16:27:36 bv Exp $
//
// PlotMuonClusterModule
//
// Begin_Html<img src="../../pedestrians.gif" align=center>
// <a href="../source_warning.html">Warning for beginners</a>.<br>
//
// This module plots clusters and fitted muon tracks to allow a
// before/after comparison.
//
// Only the Ana method has been implemented.
//
// Author:  P.S. Miyagawa 10/2000
//
// Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and
// <a href="../JOBC_Classes.html"> Job Control Classes</a> (part of
// <a href="../index.html">The MINOS Class User Guide</a>)End_Html

#include "TAxis.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TLine.h"
#include "TObjArray.h"
#include "TString.h"
#include "TStyle.h"

#include "CandData/CandRecord.h"
#include "Conventions/Munits.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "BubbleSpeak/CandDigiPairHandle.h"
#include "BubbleSpeak/CandDigiPairListHandle.h"
#include "BubbleSpeak/CandMSTClusterListHandle.h"
#include "BubbleSpeak/CandStraightClusterHandle.h"
#include "BubbleSpeak/CandThruMuonHandle.h"
#include "BubbleSpeak/CandThruMuonListHandle.h"
#include "BubbleSpeak/PlotMuonClusterModule.h"

#include <iostream.h>
#include <fstream.h>   

//......................................................................

CVSID("$Id: PlotMuonClusterModule.cxx,v 1.1 2002/06/24 16:27:36 bv Exp $");

// Declare this module to JobControll. Arguments are:
//  (1) The class name 
//  (2) The human-readable name 
//  (3) A short, human-readable description of what the module does
JOBMODULE(PlotMuonClusterModule, 
          "PlotMuonClusterModule",
          "Plots the clusters of muons from a CandThruMuonList");

//......................................................................

PlotMuonClusterModule::PlotMuonClusterModule() :
  fListClust(""), fListMuon(""),
  fAnaCanvas(0)
{
//
//  Purpose:    Default constructor.
//
//  Arguments:  n/a
//
//  Return:     n/a
//

  MSG("BubJobC", Msg::kDebug)
     << "PlotMuonClusterModule::Constructor" << endl;
  for (int i=0; i<4; i++) {
    fAnaBlanks[i] = 0;
    fAnaGraphs[i] = 0;
  }
  for (int i=0; i<2; i++) {
    fAnaFits[i] = 0;
  }
}

//......................................................................

PlotMuonClusterModule::~PlotMuonClusterModule() 
{
//
//  Purpose:    Default destructor.
//
//  Arguments:  n/a
//
//  Return:     n/a
//

  MSG("BubJobC", Msg::kDebug)
     << "PlotMuonClusterModule::Destructor" << endl;

/*
  delete fAnaCanvas;     fAnaCanvas = 0;
  for (int i=0; i<4; i++) {
cout << "loop " << i << endl;
    if (fAnaBlanks[i]) {
cout << "Deleting blank" << endl;
      delete fAnaBlanks[i];
cout << "Deleted blanks." << endl;
      fAnaBlanks[i] = 0;
    }
    if (fAnaGraphs[i]) {
      fAnaGraphs[i]->Delete();
cout << "Deleting cluster" << endl;
      delete fAnaGraphs[i];
cout << "Deleted clusters." << endl;
      fAnaGraphs[i] = 0;
    }
  }
  for (int i=0; i<2; i++) {
    if (fAnaFits[i]) {
      fAnaFits[i]->Delete();
      delete fAnaFits[i];
      fAnaFits[i] = 0;
    }
cout << "Deleted fits." << endl;
  }
*/
}

//......................................................................

JobCResult PlotMuonClusterModule::Ana(const MomNavigator *mom)
{
//
//  Purpose:  Method to plot reconstructed clusters.
//
//  Arguments:
//    mom       in    MomNavigator from which to retrieve fragment.
//
//  Return:   No decision.
//
//  Note:     This is useful for cluster parameter tuning.
//

  MSG("BubJobC", Msg::kDebug) << "PlotMuonClusterModule::Ana" << endl;

// Find PrimaryCandidateRecord fragment in MOM.
  CandRecord *candrec = dynamic_cast<CandRecord *>
             (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (candrec == 0) {
    MSG("BubJobC", Msg::kWarning) << "No PrimaryCandidateRecord in MOM."
                                                                << endl;
    return JobCResult::kError;                             // Fail event
  }

  CandMSTClusterListHandle *cllh=dynamic_cast<CandMSTClusterListHandle*>
                    (candrec->FindCandHandle("CandMSTClusterListHandle",
                                                    fListClust.Data()));
  CandThruMuonListHandle *cmlh = dynamic_cast<CandThruMuonListHandle *>
                      (candrec->FindCandHandle("CandThruMuonListHandle",
                                                     fListMuon.Data()));

  MSG("BubJobC", Msg::kDebug) << "before plot" << endl;
  PlotClusters(*cllh, *cmlh);
  MSG("BubJobC", Msg::kDebug) << "after plot" << endl;

  return JobCResult::kAOK;
}

//......................................................................

void PlotMuonClusterModule::HandleCommand(JobCommand *command) 
{
//
//  Purpose:  Method to interpret module commands.
//
//  Arguments:
//    command   in    Command to interpret.
//
//  Return:   n/a
//
//  Commands implemented:
//    Set ListClust   Set name of CandMSTClusterList to retrieve.
//    Set ListMuon    Set name of CandThruMuonList to retrieve.
//
//    Plot            Plot clusters.
//    SetPlot Zmin    Set minimum on horizontal axis.
//    SetPlot Zmax    Set maximum on horizontal axis.
//    SetPlot Umin    Set minimum on vertical U-axis.
//    SetPlot Umax    Set maximum on vertical U-axis.
//    SetPlot Vmin    Set minimum on vertical V-axis.
//    SetPlot Vmax    Set maximum on vertical V-axis.
//

   MSG("BubJobC", Msg::kDebug)
      << "PlotMuonClusterModule::HandleCommand" <<endl;

   TString cmd = command->PopCmd();

// Set reconstruction parameters.
   if (cmd == "Set") {    // Set cut threshold.
      TString opt = command->PopOpt();
      if      (opt == "ListClust") fListClust = command->PopOpt();
      else if (opt == "ListMuon")  fListMuon  = command->PopOpt();
      else {
         MSG("BubJobC", Msg::kWarning)
            << "PlotMuonClusterModule: Unrecognized option " << opt
            << endl;
      }
   }

// Plot clusters.
   else if (cmd == "Plot") RedrawPlot();

// Set cluster plot limits.
   else if (cmd == "SetPlot") {
      TString opt = command->PopOpt();
      if      (opt == "Zmin")
                      fAnaZmin    = command->PopFloatOpt() * Munits::cm;
      else if (opt == "Zmax")
                      fAnaZmax    = command->PopFloatOpt() * Munits::cm;
      else if (opt == "Umin")
                      fAnaTmin[0] = command->PopFloatOpt() * Munits::cm;
      else if (opt == "Umax")
                      fAnaTmax[0] = command->PopFloatOpt() * Munits::cm;
      else if (opt == "Vmin")
                      fAnaTmin[1] = command->PopFloatOpt() * Munits::cm;
      else if (opt == "Vmax")
                      fAnaTmax[1] = command->PopFloatOpt() * Munits::cm;
      else {
         MSG("BubJobC", Msg::kWarning)
            << "PlotMuonClusterModule: Unrecognized option " << opt
            << endl;
      }
   }

// Invalid command.
   else {
      MSG("BubJobC", Msg::kWarning)
       << "PlotMuonClusterModule: Unrecognized command " << cmd << endl;
   }
}

//......................................................................

void PlotMuonClusterModule::Help() 
{
//
//  Purpose:    Print help information for PlotMuonClusterModule.
//
//  Arguments:  n/a
//
//  Return:     n/a
//

  MSG("BubJobC", Msg::kInfo)
    << "Help for 'PlotMuonClusterModule':" << endl
    << " PlotMuonClusterModule is a module for viewing clusters" << endl
    << " and fitted muon tracks to allow a visual comparison of" << endl
    << " the results of fit rejection. The clusters are" << endl
    << " displayed by adding:" << endl
    << "     PlotMuonClusterModule::Ana" << endl
    << " to the Path to be run. Note that the colours between" << endl
    << " the two cluster graphs do not necessarily correspond," << endl
    << " but do for the muon tracks." << endl
    << endl
    << "Commands implemented:" << endl
    << " /PlotMuonClusterModule/Set ListClust <listname>" << endl
    << "     Set name of CandMSTClusterList to retrieve." << endl
    << " /PlotMuonClusterModule/Set ListMuon <listname>" << endl
    << "     Set name of CandThruMuonList to retrieve." << endl
    << endl
    << " /PlotMuonClusterModule/Plot" << endl
    << "     Draw cluster plot with current limits." << endl
    << " /PlotMuonClusterModule/SetPlot <Zmin/Zmax/Umin/Umax/Vmin/Vmax>"
    << " <value>" << endl
    << "     Set axis limit value. The command" << endl
    << "       /PlotMuonClusterModule/Plot" << endl
    << "     must be reissued afterwards to update the plot." << endl
    << endl;
}

//......................................................................

void PlotMuonClusterModule::PlotClusters(
                                   const CandMSTClusterListHandle &cllh,
                                   const CandThruMuonListHandle &cmlh)
{
//
//  Purpose:  Prepare clusters for plotting, then call method to draw
//            plot.
//
//  Arguments:
//    cllh    in     Handle to list of clusters to plot.
//    cmlh    in     Handle to list of muons with clusters to plot.
//
//  Return:   n/a
//

  // NT
 
  iEvent++;

// Clear old graphs and lines and create new ones.
  for (Int_t i=0; i<4; i++) {
    if (fAnaGraphs[i]) {
      fAnaGraphs[i]->Delete();
      delete fAnaGraphs[i];
    }
    fAnaGraphs[i] = new TObjArray();
  }
  for (Int_t i=0; i<2; i++) {
    if (fAnaFits[i]) {
      fAnaFits[i]->Delete();
      delete fAnaFits[i];
    }
    fAnaFits[i] = new TObjArray();
  }

// Initialize limits to detector extents.
  UgliGeomHandle ugh(*cllh.GetVldContext());
  ugh.GetTransverseExtent(PlaneView::kU, fAnaTmin[0], fAnaTmax[0]);
  ugh.GetTransverseExtent(PlaneView::kV, fAnaTmin[1], fAnaTmax[1]);
  ugh.GetZExtent(fAnaZmin, fAnaZmax);

/*
// Initialize limits for all graphs.
  fAnaZmin    =  1e7;
  fAnaZmax    = -1e7;
  fAnaTmin[0] =  1e7;
  fAnaTmin[1] =  1e7;
  fAnaTmax[0] = -1e7;
  fAnaTmax[1] = -1e7;
*/

// Iterate over clusters.
  Int_t npts[2] = { 0,0 };
  CandMSTClusterHandleItr clhItr(cllh.GetDaughterIterator());
  while (CandMSTClusterHandle *clh = clhItr()) {

  MSG("BubJobC", Msg::kDebug) << "before strip orientation" << endl;

// Determine strip orientation of cluster.
    Int_t clusv = -1;    // 0=U, 1=V
    switch (clh->GetPlaneView()) {
    case PlaneView::kU:
      clusv = 0;
      break;
    case PlaneView::kV:
      clusv = 1;
      break;
    case PlaneView::kA:
    case PlaneView::kB:
      MSG("BubJobC", Msg::kInfo)
         << "Plane views A and B not plotted." << endl;
      continue;  // jump to next cluster
      break;
    default:
      MSG("BubJobC", Msg::kWarning)
         << "Unrecognized strip orientation." << endl;
      continue;  // jump to next cluster
    }

  MSG("BubJobC", Msg::kDebug) << "after strip orientation" << endl;

// Fill new graph.
    MSG("BubJobC", Msg::kDebug) << "start TGraph" << endl;
    TGraph *clusg = new TGraph();
    Int_t cluspts = 0;    // Number of points in clusg.
    CandDigiPairHandleItr chhItr(clh->GetDaughterIterator());
    MSG("BubJobC", Msg::kDebug) << "after GetDaughterIterator" << endl;
//    while (CandDigiPairHandle *chh = chhItr()) {
//      MSG("BubJobC", Msg::kDebug) << "Iteration loop" << endl;
//      if(clusg&&chh){
//        clusg->SetPoint(cluspts++, chh->GetZPos(), chh->GetTPos());
//      }
//    }
    MSG("BubJobC", Msg::kDebug) << "after Iteration" << endl;
    clusg->SetMarkerStyle(8);
    clusg->SetMarkerSize(0.5);
    clusg->SetMarkerColor(npts[clusv]+2);

// Update limits.
/*
    if (clh->GetZmin() < fAnaZmin) fAnaZmin = clh->GetZmin();
    if (clh->GetZmax() > fAnaZmax) fAnaZmax = clh->GetZmax();
    if (clh->GetTmin() < fAnaTmin[clusv])
                                   fAnaTmin[clusv] = clh->GetTmin();
    if (clh->GetTmax() > fAnaTmax[clusv])
                                   fAnaTmax[clusv] = clh->GetTmax();
*/

// Store graph for future plotting and deletion.
    fAnaGraphs[2*clusv]->AddAtAndExpand(clusg, npts[clusv]++);
  }

// Iterate over muons.

  MSG("BubJobC", Msg::kDebug) << "Iterate over muons." << endl;
  Int_t nmus = 0;
  CandThruMuonHandleItr cmhItr(cmlh.GetDaughterIterator());
  while (CandThruMuonHandle *cmh = cmhItr()) {

// Retrieve U-cluster, fill new graph, and save for future plotting.
    const CandStraightClusterHandle *cchU = cmh->GetClusterU();
    TGraph *clusg = new TGraph();
    Int_t cluspts = 0;    // Number of points in clusg. 
    Int_t clusU=0;  
    Float_t chiSqU=0;  
    Float_t trueTransU=0;
    CandDigiPairHandleItr chhUItr(cchU->GetDaughterIterator());
    while (CandDigiPairHandle *chhU = chhUItr()) {
      clusg->SetPoint(cluspts++, chhU->GetZPos(), chhU->GetTPos());
// NT
      clusU++;  
      chiSqU += (chhU->GetTPos() - cmh->GetFitUFromZ(chhU->GetZPos()))*
                (chhU->GetTPos() - cmh->GetFitUFromZ(chhU->GetZPos()));
// NT    
    }

    chiSqU /= 0.0001*clusU;
    cout << "clusU, chiSqU: " << clusU << " " << chiSqU << endl;
    clusU=0;
    if(chiSqU < fMaxChiSq){
      CandDigiPairHandleItr chhUItr2(cchU->GetDaughterIterator());
      while (CandDigiPairHandle *chhU = chhUItr2()) {
        clusg->SetPoint(cluspts++, chhU->GetZPos(), chhU->GetTPos());
        clusU++;  
        ofstream out("./tracks_data.out", ios::app); 
        out << iEvent << " " << 0 << " " << clusU << " "
            << chhU->GetZPos() << " " << chhU->GetTPos() << " "
            << chhU->GetPlane() << " " << chhU->GetCharge() << endl;    
      }
    }

    clusg->SetMarkerStyle(8);
    clusg->SetMarkerSize(0.5);
    clusg->SetMarkerColor(nmus+2);
    fAnaGraphs[1]->AddAtAndExpand(clusg, nmus);

// Retrieve V-cluster, fill new graph, and save for future plotting.
    const CandStraightClusterHandle *cchV = cmh->GetClusterV();
    clusg = new TGraph();
    cluspts = 0;
    Int_t clusV=0;    
    Float_t chiSqV=0;  
    Float_t trueTransV=0;
    CandDigiPairHandleItr chhVItr(cchV->GetDaughterIterator());
    while (CandDigiPairHandle *chhV = chhVItr()) {
      clusg->SetPoint(cluspts++, chhV->GetZPos(), chhV->GetTPos()); 
// NT
      clusV++;  
      chiSqV += (chhV->GetTPos() - cmh->GetFitVFromZ(chhV->GetZPos()))*
                (chhV->GetTPos() - cmh->GetFitVFromZ(chhV->GetZPos()));      
// NT

    }    
    chiSqV /= 0.0001*clusV;
    cout << "clusV, chiSqV: " << clusV << " " << chiSqV << endl;
    clusV=0; 
    if(chiSqV < fMaxChiSq){
      CandDigiPairHandleItr chhVItr2(cchV->GetDaughterIterator());
      while (CandDigiPairHandle *chhV = chhVItr2()) {
        clusg->SetPoint(cluspts++, chhV->GetZPos(), chhV->GetTPos());
        clusV++;  
        ofstream out("./tracks_data.out", ios::app); 
        out << iEvent << " " << 1 << " " << clusV << " "
            << chhV->GetZPos() << " " << chhV->GetTPos() << " "
            << chhV->GetPlane() << " " <<  chhV->GetCharge()<< endl;    
      }
    }

    clusg->SetMarkerStyle(8);
    clusg->SetMarkerSize(0.5);
    clusg->SetMarkerColor(nmus+2);
    fAnaGraphs[3]->AddAtAndExpand(clusg, nmus);

// Set fit line endpoints.
    Float_t zbeg = cmh->GetZBeg();
    Float_t zend = cmh->GetZEnd();
    Float_t ubeg = cmh->GetFitUFromZ(zbeg);
    Float_t uend = cmh->GetFitUFromZ(zend);
    Float_t vbeg = cmh->GetFitVFromZ(zbeg);
    Float_t vend = cmh->GetFitVFromZ(zend);

// Create fit lines for each view.
    TLine *fitl = new TLine(zbeg, ubeg, zend, uend);
    fitl->SetLineColor(nmus+2);
    fAnaFits[0]->AddAtAndExpand(fitl, nmus);
    fitl = new TLine(zbeg, vbeg, zend, vend);
    fitl->SetLineColor(nmus+2);
    fAnaFits[1]->AddAtAndExpand(fitl, nmus++);

// Update plot limits.
    if (zbeg < fAnaZmin) fAnaZmin = zbeg;
    if (zend > fAnaZmax) fAnaZmax = zend;
    if (ubeg < fAnaTmin[0]) fAnaTmin[0] = ubeg;
    if (ubeg > fAnaTmax[0]) fAnaTmax[0] = ubeg;
    if (uend < fAnaTmin[0]) fAnaTmin[0] = uend;
    if (uend > fAnaTmax[0]) fAnaTmax[0] = uend;
    if (vbeg < fAnaTmin[1]) fAnaTmin[1] = vbeg;
    if (vbeg > fAnaTmax[1]) fAnaTmax[1] = vbeg;
    if (vend < fAnaTmin[1]) fAnaTmin[1] = vend;
    if (vend > fAnaTmax[1]) fAnaTmax[1] = vend;
  }

// Draw plot.
  MSG("BubJobC", Msg::kDebug) << "Draw plot" << endl;
  RedrawPlot();
}

//......................................................................

void PlotMuonClusterModule::RedrawPlot()
{
//
//  Purpose:    Draw plot of clusters.
//
//  Arguments:  n/a
//
//  Return:     n/a
//

// Turn off statistics.
  gStyle->SetOptStat(0);

// Start with clear canvas.
  if (fAnaCanvas) fAnaCanvas->Clear("D");
  else {
    fAnaCanvas = new TCanvas("PlotMuonCanvas", "Muon Cluster Canvas",
                             1200, 700);
    fAnaCanvas->SetFillColor(0);
    fAnaCanvas->Divide(2, 2);
  }

// Plot clusters.
  for (Int_t i=0; i<4; i++) {
    fAnaCanvas->cd(i+1);
    gPad->SetGrid();

// Delete old blank and create new one.
    delete fAnaBlanks[i];
    char blankname[80];
    sprintf(blankname, "blank%d", i+1);
    fAnaBlanks[i] = new TH2F(blankname, "Clusters",
                2, fAnaZmin, fAnaZmax, 2, fAnaTmin[i/2], fAnaTmax[i/2]);

    TString xtitle("z-value (");
    xtitle += Munits::base_length_name;
    xtitle += ")";
    // const_cast is only necessary until ROOT 3.00 is universal
    fAnaBlanks[i]->SetXTitle(const_cast<char*>(xtitle.Data()));

    TString ytitle((i/2 ? "v-value (" : "u-value ("));
    ytitle += Munits::base_length_name;
    ytitle += ")";
    // const_cast is only necessary until ROOT 3.00 is universal
    fAnaBlanks[i]->SetYTitle(const_cast<char*>(ytitle.Data()));

    fAnaBlanks[i]->Draw();

// Plot graphs.
    for (Int_t j=0; j<fAnaGraphs[i]->GetEntries(); j++) {
      if (TGraph *grj = (TGraph*) fAnaGraphs[i]->At(j)) grj->Draw("P");
    }
  }

// Plot fit lines.
  for (Int_t i=0; i<2; i++) {
    fAnaCanvas->cd(2*i+2);
    for (Int_t j=0; j<fAnaFits[i]->GetEntries(); j++) {
      if (TLine *fitl = dynamic_cast<TLine *>(fAnaFits[i]->At(j)))
         fitl->Draw();
    }
  }

// Update.
  fAnaCanvas->Update();
}

---