AlignmentStats.cc

Go to the documentation of this file.

#include "Alignment/CommonAlignmentMonitor/plugins/AlignmentStats.h"

#include "DataFormats/Common/interface/View.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
#include "Alignment/CommonAlignment/interface/Alignable.h"
#include "Alignment/CommonAlignment/interface/Utilities.h"

#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h"

#include "DataFormats/Alignment/interface/AlignmentClusterFlag.h"
#include "DataFormats/Alignment/interface/AliClusterValueMap.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
#include "Utilities/General/interface/ClassName.h"

using namespace std;

AlignmentStats::AlignmentStats(const edm::ParameterSet &iConfig) :
  src_(iConfig.getParameter<edm::InputTag>("src")),
  overlapAM_(iConfig.getParameter<edm::InputTag>("OverlapAssoMap")),
  keepTrackStats_(iConfig.getParameter<bool>("keepTrackStats")),
  keepHitPopulation_(iConfig.getParameter<bool>("keepHitStats")),
  statsTreeName_(iConfig.getParameter<string>("TrkStatsFileName")),
  hitsTreeName_(iConfig.getParameter<string>("HitStatsFileName")),
  prescale_(iConfig.getParameter<uint32_t>("TrkStatsPrescale")),
  lastSetup_(nullptr)
{

  //sanity checks

  //init
  outtree_=0;

}//end constructor

AlignmentStats::~AlignmentStats(){
                   //
}

void AlignmentStats::beginJob(){// const edm::EventSetup &iSetup

  //book track stats tree
  treefile_=new TFile(statsTreeName_.c_str(),"RECREATE");
  treefile_->cd();
  outtree_=new TTree("AlignmentTrackStats","Statistics of Tracks used for Alignment");
  // int nHitsinPXB[MAXTRKS_], nHitsinPXE[MAXTRKS_], nHitsinTEC[MAXTRKS_], nHitsinTIB[MAXTRKS_],nHitsinTOB[MAXTRKS_],nHitsinTID[MAXTRKS_];
  

  outtree_->Branch("Ntracks"    ,&ntracks ,"Ntracks/i");
  outtree_->Branch("Run_"        ,&run_     ,"Run_Nr/I");
  outtree_->Branch("Event"      ,&event_   ,"EventNr/I");
  outtree_->Branch("Eta"        ,&Eta     ,"Eta[Ntracks]/F");
  outtree_->Branch("Phi"        ,&Phi     ,"Phi[Ntracks]/F");
  outtree_->Branch("P"          ,&P       ,"P[Ntracks]/F");
  outtree_->Branch("Pt"         ,&Pt      ,"Pt[Ntracks]/F");
  outtree_->Branch("Chi2n"      ,&Chi2n   ,"Chi2n[Ntracks]/F");
  outtree_->Branch("Nhits"      ,&Nhits   ,"Nhits[Ntracks][7]/I");
  /*
    outtree_->Branch("NhitsPXB"       , ,);
    outtree_->Branch("NhitsPXE"       , ,);
    outtree_->Branch("NhitsTIB"       , ,);
    outtree_->Branch("NhitsTID"       , ,);
    outtree_->Branch("NhitsTOB"       , ,);
    outtree_->Branch("NhitsTOB"       , ,);
  */

  tmpPresc_=prescale_;
  firstEvent_=true;

  //load the tracker geometry from the EventSetup
  //  iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry_);
 
}//end beginJob


