/******* \class DTAnalyzer *******
 *
 * Description:
 *  
 *  detailed description
 *
 * \author : Stefano Lacaprara - INFN LNL <stefano.lacaprara@pd.infn.it>
 * $date   : 20/11/2006 16:50:57 CET $
 *
 * Modification:
 *
 *********************************/

/* This Class Header */
#include "RecoLocalMuon/DTSegment/test/DTAnalyzer.h"

/* Collaborating Class Header */
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Utilities/interface/Exception.h"

using namespace edm;

#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"

#include "DataFormats/LTCDigi/interface/LTCDigi.h"

#include "CalibMuon/DTDigiSync/interface/DTTTrigBaseSync.h"
#include "CalibMuon/DTDigiSync/interface/DTTTrigSyncFactory.h"
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
#include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
#include "RecoLocalMuon/DTSegment/test/DTMeanTimer.h"
#include "RecoLocalMuon/DTSegment/test/DTSegmentResidual.h"


#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"

// this is to handle extrapolations
#include "DataFormats/GeometrySurface/interface/Cylinder.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"


// this is to retrieve digi's
#include "DataFormats/DTDigi/interface/DTDigiCollection.h"
#include "DataFormats/MuonDetId/interface/DTWireId.h"

#include "DataFormats/DTDigi/interface/DTLocalTriggerCollection.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"

#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackReco/interface/Track.h"

/* C++ Headers */
#include <iostream>
#include <cmath>
using namespace std;

/* ====================================================================== */

