CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
ContainmentCorrectionAnalyzer Class Reference

#include <ContainmentCorrectionAnalyzer.h>

Inheritance diagram for ContainmentCorrectionAnalyzer:
edm::EDAnalyzer edm::EDConsumerBase

Public Member Functions

virtual void analyze (const edm::Event &, const edm::EventSetup &)
 
virtual void beginJob ()
 
 ContainmentCorrectionAnalyzer (const edm::ParameterSet &)
 
virtual void endJob ()
 
 ~ContainmentCorrectionAnalyzer ()
 
- Public Member Functions inherited from edm::EDAnalyzer
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
 EDAnalyzer ()
 
ModuleDescription const & moduleDescription () const
 
std::string workerType () const
 
virtual ~EDAnalyzer ()
 
- Public Member Functions inherited from edm::EDConsumerBase
 EDConsumerBase ()
 
ProductHolderIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
std::vector
< ProductHolderIndexAndSkipBit >
const & 
itemsToGetFromEvent () const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
bool registeredToConsume (ProductHolderIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
 
virtual ~EDConsumerBase ()
 

Private Member Functions

float ecalEta (float EtaParticle, float Zvertex, float plane_Radius)
 
void fillMcTruth (std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)
 
std::vector< EcalSimPhotonMCTruthfindMcTruth (std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)
 

Private Attributes

edm::InputTag BarrelSuperClusterCollection_
 
std::vector< float > e1
 
std::vector< float > e25
 
std::vector< float > e9
 
edm::InputTag EndcapSuperClusterCollection_
 
std::map< unsigned, unsigned > geantToIndex_
 
TH1F * h_EB_converted
 
TH1F * h_EB_e25EtrueReference
 
TH1F * h_EB_e9EtrueReference
 
TH1F * h_EB_eRecoEtrueReference
 
TH1F * h_EB_eTrue
 
TH1F * h_EE_converted
 
TH1F * h_EE_e25EtrueReference
 
TH1F * h_EE_e9EtrueReference
 
TH1F * h_EE_eRecoEtrueReference
 
TH1F * h_EE_eTrue
 
std::vector< int > iMC
 
std::vector< int > isConverted
 
std::vector< float > mcEnergy
 
std::vector< float > mcEta
 
std::vector< float > mcPhi
 
std::vector< float > mcPt
 
int nMCphotons
 
int nRECOphotons
 
edm::InputTag reducedBarrelRecHitCollection_
 
edm::InputTag reducedEndcapRecHitCollection_
 
std::vector< int > seedXtal
 
edm::InputTag SimTrackCollection_
 
edm::InputTag SimVertexCollection_
 
std::vector< float > superClusterEnergy
 
std::vector< float > superClusterEt
 
std::vector< float > superClusterEta
 
std::vector< float > superClusterPhi
 
std::vector< float > x_vtx
 
std::vector< float > y_vtx
 
std::vector< float > z_vtx
 

Additional Inherited Members

- Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer ModuleType
 
- Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &)
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Definition at line 47 of file ContainmentCorrectionAnalyzer.h.

Constructor & Destructor Documentation

ContainmentCorrectionAnalyzer::ContainmentCorrectionAnalyzer ( const edm::ParameterSet ps)
explicit

Definition at line 11 of file ContainmentCorrectionAnalyzer.cc.

References edm::ParameterSet::getParameter().

11  {
12 
13  BarrelSuperClusterCollection_ = ps.getParameter<InputTag>("BarrelSuperClusterCollection");
14  EndcapSuperClusterCollection_ = ps.getParameter<InputTag>("EndcapSuperClusterCollection");
15  reducedBarrelRecHitCollection_ = ps.getParameter<InputTag>("reducedBarrelRecHitCollection");
16  reducedEndcapRecHitCollection_ = ps.getParameter<InputTag>("reducedEndcapRecHitCollection");
17  SimTrackCollection_ = ps.getParameter<InputTag>("simTrackCollection");
18  SimVertexCollection_ = ps.getParameter<InputTag>("simVertexCollection");
19 }
T getParameter(std::string const &) const
ContainmentCorrectionAnalyzer::~ContainmentCorrectionAnalyzer ( )

