CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    TrackerOfflineValidation
00004 // Class:      TrackerOfflineValidation
00005 // 
00013 //
00014 // Original Author:  Erik Butz
00015 //         Created:  Tue Dec 11 14:03:05 CET 2007
00016 // $Id: TrackerOfflineValidation.cc,v 1.5 2008/05/26 16:09:24 ebutz Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 #include <map>
00024 #include <sstream>
00025 #include <math.h>
00026 
00027 // ROOT includes
00028 #include "TH1.h"
00029 #include "TH2.h"
00030 #include "TProfile.h"
00031 #include "TFile.h"
00032 
00033 // user include files
00034 #include "FWCore/Framework/interface/Frameworkfwd.h"
00035 #include "FWCore/Framework/interface/EDAnalyzer.h"
00036 #include "FWCore/Framework/interface/EventSetup.h"
00037 #include "FWCore/Framework/interface/Event.h"
00038 #include "FWCore/Framework/interface/MakerMacros.h"
00039 #include "FWCore/Framework/interface/ESHandle.h"
00040 #include "Alignment/OfflineValidation/interface/TrackerValidationVariables.h"
00041 #include "DataFormats/DetId/interface/DetId.h"
00042 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00043 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00044 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00045 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00046 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00047 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00048 #include "FWCore/ServiceRegistry/interface/Service.h"
00049 #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
00050 #include "PhysicsTools/UtilAlgos/interface/TFileDirectory.h"
00051 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00052 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
00053 #include "Alignment/CommonAlignment/interface/AlignableComposite.h"
00054 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00055 #include "Alignment/CommonAlignment/interface/Utilities.h"
00056 #include "Alignment/CommonAlignment/interface/AlignableObjectId.h"
00057 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
00058 
00059 //
00060 // class decleration
00061 //
00062 
00063 //typedef std::vector<DetId>          DetIdContainer;
00064 
00065 
00066 
00067 class TrackerOfflineValidation : public edm::EDAnalyzer {
00068 public:
00069   explicit TrackerOfflineValidation(const edm::ParameterSet&);
00070   ~TrackerOfflineValidation();
00071  
00072   
00073 private:
00074   virtual void beginJob(const edm::EventSetup&) ;
00075   virtual void analyze(const edm::Event&, const edm::EventSetup&);
00076   virtual void endJob();
00077   const edm::ParameterSet parset_;
00078   edm::ESHandle<TrackerGeometry> tkGeom_;
00079   edm::ParameterSet parameters_;
00080   
00081   struct ModuleHistos{
00082     ModuleHistos() :  ResHisto(), NormResHisto(), ResXprimeHisto(), NormResXprimeHisto() {} 
00083     TH1* ResHisto;
00084     TH1* NormResHisto;
00085     TH1* ResXprimeHisto;
00086     TH1* NormResXprimeHisto;
00087   };
00088   
00089   ModuleHistos& GetHistStructFromMap(const DetId& detid);
00090 
00091 
00092   std::vector<TH1*> vTrackHistos_;
00093   std::vector<TProfile*> vTrackProfiles_;
00094   std::vector<TH2*> vTrack2DHistos_;
00095   
00096   std::vector<TH1*> vSubdetRes_;
00097   std::vector<TH1*> vSubdetNormRes_;
00098   std::vector<TH1*> vSubdetXprimeRes_;
00099 
00100   //__gnu_cxx::hash_map<int,TrackerOfflineValidation::ModuleHistos> m_modulehistos;
00101   
00102   std::map<int,TrackerOfflineValidation::ModuleHistos> mPxbResiduals_;
00103   std::map<int,TrackerOfflineValidation::ModuleHistos> mPxeResiduals_;
00104   std::map<int,TrackerOfflineValidation::ModuleHistos> mTibResiduals_;
00105   std::map<int,TrackerOfflineValidation::ModuleHistos> mTidResiduals_;
00106   std::map<int,TrackerOfflineValidation::ModuleHistos> mTobResiduals_;
00107   std::map<int,TrackerOfflineValidation::ModuleHistos> mTecResiduals_;
00108 
00109   
00110   void bookGlobalHists(TFileDirectory &tfd);
00111   void bookDirHists(TFileDirectory &tfd, const Alignable& ali, const AlignableObjectId &aliobjid);
00112   void bookHists(TFileDirectory &tfd, const Alignable& ali, align::StructureType type, int i, const AlignableObjectId &aliobjid);
00113  
00114   void collateSummaryHists( TFileDirectory &tfd, const Alignable& ali, int i, const AlignableObjectId &aliobjid, std::vector<std::pair<TH1*,TH1*> > &v_levelProfiles);
00115   
00116   std::pair<TH1*,TH1*> bookSummaryHists(TFileDirectory &tfd, const Alignable& ali, align::StructureType type, int i, const AlignableObjectId &aliobjid); 
00117  
00118  
00119   float Fwhm(const TH1* hist);
00120 
00121   // From MillePedeAlignmentMonitor: Get Index for Arbitary vector<class> by name
00122   template <class OBJECT_TYPE>  
00123   int GetIndex(const std::vector<OBJECT_TYPE*> &vec, const TString &name);
00124   
00125  
00126   
00127   // ----------member data ---------------------------
00128 };
00129 
00130 //
00131 // constants, enums and typedefs
00132 //
00133 typedef std::vector<DetId>  DetIdContainer;
00134 //
00135 // static data member definitions
00136 //
00137 template <class OBJECT_TYPE>  
00138 int TrackerOfflineValidation::GetIndex(const std::vector<OBJECT_TYPE*> &vec, const TString &name)
00139 {
00140   int result = 0;
00141   for (typename std::vector<OBJECT_TYPE*>::const_iterator iter = vec.begin(), iterEnd = vec.end();
00142        iter != iterEnd; ++iter, ++result) {
00143     if (*iter && (*iter)->GetName() == name) return result;
00144   }
00145   edm::LogError("Alignment") << "@SUB=TrackerOfflineValidation::GetIndex" << " could not find " << name;
00146   return -1;
00147 }
00148 //
00149 // constructors and destructor
00150 //
00151 TrackerOfflineValidation::TrackerOfflineValidation(const edm::ParameterSet& iConfig)
00152   :   parset_(iConfig)
00153 {
00154    //now do what ever initialization is needed
00155  
00156 
00157 }
00158 
00159 
00160 TrackerOfflineValidation::~TrackerOfflineValidation()
00161 {
00162  
00163    // do anything here that needs to be done at desctruction time
00164    // (e.g. close files, deallocate resources etc.)
00165 
00166 }
00167 
00168 
00169 //
00170 // member functions
00171 //
00172 
00173 
00174 // ------------ method called once each job just before starting event loop  ------------
00175 void 
00176 TrackerOfflineValidation::beginJob(const edm::EventSetup& es)
00177 {
00178   es.get<TrackerDigiGeometryRecord>().get( tkGeom_ );
00179   edm::Service<TFileService> fs;    
00180   AlignableObjectId aliobjid;
00181 
00182   // construct alignable tracker to get access to alignable hierarchy 
00183   // -> no backward compatibility below 1_8_X
00184 
00185   AlignableTracker aliTracker(&(*tkGeom_));
00186   
00187 
00188   //
00189   // Book Histogramms for global track quantities
00190   TFileDirectory trackglobal = fs->mkdir("GlobalTrackVariables");
00191   
00192   this->bookGlobalHists(trackglobal);
00193   
00194   parameters_ = parset_.getParameter<edm::ParameterSet>("TH1ResModules");
00195   int32_t i_residuals_Nbins =  parameters_.getParameter<int32_t>("Nbinx");
00196   double d_residual_xmin = parameters_.getParameter<double>("xmin");
00197   double d_residual_xmax = parameters_.getParameter<double>("xmax");
00198 
00199   parameters_ = parset_.getParameter<edm::ParameterSet>("TH1NormResModules");
00200   int32_t i_normres_Nbins =  parameters_.getParameter<int32_t>("Nbinx");
00201   double d_normres_xmin = parameters_.getParameter<double>("xmin");
00202   double d_normres_xmax = parameters_.getParameter<double>("xmax");
00203   
00204 
00205   vSubdetRes_.push_back(fs->make<TH1F>("h_Residuals_pxb","Residuals in PXB;x_{pred} - x_{rec} [cm]", 
00206                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00207   vSubdetRes_.push_back(fs->make<TH1F>("h_Residuals_pxe","Residuals in PXE;x_{pred} - x_{rec} [cm]", 
00208                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00209   vSubdetRes_.push_back(fs->make<TH1F>("h_Residuals_tib","Residuals in TIB;x_{pred} - x_{rec} [cm]", 
00210                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00211   vSubdetRes_.push_back(fs->make<TH1F>("h_Residuals_tid","Residuals in TID;x_{pred} - x_{rec} [cm]", 
00212                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00213   vSubdetRes_.push_back(fs->make<TH1F>("h_Residuals_tob","Residuals in TOB;x_{pred} - x_{rec} [cm]", 
00214                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00215   vSubdetRes_.push_back(fs->make<TH1F>("h_Residuals_tec","Residuals in TEC;x_{pred} - x_{rec} [cm]", 
00216                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00217 
00218   vSubdetXprimeRes_.push_back(fs->make<TH1F>("h_XprimeResiduals_pxb","Residuals in PXB;x_{pred} - x_{rec} [cm]", 
00219                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00220   vSubdetXprimeRes_.push_back(fs->make<TH1F>("h_XprimeResiduals_pxe","Residuals in PXE;x_{pred} - x_{rec} [cm]", 
00221                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00222   vSubdetXprimeRes_.push_back(fs->make<TH1F>("h_XprimeResiduals_tib","Residuals in TIB;x_{pred} - x_{rec} [cm]", 
00223                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00224   vSubdetXprimeRes_.push_back(fs->make<TH1F>("h_XprimeResiduals_tid","Residuals in TID;x_{pred} - x_{rec} [cm]", 
00225                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00226   vSubdetXprimeRes_.push_back(fs->make<TH1F>("h_XprimeResiduals_tob","Residuals in TOB;x_{pred} - x_{rec} [cm]", 
00227                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00228   vSubdetXprimeRes_.push_back(fs->make<TH1F>("h_XprimeResiduals_tec","Residuals in TEC;x_{pred} - x_{rec} [cm]", 
00229                                        i_residuals_Nbins,d_residual_xmin,d_residual_xmax)) ;
00230 
00231 
00232   vSubdetNormRes_.push_back(fs->make<TH1F>("h_normResiduals_pxb","Normalized Residuals in PXB;(x_{pred} - x_{rec})/#sqrt{V}", 
00233                                            i_normres_Nbins,d_normres_xmin,d_normres_xmax));   
00234   vSubdetNormRes_.push_back(fs->make<TH1F>("h_normResiduals_pxe","Normalized Residuals in PXE;(x_{pred} - x_{rec})/#sqrt{V}", 
00235                                            i_normres_Nbins,d_normres_xmin,d_normres_xmax));   
00236   vSubdetNormRes_.push_back(fs->make<TH1F>("h_normResiduals_tib","Normalized Residuals in TIB;(x_{pred} - x_{rec})/#sqrt{V}", 
00237                                            i_normres_Nbins,d_normres_xmin,d_normres_xmax));   
00238   vSubdetNormRes_.push_back(fs->make<TH1F>("h_normResiduals_tid","Normalized Residuals in TID;(x_{pred} - x_{rec})/#sqrt{V}", 
00239                                            i_normres_Nbins,d_normres_xmin,d_normres_xmax));   
00240   vSubdetNormRes_.push_back(fs->make<TH1F>("h_normResiduals_tob","Normalized Residuals in TOB;(x_{pred} - x_{rec})/#sqrt{V}", 
00241                                            i_normres_Nbins,d_normres_xmin,d_normres_xmax));   
00242   vSubdetNormRes_.push_back(fs->make<TH1F>("h_normResiduals_tec","Normalized Residuals in TEC;(x_{pred} - x_{rec})/#sqrt{V}", 
00243                                            i_normres_Nbins,d_normres_xmin,d_normres_xmax));
00244 
00245 
00246   
00247 
00248   
00249   edm::LogInfo("TrackerOfflineValidation") << "There are " << (*tkGeom_).detIds().size() << " detUnits in the Geometry record" << std::endl;
00250 
00251   // recursively book histogramms on lowest level
00252   this->bookDirHists(static_cast<TFileDirectory&>(*fs), aliTracker, aliobjid);
00253   
00254  
00255 
00256 }
00257 
00258 void 
00259 TrackerOfflineValidation::bookGlobalHists(TFileDirectory &tfd )
00260 {
00261 
00262   vTrackHistos_.push_back(tfd.make<TH1F>("h_tracketa","Track #eta;#eta_{Track};Number of Tracks",90,-3.,3.));          
00263   vTrackHistos_.push_back(tfd.make<TH1F>("h_curvature","Curvature #kappa;#kappa_{Track};Number of Tracks",100,-.05,.05));
00264   vTrackHistos_.push_back(tfd.make<TH1F>("h_curvature_pos","Curvature |#kappa| Positive Tracks;|#kappa_{pos Track}|;Number of Tracks",100,.0,.05));
00265   vTrackHistos_.push_back(tfd.make<TH1F>("h_curvature_neg","Curvature |#kappa| Negative Tracks;|#kappa_{neg Track}|;Number of Tracks",100,.0,.05));
00266   vTrackHistos_.push_back(tfd.make<TH1F>("h_diff_curvature","Curvature |#kappa| Tracks Difference;|#kappa_{Track}|;# Pos Tracks - # Neg Tracks",100,.0,.05));
00267   vTrackHistos_.push_back(tfd.make<TH1F>("h_chi2","#chi^{2};#chi^{2}_{Track};Number of Tracks",500,-0.01,500.));               
00268   vTrackHistos_.push_back(tfd.make<TH1F>("h_normchi2","#chi^{2}/ndof;#chi^{2}/ndof;Number of Tracks",100,-0.01,10.));     
00269   vTrackHistos_.push_back(tfd.make<TH1F>("h_pt","p_{T};p_{T}^{track};Number of Tracks",100,0.,2500));                     
00270 
00271   vTrackProfiles_.push_back(tfd.make<TProfile>("h_d0_vs_phi","Transverse Impact Parameter vs. #phi;#phi_{Track};#LT d_{0} #GT",100,-3.15,3.15));  
00272   vTrackProfiles_.push_back(tfd.make<TProfile>("h_dz_vs_phi","Longitudinal Impact Parameter vs. #phi;#phi_{Track};#LT d_{z} #GT",100,-3.15,3.15));
00273   vTrackProfiles_.push_back(tfd.make<TProfile>("h_d0_vs_eta","Transverse Impact Parameter vs. #eta;#eta_{Track};#LT d_{0} #GT",100,-3.15,3.15));  
00274   vTrackProfiles_.push_back(tfd.make<TProfile>("h_dz_vs_eta","Longitudinal Impact Parameter vs. #eta;#eta_{Track};#LT d_{z} #GT",100,-3.15,3.15));
00275   vTrackProfiles_.push_back(tfd.make<TProfile>("h_chi2_vs_phi","#chi^{2} vs. #phi;#phi_{Track};#LT #chi^{2} #GT",100,-3.15,3.15));                
00276   vTrackProfiles_.push_back(tfd.make<TProfile>("h_normchi2_vs_phi","#chi^{2}/ndof vs. #phi;#phi_{Track};#LT #chi^{2}/ndof #GT",100,-3.15,3.15));  
00277   vTrackProfiles_.push_back(tfd.make<TProfile>("h_chi2_vs_eta","#chi^{2} vs. #eta;#eta_{Track};#LT #chi^{2} #GT",100,-3.15,3.15));                
00278   vTrackProfiles_.push_back(tfd.make<TProfile>("h_normchi2_vs_eta","#chi^{2}/ndof vs. #eta;#eta_{Track};#LT #chi^{2}/ndof #GT",100,-3.15,3.15));  
00279 
00280   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_d0_vs_phi","Transverse Impact Parameter vs. #phi;#phi_{Track};d_{0}",100, -3.15, 3.15, 100,-1.,1.) );  
00281   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_dz_vs_phi","Longitudinal Impact Parameter vs. #phi;#phi_{Track};d_{z}",100, -3.15, 3.15, 100,-100.,100.));
00282   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_d0_vs_eta","Transverse Impact Parameter vs. #eta;#eta_{Track};d_{0}",100, -3.15, 3.15, 100,-1.,1.));  
00283   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_dz_vs_eta","Longitudinal Impact Parameter vs. #eta;#eta_{Track};d_{z}",100, -3.15, 3.15, 100,-100.,100.));
00284   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2_vs_phi","#chi^{2} vs. #phi;#phi_{Track};#chi^{2}",100, -3.15, 3.15, 500, 0., 500.));                  
00285   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_normchi2_vs_phi","#chi^{2}/ndof vs. #phi;#phi_{Track};#chi^{2}/ndof",100, -3.15, 3.15, 100, 0., 10.));  
00286   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_chi2_vs_eta","#chi^{2} vs. #eta;#eta_{Track};#chi^{2}",100, -3.15, 3.15, 500, 0., 500.));                  
00287   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_normchi2_vs_eta","#chi^{2}/ndof vs. #eta;#eta_{Track};#chi^{2}/ndof",100,-3.15,3.15, 100, 0., 10.));  
00288   vTrack2DHistos_.push_back(tfd.make<TH2F>("h2_kappa_vs_phi","#kappa vs. #phi;#phi_{Track};#kappa",100,-3.15,3.15, 100, .0,.05));  
00289   
00290 
00291 
00292 
00293 }
00294 
00295 
00296 void
00297 TrackerOfflineValidation::bookDirHists( TFileDirectory &tfd, const Alignable& ali, const AlignableObjectId &aliobjid)
00298 {
00299   std::vector<Alignable*> alivec(ali.components());
00300   for(int i=0, iEnd = ali.components().size();i < iEnd; ++i) {
00301     std::string structurename  = aliobjid.typeToName((alivec)[i]->alignableObjectId());
00302     LogDebug("TrackerOfflineValidation") << "StructureName = " << structurename;
00303     std::stringstream dirname;
00304     
00305     // add no suffix counter to Strip and Pixel
00306     // just aesthetics
00307     if(structurename != "Strip" && structurename != "Pixel") {
00308       dirname << structurename << "_" << i+1;
00309     } else {
00310       dirname << structurename;
00311     }
00312     if (structurename.find("Endcap",0) != std::string::npos )   {
00313       TFileDirectory f = tfd.mkdir((dirname.str()).c_str());
00314       bookHists(f, *(alivec)[i], ali.alignableObjectId() , i, aliobjid);
00315       bookDirHists( f, *(alivec)[i], aliobjid);
00316     } else if((structurename != "Det" && structurename != "DetUnit" ) || alivec[i]->components().size() > 1) {      
00317       TFileDirectory f = tfd.mkdir((dirname.str()).c_str());
00318       bookHists(tfd, *(alivec)[i], ali.alignableObjectId() , i, aliobjid);
00319       bookDirHists( f, *(alivec)[i], aliobjid);
00320     } else {
00321       bookHists(tfd, *(alivec)[i], ali.alignableObjectId() , i, aliobjid);
00322     }
00323   }
00324 }
00325 
00326 
00327 void 
00328 TrackerOfflineValidation::bookHists(TFileDirectory &tfd, const Alignable& ali, align::StructureType type, int i, const AlignableObjectId &aliobjid)
00329 {
00330 
00331   // binnings for module level histogramms are steerable via cfg file
00332   parameters_ = parset_.getParameter<edm::ParameterSet>("TH1ResModules");
00333   int32_t i_residuals_Nbins =  parameters_.getParameter<int32_t>("Nbinx");
00334   double d_residual_xmin = parameters_.getParameter<double>("xmin");
00335   double d_residual_xmax = parameters_.getParameter<double>("xmax");
00336   parameters_ = parset_.getParameter<edm::ParameterSet>("TH1NormResModules");
00337   int32_t i_normres_Nbins =  parameters_.getParameter<int32_t>("Nbinx");
00338   double d_normres_xmin = parameters_.getParameter<double>("xmin");
00339   double d_normres_xmax = parameters_.getParameter<double>("xmax");
00340 
00341   TrackerAlignableId aliid;
00342   const DetId id = ali.id();
00343   
00344 
00345   // comparing subdetandlayer to subdetIds gives a warning at compile time
00346   // -> subdetandlayer could also be pair<uint,uint> but this has to be adapted
00347   // in AlignableObjId 
00348   std::pair<int,int> subdetandlayer = aliid.typeAndLayerFromDetId(id);
00349 
00350   align::StructureType subtype = align::invalid;
00351   
00352   // are we on or just above det, detunit level respectively?
00353   if (type == align::AlignableDetUnit )
00354     subtype = type;
00355   else if(      ali.alignableObjectId() == align::AlignableDet || 
00356                 ali.alignableObjectId() == align::AlignableDetUnit) 
00357     subtype = ali.alignableObjectId();
00358   
00359         
00360   std::stringstream histoname, histotitle, normhistoname, normhistotitle, xprimehistoname, xprimehistotitle;
00361   if( subdetandlayer.first == StripSubdetector::TID || subdetandlayer.first == StripSubdetector::TEC ||
00362       subdetandlayer.first == PixelSubdetector::PixelEndcap ) {
00363     histoname << "h_residuals_subdet_" << subdetandlayer.first 
00364               << "_wheel_" << subdetandlayer.second << "_module_" << id.rawId();
00365     xprimehistoname << "h_xprime_residuals_subdet_" << subdetandlayer.first 
00366               << "_wheel_" << subdetandlayer.second << "_module_" << id.rawId();
00367     normhistoname << "h_normresiduals_subdet_" << subdetandlayer.first 
00368                   << "_wheel_" << subdetandlayer.second << "_module_" << id.rawId();
00369     histotitle << "Residual for module " << id.rawId() << ";x_{pred} - x_{rec} [cm]";
00370     normhistotitle << "Normalized Residual for module " << id.rawId() << ";x_{pred} - x_{rec}/#sigma";
00371     xprimehistotitle << "X' Residual for module " << id.rawId() << ";x_{pred} - x_{rec} [cm]";
00372   } else if (subdetandlayer.first == StripSubdetector::TIB || subdetandlayer.first == StripSubdetector::TOB ||
00373              subdetandlayer.first == PixelSubdetector::PixelBarrel ) {
00374     histoname << "h_residuals_subdet_" << subdetandlayer.first 
00375               << "_layer_" << subdetandlayer.second << "_module_" << id.rawId();
00376     xprimehistoname << "h_xprime_residuals_subdet_" << subdetandlayer.first 
00377                     << "_layer_" << subdetandlayer.second << "_module_" << id.rawId();
00378     normhistoname << "h_normresiduals_subdet_" << subdetandlayer.first 
00379                   << "_layer_" << subdetandlayer.second << "_module_" << id.rawId();
00380     histotitle << "Residual for module " << id.rawId() << ";x_{pred} - x_{rec} [cm]";
00381     normhistotitle << "Normalized Residual for module " << id.rawId() << ";x_{pred} - x_{rec}/#sigma";
00382     xprimehistotitle << "X' Residual for module " << id.rawId() << ";x_{pred} - x_{rec} [cm]";
00383   } else {
00384     edm::LogWarning("TrackerOfflineValidation") << "@SUB=TrackerOfflineValidation::bookHists" 
00385                                                 << "Unknown subdetid: " <<  subdetandlayer.first; 
00386     
00387   }
00388   
00389   
00390   if(subtype == align::AlignableDet || subtype == align::AlignableDetUnit) {
00391     ModuleHistos &histStruct = this->GetHistStructFromMap(id);
00392     histStruct.ResHisto       =  tfd.make<TH1F>(histoname.str().c_str(),histotitle.str().c_str(),                    
00393                                          i_residuals_Nbins,d_residual_xmin,d_residual_xmax);
00394     histStruct.NormResHisto   =  tfd.make<TH1F>(normhistoname.str().c_str(),normhistotitle.str().c_str(),
00395                                                 i_normres_Nbins,d_normres_xmin,d_normres_xmax);
00396     histStruct.ResXprimeHisto =  tfd.make<TH1F>(xprimehistoname.str().c_str(),xprimehistotitle.str().c_str(),
00397                                                 i_residuals_Nbins,d_residual_xmin,d_residual_xmax);
00398   }
00399   
00400 }
00401 
00402 
00403 TrackerOfflineValidation::ModuleHistos& TrackerOfflineValidation::GetHistStructFromMap(const DetId& detid)
00404 {
00405 
00406   uint subdetid = detid.subdetId();
00407   if(subdetid == PixelSubdetector::PixelBarrel ) {
00408     return mPxbResiduals_[detid.rawId()];
00409   } else if (subdetid == PixelSubdetector::PixelEndcap) {      
00410     return mPxeResiduals_[detid.rawId()];
00411   } else if(subdetid  == StripSubdetector::TIB) {
00412     return mTibResiduals_[detid.rawId()];
00413   } else if(subdetid  == StripSubdetector::TID) {
00414     return mTidResiduals_[detid.rawId()];
00415   } else if(subdetid  == StripSubdetector::TOB) {
00416     return mTobResiduals_[detid.rawId()];
00417   } else if(subdetid  == StripSubdetector::TEC) {
00418     return mTecResiduals_[detid.rawId()];
00419   } else {
00420     throw cms::Exception("Geometry Error") 
00421       << "[TrackerOfflineValidation] Error, tried to get reference for non-tracker subdet " << subdetid 
00422       << " from detector " << detid.det();
00423     return mPxbResiduals_[0];
00424   }
00425   
00426 }
00427 
00428 
00429 // ------------ method called to for each event  ------------
00430 void
00431 TrackerOfflineValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00432 {
00433   
00434   using namespace edm;
00435   TrackerValidationVariables avalidator_(iSetup,parset_);
00436   
00437     
00438   std::vector<TrackerValidationVariables::AVTrackStruct> vTrackstruct;
00439   avalidator_.fillTrackQuantities(iEvent,vTrackstruct);
00440   std::vector<TrackerValidationVariables::AVHitStruct> v_hitstruct;
00441   avalidator_.fillHitQuantities(iEvent,v_hitstruct);
00442   
00443   
00444   for (std::vector<TrackerValidationVariables::AVTrackStruct>::const_iterator it = vTrackstruct.begin(),
00445          itEnd = vTrackstruct.end(); it != itEnd; ++it) {
00446     
00447     // Fill 1D track histos
00448     static const int etaindex = this->GetIndex(vTrackHistos_,"h_tracketa");
00449     vTrackHistos_[etaindex]->Fill(it->eta);
00450     static const int kappaindex = this->GetIndex(vTrackHistos_,"h_curvature");
00451     vTrackHistos_[kappaindex]->Fill(it->kappa);
00452     static const int kappaposindex = this->GetIndex(vTrackHistos_,"h_curvature_pos");
00453     if(it->charge > 0)
00454       vTrackHistos_[kappaposindex]->Fill(fabs(it->kappa));
00455     static const int kappanegindex = this->GetIndex(vTrackHistos_,"h_curvature_neg");
00456     if(it->charge < 0)
00457       vTrackHistos_[kappanegindex]->Fill(fabs(it->kappa));
00458     static const int normchi2index = this->GetIndex(vTrackHistos_,"h_normchi2");
00459     vTrackHistos_[normchi2index]->Fill(it->normchi2);
00460     static const int chi2index = this->GetIndex(vTrackHistos_,"h_chi2");
00461     vTrackHistos_[chi2index]->Fill(it->chi2);
00462     static const int ptindex = this->GetIndex(vTrackHistos_,"h_pt");
00463     vTrackHistos_[ptindex]->Fill(it->pt);
00464     
00465     // Fill track profiles
00466     static const int d0phiindex = this->GetIndex(vTrackProfiles_,"h_d0_vs_phi");
00467     vTrackProfiles_[d0phiindex]->Fill(it->phi,it->d0);
00468     static const int dzphiindex = this->GetIndex(vTrackProfiles_,"h_dz_vs_phi");
00469     vTrackProfiles_[dzphiindex]->Fill(it->phi,it->dz);
00470     static const int d0etaindex = this->GetIndex(vTrackProfiles_,"h_d0_vs_eta");
00471     vTrackProfiles_[d0etaindex]->Fill(it->eta,it->d0);
00472     static const int dzetaindex = this->GetIndex(vTrackProfiles_,"h_dz_vs_eta");
00473     vTrackProfiles_[dzetaindex]->Fill(it->eta,it->dz);
00474     static const int chiphiindex = this->GetIndex(vTrackProfiles_,"h_chi2_vs_phi");
00475     vTrackProfiles_[chiphiindex]->Fill(it->phi,it->chi2);
00476     static const int normchiphiindex = this->GetIndex(vTrackProfiles_,"h_normchi2_vs_phi");
00477     vTrackProfiles_[normchiphiindex]->Fill(it->phi,it->normchi2);
00478     static const int chietaindex = this->GetIndex(vTrackProfiles_,"h_chi2_vs_eta");
00479     vTrackProfiles_[chietaindex]->Fill(it->eta,it->chi2);
00480     static const int normchietaindex = this->GetIndex(vTrackProfiles_,"h_normchi2_vs_eta");
00481     vTrackProfiles_[normchietaindex]->Fill(it->eta,it->normchi2);
00482     
00483     // Fill 2D track histos
00484     static const int d0phiindex_2d = this->GetIndex(vTrack2DHistos_,"h2_d0_vs_phi");
00485     vTrack2DHistos_[d0phiindex_2d]->Fill(it->phi,it->d0);
00486     static const int dzphiindex_2d = this->GetIndex(vTrack2DHistos_,"h2_dz_vs_phi");
00487     vTrack2DHistos_[dzphiindex_2d]->Fill(it->phi,it->dz);
00488     static const int d0etaindex_2d = this->GetIndex(vTrack2DHistos_,"h2_d0_vs_eta");
00489     vTrack2DHistos_[d0etaindex_2d]->Fill(it->eta,it->d0);
00490     static const int dzetaindex_2d = this->GetIndex(vTrack2DHistos_,"h2_dz_vs_eta");
00491     vTrack2DHistos_[dzetaindex_2d]->Fill(it->eta,it->dz);
00492     static const int chiphiindex_2d = this->GetIndex(vTrack2DHistos_,"h2_chi2_vs_phi");
00493     vTrack2DHistos_[chiphiindex_2d]->Fill(it->phi,it->chi2);
00494     static const int normchiphiindex_2d = this->GetIndex(vTrack2DHistos_,"h2_normchi2_vs_phi");
00495     vTrack2DHistos_[normchiphiindex_2d]->Fill(it->phi,it->normchi2);
00496     static const int chietaindex_2d = this->GetIndex(vTrack2DHistos_,"h2_chi2_vs_eta");
00497     vTrack2DHistos_[chietaindex_2d]->Fill(it->eta,it->chi2);
00498     static const int normchietaindex_2d = this->GetIndex(vTrack2DHistos_,"h2_normchi2_vs_eta");
00499     vTrack2DHistos_[normchietaindex_2d]->Fill(it->eta,it->normchi2);
00500     static const int kappaphiindex_2d = this->GetIndex(vTrack2DHistos_,"h2_kappa_vs_phi");
00501     vTrack2DHistos_[kappaphiindex_2d]->Fill(it->phi,it->kappa);
00502      
00503   } // finish loop over track quantities
00504 
00505 
00506   // hit quantities: residuals, normalized residuals
00507   for (std::vector<TrackerValidationVariables::AVHitStruct>::const_iterator it = v_hitstruct.begin(),
00508          itEnd = v_hitstruct.end(); it != itEnd; ++it) {
00509     DetId detid(it->rawDetId);
00510     ModuleHistos &histStruct = this->GetHistStructFromMap(detid);
00511     histStruct.ResHisto->Fill(it->resX);
00512     if(it->resXprime != -999)  histStruct.ResXprimeHisto->Fill(it->resXprime);
00513     if(it->resErrX != 0)       histStruct.NormResHisto->Fill(it->resX/it->resErrX);
00514     
00515   }
00516   
00517 }
00518 
00519 
00520 
00521 // ------------ method called once each job just after ending the event loop  ------------
00522 void 
00523 TrackerOfflineValidation::endJob() {
00524   
00525   AlignableTracker aliTracker(&(*tkGeom_));
00526   edm::Service<TFileService> fs;   
00527   AlignableObjectId aliobjid;
00528 
00529 
00530   static const int kappadiffindex = this->GetIndex(vTrackHistos_,"h_diff_curvature");
00531   vTrackHistos_[kappadiffindex]->Add(vTrackHistos_[this->GetIndex(vTrackHistos_,"h_curvature_neg")],vTrackHistos_[this->GetIndex(vTrackHistos_,"h_curvature_pos")],-1,1);
00532   // Collate Information for Subdetectors
00533   // So far done by default, should be steerable
00534   for(std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator itPxb = mPxbResiduals_.begin(), 
00535         itEnd = mPxbResiduals_.end(); itPxb != itEnd;++itPxb ) {
00536     vSubdetRes_[PixelSubdetector::PixelBarrel - 1]->Add(itPxb->second.ResHisto); 
00537     vSubdetXprimeRes_[PixelSubdetector::PixelBarrel - 1]->Add(itPxb->second.ResXprimeHisto); 
00538     vSubdetNormRes_[PixelSubdetector::PixelBarrel - 1]->Add(itPxb->second.NormResHisto);
00539   }
00540   for(std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator itPxe = mPxeResiduals_.begin(), 
00541         itEnd = mPxeResiduals_.end(); itPxe != itEnd;++itPxe ) {
00542     vSubdetRes_[PixelSubdetector::PixelEndcap - 1]->Add(itPxe->second.ResHisto);
00543     vSubdetXprimeRes_[PixelSubdetector::PixelEndcap - 1]->Add(itPxe->second.ResXprimeHisto);
00544     vSubdetNormRes_[PixelSubdetector::PixelEndcap - 1]->Add(itPxe->second.NormResHisto);
00545   }
00546   for(std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator itTib = mTibResiduals_.begin(), 
00547         itEnd = mTibResiduals_.end(); itTib != itEnd;++itTib ) {
00548     vSubdetRes_[StripSubdetector::TIB - 1]->Add(itTib->second.ResHisto);
00549     vSubdetXprimeRes_[StripSubdetector::TIB - 1]->Add(itTib->second.ResXprimeHisto);
00550     vSubdetNormRes_[StripSubdetector::TIB -1]->Add(itTib->second.NormResHisto);
00551   }
00552   for(std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator itTid = mTidResiduals_.begin(), 
00553         itEnd = mTidResiduals_.end(); itTid != itEnd;++itTid ) {
00554     vSubdetRes_[StripSubdetector::TID - 1]->Add(itTid->second.ResHisto);
00555     vSubdetXprimeRes_[StripSubdetector::TID - 1]->Add(itTid->second.ResXprimeHisto);
00556     vSubdetNormRes_[StripSubdetector::TID - 1]->Add(itTid->second.NormResHisto);
00557   }
00558   for(std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator itTob = mTobResiduals_.begin(), 
00559         itEnd = mTobResiduals_.end(); itTob != itEnd;++itTob ) {
00560     vSubdetRes_[StripSubdetector::TOB - 1]->Add(itTob->second.ResHisto);
00561     vSubdetXprimeRes_[StripSubdetector::TOB - 1]->Add(itTob->second.ResXprimeHisto);
00562     vSubdetNormRes_[StripSubdetector::TOB - 1]->Add(itTob->second.NormResHisto);
00563   }
00564   for(std::map<int, TrackerOfflineValidation::ModuleHistos>::const_iterator itTec = mTecResiduals_.begin(), 
00565         itEnd = mTecResiduals_.end(); itTec != itEnd;++itTec ) {
00566     vSubdetRes_[StripSubdetector::TEC - 1]->Add(itTec->second.ResHisto);
00567     vSubdetXprimeRes_[StripSubdetector::TEC - 1]->Add(itTec->second.ResXprimeHisto);
00568     vSubdetNormRes_[StripSubdetector::TEC - 1]->Add(itTec->second.NormResHisto);
00569   }
00570   
00571   // create summary histogramms recursively
00572   std::vector<std::pair<TH1*,TH1*> > vTrackerprofiles;
00573   collateSummaryHists((*fs),(aliTracker), 0, aliobjid, vTrackerprofiles);
00574   
00575   
00576 }
00577 
00578 
00579 void
00580 TrackerOfflineValidation::collateSummaryHists( TFileDirectory &tfd, const Alignable& ali, int i, const AlignableObjectId &aliobjid, std::vector<std::pair<TH1*,TH1*> > &v_levelProfiles)
00581 {
00582   
00583   std::vector<Alignable*> alivec(ali.components());
00584 
00585   
00586   if( ((alivec)[0]->alignableObjectId() == align::AlignableDet || 
00587        (alivec)[0]->alignableObjectId() == align::AlignableDetUnit) 
00588      )  {
00589      return;
00590   }
00591 
00592   for(int iComp=0, iCompEnd = ali.components().size();iComp < iCompEnd; ++iComp) {
00593 
00594     std::vector<std::pair<TH1*,TH1*> > v_profiles;        
00595     std::string structurename  = aliobjid.typeToName((alivec)[iComp]->alignableObjectId());
00596  
00597     LogDebug("TrackerOfflineValidation") << "StructureName = " << structurename;
00598     std::stringstream dirname;
00599     
00600     // add no suffix counter to strip and pixel
00601     // just aesthetics
00602     if(structurename != "Strip" && structurename != "Pixel") {
00603       dirname << structurename << "_" << iComp+1;
00604     } else {
00605       dirname << structurename;
00606     }
00607     
00608     if( (structurename != "Det" && structurename != "DetUnit" )  || (alivec)[0]->components().size() > 1
00609        ) {
00610       TFileDirectory f = tfd.mkdir((dirname.str()).c_str());
00611       collateSummaryHists( f, *(alivec)[iComp], i, aliobjid, v_profiles);
00612     
00613       v_levelProfiles.push_back(bookSummaryHists(tfd, *(alivec[iComp]), ali.alignableObjectId(), iComp, aliobjid));
00614       
00615       for(uint n = 0; n < v_profiles.size(); ++n) {
00616         v_levelProfiles[iComp].first->SetBinContent(n+1,v_profiles[n].second->GetMean(1));
00617         v_levelProfiles[iComp].first->SetBinError(n+1,Fwhm(v_profiles[n].second)/2.);
00618         v_levelProfiles[iComp].second->Add(v_profiles[n].second);
00619       }
00620     } else {
00621       // nothing to be done for det or detunits
00622       continue;
00623     }
00624 
00625   }
00626 
00627 }
00628 
00629 std::pair<TH1*,TH1*> 
00630 TrackerOfflineValidation::bookSummaryHists(TFileDirectory &tfd, const Alignable& ali, align::StructureType type, int i, const AlignableObjectId &aliobjid)
00631 {
00632   parameters_ = parset_.getParameter<edm::ParameterSet>("TH1ResModules");
00633   int32_t i_residuals_Nbins =  parameters_.getParameter<int32_t>("Nbinx");
00634   double d_residual_xmin = parameters_.getParameter<double>("xmin");
00635   double d_residual_xmax = parameters_.getParameter<double>("xmax");
00636   uint subsize = ali.components().size();
00637   align::StructureType alitype = ali.alignableObjectId();
00638   align::StructureType subtype = ali.components()[0]->alignableObjectId();
00639   TH1 *h_thissummary = 0;
00640   
00641   if( subtype  != align::AlignableDet || (subtype  == align::AlignableDet && ali.components()[0]->components().size() == 1)
00642       ) {
00643     h_thissummary = tfd.make<TH1F>(Form("h_summary%s_%d",aliobjid.typeToName(alitype).c_str(),i), 
00644                                 Form("Summary for substructures in %s %d;%s;#LT #Delta x #GT",
00645                                      aliobjid.typeToName(alitype).c_str(),i,aliobjid.typeToName(subtype).c_str()),
00646                                 subsize,0.5,subsize+0.5)  ;
00647    
00648   } else if( subtype == align::AlignableDet && subsize > 1) {
00649     h_thissummary = tfd.make<TH1F>(Form("h_summary%s_%d",aliobjid.typeToName(alitype).c_str(),i), 
00650                                 Form("Summary for substructures in %s %d;%s;#LT #Delta x #GT",
00651                                      aliobjid.typeToName(alitype).c_str(),i,aliobjid.typeToName(subtype).c_str()),
00652                                 (2*subsize),0.5,2*subsize+0.5)  ;  
00653   } else {
00654     edm::LogWarning("TrackerOfflineValidation") << "@SUB=TrackerOfflineValidation::bookSummaryHists" 
00655                                                 << "No summary histogramm for hierarchy level" << aliobjid.typeToName(subtype);      
00656   }
00657 
00658   TH1* h_this = tfd.make<TH1F>(Form("h_%s_%d",aliobjid.typeToName(alitype).c_str(),i), 
00659                                 Form("Residual for %s %d in %s ",aliobjid.typeToName(alitype).c_str(),i,
00660                                      aliobjid.typeToName(type).c_str(),aliobjid.typeToName(subtype).c_str()),
00661                                 i_residuals_Nbins,d_residual_xmin,d_residual_xmax);
00662 
00663   
00664   
00665   // special case I: For DetUnits and Detwith  only one subcomponent start filling summary histos
00666   if( (  subtype == align::AlignableDet && ali.components()[0]->components().size() == 1) || 
00667       subtype  == align::AlignableDetUnit  
00668       ) {
00669 
00670     for(uint k=0;k<subsize;++k) {
00671       DetId detid = ali.components()[k]->id();
00672       ModuleHistos &histStruct = this->GetHistStructFromMap(detid);
00673       h_thissummary->SetBinContent(k+1, histStruct.ResHisto->GetMean());
00674       h_thissummary->SetBinError(k+1, histStruct.ResHisto->GetRMS());
00675       h_this->Add(histStruct.ResHisto);
00676     }
00677     
00678   }
00679   // special case II: Fill summary histos for dets with two detunits 
00680   else if( subtype == align::AlignableDet && subsize > 1) {
00681     for(uint k = 0; k < subsize; ++k) { 
00682       uint jEnd = ali.components()[0]->components().size();
00683       for(uint j = 0; j <  jEnd; ++j) {
00684         DetId detid = ali.components()[k]->components()[j]->id();
00685         ModuleHistos &histStruct = this->GetHistStructFromMap(detid);   
00686         h_thissummary->SetBinContent(2*k+j+1,histStruct.ResHisto->GetMean());
00687         h_thissummary->SetBinError(2*k+j+1,histStruct.ResHisto->GetRMS());
00688         h_this->Add( histStruct.ResHisto);
00689 
00690       }
00691     }
00692   }
00693   
00694 
00695   return std::make_pair(h_thissummary,h_this);
00696   
00697   
00698 }
00699 
00700 
00701 float 
00702 TrackerOfflineValidation::Fwhm (const TH1* hist) 
00703 {
00704   float fwhm = 0.;
00705   float max = hist->GetMaximum();
00706   int left = -1, right = -1;
00707   for(unsigned int i = 1, iEnd = hist->GetNbinsX(); i <= iEnd; ++i) {
00708     if(hist->GetBinContent(i) < max/2. && hist->GetBinContent(i+1) > max/2. && left == -1) {
00709       if(max/2. - hist->GetBinContent(i) < hist->GetBinContent(i+1) - max/2.) {
00710         left = i;
00711         ++i;
00712       } else {
00713         left = i+1;
00714         ++i;
00715       }
00716     }
00717     if(left != -1 && right == -1) {
00718       if(hist->GetBinContent(i) > max/2. && hist->GetBinContent(i+1) < max/2.) {
00719         if( hist->GetBinContent(i) - max/2. < max/2. - hist->GetBinContent(i+1)) {
00720           right = i;
00721         } else {
00722           right = i+1;
00723         }
00724         
00725       }
00726     }
00727   }
00728   fwhm = hist->GetXaxis()->GetBinCenter(right) - hist->GetXaxis()->GetBinCenter(left);
00729   return fwhm;
00730 }
00731 
00732 
00733 //define this as a plug-in
00734 DEFINE_FWK_MODULE(TrackerOfflineValidation);