CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/DQM/SiStripMonitorTrack/src/SiStripMonitorTrack.cc

Go to the documentation of this file.
00001 #include <FWCore/MessageLogger/interface/MessageLogger.h>
00002 
00003 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00004 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00005 
00006 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
00007 #include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
00008 #include "DataFormats/DetId/interface/DetId.h"
00009 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00010 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00011 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
00012 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00013 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00015 #include "CalibTracker/SiStripCommon/interface/SiStripDCSStatus.h"
00016 #include "CommonTools/TriggerUtils/interface/GenericTriggerEventFlag.h"
00017 
00018 #include "DQM/SiStripMonitorTrack/interface/SiStripMonitorTrack.h"
00019 
00020 #include "DQM/SiStripCommon/interface/SiStripHistoId.h"
00021 #include "TMath.h"
00022 
00023 SiStripMonitorTrack::SiStripMonitorTrack(const edm::ParameterSet& conf): 
00024   dbe(edm::Service<DQMStore>().operator->()),
00025   conf_(conf),
00026   tracksCollection_in_EventTree(true),
00027   firstEvent(-1),
00028   genTriggerEventFlag_(new GenericTriggerEventFlag(conf))
00029 {
00030   Cluster_src_   = conf.getParameter<edm::InputTag>("Cluster_src");
00031   Mod_On_        = conf.getParameter<bool>("Mod_On");
00032   Trend_On_      = conf.getParameter<bool>("Trend_On");
00033   flag_ring      = conf.getParameter<bool>("RingFlag_On");
00034   TkHistoMap_On_ = conf.getParameter<bool>("TkHistoMap_On");
00035 
00036   TrackProducer_ = conf_.getParameter<std::string>("TrackProducer");
00037   TrackLabel_ = conf_.getParameter<std::string>("TrackLabel");
00038 
00039   // cluster quality conditions 
00040   edm::ParameterSet cluster_condition = conf_.getParameter<edm::ParameterSet>("ClusterConditions");
00041   applyClusterQuality_ = cluster_condition.getParameter<bool>("On");
00042   sToNLowerLimit_      = cluster_condition.getParameter<double>("minStoN");
00043   sToNUpperLimit_      = cluster_condition.getParameter<double>("maxStoN");
00044   widthLowerLimit_     = cluster_condition.getParameter<double>("minWidth"); 
00045   widthUpperLimit_     = cluster_condition.getParameter<double>("maxWidth"); 
00046 
00047 
00048   // Create DCS Status
00049   bool checkDCS    = conf_.getParameter<bool>("UseDCSFiltering");
00050   if (checkDCS) dcsStatus_ = new SiStripDCSStatus();
00051   else dcsStatus_ = 0; 
00052 }
00053 
00054 //------------------------------------------------------------------------
00055 SiStripMonitorTrack::~SiStripMonitorTrack() { 
00056   if (dcsStatus_) delete dcsStatus_;
00057   if (genTriggerEventFlag_) delete genTriggerEventFlag_;
00058 }
00059 
00060 //------------------------------------------------------------------------
00061 void SiStripMonitorTrack::beginRun(const edm::Run& run,const edm::EventSetup& es)
00062 {
00063   //Retrieve tracker topology from geometry
00064   edm::ESHandle<TrackerTopology> tTopoHandle;
00065   es.get<IdealGeometryRecord>().get(tTopoHandle);
00066   const TrackerTopology* const tTopo = tTopoHandle.product();
00067 
00068   //get geom 
00069   es.get<TrackerDigiGeometryRecord>().get( tkgeom );
00070   LogDebug("SiStripMonitorTrack") << "[SiStripMonitorTrack::beginRun] There are "<<tkgeom->detUnits().size() <<" detectors instantiated in the geometry" << std::endl;  
00071   es.get<SiStripDetCablingRcd>().get( SiStripDetCabling_ );
00072 
00073   book(tTopo);
00074 
00075   // Initialize the GenericTriggerEventFlag
00076   if ( genTriggerEventFlag_->on() )genTriggerEventFlag_->initRun( run, es );
00077 }
00078 
00079 //------------------------------------------------------------------------
00080 void SiStripMonitorTrack::endJob(void)
00081 {
00082   if(conf_.getParameter<bool>("OutputMEsInRootFile")){
00083     dbe->showDirStructure();
00084     dbe->save(conf_.getParameter<std::string>("OutputFileName"));
00085   }
00086 }
00087 
00088 // ------------ method called to produce the data  ------------
00089 void SiStripMonitorTrack::analyze(const edm::Event& e, const edm::EventSetup& es)
00090 {
00091   // Filter out events if DCS checking is requested
00092   if (dcsStatus_ && !dcsStatus_->getStatus(e,es)) return;
00093   
00094   // Filter out events if Trigger Filtering is requested
00095   if (genTriggerEventFlag_->on()&& ! genTriggerEventFlag_->accept( e, es) ) return;
00096   
00097   //initialization of global quantities
00098   LogDebug("SiStripMonitorTrack") << "[SiStripMonitorTrack::analyse]  " << "Run " << e.id().run() << " Event " << e.id().event() << std::endl;
00099   runNb   = e.id().run();
00100   eventNb = e.id().event();
00101   vPSiStripCluster.clear();
00102   
00103   iOrbitSec = e.orbitNumber()/11223.0;
00104 
00105   // initialise # of clusters
00106   for (std::map<std::string, SubDetMEs>::iterator iSubDet = SubDetMEsMap.begin();
00107        iSubDet != SubDetMEsMap.end(); iSubDet++) {
00108     iSubDet->second.totNClustersOnTrack = 0;
00109     iSubDet->second.totNClustersOffTrack = 0;
00110   }  
00111   
00112   //Perform track study
00113   trackStudy(e, es);
00114   
00115   //Perform Cluster Study (irrespectively to tracks)
00116 
00117    AllClusters(e, es); //analyzes the off Track Clusters
00118 
00119   //Summary Counts of clusters
00120   std::map<std::string, MonitorElement*>::iterator iME;
00121   std::map<std::string, LayerMEs>::iterator        iLayerME;
00122  
00123   for (std::map<std::string, SubDetMEs>::iterator iSubDet = SubDetMEsMap.begin();
00124        iSubDet != SubDetMEsMap.end(); iSubDet++) {
00125     SubDetMEs subdet_mes = iSubDet->second;
00126     if (subdet_mes.totNClustersOnTrack > 0) {
00127       fillME(subdet_mes.nClustersOnTrack, subdet_mes.totNClustersOnTrack);
00128     }
00129     fillME(subdet_mes.nClustersOffTrack, subdet_mes.totNClustersOffTrack);
00130     if (Trend_On_) {
00131       fillME(subdet_mes.nClustersTrendOnTrack,iOrbitSec,subdet_mes.totNClustersOnTrack);
00132       fillME(subdet_mes.nClustersTrendOffTrack,iOrbitSec,subdet_mes.totNClustersOffTrack);
00133     }
00134   }  
00135 }
00136 
00137 //------------------------------------------------------------------------  
00138 void SiStripMonitorTrack::book(const TrackerTopology* tTopo)
00139 {
00140   
00141   SiStripFolderOrganizer folder_organizer;
00142   //******** TkHistoMaps
00143   if (TkHistoMap_On_) {
00144     tkhisto_StoNCorrOnTrack = new TkHistoMap("SiStrip/TkHisto" ,"TkHMap_StoNCorrOnTrack",0.0,1); 
00145     tkhisto_NumOnTrack  = new TkHistoMap("SiStrip/TkHisto", "TkHMap_NumberOfOnTrackCluster",0.0,1);
00146     tkhisto_NumOffTrack = new TkHistoMap("SiStrip/TkHisto", "TkHMap_NumberOfOfffTrackCluster",0.0,1);
00147   }
00148   //******** TkHistoMaps
00149 
00150   std::vector<uint32_t> vdetId_;
00151   SiStripDetCabling_->addActiveDetectorsRawIds(vdetId_);
00152   //Histos for each detector, layer and module
00153   for (std::vector<uint32_t>::const_iterator detid_iter=vdetId_.begin();detid_iter!=vdetId_.end();detid_iter++){  //loop on all the active detid
00154     uint32_t detid = *detid_iter;
00155     
00156     if (detid < 1){
00157       edm::LogError("SiStripMonitorTrack")<< "[" <<__PRETTY_FUNCTION__ << "] invalid detid " << detid<< std::endl;
00158       continue;
00159     }
00160 
00161     
00162     std::string name;    
00163     
00164     // book Layer and RING plots
00165     std::pair<std::string,int32_t> det_layer_pair = folder_organizer.GetSubDetAndLayer(detid,tTopo,flag_ring);
00166     
00167     SiStripHistoId hidmanager;
00168     std::string layer_id = hidmanager.getSubdetid(detid, tTopo, flag_ring);
00169     std::map<std::string, LayerMEs>::iterator iLayerME  = LayerMEsMap.find(layer_id);
00170     if(iLayerME==LayerMEsMap.end()){
00171       folder_organizer.setLayerFolder(detid, tTopo, det_layer_pair.second, flag_ring);
00172       bookLayerMEs(detid, layer_id);
00173     }
00174     // book sub-detector plots
00175     std::pair<std::string,std::string> sdet_pair = folder_organizer.getSubDetFolderAndTag(detid, tTopo);
00176     if (SubDetMEsMap.find(sdet_pair.second) == SubDetMEsMap.end()){
00177       dbe->setCurrentFolder(sdet_pair.first);
00178       bookSubDetMEs(sdet_pair.second);        
00179     }
00180     // book module plots
00181     if(Mod_On_) {
00182       folder_organizer.setDetectorFolder(detid,tTopo);
00183       bookModMEs(*detid_iter);
00184     } 
00185   }//end loop on detectors detid
00186 }
00187   
00188 //--------------------------------------------------------------------------------
00189 void SiStripMonitorTrack::bookModMEs(const uint32_t & id)//Histograms at MODULE level
00190 {
00191   std::string name = "det";
00192   SiStripHistoId hidmanager;
00193   std::string hid = hidmanager.createHistoId("",name,id);
00194   std::map<std::string, ModMEs>::iterator iModME  = ModMEsMap.find(hid);
00195   if(iModME==ModMEsMap.end()){
00196     ModMEs theModMEs; 
00197     theModMEs.ClusterStoNCorr   = 0;
00198     theModMEs.ClusterCharge     = 0;
00199     theModMEs.ClusterChargeCorr = 0;
00200     theModMEs.ClusterWidth      = 0;
00201     theModMEs.ClusterPos        = 0;
00202     theModMEs.ClusterPGV        = 0;
00203 
00204     // Cluster Width
00205     theModMEs.ClusterWidth=bookME1D("TH1ClusterWidth", hidmanager.createHistoId("ClusterWidth_OnTrack",name,id).c_str()); 
00206     dbe->tag(theModMEs.ClusterWidth,id); 
00207     // Cluster Charge
00208     theModMEs.ClusterCharge=bookME1D("TH1ClusterCharge", hidmanager.createHistoId("ClusterCharge_OnTrack",name,id).c_str());
00209     dbe->tag(theModMEs.ClusterCharge,id); 
00210     // Cluster Charge Corrected
00211     theModMEs.ClusterChargeCorr=bookME1D("TH1ClusterChargeCorr", hidmanager.createHistoId("ClusterChargeCorr_OnTrack",name,id).c_str());
00212     dbe->tag(theModMEs.ClusterChargeCorr,id); 
00213     // Cluster StoN Corrected
00214     theModMEs.ClusterStoNCorr=bookME1D("TH1ClusterStoNCorrMod", hidmanager.createHistoId("ClusterStoNCorr_OnTrack",name,id).c_str());
00215     dbe->tag(theModMEs.ClusterStoNCorr,id); 
00216     // Cluster Position
00217     short total_nr_strips = SiStripDetCabling_->nApvPairs(id) * 2 * 128;
00218     theModMEs.ClusterPos=dbe->book1D(hidmanager.createHistoId("ClusterPosition_OnTrack",name,id).c_str(),hidmanager.createHistoId("ClusterPosition_OnTrack",name,id).c_str(),total_nr_strips,0.5,total_nr_strips+0.5);
00219     dbe->tag(theModMEs.ClusterPos,id); 
00220     // Cluster PGV
00221     theModMEs.ClusterPGV=bookMEProfile("TProfileClusterPGV", hidmanager.createHistoId("PGV_OnTrack",name,id).c_str()); 
00222     dbe->tag(theModMEs.ClusterPGV,id); 
00223 
00224     ModMEsMap[hid]=theModMEs;
00225   }
00226 }
00227 //
00228 // -- Book Layer Level Histograms and Trend plots
00229 //
00230 void SiStripMonitorTrack::bookLayerMEs(const uint32_t& mod_id, std::string& layer_id)
00231 {
00232   std::string name = "layer";  
00233   std::string hname;
00234   SiStripHistoId hidmanager;
00235 
00236   LayerMEs theLayerMEs; 
00237   theLayerMEs.ClusterStoNCorrOnTrack   = 0;
00238   theLayerMEs.ClusterChargeCorrOnTrack = 0;
00239   theLayerMEs.ClusterChargeOnTrack     = 0;
00240   theLayerMEs.ClusterChargeOffTrack    = 0;
00241   theLayerMEs.ClusterNoiseOnTrack      = 0;
00242   theLayerMEs.ClusterNoiseOffTrack     = 0;
00243   theLayerMEs.ClusterWidthOnTrack      = 0;
00244   theLayerMEs.ClusterWidthOffTrack     = 0;
00245   theLayerMEs.ClusterPosOnTrack        = 0;
00246   theLayerMEs.ClusterPosOffTrack       = 0;
00247   
00248   // Cluster StoN Corrected
00249   hname = hidmanager.createHistoLayer("Summary_ClusterStoNCorr",name,layer_id,"OnTrack");
00250   theLayerMEs.ClusterStoNCorrOnTrack = bookME1D("TH1ClusterStoNCorr", hname.c_str());
00251 
00252   // Cluster Charge Corrected
00253   hname = hidmanager.createHistoLayer("Summary_ClusterChargeCorr",name,layer_id,"OnTrack");
00254   theLayerMEs.ClusterChargeCorrOnTrack = bookME1D("TH1ClusterChargeCorr", hname.c_str());
00255 
00256   // Cluster Charge (On and Off Track)
00257   hname = hidmanager.createHistoLayer("Summary_ClusterCharge",name,layer_id,"OnTrack");
00258   theLayerMEs.ClusterChargeOnTrack = bookME1D("TH1ClusterCharge", hname.c_str());
00259 
00260   hname = hidmanager.createHistoLayer("Summary_ClusterCharge",name,layer_id,"OffTrack");
00261   theLayerMEs.ClusterChargeOffTrack = bookME1D("TH1ClusterCharge", hname.c_str());
00262 
00263   // Cluster Noise (On and Off Track)
00264   hname = hidmanager.createHistoLayer("Summary_ClusterNoise",name,layer_id,"OnTrack");
00265   theLayerMEs.ClusterNoiseOnTrack = bookME1D("TH1ClusterNoise", hname.c_str()); 
00266 
00267   hname = hidmanager.createHistoLayer("Summary_ClusterNoise",name,layer_id,"OffTrack");
00268   theLayerMEs.ClusterNoiseOffTrack = bookME1D("TH1ClusterNoise", hname.c_str()); 
00269 
00270   // Cluster Width (On and Off Track)
00271   hname = hidmanager.createHistoLayer("Summary_ClusterWidth",name,layer_id,"OnTrack");
00272   theLayerMEs.ClusterWidthOnTrack = bookME1D("TH1ClusterWidth", hname.c_str()); 
00273 
00274   hname = hidmanager.createHistoLayer("Summary_ClusterWidth",name,layer_id,"OffTrack");
00275   theLayerMEs.ClusterWidthOffTrack = bookME1D("TH1ClusterWidth", hname.c_str()); 
00276 
00277   //Cluster Position
00278   short total_nr_strips = SiStripDetCabling_->nApvPairs(mod_id) * 2 * 128; 
00279   if (layer_id.find("TEC") != std::string::npos && !flag_ring)  total_nr_strips = 3 * 2 * 128;
00280   
00281   hname = hidmanager.createHistoLayer("Summary_ClusterPosition",name,layer_id,"OnTrack");
00282   theLayerMEs.ClusterPosOnTrack = dbe->book1D(hname, hname, total_nr_strips, 0.5,total_nr_strips+0.5);
00283   
00284   hname = hidmanager.createHistoLayer("Summary_ClusterPosition",name,layer_id,"OffTrack");
00285   theLayerMEs.ClusterPosOffTrack = dbe->book1D(hname, hname, total_nr_strips, 0.5,total_nr_strips+0.5);
00286   
00287   //bookeeping
00288   LayerMEsMap[layer_id]=theLayerMEs;
00289 }
00290 //
00291 // -- Book Histograms at Sub-Detector Level
00292 //
00293 void SiStripMonitorTrack::bookSubDetMEs(std::string& name){
00294 
00295   std::string subdet_tag;
00296   subdet_tag = "__" + name;
00297   std::string completeName;
00298 
00299   SubDetMEs theSubDetMEs;
00300   theSubDetMEs.totNClustersOnTrack    = 0;
00301   theSubDetMEs.totNClustersOffTrack   = 0;
00302   theSubDetMEs.nClustersOnTrack       = 0;
00303   theSubDetMEs.nClustersTrendOnTrack  = 0;
00304   theSubDetMEs.nClustersOffTrack      = 0;
00305   theSubDetMEs.nClustersTrendOffTrack = 0;
00306   theSubDetMEs.ClusterStoNCorrOnTrack = 0;
00307   theSubDetMEs.ClusterChargeOffTrack  = 0;
00308   theSubDetMEs.ClusterStoNOffTrack    = 0;
00309 
00310   // TotalNumber of Cluster OnTrack
00311   completeName = "Summary_TotalNumberOfClusters_OnTrack" + subdet_tag;
00312   theSubDetMEs.nClustersOnTrack = bookME1D("TH1nClustersOn", completeName.c_str());
00313   theSubDetMEs.nClustersOnTrack->getTH1()->StatOverflows(kTRUE);
00314   
00315   // TotalNumber of Cluster OffTrack
00316   completeName = "Summary_TotalNumberOfClusters_OffTrack" + subdet_tag;
00317   theSubDetMEs.nClustersOffTrack = bookME1D("TH1nClustersOff", completeName.c_str());
00318   theSubDetMEs.nClustersOffTrack->getTH1()->StatOverflows(kTRUE);
00319   
00320   // Cluster StoN On Track
00321   completeName = "Summary_ClusterStoNCorr_OnTrack"  + subdet_tag;
00322   theSubDetMEs.ClusterStoNCorrOnTrack = bookME1D("TH1ClusterStoNCorr", completeName.c_str());
00323   
00324   // Cluster Charge Off Track
00325   completeName = "Summary_ClusterCharge_OffTrack" + subdet_tag;
00326   theSubDetMEs.ClusterChargeOffTrack=bookME1D("TH1ClusterCharge", completeName.c_str());
00327   
00328   // Cluster Charge StoN Off Track
00329   completeName = "Summary_ClusterStoN_OffTrack"  + subdet_tag;
00330   theSubDetMEs.ClusterStoNOffTrack = bookME1D("TH1ClusterStoN", completeName.c_str());
00331   
00332   if(Trend_On_){
00333     // TotalNumber of Cluster 
00334     completeName = "Trend_TotalNumberOfClusters_OnTrack"  + subdet_tag;
00335     theSubDetMEs.nClustersTrendOnTrack = bookMETrend("TH1nClustersOn", completeName.c_str());
00336     completeName = "Trend_TotalNumberOfClusters_OffTrack"  + subdet_tag;
00337     theSubDetMEs.nClustersTrendOffTrack = bookMETrend("TH1nClustersOff", completeName.c_str());
00338   }
00339   //bookeeping
00340   SubDetMEsMap[name]=theSubDetMEs;
00341 }
00342 //--------------------------------------------------------------------------------
00343 
00344 MonitorElement* SiStripMonitorTrack::bookME1D(const char* ParameterSetLabel, const char* HistoName)
00345 {
00346   Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
00347   return dbe->book1D(HistoName,HistoName,
00348                          Parameters.getParameter<int32_t>("Nbinx"),
00349                          Parameters.getParameter<double>("xmin"),
00350                          Parameters.getParameter<double>("xmax")
00351                     );
00352 }
00353 
00354 //--------------------------------------------------------------------------------
00355 MonitorElement* SiStripMonitorTrack::bookME2D(const char* ParameterSetLabel, const char* HistoName)
00356 {
00357   Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
00358   return dbe->book2D(HistoName,HistoName,
00359                      Parameters.getParameter<int32_t>("Nbinx"),
00360                      Parameters.getParameter<double>("xmin"),
00361                      Parameters.getParameter<double>("xmax"),
00362                      Parameters.getParameter<int32_t>("Nbiny"),
00363                      Parameters.getParameter<double>("ymin"),
00364                      Parameters.getParameter<double>("ymax")
00365                      );
00366 }
00367 
00368 //--------------------------------------------------------------------------------
00369 MonitorElement* SiStripMonitorTrack::bookME3D(const char* ParameterSetLabel, const char* HistoName)
00370 {
00371   Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
00372   return dbe->book3D(HistoName,HistoName,
00373                      Parameters.getParameter<int32_t>("Nbinx"),
00374                      Parameters.getParameter<double>("xmin"),
00375                      Parameters.getParameter<double>("xmax"),
00376                      Parameters.getParameter<int32_t>("Nbiny"),
00377                      Parameters.getParameter<double>("ymin"),
00378                      Parameters.getParameter<double>("ymax"),
00379                      Parameters.getParameter<int32_t>("Nbinz"),
00380                      Parameters.getParameter<double>("zmin"),
00381                      Parameters.getParameter<double>("zmax")
00382                      );
00383 }
00384 
00385 //--------------------------------------------------------------------------------
00386 MonitorElement* SiStripMonitorTrack::bookMEProfile(const char* ParameterSetLabel, const char* HistoName)
00387 {
00388     Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
00389     return dbe->bookProfile(HistoName,HistoName,
00390                             Parameters.getParameter<int32_t>("Nbinx"),
00391                             Parameters.getParameter<double>("xmin"),
00392                             Parameters.getParameter<double>("xmax"),
00393                             Parameters.getParameter<int32_t>("Nbiny"),
00394                             Parameters.getParameter<double>("ymin"),
00395                             Parameters.getParameter<double>("ymax"),
00396                             "" );
00397 }
00398 
00399 //--------------------------------------------------------------------------------
00400 MonitorElement* SiStripMonitorTrack::bookMETrend(const char* ParameterSetLabel, const char* HistoName)
00401 {
00402   Parameters =  conf_.getParameter<edm::ParameterSet>(ParameterSetLabel);
00403   edm::ParameterSet ParametersTrend =  conf_.getParameter<edm::ParameterSet>("Trending");
00404   MonitorElement* me = dbe->bookProfile(HistoName,HistoName,
00405                                         ParametersTrend.getParameter<int32_t>("Nbins"),
00406                                         0,
00407                                         ParametersTrend.getParameter<int32_t>("Nbins"),
00408                                         100, //that parameter should not be there !?
00409                                         Parameters.getParameter<double>("xmin"),
00410                                         Parameters.getParameter<double>("xmax"),
00411                                         "" );
00412   if (me->kind() == MonitorElement::DQM_KIND_TPROFILE) me->getTH1()->SetBit(TH1::kCanRebin);
00413 
00414   if(!me) return me;
00415   me->setAxisTitle("Event Time in Seconds",1);
00416   return me;
00417 }
00418 
00419 //------------------------------------------------------------------------------------------
00420  void SiStripMonitorTrack::trackStudy(const edm::Event& ev, const edm::EventSetup& es){
00421 
00422   // track input  
00423  
00424   edm::Handle<reco::TrackCollection > trackCollectionHandle;
00425   ev.getByLabel(TrackProducer_, TrackLabel_, trackCollectionHandle);//takes the track collection
00426   if (!trackCollectionHandle.isValid()){
00427     edm::LogError("SiStripMonitorTrack")<<" Track Collection is not valid !! " << TrackLabel_<<std::endl;
00428     return;
00429   }
00430   
00431   // trajectory input
00432   edm::Handle<TrajTrackAssociationCollection> TItkAssociatorCollection;
00433   ev.getByLabel(TrackProducer_, TrackLabel_, TItkAssociatorCollection);
00434   if( !TItkAssociatorCollection.isValid()){
00435     edm::LogError("SiStripMonitorTrack")<<"Association not found "<<std::endl;
00436     return;
00437   }
00438   
00439   //Perform track study
00440   int i=0;
00441   for(TrajTrackAssociationCollection::const_iterator it =  TItkAssociatorCollection->begin();it !=  TItkAssociatorCollection->end(); ++it){
00442     const edm::Ref<std::vector<Trajectory> > traj_iterator = it->key;  
00443     // Trajectory Map, extract Trajectory for this track
00444     reco::TrackRef trackref = it->val;
00445     LogDebug("SiStripMonitorTrack")
00446       << "Track number "<< i+1 
00447       << "\n\tmomentum: " << trackref->momentum()
00448       << "\n\tPT: " << trackref->pt()
00449       << "\n\tvertex: " << trackref->vertex()
00450       << "\n\timpact parameter: " << trackref->d0()
00451       << "\n\tcharge: " << trackref->charge()
00452       << "\n\tnormalizedChi2: " << trackref->normalizedChi2() 
00453       <<"\n\tFrom EXTRA : "
00454       <<"\n\t\touter PT "<< trackref->outerPt()<<std::endl;
00455     i++;
00456 
00457     const std::vector<TrajectoryMeasurement> & measurements = traj_iterator->measurements();
00458     std::vector<TrajectoryMeasurement>::const_iterator traj_mes_iterator;
00459     int nhit=0;
00460     for(std::vector<TrajectoryMeasurement>::const_iterator traj_mes_iterator= measurements.begin();traj_mes_iterator!=measurements.end();traj_mes_iterator++){//loop on measurements
00461       //trajectory local direction and position on detector
00462       LocalPoint  stateposition;
00463       LocalVector statedirection;
00464       
00465       TrajectoryStateOnSurface  updatedtsos=traj_mes_iterator->updatedState();
00466       ConstRecHitPointer ttrh=traj_mes_iterator->recHit();
00467       if (!ttrh->isValid()) {continue;}
00468       
00469       nhit++;
00470       
00471       const ProjectedSiStripRecHit2D* phit     = dynamic_cast<const ProjectedSiStripRecHit2D*>( ttrh->hit() );
00472       const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>( ttrh->hit() );
00473       const SiStripRecHit2D* hit2D             = dynamic_cast<const SiStripRecHit2D*>( ttrh->hit() );   
00474       const SiStripRecHit1D* hit1D             = dynamic_cast<const SiStripRecHit1D*>( ttrh->hit() );   
00475       
00476       //      RecHitType type=Single;
00477 
00478       if(matchedhit){
00479         LogTrace("SiStripMonitorTrack")<<"\nMatched recHit found"<< std::endl;
00480         //      type=Matched;
00481         
00482         GluedGeomDet * gdet=(GluedGeomDet *)tkgeom->idToDet(matchedhit->geographicalId());
00483         GlobalVector gtrkdirup=gdet->toGlobal(updatedtsos.localMomentum());         
00484         //mono side
00485         const GeomDetUnit * monodet=gdet->monoDet();
00486         statedirection=monodet->toLocal(gtrkdirup);
00487         SiStripRecHit2D m = matchedhit->monoHit();
00488         if(statedirection.mag() != 0)     RecHitInfo<SiStripRecHit2D>(&m,statedirection,trackref,es);
00489         //stereo side
00490         const GeomDetUnit * stereodet=gdet->stereoDet();
00491         statedirection=stereodet->toLocal(gtrkdirup);
00492         SiStripRecHit2D s = matchedhit->stereoHit();
00493         if(statedirection.mag() != 0)     RecHitInfo<SiStripRecHit2D>(&s,statedirection,trackref,es);
00494       }
00495       else if(phit){
00496         LogTrace("SiStripMonitorTrack")<<"\nProjected recHit found"<< std::endl;
00497         //      type=Projected;
00498         GluedGeomDet * gdet=(GluedGeomDet *)tkgeom->idToDet(phit->geographicalId());
00499         
00500         GlobalVector gtrkdirup=gdet->toGlobal(updatedtsos.localMomentum());
00501         const SiStripRecHit2D&  originalhit=phit->originalHit();
00502         const GeomDetUnit * det;
00503         if(!StripSubdetector(originalhit.geographicalId().rawId()).stereo()){
00504           //mono side
00505           LogTrace("SiStripMonitorTrack")<<"\nProjected recHit found  MONO"<< std::endl;
00506           det=gdet->monoDet();
00507           statedirection=det->toLocal(gtrkdirup);
00508           if(statedirection.mag() != 0) RecHitInfo<SiStripRecHit2D>(&(phit->originalHit()),statedirection,trackref,es);
00509         }
00510         else{
00511           LogTrace("SiStripMonitorTrack")<<"\nProjected recHit found STEREO"<< std::endl;
00512           //stereo side
00513           det=gdet->stereoDet();
00514           statedirection=det->toLocal(gtrkdirup);
00515           if(statedirection.mag() != 0) RecHitInfo<SiStripRecHit2D>(&(phit->originalHit()),statedirection,trackref,es);
00516         }
00517       }else if (hit2D){
00518         statedirection=updatedtsos.localMomentum();
00519         if(statedirection.mag() != 0) RecHitInfo<SiStripRecHit2D>(hit2D,statedirection,trackref,es);
00520       } else if (hit1D) {
00521         statedirection=updatedtsos.localMomentum();
00522         if(statedirection.mag() != 0) RecHitInfo<SiStripRecHit1D>(hit1D,statedirection,trackref,es);
00523       } else {
00524         LogDebug ("SiStrioMonitorTrack") 
00525                  << " LocalMomentum: "<<statedirection
00526                  << "\nLocal x-z plane angle: "<<atan2(statedirection.x(),statedirection.z());
00527       }
00528     }
00529   }
00530 }
00531 
00532 template <class T> void SiStripMonitorTrack::RecHitInfo(const T* tkrecHit, LocalVector LV,reco::TrackRef track_ref, const edm::EventSetup& es){
00533     
00534     if(!tkrecHit->isValid()){
00535       LogTrace("SiStripMonitorTrack") <<"\t\t Invalid Hit " << std::endl;
00536       return;  
00537     }
00538     
00539     const uint32_t& detid = tkrecHit->geographicalId().rawId();
00540     if (find(ModulesToBeExcluded_.begin(),ModulesToBeExcluded_.end(),detid)!=ModulesToBeExcluded_.end()){
00541       LogTrace("SiStripMonitorTrack") << "Modules Excluded" << std::endl;
00542       return;
00543     }
00544     
00545     LogTrace("SiStripMonitorTrack")
00546       <<"\n\t\tRecHit on det "<<tkrecHit->geographicalId().rawId()
00547       <<"\n\t\tRecHit in LP "<<tkrecHit->localPosition()
00548       <<"\n\t\tRecHit in GP "<<tkgeom->idToDet(tkrecHit->geographicalId())->surface().toGlobal(tkrecHit->localPosition()) 
00549       <<"\n\t\tRecHit trackLocal vector "<<LV.x() << " " << LV.y() << " " << LV.z() <<std::endl; 
00550 
00551 
00552     //Retrieve tracker topology from geometry
00553     edm::ESHandle<TrackerTopology> tTopoHandle;
00554     es.get<IdealGeometryRecord>().get(tTopoHandle);
00555     const TrackerTopology* const tTopo = tTopoHandle.product();
00556     
00557     //Get SiStripCluster from SiStripRecHit
00558     if ( tkrecHit != NULL ){
00559       const SiStripCluster* SiStripCluster_ = &*(tkrecHit->cluster());
00560       SiStripClusterInfo SiStripClusterInfo_(*SiStripCluster_,es,detid);
00561             
00562       if ( clusterInfos(&SiStripClusterInfo_,detid, tTopo, OnTrack, LV ) ) {
00563         vPSiStripCluster.push_back(SiStripCluster_);
00564       }
00565     }else{
00566      edm::LogError("SiStripMonitorTrack") << "NULL hit" << std::endl;
00567     }     
00568   }
00569 
00570 //------------------------------------------------------------------------
00571 
00572 void SiStripMonitorTrack::AllClusters(const edm::Event& ev, const edm::EventSetup& es) 
00573 {
00574 
00575   //Retrieve tracker topology from geometry
00576   edm::ESHandle<TrackerTopology> tTopoHandle;
00577   es.get<IdealGeometryRecord>().get(tTopoHandle);
00578   const TrackerTopology* const tTopo = tTopoHandle.product();
00579     
00580   edm::Handle< edmNew::DetSetVector<SiStripCluster> > siStripClusterHandle;
00581   ev.getByLabel( Cluster_src_, siStripClusterHandle); 
00582   if (!siStripClusterHandle.isValid()){
00583     edm::LogError("SiStripMonitorTrack")<< "ClusterCollection is not valid!!" << std::endl;
00584     return;
00585   }
00586   //Loop on Dets
00587   for ( edmNew::DetSetVector<SiStripCluster>::const_iterator DSViter=siStripClusterHandle->begin(); DSViter!=siStripClusterHandle->end();DSViter++){
00588     uint32_t detid=DSViter->id();
00589     if (find(ModulesToBeExcluded_.begin(),ModulesToBeExcluded_.end(),detid)!=ModulesToBeExcluded_.end()) continue;
00590     //Loop on Clusters
00591     LogDebug("SiStripMonitorTrack") << "on detid "<< detid << " N Cluster= " << DSViter->size();
00592     edmNew::DetSet<SiStripCluster>::const_iterator ClusIter = DSViter->begin();
00593     for(; ClusIter!=DSViter->end(); ClusIter++) {
00594       if (std::find(vPSiStripCluster.begin(),vPSiStripCluster.end(),&*ClusIter) == vPSiStripCluster.end()){
00595         SiStripClusterInfo SiStripClusterInfo_(*ClusIter,es,detid);
00596         clusterInfos(&SiStripClusterInfo_,detid,tTopo,OffTrack,LV);
00597       }
00598     }
00599   }
00600 }
00601 
00602 //------------------------------------------------------------------------
00603 bool SiStripMonitorTrack::clusterInfos(SiStripClusterInfo* cluster, const uint32_t& detid, const TrackerTopology* tTopo, enum ClusterFlags flag, const LocalVector LV)
00604 {
00605   if (cluster==0) return false;
00606   // if one imposes a cut on the clusters, apply it
00607   if( (applyClusterQuality_) &&
00608       (cluster->signalOverNoise() < sToNLowerLimit_ ||
00609        cluster->signalOverNoise() > sToNUpperLimit_ ||
00610        cluster->width() < widthLowerLimit_ ||
00611        cluster->width() > widthUpperLimit_) ) return false;
00612   // start of the analysis
00613   
00614   std::pair<std::string,std::string> sdet_pair = folderOrganizer_.getSubDetFolderAndTag(detid,tTopo);
00615   std::map<std::string, SubDetMEs>::iterator iSubdet  = SubDetMEsMap.find(sdet_pair.second);
00616   if(iSubdet != SubDetMEsMap.end()){ 
00617     if (flag == OnTrack) iSubdet->second.totNClustersOnTrack++;
00618     else if (flag == OffTrack) iSubdet->second.totNClustersOffTrack++;
00619   }
00620   
00621   float cosRZ = -2;
00622   LogDebug("SiStripMonitorTrack")<< "\n\tLV " << LV.x() << " " << LV.y() << " " << LV.z() << " " << LV.mag() << std::endl;
00623   if (LV.mag()!=0){
00624     cosRZ= fabs(LV.z())/LV.mag();
00625     LogDebug("SiStripMonitorTrack")<< "\n\t cosRZ " << cosRZ << std::endl;
00626   }
00627   std::string name;
00628   
00629   // Filling SubDet/Layer Plots (on Track + off Track)
00630   fillMEs(cluster,detid,tTopo,cosRZ,flag);
00631   
00632   
00633   //******** TkHistoMaps
00634   if (TkHistoMap_On_) {
00635     uint32_t adet=cluster->detId();
00636     float noise = cluster->noiseRescaledByGain();
00637     if(flag==OnTrack){
00638       tkhisto_NumOnTrack->add(adet,1.);
00639       if(noise > 0.0) tkhisto_StoNCorrOnTrack->fill(adet,cluster->signalOverNoise()*cosRZ);
00640       if(noise == 0.0) 
00641         LogDebug("SiStripMonitorTrack") << "Module " << detid << " in Event " << eventNb << " noise " << noise << std::endl;
00642     }
00643     else if(flag==OffTrack){
00644       tkhisto_NumOffTrack->add(adet,1.);
00645       if(cluster->charge() > 250){
00646         LogDebug("SiStripMonitorTrack") << "Module firing " << detid << " in Event " << eventNb << std::endl;
00647       }
00648     }
00649   }
00650 
00651   // Module plots filled only for onTrack Clusters
00652   if(Mod_On_){
00653     if(flag==OnTrack){
00654       SiStripHistoId hidmanager2;
00655       name =hidmanager2.createHistoId("","det",detid);
00656       fillModMEs(cluster,name,cosRZ); 
00657     }
00658   }
00659   return true;
00660 }
00661 
00662 //--------------------------------------------------------------------------------
00663 void SiStripMonitorTrack::fillModMEs(SiStripClusterInfo* cluster,std::string name,float cos)
00664 {
00665   std::map<std::string, ModMEs>::iterator iModME  = ModMEsMap.find(name);
00666   if(iModME!=ModMEsMap.end()){
00667 
00668     float    StoN     = cluster->signalOverNoise();
00669     uint16_t charge   = cluster->charge();
00670     uint16_t width    = cluster->width();
00671     float    position = cluster->baryStrip(); 
00672 
00673     float noise = cluster->noiseRescaledByGain();
00674     if(noise > 0.0) fillME(iModME->second.ClusterStoNCorr ,StoN*cos);
00675     if(noise == 0.0) LogDebug("SiStripMonitorTrack") << "Module " << name << " in Event " << eventNb << " noise " << noise << std::endl;
00676     fillME(iModME->second.ClusterCharge,charge);
00677 
00678     fillME(iModME->second.ClusterChargeCorr,charge*cos);
00679 
00680     fillME(iModME->second.ClusterWidth ,width);
00681     fillME(iModME->second.ClusterPos   ,position);
00682     
00683     //fill the PGV histo
00684     float PGVmax = cluster->maxCharge();
00685     int PGVposCounter = cluster->maxIndex();
00686     for (int i= int(conf_.getParameter<edm::ParameterSet>("TProfileClusterPGV").getParameter<double>("xmin"));i<PGVposCounter;++i)
00687       fillME(iModME->second.ClusterPGV, i,0.);
00688     for (std::vector<uint8_t>::const_iterator it=cluster->stripCharges().begin();it<cluster->stripCharges().end();++it) {
00689       fillME(iModME->second.ClusterPGV, PGVposCounter++,(*it)/PGVmax);
00690     }
00691     for (int i= PGVposCounter;i<int(conf_.getParameter<edm::ParameterSet>("TProfileClusterPGV").getParameter<double>("xmax"));++i)
00692       fillME(iModME->second.ClusterPGV, i,0.);
00693     //end fill the PGV histo
00694   }
00695 }
00696 
00697 //------------------------------------------------------------------------
00698 void SiStripMonitorTrack::fillMEs(SiStripClusterInfo* cluster,uint32_t detid, const TrackerTopology* tTopo, float cos, enum ClusterFlags flag)
00699 { 
00700   std::pair<std::string,int32_t> SubDetAndLayer = folderOrganizer_.GetSubDetAndLayer(detid,tTopo,flag_ring);
00701   SiStripHistoId hidmanager1;
00702   std::string layer_id = hidmanager1.getSubdetid(detid,tTopo,flag_ring); 
00703   
00704   std::pair<std::string,std::string> sdet_pair = folderOrganizer_.getSubDetFolderAndTag(detid,tTopo);
00705   float    StoN     = cluster->signalOverNoise();
00706   float    noise    = cluster->noiseRescaledByGain();
00707   uint16_t charge   = cluster->charge();
00708   uint16_t width    = cluster->width();
00709   float    position = cluster->baryStrip(); 
00710    
00711   std::map<std::string, LayerMEs>::iterator iLayer  = LayerMEsMap.find(layer_id);
00712   if (iLayer != LayerMEsMap.end()) {
00713     if(flag==OnTrack){
00714       if(noise > 0.0) fillME(iLayer->second.ClusterStoNCorrOnTrack, StoN*cos);
00715       if(noise == 0.0) LogDebug("SiStripMonitorTrack") << "Module " << detid << " in Event " << eventNb << " noise " << cluster->noiseRescaledByGain() << std::endl;
00716       fillME(iLayer->second.ClusterChargeCorrOnTrack, charge*cos);
00717       fillME(iLayer->second.ClusterChargeOnTrack, charge);
00718       fillME(iLayer->second.ClusterNoiseOnTrack, noise);
00719       fillME(iLayer->second.ClusterWidthOnTrack, width);
00720       fillME(iLayer->second.ClusterPosOnTrack, position);
00721     } else {
00722       fillME(iLayer->second.ClusterChargeOffTrack, charge);
00723       fillME(iLayer->second.ClusterNoiseOffTrack, noise);
00724       fillME(iLayer->second.ClusterWidthOffTrack, width);
00725       fillME(iLayer->second.ClusterPosOffTrack, position);
00726     }
00727   }
00728   std::map<std::string, SubDetMEs>::iterator iSubdet  = SubDetMEsMap.find(sdet_pair.second);
00729   if(iSubdet != SubDetMEsMap.end() ){
00730     if(flag==OnTrack){
00731       if(noise > 0.0) fillME(iSubdet->second.ClusterStoNCorrOnTrack,StoN*cos);
00732     } else {
00733       fillME(iSubdet->second.ClusterChargeOffTrack,charge);
00734       if(noise > 0.0) fillME(iSubdet->second.ClusterStoNOffTrack,StoN);
00735     }
00736   }
00737 }
00738 //
00739 // -- Get Subdetector Tag from the Folder name
00740 //
00741 void SiStripMonitorTrack::getSubDetTag(std::string& folder_name, std::string& tag){
00742 
00743   tag =  folder_name.substr(folder_name.find("MechanicalView")+15);
00744   if (tag.find("side_") != std::string::npos) {
00745     tag.replace(tag.find_last_of("/"),1,"_");
00746   }
00747 }