Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:40 

0001 
0002 //-ap #include "Configuration/CSA06Skimming/interface/MCParticlePairFilter.h"
0003 
0004 #include "GeneratorInterface/GenFilters/plugins/MCParticlePairFilter.h"
0005 #include "GeneratorInterface/GenFilters/plugins/MCFilterZboostHelper.h"
0006 
0007 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0008 #include <iostream>
0009 
0010 using namespace edm;
0011 using namespace std;
0012 
0013 MCParticlePairFilter::MCParticlePairFilter(const edm::ParameterSet& iConfig)
0014     : token_(consumes<edm::HepMCProduct>(
0015           edm::InputTag(iConfig.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0016       particleCharge(iConfig.getUntrackedParameter("ParticleCharge", 0)),
0017       minInvMass(iConfig.getUntrackedParameter("MinInvMass", 0.)),
0018       maxInvMass(iConfig.getUntrackedParameter("MaxInvMass", 14000.)),
0019       minDeltaPhi(iConfig.getUntrackedParameter("MinDeltaPhi", 0.)),
0020       maxDeltaPhi(iConfig.getUntrackedParameter("MaxDeltaPhi", 6.3)),
0021       minDeltaR(iConfig.getUntrackedParameter("MinDeltaR", 0.)),
0022       maxDeltaR(iConfig.getUntrackedParameter("MaxDeltaR", 10000.)),
0023       betaBoost(iConfig.getUntrackedParameter("BetaBoost", 0.)) {
0024   //here do whatever other initialization is needed
0025   vector<int> defpid1;
0026   defpid1.push_back(0);
0027   particleID1 = iConfig.getUntrackedParameter<vector<int> >("ParticleID1", defpid1);
0028   vector<int> defpid2;
0029   defpid2.push_back(0);
0030   particleID2 = iConfig.getUntrackedParameter<vector<int> >("ParticleID2", defpid2);
0031   vector<double> defptmin;
0032   defptmin.push_back(0.);
0033   ptMin = iConfig.getUntrackedParameter<vector<double> >("MinPt", defptmin);
0034   vector<double> defpmin;
0035   defpmin.push_back(0.);
0036   pMin = iConfig.getUntrackedParameter<vector<double> >("MinP", defpmin);
0037 
0038   vector<double> defetamin;
0039   defetamin.push_back(-10.);
0040   etaMin = iConfig.getUntrackedParameter<vector<double> >("MinEta", defetamin);
0041   vector<double> defetamax;
0042   defetamax.push_back(10.);
0043   etaMax = iConfig.getUntrackedParameter<vector<double> >("MaxEta", defetamax);
0044   vector<int> defstat;
0045   defstat.push_back(0);
0046   status = iConfig.getUntrackedParameter<vector<int> >("Status", defstat);
0047 
0048   // check for correct size
0049   if (ptMin.size() != 2 || pMin.size() != 2 || etaMin.size() != 2 || etaMax.size() != 2 || status.size() != 2) {
0050     cout << "WARNING: MCParticlePairFilter : size of some vectors not matching with 2!!" << endl;
0051   }
0052 
0053   // if ptMin size smaller than 2, fill up further with defaults
0054   if (2 > ptMin.size()) {
0055     vector<double> defptmin2;
0056     for (unsigned int i = 0; i < 2; i++) {
0057       defptmin2.push_back(0.);
0058     }
0059     ptMin = defptmin2;
0060   }
0061   // if pMin size smaller than 2, fill up further with defaults
0062   if (2 > pMin.size()) {
0063     vector<double> defpmin2;
0064     for (unsigned int i = 0; i < 2; i++) {
0065       defpmin2.push_back(0.);
0066     }
0067     pMin = defpmin2;
0068   }
0069   // if etaMin size smaller than 2, fill up further with defaults
0070   if (2 > etaMin.size()) {
0071     vector<double> defetamin2;
0072     for (unsigned int i = 0; i < 2; i++) {
0073       defetamin2.push_back(-10.);
0074     }
0075     etaMin = defetamin2;
0076   }
0077   // if etaMax size smaller than 2, fill up further with defaults
0078   if (2 > etaMax.size()) {
0079     vector<double> defetamax2;
0080     for (unsigned int i = 0; i < 2; i++) {
0081       defetamax2.push_back(10.);
0082     }
0083     etaMax = defetamax2;
0084   }
0085   // if status size smaller than 2, fill up further with defaults
0086   if (2 > status.size()) {
0087     vector<int> defstat2;
0088     for (unsigned int i = 0; i < 2; i++) {
0089       defstat2.push_back(0);
0090     }
0091     status = defstat2;
0092   }
0093 }
0094 
0095 MCParticlePairFilter::~MCParticlePairFilter() {
0096   // do anything here that needs to be done at desctruction time
0097   // (e.g. close files, deallocate resources etc.)
0098 }
0099 
0100 // ------------ method called to skim the data  ------------
0101 bool MCParticlePairFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0102   using namespace edm;
0103   bool accepted = false;
0104   Handle<HepMCProduct> evt;
0105   iEvent.getByToken(token_, evt);
0106   const double pi = 3.14159;
0107 
0108   vector<HepMC::GenParticle*> typeApassed;
0109   vector<HepMC::GenParticle*> typeBpassed;
0110 
0111   const HepMC::GenEvent* myGenEvent = evt->GetEvent();
0112 
0113   for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0114        ++p) {
0115     // check for type A conditions
0116     bool gottypeAID = false;
0117     for (unsigned int j = 0; j < particleID1.size(); ++j) {
0118       if (abs((*p)->pdg_id()) == abs(particleID1[j]) || particleID1[j] == 0) {
0119         gottypeAID = true;
0120         break;
0121       }
0122     }
0123     if (gottypeAID) {
0124       HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0125       if (mom.perp() > ptMin[0] && mom.rho() > pMin[0] && mom.eta() > etaMin[0] && mom.eta() < etaMax[0] &&
0126           ((*p)->status() == status[0] || status[0] == 0)) {
0127         // passed A type conditions ...
0128         // ... now check pair-conditions with B type passed particles
0129         unsigned int i = 0;
0130         double deltaphi;
0131         double phi1 = mom.phi();
0132         double phi2;
0133         double deltaeta;
0134         double eta1 = mom.eta();
0135         double eta2;
0136         double deltaR;
0137 
0138         double tot_x = 0.;
0139         double tot_y = 0.;
0140         double tot_z = 0.;
0141         double tot_e = 0.;
0142         double invmass = 0.;
0143         int charge1 = 0;
0144         int combcharge = 0;
0145         while (!accepted && i < typeBpassed.size()) {
0146           tot_x = mom.px();
0147           tot_y = mom.py();
0148           tot_z = mom.pz();
0149           tot_e = mom.e();
0150           charge1 = charge((*p)->pdg_id());
0151           HepMC::FourVector mom_i = MCFilterZboostHelper::zboost(typeBpassed[i]->momentum(), betaBoost);
0152           tot_x += mom_i.px();
0153           tot_y += mom_i.py();
0154           tot_z += mom_i.pz();
0155           tot_e += mom_i.e();
0156           invmass = sqrt(tot_e * tot_e - tot_x * tot_x - tot_y * tot_y - tot_z * tot_z);
0157           combcharge = charge1 * charge(typeBpassed[i]->pdg_id());
0158           if (invmass > minInvMass && invmass < maxInvMass) {
0159             phi2 = mom_i.phi();
0160             deltaphi = fabs(phi1 - phi2);
0161             if (deltaphi > pi)
0162               deltaphi = 2. * pi - deltaphi;
0163             if (deltaphi > minDeltaPhi && deltaphi < maxDeltaPhi) {
0164               eta2 = mom_i.eta();
0165               deltaeta = fabs(eta1 - eta2);
0166               deltaR = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0167               if (deltaR > minDeltaR && deltaR < maxDeltaR) {
0168                 if (combcharge * particleCharge >= 0) {
0169                   accepted = true;
0170                 }
0171               }
0172             }
0173           }
0174           i++;
0175         }
0176         // if we found a matching pair quit the loop
0177         if (accepted)
0178           break;
0179         else {
0180           typeApassed.push_back(*p);  // else remember the particle to have passed type A conditions
0181         }
0182       }
0183     }
0184 
0185     // check for type B conditions
0186 
0187     bool gottypeBID = false;
0188     for (unsigned int j = 0; j < particleID2.size(); ++j) {
0189       if (abs((*p)->pdg_id()) == abs(particleID2[j]) || particleID2[j] == 0) {
0190         gottypeBID = true;
0191         break;
0192       }
0193     }
0194     if (gottypeBID) {
0195       HepMC::FourVector mom = MCFilterZboostHelper::zboost((*p)->momentum(), betaBoost);
0196       if (mom.perp() > ptMin[1] && mom.rho() > pMin[1] && mom.eta() > etaMin[1] && mom.eta() < etaMax[1] &&
0197           ((*p)->status() == status[1] || status[1] == 0)) {
0198         // passed B type conditions ...
0199         // ... now check pair-conditions with A type passed particles vector
0200         unsigned int i = 0;
0201         double deltaphi;
0202         double phi1 = mom.phi();
0203         double phi2;
0204         double deltaeta;
0205         double eta1 = mom.eta();
0206         double eta2;
0207         double deltaR;
0208 
0209         double tot_x = 0.;
0210         double tot_y = 0.;
0211         double tot_z = 0.;
0212         double tot_e = 0.;
0213         double invmass = 0.;
0214         int charge1 = 0;
0215         int combcharge = 0;
0216         while (!accepted && i < typeApassed.size()) {
0217           if ((*p) != typeApassed[i]) {
0218             tot_x = mom.px();
0219             tot_y = mom.py();
0220             tot_z = mom.pz();
0221             tot_e = mom.e();
0222             charge1 = charge((*p)->pdg_id());
0223             HepMC::FourVector mom_i = MCFilterZboostHelper::zboost(typeApassed[i]->momentum(), betaBoost);
0224             tot_x += mom_i.px();
0225             tot_y += mom_i.py();
0226             tot_z += mom_i.pz();
0227             tot_e += mom_i.e();
0228             invmass = sqrt(tot_e * tot_e - tot_x * tot_x - tot_y * tot_y - tot_z * tot_z);
0229             combcharge = charge1 * charge(typeApassed[i]->pdg_id());
0230             if (invmass > minInvMass && invmass < maxInvMass) {
0231               phi2 = mom_i.phi();
0232               deltaphi = fabs(phi1 - phi2);
0233               if (deltaphi > pi)
0234                 deltaphi = 2. * pi - deltaphi;
0235               if (deltaphi > minDeltaPhi && deltaphi < maxDeltaPhi) {
0236                 eta2 = mom_i.eta();
0237                 deltaeta = fabs(eta1 - eta2);
0238                 deltaR = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0239                 if (deltaR > minDeltaR && deltaR < maxDeltaR) {
0240                   if (combcharge * particleCharge >= 0) {
0241                     accepted = true;
0242                   }
0243                 }
0244               }
0245             }
0246           }
0247           i++;
0248         }
0249         // if we found a matching pair quit the loop
0250         if (accepted)
0251           break;
0252         else {
0253           typeBpassed.push_back(*p);  // else remember the particle to have passed type B conditions
0254         }
0255       }
0256     }
0257   }
0258 
0259   if (accepted) {
0260     return true;
0261   } else {
0262     return false;
0263   }
0264 }
0265 
0266 int MCParticlePairFilter::charge(int Id) const {
0267   //...Purpose: to give three times the charge for a particle/parton.
0268 
0269   //      ID     = particle ID
0270   //      hepchg = particle charge times 3
0271 
0272   int kqa, kq1, kq2, kq3, kqj, irt, kqx, kqn;
0273   int hepchg;
0274 
0275   constexpr const int ichg[109] = {-1, 2, -1, 2, -1, 2, -1, 2, 0,  0, -3, 0, -3, 0, -3, 0, -3, 0, 0,  0, 0,  0,
0276                                    0,  3, 0,  0, 0,  0, 0,  0, 3,  0, 3,  6, 0,  0, 3,  6, 0,  0, -1, 2, -1, 2,
0277                                    -1, 2, 0,  0, 0,  0, -3, 0, -3, 0, -3, 0, 0,  0, 0,  0, -1, 2, -1, 2, -1, 2,
0278                                    0,  0, 0,  0, -3, 0, -3, 0, -3, 0, 3,  3, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0,
0279                                    0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0,  0, 0};
0280 
0281   //...Initial values. Simple case of direct readout.
0282   hepchg = 0;
0283   kqa = abs(Id);
0284   kqn = kqa / 1000000000 % 10;
0285   kqx = kqa / 1000000 % 10;
0286   kq3 = kqa / 1000 % 10;
0287   kq2 = kqa / 100 % 10;
0288   kq1 = kqa / 10 % 10;
0289   kqj = kqa % 10;
0290   irt = kqa % 10000;
0291 
0292   //...illegal or ion
0293   //...set ion charge to zero - not enough information
0294   if (kqa == 0 || kqa >= 10000000) {
0295     if (kqn == 1) {
0296       hepchg = 0;
0297     }
0298   }
0299   //... direct translation
0300   else if (kqa <= 100) {
0301     hepchg = ichg[kqa - 1];
0302   }
0303   //... KS and KL (and undefined)
0304   else if (kqj == 0) {
0305     hepchg = 0;
0306   }
0307   //C... direct translation
0308   else if (kqx > 0 && irt < 100) {
0309     hepchg = ichg[irt - 1];
0310     if (kqa == 1000017 || kqa == 1000018) {
0311       hepchg = 0;
0312     }
0313     if (kqa == 1000034 || kqa == 1000052) {
0314       hepchg = 0;
0315     }
0316     if (kqa == 1000053 || kqa == 1000054) {
0317       hepchg = 0;
0318     }
0319     if (kqa == 5100061 || kqa == 5100062) {
0320       hepchg = 6;
0321     }
0322   }
0323   //...Construction from quark content for heavy meson, diquark, baryon.
0324   //...Mesons.
0325   else if (kq3 == 0) {
0326     hepchg = ichg[kq2 - 1] - ichg[kq1 - 1];
0327     //...Strange or beauty mesons.
0328     if ((kq2 == 3) || (kq2 == 5)) {
0329       hepchg = ichg[kq1 - 1] - ichg[kq2 - 1];
0330     }
0331   } else if (kq1 == 0) {
0332     //...Diquarks.
0333     hepchg = ichg[kq3 - 1] + ichg[kq2 - 1];
0334   }
0335 
0336   else {
0337     //...Baryons
0338     hepchg = ichg[kq3 - 1] + ichg[kq2 - 1] + ichg[kq1 - 1];
0339   }
0340 
0341   //... fix sign of charge
0342   if (Id < 0 && hepchg != 0) {
0343     hepchg = -1 * hepchg;
0344   }
0345 
0346   // cout << hepchg<< endl;
0347   return hepchg;
0348 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!