#include <Alignment/LaserAlignment/plugins/TkLasBeamFitter.cc>

Inheritance diagram for TkLasBeamFitter:

Public Member Functions
virtual void	endRun (edm::Run &run, const edm::EventSetup &setup)
virtual void	produce (edm::Event &event, const edm::EventSetup &setup)
	TkLasBeamFitter (const edm::ParameterSet &config)
	~TkLasBeamFitter ()
Private Member Functions
void	buildTrajectory (TkFittedLasBeam &beam, unsigned int &hit, vector< const GeomDetUnit * > &gd, std::vector< GlobalPoint > &globPtrack, vector< TrajectoryStateOnSurface > &tsosLas, GlobalPoint &globPref)
bool	fitBeam (TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatrix, unsigned int &hitsAtTecPlus, unsigned int &nFitParams, double &offset, double &slope, vector< GlobalPoint > &globPtrack, double &bsAngleParam, double &chi2)
void	fitter (TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatrix, unsigned int &hitsAtTecPlus, unsigned int &nFitParams, std::vector< double > &hitPhi, std::vector< double > &hitPhiError, std::vector< double > &hitZprimeError, double &zMin, double &zMax, double &bsAngleParam, double &offset, double &offsetError, double &slope, double &slopeError, double &phiAtMinusParam, double &phiAtPlusParam, double &atThetaSplitParam)
void	getLasBeams (TkFittedLasBeam &beam, vector< TrajectoryStateOnSurface > &tsosLas)
void	getLasHits (TkFittedLasBeam &beam, const SiStripLaserRecHit2D &hit, vector< const GeomDetUnit * > &gd, vector< GlobalPoint > &globHit, unsigned int &hitsAtTecPlus)
void	globalTrackPoint (TkFittedLasBeam &beam, unsigned int &hit, unsigned int &hitsAtTecPlus, double &trackPhi, double &trackPhiRef, std::vector< GlobalPoint > &globHit, std::vector< GlobalPoint > &globPtrack, GlobalPoint &globPref, std::vector< double > &hitPhiError)
void	trackPhi (TkFittedLasBeam &beam, unsigned int &hit, double &trackPhi, double &trackPhiRef, double &offset, double &slope, double &bsAngleParam, double &phiAtMinusParam, double &phiAtPlusParam, double &atThetaSplitParam, std::vector< GlobalPoint > &globHit)
Static Private Member Functions
static double	atFunction (double x, double par)
static double	tecMinusFunction (double x, double par)
static double	tecPlusFunction (double x, double par)
Private Attributes
bool	fitBeamSplitters_
edm::Service< TFileService >	fs
TH1F *	h_bsAngle
TH2F *	h_bsAngleVsBeam
TH1F *	h_chi2
TH1F *	h_chi2ndof
TH1F *	h_hitX
TH1F *	h_hitXAt
TH1F *	h_hitXTecMinus
TH1F *	h_hitXTecPlus
TH2F *	h_hitXvsZAt
TH2F *	h_hitXvsZTecMinus
TH2F *	h_hitXvsZTecPlus
TH1F *	h_pull
TH1F *	h_res
TH1F *	h_resAt
TH1F *	h_resTecMinus
TH1F *	h_resTecPlus
TH2F *	h_resVsHitAt
TH2F *	h_resVsHitTecMinus
TH2F *	h_resVsHitTecPlus
TH2F *	h_resVsZAt
TH2F *	h_resVsZTecMinus
TH2F *	h_resVsZTecPlus
unsigned int	nAtParameters_
const edm::InputTag	src_
Static Private Attributes
static vector< double >	gBarrelModuleOffset
static vector< double >	gBarrelModuleRadius
static double	gBeamR = 0.0
static double	gBeamSplitterZprime = 0.0
static double	gBeamZ0 = 0.0
static double	gBSparam = 0.0
static bool	gFitBeamSplitters = 0
static unsigned int	gHitsAtTecMinus = 0
static vector< double >	gHitZprime
static bool	gIsInnerBarrel = 0
static float	gTIBparam = 0.097614
static float	gTOBparam = 0.034949

Detailed Description

Original Authors: Gero Flucke/Kolja Kaschube Created: Wed May 6 08:43:02 CEST 2009

Id:: TkLasBeamFitter.cc,v 1.11 2011/09/16 07:43:24 mussgill Exp

Description: Fitting LAS beams with track model and providing TrajectoryStateOnSurface for hits.

Implementation:

TkLasBeamCollection read from edm::Run
all done in endRun(..) to allow a correct sequence with production of TkLasBeamCollection in LaserAlignment::endRun(..)

Definition at line 67 of file TkLasBeamFitter.cc.

Constructor & Destructor Documentation

TkLasBeamFitter::TkLasBeamFitter ( const edm::ParameterSet & config ) [explicit]

