CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes
HIPAlignmentAlgorithm Class Reference

#include <HIPAlignmentAlgorithm.h>

Inheritance diagram for HIPAlignmentAlgorithm:
AlignmentAlgorithmBase

Public Member Functions

 HIPAlignmentAlgorithm (const edm::ParameterSet &cfg)
 Constructor. More...
 
void initialize (const edm::EventSetup &setup, AlignableTracker *tracker, AlignableMuon *muon, AlignableExtras *extras, AlignmentParameterStore *store)
 Call at beginning of job. More...
 
void run (const edm::EventSetup &setup, const EventInfo &eventInfo)
 Run the algorithm. More...
 
void startNewLoop (void)
 Called at start of new loop. More...
 
void terminate (const edm::EventSetup &setup)
 Call at end of job. More...
 
 ~HIPAlignmentAlgorithm ()
 Destructor. More...
 
- Public Member Functions inherited from AlignmentAlgorithmBase
virtual bool addCalibrations (const std::vector< IntegratedCalibrationBase * > &iCals)
 
 AlignmentAlgorithmBase (const edm::ParameterSet &cfg)
 Constructor. More...
 
virtual void beginLuminosityBlock (const edm::EventSetup &setup)
 called at begin of luminosity block (no lumi block info passed yet) More...
 
virtual void beginRun (const edm::EventSetup &setup)
 called at begin of run More...
 
virtual void endLuminosityBlock (const edm::EventSetup &setup)
 called at end of luminosity block (no lumi block info passed yet) More...
 
virtual void endRun (const EndRunInfo &runInfo, const edm::EventSetup &setup)
 called at end of run - order of arguments like in EDProducer etc. More...
 
virtual bool setParametersForRunRange (const RunRange &rr)
 
virtual ~AlignmentAlgorithmBase ()
 Destructor. More...
 

Private Member Functions

void bookRoot (void)
 
double calcAPE (double *par, int iter, double function)
 
bool calcParameters (Alignable *ali)
 
void collector (void)
 
int fillEventwiseTree (const char *filename, int iter, int ierr)
 
void fillRoot (const edm::EventSetup &setup)
 
bool processHit1D (const AlignableDetOrUnitPtr &alidet, const Alignable *ali, const TrajectoryStateOnSurface &tsos, const TrackingRecHit *hit)
 
bool processHit2D (const AlignableDetOrUnitPtr &alidet, const Alignable *ali, const TrajectoryStateOnSurface &tsos, const TrackingRecHit *hit)
 
int readIterationFile (std::string filename)
 
void setAlignmentPositionError (void)
 
void writeIterationFile (std::string filename, int iter)
 

Private Attributes

int ioerr
 
bool isCollector
 
float m2_Eta
 
align::ID m2_Id
 
int m2_Layer
 
int m2_Nhit
 
align::StructureType m2_ObjId
 
float m2_Phi
 
int m2_Type
 
float m2_Xpos
 
float m2_Ypos
 
float m2_Zpos
 
align::ID m3_Id
 
align::StructureType m3_ObjId
 
float m3_par [6]
 
float m_Chi2n [MAXREC]
 
float m_d0 [MAXREC]
 
float m_dz [MAXREC]
 
float m_Eta [MAXREC]
 
int m_Nhits [MAXREC]
 
int m_nhPXB [MAXREC]
 
int m_nhPXF [MAXREC]
 
int m_Ntracks
 
float m_P [MAXREC]
 
float m_Phi [MAXREC]
 
float m_Pt [MAXREC]
 
std::string outfile
 
std::string outfile2
 
std::string outpath
 
std::string salignedfile
 
std::string siterationfile
 
std::string smisalignedfile
 
std::string sparameterfile
 
std::string ssurveyfile
 
std::string struefile
 
std::string suvarfile
 
AlignableNavigatortheAlignableDetAccessor
 
std::vector< Alignable * > theAlignables
 
AlignmentParameterStoretheAlignmentParameterStore
 
std::vector< std::pair
< std::vector< Alignable * >
, std::vector< double > > > 
theAPEParameters
 
std::vector< edm::ParameterSettheAPEParameterSet
 
bool theApplyAPE
 
int theCollectorNJobs
 
std::string theCollectorPath
 
int theCurrentPrescale
 
int theEventPrescale
 
TFile * theFile
 
TFile * theFile2
 
TFile * theFile3
 
bool theFillTrackMonitoring
 
AlignmentIORoot theIO
 
int theIteration
 
std::vector< align::StructureTypetheLevels
 
double theMaxAllowedHitPull
 
double theMaxRelParameterError
 
int theMinimumNumberOfHits
 
TTree * theTree
 
TTree * theTree2
 
TTree * theTree3
 
bool verbose
 

Static Private Attributes

static const int MAXREC = 99
 

Additional Inherited Members

- Public Types inherited from AlignmentAlgorithmBase
typedef std::pair< const
Trajectory *, const
reco::Track * > 
ConstTrajTrackPair
 
typedef std::vector
< ConstTrajTrackPair
ConstTrajTrackPairCollection
 
typedef cond::RealTimeType
< cond::runnumber >::type 
RunNumber
 
typedef std::pair< RunNumber,
RunNumber
RunRange
 

Detailed Description

Definition at line 21 of file HIPAlignmentAlgorithm.h.

Constructor & Destructor Documentation

HIPAlignmentAlgorithm::HIPAlignmentAlgorithm ( const edm::ParameterSet cfg)

Constructor.

Definition at line 37 of file HIPAlignmentAlgorithm.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), isCollector, outfile, outfile2, outpath, alignCSCRings::s, salignedfile, siterationfile, smisalignedfile, sparameterfile, ssurveyfile, AlCaHLTBitMon_QueryRunRegistry::string, AlignableObjectId::stringToId(), struefile, suvarfile, theAPEParameterSet, theApplyAPE, theCollectorNJobs, theCollectorPath, theCurrentPrescale, theEventPrescale, theFillTrackMonitoring, theLevels, theMaxAllowedHitPull, theMaxRelParameterError, and theMinimumNumberOfHits.

37  :
39 {
40 
41  // parse parameters
42 
43  verbose = cfg.getParameter<bool>("verbosity");
44 
45  outpath = cfg.getParameter<std::string>("outpath");
46  outfile = cfg.getParameter<std::string>("outfile");
47  outfile2 = cfg.getParameter<std::string>("outfile2");
48  struefile = cfg.getParameter<std::string>("trueFile");
49  smisalignedfile = cfg.getParameter<std::string>("misalignedFile");
50  salignedfile = cfg.getParameter<std::string>("alignedFile");
51  siterationfile = cfg.getParameter<std::string>("iterationFile");
52  suvarfile = cfg.getParameter<std::string>("uvarFile");
53  sparameterfile = cfg.getParameter<std::string>("parameterFile");
54  ssurveyfile = cfg.getParameter<std::string>("surveyFile");
55 
56  outfile =outpath+outfile;//Eventwise tree
57  outfile2 =outpath+outfile2;//Alignablewise tree
58  struefile =outpath+struefile;
59  smisalignedfile=outpath+smisalignedfile;
60  salignedfile =outpath+salignedfile;
61  siterationfile =outpath+siterationfile;
62  suvarfile =outpath+suvarfile;
63  sparameterfile =outpath+sparameterfile;
64  ssurveyfile =outpath+ssurveyfile;
65 
66  // parameters for APE
67  theApplyAPE = cfg.getParameter<bool>("applyAPE");
68  theAPEParameterSet = cfg.getParameter<std::vector<edm::ParameterSet> >("apeParam");
69 
70  theMaxAllowedHitPull = cfg.getParameter<double>("maxAllowedHitPull");
71  theMinimumNumberOfHits = cfg.getParameter<int>("minimumNumberOfHits");
72  theMaxRelParameterError = cfg.getParameter<double>("maxRelParameterError");
73 
74  // for collector mode (parallel processing)
75  isCollector=cfg.getParameter<bool>("collectorActive");
76  theCollectorNJobs=cfg.getParameter<int>("collectorNJobs");
77  theCollectorPath=cfg.getParameter<std::string>("collectorPath");
78  theFillTrackMonitoring=cfg.getUntrackedParameter<bool>("fillTrackMonitoring");
79 
80  if (isCollector) edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Collector mode";
81 
82  theEventPrescale = cfg.getParameter<int>("eventPrescale");
84 
85  for (std::string &s : cfg.getUntrackedParameter<std::vector<std::string> >("surveyResiduals")) {
87  }
88 
89  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] constructed.";
90 
91 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
AlignmentAlgorithmBase(const edm::ParameterSet &cfg)
Constructor.
static align::StructureType stringToId(const char *)
std::vector< edm::ParameterSet > theAPEParameterSet
std::vector< align::StructureType > theLevels
HIPAlignmentAlgorithm::~HIPAlignmentAlgorithm ( )
inline

Destructor.

Definition at line 30 of file HIPAlignmentAlgorithm.h.

30 {};

Member Function Documentation

void HIPAlignmentAlgorithm::bookRoot ( void  )
private

Definition at line 1029 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, theTree3, and cond::persistency::GLOBAL_TAG::tname.

Referenced by startNewLoop().

