SiStripLAProfileBooker.cc

Go to the documentation of this file.

/*  VI Janurary 2012 
 * This file need to be migrated to the new interface of matched hit as the mono/stero are not there anymore
 * what is returned are hits w/o localpoistion, just cluster and id
 */
#include <string>
#include <iostream>
#include <fstream>
 
#include "CalibTracker/SiStripLorentzAngle/interface/SiStripLAProfileBooker.h"
 
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
#include "Geometry/CommonTopologies/interface/StripTopology.h"
 
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
#include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
 
#include "DQM/SiStripCommon/interface/ExtractTObject.h"
#include "DQM/SiStripCommon/interface/SiStripHistoId.h"
#include "DQM/SiStripCommon/interface/SiStripFolderOrganizer.h"
#include "DQMServices/Core/interface/DQMStore.h"
 
#include "MagneticField/Engine/interface/MagneticField.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
 
#include <TProfile.h>
#include <TStyle.h>
#include <TF1.h>
 
#include <list>
 
class DetIdLess {
public:
  DetIdLess() {}
 
  bool operator()(const SiStripRecHit2D* a, const SiStripRecHit2D* b) const {
    return *(a->cluster()) < *(b->cluster());
  }
};
 
//Constructor
 
SiStripLAProfileBooker::SiStripLAProfileBooker(edm::ParameterSet const& conf)
    : conf_(conf),
      tTopoToken_(esConsumes<edm::Transition::BeginRun>()),
      tkGeomToken_(esConsumes<edm::Transition::BeginRun>()),
      magFieldToken_(esConsumes<edm::Transition::BeginRun>()),
      detCablingToken_(conf_.getParameter<bool>("UseStripCablingDB")
                           ? decltype(detCablingToken_){esConsumes<edm::Transition::BeginRun>()}
                           : decltype(detCablingToken_){}) {}
 
//BeginRun
 
