Go to the documentation of this file.00001
00002
00003
00005
00006 #include "SimG4CMS/Calo/interface/HFFibre.h"
00007
00008 #include "DetectorDescription/Core/interface/DDFilter.h"
00009 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00010 #include "DetectorDescription/Core/interface/DDValue.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012
00013 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00014 #include <iostream>
00015
00016
00017
00018 HFFibre::HFFibre(std::string & name, const DDCompactView & cpv,
00019 edm::ParameterSet const & p) {
00020
00021 edm::ParameterSet m_HF = p.getParameter<edm::ParameterSet>("HFShower");
00022 cFibre = c_light*(m_HF.getParameter<double>("CFibre"));
00023
00024 edm::LogInfo("HFShower") << "HFFibre:: Speed of light in fibre " << cFibre
00025 << " m/ns";
00026
00027 std::string attribute = "Volume";
00028 std::string value = "HF";
00029 DDSpecificsFilter filter1;
00030 DDValue ddv1(attribute,value,0);
00031 filter1.setCriteria(ddv1,DDSpecificsFilter::equals);
00032 DDFilteredView fv1(cpv);
00033 fv1.addFilter(filter1);
00034 bool dodet = fv1.firstChild();
00035
00036 if (dodet) {
00037 DDsvalues_type sv(fv1.mergedSpecifics());
00038
00039
00040 nBinAtt = -1;
00041 attL = getDDDArray("attl",sv,nBinAtt);
00042 edm::LogInfo("HFShower") << "HFFibre: " << nBinAtt << " attL ";
00043 for (int it=0; it<nBinAtt; it++)
00044 edm::LogInfo("HFShower") << "HFFibre: attL[" << it << "] = "
00045 << attL[it]*cm << "(1/cm)";
00046
00047
00048 int nb = 2;
00049 std::vector<double> nvec = getDDDArray("lambLim",sv,nb);
00050 lambLim[0] = static_cast<int>(nvec[0]);
00051 lambLim[1] = static_cast<int>(nvec[1]);
00052 edm::LogInfo("HFShower") << "HFFibre: Limits on lambda " << lambLim[0]
00053 << " and " << lambLim[1];
00054
00055
00056 nb = 0;
00057 longFL = getDDDArray("LongFL",sv,nb);
00058 edm::LogInfo("HFShower") << "HFFibre: " << nb << " Long Fibre Length";
00059 for (int it=0; it<nb; it++)
00060 edm::LogInfo("HFShower") << "HFFibre: longFL[" << it << "] = "
00061 << longFL[it]/cm << " cm";
00062 nb = 0;
00063 shortFL = getDDDArray("ShortFL",sv,nb);
00064 edm::LogInfo("HFShower") << "HFFibre: " << nb << " Short Fibre Length";
00065 for (int it=0; it<nb; it++)
00066 edm::LogInfo("HFShower") << "HFFibre: shortFL[" << it << "] = "
00067 << shortFL[it]/cm << " cm";
00068 } else {
00069 edm::LogError("HFShower") << "HFFibre: cannot get filtered "
00070 << " view for " << attribute << " matching "
00071 << name;
00072 throw cms::Exception("Unknown", "HFFibre")
00073 << "cannot match " << attribute << " to " << name <<"\n";
00074 }
00075
00076
00077 attribute = "ReadOutName";
00078 value = name;
00079 DDSpecificsFilter filter2;
00080 DDValue ddv2(attribute,value,0);
00081 filter2.setCriteria(ddv2,DDSpecificsFilter::equals);
00082 DDFilteredView fv2(cpv);
00083 fv2.addFilter(filter2);
00084 dodet = fv2.firstChild();
00085 if (dodet) {
00086 DDsvalues_type sv(fv2.mergedSpecifics());
00087
00088
00089 int nb = -1;
00090 gpar = getDDDArray("gparHF",sv,nb);
00091 edm::LogInfo("HFShower") << "HFFibre: " << nb <<" gpar (cm)";
00092 for (int i=0; i<nb; i++)
00093 edm::LogInfo("HFShower") << "HFFibre: gpar[" << i << "] = "
00094 << gpar[i]/cm << " cm";
00095
00096 nBinR = -1;
00097 radius = getDDDArray("rTable",sv,nBinR);
00098 edm::LogInfo("HFShower") << "HFFibre: " << nBinR <<" rTable (cm)";
00099 for (int i=0; i<nBinR; i++)
00100 edm::LogInfo("HFShower") << "HFFibre: radius[" << i << "] = "
00101 << radius[i]/cm << " cm";
00102 } else {
00103 edm::LogError("HFShower") << "HFFibre: cannot get filtered "
00104 << " view for " << attribute << " matching "
00105 << name;
00106 throw cms::Exception("Unknown", "HFFibre")
00107 << "cannot match " << attribute << " to " << name <<"\n";
00108 }
00109 }
00110
00111 HFFibre::~HFFibre() {}
00112
00113 double HFFibre::attLength(double lambda) {
00114
00115 int i = int(nBinAtt*(lambda - lambLim[0])/(lambLim[1]-lambLim[0]));
00116
00117 int j =i;
00118 if (i >= nBinAtt)
00119 j = nBinAtt-1;
00120 else if (i < 0)
00121 j = 0;
00122 double att = attL[j];
00123 #ifdef DebugLog
00124 edm::LogInfo("HFShower") << "HFFibre::attLength for Lambda " << lambda
00125 << " index " << i << " " << j << " Att. Length "
00126 << att;
00127 #endif
00128 return att;
00129 }
00130
00131 double HFFibre::tShift(G4ThreeVector point, int depth, int fromEndAbs) {
00132
00133 double zFibre = zShift(point, depth, fromEndAbs);
00134 double time = zFibre/cFibre;
00135 #ifdef DebugLog
00136 edm::LogInfo("HFShower") << "HFFibre::tShift for point " << point
00137 << " ( depth = " << depth <<", traversed length = "
00138 << zFibre/cm << " cm) = " << time/ns << " ns";
00139 #endif
00140 return time;
00141 }
00142
00143 double HFFibre::zShift(G4ThreeVector point, int depth, int fromEndAbs) {
00144
00145 double zFibre = 0;
00146 double hR = sqrt((point.x())*(point.x())+(point.y())*(point.y()));
00147 int ieta = 0;
00148 double length = 250*cm;
00149 if (fromEndAbs < 0) {
00150 zFibre = 0.5*gpar[1] - point.z();
00151 } else {
00152 for (int i = nBinR-1; i > 0; i--)
00153 if (hR < radius[i]) ieta = nBinR - i - 1;
00154 if (depth == 2) {
00155 if ((int)(shortFL.size()) > ieta) length = shortFL[ieta];
00156 } else {
00157 if ((int)(longFL.size()) > ieta) length = longFL[ieta];
00158 }
00159 zFibre = length;
00160 if (fromEndAbs > 0) {
00161 zFibre -= gpar[1];
00162 } else {
00163 double zz = 0.5*gpar[1] + point.z();
00164 zFibre -= zz;
00165 }
00166 if (depth == 2) zFibre += gpar[0];
00167 }
00168
00169 #ifdef DebugLog
00170 edm::LogInfo("HFShower") << "HFFibre::zShift for point " << point
00171 << " (R = " << hR/cm << " cm, Index = " << ieta
00172 << ", depth = " << depth << ", Fibre Length = "
00173 << length/cm << " cm = " << zFibre/cm
00174 << " cm)";
00175 #endif
00176 return zFibre;
00177 }
00178
00179 std::vector<double> HFFibre::getDDDArray(const std::string & str,
00180 const DDsvalues_type & sv,
00181 int & nmin) {
00182
00183 #ifdef DebugLog
00184 LogDebug("HFShower") << "HFFibre:getDDDArray called for " << str
00185 << " with nMin " << nmin;
00186 #endif
00187 DDValue value(str);
00188 if (DDfetch(&sv,value)) {
00189 #ifdef DebugLog
00190 LogDebug("HFShower") << value;
00191 #endif
00192 const std::vector<double> & fvec = value.doubles();
00193 int nval = fvec.size();
00194 if (nmin > 0) {
00195 if (nval < nmin) {
00196 edm::LogError("HFShower") << "HFFibre : # of " << str << " bins "
00197 << nval << " < " << nmin << " ==> illegal";
00198 throw cms::Exception("Unknown", "HFFibre")
00199 << "nval < nmin for array " << str <<"\n";
00200 }
00201 } else {
00202 if (nval < 1 && nmin != 0) {
00203 edm::LogError("HFShower") << "HFFibre : # of " << str << " bins "
00204 << nval << " < 1 ==> illegal (nmin="
00205 << nmin << ")";
00206 throw cms::Exception("Unknown", "HFFibre")
00207 << "nval < 1 for array " << str <<"\n";
00208 }
00209 }
00210 nmin = nval;
00211 return fvec;
00212 } else {
00213 if (nmin != 0) {
00214 edm::LogError("HFShower") << "HFFibre : cannot get array " << str;
00215 throw cms::Exception("Unknown", "HFFibre")
00216 << "cannot get array " << str <<"\n";
00217 } else {
00218 std::vector<double> fvec;
00219 return fvec;
00220 }
00221 }
00222 }