1030 {
1031  // create ROOT files
1032  theFile = new TFile(outfile.c_str(),"update");
1033  theFile->cd();
1034 
1035  // book event-wise ROOT Tree
1036 
1037  TString tname="T1";
1038  char iterString[5];
1039  snprintf(iterString, sizeof(iterString), "%i",theIteration);
1040  tname.Append("_");
1041  tname.Append(iterString);
1042 
1043  theTree = new TTree(tname,"Eventwise tree");
1044 
1045  //theTree->Branch("Run", &m_Run, "Run/I");
1046  //theTree->Branch("Event", &m_Event, "Event/I");
1047  theTree->Branch("Ntracks", &m_Ntracks, "Ntracks/I");
1048  theTree->Branch("Nhits", m_Nhits, "Nhits[Ntracks]/I");
1049  theTree->Branch("nhPXB", m_nhPXB, "nhPXB[Ntracks]/I");
1050  theTree->Branch("nhPXF", m_nhPXF, "nhPXF[Ntracks]/I");
1051  theTree->Branch("Pt", m_Pt, "Pt[Ntracks]/F");
1052  theTree->Branch("P", m_P, "P[Ntracks]/F");
1053  theTree->Branch("Eta", m_Eta, "Eta[Ntracks]/F");
1054  theTree->Branch("Phi", m_Phi, "Phi[Ntracks]/F");
1055  theTree->Branch("Chi2n", m_Chi2n, "Chi2n[Ntracks]/F");
1056  theTree->Branch("d0", m_d0, "d0[Ntracks]/F");
1057  theTree->Branch("dz", m_dz, "dz[Ntracks]/F");
1058 
1059  // book Alignable-wise ROOT Tree
1060 
1061  theFile2 = new TFile(outfile2.c_str(),"update");
1062  theFile2->cd();
1063 
1064  theTree2 = new TTree("T2","Alignablewise tree");
1065 
1066  theTree2->Branch("Nhit", &m2_Nhit, "Nhit/I");
1067  theTree2->Branch("Type", &m2_Type, "Type/I");
1068  theTree2->Branch("Layer", &m2_Layer, "Layer/I");
1069  theTree2->Branch("Xpos", &m2_Xpos, "Xpos/F");
1070  theTree2->Branch("Ypos", &m2_Ypos, "Ypos/F");
1071  theTree2->Branch("Zpos", &m2_Zpos, "Zpos/F");
1072  theTree2->Branch("Eta", &m2_Eta, "Eta/F");
1073  theTree2->Branch("Phi", &m2_Phi, "Phi/F");
1074  theTree2->Branch("Id", &m2_Id, "Id/i");
1075  theTree2->Branch("ObjId", &m2_ObjId, "ObjId/I");
1076 
1077  // book survey-wise ROOT Tree only if survey is enabled
1078  if (theLevels.size() > 0){
1079 
1080  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Survey trees booked.";
1081  theFile3 = new TFile(ssurveyfile.c_str(),"update");
1082  theFile3->cd();
1083  theTree3 = new TTree(tname, "Survey Tree");
1084  theTree3->Branch("Id", &m3_Id, "Id/i");
1085  theTree3->Branch("ObjId", &m3_ObjId, "ObjId/I");
1086  theTree3->Branch("Par", &m3_par, "Par[6]/F");
1087  }
1088 
1089  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Root trees booked.";
1090 
1091 }
static char const * tname
Definition: GTSchema.h:13
align::StructureType m2_ObjId
align::StructureType m3_ObjId
std::vector< align::StructureType > theLevels
double HIPAlignmentAlgorithm::calcAPE ( double *  par,
int  iter,
double  function 
)
private

Definition at line 1001 of file HIPAlignmentAlgorithm.cc.

References create_public_lumi_plots::exp, bookConverter::max, funct::pow(), and relval_parameters_module::step.

Referenced by setAlignmentPositionError().

1002 {
1003  double diter=(double)iter;
1004 
1005  // The following floating-point equality check is safe because this
1006  // "0." and this "1." are generated by the compiler, in the very
1007  // same file. Whatever approximization scheme it uses to turn "1."
1008  // into 0.9999999999998 in the HIPAlignmentAlgorithm::initialize is
1009  // also used here. If I'm wrong, you'll get an assertion.
1010  if (function == 0.) {
1011  return std::max(par[1],par[0]+((par[1]-par[0])/par[2])*diter);
1012  }
1013  else if (function == 1.) {
1014  return std::max(0.,par[0]*(exp(-pow(diter,par[1])/par[2])));
1015  }
1016  else if (function == 2.) {
1017  int ipar2 = (int) par[2];
1018  int step = iter/ipar2;
1019  double dstep = (double) step;
1020  return std::max(0.0, par[0] - par[1]*dstep);
1021  }
1022  else assert(false); // should have been caught in the constructor
1023 }
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
bool HIPAlignmentAlgorithm::calcParameters ( Alignable ali)
private

Definition at line 1177 of file HIPAlignmentAlgorithm.cc.

References funct::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().

1178 {
1179  // Alignment parameters
1181  // access user variables
1182  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(par->userVariables());
1183  int nhit = uservar->nhit;
1184  // The following variable is needed for the extended 1D/2D hit fix using
1185  // matrix shrinkage and expansion
1186  // int hitdim = uservar->hitdim;
1187 
1188  if (nhit < theMinimumNumberOfHits) {
1189  par->setValid(false);
1190  return false;
1191  }
1192 
1193  AlgebraicSymMatrix jtvj = uservar->jtvj;
1194  AlgebraicVector jtve = uservar->jtve;
1195 
1196  // Shrink input in case of 1D hits and 'v' selected
1197  // in alignment parameters
1198  // if (hitdim==1 && selector[1]==true) {
1199  // int iremove = 1;
1200  // if (selector[0]==false) iremove--;
1201 
1202  // AlgebraicSymMatrix tempjtvj(jtvj.num_row()-1);
1203  // int nr = 0, nc = 0;
1204  // for (int r=0;r<jtvj.num_row();r++) {
1205  // if (r==iremove) continue;
1206  // nc = 0;
1207  // for (int c=0;c<jtvj.num_col();c++) {
1208  // if (c==iremove) continue;
1209  // tempjtvj[nr][nc] = jtvj[r][c];
1210  // nc++;
1211  // }
1212  // nr++;
1213  // }
1214  // jtvj = tempjtvj;
1215 
1216  // AlgebraicVector tempjtve(jtve.num_row()-1);
1217  // nr = 0;
1218  // for (int r=0;r<jtve.num_row();r++) {
1219  // if (r==iremove) continue;
1220  // tempjtve[nr] = jtve[r];
1221  // nr++;
1222  // }
1223  // jtve = tempjtve;
1224  // }
1225 
1226  int ierr;
1227  AlgebraicSymMatrix jtvjinv = jtvj.inverse(ierr);
1228 
1229  if (ierr !=0) {
1230  edm::LogError("Alignment") << "Matrix inversion failed!";
1231  return false;
1232  }
1233 
1234  // these are the alignment corrections+covariance (for selected params)
1235  AlgebraicVector params = - (jtvjinv * jtve);
1236  AlgebraicSymMatrix cov = jtvjinv;
1237 
1238  edm::LogInfo("Alignment") << "parameters " << params;
1239 
1240  // errors of parameters
1241  int npar = params.num_row();
1242  AlgebraicVector paramerr(npar);
1243  AlgebraicVector relerr(npar);
1244  for (int i=0;i<npar;i++) {
1245  if (abs(cov[i][i])>0) paramerr[i] = sqrt(abs(cov[i][i]));
1246  else paramerr[i] = params[i];
1247  relerr[i] = abs(paramerr[i]/params[i]);
1248  if (relerr[i] >= theMaxRelParameterError) {
1249  params[i] = 0;
1250  paramerr[i]=0;
1251  }
1252  }
1253 
1254  // expand output in case of 1D hits and 'v' selected
1255  // in alignment parameters
1256  // if (hitdim==1 && selector[1]==true) {
1257  // int iremove = 1;
1258  // if (selector[0]==false) iremove--;
1259 
1260  // AlgebraicSymMatrix tempcov(cov.num_row()+1,0);
1261  // int nr = 0, nc = 0;
1262  // for (int r=0;r<cov.num_row();r++) {
1263  // if (r==iremove) nr++;
1264  // nc = 0;
1265  // for (int c=0;c<cov.num_col();c++) {
1266  // if (c==iremove) nc++;
1267  // tempcov[nr][nc] = cov[r][c];
1268  // nc++;
1269  // }
1270  // nr++;
1271  // }
1272  // cov = tempcov;
1273 
1274  // AlgebraicVector tempparams(params.num_row()+1,0);
1275  // nr = 0;
1276  // for (int r=0;r<params.num_row();r++) {
1277  // if (r==iremove) nr++;
1278  // tempparams[nr] = params[r];
1279  // nr++;
1280  // }
1281  // params = tempparams;
1282  // }
1283 
1284  // store alignment parameters
1285  AlignmentParameters* parnew = par->cloneFromSelected(params,cov);
1286  ali->setAlignmentParameters(parnew);
1287  parnew->setValid(true);
1288 
1289  return true;
1290 }
int i
Definition: DBlmapReader.cc:9
AlgebraicVector jtve
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
AlgebraicSymMatrix jtvj
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
void setAlignmentParameters(AlignmentParameters *dap)
Set the AlignmentParameters.
Definition: Alignable.cc:81
void setValid(bool v)
Set validity flag.
T sqrt(T t)
Definition: SSEVec.h:48
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
CLHEP::HepVector AlgebraicVector
CLHEP::HepSymMatrix AlgebraicSymMatrix
virtual AlignmentParameters * cloneFromSelected(const AlgebraicVector &par, const AlgebraicSymMatrix &cov) const =0
void HIPAlignmentAlgorithm::collector ( void  )
private

Definition at line 1294 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(), contentValuesCheck::ss, AlCaHLTBitMon_QueryRunRegistry::string, theAlignables, theAlignmentParameterStore, theCollectorNJobs, theCollectorPath, theFillTrackMonitoring, theIteration, and AlignmentParameters::userVariables().

Referenced by startNewLoop().

