CommonAlignmentAlgorithm/plugins/ApeSettingAlgorithm.cc

 /**
 * \file MillePedeAlignmentAlgorithm.cc
 *
 *  \author    : Gero Flucke/Ivan Reid
 *  date       : February 2009 *  $Revision: 1.7 $
 *  $Date: 2009/04/16 08:30:16 $
 *  (last update by $Author: flucke $)
 */
/*
 *# Parameters:
 *#    saveApeToASCII -- Do we write out an APE text file?
 *#    saveComposites -- Do we write APEs for composite detectors?
 *#    saveLocalNotGlobal -- Do we write the APEs in the local or global coordinates?
 *#    apeASCIISaveFile -- The name of the save-file.
 *#    readApeFromASCII -- Do we read in APEs from a text file?
 *#    readLocalNotGlobal -- Do we read APEs in the local or the global frame?
 *#    readFullLocalMatrix -- Do we read the full local matrix or just the diagonal elements?
 *#                        -- Always write full matrix
 *# Full matrix format: DetID dxx dxy dyy dxz dyz dzz
 *# Diagonal element format: DetID sqrt(dxx) sqrt(dyy) sqrt(dzz)
 *#    setComposites -- Do we set the APEs for composite detectors or just ignore them?
 *#    apeASCIIReadFile -- Input file name.
 *# Also note:
 *#    process.AlignmentProducer.saveApeToDB -- to save as an sqlite file
 *# and associated entries in _cfg.py
 */

#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
#include "Alignment/CommonAlignment/interface/AlignableModifier.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "CondFormats/Alignment/interface/AlignmentErrors.h" 
#include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
#include "CLHEP/Matrix/SymMatrix.h"

#include <fstream>
#include <string>
#include <set>

#include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h"
#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"

#include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
#include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
#include "Alignment/CommonAlignment/interface/Alignable.h"

// includes to make known that they inherit from Alignable:
#include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
#include "Alignment/MuonAlignment/interface/AlignableMuon.h"

#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"

class ApeSettingAlgorithm : public AlignmentAlgorithmBase
{
 public:
  /// Constructor
  ApeSettingAlgorithm(const edm::ParameterSet &cfg);

  /// Destructor
  virtual ~ApeSettingAlgorithm();

  /// Call at beginning of job
  virtual void initialize(const edm::EventSetup &setup, AlignableTracker *tracker,
                          AlignableMuon *muon, AlignmentParameterStore *store);

  /// Call at end of job
  virtual void terminate();

  /// Run the algorithm
  virtual void run(const edm::EventSetup &setup, const EventInfo &eventInfo);

 private:
  edm::ParameterSet         theConfig;
  AlignableNavigator       *theAlignableNavigator;
  AlignableTracker         *theTracker;
  bool                     saveApeToAscii_,readApeFromAscii_,readFullLocalMatrix_;
  bool                     readLocalNotGlobal_,saveLocalNotGlobal_;
  bool                     setComposites_,saveComposites_;
};

//____________________________________________________
//____________________________________________________
//____________________________________________________
//____________________________________________________


// Constructor ----------------------------------------------------------------
//____________________________________________________
ApeSettingAlgorithm::ApeSettingAlgorithm(const edm::ParameterSet &cfg) :
  AlignmentAlgorithmBase(cfg), theConfig(cfg),
  theAlignableNavigator(0)
{
  edm::LogInfo("Alignment") << "@SUB=ApeSettingAlgorithm" << "Start.";
  saveApeToAscii_ = theConfig.getUntrackedParameter<bool>("saveApeToASCII");
  saveComposites_ = theConfig.getUntrackedParameter<bool>("saveComposites");
  saveLocalNotGlobal_ = theConfig.getUntrackedParameter<bool>("saveLocalNotGlobal");
  readApeFromAscii_ = theConfig.getParameter<bool>("readApeFromASCII");
  readLocalNotGlobal_ = theConfig.getParameter<bool>("readLocalNotGlobal");
  readFullLocalMatrix_ = theConfig.getParameter<bool>("readFullLocalMatrix");
  setComposites_ = theConfig.getParameter<bool>("setComposites");
  
}

// Destructor ----------------------------------------------------------------
//____________________________________________________
ApeSettingAlgorithm::~ApeSettingAlgorithm()
{
  delete theAlignableNavigator;
}