void SiStripLAProfileBooker::beginRun(const edm::Run&, const edm::EventSetup& c) {
  const TrackerTopology* const tTopo = &c.getData(tTopoToken_);
  tkGeom_ = &c.getData(tkGeomToken_);
  const auto& magField = c.getData(magFieldToken_);
 
  std::vector<uint32_t> activeDets;
  if (conf_.getParameter<bool>("UseStripCablingDB")) {
    activeDets.clear();
    c.getData(detCablingToken_).addActiveDetectorsRawIds(activeDets);
  } else {
    const TrackerGeometry::DetIdContainer& Id = tkGeom_->detIds();
    TrackerGeometry::DetIdContainer::const_iterator Iditer;
    activeDets.clear();
    for (Iditer = Id.begin(); Iditer != Id.end(); Iditer++) {
      activeDets.push_back(Iditer->rawId());
    }
  }
 
  edm::InputTag TkTag = conf_.getParameter<edm::InputTag>("Tracks");
  //Get Ids;
  double ModuleRangeMin = conf_.getParameter<double>("ModuleXMin");
  double ModuleRangeMax = conf_.getParameter<double>("ModuleXMax");
  double TIBRangeMin = conf_.getParameter<double>("TIBXMin");
  double TIBRangeMax = conf_.getParameter<double>("TIBXMax");
  double TOBRangeMin = conf_.getParameter<double>("TOBXMin");
  double TOBRangeMax = conf_.getParameter<double>("TOBXMax");
  int TIB_bin = conf_.getParameter<int>("TIB_bin");
  int TOB_bin = conf_.getParameter<int>("TOB_bin");
  int SUM_bin = conf_.getParameter<int>("SUM_bin");
 
  hFile = new TFile(conf_.getUntrackedParameter<std::string>("treeName").c_str(), "RECREATE");
 
  Hit_Tree = hFile->mkdir("Hit_Tree");
  Track_Tree = hFile->mkdir("Track_Tree");
  Event_Tree = hFile->mkdir("Event_Tree");
 
  HitsTree = new TTree("HitsTree", "HitsTree");
 
  HitsTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
  HitsTree->Branch("EventNumber", &EventNumber, "EventNumber/I");
  HitsTree->Branch("TanTrackAngle", &TanTrackAngle, "TanTrackAngle/F");
  HitsTree->Branch("TanTrackAngleParallel", &TanTrackAngleParallel, "TanTrackAngleParallel/F");
  HitsTree->Branch("ClSize", &ClSize, "ClSize/I");
  HitsTree->Branch("HitCharge", &HitCharge, "HitCharge/I");
  HitsTree->Branch("Hit_Std_Dev", &hit_std_dev, "hit_std_dev/F");
  HitsTree->Branch("Type", &Type, "Type/I");
  HitsTree->Branch("Layer", &Layer, "Layer/I");
  HitsTree->Branch("Wheel", &Wheel, "Wheel/I");
  HitsTree->Branch("bw_fw", &bw_fw, "bw_fw/I");
  HitsTree->Branch("Ext_Int", &Ext_Int, "Ext_Int/I");
  HitsTree->Branch("MonoStereo", &MonoStereo, "MonoStereo/I");
  HitsTree->Branch("MagField", &MagField, "MagField/F");
  HitsTree->Branch("SignCorrection", &SignCorrection, "SignCorrection/F");
  HitsTree->Branch("XGlobal", &XGlobal, "XGlobal/F");
  HitsTree->Branch("YGlobal", &YGlobal, "YGlobal/F");
  HitsTree->Branch("ZGlobal", &ZGlobal, "ZGlobal/F");
  HitsTree->Branch("ParticleCharge", &ParticleCharge, "ParticleCharge/I");
  HitsTree->Branch("Momentum", &Momentum, "Momentum/F");
  HitsTree->Branch("pt", &pt, "pt/F");
  HitsTree->Branch("chi2norm", &chi2norm, "chi2norm/F");
  HitsTree->Branch("EtaTrack", &EtaTrack, "EtaTrack/F");
  HitsTree->Branch("PhiTrack", &PhiTrack, "PhiTrack/F");
  HitsTree->Branch("TrajSize", &trajsize, "trajsize/I");
  HitsTree->Branch("HitNr", &HitNr, "HitNr/I");
  HitsTree->Branch("HitPerTrack", &HitPerTrack, "HitPerTrack/I");
  HitsTree->Branch("id_detector", &id_detector, "id_detector/I");
  HitsTree->Branch("thick_detector", &thick_detector, "thick_detector/F");
  HitsTree->Branch("pitch_detector", &pitch_detector, "pitch_detector/F");
  HitsTree->Branch("Amplitudes", Amplitudes, "Amplitudes[ClSize]/I");
 
  HitsTree->SetDirectory(Hit_Tree);
 
  TrackTree = new TTree("TrackTree", "TrackTree");
 
  TrackTree->Branch("TrackCounter", &TrackCounter, "TrackCounter/I");
 
  TrackTree->SetDirectory(Track_Tree);
 
  EventTree = new TTree("EventTree", "EventTree");
 
  EventTree->Branch("EventCounter", &EventCounter, "EventCounter/I");
 
  EventTree->SetDirectory(Event_Tree);
 
  // use SistripHistoId for producing histogram id (and title)
  SiStripHistoId hidmanager;
 
  // create SiStripFolderOrganizer
  SiStripFolderOrganizer folder_organizer;
 
  dbe_ = edm::Service<DQMStore>().operator->();
 
  //get all detids
 
  for (std::vector<uint32_t>::const_iterator Id = activeDets.begin(); Id != activeDets.end(); Id++) {
    //  for(Iditer=Id.begin();Iditer!=Id.end();Iditer++){ //loop on detids
    DetId Iditero = DetId(*Id);
    DetId* Iditer = &Iditero;
    if ((Iditer->subdetId() == int(StripSubdetector::TIB)) ||
        (Iditer->subdetId() == int(StripSubdetector::TOB))) {  //include only barrel
 
      int module_bin = 0;
      if (Iditer->subdetId() == int(StripSubdetector::TIB)) {
        module_bin = TIB_bin;
      } else {
        module_bin = TOB_bin;
      }
 
      // create a TProfile for each module
      StripSubdetector subid(*Iditer);
      std::string hid;
      //Mono single sided detectors
      LocalPoint p;
      auto stripdet = tkGeom_->idToDet(subid);
      if (!stripdet->isLeaf())
        continue;
      const StripTopology& topol = (const StripTopology&)stripdet->topology();
      float thickness = stripdet->specificSurface().bounds().thickness();
 
      folder_organizer.setDetectorFolder(Iditer->rawId(), tTopo);
      hid = hidmanager.createHistoId(TkTag.label(), "det", Iditer->rawId());
      MonitorElement* profile = dbe_->bookProfile(hid, hid, module_bin, ModuleRangeMin, ModuleRangeMax, 20, 0, 5, "");
      detparameters* param = new detparameters;
      histos[Iditer->rawId()] = profile;
      detmap[Iditer->rawId()] = param;
      param->thickness = thickness * 10000;
      param->pitch = topol.localPitch(p) * 10000;
 
      const GlobalPoint globalp = (stripdet->surface()).toGlobal(p);
      GlobalVector globalmagdir = magField.inTesla(globalp);
      param->magfield = (stripdet->surface()).toLocal(globalmagdir);
 
      profile->setAxisTitle("tan(#theta_{t})", 1);
      profile->setAxisTitle("Cluster size", 2);
 
      // create a summary histo if it does not exist already
      std::string name;
      unsigned int layerid;
      getlayer(subid, tTopo, name, layerid);
      name += TkTag.label();
      if (summaryhisto.find(layerid) == (summaryhisto.end())) {
        folder_organizer.setSiStripFolder();
        MonitorElement* summaryprofile = nullptr;
        if (subid.subdetId() == int(StripSubdetector::TIB) || subid.subdetId() == int(StripSubdetector::TID))
          summaryprofile = dbe_->bookProfile(name, name, SUM_bin, TIBRangeMin, TIBRangeMax, 20, 0, 5, "");
        else if (subid.subdetId() == int(StripSubdetector::TOB) || subid.subdetId() == int(StripSubdetector::TEC))
          summaryprofile = dbe_->bookProfile(name, name, SUM_bin, TOBRangeMin, TOBRangeMax, 20, 0, 5, "");
        if (summaryprofile) {
          detparameters* summaryparam = new detparameters;
          summaryhisto[layerid] = summaryprofile;
          summarydetmap[layerid] = summaryparam;
          summaryparam->thickness = thickness * 10000;
          summaryparam->pitch = topol.localPitch(p) * 10000;
          summaryprofile->setAxisTitle("tan(#theta_{t})", 1);
          summaryprofile->setAxisTitle("Cluster size", 2);
        }
      }
    }
  }
 
  trackcollsize = 0;
  trajsize = 0;
  RunNumber = 0;
  EventNumber = 0;
  hitcounter = 0;
  hitcounter_2ndloop = 0;
  worse_double_hit = 0;
  better_double_hit = 0;
  eventcounter = 0;
 
  EventCounter = 1;
  TrackCounter = 1;
}
 
