#include <PhiSymmetryCalibration_step2_SM.h>

Inheritance diagram for PhiSymmetryCalibration_step2_SM:

Public Member Functions
void	analyze (const edm::Event &, const edm::EventSetup &)

void	beginJob ()

void	endJob ()

void	fillConstantsHistos ()

void	fillHistos ()

void	outResidHistos ()

	PhiSymmetryCalibration_step2_SM (const edm::ParameterSet &iConfig)

void	readEtSums ()

void	setUp (const edm::EventSetup &setup)

void	setupResidHistos ()

	~PhiSymmetryCalibration_step2_SM ()

Public Member Functions inherited from edm::EDAnalyzer
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzer ()

ModuleDescription const &	moduleDescription () const

std::string	workerType () const

virtual	~EDAnalyzer ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Private Attributes
std::vector< DetId >	barrelCells

std::vector< TH2F * >	correl_barl_histos

std::vector< TH2F * >	correl_endc_histos

EcalGeomPhiSymHelper	e_

std::vector< DetId >	endcapCells

float	epsilon_M_barl [kBarlRings][kBarlWedges][kSides]

double	epsilon_M_barl_SM_ [kBarlRings][int(kBarlWedges/nscx)][kSides]

float	epsilon_M_endc [kEndcWedgesX][kEndcWedgesY][kSides]

double	esum_barl_ [kBarlRings][kBarlWedges][kSides]

double	esum_endc_ [kEndcWedgesX][kEndcWedgesX][kSides]

double	esumMean_barl_ [kBarlRings]

double	esumMean_endc_ [kEndcEtaRings]

double	etsum_barl_ [kBarlRings][kBarlWedges][kSides]

double	etsum_barl_SM_ [kBarlRings][int(kBarlWedges/nscx)][kSides]

double	etsum_endc_ [kEndcWedgesX][kEndcWedgesX][kSides]

double	etsum_endc_uncorr [kEndcWedgesX][kEndcWedgesX][kSides]

double	etsumMean_barl_ [kBarlRings]

double	etsumMean_barl_SM_ [kBarlRings]

double	etsumMean_endc_ [kEndcEtaRings]

bool	firstpass_

bool	have_initial_miscalib_

std::string	initialmiscalibfile_

double	k_barl_ [kBarlRings]

double	k_endc_ [kEndcEtaRings]

std::vector< TH1F * >	miscal_resid_barl_histos
	res miscalib histos More...

std::vector< TH1F * >	miscal_resid_endc_histos

EcalIntercalibConstants	miscalib_
	initial miscalibration applied if any) More...

int	nBads_barl_SM_ [kBarlRings][int(kBarlWedges/nscx)][kSides]

EcalIntercalibConstants	newCalibs_
	calib constants that we are going to calculate More...

unsigned int	nhits_barl_ [kBarlRings][kBarlWedges][kSides]

unsigned int	nhits_endc_ [kEndcWedgesX][kEndcWedgesX][kSides]

std::string	oldcalibfile_

EcalIntercalibConstants	oldCalibs_
	the old calibration constants (when reiterating, the last ones derived) More...

float	rawconst_barl [kBarlRings][kBarlWedges][kSides]

float	rawconst_endc [kEndcWedgesX][kEndcWedgesX][kSides]

bool	reiteration_

int	statusThreshold_

Static Private Attributes
static const int	nscx =20

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 12 of file PhiSymmetryCalibration_step2_SM.h.

Constructor & Destructor Documentation

PhiSymmetryCalibration_step2_SM::PhiSymmetryCalibration_step2_SM ( const edm::ParameterSet & iConfig )

Definition at line 26 of file PhiSymmetryCalibration_step2_SM.cc.

References edm::ParameterSet::getUntrackedParameter(), and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                                               {
 
   statusThreshold_ =
        iConfig.getUntrackedParameter<int>("statusThreshold",0);
   have_initial_miscalib_=
        iConfig.getUntrackedParameter<bool>("haveInitialMiscalib",false);
   initialmiscalibfile_=
     iConfig.getUntrackedParameter<std::string>("initialmiscalibfile",
                            "InitialMiscalib.xml"); 
   oldcalibfile_=
     iConfig.getUntrackedParameter<std::string>("oldcalibfile",
                            "EcalIntercalibConstants.xml");
   reiteration_ = iConfig.getUntrackedParameter<bool>("reiteration",false);
   firstpass_=true;
 }

PhiSymmetryCalibration_step2_SM::~PhiSymmetryCalibration_step2_SM ( )

Definition at line 23 of file PhiSymmetryCalibration_step2_SM.cc.

23 {}

Member Function Documentation

void PhiSymmetryCalibration_step2_SM::analyze	(	const edm::Event &	ev,
		const edm::EventSetup &	se
	)

virtual

Implements edm::EDAnalyzer.

Definition at line 42 of file PhiSymmetryCalibration_step2_SM.cc.

                                                 {
 
   if (firstpass_) {
     setUp(se);
     firstpass_=false;
   }
 }

void PhiSymmetryCalibration_step2_SM::beginJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 138 of file PhiSymmetryCalibration_step2_SM.cc.

References kBarlRings, kBarlWedges, kEndcWedgesX, kEndcWedgesY, and kSides.

                                               {
   
 
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     for (int iphi=0; iphi<kBarlWedges; iphi++) {
       for (int sign=0; sign<kSides; sign++) {
     int iphi_r=int(iphi/nscx);
     
     etsum_barl_[ieta][iphi][sign]=0.;
     nhits_barl_[ieta][iphi][sign]=0;
     esum_barl_[ieta][iphi][sign]=0.;
     etsum_barl_SM_[ieta][iphi_r][sign]=0;
     nBads_barl_SM_[ieta][iphi_r][sign]=0;
     epsilon_M_barl_SM_[ieta][iphi_r][sign]=0;
       }
     }                       
     etsumMean_barl_SM_[ieta]=0.;
   }
 
   for (int ix=0; ix<kEndcWedgesX; ix++) {
     for (int iy=0; iy<kEndcWedgesY; iy++) {
       for (int sign=0; sign<kSides; sign++) {
     etsum_endc_[ix][iy][sign]=0.;
     nhits_endc_[ix][iy][sign]=0;
     esum_endc_[ix][iy][sign]=0.;
 
       }
     }
   }
 
   readEtSums();
   setupResidHistos();
 }

