/******* \class MyDTAnalyzer *******
 *
 * Description: Alignment geometry validation
 *  
 * Calculation and analysis of residuals of propagating 2D segments between 
 * SL inside each DT chamber and residuals of propagating 4D segments between 
 * consecutive chambers inside a sector.
 *
 * \author : A.Calderon <alicia.calderon.pd.infn.it> 26/10/2007
 * $date   : 12/11/2007  CET $
 *
 * Modification: 07 April 08
 *
 *********************************/

/* This Class Header */
#include "RecoLocalMuon/DTSegment/test/MyDTAnalyzer.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 "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
#include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.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"
#include "RecoLocalMuon/DTSegment/test/DTSegmentResidual.h"

/* This is to handle extrapolations and matching */
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.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/DTSLRecSegment2D.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"

#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "RecoLocalMuon/DTSegment/src/DTLinearFit.h"

/* MC analysis */
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

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

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

  // Get the debug parameter for verbose output
  debug = pset.getUntrackedParameter<bool>("debug");
  verbosity = pset.getParameter<int>("Verbosity");
  printout = pset.getParameter<bool>("PRINT");
  theRootFileName = pset.getUntrackedParameter<string>("rootFileName");

  // 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");

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

  //Is MC 
  isMC = pset.getParameter<bool>("isMC");
  theSimTrackLabel = pset.getParameter<edm::InputTag>("SimTrackLabel");

  //Cut in the Chi2 of the track 
  theChi2Cut = pset.getParameter<double>("chi2Cut");

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

  bool dirStat=TH1::AddDirectoryStatus();
  TH1::AddDirectory(kTRUE);


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

  // global histos...------------

  // segments
  new TH1F("hnSegDT","Num seg DT",50,0.,50.);
  new TH1F("hHitResidualSeg", "Hits residual wrt segments ",100,-.2,+.2);
  new TH1F("hGlobalPhi4DSegMB1","4D Segment Global Phi Angle MB1",100,0,-3);
  new TH1F("hGlobalPhi4DSegMB2","4D Segment Global Phi Angle MB2",100,0,-3);
  new TH1F("GlobalPhi4DSeg","4D Segment Global Phi Angle",100,0,-3);
  new TH1F("DiffPhi4DSegMB1MB2","Diff Global Phi Angle MB1-MB2",100,-.1,+.1);
  new TH1F("hChi24DSeg","Chi2 of 4D segments",100,0,100);
  new TH1F("hChi2NDof4DSeg","Chi2/NDof of 4D segments",100,0,100);

  // CosmicMuon reconstruction
  new TH1F("hNSA","Num SA tracks in events", 6, 0, 6);
  new TH1F("hNSegs4D","Num Seg4D in events", 12, 0, 12);
  new TH1F("hNPhiSegs","Num PhiSegs in events", 20, 0, 20);
  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("hProbChi2SA","Probability Chi2 for SA tracks", 100, 0, 1.);
  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","Num 1d Hits vs N SA", 100,-0.5,39.5, 5,-0.5, 4.5);
  new TH2F("hNSAVsNSegs2D","Num 2d Segs vs N SA", 20,-0.5,19.5, 5,-0.5, 4.5);
  new TH2F("hNSAVsNSegs4D","Num 4d Segs vs N SA", 10,-0.5,9.5, 5,-0.5, 4.5);

  new TH2F("hNHitsSAVsNHits","Num 1d Hits vs N Hits SA", 40,-0.5,39.5, 100,0,100);
  new TH2F("hNHitsSAVsNSegs2D","Num 2d Segs vs N Hits SA", 20,-0.5,19.5,100,0,100);
  new TH2F("hNHitsSAVsNSegs4D","Num 4d Segs vs N Hits SA", 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);
  
  // plots per Wheel/Sector/CHambers / SL ------------------

  for (int iwheel = 0; iwheel < 3; iwheel++) { // Wheel loop 
    for ( int sect = 1; sect < 15; sect++) {   // Sector loop
      for  ( int jst = 1; jst < 5; jst++ ) {   // station loop

  //Propagation residuals from A. Calderon
	//RPhi plane
	
	if (jst < 4) {
	  
	  stringstream hNamePropRPhi; 
	  stringstream hTitlePropRPhi;
	  hNamePropRPhi  << "hresRPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitlePropRPhi << "RPhi plane (LocalX-LocalPropX),W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1) ;
	  createTH1F(hNamePropRPhi.str(),hTitlePropRPhi.str(),"",100,-3.0, 3.0); // *** histos: hresRPhi_W0_S1_MB12, ecc... **********

	  //RZ plane
	  //if(jst != 4){
	  stringstream hNamePropZ; 
	  stringstream hTitlePropZ;
	  hNamePropZ  << "hresZ_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitlePropRPhi << "RZ plane (LocalY-LocalPropY),W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH1F(hNamePropZ.str(),hTitlePropZ.str(),"",100,-10, 10); // *** histos: hresZ_W0_S1_MB12, ecc... **********

	  //RPhi angles
	  stringstream hNameRPhi; 
	  stringstream hTitleRPhi;
	  hNameRPhi  << "hRPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleRPhi << "RPhi angle,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH1F(hNameRPhi.str(),hTitleRPhi.str(),"",100,-0.1, 0.1); // *** histos: hRPhi_W0_S1_MB12, ecc... **********

	  //RPhi angles vs Phi angle
	  stringstream hNameRPhivsPhi; 
	  stringstream hTitleRPhivsPhi;
	  hNameRPhivsPhi  << "hRPhivsPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleRPhivsPhi << "RPhi angle vs Phi,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH2F(hNameRPhivsPhi.str(),hTitleRPhivsPhi.str(),"",100,-3.0,3.0,100,-0.1, 0.1); // *** histos: hRPhi_W0_S1_MB12, ecc... **********

	  //X Resiuals vs Phi angle
	  stringstream hNameXvsPhi; 
	  stringstream hTitleXvsPhi;
	  hNameXvsPhi  << "hresXvsPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleXvsPhi << "Res X vs Phi ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH2F(hNameXvsPhi.str(),hTitleXvsPhi.str(),"",100,-3.0,3.0,100,-3.0,3.0); // *** histos: hresXvsPhi_W0_S1_MB12, ecc... **********

	  //X Resiuals / Phi angle
	  stringstream hNameXPhi; 
	  stringstream hTitleXPhi;
	  hNameXPhi  << "hresXdivPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleXPhi << "Res X / Phi ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH1F(hNameXPhi.str(),hTitleXPhi.str(),"",100,-3.0,3.0); // *** histos: hresXvsPhi_W0_S1_MB12, ecc... **********

	  //Z Resiuals vs Eta angle
	  stringstream hNameZvsEta; 
	  stringstream hTitleZvsEta;
	  hNameZvsEta  << "hresZvsEta_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleZvsEta << "Res Z vs Eta ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH2F(hNameZvsEta.str(),hTitleZvsEta.str(),"",100,-10.0,10.0,100,-3.0,3.0); // *** histos: hresZvsEta_W0_S1_MB12, ecc... **********

	  //X Resiuals vs Chi2
	  stringstream hNameXvsChi2; 
	  stringstream hTitleXvsChi2;
	  hNameXvsChi2  << "hresXvsChi2_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleXvsChi2 << "Res X vs Chi2 ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH2F(hNameXvsChi2.str(),hTitleXvsChi2.str(),"",100,0,100,100,-3.0,3.0); // *** histos: hresXvsPhi_W0_S1_MB12, ecc... **********

	  //Chi2 vs Phi
	  stringstream hNameChi2vsPhi; 
	  stringstream hTitleChi2vsPhi;
	  hNameChi2vsPhi  << "hresChi2vsPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst+1);
	  hTitleChi2vsPhi << "Chi2 vs Phi ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst+1);
	  createTH2F(hNameChi2vsPhi.str(),hTitleChi2vsPhi.str(),"",100,-3.0,3.0,100,0,100); // *** histos: hresXvsPhi_W0_S1_MB12, ecc... **********

	}

	if (jst > 1) {
	  
	  stringstream hNamePropRPhi; 
	  stringstream hTitlePropRPhi;
	  hNamePropRPhi  << "hresRPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst-1);
	  hTitlePropRPhi << "RPhi plane (LocalX-LocalPropX),W" << iwheel << "/S" << sect << "/MB" << jst << (jst-1) ;
	  createTH1F(hNamePropRPhi.str(),hTitlePropRPhi.str(),"",100,-3.0, 3.0); // *** histos: hresRPhi_W0_S1_MB12, ecc... **********
 
	  //RZ plane
	  //if(jst != 4){
	  stringstream hNamePropZ; 
	  stringstream hTitlePropZ;
	  hNamePropZ  << "hresZ_W" << iwheel << "_S" << sect << "_MB" << jst << (jst-1);
	  hTitlePropRPhi << "RZ plane (LocalY-LocalPropY),W" << iwheel << "/S" << sect << "/MB" << jst << (jst-1);
	  createTH1F(hNamePropZ.str(),hTitlePropZ.str(),"",100,-10, 10); // *** histos: hresz_W0_S1_MB12, ecc... **********

	  //RPhi angles
	  stringstream hNameRPhi; 
	  stringstream hTitleRPhi;
	  hNameRPhi  << "hRPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst-1);
	  hTitleRPhi << "RPhi angle,W" << iwheel << "/S" << sect << "/MB" << jst << (jst-1);
	  createTH1F(hNameRPhi.str(),hTitleRPhi.str(),"",100,-0.1, 0.1); // *** histos: hRPhi_W0_S1_MB12, ecc... **********

	  //X Resiuals vs Phi angle
	  stringstream hNameXvsPhi; 
	  stringstream hTitleXvsPhi;
	  hNameXvsPhi  << "hresXvsPhi_W" << iwheel << "_S" << sect << "_MB" << jst << (jst-1);
	  hTitleXvsPhi << "Res X vs Phi ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst-1);
	  createTH2F(hNameXvsPhi.str(),hTitleXvsPhi.str(),"",100,-3,3.0,100,-3.0,3.0); // *** histos: hresXvsPhi_W0_S1_MB12, ecc... **********

	  //Z Resiuals vs Eta angle
	  stringstream hNameZvsEta; 
	  stringstream hTitleZvsEta;
	  hNameZvsEta  << "hresZvsEta_W" << iwheel << "_S" << sect << "_MB" << jst << (jst-1);
	  hTitleZvsEta << "Res Z vs Eta ,W" << iwheel << "/S" << sect << "/MB" << jst << (jst-1);
	  createTH2F(hNameZvsEta.str(),hTitleZvsEta.str(),"",100,-10.0,10.0,100,-3.0,3.0); // *** histos: hresZvsEta_W0_S1_MB12, ecc... **********
	}

	//RPhi 2D segment X  between each SL
	stringstream hNameXSL; 
	stringstream hTitleXSL;
	hNameXSL  << "hXSL_W" << iwheel << "_S" << sect << "_MB" << jst ;
	hTitleXSL << "XSL X residuals,W" << iwheel << "/S" << sect << "/MB" << jst ;
	createTH1F(hNameXSL.str(),hTitleXSL.str(),"",100,-3, 3); // *** histos: hXSL_W0_S1_MB1, ecc... **********
      
	//RPhi 2D segment angles betwwwn each SL	
	stringstream hNameRPhiSL; 
	stringstream hTitleRPhiSL;
	hNameRPhiSL  << "hRPhiSL_W" << iwheel << "_S" << sect << "_MB" << jst ;
	hTitleRPhiSL << "RPhiSL angle,W" << iwheel << "/S" << sect << "/MB" << jst ;
	createTH1F(hNameRPhiSL.str(),hTitleRPhiSL.str(),"",100,-0.4, 0.4);// *** histos: hRPhiSL_W0_S1_MB1, ecc... **********

	//SLX Resiuals vs Phi angle
	stringstream hNameXSLvsPhi; 
	stringstream hTitleXSLvsPhi;
	hNameXSLvsPhi  << "hresXSLvsPhi_W" << iwheel << "_S" << sect << "_MB" << jst ;
	hTitleXSLvsPhi << "Res XSL vs Phi ,W" << iwheel << "/S" << sect << "/MB" << jst ;
	createTH2F(hNameXSLvsPhi.str(),hTitleXSLvsPhi.str(),"",100,-3.0,3.0,100,-3.0,3.0); // *** histos: hresXSLvsPhi_W0_S1_MB1, ecc... *********


	//X residuals Hit - Reco Segment vs Segment Phi angle
	stringstream hNameresHitvsPhi; 
	stringstream hTitleresHitvsPhi;
	hNameresHitvsPhi  << "hresHitvsPhi_W" << iwheel << "_S" << sect << "_MB" << jst ;
	hTitleresHitvsPhi << "Res Hit vs Phi,W" << iwheel << "/S" << sect << "/MB" << jst ;
	createTH2F(hNameresHitvsPhi.str(),hTitleresHitvsPhi.str(),"",100,-3,3,100,-1, 1); // *** histos: hresHitvsPhi_W0_S1_MB1, ecc... **********

	//X residuals Hit - Reco Segment vs X Hit
	stringstream hNameresHitvsX; 
	stringstream hTitleresHitvsX;
	hNameresHitvsX  << "hresHitvsX_W" << iwheel << "_S" << sect << "_MB" << jst ;
	hTitleresHitvsX << "Res Hit vs X,W" << iwheel << "/S" << sect << "/MB" << jst ;
	createTH2F(hNameresHitvsX.str(),hTitleresHitvsX.str(),"",100,-250,250,100,-1, 1); // *** histos: hresHitsvsX_W0_S1_MB1, ecc... *********

	//X residuals Seg - Propagated SA Track vs Segment Phi angle
	stringstream hNameresXSegSAvsPhi; 
	stringstream hTitleresXSegSAvsPhi;
	hNameresXSegSAvsPhi  << "hresXSegSAvsPhi_W" << iwheel << "_S" << sect << "_MB" << jst ;
	hTitleresXSegSAvsPhi << "X Seg-SA vs Phi,W" << iwheel << "/S" << sect << "/MB" << jst ;
	createTH2F(hNameresXSegSAvsPhi.str(),hTitleresXSegSAvsPhi.str(),"",100,-3,3,100,-3, 3); // *** histos: hresXSegSAvsPhi_W0_S1_MB1, ecc... **********

	/* // mis-alignment histos from Lacaprara
	   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 ===========================
    } // next sector ============================
  } //next wheel
  
  TH1::AddDirectory(dirStat); 

}

/* Destructor */ 
MyDTAnalyzer::~MyDTAnalyzer() {
  theFile->cd();
  theFile->Write();
  theFile->Close();
}


/* Operations */ 
// =================================================================
void MyDTAnalyzer::analyze(const Event & event,
                         const EventSetup& eventSetup) {
// =================================================================

  int nprt = event.id().event() - 100*(event.id().event()/100);

  //  cout << event.id().event() << "--" << 100*(event.id().event()/100) << endl;

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

  if (debug) cout << endl<<"--- [MyDTAnalyzer] Event analysed #Run: " <<
    event.id().run() << " #Event: " << event.id().event() << " ===========================" << endl;

  if (doHits) analyzeHits(event, eventSetup);
  if (doSegs) analyzeSegs(event, eventSetup);
  if (doSA) analyzeSA(event, eventSetup);

}

// =================================================================
void MyDTAnalyzer::analyzeHits(const Event & event,const EventSetup& eventSetup){
// =================================================================

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

  // 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 magnetic field
  ESHandle<MagneticField> theMGField;
  eventSetup.get<IdealMagneticFieldRecord>().get(theMGField);

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


  //=== MC analysis data
  if (isMC) {
    // Get the SimHit collection from the event
    Handle<PSimHitContainer> simHits;
    event.getByLabel(theSimTrackLabel, simHits);
  
    if(debug)
      cout << "   #SimHits: " << simHits->size() << endl;

    //-- Compare X SimHit position vs X RecoHit position
  }

  // -- Residual analysis per hit and station

  // loop on chambers ==============================================
  for (std::vector<DTChamber*>::const_iterator ch = chs.begin();
       ch!=chs.end() ; ++ch) {

    DTChamberId chid((*ch)->id());
    DTRecSegment4DCollection::range segsch= segs->get(chid);
     

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

    for (DTRecSegment4DCollection::const_iterator seg=segsch.first ;
         seg!=segsch.second ; ++seg ) {
      
      const DTChamber* ch = dtGeom->chamber(seg->chamberId());
      int chamber = (*ch).id().station();      
      int sector = (*ch).id().sector();
      int wheel = (*ch).id().wheel();    

      const DTChamberRecSegment2D* phiSeg= (*seg).phiSegment();

      LocalPoint lPos = (*seg).localPosition();
      LocalVector lDir = (*seg).localDirection();
      GlobalVector gDir = ch->toGlobal(lDir);
      
      if (phiSeg) {
	
	DTSegmentResidual res(phiSeg, ch);
	res.run();
	vector<DTSegmentResidual::DTResidual> deltas=res.residuals();

	//Loop over all DTResiduals vectors.
	for (vector<DTSegmentResidual::DTResidual>::const_iterator delta=deltas.begin();delta!=deltas.end(); ++delta) {
	  
	  // X Hit residuals vs Phi angle
	  stringstream hTitleresHitvsPhi;
	  hTitleresHitvsPhi << "hresHitvsPhi_W" << wheel << "_S" << sector << "_MB" << chamber;
	  histo2d(hTitleresHitvsPhi.str())->Fill(gDir.phi(),(*delta).value);	

	  // X Hit residuals vs Phi angle
	  stringstream hTitleresHitvsX;
	  hTitleresHitvsX << "hresHitvsX_W" << wheel << "_S" << sector << "_MB" << chamber;
	  histo2d(hTitleresHitvsX.str())->Fill(lPos.x(),(*delta).value);	
	}
      } // Is phiSeg
    } // End loop on 4D segments
  } // End loop on chambers
}


// =================================================================
void MyDTAnalyzer::analyzeSegs(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);
  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 magnetic field
  ESHandle<MagneticField> theMGField;
  eventSetup.get<IdealMagneticFieldRecord>().get(theMGField);

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

 
  //================================================================================================//
  //---------------  Residual 2D and 4D analysis/Alignment geometry validation ---------------------//
  //================================================================================================//
 

  //=====> Debug for printing out SL positions
  
  if (printout == true) {

    for (std::vector<DTChamber*>::const_iterator ch = chs.begin(); ch!=chs.end() ; ++ch) {
 
      DTChamberId chid((*ch)->id());

      if ( chid.wheel() == 0 && chid.sector() == 3 && chid.station() == 1) {

	cout << (*ch)->id() << endl;

	for (int il=1; il<5; il++) {        // loop on layers
	  DTLayerId layIdU = DTLayerId(chid,1,il);
	  const DTLayer* layU = (*ch)->layer(layIdU);  
	  GlobalPoint layposU =  layU->position();   
	  LocalPoint laylocalU = (*ch)->toLocal(layposU);
	  
	  cout << chid.sector() << " " << chid.station() << "  " << layIdU.superLayer() << "  " << layIdU.layer() << "  " << laylocalU.z() << endl;
	 
	}
		
	for (int il=1; il<5; il++) {        // loop on layers
	  DTLayerId layIdD = DTLayerId(chid,3,il);
	  const DTLayer* layD = (*ch)->layer(layIdD);  
	  GlobalPoint layposD =  layD->position();   
	  LocalPoint laylocalD = (*ch)->toLocal(layposD);
	  cout << chid.sector() << " " << chid.station() << "  " << layIdD.superLayer() << "  " <<  layIdD.layer() << "  " << laylocalD.z() << endl;
		
	}    
      }

      if ( chid.wheel() == 0 && chid.sector() == 1) {

	cout << chid.sector() << " " << chid.station() << "  " << (*ch)->position().x() << "   "  << (*ch)->position().y() << "   "  << (*ch)->position().z() << endl;
		
      }
    }
	    
    printout = false;
  }
 

  //=====> Residuals extrapolating 2D segments between SL inside each chamber

  // ********** loop on 2D segments *********************************
  
  const std::vector<DTSuperLayer*> & sls = dtGeom->superLayers();
  
  for (std::vector<DTSuperLayer*>::const_iterator sl = sls.begin(); sl!=sls.end() ; ++sl) {
    
    DTSuperLayerId slid((*sl)->id());
    DTRecSegment2DCollection::range segsSL= segs2d->get(slid);  
    
    const GeomDet* geomDet = theTrackingGeometry->idToDet((*sl)->geographicalId());

    for (int iwheel = 0; iwheel < 3; iwheel++) {

      DTChamberId chId  = slid.chamberId();
      const DTChamber* chamb = dtGeom->chamber(chId);
      int ist  = slid.chamberId().station();
      int isect  = slid.chamberId().sector();
      int iisect = isect;

      //RPhi SL (U == upper, D == down)
      int slU = 1;
      int slD = 3;

      DTSuperLayerId sl1Id(iwheel, ist, iisect, slU); 
      DTSuperLayerId sl2Id(iwheel, ist, iisect, slD);
     
      const DTRecSegment2DCollection::range segsSL1= segs2d->get(sl1Id);  
      const DTRecSegment2DCollection::range segsSL2= segs2d->get(sl2Id);  
	 
      int nSegsSL1 = segsSL1.second - segsSL1.first;
      int nSegsSL2 = segsSL2.second - segsSL2.first;
	 
      if ( nSegsSL1 > 0 &&  nSegsSL2 > 0) {

	const DTRecSegment2D& bestSL1Seg= getBestSegment(segsSL1);
	      
	const GeomDet* geomDet1 = theTrackingGeometry->idToDet((bestSL1Seg).geographicalId());

	if ( const DTSLRecSegment2D* dtSLseg1 = dynamic_cast< const DTSLRecSegment2D*> ((bestSL1Seg).clone()) ) {
	  
	  for (DTRecSegment2DCollection::const_iterator segSL2=segsSL2.first; segSL2!=segsSL2.second ; ++segSL2 ) {

	    const GeomDet* geomDet2 = theTrackingGeometry->idToDet((*segSL2).geographicalId());
	    unsigned int nHits1= ((bestSL1Seg).recHits()).size();
	    unsigned int nHits2= ((*segSL2).recHits()).size();
		
	    if( nHits1 == 4 && nHits2 == 4 ) { // Only take 2D segments with 4 hits each  

	      if ( const DTSLRecSegment2D* dtSLseg2 = dynamic_cast< const DTSLRecSegment2D*> ((*segSL2).clone()) ) {

		DTSuperLayerId SL2Id = dtSLseg2->superLayerId();

		//-- extrapolate to SL2Id surface:
		LocalPoint posAtMid = extrapolate(*dtSLseg1,SL2Id);		  
		LocalVector dirAtMid = bestSL1Seg.localDirection();
		//To global
		GlobalVector gdirAtMid = geomDet1->toGlobal(dirAtMid);
		//To local geomDet2
		LocalVector  ldirAtMid = geomDet2->toLocal(gdirAtMid);

		//In local gomDet2
		LocalPoint posSL2 = (*segSL2).localPosition();
		LocalVector ldirSL2 = (*segSL2).localDirection();
		// To Global
		GlobalVector gdirSL2 = geomDet2->toGlobal(ldirSL2);
	
		Double_t lPhi2 = atan(dirAtMid.z()/dirAtMid.x());

		//   if ( 1.4 < abs(lPhi2) && abs(lPhi2) < 1.7) { // Only perpendicular tracks to local chamber
		// Res in X position
		double diffX = posAtMid.x() - posSL2.x();
		
		stringstream hTitleXSL;
		hTitleXSL << "hXSL_W" << iwheel << "_S" << iisect << "_MB" << ist ;
		histo(hTitleXSL.str())->Fill( diffX );
		  
		// RPhi angle between SL 
		double diffRPhi = (atan(ldirSL2.z()/ldirSL2.x())) - (atan(ldirAtMid.z()/ldirAtMid.x()));
		stringstream hTitleRPhiSL;
		hTitleRPhiSL << "hRPhiSL_W" << iwheel << "_S" << iisect << "_MB" << ist ;
		histo(hTitleRPhiSL.str())->Fill( diffRPhi );
		 	
		// X residuaks vs Phi angle
		stringstream hTitleXSLvsPhi;
		hTitleXSLvsPhi << "hresXSLvsPhi_W" << iwheel << "_S" << iisect << "_MB" << ist ;
		histo2d(hTitleXSLvsPhi.str())->Fill(1/tan(gdirAtMid.phi()),diffX);
		  	
		//	  } // end of cut on Phi angle
	      }
	    }
	  }
	}
      }
    } // end loop on wheel
  }
  
	
  //=====> Residuals extrapolating 4D segments between consecutive chambers 
  //       MB2 to MB1, MB3 to MB2 and MB4 to MB3
 
  // loop on chambers ==============================================
  
  for (std::vector<DTChamber*>::const_iterator ch = chs.begin(); ch!=chs.end() ; ++ch) {
 
    DTChamberId chid((*ch)->id());
    //int chSector = chid.sector();
    //int chStation = 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;
 
      DTChamberId mCham((*seg).geographicalId().rawId());	
      //int station = mCham.station();

      const GeomDet* geomDet = theTrackingGeometry->idToDet((*seg).geographicalId());

      // extrapolation analysis....-----------
  
      for (int iwheel = 0; iwheel < 3; iwheel++) {
	//  int iwheel = 0;

	DTChamberId MB1Id(iwheel, ist, iisect);
     
	int ist1 = 1;
	int ist2 = 1; 
	
	if(ist == 1 || ist == 2 || ist == 3) {  
	  if(ist == 1) ist1 = 2;    
	  if(ist == 2) ist1 = 3;
	  if(ist == 3) ist1 = 4;

	  DTChamberId MB2Id(iwheel, ist1, iisect);    
	  
	  DTRecSegment4DCollection::range segsMB2= segs->get(MB2Id);
	  int nSegsMB2=segsMB2.second-segsMB2.first;
	  
	  if ( nSegsMB2 > 0 ) {

	    const DTRecSegment4D& bestMB2Seg= getBestSegment(segsMB2);
	  
	    const GeomDet* geomDet2 = theTrackingGeometry->idToDet((bestMB2Seg).geographicalId());
	    
	    const DTChamberRecSegment2D* phiSegMB2 = (bestMB2Seg).phiSegment();

	    //-- extrapolate to MB1Id surface:
	    LocalPoint posAtMid = extrapolate(bestMB2Seg, MB1Id);
	    LocalVector dirAtMid = bestMB2Seg.localDirection();
	    GlobalVector gdirAtMid = geomDet2->toGlobal(dirAtMid);
   
	    LocalPoint midSegPos= (*seg).localPosition();
	    LocalVector midSegDir = (*seg).localDirection();
	    GlobalVector gmidSegDir = geomDet->toGlobal(midSegDir);

	    const DTChamberRecSegment2D* phiSegMB1= (*seg).phiSegment();
	    
	    if (phiSegMB1 && phiSegMB2) {

	      // Cut on the number of hits 
	      unsigned int nHits1= phiSegMB1->specificRecHits().size();
	      unsigned int nHits2= phiSegMB2->specificRecHits().size();
		

	      // if( nHits1 >= 6  && nHits2 >= 6 ) { // Only take 4D segments with > 6 hits each  	    

	      LocalVector ldirMB1 = phiSegMB1->localDirection();
	      LocalVector ldirMB2 = phiSegMB2->localDirection();
	      
	      GlobalVector gdirMB1 = geomDet->toGlobal(ldirMB1);
	      GlobalVector gdirMB2 = geomDet2->toGlobal(ldirMB2);
	      
	      Double_t lPhi2 = atan(ldirMB2.z()/ldirMB2.x());
	      Double_t chi2 = bestMB2Seg.chi2();
	      Double_t nDof = bestMB2Seg.degreesOfFreedom();
	      Double_t chi2NDof = chi2/nDof;

	      histo("hChi24DSeg")->Fill(chi2);
	      histo("hChi2NDof4DSeg")->Fill(chi2NDof);

	      // if(chi2 < 15) { //Cut on Chi2 
		
	      //if ( 1.4 < abs(gdirMB2.phi()) && abs(gdirMB2.phi()) < 1.7) { // Only perpendicular tracks to the detector 

	      //if ( 1.4 < abs(lPhi2) && abs(lPhi2) < 1.7) { // Only perpendicular tracks to local chamber

	      // RPhi angle     
	      double diffRPhi = gdirMB1.phi() - gdirMB2.phi();
	      //double diffRPhi = (atan(ldirMB1.z()/ldirMB1.x())) - (atan(ldirMB2.z()/ldirMB2.x()));
	      stringstream hTitleAngle;
	      hTitleAngle << "hRPhi_W" << iwheel << "_S" << iisect << "_MB" << ist << (ist+1);
	      histo(hTitleAngle.str())->Fill( diffRPhi );

	      // RPhi angle vs Phi
	      stringstream hTitleAnglevsPhi;
	      hTitleAnglevsPhi << "hRPhivsPhi_W" << iwheel << "_S" << iisect << "_MB" << ist << (ist+1);
	      histo2d(hTitleAnglevsPhi.str())->Fill(gdirAtMid.phi(), diffRPhi );

	      // RPhi residuals
	      float distRPhi = (midSegPos-posAtMid).x();	      
	      stringstream hTitleRPhi;
	      hTitleRPhi << "hresRPhi_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo(hTitleRPhi.str())->Fill( distRPhi );

	      // RZ residuals
	      float distZ = (midSegPos-posAtMid).y();
	      stringstream hTitlez;
	      hTitlez << "hresZ_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo(hTitlez.str())->Fill( distZ );

	      //X residuals vs Phi angle
	      stringstream hTitleXvsPhi;
	      hTitleXvsPhi << "hresXvsPhi_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo2d(hTitleXvsPhi.str())->Fill(1/tan(gdirAtMid.phi()), distRPhi);
      
	      //X residuals/ Phi angle
	      stringstream hTitleXPhi;
	      hTitleXPhi << "hresXdivPhi_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo(hTitleXPhi.str())->Fill(tan(gdirAtMid.phi())*(distRPhi));
      

	      //Z residuals vs Eta angle
	      stringstream hTitleZvsEta;
	      hTitleZvsEta << "hresZvsEta_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo2d(hTitleZvsEta.str())->Fill(gdirAtMid.eta(), distZ);

	      //X residuals vs Chi2
	      stringstream hTitleXvsChi2;
	      hTitleXvsChi2 << "hresXvsChi2_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo2d(hTitleXvsChi2.str())->Fill(chi2,distRPhi);
 
	      //Chi2 vs Phi	   
	      stringstream hTitleChi2vsPhi;
	      hTitleChi2vsPhi << "hresChi2vsPhi_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist+1);
	      histo2d(hTitleChi2vsPhi.str())->Fill(gdirAtMid.phi(),chi2);
	      
	      //} // end of cut on Phi angle
	      //} // Cut on Chi2
	      //}// End cut on # recHits
	    }
	  }
	}// endif segments in MBn-1 & MBn+1   
    

	//=====> Residuals extrapolating 4D segments between consecutive chambers 
	//       MB1 to MB2, MB2 to MB3 and MB3 to MB4
 
	if(ist == 2 || ist == 3 || ist == 4) {
	  if(ist == 2) ist1 = 1;    
	  if(ist == 3) ist1 = 2;
	  if(ist == 4) ist1 = 3;

	  DTChamberId MB2Id(iwheel, ist1, iisect);    
	  
	  DTRecSegment4DCollection::range segsMB2= segs->get(MB2Id);
	  int nSegsMB2=segsMB2.second-segsMB2.first;
    
	  if ( nSegsMB2 > 0 ) {
	    //cout << " NsegMB2 : " << nSegsMB2  << endl;
	    const DTRecSegment4D& bestMB2Seg= getBestSegment(segsMB2);
	    
	    const GeomDet* geomDet2 = theTrackingGeometry->idToDet((bestMB2Seg).geographicalId());
	    
	    const DTChamberRecSegment2D* phiSegMB2 = (bestMB2Seg).phiSegment();
  
	    //    cout << bestMB2Seg << endl;

	    //-- extrapolate to MB1Id surface:
	    LocalPoint posAtMid = extrapolate(bestMB2Seg, MB1Id);
	    LocalVector dirAtMid = bestMB2Seg.localDirection();
	    GlobalVector gdirAtMid = geomDet2->toGlobal(dirAtMid);

	    //     cout << "PosAtMid " << posAtMid << endl;
	    
	    LocalPoint midSegPos= (*seg).localPosition();
	    LocalVector midSegDir = (*seg).localDirection();
	    GlobalVector gmidSegDir = geomDet->toGlobal(midSegDir);
	    
	    const DTChamberRecSegment2D* phiSegMB1= (*seg).phiSegment();
	    
	    if (phiSegMB1 && phiSegMB2) {
	
	      LocalVector ldirMB1 = phiSegMB1->localDirection();
	      LocalVector ldirMB2 = phiSegMB2->localDirection();
	      
	      GlobalVector gdirMB1 = geomDet->toGlobal(ldirMB1);
	      GlobalVector gdirMB2 = geomDet2->toGlobal(ldirMB2);
	      
	      Double_t lPhi2 = atan(ldirMB2.z()/ldirMB2.x());

	      //if ( 1.4 < abs(gdirMB2.phi()) && abs(gdirMB2.phi()) < 1.7) { // Only perpendicular tracks to the detector 

	      //if ( 1.4 < abs(lPhi2) && abs(lPhi2) < 1.7) { // Only perpendicular tracks to local chamber

	      // RPhi angle     
	      double diffRPhi = gdirMB1.phi() - gdirMB2.phi();
	      //double diffRPhi = (atan(ldirMB1.z()/ldirMB1.x())) - (atan(ldirMB2.z()/ldirMB2.x()));
	 
	      stringstream hTitleAngle;
	      hTitleAngle << "hRPhi_W" << iwheel << "_S" << iisect << "_MB" << ist << (ist-1);
	      histo(hTitleAngle.str())->Fill( diffRPhi );

	      // RPhi residuals
	      float distRPhi = (midSegPos-posAtMid).x();
	      
	      stringstream hTitleRPhi;
	      hTitleRPhi << "hresRPhi_W" << iwheel << "_S" << iisect << "_MB" << ist << (ist-1);
	      histo(hTitleRPhi.str())->Fill( distRPhi );

	      // RZ residuals
	      //if(ist1 == 1 || ist1 == 2 || ist1 == 3) {
	      float distZ = (midSegPos-posAtMid).y();
      
	      stringstream hTitlez;
	      hTitlez << "hresZ_W" << iwheel << "_S" << iisect << "_MB" << ist << (ist-1);
	      histo(hTitlez.str())->Fill( distZ );

      
	      //X residuaks vs Phi angle
	      stringstream hTitleXvsPhi;
	      hTitleXvsPhi << "hresXvsPhi_W" << iwheel << "_S" << iisect << "_MB" << ist << (ist-1);
	      histo2d(hTitleXvsPhi.str())->Fill(gdirAtMid.phi(),distRPhi);

	      //X residuaks vs Eta angle
	      stringstream hTitleZvsEta;
	      hTitleZvsEta << "hresZvsEta_W" << iwheel << "_S" << iisect << "_MB" << ist <<  (ist-1);
	      histo2d(hTitleZvsEta.str())->Fill(gdirAtMid.eta(), distZ);
	  
	      //} // end of cut on Phi angle
	    }
	  }
	}// endif segments in MBn-1 & MBn+1   
      } // end loop pn wheel
      
      
      if (debug) cout<<"---"<<endl;  

    } // end loop 4D Segments
  } // endl loop on chambers
}


