CommonAlignmentAlgorithm/plugins/ApeSettingAlgorithm.cc

 /**
 * \file MillePedeAlignmentAlgorithm.cc
 *
 *  \author    : Gero Flucke/Ivan Reid
 *  date       : February 2009 *  $Revision: 1.11 $
 *  $Date: 2011/03/22 09:49:50 $
 *  (last update by $Author: innocent $)
 */
/*
 *# 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 "Alignment/CommonAlignment/interface/AlignableExtras.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, AlignableExtras *extras,
                          AlignmentParameterStore *store);

  /// Call at end of job
  virtual void terminate(const edm::EventSetup& iSetup);

  /// 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, AlignableExtras *extras,
                                     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_)
               { AlgebraicSymMatrix33 as; 
               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();
               AlgebraicMatrix33 am;
               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();
               globErr = GlobalError(ROOT::Math::SimilarityT(am,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(const edm::EventSetup& iSetup)
{
  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.12
Committed:	Mon Jan 7 20:56:25 2013 UTC (12 years, 3 months ago) by wmtan
Content type:	text/plain
Branch:	MAIN
CVS Tags:	CMSSW_6_2_0, CMSSW_6_2_0_pre7_TS133806, CMSSW_6_2_0_pre8, V04-00-14, V04-00-13, CMSSW_6_2_0_pre7_TS132947, CMSSW_6_2_0_pre7_g496p02, CMSSW_6_2_0_pre7, V04-00-12, CMSSW_6_2_0_pre6_patch1, CMSSW_6_2_0_pre6, CMSSW_6_2_0_pre5slc6, CMSSW_6_2_0_pre5, CMSSW_6_2_0_pre4, CMSSW_6_2_0_pre3, CMSSW_6_2_0_pre2, CMSSW_6_2_0_pre1, V04-00-11, HEAD
Changes since 1.11:	+5 -5 lines
Log Message:	Get geometry from the EventSetup system, rather than having it hard coded in the data formats
#	User	Rev	Content
1	ireid	1.1	/**
2			* \file MillePedeAlignmentAlgorithm.cc
3			*
4			* \author : Gero Flucke/Ivan Reid
5	wmtan	1.12	* date : February 2009 * $Revision: 1.11 $
6			* $Date: 2011/03/22 09:49:50 $
7			* (last update by $Author: innocent $)
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	mussgill	1.10	#include "Alignment/CommonAlignment/interface/AlignableExtras.h"
53	ireid	1.1
54	ireid	1.3	#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
55	ireid	1.1
56			class ApeSettingAlgorithm : public AlignmentAlgorithmBase
57			{
58			public:
59			/// Constructor
60			ApeSettingAlgorithm(const edm::ParameterSet &cfg);
61
62			/// Destructor
63			virtual ~ApeSettingAlgorithm();
64
65			/// Call at beginning of job
66	mussgill	1.10	virtual void initialize(const edm::EventSetup &setup,
67			AlignableTracker tracker, AlignableMuon muon, AlignableExtras *extras,
68			AlignmentParameterStore *store);
69	ireid	1.1
70			/// Call at end of job
71	wmtan	1.12	virtual void terminate(const edm::EventSetup& iSetup);
72	ireid	1.1
73	flucke	1.5	/// Run the algorithm
74			virtual void run(const edm::EventSetup &setup, const EventInfo &eventInfo);
75	ireid	1.1
76			private:
77			edm::ParameterSet theConfig;
78			AlignableNavigator *theAlignableNavigator;
79			AlignableTracker *theTracker;
80	ireid	1.6	bool saveApeToAscii_,readApeFromAscii_,readFullLocalMatrix_;
81			bool readLocalNotGlobal_,saveLocalNotGlobal_;
82	ireid	1.4	bool setComposites_,saveComposites_;
83	ireid	1.1	};
84
85			//____________________________________________________
86			//____________________________________________________
87			//____________________________________________________
88			//____________________________________________________
89
90
91			// Constructor ----------------------------------------------------------------
92			//____________________________________________________
93			ApeSettingAlgorithm::ApeSettingAlgorithm(const edm::ParameterSet &cfg) :
94			AlignmentAlgorithmBase(cfg), theConfig(cfg),
95			theAlignableNavigator(0)
96			{
97			edm::LogInfo("Alignment") << "@SUB=ApeSettingAlgorithm" << "Start.";
98			saveApeToAscii_ = theConfig.getUntrackedParameter<bool>("saveApeToASCII");
99	ireid	1.4	saveComposites_ = theConfig.getUntrackedParameter<bool>("saveComposites");
100	ireid	1.6	saveLocalNotGlobal_ = theConfig.getUntrackedParameter<bool>("saveLocalNotGlobal");
101	ireid	1.1	readApeFromAscii_ = theConfig.getParameter<bool>("readApeFromASCII");
102	ireid	1.2	readLocalNotGlobal_ = theConfig.getParameter<bool>("readLocalNotGlobal");
103	ireid	1.6	readFullLocalMatrix_ = theConfig.getParameter<bool>("readFullLocalMatrix");
104	ireid	1.4	setComposites_ = theConfig.getParameter<bool>("setComposites");
105
106	ireid	1.1	}
107
108			// Destructor ----------------------------------------------------------------
109			//____________________________________________________
110			ApeSettingAlgorithm::~ApeSettingAlgorithm()
111			{
112			delete theAlignableNavigator;
113			}
114
115			// Call at beginning of job ---------------------------------------------------
116			//____________________________________________________
117			void ApeSettingAlgorithm::initialize(const edm::EventSetup &setup,
118	mussgill	1.10	AlignableTracker tracker, AlignableMuon muon, AlignableExtras *extras,
119	ireid	1.9	AlignmentParameterStore *store)
120	mussgill	1.10	{
121			theAlignableNavigator = new AlignableNavigator(tracker, muon);
122	ireid	1.2	theTracker = tracker;
123	ireid	1.9
124	ireid	1.2	if (readApeFromAscii_)
125			{ std::ifstream apeReadFile(theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath().c_str()); //requires <fstream>
126			if (!apeReadFile.good())
127	ireid	1.3	{ edm::LogInfo("Alignment") << "@SUB=initialize" <<"Problem opening APE file: skipping"
128	ireid	1.2	<< theConfig.getParameter<edm::FileInPath>("apeASCIIReadFile").fullPath();
129			return;
130			}
131			std::set<int> apeList; //To avoid duplicates
132			while (!apeReadFile.eof())
133			{ int apeId=0; double x11,x21,x22,x31,x32,x33;
134	ireid	1.6	if (!readLocalNotGlobal_ \|\| readFullLocalMatrix_)
135			{ apeReadFile>>apeId>>x11>>x21>>x22>>x31>>x32>>x33>>std::ws;}
136			else
137			{ apeReadFile>>apeId>>x11>>x22>>x33>>std::ws;}
138	ireid	1.2	//idr What sanity checks do we need to put here?
139			if (apeId != 0) //read appears valid?
140	ireid	1.9	{ if (apeList.find(apeId) == apeList.end()) //Not previously done
141	ireid	1.3	{ DetId id(apeId);
142	ireid	1.2	AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(id)); //NULL if none
143	ireid	1.4	if (alidet)
144			{ if ((alidet->components().size()<1) \|\| setComposites_) //the problem with glued dets...
145	flucke	1.7	{ GlobalError globErr;
146	ireid	1.9	if (readLocalNotGlobal_)
147	innocent	1.11	{ AlgebraicSymMatrix33 as;
148	ireid	1.6	if (readFullLocalMatrix_)
149			{ as[0][0]=x11; as[1][0]=x21; as[1][1]=x22;
150	ireid	1.9	as[2][0]=x31; as[2][1]=x32; as[2][2]=x33;
151	ireid	1.6	}
152			else
153			{ as[0][0]=x11x11; as[1][1]=x22x22; as[2][2]=x33*x33;} //local cov.
154	ireid	1.4	align::RotationType rt=alidet->globalRotation();
155	innocent	1.11	AlgebraicMatrix33 am;
156	ireid	1.4	am[0][0]=rt.xx(); am[0][1]=rt.xy(); am[0][2]=rt.xz();
157			am[1][0]=rt.yx(); am[1][1]=rt.yy(); am[1][2]=rt.yz();
158			am[2][0]=rt.zx(); am[2][1]=rt.zy(); am[2][2]=rt.zz();
159	innocent	1.11	globErr = GlobalError(ROOT::Math::SimilarityT(am,as));
160	ireid	1.4	}
161			else
162	flucke	1.7	{
163			globErr = GlobalError(x11,x21,x22,x31,x32,x33);
164	ireid	1.4	}
165	flucke	1.7	alidet->setAlignmentPositionError(globErr, false); // do not propagate down!
166	ireid	1.4	apeList.insert(apeId); //Flag it's been set
167	ireid	1.3	}
168			else
169	ireid	1.4	{ edm::LogInfo("Alignment") << "@SUB=initialize" << "Not Setting APE for Composite DetId "<<apeId;
170	ireid	1.3	}
171	ireid	1.2	}
172			}
173			else
174			{ edm::LogInfo("Alignment") << "@SUB=initialize" << "Skipping duplicate APE for DetId "<<apeId;
175			}
176	ireid	1.9	}
177	ireid	1.2	}
178			apeReadFile.close();
179			edm::LogInfo("Alignment") << "@SUB=initialize" << "Set "<<apeList.size()<<" APE values.";
180			}
181			}
182	ireid	1.4
183	ireid	1.1
184			// Call at end of job ---------------------------------------------------------
185			//____________________________________________________
186	wmtan	1.12	void ApeSettingAlgorithm::terminate(const edm::EventSetup& iSetup)
187	ireid	1.1	{
188			if (saveApeToAscii_)
189			{ AlignmentErrors* aliErr=theTracker->alignmentErrors();
190			int theSize=aliErr->m_alignError.size();
191			std::ofstream apeSaveFile(theConfig.getUntrackedParameter<std::string>("apeASCIISaveFile").c_str()); //requires <fstream>
192			for (int i=0; i < theSize; ++i)
193	ireid	1.4	{ int id= aliErr->m_alignError[i].rawId();
194			AlignableDetOrUnitPtr alidet(theAlignableNavigator->alignableFromDetId(DetId(id))); //NULL if none
195			if (alidet && ((alidet->components().size()<1) \|\| saveComposites_))
196			{ apeSaveFile<<id;
197	ireid	1.8	CLHEP::HepSymMatrix sm = aliErr->m_alignError[i].matrix();
198	ireid	1.6	if (saveLocalNotGlobal_)
199			{ align::RotationType rt=alidet->globalRotation();
200			AlgebraicMatrix am(3,3);
201			am[0][0]=rt.xx(); am[0][1]=rt.xy(); am[0][2]=rt.xz();
202			am[1][0]=rt.yx(); am[1][1]=rt.yy(); am[1][2]=rt.yz();
203			am[2][0]=rt.zx(); am[2][1]=rt.zy(); am[2][2]=rt.zz();
204	ireid	1.8	sm=sm.similarity(am); //symmetric matrix
205	ireid	1.6	} //transform to local
206	ireid	1.8	for (int j=0; j < 3; ++j)
207			for (int k=0; k <= j; ++k)
208			apeSaveFile<<" "<<sm[j][k]; //always write full matrix
209
210	ireid	1.4	apeSaveFile<<std::endl;
211			}
212	ireid	1.1	}
213			delete aliErr;
214			apeSaveFile.close();
215			}
216			// clean up at end: // FIXME: should we delete here or in destructor?
217			delete theAlignableNavigator;
218			theAlignableNavigator = 0;
219			}
220
221			// Run the algorithm on trajectories and tracks -------------------------------
222			//____________________________________________________
223	flucke	1.5	void ApeSettingAlgorithm::run(const edm::EventSetup &setup, const EventInfo &eventInfo)
224	ireid	1.1	{
225			// nothing to do here?
226			}
227
228			// Plugin definition for the algorithm
229			DEFINE_EDM_PLUGIN(AlignmentAlgorithmPluginFactory,
230			ApeSettingAlgorithm, "ApeSettingAlgorithm");
231
232
233