DeDxDiscriminatorLearner.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    DeDxDiscriminatorLearner
// Class:      DeDxDiscriminatorLearner
// 
//
// Original Author:  Loic Quertenmont(querten)
//         Created:  Mon October 20 10:09:02 CEST 2008
//


// system include files
#include <memory>

#include "CalibTracker/SiStripChannelGain/plugins/DeDxDiscriminatorLearner.h"

#include "DataFormats/TrackReco/interface/Track.h"
#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"

using namespace reco;
using namespace std;
using namespace edm;

DeDxDiscriminatorLearner::DeDxDiscriminatorLearner(const edm::ParameterSet& iConfig) : ConditionDBWriter<PhysicsTools::Calibration::HistogramD3D>(iConfig)
{
   m_tracksTag                 = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
   m_trajTrackAssociationTag   = consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation"));

   P_Min           = iConfig.getParameter<double>  ("P_Min"  );
   P_Max               = iConfig.getParameter<double>  ("P_Max"  );
   P_NBins             = iConfig.getParameter<int>     ("P_NBins");
   Path_Min            = iConfig.getParameter<double>  ("Path_Min"  );
   Path_Max            = iConfig.getParameter<double>  ("Path_Max"  );
   Path_NBins          = iConfig.getParameter<int>     ("Path_NBins");
   Charge_Min          = iConfig.getParameter<double>  ("Charge_Min"  );
   Charge_Max          = iConfig.getParameter<double>  ("Charge_Max"  );
   Charge_NBins        = iConfig.getParameter<int>     ("Charge_NBins");

   MinTrackMomentum    = iConfig.getUntrackedParameter<double>  ("minTrackMomentum"   ,  5.0);
   MaxTrackMomentum    = iConfig.getUntrackedParameter<double>  ("maxTrackMomentum"   ,  99999.0); 
   MinTrackEta         = iConfig.getUntrackedParameter<double>  ("minTrackEta"        , -5.0);
   MaxTrackEta         = iConfig.getUntrackedParameter<double>  ("maxTrackEta"        ,  5.0);
   MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  255);
   MinTrackHits        = iConfig.getUntrackedParameter<unsigned>("MinTrackHits"       ,  4);

   algoMode            = iConfig.getUntrackedParameter<string>  ("AlgoMode"           ,  "SingleJob");
   HistoFile           = iConfig.getUntrackedParameter<string>  ("HistoFile"        ,  "out.root");
   VInputFiles         = iConfig.getUntrackedParameter<vector<string> >        ("InputFiles");

   shapetest           = iConfig.getParameter<bool>("ShapeTest");
   useCalibration      = iConfig.getUntrackedParameter<bool>("UseCalibration");
   m_calibrationPath   = iConfig.getUntrackedParameter<string>("calibrationPath");
}


DeDxDiscriminatorLearner::~DeDxDiscriminatorLearner(){}

// ------------ method called once each job just before starting event loop  ------------

void  DeDxDiscriminatorLearner::algoBeginJob(const edm::EventSetup& iSetup)
{
   Charge_Vs_Path = new TH3F ("Charge_Vs_Path"     , "Charge_Vs_Path" , P_NBins, P_Min, P_Max, Path_NBins, Path_Min, Path_Max, Charge_NBins, Charge_Min, Charge_Max);

   if(useCalibration && calibGains.size()==0){
      edm::ESHandle<TrackerGeometry> tkGeom;
      iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
      m_off = tkGeom->offsetDU(GeomDetEnumerators::PixelBarrel); //index start at the first pixel

      DeDxTools::makeCalibrationMap(m_calibrationPath, *tkGeom, calibGains, m_off);
   }

   //Read the calibTree if in calibTree mode
   if(strcmp(algoMode.c_str(),"CalibTree")==0)algoAnalyzeTheTree(iSetup);
}

// ------------ method called once each job just after ending the event loop  ------------