void PhiSymmetryCalibration_step2_SM::endJob ( void )

virtual

this is the new constant, or better, the correction to be applied to the old constant

residual miscalibraition / expected precision

Reimplemented from edm::EDAnalyzer.

Definition at line 172 of file PhiSymmetryCalibration_step2_SM.cc.

References funct::abs(), EcalCondHeader::datasource_, EcalCondHeader::date_, EBDetId::ieta(), EBDetId::iphi(), EEDetId::ix(), EEDetId::iy(), kBarlRings, kBarlWedges, kEndcEtaRings, kEndcWedgesX, kEndcWedgesY, kSides, EcalCondHeader::method_, dbtoconf::out, relativeConstraints::ring, EcalCondHeader::since_, AlCaHLTBitMon_QueryRunRegistry::string, EcalCondHeader::tag_, EcalCondHeader::version_, EcalFloatCondObjectContainerXMLTranslator::writeXML(), EBDetId::zside(), and EEDetId::zside().

                                             {
 
   if (firstpass_) {
     edm::LogError("PhiSym")<< "Must process at least one event-Exiting" <<endl;
     return;
       
   }
 
   // Here the real calculation of constants happens
 
   // perform the area correction for endcap etsum
   // NOT  USED  ANYMORE
 
   /*
   for (int ix=0; ix<kEndcWedgesX; ix++) {
     for (int iy=0; iy<kEndcWedgesY; iy++) {
 
       int ring = e_.endcapRing_[ix][iy];
 
       if (ring!=-1) {
     for (int sign=0; sign<kSides; sign++) {
       etsum_endc_uncorr[ix][iy][sign] = etsum_endc_[ix][iy][sign];
       etsum_endc_[ix][iy][sign]*=meanCellArea_[ring]/cellArea_[ix][iy];
     }
       }
     }
   }
   */
 
   // ETsum histos, maps and other usefull histos (area,...)
   // are filled and saved here
   fillHistos();
 
   // write ETsum mean for all rings
   std::ofstream etsumMean_barl_out("etsumMean_barl.dat",ios::out);
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     etsumMean_barl_out << ieta << " " << etsumMean_barl_[ieta] << endl;
   }
   etsumMean_barl_out.close();
 
   std::ofstream etsumMean_endc_out("etsumMean_endc.dat",ios::out);
   for (int ring=0; ring<kEndcEtaRings; ring++) {
     etsumMean_endc_out << e_.cellPos_[ring][50].eta() << " " << etsumMean_endc_[ring] << endl;
   }
   etsumMean_endc_out.close();
   
 
   // determine barrel calibration constants
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     for (int iphi=0; iphi<kBarlWedges; iphi++) {
       for (int sign=0; sign<kSides; sign++) {
     
     
     // sc 
     int iphi_r = int(iphi/nscx);
       
       //if(nBads_barl_SM_[ieta][iphi_r][sign]>0){
       //  std::cout << "ETSUM" << etsum_barl_SM_[ieta][iphi_r][sign] << "  " <<ieta << " " << iphi_r << " " << sign << "  " << nBads_barl_SM_[ieta][iphi_r][sign]<< endl;
       //}      
         
     float epsilon_T_SM = 
        etsum_barl_SM_[ieta][iphi_r][sign] /etsumMean_barl_SM_[ieta] -1.;
     
     epsilon_M_barl_SM_[ieta][iphi_r][sign] = epsilon_T_SM/k_barl_[ieta];
     
     if(e_.goodCell_barl[ieta][iphi][sign]){
       float etsum = etsum_barl_[ieta][iphi][sign];
       float epsilon_T = (etsum/etsumMean_barl_[ieta]) - 1.;
       rawconst_barl[ieta][iphi][sign]  = epsilon_T + 1.;
       epsilon_M_barl[ieta][iphi][sign] = epsilon_T/k_barl_[ieta];
       
 
     } else {
       rawconst_barl[ieta][iphi][sign]  = 1.;
       epsilon_M_barl[ieta][iphi][sign] = 0.;
     } //if
       } //sign
     } //iphi
   } //ieta
 
     // determine endcap calibration constants
   for (int ix=0; ix<kEndcWedgesX; ix++) {
     for (int iy=0; iy<kEndcWedgesY; iy++) {
       for (int sign=0; sign<kSides; sign++) {
     int ring = e_.endcapRing_[ix][iy];
     if (ring!=-1 && e_.goodCell_endc[ix][iy][sign]) {
       float etsum = etsum_endc_[ix][iy][sign];
       float epsilon_T = (etsum/etsumMean_endc_[ring]) - 1.;
       rawconst_endc[ix][iy][sign]  = epsilon_T + 1.;
       epsilon_M_endc[ix][iy][sign] = epsilon_T/k_endc_[ring];       
     } else {
       epsilon_M_endc[ix][iy][0] = 0.;
       epsilon_M_endc[ix][iy][1] = 0.;
       rawconst_endc[ix][iy][0]  = 1.;
       rawconst_endc[ix][iy][1]  = 1.;
     } //if
       } //sign
     } //iy
   } //ix
 
 
 
   // output sc calibration
   std::fstream scfile("sccalibration.dat",std::ios::out);
   for (int ieta =0; ieta< kBarlRings; ++ieta){
     for (int iphi_r =0; iphi_r< int(kBarlWedges/nscx);++iphi_r){
       for (int sign=0; sign<kSides; sign++) {
     scfile<< ieta << " " << iphi_r << " " <<sign << " " 
           << 1/(1+epsilon_M_barl_SM_[ieta][iphi_r][sign] )<< std::endl;
       }
     }
   }
 
   std::string newcalibfile("EcalIntercalibConstants_new.xml");
 
 
 
   TFile ehistof("ehistos.root","recreate");  
 
   TH1D ebhisto("eb","eb",100, 0.,2.);
 
   std::vector<DetId>::const_iterator barrelIt=barrelCells.begin();
   for (; barrelIt!=barrelCells.end(); barrelIt++) {
     EBDetId eb(*barrelIt);
     int ieta = abs(eb.ieta())-1;
     int iphi = eb.iphi()-1;
     int sign = eb.zside()>0 ? 1 : 0;
 
     newCalibs_[eb] =  oldCalibs_[eb]/(1+epsilon_M_barl[ieta][iphi][sign]);
 
     if(e_.goodCell_barl[ieta][iphi][sign]){
 
       ebhisto.Fill(newCalibs_[eb]);
     
       miscal_resid_barl_histos[ieta]->Fill(miscalib_[eb]*newCalibs_[eb]);
       correl_barl_histos[ieta]->Fill(miscalib_[eb],newCalibs_[eb]);
     
     }
       
   }// barrelit
 
   TH1D eehisto("ee","ee",100, 0.,2.);
   std::vector<DetId>::const_iterator endcapIt=endcapCells.begin();
 
   for (; endcapIt!=endcapCells.end(); endcapIt++) {
     EEDetId ee(*endcapIt);
     int ix = ee.ix()-1;
     int iy = ee.iy()-1;
     int sign = ee.zside()>0 ? 1 : 0;
       
     newCalibs_[ee] = oldCalibs_[ee]/(1+epsilon_M_endc[ix][iy][sign]);
       
       
     if(e_.goodCell_endc[ix][iy][sign]){
  
       eehisto.Fill(newCalibs_[ee]);
       miscal_resid_endc_histos[e_.endcapRing_[ix][iy]]->Fill(miscalib_[ee]*
                                  newCalibs_[ee]);;
 
       correl_endc_histos[e_.endcapRing_[ix][iy]]->Fill(miscalib_[ee],
                                newCalibs_[ee]);
    
     }
   }//endcapit
   // Write xml file
   EcalCondHeader header;
   header.method_="phi symmetry";
   header.version_="0";
   header.datasource_="testdata";
   header.since_=1;
   header.tag_="unknown";
   header.date_="Mar 24 1973";
  
   EcalIntercalibConstantsXMLTranslator::writeXML(newcalibfile,header,
                          newCalibs_ );  
 
   eehisto.Write();
   ebhisto.Write();
   ehistof.Close();
 
   fillConstantsHistos();
   
   outResidHistos();
   
 }

