/data/refman/pasoursint/CMSSW_6_1_1/src/Alignment/OfflineValidation/plugins/ValidationMisalignedTracker.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    ValidationMisalignedTracker
// Class:      ValidationMisalignedTracker
// 
//
// Original Author:  Nicola De Filippis
//         Created:  Thu Dec 14 13:13:32 CET 2006
// $Id: ValidationMisalignedTracker.cc,v 1.7 2011/12/20 15:11:41 mussgill Exp $
//
//


#include "Alignment/OfflineValidation/plugins/ValidationMisalignedTracker.h"


// user include files

#include "DataFormats/TrackReco/interface/Track.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"

//
// constructors and destructor
//
ValidationMisalignedTracker::ValidationMisalignedTracker(const edm::ParameterSet& iConfig)
{

  //now do what ever initialization is needed
  mzmu=0.,recmzmu=0.,ptzmu=0.,recptzmu=0.,etazmu=0.,recetazmu=0., thetazmu=0.,recthetazmu=0.,phizmu=0.,recphizmu=0.;
  recenezmu=0., enezmu=0., pLzmu=0., recpLzmu=0.,yzmu=0.,recyzmu=0.,mxptmu=0.,recmxptmu=0., minptmu=0.,recminptmu=0.;
  // mzele=0.,recmzele=0.
  
  flag=0,flagrec=0,count=0,countrec=0;
  nAssoc=0;

  for (int i=0;i<2;i++){
    countpart[i]=0;
    countpartrec[i]=0; 
    for (int j=0;j<2;j++){
      ene[i][j]=0.;
      p[i][j]=0.;
      px[i][j]=0.;
      py[i][j]=0.; 
      pz[i][j]=0.;
      ptmu[i][j]=0.; 
      recene[i][j]=0.;
      recp[i][j]=0.;
      recpx[i][j]=0.;
      recpy[i][j]=0.; 
      recpz[i][j]=0.;       
      recptmu[i][j]=0.;
    }
  }


  eventCount_ = 0; 

  selection_eff    = iConfig.getUntrackedParameter<bool>("selection_eff", false);
  selection_fake   = iConfig.getUntrackedParameter<bool>("selection_fake", true);
  ZmassSelection_  = iConfig.getUntrackedParameter<bool>("ZmassSelection", false);
  simobject        = iConfig.getUntrackedParameter<std::string>("simobject","g4SimHits");
  trackassociator  = iConfig.getUntrackedParameter<std::string>("TrackAssociator","ByHits");
  associators      = iConfig.getParameter< std::vector<std::string> >("associators");
  label            = iConfig.getParameter< std::vector<edm::InputTag> >("label");
  label_tp_effic   = iConfig.getParameter< edm::InputTag >("label_tp_effic");
  label_tp_fake    = iConfig.getParameter< edm::InputTag >("label_tp_fake");

  rootfile_   = iConfig.getUntrackedParameter<std::string>("rootfile","myroot.root");
  file_ = new TFile(rootfile_.c_str(),"RECREATE");
  
  // initialize the tree
  tree_eff = new TTree("EffTracks","Efficiency Tracks Tree");

  tree_eff->Branch("Run",&irun,"irun/i");
  tree_eff->Branch("Event",&ievt,"ievt/i");

  // SimTrack
  tree_eff->Branch("TrackID",&trackType,"trackType/i");
  tree_eff->Branch("pt",&pt,"pt/F");
  tree_eff->Branch("eta",&eta,"eta/F");
  tree_eff->Branch("CotTheta",&cottheta,"cottheta/F");
  tree_eff->Branch("phi",&phi,"phi/F");
  tree_eff->Branch("d0",&d0,"d0/F");
  tree_eff->Branch("z0",&z0,"z0/F"); 
  tree_eff->Branch("nhit",&nhit,"nhit/i");
    
  // RecTrack
  tree_eff->Branch("recpt",&recpt,"recpt/F");
  tree_eff->Branch("receta",&receta,"receta/F");
  tree_eff->Branch("CotRecTheta",&reccottheta,"reccottheta/F");
  tree_eff->Branch("recphi",&recphi,"recphi/F");
  tree_eff->Branch("recd0",& recd0,"recd0/F");
  tree_eff->Branch("recz0",& recz0,"recz0/F");
  tree_eff->Branch("nAssoc",&nAssoc,"nAssoc/i");
  tree_eff->Branch("recnhit",&recnhit,"recnhit/i");
  tree_eff->Branch("CHISQ",&recchiq,"recchiq/F");

  tree_eff->Branch("reseta",&reseta,"reseta/F");
  tree_eff->Branch("respt",&respt,"respt/F");
  tree_eff->Branch("resd0",&resd0,"resd0/F");
  tree_eff->Branch("resz0",&resz0,"resz0/F");
  tree_eff->Branch("resphi",&resphi,"resphi/F");
  tree_eff->Branch("rescottheta",&rescottheta,"rescottheta/F");
  tree_eff->Branch("eff",&eff,"eff/F");

  // Invariant masses, pt of Z
  tree_eff->Branch("mzmu",&mzmu,"mzmu/F");
  tree_eff->Branch("ptzmu",&ptzmu,"ptzmu/F");
  tree_eff->Branch("pLzmu",&pLzmu,"pLzmu/F");
  tree_eff->Branch("enezmu",&enezmu,"enezmu/F");
  tree_eff->Branch("etazmu",&etazmu,"etazmu/F");
  tree_eff->Branch("thetazmu",&thetazmu,"thetazmu/F");
  tree_eff->Branch("phizmu",&phizmu,"phizmu/F");
  tree_eff->Branch("yzmu",&yzmu,"yzmu/F");
  tree_eff->Branch("mxptmu",&mxptmu,"mxptmu/F");
  tree_eff->Branch("minptmu",&minptmu,"minptmu/F");

  tree_eff->Branch("recmzmu",&recmzmu,"recmzmu/F");
  tree_eff->Branch("recptzmu",&recptzmu,"recptzmu/F");
  tree_eff->Branch("recpLzmu",&recpLzmu,"recpLzmu/F");  
  tree_eff->Branch("recenezmu",&recenezmu,"recenezmu/F");
  tree_eff->Branch("recetazmu",&recetazmu,"recetazmu/F");
  tree_eff->Branch("recthetazmu",&recthetazmu,"recthetazmu/F");
  tree_eff->Branch("recphizmu",&recphizmu,"recphizmu/F");
  tree_eff->Branch("recyzmu",&recyzmu,"recyzmu/F");
  tree_eff->Branch("recmxptmu",&recmxptmu,"recmxptmu/F");   
  tree_eff->Branch("recminptmu",&recminptmu,"recminptmu/F");       


  //tree->Branch("mzele",&ntmzele,"ntmzele/F");
  //tree->Branch("recmzele",&ntmzeleRec,"ntmzeleRec/F");
  tree_eff->Branch("chi2Associator",&recchiq,"recchiq/F");

  // Fake

  tree_fake = new TTree("FakeTracks","Fake Rate Tracks Tree");

  tree_fake->Branch("Run",&irun,"irun/i");
  tree_fake->Branch("Event",&ievt,"ievt/i");

  // SimTrack
  tree_fake->Branch("fakeTrackID",&faketrackType,"faketrackType/i");
  tree_fake->Branch("fakept",&fakept,"fakept/F");
  tree_fake->Branch("fakeeta",&fakeeta,"fakeeta/F");
  tree_fake->Branch("fakeCotTheta",&fakecottheta,"fakecottheta/F");
  tree_fake->Branch("fakephi",&fakephi,"fakephi/F");
  tree_fake->Branch("faked0",&faked0,"faked0/F");
  tree_fake->Branch("fakez0",&fakez0,"fakez0/F"); 
  tree_fake->Branch("fakenhit",&fakenhit,"fakenhit/i");
    
  // RecTrack
  tree_fake->Branch("fakerecpt",&fakerecpt,"fakerecpt/F");
  tree_fake->Branch("fakereceta",&fakereceta,"fakereceta/F");
  tree_fake->Branch("fakeCotRecTheta",&fakereccottheta,"fakereccottheta/F");
  tree_fake->Branch("fakerecphi",&fakerecphi,"fakerecphi/F");
  tree_fake->Branch("fakerecd0",& fakerecd0,"fakerecd0/F");
  tree_fake->Branch("fakerecz0",& fakerecz0,"fakerecz0/F");
  tree_fake->Branch("fakenAssoc",&fakenAssoc,"fakenAssoc/i");
  tree_fake->Branch("fakerecnhit",&fakerecnhit,"fakerecnhit/i");
  tree_fake->Branch("fakeCHISQ",&fakerecchiq,"fakerecchiq/F");

  tree_fake->Branch("fakereseta",&fakereseta,"fakereseta/F");
  tree_fake->Branch("fakerespt",&fakerespt,"fakerespt/F");
  tree_fake->Branch("fakeresd0",&fakeresd0,"fakeresd0/F");
  tree_fake->Branch("fakeresz0",&fakeresz0,"fakeresz0/F");
  tree_fake->Branch("fakeresphi",&fakeresphi,"fakeresphi/F");
  tree_fake->Branch("fakerescottheta",&fakerescottheta,"fakerescottheta/F");
  tree_fake->Branch("fake",&fake,"fake/F");

  // Invariant masses, pt of Z
  tree_fake->Branch("fakemzmu",&fakemzmu,"fakemzmu/F");
  tree_fake->Branch("fakeptzmu",&fakeptzmu,"fakeptzmu/F");
  tree_fake->Branch("fakepLzmu",&fakepLzmu,"fakepLzmu/F");
  tree_fake->Branch("fakeenezmu",&fakeenezmu,"fakeenezmu/F");
  tree_fake->Branch("fakeetazmu",&fakeetazmu,"fakeetazmu/F");
  tree_fake->Branch("fakethetazmu",&fakethetazmu,"fakethetazmu/F");
  tree_fake->Branch("fakephizmu",&fakephizmu,"fakephizmu/F");
  tree_fake->Branch("fakeyzmu",&fakeyzmu,"fakeyzmu/F");
  tree_fake->Branch("fakemxptmu",&fakemxptmu,"fakemxptmu/F");
  tree_fake->Branch("fakeminptmu",&fakeminptmu,"fakeminptmu/F");

  tree_fake->Branch("fakerecmzmu",&fakerecmzmu,"fakerecmzmu/F");
  tree_fake->Branch("fakerecptzmu",&fakerecptzmu,"fakerecptzmu/F");
  tree_fake->Branch("fakerecpLzmu",&fakerecpLzmu,"fakerecpLzmu/F");  
  tree_fake->Branch("fakerecenezmu",&fakerecenezmu,"fakerecenezmu/F");
  tree_fake->Branch("fakerecetazmu",&fakerecetazmu,"fakerecetazmu/F");
  tree_fake->Branch("fakerecthetazmu",&fakerecthetazmu,"fakerecthetazmu/F");
  tree_fake->Branch("fakerecphizmu",&fakerecphizmu,"fakerecphizmu/F");
  tree_fake->Branch("fakerecyzmu",&fakerecyzmu,"fakerecyzmu/F");
  tree_fake->Branch("fakerecmxptmu",&fakerecmxptmu,"fakerecmxptmu/F");   
  tree_fake->Branch("fakerecminptmu",&fakerecminptmu,"fakerecminptmu/F");       

  tree_fake->Branch("fakechi2Associator",&fakerecchiq,"fakerecchiq/F");

}