void DeDxDiscriminatorLearner::algoEndJob()
{
   if( strcmp(algoMode.c_str(),"MultiJob")==0){
    TFile* Output = new TFile(HistoFile.c_str(), "RECREATE");
        Charge_Vs_Path->Write();
    Output->Write();
    Output->Close();
   }else if( strcmp(algoMode.c_str(),"WriteOnDB")==0){
        TFile* Input = new TFile(HistoFile.c_str() );
    Charge_Vs_Path = (TH3F*)(Input->FindObjectAny("Charge_Vs_Path"))->Clone();  
    Input->Close();
   }else if(strcmp(algoMode.c_str(),"CalibTree")==0){
        TFile* Output = new TFile(HistoFile.c_str(), "RECREATE");
        Charge_Vs_Path->Write();
        Output->Write();
        Output->Close();
        TFile* Input = new TFile(HistoFile.c_str() );
        Charge_Vs_Path = (TH3F*)(Input->FindObjectAny("Charge_Vs_Path"))->Clone();
        Input->Close();
   }
}

void DeDxDiscriminatorLearner::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{   
   Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
   iEvent.getByToken(m_trajTrackAssociationTag, trajTrackAssociationHandle);

   edm::Handle<reco::TrackCollection> trackCollectionHandle;
   iEvent.getByToken(m_tracksTag,trackCollectionHandle);

   unsigned track_index = 0;
   for(TrajTrackAssociationCollection::const_iterator it = trajTrackAssociationHandle->begin(); it!=trajTrackAssociationHandle->end(); ++it, track_index++) {
      const Track&      track = *it->val;
      const Trajectory& traj  = *it->key;

      if(track.eta()  <MinTrackEta      || track.eta()>MaxTrackEta     ){continue;}
      if(track.pt()   <MinTrackMomentum || track.pt() >MaxTrackMomentum){continue;}
      if(track.found()<MinTrackHits                                    ){continue;}

      const vector<TrajectoryMeasurement> & measurements = traj.measurements();
      for(vector<TrajectoryMeasurement>::const_iterator it = measurements.begin(); it!=measurements.end(); it++){
         TrajectoryStateOnSurface trajState=it->updatedState();
         if( !trajState.isValid()) continue;


         const TrackingRecHit * recHit=(*it->recHit()).hit();
         if(!recHit || !recHit->isValid())continue;
         LocalVector trackDirection = trajState.localDirection();
         float cosine = trackDirection.z()/trackDirection.mag();

         processHit(recHit, trajState.localMomentum().mag(), cosine, trajState);
      }
   }

}


void DeDxDiscriminatorLearner::processHit(const TrackingRecHit* recHit, float trackMomentum, float& cosine, const TrajectoryStateOnSurface& trajState){
      auto const & thit = static_cast<BaseTrackerRecHit const&>(*recHit);
      if(!thit.isValid())return;

      auto const & clus = thit.firstClusterRef();
      if(!clus.isValid())return;

      int   NSaturating = 0;
      if(clus.isPixel()){
          return;
       }else if(clus.isStrip() && !thit.isMatched()){
          auto& detUnit     = *(recHit->detUnit());
          auto& cluster     = clus.stripCluster();
          if( cluster.amplitudes().size()>MaxNrStrips)                                            { return; }
          if( DeDxTools::IsSpanningOver2APV (cluster.firstStrip(), cluster.amplitudes().size()))  { return; }
          if(!DeDxTools::IsFarFromBorder    (trajState, &detUnit ))                               { return; }
          float pathLen     = 10.0*detUnit.surface().bounds().thickness()/fabs(cosine);
          float chargeAbs   = DeDxTools::getCharge(&cluster,NSaturating, detUnit, calibGains, m_off);
          float charge      = chargeAbs/pathLen;
          if(!shapetest || (shapetest && DeDxTools::shapeSelection(cluster))){
             Charge_Vs_Path->Fill(trackMomentum, pathLen, charge);
          }
       }else if(clus.isStrip() && thit.isMatched()){
          const SiStripMatchedRecHit2D* matchedHit=dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
          if(!matchedHit)return;
          const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(matchedHit->det());

          auto& detUnitM    = *(gdet->monoDet());
          auto& clusterM    = matchedHit->monoCluster();
          if( clusterM.amplitudes().size()>MaxNrStrips)                                             { return; }
          if( DeDxTools::IsSpanningOver2APV (clusterM.firstStrip(), clusterM.amplitudes().size()))  { return; }
          if(!DeDxTools::IsFarFromBorder    (trajState, &detUnitM ))                                { return; }
          float pathLen     = 10.0*detUnitM.surface().bounds().thickness()/fabs(cosine);
          float chargeAbs   = DeDxTools::getCharge(&clusterM,NSaturating, detUnitM, calibGains, m_off);
          float charge      = chargeAbs/pathLen;
          if(!shapetest || (shapetest && DeDxTools::shapeSelection(clusterM))){
             Charge_Vs_Path->Fill(trackMomentum, pathLen, charge);
          }

          auto& detUnitS    = *(gdet->stereoDet());
          auto& clusterS    = matchedHit->stereoCluster();
          if( clusterS.amplitudes().size()>MaxNrStrips)                                             { return; }
          if( DeDxTools::IsSpanningOver2APV (clusterS.firstStrip(), clusterS.amplitudes().size()))  { return; }
          if(!DeDxTools::IsFarFromBorder    (trajState, &detUnitS ))                                { return; }
          pathLen     = 10.0*detUnitS.surface().bounds().thickness()/fabs(cosine);
          chargeAbs   = DeDxTools::getCharge(&clusterS,NSaturating, detUnitS, calibGains, m_off);
          charge      = chargeAbs/pathLen;
          if(!shapetest || (shapetest && DeDxTools::shapeSelection(clusterS))){
             Charge_Vs_Path->Fill(trackMomentum, pathLen, charge);
          }      
       }
}