Definition at line 183 of file TkLasBeamFitter.cc.

                                                               :
  src_(iConfig.getParameter<edm::InputTag>("src")),
  fitBeamSplitters_(iConfig.getParameter<bool>("fitBeamSplitters")),
  nAtParameters_(iConfig.getParameter<unsigned int>("numberOfFittedAtParameters")),
  h_bsAngle(0), h_hitX(0), h_hitXTecPlus(0), h_hitXTecMinus(0),
  h_hitXAt(0), h_chi2(0), h_chi2ndof(0), h_pull(0), h_res(0), 
  h_resTecPlus(0), h_resTecMinus(0), h_resAt(0),
  h_bsAngleVsBeam(0), h_hitXvsZTecPlus(0), h_hitXvsZTecMinus(0),
  h_hitXvsZAt(0), h_resVsZTecPlus(0), h_resVsZTecMinus(0), h_resVsZAt(0),
  h_resVsHitTecPlus(0), h_resVsHitTecMinus(0), h_resVsHitAt(0)
{
  // declare the products to produce
  this->produces<TkFittedLasBeamCollection, edm::InRun>();
  this->produces<TsosVectorCollection, edm::InRun>();
  
  //now do what ever other initialization is needed
}

TkLasBeamFitter::~TkLasBeamFitter ( )

Definition at line 202 of file TkLasBeamFitter.cc.

{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}

Member Function Documentation

double TkLasBeamFitter::atFunction	(	double *	x,
		double *	par
	)		`[static, private]`

Definition at line 620 of file TkLasBeamFitter.cc.

References gBeamR, gBeamSplitterZprime, gBeamZ0, gBSparam, gFitBeamSplitters, gIsInnerBarrel, gTIBparam, gTOBparam, and z.

Referenced by TrackTransformerForGlobalCosmicMuons::fitter().

{
  double z = x[0]; // 'primed'? -> yes!!!
  // TECminus
  if(z < -gBeamSplitterZprime - 2.0*gBeamZ0){
    return par[0] + par[1] * z;
  }
  // BarrelMinus
  else if(-gBeamSplitterZprime - 2.0*gBeamZ0 < z && z < -gBeamZ0){
    // z value includes module offset from main beam axis
    // TOB
    if(!gIsInnerBarrel){
      return par[0] + par[1] * z + gTOBparam * (par[2] + par[4]);
    }
    // TIB
    else{
      return par[0] + par[1] * z - gTIBparam * (par[2] - par[4]);
    }
  }
  // BarrelPlus
  else if(-gBeamZ0 < z && z < gBeamSplitterZprime){
    // z value includes module offset from main beam axis
    // TOB
    if(!gIsInnerBarrel){
      return par[0] + par[1] * z + gTOBparam * (par[3] - par[4]);
    }
    // TIB
    else{
      return par[0] + par[1] * z - gTIBparam * (par[3] + par[4]);
    }
  }
  // TECplus
  else{
    if(gFitBeamSplitters){
      // par[2] = 2*tan(BeamSplitterAngle/2.0)
      return par[0] + par[1] * z - par[5] * (z - gBeamSplitterZprime)/gBeamR; // BS par: 5, 4, or 2
    }
    else{
      return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime)/gBeamR;
    }
  }
}

void TkLasBeamFitter::buildTrajectory	(	TkFittedLasBeam &	beam,
		unsigned int &	hit,
		vector< const GeomDetUnit * > &	gd,
		std::vector< GlobalPoint > &	globPtrack,
		vector< TrajectoryStateOnSurface > &	tsosLas,
		GlobalPoint &	globPref
	)		`[private]`

Definition at line 843 of file TkLasBeamFitter.cc.

References SurfaceSideDefinition::beforeSurface, fieldHandle, gBeamSplitterZprime, gBeamZ0, gHitZprime, TkLasBeam::isTecInternal(), and edm::ESHandle< T >::product().

Referenced by getLasBeams().

{
  const MagneticField* magneticField = fieldHandle.product();
  GlobalVector trajectoryState;
  
  // TECplus
  if(beam.isTecInternal(1)){
    trajectoryState = GlobalVector(globPref-globPtrack[hit]);
  }
  // TECminus
  else if(beam.isTecInternal(-1)){
    trajectoryState = GlobalVector(globPtrack[hit]-globPref);
  }
  // ATs
  else{
    // TECminus
    if(gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0){
      trajectoryState = GlobalVector(globPtrack[hit]-globPref);
    }
    // TECplus
    else if(gHitZprime[hit] > gBeamSplitterZprime){
      trajectoryState = GlobalVector(globPref-globPtrack[hit]);
    }
    // Barrel
    else{
      trajectoryState = GlobalVector(globPtrack[hit]-globPref);
    }
  }     
//   cout << "trajectory: " << trajectoryState << endl;                 
  const FreeTrajectoryState ftsLas = FreeTrajectoryState(globPtrack[hit],trajectoryState,0,magneticField);
  tsosLas.push_back(TrajectoryStateOnSurface(ftsLas,gd[hit]->surface(),
                                             SurfaceSideDefinition::beforeSurface));
}

void TkLasBeamFitter::endRun	(	edm::Run &	run,
		const edm::EventSetup &	setup
	)		`[virtual]`

Reimplemented from edm::EDProducer.

Definition at line 224 of file TkLasBeamFitter.cc.

