XANASuperClusters.cpp Source File

00001 
00002 #include <Utilities/Configuration/interface/Architecture.h>
00003 
00004 #include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
00005 #include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
00006 #include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
00007 #include <XANADOO/XANAClusters/interface/XANACellID.h>
00008 
00009 #include <ElectronPhoton/EgammaClusters/interface/EgammaBasicCluster.h>
00010 #include <ElectronPhoton/EgammaClusters/interface/EgammaBasicDynamic.h>
00011 #include <ElectronPhoton/EgammaClusters/interface/EgammaBasicIsland.h>
00012 #include <ElectronPhoton/EgammaClusters/interface/EgammaBasicHybrid.h>
00013 #include <ElectronPhoton/EgammaClusters/interface/EGBLogPosCorr.h>
00014 #include <ElectronPhoton/EgammaAnalysis/interface/EgammaHAnalyzer.h>
00015 #include <ElectronPhoton/EgammaAnalysis/interface/EgammaAlgorithm.h>
00016 #include <ElectronPhoton/EgammaClusters/interface/EGSCSeedGenerator.h>
00017 #include <ElectronPhoton/ClusterTools/interface/EgammaBremRecovery.h>
00018 #include <ElectronPhoton/ClusterTools/interface/EgammaSuperCluster.h>
00019 #include <ElectronPhoton/EgammaPreshower/interface/EgammaEndcapCluster.h>
00020 #include <ElectronPhoton/ClusterTools/interface/EgammaSCEnergyScale.h>
00021 #include <ElectronPhoton/EgammaPreshower/interface/EgammaEndcapAssociation.h>
00022 #include <ElectronPhoton/ClusterTools/interface/IsoCalculator.h>
00023 #include <ElectronPhoton/EgammaAnalysis/interface/PreshowerHits.h>
00024 //#include <ElectronPhoton/ElectronFromPixels/interface/EPBremRecovery.h>
00025 
00026 #include <Calorimetry/CaloRecHit/interface/CaloRecHit.h>
00027 
00028 #include <CARF/G3Event/interface/G3EventProxy.h>
00029 #include <CARF/Reco/interface/AutoRecCollection.h>
00030 #include <CARF/Reco/interface/RecItr.h>
00031 
00032 #include <Utilities/Notification/interface/Observer.h>
00033 #include <Utilities/Notification/interface/PackageInitializer.h>
00034 #include <Utilities/UI/interface/SimpleConfigurable.h> 
00035 #include <Utilities/Notification/interface/Singleton.h>
00036 
00037 #include <vector>
00038 
00039 struct ClusterObserver : public Observer<G3EventProxy *> {
00040 
00041   int nbClusters_ ;
00042   double energySum_ ;
00043   SimpleConfigurable<bool> writeCaloRecHits_;
00044   SimpleConfigurable<bool> produceReference_;
00045   SimpleConfigurable<std::string> clusterizer_;
00046   SimpleConfigurable<std::string> superClusterizer_;
00047   SimpleConfigurable<std::string> endcapSuperClusterizer_;
00048  
00049   IsoCalculator *ciso_; 
00050   
00051   ClusterObserver()
00052    : nbClusters_(0), energySum_(0)
00053    { 
00054     // initialize observer
00055     Observer<G3EventProxy *>::init() ;      
00056     
00057     // steering
00058     writeCaloRecHits_=SimpleConfigurable<bool>(false,"XANAClusters:ESDWithCaloRecHits");
00059     produceReference_=SimpleConfigurable<bool>(false,"XANAClusters:ProduceReference");
00060     clusterizer_=SimpleConfigurable<std::string>("EGBCluster","XANAClusters:Clusterizer");
00061     superClusterizer_=SimpleConfigurable<std::string>("EGSCluster","XANAClusters:SuperClusterizer");
00062     endcapSuperClusterizer_=SimpleConfigurable<std::string>("EGECluster","XANAClusters:EndcapSuperClusterizer");
00063 /* 
00064     // defines clustering algo  
00065     EgammaHAnalyzer<EgammaBasicCluster,EGDummy,CaloRecHit> *clustersUnit_=0; 
00066     vector<const EgammaAlgorithm<EgammaBasicCluster,EGDummy,CaloRecHit> * > algos;
00067     if (clusterizer_.value()=="EgammaBasicDynamic") {    
00068       algos.push_back(new EgammaBasicDynamic);
00069       clustersUnit_ = new EgammaHAnalyzer<EgammaBasicCluster,EGDummy,CaloRecHit>(algos,0) ;      
00070      } else if (clusterizer_.value()=="EgammaBasicIsland") { 
00071       algos.push_back(new EgammaBasicIsland); 
00072       clustersUnit_ = new EgammaHAnalyzer<EgammaBasicCluster,EGDummy,CaloRecHit>(algos,0) ;
00073      } else if (clusterizer_.value()=="EgammaBasicHybrid") { 
00074       algos.push_back(new EgammaBasicHybrid);
00075       clustersUnit_ = new EgammaHAnalyzer<EgammaBasicCluster,EGDummy,CaloRecHit>(algos,0) ;
00076      } else if (clusterizer_.value()=="EgammaDefault") { 
00077       // Egamma code uses two algorithms for clustering
00078       EgammaBasicIsland *island = new EgammaBasicIsland;
00079       EgammaBasicHybrid *hybrid = new EgammaBasicHybrid;
00080       algos.push_back(island);
00081       algos.push_back(hybrid);
00082       clustersUnit_ = new EgammaHAnalyzer<EgammaBasicCluster,EGDummy,CaloRecHit>(algos,0);
00083       // island corrections
00084       clustersUnit_->SetCalib(new EgammaCalibrator<EgammaBasicCluster>(new EGBLogPosCorr(island)));
00085       clustersUnit_->AddCalib(new EgammaCalibrator<EgammaBasicCluster>(new EGBLogPosCorr(hybrid)));
00086      } else { 
00087        throw Genexception("[ClusterObserver::ClusterObserver()], wrong cluster algo name");      
00088     }
00089     clustersUnit_->changeName("EgammaBasicCluster") ;    
00090  
00091     Reconstructor::addDefault(clustersUnit_) ;
00092 
00093     // from EgammaSuperClusters/src/EgammaReco.cc
00094     EGSCSeedGenerator *seedgen = new EGSCSeedGenerator();
00095     EgammaHAnalyzer<EgammaCluster, EGDummy, EgammaBasicCluster> *bd =
00096      new EgammaHAnalyzer<EgammaCluster, EGDummy, EgammaBasicCluster>(seedgen,0);
00097    
00098     // SuperCluster and brem recovery
00099     EgammaHAnalyzer<EgammaSuperCluster,EgammaCluster,EgammaBasicCluster> *be = 0;
00100     if (superClusterizer_.value()=="EPBremRecovery") {
00101       EPBremRecovery *brec = new EPBremRecovery(algos);
00102       be = new EgammaHAnalyzer<EgammaSuperCluster,EgammaCluster,EgammaBasicCluster>(brec,0);
00103     } else if (superClusterizer_.value()=="EgammaDefault") {
00104       EgammaBremRecovery *brec = new EgammaBremRecovery();
00105       be = new EgammaHAnalyzer<EgammaSuperCluster,EgammaCluster,EgammaBasicCluster>(brec,0);
00106     } else
00107      throw Genexception("[ClusterObserver::ClusterObserver()], wrong brem recoverer algo name");
00108   
00109     // SuperCluster calibration using Chris's escale functions
00110     be->SetCalib(new EgammaCalibrator<EgammaSuperCluster>(new EgammaSCEnergyScale()));
00111 
00112     // Endcap Preshower association/reconstruction
00113     EgammaEndcapAssociation *egea = new EgammaEndcapAssociation();
00114     EgammaHAnalyzer<EgammaEndcapCluster,EgammaSuperCluster,CaloRecHit> *eea =
00115      new EgammaHAnalyzer<EgammaEndcapCluster,EgammaSuperCluster,CaloRecHit>(egea,0);
00116 
00117     std::cout << "adding superclusters algos " << std::endl;
00118     Reconstructor::addDefault(bd);
00119     Reconstructor::addDefault(be);
00120     Reconstructor::addDefault(eea);
00121 */   
00122 
00123     // finally define isolation algorithm
00124     ciso_ = Singleton<IsoCalculator>::instance();
00125    
00126    }
00127    
00128   ~ClusterObserver()
00129    {
00130     std::cout<<"[OVAL] number of clusters : "<<nbClusters_<<std::endl ;
00131     std::cout<<"[OVAL] mean energy : "<<energySum_/nbClusters_<<std::endl ;
00132    }
00133 
00134    
00135   void upDate( G3EventProxy * ev )
00136    { 
00137     if (!ev) return ;       
00138 
00139     std::cout << "updating .." << std::endl;
00140     // from EgammaSuperClusters/src/EgammaReco.cc
00141     // Egamma code is really messy: all superclusters becomes EndcapClusters at the end !!
00142     RecCollection<EgammaEndcapCluster> endcapSuperClusters( RecQuery(endcapSuperClusterizer_.value()));
00143     RecCollection<EgammaEndcapCluster>::const_iterator endcapSuperCluster ;  
00144     // to produce reference file
00145     if (produceReference_) {
00146      // super clusters are at the end EndCapSuperClusters!!
00147      for (endcapSuperCluster=endcapSuperClusters.begin(); endcapSuperCluster!=endcapSuperClusters.end(); endcapSuperCluster++)      
00148       {
00149        // remove double reconstruction in case of EgammaDefault
00150        if (endcapSuperClusterizer_.value() == "EGECluster") {
00151          // hybrid is algo>=100   
00152          if ( (endcapSuperCluster->seed()->WhichAlgo()>99 && 
00153              fabs(endcapSuperCluster->eta())>1.4791) || 
00154             (endcapSuperCluster->seed()->WhichAlgo()<99 && 
00155              fabs(endcapSuperCluster->eta())<1.4791) ) continue;
00156        }
00157        nbClusters_++;
00158        energySum_ += endcapSuperCluster->energy() ;
00159        std::cout<<"[OVAL] new super cluster with energy "<<endcapSuperCluster->energy()<<endl ;
00160       }
00161       return ;
00162     } 
00163 
00164     // construct XANASuperClusters from SuperClusters
00165     vector<XANASuperCluster> xsclusters;
00166     XANASuperCluster *xscluster;
00167 
00168     // loop on endcaps super clusters
00169     for (endcapSuperCluster=endcapSuperClusters.begin(); endcapSuperCluster!=endcapSuperClusters.end(); endcapSuperCluster++) 
00170      {
00171       // remove double reconstruction in case of EgammaDefault
00172       if (endcapSuperClusterizer_.value() == "EGECluster") {
00173         // Egamma beautiful code: hybrid is algo>=100 !!  
00174         if ( (endcapSuperCluster->seed()->WhichAlgo()>99 && 
00175             fabs(endcapSuperCluster->eta())>1.4791) || 
00176            (endcapSuperCluster->seed()->WhichAlgo()<99 && 
00177             fabs(endcapSuperCluster->eta())<1.4791) ) continue;
00178       }
00179       nbClusters_++;
00180       energySum_ += endcapSuperCluster->energy() ;
00181       std::cout << "[update] constructing new supercluster from algo " << endcapSuperCluster->seed()->WhichAlgo()
00182       << " at eta " << endcapSuperCluster->eta()<< std::endl;
00183       xscluster = new XANASuperCluster;
00184       xscluster->setEnergy(endcapSuperCluster->energy()); 
00185       xscluster->setPosition(endcapSuperCluster->position()); 
00186       xscluster->setSum1(endcapSuperCluster->seed()->getS1()); 
00187       xscluster->setSum4(endcapSuperCluster->seed()->getS4()); 
00188       xscluster->setSum9(endcapSuperCluster->seed()->getS9()); 
00189       xscluster->setSum25(endcapSuperCluster->seed()->getS25()); 
00190       xscluster->setHadronicOverEm(endcapSuperCluster->seed()->getHoe()); 
00191       // iso not available from endcap super clusters
00192       // compute from IsoCalculator
00193       float iso = ciso_->getIso(endcapSuperCluster->eta(),endcapSuperCluster->phi());
00194       xscluster->setCaloIsolation(iso); 
00195       xscluster->setDisc1(endcapSuperCluster->seed()->getDisc1()); 
00196       xscluster->setDisc2(endcapSuperCluster->seed()->getDisc2()); 
00197       xscluster->setDisc3(endcapSuperCluster->seed()->getDisc3()); 
00198       HepSymMatrix cov(2);
00199       cov(1,1) = endcapSuperCluster->seed()->getCovEE();
00200       cov(1,2) = endcapSuperCluster->seed()->getCovEP();
00201       cov(2,2) = endcapSuperCluster->seed()->getCovPP();
00202       xscluster->setPositionCovarianceMatrix(cov);
00203       xscluster->setEnergyScaleFactor(endcapSuperCluster->energy()/endcapSuperCluster->energyUncorrected());
00204       xscluster->setAlgoName(endcapSuperClusterizer_.value());
00205       float theta = 2.*atan(exp(-endcapSuperCluster->seed()->Eta()));
00206       float energy = endcapSuperCluster->seed()->Et()/sin(theta);
00207       float x = endcapSuperCluster->seed()->Radius()*sin(theta)*cos(endcapSuperCluster->seed()->Phi());
00208       float y = endcapSuperCluster->seed()->Radius()*sin(theta)*sin(endcapSuperCluster->seed()->Phi());
00209       float z = endcapSuperCluster->seed()->Radius()*cos(theta);
00210       // Chi2 and NumberOfDigis not available from egamma seed cluster
00211       // retreiving from the basic cluster
00212       float chi2 = 0.;
00213       int numberOfDigis = 0;
00214       EgammaVSuperCluster::const_iterator icf = endcapSuperCluster->basicClusters().first;
00215       std::cout << "[update] adding seed " << std::endl;
00216       // add seed cluster
00217       XANAEmCluster *seed = new XANAEmCluster((*icf)->energy(),
00218                           (*icf)->position(),
00219                           (*icf)->chi2(),
00220                           clusterizer_.value()); 
00221       if (writeCaloRecHits_) {
00222             std::vector<const CaloRecHit *> rhits=(*icf)->rHits();
00223             for (vector<const CaloRecHit *>::const_iterator
00224              it=rhits.begin(); it != rhits.end(); it++) {
00225               seed->addHit(new XANACaloRecHit((*it)->energy(), (*it)->time(), 
00226               (*it)->chi2(), XANACellID((*it)->getMyCell().position())));
00227             } 
00228       } else
00229         seed->setNumberOfRecHits(numberOfDigis);      
00230       xscluster->setSeedCluster(seed);
00231       std::cout << "[update] adding brems " << std::endl;
00232       // now add brems clusters
00233       XANAEmCluster *brem = 0;
00234       for (EgammaVSuperCluster::const_iterator ic = (endcapSuperCluster->basicClusters().first)+1; ic != endcapSuperCluster->basicClusters().second; ic++) {
00235         // beware !! this list includes the seed cluster
00236         // for naming, stick to the RecAlgo names
00237         std::string name = "";
00238         if ((*ic)->whichAlgo() == 1) name="EGBCluster";
00239         else if ((*ic)->whichAlgo() == 2) name="EGBDynamicCluster";
00240         else if ((*ic)->whichAlgo() == 100) name="EGBCluster";  
00241         if ((*ic)->energy() != energy) {
00242           brem = new XANAEmCluster((*ic)->energy(),(*ic)->position(),(*ic)->chi2(),name);           
00243           if (writeCaloRecHits_) {
00244             std::vector<const CaloRecHit *> rhits=(*ic)->rHits();
00245             for (vector<const CaloRecHit *>::const_iterator
00246              it=rhits.begin(); it != rhits.end(); it++) {
00247               brem->addHit(new XANACaloRecHit((*it)->energy(), (*it)->time(),
00248               (*it)->chi2(), XANACellID((*it)->getMyCell().position())));
00249             } 
00250           } else
00251             brem->setNumberOfRecHits((*ic)->numberOfDigis());
00252           std::cout << "[upDate] brem energy " << brem->getEnergy() << std::endl;
00253           xscluster->addBremCluster(brem);
00254         }
00255       } 
00256       std::cout << "[update] adding presh info " << std::endl;
00257       // add preshower info when relevant (endcap)
00258       XANAPreshowerInfo *preshowerInfo = 0;
00259 //ElectronPhoton/EgammaPreshower/interface/EgammaEndcapCluster.h:
00260 //Information about presh clusters currently not stored
00261 //      if (endcapSuperCluster->energy_P1()+endcapSuperCluster->energy_P2() != 0. ) {
00262 //        cout << "[update] constructing a  XANAPreshowerInfo" << endl;
00263 //      preshowerInfo = new XANAPreshowerInfo(endcapSuperCluster->energy_P1(), 
00264 //       endcapSuperCluster->energy_P2(), endcapSuperCluster->position_P1(),
00265 //       endcapSuperCluster->position_P2(), endcapSuperCluster->calibP1(), endcapSuperCluster->calibP2(), 
00266 //       endcapSuperCluster->calibSC());
00267 //      } 
00268       xscluster->setPreshowerInfo(preshowerInfo);
00269       std::cout << "[update] all done for this supercluster " << std::endl;
00270       xsclusters.push_back(*xscluster); 
00271     }     
00272 
00273     // read XANASuperClusters 
00274     vector<XANASuperCluster>::iterator its;
00275     for (its=xsclusters.begin(); its!=xsclusters.end(); its++)  
00276      {
00277       std::cout<<"[OVAL] new super cluster with energy "<<its->getEnergy()<<std::endl ;
00278      }
00279 
00280    }
00281       
00282 } ;
00283  
00284 static PKBuilder<ClusterObserver> clusterObserver("ClusterObserver") ;
00285 
00286