00001 #include "RecoLocalCalo/CaloTowersCreator/interface/CaloTowersCreationAlgo.h"
00002 #include "Geometry/CaloTopology/interface/HcalTopology.h"
00003 #include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
00004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00006 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00007 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00008 #include "RecoLocalCalo/CaloTowersCreator/interface/HcalMaterials.h"
00009 #include "Math/Interpolator.h"
00010 #include<iostream>
00011
00012 CaloTowersCreationAlgo::CaloTowersCreationAlgo()
00013 : theEBthreshold(-1000.),
00014 theEEthreshold(-1000.),
00015
00016 theUseEtEBTresholdFlag(false),
00017 theUseEtEETresholdFlag(false),
00018 theUseSymEBTresholdFlag(false),
00019 theUseSymEETresholdFlag(false),
00020
00021
00022 theHcalThreshold(-1000.),
00023 theHBthreshold(-1000.),
00024 theHESthreshold(-1000.),
00025 theHEDthreshold(-1000.),
00026 theHOthreshold0(-1000.),
00027 theHOthresholdPlus1(-1000.),
00028 theHOthresholdMinus1(-1000.),
00029 theHOthresholdPlus2(-1000.),
00030 theHOthresholdMinus2(-1000.),
00031 theHF1threshold(-1000.),
00032 theHF2threshold(-1000.),
00033 theEBGrid(std::vector<double>(5,10.)),
00034 theEBWeights(std::vector<double>(5,1.)),
00035 theEEGrid(std::vector<double>(5,10.)),
00036 theEEWeights(std::vector<double>(5,1.)),
00037 theHBGrid(std::vector<double>(5,10.)),
00038 theHBWeights(std::vector<double>(5,1.)),
00039 theHESGrid(std::vector<double>(5,10.)),
00040 theHESWeights(std::vector<double>(5,1.)),
00041 theHEDGrid(std::vector<double>(5,10.)),
00042 theHEDWeights(std::vector<double>(5,1.)),
00043 theHOGrid(std::vector<double>(5,10.)),
00044 theHOWeights(std::vector<double>(5,1.)),
00045 theHF1Grid(std::vector<double>(5,10.)),
00046 theHF1Weights(std::vector<double>(5,1.)),
00047 theHF2Grid(std::vector<double>(5,10.)),
00048 theHF2Weights(std::vector<double>(5,1.)),
00049 theEBweight(1.),
00050 theEEweight(1.),
00051 theHBweight(1.),
00052 theHESweight(1.),
00053 theHEDweight(1.),
00054 theHOweight(1.),
00055 theHF1weight(1.),
00056 theHF2weight(1.),
00057 theEcutTower(-1000.),
00058 theEBSumThreshold(-1000.),
00059 theEESumThreshold(-1000.),
00060 theHcalTopology(0),
00061 theGeometry(0),
00062 theTowerConstituentsMap(0),
00063 theHOIsUsed(true),
00064
00065 theMomConstrMethod(0),
00066 theMomHBDepth(0.),
00067 theMomHEDepth(0.),
00068 theMomEBDepth(0.),
00069 theMomEEDepth(0.)
00070 {
00071 }
00072
00073 CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold, double EEthreshold,
00074
00075 bool useEtEBTreshold,
00076 bool useEtEETreshold,
00077 bool useSymEBTreshold,
00078 bool useSymEETreshold,
00079
00080 double HcalThreshold,
00081 double HBthreshold, double HESthreshold, double HEDthreshold,
00082 double HOthreshold0, double HOthresholdPlus1, double HOthresholdMinus1,
00083 double HOthresholdPlus2, double HOthresholdMinus2,
00084 double HF1threshold, double HF2threshold,
00085 double EBweight, double EEweight,
00086 double HBweight, double HESweight, double HEDweight,
00087 double HOweight, double HF1weight, double HF2weight,
00088 double EcutTower, double EBSumThreshold, double EESumThreshold,
00089 bool useHO,
00090
00091 int momConstrMethod,
00092 double momHBDepth,
00093 double momHEDepth,
00094 double momEBDepth,
00095 double momEEDepth)
00096
00097 : theEBthreshold(EBthreshold),
00098 theEEthreshold(EEthreshold),
00099
00100 theUseEtEBTresholdFlag(useEtEBTreshold),
00101 theUseEtEETresholdFlag(useEtEETreshold),
00102 theUseSymEBTresholdFlag(useSymEBTreshold),
00103 theUseSymEETresholdFlag(useSymEETreshold),
00104
00105 theHcalThreshold(HcalThreshold),
00106 theHBthreshold(HBthreshold),
00107 theHESthreshold(HESthreshold),
00108 theHEDthreshold(HEDthreshold),
00109 theHOthreshold0(HOthreshold0),
00110 theHOthresholdPlus1(HOthresholdPlus1),
00111 theHOthresholdMinus1(HOthresholdMinus1),
00112 theHOthresholdPlus2(HOthresholdPlus2),
00113 theHOthresholdMinus2(HOthresholdMinus2),
00114 theHF1threshold(HF1threshold),
00115 theHF2threshold(HF2threshold),
00116 theEBGrid(std::vector<double>(5,10.)),
00117 theEBWeights(std::vector<double>(5,1.)),
00118 theEEGrid(std::vector<double>(5,10.)),
00119 theEEWeights(std::vector<double>(5,1.)),
00120 theHBGrid(std::vector<double>(5,10.)),
00121 theHBWeights(std::vector<double>(5,1.)),
00122 theHESGrid(std::vector<double>(5,10.)),
00123 theHESWeights(std::vector<double>(5,1.)),
00124 theHEDGrid(std::vector<double>(5,10.)),
00125 theHEDWeights(std::vector<double>(5,1.)),
00126 theHOGrid(std::vector<double>(5,10.)),
00127 theHOWeights(std::vector<double>(5,1.)),
00128 theHF1Grid(std::vector<double>(5,10.)),
00129 theHF1Weights(std::vector<double>(5,1.)),
00130 theHF2Grid(std::vector<double>(5,10.)),
00131 theHF2Weights(std::vector<double>(5,1.)),
00132 theEBweight(EBweight),
00133 theEEweight(EEweight),
00134 theHBweight(HBweight),
00135 theHESweight(HESweight),
00136 theHEDweight(HEDweight),
00137 theHOweight(HOweight),
00138 theHF1weight(HF1weight),
00139 theHF2weight(HF2weight),
00140 theEcutTower(EcutTower),
00141 theEBSumThreshold(EBSumThreshold),
00142 theEESumThreshold(EESumThreshold),
00143 theHOIsUsed(useHO),
00144
00145 theMomConstrMethod(momConstrMethod),
00146 theMomHBDepth(momHBDepth),
00147 theMomHEDepth(momHEDepth),
00148 theMomEBDepth(momEBDepth),
00149 theMomEEDepth(momEEDepth)
00150
00151 {
00152 }
00153
00154 CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold, double EEthreshold,
00155
00156 bool useEtEBTreshold,
00157 bool useEtEETreshold,
00158 bool useSymEBTreshold,
00159 bool useSymEETreshold,
00160
00161 double HcalThreshold,
00162 double HBthreshold, double HESthreshold, double HEDthreshold,
00163 double HOthreshold0, double HOthresholdPlus1, double HOthresholdMinus1,
00164 double HOthresholdPlus2, double HOthresholdMinus2,
00165 double HF1threshold, double HF2threshold,
00166 std::vector<double> EBGrid, std::vector<double> EBWeights,
00167 std::vector<double> EEGrid, std::vector<double> EEWeights,
00168 std::vector<double> HBGrid, std::vector<double> HBWeights,
00169 std::vector<double> HESGrid, std::vector<double> HESWeights,
00170 std::vector<double> HEDGrid, std::vector<double> HEDWeights,
00171 std::vector<double> HOGrid, std::vector<double> HOWeights,
00172 std::vector<double> HF1Grid, std::vector<double> HF1Weights,
00173 std::vector<double> HF2Grid, std::vector<double> HF2Weights,
00174 double EBweight, double EEweight,
00175 double HBweight, double HESweight, double HEDweight,
00176 double HOweight, double HF1weight, double HF2weight,
00177 double EcutTower, double EBSumThreshold, double EESumThreshold,
00178 bool useHO,
00179
00180 int momConstrMethod,
00181 double momHBDepth,
00182 double momHEDepth,
00183 double momEBDepth,
00184 double momEEDepth)
00185
00186 : theEBthreshold(EBthreshold),
00187 theEEthreshold(EEthreshold),
00188
00189 theUseEtEBTresholdFlag(useEtEBTreshold),
00190 theUseEtEETresholdFlag(useEtEETreshold),
00191 theUseSymEBTresholdFlag(useSymEBTreshold),
00192 theUseSymEETresholdFlag(useSymEETreshold),
00193
00194 theHcalThreshold(HcalThreshold),
00195 theHBthreshold(HBthreshold),
00196 theHESthreshold(HESthreshold),
00197 theHEDthreshold(HEDthreshold),
00198 theHOthreshold0(HOthreshold0),
00199 theHOthresholdPlus1(HOthresholdPlus1),
00200 theHOthresholdMinus1(HOthresholdMinus1),
00201 theHOthresholdPlus2(HOthresholdPlus2),
00202 theHOthresholdMinus2(HOthresholdMinus2),
00203 theHF1threshold(HF1threshold),
00204 theHF2threshold(HF2threshold),
00205 theEBGrid(EBGrid),
00206 theEBWeights(EBWeights),
00207 theEEGrid(EEGrid),
00208 theEEWeights(EEWeights),
00209 theHBGrid(HBGrid),
00210 theHBWeights(HBWeights),
00211 theHESGrid(HESGrid),
00212 theHESWeights(HESWeights),
00213 theHEDGrid(HEDGrid),
00214 theHEDWeights(HEDWeights),
00215 theHOGrid(HOGrid),
00216 theHOWeights(HOWeights),
00217 theHF1Grid(HF1Grid),
00218 theHF1Weights(HF1Weights),
00219 theHF2Grid(HF2Grid),
00220 theHF2Weights(HF2Weights),
00221 theEBweight(EBweight),
00222 theEEweight(EEweight),
00223 theHBweight(HBweight),
00224 theHESweight(HESweight),
00225 theHEDweight(HEDweight),
00226 theHOweight(HOweight),
00227 theHF1weight(HF1weight),
00228 theHF2weight(HF2weight),
00229 theEcutTower(EcutTower),
00230 theEBSumThreshold(EBSumThreshold),
00231 theEESumThreshold(EESumThreshold),
00232 theHOIsUsed(useHO),
00233
00234 theMomConstrMethod(momConstrMethod),
00235 theMomHBDepth(momHBDepth),
00236 theMomHEDepth(momHEDepth),
00237 theMomEBDepth(momEBDepth),
00238 theMomEEDepth(momEEDepth)
00239
00240 {
00241 }
00242
00243
00244 void CaloTowersCreationAlgo::setGeometry(const CaloTowerConstituentsMap* ctt, const HcalTopology* topo, const CaloGeometry* geo) {
00245 theTowerConstituentsMap=ctt;
00246 theHcalTopology = topo;
00247 theGeometry = geo;
00248 theTowerGeometry=geo->getSubdetectorGeometry(DetId::Calo,CaloTowerDetId::SubdetId);
00249 }
00250
00251 void CaloTowersCreationAlgo::begin() {
00252 theTowerMap.clear();
00253 hcalDropChMap.clear();
00254 }
00255
00256 void CaloTowersCreationAlgo::process(const HBHERecHitCollection& hbhe) {
00257 for(HBHERecHitCollection::const_iterator hbheItr = hbhe.begin();
00258 hbheItr != hbhe.end(); ++hbheItr)
00259 assignHit(&(*hbheItr));
00260 }
00261
00262 void CaloTowersCreationAlgo::process(const HORecHitCollection& ho) {
00263 for(HORecHitCollection::const_iterator hoItr = ho.begin();
00264 hoItr != ho.end(); ++hoItr)
00265 assignHit(&(*hoItr));
00266 }
00267
00268 void CaloTowersCreationAlgo::process(const HFRecHitCollection& hf) {
00269 for(HFRecHitCollection::const_iterator hfItr = hf.begin();
00270 hfItr != hf.end(); ++hfItr)
00271 assignHit(&(*hfItr));
00272 }
00273
00274 void CaloTowersCreationAlgo::process(const EcalRecHitCollection& ec) {
00275 for(EcalRecHitCollection::const_iterator ecItr = ec.begin();
00276 ecItr != ec.end(); ++ecItr)
00277 assignHit(&(*ecItr));
00278 }
00279
00280
00281
00282
00283
00284 void CaloTowersCreationAlgo::process(const CaloTowerCollection& ctc) {
00285 for(CaloTowerCollection::const_iterator ctcItr = ctc.begin();
00286 ctcItr != ctc.end(); ++ctcItr) {
00287 rescale(&(*ctcItr));
00288 }
00289 }
00290
00291
00292
00293 void CaloTowersCreationAlgo::finish(CaloTowerCollection& result) {
00294
00295 for(MetaTowerMap::const_iterator mapItr = theTowerMap.begin();
00296 mapItr != theTowerMap.end(); ++mapItr) {
00297
00298
00299
00300 if ( (mapItr->second).metaConstituents.size()<1) continue;
00301
00302 CaloTower ct=convert(mapItr->first,mapItr->second);
00303 if (ct.constituentsSize()>0 && ct.energy()>theEcutTower) {
00304 result.push_back(ct);
00305 }
00306 }
00307 theTowerMap.clear();
00308 }
00309
00310
00311 void CaloTowersCreationAlgo::rescaleTowers(const CaloTowerCollection& ctc, CaloTowerCollection& ctcResult) {
00312
00313 for (CaloTowerCollection::const_iterator ctcItr = ctc.begin();
00314 ctcItr != ctc.end(); ++ctcItr) {
00315
00316 CaloTowerDetId twrId = ctcItr->id();
00317 double newE_em = ctcItr->emEnergy();
00318 double newE_had = ctcItr->hadEnergy();
00319 double newE_outer = ctcItr->outerEnergy();
00320
00321 double threshold = 0.0;
00322 double weight = 1.0;
00323
00324
00325 if (ctcItr->ietaAbs()>=30) {
00326 double E_short = 0.5 * newE_had;
00327 double E_long = newE_em + 0.5 * newE_had;
00328
00329 E_long *= theHF1weight;
00330 E_short *= theHF2weight;
00331
00332 newE_em = E_long - E_short;
00333 newE_had = 2.0 * E_short;
00334 }
00335
00336 else {
00337
00338
00339 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
00340 DetId constId = ctcItr->constituent(iConst);
00341 if (constId.det()!=DetId::Ecal) continue;
00342 getThresholdAndWeight(constId, threshold, weight);
00343 newE_em *= weight;
00344 break;
00345 }
00346
00347 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
00348 DetId constId = ctcItr->constituent(iConst);
00349 if (constId.det()!=DetId::Hcal) continue;
00350 if (HcalDetId(constId).subdet()!=HcalOuter) continue;
00351 getThresholdAndWeight(constId, threshold, weight);
00352 newE_outer *= weight;
00353 break;
00354 }
00355
00356 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
00357 DetId constId = ctcItr->constituent(iConst);
00358 if (constId.det()!=DetId::Hcal) continue;
00359 if (HcalDetId(constId).subdet()==HcalOuter) continue;
00360 getThresholdAndWeight(constId, threshold, weight);
00361 newE_had *= weight;
00362 if (ctcItr->ietaAbs()>16) newE_outer *= weight;
00363 break;
00364 }
00365
00366 }
00367
00368
00369
00370 double newE_hadTot = (theHOIsUsed && twrId.ietaAbs()<16)? newE_had+newE_outer : newE_had;
00371
00372 GlobalPoint emPoint = ctcItr->emPosition();
00373 GlobalPoint hadPoint = ctcItr->emPosition();
00374
00375 double f_em = 1.0/cosh(emPoint.eta());
00376 double f_had = 1.0/cosh(hadPoint.eta());
00377
00378 CaloTower::PolarLorentzVector towerP4;
00379
00380 if (ctcItr->ietaAbs()<30) {
00381 if (newE_em>0) towerP4 += CaloTower::PolarLorentzVector(newE_em*f_em, emPoint.eta(), emPoint.phi(), 0);
00382 if (newE_hadTot>0) towerP4 += CaloTower::PolarLorentzVector(newE_hadTot*f_had, hadPoint.eta(), hadPoint.phi(), 0);
00383 }
00384 else {
00385 double newE_tot = newE_em + newE_had;
00386
00387 if (newE_tot>0) towerP4 += CaloTower::PolarLorentzVector(newE_tot*f_had, hadPoint.eta(), hadPoint.phi(), 0);
00388 }
00389
00390
00391
00392 CaloTower rescaledTower(twrId, newE_em, newE_had, newE_outer, -1, -1, towerP4, emPoint, hadPoint);
00393
00394 rescaledTower.setEcalTime( int(ctcItr->ecalTime()*100.0 + 0.5) );
00395 rescaledTower.setHcalTime( int(ctcItr->hcalTime()*100.0 + 0.5) );
00396
00397 std::vector<DetId> contains;
00398 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
00399 contains.push_back(ctcItr->constituent(iConst));
00400 }
00401 rescaledTower.addConstituents(contains);
00402
00403 rescaledTower.setCaloTowerStatus(ctcItr->towerStatusWord());
00404
00405 ctcResult.push_back(rescaledTower);
00406
00407 }
00408
00409
00410 }
00411
00412
00413
00414 void CaloTowersCreationAlgo::assignHit(const CaloRecHit * recHit) {
00415 DetId detId = recHit->detid();
00416
00417 unsigned int chStatusForCT = (detId.det()==DetId::Hcal)?
00418 hcalChanStatusForCaloTower(recHit) :
00419 ecalChanStatusForCaloTower(recHit);
00420
00421
00422
00423 if (chStatusForCT==CaloTowersCreationAlgo::IgnoredChan) return;
00424
00425 double threshold, weight;
00426 getThresholdAndWeight(detId, threshold, weight);
00427
00428 double energy = recHit->energy();
00429 double e = energy * weight;
00430
00431
00432
00433 if (detId.det()==DetId::Hcal &&
00434 HcalDetId(detId).subdet()==HcalEndcap &&
00435 HcalDetId(detId).depth()==3 &&
00436 HcalDetId(detId).ietaAbs()==28) {
00437
00438 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
00439 if (towerDetId.null()) return;
00440 MetaTower & tower28 = find(towerDetId);
00441 CaloTowerDetId towerDetId29 = CaloTowerDetId(towerDetId.ieta()+
00442 towerDetId.zside(),
00443 towerDetId.iphi());
00444 MetaTower & tower29 = find(towerDetId29);
00445
00447
00448
00449
00450 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
00451 tower28.numBadHcalCells += 1;
00452 tower29.numBadHcalCells += 1;
00453 }
00454
00455 else if (0.5*energy >= threshold) {
00456
00457 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
00458 tower28.numRecHcalCells += 1;
00459 tower29.numRecHcalCells += 1;
00460 }
00461 else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
00462 tower28.numProbHcalCells += 1;
00463 tower29.numProbHcalCells += 1;
00464 }
00465
00466
00467 double e28 = 0.5 * e;
00468 double e29 = 0.5 * e;
00469
00470 tower28.E_had += e28;
00471 tower28.E += e28;
00472 std::pair<DetId,double> mc(detId,e28);
00473 tower28.metaConstituents.push_back(mc);
00474
00475 tower29.E_had += e29;
00476 tower29.E += e29;
00477 tower29.metaConstituents.push_back(mc);
00478
00479
00480
00481
00482 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
00483 tower28.hadSumTimeTimesE += ( e28 * recHit->time() );
00484 tower28.hadSumEForTime += e28;
00485 tower29.hadSumTimeTimesE += ( e29 * recHit->time() );
00486 tower29.hadSumEForTime += e29;
00487 }
00488
00489
00490 tower28.E_outer += e28;
00491 tower29.E_outer += e29;
00492 }
00493
00494 }
00495
00496 else {
00497
00498 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
00499
00500 if (towerDetId.null()) return;
00501 MetaTower & tower = find(towerDetId);
00502
00503 DetId::Detector det = detId.det();
00504
00505 if (det == DetId::Ecal) {
00506
00508
00509
00510
00511
00512
00513
00514 bool passEmThreshold = false;
00515
00516 if (detId.subdetId() == EcalBarrel) {
00517 if (theUseEtEBTresholdFlag) energy /= cosh( (theGeometry->getGeometry(detId)->getPosition()).eta() ) ;
00518 if (theUseSymEBTresholdFlag) passEmThreshold = (fabs(energy) >= threshold);
00519 else passEmThreshold = (energy >= threshold);
00520
00521 }
00522 else if (detId.subdetId() == EcalEndcap) {
00523 if (theUseEtEETresholdFlag) energy /= cosh( (theGeometry->getGeometry(detId)->getPosition()).eta() ) ;
00524 if (theUseSymEETresholdFlag) passEmThreshold = (fabs(energy) >= threshold);
00525 else passEmThreshold = (energy >= threshold);
00526 }
00527
00528
00529
00530 if (chStatusForCT != CaloTowersCreationAlgo::BadChan && passEmThreshold) {
00531 tower.E_em += e;
00532 tower.E += e;
00533
00534 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
00535 tower.numRecEcalCells += 1;
00536 }
00537 else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
00538 tower.numProbEcalCells += 1;
00539 }
00540
00541
00542
00543
00544
00545 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan && e>0 ) {
00546 tower.emSumTimeTimesE += ( e * recHit->time() );
00547 tower.emSumEForTime += e;
00548 }
00549
00550 std::pair<DetId,double> mc(detId,e);
00551 tower.metaConstituents.push_back(mc);
00552 }
00553
00554 }
00555
00556
00557 else {
00558 HcalDetId hcalDetId(detId);
00559
00561
00562 if(hcalDetId.subdet() == HcalOuter) {
00563
00564 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
00565 if (theHOIsUsed) tower.numBadHcalCells += 1;
00566 }
00567
00568 else if (energy >= threshold) {
00569 tower.E_outer += e;
00570
00571 if(theHOIsUsed) {
00572 tower.E += e;
00573
00574 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
00575 tower.numRecHcalCells += 1;
00576 }
00577 else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
00578 tower.numProbHcalCells += 1;
00579 }
00580 }
00581
00582
00583
00584 std::pair<DetId,double> mc(detId,e);
00585 tower.metaConstituents.push_back(mc);
00586
00587 }
00588
00589 }
00590
00591
00592 else if(hcalDetId.subdet() == HcalForward) {
00593
00594 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
00595 tower.numBadHcalCells += 1;
00596 }
00597
00598 else if (energy >= threshold) {
00599 if (hcalDetId.depth() == 1) {
00600
00601 tower.E_em += e;
00602 }
00603 else {
00604
00605 tower.E_em -= e;
00606 tower.E_had += 2. * e;
00607 }
00608 tower.E += e;
00609 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
00610 tower.numRecHcalCells += 1;
00611 }
00612 else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
00613 tower.numProbHcalCells += 1;
00614 }
00615
00616
00617
00618 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
00619 tower.hadSumTimeTimesE += ( e * recHit->time() );
00620 tower.hadSumEForTime += e;
00621 }
00622
00623 std::pair<DetId,double> mc(detId,e);
00624 tower.metaConstituents.push_back(mc);
00625
00626 }
00627
00628 }
00629
00630 else {
00631
00632
00633 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
00634 tower.numBadHcalCells += 1;
00635 }
00636 else if (energy >= threshold) {
00637 tower.E_had += e;
00638 tower.E += e;
00639 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
00640 tower.numRecHcalCells += 1;
00641 }
00642 else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
00643 tower.numProbHcalCells += 1;
00644 }
00645
00646
00647
00648 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
00649 tower.hadSumTimeTimesE += ( e * recHit->time() );
00650 tower.hadSumEForTime += e;
00651 }
00652
00653
00654 if (HcalDetId(detId).subdet()==HcalEndcap) {
00655 if ( (HcalDetId(detId).depth()==2 && HcalDetId(detId).ietaAbs()>=18 && HcalDetId(detId).ietaAbs()<27) ||
00656 (HcalDetId(detId).depth()==3 && HcalDetId(detId).ietaAbs()==27) ||
00657 (HcalDetId(detId).depth()==3 && HcalDetId(detId).ietaAbs()==16) ) {
00658 tower.E_outer += e;
00659 }
00660 }
00661
00662 std::pair<DetId,double> mc(detId,e);
00663 tower.metaConstituents.push_back(mc);
00664
00665 }
00666
00667 }
00668
00669 }
00670
00671 }
00672
00673 }
00674
00675
00676
00677
00678
00679
00680
00681
00682 void CaloTowersCreationAlgo::rescale(const CaloTower * ct) {
00683 double threshold, weight;
00684
00685 CaloTowerDetId towerDetId = ct->id();
00686
00687 MetaTower & tower = find(towerDetId);
00688
00689 tower.E_em = 0.;
00690 tower.E_had = 0.;
00691 tower.E_outer = 0.;
00692 for (unsigned int i=0; i<ct->constituentsSize(); i++) {
00693 DetId detId = ct->constituent(i);
00694 getThresholdAndWeight(detId, threshold, weight);
00695 DetId::Detector det = detId.det();
00696 if(det == DetId::Ecal) {
00697 tower.E_em = ct->emEnergy()*weight;
00698 }
00699 else {
00700 HcalDetId hcalDetId(detId);
00701 if(hcalDetId.subdet() == HcalForward) {
00702 if (hcalDetId.depth()==1) tower.E_em = ct->emEnergy()*weight;
00703 if (hcalDetId.depth()==2) tower.E_had = ct->hadEnergy()*weight;
00704 }
00705 else if(hcalDetId.subdet() == HcalOuter) {
00706 tower.E_outer = ct->outerEnergy()*weight;
00707 }
00708 else {
00709 tower.E_had = ct->hadEnergy()*weight;
00710 }
00711 }
00712 tower.E = tower.E_had+tower.E_em+tower.E_outer;
00713
00714
00715
00716 std::pair<DetId, double> mc(detId, 0);
00717 tower.metaConstituents.push_back(mc);
00718 }
00719
00720
00721 tower.emSumTimeTimesE = ct->ecalTime();
00722 tower.hadSumTimeTimesE = ct->hcalTime();
00723 tower.emSumEForTime = 1.0;
00724 tower.hadSumEForTime = 1.0;
00725 }
00726
00727 CaloTowersCreationAlgo::MetaTower::MetaTower() :
00728 E(0),E_em(0),E_had(0),E_outer(0), emSumTimeTimesE(0), hadSumTimeTimesE(0), emSumEForTime(0), hadSumEForTime(0),
00729 numBadEcalCells(0), numRecEcalCells(0), numProbEcalCells(0), numBadHcalCells(0), numRecHcalCells(0), numProbHcalCells(0) { }
00730
00731
00732 CaloTowersCreationAlgo::MetaTower & CaloTowersCreationAlgo::find(const CaloTowerDetId & detId) {
00733 MetaTowerMap::iterator itr = theTowerMap.find(detId);
00734 if(itr == theTowerMap.end()) {
00735
00736 MetaTower t;
00737
00738
00739 theTowerMap.insert(std::pair<CaloTowerDetId, CaloTowersCreationAlgo::MetaTower>(detId, t));
00740 itr = theTowerMap.find(detId);
00741
00742
00743 }
00744 return itr->second;
00745 }
00746
00747 CaloTower CaloTowersCreationAlgo::convert(const CaloTowerDetId& id, const MetaTower& mt) {
00748
00749 double ecalThres=(id.ietaAbs()<=17)?(theEBSumThreshold):(theEESumThreshold);
00750 double E=mt.E;
00751 double E_em=mt.E_em;
00752 double E_had=mt.E_had;
00753 double E_outer=mt.E_outer;
00754
00755
00756
00757
00758 if (id.ietaAbs()<=17) {
00759 theMomHadDepth = theMomHBDepth;
00760 theMomEmDepth = theMomEBDepth;
00761 }
00762 else {
00763 theMomHadDepth = theMomHEDepth;
00764 theMomEmDepth = theMomEEDepth;
00765 }
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 float ecalTime = (mt.emSumEForTime>0)? mt.emSumTimeTimesE/mt.emSumEForTime : -9999;
00777 float hcalTime = (mt.hadSumEForTime>0)? mt.hadSumTimeTimesE/mt.hadSumEForTime : -9999;
00778
00779 std::vector<std::pair<DetId,double> > metaContains=mt.metaConstituents;
00780
00781 if (id.ietaAbs()<=29 && E_em<ecalThres) {
00782 E-=E_em;
00783 E_em=0;
00784 std::vector<std::pair<DetId,double> > metaContains_noecal;
00785
00786 for (std::vector<std::pair<DetId,double> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i)
00787 if (i->first.det()!=DetId::Ecal) metaContains_noecal.push_back(*i);
00788 metaContains.swap(metaContains_noecal);
00789 }
00790
00791 if (id.ietaAbs()<=29 && E_had<theHcalThreshold) {
00792 E-=E_had;
00793
00794 if (theHOIsUsed && id.ietaAbs()<16) E-=E_outer;
00795
00796 E_had=0;
00797 E_outer=0;
00798 std::vector<std::pair<DetId,double> > metaContains_nohcal;
00799
00800 for (std::vector<std::pair<DetId,double> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i)
00801 if (i->first.det()!=DetId::Hcal) metaContains_nohcal.push_back(*i);
00802 metaContains.swap(metaContains_nohcal);
00803 }
00804
00805
00806 double E_had_tot = (theHOIsUsed && id.ietaAbs()<16)? E_had+E_outer : E_had;
00807
00808
00809
00810
00811 GlobalPoint emPoint, hadPoint;
00812
00813 CaloTower::PolarLorentzVector towerP4;
00814
00815
00816 switch (theMomConstrMethod) {
00817
00818 case 0 :
00819 {
00820 GlobalPoint p=theTowerGeometry->getGeometry(id)->getPosition();
00821
00822 double pf=1.0/cosh(p.eta());
00823 if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);
00824
00825 emPoint = p;
00826 hadPoint = p;
00827 }
00828 break;
00829
00830 case 1 :
00831 {
00832 if (id.ietaAbs()<=29) {
00833 if (E_em>0) {
00834 emPoint = emShwrPos(metaContains, theMomEmDepth, E_em);
00835 double emPf = 1.0/cosh(emPoint.eta());
00836 towerP4 += CaloTower::PolarLorentzVector(E_em*emPf, emPoint.eta(), emPoint.phi(), 0);
00837 }
00838 if ( (E_had + E_outer) >0) {
00839 hadPoint = hadShwrPos(id, theMomHadDepth);
00840 double hadPf = 1.0/cosh(hadPoint.eta());
00841
00842 if (E_had_tot>0) {
00843 towerP4 += CaloTower::PolarLorentzVector(E_had_tot*hadPf, hadPoint.eta(), hadPoint.phi(), 0);
00844 }
00845 }
00846 }
00847 else {
00848 GlobalPoint p=theTowerGeometry->getGeometry(id)->getPosition();
00849 double pf=1.0/cosh(p.eta());
00850 if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);
00851 emPoint = p;
00852 hadPoint = p;
00853 }
00854 }
00855 break;
00856
00857 case 2:
00858 {
00859 if (id.ietaAbs()<=29) {
00860 if (E_em>0) emPoint = emShwrLogWeightPos(metaContains, theMomEmDepth, E_em);
00861 else emPoint = theTowerGeometry->getGeometry(id)->getPosition();
00862
00863 double sumPf = 1.0/cosh(emPoint.eta());
00864 if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*sumPf, emPoint.eta(), emPoint.phi(), 0);
00865
00866 hadPoint = emPoint;
00867 }
00868 else {
00869 GlobalPoint p=theTowerGeometry->getGeometry(id)->getPosition();
00870 double pf=1.0/cosh(p.eta());
00871 if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);
00872 emPoint = p;
00873 hadPoint = p;
00874 }
00875 }
00876 break;
00877
00878 }
00879
00880
00881 CaloTower retval(id, E_em, E_had, E_outer, -1, -1, towerP4, emPoint, hadPoint);
00882
00883
00884 retval.setEcalTime(compactTime(ecalTime));
00885 retval.setHcalTime(compactTime(hcalTime));
00886
00887
00888
00889
00890
00891
00892 unsigned int numBadHcalChan = mt.numBadHcalCells;
00893
00894 unsigned int numBadEcalChan = 0;
00895
00896 unsigned int numRecHcalChan = mt.numRecHcalCells;
00897 unsigned int numRecEcalChan = mt.numRecEcalCells;
00898 unsigned int numProbHcalChan = mt.numProbHcalCells;
00899 unsigned int numProbEcalChan = mt.numProbEcalCells;
00900
00901
00902 if (hcalDropChMap.find(id) != hcalDropChMap.end()) numBadHcalChan += hcalDropChMap[id];
00903
00904
00905
00906
00907
00908 std::vector<DetId> allConstituents = theTowerConstituentsMap->constituentsOf(id);
00909
00910 for (std::vector<DetId>::iterator ac_it=allConstituents.begin();
00911 ac_it!=allConstituents.end(); ++ac_it) {
00912
00913 if (ac_it->det()!=DetId::Ecal) continue;
00914
00915 int thisEcalSevLvl = -999;
00916
00917 if (ac_it->subdetId() == EcalBarrel && theEbHandle.isValid()) {
00918 thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEbHandle, *theEcalChStatus);
00919 }
00920 else if (ac_it->subdetId() == EcalEndcap && theEeHandle.isValid()) {
00921 thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEeHandle, *theEcalChStatus);
00922 }
00923
00924
00925
00926 if (thisEcalSevLvl>int(theEcalAcceptSeverityLevel) ||
00927 ( thisEcalSevLvl==EcalSeverityLevelAlgo::kRecovered && !theRecoveredEcalHitsAreUsed) ) {
00928
00929 ++numBadEcalChan;
00930 }
00931
00932 }
00933
00934
00935 retval.setCaloTowerStatus(numBadHcalChan, numBadEcalChan,
00936 numRecHcalChan, numRecEcalChan,
00937 numProbHcalChan, numProbEcalChan);
00938
00939 double maxCellE = -999.0;
00940
00941 std::vector<DetId> contains;
00942 for (std::vector<std::pair<DetId,double> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i) {
00943
00944 contains.push_back(i->first);
00945
00946 if (maxCellE < i->second) {
00947
00948
00949
00950
00951
00952
00953 if (i->first.det()==DetId::Ecal) {
00954 maxCellE = i->second;
00955 }
00956 else {
00957 if (HcalDetId(i->first).subdet() != HcalOuter)
00958 maxCellE = i->second;
00959 else if (theHOIsUsed) maxCellE = i->second;
00960 }
00961
00962 }
00963
00964 }
00965
00966 retval.addConstituents(contains);
00967 retval.setHottestCellE(maxCellE);
00968
00969 return retval;
00970 }
00971
00972
00973
00974 void CaloTowersCreationAlgo::getThresholdAndWeight(const DetId & detId, double & threshold, double & weight) const {
00975 DetId::Detector det = detId.det();
00976 weight=0;
00977
00978 if(det == DetId::Ecal) {
00979
00980
00981 EcalSubdetector subdet = (EcalSubdetector)(detId.subdetId());
00982 if(subdet == EcalBarrel) {
00983 threshold = theEBthreshold;
00984 weight = theEBweight;
00985 if (weight <= 0.) {
00986 ROOT::Math::Interpolator my(theEBGrid,theEBWeights,ROOT::Math::Interpolation::kAKIMA);
00987 weight = my.Eval(theEBEScale);
00988 }
00989 }
00990 else if(subdet == EcalEndcap) {
00991 threshold = theEEthreshold;
00992 weight = theEEweight;
00993 if (weight <= 0.) {
00994 ROOT::Math::Interpolator my(theEEGrid,theEEWeights,ROOT::Math::Interpolation::kAKIMA);
00995 weight = my.Eval(theEEEScale);
00996 }
00997 }
00998 }
00999 else if(det == DetId::Hcal) {
01000 HcalDetId hcalDetId(detId);
01001 HcalSubdetector subdet = hcalDetId.subdet();
01002
01003 if(subdet == HcalBarrel) {
01004 threshold = theHBthreshold;
01005 weight = theHBweight;
01006 if (weight <= 0.) {
01007 ROOT::Math::Interpolator my(theHBGrid,theHBWeights,ROOT::Math::Interpolation::kAKIMA);
01008 weight = my.Eval(theHBEScale);
01009 }
01010 }
01011
01012 else if(subdet == HcalEndcap) {
01013
01014 if(hcalDetId.ietaAbs() < theHcalTopology->firstHEDoublePhiRing()) {
01015 threshold = theHESthreshold;
01016 weight = theHESweight;
01017 if (weight <= 0.) {
01018 ROOT::Math::Interpolator my(theHESGrid,theHESWeights,ROOT::Math::Interpolation::kAKIMA);
01019 weight = my.Eval(theHESEScale);
01020 }
01021 }
01022 else {
01023 threshold = theHEDthreshold;
01024 weight = theHEDweight;
01025 if (weight <= 0.) {
01026 ROOT::Math::Interpolator my(theHEDGrid,theHEDWeights,ROOT::Math::Interpolation::kAKIMA);
01027 weight = my.Eval(theHEDEScale);
01028 }
01029 }
01030 }
01031
01032 else if(subdet == HcalOuter) {
01033
01034 if(hcalDetId.ietaAbs() <= 4) threshold = theHOthreshold0;
01035 else if(hcalDetId.ieta() < 0) {
01036
01037 threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdMinus1 : theHOthresholdMinus2;
01038 } else {
01039
01040 threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdPlus1 : theHOthresholdPlus2;
01041 }
01042 weight = theHOweight;
01043 if (weight <= 0.) {
01044 ROOT::Math::Interpolator my(theHOGrid,theHOWeights,ROOT::Math::Interpolation::kAKIMA);
01045 weight = my.Eval(theHOEScale);
01046 }
01047 }
01048
01049 else if(subdet == HcalForward) {
01050 if(hcalDetId.depth() == 1) {
01051 threshold = theHF1threshold;
01052 weight = theHF1weight;
01053 if (weight <= 0.) {
01054 ROOT::Math::Interpolator my(theHF1Grid,theHF1Weights,ROOT::Math::Interpolation::kAKIMA);
01055 weight = my.Eval(theHF1EScale);
01056 }
01057 } else {
01058 threshold = theHF2threshold;
01059 weight = theHF2weight;
01060 if (weight <= 0.) {
01061 ROOT::Math::Interpolator my(theHF2Grid,theHF2Weights,ROOT::Math::Interpolation::kAKIMA);
01062 weight = my.Eval(theHF2EScale);
01063 }
01064 }
01065 }
01066 }
01067 else {
01068 std::cout << "BAD CELL det " << det << std::endl;
01069 }
01070 }
01071
01072 void CaloTowersCreationAlgo::setEBEScale(double scale){
01073 if (scale>0.00001) *&theEBEScale = scale;
01074 else *&theEBEScale = 50.;
01075 }
01076
01077 void CaloTowersCreationAlgo::setEEEScale(double scale){
01078 if (scale>0.00001) *&theEEEScale = scale;
01079 else *&theEEEScale = 50.;
01080 }
01081
01082 void CaloTowersCreationAlgo::setHBEScale(double scale){
01083 if (scale>0.00001) *&theHBEScale = scale;
01084 else *&theHBEScale = 50.;
01085 }
01086
01087 void CaloTowersCreationAlgo::setHESEScale(double scale){
01088 if (scale>0.00001) *&theHESEScale = scale;
01089 else *&theHESEScale = 50.;
01090 }
01091
01092 void CaloTowersCreationAlgo::setHEDEScale(double scale){
01093 if (scale>0.00001) *&theHEDEScale = scale;
01094 else *&theHEDEScale = 50.;
01095 }
01096
01097 void CaloTowersCreationAlgo::setHOEScale(double scale){
01098 if (scale>0.00001) *&theHOEScale = scale;
01099 else *&theHOEScale = 50.;
01100 }
01101
01102 void CaloTowersCreationAlgo::setHF1EScale(double scale){
01103 if (scale>0.00001) *&theHF1EScale = scale;
01104 else *&theHF1EScale = 50.;
01105 }
01106
01107 void CaloTowersCreationAlgo::setHF2EScale(double scale){
01108 if (scale>0.00001) *&theHF2EScale = scale;
01109 else *&theHF2EScale = 50.;
01110 }
01111
01112
01113 GlobalPoint CaloTowersCreationAlgo::emCrystalShwrPos(DetId detId, float fracDepth) {
01114 const CaloCellGeometry* cellGeometry = theGeometry->getGeometry(detId);
01115 GlobalPoint point = cellGeometry->getPosition();
01116
01117 if (fracDepth<0) fracDepth=0;
01118 else if (fracDepth>1) fracDepth=1;
01119
01120 if (fracDepth>0.0) {
01121 CaloCellGeometry::CornersVec cv = cellGeometry->getCorners();
01122 GlobalPoint backPoint = GlobalPoint( 0.25*( cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x() ),
01123 0.25*( cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y() ),
01124 0.25*( cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z() ) );
01125 point += fracDepth * (backPoint-point);
01126 }
01127
01128 return point;
01129 }
01130
01131 GlobalPoint CaloTowersCreationAlgo::hadSegmentShwrPos(DetId detId, float fracDepth) {
01132 const CaloCellGeometry* cellGeometry = theGeometry->getGeometry(detId);
01133 GlobalPoint point = cellGeometry->getPosition();
01134
01135 if (fracDepth<0) fracDepth=0;
01136 else if (fracDepth>1) fracDepth=1;
01137
01138 if (fracDepth>0.0) {
01139 CaloCellGeometry::CornersVec cv = cellGeometry->getCorners();
01140 GlobalPoint backPoint = GlobalPoint( 0.25*( cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x() ),
01141 0.25*( cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y() ),
01142 0.25*( cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z() ) );
01143 point += fracDepth * (backPoint-point);
01144 }
01145
01146 return point;
01147 }
01148
01149
01150 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(std::vector<std::pair<DetId,double> >& metaContains,
01151 float fracDepth, double hadE) {
01152
01153
01154
01155 if (hadE<=0) return GlobalPoint(0,0,0);
01156
01157 double hadX = 0.0;
01158 double hadY = 0.0;
01159 double hadZ = 0.0;
01160
01161 int nConst = 0;
01162
01163 std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
01164 for (; mc_it!=metaContains.end(); ++mc_it) {
01165 if (mc_it->first.det() != DetId::Hcal) continue;
01166
01167 if (HcalDetId(mc_it->first).subdet() == HcalOuter) continue;
01168 ++nConst;
01169
01170 GlobalPoint p = hadSegmentShwrPos(mc_it->first, fracDepth);
01171
01172
01173
01174 hadX += p.x();
01175 hadY += p.y();
01176 hadZ += p.z();
01177 }
01178
01179 return GlobalPoint(hadX/nConst, hadY/nConst, hadZ/nConst);
01180 }
01181
01182
01183 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(CaloTowerDetId towerId, float fracDepth) {
01184
01185
01186
01187
01188
01189
01190 if (fracDepth < 0) fracDepth = 0;
01191 else if (fracDepth > 1) fracDepth = 1;
01192
01193 GlobalPoint point(0,0,0);
01194
01195 int iEta = towerId.ieta();
01196 int iPhi = towerId.iphi();
01197
01198 HcalDetId frontCellId, backCellId;
01199
01200 if (towerId.ietaAbs() <= 14) {
01201
01202 frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01203 backCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01204 }
01205 else if (towerId.ietaAbs() == 15) {
01206
01207 frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01208 backCellId = HcalDetId(HcalBarrel, iEta, iPhi, 2);
01209 }
01210 else if (towerId.ietaAbs() == 16) {
01211
01212 frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01213 backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 3);
01214 }
01215 else if (towerId.ietaAbs() == 17) {
01216
01217 frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
01218 backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
01219 }
01220 else if (towerId.ietaAbs() >= 18 && towerId.ietaAbs() <= 26) {
01221
01222 frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
01223 backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 2);
01224 }
01225 else if (towerId.ietaAbs() <= 29) {
01226
01227 frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
01228
01229 if (iEta == 29) iEta = 28;
01230 if (iEta == -29) iEta = -28;
01231 backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 3);
01232 }
01233 else if (towerId.ietaAbs() >= 30) {
01234
01235 frontCellId = HcalDetId(HcalForward, iEta, iPhi, 1);
01236 backCellId = HcalDetId(HcalForward, iEta, iPhi, 1);
01237 }
01238 else {
01239
01240 return point;
01241 }
01242
01243 point = hadShwPosFromCells(DetId(frontCellId), DetId(backCellId), fracDepth);
01244
01245 return point;
01246 }
01247
01248 GlobalPoint CaloTowersCreationAlgo::hadShwPosFromCells(DetId frontCellId, DetId backCellId, float fracDepth) {
01249
01250
01251
01252
01253 const CaloCellGeometry* frontCellGeometry = theGeometry->getGeometry(DetId(frontCellId));
01254 const CaloCellGeometry* backCellGeometry = theGeometry->getGeometry(DetId(backCellId));
01255
01256 GlobalPoint point = frontCellGeometry->getPosition();
01257
01258 CaloCellGeometry::CornersVec cv = backCellGeometry->getCorners();
01259
01260 GlobalPoint backPoint = GlobalPoint(0.25 * (cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x()),
01261 0.25 * (cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y()),
01262 0.25 * (cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z()));
01263
01264 point += fracDepth * (backPoint - point);
01265
01266 return point;
01267 }
01268
01269
01270 GlobalPoint CaloTowersCreationAlgo::emShwrPos(std::vector<std::pair<DetId,double> >& metaContains,
01271 float fracDepth, double emE) {
01272
01273 if (emE<=0) return GlobalPoint(0,0,0);
01274
01275 double emX = 0.0;
01276 double emY = 0.0;
01277 double emZ = 0.0;
01278
01279 double eSum = 0;
01280
01281 std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
01282 for (; mc_it!=metaContains.end(); ++mc_it) {
01283 if (mc_it->first.det() != DetId::Ecal) continue;
01284 GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
01285 double e = mc_it->second;
01286
01287 if (e>0) {
01288 emX += p.x() * e;
01289 emY += p.y() * e;
01290 emZ += p.z() * e;
01291 eSum += e;
01292 }
01293
01294 }
01295
01296 return GlobalPoint(emX/eSum, emY/eSum, emZ/eSum);
01297 }
01298
01299
01300 GlobalPoint CaloTowersCreationAlgo::emShwrLogWeightPos(std::vector<std::pair<DetId,double> >& metaContains,
01301 float fracDepth, double emE) {
01302
01303 double emX = 0.0;
01304 double emY = 0.0;
01305 double emZ = 0.0;
01306
01307 double weight = 0;
01308 double sumWeights = 0;
01309 double sumEmE = 0;
01310 double crystalThresh = 0.015 * emE;
01311
01312 std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
01313 for (; mc_it!=metaContains.end(); ++mc_it) {
01314 if (mc_it->second < 0) continue;
01315 if (mc_it->first.det() == DetId::Ecal && mc_it->second > crystalThresh) sumEmE += mc_it->second;
01316 }
01317
01318 for (mc_it = metaContains.begin(); mc_it!=metaContains.end(); ++mc_it) {
01319
01320 if (mc_it->first.det() != DetId::Ecal || mc_it->second < crystalThresh) continue;
01321
01322 GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
01323
01324 weight = 4.2 + log(mc_it->second/sumEmE);
01325 sumWeights += weight;
01326
01327 emX += p.x() * weight;
01328 emY += p.y() * weight;
01329 emZ += p.z() * weight;
01330 }
01331
01332 return GlobalPoint(emX/sumWeights, emY/sumWeights, emZ/sumWeights);
01333 }
01334
01335
01336
01337
01338
01339 int CaloTowersCreationAlgo::compactTime(float time) {
01340
01341 const float timeUnit = 0.01;
01342
01343 if (time> 300.0) return 30000;
01344 if (time< -300.0) return -30000;
01345
01346 return int(time/timeUnit + 0.5);
01347
01348 }
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359 void CaloTowersCreationAlgo::makeHcalDropChMap() {
01360
01361
01362
01363
01364
01365 std::vector<DetId> allChanInStatusCont = theHcalChStatus->getAllChannels();
01366
01367 for (std::vector<DetId>::iterator it = allChanInStatusCont.begin(); it!=allChanInStatusCont.end(); ++it) {
01368
01369 const uint32_t dbStatusFlag = theHcalChStatus->getValues(*it)->getValue();
01370
01371 if (theHcalSevLvlComputer->dropChannel(dbStatusFlag)) {
01372
01373 CaloTowerDetId twrId = theTowerConstituentsMap->towerOf(*it);
01374
01375 hcalDropChMap[twrId] +=1;
01376
01377
01378 if (HcalDetId(*it).subdet()==HcalEndcap &&
01379 HcalDetId(*it).depth()==3 &&
01380 HcalDetId(*it).ietaAbs()==28) {
01381
01382 CaloTowerDetId twrId29 = CaloTowerDetId(twrId.ieta()+twrId.zside(), twrId.iphi());
01383 hcalDropChMap[twrId29] +=1;
01384 }
01385
01386 }
01387
01388 }
01389
01390 return;
01391 }
01392
01393
01394
01396
01397 unsigned int CaloTowersCreationAlgo::hcalChanStatusForCaloTower(const CaloRecHit* hit) {
01398
01399 const DetId id = hit->detid();
01400
01401 const uint32_t recHitFlag = hit->flags();
01402 const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
01403
01404 int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
01405 bool isRecovered = theHcalSevLvlComputer->recoveredRecHit(id, recHitFlag);
01406
01407
01408
01409 if (useRejectedHitsOnly) {
01410
01411 if (!isRecovered) {
01412 if (severityLevel <= int(theHcalAcceptSeverityLevel) ||
01413 severityLevel > int(theHcalAcceptSeverityLevelForRejectedHit)) return CaloTowersCreationAlgo::IgnoredChan;
01414
01415 }
01416 else {
01417
01418 if (theRecoveredHcalHitsAreUsed || !useRejectedRecoveredHcalHits) {
01419
01420 return CaloTowersCreationAlgo::IgnoredChan;
01421 }
01422 else if (useRejectedRecoveredHcalHits) {
01423 return CaloTowersCreationAlgo::RecoveredChan;
01424 }
01425
01426 }
01427
01428
01429
01430 return CaloTowersCreationAlgo::ProblematicChan;
01431
01432 }
01433
01434
01435
01436
01437
01438
01439 if (severityLevel == 0) return CaloTowersCreationAlgo::GoodChan;
01440
01441 if (isRecovered) {
01442 return (theRecoveredHcalHitsAreUsed) ?
01443 CaloTowersCreationAlgo::RecoveredChan : CaloTowersCreationAlgo::BadChan;
01444 }
01445 else {
01446 if (severityLevel > int(theHcalAcceptSeverityLevel)) {
01447 return CaloTowersCreationAlgo::BadChan;
01448 }
01449 else {
01450 return CaloTowersCreationAlgo::ProblematicChan;
01451 }
01452 }
01453
01454 }
01455
01456
01457
01458 unsigned int CaloTowersCreationAlgo::ecalChanStatusForCaloTower(const CaloRecHit* hit) {
01459
01460 const DetId id = hit->detid();
01461
01462
01463
01464
01465
01466
01467
01468
01469 int severityLevel = 999;
01470
01471 if (id.subdetId() == EcalBarrel) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEbHandle, *theEcalChStatus);
01472 else if (id.subdetId() == EcalEndcap) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEeHandle, *theEcalChStatus);
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486 bool isRecovered = (severityLevel == EcalSeverityLevelAlgo::kRecovered);
01487
01488
01489
01490 if (useRejectedHitsOnly) {
01491
01492 if (!isRecovered) {
01493 if (severityLevel <= int(theEcalAcceptSeverityLevel) ||
01494 severityLevel > int(theEcalAcceptSeverityLevelForRejectedHit)) return CaloTowersCreationAlgo::IgnoredChan;
01495
01496 }
01497 else {
01498
01499 if (theRecoveredEcalHitsAreUsed || !useRejectedRecoveredEcalHits) {
01500
01501 return CaloTowersCreationAlgo::IgnoredChan;
01502 }
01503 else if (useRejectedRecoveredEcalHits) {
01504 return CaloTowersCreationAlgo::RecoveredChan;
01505 }
01506
01507 }
01508
01509
01510 return CaloTowersCreationAlgo::ProblematicChan;
01511
01512 }
01513
01514
01515
01516
01517
01518
01519 if (severityLevel == EcalSeverityLevelAlgo::kGood) return CaloTowersCreationAlgo::GoodChan;
01520
01521 if (isRecovered) {
01522 return (theRecoveredEcalHitsAreUsed) ?
01523 CaloTowersCreationAlgo::RecoveredChan : CaloTowersCreationAlgo::BadChan;
01524 }
01525 else {
01526 if (severityLevel > int(theEcalAcceptSeverityLevel)) {
01527 return CaloTowersCreationAlgo::BadChan;
01528 }
01529 else {
01530 return CaloTowersCreationAlgo::ProblematicChan;
01531 }
01532 }
01533
01534
01535 }
01536