Definition at line 21 of file ContainmentCorrectionAnalyzer.cc.

21 { }

Member Function Documentation

void ContainmentCorrectionAnalyzer::analyze ( const edm::Event evt,
const edm::EventSetup es 
)
virtual

Implements edm::EDAnalyzer.

Definition at line 40 of file ContainmentCorrectionAnalyzer.cc.

References deltaR(), alignCSCRings::e, reco::e1, jptDQMConfig_cff::eMax, eta(), HcalObjRepresent::Fill(), edm::EventSetup::get(), edm::Event::getByLabel(), i, edm::EventBase::id(), edm::ESHandleBase::isValid(), edm::HandleBase::isValid(), interactiveExample::photons, pi, reco::CaloCluster::position(), funct::pow(), reco::SuperCluster::preshowerEnergy(), edm::Handle< T >::product(), edm::ESHandle< class >::product(), reco::SuperCluster::rawEnergy(), reco::SuperCluster::seed(), and mathSSE::sqrt().

40  {
41 
42  LogInfo("ContainmentCorrectionAnalyzer") << "Analyzing event " << evt.id() << "\n";
43 
44  // taking the needed collections
45  std::vector<SimTrack> theSimTracks;
48  if (SimTk.isValid() ) theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
49  else { LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << SimTrackCollection_.label(); }
50 
51  std::vector<SimVertex> theSimVertexes;
53  evt.getByLabel(SimVertexCollection_,SimVtx);
54  if (SimVtx.isValid()) theSimVertexes.insert(theSimVertexes.end(),SimVtx->begin(),SimVtx->end());
55  else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << SimVertexCollection_.label(); }
56 
57  const reco::SuperClusterCollection* BarrelSuperClusters = 0;
58  Handle<reco::SuperClusterCollection> pHybridBarrelSuperClusters;
59  evt.getByLabel(BarrelSuperClusterCollection_, pHybridBarrelSuperClusters);
60  if (pHybridBarrelSuperClusters.isValid()) { BarrelSuperClusters = pHybridBarrelSuperClusters.product(); }
61  else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << BarrelSuperClusterCollection_.label(); }
62 
63  const reco::SuperClusterCollection* EndcapSuperClusters = 0;
64  Handle<reco::SuperClusterCollection> pMulti5x5EndcapSuperClusters;
65  evt.getByLabel(EndcapSuperClusterCollection_, pMulti5x5EndcapSuperClusters);
66  if (pMulti5x5EndcapSuperClusters.isValid()) EndcapSuperClusters = pMulti5x5EndcapSuperClusters.product();
67  else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << EndcapSuperClusterCollection_.label(); }
68 
69  const EcalRecHitCollection *ebRecHits = 0;
71  evt.getByLabel( reducedBarrelRecHitCollection_, pEBRecHits );
72  if (pEBRecHits.isValid()) ebRecHits = pEBRecHits.product();
73  else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << reducedBarrelRecHitCollection_.label(); }
74 
75  const EcalRecHitCollection *eeRecHits = 0;
77  evt.getByLabel( reducedEndcapRecHitCollection_, pEERecHits );
78  if (pEERecHits.isValid()) eeRecHits = pEERecHits.product();
79  else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << reducedEndcapRecHitCollection_.label(); }
80 
81  const CaloTopology *topology = 0;
82  ESHandle<CaloTopology> pTopology;
83  es.get<CaloTopologyRecord>().get(pTopology);
84  if(pTopology.isValid()) topology = pTopology.product();
85 
86  std::vector<EcalSimPhotonMCTruth> photons=findMcTruth(theSimTracks,theSimVertexes);
87 
88  nMCphotons = 0;
89  nRECOphotons = 0;
90 
91  int mc_size = photons.size();
92  mcEnergy.resize(mc_size);
93  mcEta.resize(mc_size);
94  mcPhi.resize(mc_size);
95  mcPt.resize(mc_size);
96  isConverted.resize(mc_size);
97  x_vtx.resize(mc_size);
98  y_vtx.resize(mc_size);
99  z_vtx.resize(mc_size);
100 
101  superClusterEnergy.resize(mc_size);
102  superClusterEta.resize(mc_size);
103  superClusterPhi.resize(mc_size);
104  superClusterEt.resize(mc_size);
105  e1.resize(mc_size);
106  e9.resize(mc_size);
107  e25.resize(mc_size);
108  seedXtal.resize(mc_size);
109  iMC.resize(mc_size);
110 
111 
112  // loop over MC truth photons
113  for (unsigned int ipho=0;ipho<photons.size();ipho++) {
114 
115  math::XYZTLorentzVectorD vtx = photons[ipho].primaryVertex();
116  double phiTrue = photons[ipho].fourMomentum().phi();
117  double vtxPerp = sqrt(vtx.x()*vtx.x() + vtx.y()*vtx.y());
118  double etaTrue = ecalEta(photons[ipho].fourMomentum().eta(), vtx.z(), vtxPerp);
119  double etTrue = photons[ipho].fourMomentum().e()/cosh(etaTrue);
120  nMCphotons++;
121  mcEnergy[nMCphotons-1]=photons[ipho].fourMomentum().e();
122  mcEta[nMCphotons-1]=etaTrue;
123  mcPhi[nMCphotons-1]=phiTrue;
124  mcPt[nMCphotons-1]=etTrue;
125  isConverted[nMCphotons-1]=photons[ipho].isAConversion();
126  x_vtx[nMCphotons-1]=vtx.x();
127  y_vtx[nMCphotons-1]=vtx.y();
128  z_vtx[nMCphotons-1]=vtx.z();
129 
130  // check histos for MC truth
131  if(std::fabs(etaTrue) < 1.479) {
132  h_EB_eTrue ->Fill(photons[ipho].fourMomentum().e());
133  h_EB_converted ->Fill(photons[ipho].isAConversion());
134  }
135  if(std::fabs(etaTrue) >= 1.6) {
136  h_EE_eTrue ->Fill(photons[ipho].fourMomentum().e());
137  h_EE_converted ->Fill(photons[ipho].isAConversion());
138  }
139 
140  // barrel
141  if(std::fabs(etaTrue) < 1.479) {
142  double etaCurrent; // , etaFound = 0; // UNUSED
143  double phiCurrent;
144  // double etCurrent, etFound = 0; // UNUSED
145  const reco::SuperCluster* nearSC = 0;
146 
147  double closestParticleDistance = 999;
148  for(reco::SuperClusterCollection::const_iterator aClus = BarrelSuperClusters->begin();
149  aClus != BarrelSuperClusters->end(); aClus++) {
150  etaCurrent = aClus->position().eta();
151  phiCurrent = aClus->position().phi();
152  // etCurrent = aClus->energy()/std::cosh(etaCurrent); // UNUSED
153  double deltaR = std::sqrt(std::pow(etaCurrent-etaTrue,2)+std::pow(phiCurrent-phiTrue,2));
154  if(deltaR < closestParticleDistance) {
155  // etFound = etCurrent; // UNUSED
156  // etaFound = etaCurrent; // UNUSED
157  closestParticleDistance = deltaR;
158  nearSC=&(*aClus);
159  }
160  }
161 
162  if(closestParticleDistance < 0.3) {
163  // Is a matched particle dumping informations
164  nRECOphotons++;
166  superClusterEta[nRECOphotons-1]=nearSC->position().eta();
167  superClusterPhi[nRECOphotons-1]=nearSC->position().phi();
168  superClusterEt[nRECOphotons-1]=nearSC->rawEnergy()/std::cosh(nearSC->position().eta());
170 
171  reco::CaloClusterPtr theSeed = nearSC->seed();
172  seedXtal[nRECOphotons-1] = EcalClusterTools::getMaximum(*theSeed, ebRecHits).first;
173  e1[nRECOphotons-1] = EcalClusterTools::eMax(*theSeed, ebRecHits);
174  e9[nRECOphotons-1] = EcalClusterTools::e3x3(*theSeed, ebRecHits, topology );
175  e25[nRECOphotons-1] = EcalClusterTools::e5x5(*theSeed, ebRecHits, topology );
176  }
177  }
178 
179  // endcap
180  if(std::fabs(etaTrue) >= 1.6) {
181  double etaCurrent; // , etaFound = 0; // UNUSED
182  double phiCurrent;
183  // double etCurrent, etFound = 0; // UNUSED
184  const reco::SuperCluster* nearSC = 0;
185 
186  double closestParticleDistance = 999;
187  for(reco::SuperClusterCollection::const_iterator aClus = EndcapSuperClusters->begin();
188  aClus != EndcapSuperClusters->end(); aClus++) {
189  etaCurrent = aClus->position().eta();
190  phiCurrent = aClus->position().phi();
191  // etCurrent = aClus->energy()/std::cosh(etaCurrent);
192  double deltaR = std::sqrt(std::pow(etaCurrent-etaTrue,2)+std::pow(phiCurrent-phiTrue,2));
193  if(deltaR < closestParticleDistance) {
194  // etFound = etCurrent; // UNUSED
195  // etaFound = etaCurrent; // UNUSED
196  closestParticleDistance = deltaR;
197  nearSC=&(*aClus);
198  }
199  }
200 
201  if(closestParticleDistance < 0.3) {
202  //Is a matched particle dumping informations
203  nRECOphotons++;
204  float psEnergy = nearSC->preshowerEnergy();
205  superClusterEnergy[nRECOphotons-1]=nearSC->rawEnergy() + psEnergy;
206  superClusterEta[nRECOphotons-1]=nearSC->position().eta();
207  superClusterPhi[nRECOphotons-1]=nearSC->position().phi();
208  superClusterEt[nRECOphotons-1]=(nearSC->rawEnergy() + psEnergy)/std::cosh(nearSC->position().eta());
210 
211  reco::CaloClusterPtr theSeed = nearSC->seed();
212  seedXtal[nRECOphotons-1] = EcalClusterTools::getMaximum(*theSeed, eeRecHits).first;
213  e1[nRECOphotons-1] = EcalClusterTools::eMax(*theSeed, eeRecHits) + psEnergy;
214  e9[nRECOphotons-1] = EcalClusterTools::e3x3(*theSeed, eeRecHits, topology ) + psEnergy;
215  e25[nRECOphotons-1] = EcalClusterTools::e5x5(*theSeed, eeRecHits, topology ) + psEnergy;
216  }
217  }
218  }
219 
220 
221  // containment analysis for unconverted photons in the reference region only
222  for (int i=0;i<nRECOphotons;i++) {
223 
224  // barrel
225  if (fabs(superClusterEta[i]) < 1.479 ) {
226  if (isConverted[iMC[i]] != 1){
227  int ietaAbs = (seedXtal[i]>>9)&0x7F;
228  int iphi = seedXtal[i]&0x1FF;
229  if (ietaAbs > 5 && ietaAbs < 21 && ((iphi % 20) > 5) && ((iphi % 20) < 16) ) {
233  }
234  }
235  }
236 
237  // endcap
238  if (fabs(superClusterEta[i]) > 1.6 ) {
239  if (isConverted[iMC[i]] != 1) {
240  if ( fabs(superClusterEta[i]) > 1.7 && fabs(superClusterEta[i] < 2.3) &&
241  ((superClusterPhi[i] > -CLHEP::pi/2. + 0.1 && superClusterPhi[i] < CLHEP::pi/2. - 0.1) ||
242  (superClusterPhi[i] > CLHEP::pi/2. + 0.1) ||
243  (superClusterPhi[i] < -CLHEP::pi/2. - 0.1)
244  ))
245  {
249  }
250  }
251  }
252 
253  } // loop over reco photons
254 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:124
int i
Definition: DBlmapReader.cc:9
std::vector< EcalSimPhotonMCTruth > findMcTruth(std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
T eta() const
float ecalEta(float EtaParticle, float Zvertex, float plane_Radius)
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
T sqrt(T t)
Definition: SSEVec.h:48
bool isValid() const
Definition: HandleBase.h:76
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:390
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
double rawEnergy() const
raw uncorrected energy (sum of energies of component BasicClusters)
Definition: SuperCluster.h:47
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
T const * product() const
Definition: Handle.h:81
std::string const & label() const
Definition: InputTag.h:42
edm::EventID id() const
Definition: EventBase.h:56
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:66
double preshowerEnergy() const
energy deposited in preshower
Definition: SuperCluster.h:50
double pi
bool isValid() const
Definition: ESHandle.h:37
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void ContainmentCorrectionAnalyzer::beginJob ( void  )
virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 23 of file ContainmentCorrectionAnalyzer.cc.

References TFileService::make().

23  {
24 
26 
27  // Define reference histograms
28  h_EB_eRecoEtrueReference = fs->make<TH1F>("EB_eRecoEtrueReference","EB_eRecoEtrueReference",440,0.,1.1);
29  h_EB_e9EtrueReference = fs->make<TH1F>("EB_e9EtrueReference", "EB_e9EtrueReference", 440,0.,1.1);
30  h_EB_e25EtrueReference = fs->make<TH1F>("EB_e25EtrueReference", "EB_e25EtrueReference", 440,0.,1.1);
31  h_EE_eRecoEtrueReference = fs->make<TH1F>("EE_eRecoEtrueReference","EE_eRecoEtrueReference",440,0.,1.1);
32  h_EE_e9EtrueReference = fs->make<TH1F>("EE_e9EtrueReference", "EE_e9EtrueReference", 440,0.,1.1);
33  h_EE_e25EtrueReference = fs->make<TH1F>("EE_e25EtrueReference", "EE_e25EtrueReference", 440,0.,1.1);
34  h_EB_eTrue = fs->make<TH1F>("EB_eTrue", "EB_eTrue", 41,40.,60.);
35  h_EE_eTrue = fs->make<TH1F>("EE_eTrue", "EE_eTrue", 41,40.,60.);
36  h_EB_converted = fs->make<TH1F>("EB_converted", "EB_converted", 2,0.,2.);
37  h_EE_converted = fs->make<TH1F>("EE_converted", "EE_converted", 2,0.,2.);
38 }
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
float ContainmentCorrectionAnalyzer::ecalEta ( float  EtaParticle,
float  Zvertex,
float  plane_Radius 
)
private

Definition at line 258 of file ContainmentCorrectionAnalyzer.cc.

References ETA, etaBarrelEndcap, create_public_lumi_plots::log, Geom::pi(), R_ECAL, funct::tan(), and Z_Endcap.

258  {
259 
260  const float R_ECAL = 136.5;
261  const float Z_Endcap = 328.0;
262  const float etaBarrelEndcap = 1.479;
263 
264  if(EtaParticle != 0.) {
265  float Theta = 0.0 ;
266  float ZEcal = (R_ECAL-plane_Radius)*sinh(EtaParticle)+Zvertex;
267 
268  if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
269  if(Theta<0.0) Theta = Theta+Geom::pi() ;
270 
271  float ETA = - log(tan(0.5*Theta));
272 
273  if( fabs(ETA) > etaBarrelEndcap ) {
274  float Zend = Z_Endcap ;
275  if(EtaParticle<0.0 ) Zend = -Zend ;
276  float Zlen = Zend - Zvertex ;
277  float RR = Zlen/sinh(EtaParticle);
278  Theta = atan((RR+plane_Radius)/Zend);
279  if(Theta<0.0) Theta = Theta+Geom::pi() ;
280  ETA = - log(tan(0.5*Theta));
281  }
282 
283  return ETA;
284  }
285  else {
286  LogWarning("") << "[ContainmentCorrectionAnalyzer::ecalEta] Warning: Eta equals to zero, not correcting" ;
287  return EtaParticle;
288  }
289 }
#define ETA
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
static const float etaBarrelEndcap
static const float Z_Endcap
double pi()
Definition: Pi.h:31
static const float R_ECAL
void ContainmentCorrectionAnalyzer::endJob ( void  )
virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 256 of file ContainmentCorrectionAnalyzer.cc.

256 { }
void ContainmentCorrectionAnalyzer::fillMcTruth ( std::vector< SimTrack > &  theSimTracks,
std::vector< SimVertex > &  theSimVertices 
)
private

Definition at line 424 of file ContainmentCorrectionAnalyzer.cc.

424  {
425 
426  unsigned nVtx = simVertices.size();
427  unsigned nTks = simTracks.size();
428  if ( nVtx == 0 ) return;
429  // create a map associating geant particle id and position in the event SimTrack vector
430  for( unsigned it=0; it<nTks; ++it ) {
431  geantToIndex_[ simTracks[it].trackId() ] = it;
432  }
433 }
std::map< unsigned, unsigned > geantToIndex_
std::vector< EcalSimPhotonMCTruth > ContainmentCorrectionAnalyzer::findMcTruth ( std::vector< SimTrack > &  theSimTracks,
std::vector< SimVertex > &  theSimVertices 
)
private

store this electron since it's from a converted photon

Definition at line 292 of file ContainmentCorrectionAnalyzer.cc.

References funct::abs(), reco::if(), SimVertex::parentIndex(), CoreSimVertex::position(), and query::result.

292  {
293 
294  std::vector<EcalSimPhotonMCTruth> result;
295 
296  geantToIndex_.clear();
297  // int idTrk1_[10]; // UNUSED
298  // int idTrk2_[10]; // UNUSED
299 
300  // Local variables
301  // const int SINGLE=1; // UNUSED
302  // const int DOUBLE=2; // UNUSED
303  // const int PYTHIA=3; // UNUSED
304  const int ELECTRON_FLAV=1;
305  const int PIZERO_FLAV=2;
306  const int PHOTON_FLAV=3;
307 
308  // int ievtype=0; // UNUSED
309  int ievflav=0;
310  std::vector<SimTrack*> photonTracks;
311  std::vector<SimTrack*> pizeroTracks;
312  std::vector<const SimTrack *> trkFromConversion;
313  SimVertex primVtx;
314  std::vector<int> convInd;
315 
316  fillMcTruth(theSimTracks, theSimVertices);
317  int iPV=-1;
318  int partType1=0;
319  int partType2=0;
320  std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin();
321  if ( !(*iFirstSimTk).noVertex() ) {
322  iPV = (*iFirstSimTk).vertIndex();
323  int vtxId = (*iFirstSimTk).vertIndex();
324  primVtx = theSimVertices[vtxId];
325  partType1 = (*iFirstSimTk).type();
326  }
327 
328  // Look at a second track
329  iFirstSimTk++;
330  if ( iFirstSimTk!= theSimTracks.end() ) {
331  if ( (*iFirstSimTk).vertIndex() == iPV) {
332  partType2 = (*iFirstSimTk).type();
333  }
334  }
335  int npv=0;
336  int iPho=0;
337  for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
338  if ( (*iSimTk).noVertex() ) continue;
339  // int vertexId = (*iSimTk).vertIndex(); // UNUSED
340  // SimVertex vertex = theSimVertices[vertexId]; // UNUSED
341  if ( (*iSimTk).vertIndex() == iPV ) {
342  npv++;
343  if ( (*iSimTk).type() == 22) {
344  convInd.push_back(0);
345  photonTracks.push_back( &(*iSimTk) );
346  // math::XYZTLorentzVectorD momentum = (*iSimTk).momentum(); // UNUSED
347  }
348  }
349  }
350 
351  if(npv > 4) { // ievtype = PYTHIA; // UNUSED
352  } else if(npv == 1) {
353  if( abs(partType1) == 11 ) { /* ievtype = SINGLE; ==UNUSED== */ ievflav = ELECTRON_FLAV; }
354  else if(partType1 == 111) { /* ievtype = SINGLE; ==UNUSED== */ ievflav = PIZERO_FLAV; }
355  else if(partType1 == 22) { /* ievtype = SINGLE; ==UNUSED== */ ievflav = PHOTON_FLAV; }
356  } else if(npv == 2) {
357  if ( abs(partType1) == 11 && abs(partType2) == 11 ) { /* ievtype = DOUBLE; ==UNUSED== */ ievflav = ELECTRON_FLAV; }
358  else if(partType1 == 111 && partType2 == 111) { /* ievtype = DOUBLE; ==UNUSED== */ ievflav = PIZERO_FLAV; }
359  else if(partType1 == 22 && partType2 == 22) { /* ievtype = DOUBLE; ==UNUSED== */ ievflav = PHOTON_FLAV; }
360  }
361 
362  // Look into converted photons
363  int isAconversion=0;
364  if(ievflav == PHOTON_FLAV) {
365 
366  int nConv=0;
367  int iConv=0;
368  iPho=0;
369  for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
370  trkFromConversion.clear();
371  for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
372  if ( (*iSimTk).noVertex() ) continue;
373  if ( (*iSimTk).vertIndex() == iPV ) continue;
374  if ( abs((*iSimTk).type()) != 11 ) continue;
375  int vertexId = (*iSimTk).vertIndex();
376  SimVertex vertex = theSimVertices[vertexId];
377  int motherId=-1;
378  if ( vertex.parentIndex() ) {
379  unsigned motherGeantId = vertex.parentIndex();
380  std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
381  if(association != geantToIndex_.end() )
382  motherId = association->second;
383  //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
384 
385  if ( theSimTracks[motherId].trackId() == (*iPhoTk)->trackId() ) {
387  trkFromConversion.push_back(&(*iSimTk ) );
388  }
389  }
390  } // loop over the SimTracks
391 
392  if ( trkFromConversion.size() > 0 ) {
393  isAconversion=1;
394  nConv++;
395  convInd[iPho]=nConv;
396  int convVtxId = trkFromConversion[0]->vertIndex();
397  SimVertex convVtx = theSimVertices[convVtxId];
398  math::XYZTLorentzVectorD vtxPosition = convVtx.position();
399  // math::XYZTLorentzVectorD momentum = (*iPhoTk)->momentum(); // UNUSED
400  if ( nConv <= 10) {
401  if ( trkFromConversion.size() > 1) {
402  // idTrk1_[iConv]= trkFromConversion[0]->trackId(); // UNUSED
403  // idTrk2_[iConv]= trkFromConversion[1]->trackId(); // UNUSED
404  } else {
405  // idTrk1_[iConv]=trkFromConversion[0]->trackId(); // UNUSED
406  // idTrk2_[iConv]=-1; // UNUSED
407  }
408  }
409  iConv++;
410 
411  result.push_back( EcalSimPhotonMCTruth(isAconversion, (*iPhoTk)->momentum(), vtxPosition.pt(), vtxPosition.z() , vtxPosition, primVtx.position(), trkFromConversion ));
412  } else {
413  isAconversion=0;
414  math::XYZTLorentzVectorD vtxPosition(0.,0.,0.,0.);
415  result.push_back( EcalSimPhotonMCTruth(isAconversion, (*iPhoTk)->momentum(), vtxPosition.pt(), vtxPosition.z() , vtxPosition, primVtx.position(), trkFromConversion ));
416  }
417  iPho++;
418  } // loop over the primary photons
419  } // Event with one or two photons
420 
421  return result;
422 }
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
std::map< unsigned, unsigned > geantToIndex_
int parentIndex() const
Definition: SimVertex.h:33
tuple result
Definition: query.py:137
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const math::XYZTLorentzVectorD & position() const
Definition: CoreSimVertex.h:26
if(dp >Float(M_PI)) dp-
void fillMcTruth(std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)