void PhiSymmetryCalibration_step2_SM::fillConstantsHistos ( )

Definition at line 364 of file PhiSymmetryCalibration_step2_SM.cc.

References f, kBarlRings, kBarlWedges, kEndcWedgesX, kEndcWedgesY, kSides, and EEDetId::validDetId().

                                                           {
   
   TFile f("CalibHistos.root","recreate");  
 
   TH2F barreletamap("barreletamap","barreletamap",171, -85,86,100,0.,2.);
   TH2F barreletamapraw("barreletamapraw","barreletamapraw",171, -85,86,100,0.,2.);
 
   TH2F barrelmapold("barrelmapold","barrelmapold",360,1.,361.,171,-85.,86.);
   TH2F barrelmapnew("barrelmapnew","barrelmapnew",360,1.,361.,171,-85.,86.);
   TH2F barrelmapratio("barrelmapratio","barrelmapratio",360,1.,361.,171,-85.,86.);
 
   TH1F rawconst_endc_h("rawconst_endc","rawconst_endc",100,0.,2.);
   TH1F const_endc_h("const_endc","const_endc",100,0.,2.);
 
   TH1F oldconst_endc_h("oldconst_endc","oldconst_endc;oldCalib;",200,0,2);
   TH2F newvsraw_endc_h("newvsraw_endc","newvsraw_endc;rawConst;newCalib",200,0,2,200,0,2);
 
   TH2F endcapmapold_plus("endcapmapold_plus","endcapmapold_plus",100,1.,101.,100,1.,101.);
   TH2F endcapmapnew_plus("endcapmapnew_plus","endcapmapnew_plus",100,1.,101.,100,1.,101.);
   TH2F endcapmapratio_plus("endcapmapratio_plus","endcapmapratio_plus",100,1.,101.,100,1.,101.);
 
   TH2F endcapmapold_minus("endcapmapold_minus","endcapmapold_minus",100,1.,101.,100,1.,101.);
   TH2F endcapmapnew_minus("endcapmapnew_minus","endcapmapnew_minus",100,1.,101.,100,1.,101.);
   TH2F endcapmapratio_minus("endcapmapratio_minus","endcapmapratio_minus",100,1.,101.,100,1.,101.);
 
 
   for (int sign=0; sign<kSides; sign++) {
 
     int thesign = sign==1 ? 1:-1;
 
     for (int ieta=0; ieta<kBarlRings; ieta++) {
       for (int iphi=0; iphi<kBarlWedges; iphi++) {
     if(e_.goodCell_barl[ieta][iphi][sign]){
 
       EBDetId eb(thesign*( ieta+1 ), iphi+1);
       //int mod20= (iphi+1)%20;
       //if (mod20==0 || mod20==1 ||mod20==2) continue;  // exclude SM boundaries
       barreletamap.Fill(ieta*thesign + thesign,newCalibs_[eb]);
       barreletamapraw.Fill(ieta*thesign + thesign,rawconst_barl[ieta][iphi][sign]);
       
       barrelmapold.Fill(iphi+1,ieta*thesign + thesign, oldCalibs_[eb]);
       barrelmapnew.Fill(iphi+1,ieta*thesign + thesign, newCalibs_[eb]);
       barrelmapratio.Fill(iphi+1,ieta*thesign + thesign, newCalibs_[eb]/oldCalibs_[eb]);
     }//if
       }//iphi
     }//ieta
 
     for (int ix=0; ix<kEndcWedgesX; ix++) {
       for (int iy=0; iy<kEndcWedgesY; iy++) {
     if (e_.goodCell_endc[ix][iy][sign]){
       if (! EEDetId::validDetId(ix+1, iy+1,thesign)) continue;
       EEDetId ee(ix+1, iy+1,thesign);
 
       rawconst_endc_h.Fill(rawconst_endc[ix][iy][sign]);
       const_endc_h.Fill(newCalibs_[ee]);
       oldconst_endc_h.Fill(oldCalibs_[ee]);
       newvsraw_endc_h.Fill(rawconst_endc[ix][iy][sign],newCalibs_[ee]);
 
       if(sign==1){
         endcapmapold_plus.Fill(ix+1,iy+1,oldCalibs_[ee]);
         endcapmapnew_plus.Fill(ix+1,iy+1,newCalibs_[ee]);
         endcapmapratio_plus.Fill(ix+1,iy+1,newCalibs_[ee]/oldCalibs_[ee]);
       }
       else{
         endcapmapold_minus.Fill(ix+1,iy+1,oldCalibs_[ee]);
         endcapmapnew_minus.Fill(ix+1,iy+1,newCalibs_[ee]);
         endcapmapratio_minus.Fill(ix+1,iy+1,newCalibs_[ee]/oldCalibs_[ee]);
       }
 
     }//if
       }//iy
     }//ix
     
   } // sides
 
   barreletamap.Write();
   barreletamapraw.Write();
   rawconst_endc_h.Write();
   const_endc_h.Write();
   oldconst_endc_h.Write();
   newvsraw_endc_h.Write();
   barrelmapold.Write();
   barrelmapnew.Write();
   barrelmapratio.Write();
   endcapmapold_plus.Write();
   endcapmapnew_plus.Write();
   endcapmapratio_plus.Write();
   endcapmapold_minus.Write();
   endcapmapnew_minus.Write();
   endcapmapratio_minus.Write();
 
   f.Close();
 }