/* Constructor */ 
DTAnalyzer::DTAnalyzer(const ParameterSet& pset) : _ev(0){

  theSync =
    DTTTrigSyncFactory::get()->create(pset.getUntrackedParameter<string>("tTrigMode"),
                                      pset.getUntrackedParameter<ParameterSet>("tTrigModeConfig"));
  // Get the debug parameter for verbose output
  debug = pset.getUntrackedParameter<bool>("debug");
  theRootFileName = pset.getUntrackedParameter<string>("rootFileName");

  // the name of the digis collection
  theDTLocalTriggerLabel = pset.getParameter<string>("DTLocalTriggerLabel");

  // the name of the 1D rec hits collection
  theRecHits1DLabel = pset.getParameter<string>("recHits1DLabel");

  // the name of the 2D rec hits collection
  theRecHits2DLabel = pset.getParameter<string>("recHits2DLabel");

  // the name of the 4D rec hits collection
  theRecHits4DLabel = pset.getParameter<string>("recHits4DLabel");

  // the name of the 4D rec hits collection
  theSTAMuonLabel = pset.getParameter<string>("SALabel");

  // the name of the propagator
  thePropagatorName = pset.getParameter<std::string>("PropagatorName");
  thePropagator = 0;

  // if MC
  mc = pset.getParameter<bool>("isMC");

  LCT_RPC = pset.getParameter<bool>("LCT_RPC");
  LCT_DT = pset.getParameter<bool>("LCT_DT");
  LCT_CSC = pset.getParameter<bool>("LCT_CSC");


  doTrig = pset.getParameter<bool>("doTrig");
  doHits = pset.getParameter<bool>("doHits");
  doSegs = pset.getParameter<bool>("doSegs");
  doSA = pset.getParameter<bool>("doSA");


  // Create the root file
  theFile = new TFile(theRootFileName.c_str(), "RECREATE");

// ==================
// histogram booking
// ==================

// global histos...------------
// 1d recHits (in whole detector)
  new TH1F("hnHitDT","Num 1d hits DT",200,0.,200.);
// segments (in whole detector)
  new TH1F("hnSegDT","Num seg DT",50,0.,50.);  
// Trigger:
  new TH1F ("TriggerInclusive","Trigger Inclusive",60,0.5,60.5);
  for (int ise=0; ise<12; ise++) sects[ise]=0;
  new TH1F ("SectorTriggerMatrix","SectorTrigger Matrix",21,-0.5,20.5);
  
// plots per Sector/CHambers / SL ------------------

   for ( int sect = 1; sect < 15; sect++) {   // Sector loop

//   occupancy digi scatter plots
     stringstream hNameOcc; 
     stringstream hTitleOcc; 
     hNameOcc  << "hDigiXY_S" << sect ;
     hTitleOcc << "Digis in Sect" << sect ;
     createTH2F(hNameOcc.str(),hTitleOcc.str(),"",100, 0.,100.,80,0.,80.); // ** histos hDigiXY_S1, ecc... ********

//   sectors trigger matrix
     stringstream hTitleTrigMat;
     stringstream hNameTrigMat;
     hNameTrigMat << "TriggerMatrix" << sect ;
     hTitleTrigMat << "Trigger Matrix, Sector " << sect ;
     createTH1F (hNameTrigMat.str(),hTitleTrigMat.str(),"",20,-0.5,19.5);

     for  ( int jst = 1; jst < 5; jst++ ) {   // station loop

//   trigger Histos: quality & bxID ----------------------
     stringstream hNameTrigg; 
     stringstream hTitleTrigg; 
     hNameTrigg  << "hTrigBits_S" << sect << "_MB" << jst ;
     hTitleTrigg << "Highest trigg qual,S" << sect << "/MB" << jst ;
     createTH1F(hNameTrigg.str(),hTitleTrigg.str(),"",11,-1.,10.); // *** histos: hTrigBits_S1_MB1, ecc... **********
     for  ( int iqual = 1; iqual < 7; iqual++ ) { 
       stringstream hNameTriggBX; 
       stringstream hTitleTriggBX; 
       if (iqual == 1 || iqual > 3) {  // book histos for qual=1 (== all uncorrelated trigger), 4=LL, 5=HL, 6=HH
        hNameTriggBX  << "hTrigBX_S" << sect << "_MB" << jst << "_qual" << iqual ;
        hTitleTriggBX << "Highest qual BXid,S" << sect << "/MB" << jst  << " qual=" << iqual;
        createTH1F(hNameTriggBX.str(),hTitleTriggBX.str(),"",30,0.,30.); // *** histos: hTrigBX_S1_MB1_qual1, ecc... **********
       }
     } // next quality flag 

        
//   Digi histos: ------------------------
//   time boxes per chamber....
//     float tmin = 3500., tmax = 4500.;  // full YB0, local Runs 25946,25602,...
//     float tmin = 2300., tmax = 3300.;  // Sept Global Runs
//    float tmin = 2800., tmax = 3800.;  // GREN....
     float tmin = -200., tmax = 800.;  // for tbox ttrigger-subtracted
     stringstream hNameDigi; 
     stringstream hTitleDigi; 
     hNameDigi << "htimeYB0_S" << sect << "_MB" << jst  ;
     hTitleDigi << "Tbox S" << sect << "/MB" << jst  ;
     createTH1F(hNameDigi.str(),hTitleDigi.str(),"",100,tmin, tmax); // *** histos: htimeYB0_S1_MB1, ecc... **********
     for (int sl = 1; sl < 4; sl++ ) {  // histograms per SL
//     time boxes per SL....
       stringstream hTitleDigiSL; 
       hTitleDigiSL << "htimeYB0_S" << sect << "_MB" << jst << "_SL" << sl  ;
       createTH1F(hTitleDigiSL.str(),"Digi time per SL","",100,tmin, tmax); // *** histos: htimeYB0_S1_MB1_SL1, ecc... **********
//     tmax plots per SL...
       float t1 = 300., t2 = 450. ;
       stringstream hNameTmax; 
       stringstream hTitleTmax; 
       hNameTmax  << "htmaxYB0_S" << sect << "_MB" << jst << "_SL" << sl  ;
       hTitleTmax << "Tmax S" << sect << "/MB" << jst << "/SL" << sl  ;
       createTH1F(hNameTmax.str(),hTitleTmax.str(),"",100,t1, t2); // *** histos: htmaxYB0_S1_MB1_SL1, ecc... **********
     } // next SL

//   hit residuals in segments -------------
     float xmin = -0.25, xmax = 0.25 ; // cm
     stringstream hNameRes; 
     stringstream hTitleRes; 
     hNameRes << "hResX_S" << sect << "_MB" << jst ;
     hTitleRes << "Hit residuals S" << sect << "/MB" << jst ;
     createTH1F(hNameRes.str(),hTitleRes.str(),"",100,xmin, xmax); // *** histos: hResX_S1_MB1, ecc... **********

//   philocal of rec segments
     float phimin = -50., phimax = 50.;
     stringstream hNamePhi; 
     stringstream hTitlePhi; 
     hNamePhi  << "hPhi_S" << sect << "_MB" << jst ;
     hTitlePhi << "Phi local S" << sect << "/MB" << jst ;
     createTH1F(hNamePhi.str(),hTitlePhi.str(),"",100,phimin, phimax); // *** histos: hPhi_S1_MB1, ecc... **********

//    hits in track segments
     stringstream hNameHits; 
     stringstream hTitleHits; 
     hNameHits  << "hNhits_S" << sect << "_MB" << jst ;
     hTitleHits << "Nr.hits in segment, S" << sect << "/MB" << jst ;
     createTH1F(hNameHits.str(),hTitleHits.str(),"",15,0, 15); // *** histos: hNhits_S1_MB1, ecc... **********

//   mis-alignment histos
     stringstream hNameAliX; 
     stringstream hTitleAliX; 
     hNameAliX  << "halignX_S" << sect << "_MB" << jst ;
     hTitleAliX << "X-Misalignment, S" << sect << "/MB" << jst ;
     createTH1F(hNameAliX.str(),hTitleAliX.str(),"",100,-5.0, 5.0); // *** histos: halignX_S1_MB1, ecc... **********
     stringstream hNameAliZ; 
     stringstream hTitleAliZ; 
     hNameAliZ  << "halignZ_S" << sect << "_MB" << jst ;
     hTitleAliZ << "X-Misalignment, S" << sect << "/MB" << jst ;
     createTH1F(hNameAliZ.str(),hTitleAliZ.str(),"",100,-5.0, 5.0); // *** histos: halignZ_S1_MB1, ecc... **********


   }   // next station ===========================

// segments in Sector
   stringstream hNameNsegm; 
   stringstream hTitleNsegm; 
   hNameNsegm << "hNsegs_S" << sect  ;
   hTitleNsegm << "Segments in sect " << sect  ;
   createTH1F(hNameNsegm.str(),hTitleNsegm.str(),"",10, 0, 10); // *** histos: hNsegs_S1 ecc... **********

// associated hits to global track per sector
   stringstream hNameNassH; 
   stringstream hTitleNassH; 
   hNameNassH<< "hNhass_S" << sect  ;
   hTitleNassH<< "Associated hits in sect." << sect  ;
   createTH1F(hNameNassH.str(),hTitleNassH.str(),"",50, 0, 50); // *** histos: hNhass_S1 ecc... **********
 
  } // next sector ============================
          

	  
  // STA Muon plots
  new TH1F("hNSA","Num SA tracks in events", 6, -0.5, 5.5);
  new TH1F("hNhitsSA","Num hits in SA tracks", 50, .5, 50.5);
  new TH1F("hChi2SA","#chi^{2}/NDoF for SA tracks", 100, 0, 100.);
  new TH1F("hSAValidHits","#Valid Hits for SA tracks", 100, 0, 100.);
  new TH1F("hSALostHits","#Lost Hits for SA tracks", 100, 0, 100.);
  new TH1F("hPIPSA","p for SA tracks @ IP", 100, 0., 100);
  new TH1F("hPtIPSA","pt for SA tracks @ IP", 100, 0., 100);
//  new TH1F("hPhiIPSA","#phi for SA tracks @ IP", 100, -M_PI_2, M_PI_2);
  new TH1F("hPhiIPSA","#phi for SA tracks @ IP", 100, -3.14, 3.14);
  new TH1F("hEtaIPSA","#eta for SA tracks @ IP", 100, -2.5, 2.5);
  new TH1F("hrIPSA"," r at IP",100,0.,500.);
  new TH1F("hzIPSA"," z at IP",100,-750.,750.);
  new TH2F("hrVszIPSA"," r vs z at IP",100,-750.,750.,100,0.,500.);

  new TH1F("hPInnerTSOSSA","p for SA tracks @ innermost TSOS", 100, 0., 100);
  new TH1F("hPtInnerTSOSSA","pt for SA tracks @ innermost TSOS", 100, 0., 100);
//  new TH1F("hPhiInnerTSOSSA","#phi for SA tracks @ innermost TSOS", 100, -M_PI_2, M_PI_2);
  new TH1F("hPhiInnerTSOSSA","#phi for SA tracks @ innermost TSOS", 100, -3.14, 3.14);
  new TH1F("hEtaInnerTSOSSA","#eta for SA tracks @ innermost TSOS", 100, -2.5, 2.5);
  new TH1F("hInnerRSA","Radius of innermost TSOS for SA tracks", 100, 400, 1000.);
  new TH1F("hOuterRSA","Radius of outermost TSOS for SA tracks", 100, 400, 1000.);
  new TH2F("hInnerOuterRSA","Radius of outermost TSOS vs innermost one for SA tracks",
           40, 400, 1000.,40, 400, 1000.);

                  
  new TH2F("hNSAVsNHits",  "Nr.SA tracks vs Nr.1d Hits", 100,-0.5,39.5, 5,-0.5, 4.5);
  new TH2F("hNSAVsNSegs2D","Nr.SA tracks vs Nr.2d Segs", 20,-0.5,19.5, 5,-0.5, 4.5);
  new TH2F("hNSAVsNSegs4D","Nr.SA tracks vs Nr.4d Segs", 10,-0.5,9.5, 5,-0.5, 4.5);

  new TH2F("hNHitsSAVsNHits",  "N Hits SA vs Nr. 1d Hits", 40,-0.5,39.5,100,0,100);
  new TH2F("hNHitsSAVsNSegs2D","N Hits SA vs Nr. 2d Segs", 20,-0.5,19.5,100,0,100);
  new TH2F("hNHitsSAVsNSegs4D","N Hits SA vs Nr. 4d Segs", 10,-0.5,9.5, 100,0,100);

  new TH2F("hHitsPosXYSA","Hits position (x,y) SA",150,-750,750,150,-750,750);
  new TH2F("hHitsPosXYSA_1","Hits position (x,y) SA",150,0,750,150,0,750);
  new TH2F("hHitsPosXYSA_2","Hits position (x,y) SA",150,-750,0,150,0,750);
  new TH2F("hHitsPosXYSA_3","Hits position (x,y) SA",150,-750,0,150,-750,0);
  new TH2F("hHitsPosXYSA_4","Hits position (x,y) SA",150,0,750,150,-750,0);
  new TH2F("hHitsPosXZSA","Hits position (x,z) SA",100,-700,700,100,-800,800);
  new TH2F("hHitsPosYZSA","Hits position (y,z) SA",100,-700,700,100,-800,800);

// Global Phi of associated rechits  
  for  ( int jst = 1; jst < 5; jst++ ) {   // station loop
     stringstream hTitlePhihit;
     stringstream hNamePhihit;
     hNamePhihit << "hPhiHit_MB" << jst ;
     hTitlePhihit << "Phi of ass.hit, MB" << jst ;
     createTH1F(hNamePhihit.str(),hTitlePhihit.str(),"",180,-180., 180.); // *** histos: hPhiHit_MB1 ecc... **********
     stringstream hNamePhiGlob;
     stringstream hTitlePhiGlob;
     hNamePhiGlob << "hPhiGlob_MB" << jst;
     hTitlePhiGlob << "Phi-Trigger eff, MB" << jst;
     createTH1F(hNamePhiGlob.str(),hTitlePhiGlob.str(),"",180,-180., 180.); // *** histos: hPhiGlob_MB1 ecc... **********
     stringstream hNamePhiTrigg;
     stringstream hTitlePhiTrigg;
     hNamePhiTrigg << "hPhiTrigg_MB" << jst;
     hTitlePhiTrigg << "Phi-Trigger eff, MB"  << jst;
     createTH1F(hNamePhiTrigg.str(),hTitlePhiTrigg.str(),"",180,-180., 180.); // *** histos: hPhiTrigg_MB1 ecc... **********
  } // nest station
  
//  extrapolations to ECAL, HCAL  
  float zmin = -400., zmax = 400.;
  new TH2F("hphivszECAL","Extrap.position phi-vs-z to ECAL ",100,zmin,zmax,90,-180,180);
  new TH2F("hphivszHCAL","Extrap.position phi-vs-z to HCAL ",100,zmin,zmax,90,-180,180);

} // end constructor =================================================================

/* Destructor */ 
DTAnalyzer::~DTAnalyzer() {
  theFile->cd();
  if (doTrig) LabelTriggerMatrix();
  theFile->Write();
  theFile->Close();
} // end destructor ================================================================

