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 "CommonTools/UtilAlgos/interface/TFileService.h"
00011
00012 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00013 #include "CLHEP/Units/GlobalSystemOfUnits.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 fillHistos = p.getUntrackedParameter<bool>("FillHisto", true);
00023 printSum = p.getUntrackedParameter<bool>("PrintSummary", false);
00024 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: FillHisto : "
00025 << fillHistos << " PrintSummary " << printSum
00026 << " == Eta plot: NX " << binEta << " Range "
00027 << -maxEta << ":" << maxEta <<" Phi plot: NX "
00028 << binPhi << " Range " << -pi << ":" << pi
00029 << " (Eta limit " << etaLow << ":" << etaHigh
00030 <<")";
00031 if (fillHistos) book();
00032
00033 }
00034
00035 void MaterialBudgetHcalHistos::fillBeginJob(const DDCompactView & cpv) {
00036
00037 if (fillHistos) {
00038 std::string attribute = "ReadOutName";
00039 std::string value = "HcalHits";
00040 DDSpecificsFilter filter1;
00041 DDValue ddv1(attribute,value,0);
00042 filter1.setCriteria(ddv1,DDSpecificsFilter::equals);
00043 DDFilteredView fv1(cpv);
00044 fv1.addFilter(filter1);
00045 sensitives = getNames(fv1);
00046 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
00047 << "tested for " << attribute << " = "
00048 << value << " has " << sensitives.size()
00049 << " elements";
00050 for (unsigned int i=0; i<sensitives.size(); i++)
00051 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: sensitives["
00052 << i << "] = " << sensitives[i];
00053
00054 attribute = "Volume";
00055 value = "HF";
00056 DDSpecificsFilter filter2;
00057 DDValue ddv2(attribute,value,0);
00058 filter2.setCriteria(ddv2,DDSpecificsFilter::equals);
00059 DDFilteredView fv2(cpv);
00060 fv2.addFilter(filter2);
00061 hfNames = getNames(fv2);
00062 fv2.firstChild();
00063 DDsvalues_type sv(fv2.mergedSpecifics());
00064 std::vector<double> temp = getDDDArray("Levels",sv);
00065 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be "
00066 << "tested for " << attribute << " = "
00067 << value << " has " << hfNames.size()
00068 << " elements";
00069 for (unsigned int i=0; i<hfNames.size(); i++) {
00070 int level = static_cast<int>(temp[i]);
00071 hfLevels.push_back(level);
00072 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: HF[" << i
00073 << "] = " << hfNames[i] << " at level "
00074 << hfLevels[i];
00075 }
00076
00077 std::string ecalRO[2] = {"EcalHitsEB", "EcalHitsEE"};
00078 attribute = "ReadOutName";
00079 for (int k=0; k<2; k++) {
00080 value = ecalRO[k];
00081 DDSpecificsFilter filter3;
00082 DDValue ddv3(attribute,value,0);
00083 filter3.setCriteria(ddv3,DDSpecificsFilter::equals);
00084 DDFilteredView fv3(cpv);
00085 fv3.addFilter(filter3);
00086 std::vector<std::string> senstmp = getNames(fv3);
00087 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Names to be"
00088 << " tested for " << attribute << " = "
00089 << value << " has " << senstmp.size()
00090 << " elements";
00091 for (unsigned int i=0; i<senstmp.size(); i++)
00092 sensitiveEC.push_back(senstmp[i]);
00093 }
00094 for (unsigned int i=0; i<sensitiveEC.size(); i++)
00095 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos:sensitiveEC["
00096 << i << "] = " << sensitiveEC[i];
00097 }
00098 }
00099
00100 void MaterialBudgetHcalHistos::fillStartTrack(const G4Track* aTrack) {
00101
00102 id = layer = steps = 0;
00103 radLen = intLen = stepLen = 0;
00104 nlayHB = nlayHE = nlayHF = nlayHO = 0;
00105
00106 const G4ThreeVector& dir = aTrack->GetMomentum() ;
00107 if (dir.theta() != 0 ) {
00108 eta = dir.eta();
00109 } else {
00110 eta = -99;
00111 }
00112 phi = dir.phi();
00113 double theEnergy = aTrack->GetTotalEnergy();
00114 int theID = (int)(aTrack->GetDefinition()->GetPDGEncoding());
00115
00116 if (printSum) {
00117 matList.clear();
00118 stepLength.clear();
00119 radLength.clear();
00120 intLength.clear();
00121 }
00122
00123 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Track "
00124 << aTrack->GetTrackID() << " Code " << theID
00125 << " Energy " << theEnergy/GeV << " GeV; Eta "
00126 << eta << " Phi " << phi/deg << " PT "
00127 << dir.perp()/GeV << " GeV *****";
00128 }
00129
00130
00131 void MaterialBudgetHcalHistos::fillPerStep(const G4Step* aStep) {
00132
00133 G4Material * material = aStep->GetPreStepPoint()->GetMaterial();
00134 double step = aStep->GetStepLength();
00135 double radl = material->GetRadlen();
00136 double intl = material->GetNuclearInterLength();
00137 double density = material->GetDensity() / (g/cm3);
00138
00139 int idOld = id;
00140 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
00141 std::string name = touch->GetVolume(0)->GetName();
00142 std::string matName = material->GetName();
00143 if (printSum) {
00144 bool found = false;
00145 for (unsigned int ii=0; ii<matList.size(); ii++) {
00146 if (matList[ii] == matName) {
00147 stepLength[ii] += step;
00148 radLength[ii] += (step/radl);
00149 intLength[ii] += (step/intl);
00150 found = true;
00151 break;
00152 }
00153 }
00154 if (!found) {
00155 matList.push_back(matName);
00156 stepLength.push_back(step);
00157 radLength.push_back(step/radl);
00158 intLength.push_back(step/intl);
00159 }
00160 edm::LogInfo("MaterialBudget") << name << " " << step << " " << matName
00161 << " " << stepLen << " " << step/radl << " "
00162 << radLen << " " <<step/intl << " " <<intLen;
00163 } else {
00164 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Step at "
00165 << name << " Length " << step << " in "
00166 << matName << " of density " << density
00167 << " g/cc; Radiation Length " <<radl <<" mm;"
00168 << " Interaction Length " << intl << " mm\n"
00169 << " Position "
00170 << aStep->GetPreStepPoint()->GetPosition()
00171 << " Cylindrical R "
00172 <<aStep->GetPreStepPoint()->GetPosition().perp()
00173 << " Length (so far) " << stepLen << " L/X0 "
00174 << step/radl << "/" << radLen << " L/Lambda "
00175 << step/intl << "/" << intLen;
00176 }
00177
00178 int det=0, lay=0;
00179 if (fillHistos) {
00180 if (isItEC(name)) {
00181 det = 1;
00182 lay = 1;
00183 } else {
00184 if (isSensitive(name)) {
00185 if (isItHF(touch)) {
00186 det = 5;
00187 lay = 21;
00188 if (lay != layer) nlayHF++;
00189 } else {
00190 det = (touch->GetReplicaNumber(1))/1000;
00191 lay = (touch->GetReplicaNumber(0)/10)%100 + 3;
00192 if (det == 4) {
00193 double abeta = std::abs(eta);
00194 if (abeta < 1.479) lay = layer + 1;
00195 else lay--;
00196 if (lay < 3) lay = 3;
00197 if (lay == layer) lay++;
00198 if (lay > 20) lay = 20;
00199 if (lay != layer) nlayHE++;
00200 } else if (lay != layer) {
00201 if (lay >= 20) nlayHO++;
00202 else nlayHB++;
00203 }
00204 }
00205 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: Det " << det
00206 << " Layer " << lay << " Eta " << eta
00207 << " Phi " << phi/deg;
00208 } else if (layer == 1) {
00209 det = -1;
00210 lay = 2;
00211 }
00212 }
00213 if (det != 0) {
00214 if (lay != layer) {
00215 id = lay;
00216 layer = lay;
00217 }
00218 }
00219
00220 if (id > idOld) {
00221
00222 fillHisto(id-1);
00223 }
00224 }
00225
00226 stepLen += step;
00227 radLen += step/radl;
00228 intLen += step/intl;
00229 if (fillHistos) {
00230 if (layer == 21 && det == 5) {
00231 if (!isItHF(aStep->GetPostStepPoint()->GetTouchable())) {
00232 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos: After HF in "
00233 << aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)->GetName();
00234 fillHisto(id);
00235 id++;
00236 layer = 0;
00237 }
00238 }
00239 }
00240 }
00241
00242
00243 void MaterialBudgetHcalHistos::fillEndTrack() {
00244 edm::LogInfo("MaterialBudget") << "Number of layers hit in HB:" << nlayHB
00245 << " HE:" << nlayHE << " HO:" << nlayHO
00246 << " HF:" << nlayHF;
00247 if (fillHistos) {
00248 fillHisto(maxSet-1);
00249 fillLayer();
00250 }
00251 if (printSum) {
00252 for (unsigned int ii=0; ii<matList.size(); ii++) {
00253 edm::LogInfo("MaterialBudget") << matList[ii] << "\t" << stepLength[ii]
00254 << "\t" << radLength[ii] << "\t"
00255 << intLength[ii];
00256 }
00257 }
00258 }
00259
00260 void MaterialBudgetHcalHistos::book() {
00261
00262
00263 edm::Service<TFileService> tfile;
00264
00265 if ( !tfile.isAvailable() )
00266 throw cms::Exception("BadConfig") << "TFileService unavailable: "
00267 << "please add it to config file";
00268
00269 double maxPhi=pi;
00270 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user "
00271 << "histos === with " << binEta << " bins "
00272 << "in eta from " << -maxEta << " to "
00273 << maxEta << " and " << binPhi << " bins "
00274 << "in phi from " << -maxPhi << " to "
00275 << maxPhi;
00276
00277 char name[10], title[40];
00278
00279 for (int i=0; i<maxSet; i++) {
00280 sprintf(name, "%d", i+100);
00281 sprintf(title, "MB(X0) prof Eta in region %d", i);
00282 me100[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00283 sprintf(name, "%d", i+200);
00284 sprintf(title, "MB(L0) prof Eta in region %d", i);
00285 me200[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00286 sprintf(name, "%d", i+300);
00287 sprintf(title, "MB(Step) prof Eta in region %d", i);
00288 me300[i] = tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00289 sprintf(name, "%d", i+400);
00290 sprintf(title, "Eta in region %d", i);
00291 me400[i] = tfile->make<TH1F>(name, title, binEta, -maxEta, maxEta);
00292 sprintf(name, "%d", i+500);
00293 sprintf(title, "MB(X0) prof Ph in region %d", i);
00294 me500[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00295 sprintf(name, "%d", i+600);
00296 sprintf(title, "MB(L0) prof Ph in region %d", i);
00297 me600[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00298 sprintf(name, "%d", i+700);
00299 sprintf(title, "MB(Step) prof Ph in region %d", i);
00300 me700[i] = tfile->make<TProfile>(name, title, binPhi, -maxPhi, maxPhi);
00301 sprintf(name, "%d", i+800);
00302 sprintf(title, "Phi in region %d", i);
00303 me800[i] = tfile->make<TH1F>(name, title, binPhi, -maxPhi, maxPhi);
00304 sprintf(name, "%d", i+900);
00305 sprintf(title, "MB(X0) prof Eta Phi in region %d", i);
00306 me900[i] = tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
00307 binPhi/2, -maxPhi, maxPhi);
00308 sprintf(name, "%d", i+1000);
00309 sprintf(title, "MB(L0) prof Eta Phi in region %d", i);
00310 me1000[i]= tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
00311 binPhi/2, -maxPhi, maxPhi);
00312 sprintf(name, "%d", i+1100);
00313 sprintf(title, "MB(Step) prof Eta Phi in region %d", i);
00314 me1100[i]= tfile->make<TProfile2D>(name, title, binEta/2, -maxEta, maxEta,
00315 binPhi/2, -maxPhi, maxPhi);
00316 sprintf(name, "%d", i+1200);
00317 sprintf(title, "Eta vs Phi in region %d", i);
00318 me1200[i]= tfile->make<TH2F>(name, title, binEta/2, -maxEta, maxEta,
00319 binPhi/2, -maxPhi, maxPhi);
00320 }
00321 for (int i=0; i<maxSet2; i++) {
00322 sprintf(name, "%d", i+1300);
00323 sprintf(title, "Events with layers Hit (0 all, 1 HB, ..) for %d", i);
00324 me1300[i]= tfile->make<TH1F>(name, title, binEta, -maxEta, maxEta);
00325 sprintf(name, "%d", i+1400);
00326 sprintf(title, "Eta vs Phi for layers hit in %d", i);
00327 me1400[i]= tfile->make<TH2F>(name, title, binEta/2, -maxEta, maxEta,
00328 binPhi/2, -maxPhi, maxPhi);
00329 sprintf(name, "%d", i+1500);
00330 sprintf(title, "Number of layers crossed (0 all, 1 HB, ..) for %d", i);
00331 me1500[i]= tfile->make<TProfile>(name, title, binEta, -maxEta, maxEta);
00332 }
00333
00334 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Booking user "
00335 << "histos done ===";
00336
00337 }
00338
00339 void MaterialBudgetHcalHistos::fillHisto(int ii) {
00340
00341 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos:FillHisto called "
00342 << "with index " << ii << " integrated step "
00343 << stepLen << " X0 " << radLen << " Lamda "
00344 << intLen;
00345
00346 if (ii >=0 && ii < maxSet) {
00347 me100[ii]->Fill(eta, radLen);
00348 me200[ii]->Fill(eta, intLen);
00349 me300[ii]->Fill(eta, stepLen);
00350 me400[ii]->Fill(eta);
00351
00352 if (eta >= etaLow && eta <= etaHigh) {
00353 me500[ii]->Fill(phi, radLen);
00354 me600[ii]->Fill(phi, intLen);
00355 me700[ii]->Fill(phi, stepLen);
00356 me800[ii]->Fill(phi);
00357 }
00358
00359 me900[ii]->Fill(eta, phi, radLen);
00360 me1000[ii]->Fill(eta, phi, intLen);
00361 me1100[ii]->Fill(eta, phi, stepLen);
00362 me1200[ii]->Fill(eta, phi);
00363
00364 }
00365 }
00366
00367 void MaterialBudgetHcalHistos::fillLayer() {
00368
00369 me1300[0]->Fill(eta);
00370 me1400[0]->Fill(eta,phi);
00371 if (nlayHB > 0) {
00372 me1300[1]->Fill(eta);
00373 me1400[1]->Fill(eta,phi);
00374 }
00375 if (nlayHB >= 16) {
00376 me1300[2]->Fill(eta);
00377 me1400[2]->Fill(eta,phi);
00378 }
00379 if (nlayHE > 0) {
00380 me1300[3]->Fill(eta);
00381 me1400[3]->Fill(eta,phi);
00382 }
00383 if (nlayHE >= 16) {
00384 me1300[4]->Fill(eta);
00385 me1400[4]->Fill(eta,phi);
00386 }
00387 if (nlayHO > 0) {
00388 me1300[5]->Fill(eta);
00389 me1400[5]->Fill(eta,phi);
00390 }
00391 if (nlayHO >= 2) {
00392 me1300[6]->Fill(eta);
00393 me1400[6]->Fill(eta,phi);
00394 }
00395 if (nlayHF > 0) {
00396 me1300[7]->Fill(eta);
00397 me1400[7]->Fill(eta,phi);
00398 }
00399 if (nlayHB > 0 || nlayHE > 0 || (nlayHF > 0 && std::abs(eta) > 3.0)) {
00400 me1300[8]->Fill(eta);
00401 me1400[8]->Fill(eta,phi);
00402 }
00403 me1500[0]->Fill(eta,(double)(nlayHB+nlayHO+nlayHE+nlayHF));
00404 me1500[1]->Fill(eta,(double)(nlayHB));
00405 me1500[2]->Fill(eta,(double)(nlayHE));
00406 me1500[4]->Fill(eta,(double)(nlayHF));
00407 }
00408
00409 void MaterialBudgetHcalHistos::hend() {
00410 edm::LogInfo("MaterialBudget") << "MaterialBudgetHcalHistos: Save user "
00411 << "histos ===";
00412 }
00413
00414 std::vector<std::string> MaterialBudgetHcalHistos::getNames(DDFilteredView& fv) {
00415
00416 std::vector<std::string> tmp;
00417 bool dodet = fv.firstChild();
00418 while (dodet) {
00419 const DDLogicalPart & log = fv.logicalPart();
00420 std::string namx = log.name().name();
00421 bool ok = true;
00422 for (unsigned int i=0; i<tmp.size(); i++)
00423 if (namx == tmp[i]) ok = false;
00424 if (ok) tmp.push_back(namx);
00425 dodet = fv.next();
00426 }
00427 return tmp;
00428 }
00429
00430 std::vector<double> MaterialBudgetHcalHistos::getDDDArray(const std::string & str,
00431 const DDsvalues_type & sv) {
00432
00433 LogDebug("MaterialBudget") << "MaterialBudgetHcalHistos:getDDDArray called "
00434 << "for " << str;
00435 DDValue value(str);
00436 if (DDfetch(&sv,value)) {
00437 LogDebug("MaterialBudget") << value;
00438 const std::vector<double> & fvec = value.doubles();
00439 int nval = fvec.size();
00440 if (nval < 1) {
00441 edm::LogError("MaterialBudget") << "MaterialBudgetHcalHistos : # of "
00442 << str << " bins " << nval
00443 << " < 1 ==> illegal";
00444 throw cms::Exception("Unknown", "MaterialBudgetHcalHistos")
00445 << "nval < 1 for array " << str << "\n";
00446 }
00447
00448 return fvec;
00449 } else {
00450 edm::LogError("MaterialBudget") << "MaterialBudgetHcalHistos : cannot get "
00451 << "array " << str;
00452 throw cms::Exception("Unknown", "MaterialBudgetHcalHistos")
00453 << "cannot get array " << str << "\n";
00454 }
00455 }
00456
00457 bool MaterialBudgetHcalHistos::isSensitive (std::string name) {
00458
00459 std::vector<std::string>::const_iterator it = sensitives.begin();
00460 std::vector<std::string>::const_iterator itEnd = sensitives.end();
00461 for (; it != itEnd; ++it)
00462 if (name == *it) return true;
00463 return false;
00464 }
00465
00466 bool MaterialBudgetHcalHistos::isItHF (const G4VTouchable* touch) {
00467
00468
00469 int levels = ((touch->GetHistoryDepth())+1);
00470 for (unsigned int it=0; it < hfNames.size(); it++) {
00471 if (levels >= hfLevels[it]) {
00472 std::string name = touch->GetVolume(levels-hfLevels[it])->GetName();
00473 if (name == hfNames[it]) return true;
00474 }
00475 }
00476 return false;
00477 }
00478
00479 bool MaterialBudgetHcalHistos::isItEC (std::string name) {
00480
00481 std::vector<std::string>::const_iterator it = sensitiveEC.begin();
00482 std::vector<std::string>::const_iterator itEnd = sensitiveEC.end();
00483 for (; it != itEnd; ++it)
00484 if (name == *it) return true;
00485 return false;
00486 }
00487
1.7.3