MuonDTDigis.cc

Go to the documentation of this file.

#include "MuonDTDigis.h"

#include <iostream>
#include <string>

#include "TFile.h"

#include "SimMuon/DTDigitizer/test/Histograms.h"

using namespace edm;
using namespace std;

MuonDTDigis::MuonDTDigis(const ParameterSet& pset){

  // ----------------------
  // Get the debug parameter for verbose output
  verbose_ = pset.getUntrackedParameter<bool>("verbose",false);

  // the name of the Digi collection
  SimHitToken_ = consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("SimHitLabel"));

  // the name of the Digi collection
  DigiToken_ = consumes<DTDigiCollection> (pset.getParameter<edm::InputTag>(("DigiLabel")));

  hDigis_global = std::make_unique<hDigis>("Global");
  hDigis_W0 = std::make_unique<hDigis>("Wheel0");
  hDigis_W1 = std::make_unique<hDigis>("Wheel1");
  hDigis_W2 = std::make_unique<hDigis>("Wheel2");
  hAllHits = std::make_unique<hHits>("AllHits");
}

MuonDTDigis::~MuonDTDigis(){
  if(verbose_)
    cout << "[MuonDTDigis] Destructor called" << endl;
}

void MuonDTDigis::bookHistograms(DQMStore::IBooker & iBooker, edm::Run const & iRun, edm::EventSetup const & /* iSetup */)
{
  meDigiTimeBox_          = 0;
  meDigiTimeBox_wheel2m_  = 0;
  meDigiTimeBox_wheel1m_  = 0;
  meDigiTimeBox_wheel0_   = 0;
  meDigiTimeBox_wheel1p_  = 0;
  meDigiTimeBox_wheel2p_  = 0;
  meDigiEfficiency_       = 0;
  meDigiEfficiencyMu_     = 0;
  meDoubleDigi_           = 0;
  meSimvsDigi_            = 0;
  meWire_DoubleDigi_      = 0;

  meMB1_sim_occup_        = 0;
  meMB1_digi_occup_       = 0;
  meMB2_sim_occup_        = 0;
  meMB2_digi_occup_       = 0;
  meMB3_sim_occup_        = 0;
  meMB3_digi_occup_       = 0;
  meMB4_sim_occup_        = 0;
  meMB4_digi_occup_       = 0;

  meDigiHisto_            = 0;

  // ----------------------
  // We go
  // ----------------------

  Char_t histo_n[100];
  Char_t histo_t[100];

  iBooker.setCurrentFolder("MuonDTDigisV/DTDigiValidationTask");

  sprintf (histo_n, "DigiTimeBox" );
  sprintf (histo_t, "Digi Time Box" );
  meDigiTimeBox_ = iBooker.book1D(histo_n, histo_t, 1536,0,1200);

  sprintf (histo_n, "DigiTimeBox_wheel2m" );
  sprintf (histo_t, "Digi Time Box wheel -2" );
  meDigiTimeBox_wheel2m_ = iBooker.book1D(histo_n, histo_t, 384,0,1200);

  sprintf (histo_n, "DigiTimeBox_wheel1m" );
  sprintf (histo_t, "Digi Time Box wheel -1" );
  meDigiTimeBox_wheel1m_ = iBooker.book1D(histo_n, histo_t, 384,0,1200);

  sprintf (histo_n, "DigiTimeBox_wheel0" );
  sprintf (histo_t, "Digi Time Box wheel 0" );
  meDigiTimeBox_wheel0_ = iBooker.book1D(histo_n, histo_t, 384,0,1200);

  sprintf (histo_n, "DigiTimeBox_wheel1p" );
  sprintf (histo_t, "Digi Time Box wheel 1" );
  meDigiTimeBox_wheel1p_ = iBooker.book1D(histo_n, histo_t, 384,0,1200);

  sprintf (histo_n, "DigiTimeBox_wheel2p" );
  sprintf (histo_t, "Digi Time Box wheel 2" );
  meDigiTimeBox_wheel2p_ = iBooker.book1D(histo_n, histo_t, 384,0,1200);

  sprintf (histo_n, "DigiEfficiencyMu" );
  sprintf (histo_t, "Ratio (#Digis Mu)/(#SimHits Mu)" );
  meDigiEfficiencyMu_ = iBooker.book1D(histo_n, histo_t, 100, 0., 5.);

  sprintf (histo_n, "DigiEfficiency" );
  sprintf (histo_t, "Ratio (#Digis)/(#SimHits)" );
  meDigiEfficiency_ = iBooker.book1D(histo_n, histo_t, 100, 0., 5.);

  sprintf (histo_n, "Number_Digi_per_layer" );
  sprintf (histo_t, "Number_Digi_per_layer" );
  meDoubleDigi_ = iBooker.book1D(histo_n, histo_t, 10,0.,10.);

  sprintf (histo_n, "Number_simhit_vs_digi" );
  sprintf (histo_t, "Number_simhit_vs_digi" );
  meSimvsDigi_ = iBooker.book2D(histo_n, histo_t, 100, 0., 140., 100, 0., 140.);

  sprintf (histo_n, "Wire_Number_with_double_Digi" );
  sprintf (histo_t, "Wire_Number_with_double_Digi" );
  meWire_DoubleDigi_ = iBooker.book1D(histo_n, histo_t, 100,0.,100.);

  sprintf (histo_n, "Simhit_occupancy_MB1" );
  sprintf (histo_t, "Simhit_occupancy_MB1" );
  meMB1_sim_occup_ = iBooker.book1D(histo_n, histo_t, 55, 0., 55. );

  sprintf (histo_n, "Digi_occupancy_MB1" );
  sprintf (histo_t, "Digi_occupancy_MB1" );
  meMB1_digi_occup_ = iBooker.book1D(histo_n, histo_t, 55, 0., 55. );

  sprintf (histo_n, "Simhit_occupancy_MB2" );
  sprintf (histo_t, "Simhit_occupancy_MB2" );
  meMB2_sim_occup_ = iBooker.book1D(histo_n, histo_t, 63, 0., 63. );

  sprintf (histo_n, "Digi_occupancy_MB2" );
  sprintf (histo_t, "Digi_occupancy_MB2" );
  meMB2_digi_occup_ = iBooker.book1D(histo_n, histo_t, 63, 0., 63. );

  sprintf (histo_n, "Simhit_occupancy_MB3" );
  sprintf (histo_t, "Simhit_occupancy_MB3" );
  meMB3_sim_occup_ = iBooker.book1D(histo_n, histo_t, 75, 0., 75. );

  sprintf (histo_n, "Digi_occupancy_MB3" );
  sprintf (histo_t, "Digi_occupancy_MB3" );
  meMB3_digi_occup_ = iBooker.book1D(histo_n, histo_t, 75, 0., 75. );

  sprintf (histo_n, "Simhit_occupancy_MB4" );
  sprintf (histo_t, "Simhit_occupancy_MB4" );
  meMB4_sim_occup_ = iBooker.book1D(histo_n, histo_t, 99, 0., 99. );

  sprintf (histo_n, "Digi_occupancy_MB4" );
  sprintf (histo_t, "Digi_occupancy_MB4" );
  meMB4_digi_occup_ = iBooker.book1D(histo_n, histo_t, 99, 0., 99. );

  // Begona's option
  char stringcham[40];
  for ( int slnum = 1; slnum < 62; ++slnum ) {
    sprintf(stringcham, "DigiTimeBox_slid_%d", slnum) ;
    meDigiHisto_ =  iBooker.book1D(stringcham, stringcham, 100,0,1200);
    meDigiTimeBox_SL_.push_back(meDigiHisto_);
  }  
}