// =================================================================
void MyDTAnalyzer::analyzeSA(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); 
 
  // Get the magnetic field
  ESHandle<MagneticField> theMGField;
  eventSetup.get<IdealMagneticFieldRecord>().get(theMGField);

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

  if (isMC) {
    Handle<edm::SimTrackContainer> simTracks;
    event.getByLabel(theSimTrackLabel,simTracks);
    cout << "simTracks: " <<simTracks->size() <<endl;
  }

  reco::TrackCollection::const_iterator staTrack;

  double recPt=0.;

  if (debug && staTracks->size() ) 
    cout << endl<<"Run/Event " << event.id().run() << ":" << event.id().event() << " SA tracks " << staTracks->size() << endl;  
  
  // -- Create the propagator
  thePropagator1 = new SteppingHelixPropagator(&*theMGField, alongMomentum);

  // ********* loop on SA tracks ************************************************
  
  for (staTrack = staTracks->begin(); staTrack != staTracks->end(); ++staTrack){
    
    reco::TransientTrack track(*staTrack,&*theMGField,theTrackingGeometry);
    
    TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();

    LocalVector ldirSA = innerTSOS.localDirection();
    
    Double_t lphiSA = atan(ldirSA.z()/ldirSA.x());

    //if (1.4 < abs(lphiSA) && abs(lphiSA) < 1.7) { //only tracks perpendiculaar to each sector

    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());
    histo("hNSegs4D")->Fill(segs->size());

    if (track.normalizedChi2() < theChi2Cut ){
      
      // -- Some Information at IP
      if (debug) {
	cout << " Track state at closest IP :" << endl;
	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;
      }


      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 );
      //Double_t probChi2 = 100*(TMath::Prob(track.chi2(), track.ndof()));
    
      histo("hChi2SA")->Fill(track.normalizedChi2());
      //histo("hProbChi2SA")->Fill(probChi2);
      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);

    // -- Some Information at innerTSOS
        
    //TrajectoryStateOnSurface innerTSOS = track.innermostMeasurementState();
	
	if (debug) {
	  cout << "InnerMost TSOS  (Track State On Surface):"<<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("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());
      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 rechit
      for(trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit){
	
	const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());

	//From Local Hit position to Global Hit position
	GlobalPoint gpos=geomDet->toGlobal((*recHit)->localPosition());
	
	//    cout << (*recHit)->geographicalId() << endl;
	
	//      const TrackingRecHit* trh = (*recHit)->recHits();
	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 (debug)  cout<< wireId << " r: "<< gpos.perp() <<" x,y,z: "<< gpos.x() << " " << gpos.y() << " " << gpos.z() << " " << "phi" << 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();
      } // next associated rechit


      //=====> Residuals of SA Track and the hits

      // -------------  Try to extrapolate using thePropagator segment to segment 
 
      if (debug) 
	cout << "Begin the propagator" << endl; 

      //If the state is valid
      if(innerTSOS.isValid()) {

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

	  // ********** loop on 4D segments *********************************
	  for (DTRecSegment4DCollection::const_iterator seg=segsch.first; seg!=segsch.second; ++seg ) {

	    const GeomDet* geomDet = theTrackingGeometry->idToDet((*seg).geographicalId());
	    int id = (*seg).geographicalId().subdetId();
	    DTChamberId mCham((*seg).geographicalId().rawId());	
	    int station = mCham.station();
	    int sector = mCham.sector();
	    int wheel = mCham.wheel();

	    // -- This try block is to catch some exceptions given by the propagator
	    TrajectoryStateOnSurface destiny;
	    try {
	      
	      destiny = thePropagator1->propagate(*(innerTSOS.freeState()),geomDet->surface());

	      if(!destiny.isValid() || !destiny.hasError()){ 
		continue;
	      }
	    }catch(...){ 
	      
	      cout << "<MyDTAnalyzer> Exception in propagator catched" << endl;
	      continue;
	    }
  
	    LocalPoint lsegPos = (*seg).localPosition();
	    LocalVector lsegDir = (*seg).localDirection();
	    LocalPoint lpredPos = (destiny).localPosition();
	    LocalVector lpredDir = (destiny).localDirection();
	    GlobalVector gsegDir = geomDet->toGlobal(lsegDir);
	    GlobalVector gpredDir = geomDet->toGlobal(lpredDir);

	    Double_t diffRPhi = lsegPos.x() - lpredPos.x();

	    // X Seg - Propagated SA  vs Phi angle
	    stringstream hTitleresXSegSAvsPhi;
	    hTitleresXSegSAvsPhi << "hresXSegSAvsPhi_W" << wheel << "_S" << sector << "_MB" << station;
	    histo2d(hTitleresXSegSAvsPhi.str())->Fill(gpredDir.phi(),diffRPhi);

	  } // End loop on 4D segments 
	} // End loop on chambers
      } 
    } //Cut on Chi2
    //} // tracks perpendicular to each sector
  } // next SA track ============================================
}


