CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
KineParticleFilter.cc
Go to the documentation of this file.
2 
3 //Framework Headers
5 
6 #include <vector>
7 #include <iterator>
8 
11 {
12 
13  // Set the kinematic cuts
14  // Upper abs(eta) bound
15  etaMax = kine.getParameter<double>("etaMax");
16  // Lower pT bound (charged, in GeV/c)
17  pTMin = kine.getParameter<double>("pTMin");
18  // Lower E bound - reject (all, in GeV)
19  EMin = kine.getParameter<double>("EMin");
20  // Lower E bound - accept (all, in GeV)
21  EMax = kine.getParameter<double>("EProton");
22 
23  // pdg codes of the particles to be removed from the events
24  // ParameterSet cannot handle sets, only vectors
25  std::vector<int> tmpcodes
26  = kine.getUntrackedParameter< std::vector<int> >
27  ("forbiddenPdgCodes", std::vector<int>() );
28 
29  std::copy(tmpcodes.begin(),
30  tmpcodes.end(),
31  std::insert_iterator< std::set<int> >(forbiddenPdgCodes,
32  forbiddenPdgCodes.begin() ));
33 
34  if( !forbiddenPdgCodes.empty() ) {
35  std::cout<<"KineParticleFilter : Forbidden PDG codes : ";
36  copy(forbiddenPdgCodes.begin(), forbiddenPdgCodes.end(),
37  std::ostream_iterator<int>(std::cout, " "));
38  }
39 
40  // Change eta cuts to cos**2(theta) cuts (less CPU consuming)
41  if ( etaMax > 20. ) etaMax = 20.; // Protection against paranoid people.
42  double cosMax = (std::exp(2.*etaMax)-1.) / (std::exp(2.*etaMax)+1.);
43  cos2Max = cosMax*cosMax;
44 
45  double etaPreshMin = 1.479;
46  double etaPreshMax = 1.594;
47  double cosPreshMin = (std::exp(2.*etaPreshMin)-1.) / (std::exp(2.*etaPreshMin)+1.);
48  double cosPreshMax = (std::exp(2.*etaPreshMax)-1.) / (std::exp(2.*etaPreshMax)+1.);
49  cos2PreshMin = cosPreshMin*cosPreshMin;
50  cos2PreshMax = cosPreshMax*cosPreshMax;
51 
52  // Change pt cut to pt**2 cut (less CPU consuming)
53  pTMin *= pTMin;
54 
55 }
56 
58 {
59 
60  // Do not consider quarks, gluons, Z, W, strings, diquarks
61  // ... and supesymmetric particles
62  int pId = abs(p->pid());
63 
64  // Vertices are coming with pId = 0
65  if ( pId != 0 ) {
66  bool particleCut = ( pId > 10 && pId != 12 && pId != 14 &&
67  pId != 16 && pId != 18 && pId != 21 &&
68  (pId < 23 || pId > 40 ) &&
69  (pId < 81 || pId > 100 ) && pId != 2101 &&
70  pId != 3101 && pId != 3201 && pId != 1103 &&
71  pId != 2103 && pId != 2203 && pId != 3103 &&
72  pId != 3203 && pId != 3303 );
73  // particleCut = particleCut || pId == 0;
74 
75 
76  if ( !particleCut ) return false;
77 
78  // Keep protons with energy in excess of 5 TeV
79  bool protonTaggers = (pId == 2212 && p->E() >= EMax) ;
80  if ( protonTaggers ) return true;
81 
82  std::set<int>::iterator is = forbiddenPdgCodes.find(pId);
83  if( is != forbiddenPdgCodes.end() ) return false;
84 
85  // bool kineCut = pId == 0;
86  // Cut on kinematic properties
87  // Cut on the energy of all particles
88  bool eneCut = p->E() >= EMin;
89  if (!eneCut) return false;
90 
91  // Cut on the transverse momentum of charged particles
92  bool pTCut = p->charge()==0 || p->Perp2()>=pTMin;
93  if (!pTCut) return false;
94 
95  // Cut on eta if the origin vertex is close to the beam
96  // bool etaCut = (p->vertex()-mainVertex).perp()>5. || fabs(p->eta())<=etaMax;
97  bool etaCut = (p->vertex()-mainVertex).Perp2()>25. || p->cos2Theta()<= cos2Max;
98 
99  /*
100  if ( etaCut != etaCut2 )
101  cout << "WANRNING ! etaCut != etaCut2 "
102  << etaCut << " "
103  << etaCut2 << " "
104  << (p->eta()) << " " << etaMax << " "
105  << p->vect().cos2Theta() << " " << cos2Max << endl;
106  */
107  if (!etaCut) return false;
108 
109  // Cut on the origin vertex position (prior to the ECAL for all
110  // particles, except for muons ! Just modified: Muons included as well !
111  double radius2 = p->R2();
112  double zed = fabs(p->Z());
113  double cos2Tet = p->cos2ThetaV();
114  // Ecal entrance
115  bool ecalAcc = ( (radius2<129.01*129.01 && zed<317.01) ||
116  (cos2Tet>cos2PreshMin && cos2Tet<cos2PreshMax
117  && radius2<171.11*171.11 && zed<317.01) );
118 
119  return ecalAcc;
120 
121  } else {
122  // Cut for vertices
123  double radius2 = p->Perp2();
124  double zed = fabs(p->Pz());
125  double cos2Tet = p->cos2Theta();
126 
127  // Vertices must be before the Ecal entrance
128  bool ecalAcc = ( (radius2<129.01*129.01 && zed<317.01) ||
129  (cos2Tet>cos2PreshMin && cos2Tet<cos2PreshMax
130  && radius2<171.11*171.11 && zed<317.01) );
131 
132  return ecalAcc;
133 
134  }
135 
136 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
#define abs(x)
Definition: mlp_lapack.h:159
int pid() const
get the HEP particle ID number
Definition: RawParticle.h:264
double R2() const
vertex radius**2
Definition: RawParticle.h:278
double cos2ThetaV() const
Definition: RawParticle.h:268
double Z() const
z of vertex
Definition: RawParticle.h:275
double charge() const
get the MEASURED charge
Definition: RawParticle.h:281
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:284
std::set< int > forbiddenPdgCodes
KineParticleFilter(const edm::ParameterSet &kine)
virtual bool isOKForMe(const RawParticle *p) const
the real selection is done here
tuple cout
Definition: gather_cfg.py:121
XYZTLorentzVector mainVertex
double cos2Theta() const
Cos**2(theta) is faster to determine than eta.
Definition: RawParticle.h:267