SWGuideAlignmentMonitors < CMSPublic

---+!! Creating a Monitoring Module for =AlignmentProducer=
%COMPLETE5%

%TOC{title="Contents:"}%

---++ Goal of this page

See [[https://twiki.cern.ch/twiki/bin/view/CMS/SWGuideAlignmentAlgorithms#Monitoring][Alignment algorithms::Monitoring]] for an introduction to the =CommonAlignmentMonitor= package.  This page explains how to make your own histogram module.

---++ A consistent release

You need to start with a consistent set of tags to compile and run the =AlignmentProducer=, because (a) you'll need it to test your module and (b) what else would you be monitoring?  _(Need to update if this package is ever applied to calibration...)_

*Update:* version 1_5_0 of !AlignmentProducer works with the latest version of !CommonAlignmentMonitor (V00-02-02).  Use these or a later version.  If you need to see the recipe for 1_3_X and 1_4_X, look in the history of this twiki page (r5).

---++ Your new monitor

You will make a new subclass of =AlignmentMonitorBase=, following =AlignmentMoniterHIP= as a model.  The name should conform to =AlignmentMonitorXXX= to keep all histogram modules in the same corner of =SEAL= namespace.  Change directories to =Alignment/CommonAlignmentMonitor/plugins/= and create a new file named =AlignmentMonitorXXX.cc=.  (If using the old plugin manager, put the file in =/src/=.)  Fill the new file with
<verbatim>
// -*- C++ -*-
//
// Package:     CommonAlignmentProducer
// Class  :     AlignmentMonitorXXX
// 
// Implementation:
//     <Notes on implementation>
//
// Original Author:  Who?
//         Created:  When?
// $Id: SWGuideAlignmentMonitors.txt,v 1.6 2007/07/09 13:52:15 pivarski Exp $
//

#include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
#include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"

class AlignmentMonitorXXX: public AlignmentMonitorBase {
   public:
      AlignmentMonitorXXX(const edm::ParameterSet& cfg): AlignmentMonitorBase(cfg) { };
      ~AlignmentMonitorXXX() {};

      void book();
      void event(const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& iTrajTracks);
      void afterAlignment(const edm::EventSetup &iSetup);

   private:
};

void AlignmentMonitorHIP::book() {
}

void AlignmentMonitorHIP::event(const edm::EventSetup &iSetup, const ConstTrajTrackPairCollection& tracks) {
}

void AlignmentMonitorHIP::afterAlignment(const edm::EventSetup &iSetup) {
}

// DEFINE_SEAL_MODULE();
// DEFINE_SEAL_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorXXX, "AlignmentMonitorXXX");
DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorXXX, "AlignmentMonitorXXX");
</verbatim>

If using the old plugin manager, you need to use the =DEFINE_SEAL_MODULE()=, =DEFINE_SEAL_PLUGIN()= sequence; if using the new plugin manager (=1_5_0_pre1= or later), use =DEFINE_EDM_PLUGIN()=.  Also, if using the new plugin manager, add
<verbatim>
<library name="AlignmentMonitorXXX" file="AlignmentMonitorXXX.cc">
 <use name=Alignment/CommonAlignmentMonitor>
 <use name=TrackingTools/TrackFitters>
 <flags EDM_PLUGIN=1>
</library>
</verbatim>
to =Alignment/CommonAlignmentMonitors/plugins/BuildFile=.

_One last thing._  (Sorry about this mess.)  If you're using the old plugin manager, you can't compile multiple monitors because the =DEFINE_SEAL_MODULE()= line should appear only once.  The proper way to handle this would have been to make a =SealModule.cc= file, and split all the monitors into =.cc= and =.h= files.  To get a local copy working, just delete =AlignmentMonitorHIP.cc=.  We'll have to adapt everything to the new plugin manager to put it in the database anyway.

Now it should compile.

I've found that sometimes new monitors are not automatically loaded into =SEAL=.  To force it to be loaded, run
<verbatim>
SealPluginsDump | grep AlignmentMonitor
</verbatim>
This command is available after you run =eval `scramv1 run -sh`= or =-csh=.

---+++ Adding histograms

*1. Declare your histogram* in the =private:= field of your class definition, e.g.
<verbatim>
   private:
      TH1F *m_hist;
</verbatim>
The histogram may be anything that descends from =TH1= (that is, any ROOT histogram, including 2D histograms and profile plots).

*2. Book it* in the =book()= method like this:
<verbatim>
void AlignmentMonitorHIP::book() {
    m_hist = (TH1F*)(add("/", new TH1F("hist", "normal constructor", 100, 0., 1.)));
}
</verbatim>

The =add()= function manages the ROOT file, handling iterations and job collection automatically.  Let's step through the program flow:
   * First, the ROOT constructor (in this case, =TH1F=) does its thing in the normal way.
   * This is passed to =add()= with a string ="/"=.  The string tells =add()= where in the ROOT file the histogram belongs.
   * If this is a new ROOT file (first iteration), =add()= will create whatever parent directories are needed and place the histogram at that point in the file.
   * If this is an old ROOT file (second or later iteration), =add()= will find the old histogram in the file and throw away the new one you just made.
   * The =add()= function returns a histogram, either the new one or the old one.  (Unfortunately, you need to cast it.  But that isn't _so_ bad, because the cast isn't very far away from the constructor, so it's easy to check for incorrect casts.)
On the first iteration, =book()= books the histograms and saves them in the normal way, but on subsequent iterations, =book()= is actually resetting your =m_hist= pointers to the histograms they belong to: the newly constructed histogram becomes a look-up specifier that is discarded when the search is complete.  This is convenient when an iterative algorithm is executed in several invocations of =cmsRun=, and we want histograms that use the iteration number as data (e.g. residual RMS versus iteration number).

Most histograms, however, should be reset with each iteration.  Better still, they should be booked anew so that iteration N results can be compared with iteration N+1.  To do this, replace the ="/"= string with ="/iterN/"= (both slashes are necessary).  This is code for ="/iter1/"= on the first iteration, ="/iter2/"= on the second iteration, etc.  The histogram specified by the constructor will not be thrown away; it will be placed into the directory appropriate to the current iteration.  You can make arbitrarily deep directory specifications, such as "/iterN/residprofiles/" and "/iterN/overlapplots/station11/": the whole subdirectory structure will be copied from iteration to iteration.

Directory names always begin and end with slashes.  Do not use the standard ROOT =TFile= and =TDirectory= tools: make all directories through the =add()= function.  Don't make a directory named "/iter2/" or something--- that's just asking for trouble.

The =add()= function also merges histograms in a collection job after parallel processing.  In this case, newly-constructed histograms are saved into the grand total ROOT file, and =add()= uses the names and directories to find all the subjob histograms that need to be merged into the grand total.  All of that happens automatically.

You may be wondering, I see a "new", where's the "delete"?  It happens at the end of the iteration in !AlignmentMonitorBase, but only for histograms that _should_ be deleted at the end of each iteration (histograms in the "iterN" directories).  That is to say, your plugin creates references to histograms, and the !AlignmentMonitor framework "steals" those references; you are no longer responsible for deleting them.  So don't delete them, or you'll cause a segmentation fault!

*3. Fill your histograms.*  You may do this in =event()= or =afterAlignment()=, whichever is appropriate.

The =event(const edm::EventSetup &iSetup, const= =ConstTrajTrackPairCollection= =&tracks)= function is called in the event loop, and it supplies a list of trajectory-track pairs.  Here's how to iterate over them:
<verbatim>
for (ConstTrajTrackPairCollection::const_iterator iter = tracks.begin();  iter != tracks.end();  ++iter) {
    const Trajectory *traj = iter->first;
    const reco::Track *track = iter->second;

    std::vector<TrajectoryMeasurement> trajectoryIntersections = traj->measurements();
    for (std::vector<TrajectoryMeasurement>::const_iterator interIter = trajectoryIntersections.begin();  interIter != trajectoryIntersections.end();  ++interIter) {
        const TrajectoryMeasurement trajectoryIntersection = *interIter;
        const TransientTrackingRecHit *hit = &(*interIter.recHit());
        ...
    }
}
</verbatim>

The =afterAlignment(const edm::EventSetup &iSetup)= function is called at the end of an iteration, after updating the alignable geometry (=AlignableTracker= and =AlignableMuon=) but before updating the production geometry (=TrackerGeometry=, =DTGeometry=, and =CSCGeometry=).  You can walk through the =AlignableTracker=, =AlignableMuon=, or =AlignableParameterStore= hierarchies to get the new orientations.

You have access to the current =AlignableTracker= through =pTracker()=, the =AlignableMuon= through =pMuon()=, the =AlignmentParameterStore= through =pStore()=, and the =AlignableNavigator= through =pNavigator()= (all are pointers, NULL if not defined).  You can also access the current iteration number through =iteration()= (it starts with 1).

---+++ Ntuples or =TTrees=

Before creating thousands of static histograms, it's good to explore the data with an ntuple, at least to get the binning right.  The =add()= function described above places any =TObject= in the file, including =TTrees=, so you can do the following in =book()=:
<verbatim>
m_tree = (TTree*)(add("/iterN/", new TTree("tree", "tree")));
</verbatim>
_However,_ we still need to add branches, and =AlignmentMonitorBase= has no equivalent of =add()= for branches.  This means that if you simply type
<verbatim>
m_tree->Branch("x", &m_x, "x/F");
m_tree->Branch("y", &m_y, "y/F");
m_tree->Branch("z", &m_z, "z/F");
</verbatim>
then new branches will be added and re-added with every iteration.  You can get around this by putting the above in an 
<verbatim>
if (iteration() == 1) { }
</verbatim>
block, but then on subsequent iterations, you'll have no way to access ="x/F"=, ="y/F"=, and ="z/F"=.  Therefore, if your ntuple is not in an ="/iterN/"= directory, be sure that you only =Fill()= it in iteration 1 (or you will fill it with zeros, or NaNs, or random addresses in memory).  If the ntuple is in ="/iterN/"=, adding branches will work, since it's a new ntuple with each iteration.

*Here are some hard-and-fast rules:* if you want to create a new ntuple with each iteration, do so like this:
<verbatim>
m_tree = (TTree*)(add("/iterN/", new TTree("tree", "tree"))
m_tree->Branch("x", &m_x, "x/F")
</verbatim>
If you want to create one ntuple, for use in all iterations, do so like this:
<verbatim>
m_tree = ((TTree*)(add("/", new TTree("tree", "tree"))
if (iteration() == 1) {
    m_tree->Branch("x", &m_x, "x/F")
}
</verbatim>
and then you _must_ use
<verbatim>
replace AlignmentProducer.maxLoops = 22
</verbatim>
to iterate, you cannot call =cmsRun= multiple times.

Ntuples are not merged in a collection job.  It's unclear whether we'd want to copy all the ntuple data or merely link it up with a =TChain=.  Anyway, we'll only be parallel processing when we understand our algorithms well enough to have static histograms, anyway.

---+++ Configuring =AlignmentProducer= to use your new monitoring package

That was described in [[https://twiki.cern.ch/twiki/bin/view/CMS/SWGuideAlignmentAlgorithms#Monitoring][Alignment algorithms::Monitoring]].

---++ Recipes

Histogram code will involve a lot of similar constructions, such as calculating residuals.  Here are some code snippets that you can use to quickly make your plots.  I encourage others to add to this list (that's what wikis are for)!

---+++ Calculating residuals

<verbatim>
#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
</verbatim>
and in your loop over hits,

<verbatim>
const DetId id = hit->geographicalId();
if (hit->isValid()  &&  pNavigator()->detAndSubdetInMap(id)) {
    TrajectoryStateOnSurface combinedTrajInter = tsoscomb.combine(trajectoryIntersection.forwardPredictedState(), trajectoryIntersection.backwardPredictedState());
    double residual = combinedTrajInter.localPosition().x() - hit->localPosition().x());
}
</verbatim>

---+++ Booking histograms from selected =Alignables=

In the class declaration,
<verbatim>
      std::map<Alignable*, TH1F*> m_residuals;
</verbatim>

and in =book()=,
<verbatim>
   std::vector<Alignable*> alignables = pStore()->alignables();
   for (std::vector<Alignable*>::const_iterator it = alignables.begin();  it != alignables.end();  ++it) {
      char name[256], title[256];
      sprintf(name,  "xresid%d", (*it)->geomDetId().rawId());
      sprintf(title, "x residual for DetId %d (cm)", (*it)->geomDetId().rawId());

      m_residuals[*it] = (TH1F*)(add("/iterN/", new TH1F(name, title, 100, -5., 5.)));
   }
</verbatim>

and then in =event()=,
<verbatim>
const DetId id = hit->geographicalId();
Alignable *alignable = pNavigator()->alignableFromDetId(id);
std::map<Alignable*, TH1F*>::const_iterator search = m_residuals.find(alignable);
while (search == m_residuals.end()  &&  (alignable = alignable->mother())) search = m_residuals.find(alignable);
if (search != m_residuals.end()) {
    search->second->FIll(residual);
}
</verbatim>

#ReviewStatus
---++!! Review Status

| *Reviewer/Editor and Date* | *Comments* |
| Main.pivarski - 07 May 2007 | Page created |

%RESPONSIBLE% Main.pivarski %BR%
%REVIEW% Reviewer
This topic: CMSPublic > DefaultWeb > WebHome > SWGuide > SWGuideCalAli > SWGuideTrackAlignment > SWGuideAlignmentAlgorithms > SWGuideAlignmentMonitors
Topic revision: r7 - 2007-07-16 - JennyWilliams
Copyright &© 2008-2022 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
or Ideas, requests, problems regarding TWiki? use Discourse or Send feedback