void PhiSymmetryCalibration_step2_SM::fillHistos ( )

Definition at line 468 of file PhiSymmetryCalibration_step2_SM.cc.

References f, kBarlRings, kBarlWedges, kEndcEtaRings, kEndcWedgesX, kEndcWedgesY, kSides, relativeConstraints::ring, and edmStreamStallGrapher::t.

 {
   TFile f("PhiSymmetryCalibration.root","recreate");
 
   std::vector<TH1F*> etsum_barl_histos(kBarlRings);
   std::vector<TH1F*> esum_barl_histos(kBarlRings);
   
   // determine ranges of ET sums to get histo bounds and book histos (barrel)
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     float low=999999.;
     float high=0.;
     float low_e=999999.;
     float high_e=0.;
 
     for (int iphi=0; iphi<kBarlWedges; iphi++) {
       for (int sign=0; sign<kSides; sign++) {
     float etsum = etsum_barl_[ieta][iphi][sign];
     if (etsum<low && etsum!=0.) low=etsum;
     if (etsum>high) high=etsum;
 
     float esum = esum_barl_[ieta][iphi][sign];
     if (esum<low_e && esum!=0.) low_e=esum;
     if (esum>high_e) high_e=esum;
       }
     }
     
     ostringstream t;
     t << "etsum_barl_" << ieta+1;
     etsum_barl_histos[ieta]=new TH1F(t.str().c_str(),"",50,low-.2*low,high+.1*high);
     t.str("");
 
     t << "esum_barl_" << ieta+1;
     esum_barl_histos[ieta]=new TH1F(t.str().c_str(),"",50,low_e-.2*low_e,high_e+.1*high_e);
     t.str("");
 
     // fill barrel ET sum histos
     etsumMean_barl_[ieta]=0.;
     esumMean_barl_[ieta]=0.;
     
     for (int iphi=0; iphi<kBarlWedges; iphi++) {
       for (int sign=0; sign<kSides; sign++) {
     
     // mean for the SC
     int iphi_r = int(iphi/nscx);
     
     if( !(iphi%nscx)){ 
      // bad channel correction
          etsum_barl_SM_[ieta][iphi_r][sign] = etsum_barl_SM_[ieta][iphi_r][sign]*nscx/(nscx- nBads_barl_SM_[ieta][iphi_r][sign]);      
 //         std::cout << "ETSUM M " << ieta << " " << iphi_r << " " << 
 //       sign << " " << etsum_barl_SM_[ieta][iphi_r][sign] << "  " << nBads_barl_SM_[ieta][iphi_r][sign]<< endl;
      etsumMean_barl_SM_[ieta] += 
         etsum_barl_SM_[ieta][iphi_r][sign]/(2*int(kBarlWedges/nscx));
     }
     if(e_.goodCell_barl[ieta][iphi][sign]){
       float etsum = etsum_barl_[ieta][iphi][sign];
       float esum  = esum_barl_[ieta][iphi][sign];
       etsum_barl_histos[ieta]->Fill(etsum);
       esum_barl_histos[ieta]->Fill(esum);
       etsumMean_barl_[ieta]+=etsum;
 
 
       esumMean_barl_[ieta]+=esum;
     }
       }
     }
 
     etsum_barl_histos[ieta]->Write();
     esum_barl_histos[ieta]->Write();
     etsumMean_barl_[ieta]/=(720.-e_.nBads_barl[ieta]);
     esumMean_barl_[ieta]/=(720.-e_.nBads_barl[ieta]);
     
     delete etsum_barl_histos[ieta];
     delete esum_barl_histos[ieta]; //VS
   }
 
   std::vector<TH1F*> etsum_endc_histos(kEndcEtaRings);
   std::vector<TH1F*> etsum_endc_uncorr_histos(kEndcEtaRings);
   std::vector<TH1F*> esum_endc_histos(kEndcEtaRings);
 
   std::vector<TH2F*> etsumvsarea_endc_histos(kEndcEtaRings);
   std::vector<TH2F*> esumvsarea_endc_histos(kEndcEtaRings);
 
   // determine ranges of ET sums to get histo bounds and book histos (endcap)
   for (int ring=0; ring<kEndcEtaRings; ring++) {
 
     float low=FLT_MAX;
     float low_uncorr=FLT_MAX;
     float high=0.;
     float high_uncorr=0;
     float low_e=FLT_MAX;
     float high_e=0.;
     float low_a=1.;
     float high_a=0.;
     for (int ix=0; ix<kEndcWedgesX; ix++) {
       for (int iy=0; iy<kEndcWedgesY; iy++) {
     if (e_.endcapRing_[ix][iy]==ring) {
       for (int sign=0; sign<kSides; sign++) {
         float etsum = etsum_endc_[ix][iy][sign];
         if (etsum<low && etsum!=0.) low=etsum;
         if (etsum>high) high=etsum;
 
         float etsum_uncorr = etsum_endc_uncorr[ix][iy][sign];
         if (etsum_uncorr<low_uncorr && etsum_uncorr!=0.) low_uncorr=etsum_uncorr;
         if (etsum_uncorr>high_uncorr) high_uncorr=etsum_uncorr;
 
         float esum = esum_endc_[ix][iy][sign];
         if (esum<low_e && esum!=0.) low_e=esum;
         if (esum>high_e) high_e=esum;
 
         float area = e_.cellArea_[ix][iy];
         if (area<low_a) low_a=area;
         if (area>high_a) high_a=area;
       }
     }
       }
     }
 
     ostringstream t;
     t<<"etsum_endc_" << ring+1;
     etsum_endc_histos[ring]= new TH1F(t.str().c_str(),"",50,low-.2*low,high+.1*high);
     t.str("");
 
     t<<"etsum_endc_uncorr_" << ring+1;
     etsum_endc_uncorr_histos[ring]= new TH1F(t.str().c_str(),"",50,low_uncorr-.2*low_uncorr,high_uncorr+.1*high_uncorr);
     t.str("");
 
     t<<"esum_endc_" << ring+1;
     esum_endc_histos[ring]= new TH1F(t.str().c_str(),"",50,low_e-.2*low_e,high_e+.1*high_e);
     t.str("");
 
     t<<"etsumvsarea_endc_" << ring+1;
     etsumvsarea_endc_histos[ring]= new TH2F(t.str().c_str(),";A_{#eta#phi};#Sigma E_{T}",50,low_a,high_a,50,low,high);
     t.str("");
 
     t<<"esumvsarea_endc_" << ring+1;
     esumvsarea_endc_histos[ring]= new TH2F(t.str().c_str(),";A_{#eta#phi};#Sigma E",50,low_a,high_a,50,low_e,high_e);
     t.str("");
 
     // fill endcap ET sum histos
     etsumMean_endc_[ring]=0.;
     esumMean_endc_[ring]=0.;
     for (int ix=0; ix<kEndcWedgesX; ix++) {
       for (int iy=0; iy<kEndcWedgesY; iy++) {
     if (e_.endcapRing_[ix][iy]==ring) {
       for (int sign=0; sign<kSides; sign++) {
         if(e_.goodCell_endc[ix][iy][sign]){
           float etsum = etsum_endc_[ix][iy][sign];
           float esum  = esum_endc_[ix][iy][sign];
           float etsum_uncorr = etsum_endc_uncorr[ix][iy][sign];
           etsum_endc_histos[ring]->Fill(etsum);
           etsum_endc_uncorr_histos[ring]->Fill(etsum_uncorr);
           esum_endc_histos[ring]->Fill(esum);
           
           float area = e_.cellArea_[ix][iy];
           etsumvsarea_endc_histos[ring]->Fill(area,etsum);
           esumvsarea_endc_histos[ring]->Fill(area,esum);
           
           etsumMean_endc_[ring]+=etsum;
           esumMean_endc_[ring]+=esum;
         }
       }
     }
       }
     }
 
     etsum_endc_histos[ring]->Write();
     etsum_endc_uncorr_histos[ring]->Write();
     esum_endc_histos[ring]->Write();
     etsumMean_endc_[ring]/=(float(e_.nRing_[ring]*2-e_.nBads_endc[ring]));
     esumMean_endc_[ring]/=(float(e_.nRing_[ring]*2-e_.nBads_endc[ring]));
     etsumvsarea_endc_histos[ring]->Write();
     esumvsarea_endc_histos[ring]->Write();
 
     delete etsum_endc_histos[ring];
     delete etsum_endc_uncorr_histos[ring];
     delete esum_endc_histos[ring];
     delete etsumvsarea_endc_histos[ring];
     delete esumvsarea_endc_histos[ring];
   }//ring
 
 
   // Maps of etsum in EB and EE
   TH2F barreletamap("barreletamap","barreletamap",171, -85,86,100,0,2);
   TH2F barrelmap("barrelmap","barrelmap - #frac{#Sigma E_{T}}{<#Sigma E_{T}>_{0}}",360,1,360, 171, -85,86);
   TH2F barrelmap_e("barrelmape","barrelmape - #frac{#Sigma E}{<#Sigma E>_{0}}",360,1,360, 171, -85,86);
   TH2F barrelmap_divided("barrelmapdiv","barrelmapdivided - #frac{#Sigma E_{T}}{hits}",360,1,360,171,-85,86);
   TH2F barrelmap_e_divided("barrelmapediv","barrelmapedivided - #frac{#Sigma E}{hits}",360,1,360,171,-85,86);
   TH2F endcmap_plus_corr("endcapmapplus_corrected","endcapmapplus - #frac{#Sigma E_{T}}{<#Sigma E_{T}>_{38}}",100,1,101,100,1,101);
   TH2F endcmap_minus_corr("endcapmapminus_corrected","endcapmapminus - #frac{#Sigma E_{T}}{<#Sigma E_{T}>_{38}}",100,1,101,100,1,101);
   TH2F endcmap_plus_uncorr("endcapmapplus_uncorrected","endcapmapplus_uncor - #frac{#Sigma E_{T}}{<#Sigma E_{T}>_{38}}",100,1,101,100,1,101);
   TH2F endcmap_minus_uncorr("endcapmapminus_uncorrected","endcapmapminus_uncor - #frac{#Sigma E_{T}}{<#Sigma E_{T}>_{38}}",100,1,101,100,1,101);
   TH2F endcmap_e_plus("endcapmapeplus","endcapmapeplus - #frac{#Sigma E}{<#Sigma E>_{38}}",100,1,101,100,1,101);
   TH2F endcmap_e_minus("endcapmapeminus","endcapmapeminus - #frac{#Sigma E}{<#Sigma E>_{38}}",100,1,101,100,1,101);
 
   for (int sign=0; sign<kSides; sign++) {
 
     int thesign = sign==1 ? 1:-1;
 
     for (int ieta=0; ieta<kBarlRings; ieta++) {
       for (int iphi=0; iphi<kBarlWedges; iphi++) {
     if(e_.goodCell_barl[ieta][iphi][sign]){
       barrelmap.Fill(iphi+1,ieta*thesign + thesign, etsum_barl_[ieta][iphi][sign]/etsumMean_barl_[0]);
       barrelmap_e.Fill(iphi+1,ieta*thesign + thesign, esum_barl_[ieta][iphi][sign]/esumMean_barl_[0]); //VS
       if (!nhits_barl_[ieta][iphi][sign]) nhits_barl_[ieta][iphi][sign] =1;
       barrelmap_divided.Fill( iphi+1,ieta*thesign + thesign, etsum_barl_[ieta][iphi][sign]/nhits_barl_[ieta][iphi][sign]);
       barrelmap_e_divided.Fill( iphi+1,ieta*thesign + thesign, esum_barl_[ieta][iphi][sign]/nhits_barl_[ieta][iphi][sign]); //VS
       //int mod20= (iphi+1)%20;
       //if (mod20==0 || mod20==1 ||mod20==2) continue;  // exclude SM boundaries
       barreletamap.Fill(ieta*thesign + thesign,etsum_barl_[ieta][iphi][sign]/etsumMean_barl_[0]);
     }//if
       }//iphi
     }//ieta
 
     for (int ix=0; ix<kEndcWedgesX; ix++) {
       for (int iy=0; iy<kEndcWedgesY; iy++) {
     if (sign==1) {
       endcmap_plus_corr.Fill(ix+1,iy+1,etsum_endc_[ix][iy][sign]/etsumMean_endc_[38]);
       endcmap_plus_uncorr.Fill(ix+1,iy+1,etsum_endc_uncorr[ix][iy][sign]/etsumMean_endc_[38]);
       endcmap_e_plus.Fill(ix+1,iy+1,esum_endc_[ix][iy][sign]/esumMean_endc_[38]);
     }
     else{ 
       endcmap_minus_corr.Fill(ix+1,iy+1,etsum_endc_[ix][iy][sign]/etsumMean_endc_[38]);
       endcmap_minus_uncorr.Fill(ix+1,iy+1,etsum_endc_uncorr[ix][iy][sign]/etsumMean_endc_[38]);
       endcmap_e_minus.Fill(ix+1,iy+1,esum_endc_[ix][iy][sign]/esumMean_endc_[38]);
     }
       }//iy
     }//ix
 
   }  //sign
   
  
 
   barreletamap.Write();
   barrelmap_divided.Write();
   barrelmap.Write();
   barrelmap_e_divided.Write();
   barrelmap_e.Write();
   endcmap_plus_corr.Write();
   endcmap_minus_corr.Write();
   endcmap_plus_uncorr.Write();
   endcmap_minus_uncorr.Write();
   endcmap_e_plus.Write();
   endcmap_e_minus.Write();
 
 
   vector<TH1F*> etavsphi_endc(kEndcEtaRings);
   vector<TH1F*> areavsphi_endc(kEndcEtaRings);
   vector<TH1F*> etsumvsphi_endcp_corr(kEndcEtaRings);
   vector<TH1F*> etsumvsphi_endcm_corr(kEndcEtaRings);
   vector<TH1F*> etsumvsphi_endcp_uncorr(kEndcEtaRings);
   vector<TH1F*> etsumvsphi_endcm_uncorr(kEndcEtaRings);
   vector<TH1F*> esumvsphi_endcp(kEndcEtaRings);
   vector<TH1F*> esumvsphi_endcm(kEndcEtaRings);
 
   std::vector<TH1F*> deltaeta_histos(kEndcEtaRings);
   std::vector<TH1F*> deltaphi_histos(kEndcEtaRings);
 
   for(int ring =0; ring<kEndcEtaRings;++ring){
     
     ostringstream t;
     t<< "etavsphi_endc_"<<ring;
     etavsphi_endc[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t<< "areavsphi_endc_"<<ring;
     areavsphi_endc[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t<< "etsumvsphi_endcp_corr_"<<ring;
     etsumvsphi_endcp_corr[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t << "etsumvsphi_endcm_corr_"<<ring;
     etsumvsphi_endcm_corr[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t << "etsumvsphi_endcp_uncorr_"<<ring;
     etsumvsphi_endcp_uncorr[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t << "etsumvsphi_endcm_uncorr_"<<ring;
     etsumvsphi_endcm_uncorr[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
     
     t << "esumvsphi_endcp_"<<ring;
     esumvsphi_endcp[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t << "esumvsphi_endcm_"<<ring;
     esumvsphi_endcm[ring] = new TH1F(t.str().c_str(), t.str().c_str(),e_.nRing_[ring],0,e_.nRing_[ring]);
     t.str("");
 
     t << "deltaeta_" << ring;
     deltaeta_histos[ring]= new TH1F(t.str().c_str(),"",50,-.1,.1);
     t.str("");
     t << "deltaphi_" << ring;
     deltaphi_histos[ring]= new TH1F(t.str().c_str(),"",50,-.1,.1);
     t.str("");
   }
 
   for (int ix=0; ix<kEndcWedgesX; ix++) {
     for (int iy=0; iy<kEndcWedgesY; iy++) {
 
       int ring = e_.endcapRing_[ix][iy];
       if (ring!=-1) {
     int iphi_endc=-1;
     for (int ip=0; ip<e_.nRing_[ring]; ip++) {
       if (e_.cellPhi_[ix][iy]==e_.phi_endc_[ip][ring]) iphi_endc=ip;
     }
 
     if(iphi_endc!=-1){
       for (int sign=0; sign<kSides; sign++) {
         if(e_.goodCell_endc[ix][iy][sign]){
           if (sign==1){
         etsumvsphi_endcp_corr[ring]->Fill(iphi_endc,etsum_endc_[ix][iy][sign]);
         etsumvsphi_endcp_uncorr[ring]->Fill(iphi_endc,etsum_endc_uncorr[ix][iy][sign]);
         esumvsphi_endcp[ring]->Fill(iphi_endc,esum_endc_[ix][iy][sign]);
           } else {
         etsumvsphi_endcm_corr[ring]->Fill(iphi_endc,etsum_endc_[ix][iy][sign]);
         etsumvsphi_endcm_uncorr[ring]->Fill(iphi_endc,etsum_endc_uncorr[ix][iy][sign]);
         esumvsphi_endcm[ring]->Fill(iphi_endc,esum_endc_[ix][iy][sign]);
           }
         }//if
       }//sign
       etavsphi_endc[ring]->Fill(iphi_endc,e_.cellPos_[ix][iy].eta());
       areavsphi_endc[ring]->Fill(iphi_endc,e_.cellArea_[ix][iy]);
     } //if iphi_endc
 
       }//if ring
     }//iy
   } //ix
 
 
 
   for(int ring =0; ring<kEndcEtaRings;++ring){
     
     etavsphi_endc[ring]->Write();
     areavsphi_endc[ring]->Write();
     etsumvsphi_endcp_corr[ring]->Write();
     etsumvsphi_endcm_corr[ring]->Write();
     etsumvsphi_endcp_uncorr[ring]->Write();
     etsumvsphi_endcm_uncorr[ring]->Write();
     esumvsphi_endcp[ring]->Write();
     esumvsphi_endcm[ring]->Write();
     deltaeta_histos[ring]->Write();
     deltaphi_histos[ring]->Write();
 
     
     delete etsumvsphi_endcp_corr[ring];
     delete etsumvsphi_endcm_corr[ring];
     delete etsumvsphi_endcp_uncorr[ring];
     delete etsumvsphi_endcm_uncorr[ring];
     delete etavsphi_endc[ring];
     delete areavsphi_endc[ring];
     delete esumvsphi_endcp[ring];
     delete esumvsphi_endcm[ring];
     delete deltaeta_histos[ring];
     delete deltaphi_histos[ring];
   }
 
 
   f.Close();
 }

void PhiSymmetryCalibration_step2_SM::outResidHistos ( )

Definition at line 912 of file PhiSymmetryCalibration_step2_SM.cc.

References f, and relativeConstraints::ring.

                                                      {
 
   // output histograms of residual miscalibrations
   TFile f("PhiSymmetryCalibration_miscal_resid.root","recreate");
   for (int ieta=0; ieta<85; ieta++) {
     miscal_resid_barl_histos[ieta]->Write();
     correl_barl_histos[ieta]->Write();
 
     delete miscal_resid_barl_histos[ieta];
     delete correl_barl_histos[ieta];
   }
   
   for (int ring=0; ring<39; ring++) {
     miscal_resid_endc_histos[ring]->Write();
     correl_endc_histos[ring]->Write();
     
     delete  miscal_resid_endc_histos[ring];
     delete  correl_endc_histos[ring];
 
   }
   f.Close(); 
 }

void PhiSymmetryCalibration_step2_SM::readEtSums ( )

Definition at line 833 of file PhiSymmetryCalibration_step2_SM.cc.

References recoMuon::in, kBarlRings, kEndcEtaRings, and relativeConstraints::ring.

                                                 {
 
 
   //read in ET sums
   
   int ieta,iphi,sign,ix,iy,dummy;
   double etsum;
   unsigned int nhits;
   std::ifstream etsum_barl_in("etsum_barl.dat", ios::in);
   while ( etsum_barl_in >> dummy >> ieta >> iphi >> sign >> etsum >> nhits ) {
     etsum_barl_[ieta][iphi][sign]+=etsum;
     nhits_barl_[ieta][iphi][sign]+=nhits;
     
     // fill etsums for the SM calibration
     int iphi_r = int(iphi/nscx);
     etsum_barl_SM_[ieta][iphi_r][sign]+= etsum;
 //      etsum*nscx/(nscx- nBads_barl_SM_[ieta][iphi_r][sign]);
  //     if(nBads_barl_SM_[ieta][iphi_r][sign]>0){
  //       std::cout << "ETSUM" << etsum_barl_SM_[ieta][iphi_r][sign] << "  " << nscx << "  " << nBads_barl_SM_[ieta][iphi_r][sign]<< endl;
  //     }      
   }
  
   std::ifstream etsum_endc_in("etsum_endc.dat", ios::in);
   while ( etsum_endc_in >> dummy >> ix >> iy >> sign >> etsum >> nhits>>dummy ) {
     etsum_endc_[ix][iy][sign]+=etsum;
     nhits_endc_[ix][iy][sign]+=nhits;
 
   }
  
   std::ifstream k_barl_in("k_barl.dat", ios::in);
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     k_barl_in >> dummy >> k_barl_[ieta];
   }
  
   std::ifstream k_endc_in("k_endc.dat", ios::in);
   for (int ring=0; ring<kEndcEtaRings; ring++) {
     k_endc_in >> dummy >> k_endc_[ring];
   }
   
 
 }

void PhiSymmetryCalibration_step2_SM::setUp ( const edm::EventSetup & setup )

if a miscalibration was applied, load it, if not put it to 1

Definition at line 51 of file PhiSymmetryCalibration_step2_SM.cc.

References DetId::Ecal, EcalBarrel, EcalEndcap, edm::EventSetup::get(), h, kBarlRings, kBarlWedges, kSides, AlCaHLTBitMon_ParallelJobs::p, cmsHarvester::path, EcalFloatCondObjectContainerXMLTranslator::readXML(), and run_regression::ret.

                                                                   {
 
   edm::ESHandle<EcalChannelStatus> chStatus;
   se.get<EcalChannelStatusRcd>().get(chStatus);
 
   edm::ESHandle<CaloGeometry> geoHandle;
   se.get<CaloGeometryRecord>().get(geoHandle);
 
   barrelCells = geoHandle->getValidDetIds(DetId::Ecal, EcalBarrel);
   endcapCells = geoHandle->getValidDetIds(DetId::Ecal, EcalEndcap);
 
   e_.setup(&(*geoHandle), &(*chStatus), statusThreshold_);
 
 
 
   for (int sign=0; sign<kSides; sign++) {
     for (int ieta=0; ieta<kBarlRings; ieta++) {
       for (int iphi=0; iphi<kBarlWedges; iphi++) {
     int iphi_r=int(iphi/nscx);
     if( !e_.goodCell_barl[ieta][iphi][sign] ){
       nBads_barl_SM_[ieta][iphi_r][sign]++;
       // std::cout << "N BAD CELL " << nBads_barl_SM_[ieta][iphi_r][sign] << endl; 
          }
       }
     }
   }
   
 
 
   if (have_initial_miscalib_){
 
     EcalCondHeader h;
     namespace fs = boost::filesystem;
     fs::path p(initialmiscalibfile_.c_str());
     if (!fs::exists(p)) edm::LogError("PhiSym") << "File not found: " 
                         << initialmiscalibfile_ <<endl;  
     
     int ret=
       EcalIntercalibConstantsXMLTranslator::readXML(initialmiscalibfile_,h,miscalib_);    
     if (ret) edm::LogError("PhiSym")<<"Error reading XML files"<<endl;
   } else {
 
     for (vector<DetId>::iterator it=barrelCells.begin(); it!=barrelCells.end(); ++it){
       miscalib_[*it]=1;
      }
 
     for (vector<DetId>::iterator it=endcapCells.begin(); it!=endcapCells.end(); ++it){
       miscalib_[*it]=1;
  
     }
   }
     
   // if we are reiterating, read constants from previous iter
   // if not put them to one
   if (reiteration_){
     
      
     EcalCondHeader h;
     namespace fs = boost::filesystem;
     fs::path p(oldcalibfile_.c_str());
     if (!fs::exists(p)) edm::LogError("PhiSym") << "File not found: " 
                         << oldcalibfile_ <<endl;  
     
     int ret=
       EcalIntercalibConstantsXMLTranslator::readXML(oldcalibfile_,h,
                             oldCalibs_); 
    
     if (ret) edm::LogError("PhiSym")<<"Error reading XML files"<<endl;;
     
   } else {
     
     for (vector<DetId>::iterator it=barrelCells.begin(); 
      it!=barrelCells.end(); ++it)
       oldCalibs_[*it]=1;
      
 
     for (vector<DetId>::iterator it=endcapCells.begin(); 
      it!=endcapCells.end(); ++it)
       oldCalibs_[*it]=1;
       
    
   } // else 
   
 }

void PhiSymmetryCalibration_step2_SM::setupResidHistos ( )

Definition at line 877 of file PhiSymmetryCalibration_step2_SM.cc.

References kBarlRings, kEndcEtaRings, relativeConstraints::ring, and reco::t2.

                                                       {
 
   miscal_resid_barl_histos.resize(kBarlRings);
   correl_barl_histos.resize(kBarlRings);  
  
 
   for (int ieta=0; ieta<kBarlRings; ieta++) {
     ostringstream t1; 
     t1<<"mr_barl_"<<ieta+1;
     miscal_resid_barl_histos[ieta] = new TH1F(t1.str().c_str(),"",100,0.,2.);
     ostringstream t2;
     t2<<"co_barl_"<<ieta+1;
     correl_barl_histos[ieta] = new TH2F(t2.str().c_str(),"",50,.5,1.5,50,.5,1.5);
   }
 
   miscal_resid_endc_histos.resize(kEndcEtaRings);
   correl_endc_histos.resize(kEndcEtaRings);
 
 
   for (int ring=0; ring<kEndcEtaRings; ring++) {
     ostringstream t1;
     t1<<"mr_endc_"<< ring+1;
     miscal_resid_endc_histos[ring] = new TH1F(t1.str().c_str(),"",100,0.,2.);
     ostringstream t2;
     t2<<"co_endc_"<<ring+1;
     correl_endc_histos[ring] = new TH2F(t2.str().c_str(),"",50,.5,1.5,50,.5,1.5);
   }
 
  
 
   
 
 }