CMSSW Reference Manual

00001 // -*- C++ -*-
00002 // Package:    SiStripChannelGain
00003 // Class:      SiStripGainCosmicCalculator
00004 // Original Author:  G. Bruno, D. Kcira
00005 //         Created:  Mon May 20 10:04:31 CET 2007
00006 // $Id: SiStripGainCosmicCalculator.cc,v 1.8 2008/04/23 08:47:10 querten Exp $
00007 #include "CalibTracker/SiStripChannelGain/plugins/SiStripGainCosmicCalculator.h"
00008 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
00009 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "FWCore/Framework/interface/ESHandle.h"
00012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00013 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00014 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
00015 #include "Geometry/CommonTopologies/interface/StripTopology.h"
00016 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
00017 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
00018 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
00019 #include "CLHEP/Random/RandFlat.h"
00020 #include "CLHEP/Random/RandGauss.h"
00021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00023 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h"
00024 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
00025 #include "DataFormats/TrackReco/interface/Track.h"
00026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00027 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00028 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00029 #include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
00030 #include "DataFormats/DetId/interface/DetId.h"
00031 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00032 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00033 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00034 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00035 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00036 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00037 //#include "DQM/SiStripCommon/interface/SiStripGenerateKey.h"
00038 
00039 //---------------------------------------------------------------------------------------------------------
00040 SiStripGainCosmicCalculator::SiStripGainCosmicCalculator(const edm::ParameterSet& iConfig) : ConditionDBWriter<SiStripApvGain>::ConditionDBWriter<SiStripApvGain>(iConfig){
00041   edm::LogInfo("SiStripGainCosmicCalculator::SiStripGainCosmicCalculator");
00042   ExpectedChargeDeposition = 200.;
00043   edm::LogInfo("SiStripApvGainCalculator::SiStripApvGainCalculator")<<"ExpectedChargeDeposition="<<ExpectedChargeDeposition;
00044 
00045   TrackProducer =  iConfig.getParameter<std::string>("TrackProducer");
00046   TrackLabel    =  iConfig.getParameter<std::string>("TrackLabel");
00047 
00048   detModulesToBeExcluded.clear(); detModulesToBeExcluded = iConfig.getParameter< std::vector<unsigned> >("detModulesToBeExcluded");
00049   MinNrEntries = iConfig.getUntrackedParameter<unsigned>("minNrEntries", 20);
00050   MaxChi2OverNDF = iConfig.getUntrackedParameter<double>("maxChi2OverNDF", 5.);
00051 
00052   outputHistogramsInRootFile = iConfig.getParameter<bool>("OutputHistogramsInRootFile");
00053   outputFileName = iConfig.getParameter<std::string>("OutputFileName");
00054 
00055   edm::LogInfo("SiStripApvGainCalculator")<<"Clusters from "<<detModulesToBeExcluded.size()<<" modules will be ignored in the calibration:";
00056   edm::LogInfo("SiStripApvGainCalculator")<<"The calibration for these DetIds will be set to a default value";
00057   for( std::vector<uint32_t>::const_iterator imod = detModulesToBeExcluded.begin(); imod != detModulesToBeExcluded.end(); imod++){
00058     edm::LogInfo("SiStripApvGainCalculator")<<"exclude detid = "<< *imod;
00059   }
00060 
00061   printdebug_ = iConfig.getUntrackedParameter<bool>("printDebug", false);
00062 }
00063 
00064 
00065 SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator(){
00066   edm::LogInfo("SiStripGainCosmicCalculator::~SiStripGainCosmicCalculator");
00067 }
00068 
00069 void SiStripGainCosmicCalculator::algoEndJob(){
00070 }
00071 
00072 void SiStripGainCosmicCalculator::algoBeginJob(const edm::EventSetup& iSetup)
00073 {
00074    eventSetupCopy_ = &iSetup;
00075    std::cout<<"SiStripGainCosmicCalculator::algoBeginJob called"<<std::endl;
00076    total_nr_of_events = 0;
00077    HlistAPVPairs = new TObjArray(); HlistOtherHistos = new TObjArray();
00078    //
00079    HlistOtherHistos->Add(new TH1F( Form("APVPairCorrections"), Form("APVPairCorrections"), 50,-1.,4.));
00080    HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB1mono"),Form("APVPairCorrectionsTIB1mono"),50,-1.,4.));
00081    HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB1stereo"),Form("APVPairCorrectionsTIB1stereo"),50,-1.,4.));
00082    HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTIB2"),Form("APVPairCorrectionsTIB2"),50,-1.,4.));
00083    HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTOB1"),Form("APVPairCorrectionsTOB1"),50,-1.,4.));
00084    HlistOtherHistos->Add(new TH1F(Form("APVPairCorrectionsTOB2"),Form("APVPairCorrectionsTOB2"),50,-1.,4.));
00085    HlistOtherHistos->Add(new TH1F(Form("LocalAngle"),Form("LocalAngle"),70,-0.1,3.4));
00086    HlistOtherHistos->Add(new TH1F(Form("LocalAngleAbsoluteCosine"),Form("LocalAngleAbsoluteCosine"),48,-0.1,1.1));
00087    HlistOtherHistos->Add(new TH1F(Form("LocalPosition_cm"),Form("LocalPosition_cm"),100,-5.,5.));
00088    HlistOtherHistos->Add(new TH1F(Form("LocalPosition_normalized"),Form("LocalPosition_normalized"),100,-1.1,1.1));
00089    TH1F* local_histo = new TH1F(Form("SiStripRecHitType"),Form("SiStripRecHitType"),2,0.5,2.5); HlistOtherHistos->Add(local_histo);
00090    local_histo->GetXaxis()->SetBinLabel(1,"simple"); local_histo->GetXaxis()->SetBinLabel(2,"matched");
00091 
00092    // get cabling and find out list of active detectors
00093    edm::ESHandle<SiStripDetCabling> siStripDetCabling; iSetup.get<SiStripDetCablingRcd>().get(siStripDetCabling);
00094    std::vector<uint32_t> activeDets; activeDets.clear();
00095    SelectedDetIds.clear();
00096    siStripDetCabling->addActiveDetectorsRawIds(activeDets);
00097 //    SelectedDetIds = activeDets; // all active detector modules
00098    // use SiStripSubStructure for selecting certain regions
00099    SiStripSubStructure substructure;
00100    substructure.getTIBDetectors(activeDets, SelectedDetIds, 0, 0, 0, 0); // this adds rawDetIds to SelectedDetIds
00101    substructure.getTOBDetectors(activeDets, SelectedDetIds, 0, 0, 0);    // this adds rawDetIds to SelectedDetIds
00102    // get tracker geometry and find nr. of apv pairs for each active detector 
00103    edm::ESHandle<TrackerGeometry> tkGeom; iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );     
00104    for(TrackerGeometry::DetContainer::const_iterator it = tkGeom->dets().begin(); it != tkGeom->dets().end(); it++){ // loop over detector modules
00105      if( dynamic_cast<StripGeomDetUnit*>((*it))!=0){
00106        uint32_t detid=((*it)->geographicalId()).rawId();
00107        // get thickness for all detector modules, not just for active, this is strange 
00108        double module_thickness = (*it)->surface().bounds().thickness(); // get thickness of detector from GeomDet (DetContainer == vector<GeomDet*>)
00109        thickness_map.insert(std::make_pair(detid,module_thickness));
00110        //
00111        bool is_active_detector = false;
00112        for(std::vector<uint32_t>::iterator iactive = SelectedDetIds.begin(); iactive != SelectedDetIds.end(); iactive++){
00113          if( *iactive == detid ){
00114            is_active_detector = true;
00115            break; // leave for loop if found matching detid
00116          }
00117        }
00118        //
00119        bool exclude_this_detid = false;
00120        for( std::vector<uint32_t>::const_iterator imod = detModulesToBeExcluded.begin(); imod != detModulesToBeExcluded.end(); imod++ ){
00121            if(*imod == detid) exclude_this_detid = true; // found in exclusion list
00122            break;
00123        }
00124        //
00125        if(is_active_detector && (!exclude_this_detid)){ // check whether is active detector and that should not be excluded
00126          const StripTopology& p = dynamic_cast<StripGeomDetUnit*>((*it))->specificTopology();
00127          unsigned short NAPVPairs = p.nstrips()/256;
00128          if( NAPVPairs<2 || NAPVPairs>3 ) {
00129            edm::LogError("SiStripGainCosmicCalculator")<<"Problem with Number of strips in detector: "<<p.nstrips()<<" Exiting program";
00130            exit(1);
00131          }
00132          for(int iapp = 0; iapp<NAPVPairs; iapp++){
00133            TString hid = Form("ChargeAPVPair_%i_%i",detid,iapp);
00134            HlistAPVPairs->Add(new TH1F(hid,hid,45,0.,1350.)); // multiply by 3 to take into account division by width
00135          }
00136        }
00137      }
00138    }
00139 }
00140 
00141 //---------------------------------------------------------------------------------------------------------
00142 void SiStripGainCosmicCalculator::algoAnalyze(const edm::Event & iEvent, const edm::EventSetup& iSetup){
00143   using namespace edm;
00144   total_nr_of_events++;
00145 
00146   //TO BE RESTORED
00147   //  anglefinder_->init(event,iSetup);
00148 
00149 
00150   // get seeds
00151 //  edm::Handle<TrajectorySeedCollection> seedcoll;
00152 //  event.getByType(seedcoll);
00153   // get tracks
00154   Handle<reco::TrackCollection> trackCollection; iEvent.getByLabel(TrackProducer, TrackLabel, trackCollection);
00155   const reco::TrackCollection *tracks=trackCollection.product();
00156 
00157 //  // get magnetic field
00158 //  edm::ESHandle<MagneticField> esmagfield;
00159 //  es.get<IdealMagneticFieldRecord>().get(esmagfield);
00160 //  magfield=&(*esmagfield);
00161   // loop over tracks
00162   for(reco::TrackCollection::const_iterator itr = tracks->begin(); itr != tracks->end(); itr++){ // looping over tracks
00163 
00164     //TO BE RESTORED
00165     //    std::vector<std::pair<const TrackingRecHit *,float> >hitangle =anglefinder_->findtrackangle((*(*seedcoll).begin()),*itr);
00166     std::vector<std::pair<const TrackingRecHit *,float> >hitangle;// =anglefinder_->findtrackangle((*(*seedcoll).begin()),*itr);
00167 
00168     for(std::vector<std::pair<const TrackingRecHit *,float> >::const_iterator hitangle_iter=hitangle.begin();hitangle_iter!=hitangle.end();hitangle_iter++){
00169       const TrackingRecHit * trechit = hitangle_iter->first;
00170       float local_angle=hitangle_iter->second;
00171       LocalPoint local_position= trechit->localPosition();
00172       const SiStripRecHit2D* sistripsimplehit=dynamic_cast<const SiStripRecHit2D*>(trechit);
00173       const SiStripMatchedRecHit2D* sistripmatchedhit=dynamic_cast<const SiStripMatchedRecHit2D*>(trechit);
00174 //      std::cout<<" hit/matched "<<std::ios::hex<<sistripsimplehit<<" "<<sistripmatchedhit<<std::endl;
00175       ((TH1F*) HlistOtherHistos->FindObject("LocalAngle"))->Fill(local_angle);
00176       ((TH1F*) HlistOtherHistos->FindObject("LocalAngleAbsoluteCosine"))->Fill(fabs(cos(local_angle)));
00177       if(sistripsimplehit){
00178         ((TH1F*) HlistOtherHistos->FindObject("SiStripRecHitType"))->Fill(1.);
00179         const SiStripRecHit2D::ClusterRef & cluster=sistripsimplehit->cluster();
00180         const std::vector<uint8_t>& ampls = cluster->amplitudes();
00181 //        const std::vector<uint16_t>& ampls = cluster->amplitudes();
00182         uint32_t thedetid  = cluster->geographicalId();
00183         double module_width = moduleWidth(thedetid, &iSetup);
00184         ((TH1F*) HlistOtherHistos->FindObject("LocalPosition_cm"))->Fill(local_position.x());
00185         ((TH1F*) HlistOtherHistos->FindObject("LocalPosition_normalized"))->Fill(local_position.x()/module_width);
00186         double module_thickness = moduleThickness(thedetid, &iSetup);
00187 //        int ifirststrip= cluster->firstStrip();
00188 //        int theapvpairid = int(float(ifirststrip)/256.);
00189         TH1F* histopointer = (TH1F*) HlistAPVPairs->FindObject(Form("ChargeAPVPair_%i_%i,thedetid,theapvpairid"));
00190         if( histopointer ){
00191           short cCharge = 0;
00192           for(unsigned int iampl = 0; iampl<ampls.size(); iampl++){
00193           cCharge += ampls[iampl];
00194           }
00195           double cluster_charge_over_path = ((double)cCharge) * fabs(cos(local_angle)) / ( 10. * module_thickness);
00196           histopointer->Fill(cluster_charge_over_path);
00197         }
00198       }else{
00199            if(sistripmatchedhit) ((TH1F*) HlistOtherHistos->FindObject("SiStripRecHitType"))->Fill(2.);
00200       }
00201     }
00202   }
00203 }
00204 
00205 
00206 //---------------------------------------------------------------------------------------------------------
00207 std::pair<double,double> SiStripGainCosmicCalculator::getPeakOfLandau( TH1F * inputHisto){ // automated fitting with finding of the appropriate nr. of ADCs
00208   // set some default dummy value and return if no entries
00209   double adcs = -0.5; double error = 0.; double nr_of_entries = inputHisto->GetEntries();
00210   if(nr_of_entries < MinNrEntries){
00211     return std::make_pair(adcs,error);
00212   }
00213 //
00214 //  // fit with initial setting of  parameter values
00215 //  double rms_of_histogram = inputHisto->GetRMS();
00216 //  TF1 *landaufit = new TF1("landaufit","landau",0.,450.);
00217 //  landaufit->SetParameters(nr_of_entries,mean_of_histogram,rms_of_histogram);
00218 //  inputHisto->Fit("landaufit","0Q+");
00219 //  delete landaufit;
00220 //
00221   // perform fit with standard landau
00222   inputHisto->Fit("landau","0Q");
00223   TF1 * fitfunction = (TF1*) inputHisto->GetListOfFunctions()->First();
00224   adcs = fitfunction->GetParameter("MPV");
00225   error = fitfunction->GetParError(1); // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
00226   double chi2 = fitfunction->GetChisquare();
00227   double ndf = fitfunction->GetNDF();
00228   double chi2overndf = chi2 / ndf;
00229   // in case things went wrong, try to refit in smaller range
00230   if(adcs< 2. || (error/adcs)>1.8 ){
00231      inputHisto->Fit("landau","0Q",0,0.,400.);
00232      TF1 * fitfunction2 = (TF1*) inputHisto->GetListOfFunctions()->First();
00233      std::cout<<"refitting landau for histogram "<<inputHisto->GetTitle()<<std::endl;
00234      std::cout<<"initial error/adcs ="<<error<<" / "<<adcs<<std::endl;
00235      std::cout<<"new     error/adcs ="<<fitfunction2->GetParError(1)<<" / "<<fitfunction2->GetParameter("MPV")<<std::endl;
00236      adcs = fitfunction2->GetParameter("MPV");
00237      error = fitfunction2->GetParError(1); // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
00238      chi2 = fitfunction2->GetChisquare();
00239      ndf = fitfunction2->GetNDF();
00240      chi2overndf = chi2 / ndf;
00241    }
00242    // if still wrong, give up
00243    if(adcs<2. || chi2overndf>MaxChi2OverNDF){
00244      adcs = -0.5; error = 0.;
00245    }
00246   return std::make_pair(adcs,error);
00247 }
00248 
00249 //---------------------------------------------------------------------------------------------------------
00250 double SiStripGainCosmicCalculator::moduleWidth(const uint32_t detid, const edm::EventSetup* iSetup) // get width of the module detid
00251 { //dk: copied from A. Giammanco and hacked,  module_width values : 10.49 12.03 6.144 7.14 9.3696
00252   edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );     
00253   double module_width=0.;
00254   const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
00255   if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
00256     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
00257   }else{
00258     module_width = it->surface().bounds().width();
00259   }
00260   return module_width;
00261 }
00262 
00263 //---------------------------------------------------------------------------------------------------------
00264 double SiStripGainCosmicCalculator::moduleThickness(const uint32_t detid, const edm::EventSetup* iSetup) // get thickness of the module detid
00265 { //dk: copied from A. Giammanco and hacked
00266   edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );
00267   double module_thickness=0.;
00268   const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
00269   if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
00270     std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
00271   }else{
00272     module_thickness = it->surface().bounds().thickness();
00273   }
00274   return module_thickness;
00275 }
00276 
00277 //---------------------------------------------------------------------------------------------------------
00278 SiStripApvGain * SiStripGainCosmicCalculator::getNewObject() {
00279   std::cout<<"SiStripGainCosmicCalculator::getNewObject called"<<std::endl;
00280 
00281   std::cout<<"total_nr_of_events="<<total_nr_of_events<<std::endl;
00282   // book some more histograms
00283   TH1F *ChargeOfEachAPVPair = new TH1F("ChargeOfEachAPVPair","ChargeOfEachAPVPair",1,0,1); ChargeOfEachAPVPair->SetBit(TH1::kCanRebin);
00284   TH1F *EntriesApvPairs = new TH1F("EntriesApvPairs","EntriesApvPairs",1,0,1); EntriesApvPairs->SetBit(TH1::kCanRebin);
00285   TH1F * NrOfEntries = new TH1F("NrOfEntries","NrOfEntries",351,-0.5,350.5);// NrOfEntries->SetBit(TH1::kCanRebin);
00286   TH1F * ModuleThickness = new TH1F("ModuleThickness","ModuleThickness",2,0.5,2.5); HlistOtherHistos->Add(ModuleThickness);
00287   ModuleThickness->GetXaxis()->SetBinLabel(1,"320mu"); ModuleThickness->GetXaxis()->SetBinLabel(2,"500mu"); ModuleThickness->SetYTitle("Nr APVPairs");
00288   TH1F * ModuleWidth = new TH1F("ModuleWidth","ModuleWidth",5,0.5,5.5); HlistOtherHistos->Add(ModuleWidth);
00289   ModuleWidth->GetXaxis()->SetBinLabel(1,"6.144cm"); ModuleWidth->GetXaxis()->SetBinLabel(2,"7.14cm");
00290   ModuleWidth->GetXaxis()->SetBinLabel(3,"9.3696cm"); ModuleWidth->GetXaxis()->SetBinLabel(4,"10.49cm");
00291   ModuleWidth->GetXaxis()->SetBinLabel(5,"12.03cm");
00292   ModuleWidth->SetYTitle("Nr APVPairs");
00293   // loop over single histograms and extract peak value of charge
00294   HlistAPVPairs->Sort(); // sort alfabetically
00295   TIter hiterator(HlistAPVPairs);
00296   double MeanCharge = 0.;
00297   double NrOfApvPairs = 0.;
00298   TH1F *MyHisto = (TH1F*)hiterator();
00299   while( MyHisto ){
00300     TString histo_title = MyHisto->GetTitle();
00301     if(histo_title.Contains("ChargeAPVPair_")){
00302       std::pair<double,double> two_values = getPeakOfLandau(MyHisto);
00303       double local_nrofadcs = two_values.first;
00304       double local_sigma = two_values.second;
00305       ChargeOfEachAPVPair->Fill(histo_title, local_nrofadcs);
00306       int ichbin  = ChargeOfEachAPVPair->GetXaxis()->FindBin(histo_title.Data());
00307       ChargeOfEachAPVPair->SetBinError(ichbin,local_sigma);
00308       EntriesApvPairs->Fill(histo_title, MyHisto->GetEntries());
00309       NrOfEntries->Fill(MyHisto->GetEntries());
00310       if(local_nrofadcs > 0){ // if nr of adcs is negative, the fitting routine could not extract meaningfull numbers
00311        MeanCharge += local_nrofadcs;
00312        NrOfApvPairs += 1.; // count nr of apv pairs since do not know whether nr of bins of histogram is the same
00313       }
00314     }
00315     MyHisto = (TH1F*)hiterator();
00316   }
00317   ChargeOfEachAPVPair->LabelsDeflate("X"); EntriesApvPairs->LabelsDeflate("X"); // trim nr. of bins to match active labels
00318   HlistOtherHistos->Add(ChargeOfEachAPVPair);
00319   HlistOtherHistos->Add(EntriesApvPairs);
00320   HlistOtherHistos->Add(NrOfEntries); 
00321   MeanCharge = MeanCharge / NrOfApvPairs;
00322   // calculate correction
00323   TH1F* CorrectionOfEachAPVPair = (TH1F*) ChargeOfEachAPVPair->Clone("CorrectionOfEachAPVPair");
00324   TH1F *ChargeOfEachAPVPairControlView = new TH1F("ChargeOfEachAPVPairControlView","ChargeOfEachAPVPairControlView",1,0,1); ChargeOfEachAPVPairControlView->SetBit(TH1::kCanRebin);
00325 TH1F *CorrectionOfEachAPVPairControlView = new TH1F("CorrectionOfEachAPVPairControlView","CorrectionOfEachAPVPairControlView",1,0,1); CorrectionOfEachAPVPairControlView->SetBit(TH1::kCanRebin);
00326   std::ofstream APVPairTextOutput("apvpair_corrections.txt");
00327   APVPairTextOutput<<"# MeanCharge = "<<MeanCharge<<std::endl;
00328   APVPairTextOutput<<"# Nr. of APVPairs = "<<NrOfApvPairs<<std::endl;
00329   for(int ibin=1; ibin <= ChargeOfEachAPVPair->GetNbinsX(); ibin++){
00330      TString local_bin_label = ChargeOfEachAPVPair->GetXaxis()->GetBinLabel(ibin);
00331      double local_charge_over_path = ChargeOfEachAPVPair->GetBinContent(ibin);
00332      if(local_bin_label.Contains("ChargeAPVPair_") && local_charge_over_path > 0.0000001){ // calculate correction only for meaningful numbers
00333        uint32_t extracted_detid; std::istringstream read_label((local_bin_label(14,9)).Data()); read_label >> extracted_detid; 
00334        unsigned short extracted_apvpairid; std::istringstream read_apvpair((local_bin_label(24,1)).Data()); read_apvpair >> extracted_apvpairid; 
00335        double local_error_of_charge = ChargeOfEachAPVPair->GetBinError(ibin);
00336        double local_correction = -0.5;
00337        double local_error_correction = 0.;
00338        local_correction = MeanCharge / local_charge_over_path; // later use ExpectedChargeDeposition instead of MeanCharge
00339        local_error_correction = local_correction * local_error_of_charge / local_charge_over_path;
00340        if(local_error_correction>1.8){ // understand why error too large sometimes
00341          std::cout<<"too large error "<<local_error_correction<<" for histogram "<<local_bin_label<<std::endl;
00342        }
00343        double nr_of_entries = EntriesApvPairs->GetBinContent(ibin);
00344        APVPairTextOutput<<local_bin_label<<" "<<local_correction<<" "<<local_charge_over_path<<" "<<nr_of_entries<<std::endl;
00345        CorrectionOfEachAPVPair->SetBinContent(ibin, local_correction);
00346        CorrectionOfEachAPVPair->SetBinError(ibin, local_error_correction);
00347        ((TH1F*) HlistOtherHistos->FindObject("APVPairCorrections"))->Fill(local_correction);
00348        DetId thedetId = DetId(extracted_detid);
00349        uint generalized_layer = 0;
00350        // calculate generalized_layer:  31,32 = TIB1, 33 = TIB2, 33 = TIB3, 51 = TOB1, 52 = TOB2, 60 = TEC
00351        if(thedetId.subdetId()==StripSubdetector::TIB){
00352           TIBDetId ptib = TIBDetId(thedetId.rawId());
00353           generalized_layer = 10*thedetId.subdetId() + ptib.layer() + ptib.stereo();
00354           if(ptib.layer()==2){
00355             generalized_layer++;
00356             if (ptib.glued()) edm::LogError("ClusterMTCCFilter")<<"WRONGGGG"<<std::endl;
00357           }
00358         }else{
00359           generalized_layer = 10*thedetId.subdetId();
00360           if(thedetId.subdetId()==StripSubdetector::TOB){
00361             TOBDetId ptob = TOBDetId(thedetId.rawId());
00362             generalized_layer += ptob.layer();
00363           }
00364         }
00365        if(generalized_layer==31){
00366          ((TH1F*) HlistOtherHistos->FindObject("APVPairCorrectionsTIB1mono"))->Fill(local_correction);
00367        }
00368        if(generalized_layer==32){
00369          ((TH1F*) HlistOtherHistos->FindObject("APVPairCorrectionsTIB1stereo"))->Fill(local_correction);
00370        }
00371        if(generalized_layer==33){
00372         ((TH1F*) HlistOtherHistos->FindObject("APVPairCorrectionsTIB2"))->Fill(local_correction);
00373        }
00374        if(generalized_layer==51){
00375         ((TH1F*) HlistOtherHistos->FindObject("APVPairCorrectionsTOB1"))->Fill(local_correction);
00376        }
00377        if(generalized_layer==52){
00378         ((TH1F*) HlistOtherHistos->FindObject("APVPairCorrectionsTOB2"))->Fill(local_correction);
00379        }
00380        // control view
00381        edm::ESHandle<SiStripDetCabling> siStripDetCabling; eventSetupCopy_->get<SiStripDetCablingRcd>().get(siStripDetCabling);
00382        const FedChannelConnection& fedchannelconnection = siStripDetCabling->getConnection( extracted_detid, extracted_apvpairid );
00383        std::ostringstream local_key;
00384        // in S. Mersi's analysis the APVPair id seems to be used instead of the lldChannel, hence use the same here
00385        local_key<<"fecCrate"<<fedchannelconnection.fecCrate()<<"_fecSlot"<<fedchannelconnection.fecSlot()<<"_fecRing"<<fedchannelconnection.fecRing()<<"_ccuAddr"<<fedchannelconnection.ccuAddr()<<"_ccuChan"<<fedchannelconnection.ccuChan()<<"_apvPair"<<extracted_apvpairid;
00386        TString control_key = local_key.str();
00387        ChargeOfEachAPVPairControlView->Fill(control_key,local_charge_over_path);
00388        int ibin1  = ChargeOfEachAPVPairControlView->GetXaxis()->FindBin(control_key);
00389        ChargeOfEachAPVPairControlView->SetBinError(ibin1,local_error_of_charge);
00390        CorrectionOfEachAPVPairControlView->Fill(control_key, local_correction);
00391        int ibin2  = CorrectionOfEachAPVPairControlView->GetXaxis()->FindBin(control_key);
00392        CorrectionOfEachAPVPairControlView->SetBinError(ibin2, local_error_correction);
00393        // thickness of each module
00394        double module_thickness = moduleThickness(extracted_detid, eventSetupCopy_);
00395        if( fabs(module_thickness - 0.032)<0.001 ) ModuleThickness->Fill(1);
00396        if( fabs(module_thickness - 0.05)<0.001 )  ModuleThickness->Fill(2);
00397        // width of each module     
00398        double module_width = moduleWidth(extracted_detid, eventSetupCopy_);
00399        if(fabs(module_width-6.144)<0.01) ModuleWidth->Fill(1);
00400        if(fabs(module_width-7.14)<0.01) ModuleWidth->Fill(2);
00401        if(fabs(module_width-9.3696)<0.01) ModuleWidth->Fill(3);
00402        if(fabs(module_width-10.49)<0.01) ModuleWidth->Fill(4);
00403        if(fabs(module_width-12.03)<0.01) ModuleWidth->Fill(5);
00404      }
00405   }
00406   HlistOtherHistos->Add(CorrectionOfEachAPVPair);
00407   ChargeOfEachAPVPairControlView->LabelsDeflate("X");
00408   CorrectionOfEachAPVPairControlView->LabelsDeflate("X");
00409   HlistOtherHistos->Add(ChargeOfEachAPVPairControlView);
00410   HlistOtherHistos->Add(CorrectionOfEachAPVPairControlView);
00411   // output histograms to file
00412 
00413 
00414   if(outputHistogramsInRootFile){
00415     TFile *outputfile = new TFile(outputFileName,"RECREATE");
00416     HlistAPVPairs->Write();
00417     HlistOtherHistos->Write();
00418     outputfile->Close();
00419   }
00420 
00421   SiStripApvGain * obj = new SiStripApvGain();
00422 
00423 //   for(std::map<uint32_t,OptoScanAnalysis*>::const_iterator it = analyses.begin(); it != analyses.end(); it++){
00424 //     //Generate Gain for det detid
00425 //     std::vector<float> theSiStripVector;
00426 //     for(unsigned short j=0; j<it->second; j++){
00427 //       float gain;
00428 
00429 //       //      if(sigmaGain_/meanGain_ < 0.00001) gain = meanGain_;
00430 //       //      else{
00431 //       gain = RandGauss::shoot(meanGain_, sigmaGain_);
00432 //       if(gain<=minimumPosValue_) gain=minimumPosValue_;
00433 //       //      }
00434 
00435 //       if (printdebug_)
00436 //      edm::LogInfo("SiStripGainCalculator") << "detid " << it->first << " \t"
00437 //                                            << " apv " << j << " \t"
00438 //                                            << gain    << " \t" 
00439 //                                            << std::endl;         
00440 //       theSiStripVector.push_back(gain);
00441 //     }
00442 //     SiStripApvGain::Range range(theSiStripVector.begin(),theSiStripVector.end());
00443 //     if ( ! obj->put(it->first,range) )
00444 //       edm::LogError("SiStripGainCalculator")<<"[SiStripGainCalculator::beginJob] detid already exists"<<std::endl;
00445 //   }
00446   
00447   return obj;
00448 }
00449