CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoEgamma/EgammaTools/src/ggPFClusters.cc

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaTools/interface/ggPFClusters.h"
00002 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00003 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00004 #include "DataFormats/Math/interface/deltaR.h"
00005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00006 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
00007 #include "TMath.h"
00008 /*
00009 class by Rishi Patel rpatel@ cern.ch
00010 */
00011 
00012 
00013 ggPFClusters::ggPFClusters(
00014                            //   reco::Photon PFPhoton,
00015                            edm::Handle<EcalRecHitCollection>& EBReducedRecHits,
00016                            edm::Handle<EcalRecHitCollection>& EEReducedRecHits,
00017                            const CaloSubdetectorGeometry* geomBar,
00018                            const CaloSubdetectorGeometry* geomEnd
00019                            
00020                            ):
00021   EBReducedRecHits_(EBReducedRecHits),
00022   EEReducedRecHits_(EEReducedRecHits),
00023   geomBar_(geomBar),
00024   geomEnd_(geomEnd)
00025 {
00026   
00027 }
00028 
00029 
00030 ggPFClusters::~ggPFClusters(){
00031   
00032 }
00033 
00034 
00035 vector<reco::CaloCluster> ggPFClusters::getPFClusters(reco::SuperCluster sc){ 
00036   vector<reco::CaloCluster> PFClust;
00037   reco::CaloCluster_iterator it=sc.clustersBegin();
00038   //get PFCluster Position from Basic Clusters, get Energy from RecHits
00039   for(;it!=sc.clustersEnd();++it){
00040     std::vector< std::pair<DetId, float> >bcCells=(*it)->hitsAndFractions();
00041     DetId seedXtalId = bcCells[0].first ;
00042     int detector = seedXtalId.subdetId(); //use Seed to check if EB or EE
00043     bool isEb;
00044     if(detector==1)isEb=true;
00045     else isEb=false;
00046     
00047     float ClusSum=SumPFRecHits(bcCells, isEb); //return PFCluster Energy by Matching to RecHit and using the fractions from the bcCells
00048     CaloCluster calo(ClusSum, (*it)->position());//store in CaloCluster (parent of PFCluster)
00049     for(unsigned int i=0; i<bcCells.size(); ++i){
00050       calo.addHitAndFraction(bcCells[i].first, bcCells[i].second);//Store DetIds and Fractions
00051     }
00052     PFClust.push_back(calo);//store in the Vector
00053   }
00054   return PFClust;
00055 }
00056 
00057 float ggPFClusters::SumPFRecHits(std::vector< std::pair<DetId, float> >& bcCells, bool isEB){
00058   float ClustSum=0;
00059   for(unsigned int i=0; i<bcCells.size(); ++i){//loop over the Basic Cluster Cells
00060     EcalRecHitCollection::const_iterator eb;
00061     EcalRecHitCollection::const_iterator ee;
00062     
00063     if(isEB){
00064       for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){//loop over RecHits           
00065         if(eb->id().rawId()==bcCells[i].first.rawId()){//match
00066           float cellE=eb->energy()* bcCells[i].second;//Energy times fraction
00067           ClustSum=ClustSum+cellE;
00068         }
00069       }
00070     }  
00071     else{
00072       for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){//loop over RecHits           
00073         if(ee->id().rawId()==bcCells[i].first.rawId()){//match
00074           float cellE=ee->energy()* bcCells[i].second;//Energy times fraction
00075           ClustSum=ClustSum+cellE;
00076           break;
00077         }
00078       }
00079     }
00080   }    
00081   return ClustSum;  
00082 }
00083 
00084 float ggPFClusters::getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::SuperCluster sc){
00085   float overlap=0;
00086   reco::CaloCluster_iterator pit=sc.clustersBegin();
00087   std::vector< std::pair<DetId, float> >bcCellsPF=PFClust.hitsAndFractions();
00088   
00089   DetId seedXtalId = bcCellsPF[0].first ;
00090   int detector = seedXtalId.subdetId() ;
00091   bool isEb;
00092   std::vector< std::pair<DetId, float> >bcCellsreco;
00093   if(detector==1)isEb=true;
00094   else isEb=false;
00095   for(;pit!=sc.clustersEnd();++pit){//fill vector of basic Clusters from SuperCluster
00096     std::vector< std::pair<DetId, float> >bcCells2=(*pit)->hitsAndFractions();
00097     for(unsigned int h=0; h<bcCells2.size(); ++h)bcCellsreco.push_back(bcCells2[h]);
00098   }
00099   float clustOverlap=0;
00100   clustOverlap=PFRecHitsSCOverlap(//find overlap of a PFCluster with SuperCluster
00101                                   bcCellsPF, 
00102                                   bcCellsreco,
00103                                   isEb);
00104   overlap=clustOverlap;
00105   return overlap;
00106 
00107 
00108 }
00109 
00110 
00111 float ggPFClusters::getPFSuperclusterOverlap(reco::CaloCluster PFClust, reco::Photon phot){
00112   float overlap=0;
00113   SuperClusterRef recoSC=phot.superCluster();
00114   reco::CaloCluster_iterator pit=recoSC->clustersBegin();
00115   
00116   std::vector< std::pair<DetId, float> >bcCellsPF=PFClust.hitsAndFractions();
00117   
00118   DetId seedXtalId = bcCellsPF[0].first ;
00119   int detector = seedXtalId.subdetId() ;
00120   bool isEb;
00121   std::vector< std::pair<DetId, float> >bcCellsreco;
00122   if(detector==1)isEb=true;
00123   else isEb=false;
00124   for(;pit!=recoSC->clustersEnd();++pit){//fill vector of basic Clusters from SuperCluster
00125     std::vector< std::pair<DetId, float> >bcCells2=(*pit)->hitsAndFractions();
00126     for(unsigned int h=0; h<bcCells2.size(); ++h)bcCellsreco.push_back(bcCells2[h]);
00127     }
00128   float clustOverlap=0;
00129   clustOverlap=PFRecHitsSCOverlap(//find overlap of a PFCluster with SuperCluster
00130                                   bcCellsPF, 
00131                                   bcCellsreco,
00132                                   isEb);
00133   overlap=clustOverlap;
00134   return overlap;
00135 }
00136 
00137 float ggPFClusters::PFRecHitsSCOverlap(
00138                                        std::vector< std::pair<DetId, float> >& bcCells1, 
00139                                        std::vector< std::pair<DetId, float> >& bcCells2,
00140                                        bool isEB){
00141   float OverlapSum=0;
00142   multimap <DetId, float> pfDID;
00143   multimap <DetId, float> recoDID;
00144   multimap<DetId, float>::iterator pfit;
00145   multimap<DetId, float>::iterator pit;
00146   vector<DetId>matched;
00147   vector<float>matchedfrac;
00148   //fill Multimap of DetId, fraction for PFCluster
00149   for(unsigned int i=0; i<bcCells1.size(); ++i){ 
00150     pfDID.insert(make_pair(bcCells1[i].first, bcCells1[i].second));    
00151   }
00152    //fill Multimap of DetId, fraction for SuperCluster
00153   for(unsigned int i=0; i<bcCells2.size(); ++i){ 
00154     recoDID.insert(make_pair(bcCells2[i].first, bcCells2[i].second));    
00155   }
00156   pit=recoDID.begin();
00157   pfit=pfDID.begin();
00158   
00159   for(; pfit!=pfDID.end(); ++pfit){
00160     pit=recoDID.find(pfit->first);//match DetId from PFCluster to Supercluster
00161     if(pit!=recoDID.end()){
00162       // cout<<"Found Match "<<endl; 
00163       matched.push_back(pfit->first);//store detId
00164       matchedfrac.push_back(pfit->second);//store fraction
00165     }
00166   }
00167   
00168   for(unsigned int m=0; m<matched.size(); ++m){ //loop over matched cells
00169     DetId det=matched[m];
00170     EcalRecHitCollection::const_iterator eb;
00171     EcalRecHitCollection::const_iterator ee;
00172     if(isEB){
00173       for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){              
00174         if(eb->id().rawId()==det.rawId()){
00175           float cellE=eb->energy()* matchedfrac[m]; //compute overlap  
00176           OverlapSum=OverlapSum+cellE;
00177           break;
00178         }
00179       }
00180     }
00181     else{
00182       for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){              
00183         if(ee->id().rawId()==det.rawId()){
00184           float cellE=ee->energy() *  matchedfrac[m];//compute overlap  
00185           OverlapSum=OverlapSum+cellE;
00186           break;
00187         }
00188       }
00189     }
00190     
00191   }
00192   return OverlapSum;
00193 }
00194 
00195 void ggPFClusters::localCoordsEB( reco::CaloCluster clus, float &etacry, float &phicry, int &ieta, int &iphi, float &thetatilt, float &phitilt){
00196   const math::XYZPoint position_ = clus.position(); 
00197   double Theta = -position_.theta()+0.5*TMath::Pi();
00198   double Eta = position_.eta();
00199   double Phi = TVector2::Phi_mpi_pi(position_.phi());
00200   
00201   //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
00202   // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
00203   const float X0 = 0.89; const float T0 = 7.4;
00204   double depth = X0 * (T0 + log(clus.energy()));
00205   
00206   //find max energy crystal
00207   std::vector< std::pair<DetId, float> > crystals_vector = clus.hitsAndFractions();
00208   float drmin = 999.;
00209   EBDetId crystalseed(crystals_vector[0].first);
00210   //printf("starting loop over crystals, etot = %5f:\n",clus.energy());
00211   for (unsigned int icry=0; icry!=crystals_vector.size(); ++icry) {    
00212     
00213     EBDetId crystal(crystals_vector[icry].first);
00214     
00215     const CaloCellGeometry* cell=geomBar_->getGeometry(crystal);
00216     GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00217     double EtaCentr = center_pos.eta();
00218     double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00219     
00220     float dr = reco::deltaR(Eta,Phi,EtaCentr,PhiCentr);
00221     if (dr<drmin) {
00222       drmin = dr;
00223       crystalseed = crystal;
00224     }
00225     
00226   }
00227   
00228   ieta = crystalseed.ieta();
00229   iphi = crystalseed.iphi();
00230   
00231   // Get center cell position from shower depth
00232   const CaloCellGeometry* cell=geomBar_->getGeometry(crystalseed);
00233   const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid*>(cell);
00234   
00235    thetatilt = cpyr->getThetaAxis();
00236    phitilt = cpyr->getPhiAxis();
00237   
00238   GlobalPoint center_pos = cpyr->getPosition(depth);
00239   
00240   double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00241   double PhiWidth = (TMath::Pi()/180.);
00242   phicry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
00243   //Some flips to take into account ECAL barrel symmetries:
00244   if (ieta<0) phicry *= -1.;  
00245   
00246   double ThetaCentr = -center_pos.theta()+0.5*TMath::Pi();
00247   double ThetaWidth = (TMath::Pi()/180.)*TMath::Cos(ThetaCentr);
00248   etacry = (Theta-ThetaCentr)/ThetaWidth;
00249   //cout<<"eta, phi raw w/o widths "<<(Theta-ThetaCentr)<<", "<<(TVector2::Phi_mpi_pi(Phi-PhiCentr))<<endl;
00250   //flip to take into account ECAL barrel symmetries:
00251   if (ieta<0) etacry *= -1.;  
00252   return;
00253   
00254 }
00255 
00256 void ggPFClusters::localCoordsEE(reco::CaloCluster clus, float &xcry, float &ycry, int &ix, int &iy, float &thetatilt, float &phitilt){
00257   const math::XYZPoint position_ = clus.position(); 
00258   //double Theta = -position_.theta()+0.5*TMath::Pi();
00259   double Eta = position_.eta();
00260   double Phi = TVector2::Phi_mpi_pi(position_.phi());
00261   double X = position_.x();
00262   double Y = position_.y();
00263   
00264   //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
00265   // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
00266   const float X0 = 0.89; float T0 = 1.2;
00267   //different T0 value if outside of preshower coverage
00268   if (TMath::Abs(clus.eta())<1.653) T0 = 3.1;
00269   
00270   double depth = X0 * (T0 + log(clus.energy()));
00271   
00272   //find max energy crystal
00273   std::vector< std::pair<DetId, float> > crystals_vector = clus.hitsAndFractions();
00274   float drmin = 999.;
00275   EEDetId crystalseed(crystals_vector[0].first);
00276   //printf("starting loop over crystals, etot = %5f:\n",bclus.energy());
00277   for (unsigned int icry=0; icry!=crystals_vector.size(); ++icry) {    
00278     
00279     EEDetId crystal(crystals_vector[icry].first);
00280         
00281     const CaloCellGeometry* cell=geomEnd_->getGeometry(crystal);
00282     GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00283     double EtaCentr = center_pos.eta();
00284     double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00285 
00286     float dr = reco::deltaR(Eta,Phi,EtaCentr,PhiCentr);
00287     if (dr<drmin) {
00288       drmin = dr;
00289       crystalseed = crystal;
00290     }
00291     
00292   }
00293   
00294   ix = crystalseed.ix();
00295   iy = crystalseed.iy();
00296   
00297   // Get center cell position from shower depth
00298   const CaloCellGeometry* cell=geomEnd_->getGeometry(crystalseed);
00299   const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid*>(cell);
00300 
00301   thetatilt = cpyr->getThetaAxis();
00302   phitilt = cpyr->getPhiAxis();
00303 
00304   GlobalPoint center_pos = cpyr->getPosition(depth);
00305   
00306   double XCentr = center_pos.x();
00307   double XWidth = 2.59;
00308   xcry = (X-XCentr)/XWidth;
00309  
00310   
00311   double YCentr = center_pos.y();
00312   double YWidth = 2.59;
00313   ycry = (Y-YCentr)/YWidth;  
00314   //cout<<"x, y raw w/o widths "<<(X-XCentr)<<", "<<(Y-YCentr)<<endl;
00315 }
00316 
00317 float ggPFClusters::get5x5Element(int i, int j,
00318                                   std::vector< std::pair<DetId, float> >& 
00319                                   bcCells, 
00320                                   bool isEB){
00321   
00322   Fill5x5Map(bcCells,isEB);
00323   if(abs(i)>2 ||abs(j)>2)return 0; //outside 5x5
00324   int ind1=i+2;
00325   int ind2=j+2;
00326   float E=e5x5_[ind1][ind2]; //return element from 5x5 cells 
00327   return E;
00328 }
00329 
00330 void ggPFClusters::Fill5x5Map(std::vector< std::pair<DetId, float> >& bcCells, 
00331                               bool isEB){
00332   
00333   for(int i=0; i<5; ++i)
00334     for(int j=0; j<5; ++j)e5x5_[i][j]=0;
00335   //cout<<"here 2 "<<endl;
00336   EcalRecHitCollection::const_iterator eb;
00337   EcalRecHitCollection::const_iterator ee;
00338   
00339   DetId idseed=FindSeed(bcCells, isEB);//return seed crystal
00340   
00341   for(unsigned int i=0; i<bcCells.size(); ++i){
00342     DetId id=bcCells[i].first;
00343     if(isEB){
00344       EBDetId EBidSeed=EBDetId(idseed.rawId());
00345       int deta=EBDetId::distanceEta(id,idseed);
00346       int dphi=EBDetId::distancePhi(id,idseed);
00347       if(abs(dphi)<=2 && abs(deta)<=2){
00348         for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00349           if(eb->id().rawId()==bcCells[i].first.rawId()){
00350             EBDetId EBid=EBDetId(id.rawId());
00351             int ind1=EBidSeed.ieta()-EBid.ieta();
00352             int ind2=EBid.iphi()-EBidSeed.iphi();
00353             if(EBidSeed.ieta() * EBid.ieta() > 0){
00354               ind1=EBid.ieta()-EBidSeed.ieta();
00355             }
00356             else{ //near EB+ EB-
00357               //bugged
00358               //ind1=(1-(EBidSeed.ieta()-EBid.ieta())); 
00359               ind1 =  (EBid.ieta()-EBidSeed.ieta())*(abs(EBid.ieta()-EBidSeed.ieta())-1)/abs(EBid.ieta()-EBidSeed.ieta());
00360             }
00361             int iEta=ind1+2;
00362             int iPhi=ind2+2;
00363             e5x5_[iEta][iPhi]=eb->energy()* bcCells[i].second;
00364           }       
00365         }
00366       }
00367     }
00368     else{
00369       EEDetId EEidSeed=EEDetId(idseed.rawId());
00370       int dx=EEDetId::distanceX(id,idseed);
00371       int dy=EEDetId::distanceY(id,idseed);
00372       if(abs(dx)<=2 && abs(dy)<=2){
00373         for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00374           if(ee->id().rawId()==bcCells[i].first.rawId()){
00375             EEDetId EEid=EEDetId(id.rawId());    
00376             int ind1=EEid.ix()-EEidSeed.ix();
00377             int ind2=EEid.iy()-EEidSeed.iy();
00378             int ix=ind1+2;
00379             int iy=ind2+2;
00380             e5x5_[ix][iy]=ee->energy()* bcCells[i].second;
00381           }
00382         }
00383       }   
00384     }
00385   }
00386 }
00387 
00388 DetId ggPFClusters::FindSeed(std::vector< std::pair<DetId, float> >& bcCells, bool isEB){
00389   //first find seed:
00390   EcalRecHitCollection::const_iterator eb;
00391   EcalRecHitCollection::const_iterator ee;
00392   DetId idseed=bcCells[0].first;
00393   float PFSeedE=0;
00394   //find seed by largest energy matching
00395   for(unsigned int i=0; i<bcCells.size(); ++i){
00396     if(isEB){      
00397       for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){              
00398         if(eb->id().rawId()==bcCells[i].first.rawId()){
00399           float cellE=eb->energy()* bcCells[i].second;
00400           if(PFSeedE<cellE){                  
00401             PFSeedE=cellE;
00402             idseed=bcCells[i].first;
00403           }
00404           break;
00405         }
00406       }                 
00407     }
00408     else{
00409       for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00410         
00411         if(ee->id().rawId()==bcCells[i].first.rawId()){
00412           float cellE=ee->energy()* bcCells[i].second;
00413           if(PFSeedE<cellE){
00414             PFSeedE=cellE;
00415             idseed=bcCells[i].first;
00416           }
00417           break;
00418         }
00419       }
00420     }
00421   }
00422   return idseed;
00423 }
00424 
00425 std::pair<float, float> ggPFClusters::ClusterWidth(
00426                                                    vector<reco::CaloCluster>& PFClust){
00427   pair<float, float> widths(0,0);
00428   multimap<float, unsigned int>ClusterMap;
00429   //fill in Multimap to order by Energy (ascending order
00430   for(unsigned int i=0; i<PFClust.size();++i)ClusterMap.insert(make_pair(PFClust[i].energy(), i));
00431   //reverse iterator to get cluster with largest first 
00432   std::multimap<float, unsigned int>::reverse_iterator it;
00433   it=ClusterMap.rbegin();
00434   unsigned int max_c=(*it).second;
00435   std::vector< std::pair<DetId, float> >bcCells=PFClust[max_c].hitsAndFractions();
00436   EcalRecHitCollection::const_iterator eb;
00437   EcalRecHitCollection::const_iterator ee;
00438   float numeratorEtaWidth=0;
00439   float numeratorPhiWidth=0;
00440   float denominator=0;
00441   
00442   DetId seedXtalId = bcCells[0].first ;
00443   int detector = seedXtalId.subdetId() ;
00444   bool isEb;
00445   if(detector==1)isEb=true;
00446   else isEb=false;
00447   
00448   
00449   if(isEb){
00450     for(unsigned int i=0; i<bcCells.size(); ++i){
00451       for(eb=EBReducedRecHits_->begin();eb!= EBReducedRecHits_->end();++eb){
00452         if(eb->id().rawId()==bcCells[i].first.rawId()){
00453           
00454           double energyHit = eb->energy()*bcCells[i].second;
00455           DetId id=bcCells[i].first;
00456           float eta=geomBar_->getGeometry(id)->getPosition().eta();
00457           float phi=geomBar_->getGeometry(id)->getPosition().phi();
00458           float dEta = eta - PFClust[max_c].eta();
00459           float dPhi = phi - PFClust[max_c].phi();
00460           if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
00461           if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
00462           numeratorEtaWidth=dEta*dEta*energyHit+numeratorEtaWidth;
00463           numeratorPhiWidth=dPhi*dPhi*energyHit+numeratorPhiWidth;
00464           denominator=energyHit+denominator;
00465           break;
00466         }
00467       }
00468       
00469     }
00470   }
00471   else{
00472     for(unsigned int i=0; i<bcCells.size(); ++i){
00473       for(ee=EEReducedRecHits_->begin();ee!= EEReducedRecHits_->end();++ee){
00474         if(ee->id().rawId()==bcCells[i].first.rawId()){   
00475           double energyHit = ee->energy()*bcCells[i].second;
00476           DetId id=bcCells[i].first;
00477           float eta=geomEnd_->getGeometry(id)->getPosition().eta();
00478           float phi=geomEnd_->getGeometry(id)->getPosition().phi();
00479           float dEta = eta - PFClust[max_c].eta();
00480           float dPhi = phi - PFClust[max_c].phi();
00481           if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
00482           if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
00483           numeratorEtaWidth=dEta*dEta*energyHit+numeratorEtaWidth;
00484           numeratorPhiWidth=dPhi*dPhi*energyHit+numeratorPhiWidth;
00485           denominator=energyHit+denominator;
00486           break;          
00487         }
00488       }   
00489     }
00490   }
00491   float etaWidth=sqrt(numeratorEtaWidth/denominator);
00492   float phiWidth=sqrt(numeratorPhiWidth/denominator);
00493   widths=make_pair(etaWidth, phiWidth);
00494   return widths;
00495 }
00496 
00497 double ggPFClusters::LocalEnergyCorrection(const GBRForest *ReaderLCEB, const GBRForest *ReaderLCEE,reco::CaloCluster PFClust,float beamspotZ ){
00498   
00499   //double Coor=1;
00500   double PFClustCorr=PFClust.energy();
00501   float ClusEta=PFClust.eta();
00502   float ClusE=PFClust.energy();
00503   float ClusPhi=PFClust.phi();
00504   
00505   std::vector<float>inputs;
00506   std::vector< std::pair<DetId, float> >bcCells=PFClust.hitsAndFractions();
00507   bool isEb;
00508   DetId seedXtalId = bcCells[0].first ;
00509   int detector = seedXtalId.subdetId();
00510   if(detector==1)isEb=true;
00511   else isEb=false;
00512   //Shower Shape Variables:
00513   Fill5x5Map(bcCells, isEb);
00514   float Eseed=get5x5Element(0, 0, bcCells, isEb);
00515   
00516   float Etop=get5x5Element(0, 1, bcCells, isEb);
00517   float Ebottom=get5x5Element(0, -1, bcCells, isEb);
00518   float Eleft=get5x5Element(-1, 0, bcCells, isEb);
00519   float Eright=get5x5Element(1, 0, bcCells, isEb);
00520   float E5x5=0; float E3x3=0;
00521   for(int i=-2; i<3; ++i)
00522     for(int j=-2; j<3; ++j){
00523       float e=get5x5Element(i, j, bcCells, isEb);
00524       if(abs(i)<2)E3x3=E3x3+e;
00525       E5x5=e+E5x5;
00526     }
00527   //Fill5x5Map(bcCells, isEb);
00528   float E3x1=e5x5_[2][2]+e5x5_[1][2]+e5x5_[3][2];
00529   float E1x3=e5x5_[2][2]+e5x5_[2][1]+e5x5_[2][3];
00530   float E1x5=e5x5_[2][2]+e5x5_[2][0]+e5x5_[2][1]+e5x5_[2][3]+e5x5_[2][4];
00531   
00532   float E2x5Top=e5x5_[0][4]+e5x5_[1][4]+e5x5_[2][4]+e5x5_[3][4]+e5x5_[4][4]
00533     +e5x5_[0][3]+e5x5_[1][3]+e5x5_[2][3]+e5x5_[3][3]+e5x5_[4][3];
00534   //add up bottom edge of 5x5 2 rows
00535   float E2x5Bottom=e5x5_[0][0]+e5x5_[1][0]+e5x5_[2][0]+e5x5_[3][0]+e5x5_[4][0]
00536     +e5x5_[0][1]+e5x5_[1][1]+e5x5_[2][1]+e5x5_[3][1]+e5x5_[4][1];
00537   //add up left edge of 5x5 2 rows
00538   float E2x5Left=e5x5_[0][1]+e5x5_[0][1]+e5x5_[0][2]+e5x5_[0][3]+e5x5_[0][4]
00539     +e5x5_[1][0]+e5x5_[1][1]+e5x5_[1][2]+e5x5_[1][3]+e5x5_[1][4];
00540   //add up right edge of 5x5 2 rows
00541   float E2x5Right=e5x5_[4][0]+e5x5_[4][1]+e5x5_[4][2]+e5x5_[4][3]+e5x5_[4][4]
00542     +e5x5_[3][0]+e5x5_[3][1]+e5x5_[3][2]+e5x5_[3][3]+e5x5_[3][4];
00543   //find max 2x5 from the center
00544   float centerstrip=e5x5_[2][2]+e5x5_[2][0]+e5x5_[2][1]+e5x5_[2][3]+e5x5_[2][4];
00545   float rightstrip=e5x5_[3][2]+e5x5_[3][0]+e5x5_[3][1]+e5x5_[3][3]+e5x5_[3][4];
00546   float leftstrip=e5x5_[1][2]+e5x5_[1][0]+e5x5_[1][1]+e5x5_[1][3]+e5x5_[1][4];
00547   float E2x5Max=0;
00548   if(rightstrip>leftstrip)E2x5Max=rightstrip+centerstrip;
00549   else E2x5Max=leftstrip+centerstrip;
00550   //get Local Coordinates
00551   if(isEb){
00552     float etacry; float phicry; int ieta; int iphi; float thetatilt; float phitilt;
00553     int iEtaCrack=3;
00554     if(abs(ieta)==1 || abs(ieta)==2 )
00555       iEtaCrack=abs(ieta);
00556     if(abs(ieta)>2 && abs(ieta)<24)
00557       iEtaCrack=3;
00558     if(abs(ieta)==24)
00559       iEtaCrack=4;
00560     if(abs(ieta)==25)
00561       iEtaCrack=5;
00562     if(abs(ieta)==26)
00563       iEtaCrack=6;
00564     if(abs(ieta)==27)
00565       iEtaCrack=7;
00566     if(abs(ieta)>27 &&  abs(ieta)<44)
00567       iEtaCrack=8;
00568     if(abs(ieta)==44)
00569       iEtaCrack=9;
00570     if(abs(ieta)==45)
00571       iEtaCrack=10;
00572     if(abs(ieta)==46)
00573       iEtaCrack=11;
00574     if(abs(ieta)==47)
00575       iEtaCrack=12;
00576     if(abs(ieta)>47 &&  abs(ieta)<64)
00577       iEtaCrack=13;
00578     if(abs(ieta)==64)
00579       iEtaCrack=14;
00580     if(abs(ieta)==65)
00581       iEtaCrack=15;
00582     if(abs(ieta)==66)
00583       iEtaCrack=16;
00584     if(abs(ieta)==67)
00585       iEtaCrack=17;
00586     if(abs(ieta)>67 &&  abs(ieta)<84)
00587       iEtaCrack=18;
00588     if(abs(ieta)==84)
00589       iEtaCrack=19;
00590     if(abs(ieta)==85)
00591       iEtaCrack=20;
00592     localCoordsEB(PFClust, etacry, phicry, ieta, iphi, thetatilt, phitilt);
00593     inputs.push_back(beamspotZ);
00594     inputs.push_back(ClusEta/fabs(ClusEta));
00595     inputs.push_back(fabs(ClusEta));
00596     inputs.push_back(fabs(ClusPhi));
00597     inputs.push_back(log(ClusE));
00598     inputs.push_back((Eseed/ClusE));
00599     inputs.push_back((Etop/ClusE));
00600     inputs.push_back((Ebottom/ClusE));
00601     inputs.push_back((Eleft/ClusE));
00602     inputs.push_back((Eright/ClusE));
00603     inputs.push_back(E3x3/ClusE);
00604     inputs.push_back(E1x3/ClusE);
00605     inputs.push_back(E3x1/ClusE);
00606     inputs.push_back(E5x5/ClusE);
00607     inputs.push_back(E1x5/ClusE);
00608     inputs.push_back(E2x5Max/ClusE);
00609     inputs.push_back(E2x5Top/ClusE);
00610     inputs.push_back(E2x5Bottom/ClusE);
00611     inputs.push_back(E2x5Left/ClusE);
00612     inputs.push_back(E2x5Right/ClusE);
00613     inputs.push_back(etacry);
00614     inputs.push_back(phicry);
00615     inputs.push_back(iphi%2);
00616     inputs.push_back(ieta%5);
00617     inputs.push_back(iphi%20);
00618     inputs.push_back(iEtaCrack);
00619     int size=inputs.size();
00620     float PFInputs[26];
00621     for(int i=0; i<size; ++i)PFInputs[i]=inputs[i];
00622     PFClustCorr= ReaderLCEB->GetResponse(PFInputs)*ClusE;
00623   }
00624   else{    
00625     float xcry; float ycry; int ix; int iy; float thetatilt; float phitilt;
00626     localCoordsEE(PFClust, xcry, ycry, ix, iy, thetatilt, phitilt);
00627     inputs.push_back(beamspotZ);
00628     inputs.push_back(ClusEta/fabs(ClusEta));
00629     inputs.push_back(fabs(ClusEta));
00630     inputs.push_back(fabs(ClusPhi));
00631     inputs.push_back(log(ClusE));
00632     inputs.push_back((Eseed/ClusE));
00633     inputs.push_back((Etop/ClusE));
00634     inputs.push_back((Ebottom/ClusE));
00635     inputs.push_back((Eleft/ClusE));
00636     inputs.push_back((Eright/ClusE));
00637     inputs.push_back(E3x3/ClusE);
00638     inputs.push_back(E1x3/ClusE);
00639     inputs.push_back(E3x1/ClusE);
00640     inputs.push_back(E5x5/ClusE);
00641     inputs.push_back(E1x5/ClusE);
00642     inputs.push_back(E2x5Max/ClusE);
00643     inputs.push_back(E2x5Top/ClusE);
00644     inputs.push_back(E2x5Bottom/ClusE);
00645     inputs.push_back(E2x5Left/ClusE);
00646     inputs.push_back(E2x5Right/ClusE);
00647     inputs.push_back(xcry);
00648     inputs.push_back(ycry);
00649     int size=inputs.size();
00650     float PFInputs[22];
00651     for(int i=0; i<size; ++i)PFInputs[i]=inputs[i];
00652     PFClustCorr= ReaderLCEE->GetResponse(PFInputs) *ClusE;
00653   }
00654   
00655   return PFClustCorr;
00656 }
00657 void ggPFClusters::BasicClusterPFCandLink(           
00658                                           reco::SuperCluster sc, 
00659                                           std::vector<reco::PFCandidatePtr>&insideBox,
00660                                           std::vector<DetId>& MatchedRH
00661                                           ){
00662   std::vector<reco::PFCandidatePtr>Linked;
00663   for(unsigned int p=0; p<insideBox.size(); ++p){
00664     math::XYZPointF position_ = insideBox[p]->positionAtECALEntrance();
00665     //math::XYZPointF positionvtx(position_.x()+insideBox[p]->vx(), 
00666     //                  position_.y()+insideBox[p]->vy(),
00667     //                  position_.z()+insideBox[p]->vz());
00668     //position_=positionvtx;
00669     //math::XYZVector position_=insideBox[p]->momentum();
00670     if(insideBox[p]->pdgId()==22){
00671       double Theta = -position_.theta()+0.5*TMath::Pi();
00672       double Eta = position_.eta();
00673       double Phi = TVector2::Phi_mpi_pi(position_.phi());
00674       double X = position_.x();
00675       double Y = position_.y();
00676       reco::CaloCluster_iterator cit=sc.clustersBegin();
00677       std::vector< std::pair<DetId, float> > crystals_vector = (*cit)->hitsAndFractions();
00678       DetId seedXtalId = crystals_vector[0].first;
00679       int detector = seedXtalId.subdetId();
00680       bool isEb=false;
00681       float X0 = 0.89; float T0 = 7.4;
00682       double depth = X0 * (T0 + log((*cit)->energy()));
00683       if(detector==1){
00684         X0 = 0.89;  T0 = 7.4;
00685         depth = X0 * (T0 + log((*cit)->energy()));
00686         isEb=true;
00687       }
00688       else{
00689         X0 = 0.89; T0 = 1.2;
00690         if(fabs(Eta)<1.653)T0=3.1;
00691         depth = X0 * (T0 + log((*cit)->energy()));
00692         
00693         isEb=false;
00694       }
00695       crystals_vector.clear();
00696       for(; cit!=sc.clustersEnd(); ++cit){
00697         bool matchBC=false;
00698         crystals_vector = (*cit)->hitsAndFractions();
00699         
00700         
00701         for (unsigned int icry=0; icry<crystals_vector.size(); ++icry){
00702           
00703           if(isEb){
00704             
00705             EBDetId crystal(crystals_vector[icry].first);
00706             const CaloCellGeometry* cell=geomBar_->getGeometry(crystal);
00707             
00708             GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00709             //GlobalPoint vtx(insideBox[p]->vx(), insideBox[p]->vy(), insideBox[p]->vz());
00710             //GlobalPoint center_pos(oldcenter_pos.x()-vtx.x(), 
00711             //             oldcenter_pos.y()-vtx.y(),
00712             //             oldcenter_pos.z()-vtx.z());
00713             //double EtaCentr = center_pos.eta();
00714              double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
00715              double PhiWidth = (TMath::Pi()/180.);
00716              double ThetaCentr = -center_pos.theta()+0.5*TMath::Pi();
00717              double ThetaWidth = (TMath::Pi()/180.)*TMath::Cos(ThetaCentr);
00718              double phicry = (TVector2::Phi_mpi_pi(Phi-PhiCentr))/PhiWidth;
00719              double etacry=(Theta-ThetaCentr)/ThetaWidth;
00720              if(fabs(etacry)<0.6 && fabs(phicry)<0.6){
00721                Linked.push_back(insideBox[p]);
00722                matchBC=true;
00723                MatchedRH.push_back(crystals_vector[icry].first);
00724                break;
00725                
00726              }
00727              
00728           }//Barrel
00729           else{
00730              
00731             EEDetId crystal(crystals_vector[icry].first);
00732             const CaloCellGeometry* cell=geomEnd_->getGeometry(crystal);
00733             
00734             GlobalPoint center_pos = (dynamic_cast<const TruncatedPyramid*>(cell))->getPosition(depth);
00735             //GlobalPoint vtx(insideBox[p]->vx(), insideBox[p]->vy(), insideBox[p]->vz());
00736             //GlobalPoint center_pos(oldcenter_pos.x()-vtx.x(), oldcenter_pos.y()-vtx.y(), oldcenter_pos.z()-vtx.z());
00737             double XCentr = center_pos.x();
00738             double XWidth = 2.59;
00739             double xcry = (X-XCentr)/XWidth;
00740             double YCentr = center_pos.y();
00741             double YWidth = 2.59;
00742             double ycry = (Y-YCentr)/YWidth;  
00743             if(fabs(xcry)<0.6 && fabs(ycry)<0.6){
00744               Linked.push_back(insideBox[p]);
00745               MatchedRH.push_back(crystals_vector[icry].first);
00746               matchBC=true;
00747               break;
00748             }
00749             
00750           }//Endcap
00751           
00752         }
00753         if(matchBC)break;
00754       }      
00755     }
00756     if(abs(insideBox[p]->pdgId())==211){
00757       float drmin = 999.;
00758       float deta=999;
00759       float dphi=999;
00760       reco::CaloCluster_iterator cit=sc.clustersBegin();
00761       DetId seedCrystal=(*cit)->hitsAndFractions()[0].first;
00762       for(; cit!=sc.clustersEnd(); ++cit){
00763         DetId crystals_vector_seed=(*cit)->hitsAndFractions()[0].first;
00764         math::XYZVector photon_directionWrtVtx((*cit)->position().x()-insideBox[p]->vx(), 
00765                                                (*cit)->position().y()-insideBox[p]->vy(),
00766                                                (*cit)->position().z()-insideBox[p]->vz()
00767                                                );
00768         float dR=deltaR(photon_directionWrtVtx.eta(), photon_directionWrtVtx.phi(), position_.eta(), position_.phi());
00769         //float dR=deltaR((*cit)->eta(), (*cit)->phi(), position_.eta(), position_.phi());
00770         if(dR<drmin){
00771           drmin=dR;
00772           deta=fabs(photon_directionWrtVtx.eta()-position_.eta());
00773           dphi=acos(cos(photon_directionWrtVtx.phi()- position_.phi()));
00774           seedCrystal=crystals_vector_seed;
00775         }
00776       }
00777       if(deta<0.05 && dphi<0.07){
00778         Linked.push_back(insideBox[p]);
00779         MatchedRH.push_back(seedCrystal);
00780       }
00781     }
00782     
00783   }
00784   insideBox.clear();
00785   insideBox=Linked;
00786 }