//this function is only used when we run over a calibTree instead of running over EDM files
void DeDxDiscriminatorLearner::algoAnalyzeTheTree(const edm::EventSetup& iSetup)
{
   edm::ESHandle<TrackerGeometry> tkGeom;
   iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );

   unsigned int NEvent = 0;
   for(unsigned int i=0;i<VInputFiles.size();i++){
      printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
      TChain* tree = new TChain("gainCalibrationTree/tree");
      tree->Add(VInputFiles[i].c_str());

      TString EventPrefix("");
      TString EventSuffix("");

      TString TrackPrefix("track");
      TString TrackSuffix("");

      TString CalibPrefix("GainCalibration");
      TString CalibSuffix("");

      unsigned int                 eventnumber    = 0;    tree->SetBranchAddress(EventPrefix + "event"          + EventSuffix, &eventnumber   , NULL);
      unsigned int                 runnumber      = 0;    tree->SetBranchAddress(EventPrefix + "run"            + EventSuffix, &runnumber     , NULL);
      std::vector<bool>*           TrigTech       = 0;    tree->SetBranchAddress(EventPrefix + "TrigTech"       + EventSuffix, &TrigTech   , NULL);

      std::vector<double>*         trackchi2ndof  = 0;    tree->SetBranchAddress(TrackPrefix + "chi2ndof"       + TrackSuffix, &trackchi2ndof , NULL);
      std::vector<float>*          trackp         = 0;    tree->SetBranchAddress(TrackPrefix + "momentum"       + TrackSuffix, &trackp        , NULL);
      std::vector<float>*          trackpt        = 0;    tree->SetBranchAddress(TrackPrefix + "pt"             + TrackSuffix, &trackpt       , NULL);
      std::vector<double>*         tracketa       = 0;    tree->SetBranchAddress(TrackPrefix + "eta"            + TrackSuffix, &tracketa      , NULL);
      std::vector<double>*         trackphi       = 0;    tree->SetBranchAddress(TrackPrefix + "phi"            + TrackSuffix, &trackphi      , NULL);
      std::vector<unsigned int>*   trackhitsvalid = 0;    tree->SetBranchAddress(TrackPrefix + "hitsvalid"      + TrackSuffix, &trackhitsvalid, NULL);

      std::vector<int>*            trackindex     = 0;    tree->SetBranchAddress(CalibPrefix + "trackindex"     + CalibSuffix, &trackindex    , NULL);
      std::vector<unsigned int>*   rawid          = 0;    tree->SetBranchAddress(CalibPrefix + "rawid"          + CalibSuffix, &rawid         , NULL);
      std::vector<unsigned short>* firststrip     = 0;    tree->SetBranchAddress(CalibPrefix + "firststrip"     + CalibSuffix, &firststrip    , NULL);
      std::vector<unsigned short>* nstrips        = 0;    tree->SetBranchAddress(CalibPrefix + "nstrips"        + CalibSuffix, &nstrips       , NULL);
      std::vector<unsigned int>*   charge         = 0;    tree->SetBranchAddress(CalibPrefix + "charge"         + CalibSuffix, &charge        , NULL);
      std::vector<float>*          path           = 0;    tree->SetBranchAddress(CalibPrefix + "path"           + CalibSuffix, &path          , NULL);
      std::vector<unsigned char>*  amplitude      = 0;    tree->SetBranchAddress(CalibPrefix + "amplitude"      + CalibSuffix, &amplitude     , NULL);
      std::vector<double>*         gainused       = 0;    tree->SetBranchAddress(CalibPrefix + "gainused"       + CalibSuffix, &gainused      , NULL);

      printf("Number of Events = %i + %i = %i\n",NEvent,(unsigned int)tree->GetEntries(),(unsigned int)(NEvent+tree->GetEntries()));NEvent+=tree->GetEntries();
      printf("Progressing Bar              :0%%       20%%       40%%       60%%       80%%       100%%\n");
      printf("Looping on the Tree          :");
      int TreeStep = tree->GetEntries()/50;if(TreeStep<=1)TreeStep=1;
      for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
      if(ientry%TreeStep==0){printf(".");fflush(stdout);}
         tree->GetEntry(ientry);

         int FirstAmplitude = 0;
         for(unsigned int c=0;c<(*path).size();c++){
            FirstAmplitude+=(*nstrips)[c];
            int t = (*trackindex)[c];
            if((*trackpt)[t]<5)continue;
            if((*trackhitsvalid)[t]<5)continue;

            int Charge = 0; 
            if(useCalibration){
               auto & gains     = calibGains[tkGeom->idToDetUnit(DetId((*rawid)[c]))->index()-m_off];
               auto & gain      = gains[(*firststrip)[c]/128];
               for(unsigned int s=0;s<(*nstrips)[c];s++){
                 int StripCharge =  (*amplitude)[FirstAmplitude-(*nstrips)[c]+s];
                 if(StripCharge<254){
                    StripCharge=(int)(StripCharge/gain);
                    if(StripCharge>=1024){
                       StripCharge = 255;
                    }else if(StripCharge>=254){
                       StripCharge = 254;
                    }
                 }
                 Charge += StripCharge;
               }
            }else{
               Charge = (*charge)[c];
            }

//          printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[c],Charge,Gains[(*rawid)[c]]);
            double ClusterChargeOverPath   =  ( (double) Charge )/(*path)[c] ;       
            Charge_Vs_Path->Fill((*trackp)[t],(*path)[c],ClusterChargeOverPath);
         }
      }printf("\n");
  }
}