1295 {
1296  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::collector] called for iteration "
1297  << theIteration << std::endl;
1298 
1299  HIPUserVariablesIORoot HIPIO;
1300 
1301  for (int ijob=1;ijob<=theCollectorNJobs;ijob++) {
1302 
1303  edm::LogWarning("Alignment") << "reading uservar for job " << ijob;
1304 
1305  std::stringstream ss;
1306  std::string str;
1307  ss << ijob;
1308  ss >> str;
1309  std::string uvfile = theCollectorPath+"/job"+str+"/IOUserVariables.root";
1310 
1311  std::vector<AlignmentUserVariables*> uvarvec =
1312  HIPIO.readHIPUserVariables(theAlignables, uvfile.c_str(),
1313  theIteration, ioerr);
1314 
1315  if (ioerr!=0) {
1316  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::collector] could not read user variable files for job "
1317  << ijob;
1318  continue;
1319  }
1320 
1321  // add
1322  std::vector<AlignmentUserVariables*> uvarvecadd;
1323  std::vector<AlignmentUserVariables*>::const_iterator iuvarnew=uvarvec.begin();
1324  for (std::vector<Alignable*>::const_iterator it=theAlignables.begin();
1325  it!=theAlignables.end();
1326  ++it) {
1327  Alignable* ali = *it;
1329 
1330  HIPUserVariables* uvarold = dynamic_cast<HIPUserVariables*>(ap->userVariables());
1331  HIPUserVariables* uvarnew = dynamic_cast<HIPUserVariables*>(*iuvarnew);
1332 
1333  HIPUserVariables* uvar = uvarold->clone();
1334  if (uvarnew!=0) {
1335  uvar->nhit = (uvarold->nhit)+(uvarnew->nhit);
1336  uvar->jtvj = (uvarold->jtvj)+(uvarnew->jtvj);
1337  uvar->jtve = (uvarold->jtve)+(uvarnew->jtve);
1338  uvar->alichi2 = (uvarold->alichi2)+(uvarnew->alichi2);
1339  uvar->alindof = (uvarold->alindof)+(uvarnew->alindof);
1340  delete uvarnew;
1341  }
1342 
1343  uvarvecadd.push_back(uvar);
1344  iuvarnew++;
1345  }
1346 
1348 
1349  //fill Eventwise Tree
1350  if (theFillTrackMonitoring) {
1351  uvfile = theCollectorPath+"/job"+str+"/HIPAlignmentEvents.root";
1352  edm::LogWarning("Alignment") << "Added to the tree "
1353  << fillEventwiseTree(uvfile.c_str(), theIteration, ioerr)
1354  << "tracks";
1355  }
1356 
1357  }//end loop on jobs
1358 }
void attachUserVariables(const align::Alignables &alivec, const std::vector< AlignmentUserVariables * > &uvarvec, int &ierr)
Attach User Variables to given alignables.
AlignmentParameterStore * theAlignmentParameterStore
int fillEventwiseTree(const char *filename, int iter, int ierr)
AlgebraicVector jtve
HIPUserVariables * clone(void) const
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
AlgebraicSymMatrix jtvj
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
std::vector< AlignmentUserVariables * > readHIPUserVariables(const Alignables &alivec, const char *filename, int iter, int &ierr)
std::vector< Alignable * > theAlignables
int HIPAlignmentAlgorithm::fillEventwiseTree ( const char *  filename,
int  iter,
int  ierr 
)
private

Definition at line 1362 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().

1363 {
1364  int totntrk = 0;
1365  char treeName[64];
1366  snprintf(treeName, sizeof(treeName), "T1_%d", iter);
1367  //open the file "HIPAlignmentEvents.root" in the job directory
1368  TFile *jobfile = new TFile(filename, "READ");
1369  //grab the tree corresponding to this iteration
1370  TTree *jobtree = (TTree*)jobfile->Get(treeName);
1371  //address and read the variables
1372  static const int nmaxtrackperevent = 1000;
1373  int jobNtracks, jobNhitspertrack[nmaxtrackperevent], jobnhPXB[nmaxtrackperevent], jobnhPXF[nmaxtrackperevent];
1374  float jobP[nmaxtrackperevent], jobPt[nmaxtrackperevent], jobEta[nmaxtrackperevent] , jobPhi[nmaxtrackperevent];
1375  float jobd0[nmaxtrackperevent], jobdz[nmaxtrackperevent] , jobChi2n[nmaxtrackperevent];
1376 
1377  jobtree->SetBranchAddress("Ntracks", &jobNtracks);
1378  jobtree->SetBranchAddress("Nhits", jobNhitspertrack);
1379  jobtree->SetBranchAddress("nhPXB", jobnhPXB);
1380  jobtree->SetBranchAddress("nhPXF", jobnhPXF);
1381  jobtree->SetBranchAddress("Pt", jobPt);
1382  jobtree->SetBranchAddress("P", jobP);
1383  jobtree->SetBranchAddress("d0", jobd0);
1384  jobtree->SetBranchAddress("dz", jobdz);
1385  jobtree->SetBranchAddress("Eta", jobEta);
1386  jobtree->SetBranchAddress("Phi", jobPhi);
1387  jobtree->SetBranchAddress("Chi2n", jobChi2n);
1388  int ievent = 0;
1389  for (ievent=0;ievent<jobtree->GetEntries();++ievent) {
1390  jobtree->GetEntry(ievent);
1391 
1392  //fill the collector tree with them
1393 
1394  // TO BE IMPLEMENTED: a prescale factor like in run()
1395  m_Ntracks = jobNtracks;
1396  int ntrk = 0;
1397  while (ntrk<m_Ntracks) {
1398  if (ntrk<MAXREC) {
1399  totntrk = ntrk+1;
1400  m_Nhits[ntrk] = jobNhitspertrack[ntrk];
1401  m_Pt[ntrk] = jobPt[ntrk];
1402  m_P[ntrk] = jobP[ntrk];
1403  m_nhPXB[ntrk] = jobnhPXB[ntrk];
1404  m_nhPXF[ntrk] = jobnhPXF[ntrk];
1405  m_Eta[ntrk] = jobEta[ntrk];
1406  m_Phi[ntrk] = jobPhi[ntrk];
1407  m_Chi2n[ntrk] = jobChi2n[ntrk];
1408  m_d0[ntrk] = jobd0[ntrk];
1409  m_dz[ntrk] = jobdz[ntrk];
1410  }//end if j<MAXREC
1411  else{
1412  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::fillEventwiseTree] Number of tracks in Eventwise tree exceeds MAXREC: "
1413  << m_Ntracks << " Skipping exceeding tracks.";
1414  ntrk = m_Ntracks+1;
1415  }
1416  ++ntrk;
1417  }//end while loop
1418  theTree->Fill();
1419  }//end loop on i - entries in the job tree
1420 
1421  //clean up
1422  delete jobtree;
1423  delete jobfile;
1424 
1425  return totntrk;
1426 }//end fillEventwiseTree
tuple filename
Definition: lut2db_cfg.py:20
void HIPAlignmentAlgorithm::fillRoot ( const edm::EventSetup setup)
private

Definition at line 1096 of file HIPAlignmentAlgorithm.cc.

References Alignable::alignableObjectId(), Alignable::alignmentParameters(), PV3DBase< T, PVType, FrameType >::eta(), edm::EventSetup::get(), 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(), GloballyPositioned< T >::position(), edm::ESHandle< class >::product(), 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().

1097 {
1098  using std::setw;
1099  theFile2->cd();
1100 
1101  int naligned=0;
1102 
1103  //Retrieve tracker topology from geometry
1104  edm::ESHandle<TrackerTopology> tTopoHandle;
1105  iSetup.get<IdealGeometryRecord>().get(tTopoHandle);
1106  const TrackerTopology* const tTopo = tTopoHandle.product();
1107 
1108  for (std::vector<Alignable*>::const_iterator it=theAlignables.begin();
1109  it!=theAlignables.end();
1110  ++it) {
1111  Alignable* ali = (*it);
1113 
1114  // consider only those parameters classified as 'valid'
1115  if (dap->isValid()) {
1116 
1117  // get number of hits from user variable
1118  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(dap->userVariables());
1119  m2_Nhit = uservar->nhit;
1120 
1121  // get type/layer
1122  std::pair<int,int> tl = theAlignmentParameterStore->typeAndLayer(ali, tTopo);
1123  m2_Type = tl.first;
1124  m2_Layer = tl.second;
1125 
1126  // get identifier (as for IO)
1127  m2_Id = ali->id();
1128  m2_ObjId = ali->alignableObjectId();
1129 
1130  // get position
1131  GlobalPoint pos = ali->surface().position();
1132  m2_Xpos = pos.x();
1133  m2_Ypos = pos.y();
1134  m2_Zpos = pos.z();
1135  m2_Eta = pos.eta();
1136  m2_Phi = pos.phi();
1137 
1138  AlgebraicVector pars = dap->parameters();
1139 
1140  if (verbose) {
1141  edm::LogVerbatim("Alignment")
1142  << "------------------------------------------------------------------------\n"
1143  << " ALIGNABLE: " << setw(6) << naligned
1144  << '\n'
1145  << "hits: " << setw(4) << m2_Nhit
1146  << " type: " << setw(4) << m2_Type
1147  << " layer: " << setw(4) << m2_Layer
1148  << " id: " << setw(4) << m2_Id
1149  << " objId: " << setw(4) << m2_ObjId
1150  << '\n'
1151  << std::fixed << std::setprecision(5)
1152  << "x,y,z: "
1153  << setw(12) << m2_Xpos
1154  << setw(12) << m2_Ypos
1155  << setw(12) << m2_Zpos
1156  << " eta,phi: "
1157  << setw(12) << m2_Eta
1158  << setw(12) << m2_Phi
1159  << '\n'
1160  << "params: "
1161  << setw(12) << pars[0]
1162  << setw(12) << pars[1]
1163  << setw(12) << pars[2]
1164  << setw(12) << pars[3]
1165  << setw(12) << pars[4]
1166  << setw(12) << pars[5];
1167  }
1168 
1169  naligned++;
1170  theTree2->Fill();
1171  }
1172  }
1173 }
align::ID id() const
Return the ID of Alignable, i.e. DetId of &#39;first&#39; component GeomDet(Unit).
Definition: Alignable.h:180
AlignmentParameterStore * theAlignmentParameterStore
std::pair< int, int > typeAndLayer(const Alignable *ali, const TrackerTopology *tTopo) const
Obtain type and layer from Alignable.
align::StructureType m2_ObjId
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
const AlgebraicVector & parameters(void) const
Get alignment parameters.
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
T z() const
Definition: PV3DBase.h:64
virtual StructureType alignableObjectId() const =0
Return the alignable type identifier.
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:126
CLHEP::HepVector AlgebraicVector
T const * product() const
Definition: ESHandle.h:62
T eta() const
Definition: PV3DBase.h:76
bool isValid(void) const
Get validity flag.
T x() const
Definition: PV3DBase.h:62
const PositionType & position() const
std::vector< Alignable * > theAlignables
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(), edm::hlt::Exception, edm::ParameterSet::getParameter(), i, AlignmentParameterSelector::selectedAlignables(), AlCaHLTBitMon_QueryRunRegistry::string, theAlignableDetAccessor, theAlignables, theAlignmentParameterStore, theAPEParameters, theAPEParameterSet, and theApplyAPE.