void  MuonDTDigis::analyze(const Event & event, const EventSetup& eventSetup){

  if(verbose_)
    cout << "--- [MuonDTDigis] Analysing Event: #Run: " << event.id().run()
         << " #Event: " << event.id().event() << endl;

  // Get the DT Geometry
  ESHandle<DTGeometry> muonGeom;
  eventSetup.get<MuonGeometryRecord>().get(muonGeom);

  // Get the Digi collection from the event
  Handle<DTDigiCollection> dtDigis;
  event.getByToken(DigiToken_, dtDigis);

  // Get simhits
  Handle<PSimHitContainer> simHits;
  event.getByToken(SimHitToken_, simHits);


  int num_mudigis;
  int num_musimhits;
  int num_digis;
  int num_digis_layer;
  int cham_num ;
  int wire_touched;
  num_digis = 0;
  num_mudigis = 0;
  num_musimhits = 0;
  DTWireIdMap wireMap;

  for(vector<PSimHit>::const_iterator hit = simHits->begin();
      hit != simHits->end(); hit++){
    // Create the id of the wire, the simHits in the DT known also the wireId
    DTWireId wireId(hit->detUnitId());
    //   cout << " PSimHits wire id " << wireId << " part type " << hit->particleType() << endl;

    // Fill the map
    wireMap[wireId].push_back(&(*hit));

    LocalPoint entryP = hit->entryPoint();
    LocalPoint exitP = hit->exitPoint();
    int partType = hit->particleType();
    if ( abs(partType) == 13 ) num_musimhits++;

    if ( wireId.station() == 1 && abs(partType) == 13 ) meMB1_sim_occup_->Fill(wireId.wire());
    if ( wireId.station() == 2 && abs(partType) == 13 ) meMB2_sim_occup_->Fill(wireId.wire());
    if ( wireId.station() == 3 && abs(partType) == 13 ) meMB3_sim_occup_->Fill(wireId.wire());
    if ( wireId.station() == 4 && abs(partType) == 13 ) meMB4_sim_occup_->Fill(wireId.wire());


    float path = (exitP-entryP).mag();
    float path_x = fabs((exitP-entryP).x());

    hAllHits->Fill(entryP.x(),exitP.x(),
           entryP.y(),exitP.y(),
           entryP.z(),exitP.z(),
           path , path_x,
           partType, hit->processType(),
                   hit->pabs());

  }

  //  cout << "num muon simhits " << num_musimhits << endl;

  DTDigiCollection::DigiRangeIterator detUnitIt;
  for (detUnitIt=dtDigis->begin();
       detUnitIt!=dtDigis->end();
       ++detUnitIt){

    const DTLayerId& id = (*detUnitIt).first;
    const DTDigiCollection::Range& range = (*detUnitIt).second;

    num_digis_layer = 0 ;
    cham_num = 0 ;
    wire_touched = 0;

    // Loop over the digis of this DetUnit
    for (DTDigiCollection::const_iterator digiIt = range.first;
     digiIt!=range.second;
     ++digiIt){
      //   cout<<" Wire: "<<(*digiIt).wire()<<endl
      //  <<" digi time (ns): "<<(*digiIt).time()<<endl;

      num_digis++;
      num_digis_layer++;
      if (num_digis_layer > 1 )
      {
        if ( (*digiIt).wire() == wire_touched )
        {
          meWire_DoubleDigi_->Fill((*digiIt).wire());
          //      cout << "old & new wire " << wire_touched << " " << (*digiIt).wire() << endl;
        }
      }
      wire_touched = (*digiIt).wire();

      meDigiTimeBox_->Fill((*digiIt).time());
      if (id.wheel() == -2 ) meDigiTimeBox_wheel2m_->Fill((*digiIt).time());
      if (id.wheel() == -1 ) meDigiTimeBox_wheel1m_->Fill((*digiIt).time());
      if (id.wheel() == 0 )  meDigiTimeBox_wheel0_ ->Fill((*digiIt).time());
      if (id.wheel() == 1 )  meDigiTimeBox_wheel1p_->Fill((*digiIt).time());
      if (id.wheel() == 2 )  meDigiTimeBox_wheel2p_->Fill((*digiIt).time());

      //   Superlayer number and fill histo with digi timebox

      cham_num = (id.wheel() +2)*12 + (id.station() -1)*3 + id.superlayer();
      //   cout << " Histo number " << cham_num << endl;

      meDigiTimeBox_SL_[cham_num]->Fill((*digiIt).time());

      //    cout << " size de digis " << (*digiIt).size() << endl;

      DTWireId wireId(id,(*digiIt).wire());
      if (wireId.station() == 1 ) meMB1_digi_occup_->Fill((*digiIt).wire());
      if (wireId.station() == 2 ) meMB2_digi_occup_->Fill((*digiIt).wire());
      if (wireId.station() == 3 ) meMB3_digi_occup_->Fill((*digiIt).wire());
      if (wireId.station() == 4 ) meMB4_digi_occup_->Fill((*digiIt).wire());

      int mu=0;
      float theta = 0;

      for(vector<const PSimHit*>::iterator hit = wireMap[wireId].begin();
      hit != wireMap[wireId].end(); hit++)
    if( abs((*hit)->particleType()) == 13){
      theta = atan( (*hit)->momentumAtEntry().x()/ (-(*hit)->momentumAtEntry().z()) )*180/M_PI;
          //      cout<<"momentum x: "<<(*hit)->momentumAtEntry().x()<<endl
          //          <<"momentum z: "<<(*hit)->momentumAtEntry().z()<<endl
          //          <<"atan: "<<theta<<endl;
      mu++;
    }

      if( mu ) num_mudigis++;

      if(mu && theta){
    hDigis_global->Fill((*digiIt).time(),theta,id.superlayer());
    //filling digi histos for wheel and for RZ and RPhi
    WheelHistos(id.wheel())->Fill((*digiIt).time(),theta,id.superlayer());
      }

    }// for digis in layer

    meDoubleDigi_->Fill( (float)num_digis_layer );

  }// for layers

  //cout << "num_digis " << num_digis << "mu digis " << num_mudigis << endl;

  if (num_musimhits != 0) {
    meDigiEfficiencyMu_->Fill( (float)num_mudigis/(float)num_musimhits );
    meDigiEfficiency_->Fill( (float)num_digis/(float)num_musimhits );
  }

  meSimvsDigi_->Fill( (float)num_musimhits, (float)num_digis ) ;
  //  cout<<"--------------"<<endl;

}

hDigis* MuonDTDigis::WheelHistos(int wheel){
  switch(abs(wheel)){

    case 0: return  hDigis_W0.get();

    case 1: return  hDigis_W1.get();

    case 2: return  hDigis_W2.get();

    default: return NULL;
  }
}
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "DQMServices/Core/interface/DQMStore.h"

DEFINE_FWK_MODULE(MuonDTDigis);