Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /*  VI Janurary 2012 
0002  * This file need to be migrated to the new interface of matched hit as the mono/stero are not there anymore
0003  * what is returned are hits w/o localpoistion, just cluster and id
0004  */
0005 #include <string>
0006 #include <iostream>
0007 #include <fstream>
0008 
0009 #include "CalibTracker/SiStripLorentzAngle/interface/SiStripLAProfileBooker.h"
0010 
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0017 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0018 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0019 
0020 #include "DataFormats/Common/interface/Handle.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0023 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0024 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0025 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0027 
0028 #include "DQM/SiStripCommon/interface/ExtractTObject.h"
0029 #include "DQM/SiStripCommon/interface/SiStripHistoId.h"
0030 #include "DQM/SiStripCommon/interface/SiStripFolderOrganizer.h"
0031 #include "DQMServices/Core/interface/DQMStore.h"
0032 
0033 #include "MagneticField/Engine/interface/MagneticField.h"
0034 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0035 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0036 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0037 
0038 #include <TProfile.h>
0039 #include <TStyle.h>
0040 #include <TF1.h>
0041 
0042 #include <list>
0043 
0044 class DetIdLess {
0045 public:
0046   DetIdLess() {}
0047 
0048   bool operator()(const SiStripRecHit2D* a, const SiStripRecHit2D* b) const {
0049     return *(a->cluster()) < *(b->cluster());
0050   }
0051 };
0052 
0053 //Constructor
0054 
0055 SiStripLAProfileBooker::SiStripLAProfileBooker(edm::ParameterSet const& conf)
0056     : conf_(conf),
0057       tTopoToken_(esConsumes<edm::Transition::BeginRun>()),
0058       tkGeomToken_(esConsumes<edm::Transition::BeginRun>()),
0059       magFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0060       detCablingToken_(conf_.getParameter<bool>("UseStripCablingDB")
0061                            ? decltype(detCablingToken_){esConsumes<edm::Transition::BeginRun>()}
0062                            : decltype(detCablingToken_){}) {}
0063 
0064 //BeginRun
0065 
0066 void SiStripLAProfileBooker::beginRun(const edm::Run&, const edm::EventSetup& c) {
0067   const TrackerTopology* const tTopo = &c.getData(tTopoToken_);
0068   tkGeom_ = &c.getData(tkGeomToken_);
0069   const auto& magField = c.getData(magFieldToken_);
0070 
0071   std::vector<uint32_t> activeDets;
0072   if (conf_.getParameter<bool>("UseStripCablingDB")) {
0073     activeDets.clear();
0074     c.getData(detCablingToken_).addActiveDetectorsRawIds(activeDets);
0075   } else {
0076     const TrackerGeometry::DetIdContainer& Id = tkGeom_->detIds();
0077     TrackerGeometry::DetIdContainer::const_iterator Iditer;
0078     activeDets.clear();
0079     for (Iditer = Id.begin(); Iditer != Id.end(); Iditer++) {
0080       activeDets.push_back(Iditer->rawId());
0081     }
0082   }
0083 
0084   edm::InputTag TkTag = conf_.getParameter<edm::InputTag>("Tracks");
0085   //Get Ids;
0086   double ModuleRangeMin = conf_.getParameter<double>("ModuleXMin");
0087   double ModuleRangeMax = conf_.getParameter<double>("ModuleXMax");
0088   double TIBRangeMin = conf_.getParameter<double>("TIBXMin");
0089   double TIBRangeMax = conf_.getParameter<double>("TIBXMax");
0090   double TOBRangeMin = conf_.getParameter<double>("TOBXMin");
0091   double TOBRangeMax = conf_.getParameter<double>("TOBXMax");
0092   int TIB_bin = conf_.getParameter<int>("TIB_bin");
0093   int TOB_bin = conf_.getParameter<int>("TOB_bin");
0094   int SUM_bin = conf_.getParameter<int>("SUM_bin");
0095 
0096   hFile = new TFile(conf_.getUntrackedParameter<std::string>("treeName").c_str(), "RECREATE");
0097 
0098   Hit_Tree = hFile->mkdir("Hit_Tree");
0099   Track_Tree = hFile->mkdir("Track_Tree");
0100   Event_Tree = hFile->mkdir("Event_Tree");
0101 
0102   HitsTree = new TTree("HitsTree", "HitsTree");
0103 
0104   HitsTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
0105   HitsTree->Branch("EventNumber", &EventNumber, "EventNumber/I");
0106   HitsTree->Branch("TanTrackAngle", &TanTrackAngle, "TanTrackAngle/F");
0107   HitsTree->Branch("TanTrackAngleParallel", &TanTrackAngleParallel, "TanTrackAngleParallel/F");
0108   HitsTree->Branch("ClSize", &ClSize, "ClSize/I");
0109   HitsTree->Branch("HitCharge", &HitCharge, "HitCharge/I");
0110   HitsTree->Branch("Hit_Std_Dev", &hit_std_dev, "hit_std_dev/F");
0111   HitsTree->Branch("Type", &Type, "Type/I");
0112   HitsTree->Branch("Layer", &Layer, "Layer/I");
0113   HitsTree->Branch("Wheel", &Wheel, "Wheel/I");
0114   HitsTree->Branch("bw_fw", &bw_fw, "bw_fw/I");
0115   HitsTree->Branch("Ext_Int", &Ext_Int, "Ext_Int/I");
0116   HitsTree->Branch("MonoStereo", &MonoStereo, "MonoStereo/I");
0117   HitsTree->Branch("MagField", &MagField, "MagField/F");
0118   HitsTree->Branch("SignCorrection", &SignCorrection, "SignCorrection/F");
0119   HitsTree->Branch("XGlobal", &XGlobal, "XGlobal/F");
0120   HitsTree->Branch("YGlobal", &YGlobal, "YGlobal/F");
0121   HitsTree->Branch("ZGlobal", &ZGlobal, "ZGlobal/F");
0122   HitsTree->Branch("ParticleCharge", &ParticleCharge, "ParticleCharge/I");
0123   HitsTree->Branch("Momentum", &Momentum, "Momentum/F");
0124   HitsTree->Branch("pt", &pt, "pt/F");
0125   HitsTree->Branch("chi2norm", &chi2norm, "chi2norm/F");
0126   HitsTree->Branch("EtaTrack", &EtaTrack, "EtaTrack/F");
0127   HitsTree->Branch("PhiTrack", &PhiTrack, "PhiTrack/F");
0128   HitsTree->Branch("TrajSize", &trajsize, "trajsize/I");
0129   HitsTree->Branch("HitNr", &HitNr, "HitNr/I");
0130   HitsTree->Branch("HitPerTrack", &HitPerTrack, "HitPerTrack/I");
0131   HitsTree->Branch("id_detector", &id_detector, "id_detector/I");
0132   HitsTree->Branch("thick_detector", &thick_detector, "thick_detector/F");
0133   HitsTree->Branch("pitch_detector", &pitch_detector, "pitch_detector/F");
0134   HitsTree->Branch("Amplitudes", Amplitudes, "Amplitudes[ClSize]/I");
0135 
0136   HitsTree->SetDirectory(Hit_Tree);
0137 
0138   TrackTree = new TTree("TrackTree", "TrackTree");
0139 
0140   TrackTree->Branch("TrackCounter", &TrackCounter, "TrackCounter/I");
0141 
0142   TrackTree->SetDirectory(Track_Tree);
0143 
0144   EventTree = new TTree("EventTree", "EventTree");
0145 
0146   EventTree->Branch("EventCounter", &EventCounter, "EventCounter/I");
0147 
0148   EventTree->SetDirectory(Event_Tree);
0149 
0150   // use SistripHistoId for producing histogram id (and title)
0151   SiStripHistoId hidmanager;
0152 
0153   // create SiStripFolderOrganizer
0154   SiStripFolderOrganizer folder_organizer;
0155 
0156   dbe_ = edm::Service<DQMStore>().operator->();
0157 
0158   //get all detids
0159 
0160   for (std::vector<uint32_t>::const_iterator Id = activeDets.begin(); Id != activeDets.end(); Id++) {
0161     //  for(Iditer=Id.begin();Iditer!=Id.end();Iditer++){ //loop on detids
0162     DetId Iditero = DetId(*Id);
0163     DetId* Iditer = &Iditero;
0164     if ((Iditer->subdetId() == int(StripSubdetector::TIB)) ||
0165         (Iditer->subdetId() == int(StripSubdetector::TOB))) {  //include only barrel
0166 
0167       int module_bin = 0;
0168       if (Iditer->subdetId() == int(StripSubdetector::TIB)) {
0169         module_bin = TIB_bin;
0170       } else {
0171         module_bin = TOB_bin;
0172       }
0173 
0174       // create a TProfile for each module
0175       StripSubdetector subid(*Iditer);
0176       std::string hid;
0177       //Mono single sided detectors
0178       LocalPoint p;
0179       auto stripdet = tkGeom_->idToDet(subid);
0180       if (!stripdet->isLeaf())
0181         continue;
0182       const StripTopology& topol = (const StripTopology&)stripdet->topology();
0183       float thickness = stripdet->specificSurface().bounds().thickness();
0184 
0185       folder_organizer.setDetectorFolder(Iditer->rawId(), tTopo);
0186       hid = hidmanager.createHistoId(TkTag.label(), "det", Iditer->rawId());
0187       MonitorElement* profile = dbe_->bookProfile(hid, hid, module_bin, ModuleRangeMin, ModuleRangeMax, 20, 0, 5, "");
0188       detparameters* param = new detparameters;
0189       histos[Iditer->rawId()] = profile;
0190       detmap[Iditer->rawId()] = param;
0191       param->thickness = thickness * 10000;
0192       param->pitch = topol.localPitch(p) * 10000;
0193 
0194       const GlobalPoint globalp = (stripdet->surface()).toGlobal(p);
0195       GlobalVector globalmagdir = magField.inTesla(globalp);
0196       param->magfield = (stripdet->surface()).toLocal(globalmagdir);
0197 
0198       profile->setAxisTitle("tan(#theta_{t})", 1);
0199       profile->setAxisTitle("Cluster size", 2);
0200 
0201       // create a summary histo if it does not exist already
0202       std::string name;
0203       unsigned int layerid;
0204       getlayer(subid, tTopo, name, layerid);
0205       name += TkTag.label();
0206       if (summaryhisto.find(layerid) == (summaryhisto.end())) {
0207         folder_organizer.setSiStripFolder();
0208         MonitorElement* summaryprofile = nullptr;
0209         if (subid.subdetId() == int(StripSubdetector::TIB) || subid.subdetId() == int(StripSubdetector::TID))
0210           summaryprofile = dbe_->bookProfile(name, name, SUM_bin, TIBRangeMin, TIBRangeMax, 20, 0, 5, "");
0211         else if (subid.subdetId() == int(StripSubdetector::TOB) || subid.subdetId() == int(StripSubdetector::TEC))
0212           summaryprofile = dbe_->bookProfile(name, name, SUM_bin, TOBRangeMin, TOBRangeMax, 20, 0, 5, "");
0213         if (summaryprofile) {
0214           detparameters* summaryparam = new detparameters;
0215           summaryhisto[layerid] = summaryprofile;
0216           summarydetmap[layerid] = summaryparam;
0217           summaryparam->thickness = thickness * 10000;
0218           summaryparam->pitch = topol.localPitch(p) * 10000;
0219           summaryprofile->setAxisTitle("tan(#theta_{t})", 1);
0220           summaryprofile->setAxisTitle("Cluster size", 2);
0221         }
0222       }
0223     }
0224   }
0225 
0226   trackcollsize = 0;
0227   trajsize = 0;
0228   RunNumber = 0;
0229   EventNumber = 0;
0230   hitcounter = 0;
0231   hitcounter_2ndloop = 0;
0232   worse_double_hit = 0;
0233   better_double_hit = 0;
0234   eventcounter = 0;
0235 
0236   EventCounter = 1;
0237   TrackCounter = 1;
0238 }
0239 
0240 void SiStripLAProfileBooker::endRun(const edm::Run&, const edm::EventSetup& c) {}
0241 
0242 SiStripLAProfileBooker::~SiStripLAProfileBooker() {
0243   detparmap::iterator detpariter;
0244   for (detpariter = detmap.begin(); detpariter != detmap.end(); ++detpariter)
0245     delete detpariter->second;
0246   for (detpariter = summarydetmap.begin(); detpariter != summarydetmap.end(); ++detpariter)
0247     delete detpariter->second;
0248   delete hFile;
0249 }
0250 
0251 // Analyzer: Functions that gets called by framework every event
0252 
0253 void SiStripLAProfileBooker::analyze(const edm::Event& e, const edm::EventSetup& es) {
0254   const TrackerTopology* const tTopo = &es.getData(tTopoToken_);
0255 
0256   RunNumber = e.id().run();
0257   EventNumber = e.id().event();
0258 
0259   eventcounter++;
0260 
0261   EventTree->Fill();
0262 
0263   //Analysis of Trajectory-RecHits
0264 
0265   edm::InputTag TkTag = conf_.getParameter<edm::InputTag>("Tracks");
0266 
0267   edm::Handle<reco::TrackCollection> trackCollection;
0268   e.getByLabel(TkTag, trackCollection);
0269 
0270   edm::Handle<std::vector<Trajectory> > TrajectoryCollection;
0271   e.getByLabel(TkTag, TrajectoryCollection);
0272 
0273   edm::Handle<TrajTrackAssociationCollection> TrajTrackMap;
0274   e.getByLabel(TkTag, TrajTrackMap);
0275 
0276   const reco::TrackCollection* tracks = trackCollection.product();
0277 
0278   // FIXME this has to be changed to use pointers to clusters...
0279   std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess> hitangleassociation;
0280   std::list<SiStripRecHit2D> cache;  // ugly, inefficient, effective in making the above working
0281 
0282   trackcollsize = tracks->size();
0283   trajsize = TrajectoryCollection->size();
0284 
0285   edm::LogInfo("SiStripLAProfileBooker::analyze") << " Number of tracks in event = " << trackcollsize << "\n";
0286   edm::LogInfo("SiStripLAProfileBooker::analyze") << " Number of trajectories in event = " << trajsize << "\n";
0287 
0288   TrajTrackAssociationCollection::const_iterator TrajTrackIter;
0289 
0290   for (TrajTrackIter = TrajTrackMap->begin(); TrajTrackIter != TrajTrackMap->end();
0291        TrajTrackIter++) {  //loop on trajectories
0292 
0293     if (TrajTrackIter->key->foundHits() >= 5) {
0294       TrackTree->Fill();
0295 
0296       ParticleCharge = -99;
0297       Momentum = -99;
0298       pt = -99;
0299       chi2norm = -99;
0300       HitPerTrack = -99;
0301       EtaTrack = -99;
0302       PhiTrack = -99;
0303 
0304       ParticleCharge = TrajTrackIter->val->charge();
0305       pt = TrajTrackIter->val->pt();
0306       Momentum = TrajTrackIter->val->p();
0307       chi2norm = TrajTrackIter->val->normalizedChi2();
0308       EtaTrack = TrajTrackIter->val->eta();
0309       PhiTrack = TrajTrackIter->val->phi();
0310       HitPerTrack = TrajTrackIter->key->foundHits();
0311 
0312       std::vector<TrajectoryMeasurement> TMeas = TrajTrackIter->key->measurements();
0313       std::vector<TrajectoryMeasurement>::iterator itm;
0314 
0315       for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {  //loop on hits
0316 
0317         int i;
0318         for (i = 0; i < 100; i++) {
0319           Amplitudes[i] = 0;
0320         }
0321 
0322         TanTrackAngle = -99;
0323         TanTrackAngleParallel = -99;
0324         ClSize = -99;
0325         HitCharge = 0;
0326         Type = -99;
0327         Layer = -99;
0328         Wheel = -99;
0329         bw_fw = -99;
0330         Ext_Int = -99;
0331         MonoStereo = -99;
0332         MagField = -99;
0333         SignCorrection = -99;
0334         XGlobal = -99;
0335         YGlobal = -99;
0336         ZGlobal = -99;
0337         barycenter = -99;
0338         hit_std_dev = -99;
0339         sumx = 0;
0340         id_detector = -1;
0341         thick_detector = -1;
0342         pitch_detector = -1;
0343         HitNr = 1;
0344 
0345         SiStripRecHit2D lhit;
0346         TrajectoryStateOnSurface tsos = itm->updatedState();
0347         const TransientTrackingRecHit::ConstRecHitPointer thit = itm->recHit();
0348         if ((thit->geographicalId().subdetId() == int(StripSubdetector::TIB)) ||
0349             thit->geographicalId().subdetId() == int(StripSubdetector::TOB)) {  //include only barrel
0350           const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>((*thit).hit());
0351           const ProjectedSiStripRecHit2D* phit = dynamic_cast<const ProjectedSiStripRecHit2D*>((*thit).hit());
0352           const SiStripRecHit2D* hit = dynamic_cast<const SiStripRecHit2D*>((*thit).hit());
0353           if (phit) {
0354             lhit = phit->originalHit();
0355             hit = &lhit;
0356           }
0357 
0358           LocalVector trackdirection = tsos.localDirection();
0359 
0360           if (matchedhit) {  //if matched hit...
0361 
0362             GluedGeomDet* gdet = (GluedGeomDet*)tkGeom_->idToDet(matchedhit->geographicalId());
0363 
0364             GlobalVector gtrkdir = gdet->toGlobal(trackdirection);
0365 
0366             // THIS THE POINTER TO THE MONO HIT OF A MATCHED HIT
0367 
0368             // top be migrated to the more direct interface of matchedhit
0369             cache.push_back(matchedhit->monoHit());
0370             const SiStripRecHit2D* monohit = &cache.back();
0371             const SiStripRecHit2D::ClusterRef& monocluster = monohit->cluster();
0372             const GeomDetUnit* monodet = gdet->monoDet();
0373             // this does not exists anymore! either project the matched or use CPE
0374             const LocalPoint monoposition = monohit->localPosition();
0375 
0376             StripSubdetector detid = (StripSubdetector)monohit->geographicalId();
0377             id_detector = detid.rawId();
0378             thick_detector = monodet->specificSurface().bounds().thickness();
0379             const StripTopology& mtopol = (StripTopology&)monodet->topology();
0380             pitch_detector = mtopol.localPitch(monoposition);
0381             const GlobalPoint monogposition = (monodet->surface()).toGlobal(monoposition);
0382             ClSize = (monocluster->amplitudes()).size();
0383 
0384             const auto& amplitudes = monocluster->amplitudes();
0385 
0386             barycenter = monocluster->barycenter() - 0.5;
0387             uint16_t FirstStrip = monocluster->firstStrip();
0388             auto begin = amplitudes.begin();
0389             nstrip = 0;
0390             for (auto idigi = begin; idigi != amplitudes.end(); idigi++) {
0391               Amplitudes[nstrip] = *idigi;
0392               sumx += pow(((FirstStrip + idigi - begin) - barycenter), 2) * (*idigi);
0393               HitCharge += *idigi;
0394             }
0395             hit_std_dev = sqrt(sumx / HitCharge);
0396 
0397             XGlobal = monogposition.x();
0398             YGlobal = monogposition.y();
0399             ZGlobal = monogposition.z();
0400 
0401             Type = detid.subdetId();
0402             MonoStereo = detid.stereo();
0403 
0404             if (detid.subdetId() == int(StripSubdetector::TIB)) {
0405               Layer = tTopo->tibLayer(detid);
0406               bw_fw = tTopo->tibStringInfo(detid)[0];
0407               Ext_Int = tTopo->tibStringInfo(detid)[1];
0408             }
0409             if (detid.subdetId() == int(StripSubdetector::TOB)) {
0410               Layer = tTopo->tobLayer(detid);
0411               bw_fw = tTopo->tobRodInfo(detid)[0];
0412             }
0413             if (detid.subdetId() == int(StripSubdetector::TID)) {
0414               Wheel = tTopo->tidWheel(detid);
0415               bw_fw = tTopo->tidModuleInfo(detid)[0];
0416             }
0417             if (detid.subdetId() == int(StripSubdetector::TEC)) {
0418               Wheel = tTopo->tecWheel(detid);
0419               bw_fw = tTopo->tecPetalInfo(detid)[0];
0420             }
0421 
0422             LocalVector monotkdir = monodet->toLocal(gtrkdir);
0423 
0424             if (monotkdir.z() != 0) {
0425               // THE LOCAL ANGLE (MONO)
0426               float tanangle = monotkdir.x() / monotkdir.z();
0427               TanTrackAngleParallel = monotkdir.y() / monotkdir.z();
0428               TanTrackAngle = tanangle;
0429               detparmap::iterator TheDet = detmap.find(detid.rawId());
0430               LocalVector localmagdir;
0431               if (TheDet != detmap.end())
0432                 localmagdir = TheDet->second->magfield;
0433               MagField = localmagdir.mag();
0434               if (MagField != 0.) {
0435                 LocalVector monoylocal(0, 1, 0);
0436                 float signcorrection = (localmagdir * monoylocal) / (MagField);
0437                 if (signcorrection != 0)
0438                   SignCorrection = 1 / signcorrection;
0439               }
0440 
0441               std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator alreadystored =
0442                   hitangleassociation.find(monohit);
0443 
0444               if (alreadystored != hitangleassociation.end()) {  //decide which hit take
0445                 if (itm->estimate() > alreadystored->second.first) {
0446                   worse_double_hit++;
0447                 }
0448                 if (itm->estimate() < alreadystored->second.first) {
0449                   better_double_hit++;
0450                   hitangleassociation.insert(std::make_pair(monohit, std::make_pair(itm->estimate(), tanangle)));
0451                 }
0452               } else {
0453                 hitangleassociation.insert(make_pair(monohit, std::make_pair(itm->estimate(), tanangle)));
0454                 HitsTree->Fill();
0455                 hitcounter++;
0456               }
0457 
0458               // THIS THE POINTER TO THE STEREO HIT OF A MATCHED HIT
0459 
0460               // top be migrated to the more direct interface of matchedhit
0461               cache.push_back(matchedhit->stereoHit());
0462               const SiStripRecHit2D* stereohit = &cache.back();
0463               const SiStripRecHit2D::ClusterRef& stereocluster = stereohit->cluster();
0464               const GeomDetUnit* stereodet = gdet->stereoDet();
0465               // this does not exists anymore! either project the matched or use CPE
0466               const LocalPoint stereoposition = stereohit->localPosition();
0467               StripSubdetector detid = (StripSubdetector)stereohit->geographicalId();
0468               id_detector = detid.rawId();
0469               thick_detector = stereodet->specificSurface().bounds().thickness();
0470               const StripTopology& stopol = (StripTopology&)stereodet->topology();
0471               pitch_detector = stopol.localPitch(stereoposition);
0472               const GlobalPoint stereogposition = (stereodet->surface()).toGlobal(stereoposition);
0473 
0474               ClSize = (stereocluster->amplitudes()).size();
0475 
0476               const auto& amplitudes = stereocluster->amplitudes();
0477 
0478               barycenter = stereocluster->barycenter() - 0.5;
0479               uint16_t FirstStrip = stereocluster->firstStrip();
0480               auto begin = amplitudes.begin();
0481               nstrip = 0;
0482               for (auto idigi = begin; idigi != amplitudes.end(); idigi++) {
0483                 Amplitudes[nstrip] = *idigi;
0484                 sumx += pow(((FirstStrip + idigi - begin) - barycenter), 2) * (*idigi);
0485                 HitCharge += *idigi;
0486               }
0487               hit_std_dev = sqrt(sumx / HitCharge);
0488 
0489               XGlobal = stereogposition.x();
0490               YGlobal = stereogposition.y();
0491               ZGlobal = stereogposition.z();
0492 
0493               Type = detid.subdetId();
0494               MonoStereo = detid.stereo();
0495 
0496               if (detid.subdetId() == int(StripSubdetector::TIB)) {
0497                 Layer = tTopo->tibLayer(detid);
0498                 bw_fw = tTopo->tibStringInfo(detid)[0];
0499                 Ext_Int = tTopo->tibStringInfo(detid)[1];
0500               }
0501               if (detid.subdetId() == int(StripSubdetector::TOB)) {
0502                 Layer = tTopo->tobLayer(detid);
0503                 bw_fw = tTopo->tobRodInfo(detid)[0];
0504               }
0505               if (detid.subdetId() == int(StripSubdetector::TID)) {
0506                 Wheel = tTopo->tidWheel(detid);
0507                 bw_fw = tTopo->tidModuleInfo(detid)[0];
0508               }
0509               if (detid.subdetId() == int(StripSubdetector::TEC)) {
0510                 Wheel = tTopo->tecWheel(detid);
0511                 bw_fw = tTopo->tecPetalInfo(detid)[0];
0512               }
0513 
0514               LocalVector stereotkdir = stereodet->toLocal(gtrkdir);
0515 
0516               if (stereotkdir.z() != 0) {
0517                 // THE LOCAL ANGLE (STEREO)
0518                 float tanangle = stereotkdir.x() / stereotkdir.z();
0519                 TanTrackAngleParallel = stereotkdir.y() / stereotkdir.z();
0520                 TanTrackAngle = tanangle;
0521                 detparmap::iterator TheDet = detmap.find(detid.rawId());
0522                 LocalVector localmagdir;
0523                 if (TheDet != detmap.end())
0524                   localmagdir = TheDet->second->magfield;
0525                 MagField = localmagdir.mag();
0526                 LocalVector stereoylocal(0, 1, 0);
0527                 if (MagField != 0.) {
0528                   float signcorrection = (localmagdir * stereoylocal) / (MagField);
0529                   if (signcorrection != 0)
0530                     SignCorrection = 1 / signcorrection;
0531                 }
0532 
0533                 std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator alreadystored =
0534                     hitangleassociation.find(stereohit);
0535 
0536                 if (alreadystored != hitangleassociation.end()) {  //decide which hit take
0537                   if (itm->estimate() > alreadystored->second.first) {
0538                     worse_double_hit++;
0539                   }
0540                   if (itm->estimate() < alreadystored->second.first) {
0541                     better_double_hit++;
0542                     hitangleassociation.insert(std::make_pair(stereohit, std::make_pair(itm->estimate(), tanangle)));
0543                   }
0544                 } else {
0545                   hitangleassociation.insert(std::make_pair(stereohit, std::make_pair(itm->estimate(), tanangle)));
0546                   HitsTree->Fill();
0547                   hitcounter++;
0548                 }
0549               }
0550             }
0551           } else if (hit) {
0552             //  hit= POINTER TO THE RECHIT
0553 
0554             const SiStripRecHit2D::ClusterRef& cluster = hit->cluster();
0555 
0556             GeomDetUnit* gdet = (GeomDetUnit*)tkGeom_->idToDet(hit->geographicalId());
0557             const LocalPoint position = hit->localPosition();
0558             StripSubdetector detid = (StripSubdetector)hit->geographicalId();
0559             id_detector = detid.rawId();
0560             thick_detector = gdet->specificSurface().bounds().thickness();
0561             const StripTopology& topol = (StripTopology&)gdet->topology();
0562             pitch_detector = topol.localPitch(position);
0563             const GlobalPoint gposition = (gdet->surface()).toGlobal(position);
0564 
0565             ClSize = (cluster->amplitudes()).size();
0566 
0567             const auto& amplitudes = cluster->amplitudes();
0568 
0569             barycenter = cluster->barycenter() - 0.5;
0570             uint16_t FirstStrip = cluster->firstStrip();
0571             nstrip = 0;
0572             auto begin = amplitudes.begin();
0573             for (auto idigi = amplitudes.begin(); idigi != amplitudes.end(); idigi++) {
0574               Amplitudes[nstrip] = *idigi;
0575               sumx += pow(((FirstStrip + idigi - begin) - barycenter), 2) * (*idigi);
0576               HitCharge += *idigi;
0577             }
0578             hit_std_dev = sqrt(sumx / HitCharge);
0579 
0580             XGlobal = gposition.x();
0581             YGlobal = gposition.y();
0582             ZGlobal = gposition.z();
0583 
0584             Type = detid.subdetId();
0585             MonoStereo = detid.stereo();
0586 
0587             if (detid.subdetId() == int(StripSubdetector::TIB)) {
0588               Layer = tTopo->tibLayer(detid);
0589               bw_fw = tTopo->tibStringInfo(detid)[0];
0590               Ext_Int = tTopo->tibStringInfo(detid)[1];
0591             }
0592             if (detid.subdetId() == int(StripSubdetector::TOB)) {
0593               Layer = tTopo->tobLayer(detid);
0594               bw_fw = tTopo->tobRodInfo(detid)[0];
0595             }
0596             if (detid.subdetId() == int(StripSubdetector::TID)) {
0597               Wheel = tTopo->tidWheel(detid);
0598               bw_fw = tTopo->tidModuleInfo(detid)[0];
0599             }
0600             if (detid.subdetId() == int(StripSubdetector::TEC)) {
0601               Wheel = tTopo->tecWheel(detid);
0602               bw_fw = tTopo->tecPetalInfo(detid)[0];
0603             }
0604 
0605             if (trackdirection.z() != 0) {
0606               // THE LOCAL ANGLE
0607               float tanangle = trackdirection.x() / trackdirection.z();
0608               TanTrackAngleParallel = trackdirection.y() / trackdirection.z();
0609               TanTrackAngle = tanangle;
0610               detparmap::iterator TheDet = detmap.find(detid.rawId());
0611               LocalVector localmagdir;
0612               if (TheDet != detmap.end())
0613                 localmagdir = TheDet->second->magfield;
0614               MagField = localmagdir.mag();
0615               if (MagField != 0.) {
0616                 LocalVector ylocal(0, 1, 0);
0617                 float signcorrection = (localmagdir * ylocal) / (MagField);
0618                 if (signcorrection != 0)
0619                   SignCorrection = 1 / signcorrection;
0620               }
0621 
0622               std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator alreadystored =
0623                   hitangleassociation.find(hit);
0624 
0625               if (alreadystored != hitangleassociation.end()) {  //decide which hit take
0626                 if (itm->estimate() > alreadystored->second.first) {
0627                   worse_double_hit++;
0628                 }
0629                 if (itm->estimate() < alreadystored->second.first) {
0630                   better_double_hit++;
0631                   hitangleassociation.insert(std::make_pair(hit, std::make_pair(itm->estimate(), tanangle)));
0632                 }
0633               } else {
0634                 hitangleassociation.insert(std::make_pair(hit, std::make_pair(itm->estimate(), tanangle)));
0635                 HitsTree->Fill();
0636                 hitcounter++;
0637               }
0638             }
0639           }
0640         }
0641       }
0642     }
0643   }
0644   std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator hitsiter;
0645 
0646   for (hitsiter = hitangleassociation.begin(); hitsiter != hitangleassociation.end(); hitsiter++) {
0647     hitcounter_2ndloop++;
0648 
0649     const SiStripRecHit2D* hit = hitsiter->first;
0650     const SiStripRecHit2D::ClusterRef& cluster = hit->cluster();
0651 
0652     size = (cluster->amplitudes()).size();
0653 
0654     StripSubdetector detid = (StripSubdetector)hit->geographicalId();
0655 
0656     float tangent = hitsiter->second.second;
0657 
0658     //Sign and XZ plane projection correction applied in TrackLocalAngle (TIB|TOB layers)
0659 
0660     detparmap::iterator thedet = detmap.find(detid.rawId());
0661     LocalVector localmagdir;
0662     if (thedet != detmap.end())
0663       localmagdir = thedet->second->magfield;
0664     float localmagfield = localmagdir.mag();
0665 
0666     if (localmagfield != 0.) {
0667       LocalVector ylocal(0, 1, 0);
0668 
0669       float normprojection = (localmagdir * ylocal) / (localmagfield);
0670 
0671       if (normprojection == 0.)
0672         LogDebug("SiStripLAProfileBooker::analyze") << "Error: YBprojection = 0";
0673 
0674       else {
0675         float signprojcorrection = 1 / normprojection;
0676         tangent *= signprojcorrection;
0677       }
0678     }
0679 
0680     //Filling histograms
0681 
0682     histomap::iterator thehisto = histos.find(detid.rawId());
0683 
0684     if (thehisto == histos.end())
0685       edm::LogError("SiStripLAProfileBooker::analyze")
0686           << "Error: the profile associated to" << detid.rawId() << "does not exist! ";
0687     else
0688       thehisto->second->Fill(tangent, size);
0689 
0690     //Summary histograms
0691     std::string name;
0692     unsigned int layerid;
0693     getlayer(detid, tTopo, name, layerid);
0694     histomap::iterator thesummaryhisto = summaryhisto.find(layerid);
0695     if (thesummaryhisto == summaryhisto.end())
0696       edm::LogError("SiStripLAProfileBooker::analyze")
0697           << "Error: the profile associated to subdet " << name << "does not exist! ";
0698     else
0699       thesummaryhisto->second->Fill(tangent, size);
0700   }
0701 }
0702 
0703 //Makename function
0704 
0705 void SiStripLAProfileBooker::getlayer(const DetId& detid,
0706                                       const TrackerTopology* tTopo,
0707                                       std::string& name,
0708                                       unsigned int& layerid) {
0709   int layer = 0;
0710   std::stringstream layernum;
0711 
0712   if (detid.subdetId() == int(StripSubdetector::TIB)) {
0713     name += "TIB_Layer_";
0714     layer = tTopo->tibLayer(detid);
0715   }
0716 
0717   else if (detid.subdetId() == int(StripSubdetector::TID)) {
0718     name += "TID_Ring_";
0719     layer = tTopo->tidRing(detid);
0720   }
0721 
0722   else if (detid.subdetId() == int(StripSubdetector::TOB)) {
0723     name += "TOB_Layer_";
0724     layer = tTopo->tobLayer(detid);
0725 
0726   }
0727 
0728   else if (detid.subdetId() == int(StripSubdetector::TEC)) {
0729     name += "TEC_Ring_";
0730     layer = tTopo->tecRing(detid);
0731   }
0732   layernum << layer;
0733   name += layernum.str();
0734   layerid = detid.subdetId() * 10 + layer;
0735 }
0736 
0737 void SiStripLAProfileBooker::endJob() {
0738   std::string outputFile_ = conf_.getUntrackedParameter<std::string>("fileName", "LorentzAngle.root");
0739   dbe_->save(outputFile_);
0740 
0741   hFile->Write();
0742   hFile->Close();
0743 }
[ 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!