// tools for trigger analysis.....==============================================
void DTAnalyzer::LabelTriggerMatrix() {

  //Labels for SectorTriggerMatrix:

  char A[3]; char B[3]; char C[3]; char D[3];
  sprintf (A,"%u",sects[0]);  sprintf (B,"%u",sects[1]);  sprintf (C,"%u",sects[2]); sprintf (D,"%u",sects[3]);
  std::string lab1[] = {" "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "," "};
  
  lab1[0]="no trig";
  lab1[1]="other sects";
  lab1[3]= lab1[3] + "Sect" + A;
  if (sects[1]) {
   lab1[4]= lab1[4] + "Sect" + B;
   lab1[8]=lab1[8] + A + " & " + B;
   if (sects[2]) {
    lab1[5]=lab1[5] + "Sect" + C;
    lab1[9]=lab1[9] + A + " & " + C;
    lab1[11]=lab1[11] + B +" & " + C;
    lab1[15]=lab1[15] + A + " & " + B +" & " + C;
    if (sects[3]) {
     lab1[6]=lab1[6] + "Sect" + D;
     lab1[10]=lab1[10] + A + " & " + D;
     lab1[12]=lab1[12] + B + " & " + D;
     lab1[13]=lab1[13] + C + " & " + D;

     lab1[16]=lab1[16] + A + " & " + B + " & " + D;
     lab1[17]=lab1[17] + A + " & " + C + " & " + D;
     lab1[18]=lab1[18] + B + " & " + C + " & " + D;

     lab1[20]=lab1[20] + A + " & " + B + " & " + C + " & " + D;
    }
   }
  }
  for (int ibin=1; ibin<21; ibin++) histo("SectorTriggerMatrix")->GetXaxis()->SetBinLabel(ibin,lab1[ibin-1].c_str()); 

// Labels for TriggerMatrix:

  char* lab[]={"no trig"," ","MB1","MB2","MB3","MB4"," ","1 & 2","1 & 3","1 & 4","2 & 3","2 & 4",
               "3 & 4"," ","1 & 2 & 3","1 & 2 & 4","1 & 3 & 4","2 & 3 & 4"," ","1 & 2 & 3 & 4"};
  for (int ise=1; ise<13; ise++) {
   stringstream hName; 
   hName  << "TriggerMatrix" << ise ;
   for (int ibin=1; ibin<21; ibin++) histo(hName.str())->GetXaxis()->SetBinLabel(ibin,lab[ibin-1]); 
  }

// Labels for TriggerInclusive:
  char label[10];
  for (int ise=1; ise<13; ise++) {
   for (int ibin=1; ibin<5; ibin++) {
    sprintf (label,"S%uMB%u",ise,ibin);
    histo("TriggerInclusive")->GetXaxis()->SetBinLabel((ise-1)*5+ibin,label);
   } 
   histo("TriggerInclusive")->GetXaxis()->SetBinLabel(ise*5," ");
  }

}

// Operations : main module ======================================================================
void DTAnalyzer::analyze(const Event & event,
                         const EventSetup& eventSetup) {
// ===============================================================================================

  theSync->setES(eventSetup);

  int nprt = event.id().event() - 100*(event.id().event()/100);
  if( nprt == 1)
    cout << "Run:Event analyzed " << event.id().run() << ":" << event.id().event() <<
    " Num " << _ev++ << endl;

  if (debug) cout << endl<<"--- [DTAnalyzer] Event analysed #Run: " <<
    event.id().run() << " #Event: " << event.id().event() << " =====================================" <<  endl;
                      
// go on with the rest of analysis: DIGI/RecHits, Segments, STAmuons... ============
  if (selectEvent()) {
    if (doTrig) analyzeTrigger(event, eventSetup);
    if (doHits) analyzeDTHits(event, eventSetup);
    if (doSegs) analyzeDTSegments(event, eventSetup);
    if (doSA) analyzeSATrack(event, eventSetup);
  }

} // end main =========================================================================================== 

bool DTAnalyzer::selectEvent() const {
  bool trigger=true;
//  if( qual[0] > 3 ) trigger = true ; // HH,HL,LL trigger in MB1
  return trigger;
}

// ============================================================================================================
// Hits & RecHits analysis 
void DTAnalyzer::analyzeDTHits(const Event & event,
                               const EventSetup& eventSetup){
// ============================================================================================================

   int Nnoisy = 90;
//  100000*wireId.sector() + 10000*wireId.station()+ 1000*wireId.superlayer()+ 100*wireId.layer()+wireId.wire();
    int noisy_ch[] = {
112254,112255,112256, 112257, 112307, 112353, 112354, 112355, 112356, 112357,
122210,122211,122355, 221106, 221201, 221202, 221203, 221219, 221220, 221303,
222005,222006,222007, 223201, 223203, 231209, 231212, 233324, 233428, 233430,
233431,311348,311349, 321259, 321260, 413349, 422101, 422201, 521303, 611149,
611209,611210,611349, 613201, 613203, 613204, 613301, 613302, 621139, 621140,
621240,621360,622109, 622110, 622112, 622125, 622201, 622209, 622210, 622211,
622212,622226,622302, 631355, 631356, 632056, 643104, 711108, 711208, 721109,
721209,721210,721211, 721360, 812354, 812355, 812356, 911349, 912355, 912357,
921158,921159,943348,1122101,1122357,1132046,1211349,1213248,1213349,1232357 } ;
 

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


   // Get the digis from the event ======================
  Handle<DTDigiCollection> digis; 
//  event.getByLabel(digiLabel, digis);
  event.getByLabel(theDTLocalTriggerLabel, digis);
  
 // Iterate through all digi collections ordered by LayerId  -------------------  

  DTDigiCollection::DigiRangeIterator dtLayerIt;
  for (dtLayerIt = digis->begin();
       dtLayerIt != digis->end();  ++dtLayerIt){   // loop on layerID -------
    // The layerId
    const DTLayerId layerId = (*dtLayerIt).first;
//    const DTSuperLayerId slId = layerId.superlayerId();
    // Get the iterators over the digis associated with this LayerId
    const DTDigiCollection::Range& digiRange = (*dtLayerIt).second;
    // Loop over all digis in the given range
    for (DTDigiCollection::const_iterator digi = digiRange.first;   // loop on digi --------
	 digi != digiRange.second;
	 digi++) {                                        
      const DTWireId wireId(layerId, (*digi).wire());
/*    cout << "Wire: " << wireId << endl;
      cout << "  wheel   " << wireId.wheel() << endl;
      cout << "  sector  " << wireId.sector() << endl;
      cout << "  chamber " << wireId.station() << endl;
      cout << "  SL      " << wireId.superlayer() << endl;
      cout << "  layer   " << wireId.layer() << endl;
      cout << "  wire    " << wireId.wire() << endl;  */
//      int ix = 100*wireId.sector() + wireId.wire();
      int ix = wireId.wire();
      int iy = 20*(wireId.station()-1) + 4*(wireId.superlayer()-1) + wireId.layer();
      int se = wireId.sector();
//    occupancy scatter plot -----------------------------------
      stringstream hTitleDigi; 
      hTitleDigi << "hDigiXY_S" << se ;
      histo2d(hTitleDigi.str())->Fill( ix,iy ); // *** histos: hDigiXY_S1, ecc... **********

// the ttrigg used for this channel
   float ttrig = theSync->offset(wireId);
                                 
//   float TDCtime = (*digi).time()/ 0.78125 ;
//   float TDCtime = (*digi).time() ;
   float TDCtime = (*digi).time() - ttrig ; // plot the TDC times ttrigg-subtracted....


// check noisy channel
   int iskip = 0;   
// flag as defined in privbate macro Noise_ana.C
   int ichannel = 100000*wireId.sector() + 10000*wireId.station()+ 1000*wireId.superlayer()+ 100*wireId.layer()+wireId.wire();
   for (int i = 0; i <= Nnoisy; i++ ) {
     if (ichannel == noisy_ch[i] ) iskip = 1;
   }
   if (iskip == 0)  {   // apply noise filter
//   if ( wireId.wheel() == 0 ) {  // YB0 =========================
        stringstream hTitle, hTitleSL; 
        hTitle << "htimeYB0_S" << wireId.sector() << "_MB" << wireId.station() ;
        histo(hTitle.str())->Fill( TDCtime ); // *** histos: htimeYB0_1_MB1, ecc... **********
        hTitleSL << "htimeYB0_S" << wireId.sector() << "_MB" << wireId.station() << "_SL"<< wireId.superlayer() ;
        histo(hTitleSL.str())->Fill( TDCtime ); // *** histos: htimeYB0_S1_MB1_SL1, ecc... **********		
//    } // endif YB0 ==============================================
   } // endif noise filter

    
   } // next digi --------------------
  } // next LayerId -----------------


//  get the tracking geometry
  ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
  eventSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);


  // Get the 1D rechits from the event --------------
  Handle<DTRecHitCollection> dtRecHits; 
  event.getByLabel(theRecHits1DLabel, dtRecHits);

  int nHitDT = dtRecHits->size();
  histo("hnHitDT")->Fill(nHitDT);

// ********* recHits loop ***********************
/*
  for (DTRecHitCollection::const_iterator hit = dtRecHits->begin();
       hit!=dtRecHits->end();  ++hit) {
    // Get the wireId of the rechit
    DTWireId wireId = (*hit).wireId();
    cout << "Wire: " << wireId << endl;
      cout << "  wheel   " << wireId.wheel() << endl;
      cout << "  sector  " << wireId.sector() << endl;
      cout << "  chamber " << wireId.station() << endl;
      cout << "  SL      " << wireId.superlayer() << endl;
      cout << "  layer   " << wireId.layer() << endl;
      cout << "  wire    " << wireId.wire() << endl;  

    float ttrig = theSync->offset(wireId);
    //cout << "TTrig " << ttrig << endl;
    float time = (*hit).digiTime()  - ttrig ;
    double xLeft = (*hit).localPosition(DTEnums::Left).x();
    double xRight = (*hit).localPosition(DTEnums::Right).x();
  } // ********* end recHits loop  *******************
*/

} // end DTHits analyzer ===============================================