SiStripLAProfileBooker::~SiStripLAProfileBooker() {
  detparmap::iterator detpariter;
  for (detpariter = detmap.begin(); detpariter != detmap.end(); ++detpariter)
    delete detpariter->second;
  for (detpariter = summarydetmap.begin(); detpariter != summarydetmap.end(); ++detpariter)
    delete detpariter->second;
  delete hFile;
}
 
// Analyzer: Functions that gets called by framework every event
 
void SiStripLAProfileBooker::analyze(const edm::Event& e, const edm::EventSetup& es) {
  const TrackerTopology* const tTopo = &es.getData(tTopoToken_);
 
  RunNumber = e.id().run();
  EventNumber = e.id().event();
 
  eventcounter++;
 
  EventTree->Fill();
 
  //Analysis of Trajectory-RecHits
 
  edm::InputTag TkTag = conf_.getParameter<edm::InputTag>("Tracks");
 
  edm::Handle<reco::TrackCollection> trackCollection;
  e.getByLabel(TkTag, trackCollection);
 
  edm::Handle<std::vector<Trajectory> > TrajectoryCollection;
  e.getByLabel(TkTag, TrajectoryCollection);
 
  edm::Handle<TrajTrackAssociationCollection> TrajTrackMap;
  e.getByLabel(TkTag, TrajTrackMap);
 
  const reco::TrackCollection* tracks = trackCollection.product();
 
  // FIXME this has to be changed to use pointers to clusters...
  std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess> hitangleassociation;
  std::list<SiStripRecHit2D> cache;  // ugly, inefficient, effective in making the above working
 
  trackcollsize = tracks->size();
  trajsize = TrajectoryCollection->size();
 
  edm::LogInfo("SiStripLAProfileBooker::analyze") << " Number of tracks in event = " << trackcollsize << "\n";
  edm::LogInfo("SiStripLAProfileBooker::analyze") << " Number of trajectories in event = " << trajsize << "\n";
 
  TrajTrackAssociationCollection::const_iterator TrajTrackIter;
 
  for (TrajTrackIter = TrajTrackMap->begin(); TrajTrackIter != TrajTrackMap->end();
       TrajTrackIter++) {  //loop on trajectories
 
    if (TrajTrackIter->key->foundHits() >= 5) {
      TrackTree->Fill();
 
      ParticleCharge = -99;
      Momentum = -99;
      pt = -99;
      chi2norm = -99;
      HitPerTrack = -99;
      EtaTrack = -99;
      PhiTrack = -99;
 
      ParticleCharge = TrajTrackIter->val->charge();
      pt = TrajTrackIter->val->pt();
      Momentum = TrajTrackIter->val->p();
      chi2norm = TrajTrackIter->val->normalizedChi2();
      EtaTrack = TrajTrackIter->val->eta();
      PhiTrack = TrajTrackIter->val->phi();
      HitPerTrack = TrajTrackIter->key->foundHits();
 
      std::vector<TrajectoryMeasurement> TMeas = TrajTrackIter->key->measurements();
      std::vector<TrajectoryMeasurement>::iterator itm;
 
      for (itm = TMeas.begin(); itm != TMeas.end(); itm++) {  //loop on hits
 
        int i;
        for (i = 0; i < 100; i++) {
          Amplitudes[i] = 0;
        }
 
        TanTrackAngle = -99;
        TanTrackAngleParallel = -99;
        ClSize = -99;
        HitCharge = 0;
        Type = -99;
        Layer = -99;
        Wheel = -99;
        bw_fw = -99;
        Ext_Int = -99;
        MonoStereo = -99;
        MagField = -99;
        SignCorrection = -99;
        XGlobal = -99;
        YGlobal = -99;
        ZGlobal = -99;
        barycenter = -99;
        hit_std_dev = -99;
        sumx = 0;
        id_detector = -1;
        thick_detector = -1;
        pitch_detector = -1;
        HitNr = 1;
 
        SiStripRecHit2D lhit;
        TrajectoryStateOnSurface tsos = itm->updatedState();
        const TransientTrackingRecHit::ConstRecHitPointer thit = itm->recHit();
        if ((thit->geographicalId().subdetId() == int(StripSubdetector::TIB)) ||
            thit->geographicalId().subdetId() == int(StripSubdetector::TOB)) {  //include only barrel
          const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>((*thit).hit());
          const ProjectedSiStripRecHit2D* phit = dynamic_cast<const ProjectedSiStripRecHit2D*>((*thit).hit());
          const SiStripRecHit2D* hit = dynamic_cast<const SiStripRecHit2D*>((*thit).hit());
          if (phit) {
            lhit = phit->originalHit();
            hit = &lhit;
          }
 
          LocalVector trackdirection = tsos.localDirection();
 
          if (matchedhit) {  //if matched hit...
 
            GluedGeomDet* gdet = (GluedGeomDet*)tkGeom_->idToDet(matchedhit->geographicalId());
 
            GlobalVector gtrkdir = gdet->toGlobal(trackdirection);
 
            // THIS THE POINTER TO THE MONO HIT OF A MATCHED HIT
 
            // top be migrated to the more direct interface of matchedhit
            cache.push_back(matchedhit->monoHit());
            const SiStripRecHit2D* monohit = &cache.back();
            const SiStripRecHit2D::ClusterRef& monocluster = monohit->cluster();
            const GeomDetUnit* monodet = gdet->monoDet();
            // this does not exists anymore! either project the matched or use CPE
            const LocalPoint monoposition = monohit->localPosition();
 
            StripSubdetector detid = (StripSubdetector)monohit->geographicalId();
            id_detector = detid.rawId();
            thick_detector = monodet->specificSurface().bounds().thickness();
            const StripTopology& mtopol = (StripTopology&)monodet->topology();
            pitch_detector = mtopol.localPitch(monoposition);
            const GlobalPoint monogposition = (monodet->surface()).toGlobal(monoposition);
            ClSize = (monocluster->amplitudes()).size();
 
            const auto& amplitudes = monocluster->amplitudes();
 
            barycenter = monocluster->barycenter() - 0.5;
            uint16_t FirstStrip = monocluster->firstStrip();
            auto begin = amplitudes.begin();
            nstrip = 0;
            for (auto idigi = begin; idigi != amplitudes.end(); idigi++) {
              Amplitudes[nstrip] = *idigi;
              sumx += pow(((FirstStrip + idigi - begin) - barycenter), 2) * (*idigi);
              HitCharge += *idigi;
            }
            hit_std_dev = sqrt(sumx / HitCharge);
 
            XGlobal = monogposition.x();
            YGlobal = monogposition.y();
            ZGlobal = monogposition.z();
 
            Type = detid.subdetId();
            MonoStereo = detid.stereo();
 
            if (detid.subdetId() == int(StripSubdetector::TIB)) {
              Layer = tTopo->tibLayer(detid);
              bw_fw = tTopo->tibStringInfo(detid)[0];
              Ext_Int = tTopo->tibStringInfo(detid)[1];
            }
            if (detid.subdetId() == int(StripSubdetector::TOB)) {
              Layer = tTopo->tobLayer(detid);
              bw_fw = tTopo->tobRodInfo(detid)[0];
            }
            if (detid.subdetId() == int(StripSubdetector::TID)) {
              Wheel = tTopo->tidWheel(detid);
              bw_fw = tTopo->tidModuleInfo(detid)[0];
            }
            if (detid.subdetId() == int(StripSubdetector::TEC)) {
              Wheel = tTopo->tecWheel(detid);
              bw_fw = tTopo->tecPetalInfo(detid)[0];
            }
 
            LocalVector monotkdir = monodet->toLocal(gtrkdir);
 
            if (monotkdir.z() != 0) {
              // THE LOCAL ANGLE (MONO)
              float tanangle = monotkdir.x() / monotkdir.z();
              TanTrackAngleParallel = monotkdir.y() / monotkdir.z();
              TanTrackAngle = tanangle;
              detparmap::iterator TheDet = detmap.find(detid.rawId());
              LocalVector localmagdir;
              if (TheDet != detmap.end())
                localmagdir = TheDet->second->magfield;
              MagField = localmagdir.mag();
              if (MagField != 0.) {
                LocalVector monoylocal(0, 1, 0);
                float signcorrection = (localmagdir * monoylocal) / (MagField);
                if (signcorrection != 0)
                  SignCorrection = 1 / signcorrection;
              }
 
              std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator alreadystored =
                  hitangleassociation.find(monohit);
 
              if (alreadystored != hitangleassociation.end()) {  //decide which hit take
                if (itm->estimate() > alreadystored->second.first) {
                  worse_double_hit++;
                }
                if (itm->estimate() < alreadystored->second.first) {
                  better_double_hit++;
                  hitangleassociation.insert(std::make_pair(monohit, std::make_pair(itm->estimate(), tanangle)));
                }
              } else {
                hitangleassociation.insert(make_pair(monohit, std::make_pair(itm->estimate(), tanangle)));
                HitsTree->Fill();
                hitcounter++;
              }
 
              // THIS THE POINTER TO THE STEREO HIT OF A MATCHED HIT
 
              // top be migrated to the more direct interface of matchedhit
              cache.push_back(matchedhit->stereoHit());
              const SiStripRecHit2D* stereohit = &cache.back();
              const SiStripRecHit2D::ClusterRef& stereocluster = stereohit->cluster();
              const GeomDetUnit* stereodet = gdet->stereoDet();
              // this does not exists anymore! either project the matched or use CPE
              const LocalPoint stereoposition = stereohit->localPosition();
              StripSubdetector detid = (StripSubdetector)stereohit->geographicalId();
              id_detector = detid.rawId();
              thick_detector = stereodet->specificSurface().bounds().thickness();
              const StripTopology& stopol = (StripTopology&)stereodet->topology();
              pitch_detector = stopol.localPitch(stereoposition);
              const GlobalPoint stereogposition = (stereodet->surface()).toGlobal(stereoposition);
 
              ClSize = (stereocluster->amplitudes()).size();
 
              const auto& amplitudes = stereocluster->amplitudes();
 
              barycenter = stereocluster->barycenter() - 0.5;
              uint16_t FirstStrip = stereocluster->firstStrip();
              auto begin = amplitudes.begin();
              nstrip = 0;
              for (auto idigi = begin; idigi != amplitudes.end(); idigi++) {
                Amplitudes[nstrip] = *idigi;
                sumx += pow(((FirstStrip + idigi - begin) - barycenter), 2) * (*idigi);
                HitCharge += *idigi;
              }
              hit_std_dev = sqrt(sumx / HitCharge);
 
              XGlobal = stereogposition.x();
              YGlobal = stereogposition.y();
              ZGlobal = stereogposition.z();
 
              Type = detid.subdetId();
              MonoStereo = detid.stereo();
 
              if (detid.subdetId() == int(StripSubdetector::TIB)) {
                Layer = tTopo->tibLayer(detid);
                bw_fw = tTopo->tibStringInfo(detid)[0];
                Ext_Int = tTopo->tibStringInfo(detid)[1];
              }
              if (detid.subdetId() == int(StripSubdetector::TOB)) {
                Layer = tTopo->tobLayer(detid);
                bw_fw = tTopo->tobRodInfo(detid)[0];
              }
              if (detid.subdetId() == int(StripSubdetector::TID)) {
                Wheel = tTopo->tidWheel(detid);
                bw_fw = tTopo->tidModuleInfo(detid)[0];
              }
              if (detid.subdetId() == int(StripSubdetector::TEC)) {
                Wheel = tTopo->tecWheel(detid);
                bw_fw = tTopo->tecPetalInfo(detid)[0];
              }
 
              LocalVector stereotkdir = stereodet->toLocal(gtrkdir);
 
              if (stereotkdir.z() != 0) {
                // THE LOCAL ANGLE (STEREO)
                float tanangle = stereotkdir.x() / stereotkdir.z();
                TanTrackAngleParallel = stereotkdir.y() / stereotkdir.z();
                TanTrackAngle = tanangle;
                detparmap::iterator TheDet = detmap.find(detid.rawId());
                LocalVector localmagdir;
                if (TheDet != detmap.end())
                  localmagdir = TheDet->second->magfield;
                MagField = localmagdir.mag();
                LocalVector stereoylocal(0, 1, 0);
                if (MagField != 0.) {
                  float signcorrection = (localmagdir * stereoylocal) / (MagField);
                  if (signcorrection != 0)
                    SignCorrection = 1 / signcorrection;
                }
 
                std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator alreadystored =
                    hitangleassociation.find(stereohit);
 
                if (alreadystored != hitangleassociation.end()) {  //decide which hit take
                  if (itm->estimate() > alreadystored->second.first) {
                    worse_double_hit++;
                  }
                  if (itm->estimate() < alreadystored->second.first) {
                    better_double_hit++;
                    hitangleassociation.insert(std::make_pair(stereohit, std::make_pair(itm->estimate(), tanangle)));
                  }
                } else {
                  hitangleassociation.insert(std::make_pair(stereohit, std::make_pair(itm->estimate(), tanangle)));
                  HitsTree->Fill();
                  hitcounter++;
                }
              }
            }
          } else if (hit) {
            //  hit= POINTER TO THE RECHIT
 
            const SiStripRecHit2D::ClusterRef& cluster = hit->cluster();
 
            GeomDetUnit* gdet = (GeomDetUnit*)tkGeom_->idToDet(hit->geographicalId());
            const LocalPoint position = hit->localPosition();
            StripSubdetector detid = (StripSubdetector)hit->geographicalId();
            id_detector = detid.rawId();
            thick_detector = gdet->specificSurface().bounds().thickness();
            const StripTopology& topol = (StripTopology&)gdet->topology();
            pitch_detector = topol.localPitch(position);
            const GlobalPoint gposition = (gdet->surface()).toGlobal(position);
 
            ClSize = (cluster->amplitudes()).size();
 
            const auto& amplitudes = cluster->amplitudes();
 
            barycenter = cluster->barycenter() - 0.5;
            uint16_t FirstStrip = cluster->firstStrip();
            nstrip = 0;
            auto begin = amplitudes.begin();
            for (auto idigi = amplitudes.begin(); idigi != amplitudes.end(); idigi++) {
              Amplitudes[nstrip] = *idigi;
              sumx += pow(((FirstStrip + idigi - begin) - barycenter), 2) * (*idigi);
              HitCharge += *idigi;
            }
            hit_std_dev = sqrt(sumx / HitCharge);
 
            XGlobal = gposition.x();
            YGlobal = gposition.y();
            ZGlobal = gposition.z();
 
            Type = detid.subdetId();
            MonoStereo = detid.stereo();
 
            if (detid.subdetId() == int(StripSubdetector::TIB)) {
              Layer = tTopo->tibLayer(detid);
              bw_fw = tTopo->tibStringInfo(detid)[0];
              Ext_Int = tTopo->tibStringInfo(detid)[1];
            }
            if (detid.subdetId() == int(StripSubdetector::TOB)) {
              Layer = tTopo->tobLayer(detid);
              bw_fw = tTopo->tobRodInfo(detid)[0];
            }
            if (detid.subdetId() == int(StripSubdetector::TID)) {
              Wheel = tTopo->tidWheel(detid);
              bw_fw = tTopo->tidModuleInfo(detid)[0];
            }
            if (detid.subdetId() == int(StripSubdetector::TEC)) {
              Wheel = tTopo->tecWheel(detid);
              bw_fw = tTopo->tecPetalInfo(detid)[0];
            }
 
            if (trackdirection.z() != 0) {
              // THE LOCAL ANGLE
              float tanangle = trackdirection.x() / trackdirection.z();
              TanTrackAngleParallel = trackdirection.y() / trackdirection.z();
              TanTrackAngle = tanangle;
              detparmap::iterator TheDet = detmap.find(detid.rawId());
              LocalVector localmagdir;
              if (TheDet != detmap.end())
                localmagdir = TheDet->second->magfield;
              MagField = localmagdir.mag();
              if (MagField != 0.) {
                LocalVector ylocal(0, 1, 0);
                float signcorrection = (localmagdir * ylocal) / (MagField);
                if (signcorrection != 0)
                  SignCorrection = 1 / signcorrection;
              }
 
              std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator alreadystored =
                  hitangleassociation.find(hit);
 
              if (alreadystored != hitangleassociation.end()) {  //decide which hit take
                if (itm->estimate() > alreadystored->second.first) {
                  worse_double_hit++;
                }
                if (itm->estimate() < alreadystored->second.first) {
                  better_double_hit++;
                  hitangleassociation.insert(std::make_pair(hit, std::make_pair(itm->estimate(), tanangle)));
                }
              } else {
                hitangleassociation.insert(std::make_pair(hit, std::make_pair(itm->estimate(), tanangle)));
                HitsTree->Fill();
                hitcounter++;
              }
            }
          }
        }
      }
    }
  }
  std::map<const SiStripRecHit2D*, std::pair<float, float>, DetIdLess>::iterator hitsiter;
 
  for (hitsiter = hitangleassociation.begin(); hitsiter != hitangleassociation.end(); hitsiter++) {
    hitcounter_2ndloop++;
 
    const SiStripRecHit2D* hit = hitsiter->first;
    const SiStripRecHit2D::ClusterRef& cluster = hit->cluster();
 
    size = (cluster->amplitudes()).size();
 
    StripSubdetector detid = (StripSubdetector)hit->geographicalId();
 
    float tangent = hitsiter->second.second;
 
    //Sign and XZ plane projection correction applied in TrackLocalAngle (TIB|TOB layers)
 
    detparmap::iterator thedet = detmap.find(detid.rawId());
    LocalVector localmagdir;
    if (thedet != detmap.end())
      localmagdir = thedet->second->magfield;
    float localmagfield = localmagdir.mag();
 
    if (localmagfield != 0.) {
      LocalVector ylocal(0, 1, 0);
 
      float normprojection = (localmagdir * ylocal) / (localmagfield);
 
      if (normprojection == 0.)
        LogDebug("SiStripLAProfileBooker::analyze") << "Error: YBprojection = 0";
 
      else {
        float signprojcorrection = 1 / normprojection;
        tangent *= signprojcorrection;
      }
    }
 
    //Filling histograms
 
    histomap::iterator thehisto = histos.find(detid.rawId());
 
    if (thehisto == histos.end())
      edm::LogError("SiStripLAProfileBooker::analyze")
          << "Error: the profile associated to" << detid.rawId() << "does not exist! ";
    else
      thehisto->second->Fill(tangent, size);
 
    //Summary histograms
    std::string name;
    unsigned int layerid;
    getlayer(detid, tTopo, name, layerid);
    histomap::iterator thesummaryhisto = summaryhisto.find(layerid);
    if (thesummaryhisto == summaryhisto.end())
      edm::LogError("SiStripLAProfileBooker::analyze")
          << "Error: the profile associated to subdet " << name << "does not exist! ";
    else
      thesummaryhisto->second->Fill(tangent, size);
  }
}
 