Member Data Documentation

edm::InputTag ContainmentCorrectionAnalyzer::BarrelSuperClusterCollection_
private

Definition at line 64 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::e1
private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::e25
private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::e9
private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

edm::InputTag ContainmentCorrectionAnalyzer::EndcapSuperClusterCollection_
private

Definition at line 65 of file ContainmentCorrectionAnalyzer.h.

std::map<unsigned, unsigned> ContainmentCorrectionAnalyzer::geantToIndex_
private

Definition at line 79 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_converted
private

Definition at line 89 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_e25EtrueReference
private

Definition at line 83 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_e9EtrueReference
private

Definition at line 82 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_eRecoEtrueReference
private

Definition at line 81 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_eTrue
private

Definition at line 87 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_converted
private

Definition at line 90 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_e25EtrueReference
private

Definition at line 86 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_e9EtrueReference
private

Definition at line 85 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_eRecoEtrueReference
private

Definition at line 84 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_eTrue
private

Definition at line 88 of file ContainmentCorrectionAnalyzer.h.

std::vector<int> ContainmentCorrectionAnalyzer::iMC
private

Definition at line 77 of file ContainmentCorrectionAnalyzer.h.

std::vector<int> ContainmentCorrectionAnalyzer::isConverted
private

Definition at line 71 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcEnergy
private

