Reference Manual

Go to the documentation of this file.
00001 
00002 //
00003 // Original Author:  Gero FLUCKE
00004 //         Created:  Mon Mar  5 16:32:01 CET 2007
00005 // $Id: CosmicGenFilterHelix.cc,v 1.12 2010/01/05 13:49:07 hegner Exp $
00006 
00007 #include "GeneratorInterface/GenFilters/interface/CosmicGenFilterHelix.h"
00008 
00009 // user include files
00010 #include "FWCore/Framework/interface/ESHandle.h"
00011 #include "DataFormats/Common/interface/Handle.h"
00012 
00013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00014 
00015 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00016 
00017 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00018 
00019 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00020 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00021 
00022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00023 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00024 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00025 
00026 
00027 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00028 #include "CommonTools/Utils/interface/TFileDirectory.h"
00029 #include "FWCore/ServiceRegistry/interface/Service.h"
00030 
00031 #include <TMath.h>
00032 #include <TH1.h> // include checker: don't touch!
00033 #include <TH2.h> // include checker: don't touch!
00034 
00035 #include <utility> // for std::pair
00036 #include <string>
00037 
00038 
00039 //
00040 // constructors and destructor
00041 //
00042 
00044 CosmicGenFilterHelix::CosmicGenFilterHelix(const edm::ParameterSet& cfg) 
00045   : theSrc(cfg.getParameter<edm::InputTag>("src")),
00046     theIds(cfg.getParameter<std::vector<int> >("pdgIds")),
00047     theCharges(cfg.getParameter<std::vector<int> >("charges")),
00048     thePropagatorName(cfg.getParameter<std::string>("propagator")),
00049     theMinP2(cfg.getParameter<double>("minP")*cfg.getParameter<double>("minP")),
00050     theMinPt2(cfg.getParameter<double>("minPt")*cfg.getParameter<double>("minPt")),
00051     theDoMonitor(cfg.getUntrackedParameter<bool>("doMonitor"))
00052 {
00053   if (theIds.size() != theCharges.size()) {
00054     throw cms::Exception("BadConfig") << "CosmicGenFilterHelix: "
00055                                       << "'pdgIds' and 'charges' need same length.";
00056   }
00057   Surface::Scalar radius = cfg.getParameter<double>("radius");
00058   Surface::Scalar maxZ   = cfg.getParameter<double>("maxZ");
00059   Surface::Scalar minZ   = cfg.getParameter<double>("minZ");
00060   
00061   if (maxZ < minZ) {
00062     throw cms::Exception("BadConfig") << "CosmicGenFilterHelix: maxZ (" << maxZ
00063                                       << ") smaller than minZ (" << minZ << ").";
00064   }
00065   
00066   const Surface::RotationType dummyRot;
00067   theTargetCylinder = Cylinder::build(Surface::PositionType(0.,0.,0.), dummyRot, radius);
00068   theTargetPlaneMin = Plane::build(Surface::PositionType(0.,0.,minZ), dummyRot);
00069   theTargetPlaneMax = Plane::build(Surface::PositionType(0.,0.,maxZ), dummyRot);
00070 }
00071 
00073 CosmicGenFilterHelix::~CosmicGenFilterHelix()
00074 {
00075 }
00076 
00077 
00078 //
00079 // member functions
00080 //
00081 
00083 bool CosmicGenFilterHelix::filter(edm::Event &iEvent, const edm::EventSetup &iSetup)
00084 {
00085   edm::Handle<edm::HepMCProduct> hepMCEvt;
00086   iEvent.getByLabel(theSrc, hepMCEvt);
00087   const HepMC::GenEvent *mCEvt = hepMCEvt->GetEvent();
00088   const MagneticField *bField = this->getMagneticField(iSetup); // should be fast (?)
00089   const Propagator *propagator = this->getPropagator(iSetup);
00090 
00091   ++theNumTotal;
00092   bool result = false;
00093   for (HepMC::GenEvent::particle_const_iterator iPart = mCEvt->particles_begin(),
00094          endPart = mCEvt->particles_end(); iPart != endPart; ++iPart) {
00095     int charge = 0; // there is no method providing charge in GenParticle :-(
00096     if ((*iPart)->status() != 1) continue; // look only at stable particles
00097     if (!this->charge((*iPart)->pdg_id(), charge)) continue;
00098 
00099     // Get the position and momentum
00100     const HepMC::ThreeVector hepVertex((*iPart)->production_vertex()->point3d());
00101     const GlobalPoint vert(hepVertex.x()/10., hepVertex.y()/10., hepVertex.z()/10.); // to cm
00102     const HepMC::FourVector hepMomentum((*iPart)->momentum());
00103     const GlobalVector mom(hepMomentum.x(), hepMomentum.y(), hepMomentum.z());
00104 
00105     if (theDoMonitor) this->monitorStart(vert, mom, charge, theHistsBefore);
00106 
00107     if (this->propagateToCutCylinder(vert, mom, charge, bField, propagator)) {
00108       result = true;
00109     }
00110   }
00111 
00112   if (result) ++theNumPass;
00113   return result;
00114 }
00115 
00116 //_________________________________________________________________________________________________
00117 bool CosmicGenFilterHelix::propagateToCutCylinder(const GlobalPoint &vertStart,
00118                                                   const GlobalVector &momStart,
00119                                                   int charge, const MagneticField *field,
00120                                                   const Propagator *propagator)
00121 {
00122   typedef std::pair<TrajectoryStateOnSurface, double> TsosPath;
00123 
00124   const FreeTrajectoryState fts(GlobalTrajectoryParameters(vertStart, momStart, charge, field));
00125 
00126   bool result = true;
00127   TsosPath aTsosPath(propagator->propagateWithPath(fts, *theTargetCylinder));
00128   if (!aTsosPath.first.isValid()) {
00129     result = false;
00130   } else if (aTsosPath.first.globalPosition().z() < theTargetPlaneMin->position().z()) {
00131     // If on cylinder, but outside minimum z, try minimum z-plane:
00132     // (Would it be possible to miss radius on plane, but reach cylinder afterwards in z-range?
00133     //  No, at least not in B-field parallel to z-axis which is cylinder axis.)
00134     aTsosPath = propagator->propagateWithPath(fts, *theTargetPlaneMin);
00135     if (!aTsosPath.first.isValid()
00136         || aTsosPath.first.globalPosition().perp() > theTargetCylinder->radius()) {
00137       result = false;
00138     }
00139   } else if (aTsosPath.first.globalPosition().z() > theTargetPlaneMax->position().z()) {
00140     // Analog for outside maximum z:
00141     aTsosPath = propagator->propagateWithPath(fts, *theTargetPlaneMax);
00142     if (!aTsosPath.first.isValid()
00143         || aTsosPath.first.globalPosition().perp() > theTargetCylinder->radius()) {
00144       result = false;
00145     }
00146   }
00147 
00148   if (result) {
00149     const GlobalVector momEnd(aTsosPath.first.globalMomentum());
00150     if (momEnd.perp2() < theMinPt2 || momEnd.mag2() < theMinP2) {
00151       result = false;
00152     } else if (theDoMonitor) {
00153       const GlobalPoint vertEnd(aTsosPath.first.globalPosition());
00154       this->monitorStart(vertStart, momStart, charge, theHistsAfter);
00155       this->monitorEnd(vertEnd, momEnd, vertStart, momStart, aTsosPath.second, theHistsAfter);
00156     }
00157   }
00158 
00159   return result;
00160 }
00161 
00162 
00163 // ------------ method called once each job just before starting event loop  ------------
00164 void CosmicGenFilterHelix::beginJob()
00165 {
00166   if (theDoMonitor) {
00167     this->createHistsStart("start", theHistsBefore);
00168     this->createHistsStart("startAfter", theHistsAfter);
00169     // must be after the line above: hist indices are static in monitorStart(...)
00170     this->createHistsEnd("end", theHistsAfter);
00171   }
00172 
00173   theNumTotal = theNumPass = 0;
00174 }
00175 
00177 void CosmicGenFilterHelix::createHistsStart(const char *dirName, TObjArray &hists)
00178 {
00179   edm::Service<TFileService> fs;
00180   TFileDirectory fd(fs->mkdir(dirName, dirName));
00181   
00182   hists.Add(fd.make<TH1F>("momentumP", "|p(#mu^{+})| (start);|p| [GeV]",100, 0., 1000.));
00183   hists.Add(fd.make<TH1F>("momentumM", "|p(#mu^{-})| (start);|p| [GeV]",100, 0., 1000.));
00184   hists.Add(fd.make<TH1F>("momentum2", "|p(#mu)| (start);|p| [GeV]",100, 0., 25.));
00185   const int kNumBins = 50;
00186   double pBinsLog[kNumBins+1] = {0.}; // fully initialised with 0.
00187   this->equidistLogBins(pBinsLog, kNumBins, 1., 4000.);
00188   hists.Add(fd.make<TH1F>("momentumLog", "|p(#mu)| (start);|p| [GeV]", kNumBins, pBinsLog));
00189   hists.Add(fd.make<TH1F>("phi", "start p_{#phi(#mu)};#phi", 100, -TMath::Pi(), TMath::Pi()));
00190   hists.Add(fd.make<TH1F>("cosPhi", "cos(p_{#phi(#mu)}) (start);cos(#phi)", 100, -1., 1.));
00191   hists.Add(fd.make<TH1F>("phiXz", "start p_{#phi_{xz}(#mu)};#phi_{xz}",
00192                           100, -TMath::Pi(), TMath::Pi()));
00193   hists.Add(fd.make<TH1F>("theta", "#theta(#mu) (start);#theta", 100, 0., TMath::Pi()));
00194   hists.Add(fd.make<TH1F>("thetaY", "#theta_{y}(#mu) (start);#theta_{y}", 100,0.,TMath::Pi()/2.));
00195   
00196   hists.Add(fd.make<TH2F>("momVsPhi", "|p(#mu)| vs #phi (start);#phi;|p| [GeV]",
00197                           50, -TMath::Pi(), TMath::Pi(), 50, 1.5, 1000.));
00198   hists.Add(fd.make<TH2F>("momVsTheta", "|p(#mu)| vs #theta (start);#theta;|p| [GeV]",
00199                           50, 0., TMath::Pi(), 50, 1.5, 1000.));
00200   hists.Add(fd.make<TH2F>("momVsThetaY", "|p(#mu)| vs #theta_{y} (start);#theta_{y};|p| [GeV]",
00201                           50, 0., TMath::Pi()/2., 50, 1.5, 1000.));
00202   hists.Add(fd.make<TH2F>("momVsZ", "|p(#mu)| vs z (start);z [cm];|p| [GeV]",
00203                           50, -1600., 1600., 50, 1.5, 1000.));
00204   hists.Add(fd.make<TH2F>("thetaVsZ", "#theta vs z (start);z [cm];#theta",
00205                           50, -1600., 1600., 50, 0., TMath::Pi()));
00206   hists.Add(fd.make<TH2F>("thetaYvsZ", "#theta_{y} vs z (start);z [cm];#theta",
00207                           50, -1600., 1600., 50, 0., TMath::PiOver2()));
00208   hists.Add(fd.make<TH2F>("yVsThetaY", "#theta_{y}(#mu) vs y (start);#theta_{y};y [cm]",
00209                           50, 0., TMath::Pi()/2., 50, -1000., 1000.));
00210   hists.Add(fd.make<TH2F>("yVsThetaYnoR", "#theta_{y}(#mu) vs y (start, barrel);#theta_{y};y [cm]",
00211                           50, 0., TMath::Pi()/2., 50, -1000., 1000.));
00212   
00213   hists.Add(fd.make<TH1F>("radius", "start radius;r [cm]", 100, 0., 1000.));
00214   hists.Add(fd.make<TH1F>("z", "start z;z [cm]", 100, -1600., 1600.));
00215   hists.Add(fd.make<TH2F>("xyPlane", "start xy;x [cm];y [cm]", 50, -1000., 1000.,
00216                           50, -1000., 1000.));
00217   hists.Add(fd.make<TH2F>("rzPlane",
00218                           "start rz (y < 0 #Rightarrow r_{#pm} = -r);z [cm];r_{#pm} [cm]",
00219                           50, -1600., 1600., 50, -1000., 1000.));
00220 }
00221 
00223 void CosmicGenFilterHelix::createHistsEnd(const char *dirName, TObjArray &hists)
00224 {
00225   edm::Service<TFileService> fs;
00226   TFileDirectory fd(fs->mkdir(dirName, dirName));
00227 
00228   const int kNumBins = 50;
00229   double pBinsLog[kNumBins+1] = {0.}; // fully initialised with 0.
00230   this->equidistLogBins(pBinsLog, kNumBins, 1., 4000.);
00231 
00232   // take care: hist names must differ from those in createHistsStart!
00233   hists.Add(fd.make<TH1F>("pathEnd", "path until cylinder;s [cm]", 100, 0., 2000.));
00234   hists.Add(fd.make<TH1F>("momEnd", "|p_{end}|;p [GeV]", 100, 0., 1000.));
00235   hists.Add(fd.make<TH1F>("momEndLog", "|p_{end}|;p [GeV]", kNumBins, pBinsLog));
00236   hists.Add(fd.make<TH1F>("ptEnd", "p_{t} (end);p_{t} [GeV]", 100, 0., 750.));
00237   hists.Add(fd.make<TH1F>("ptEndLog", "p_{t} (end);p_{t} [GeV]", kNumBins, pBinsLog));
00238   hists.Add(fd.make<TH1F>("phiXzEnd", "#phi_{xz} (end);#phi_{xz}", 100,-TMath::Pi(),TMath::Pi()));
00239   hists.Add(fd.make<TH1F>("thetaYEnd","#theta_{y} (end);#theta_{y}", 100, 0., TMath::Pi()));
00240   
00241   hists.Add(fd.make<TH1F>("momStartEnd", "|p_{start}|-|p_{end}|;#Deltap [GeV]",100,0.,15.));
00242   hists.Add(fd.make<TH1F>("momStartEndRel", "(p_{start}-p_{end})/p_{start};#Deltap_{rel}",
00243                           100,.0,1.));
00244   hists.Add(fd.make<TH1F>("phiXzStartEnd", "#phi_{xz,start}-#phi_{xz,end};#Delta#phi_{xz}",
00245                           100,-1.,1.));
00246   hists.Add(fd.make<TH1F>("thetaYStartEnd","#theta_{y,start}-#theta_{y,end};#Delta#theta_{y}",
00247                           100,-1.,1.));
00248   
00249   hists.Add(fd.make<TH2F>("phiXzStartVsEnd",
00250                           "#phi_{xz} start vs end;#phi_{xz}^{end};#phi_{xz}^{start}",
00251                           50, -TMath::Pi(), TMath::Pi(), 50, -TMath::Pi(), TMath::Pi()));
00252   hists.Add(fd.make<TH2F>("thetaYStartVsEnd",
00253                           "#theta_{y} start vs end;#theta_{y}^{end};#theta_{y}^{start}",
00254                           50, 0., TMath::Pi(), 50, 0., TMath::Pi()/2.));
00255   
00256   hists.Add(fd.make<TH2F>("momStartEndRelVsZ",
00257                           "(p_{start}-p_{end})/p_{start} vs z_{start};z [cm];#Deltap_{rel}",
00258                           50, -1600., 1600., 50,.0,.8));
00259   hists.Add(fd.make<TH2F>("phiXzStartEndVsZ", 
00260                           "#phi_{xz,start}-#phi_{xz,end} vs z_{start};z [cm];#Delta#phi_{xz}",
00261                           50, -1600., 1600., 50,-1., 1.));
00262   hists.Add(fd.make<TH2F>("thetaYStartEndVsZ",
00263                           "#theta_{y,start}-#theta_{y,end} vs z_{start};z [cm];#Delta#theta_{y}",
00264                           50, -1600., 1600., 50,-.5,.5));
00265   hists.Add(fd.make<TH2F>("momStartEndRelVsP",
00266                           "(p_{start}-p_{end})/p_{start} vs p_{start};p [GeV];#Deltap_{rel}",
00267                           kNumBins, pBinsLog, 50, .0, .8));
00268   hists.Add(fd.make<TH2F>("phiXzStartEndVsP", 
00269                           "#phi_{xz,start}-#phi_{xz,end} vs |p|_{start};p [GeV];#Delta#phi_{xz}",
00270                           kNumBins, pBinsLog, 100,-1.5, 1.5));
00271   hists.Add(fd.make<TH2F>("thetaYStartEndVsP",
00272                           "#theta_{y,start}-#theta_{y,end} vs |p|_{start};p [GeV];#Delta#theta_{y}",
00273                           kNumBins, pBinsLog, 100,-1.,1.));
00274   
00275   const double maxR = theTargetCylinder->radius() * 1.1;
00276   hists.Add(fd.make<TH1F>("radiusEnd", "end radius;r [cm]", 100, 0., maxR));
00277   double minZ = theTargetPlaneMin->position().z();
00278   minZ -= TMath::Abs(minZ) * 0.1;
00279   double maxZ = theTargetPlaneMax->position().z();
00280   maxZ += TMath::Abs(maxZ) * 0.1;
00281   hists.Add(fd.make<TH1F>("zEnd", "end z;z [cm]", 100, minZ, maxZ));
00282   hists.Add(fd.make<TH1F>("zDiff", "z_{start}-z_{end};#Deltaz [cm]", 100, -1000., 1000.));
00283   hists.Add(fd.make<TH2F>("xyPlaneEnd", "end xy;x [cm];y [cm]", 100, -maxR, maxR, 100,-maxR,maxR));
00284   
00285   hists.Add(fd.make<TH2F>("rzPlaneEnd", "end rz (y<0 #Rightarrow r_{#pm}=-r);z [cm];r_{#pm} [cm]",
00286                           50, minZ, maxZ, 50, -maxR, maxR));
00287   hists.Add(fd.make<TH2F>("thetaVsZend", "#theta vs z (end);z [cm];#theta",
00288                           50, minZ, maxZ, 50, 0., TMath::Pi()));
00289 }
00290 
00291 // ------------ method called once each job just after ending the event loop  ------------
00292 void CosmicGenFilterHelix::endJob()
00293 {
00294   const char *border = "////////////////////////////////////////////////////////";
00295   const char *line = "\n// ";
00296   edm::LogInfo("Filter") << "@SUB=CosmicGenFilterHelix::endJob"
00297                          << border << line
00298                          << theNumPass << " events out of " << theNumTotal
00299                          << ", i.e. " << theNumPass*100./theNumTotal << "%, "
00300                          << "reached target cylinder," 
00301                          << line << "defined by r < " 
00302                          << theTargetCylinder->radius() << " cm and " 
00303                          << theTargetPlaneMin->position().z() << " < z < " 
00304                          << theTargetPlaneMax->position().z() << " cm."
00305                          << line << "Minimal required (transverse) momentum was "
00306                          << TMath::Sqrt(theMinP2) << " (" << TMath::Sqrt(theMinPt2) << ") GeV."
00307                          << "\n" << border;
00308 }
00309 
00311 bool CosmicGenFilterHelix::charge(int id, int &charge) const
00312 {
00313   std::vector<int>::const_iterator iC = theCharges.begin();
00314   for (std::vector<int>::const_iterator i = theIds.begin(), end = theIds.end(); i != end;
00315        ++i,++iC) {
00316     if (*i == id) {
00317       charge = *iC;
00318       return true;
00319     }
00320   }
00321 
00322   return false;
00323 }
00324 
00325 //_________________________________________________________________________________________________
00326 const MagneticField* CosmicGenFilterHelix::getMagneticField(const edm::EventSetup &setup) const
00327 {
00328   edm::ESHandle<MagneticField> fieldHandle;
00329   setup.get<IdealMagneticFieldRecord>().get(fieldHandle); 
00330 
00331   return fieldHandle.product();
00332 }
00333 
00334 //_________________________________________________________________________________________________
00335 const Propagator* CosmicGenFilterHelix::getPropagator(const edm::EventSetup &setup) const
00336 {
00337   edm::ESHandle<Propagator> propHandle;
00338   setup.get<TrackingComponentsRecord>().get(thePropagatorName, propHandle);
00339 
00340   const Propagator *prop = propHandle.product();
00341   if (!dynamic_cast<const SteppingHelixPropagator*>(prop)) {
00342     edm::LogWarning("BadConfig") << "@SUB=CosmicGenFilterHelix::getPropagator"
00343                                  << "Not a SteppingHelixPropagator!";
00344 
00345   }
00346   return prop;
00347 }
00348 
00349 //_________________________________________________________________________________________________
00350 void CosmicGenFilterHelix::monitorStart(const GlobalPoint &vert, const GlobalVector &mom,
00351                                         int charge, TObjArray &hists)
00352 {
00353   const double scalarMom = mom.mag();
00354   const double phi = mom.phi();
00355   const double phiXz = TMath::ATan2(mom.z(), mom.x());
00356   const double theta = mom.theta();
00357   const double thetaY = TMath::ATan2(TMath::Sqrt(mom.x()*mom.x()+mom.z()*mom.z()), -mom.y());
00358 
00359   const double z = vert.z();
00360   const double r = vert.perp();
00361 
00362   static int iMomP = hists.IndexOf(hists.FindObject("momentumP"));
00363   static int iMomM = hists.IndexOf(hists.FindObject("momentumM"));
00364   if (charge > 0) static_cast<TH1*>(hists[iMomP])->Fill(scalarMom);
00365   else            static_cast<TH1*>(hists[iMomM])->Fill(scalarMom);
00366   static int iMom2 = hists.IndexOf(hists.FindObject("momentum2"));
00367   static_cast<TH1*>(hists[iMom2])->Fill(scalarMom);
00368   static int iMomLog = hists.IndexOf(hists.FindObject("momentumLog"));
00369   static_cast<TH1*>(hists[iMomLog])->Fill(scalarMom);
00370   static int iPhi = hists.IndexOf(hists.FindObject("phi"));
00371   static_cast<TH1*>(hists[iPhi])->Fill(phi);
00372   static int iCosPhi = hists.IndexOf(hists.FindObject("cosPhi"));
00373   static_cast<TH1*>(hists[iCosPhi])->Fill(TMath::Cos(phi));
00374   static int iPhiXz = hists.IndexOf(hists.FindObject("phiXz"));
00375   static_cast<TH1*>(hists[iPhiXz])->Fill(phiXz);
00376   static int iTheta = hists.IndexOf(hists.FindObject("theta"));
00377   static_cast<TH1*>(hists[iTheta])->Fill(theta);
00378   static int iThetaY = hists.IndexOf(hists.FindObject("thetaY"));
00379   static_cast<TH1*>(hists[iThetaY])->Fill(thetaY);
00380 
00381   static int iMomVsTheta = hists.IndexOf(hists.FindObject("momVsTheta"));
00382   static_cast<TH2*>(hists[iMomVsTheta])->Fill(theta, scalarMom);
00383   static int iMomVsThetaY = hists.IndexOf(hists.FindObject("momVsThetaY"));
00384   static_cast<TH2*>(hists[iMomVsThetaY])->Fill(thetaY, scalarMom);
00385   static int iMomVsPhi = hists.IndexOf(hists.FindObject("momVsPhi"));
00386   static_cast<TH2*>(hists[iMomVsPhi])->Fill(phi, scalarMom);
00387   static int iMomVsZ = hists.IndexOf(hists.FindObject("momVsZ"));
00388   static_cast<TH2*>(hists[iMomVsZ])->Fill(z, scalarMom);
00389   static int iThetaVsZ = hists.IndexOf(hists.FindObject("thetaVsZ"));
00390   static_cast<TH2*>(hists[iThetaVsZ])->Fill(z, theta);
00391   static int iThetaYvsZ = hists.IndexOf(hists.FindObject("thetaYvsZ"));
00392   static_cast<TH2*>(hists[iThetaYvsZ])->Fill(z, thetaY);
00393   static int iYvsThetaY = hists.IndexOf(hists.FindObject("yVsThetaY"));
00394   static_cast<TH2*>(hists[iYvsThetaY])->Fill(thetaY, vert.y());
00395   static int iYvsThetaYnoR = hists.IndexOf(hists.FindObject("yVsThetaYnoR"));
00396   if (z > -1400. && z < 1400.) {
00397     static_cast<TH2*>(hists[iYvsThetaYnoR])->Fill(thetaY, vert.y());
00398   }
00399 
00400   static int iRadius = hists.IndexOf(hists.FindObject("radius"));
00401   static_cast<TH1*>(hists[iRadius])->Fill(r);
00402   static int iZ = hists.IndexOf(hists.FindObject("z"));
00403   static_cast<TH1*>(hists[iZ])->Fill(z);
00404   static int iXy = hists.IndexOf(hists.FindObject("xyPlane"));
00405   static_cast<TH1*>(hists[iXy])->Fill(vert.x(), vert.y());
00406   static int iRz = hists.IndexOf(hists.FindObject("rzPlane"));
00407   static_cast<TH1*>(hists[iRz])->Fill(z, (vert.y() > 0 ? r : -r));
00408 }
00409 
00410 //_________________________________________________________________________________________________
00411 void CosmicGenFilterHelix::monitorEnd(const GlobalPoint &endVert, const GlobalVector &endMom,
00412                                       const GlobalPoint &vert, const GlobalVector &mom,
00413                                       double path, TObjArray &hists)
00414 {
00415   const double scalarMomStart = mom.mag();
00416   const double phiXzStart = TMath::ATan2(mom.z(), mom.x());
00417   const double thetaYStart = TMath::ATan2(TMath::Sqrt(mom.x()*mom.x()+mom.z()*mom.z()), -mom.y());
00418   const double scalarMomEnd = endMom.mag();
00419   const double ptEnd = endMom.perp();
00420   const double phiXzEnd = TMath::ATan2(endMom.z(), endMom.x());
00421   const double thetaYEnd = TMath::ATan2(TMath::Sqrt(endMom.x()*endMom.x()+endMom.z()*endMom.z()),
00422                                         -endMom.y());
00423   const double thetaEnd = endMom.theta();
00424 
00425   const double diffMomRel = (scalarMomStart-scalarMomEnd)/scalarMomStart;
00426 
00427   const double zEnd = endVert.z();
00428   const double rEnd = endVert.perp();
00429   const double diffZ = zEnd - vert.z();
00430 
00431   static int iPathEnd = hists.IndexOf(hists.FindObject("pathEnd"));
00432   static_cast<TH1*>(hists[iPathEnd])->Fill(path);
00433   static int iMomEnd = hists.IndexOf(hists.FindObject("momEnd"));
00434   static_cast<TH1*>(hists[iMomEnd])->Fill(scalarMomEnd);
00435   static int iMomEndLog = hists.IndexOf(hists.FindObject("momEndLog"));
00436   static_cast<TH1*>(hists[iMomEndLog])->Fill(scalarMomEnd);
00437   static int iPtEnd = hists.IndexOf(hists.FindObject("ptEnd"));
00438   static_cast<TH1*>(hists[iPtEnd])->Fill(ptEnd);
00439   static int iPtEndLog = hists.IndexOf(hists.FindObject("ptEndLog"));
00440   static_cast<TH1*>(hists[iPtEndLog])->Fill(ptEnd);
00441   static int iPhiXzEnd = hists.IndexOf(hists.FindObject("phiXzEnd"));
00442   static_cast<TH1*>(hists[iPhiXzEnd])->Fill(phiXzEnd);
00443   static int iThetaYEnd = hists.IndexOf(hists.FindObject("thetaYEnd"));
00444   static_cast<TH1*>(hists[iThetaYEnd])->Fill(thetaYEnd);
00445 
00446   static int iMomStartEnd = hists.IndexOf(hists.FindObject("momStartEnd"));
00447   static_cast<TH1*>(hists[iMomStartEnd])->Fill(scalarMomStart-scalarMomEnd);
00448   static int iMomStartEndRel = hists.IndexOf(hists.FindObject("momStartEndRel"));
00449   static_cast<TH1*>(hists[iMomStartEndRel])->Fill(diffMomRel);
00450   static int iPhiStartEnd = hists.IndexOf(hists.FindObject("phiXzStartEnd"));
00451   static_cast<TH1*>(hists[iPhiStartEnd])->Fill(phiXzStart-phiXzEnd);
00452   static int iThetaStartEnd = hists.IndexOf(hists.FindObject("thetaYStartEnd"));
00453   static_cast<TH1*>(hists[iThetaStartEnd])->Fill(thetaYStart-thetaYEnd);
00454 
00455   static int iPhiStartVsEnd = hists.IndexOf(hists.FindObject("phiXzStartVsEnd"));
00456   static_cast<TH2*>(hists[iPhiStartVsEnd])->Fill(phiXzEnd, phiXzStart);
00457   static int iThetaStartVsEnd = hists.IndexOf(hists.FindObject("thetaYStartVsEnd"));
00458   static_cast<TH2*>(hists[iThetaStartVsEnd])->Fill(thetaYEnd, thetaYStart);
00459 
00460   static int iMomStartEndRelVsZ = hists.IndexOf(hists.FindObject("momStartEndRelVsZ"));
00461   static_cast<TH2*>(hists[iMomStartEndRelVsZ])->Fill(vert.z(), diffMomRel);
00462   static int iPhiStartEndVsZ = hists.IndexOf(hists.FindObject("phiXzStartEndVsZ"));
00463   static_cast<TH2*>(hists[iPhiStartEndVsZ])->Fill(vert.z(), phiXzStart-phiXzEnd);
00464   static int iThetaStartEndVsZ = hists.IndexOf(hists.FindObject("thetaYStartEndVsZ"));
00465   static_cast<TH2*>(hists[iThetaStartEndVsZ])->Fill(vert.z(), thetaYStart-thetaYEnd);
00466   static int iMomStartEndRelVsP = hists.IndexOf(hists.FindObject("momStartEndRelVsP"));
00467   static_cast<TH2*>(hists[iMomStartEndRelVsP])->Fill(scalarMomStart, diffMomRel);
00468   static int iPhiStartEndVsP = hists.IndexOf(hists.FindObject("phiXzStartEndVsP"));
00469   static_cast<TH2*>(hists[iPhiStartEndVsP])->Fill(scalarMomStart, phiXzStart-phiXzEnd);
00470   static int iThetaStartEndVsP = hists.IndexOf(hists.FindObject("thetaYStartEndVsP"));
00471   static_cast<TH2*>(hists[iThetaStartEndVsP])->Fill(scalarMomStart, thetaYStart-thetaYEnd);
00472 
00473   static int iRadiusEnd = hists.IndexOf(hists.FindObject("radiusEnd"));
00474   static_cast<TH1*>(hists[iRadiusEnd])->Fill(rEnd);
00475   static int iZend = hists.IndexOf(hists.FindObject("zEnd"));
00476   static_cast<TH1*>(hists[iZend])->Fill(zEnd);
00477   static int iZdiff = hists.IndexOf(hists.FindObject("zDiff"));
00478   static_cast<TH1*>(hists[iZdiff])->Fill(diffZ);
00479   static int iXyPlaneEnd = hists.IndexOf(hists.FindObject("xyPlaneEnd"));
00480   static_cast<TH1*>(hists[iXyPlaneEnd])->Fill(endVert.x(), endVert.y());
00481   static int iRzPlaneEnd = hists.IndexOf(hists.FindObject("rzPlaneEnd"));
00482   static_cast<TH1*>(hists[iRzPlaneEnd])->Fill(zEnd, (endVert.y() > 0 ? rEnd : -rEnd));
00483   static int iThetaVsZend = hists.IndexOf(hists.FindObject("thetaVsZend"));
00484   static_cast<TH2*>(hists[iThetaVsZend])->Fill(zEnd, thetaEnd);
00485 }
00486 
00487 //_________________________________________________________________________________________________
00488 bool CosmicGenFilterHelix::equidistLogBins(double* bins, int nBins, 
00489                                            double first, double last) const
00490 {
00491   // Filling 'bins' with borders of 'nBins' bins between 'first' and 'last'
00492   // that are equidistant when viewed in log scale,
00493   // so 'bins' must have length nBins+1;
00494   // If 'first', 'last' or 'nBins' are not positive, failure is reported.
00495 
00496   if (nBins < 1 || first <= 0. || last <= 0.) return false;
00497 
00498   bins[0] = first;
00499   bins[nBins] = last;
00500   const double firstLog = TMath::Log10(bins[0]);
00501   const double lastLog  = TMath::Log10(bins[nBins]);
00502   for (int i = 1; i < nBins; ++i) {
00503     bins[i] = TMath::Power(10., firstLog + i*(lastLog-firstLog)/(nBins));
00504   }
00505 
00506   return true;
00507 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/GeneratorInterface/GenFilters/src/CosmicGenFilterHelix.cc