// Call at beginning of job ---------------------------------------------------
//____________________________________________________
void ApeSettingAlgorithm::initialize(const edm::EventSetup &setup, 
                                             AlignableTracker *tracker, AlignableMuon *muon,
                                             AlignmentParameterStore *store)
{ theAlignableNavigator = new AlignableNavigator(tracker, muon);
 theTracker = tracker;

 if (readApeFromAscii_)
   { std::ifstream apeReadFile(theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath().c_str()); //requires <fstream>
   if (!apeReadFile.good())
     { edm::LogInfo("Alignment") << "@SUB=initialize" <<"Problem opening APE file: skipping"
                                 << theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath();
     return;
     }
   std::set<int> apeList; //To avoid duplicates
   while (!apeReadFile.eof())
     { int apeId=0; double x11,x21,x22,x31,x32,x33;
     if (!readLocalNotGlobal_ || readFullLocalMatrix_) 
       { apeReadFile>>apeId>>x11>>x21>>x22>>x31>>x32>>x33>>std::ws;}
     else
       { apeReadFile>>apeId>>x11>>x22>>x33>>std::ws;}
     //idr What sanity checks do we need to put here?
     if (apeId != 0) //read appears valid?
       if (apeList.find(apeId) == apeList.end()) //Not previously done
         {  DetId id(apeId);
         AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(id)); //NULL if none
         if (alidet)
           { if ((alidet->components().size()<1) || setComposites_) //the problem with glued dets...
             { GlobalError globErr;
               if (readLocalNotGlobal_)
               { AlgebraicSymMatrix as(3,0); 
               if (readFullLocalMatrix_)
                 { as[0][0]=x11; as[1][0]=x21; as[1][1]=x22;
                   as[2][0]=x31; as[2][1]=x32; as[2][2]=x33;
                 }
               else
                 { as[0][0]=x11*x11; as[1][1]=x22*x22; as[2][2]=x33*x33;} //local cov.
               align::RotationType rt=alidet->globalRotation();
               AlgebraicMatrix am(3,3);
               am[0][0]=rt.xx(); am[0][1]=rt.xy(); am[0][2]=rt.xz();
               am[1][0]=rt.yx(); am[1][1]=rt.yy(); am[1][2]=rt.yz();
               am[2][0]=rt.zx(); am[2][1]=rt.zy(); am[2][2]=rt.zz();
               as=as.similarityT(am); //symmetric matrix
               globErr = GlobalError( as );
               }
             else
               {
                 globErr = GlobalError(x11,x21,x22,x31,x32,x33);
               }
             alidet->setAlignmentPositionError(globErr, false); // do not propagate down!
             apeList.insert(apeId); //Flag it's been set
             }
           else
             { edm::LogInfo("Alignment") << "@SUB=initialize" << "Not Setting APE for Composite DetId "<<apeId;
             }
           }
         }
       else
         { edm::LogInfo("Alignment") << "@SUB=initialize" << "Skipping duplicate APE for DetId "<<apeId;
         }
     }
   apeReadFile.close();
   edm::LogInfo("Alignment") << "@SUB=initialize" << "Set "<<apeList.size()<<" APE values.";
   }
}
 

// Call at end of job ---------------------------------------------------------
//____________________________________________________
void ApeSettingAlgorithm::terminate()
{
  if (saveApeToAscii_)
    { AlignmentErrors* aliErr=theTracker->alignmentErrors();
    int theSize=aliErr->m_alignError.size();
    std::ofstream apeSaveFile(theConfig.getUntrackedParameter<std::string>("apeASCIISaveFile").c_str()); //requires <fstream>
    for (int i=0; i < theSize; ++i)
      { int id= aliErr->m_alignError[i].rawId();
      AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(DetId(id))); //NULL if none
      if (alidet && ((alidet->components().size()<1) || saveComposites_))
        { apeSaveFile<<id;
        CLHEP::HepSymMatrix sm = aliErr->m_alignError[i].matrix();
        if (saveLocalNotGlobal_)
          { align::RotationType rt=alidet->globalRotation();
          AlgebraicMatrix am(3,3);
          am[0][0]=rt.xx(); am[0][1]=rt.xy(); am[0][2]=rt.xz();
          am[1][0]=rt.yx(); am[1][1]=rt.yy(); am[1][2]=rt.yz();
          am[2][0]=rt.zx(); am[2][1]=rt.zy(); am[2][2]=rt.zz();
          sm=sm.similarity(am); //symmetric matrix
          } //transform to local
        for (int j=0; j < 3; ++j)
          for (int k=0; k <= j; ++k)
            apeSaveFile<<"  "<<sm[j][k]; //always write full matrix
        
        apeSaveFile<<std::endl;
        }
      }
    delete aliErr;
    apeSaveFile.close();
    }
  // clean up at end:  // FIXME: should we delete here or in destructor?
  delete theAlignableNavigator;
  theAlignableNavigator = 0;
}

