Reference Manual

Go to the documentation of this file.
00001 #include <iostream>
00002 //
00003 
00004 #include "DQMOffline/EGamma/plugins/PiZeroAnalyzer.h"
00005 
00006 
00007 //#define TWOPI 6.283185308
00008 //
00009 
00023 using namespace std;
00024 
00025 
00026 PiZeroAnalyzer::PiZeroAnalyzer( const edm::ParameterSet& pset )
00027 {
00028 
00029     fName_              = pset.getUntrackedParameter<std::string>("Name");
00030     verbosity_          = pset.getUntrackedParameter<int>("Verbosity");
00031 
00032     prescaleFactor_     = pset.getUntrackedParameter<int>("prescaleFactor",1);
00033 
00034 
00035 
00036     barrelEcalHits_     = pset.getParameter<edm::InputTag>("barrelEcalHits");
00037     endcapEcalHits_     = pset.getParameter<edm::InputTag>("endcapEcalHits");
00038 
00039 
00040     standAlone_         = pset.getParameter<bool>("standAlone");
00041 
00042 
00043 
00044 
00045 
00046     // parameters for Pizero finding
00047     seleXtalMinEnergy_    = pset.getParameter<double> ("seleXtalMinEnergy");
00048     clusSeedThr_          = pset.getParameter<double> ("clusSeedThr");
00049     clusEtaSize_          = pset.getParameter<int> ("clusEtaSize");
00050     clusPhiSize_          = pset.getParameter<int> ("clusPhiSize");
00051     ParameterLogWeighted_ = pset.getParameter<bool> ("ParameterLogWeighted");
00052     ParameterX0_          = pset.getParameter<double> ("ParameterX0");
00053     ParameterT0_barl_     = pset.getParameter<double> ("ParameterT0_barl");
00054     ParameterW0_          = pset.getParameter<double> ("ParameterW0");
00055 
00056     selePtGammaOne_       = pset.getParameter<double> ("selePtGammaOne");
00057     selePtGammaTwo_       = pset.getParameter<double> ("selePtGammaTwo");
00058     seleS4S9GammaOne_     = pset.getParameter<double> ("seleS4S9GammaOne");
00059     seleS4S9GammaTwo_     = pset.getParameter<double> ("seleS4S9GammaTwo");
00060     selePtPi0_            = pset.getParameter<double> ("selePtPi0");
00061     selePi0Iso_           = pset.getParameter<double> ("selePi0Iso");
00062     selePi0BeltDR_        = pset.getParameter<double> ("selePi0BeltDR");
00063     selePi0BeltDeta_      = pset.getParameter<double> ("selePi0BeltDeta");
00064     seleMinvMaxPi0_       = pset.getParameter<double> ("seleMinvMaxPi0");
00065     seleMinvMinPi0_       = pset.getParameter<double> ("seleMinvMinPi0");
00066 
00067     parameters_ = pset;
00068 
00069 
00070 }
00071 
00072 
00073 
00074 PiZeroAnalyzer::~PiZeroAnalyzer() {
00075 
00076 
00077 
00078 
00079 }
00080 
00081 
00082 void PiZeroAnalyzer::beginJob()
00083 {
00084 
00085 
00086   nEvt_=0;
00087   nEntry_=0;
00088 
00089   dbe_ = 0;
00090   dbe_ = edm::Service<DQMStore>().operator->();
00091 
00092 
00093 
00094  if (dbe_) {
00095     if (verbosity_ > 0 ) {
00096       dbe_->setVerbose(1);
00097     } else {
00098       dbe_->setVerbose(0);
00099     }
00100   }
00101   if (dbe_) {
00102     if (verbosity_ > 0 ) dbe_->showDirStructure();
00103   }
00104 
00105 
00106 
00107 
00108   //booking all histograms
00109 
00110   if (dbe_) {
00111 
00112     currentFolder_.str("");
00113     currentFolder_ << "Egamma/PiZeroAnalyzer/";
00114     dbe_->setCurrentFolder(currentFolder_.str());
00115 
00116 
00117 
00118 
00119     hMinvPi0EB_ = dbe_->book1D("Pi0InvmassEB","Pi0 Invariant Mass in EB",100,0.,0.5);
00120     hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ",1);
00121 
00122     hPt1Pi0EB_ = dbe_->book1D("Pt1Pi0EB","Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
00123     hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ",1);
00124 
00125     hPt2Pi0EB_ = dbe_->book1D("Pt2Pi0EB","Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
00126     hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ",1);
00127 
00128     hPtPi0EB_ = dbe_->book1D("PtPi0EB","Pi0 Pt in EB",100,0.,20.);
00129     hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ",1);
00130 
00131     hIsoPi0EB_ = dbe_->book1D("IsoPi0EB","Pi0 Iso in EB",50,0.,1.);
00132     hIsoPi0EB_->setAxisTitle("Pi0 Iso",1);
00133 
00134 
00135   }
00136 
00137 }
00138 
00139 
00140 
00141 
00142 
00143 
00144 void PiZeroAnalyzer::analyze( const edm::Event& e, const edm::EventSetup& esup )
00145 {
00146 
00147   using namespace edm;
00148 
00149   if (nEvt_% prescaleFactor_ ) return;
00150   nEvt_++;
00151   LogInfo("PiZeroAnalyzer") << "PiZeroAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
00152 
00153 
00154   // Get EcalRecHits
00155   bool validEcalRecHits=true;
00156   Handle<EcalRecHitCollection> barrelHitHandle;
00157   EcalRecHitCollection barrelRecHits;
00158   e.getByLabel(barrelEcalHits_, barrelHitHandle);
00159   if (!barrelHitHandle.isValid()) {
00160     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<barrelEcalHits_.label();
00161     validEcalRecHits=false;
00162   }
00163 
00164   Handle<EcalRecHitCollection> endcapHitHandle;
00165   e.getByLabel(endcapEcalHits_, endcapHitHandle);
00166   EcalRecHitCollection endcapRecHits;
00167   if (!endcapHitHandle.isValid()) {
00168     edm::LogError("PhotonProducer") << "Error! Can't get the product "<<endcapEcalHits_.label();
00169     validEcalRecHits=false;
00170   }
00171 
00172   if (validEcalRecHits) makePizero(esup,  barrelHitHandle, endcapHitHandle);
00173 
00174 
00175 
00176 }
00177 
00178 void PiZeroAnalyzer::makePizero ( const edm::EventSetup& es, const edm::Handle<EcalRecHitCollection> rhEB,  const edm::Handle<EcalRecHitCollection> rhEE ) {
00179 
00180   const EcalRecHitCollection *hitCollection_p = rhEB.product();
00181 
00182   edm::ESHandle<CaloGeometry> geoHandle;
00183   es.get<CaloGeometryRecord>().get(geoHandle);
00184 
00185   edm::ESHandle<CaloTopology> theCaloTopology;
00186   es.get<CaloTopologyRecord>().get(theCaloTopology);
00187 
00188 
00189   const CaloSubdetectorGeometry *geometry_p;
00190   const CaloSubdetectorTopology *topology_p;
00191   const CaloSubdetectorGeometry *geometryES_p;
00192   geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00193   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00194 
00195   // Parameters for the position calculation:
00196   edm::ParameterSet posCalcParameters =
00197     parameters_.getParameter<edm::ParameterSet>("posCalcParameters");
00198   PositionCalc posCalculator_ = PositionCalc(posCalcParameters);
00199   //
00200   std::map<DetId, EcalRecHit> recHitsEB_map;
00201   //
00202   std::vector<EcalRecHit> seeds;
00203 
00204   seeds.clear();
00205   //
00206   vector<EBDetId> usedXtals;
00207   usedXtals.clear();
00208   //
00209   EcalRecHitCollection::const_iterator itb;
00210   //
00211   static const int MAXCLUS = 2000;
00212   int nClus=0;
00213   vector<float> eClus;
00214   vector<float> etClus;
00215   vector<float> etaClus;
00216   vector<float> phiClus;
00217   vector<EBDetId> max_hit;
00218   vector< vector<EcalRecHit> > RecHitsCluster;
00219   vector<float> s4s9Clus;
00220 
00221   // find cluster seeds in EB
00222   for(itb=rhEB->begin(); itb!=rhEB->end(); ++itb){
00223     EBDetId id(itb->id());
00224     double energy = itb->energy();
00225     if (energy > seleXtalMinEnergy_) {
00226       std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
00227       recHitsEB_map.insert(map_entry);
00228     }
00229     if (energy > clusSeedThr_) seeds.push_back(*itb);
00230   } // Eb rechits
00231 
00232   sort(seeds.begin(), seeds.end(), ecalRecHitLess());
00233   for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
00234     EBDetId seed_id = itseed->id();
00235     std::vector<EBDetId>::const_iterator usedIds;
00236 
00237     bool seedAlreadyUsed=false;
00238     for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
00239       if(*usedIds==seed_id){
00240         seedAlreadyUsed=true;
00241         //cout<< " Seed with energy "<<itseed->energy()<<" was used !"<<endl;
00242         break;
00243       }
00244     }
00245     if(seedAlreadyUsed)continue;
00246     topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal,EcalBarrel);
00247     std::vector<DetId> clus_v = topology_p->getWindow(seed_id,clusEtaSize_,clusPhiSize_);
00248     //std::vector<DetId> clus_used;
00249     std::vector<std::pair<DetId, float> > clus_used;
00250 
00251     vector<EcalRecHit> RecHitsInWindow;
00252 
00253     double simple_energy = 0;
00254 
00255     for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
00256       // EBDetId EBdet = *det;
00257       //      cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi "<<EBdet.iphi()<<endl;
00258       bool  HitAlreadyUsed=false;
00259       for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
00260         if(*usedIds==*det){
00261           HitAlreadyUsed=true;
00262           break;
00263         }
00264       }
00265       if(HitAlreadyUsed)continue;
00266       if (recHitsEB_map.find(*det) != recHitsEB_map.end()){
00267         //      cout<<" Used det "<< EBdet<<endl;
00268         std::map<DetId, EcalRecHit>::iterator aHit;
00269         aHit = recHitsEB_map.find(*det);
00270         usedXtals.push_back(*det);
00271         RecHitsInWindow.push_back(aHit->second);
00272         clus_used.push_back( std::pair<DetId, float>(*det, 1.) );
00273         simple_energy = simple_energy + aHit->second.energy();
00274       }
00275     }
00276 
00277     math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_p,geometry_p,geometryES_p);
00278     float theta_s = 2. * atan(exp(-clus_pos.eta()));
00279     float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
00280     float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
00281     //      float p0z_s = simple_energy * cos(theta_s);
00282     float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
00283 
00284     //cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<" "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
00285 
00286     eClus.push_back(simple_energy);
00287     etClus.push_back(et_s);
00288     etaClus.push_back(clus_pos.eta());
00289     phiClus.push_back(clus_pos.phi());
00290     max_hit.push_back(seed_id);
00291     RecHitsCluster.push_back(RecHitsInWindow);
00292     //Compute S4/S9 variable
00293     //We are not sure to have 9 RecHits so need to check eta and phi:
00294     float s4s9_[4];
00295     for(int i=0;i<4;i++)s4s9_[i]= itseed->energy();
00296     for(unsigned int j=0; j<RecHitsInWindow.size();j++){
00297       //cout << " Simple cluster rh, ieta, iphi : "<<((EBDetId)RecHitsInWindow[j].id()).ieta()<<" "<<((EBDetId)RecHitsInWindow[j].id()).iphi()<<endl;
00298       if((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()-1 && seed_id.ieta()!=1 ) || ( seed_id.ieta()==1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()-2))){
00299         if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()-1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()-1 ){
00300           s4s9_[0]+=RecHitsInWindow[j].energy();
00301         }else{
00302           if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()){
00303             s4s9_[0]+=RecHitsInWindow[j].energy();
00304             s4s9_[1]+=RecHitsInWindow[j].energy();
00305           }else{
00306             if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()+1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()+1 ){
00307               s4s9_[1]+=RecHitsInWindow[j].energy();
00308             }
00309           }
00310         }
00311       }else{
00312         if(((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()){
00313           if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()-1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()-1 ){
00314             s4s9_[0]+=RecHitsInWindow[j].energy();
00315             s4s9_[3]+=RecHitsInWindow[j].energy();
00316           }else{
00317             if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()+1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()+1 ){
00318               s4s9_[1]+=RecHitsInWindow[j].energy();
00319               s4s9_[2]+=RecHitsInWindow[j].energy();
00320             }
00321           }
00322         }else{
00323           if((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()+1 && seed_id.ieta()!=-1 ) || ( seed_id.ieta()==-1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()+2))){
00324             if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()-1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()-1 ){
00325               s4s9_[3]+=RecHitsInWindow[j].energy();
00326             }else{
00327               if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()){
00328                 s4s9_[2]+=RecHitsInWindow[j].energy();
00329                 s4s9_[3]+=RecHitsInWindow[j].energy();
00330               }else{
00331                 if(((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()+1 ||((EBDetId)RecHitsInWindow[j].id()).iphi()-360 == seed_id.iphi()+1 ){
00332                   s4s9_[2]+=RecHitsInWindow[j].energy();
00333                 }
00334               }
00335             }
00336           }else{
00337             cout<<" (EBDetId)RecHitsInWindow[j].id()).ieta() "<<((EBDetId)RecHitsInWindow[j].id()).ieta()<<" seed_id.ieta() "<<seed_id.ieta()<<endl;
00338             cout<<" Problem with S4 calculation "<<endl;return;
00339           }
00340         }
00341       }
00342     }
00343     s4s9Clus.push_back(*max_element( s4s9_,s4s9_+4)/simple_energy);
00344     //    cout<<" s4s9Clus[0] "<<s4s9_[0]/simple_energy<<" s4s9Clus[1] "<<s4s9_[1]/simple_energy<<" s4s9Clus[2] "<<s4s9_[2]/simple_energy<<" s4s9Clus[3] "<<s4s9_[3]/simple_energy<<endl;
00345     //    cout<<" Max "<<*max_element( s4s9_,s4s9_+4)/simple_energy<<endl;
00346     nClus++;
00347     if (nClus == MAXCLUS) return;
00348   }  //  End loop over seed clusters
00349 
00350   // cout<< " Pi0 clusters: "<<nClus<<endl;
00351 
00352   // Selection, based on Simple clustering
00353   //pi0 candidates
00354   static const int MAXPI0S = 200;
00355   int npi0_s=0;
00356 
00357   vector<EBDetId> scXtals;
00358   scXtals.clear();
00359 
00360   if (nClus <= 1) return;
00361   for(Int_t i=0 ; i<nClus ; i++){
00362     for(Int_t j=i+1 ; j<nClus ; j++){
00363       //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"  etClus[j] "<<etClus[j]<<endl;
00364       if( etClus[i]>selePtGammaOne_ && etClus[j]>selePtGammaTwo_ && s4s9Clus[i]>seleS4S9GammaOne_ && s4s9Clus[j]>seleS4S9GammaTwo_){
00365         float theta_0 = 2. * atan(exp(-etaClus[i]));
00366         float theta_1 = 2. * atan(exp(-etaClus[j]));
00367 
00368         float p0x = eClus[i] * sin(theta_0) * cos(phiClus[i]);
00369         float p1x = eClus[j] * sin(theta_1) * cos(phiClus[j]);
00370         float p0y = eClus[i] * sin(theta_0) * sin(phiClus[i]);
00371         float p1y = eClus[j] * sin(theta_1) * sin(phiClus[j]);
00372         float p0z = eClus[i] * cos(theta_0);
00373         float p1z = eClus[j] * cos(theta_1);
00374 
00375         float pt_pi0 = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
00376         //      cout<<" pt_pi0 "<<pt_pi0<<endl;
00377         if (pt_pi0 < selePtPi0_)continue;
00378         float m_inv = sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
00379         if ( (m_inv<seleMinvMaxPi0_) && (m_inv>seleMinvMinPi0_) ){
00380 
00381           //New Loop on cluster to measure isolation:
00382           vector<int> IsoClus;
00383           IsoClus.clear();
00384           float Iso = 0;
00385           TVector3 pi0vect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
00386           for(Int_t k=0 ; k<nClus ; k++){
00387             if(k==i || k==j)continue;
00388             TVector3 Clusvect = TVector3(eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * cos(phiClus[k]), eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * sin(phiClus[k]) , eClus[k] * cos(2. * atan(exp(-etaClus[k]))));
00389             float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
00390             float drclpi0 = Clusvect.DeltaR(pi0vect);
00391 
00392             if((drclpi0<selePi0BeltDR_) && (dretaclpi0<selePi0BeltDeta_) ){
00393 
00394               Iso = Iso + etClus[k];
00395               IsoClus.push_back(k);
00396             }
00397           }
00398 
00399 
00400           if(Iso/pt_pi0<selePi0Iso_){
00401 
00402             hMinvPi0EB_->Fill(m_inv);
00403             hPt1Pi0EB_->Fill(etClus[i]);
00404             hPt2Pi0EB_->Fill(etClus[j]);
00405             hPtPi0EB_->Fill(pt_pi0);
00406             hIsoPi0EB_->Fill(Iso/pt_pi0);
00407 
00408 
00409             npi0_s++;
00410           }
00411 
00412           if(npi0_s == MAXPI0S) return;
00413         }
00414       }
00415     }
00416   }
00417 
00418 }
00419 
00420 
00421 
00422 void PiZeroAnalyzer::endJob()
00423 {
00424 
00425 
00426 
00427   bool outputMEsInRootFile = parameters_.getParameter<bool>("OutputMEsInRootFile");
00428   std::string outputFileName = parameters_.getParameter<std::string>("OutputFileName");
00429   if(outputMEsInRootFile){
00430     dbe_->save(outputFileName);
00431   }
00432 
00433   edm::LogInfo("PiZeroAnalyzer") << "Analyzed " << nEvt_  << "\n";
00434   return ;
00435 }
00436 
00437
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DQMOffline/EGamma/plugins/PiZeroAnalyzer.cc
