CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/Alignment/OfflineValidation/plugins/PrimaryVertexValidation.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    PrimaryVertexValidation
00004 // Class:      PrimaryVertexValidation
00005 // 
00013 //
00014 // Original Author:  Marco Musich
00015 //         Created:  Tue Mar 02 10:39:34 CDT 2010
00016 //
00017 
00018 // system include files
00019 #include <memory>
00020 
00021 // user include files
00022 #include "Alignment/OfflineValidation/plugins/PrimaryVertexValidation.h"
00023 #include "FWCore/Framework/interface/Frameworkfwd.h"
00024 #include "FWCore/Framework/interface/EDAnalyzer.h"
00025 
00026 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
00027 #include "FWCore/Framework/interface/Event.h"
00028 #include <FWCore/Framework/interface/ESHandle.h>
00029 #include "FWCore/Framework/interface/MakerMacros.h"
00030 
00031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00032 
00033 #include <SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h>
00034 #include <SimDataFormats/TrackingHit/interface/PSimHit.h>
00035 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00036 
00037 #include "TH1F.h"
00038 #include "TH2F.h"
00039 #include "TF1.h"
00040 #include "TVector3.h"
00041 #include "TFile.h"
00042 #include "TROOT.h"
00043 #include "TChain.h"
00044 #include "TNtuple.h"
00045 #include <TMatrixD.h>
00046 #include <TVectorD.h>
00047 #include "FWCore/ServiceRegistry/interface/Service.h"
00048 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00049 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00050 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00051 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00052 #include <Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h>
00053 #include <DataFormats/GeometrySurface/interface/Surface.h>
00054 #include <DataFormats/GeometrySurface/interface/GloballyPositioned.h>
00055 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00056 #include "MagneticField/Engine/interface/MagneticField.h" 
00057 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" 
00058 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
00059 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00060 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00061 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00062 
00063 #include "DataFormats/TrackReco/interface/Track.h"
00064 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00065 #include "DataFormats/VertexReco/interface/Vertex.h"
00066 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00067 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h"
00068 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h"
00069 
00070  const int kBPIX = PixelSubdetector::PixelBarrel;
00071  const int kFPIX = PixelSubdetector::PixelEndcap;
00072 
00073 // Constructor
00074 
00075 PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfig):
00076   storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
00077   lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
00078   useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
00079   askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
00080   ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt",0.)),
00081   etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta",2.4)),
00082   nBins_(iConfig.getUntrackedParameter<int>("numberOfBins",24)),
00083   debug_(iConfig.getParameter<bool>("Debug")),
00084   TrackCollectionTag_(iConfig.getParameter<edm::InputTag>("TrackCollectionTag"))
00085 {
00086   
00087   // now do what ever initialization is needed
00088   // initialize phase space boundaries
00089   
00090   phipitch_ = (2*TMath::Pi())/nBins_;
00091   etapitch_ = 5./nBins_;
00092 
00093   // old version
00094   // theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
00095 
00096   // select and configure the track filter 
00097   theTrackFilter_= new TrackFilterForPVFinding(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters") );
00098   // select and configure the track clusterizer  
00099   std::string clusteringAlgorithm=iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
00100   if (clusteringAlgorithm=="gap"){
00101     theTrackClusterizer_ = new GapClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkGapClusParameters"));
00102   }else if(clusteringAlgorithm=="DA"){
00103     theTrackClusterizer_ = new DAClusterizerInZ(iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<edm::ParameterSet>("TkDAClusParameters"));
00104   }else{
00105     throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);  
00106   }
00107 }
00108    
00109 // Destructor
00110 PrimaryVertexValidation::~PrimaryVertexValidation()
00111 {
00112    // do anything here that needs to be done at desctruction time
00113    // (e.g. close files, deallocate resources etc.)
00114 }
00115 
00116 
00117 //
00118 // member functions
00119 //
00120 
00121 // ------------ method called to for each event  ------------
00122 void
00123 PrimaryVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00124 {
00125 
00126   using namespace std;
00127   using namespace IPTools;
00128 
00129   Nevt_++;  
00130 
00131   //=======================================================
00132   // Initialize Root-tuple variables
00133   //=======================================================
00134 
00135   SetVarToZero();
00136  
00137   //=======================================================
00138   // Retrieve the Magnetic Field information
00139   //=======================================================
00140 
00141   edm::ESHandle<MagneticField> theMGField;
00142   iSetup.get<IdealMagneticFieldRecord>().get( theMGField );
00143 
00144   //=======================================================
00145   // Retrieve the Tracking Geometry information
00146   //=======================================================
00147 
00148   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
00149   iSetup.get<GlobalTrackingGeometryRecord>().get( theTrackingGeometry );
00150 
00151   //=======================================================
00152   // Retrieve the Track information
00153   //=======================================================
00154   
00155   edm::Handle<reco::TrackCollection>  trackCollectionHandle;
00156   iEvent.getByLabel(TrackCollectionTag_, trackCollectionHandle);
00157   
00158   //=======================================================
00159   // Retrieve offline vartex information (only for reco)
00160   //=======================================================
00161  
00162   edm::Handle<reco::VertexCollection> vertices;
00163   try {
00164     iEvent.getByLabel("offlinePrimaryVertices", vertices);
00165   } catch (...) {
00166     if(debug_)
00167       cout << "No offlinePrimaryVertices found!" << endl;
00168   }
00169   if ( vertices.isValid() ) {
00170     xOfflineVertex_ = (*vertices)[0].x();
00171     yOfflineVertex_ = (*vertices)[0].y();
00172     zOfflineVertex_ = (*vertices)[0].z();
00173   }
00174 
00175   unsigned int vertexCollectionSize = vertices.product()->size();
00176   int nvvertex = 0;
00177   
00178   for (unsigned int i=0; i<vertexCollectionSize; i++) {
00179     const reco::Vertex& vertex = vertices->at(i);
00180     if (vertex.isValid()) nvvertex++;
00181   }
00182 
00183   nOfflineVertices_ = nvvertex;
00184  
00185 
00186   if ( vertices->size() && useTracksFromRecoVtx_ ) {
00187    
00188     double sumpt    = 0;
00189     size_t ntracks  = 0;
00190     double chi2ndf  = 0.; 
00191     double chi2prob = 0.;
00192 
00193     if (!vertices->at(0).isFake()) {
00194       
00195       reco::Vertex pv = vertices->at(0);
00196       
00197       ntracks  = pv.tracksSize();
00198       chi2ndf  = pv.normalizedChi2();
00199       chi2prob = TMath::Prob(pv.chi2(),(int)pv.ndof());
00200       
00201       h_recoVtxNtracks_->Fill(ntracks);   
00202       h_recoVtxChi2ndf_->Fill(chi2ndf);    
00203       h_recoVtxChi2Prob_->Fill(chi2prob);   
00204       
00205       for (reco::Vertex::trackRef_iterator itrk = pv.tracks_begin();itrk != pv.tracks_end(); ++itrk) {
00206         double pt = (**itrk).pt();
00207         sumpt += pt*pt;
00208         
00209         const math::XYZPoint myVertex(pv.position().x(),pv.position().y(),pv.position().z());
00210         
00211         double dxyRes = (**itrk).dxy(myVertex);
00212         double dzRes  = (**itrk).dz(myVertex);
00213         
00214         double dxy_err = (**itrk).dxyError();
00215         double dz_err  = (**itrk).dzError();
00216         
00217         float_t trackphi = ((**itrk).phi())*(180/TMath::Pi());
00218         float_t tracketa = (**itrk).eta();
00219         
00220         for(Int_t i=0; i<nBins_; i++){
00221           
00222           float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
00223           float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
00224           
00225           float etaF=-2.5+i*etapitch_;
00226           float etaL=-2.5+(i+1)*etapitch_;
00227           
00228           if(trackphi >= phiF && trackphi < phiL ){
00229             
00230             a_dxyPhiBiasResiduals[i]->Fill(dxyRes*cmToum);
00231             a_dzPhiBiasResiduals[i]->Fill(dzRes*cmToum); 
00232             n_dxyPhiBiasResiduals[i]->Fill((dxyRes)/dxy_err);
00233             n_dzPhiBiasResiduals[i]->Fill((dzRes)/dz_err); 
00234             
00235             for(Int_t j=0; j<nBins_; j++){
00236               
00237               float etaJ=-2.5+j*etapitch_;
00238               float etaK=-2.5+(j+1)*etapitch_;
00239               
00240               if(tracketa >= etaJ && tracketa < etaK ){
00241                 
00242                 a_dxyBiasResidualsMap[i][j]->Fill(dxyRes*cmToum); 
00243                 a_dzBiasResidualsMap[i][j]->Fill(dzRes*cmToum);   
00244                 
00245                 n_dxyBiasResidualsMap[i][j]->Fill((dxyRes)/dxy_err); 
00246                 n_dzBiasResidualsMap[i][j]->Fill((dzRes)/dz_err);  
00247                 
00248               }
00249             }
00250           }             
00251           
00252           if(tracketa >= etaF && tracketa < etaL ){
00253             a_dxyEtaBiasResiduals[i]->Fill(dxyRes*cmToum);
00254             a_dzEtaBiasResiduals[i]->Fill(dzRes*cmToum); 
00255             n_dxyEtaBiasResiduals[i]->Fill((dxyRes)/dxy_err);
00256             n_dzEtaBiasResiduals[i]->Fill((dzRes)/dz_err);
00257           }
00258         }
00259       }
00260       
00261       h_recoVtxSumPt_->Fill(sumpt);   
00262 
00263     }
00264   }
00265 
00266   //=======================================================
00267   // Retrieve Beamspot information
00268   //=======================================================
00269 
00270   reco::BeamSpot beamSpot;
00271   edm::Handle<reco::BeamSpot> beamSpotHandle;
00272   iEvent.getByLabel("offlineBeamSpot", beamSpotHandle);
00273     
00274   if ( beamSpotHandle.isValid() )
00275     {
00276       beamSpot = *beamSpotHandle;
00277       BSx0_ = beamSpot.x0();
00278       BSy0_ = beamSpot.y0();
00279       BSz0_ = beamSpot.z0();
00280       Beamsigmaz_ = beamSpot.sigmaZ();    
00281       Beamdxdz_ = beamSpot.dxdz();           
00282       BeamWidthX_ = beamSpot.BeamWidthX();
00283       BeamWidthY_ = beamSpot.BeamWidthY();
00284     } else
00285     {
00286       if(debug_)
00287         cout << "No BeamSpot found!" << endl;
00288     }
00289   
00290   if(debug_)
00291     std::cout<<"Beamspot x:"<<BSx0_<<" y:"<<BSy0_<<" z:"<<BSz0_<<std::endl; 
00292   
00293   //double sigmaz = beamSpot.sigmaZ();
00294   //double dxdz = beamSpot.dxdz();
00295   //double BeamWidth = beamSpot.BeamWidth();
00296   
00297   //=======================================================
00298   // Starts here ananlysis
00299   //=======================================================
00300   
00301   RunNumber_=iEvent.eventAuxiliary().run();
00302   LuminosityBlockNumber_=iEvent.eventAuxiliary().luminosityBlock();
00303   EventNumber_=iEvent.eventAuxiliary().id().event();
00304     
00305   if(debug_)
00306     std::cout<<"PrimaryVertexValidation::analyze() looping over "<<trackCollectionHandle->size()<< "tracks." <<std::endl;       
00307 
00308   //======================================================
00309   // Interface RECO tracks to vertex reconstruction
00310   //======================================================
00311  
00312   std::vector<reco::TransientTrack> t_tks;
00313   unsigned int k = 0;   
00314   for(reco::TrackCollection::const_iterator track = trackCollectionHandle->begin(); track!= trackCollectionHandle->end(); ++track, ++k){
00315   
00316     reco::TrackRef trackref(trackCollectionHandle,k);
00317     reco::TransientTrack theTTRef = reco::TransientTrack(trackref, &*theMGField, theTrackingGeometry );
00318     t_tks.push_back(theTTRef);
00319   
00320   }
00321   
00322   if(debug_) {cout << "PrimaryVertexValidation"
00323                   << "Found: " << t_tks.size() << " reconstructed tracks" << "\n";
00324   }
00325   
00326   //======================================================
00327   // select the tracks
00328   //======================================================
00329 
00330   std::vector<reco::TransientTrack> seltks = theTrackFilter_->select(t_tks);
00331     
00332   //======================================================
00333   // clusterize tracks in Z
00334   //======================================================
00335 
00336   vector< vector<reco::TransientTrack> > clusters = theTrackClusterizer_->clusterize(seltks);
00337   
00338   if (debug_){
00339     cout <<  " clustering returned  "<< clusters.size() << " clusters  from " << t_tks.size() << " selected tracks" <<endl;
00340   }
00341   
00342   nClus_=clusters.size();  
00343 
00344   //======================================================
00345   // Starts loop on clusters 
00346   //======================================================
00347 
00348   for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin(); iclus != clusters.end(); iclus++) {
00349 
00350     nTracksPerClus_=0;
00351 
00352     unsigned int i = 0;   
00353     for(vector<reco::TransientTrack>::const_iterator theTrack = iclus->begin(); theTrack!= iclus->end(); ++theTrack, ++i)
00354       {
00355         if ( nTracks_ >= nMaxtracks_ ) {
00356           std::cout << " PrimaryVertexValidation::analyze() : Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_ << std::endl;
00357           continue;
00358         }
00359         
00360         pt_[nTracks_]       = theTrack->track().pt();
00361         p_[nTracks_]        = theTrack->track().p();
00362         nhits_[nTracks_]    = theTrack->track().numberOfValidHits();
00363         eta_[nTracks_]      = theTrack->track().eta();
00364         theta_[nTracks_]    = theTrack->track().theta();
00365         phi_[nTracks_]      = theTrack->track().phi();
00366         chi2_[nTracks_]     = theTrack->track().chi2();
00367         chi2ndof_[nTracks_] = theTrack->track().normalizedChi2();
00368         charge_[nTracks_]   = theTrack->track().charge();
00369         qoverp_[nTracks_]   = theTrack->track().qoverp();
00370         dz_[nTracks_]       = theTrack->track().dz();
00371         dxy_[nTracks_]      = theTrack->track().dxy();
00372         
00373         reco::TrackBase::TrackQuality _trackQuality = reco::TrackBase::qualityByName("highPurity");
00374         isHighPurity_[nTracks_] = theTrack->track().quality(_trackQuality);
00375         
00376         math::XYZPoint point(beamSpot.x0(),beamSpot.y0(), beamSpot.z0());
00377         dxyBs_[nTracks_]    = theTrack->track().dxy(point);
00378         dzBs_[nTracks_]     = theTrack->track().dz(point);
00379 
00380         xPCA_[nTracks_]     = theTrack->track().vertex().x();
00381         yPCA_[nTracks_]     = theTrack->track().vertex().y();
00382         zPCA_[nTracks_]     = theTrack->track().vertex().z(); 
00383         
00384         //=======================================================
00385         // Retrieve rechit information
00386         //=======================================================  
00387         
00388         int nRecHit1D =0;
00389         int nRecHit2D =0;
00390         int nhitinTIB =0; 
00391         int nhitinTOB =0; 
00392         int nhitinTID =0; 
00393         int nhitinTEC =0; 
00394         int nhitinBPIX=0;
00395         int nhitinFPIX=0; 
00396         
00397         for (trackingRecHit_iterator iHit = theTrack->recHitsBegin(); iHit != theTrack->recHitsEnd(); ++iHit) {
00398           if((*iHit)->isValid()) {      
00399             
00400             if (this->isHit2D(**iHit)) {++nRecHit2D;}
00401             else {++nRecHit1D; }
00402             
00403             int type =(*iHit)->geographicalId().subdetId();
00404             
00405             if(type==int(StripSubdetector::TIB)){++nhitinTIB;}
00406             if(type==int(StripSubdetector::TOB)){++nhitinTOB;}
00407             if(type==int(StripSubdetector::TID)){++nhitinTID;}
00408             if(type==int(StripSubdetector::TEC)){++nhitinTEC;}
00409             if(type==int(                kBPIX)){++nhitinBPIX;}
00410             if(type==int(                kFPIX)){++nhitinFPIX;}
00411           }
00412         }      
00413 
00414         nhits1D_[nTracks_]     = nRecHit1D;
00415         nhits2D_[nTracks_]     = nRecHit2D;
00416         nhitsBPIX_[nTracks_]   = nhitinBPIX;
00417         nhitsFPIX_[nTracks_]   = nhitinFPIX;
00418         nhitsTIB_[nTracks_]    = nhitinTIB;
00419         nhitsTID_[nTracks_]    = nhitinTID;
00420         nhitsTOB_[nTracks_]    = nhitinTOB;
00421         nhitsTEC_[nTracks_]    = nhitinTEC;
00422         
00423         //=======================================================
00424         // Good tracks for vertexing selection
00425         //=======================================================  
00426 
00427         bool pass = true;
00428         if(askFirstLayerHit_) pass = this->hasFirstLayerPixelHits((*theTrack));
00429         if (pass){
00430           isGoodTrack_[nTracks_]=1;
00431         }
00432       
00433         //=======================================================
00434         // Fit unbiased vertex
00435         //=======================================================
00436         
00437         vector<reco::TransientTrack> theFinalTracks;
00438         theFinalTracks.clear();
00439 
00440         for(vector<reco::TransientTrack>::const_iterator tk = iclus->begin(); tk!= iclus->end(); ++tk){
00441           
00442           pass = this->hasFirstLayerPixelHits((*tk));
00443           if (pass){
00444             if( tk == theTrack ) continue;
00445             else {
00446               theFinalTracks.push_back((*tk));
00447             }
00448           }
00449         }
00450         
00451         if(theFinalTracks.size() > 2){
00452             
00453           if(debug_)
00454             std::cout <<"PrimaryVertexValidation::analyze() :Transient Track Collection size: "<<theFinalTracks.size()<<std::endl;
00455           
00456           try{
00457               
00458             VertexFitter<5>* fitter = new AdaptiveVertexFitter;
00459             TransientVertex theFittedVertex = fitter->vertex(theFinalTracks);
00460             
00461             if(theFittedVertex.isValid ()){
00462               
00463               if(theFittedVertex.hasTrackWeight()){
00464                 for(size_t rtracks= 0; rtracks < theFinalTracks.size(); rtracks++){
00465                   sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTracks[rtracks]);
00466                 }
00467               }
00468               
00469               const math::XYZPoint theRecoVertex(xOfflineVertex_,yOfflineVertex_,zOfflineVertex_);
00470               
00471               const math::XYZPoint myVertex(theFittedVertex.position().x(),theFittedVertex.position().y(),theFittedVertex.position().z());
00472               hasRecVertex_[nTracks_]    = 1;
00473               xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
00474               yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
00475               zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
00476               
00477               chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
00478               chi2UnbiasedVertex_[nTracks_]     = theFittedVertex.totalChiSquared();
00479               DOFUnbiasedVertex_[nTracks_]      = theFittedVertex.degreesOfFreedom();   
00480               tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
00481 
00482               h_fitVtxNtracks_->Fill(theFinalTracks.size());        
00483               h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
00484               h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
00485               h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());        
00486               h_fitVtxChi2Prob_->Fill(TMath::Prob(theFittedVertex.totalChiSquared(),(int)theFittedVertex.degreesOfFreedom()));       
00487                        
00488               dszFromMyVertex_[nTracks_] = theTrack->track().dsz(myVertex);
00489               dxyFromMyVertex_[nTracks_] = theTrack->track().dxy(myVertex);
00490               dzFromMyVertex_[nTracks_]  = theTrack->track().dz(myVertex);
00491               
00492               double dz_err = hypot(theTrack->track().dzError(),theFittedVertex.positionError().czz());
00493 
00494               // PV2D 
00495 
00496               // std::pair<bool,Measurement1D> ip2dpv = absoluteTransverseImpactParameter(*theTrack,theFittedVertex);
00497               // double ip2d_corr = ip2dpv.second.value();
00498               // double ip2d_err  = ip2dpv.second.error();
00499 
00500               std::pair<bool,Measurement1D> s_ip2dpv =
00501                 signedTransverseImpactParameter(*theTrack,GlobalVector(theTrack->track().px(),
00502                                                                        theTrack->track().py(),
00503                                                                        theTrack->track().pz()),
00504                                                 theFittedVertex); 
00505               
00506               // double s_ip2dpv_corr = s_ip2dpv.second.value();
00507               double s_ip2dpv_err  = s_ip2dpv.second.error();
00508               
00509               // PV3D
00510 
00511               // std::pair<bool,Measurement1D> ip3dpv = absoluteImpactParameter3D(*theTrack,theFittedVertex);
00512               // double ip3d_corr = ip3dpv.second.value(); 
00513               // double ip3d_err  = ip3dpv.second.error(); 
00514 
00515               // std::pair<bool,Measurement1D> s_ip3dpv = 
00516               //        signedImpactParameter3D(*theTrack,GlobalVector(theTrack->track().px(),
00517               //                                                       theTrack->track().py(),
00518               //                                                       theTrack->track().pz()),
00519               //                        theFittedVertex);
00520               
00521               // double s_ip3dpv_corr = s_ip3dpv.second.value();
00522               // double s_ip3dpv_err  = s_ip3dpv.second.error();
00523               
00524               dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
00525               dzErrorFromMyVertex_[nTracks_]  = dz_err;
00526               
00527               IPTsigFromMyVertex_[nTracks_]   = (theTrack->track().dxy(myVertex))/s_ip2dpv_err;
00528               IPLsigFromMyVertex_[nTracks_]   = (theTrack->track().dz(myVertex))/dz_err;
00529               
00530               // fill directly the histograms of residuals
00531               
00532               float_t trackphi = (theTrack->track().phi())*(180/TMath::Pi());
00533               float_t tracketa = theTrack->track().eta();
00534               float_t trackpt  = theTrack->track().pt();
00535 
00536               // checks on the probe track quality
00537               if(trackpt >= ptOfProbe_ && fabs(tracketa)<= etaOfProbe_){
00538 
00539                 h_probePt_->Fill(theTrack->track().pt());
00540                 h_probeEta_->Fill(theTrack->track().eta());
00541                 h_probePhi_->Fill(theTrack->track().phi());
00542                 h_probeChi2_->Fill(theTrack->track().chi2());
00543                 h_probeNormChi2_->Fill(theTrack->track().normalizedChi2());
00544                 h_probeCharge_->Fill(theTrack->track().charge());
00545                 h_probeQoverP_->Fill(theTrack->track().qoverp());
00546                 h_probedz_->Fill(theTrack->track().dz(theRecoVertex));
00547                 h_probedxy_->Fill((theTrack->track().dxy(theRecoVertex)));
00548 
00549                 h_probeHits_->Fill(theTrack->track().numberOfValidHits());       
00550                 h_probeHits1D_->Fill(nRecHit1D);
00551                 h_probeHits2D_->Fill(nRecHit2D);
00552                 h_probeHitsInTIB_->Fill(nhitinBPIX);
00553                 h_probeHitsInTOB_->Fill(nhitinFPIX);
00554                 h_probeHitsInTID_->Fill(nhitinTIB);
00555                 h_probeHitsInTEC_->Fill(nhitinTID);
00556                 h_probeHitsInBPIX_->Fill(nhitinTOB);
00557                 h_probeHitsInFPIX_->Fill(nhitinTEC);
00558                 
00559                 for(Int_t i=0; i<nBins_; i++){
00560                   
00561                   float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
00562                   float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
00563                   
00564                   float etaF=-2.5+i*etapitch_;
00565                   float etaL=-2.5+(i+1)*etapitch_;
00566                           
00567                   if(trackphi >= phiF && trackphi < phiL ){
00568                     a_dxyPhiResiduals[i]->Fill(theTrack->track().dxy(myVertex)*cmToum);
00569                     a_dzPhiResiduals[i]->Fill(theTrack->track().dz(myVertex)*cmToum); 
00570                     n_dxyPhiResiduals[i]->Fill((theTrack->track().dxy(myVertex))/s_ip2dpv_err);
00571                     n_dzPhiResiduals[i]->Fill((theTrack->track().dz(myVertex))/dz_err); 
00572 
00573                     for(Int_t j=0; j<nBins_; j++){
00574 
00575                       float etaJ=-2.5+j*etapitch_;
00576                       float etaK=-2.5+(j+1)*etapitch_;
00577 
00578                       if(tracketa >= etaJ && tracketa < etaK ){
00579                         
00580                         a_dxyResidualsMap[i][j]->Fill(theTrack->track().dxy(myVertex)*cmToum); 
00581                         a_dzResidualsMap[i][j]->Fill(theTrack->track().dz(myVertex)*cmToum);   
00582                         
00583                         n_dxyResidualsMap[i][j]->Fill((theTrack->track().dxy(myVertex))/s_ip2dpv_err); 
00584                         n_dzResidualsMap[i][j]->Fill((theTrack->track().dz(myVertex))/dz_err);  
00585                         
00586                       }
00587                     }
00588                   }             
00589                   
00590                   if(tracketa >= etaF && tracketa < etaL ){
00591                     a_dxyEtaResiduals[i]->Fill(theTrack->track().dxy(myVertex)*cmToum);
00592                     a_dzEtaResiduals[i]->Fill(theTrack->track().dz(myVertex)*cmToum); 
00593                     n_dxyEtaResiduals[i]->Fill((theTrack->track().dxy(myVertex))/s_ip2dpv_err);
00594                     n_dzEtaResiduals[i]->Fill((theTrack->track().dz(myVertex))/dz_err);
00595                   }
00596                 }
00597               }
00598                   
00599               if(debug_){
00600                 std::cout<<"PrimaryVertexValidation::analyze() : myVertex.x()= "<<myVertex.x()<<" myVertex.y()= "<<myVertex.y()<<" theFittedVertex.z()= "<<myVertex.z()<<std::endl;       
00601                 std::cout<<"PrimaryVertexValidation::analyze() : theTrack->track().dz(myVertex)= "<<theTrack->track().dz(myVertex)<<std::endl;
00602                 std::cout<<"PrimaryVertexValidation::analyze() : zPCA -myVertex.z() = "<<(theTrack->track().vertex().z() -myVertex.z() )<<std::endl; 
00603               }// ends if debug_
00604             } // ends if the fitted vertex is Valid
00605 
00606             delete fitter;
00607 
00608           }  catch ( cms::Exception& er ) {
00609             LogTrace("PrimaryVertexValidation::analyze RECO")<<"caught std::exception "<<er.what()<<std::endl;
00610           }
00611                 
00612         } //ends if theFinalTracks.size() > 2
00613         
00614         else {
00615           if(debug_)
00616               std::cout << "PrimaryVertexValidation::analyze() :Not enough tracks to make a vertex.  Returns no vertex info" << std::endl;
00617         }
00618           
00619         ++nTracks_;  
00620         ++nTracksPerClus_;
00621 
00622         if(debug_)
00623           cout<< "Track "<<i<<" : pT = "<<theTrack->track().pt()<<endl;
00624         
00625       }// for loop on tracks
00626 
00627   } // for loop on track clusters
00628   
00629 
00630   // Fill the TTree if needed
00631 
00632   if(storeNtuple_){  
00633     rootTree_->Fill();
00634   }
00635 
00636   
00637 } 
00638 
00639 // ------------ method called to discriminate 1D from 2D hits  ------------
00640 bool PrimaryVertexValidation::isHit2D(const TrackingRecHit &hit) const
00641 {
00642   if (hit.dimension() < 2) {
00643     return false; // some (muon...) stuff really has RecHit1D
00644   } else {
00645     const DetId detId(hit.geographicalId());
00646     if (detId.det() == DetId::Tracker) {
00647       if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
00648         return true; // pixel is always 2D
00649       } else { // should be SiStrip now
00650         if (dynamic_cast<const SiStripRecHit2D*>(&hit)) return false; // normal hit
00651         else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit)) return true; // matched is 2D
00652         else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) return false; // crazy hit...
00653         else {
00654           edm::LogError("UnkownType") << "@SUB=AlignmentTrackSelector::isHit2D"
00655                                       << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
00656                                       << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
00657           return false;
00658         }
00659       }
00660     } else { // not tracker??
00661       edm::LogWarning("DetectorMismatch") << "@SUB=AlignmentTrackSelector::isHit2D"
00662                                           << "Hit not in tracker with 'official' dimension >=2.";
00663       return true; // dimension() >= 2 so accept that...
00664     }
00665   }
00666   // never reached...
00667 }
00668 
00669 // ------------ method to check the presence of pixel hits  ------------
00670 bool PrimaryVertexValidation::hasFirstLayerPixelHits(const reco::TransientTrack track)
00671 {
00672   bool accepted = false;
00673   // hit pattern of the track
00674   const reco::HitPattern& p = track.hitPattern();      
00675   for (int i=0; i<p.numberOfHits(); i++) {
00676     uint32_t pattern = p.getHitPattern(i);   
00677     if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern) ) {
00678       if (p.getLayer(pattern) == 1) {
00679         if (p.validHitFilter(pattern)) {
00680           accepted = true;
00681         }
00682       }
00683      }
00684   }
00685   return accepted;
00686 } 
00687 
00688 // ------------ method called once each job before begining the event loop  ------------
00689 void PrimaryVertexValidation::beginJob()
00690 {
00691   edm::LogInfo("beginJob") << "Begin Job" << std::endl;
00692   // Define TTree for output
00693   Nevt_    = 0;
00694   
00695   //  rootFile_ = new TFile(filename_.c_str(),"recreate");
00696   edm::Service<TFileService> fs;
00697   rootTree_ = fs->make<TTree>("tree","PV Validation tree");
00698  
00699   // Track Paramters 
00700 
00701   if(lightNtupleSwitch_){ 
00702 
00703     rootTree_->Branch("EventNumber",&EventNumber_,"EventNumber/i");
00704     rootTree_->Branch("RunNumber",&RunNumber_,"RunNumber/i");
00705     rootTree_->Branch("LuminosityBlockNumber",&LuminosityBlockNumber_,"LuminosityBlockNumber/i");
00706     rootTree_->Branch("nOfflineVertices",&nOfflineVertices_,"nOfflineVertices/I");
00707     rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
00708     rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
00709     rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
00710     rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
00711     rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
00712     rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D"); 
00713     rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");    
00714     rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D"); 
00715     rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
00716     rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");
00717     rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity_[nTracks]/I");
00718     
00719   } else {
00720     
00721     rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
00722     rootTree_->Branch("nTracksPerClus",&nTracksPerClus_,"nTracksPerClus/I");
00723     rootTree_->Branch("nClus",&nClus_,"nClus/I");
00724     rootTree_->Branch("xOfflineVertex",&xOfflineVertex_,"xOfflineVertex/D");
00725     rootTree_->Branch("yOfflineVertex",&yOfflineVertex_,"yOfflineVertex/D");
00726     rootTree_->Branch("zOfflineVertex",&zOfflineVertex_,"zOfflineVertex/D");
00727     rootTree_->Branch("BSx0",&BSx0_,"BSx0/D");
00728     rootTree_->Branch("BSy0",&BSy0_,"BSy0/D");
00729     rootTree_->Branch("BSz0",&BSz0_,"BSz0/D");
00730     rootTree_->Branch("Beamsigmaz",&Beamsigmaz_,"Beamsigmaz/D");
00731     rootTree_->Branch("Beamdxdz",&Beamdxdz_,"Beamdxdz/D");
00732     rootTree_->Branch("BeamWidthX",&BeamWidthX_,"BeamWidthX/D");
00733     rootTree_->Branch("BeamWidthY",&BeamWidthY_,"BeamWidthY/D");
00734     rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
00735     rootTree_->Branch("p",&p_,"p[nTracks]/D");
00736     rootTree_->Branch("nhits",&nhits_,"nhits[nTracks]/I");
00737     rootTree_->Branch("nhits1D",&nhits1D_,"nhits1D[nTracks]/I");
00738     rootTree_->Branch("nhits2D",&nhits2D_,"nhits2D[nTracks]/I");
00739     rootTree_->Branch("nhitsBPIX",&nhitsBPIX_,"nhitsBPIX[nTracks]/I");
00740     rootTree_->Branch("nhitsFPIX",&nhitsFPIX_,"nhitsFPIX[nTracks]/I");
00741     rootTree_->Branch("nhitsTIB",&nhitsTIB_,"nhitsTIB[nTracks]/I");
00742     rootTree_->Branch("nhitsTID",&nhitsTID_,"nhitsTID[nTracks]/I");
00743     rootTree_->Branch("nhitsTOB",&nhitsTOB_,"nhitsTOB[nTracks]/I");
00744     rootTree_->Branch("nhitsTEC",&nhitsTEC_,"nhitsTEC[nTracks]/I");
00745     rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
00746     rootTree_->Branch("theta",&theta_,"theta[nTracks]/D");
00747     rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
00748     rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
00749     rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
00750     rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
00751     rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
00752     rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
00753     rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
00754     rootTree_->Branch("dzBs",&dzBs_,"dzBs[nTracks]/D");
00755     rootTree_->Branch("dxyBs",&dxyBs_,"dxyBs[nTracks]/D");
00756     rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
00757     rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
00758     rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
00759     rootTree_->Branch("xUnbiasedVertex",&xUnbiasedVertex_,"xUnbiasedVertex[nTracks]/D");
00760     rootTree_->Branch("yUnbiasedVertex",&yUnbiasedVertex_,"yUnbiasedVertex[nTracks]/D");
00761     rootTree_->Branch("zUnbiasedVertex",&zUnbiasedVertex_,"zUnbiasedVertex[nTracks]/D");
00762     rootTree_->Branch("chi2normUnbiasedVertex",&chi2normUnbiasedVertex_,"chi2normUnbiasedVertex[nTracks]/F");
00763     rootTree_->Branch("chi2UnbiasedVertex",&chi2UnbiasedVertex_,"chi2UnbiasedVertex[nTracks]/F");
00764     rootTree_->Branch("DOFUnbiasedVertex",&DOFUnbiasedVertex_," DOFUnbiasedVertex[nTracks]/F");
00765     rootTree_->Branch("sumOfWeightsUnbiasedVertex",&sumOfWeightsUnbiasedVertex_,"sumOfWeightsUnbiasedVertex[nTracks]/F");
00766     rootTree_->Branch("tracksUsedForVertexing",&tracksUsedForVertexing_,"tracksUsedForVertexing[nTracks]/I");
00767     rootTree_->Branch("dxyFromMyVertex",&dxyFromMyVertex_,"dxyFromMyVertex[nTracks]/D");
00768     rootTree_->Branch("dzFromMyVertex",&dzFromMyVertex_,"dzFromMyVertex[nTracks]/D");
00769     rootTree_->Branch("dszFromMyVertex",&dszFromMyVertex_,"dszFromMyVertex[nTracks]/D");
00770     rootTree_->Branch("dxyErrorFromMyVertex",&dxyErrorFromMyVertex_,"dxyErrorFromMyVertex_[nTracks]/D"); 
00771     rootTree_->Branch("dzErrorFromMyVertex",&dzErrorFromMyVertex_,"dzErrorFromMyVertex_[nTracks]/D");  
00772     rootTree_->Branch("IPTsigFromMyVertex",&IPTsigFromMyVertex_,"IPTsigFromMyVertex_[nTracks]/D");    
00773     rootTree_->Branch("IPLsigFromMyVertex",&IPLsigFromMyVertex_,"IPLsigFromMyVertex_[nTracks]/D"); 
00774     rootTree_->Branch("hasRecVertex",&hasRecVertex_,"hasRecVertex[nTracks]/I");
00775     rootTree_->Branch("isGoodTrack",&isGoodTrack_,"isGoodTrack[nTracks]/I");   
00776   }
00777 
00778   // probe track histograms
00779 
00780   TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
00781 
00782   h_probePt_         = ProbeFeatures.make<TH1F>("h_probePt","p_{T} of probe track;track p_{T} (GeV); tracks",100,0.,50.);   
00783   h_probeEta_        = ProbeFeatures.make<TH1F>("h_probeEta","#eta of probe track;track #eta; tracks",54,-2.7,2.7);  
00784   h_probePhi_        = ProbeFeatures.make<TH1F>("h_probePhi","#phi of probe track;track #phi [rad]; tracks",100,-3.15,3.15);  
00785   h_probeChi2_       = ProbeFeatures.make<TH1F>("h_probeChi2","#chi^{2} of probe track;track #chi^{2}; tracks",100,0.,100.); 
00786   h_probeNormChi2_   = ProbeFeatures.make<TH1F>("h_probeNormChi2"," normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks",100,0.,10.);
00787   h_probeCharge_     = ProbeFeatures.make<TH1F>("h_probeCharge","charge of profe track;track charge Q;tracks",3,-1.5,1.5);
00788   h_probeQoverP_     = ProbeFeatures.make<TH1F>("h_probeQoverP","q/p of probe track; track Q/p (GeV^{-1})",200,-1.,1.);
00789   h_probedz_         = ProbeFeatures.make<TH1F>("h_probedz","d_{z} of probe track;track d_{z} (cm);tracks",100,-5.,5.);  
00790   h_probedxy_        = ProbeFeatures.make<TH1F>("h_probedxy","d_{xy} of probe track;track d_{xy} (#mum);tracks",200,-1.,1.);      
00791 
00792   h_probeHits_       = ProbeFeatures.make<TH1F>("h_probeNRechits"    ,"N_{hits}     ;N_{hits}    ;tracks",40,-0.5,39.5);
00793   h_probeHits1D_     = ProbeFeatures.make<TH1F>("h_probeNRechits1D"  ,"N_{hits} 1D  ;N_{hits} 1D ;tracks",40,-0.5,39.5);
00794   h_probeHits2D_     = ProbeFeatures.make<TH1F>("h_probeNRechits2D"  ,"N_{hits} 2D  ;N_{hits} 2D ;tracks",40,-0.5,39.5);
00795   h_probeHitsInTIB_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTIB" ,"N_{hits} TIB ;N_{hits} TIB;tracks",40,-0.5,39.5);
00796   h_probeHitsInTOB_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTOB" ,"N_{hits} TOB ;N_{hits} TOB;tracks",40,-0.5,39.5);
00797   h_probeHitsInTID_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTID" ,"N_{hits} TID ;N_{hits} TID;tracks",40,-0.5,39.5);
00798   h_probeHitsInTEC_  = ProbeFeatures.make<TH1F>("h_probeNRechitsTEC" ,"N_{hits} TEC ;N_{hits} TEC;tracks",40,-0.5,39.5);
00799   h_probeHitsInBPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX","N_{hits} BPIX;N_{hits} BPIX;tracks",40,-0.5,39.5);
00800   h_probeHitsInFPIX_ = ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX","N_{hits} FPIX;N_{hits} FPIX;tracks",40,-0.5,39.5);
00801 
00802   TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
00803   h_fitVtxNtracks_          = RefitVertexFeatures.make<TH1F>("h_fitVtxNtracks"  ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices"        ,100,-0.5,99.5);
00804   h_fitVtxNdof_             = RefitVertexFeatures.make<TH1F>("h_fitVtxNdof"     ,"N^{vtx}_{DOF};N^{vtx}_{DOF};vertices"          ,100,-0.5,99.5);
00805   h_fitVtxChi2_             = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2"     ,"#chi^{2} vtx;#chi^{2} vtx;vertices"            ,100,-0.5,99.5);
00806   h_fitVtxChi2ndf_          = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2ndf"  ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices"    ,100,-0.5,9.5);
00807   h_fitVtxChi2Prob_         = RefitVertexFeatures.make<TH1F>("h_fitVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
00808 
00809   if(useTracksFromRecoVtx_) {
00810 
00811     TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
00812     h_recoVtxNtracks_          = RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks"  ,"N^{vtx}_{trks};N^{vtx}_{trks};vertices"        ,100,-0.5,99.5);
00813     h_recoVtxChi2ndf_          = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf"  ,"#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices"    ,10,-0.5,9.5);
00814     h_recoVtxChi2Prob_         = RecoVertexFeatures.make<TH1F>("h_recoVtxChi2Prob" ,"Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices",40,0.,1.);
00815     h_recoVtxSumPt_            = RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt"    ,"Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices"  ,100,0.,200.);
00816     
00817   }
00818 
00819   // initialize the residuals histograms 
00820 
00821   float dxymax_phi = 2000; 
00822   float dzmax_phi  = 2000; 
00823   float dxymax_eta = 3000; 
00824   float dzmax_eta  = 3000;
00825 
00826   const Int_t mybins_ = 500;
00827 
00829   //  
00830   // Usual plots from refitting the vertex
00831   // The vertex is refit without the probe track
00832   //
00834 
00835   TFileDirectory AbsTransPhiRes  = fs->mkdir("Abs_Transv_Phi_Residuals");
00836   TFileDirectory AbsTransEtaRes  = fs->mkdir("Abs_Transv_Eta_Residuals");
00837                                                           
00838   TFileDirectory AbsLongPhiRes   = fs->mkdir("Abs_Long_Phi_Residuals");
00839   TFileDirectory AbsLongEtaRes   = fs->mkdir("Abs_Long_Eta_Residuals");
00840                                   
00841   TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
00842   TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
00843                                                           
00844   TFileDirectory NormLongPhiRes  = fs->mkdir("Norm_Long_Phi_Residuals");
00845   TFileDirectory NormLongEtaRes  = fs->mkdir("Norm_Long_Eta_Residuals");
00846 
00847   TFileDirectory AbsDoubleDiffRes   = fs->mkdir("Abs_DoubleDiffResiduals");
00848   TFileDirectory NormDoubleDiffRes  = fs->mkdir("Norm_DoubleDiffResiduals");
00849 
00850   for ( int i=0; i<nBins_; ++i ) {
00851 
00852     float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
00853     float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
00854     
00855     float etaF=-2.5+i*etapitch_;
00856     float etaL=-2.5+(i+1)*etapitch_;
00857     
00858     a_dxyPhiResiduals[i] = AbsTransPhiRes.make<TH1F>(Form("histo_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi,dxymax_phi);
00859     a_dxyEtaResiduals[i] = AbsTransEtaRes.make<TH1F>(Form("histo_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta,dxymax_eta);
00860 
00861     a_dzPhiResiduals[i]  = AbsLongPhiRes.make<TH1F>(Form("histo_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f #circ;d_{z} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi,dzmax_phi);
00862     a_dzEtaResiduals[i]  = AbsLongEtaRes.make<TH1F>(Form("histo_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta,dzmax_eta);
00863                                                         
00864     n_dxyPhiResiduals[i] = NormTransPhiRes.make<TH1F>(Form("histo_norm_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi/100.,dxymax_phi/100.);
00865     n_dxyEtaResiduals[i] = NormTransEtaRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta/100.,dxymax_eta/100.);
00866 
00867     n_dzPhiResiduals[i]  = NormLongPhiRes.make<TH1F>(Form("histo_norm_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi/100.,dzmax_phi/100.);
00868     n_dzEtaResiduals[i]  = NormLongEtaRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta/100.,dzmax_eta/100.);
00869 
00870     for ( int j=0; j<nBins_; ++j ) {
00871       
00872       a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
00873       a_dzResidualsMap[i][j]  = AbsDoubleDiffRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
00874                                        
00875       n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
00876       n_dzResidualsMap[i][j]  = NormDoubleDiffRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
00877     }
00878 
00879   }
00880 
00881   // declaration of the directories
00882 
00883   TFileDirectory MeanTrendsDir   = fs->mkdir("MeanTrends");
00884   TFileDirectory WidthTrendsDir  = fs->mkdir("WidthTrends");
00885   TFileDirectory MedianTrendsDir = fs->mkdir("MedianTrends");
00886   TFileDirectory MADTrendsDir    = fs->mkdir("MADTrends");
00887 
00888   TFileDirectory Mean2DMapsDir   = fs->mkdir("MeanMaps");
00889   TFileDirectory Width2DMapsDir  = fs->mkdir("WidthMaps");
00890 
00891   Double_t highedge=nBins_-0.5;
00892   Double_t lowedge=-0.5;
00893 
00894   // means and widths from the fit
00895 
00896   a_dxyPhiMeanTrend  = MeanTrendsDir.make<TH1F> ("means_dxy_phi","#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge); 
00897   a_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_phi","#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
00898   a_dzPhiMeanTrend   = MeanTrendsDir.make<TH1F> ("means_dz_phi","#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge); 
00899   a_dzPhiWidthTrend  = WidthTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
00900                           
00901   a_dxyEtaMeanTrend  = MeanTrendsDir.make<TH1F> ("means_dxy_eta","#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge);
00902   a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
00903   a_dzEtaMeanTrend   = MeanTrendsDir.make<TH1F> ("means_dz_eta","#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge); 
00904   a_dzEtaWidthTrend  = WidthTrendsDir.make<TH1F>("widths_dz_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
00905                           
00906   n_dxyPhiMeanTrend  = MeanTrendsDir.make<TH1F> ("norm_means_dxy_phi","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",nBins_,lowedge,highedge);
00907   n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_phi","width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
00908   n_dzPhiMeanTrend   = MeanTrendsDir.make<TH1F> ("norm_means_dz_phi","#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge); 
00909   n_dzPhiWidthTrend  = WidthTrendsDir.make<TH1F>("norm_widths_dz_phi","width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
00910                                                                                         
00911   n_dxyEtaMeanTrend  = MeanTrendsDir.make<TH1F> ("norm_means_dxy_eta","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
00912   n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("norm_widths_dxy_eta","width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
00913   n_dzEtaMeanTrend   = MeanTrendsDir.make<TH1F> ("norm_means_dz_eta","#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);  
00914   n_dzEtaWidthTrend  = WidthTrendsDir.make<TH1F>("norm_widths_dz_eta","width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                        
00915   
00916   // 2D maps
00917 
00918   a_dxyMeanMap       =  Mean2DMapsDir.make<TH2F>  ("means_dxy_map","#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00919   a_dzMeanMap        =  Mean2DMapsDir.make<TH2F>  ("means_dz_map","#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00920                      
00921   n_dxyMeanMap       =  Mean2DMapsDir.make<TH2F>  ("norm_means_dxy_map","#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00922   n_dzMeanMap        =  Mean2DMapsDir.make<TH2F>  ("norm_means_dz_map","#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00923                      
00924   a_dxyWidthMap      =  Width2DMapsDir.make<TH2F> ("widths_dxy_map","#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00925   a_dzWidthMap       =  Width2DMapsDir.make<TH2F> ("widths_dz_map","#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00926                      
00927   n_dxyWidthMap      =  Width2DMapsDir.make<TH2F> ("norm_widths_dxy_map","width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00928   n_dzWidthMap       =  Width2DMapsDir.make<TH2F> ("norm_widths_dz_map","width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
00929 
00930   // medians and MADs
00931 
00932   a_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_phi","Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);                         
00933   a_dxyPhiMADTrend    = MADTrendsDir.make<TH1F>   ("MADs_dxy_phi","Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);                              
00934   a_dzPhiMedianTrend  = MedianTrendsDir.make<TH1F>("medians_dz_phi","Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);                                   
00935   a_dzPhiMADTrend     = MADTrendsDir.make<TH1F>   ("MADs_dz_phi","Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);                                 
00936   
00937   a_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("medians_dxy_eta","Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);                                       
00938   a_dxyEtaMADTrend    = MADTrendsDir.make<TH1F>   ("MADs_dxy_eta","Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);                                           
00939   a_dzEtaMedianTrend  = MedianTrendsDir.make<TH1F>("medians_dz_eta","Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);                                          
00940   a_dzEtaMADTrend     = MADTrendsDir.make<TH1F>   ("MADs_dz_eta","Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);                                    
00941   
00942   n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_phi","Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge); 
00943   n_dxyPhiMADTrend    = MADTrendsDir.make<TH1F>   ("norm_MADs_dxy_phi","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);   
00944   n_dzPhiMedianTrend  = MedianTrendsDir.make<TH1F>("norm_medians_dz_phi","Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);      
00945   n_dzPhiMADTrend     = MADTrendsDir.make<TH1F>   ("norm_MADs_dz_phi","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);     
00946   
00947   n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>("norm_medians_dxy_eta","Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                
00948   n_dxyEtaMADTrend    = MADTrendsDir.make<TH1F>   ("norm_MADs_dxy_eta","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);              
00949   n_dzEtaMedianTrend  = MedianTrendsDir.make<TH1F>("norm_medians_dz_eta","Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                    
00950   n_dzEtaMADTrend     = MADTrendsDir.make<TH1F>   ("norm_MADs_dz_eta","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                      
00951 
00952 
00954   //  
00955   // plots of biased residuals
00956   // The vertex is refit without the probe track
00957   //
00959 
00960   if (useTracksFromRecoVtx_){
00961 
00962     TFileDirectory AbsTransPhiBiasRes  = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
00963     TFileDirectory AbsTransEtaBiasRes  = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
00964     
00965     TFileDirectory AbsLongPhiBiasRes   = fs->mkdir("Abs_Long_Phi_BiasResiduals");
00966     TFileDirectory AbsLongEtaBiasRes   = fs->mkdir("Abs_Long_Eta_BiasResiduals");
00967     
00968     TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
00969     TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
00970     
00971     TFileDirectory NormLongPhiBiasRes  = fs->mkdir("Norm_Long_Phi_BiasResiduals");
00972     TFileDirectory NormLongEtaBiasRes  = fs->mkdir("Norm_Long_Eta_BiasResiduals");
00973     
00974     TFileDirectory AbsDoubleDiffBiasRes   = fs->mkdir("Abs_DoubleDiffBiasResiduals");
00975     TFileDirectory NormDoubleDiffBiasRes  = fs->mkdir("Norm_DoubleDiffBiasResiduals");
00976     
00977     for ( int i=0; i<nBins_; ++i ) {
00978       
00979       float phiF = (-TMath::Pi()+i*phipitch_)*(180/TMath::Pi());
00980       float phiL = (-TMath::Pi()+(i+1)*phipitch_)*(180/TMath::Pi());
00981       
00982       float etaF=-2.5+i*etapitch_;
00983       float etaL=-2.5+(i+1)*etapitch_;
00984       
00985       a_dxyPhiBiasResiduals[i] = AbsTransPhiBiasRes.make<TH1F>(Form("histo_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi,dxymax_phi);
00986       a_dxyEtaBiasResiduals[i] = AbsTransEtaBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta,dxymax_eta);
00987       
00988       a_dzPhiBiasResiduals[i]  = AbsLongPhiBiasRes.make<TH1F>(Form("histo_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f #circ;d_{z} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi,dzmax_phi);
00989       a_dzEtaBiasResiduals[i]  = AbsLongEtaBiasRes.make<TH1F>(Form("histo_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta,dzmax_eta);
00990       
00991       n_dxyPhiBiasResiduals[i] = NormTransPhiBiasRes.make<TH1F>(Form("histo_norm_dxy_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",phiF,phiL),mybins_,-dxymax_phi/100.,dxymax_phi/100.);
00992       n_dxyEtaBiasResiduals[i] = NormTransEtaBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL),mybins_,-dxymax_eta/100.,dxymax_eta/100.);
00993       
00994       n_dzPhiBiasResiduals[i]  = NormLongPhiBiasRes.make<TH1F>(Form("histo_norm_dz_phi_plot%i",i),Form("%.2f#circ<#varphi^{probe}_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",phiF,phiL),mybins_,-dzmax_phi/100.,dzmax_phi/100.);
00995       n_dzEtaBiasResiduals[i]  = NormLongEtaBiasRes.make<TH1F>(Form("histo_norm_dz_eta_plot%i",i),Form("%.2f<#eta^{probe}_{tk}<%.2f;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL),mybins_,-dzmax_eta/100.,dzmax_eta/100.);
00996       
00997       for ( int j=0; j<nBins_; ++j ) {
00998         
00999       a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
01000       a_dzBiasResidualsMap[i][j]  = AbsDoubleDiffBiasRes.make<TH1F>(Form("histo_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta,dzmax_eta);
01001       
01002       n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
01003       n_dzBiasResidualsMap[i][j]  = NormDoubleDiffBiasRes.make<TH1F>(Form("histo_norm_dxy_eta_plot%i_phi_plot%i",i,j),Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}} [#mum];tracks",etaF,etaL,phiF,phiL),mybins_,-dzmax_eta/100,dzmax_eta/100);
01004       }
01005       
01006     }
01007     
01008     // declaration of the directories
01009 
01010     TFileDirectory MeanBiasTrendsDir   = fs->mkdir("MeanBiasTrends");
01011     TFileDirectory WidthBiasTrendsDir  = fs->mkdir("WidthBiasTrends");
01012     TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
01013     TFileDirectory MADBiasTrendsDir    = fs->mkdir("MADBiasTrends");
01014     
01015     TFileDirectory Mean2DBiasMapsDir   = fs->mkdir("MeanBiasMaps");
01016     TFileDirectory Width2DBiasMapsDir  = fs->mkdir("WidthBiasMaps");
01017     
01018     // means and widths from the fit
01019     
01020     a_dxyPhiMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("means_dxy_phi","#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge); 
01021     a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_phi","#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
01022     a_dzPhiMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("means_dz_phi","#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge); 
01023     a_dzPhiWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("widths_dz_phi","#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
01024     
01025     a_dxyEtaMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("means_dxy_eta","#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]",nBins_,lowedge,highedge);
01026     a_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",nBins_,lowedge,highedge);
01027     a_dzEtaMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("means_dz_eta","#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]",nBins_,lowedge,highedge); 
01028     a_dzEtaWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("widths_dz_eta","#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",nBins_,lowedge,highedge);
01029     
01030     n_dxyPhiMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_phi","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",nBins_,lowedge,highedge);
01031     n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_phi","width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);
01032     n_dzPhiMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_phi","#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge); 
01033     n_dzPhiWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_phi","width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
01034     
01035     n_dxyEtaMeanBiasTrend  = MeanBiasTrendsDir.make<TH1F> ("norm_means_dxy_eta","#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);
01036     n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>("norm_widths_dxy_eta","width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);
01037     n_dzEtaMeanBiasTrend   = MeanBiasTrendsDir.make<TH1F> ("norm_means_dz_eta","#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",nBins_,lowedge,highedge);  
01038     n_dzEtaWidthBiasTrend  = WidthBiasTrendsDir.make<TH1F>("norm_widths_dz_eta","width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                        
01039     
01040     // 2D maps
01041     
01042     a_dxyMeanBiasMap       =  Mean2DBiasMapsDir.make<TH2F>  ("means_dxy_map","#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01043     a_dzMeanBiasMap        =  Mean2DBiasMapsDir.make<TH2F>  ("means_dz_map","#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01044     
01045     n_dxyMeanBiasMap       =  Mean2DBiasMapsDir.make<TH2F>  ("norm_means_dxy_map","#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01046     n_dzMeanBiasMap        =  Mean2DBiasMapsDir.make<TH2F>  ("norm_means_dz_map","#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01047     
01048     a_dxyWidthBiasMap      =  Width2DBiasMapsDir.make<TH2F> ("widths_dxy_map","#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01049     a_dzWidthBiasMap       =  Width2DBiasMapsDir.make<TH2F> ("widths_dz_map","#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01050     
01051     n_dxyWidthBiasMap      =  Width2DBiasMapsDir.make<TH2F> ("norm_widths_dxy_map","width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01052     n_dzWidthBiasMap       =  Width2DBiasMapsDir.make<TH2F> ("norm_widths_dz_map","width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",nBins_,lowedge,highedge,nBins_,lowedge,highedge);
01053     
01054     // medians and MADs
01055     
01056     a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_phi","Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);                       
01057     a_dxyPhiMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("MADs_dxy_phi","Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);                                    
01058     a_dzPhiMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("medians_dz_phi","Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);                                 
01059     a_dzPhiMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("MADs_dz_phi","Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);                               
01060     
01061     a_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta","Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",nBins_,lowedge,highedge);                                             
01062     a_dxyEtaMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("MADs_dxy_eta","Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",nBins_,lowedge,highedge);                                         
01063     a_dzEtaMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("medians_dz_eta","Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",nBins_,lowedge,highedge);                                        
01064     a_dzEtaMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("MADs_dz_eta","Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",nBins_,lowedge,highedge);                                          
01065     
01066     n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_phi","Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge); 
01067     n_dxyPhiMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("norm_MADs_dxy_phi","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",nBins_,lowedge,highedge);   
01068     n_dzPhiMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_phi","Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);      
01069     n_dzPhiMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("norm_MADs_dz_phi","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);           
01070     
01071     n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>("norm_medians_dxy_eta","Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);              
01072     n_dxyEtaMADBiasTrend    = MADBiasTrendsDir.make<TH1F>   ("norm_MADs_dxy_eta","Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                    
01073     n_dzEtaMedianBiasTrend  = MedianBiasTrendsDir.make<TH1F>("norm_medians_dz_eta","Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                  
01074     n_dzEtaMADBiasTrend     = MADBiasTrendsDir.make<TH1F>   ("norm_MADs_dz_eta","Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",nBins_,lowedge,highedge);                      
01075     
01076   }
01077 }
01078 // ------------ method called once each job just after ending the event loop  ------------
01079 void PrimaryVertexValidation::endJob() 
01080 {
01081 
01082   std::cout<<"######################################"<<std::endl;
01083   std::cout<<"Number of analyzed events: "<<Nevt_<<std::endl;
01084   std::cout<<"######################################"<<std::endl;
01085 
01086   if(useTracksFromRecoVtx_){
01087 
01088     FillTrendPlot(a_dxyPhiMeanBiasTrend ,a_dxyPhiBiasResiduals,"mean","phi");  
01089     FillTrendPlot(a_dxyPhiWidthBiasTrend,a_dxyPhiBiasResiduals,"width","phi");
01090     FillTrendPlot(a_dzPhiMeanBiasTrend  ,a_dzPhiBiasResiduals ,"mean","phi");   
01091     FillTrendPlot(a_dzPhiWidthBiasTrend ,a_dzPhiBiasResiduals ,"width","phi");  
01092     
01093     FillTrendPlot(a_dxyEtaMeanBiasTrend ,a_dxyEtaBiasResiduals,"mean","eta"); 
01094     FillTrendPlot(a_dxyEtaWidthBiasTrend,a_dxyEtaBiasResiduals,"width","eta");
01095     FillTrendPlot(a_dzEtaMeanBiasTrend  ,a_dzEtaBiasResiduals ,"mean","eta"); 
01096     FillTrendPlot(a_dzEtaWidthBiasTrend ,a_dzEtaBiasResiduals ,"width","eta");
01097     
01098     FillTrendPlot(n_dxyPhiMeanBiasTrend ,n_dxyPhiBiasResiduals,"mean","phi"); 
01099     FillTrendPlot(n_dxyPhiWidthBiasTrend,n_dxyPhiBiasResiduals,"width","phi");
01100     FillTrendPlot(n_dzPhiMeanBiasTrend  ,n_dzPhiBiasResiduals ,"mean","phi"); 
01101     FillTrendPlot(n_dzPhiWidthBiasTrend ,n_dzPhiBiasResiduals ,"width","phi");
01102     
01103     FillTrendPlot(n_dxyEtaMeanBiasTrend ,n_dxyEtaBiasResiduals,"mean","eta"); 
01104     FillTrendPlot(n_dxyEtaWidthBiasTrend,n_dxyEtaBiasResiduals,"width","eta");
01105     FillTrendPlot(n_dzEtaMeanBiasTrend  ,n_dzEtaBiasResiduals ,"mean","eta"); 
01106     FillTrendPlot(n_dzEtaWidthBiasTrend ,n_dzEtaBiasResiduals ,"width","eta");
01107     
01108     // medians and MADs   
01109     
01110     FillTrendPlot(a_dxyPhiMedianBiasTrend,a_dxyPhiBiasResiduals,"median","phi");  
01111     FillTrendPlot(a_dxyPhiMADBiasTrend   ,a_dxyPhiBiasResiduals,"mad","phi"); 
01112     FillTrendPlot(a_dzPhiMedianBiasTrend ,a_dzPhiBiasResiduals ,"median","phi");  
01113     FillTrendPlot(a_dzPhiMADBiasTrend    ,a_dzPhiBiasResiduals ,"mad","phi"); 
01114     
01115     FillTrendPlot(a_dxyEtaMedianBiasTrend,a_dxyEtaBiasResiduals,"median","eta");  
01116     FillTrendPlot(a_dxyEtaMADBiasTrend   ,a_dxyEtaBiasResiduals,"mad","eta"); 
01117     FillTrendPlot(a_dzEtaMedianBiasTrend ,a_dzEtaBiasResiduals ,"median","eta");  
01118     FillTrendPlot(a_dzEtaMADBiasTrend    ,a_dzEtaBiasResiduals ,"mad","eta"); 
01119     
01120     FillTrendPlot(n_dxyPhiMedianBiasTrend,n_dxyPhiBiasResiduals,"median","phi");  
01121     FillTrendPlot(n_dxyPhiMADBiasTrend   ,n_dxyPhiBiasResiduals,"mad","phi"); 
01122     FillTrendPlot(n_dzPhiMedianBiasTrend ,n_dzPhiBiasResiduals ,"median","phi");  
01123     FillTrendPlot(n_dzPhiMADBiasTrend    ,n_dzPhiBiasResiduals ,"mad","phi"); 
01124     
01125     FillTrendPlot(n_dxyEtaMedianBiasTrend,n_dxyEtaBiasResiduals,"median","eta");  
01126     FillTrendPlot(n_dxyEtaMADBiasTrend   ,n_dxyEtaBiasResiduals,"mad","eta"); 
01127     FillTrendPlot(n_dzEtaMedianBiasTrend ,n_dzEtaBiasResiduals ,"median","eta");  
01128     FillTrendPlot(n_dzEtaMADBiasTrend    ,n_dzEtaBiasResiduals ,"mad","eta"); 
01129     
01130   }
01131 
01132   FillTrendPlot(a_dxyPhiMeanTrend ,a_dxyPhiResiduals,"mean","phi");  
01133   FillTrendPlot(a_dxyPhiWidthTrend,a_dxyPhiResiduals,"width","phi");
01134   FillTrendPlot(a_dzPhiMeanTrend  ,a_dzPhiResiduals ,"mean","phi");   
01135   FillTrendPlot(a_dzPhiWidthTrend ,a_dzPhiResiduals ,"width","phi");  
01136   
01137   FillTrendPlot(a_dxyEtaMeanTrend ,a_dxyEtaResiduals,"mean","eta"); 
01138   FillTrendPlot(a_dxyEtaWidthTrend,a_dxyEtaResiduals,"width","eta");
01139   FillTrendPlot(a_dzEtaMeanTrend  ,a_dzEtaResiduals ,"mean","eta"); 
01140   FillTrendPlot(a_dzEtaWidthTrend ,a_dzEtaResiduals ,"width","eta");
01141   
01142   FillTrendPlot(n_dxyPhiMeanTrend ,n_dxyPhiResiduals,"mean","phi"); 
01143   FillTrendPlot(n_dxyPhiWidthTrend,n_dxyPhiResiduals,"width","phi");
01144   FillTrendPlot(n_dzPhiMeanTrend  ,n_dzPhiResiduals ,"mean","phi"); 
01145   FillTrendPlot(n_dzPhiWidthTrend ,n_dzPhiResiduals ,"width","phi");
01146   
01147   FillTrendPlot(n_dxyEtaMeanTrend ,n_dxyEtaResiduals,"mean","eta"); 
01148   FillTrendPlot(n_dxyEtaWidthTrend,n_dxyEtaResiduals,"width","eta");
01149   FillTrendPlot(n_dzEtaMeanTrend  ,n_dzEtaResiduals ,"mean","eta"); 
01150   FillTrendPlot(n_dzEtaWidthTrend ,n_dzEtaResiduals ,"width","eta");
01151     
01152   // medians and MADs     
01153   
01154   FillTrendPlot(a_dxyPhiMedianTrend,a_dxyPhiResiduals,"median","phi");  
01155   FillTrendPlot(a_dxyPhiMADTrend   ,a_dxyPhiResiduals,"mad","phi"); 
01156   FillTrendPlot(a_dzPhiMedianTrend ,a_dzPhiResiduals ,"median","phi");  
01157   FillTrendPlot(a_dzPhiMADTrend    ,a_dzPhiResiduals ,"mad","phi"); 
01158   
01159   FillTrendPlot(a_dxyEtaMedianTrend,a_dxyEtaResiduals,"median","eta");  
01160   FillTrendPlot(a_dxyEtaMADTrend   ,a_dxyEtaResiduals,"mad","eta"); 
01161   FillTrendPlot(a_dzEtaMedianTrend ,a_dzEtaResiduals ,"median","eta");  
01162   FillTrendPlot(a_dzEtaMADTrend    ,a_dzEtaResiduals ,"mad","eta"); 
01163   
01164   FillTrendPlot(n_dxyPhiMedianTrend,n_dxyPhiResiduals,"median","phi");  
01165   FillTrendPlot(n_dxyPhiMADTrend   ,n_dxyPhiResiduals,"mad","phi"); 
01166   FillTrendPlot(n_dzPhiMedianTrend ,n_dzPhiResiduals ,"median","phi");  
01167   FillTrendPlot(n_dzPhiMADTrend    ,n_dzPhiResiduals ,"mad","phi"); 
01168   
01169   FillTrendPlot(n_dxyEtaMedianTrend,n_dxyEtaResiduals,"median","eta");  
01170   FillTrendPlot(n_dxyEtaMADTrend   ,n_dxyEtaResiduals,"mad","eta"); 
01171   FillTrendPlot(n_dzEtaMedianTrend ,n_dzEtaResiduals ,"median","eta");  
01172   FillTrendPlot(n_dzEtaMADTrend    ,n_dzEtaResiduals ,"mad","eta"); 
01173     
01174 }
01175 
01176 //*************************************************************
01177 void PrimaryVertexValidation::SetVarToZero() 
01178 //*************************************************************
01179 {
01180   nTracks_ = 0;
01181   nClus_ = 0;
01182   nOfflineVertices_=0;
01183   RunNumber_ =0;
01184   LuminosityBlockNumber_=0;
01185   xOfflineVertex_ =-999.;
01186   yOfflineVertex_ =-999.;
01187   zOfflineVertex_ =-999.;
01188   BSx0_ = -999.;
01189   BSy0_ = -999.;
01190   BSz0_ = -999.;
01191   Beamsigmaz_=-999.;
01192   Beamdxdz_=-999.;   
01193   BeamWidthX_=-999.;
01194   BeamWidthY_=-999.;
01195 
01196   for ( int i=0; i<nMaxtracks_; ++i ) {
01197     
01198     pt_[i]        = 0;
01199     p_[i]         = 0;
01200     nhits_[i]     = 0;
01201     nhits1D_[i]   = 0;
01202     nhits2D_[i]   = 0;
01203     nhitsBPIX_[i] = 0;
01204     nhitsFPIX_[i] = 0;
01205     nhitsTIB_[i]  = 0;
01206     nhitsTID_[i]  = 0;
01207     nhitsTOB_[i]  = 0;
01208     nhitsTEC_[i]  = 0;
01209     isHighPurity_[i] = 0;
01210     eta_[i]       = 0;
01211     theta_[i]     = 0;
01212     phi_[i]       = 0;
01213     chi2_[i]      = 0;
01214     chi2ndof_[i]  = 0;
01215     charge_[i]    = 0;
01216     qoverp_[i]    = 0;
01217     dz_[i]        = 0;
01218     dxy_[i]       = 0;
01219     dzBs_[i]      = 0;
01220     dxyBs_[i]     = 0;
01221     xPCA_[i]      = 0;
01222     yPCA_[i]      = 0;
01223     zPCA_[i]      = 0;
01224     xUnbiasedVertex_[i] =0;    
01225     yUnbiasedVertex_[i] =0;
01226     zUnbiasedVertex_[i] =0;
01227     chi2normUnbiasedVertex_[i]=0;
01228     chi2UnbiasedVertex_[i]=0;
01229     DOFUnbiasedVertex_[i]=0;
01230     sumOfWeightsUnbiasedVertex_[i]=0;
01231     tracksUsedForVertexing_[i]=0;
01232     dxyFromMyVertex_[i]=0;
01233     dzFromMyVertex_[i]=0;
01234     dszFromMyVertex_[i]=0;
01235     dxyErrorFromMyVertex_[i]=0; 
01236     dzErrorFromMyVertex_[i]=0; 
01237     IPTsigFromMyVertex_[i]=0;   
01238     IPLsigFromMyVertex_[i]=0;    
01239     hasRecVertex_[i] = 0;
01240     isGoodTrack_[i]  = 0;
01241   } 
01242 }
01243 
01244 //*************************************************************
01245 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMedian(TH1F *histo)
01246 //*************************************************************
01247 {
01248   Double_t median = 999;
01249   int nbins = histo->GetNbinsX();
01250 
01251   //extract median from histogram
01252   double *x = new double[nbins];
01253   double *y = new double[nbins];
01254   for (int j = 0; j < nbins; j++) {
01255     x[j] = histo->GetBinCenter(j+1);
01256     y[j] = histo->GetBinContent(j+1);
01257   }
01258   median = TMath::Median(nbins, x, y);
01259   
01260   delete[] x; x = 0;
01261   delete [] y; y = 0;  
01262 
01263   std::pair<Double_t,Double_t> result;
01264   result = std::make_pair(median,median/TMath::Sqrt(histo->GetEntries()));
01265 
01266   return result;
01267 
01268 }
01269 
01270 //*************************************************************
01271 std::pair<Double_t,Double_t> PrimaryVertexValidation::getMAD(TH1F *histo)
01272 //*************************************************************
01273 {
01274 
01275   int nbins = histo->GetNbinsX();
01276   Double_t median = getMedian(histo).first;
01277   Double_t x_lastBin = histo->GetBinLowEdge(nbins+1);
01278   const char *HistoName =histo->GetName();
01279   TString Finalname = Form("resMed%s",HistoName);
01280   TH1F *newHisto = new TH1F(Finalname,Finalname,nbins,0.,x_lastBin);
01281   Double_t *residuals = new Double_t[nbins];
01282   Double_t *weights = new Double_t[nbins];
01283 
01284   for (int j = 0; j < nbins; j++) {
01285     residuals[j] = TMath::Abs(median - histo->GetBinCenter(j+1));
01286     weights[j]=histo->GetBinContent(j+1);
01287     newHisto->Fill(residuals[j],weights[j]);
01288   }
01289   
01290   Double_t theMAD = (getMedian(newHisto).first)*1.4826;
01291   newHisto->Delete("");
01292   
01293   std::pair<Double_t,Double_t> result;
01294   result = std::make_pair(theMAD,theMAD/histo->GetEntries());
01295 
01296   return result;
01297 
01298 }
01299 
01300 //*************************************************************
01301 std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t>  > PrimaryVertexValidation::fitResiduals(TH1 *hist)
01302 //*************************************************************
01303 {
01304   //float fitResult(9999);
01305   //if (hist->GetEntries() < 20) return ;
01306   
01307   float mean  = hist->GetMean();
01308   float sigma = hist->GetRMS();
01309   
01310   TF1 func("tmp", "gaus", mean - 1.5*sigma, mean + 1.5*sigma); 
01311   if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
01312     mean  = func.GetParameter(1);
01313     sigma = func.GetParameter(2);
01314     // second fit: three sigma of first fit around mean of first fit
01315     func.SetRange(mean - 2*sigma, mean + 2*sigma);
01316       // I: integral gives more correct results if binning is too wide
01317       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
01318     if (0 == hist->Fit(&func, "Q0LR")) {
01319       if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
01320         hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
01321       }
01322     }
01323   }
01324 
01325   float res_mean  = func.GetParameter(1);
01326   float res_width = func.GetParameter(2);
01327   
01328   float res_mean_err  = func.GetParError(1);
01329   float res_width_err = func.GetParError(2);
01330 
01331   std::pair<Double_t,Double_t> resultM;
01332   std::pair<Double_t,Double_t> resultW;
01333 
01334   resultM = std::make_pair(res_mean,res_mean_err);
01335   resultW = std::make_pair(res_width,res_width_err);
01336 
01337   std::pair<std::pair<Double_t,Double_t>, std::pair<Double_t,Double_t>  > result;
01338   
01339   result = std::make_pair(resultM,resultW);
01340   return result;
01341 }
01342 
01343 //*************************************************************
01344 void PrimaryVertexValidation::FillTrendPlot(TH1F* trendPlot, TH1F* residualsPlot[100], TString fitPar_, TString var_)
01345 //*************************************************************
01346 {
01347  
01348   float phiInterval = (360.)/nBins_;
01349   float etaInterval = 5./nBins_;
01350   
01351   for ( int i=0; i<nBins_; ++i ) {
01352     
01353     char phipositionString[129];
01354     float phiposition = (-180+i*phiInterval)+(phiInterval/2);
01355     sprintf(phipositionString,"%.f",phiposition);
01356     
01357     char etapositionString[129];
01358     float etaposition = (-2.5+i*etaInterval)+(etaInterval/2);
01359     sprintf(etapositionString,"%.1f",etaposition);
01360     
01361     if(fitPar_=="mean"){
01362       float mean_      = fitResiduals(residualsPlot[i]).first.first;
01363       float meanErr_   = fitResiduals(residualsPlot[i]).first.second;
01364       trendPlot->SetBinContent(i+1,mean_);
01365       trendPlot->SetBinError(i+1,meanErr_);
01366     } else if (fitPar_=="width"){
01367       float width_     = fitResiduals(residualsPlot[i]).second.first;
01368       float widthErr_  = fitResiduals(residualsPlot[i]).second.second;
01369       trendPlot->SetBinContent(i+1,width_);
01370       trendPlot->SetBinError(i+1,widthErr_);
01371     } else if (fitPar_=="median"){
01372       float median_    = getMedian(residualsPlot[i]).first;
01373       float medianErr_ = getMedian(residualsPlot[i]).second;
01374       trendPlot->SetBinContent(i+1,median_);
01375       trendPlot->SetBinError(i+1,medianErr_);
01376     } else if (fitPar_=="mad"){
01377       float mad_       = getMAD(residualsPlot[i]).first; 
01378       float madErr_    = getMAD(residualsPlot[i]).second;
01379       trendPlot->SetBinContent(i+1,mad_);
01380       trendPlot->SetBinError(i+1,madErr_);
01381     } else {
01382       std::cout<<"PrimaryVertexValidation::FillTrendPlot() "<<fitPar_<<" unknown estimator!"<<std::endl;
01383     }
01384 
01385     if(var_=="eta"){
01386       trendPlot->GetXaxis()->SetBinLabel(i+1,phipositionString); 
01387     } else if(var_=="phi"){
01388       trendPlot->GetXaxis()->SetBinLabel(i+1,etapositionString); 
01389     } else {
01390       std::cout<<"PrimaryVertexValidation::FillTrendPlot() "<<var_<<" unknown track parameter!"<<std::endl;
01391     }
01392   }
01393 }
01394 
01395 
01396 
01397 //define this as a plug-in
01398 DEFINE_FWK_MODULE(PrimaryVertexValidation);
Generated for CMSSW Reference Manual by doxygen 1.7.3