//=== Utilities ===//

//=================================================================
double MyDTAnalyzer::extrapolate(char *option_, const DTRecSegment4D seg1_, const DTRecSegment4D seg2_, ESHandle<GlobalTrackingGeometry> &theTrackingGeometry) {
// =================================================================

const GeomDet* geomDet1_ = theTrackingGeometry->idToDet((seg1_).geographicalId());
const GeomDet* geomDet2_ = theTrackingGeometry->idToDet((seg2_).geographicalId());

 double residual = 99999;

//Get GlobalPoition of seg1 in geomDet1 
  GlobalPoint gpos1 = geomDet1_->toGlobal((seg1_).localPosition());
  GlobalVector gdir1 = geomDet1_->toGlobal((seg1_).localDirection());
  
//get all in the frame of geomDet2
  LocalPoint pos1 = (geomDet2_)->toLocal(gpos1);
  LocalVector dir1 = (geomDet2_)->toLocal(gdir1);
  
  LocalPoint pos2 = (seg2_).localPosition();
  LocalVector dir2 = (seg2_).localDirection();

// Extrapolate
  LocalVector unitDir = dir1.unit();
  double lambda = (pos2.z()-pos1.z())/(unitDir.z()); //  >0

  LocalPoint pos1to2 = pos1 + (lambda*unitDir);


  //Calculate RPhi residuals
  if(option_ == "rphi"){
    if(seg1_.hasPhi() && seg2_.hasPhi()){
      residual = (pos2 - pos1to2).x();
    }
  }

 //Calculate Z residuals
  if(option_ == "zed"){
    if(seg1_.hasZed() && seg2_.hasZed()){
      residual = (pos2 - pos1to2).y();
    }
  }
  
//Calculate RPhi angle residuals
  
  if(option_ == "phiangle"){
    if(seg1_.hasPhi() && seg2_.hasPhi()){
      residual = (dir2 - dir1).phi();
    }
  }



  return residual;

}
  


