FP420RecoMain.cc

Go to the documentation of this file.

// File: FP420RecoMain.cc
// Date: 12.2006
// Description: FP420RecoMain for FP420
// Modifications:
#include <vector>
#include <iostream>
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
 
#include "RecoRomanPot/RecoFP420/interface/FP420RecoMain.h"
 
#include "DataFormats/FP420Cluster/interface/TrackFP420.h"
#include "DataFormats/FP420Cluster/interface/TrackCollectionFP420.h"
#include "DataFormats/FP420Cluster/interface/RecoFP420.h"
#include "DataFormats/FP420Cluster/interface/RecoCollectionFP420.h"
#include "RecoRomanPot/RecoFP420/interface/RecoProducerFP420.h"
 
// #include "CLHEP/Vector/LorentzVector.h"
// #include "CLHEP/Random/RandFlat.h"
#include <cmath>
 
using namespace std;
 
FP420RecoMain::FP420RecoMain(const edm::ParameterSet &conf) : conf_(conf) {
  verbosity = conf_.getUntrackedParameter<int>("VerbosityLevel");
  m_rpp420_f = conf_.getParameter<double>("RP420f");  //mm
  m_rpp420_b = conf_.getParameter<double>("RP420b");  //mm
  m_zreff = conf_.getParameter<double>("zreff");      //mm
  m_zrefb = conf_.getParameter<double>("zrefb");      //mm
  dn0 = conf_.getParameter<int>("NumberFP420Detectors");
 
  if (verbosity > 0) {
    std::cout << "FP420RecoMain constructor::" << std::endl;
    std::cout << "m_rpp420_f=" << m_rpp420_f << " m_rpp420_b=" << m_rpp420_b << std::endl;
    std::cout << "m_zreff=" << m_zreff << " m_zrefb=" << m_zrefb << std::endl;
  }
 
  double eee1 = 11.;
  double eee2 = 12.;
  //  zinibeg_ = (eee1-eee2)/2.;
  zinibeg_ = 0.;
  //
  if (verbosity > 1) {
    std::cout << "FP420RecoMain constructor::" << std::endl;
    std::cout << " eee1=" << eee1 << " eee2=" << eee2 << " zinibeg =" << zinibeg_ << std::endl;
  }
  finderParameters_ = new RecoProducerFP420(conf_);
}
 
FP420RecoMain::~FP420RecoMain() {
  if (finderParameters_ != nullptr) {
    delete finderParameters_;
  }
}
 
