00001
00002
00003
00005 #include "Validation/HcalHits/interface/SimG4HcalValidation.h"
00006
00007 #include "SimG4Core/Notification/interface/BeginOfJob.h"
00008 #include "SimG4Core/Notification/interface/BeginOfRun.h"
00009 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
00010 #include "SimG4Core/Notification/interface/EndOfEvent.h"
00011
00012
00013 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00014 #include "SimG4CMS/Calo/interface/HCalSD.h"
00015 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
00016 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
00017 #include "DataFormats/Math/interface/Point3D.h"
00018
00019 #include "FWCore/Framework/interface/Event.h"
00020 #include "FWCore/Framework/interface/EventSetup.h"
00021 #include "FWCore/Framework/interface/ESTransientHandle.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023
00024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00025 #include "DetectorDescription/Core/interface/DDCompactView.h"
00026
00027 #include "G4SDManager.hh"
00028 #include "G4Step.hh"
00029 #include "G4Track.hh"
00030 #include "G4VProcess.hh"
00031 #include "G4HCofThisEvent.hh"
00032 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00033 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00034 #include <cmath>
00035 #include <iostream>
00036 #include <iomanip>
00037
00038 SimG4HcalValidation::SimG4HcalValidation(const edm::ParameterSet &p):
00039 jetf(0), numberingFromDDD(0), org(0) {
00040
00041 edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("SimG4HcalValidation");
00042 infolevel = m_Anal.getParameter<int>("InfoLevel");
00043 hcalOnly = m_Anal.getParameter<bool>("HcalClusterOnly");
00044 applySampling = m_Anal.getParameter<bool>("HcalSampling");
00045 coneSize = m_Anal.getParameter<double>("ConeSize");
00046 ehitThreshold = m_Anal.getParameter<double>("EcalHitThreshold");
00047 hhitThreshold = m_Anal.getParameter<double>("HcalHitThreshold");
00048 timeLowlim = m_Anal.getParameter<double>("TimeLowLimit");
00049 timeUplim = m_Anal.getParameter<double>("TimeUpLimit");
00050 jetThreshold = m_Anal.getParameter<double>("JetThreshold");
00051 eta0 = m_Anal.getParameter<double>("Eta0");
00052 phi0 = m_Anal.getParameter<double>("Phi0");
00053 names = m_Anal.getParameter<std::vector<std::string> >("Names");
00054 labelLayer = m_Anal.getUntrackedParameter<std::string>("LabelLayerInfo","HcalInfoLayer");
00055 labelNxN = m_Anal.getUntrackedParameter<std::string>("LabelNxNInfo","HcalInfoNxN");
00056 labelJets = m_Anal.getUntrackedParameter<std::string>("LabelJetsInfo","HcalInfoJets");
00057
00058 produces<PHcalValidInfoLayer>(labelLayer);
00059 if (infolevel > 0) produces<PHcalValidInfoNxN>(labelNxN);
00060 if (infolevel > 1) produces<PHcalValidInfoJets>(labelJets);
00061
00062 edm::LogInfo("ValidHcal") << "HcalTestAnalysis:: Initialised as observer of "
00063 << "begin/end events and of G4step with Parameter "
00064 << "values: \n\tInfoLevel = " << infolevel
00065 << "\n\thcalOnly = " << hcalOnly
00066 << "\n\tapplySampling = " << applySampling
00067 << "\n\tconeSize = " << coneSize
00068 << "\n\tehitThreshold = " << ehitThreshold
00069 << "\n\thhitThreshold = " << hhitThreshold
00070 << "\n\tttimeLowlim = " << timeLowlim
00071 << "\n\tttimeUplim = " << timeUplim
00072 << "\n\teta0 = " << eta0
00073 << "\n\tphi0 = " << phi0
00074 << "\nLabels (Layer): " << labelLayer
00075 << " (NxN): " << labelNxN << " (Jets): "
00076 << labelJets;
00077
00078 init();
00079 }
00080
00081 SimG4HcalValidation::~SimG4HcalValidation() {
00082
00083 edm::LogInfo("ValidHcal") << "\n --------> Total number of selected entries"
00084 << " : " << count << "\nPointers:: JettFinder "
00085 << jetf << ", Numbering Scheme " << org
00086 << " and FromDDD " << numberingFromDDD;
00087 if (jetf) {
00088 edm::LogInfo("ValidHcal") << "Delete Jetfinder";
00089 delete jetf;
00090 jetf = 0;
00091 }
00092 if (numberingFromDDD) {
00093 edm::LogInfo("ValidHcal") << "Delete HcalNumberingFromDDD";
00094 delete numberingFromDDD;
00095 numberingFromDDD = 0;
00096 }
00097 }
00098
00099 void SimG4HcalValidation::produce(edm::Event& e, const edm::EventSetup&) {
00100
00101 std::auto_ptr<PHcalValidInfoLayer> productLayer(new PHcalValidInfoLayer);
00102 layerAnalysis(*productLayer);
00103 e.put(productLayer,labelLayer);
00104
00105 if (infolevel > 0) {
00106 std::auto_ptr<PHcalValidInfoNxN> productNxN(new PHcalValidInfoNxN);
00107 nxNAnalysis(*productNxN);
00108 e.put(productNxN,labelNxN);
00109 }
00110
00111 if (infolevel > 1) {
00112 std::auto_ptr<PHcalValidInfoJets> productJets(new PHcalValidInfoJets);
00113 jetAnalysis(*productJets);
00114 e.put(productJets,labelJets);
00115 }
00116 }
00117
00118
00119 void SimG4HcalValidation::init() {
00120
00121 float sHB[4] = { 117., 117., 178., 217.};
00122 float sHE[3] = { 178., 178., 178.};
00123 float sHF[3] = { 2.84, 2.09, 0.};
00124
00125 float deta[4] = { 0.0435, 0.1305, 0.2175, 0.3045 };
00126 float dphi[4] = { 0.0436, 0.1309, 0.2182, 0.3054 };
00127
00128 int i=0;
00129 for (i = 0; i < 4; i++) {
00130 scaleHB.push_back(sHB[i]);
00131 }
00132 for (i = 0; i < 3; i++) {
00133 scaleHE.push_back(sHE[i]);
00134 }
00135 for (i = 0; i < 3; i++) {
00136 scaleHF.push_back(sHF[i]);
00137 }
00138
00139
00140 for(i = 0; i < 4; i++) {
00141 dEta.push_back(deta[i]);
00142 dPhi.push_back(dphi[i]);
00143 }
00144
00145
00146 jetf = new SimG4HcalHitJetFinder(coneSize);
00147
00148
00149 count = 0;
00150
00151 }
00152
00153 void SimG4HcalValidation::update(const BeginOfJob * job) {
00154
00155
00156 edm::ESTransientHandle<DDCompactView> pDD;
00157 (*job)()->get<IdealGeometryRecord>().get(pDD);
00158 edm::LogInfo("ValidHcal") << "HcalTestAnalysis:: Initialise "
00159 << "HcalNumberingFromDDD for " << names[0];
00160 numberingFromDDD = new HcalNumberingFromDDD(names[0], (*pDD));
00161
00162
00163 org = new HcalTestNumberingScheme(false);
00164
00165 }
00166
00167 void SimG4HcalValidation::update(const BeginOfRun * run) {
00168
00169 int irun = (*run)()->GetRunID();
00170
00171 edm::LogInfo("ValidHcal") <<" =====> Begin of Run = " << irun;
00172
00173 std::string sdname = names[0];
00174 G4SDManager* sd = G4SDManager::GetSDMpointerIfExist();
00175 if (sd != 0) {
00176 G4VSensitiveDetector* aSD = sd->FindSensitiveDetector(sdname);
00177 if (aSD==0) {
00178 edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: No SD"
00179 << " with name " << sdname << " in this "
00180 << "Setup";
00181 } else {
00182 HCalSD* theCaloSD = dynamic_cast<HCalSD*>(aSD);
00183 edm::LogInfo("ValidHcal") << "SimG4HcalValidation::beginOfRun: Finds SD"
00184 << "with name " << theCaloSD->GetName()
00185 << " in this Setup";
00186 if (org) {
00187 theCaloSD->setNumberingScheme(org);
00188 edm::LogInfo("ValidHcal") << "SimG4HcalValidation::beginOfRun: set a "
00189 << "new numbering scheme";
00190 }
00191 }
00192 } else {
00193 edm::LogWarning("ValidHcal") << "SimG4HcalValidation::beginOfRun: Could "
00194 << "not get SD Manager!";
00195 }
00196
00197 }
00198
00199
00200 void SimG4HcalValidation::update(const BeginOfEvent * evt) {
00201
00202 int i = 0;
00203 edepEB = edepEE = edepHB = edepHE = edepHO = 0.;
00204 for (i = 0; i < 5; i++) edepd[i] = 0.;
00205 for (i = 0; i < 20; i++) edepl[i] = 0.;
00206 vhitec = vhithc = enEcal = enHcal = 0;
00207
00208 clear();
00209
00210 int iev = (*evt)()->GetEventID();
00211 LogDebug("ValidHcal") << "SimG4HcalValidation: =====> Begin of event = "
00212 << iev;
00213 }
00214
00215
00216 void SimG4HcalValidation::update(const G4Step * aStep) {
00217
00218 if (aStep != NULL) {
00219 G4VPhysicalVolume* curPV = aStep->GetPreStepPoint()->GetPhysicalVolume();
00220 G4String name = curPV->GetName();
00221 name.assign(name,0,3);
00222 double edeposit = aStep->GetTotalEnergyDeposit();
00223 int layer=-1, depth=-1;
00224 if (name == "EBR") {
00225 depth = 0;
00226 edepEB += edeposit;
00227 } else if (name == "EFR") {
00228 depth = 0;
00229 edepEE += edeposit;
00230 } else if (name == "HBS") {
00231 layer = (curPV->GetCopyNo()/10)%100;
00232 depth = (curPV->GetCopyNo())%10 + 1;
00233 if (depth > 0 && depth < 4 && layer >= 0 && layer < 17) {
00234 edepHB += edeposit;
00235 } else {
00236 edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error "
00237 << curPV->GetName() << curPV->GetCopyNo();
00238 depth = -1; layer = -1;
00239 }
00240 } else if (name == "HES") {
00241 layer = (curPV->GetCopyNo()/10)%100;
00242 depth = (curPV->GetCopyNo())%10 + 1;
00243 if (depth > 0 && depth < 3 && layer >= 0 && layer < 19) {
00244 edepHE += edeposit;
00245 } else {
00246 edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error "
00247 << curPV->GetName() << curPV->GetCopyNo();
00248 depth = -1; layer = -1;
00249 }
00250 } else if (name == "HTS") {
00251 layer = (curPV->GetCopyNo()/10)%100;
00252 depth = (curPV->GetCopyNo())%10 + 1;
00253 if (depth > 3 && depth < 5 && layer >= 17 && layer < 20) {
00254 edepHO += edeposit;
00255 } else {
00256 edm::LogWarning("ValidHcal") << "SimG4HcalValidation:Error "
00257 << curPV->GetName() <<curPV->GetCopyNo();
00258 depth = -1; layer = -1;
00259 }
00260 }
00261 if (depth >= 0 && depth < 5) edepd[depth] += edeposit;
00262 if (layer >= 0 && layer < 20) edepl[layer] += edeposit;
00263
00264 if (layer >= 0 && layer < 20) {
00265 LogDebug("ValidHcal") << "SimG4HcalValidation:: G4Step: " << name
00266 << " Layer " << std::setw(3) << layer << " Depth "
00267 << std::setw(2) << depth << " Edep " <<std::setw(6)
00268 << edeposit/MeV << " MeV";
00269 }
00270 }
00271 }
00272
00273
00274 void SimG4HcalValidation::update(const EndOfEvent * evt) {
00275
00276 count++;
00277
00278
00279 fill(evt);
00280 LogDebug("ValidHcal") << "SimG4HcalValidation:: --- after Fill ";
00281 }
00282
00283
00284 void SimG4HcalValidation::fill(const EndOfEvent * evt) {
00285
00286 LogDebug("ValidHcal") << "SimG4HcalValidation:Fill event "
00287 << (*evt)()->GetEventID();
00288
00289
00290 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
00291
00292 int nhc = 0, j = 0;
00293
00294
00295 int HCHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[0]);
00296 CaloG4HitCollection* theHCHC = (CaloG4HitCollection*) allHC->GetHC(HCHCid);
00297 LogDebug("ValidHcal") << "SimG4HcalValidation :: Hit Collection for "
00298 << names[0] << " of ID " << HCHCid <<" is obtained at "
00299 << theHCHC;
00300 if (HCHCid >= 0 && theHCHC > 0) {
00301 for (j = 0; j < theHCHC->entries(); j++) {
00302
00303 CaloG4Hit* aHit = (*theHCHC)[j];
00304
00305 double e = aHit->getEnergyDeposit()/GeV;
00306 double time = aHit->getTimeSlice();
00307
00308 math::XYZPoint pos = aHit->getPosition();
00309 double theta = pos.theta();
00310 double eta = -log(tan(theta * 0.5));
00311 double phi = pos.phi();
00312
00313 uint32_t unitID = aHit->getUnitID();
00314 int subdet, zside, layer, etaIndex, phiIndex, lay;
00315 org->unpackHcalIndex(unitID,subdet,zside,layer,etaIndex,phiIndex,lay);
00316
00317
00318 layer--;
00319 std::string det = "HB";
00320 if (subdet == static_cast<int>(HcalForward)) {
00321 det = "HF";
00322 uint16_t depth = aHit->getDepth();
00323 if (depth != 0) { layer += 2; lay += 2; }
00324 if (layer == 1) vhithc += e;
00325 else if (layer == 0) vhitec += e;
00326 else
00327 edm::LogInfo("ValidHcal") << "SimG4HcalValidation::HitPMT "
00328 << subdet << " " << (2*zside-1)*etaIndex
00329 << " " << phiIndex << " " << layer+1
00330 << " R " << pos.rho() << " Phi " << phi/deg
00331 << " Edep " << e << " Time " << time;
00332 } else if (subdet == static_cast<int>(HcalEndcap)) {
00333 if (etaIndex <= 20) {
00334 det = "HES";
00335 } else {
00336 det = "HED";
00337 }
00338 }
00339 LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Hcal "
00340 << det << " layer " << std::setw(2) << layer
00341 << " lay " << std::setw(2) << lay << " time "
00342 << std::setw(6) << time << " theta "
00343 << std::setw(8) << theta << " eta " << std::setw(8)
00344 << eta << " phi " << std::setw(8) << phi << " e "
00345 << std::setw(8) << e << " ID 0x" << std::hex
00346 << unitID << " ID dec " << std::dec << (int)unitID;
00347
00348
00349 if (applySampling) e *= getHcalScale(det,layer);
00350
00351
00352 if (time >= timeLowlim && time <= timeUplim && e > hhitThreshold ) {
00353 enHcal += e;
00354 CaloHit ahit(subdet,lay,e,eta,phi,time,unitID);
00355 hitcache.push_back(ahit);
00356 ++nhc;
00357 }
00358 }
00359 }
00360 LogDebug("ValidHcal") << "SimG4HcalValidation:: HCAL hits : " << nhc;
00361
00362 if (!hcalOnly) {
00363 int ndets = names.size();
00364 if (ndets > 3) ndets = 3;
00365 for (int idty = 1; idty < ndets; idty++) {
00366 std::string det = "EB";
00367 if (idty == 2) det = "EE";
00368 else if (idty > 2) det = "ES";
00369
00370 int nec = 0;
00371 int ECHCid = G4SDManager::GetSDMpointer()->GetCollectionID(names[idty]);
00372 CaloG4HitCollection* theECHC =(CaloG4HitCollection*)allHC->GetHC(ECHCid);
00373 LogDebug("ValidHcal") << "SimG4HcalValidation:: Hit Collection for "
00374 << names[idty] << " of ID " << ECHCid
00375 << " is obtained at " << theECHC;
00376 if (ECHCid >= 0 && theECHC > 0) {
00377 for (j = 0; j < theECHC->entries(); j++) {
00378
00379 CaloG4Hit* aHit = (*theECHC)[j];
00380
00381 double e = aHit->getEnergyDeposit()/GeV;
00382 double time = aHit->getTimeSlice();
00383
00384 math::XYZPoint pos = aHit->getPosition();
00385 double theta = pos.theta();
00386 double eta = -log(tan(theta/2.));
00387 double phi = pos.phi();
00388 HcalNumberingFromDDD::HcalID id = numberingFromDDD->unitID(eta,phi,1,1);
00389 uint32_t unitID = org->getUnitID(id);
00390 int subdet, zside, layer, ieta, iphi, lay;
00391 org->unpackHcalIndex(unitID,subdet,zside,layer,ieta,iphi,lay);
00392 subdet = idty+9;
00393 layer = 0;
00394 unitID = org->packHcalIndex(subdet,zside,layer,ieta,iphi,lay);
00395
00396
00397 if (time >= timeLowlim && time <= timeUplim && e > ehitThreshold ) {
00398 enEcal += e;
00399 CaloHit ahit(subdet,lay,e,eta,phi,time,unitID);
00400 hitcache.push_back(ahit);
00401 ++nec;
00402 }
00403
00404 LogDebug("ValidHcal") << "SimG4HcalValidation::debugFill Ecal "
00405 << det << " layer " << std::setw(2) << layer
00406 << " lay " << std::setw(2) << lay << " time "
00407 << std::setw(6) << time << " theta "
00408 << std::setw(8) << theta << " eta "
00409 << std::setw(8) << eta << " phi "
00410 << std::setw(8) << phi << " e " << std::setw(8)
00411 << e << " ID 0x" << std::hex << unitID
00412 << " ID dec " << std::dec << (int)unitID;
00413
00414 }
00415 }
00416
00417 LogDebug("ValidHcal") << "SimG4HcalValidation:: " << det << " hits : "
00418 << nec;
00419 }
00420 }
00421
00422 }
00423
00424 void SimG4HcalValidation::layerAnalysis(PHcalValidInfoLayer& product) {
00425
00426 int i = 0;
00427 LogDebug("ValidHcal") << "\n ===>>> SimG4HcalValidation: Energy deposit "
00428 << "in MeV\n at EB : " << std::setw(6) << edepEB/MeV
00429 << "\n at EE : " << std::setw(6) << edepEE/MeV
00430 << "\n at HB : " << std::setw(6) << edepHB/MeV
00431 << "\n at HE : " << std::setw(6) << edepHE/MeV
00432 << "\n at HO : " << std::setw(6) << edepHO/MeV
00433 << "\n ---- SimG4HcalValidation: Energy deposit in";
00434 for (i = 0; i < 5; i++)
00435 LogDebug("ValidHcal") << " Depth " << std::setw(2) << i << " E "
00436 << std::setw(8) << edepd[i]/MeV << " MeV";
00437 LogDebug("ValidHcal") << " ---- SimG4HcalValidation: Energy deposit in"
00438 << "layers";
00439 for (i = 0; i < 20; i++)
00440 LogDebug("ValidHcal") << " Layer " << std::setw(2) << i << " E "
00441 << std::setw(8) << edepl[i]/MeV << " MeV";
00442
00443 product.fillLayers(edepl, edepd, edepHO, edepHB+edepHE, edepEB+edepEE);
00444
00445
00446 product.fillHF(vhitec, vhithc, enEcal, enHcal);
00447 LogDebug("ValidHcal") << "SimG4HcalValidation::HF hits " << vhitec
00448 << " in EC and " << vhithc << " in HC\n"
00449 << " HB/HE " << enEcal
00450 << " in EC and " << enHcal << " in HC";
00451
00452
00453 fetchHits(product);
00454
00455 LogDebug("ValidHcal") << " layerAnalysis ===> after fetchHits";
00456
00457 }
00458
00459
00460 void SimG4HcalValidation::nxNAnalysis(PHcalValidInfoNxN& product) {
00461
00462 std::vector<CaloHit> * hits = &hitcache;
00463 std::vector<CaloHit>::iterator hit_itr;
00464
00465 LogDebug("ValidHcal") << "SimG4HcalValidation::NxNAnalysis : entrance ";
00466
00467 collectEnergyRdir(eta0,phi0);
00468
00469
00470 LogDebug("ValidHcal") << " NxNAnalysis : coolectEnergyRdir - Ecal " << een
00471 << " Hcal " << hen;
00472
00473 double etot = 0.;
00474 double ee = 0.;
00475 double he = 0.;
00476 double hoe = 0.;
00477
00478 int max = dEta.size();
00479
00480 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
00481
00482 double e = hit_itr->e();
00483 double t = hit_itr->t();
00484 double eta = hit_itr->eta();
00485 double phi = hit_itr->phi();
00486 int type = hit_itr->det();
00487 int layer = hit_itr->layer();
00488
00489
00490
00491 if (fabs(eta0-eta) <= dEta[max-1] && fabs(phi0-phi) <= dPhi[max-1]) {
00492 etot += e;
00493 if (type == 10 || type == 11 || type == 12) ee += e;
00494 if (type == static_cast<int>(HcalBarrel) ||
00495 type == static_cast<int>(HcalEndcap) ||
00496 type == static_cast<int>(HcalForward)) {
00497 he += e;
00498 if (type == static_cast<int>(HcalBarrel) && layer > 17) hoe += e;
00499
00500
00501 for (int i = 0; i < max; i++ ) {
00502 if ((eta0-eta) <= dEta[i] && (eta0-eta) >= -dEta[i] &&
00503 (phi0-phi) <= dPhi[i] && (phi0-phi) >= -dPhi[i]) {
00504
00505 LogDebug("ValidHcal") << "SimG4HcalValidation:: nxNAnalysis eta0,"
00506 << " phi0 = " << eta0 << " " << phi0
00507 << " type, layer = " << type << ","
00508 << layer << " eta, phi = " << eta << " "
00509 << phi;
00510
00511 product.fillTProfileNxN(e, i, t);
00512 break;
00513 }
00514 }
00515
00516 }
00517 }
00518 }
00519
00520 product.fillEcollectNxN(ee, he, hoe, etot);
00521 product.fillHvsE(een, hen, hoen, een+hen);
00522
00523 LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
00524
00525
00526 }
00527
00528
00529
00530 void SimG4HcalValidation::jetAnalysis(PHcalValidInfoJets& product) {
00531
00532 std::vector<CaloHit> * hhit = &hitcache;
00533
00534 jetf->setInput(hhit);
00535
00536
00537 std::vector<SimG4HcalHitCluster> * result = jetf->getClusters(hcalOnly);
00538
00539 std::vector<SimG4HcalHitCluster>::iterator clus_itr;
00540
00541 LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size()
00542 << " clusters ---------------";
00543 for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
00544 LogDebug("ValidHcal") << (*clus_itr);
00545
00546 std::vector<double> enevec, etavec, phivec;
00547
00548 if ((*result).size() > 0) {
00549
00550 sort((*result).begin(),(*result).end());
00551
00552 clus_itr = result->begin();
00553 double etac = clus_itr->eta();
00554 double phic = clus_itr->phi();
00555
00556 double ecal_collect = 0.;
00557 if (!hcalOnly) {
00558 ecal_collect = clus_itr->collectEcalEnergyR();}
00559 else {
00560 collectEnergyRdir(etac, phic);
00561 ecal_collect = een;
00562 }
00563 LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect = "
00564 << ecal_collect;
00565
00566
00567 double dist = jetf->rDist(eta0, phi0, etac, phic);
00568 LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation = " << dist;
00569 product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
00570
00571 LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
00572
00573 std::vector<CaloHit> * hits = clus_itr->getHits() ;
00574 std::vector<CaloHit>::iterator hit_itr;
00575
00576 double ee = 0., he = 0., hoe = 0., etot = 0.;
00577
00578
00579 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
00580 double e = hit_itr->e();
00581 double t = hit_itr->t();
00582 double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
00583
00584
00585 etot += e;
00586 if (hit_itr->det() == 10 || hit_itr->det() == 11 ||
00587 hit_itr->det() == 12) ee += e;
00588 if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
00589 hit_itr->det() == static_cast<int>(HcalEndcap) ||
00590 hit_itr->det() == static_cast<int>(HcalForward)) {
00591 he += e;
00592 if (hit_itr->det() == static_cast<int>(HcalBarrel) &&
00593 hit_itr->layer() > 17)
00594 hoe += e;
00595 }
00596
00597 if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
00598 hit_itr->det() == static_cast<int>(HcalEndcap) ||
00599 hit_itr->det() == static_cast<int>(HcalForward)) {
00600 product.fillTProfileJet(he, r, t);
00601 }
00602 }
00603
00604 product.fillEcollectJet(ee, he, hoe, etot);
00605
00606 LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
00607 << "ee/he/hoe/etot " << ee << "/" << he << "/"
00608 << hoe << "/" << etot;
00609
00610
00611 for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
00612 if (clus_itr->e() > jetThreshold) {
00613 enevec.push_back(clus_itr->e());
00614 etavec.push_back(clus_itr->eta());
00615 phivec.push_back(clus_itr->phi());
00616 }
00617 }
00618 product.fillJets(enevec, etavec, phivec);
00619
00620 LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
00621 << " JetAnalysis ===> (*result).size() "
00622 << (*result).size();
00623
00624
00625 if (etavec.size() > 1) {
00626 if (etavec[0] > -2.5 && etavec[0] < 2.5 &&
00627 etavec[1] > -2.5 && etavec[1] < 2.5) {
00628
00629 LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
00630 << " JetAnalysis ===> Di-jet vectors ";
00631 for (unsigned int i = 0; i < enevec.size(); i++)
00632 LogDebug("ValidHcal") << " e, eta, phi = " << enevec[i] << " "
00633 << etavec[i] << " " << phivec[i];
00634
00635 double et0 = enevec[0] / cosh(etavec[0]);
00636 double px0 = et0 * cos(phivec[0]);
00637 double py0 = et0 * sin(phivec[0]);
00638 double pz0 = et0 * sinh(etavec[0]);
00639 double et1 = enevec[1] / cosh(etavec[1]);
00640 double px1 = et1 * cos(phivec[1]);
00641 double py1 = et1 * sin(phivec[1]);
00642 double pz1 = et1 * sinh(etavec[1]);
00643
00644 double dijetmass2 = (enevec[0]+enevec[1])*(enevec[0]+enevec[1])
00645 - (px1+px0)*(px1+px0) - (py1+py0)*(py1+py0) - (pz1+pz0)*(pz1+pz0);
00646
00647 LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ "
00648 << dijetmass2;
00649
00650 double dijetmass;
00651 if(dijetmass2 >= 0.) dijetmass = sqrt(dijetmass2);
00652 else dijetmass = -sqrt(-1. * dijetmass2);
00653
00654 product.fillDiJets(dijetmass);
00655
00656 LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
00657 }
00658 }
00659 }
00660 }
00661
00662
00663 void SimG4HcalValidation::fetchHits(PHcalValidInfoLayer& product) {
00664
00665 LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with "
00666 << hitcache.size() << " hits";
00667 int nHit = hitcache.size();
00668 int hit = 0;
00669 int i;
00670 std::vector<CaloHit>::iterator itr;
00671 std::vector<CaloHit*> lhits(nHit);
00672 for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
00673 uint32_t unitID=itr->id();
00674 int subdet, zside, group, ieta, iphi, lay;
00675 HcalTestNumbering::unpackHcalIndex(unitID,subdet,zside,group,
00676 ieta,iphi,lay);
00677 subdet = itr->det();
00678 lay = itr->layer();
00679 group = (subdet&15)<<20;
00680 group += ((lay-1)&31)<<15;
00681 group += (zside&1)<<14;
00682 group += (ieta&127)<<7;
00683 group += (iphi&127);
00684 itr->setId(group);
00685 lhits[i] = &hitcache[i];
00686 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i
00687 << " " << hitcache[i] << "\n"
00688 << "SimG4HcalValidation::fetchHits:Copied " << i
00689 << " " << *lhits[i];
00690 }
00691 sort(lhits.begin(),lhits.end(),CaloHitIdMore());
00692 std::vector<CaloHit*>::iterator k1, k2;
00693 for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
00694 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i
00695 << " " << **k1;
00696 int nHits = 0;
00697 for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
00698 double ehit = (**k1).e();
00699 double t = (**k1).t();
00700 uint32_t unitID= (**k1).id();
00701 int jump = 0;
00702 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i
00703 << " U/T/E" << " 0x" << std::hex << unitID
00704 << std::dec << " " << t << " " << ehit;
00705 for (k2 = k1+1; k2 != lhits.end() && (t-(**k2).t()) < 1 &&
00706 (t-(**k2).t()) > -1 && unitID == (**k2).id(); k2++) {
00707 ehit += (**k2).e();
00708 LogDebug("ValidHcal") << "\t + " << (**k2).e();
00709 jump++;
00710 }
00711 LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
00712
00713 double eta = (*k1)->eta();
00714 double phi = (*k1)->phi();
00715 int lay = ((unitID>>15)&31) + 1;
00716 int subdet = (unitID>>20)&15;
00717 int zside = (unitID>>14)&1;
00718 int ieta = (unitID>>7)&127;
00719 int iphi = (unitID)&127;
00720
00721
00722 product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
00723 nHits++;
00724
00725 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits
00726 << " " << i << " ID 0x" << std::hex << unitID
00727 << " det " << std::dec << subdet << " " << lay
00728 << " " << zside << " " << ieta << " " << iphi
00729 << " Time " << t << " E " << ehit;
00730
00731 i += jump;
00732 k1 += jump;
00733 }
00734
00735 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with "
00736 << nHit << " hits" << " and writes out " << nHits
00737 << '(' << hit << ") hits";
00738
00739 }
00740
00741 void SimG4HcalValidation::clear(){
00742 hitcache.erase( hitcache.begin(), hitcache.end());
00743 }
00744
00745
00746 void SimG4HcalValidation::collectEnergyRdir (const double eta0,
00747 const double phi0) {
00748
00749 std::vector<CaloHit> * hits = &hitcache;
00750 std::vector<CaloHit>::iterator hit_itr;
00751
00752 double sume = 0., sumh = 0., sumho = 0.;
00753
00754 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
00755
00756 double e = hit_itr->e();
00757 double eta = hit_itr->eta();
00758 double phi = hit_itr->phi();
00759
00760 int type = hit_itr->det();
00761
00762 double r = jetf->rDist(eta0, phi0, eta, phi);
00763 if (r < coneSize) {
00764 if (type == 10 || type == 11 || type == 12) {
00765 sume += e;
00766 } else {
00767 sumh += e;
00768 if (type == static_cast<int>(HcalBarrel) &&
00769 hit_itr->layer() > 17) sumho += e;
00770 }
00771 }
00772 }
00773
00774 een = sume;
00775 hen = sumh;
00776 hoen = sumho;
00777 }
00778
00779
00780 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
00781
00782 double tmp = 0.;
00783
00784 if (det == "HB") {
00785 tmp = scaleHB[layer];
00786 } else if (det == "HES" || det == "HED") {
00787 tmp = scaleHE[layer];
00788 } else if (det == "HF") {
00789 tmp = scaleHF[layer];
00790 }
00791
00792 return tmp;
00793 }
1.7.3