00001 #include "Calibration/Tools/interface/CalibElectron.h"
00002 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00003 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00004 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00005 #include "DataFormats/EgammaCandidates/interface/Electron.h"
00006 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00007 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
00008 #include "Calibration/Tools/interface/EcalIndexingTools.h"
00009 #include <iostream>
00010
00011 using namespace calib;
00012 using namespace std;
00013
00014
00015 CalibElectron::CalibElectron() : theElectron_(0), theHits_(0), theEEHits_(0)
00016 {
00017 }
00018
00019
00020 std::vector< std::pair<int,float> > CalibElectron::getCalibModulesWeights(TString calibtype)
00021 {
00022 std::vector< std::pair<int,float> > theWeights;
00023
00024 if (calibtype == "RING")
00025 {
00026 float w_ring[EcalRingCalibrationTools::N_RING_TOTAL];
00027
00028 for (int i=0;i<EcalRingCalibrationTools::N_RING_TOTAL;++i)
00029 w_ring[i]=0.;
00030
00031 std::vector<DetId> scDetIds = theElectron_->superCluster()->getHitsByDetId();
00032
00033
00034 for(std::vector<DetId>::const_iterator idIt=scDetIds.begin(); idIt!=scDetIds.end(); ++idIt){
00035
00036 const EcalRecHit* rh=0;
00037 if ( (*idIt).subdetId() == EcalBarrel)
00038 rh = &*(theHits_->find(*idIt));
00039 else if ( (*idIt).subdetId() == EcalEndcap)
00040 rh = &*(theEEHits_->find(*idIt));
00041 if (!rh)
00042 std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
00043 w_ring[EcalRingCalibrationTools::getRingIndex(*idIt)]+=rh->energy();
00044
00045 }
00046
00047 for (int i=0;i<EcalRingCalibrationTools::N_RING_TOTAL;++i)
00048 if (w_ring[i]!=0.)
00049 theWeights.push_back(std::pair<int,float>(i,w_ring[i]));
00050
00051 }
00052
00053 else if(calibtype == "MODULE")
00054 {
00055
00056 float w_ring[EcalRingCalibrationTools::N_MODULES_BARREL];
00057
00058 for (int i=0;i<EcalRingCalibrationTools::N_MODULES_BARREL;++i)
00059 w_ring[i]=0.;
00060
00061 std::vector<DetId> scDetIds = theElectron_->superCluster()->getHitsByDetId();
00062
00063
00064 for(std::vector<DetId>::const_iterator idIt=scDetIds.begin(); idIt!=scDetIds.end(); ++idIt){
00065
00066 const EcalRecHit* rh=0;
00067 if ( (*idIt).subdetId() == EcalBarrel)
00068 rh = &*(theHits_->find(*idIt));
00069 else if ( (*idIt).subdetId() == EcalEndcap)
00070 rh = &*(theEEHits_->find(*idIt));
00071 if (!rh)
00072 std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
00073 w_ring[EcalRingCalibrationTools::getModuleIndex(*idIt)]+=rh->energy();
00074
00075 }
00076
00077 for (int i=0;i<EcalRingCalibrationTools::N_MODULES_BARREL;++i)
00078 if (w_ring[i]!=0.)
00079 theWeights.push_back(std::pair<int,float>(i,w_ring[i]));
00080
00081
00082 }
00083
00084 else if(calibtype == "ABS_SCALE")
00085 {
00086 std::cout<< "ENTERING CalibElectron, ABS SCALE mode"<<std::endl;
00087
00088 float w_ring(0.);
00089
00090 std::vector<DetId> scDetIds = theElectron_->superCluster()->getHitsByDetId();
00091
00092
00093 for(std::vector<DetId>::const_iterator idIt=scDetIds.begin(); idIt!=scDetIds.end(); ++idIt){
00094
00095 const EcalRecHit* rh=0;
00096 if ( (*idIt).subdetId() == EcalBarrel)
00097 rh = &*(theHits_->find(*idIt));
00098 else if ( (*idIt).subdetId() == EcalEndcap)
00099 rh = &*(theEEHits_->find(*idIt));
00100 if (!rh)
00101 std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
00102
00103 w_ring += rh->energy();
00104
00105 }
00106
00107 if(w_ring != 0.)
00108 theWeights.push_back(std::pair<int,float>(0,w_ring));
00109 cout << " ABS SCALE - energy sum " << w_ring << endl;
00110
00111 }
00112
00113 else if(calibtype == "ETA_ET_MODE")
00114 {
00115
00116 float w_ring[ 200 ];
00117
00118 for (int i=0; i< EcalIndexingTools::getInstance()->getNumberOfChannels(); ++i)
00119 w_ring[i]=0.;
00120
00121 std::vector<DetId> scDetIds = theElectron_->superCluster()->getHitsByDetId();
00122
00123
00124 for(std::vector<DetId>::const_iterator idIt=scDetIds.begin(); idIt!=scDetIds.end(); ++idIt){
00125
00126 const EcalRecHit* rh=0;
00127 if ( (*idIt).subdetId() == EcalBarrel)
00128 rh = &*(theHits_->find(*idIt));
00129 else if ( (*idIt).subdetId() == EcalEndcap)
00130 rh = &*(theEEHits_->find(*idIt));
00131 if (!rh)
00132 std::cout << "CalibElectron::BIG ERROR::RecHit NOT FOUND" << std::endl;
00133
00134 float eta = fabs( theElectron_->eta() );
00135 float theta = 2. * atan( exp(- eta) );
00136 float et = theElectron_->superCluster()->energy() * sin(theta);
00137
00138 int in = EcalIndexingTools::getInstance()->getProgressiveIndex(eta, et);
00139
00140 w_ring[in]+=rh->energy();
00141
00142
00143 std::cout << "CalibElectron::filling channel " << in << " with value " << theElectron_->superCluster()->energy() << std::endl;
00144 }
00145
00146 for (int i=0; i< EcalIndexingTools::getInstance()->getNumberOfChannels(); ++i){
00147 if (w_ring[i]!=0.){
00148 theWeights.push_back(std::pair<int,float>(i,w_ring[i]));
00149 cout << " ring " << i << " - energy sum " << w_ring[i] << endl;
00150 }
00151
00152 }
00153
00154 }
00155
00156 else
00157 {
00158 cout << "CalibType not yet implemented" << endl;
00159
00160 }
00161
00162 return theWeights;
00163
00164 }
00165