void FP420RecoMain::run(
    edm::Handle<TrackCollectionFP420> &input, RecoCollectionFP420 *toutput, double VtxX, double VtxY, double VtxZ) {
  // initialization
  bool first = true;
  //  finderParameters_->clear();
  // finderParameters_->setIP( 0., 0., 0. );
  std::vector<TrackFP420> rhits;
  int restracks = 10;  // max # tracks
  rhits.reserve(restracks);
  rhits.clear();
 
  // loop over detunits:
  for (int number_detunits = 1; number_detunits < dn0; number_detunits++) {
    unsigned int StID = number_detunits;
    std::vector<RecoFP420> rcollector;
    int restracks = 10;  // max # tracks
    rcollector.reserve(restracks);
    rcollector.clear();
 
    std::vector<TrackFP420> collector;
    collector.clear();
    TrackCollectionFP420::Range outputRange;
    unsigned int StIDTrack = 1111;
    double z420 = m_rpp420_f;
    double zref1 = m_zreff;
    double zinibeg = zinibeg_;
    double VtxXcur = VtxX;
    double VtxYcur = VtxY;
    double VtxZcur = VtxZ;
    if (StID == 2) {
      StIDTrack = 2222;
      z420 = -m_rpp420_b;
      zref1 = -m_zrefb;
      zinibeg = -zinibeg_;
      //  VtxXcur = -VtxX;
      // VtxYcur = -VtxY;
      // VtxZcur = -VtxZ;
    }
    double z1 = z420 + zinibeg - VtxZcur;
    double z2 = z420 + zinibeg + zref1 - VtxZcur;
    if (verbosity > 1) {
      std::cout << "FP420RecoMain: StIDTrack=" << StIDTrack << std::endl;
    }
    outputRange = input->get(StIDTrack);
    //
    // fill output in collector vector (for may be sorting? or other checks)
    //
    TrackCollectionFP420::ContainerIterator sort_begin = outputRange.first;
    TrackCollectionFP420::ContainerIterator sort_end = outputRange.second;
    //
    for (; sort_begin != sort_end; ++sort_begin) {
      collector.push_back(*sort_begin);
    }  // for  sort_begin
    if (verbosity > 1) {
      std::cout << "FP420RecoMain: track collector.size=" << collector.size() << std::endl;
    }
    std::vector<TrackFP420>::const_iterator simHitIter = collector.begin();
    std::vector<TrackFP420>::const_iterator simHitIterEnd = collector.end();
    for (; simHitIter != simHitIterEnd; ++simHitIter) {
      const TrackFP420 itrack = *simHitIter;
      double x1 = (itrack.bx() * z1 + (itrack.ax() - VtxXcur)) * 1000.;  //um
      double y1 = (itrack.by() * z1 + (itrack.ay() - VtxYcur)) * 1000.;  //um
      double x2 = (itrack.bx() * z2 + (itrack.ax() - VtxXcur)) * 1000.;  //um
      double y2 = (itrack.by() * z2 + (itrack.ay() - VtxYcur)) * 1000.;  //um
      if (verbosity == -49) {
        std::cout << "==================================================================== " << std::endl;
        std::cout << "FP420RecoMain: StID= " << StID << std::endl;
        std::cout << "input coord. in mm:  z1= " << z1 << std::endl;
        std::cout << "input coord. in mm:  z2= " << z2 << std::endl;
        std::cout << "input:  itrack.bx()= " << itrack.bx() << std::endl;
        std::cout << "input:  itrack.ax()= " << itrack.ax() << std::endl;
        std::cout << "input:  itrack.by()= " << itrack.by() << std::endl;
        std::cout << "input:  itrack.ay()= " << itrack.ay() << std::endl;
 
        std::cout << "input: in um X1noVtx= " << (itrack.bx() * (z420 + zinibeg) + itrack.ax()) * 1000. << std::endl;
        std::cout << "input: in um Y1noVtx= " << (itrack.by() * (z420 + zinibeg) + itrack.ay()) * 1000. << std::endl;
        std::cout << "input: in um X2noVtx= " << (itrack.bx() * (z420 + zinibeg + zref1) + itrack.ax()) * 1000.
                  << std::endl;
        std::cout << "input: in um Y2noVtx= " << (itrack.by() * (z420 + zinibeg + zref1) + itrack.ay()) * 1000.
                  << std::endl;
 
        std::cout << "input:  in mm VtxXcur= " << VtxXcur << std::endl;
        std::cout << "input:  in mm VtxYcur= " << VtxYcur << std::endl;
        std::cout << "input:  in mm VtxZcur= " << VtxZcur << std::endl;
        std::cout << "input coord. in um:  x1= " << x1 << std::endl;
        std::cout << "input coord. in um:  y1= " << y1 << std::endl;
        std::cout << "input coord. in um:  x2= " << x2 << std::endl;
        std::cout << "input coord. in um:  y2= " << y2 << std::endl;
      }
      double zz1 = fabs(z1);
      double zz2 = fabs(z2);
      rcollector = finderParameters_->reconstruct(
          StID, x1, y1, x2, y2, zz1, zz2);  // x1,y1,x2,y2 input coord. in um; z1, z2 in mm
    }
 
    if (verbosity > 1) {
      std::cout << "FP420RecoMain: track rcollector.size=" << rcollector.size() << std::endl;
    }
    if (!rcollector.empty()) {
      RecoCollectionFP420::Range rinputRange;
      rinputRange.first = rcollector.begin();
      rinputRange.second = rcollector.end();
 
      if (first) {
        // use it only if RecoCollectionFP420 is the RecoCollection of one event, otherwise, do not use (loose 1st cl. of 1st event only)
        first = false;
        unsigned int StID0 = 0;
        toutput->put(rinputRange, StID0);  // !!! put into adress 0 for detID which will not be used never
        if (verbosity > 1) {
          std::cout << "FP420RecoMain: put(rinputRange,StID0)" << std::endl;
        }
      }  //if ( first )
 
      // !!! put                                        !!! put
      toutput->put(rinputRange, StID);
      if (verbosity > 1) {
        std::cout << "FP420RecoMain: put(rinputRange,StID)" << std::endl;
      }
 
    }  // if rcollector.size
 
  }  // for loop over detunits
 
  if (verbosity > 1) {
    //     check of access to the zcollector:
    for (int number_detunits = 1; number_detunits < dn0; number_detunits++) {
      int StID = number_detunits;
      if (verbosity > 1) {
        std::cout << " ===" << std::endl;
        std::cout << " ===" << std::endl;
        std::cout << "FP420RecoMain:   re-new  StID= " << StID << std::endl;
      }
      std::vector<RecoFP420> zcollector;
      zcollector.clear();
      RecoCollectionFP420::Range zoutputRange;
      zoutputRange = toutput->get(StID);
      // fill output in zcollector vector (for may be sorting? or other checks)
      RecoCollectionFP420::ContainerIterator sort_begin = zoutputRange.first;
      RecoCollectionFP420::ContainerIterator sort_end = zoutputRange.second;
      for (; sort_begin != sort_end; ++sort_begin) {
        zcollector.push_back(*sort_begin);
      }  // for
      std::cout << "=======FP420RecoMain:check of re-new zcollector size = " << zcollector.size() << std::endl;
      std::cout << " ===" << std::endl;
      std::cout << " ===" << std::endl;
      std::vector<RecoFP420>::const_iterator simHitIter = zcollector.begin();
      std::vector<RecoFP420>::const_iterator simHitIterEnd = zcollector.end();
      // loop in recoess
      for (; simHitIter != simHitIterEnd; ++simHitIter) {
        const RecoFP420 itrack = *simHitIter;
        //double e0, double x0, double y0, double tx0, double ty0, double q20, int direction
        std::cout << "FP420RecoMain:check: direction = " << itrack.direction() << std::endl;
        std::cout << "  e0 = " << itrack.e0() << "  q20 = " << itrack.q20() << std::endl;
        std::cout << "  tx0 = " << itrack.tx0() << "  ty0 = " << itrack.ty0() << std::endl;
        std::cout << " x0= " << itrack.x0() << " y0= " << itrack.y0() << std::endl;
        std::cout << " ===" << std::endl;
        std::cout << " ===" << std::endl;
        std::cout << " =======================" << std::endl;
      }
    }
 
    //==================================
 
    //     end of check of access to the strip collection
    std::cout << "=======            FP420RecoMain:                    end of check     " << std::endl;
 
  }  //if (verbos
}