file ExCompMuon.cc
====================

// -----------------------------------------------------------------------------
//  Example program trying to match muon barrel tracks and generator level muons
// -----------------------------------------------------------------------------
//  $Date: 2001/03/22 13:21:56 $
//  $Revision: 1.13 $
// -----------------------------------------------------------------------------
//  Author: Stephan Wynhoff - RWTH-Aachen
// -----------------------------------------------------------------------------

#include "Utilities/Configuration/interface/Architecture.h"

#include "Utilities/Notification/interface/PackageInitializer.h"

#include "Examples/CompGenRec/interface/ExCompMuonInfo.h"

static PKBuilder<ExCompMuonInfo> mybuilder("ExCompMuonInfoBuilder");


file  interface/ExCompMuonInfo.h
================================


#ifndef EXCOMPMUONINFO_H
#define EXCOMPMUONINFO_H
/** \file Examples/CompGenRec/interface/ExCompMuonInfo.h
 *  Compares generated muons with reconstructed tracker tracks, 
 *  calorimetric clusters and muon barreltracks. */

#include "Utilities/Notification/interface/Observer.h" 

#include "CLHEP/Vector/LorentzVector.h"

class ExCompMuonHistos;
class G3EventProxy;
class TrackReconstructor;

/** Compares generated muons with reconstructed tracker tracks, 
 *  calorimetric clusters and muon barreltracks.
 *  Histograms are filled using class ExCompMuonHistos .
 * \author Stephan Wynhoff, CERN  */
class ExCompMuonInfo : private ActiveObserver<G3EventProxy *> {
public:
  ExCompMuonInfo();  //!< The constructor
  ~ExCompMuonInfo(); //!< The destructor
  /// this method will do the user analysis
  void analysis(G3EventProxy * ev);
private:
  /// don't change the name "upDate" - this method is mandatory. 
  void upDate(G3EventProxy * ev) { if (ev!=0) analysis(ev);}
  void getGeneratorParticles(G3EventProxy * ev);
  void getCalorimeterClusters(G3EventProxy * ev);
  void getTrackerTracks(G3EventProxy * ev);
  void getMuonTracks(G3EventProxy * ev);
private:
  ExCompMuonHistos * UserHists; //!< to access the histograms
  TrackReconstructor *  myTrackFinder;
  int ncluster, ntracks, nmuons;
  int matchmuon,matchcluster,matchtrack;
  float EMuRes, ETkRes;
  vector <Hep3Vector> GoodGenMuon, GoodRecMuon, GoodCluster, GoodTracksAtVertex, GoodTracksAtEcal;
};

#endif



file ExCompMuonInfo.cc
=======================


#include "Utilities/Configuration/interface/Architecture.h"

#include "CARF/G3Event/interface/G3EventProxy.h"
#include "CARF/Reco/interface/RecItr.h"

#include "GeneratorInterface/HepEvent/interface/RawHepEvent.h"
#include "GeneratorInterface/HepEvent/interface/HepEventG3EventProxyReader.h"
#include "GeneratorInterface/Particle/interface/RawParticle.h"
#include "GeneratorInterface/Particle/interface/RawHepEventParticle.h"

#include "Calorimetry/CaloCluster/interface/CaloCluster.h"

#include "CommonDet/PatternTools/interface/TrackReconstructor.h"
#include "CommonDet/PatternTools/interface/TTrack.h"
#include "CommonDet/PatternPrimitives/interface/GeaneWrapper.h"

#include "Tracker/GtfPattern/interface/CombinatorialTrackFinder.h"

#include "Examples/CompGenRec/interface/ExCompMuonHistos.h"
#include "Examples/CompGenRec/interface/ExCompMuonInfo.h"

#include <iostream>
#include <iomanip>

ExCompMuonInfo::ExCompMuonInfo() {
  cout << "Construct ExCompMuonInfo" << endl;
  UserHists = new ExCompMuonHistos();    // instantiate the user histogram 
  myTrackFinder=0;
}

ExCompMuonInfo::~ExCompMuonInfo() {
  cout << "Destruct ExCompMuonInfo" << endl;
  delete UserHists;
}