References gather_cfg::cout, fieldHandle, fitBeamSplitters_, fs, gBSparam, geometry, edm::EventSetup::get(), edm::Run::getByLabel(), getLasBeams(), gFitBeamSplitters, h_bsAngle, h_bsAngleVsBeam, h_chi2, h_chi2ndof, h_hitX, h_hitXAt, h_hitXTecMinus, h_hitXTecPlus, h_hitXvsZAt, h_hitXvsZTecMinus, h_hitXvsZTecPlus, h_pull, h_res, h_resAt, h_resTecMinus, h_resTecPlus, h_resVsHitAt, h_resVsHitTecMinus, h_resVsHitTecPlus, h_resVsZAt, h_resVsZTecMinus, h_resVsZTecPlus, laserBeams, and edm::Run::put().

{
// }
// // FIXME!
// // Indeed, that should be in endRun(..) - as soon as AlignmentProducer can call
// // the algorithm's endRun correctly!
//
//
// void TkLasBeamFitter::beginRun(edm::Run &run, const edm::EventSetup &setup)
// {

  // book histograms
  h_hitX = fs->make<TH1F>("hitX","local x of LAS hits;local x [cm];N",100,-0.5,0.5);
  h_hitXTecPlus = fs->make<TH1F>("hitXTecPlus","local x of LAS hits in TECplus;local x [cm];N",100,-0.5,0.5);
  h_hitXTecMinus = fs->make<TH1F>("hitXTecMinus","local x of LAS hits in TECminus;local x [cm];N",100,-0.5,0.5);
  h_hitXAt = fs->make<TH1F>("hitXAt","local x of LAS hits in ATs;local x [cm];N",100,-2.5,2.5);
  h_hitXvsZTecPlus = fs->make<TH2F>("hitXvsZTecPlus","local x vs z in TECplus;z [cm];local x [cm]",80,120,280,100,-0.5,0.5);
  h_hitXvsZTecMinus = fs->make<TH2F>("hitXvsZTecMinus","local x vs z in TECMinus;z [cm];local x [cm]",80,-280,-120,100,-0.5,0.5);
  h_hitXvsZAt = fs->make<TH2F>("hitXvsZAt","local x vs z in ATs;z [cm];local x [cm]",200,-200,200,100,-0.5,0.5);
  h_chi2 = fs->make<TH1F>("chi2","#chi^{2};#chi^{2};N",100,0,2000);
  h_chi2ndof = fs->make<TH1F>("chi2ndof","#chi^{2} per degree of freedom;#chi^{2}/N_{dof};N",100,0,300);
  h_pull = fs->make<TH1F>("pull","pulls of #phi residuals;pull;N",50,-10,10);
  h_res = fs->make<TH1F>("res","#phi residuals;#phi_{track} - #phi_{hit} [rad];N",60,-0.0015,0.0015);
  h_resTecPlus = fs->make<TH1F>("resTecPlus","#phi residuals in TECplus;#phi_{track} - #phi_{hit} [rad];N",30,-0.0015,0.0015);
  h_resTecMinus = fs->make<TH1F>("resTecMinus","#phi residuals in TECminus;#phi_{track} - #phi_{hit} [rad];N",60,-0.0015,0.0015);
  h_resAt = fs->make<TH1F>("resAt","#phi residuals in ATs;#phi_{track} - #phi_{hit} [rad];N",30,-0.0015,0.0015);
  h_resVsZTecPlus = fs->make<TH2F>("resVsZTecPlus","phi residuals vs. z in TECplus;z [cm];#phi_{track} - #phi_{hit} [rad]",
                         80,120,280,100,-0.0015,0.0015);
  h_resVsZTecMinus = fs->make<TH2F>("resVsZTecMinus","phi residuals vs. z in TECminus;z [cm];#phi_{track} - #phi_{hit} [rad]",
                          80,-280,-120,100,-0.0015,0.0015);
  h_resVsZAt = fs->make<TH2F>("resVsZAt","#phi residuals vs. z in ATs;N;#phi_{track} - #phi_{hit} [rad]",
                    200,-200,200,100,-0.0015,0.0015);
  h_resVsHitTecPlus = fs->make<TH2F>("resVsHitTecPlus","#phi residuals vs. hits in TECplus;hit no.;#phi_{track} - #phi_{hit} [rad]",
                           144,0,144,100,-0.0015,0.0015);
  h_resVsHitTecMinus = fs->make<TH2F>("resVsHitTecMinus","#phi residuals vs. hits in TECminus;hit no.;#phi_{track} - #phi_{hit} [rad]",
                            144,0,144,100,-0.0015,0.0015);
  h_resVsHitAt = fs->make<TH2F>("resVsHitAt","#phi residuals vs. hits in ATs;hit no.;#phi_{track} - #phi_{hit} [rad]",
                      176,0,176,100,-0.0015,0.0015);
  h_bsAngle = fs->make<TH1F>("bsAngle","fitted beam splitter angle;BS angle [rad];N",40,-0.004,0.004);
  h_bsAngleVsBeam = fs->make<TH2F>("bsAngleVsBeam","fitted beam splitter angle per beam;Beam no.;BS angle [rad]",
                         40,0,300,100,-0.004,0.004);

  // Create output collections - they are parallel.
  // (edm::Ref etc. and thus edm::AssociationVector are not supported for edm::Run...)
  std::auto_ptr<TkFittedLasBeamCollection> fittedBeams(new TkFittedLasBeamCollection);
  // One std::vector<TSOS> for each TkFittedLasBeam:
  std::auto_ptr<TsosVectorCollection> tsosesVec(new TsosVectorCollection);

  // get TkLasBeams, Tracker geometry, magnetic field
  run.getByLabel( "LaserAlignment", "tkLaserBeams", laserBeams );
  setup.get<TrackerDigiGeometryRecord>().get(geometry);
  setup.get<IdealMagneticFieldRecord>().get(fieldHandle);

  // hack for fixed BSparams (ugly!)
//   double bsParams[34] = {-0.000266,-0.000956,-0.001205,-0.000018,-0.000759,0.002554,
//                       0.000465,0.000975,0.001006,0.002027,-0.001263,-0.000763,
//                       -0.001702,0.000906,-0.002120,0.001594,0.000661,-0.000457,
//                       -0.000447,0.000347,-0.002266,-0.000446,0.000659,0.000018,
//                       -0.001630,-0.000324,
//                       // ATs
//                       -999.,-0.001709,-0.002091,-999.,
//                       -0.001640,-999.,-0.002444,-0.002345};

  double bsParams[40] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};

  // beam counter
  unsigned int beamNo(0);
  // fit BS? If false, values from bsParams are taken
  gFitBeamSplitters = fitBeamSplitters_;
  if(fitBeamSplitters_) cout << "Fitting BS!" << endl;
  else cout << "BS fixed, not fitted!" << endl;

  // loop over LAS beams
  for(TkLasBeamCollection::const_iterator iBeam = laserBeams->begin(), iEnd = laserBeams->end();
       iBeam != iEnd; ++iBeam){

    TkFittedLasBeam beam(*iBeam);
    vector<TrajectoryStateOnSurface> tsosLas;

    // set BS param for fit
    gBSparam = bsParams[beamNo];

    // call main function; all other functions are called inside getLasBeams(..)
    this->getLasBeams(beam, tsosLas);
    
    // fill output products
    fittedBeams->push_back(beam);
    tsosesVec->push_back(tsosLas);

//     if(!this->fitBeam(fittedBeams->back(), tsosesVec->back())){
//       edm::LogError("BadFit") 
//       << "Problems fitting TkLasBeam, id " << fittedBeams->back().getBeamId() << ".";
//       fittedBeams->pop_back(); // remove last entry added just before
//       tsosesVec->pop_back();   // dito
//     }

    beamNo++;
  }
  
  // finally, put fitted beams and TSOS vectors into run
  run.put(fittedBeams);
  run.put(tsosesVec);
}

