00001 #include "Validation/Geometry/interface/MaterialBudgetHcalHistos.h"
00002
00003 #include "DetectorDescription/Core/interface/DDFilter.h"
00004 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
00005 #include "DetectorDescription/Core/interface/DDSplit.h"
00006 #include "DetectorDescription/Core/interface/DDValue.h"
00007 #include "FWCore/Utilities/interface/Exception.h"
00008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
00011
00012 #include "CLHEP/Units/PhysicalConstants.h"
00013 #include "CLHEP/Units/SystemOfUnits.h"
00014
00015 MaterialBudgetHcalHistos::MaterialBudgetHcalHistos(const edm::ParameterSet &p){
00016
00017 binEta = p.getUntrackedParameter<int>("NBinEta", 260);
00018 binPhi = p.getUntrackedParameter<int>("NBinPhi", 180);
00019 maxEta = p.getUntrackedParameter<double>("MaxEta", 5.2);
00020 etaLow = p.getUntrackedParameter<double>("EtaLow", -5.2);
00021 etaHigh = p.getUntrackedParameter<double>("EtaHigh", 5.2);
00022 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Eta plot: NX "
00023 << binEta << " Range " << -maxEta << ":"
00024 << maxEta << " Phi plot: NX " << binPhi
00025 << " Range " << -pi << ":" << pi << " (Eta "
00026 << "limit " << etaLow << ":" << etaHigh <<")";
00027 book();
00028
00029 }
00030
00031 void MaterialBudgetHcalHistos::fillBeginJob(const DDCompactView & cpv) {
00032
00033 std::string attribute = "ReadOutName";
00034 std::string value = "HcalHits";
00035 DDSpecificsFilter filter1;
00036 DDValue ddv1(attribute,value,0);
00037 filter1.setCriteria(ddv1,DDSpecificsFilter::equals);
00038 DDFilteredView fv1(cpv);
00039 fv1.addFilter(filter1);
00040 sensitives = getNames(fv1);
00041 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
00042 << "tested for " << attribute << " = "
00043 << value << " has " << sensitives.size()
00044 << " elements";
00045 for (unsigned int i=0; i<sensitives.size(); i++)
00046 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: sensitives["
00047 << i << "] = " << sensitives[i];
00048
00049 attribute = "Volume";
00050 value = "HF";
00051 DDSpecificsFilter filter2;
00052 DDValue ddv2(attribute,value,0);
00053 filter2.setCriteria(ddv2,DDSpecificsFilter::equals);
00054 DDFilteredView fv2(cpv);
00055 fv2.addFilter(filter2);
00056 hfNames = getNames(fv2);
00057 fv2.firstChild();
00058 DDsvalues_type sv(fv2.mergedSpecifics());
00059 std::vector<double> temp = getDDDArray("Levels",sv);
00060 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
00061 << "tested for " << attribute << " = "
00062 << value << " has " << hfNames.size()
00063 << " elements";
00064 for (unsigned int i=0; i<hfNames.size(); i++) {
00065 int level = static_cast<int>(temp[i]);
00066 hfLevels.push_back(level);
00067 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: HF[" << i
00068 << "] = " << hfNames[i] << " at level "
00069 << hfLevels[i];
00070 }
00071
00072 std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
00073 attribute = "ReadOutName";
00074 for (int k=0; k<2; k++) {
00075 value = ecalRO[k];
00076 DDSpecificsFilter filter3;
00077 DDValue ddv3(attribute,value,0);
00078 filter3.setCriteria(ddv3,DDSpecificsFilter::equals);
00079 DDFilteredView fv3(cpv);
00080 fv3.addFilter(filter3);
00081 std::vector<std::string> senstmp = getNames(fv3);
00082 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
00083 << "tested for " << attribute << " = "
00084 << value << " has " << senstmp.size()
00085 << " elements";
00086 for (unsigned int i=0; i<senstmp.size(); i++)
00087 sensitiveEC.push_back(senstmp[i]);
00088 }
00089 for (unsigned int i=0; i<sensitiveEC.size(); i++)
00090 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: sensitiveEC["
00091 << i << "] = " << sensitiveEC[i];
00092 }
00093
00094 void MaterialBudgetHcalHistos::fillStartTrack(const G4Track* aTrack) {
00095
00096 id = layer = steps = 0;
00097 radLen = intLen = stepLen = 0;
00098
00099 const G4ThreeVector& dir = aTrack->GetMomentum() ;
00100 if (dir.theta() != 0 ) {
00101 eta = dir.eta();
00102 } else {
00103 eta = -99;
00104 }
00105 phi = dir.phi();
00106 double theEnergy = aTrack->GetTotalEnergy();
00107 int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
00108
00109 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: Track "
00110 << aTrack->GetTrackID() << " Code " << theID
00111 << " Energy " << theEnergy/GeV << " GeV; Eta "
00112 << eta << " Phi " << phi/deg << " PT "
00113 << dir.perp()/GeV << " GeV *****";
00114 }
00115
00116
00117 void MaterialBudgetHcalHistos::fillPerStep(const G4Step* aStep) {
00118
00119 G4Material * material = aStep->GetPreStepPoint()->GetMaterial();
00120 double step = aStep->GetStepLength();
00121 double radl = material->GetRadlen();
00122 double intl = material->GetNuclearInterLength();
00123 double density = material->GetDensity() / (g/cm3);
00124
00125 int idOld = id;
00126 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
00127 std::string name = touch->GetVolume(0)->GetName();
00128 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: Step at " << name
00129 << " Length " << step << " in "
00130 << material->GetName() << " of density "
00131 << density << " g/cc; Radiation Length "
00132 << radl << " mm; Interaction Length " << intl
00133 << " mm\n Position "
00134 << aStep->GetPreStepPoint()->GetPosition()
00135 << " Length (so far) " << stepLen << " L/X0 "
00136 << step/radl << "/" << radLen << " L/Lambda "
00137 << step/intl << "/" << intLen;
00138
00139 int det=0, lay=0;
00140 if (isItEC(name)) {
00141 det = 1;
00142 lay = 1;
00143 } else {
00144 if (isSensitive(name)) {
00145 if (isItHF(touch)) {
00146 det = 5;
00147 lay = 21;
00148 } else {
00149 det = (touch->GetReplicaNumber(1))/1000;
00150 lay = (touch->GetReplicaNumber(0)/10)%100 + 3;
00151 if (det == 4) {
00152 double abeta = std::abs(eta);
00153 if (abeta < 1.479) lay = layer + 1;
00154 else lay--;
00155 if (lay < 3) lay = 3;
00156 if (lay == layer) lay++;
00157 if (lay > 20) lay = 20;
00158 }
00159 }
00160 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: Det " << det
00161 << " Layer " << lay << " Eta " << eta
00162 << " Phi " << phi/deg;
00163 } else if (layer == 1) {
00164 det = -1;
00165 lay = 2;
00166 }
00167 }
00168 if (det != 0) {
00169 if (lay != layer) {
00170 id = lay;
00171 layer = lay;
00172 }
00173 }
00174
00175 if (id > idOld) {
00176
00177 fillHisto(id-1);
00178 }
00179
00180 stepLen += step;
00181 radLen += step/radl;
00182 intLen += step/intl;
00183 if (layer == 21 && det == 5) {
00184 if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
00185 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: After HF in "
00186 << aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName();
00187 fillHisto(id);
00188 id++;
00189 layer = 0;
00190 }
00191 }
00192 }
00193
00194
00195 void MaterialBudgetHcalHistos::fillEndTrack() {
00196 fillHisto(maxSet-1);
00197 }
00198
00199 void MaterialBudgetHcalHistos::book() {
00200
00201
00202 edm::Service<TFileService> tfile;
00203
00204 if ( !tfile.isAvailable() )
00205 throw cms::Exception("BadConfig") << "TFileService unavailable: "
00206 << "please add it to config file";
00207
00208 double maxPhi=pi;
00209 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user "
00210 << "histos === with " << binEta << " bins "
00211 << "in eta from " << -maxEta << " to "
00212 << maxEta << " and " << binPhi << " bins "
00213 << "in phi from " << -maxPhi << " to "
00214 << maxPhi;
00215
00216 char name[10], title[40];
00217
00218 for (int i=0; i<maxSet; i++) {
00219 sprintf(name, "%d", i+100);
00220 sprintf(title, "MB(X0) prof Eta in region %d", i);
00221 me100[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00222 sprintf(name, "%d", i+200);
00223 sprintf(title, "MB(L0) prof Eta in region %d", i);
00224 me200[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00225 sprintf(name, "%d", i+300);
00226 sprintf(title, "MB(Step) prof Eta in region %d", i);
00227 me300[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00228 sprintf(name, "%d", i+400);
00229 sprintf(title, "Eta in region %d", i);
00230 me400[i] = tfile->make<TH1F>(name, title, binEta, -maxEta, maxEta);
00231 sprintf(name, "%d", i+500);
00232 sprintf(title, "MB(X0) prof Ph in region %d", i);
00233 me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00234 sprintf(name, "%d", i+600);
00235 sprintf(title, "MB(L0) prof Ph in region %d", i);
00236 me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00237 sprintf(name, "%d", i+700);
00238 sprintf(title, "MB(Step) prof Ph in region %d", i);
00239 me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00240 sprintf(name, "%d", i+800);
00241 sprintf(title, "Phi in region %d", i);
00242 me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
00243 sprintf(name, "%d", i+900);
00244 sprintf(title, "MB(X0) prof Eta Phi in region %d", i);
00245 me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
00246 binPhi/2, -maxPhi, maxPhi);
00247 sprintf(name, "%d", i+1000);
00248 sprintf(title, "MB(L0) prof Eta Phi in region %d", i);
00249 me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
00250 binPhi/2, -maxPhi, maxPhi);
00251 sprintf(name, "%d", i+1100);
00252 sprintf(title, "MB(Step) prof Eta Phi in region %d", i);
00253 me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
00254 binPhi/2, -maxPhi, maxPhi);
00255 sprintf(name, "%d", i+1200);
00256 sprintf(title, "Eta vs Phi in region %d", i);
00257 me1200[i]= tfile->make<TH2F>(name, title, binEta/2, -maxEta, maxEta,
00258 binPhi/2, -maxPhi, maxPhi);
00259 }
00260
00261 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user "
00262 << "histos done ===";
00263
00264 }
00265
00266 void MaterialBudgetHcalHistos::fillHisto(int ii) {
00267
00268 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos:FillHisto called "
00269 << "with index " << ii << " integrated step "
00270 << stepLen << " X0 " << radLen << " Lamda "
00271 << intLen;
00272
00273 if (ii >=0 && ii < maxSet) {
00274 me100[ii]->Fill(eta, radLen);
00275 me200[ii]->Fill(eta, intLen);
00276 me300[ii]->Fill(eta, stepLen);
00277 me400[ii]->Fill(eta);
00278
00279 if (eta >= etaLow && eta <= etaHigh) {
00280 me500[ii]->Fill(phi, radLen);
00281 me600[ii]->Fill(phi, intLen);
00282 me700[ii]->Fill(phi, stepLen);
00283 me800[ii]->Fill(phi);
00284 }
00285
00286 me900[ii]->Fill(eta, phi, radLen);
00287 me1000[ii]->Fill(eta, phi, intLen);
00288 me1100[ii]->Fill(eta, phi, stepLen);
00289 me1200[ii]->Fill(eta, phi);
00290
00291 }
00292 }
00293
00294 void MaterialBudgetHcalHistos::hend() {
00295 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Save user "
00296 << "histos ===";
00297 }
00298
00299 std::vector<std::string> MaterialBudgetHcalHistos::getNames(DDFilteredView& fv) {
00300
00301 std::vector<std::string> tmp;
00302 bool dodet = fv.firstChild();
00303 while (dodet) {
00304 const DDLogicalPart & log = fv.logicalPart();
00305 std::string namx = DDSplit(log.name()).first;
00306 bool ok = true;
00307 for (unsigned int i=0; i<tmp.size(); i++)
00308 if (namx == tmp[i]) ok = false;
00309 if (ok) tmp.push_back(namx);
00310 dodet = fv.next();
00311 }
00312 return tmp;
00313 }
00314
00315 std::vector<double> MaterialBudgetHcalHistos::getDDDArray(const std::string & str,
00316 const DDsvalues_type & sv) {
00317
00318 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos:getDDDArray called "
00319 << "for " << str;
00320 DDValue value(str);
00321 if (DDfetch(&sv,value)) {
00322 LogDebug("MaterialBudget") << value;
00323 const std::vector<double> & fvec = value.doubles();
00324 int nval = fvec.size();
00325 if (nval < 1) {
00326 edm::LogError("MaterialBudget") << "MaterialBudgetHcalHistos : # of "
00327 << str << " bins " << nval
00328 << " < 1 ==> illegal";
00329 throw cms::Exception("Unknown", "MaterialBudgetHcalHistos")
00330 << "nval < 1 for array " << str << "\n";
00331 }
00332
00333 return fvec;
00334 } else {
00335 edm::LogError("MaterialBudget") << "MaterialBudgetHcalHistos : cannot get "
00336 << "array " << str;
00337 throw cms::Exception("Unknown", "MaterialBudgetHcalHistos")
00338 << "cannot get array " << str << "\n";
00339 }
00340 }
00341
00342 bool MaterialBudgetHcalHistos::isSensitive (std::string name) {
00343
00344 std::vector<std::string>::const_iterator it = sensitives.begin();
00345 for (; it != sensitives.end(); it++)
00346 if (name == *it) return true;
00347 return false;
00348 }
00349
00350 bool MaterialBudgetHcalHistos::isItHF (const G4VTouchable* touch) {
00351
00352 std::vector<std::string>::const_iterator it = hfNames.begin();
00353 int levels = ((touch->GetHistoryDepth())+1);
00354 for (unsigned int it=0; it < hfNames.size(); it++) {
00355 if (levels >= hfLevels[it]) {
00356 std::string name = touch->GetVolume(levels-hfLevels[it])->GetName();
00357 if (name == hfNames[it]) return true;
00358 }
00359 }
00360 return false;
00361 }
00362
00363 bool MaterialBudgetHcalHistos::isItEC (std::string name) {
00364
00365 std::vector<std::string>::const_iterator it = sensitiveEC.begin();
00366 for (; it != sensitiveEC.end(); it++)
00367 if (name == *it) return true;
00368 return false;
00369 }
00370
1.5.4