CommonAlignmentAlgorithm/plugins/SiPixelLorentzAngleCalibration.cc

/// \class SiPixelLorentzAngleCalibration
///
/// Calibration of Lorentz angle for the pixel,
/// integrated in the alignment algorithms.
///
/// Note that not all algorithms support this...
///
///  \author    : Gero Flucke
///  date       : September 2012
///  $Revision: 1.1 $
///  $Date: 2012/09/05 07:56:27 $
///  (last update by $Author: flucke $)

#include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"

#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
#include "CondFormats/DataRecord/interface/SiPixelLorentzAngleRcd.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"

#include "TTree.h"
#include "TFile.h"
#include "TString.h"

// #include <iostream>
#include <vector>
#include <sstream>

class SiPixelLorentzAngleCalibration : public IntegratedCalibrationBase
{
public:
  /// Constructor
  explicit SiPixelLorentzAngleCalibration(const edm::ParameterSet &cfg);
  
  /// Destructor
  virtual ~SiPixelLorentzAngleCalibration();

  /// How many parameters does this calibration define?
  virtual unsigned int numParameters() const;

  // /// Return all derivatives,
  // /// default implementation uses other derivatives(..) method,
  // /// but can be overwritten in derived class for efficiency.
  // virtual std::vector<double> derivatives(const TransientTrackingRecHit &hit,
  //                                      const TrajectoryStateOnSurface &tsos,
  //                                      const edm::EventSetup &setup,
  //                                      const EventInfo &eventInfo) const;

  /// Return non-zero derivatives for x- and y-measurements with their indices by reference.
  /// Return value is their number.
  virtual unsigned int derivatives(std::vector<ValuesIndexPair> &outDerivInds,
                                   const TransientTrackingRecHit &hit,
                                   const TrajectoryStateOnSurface &tsos,
                                   const edm::EventSetup &setup,
                                   const EventInfo &eventInfo) const;

  /// Setting the determined parameter identified by index,
  /// should return false if out-of-bounds, true otherwise.
  virtual bool setParameter(unsigned int index, double value);

  /// Setting the determined parameter uncertainty identified by index,
  /// returns false if out-of-bounds, true otherwise.
  virtual bool setParameterError(unsigned int index, double error);

  /// Return current value of parameter identified by index.
  /// Should return 0. if index out-of-bounds.
  virtual double getParameter(unsigned int index) const;

  /// Return current value of parameter identified by index.
  /// Returns 0. if index out-of-bounds or if errors undetermined.
  virtual double getParameterError(unsigned int index) const;

  // /// Call at beginning of job:
  // virtual void beginOfJob(const AlignableTracker *tracker,
  //                      const AlignableMuon *muon,
  //                      const AlignableExtras *extras);

  /// Called at end of a the job of the AlignmentProducer.
  /// Write out determined parameters.
  virtual void endOfJob();

private:
  /// If called the first time, fill 'siPixelLorentzAngleInput_',
  /// later check that LorentzAngle has not changed.
  bool checkLorentzAngleInput(const edm::EventSetup &setup, const EventInfo &eventInfo);
  /// Input LorentzAngle values:
  /// - either from EventSetup of first call to derivatives(..)
  /// - or created from files of passed by configuration (i.e. from parallel processing)
  const SiPixelLorentzAngle* getLorentzAnglesInput();

  /// Determined parameter value for this detId (detId not treated => 0.)
  /// and the given run.
  double getParameterForDetId(unsigned int detId, edm::RunNumber_t run) const;
  /// Index of parameter for given detId (detId not treated => < 0)
  /// and the given run.
  int getParameterIndexFromDetId(unsigned int detId, edm::RunNumber_t run) const;
  /// Total number of IOVs.
  unsigned int numIovs() const;
  /// First run of iov (0 if iovNum not treated).
  edm::RunNumber_t firstRunOfIOV(unsigned int iovNum) const;

  void writeTree(const SiPixelLorentzAngle *lorentzAngle, const char *treeName) const;
  SiPixelLorentzAngle* createFromTree(const char *fileName, const char *treeName) const;

  const bool saveToDB_;
  const std::string recordNameDBwrite_;
  const std::string outFileName_;
  const std::vector<std::string> mergeFileNames_;

  edm::ESWatcher<SiPixelLorentzAngleRcd> watchLorentzAngleRcd_;

  // const AlignableTracker *alignableTracker_;
  SiPixelLorentzAngle *siPixelLorentzAngleInput_;
  std::vector<double> parameters_;
  std::vector<double> paramUncertainties_;
};

//======================================================================
//======================================================================
//======================================================================

