CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoRomanPot/RecoFP420/src/RecoProducerFP420.cc

Go to the documentation of this file.
00001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00002 #include "FWCore/ParameterSet/interface/FileInPath.h"
00003 
00004 #include "RecoRomanPot/RecoFP420/interface/RecoProducerFP420.h"
00005 #include "DataFormats/FP420Cluster/interface/RecoFP420.h"
00006 #include "DataFormats/FP420Cluster/interface/RecoCollectionFP420.h"
00007 
00008 #include "CLHEP/Vector/LorentzVector.h"
00009 
00010 #include <math.h>
00011 
00012 //#include <iostream>
00013 
00014 RecoProducerFP420::RecoProducerFP420(const edm::ParameterSet& conf):conf_(conf)  { 
00015 
00016   // Create LHC beam line
00017   //  double length  = param.getParameter<double>("BeamLineLength");
00018   //  std::cout << " BeamLineLength = " << length << std::endl;
00019 
00020   length         = conf_.getParameter<double>("BeamLineLength" );//m
00021   verbosity      = conf_.getUntrackedParameter<int>("VerbosityLevel");
00022   beam1filename  = conf_.getParameter<std::string>("Beam1");
00023   beam2filename  = conf_.getParameter<std::string>("Beam2");  
00024 
00025   edm::LogInfo ("RecoProducerFP420") << "RecoProducerFP420 parameters: \n" 
00026                                      << "   Verbosity: " << verbosity << "\n"
00027                                      << "   length:    " << length << "\n";
00028   if (verbosity > 1) {
00029     std::cout << "  RecoProducerFP420: constructor    " << std::endl;
00030     std::cout << "   BeamLineLength:    " << length << std::endl;
00031   }
00032   
00033   //  edm::FileInPath b1("SimTransport/HectorData/twiss_ip5_b1_v6.5.txt");
00034   //  edm::FileInPath b2("SimTransport/HectorData/twiss_ip5_b2_v6.5.txt");
00035   
00036   edm::FileInPath b1(beam1filename.c_str());
00037   edm::FileInPath b2(beam2filename.c_str());
00038   
00039   
00040   m_beamline1 = new H_BeamLine(  1, length + 0.1 ); // (direction, length)
00041   m_beamline2 = new H_BeamLine( -1, length + 0.1 ); //
00042   
00043   try {
00044     m_beamline1->fill( b1.fullPath(),  1, "IP5" );
00045     m_beamline2->fill( b2.fullPath(), -1, "IP5" );
00046   } catch ( const edm::Exception& e ) {
00047     std::string msg = e.what();
00048     msg += " caught in RecoProducerFP420... \nERROR: Could not locate SimTransport/HectorData data files.";
00049     edm::LogError ("DataNotFound") << msg;
00050   }
00051   
00052   m_beamline1->offsetElements( 120, -0.097 );
00053   m_beamline2->offsetElements( 120, +0.097 );
00054   
00055   m_beamline1->calcMatrix();
00056   m_beamline2->calcMatrix();
00057 
00058   edm::LogInfo ("RecoProducerFP420") << "==============================\n";
00059 
00060 }
00061 
00062 RecoProducerFP420::~RecoProducerFP420(){}
00063 
00064 std::vector<RecoFP420>  RecoProducerFP420::reconstruct(int direction, double x1_420, double y1_420, double x2_420, double y2_420, double z1_420, double z2_420){
00065   // ==================================================================================  
00066   // ==================================================================================  
00067   std::vector<RecoFP420> rhits;
00068   int restracks = 10;// max # tracks
00069   rhits.reserve(restracks); 
00070   rhits.clear();
00071   // ==================================================================================  
00072 // trivial (TM), angle compensation (AM) and position compensation (PM) methods
00073 // #define TM 1    #define AM 2   #define PM 3
00074   m_tx0=-100., m_ty0=-100., m_x0=-100., m_y0=-100.;
00075   if ( direction == 1 ) {
00076     m_rp420_f = new H_RecRPObject( z1_420*0.001, z2_420*0.001, *m_beamline1 );// m
00077     if ( m_rp420_f ) {
00078       if (verbosity >1) {
00079         std::cout << "  RecoProducerFP420: input coord. in um   " << std::endl;
00080         std::cout << "   x1_420:    " << x1_420 << "   y1_420:    " << y1_420 << std::endl;
00081         std::cout << "   x2_420:    " << x2_420 << "   y2_420:    " << y2_420 << std::endl;
00082       }
00083       m_rp420_f->setPositions( x1_420, y1_420, x2_420 ,y2_420 );//input coord. in um
00084       m_e   = m_rp420_f->getE( AM );// GeV
00085       //  std::cout << "   m_e1:    " << m_rp420_f->getE( TM ) << std::endl;
00086       //  std::cout << "   m_e2:    " << m_rp420_f->getE( AM ) << std::endl;
00087       m_tx0 = m_rp420_f->getTXIP();// urad
00088       m_ty0 = m_rp420_f->getTYIP();// urad
00089       m_x0  = m_rp420_f->getX0();// um
00090       m_y0  = m_rp420_f->getY0();// um
00091       m_q2  = m_rp420_f->getQ2();// GeV^2
00092     }// if ( m_rp420_f
00093   }// if ( dire
00094   else if ( direction == 2 ) {
00095     m_rp420_b = new H_RecRPObject( z1_420*0.001, z2_420*0.001, *m_beamline2 );// m
00096     if ( m_rp420_b ) {
00097       m_rp420_b->setPositions( x1_420, y1_420, x2_420 ,y2_420 );// input coord. in um
00098       m_e   = m_rp420_b->getE( AM );// GeV
00099       m_tx0 = m_rp420_b->getTXIP();// urad
00100       m_ty0 = m_rp420_b->getTYIP();// urad
00101       m_x0  = m_rp420_b->getX0();// um
00102       m_y0  = m_rp420_b->getY0();// um
00103       m_q2  = m_rp420_b->getQ2();// GeV^2
00104     }// if ( m_rp420_b
00105   }
00106   else{
00107     return rhits;
00108   }
00109   
00110   // ==============================  
00111   if (verbosity > 1) {
00112     std::cout << "  RecoProducerFP420: rhits.push_back    " << std::endl;
00113     std::cout << "   m_e:    " << m_e << std::endl;
00114     std::cout << "   m_x0:    " << m_x0 << std::endl;
00115     std::cout << "   m_y0:    " << m_y0 << std::endl;
00116     std::cout << "   m_tx0:    " << m_tx0  << std::endl;
00117     std::cout << "   m_ty0:    " << m_ty0  << std::endl;
00118     std::cout << "   m_q2:    " << m_q2  << std::endl;
00119     std::cout << "   direction:    " << direction  << std::endl;
00120   }
00121   rhits.push_back( RecoFP420(m_e,m_x0,m_y0,m_tx0,m_ty0,m_q2,direction) );
00122   // ==============================  
00124     return rhits;
00125     //============
00126 }
00127 
00128 
00129