Pi0FixedMassWindowCalibration.cc

Go to the documentation of this file.

#include "Calibration/EcalCalibAlgos/interface/Pi0FixedMassWindowCalibration.h"

// System include files

// Framework


// Conditions database

#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
#include "Calibration/Tools/interface/Pi0CalibXMLwriter.h"

// Reconstruction Classes
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
// Geometry
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
#include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"

// EgammaCoreTools
#include "DataFormats/EgammaReco/interface/ClusterShape.h"
#include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"

//const double Pi0Calibration::PDGPi0Mass =  0.1349766;

using namespace std;

//_____________________________________________________________________________

Pi0FixedMassWindowCalibration::Pi0FixedMassWindowCalibration(const edm::ParameterSet& iConfig) :
  theMaxLoops( iConfig.getUntrackedParameter<unsigned int>("maxLoops",0) ),
  ecalHitsProducer_( iConfig.getParameter< std::string > ("ecalRecHitsProducer") ),
  barrelHits_( iConfig.getParameter< std::string > ("barrelHitCollection") )
{

  std::cout << "[Pi0FixedMassWindowCalibration] Constructor "<<std::endl;
  // The verbosity level
  std::string verbosityString = iConfig.getParameter<std::string>("VerbosityLevel");
  if      (verbosityString == "DEBUG")   verbosity = IslandClusterAlgo::pDEBUG;
  else if (verbosityString == "WARNING") verbosity = IslandClusterAlgo::pWARNING;
  else if (verbosityString == "INFO")    verbosity = IslandClusterAlgo::pINFO;
  else                                   verbosity = IslandClusterAlgo::pERROR;

  // The names of the produced cluster collections
  barrelClusterCollection_  = iConfig.getParameter<std::string>("barrelClusterCollection");

  // Island algorithm parameters
  double barrelSeedThreshold = iConfig.getParameter<double>("IslandBarrelSeedThr");
  double endcapSeedThreshold = iConfig.getParameter<double>("IslandEndcapSeedThr");


  // Selection algorithm parameters
  selePi0PtGammaOneMin_ = iConfig.getParameter<double>("selePi0PtGammaOneMin");
  selePi0PtGammaTwoMin_ = iConfig.getParameter<double>("selePi0PtGammaTwoMin");

  selePi0DRBelt_ = iConfig.getParameter<double>("selePi0DRBelt");
  selePi0DetaBelt_ = iConfig.getParameter<double>("selePi0DetaBelt");

  selePi0PtPi0Min_ = iConfig.getParameter<double>("selePi0PtPi0Min");

  selePi0S4S9GammaOneMin_ = iConfig.getParameter<double>("selePi0S4S9GammaOneMin");
  selePi0S4S9GammaTwoMin_ = iConfig.getParameter<double>("selePi0S4S9GammaTwoMin");
  selePi0S9S25GammaOneMin_ = iConfig.getParameter<double>("selePi0S9S25GammaOneMin");
  selePi0S9S25GammaTwoMin_ = iConfig.getParameter<double>("selePi0S9S25GammaTwoMin");
 
  selePi0EtBeltIsoRatioMax_ = iConfig.getParameter<double>("selePi0EtBeltIsoRatioMax");

  selePi0MinvMeanFixed_ = iConfig.getParameter<double>("selePi0MinvMeanFixed");
  selePi0MinvSigmaFixed_ = iConfig.getParameter<double>("selePi0MinvSigmaFixed");




  // Parameters for the position calculation:
  edm::ParameterSet posCalcParameters = 
    iConfig.getParameter<edm::ParameterSet>("posCalcParameters");
  posCalculator_ = PositionCalc(posCalcParameters); 
  shapeAlgo_ = ClusterShapeAlgo(posCalcParameters);
  clustershapecollectionEB_ = iConfig.getParameter<std::string>("clustershapecollectionEB");

  //AssociationMap
  barrelClusterShapeAssociation_ = iConfig.getParameter<std::string>("barrelShapeAssociation");

  island_p = new IslandClusterAlgo(barrelSeedThreshold, endcapSeedThreshold, posCalculator_,verbosity);

  theParameterSet=iConfig;


}


