35 theMaxLoops( iConfig.getUntrackedParameter<unsigned int>(
"maxLoops",0) ),
36 ecalHitsProducer_( iConfig.getParameter< std::string > (
"ecalRecHitsProducer") ),
37 barrelHits_( iConfig.getParameter< std::string > (
"barrelHitCollection") )
40 std::cout <<
"[Pi0FixedMassWindowCalibration] Constructor "<<std::endl;
42 std::string verbosityString = iConfig.
getParameter<std::string>(
"VerbosityLevel");
52 double barrelSeedThreshold = iConfig.
getParameter<
double>(
"IslandBarrelSeedThr");
53 double endcapSeedThreshold = iConfig.
getParameter<
double>(
"IslandEndcapSeedThr");
122 std::cout <<
"[Pi0FixedMassWindowCalibration] endOfJob"<<endl;
128 std::vector<DetId>::const_iterator barrelIt=
barrelCells.begin();
131 int ieta = eb.
ieta();
132 int iphi = eb.
iphi();
133 int sign = eb.
zside()>0 ? 1 : 0;
144 for (
int sign=0; sign<2; sign++) {
145 for (
int ieta=0; ieta<85; ieta++) {
146 for (
int iphi=0; iphi<360; iphi++) {
147 wxtals[ieta][iphi][sign]=0.;
152 std::cout <<
"[Pi0FixedMassWindowCalibration] Starting loop number " << iLoop<<std::endl;
163 std::cout <<
"[Pi0FixedMassWindowCalibration] Ending loop " << iLoop<<std::endl;
165 for (
int sign=0; sign<2; sign++) {
166 for (
int ieta=0; ieta<85; ieta++) {
167 for (
int iphi=0; iphi<360; iphi++) {
169 if (
wxtals[ieta][iphi][sign] == 0 )
175 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;
183 std::vector<DetId>::const_iterator barrelIt=
barrelCells.begin();
186 int ieta = eb.
ieta();
187 int iphi = eb.
iphi();
188 int sign = eb.
zside()>0 ? 1 : 0;
191 std::cout <<
"Calib constant for barrel crystal "
192 <<
" (" << ieta <<
"," << iphi <<
") changed from "
200 for (
int sign=0; sign<2; sign++) {
201 for (
int ieta=0; ieta<85; ieta++) {
202 for (
int iphi=0; iphi<360; iphi++) {
230 for (
int sign=0; sign<2; sign++) {
231 for (
int ieta=0; ieta<85; ieta++) {
232 for (
int iphi=0; iphi<360; iphi++) {
235 wxtals[ieta][iphi][sign]=0.;
248 std::cout <<
"Taken EcalIntercalibConstants" << std::endl;
249 imap = pIcal.
product()->getMap();
252 std::cerr <<
"Error! can't get EcalIntercalibConstants " << std::endl;
263 std::vector<DetId>::const_iterator barrelIt;
269 if ( itcalib == imap.
end() ) {
273 int sign = eb.
zside()>0 ? 1 : 0;
275 if (eb.
iphi()==1)
std::cout <<
"Read old constant for crystal "
276 <<
" (" << eb.
ieta() <<
"," << eb.
iphi()
277 <<
") : " << calib << std::endl;
294 std::cout <<
"[Pi0FixedMassWindowCalibration] Events processed: "<<
nevent<<std::endl;
309 EBDetId ebrhdetid = aRecHitEB->detid();
314 int sign = ebrhdetid.
zside()>0 ? 1 : 0;
316 recalibEcalRecHitCollection->push_back(aHit);
332 std::pair<DetId, EcalRecHit> map_entry(aRecHitEB->id(), *aRecHitEB);
348 std::string clustershapetag;
385 std::vector <reco::ClusterShape> ClusVec_recalib;
386 for (
int erg=0;erg<int(clusters_recalib.size());++erg){
388 ClusVec_recalib.push_back(TestShape_recalib);
393 clustersshapes_p_recalib->assign(ClusVec_recalib.begin(), ClusVec_recalib.end());
395 cout<<
"[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster collection size: "<<clusters_recalib.size()<<endl;
396 cout<<
"[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster Shape Collection size: "<<clustersshapes_p_recalib->size()<<endl;
403 static const int MAXBCEB = 200;
404 static const int MAXBCEBRH = 200;
406 float eIslandBCEB[MAXBCEB];
407 float etIslandBCEB[MAXBCEB];
408 float etaIslandBCEB[MAXBCEB];
409 float phiIslandBCEB[MAXBCEB];
410 int ietaIslandBCEB[MAXBCEB];
411 int iphiIslandBCEB[MAXBCEB];
412 float e2x2IslandBCEB[MAXBCEB];
413 float e3x3IslandBCEB[MAXBCEB];
414 float e5x5IslandBCEB[MAXBCEB];
417 int nIslandBCEBRecHits[MAXBCEB];
419 int ietaIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
420 int iphiIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
421 int zsideIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
422 float eIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
425 for(
int i=0;
i<MAXBCEB;
i++){
428 etaIslandBCEB[
i] = 0;
429 phiIslandBCEB[
i] = 0;
430 ietaIslandBCEB[
i] = 0;
431 iphiIslandBCEB[
i] = 0;
432 e2x2IslandBCEB[
i] = 0;
433 e3x3IslandBCEB[
i] = 0;
434 e5x5IslandBCEB[
i] = 0;
435 nIslandBCEBRecHits[
i] = 0;
436 for(
int j=0;
j<MAXBCEBRH;
j++){
438 ietaIslandBCEBRecHits[
i][
j] = 0;
439 iphiIslandBCEBRecHits[
i][
j] = 0;
440 zsideIslandBCEBRecHits[
i][
j] = 0;
441 eIslandBCEBRecHits[
i][
j] = 0;
447 for(reco::BasicClusterCollection::const_iterator aClus = clusters_recalib.begin(); aClus != clusters_recalib.end(); aClus++) {
448 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;
450 eIslandBCEB[nIslandBCEB] = aClus->energy();
451 etIslandBCEB[nIslandBCEB] = aClus->energy()*
sin(aClus->position().theta());
452 etaIslandBCEB[nIslandBCEB] = aClus->position().eta();
453 phiIslandBCEB[nIslandBCEB] = aClus->position().phi();
455 e2x2IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e2x2();
456 e3x3IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e3x3();
457 e5x5IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e5x5();
459 nIslandBCEBRecHits[nIslandBCEB] = aClus->size();
461 std::vector<std::pair< DetId,float> > hits = aClus->hitsAndFractions();
462 std::vector<std::pair< DetId,float> >::iterator
hit;
463 std::map<DetId, EcalRecHit>::iterator aHit;
466 for(hit = hits.begin(); hit != hits.end(); hit++)
472 EBDetId sel_rh = aHit->second.detid();
477 ietaIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.
ieta();
478 iphiIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.
iphi();
479 zsideIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.
zside();
480 eIslandBCEBRecHits[nIslandBCEB][irhcount] = aHit->second.energy();
494 for(
int i=0 ;
i<nIslandBCEB ;
i++)
496 for(
int j=
i+1 ;
j<nIslandBCEB ;
j++)
502 float theta_0 = 2. * atan(
exp(-etaIslandBCEB[
i]));
503 float theta_1 = 2. * atan(
exp(-etaIslandBCEB[
j]));
505 float p0x = eIslandBCEB[
i] *
sin(theta_0) *
cos(phiIslandBCEB[i]);
506 float p1x = eIslandBCEB[
j] *
sin(theta_1) *
cos(phiIslandBCEB[j]);
508 float p0y = eIslandBCEB[
i] *
sin(theta_0) *
sin(phiIslandBCEB[i]);
509 float p1y = eIslandBCEB[
j] *
sin(theta_1) *
sin(phiIslandBCEB[j]);
511 float p0z = eIslandBCEB[
i] *
cos(theta_0);
512 float p1z = eIslandBCEB[
j] *
cos(theta_1);
514 float pi0_px = p0x + p1x;
515 float pi0_py = p0y + p1y;
516 float pi0_pz = p0z + p1z;
518 float pi0_ptot =
sqrt (pi0_px*pi0_px + pi0_py*pi0_py + pi0_pz*pi0_pz);
520 float pi0_theta = acos(pi0_pz/pi0_ptot);
521 float pi0_eta = -
log(
tan(pi0_theta/2));
522 float pi0_phi = atan(pi0_py/pi0_px);
528 for(Int_t
k=0 ;
k<nIslandBCEB ;
k++)
530 if( (
k != i) && (
k != j) )
532 float dr_pi0_k =
sqrt( (etaIslandBCEB[
k]-pi0_eta)*(etaIslandBCEB[
k]-pi0_eta) + (phiIslandBCEB[
k]-pi0_phi)*(phiIslandBCEB[
k]-pi0_phi) );
533 float deta_pi0_k = fabs(etaIslandBCEB[
k]-pi0_eta);
539 float pt_pi0 =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
545 float m_inv =
sqrt ( (eIslandBCEB[i] + eIslandBCEB[j])*(eIslandBCEB[i] + eIslandBCEB[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
546 cout <<
" pi0 (pt>2.5 GeV) m_inv = "<<m_inv<<endl;
548 float s4s9_1 = e2x2IslandBCEB[
i]/e3x3IslandBCEB[
i];
549 float s4s9_2 = e2x2IslandBCEB[
j]/e3x3IslandBCEB[
j];
551 float s9s25_1 = e3x3IslandBCEB[
i]/e5x5IslandBCEB[
i];
552 float s9s25_2 = e3x3IslandBCEB[
j]/e5x5IslandBCEB[
j];
563 cout<<
" Pi0 Good candidate : minv = "<<m_inv<<endl;
564 for(
int kk=0 ; kk<nIslandBCEBRecHits[
i] ; kk++)
566 int ieta_xtal = ietaIslandBCEBRecHits[
i][kk];
567 int iphi_xtal = iphiIslandBCEBRecHits[
i][kk];
568 int sign = zsideIslandBCEBRecHits[
i][kk]>0 ? 1 : 0;
569 wxtals[
abs(ieta_xtal)-1][iphi_xtal-1][sign] =
wxtals[
abs(ieta_xtal)-1][iphi_xtal-1][sign] + eIslandBCEBRecHits[
i][kk]/e3x3IslandBCEB[
i];
571 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;
574 for(
int kk=0 ; kk<nIslandBCEBRecHits[
j] ; kk++)
576 int ieta_xtal = ietaIslandBCEBRecHits[
j][kk];
577 int iphi_xtal = iphiIslandBCEBRecHits[
j][kk];
578 int sign = zsideIslandBCEBRecHits[
j][kk]>0 ? 1 : 0;
579 wxtals[
abs(ieta_xtal)-1][iphi_xtal-1][sign] =
wxtals[
abs(ieta_xtal)-1][iphi_xtal-1][sign] + eIslandBCEBRecHits[
j][kk]/e3x3IslandBCEB[
j];
581 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;
595 cout<<
" Not enough ECAL Barrel Basic Clusters: "<<nIslandBCEB<<endl;
double selePi0PtGammaTwoMin_
virtual void beginOfJob()
Called at beginning of job.
T getParameter(std::string const &) const
edm::ParameterSet theParameterSet
std::string clustershapecollectionEB_
double selePi0S9S25GammaOneMin_
virtual void startingNewLoop(unsigned int iLoop)
Called at beginning of loop.
void writeLine(EBDetId const &, float)
Sin< T >::type sin(const T &t)
std::vector< DetId > barrelCells
double selePi0S4S9GammaOneMin_
std::vector< T >::const_iterator const_iterator
std::vector< ClusterShape > ClusterShapeCollection
collection of ClusterShape objects
Exp< T >::type exp(const T &t)
IslandClusterAlgo::VerbosityLevel verbosity
PositionCalc posCalculator_
int iphi() const
get the crystal iphi
uint32_t rawId() const
get the raw id
double selePi0PtGammaOneMin_
virtual Status duringLoop(const edm::Event &, const edm::EventSetup &)
Called at each event.
reco::ClusterShape Calculate(const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry, const CaloSubdetectorTopology *topology)
MVATrainerComputer * calib
double selePi0MinvMeanFixed_
virtual Status endOfLoop(const edm::EventSetup &, unsigned int iLoop)
Called at end of loop.
std::string barrelClusterCollection_
Cos< T >::type cos(const T &t)
std::vector< reco::BasicCluster > makeClusters(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry, const CaloSubdetectorTopology *topology_p, const CaloSubdetectorGeometry *geometryES_p, EcalPart ecalPart, bool regional=false, const std::vector< EcalEtaPhiRegion > ®ions=std::vector< EcalEtaPhiRegion >())
Tan< T >::type tan(const T &t)
double selePi0MinvSigmaFixed_
double selePi0S9S25GammaTwoMin_
int ieta() const
get the crystal ieta
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
const EcalRecHitCollection * ecalRecHitBarrelCollection
const_iterator end() const
double newCalibs_barl[85][360][2]
double selePi0EtBeltIsoRatioMax_
Log< T >::type log(const T &t)
virtual void endOfJob()
Called at end of job.
~Pi0FixedMassWindowCalibration()
Destructor.
std::vector< Item >::const_iterator const_iterator
T const * product() const
double selePi0S4S9GammaTwoMin_
std::vector< BasicCluster > BasicClusterCollection
collection of BasicCluster objects
std::string barrelClusterShapeAssociation_
IslandClusterAlgo * island_p
std::vector< DetId > getValidDetIds() const
Get the list of all valid detector ids.
ESHandle< TrackerGeometry > geometry
double mwxtals[85][360][2]
std::string ecalHitsProducer_
const_iterator find(uint32_t rawId) const
double oldCalibs_barl[85][360][2]
const_iterator end() const
std::map< DetId, EcalRecHit > * recHitsEB_map
Pi0FixedMassWindowCalibration(const edm::ParameterSet &iConfig)
Constructor.
const EcalRecHitCollection * recalibEcalRecHitCollection
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
const_iterator begin() const
ClusterShapeAlgo shapeAlgo_
int zside() const
get the z-side of the crystal (1/-1)
float EcalIntercalibConstant
double wxtals[85][360][2]