99 {
100  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Initializing...";
101 
102  // accessor Det->AlignableDet
103  if ( !muon )
105  else if ( !tracker )
107  else
108  theAlignableDetAccessor = new AlignableNavigator(tracker, muon);
109 
110  // set alignmentParameterStore
112 
113  // get alignables
115 
116  // clear theAPEParameters, if necessary
117  theAPEParameters.clear();
118 
119  // get APE parameters
120  if(theApplyAPE){
121  AlignmentParameterSelector selector(tracker, muon);
122  for (std::vector<edm::ParameterSet>::const_iterator setiter = theAPEParameterSet.begin();
123  setiter != theAPEParameterSet.end();
124  ++setiter) {
125  std::vector<Alignable*> alignables;
126 
127  selector.clear();
128  edm::ParameterSet selectorPSet = setiter->getParameter<edm::ParameterSet>("Selector");
129  std::vector<std::string> alignParams = selectorPSet.getParameter<std::vector<std::string> >("alignParams");
130  if (alignParams.size() == 1 && alignParams[0] == std::string("selected")) {
131  alignables = theAlignables;
132  }
133  else {
134  selector.addSelections(selectorPSet);
135  alignables = selector.selectedAlignables();
136  }
137 
138  std::vector<double> apeSPar = setiter->getParameter<std::vector<double> >("apeSPar");
139  std::vector<double> apeRPar = setiter->getParameter<std::vector<double> >("apeRPar");
140  std::string function = setiter->getParameter<std::string>("function");
141 
142  if (apeSPar.size() != 3 || apeRPar.size() != 3)
143  throw cms::Exception("BadConfig") << "apeSPar and apeRPar must have 3 values each" << std::endl;
144 
145  for (std::vector<double>::const_iterator i = apeRPar.begin(); i != apeRPar.end(); ++i) {
146  apeSPar.push_back(*i);
147  }
148 
149  if (function == std::string("linear")) {
150  apeSPar.push_back(0.); // c.f. note in calcAPE
151  }
152  else if (function == std::string("exponential")) {
153  apeSPar.push_back(1.); // c.f. note in calcAPE
154  }
155  else if (function == std::string("step")) {
156  apeSPar.push_back(2.); // c.f. note in calcAPE
157  }
158  else {
159  throw cms::Exception("BadConfig") << "APE function must be \"linear\" or \"exponential\"." << std::endl;
160  }
161 
162  theAPEParameters.push_back(std::pair<std::vector<Alignable*>, std::vector<double> >(alignables, apeSPar));
163  }
164  }
165 }
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
AlignmentParameterStore * theAlignmentParameterStore
std::vector< std::pair< std::vector< Alignable * >, std::vector< double > > > theAPEParameters
AlignableNavigator * theAlignableDetAccessor
std::vector< edm::ParameterSet > theAPEParameterSet
std::vector< Alignable * > theAlignables
const align::Alignables & alignables(void) const
get all alignables
bool HIPAlignmentAlgorithm::processHit1D ( const AlignableDetOrUnitPtr alidet,
const Alignable ali,
const TrajectoryStateOnSurface tsos,
const TrackingRecHit hit 
)
private

Definition at line 434 of file HIPAlignmentAlgorithm.cc.

References HIPUserVariables::alichi2, Alignable::alignmentParameters(), HIPUserVariables::alindof, HIPUserVariables::jtve, HIPUserVariables::jtvj, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localPosition(), TrackingRecHit::localPosition(), TrackingRecHit::localPositionError(), HIPUserVariables::nhit, LocalTrajectoryError::positionError(), AlignmentParameters::selectedDerivatives(), mathSSE::sqrt(), theMaxAllowedHitPull, AlignmentParameters::userVariables(), PV3DBase< T, PVType, FrameType >::x(), and LocalError::xx().

Referenced by run().

438 {
439  static const unsigned int hitDim = 1;
440 
441  // get trajectory impact point
442  LocalPoint alvec = tsos.localPosition();
443  AlgebraicVector pos(hitDim);
444  pos[0] = alvec.x();
445 
446  // get impact point covariance
447  AlgebraicSymMatrix ipcovmat(hitDim);
448  ipcovmat[0][0] = tsos.localError().positionError().xx();
449 
450  // get hit local position and covariance
451  AlgebraicVector coor(hitDim);
452  coor[0] = hit->localPosition().x();
453 
454  AlgebraicSymMatrix covmat(hitDim);
455  covmat[0][0] = hit->localPositionError().xx();
456 
457  // add hit and impact point covariance matrices
458  covmat = covmat + ipcovmat;
459 
460  // calculate the x pull of this hit
461  double xpull = 0.;
462  if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0]));
463 
464  // get Alignment Parameters
465  AlignmentParameters* params = ali->alignmentParameters();
466  // get derivatives
467  AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
468  // calculate user parameters
469  int npar = derivs2D.num_row();
470  // std::cout << "Dumping derivsSEL: \n" << derivs2D <<std::endl;
471  AlgebraicMatrix derivs(npar, hitDim, 0);
472 
473  for (int ipar=0;ipar<npar;ipar++) {
474  for (unsigned int idim=0;idim<hitDim;idim++) {
475  derivs[ipar][idim] = derivs2D[ipar][idim];
476  }
477  }
478  // std::cout << "Dumping derivs: \n" << derivs << std::endl;
479  // invert covariance matrix
480  int ierr;
481  // if (covmat[1][1]>4.0) nhitDim = hitDim;
482 
483  covmat.invert(ierr);
484  if (ierr != 0) {
485  edm::LogError("Alignment") << "Matrix inversion failed!";
486  return false;
487  }
488 
489  bool useThisHit = (theMaxAllowedHitPull <= 0.);
490 
491  useThisHit |= (fabs(xpull) < theMaxAllowedHitPull);
492 
493  // bailing out
494  if (!useThisHit) return false;
495 
496  // std::cout << "We use this hit in " << subDet << std::endl;
497  // std::cout << "Npars= " << npar << " NhitDim=1 Size of derivs: "
498  // << derivs.num_row() << " x " << derivs.num_col() << std::endl;
499 
500  // AlgebraicMatrix jtvjtmp();
501  AlgebraicMatrix covtmp(covmat);
502  // std::cout << "Preapring JTVJTMP -> " << std::flush;
503  AlgebraicMatrix jtvjtmp(derivs * covtmp *derivs.T());
504  // std::cout << "TMP JTVJ= \n" << jtvjtmp << std::endl;
505  AlgebraicSymMatrix thisjtvj(npar);
506  AlgebraicVector thisjtve(npar);
507  // std::cout << "Preparing ThisJtVJ" << std::flush;
508  // thisjtvj = covmat.similarity(derivs);
509  thisjtvj.assign(jtvjtmp);
510  // std::cout<<" ThisJtVE"<<std::endl;
511  thisjtve=derivs * covmat * (pos-coor);
512 
513  AlgebraicVector hitresidual(hitDim);
514  hitresidual[0] = (pos[0] - coor[0]);
515 
516  AlgebraicMatrix hitresidualT;
517  hitresidualT = hitresidual.T();
518  // std::cout << "HitResidualT = \n" << hitresidualT << std::endl;
519 
520  // access user variables (via AlignmentParameters)
521  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
522  uservar->jtvj += thisjtvj;
523  uservar->jtve += thisjtve;
524  uservar->nhit++;
525  // The following variable is needed for the extended 1D/2D hit fix using
526  // matrix shrinkage and expansion
527  // uservar->hitdim = hitDim;
528 
529  //for alignable chi squared
530  float thischi2;
531  thischi2 = (hitresidualT *covmat *hitresidual)[0];
532 
533  if ( verbose && ((thischi2/ static_cast<float>(uservar->nhit)) >10.0) ) {
534  edm::LogWarning("Alignment") << "Added to Chi2 the number " << thischi2 <<" having "
535  << uservar->nhit << " dof " << std::endl << "X-resid "
536  << hitresidual[0] << " Y-resid "
537  << hitresidual[1] << std::endl << " Cov^-1 matr (covmat): [0][0]= "
538  << covmat[0][0] << " [0][1]= "
539  << covmat[0][1] << " [1][0]= "
540  << covmat[1][0] << " [1][1]= "
541  << covmat[1][1] << std::endl;
542  }
543 
544  uservar->alichi2 += thischi2; // a bit weird, but vector.transposed * matrix * vector doesn't give a double in CMSSW's opinion
545  uservar->alindof += hitDim; // 2D hits contribute twice to the ndofs
546 
547  return true;
548 }
float xx() const
Definition: LocalError.h:24
AlgebraicVector jtve
LocalError positionError() const
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
AlgebraicSymMatrix jtvj
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
CLHEP::HepMatrix AlgebraicMatrix
T sqrt(T t)
Definition: SSEVec.h:48
const LocalTrajectoryError & localError() const
CLHEP::HepVector AlgebraicVector
virtual LocalError localPositionError() const =0
CLHEP::HepSymMatrix AlgebraicSymMatrix
T x() const
Definition: PV3DBase.h:62
virtual LocalPoint localPosition() const =0
virtual AlgebraicMatrix selectedDerivatives(const TrajectoryStateOnSurface &tsos, const AlignableDetOrUnitPtr &alidet) const
bool HIPAlignmentAlgorithm::processHit2D ( const AlignableDetOrUnitPtr alidet,
const Alignable ali,
const TrajectoryStateOnSurface tsos,
const TrackingRecHit hit 
)
private