SiPixelLorentzAngleCalibration::SiPixelLorentzAngleCalibration(const edm::ParameterSet &cfg)
  : IntegratedCalibrationBase(cfg),
    saveToDB_(cfg.getParameter<bool>("saveToDB")),
    recordNameDBwrite_(cfg.getParameter<std::string>("recordNameDBwrite")),
    outFileName_(cfg.getParameter<std::string>("treeFile")),
    mergeFileNames_(cfg.getParameter<std::vector<std::string> >("mergeTreeFiles")),
    //    alignableTracker_(0),
    siPixelLorentzAngleInput_(0)
{
  // FIXME: Which granularity, leading to how many parameters?
  parameters_.resize(2, 0.); // currently two parameters (ring1-4, 5-8), start value 0.
  paramUncertainties_.resize(2, 0.); // dito for errors

  edm::LogInfo("Alignment") << "@SUB=SiPixelLorentzAngleCalibration" << "Created with name "
                            << this->name() << "',\n" << this->numParameters() << " parameters to be determined,"
                            << "\nsaveToDB = " << saveToDB_
                            << "\n outFileName = " << outFileName_
                            << "\n N(merge files) = " << mergeFileNames_.size();
  if (mergeFileNames_.size()) {
    edm::LogInfo("Alignment") << "@SUB=SiPixelLorentzAngleCalibration"
                              << "First file to merge: " << mergeFileNames_[0];
  }
}
  
//======================================================================
SiPixelLorentzAngleCalibration::~SiPixelLorentzAngleCalibration()
{
  //  std::cout << "Destroy SiPixelLorentzAngleCalibration named " << this->name() << std::endl;
  delete siPixelLorentzAngleInput_;
}

//======================================================================
unsigned int SiPixelLorentzAngleCalibration::numParameters() const
{
  return parameters_.size();
}

//======================================================================
unsigned int
SiPixelLorentzAngleCalibration::derivatives(std::vector<ValuesIndexPair> &outDerivInds,
                                            const TransientTrackingRecHit &hit,
                                            const TrajectoryStateOnSurface &tsos,
                                            const edm::EventSetup &setup,
                                            const EventInfo &eventInfo) const
{
  // ugly const-cast:
  // But it is either only first initialisation or throwing an exception...
  const_cast<SiPixelLorentzAngleCalibration*>(this)->checkLorentzAngleInput(setup, eventInfo);

  outDerivInds.clear();

  if (hit.det()) { // otherwise 'constraint hit' or whatever
    
    const int index = this->getParameterIndexFromDetId(hit.det()->geographicalId(),
                                                       eventInfo.eventId_.run());
    if (index >= 0) { // otherwise not treated
      edm::ESHandle<MagneticField> magneticField;
      setup.get<IdealMagneticFieldRecord>().get(magneticField);
      const GlobalVector bField(magneticField->inTesla(hit.det()->surface().position()));
      const LocalVector bFieldLocal(hit.det()->surface().toLocal(bField));
      const double dZ = hit.det()->surface().bounds().thickness(); // it is a float only...
      // shift due to LA: dx = tan(LA) * dz = mobility * B_y * dz
      const double xDerivative = bFieldLocal.y() * dZ; // parameter is mobility!
      if (xDerivative) { // If field is zero, this is zero: do not return it
        const Values derivs(xDerivative, 0.); // yDerivative = 0.
        outDerivInds.push_back(ValuesIndexPair(derivs, index));
      }
    }
  } else {
    edm::LogWarning("Alignment") << "@SUB=SiPixelLorentzAngleCalibration::derivatives2"
                                 << "Hit without GeomDet, skip!";
  }
  
  return outDerivInds.size();
}

//======================================================================
bool SiPixelLorentzAngleCalibration::setParameter(unsigned int index, double value)
{
  if (index >= parameters_.size()) {
    return false;
  } else {
    parameters_[index] = value;
    return true;
  }
}

//======================================================================
bool SiPixelLorentzAngleCalibration::setParameterError(unsigned int index, double error)
{
  if (index >= paramUncertainties_.size()) {
    return false;
  } else {
    paramUncertainties_[index] = error;
    return true;
  }
}

//======================================================================
double SiPixelLorentzAngleCalibration::getParameter(unsigned int index) const
{
  //   if (index >= parameters_.size()) {
  //     return 0.;
  //   } else {
  //     return parameters_[index];
  //   }
  return (index >= parameters_.size() ? 0. : parameters_[index]);
}

//======================================================================
double SiPixelLorentzAngleCalibration::getParameterError(unsigned int index) const
{
  //   if (index >= paramUncertainties_.size()) {
  //     return 0.;
  //   } else {
  //     return paramUncertainties_[index];
  //   }
  return (index >= paramUncertainties_.size() ? 0. : paramUncertainties_[index]);
}

// //======================================================================
// void SiPixelLorentzAngleCalibration::beginOfJob(const AlignableTracker *tracker,
//                                                 const AlignableMuon */*muon*/,
//                                                 const AlignableExtras */*extras*/)
// {
//   alignableTracker_ = tracker;
// }


