|
|
|
#include <Calibration/EcalCalibAlgos/interface/PhiSymmetryCalibration.h>
Definition at line 37 of file PhiSymmetryCalibration.h.
| PhiSymmetryCalibration::PhiSymmetryCalibration | ( | const edm::ParameterSet & | iConfig | ) |
Constructor.
Definition at line 30 of file PhiSymmetryCalibration.cc.
00030 : 00031 nevent(0), 00032 ecalHitsProducer_( iConfig.getParameter< std::string > ("ecalRecHitsProducer") ), 00033 barrelHits_( iConfig.getParameter< std::string > ("barrelHitCollection") ), 00034 endcapHits_( iConfig.getParameter< std::string > ("endcapHitCollection") ), 00035 eCut_barl_( iConfig.getParameter< double > ("eCut_barrel") ), 00036 eCut_endc_( iConfig.getParameter< double > ("eCut_endcap") ), 00037 eventSet_( iConfig.getParameter< int > ("eventSet") ) 00038 { 00039 00040 edm::LogInfo("Calibration") << "[PhiSymmetryCalibration] Constructor called ..."; 00041 00042 00043 }
| PhiSymmetryCalibration::~PhiSymmetryCalibration | ( | ) |
| void PhiSymmetryCalibration::analyze | ( | const edm::Event & | event, | |
| const edm::EventSetup & | setup | |||
| ) | [virtual] |
Called at each event.
Implements edm::EDAnalyzer.
Definition at line 533 of file PhiSymmetryCalibration.cc.
References funct::abs(), barrelHits_, cellArea_, cellEta_, cellPos_, ecalHitsProducer_, eCut_barl_, eCut_endc_, endcapHits_, endcapRing_, lat::endl(), eta, PV3DBase< T, PVType, FrameType >::eta(), etsum_barl_, etsum_barl_miscal_, etsum_endc_, etsum_endc_miscal_, eventSet_, EBDetId::ieta(), EBDetId::iphi(), EEDetId::ix(), EEDetId::iy(), kNMiscalBins, LogDebug, meanCellArea_, miscal_, nevent, std, and EEDetId::zside().
00535 { 00536 using namespace edm; 00537 using namespace std; 00538 00539 nevent++; 00540 00541 edm::LogInfo("Calibration") << "[PhiSymmetryCalibration] New Event --------------------------------------------------------------"; 00542 00543 if ((nevent<100 && nevent%10==0) 00544 ||(nevent<1000 && nevent%100==0) 00545 ||(nevent<10000 && nevent%100==0) 00546 ||(nevent<100000 && nevent%1000==0) 00547 ||(nevent<10000000 && nevent%1000==0)) 00548 edm::LogWarning("Calibration") << "[PhiSymmetryCalibration] Events processed: "<<nevent; 00549 00550 00551 Handle<EBRecHitCollection> barrelRecHitsHandle; 00552 Handle<EERecHitCollection> endcapRecHitsHandle; 00553 00554 event.getByLabel(ecalHitsProducer_,barrelHits_,barrelRecHitsHandle); 00555 if (!barrelRecHitsHandle.isValid()) { 00556 LogDebug("") << "PhiSymmetryCalibration: Error! can't get product!" << std::endl; 00557 } 00558 00559 event.getByLabel(ecalHitsProducer_,endcapHits_,endcapRecHitsHandle); 00560 if (!endcapRecHitsHandle.isValid()) { 00561 LogDebug("") << "PhiSymmetryCalibration: Error! can't get product!" << std::endl; 00562 } 00563 00564 //Select interesting EcalRecHits (barrel) 00565 EBRecHitCollection::const_iterator itb; 00566 for (itb=barrelRecHitsHandle->begin(); itb!=barrelRecHitsHandle->end(); itb++) { 00567 EBDetId hit = EBDetId(itb->id()); 00568 float eta = cellEta_[abs(hit.ieta())-1]; 00569 float et = itb->energy()/cosh(eta); 00570 float et_thr = eCut_barl_/cosh(eta); 00571 et_thr*=1.05; 00572 if (et > et_thr && et < et_thr+0.8) { 00573 int sign = hit.ieta()>0 ? 1 : 0; 00574 etsum_barl_[abs(hit.ieta())-1][hit.iphi()-1][sign] += et; 00575 } 00576 00577 if (eventSet_==1) { 00578 //Apply a miscalibration to all crystals and increment the 00579 //ET sum, combined for all crystals 00580 for (int imiscal=0; imiscal<kNMiscalBins; imiscal++) { 00581 if (miscal_[imiscal]*et > et_thr && miscal_[imiscal]*et < et_thr+0.8) { 00582 etsum_barl_miscal_[imiscal][abs(hit.ieta())-1] += miscal_[imiscal]*et; 00583 } 00584 } 00585 } 00586 00587 } 00588 00589 //Select interesting EcalRecHits (endcaps) 00590 EERecHitCollection::const_iterator ite; 00591 for (ite=endcapRecHitsHandle->begin(); ite!=endcapRecHitsHandle->end(); ite++) { 00592 EEDetId hit = EEDetId(ite->id()); 00593 float eta = cellPos_[hit.ix()-1][hit.iy()-1].eta(); 00594 float et = ite->energy()/cosh(eta); 00595 float et_thr = eCut_endc_/cosh(eta); 00596 et_thr*=1.05; 00597 if (et > et_thr && et < et_thr+0.8) { 00598 int sign = hit.zside()>0 ? 1 : 0; 00599 etsum_endc_[hit.ix()-1][hit.iy()-1][sign] += et; 00600 } 00601 00602 if (eventSet_==1) { 00603 //Apply a miscalibration to all crystals and increment the 00604 //ET sum, combined for all crystals 00605 for (int imiscal=0; imiscal<kNMiscalBins; imiscal++) { 00606 if (miscal_[imiscal]*et > et_thr && miscal_[imiscal]*et < et_thr+0.8) { 00607 int ring = endcapRing_[hit.ix()-1][hit.iy()-1]; 00608 etsum_endc_miscal_[imiscal][ring] += miscal_[imiscal]*et*meanCellArea_[ring]/cellArea_[hit.ix()-1][hit.iy()-1]; 00609 } 00610 } 00611 } 00612 00613 } 00614 00615 }
| void PhiSymmetryCalibration::beginJob | ( | const edm::EventSetup & | iSetup | ) | [virtual] |
Called at beginning of job.
Reimplemented from edm::EDAnalyzer.
Definition at line 57 of file PhiSymmetryCalibration.cc.
References funct::abs(), barrelCells, calib, cellArea_, cellEta_, cellPhi_, cellPos_, TestMuL1L2Filter_cff::cerr, GenMuonPlsPt100GeV_cfg::cout, DetId::Ecal, EcalBarrel, EcalEndcap, EcalCondObjectContainer< T >::end(), endcapCells, endcapRing_, lat::endl(), eta, PV3DBase< T, PVType, FrameType >::eta(), etaBoundary_, etsum_barl_, etsum_barl_miscal_, etsum_endc_, etsum_endc_miscal_, eventSet_, exception, EcalCondObjectContainer< T >::find(), edm::EventSetup::get(), CaloCellGeometry::getCorners(), CaloSubdetectorGeometry::getGeometry(), EcalCondObjectContainer< T >::getMap(), CaloCellGeometry::getPosition(), CaloGeometry::getSubdetectorGeometry(), CaloGeometry::getValidDetIds(), i, EBDetId::ieta(), EBDetId::iphi(), EEDetId::ix(), EEDetId::iy(), kBarlRings, kBarlWedges, kEndcEtaRings, kEndcWedgesX, kEndcWedgesY, kMaxEndciPhi, kNMiscalBins, kSides, meanCellArea_, miscal_, nRing_, oldCalibs_barl, oldCalibs_endc, PV3DBase< T, PVType, FrameType >::phi(), phi, phi_endc, edm::ESHandle< T >::product(), DetId::rawId(), EcalCondObjectContainer< T >::size(), EEDetId::zside(), and EBDetId::zside().
00058 { 00059 00060 edm::LogInfo("Calibration") << "[PhiSymmetryCalibration] At begin job ..."; 00061 00062 // initialize arrays 00063 00064 for (int sign=0; sign<kSides; sign++) { 00065 for (int ieta=0; ieta<kBarlRings; ieta++) { 00066 for (int iphi=0; iphi<kBarlWedges; iphi++) { 00067 etsum_barl_[ieta][iphi][sign]=0.; 00068 } 00069 } 00070 for (int ix=0; ix<kEndcWedgesX; ix++) { 00071 for (int iy=0; iy<kEndcWedgesY; iy++) { 00072 etsum_endc_[ix][iy][sign]=0.; 00073 } 00074 } 00075 } 00076 00077 00078 for (int ieta=0; ieta<kBarlRings; ieta++) cellEta_[ieta]=0.; 00079 00080 for (int ix=0; ix<kEndcWedgesX; ix++) { 00081 for (int iy=0; iy<kEndcWedgesY; iy++) { 00082 cellPos_[ix][iy] = GlobalPoint(0.,0.,0.); 00083 cellPhi_[ix][iy]=0.; 00084 cellArea_[ix][iy]=0.; 00085 endcapRing_[ix][iy]=-1; 00086 } 00087 } 00088 00089 for (int imiscal=0; imiscal<kNMiscalBins; imiscal++) { 00090 miscal_[imiscal]=.95+float(imiscal)/200.; 00091 for (int ieta=0; ieta<kBarlRings; ieta++) etsum_barl_miscal_[imiscal][ieta]=0.; 00092 for (int ring=0; ring<kEndcEtaRings; ring++) etsum_endc_miscal_[imiscal][ring]=0.; 00093 } 00094 00095 // get initial constants out of DB 00096 EcalIntercalibConstantMap imap; 00097 if (eventSet_==0) { 00098 edm::ESHandle<EcalIntercalibConstants> pIcal; 00099 try { 00100 iSetup.get<EcalIntercalibConstantsRcd>().get(pIcal); 00101 std::cout << "Taken EcalIntercalibConstants" << std::endl; 00102 imap = pIcal.product()->getMap(); 00103 std::cout << "imap.size() = " << imap.size() << std::endl; 00104 } catch ( std::exception& ex ) { 00105 std::cerr << "Error! can't get EcalIntercalibConstants " << std::endl; 00106 } 00107 } 00108 00109 // get the ecal geometry: 00110 edm::ESHandle<CaloGeometry> geoHandle; 00111 iSetup.get<CaloGeometryRecord>().get(geoHandle); 00112 const CaloGeometry& geometry = *geoHandle; 00113 const CaloSubdetectorGeometry *barrelGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel); 00114 const CaloSubdetectorGeometry *endcapGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);; 00115 00116 00117 // loop over all barrel crystals 00118 barrelCells = geometry.getValidDetIds(DetId::Ecal, EcalBarrel); 00119 std::vector<DetId>::const_iterator barrelIt; 00120 for (barrelIt=barrelCells.begin(); barrelIt!=barrelCells.end(); barrelIt++) { 00121 EBDetId eb(*barrelIt); 00122 00123 if (eventSet_==0) { 00124 // get the initial calibration constants 00125 EcalIntercalibConstantMap::const_iterator itcalib = imap.find(eb.rawId()); 00126 if ( itcalib == imap.end() ) { 00127 // FIXME -- throw error 00128 } 00129 EcalIntercalibConstant calib = (*itcalib); 00130 int sign = eb.zside()>0 ? 1 : 0; 00131 oldCalibs_barl[abs(eb.ieta())-1][eb.iphi()-1][sign] = calib; 00132 if (eb.iphi()==1) std::cout << "Read old constant for crystal " 00133 << " (" << eb.ieta() << "," << eb.iphi() 00134 << ") : " << calib << std::endl; 00135 } 00136 00137 // store eta value for each ring 00138 if (eb.ieta()>0 &&eb.iphi()==1) { 00139 const CaloCellGeometry *cellGeometry = barrelGeometry->getGeometry(*barrelIt); 00140 cellEta_[eb.ieta()-1] = cellGeometry->getPosition().eta(); 00141 } 00142 } 00143 00144 // loop over all endcap crystals 00145 endcapCells = geometry.getValidDetIds(DetId::Ecal, EcalEndcap); 00146 std::vector<DetId>::const_iterator endcapIt; 00147 for (endcapIt=endcapCells.begin(); endcapIt!=endcapCells.end(); endcapIt++) { 00148 const CaloCellGeometry *cellGeometry = endcapGeometry->getGeometry(*endcapIt); 00149 EEDetId ee(*endcapIt); 00150 int ix=ee.ix()-1; 00151 int iy=ee.iy()-1; 00152 00153 if (eventSet_==0) { 00154 // get the initial calibration constants 00155 EcalIntercalibConstantMap::const_iterator itcalib = imap.find(ee.rawId()); 00156 if ( itcalib == imap.end() ) { 00157 // FIXME -- throw error 00158 } 00159 EcalIntercalibConstant calib = (*itcalib); 00160 int sign = ee.zside()>0 ? 1 : 0; 00161 oldCalibs_endc[ix][iy][sign] = calib; 00162 if (ix==49) std::cout << "Read old constant for xcrystal " 00163 << " (" << ix << "," << iy 00164 << ") : " << calib << std::endl; 00165 } 00166 00167 // store all crystal positions 00168 cellPos_[ix][iy] = cellGeometry->getPosition(); 00169 cellPhi_[ix][iy] = cellGeometry->getPosition().phi(); 00170 00171 00172 00173 // calculate and store eta-phi area for each crystal front face Shoelace formuls 00174 const CaloCellGeometry::CornersVec& cellCorners (cellGeometry->getCorners()) ; 00175 cellArea_[ix][iy]=0.; 00176 for (int i=0; i<4; i++) { 00177 int iplus1 = i==3 ? 0 : i+1; 00178 cellArea_[ix][iy] += 00179 cellCorners[i].eta()*float(cellCorners[iplus1].phi()) - 00180 cellCorners[iplus1].eta()*float(cellCorners[i].phi()); 00181 } 00182 cellArea_[ix][iy] = cellArea_[ix][iy]/2.; 00183 00184 } 00185 00186 00187 // get eta for each endcap ring 00188 float eta_ring[kEndcEtaRings]; 00189 for (int ring=0; ring<kEndcEtaRings; ring++) { 00190 eta_ring[ring]=cellPos_[ring][50].eta(); 00191 } 00192 00193 // get eta boundaries for each endcap ring 00194 etaBoundary_[0]=1.479; 00195 etaBoundary_[kEndcEtaRings]=4.; 00196 for (int ring=1; ring<kEndcEtaRings; ring++) { 00197 etaBoundary_[ring]=(eta_ring[ring]+eta_ring[ring-1])/2.; 00198 } 00199 00200 // determine to which ring each endcap crystal belongs, 00201 // the number of crystals in each ring, 00202 // and the mean eta-phi area of the crystals in each ring 00203 for (int ring=0; ring<kEndcEtaRings; ring++) { 00204 nRing_[ring]=0; 00205 meanCellArea_[ring]=0.; 00206 for (int ix=0; ix<kEndcWedgesX; ix++) { 00207 for (int iy=0; iy<kEndcWedgesY; iy++) { 00208 if (fabs(cellPos_[ix][iy].eta())>etaBoundary_[ring] && 00209 fabs(cellPos_[ix][iy].eta())<etaBoundary_[ring+1]) { 00210 meanCellArea_[ring]+=cellArea_[ix][iy]; 00211 endcapRing_[ix][iy]=ring; 00212 nRing_[ring]++; 00213 } 00214 } 00215 } 00216 00217 meanCellArea_[ring]/=nRing_[ring]; 00218 00219 std::cout << nRing_[ring] << " crystals with mean area " << meanCellArea_[ring] << " in endcap ring " << ring << " (" << etaBoundary_[ring] << "<eta<" << etaBoundary_[ring+1] << ")" << std::endl; 00220 } 00221 00222 00223 00224 for (int ring=0; ring<kEndcEtaRings; ring++) { 00225 for (int i=0; i<kMaxEndciPhi; i++) phi_endc[i][ring]=0.; 00226 00227 float philast=-999.; 00228 for (int ip=0; ip<nRing_[ring]; ip++) { 00229 float phimin=999.; 00230 for (int ix=0; ix<kEndcWedgesX; ix++) { 00231 for (int iy=0; iy<kEndcWedgesY; iy++) { 00232 if (endcapRing_[ix][iy]==ring) { 00233 if (cellPhi_[ix][iy]<phimin && cellPhi_[ix][iy]>philast) { 00234 phimin=cellPhi_[ix][iy]; 00235 } 00236 } 00237 } 00238 } 00239 00240 phi_endc[ip][ring]=phimin; 00241 philast=phimin; 00242 } 00243 } 00244 00245 00246 }
Called at end of job.
Reimplemented from edm::EDAnalyzer.
Definition at line 251 of file PhiSymmetryCalibration.cc.
References funct::abs(), barrelCells, cellArea_, cellEta_, cellPos_, GenMuonPlsPt100GeV_cfg::cout, dummy, EcalBarrel, EcalEndcap, endcapCells, endcapRing_, lat::endl(), PV3DBase< T, PVType, FrameType >::eta(), etsum_barl_, etsum_endc_, etsumMean_barl_, etsumMean_endc_, eventSet_, f, fillHistos(), getKfactors(), EBDetId::ieta(), in, EBDetId::iphi(), EEDetId::ix(), EEDetId::iy(), k_barl_, k_endc_, kBarlRings, kBarlWedges, kEndcEtaRings, kEndcWedgesX, kEndcWedgesY, kSides, meanCellArea_, oldCalibs_barl, oldCalibs_endc, out, calibXMLwriter::writeLine(), EEDetId::zside(), and EBDetId::zside().
00252 { 00253 00254 edm::LogWarning("Calibration") << "[PhiSymmetryCalibration] At end of job"; 00255 return; 00256 00257 if (eventSet_==1) { 00258 // calculate factors to convert from fractional deviation of ET sum from 00259 // the mean to the estimate of the miscalibration factor 00260 getKfactors(); 00261 00262 00263 std::ofstream k_barl_out("k_barl.dat", ios::out); 00264 for (int ieta=0; ieta<kBarlRings; ieta++) 00265 k_barl_out << ieta << " " << k_barl_[ieta] << endl; 00266 k_barl_out.close(); 00267 00268 std::ofstream k_endc_out("k_endc.dat", ios::out); 00269 for (int ring=0; ring<kEndcEtaRings; ring++) 00270 k_endc_out << ring << " " << k_endc_[ring] << endl; 00271 k_endc_out.close(); 00272 } 00273 00274 if (eventSet_!=0) { 00275 00276 //output ET sums 00277 stringstream etsum_file_barl; 00278 etsum_file_barl << "etsum_barl_"<<eventSet_<<".dat"; 00279 00280 std::ofstream etsum_barl_out(etsum_file_barl.str().c_str(),ios::out); 00281 00282 for (int ieta=0; ieta<kBarlRings; ieta++) { 00283 for (int iphi=0; iphi<kBarlWedges; iphi++) { 00284 for (int sign=0; sign<kSides; sign++) { 00285 etsum_barl_out << eventSet_ << " " << ieta << " " << iphi << " " << sign << " "<< etsum_barl_[ieta][iphi][sign] << endl; 00286 } 00287 } 00288 } 00289 etsum_barl_out.close(); 00290 00291 stringstream etsum_file_endc; 00292 etsum_file_endc << "etsum_endc_"<<eventSet_<<".dat"; 00293 00294 std::ofstream etsum_endc_out(etsum_file_endc.str().c_str(),ios::out); 00295 for (int ix=0; ix<kEndcWedgesX; ix++) { 00296 for (int iy=0; iy<kEndcWedgesY; iy++) { 00297 int ring = endcapRing_[ix][iy]; 00298 if (ring!=-1) { 00299 for (int sign=0; sign<kSides; sign++) { 00300 etsum_endc_out << eventSet_ << " " << ix << " " << iy << " " << sign << " " << etsum_endc_[ix][iy][sign] << endl; 00301 } 00302 } 00303 } 00304 } 00305 etsum_endc_out.close(); 00306 00307 } else { // eventset =0 00308 00309 for (int ieta=0; ieta<kBarlRings; ieta++) { 00310 for (int iphi=0; iphi<kBarlWedges; iphi++) { 00311 for (int sign=0; sign<kSides; sign++) { 00312 etsum_barl_[ieta][iphi][sign]=0.; 00313 } 00314 } 00315 } 00316 00317 for (int ix=0; ix<kEndcWedgesX; ix++) { 00318 for (int iy=0; iy<kEndcWedgesY; iy++) { 00319 for (int sign=0; sign<kSides; sign++) { 00320 etsum_endc_[ix][iy][sign]=0.; 00321 } 00322 } 00323 } 00324 00325 //read in ET sums 00326 00327 int ieta,iphi,sign,ix,iy,dummy; 00328 double etsum; 00329 std::ifstream etsum_barl_in("etsum_barl.dat", ios::in); 00330 while ( etsum_barl_in >> dummy >> ieta >> iphi >> sign >> etsum) { 00331 etsum_barl_[ieta][iphi][sign]+=etsum; 00332 } 00333 etsum_barl_in.close(); 00334 00335 std::ifstream etsum_endc_in("etsum_endc.dat", ios::in); 00336 while (etsum_endc_in >> dummy >> ix >> iy >> sign >> etsum) { 00337 etsum_endc_[ix][iy][sign]+=etsum; 00338 } 00339 etsum_endc_in.close(); 00340 00341 std::ifstream k_barl_in("k_barl.dat", ios::in); 00342 for (int ieta=0; ieta<kBarlRings; ieta++) { 00343 k_barl_in >> dummy >> k_barl_[ieta]; 00344 } 00345 k_barl_in.close(); 00346 00347 std::ifstream k_endc_in("k_endc.dat", ios::in); 00348 for (int ring=0; ring<kEndcEtaRings; ring++) { 00349 k_endc_in >> dummy >> k_endc_[ring]; 00350 } 00351 k_endc_in.close(); 00352 00353 00354 00355 for (int ix=0; ix<kEndcWedgesX; ix++) { 00356 for (int iy=0; iy<kEndcWedgesY; iy++) { 00357 00358 int ring = endcapRing_[ix][iy]; 00359 00360 if (ring!=-1) { 00361 for (int sign=0; sign<kSides; sign++) { 00362 00363 etsum_endc_[ix][iy][sign]*=meanCellArea_[ring]/cellArea_[ix][iy]; 00364 } 00365 } 00366 } 00367 } 00368 00369 // fill output ET sum histograms 00370 fillHistos(); 00371 00372 std::ofstream etsumMean_barl_out("etsumMean_barl.dat",ios::out); 00373 for (int ieta=0; ieta<kBarlRings; ieta++) { 00374 etsumMean_barl_out << cellEta_[ieta] << " " << etsumMean_barl_[ieta] << endl; 00375 } 00376 etsumMean_barl_out.close(); 00377 00378 std::ofstream etsumMean_endc_out("etsumMean_endc.dat",ios::out); 00379 for (int ring=0; ring<kEndcEtaRings; ring++) { 00380 etsumMean_endc_out << cellPos_[ring][50].eta() << " " << etsumMean_endc_[ring] << endl; 00381 } 00382 etsumMean_endc_out.close(); 00383 00384 std::ofstream area_out("area.dat",ios::out); 00385 for (int ring=0; ring<kEndcEtaRings; ring++) { 00386 area_out << meanCellArea_[ring] << endl; 00387 } 00388 area_out.close(); 00389 00390 // Determine barrel calibration constants 00391 float epsilon_M_barl[kBarlRings][kBarlWedges][kSides]; 00392 for (int ieta=0; ieta<kBarlRings; ieta++) { 00393 for (int iphi=0; iphi<kBarlWedges; iphi++) { 00394 for (int sign=0; sign<kSides; sign++) { 00395 float etsum = etsum_barl_[ieta][iphi][sign]; 00396 float epsilon_T = (etsum/etsumMean_barl_[ieta])-1; 00397 epsilon_M_barl[ieta][iphi][sign] = epsilon_T/k_barl_[ieta]; 00398 } 00399 } 00400 } 00401 00402 // Determine endcap calibration constants 00403 float epsilon_M_endc[kEndcWedgesX][kEndcWedgesY][2]; 00404 for (int ix=0; ix<kEndcWedgesX; ix++) { 00405 for (int iy=0; iy<kEndcWedgesY; iy++) { 00406 int ring = endcapRing_[ix][iy]; 00407 if (ring!=-1) { 00408 for (int sign=0; sign<kSides; sign++) { 00409 float etsum = etsum_endc_[ix][iy][sign]; 00410 float epsilon_T = (etsum/etsumMean_endc_[ring])-1; 00411 epsilon_M_endc[ix][iy][sign] = epsilon_T/k_endc_[ring]; 00412 } 00413 } else { 00414 epsilon_M_endc[ix][iy][0] = 0.; 00415 epsilon_M_endc[ix][iy][1] = 0.; 00416 } 00417 } 00418 } 00419 00420 // Write new calibration constants 00421 00422 calibXMLwriter barrelWriter(EcalBarrel); 00423 calibXMLwriter endcapWriter(EcalEndcap); 00424 00425 double newCalibs_barl[kEndcWedgesX][kEndcWedgesX][kSides]; 00426 double newCalibs_endc[kEndcWedgesX][kEndcWedgesX][kSides]; 00427 00428 std::vector<TH1F*> miscal_resid_barl_histos(kBarlRings); 00429 std::vector<TH2F*> correl_barl_histos(kBarlRings); 00430 00431 for (int ieta=0; ieta<kBarlRings; ieta++) { 00432 00433 ostringstream t1; 00434 t1<<"mr_barl_"<<ieta+1; 00435 miscal_resid_barl_histos[ieta] = new TH1F(t1.str().c_str(),"",50,.8,1.2); 00436 ostringstream t2; 00437 t2<<"co_barl_"<<ieta+1; 00438 correl_barl_histos[ieta] = new TH2F(t1.str().c_str(),"",50,.8,1.2,50,.8,1.2); 00439 } 00440 00441 std::vector<TH1F*> miscal_resid_endc_histos(kEndcEtaRings); 00442 std::vector<TH2F*> correl_endc_histos(kEndcEtaRings); 00443 00444 for (int ring=0; ring<kEndcEtaRings; ring++) { 00445 ostringstream t1; 00446 t1<<"mr_endc_"<< ring+1; 00447 miscal_resid_endc_histos[ring] = new TH1F(t1.str().c_str(),"",50,.8,1.2); 00448 ostringstream t2; 00449 t2<<"co_endc_"<<ring+1; 00450 correl_endc_histos[ring] = new TH2F(t2.str().c_str(),"",50,.8,1.2,50,.8,1.2); 00451 } 00452 00453 std::vector<DetId>::const_iterator barrelIt=barrelCells.begin(); 00454 TFile eehistof("eehisto.root","recreate"); 00455 TH1D ebhisto("eb","eb",100, -2,2); 00456 for (; barrelIt!=barrelCells.end(); barrelIt++) { 00457 EBDetId eb(*barrelIt); 00458 int ieta = abs(eb.ieta())-1; 00459 int iphi = eb.iphi()-1; 00460 int sign = eb.zside()>0 ? 1 : 0; 00461 newCalibs_barl[ieta][iphi][sign] = 00462 oldCalibs_barl[ieta][iphi][sign]/ 00463 (1+epsilon_M_barl[ieta][iphi][sign]); 00464 ebhisto.Fill(newCalibs_barl[ieta][iphi][sign]); 00465 barrelWriter.writeLine(eb,newCalibs_barl[ieta][iphi][sign]); 00466 miscal_resid_barl_histos[ieta]->Fill(newCalibs_barl[ieta][iphi][sign]); 00467 correl_barl_histos[ieta]->Fill(oldCalibs_barl[ieta][iphi][sign],1+epsilon_M_barl[ieta][iphi][sign]); 00468 if (iphi==1) { 00469 std::cout << "Calib constant for barrel crystal " 00470 << " (" << eb.ieta() << "," << eb.iphi() << ") changed from " 00471 << oldCalibs_barl[ieta][iphi][sign] << " to " 00472 << newCalibs_barl[ieta][iphi][sign] << std::endl; 00473 } 00474 } 00475 00476 00477 TH1D eehisto("ee","ee",100, -2,2); 00478 std::vector<DetId>::const_iterator endcapIt=endcapCells.begin(); 00479 for (; endcapIt!=endcapCells.end(); endcapIt++) { 00480 EEDetId ee(*endcapIt); 00481 int ix = ee.ix()-1; 00482 int iy = ee.iy()-1; 00483 int sign = ee.zside()>0 ? 1 : 0; 00484 newCalibs_endc[ix][iy][sign] = 00485 oldCalibs_endc[ix][iy][sign]/ 00486 (1+epsilon_M_endc[ix][iy][sign])/1.03; 00487 eehisto.Fill(newCalibs_endc[ix][iy][sign]); 00488 endcapWriter.writeLine(ee,newCalibs_endc[ix][iy][sign]); 00489 miscal_resid_endc_histos[endcapRing_[ix][iy]]->Fill(newCalibs_endc[ix][iy][sign]); 00490 correl_endc_histos[endcapRing_[ix][iy]]->Fill(oldCalibs_endc[ix][iy][sign],1+epsilon_M_endc[ix][iy][sign]); 00491 if (ix==50) { 00492 std::cout << "Calib constant for endcap crystal " 00493 << " (" << ix << "," << iy << "," << sign << ") changed from " 00494 << oldCalibs_endc[ix][iy][sign] << " to " 00495 << newCalibs_endc[ix][iy][sign] << std::endl; 00496 } 00497 } 00498 00499 eehisto.Write(); 00500 ebhisto.Write(); 00501 eehistof.Close(); 00502 00503 00504 // Output histograms of residual miscalibrations 00505 TFile f("PhiSymmetryCalibration_miscal_resid.root","recreate"); 00506 for (int ieta=0; ieta<kBarlRings; ieta++) { 00507 miscal_resid_barl_histos[ieta]->Fit("gaus"); 00508 miscal_resid_barl_histos[ieta]->Write(); 00509 correl_barl_histos[ieta]->Write(); 00510 00511 delete miscal_resid_barl_histos[ieta]; 00512 delete correl_barl_histos[ieta]; 00513 } 00514 for (int ring=0; ring<kEndcEtaRings; ring++) { 00515 miscal_resid_endc_histos[ring]->Fit("gaus"); 00516 miscal_resid_endc_histos[ring]->Write(); 00517 correl_endc_histos[ring]->Write(); 00518 00519 delete miscal_resid_endc_histos[ring]; 00520 delete correl_endc_histos[ring]; 00521 00522 } 00523 f.Close(); 00524 00525 00526 } 00527 00528 }
| void PhiSymmetryCalibration::fillHistos | ( | ) | [private] |
Definition at line 704 of file PhiSymmetryCalibration.cc.
References endcapRing_, etsum_barl_, etsum_endc_, etsumMean_barl_, etsumMean_endc_, f, WenuSkim_TriggerBit_cff::high, kBarlRings, kBarlWedges, kEndcEtaRings, kEndcWedgesX, kEndcWedgesY, kSides, nRing_, and t.
Referenced by endJob().
00705 { 00706 00707 00708 TFile f("PhiSymmetryCalibration.root","recreate"); 00709 00710 std::vector<TH1F*> etsum_barl_histos(kBarlRings); 00711 00712 for (int ieta=0; ieta<kBarlRings; ieta++) { 00713 00714 // Determine ranges of ET sums to get histo bounds and book histos (barrel) 00715 float low=999999.; 00716 float high=0.; 00717 for (int iphi=0; iphi<kBarlWedges; iphi++) { 00718 for (int sign=0; sign<kSides; sign++) { 00719 float etsum = etsum_barl_[ieta][iphi][sign]; 00720 if (etsum<low) low=etsum; 00721 if (etsum>high) high=etsum; 00722 } 00723 } 00724 00725 ostringstream t; 00726 t<< "etsum_barl_" << ieta+1; 00727 00728 etsum_barl_histos[ieta]=new TH1F(t.str().c_str(),"",50,low-.2*low,high+.1*high); 00729 00730 // Fill barrel ET sum histos 00731 etsumMean_barl_[ieta]=0.; 00732 for (int iphi=0; iphi<kBarlWedges; iphi++) { 00733 for (int sign=0; sign<kSides; sign++) { 00734 float etsum = etsum_barl_[ieta][iphi][sign]; 00735 etsum_barl_histos[ieta]->Fill(etsum); 00736 etsumMean_barl_[ieta]+=etsum; 00737 } 00738 } 00739 etsum_barl_histos[ieta]->Fit("gaus"); 00740 etsum_barl_histos[ieta]->Write(); 00741 etsumMean_barl_[ieta]/=720.; 00742 delete etsum_barl_histos[ieta]; 00743 } 00744 00745 std::vector<TH1F*> etsum_endc_histos(kEndcEtaRings); 00746 00747 for (int ring=0; ring<kEndcEtaRings; ring++) { 00748 00749 // Determine ranges of ET sums to get histo bounds and book histos (endcap) 00750 float low=999999.; 00751 float high=0.; 00752 for (int ix=0; ix<kEndcWedgesX; ix++) { 00753 for (int iy=0; iy<kEndcWedgesY; iy++) { 00754 if (endcapRing_[ix][iy]==ring) { 00755 for (int sign=0; sign<kSides; sign++) { 00756 float etsum = etsum_endc_[ix][iy][sign]; 00757 if (etsum<low) low=etsum; 00758 if (etsum>high) high=etsum; 00759 } 00760 } 00761 } 00762 } 00763 ostringstream t; 00764 t<<"etsum_endc_" << ring+1; 00765 etsum_endc_histos[ring]= new TH1F(t.str().c_str(),"",50,low-.2*low,high+.1*high); 00766 00767 // Fill endcap ET sum histos* 00768 etsumMean_endc_[ring]=0.; 00769 for (int ix=0; ix<kEndcWedgesX; ix++) { 00770 for (int iy=0; iy<kEndcWedgesY; iy++) { 00771 if (endcapRing_[ix][iy]==ring) { 00772 for (int sign=0; sign<kSides; sign++) { 00773 float etsum = etsum_endc_[ix][iy][sign]; 00774 etsum_endc_histos[ring]->Fill(etsum); 00775 etsumMean_endc_[ring]+=etsum; 00776 } 00777 } 00778 } 00779 } 00780 etsum_endc_histos[ring]->Fit("gaus"); 00781 etsum_endc_histos[ring]->Write(); 00782 etsumMean_endc_[ring]/=(float(nRing_[ring]*2)); 00783 delete etsum_endc_histos[ring]; 00784 } 00785 00786 f.Close(); 00787 00788 00789 }
| void PhiSymmetryCalibration::getKfactors | ( | ) | [private] |
Definition at line 619 of file PhiSymmetryCalibration.cc.
References GenMuonPlsPt100GeV_cfg::cout, lat::endl(), etsum_barl_miscal_, etsum_endc_miscal_, f, k_barl_, k_endc_, kBarlRings, kEndcEtaRings, kNMiscalBins, miscal_, and t.
Referenced by endJob().
00620 { 00621 00622 float epsilon_T[kNMiscalBins]; 00623 float epsilon_M[kNMiscalBins]; 00624 00625 std::vector<TGraph*> k_barl_graph(kBarlRings); 00626 std::vector<TCanvas*> k_barl_plot(kBarlRings); 00627 00628 for (int ieta=0; ieta<kBarlRings; ieta++) { 00629 for (int imiscal=0; imiscal<kNMiscalBins; imiscal++) { 00630 epsilon_T[imiscal] = etsum_barl_miscal_[imiscal][ieta]/etsum_barl_miscal_[10][ieta] - 1.; 00631 epsilon_M[imiscal]=miscal_[imiscal]-1.; 00632 } 00633 00634 k_barl_graph[ieta] = new TGraph (kNMiscalBins,epsilon_M,epsilon_T); 00635 k_barl_graph[ieta]->Fit("pol1"); 00636 00637 00638 ostringstream t; 00639 t<< "k_barl_" << ieta+1; 00640 k_barl_plot[ieta] = new TCanvas(t.str().c_str(),""); 00641 k_barl_plot[ieta]->SetFillColor(10); 00642 k_barl_plot[ieta]->SetGrid(); 00643 k_barl_graph[ieta]->SetMarkerSize(1.); 00644 k_barl_graph[ieta]->SetMarkerColor(4); 00645 k_barl_graph[ieta]->SetMarkerStyle(20); 00646 k_barl_graph[ieta]->GetXaxis()->SetLimits(-.06,.06); 00647 k_barl_graph[ieta]->GetXaxis()->SetTitleSize(.05); 00648 k_barl_graph[ieta]->GetYaxis()->SetTitleSize(.05); 00649 k_barl_graph[ieta]->GetXaxis()->SetTitle("#epsilon_{M}"); 00650 k_barl_graph[ieta]->GetYaxis()->SetTitle("#epsilon_{T}"); 00651 k_barl_graph[ieta]->Draw("AP"); 00652 00653 00654 k_barl_[ieta] = k_barl_graph[ieta]->GetFunction("pol1")->GetParameter(1); 00655 00656 } 00657 00658 00659 std::vector<TGraph*> k_endc_graph(kBarlRings); 00660 std::vector<TCanvas*> k_endc_plot(kBarlRings); 00661 00662 for (int ring=0; ring<kEndcEtaRings; ring++) { 00663 for (int imiscal=0; imiscal<kNMiscalBins; imiscal++) { 00664 epsilon_T[imiscal] = etsum_endc_miscal_[imiscal][ring]/etsum_endc_miscal_[10][ring] - 1.; 00665 epsilon_M[imiscal]=miscal_[imiscal]-1.; 00666 } 00667 k_endc_graph[ring] = new TGraph (kNMiscalBins,epsilon_M,epsilon_T); 00668 k_endc_graph[ring]->Fit("pol1"); 00669 00670 ostringstream t; 00671 t<< "k_endc_"<< ring+1; 00672 k_endc_plot[ring] = new TCanvas(t.str().c_str(),""); 00673 k_endc_plot[ring]->SetFillColor(10); 00674 k_endc_plot[ring]->SetGrid(); 00675 k_endc_graph[ring]->SetMarkerSize(1.); 00676 k_endc_graph[ring]->SetMarkerColor(4); 00677 k_endc_graph[ring]->SetMarkerStyle(20); 00678 k_endc_graph[ring]->GetXaxis()->SetLimits(-.06,.06); 00679 k_endc_graph[ring]->GetXaxis()->SetTitleSize(.05); 00680 k_endc_graph[ring]->GetYaxis()->SetTitleSize(.05); 00681 k_endc_graph[ring]->GetXaxis()->SetTitle("#epsilon_{M}"); 00682 k_endc_graph[ring]->GetYaxis()->SetTitle("#epsilon_{T}"); 00683 k_endc_graph[ring]->Draw("AP"); 00684 00685 k_endc_[ring] = k_endc_graph[ring]->GetFunction("pol1")->GetParameter(1); 00686 std::cout << "k_endc_[" << ring << "]=" << k_endc_[ring] << std::endl; 00687 } 00688 00689 TFile f("PhiSymmetryCalibration_kFactors.root","recreate"); 00690 for (int ieta=0; ieta<kBarlRings; ieta++) { 00691 k_barl_plot[ieta]->Write(); 00692 delete k_barl_plot[ieta]; 00693 delete k_barl_graph[ieta]; 00694 } 00695 for (int ring=0; ring<kEndcEtaRings; ring++) { 00696 k_endc_plot[ring]->Write(); 00697 delete k_endc_plot[ring]; 00698 delete k_endc_graph[ring]; 00699 } 00700 f.Close(); 00701 00702 }
std::vector<DetId> PhiSymmetryCalibration::barrelCells [private] |
std::string PhiSymmetryCalibration::barrelHits_ [private] |
double PhiSymmetryCalibration::cellArea_[kEndcWedgesX][kEndcWedgesY] [private] |
Definition at line 94 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and endJob().
double PhiSymmetryCalibration::cellEta_[kBarlRings] [private] |
Definition at line 90 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and endJob().
double PhiSymmetryCalibration::cellPhi_[kEndcWedgesX][kEndcWedgesY] [private] |
Definition at line 92 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and endJob().
std::string PhiSymmetryCalibration::ecalHitsProducer_ [private] |
double PhiSymmetryCalibration::eCut_barl_ [private] |
double PhiSymmetryCalibration::eCut_endc_ [private] |
std::vector<DetId> PhiSymmetryCalibration::endcapCells [private] |
std::string PhiSymmetryCalibration::endcapHits_ [private] |
Definition at line 97 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), endJob(), and fillHistos().
double PhiSymmetryCalibration::etaBoundary_[kEndcEtaRings+1] [private] |
double PhiSymmetryCalibration::etsum_barl_[kBarlRings][kBarlWedges][kSides] [private] |
Definition at line 80 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), endJob(), and fillHistos().
double PhiSymmetryCalibration::etsum_barl_miscal_[kNMiscalBins][kBarlRings] [private] |
Definition at line 86 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and getKfactors().
double PhiSymmetryCalibration::etsum_endc_[kEndcWedgesX][kEndcWedgesX][kSides] [private] |
Definition at line 81 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), endJob(), and fillHistos().
double PhiSymmetryCalibration::etsum_endc_miscal_[kNMiscalBins][kEndcEtaRings] [private] |
Definition at line 87 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and getKfactors().
double PhiSymmetryCalibration::etsumMean_barl_[kBarlRings] [private] |
double PhiSymmetryCalibration::etsumMean_endc_[kEndcEtaRings] [private] |
int PhiSymmetryCalibration::eventSet_ [private] |
Definition at line 119 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and endJob().
double PhiSymmetryCalibration::k_barl_[kBarlRings] [private] |
double PhiSymmetryCalibration::k_endc_[kEndcEtaRings] [private] |
const int PhiSymmetryCalibration::kBarlRings = 85 [static, private] |
Definition at line 69 of file PhiSymmetryCalibration.h.
Referenced by beginJob(), endJob(), fillHistos(), and getKfactors().
const int PhiSymmetryCalibration::kBarlWedges = 360 [static, private] |
Definition at line 70 of file PhiSymmetryCalibration.h.
Referenced by beginJob(), endJob(), and fillHistos().
const int PhiSymmetryCalibration::kEndcEtaRings = 39 [static, private] |
Definition at line 76 of file PhiSymmetryCalibration.h.
Referenced by beginJob(), endJob(), fillHistos(), and getKfactors().
const int PhiSymmetryCalibration::kEndcWedgesX = 100 [static, private] |
Definition at line 73 of file PhiSymmetryCalibration.h.
Referenced by beginJob(), endJob(), and fillHistos().
const int PhiSymmetryCalibration::kEndcWedgesY = 100 [static, private] |
Definition at line 74 of file PhiSymmetryCalibration.h.
Referenced by beginJob(), endJob(), and fillHistos().
const int PhiSymmetryCalibration::kMaxEndciPhi = 360 [static, private] |
const int PhiSymmetryCalibration::kNMiscalBins = 21 [static, private] |
Definition at line 77 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and getKfactors().
const int PhiSymmetryCalibration::kSides = 2 [static, private] |
Definition at line 71 of file PhiSymmetryCalibration.h.
Referenced by beginJob(), endJob(), and fillHistos().
double PhiSymmetryCalibration::meanCellArea_[kEndcEtaRings] [private] |
Definition at line 95 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and endJob().
double PhiSymmetryCalibration::miscal_[kNMiscalBins] [private] |
Definition at line 103 of file PhiSymmetryCalibration.h.
Referenced by analyze(), beginJob(), and getKfactors().
int PhiSymmetryCalibration::nevent [private] |
int PhiSymmetryCalibration::nRing_[kEndcEtaRings] [private] |
double PhiSymmetryCalibration::oldCalibs_barl[kBarlRings][kBarlWedges][kSides] [private] |
double PhiSymmetryCalibration::oldCalibs_endc[kEndcWedgesX][kEndcWedgesY][kSides] [private] |
float PhiSymmetryCalibration::phi_endc[kMaxEndciPhi][kEndcEtaRings] [private] |
1.5.4