Definition at line 550 of file HIPAlignmentAlgorithm.cc.

References HIPUserVariables::alichi2, Alignable::alignmentParameters(), HIPUserVariables::alindof, HIPUserVariables::jtve, HIPUserVariables::jtvj, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localPosition(), TrackingRecHit::localPosition(), TrackingRecHit::localPositionError(), HIPUserVariables::nhit, 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().

554 {
555  static const unsigned int hitDim = 2;
556 
557  // get trajectory impact point
558  LocalPoint alvec = tsos.localPosition();
559  AlgebraicVector pos(hitDim);
560  pos[0] = alvec.x();
561  pos[1] = alvec.y();
562 
563  // get impact point covariance
564  AlgebraicSymMatrix ipcovmat(hitDim);
565  ipcovmat[0][0] = tsos.localError().positionError().xx();
566  ipcovmat[1][1] = tsos.localError().positionError().yy();
567  ipcovmat[0][1] = tsos.localError().positionError().xy();
568 
569  // get hit local position and covariance
570  AlgebraicVector coor(hitDim);
571  coor[0] = hit->localPosition().x();
572  coor[1] = hit->localPosition().y();
573 
574  AlgebraicSymMatrix covmat(hitDim);
575  covmat[0][0] = hit->localPositionError().xx();
576  covmat[1][1] = hit->localPositionError().yy();
577  covmat[0][1] = hit->localPositionError().xy();
578 
579  // add hit and impact point covariance matrices
580  covmat = covmat + ipcovmat;
581 
582  // calculate the x pull and y pull of this hit
583  double xpull = 0.;
584  double ypull = 0.;
585  if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0]));
586  if (covmat[1][1] != 0.) ypull = (pos[1] - coor[1])/sqrt(fabs(covmat[1][1]));
587 
588  // get Alignment Parameters
589  AlignmentParameters* params = ali->alignmentParameters();
590  // get derivatives
591  AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
592  // calculate user parameters
593  int npar = derivs2D.num_row();
594  // std::cout << "Dumping derivsSEL: \n"<< derivs2D <<std::endl;
595  AlgebraicMatrix derivs(npar, hitDim, 0);
596 
597  for (int ipar=0;ipar<npar;ipar++) {
598  for (unsigned int idim=0;idim<hitDim;idim++) {
599  derivs[ipar][idim] = derivs2D[ipar][idim];
600  }
601  }
602  // std::cout << "Dumping derivs: \n" << derivs << std::endl;
603  // invert covariance matrix
604  int ierr;
605  // if (covmat[1][1]>4.0) nhitDim = hitDim;
606 
607  covmat.invert(ierr);
608  if (ierr != 0) {
609  edm::LogError("Alignment") << "Matrix inversion failed!";
610  return false;
611  }
612 
613  bool useThisHit = (theMaxAllowedHitPull <= 0.);
614 
615  useThisHit |= (fabs(xpull) < theMaxAllowedHitPull && fabs(ypull) < theMaxAllowedHitPull);
616 
617  // bailing out
618  if (!useThisHit) return false;
619 
620  // std::cout << "We use this hit in " << subDet << std::endl;
621  // std::cout << "Npars= " << npar << " NhitDim=2 Size of derivs: "
622  // << derivs.num_row() << " x " << derivs.num_col() << std::endl;
623 
624  // AlgebraicMatrix jtvjtmp();
625  AlgebraicMatrix covtmp(covmat);
626  // std::cout << "Preapring JTVJTMP -> " << std::flush;
627  AlgebraicMatrix jtvjtmp(derivs * covtmp *derivs.T());
628  // std::cout << "TMP JTVJ= \n" << jtvjtmp << std::endl;
629  AlgebraicSymMatrix thisjtvj(npar);
630  AlgebraicVector thisjtve(npar);
631  // std::cout << "Preparing ThisJtVJ" << std::flush;
632  // thisjtvj = covmat.similarity(derivs);
633  thisjtvj.assign(jtvjtmp);
634  // std::cout<<" ThisJtVE"<<std::endl;
635  thisjtve=derivs * covmat * (pos-coor);
636 
637  AlgebraicVector hitresidual(hitDim);
638  hitresidual[0] = (pos[0] - coor[0]);
639  hitresidual[1] = (pos[1] - coor[1]);
640  // if(nhitDim>1) {
641  // hitresidual[1] =0.0;
642  // nhitDim=1;
643  // }
644 
645  AlgebraicMatrix hitresidualT;
646  hitresidualT = hitresidual.T();
647  // std::cout << "HitResidualT = \n" << hitresidualT << std::endl;
648  // access user variables (via AlignmentParameters)
649  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
650  uservar->jtvj += thisjtvj;
651  uservar->jtve += thisjtve;
652  uservar->nhit++;
653  // The following variable is needed for the extended 1D/2D hit fix using
654  // matrix shrinkage and expansion
655  // uservar->hitdim = hitDim;
656 
657  //for alignable chi squared
658  float thischi2;
659  thischi2 = (hitresidualT *covmat *hitresidual)[0];
660 
661  if ( verbose && ((thischi2/ static_cast<float>(uservar->nhit)) >10.0) ) {
662  edm::LogWarning("Alignment") << "Added to Chi2 the number " << thischi2 <<" having "
663  << uservar->nhit << " dof " << std::endl << "X-resid "
664  << hitresidual[0] << " Y-resid "
665  << hitresidual[1] << std::endl << " Cov^-1 matr (covmat): [0][0]= "
666  << covmat[0][0] << " [0][1]= "
667  << covmat[0][1] << " [1][0]= "
668  << covmat[1][0] << " [1][1]= "
669  << covmat[1][1] << std::endl;
670  }
671 
672  uservar->alichi2 += thischi2; // a bit weird, but vector.transposed * matrix * vector doesn't give a double in CMSSW's opinion
673  uservar->alindof += hitDim; // 2D hits contribute twice to the ndofs
674 
675  return true;
676 }
float xx() const
Definition: LocalError.h:24
T y() const
Definition: PV3DBase.h:63
AlgebraicVector jtve
LocalError positionError() const
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
AlgebraicSymMatrix jtvj
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
float xy() const
Definition: LocalError.h:25
CLHEP::HepMatrix AlgebraicMatrix
float yy() const
Definition: LocalError.h:26
T sqrt(T t)
Definition: SSEVec.h:48
const LocalTrajectoryError & localError() const
CLHEP::HepVector AlgebraicVector
virtual LocalError localPositionError() const =0
CLHEP::HepSymMatrix AlgebraicSymMatrix
T x() const
Definition: PV3DBase.h:62
virtual LocalPoint localPosition() const =0
virtual AlgebraicMatrix selectedDerivatives(const TrajectoryStateOnSurface &tsos, const AlignableDetOrUnitPtr &alidet) const
int HIPAlignmentAlgorithm::readIterationFile ( std::string  filename)
private

Definition at line 909 of file HIPAlignmentAlgorithm.cc.

References recoMuon::in, and query::result.

Referenced by startNewLoop().

910 {
911  int result;
912 
913  std::ifstream inIterFile(filename.c_str(), std::ios::in);
914  if (!inIterFile) {
915  edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] ERROR! "
916  << "Unable to open Iteration file";
917  result = -1;
918  }
919  else {
920  inIterFile >> result;
921  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] "
922  << "Read last iteration number from file: " << result;
923  }
924  inIterFile.close();
925 
926  return result;
927 }
tuple result
Definition: query.py:137
tuple filename
Definition: lut2db_cfg.py:20
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(), reco::TrackBase::d0(), AlignableNavigator::detAndSubdetInMap(), reco::TrackBase::dz(), reco::TrackBase::eta(), eta(), edm::hlt::Exception, TrajectoryMeasurement::forwardPredictedState(), TrackingRecHit::geographicalId(), TrackingRecHit::hit(), reco::TrackBase::hitPattern(), isCollector, AlignmentClusterFlag::isTaken(), TrajectoryStateOnSurface::isValid(), TrackingRecHit::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::HitPattern::numberOfValidPixelBarrelHits(), reco::HitPattern::numberOfValidPixelEndcapHits(), reco::TrackBase::p(), AlCaHLTBitMon_ParallelJobs::p, phi, reco::TrackBase::phi(), processHit1D(), processHit2D(), RecoTauCleanerPlugins::pt, reco::TrackBase::pt(), TrajectoryMeasurement::recHit(), AlignmentParameterStore::selectParameters(), DetId::subdetId(), theAlignableDetAccessor, theAlignmentParameterStore, theCurrentPrescale, theEventPrescale, theFile, theFillTrackMonitoring, theTree, testEve_cfg::tracks, and AlignmentAlgorithmBase::EventInfo::trajTrackPairs().

Referenced by Types.EventID::cppID(), and Types.LuminosityBlockID::cppID().

