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

Inheritance diagram for TkLasBeamFitter:

Public Member Functions
void	endRunProduce (edm::Run &run, const edm::EventSetup &setup) override

void	produce (edm::Event &event, const edm::EventSetup &setup) override

	TkLasBeamFitter (const edm::ParameterSet &config)

	~TkLasBeamFitter () override

Public Member Functions inherited from edm::one::EDProducer< edm::EndRunProducer >
	EDProducer ()=default

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	hasAbilityToProduceInLumis () const final

bool	hasAbilityToProduceInRuns () const final

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

Public Member Functions inherited from edm::one::EDProducerBase
	EDProducerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDProducerBase () override

Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

std::vector< edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

std::vector< edm::ProductResolverIndex > const &	putTokenIndexToProductResolverIndex () const

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

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

void	resolvePutIndicies (BranchType iBranchType, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)

	~ProducerBase () noexcept(false) override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

virtual	~EDConsumerBase () noexcept(false)

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 = false

static unsigned int	gHitsAtTecMinus = 0

static vector< double >	gHitZprime

static bool	gIsInnerBarrel = false

static float	gTIBparam = 0.097614

static float	gTOBparam = 0.034949

Additional Inherited Members
Public Types inherited from edm::one::EDProducerBase
typedef EDProducerBase	ModuleType

Public Types inherited from edm::ProducerBase
using	ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const , const char , edm::ProductResolverIndex >>

typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

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

Id: TkLasBeamFitter.cc,v 1.12 2012/01/25 17:30:30 innocent 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 65 of file TkLasBeamFitter.cc.

Constructor & Destructor Documentation

TkLasBeamFitter::TkLasBeamFitter ( const edm::ParameterSet & config )

explicit

Definition at line 181 of file TkLasBeamFitter.cc.

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

TkLasBeamFitter::~TkLasBeamFitter ( )

override

Definition at line 200 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 618 of file TkLasBeamFitter.cc.

References gBeamR, gBeamSplitterZprime, gBeamZ0, gBSparam, gFitBeamSplitters, gIsInnerBarrel, gTIBparam, gTOBparam, and 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 841 of file TkLasBeamFitter.cc.

References SurfaceSideDefinition::beforeSurface, gBeamSplitterZprime, gBeamZ0, gHitZprime, TkLasBeam::isTecInternal(), seedCreatorFromRegionConsecutiveHitsEDProducer_cff::magneticField, 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::endRunProduce	(	edm::Run &	run,
		const edm::EventSetup &	setup
	)

override

Definition at line 222 of file TkLasBeamFitter.cc.

References gather_cfg::cout, fitBeamSplitters_, fs, gBSparam, 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, TFileService::make(), eostools::move(), 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...)
   auto fittedBeams = std::make_unique<TkFittedLasBeamCollection>();
   // One std::vector<TSOS> for each TkFittedLasBeam:
   auto tsosesVec = std::make_unique<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(std::move(fittedBeams));
   run.put(std::move(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 879 of file TkLasBeamFitter.cc.

References DEFINE_FWK_MODULE, fitBeamSplitters_, gBarrelModuleOffset, gBeamSplitterZprime, gBeamZ0, gHitsAtTecMinus, gHitZprime, TkLasBeam::isAlignmentTube(), TkLasBeam::isTecInternal(), PFRecoTauDiscriminationByIsolation_cfi::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 663 of file TkLasBeamFitter.cc.

References atFunction(), cuy::col, fitBeamSplitters_, gBSparam, gHitZprime, 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 330 of file TkLasBeamFitter.cc.

References funct::abs(), TkLasBeam::begin(), buildTrajectory(), vertices_cff::chi2, 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_, PFRecoTauDiscriminationByIsolation_cfi::offset, perp(), phi, slope, trackPhi(), z, hfnoseSimHitStudy_cfi::zMax, and hfnoseSimHitStudy_cfi::zMin.

Referenced by endRunProduce().

 {
   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 = true;
       }
       else{
     gIsInnerBarrel = false;
       }
       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 560 of file TkLasBeamFitter.cc.

References funct::abs(), SiStripLaserRecHit2D::getDetId(), gHitsAtTecMinus, TrackerGeometry::idToDetUnit(), 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 808 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
	)

override

Definition at line 215 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 600 of file TkLasBeamFitter.cc.

References gBeamR, gBeamSplitterZprime, gBSparam, gFitBeamSplitters, and 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 581 of file TkLasBeamFitter.cc.

References gBeamR, gBeamSplitterZprime, gBSparam, gFitBeamSplitters, and 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 737 of file TkLasBeamFitter.cc.

References funct::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;
     }
   }
 }