bool TkLasBeamFitter::fitBeam	(	TkFittedLasBeam &	beam,
		AlgebraicSymMatrix &	covMatrix,
		unsigned int &	hitsAtTecPlus,
		unsigned int &	nFitParams,
		double &	offset,
		double &	slope,
		vector< GlobalPoint > &	globPtrack,
		double &	bsAngleParam,
		double &	chi2
	)		`[private]`

Definition at line 881 of file TkLasBeamFitter.cc.

References fitBeamSplitters_, gBarrelModuleOffset, gBeamSplitterZprime, gBeamZ0, gHitsAtTecMinus, gHitZprime, TkLasBeam::isAlignmentTube(), TkLasBeam::isTecInternal(), evf::evtn::offset(), TkFittedLasBeam::setParameters(), and slope.

Referenced by getLasBeams().

{
  // set beam parameters for beam output
  unsigned int paramType(0);
  if(!fitBeamSplitters_) paramType = 1;
  if(beam.isAlignmentTube() && hitsAtTecPlus == 0) paramType = 0;
//   const unsigned int nPedeParams = nFitParams + paramType;

  // test without BS params
  const unsigned int nPedeParams(nFitParams);
//   cout << "number of Pede parameters: " << nPedeParams << endl;

  std::vector<TkFittedLasBeam::Scalar> params(nPedeParams);
  params[0] = offset;
  params[1] = slope;
  // no BS parameter for AT beams without TECplus hits
//   if(beam.isTecInternal() || hitsAtTecPlus > 0) params[2] = bsAngleParam;

  AlgebraicMatrix derivatives(gHitZprime.size(), nPedeParams);
  // fill derivatives matrix with local track derivatives
  for(unsigned int hit = 0; hit < gHitZprime.size(); ++hit){
    
    // d(delta phi)/d(offset) is identical for every hit
    derivatives[hit][0] = 1.0;

    // d(delta phi)/d(slope) and d(delta phi)/d(bsAngleParam) depend on parametrizations
    // TECplus
    if(beam.isTecInternal(1)){
      derivatives[hit][1] = globPtrack[hit].z();
//       if(gHitZprime[hit] < gBeamSplitterZprime){
//      derivatives[hit][2] = 0.0;
//       }
//       else{
//      derivatives[hit][2] = - (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
//       }
    }
    // TECminus
    else if(beam.isTecInternal(-1)){
      derivatives[hit][1] = globPtrack[hit].z();
//       if(gHitZprime[hit] > gBeamSplitterZprime){
//      derivatives[hit][2] = 0.0;
//       }
//       else{
//      derivatives[hit][2] = (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
//       }
    }
    // ATs
    else{
      // TECminus
      if(gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0){
        derivatives[hit][1] = globPtrack[hit].z();
//      if(hitsAtTecPlus > 0){
//        derivatives[hit][2] = 0.0;
//      }
      }
      // TECplus
      else if(gHitZprime[hit] > gBeamSplitterZprime){
        derivatives[hit][1] = globPtrack[hit].z();
//      if(hitsAtTecPlus > 0){
//        derivatives[hit][2] = - (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
//      }
      }
      // Barrel
      else{
        derivatives[hit][1] = globPtrack[hit].z() - gBarrelModuleOffset[hit-gHitsAtTecMinus];
//      if(hitsAtTecPlus > 0){
//        derivatives[hit][2] = 0.0;
//      }
      }
    }
  }

  unsigned int firstFixedParam(covMatrix.num_col()); // FIXME! --> no, is fine!!!
//   unsigned int firstFixedParam = nPedeParams - 1;
//   if(beam.isAlignmentTube() && hitsAtTecPlus == 0) firstFixedParam = nPedeParams;
//   cout << "first fixed parameter: " << firstFixedParam << endl;
  // set fit results
  beam.setParameters(paramType, params, covMatrix, derivatives, firstFixedParam, chi2);

  return true; // return false in case of problems
}