void ExCompMuonInfo::analysis(G3EventProxy * ev) { 

  cout << "=== Event #"<< ev->simTrigger()->id().eventInRun() 
       << "==================================="<< endl;

  ncluster=0; ntracks=0; nmuons=0;
  this->getGeneratorParticles(ev);
  this->getCalorimeterClusters(ev);
  this->getMuonTracks(ev);
  this->getTrackerTracks(ev);

  cout << "Start matching -------------------------------------" << endl;
  float delta, tmpdelta;
  int ig=0, ir=0, nmatches=0;
  unsigned int i1,i2;
  EMuRes=-99.;
  ETkRes=-9.9;
  cout << "Gen/Rec/Clu/Tk: " << GoodGenMuon.size() 
       <<"/"<<GoodRecMuon.size() 
       <<"/"<<GoodCluster.size() 
       <<"/"<<GoodTracksAtVertex.size() 
       << endl;
  matchmuon=matchcluster=matchtrack=0;
  for(i1=0; i1<GoodGenMuon.size(); i1++) {
    delta=999. ;
    cout << "GM=" << setw(2) << i1 
	 << " theta=" << setw(5) << GoodGenMuon[i1].theta()
	 << " phi=" << setw(5) << GoodGenMuon[i1].phi() << endl;;
    for (i2=0; i2<GoodRecMuon.size(); i2++) {
      cout << "RM=" << setw(2) << i2 
	   << " theta=" << setw(5) << GoodRecMuon[i2].theta()
	   << " phi=" << setw(5) << GoodRecMuon[i2].phi() ;
      tmpdelta=GoodGenMuon[i1].angle(GoodRecMuon[i2]);
      cout << " Delta=" << setw(8) << tmpdelta << "=" 
	   << setw(8) << tmpdelta*180./M_PI ;
      if (tmpdelta < delta) {
	cout << " better" << endl;
	delta = tmpdelta;
	ig=i1; ir=i2;
      }
      else {
	cout << endl;
      }
    }
    float dphimagneticfield;
    // matching range is 15 degrees
    if (delta < 150.*M_PI/180. && GoodRecMuon.size() != 0) {
      dphimagneticfield=0.;
      cout << "Match: GenMuon " << setw(2) << i1+1
	   << " with RecMuon " << setw(2) << ir+1 
	   << " dTheta=" << setw(8) << 180./M_PI*fabs(GoodGenMuon[ig].theta()-
						      GoodRecMuon[ir].theta())
	   << " dPhi=" << setw(8) << 180./M_PI*fabs(GoodGenMuon[ig].phi()-
						    GoodRecMuon[ir].phi()) 
	   << "(" << GoodGenMuon[ig].phi() << "/" << GoodRecMuon[ir].phi() << ")"
	   << " dPsi=" << setw(8) << delta-dphimagneticfield*M_PI/180.
	   << endl;
      nmatches++;
      matchmuon++;
      if (GoodGenMuon[ig].mag() != 0) 
	EMuRes=(GoodGenMuon[ig].mag()-GoodRecMuon[ir].mag())/GoodGenMuon[ig].mag();
      UserHists->FillMatching(180./M_PI*fabs(GoodGenMuon[ig].phi()-
					     GoodRecMuon[ir].phi())-dphimagneticfield,
			      180./M_PI*fabs(GoodGenMuon[ig].theta()-
					     GoodRecMuon[ir].theta())-dphimagneticfield,
			      delta*180./M_PI-dphimagneticfield,1);
    }
    delta=999.;
    for ( i2=0; i2<GoodCluster.size(); i2++) {
      cout << "CL=" << setw(2) << i2 
	   << " theta=" << setw(5) << GoodCluster[i2].theta()
	   << " phi=" << setw(5) << GoodCluster[i2].phi() ;
      tmpdelta=GoodGenMuon[i1].angle(GoodCluster[i2]);
      cout << " Delta=" << setw(8) << tmpdelta << "=" 
	   << setw(8) << tmpdelta*180./M_PI ;
      if (tmpdelta < delta) {
	cout << " better" << endl;
	delta = tmpdelta;
	ig=i1; ir=i2;
      }
      else {
	cout << endl;
      }
    }
    // matching range is 15 degrees
    if (delta < 150.*M_PI/180. && GoodCluster.size() != 0) {
      dphimagneticfield=0.;
      cout << "Match: GenMuon " << setw(2) << i1+1
	   << " with Cluster " << setw(2) << ir+1 
	   << " dTheta=" << setw(8) << 180./M_PI*fabs(GoodGenMuon[ig].theta()-
						      GoodCluster[ir].theta())
	   << " dPhi=" << setw(8) << 180./M_PI*fabs(GoodGenMuon[ig].phi()-
						    GoodCluster[ir].phi()) 
	   << "(" << GoodGenMuon[ig].phi() << "/" << GoodCluster[ir].phi() << ")"
	   << " dPsi=" << setw(8) << delta-dphimagneticfield*M_PI/180.
	   << endl;
      matchcluster++;
      UserHists->FillMatching(180./M_PI*fabs(GoodGenMuon[ig].phi()-
					     GoodCluster[ir].phi())-dphimagneticfield,
			      0.,
			      delta*180./M_PI-dphimagneticfield,2);
    }
    // match Tracks with Muons at Ecal radius
    delta=999.;
    for ( i2=0; i2<GoodTracksAtVertex.size(); i2++) {
      cout << "Tk=" << setw(2) << i2 
	   << " theta=" << setw(5) << GoodTracksAtVertex[i2].theta()
	   << " phi=" << setw(5) << GoodTracksAtVertex[i2].phi() ;
      tmpdelta=GoodGenMuon[i1].angle(GoodTracksAtEcal[i2]);
      cout << " Delta=" << setw(8) << tmpdelta << "=" 
	   << setw(8) << tmpdelta*180./M_PI ;
      if (tmpdelta < delta) {
	cout << " better" << endl;
	delta = tmpdelta;
	ig=i1; ir=i2;
      }
      else {
	cout << endl;
      }
    }
    // matching range is 15 degrees
    if (delta < 150.*M_PI/180. && GoodTracksAtVertex.size() != 0) {
      // recalculate difference to Generator Muon at vertex
      delta = GoodGenMuon[ig].angle(GoodTracksAtVertex[ir]);
      cout << "Match: GenMuon " << setw(2) << i1+1
	   << " with TracksAtVertex " << setw(2) << ir+1 
	   << " dTheta=" << setw(8) << 180./M_PI*fabs(GoodGenMuon[ig].theta()-
						      GoodTracksAtVertex[ir].theta())
	   << " dPhi=" << setw(8) << 180./M_PI*fabs(GoodGenMuon[ig].phi()-
						    GoodTracksAtVertex[ir].phi()) 
	   << "(" << GoodGenMuon[ig].phi() << "/" << GoodTracksAtVertex[ir].phi() << ")"
	   << " dPsi=" << setw(8) << delta*180./M_PI
	   << endl;
      matchtrack++;
      if (GoodGenMuon[ig].mag() != 0) 
	ETkRes=(GoodGenMuon[ig].mag()-GoodTracksAtVertex[ir].mag())/GoodGenMuon[ig].mag();
      UserHists->FillMatching(180./M_PI*fabs(GoodGenMuon[ig].phi()-
					     GoodTracksAtVertex[ir].phi()),
			      180./M_PI*fabs(GoodGenMuon[ig].theta()-
					     GoodTracksAtVertex[ir].theta()),
			      delta*180./M_PI,3);
    }
    UserHists->FillEres(EMuRes,ETkRes);
  }
  UserHists->FillStatistics(GoodGenMuon.size(),GoodCluster.size(),GoodTracksAtVertex.size(),
			    matchmuon,matchcluster,matchtrack);

  cout << "--- done ----------------------------------" << endl;
}