//======================================================================
void SiPixelLorentzAngleCalibration::endOfJob()
{
  // loginfo output
  std::ostringstream out;
  out << "Parameter results\n";
  for (unsigned int iPar = 0; iPar < parameters_.size(); ++iPar) {
    out << iPar << ": " << parameters_[iPar] << " +- " << paramUncertainties_[iPar] << "\n";
  }
  edm::LogInfo("Alignment") << "@SUB=SiPixelLorentzAngleCalibration::endOfJob" << out.str();

  // now write 'input' tree
  const SiPixelLorentzAngle *input = this->getLorentzAnglesInput(); // never NULL
  const std::string treeName(this->name() + '_');
  this->writeTree(input, (treeName + "input").c_str());

  for (unsigned int iIOV = 0; iIOV < this->numIovs(); ++iIOV) {
    cond::Time_t firstRunOfIOV = this->firstRunOfIOV(iIOV);
    SiPixelLorentzAngle *output = new SiPixelLorentzAngle;
    // Loop on map of values from input and add (possible) parameter results
    for (auto iterIdValue = input->getLorentzAngles().begin();
         iterIdValue != input->getLorentzAngles().end(); ++iterIdValue) {
      // type of (*iterIdValue) is pair<unsigned int, float>
      const unsigned int detId = iterIdValue->first; // key of map is DetId
      // Nasty: putLorentzAngle(..) takes float by reference - not even const reference!
      float value = iterIdValue->second + this->getParameterForDetId(detId, firstRunOfIOV);
      output->putLorentzAngle(detId, value); // put result in output
    }

    // Write this even for mille jobs?
    this->writeTree(output, (treeName + Form("result_%lld", firstRunOfIOV)).c_str());

    if (saveToDB_) { // If requested, write out to DB 
      edm::Service<cond::service::PoolDBOutputService> dbService;
      if (dbService.isAvailable()) {
        dbService->writeOne(output, firstRunOfIOV, recordNameDBwrite_.c_str());
        // no 'delete output;': writeOne(..) took over ownership
      } else {
        delete output;
        edm::LogError("BadConfig") << "@SUB=SiPixelLorentzAngleCalibration::endOfJob"
                                   << "No PoolDBOutputService available, but saveToDB true!";
      }
    } else {
      delete output;
    }
  } // end loop on IOVs

}

//======================================================================
bool SiPixelLorentzAngleCalibration::checkLorentzAngleInput(const edm::EventSetup &setup,
                                                            const EventInfo &eventInfo)
{
  edm::ESHandle<SiPixelLorentzAngle> lorentzAngleHandle;
  if (!siPixelLorentzAngleInput_) {
    setup.get<SiPixelLorentzAngleRcd>().get(lorentzAngleHandle);
    siPixelLorentzAngleInput_ = new SiPixelLorentzAngle(*lorentzAngleHandle);
  } else {
    if (watchLorentzAngleRcd_.check(setup)) { // new IOV of input
      setup.get<SiPixelLorentzAngleRcd>().get(lorentzAngleHandle);
      if (lorentzAngleHandle->getLorentzAngles() // but only bad if non-identical values
          != siPixelLorentzAngleInput_->getLorentzAngles()) { // (comparing maps)
        // FIXME: Could different maps have same content, but different order?
        //        Or 'floating point comparison' problem?
        throw cms::Exception("BadInput")
          << "SiPixelLorentzAngleCalibration::checkLorentzAngleInput:\n"
          << "Content of SiPixelLorentzAngle changed at run " << eventInfo.eventId_.run()
          << ", but algorithm expects constant input!\n";
        return false; // not reached...
      }
    }
  }
  
  return true;
}

//======================================================================
const SiPixelLorentzAngle* SiPixelLorentzAngleCalibration::getLorentzAnglesInput()
{
  // For parallel processing in Millepede II, create SiPixelLorentzAngle
  // from info stored in files of parallel jobs and check that they are identical.
  // If this job has run on data, still check that LA is identical to the ones
  // from mergeFileNames_.
  const std::string treeName(this->name() + "_input");
  for (auto iFile = mergeFileNames_.begin(); iFile != mergeFileNames_.end(); ++iFile) {
    SiPixelLorentzAngle* la = this->createFromTree(iFile->c_str(), treeName.c_str());
    // siPixelLorentzAngleInput_ could be non-null from previous file of this loop
    // or from checkLorentzAngleInput(..) when running on data in this job as well
    if (!siPixelLorentzAngleInput_ || siPixelLorentzAngleInput_->getLorentzAngles().empty()) {
      delete siPixelLorentzAngleInput_; // NULL or empty
      siPixelLorentzAngleInput_ = la;
    } else {
      // FIXME: about comparison of maps see comments in checkLorentzAngleInput
      if (!la->getLorentzAngles().empty() && // single job might not have got events
          la->getLorentzAngles() != siPixelLorentzAngleInput_->getLorentzAngles()) {
        // Throw exception instead of error?
        edm::LogError("NoInput") << "@SUB=SiPixelLorentzAngleCalibration::getLorentzAnglesInput"
                                 << "Different input values from tree " << treeName
                                 << " in file " << *iFile << ".";
        
      }
      delete la;
    }
  }

  if (!siPixelLorentzAngleInput_) { // no files nor ran on events
    siPixelLorentzAngleInput_ = new SiPixelLorentzAngle;
    edm::LogError("NoInput") << "@SUB=SiPixelLorentzAngleCalibration::getLorentzAnglesInput"
                             << "No input, create an empty one!";
  } else if (siPixelLorentzAngleInput_->getLorentzAngles().empty()) {
    edm::LogError("NoInput") << "@SUB=SiPixelLorentzAngleCalibration::getLorentzAnglesInput"
                             << "Empty result!";
  }

  return siPixelLorentzAngleInput_;
}

