CMSSW Reference Manual

AlCaPi0BasicClusterRecHitsProducer::AlCaPi0BasicClusterRecHitsProducer ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 8 of file AlCaPi0BasicClusterRecHitsProducer.cc.

References barrelHits_, ecalHitsProducer_, gammaCandEtaSize_, gammaCandPhiSize_, edm::ParameterSet::getParameter(), islandBCColl_, islandBCProd_, pi0BarrelHits_, seleMinvMaxPi0_, seleMinvMinPi0_, selePtGammaOne_, selePtGammaTwo_, and selePtPi0_.

00009 {
00010   ecalHitsProducer_ = iConfig.getParameter< std::string > ("ecalRecHitsProducer");
00011   barrelHits_ = iConfig.getParameter< std::string > ("barrelHitCollection");
00012   pi0BarrelHits_ = iConfig.getParameter< std::string > ("pi0BarrelHitCollection");
00013   islandBCProd_ = iConfig.getParameter< std::string > ("islandBCProd");
00014   islandBCColl_ = iConfig.getParameter< std::string > ("islandBCColl");
00015 
00016   gammaCandEtaSize_ = iConfig.getParameter<int> ("gammaCandEtaSize");
00017   gammaCandPhiSize_ = iConfig.getParameter<int> ("gammaCandPhiSize");
00018   if ( gammaCandPhiSize_ % 2 == 0 ||  gammaCandEtaSize_ % 2 == 0)
00019     edm::LogError("AlCaPi0BasicClusterRecHitsProducerError") << "Size of eta/phi should be odd numbers";
00020 
00021   selePtGammaOne_ = iConfig.getParameter<double> ("selePtGammaOne");  
00022   selePtGammaTwo_ = iConfig.getParameter<double> ("selePtGammaTwo");  
00023   selePtPi0_ = iConfig.getParameter<double> ("selePtPi0");  
00024   seleMinvMaxPi0_ = iConfig.getParameter<double> ("seleMinvMaxPi0");  
00025   seleMinvMinPi0_ = iConfig.getParameter<double> ("seleMinvMinPi0");  
00026 
00027 
00028   //register your products
00029   produces< EBRecHitCollection >(pi0BarrelHits_);
00030 }

AlCaPi0BasicClusterRecHitsProducer::~AlCaPi0BasicClusterRecHitsProducer ( )

Definition at line 33 of file AlCaPi0BasicClusterRecHitsProducer.cc.

References GenMuonPlsPt100GeV_cfg::cout, TimingReport::current(), and TimingReport::dump().

00034 {
00035  
00036   TimingReport::current()->dump(std::cout);
00037 
00038 }

void AlCaPi0BasicClusterRecHitsProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)			`[virtual]`

Implements edm::EDProducer.

Definition at line 43 of file AlCaPi0BasicClusterRecHitsProducer.cc.

References barrelHits_, category, TimerStack::clear_stack(), funct::cos(), GenMuonPlsPt100GeV_cfg::cout, DetId::Ecal, EcalBarrel, ecalHitsProducer_, lat::endl(), EBDetId::ETAPHIMODE, funct::exp(), gammaCandEtaSize_, gammaCandPhiSize_, edm::Event::getByLabel(), i, EBDetId::ieta(), EBDetId::iphi(), islandBCColl_, islandBCProd_, j, LogDebug, DetId::null(), pi0BarrelHits_, TimerStack::pop_and_push(), TimerStack::push(), edm::Event::put(), recHitsEB_map, row, edm::second(), seleMinvMaxPi0_, seleMinvMinPi0_, selePtGammaOne_, selePtGammaTwo_, selePtPi0_, funct::sin(), funct::sqrt(), std, and theta.