void TkLasBeamFitter::fitter	(	TkFittedLasBeam &	beam,
		AlgebraicSymMatrix &	covMatrix,
		unsigned int &	hitsAtTecPlus,
		unsigned int &	nFitParams,
		std::vector< double > &	hitPhi,
		std::vector< double > &	hitPhiError,
		std::vector< double > &	hitZprimeError,
		double &	zMin,
		double &	zMax,
		double &	bsAngleParam,
		double &	offset,
		double &	offsetError,
		double &	slope,
		double &	slopeError,
		double &	phiAtMinusParam,
		double &	phiAtPlusParam,
		double &	atThetaSplitParam
	)		`[private]`

Referenced by getLasBeams().

void TkLasBeamFitter::getLasBeams	(	TkFittedLasBeam &	beam,
		vector< TrajectoryStateOnSurface > &	tsosLas
	)		`[private]`

Fit 'beam' using info from its base class TkLasBeam and set its parameters. Also fill 'tsoses' with TSOS for each LAS hit.

Definition at line 332 of file TkLasBeamFitter.cc.

References abs, TkLasBeam::begin(), buildTrajectory(), gather_cfg::cout, TkLasBeam::end(), fitBeam(), fitBeamSplitters_, fitter(), gBarrelModuleOffset, gBarrelModuleRadius, gBeamR, gBeamSplitterZprime, gBeamZ0, TkLasBeam::getBeamId(), getLasHits(), gHitsAtTecMinus, gHitZprime, gIsInnerBarrel, globalTrackPoint(), h_bsAngle, h_bsAngleVsBeam, h_chi2, h_chi2ndof, h_hitX, h_hitXAt, h_hitXTecMinus, h_hitXTecPlus, h_hitXvsZAt, h_hitXvsZTecMinus, h_hitXvsZTecPlus, h_pull, h_res, h_resAt, h_resTecMinus, h_resTecPlus, h_resVsHitAt, h_resVsHitTecMinus, h_resVsHitTecPlus, h_resVsZAt, h_resVsZTecMinus, h_resVsZTecPlus, TkLasBeam::isAlignmentTube(), TkLasBeam::isRing6(), TkLasBeam::isTecInternal(), SiStripLaserRecHit2D::localPosition(), nAtParameters_, evf::evtn::offset(), perp(), phi, slope, trackPhi(), PV3DBase< T, PVType, FrameType >::x(), and z.

Referenced by endRun().

