|
|
|
#include <Alignment/HIPAlignmentAlgorithm/interface/HIPAlignmentAlgorithm.h>
Definition at line 12 of file HIPAlignmentAlgorithm.h.
| HIPAlignmentAlgorithm::HIPAlignmentAlgorithm | ( | const edm::ParameterSet & | cfg | ) |
Constructor.
Definition at line 32 of file HIPAlignmentAlgorithm.cc.
References dummy, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), isCollector, edm::es::l(), AlignableObjectId::nameToType(), outfile, outfile2, outpath, salignedfile, siterationfile, smisalignedfile, sparameterfile, struefile, suvarfile, theAPEParameterSet, theApplyAPE, theCollectorNJobs, theCollectorPath, theCurrentPrescale, theEventPrescale, theLevels, theMaxAllowedHitPull, theMaxRelParameterError, and theMinimumNumberOfHits.
00032 : 00033 AlignmentAlgorithmBase( cfg ) 00034 { 00035 00036 // parse parameters 00037 00038 verbose = cfg.getParameter<bool>("verbosity"); 00039 00040 outpath = cfg.getParameter<string>("outpath"); 00041 outfile = cfg.getParameter<string>("outfile"); 00042 outfile2 = cfg.getParameter<string>("outfile2"); 00043 struefile = cfg.getParameter<string>("trueFile"); 00044 smisalignedfile = cfg.getParameter<string>("misalignedFile"); 00045 salignedfile = cfg.getParameter<string>("alignedFile"); 00046 siterationfile = cfg.getParameter<string>("iterationFile"); 00047 suvarfile = cfg.getParameter<string>("uvarFile"); 00048 sparameterfile = cfg.getParameter<string>("parameterFile"); 00049 00050 outfile =outpath+outfile; 00051 outfile2 =outpath+outfile2; 00052 struefile =outpath+struefile; 00053 smisalignedfile=outpath+smisalignedfile; 00054 salignedfile =outpath+salignedfile; 00055 siterationfile =outpath+siterationfile; 00056 suvarfile =outpath+suvarfile; 00057 sparameterfile =outpath+sparameterfile; 00058 00059 // parameters for APE 00060 theApplyAPE = cfg.getParameter<bool>("applyAPE"); 00061 theAPEParameterSet = cfg.getParameter<std::vector<edm::ParameterSet> >("apeParam"); 00062 00063 theMaxAllowedHitPull = cfg.getParameter<double>("maxAllowedHitPull"); 00064 theMinimumNumberOfHits = cfg.getParameter<int>("minimumNumberOfHits"); 00065 theMaxRelParameterError = cfg.getParameter<double>("maxRelParameterError"); 00066 00067 // for collector mode (parallel processing) 00068 isCollector=cfg.getParameter<bool>("collectorActive"); 00069 theCollectorNJobs=cfg.getParameter<int>("collectorNJobs"); 00070 theCollectorPath=cfg.getParameter<string>("collectorPath"); 00071 if (isCollector) edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Collector mode"; 00072 00073 theEventPrescale = cfg.getParameter<int>("eventPrescale"); 00074 theCurrentPrescale = theEventPrescale; 00075 00076 AlignableObjectId dummy; 00077 00078 const std::vector<std::string>& levels = cfg.getUntrackedParameter< std::vector<std::string> >("surveyResiduals"); 00079 00080 for (unsigned int l = 0; l < levels.size(); ++l) 00081 { 00082 theLevels.push_back( dummy.nameToType(levels[l]) ); 00083 } 00084 00085 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] constructed."; 00086 00087 }
| HIPAlignmentAlgorithm::~HIPAlignmentAlgorithm | ( | ) | [inline] |
Definition at line 644 of file HIPAlignmentAlgorithm.cc.
References m2_Eta, m2_Id, m2_Layer, m2_Nhit, m2_ObjId, m2_Phi, m2_Type, m2_Xpos, m2_Ypos, m2_Zpos, m_Chi2n, m_Eta, m_Nhits, m_Ntracks, m_Phi, m_Pt, outfile, outfile2, theFile, theFile2, theIteration, theTree, and theTree2.
Referenced by startNewLoop().
00645 { 00646 // create ROOT files 00647 theFile = new TFile(outfile.c_str(),"update"); 00648 theFile->cd(); 00649 00650 // book event-wise ROOT Tree 00651 00652 TString tname="T1"; 00653 char iterString[5]; 00654 sprintf(iterString, "%i",theIteration); 00655 tname.Append("_"); 00656 tname.Append(iterString); 00657 00658 theTree = new TTree(tname,"Eventwise tree"); 00659 00660 //theTree->Branch("Run", &m_Run, "Run/I"); 00661 //theTree->Branch("Event", &m_Event, "Event/I"); 00662 theTree->Branch("Ntracks", &m_Ntracks, "Ntracks/I"); 00663 theTree->Branch("Nhits", m_Nhits, "Nhits[Ntracks]/I"); 00664 theTree->Branch("Pt", m_Pt, "Pt[Ntracks]/F"); 00665 theTree->Branch("Eta", m_Eta, "Eta[Ntracks]/F"); 00666 theTree->Branch("Phi", m_Phi, "Phi[Ntracks]/F"); 00667 theTree->Branch("Chi2n", m_Chi2n, "Chi2n[Ntracks]/F"); 00668 00669 // book Alignable-wise ROOT Tree 00670 00671 theFile2 = new TFile(outfile2.c_str(),"update"); 00672 theFile2->cd(); 00673 00674 theTree2 = new TTree("T2","Alignablewise tree"); 00675 00676 theTree2->Branch("Nhit", &m2_Nhit, "Nhit/I"); 00677 theTree2->Branch("Type", &m2_Type, "Type/I"); 00678 theTree2->Branch("Layer", &m2_Layer, "Layer/I"); 00679 theTree2->Branch("Xpos", &m2_Xpos, "Xpos/F"); 00680 theTree2->Branch("Ypos", &m2_Ypos, "Ypos/F"); 00681 theTree2->Branch("Zpos", &m2_Zpos, "Zpos/F"); 00682 theTree2->Branch("Eta", &m2_Eta, "Eta/F"); 00683 theTree2->Branch("Phi", &m2_Phi, "Phi/F"); 00684 theTree2->Branch("Id", &m2_Id, "Id/i"); 00685 theTree2->Branch("ObjId", &m2_ObjId, "ObjId/I"); 00686 00687 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Root trees booked."; 00688 00689 }
| double HIPAlignmentAlgorithm::calcAPE | ( | double * | par, | |
| int | iter, | |||
| double | function | |||
| ) | [private] |
Definition at line 622 of file HIPAlignmentAlgorithm.cc.
References funct::exp(), max, and funct::pow().
Referenced by setAlignmentPositionError().
00623 { 00624 double diter=(double)iter; 00625 00626 // The following floating-point equality check is safe because this 00627 // "0." and this "1." are generated by the compiler, in the very 00628 // same file. Whatever approximization scheme it uses to turn "1." 00629 // into 0.9999999999998 in the HIPAlignmentAlgorithm::initialize is 00630 // also used here. If I'm wrong, you'll get an assertion. 00631 if (function == 0.) { 00632 return max(par[1],par[0]+((par[1]-par[0])/par[2])*diter); 00633 } 00634 else if (function == 1.) { 00635 return max(0.,par[0]*(exp(-pow(diter,par[1])/par[2]))); 00636 } 00637 else assert(false); // should have been caught in the constructor 00638 }
Definition at line 770 of file HIPAlignmentAlgorithm.cc.
References funct::abs(), Alignable::alignmentParameters(), AlignmentParameters::cloneFromSelected(), i, HIPUserVariables::jtve, HIPUserVariables::jtvj, HIPUserVariables::nhit, params, Alignable::setAlignmentParameters(), AlignmentParameters::setValid(), funct::sqrt(), theMaxRelParameterError, theMinimumNumberOfHits, and AlignmentParameters::userVariables().
Referenced by terminate().
00771 { 00772 00773 // Alignment parameters 00774 AlignmentParameters* par = ali->alignmentParameters(); 00775 // access user variables 00776 HIPUserVariables* uservar = 00777 dynamic_cast<HIPUserVariables*>(par->userVariables()); 00778 int nhit = uservar->nhit; 00779 00780 if (nhit < theMinimumNumberOfHits) { 00781 par->setValid(false); 00782 return false; 00783 } 00784 00785 AlgebraicSymMatrix jtvj = uservar->jtvj; 00786 AlgebraicVector jtve = uservar->jtve; 00787 int ierr; 00788 AlgebraicSymMatrix jtvjinv=jtvj.inverse(ierr); 00789 if (ierr !=0) { 00790 edm::LogError("Alignment") << "Matrix inversion failed!"; 00791 return false; 00792 } 00793 00794 // these are the alignment corrections+covariance (for selected params) 00795 AlgebraicVector params = - (jtvjinv * jtve); 00796 AlgebraicSymMatrix cov = jtvjinv; 00797 00798 edm::LogInfo("Alignment") << "parameters " << params; 00799 00800 // errors of parameters 00801 int npar=params.num_row(); 00802 AlgebraicVector paramerr(npar); 00803 AlgebraicVector relerr(npar); 00804 for (int i=0;i<npar;i++) { 00805 if (abs(cov[i][i])>0) paramerr[i]=sqrt(abs(cov[i][i])); 00806 else paramerr[i]=params[i]; 00807 relerr[i] = abs(paramerr[i]/params[i]); 00808 if (relerr[i] >= theMaxRelParameterError) { 00809 params[i]=0; 00810 paramerr[i]=0; 00811 } 00812 } 00813 00814 // store alignment parameters 00815 AlignmentParameters* parnew = par->cloneFromSelected(params,cov); 00816 ali->setAlignmentParameters(parnew); 00817 parnew->setValid(true); 00818 return true; 00819 }
Definition at line 823 of file HIPAlignmentAlgorithm.cc.
References Alignable::alignmentParameters(), AlignmentParameterStore::attachUserVariables(), HIPUserVariables::clone(), lat::endl(), ioerr, it, HIPUserVariables::jtve, HIPUserVariables::jtvj, HIPUserVariables::nhit, HIPUserVariablesIORoot::readHIPUserVariables(), ss, theAlignables, theAlignmentParameterStore, theCollectorNJobs, theCollectorPath, theIteration, and AlignmentParameters::userVariables().
Referenced by startNewLoop().
00824 { 00825 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::collector] called for iteration " << theIteration <<endl; 00826 00827 HIPUserVariablesIORoot HIPIO; 00828 00829 for (int ijob=1; ijob<=theCollectorNJobs; ijob++) { 00830 00831 edm::LogWarning("Alignment") <<"reading uservar for job " << ijob; 00832 00833 stringstream ss; 00834 string str; 00835 ss << ijob; 00836 ss >> str; 00837 string uvfile = theCollectorPath+"/job"+str+"/IOUserVariables.root"; 00838 00839 vector<AlignmentUserVariables*> uvarvec = 00840 HIPIO.readHIPUserVariables (theAlignables,(char*)uvfile.c_str(), 00841 theIteration,ioerr); 00842 00843 if (ioerr!=0) { 00844 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::collector] could not read user variable files for job" << ijob; 00845 continue; 00846 } 00847 00848 // add 00849 vector<AlignmentUserVariables*> uvarvecadd; 00850 vector<AlignmentUserVariables*>::const_iterator iuvarnew=uvarvec.begin(); 00851 for (vector<Alignable*>::const_iterator it=theAlignables.begin(); 00852 it!=theAlignables.end(); it++) { 00853 Alignable* ali = *it; 00854 AlignmentParameters* ap = ali->alignmentParameters(); 00855 00856 HIPUserVariables* uvarold = 00857 dynamic_cast<HIPUserVariables*>(ap->userVariables()); 00858 HIPUserVariables* uvarnew = 00859 dynamic_cast<HIPUserVariables*>(*iuvarnew); 00860 00861 HIPUserVariables* uvar = uvarold->clone(); 00862 if (uvarnew!=0) { 00863 uvar->nhit=(uvarold->nhit)+(uvarnew->nhit); 00864 uvar->jtvj=(uvarold->jtvj)+(uvarnew->jtvj); 00865 uvar->jtve=(uvarold->jtve)+(uvarnew->jtve); 00866 delete uvarnew; 00867 } 00868 00869 uvarvecadd.push_back(uvar); 00870 iuvarnew++; 00871 } 00872 00873 theAlignmentParameterStore->attachUserVariables(theAlignables, 00874 uvarvecadd,ioerr); 00875 00876 } 00877 00878 }
Definition at line 694 of file HIPAlignmentAlgorithm.cc.
References Alignable::alignableObjectId(), Alignable::alignmentParameters(), PV3DBase< T, PVType, FrameType >::eta(), Alignable::id(), AlignmentParameters::isValid(), it, m2_Eta, m2_Id, m2_Layer, m2_Nhit, m2_ObjId, m2_Phi, m2_Type, m2_Xpos, m2_Ypos, m2_Zpos, HIPUserVariables::nhit, AlignmentParameters::parameters(), pars, PV3DBase< T, PVType, FrameType >::phi(), GloballyPositioned< T >::position(), Alignable::surface(), theAlignables, theAlignmentParameterStore, theFile2, theTree2, AlignmentParameterStore::typeAndLayer(), AlignmentParameters::userVariables(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by terminate().
00695 { 00696 theFile2->cd(); 00697 00698 int naligned=0; 00699 00700 for(vector<Alignable*>::const_iterator 00701 it=theAlignables.begin(); it!=theAlignables.end(); it++) { 00702 Alignable* ali=(*it); 00703 AlignmentParameters* dap = ali->alignmentParameters(); 00704 00705 // consider only those parameters classified as 'valid' 00706 if (dap->isValid()) { 00707 00708 // get number of hits from user variable 00709 HIPUserVariables* uservar = 00710 dynamic_cast<HIPUserVariables*>(dap->userVariables()); 00711 m2_Nhit = uservar->nhit; 00712 00713 // get type/layer 00714 std::pair<int,int> tl=theAlignmentParameterStore->typeAndLayer(ali); 00715 m2_Type=tl.first; 00716 m2_Layer=tl.second; 00717 00718 // get identifier (as for IO) 00719 m2_Id = ali->id(); 00720 m2_ObjId = ali->alignableObjectId(); 00721 00722 // get position 00723 GlobalPoint pos=ali->surface().position(); 00724 m2_Xpos=pos.x(); 00725 m2_Ypos=pos.y(); 00726 m2_Zpos=pos.z(); 00727 m2_Eta=pos.eta(); 00728 m2_Phi=pos.phi(); 00729 00730 AlgebraicVector pars=dap->parameters(); 00731 00732 if (verbose) 00733 { 00734 edm::LogVerbatim("Alignment") 00735 << "------------------------------------------------------------------------\n" 00736 << " ALIGNABLE: " << setw(6) << naligned 00737 << '\n' 00738 << "hits: " << setw(4) << m2_Nhit 00739 << " type: " << setw(4) << m2_Type 00740 << " layer: " << setw(4) << m2_Layer 00741 << " id: " << setw(4) << m2_Id 00742 << " objId: " << setw(4) << m2_ObjId 00743 << '\n' 00744 << fixed << setprecision(5) 00745 << "x,y,z: " 00746 << setw(12) << m2_Xpos 00747 << setw(12) << m2_Ypos 00748 << setw(12) << m2_Zpos 00749 << " eta,phi: " 00750 << setw(12) << m2_Eta 00751 << setw(12) << m2_Phi 00752 << '\n' 00753 << "params: " 00754 << setw(12) << pars[0] 00755 << setw(12) << pars[1] 00756 << setw(12) << pars[2] 00757 << setw(12) << pars[3] 00758 << setw(12) << pars[4] 00759 << setw(12) << pars[5]; 00760 } 00761 00762 naligned++; 00763 theTree2->Fill(); 00764 } 00765 } 00766 }
| void HIPAlignmentAlgorithm::initialize | ( | const edm::EventSetup & | setup, | |
| AlignableTracker * | tracker, | |||
| AlignableMuon * | muon, | |||
| AlignmentParameterStore * | store | |||
| ) | [virtual] |
Call at beginning of job.
Implements AlignmentAlgorithmBase.
Definition at line 92 of file HIPAlignmentAlgorithm.cc.
References AlignmentParameterSelector::addSelections(), AlignmentParameterStore::alignables(), AlignmentParameterSelector::clear(), lat::endl(), Exception, edm::ParameterSet::getParameter(), i, AlignmentParameterSelector::selectedAlignables(), theAlignableDetAccessor, theAlignables, theAlignmentParameterStore, theAPEParameters, theAPEParameterSet, and theApplyAPE.
00095 { 00096 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Initializing..."; 00097 00098 // accessor Det->AlignableDet 00099 if ( !muon ) 00100 theAlignableDetAccessor = new AlignableNavigator(tracker); 00101 else if ( !tracker ) 00102 theAlignableDetAccessor = new AlignableNavigator(muon); 00103 else 00104 theAlignableDetAccessor = new AlignableNavigator(tracker,muon); 00105 00106 // set alignmentParameterStore 00107 theAlignmentParameterStore=store; 00108 00109 // get alignables 00110 theAlignables = theAlignmentParameterStore->alignables(); 00111 00112 // clear theAPEParameters, if necessary 00113 theAPEParameters.clear(); 00114 00115 // get APE parameters 00116 if (theApplyAPE) { 00117 AlignmentParameterSelector selector(tracker, muon); 00118 for (std::vector<edm::ParameterSet>::const_iterator setiter = theAPEParameterSet.begin(); setiter != theAPEParameterSet.end(); ++setiter) { 00119 std::vector<Alignable*> alignables; 00120 00121 selector.clear(); 00122 edm::ParameterSet selectorPSet = setiter->getParameter<edm::ParameterSet>("Selector"); 00123 std::vector<std::string> alignParams = selectorPSet.getParameter<std::vector<std::string> >("alignParams"); 00124 if (alignParams.size() == 1 && alignParams[0] == std::string("selected")) { 00125 alignables = theAlignables; 00126 } 00127 else { 00128 selector.addSelections(selectorPSet); 00129 alignables = selector.selectedAlignables(); 00130 } 00131 00132 std::vector<double> apeSPar = setiter->getParameter<std::vector<double> >("apeSPar"); 00133 std::vector<double> apeRPar = setiter->getParameter<std::vector<double> >("apeRPar"); 00134 std::string function = setiter->getParameter<std::string>("function"); 00135 00136 if (apeSPar.size() != 3 || apeRPar.size() != 3) 00137 throw cms::Exception("BadConfig") << "apeSPar and apeRPar must have 3 values each" << std::endl; 00138 00139 for (std::vector<double>::const_iterator i = apeRPar.begin(); i != apeRPar.end(); ++i) { 00140 apeSPar.push_back(*i); 00141 } 00142 00143 if (function == std::string("linear")) { 00144 apeSPar.push_back(0.); // c.f. note in calcAPE 00145 } 00146 else if (function == std::string("exponential")) { 00147 apeSPar.push_back(1.); // c.f. note in calcAPE 00148 } 00149 else { 00150 throw cms::Exception("BadConfig") << "APE function must be \"linear\" or \"exponential\"." << std::endl; 00151 } 00152 00153 theAPEParameters.push_back(std::pair<std::vector<Alignable*>, std::vector<double> >(alignables, apeSPar)); 00154 } 00155 } 00156 }
| int HIPAlignmentAlgorithm::readIterationFile | ( | std::string | filename | ) | [private] |
Referenced by startNewLoop().
| void HIPAlignmentAlgorithm::run | ( | const edm::EventSetup & | setup, | |
| const ConstTrajTrackPairCollection & | tracks | |||
| ) |
Run the algorithm on trajectories and tracks.
Definition at line 353 of file HIPAlignmentAlgorithm.cc.
References CompositeAlignmentParameters::alignableFromAlignableDet(), AlignableNavigator::alignableFromGeomDet(), AlignableNavigator::alignablesFromHits(), Alignable::alignmentParameters(), TrajectoryMeasurement::backwardPredictedState(), TrajectoryStateCombiner::combine(), AlignableNavigator::detAndSubdetInMap(), eta, reco::TrackBase::eta(), TrajectoryMeasurement::forwardPredictedState(), TrackingRecHit::geographicalId(), isCollector, TrackingRecHit::isValid(), TrajectoryStateOnSurface::isValid(), it, m_Chi2n, m_Eta, m_Nhits, m_Ntracks, m_Phi, m_Pt, MAXREC, Trajectory::measurements(), reco::TrackBase::normalizedChi2(), reco::TrackBase::numberOfValidHits(), params, phi, reco::TrackBase::phi(), reco::TrackBase::pt(), TrajectoryMeasurement::recHit(), AlignmentParameters::selectedDerivatives(), AlignmentParameterStore::selectParameters(), funct::sqrt(), theAlignableDetAccessor, theAlignmentParameterStore, theCurrentPrescale, theEventPrescale, theFile, theMaxAllowedHitPull, theTree, track, AlignmentParameters::userVariables(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().
00355 { 00356 if (isCollector) return; 00357 00358 TrajectoryStateCombiner tsoscomb; 00359 00360 int itr=0; 00361 m_Ntracks=0; 00362 00363 theFile->cd(); 00364 00365 // loop over tracks 00366 for( ConstTrajTrackPairCollection::const_iterator it=tracks.begin(); 00367 it!=tracks.end();it++) { 00368 00369 const Trajectory* traj = (*it).first; 00370 const reco::Track* track = (*it).second; 00371 00372 float pt = track->pt(); 00373 float eta = track->eta(); 00374 float phi = track->phi(); 00375 float chi2n = track->normalizedChi2(); 00376 int nhit = track->numberOfValidHits(); 00377 00378 if (verbose) edm::LogInfo("Alignment") << "New track pt,eta,phi,chi2n,hits: " << pt <<","<< eta <<","<< phi <<","<< chi2n << ","<<nhit; 00379 00380 // fill track parameters in root tree 00381 if (itr<MAXREC) { 00382 m_Nhits[itr]=nhit; 00383 m_Pt[itr]=pt; 00384 m_Eta[itr]=eta; 00385 m_Phi[itr]=phi; 00386 m_Chi2n[itr]=chi2n; 00387 itr++; 00388 m_Ntracks=itr; 00389 } 00390 00391 vector<const TransientTrackingRecHit*> hitvec; 00392 vector<TrajectoryStateOnSurface> tsosvec; 00393 00394 // loop over measurements 00395 vector<TrajectoryMeasurement> measurements = traj->measurements(); 00396 for (vector<TrajectoryMeasurement>::iterator im=measurements.begin(); 00397 im!=measurements.end(); im++) { 00398 TrajectoryMeasurement meas = *im; 00399 const TransientTrackingRecHit* hit = &(*meas.recHit()); 00400 if (hit->isValid() && theAlignableDetAccessor->detAndSubdetInMap( hit->geographicalId() )) { 00401 // this is the updated state (including the current hit) 00402 //TrajectoryStateOnSurface tsos=meas.updatedState(); 00403 // combine fwd and bwd predicted state to get state 00404 // which excludes current hit 00405 TrajectoryStateOnSurface tsos = 00406 tsoscomb.combine(meas.forwardPredictedState(), 00407 meas.backwardPredictedState()); 00408 if (tsos.isValid()) 00409 { 00410 hitvec.push_back(hit); 00411 tsosvec.push_back(tsos); 00412 } 00413 } 00414 } 00415 00416 // transform RecHit vector to AlignableDet vector 00417 vector <AlignableDetOrUnitPtr> alidetvec = 00418 theAlignableDetAccessor->alignablesFromHits(hitvec); 00419 00420 // get concatenated alignment parameters for list of alignables 00421 CompositeAlignmentParameters aap = 00422 theAlignmentParameterStore->selectParameters(alidetvec); 00423 00424 vector<TrajectoryStateOnSurface>::const_iterator itsos=tsosvec.begin(); 00425 vector<const TransientTrackingRecHit*>::const_iterator ihit=hitvec.begin(); 00426 00427 // loop over vectors(hit,tsos) 00428 while (itsos != tsosvec.end()) 00429 { 00430 // get AlignableDet for this hit 00431 const GeomDet* det=(*ihit)->det(); 00432 AlignableDetOrUnitPtr alidet = 00433 theAlignableDetAccessor->alignableFromGeomDet(det); 00434 00435 // get relevant Alignable 00436 Alignable* ali=aap.alignableFromAlignableDet(alidet); 00437 00438 if (ali!=0) { 00439 // get trajectory impact point 00440 LocalPoint alvec = (*itsos).localPosition(); 00441 AlgebraicVector pos(2); 00442 pos[0]=alvec.x(); // local x 00443 pos[1]=alvec.y(); // local y 00444 00445 // get impact point covariance 00446 AlgebraicSymMatrix ipcovmat(2); 00447 ipcovmat[0][0] = (*itsos).localError().positionError().xx(); 00448 ipcovmat[1][1] = (*itsos).localError().positionError().yy(); 00449 ipcovmat[0][1] = (*itsos).localError().positionError().xy(); 00450 00451 // get hit local position and covariance 00452 AlgebraicVector coor(2); 00453 coor[0] = (*ihit)->localPosition().x(); 00454 coor[1] = (*ihit)->localPosition().y(); 00455 00456 AlgebraicSymMatrix covmat(2); 00457 covmat[0][0] = (*ihit)->localPositionError().xx(); 00458 covmat[1][1] = (*ihit)->localPositionError().yy(); 00459 covmat[0][1] = (*ihit)->localPositionError().xy(); 00460 00461 // add hit and impact point covariance matrices 00462 covmat = covmat + ipcovmat; 00463 00464 // calculate the x pull and y pull of this hit 00465 double xpull = 0.; 00466 double ypull = 0.; 00467 if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0])); 00468 if (covmat[1][1] != 0.) ypull = (pos[1] - coor[1])/sqrt(fabs(covmat[1][1])); 00469 00470 // get Alignment Parameters 00471 AlignmentParameters* params = ali->alignmentParameters(); 00472 // get derivatives 00473 AlgebraicMatrix derivs=params->selectedDerivatives(*itsos,alidet); 00474 00475 // invert covariance matrix 00476 int ierr; 00477 covmat.invert(ierr); 00478 if (ierr != 0) { 00479 edm::LogError("Alignment") << "Matrix inversion failed!"; 00480 return; 00481 } 00482 00483 bool useThisHit = (theMaxAllowedHitPull <= 0.); 00484 00485 // ignore track minus center-of-chamber "residual" from 1d hits (only muon drift tubes) 00486 if ((*ihit)->dimension() == 1) { 00487 covmat[1][1] = 0.; 00488 covmat[0][1] = 0.; 00489 00490 useThisHit = useThisHit || (fabs(xpull) < theMaxAllowedHitPull); 00491 } 00492 else { 00493 useThisHit = useThisHit || (fabs(xpull) < theMaxAllowedHitPull && fabs(ypull) < theMaxAllowedHitPull); 00494 } 00495 00496 if (useThisHit) { 00497 // calculate user parameters 00498 int npar=derivs.num_row(); 00499 AlgebraicSymMatrix thisjtvj(npar); 00500 AlgebraicVector thisjtve(npar); 00501 thisjtvj=covmat.similarity(derivs); 00502 thisjtve=derivs * covmat * (pos-coor); 00503 00504 // access user variables (via AlignmentParameters) 00505 HIPUserVariables* uservar = 00506 dynamic_cast<HIPUserVariables*>(params->userVariables()); 00507 uservar->jtvj += thisjtvj; 00508 uservar->jtve += thisjtve; 00509 uservar->nhit ++; 00510 } 00511 } 00512 00513 itsos++; 00514 ihit++; 00515 } 00516 00517 } // end of track loop 00518 00519 // fill eventwise root tree (with prescale defined in pset) 00520 theCurrentPrescale--; 00521 if (theCurrentPrescale<=0) { 00522 theTree->Fill(); 00523 theCurrentPrescale=theEventPrescale; 00524 } 00525 00526 }
Definition at line 569 of file HIPAlignmentAlgorithm.cc.
References calcAPE(), i, pars, AlignmentParameterStore::setAlignmentPositionError(), theAlignmentParameterStore, theAPEParameters, theApplyAPE, and theIteration.
Referenced by startNewLoop().
00570 { 00571 00572 00573 // Check if user wants to override APE 00574 if ( !theApplyAPE ) 00575 { 00576 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] No APE applied"; 00577 return; // NO APE APPLIED 00578 } 00579 00580 00581 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] Apply APE!"; 00582 00583 double apeSPar[3], apeRPar[3]; 00584 for (std::vector<std::pair<std::vector<Alignable*>, std::vector<double> > >::const_iterator alipars = theAPEParameters.begin(); 00585 alipars != theAPEParameters.end(); 00586 ++alipars) { 00587 const std::vector<Alignable*> &alignables = alipars->first; 00588 const std::vector<double> &pars = alipars->second; 00589 00590 apeSPar[0] = pars[0]; 00591 apeSPar[1] = pars[1]; 00592 apeSPar[2] = pars[2]; 00593 apeRPar[0] = pars[3]; 00594 apeRPar[1] = pars[4]; 00595 apeRPar[2] = pars[5]; 00596 00597 double function = pars[6]; 00598 00599 // Printout for debug 00600 printf("APE: %d alignables\n", alignables.size()); 00601 for ( int i=0; i<21; ++i ) { 00602 double apelinstest=calcAPE(apeSPar,i,0.); 00603 double apeexpstest=calcAPE(apeSPar,i,1.); 00604 double apelinrtest=calcAPE(apeRPar,i,0.); 00605 double apeexprtest=calcAPE(apeRPar,i,1.); 00606 printf("APE: iter slin sexp rlin rexp: %5d %12.5f %12.5f %12.5f %12.5f\n", 00607 i,apelinstest,apeexpstest,apelinrtest,apeexprtest); 00608 } 00609 00610 // set APE 00611 double apeshift=calcAPE(apeSPar,theIteration,function); 00612 double aperot =calcAPE(apeRPar,theIteration,function); 00613 theAlignmentParameterStore->setAlignmentPositionError( alignables, apeshift, aperot ); 00614 } 00615 00616 }
Called at start of new loop.
Reimplemented from AlignmentAlgorithmBase.
Definition at line 159 of file HIPAlignmentAlgorithm.cc.
References Alignable::alignmentParameters(), AlignmentParameterStore::applyAlignableAbsolutePositions(), bookRoot(), collector(), i, SurveyResidual::inverseCovariance(), ioerr, isCollector, it, HIPUserVariables::jtve, HIPUserVariables::jtvj, edm::es::l(), AlignmentParameters::numSelected(), AlignmentIORoot::readAlignableAbsolutePositions(), readIterationFile(), res, salignedfile, SurveyResidual::sensorResidual(), setAlignmentPositionError(), AlignmentParameters::setUserVariables(), siterationfile, smisalignedfile, struefile, theAlignables, theAlignmentParameterStore, theIO, theIteration, theLevels, AlignmentParameters::userVariables(), SurveyResidual::valid(), AlignmentIORoot::writeAlignableAbsolutePositions(), and AlignmentIORoot::writeAlignableOriginalPositions().
00160 { 00161 00162 // iterate over all alignables and attach user variables 00163 for( vector<Alignable*>::const_iterator it=theAlignables.begin(); 00164 it!=theAlignables.end(); it++ ) 00165 { 00166 AlignmentParameters* ap = (*it)->alignmentParameters(); 00167 int npar=ap->numSelected(); 00168 HIPUserVariables* userpar = new HIPUserVariables(npar); 00169 ap->setUserVariables(userpar); 00170 } 00171 00172 // try to read in alignment parameters from a previous iteration 00173 AlignablePositions theAlignablePositionsFromFile = 00174 theIO.readAlignableAbsolutePositions(theAlignables, 00175 (char*)salignedfile.c_str(),-1,ioerr); 00176 00177 int numAlignablesFromFile = theAlignablePositionsFromFile.size(); 00178 00179 if (numAlignablesFromFile==0) { // file not there: first iteration 00180 00181 // set iteration number to 1 00182 if (isCollector) theIteration=0; 00183 else theIteration=1; 00184 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] File not found => iteration "<<theIteration; 00185 00186 // get true (de-misaligned positions) and write to root file 00187 // hardcoded iteration=1 00188 theIO.writeAlignableOriginalPositions(theAlignables, 00189 (char*)struefile.c_str(),1,false,ioerr); 00190 00191 // get misaligned positions and write to root file 00192 // hardcoded iteration=1 00193 theIO.writeAlignableAbsolutePositions(theAlignables, 00194 (char*)smisalignedfile.c_str(),1,false,ioerr); 00195 00196 } 00197 00198 else { // there have been previous iterations 00199 00200 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Alignables Read " 00201 << numAlignablesFromFile; 00202 00203 // get iteration number from file 00204 theIteration = readIterationFile(siterationfile); 00205 00206 // increase iteration 00207 theIteration++; 00208 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Iteration increased by one!"; 00209 00210 // now apply psotions of file from prev iteration 00211 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Apply positions from file ..."; 00212 theAlignmentParameterStore->applyAlignableAbsolutePositions(theAlignables, 00213 theAlignablePositionsFromFile,ioerr); 00214 00215 } 00216 00217 edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Current Iteration number: " 00218 << theIteration; 00219 00220 if (theLevels.size() > 0) 00221 { 00222 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint"; 00223 00224 unsigned int nAlignable = theAlignables.size(); 00225 00226 for (unsigned int i = 0; i < nAlignable; ++i) 00227 { 00228 const Alignable* ali = theAlignables[i]; 00229 00230 AlignmentParameters* ap = ali->alignmentParameters(); 00231 00232 HIPUserVariables* uservar = 00233 dynamic_cast<HIPUserVariables*>(ap->userVariables()); 00234 00235 for (unsigned int l = 0; l < theLevels.size(); ++l) 00236 { 00237 SurveyResidual res(*ali, theLevels[l], true); 00238 00239 if ( res.valid() ) 00240 { 00241 AlgebraicSymMatrix invCov = res.inverseCovariance(); 00242 00243 uservar->jtvj += invCov; 00244 uservar->jtve += invCov * res.sensorResidual(); 00245 } 00246 } 00247 00248 // align::LocalVectors residuals = res1.pointsResidual(); 00249 00250 // unsigned int nPoints = residuals.size(); 00251 00252 // for (unsigned int k = 0; k < nPoints; ++k) 00253 // { 00254 // AlgebraicMatrix J = term->survey()->derivatives(k); 00255 // AlgebraicVector e(3); // local residual 00256 00257 // const align::LocalVector& lr = residuals[k]; 00258 00259 // e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z(); 00260 00261 // uservar->jtvj += invCov1.similarity(J); 00262 // uservar->jtve += J * (invCov1 * e); 00263 // } 00264 00265 } 00266 } 00267 00268 // set alignment position error 00269 setAlignmentPositionError(); 00270 00271 // book root trees 00272 bookRoot(); 00273 00274 // run collector job if we are in parallel mode 00275 if (isCollector) collector(); 00276 00277 }
Call at end of job.
Implements AlignmentAlgorithmBase.
Definition at line 281 of file HIPAlignmentAlgorithm.cc.
References Alignable::alignmentParameters(), AlignmentParameterStore::applyParameters(), calcParameters(), fillRoot(), ioerr, it, salignedfile, AlignmentParameters::setValid(), siterationfile, sparameterfile, suvarfile, parsecf::pyparsing::test(), theAlignables, theAlignmentParameterStore, theFile, theFile2, theIO, theIteration, theTree, theTree2, AlignmentIORoot::writeAlignableAbsolutePositions(), AlignmentIORoot::writeAlignmentParameters(), HIPUserVariablesIORoot::writeHIPUserVariables(), and writeIterationFile().
00282 { 00283 00284 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Terminating"; 00285 00286 // write user variables 00287 HIPUserVariablesIORoot HIPIO; 00288 HIPIO.writeHIPUserVariables (theAlignables,(char*)suvarfile.c_str(), 00289 theIteration,false,ioerr); 00290 00291 // now calculate alignment corrections ... 00292 int ialigned=0; 00293 // iterate over alignment parameters 00294 for(vector<Alignable*>::const_iterator 00295 it=theAlignables.begin(); it!=theAlignables.end(); it++) { 00296 Alignable* ali=(*it); 00297 // Alignment parameters 00298 AlignmentParameters* par = ali->alignmentParameters(); 00299 // try to calculate parameters 00300 bool test = calcParameters(ali); 00301 // if successful, apply parameters 00302 if (test) { 00303 edm::LogInfo("Alignment") << "now apply params"; 00304 theAlignmentParameterStore->applyParameters(ali); 00305 // set these parameters 'valid' 00306 ali->alignmentParameters()->setValid(true); 00307 // increase counter 00308 ialigned++; 00309 } 00310 else par->setValid(false); 00311 } 00312 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::terminate] Aligned units: " << ialigned; 00313 00314 // fill alignable wise root tree 00315 fillRoot(); 00316 00317 edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Writing aligned parameters to file: " << theAlignables.size(); 00318 00319 // write new absolute positions to disk 00320 theIO.writeAlignableAbsolutePositions(theAlignables, 00321 (char*)salignedfile.c_str(),theIteration,false,ioerr); 00322 00323 // write alignment parameters to disk 00324 theIO.writeAlignmentParameters(theAlignables, 00325 (char*)sparameterfile.c_str(),theIteration,false,ioerr); 00326 00327 // write iteration number to file 00328 writeIterationFile(siterationfile,theIteration); 00329 00330 00331 // write out trees and close root file 00332 00333 // eventwise tree 00334 theFile->cd(); 00335 theTree->Write(); 00336 theFile->Close(); 00337 delete theFile; 00338 00339 // alignable-wise tree is only filled once 00340 //if ((!isCollector && theIteration==1)|| 00341 // ( isCollector && theIteration==2)) { 00342 if (theIteration==1) { // only for 1st iteration 00343 theFile2->cd(); 00344 theTree2->Write(); 00345 theFile2->Close(); 00346 delete theFile2; 00347 } 00348 00349 }
Referenced by terminate().
int HIPAlignmentAlgorithm::ioerr [private] |
Definition at line 57 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), startNewLoop(), and terminate().
bool HIPAlignmentAlgorithm::isCollector [private] |
Definition at line 78 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), run(), and startNewLoop().
float HIPAlignmentAlgorithm::m2_Eta [private] |
align::ID HIPAlignmentAlgorithm::m2_Id [private] |
int HIPAlignmentAlgorithm::m2_Layer [private] |
int HIPAlignmentAlgorithm::m2_Nhit [private] |
float HIPAlignmentAlgorithm::m2_Phi [private] |
int HIPAlignmentAlgorithm::m2_Type [private] |
float HIPAlignmentAlgorithm::m2_Xpos [private] |
float HIPAlignmentAlgorithm::m2_Ypos [private] |
float HIPAlignmentAlgorithm::m2_Zpos [private] |
float HIPAlignmentAlgorithm::m_Chi2n[MAXREC] [private] |
float HIPAlignmentAlgorithm::m_Eta[MAXREC] [private] |
int HIPAlignmentAlgorithm::m_Nhits[MAXREC] [private] |
int HIPAlignmentAlgorithm::m_Ntracks [private] |
float HIPAlignmentAlgorithm::m_Phi[MAXREC] [private] |
float HIPAlignmentAlgorithm::m_Pt[MAXREC] [private] |
const int HIPAlignmentAlgorithm::MAXREC = 99 [static, private] |
std::string HIPAlignmentAlgorithm::outfile [private] |
Definition at line 65 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::outfile2 [private] |
Definition at line 65 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::outpath [private] |
std::string HIPAlignmentAlgorithm::salignedfile [private] |
Definition at line 66 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), startNewLoop(), and terminate().
std::string HIPAlignmentAlgorithm::siterationfile [private] |
Definition at line 66 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), startNewLoop(), and terminate().
std::string HIPAlignmentAlgorithm::smisalignedfile [private] |
Definition at line 66 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and startNewLoop().
std::string HIPAlignmentAlgorithm::sparameterfile [private] |
Definition at line 65 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and terminate().
std::string HIPAlignmentAlgorithm::struefile [private] |
Definition at line 66 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and startNewLoop().
std::string HIPAlignmentAlgorithm::suvarfile [private] |
Definition at line 65 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and terminate().
std::vector<Alignable*> HIPAlignmentAlgorithm::theAlignables [private] |
Definition at line 53 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), fillRoot(), initialize(), startNewLoop(), and terminate().
Definition at line 52 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), fillRoot(), initialize(), run(), setAlignmentPositionError(), startNewLoop(), and terminate().
std::vector<std::pair<std::vector<Alignable*>, std::vector<double> > > HIPAlignmentAlgorithm::theAPEParameters [private] |
Definition at line 70 of file HIPAlignmentAlgorithm.h.
Referenced by initialize(), and setAlignmentPositionError().
std::vector<edm::ParameterSet> HIPAlignmentAlgorithm::theAPEParameterSet [private] |
Definition at line 69 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and initialize().
bool HIPAlignmentAlgorithm::theApplyAPE [private] |
Definition at line 68 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), initialize(), and setAlignmentPositionError().
int HIPAlignmentAlgorithm::theCollectorNJobs [private] |
Definition at line 79 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::theCollectorPath [private] |
Definition at line 80 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), and HIPAlignmentAlgorithm().
int HIPAlignmentAlgorithm::theCurrentPrescale [private] |
Definition at line 81 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and run().
int HIPAlignmentAlgorithm::theEventPrescale [private] |
Definition at line 81 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and run().
TFile* HIPAlignmentAlgorithm::theFile [private] |
Definition at line 86 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), run(), and terminate().
TFile* HIPAlignmentAlgorithm::theFile2 [private] |
Definition at line 88 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillRoot(), and terminate().
AlignmentIORoot HIPAlignmentAlgorithm::theIO [private] |
Definition at line 56 of file HIPAlignmentAlgorithm.h.
Referenced by startNewLoop(), and terminate().
int HIPAlignmentAlgorithm::theIteration [private] |
Definition at line 58 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), collector(), setAlignmentPositionError(), startNewLoop(), and terminate().
std::vector<align::StructureType> HIPAlignmentAlgorithm::theLevels [private] |
Definition at line 83 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and startNewLoop().
double HIPAlignmentAlgorithm::theMaxAllowedHitPull [private] |
Definition at line 72 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and run().
double HIPAlignmentAlgorithm::theMaxRelParameterError [private] |
Definition at line 76 of file HIPAlignmentAlgorithm.h.
Referenced by calcParameters(), and HIPAlignmentAlgorithm().
Definition at line 74 of file HIPAlignmentAlgorithm.h.
Referenced by calcParameters(), and HIPAlignmentAlgorithm().
TTree* HIPAlignmentAlgorithm::theTree [private] |
Definition at line 87 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), run(), and terminate().
TTree* HIPAlignmentAlgorithm::theTree2 [private] |
Definition at line 89 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillRoot(), and terminate().
bool HIPAlignmentAlgorithm::verbose [private] |
Definition at line 63 of file HIPAlignmentAlgorithm.h.
1.5.4