680 {
681  if (isCollector) return;
682 
683  TrajectoryStateCombiner tsoscomb;
684 
685  // AM: not really needed
686  // AM: m_Ntracks = 0 should be sufficient
687  int itr=0;
688  m_Ntracks=0;
689  for(itr=0;itr<MAXREC;++itr){
690  m_Nhits[itr]=0;
691  m_Pt[itr]=-5.0;
692  m_P[itr]=-5.0;
693  m_nhPXB[itr]=0;
694  m_nhPXF[itr]=0;
695  m_Eta[itr]=-99.0;
696  m_Phi[itr]=-4.0;
697  m_Chi2n[itr]=-11.0;
698  m_d0[itr]=-999;
699  m_dz[itr]=-999;
700  }
701  itr=0;
702 
703  // AM: what is this needed for?
704  theFile->cd();
705 
706  // loop over tracks
707  const ConstTrajTrackPairCollection &tracks = eventInfo.trajTrackPairs();
708  for (ConstTrajTrackPairCollection::const_iterator it=tracks.begin();
709  it!=tracks.end();
710  ++it) {
711 
712  const Trajectory* traj = (*it).first;
713  const reco::Track* track = (*it).second;
714 
715  float pt = track->pt();
716  float eta = track->eta();
717  float phi = track->phi();
718  float p = track->p();
719  float chi2n = track->normalizedChi2();
720  int nhit = track->numberOfValidHits();
721  float d0 = track->d0();
722  float dz = track->dz();
723 
724  int nhpxb = track->hitPattern().numberOfValidPixelBarrelHits();
725  int nhpxf = track->hitPattern().numberOfValidPixelEndcapHits();
726 
727  if (verbose) edm::LogInfo("Alignment") << "New track pt,eta,phi,chi2n,hits: "
728  << pt << ","
729  << eta << ","
730  << phi << ","
731  << chi2n << ","
732  << nhit;
733  // edm::LogWarning("Alignment") << "New track pt,eta,phi,chi2n,hits: "
734  // << pt << ","
735  // << eta << ","
736  // << phi << ","
737  // << chi2n << ","
738  // << nhit;
739 
740  // fill track parameters in root tree
741  if (itr<MAXREC) {
742  m_Nhits[itr]=nhit;
743  m_Pt[itr]=pt;
744  m_P[itr]=p;
745  m_Eta[itr]=eta;
746  m_Phi[itr]=phi;
747  m_Chi2n[itr]=chi2n;
748  m_nhPXB[itr]=nhpxb;
749  m_nhPXF[itr]=nhpxf;
750  m_d0[itr]=d0;
751  m_dz[itr]=dz;
752  itr++;
753  m_Ntracks=itr;
754  }
755  // AM: Can be simplified
756 
757  std::vector<const TransientTrackingRecHit*> hitvec;
758  std::vector<TrajectoryStateOnSurface> tsosvec;
759 
760  // loop over measurements
761  std::vector<TrajectoryMeasurement> measurements = traj->measurements();
762  for (std::vector<TrajectoryMeasurement>::iterator im=measurements.begin();
763  im!=measurements.end();
764  ++im) {
765 
766  TrajectoryMeasurement meas = *im;
767 
768  // const TransientTrackingRecHit* ttrhit = &(*meas.recHit());
769  // const TrackingRecHit *hit = ttrhit->hit();
770  const TransientTrackingRecHit* hit = &(*meas.recHit());
771 
773 
774  // this is the updated state (including the current hit)
775  // TrajectoryStateOnSurface tsos=meas.updatedState();
776  // combine fwd and bwd predicted state to get state
777  // which excludes current hit
778 
780  bool skiphit = false;
781  if (eventInfo.clusterValueMap()) {
782  // check from the PrescalingMap if the hit was taken.
783  // If not skip to the next TM
784  // bool hitTaken=false;
785  AlignmentClusterFlag myflag;
786 
787  int subDet = hit->geographicalId().subdetId();
788  //take the actual RecHit out of the Transient one
789  const TrackingRecHit *rechit=hit->hit();
790  if (subDet>2) { // AM: if possible use enum instead of hard-coded value
791  const std::type_info &type = typeid(*rechit);
792 
793  if (type == typeid(SiStripRecHit1D)) {
794 
795  const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(rechit);
796  if (stripHit1D) {
797  SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
798  // myflag=PrescMap[stripclust];
799  myflag = (*eventInfo.clusterValueMap())[stripclust];
800  } else {
801  edm::LogError("HIPAlignmentAlgorithm")
802  << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
803  << "TypeId of the RecHit: " << className(*hit) <<std::endl;
804  }
805 
806  }//end if type = SiStripRecHit1D
807  else if(type == typeid(SiStripRecHit2D)){
808 
809  const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(rechit);
810  if (stripHit2D) {
811  SiStripRecHit2D::ClusterRef stripclust(stripHit2D->cluster());
812  // myflag=PrescMap[stripclust];
813  myflag = (*eventInfo.clusterValueMap())[stripclust];
814  } else {
815  edm::LogError("HIPAlignmentAlgorithm")
816  << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
817  // << "TypeId of the TTRH: " << className(*ttrhit) << std::endl;
818  << "TypeId of the TTRH: " << className(*hit) << std::endl;
819  }
820  } //end if type == SiStripRecHit2D
821  } //end if hit from strips
822  else {
823  const SiPixelRecHit* pixelhit= dynamic_cast<const SiPixelRecHit*>(rechit);
824  if (pixelhit) {
825  SiPixelClusterRefNew pixelclust(pixelhit->cluster());
826  // myflag=PrescMap[pixelclust];
827  myflag = (*eventInfo.clusterValueMap())[pixelclust];
828  }
829  else {
830  edm::LogError("HIPAlignmentAlgorithm")
831  << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Pixel RecHit failed! "
832  // << "TypeId of the TTRH: " << className(*ttrhit) << std::endl;
833  << "TypeId of the TTRH: " << className(*hit) << std::endl;
834  }
835  } //end 'else' it is a pixel hit
836  // bool hitTaken=myflag.isTaken();
837  if (!myflag.isTaken()) {
838  skiphit=true;
839  //std::cout<<"Hit from subdet "<<rechit->geographicalId().subdetId()<<" prescaled out."<<std::endl;
840  continue;
841  }
842  }//end if Prescaled Hits
844  if (skiphit) {
845  throw cms::Exception("LogicError")
846  << "ERROR in <HIPAlignmentAlgorithm::run>: this hit should have been skipped!"
847  << std::endl;
848  }
849 
851  meas.backwardPredictedState());
852 
853  if(tsos.isValid()){
854  // hitvec.push_back(ttrhit);
855  hitvec.push_back(hit);
856  tsosvec.push_back(tsos);
857  }
858 
859  } //hit valid
860  }
861 
862  // transform RecHit vector to AlignableDet vector
863  std::vector <AlignableDetOrUnitPtr> alidetvec = theAlignableDetAccessor->alignablesFromHits(hitvec);
864 
865  // get concatenated alignment parameters for list of alignables
867 
868  std::vector<TrajectoryStateOnSurface>::const_iterator itsos=tsosvec.begin();
869  std::vector<const TransientTrackingRecHit*>::const_iterator ihit=hitvec.begin();
870 
871  // loop over vectors(hit,tsos)
872  while (itsos != tsosvec.end()) {
873 
874  // get AlignableDet for this hit
875  const GeomDet* det = (*ihit)->det();
876  // int subDet= (*ihit)->geographicalId().subdetId();
877  uint32_t nhitDim = (*ihit)->dimension();
878 
880 
881  // get relevant Alignable
882  Alignable* ali = aap.alignableFromAlignableDet(alidet);
883 
884  if (ali!=0) {
885  if (nhitDim==1) {
886  processHit1D(alidet, ali, *itsos, *ihit);
887  } else if (nhitDim==2) {
888  processHit2D(alidet, ali, *itsos, *ihit);
889  }
890  }
891 
892  itsos++;
893  ihit++;
894  }
895  } // end of track loop
896 
897  // fill eventwise root tree (with prescale defined in pset)
900  if (theCurrentPrescale<=0) {
901  theTree->Fill();
903  }
904  }
905 }
ClusterRef cluster() const
double p() const
momentum vector magnitude
Definition: TrackBase.h:127
type
Definition: HCALResponse.h:21
ConstRecHitPointer const & recHit() const
double d0() const
dxy parameter in perigee convention (d0 = - dxy)
Definition: TrackBase.h:121
AlignableDetOrUnitPtr alignableFromGeomDet(const GeomDet *geomDet)
Returns AlignableDetOrUnitPtr corresponding to given GeomDet.
AlignmentParameterStore * theAlignmentParameterStore
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:109
TSOS combine(const TSOS &pTsos1, const TSOS &pTsos2) const
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:137
T eta() const
Alignable * alignableFromAlignableDet(const AlignableDetOrUnitPtr &adet) const
Get relevant Alignable from AlignableDet.
DataContainer const & measurements() const
Definition: Trajectory.h:203
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:139
int numberOfValidPixelBarrelHits() const
Definition: HitPattern.h:605
double pt() const
track transverse momentum
Definition: TrackBase.h:129
bool processHit2D(const AlignableDetOrUnitPtr &alidet, const Alignable *ali, const TrajectoryStateOnSurface &tsos, const TrackingRecHit *hit)
AlignableNavigator * theAlignableDetAccessor
ClusterRef cluster() const
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:230
bool processHit1D(const AlignableDetOrUnitPtr &alidet, const Alignable *ali, const TrajectoryStateOnSurface &tsos, const TrackingRecHit *hit)
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:221
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
TrajectoryStateOnSurface const & forwardPredictedState() const
Access to forward predicted state (from fitter or builder)
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:125
virtual TrackingRecHit const * hit() const
tuple tracks
Definition: testEve_cfg.py:39
ClusterRef cluster() const
Definition: SiPixelRecHit.h:49
CompositeAlignmentParameters selectParameters(const std::vector< AlignableDet * > &alignabledets) const
bool isValid() const
int numberOfValidPixelEndcapHits() const
Definition: HitPattern.h:609
std::vector< AlignableDetOrUnitPtr > alignablesFromHits(const std::vector< const TransientTrackingRecHit * > &hitvec)
Returns vector AlignableDetOrUnitPtr for given vector of Hits.
DetId geographicalId() const
bool detAndSubdetInMap(const DetId &detid) const
Given a DetId, returns true if DetIds with this detector and subdetector id are in the map (not neces...
std::vector< ConstTrajTrackPair > ConstTrajTrackPairCollection
std::string className(const T &t)
Definition: ClassName.h:30
TrajectoryStateOnSurface const & backwardPredictedState() const
Access to backward predicted state (from smoother)
Our base class.
Definition: SiPixelRecHit.h:23
Definition: DDAxes.h:10
void HIPAlignmentAlgorithm::setAlignmentPositionError ( void  )
private

Definition at line 948 of file HIPAlignmentAlgorithm.cc.

References calcAPE(), i, AlignmentParameterStore::setAlignmentPositionError(), theAlignmentParameterStore, theAPEParameters, theApplyAPE, and theIteration.

Referenced by startNewLoop().

949 {
950 
951 
952  // Check if user wants to override APE
953  if ( !theApplyAPE )
954  {
955  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] No APE applied";
956  return; // NO APE APPLIED
957  }
958 
959 
960  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] Apply APE!";
961 
962  double apeSPar[3], apeRPar[3];
963  for (std::vector<std::pair<std::vector<Alignable*>, std::vector<double> > >::const_iterator alipars = theAPEParameters.begin();
964  alipars != theAPEParameters.end();
965  ++alipars) {
966  const std::vector<Alignable*> &alignables = alipars->first;
967  const std::vector<double> &pars = alipars->second;
968 
969  apeSPar[0] = pars[0];
970  apeSPar[1] = pars[1];
971  apeSPar[2] = pars[2];
972  apeRPar[0] = pars[3];
973  apeRPar[1] = pars[4];
974  apeRPar[2] = pars[5];
975 
976  double function = pars[6];
977 
978  // Printout for debug
979  printf("APE: %u alignables\n", (unsigned int)alignables.size());
980  for ( int i=0; i<21; ++i ) {
981  double apelinstest=calcAPE(apeSPar,i,0.);
982  double apeexpstest=calcAPE(apeSPar,i,1.);
983  double apelinrtest=calcAPE(apeRPar,i,0.);
984  double apeexprtest=calcAPE(apeRPar,i,1.);
985  printf("APE: iter slin sexp rlin rexp: %5d %12.5f %12.5f %12.5f %12.5f\n",
986  i,apelinstest,apeexpstest,apelinrtest,apeexprtest);
987  }
988 
989  // set APE
990  double apeshift=calcAPE(apeSPar,theIteration,function);
991  double aperot =calcAPE(apeRPar,theIteration,function);
992  theAlignmentParameterStore->setAlignmentPositionError( alignables, apeshift, aperot );
993  }
994 
995 }
int i
Definition: DBlmapReader.cc:9
AlignmentParameterStore * theAlignmentParameterStore
void setAlignmentPositionError(const align::Alignables &alivec, double valshift, double valrot)
Set Alignment position error.
std::vector< std::pair< std::vector< Alignable * >, std::vector< double > > > theAPEParameters
double calcAPE(double *par, int iter, double function)
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().