//======================================================================
double SiPixelLorentzAngleCalibration::getParameterForDetId(unsigned int detId,
                                                            edm::RunNumber_t run) const
{
  const int index = this->getParameterIndexFromDetId(detId, run);

  return (index < 0 ? 0. : parameters_[index]);
}

//======================================================================
int SiPixelLorentzAngleCalibration::getParameterIndexFromDetId(unsigned int detId,
                                                               edm::RunNumber_t run) const
{
  // Return the index of the parameter that is used for this DetId.
  // If this DetId is not treated, return values < 0.
  
  // FIXME: Extend to configurable granularity? 
  //        Including treatment of run dependence?
  const PXBDetId id(detId);
  if (id.det() == DetId::Tracker && id.subdetId() == PixelSubdetector::PixelBarrel) {
    if (id.module() >= 1 && id.module() <= 4) return 0;
    if (id.module() >= 5 && id.module() <= 8) return 1;
    edm::LogWarning("Alignment")
        << "@SUB=SiPixelLorentzAngleCalibration::getParameterIndexFromDetId"
        << "Module should be 1-8, but is " << id.module() << " => skip!";
  }

  return -1;
}

//======================================================================
unsigned int SiPixelLorentzAngleCalibration::numIovs() const
{
  // FIXME: Needed to include treatment of run dependence!
  return 1; 
}

//======================================================================
edm::RunNumber_t SiPixelLorentzAngleCalibration::firstRunOfIOV(unsigned int iovNum) const
{
  // FIXME: Needed to include treatment of run dependence!
  if (iovNum < this->numIovs()) return 1;
  else return 0;
}


//======================================================================
void SiPixelLorentzAngleCalibration::writeTree(const SiPixelLorentzAngle *lorentzAngle,
                                               const char *treeName) const
{
  if (!lorentzAngle) return;

  TFile* file = TFile::Open(outFileName_.c_str(), "UPDATE");
  if (!file) {
    edm::LogError("BadConfig") << "@SUB=SiPixelLorentzAngleCalibration::writeTree"
                               << "Could not open file '" << outFileName_ << "'.";
    return;
  }

  TTree *tree = new TTree(treeName, treeName);
  unsigned int id = 0;
  float value = 0.;
  tree->Branch("detId", &id, "detId/i");
  tree->Branch("value", &value, "value/F");

  for (auto iterIdValue = lorentzAngle->getLorentzAngles().begin();
       iterIdValue != lorentzAngle->getLorentzAngles().end(); ++iterIdValue) {
    // type of (*iterIdValue) is pair<unsigned int, float>
    id = iterIdValue->first; // key of map is DetId
    value = iterIdValue->second;
    tree->Fill();
  }
  tree->Write();
  delete file; // tree vanishes with the file... (?)

}

//======================================================================
SiPixelLorentzAngle* 
SiPixelLorentzAngleCalibration::createFromTree(const char *fileName, const char *treeName) const
{
  TTree *tree = 0;
  TFile* file = TFile::Open(fileName, "READ");
  if (file) file->GetObject(treeName, tree);

  SiPixelLorentzAngle *result = 0;
  if (tree) {
    unsigned int id = 0;
    float value = 0.;
    tree->SetBranchAddress("detId", &id);
    tree->SetBranchAddress("value", &value);

    result = new SiPixelLorentzAngle;
    const Long64_t nEntries = tree->GetEntries();
    for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
      tree->GetEntry(iEntry);
      result->putLorentzAngle(id, value);
    }
  } else {
    edm::LogError("Alignment") << "@SUB=SiPixelLorentzAngleCalibration::createFromTree"
                               << "Could not get TTree '" << treeName << "' from file '"
                               << fileName << "'.";
  }

  delete file; // tree will vanish with file
  return result;
}


//======================================================================
//======================================================================
// Plugin definition

#include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationPluginFactory.h"

DEFINE_EDM_PLUGIN(IntegratedCalibrationPluginFactory,
                   SiPixelLorentzAngleCalibration, "SiPixelLorentzAngleCalibration");