//Makename function
 
void SiStripLAProfileBooker::getlayer(const DetId& detid,
                                      const TrackerTopology* tTopo,
                                      std::string& name,
                                      unsigned int& layerid) {
  int layer = 0;
  std::stringstream layernum;
 
  if (detid.subdetId() == int(StripSubdetector::TIB)) {
    name += "TIB_Layer_";
    layer = tTopo->tibLayer(detid);
  }
 
  else if (detid.subdetId() == int(StripSubdetector::TID)) {
    name += "TID_Ring_";
    layer = tTopo->tidRing(detid);
  }
 
  else if (detid.subdetId() == int(StripSubdetector::TOB)) {
    name += "TOB_Layer_";
    layer = tTopo->tobLayer(detid);
 
  }
 
  else if (detid.subdetId() == int(StripSubdetector::TEC)) {
    name += "TEC_Ring_";
    layer = tTopo->tecRing(detid);
  }
  layernum << layer;
  name += layernum.str();
  layerid = detid.subdetId() * 10 + layer;
}
 
void SiStripLAProfileBooker::endJob() {
  std::string outputFile_ = conf_.getUntrackedParameter<std::string>("fileName", "LorentzAngle.root");
  dbe_->save(outputFile_);
 
  hFile->Write();
  hFile->Close();
}