void ExCompMuonInfo::getGeneratorParticles(G3EventProxy * ev) {
  Hep3Vector vtmp(1.,1.,1.);
  cout << "Generator level  ---------------------------------" << endl;
  const RawParticle * aStableParticle;
  HepEventG3EventProxyReader EventReader(ev);
  RawHepEvent MyHepEvent;
  MyHepEvent.read(&EventReader);
  cout << "This Event has " << MyHepEvent.nStable() << " particles." << endl;
  GoodGenMuon.clear();
  for(unsigned int i=1; i<=MyHepEvent.nStable(); i++) {
    aStableParticle = MyHepEvent.getStableParticle(i);
    //      cout << "Particle: " << (*aStableParticle) << endl;
    if ( aStableParticle != 0 
	 && abs(aStableParticle->pid()) == 13) {
      if ( fabs(aStableParticle->eta())<1.) 
	cout << "A barrel muon:  E=" ;
      else
	cout << "An endcap muon: E=" ;
      cout.setf(ios::showpoint);
      cout << setw(8) << setprecision(3) << aStableParticle->e()
	   << ", theta="
	   << setw(8) << setprecision(4) << aStableParticle->theta()
	   << ", phi="
	   << setw(8) << setprecision(4) << aStableParticle->phi()
	   << "="<< aStableParticle->phi()*180/3.1416  
	   << endl;
      vtmp.setMag(   (HepDouble) aStableParticle->e() );
      vtmp.setTheta( (HepDouble) aStableParticle->theta() );
      vtmp.setPhi(   (HepDouble) aStableParticle->phi() );
      GoodGenMuon.push_back(vtmp);
    }
  }

}

