CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoLocalCalo/CaloTowersCreator/src/CaloTowersCreationAlgo.cc

Go to the documentation of this file.
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    // (for momentum reconstruction algorithm)
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                                                // (momentum reconstruction algorithm)
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     // (momentum reconstruction algorithm)
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                                                // (for the momentum construction algorithm)
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     // (momentum reconstruction algorithm)
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 // this method should not be used any more as the towers in the changed format
00281 // can not be properly rescaled with the "rescale" method.
00282 // "rescale was replaced by "rescaleTowers"
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   // now copy this map into the final collection
00295   for(MetaTowerMap::const_iterator mapItr = theTowerMap.begin();
00296       mapItr != theTowerMap.end(); ++mapItr) {
00297 
00298     // Convert only if there is at least one constituent in the metatower. 
00299     // The check of constituents size in the coverted tower is still needed!
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(); // save the memory
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; // not used: we do not change thresholds
00322       double weight    = 1.0;
00323 
00324       // HF
00325       if (ctcItr->ietaAbs()>=30) {
00326         double E_short = 0.5 * newE_had;             // from the definitions for HF
00327         double E_long  = newE_em + 0.5 * newE_had;   //
00328         // scale
00329         E_long  *= theHF1weight;
00330         E_short *= theHF2weight;
00331         // convert
00332         newE_em  = E_long - E_short;
00333         newE_had = 2.0 * E_short;
00334       }
00335 
00336       else {   // barrel/endcap
00337 
00338         // find if its in EB, or EE; determine from first ecal constituent found
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         // HO
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         // HB/HE
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     }   // barrel/endcap region
00367 
00368     // now make the new tower
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       // for HF we use common point for ecal, hcal shower positions regardless of the method
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     // copy the timings, have to convert back to int, 1 unit = 0.01 ns
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     } // end of loop over towers
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   // this is for skipping channls: mostly needed for the creation of
00422   // bad towers from hits i the bad channel collections.
00423   if (chStatusForCT==CaloTowersCreationAlgo::IgnoredChan) return;
00424 
00425   double threshold, weight;
00426   getThresholdAndWeight(detId, threshold, weight);
00427 
00428   double energy = recHit->energy();  // original RecHit energy is used to apply thresholds  
00429   double e = energy * weight;        // energies scaled by user weight: used in energy assignments
00430         
00431         
00432   // SPECIAL handling of tower 28/depth 3 --> half into tower 28 and half into tower 29
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     // bad channels are counted regardless of energy threshold
00449 
00450     if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
00451         tower28.numBadHcalCells += 1;
00452         tower29.numBadHcalCells += 1;
00453     }
00454 
00455     else if (0.5*energy >= threshold) {  // not bad channel: use energy if above 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       // NOTE DIVIDE BY 2!!!
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       // time info: do not use in averaging if timing error is found: need 
00480       // full set of status info to implement: use only "good" channels for now
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       // store the energy in layer 3 also in E_outer
00490       tower28.E_outer += e28;
00491       tower29.E_outer += e29;
00492     } // not a "bad" hit
00493       
00494   }  // end of special case 
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       // For ECAL we count all bad channels after the metatower is complete 
00510 
00511       // Include options for symmetric thresholds and cut on Et
00512       // for ECAL RecHits
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       //      if (chStatusForCT != CaloTowersCreationAlgo::BadChan && energy >= threshold) {
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         // change when full status info is available
00542         // for now use only good channels
00543         
00544         // add e>0 check (new options allow e<0)
00545         if (chStatusForCT == CaloTowersCreationAlgo::GoodChan && e>0 ) {
00546           tower.emSumTimeTimesE += ( e * recHit->time() );
00547           tower.emSumEForTime   += e;  // see above
00548         }
00549 
00550         std::pair<DetId,double> mc(detId,e);
00551         tower.metaConstituents.push_back(mc);
00552       }
00553     
00554     }  // end of ECAL
00555     
00556     // HCAL
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; // store HO energy even if HO is not used
00570           // add energy of the tower and/or flag if theHOIsUsed
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           } // HO is used
00581         
00582           
00583           // add HO to constituents even if it is not used: JetMET wants to keep these towers
00584           std::pair<DetId,double> mc(detId,e);
00585           tower.metaConstituents.push_back(mc);   
00586 
00587         } // not a bad channel, energy above threshold
00588       
00589       }  // HO hit 
00590       
00591       // HF calculates EM fraction differently
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             // long fiber, so E_EM = E(Long) - E(Short)
00601             tower.E_em += e;
00602           } 
00603           else {
00604             // short fiber, EHAD = 2 * E(Short)
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           // put the timing in HCAL -> have to check timing errors when available
00617           // for now use only good channels
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         } // not a bad HF channel, energy above threshold
00627       
00628       } // HF hit
00629       
00630       else {
00631         // HCAL situation normal in HB/HE
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           // Timing information: need specific accessors
00647           // for now use only good channels
00648           if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
00649             tower.hadSumTimeTimesE += ( e * recHit->time() );
00650             tower.hadSumEForTime   += e;
00651           }
00652           // store energy in highest depth for towers 18-27 (for electron,photon ID in endcap)
00653           // also, store energy in HE part of tower 16 (for JetMET cleanup)
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         }   // not a "bad" channel, energy above threshold
00666       
00667       } // channel in HBHE (excluding twrs 28,29)
00668     
00669     }
00670     
00671   }  // recHit normal case (not in HE towers 28,29)
00672 
00673 }  // end of assignHit method
00674 
00675 
00676 
00677 
00678 // This method is not flexible enough for the new CaloTower format. 
00679 // For now make a quick compatibility "fix" : WILL NOT WORK CORRECTLY with anything 
00680 // except the default simple p4 assignment!!!
00681 // Must be rewritten for full functionality.
00682 void CaloTowersCreationAlgo::rescale(const CaloTower * ct) {
00683   double threshold, weight;
00684 
00685   CaloTowerDetId towerDetId = ct->id();
00686 //  if (towerDetId.null()) return;    
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     // this is to be compliant with the new MetaTower setup
00715     // used only for the default simple vector assignment
00716     std::pair<DetId, double> mc(detId, 0);
00717     tower.metaConstituents.push_back(mc);
00718   }
00719 
00720   // preserve time inforamtion
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     // need to build a new tower
00736     MetaTower t;
00737 
00738     // store it in the map
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     // conditional assignment of depths for barrel/endcap
00756     // Some additional tuning may be required in the transitional region
00757     // 14<|iEta|<19
00758     if (id.ietaAbs()<=17) {
00759       theMomHadDepth = theMomHBDepth;
00760       theMomEmDepth  = theMomEBDepth;
00761     }
00762     else {
00763       theMomHadDepth = theMomHEDepth;
00764       theMomEmDepth  = theMomEEDepth;
00765     }
00766 
00767 
00768     // Note: E_outer is used to save HO energy OR energy in the outermost depths in endcap region
00769     // In the methods with separate treatment of EM and HAD components:
00770     //  - HO is not used to determine direction, however HO energy is added to get "total had energy"
00771     //  => Check if the tower is within HO coverage before adding E_outer to the "total had" energy
00772     //     else the energy will be double counted
00773     // When summing up the energy of the tower these checks are performed in the loops over RecHits
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) { // ignore EM threshold in HF
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; // not subtracted before, think it should be done
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     // create CaloTower using the selected algorithm
00810 
00811     GlobalPoint emPoint, hadPoint;
00812 
00813     CaloTower::PolarLorentzVector towerP4;
00814     
00815 
00816   switch (theMomConstrMethod) {
00817 
00818   case 0 :
00819     {  // Simple 4-momentum assignment
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     }  // end case 0
00828     break;
00829 
00830   case 1 :
00831     {   // separate 4-vectors for ECAL, HCAL, add to get the 4-vector of the tower (=>tower has mass!)
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 {  // forward detector: use the CaloTower position 
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);  // simple momentum assignment, same position
00851         emPoint  = p;   
00852         hadPoint = p;
00853       }
00854     }  // end case 1
00855     break;
00856 
00857   case 2:
00858     {   // use ECAL position for the tower (when E_cal>0), else default CaloTower position (massless tower)
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 {  // forward detector: use the CaloTower position 
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);  // simple momentum assignment, same position
00872         emPoint  = p;   
00873         hadPoint = p;
00874       }
00875     }   // end case 2
00876     break;
00877 
00878   }  // end of decision on p4 reconstruction method
00879 
00880 
00881     CaloTower retval(id, E_em, E_had, E_outer, -1, -1, towerP4, emPoint, hadPoint);
00882 
00883     // set the timings
00884     retval.setEcalTime(compactTime(ecalTime));
00885     retval.setHcalTime(compactTime(hcalTime));
00886 
00887     // set the CaloTower status word =====================================
00888     // Channels must be counter exclusively in the defined cathegories
00889     // "Bad" channels (not used in energy assignment) can be flagged during
00890     // CaloTower creation only if specified in the configuration file
00891 
00892     unsigned int numBadHcalChan  = mt.numBadHcalCells;
00893     //    unsigned int numBadEcalChan  = mt.numBadEcalCells;
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     // now add dead/off/... channels not used in RecHit reconstruction for HCAL 
00902     if (hcalDropChMap.find(id) != hcalDropChMap.end()) numBadHcalChan += hcalDropChMap[id];
00903     
00904 
00905     // for ECAL the number of all bad channels is obtained here -----------------------
00906 
00907     // get all possible constituents of the tower
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       // "bad" include "recovered" if the flag not to use them was set 
00925       // this is consistent with the treatment of what channels are accepted
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; // for storing the hottest cell E in the calotower
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         // need an extra check because of the funny towers that are empty except for the presence of an HO
00948         // hit in the constituents (JetMET wanted them saved)
00949         // This constituent is only used for storing the tower, but should not be concidered as a hot cell canditate for
00950         // configurations with useHO = false
00951 
00952 
00953         if (i->first.det()==DetId::Ecal) {  // ECAL
00954           maxCellE = i->second;
00955         }       
00956         else {  // HCAL
00957           if (HcalDetId(i->first).subdet() != HcalOuter) 
00958             maxCellE = i->second;
00959           else if (theHOIsUsed) maxCellE = i->second;
00960         }
00961 
00962       } // found higher E cell
00963 
00964     } // loop over matacontains
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; // in case the hit is not identified
00977 
00978   if(det == DetId::Ecal) {
00979     // may or may not be EB.  We'll find out.
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       // check if it's single or double tower
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       //check if it's ring 0 or +1 or +2 or -1 or -2
01034       if(hcalDetId.ietaAbs() <= 4) threshold = theHOthreshold0;
01035       else if(hcalDetId.ieta() < 0) {
01036         // set threshold for ring -1 or -2
01037         threshold = (hcalDetId.ietaAbs() <= 10) ?  theHOthresholdMinus1 : theHOthresholdMinus2;
01038       } else {
01039         // set threshold for ring +1 or +2
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();  // face of the cell
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();  // face of the cell
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   // this is based on available RecHits, can lead to different actual depths if
01154   // hits in multi-depth towers are not all there
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     // do not use HO for deirection calculations for now
01167     if (HcalDetId(mc_it->first).subdet() == HcalOuter) continue;
01168     ++nConst;
01169 
01170     GlobalPoint p = hadSegmentShwrPos(mc_it->first, fracDepth);
01171 
01172     // longitudinal segmentation: do not weight by energy,
01173     // get the geometrical position
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   // set depth using geometry of cells that are associated with the
01186   // tower (regardless if they have non-zero energies)
01187 
01188 //  if (hadE <= 0) return GlobalPoint(0, 0, 0);
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     // barrel, one depth only
01202     frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01203     backCellId  = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01204   }
01205   else if (towerId.ietaAbs() == 15) {
01206     // barrel, two depths
01207     frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01208     backCellId  = HcalDetId(HcalBarrel, iEta, iPhi, 2);
01209   }
01210   else if (towerId.ietaAbs() == 16) {
01211     // barrel and endcap: two depths HB, one depth HE 
01212     frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
01213     backCellId  = HcalDetId(HcalEndcap, iEta, iPhi, 3);  // this cell is in endcap!
01214   }
01215   else if (towerId.ietaAbs() == 17) {
01216     // endcap, one depth only
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   // endcap: two depths
01222   frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
01223   backCellId  = HcalDetId(HcalEndcap, iEta, iPhi, 2);
01224   }
01225   else if (towerId.ietaAbs() <= 29) {
01226   // endcap: three depths
01227   frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
01228   // there is no iEta=29 for depth 3
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   // forward, take the goemetry for long fibers
01235   frontCellId = HcalDetId(HcalForward, iEta, iPhi, 1);
01236   backCellId  = HcalDetId(HcalForward, iEta, iPhi, 1);
01237   }
01238   else {
01239     // should not get here
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    // uses the "front" and "back" cells
01251    // to determine the axis. point set by the predefined depth.
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;  // add crystals with E/E_EM > 1.5%
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; // discretization (ns)
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 // Bad/anomolous cell handling 
01355 
01356 
01357 
01358 
01359 void CaloTowersCreationAlgo::makeHcalDropChMap() {
01360 
01361   // This method fills the map of number of dead channels for the calotower,
01362   // The key of the map is CaloTowerDetId.
01363   // By definition these channels are not going to be in the RecHit collections.
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       // special case for tower 29: if HCAL hit is in depth 3 add to twr 29 as well
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   // For use with hits rejected in the default reconstruction
01409   if (useRejectedHitsOnly) {
01410     
01411     if (!isRecovered) {
01412       if (severityLevel <= int(theHcalAcceptSeverityLevel) ||
01413           severityLevel > int(theHcalAcceptSeverityLevelForRejectedHit)) return CaloTowersCreationAlgo::IgnoredChan;
01414       // this hit was either already accepted or is worse than 
01415     }    
01416     else {
01417       
01418       if (theRecoveredHcalHitsAreUsed || !useRejectedRecoveredHcalHits) {
01419         // skip recovered hits either because they were already used or because there was an explicit instruction
01420         return CaloTowersCreationAlgo::IgnoredChan;
01421       }
01422       else if (useRejectedRecoveredHcalHits) {
01423         return CaloTowersCreationAlgo::RecoveredChan;
01424       }  
01425   
01426     }  // recovered channels
01427 
01428     // clasify channels as problematic: no good hits are supposed to be present in the
01429     // extra rechit collections
01430     return CaloTowersCreationAlgo::ProblematicChan;
01431 
01432   }  // treatment of rejected hits
01433 
01434 
01435 
01436 
01437   // this is for the regular reconstruction sequence
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   //  uint16_t dbStatus = theEcalChStatus->find(id)->getStatusCode();
01463   //  uint32_t rhFlags  = hit->flags();
01464   //  int severityLevel = theEcalSevLvlAlgo->severityLevel(rhFlags, dbStatus);
01465   // The methods above will become private and cannot be usef for flagging ecal spikes.
01466   // Use the recommended interface - we leave the parameters for spilke removal to be specified by ECAL.
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   // there should be no other ECAL types used in this reconstruction
01475 
01476   // The definition of ECAL severity levels uses categories that
01477   // are similar to the defined for CaloTower. (However, the categorization
01478   // for CaloTowers depends on the specified maximum acceptabel severity and therefore cannnot
01479   // be exact correspondence between the two. ECAL has additional categories describing modes of failure.)
01480   // This approach is different from the initial idea and from
01481   // the implementation for HCAL. Still make the logic similar to HCAL so that one has the ability to 
01482   // exclude problematic channels as defined by ECAL.
01483   // For definitions of ECAL severity levels see RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h
01484 
01485 
01486   bool isRecovered = (severityLevel == EcalSeverityLevelAlgo::kRecovered);
01487 
01488 
01489   // for use with rejected by the regular reconstruction hits
01490   if (useRejectedHitsOnly) {
01491     
01492     if (!isRecovered) {
01493       if (severityLevel <= int(theEcalAcceptSeverityLevel) ||
01494           severityLevel > int(theEcalAcceptSeverityLevelForRejectedHit)) return CaloTowersCreationAlgo::IgnoredChan;
01495       // this hit was either accepted already, or counted as bad under the conditions
01496     }    
01497     else {
01498       
01499       if (theRecoveredEcalHitsAreUsed || !useRejectedRecoveredEcalHits) {
01500         // skip recovered hits either because they were already used or because there was an explicit instruction
01501         return CaloTowersCreationAlgo::IgnoredChan;
01502       }
01503       else if (useRejectedRecoveredEcalHits) {
01504         return CaloTowersCreationAlgo::RecoveredChan;
01505       }  
01506   
01507     }  // recovered channels
01508 
01509     // clasify channels as problematic
01510     return CaloTowersCreationAlgo::ProblematicChan;
01511 
01512   }  // treatment of rejected hits
01513 
01514 
01515 
01516 
01517 
01518   // for normal reconstruction
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