/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/GeneratorInterface/GenFilters/src/BHFilter.cc

Go to the documentation of this file.

// livio.fano@cern.ch

#include "GeneratorInterface/GenFilters/interface/BHFilter.h"
#include "DataFormats/Common/interface/Handle.h"

using namespace std;
namespace cms

{
BHFilter::BHFilter(const edm::ParameterSet& conf):    conf_(conf)
{
  rBounds = conf_.getParameter< vector<double> >("radii");
  zBounds = conf_.getParameter< vector<double> >("zeds");
  bFields = conf_.getParameter< vector<double> >("bfiel");
  bReduction = conf_.getParameter< double >("factor");
  trig_ = conf_.getParameter< int >("trig_type");
  trig2_ = conf_.getParameter< int >("scintillators_type");

  for ( unsigned i=0; i<bFields.size(); ++i ) { 
    bFields[i] *= bReduction;
    cout << "r/z/b = " << rBounds[i] << " " << zBounds[i] << " " << bFields[i] << endl;
  }

  if(trig_==0){ cout <<endl << "trigger is both + and - BSC " << endl;}
  if(trig_==1){ cout <<endl << "trigger is + side BSC " << endl;}
  if(trig_==-1){ cout <<endl << "trigger is - side BSC " << endl;}
  
  cout << endl;
  
  if(trig2_==0){ cout <<endl << "trigger is both PADs and DISKs " << endl;}
  if(trig2_==-1){ cout <<endl << "trigger is only PADs " << endl;}
  if(trig2_==1){ cout <<endl << "trigger is only DISKs " << endl;}

  if(trig_!=0 && trig_!=-1 && trig_!=1 && trig2_!=0 && trig2_!=-1 && trig2_!=1)
    {
      cout <<endl << endl <<endl<< "WARNING!! BSC trigger/scintillator type not properly defined " << endl;
      cout <<endl << endl <<endl<< "WARNING!! BSC trigger/scintillator type not properly defined " << endl << endl << endl;
    }

}

bool BHFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup)
{


  edm::Handle<edm::HepMCProduct>HepMCEvt;
  iEvent.getByLabel("generator","",HepMCEvt);
  const HepMC::GenEvent* MCEvt = HepMCEvt->GetEvent();

  BaseParticlePropagator PP;
  map<unsigned,XYZTLorentzVector> myHits;



  for(HepMC::GenEvent::particle_const_iterator i=MCEvt->particles_begin(); i != MCEvt->particles_end();++i)
    {



      pad_plus = false;
      pad_minus = false;
      circ_plus = false;
      circ_minus = false;

      int myId = (*i)->pdg_id();


        const HepMC::GenVertex * vertex_=(*i)->production_vertex();
        double xv = vertex_->position().x();
        double yv = vertex_->position().y();
        double zv = vertex_->position().z();
        double tv = vertex_->position().t();
        XYZTLorentzVector vertex = XYZTLorentzVector(xv,yv,zv,tv);

        //vertex_->position();

        //      HepMC::FourVector 
        XYZTLorentzVector momentum = XYZTLorentzVector(((*i)->momentum()).x(), ((*i)->momentum()).y(), ((*i)->momentum()).z(), ((*i)->momentum()).t());

        RawParticle myMuon(-momentum, vertex/10.);

          if ( myId < 0 ) 
            myMuon.setCharge(-1.);
          else
            myMuon.setCharge(+1.);

          BaseParticlePropagator PP(myMuon,0.,0.,0.);
          

          float x0= vertex.x();
          float y0= vertex.y();
          float z0= vertex.z();
          float px0=momentum.x();
          float py0=momentum.y();
          float pz0=momentum.z();

          // beam 1 or 2 ?
          // propagator -> need to be implemented 
          

              // intersection point - particle/BSC1
              float zs = 10860.;
              float ys = (zs-z0)*(py0/pz0)+y0;
              float xs = (ys-y0)*(px0/py0)+x0;
              
              
              if(xs<0 && ys>0) {xs=-xs;}
              if(xs>0 && ys<0) {ys=-ys;}
              if(xs<0 && ys<0) {xs=-xs; ys=-ys;}
              
              
              // scintillator pads
              float A[2]={732.7, 24.7};
              float B[2]={895.3, 187.3};
              float C[2]={144.0, 850.2};
              float D[2]={69.8, 776.0};
              float m1=(B[1]-A[1])/(B[0]-A[0]);
              float m2=(C[1]-B[1])/(C[0]-B[0]);
              float m3=(D[1]-C[1])/(D[0]-C[0]);
              float m4=(A[1]-D[1])/(A[0]-D[0]);
              float y1=m1*(xs-A[0])+A[1];
              float y2=m2*(xs-B[0])+B[1];
              float y3=m3*(xs-C[0])+C[1];
              float y4=m4*(xs-D[0])+D[1];
              
              
              
              
              // trigger conditions
              
              if(ys>y1 && ys<y2 && ys<y3 && ys>y4)
                {
                  pad_plus = true;
                  //  cout << " trig1+" << endl;
                }
              
              if((ys<sqrt(450*450-xs*xs) && ys>sqrt(208*208-xs*xs)) || (ys<sqrt(450*450-xs*xs) && xs>208)) 
                {
                  circ_plus = true;
                  //  cout << " trig2+" << endl;
                }
              
              
              // intersection point - particle/BSC2
              zs = -10860.;
              ys = (zs-z0)*(py0/pz0)+y0;
              xs = (ys-y0)*(px0/py0)+x0;
              
              //  cout << endl << " xs ys zs " << xs << " " << ys << " " << zs; 
              
              
              if(xs<0 && ys>0) {xs=-xs;}
              if(xs>0 && ys<0) {ys=-ys;}
              if(xs<0 && ys<0) {xs=-xs; ys=-ys;}
              
              
              // scintillator pads
              y1=m1*(xs-A[0])+A[1];
              y2=m2*(xs-B[0])+B[1];
              y3=m3*(xs-C[0])+C[1];
              y4=m4*(xs-D[0])+D[1];
              
              
              
              // trigger conditions
              
              if(ys>y1 && ys<y2 && ys<y3 && ys>y4)
                {
                  pad_minus = true;
                  // cout << " trig1-" << endl;
                }
              
              if((ys<sqrt(450*450-xs*xs) && ys>sqrt(208*208-xs*xs)) || (ys<sqrt(450*450-xs*xs) && xs>208)) 
                {
                  circ_minus = true;
                  //   cout << " trig2-" << endl;
                }
              
              
              // final selection

              if(trig2_==-1) 
                {
                  pad_plus = false;
                  pad_minus = false;
                }
              if(trig2_==1) 
                {
                  circ_plus = false;
                  circ_minus = false; 
                }
              
              if(trig_==0 && (pad_plus || circ_plus) && (pad_minus || circ_minus) )
                {
                  //              cout << "triggg 0 " << endl;
                  return true;
                }
              if(trig_==1 && (pad_plus || circ_plus)) 
                {
                  //              cout << "triggg 1 " << endl;
                  return true;
                }
              if(trig_==-1 && (pad_minus || circ_minus)) 
                {
                  //              cout << "triggg -1 " << endl;
                  return true;
                }

    }
          
          return false;

}



}