PhysicsTools::Calibration::HistogramD3D* DeDxDiscriminatorLearner::getNewObject()
{
   PhysicsTools::Calibration::HistogramD3D* obj;
   obj = new PhysicsTools::Calibration::HistogramD3D(
                Charge_Vs_Path->GetNbinsX(), Charge_Vs_Path->GetXaxis()->GetXmin(),  Charge_Vs_Path->GetXaxis()->GetXmax(),
                Charge_Vs_Path->GetNbinsY(), Charge_Vs_Path->GetYaxis()->GetXmin(),  Charge_Vs_Path->GetYaxis()->GetXmax(),
            Charge_Vs_Path->GetNbinsZ(), Charge_Vs_Path->GetZaxis()->GetXmin(),  Charge_Vs_Path->GetZaxis()->GetXmax());

   for(int ix=0; ix<=Charge_Vs_Path->GetNbinsX()+1; ix++){
      for(int iy=0; iy<=Charge_Vs_Path->GetNbinsY()+1; iy++){
         for(int iz=0; iz<=Charge_Vs_Path->GetNbinsZ()+1; iz++){
            obj->setBinContent(ix, iy, iz, Charge_Vs_Path->GetBinContent(ix,iy, iz) );       
//          if(Charge_Vs_Path->GetBinContent(ix,iy)!=0)printf("%i %i %i --> %f\n",ix,iy, iz, Charge_Vs_Path->GetBinContent(ix,iy,iz)); 
         }
      }
   }

   return obj;
}


//define this as a plug-in
DEFINE_FWK_MODULE(DeDxDiscriminatorLearner);