void ExCompMuonInfo::getCalorimeterClusters(G3EventProxy * ev) {
  Hep3Vector vtmp(1.,1.,1.);
  cout << "ECAL cluster search ---------------------------------" << endl;
  RecItr<CaloCluster> MyCluster(ev->recEvent());
  GoodCluster.clear();
  while (MyCluster.next()) {
    ncluster++;
    cout << setw(3) << ncluster <<": E=" ;
    cout.setf(ios::showpoint);
    cout << setw(8) << setprecision(3) << MyCluster->Energy()
	 << ", theta="
	 << setw(8) << setprecision(4) << MyCluster->Theta()
	 << ", phi="
	 << setw(8) << setprecision(4) << MyCluster->Phi()  
	 << "="<< MyCluster->Phi()*180/3.1416  
	 << endl;
    if (fabs(MyCluster->Eta()) < 1.156) {
      vtmp.setMag(   (HepDouble) MyCluster->Energy() );
      vtmp.setTheta( (HepDouble) MyCluster->Theta() );
      vtmp.setPhi(   (HepDouble) MyCluster->Phi() );
      GoodCluster.push_back(vtmp);
    }
  }
}

void ExCompMuonInfo::getTrackerTracks(G3EventProxy * ev) {
  cout << "Tracker track search ------------------------------" << endl;
  // instantiate the TrackFinder
  if (myTrackFinder==0) {
    cout << "Create Track reconstructor" << endl;
    myTrackFinder = new TrackReconstructor(new CombinatorialTrackFinder, "TrackerReconstructor");
  }
  Hep3Vector vtmp(1.,1.,1.);
  RecItr<TTrack> MyTrack(ev->recEvent(),"TrackerReconstructor");
  ntracks=0;
  GoodTracksAtVertex.clear();
  GoodTracksAtEcal.clear();
  string volumeName = "EBRY";
  GeaneWrapper geaneWrapper;
  while (MyTrack.next()) {
    ntracks++;
    //FreeTrajectoryState ip = (*MyTrack).impactPointState();
    //cout << "After ip= "  << endl << flush;
    TrajectoryStateOnSurface ipos = (*MyTrack).impactPointStateOnSurface();
    //    if (ipos.isValsimTrigger()->id()) {
    if (ipos.isValid()) {
      FreeTrajectoryState ip = *(ipos.freeTrajectoryState());
      //      FreeTrajectoryState ts = (*MyTrack).stateOnVolume(volumeName);
      FreeTrajectoryState ts = geaneWrapper.propagate(*((*MyTrack).outermostState().freeTrajectoryState()),volumeName);
      cout << setw(3) << ntracks <<": p=" ;
      cout.setf(ios::showpoint);
      cout << setw(8) << setprecision(3) << ip.momentum().mag()
	   << ", theta="
	   << setw(8) << setprecision(4) << ip.momentum().theta()
	   << ", phi="
	   << setw(8) << setprecision(4) << ip.momentum().phi()
	   << "="<< ip.momentum().phi()*180/3.1416  
	   << ", eta="
	   << setw(8) << setprecision(4) << -log(tan(ip.momentum().theta()/2.))
	   << " // ";
      if (ts.momentum().mag() > 0.00001) {
	cout << setw(8) << setprecision(3) << ts.momentum().mag()
	     << ", theta="
	     << setw(8) << setprecision(4) << ts.momentum().theta()
	     << ", phi="
	     << setw(8) << setprecision(4) << ts.momentum().phi()
	     << "="<< ts.momentum().phi()*180/3.1416  
	     << ", eta="
	     << setw(8) << setprecision(4) << -log(tan(ts.momentum().theta()/2.));
	if (abs(log(tan(ts.momentum().theta()/2.)))<1.156) {
	  cout << " is good";
	  vtmp.setMag(ip.momentum().mag());
	  vtmp.setPhi(ip.momentum().phi());
	  vtmp.setTheta(ip.momentum().theta());
	  GoodTracksAtVertex.push_back(vtmp);
	  vtmp.setMag(ts.momentum().mag());
	  vtmp.setPhi(ts.momentum().phi());
	  vtmp.setTheta(ts.momentum().theta());
	  GoodTracksAtEcal.push_back(vtmp);
	}
      }
      else {
	cout << " is bad";
      }
      cout << endl;
    }
  }
}