void AlignmentStats::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup){

  lastSetup_ = &iSetup;

  //load list of detunits needed then in endJob
 
  if(firstEvent_){
    edm::ESHandle<TrackerGeometry> tmpTkGeometry;
    iSetup.get<TrackerDigiGeometryRecord>().get(tmpTkGeometry);
    trackerGeometry_=&(*tmpTkGeometry);
    firstEvent_=false;
  }


  //take trajectories and tracks to loop on
  // edm::Handle<TrajTrackAssociationCollection> TrackAssoMap;
  edm::Handle<reco::TrackCollection> Tracks;
  iEvent.getByLabel(src_, Tracks);

  //take overlap HitAssomap
  edm::Handle<AliClusterValueMap> hMap;
  iEvent.getByLabel(overlapAM_, hMap);
  const AliClusterValueMap &OverlapMap=*hMap;

  // Initialise
  run_=1;
  event_=1;
  ntracks=0;
  run_=iEvent.id().run();
  event_=iEvent.id().event();
  ntracks=Tracks->size();
  //  if(ntracks>1) std::cout<<"~~~~~~~~~~~~\n For this event processing "<<ntracks<<" tracks"<<std::endl;

  unsigned int trk_cnt=0;
  
  for(int j=0;j<MAXTRKS_;j++){
    Eta[j]=-9999.0;
    Phi[j]=-8888.0;
    P[j]  =-7777.0;
    Pt[j] =-6666.0;
    Chi2n[j]=-2222.0;
    for(int k=0;k<7;k++){
      Nhits[j][k]=0;
    }
  }
  
  // int npxbhits=0;

  //loop on tracks
  for(std::vector<reco::Track>::const_iterator ittrk = Tracks->begin(), edtrk = Tracks->end(); ittrk != edtrk; ++ittrk){
    Eta[trk_cnt]=ittrk->eta();
    Phi[trk_cnt]=ittrk->phi();
    Chi2n[trk_cnt]=ittrk->normalizedChi2();
    P[trk_cnt]=ittrk->p();
    Pt[trk_cnt]=ittrk->pt();
    Nhits[trk_cnt][0]=ittrk->numberOfValidHits();
 
    //   if(ntracks>1)std::cout<<"Track #"<<trk_cnt+1<<" params:    Eta="<< Eta[trk_cnt]<<"  Phi="<< Phi[trk_cnt]<<"  P="<<P[trk_cnt]<<"   Nhits="<<Nhits[trk_cnt][0]<<std::endl;
   
    int nhit=0;
    //loop on tracking rechits
    //std::cout << "   loop on hits of track #" << (itt - tracks->begin()) << std::endl;
    for (trackingRecHit_iterator ith = ittrk->recHitsBegin(), edh = ittrk->recHitsEnd(); ith != edh; ++ith) {

      const TrackingRecHit *hit = *ith; // ith is an iterator on edm::Ref to rechit
      if(! hit->isValid())continue;
      DetId detid = hit->geographicalId();
      int subDet = detid.subdetId();
      uint32_t rawId = hit->geographicalId().rawId();

      //  if(subDet==1)npxbhits++;

      //look if you find this detid in the map
      DetHitMap::iterator mapiter;
      mapiter= hitmap_.find(rawId);
      if(mapiter!=hitmap_.end() ){//present, increase its value by one
    //  hitmap_[rawId]=hitmap_[rawId]+1;
    ++(hitmap_[rawId]);
      }
      else{//not present, let's add this key to the map with value=1
    hitmap_.insert(pair<uint32_t, uint32_t>(rawId, 1));
      }
      
      AlignmentClusterFlag inval;

      bool hitInPixel= (subDet==PixelSubdetector::PixelBarrel)||(subDet==PixelSubdetector::PixelEndcap);
      bool hitInStrip=(subDet==SiStripDetId::TIB) || (subDet==SiStripDetId::TID) ||(subDet==SiStripDetId::TOB) ||(subDet==SiStripDetId::TEC);

      if(!(hitInPixel||hitInStrip)){
    //skip only this hit, don't stop everything throwing an exception
    edm::LogError("AlignmentStats")<<"Hit not belonging neither to pixels nor to strips! Skipping it. SubDet= "<<subDet;
    continue;
      }

      //check also if the hit is an overlap. If yes fill a dedicated hitmap
      if (hitInStrip){
    const std::type_info &type = typeid(*hit);

    if (type == typeid(SiStripRecHit1D)) {
      //Notice the difference respect to when one loops on Trajectories: the recHit is a TrackingRecHit and not a TransientTrackingRecHit
      const SiStripRecHit1D* striphit=dynamic_cast<const  SiStripRecHit1D*>(hit);
      if(striphit!=0){
        SiStripRecHit1D::ClusterRef stripclust(striphit->cluster());
        inval = OverlapMap[stripclust];
      }
      else{
        //  edm::LogError("AlignmentStats")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit1D failed!   TypeId of the RecHit: "<<className(*hit);
        throw cms::Exception("NullPointerError")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit1D failed!   TypeId of the RecHit: "<<className(*hit);
      }
    }//end if sistriprechit1D
    if (type == typeid(SiStripRecHit2D)) {
      const SiStripRecHit2D* striphit=dynamic_cast<const  SiStripRecHit2D*>(hit);
      if(striphit!=0){
        SiStripRecHit2D::ClusterRef stripclust(striphit->cluster());
        inval = OverlapMap[stripclust];
        //cout<<"Taken the Strip Cluster with ProdId "<<stripclust.id() <<"; the Value in the map is "<<inval<<"  (DetId is "<<hit->geographicalId().rawId()<<")"<<endl;
      }
      else{
        throw cms::Exception("NullPointerError")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit2D failed!   TypeId of the RecHit: "<<className(*hit);
        //    edm::LogError("AlignmentStats")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Strip RecHit2D failed!   TypeId of the RecHit: "<<className(*hit);
      }
    }//end if sistriprechit2D

      }//end if hit in Strips
      else{
    const SiPixelRecHit* pixelhit= dynamic_cast<const SiPixelRecHit*>(hit);
    if(pixelhit!=0){
      SiPixelClusterRefNew pixclust(pixelhit->cluster());
      inval = OverlapMap[pixclust];
      //cout<<"Taken the Pixel Cluster with ProdId "<<pixclust.id() <<"; the Value in the map is "<<inval<<"  (DetId is "<<hit->geographicalId().rawId()<<")"<<endl;
    }
    else{
      edm::LogError("AlignmentStats")<<"ERROR in <AlignmentStats::analyze>: Dynamic cast of Pixel RecHit failed!   TypeId of the RecHit: "<<className(*hit);
    }
      }//end else hit is in Pixel


      bool isOverlapHit(inval.isOverlap());
      
      if( isOverlapHit ){
    //cout<<"This hit is an overlap !"<<endl;
    DetHitMap::iterator overlapiter;
    overlapiter=overlapmap_.find(rawId);
    
    if(overlapiter!=overlapmap_.end() ){//the det already collected at least an overlap, increase its value by one
        overlapmap_[rawId]=overlapmap_[rawId]+1;
    }
    else{//first overlap on det unit, let's add it to the map
      overlapmap_.insert(pair<uint32_t, uint32_t>(rawId, 1));
    }
      }//end if the hit is an overlap
      



      int subdethit= static_cast<int>(hit->geographicalId().subdetId());
      // if(ntracks>1)std::cout<<"Hit in SubDet="<<subdethit<<std::endl;
      Nhits[trk_cnt][subdethit]=Nhits[trk_cnt][subdethit]+1;
      nhit++;
    }//end loop on trackingrechits
    trk_cnt++;

  }//end loop on tracks

  //  cout<<"Total number of pixel hits is "<<npxbhits<<endl;

  tmpPresc_--;
  if(tmpPresc_<1){
    outtree_->Fill();
    tmpPresc_=prescale_;
  }
  if(trk_cnt!=ntracks)edm::LogError("AlignmentStats")<<"\nERROR! trk_cnt="<<trk_cnt<<"   ntracks="<<ntracks;
  trk_cnt=0;

  return;
}

