Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:39:17 

0001 /**\class TkLasBeamFitter TkLasBeamFitter.cc Alignment/LaserAlignment/plugins/TkLasBeamFitter.cc
0002 
0003   Original Authors:  Gero Flucke/Kolja Kaschube
0004            Created:  Wed May  6 08:43:02 CEST 2009
0005            $Id: TkLasBeamFitter.cc,v 1.12 2012/01/25 17:30:30 innocent Exp $
0006 
0007  Description: Fitting LAS beams with track model and providing TrajectoryStateOnSurface for hits.
0008 
0009  Implementation:
0010     - TkLasBeamCollection read from edm::Run
0011     - all done in endRun(..) to allow a correct sequence with 
0012       production of TkLasBeamCollection in LaserAlignment::endRun(..)
0013 */
0014 
0015 // framework include files
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/one/EDProducer.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/Framework/interface/Run.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0027 
0028 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0029 #include <Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h>
0030 #include <Geometry/Records/interface/TrackerDigiGeometryRecord.h>
0031 #include <MagneticField/Engine/interface/MagneticField.h>
0032 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0033 
0034 // data formats
0035 // for edm::InRun
0036 #include "DataFormats/Provenance/interface/BranchType.h"
0037 
0038 // laser data formats
0039 #include "DataFormats/Alignment/interface/TkLasBeam.h"
0040 #include "DataFormats/Alignment/interface/TkFittedLasBeam.h"
0041 
0042 // further includes
0043 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0044 #include "Alignment/LaserAlignment/interface/TsosVectorCollection.h"
0045 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0046 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0047 
0048 #include <iostream>
0049 #include "TMinuit.h"
0050 #include "TGraphErrors.h"
0051 #include "TF1.h"
0052 #include "TH1.h"
0053 #include "TH2.h"
0054 
0055 using namespace edm;
0056 using namespace std;
0057 
0058 //
0059 // class declaration
0060 //
0061 
0062 class TkLasBeamFitter : public edm::one::EDProducer<edm::EndRunProducer> {
0063 public:
0064   explicit TkLasBeamFitter(const edm::ParameterSet &config);
0065   ~TkLasBeamFitter() override;
0066 
0067   //virtual void beginJob(const edm::EventSetup& /*access deprecated*/) {}
0068   void produce(edm::Event &event, const edm::EventSetup &setup) override;
0069   // virtual void beginRun(edm::Run &run, const edm::EventSetup &setup);
0070   void endRunProduce(edm::Run &run, const edm::EventSetup &setup) override;
0071   //virtual void endJob() {}
0072 
0073 private:
0074   /// Fit 'beam' using info from its base class TkLasBeam and set its parameters.
0075   /// Also fill 'tsoses' with TSOS for each LAS hit.
0076   void getLasBeams(TkFittedLasBeam &beam, vector<TrajectoryStateOnSurface> &tsosLas);
0077   void getLasHits(TkFittedLasBeam &beam,
0078                   const SiStripLaserRecHit2D &hit,
0079                   vector<const GeomDetUnit *> &gd,
0080                   vector<GlobalPoint> &globHit,
0081                   unsigned int &hitsAtTecPlus);
0082   //   void fillVectors(TkFittedLasBeam &beam);
0083 
0084   // need static functions to be used in fitter(..);
0085   // all parameters used therein have to be static, as well (see below)
0086   static double tecPlusFunction(double *x, double *par);
0087   static double tecMinusFunction(double *x, double *par);
0088   static double atFunction(double *x, double *par);
0089 
0090   void fitter(TkFittedLasBeam &beam,
0091               AlgebraicSymMatrix &covMatrix,
0092               unsigned int &hitsAtTecPlus,
0093               unsigned int &nFitParams,
0094               std::vector<double> &hitPhi,
0095               std::vector<double> &hitPhiError,
0096               std::vector<double> &hitZprimeError,
0097               double &zMin,
0098               double &zMax,
0099               double &bsAngleParam,
0100               double &offset,
0101               double &offsetError,
0102               double &slope,
0103               double &slopeError,
0104               double &phiAtMinusParam,
0105               double &phiAtPlusParam,
0106               double &atThetaSplitParam);
0107 
0108   void trackPhi(TkFittedLasBeam &beam,
0109                 unsigned int &hit,
0110                 double &trackPhi,
0111                 double &trackPhiRef,
0112                 double &offset,
0113                 double &slope,
0114                 double &bsAngleParam,
0115                 double &phiAtMinusParam,
0116                 double &phiAtPlusParam,
0117                 double &atThetaSplitParam,
0118                 std::vector<GlobalPoint> &globHit);
0119 
0120   void globalTrackPoint(TkFittedLasBeam &beam,
0121                         unsigned int &hit,
0122                         unsigned int &hitsAtTecPlus,
0123                         double &trackPhi,
0124                         double &trackPhiRef,
0125                         std::vector<GlobalPoint> &globHit,
0126                         std::vector<GlobalPoint> &globPtrack,
0127                         GlobalPoint &globPref,
0128                         std::vector<double> &hitPhiError);
0129 
0130   void buildTrajectory(TkFittedLasBeam &beam,
0131                        unsigned int &hit,
0132                        vector<const GeomDetUnit *> &gd,
0133                        std::vector<GlobalPoint> &globPtrack,
0134                        vector<TrajectoryStateOnSurface> &tsosLas,
0135                        GlobalPoint &globPref);
0136 
0137   bool fitBeam(TkFittedLasBeam &beam,
0138                AlgebraicSymMatrix &covMatrix,
0139                unsigned int &hitsAtTecPlus,
0140                unsigned int &nFitParams,
0141                double &offset,
0142                double &slope,
0143                vector<GlobalPoint> &globPtrack,
0144                double &bsAngleParam,
0145                double &chi2);
0146 
0147   // ----------member data ---------------------------
0148 
0149   // ES token
0150   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0151   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0152 
0153   // handles
0154   edm::Handle<TkLasBeamCollection> laserBeams;
0155   edm::ESHandle<MagneticField> fieldHandle;
0156   edm::ESHandle<TrackerGeometry> geometry;
0157 
0158   const edm::InputTag src_;
0159   bool fitBeamSplitters_;
0160   unsigned int nAtParameters_;
0161 
0162   edm::Service<TFileService> fs;
0163 
0164   // static parameters used in static parametrization functions
0165   static vector<double> gHitZprime;
0166   static vector<double> gBarrelModuleRadius;
0167   static vector<double> gBarrelModuleOffset;
0168   static float gTIBparam;
0169   static float gTOBparam;
0170   static double gBeamR;
0171   static double gBeamZ0;
0172   static double gBeamSplitterZprime;
0173   static unsigned int gHitsAtTecMinus;
0174   static double gBSparam;
0175   static bool gFitBeamSplitters;
0176   static bool gIsInnerBarrel;
0177 
0178   // histograms
0179   TH1F *h_bsAngle, *h_hitX, *h_hitXTecPlus, *h_hitXTecMinus, *h_hitXAt, *h_chi2, *h_chi2ndof, *h_pull, *h_res,
0180       *h_resTecPlus, *h_resTecMinus, *h_resAt;
0181   TH2F *h_bsAngleVsBeam, *h_hitXvsZTecPlus, *h_hitXvsZTecMinus, *h_hitXvsZAt, *h_resVsZTecPlus, *h_resVsZTecMinus,
0182       *h_resVsZAt, *h_resVsHitTecPlus, *h_resVsHitTecMinus, *h_resVsHitAt;
0183 };
0184 
0185 //
0186 // constants, enums and typedefs
0187 //
0188 
0189 //
0190 // static data member definitions
0191 //
0192 
0193 // static parameters used in parametrization functions
0194 vector<double> TkLasBeamFitter::gHitZprime;
0195 vector<double> TkLasBeamFitter::gBarrelModuleRadius;
0196 vector<double> TkLasBeamFitter::gBarrelModuleOffset;
0197 float TkLasBeamFitter::gTIBparam = 0.097614;  // = abs(r_offset/r_module) (nominal!)
0198 float TkLasBeamFitter::gTOBparam = 0.034949;  // = abs(r_offset/r_module) (nominal!)
0199 double TkLasBeamFitter::gBeamR = 0.0;
0200 double TkLasBeamFitter::gBeamZ0 = 0.0;
0201 double TkLasBeamFitter::gBeamSplitterZprime = 0.0;
0202 unsigned int TkLasBeamFitter::gHitsAtTecMinus = 0;
0203 double TkLasBeamFitter::gBSparam = 0.0;
0204 bool TkLasBeamFitter::gFitBeamSplitters = false;
0205 bool TkLasBeamFitter::gIsInnerBarrel = false;
0206 
0207 //
0208 // constructors and destructor
0209 //
0210 TkLasBeamFitter::TkLasBeamFitter(const edm::ParameterSet &iConfig)
0211     : magFieldToken_(esConsumes<edm::Transition::EndRun>()),
0212       geomToken_(esConsumes<edm::Transition::EndRun>()),
0213       src_(iConfig.getParameter<edm::InputTag>("src")),
0214       fitBeamSplitters_(iConfig.getParameter<bool>("fitBeamSplitters")),
0215       nAtParameters_(iConfig.getParameter<unsigned int>("numberOfFittedAtParameters")),
0216       h_bsAngle(nullptr),
0217       h_hitX(nullptr),
0218       h_hitXTecPlus(nullptr),
0219       h_hitXTecMinus(nullptr),
0220       h_hitXAt(nullptr),
0221       h_chi2(nullptr),
0222       h_chi2ndof(nullptr),
0223       h_pull(nullptr),
0224       h_res(nullptr),
0225       h_resTecPlus(nullptr),
0226       h_resTecMinus(nullptr),
0227       h_resAt(nullptr),
0228       h_bsAngleVsBeam(nullptr),
0229       h_hitXvsZTecPlus(nullptr),
0230       h_hitXvsZTecMinus(nullptr),
0231       h_hitXvsZAt(nullptr),
0232       h_resVsZTecPlus(nullptr),
0233       h_resVsZTecMinus(nullptr),
0234       h_resVsZAt(nullptr),
0235       h_resVsHitTecPlus(nullptr),
0236       h_resVsHitTecMinus(nullptr),
0237       h_resVsHitAt(nullptr) {
0238   // declare the products to produce
0239   this->produces<TkFittedLasBeamCollection, edm::Transition::EndRun>();
0240   this->produces<TsosVectorCollection, edm::Transition::EndRun>();
0241 
0242   //now do what ever other initialization is needed
0243 }
0244 
0245 //---------------------------------------------------------------------------------------
0246 TkLasBeamFitter::~TkLasBeamFitter() {
0247   // do anything here that needs to be done at desctruction time
0248   // (e.g. close files, deallocate resources etc.)
0249 }
0250 
0251 //
0252 // member functions
0253 //
0254 
0255 //---------------------------------------------------------------------------------------
0256 // ------------ method called to produce the data  ------------
0257 void TkLasBeamFitter::produce(edm::Event &iEvent, const edm::EventSetup &setup) {
0258   // Nothing per event!
0259 }
0260 
0261 //---------------------------------------------------------------------------------------
0262 // ------------ method called at end of each run  ---------------------------------------
0263 void TkLasBeamFitter::endRunProduce(edm::Run &run, const edm::EventSetup &setup) {
0264   // }
0265   // // FIXME!
0266   // // Indeed, that should be in endRun(..) - as soon as AlignmentProducer can call
0267   // // the algorithm's endRun correctly!
0268   //
0269   //
0270   // void TkLasBeamFitter::beginRun(edm::Run &run, const edm::EventSetup &setup)
0271   // {
0272 
0273   // book histograms
0274   h_hitX = fs->make<TH1F>("hitX", "local x of LAS hits;local x [cm];N", 100, -0.5, 0.5);
0275   h_hitXTecPlus = fs->make<TH1F>("hitXTecPlus", "local x of LAS hits in TECplus;local x [cm];N", 100, -0.5, 0.5);
0276   h_hitXTecMinus = fs->make<TH1F>("hitXTecMinus", "local x of LAS hits in TECminus;local x [cm];N", 100, -0.5, 0.5);
0277   h_hitXAt = fs->make<TH1F>("hitXAt", "local x of LAS hits in ATs;local x [cm];N", 100, -2.5, 2.5);
0278   h_hitXvsZTecPlus =
0279       fs->make<TH2F>("hitXvsZTecPlus", "local x vs z in TECplus;z [cm];local x [cm]", 80, 120, 280, 100, -0.5, 0.5);
0280   h_hitXvsZTecMinus =
0281       fs->make<TH2F>("hitXvsZTecMinus", "local x vs z in TECMinus;z [cm];local x [cm]", 80, -280, -120, 100, -0.5, 0.5);
0282   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);
0283   h_chi2 = fs->make<TH1F>("chi2", "#chi^{2};#chi^{2};N", 100, 0, 2000);
0284   h_chi2ndof = fs->make<TH1F>("chi2ndof", "#chi^{2} per degree of freedom;#chi^{2}/N_{dof};N", 100, 0, 300);
0285   h_pull = fs->make<TH1F>("pull", "pulls of #phi residuals;pull;N", 50, -10, 10);
0286   h_res = fs->make<TH1F>("res", "#phi residuals;#phi_{track} - #phi_{hit} [rad];N", 60, -0.0015, 0.0015);
0287   h_resTecPlus =
0288       fs->make<TH1F>("resTecPlus", "#phi residuals in TECplus;#phi_{track} - #phi_{hit} [rad];N", 30, -0.0015, 0.0015);
0289   h_resTecMinus = fs->make<TH1F>(
0290       "resTecMinus", "#phi residuals in TECminus;#phi_{track} - #phi_{hit} [rad];N", 60, -0.0015, 0.0015);
0291   h_resAt = fs->make<TH1F>("resAt", "#phi residuals in ATs;#phi_{track} - #phi_{hit} [rad];N", 30, -0.0015, 0.0015);
0292   h_resVsZTecPlus = fs->make<TH2F>("resVsZTecPlus",
0293                                    "phi residuals vs. z in TECplus;z [cm];#phi_{track} - #phi_{hit} [rad]",
0294                                    80,
0295                                    120,
0296                                    280,
0297                                    100,
0298                                    -0.0015,
0299                                    0.0015);
0300   h_resVsZTecMinus = fs->make<TH2F>("resVsZTecMinus",
0301                                     "phi residuals vs. z in TECminus;z [cm];#phi_{track} - #phi_{hit} [rad]",
0302                                     80,
0303                                     -280,
0304                                     -120,
0305                                     100,
0306                                     -0.0015,
0307                                     0.0015);
0308   h_resVsZAt = fs->make<TH2F>(
0309       "resVsZAt", "#phi residuals vs. z in ATs;N;#phi_{track} - #phi_{hit} [rad]", 200, -200, 200, 100, -0.0015, 0.0015);
0310   h_resVsHitTecPlus = fs->make<TH2F>("resVsHitTecPlus",
0311                                      "#phi residuals vs. hits in TECplus;hit no.;#phi_{track} - #phi_{hit} [rad]",
0312                                      144,
0313                                      0,
0314                                      144,
0315                                      100,
0316                                      -0.0015,
0317                                      0.0015);
0318   h_resVsHitTecMinus = fs->make<TH2F>("resVsHitTecMinus",
0319                                       "#phi residuals vs. hits in TECminus;hit no.;#phi_{track} - #phi_{hit} [rad]",
0320                                       144,
0321                                       0,
0322                                       144,
0323                                       100,
0324                                       -0.0015,
0325                                       0.0015);
0326   h_resVsHitAt = fs->make<TH2F>("resVsHitAt",
0327                                 "#phi residuals vs. hits in ATs;hit no.;#phi_{track} - #phi_{hit} [rad]",
0328                                 176,
0329                                 0,
0330                                 176,
0331                                 100,
0332                                 -0.0015,
0333                                 0.0015);
0334   h_bsAngle = fs->make<TH1F>("bsAngle", "fitted beam splitter angle;BS angle [rad];N", 40, -0.004, 0.004);
0335   h_bsAngleVsBeam = fs->make<TH2F>(
0336       "bsAngleVsBeam", "fitted beam splitter angle per beam;Beam no.;BS angle [rad]", 40, 0, 300, 100, -0.004, 0.004);
0337 
0338   // Create output collections - they are parallel.
0339   // (edm::Ref etc. and thus edm::AssociationVector are not supported for edm::Run...)
0340   auto fittedBeams = std::make_unique<TkFittedLasBeamCollection>();
0341   // One std::vector<TSOS> for each TkFittedLasBeam:
0342   auto tsosesVec = std::make_unique<TsosVectorCollection>();
0343 
0344   // get TkLasBeams, Tracker geometry, magnetic field
0345   run.getByLabel("LaserAlignment", "tkLaserBeams", laserBeams);
0346   geometry = setup.getHandle(geomToken_);
0347   fieldHandle = setup.getHandle(magFieldToken_);
0348 
0349   // hack for fixed BSparams (ugly!)
0350   //   double bsParams[34] = {-0.000266,-0.000956,-0.001205,-0.000018,-0.000759,0.002554,
0351   //             0.000465,0.000975,0.001006,0.002027,-0.001263,-0.000763,
0352   //             -0.001702,0.000906,-0.002120,0.001594,0.000661,-0.000457,
0353   //             -0.000447,0.000347,-0.002266,-0.000446,0.000659,0.000018,
0354   //             -0.001630,-0.000324,
0355   //             // ATs
0356   //             -999.,-0.001709,-0.002091,-999.,
0357   //             -0.001640,-999.,-0.002444,-0.002345};
0358 
0359   double bsParams[40] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0360                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
0361 
0362   // beam counter
0363   unsigned int beamNo(0);
0364   // fit BS? If false, values from bsParams are taken
0365   gFitBeamSplitters = fitBeamSplitters_;
0366   if (fitBeamSplitters_)
0367     cout << "Fitting BS!" << endl;
0368   else
0369     cout << "BS fixed, not fitted!" << endl;
0370 
0371   // loop over LAS beams
0372   for (TkLasBeamCollection::const_iterator iBeam = laserBeams->begin(), iEnd = laserBeams->end(); iBeam != iEnd;
0373        ++iBeam) {
0374     TkFittedLasBeam beam(*iBeam);
0375     vector<TrajectoryStateOnSurface> tsosLas;
0376 
0377     // set BS param for fit
0378     gBSparam = bsParams[beamNo];
0379 
0380     // call main function; all other functions are called inside getLasBeams(..)
0381     this->getLasBeams(beam, tsosLas);
0382 
0383     // fill output products
0384     fittedBeams->push_back(beam);
0385     tsosesVec->push_back(tsosLas);
0386 
0387     //     if(!this->fitBeam(fittedBeams->back(), tsosesVec->back())){
0388     //       edm::LogError("BadFit")
0389     //   << "Problems fitting TkLasBeam, id " << fittedBeams->back().getBeamId() << ".";
0390     //       fittedBeams->pop_back(); // remove last entry added just before
0391     //       tsosesVec->pop_back();   // dito
0392     //     }
0393 
0394     beamNo++;
0395   }
0396 
0397   // finally, put fitted beams and TSOS vectors into run
0398   run.put(std::move(fittedBeams));
0399   run.put(std::move(tsosesVec));
0400 }
0401 
0402 // methods for las data processing
0403 
0404 // -------------- loop over beams, call functions ----------------------------
0405 void TkLasBeamFitter::getLasBeams(TkFittedLasBeam &beam, vector<TrajectoryStateOnSurface> &tsosLas) {
0406   cout << "---------------------------------------" << endl;
0407   cout << "beam id: " << beam.getBeamId()  // << " isTec: " << (beam.isTecInternal() ? "Y" : "N")
0408        << " isTec+: " << (beam.isTecInternal(1) ? "Y" : "N") << " isTec-: " << (beam.isTecInternal(-1) ? "Y" : "N")
0409        << " isAt: " << (beam.isAlignmentTube() ? "Y" : "N") << " isR6: " << (beam.isRing6() ? "Y" : "N") << endl;
0410 
0411   // reset static variables
0412   gHitsAtTecMinus = 0;
0413   gHitZprime.clear();
0414   gBarrelModuleRadius.clear();
0415   gBarrelModuleOffset.clear();
0416 
0417   // set right beam radius
0418   gBeamR = beam.isRing6() ? 84.0 : 56.4;
0419 
0420   vector<const GeomDetUnit *> gd;
0421   vector<GlobalPoint> globHit;
0422   unsigned int hitsAtTecPlus(0);
0423   double sumZ(0.);
0424 
0425   // loop over hits
0426   for (TkLasBeam::const_iterator iHit = beam.begin(); iHit < beam.end(); ++iHit) {
0427     // iHit is a SiStripLaserRecHit2D
0428 
0429     const SiStripLaserRecHit2D &hit(*iHit);
0430 
0431     this->getLasHits(beam, hit, gd, globHit, hitsAtTecPlus);
0432     sumZ += globHit.back().z();
0433 
0434     // fill histos
0435     h_hitX->Fill(hit.localPosition().x());
0436     // TECplus
0437     if (beam.isTecInternal(1)) {
0438       h_hitXTecPlus->Fill(hit.localPosition().x());
0439       h_hitXvsZTecPlus->Fill(globHit.back().z(), hit.localPosition().x());
0440     }
0441     // TECminus
0442     else if (beam.isTecInternal(-1)) {
0443       h_hitXTecMinus->Fill(hit.localPosition().x());
0444       h_hitXvsZTecMinus->Fill(globHit.back().z(), hit.localPosition().x());
0445     }
0446     // ATs
0447     else {
0448       h_hitXAt->Fill(hit.localPosition().x());
0449       h_hitXvsZAt->Fill(globHit.back().z(), hit.localPosition().x());
0450     }
0451   }
0452 
0453   gBeamZ0 = sumZ / globHit.size();
0454   double zMin(0.), zMax(0.);
0455   // TECplus
0456   if (beam.isTecInternal(1)) {
0457     gBeamSplitterZprime = 205.75 - gBeamZ0;
0458     zMin = 120.0 - gBeamZ0;
0459     zMax = 280.0 - gBeamZ0;
0460   }
0461   // TECminus
0462   else if (beam.isTecInternal(-1)) {
0463     gBeamSplitterZprime = -205.75 - gBeamZ0;
0464     zMin = -280.0 - gBeamZ0;
0465     zMax = -120.0 - gBeamZ0;
0466   }
0467   // AT
0468   else {
0469     gBeamSplitterZprime = 112.3 - gBeamZ0;
0470     zMin = -200.0 - gBeamZ0;
0471     zMax = 200.0 - gBeamZ0;
0472   }
0473 
0474   // fill vectors for fitted quantities
0475   vector<double> hitPhi, hitPhiError, hitZprimeError;
0476 
0477   for (unsigned int hit = 0; hit < globHit.size(); ++hit) {
0478     hitPhi.push_back(static_cast<double>(globHit[hit].phi()));
0479     // localPositionError[hit] or assume 0.003, 0.006
0480     hitPhiError.push_back(0.003 / globHit[hit].perp());
0481     // no errors on z, fill with zeros
0482     hitZprimeError.push_back(0.0);
0483     // barrel-specific values
0484     if (beam.isAlignmentTube() && abs(globHit[hit].z()) < 112.3) {
0485       gBarrelModuleRadius.push_back(globHit[hit].perp());
0486       gBarrelModuleOffset.push_back(gBarrelModuleRadius.back() - gBeamR);
0487       // TIB/TOB flag
0488       if (gBarrelModuleOffset.back() < 0.0) {
0489         gIsInnerBarrel = true;
0490       } else {
0491         gIsInnerBarrel = false;
0492       }
0493       gHitZprime.push_back(globHit[hit].z() - gBeamZ0 - abs(gBarrelModuleOffset.back()));
0494     }
0495     // non-barrel z'
0496     else {
0497       gHitZprime.push_back(globHit[hit].z() - gBeamZ0);
0498     }
0499   }
0500 
0501   // number of fit parameters, 3 for TECs (always!); 3, 5, or 6 for ATs
0502   unsigned int tecParams(3), atParams(0);
0503   if (nAtParameters_ == 3)
0504     atParams = 3;
0505   else if (nAtParameters_ == 5)
0506     atParams = 5;
0507   else
0508     atParams = 6;  // <-- default value, recommended
0509   unsigned int nFitParams(0);
0510   if (!fitBeamSplitters_ || (hitsAtTecPlus == 0 && beam.isAlignmentTube())) {
0511     tecParams = tecParams - 1;
0512     atParams = atParams - 1;
0513   }
0514   if (beam.isTecInternal()) {
0515     nFitParams = tecParams;
0516   } else {
0517     nFitParams = atParams;
0518   }
0519 
0520   // fit parameter definitions
0521   double offset(0.), offsetError(0.), slope(0.), slopeError(0.), bsAngleParam(0.), phiAtMinusParam(0.),
0522       phiAtPlusParam(0.), atThetaSplitParam(0.);
0523   AlgebraicSymMatrix covMatrix;
0524   if (!fitBeamSplitters_ || (beam.isAlignmentTube() && hitsAtTecPlus == 0)) {
0525     covMatrix = AlgebraicSymMatrix(nFitParams, 1);
0526   } else {
0527     covMatrix = AlgebraicSymMatrix(nFitParams - 1, 1);
0528   }
0529 
0530   this->fitter(beam,
0531                covMatrix,
0532                hitsAtTecPlus,
0533                nFitParams,
0534                hitPhi,
0535                hitPhiError,
0536                hitZprimeError,
0537                zMin,
0538                zMax,
0539                bsAngleParam,
0540                offset,
0541                offsetError,
0542                slope,
0543                slopeError,
0544                phiAtMinusParam,
0545                phiAtPlusParam,
0546                atThetaSplitParam);
0547 
0548   vector<GlobalPoint> globPtrack;
0549   GlobalPoint globPref;
0550   double chi2(0.);
0551 
0552   for (unsigned int hit = 0; hit < gHitZprime.size(); ++hit) {
0553     // additional phi value (trackPhiRef) for trajectory calculation
0554     double trackPhi(0.), trackPhiRef(0.);
0555 
0556     this->trackPhi(beam,
0557                    hit,
0558                    trackPhi,
0559                    trackPhiRef,
0560                    offset,
0561                    slope,
0562                    bsAngleParam,
0563                    phiAtMinusParam,
0564                    phiAtPlusParam,
0565                    atThetaSplitParam,
0566                    globHit);
0567 
0568     cout << "track phi = " << trackPhi << ", hit phi = " << hitPhi[hit] << ", zPrime = " << gHitZprime[hit]
0569          << " r = " << globHit[hit].perp() << endl;
0570 
0571     this->globalTrackPoint(beam, hit, hitsAtTecPlus, trackPhi, trackPhiRef, globHit, globPtrack, globPref, hitPhiError);
0572 
0573     // calculate residuals = pred - hit (in global phi)
0574     const double phiResidual = globPtrack[hit].phi() - globHit[hit].phi();
0575     // pull calculation (FIX!!!)
0576     const double phiResidualPull = phiResidual / hitPhiError[hit];
0577     //       sqrt(hitPhiError[hit]*hitPhiError[hit] +
0578     //     (offsetError*offsetError + globPtrack[hit].z()*globPtrack[hit].z() * slopeError*slopeError));
0579 
0580     // calculate chi2
0581     chi2 += phiResidual * phiResidual / (hitPhiError[hit] * hitPhiError[hit]);
0582 
0583     // fill histos
0584     h_res->Fill(phiResidual);
0585     // TECplus
0586     if (beam.isTecInternal(1)) {
0587       h_pull->Fill(phiResidualPull);
0588       h_resTecPlus->Fill(phiResidual);
0589       h_resVsZTecPlus->Fill(globPtrack[hit].z(), phiResidual);
0590       // Ring 6
0591       if (beam.isRing6()) {
0592         h_resVsHitTecPlus->Fill(hit + (beam.getBeamId() - 1) / 10 * 9 + 72, phiResidual);
0593       }
0594       // Ring 4
0595       else {
0596         h_resVsHitTecPlus->Fill(hit + beam.getBeamId() / 10 * 9, phiResidual);
0597       }
0598     }
0599     // TECminus
0600     else if (beam.isTecInternal(-1)) {
0601       h_pull->Fill(phiResidualPull);
0602       h_resTecMinus->Fill(phiResidual);
0603       h_resVsZTecMinus->Fill(globPtrack[hit].z(), phiResidual);
0604       // Ring 6
0605       if (beam.isRing6()) {
0606         h_resVsHitTecMinus->Fill(hit + (beam.getBeamId() - 101) / 10 * 9 + 72, phiResidual);
0607       }
0608       // Ring 4
0609       else {
0610         h_resVsHitTecMinus->Fill(hit + (beam.getBeamId() - 100) / 10 * 9, phiResidual);
0611       }
0612     }
0613     // ATs
0614     else {
0615       h_pull->Fill(phiResidualPull);
0616       h_resAt->Fill(phiResidual);
0617       h_resVsZAt->Fill(globPtrack[hit].z(), phiResidual);
0618       h_resVsHitAt->Fill(hit + (beam.getBeamId() - 200) / 10 * 22, phiResidual);
0619     }
0620 
0621     this->buildTrajectory(beam, hit, gd, globPtrack, tsosLas, globPref);
0622   }
0623 
0624   cout << "chi^2 = " << chi2 << ", chi^2/ndof = " << chi2 / (gHitZprime.size() - nFitParams) << endl;
0625   this->fitBeam(beam, covMatrix, hitsAtTecPlus, nFitParams, offset, slope, globPtrack, bsAngleParam, chi2);
0626 
0627   cout << "bsAngleParam = " << bsAngleParam << endl;
0628 
0629   // fill histos
0630   // include slope, offset, covariance plots here
0631   h_chi2->Fill(chi2);
0632   h_chi2ndof->Fill(chi2 / (gHitZprime.size() - nFitParams));
0633   if (bsAngleParam != 0.0) {
0634     h_bsAngle->Fill(2.0 * atan(0.5 * bsAngleParam));
0635     h_bsAngleVsBeam->Fill(beam.getBeamId(), 2.0 * atan(0.5 * bsAngleParam));
0636   }
0637 }
0638 
0639 // --------- get hits, convert to global coordinates ---------------------------
0640 void TkLasBeamFitter::getLasHits(TkFittedLasBeam &beam,
0641                                  const SiStripLaserRecHit2D &hit,
0642                                  vector<const GeomDetUnit *> &gd,
0643                                  vector<GlobalPoint> &globHit,
0644                                  unsigned int &hitsAtTecPlus) {
0645   // get global position of LAS hits
0646   gd.push_back(geometry->idToDetUnit(hit.getDetId()));
0647   GlobalPoint globPtemp(gd.back()->toGlobal(hit.localPosition()));
0648 
0649   // testing: globPtemp should be right
0650   globHit.push_back(globPtemp);
0651 
0652   if (beam.isAlignmentTube()) {
0653     if (abs(globPtemp.z()) > 112.3) {
0654       if (globPtemp.z() < 112.3)
0655         gHitsAtTecMinus++;
0656       else
0657         hitsAtTecPlus++;
0658     }
0659   }
0660 }
0661 
0662 // ------------ parametrization functions for las beam fits ------------
0663 double TkLasBeamFitter::tecPlusFunction(double *x, double *par) {
0664   double z = x[0];  // 'primed'? -> yes!!!
0665 
0666   if (z < gBeamSplitterZprime) {
0667     return par[0] + par[1] * z;
0668   } else {
0669     if (gFitBeamSplitters) {
0670       // par[2] = 2*tan(BeamSplitterAngle/2.0)
0671       return par[0] + par[1] * z - par[2] * (z - gBeamSplitterZprime) / gBeamR;
0672     } else {
0673       return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime) / gBeamR;
0674     }
0675   }
0676 }
0677 
0678 double TkLasBeamFitter::tecMinusFunction(double *x, double *par) {
0679   double z = x[0];  // 'primed'? -> yes!!!
0680 
0681   if (z > gBeamSplitterZprime) {
0682     return par[0] + par[1] * z;
0683   } else {
0684     if (gFitBeamSplitters) {
0685       // par[2] = 2*tan(BeamSplitterAngle/2.0)
0686       return par[0] + par[1] * z + par[2] * (z - gBeamSplitterZprime) / gBeamR;
0687     } else {
0688       return par[0] + par[1] * z + gBSparam * (z - gBeamSplitterZprime) / gBeamR;
0689     }
0690   }
0691 }
0692 
0693 double TkLasBeamFitter::atFunction(double *x, double *par) {
0694   double z = x[0];  // 'primed'? -> yes!!!
0695   // TECminus
0696   if (z < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
0697     return par[0] + par[1] * z;
0698   }
0699   // BarrelMinus
0700   else if (-gBeamSplitterZprime - 2.0 * gBeamZ0 < z && z < -gBeamZ0) {
0701     // z value includes module offset from main beam axis
0702     // TOB
0703     if (!gIsInnerBarrel) {
0704       return par[0] + par[1] * z + gTOBparam * (par[2] + par[4]);
0705     }
0706     // TIB
0707     else {
0708       return par[0] + par[1] * z - gTIBparam * (par[2] - par[4]);
0709     }
0710   }
0711   // BarrelPlus
0712   else if (-gBeamZ0 < z && z < gBeamSplitterZprime) {
0713     // z value includes module offset from main beam axis
0714     // TOB
0715     if (!gIsInnerBarrel) {
0716       return par[0] + par[1] * z + gTOBparam * (par[3] - par[4]);
0717     }
0718     // TIB
0719     else {
0720       return par[0] + par[1] * z - gTIBparam * (par[3] + par[4]);
0721     }
0722   }
0723   // TECplus
0724   else {
0725     if (gFitBeamSplitters) {
0726       // par[2] = 2*tan(BeamSplitterAngle/2.0)
0727       return par[0] + par[1] * z - par[5] * (z - gBeamSplitterZprime) / gBeamR;  // BS par: 5, 4, or 2
0728     } else {
0729       return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime) / gBeamR;
0730     }
0731   }
0732 }
0733 
0734 // ------------ perform fit of beams ------------------------------------
0735 void TkLasBeamFitter::fitter(TkFittedLasBeam &beam,
0736                              AlgebraicSymMatrix &covMatrix,
0737                              unsigned int &hitsAtTecPlus,
0738                              unsigned int &nFitParams,
0739                              vector<double> &hitPhi,
0740                              vector<double> &hitPhiError,
0741                              vector<double> &hitZprimeError,
0742                              double &zMin,
0743                              double &zMax,
0744                              double &bsAngleParam,
0745                              double &offset,
0746                              double &offsetError,
0747                              double &slope,
0748                              double &slopeError,
0749                              double &phiAtMinusParam,
0750                              double &phiAtPlusParam,
0751                              double &atThetaSplitParam) {
0752   TGraphErrors *lasData =
0753       new TGraphErrors(gHitZprime.size(), &(gHitZprime[0]), &(hitPhi[0]), &(hitZprimeError[0]), &(hitPhiError[0]));
0754 
0755   // do fit (R = entire range)
0756   if (beam.isTecInternal(1)) {
0757     TF1 tecPlus("tecPlus", tecPlusFunction, zMin, zMax, nFitParams);
0758     tecPlus.SetParameter(1, 0);               // slope
0759     tecPlus.SetParameter(nFitParams - 1, 0);  // BS
0760     lasData->Fit(&tecPlus, "R");              // "R", "RV" or "RQ"
0761   } else if (beam.isTecInternal(-1)) {
0762     TF1 tecMinus("tecMinus", tecMinusFunction, zMin, zMax, nFitParams);
0763     tecMinus.SetParameter(1, 0);               // slope
0764     tecMinus.SetParameter(nFitParams - 1, 0);  // BS
0765     lasData->Fit(&tecMinus, "R");
0766   } else {
0767     TF1 at("at", atFunction, zMin, zMax, nFitParams);
0768     at.SetParameter(1, 0);               // slope
0769     at.SetParameter(nFitParams - 1, 0);  // BS
0770     lasData->Fit(&at, "R");
0771   }
0772 
0773   // get values and errors for offset and slope
0774   gMinuit->GetParameter(0, offset, offsetError);
0775   gMinuit->GetParameter(1, slope, slopeError);
0776 
0777   // additional AT parameters
0778   // define param errors that are not used later
0779   double bsAngleParamError(0.), phiAtMinusParamError(0.), phiAtPlusParamError(0.), atThetaSplitParamError(0.);
0780 
0781   if (beam.isAlignmentTube()) {
0782     gMinuit->GetParameter(2, phiAtMinusParam, phiAtMinusParamError);
0783     gMinuit->GetParameter(3, phiAtPlusParam, phiAtPlusParamError);
0784     gMinuit->GetParameter(4, atThetaSplitParam, atThetaSplitParamError);
0785   }
0786   // get Beam Splitter parameters
0787   if (fitBeamSplitters_) {
0788     if (beam.isAlignmentTube() && hitsAtTecPlus == 0) {
0789       bsAngleParam = gBSparam;
0790     } else {
0791       gMinuit->GetParameter(nFitParams - 1, bsAngleParam, bsAngleParamError);
0792     }
0793   } else {
0794     bsAngleParam = gBSparam;
0795   }
0796 
0797   // fill covariance matrix
0798   vector<double> vec(covMatrix.num_col() * covMatrix.num_col());
0799   gMinuit->mnemat(&vec[0], covMatrix.num_col());
0800   for (int col = 0; col < covMatrix.num_col(); col++) {
0801     for (int row = 0; row < covMatrix.num_col(); row++) {
0802       covMatrix[col][row] = vec[row + covMatrix.num_col() * col];
0803     }
0804   }
0805   // compute correlation between parameters
0806   //   double corr01 = covMatrix[1][0]/(offsetError*slopeError);
0807 
0808   delete lasData;
0809 }
0810 
0811 // -------------- calculate track phi value ----------------------------------
0812 void TkLasBeamFitter::trackPhi(TkFittedLasBeam &beam,
0813                                unsigned int &hit,
0814                                double &trackPhi,
0815                                double &trackPhiRef,
0816                                double &offset,
0817                                double &slope,
0818                                double &bsAngleParam,
0819                                double &phiAtMinusParam,
0820                                double &phiAtPlusParam,
0821                                double &atThetaSplitParam,
0822                                vector<GlobalPoint> &globHit) {
0823   // TECplus
0824   if (beam.isTecInternal(1)) {
0825     if (gHitZprime[hit] < gBeamSplitterZprime) {
0826       trackPhi = offset + slope * gHitZprime[hit];
0827       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
0828     } else {
0829       trackPhi = offset + slope * gHitZprime[hit] - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime) / gBeamR;
0830       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0) -
0831                     bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime) / gBeamR;
0832     }
0833   }
0834   // TECminus
0835   else if (beam.isTecInternal(-1)) {
0836     if (gHitZprime[hit] > gBeamSplitterZprime) {
0837       trackPhi = offset + slope * gHitZprime[hit];
0838       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
0839     } else {
0840       trackPhi = offset + slope * gHitZprime[hit] + bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime) / gBeamR;
0841       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0) +
0842                     bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime) / gBeamR;
0843     }
0844   }
0845   // ATs
0846   else {
0847     // TECminus
0848     if (gHitZprime[hit] < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
0849       trackPhi = offset + slope * gHitZprime[hit];
0850       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
0851     }
0852     // BarrelMinus
0853     else if (-gBeamSplitterZprime - 2.0 * gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < -gBeamZ0) {
0854       if (!gIsInnerBarrel) {
0855         trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtMinusParam + atThetaSplitParam);
0856       } else {
0857         trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtMinusParam - atThetaSplitParam);
0858       }
0859       trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit - gHitsAtTecMinus]));
0860     }
0861     // BarrelPlus
0862     else if (-gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < gBeamSplitterZprime) {
0863       if (!gIsInnerBarrel) {
0864         trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtPlusParam - atThetaSplitParam);
0865       } else {
0866         trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtPlusParam + atThetaSplitParam);
0867       }
0868       trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit - gHitsAtTecMinus]));
0869     }
0870     // TECplus
0871     else {
0872       trackPhi = offset + slope * gHitZprime[hit] - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime) / gBeamR;
0873       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0) -
0874                     bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime) / gBeamR;
0875     }
0876   }
0877 }
0878 
0879 // -------------- calculate global track points, hit residuals, chi2 ----------------------------------
0880 void TkLasBeamFitter::globalTrackPoint(TkFittedLasBeam &beam,
0881                                        unsigned int &hit,
0882                                        unsigned int &hitsAtTecPlus,
0883                                        double &trackPhi,
0884                                        double &trackPhiRef,
0885                                        vector<GlobalPoint> &globHit,
0886                                        vector<GlobalPoint> &globPtrack,
0887                                        GlobalPoint &globPref,
0888                                        vector<double> &hitPhiError) {
0889   // TECs
0890   if (beam.isTecInternal(0)) {
0891     globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
0892     globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
0893   }
0894   // ATs
0895   else {
0896     // TECminus
0897     if (hit < gHitsAtTecMinus) {  // gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0
0898       globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
0899       globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
0900     }
0901     // TECplus
0902     else if (hit > gHitZprime.size() - hitsAtTecPlus - 1) {  // gHitZprime[hit] > gBeamSplitterZprime
0903       globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
0904       globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
0905     }
0906     // Barrel
0907     else {
0908       globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(globHit[hit].perp(), trackPhi, globHit[hit].z())));
0909       globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z()));
0910     }
0911   }
0912 }
0913 
0914 // ----------- create TrajectoryStateOnSurface for each track hit ----------------------------------------------
0915 void TkLasBeamFitter::buildTrajectory(TkFittedLasBeam &beam,
0916                                       unsigned int &hit,
0917                                       vector<const GeomDetUnit *> &gd,
0918                                       vector<GlobalPoint> &globPtrack,
0919                                       vector<TrajectoryStateOnSurface> &tsosLas,
0920                                       GlobalPoint &globPref) {
0921   const MagneticField *magneticField = fieldHandle.product();
0922   GlobalVector trajectoryState;
0923 
0924   // TECplus
0925   if (beam.isTecInternal(1)) {
0926     trajectoryState = GlobalVector(globPref - globPtrack[hit]);
0927   }
0928   // TECminus
0929   else if (beam.isTecInternal(-1)) {
0930     trajectoryState = GlobalVector(globPtrack[hit] - globPref);
0931   }
0932   // ATs
0933   else {
0934     // TECminus
0935     if (gHitZprime[hit] < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
0936       trajectoryState = GlobalVector(globPtrack[hit] - globPref);
0937     }
0938     // TECplus
0939     else if (gHitZprime[hit] > gBeamSplitterZprime) {
0940       trajectoryState = GlobalVector(globPref - globPtrack[hit]);
0941     }
0942     // Barrel
0943     else {
0944       trajectoryState = GlobalVector(globPtrack[hit] - globPref);
0945     }
0946   }
0947   //   cout << "trajectory: " << trajectoryState << endl;
0948   const FreeTrajectoryState ftsLas = FreeTrajectoryState(globPtrack[hit], trajectoryState, 0, magneticField);
0949   tsosLas.push_back(TrajectoryStateOnSurface(ftsLas, gd[hit]->surface(), SurfaceSideDefinition::beforeSurface));
0950 }
0951 
0952 //---------------------- set beam parameters for fittedBeams ---------------------------------
0953 bool TkLasBeamFitter::fitBeam(TkFittedLasBeam &beam,
0954                               AlgebraicSymMatrix &covMatrix,
0955                               unsigned int &hitsAtTecPlus,
0956                               unsigned int &nFitParams,
0957                               double &offset,
0958                               double &slope,
0959                               vector<GlobalPoint> &globPtrack,
0960                               double &bsAngleParam,
0961                               double &chi2) {
0962   // set beam parameters for beam output
0963   unsigned int paramType(0);
0964   if (!fitBeamSplitters_)
0965     paramType = 1;
0966   if (beam.isAlignmentTube() && hitsAtTecPlus == 0)
0967     paramType = 0;
0968   //   const unsigned int nPedeParams = nFitParams + paramType;
0969 
0970   // test without BS params
0971   const unsigned int nPedeParams(nFitParams);
0972   //   cout << "number of Pede parameters: " << nPedeParams << endl;
0973 
0974   std::vector<TkFittedLasBeam::Scalar> params(nPedeParams);
0975   params[0] = offset;
0976   params[1] = slope;
0977   // no BS parameter for AT beams without TECplus hits
0978   //   if(beam.isTecInternal() || hitsAtTecPlus > 0) params[2] = bsAngleParam;
0979 
0980   AlgebraicMatrix derivatives(gHitZprime.size(), nPedeParams);
0981   // fill derivatives matrix with local track derivatives
0982   for (unsigned int hit = 0; hit < gHitZprime.size(); ++hit) {
0983     // d(delta phi)/d(offset) is identical for every hit
0984     derivatives[hit][0] = 1.0;
0985 
0986     // d(delta phi)/d(slope) and d(delta phi)/d(bsAngleParam) depend on parametrizations
0987     // TECplus
0988     if (beam.isTecInternal(1)) {
0989       derivatives[hit][1] = globPtrack[hit].z();
0990       //       if(gHitZprime[hit] < gBeamSplitterZprime){
0991       //    derivatives[hit][2] = 0.0;
0992       //       }
0993       //       else{
0994       //    derivatives[hit][2] = - (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
0995       //       }
0996     }
0997     // TECminus
0998     else if (beam.isTecInternal(-1)) {
0999       derivatives[hit][1] = globPtrack[hit].z();
1000       //       if(gHitZprime[hit] > gBeamSplitterZprime){
1001       //    derivatives[hit][2] = 0.0;
1002       //       }
1003       //       else{
1004       //    derivatives[hit][2] = (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
1005       //       }
1006     }
1007     // ATs
1008     else {
1009       // TECminus
1010       if (gHitZprime[hit] < -gBeamSplitterZprime - 2.0 * gBeamZ0) {
1011         derivatives[hit][1] = globPtrack[hit].z();
1012         //  if(hitsAtTecPlus > 0){
1013         //    derivatives[hit][2] = 0.0;
1014         //  }
1015       }
1016       // TECplus
1017       else if (gHitZprime[hit] > gBeamSplitterZprime) {
1018         derivatives[hit][1] = globPtrack[hit].z();
1019         //  if(hitsAtTecPlus > 0){
1020         //    derivatives[hit][2] = - (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
1021         //  }
1022       }
1023       // Barrel
1024       else {
1025         derivatives[hit][1] = globPtrack[hit].z() - gBarrelModuleOffset[hit - gHitsAtTecMinus];
1026         //  if(hitsAtTecPlus > 0){
1027         //    derivatives[hit][2] = 0.0;
1028         //  }
1029       }
1030     }
1031   }
1032 
1033   unsigned int firstFixedParam(covMatrix.num_col());  // FIXME! --> no, is fine!!!
1034                                                       //   unsigned int firstFixedParam = nPedeParams - 1;
1035   //   if(beam.isAlignmentTube() && hitsAtTecPlus == 0) firstFixedParam = nPedeParams;
1036   //   cout << "first fixed parameter: " << firstFixedParam << endl;
1037   // set fit results
1038   beam.setParameters(paramType, params, covMatrix, derivatives, firstFixedParam, chi2);
1039 
1040   return true;  // return false in case of problems
1041 }
1042 
1043 //---------------------------------------------------------------------------------------
1044 //define this as a plug-in
1045 DEFINE_FWK_MODULE(TkLasBeamFitter);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!