Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:35 

0001 
0002 #include <memory>
0003 #include <string>
0004 #include <iostream>
0005 #include <fstream>
0006 #include <TMath.h>
0007 
0008 #include "FWCore/PluginManager/interface/ModuleDef.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 
0011 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0013 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0014 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0015 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0016 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0017 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0018 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0019 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0020 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0021 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0022 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0023 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0024 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0025 
0026 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0027 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0028 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0029 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0030 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0031 #include "SiPixelLorentzAngle.h"
0032 
0033 using namespace std;
0034 using namespace edm;
0035 using namespace reco;
0036 using namespace analyzer;
0037 
0038 SiPixelLorentzAngle::SiPixelLorentzAngle(edm::ParameterSet const& conf)
0039     : filename_(conf.getParameter<std::string>("fileName")),
0040       filenameFit_(conf.getParameter<std::string>("fileNameFit")),
0041       ptmin_(conf.getParameter<double>("ptMin")),
0042       simData_(conf.getParameter<bool>("simData")),
0043       normChi2Max_(conf.getParameter<double>("normChi2Max")),
0044       clustSizeYMin_(conf.getParameter<int>("clustSizeYMin")),
0045       residualMax_(conf.getParameter<double>("residualMax")),
0046       clustChargeMax_(conf.getParameter<double>("clustChargeMax")),
0047       hist_depth_(conf.getParameter<int>("binsDepth")),
0048       hist_drift_(conf.getParameter<int>("binsDrift")),
0049       trackerHitAssociatorConfig_(consumesCollector()) {
0050   //    anglefinder_=new  TrackLocalAngle(conf);
0051   hist_x_ = 50;
0052   hist_y_ = 100;
0053   min_x_ = -500.;
0054   max_x_ = 500.;
0055   min_y_ = -1500.;
0056   max_y_ = 500.;
0057   width_ = 0.0285;
0058   min_depth_ = -100.;
0059   max_depth_ = 400.;
0060   min_drift_ = -1000.;  //-200.;(conf.getParameter<double>("residualMax"))
0061   max_drift_ = 1000.;   //400.;
0062 
0063   t_trajTrack = consumes<TrajTrackAssociationCollection>(conf.getParameter<edm::InputTag>("src"));
0064   trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0065   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0066 }
0067 
0068 // Virtual destructor needed.
0069 SiPixelLorentzAngle::~SiPixelLorentzAngle() = default;
0070 
0071 void SiPixelLorentzAngle::beginJob() {
0072   //    cout << "started SiPixelLorentzAngle" << endl;
0073   hFile_ = new TFile(filename_.c_str(), "RECREATE");
0074   int bufsize = 64000;
0075   // create tree structure
0076   SiPixelLorentzAngleTree_ = new TTree("SiPixelLorentzAngleTree_", "SiPixel LorentzAngle tree", bufsize);
0077   SiPixelLorentzAngleTree_->Branch("run", &run_, "run/I", bufsize);
0078   SiPixelLorentzAngleTree_->Branch("event", &event_, "event/I", bufsize);
0079   SiPixelLorentzAngleTree_->Branch("module", &module_, "module/I", bufsize);
0080   SiPixelLorentzAngleTree_->Branch("ladder", &ladder_, "ladder/I", bufsize);
0081   SiPixelLorentzAngleTree_->Branch("layer", &layer_, "layer/I", bufsize);
0082   SiPixelLorentzAngleTree_->Branch("isflipped", &isflipped_, "isflipped/I", bufsize);
0083   SiPixelLorentzAngleTree_->Branch("pt", &pt_, "pt/F", bufsize);
0084   SiPixelLorentzAngleTree_->Branch("eta", &eta_, "eta/F", bufsize);
0085   SiPixelLorentzAngleTree_->Branch("phi", &phi_, "phi/F", bufsize);
0086   SiPixelLorentzAngleTree_->Branch("chi2", &chi2_, "chi2/D", bufsize);
0087   SiPixelLorentzAngleTree_->Branch("ndof", &ndof_, "ndof/D", bufsize);
0088   SiPixelLorentzAngleTree_->Branch("trackhit", &trackhit_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
0089   SiPixelLorentzAngleTree_->Branch("simhit", &simhit_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
0090   SiPixelLorentzAngleTree_->Branch("npix", &pixinfo_.npix, "npix/I", bufsize);
0091   SiPixelLorentzAngleTree_->Branch("rowpix", pixinfo_.row, "row[npix]/F", bufsize);
0092   SiPixelLorentzAngleTree_->Branch("colpix", pixinfo_.col, "col[npix]/F", bufsize);
0093   SiPixelLorentzAngleTree_->Branch("adc", pixinfo_.adc, "adc[npix]/F", bufsize);
0094   SiPixelLorentzAngleTree_->Branch("xpix", pixinfo_.x, "x[npix]/F", bufsize);
0095   SiPixelLorentzAngleTree_->Branch("ypix", pixinfo_.y, "y[npix]/F", bufsize);
0096   SiPixelLorentzAngleTree_->Branch(
0097       "clust",
0098       &clust_,
0099       "x/F:y/F:charge/F:size_x/I:size_y/I:maxPixelCol/I:maxPixelRow:minPixelCol/I:minPixelRow/I",
0100       bufsize);
0101   SiPixelLorentzAngleTree_->Branch("rechit", &rechit_, "x/F:y/F", bufsize);
0102 
0103   SiPixelLorentzAngleTreeForward_ =
0104       new TTree("SiPixelLorentzAngleTreeForward_", "SiPixel LorentzAngle tree forward", bufsize);
0105   SiPixelLorentzAngleTreeForward_->Branch("run", &run_, "run/I", bufsize);
0106   SiPixelLorentzAngleTreeForward_->Branch("event", &event_, "event/I", bufsize);
0107   SiPixelLorentzAngleTreeForward_->Branch("side", &sideF_, "side/I", bufsize);
0108   SiPixelLorentzAngleTreeForward_->Branch("disk", &diskF_, "disk/I", bufsize);
0109   SiPixelLorentzAngleTreeForward_->Branch("blade", &bladeF_, "blade/I", bufsize);
0110   SiPixelLorentzAngleTreeForward_->Branch("panel", &panelF_, "panel/I", bufsize);
0111   SiPixelLorentzAngleTreeForward_->Branch("module", &moduleF_, "module/I", bufsize);
0112   SiPixelLorentzAngleTreeForward_->Branch("pt", &pt_, "pt/F", bufsize);
0113   SiPixelLorentzAngleTreeForward_->Branch("eta", &eta_, "eta/F", bufsize);
0114   SiPixelLorentzAngleTreeForward_->Branch("phi", &phi_, "phi/F", bufsize);
0115   SiPixelLorentzAngleTreeForward_->Branch("chi2", &chi2_, "chi2/D", bufsize);
0116   SiPixelLorentzAngleTreeForward_->Branch("ndof", &ndof_, "ndof/D", bufsize);
0117   SiPixelLorentzAngleTreeForward_->Branch("trackhit", &trackhitF_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
0118   SiPixelLorentzAngleTreeForward_->Branch("simhit", &simhitF_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
0119   SiPixelLorentzAngleTreeForward_->Branch("npix", &pixinfoF_.npix, "npix/I", bufsize);
0120   SiPixelLorentzAngleTreeForward_->Branch("rowpix", pixinfoF_.row, "row[npix]/F", bufsize);
0121   SiPixelLorentzAngleTreeForward_->Branch("colpix", pixinfoF_.col, "col[npix]/F", bufsize);
0122   SiPixelLorentzAngleTreeForward_->Branch("adc", pixinfoF_.adc, "adc[npix]/F", bufsize);
0123   SiPixelLorentzAngleTreeForward_->Branch("xpix", pixinfoF_.x, "x[npix]/F", bufsize);
0124   SiPixelLorentzAngleTreeForward_->Branch("ypix", pixinfoF_.y, "y[npix]/F", bufsize);
0125   SiPixelLorentzAngleTreeForward_->Branch(
0126       "clust",
0127       &clustF_,
0128       "x/F:y/F:charge/F:size_x/I:size_y/I:maxPixelCol/I:maxPixelRow:minPixelCol/I:minPixelRow/I",
0129       bufsize);
0130   SiPixelLorentzAngleTreeForward_->Branch("rechit", &rechitF_, "x/F:y/F", bufsize);
0131 
0132   //book histograms
0133   char name[128];
0134   for (int i_module = 1; i_module <= 8; i_module++) {
0135     for (int i_layer = 1; i_layer <= 3; i_layer++) {
0136       sprintf(name, "h_drift_depth_adc_layer%i_module%i", i_layer, i_module);
0137       _h_drift_depth_adc_[i_module + (i_layer - 1) * 8] =
0138           new TH2F(name, name, hist_drift_, min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0139       sprintf(name, "h_drift_depth_adc2_layer%i_module%i", i_layer, i_module);
0140       _h_drift_depth_adc2_[i_module + (i_layer - 1) * 8] =
0141           new TH2F(name, name, hist_drift_, min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0142       sprintf(name, "h_drift_depth_noadc_layer%i_module%i", i_layer, i_module);
0143       _h_drift_depth_noadc_[i_module + (i_layer - 1) * 8] =
0144           new TH2F(name, name, hist_drift_, min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0145       sprintf(name, "h_drift_depth_layer%i_module%i", i_layer, i_module);
0146       _h_drift_depth_[i_module + (i_layer - 1) * 8] =
0147           new TH2F(name, name, hist_drift_, min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0148       sprintf(name, "h_mean_layer%i_module%i", i_layer, i_module);
0149       _h_mean_[i_module + (i_layer - 1) * 8] = new TH1F(name, name, hist_depth_, min_depth_, max_depth_);
0150     }
0151   }
0152 
0153   // just for some expaining plots
0154   h_cluster_shape_adc_ = new TH2F(
0155       "h_cluster_shape_adc", "cluster shape with adc weight", hist_x_, min_x_, max_x_, hist_y_, min_y_, max_y_);
0156   h_cluster_shape_noadc_ = new TH2F(
0157       "h_cluster_shape_noadc", "cluster shape without adc weight", hist_x_, min_x_, max_x_, hist_y_, min_y_, max_y_);
0158   h_cluster_shape_ = new TH2F("h_cluster_shape", "cluster shape", hist_x_, min_x_, max_x_, hist_y_, min_y_, max_y_);
0159   h_cluster_shape_adc_rot_ = new TH2F(
0160       "h_cluster_shape_adc_rot", "cluster shape with adc weight", hist_x_, min_x_, max_x_, hist_y_, -max_y_, -min_y_);
0161   h_cluster_shape_noadc_rot_ = new TH2F("h_cluster_shape_noadc_rot",
0162                                         "cluster shape without adc weight",
0163                                         hist_x_,
0164                                         min_x_,
0165                                         max_x_,
0166                                         hist_y_,
0167                                         -max_y_,
0168                                         -min_y_);
0169   h_cluster_shape_rot_ =
0170       new TH2F("h_cluster_shape_rot", "cluster shape", hist_x_, min_x_, max_x_, hist_y_, -max_y_, -min_y_);
0171   h_tracks_ = new TH1F("h_tracks", "h_tracks", 2, 0., 2.);
0172   event_counter_ = 0;
0173   trackEventsCounter_ = 0;
0174   //    trackcounter_ = 0;
0175   hitCounter_ = 0;
0176   usedHitCounter_ = 0;
0177   pixelTracksCounter_ = 0;
0178 }
0179 
0180 // Functions that gets called by framework every event
0181 void SiPixelLorentzAngle::analyze(const edm::Event& e, const edm::EventSetup& es) {
0182   //Retrieve tracker topology from geometry
0183   edm::ESHandle<TrackerTopology> tTopoHandle = es.getHandle(trackerTopoToken_);
0184   const TrackerTopology* const tTopo = tTopoHandle.product();
0185 
0186   event_counter_++;
0187   //    if(event_counter_ % 500 == 0) cout << "event number " << event_counter_ << endl;
0188   cout << "event number " << event_counter_ << endl;
0189 
0190   edm::ESHandle<TrackerGeometry> estracker = es.getHandle(trackerGeomToken_);
0191   const TrackerGeometry* tracker = &(*estracker);
0192 
0193   std::unique_ptr<TrackerHitAssociator> associate;
0194   if (simData_)
0195     associate = std::make_unique<TrackerHitAssociator>(e, trackerHitAssociatorConfig_);
0196   // restet values
0197   module_ = -1;
0198   layer_ = -1;
0199   ladder_ = -1;
0200   isflipped_ = -1;
0201   pt_ = -999;
0202   eta_ = 999;
0203   phi_ = 999;
0204   pixinfo_.npix = 0;
0205 
0206   run_ = e.id().run();
0207   event_ = e.id().event();
0208 
0209   // get the association map between tracks and trajectories
0210   edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
0211   e.getByToken(t_trajTrack, trajTrackCollectionHandle);
0212   if (!trajTrackCollectionHandle->empty()) {
0213     trackEventsCounter_++;
0214     for (TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin();
0215          it != trajTrackCollectionHandle->end();
0216          ++it) {
0217       const Track& track = *it->val;
0218       const Trajectory& traj = *it->key;
0219 
0220       // get the trajectory measurements
0221       std::vector<TrajectoryMeasurement> tmColl = traj.measurements();
0222       //            TrajectoryStateOnSurface tsos = tsoscomb( itTraj->forwardPredictedState(), itTraj->backwardPredictedState() );
0223       pt_ = track.pt();
0224       eta_ = track.eta();
0225       phi_ = track.phi();
0226       chi2_ = traj.chiSquared();
0227       ndof_ = traj.ndof();
0228       if (pt_ < ptmin_)
0229         continue;
0230       // iterate over trajectory measurements
0231       std::vector<PSimHit> matched;
0232       h_tracks_->Fill(0);
0233       bool pixeltrack = false;
0234       for (std::vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj != tmColl.end();
0235            itTraj++) {
0236         if (!itTraj->updatedState().isValid())
0237           continue;
0238         TransientTrackingRecHit::ConstRecHitPointer recHit = itTraj->recHit();
0239         if (!recHit->isValid() || recHit->geographicalId().det() != DetId::Tracker)
0240           continue;
0241         unsigned int subDetID = (recHit->geographicalId().subdetId());
0242         if (subDetID == PixelSubdetector::PixelBarrel || subDetID == PixelSubdetector::PixelEndcap) {
0243           if (!pixeltrack) {
0244             h_tracks_->Fill(1);
0245             pixelTracksCounter_++;
0246           }
0247           pixeltrack = true;
0248         }
0249         if (subDetID == PixelSubdetector::PixelBarrel) {
0250           hitCounter_++;
0251 
0252           DetId detIdObj = recHit->geographicalId();
0253           const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(tracker->idToDet(detIdObj));
0254           if (!theGeomDet)
0255             continue;
0256 
0257           const PixelTopology* topol = &(theGeomDet->specificTopology());
0258 
0259           if (!topol)
0260             continue;
0261 
0262           layer_ = tTopo->pxbLayer(detIdObj);
0263           ladder_ = tTopo->pxbLadder(detIdObj);
0264           module_ = tTopo->pxbModule(detIdObj);
0265           float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 0.)).perp();
0266           float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 1.)).perp();
0267           if (tmp2 < tmp1)
0268             isflipped_ = 1;
0269           else
0270             isflipped_ = 0;
0271           const SiPixelRecHit* recHitPix = dynamic_cast<const SiPixelRecHit*>((*recHit).hit());
0272           if (!recHitPix)
0273             continue;
0274           rechit_.x = recHitPix->localPosition().x();
0275           rechit_.y = recHitPix->localPosition().y();
0276           SiPixelRecHit::ClusterRef const& cluster = recHitPix->cluster();
0277 
0278           // fill entries in clust_
0279           clust_.x = (cluster)->x();
0280           clust_.y = (cluster)->y();
0281           clust_.charge = (cluster->charge()) / 1000.;
0282           clust_.size_x = cluster->sizeX();
0283           clust_.size_y = cluster->sizeY();
0284           clust_.maxPixelCol = cluster->maxPixelCol();
0285           clust_.maxPixelRow = cluster->maxPixelRow();
0286           clust_.minPixelCol = cluster->minPixelCol();
0287           clust_.minPixelRow = cluster->minPixelRow();
0288           // fill entries in pixinfo_:
0289           fillPix(*cluster, topol, pixinfo_);
0290           // fill the trackhit info
0291           TrajectoryStateOnSurface tsos = itTraj->updatedState();
0292           if (!tsos.isValid()) {
0293             cout << "tsos not valid" << endl;
0294             continue;
0295           }
0296           LocalVector trackdirection = tsos.localDirection();
0297           LocalPoint trackposition = tsos.localPosition();
0298 
0299           if (trackdirection.z() == 0)
0300             continue;
0301           // the local position and direction
0302           trackhit_.alpha = atan2(trackdirection.z(), trackdirection.x());
0303           trackhit_.beta = atan2(trackdirection.z(), trackdirection.y());
0304           trackhit_.gamma = atan2(trackdirection.x(), trackdirection.y());
0305           trackhit_.x = trackposition.x();
0306           trackhit_.y = trackposition.y();
0307 
0308           // fill entries in simhit_:
0309           if (simData_) {
0310             matched.clear();
0311             matched = associate->associateHit((*recHitPix));
0312             float dr_start = 9999.;
0313             for (std::vector<PSimHit>::iterator isim = matched.begin(); isim != matched.end(); ++isim) {
0314               DetId simdetIdObj((*isim).detUnitId());
0315               if (simdetIdObj == detIdObj) {
0316                 float sim_x1 = (*isim).entryPoint().x();  // width (row index, in col direction)
0317                 float sim_y1 = (*isim).entryPoint().y();  // length (col index, in row direction)
0318                 float sim_x2 = (*isim).exitPoint().x();
0319                 float sim_y2 = (*isim).exitPoint().y();
0320                 float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0321                 float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0322                 float sim_px = (*isim).momentumAtEntry().x();
0323                 float sim_py = (*isim).momentumAtEntry().y();
0324                 float sim_pz = (*isim).momentumAtEntry().z();
0325 
0326                 float dr = (sim_xpos - (recHitPix->localPosition().x())) * (sim_xpos - recHitPix->localPosition().x()) +
0327                            (sim_ypos - recHitPix->localPosition().y()) * (sim_ypos - recHitPix->localPosition().y());
0328                 if (dr < dr_start) {
0329                   simhit_.x = sim_xpos;
0330                   simhit_.y = sim_ypos;
0331                   simhit_.alpha = atan2(sim_pz, sim_px);
0332                   simhit_.beta = atan2(sim_pz, sim_py);
0333                   simhit_.gamma = atan2(sim_px, sim_py);
0334                   dr_start = dr;
0335                 }
0336               }
0337             }  // end of filling simhit_
0338           }
0339           // is one pixel in cluster a large pixel ? (hit will be excluded)
0340           bool large_pix = false;
0341           for (int j = 0; j < pixinfo_.npix; j++) {
0342             int colpos = static_cast<int>(pixinfo_.col[j]);
0343             if (pixinfo_.row[j] == 0 || pixinfo_.row[j] == 79 || pixinfo_.row[j] == 80 || pixinfo_.row[j] == 159 ||
0344                 colpos % 52 == 0 || colpos % 52 == 51) {
0345               large_pix = true;
0346             }
0347           }
0348 
0349           double residual = TMath::Sqrt((trackhit_.x - rechit_.x) * (trackhit_.x - rechit_.x) +
0350                                         (trackhit_.y - rechit_.y) * (trackhit_.y - rechit_.y));
0351 
0352           SiPixelLorentzAngleTree_->Fill();
0353           if (!large_pix && (chi2_ / ndof_) < normChi2Max_ && cluster->sizeY() >= clustSizeYMin_ &&
0354               residual < residualMax_ && (cluster->charge() < clustChargeMax_)) {
0355             usedHitCounter_++;
0356             // iterate over pixels in hit
0357             for (int j = 0; j < pixinfo_.npix; j++) {
0358               // use trackhits
0359               float dx = (pixinfo_.x[j] - (trackhit_.x - width_ / 2. / TMath::Tan(trackhit_.alpha))) * 10000.;
0360               float dy = (pixinfo_.y[j] - (trackhit_.y - width_ / 2. / TMath::Tan(trackhit_.beta))) * 10000.;
0361               float depth = dy * tan(trackhit_.beta);
0362               float drift = dx - dy * tan(trackhit_.gamma);
0363               _h_drift_depth_adc_[module_ + (layer_ - 1) * 8]->Fill(drift, depth, pixinfo_.adc[j]);
0364               _h_drift_depth_adc2_[module_ + (layer_ - 1) * 8]->Fill(drift, depth, pixinfo_.adc[j] * pixinfo_.adc[j]);
0365               _h_drift_depth_noadc_[module_ + (layer_ - 1) * 8]->Fill(drift, depth);
0366               if (layer_ == 3 && module_ == 1 && isflipped_) {
0367                 float dx_rot = dx * TMath::Cos(trackhit_.gamma) + dy * TMath::Sin(trackhit_.gamma);
0368                 float dy_rot = dy * TMath::Cos(trackhit_.gamma) - dx * TMath::Sin(trackhit_.gamma);
0369                 h_cluster_shape_adc_->Fill(dx, dy, pixinfo_.adc[j]);
0370                 h_cluster_shape_noadc_->Fill(dx, dy);
0371                 h_cluster_shape_adc_rot_->Fill(dx_rot, dy_rot, pixinfo_.adc[j]);
0372                 h_cluster_shape_noadc_rot_->Fill(dx_rot, dy_rot);
0373               }
0374             }
0375           }
0376         } else if (subDetID == PixelSubdetector::PixelEndcap) {
0377           DetId detIdObj = recHit->geographicalId();
0378           const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(tracker->idToDet(detIdObj));
0379           if (!theGeomDet)
0380             continue;
0381 
0382           const PixelTopology* topol = &(theGeomDet->specificTopology());
0383 
0384           if (!topol)
0385             continue;
0386 
0387           sideF_ = tTopo->pxfSide(detIdObj);
0388           diskF_ = tTopo->pxfDisk(detIdObj);
0389           bladeF_ = tTopo->pxfBlade(detIdObj);
0390           panelF_ = tTopo->pxfPanel(detIdObj);
0391           moduleF_ = tTopo->pxfModule(detIdObj);
0392           //float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
0393           //float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
0394           //if ( tmp2<tmp1 ) isflipped_ = 1;
0395           //else isflipped_ = 0;
0396           const SiPixelRecHit* recHitPix = dynamic_cast<const SiPixelRecHit*>((*recHit).hit());
0397           if (!recHitPix)
0398             continue;
0399           rechitF_.x = recHitPix->localPosition().x();
0400           rechitF_.y = recHitPix->localPosition().y();
0401           SiPixelRecHit::ClusterRef const& cluster = recHitPix->cluster();
0402 
0403           // fill entries in clust_
0404           clustF_.x = (cluster)->x();
0405           clustF_.y = (cluster)->y();
0406           clustF_.charge = (cluster->charge()) / 1000.;
0407           clustF_.size_x = cluster->sizeX();
0408           clustF_.size_y = cluster->sizeY();
0409           clustF_.maxPixelCol = cluster->maxPixelCol();
0410           clustF_.maxPixelRow = cluster->maxPixelRow();
0411           clustF_.minPixelCol = cluster->minPixelCol();
0412           clustF_.minPixelRow = cluster->minPixelRow();
0413           // fill entries in pixinfo_:
0414           fillPix(*cluster, topol, pixinfoF_);
0415           // fill the trackhit info
0416           TrajectoryStateOnSurface tsos = itTraj->updatedState();
0417           if (!tsos.isValid()) {
0418             cout << "tsos not valid" << endl;
0419             continue;
0420           }
0421           LocalVector trackdirection = tsos.localDirection();
0422           LocalPoint trackposition = tsos.localPosition();
0423 
0424           if (trackdirection.z() == 0)
0425             continue;
0426           // the local position and direction
0427           trackhitF_.alpha = atan2(trackdirection.z(), trackdirection.x());
0428           trackhitF_.beta = atan2(trackdirection.z(), trackdirection.y());
0429           trackhitF_.gamma = atan2(trackdirection.x(), trackdirection.y());
0430           trackhitF_.x = trackposition.x();
0431           trackhitF_.y = trackposition.y();
0432 
0433           // fill entries in simhit_:
0434           if (simData_) {
0435             matched.clear();
0436             matched = associate->associateHit((*recHitPix));
0437             float dr_start = 9999.;
0438             for (std::vector<PSimHit>::iterator isim = matched.begin(); isim != matched.end(); ++isim) {
0439               DetId simdetIdObj((*isim).detUnitId());
0440               if (simdetIdObj == detIdObj) {
0441                 float sim_x1 = (*isim).entryPoint().x();  // width (row index, in col direction)
0442                 float sim_y1 = (*isim).entryPoint().y();  // length (col index, in row direction)
0443                 float sim_x2 = (*isim).exitPoint().x();
0444                 float sim_y2 = (*isim).exitPoint().y();
0445                 float sim_xpos = 0.5 * (sim_x1 + sim_x2);
0446                 float sim_ypos = 0.5 * (sim_y1 + sim_y2);
0447                 float sim_px = (*isim).momentumAtEntry().x();
0448                 float sim_py = (*isim).momentumAtEntry().y();
0449                 float sim_pz = (*isim).momentumAtEntry().z();
0450 
0451                 float dr = (sim_xpos - (recHitPix->localPosition().x())) * (sim_xpos - recHitPix->localPosition().x()) +
0452                            (sim_ypos - recHitPix->localPosition().y()) * (sim_ypos - recHitPix->localPosition().y());
0453                 if (dr < dr_start) {
0454                   simhitF_.x = sim_xpos;
0455                   simhitF_.y = sim_ypos;
0456                   simhitF_.alpha = atan2(sim_pz, sim_px);
0457                   simhitF_.beta = atan2(sim_pz, sim_py);
0458                   simhitF_.gamma = atan2(sim_px, sim_py);
0459                   dr_start = dr;
0460                 }
0461               }
0462             }  // end of filling simhit_
0463           }
0464           SiPixelLorentzAngleTreeForward_->Fill();
0465         }
0466       }  //end iteration over trajectory measurements
0467     }    //end iteration over trajectories
0468   }
0469 }
0470 
0471 void SiPixelLorentzAngle::endJob() {
0472   // produce histograms with the average adc counts
0473   for (int i_ring = 1; i_ring <= 24; i_ring++) {
0474     _h_drift_depth_[i_ring]->Divide(_h_drift_depth_adc_[i_ring], _h_drift_depth_noadc_[i_ring]);
0475   }
0476 
0477   h_drift_depth_adc_slice_ =
0478       new TH1F("h_drift_depth_adc_slice", "slice of adc histogram", hist_drift_, min_drift_, max_drift_);
0479 
0480   TF1* f1 = new TF1("f1", "[0] + [1]*x", 50., 235.);
0481   f1->SetParName(0, "p0");
0482   f1->SetParName(1, "p1");
0483   f1->SetParameter(0, 0);
0484   f1->SetParameter(1, 0.4);
0485   ofstream fLorentzFit(filenameFit_.c_str(), ios::trunc);
0486   cout.precision(4);
0487   fLorentzFit << "module"
0488               << "\t"
0489               << "layer"
0490               << "\t"
0491               << "offset"
0492               << "\t"
0493               << "error"
0494               << "\t"
0495               << "slope"
0496               << "\t"
0497               << "error"
0498               << "\t"
0499                  "rel.err"
0500               << "\t"
0501                  "pull"
0502               << "\t"
0503               << "chi2"
0504               << "\t"
0505               << "prob" << endl;
0506   //loop over modlues and layers to fit the lorentz angle
0507   for (int i_layer = 1; i_layer <= 3; i_layer++) {
0508     for (int i_module = 1; i_module <= 8; i_module++) {
0509       //loop over bins in depth (z-local-coordinate) (in order to fit slices)
0510       for (int i = 1; i <= hist_depth_; i++) {
0511         findMean(i, (i_module + (i_layer - 1) * 8));
0512       }  // end loop over bins in depth
0513       _h_mean_[i_module + (i_layer - 1) * 8]->Fit(f1, "ERQ");
0514       double p0 = f1->GetParameter(0);
0515       double e0 = f1->GetParError(0);
0516       double p1 = f1->GetParameter(1);
0517       double e1 = f1->GetParError(1);
0518       double chi2 = f1->GetChisquare();
0519       double prob = f1->GetProb();
0520       fLorentzFit << setprecision(4) << i_module << "\t" << i_layer << "\t" << p0 << "\t" << e0 << "\t" << p1
0521                   << setprecision(3) << "\t" << e1 << "\t" << e1 / p1 * 100. << "\t" << (p1 - 0.424) / e1 << "\t"
0522                   << chi2 << "\t" << prob << endl;
0523     }
0524   }  // end loop over modules and layers
0525   fLorentzFit.close();
0526   hFile_->cd();
0527   for (int i_module = 1; i_module <= 8; i_module++) {
0528     for (int i_layer = 1; i_layer <= 3; i_layer++) {
0529       _h_drift_depth_adc_[i_module + (i_layer - 1) * 8]->Write();
0530       _h_drift_depth_adc2_[i_module + (i_layer - 1) * 8]->Write();
0531       _h_drift_depth_noadc_[i_module + (i_layer - 1) * 8]->Write();
0532       _h_drift_depth_[i_module + (i_layer - 1) * 8]->Write();
0533       _h_mean_[i_module + (i_layer - 1) * 8]->Write();
0534     }
0535   }
0536   h_cluster_shape_adc_->Write();
0537   h_cluster_shape_noadc_->Write();
0538   h_cluster_shape_adc_rot_->Write();
0539   h_cluster_shape_noadc_rot_->Write();
0540   h_tracks_->Write();
0541 
0542   hFile_->Write();
0543   hFile_->Close();
0544   cout << "events: " << event_counter_ << endl;
0545   cout << "events with tracks: " << trackEventsCounter_ << endl;
0546   cout << "events with pixeltracks: " << pixelTracksCounter_ << endl;
0547   cout << "hits in the pixel: " << hitCounter_ << endl;
0548   cout << "number of used Hits: " << usedHitCounter_ << endl;
0549 }
0550 
0551 inline void SiPixelLorentzAngle::fillPix(const SiPixelCluster& LocPix, const PixelTopology* topol, Pixinfo& pixinfo)
0552 
0553 {
0554   const std::vector<SiPixelCluster::Pixel>& pixvector = LocPix.pixels();
0555   pixinfo.npix = 0;
0556   for (std::vector<SiPixelCluster::Pixel>::const_iterator itPix = pixvector.begin(); itPix != pixvector.end();
0557        itPix++) {
0558     //  for(pixinfo.npix = 0; pixinfo.npix < static_cast<int>(pixvector.size()); ++pixinfo.npix) {
0559     pixinfo.row[pixinfo.npix] = itPix->x;
0560     pixinfo.col[pixinfo.npix] = itPix->y;
0561     pixinfo.adc[pixinfo.npix] = itPix->adc;
0562     LocalPoint lp = topol->localPosition(MeasurementPoint(itPix->x + 0.5, itPix->y + 0.5));
0563     pixinfo.x[pixinfo.npix] = lp.x();
0564     pixinfo.y[pixinfo.npix] = lp.y();
0565     pixinfo.npix++;
0566   }
0567 }
0568 
0569 void SiPixelLorentzAngle::findMean(int i, int i_ring) {
0570   double nentries = 0;
0571 
0572   h_drift_depth_adc_slice_->Reset("ICE");
0573 
0574   // determine sigma and sigma^2 of the adc counts and average adc counts
0575   //loop over bins in drift width
0576   for (int j = 1; j <= hist_drift_; j++) {
0577     if (_h_drift_depth_noadc_[i_ring]->GetBinContent(j, i) >= 1) {
0578       double adc_error2 =
0579           (_h_drift_depth_adc2_[i_ring]->GetBinContent(j, i) - _h_drift_depth_adc_[i_ring]->GetBinContent(j, i) *
0580                                                                    _h_drift_depth_adc_[i_ring]->GetBinContent(j, i) /
0581                                                                    _h_drift_depth_noadc_[i_ring]->GetBinContent(j, i)) /
0582           _h_drift_depth_noadc_[i_ring]->GetBinContent(j, i);
0583       _h_drift_depth_adc_[i_ring]->SetBinError(j, i, sqrt(adc_error2));
0584       double error2 = adc_error2 / (_h_drift_depth_noadc_[i_ring]->GetBinContent(j, i) - 1.);
0585       _h_drift_depth_[i_ring]->SetBinError(j, i, sqrt(error2));
0586     } else {
0587       _h_drift_depth_[i_ring]->SetBinError(j, i, 0);
0588       _h_drift_depth_adc_[i_ring]->SetBinError(j, i, 0);
0589     }
0590     h_drift_depth_adc_slice_->SetBinContent(j, _h_drift_depth_adc_[i_ring]->GetBinContent(j, i));
0591     h_drift_depth_adc_slice_->SetBinError(j, _h_drift_depth_adc_[i_ring]->GetBinError(j, i));
0592     nentries += _h_drift_depth_noadc_[i_ring]->GetBinContent(j, i);
0593   }  // end loop over bins in drift width
0594 
0595   double mean = h_drift_depth_adc_slice_->GetMean(1);
0596   double error = 0;
0597   if (nentries != 0) {
0598     error = h_drift_depth_adc_slice_->GetRMS(1) / sqrt(nentries);
0599   }
0600 
0601   _h_mean_[i_ring]->SetBinContent(i, mean);
0602   _h_mean_[i_ring]->SetBinError(i, error);
0603 }
0604 
0605 DEFINE_FWK_MODULE(SiPixelLorentzAngle);
[ 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!