00044 {
00045   using namespace edm;
00046   using namespace std;
00047 
00048   // Timer
00049   const std::string category = "AlCaPi0BasicClusterRecHitsProducer";
00050   TimerStack timers;
00051   string timerName = category + "::Total";
00052   timers.push(timerName);
00053 
00054 
00055 
00056   Handle<EBRecHitCollection> barrelRecHitsHandle;
00057 
00058   iEvent.getByLabel(ecalHitsProducer_,barrelHits_,barrelRecHitsHandle);
00059   if (!barrelRecHitsHandle.isValid()) {
00060     LogDebug("") << "AlCaPi0RecHitsProducer: Error! can't get product!" << std::endl;
00061   }
00062 
00063   recHitsEB_map = new std::map<DetId, EcalRecHit>();
00064 
00065 
00066   //Create empty output collections
00067   std::auto_ptr< EBRecHitCollection > pi0EBRecHitCollection( new EBRecHitCollection );
00068 
00069   //Select interesting EcalRecHits (barrel)
00070   EBRecHitCollection::const_iterator itb;
00071   cout<< "   EB RecHits #: "<<barrelRecHitsHandle->size()<<endl;
00072   for (itb=barrelRecHitsHandle->begin(); itb!=barrelRecHitsHandle->end(); itb++) {
00073 
00074     std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
00075     recHitsEB_map->insert(map_entry);
00076 
00077   }
00078 
00079   timerName = category + "::readEBRecHitsCollection";
00080   timers.push(timerName);
00081 
00082 
00083 
00084   // Get ECAL Barrel Island Basic Clusters collection
00085   // ECAL Barrel Island Basic Clusters 
00086      static const int MAXBCEB = 200;
00087      int nIslandBCEB;
00088      float eIslandBCEB[MAXBCEB];
00089      float etIslandBCEB[MAXBCEB];
00090      float etaIslandBCEB[MAXBCEB];
00091      float phiIslandBCEB[MAXBCEB];
00092 
00093   nIslandBCEB=0;
00094   for(int i=0; i<MAXBCEB; i++){
00095     eIslandBCEB[i] = 0;
00096     etIslandBCEB[i] = 0;
00097     etaIslandBCEB[i] = 0;
00098     phiIslandBCEB[i] = 0;
00099   }
00100   //pi0 candidates
00101   static const int MAXPI0 = 20;
00102   int npi0=0;
00103   int BC1[MAXPI0],BC2[MAXPI0];
00104   for(int i=0; i<MAXPI0; i++){
00105     BC1[i]=0;
00106     BC2[i]=0;
00107   }
00108 
00109   Handle<reco::BasicClusterCollection> pIslandBarrelBasicClusters;
00110   iEvent.getByLabel(islandBCProd_, islandBCColl_, pIslandBarrelBasicClusters);
00111   const reco::BasicClusterCollection* islandBarrelBasicClusters = pIslandBarrelBasicClusters.product();
00112   cout << " island Basic Clusters # "<< islandBarrelBasicClusters->size() <<endl;
00113 
00114   for(reco::BasicClusterCollection::const_iterator aClus = islandBarrelBasicClusters->begin();
00115       aClus != islandBarrelBasicClusters->end(); aClus++) {
00116 
00117     //cout << " island Basic Cluster E, #xtals: "<<  aClus->energy()<<" "<< aClus->getHitsByDetId().size() <<endl;
00118     //cout << " island Basic Cluster Position (x,y,z): "<<aClus->position().X()<<" "<<aClus->position().Y()<<" "<<aClus->position().Z()<<" "<<endl;
00119     //cout << " island Basic Cluster Position (eta,phi): "<<aClus->position().eta()<<" "<<aClus->position().phi()<<endl;
00120 
00121     float theta = 2. * atan(exp(-aClus->position().eta()));
00122     float p0x = aClus->energy() * sin(theta) * cos(aClus->position().phi());
00123     float p0y = aClus->energy() * sin(theta) * sin(aClus->position().phi());
00124     float p0z = aClus->energy() * cos(theta);
00125     float et = sqrt( p0x*p0x + p0y*p0y);
00126 
00127     //cout << " island Basic Cluster E,Et,px,py,pz: "<<aClus->energy()<<" "<<et<<" "<<p0x<<" "<<p0y<<" "<<p0z<<endl; 
00128 
00129     eIslandBCEB[nIslandBCEB] = aClus->energy();
00130     etIslandBCEB[nIslandBCEB] = et;
00131     etaIslandBCEB[nIslandBCEB] = aClus->position().eta();
00132     phiIslandBCEB[nIslandBCEB] = aClus->position().phi();
00133      
00134     nIslandBCEB++;
00135   }
00136 
00137   timerName = category + "::readIslandBasicClustersCollection";
00138   timers.push(timerName);
00139 
00140 
00141 
00142   // Selection, based on ECAL Barrel Basic Clusters
00143 
00144   if (nIslandBCEB > 1) 
00145     {
00146       for(Int_t i=0 ; i<nIslandBCEB ; i++)
00147         {
00148           for(Int_t j=i+1 ; j<nIslandBCEB ; j++)
00149             {
00150 
00151               if( etIslandBCEB[i]>selePtGammaOne_ && etIslandBCEB[j]>selePtGammaTwo_) 
00152                 {
00153                 
00154                   float theta_0 = 2. * atan(exp(-etaIslandBCEB[i]));
00155                   float theta_1 = 2. * atan(exp(-etaIslandBCEB[j]));
00156                 
00157                   float p0x = eIslandBCEB[i] * sin(theta_0) * cos(phiIslandBCEB[i]);
00158                   float p1x = eIslandBCEB[j] * sin(theta_1) * cos(phiIslandBCEB[j]);
00159                 
00160                   float p0y = eIslandBCEB[i] * sin(theta_0) * sin(phiIslandBCEB[i]);
00161                   float p1y = eIslandBCEB[j] * sin(theta_1) * sin(phiIslandBCEB[j]);
00162                   float p0z = eIslandBCEB[i] * cos(theta_0);
00163                   float p1z = eIslandBCEB[j] * cos(theta_1);
00164                 
00165                   //
00166                 
00167                   float pt_pi0 = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
00168                   //float dr_pi0 = sqrt ( (etaIslandBCEB[i]-etaIslandBCEB[j])*(etaIslandBCEB[i]-etaIslandBCEB[j]) + (phiIslandBCEB[i]-phiIslandBCEB[j])*(phiIslandBCEB[i]-phiIslandBCEB[j]) );
00169                   if (pt_pi0 > selePtPi0_) 
00170                     {
00171                       float m_inv = sqrt ( (eIslandBCEB[i] + eIslandBCEB[j])*(eIslandBCEB[i] + eIslandBCEB[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
00172                       if ( (m_inv<seleMinvMaxPi0_) && (m_inv>seleMinvMinPi0_) )
00173                         {
00174                           cout <<" pi0 Candidate (pt>2.5 GeV, m_inv<0.2) pt,m_inv,i,j :   "<<pt_pi0<<" "<<m_inv<<" "<<i<<" "<<j<<" "<<endl;  
00175                           BC1[npi0]=i;
00176                           BC2[npi0]=j;
00177                           npi0++;
00178                         }
00179 
00180                     }
00181                 
00182                 }
00183             } // End of the "j" loop over BCEB
00184         } // End of the "i" loop over BCEB
00185 
00186 
00187       cout<<" "<<endl;
00188       cout<<"  Pi0 candidates #: "<<npi0<<endl;
00189 
00190       timerName = category + "::makePi0Cand";
00191       timers.pop_and_push(timerName);
00192 
00193 
00194 
00195       vector<EBDetId> scXtals;
00196       scXtals.clear();
00197       for(Int_t i=0 ; i<npi0 ; i++)
00198         {
00199           // cout<<"     Pi0 i,Bc1, Bc2 "<<i<<" "<<BC1[i]<<" "<<BC2[i]<<endl; 
00200 
00201           int intbc=0;
00202 
00203           for(reco::BasicClusterCollection::const_iterator aClus = islandBarrelBasicClusters->begin();
00204               aClus != islandBarrelBasicClusters->end(); aClus++) {
00205 
00206             //   cout << " intbc, BC1[i], BC2[i] "<<intbc<<" "<<BC1[i]<<" "<<BC2[i]<<endl; 
00207             if((intbc==BC1[i]) || (intbc==BC2[i]))
00208               {
00209 
00210                 std::vector<DetId> hits = aClus->getHitsByDetId();
00211                 //std::vector<DetId>::iterator hit;
00212                 std::map<DetId, EcalRecHit>::iterator aHit;
00213 
00214                 // New addition - all existed xtals inside a sliding window will be added 
00215 
00216                 double currEnergy = 0.;
00217                 EBDetId maxHit(0);
00218 
00219                 for( std::vector<DetId>::const_iterator idsIt = hits.begin(); idsIt != hits.end(); ++idsIt) {
00220           
00221                   if((*idsIt).subdetId()!=EcalBarrel || (*idsIt).det()!= DetId::Ecal) continue;
00222           
00223                   if(((recHitsEB_map->find(*idsIt))->second).energy() > currEnergy) {
00224                     currEnergy=((recHitsEB_map->find(*idsIt))->second).energy();
00225                     maxHit=*idsIt;
00226                   }
00227                   aHit = recHitsEB_map->find(*idsIt);
00228                   pi0EBRecHitCollection->push_back(aHit->second);
00229                   scXtals.push_back(*idsIt);
00230 
00231                   //EBDetId sel_rh_bc = aHit->second.detid();
00232                   // cout << "       RecHit Ok, belongs to cluster # "<< intbc<<" : z,ieta,iphi "<<sel_rh_bc.zside()<<" "<<sel_rh_bc.ieta()<<" "<<sel_rh_bc.iphi()<<endl;    
00233                   //cout << "       RecHit Ok, belongs to cluster # "<< intbc<<" : tower_ieta,tower_iphi "<<sel_rh_bc.tower_ieta()<<" "<<sel_rh_bc.tower_iphi()<<endl;   
00234                   //cout << "       RecHit Ok, belongs to  cluster # "<<intbc <<" : iSM, ic "<<sel_rh_bc.ism()<<" "<<sel_rh_bc.ic()<<endl;
00235 
00236                 }
00237 
00238                 if (!maxHit.null())
00239                   for (int icry=0;icry< gammaCandEtaSize_*gammaCandPhiSize_;icry++)
00240                     {
00241               
00242                       int row = icry / gammaCandEtaSize_;
00243                       int column= icry % gammaCandEtaSize_;
00244                       int curr_eta=maxHit.ieta() + column - (gammaCandEtaSize_/2);
00245                       int curr_phi=maxHit.iphi() + row - (gammaCandPhiSize_/2);
00246               
00247                       if (curr_eta * maxHit.ieta() <= 0) {if (maxHit.ieta() > 0) curr_eta--; else curr_eta++; }  // JUMP over 0
00248                       if (curr_phi < 1) curr_phi += 360;
00249                       if (curr_phi > 360) curr_phi -= 360;
00250               
00251                       try
00252                         {
00253                           EBDetId det = EBDetId(curr_eta,curr_phi,EBDetId::ETAPHIMODE);
00254                           std::vector<EBDetId>::const_iterator usedIds;
00255                           
00256                           bool HitAlreadyUsed=false;
00257                           for(usedIds=scXtals.begin(); usedIds!=scXtals.end(); usedIds++)
00258                             if(*usedIds==det)
00259                               {
00260                                 HitAlreadyUsed=true;
00261                                 break;
00262                               }
00263                           
00264                           if(!HitAlreadyUsed)
00265                             if (recHitsEB_map->find(det) != recHitsEB_map->end())
00266                               {
00267                                 aHit = recHitsEB_map->find(det);
00268                                 pi0EBRecHitCollection->push_back(aHit->second);
00269                                 scXtals.push_back(det);
00270 
00271                                 //EBDetId sel_rh = aHit->second.detid();
00272                                 //cout << "       RecHit Ok 20x20 matrix outside cluster # "<<intbc<<" : z,ieta,iphi "<<sel_rh.zside()<<" "<<sel_rh.ieta()<<" "<<sel_rh.iphi()<<endl;    
00273                                 //cout << "       RecHit Ok 20x20 matrix outside cluster # "<<intbc<<" : tower_ieta,tower_iphi "<<sel_rh.tower_ieta()<<" "<<sel_rh.tower_iphi()<<endl;   
00274                                 //cout << "       RecHit Ok 20x20 matrix outside cluster # "<<intbc<<" : iSM, ic "<<sel_rh.ism()<<" "<<sel_rh.ic()<<endl;
00275 
00276                               }
00277                         }
00278                       catch (...)
00279                         {
00280                         }
00281                     }
00282 
00283               }
00284             intbc++;
00285                 
00286           }
00287         }
00288     }
00289 
00290   timerName = category + "::preparePi0RecHitsCollection";
00291   timers.pop_and_push(timerName);
00292 
00293 
00294       //Put selected information in the event
00295       //if (npi0>0) iEvent.put( pi0EBRecHitCollection, pi0BarrelHits_);
00296 
00297   cout<< "   EB RecHits # in Collection: "<<pi0EBRecHitCollection->size()<<endl;
00298   iEvent.put( pi0EBRecHitCollection, pi0BarrelHits_);
00299 
00300   timerName = category + "::storePi0RecHitsCollection";
00301   timers.pop_and_push(timerName);
00302   
00303   timers.clear_stack();
00304   delete recHitsEB_map;
00305 
00306 }

std::string AlCaPi0BasicClusterRecHitsProducer::barrelHits_ [private]

Definition at line 45 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

std::string AlCaPi0BasicClusterRecHitsProducer::ecalHitsProducer_ [private]

Definition at line 44 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

int AlCaPi0BasicClusterRecHitsProducer::gammaCandEtaSize_ [private]

Definition at line 50 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

int AlCaPi0BasicClusterRecHitsProducer::gammaCandPhiSize_ [private]

Definition at line 51 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

std::string AlCaPi0BasicClusterRecHitsProducer::islandBCColl_ [private]

Definition at line 48 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

std::string AlCaPi0BasicClusterRecHitsProducer::islandBCProd_ [private]

Definition at line 47 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

std::string AlCaPi0BasicClusterRecHitsProducer::pi0BarrelHits_ [private]

Definition at line 46 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

std::map<DetId, EcalRecHit>* AlCaPi0BasicClusterRecHitsProducer::recHitsEB_map [private]

Definition at line 59 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by produce().

double AlCaPi0BasicClusterRecHitsProducer::seleMinvMaxPi0_ [private]

Definition at line 56 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

double AlCaPi0BasicClusterRecHitsProducer::seleMinvMinPi0_ [private]

Definition at line 57 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

double AlCaPi0BasicClusterRecHitsProducer::selePtGammaOne_ [private]

Definition at line 53 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

double AlCaPi0BasicClusterRecHitsProducer::selePtGammaTwo_ [private]

Definition at line 54 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

double AlCaPi0BasicClusterRecHitsProducer::selePtPi0_ [private]

Definition at line 55 of file AlCaPi0BasicClusterRecHitsProducer.h.

Referenced by AlCaPi0BasicClusterRecHitsProducer(), and produce().

AlCaPi0BasicClusterRecHitsProducer Class Reference

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	AlCaPi0BasicClusterRecHitsProducer (const edm::ParameterSet &)
virtual void	produce (edm::Event &, const edm::EventSetup &)
	~AlCaPi0BasicClusterRecHitsProducer ()
Private Attributes
std::string	barrelHits_
std::string	ecalHitsProducer_
int	gammaCandEtaSize_
int	gammaCandPhiSize_
std::string	islandBCColl_
std::string	islandBCProd_
std::string	pi0BarrelHits_
std::map< DetId, EcalRecHit > *	recHitsEB_map
double	seleMinvMaxPi0_
double	seleMinvMinPi0_
double	selePtGammaOne_
double	selePtGammaTwo_
double	selePtPi0_