ValidationMisalignedTracker::~ValidationMisalignedTracker()
{
 
 
  std::cout << "ValidationMisalignedTracker::endJob Processed " << eventCount_
            << " events" << std::endl;
                                                                                                                   
  // store the tree in the output file
  file_->Write();
                                                                                                                   
  
  // Closing the file deletes the tree.
  file_->Close();
  tree_eff=0;
  tree_fake=0;
}


//
// member functions
//

// ------------ method called to for each event  ------------
void
ValidationMisalignedTracker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   if (watchTrackAssociatorRecord_.check(iSetup)) {
     associatore.clear();
     edm::ESHandle<TrackAssociatorBase> theAssociator;
     for (unsigned int w=0;w<associators.size();w++) {
       iSetup.get<TrackAssociatorRecord>().get(associators[w], theAssociator);
       associatore.push_back( theAssociator.product() );
     }
   }
   
   edm::LogInfo("Tracker Misalignment Validation") << "\n Starting!";

   // Monte Carlo Z selection
   skip=false;
   std::vector<int> indmu;
   
   if ( selection_eff && ZmassSelection_){ 
     edm::Handle<edm::HepMCProduct> evt;
     iEvent.getByLabel("source", evt);
     bool accepted = false;
     bool foundmuons=false;
     HepMC::GenEvent * myGenEvent = new  HepMC::GenEvent(*(evt->GetEvent()));
     
     for ( HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p ) { 
       if ( !accepted && ( (*p)->pdg_id() == 23 ) && (*p)->status() == 3 ) { 
         accepted=true;
         for( HepMC::GenVertex::particle_iterator aDaughter=(*p)->end_vertex()->particles_begin(HepMC::descendants); aDaughter !=(*p)->end_vertex()->particles_end(HepMC::descendants);aDaughter++){
           if ( abs((*aDaughter)->pdg_id())==13) {
             foundmuons=true;
             if ((*aDaughter)->status()!=1 ) {
               for( HepMC::GenVertex::particle_iterator byaDaughter=(*aDaughter)->end_vertex()->particles_begin(HepMC::descendants); byaDaughter !=(*aDaughter)->end_vertex()->particles_end(HepMC::descendants);byaDaughter++){
                 if ((*byaDaughter)->status()==1 && abs((*byaDaughter)->pdg_id())==13) {                   
                   indmu.push_back((*byaDaughter)->barcode());
                   std::cout << "Stable muon from Z with charge " <<  (*byaDaughter)->pdg_id()  << " and index "<< (*byaDaughter)->barcode() << std::endl;
                 }
               }
             }
             else {            
               indmu.push_back((*aDaughter)->barcode());
               std::cout << "Stable muon from Z with charge " <<  (*aDaughter)->pdg_id() << " and index "<< (*aDaughter)->barcode() << std::endl; 
             }       
           }               
         }
         if (!foundmuons){
           std::cout << "No muons from Z ...skip event" << std::endl;
           skip=true;
         } 
       }
     }     
     if ( !accepted) {
       std::cout << "No Z particles in the event ...skip event" << std::endl;
       skip=true;
     }   
   }
   else {
     skip=false;
   }
   
   //
   // Retrieve tracker geometry from event setup
   //
   edm::ESHandle<TrackerGeometry> trackerGeometry;
   iSetup.get<TrackerDigiGeometryRecord>().get( trackerGeometry );
   GeomDet* testGeomDet = trackerGeometry->detsTOB().front();
   std::cout << testGeomDet->position() << std::endl;
   
   
   //Dump Run and Event
   irun=iEvent.id().run();
   ievt=iEvent.id().event();

   // Reset tree variables
   int countpart[2]={0,0},countpartrec[2]={0,0},flag=0,flagrec=0,count=0,countrec=0;
   //int countsim=0;
   float ene[2][2],px[2][2],py[2][2],pz[2][2],ptmu[2][2];
   float recene[2][2],recp[2][2],recpx[2][2],recpy[2][2],recpz[2][2],recptmu[2][2];
   
   for (int i=0;i<2;i++){
     for (int j=0;j<2;j++){
       ene[i][j]=0.;
       px[i][j]=0.;
       py[i][j]=0.; 
       pz[i][j]=0.; 
       ptmu[i][j]=0.;
       recene[i][j]=0.;
       recp[i][j]=0.;
       recpx[i][j]=0.;
       recpy[i][j]=0.; 
       recpz[i][j]=0.;       
       recptmu[i][j]=0.;
     }
   }
   
   
   edm::Handle<TrackingParticleCollection>  TPCollectionHeff ;
   iEvent.getByLabel(label_tp_effic,TPCollectionHeff);
   const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
   
   edm::Handle<TrackingParticleCollection>  TPCollectionHfake ;
   iEvent.getByLabel(label_tp_fake,TPCollectionHfake);
   const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
   
               
   int w=0;
   for (unsigned int ww=0;ww<associators.size();ww++){
      //
      //get collections from the event
      //

      edm::InputTag algo = label[0];
     
      edm::Handle<edm::View<reco::Track> > trackCollection;
      iEvent.getByLabel(algo, trackCollection);
      const edm::View<reco::Track> tC = *(trackCollection.product());
            
           
      //associate tracks
      LogTrace("TrackValidator") << "Calling associateRecoToSim method" << "\n";
           reco::RecoToSimCollection recSimColl=associatore[ww]->associateRecoToSim(trackCollection,
                                                              TPCollectionHfake,
                                                                      &iEvent);
      
      LogTrace("TrackValidator") << "Calling associateSimToReco method" << "\n";
      reco::SimToRecoCollection simRecColl=associatore[ww]->associateSimToReco(trackCollection,
                                                                      TPCollectionHeff, 
                                                                      &iEvent);

   

      //
      //compute number of tracks per eta interval
      //

      if (selection_eff && !skip ) {
        std::cout << "Computing Efficiency" << std::endl;

        edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles (before cuts): " << tPCeff.size() << "\n";
        int ats = 0;
        int st=0;
        for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
          
          // Initialize variables
          eta = 0.,theta=0.,phi=0.,pt=0.,cottheta=0.,costheta=0.;
          d0=0.,z0=0.;
          nhit=0;
          receta = 0.,rectheta = 0.,recphi = 0.,recpt = 0.,reccottheta=0.,recd0=0.,recz0=0.;
          respt = 0.,resd0 = 0.,resz0 = 0.,reseta = 0.,resphi=0.,rescottheta=0.;
          recchiq = 0.;
          recnhit = 0;
          trackType = 0;
          eff=0;        
          
          // typedef edm::Ref<TrackingParticleCollection> TrackingParticleRef;
          TrackingParticleRef tp(TPCollectionHeff, i);
          if (tp->charge()==0) continue;
          st++;
          //pt=sqrt(tp->momentum().perp2());
          //eta=tp->momentum().eta();
          //vpos=tp->vertex().perp2()));
          
          const SimTrack * simulatedTrack = &(*tp->g4Track_begin());
        
          edm::ESHandle<MagneticField> theMF;
          iSetup.get<IdealMagneticFieldRecord>().get(theMF);
          FreeTrajectoryState 
            ftsAtProduction(GlobalPoint(tp->vertex().x(),tp->vertex().y(),tp->vertex().z()),
                            GlobalVector(simulatedTrack->momentum().x(),simulatedTrack->momentum().y(),simulatedTrack->momentum().z()),
                            TrackCharge(tp->charge()),
                            theMF.product());
          TSCPBuilderNoMaterial tscpBuilder;
          TrajectoryStateClosestToPoint tsAtClosestApproach 
            = tscpBuilder(ftsAtProduction,GlobalPoint(0,0,0));//as in TrackProducerAlgorithm
          GlobalPoint v = tsAtClosestApproach.theState().position();
          GlobalVector p = tsAtClosestApproach.theState().momentum();

          //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
          //  double lambdaSim = M_PI/2-p.theta();
          //  double phiSim    = p.phi();
          double dxySim    = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
          double dszSim    = v.z()*p.perp()/p.mag() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.mag();
          d0     = float(-dxySim);
          z0     = float(dszSim*p.mag()/p.perp());

  
          if (abs(simulatedTrack->type())==13 && simulatedTrack->genpartIndex() != -1 ) {
            std::cout << " TRACCIA SIM DI MUONI " << std::endl;  
            std::cout << "Gen part " << simulatedTrack->genpartIndex()<< std::endl;
            trackType=simulatedTrack->type();
            theta=simulatedTrack->momentum().theta();
            costheta=cos(theta);
            cottheta=1./tan(theta);
            
            eta=simulatedTrack->momentum().eta();
            phi=simulatedTrack->momentum().phi();
            pt=simulatedTrack->momentum().pt();
            nhit=tp->matchedHit();
            

            std::cout << "3) Before assoc: SimTrack of type = " << simulatedTrack->type() 
                      << " ,at eta = " << eta 
                      << " ,with pt at vertex = " << simulatedTrack->momentum().pt() << " GeV/c"
                      << " ,d0 =" << d0 
                      << " ,z0 =" << z0 
                      << " ,nhit=" << nhit
                      << std::endl;

            if ( ZmassSelection_ ){
              if (abs(trackType)==13 && (simulatedTrack->genpartIndex()==indmu[0] || simulatedTrack->genpartIndex()==indmu[1] )) { 
                std::cout << " TRACK sim of muons from Z " << std::endl;  
                flag=0;
                count=countpart[0];
                countpart[0]++;
              }
              else if (abs(trackType)==11) {
                //std::cout << " TRACCIA SIM DI ELETTRONI " << std::endl;  
                flag=1;
                count=countpart[1];
                countpart[1]++;
              }
              
              
              px[flag][count]=simulatedTrack->momentum().x();   
              py[flag][count]=simulatedTrack->momentum().y();
              pz[flag][count]=simulatedTrack->momentum().z();
              ptmu[flag][count]=simulatedTrack->momentum().pt(); 
              ene[flag][count]=simulatedTrack->momentum().e();
            }
            
            
            std::vector<std::pair<edm::RefToBase<reco::Track>, double> > rt;
            if(simRecColl.find(tp) != simRecColl.end()){
              
              rt = simRecColl[tp];
              if (rt.size()!=0) {
                
                edm::RefToBase<reco::Track> t = rt.begin()->first;
                ats++;

                // bool flagptused=false;
                // for (unsigned int j=0;j<ptused.size();j++){
                //   if (fabs(t->pt()-ptused[j])<0.001) {
                //     flagptused=true;
                //   }
                // }

                edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st << " with pt=" << t->pt() 
                                                   << " associated with quality:" << rt.begin()->second <<"\n";
                std::cout << "Reconstructed Track:" << t->pt()<< std::endl;
                std::cout << "\tpT: " << t->pt()<< std::endl;
                std::cout << "\timpact parameter:d0: " << t->d0()<< std::endl;
                std::cout << "\timpact parameter:z0: " << t->dz()<< std::endl;
                std::cout << "\tAzimuthal angle of point of closest approach:" << t->phi()<< std::endl;
                std::cout << "\tcharge: " << t->charge()<< std::endl;
                std::cout << "\teta: " << t->eta()<< std::endl;
                std::cout << "\tnormalizedChi2: " << t->normalizedChi2()<< std::endl;

                recnhit=t->numberOfValidHits();
                recchiq=t->normalizedChi2();
                rectheta=t->theta();
                reccottheta=1./tan(rectheta);
                //receta=-log(tan(rectheta/2.));
                receta=t->momentum().eta();
                //         reccostheta=cos(matchedrectrack->momentum().theta());
                recphi=t->phi(); 
                recpt=t->pt();
                ptused.push_back(recpt);
                recd0=t->d0();
                recz0=t->dz();

                std::cout << "5) After call to associator: the best match has " 
                          << recnhit << " hits, Chi2 = " 
                          << recchiq << ", pt at vertex = " 
                          << recpt << " GeV/c, " 
                          << ", recd0 = " << recd0 
                          << ", recz0= " << recz0
                          << std::endl;


                respt=recpt - pt;
                resd0=recd0-d0;
                resz0=recz0-z0;
                reseta=receta-eta;
                resphi=recphi-phi;
                rescottheta=reccottheta-cottheta;
                eff=1;

                std::cout << "6) Transverse momentum residual=" << respt 
                          << " ,d0 residual=" << resd0 
                          << " ,z0 residual=" << resz0 
                          << " with eff=" << eff << std::endl;
                
                if ( ZmassSelection_ ){
                  
                  if (abs(trackType)==13) { 
                    std::cout << " TRACCIA RECO DI MUONI " << std::endl;  
                    flagrec=0;
                    countrec=countpartrec[0];
                    countpartrec[0]++;
                  }
                  else if (abs(trackType)==11) {
                    std::cout << " TRACCIA RECO DI ELETTRONI " << std::endl;  
                    flagrec=1;
                    countrec=countpartrec[1];
                    countpartrec[1]++;
                  }
                  
                  recp[flagrec][countrec]=sqrt(t->momentum().mag2());
                  recpx[flagrec][countrec]=t->momentum().x();   
                  recpy[flagrec][countrec]=t->momentum().y();
                  recpz[flagrec][countrec]=t->momentum().z();
                  recptmu[flagrec][countrec]=sqrt( (t->momentum().x()*t->momentum().x()) + (t->momentum().y()*t->momentum().y()) );
                  if (abs(trackType)==13) recene[flagrec][countrec]=sqrt(recp[flagrec][countrec]*recp[flagrec][countrec]+0.105*0.105);
                  if (abs(trackType)==11) recene[flagrec][countrec]=sqrt(recp[flagrec][countrec]*recp[flagrec][countrec]+0.0005*0.0005);
                }
                
                std::cout << "7) Transverse momentum reconstructed =" << recpt 
                          << " at  eta= " << receta 
                          << " and phi= " << recphi 
                          << std::endl;
                
              }
            }
            else{
              edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st
                                                 << " with pt=" << sqrt(tp->momentum().perp2())
                                                 << " NOT associated to any reco::Track" << "\n";
              receta =-100.;
              recphi =-100.;
              recpt  =-100.;
              recd0  =-100.;
              recz0  =-100;
              respt  =-100.;
              resd0  =-100.;
              resz0  =-100.;
              resphi =-100.;
              reseta =-100.;
              rescottheta=-100.;
              recnhit=100;
              recchiq=-100;
              eff=0;
              flagrec=100;            
            }
            
            std::cout << "Eff=" << eff << std::endl;
               
               // simulated muons
               
               std::cout <<"Flag is" << flag << std::endl;
               std::cout <<"RecFlag is" << flagrec << std::endl;
               
               if (countpart[0]==2 && flag==0) {
                 mzmu=sqrt(
                           (ene[0][0]+ene[0][1])*(ene[0][0]+ene[0][1])-
                           (px[0][0]+px[0][1])*(px[0][0]+px[0][1])-
                           (py[0][0]+py[0][1])*(py[0][0]+py[0][1])-
                           (pz[0][0]+pz[0][1])*(pz[0][0]+pz[0][1])
                           );
                 std::cout << "Mzmu " << mzmu << std::endl;
                 ptzmu=sqrt(
                           (px[0][0]+px[0][1])*(px[0][0]+px[0][1])+
                           (py[0][0]+py[0][1])*(py[0][0]+py[0][1])
                           );
                 
                 pLzmu=pz[0][0]+pz[0][1];
                 enezmu=ene[0][0]+ene[0][1];
                 phizmu=atan2((py[0][0]+py[0][1]),(px[0][0]+px[0][1]));
                 thetazmu=atan2(ptzmu,(pz[0][0]+pz[0][1]));
                 etazmu=-log(tan(thetazmu*3.14/360.));
                 yzmu=0.5*log((enezmu+pLzmu)/(enezmu-pLzmu));
                 mxptmu=std::max( ptmu[0][0], ptmu[0][1]);
                 minptmu=std::min( ptmu[0][0], ptmu[0][1]);
               }
               else {
                 mzmu=-100.;
                 ptzmu=-100.;
                 pLzmu=-100.;
                 enezmu=-100.;
                 etazmu=-100.;
                 phizmu=-100.;
                 thetazmu=-100.;
                 yzmu=-100.;
                 mxptmu=-100.;
                 minptmu=-100.;
               }      

               // reconstructed muons
               if (countpartrec[0]==2 && flagrec==0 ){
                 recmzmu=sqrt(
                              (recene[0][0]+recene[0][1])*(recene[0][0]+recene[0][1])-
                              (recpx[0][0]+recpx[0][1])*(recpx[0][0]+recpx[0][1])-
                              (recpy[0][0]+recpy[0][1])*(recpy[0][0]+recpy[0][1])-
                              (recpz[0][0]+recpz[0][1])*(recpz[0][0]+recpz[0][1])
                              );
                 std::cout << "RecMzmu " << recmzmu << std::endl;
                 recptzmu=sqrt(
                               (recpx[0][0]+recpx[0][1])*(recpx[0][0]+recpx[0][1])+
                               (recpy[0][0]+recpy[0][1])*(recpy[0][0]+recpy[0][1])
                               );
                 
                 recpLzmu=recpz[0][0]+recpz[0][1];      
                 recenezmu=recene[0][0]+recene[0][1];
                 recphizmu=atan2((recpy[0][0]+recpy[0][1]),(recpx[0][0]+recpx[0][1]));
                 recthetazmu=atan2(recptzmu,(recpz[0][0]+recpz[0][1]));
                 recetazmu=-log(tan(recthetazmu*3.14/360.));
                 recyzmu=0.5*log((recenezmu+recpLzmu)/(recenezmu-recpLzmu));
                 recmxptmu=std::max(recptmu[0][0], recptmu[0][1]);
                 recminptmu=std::min( recptmu[0][0], recptmu[0][1]);
               }
               else {
                 recmzmu=-100.;
                 recptzmu=-100.; 
                 recpLzmu=-100.;                 
                 recenezmu=-100.;       
                 recetazmu=-100.;
                 recphizmu=-100.;
                 recthetazmu=-100.;
                 recyzmu=-100.;
                 recmxptmu=-100;
                 recminptmu=-100.;
               }               
               
               tree_eff->Fill();

          } // end of loop on muons
        } // end of loop for tracking particle
      } // end of loop for efficiency

      //
      // Fake Rate
      // 
      if (selection_fake ) {
        std::cout << "Computing Fake Rate" << std::endl;

        fakeeta = 0.,faketheta=0.,fakephi=0.,fakept=0.,fakecottheta=0.,fakecostheta=0.;
        faked0=0.,fakez0=0.;
        fakenhit=0;
        fakereceta = 0.,fakerectheta = 0.,fakerecphi = 0.,fakerecpt = 0.,fakereccottheta=0.,fakerecd0=0.,fakerecz0=0.;
        fakerespt = 0.,fakeresd0 = 0.,fakeresz0 = 0.,fakereseta = 0.,fakeresphi=0.,fakerescottheta=0.;
        fakerecchiq = 0.;
        fakerecnhit = 0;
        faketrackType = 0;
        fake=0;

        
        //      int at=0;
        int rT=0;
        for(reco::TrackCollection::size_type i=0; i<tC.size(); ++i){
          edm::RefToBase<reco::Track> track(trackCollection, i);
          rT++;

          fakeeta = 0.,faketheta=0.,fakephi=0.,fakept=0.,fakecottheta=0.,fakecostheta=0.;
          faked0=0.,fakez0=0.;
          fakenhit=0;
          fakereceta = 0.,fakerectheta = 0.,fakerecphi = 0.,fakerecpt = 0.,fakereccottheta=0.,fakerecd0=0.,fakerecz0=0.;
          fakerespt = 0.,fakeresd0 = 0.,fakeresz0 = 0.,fakereseta = 0.,fakeresphi=0.,fakerescottheta=0.;
          fakerecchiq = 0.;
          fakerecnhit = 0;
          faketrackType = 0;
          fake=0;
          
          fakerecnhit=track->numberOfValidHits();
          fakerecchiq=track->normalizedChi2();
          fakerectheta=track->theta();
          fakereccottheta=1./tan(rectheta);
          //fakereceta=-log(tan(rectheta/2.));
          fakereceta=track->momentum().eta();
          //       fakereccostheta=cos(track->momentum().theta());
          fakerecphi=track->phi(); 
          fakerecpt=track->pt();
          fakerecd0=track->d0();
          fakerecz0=track->dz();
                  
          std::cout << "1) Before assoc: TkRecTrack at eta = " << fakereceta << std::endl;
          std::cout << "Track number "<< i << std::endl ;
          std::cout << "\tPT: " << track->pt()<< std::endl;
          std::cout << "\timpact parameter:d0: " << track->d0()<< std::endl;
          std::cout << "\timpact parameter:z0: " << track->dz()<< std::endl;
          std::cout << "\tAzimuthal angle of point of closest approach:" << track->phi()<< std::endl;
          std::cout << "\tcharge: " << track->charge()<< std::endl;
          std::cout << "\teta: " << track->eta()<< std::endl;
          std::cout << "\tnormalizedChi2: " << track->normalizedChi2()<< std::endl;
           
           
          std::vector<std::pair<TrackingParticleRef, double> > tp;
          
          //Compute fake rate vs eta
          if(recSimColl.find(track) != recSimColl.end()){
            tp = recSimColl[track];
            if (tp.size()!=0) {
              edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt() 
                                                 << " associated with quality:" << tp.begin()->second <<"\n";
              
              
              TrackingParticleRef tpr = tp.begin()->first;
              const SimTrack * fakeassocTrack = &(*tpr->g4Track_begin());

              edm::ESHandle<MagneticField> theMF;
              iSetup.get<IdealMagneticFieldRecord>().get(theMF);
              FreeTrajectoryState 
                ftsAtProduction(GlobalPoint(tpr->vertex().x(),tpr->vertex().y(),tpr->vertex().z()),
                                GlobalVector(fakeassocTrack->momentum().x(),fakeassocTrack->momentum().y(),fakeassocTrack->momentum().z()),
                                TrackCharge(tpr->charge()),
                                theMF.product());
              TSCPBuilderNoMaterial tscpBuilder;
              TrajectoryStateClosestToPoint tsAtClosestApproach 
                = tscpBuilder(ftsAtProduction,GlobalPoint(0,0,0));//as in TrackProducerAlgorithm
              GlobalPoint v = tsAtClosestApproach.theState().position();
              GlobalVector p = tsAtClosestApproach.theState().momentum();
              
              //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
              //  double lambdaSim = M_PI/2-p.theta();
              //  double phiSim    = p.phi();
              double dxySim    = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
              double dszSim    = v.z()*p.perp()/p.mag() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.mag();
              faked0     = float(-dxySim);
              fakez0     = float(dszSim*p.mag()/p.perp());


              faketrackType=fakeassocTrack->type();
              faketheta=fakeassocTrack->momentum().theta();
              fakecottheta=1./tan(faketheta);
              fakeeta=fakeassocTrack->momentum().eta();
              fakephi=fakeassocTrack->momentum().phi();
              fakept=fakeassocTrack->momentum().pt();
              fakenhit=tpr->matchedHit();

              std::cout << "4) After call to associator: the best SimTrack match is of type" << fakeassocTrack->type() 
                        << " ,at eta = " << fakeeta 
                        << " and phi = " << fakephi  
                        << " ,with pt at vertex = " << fakept << " GeV/c" 
                        << " ,d0 global = " << faked0 
                        << " ,z0 = " << fakez0
                        << std::endl;
              fake=1;

              fakerespt=fakerecpt-fakept;
              fakeresd0=fakerecd0-faked0;
              fakeresz0=fakerecz0-fakez0;
              fakereseta=-log(tan(fakerectheta/2.))-(-log(tan(faketheta/2.)));
              fakeresphi=fakerecphi-fakephi;
              fakerescottheta=fakereccottheta-fakecottheta;
              
            }
          }
          else{
            edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
                                               << " NOT associated to any TrackingParticle" << "\n";
            
            fakeeta =-100.;
            faketheta=-100;
            fakephi =-100.;
            fakept  =-100.;
            faked0  =-100.;
            fakez0  =-100;
            fakerespt  =-100.;
            fakeresd0  =-100.;
            fakeresz0  =-100.;
            fakeresphi =-100.;
            fakereseta =-100.;
            fakerescottheta=-100.;
            fakenhit=100;
            fake=0;
          }
          
          tree_fake->Fill();
        }       
        
      } // End of loop on fakerate       
      
      w++;
      
   } // End of loop on associators
}

// ------------ method called once each job just after ending the event loop  ------------
void ValidationMisalignedTracker::endJob() {

  std::cout << "\t Misalignment analysis completed \n" << std::endl;  

}

DEFINE_FWK_MODULE(ValidationMisalignedTracker);