void ExCompMuonInfo::getMuonTracks(G3EventProxy * ev) {
  Hep3Vector vtmp(1.,1.,1.);
  cout << "Muon track search ---------------------------------" << endl;
  RecItr<TTrack> MyMuonTrack(ev->recEvent(),"MuonReconstructor");
  GoodRecMuon.clear();
  while (MyMuonTrack.next()) {
    nmuons++;
    cout << setw(3) << nmuons  ;
    TrajectoryStateOnSurface ipos = (*MyMuonTrack).impactPointStateOnSurface();
    if (ipos.isValid()) {
      FreeTrajectoryState ip = *(ipos.freeTrajectoryState());
      cout << setw(3) << ntracks <<": p=" ;
      cout.setf(ios::showpoint);
      cout << setw(8) << setprecision(3) << ip.momentum().mag()
	   << ", theta="
	   << setw(8) << setprecision(4) << ip.momentum().theta()
	   << ", phi="
	   << setw(8) << setprecision(4) << ip.momentum().phi()
	   << "="<< ip.momentum().phi()*180/3.1416  
	   << ", eta="
	   << setw(8) << setprecision(4) << -log(tan(ip.momentum().theta()/2.))
	   << " // ";
      vtmp.setMag( (HepDouble)  ip.momentum().mag() );
      vtmp.setTheta( (HepDouble) ip.momentum().theta() ) ;
      vtmp.setPhi( (HepDouble) ip.momentum().phi() );
      GoodRecMuon.push_back(vtmp);
      cout << endl;
    }
    cout << endl;
  }    
}


file  interface/ExCompMuonHistos.h
==================================



#ifndef ExCompMuonHistos_h
#define ExCompMuonHistos_h
/** \file Examples/CompGenRec/interface/ExCompMuonHistos.h
 * Histograms for class ExCompMuonInfo */

#include "Utilities/CHBook4/interface/HB4Histogrammer.h"
#include "Utilities/CHBook4/interface/CHBookHisto.h"

#include <memory>

/** Class to handle histograms for ExCompMuonInfo .
 * \author Stephan Wynhoff, CERN  */
class ExCompMuonHistos : private HB4Histogrammer {
  typedef auto_ptr<CHBookHisto> HP;
public:
  /** The constructor with the filename as argument */
  ExCompMuonHistos() : HB4Histogrammer("comp_muon.hbook")
  { BaseInit(); }
  /** The destructor save the histograms to disk */
  ~ExCompMuonHistos() {  BaseEnd();}
  // the public functions to access the user histograms
  /** For successfull matches fill histograms with difference \a deta in 
   * pseudorapidity, \a dphi in azimuthal angle and \a dpsi in 3D. \a which
   * determins whether tracker tracks (\a which=1), cluster (\a which=2)
   * or muon tracks (\a which=3) are compared with the generator muons. */
  void FillMatching(float dphi, float deta, float dpsi, int which);
  /** Do some statistics on the number of found and matched muon tracks 
   *  (\a hasmuon ,\a matchmuon ), calorimetric clusters (\a hascluster ,
   *  \a matchcluster ) and tracker tracks (\a hastrack ,\a matchtrack ).
   */
  void FillStatistics(int hasmuon, int hascluster, int hastrack, 
		      int matchmuon, int matchcluster, int matchtrack);
  /** Fill histograms for the energy resolution comparing generator information 
   *  with tracks from muon chambers (\a mures ) and from central tracker
   *  (\a tkres ). */
  void FillEres(float mures, float tkres);
private:
  virtual void book();  //<! HBOOK user booking
  virtual void hend();  //<! HBOOK user ending
private:
  // these are accessable only via public methods
  HP hStatistics;
  HP hEMuRes, hETkRes;
  HP hMgMrDeltaPhi, hMgMrDeltaEta, hMgMrDeltaPsi;
  HP hMgClDeltaPhi, hMgClDeltaEta, hMgClDeltaPsi;
  HP hMgTkDeltaPhi, hMgTkDeltaEta, hMgTkDeltaPsi;
};

