45 edm::LogInfo(
"Alignment") <<
"[MuonDTLocalMillepedeAlgorithm] constructed.";
66 if(workingmode == 0) {
67 edm::LogInfo(
"Alignment") <<
"[MuonDTLocalMillepedeAlgorithm] Running on production mode.";
69 snprintf(nameOfFile,
sizeof(nameOfFile),
"%s/MyNtupleResidual.root",
ntuplePath.c_str());
70 f =
new TFile(nameOfFile,
"RECREATE");
73 }
else if (workingmode == 1) {
74 edm::LogInfo(
"Alignment") <<
"[MuonDTLocalMillepedeAlgorithm] Running SLToSL algorithm.";
76 edm::LogInfo(
"Alignment") <<
"[MuonDTLocalMillepedeAlgorithm] Running Local Millepede algorithm.";
115 edm::LogInfo(
"Alignment") <<
"[MuonDTLocalMillepedeAlgorithm] Starting SLToSL algorithm";
118 edm::LogInfo(
"Alignment") <<
"[MuonDTLocalMillepedeAlgorithm] Starting local MuonMillepede algorithm";
140 for( ConstTrajTrackPairCollection::const_iterator it=tracks.begin();
141 it!=tracks.end();it++) {
153 if(pt < ptMin || pt >
ptMax)
continue;
155 vector<const TransientTrackingRecHit*> hitvec;
156 vector<TrajectoryMeasurement> measurements = traj->
measurements();
160 for (vector<TrajectoryMeasurement>::iterator im=measurements.begin();
161 im!=measurements.end(); im++) {
167 hitvec.push_back(hit);
174 vector<const TransientTrackingRecHit*>::const_iterator ihit=hitvec.begin();
178 while (ihit != hitvec.end())
180 const GeomDet* det=(*ihit)->det();
184 DTChamberId mChamber(mLayer.wheel(), mLayer.station(), mLayer.sector());
191 myTrack1D.
erx[ch_counter] = (*ihit)->localPositionError().xx();
216 bool saveThis =
false;
220 int numlayer[20][12];
221 for(
int s = 0;
s < 20; ++
s) {
222 for(
int k = 0;
k < 5; ++
k)
id[
s][
k] = 0;
223 for(
int k = 0;
k < 12; ++
k) numlayer[
s][
k] = 0;
230 for(
int counterCham = 0; counterCham < nChambers; counterCham++) {
235 id[counterCham][4]++;
237 id[counterCham][3]++;
239 for (
int ila = 1; ila<=4; ila++)
242 numlayer[counterCham][jla]++;
256 for (
int ila = 1; ila<=4; ila++)
259 numlayer[nChambers][jla]++;
267 xc[iseg][ihit] = -250.;
268 yc[iseg][ihit] = -250.;
269 zc[iseg][ihit] = -250.;
270 ex[iseg][ihit] = -250.;
271 xcp[iseg][ihit] = -250.;
272 ycp[iseg][ihit] = -250.;
273 excp[iseg][ihit] = -250.;
274 eycp[iseg][ihit] = -250.;
282 bool GoodPhiChamber =
true, GoodThetaChamber =
true;
283 for (
int ila = 1; ila<=12; ila++) {
284 if (numlayer[
counter][ila-1]!=1 && (ila<5 || ila>8)) GoodPhiChamber =
false;
285 if (numlayer[
counter][ila-1]!=1 && (ila<9 || ila>4) &&
id[
counter][1]!=4) GoodThetaChamber =
false;
289 GoodPhiChamber && GoodThetaChamber) {
291 TMatrixD phiProjection(2,2);
292 TMatrixD thetaProjection(2,2);
293 TMatrixD bphiProjection(2,1);
294 TMatrixD bthetaProjection(2,1);
296 TMatrixD phiProjectionSL1(2,2);
297 TMatrixD bphiProjectionSL1(2,1);
298 TMatrixD phiProjectionSL3(2,2);
299 TMatrixD bphiProjectionSL3(2,1);
304 int numh1 = 0, numh2 = 0, numh3 = 0;
349 if (phiProjection(0,0) != 0 && phiProjectionSL1(0,0) != 0 && phiProjectionSL3(0,0) != 0 &&
350 (thetaProjection(0,0) != 0 ||
id[
counter][1] == 4)) {
356 if(thetaProjection(0,0) != 0 &&
id[
counter][1] != 4) {
357 thetaProjection.Invert();
358 TMatrixD solution = thetaProjection*bthetaProjection;
365 phiProjection.Invert();
366 phiProjectionSL1.Invert();
367 phiProjectionSL3.Invert();
368 TMatrixD solution = phiProjection*bphiProjection;
369 TMatrixD solutionSL1 = phiProjectionSL1*bphiProjectionSL1;
370 TMatrixD solutionSL3 = phiProjectionSL3*bphiProjectionSL3;
381 xSL1SL3[
nseg] = solutionSL1(0,0) + SL3_z_ave * solutionSL1(1,0);
387 xSL3SL1[
nseg] = solutionSL3(0,0) + SL1_z_ave * solutionSL3(1,0);
double p() const
momentum vector magnitude
T getParameter(std::string const &) const
ConstRecHitPointer const & recHit() const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
double phi() const
azimuthal angle of momentum vector
std::vector< Alignable * > theAlignables
const ConstTrajTrackPairCollection & trajTrackPairs() const
define event information passed to algorithms
uint32_t rawId() const
get the raw id
DataContainer const & measurements() const
double eta() const
pseudorapidity of momentum vector
const GeomDet * det() const
void initialize(const edm::EventSetup &setup, AlignableTracker *tracker, AlignableMuon *muon, AlignmentParameterStore *store)
Call at beginning of job.
double pt() const
track transverse momentum
DetId geographicalId() const
The label of this GeomDet.
MuonDTLocalMillepedeAlgorithm(const edm::ParameterSet &cfg)
Constructor.
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
void run(const edm::EventSetup &setup, const EventInfo &eventInfo)
Run the algorithm on trajectories and tracks.
static std::atomic< unsigned int > counter
int charge() const
track electric charge
#define DEFINE_EDM_PLUGIN(factory, type, name)
Detector det() const
get the detector field from this detid
AlignableNavigator * theAlignableDetAccessor
Constructor of the full muon geometry.
AlignmentParameterStore * theAlignmentParameterStore
const align::Alignables & alignables(void) const
get all alignables
std::vector< ConstTrajTrackPair > ConstTrajTrackPairCollection
void terminate(void)
Call at end of job.