CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Alignment/LaserAlignment/plugins/TkLasBeamFitter.cc

Go to the documentation of this file.
00001 
00016 // system include files
00017 #include <memory>
00018 
00019 // framework include files
00020 #include "FWCore/Framework/interface/Frameworkfwd.h"
00021 #include "FWCore/Framework/interface/EDProducer.h"
00022 #include "FWCore/Framework/interface/ESHandle.h"
00023 #include "FWCore/Framework/interface/Event.h"
00024 #include "FWCore/Framework/interface/EventSetup.h"
00025 #include "FWCore/Framework/interface/MakerMacros.h"
00026 #include "FWCore/Framework/interface/Run.h"
00027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00028 #include "FWCore/Utilities/interface/InputTag.h"
00029 #include "FWCore/ServiceRegistry/interface/Service.h"
00030 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00031 
00032 #include <Geometry/CommonDetUnit/interface/GeomDetUnit.h>
00033 #include <Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h>
00034 #include <Geometry/Records/interface/TrackerDigiGeometryRecord.h>
00035 #include <MagneticField/Engine/interface/MagneticField.h>
00036 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00037 
00038 // data formats
00039 // for edm::InRun
00040 #include "DataFormats/Provenance/interface/BranchType.h"
00041 
00042 // laser data formats
00043 #include "DataFormats/Alignment/interface/TkLasBeam.h"
00044 #include "DataFormats/Alignment/interface/TkFittedLasBeam.h"
00045 
00046 // further includes
00047 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00048 #include "Alignment/LaserAlignment/interface/TsosVectorCollection.h"
00049 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00050 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00051 
00052 #include <iostream>
00053 #include "TMinuit.h"
00054 #include "TGraphErrors.h"
00055 #include "TF1.h"
00056 #include "TH1.h"
00057 #include "TH2.h"
00058 
00059 using namespace edm;
00060 using namespace std;
00061 
00062 
00063 //
00064 // class declaration
00065 //
00066 
00067 class TkLasBeamFitter : public edm::EDProducer {
00068 public:
00069   explicit TkLasBeamFitter(const edm::ParameterSet &config);
00070   ~TkLasBeamFitter();
00071   
00072   //virtual void beginJob(const edm::EventSetup& /*access deprecated*/) {}
00073   virtual void produce(edm::Event &event, const edm::EventSetup &setup);
00074   // virtual void beginRun(edm::Run &run, const edm::EventSetup &setup);
00075   virtual void endRun(edm::Run &run, const edm::EventSetup &setup);
00076   //virtual void endJob() {}
00077 
00078 private:
00081   void getLasBeams(TkFittedLasBeam &beam,vector<TrajectoryStateOnSurface> &tsosLas);
00082   void getLasHits(TkFittedLasBeam &beam, const SiStripLaserRecHit2D &hit, 
00083                   vector<const GeomDetUnit*> &gd, vector<GlobalPoint> &globHit,
00084                   unsigned int &hitsAtTecPlus); 
00085 //   void fillVectors(TkFittedLasBeam &beam);
00086 
00087   // need static functions to be used in fitter(..);
00088   // all parameters used therein have to be static, as well (see below)
00089   static double tecPlusFunction(double *x, double *par);
00090   static double tecMinusFunction(double *x, double *par);
00091   static double atFunction(double *x, double *par);
00092 
00093   void fitter(TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatrix, 
00094               unsigned int &hitsAtTecPlus, unsigned int &nFitParams,
00095               std::vector<double> &hitPhi, std::vector<double> &hitPhiError, std::vector<double> &hitZprimeError,
00096               double &zMin, double &zMax, double &bsAngleParam,
00097               double &offset, double &offsetError, double &slope, double &slopeError,
00098               double &phiAtMinusParam, double &phiAtPlusParam,
00099               double &atThetaSplitParam);
00100 
00101   void trackPhi(TkFittedLasBeam &beam, unsigned int &hit,
00102                 double &trackPhi, double &trackPhiRef,
00103                 double &offset, double &slope, double &bsAngleParam, 
00104                 double &phiAtMinusParam, double &phiAtPlusParam,
00105                 double &atThetaSplitParam, std::vector<GlobalPoint> &globHit);
00106 
00107   void globalTrackPoint(TkFittedLasBeam &beam, 
00108                         unsigned int &hit, unsigned int &hitsAtTecPlus,
00109                         double &trackPhi, double &trackPhiRef, 
00110                         std::vector<GlobalPoint> &globHit, std::vector<GlobalPoint> &globPtrack,
00111                         GlobalPoint &globPref, std::vector<double> &hitPhiError);
00112 
00113   void buildTrajectory(TkFittedLasBeam &beam, unsigned int &hit,
00114                        vector<const GeomDetUnit*> &gd, std::vector<GlobalPoint> &globPtrack,
00115                        vector<TrajectoryStateOnSurface> &tsosLas, GlobalPoint &globPref);
00116 
00117   bool fitBeam(TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatrix, 
00118                unsigned int &hitsAtTecPlus, unsigned int &nFitParams,
00119                double &offset, double &slope, vector<GlobalPoint> &globPtrack, 
00120                double &bsAngleParam, double &chi2);
00121 
00122   // ----------member data ---------------------------
00123   const edm::InputTag src_;
00124   bool fitBeamSplitters_;
00125   unsigned int nAtParameters_;
00126 
00127   edm::Service<TFileService> fs;
00128 
00129   // static parameters used in static parametrization functions
00130   static vector<double> gHitZprime;
00131   static vector<double> gBarrelModuleRadius;
00132   static vector<double> gBarrelModuleOffset;
00133   static float gTIBparam;
00134   static float gTOBparam;
00135   static double gBeamR;
00136   static double gBeamZ0;
00137   static double gBeamSplitterZprime;
00138   static unsigned int gHitsAtTecMinus;
00139   static double gBSparam;
00140   static bool gFitBeamSplitters;
00141   static bool gIsInnerBarrel;
00142 
00143   // histograms
00144   TH1F *h_bsAngle, *h_hitX, *h_hitXTecPlus, *h_hitXTecMinus,
00145     *h_hitXAt, *h_chi2, *h_chi2ndof, *h_pull, *h_res, 
00146     *h_resTecPlus, *h_resTecMinus, *h_resAt;
00147   TH2F *h_bsAngleVsBeam, *h_hitXvsZTecPlus, *h_hitXvsZTecMinus,
00148     *h_hitXvsZAt, *h_resVsZTecPlus, *h_resVsZTecMinus, *h_resVsZAt,
00149     *h_resVsHitTecPlus, *h_resVsHitTecMinus, *h_resVsHitAt;
00150 }; 
00151 
00152 //
00153 // constants, enums and typedefs
00154 //
00155 
00156 
00157 //
00158 // static data member definitions
00159 //
00160 
00161 // static parameters used in parametrization functions
00162 vector<double> TkLasBeamFitter::gHitZprime;
00163 vector<double> TkLasBeamFitter::gBarrelModuleRadius;
00164 vector<double> TkLasBeamFitter::gBarrelModuleOffset;
00165 float TkLasBeamFitter::gTIBparam = 0.097614; // = abs(r_offset/r_module) (nominal!)
00166 float TkLasBeamFitter::gTOBparam = 0.034949; // = abs(r_offset/r_module) (nominal!)
00167 double TkLasBeamFitter::gBeamR = 0.0;
00168 double TkLasBeamFitter::gBeamZ0 = 0.0; 
00169 double TkLasBeamFitter::gBeamSplitterZprime = 0.0;
00170 unsigned int TkLasBeamFitter::gHitsAtTecMinus = 0;
00171 double TkLasBeamFitter::gBSparam = 0.0;
00172 bool TkLasBeamFitter::gFitBeamSplitters = 0;
00173 bool TkLasBeamFitter::gIsInnerBarrel = 0;
00174 
00175 // handles
00176 Handle<TkLasBeamCollection> laserBeams;
00177 ESHandle<MagneticField> fieldHandle;
00178 ESHandle<TrackerGeometry> geometry;
00179 
00180 //
00181 // constructors and destructor
00182 //
00183 TkLasBeamFitter::TkLasBeamFitter(const edm::ParameterSet &iConfig) :
00184   src_(iConfig.getParameter<edm::InputTag>("src")),
00185   fitBeamSplitters_(iConfig.getParameter<bool>("fitBeamSplitters")),
00186   nAtParameters_(iConfig.getParameter<unsigned int>("numberOfFittedAtParameters"))
00187 {
00188   // declare the products to produce
00189   this->produces<TkFittedLasBeamCollection, edm::InRun>();
00190   this->produces<TsosVectorCollection, edm::InRun>();
00191   
00192   //now do what ever other initialization is needed
00193 }
00194 
00195 //---------------------------------------------------------------------------------------
00196 TkLasBeamFitter::~TkLasBeamFitter()
00197 {
00198  
00199    // do anything here that needs to be done at desctruction time
00200    // (e.g. close files, deallocate resources etc.)
00201 
00202 }
00203 
00204 
00205 //
00206 // member functions
00207 //
00208 
00209 //---------------------------------------------------------------------------------------
00210 // ------------ method called to produce the data  ------------
00211 void TkLasBeamFitter::produce(edm::Event &iEvent, const edm::EventSetup &setup)
00212 {
00213   // Nothing per event!
00214 }
00215 
00216 //---------------------------------------------------------------------------------------
00217 // ------------ method called at end of each run  ---------------------------------------
00218 void TkLasBeamFitter::endRun(edm::Run &run, const edm::EventSetup &setup)
00219 {
00220 // }
00221 // // FIXME!
00222 // // Indeed, that should be in endRun(..) - as soon as AlignmentProducer can call
00223 // // the algorithm's endRun correctly!
00224 //
00225 //
00226 // void TkLasBeamFitter::beginRun(edm::Run &run, const edm::EventSetup &setup)
00227 // {
00228 
00229   // book histograms
00230   h_hitX = fs->make<TH1F>("hitX","local x of LAS hits;local x [cm];N",100,-0.5,0.5);
00231   h_hitXTecPlus = fs->make<TH1F>("hitXTecPlus","local x of LAS hits in TECplus;local x [cm];N",100,-0.5,0.5);
00232   h_hitXTecMinus = fs->make<TH1F>("hitXTecMinus","local x of LAS hits in TECminus;local x [cm];N",100,-0.5,0.5);
00233   h_hitXAt = fs->make<TH1F>("hitXAt","local x of LAS hits in ATs;local x [cm];N",100,-2.5,2.5);
00234   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);
00235   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);
00236   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);
00237   h_chi2 = fs->make<TH1F>("chi2","#chi^{2};#chi^{2};N",100,0,2000);
00238   h_chi2ndof = fs->make<TH1F>("chi2ndof","#chi^{2} per degree of freedom;#chi^{2}/N_{dof};N",100,0,300);
00239   h_pull = fs->make<TH1F>("pull","pulls of #phi residuals;pull;N",50,-10,10);
00240   h_res = fs->make<TH1F>("res","#phi residuals;#phi_{track} - #phi_{hit} [rad];N",60,-0.0015,0.0015);
00241   h_resTecPlus = fs->make<TH1F>("resTecPlus","#phi residuals in TECplus;#phi_{track} - #phi_{hit} [rad];N",30,-0.0015,0.0015);
00242   h_resTecMinus = fs->make<TH1F>("resTecMinus","#phi residuals in TECminus;#phi_{track} - #phi_{hit} [rad];N",60,-0.0015,0.0015);
00243   h_resAt = fs->make<TH1F>("resAt","#phi residuals in ATs;#phi_{track} - #phi_{hit} [rad];N",30,-0.0015,0.0015);
00244   h_resVsZTecPlus = fs->make<TH2F>("resVsZTecPlus","phi residuals vs. z in TECplus;z [cm];#phi_{track} - #phi_{hit} [rad]",
00245                          80,120,280,100,-0.0015,0.0015);
00246   h_resVsZTecMinus = fs->make<TH2F>("resVsZTecMinus","phi residuals vs. z in TECminus;z [cm];#phi_{track} - #phi_{hit} [rad]",
00247                           80,-280,-120,100,-0.0015,0.0015);
00248   h_resVsZAt = fs->make<TH2F>("resVsZAt","#phi residuals vs. z in ATs;N;#phi_{track} - #phi_{hit} [rad]",
00249                     200,-200,200,100,-0.0015,0.0015);
00250   h_resVsHitTecPlus = fs->make<TH2F>("resVsHitTecPlus","#phi residuals vs. hits in TECplus;hit no.;#phi_{track} - #phi_{hit} [rad]",
00251                            144,0,144,100,-0.0015,0.0015);
00252   h_resVsHitTecMinus = fs->make<TH2F>("resVsHitTecMinus","#phi residuals vs. hits in TECminus;hit no.;#phi_{track} - #phi_{hit} [rad]",
00253                             144,0,144,100,-0.0015,0.0015);
00254   h_resVsHitAt = fs->make<TH2F>("resVsHitAt","#phi residuals vs. hits in ATs;hit no.;#phi_{track} - #phi_{hit} [rad]",
00255                       176,0,176,100,-0.0015,0.0015);
00256   h_bsAngle = fs->make<TH1F>("bsAngle","fitted beam splitter angle;BS angle [rad];N",40,-0.004,0.004);
00257   h_bsAngleVsBeam = fs->make<TH2F>("bsAngleVsBeam","fitted beam splitter angle per beam;Beam no.;BS angle [rad]",
00258                          40,0,300,100,-0.004,0.004);
00259 
00260   // Create output collections - they are parallel.
00261   // (edm::Ref etc. and thus edm::AssociationVector are not supported for edm::Run...)
00262   std::auto_ptr<TkFittedLasBeamCollection> fittedBeams(new TkFittedLasBeamCollection);
00263   // One std::vector<TSOS> for each TkFittedLasBeam:
00264   std::auto_ptr<TsosVectorCollection> tsosesVec(new TsosVectorCollection);
00265 
00266   // get TkLasBeams, Tracker geometry, magnetic field
00267   run.getByLabel( "LaserAlignment", "tkLaserBeams", laserBeams );
00268   setup.get<TrackerDigiGeometryRecord>().get(geometry);
00269   setup.get<IdealMagneticFieldRecord>().get(fieldHandle);
00270 
00271   // hack for fixed BSparams (ugly!)
00272 //   double bsParams[34] = {-0.000266,-0.000956,-0.001205,-0.000018,-0.000759,0.002554,
00273 //                       0.000465,0.000975,0.001006,0.002027,-0.001263,-0.000763,
00274 //                       -0.001702,0.000906,-0.002120,0.001594,0.000661,-0.000457,
00275 //                       -0.000447,0.000347,-0.002266,-0.000446,0.000659,0.000018,
00276 //                       -0.001630,-0.000324,
00277 //                       // ATs
00278 //                       -999.,-0.001709,-0.002091,-999.,
00279 //                       -0.001640,-999.,-0.002444,-0.002345};
00280 
00281   double bsParams[40] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
00282                          0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
00283 
00284   // beam counter
00285   unsigned int beamNo(0);
00286   // fit BS? If false, values from bsParams are taken
00287   gFitBeamSplitters = fitBeamSplitters_;
00288   if(fitBeamSplitters_) cout << "Fitting BS!" << endl;
00289   else cout << "BS fixed, not fitted!" << endl;
00290 
00291   // loop over LAS beams
00292   for(TkLasBeamCollection::const_iterator iBeam = laserBeams->begin(), iEnd = laserBeams->end();
00293        iBeam != iEnd; ++iBeam){
00294 
00295     TkFittedLasBeam beam(*iBeam);
00296     vector<TrajectoryStateOnSurface> tsosLas;
00297 
00298     // set BS param for fit
00299     gBSparam = bsParams[beamNo];
00300 
00301     // call main function; all other functions are called inside getLasBeams(..)
00302     this->getLasBeams(beam, tsosLas);
00303     
00304     // fill output products
00305     fittedBeams->push_back(beam);
00306     tsosesVec->push_back(tsosLas);
00307 
00308 //     if(!this->fitBeam(fittedBeams->back(), tsosesVec->back())){
00309 //       edm::LogError("BadFit") 
00310 //       << "Problems fitting TkLasBeam, id " << fittedBeams->back().getBeamId() << ".";
00311 //       fittedBeams->pop_back(); // remove last entry added just before
00312 //       tsosesVec->pop_back();   // dito
00313 //     }
00314 
00315     beamNo++;
00316   }
00317   
00318   // finally, put fitted beams and TSOS vectors into run
00319   run.put(fittedBeams);
00320   run.put(tsosesVec);
00321 }
00322 
00323 // methods for las data processing
00324 
00325 // -------------- loop over beams, call functions ----------------------------
00326 void TkLasBeamFitter::getLasBeams(TkFittedLasBeam &beam, vector<TrajectoryStateOnSurface> &tsosLas)
00327 {
00328   cout << "---------------------------------------" << endl;
00329   cout << "beam id: " << beam.getBeamId() // << " isTec: " << (beam.isTecInternal() ? "Y" : "N") 
00330        << " isTec+: " << (beam.isTecInternal(1) ? "Y" : "N") << " isTec-: " << (beam.isTecInternal(-1) ? "Y" : "N")
00331        << " isAt: " << (beam.isAlignmentTube() ? "Y" : "N") << " isR6: " << (beam.isRing6() ? "Y" : "N")
00332        << endl;
00333  
00334   // reset static variables 
00335   gHitsAtTecMinus = 0;
00336   gHitZprime.clear();
00337   gBarrelModuleRadius.clear();
00338   gBarrelModuleOffset.clear();
00339 
00340   // set right beam radius
00341   gBeamR = beam.isRing6() ? 84.0 : 56.4;
00342   
00343   vector<const GeomDetUnit*> gd;
00344   vector<GlobalPoint> globHit;
00345   unsigned int hitsAtTecPlus(0);
00346   double sumZ(0.);
00347 
00348   // loop over hits
00349   for( TkLasBeam::const_iterator iHit = beam.begin(); iHit < beam.end(); ++iHit ){
00350     // iHit is a SiStripLaserRecHit2D
00351 
00352     const SiStripLaserRecHit2D hit(*iHit);
00353 
00354     this->getLasHits(beam, hit, gd, globHit, hitsAtTecPlus);
00355     sumZ += globHit.back().z();
00356 
00357     // fill histos
00358     h_hitX->Fill(hit.localPosition().x());
00359     // TECplus
00360     if(beam.isTecInternal(1)){
00361       h_hitXTecPlus->Fill(hit.localPosition().x());
00362       h_hitXvsZTecPlus->Fill(globHit.back().z(),hit.localPosition().x());
00363     }
00364     // TECminus
00365     else if(beam.isTecInternal(-1)){
00366       h_hitXTecMinus->Fill(hit.localPosition().x());
00367       h_hitXvsZTecMinus->Fill(globHit.back().z(),hit.localPosition().x());
00368     }
00369     // ATs
00370     else{
00371       h_hitXAt->Fill(hit.localPosition().x());
00372       h_hitXvsZAt->Fill(globHit.back().z(),hit.localPosition().x());
00373     }  
00374   }
00375 
00376   gBeamZ0 = sumZ / globHit.size();
00377   double zMin(0.), zMax(0.);
00378   // TECplus
00379   if(beam.isTecInternal(1)){
00380     gBeamSplitterZprime = 205.75 - gBeamZ0;
00381     zMin = 120.0 - gBeamZ0;
00382     zMax = 280.0 - gBeamZ0;
00383   }
00384   // TECminus
00385   else if(beam.isTecInternal(-1)){
00386     gBeamSplitterZprime = -205.75 - gBeamZ0;
00387     zMin = -280.0 - gBeamZ0;
00388     zMax = -120.0 - gBeamZ0;
00389   }
00390   // AT
00391   else{
00392     gBeamSplitterZprime = 112.3 - gBeamZ0;
00393     zMin = -200.0 - gBeamZ0;
00394     zMax = 200.0 - gBeamZ0;
00395   }
00396 
00397   // fill vectors for fitted quantities
00398   vector<double> hitPhi, hitPhiError, hitZprimeError;
00399 
00400   for(unsigned int hit = 0; hit < globHit.size(); ++hit){
00401     hitPhi.push_back(static_cast<double>(globHit[hit].phi()));
00402     // localPositionError[hit] or assume 0.003, 0.006
00403     hitPhiError.push_back( 0.003 / globHit[hit].perp());
00404     // no errors on z, fill with zeros
00405     hitZprimeError.push_back(0.0);
00406     // barrel-specific values
00407     if(beam.isAlignmentTube() && abs(globHit[hit].z()) < 112.3){
00408       gBarrelModuleRadius.push_back(globHit[hit].perp());
00409       gBarrelModuleOffset.push_back(gBarrelModuleRadius.back() - gBeamR);
00410       // TIB/TOB flag
00411       if(gBarrelModuleOffset.back() < 0.0){
00412         gIsInnerBarrel = 1;
00413       }
00414       else{
00415         gIsInnerBarrel = 0;
00416       }
00417       gHitZprime.push_back(globHit[hit].z() - gBeamZ0 - abs(gBarrelModuleOffset.back())); 
00418     }
00419     // non-barrel z'
00420     else{
00421       gHitZprime.push_back(globHit[hit].z() - gBeamZ0);
00422     }
00423   }
00424 
00425   // number of fit parameters, 3 for TECs (always!); 3, 5, or 6 for ATs
00426   unsigned int tecParams(3), atParams(0);
00427   if(nAtParameters_ == 3) atParams = 3;
00428   else if(nAtParameters_ == 5) atParams = 5;
00429   else atParams = 6;                            // <-- default value, recommended
00430   unsigned int nFitParams(0);
00431   if(!fitBeamSplitters_ || 
00432      (hitsAtTecPlus == 0 && beam.isAlignmentTube() ) ){
00433     tecParams = tecParams - 1;
00434     atParams = atParams - 1;
00435   }
00436   if(beam.isTecInternal()){
00437     nFitParams = tecParams;
00438   }
00439   else{
00440     nFitParams = atParams;
00441   }
00442   
00443   // fit parameter definitions
00444   double offset(0.), offsetError(0.), slope(0.), slopeError(0.),
00445     bsAngleParam(0.), phiAtMinusParam(0.), phiAtPlusParam(0.),
00446     atThetaSplitParam(0.);
00447   AlgebraicSymMatrix covMatrix;
00448   if(!fitBeamSplitters_ || (beam.isAlignmentTube() && hitsAtTecPlus == 0)){
00449     covMatrix = AlgebraicSymMatrix(nFitParams, 1);
00450   }
00451   else{
00452     covMatrix = AlgebraicSymMatrix(nFitParams - 1, 1);
00453   }
00454 
00455   this->fitter(beam, covMatrix, 
00456                hitsAtTecPlus, nFitParams, 
00457                hitPhi, hitPhiError, hitZprimeError, 
00458                zMin, zMax, bsAngleParam,
00459                offset, offsetError, slope, slopeError,
00460                phiAtMinusParam, phiAtPlusParam,
00461                atThetaSplitParam); 
00462   
00463   vector<GlobalPoint> globPtrack;
00464   GlobalPoint globPref;
00465   double chi2(0.);
00466 
00467   for(unsigned int hit = 0; hit < gHitZprime.size(); ++hit){
00468     
00469     // additional phi value (trackPhiRef) for trajectory calculation
00470     double trackPhi(0.), trackPhiRef(0.); 
00471     
00472     this->trackPhi(beam, hit, trackPhi, trackPhiRef,
00473                    offset, slope, bsAngleParam, 
00474                    phiAtMinusParam, phiAtPlusParam,
00475                    atThetaSplitParam, globHit);
00476 
00477     cout << "track phi = " << trackPhi 
00478          << ", hit phi = " << hitPhi[hit] 
00479          << ", zPrime = " << gHitZprime[hit] 
00480          << " r = " << globHit[hit].perp() << endl;
00481 
00482     this->globalTrackPoint(beam, hit, hitsAtTecPlus, 
00483                            trackPhi, trackPhiRef, 
00484                            globHit, globPtrack, globPref, 
00485                            hitPhiError);
00486 
00487     // calculate residuals = pred - hit (in global phi)
00488     const double phiResidual = globPtrack[hit].phi() - globHit[hit].phi();
00489     // pull calculation (FIX!!!)
00490     const double phiResidualPull = phiResidual / hitPhiError[hit];
00491     //       sqrt(hitPhiError[hit]*hitPhiError[hit] + 
00492     //     (offsetError*offsetError + globPtrack[hit].z()*globPtrack[hit].z() * slopeError*slopeError));
00493 
00494     // calculate chi2
00495     chi2 += phiResidual*phiResidual / (hitPhiError[hit]*hitPhiError[hit]);
00496     
00497     // fill histos
00498     h_res->Fill(phiResidual);
00499     // TECplus
00500     if(beam.isTecInternal(1)){
00501       h_pull->Fill(phiResidualPull);
00502       h_resTecPlus->Fill(phiResidual);
00503       h_resVsZTecPlus->Fill(globPtrack[hit].z(), phiResidual);
00504       // Ring 6
00505       if(beam.isRing6()){
00506         h_resVsHitTecPlus->Fill(hit+(beam.getBeamId()-1)/10*9+72, phiResidual);
00507       }
00508       // Ring 4
00509       else{
00510         h_resVsHitTecPlus->Fill(hit+beam.getBeamId()/10*9, phiResidual);
00511       }
00512     }
00513     // TECminus
00514     else if(beam.isTecInternal(-1)){
00515       h_pull->Fill(phiResidualPull);
00516       h_resTecMinus->Fill(phiResidual);
00517       h_resVsZTecMinus->Fill(globPtrack[hit].z(), phiResidual);
00518       // Ring 6
00519       if(beam.isRing6()){
00520         h_resVsHitTecMinus->Fill(hit+(beam.getBeamId()-101)/10*9+72, phiResidual);
00521       }
00522       // Ring 4
00523       else{
00524         h_resVsHitTecMinus->Fill(hit+(beam.getBeamId()-100)/10*9, phiResidual);
00525       }
00526     }
00527     // ATs
00528     else{
00529       h_pull->Fill(phiResidualPull);
00530       h_resAt->Fill(phiResidual);
00531       h_resVsZAt->Fill(globPtrack[hit].z(), phiResidual);
00532       h_resVsHitAt->Fill(hit+(beam.getBeamId()-200)/10*22, phiResidual);
00533     }
00534 
00535     this->buildTrajectory(beam, hit, gd, globPtrack, tsosLas, globPref);
00536   }
00537 
00538   cout << "chi^2 = " << chi2 << ", chi^2/ndof = " << chi2/(gHitZprime.size() - nFitParams) << endl;
00539   this->fitBeam(beam, covMatrix, hitsAtTecPlus, nFitParams,
00540                 offset, slope, globPtrack, bsAngleParam, chi2);
00541   
00542   cout << "bsAngleParam = " << bsAngleParam << endl;
00543 
00544   // fill histos
00545   // include slope, offset, covariance plots here
00546   h_chi2->Fill(chi2);
00547   h_chi2ndof->Fill(chi2/(gHitZprime.size() - nFitParams));
00548   if(bsAngleParam != 0.0){
00549     h_bsAngle->Fill(2.0*atan(0.5*bsAngleParam));
00550     h_bsAngleVsBeam->Fill(beam.getBeamId(), 2.0*atan(0.5*bsAngleParam));
00551   }
00552 }
00553 
00554 
00555 // --------- get hits, convert to global coordinates ---------------------------
00556 void TkLasBeamFitter::getLasHits(TkFittedLasBeam &beam, const SiStripLaserRecHit2D &hit, 
00557                                  vector<const GeomDetUnit*> &gd, vector<GlobalPoint> &globHit,
00558                                  unsigned int &hitsAtTecPlus)
00559 { 
00560   // get global position of LAS hits
00561   gd.push_back(geometry->idToDetUnit(hit.getDetId()));
00562   GlobalPoint globPtemp(gd.back()->toGlobal(hit.localPosition()));
00563 
00564   // testing: globPtemp should be right
00565   globHit.push_back(globPtemp);
00566 
00567   if(beam.isAlignmentTube()){
00568     if(abs(globPtemp.z()) > 112.3){
00569       if(globPtemp.z() < 112.3) gHitsAtTecMinus++ ;
00570       else hitsAtTecPlus++ ;
00571     }
00572   }
00573 }
00574 
00575 
00576 // ------------ parametrization functions for las beam fits ------------
00577 double TkLasBeamFitter::tecPlusFunction(double *x, double *par)
00578 {
00579   double z = x[0]; // 'primed'? -> yes!!!
00580 
00581   if(z < gBeamSplitterZprime){
00582     return par[0] + par[1] * z;
00583   } 
00584   else{
00585     if(gFitBeamSplitters){
00586       // par[2] = 2*tan(BeamSplitterAngle/2.0)
00587       return par[0] + par[1] * z - par[2] * (z - gBeamSplitterZprime)/gBeamR; 
00588     }
00589     else{
00590       return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime)/gBeamR; 
00591     }
00592   }
00593 }
00594 
00595 
00596 double TkLasBeamFitter::tecMinusFunction(double *x, double *par)
00597 {
00598   double z = x[0]; // 'primed'? -> yes!!!
00599 
00600   if(z > gBeamSplitterZprime){
00601     return par[0] + par[1] * z;
00602   } 
00603   else{
00604     if(gFitBeamSplitters){
00605       // par[2] = 2*tan(BeamSplitterAngle/2.0)
00606       return par[0] + par[1] * z + par[2] * (z - gBeamSplitterZprime)/gBeamR; 
00607     }
00608     else{
00609       return par[0] + par[1] * z + gBSparam * (z - gBeamSplitterZprime)/gBeamR; 
00610     }
00611   }
00612 }
00613  
00614 double TkLasBeamFitter::atFunction(double *x, double *par)
00615 {
00616   double z = x[0]; // 'primed'? -> yes!!!
00617   // TECminus
00618   if(z < -gBeamSplitterZprime - 2.0*gBeamZ0){
00619     return par[0] + par[1] * z;
00620   }
00621   // BarrelMinus
00622   else if(-gBeamSplitterZprime - 2.0*gBeamZ0 < z && z < -gBeamZ0){
00623     // z value includes module offset from main beam axis
00624     // TOB
00625     if(!gIsInnerBarrel){
00626       return par[0] + par[1] * z + gTOBparam * (par[2] + par[4]);
00627     }
00628     // TIB
00629     else{
00630       return par[0] + par[1] * z - gTIBparam * (par[2] - par[4]);
00631     }
00632   }
00633   // BarrelPlus
00634   else if(-gBeamZ0 < z && z < gBeamSplitterZprime){
00635     // z value includes module offset from main beam axis
00636     // TOB
00637     if(!gIsInnerBarrel){
00638       return par[0] + par[1] * z + gTOBparam * (par[3] - par[4]);
00639     }
00640     // TIB
00641     else{
00642       return par[0] + par[1] * z - gTIBparam * (par[3] + par[4]);
00643     }
00644   }
00645   // TECplus
00646   else{
00647     if(gFitBeamSplitters){
00648       // par[2] = 2*tan(BeamSplitterAngle/2.0)
00649       return par[0] + par[1] * z - par[5] * (z - gBeamSplitterZprime)/gBeamR; // BS par: 5, 4, or 2
00650     }
00651     else{
00652       return par[0] + par[1] * z - gBSparam * (z - gBeamSplitterZprime)/gBeamR;
00653     }
00654   }
00655 }
00656 
00657 
00658 // ------------ perform fit of beams ------------------------------------
00659 void TkLasBeamFitter::fitter(TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatrix, 
00660                              unsigned int &hitsAtTecPlus, unsigned int &nFitParams,
00661                              vector<double> &hitPhi, vector<double> &hitPhiError, vector<double> &hitZprimeError,
00662                              double &zMin, double &zMax, double &bsAngleParam, 
00663                              double &offset, double &offsetError, double &slope, double &slopeError,
00664                              double &phiAtMinusParam, double &phiAtPlusParam, double &atThetaSplitParam)
00665 {
00666   TGraphErrors *lasData = new TGraphErrors(gHitZprime.size(), 
00667                                            &(gHitZprime[0]), &(hitPhi[0]), 
00668                                            &(hitZprimeError[0]), &(hitPhiError[0]));
00669   
00670   // do fit (R = entire range)
00671   if(beam.isTecInternal(1)){
00672     TF1 tecPlus("tecPlus", tecPlusFunction, zMin, zMax, nFitParams );
00673     tecPlus.SetParameter( 1, 0 ); // slope
00674     tecPlus.SetParameter( nFitParams - 1, 0 ); // BS 
00675     lasData->Fit("tecPlus", "R"); // "R", "RV" or "RQ"
00676   }
00677   else if(beam.isTecInternal(-1)){
00678     TF1 tecMinus("tecMinus", tecMinusFunction, zMin, zMax, nFitParams );
00679     tecMinus.SetParameter( 1, 0 ); // slope
00680     tecMinus.SetParameter( nFitParams - 1, 0 ); // BS 
00681     lasData->Fit("tecMinus", "R");
00682   }
00683   else{
00684     TF1 at("at", atFunction, zMin, zMax, nFitParams );
00685     at.SetParameter( 1, 0 ); // slope
00686     at.SetParameter( nFitParams - 1, 0 ); // BS 
00687     lasData->Fit("at","R");
00688   }
00689   
00690   // get values and errors for offset and slope
00691   gMinuit->GetParameter(0, offset, offsetError);
00692   gMinuit->GetParameter(1, slope, slopeError);
00693 
00694   // additional AT parameters
00695   // define param errors that are not used later
00696   double bsAngleParamError(0.), phiAtMinusParamError(0.),
00697     phiAtPlusParamError(0.), atThetaSplitParamError(0.);
00698 
00699   if(beam.isAlignmentTube()){
00700     gMinuit->GetParameter(2, phiAtMinusParam, phiAtMinusParamError);
00701     gMinuit->GetParameter(3, phiAtPlusParam, phiAtPlusParamError);
00702     gMinuit->GetParameter(4, atThetaSplitParam, atThetaSplitParamError);
00703   }
00704   // get Beam Splitter parameters
00705   if(fitBeamSplitters_){
00706     if(beam.isAlignmentTube() && hitsAtTecPlus == 0){
00707       bsAngleParam = gBSparam;
00708     }
00709     else{
00710       gMinuit->GetParameter( nFitParams - 1 , bsAngleParam, bsAngleParamError);
00711     }
00712   }
00713   else{
00714     bsAngleParam = gBSparam;
00715   }
00716 
00717   // fill covariance matrix
00718   vector<double> vec( covMatrix.num_col() * covMatrix.num_col() );
00719   gMinuit->mnemat( &vec[0], covMatrix.num_col() );
00720   for(int col = 0; col < covMatrix.num_col(); col++){
00721     for(int row = 0; row < covMatrix.num_col(); row++){
00722       covMatrix[col][row] = vec[row + covMatrix.num_col()*col];
00723     }
00724   }
00725   // compute correlation between parameters
00726 //   double corr01 = covMatrix[1][0]/(offsetError*slopeError);
00727 
00728   delete lasData;
00729 }
00730 
00731 
00732 // -------------- calculate track phi value ----------------------------------
00733 void TkLasBeamFitter::trackPhi(TkFittedLasBeam &beam, unsigned int &hit, 
00734                                double &trackPhi, double &trackPhiRef,
00735                                double &offset, double &slope, double &bsAngleParam,
00736                                double &phiAtMinusParam, double &phiAtPlusParam,
00737                                double &atThetaSplitParam, vector<GlobalPoint> &globHit)
00738 {
00739   // TECplus
00740   if(beam.isTecInternal(1)){
00741     if(gHitZprime[hit] < gBeamSplitterZprime){
00742       trackPhi = offset + slope * gHitZprime[hit];
00743       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
00744     }
00745     else{
00746       trackPhi = offset + slope * gHitZprime[hit] 
00747         - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime)/gBeamR;
00748       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0)
00749         - bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime)/gBeamR;
00750     }
00751   }
00752   // TECminus
00753   else if(beam.isTecInternal(-1)){
00754     if(gHitZprime[hit] > gBeamSplitterZprime){
00755       trackPhi = offset + slope * gHitZprime[hit];
00756       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
00757     }
00758     else{
00759       trackPhi = offset + slope * gHitZprime[hit] 
00760         + bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime)/gBeamR;
00761       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0)
00762         + bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime)/gBeamR;
00763     }
00764   }
00765   // ATs
00766   else{
00767     // TECminus
00768     if(gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0){
00769       trackPhi = offset + slope * gHitZprime[hit];
00770       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0);
00771     }
00772     // BarrelMinus
00773     else if(-gBeamSplitterZprime - 2.0*gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < -gBeamZ0){
00774       if(!gIsInnerBarrel){
00775         trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtMinusParam + atThetaSplitParam); 
00776       }
00777       else{
00778         trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtMinusParam - atThetaSplitParam);
00779       }
00780       trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit-gHitsAtTecMinus]));
00781     }
00782     // BarrelPlus
00783     else if(-gBeamZ0 < gHitZprime[hit] && gHitZprime[hit] < gBeamSplitterZprime){
00784       if(!gIsInnerBarrel){
00785         trackPhi = offset + slope * gHitZprime[hit] + gTOBparam * (phiAtPlusParam - atThetaSplitParam);
00786       }
00787       else{
00788         trackPhi = offset + slope * gHitZprime[hit] - gTIBparam * (phiAtPlusParam + atThetaSplitParam);
00789       }
00790       trackPhiRef = offset + slope * (gHitZprime[hit] + abs(gBarrelModuleOffset[hit-gHitsAtTecMinus]));
00791     }
00792     // TECplus
00793     else{
00794       trackPhi = offset + slope * gHitZprime[hit] 
00795         - bsAngleParam * (gHitZprime[hit] - gBeamSplitterZprime)/gBeamR;
00796       trackPhiRef = offset + slope * (gHitZprime[hit] + 1.0)
00797         - bsAngleParam * ((gHitZprime[hit] + 1.0) - gBeamSplitterZprime)/gBeamR;
00798     }
00799   }
00800 }
00801 
00802 
00803 // -------------- calculate global track points, hit residuals, chi2 ----------------------------------
00804 void TkLasBeamFitter::globalTrackPoint(TkFittedLasBeam &beam, 
00805                                        unsigned int &hit, unsigned int &hitsAtTecPlus,
00806                                        double &trackPhi, double &trackPhiRef, 
00807                                        vector<GlobalPoint> &globHit, vector<GlobalPoint> &globPtrack,
00808                                        GlobalPoint &globPref, vector<double> &hitPhiError)
00809 {
00810   // TECs
00811   if(beam.isTecInternal(0)){
00812     globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
00813     globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
00814   }
00815   // ATs
00816   else{
00817     // TECminus
00818     if(hit < gHitsAtTecMinus){ // gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0
00819       globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
00820       globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
00821     }
00822     // TECplus
00823     else if(hit > gHitZprime.size() - hitsAtTecPlus - 1){ // gHitZprime[hit] > gBeamSplitterZprime
00824       globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhi, globHit[hit].z())));
00825       globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z() + 1.0));
00826     }
00827     // Barrel
00828     else{
00829       globPtrack.push_back(GlobalPoint(GlobalPoint::Cylindrical(globHit[hit].perp(), trackPhi, globHit[hit].z())));
00830       globPref = GlobalPoint(GlobalPoint::Cylindrical(gBeamR, trackPhiRef, globHit[hit].z()));
00831     }
00832   }
00833 }
00834 
00835 
00836 // ----------- create TrajectoryStateOnSurface for each track hit ----------------------------------------------
00837 void TkLasBeamFitter::buildTrajectory(TkFittedLasBeam &beam, unsigned int &hit, 
00838                                       vector<const GeomDetUnit*> &gd, vector<GlobalPoint> &globPtrack,
00839                                       vector<TrajectoryStateOnSurface> &tsosLas, GlobalPoint &globPref)
00840 {
00841   const MagneticField* magneticField = fieldHandle.product();
00842   GlobalVector trajectoryState;
00843   
00844   // TECplus
00845   if(beam.isTecInternal(1)){
00846     trajectoryState = GlobalVector(globPref-globPtrack[hit]);
00847   }
00848   // TECminus
00849   else if(beam.isTecInternal(-1)){
00850     trajectoryState = GlobalVector(globPtrack[hit]-globPref);
00851   }
00852   // ATs
00853   else{
00854     // TECminus
00855     if(gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0){
00856       trajectoryState = GlobalVector(globPtrack[hit]-globPref);
00857     }
00858     // TECplus
00859     else if(gHitZprime[hit] > gBeamSplitterZprime){
00860       trajectoryState = GlobalVector(globPref-globPtrack[hit]);
00861     }
00862     // Barrel
00863     else{
00864       trajectoryState = GlobalVector(globPtrack[hit]-globPref);
00865     }
00866   }     
00867 //   cout << "trajectory: " << trajectoryState << endl;                 
00868   const FreeTrajectoryState ftsLas = FreeTrajectoryState(globPtrack[hit],trajectoryState,0,magneticField);
00869   tsosLas.push_back(TrajectoryStateOnSurface(ftsLas,gd[hit]->surface(),
00870                                              SurfaceSideDefinition::beforeSurface));
00871 }
00872 
00873 
00874 //---------------------- set beam parameters for fittedBeams ---------------------------------
00875 bool TkLasBeamFitter::fitBeam(TkFittedLasBeam &beam, AlgebraicSymMatrix &covMatrix, 
00876                               unsigned int &hitsAtTecPlus, unsigned int &nFitParams,
00877                               double &offset, double &slope, vector<GlobalPoint> &globPtrack,
00878                               double &bsAngleParam, double &chi2)
00879 {
00880   // set beam parameters for beam output
00881   unsigned int paramType(0);
00882   if(!fitBeamSplitters_) paramType = 1;
00883   if(beam.isAlignmentTube() && hitsAtTecPlus == 0) paramType = 0;
00884 //   const unsigned int nPedeParams = nFitParams + paramType;
00885 
00886   // test without BS params
00887   const unsigned int nPedeParams(nFitParams);
00888 //   cout << "number of Pede parameters: " << nPedeParams << endl;
00889 
00890   std::vector<TkFittedLasBeam::Scalar> params(nPedeParams);
00891   params[0] = offset;
00892   params[1] = slope;
00893   // no BS parameter for AT beams without TECplus hits
00894 //   if(beam.isTecInternal() || hitsAtTecPlus > 0) params[2] = bsAngleParam;
00895 
00896   AlgebraicMatrix derivatives(gHitZprime.size(), nPedeParams);
00897   // fill derivatives matrix with local track derivatives
00898   for(unsigned int hit = 0; hit < gHitZprime.size(); ++hit){
00899     
00900     // d(delta phi)/d(offset) is identical for every hit
00901     derivatives[hit][0] = 1.0;
00902 
00903     // d(delta phi)/d(slope) and d(delta phi)/d(bsAngleParam) depend on parametrizations
00904     // TECplus
00905     if(beam.isTecInternal(1)){
00906       derivatives[hit][1] = globPtrack[hit].z();
00907 //       if(gHitZprime[hit] < gBeamSplitterZprime){
00908 //      derivatives[hit][2] = 0.0;
00909 //       }
00910 //       else{
00911 //      derivatives[hit][2] = - (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
00912 //       }
00913     }
00914     // TECminus
00915     else if(beam.isTecInternal(-1)){
00916       derivatives[hit][1] = globPtrack[hit].z();
00917 //       if(gHitZprime[hit] > gBeamSplitterZprime){
00918 //      derivatives[hit][2] = 0.0;
00919 //       }
00920 //       else{
00921 //      derivatives[hit][2] = (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
00922 //       }
00923     }
00924     // ATs
00925     else{
00926       // TECminus
00927       if(gHitZprime[hit] < -gBeamSplitterZprime - 2.0*gBeamZ0){
00928         derivatives[hit][1] = globPtrack[hit].z();
00929 //      if(hitsAtTecPlus > 0){
00930 //        derivatives[hit][2] = 0.0;
00931 //      }
00932       }
00933       // TECplus
00934       else if(gHitZprime[hit] > gBeamSplitterZprime){
00935         derivatives[hit][1] = globPtrack[hit].z();
00936 //      if(hitsAtTecPlus > 0){
00937 //        derivatives[hit][2] = - (globPtrack[hit].z() - gBeamSplitterZprime) / gBeamR;
00938 //      }
00939       }
00940       // Barrel
00941       else{
00942         derivatives[hit][1] = globPtrack[hit].z() - gBarrelModuleOffset[hit-gHitsAtTecMinus];
00943 //      if(hitsAtTecPlus > 0){
00944 //        derivatives[hit][2] = 0.0;
00945 //      }
00946       }
00947     }
00948   }
00949 
00950   unsigned int firstFixedParam(covMatrix.num_col()); // FIXME! --> no, is fine!!!
00951 //   unsigned int firstFixedParam = nPedeParams - 1;
00952 //   if(beam.isAlignmentTube() && hitsAtTecPlus == 0) firstFixedParam = nPedeParams;
00953 //   cout << "first fixed parameter: " << firstFixedParam << endl;
00954   // set fit results
00955   beam.setParameters(paramType, params, covMatrix, derivatives, firstFixedParam, chi2);
00956 
00957   return true; // return false in case of problems
00958 }
00959 
00960 
00961 //---------------------------------------------------------------------------------------
00962 //define this as a plug-in
00963 DEFINE_FWK_MODULE(TkLasBeamFitter);