TrackerFilters/plugins/CosmicTIFTrigFilter.cc

0001 // livio.fano@cern.ch
0002
0003 #include "CosmicTIFTrigFilter.h"
0004
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "HepMC/GenVertex.h"
0009 #include <map>
0010 #include <vector>
0011
0012 using namespace std;
0013 namespace cms
0014
0015 {
0016   CosmicTIFTrigFilter::CosmicTIFTrigFilter(const edm::ParameterSet &conf)
0017       : m_Token(consumes<edm::HepMCProduct>(
0018             conf.getParameter<edm::ParameterSet>("Generator").getParameter<std::string>("HepMCProductLabel"))) {
0019     trigconf = conf.getParameter<int>("trig_conf");
0020     trigS1 = conf.getParameter<std::vector<double>>("PosScint1");
0021     trigS2 = conf.getParameter<std::vector<double>>("PosScint2");
0022     trigS3 = conf.getParameter<std::vector<double>>("PosScint3");
0023     trigS4 = conf.getParameter<std::vector<double>>("PosScint4");
0024
0025     /*
0026   std::cout << "S1 = " << trigS1[0] << ", " <<  trigS1[1] << ", " <<trigS1[2]
0027             << "\nS2 = " << trigS2[0] << ", " <<  trigS2[1] << ", " <<trigS2[2]
0028             << "\nS3 = " << trigS3[0] << ", " <<  trigS3[1] << ", " <<trigS3[2]
0029             << "\nS4 = " << trigS4[0] << ", " <<  trigS4[1] << ", " <<trigS4[2]
0030             << std::endl;
0031   */
0032   }
0033
0034   bool CosmicTIFTrigFilter::filter(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0035     edm::Handle<edm::HepMCProduct> HepMCEvt;
0036     iEvent.getByToken(m_Token, HepMCEvt);
0037
0038     const HepMC::GenEvent *MCEvt = HepMCEvt->GetEvent();
0039
0040     bool hit1 = false;
0041     bool hit2 = false;
0042     bool hit3 = false;
0043     bool hit4 = false;
0044
0045     for (HepMC::GenEvent::particle_const_iterator i = MCEvt->particles_begin(); i != MCEvt->particles_end(); ++i) {
0046       int myId = (*i)->pdg_id();
0047       if (abs(myId) == 13) {
0048         // Get the muon position and momentum
0049         HepMC::GenVertex *pv = (*i)->production_vertex();
0050         const HepMC::FourVector &vertex = pv->position();
0051
0052         HepMC::FourVector momentum = (*i)->momentum();
0053
0054         // std::cout << "\t vertex for cut = " << vertex << std::endl;
0055         // std::cout << "\t momentum  = " << momentum << std::endl;
0056
0057         if (trigconf == 1) {
0058           HepMC::FourVector S1(trigS1[0], trigS1[1], trigS1[2], 0.);
0059           HepMC::FourVector S2(trigS2[0], trigS2[1], trigS2[2], 0.);
0060           HepMC::FourVector S3(trigS3[0], trigS3[1], trigS3[2], 0.);
0061
0062           hit1 = Sci_trig(vertex, momentum, S1);
0063           hit2 = Sci_trig(vertex, momentum, S2);
0064           hit3 = Sci_trig(vertex, momentum, S3);
0065
0066           // trigger conditions
0067
0068           if ((hit1 && hit2) || (hit3 && hit2)) {
0069             /*
0070           cout << "\tGot a trigger in configuration A " << endl;
0071           if(hit1)cout << "hit1 " << endl;
0072           if(hit2)cout << "hit2 " << endl;
0073           if(hit3)cout << "hit3 " << endl;
0074           */
0075             trig1++;
0076             return true;
0077           }
0078         } else if (trigconf == 2) {
0079           HepMC::FourVector S1(trigS1[0], trigS1[1], trigS1[2], 0.);
0080           HepMC::FourVector S2(trigS2[0], trigS2[1], trigS2[2], 0.);
0081           HepMC::FourVector S3(trigS3[0], trigS3[1], trigS3[2], 0.);
0082
0083           hit1 = Sci_trig(vertex, momentum, S1);
0084           hit2 = Sci_trig(vertex, momentum, S2);
0085           hit3 = Sci_trig(vertex, momentum, S3);
0086
0087           // trigger conditions
0088
0089           if ((hit1 && hit2) || (hit3 && hit2)) {
0090             /*
0091           cout << "\tGot a trigger in configuration B " << endl;
0092           if(hit1)cout << "hit1 " << endl;
0093           if(hit2)cout << "hit2 " << endl;
0094           if(hit3)cout << "hit3 " << endl;
0095           */
0096             trig2++;
0097             return true;
0098           }
0099
0100         } else if (trigconf == 3) {
0101           HepMC::FourVector S1(trigS1[0], trigS1[1], trigS1[2], 0.);
0102           HepMC::FourVector S2(trigS2[0], trigS2[1], trigS2[2], 0.);
0103           HepMC::FourVector S3(trigS3[0], trigS3[1], trigS3[2], 0.);
0104           HepMC::FourVector S4(trigS4[0], trigS4[1], trigS4[2], 0.);
0105
0106           /*          std::cout << "S1 = " << S1.x() << "," << S1.y() << ", " <<
0107            S1.z()
0108                   << "\nS2 = " << S2.x() << "," << S2.y() << ", " << S2.z()
0109                   << "\nS3 = " << S3.x() << "," << S3.y() << ", " << S3.z()
0110                   << "\nS4 = " << S4.x() << "," << S4.y() << ", " << S4.z()
0111                   << std::endl;
0112         */
0113
0114           hit1 = Sci_trig(vertex, momentum, S1);
0115           hit2 = Sci_trig(vertex, momentum, S2);
0116           hit3 = Sci_trig(vertex, momentum, S3);
0117           hit4 = Sci_trig(vertex, momentum, S4);
0118
0119           // trigger conditions
0120           if ((hit1 && hit2) || (hit3 && hit2) || (hit1 && hit4) || (hit3 && hit4)) {
0121             /*
0122           cout << "\tGot a trigger in configuration C " << endl;
0123           if(hit1)cout << "hit1 " << endl;
0124           if(hit2)cout << "hit2 " << endl;
0125           if(hit3)cout << "hit3 " << endl;
0126           if(hit4)cout << "hit4 " << endl;
0127           */
0128
0129             trig3++;
0130             return true;
0131           }
0132         }
0133       }
0134     }
0135
0136     return false;
0137   }
0138
0139   bool CosmicTIFTrigFilter::Sci_trig(const HepMC::FourVector &vertex,
0140                                      const HepMC::FourVector &momentum,
0141                                      const HepMC::FourVector &S) {
0142     float x0 = vertex.x();
0143     float y0 = vertex.y();
0144     float z0 = vertex.z();
0145     float px0 = momentum.x();
0146     float py0 = momentum.y();
0147     float pz0 = momentum.z();
0148     float Sx = S.x();
0149     float Sy = S.y();
0150     float Sz = S.z();
0151
0152     float zs = (Sy - y0) * (pz0 / py0) + z0;
0153     float xs = (Sy - y0) * (px0 / py0) + x0;
0154
0155     //    std::cout << Sx << " " << Sz << " " << xs << " " << zs << std::endl;
0156     // std::cout << x0 << " " << z0 << " " << px0 << " " << py0 << " " << pz0 <<
0157     // endl;
0158
0159     if ((xs < Sx + 500 && xs > Sx - 500) && (zs < Sz + 500 && zs > Sz - 500)) {
0160       // std::cout << "PASSED" << std::endl;
0161       return true;
0162     } else {
0163       return false;
0164     }
0165   }
0166 }  // namespace cms