//_____________________________________________________________________________
// Close files, etc.

Pi0FixedMassWindowCalibration::~Pi0FixedMassWindowCalibration()
{

}

//_____________________________________________________________________________
// Initialize algorithm

void Pi0FixedMassWindowCalibration::beginOfJob( )
{

  //std::cout << "[Pi0FixedMassWindowCalibration] beginOfJob "<<std::endl;

  isfirstcall_=true;
  
 

}


void Pi0FixedMassWindowCalibration::endOfJob()
{
  
  std::cout << "[Pi0FixedMassWindowCalibration] endOfJob"<<endl;

  // Write new calibration constants

  Pi0CalibXMLwriter barrelWriter(EcalBarrel,99);

  std::vector<DetId>::const_iterator barrelIt=barrelCells.begin();
  for (; barrelIt!=barrelCells.end(); barrelIt++) {
    EBDetId eb(*barrelIt);
    int ieta = eb.ieta();
    int iphi = eb.iphi();
    int sign = eb.zside()>0 ? 1 : 0;
    barrelWriter.writeLine(eb,newCalibs_barl[abs(ieta)-1][iphi-1][sign]);
  }

}

//_____________________________________________________________________________
// Called at beginning of loop
void Pi0FixedMassWindowCalibration::startingNewLoop(unsigned int iLoop )
{

  for (int sign=0; sign<2; sign++) {
    for (int ieta=0; ieta<85; ieta++) {
      for (int iphi=0; iphi<360; iphi++) {
    wxtals[ieta][iphi][sign]=0.;
    mwxtals[ieta][iphi][sign]=0.;
      }
    }
  }
  std::cout << "[Pi0FixedMassWindowCalibration] Starting loop number " << iLoop<<std::endl;

}


//_____________________________________________________________________________
// Called at end of loop

edm::EDLooper::Status 
Pi0FixedMassWindowCalibration::endOfLoop(const edm::EventSetup& iSetup, unsigned int iLoop)
{
  std::cout << "[Pi0FixedMassWindowCalibration] Ending loop " << iLoop<<std::endl;

  for (int sign=0; sign<2; sign++) {
    for (int ieta=0; ieta<85; ieta++) {
      for (int iphi=0; iphi<360; iphi++) {

        if (wxtals[ieta][iphi][sign] == 0 ) 
      {
        newCalibs_barl[ieta][iphi][sign]=oldCalibs_barl[ieta][iphi][sign];
      } else {
        newCalibs_barl[ieta][iphi][sign]=oldCalibs_barl[ieta][iphi][sign]*(mwxtals[ieta][iphi][sign]/wxtals[ieta][iphi][sign]);
      }
    cout<< " New calibration constant: ieta iphi sign - old,mwxtals,wxtals,new: "<<ieta<<" "<<iphi<<" "<<sign<<" - "<<oldCalibs_barl[ieta][iphi][sign]<<" "<<mwxtals[ieta][iphi][sign]<<" "<<wxtals[ieta][iphi][sign]<<" "<<newCalibs_barl[ieta][iphi][sign]<<endl;

      }
    }
  }
    
  Pi0CalibXMLwriter barrelWriter(EcalBarrel,iLoop+1);

  std::vector<DetId>::const_iterator barrelIt=barrelCells.begin();
  for (; barrelIt!=barrelCells.end(); barrelIt++) {
    EBDetId eb(*barrelIt);
    int ieta = eb.ieta();
    int iphi = eb.iphi();
    int sign = eb.zside()>0 ? 1 : 0;
    barrelWriter.writeLine(eb,newCalibs_barl[abs(ieta)-1][iphi-1][sign]);
    if (iphi==1) {
      std::cout << "Calib constant for barrel crystal "
        << " (" << ieta << "," << iphi << ") changed from "
        << oldCalibs_barl[abs(ieta)-1][iphi-1][sign] << " to "
        << newCalibs_barl[abs(ieta)-1][iphi-1][sign] << std::endl;
    }
  }

  // replace old calibration constants with new one

  for (int sign=0; sign<2; sign++) {
    for (int ieta=0; ieta<85; ieta++) {
      for (int iphi=0; iphi<360; iphi++) {
        oldCalibs_barl[ieta][iphi][sign]=newCalibs_barl[ieta][iphi][sign];
    newCalibs_barl[ieta][iphi][sign]=0;
      }
    }
  }


  if ( iLoop == theMaxLoops-1 || iLoop >= theMaxLoops ) return kStop;
  else return kContinue;
}

