Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /// \file CosmicGenFilterHelix.cc
0002 //
0003 // Original Author:  Gero FLUCKE
0004 //         Created:  Mon Mar  5 16:32:01 CET 2007
0005 
0006 #include "GeneratorInterface/GenFilters/plugins/CosmicGenFilterHelix.h"
0007 
0008 // user include files
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0015 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0016 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0017 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0018 #include "CommonTools/Utils/interface/TFileDirectory.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 
0021 #include <TMath.h>
0022 #include <TH1.h>  // include checker: don't touch!
0023 #include <TH2.h>  // include checker: don't touch!
0024 
0025 #include <utility>  // for std::pair
0026 #include <string>
0027 
0028 //
0029 // constructors and destructor
0030 //
0031 
0032 ///////////////////////////////////////////////////////////////////////////////////////////////////
0033 CosmicGenFilterHelix::CosmicGenFilterHelix(const edm::ParameterSet &cfg)
0034     : theMFToken(esConsumes()),
0035       thePropToken(esConsumes(edm::ESInputTag("", cfg.getParameter<std::string>("propagator")))),
0036       theIds(cfg.getParameter<std::vector<int> >("pdgIds")),
0037       theCharges(cfg.getParameter<std::vector<int> >("charges")),
0038       theMinP2(cfg.getParameter<double>("minP") * cfg.getParameter<double>("minP")),
0039       theMinPt2(cfg.getParameter<double>("minPt") * cfg.getParameter<double>("minPt")),
0040       theDoMonitor(cfg.getUntrackedParameter<bool>("doMonitor")) {
0041   theSrcToken = consumes<edm::HepMCProduct>(cfg.getParameter<edm::InputTag>("src"));
0042 
0043   usesResource(TFileService::kSharedResource);
0044 
0045   if (theIds.size() != theCharges.size()) {
0046     throw cms::Exception("BadConfig") << "CosmicGenFilterHelix: "
0047                                       << "'pdgIds' and 'charges' need same length.";
0048   }
0049   Surface::Scalar radius = cfg.getParameter<double>("radius");
0050   Surface::Scalar maxZ = cfg.getParameter<double>("maxZ");
0051   Surface::Scalar minZ = cfg.getParameter<double>("minZ");
0052 
0053   if (maxZ < minZ) {
0054     throw cms::Exception("BadConfig") << "CosmicGenFilterHelix: maxZ (" << maxZ << ") smaller than minZ (" << minZ
0055                                       << ").";
0056   }
0057 
0058   const Surface::RotationType dummyRot;
0059   theTargetCylinder = Cylinder::build(Surface::PositionType(0., 0., 0.), dummyRot, radius);
0060   theTargetPlaneMin = Plane::build(Surface::PositionType(0., 0., minZ), dummyRot);
0061   theTargetPlaneMax = Plane::build(Surface::PositionType(0., 0., maxZ), dummyRot);
0062 }
0063 
0064 ///////////////////////////////////////////////////////////////////////////////////////////////////
0065 CosmicGenFilterHelix::~CosmicGenFilterHelix() {}
0066 
0067 //
0068 // member functions
0069 //
0070 
0071 ///////////////////////////////////////////////////////////////////////////////////////////////////
0072 bool CosmicGenFilterHelix::filter(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0073   edm::Handle<edm::HepMCProduct> hepMCEvt;
0074   iEvent.getByToken(theSrcToken, hepMCEvt);
0075   const HepMC::GenEvent *mCEvt = hepMCEvt->GetEvent();
0076   const MagneticField *bField = this->getMagneticField(iSetup);  // should be fast (?)
0077   const Propagator *propagator = this->getPropagator(iSetup);
0078 
0079   ++theNumTotal;
0080   bool result = false;
0081   for (HepMC::GenEvent::particle_const_iterator iPart = mCEvt->particles_begin(), endPart = mCEvt->particles_end();
0082        iPart != endPart;
0083        ++iPart) {
0084     int charge = 0;  // there is no method providing charge in GenParticle :-(
0085     if ((*iPart)->status() != 1)
0086       continue;  // look only at stable particles
0087     if (!this->charge((*iPart)->pdg_id(), charge))
0088       continue;
0089 
0090     // Get the position and momentum
0091     const HepMC::ThreeVector hepVertex((*iPart)->production_vertex()->point3d());
0092     const GlobalPoint vert(hepVertex.x() / 10., hepVertex.y() / 10., hepVertex.z() / 10.);  // to cm
0093     const HepMC::FourVector hepMomentum((*iPart)->momentum());
0094     const GlobalVector mom(hepMomentum.x(), hepMomentum.y(), hepMomentum.z());
0095 
0096     if (theDoMonitor)
0097       this->monitorStart(vert, mom, charge, theHistsBefore);
0098 
0099     if (this->propagateToCutCylinder(vert, mom, charge, bField, propagator)) {
0100       result = true;
0101     }
0102   }
0103 
0104   if (result)
0105     ++theNumPass;
0106   return result;
0107 }
0108 
0109 //_________________________________________________________________________________________________
0110 bool CosmicGenFilterHelix::propagateToCutCylinder(const GlobalPoint &vertStart,
0111                                                   const GlobalVector &momStart,
0112                                                   int charge,
0113                                                   const MagneticField *field,
0114                                                   const Propagator *propagator) {
0115   typedef std::pair<TrajectoryStateOnSurface, double> TsosPath;
0116 
0117   const FreeTrajectoryState fts(GlobalTrajectoryParameters(vertStart, momStart, charge, field));
0118 
0119   bool result = true;
0120   TsosPath aTsosPath(propagator->propagateWithPath(fts, *theTargetCylinder));
0121   if (!aTsosPath.first.isValid()) {
0122     result = false;
0123   } else if (aTsosPath.first.globalPosition().z() < theTargetPlaneMin->position().z()) {
0124     // If on cylinder, but outside minimum z, try minimum z-plane:
0125     // (Would it be possible to miss radius on plane, but reach cylinder afterwards in z-range?
0126     //  No, at least not in B-field parallel to z-axis which is cylinder axis.)
0127     aTsosPath = propagator->propagateWithPath(fts, *theTargetPlaneMin);
0128     if (!aTsosPath.first.isValid() || aTsosPath.first.globalPosition().perp() > theTargetCylinder->radius()) {
0129       result = false;
0130     }
0131   } else if (aTsosPath.first.globalPosition().z() > theTargetPlaneMax->position().z()) {
0132     // Analog for outside maximum z:
0133     aTsosPath = propagator->propagateWithPath(fts, *theTargetPlaneMax);
0134     if (!aTsosPath.first.isValid() || aTsosPath.first.globalPosition().perp() > theTargetCylinder->radius()) {
0135       result = false;
0136     }
0137   }
0138 
0139   if (result) {
0140     const GlobalVector momEnd(aTsosPath.first.globalMomentum());
0141     if (momEnd.perp2() < theMinPt2 || momEnd.mag2() < theMinP2) {
0142       result = false;
0143     } else if (theDoMonitor) {
0144       const GlobalPoint vertEnd(aTsosPath.first.globalPosition());
0145       this->monitorStart(vertStart, momStart, charge, theHistsAfter);
0146       this->monitorEnd(vertEnd, momEnd, vertStart, momStart, aTsosPath.second, theHistsAfter);
0147     }
0148   }
0149 
0150   return result;
0151 }
0152 
0153 // ------------ method called once each job just before starting event loop  ------------
0154 void CosmicGenFilterHelix::beginJob() {
0155   if (theDoMonitor) {
0156     this->createHistsStart("start", theHistsBefore);
0157     this->createHistsStart("startAfter", theHistsAfter);
0158     // must be after the line above: hist indices are static in monitorStart(...)
0159     this->createHistsEnd("end", theHistsAfter);
0160   }
0161 
0162   theNumTotal = theNumPass = 0;
0163 }
0164 
0165 ///////////////////////////////////////////////////////////////////////////////////////////////////
0166 void CosmicGenFilterHelix::createHistsStart(const char *dirName, TObjArray &hists) {
0167   edm::Service<TFileService> fs;
0168   TFileDirectory fd(fs->mkdir(dirName, dirName));
0169 
0170   hists.Add(fd.make<TH1F>("momentumP", "|p(#mu^{+})| (start);|p| [GeV]", 100, 0., 1000.));
0171   hists.Add(fd.make<TH1F>("momentumM", "|p(#mu^{-})| (start);|p| [GeV]", 100, 0., 1000.));
0172   hists.Add(fd.make<TH1F>("momentum2", "|p(#mu)| (start);|p| [GeV]", 100, 0., 25.));
0173   const int kNumBins = 50;
0174   double pBinsLog[kNumBins + 1] = {0.};  // fully initialised with 0.
0175   this->equidistLogBins(pBinsLog, kNumBins, 1., 4000.);
0176   hists.Add(fd.make<TH1F>("momentumLog", "|p(#mu)| (start);|p| [GeV]", kNumBins, pBinsLog));
0177   hists.Add(fd.make<TH1F>("phi", "start p_{#phi(#mu)};#phi", 100, -TMath::Pi(), TMath::Pi()));
0178   hists.Add(fd.make<TH1F>("cosPhi", "cos(p_{#phi(#mu)}) (start);cos(#phi)", 100, -1., 1.));
0179   hists.Add(fd.make<TH1F>("phiXz", "start p_{#phi_{xz}(#mu)};#phi_{xz}", 100, -TMath::Pi(), TMath::Pi()));
0180   hists.Add(fd.make<TH1F>("theta", "#theta(#mu) (start);#theta", 100, 0., TMath::Pi()));
0181   hists.Add(fd.make<TH1F>("thetaY", "#theta_{y}(#mu) (start);#theta_{y}", 100, 0., TMath::Pi() / 2.));
0182 
0183   hists.Add(fd.make<TH2F>(
0184       "momVsPhi", "|p(#mu)| vs #phi (start);#phi;|p| [GeV]", 50, -TMath::Pi(), TMath::Pi(), 50, 1.5, 1000.));
0185   hists.Add(
0186       fd.make<TH2F>("momVsTheta", "|p(#mu)| vs #theta (start);#theta;|p| [GeV]", 50, 0., TMath::Pi(), 50, 1.5, 1000.));
0187   hists.Add(fd.make<TH2F>(
0188       "momVsThetaY", "|p(#mu)| vs #theta_{y} (start);#theta_{y};|p| [GeV]", 50, 0., TMath::Pi() / 2., 50, 1.5, 1000.));
0189   hists.Add(fd.make<TH2F>("momVsZ", "|p(#mu)| vs z (start);z [cm];|p| [GeV]", 50, -1600., 1600., 50, 1.5, 1000.));
0190   hists.Add(fd.make<TH2F>("thetaVsZ", "#theta vs z (start);z [cm];#theta", 50, -1600., 1600., 50, 0., TMath::Pi()));
0191   hists.Add(
0192       fd.make<TH2F>("thetaYvsZ", "#theta_{y} vs z (start);z [cm];#theta", 50, -1600., 1600., 50, 0., TMath::PiOver2()));
0193   hists.Add(fd.make<TH2F>(
0194       "yVsThetaY", "#theta_{y}(#mu) vs y (start);#theta_{y};y [cm]", 50, 0., TMath::Pi() / 2., 50, -1000., 1000.));
0195   hists.Add(fd.make<TH2F>("yVsThetaYnoR",
0196                           "#theta_{y}(#mu) vs y (start, barrel);#theta_{y};y [cm]",
0197                           50,
0198                           0.,
0199                           TMath::Pi() / 2.,
0200                           50,
0201                           -1000.,
0202                           1000.));
0203 
0204   hists.Add(fd.make<TH1F>("radius", "start radius;r [cm]", 100, 0., 1000.));
0205   hists.Add(fd.make<TH1F>("z", "start z;z [cm]", 100, -1600., 1600.));
0206   hists.Add(fd.make<TH2F>("xyPlane", "start xy;x [cm];y [cm]", 50, -1000., 1000., 50, -1000., 1000.));
0207   hists.Add(fd.make<TH2F>(
0208       "rzPlane", "start rz (y < 0 #Rightarrow r_{#pm} = -r);z [cm];r_{#pm} [cm]", 50, -1600., 1600., 50, -1000., 1000.));
0209 }
0210 
0211 ///////////////////////////////////////////////////////////////////////////////////////////////////
0212 void CosmicGenFilterHelix::createHistsEnd(const char *dirName, TObjArray &hists) {
0213   edm::Service<TFileService> fs;
0214   TFileDirectory fd(fs->mkdir(dirName, dirName));
0215 
0216   const int kNumBins = 50;
0217   double pBinsLog[kNumBins + 1] = {0.};  // fully initialised with 0.
0218   this->equidistLogBins(pBinsLog, kNumBins, 1., 4000.);
0219 
0220   // take care: hist names must differ from those in createHistsStart!
0221   hists.Add(fd.make<TH1F>("pathEnd", "path until cylinder;s [cm]", 100, 0., 2000.));
0222   hists.Add(fd.make<TH1F>("momEnd", "|p_{end}|;p [GeV]", 100, 0., 1000.));
0223   hists.Add(fd.make<TH1F>("momEndLog", "|p_{end}|;p [GeV]", kNumBins, pBinsLog));
0224   hists.Add(fd.make<TH1F>("ptEnd", "p_{t} (end);p_{t} [GeV]", 100, 0., 750.));
0225   hists.Add(fd.make<TH1F>("ptEndLog", "p_{t} (end);p_{t} [GeV]", kNumBins, pBinsLog));
0226   hists.Add(fd.make<TH1F>("phiXzEnd", "#phi_{xz} (end);#phi_{xz}", 100, -TMath::Pi(), TMath::Pi()));
0227   hists.Add(fd.make<TH1F>("thetaYEnd", "#theta_{y} (end);#theta_{y}", 100, 0., TMath::Pi()));
0228 
0229   hists.Add(fd.make<TH1F>("momStartEnd", "|p_{start}|-|p_{end}|;#Deltap [GeV]", 100, 0., 15.));
0230   hists.Add(fd.make<TH1F>("momStartEndRel", "(p_{start}-p_{end})/p_{start};#Deltap_{rel}", 100, .0, 1.));
0231   hists.Add(fd.make<TH1F>("phiXzStartEnd", "#phi_{xz,start}-#phi_{xz,end};#Delta#phi_{xz}", 100, -1., 1.));
0232   hists.Add(fd.make<TH1F>("thetaYStartEnd", "#theta_{y,start}-#theta_{y,end};#Delta#theta_{y}", 100, -1., 1.));
0233 
0234   hists.Add(fd.make<TH2F>("phiXzStartVsEnd",
0235                           "#phi_{xz} start vs end;#phi_{xz}^{end};#phi_{xz}^{start}",
0236                           50,
0237                           -TMath::Pi(),
0238                           TMath::Pi(),
0239                           50,
0240                           -TMath::Pi(),
0241                           TMath::Pi()));
0242   hists.Add(fd.make<TH2F>("thetaYStartVsEnd",
0243                           "#theta_{y} start vs end;#theta_{y}^{end};#theta_{y}^{start}",
0244                           50,
0245                           0.,
0246                           TMath::Pi(),
0247                           50,
0248                           0.,
0249                           TMath::Pi() / 2.));
0250 
0251   hists.Add(fd.make<TH2F>("momStartEndRelVsZ",
0252                           "(p_{start}-p_{end})/p_{start} vs z_{start};z [cm];#Deltap_{rel}",
0253                           50,
0254                           -1600.,
0255                           1600.,
0256                           50,
0257                           .0,
0258                           .8));
0259   hists.Add(fd.make<TH2F>("phiXzStartEndVsZ",
0260                           "#phi_{xz,start}-#phi_{xz,end} vs z_{start};z [cm];#Delta#phi_{xz}",
0261                           50,
0262                           -1600.,
0263                           1600.,
0264                           50,
0265                           -1.,
0266                           1.));
0267   hists.Add(fd.make<TH2F>("thetaYStartEndVsZ",
0268                           "#theta_{y,start}-#theta_{y,end} vs z_{start};z [cm];#Delta#theta_{y}",
0269                           50,
0270                           -1600.,
0271                           1600.,
0272                           50,
0273                           -.5,
0274                           .5));
0275   hists.Add(fd.make<TH2F>("momStartEndRelVsP",
0276                           "(p_{start}-p_{end})/p_{start} vs p_{start};p [GeV];#Deltap_{rel}",
0277                           kNumBins,
0278                           pBinsLog,
0279                           50,
0280                           .0,
0281                           .8));
0282   hists.Add(fd.make<TH2F>("phiXzStartEndVsP",
0283                           "#phi_{xz,start}-#phi_{xz,end} vs |p|_{start};p [GeV];#Delta#phi_{xz}",
0284                           kNumBins,
0285                           pBinsLog,
0286                           100,
0287                           -1.5,
0288                           1.5));
0289   hists.Add(fd.make<TH2F>("thetaYStartEndVsP",
0290                           "#theta_{y,start}-#theta_{y,end} vs |p|_{start};p [GeV];#Delta#theta_{y}",
0291                           kNumBins,
0292                           pBinsLog,
0293                           100,
0294                           -1.,
0295                           1.));
0296 
0297   const double maxR = theTargetCylinder->radius() * 1.1;
0298   hists.Add(fd.make<TH1F>("radiusEnd", "end radius;r [cm]", 100, 0., maxR));
0299   double minZ = theTargetPlaneMin->position().z();
0300   minZ -= TMath::Abs(minZ) * 0.1;
0301   double maxZ = theTargetPlaneMax->position().z();
0302   maxZ += TMath::Abs(maxZ) * 0.1;
0303   hists.Add(fd.make<TH1F>("zEnd", "end z;z [cm]", 100, minZ, maxZ));
0304   hists.Add(fd.make<TH1F>("zDiff", "z_{start}-z_{end};#Deltaz [cm]", 100, -1000., 1000.));
0305   hists.Add(fd.make<TH2F>("xyPlaneEnd", "end xy;x [cm];y [cm]", 100, -maxR, maxR, 100, -maxR, maxR));
0306 
0307   hists.Add(fd.make<TH2F>(
0308       "rzPlaneEnd", "end rz (y<0 #Rightarrow r_{#pm}=-r);z [cm];r_{#pm} [cm]", 50, minZ, maxZ, 50, -maxR, maxR));
0309   hists.Add(fd.make<TH2F>("thetaVsZend", "#theta vs z (end);z [cm];#theta", 50, minZ, maxZ, 50, 0., TMath::Pi()));
0310 }
0311 
0312 // ------------ method called once each job just after ending the event loop  ------------
0313 void CosmicGenFilterHelix::endJob() {
0314   const char *border = "////////////////////////////////////////////////////////";
0315   const char *line = "\n// ";
0316   edm::LogInfo("Filter") << "@SUB=CosmicGenFilterHelix::endJob" << border << line << theNumPass << " events out of "
0317                          << theNumTotal << ", i.e. " << theNumPass * 100. / theNumTotal << "%, "
0318                          << "reached target cylinder," << line << "defined by r < " << theTargetCylinder->radius()
0319                          << " cm and " << theTargetPlaneMin->position().z() << " < z < "
0320                          << theTargetPlaneMax->position().z() << " cm." << line
0321                          << "Minimal required (transverse) momentum was " << TMath::Sqrt(theMinP2) << " ("
0322                          << TMath::Sqrt(theMinPt2) << ") GeV."
0323                          << "\n"
0324                          << border;
0325 }
0326 
0327 ///////////////////////////////////////////////////////////////////////////////////////////////////
0328 bool CosmicGenFilterHelix::charge(int id, int &charge) const {
0329   std::vector<int>::const_iterator iC = theCharges.begin();
0330   for (std::vector<int>::const_iterator i = theIds.begin(), end = theIds.end(); i != end; ++i, ++iC) {
0331     if (*i == id) {
0332       charge = *iC;
0333       return true;
0334     }
0335   }
0336 
0337   return false;
0338 }
0339 
0340 //_________________________________________________________________________________________________
0341 const MagneticField *CosmicGenFilterHelix::getMagneticField(const edm::EventSetup &setup) const {
0342   return &setup.getData(theMFToken);
0343 }
0344 
0345 //_________________________________________________________________________________________________
0346 const Propagator *CosmicGenFilterHelix::getPropagator(const edm::EventSetup &setup) const {
0347   const Propagator *prop = &setup.getData(thePropToken);
0348   if (!dynamic_cast<const SteppingHelixPropagator *>(prop)) {
0349     edm::LogWarning("BadConfig") << "@SUB=CosmicGenFilterHelix::getPropagator"
0350                                  << "Not a SteppingHelixPropagator!";
0351   }
0352   return prop;
0353 }
0354 
0355 //_________________________________________________________________________________________________
0356 void CosmicGenFilterHelix::monitorStart(const GlobalPoint &vert, const GlobalVector &mom, int charge, TObjArray &hists) {
0357   const double scalarMom = mom.mag();
0358   const double phi = mom.phi();
0359   const double phiXz = TMath::ATan2(mom.z(), mom.x());
0360   const double theta = mom.theta();
0361   const double thetaY = TMath::ATan2(TMath::Sqrt(mom.x() * mom.x() + mom.z() * mom.z()), -mom.y());
0362 
0363   const double z = vert.z();
0364   const double r = vert.perp();
0365 
0366   const static int iMomP = hists.IndexOf(hists.FindObject("momentumP"));
0367   const static int iMomM = hists.IndexOf(hists.FindObject("momentumM"));
0368   if (charge > 0)
0369     static_cast<TH1 *>(hists[iMomP])->Fill(scalarMom);
0370   else
0371     static_cast<TH1 *>(hists[iMomM])->Fill(scalarMom);
0372   const static int iMom2 = hists.IndexOf(hists.FindObject("momentum2"));
0373   static_cast<TH1 *>(hists[iMom2])->Fill(scalarMom);
0374   const static int iMomLog = hists.IndexOf(hists.FindObject("momentumLog"));
0375   static_cast<TH1 *>(hists[iMomLog])->Fill(scalarMom);
0376   const static int iPhi = hists.IndexOf(hists.FindObject("phi"));
0377   static_cast<TH1 *>(hists[iPhi])->Fill(phi);
0378   const static int iCosPhi = hists.IndexOf(hists.FindObject("cosPhi"));
0379   static_cast<TH1 *>(hists[iCosPhi])->Fill(TMath::Cos(phi));
0380   const static int iPhiXz = hists.IndexOf(hists.FindObject("phiXz"));
0381   static_cast<TH1 *>(hists[iPhiXz])->Fill(phiXz);
0382   const static int iTheta = hists.IndexOf(hists.FindObject("theta"));
0383   static_cast<TH1 *>(hists[iTheta])->Fill(theta);
0384   const static int iThetaY = hists.IndexOf(hists.FindObject("thetaY"));
0385   static_cast<TH1 *>(hists[iThetaY])->Fill(thetaY);
0386 
0387   const static int iMomVsTheta = hists.IndexOf(hists.FindObject("momVsTheta"));
0388   static_cast<TH2 *>(hists[iMomVsTheta])->Fill(theta, scalarMom);
0389   const static int iMomVsThetaY = hists.IndexOf(hists.FindObject("momVsThetaY"));
0390   static_cast<TH2 *>(hists[iMomVsThetaY])->Fill(thetaY, scalarMom);
0391   const static int iMomVsPhi = hists.IndexOf(hists.FindObject("momVsPhi"));
0392   static_cast<TH2 *>(hists[iMomVsPhi])->Fill(phi, scalarMom);
0393   const static int iMomVsZ = hists.IndexOf(hists.FindObject("momVsZ"));
0394   static_cast<TH2 *>(hists[iMomVsZ])->Fill(z, scalarMom);
0395   const static int iThetaVsZ = hists.IndexOf(hists.FindObject("thetaVsZ"));
0396   static_cast<TH2 *>(hists[iThetaVsZ])->Fill(z, theta);
0397   const static int iThetaYvsZ = hists.IndexOf(hists.FindObject("thetaYvsZ"));
0398   static_cast<TH2 *>(hists[iThetaYvsZ])->Fill(z, thetaY);
0399   const static int iYvsThetaY = hists.IndexOf(hists.FindObject("yVsThetaY"));
0400   static_cast<TH2 *>(hists[iYvsThetaY])->Fill(thetaY, vert.y());
0401   const static int iYvsThetaYnoR = hists.IndexOf(hists.FindObject("yVsThetaYnoR"));
0402   if (z > -1400. && z < 1400.) {
0403     static_cast<TH2 *>(hists[iYvsThetaYnoR])->Fill(thetaY, vert.y());
0404   }
0405 
0406   const static int iRadius = hists.IndexOf(hists.FindObject("radius"));
0407   static_cast<TH1 *>(hists[iRadius])->Fill(r);
0408   const static int iZ = hists.IndexOf(hists.FindObject("z"));
0409   static_cast<TH1 *>(hists[iZ])->Fill(z);
0410   const static int iXy = hists.IndexOf(hists.FindObject("xyPlane"));
0411   static_cast<TH1 *>(hists[iXy])->Fill(vert.x(), vert.y());
0412   const static int iRz = hists.IndexOf(hists.FindObject("rzPlane"));
0413   static_cast<TH1 *>(hists[iRz])->Fill(z, (vert.y() > 0 ? r : -r));
0414 }
0415 
0416 //_________________________________________________________________________________________________
0417 void CosmicGenFilterHelix::monitorEnd(const GlobalPoint &endVert,
0418                                       const GlobalVector &endMom,
0419                                       const GlobalPoint &vert,
0420                                       const GlobalVector &mom,
0421                                       double path,
0422                                       TObjArray &hists) {
0423   const double scalarMomStart = mom.mag();
0424   const double phiXzStart = TMath::ATan2(mom.z(), mom.x());
0425   const double thetaYStart = TMath::ATan2(TMath::Sqrt(mom.x() * mom.x() + mom.z() * mom.z()), -mom.y());
0426   const double scalarMomEnd = endMom.mag();
0427   const double ptEnd = endMom.perp();
0428   const double phiXzEnd = TMath::ATan2(endMom.z(), endMom.x());
0429   const double thetaYEnd = TMath::ATan2(TMath::Sqrt(endMom.x() * endMom.x() + endMom.z() * endMom.z()), -endMom.y());
0430   const double thetaEnd = endMom.theta();
0431 
0432   const double diffMomRel = (scalarMomStart - scalarMomEnd) / scalarMomStart;
0433 
0434   const double zEnd = endVert.z();
0435   const double rEnd = endVert.perp();
0436   const double diffZ = zEnd - vert.z();
0437 
0438   const static int iPathEnd = hists.IndexOf(hists.FindObject("pathEnd"));
0439   static_cast<TH1 *>(hists[iPathEnd])->Fill(path);
0440   const static int iMomEnd = hists.IndexOf(hists.FindObject("momEnd"));
0441   static_cast<TH1 *>(hists[iMomEnd])->Fill(scalarMomEnd);
0442   const static int iMomEndLog = hists.IndexOf(hists.FindObject("momEndLog"));
0443   static_cast<TH1 *>(hists[iMomEndLog])->Fill(scalarMomEnd);
0444   const static int iPtEnd = hists.IndexOf(hists.FindObject("ptEnd"));
0445   static_cast<TH1 *>(hists[iPtEnd])->Fill(ptEnd);
0446   const static int iPtEndLog = hists.IndexOf(hists.FindObject("ptEndLog"));
0447   static_cast<TH1 *>(hists[iPtEndLog])->Fill(ptEnd);
0448   const static int iPhiXzEnd = hists.IndexOf(hists.FindObject("phiXzEnd"));
0449   static_cast<TH1 *>(hists[iPhiXzEnd])->Fill(phiXzEnd);
0450   const static int iThetaYEnd = hists.IndexOf(hists.FindObject("thetaYEnd"));
0451   static_cast<TH1 *>(hists[iThetaYEnd])->Fill(thetaYEnd);
0452 
0453   const static int iMomStartEnd = hists.IndexOf(hists.FindObject("momStartEnd"));
0454   static_cast<TH1 *>(hists[iMomStartEnd])->Fill(scalarMomStart - scalarMomEnd);
0455   const static int iMomStartEndRel = hists.IndexOf(hists.FindObject("momStartEndRel"));
0456   static_cast<TH1 *>(hists[iMomStartEndRel])->Fill(diffMomRel);
0457   const static int iPhiStartEnd = hists.IndexOf(hists.FindObject("phiXzStartEnd"));
0458   static_cast<TH1 *>(hists[iPhiStartEnd])->Fill(phiXzStart - phiXzEnd);
0459   const static int iThetaStartEnd = hists.IndexOf(hists.FindObject("thetaYStartEnd"));
0460   static_cast<TH1 *>(hists[iThetaStartEnd])->Fill(thetaYStart - thetaYEnd);
0461 
0462   const static int iPhiStartVsEnd = hists.IndexOf(hists.FindObject("phiXzStartVsEnd"));
0463   static_cast<TH2 *>(hists[iPhiStartVsEnd])->Fill(phiXzEnd, phiXzStart);
0464   const static int iThetaStartVsEnd = hists.IndexOf(hists.FindObject("thetaYStartVsEnd"));
0465   static_cast<TH2 *>(hists[iThetaStartVsEnd])->Fill(thetaYEnd, thetaYStart);
0466 
0467   const static int iMomStartEndRelVsZ = hists.IndexOf(hists.FindObject("momStartEndRelVsZ"));
0468   static_cast<TH2 *>(hists[iMomStartEndRelVsZ])->Fill(vert.z(), diffMomRel);
0469   const static int iPhiStartEndVsZ = hists.IndexOf(hists.FindObject("phiXzStartEndVsZ"));
0470   static_cast<TH2 *>(hists[iPhiStartEndVsZ])->Fill(vert.z(), phiXzStart - phiXzEnd);
0471   const static int iThetaStartEndVsZ = hists.IndexOf(hists.FindObject("thetaYStartEndVsZ"));
0472   static_cast<TH2 *>(hists[iThetaStartEndVsZ])->Fill(vert.z(), thetaYStart - thetaYEnd);
0473   const static int iMomStartEndRelVsP = hists.IndexOf(hists.FindObject("momStartEndRelVsP"));
0474   static_cast<TH2 *>(hists[iMomStartEndRelVsP])->Fill(scalarMomStart, diffMomRel);
0475   const static int iPhiStartEndVsP = hists.IndexOf(hists.FindObject("phiXzStartEndVsP"));
0476   static_cast<TH2 *>(hists[iPhiStartEndVsP])->Fill(scalarMomStart, phiXzStart - phiXzEnd);
0477   const static int iThetaStartEndVsP = hists.IndexOf(hists.FindObject("thetaYStartEndVsP"));
0478   static_cast<TH2 *>(hists[iThetaStartEndVsP])->Fill(scalarMomStart, thetaYStart - thetaYEnd);
0479 
0480   const static int iRadiusEnd = hists.IndexOf(hists.FindObject("radiusEnd"));
0481   static_cast<TH1 *>(hists[iRadiusEnd])->Fill(rEnd);
0482   const static int iZend = hists.IndexOf(hists.FindObject("zEnd"));
0483   static_cast<TH1 *>(hists[iZend])->Fill(zEnd);
0484   const static int iZdiff = hists.IndexOf(hists.FindObject("zDiff"));
0485   static_cast<TH1 *>(hists[iZdiff])->Fill(diffZ);
0486   const static int iXyPlaneEnd = hists.IndexOf(hists.FindObject("xyPlaneEnd"));
0487   static_cast<TH1 *>(hists[iXyPlaneEnd])->Fill(endVert.x(), endVert.y());
0488   const static int iRzPlaneEnd = hists.IndexOf(hists.FindObject("rzPlaneEnd"));
0489   static_cast<TH1 *>(hists[iRzPlaneEnd])->Fill(zEnd, (endVert.y() > 0 ? rEnd : -rEnd));
0490   const static int iThetaVsZend = hists.IndexOf(hists.FindObject("thetaVsZend"));
0491   static_cast<TH2 *>(hists[iThetaVsZend])->Fill(zEnd, thetaEnd);
0492 }
0493 
0494 //_________________________________________________________________________________________________
0495 bool CosmicGenFilterHelix::equidistLogBins(double *bins, int nBins, double first, double last) const {
0496   // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last'
0497   // that are equidistant when viewed in log scale,
0498   // so 'bins' must have length nBins+1;
0499   // If 'first', 'last' or 'nBins' are not positive, failure is reported.
0500 
0501   if (nBins < 1 || first <= 0. || last <= 0.)
0502     return false;
0503 
0504   bins[0] = first;
0505   bins[nBins] = last;
0506   const double firstLog = TMath::Log10(bins[0]);
0507   const double lastLog = TMath::Log10(bins[nBins]);
0508   for (int i = 1; i < nBins; ++i) {
0509     bins[i] = TMath::Power(10., firstLog + i * (lastLog - firstLog) / (nBins));
0510   }
0511 
0512   return true;
0513 }
[ 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!