// ==========================================================================================================
// Segment analysis:
void DTAnalyzer::analyzeDTSegments(const Event & event,
                                   const EventSetup& eventSetup){
// ==========================================================================================================

  if (debug) cout << " analyzeDTSegments: " << " Run/Event " << event.id().run() << ":" << event.id().event() << endl;


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

  // Get the digis from the event
  Handle<DTLocalTriggerCollection> digis; 
  if (!mc) event.getByLabel(theDTLocalTriggerLabel, digis);

  // Get the 4D rechit collection from the event -------------------
  edm::Handle<DTRecSegment4DCollection> segs;
  event.getByLabel(theRecHits4DLabel, segs);
  if (debug) cout << "4d segments: " << segs->size() << endl;

  // Get the 2D rechit collection from the event -------------------
  edm::Handle<DTRecSegment2DCollection> segs2d;
  event.getByLabel(theRecHits2DLabel, segs2d);
  if (debug) cout << "2d segments:" << segs2d->size() << endl;

  // Get the 1D rechits from the event --------------
  Handle<DTRecHitCollection> dtRecHits; 
  event.getByLabel(theRecHits1DLabel, dtRecHits);
  if (debug) cout << "1d rechits: " << dtRecHits->size() << endl;

  int nsegs = segs->size();
  histo("hnSegDT")->Fill(nsegs);
  const std::vector<DTChamber*> & chs = dtGeom->chambers();


// =========================================================
//  local trigger analysis for trigg. efficiency computation
//==========================================================
  // get the DT local trigger collection =======================
  edm::Handle<DTLocalTriggerCollection> allLocalTriggers;
  event.getByLabel("dtunpacker", allLocalTriggers); 
 DTLocalTriggerCollection::DigiRangeIterator chambIt;
// Loop over chambers present in DTLocalTriggerCollection
  bool hasTr[]={false,false,false,false,false};
  bool hasTrOut[]={false,false,false,false,false};
  int bx[4][5]; int qual[4][14];
  for (int sec=1; sec<15; ++sec) { // section 1 to 14
    for (int ist=0; ist<5; ++ist) { 
      qual[ist-1][sec-1] = -1;
    }    
  } 
  int SCsect=0; int SCst=0;
  for (chambIt=allLocalTriggers->begin();chambIt!=allLocalTriggers->end();++chambIt){ // loop on chambers ------
    const DTChamberId& id = (*chambIt).first;
    SCst=id.station();
    SCsect=id.sector();
    bx[SCst-1][0]=0;     bx[SCst-1][1]=0;
    const DTLocalTriggerCollection::Range& range = (*chambIt).second;    
    int ntrCh=0;
//  loop over triggers of this chamber
    for (DTLocalTriggerCollection::const_iterator trigtrack = range.first;trigtrack!=range.second;++trigtrack){
      if (trigtrack->trOut()) hasTrOut[SCst-1]=true;
      if (trigtrack->quality()<7) {
       hasTr[SCst-1]=true; // SCst = numero stazione : la stazione ha un trigger
       bx[SCst-1][ntrCh]=trigtrack->bx();  // dice a che bx c'e' stato trigger ..  
 //      float ibx=trigtrack->bx();
       ntrCh++; 
//       cout << " SCsect: " << SCsect << " SCst: " << SCst << " quality " << trigtrack->quality() << endl;
       if (ntrCh>4) break;
       if (qual[SCst-1][SCsect-1]!=-2) { 
         int iQual=trigtrack->quality();
         if( iQual > qual[SCst-1][SCsect-1] ) qual[SCst-1][SCsect-1]= iQual;  // store the best trigger quality for this station
       } //end if qual != -2
      } //  quality <7 
    } // end loop on triggers in this chambers -------------    
  } //end loop on chambers -------  
// ===== end local trigger analysis =========================================   

// counters for segments in sectors
   int nSegSect[15];
    for (int sec=1; sec<=15; ++sec) { // section 1 to 14
      nSegSect[sec] = 0;
    }


// loop on chambers ==============================================
  for (std::vector<DTChamber*>::const_iterator ch = chs.begin();
       ch!=chs.end() ; ++ch) {
    DTChamberId chid((*ch)->id());
    //cout << "chid " << chid << endl;
    //int w= chid.wheel();
   // int se= chid.sector();
    //int st= chid.station();
    DTRecSegment4DCollection::range segsch= segs->get(chid);


// ********** loop on 4D segments *********************************

    for (DTRecSegment4DCollection::const_iterator seg=segsch.first ;
         seg!=segsch.second ; ++seg ) {

      int ist    = seg->chamberId().station();
 //     int iwheel = seg->chamberId().wheel() + 3;  // Wheel -2,.+2  => 1,...5
      int isect  = seg->chamberId().sector();
      int iisect = isect;
      if (isect == 13 ) isect = 4;
      if (isect == 14 ) isect = 10;
      
            const DTChamber* ch = dtGeom->chamber(seg->chamberId());
            GlobalPoint   glbPoint = ch->toGlobal((*seg).localPosition());
            GlobalVector  glbDir   = ch->toGlobal((*seg).localDirection());

      // first the the two projection separately -----------------
      // phi segment
      if(debug) cout  << seg->chamberId() << endl;
//      if(debug) cout  << " 4D  segm: " << *seg << endl;
      if(debug) cout  << " 4D  segm: x " << glbPoint.x() << " y " << glbPoint.y() << " z " << glbPoint.z() << endl;
      const DTChamberRecSegment2D* phiSeg= (*seg).phiSegment();
//      cout << (*seg).localDirection() << endl;
//      cout <<  -atan(glbDir.x()/ glbDir.y()) << endl;
      float localPhi = (*seg).localDirection().x();
      float radtodeg = 57.296 ;
      float localPhideg = radtodeg * atan(localPhi);
//      vector<DTRecHit1D> phiHits;   
      std::vector<DTRecHit1D> phiHits; 
      int NtkHit = 0;    
      if (phiSeg) {
//        if(debug) cout << "Phi segm: " << *phiSeg << endl;
        phiHits = phiSeg->specificRecHits();
        if(debug) cout << " Nhits in Phi-segment " << phiHits.size() << endl;
        NtkHit = phiHits.size();
        if (phiHits.size()>=6)      { // Nhits > 5
//	 cout << " Event " << event.id().event() << " Sect " << isect << " MB" << ist << " Philocal " << localPhideg << endl;
        stringstream hTitlePhi;
        hTitlePhi << "hPhi_S" << isect << "_MB" << ist;
        histo(hTitlePhi.str())->Fill( localPhideg ); // *** histos: hPhiS1_MB1, ecc... **********

	
//        ---- trigger efficiency analysis -----------
        if (abs(localPhideg)< 30. )   { //  phi < 30 deg
//        set flag for trigger efficiency study: 1= use it (event triggered by other stations), 0= don't use it
          int MBeff = 0;  // don't use this station [ist,isec ]
          for (int sec1=1; sec1<15; ++sec1)       {
           for ( int ist1=1; ist1 < 5; ++ist1) { 
              if( ist1 != ist  || sec1 != iisect ) {
                if ( qual[ist1-1][sec1-1] > 4 ) MBeff = 1; // event triggered also by other stations; use this station
 	      }  
            } // next station
          } // next sector
          if( MBeff == 1) { // event triggered by other stations
            float phiglobal =  glbPoint.phi(); // phi position of reco segments
	    if (phiHits.size()>=7) {  // HH & HL potential candidates
              stringstream hTitlePhiGlob;
              hTitlePhiGlob << "hPhiGlob_MB" << ist;
              stringstream hTitlePhiTrigg;
              hTitlePhiTrigg << "hPhiTrigg_MB" << ist;
              histo(hTitlePhiGlob.str())->Fill( phiglobal*radtodeg  ); // *** histos: hPhiGlob_MB1, ecc... **********
              if( qual[ist-1][iisect-1] > 4) histo(hTitlePhiTrigg.str())->Fill( phiglobal*radtodeg  ); // *** hPhiTrigg_MB1, ...
            }
	  } // endif phi < 30 deg
//       end trigger efficiency analysis ------------

/*
// interpolation analysis....-----------
  // get position at Mid by interpolating the position (not direction) of best
  // segment in Bot and Top to Mid surface
    int iwheel = 0;
    DTChamberId MB2Id(iwheel, ist, iisect);
  // Get segments in MB1 (if any)
    int ist1 = ist-1;
    if(ist == 1 || ist == 4 ) ist1 =2;    
    DTChamberId MB1Id(iwheel, ist1, iisect);
    DTRecSegment4DCollection::range segsMB1= segs->get(MB1Id);
    int nSegsMB1=segsMB1.second-segsMB1.first;
  // Get segments in MB3 (if any)
    int ist2 = ist+1;
    if(ist == 1 || ist == 4 ) ist2 =3;
    DTChamberId MB3Id(iwheel, ist2, iisect);
    DTRecSegment4DCollection::range segsMB3= segs->get(MB3Id);
    int nSegsMB3=segsMB3.second-segsMB3.first;
   // check if any segments in MB1 & MB3
    if (nSegsMB1 > 0 && nSegsMB3 > 0) {
 //     cout << " NsegMB1 : " << nSegsMB1 << " NsegMB3 : " << nSegsMB3 << endl;
      const DTRecSegment4D& bestMB1Seg= getBestSegment(segsMB1);
  //    cout << bestMB1Seg << endl;
      const DTRecSegment4D& bestMB3Seg= getBestSegment(segsMB3);
  //    cout << bestMB3Seg << endl;
      LocalPoint posAtMid = interpolate(bestMB3Seg, bestMB1Seg, MB2Id);
 //     cout << "PosAtMid " << posAtMid << endl;
      LocalPoint midSegPos= (*seg).localPosition();
      float distx = (midSegPos-posAtMid).x();
      float disty = (midSegPos-posAtMid).y();
      stringstream hTitlex;
      hTitlex << "halignX_S" << isect << "_MB" << ist ;
      histo(hTitlex.str())->Fill( distx );
      stringstream hTitlez;
      hTitlez << "halignZ_S" << isect << "_MB" << ist ;
      histo(hTitlez.str())->Fill( disty );
    } // endif segments in MBn-1 & MBn+1   
*/

// extrapolation analysis....-----------
  // get position at Mid by extrapolating the position + direction of best
  // segment in MB1
    int iwheel = 0;
    DTChamberId MB2Id(iwheel, ist, iisect);
  // Get segments in MB1 (if any) when looking to MB2/3/4
  // get segments in MB2 when looking to MB1
    int ist1 = 1;
    if(ist == 1  ) ist1 =2;    
    DTChamberId MB1Id(iwheel, ist1, iisect);
    DTRecSegment4DCollection::range segsMB1= segs->get(MB1Id);
    int nSegsMB1=segsMB1.second-segsMB1.first;
    if (nSegsMB1 > 0 ) {
 //     cout << " NsegMB1 : " << nSegsMB1  << endl;
      const DTRecSegment4D& bestMB1Seg= getBestSegment(segsMB1);
//    cout << bestMB1Seg << endl;
// extrapolate to MB2Id surface:
      LocalPoint posAtMid = extrapolate(bestMB1Seg, MB2Id);
 //     cout << "PosAtMid " << posAtMid << endl;
      LocalPoint midSegPos= (*seg).localPosition();
      float distx = (midSegPos-posAtMid).x();
      float disty = (midSegPos-posAtMid).y();
      stringstream hTitlex;
      hTitlex << "halignX_S" << isect << "_MB" << ist ;
      histo(hTitlex.str())->Fill( distx );
      stringstream hTitlez;
      hTitlez << "halignZ_S" << isect << "_MB" << ist ;
      histo(hTitlez.str())->Fill( disty );
    } // endif segments in MBn-1 & MBn+1   

// ----------------------------------------------------  



//  prepare vectors for re-fitting-----------------------
         float xfit[12], yfit [12];
         float sfit[12]; // flag (-1 or +1) for left/right semicell, input to fitline_t0 function 
//        for (unsigned int sphit = 0; sphit <  phiHits.size(); ++sphit){  // loop on segment hits
//            if(debug) cout <<(phiHits[sphit]).wireId() << endl;	   
//        } // next hit in 2D phi segment        
          DTSuperLayerId slid1(phiSeg->chamberId(),1);
          /// Mean timer analysis
          DTMeanTimer meanTimer1(dtGeom->superLayer(slid1), phiHits, eventSetup, theSync);
          vector<double> tMaxs1=meanTimer1.run();	
//        refit: prepare the input vectors	  
          DTChamberId chId = phiSeg->chamberId();
          const DTChamber* chamb = dtGeom->chamber(chId);
          int nhitfit = 0;  // nr. of hits in the re-fit
	  float tmax123[] = {0.,0.,0.,0.};
          for (int isl=1; isl<4; isl=isl+2) {  // loop on SL 1 & 3
            for (int il=1; il<5; il++) {        // loop on layers
              DTLayerId layId = DTLayerId(chId,isl,il);
              const DTLayer* lay = chamb->layer(layId);  
              GlobalPoint laypos =  lay->position();   
              LocalPoint laylocal = chamb->toLocal(laypos);
              double xextr = phiSeg->localPosition().x() +   // segment extrap. to this layer
     	                   ( phiSeg->localDirection().x()/phiSeg->localDirection().z() ) * laylocal.z();		     
              std::vector< DTRecHit1D > specificHits =  phiSeg->specificRecHits();
              for (unsigned int sphit = 0; sphit <  specificHits.size(); ++sphit){  // loop on segment hits
               if ((specificHits[sphit]).wireId().layerId().superlayerId().superLayer() == isl &&
                  (specificHits[sphit]).wireId().layer() == il ) {
                 float xhit = (specificHits[sphit]).localPosition().x() + laylocal.x() ; // localPosition is in the layer (not in the cell)
                 DTWireId wireId = (specificHits[sphit]).wireId();
//		 cout << " digiTime " << (specificHits[sphit]).digiTime() << " xlocal " << (specificHits[sphit]).localPosition().x() << endl;
                 float ttrig = theSync->offset(wireId);
//		 cout << "Sect " << isect << " MB" << ist << " SL" << isl << 
//		         " digiTime " << (specificHits[sphit]).digiTime() << " ttrigg " << ttrig << endl;
	         int ill = il;
	         if ( isl > 1) ill = il + 4;
	         if (il == 1 || il == 3 ) tmax123[isl] = tmax123[isl] + 0.5*( (specificHits[sphit]).digiTime() - ttrig );
                 if (il == 2)             tmax123[isl] = tmax123[isl] + (specificHits[sphit]).digiTime() - ttrig;
	         xfit[nhitfit] = laylocal.z(); // input to fitline function
	         yfit[nhitfit] = xhit;         // input to fitline function
                 sfit[nhitfit] = 1.;         // mark left/right semi-cell (1=Right) input to fitline_t0 function )
                 if( (specificHits[sphit]).lrSide() == 2 ) sfit[nhitfit] = -1.; // -1 = Left
	         nhitfit++; 
// ============== histograms of residuals ======================
  	         stringstream hTitle;
                 hTitle << "hResX_S" << isect << "_MB" << ist ;
                 histo(hTitle.str())->Fill( xhit-xextr); // *** histos: hResX_S1_MB1, ecc... **********		    
// ==============================================================
	       }  // endif hitlayer == layer 
             } // end loop on hits
           } // end loop on layers	      
         } // end loop on SL
// ======= tmax plots =======================================
          stringstream hTitle1, hTitle2, hTitle3;
          hTitle1 << "htmaxYB0_S" << isect << "_MB" << ist << "_SL1";
          histo(hTitle1.str())->Fill( tmax123[1]);
          hTitle2 << "htmaxYB0_S" << isect << "_MB" << ist << "_SL3";
          histo(hTitle2.str())->Fill( tmax123[3]);
//          hTitle3 << "hdtmaxMB" << ist << "_" << isect;
//          histo(hTitle3.str())->Fill( tmax123[3]-tmax123[1]);
// ===========================================================
	} // endif abs(phiLocal) < 35 deg	  
	} // endif phiHts.size() >= 6
      } // endif phiseg ---------------------------
      // zed segment
      const DTSLRecSegment2D* zedSeg =(*seg).zSegment();
      //cout << "zedSeg " << zedSeg << endl;
      vector<DTRecHit1D> zedHits;
      if (zedSeg) {
//        if(debug) cout << "Zed segm: " << *zedSeg << endl;  
        zedHits = zedSeg->specificRecHits();
        if(debug) cout << " Nhits in z-segment " << zedHits.size() << endl;
        NtkHit = NtkHit + zedHits.size();
//        for (unsigned int sphit = 0; sphit <  zedHits.size(); ++sphit){  // loop on segment hits
//            if(debug) cout <<(zedHits[sphit]).wireId() << endl;	   
//        } // next hit in 2D zed segment  
      } // endif zedSeg


// ==== fill  hit multiplicity histos per sector ===============
      stringstream hTitle;
      hTitle << "hNhits_S" << isect << "_MB" << ist ;
      histo(hTitle.str())->Fill( NtkHit ); // hNhits_S1_MB1, ....ecc
	  
      if (phiSeg && zedSeg) nSegSect[isect]++;
           
    } // end loop 4D segment in chamber ========================

  } // end loop on all chambers =========================


// ==== fill segment multiplicity histos per sector ===============
    for (int sec=1; sec<=12; ++sec) { // section 1 to 12
         stringstream hTitle;
          hTitle << "hNsegs_S" << sec;
         if( nSegSect[sec] > 0 )  histo(hTitle.str())->Fill( nSegSect[sec]); // hNsegs_Sect1, ....ecc
    }
// ================================================================ 

        
} // end DTsegments analyzer =============================================

  
// ==========================================================================================================
// StanAlone Track muon analysis
void DTAnalyzer::analyzeSATrack(const Event & event,
                                const EventSetup& eventSetup){
// ==========================================================================================================

  if (debug) cout << " analyzeSATracks : ------------------------------------------" << endl;

  MuonPatternRecoDumper muonDumper;

  // Get the DT Geometry
  ESHandle<DTGeometry> dtGeom;
  eventSetup.get<MuonGeometryRecord>().get(dtGeom);
//  const std::vector<DTChamber*> & chs = dtGeom->chambers();

  
  // Get the 4D rechit collection from the event -------------------
  edm::Handle<DTRecSegment4DCollection> segs;
  event.getByLabel(theRecHits4DLabel, segs);

  // Get the 2D rechit collection from the event -------------------
  edm::Handle<DTRecSegment2DCollection> segs2d;
  event.getByLabel(theRecHits2DLabel, segs2d);

  // Get the 1D rechits from the event --------------
  Handle<DTRecHitCollection> dtRecHits; 
  event.getByLabel(theRecHits1DLabel, dtRecHits);

  // Get the RecTrack collection from the event
  Handle<reco::TrackCollection> staTracks;
  event.getByLabel(theSTAMuonLabel, staTracks);

  ESHandle<MagneticField> theMGField;
  eventSetup.get<IdealMagneticFieldRecord>().get(theMGField);

  ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
  eventSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);


  reco::TrackCollection::const_iterator staTrack;

  double recPt=0.;
  histo("hNSA")->Fill(staTracks->size());
  histo2d("hNSAVsNHits")->Fill(dtRecHits->size(),staTracks->size());
  histo2d("hNSAVsNSegs2D")->Fill(segs2d->size(),staTracks->size());
  histo2d("hNSAVsNSegs4D")->Fill(segs->size(),staTracks->size());

  if (debug && staTracks->size() ) 
    cout << endl<<"Run/Event " << event.id().run() << ":" << event.id().event() << 
      " SA tracks " << staTracks->size() << endl;

