|
|
|
#include <HIPAlignmentAlgorithm.h>
Definition at line 24 of file HIPAlignmentAlgorithm.h.
| HIPAlignmentAlgorithm::HIPAlignmentAlgorithm | ( | const edm::ParameterSet & | cfg | ) |
Constructor.
Definition at line 33 of file HIPAlignmentAlgorithm.cc.
References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), isCollector, prof2calltree::l, AlignableObjectId::nameToType(), outfile, outfile2, outpath, salignedfile, siterationfile, smisalignedfile, sparameterfile, ssurveyfile, struefile, suvarfile, theAPEParameterSet, theApplyAPE, theCollectorNJobs, theCollectorPath, theCurrentPrescale, theEventPrescale, theFillTrackMonitoring, theLevels, theMaxAllowedHitPull, theMaxRelParameterError, and theMinimumNumberOfHits.
: AlignmentAlgorithmBase( cfg ) { // parse parameters verbose = cfg.getParameter<bool>("verbosity"); outpath = cfg.getParameter<string>("outpath"); outfile = cfg.getParameter<string>("outfile"); outfile2 = cfg.getParameter<string>("outfile2"); struefile = cfg.getParameter<string>("trueFile"); smisalignedfile = cfg.getParameter<string>("misalignedFile"); salignedfile = cfg.getParameter<string>("alignedFile"); siterationfile = cfg.getParameter<string>("iterationFile"); suvarfile = cfg.getParameter<string>("uvarFile"); sparameterfile = cfg.getParameter<string>("parameterFile"); ssurveyfile = cfg.getParameter<string>("surveyFile"); outfile =outpath+outfile;//Eventwise tree outfile2 =outpath+outfile2;//Alignablewise tree struefile =outpath+struefile; smisalignedfile=outpath+smisalignedfile; salignedfile =outpath+salignedfile; siterationfile =outpath+siterationfile; suvarfile =outpath+suvarfile; sparameterfile =outpath+sparameterfile; ssurveyfile =outpath+ssurveyfile; // parameters for APE theApplyAPE = cfg.getParameter<bool>("applyAPE"); theAPEParameterSet = cfg.getParameter<std::vector<edm::ParameterSet> >("apeParam"); theMaxAllowedHitPull = cfg.getParameter<double>("maxAllowedHitPull"); theMinimumNumberOfHits = cfg.getParameter<int>("minimumNumberOfHits"); theMaxRelParameterError = cfg.getParameter<double>("maxRelParameterError"); // for collector mode (parallel processing) isCollector=cfg.getParameter<bool>("collectorActive"); theCollectorNJobs=cfg.getParameter<int>("collectorNJobs"); theCollectorPath=cfg.getParameter<string>("collectorPath"); theFillTrackMonitoring=cfg.getUntrackedParameter<bool>("fillTrackMonitoring"); if (isCollector) edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Collector mode"; theEventPrescale = cfg.getParameter<int>("eventPrescale"); theCurrentPrescale = theEventPrescale; AlignableObjectId dummy; const std::vector<std::string>& levels = cfg.getUntrackedParameter<std::vector<std::string> >("surveyResiduals"); for (unsigned int l = 0; l < levels.size(); ++l) { theLevels.push_back( dummy.nameToType(levels[l]) ); } edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] constructed."; }
| HIPAlignmentAlgorithm::~HIPAlignmentAlgorithm | ( | ) | [inline] |
| void HIPAlignmentAlgorithm::bookRoot | ( | void | ) | [private] |
Definition at line 1030 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, m3_Id, m3_ObjId, m3_par, m_Chi2n, m_d0, m_dz, m_Eta, m_Nhits, m_nhPXB, m_nhPXF, m_Ntracks, m_P, m_Phi, m_Pt, outfile, outfile2, ssurveyfile, theFile, theFile2, theFile3, theIteration, theLevels, theTree, theTree2, and theTree3.
Referenced by startNewLoop().
{
// create ROOT files
theFile = new TFile(outfile.c_str(),"update");
theFile->cd();
// book event-wise ROOT Tree
TString tname="T1";
char iterString[5];
sprintf(iterString, "%i",theIteration);
tname.Append("_");
tname.Append(iterString);
theTree = new TTree(tname,"Eventwise tree");
//theTree->Branch("Run", &m_Run, "Run/I");
//theTree->Branch("Event", &m_Event, "Event/I");
theTree->Branch("Ntracks", &m_Ntracks, "Ntracks/I");
theTree->Branch("Nhits", m_Nhits, "Nhits[Ntracks]/I");
theTree->Branch("nhPXB", m_nhPXB, "nhPXB[Ntracks]/I");
theTree->Branch("nhPXF", m_nhPXF, "nhPXF[Ntracks]/I");
theTree->Branch("Pt", m_Pt, "Pt[Ntracks]/F");
theTree->Branch("P", m_P, "P[Ntracks]/F");
theTree->Branch("Eta", m_Eta, "Eta[Ntracks]/F");
theTree->Branch("Phi", m_Phi, "Phi[Ntracks]/F");
theTree->Branch("Chi2n", m_Chi2n, "Chi2n[Ntracks]/F");
theTree->Branch("d0", m_d0, "d0[Ntracks]/F");
theTree->Branch("dz", m_dz, "dz[Ntracks]/F");
// book Alignable-wise ROOT Tree
theFile2 = new TFile(outfile2.c_str(),"update");
theFile2->cd();
theTree2 = new TTree("T2","Alignablewise tree");
theTree2->Branch("Nhit", &m2_Nhit, "Nhit/I");
theTree2->Branch("Type", &m2_Type, "Type/I");
theTree2->Branch("Layer", &m2_Layer, "Layer/I");
theTree2->Branch("Xpos", &m2_Xpos, "Xpos/F");
theTree2->Branch("Ypos", &m2_Ypos, "Ypos/F");
theTree2->Branch("Zpos", &m2_Zpos, "Zpos/F");
theTree2->Branch("Eta", &m2_Eta, "Eta/F");
theTree2->Branch("Phi", &m2_Phi, "Phi/F");
theTree2->Branch("Id", &m2_Id, "Id/i");
theTree2->Branch("ObjId", &m2_ObjId, "ObjId/I");
// book survey-wise ROOT Tree only if survey is enabled
if (theLevels.size() > 0){
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Survey trees booked.";
theFile3 = new TFile(ssurveyfile.c_str(),"update");
theFile3->cd();
theTree3 = new TTree(tname, "Survey Tree");
theTree3->Branch("Id", &m3_Id, "Id/i");
theTree3->Branch("ObjId", &m3_ObjId, "ObjId/I");
theTree3->Branch("Par", &m3_par, "Par[6]/F");
}
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Root trees booked.";
}
| double HIPAlignmentAlgorithm::calcAPE | ( | double * | par, |
| int | iter, | ||
| double | function | ||
| ) | [private] |
Definition at line 1002 of file HIPAlignmentAlgorithm.cc.
References funct::exp(), max(), funct::pow(), and launcher::step.
Referenced by setAlignmentPositionError().
{
double diter=(double)iter;
// The following floating-point equality check is safe because this
// "0." and this "1." are generated by the compiler, in the very
// same file. Whatever approximization scheme it uses to turn "1."
// into 0.9999999999998 in the HIPAlignmentAlgorithm::initialize is
// also used here. If I'm wrong, you'll get an assertion.
if (function == 0.) {
return max(par[1],par[0]+((par[1]-par[0])/par[2])*diter);
}
else if (function == 1.) {
return max(0.,par[0]*(exp(-pow(diter,par[1])/par[2])));
}
else if (function == 2.) {
int ipar2 = (int) par[2];
int step = iter/ipar2;
double dstep = (double) step;
return max(0.0, par[0] - par[1]*dstep);
}
else assert(false); // should have been caught in the constructor
}
| bool HIPAlignmentAlgorithm::calcParameters | ( | Alignable * | ali | ) | [private] |
Definition at line 1172 of file HIPAlignmentAlgorithm.cc.
References abs, Alignable::alignmentParameters(), AlignmentParameters::cloneFromSelected(), i, HIPUserVariables::jtve, HIPUserVariables::jtvj, HIPUserVariables::nhit, Alignable::setAlignmentParameters(), AlignmentParameters::setValid(), mathSSE::sqrt(), theMaxRelParameterError, theMinimumNumberOfHits, and AlignmentParameters::userVariables().
Referenced by terminate().
{
// Alignment parameters
AlignmentParameters* par = ali->alignmentParameters();
// access user variables
HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(par->userVariables());
int nhit = uservar->nhit;
// The following variable is needed for the extended 1D/2D hit fix using
// matrix shrinkage and expansion
// int hitdim = uservar->hitdim;
if (nhit < theMinimumNumberOfHits) {
par->setValid(false);
return false;
}
AlgebraicSymMatrix jtvj = uservar->jtvj;
AlgebraicVector jtve = uservar->jtve;
// Shrink input in case of 1D hits and 'v' selected
// in alignment parameters
// if (hitdim==1 && selector[1]==true) {
// int iremove = 1;
// if (selector[0]==false) iremove--;
// AlgebraicSymMatrix tempjtvj(jtvj.num_row()-1);
// int nr = 0, nc = 0;
// for (int r=0;r<jtvj.num_row();r++) {
// if (r==iremove) continue;
// nc = 0;
// for (int c=0;c<jtvj.num_col();c++) {
// if (c==iremove) continue;
// tempjtvj[nr][nc] = jtvj[r][c];
// nc++;
// }
// nr++;
// }
// jtvj = tempjtvj;
// AlgebraicVector tempjtve(jtve.num_row()-1);
// nr = 0;
// for (int r=0;r<jtve.num_row();r++) {
// if (r==iremove) continue;
// tempjtve[nr] = jtve[r];
// nr++;
// }
// jtve = tempjtve;
// }
int ierr;
AlgebraicSymMatrix jtvjinv = jtvj.inverse(ierr);
if (ierr !=0) {
edm::LogError("Alignment") << "Matrix inversion failed!";
return false;
}
// these are the alignment corrections+covariance (for selected params)
AlgebraicVector params = - (jtvjinv * jtve);
AlgebraicSymMatrix cov = jtvjinv;
edm::LogInfo("Alignment") << "parameters " << params;
// errors of parameters
int npar = params.num_row();
AlgebraicVector paramerr(npar);
AlgebraicVector relerr(npar);
for (int i=0;i<npar;i++) {
if (abs(cov[i][i])>0) paramerr[i] = sqrt(abs(cov[i][i]));
else paramerr[i] = params[i];
relerr[i] = abs(paramerr[i]/params[i]);
if (relerr[i] >= theMaxRelParameterError) {
params[i] = 0;
paramerr[i]=0;
}
}
// expand output in case of 1D hits and 'v' selected
// in alignment parameters
// if (hitdim==1 && selector[1]==true) {
// int iremove = 1;
// if (selector[0]==false) iremove--;
// AlgebraicSymMatrix tempcov(cov.num_row()+1,0);
// int nr = 0, nc = 0;
// for (int r=0;r<cov.num_row();r++) {
// if (r==iremove) nr++;
// nc = 0;
// for (int c=0;c<cov.num_col();c++) {
// if (c==iremove) nc++;
// tempcov[nr][nc] = cov[r][c];
// nc++;
// }
// nr++;
// }
// cov = tempcov;
// AlgebraicVector tempparams(params.num_row()+1,0);
// nr = 0;
// for (int r=0;r<params.num_row();r++) {
// if (r==iremove) nr++;
// tempparams[nr] = params[r];
// nr++;
// }
// params = tempparams;
// }
// store alignment parameters
AlignmentParameters* parnew = par->cloneFromSelected(params,cov);
ali->setAlignmentParameters(parnew);
parnew->setValid(true);
return true;
}
| void HIPAlignmentAlgorithm::collector | ( | void | ) | [private] |
Definition at line 1289 of file HIPAlignmentAlgorithm.cc.
References HIPUserVariables::alichi2, Alignable::alignmentParameters(), HIPUserVariables::alindof, AlignmentParameterStore::attachUserVariables(), HIPUserVariables::clone(), fillEventwiseTree(), ioerr, HIPUserVariables::jtve, HIPUserVariables::jtvj, HIPUserVariables::nhit, HIPUserVariablesIORoot::readHIPUserVariables(), theAlignables, theAlignmentParameterStore, theCollectorNJobs, theCollectorPath, theFillTrackMonitoring, theIteration, and AlignmentParameters::userVariables().
Referenced by startNewLoop().
{
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::collector] called for iteration "
<< theIteration << std::endl;
HIPUserVariablesIORoot HIPIO;
for (int ijob=1;ijob<=theCollectorNJobs;ijob++) {
edm::LogWarning("Alignment") << "reading uservar for job " << ijob;
stringstream ss;
string str;
ss << ijob;
ss >> str;
string uvfile = theCollectorPath+"/job"+str+"/IOUserVariables.root";
vector<AlignmentUserVariables*> uvarvec =
HIPIO.readHIPUserVariables(theAlignables, (char*)uvfile.c_str(),
theIteration, ioerr);
if (ioerr!=0) {
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::collector] could not read user variable files for job "
<< ijob;
continue;
}
// add
vector<AlignmentUserVariables*> uvarvecadd;
vector<AlignmentUserVariables*>::const_iterator iuvarnew=uvarvec.begin();
for (vector<Alignable*>::const_iterator it=theAlignables.begin();
it!=theAlignables.end();
++it) {
Alignable* ali = *it;
AlignmentParameters* ap = ali->alignmentParameters();
HIPUserVariables* uvarold = dynamic_cast<HIPUserVariables*>(ap->userVariables());
HIPUserVariables* uvarnew = dynamic_cast<HIPUserVariables*>(*iuvarnew);
HIPUserVariables* uvar = uvarold->clone();
if (uvarnew!=0) {
uvar->nhit = (uvarold->nhit)+(uvarnew->nhit);
uvar->jtvj = (uvarold->jtvj)+(uvarnew->jtvj);
uvar->jtve = (uvarold->jtve)+(uvarnew->jtve);
uvar->alichi2 = (uvarold->alichi2)+(uvarnew->alichi2);
uvar->alindof = (uvarold->alindof)+(uvarnew->alindof);
delete uvarnew;
}
uvarvecadd.push_back(uvar);
iuvarnew++;
}
theAlignmentParameterStore->attachUserVariables(theAlignables, uvarvecadd, ioerr);
//fill Eventwise Tree
if (theFillTrackMonitoring) {
uvfile = theCollectorPath+"/job"+str+"/HIPAlignmentEvents.root";
edm::LogWarning("Alignment") << "Added to the tree "
<< fillEventwiseTree(uvfile.c_str(), theIteration, ioerr)
<< "tracks";
}
}//end loop on jobs
}
| int HIPAlignmentAlgorithm::fillEventwiseTree | ( | const char * | filename, |
| int | iter, | ||
| int | ierr | ||
| ) | [private] |
Definition at line 1357 of file HIPAlignmentAlgorithm.cc.
References m_Chi2n, m_d0, m_dz, m_Eta, m_Nhits, m_nhPXB, m_nhPXF, m_Ntracks, m_P, m_Phi, m_Pt, MAXREC, and theTree.
Referenced by collector().
{
int totntrk = 0;
char treeName[64];
sprintf(treeName, "T1_%d", iter);
//open the file "HIPAlignmentEvents.root" in the job directory
TFile *jobfile = new TFile(filename, "READ");
//grab the tree corresponding to this iteration
TTree *jobtree = (TTree*)jobfile->Get(treeName);
//address and read the variables
static const int nmaxtrackperevent = 1000;
int jobNtracks, jobNhitspertrack[nmaxtrackperevent], jobnhPXB[nmaxtrackperevent], jobnhPXF[nmaxtrackperevent];
float jobP[nmaxtrackperevent], jobPt[nmaxtrackperevent], jobEta[nmaxtrackperevent] , jobPhi[nmaxtrackperevent];
float jobd0[nmaxtrackperevent], jobdz[nmaxtrackperevent] , jobChi2n[nmaxtrackperevent];
jobtree->SetBranchAddress("Ntracks", &jobNtracks);
jobtree->SetBranchAddress("Nhits", jobNhitspertrack);
jobtree->SetBranchAddress("nhPXB", jobnhPXB);
jobtree->SetBranchAddress("nhPXF", jobnhPXF);
jobtree->SetBranchAddress("Pt", jobPt);
jobtree->SetBranchAddress("P", jobP);
jobtree->SetBranchAddress("d0", jobd0);
jobtree->SetBranchAddress("dz", jobdz);
jobtree->SetBranchAddress("Eta", jobEta);
jobtree->SetBranchAddress("Phi", jobPhi);
jobtree->SetBranchAddress("Chi2n", jobChi2n);
int ievent = 0;
for (ievent=0;ievent<jobtree->GetEntries();++ievent) {
jobtree->GetEntry(ievent);
//fill the collector tree with them
// TO BE IMPLEMENTED: a prescale factor like in run()
m_Ntracks = jobNtracks;
int ntrk = 0;
while (ntrk<m_Ntracks) {
if (ntrk<MAXREC) {
totntrk = ntrk+1;
m_Nhits[ntrk] = jobNhitspertrack[ntrk];
m_Pt[ntrk] = jobPt[ntrk];
m_P[ntrk] = jobP[ntrk];
m_nhPXB[ntrk] = jobnhPXB[ntrk];
m_nhPXF[ntrk] = jobnhPXF[ntrk];
m_Eta[ntrk] = jobEta[ntrk];
m_Phi[ntrk] = jobPhi[ntrk];
m_Chi2n[ntrk] = jobChi2n[ntrk];
m_d0[ntrk] = jobd0[ntrk];
m_dz[ntrk] = jobdz[ntrk];
}//end if j<MAXREC
else{
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::fillEventwiseTree] Number of tracks in Eventwise tree exceeds MAXREC: "
<< m_Ntracks << " Skipping exceeding tracks.";
ntrk = m_Ntracks+1;
}
++ntrk;
}//end while loop
theTree->Fill();
}//end loop on i - entries in the job tree
//clean up
delete jobtree;
delete jobfile;
return totntrk;
}//end fillEventwiseTree
| void HIPAlignmentAlgorithm::fillRoot | ( | void | ) | [private] |
Definition at line 1097 of file HIPAlignmentAlgorithm.cc.
References Alignable::alignableObjectId(), Alignable::alignmentParameters(), PV3DBase< T, PVType, FrameType >::eta(), Alignable::id(), AlignmentParameters::isValid(), m2_Eta, m2_Id, m2_Layer, m2_Nhit, m2_ObjId, m2_Phi, m2_Type, m2_Xpos, m2_Ypos, m2_Zpos, HIPUserVariables::nhit, AlignmentParameters::parameters(), PV3DBase< T, PVType, FrameType >::phi(), pos, 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().
{
theFile2->cd();
int naligned=0;
for (vector<Alignable*>::const_iterator it=theAlignables.begin();
it!=theAlignables.end();
++it) {
Alignable* ali = (*it);
AlignmentParameters* dap = ali->alignmentParameters();
// consider only those parameters classified as 'valid'
if (dap->isValid()) {
// get number of hits from user variable
HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(dap->userVariables());
m2_Nhit = uservar->nhit;
// get type/layer
std::pair<int,int> tl = theAlignmentParameterStore->typeAndLayer(ali);
m2_Type = tl.first;
m2_Layer = tl.second;
// get identifier (as for IO)
m2_Id = ali->id();
m2_ObjId = ali->alignableObjectId();
// get position
GlobalPoint pos = ali->surface().position();
m2_Xpos = pos.x();
m2_Ypos = pos.y();
m2_Zpos = pos.z();
m2_Eta = pos.eta();
m2_Phi = pos.phi();
AlgebraicVector pars = dap->parameters();
if (verbose) {
edm::LogVerbatim("Alignment")
<< "------------------------------------------------------------------------\n"
<< " ALIGNABLE: " << setw(6) << naligned
<< '\n'
<< "hits: " << setw(4) << m2_Nhit
<< " type: " << setw(4) << m2_Type
<< " layer: " << setw(4) << m2_Layer
<< " id: " << setw(4) << m2_Id
<< " objId: " << setw(4) << m2_ObjId
<< '\n'
<< fixed << setprecision(5)
<< "x,y,z: "
<< setw(12) << m2_Xpos
<< setw(12) << m2_Ypos
<< setw(12) << m2_Zpos
<< " eta,phi: "
<< setw(12) << m2_Eta
<< setw(12) << m2_Phi
<< '\n'
<< "params: "
<< setw(12) << pars[0]
<< setw(12) << pars[1]
<< setw(12) << pars[2]
<< setw(12) << pars[3]
<< setw(12) << pars[4]
<< setw(12) << pars[5];
}
naligned++;
theTree2->Fill();
}
}
}
| void HIPAlignmentAlgorithm::initialize | ( | const edm::EventSetup & | setup, |
| AlignableTracker * | tracker, | ||
| AlignableMuon * | muon, | ||
| AlignableExtras * | extras, | ||
| AlignmentParameterStore * | store | ||
| ) | [virtual] |
Call at beginning of job.
Implements AlignmentAlgorithmBase.
Definition at line 96 of file HIPAlignmentAlgorithm.cc.
References AlignmentParameterSelector::addSelections(), AlignmentParameterStore::alignables(), MuonAlignmentFromReference_cff::alignParams, AlignmentParameterSelector::clear(), Exception, edm::ParameterSet::getParameter(), i, AlignmentParameterSelector::selectedAlignables(), theAlignableDetAccessor, theAlignables, theAlignmentParameterStore, theAPEParameters, theAPEParameterSet, and theApplyAPE.
{
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Initializing...";
// accessor Det->AlignableDet
if ( !muon )
theAlignableDetAccessor = new AlignableNavigator(tracker);
else if ( !tracker )
theAlignableDetAccessor = new AlignableNavigator(muon);
else
theAlignableDetAccessor = new AlignableNavigator(tracker, muon);
// set alignmentParameterStore
theAlignmentParameterStore=store;
// get alignables
theAlignables = theAlignmentParameterStore->alignables();
// clear theAPEParameters, if necessary
theAPEParameters.clear();
// get APE parameters
if(theApplyAPE){
AlignmentParameterSelector selector(tracker, muon);
for (std::vector<edm::ParameterSet>::const_iterator setiter = theAPEParameterSet.begin();
setiter != theAPEParameterSet.end();
++setiter) {
std::vector<Alignable*> alignables;
selector.clear();
edm::ParameterSet selectorPSet = setiter->getParameter<edm::ParameterSet>("Selector");
std::vector<std::string> alignParams = selectorPSet.getParameter<std::vector<std::string> >("alignParams");
if (alignParams.size() == 1 && alignParams[0] == std::string("selected")) {
alignables = theAlignables;
}
else {
selector.addSelections(selectorPSet);
alignables = selector.selectedAlignables();
}
std::vector<double> apeSPar = setiter->getParameter<std::vector<double> >("apeSPar");
std::vector<double> apeRPar = setiter->getParameter<std::vector<double> >("apeRPar");
std::string function = setiter->getParameter<std::string>("function");
if (apeSPar.size() != 3 || apeRPar.size() != 3)
throw cms::Exception("BadConfig") << "apeSPar and apeRPar must have 3 values each" << std::endl;
for (std::vector<double>::const_iterator i = apeRPar.begin(); i != apeRPar.end(); ++i) {
apeSPar.push_back(*i);
}
if (function == std::string("linear")) {
apeSPar.push_back(0.); // c.f. note in calcAPE
}
else if (function == std::string("exponential")) {
apeSPar.push_back(1.); // c.f. note in calcAPE
}
else if (function == std::string("step")) {
apeSPar.push_back(2.); // c.f. note in calcAPE
}
else {
throw cms::Exception("BadConfig") << "APE function must be \"linear\" or \"exponential\"." << std::endl;
}
theAPEParameters.push_back(std::pair<std::vector<Alignable*>, std::vector<double> >(alignables, apeSPar));
}
}
}
| bool HIPAlignmentAlgorithm::processHit1D | ( | const AlignableDetOrUnitPtr | alidet, |
| const Alignable * | ali, | ||
| const TrajectoryStateOnSurface & | tsos, | ||
| const TransientTrackingRecHit * | hit | ||
| ) | [private] |
Definition at line 434 of file HIPAlignmentAlgorithm.cc.
References HIPUserVariables::alichi2, Alignable::alignmentParameters(), HIPUserVariables::alindof, HIPUserVariables::jtve, HIPUserVariables::jtvj, TrajectoryStateOnSurface::localError(), TrackingRecHit::localPosition(), TrajectoryStateOnSurface::localPosition(), TrackingRecHit::localPositionError(), HIPUserVariables::nhit, pos, LocalTrajectoryError::positionError(), AlignmentParameters::selectedDerivatives(), mathSSE::sqrt(), theMaxAllowedHitPull, AlignmentParameters::userVariables(), PV3DBase< T, PVType, FrameType >::x(), and LocalError::xx().
Referenced by run().
{
static const unsigned int hitDim = 1;
// get trajectory impact point
LocalPoint alvec = tsos.localPosition();
AlgebraicVector pos(hitDim);
pos[0] = alvec.x();
// get impact point covariance
AlgebraicSymMatrix ipcovmat(hitDim);
ipcovmat[0][0] = tsos.localError().positionError().xx();
// get hit local position and covariance
AlgebraicVector coor(hitDim);
coor[0] = hit->localPosition().x();
AlgebraicSymMatrix covmat(hitDim);
covmat[0][0] = hit->localPositionError().xx();
// add hit and impact point covariance matrices
covmat = covmat + ipcovmat;
// calculate the x pull of this hit
double xpull = 0.;
if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0]));
// get Alignment Parameters
AlignmentParameters* params = ali->alignmentParameters();
// get derivatives
AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
// calculate user parameters
int npar = derivs2D.num_row();
// std::cout << "Dumping derivsSEL: \n" << derivs2D <<std::endl;
AlgebraicMatrix derivs(npar, hitDim, 0);
for (int ipar=0;ipar<npar;ipar++) {
for (unsigned int idim=0;idim<hitDim;idim++) {
derivs[ipar][idim] = derivs2D[ipar][idim];
}
}
// std::cout << "Dumping derivs: \n" << derivs << std::endl;
// invert covariance matrix
int ierr;
// if (covmat[1][1]>4.0) nhitDim = hitDim;
covmat.invert(ierr);
if (ierr != 0) {
edm::LogError("Alignment") << "Matrix inversion failed!";
return false;
}
bool useThisHit = (theMaxAllowedHitPull <= 0.);
useThisHit |= (fabs(xpull) < theMaxAllowedHitPull);
// bailing out
if (!useThisHit) return false;
// std::cout << "We use this hit in " << subDet << std::endl;
// std::cout << "Npars= " << npar << " NhitDim=1 Size of derivs: "
// << derivs.num_row() << " x " << derivs.num_col() << std::endl;
// AlgebraicMatrix jtvjtmp();
AlgebraicMatrix covtmp(covmat);
// std::cout << "Preapring JTVJTMP -> " << std::flush;
AlgebraicMatrix jtvjtmp(derivs * covtmp *derivs.T());
// std::cout << "TMP JTVJ= \n" << jtvjtmp << std::endl;
AlgebraicSymMatrix thisjtvj(npar);
AlgebraicVector thisjtve(npar);
// std::cout << "Preparing ThisJtVJ" << std::flush;
// thisjtvj = covmat.similarity(derivs);
thisjtvj.assign(jtvjtmp);
// cout<<" ThisJtVE"<<std::endl;
thisjtve=derivs * covmat * (pos-coor);
AlgebraicVector hitresidual(hitDim);
hitresidual[0] = (pos[0] - coor[0]);
AlgebraicMatrix hitresidualT;
hitresidualT = hitresidual.T();
// std::cout << "HitResidualT = \n" << hitresidualT << std::endl;
// access user variables (via AlignmentParameters)
HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
uservar->jtvj += thisjtvj;
uservar->jtve += thisjtve;
uservar->nhit++;
// The following variable is needed for the extended 1D/2D hit fix using
// matrix shrinkage and expansion
// uservar->hitdim = hitDim;
//for alignable chi squared
float thischi2;
thischi2 = (hitresidualT *covmat *hitresidual)[0];
if ( verbose && ((thischi2/ static_cast<float>(uservar->nhit)) >10.0) ) {
edm::LogWarning("Alignment") << "Added to Chi2 the number " << thischi2 <<" having "
<< uservar->nhit << " dof " << std::endl << "X-resid "
<< hitresidual[0] << " Y-resid "
<< hitresidual[1] << std::endl << " Cov^-1 matr (covmat): [0][0]= "
<< covmat[0][0] << " [0][1]= "
<< covmat[0][1] << " [1][0]= "
<< covmat[1][0] << " [1][1]= "
<< covmat[1][1] << std::endl;
}
uservar->alichi2 += thischi2; // a bit weird, but vector.transposed * matrix * vector doesn't give a double in CMSSW's opinion
uservar->alindof += hitDim; // 2D hits contribute twice to the ndofs
return true;
}
| bool HIPAlignmentAlgorithm::processHit2D | ( | const AlignableDetOrUnitPtr | alidet, |
| const Alignable * | ali, | ||
| const TrajectoryStateOnSurface & | tsos, | ||
| const TransientTrackingRecHit * | hit | ||
| ) | [private] |
Definition at line 550 of file HIPAlignmentAlgorithm.cc.
References HIPUserVariables::alichi2, Alignable::alignmentParameters(), HIPUserVariables::alindof, HIPUserVariables::jtve, HIPUserVariables::jtvj, TrajectoryStateOnSurface::localError(), TrackingRecHit::localPosition(), TrajectoryStateOnSurface::localPosition(), TrackingRecHit::localPositionError(), HIPUserVariables::nhit, pos, LocalTrajectoryError::positionError(), AlignmentParameters::selectedDerivatives(), mathSSE::sqrt(), theMaxAllowedHitPull, AlignmentParameters::userVariables(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), LocalError::xy(), PV3DBase< T, PVType, FrameType >::y(), and LocalError::yy().
Referenced by run().
{
static const unsigned int hitDim = 2;
// get trajectory impact point
LocalPoint alvec = tsos.localPosition();
AlgebraicVector pos(hitDim);
pos[0] = alvec.x();
pos[1] = alvec.y();
// get impact point covariance
AlgebraicSymMatrix ipcovmat(hitDim);
ipcovmat[0][0] = tsos.localError().positionError().xx();
ipcovmat[1][1] = tsos.localError().positionError().yy();
ipcovmat[0][1] = tsos.localError().positionError().xy();
// get hit local position and covariance
AlgebraicVector coor(hitDim);
coor[0] = hit->localPosition().x();
coor[1] = hit->localPosition().y();
AlgebraicSymMatrix covmat(hitDim);
covmat[0][0] = hit->localPositionError().xx();
covmat[1][1] = hit->localPositionError().yy();
covmat[0][1] = hit->localPositionError().xy();
// add hit and impact point covariance matrices
covmat = covmat + ipcovmat;
// calculate the x pull and y pull of this hit
double xpull = 0.;
double ypull = 0.;
if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0]));
if (covmat[1][1] != 0.) ypull = (pos[1] - coor[1])/sqrt(fabs(covmat[1][1]));
// get Alignment Parameters
AlignmentParameters* params = ali->alignmentParameters();
// get derivatives
AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
// calculate user parameters
int npar = derivs2D.num_row();
// std::cout << "Dumping derivsSEL: \n"<< derivs2D <<std::endl;
AlgebraicMatrix derivs(npar, hitDim, 0);
for (int ipar=0;ipar<npar;ipar++) {
for (unsigned int idim=0;idim<hitDim;idim++) {
derivs[ipar][idim] = derivs2D[ipar][idim];
}
}
// std::cout << "Dumping derivs: \n" << derivs << std::endl;
// invert covariance matrix
int ierr;
// if (covmat[1][1]>4.0) nhitDim = hitDim;
covmat.invert(ierr);
if (ierr != 0) {
edm::LogError("Alignment") << "Matrix inversion failed!";
return false;
}
bool useThisHit = (theMaxAllowedHitPull <= 0.);
useThisHit |= (fabs(xpull) < theMaxAllowedHitPull && fabs(ypull) < theMaxAllowedHitPull);
// bailing out
if (!useThisHit) return false;
// std::cout << "We use this hit in " << subDet << std::endl;
// std::cout << "Npars= " << npar << " NhitDim=2 Size of derivs: "
// << derivs.num_row() << " x " << derivs.num_col() << std::endl;
// AlgebraicMatrix jtvjtmp();
AlgebraicMatrix covtmp(covmat);
// std::cout << "Preapring JTVJTMP -> " << std::flush;
AlgebraicMatrix jtvjtmp(derivs * covtmp *derivs.T());
// std::cout << "TMP JTVJ= \n" << jtvjtmp << std::endl;
AlgebraicSymMatrix thisjtvj(npar);
AlgebraicVector thisjtve(npar);
// std::cout << "Preparing ThisJtVJ" << std::flush;
// thisjtvj = covmat.similarity(derivs);
thisjtvj.assign(jtvjtmp);
// cout<<" ThisJtVE"<<std::endl;
thisjtve=derivs * covmat * (pos-coor);
AlgebraicVector hitresidual(hitDim);
hitresidual[0] = (pos[0] - coor[0]);
hitresidual[1] = (pos[1] - coor[1]);
// if(nhitDim>1) {
// hitresidual[1] =0.0;
// nhitDim=1;
// }
AlgebraicMatrix hitresidualT;
hitresidualT = hitresidual.T();
// std::cout << "HitResidualT = \n" << hitresidualT << std::endl;
// access user variables (via AlignmentParameters)
HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
uservar->jtvj += thisjtvj;
uservar->jtve += thisjtve;
uservar->nhit++;
// The following variable is needed for the extended 1D/2D hit fix using
// matrix shrinkage and expansion
// uservar->hitdim = hitDim;
//for alignable chi squared
float thischi2;
thischi2 = (hitresidualT *covmat *hitresidual)[0];
if ( verbose && ((thischi2/ static_cast<float>(uservar->nhit)) >10.0) ) {
edm::LogWarning("Alignment") << "Added to Chi2 the number " << thischi2 <<" having "
<< uservar->nhit << " dof " << std::endl << "X-resid "
<< hitresidual[0] << " Y-resid "
<< hitresidual[1] << std::endl << " Cov^-1 matr (covmat): [0][0]= "
<< covmat[0][0] << " [0][1]= "
<< covmat[0][1] << " [1][0]= "
<< covmat[1][0] << " [1][1]= "
<< covmat[1][1] << std::endl;
}
uservar->alichi2 += thischi2; // a bit weird, but vector.transposed * matrix * vector doesn't give a double in CMSSW's opinion
uservar->alindof += hitDim; // 2D hits contribute twice to the ndofs
return true;
}
| int HIPAlignmentAlgorithm::readIterationFile | ( | std::string | filename | ) | [private] |
Definition at line 910 of file HIPAlignmentAlgorithm.cc.
References recoMuon::in, and query::result.
Referenced by startNewLoop().
{
int result;
ifstream inIterFile((char*)filename.c_str(), ios::in);
if (!inIterFile) {
edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] ERROR! "
<< "Unable to open Iteration file";
result = -1;
}
else {
inIterFile >> result;
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] "
<< "Read last iteration number from file: " << result;
}
inIterFile.close();
return result;
}
| void HIPAlignmentAlgorithm::run | ( | const edm::EventSetup & | setup, |
| const EventInfo & | eventInfo | ||
| ) | [virtual] |
Run the algorithm.
Implements AlignmentAlgorithmBase.
Definition at line 679 of file HIPAlignmentAlgorithm.cc.
References CompositeAlignmentParameters::alignableFromAlignableDet(), AlignableNavigator::alignableFromGeomDet(), AlignableNavigator::alignablesFromHits(), TrajectoryMeasurement::backwardPredictedState(), className(), SiStripRecHit1D::cluster(), SiStripRecHit2D::cluster(), SiPixelRecHit::cluster(), AlignmentAlgorithmBase::EventInfo::clusterValueMap_, TrajectoryStateCombiner::combine(), debug_cff::d0, reco::TrackBase::d0(), AlignableNavigator::detAndSubdetInMap(), reco::TrackBase::dz(), eta(), reco::TrackBase::eta(), Exception, TrajectoryMeasurement::forwardPredictedState(), TrackingRecHit::geographicalId(), TransientTrackingRecHit::hit(), reco::TrackBase::hitPattern(), isCollector, AlignmentClusterFlag::isTaken(), TrackingRecHit::isValid(), TrajectoryStateOnSurface::isValid(), m_Chi2n, m_d0, m_dz, m_Eta, m_Nhits, m_nhPXB, m_nhPXF, m_Ntracks, m_P, m_Phi, m_Pt, MAXREC, Trajectory::measurements(), reco::TrackBase::normalizedChi2(), reco::TrackBase::numberOfValidHits(), reco::TrackBase::p(), AlCaHLTBitMon_ParallelJobs::p, phi, reco::TrackBase::phi(), processHit1D(), processHit2D(), reco::TrackBase::pt(), TrajectoryMeasurement::recHit(), AlignmentParameterStore::selectParameters(), DetId::subdetId(), theAlignableDetAccessor, theAlignmentParameterStore, theCurrentPrescale, theEventPrescale, theFile, theFillTrackMonitoring, theTree, testEve_cfg::tracks, and AlignmentAlgorithmBase::EventInfo::trajTrackPairs_.
{
if (isCollector) return;
TrajectoryStateCombiner tsoscomb;
// AM: not really needed
// AM: m_Ntracks = 0 should be sufficient
int itr=0;
m_Ntracks=0;
for(itr=0;itr<MAXREC;++itr){
m_Nhits[itr]=0;
m_Pt[itr]=-5.0;
m_P[itr]=-5.0;
m_nhPXB[itr]=0;
m_nhPXF[itr]=0;
m_Eta[itr]=-99.0;
m_Phi[itr]=-4.0;
m_Chi2n[itr]=-11.0;
m_d0[itr]=-999;
m_dz[itr]=-999;
}
itr=0;
// AM: what is this needed for?
theFile->cd();
// loop over tracks
const ConstTrajTrackPairCollection &tracks = eventInfo.trajTrackPairs_;
for (ConstTrajTrackPairCollection::const_iterator it=tracks.begin();
it!=tracks.end();
++it) {
const Trajectory* traj = (*it).first;
const reco::Track* track = (*it).second;
float pt = track->pt();
float eta = track->eta();
float phi = track->phi();
float p = track->p();
float chi2n = track->normalizedChi2();
int nhit = track->numberOfValidHits();
float d0 = track->d0();
float dz = track->dz();
int nhpxb = track->hitPattern().numberOfValidPixelBarrelHits();
int nhpxf = track->hitPattern().numberOfValidPixelEndcapHits();
if (verbose) edm::LogInfo("Alignment") << "New track pt,eta,phi,chi2n,hits: "
<< pt << ","
<< eta << ","
<< phi << ","
<< chi2n << ","
<< nhit;
// edm::LogWarning("Alignment") << "New track pt,eta,phi,chi2n,hits: "
// << pt << ","
// << eta << ","
// << phi << ","
// << chi2n << ","
// << nhit;
// fill track parameters in root tree
if (itr<MAXREC) {
m_Nhits[itr]=nhit;
m_Pt[itr]=pt;
m_P[itr]=p;
m_Eta[itr]=eta;
m_Phi[itr]=phi;
m_Chi2n[itr]=chi2n;
m_nhPXB[itr]=nhpxb;
m_nhPXF[itr]=nhpxf;
m_d0[itr]=d0;
m_dz[itr]=dz;
itr++;
m_Ntracks=itr;
}
// AM: Can be simplified
vector<const TransientTrackingRecHit*> hitvec;
vector<TrajectoryStateOnSurface> tsosvec;
// loop over measurements
vector<TrajectoryMeasurement> measurements = traj->measurements();
for (vector<TrajectoryMeasurement>::iterator im=measurements.begin();
im!=measurements.end();
++im) {
TrajectoryMeasurement meas = *im;
// const TransientTrackingRecHit* ttrhit = &(*meas.recHit());
// const TrackingRecHit *hit = ttrhit->hit();
const TransientTrackingRecHit* hit = &(*meas.recHit());
if (hit->isValid() && theAlignableDetAccessor->detAndSubdetInMap( hit->geographicalId() )) {
// this is the updated state (including the current hit)
// TrajectoryStateOnSurface tsos=meas.updatedState();
// combine fwd and bwd predicted state to get state
// which excludes current hit
bool skiphit = false;
if (eventInfo.clusterValueMap_) {
// check from the PrescalingMap if the hit was taken.
// If not skip to the next TM
// bool hitTaken=false;
AlignmentClusterFlag myflag;
int subDet = hit->geographicalId().subdetId();
//take the actual RecHit out of the Transient one
const TrackingRecHit *rechit=hit->hit();
if (subDet>2) { // AM: if possible use enum instead of hard-coded value
const std::type_info &type = typeid(*rechit);
if (type == typeid(SiStripRecHit1D)) {
const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(rechit);
if (stripHit1D) {
SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
// myflag=PrescMap[stripclust];
myflag = (*eventInfo.clusterValueMap_)[stripclust];
} else {
edm::LogError("HIPAlignmentAlgorithm")
<< "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
<< "TypeId of the RecHit: " << className(*hit) <<endl;
}
}//end if type = SiStripRecHit1D
else if(type == typeid(SiStripRecHit2D)){
const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(rechit);
if (stripHit2D) {
SiStripRecHit2D::ClusterRef stripclust(stripHit2D->cluster());
// myflag=PrescMap[stripclust];
myflag = (*eventInfo.clusterValueMap_)[stripclust];
} else {
edm::LogError("HIPAlignmentAlgorithm")
<< "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
// << "TypeId of the TTRH: " << className(*ttrhit) << endl;
<< "TypeId of the TTRH: " << className(*hit) << endl;
}
} //end if type == SiStripRecHit2D
} //end if hit from strips
else {
const SiPixelRecHit* pixelhit= dynamic_cast<const SiPixelRecHit*>(rechit);
if (pixelhit) {
SiPixelClusterRefNew pixelclust(pixelhit->cluster());
// myflag=PrescMap[pixelclust];
myflag = (*eventInfo.clusterValueMap_)[pixelclust];
}
else {
edm::LogError("HIPAlignmentAlgorithm")
<< "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Pixel RecHit failed! "
// << "TypeId of the TTRH: " << className(*ttrhit) << endl;
<< "TypeId of the TTRH: " << className(*hit) << endl;
}
} //end 'else' it is a pixel hit
// bool hitTaken=myflag.isTaken();
if (!myflag.isTaken()) {
skiphit=true;
//cout<<"Hit from subdet "<<rechit->geographicalId().subdetId()<<" prescaled out."<<endl;
continue;
}
}//end if Prescaled Hits
if (skiphit) {
throw cms::Exception("LogicError")
<< "ERROR in <HIPAlignmentAlgorithm::run>: this hit should have been skipped!"
<< endl;
}
TrajectoryStateOnSurface tsos = tsoscomb.combine(meas.forwardPredictedState(),
meas.backwardPredictedState());
if(tsos.isValid()){
// hitvec.push_back(ttrhit);
hitvec.push_back(hit);
tsosvec.push_back(tsos);
}
} //hit valid
}
// transform RecHit vector to AlignableDet vector
vector <AlignableDetOrUnitPtr> alidetvec = theAlignableDetAccessor->alignablesFromHits(hitvec);
// get concatenated alignment parameters for list of alignables
CompositeAlignmentParameters aap = theAlignmentParameterStore->selectParameters(alidetvec);
vector<TrajectoryStateOnSurface>::const_iterator itsos=tsosvec.begin();
vector<const TransientTrackingRecHit*>::const_iterator ihit=hitvec.begin();
// loop over vectors(hit,tsos)
while (itsos != tsosvec.end()) {
// get AlignableDet for this hit
const GeomDet* det = (*ihit)->det();
// int subDet= (*ihit)->geographicalId().subdetId();
uint32_t nhitDim = (*ihit)->dimension();
AlignableDetOrUnitPtr alidet = theAlignableDetAccessor->alignableFromGeomDet(det);
// get relevant Alignable
Alignable* ali = aap.alignableFromAlignableDet(alidet);
if (ali!=0) {
bool processHitReturnValue;
if (nhitDim==1) {
processHitReturnValue = processHit1D(alidet, ali, *itsos, *ihit);
} else if (nhitDim==2) {
processHitReturnValue = processHit2D(alidet, ali, *itsos, *ihit);
}
}
itsos++;
ihit++;
}
} // end of track loop
// fill eventwise root tree (with prescale defined in pset)
if (theFillTrackMonitoring) {
theCurrentPrescale--;
if (theCurrentPrescale<=0) {
theTree->Fill();
theCurrentPrescale = theEventPrescale;
}
}
}
| void HIPAlignmentAlgorithm::setAlignmentPositionError | ( | void | ) | [private] |
Definition at line 949 of file HIPAlignmentAlgorithm.cc.
References calcAPE(), i, AlignmentParameterStore::setAlignmentPositionError(), theAlignmentParameterStore, theAPEParameters, theApplyAPE, and theIteration.
Referenced by startNewLoop().
{
// Check if user wants to override APE
if ( !theApplyAPE )
{
edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] No APE applied";
return; // NO APE APPLIED
}
edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] Apply APE!";
double apeSPar[3], apeRPar[3];
for (std::vector<std::pair<std::vector<Alignable*>, std::vector<double> > >::const_iterator alipars = theAPEParameters.begin();
alipars != theAPEParameters.end();
++alipars) {
const std::vector<Alignable*> &alignables = alipars->first;
const std::vector<double> &pars = alipars->second;
apeSPar[0] = pars[0];
apeSPar[1] = pars[1];
apeSPar[2] = pars[2];
apeRPar[0] = pars[3];
apeRPar[1] = pars[4];
apeRPar[2] = pars[5];
double function = pars[6];
// Printout for debug
printf("APE: %u alignables\n", (unsigned int)alignables.size());
for ( int i=0; i<21; ++i ) {
double apelinstest=calcAPE(apeSPar,i,0.);
double apeexpstest=calcAPE(apeSPar,i,1.);
double apelinrtest=calcAPE(apeRPar,i,0.);
double apeexprtest=calcAPE(apeRPar,i,1.);
printf("APE: iter slin sexp rlin rexp: %5d %12.5f %12.5f %12.5f %12.5f\n",
i,apelinstest,apeexpstest,apelinrtest,apeexprtest);
}
// set APE
double apeshift=calcAPE(apeSPar,theIteration,function);
double aperot =calcAPE(apeRPar,theIteration,function);
theAlignmentParameterStore->setAlignmentPositionError( alignables, apeshift, aperot );
}
}
| void HIPAlignmentAlgorithm::startNewLoop | ( | void | ) | [virtual] |
Called at start of new loop.
Reimplemented from AlignmentAlgorithmBase.
Definition at line 168 of file HIPAlignmentAlgorithm.cc.
References AlignmentParameterStore::applyAlignableAbsolutePositions(), bookRoot(), collector(), ioerr, isCollector, AlignmentParameters::numSelected(), AlignmentIORoot::readAlignableAbsolutePositions(), readIterationFile(), salignedfile, setAlignmentPositionError(), AlignmentParameters::setUserVariables(), siterationfile, smisalignedfile, struefile, theAlignables, theAlignmentParameterStore, theIO, theIteration, AlignmentIORoot::writeAlignableAbsolutePositions(), and AlignmentIORoot::writeAlignableOriginalPositions().
{
// iterate over all alignables and attach user variables
for( vector<Alignable*>::const_iterator it=theAlignables.begin();
it!=theAlignables.end(); it++ )
{
AlignmentParameters* ap = (*it)->alignmentParameters();
int npar=ap->numSelected();
HIPUserVariables* userpar = new HIPUserVariables(npar);
ap->setUserVariables(userpar);
}
// try to read in alignment parameters from a previous iteration
AlignablePositions theAlignablePositionsFromFile =
theIO.readAlignableAbsolutePositions(theAlignables,
(char*)salignedfile.c_str(),-1,ioerr);
int numAlignablesFromFile = theAlignablePositionsFromFile.size();
if (numAlignablesFromFile==0) { // file not there: first iteration
// set iteration number to 1
if (isCollector) theIteration=0;
else theIteration=1;
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] File not found => iteration "<<theIteration;
// get true (de-misaligned positions) and write to root file
// hardcoded iteration=1
theIO.writeAlignableOriginalPositions(theAlignables,
(char*)struefile.c_str(),1,false,ioerr);
// get misaligned positions and write to root file
// hardcoded iteration=1
theIO.writeAlignableAbsolutePositions(theAlignables,
(char*)smisalignedfile.c_str(),1,false,ioerr);
}
else { // there have been previous iterations
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Alignables Read "
<< numAlignablesFromFile;
// get iteration number from file
theIteration = readIterationFile(siterationfile);
// increase iteration
theIteration++;
edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Iteration increased by one!";
// now apply psotions of file from prev iteration
edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Apply positions from file ...";
theAlignmentParameterStore->applyAlignableAbsolutePositions(theAlignables,
theAlignablePositionsFromFile,ioerr);
}
edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Current Iteration number: "
<< theIteration;
// book root trees
bookRoot();
/*---------------------moved to terminate------------------------------
if (theLevels.size() > 0)
{
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
unsigned int nAlignable = theAlignables.size();
for (unsigned int i = 0; i < nAlignable; ++i)
{
const Alignable* ali = theAlignables[i];
AlignmentParameters* ap = ali->alignmentParameters();
HIPUserVariables* uservar =
dynamic_cast<HIPUserVariables*>(ap->userVariables());
for (unsigned int l = 0; l < theLevels.size(); ++l)
{
SurveyResidual res(*ali, theLevels[l], true);
if ( res.valid() )
{
AlgebraicSymMatrix invCov = res.inverseCovariance();
// variable for tree
AlgebraicVector sensResid = res.sensorResidual();
m3_Id = ali->id();
m3_ObjId = theLevels[l];
m3_par[0] = sensResid[0]; m3_par[1] = sensResid[1]; m3_par[2] = sensResid[2];
m3_par[3] = sensResid[3]; m3_par[4] = sensResid[4]; m3_par[5] = sensResid[5];
uservar->jtvj += invCov;
uservar->jtve += invCov * sensResid;
theTree3->Fill();
}
}
// align::LocalVectors residuals = res1.pointsResidual();
// unsigned int nPoints = residuals.size();
// for (unsigned int k = 0; k < nPoints; ++k)
// {
// AlgebraicMatrix J = term->survey()->derivatives(k);
// AlgebraicVector e(3); // local residual
// const align::LocalVector& lr = residuals[k];
// e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z();
// uservar->jtvj += invCov1.similarity(J);
// uservar->jtve += J * (invCov1 * e);
// }
}
}
//------------------------------------------------------*/
// set alignment position error
setAlignmentPositionError();
// run collector job if we are in parallel mode
if (isCollector) collector();
}
| void HIPAlignmentAlgorithm::terminate | ( | void | ) | [virtual] |
Call at end of job.
Implements AlignmentAlgorithmBase.
Definition at line 303 of file HIPAlignmentAlgorithm.cc.
References Alignable::alignmentParameters(), AlignmentParameterStore::applyParameters(), calcParameters(), fillRoot(), i, Alignable::id(), SurveyResidual::inverseCovariance(), ioerr, HIPUserVariables::jtve, HIPUserVariables::jtvj, prof2calltree::l, m3_Id, m3_ObjId, m3_par, salignedfile, SurveyResidual::sensorResidual(), AlignmentParameters::setValid(), siterationfile, sparameterfile, suvarfile, MultipleCompare::test, theAlignables, theAlignmentParameterStore, theFile, theFile2, theFile3, theIO, theIteration, theLevels, theTree, theTree2, theTree3, AlignmentParameters::userVariables(), SurveyResidual::valid(), AlignmentIORoot::writeAlignableAbsolutePositions(), AlignmentIORoot::writeAlignmentParameters(), HIPUserVariablesIORoot::writeHIPUserVariables(), and writeIterationFile().
{
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Terminating";
// calculating survey residuals
if (theLevels.size() > 0)
{
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
unsigned int nAlignable = theAlignables.size();
for (unsigned int i = 0; i < nAlignable; ++i)
{
const Alignable* ali = theAlignables[i];
AlignmentParameters* ap = ali->alignmentParameters();
HIPUserVariables* uservar =
dynamic_cast<HIPUserVariables*>(ap->userVariables());
for (unsigned int l = 0; l < theLevels.size(); ++l)
{
SurveyResidual res(*ali, theLevels[l], true);
if ( res.valid() )
{
AlgebraicSymMatrix invCov = res.inverseCovariance();
// variable for tree
AlgebraicVector sensResid = res.sensorResidual();
m3_Id = ali->id();
m3_ObjId = theLevels[l];
m3_par[0] = sensResid[0]; m3_par[1] = sensResid[1]; m3_par[2] = sensResid[2];
m3_par[3] = sensResid[3]; m3_par[4] = sensResid[4]; m3_par[5] = sensResid[5];
uservar->jtvj += invCov;
uservar->jtve += invCov * sensResid;
theTree3->Fill();
}
}
// align::LocalVectors residuals = res1.pointsResidual();
// unsigned int nPoints = residuals.size();
// for (unsigned int k = 0; k < nPoints; ++k)
// {
// AlgebraicMatrix J = term->survey()->derivatives(k);
// AlgebraicVector e(3); // local residual
// const align::LocalVector& lr = residuals[k];
// e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z();
// uservar->jtvj += invCov1.similarity(J);
// uservar->jtve += J * (invCov1 * e);
// }
}
}
// write user variables
HIPUserVariablesIORoot HIPIO;
HIPIO.writeHIPUserVariables (theAlignables,(char*)suvarfile.c_str(),
theIteration,false,ioerr);
// now calculate alignment corrections ...
int ialigned=0;
// iterate over alignment parameters
for(vector<Alignable*>::const_iterator
it=theAlignables.begin(); it!=theAlignables.end(); it++) {
Alignable* ali=(*it);
// Alignment parameters
AlignmentParameters* par = ali->alignmentParameters();
// try to calculate parameters
bool test = calcParameters(ali);
// if successful, apply parameters
if (test) {
edm::LogInfo("Alignment") << "now apply params";
theAlignmentParameterStore->applyParameters(ali);
// set these parameters 'valid'
ali->alignmentParameters()->setValid(true);
// increase counter
ialigned++;
}
else par->setValid(false);
}
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::terminate] Aligned units: " << ialigned;
// fill alignable wise root tree
fillRoot();
edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Writing aligned parameters to file: " << theAlignables.size();
// write new absolute positions to disk
theIO.writeAlignableAbsolutePositions(theAlignables,
(char*)salignedfile.c_str(),theIteration,false,ioerr);
// write alignment parameters to disk
theIO.writeAlignmentParameters(theAlignables,
(char*)sparameterfile.c_str(),theIteration,false,ioerr);
// write iteration number to file
writeIterationFile(siterationfile,theIteration);
// write out trees and close root file
// eventwise tree
theFile->cd();
theTree->Write();
delete theFile;
if (theLevels.size() > 0){
theFile3->cd();
theTree3->Write();
delete theFile3;
}
// alignable-wise tree is only filled once
if (theIteration==1) { // only for 1st iteration
theFile2->cd();
theTree2->Write();
delete theFile2;
}
}
| void HIPAlignmentAlgorithm::writeIterationFile | ( | std::string | filename, |
| int | iter | ||
| ) | [private] |
Definition at line 932 of file HIPAlignmentAlgorithm.cc.
References dbtoconf::out.
Referenced by terminate().
{
ofstream outIterFile((char*)(filename.c_str()), ios::out);
if (!outIterFile) {
edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::writeIterationFile] ERROR: Unable to write Iteration file";
}
else {
outIterFile << iter;
edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::writeIterationFile] writing iteration number to file: " << iter;
}
outIterFile.close();
}
int HIPAlignmentAlgorithm::ioerr [private] |
Definition at line 79 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), startNewLoop(), and terminate().
bool HIPAlignmentAlgorithm::isCollector [private] |
Definition at line 101 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), run(), and startNewLoop().
float HIPAlignmentAlgorithm::m2_Eta [private] |
Definition at line 125 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
align::ID HIPAlignmentAlgorithm::m2_Id [private] |
Definition at line 126 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
int HIPAlignmentAlgorithm::m2_Layer [private] |
Definition at line 124 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
int HIPAlignmentAlgorithm::m2_Nhit [private] |
Definition at line 124 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
Definition at line 127 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
float HIPAlignmentAlgorithm::m2_Phi [private] |
Definition at line 125 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
int HIPAlignmentAlgorithm::m2_Type [private] |
Definition at line 124 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
float HIPAlignmentAlgorithm::m2_Xpos [private] |
Definition at line 125 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
float HIPAlignmentAlgorithm::m2_Ypos [private] |
Definition at line 125 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
float HIPAlignmentAlgorithm::m2_Zpos [private] |
Definition at line 125 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and fillRoot().
align::ID HIPAlignmentAlgorithm::m3_Id [private] |
Definition at line 130 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and terminate().
Definition at line 131 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and terminate().
float HIPAlignmentAlgorithm::m3_par[6] [private] |
Definition at line 132 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and terminate().
float HIPAlignmentAlgorithm::m_Chi2n[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
float HIPAlignmentAlgorithm::m_d0[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
float HIPAlignmentAlgorithm::m_dz[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
float HIPAlignmentAlgorithm::m_Eta[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
int HIPAlignmentAlgorithm::m_Nhits[MAXREC] [private] |
Definition at line 120 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
int HIPAlignmentAlgorithm::m_nhPXB[MAXREC] [private] |
Definition at line 120 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
int HIPAlignmentAlgorithm::m_nhPXF[MAXREC] [private] |
Definition at line 120 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
int HIPAlignmentAlgorithm::m_Ntracks [private] |
Definition at line 120 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
float HIPAlignmentAlgorithm::m_P[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
float HIPAlignmentAlgorithm::m_Phi[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
float HIPAlignmentAlgorithm::m_Pt[MAXREC] [private] |
Definition at line 121 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), and run().
const int HIPAlignmentAlgorithm::MAXREC = 99 [static, private] |
Definition at line 118 of file HIPAlignmentAlgorithm.h.
Referenced by fillEventwiseTree(), and run().
std::string HIPAlignmentAlgorithm::outfile [private] |
Definition at line 87 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::outfile2 [private] |
Definition at line 87 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::outpath [private] |
Definition at line 87 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::salignedfile [private] |
Definition at line 88 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), startNewLoop(), and terminate().
std::string HIPAlignmentAlgorithm::siterationfile [private] |
Definition at line 88 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), startNewLoop(), and terminate().
std::string HIPAlignmentAlgorithm::smisalignedfile [private] |
Definition at line 88 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and startNewLoop().
std::string HIPAlignmentAlgorithm::sparameterfile [private] |
Definition at line 87 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and terminate().
std::string HIPAlignmentAlgorithm::ssurveyfile [private] |
Definition at line 88 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::struefile [private] |
Definition at line 88 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and startNewLoop().
std::string HIPAlignmentAlgorithm::suvarfile [private] |
Definition at line 87 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and terminate().
Definition at line 76 of file HIPAlignmentAlgorithm.h.
Referenced by initialize(), and run().
std::vector<Alignable*> HIPAlignmentAlgorithm::theAlignables [private] |
Definition at line 75 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), fillRoot(), initialize(), startNewLoop(), and terminate().
Definition at line 74 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 93 of file HIPAlignmentAlgorithm.h.
Referenced by initialize(), and setAlignmentPositionError().
std::vector<edm::ParameterSet> HIPAlignmentAlgorithm::theAPEParameterSet [private] |
Definition at line 92 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and initialize().
bool HIPAlignmentAlgorithm::theApplyAPE [private] |
Definition at line 91 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), initialize(), and setAlignmentPositionError().
int HIPAlignmentAlgorithm::theCollectorNJobs [private] |
Definition at line 102 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), and HIPAlignmentAlgorithm().
std::string HIPAlignmentAlgorithm::theCollectorPath [private] |
Definition at line 103 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), and HIPAlignmentAlgorithm().
int HIPAlignmentAlgorithm::theCurrentPrescale [private] |
Definition at line 104 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and run().
int HIPAlignmentAlgorithm::theEventPrescale [private] |
Definition at line 104 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), and run().
TFile* HIPAlignmentAlgorithm::theFile [private] |
Definition at line 110 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), run(), and terminate().
TFile* HIPAlignmentAlgorithm::theFile2 [private] |
Definition at line 112 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillRoot(), and terminate().
TFile* HIPAlignmentAlgorithm::theFile3 [private] |
Definition at line 114 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and terminate().
bool HIPAlignmentAlgorithm::theFillTrackMonitoring [private] |
Definition at line 105 of file HIPAlignmentAlgorithm.h.
Referenced by collector(), HIPAlignmentAlgorithm(), and run().
AlignmentIORoot HIPAlignmentAlgorithm::theIO [private] |
Definition at line 78 of file HIPAlignmentAlgorithm.h.
Referenced by startNewLoop(), and terminate().
int HIPAlignmentAlgorithm::theIteration [private] |
Definition at line 80 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), collector(), setAlignmentPositionError(), startNewLoop(), and terminate().
std::vector<align::StructureType> HIPAlignmentAlgorithm::theLevels [private] |
Definition at line 107 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), HIPAlignmentAlgorithm(), and terminate().
double HIPAlignmentAlgorithm::theMaxAllowedHitPull [private] |
Definition at line 95 of file HIPAlignmentAlgorithm.h.
Referenced by HIPAlignmentAlgorithm(), processHit1D(), and processHit2D().
double HIPAlignmentAlgorithm::theMaxRelParameterError [private] |
Definition at line 99 of file HIPAlignmentAlgorithm.h.
Referenced by calcParameters(), and HIPAlignmentAlgorithm().
int HIPAlignmentAlgorithm::theMinimumNumberOfHits [private] |
Definition at line 97 of file HIPAlignmentAlgorithm.h.
Referenced by calcParameters(), and HIPAlignmentAlgorithm().
TTree* HIPAlignmentAlgorithm::theTree [private] |
Definition at line 111 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillEventwiseTree(), run(), and terminate().
TTree* HIPAlignmentAlgorithm::theTree2 [private] |
Definition at line 113 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), fillRoot(), and terminate().
TTree* HIPAlignmentAlgorithm::theTree3 [private] |
Definition at line 115 of file HIPAlignmentAlgorithm.h.
Referenced by bookRoot(), and terminate().
bool HIPAlignmentAlgorithm::verbose [private] |
Definition at line 85 of file HIPAlignmentAlgorithm.h.
1.7.3