// Run the algorithm on trajectories and tracks -------------------------------
//____________________________________________________
void ApeSettingAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo)
{
  // nothing to do here?
}

// Plugin definition for the algorithm
DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory,
                   ApeSettingAlgorithm, "ApeSettingAlgorithm");


Revision:	1.8
Committed:	Fri Apr 17 12:55:40 2009 UTC (16 years ago) by ireid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	CMSSW_3_3_6_patch4, CMSSW_3_3_6_patch6, CMSSW_3_3_6_patch5, CMSSW_3_3_6_patch3, CMSSW_3_3_3_TSG, CMSSW_3_3_6_patch2, CMSSW_3_3_6_patch1, CMSSW_3_3_6, CMSSW_3_3_5, CMSSW_3_3_4, CMSSW_3_3_3_patch1, CMSSW_3_3_3, CMSSW_3_3_2, CMSSW_3_3_1, CMSSW_3_2_8, CMSSW_3_2_7, CMSSW_3_3_0, CMSSW_3_1_4, CMSSW_3_3_0_pre6, CMSSW_3_3_0_pre5, CMSSW_3_1_3, CMSSW_3_3_0_pre4, CMSSW_3_3_0_pre3, CMSSW_3_3_0_pre2, CMSSW_3_2_6, CMSSW_3_3_0_pre1, CMSSW_3_2_5, CMSSW_3_2_4, CMSSW_3_2_3, V03-01-02, CMSSW_3_2_2_patch2, CMSSW_3_2_2_patch1, CMSSW_3_2_2, CMSSW_3_1_2, CMSSW_3_2_1, CMSSW_3_2_0, CMSSW_3_1_1, CMSSW_3_1_0_patch1, CMSSW_3_1_0, CMSSW_3_1_0_pre11, V03-01-01, CMSSW_3_1_0_pre10, CMSSW_3_1_0_pre9, CMSSW_3_1_0_pre8, CMSSW_3_1_0_pre7, V03-01-00, V03-00-05, V03-00-03
Changes since 1.7:	+13 -17 lines
Log Message:	Insert similarity transforms into code after mods to stop propagation of APEs
#	User	Rev	Content
1	ireid	1.1	/**
2			* \file MillePedeAlignmentAlgorithm.cc
3			*
4			* \author : Gero Flucke/Ivan Reid
5	ireid	1.8	* date : February 2009 * $Revision: 1.7 $
6			* $Date: 2009/04/16 08:30:16 $
7			* (last update by $Author: flucke $)
8	ireid	1.6	*/
9			/*
10			*# Parameters:
11			*# saveApeToASCII -- Do we write out an APE text file?
12			*# saveComposites -- Do we write APEs for composite detectors?
13			*# saveLocalNotGlobal -- Do we write the APEs in the local or global coordinates?
14			*# apeASCIISaveFile -- The name of the save-file.
15			*# readApeFromASCII -- Do we read in APEs from a text file?
16			*# readLocalNotGlobal -- Do we read APEs in the local or the global frame?
17			*# readFullLocalMatrix -- Do we read the full local matrix or just the diagonal elements?
18	ireid	1.8	*# -- Always write full matrix
19	ireid	1.6	*# Full matrix format: DetID dxx dxy dyy dxz dyz dzz
20			*# Diagonal element format: DetID sqrt(dxx) sqrt(dyy) sqrt(dzz)
21			*# setComposites -- Do we set the APEs for composite detectors or just ignore them?
22			*# apeASCIIReadFile -- Input file name.
23			*# Also note:
24			*# process.AlignmentProducer.saveApeToDB -- to save as an sqlite file
25			*# and associated entries in _cfg.py
26	ireid	1.1	*/
27
28			#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmPluginFactory.h"
29			#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentAlgorithmBase.h"
30	ireid	1.2	#include "Alignment/CommonAlignment/interface/AlignableModifier.h"
31	ireid	1.1
32			#include "FWCore/ParameterSet/interface/ParameterSet.h"
33
34			#include "CondFormats/Alignment/interface/AlignmentErrors.h"
35	ireid	1.3	#include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
36	ireid	1.1	#include "CLHEP/Matrix/SymMatrix.h"
37
38			#include <fstream>
39			#include <string>
40			#include <set>
41
42			#include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h"
43			#include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentParameterStore.h"
44
45			#include "Alignment/CommonAlignment/interface/AlignableNavigator.h"
46			#include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
47			#include "Alignment/CommonAlignment/interface/Alignable.h"
48
49			// includes to make known that they inherit from Alignable:
50			#include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
51			#include "Alignment/MuonAlignment/interface/AlignableMuon.h"
52
53	ireid	1.3	#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
54	ireid	1.1
55			class ApeSettingAlgorithm : public AlignmentAlgorithmBase
56			{
57			public:
58			/// Constructor
59			ApeSettingAlgorithm(const edm::ParameterSet &cfg);
60
61			/// Destructor
62			virtual ~ApeSettingAlgorithm();
63
64			/// Call at beginning of job
65			virtual void initialize(const edm::EventSetup &setup, AlignableTracker *tracker,
66			AlignableMuon muon, AlignmentParameterStore store);
67
68			/// Call at end of job
69			virtual void terminate();
70
71	flucke	1.5	/// Run the algorithm
72			virtual void run(const edm::EventSetup &setup, const EventInfo &eventInfo);
73	ireid	1.1
74			private:
75			edm::ParameterSet theConfig;
76			AlignableNavigator *theAlignableNavigator;
77			AlignableTracker *theTracker;
78	ireid	1.6	bool saveApeToAscii_,readApeFromAscii_,readFullLocalMatrix_;
79			bool readLocalNotGlobal_,saveLocalNotGlobal_;
80	ireid	1.4	bool setComposites_,saveComposites_;
81	ireid	1.1	};
82
83			//____________________________________________________
84			//____________________________________________________
85			//____________________________________________________
86			//____________________________________________________
87
88
89			// Constructor ----------------------------------------------------------------
90			//____________________________________________________
91			ApeSettingAlgorithm::ApeSettingAlgorithm(const edm::ParameterSet &cfg) :
92			AlignmentAlgorithmBase(cfg), theConfig(cfg),
93			theAlignableNavigator(0)
94			{
95			edm::LogInfo("Alignment") << "@SUB=ApeSettingAlgorithm" << "Start.";
96			saveApeToAscii_ = theConfig.getUntrackedParameter<bool>("saveApeToASCII");
97	ireid	1.4	saveComposites_ = theConfig.getUntrackedParameter<bool>("saveComposites");
98	ireid	1.6	saveLocalNotGlobal_ = theConfig.getUntrackedParameter<bool>("saveLocalNotGlobal");
99	ireid	1.1	readApeFromAscii_ = theConfig.getParameter<bool>("readApeFromASCII");
100	ireid	1.2	readLocalNotGlobal_ = theConfig.getParameter<bool>("readLocalNotGlobal");
101	ireid	1.6	readFullLocalMatrix_ = theConfig.getParameter<bool>("readFullLocalMatrix");
102	ireid	1.4	setComposites_ = theConfig.getParameter<bool>("setComposites");
103
104	ireid	1.1	}
105
106			// Destructor ----------------------------------------------------------------
107			//____________________________________________________
108			ApeSettingAlgorithm::~ApeSettingAlgorithm()
109			{
110			delete theAlignableNavigator;
111			}
112
113			// Call at beginning of job ---------------------------------------------------
114			//____________________________________________________
115			void ApeSettingAlgorithm::initialize(const edm::EventSetup &setup,
116			AlignableTracker tracker, AlignableMuon muon,
117			AlignmentParameterStore *store)
118	ireid	1.2	{ theAlignableNavigator = new AlignableNavigator(tracker, muon);
119			theTracker = tracker;
120	ireid	1.8
121	ireid	1.2	if (readApeFromAscii_)
122			{ std::ifstream apeReadFile(theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath().c_str()); //requires <fstream>
123			if (!apeReadFile.good())
124	ireid	1.3	{ edm::LogInfo("Alignment") << "@SUB=initialize" <<"Problem opening APE file: skipping"
125	ireid	1.2	<< theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath();
126			return;
127			}
128			std::set<int> apeList; //To avoid duplicates
129			while (!apeReadFile.eof())
130			{ int apeId=0; double x11,x21,x22,x31,x32,x33;
131	ireid	1.6	if (!readLocalNotGlobal_ \|\| readFullLocalMatrix_)
132			{ apeReadFile>>apeId>>x11>>x21>>x22>>x31>>x32>>x33>>std::ws;}
133			else
134			{ apeReadFile>>apeId>>x11>>x22>>x33>>std::ws;}
135	ireid	1.2	//idr What sanity checks do we need to put here?
136			if (apeId != 0) //read appears valid?
137			if (apeList.find(apeId) == apeList.end()) //Not previously done
138	ireid	1.3	{ DetId id(apeId);
139	ireid	1.2	AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(id)); //NULL if none
140	ireid	1.4	if (alidet)
141			{ if ((alidet->components().size()<1) \|\| setComposites_) //the problem with glued dets...
142	flucke	1.7	{ GlobalError globErr;
143			if (readLocalNotGlobal_)
144	ireid	1.4	{ AlgebraicSymMatrix as(3,0);
145	ireid	1.6	if (readFullLocalMatrix_)
146			{ as[0][0]=x11; as[1][0]=x21; as[1][1]=x22;
147			as[2][0]=x31; as[2][1]=x32; as[2][2]=x33;
148			}
149			else
150			{ as[0][0]=x11x11; as[1][1]=x22x22; as[2][2]=x33*x33;} //local cov.
151	ireid	1.4	align::RotationType rt=alidet->globalRotation();
152			AlgebraicMatrix am(3,3);
153			am[0][0]=rt.xx(); am[0][1]=rt.xy(); am[0][2]=rt.xz();
154			am[1][0]=rt.yx(); am[1][1]=rt.yy(); am[1][2]=rt.yz();
155			am[2][0]=rt.zx(); am[2][1]=rt.zy(); am[2][2]=rt.zz();
156	ireid	1.8	as=as.similarityT(am); //symmetric matrix
157			globErr = GlobalError( as );
158	ireid	1.4	}
159			else
160	flucke	1.7	{
161			globErr = GlobalError(x11,x21,x22,x31,x32,x33);
162	ireid	1.4	}
163	flucke	1.7	alidet->setAlignmentPositionError(globErr, false); // do not propagate down!
164	ireid	1.4	apeList.insert(apeId); //Flag it's been set
165	ireid	1.3	}
166			else
167	ireid	1.4	{ edm::LogInfo("Alignment") << "@SUB=initialize" << "Not Setting APE for Composite DetId "<<apeId;
168	ireid	1.3	}
169	ireid	1.2	}
170			}
171			else
172			{ edm::LogInfo("Alignment") << "@SUB=initialize" << "Skipping duplicate APE for DetId "<<apeId;
173			}
174			}
175			apeReadFile.close();
176			edm::LogInfo("Alignment") << "@SUB=initialize" << "Set "<<apeList.size()<<" APE values.";
177			}
178			}
179	ireid	1.4
180	ireid	1.1
181			// Call at end of job ---------------------------------------------------------
182			//____________________________________________________
183			void ApeSettingAlgorithm::terminate()
184			{
185			if (saveApeToAscii_)
186			{ AlignmentErrors* aliErr=theTracker->alignmentErrors();
187			int theSize=aliErr->m_alignError.size();
188			std::ofstream apeSaveFile(theConfig.getUntrackedParameter<std::string>("apeASCIISaveFile").c_str()); //requires <fstream>
189			for (int i=0; i < theSize; ++i)
190	ireid	1.4	{ int id= aliErr->m_alignError[i].rawId();
191			AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(DetId(id))); //NULL if none
192			if (alidet && ((alidet->components().size()<1) \|\| saveComposites_))
193			{ apeSaveFile<<id;
194	ireid	1.8	CLHEP::HepSymMatrix sm = aliErr->m_alignError[i].matrix();
195	ireid	1.6	if (saveLocalNotGlobal_)
196			{ align::RotationType rt=alidet->globalRotation();
197			AlgebraicMatrix am(3,3);
198			am[0][0]=rt.xx(); am[0][1]=rt.xy(); am[0][2]=rt.xz();
199			am[1][0]=rt.yx(); am[1][1]=rt.yy(); am[1][2]=rt.yz();
200			am[2][0]=rt.zx(); am[2][1]=rt.zy(); am[2][2]=rt.zz();
201	ireid	1.8	sm=sm.similarity(am); //symmetric matrix
202	ireid	1.6	} //transform to local
203	ireid	1.8	for (int j=0; j < 3; ++j)
204			for (int k=0; k <= j; ++k)
205			apeSaveFile<<" "<<sm[j][k]; //always write full matrix
206
207	ireid	1.4	apeSaveFile<<std::endl;
208			}
209	ireid	1.1	}
210			delete aliErr;
211			apeSaveFile.close();
212			}
213			// clean up at end: // FIXME: should we delete here or in destructor?
214			delete theAlignableNavigator;
215			theAlignableNavigator = 0;
216			}
217
218			// Run the algorithm on trajectories and tracks -------------------------------
219			//____________________________________________________
220	flucke	1.5	void ApeSettingAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo)
221	ireid	1.1	{
222			// nothing to do here?
223			}
224
225			// Plugin definition for the algorithm
226			DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory,
227			ApeSettingAlgorithm, "ApeSettingAlgorithm");
228
229
230