// ********* loop on SA tracks ************************************************

  for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack){
    
    reco::TransientTrack track(*staTrack,&*theMGField,theTrackingGeometry); 

    if (debug) {
      cout << muonDumper.dumpFTS(track.impactPointTSCP().theState());

//      double posx = track.impactPointTSCP().position().x();
      recPt = track.impactPointTSCP().momentum().perp(); 
//      cout<<" x = " << posx << endl;   
      cout<<" p: "<<track.impactPointTSCP().momentum().mag()<< " pT: "<<recPt<<endl;
      cout<<" normalized chi2: "<<track.normalizedChi2()<<endl;
      //cout << " Z magentic Field: "<<track.field()->inTesla().z()<<endl;
    }

    double posx = track.impactPointTSCP().position().x();
    double posy = track.impactPointTSCP().position().y();
    double posz = track.impactPointTSCP().position().z();
    double radius = sqrt( posx*posx + posy*posy );

//  total hits in STA
    histo("hNhitsSA")->Fill(staTrack->recHitsSize());
    histo2d("hNHitsSAVsNHits")->Fill(dtRecHits->size(),staTrack->recHitsSize());
    histo2d("hNHitsSAVsNSegs2D")->Fill(segs2d->size(),staTrack->recHitsSize());
    histo2d("hNHitsSAVsNSegs4D")->Fill(segs->size(),staTrack->recHitsSize());

//  chi2...
    histo("hChi2SA")->Fill(staTrack->normalizedChi2());
    histo("hSAValidHits")->Fill(staTrack->numberOfValidHits());
    histo("hSALostHits")->Fill(staTrack->numberOfLostHits());

//   plots at Impact Point
    histo("hPIPSA")->Fill(track.impactPointTSCP().momentum().mag());
    histo("hPtIPSA")->Fill(recPt);
    histo("hPhiIPSA")->Fill(track.impactPointTSCP().momentum().phi());
    histo("hEtaIPSA")->Fill(track.impactPointTSCP().momentum().eta());
    histo("hrIPSA")->Fill(radius);
    histo("hzIPSA")->Fill(posz);
    histo2d("hrVszIPSA")->Fill(posz,radius);
    
    TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
    if (debug) {
      cout << "Inner TSOS:"<<endl;
      cout << muonDumper.dumpTSOS(innerTSOS);
      cout<<" p: "<<innerTSOS.globalMomentum().mag()<< " pT: "<<innerTSOS.globalMomentum().perp()<<endl;
    }

    histo("hPInnerTSOSSA")->Fill(innerTSOS.globalMomentum().mag());
    histo("hPtInnerTSOSSA")->Fill(innerTSOS.globalMomentum().perp());
    histo("hPhiInnerTSOSSA")->Fill(innerTSOS.globalMomentum().phi());
    histo("hEtaInnerTSOSSA")->Fill(innerTSOS.globalMomentum().eta());

    histo("hInnerRSA")->Fill(sqrt(staTrack->innerPosition().perp2()));
    histo("hOuterRSA")->Fill(sqrt(staTrack->outerPosition().perp2()));
    histo2d("hInnerOuterRSA")->Fill(sqrt(staTrack->innerPosition().perp2()),
                                    sqrt(staTrack->outerPosition().perp2()));


    if (debug) cout<<" Associated RecHits:"<<endl;

    trackingRecHit_iterator rhbegin = staTrack->recHitsBegin();
    trackingRecHit_iterator rhend = staTrack->recHitsEnd();

//  loop on associated rechits ----------------------------------------------------------
    int NHass = 0;
    int sect = 0;
    int sect1 = 0;
    int wheel1 = 0;
    int sect2 = 0;
    int wheel2 = 0;

    for(trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit){
      const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
      GlobalPoint gpos=geomDet->toGlobal((*recHit)->localPosition());
//    cout << (*recHit)->geographicalId() << endl;
//      const TrackingRecHit* trh = (*recHit)->recHits();
      NHass++;
      DTWireId wireId ;
      if ( const DTRecHit1D* dthit = dynamic_cast< const DTRecHit1D*> ((*recHit)->clone()) ) {
       wireId = dthit ->wireId();
/*       cout << "  Wire:   " << wireId << endl;
       cout << "  wheel   " << wireId.wheel() << endl;
       cout << "  sector  " << wireId.sector() << endl;
       cout << "  chamber " << wireId.station() << endl;
       cout << "  SL      " << wireId.superlayer() << endl;
       cout << "  layer   " << wireId.layer() << endl;
       cout << "  wire    " << wireId.wire() << endl; */
     }

      if (NHass == 1) sect = wireId.sector();
      //if (NHass == 1) {
      //sect1 = wireId.sector();
      //wheel1 = wireId.wheel();
      //}

      //if (NHass > 1) {
      //sect2 = wireId.sector();
      //wheel2 = wireId.wheel();

      //if ( sect2 == sect1 && wheel2 == wheel1) continue;
      //else if(sect2 != sect1 || wheel2 != wheel1) {
      //  NHass = 0;
      //  break; 

      //}
      //}

      //if(NHass != 0) {

      //histo("hChi2SA")->Fill(staTrack->normalizedChi2());
      //}


      if (debug)  cout<< wireId << " r: "<< gpos.perp() <<" x,y,z: "<< gpos.x() << " " << gpos.y() << " " << gpos.z() << " " << gpos.phi() << endl;
      histo2d("hHitsPosXYSA")->Fill(gpos.x(),gpos.y());
      histo2d("hHitsPosXYSA_1")->Fill(gpos.x(),gpos.y());
      histo2d("hHitsPosXYSA_2")->Fill(gpos.x(),gpos.y());
      histo2d("hHitsPosXYSA_3")->Fill(gpos.x(),gpos.y());
      histo2d("hHitsPosXYSA_4")->Fill(gpos.x(),gpos.y());
      histo2d("hHitsPosXZSA")->Fill(gpos.z(),gpos.x());
      histo2d("hHitsPosYZSA")->Fill(gpos.z(),gpos.y());
      float radtodeg = 57.296 ;
      float phi =  gpos.phi();
      stringstream hTitlePhihit;
      hTitlePhihit << "hPhiHit_MB" << wireId.station() ;
      if(wireId.superlayer() != 2) histo(hTitlePhihit.str())->Fill( radtodeg*phi); // hPhiHit_MB1, ....ecc                 
    } // next associated rechit --------------------------------
  
// ==== fill ass hit multiplicity histos per sector ===============
     stringstream hTitle;
     hTitle << "hNhass_S" << sect ;
     histo(hTitle.str())->Fill( NHass); // hNhass_S1_MB1, ....ecc


// ---- start extrapolation analysis ------------------
    // try to extrapolate using thePropagator ------
  if (!thePropagator){
    ESHandle<Propagator> prop;
    eventSetup.get<TrackingComponentsRecord>().get(thePropagatorName, prop);
    thePropagator = prop->clone();
    thePropagator->setPropagationDirection(anyDirection);
  }

    float radtodeg = 57.296 ;
    // Get a surface (here a cylinder of radius 1290mm) : ECAL
    float radiusECAL = 129.0; // radius in centimeter
    Cylinder::PositionType pos0;
    Cylinder::RotationType rot0;
    const Cylinder::CylinderPointer ecal = Cylinder::build(pos0, rot0,radiusECAL);
    //cout << "Cyl " << ecal->radius() << endl;
    TrajectoryStateOnSurface tsosAtEcal =
      thePropagator->propagate(*innerTSOS.freeState(), *ecal);
    if (tsosAtEcal.isValid()) {
//      cout << "extrap to ECAL : " << tsosAtEcal.globalPosition() << " r " <<
//        tsosAtEcal.globalPosition().perp() << endl;
	float phiEcal = -radtodeg*acos(tsosAtEcal.globalPosition().x()/tsosAtEcal.globalPosition().perp());
	float zEcal = tsosAtEcal.globalPosition().z();
        if (abs(zEcal) < 290. ) histo2d("hphivszECAL")->Fill(zEcal,phiEcal);
      }
//    else 
//      cout << "Extrapolation to ECAL failed" << endl;

    // Get a surface (here a cylinder of radius 1811 mm): HCAL
    float radiusHCAL = 181.1; // radius in centimeter
    const Cylinder::CylinderPointer hcal = Cylinder::build(pos0, rot0,radiusHCAL);
    //cout << "Cyl " << hcal->radius() << endl;
    TrajectoryStateOnSurface tsosAtHcal =
      thePropagator->propagate(*innerTSOS.freeState(), *hcal);
    if (tsosAtHcal.isValid()) {
//      cout << "extrap to HCAL : " << tsosAtHcal.globalPosition() << " r " <<
//        tsosAtHcal.globalPosition().perp() << endl;
	float phiHcal = -radtodeg*acos(tsosAtHcal.globalPosition().x()/tsosAtHcal.globalPosition().perp());
	float zHcal = tsosAtHcal.globalPosition().z();
        if (abs(zHcal) < 390. ) histo2d("hphivszHCAL")->Fill(zHcal,phiHcal);
       }
//    else 
//      cout << "Extrapolation to HCAL failed" << endl;
// --- end extraploation analysis -------------------------

        
  } // next SA track ============================================
  
  if (debug) cout<<"---  end event ------------------------------"<<endl;  
}