169 {
170 
171  // iterate over all alignables and attach user variables
172  for( std::vector<Alignable*>::const_iterator it=theAlignables.begin();
173  it!=theAlignables.end(); it++ )
174  {
175  AlignmentParameters* ap = (*it)->alignmentParameters();
176  int npar=ap->numSelected();
177  HIPUserVariables* userpar = new HIPUserVariables(npar);
178  ap->setUserVariables(userpar);
179  }
180 
181  // try to read in alignment parameters from a previous iteration
182  AlignablePositions theAlignablePositionsFromFile =
184  salignedfile.c_str(),-1,ioerr);
185 
186  int numAlignablesFromFile = theAlignablePositionsFromFile.size();
187 
188  if (numAlignablesFromFile==0) { // file not there: first iteration
189 
190  // set iteration number to 1
191  if (isCollector) theIteration=0;
192  else theIteration=1;
193  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] File not found => iteration "<<theIteration;
194 
195  // get true (de-misaligned positions) and write to root file
196  // hardcoded iteration=1
198  struefile.c_str(),1,false,ioerr);
199 
200  // get misaligned positions and write to root file
201  // hardcoded iteration=1
203  smisalignedfile.c_str(),1,false,ioerr);
204 
205  }
206 
207  else { // there have been previous iterations
208 
209  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Alignables Read "
210  << numAlignablesFromFile;
211 
212  // get iteration number from file
214 
215  // increase iteration
216  theIteration++;
217  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Iteration increased by one!";
218 
219  // now apply psotions of file from prev iteration
220  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Apply positions from file ...";
222  theAlignablePositionsFromFile,ioerr);
223 
224  }
225 
226  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Current Iteration number: "
227  << theIteration;
228 
229 
230  // book root trees
231  bookRoot();
232 
233 
234  /*---------------------moved to terminate------------------------------
235  if (theLevels.size() > 0)
236  {
237  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
238 
239  unsigned int nAlignable = theAlignables.size();
240 
241  for (unsigned int i = 0; i < nAlignable; ++i)
242  {
243  const Alignable* ali = theAlignables[i];
244 
245  AlignmentParameters* ap = ali->alignmentParameters();
246 
247  HIPUserVariables* uservar =
248  dynamic_cast<HIPUserVariables*>(ap->userVariables());
249 
250  for (unsigned int l = 0; l < theLevels.size(); ++l)
251  {
252  SurveyResidual res(*ali, theLevels[l], true);
253 
254  if ( res.valid() )
255  {
256  AlgebraicSymMatrix invCov = res.inverseCovariance();
257 
258  // variable for tree
259  AlgebraicVector sensResid = res.sensorResidual();
260  m3_Id = ali->id();
261  m3_ObjId = theLevels[l];
262  m3_par[0] = sensResid[0]; m3_par[1] = sensResid[1]; m3_par[2] = sensResid[2];
263  m3_par[3] = sensResid[3]; m3_par[4] = sensResid[4]; m3_par[5] = sensResid[5];
264 
265  uservar->jtvj += invCov;
266  uservar->jtve += invCov * sensResid;
267 
268  theTree3->Fill();
269  }
270  }
271 
272  // align::LocalVectors residuals = res1.pointsResidual();
273 
274  // unsigned int nPoints = residuals.size();
275 
276  // for (unsigned int k = 0; k < nPoints; ++k)
277  // {
278  // AlgebraicMatrix J = term->survey()->derivatives(k);
279  // AlgebraicVector e(3); // local residual
280 
281  // const align::LocalVector& lr = residuals[k];
282 
283  // e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z();
284 
285  // uservar->jtvj += invCov1.similarity(J);
286  // uservar->jtve += J * (invCov1 * e);
287  // }
288 
289  }
290  }
291  //------------------------------------------------------*/
292 
293  // set alignment position error
295 
296  // run collector job if we are in parallel mode
297  if (isCollector) collector();
298 
299 }
AlignmentParameterStore * theAlignmentParameterStore
void writeAlignableOriginalPositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write Alignable original (before misalignment) absolute positions
void applyAlignableAbsolutePositions(const align::Alignables &alis, const AlignablePositions &newpos, int &ierr)
apply absolute positions to alignables
AlignablePositions readAlignableAbsolutePositions(const align::Alignables &alivec, const char *filename, int iter, int &ierr)
read Alignable current absolute positions
int readIterationFile(std::string filename)
int numSelected(void) const
Get number of selected parameters.
void setUserVariables(AlignmentUserVariables *auv)
Set pointer to user variables.
std::vector< AlignableAbsData > AlignablePositions
Definition: AlignableData.h:51
void writeAlignableAbsolutePositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write Alignable current absolute positions
std::vector< Alignable * > theAlignables
void HIPAlignmentAlgorithm::terminate ( const edm::EventSetup setup)
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, run_regression::test, theAlignables, theAlignmentParameterStore, theFile, theFile2, theFile3, theIO, theIteration, theLevels, theTree, theTree2, theTree3, AlignmentParameters::userVariables(), SurveyResidual::valid(), AlignmentIORoot::writeAlignableAbsolutePositions(), AlignmentIORoot::writeAlignmentParameters(), HIPUserVariablesIORoot::writeHIPUserVariables(), and writeIterationFile().

304 {
305 
306  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Terminating";
307 
308  // calculating survey residuals
309  if (theLevels.size() > 0)
310  {
311  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
312 
313  unsigned int nAlignable = theAlignables.size();
314 
315  for (unsigned int i = 0; i < nAlignable; ++i)
316  {
317  const Alignable* ali = theAlignables[i];
318 
320 
321  HIPUserVariables* uservar =
322  dynamic_cast<HIPUserVariables*>(ap->userVariables());
323 
324  for (unsigned int l = 0; l < theLevels.size(); ++l)
325  {
326  SurveyResidual res(*ali, theLevels[l], true);
327 
328  if ( res.valid() )
329  {
330  AlgebraicSymMatrix invCov = res.inverseCovariance();
331 
332  // variable for tree
333  AlgebraicVector sensResid = res.sensorResidual();
334  m3_Id = ali->id();
335  m3_ObjId = theLevels[l];
336  m3_par[0] = sensResid[0]; m3_par[1] = sensResid[1]; m3_par[2] = sensResid[2];
337  m3_par[3] = sensResid[3]; m3_par[4] = sensResid[4]; m3_par[5] = sensResid[5];
338 
339  uservar->jtvj += invCov;
340  uservar->jtve += invCov * sensResid;
341 
342  theTree3->Fill();
343  }
344  }
345 
346  // align::LocalVectors residuals = res1.pointsResidual();
347 
348  // unsigned int nPoints = residuals.size();
349 
350  // for (unsigned int k = 0; k < nPoints; ++k)
351  // {
352  // AlgebraicMatrix J = term->survey()->derivatives(k);
353  // AlgebraicVector e(3); // local residual
354 
355  // const align::LocalVector& lr = residuals[k];
356 
357  // e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z();
358 
359  // uservar->jtvj += invCov1.similarity(J);
360  // uservar->jtve += J * (invCov1 * e);
361  // }
362 
363  }
364  }
365 
366  // write user variables
369  theIteration,false,ioerr);
370 
371  // now calculate alignment corrections ...
372  int ialigned=0;
373  // iterate over alignment parameters
374  for(std::vector<Alignable*>::const_iterator
375  it=theAlignables.begin(); it!=theAlignables.end(); it++) {
376  Alignable* ali=(*it);
377  // Alignment parameters
379 
380  // try to calculate parameters
381  bool test = calcParameters(ali);
382 
383  // if successful, apply parameters
384  if (test) {
385  edm::LogInfo("Alignment") << "now apply params";
387  // set these parameters 'valid'
388  ali->alignmentParameters()->setValid(true);
389  // increase counter
390  ialigned++;
391  }
392  else par->setValid(false);
393  }
394  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::terminate] Aligned units: " << ialigned;
395 
396  // fill alignable wise root tree
397  fillRoot(iSetup);
398 
399  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Writing aligned parameters to file: " << theAlignables.size();
400 
401  // write new absolute positions to disk
403  salignedfile.c_str(),theIteration,false,ioerr);
404 
405  // write alignment parameters to disk
407  sparameterfile.c_str(),theIteration,false,ioerr);
408 
409  // write iteration number to file
411 
412  // write out trees and close root file
413 
414  // eventwise tree
415  theFile->cd();
416  theTree->Write();
417  delete theFile;
418 
419  if (theLevels.size() > 0){
420  theFile3->cd();
421  theTree3->Write();
422  delete theFile3;
423  }
424 
425  // alignable-wise tree is only filled once
426  if (theIteration==1) { // only for 1st iteration
427  theFile2->cd();
428  theTree2->Write();
429  delete theFile2;
430  }
431 
432 }
align::ID id() const
Return the ID of Alignable, i.e. DetId of &#39;first&#39; component GeomDet(Unit).
Definition: Alignable.h:180
int i
Definition: DBlmapReader.cc:9
AlignmentParameterStore * theAlignmentParameterStore
bool calcParameters(Alignable *ali)
void writeIterationFile(std::string filename, int iter)
void applyParameters(void)
Obsolete: Use AlignableNavigator::alignableDetFromDetId and alignableFromAlignableDet.
AlgebraicVector jtve
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
void writeHIPUserVariables(const Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
AlgebraicSymMatrix jtvj
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
void setValid(bool v)
Set validity flag.
void fillRoot(const edm::EventSetup &setup)
CLHEP::HepVector AlgebraicVector
void writeAlignmentParameters(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write AlignmentParameters
align::StructureType m3_ObjId
CLHEP::HepSymMatrix AlgebraicSymMatrix
void writeAlignableAbsolutePositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write Alignable current absolute positions
std::vector< align::StructureType > theLevels
std::vector< Alignable * > theAlignables
void HIPAlignmentAlgorithm::writeIterationFile ( std::string  filename,
int  iter 
)
private

Definition at line 931 of file HIPAlignmentAlgorithm.cc.

References getDQMSummary::iter, and dbtoconf::out.

Referenced by terminate().

932 {
933  std::ofstream outIterFile((filename.c_str()), std::ios::out);
934  if (!outIterFile) {
935  edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::writeIterationFile] ERROR: Unable to write Iteration file";
936  }
937  else {
938  outIterFile << iter;
939  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::writeIterationFile] writing iteration number to file: " << iter;
940  }
941  outIterFile.close();
942 }
tuple out
Definition: dbtoconf.py:99
tuple filename
Definition: lut2db_cfg.py:20