//_____________________________________________________________________________
// Called at each event

edm::EDLooper::Status 
Pi0FixedMassWindowCalibration::duringLoop(const edm::Event& event, 
  const edm::EventSetup& setup )
{
  using namespace edm;
  using namespace std;

  // this chunk used to belong to beginJob(isetup). Moved here
  // with the beginJob without arguments migration
  
  if (isfirstcall_){
   // initialize arrays

    for (int sign=0; sign<2; sign++) {
      for (int ieta=0; ieta<85; ieta++) {
    for (int iphi=0; iphi<360; iphi++) {
      oldCalibs_barl[ieta][iphi][sign]=0.;
      newCalibs_barl[ieta][iphi][sign]=0.;
      wxtals[ieta][iphi][sign]=0.;
      mwxtals[ieta][iphi][sign]=0.;
    }
      }
    }
    
    // get initial constants out of DB
    
    edm::ESHandle<EcalIntercalibConstants> pIcal;
    EcalIntercalibConstantMap imap;
    
    try {
      setup.get<EcalIntercalibConstantsRcd>().get(pIcal);
      std::cout << "Taken EcalIntercalibConstants" << std::endl;
      imap = pIcal.product()->getMap();
      std::cout << "imap.size() = " << imap.size() << std::endl;
    } catch ( std::exception& ex ) {     
      std::cerr << "Error! can't get EcalIntercalibConstants " << std::endl;
    } 
    
    // get the ecal geometry:
    edm::ESHandle<CaloGeometry> geoHandle;
    setup.get<CaloGeometryRecord>().get(geoHandle);
    const CaloGeometry& geometry = *geoHandle;
    //const CaloSubdetectorGeometry *barrelGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
    
    // loop over all barrel crystals
    barrelCells = geometry.getValidDetIds(DetId::Ecal, EcalBarrel);
    std::vector<DetId>::const_iterator barrelIt;
    for (barrelIt=barrelCells.begin(); barrelIt!=barrelCells.end(); barrelIt++) {
      EBDetId eb(*barrelIt);
      
      // get the initial calibration constants
      EcalIntercalibConstantMap::const_iterator itcalib = imap.find(eb.rawId());
      if ( itcalib == imap.end() ) {
    // FIXME -- throw error
      }
      EcalIntercalibConstant calib = (*itcalib);
      int sign = eb.zside()>0 ? 1 : 0;
      oldCalibs_barl[abs(eb.ieta())-1][eb.iphi()-1][sign] = calib;
      if (eb.iphi()==1) std::cout << "Read old constant for crystal "
                  << " (" << eb.ieta() << "," << eb.iphi()
                  << ") : " << calib << std::endl;
      
      
      
    }
    isfirstcall_=false;
  }
  


  nevent++;

  if ((nevent<100 && nevent%10==0) 
      ||(nevent<1000 && nevent%100==0) 
      ||(nevent<10000 && nevent%100==0) 
      ||(nevent<100000 && nevent%1000==0) 
      ||(nevent<10000000 && nevent%1000==0))
    std::cout << "[Pi0FixedMassWindowCalibration] Events processed: "<<nevent<<std::endl;

  recHitsEB_map = new std::map<DetId, EcalRecHit>();

  EcalRecHitCollection* recalibEcalRecHitCollection( new EcalRecHitCollection );

  int nRecHitsEB=0;
  Handle<EcalRecHitCollection> pEcalRecHitBarrelCollection;
  event.getByLabel(ecalHitsProducer_, barrelHits_, pEcalRecHitBarrelCollection);
  const EcalRecHitCollection* ecalRecHitBarrelCollection = pEcalRecHitBarrelCollection.product();
  cout << " ECAL Barrel RecHits # "<< ecalRecHitBarrelCollection->size() <<endl;
  for(EcalRecHitCollection::const_iterator aRecHitEB = ecalRecHitBarrelCollection->begin(); aRecHitEB != ecalRecHitBarrelCollection->end(); aRecHitEB++) {
    //cout << " ECAL Barrel RecHit #,E,time,det,subdetid: "<<nRecHitsEB<<" "<<aRecHitEB->energy()<<" "<<aRecHitEB->time()<<" "<<aRecHitEB->detid().det()<<" "<<aRecHitEB->detid().subdetId()<<endl;


    EBDetId ebrhdetid = aRecHitEB->detid();
    //cout << " EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
    //cout << " EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
    //cout << " EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;

    int sign = ebrhdetid.zside()>0 ? 1 : 0;
    EcalRecHit aHit(aRecHitEB->id(),aRecHitEB->energy()*oldCalibs_barl[abs(ebrhdetid.ieta())-1][ebrhdetid.iphi()-1][sign],aRecHitEB->time());
    recalibEcalRecHitCollection->push_back(aHit);

    nRecHitsEB++;

  }

  //  cout<<" Recalib size: "<<recalibEcalRecHitCollection->size()<<endl;
  int irecalib=0;
  for(EcalRecHitCollection::const_iterator aRecHitEB = recalibEcalRecHitCollection->begin(); aRecHitEB != recalibEcalRecHitCollection->end(); aRecHitEB++) {
    //cout << " [recalibrated] ECAL Barrel RecHit #,E,time,det,subdetid: "<<irecalib<<" "<<aRecHitEB->energy()<<" "<<aRecHitEB->time()<<" "<<aRecHitEB->detid().det()<<" "<<aRecHitEB->detid().subdetId()<<endl;

    //    EBDetId ebrhdetid = aRecHitEB->detid();
    //cout << " [recalibrated] EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
    //cout << " [recalibrated] EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
    //cout << " [recalibrated] EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;

    std::pair<DetId, EcalRecHit> map_entry(aRecHitEB->id(), *aRecHitEB);
    recHitsEB_map->insert(map_entry);


    irecalib++;

  }

  
  // get the geometry and topology from the event setup:
  edm::ESHandle<CaloGeometry> geoHandle;
  setup.get<CaloGeometryRecord>().get(geoHandle);

  const CaloSubdetectorGeometry *geometry_p;
  CaloSubdetectorTopology *topology_p;

  std::string clustershapetag;
  geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
  topology_p = new EcalBarrelTopology(geoHandle);

  const CaloSubdetectorGeometry *geometryES_p;
  geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);

  /*
  reco::BasicClusterCollection clusters;
  clusters = island_p->makeClusters(ecalRecHitBarrelCollection,geometry_p,topology_p,geometryES_p,IslandClusterAlgo::barrel);
  
  //Create associated ClusterShape objects.
  std::vector <reco::ClusterShape> ClusVec;
  for (int erg=0;erg<int(clusters.size());++erg){
    reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg],ecalRecHitBarrelCollection,geometry_p,topology_p);
    ClusVec.push_back(TestShape);
  }

  //Put clustershapes in event, but retain a Handle on them.
  std::auto_ptr<reco::ClusterShapeCollection> clustersshapes_p(new reco::ClusterShapeCollection);
  clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());

  cout<<"[Pi0Calibration] Basic Cluster collection size: "<<clusters.size()<<endl;
  cout<<"[Pi0Calibration] Basic Cluster Shape Collection size: "<<clustersshapes_p->size()<<endl;

  int iClus=0;
  for(reco::BasicClusterCollection::const_iterator aClus = clusters.begin(); aClus != clusters.end(); aClus++) {
    cout<<" CLUSTER : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: "<<iClus<<" "<<aClus->getHitsByDetId().size()<<" "<<aClus->energy()<<" "<<aClus->energy()*sin(aClus->position().theta())<<" "<<aClus->position().eta()<<" "<<aClus->position().phi()<<" "<<(*clustersshapes_p)[iClus].e2x2()<<" "<<(*clustersshapes_p)[iClus].e3x3()<<" "<<(*clustersshapes_p)[iClus].e5x5()<<endl; 
    iClus++;
  }
  */

  // recalibrated clusters
  reco::BasicClusterCollection clusters_recalib;
  clusters_recalib = island_p->makeClusters(recalibEcalRecHitCollection,geometry_p,topology_p,geometryES_p,IslandClusterAlgo::barrel);
  
  //Create associated ClusterShape objects.
  std::vector <reco::ClusterShape> ClusVec_recalib;
  for (int erg=0;erg<int(clusters_recalib.size());++erg){
    reco::ClusterShape TestShape_recalib = shapeAlgo_.Calculate(clusters_recalib[erg],recalibEcalRecHitCollection,geometry_p,topology_p);
    ClusVec_recalib.push_back(TestShape_recalib);
  }

  //Put clustershapes in event, but retain a Handle on them.
  std::auto_ptr<reco::ClusterShapeCollection> clustersshapes_p_recalib(new reco::ClusterShapeCollection);
  clustersshapes_p_recalib->assign(ClusVec_recalib.begin(), ClusVec_recalib.end());

  cout<<"[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster collection size: "<<clusters_recalib.size()<<endl;
  cout<<"[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster Shape Collection size: "<<clustersshapes_p_recalib->size()<<endl;


  // pizero selection

  // Get ECAL Barrel Island Basic Clusters collection
  // ECAL Barrel Island Basic Clusters 
  static const int MAXBCEB = 200;
  static const int MAXBCEBRH = 200;
  int nIslandBCEB;
  float eIslandBCEB[MAXBCEB];
  float etIslandBCEB[MAXBCEB];
  float etaIslandBCEB[MAXBCEB];
  float phiIslandBCEB[MAXBCEB];
  float e2x2IslandBCEB[MAXBCEB];
  float e3x3IslandBCEB[MAXBCEB];
  float e5x5IslandBCEB[MAXBCEB];
  // indexes to the RecHits assiciated with 
  // ECAL Barrel Island Basic Clusters
  int nIslandBCEBRecHits[MAXBCEB];
  //  int indexIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
  int ietaIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
  int iphiIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
  int zsideIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
  float eIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];

  nIslandBCEB=0;
  for(int i=0; i<MAXBCEB; i++){
    eIslandBCEB[i] = 0;
    etIslandBCEB[i] = 0;
    etaIslandBCEB[i] = 0;
    phiIslandBCEB[i] = 0;
    e2x2IslandBCEB[i] = 0;
    e3x3IslandBCEB[i] = 0;
    e5x5IslandBCEB[i] = 0;
     nIslandBCEBRecHits[i] = 0;
     for(int j=0;j<MAXBCEBRH;j++){
       //       indexIslandBCEBRecHits[i][j] = 0;
       ietaIslandBCEBRecHits[i][j] = 0;
       iphiIslandBCEBRecHits[i][j] = 0;
       zsideIslandBCEBRecHits[i][j] = 0;
       eIslandBCEBRecHits[i][j] = 0;
     }
  }


  int iClus_recalib=0;
  for(reco::BasicClusterCollection::const_iterator aClus = clusters_recalib.begin(); aClus != clusters_recalib.end(); aClus++) {
    cout<<" CLUSTER [recalibration] : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: "<<iClus_recalib<<" "<<aClus->size()<<" "<<aClus->energy()<<" "<<aClus->energy()*sin(aClus->position().theta())<<" "<<aClus->position().eta()<<" "<<aClus->position().phi()<<" "<<(*clustersshapes_p_recalib)[iClus_recalib].e2x2()<<" "<<(*clustersshapes_p_recalib)[iClus_recalib].e3x3()<<" "<<(*clustersshapes_p_recalib)[iClus_recalib].e5x5()<<endl; 

    eIslandBCEB[nIslandBCEB] = aClus->energy();
    etIslandBCEB[nIslandBCEB] = aClus->energy()*sin(aClus->position().theta());
    etaIslandBCEB[nIslandBCEB] = aClus->position().eta();
    phiIslandBCEB[nIslandBCEB] = aClus->position().phi();
    
    e2x2IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e2x2();    
    e3x3IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e3x3();    
    e5x5IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e5x5();    

    nIslandBCEBRecHits[nIslandBCEB] = aClus->size();

    std::vector<std::pair< DetId,float> > hits = aClus->hitsAndFractions();
    std::vector<std::pair< DetId,float> >::iterator hit;
    std::map<DetId, EcalRecHit>::iterator aHit;

    int irhcount=0;
    for(hit = hits.begin(); hit != hits.end(); hit++)
      {
        // need to get hit by DetID in order to get energy
        aHit = recHitsEB_map->find((*hit).first);
        //cout << "       RecHit #: "<<irhcount<<" from Basic Cluster with Energy: "<<aHit->second.energy()<<endl; 

        EBDetId sel_rh = aHit->second.detid();
        //cout << "       RecHit: z,ieta,iphi "<<sel_rh.zside()<<" "<<sel_rh.ieta()<<" "<<sel_rh.iphi()<<endl;
        //cout << "       RecHit: tower_ieta,tower_iphi "<<sel_rh.tower_ieta()<<" "<<sel_rh.tower_iphi()<<endl;
        //cout << "       RecHit: iSM, ic "<<sel_rh.ism()<<" "<<sel_rh.ic()<<endl;

    ietaIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.ieta();
    iphiIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.iphi();
    zsideIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.zside();
    eIslandBCEBRecHits[nIslandBCEB][irhcount] = aHit->second.energy();

        irhcount++;
      }
    nIslandBCEB++;
    iClus_recalib++;


  }

  // Selection, based on ECAL Barrel Basic Clusters

  if (nIslandBCEB > 1) 
    {
      for(int i=0 ; i<nIslandBCEB ; i++)
        {
          for(int j=i+1 ; j<nIslandBCEB ; j++)
            {

              if( etIslandBCEB[i]>selePi0PtGammaOneMin_ && etIslandBCEB[j]>selePi0PtGammaOneMin_) 
        {
                
          float theta_0 = 2. * atan(exp(-etaIslandBCEB[i]));
          float theta_1 = 2. * atan(exp(-etaIslandBCEB[j]));
                
          float p0x = eIslandBCEB[i] * sin(theta_0) * cos(phiIslandBCEB[i]);
          float p1x = eIslandBCEB[j] * sin(theta_1) * cos(phiIslandBCEB[j]);
                
          float p0y = eIslandBCEB[i] * sin(theta_0) * sin(phiIslandBCEB[i]);
          float p1y = eIslandBCEB[j] * sin(theta_1) * sin(phiIslandBCEB[j]);
                
          float p0z = eIslandBCEB[i] * cos(theta_0);
          float p1z = eIslandBCEB[j] * cos(theta_1);

          float pi0_px = p0x + p1x;
          float pi0_py = p0y + p1y;
          float pi0_pz = p0z + p1z;
                
          float pi0_ptot = sqrt (pi0_px*pi0_px + pi0_py*pi0_py + pi0_pz*pi0_pz); 

          float pi0_theta = acos(pi0_pz/pi0_ptot);
          float pi0_eta = -log(tan(pi0_theta/2));
          float pi0_phi = atan(pi0_py/pi0_px);
          //cout << " pi0_theta, pi0_eta, pi0_phi "<<pi0_theta<<" "<<pi0_eta<<" "<<pi0_phi<<endl;

          // belt isolation
                
          float et_belt=0;
          for(Int_t k=0 ; k<nIslandBCEB ; k++)
            {
              if( (k != i) && (k != j) ) 
            {
              float dr_pi0_k = sqrt( (etaIslandBCEB[k]-pi0_eta)*(etaIslandBCEB[k]-pi0_eta) + (phiIslandBCEB[k]-pi0_phi)*(phiIslandBCEB[k]-pi0_phi) );
              float deta_pi0_k = fabs(etaIslandBCEB[k]-pi0_eta);
              if ( (dr_pi0_k<selePi0DRBelt_) && (deta_pi0_k<selePi0DetaBelt_) ) et_belt = et_belt + etIslandBCEB[k];
            }
            }

                
          float pt_pi0 = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
          //float dr_pi0 = sqrt ( (etaIslandBCEB[i]-etaIslandBCEB[j])*(etaIslandBCEB[i]-etaIslandBCEB[j]) + (phiIslandBCEB[i]-phiIslandBCEB[j])*(phiIslandBCEB[i]-phiIslandBCEB[j]) );

          //cout <<" pi0 pt,dr:  "<<pt_pi0<<" "<<dr_pi0<<endl; 
          if (pt_pi0 > selePi0PtPi0Min_) 
            {
              float m_inv = sqrt ( (eIslandBCEB[i] + eIslandBCEB[j])*(eIslandBCEB[i] + eIslandBCEB[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
              cout <<" pi0 (pt>2.5 GeV) m_inv = "<<m_inv<<endl;  

              float s4s9_1 = e2x2IslandBCEB[i]/e3x3IslandBCEB[i];
              float s4s9_2 = e2x2IslandBCEB[j]/e3x3IslandBCEB[j];

              float s9s25_1 = e3x3IslandBCEB[i]/e5x5IslandBCEB[i];
              float s9s25_2 = e3x3IslandBCEB[j]/e5x5IslandBCEB[j];

              //float s9Esc_1 = e3x3IslandBCEB[i]/eIslandBCEB[i];
              //float s9Esc_2 = e3x3IslandBCEB[j]/eIslandBCEB[j];

              if (s4s9_1>selePi0S4S9GammaOneMin_ && s4s9_2>selePi0S4S9GammaTwoMin_ && s9s25_1>selePi0S9S25GammaOneMin_ && s9s25_2>selePi0S9S25GammaTwoMin_ && (et_belt/pt_pi0 < selePi0EtBeltIsoRatioMax_) ) 
            {
              //good pizero candidate
              if ( m_inv>(selePi0MinvMeanFixed_ - 2*selePi0MinvSigmaFixed_) && m_inv<(selePi0MinvMeanFixed_ + 2*selePi0MinvSigmaFixed_) ) 
                {
                  //fill wxtals and mwxtals weights
                  cout<<" Pi0 Good candidate : minv = "<<m_inv<<endl;
                  for(int kk=0 ; kk<nIslandBCEBRecHits[i] ; kk++)
                {
                  int ieta_xtal = ietaIslandBCEBRecHits[i][kk];
                  int iphi_xtal = iphiIslandBCEBRecHits[i][kk];
                  int sign = zsideIslandBCEBRecHits[i][kk]>0 ? 1 : 0;
                  wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + eIslandBCEBRecHits[i][kk]/e3x3IslandBCEB[i];
                  mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + (selePi0MinvMeanFixed_/m_inv)*(selePi0MinvMeanFixed_/m_inv)*(eIslandBCEBRecHits[i][kk]/e3x3IslandBCEB[i]);
                  cout<< "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: "<<ieta_xtal<<" "<<iphi_xtal<<" "<<sign<<" "<<eIslandBCEBRecHits[i][kk]<<" "<<e3x3IslandBCEB[i]<<" "<<wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<" "<<mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<endl;
                }

                  for(int kk=0 ; kk<nIslandBCEBRecHits[j] ; kk++)
                {
                  int ieta_xtal = ietaIslandBCEBRecHits[j][kk];
                  int iphi_xtal = iphiIslandBCEBRecHits[j][kk];
                  int sign = zsideIslandBCEBRecHits[j][kk]>0 ? 1 : 0;
                  wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + eIslandBCEBRecHits[j][kk]/e3x3IslandBCEB[j];
                  mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + (selePi0MinvMeanFixed_/m_inv)*(selePi0MinvMeanFixed_/m_inv)*(eIslandBCEBRecHits[j][kk]/e3x3IslandBCEB[j]);
                  cout<< "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: "<<ieta_xtal<<" "<<iphi_xtal<<" "<<sign<<" "<<eIslandBCEBRecHits[j][kk]<<" "<<e3x3IslandBCEB[j]<<" "<<wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<" "<<mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<endl;
                }

                }
            }
            }
        }
                
          
        } // End of the "j" loop over BCEB
    } // End of the "i" loop over BCEB


    } else {
      cout<< " Not enough ECAL Barrel Basic Clusters: "<<nIslandBCEB<<endl;
    }






  return kContinue;
}

// ----------------------------------------------------------------------------