{
  cout << "---------------------------------------" << endl;
  cout << "beam id: " << beam.getBeamId() // << " isTec: " << (beam.isTecInternal() ? "Y" : "N") 
       << " isTec+: " << (beam.isTecInternal(1) ? "Y" : "N") << " isTec-: " << (beam.isTecInternal(-1) ? "Y" : "N")
       << " isAt: " << (beam.isAlignmentTube() ? "Y" : "N") << " isR6: " << (beam.isRing6() ? "Y" : "N")
       << endl;
 
  // reset static variables 
  gHitsAtTecMinus = 0;
  gHitZprime.clear();
  gBarrelModuleRadius.clear();
  gBarrelModuleOffset.clear();

  // set right beam radius
  gBeamR = beam.isRing6() ? 84.0 : 56.4;
  
  vector<const GeomDetUnit*> gd;
  vector<GlobalPoint> globHit;
  unsigned int hitsAtTecPlus(0);
  double sumZ(0.);

  // loop over hits
  for( TkLasBeam::const_iterator iHit = beam.begin(); iHit < beam.end(); ++iHit ){
    // iHit is a SiStripLaserRecHit2D

    const SiStripLaserRecHit2D hit(*iHit);

    this->getLasHits(beam, hit, gd, globHit, hitsAtTecPlus);
    sumZ += globHit.back().z();

    // fill histos
    h_hitX->Fill(hit.localPosition().x());
    // TECplus
    if(beam.isTecInternal(1)){
      h_hitXTecPlus->Fill(hit.localPosition().x());
      h_hitXvsZTecPlus->Fill(globHit.back().z(),hit.localPosition().x());
    }
    // TECminus
    else if(beam.isTecInternal(-1)){
      h_hitXTecMinus->Fill(hit.localPosition().x());
      h_hitXvsZTecMinus->Fill(globHit.back().z(),hit.localPosition().x());
    }
    // ATs
    else{
      h_hitXAt->Fill(hit.localPosition().x());
      h_hitXvsZAt->Fill(globHit.back().z(),hit.localPosition().x());
    }  
  }

  gBeamZ0 = sumZ / globHit.size();
  double zMin(0.), zMax(0.);
  // TECplus
  if(beam.isTecInternal(1)){
    gBeamSplitterZprime = 205.75 - gBeamZ0;
    zMin = 120.0 - gBeamZ0;
    zMax = 280.0 - gBeamZ0;
  }
  // TECminus
  else if(beam.isTecInternal(-1)){
    gBeamSplitterZprime = -205.75 - gBeamZ0;
    zMin = -280.0 - gBeamZ0;
    zMax = -120.0 - gBeamZ0;
  }
  // AT
  else{
    gBeamSplitterZprime = 112.3 - gBeamZ0;
    zMin = -200.0 - gBeamZ0;
    zMax = 200.0 - gBeamZ0;
  }

  // fill vectors for fitted quantities
  vector<double> hitPhi, hitPhiError, hitZprimeError;

  for(unsigned int hit = 0; hit < globHit.size(); ++hit){
    hitPhi.push_back(static_cast<double>(globHit[hit].phi()));
    // localPositionError[hit] or assume 0.003, 0.006
    hitPhiError.push_back( 0.003 / globHit[hit].perp());
    // no errors on z, fill with zeros
    hitZprimeError.push_back(0.0);
    // barrel-specific values
    if(beam.isAlignmentTube() && abs(globHit[hit].z()) < 112.3){
      gBarrelModuleRadius.push_back(globHit[hit].perp());
      gBarrelModuleOffset.push_back(gBarrelModuleRadius.back() - gBeamR);
      // TIB/TOB flag
      if(gBarrelModuleOffset.back() < 0.0){
        gIsInnerBarrel = 1;
      }
      else{
        gIsInnerBarrel = 0;
      }
      gHitZprime.push_back(globHit[hit].z() - gBeamZ0 - abs(gBarrelModuleOffset.back())); 
    }
    // non-barrel z'
    else{
      gHitZprime.push_back(globHit[hit].z() - gBeamZ0);
    }
  }

  // number of fit parameters, 3 for TECs (always!); 3, 5, or 6 for ATs
  unsigned int tecParams(3), atParams(0);
  if(nAtParameters_ == 3) atParams = 3;
  else if(nAtParameters_ == 5) atParams = 5;
  else atParams = 6;                            // <-- default value, recommended
  unsigned int nFitParams(0);
  if(!fitBeamSplitters_ || 
     (hitsAtTecPlus == 0 && beam.isAlignmentTube() ) ){
    tecParams = tecParams - 1;
    atParams = atParams - 1;
  }
  if(beam.isTecInternal()){
    nFitParams = tecParams;
  }
  else{
    nFitParams = atParams;
  }
  
  // fit parameter definitions
  double offset(0.), offsetError(0.), slope(0.), slopeError(0.),
    bsAngleParam(0.), phiAtMinusParam(0.), phiAtPlusParam(0.),
    atThetaSplitParam(0.);
  AlgebraicSymMatrix covMatrix;
  if(!fitBeamSplitters_ || (beam.isAlignmentTube() && hitsAtTecPlus == 0)){
    covMatrix = AlgebraicSymMatrix(nFitParams, 1);
  }
  else{
    covMatrix = AlgebraicSymMatrix(nFitParams - 1, 1);
  }

  this->fitter(beam, covMatrix, 
               hitsAtTecPlus, nFitParams, 
               hitPhi, hitPhiError, hitZprimeError, 
               zMin, zMax, bsAngleParam,
               offset, offsetError, slope, slopeError,
               phiAtMinusParam, phiAtPlusParam,
               atThetaSplitParam); 
  
  vector<GlobalPoint> globPtrack;
  GlobalPoint globPref;
  double chi2(0.);

  for(unsigned int hit = 0; hit < gHitZprime.size(); ++hit){
    
    // additional phi value (trackPhiRef) for trajectory calculation
    double trackPhi(0.), trackPhiRef(0.); 
    
    this->trackPhi(beam, hit, trackPhi, trackPhiRef,
                   offset, slope, bsAngleParam, 
                   phiAtMinusParam, phiAtPlusParam,
                   atThetaSplitParam, globHit);

    cout << "track phi = " << trackPhi 
         << ", hit phi = " << hitPhi[hit] 
         << ", zPrime = " << gHitZprime[hit] 
         << " r = " << globHit[hit].perp() << endl;

    this->globalTrackPoint(beam, hit, hitsAtTecPlus, 
                           trackPhi, trackPhiRef, 
                           globHit, globPtrack, globPref, 
                           hitPhiError);

    // calculate residuals = pred - hit (in global phi)
    const double phiResidual = globPtrack[hit].phi() - globHit[hit].phi();
    // pull calculation (FIX!!!)
    const double phiResidualPull = phiResidual / hitPhiError[hit];
    //       sqrt(hitPhiError[hit]*hitPhiError[hit] + 
    //     (offsetError*offsetError + globPtrack[hit].z()*globPtrack[hit].z() * slopeError*slopeError));

    // calculate chi2
    chi2 += phiResidual*phiResidual / (hitPhiError[hit]*hitPhiError[hit]);
    
    // fill histos
    h_res->Fill(phiResidual);
    // TECplus
    if(beam.isTecInternal(1)){
      h_pull->Fill(phiResidualPull);
      h_resTecPlus->Fill(phiResidual);
      h_resVsZTecPlus->Fill(globPtrack[hit].z(), phiResidual);
      // Ring 6
      if(beam.isRing6()){
        h_resVsHitTecPlus->Fill(hit+(beam.getBeamId()-1)/10*9+72, phiResidual);
      }
      // Ring 4
      else{
        h_resVsHitTecPlus->Fill(hit+beam.getBeamId()/10*9, phiResidual);
      }
    }
    // TECminus
    else if(beam.isTecInternal(-1)){
      h_pull->Fill(phiResidualPull);
      h_resTecMinus->Fill(phiResidual);
      h_resVsZTecMinus->Fill(globPtrack[hit].z(), phiResidual);
      // Ring 6
      if(beam.isRing6()){
        h_resVsHitTecMinus->Fill(hit+(beam.getBeamId()-101)/10*9+72, phiResidual);
      }
      // Ring 4
      else{
        h_resVsHitTecMinus->Fill(hit+(beam.getBeamId()-100)/10*9, phiResidual);
      }
    }
    // ATs
    else{
      h_pull->Fill(phiResidualPull);
      h_resAt->Fill(phiResidual);
      h_resVsZAt->Fill(globPtrack[hit].z(), phiResidual);
      h_resVsHitAt->Fill(hit+(beam.getBeamId()-200)/10*22, phiResidual);
    }

    this->buildTrajectory(beam, hit, gd, globPtrack, tsosLas, globPref);
  }

  cout << "chi^2 = " << chi2 << ", chi^2/ndof = " << chi2/(gHitZprime.size() - nFitParams) << endl;
  this->fitBeam(beam, covMatrix, hitsAtTecPlus, nFitParams,
                offset, slope, globPtrack, bsAngleParam, chi2);
  
  cout << "bsAngleParam = " << bsAngleParam << endl;

  // fill histos
  // include slope, offset, covariance plots here
  h_chi2->Fill(chi2);
  h_chi2ndof->Fill(chi2/(gHitZprime.size() - nFitParams));
  if(bsAngleParam != 0.0){
    h_bsAngle->Fill(2.0*atan(0.5*bsAngleParam));
    h_bsAngleVsBeam->Fill(beam.getBeamId(), 2.0*atan(0.5*bsAngleParam));
  }
}