// =================================================================
void MyDTAnalyzer::calculateResiduals(long rawId_, int subDet_M, int station_M,	const GeomDet* geomDet_, DTRecSegment4D seg_, TrajectoryStateOnSurface destiny_){
// =================================================================
  
  if (debug){
	cout << "- Unsing calculateResiduals -" << endl;}


  //The rawId of the detector
  theRawId = rawId_;
  theSubDet = subDet_M;
  theStation = station_M;
  
  //Position and orientation vectors
  theGlobalPosition = geomDet_->toGlobal((seg_).localPosition());
  theGlobalDirection = geomDet_->toGlobal((seg_).localDirection());
  theLocalPosition = (seg_).localPosition();
  theLocalPositionError = (seg_).localPositionError();
  theLocalDirection = (seg_).localDirection();
  theLocalDirectionError = (seg_).localDirectionError();

  //Position and orientation predicted vectors
  thePredictedGlobalPosition = (destiny_).freeState()->position();
  thePredictedGlobalDirection = (destiny_).globalDirection();
  thePredictedLocalPosition = (destiny_).localPosition();
  thePredictedLocalDirection = (destiny_).localDirection();

  //Calculation of residuals
  if(theStation != 4 && theSubDet == 1) {
    //DT internal survey reference system
    theResidualRPhi = theLocalPosition.x() - thePredictedLocalPosition.x();
    theResidualZ = theLocalPosition.y()- thePredictedLocalPosition.y();
    theResidual = sqrt(theLocalPosition.x()*theLocalPosition.x()+theLocalPosition.y()*theLocalPosition.y()) - sqrt(thePredictedLocalPosition.x()*thePredictedLocalPosition.x()+thePredictedLocalPosition.y()*thePredictedLocalPosition.y());
    theResidualPhi = TMath::ATan2( theLocalDirection.x(), theLocalDirection.z() )- TMath::ATan2( thePredictedLocalDirection.x(), thePredictedLocalDirection.z() );
    theResidualTheta = TMath::ATan2( theLocalDirection.z(), theLocalDirection.y() )- TMath::ATan2( thePredictedLocalDirection.z(), thePredictedLocalDirection.y() );
  } else {
    theResidualRPhi = theLocalPosition.x() - thePredictedLocalPosition.x();
    theResidualZ = 0;
    theResidual = theResidualRPhi;
    theResidualPhi = TMath::ATan2( theLocalDirection.x(), theLocalDirection.z() )- TMath::ATan2( thePredictedLocalDirection.x(), thePredictedLocalDirection.z() );
    theResidualTheta = 0;
  }

  theResiduals.ResizeTo(NDOFCoor,1);
  theResiduals(0,0) = theResidualRPhi;
  theResiduals(1,0) = theResidualZ;
  /*theResiduals(2,0) = theResidualPhi;
    theResiduals(3,0) = theResidualTheta;*/
  
   if(TMath::Abs(theResiduals(0,0)) < 5.0 && TMath::Abs(theResiduals(1,0)) < 8.0) {
    validPoint = true;
  } else {
    validPoint = false;
  } 
}



