#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 ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

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

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

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
	)

staticprivate

Definition at line 620 of file TkLasBeamFitter.cc.

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

Referenced by 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< class >::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, TFileDirectory::make(), 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

Definition at line 665 of file TkLasBeamFitter.cc.

References asciidump::at, atFunction(), fitBeamSplitters_, gBSparam, gHitZprime, fitWZ::gMinuit, TkLasBeam::isAlignmentTube(), TkLasBeam::isTecInternal(), tecMinusFunction(), and tecPlusFunction().

Referenced by getLasBeams().

 {
   TGraphErrors *lasData = new TGraphErrors(gHitZprime.size(), 
                        &(gHitZprime[0]), &(hitPhi[0]), 
                        &(hitZprimeError[0]), &(hitPhiError[0]));
   
   // do fit (R = entire range)
   if(beam.isTecInternal(1)){
     TF1 tecPlus("tecPlus", tecPlusFunction, zMin, zMax, nFitParams );
     tecPlus.SetParameter( 1, 0 ); // slope
     tecPlus.SetParameter( nFitParams - 1, 0 ); // BS 
     lasData->Fit("tecPlus", "R"); // "R", "RV" or "RQ"
   }
   else if(beam.isTecInternal(-1)){
     TF1 tecMinus("tecMinus", tecMinusFunction, zMin, zMax, nFitParams );
     tecMinus.SetParameter( 1, 0 ); // slope
     tecMinus.SetParameter( nFitParams - 1, 0 ); // BS 
     lasData->Fit("tecMinus", "R");
   }
   else{
     TF1 at("at", atFunction, zMin, zMax, nFitParams );
     at.SetParameter( 1, 0 ); // slope
     at.SetParameter( nFitParams - 1, 0 ); // BS 
     lasData->Fit("at","R");
   }
   
   // get values and errors for offset and slope
   gMinuit->GetParameter(0, offset, offsetError);
   gMinuit->GetParameter(1, slope, slopeError);
 
   // additional AT parameters
   // define param errors that are not used later
   double bsAngleParamError(0.), phiAtMinusParamError(0.),
     phiAtPlusParamError(0.), atThetaSplitParamError(0.);
 
   if(beam.isAlignmentTube()){
     gMinuit->GetParameter(2, phiAtMinusParam, phiAtMinusParamError);
     gMinuit->GetParameter(3, phiAtPlusParam, phiAtPlusParamError);
     gMinuit->GetParameter(4, atThetaSplitParam, atThetaSplitParamError);
   }
   // get Beam Splitter parameters
   if(fitBeamSplitters_){
     if(beam.isAlignmentTube() && hitsAtTecPlus == 0){
       bsAngleParam = gBSparam;
     }
     else{
       gMinuit->GetParameter( nFitParams - 1 , bsAngleParam, bsAngleParamError);
     }
   }
   else{
     bsAngleParam = gBSparam;
   }
 
   // fill covariance matrix
   vector<double> vec( covMatrix.num_col() * covMatrix.num_col() );
   gMinuit->mnemat( &vec[0], covMatrix.num_col() );
   for(int col = 0; col < covMatrix.num_col(); col++){
     for(int row = 0; row < covMatrix.num_col(); row++){
       covMatrix[col][row] = vec[row + covMatrix.num_col()*col];
     }
   }
   // compute correlation between parameters
 //   double corr01 = covMatrix[1][0]/(offsetError*slopeError);
 
   delete lasData;
 }

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 detailsBasic3DVector::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 detailsBasic3DVector::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.

Referenced by JSONExport.JsonExport::export(), HTMLExport.HTMLExport::export(), and HTMLExport.HTMLExportStatic::export().

 {
   // Nothing per event!
 }

double TkLasBeamFitter::tecMinusFunction	(	double *	x,
		double *	par
	)

staticprivate

Definition at line 602 of file TkLasBeamFitter.cc.

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

Referenced by 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
	)

staticprivate

Definition at line 583 of file TkLasBeamFitter.cc.

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

Referenced by 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;
     }
   }
 }