Definition at line 70 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcEta
private

Definition at line 70 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcPhi
private

Definition at line 70 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcPt
private

Definition at line 70 of file ContainmentCorrectionAnalyzer.h.

int ContainmentCorrectionAnalyzer::nMCphotons
private

Definition at line 69 of file ContainmentCorrectionAnalyzer.h.

int ContainmentCorrectionAnalyzer::nRECOphotons
private

Definition at line 73 of file ContainmentCorrectionAnalyzer.h.

edm::InputTag ContainmentCorrectionAnalyzer::reducedBarrelRecHitCollection_
private

Definition at line 66 of file ContainmentCorrectionAnalyzer.h.

edm::InputTag ContainmentCorrectionAnalyzer::reducedEndcapRecHitCollection_
private

Definition at line 67 of file ContainmentCorrectionAnalyzer.h.

std::vector<int> ContainmentCorrectionAnalyzer::seedXtal
private

Definition at line 76 of file ContainmentCorrectionAnalyzer.h.

edm::InputTag ContainmentCorrectionAnalyzer::SimTrackCollection_
private

Definition at line 62 of file ContainmentCorrectionAnalyzer.h.

edm::InputTag ContainmentCorrectionAnalyzer::SimVertexCollection_
private

Definition at line 63 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterEnergy
private

Definition at line 74 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterEt
private

Definition at line 74 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterEta
private

Definition at line 74 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterPhi
private

Definition at line 74 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::x_vtx
private

Definition at line 72 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::y_vtx
private

Definition at line 72 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::z_vtx
private

Definition at line 72 of file ContainmentCorrectionAnalyzer.h.