// ================================================================= 
DTRecSegment4D MyDTAnalyzer::matchSegs(Propagator *thePropagator_,TrajectoryStateOnSurface innerTSOS_, edm::ESHandle<GlobalTrackingGeometry> &theTrackingGeometry,  std::vector<DTRecSegment4D> segsch_ ){
// =================================================================

  if(debug) {
      cout << "- Unsing matchSegs -" << endl;}

  
     DTRecSegment4D segSATrack_; 
    
     //double resRPhi = 9999;
     //double resZ = 9999;
     double res = 9999;

    if(innerTSOS_.isValid()) {

  if (segsch_.size() != 0){ 
      
    for (int seg_=0; seg_ < segsch_.size(); seg_++) {

      const GeomDet* geomDet_ = theTrackingGeometry->idToDet((segsch_.at(seg_)).geographicalId());
	  
	 unsigned int id_ = (segsch_.at(seg_)).geographicalId().subdetId();

	 DTChamberId mCham_((segsch_.at(seg_)).geographicalId().rawId());	

// -- Try to associated hits from 4DSegment and SATrack
 
  if (debug) cout<<" -- SA associated RecHits wrt 4D Segment: -- "<<endl;

// -- This try block is to catch some exceptions given by the propagator
     
		     try {
	      
	TrajectoryStateOnSurface destiny_ = thePropagator_->propagate(*(innerTSOS_.freeState()),geomDet_->surface());

	 if(!destiny_.isValid() || !destiny_.hasError()){ 
	   continue;
	 }

	      // For each segment calculate residuals 

      calculateResiduals(mCham_, id_, mCham_.station(), geomDet_, (segsch_.at(seg_)), destiny_);

   
      //if(validPoint && theResidualRPhi < resRPhi ){
	   if(validPoint && theResidual < res ){
	       
	     //resRPhi=theResidualRPhi;
	     res = theResidual;
	     segSATrack_ = segsch_.at(seg_);
	     
	       }
			 
		     }catch(...){ 
	  
     cout << "<MyDTAnalyzer> Exception in propagator catched" << endl;
     continue;
		     }
		     
       } //end loop on segments
     }
    }
	
    return segSATrack_;					     
}		 

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

