CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10_patch2/src/RecoLocalTracker/SiStripZeroSuppression/plugins/SiStripBaselineAnalyzer.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    SiStripBaselineAnalyzer
00004 // Class:      SiStripBaselineAnalyzer
00005 // 
00013 //
00014 // Original Author:  Ivan Amos Cali
00015 //         Created:  Mon Jul 28 14:10:52 CEST 2008
00016 // $Id: SiStripBaselineAnalyzer.cc,v 1.8 2011/10/05 09:08:20 icali Exp $
00017 //
00018 //
00019  
00020 
00021 // system include files
00022 #include <memory>
00023 #include <iostream>
00024 
00025 // user include files
00026 #include "FWCore/Framework/interface/Frameworkfwd.h"
00027 #include "FWCore/Framework/interface/EDAnalyzer.h"
00028 #include "FWCore/Framework/interface/MakerMacros.h"
00029 
00030 #include "FWCore/Framework/interface/Event.h"
00031 #include "FWCore/Framework/interface/MakerMacros.h"
00032 #include "FWCore/Framework/interface/ESHandle.h"
00033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00034 #include "DataFormats/Common/interface/DetSet.h"
00035 #include "DataFormats/Common/interface/DetSetVector.h"
00036 #include "DataFormats/Common/interface/DetSetVectorNew.h"
00037 
00038 #include "DataFormats/SiStripDigi/interface/SiStripProcessedRawDigi.h"
00039 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
00040 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
00041 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
00042 
00043 #include "CondFormats/SiStripObjects/interface/SiStripPedestals.h"
00044 #include "CondFormats/DataRecord/interface/SiStripPedestalsRcd.h"
00045 
00046 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripPedestalsSubtractor.h"
00047 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripCommonModeNoiseSubtractor.h"
00048 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripRawProcessingFactory.h"
00049 
00050 #include "FWCore/ServiceRegistry/interface/Service.h"
00051 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00052 #include "CommonTools/Utils/interface/TFileDirectory.h"
00053 
00054 //ROOT inclusion
00055 #include "TROOT.h"
00056 #include "TFile.h"
00057 #include "TNtuple.h"
00058 #include "TMath.h"
00059 #include "TCanvas.h"
00060 #include "TH1F.h"
00061 #include "TH2F.h"
00062 #include "TProfile.h"
00063 #include "TList.h"
00064 #include "TString.h"
00065 #include "TStyle.h"
00066 #include "TGraph.h"
00067 #include "TMultiGraph.h"
00068 #include "THStack.h"
00069 
00070 
00071 //
00072 // class decleration
00073 //
00074 
00075 class SiStripBaselineAnalyzer : public edm::EDAnalyzer {
00076    public:
00077       explicit SiStripBaselineAnalyzer(const edm::ParameterSet&);
00078       ~SiStripBaselineAnalyzer();
00079 
00080 
00081    private:
00082       virtual void beginJob() ;
00083       virtual void analyze(const edm::Event&, const edm::EventSetup&);
00084       virtual void endJob() ;
00085       
00086           std::auto_ptr<SiStripPedestalsSubtractor>   subtractorPed_;
00087           edm::ESHandle<SiStripPedestals> pedestalsHandle;
00088           std::vector<int> pedestals;
00089           uint32_t peds_cache_id;
00090 
00091           bool plotClusters_;
00092           bool plotBaseline_;
00093           bool plotBaselinePoints_;
00094           bool plotRawDigi_;
00095           bool plotAPVCM_;
00096           bool plotPedestals_;
00097           
00098           edm::InputTag srcBaseline_;
00099           edm::InputTag srcBaselinePoints_;
00100           edm::InputTag srcAPVCM_;
00101           edm::InputTag srcProcessedRawDigi_;
00102       
00103           edm::Service<TFileService> fs_;
00104   
00105           TH1F* h1BadAPVperEvent_;
00106           
00107           TH1F* h1ProcessedRawDigis_;
00108           TH1F* h1Baseline_;
00109           TH1F* h1Clusters_;
00110           TH1F* h1APVCM_;
00111           TH1F* h1Pedestals_;     
00112           
00113           TCanvas* Canvas_;
00114           std::vector<TH1F> vProcessedRawDigiHisto_;
00115           std::vector<TH1F> vBaselineHisto_;
00116           std::vector<TH1F> vBaselinePointsHisto_;
00117           std::vector<TH1F> vClusterHisto_;
00118           
00119           uint16_t nModuletoDisplay_;
00120           uint16_t actualModule_;
00121 };
00122 
00123 
00124 SiStripBaselineAnalyzer::SiStripBaselineAnalyzer(const edm::ParameterSet& conf){
00125    
00126   srcBaseline_ =  conf.getParameter<edm::InputTag>( "srcBaseline" );
00127   srcBaselinePoints_ = conf.getParameter<edm::InputTag>( "srcBaselinePoints" );
00128   srcProcessedRawDigi_ =  conf.getParameter<edm::InputTag>( "srcProcessedRawDigi" );
00129   srcAPVCM_ =  conf.getParameter<edm::InputTag>( "srcAPVCM" );
00130   subtractorPed_ = SiStripRawProcessingFactory::create_SubtractorPed(conf.getParameter<edm::ParameterSet>("Algorithms"));
00131   nModuletoDisplay_ = conf.getParameter<uint32_t>( "nModuletoDisplay" );
00132   plotClusters_ = conf.getParameter<bool>( "plotClusters" );
00133   plotBaseline_ = conf.getParameter<bool>( "plotBaseline" );
00134   plotBaselinePoints_ = conf.getParameter<bool>( "plotBaselinePoints" );
00135   plotRawDigi_ = conf.getParameter<bool>( "plotRawDigi" );
00136   plotAPVCM_ = conf.getParameter<bool>( "plotAPVCM" );
00137   plotPedestals_ = conf.getParameter<bool>( "plotPedestals" );
00138 
00139   h1BadAPVperEvent_ = fs_->make<TH1F>("BadAPV/Event","BadAPV/Event", 2001, -0.5, 2000.5);
00140   h1BadAPVperEvent_->SetXTitle("# Modules with Bad APVs");
00141   h1BadAPVperEvent_->SetYTitle("Entries");
00142   h1BadAPVperEvent_->SetLineWidth(2);
00143   h1BadAPVperEvent_->SetLineStyle(2);
00144 
00145   h1APVCM_ = fs_->make<TH1F>("APV CM","APV CM", 2048, -1023.5, 1023.5);
00146   h1APVCM_->SetXTitle("APV CM [adc]");
00147   h1APVCM_->SetYTitle("Entries");
00148   h1APVCM_->SetLineWidth(2);
00149   h1APVCM_->SetLineStyle(2);
00150 
00151   h1Pedestals_ = fs_->make<TH1F>("Pedestals","Pedestals", 2048, -1023.5, 1023.5);
00152   h1Pedestals_->SetXTitle("Pedestals [adc]");
00153   h1Pedestals_->SetYTitle("Entries");
00154   h1Pedestals_->SetLineWidth(2);
00155   h1Pedestals_->SetLineStyle(2);
00156 
00157   
00158  
00159 }
00160 
00161 
00162 SiStripBaselineAnalyzer::~SiStripBaselineAnalyzer()
00163 {
00164  
00165    
00166 
00167 }
00168 
00169 void
00170 SiStripBaselineAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& es)
00171 {
00172    using namespace edm;
00173    if(plotPedestals_&&actualModule_ ==0){
00174       uint32_t p_cache_id = es.get<SiStripPedestalsRcd>().cacheIdentifier();
00175       if(p_cache_id != peds_cache_id) {
00176         es.get<SiStripPedestalsRcd>().get(pedestalsHandle);
00177         peds_cache_id = p_cache_id;
00178       }
00179       
00180       
00181       std::vector<uint32_t> detIdV;
00182       pedestalsHandle->getDetIds(detIdV);
00183       
00184       for(uint32_t i=0; i < detIdV.size(); ++i){
00185         pedestals.clear();
00186         SiStripPedestals::Range pedestalsRange = pedestalsHandle->getRange(detIdV[i]);
00187         pedestals.resize((pedestalsRange.second- pedestalsRange.first)*8/10);
00188         pedestalsHandle->allPeds(pedestals, pedestalsRange);
00189         for(uint32_t it=0; it < pedestals.size(); ++it) h1Pedestals_->Fill(pedestals[it]);
00190       }
00191    }
00192 
00193    if(plotAPVCM_){
00194      edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi> > moduleCM;
00195      edm::InputTag CMLabel("siStripZeroSuppression:APVCM");
00196      e.getByLabel(srcAPVCM_,moduleCM);
00197 
00198      edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator itCMDetSetV =moduleCM->begin();
00199      for (; itCMDetSetV != moduleCM->end(); ++itCMDetSetV){  
00200        edm::DetSet<SiStripProcessedRawDigi>::const_iterator  itCM= itCMDetSetV->begin();
00201        for(;itCM != itCMDetSetV->end(); ++itCM) h1APVCM_->Fill(itCM->adc());
00202      }
00203    }
00204    
00205    if(!plotRawDigi_) return;
00206    subtractorPed_->init(es);
00207    
00208    
00209  
00210    edm::Handle< edm::DetSetVector<SiStripRawDigi> > moduleRawDigi;
00211    e.getByLabel(srcProcessedRawDigi_,moduleRawDigi);
00212  
00213    edm::Handle<edm::DetSetVector<SiStripProcessedRawDigi> > moduleBaseline;
00214    if(plotBaseline_) e.getByLabel(srcBaseline_, moduleBaseline); 
00215 
00216    edm::Handle<edm::DetSetVector<SiStripDigi> > moduleBaselinePoints;
00217    if(plotBaselinePoints_) e.getByLabel(srcBaseline_, moduleBaselinePoints); 
00218    
00219    edm::Handle<edmNew::DetSetVector<SiStripCluster> > clusters;
00220    if(plotClusters_){
00221         edm::InputTag clusLabel("siStripClusters");
00222         e.getByLabel(clusLabel, clusters);
00223    }
00224    
00225    char detIds[20];
00226    char evs[20];
00227    char runs[20];    
00228    
00229 
00230    TFileDirectory sdProcessedRawDigis_= fs_->mkdir("ProcessedRawDigis");
00231    TFileDirectory sdBaseline_= fs_->mkdir("Baseline");
00232    TFileDirectory sdBaselinePoints_= fs_->mkdir("BaselinePoints");
00233    TFileDirectory sdClusters_= fs_->mkdir("Clusters");
00234    
00235 
00236    edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator itDSBaseline;
00237    if(plotBaseline_) itDSBaseline = moduleBaseline->begin();
00238    edm::DetSetVector<SiStripRawDigi>::const_iterator itRawDigis = moduleRawDigi->begin();
00239    
00240    uint32_t NBabAPVs = moduleRawDigi->size();     
00241    std::cout<< "Number of module with HIP in this event: " << NBabAPVs << std::endl;
00242    h1BadAPVperEvent_->Fill(NBabAPVs);
00243    
00244    for (; itRawDigis != moduleRawDigi->end(); ++itRawDigis) {
00245       if(actualModule_ > nModuletoDisplay_) return;
00246       uint32_t detId = itRawDigis->id;
00247           
00248       if(plotBaseline_){
00249         //std::cout << "bas id: " << itDSBaseline->id << " raw id: " << detId << std::endl;
00250         if(itDSBaseline->id != detId){
00251                 std::cout << "Collections out of Synch. Something of fishy is going on ;-)" << std::endl;
00252                 return;
00253         }         
00254       }
00255       
00256     
00257       actualModule_++;
00258       uint32_t event = e.id().event();
00259       uint32_t run = e.id().run();
00260       //std::cout << "processing module N: " << actualModule_<< " detId: " << detId << " event: "<< event << std::endl; 
00261          
00262 
00263           
00264       edm::DetSet<SiStripRawDigi>::const_iterator itRaw = itRawDigis->begin(); 
00265       bool restAPV[6] = {0,0,0,0,0,0};
00266       int strip =0, totADC=0;
00267       int minAPVRes = 7, maxAPVRes = -1;
00268       for(;itRaw != itRawDigis->end(); ++itRaw, ++strip){
00269             float adc = itRaw->adc();
00270             totADC+= adc;
00271             if(strip%127 ==0){
00272                 //std::cout << "totADC " << totADC << std::endl;
00273               int APV = strip/128;
00274               if(totADC!= 0){
00275                 restAPV[APV] = true;
00276                         totADC =0;
00277                         if(APV>maxAPVRes) maxAPVRes = APV;
00278                         if(APV<minAPVRes) minAPVRes = APV;
00279               }
00280             }
00281       }
00282 
00283       uint16_t bins =768;
00284       float minx = -0.5, maxx=767.5;
00285       if(minAPVRes !=7){
00286         minx = minAPVRes * 128 -0.5;
00287         maxx = maxAPVRes * 128 + 127.5;
00288         bins = maxx-minx;
00289       }
00290       
00291       sprintf(detIds,"%ul", detId);
00292       sprintf(evs,"%ul", event);
00293       sprintf(runs,"%ul", run);
00294       char* dHistoName = Form("Id:%s_run:%s_ev:%s",detIds, runs, evs);
00295       h1ProcessedRawDigis_ = sdProcessedRawDigis_.make<TH1F>(dHistoName,dHistoName, bins, minx, maxx); 
00296       
00297       if(plotBaseline_){
00298         h1Baseline_ = sdBaseline_.make<TH1F>(dHistoName,dHistoName, bins, minx, maxx); 
00299         h1Baseline_->SetXTitle("strip#");
00300         h1Baseline_->SetYTitle("ADC");
00301         h1Baseline_->SetMaximum(1024.);
00302         h1Baseline_->SetMinimum(-300.);
00303         h1Baseline_->SetLineWidth(2);
00304         h1Baseline_->SetLineStyle(2);
00305         h1Baseline_->SetLineColor(2);
00306       }
00307 
00308       if(plotClusters_){
00309         h1Clusters_ = sdClusters_.make<TH1F>(dHistoName,dHistoName, bins, minx, maxx);
00310           
00311         h1Clusters_->SetXTitle("strip#");
00312         h1Clusters_->SetYTitle("ADC");
00313         h1Clusters_->SetMaximum(1024.);
00314         h1Clusters_->SetMinimum(-300.);
00315         h1Clusters_->SetLineWidth(2);
00316         h1Clusters_->SetLineStyle(2);
00317         h1Clusters_->SetLineColor(3);
00318       }
00319 
00320       h1ProcessedRawDigis_->SetXTitle("strip#");  
00321       h1ProcessedRawDigis_->SetYTitle("ADC");
00322       h1ProcessedRawDigis_->SetMaximum(1024.);
00323       h1ProcessedRawDigis_->SetMinimum(-300.);
00324       h1ProcessedRawDigis_->SetLineWidth(2);
00325       
00326       std::vector<int16_t> ProcessedRawDigis(itRawDigis->size());
00327       subtractorPed_->subtract( *itRawDigis, ProcessedRawDigis);
00328 
00329       edm::DetSet<SiStripProcessedRawDigi>::const_iterator  itBaseline;
00330       if(plotBaseline_) itBaseline = itDSBaseline->begin(); 
00331       std::vector<int16_t>::const_iterator itProcessedRawDigis;
00332                   
00333       strip =0;      
00334       for(itProcessedRawDigis = ProcessedRawDigis.begin();itProcessedRawDigis != ProcessedRawDigis.end(); ++itProcessedRawDigis){
00335         if(restAPV[strip/128]){
00336           float adc = *itProcessedRawDigis;       
00337           h1ProcessedRawDigis_->Fill(strip, adc);
00338           if(plotBaseline_){
00339             h1Baseline_->Fill(strip, itBaseline->adc()); 
00340             ++itBaseline;
00341           }
00342          }
00343         ++strip;
00344        }          
00345        
00346       if(plotBaseline_) ++itDSBaseline; 
00347       if(plotClusters_){
00348           edmNew::DetSetVector<SiStripCluster>::const_iterator itClusters = clusters->begin();
00349           for ( ; itClusters != clusters->end(); ++itClusters ){
00350                 for ( edmNew::DetSet<SiStripCluster>::const_iterator clus =     itClusters->begin(); clus != itClusters->end(); ++clus){
00351                   if(itClusters->id() == detId){
00352                     int firststrip = clus->firstStrip();
00353                     //std::cout << "Found cluster in detId " << detId << " " << firststrip << " " << clus->amplitudes().size() << " -----------------------------------------------" << std::endl;              
00354                     strip=0;
00355                     for( std::vector<uint8_t>::const_iterator itAmpl = clus->amplitudes().begin(); itAmpl != clus->amplitudes().end(); ++itAmpl){
00356                       h1Clusters_->Fill(firststrip+strip, *itAmpl);
00357                       ++strip;
00358                     }
00359                   }              
00360                 }
00361           }
00362       }
00363    }            
00364 
00365 }
00366 
00367 
00368 // ------------ method called once each job just before starting event loop  ------------
00369 void SiStripBaselineAnalyzer::beginJob()
00370 {
00371   
00372   
00373 actualModule_ =0;  
00374    
00375  
00376 }
00377 
00378 // ------------ method called once each job just after ending the event loop  ------------
00379 void 
00380 SiStripBaselineAnalyzer::endJob() {
00381      
00382 }
00383 
00384 //define this as a plug-in
00385 DEFINE_FWK_MODULE(SiStripBaselineAnalyzer);
00386