void AlignmentStats::endJob(){

  treefile_->cd();
  edm::LogInfo("AlignmentStats")<<"Writing out the TrackStatistics in "<<gDirectory->GetPath()<<std::endl;
  outtree_->Write();
  delete outtree_;

  //create tree with hit maps (hitstree)
  //book track stats tree
  TFile *hitsfile=new TFile(hitsTreeName_.c_str(),"RECREATE");
  hitsfile->cd();
  TTree *hitstree=new TTree("AlignmentHitMap","Maps of Hits used for Alignment");

  unsigned int id=0,nhits=0,noverlaps=0;
  float posX(-99999.0),posY(-77777.0),posZ(-88888.0);
  float posEta(-6666.0),posPhi(-5555.0),posR(-4444.0);
  int subdet=0;
  unsigned int layer=0; 
  bool is2D=false,isStereo=false;
  hitstree->Branch("DetId",    &id ,      "DetId/i");
  hitstree->Branch("Nhits",    &nhits ,   "Nhits/i");
  hitstree->Branch("Noverlaps",&noverlaps,"Noverlaps/i");
  hitstree->Branch("SubDet",   &subdet,   "SubDet/I");
  hitstree->Branch("Layer",    &layer,    "Layer/i");
  hitstree->Branch("is2D" ,    &is2D,     "is2D/B");
  hitstree->Branch("isStereo", &isStereo, "isStereo/B");
  hitstree->Branch("posX",     &posX,     "posX/F");
  hitstree->Branch("posY",     &posY,     "posY/F");
  hitstree->Branch("posZ",     &posZ,     "posZ/F");
  hitstree->Branch("posR",     &posR,     "posR/F");
  hitstree->Branch("posEta",   &posEta,   "posEta/F");
  hitstree->Branch("posPhi",   &posPhi,   "posPhi/F");

  /*
  TTree *overlapstree=new TTree("OverlapHitMap","Maps of Overlaps used for Alignment");
  hitstree->Branch("DetId",   &id ,     "DetId/i");
  hitstree->Branch("NOverlaps",   &nhits ,  "Nhits/i");
  hitstree->Branch("SubDet",  &subdet,  "SubDet/I");
  hitstree->Branch("Layer",   &layer,   "Layer/i");
  hitstree->Branch("is2D" ,   &is2D,    "is2D/B");
  hitstree->Branch("isStereo",&isStereo,"isStereo/B");
  hitstree->Branch("posX",    &posX,    "posX/F");
  hitstree->Branch("posY",    &posY,    "posY/F");
  hitstree->Branch("posZ",    &posZ,    "posZ/F");
  hitstree->Branch("posR",    &posR,    "posR/F");
  hitstree->Branch("posEta",  &posEta,  "posEta/F");
  hitstree->Branch("posPhi",  &posPhi,  "posPhi/F");
  */

  //Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHandle;
  lastSetup_->get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology* const tTopo = tTopoHandle.product();

  AlignableTracker* theAliTracker=new AlignableTracker(&(*trackerGeometry_), tTopo);
  const std::vector<Alignable*>& Detunitslist=theAliTracker->deepComponents();
  int ndetunits=Detunitslist.size();
  edm::LogInfo("AlignmentStats")<<"Number of DetUnits in the AlignableTracker: "<< ndetunits<<std::endl;

  for(int det_cnt=0;det_cnt< ndetunits;++det_cnt){
    
    //re-initialize for safety
    id=0;
    nhits=0;
    noverlaps=0;
    posX=-99999.0;
    posY=-77777.0;
    posZ=-88888.0;
    posEta=-6666.0;
    posPhi=-5555.0;
    posR=-4444.0;
    subdet=0;
    layer=0; 
    is2D=false;
    isStereo=false;

    //if detunit in vector is found also in the map, look for how many hits were collected 
    //and save in the tree this number
    id=static_cast <uint32_t>( Detunitslist[det_cnt]->id() );
    if( hitmap_.find(id) != hitmap_.end() ){
      nhits=hitmap_[id];
    }
    //if not, save nhits=0
    else{
      nhits=0;
      hitmap_.insert(pair<uint32_t, uint32_t>(id, 0));
    }

    if( overlapmap_.find(id) != overlapmap_.end() ){
      noverlaps=overlapmap_[id];
    }
    //if not, save nhits=0
    else{
      noverlaps=0;
      overlapmap_.insert(pair<uint32_t, uint32_t>(id, 0));
    }

  //take other geometrical infos from the det
   posX= Detunitslist[det_cnt]->globalPosition().x();
   posY= Detunitslist[det_cnt]->globalPosition().y();
   posZ= Detunitslist[det_cnt]->globalPosition().z();

   align::GlobalVector vec(posX,posY,posZ);
   posR = vec.perp();
   posPhi = vec.phi();
   posEta = vec.eta();
   //   posPhi = atan2(posY,posX);  

   DetId detid(id);
   subdet=detid.subdetId();

   //get layers, petals, etc...
   if(subdet==PixelSubdetector::PixelBarrel){//PXB
     
     layer=tTopo->pxbLayer(id);
     is2D=true;
     isStereo=false;
   }
   else if(subdet==PixelSubdetector::PixelEndcap){
     
     layer=tTopo->pxfDisk(id);
     is2D=true;
     isStereo=false;
   }
   else if(subdet==SiStripDetId::TIB){
     
     layer=tTopo->tibLayer(id);
     is2D=tTopo->tibIsDoubleSide(id);
     isStereo=tTopo->tibIsStereo(id);
   }
   else if(subdet==SiStripDetId::TID){
     
     layer=tTopo->tidWheel(id);
     is2D=tTopo->tidIsDoubleSide(id);
     isStereo=tTopo->tidIsStereo(id);
   }
   else if(subdet==SiStripDetId::TOB){
     
     layer=tTopo->tobLayer(id);
     is2D=tTopo->tobIsDoubleSide(id);
     isStereo=tTopo->tobIsStereo(id);
   }
   else if(subdet==SiStripDetId::TEC){
     
     layer=tTopo->tecWheel(id);
     is2D=tTopo->tecIsDoubleSide(id);
     isStereo=tTopo->tecIsStereo(id);
   }
   else{
     edm::LogError("AlignmentStats")<<"Detector not belonging neither to pixels nor to strips! Skipping it. SubDet= "<<subdet;
   }

  //write in the hitstree
    hitstree->Fill();
  }//end loop over detunits


  //save hitstree
  hitstree->Write();
  delete hitstree;
  //delete Detunitslist;
  hitmap_.clear();
  overlapmap_.clear();
  delete hitsfile;
}
// ========= MODULE DEF ==============
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(AlignmentStats);