LocalPoint MyDTAnalyzer::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;
}

// extrapolation between SL  ------------------

LocalPoint MyDTAnalyzer::extrapolate(const DTSLRecSegment2D& segSL1, 
                                   const DTSuperLayerId& SL2Id) const {


  // Get GlobalPosition & direction of 2DSeg in SL1

  GlobalPoint gpos1=(dtGeom->superLayer(segSL1.superLayerId()))->toGlobal(segSL1.localPosition());
  GlobalVector gdir1=(dtGeom->superLayer(segSL1.superLayerId()))->toGlobal(segSL1.localDirection());

     // extrapolate....:
  // get all in SL2 frame
  LocalPoint     pos1=(dtGeom->superLayer(SL2Id))->toLocal(gpos1);
  LocalVector    dir1=(dtGeom->superLayer(SL2Id))->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 -----------from S. Lacaprara--------------------------
// as usual max number of hits and min chi2
const DTRecSegment4D& MyDTAnalyzer::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;
}

 
// best segment utilities -----------from S. Lacaprara--------------------------
// as usual max number of hits and min chi2
const DTRecSegment2D& MyDTAnalyzer::getBestSegment(const DTRecSegment2DCollection::range& segsSL) const{

  DTRecSegment2DCollection::const_iterator bestIter;
  unsigned int nHitBest=0;
  double chi2Best=99999.;
  for (DTRecSegment2DCollection::const_iterator segSL=segsSL.first ;
       segSL!=segsSL.second ; ++segSL ) {
    unsigned int nHits= ((*segSL).recHits()).size();
    
    // we are going to take only segments with 4 hits.

    if (nHits == nHitBest) {

      if ((*segSL).chi2()/(*segSL).degreesOfFreedom() < chi2Best ) {
        chi2Best=(*segSL).chi2()/(*segSL).degreesOfFreedom();
        bestIter = segSL;
      }
    }

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

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


string MyDTAnalyzer::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 MyDTAnalyzer::toString(const DTSuperLayerId& id) const {
  stringstream result;
  result << "_Wh" << id.wheel() 
    << "_Sec" << id.sector() 
    << "_St" << id.station()
    << "_Sl" << id.superLayer();
  return result.str();
}

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

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


void MyDTAnalyzer::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 MyDTAnalyzer::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);
}

DEFINE_FWK_MODULE(MyDTAnalyzer);