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 int save_i = 0;
00492 if (fabs(eta0-eta) <= dEta[max-1] && fabs(phi0-phi) <= dPhi[max-1]) {
00493 etot += e;
00494 if (type == 10 || type == 11 || type == 12) ee += e;
00495 if (type == static_cast<int>(HcalBarrel) ||
00496 type == static_cast<int>(HcalEndcap) ||
00497 type == static_cast<int>(HcalForward)) {
00498 he += e;
00499 if (type == static_cast<int>(HcalBarrel) && layer > 17) hoe += e;
00500
00501
00502 for (int i = 0; i < max; i++ ) {
00503 if ((eta0-eta) <= dEta[i] && (eta0-eta) >= -dEta[i] &&
00504 (phi0-phi) <= dPhi[i] && (phi0-phi) >= -dPhi[i]) {
00505
00506 LogDebug("ValidHcal") << "SimG4HcalValidation:: nxNAnalysis eta0,"
00507 << " phi0 = " << eta0 << " " << phi0
00508 << " type, layer = " << type << ","
00509 << layer << " eta, phi = " << eta << " "
00510 << phi;
00511
00512 product.fillTProfileNxN(e, i, t);
00513 save_i = i;
00514 break;
00515 }
00516 }
00517
00518 }
00519 }
00520 }
00521
00522 product.fillEcollectNxN(ee, he, hoe, etot);
00523 product.fillHvsE(een, hen, hoen, een+hen);
00524
00525 LogDebug("ValidHcal") << " nxNAnalysis ===> after fillHvsE";
00526
00527
00528 }
00529
00530
00531
00532 void SimG4HcalValidation::jetAnalysis(PHcalValidInfoJets& product) {
00533
00534 std::vector<CaloHit> * hhit = &hitcache;
00535
00536 jetf->setInput(hhit);
00537
00538
00539 std::vector<SimG4HcalHitCluster> * result = jetf->getClusters(hcalOnly);
00540
00541 std::vector<SimG4HcalHitCluster>::iterator clus_itr;
00542
00543 LogDebug("ValidHcal") << "\n ---------- Final list of " << (*result).size()
00544 << " clusters ---------------";
00545 for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++)
00546 LogDebug("ValidHcal") << (*clus_itr);
00547
00548 std::vector<double> enevec, etavec, phivec;
00549
00550 if ((*result).size() > 0) {
00551
00552 sort((*result).begin(),(*result).end());
00553
00554 clus_itr = result->begin();
00555 double etac = clus_itr->eta();
00556 double phic = clus_itr->phi();
00557
00558 double ecal_collect = 0.;
00559 if (!hcalOnly) {
00560 ecal_collect = clus_itr->collectEcalEnergyR();}
00561 else {
00562 collectEnergyRdir(etac, phic);
00563 ecal_collect = een;
00564 }
00565 LogDebug("ValidHcal") << " JetAnalysis ===> ecal_collect = "
00566 << ecal_collect;
00567
00568
00569 double dist = jetf->rDist(eta0, phi0, etac, phic);
00570 LogDebug("ValidHcal") << " JetAnalysis ===> eta phi deviation = " << dist;
00571 product.fillEtaPhiProfileJet(eta0, phi0, etac, phic, dist);
00572
00573 LogDebug("ValidHcal") << " JetAnalysis ===> after fillEtaPhiProfileJet";
00574
00575 std::vector<CaloHit> * hits = clus_itr->getHits() ;
00576 std::vector<CaloHit>::iterator hit_itr;
00577
00578 double ee = 0., he = 0., hoe = 0., etot = 0.;
00579
00580
00581 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
00582 double e = hit_itr->e();
00583 double t = hit_itr->t();
00584 double r = jetf->rDist(&(*clus_itr), &(*hit_itr));
00585
00586
00587 etot += e;
00588 if (hit_itr->det() == 10 || hit_itr->det() == 11 ||
00589 hit_itr->det() == 12) ee += e;
00590 if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
00591 hit_itr->det() == static_cast<int>(HcalEndcap) ||
00592 hit_itr->det() == static_cast<int>(HcalForward)) {
00593 he += e;
00594 if (hit_itr->det() == static_cast<int>(HcalBarrel) &&
00595 hit_itr->layer() > 17)
00596 hoe += e;
00597 }
00598
00599 if (hit_itr->det() == static_cast<int>(HcalBarrel) ||
00600 hit_itr->det() == static_cast<int>(HcalEndcap) ||
00601 hit_itr->det() == static_cast<int>(HcalForward)) {
00602 product.fillTProfileJet(he, r, t);
00603 }
00604 }
00605
00606 product.fillEcollectJet(ee, he, hoe, etot);
00607
00608 LogDebug("ValidHcal") << " JetAnalysis ===> after fillEcollectJet: "
00609 << "ee/he/hoe/etot " << ee << "/" << he << "/"
00610 << hoe << "/" << etot;
00611
00612
00613 for (clus_itr = result->begin(); clus_itr < result->end(); clus_itr++) {
00614 if (clus_itr->e() > jetThreshold) {
00615 enevec.push_back(clus_itr->e());
00616 etavec.push_back(clus_itr->eta());
00617 phivec.push_back(clus_itr->phi());
00618 }
00619 }
00620 product.fillJets(enevec, etavec, phivec);
00621
00622 LogDebug("ValidHcal") << " JetAnalysis ===> after fillJets\n"
00623 << " JetAnalysis ===> (*result).size() "
00624 << (*result).size();
00625
00626
00627 if (etavec.size() > 1) {
00628 if (etavec[0] > -2.5 && etavec[0] < 2.5 &&
00629 etavec[1] > -2.5 && etavec[1] < 2.5) {
00630
00631 LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet mass enter\n"
00632 << " JetAnalysis ===> Di-jet vectors ";
00633 for (unsigned int i = 0; i < enevec.size(); i++)
00634 LogDebug("ValidHcal") << " e, eta, phi = " << enevec[i] << " "
00635 << etavec[i] << " " << phivec[i];
00636
00637 double et0 = enevec[0] / cosh(etavec[0]);
00638 double px0 = et0 * cos(phivec[0]);
00639 double py0 = et0 * sin(phivec[0]);
00640 double pz0 = et0 * sinh(etavec[0]);
00641 double et1 = enevec[1] / cosh(etavec[1]);
00642 double px1 = et1 * cos(phivec[1]);
00643 double py1 = et1 * sin(phivec[1]);
00644 double pz1 = et1 * sinh(etavec[1]);
00645
00646 double dijetmass2 = (enevec[0]+enevec[1])*(enevec[0]+enevec[1])
00647 - (px1+px0)*(px1+px0) - (py1+py0)*(py1+py0) - (pz1+pz0)*(pz1+pz0);
00648
00649 LogDebug("ValidHcal") << " JetAnalysis ===> Di-jet massSQ "
00650 << dijetmass2;
00651
00652 double dijetmass;
00653 if(dijetmass2 >= 0.) dijetmass = sqrt(dijetmass2);
00654 else dijetmass = -sqrt(-1. * dijetmass2);
00655
00656 product.fillDiJets(dijetmass);
00657
00658 LogDebug("ValidHcal") << " JetAnalysis ===> after fillDiJets";
00659 }
00660 }
00661 }
00662 }
00663
00664
00665 void SimG4HcalValidation::fetchHits(PHcalValidInfoLayer& product) {
00666
00667 LogDebug("ValidHcal") << "Enter SimG4HcalValidation::fetchHits with "
00668 << hitcache.size() << " hits";
00669 int nHit = hitcache.size();
00670 int hit = 0;
00671 int i;
00672 std::vector<CaloHit>::iterator itr;
00673 std::vector<CaloHit*> lhits(nHit);
00674 for (i = 0, itr = hitcache.begin(); itr != hitcache.end(); i++, itr++) {
00675 uint32_t unitID=itr->id();
00676 int subdet, zside, group, ieta, iphi, lay;
00677 HcalTestNumbering::unpackHcalIndex(unitID,subdet,zside,group,
00678 ieta,iphi,lay);
00679 subdet = itr->det();
00680 lay = itr->layer();
00681 group = (subdet&15)<<20;
00682 group += ((lay-1)&31)<<15;
00683 group += (zside&1)<<14;
00684 group += (ieta&127)<<7;
00685 group += (iphi&127);
00686 itr->setId(group);
00687 lhits[i] = &hitcache[i];
00688 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Original " << i
00689 << " " << hitcache[i] << "\n"
00690 << "SimG4HcalValidation::fetchHits:Copied " << i
00691 << " " << *lhits[i];
00692 }
00693 sort(lhits.begin(),lhits.end(),CaloHitIdMore());
00694 std::vector<CaloHit*>::iterator k1, k2;
00695 for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++)
00696 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Sorted " << i
00697 << " " << **k1;
00698 int nHits = 0;
00699 for (i = 0, k1 = lhits.begin(); k1 != lhits.end(); i++, k1++) {
00700 double ehit = (**k1).e();
00701 double t = (**k1).t();
00702 uint32_t unitID= (**k1).id();
00703 int jump = 0;
00704 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Start " << i
00705 << " U/T/E" << " 0x" << std::hex << unitID
00706 << std::dec << " " << t << " " << ehit;
00707 for (k2 = k1+1; k2 != lhits.end() && (t-(**k2).t()) < 1 &&
00708 (t-(**k2).t()) > -1 && unitID == (**k2).id(); k2++) {
00709 ehit += (**k2).e();
00710 LogDebug("ValidHcal") << "\t + " << (**k2).e();
00711 jump++;
00712 }
00713 LogDebug("ValidHcal") << "\t = " << ehit << " in " << jump;
00714
00715 double eta = (*k1)->eta();
00716 double phi = (*k1)->phi();
00717 int lay = ((unitID>>15)&31) + 1;
00718 int subdet = (unitID>>20)&15;
00719 int zside = (unitID>>14)&1;
00720 int ieta = (unitID>>7)&127;
00721 int iphi = (unitID)&127;
00722
00723
00724 product.fillHits(nHits, lay, subdet, eta, phi, ehit, t);
00725 nHits++;
00726
00727 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits:Hit " << nHits
00728 << " " << i << " ID 0x" << std::hex << unitID
00729 << " det " << std::dec << subdet << " " << lay
00730 << " " << zside << " " << ieta << " " << iphi
00731 << " Time " << t << " E " << ehit;
00732
00733 i += jump;
00734 k1 += jump;
00735 }
00736
00737 LogDebug("ValidHcal") << "SimG4HcalValidation::fetchHits called with "
00738 << nHit << " hits" << " and writes out " << nHits
00739 << '(' << hit << ") hits";
00740
00741 }
00742
00743 void SimG4HcalValidation::clear(){
00744 hitcache.erase( hitcache.begin(), hitcache.end());
00745 }
00746
00747
00748 void SimG4HcalValidation::collectEnergyRdir (const double eta0,
00749 const double phi0) {
00750
00751 std::vector<CaloHit> * hits = &hitcache;
00752 std::vector<CaloHit>::iterator hit_itr;
00753
00754 double sume = 0., sumh = 0., sumho = 0.;
00755
00756 for (hit_itr = hits->begin(); hit_itr < hits->end(); hit_itr++) {
00757
00758 double e = hit_itr->e();
00759 double eta = hit_itr->eta();
00760 double phi = hit_itr->phi();
00761
00762 int type = hit_itr->det();
00763
00764 double r = jetf->rDist(eta0, phi0, eta, phi);
00765 if (r < coneSize) {
00766 if (type == 10 || type == 11 || type == 12) {
00767 sume += e;
00768 } else {
00769 sumh += e;
00770 if (type == static_cast<int>(HcalBarrel) &&
00771 hit_itr->layer() > 17) sumho += e;
00772 }
00773 }
00774 }
00775
00776 een = sume;
00777 hen = sumh;
00778 hoen = sumho;
00779 }
00780
00781
00782 double SimG4HcalValidation::getHcalScale(std::string det, int layer) const {
00783
00784 double tmp = 0.;
00785
00786 if (det == "HB") {
00787 tmp = scaleHB[layer];
00788 } else if (det == "HES" || det == "HED") {
00789 tmp = scaleHE[layer];
00790 } else if (det == "HF") {
00791 tmp = scaleHF[layer];
00792 }
00793
00794 return tmp;
00795 }
1.7.3