Revision:	1.1
Committed:	Wed Sep 5 11:50:22 2012 UTC (12 years, 8 months ago) by flucke
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V04-00-05
Log Message:	add SiPixelLorentzAngleCalibration: basically copy paste of SiStripLorentzAngleCalibration - replacing SiStripLorentzAngle by SiPixelLorentzAngle - removing peak vs deco mode - currently foresee BPIX, separated for ring 1-4 and 5-8
#	User	Rev	Content
1	flucke	1.1	/// \class SiPixelLorentzAngleCalibration
2			///
3			/// Calibration of Lorentz angle for the pixel,
4			/// integrated in the alignment algorithms.
5			///
6			/// Note that not all algorithms support this...
7			///
8			/// \author : Gero Flucke
9			/// date : September 2012
10			/// $Revision: 1.1 $
11			/// $Date: 2012/09/05 07:56:27 $
12			/// (last update by $Author: flucke $)
13
14			#include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationBase.h"
15
16			#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
17			#include "CondFormats/DataRecord/interface/SiPixelLorentzAngleRcd.h"
18			#include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
19
20			#include "DataFormats/DetId/interface/DetId.h"
21			#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
22
23			#include "FWCore/Framework/interface/ESHandle.h"
24			#include "FWCore/Framework/interface/ESWatcher.h"
25			#include "FWCore/MessageLogger/interface/MessageLogger.h"
26			#include "FWCore/ServiceRegistry/interface/Service.h"
27
28			#include "Geometry/CommonDetUnit/interface/GeomDet.h"
29			#include "MagneticField/Engine/interface/MagneticField.h"
30			#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
31			#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
32
33			#include "TTree.h"
34			#include "TFile.h"
35			#include "TString.h"
36
37			// #include <iostream>
38			#include <vector>
39			#include <sstream>
40
41			class SiPixelLorentzAngleCalibration : public IntegratedCalibrationBase
42			{
43			public:
44			/// Constructor
45			explicit SiPixelLorentzAngleCalibration(const edm::ParameterSet &cfg);
46
47			/// Destructor
48			virtual ~SiPixelLorentzAngleCalibration();
49
50			/// How many parameters does this calibration define?
51			virtual unsigned int numParameters() const;
52
53			// /// Return all derivatives,
54			// /// default implementation uses other derivatives(..) method,
55			// /// but can be overwritten in derived class for efficiency.
56			// virtual std::vector<double> derivatives(const TransientTrackingRecHit &hit,
57			// const TrajectoryStateOnSurface &tsos,
58			// const edm::EventSetup &setup,
59			// const EventInfo &eventInfo) const;
60
61			/// Return non-zero derivatives for x- and y-measurements with their indices by reference.
62			/// Return value is their number.
63			virtual unsigned int derivatives(std::vector<ValuesIndexPair> &outDerivInds,
64			const TransientTrackingRecHit &hit,
65			const TrajectoryStateOnSurface &tsos,
66			const edm::EventSetup &setup,
67			const EventInfo &eventInfo) const;
68
69			/// Setting the determined parameter identified by index,
70			/// should return false if out-of-bounds, true otherwise.
71			virtual bool setParameter(unsigned int index, double value);
72
73			/// Setting the determined parameter uncertainty identified by index,
74			/// returns false if out-of-bounds, true otherwise.
75			virtual bool setParameterError(unsigned int index, double error);
76
77			/// Return current value of parameter identified by index.
78			/// Should return 0. if index out-of-bounds.
79			virtual double getParameter(unsigned int index) const;
80
81			/// Return current value of parameter identified by index.
82			/// Returns 0. if index out-of-bounds or if errors undetermined.
83			virtual double getParameterError(unsigned int index) const;
84
85			// /// Call at beginning of job:
86			// virtual void beginOfJob(const AlignableTracker *tracker,
87			// const AlignableMuon *muon,
88			// const AlignableExtras *extras);
89
90			/// Called at end of a the job of the AlignmentProducer.
91			/// Write out determined parameters.
92			virtual void endOfJob();
93
94			private:
95			/// If called the first time, fill 'siPixelLorentzAngleInput_',
96			/// later check that LorentzAngle has not changed.
97			bool checkLorentzAngleInput(const edm::EventSetup &setup, const EventInfo &eventInfo);
98			/// Input LorentzAngle values:
99			/// - either from EventSetup of first call to derivatives(..)
100			/// - or created from files of passed by configuration (i.e. from parallel processing)
101			const SiPixelLorentzAngle* getLorentzAnglesInput();
102
103			/// Determined parameter value for this detId (detId not treated => 0.)
104			/// and the given run.
105			double getParameterForDetId(unsigned int detId, edm::RunNumber_t run) const;
106			/// Index of parameter for given detId (detId not treated => < 0)
107			/// and the given run.
108			int getParameterIndexFromDetId(unsigned int detId, edm::RunNumber_t run) const;
109			/// Total number of IOVs.
110			unsigned int numIovs() const;
111			/// First run of iov (0 if iovNum not treated).
112			edm::RunNumber_t firstRunOfIOV(unsigned int iovNum) const;
113
114			void writeTree(const SiPixelLorentzAngle lorentzAngle, const char treeName) const;
115			SiPixelLorentzAngle* createFromTree(const char fileName, const char treeName) const;
116
117			const bool saveToDB_;
118			const std::string recordNameDBwrite_;
119			const std::string outFileName_;
120			const std::vector<std::string> mergeFileNames_;
121
122			edm::ESWatcher<SiPixelLorentzAngleRcd> watchLorentzAngleRcd_;
123
124			// const AlignableTracker *alignableTracker_;
125			SiPixelLorentzAngle *siPixelLorentzAngleInput_;
126			std::vector<double> parameters_;
127			std::vector<double> paramUncertainties_;
128			};
129
130			//======================================================================
131			//======================================================================
132			//======================================================================
133
134			SiPixelLorentzAngleCalibration::SiPixelLorentzAngleCalibration(const edm::ParameterSet &cfg)
135			: IntegratedCalibrationBase(cfg),
136			saveToDB_(cfg.getParameter<bool>("saveToDB")),
137			recordNameDBwrite_(cfg.getParameter<std::string>("recordNameDBwrite")),
138			outFileName_(cfg.getParameter<std::string>("treeFile")),
139			mergeFileNames_(cfg.getParameter<std::vector<std::string> >("mergeTreeFiles")),
140			// alignableTracker_(0),
141			siPixelLorentzAngleInput_(0)
142			{
143			// FIXME: Which granularity, leading to how many parameters?
144			parameters_.resize(2, 0.); // currently two parameters (ring1-4, 5-8), start value 0.
145			paramUncertainties_.resize(2, 0.); // dito for errors
146
147			edm::LogInfo("Alignment") << "@SUB=SiPixelLorentzAngleCalibration" << "Created with name "
148			<< this->name() << "',\n" << this->numParameters() << " parameters to be determined,"
149			<< "\nsaveToDB = " << saveToDB_
150			<< "\n outFileName = " << outFileName_
151			<< "\n N(merge files) = " << mergeFileNames_.size();
152			if (mergeFileNames_.size()) {
153			edm::LogInfo("Alignment") << "@SUB=SiPixelLorentzAngleCalibration"
154			<< "First file to merge: " << mergeFileNames_[0];
155			}
156			}
157
158			//======================================================================
159			SiPixelLorentzAngleCalibration::~SiPixelLorentzAngleCalibration()
160			{
161			// std::cout << "Destroy SiPixelLorentzAngleCalibration named " << this->name() << std::endl;
162			delete siPixelLorentzAngleInput_;
163			}
164
165			//======================================================================
166			unsigned int SiPixelLorentzAngleCalibration::numParameters() const
167			{
168			return parameters_.size();
169			}
170
171			//======================================================================
172			unsigned int
173			SiPixelLorentzAngleCalibration::derivatives(std::vector<ValuesIndexPair> &outDerivInds,
174			const TransientTrackingRecHit &hit,
175			const TrajectoryStateOnSurface &tsos,
176			const edm::EventSetup &setup,
177			const EventInfo &eventInfo) const
178			{
179			// ugly const-cast:
180			// But it is either only first initialisation or throwing an exception...
181			const_cast<SiPixelLorentzAngleCalibration*>(this)->checkLorentzAngleInput(setup, eventInfo);
182
183			outDerivInds.clear();
184
185			if (hit.det()) { // otherwise 'constraint hit' or whatever
186
187			const int index = this->getParameterIndexFromDetId(hit.det()->geographicalId(),
188			eventInfo.eventId_.run());
189			if (index >= 0) { // otherwise not treated
190			edm::ESHandle<MagneticField> magneticField;
191			setup.get<IdealMagneticFieldRecord>().get(magneticField);
192			const GlobalVector bField(magneticField->inTesla(hit.det()->surface().position()));
193			const LocalVector bFieldLocal(hit.det()->surface().toLocal(bField));
194			const double dZ = hit.det()->surface().bounds().thickness(); // it is a float only...
195			// shift due to LA: dx = tan(LA) * dz = mobility * B_y * dz
196			const double xDerivative = bFieldLocal.y() * dZ; // parameter is mobility!
197			if (xDerivative) { // If field is zero, this is zero: do not return it
198			const Values derivs(xDerivative, 0.); // yDerivative = 0.
199			outDerivInds.push_back(ValuesIndexPair(derivs, index));
200			}
201			}
202			} else {
203			edm::LogWarning("Alignment") << "@SUB=SiPixelLorentzAngleCalibration::derivatives2"
204			<< "Hit without GeomDet, skip!";
205			}
206
207			return outDerivInds.size();
208			}
209
210			//======================================================================
211			bool SiPixelLorentzAngleCalibration::setParameter(unsigned int index, double value)
212			{
213			if (index >= parameters_.size()) {
214			return false;
215			} else {
216			parameters_[index] = value;
217			return true;
218			}
219			}
220
221			//======================================================================
222			bool SiPixelLorentzAngleCalibration::setParameterError(unsigned int index, double error)
223			{
224			if (index >= paramUncertainties_.size()) {
225			return false;
226			} else {
227			paramUncertainties_[index] = error;
228			return true;
229			}
230			}
231
232			//======================================================================
233			double SiPixelLorentzAngleCalibration::getParameter(unsigned int index) const
234			{
235			// if (index >= parameters_.size()) {
236			// return 0.;
237			// } else {
238			// return parameters_[index];
239			// }
240			return (index >= parameters_.size() ? 0. : parameters_[index]);
241			}
242
243			//======================================================================
244			double SiPixelLorentzAngleCalibration::getParameterError(unsigned int index) const
245			{
246			// if (index >= paramUncertainties_.size()) {
247			// return 0.;
248			// } else {
249			// return paramUncertainties_[index];
250			// }
251			return (index >= paramUncertainties_.size() ? 0. : paramUncertainties_[index]);
252			}
253
254			// //======================================================================
255			// void SiPixelLorentzAngleCalibration::beginOfJob(const AlignableTracker *tracker,
256			// const AlignableMuon /muon*/,
257			// const AlignableExtras /extras*/)
258			// {
259			// alignableTracker_ = tracker;
260			// }
261
262
263			//======================================================================
264			void SiPixelLorentzAngleCalibration::endOfJob()
265			{
266			// loginfo output
267			std::ostringstream out;
268			out << "Parameter results\n";
269			for (unsigned int iPar = 0; iPar < parameters_.size(); ++iPar) {
270			out << iPar << ": " << parameters_[iPar] << " +- " << paramUncertainties_[iPar] << "\n";
271			}
272			edm::LogInfo("Alignment") << "@SUB=SiPixelLorentzAngleCalibration::endOfJob" << out.str();
273
274			// now write 'input' tree
275			const SiPixelLorentzAngle *input = this->getLorentzAnglesInput(); // never NULL
276			const std::string treeName(this->name() + '_');
277			this->writeTree(input, (treeName + "input").c_str());
278
279			for (unsigned int iIOV = 0; iIOV < this->numIovs(); ++iIOV) {
280			cond::Time_t firstRunOfIOV = this->firstRunOfIOV(iIOV);
281			SiPixelLorentzAngle *output = new SiPixelLorentzAngle;
282			// Loop on map of values from input and add (possible) parameter results
283			for (auto iterIdValue = input->getLorentzAngles().begin();
284			iterIdValue != input->getLorentzAngles().end(); ++iterIdValue) {
285			// type of (*iterIdValue) is pair<unsigned int, float>
286			const unsigned int detId = iterIdValue->first; // key of map is DetId
287			// Nasty: putLorentzAngle(..) takes float by reference - not even const reference!
288			float value = iterIdValue->second + this->getParameterForDetId(detId, firstRunOfIOV);
289			output->putLorentzAngle(detId, value); // put result in output
290			}
291
292			// Write this even for mille jobs?
293			this->writeTree(output, (treeName + Form("result_%lld", firstRunOfIOV)).c_str());
294
295			if (saveToDB_) { // If requested, write out to DB
296			edm::Service<cond::service::PoolDBOutputService> dbService;
297			if (dbService.isAvailable()) {
298			dbService->writeOne(output, firstRunOfIOV, recordNameDBwrite_.c_str());
299			// no 'delete output;': writeOne(..) took over ownership
300			} else {
301			delete output;
302			edm::LogError("BadConfig") << "@SUB=SiPixelLorentzAngleCalibration::endOfJob"
303			<< "No PoolDBOutputService available, but saveToDB true!";
304			}
305			} else {
306			delete output;
307			}
308			} // end loop on IOVs
309
310			}
311
312			//======================================================================
313			bool SiPixelLorentzAngleCalibration::checkLorentzAngleInput(const edm::EventSetup &setup,
314			const EventInfo &eventInfo)
315			{
316			edm::ESHandle<SiPixelLorentzAngle> lorentzAngleHandle;
317			if (!siPixelLorentzAngleInput_) {
318			setup.get<SiPixelLorentzAngleRcd>().get(lorentzAngleHandle);
319			siPixelLorentzAngleInput_ = new SiPixelLorentzAngle(*lorentzAngleHandle);
320			} else {
321			if (watchLorentzAngleRcd_.check(setup)) { // new IOV of input
322			setup.get<SiPixelLorentzAngleRcd>().get(lorentzAngleHandle);
323			if (lorentzAngleHandle->getLorentzAngles() // but only bad if non-identical values
324			!= siPixelLorentzAngleInput_->getLorentzAngles()) { // (comparing maps)
325			// FIXME: Could different maps have same content, but different order?
326			// Or 'floating point comparison' problem?
327			throw cms::Exception("BadInput")
328			<< "SiPixelLorentzAngleCalibration::checkLorentzAngleInput:\n"
329			<< "Content of SiPixelLorentzAngle changed at run " << eventInfo.eventId_.run()
330			<< ", but algorithm expects constant input!\n";
331			return false; // not reached...
332			}
333			}
334			}
335
336			return true;
337			}
338
339			//======================================================================
340			const SiPixelLorentzAngle* SiPixelLorentzAngleCalibration::getLorentzAnglesInput()
341			{
342			// For parallel processing in Millepede II, create SiPixelLorentzAngle
343			// from info stored in files of parallel jobs and check that they are identical.
344			// If this job has run on data, still check that LA is identical to the ones
345			// from mergeFileNames_.
346			const std::string treeName(this->name() + "_input");
347			for (auto iFile = mergeFileNames_.begin(); iFile != mergeFileNames_.end(); ++iFile) {
348			SiPixelLorentzAngle* la = this->createFromTree(iFile->c_str(), treeName.c_str());
349			// siPixelLorentzAngleInput_ could be non-null from previous file of this loop
350			// or from checkLorentzAngleInput(..) when running on data in this job as well
351			if (!siPixelLorentzAngleInput_ \|\| siPixelLorentzAngleInput_->getLorentzAngles().empty()) {
352			delete siPixelLorentzAngleInput_; // NULL or empty
353			siPixelLorentzAngleInput_ = la;
354			} else {
355			// FIXME: about comparison of maps see comments in checkLorentzAngleInput
356			if (!la->getLorentzAngles().empty() && // single job might not have got events
357			la->getLorentzAngles() != siPixelLorentzAngleInput_->getLorentzAngles()) {
358			// Throw exception instead of error?
359			edm::LogError("NoInput") << "@SUB=SiPixelLorentzAngleCalibration::getLorentzAnglesInput"
360			<< "Different input values from tree " << treeName
361			<< " in file " << *iFile << ".";
362
363			}
364			delete la;
365			}
366			}
367
368			if (!siPixelLorentzAngleInput_) { // no files nor ran on events
369			siPixelLorentzAngleInput_ = new SiPixelLorentzAngle;
370			edm::LogError("NoInput") << "@SUB=SiPixelLorentzAngleCalibration::getLorentzAnglesInput"
371			<< "No input, create an empty one!";
372			} else if (siPixelLorentzAngleInput_->getLorentzAngles().empty()) {
373			edm::LogError("NoInput") << "@SUB=SiPixelLorentzAngleCalibration::getLorentzAnglesInput"
374			<< "Empty result!";
375			}
376
377			return siPixelLorentzAngleInput_;
378			}
379
380			//======================================================================
381			double SiPixelLorentzAngleCalibration::getParameterForDetId(unsigned int detId,
382			edm::RunNumber_t run) const
383			{
384			const int index = this->getParameterIndexFromDetId(detId, run);
385
386			return (index < 0 ? 0. : parameters_[index]);
387			}
388
389			//======================================================================
390			int SiPixelLorentzAngleCalibration::getParameterIndexFromDetId(unsigned int detId,
391			edm::RunNumber_t run) const
392			{
393			// Return the index of the parameter that is used for this DetId.
394			// If this DetId is not treated, return values < 0.
395
396			// FIXME: Extend to configurable granularity?
397			// Including treatment of run dependence?
398			const PXBDetId id(detId);
399			if (id.det() == DetId::Tracker && id.subdetId() == PixelSubdetector::PixelBarrel) {
400			if (id.module() >= 1 && id.module() <= 4) return 0;
401			if (id.module() >= 5 && id.module() <= 8) return 1;
402			edm::LogWarning("Alignment")
403			<< "@SUB=SiPixelLorentzAngleCalibration::getParameterIndexFromDetId"
404			<< "Module should be 1-8, but is " << id.module() << " => skip!";
405			}
406
407			return -1;
408			}
409
410			//======================================================================
411			unsigned int SiPixelLorentzAngleCalibration::numIovs() const
412			{
413			// FIXME: Needed to include treatment of run dependence!
414			return 1;
415			}
416
417			//======================================================================
418			edm::RunNumber_t SiPixelLorentzAngleCalibration::firstRunOfIOV(unsigned int iovNum) const
419			{
420			// FIXME: Needed to include treatment of run dependence!
421			if (iovNum < this->numIovs()) return 1;
422			else return 0;
423			}
424
425
426			//======================================================================
427			void SiPixelLorentzAngleCalibration::writeTree(const SiPixelLorentzAngle *lorentzAngle,
428			const char *treeName) const
429			{
430			if (!lorentzAngle) return;
431
432			TFile* file = TFile::Open(outFileName_.c_str(), "UPDATE");
433			if (!file) {
434			edm::LogError("BadConfig") << "@SUB=SiPixelLorentzAngleCalibration::writeTree"
435			<< "Could not open file '" << outFileName_ << "'.";
436			return;
437			}
438
439			TTree *tree = new TTree(treeName, treeName);
440			unsigned int id = 0;
441			float value = 0.;
442			tree->Branch("detId", &id, "detId/i");
443			tree->Branch("value", &value, "value/F");
444
445			for (auto iterIdValue = lorentzAngle->getLorentzAngles().begin();
446			iterIdValue != lorentzAngle->getLorentzAngles().end(); ++iterIdValue) {
447			// type of (*iterIdValue) is pair<unsigned int, float>
448			id = iterIdValue->first; // key of map is DetId
449			value = iterIdValue->second;
450			tree->Fill();
451			}
452			tree->Write();
453			delete file; // tree vanishes with the file... (?)
454
455			}
456
457			//======================================================================
458			SiPixelLorentzAngle*
459			SiPixelLorentzAngleCalibration::createFromTree(const char fileName, const char treeName) const
460			{
461			TTree *tree = 0;
462			TFile* file = TFile::Open(fileName, "READ");
463			if (file) file->GetObject(treeName, tree);
464
465			SiPixelLorentzAngle *result = 0;
466			if (tree) {
467			unsigned int id = 0;
468			float value = 0.;
469			tree->SetBranchAddress("detId", &id);
470			tree->SetBranchAddress("value", &value);
471
472			result = new SiPixelLorentzAngle;
473			const Long64_t nEntries = tree->GetEntries();
474			for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
475			tree->GetEntry(iEntry);
476			result->putLorentzAngle(id, value);
477			}
478			} else {
479			edm::LogError("Alignment") << "@SUB=SiPixelLorentzAngleCalibration::createFromTree"
480			<< "Could not get TTree '" << treeName << "' from file '"
481			<< fileName << "'.";
482			}
483
484			delete file; // tree will vanish with file
485			return result;
486			}
487
488
489			//======================================================================
490			//======================================================================
491			// Plugin definition
492
493			#include "Alignment/CommonAlignmentAlgorithm/interface/IntegratedCalibrationPluginFactory.h"
494
495			DEFINE_EDM_PLUGIN(IntegratedCalibrationPluginFactory,
496			SiPixelLorentzAngleCalibration, "SiPixelLorentzAngleCalibration");