void TkLasBeamFitter::getLasHits	(	TkFittedLasBeam &	beam,
		const SiStripLaserRecHit2D &	hit,
		vector< const GeomDetUnit * > &	gd,
		vector< GlobalPoint > &	globHit,
		unsigned int &	hitsAtTecPlus
	)		`[private]`

Definition at line 562 of file TkLasBeamFitter.cc.

References abs, geometry, SiStripLaserRecHit2D::getDetId(), gHitsAtTecMinus, TkLasBeam::isAlignmentTube(), SiStripLaserRecHit2D::localPosition(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by getLasBeams().

{ 
  // get global position of LAS hits
  gd.push_back(geometry->idToDetUnit(hit.getDetId()));
  GlobalPoint globPtemp(gd.back()->toGlobal(hit.localPosition()));

  // testing: globPtemp should be right
  globHit.push_back(globPtemp);

  if(beam.isAlignmentTube()){
    if(abs(globPtemp.z()) > 112.3){
      if(globPtemp.z() < 112.3) gHitsAtTecMinus++ ;
      else hitsAtTecPlus++ ;
    }
  }
}

void TkLasBeamFitter::globalTrackPoint	(	TkFittedLasBeam &	beam,
		unsigned int &	hit,
		unsigned int &	hitsAtTecPlus,
		double &	trackPhi,
		double &	trackPhiRef,
		std::vector< GlobalPoint > &	globHit,
		std::vector< GlobalPoint > &	globPtrack,
		GlobalPoint &	globPref,
		std::vector< double > &	hitPhiError
	)		`[private]`

Definition at line 810 of file TkLasBeamFitter.cc.

References gBeamR, gHitsAtTecMinus, gHitZprime, TkLasBeam::isTecInternal(), perp(), and z.

Referenced by getLasBeams().

{
  // TECs
  if(beam.isTecInternal(0)){
    globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
    globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
  }
  // ATs
  else{
    // TECminus
    if(hit < gHitsAtTecMinus){ // gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0
      globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
      globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
    }
    // TECplus
    else if(hit > gHitZprime.size() - hitsAtTecPlus - 1){ // gHitZprime[hit] > gBeamSplitterZprime
      globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
      globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
    }
    // Barrel
    else{
      globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(globHit[hit].perp(), trackPhi, globHit[hit].z())));
      globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z()));
    }
  }
}