// ==========================================================================================================
// trigger analysis
void DTAnalyzer::analyzeTrigger(const Event & event,
                                const EventSetup& eventSetup){
// ==========================================================================================================

  if (debug) cout << " analyzeTrigger : ------------------------------------------" << endl;

  // Get the DT Geometry
  ESHandle<DTGeometry> dtGeom;
  eventSetup.get<MuonGeometryRecord>().get(dtGeom);
//  const std::vector<DTChamber*> & chs = dtGeom->chambers();

  // get the DT local trigger collection =======================
  edm::Handle<DTLocalTriggerCollection> allLocalTriggers;
  event.getByLabel("dtunpacker", allLocalTriggers); 

  bool hasTr[12][4];
  for (int ise=0;ise<12;ise++) for(int ist=0;ist<4;ist++) hasTr[ise][ist]=false;
  bool SecthasTr[]={false,false,false,false,false,false,false,false,false,false,false,false};

  int SCsect=0; int SCst=0;

  int bx[4][5]; 
  float bxbest[4][14];
//  float qual[]={-1,-1,-1,-1};  
  int qual[4][14];
  for (int sec=1; sec<15; ++sec) { // section 1 to 14
    for (int ist=0; ist<5; ++ist) { 
      qual[ist-1][sec-1] = -1;
      bxbest[ist-1][sec-1] = -1;
    }    
  } 
  
  DTLocalTriggerCollection::DigiRangeIterator chambIt;

  for (chambIt=allLocalTriggers->begin();chambIt!=allLocalTriggers->end();++chambIt){ // loop on chambers ------

    const DTChamberId& id = (*chambIt).first;
    SCst=id.station();
    SCsect=id.sector();

    bx[SCst-1][0]=0;     bx[SCst-1][1]=0;

    int countsec=-1;
    for (int ns=0; ns<12; ns++) {
     if (SCsect == sects[ns] ) {
       countsec=ns;
       break;
     }
     else if (sects[ns]==0) {
      sects[ns]=SCsect;
      countsec=ns;
      break;
     }
    }

    int ntrCh=0;

    const DTLocalTriggerCollection::Range& range = (*chambIt).second;
  
   // Loop over triggers of this chamber -----------------------------------------------
    for (DTLocalTriggerCollection::const_iterator trigtrack = range.first;trigtrack!=range.second;++trigtrack){
       if (trigtrack->quality()<7 && trigtrack->quality()>1) {
       SecthasTr[countsec]=true;    
       hasTr[SCsect-1][SCst-1]=true;
      }
      if (trigtrack->quality()<7) {
       bx[SCst-1][ntrCh]=trigtrack->bx();  // dice a che bx c'e' stato trigger  
       float ibx=trigtrack->bx();
       ntrCh++; 
       if (ntrCh>4) break;
          int iQual=trigtrack->quality();
          if( iQual > qual[SCst-1][SCsect-1] ) {
	   qual[SCst-1][SCsect-1]= iQual;  // store the quality for this station
	   bxbest[SCst-1][SCsect-1] = ibx;
	  } 
      } //  endif quality <7 
    } // end loop on triggers in this chamber --------------------------

//    plot the highest quality phi-trigger
      stringstream hTitleTrigg; 
      hTitleTrigg << "hTrigBits_S" << SCsect << "_MB" << SCst ;
      if(SCsect != 0) histo(hTitleTrigg.str())->Fill( qual[SCst-1][SCsect-1] ); // *** histos: hTrigBits_S1_MB1, ecc... **********
//    plot the BXid of highest quality
      int iqual = qual[SCst-1][SCsect-1];
      if(iqual < 4) iqual = 1; // flag for all uncorrelated phi trigger
      stringstream hTitleTriggBX; 
      hTitleTriggBX << "hTrigBX_S" << SCsect << "_MB" << SCst  << "_qual" << iqual;
      if(SCsect != 0) histo(hTitleTriggBX.str())->Fill( bxbest[SCst-1][SCsect-1] ); // *** histos: hTrigBX_S1_MB1_qual1, ecc... **********

  } // end loop on chambers ----------------------------------------------

  // Now let's Fill summary plots !

  int SectrigMFill=0;
  if (sects[4]) SectrigMFill=1;

  if (SecthasTr[0] && !SecthasTr[1] && !SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 3;
  if (!SecthasTr[0] && SecthasTr[1] && !SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 4;
  if (!SecthasTr[0] && !SecthasTr[1] && SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 5;
  if (!SecthasTr[0] && !SecthasTr[1] && !SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 6;


  if (SecthasTr[0] && SecthasTr[1] && !SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 8;
  if (SecthasTr[0] && !SecthasTr[1] && SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 9;
  if (SecthasTr[0] && !SecthasTr[1] && !SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 10;
  if (!SecthasTr[0] && SecthasTr[1] && SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 11;
  if (!SecthasTr[0] && SecthasTr[1] && !SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 12;
  if (!SecthasTr[0] && !SecthasTr[1] && SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 13;

  if (SecthasTr[0] && SecthasTr[1] && SecthasTr[2] && !SecthasTr[3] ) SectrigMFill = 15;
  if (SecthasTr[0] && SecthasTr[1] && !SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 16;
  if (SecthasTr[0] && !SecthasTr[1] && SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 17;
  if (!SecthasTr[0] && SecthasTr[1] && SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 18;

  if (SecthasTr[0] && SecthasTr[1] && SecthasTr[2] && SecthasTr[3] ) SectrigMFill = 20;

  histo("SectorTriggerMatrix")->Fill(SectrigMFill);

  for (int ns=0; ns<12; ns++) {

   int sector=sects[ns];

   if (!sector) continue;

   if (SecthasTr[ns]){

    for (int ist=0; ist<4; ist++){
     if (hasTr[sector-1][ist]) histo("TriggerInclusive")->Fill((sector-1)*5+ist+1);
    }
 
    int trigMFill=0;

    if (hasTr[sector-1][0] && !hasTr[sector-1][1] && !hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 2;
    if (!hasTr[sector-1][0] && hasTr[sector-1][1] && !hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 3;
    if (!hasTr[sector-1][0] && !hasTr[sector-1][1] && hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 4;
    if (!hasTr[sector-1][0] && !hasTr[sector-1][1] && !hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 5;
    if (hasTr[sector-1][0] && hasTr[sector-1][1] && !hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 7;
    if (hasTr[sector-1][0] && !hasTr[sector-1][1] && hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 8;
    if (hasTr[sector-1][0] && !hasTr[sector-1][1] && !hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 9;
    if (!hasTr[sector-1][0] && hasTr[sector-1][1] && hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 10;
    if (!hasTr[sector-1][0] && hasTr[sector-1][1] && !hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 11;
    if (!hasTr[sector-1][0] && !hasTr[sector-1][1] && hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 12;
    if (hasTr[sector-1][0] && hasTr[sector-1][1] && hasTr[sector-1][2] && !hasTr[sector-1][3]) trigMFill = 14;
    if (hasTr[sector-1][0] && hasTr[sector-1][1] && !hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 15;
    if (hasTr[sector-1][0] && !hasTr[sector-1][1] && hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 16;
    if (!hasTr[sector-1][0] && hasTr[sector-1][1] && hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 17;
    if (hasTr[sector-1][0] && hasTr[sector-1][1] && hasTr[sector-1][2] && hasTr[sector-1][3]) trigMFill = 19;

    stringstream hNameTrigMat; 
    hNameTrigMat << "TriggerMatrix" << sector ;
    histo(hNameTrigMat.str())->Fill(trigMFill);

   } // if the sector has trigger
  } // loop on sectors
}


//  utilities ###########################################

// interpolation between stations -----------------------------

LocalPoint DTAnalyzer::interpolate(const DTRecSegment4D& seg1,
                                      const DTRecSegment4D& seg3,
                                      const DTChamberId& id2) const {
  // Get GlobalPoition of Seg in MB1
  GlobalPoint
    gpos1=(dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());

  // Get GlobalPoition of Seg in MB3
  GlobalPoint
    gpos3=(dtGeom->chamber(seg3.chamberId()))->toGlobal(seg3.localPosition());

   
  // interpolate
  // get all in MB2 frame
  LocalPoint pos1=(dtGeom->chamber(id2))->toLocal(gpos1);
  LocalPoint pos3=(dtGeom->chamber(id2))->toLocal(gpos3);


  // case 1: 1 and 3 has both projection. No problem

  // case 2: one projection is missing for one of the segments. Keep the other's
  // segment position
  if (!seg1.hasZed()) pos1=LocalPoint(pos1.x(),pos3.y(),pos1.z());
  if (!seg3.hasZed()) pos3=LocalPoint(pos3.x(),pos1.y(),pos3.z());

  if (!seg1.hasPhi()) pos1=LocalPoint(pos3.x(),pos1.y(),pos1.z());
  if (!seg3.hasPhi()) pos3=LocalPoint(pos1.x(),pos3.y(),pos3.z());

  // cout << "pos1 " << pos1 << endl;
  // cout << "pos3 " << pos3 << endl;
  // direction
  LocalVector dir = (pos3-pos1).unit(); // z points inward!
  // cout << "dir " << dir << endl;
  LocalPoint pos2 = pos1+dir*pos1.z()/(-dir.z());
  // cout << "pos2 " << pos2 << endl;


  return pos2;
}

// extrapolation from one station to another ------------------

LocalPoint DTAnalyzer::extrapolate(const DTRecSegment4D& seg1, 
                                   const DTChamberId& id2) const {

  // Get GlobalPosition & direction of Seg in MB1
  GlobalPoint gpos1=(dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localPosition());
  GlobalVector gdir1=(dtGeom->chamber(seg1.chamberId()))->toGlobal(seg1.localDirection());

     // extrapolate....:
  // get all in MB2 frame
  LocalPoint     pos1=(dtGeom->chamber(id2))->toLocal(gpos1);
  LocalVector    dir1=(dtGeom->chamber(id2))->toLocal(gdir1);


  // direction
//  LocalVector dir = (pos3-pos1).unit(); // z points inward!
  // cout << "dir " << dir << endl;
//  LocalPoint pos2 = pos1+dir*pos1.z()/(-dir.z());
  LocalPoint pos2 = pos1+dir1*pos1.z()/(-dir1.z());
  // cout << "pos2 " << pos2 << endl;

  return pos2;
}

// best segment utilities --------------------------------------

// as usual max number of hits and min chi2
const DTRecSegment4D& DTAnalyzer::getBestSegment(const DTRecSegment4DCollection::range& segs) const{
  DTRecSegment4DCollection::const_iterator bestIter;
  unsigned int nHitBest=0;
  double chi2Best=99999.;
  for (DTRecSegment4DCollection::const_iterator seg=segs.first ;
       seg!=segs.second ; ++seg ) {
    unsigned int nHits= ((*seg).hasPhi() ? (*seg).phiSegment()->recHits().size() : 0 ) ;
    nHits+= ((*seg).hasZed() ?  (*seg).zSegment()->recHits().size() : 0 );

    if (nHits==nHitBest) {
      if ((*seg).chi2()/(*seg).degreesOfFreedom() < chi2Best ) {
        chi2Best=(*seg).chi2()/(*seg).degreesOfFreedom();
        bestIter = seg;
      }
    }
    else if (nHits>nHitBest) {
      nHitBest=nHits;
      bestIter = seg;
    }
  }
  return *bestIter;
}


const DTRecSegment4D* DTAnalyzer::getBestSegment(const DTRecSegment4D* s1,
                                                    const DTRecSegment4D* s2) const{

  if (!s1) return s2;
  if (!s2) return s1;
  unsigned int nHits1= (s1->hasPhi() ? s1->phiSegment()->recHits().size() : 0 ) ;
  nHits1+= (s1->hasZed() ?  s1->zSegment()->recHits().size() : 0 );

  unsigned int nHits2= (s2->hasPhi() ? s2->phiSegment()->recHits().size() : 0 ) ;
  nHits2+= (s2->hasZed() ?  s2->zSegment()->recHits().size() : 0 );

  if (nHits1==nHits2) {
    if (s1->chi2()/s1->degreesOfFreedom() < s2->chi2()/s2->degreesOfFreedom() ) 
      return s1;
    else
      return s2;
  }
  else if (nHits1>nHits2) return s1;
  return s2;
}


// ---------------------------------------

TH1F* DTAnalyzer::histo(const string& name) const{
  if (TH1F* h =  dynamic_cast<TH1F*>(theFile->Get(name.c_str())) ) return h;
  else throw cms::Exception("DTAnalyzer") << " Not a TH1F " << name;
}

TH2F* DTAnalyzer::histo2d(const string& name) const{
  if (TH2F* h =  dynamic_cast<TH2F*>(theFile->Get(name.c_str())) ) return h;
  else throw  cms::Exception("DTAnalyzer") << " Not a TH2F " << name;
}

bool DTAnalyzer::getLCT(LCTType t) const {
  return LCT.test(t);
}

string DTAnalyzer::toString(const DTLayerId& id) const {
  stringstream result;
  result << "_Wh" << id.wheel() 
    << "_Sec" << id.sector() 
    << "_St" << id.station()
    << "_Sl" << id.superLayer()
    << "_Lay" << id.layer();
  return result.str();
}

string DTAnalyzer::toString(const DTSuperLayerId& id) const {
  stringstream result;
  result << "_Wh" << id.wheel() 
    << "_Sec" << id.sector() 
    << "_St" << id.station()
    << "_Sl" << id.superLayer();
  return result.str();
}

string DTAnalyzer::toString(const DTChamberId& id) const {
  stringstream result;
  result << "_Wh" << id.wheel() 
    << "_Sec" << id.sector() 
    << "_St" << id.station();
  return result.str();
}

template<class T> string DTAnalyzer::hName(const string& s, const T& id) const {
  string name(toString(id));
  stringstream hName;
  hName << s << name;
  return hName.str();
}

void DTAnalyzer::createTH1F(const std::string& name,
                            const std::string& title,
                            const std::string& suffix,
                            int nbin,
                            const double& binMin,
                            const double& binMax) const {
  stringstream hName;
  stringstream hTitle;
  hName << name << suffix;
  hTitle << title << suffix;
  new TH1F(hName.str().c_str(), hTitle.str().c_str(), nbin,binMin,binMax);
}

void DTAnalyzer::createTH2F(const std::string& name,
                            const std::string& title,
                            const std::string& suffix,
                            int nBinX,
                            const double& binXMin,
                            const double& binXMax,
                            int nBinY,
                            const double& binYMin,
                            const double& binYMax) const {
  stringstream hName;
  stringstream hTitle;
  hName << name << suffix;
  hTitle << title << suffix;
  new TH2F(hName.str().c_str(), hTitle.str().c_str(), nBinX,binXMin,binXMax, nBinY,binYMin,binYMax);
}

    
//  end utilities ...#######################################################   

DEFINE_FWK_MODULE(DTAnalyzer);