Member Data Documentation

int HIPAlignmentAlgorithm::ioerr
private

Definition at line 76 of file HIPAlignmentAlgorithm.h.

Referenced by collector(), startNewLoop(), and terminate().

bool HIPAlignmentAlgorithm::isCollector
private

Definition at line 98 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), run(), and startNewLoop().

float HIPAlignmentAlgorithm::m2_Eta
private

Definition at line 122 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

align::ID HIPAlignmentAlgorithm::m2_Id
private

Definition at line 123 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

int HIPAlignmentAlgorithm::m2_Layer
private

Definition at line 121 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

int HIPAlignmentAlgorithm::m2_Nhit
private

Definition at line 121 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

align::StructureType HIPAlignmentAlgorithm::m2_ObjId
private

Definition at line 124 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

float HIPAlignmentAlgorithm::m2_Phi
private

Definition at line 122 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

int HIPAlignmentAlgorithm::m2_Type
private

Definition at line 121 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

float HIPAlignmentAlgorithm::m2_Xpos
private

Definition at line 122 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

float HIPAlignmentAlgorithm::m2_Ypos
private

Definition at line 122 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

float HIPAlignmentAlgorithm::m2_Zpos
private

Definition at line 122 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and fillRoot().

align::ID HIPAlignmentAlgorithm::m3_Id
private

Definition at line 127 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and terminate().

align::StructureType HIPAlignmentAlgorithm::m3_ObjId
private

Definition at line 128 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and terminate().

float HIPAlignmentAlgorithm::m3_par[6]
private

Definition at line 129 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and terminate().

float HIPAlignmentAlgorithm::m_Chi2n[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

float HIPAlignmentAlgorithm::m_d0[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

float HIPAlignmentAlgorithm::m_dz[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

float HIPAlignmentAlgorithm::m_Eta[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

int HIPAlignmentAlgorithm::m_Nhits[MAXREC]
private

Definition at line 117 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

int HIPAlignmentAlgorithm::m_nhPXB[MAXREC]
private

Definition at line 117 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

int HIPAlignmentAlgorithm::m_nhPXF[MAXREC]
private

Definition at line 117 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

int HIPAlignmentAlgorithm::m_Ntracks
private

Definition at line 117 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

float HIPAlignmentAlgorithm::m_P[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

float HIPAlignmentAlgorithm::m_Phi[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

float HIPAlignmentAlgorithm::m_Pt[MAXREC]
private

Definition at line 118 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), and run().

const int HIPAlignmentAlgorithm::MAXREC = 99
staticprivate

Definition at line 115 of file HIPAlignmentAlgorithm.h.

Referenced by fillEventwiseTree(), and run().

std::string HIPAlignmentAlgorithm::outfile
private

Definition at line 84 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and HIPAlignmentAlgorithm().

std::string HIPAlignmentAlgorithm::outfile2
private

Definition at line 84 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and HIPAlignmentAlgorithm().

std::string HIPAlignmentAlgorithm::outpath
private

Definition at line 84 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm().

std::string HIPAlignmentAlgorithm::salignedfile
private

Definition at line 85 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), startNewLoop(), and terminate().

std::string HIPAlignmentAlgorithm::siterationfile
private

Definition at line 85 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), startNewLoop(), and terminate().

std::string HIPAlignmentAlgorithm::smisalignedfile
private

Definition at line 85 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and startNewLoop().

std::string HIPAlignmentAlgorithm::sparameterfile
private

Definition at line 84 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and terminate().

std::string HIPAlignmentAlgorithm::ssurveyfile
private

Definition at line 85 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and HIPAlignmentAlgorithm().

std::string HIPAlignmentAlgorithm::struefile
private

Definition at line 85 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and startNewLoop().

std::string HIPAlignmentAlgorithm::suvarfile
private

Definition at line 84 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and terminate().

AlignableNavigator* HIPAlignmentAlgorithm::theAlignableDetAccessor
private

Definition at line 73 of file HIPAlignmentAlgorithm.h.

Referenced by initialize(), and run().

std::vector<Alignable*> HIPAlignmentAlgorithm::theAlignables
private

Definition at line 72 of file HIPAlignmentAlgorithm.h.

Referenced by collector(), fillRoot(), initialize(), startNewLoop(), and terminate().

AlignmentParameterStore* HIPAlignmentAlgorithm::theAlignmentParameterStore
private
std::vector<std::pair<std::vector<Alignable*>, std::vector<double> > > HIPAlignmentAlgorithm::theAPEParameters
private

Definition at line 90 of file HIPAlignmentAlgorithm.h.

Referenced by initialize(), and setAlignmentPositionError().

std::vector<edm::ParameterSet> HIPAlignmentAlgorithm::theAPEParameterSet
private

Definition at line 89 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and initialize().

bool HIPAlignmentAlgorithm::theApplyAPE
private
int HIPAlignmentAlgorithm::theCollectorNJobs
private

Definition at line 99 of file HIPAlignmentAlgorithm.h.

Referenced by collector(), and HIPAlignmentAlgorithm().

std::string HIPAlignmentAlgorithm::theCollectorPath
private

Definition at line 100 of file HIPAlignmentAlgorithm.h.

Referenced by collector(), and HIPAlignmentAlgorithm().

int HIPAlignmentAlgorithm::theCurrentPrescale
private

Definition at line 101 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and run().

int HIPAlignmentAlgorithm::theEventPrescale
private

Definition at line 101 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), and run().

TFile* HIPAlignmentAlgorithm::theFile
private

Definition at line 107 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), run(), and terminate().

TFile* HIPAlignmentAlgorithm::theFile2
private

Definition at line 109 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillRoot(), and terminate().

TFile* HIPAlignmentAlgorithm::theFile3
private

Definition at line 111 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and terminate().

bool HIPAlignmentAlgorithm::theFillTrackMonitoring
private

Definition at line 102 of file HIPAlignmentAlgorithm.h.

Referenced by collector(), HIPAlignmentAlgorithm(), and run().

AlignmentIORoot HIPAlignmentAlgorithm::theIO
private

Definition at line 75 of file HIPAlignmentAlgorithm.h.

Referenced by startNewLoop(), and terminate().

int HIPAlignmentAlgorithm::theIteration
private
std::vector<align::StructureType> HIPAlignmentAlgorithm::theLevels
private

Definition at line 104 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), HIPAlignmentAlgorithm(), and terminate().

double HIPAlignmentAlgorithm::theMaxAllowedHitPull
private

Definition at line 92 of file HIPAlignmentAlgorithm.h.

Referenced by HIPAlignmentAlgorithm(), processHit1D(), and processHit2D().

double HIPAlignmentAlgorithm::theMaxRelParameterError
private

Definition at line 96 of file HIPAlignmentAlgorithm.h.

Referenced by calcParameters(), and HIPAlignmentAlgorithm().

int HIPAlignmentAlgorithm::theMinimumNumberOfHits
private

Definition at line 94 of file HIPAlignmentAlgorithm.h.

Referenced by calcParameters(), and HIPAlignmentAlgorithm().

TTree* HIPAlignmentAlgorithm::theTree
private

Definition at line 108 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillEventwiseTree(), run(), and terminate().

TTree* HIPAlignmentAlgorithm::theTree2
private

Definition at line 110 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), fillRoot(), and terminate().

TTree* HIPAlignmentAlgorithm::theTree3
private

Definition at line 112 of file HIPAlignmentAlgorithm.h.

Referenced by bookRoot(), and terminate().

bool HIPAlignmentAlgorithm::verbose
private