void TkLasBeamFitter::produce	(	edm::Event &	event,
		const edm::EventSetup &	setup
	)		`[virtual]`

Implements edm::EDProducer.

Definition at line 217 of file TkLasBeamFitter.cc.

{
  // Nothing per event!
}

double TkLasBeamFitter::tecMinusFunction	(	double *	x,
		double *	par
	)		`[static, private]`

Definition at line 602 of file TkLasBeamFitter.cc.

References gBeamR, gBeamSplitterZprime, gBSparam, gFitBeamSplitters, and z.

Referenced by TrackTransformerForGlobalCosmicMuons::fitter().

{
  double z = x[0]; // 'primed'? -> yes!!!

  if(z > gBeamSplitterZprime){
    return par[0] + par[1] * z;
  } 
  else{
    if(gFitBeamSplitters){
      // par[2] = 2*tan(BeamSplitterAngle/2.0)
      return par[0] + par[1] * z + par[2] * (z - gBeamSplitterZprime)/gBeamR; 
    }
    else{
      return par[0] + par[1] * z + gBSparam * (z - gBeamSplitterZprime)/gBeamR; 
    }
  }
}

double TkLasBeamFitter::tecPlusFunction	(	double *	x,
		double *	par
	)		`[static, private]`

Definition at line 583 of file TkLasBeamFitter.cc.

References gBeamR, gBeamSplitterZprime, gBSparam, gFitBeamSplitters, and z.

Referenced by TrackTransformerForGlobalCosmicMuons::fitter().

{
  double z = x[0]; // 'primed'? -> yes!!!

  if(z < gBeamSplitterZprime){
    return par[0] + par[1] * z;
  } 
  else{
    if(gFitBeamSplitters){
      // par[2] = 2*tan(BeamSplitterAngle/2.0)
      return par[0] + par[1] * z - par[2] * (z - gBeamSplitterZprime)/gBeamR; 
    }
    else{
      return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime)/gBeamR; 
    }
  }
}

void TkLasBeamFitter::trackPhi	(	TkFittedLasBeam &	beam,
		unsigned int &	hit,
		double &	trackPhi,
		double &	trackPhiRef,
		double &	offset,
		double &	slope,
		double &	bsAngleParam,
		double &	phiAtMinusParam,
		double &	phiAtPlusParam,
		double &	atThetaSplitParam,
		std::vector< GlobalPoint > &	globHit
	)		`[private]`

Definition at line 739 of file TkLasBeamFitter.cc.

References abs, gBarrelModuleOffset, gBeamR, gBeamSplitterZprime, gBeamZ0, gHitsAtTecMinus, gHitZprime, gIsInnerBarrel, gTIBparam, gTOBparam, and TkLasBeam::isTecInternal().

Referenced by getLasBeams().

{
  // TECplus
  if(beam.isTecInternal(1)){
    if(gHitZprime[hit] < gBeamSplitterZprime){
      trackPhi = offset + slope * gHitZprime[hit];
      trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
    }
    else{
      trackPhi = offset + slope * gHitZprime[hit] 
        - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime)/gBeamR;
      trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0)
        - bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime)/gBeamR;
    }
  }
  // TECminus
  else if(beam.isTecInternal(-1)){
    if(gHitZprime[hit] > gBeamSplitterZprime){
      trackPhi = offset + slope * gHitZprime[hit];
      trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
    }
    else{
      trackPhi = offset + slope * gHitZprime[hit] 
        + bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime)/gBeamR;
      trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0)
        + bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime)/gBeamR;
    }
  }
  // ATs
  else{
    // TECminus
    if(gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0){
      trackPhi = offset + slope * gHitZprime[hit];
      trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
    }
    // BarrelMinus
    else if(-gBeamSplitterZprime - 2.0*gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < -gBeamZ0){
      if(!gIsInnerBarrel){
        trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtMinusParam + atThetaSplitParam); 
      }
      else{
        trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtMinusParam - atThetaSplitParam);
      }
      trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit-gHitsAtTecMinus]));
    }
    // BarrelPlus
    else if(-gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < gBeamSplitterZprime){
      if(!gIsInnerBarrel){
        trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtPlusParam - atThetaSplitParam);
      }
      else{
        trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtPlusParam + atThetaSplitParam);
      }
      trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit-gHitsAtTecMinus]));
    }
    // TECplus
    else{
      trackPhi = offset + slope * gHitZprime[hit] 
        - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime)/gBeamR;
      trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0)
        - bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime)/gBeamR;
    }
  }
}