#endif

file  ExCompMuonHistos.cc
=========================




#include "Utilities/Configuration/interface/Architecture.h"

#include "Examples/CompGenRec/interface/ExCompMuonHistos.h"
#include <iostream>

void ExCompMuonHistos::FillMatching(float dphi, float deta, float dpsi, int which) {
  switch (which) {
  case 1: hMgMrDeltaPhi->fill(dphi);
    hMgMrDeltaEta->fill(deta);
    hMgMrDeltaPsi->fill(dpsi);
    break;
  case 2: hMgClDeltaPhi->fill(dphi);
    hMgClDeltaEta->fill(deta);
    hMgClDeltaPsi->fill(dpsi);
    break;
  case 3: hMgTkDeltaPhi->fill(dphi);
    hMgTkDeltaEta->fill(deta);
    hMgTkDeltaPsi->fill(dpsi);
    break;
  default:
    break;
  }
}

void ExCompMuonHistos::FillStatistics(int hasmuon, int hascluster, int hastrack, 
		      int matchmuon, int matchcluster, int matchtrack) {
  hStatistics->fill(1.,1.);
  hStatistics->fill(2.,hasmuon);
  hStatistics->fill(3.,matchmuon);
  hStatistics->fill(4.,hascluster);
  hStatistics->fill(5.,matchcluster);
  hStatistics->fill(6.,hastrack);
  hStatistics->fill(7.,matchtrack);
}

void ExCompMuonHistos::FillEres(float mures, float tkres) {
  hEMuRes->fill(100.*mures);
  hETkRes->fill(100.*tkres);
}

// here the user histograms are booked
void ExCompMuonHistos::book() {
  cout << "=== booking user histograms ===" << endl;
  hStatistics = HP(new CHBookHisto(1,"Statistics",10,0.5,10.5));

  hEMuRes = HP(new CHBookHisto(2,"Mu(Egen-Erec)/Egen \"MYN#",300,-100.,200.));
  hETkRes = HP(new CHBookHisto(3,"Tk(Egen-Erec)/Egen \"MYN#",300,-10.,20.));

  hMgMrDeltaPhi = HP(new CHBookHisto(10,"Muon dPhi(Gen,Rec) \"M#degree\"N#",250,-5.,50.));
  hMgMrDeltaEta = HP(new CHBookHisto(20,"Muon dTheta(Gen,Rec) \"M#degree\"N#",250,-5.,50.));
  hMgMrDeltaPsi = HP(new CHBookHisto(30,"Muon dPsi(Gen,Rec) \"M#degree\"N#",125,-5.,20.));
  hMgClDeltaPhi = HP(new CHBookHisto(11,"dPhi(GenMu,RecCl) \"M#degree\"N#",100,-2.,25.));
  hMgClDeltaEta = HP(new CHBookHisto(21,"dEta(GenMu,RecCl)",100,-1.5,1.5));
  hMgClDeltaPsi = HP(new CHBookHisto(31,"dPsi(GenMu,RecCl) \"M#degree\"N#",100,-1.,25.));
  hMgTkDeltaPhi = HP(new CHBookHisto(12,"dPhi(GenMu,RecTk) \"M#degree\"N#",90,-.05,.25));
  hMgTkDeltaEta = HP(new CHBookHisto(22,"dTheta(GenMu,RecTk)",90,-.05,.25));
  hMgTkDeltaPsi = HP(new CHBookHisto(32,"dPsi(GenMu,RecTk) \"M#degree\"MN#",125,-.05,.2));
}

// here one can print the histograms or 
// manipulate them before they are written to file
void ExCompMuonHistos::hend() {
  cout << "=== save user histograms ===" << endl;
}