00001 #include "Validation/HcalHits/interface/HcalSimHitStudy.h"
00002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00003
00004 #include "FWCore/Utilities/interface/Exception.h"
00005
00006 HcalSimHitStudy::HcalSimHitStudy(const edm::ParameterSet& ps) {
00007
00008 g4Label = ps.getUntrackedParameter<std::string>("moduleLabel","g4SimHits");
00009 hcalHits = ps.getUntrackedParameter<std::string>("HitCollection","HcalHits");
00010 outFile_ = ps.getUntrackedParameter<std::string>("outputFile", "hcHit.root");
00011 verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
00012 checkHit_= true;
00013
00014 edm::LogInfo("HcalSim") << "Module Label: " << g4Label << " Hits: "
00015 << hcalHits << " / "<< checkHit_
00016 << " Output: " << outFile_;
00017
00018 dbe_ = edm::Service<DQMStore>().operator->();
00019 if (dbe_) {
00020 if (verbose_) {
00021 dbe_->setVerbose(1);
00022 sleep (3);
00023 dbe_->showDirStructure();
00024 } else {
00025 dbe_->setVerbose(0);
00026 }
00027 }
00028 }
00029
00030 HcalSimHitStudy::~HcalSimHitStudy() {}
00031
00032 void HcalSimHitStudy::beginJob() {
00033
00034 if (dbe_) {
00035 dbe_->setCurrentFolder("HcalHitsV/HcalSimHitsTask");
00036
00037
00038 if (checkHit_) {
00039 meAllNHit_ = dbe_->book1D("Hit01","Number of Hits in HCal",1000,0.,1000.);
00040 meBadDetHit_= dbe_->book1D("Hit02","Hits with wrong Det", 100,0.,100.);
00041 meBadSubHit_= dbe_->book1D("Hit03","Hits with wrong Subdet",100,0.,100.);
00042 meBadIdHit_ = dbe_->book1D("Hit04","Hits with wrong ID", 100,0.,100.);
00043 meHBNHit_ = dbe_->book1D("Hit05","Number of Hits in HB",1000,0.,1000.);
00044 meHENHit_ = dbe_->book1D("Hit06","Number of Hits in HE",1000,0.,1000.);
00045 meHONHit_ = dbe_->book1D("Hit07","Number of Hits in HO",1000,0.,1000.);
00046 meHFNHit_ = dbe_->book1D("Hit08","Number of Hits in HF",1000,0.,1000.);
00047 meDetectHit_= dbe_->book1D("Hit09","Detector ID", 50,0.,50.);
00048 meSubdetHit_= dbe_->book1D("Hit10","Subdetectors in HCal", 50,0.,50.);
00049 meDepthHit_ = dbe_->book1D("Hit11","Depths in HCal", 20,0.,20.);
00050 meEtaHit_ = dbe_->book1D("Hit12","Eta in HCal", 100,-50.,50.);
00051 mePhiHit_ = dbe_->book1D("Hit13","Phi in HCal", 100,0.,100.);
00052 meEnergyHit_= dbe_->book1D("Hit14","Energy in HCal", 100,0.,1.);
00053 meTimeHit_ = dbe_->book1D("Hit15","Time in HCal", 100,0.,400.);
00054 meTimeWHit_ = dbe_->book1D("Hit16","Time in HCal (E wtd)", 100,0.,400.);
00055 meHBDepHit_ = dbe_->book1D("Hit17","Depths in HB", 20,0.,20.);
00056 meHEDepHit_ = dbe_->book1D("Hit18","Depths in HE", 20,0.,20.);
00057 meHODepHit_ = dbe_->book1D("Hit19","Depths in HO", 20,0.,20.);
00058 meHFDepHit_ = dbe_->book1D("Hit20","Depths in HF", 20,0.,20.);
00059 meHBEtaHit_ = dbe_->book1D("Hit21","Eta in HB", 100,-50.,50.);
00060 meHEEtaHit_ = dbe_->book1D("Hit22","Eta in HE", 100,-50.,50.);
00061 meHOEtaHit_ = dbe_->book1D("Hit23","Eta in HO", 100,-50.,50.);
00062 meHFEtaHit_ = dbe_->book1D("Hit24","Eta in HF", 100,-50.,50.);
00063 meHBPhiHit_ = dbe_->book1D("Hit25","Phi in HB", 100,0.,100.);
00064 meHEPhiHit_ = dbe_->book1D("Hit26","Phi in HE", 100,0.,100.);
00065 meHOPhiHit_ = dbe_->book1D("Hit27","Phi in HO", 100,0.,100.);
00066 meHFPhiHit_ = dbe_->book1D("Hit28","Phi in HF", 100,0.,100.);
00067 meHBEneHit_ = dbe_->book1D("Hit29","Energy in HB", 100,0.,1.);
00068 meHEEneHit_ = dbe_->book1D("Hit30","Energy in HE", 100,0.,1.);
00069 meHOEneHit_ = dbe_->book1D("Hit31","Energy in HO", 100,0.,1.);
00070 meHFEneHit_ = dbe_->book1D("Hit32","Energy in HF", 100,0.,100.);
00071 meHBTimHit_ = dbe_->book1D("Hit33","Time in HB", 100,0.,400.);
00072 meHETimHit_ = dbe_->book1D("Hit34","Time in HE", 100,0.,400.);
00073 meHOTimHit_ = dbe_->book1D("Hit35","Time in HO", 100,0.,400.);
00074 meHFTimHit_ = dbe_->book1D("Hit36","Time in HF", 100,0.,400.);
00075 meHBEneHit2_ = dbe_->book1D("Hit37","Energy in HB 2", 100,0.,0.0001);
00076 meHEEneHit2_ = dbe_->book1D("Hit38","Energy in HE 2", 100,0.,0.0001);
00077 meHOEneHit2_ = dbe_->book1D("Hit39","Energy in HO 2", 100,0.,0.0001);
00078 meHFEneHit2_ = dbe_->book1D("Hit40","Energy in HF 2", 100,0.,0.0001);
00079 meHBL10Ene_ = dbe_->book1D("Hit41","Log10Energy in HB", 140, -10., 4. );
00080 meHEL10Ene_ = dbe_->book1D("Hit42","Log10Energy in HE", 140, -10., 4. );
00081 meHFL10Ene_ = dbe_->book1D("Hit43","Log10Energy in HF", 140, -10., 4. );
00082 meHOL10Ene_ = dbe_->book1D("Hit44","Log10Energy in HO", 140, -10., 4. );
00083 meHBL10EneP_ = dbe_->bookProfile("Hit45","Log10Energy in HB vs Hit contribution", 140, -10., 4., 100, 0., 1. );
00084 meHEL10EneP_ = dbe_->bookProfile("Hit46","Log10Energy in HE vs Hit contribution", 140, -10., 4., 100, 0., 1. );
00085 meHFL10EneP_ = dbe_->bookProfile("Hit47","Log10Energy in HF vs Hit contribution", 140, -10., 4., 100, 0., 1. );
00086 meHOL10EneP_ = dbe_->bookProfile("Hit48","Log10Energy in HO vs Hit contribution", 140, -10., 4., 100, 0., 1. );
00087 }
00088 }
00089 }
00090
00091 void HcalSimHitStudy::endJob() {
00092 if (dbe_ && outFile_.size() > 0) dbe_->save(outFile_);
00093 }
00094
00095 void HcalSimHitStudy::analyze(const edm::Event& e, const edm::EventSetup& ) {
00096
00097 LogDebug("HcalSim") << "Run = " << e.id().run() << " Event = "
00098 << e.id().event();
00099
00100 std::vector<PCaloHit> caloHits;
00101 edm::Handle<edm::PCaloHitContainer> hitsHcal;
00102
00103 bool getHits = false;
00104 if (checkHit_) {
00105 e.getByLabel(g4Label,hcalHits,hitsHcal);
00106 if (hitsHcal.isValid()) getHits = true;
00107 }
00108
00109 LogDebug("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
00110
00111 if (getHits) {
00112 caloHits.insert(caloHits.end(),hitsHcal->begin(),hitsHcal->end());
00113 LogDebug("HcalSim") << "HcalValidation: Hit buffer "
00114 << caloHits.size();
00115 analyzeHits (caloHits);
00116 }
00117 }
00118
00119 void HcalSimHitStudy::analyzeHits (std::vector<PCaloHit>& hits) {
00120
00121 int nHit = hits.size();
00122 int nHB=0, nHE=0, nHO=0, nHF=0, nBad1=0, nBad2=0, nBad=0;
00123 std::vector<double> encontHB(140, 0.);
00124 std::vector<double> encontHE(140, 0.);
00125 std::vector<double> encontHF(140, 0.);
00126 std::vector<double> encontHO(140, 0.);
00127 double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0;
00128
00129 for (int i=0; i<nHit; i++) {
00130 double energy = hits[i].energy();
00131 double log10en = log10(energy);
00132 int log10i = int( (log10en+10.)*10. );
00133 double time = hits[i].time();
00134 unsigned int id_ = hits[i].id();
00135 HcalDetId id = HcalDetId(id_);
00136 int det = id.det();
00137 int subdet = id.subdet();
00138 int depth = id.depth();
00139 int eta = id.ieta();
00140 int phi = id.iphi();
00141 LogDebug("HcalSim") << "Hit[" << i << "] ID " << std::hex << id_
00142 << std::dec << " Det " << det << " Sub "
00143 << subdet << " depth " << depth << " Eta " << eta
00144 << " Phi " << phi << " E " << energy << " time "
00145 << time;
00146 if (det == 4) {
00147 if (subdet == static_cast<int>(HcalBarrel)) nHB++;
00148 else if (subdet == static_cast<int>(HcalEndcap)) nHE++;
00149 else if (subdet == static_cast<int>(HcalOuter)) nHO++;
00150 else if (subdet == static_cast<int>(HcalForward)) nHF++;
00151 else { nBad++; nBad2++;}
00152 } else { nBad++; nBad1++;}
00153 if (dbe_) {
00154 meDetectHit_->Fill(double(det));
00155 if (det == 4) {
00156 meSubdetHit_->Fill(double(subdet));
00157 meDepthHit_->Fill(double(depth));
00158 meEtaHit_->Fill(double(eta));
00159 mePhiHit_->Fill(double(phi));
00160 meEnergyHit_->Fill(energy);
00161 meTimeHit_->Fill(time);
00162 meTimeWHit_->Fill(double(time),energy);
00163 if (subdet == static_cast<int>(HcalBarrel)) {
00164 meHBDepHit_->Fill(double(depth));
00165 meHBEtaHit_->Fill(double(eta));
00166 meHBPhiHit_->Fill(double(phi));
00167 meHBEneHit_->Fill(energy);
00168 meHBEneHit2_->Fill(energy);
00169 meHBTimHit_->Fill(time);
00170 meHBL10Ene_->Fill(log10en);
00171 if( log10i >=0 && log10i < 140 ) encontHB[log10i] += energy;
00172 entotHB += energy;
00173 } else if (subdet == static_cast<int>(HcalEndcap)) {
00174 meHEDepHit_->Fill(double(depth));
00175 meHEEtaHit_->Fill(double(eta));
00176 meHEPhiHit_->Fill(double(phi));
00177 meHEEneHit_->Fill(energy);
00178 meHEEneHit2_->Fill(energy);
00179 meHETimHit_->Fill(time);
00180 meHEL10Ene_->Fill(log10en);
00181 if( log10i >=0 && log10i < 140 ) encontHE[log10i] += energy;
00182 entotHE += energy;
00183 } else if (subdet == static_cast<int>(HcalOuter)) {
00184 meHODepHit_->Fill(double(depth));
00185 meHOEtaHit_->Fill(double(eta));
00186 meHOPhiHit_->Fill(double(phi));
00187 meHOEneHit_->Fill(energy);
00188 meHOEneHit2_->Fill(energy);
00189 meHOTimHit_->Fill(time);
00190 meHOL10Ene_->Fill(log10en);
00191 if( log10i >=0 && log10i < 140 ) encontHO[log10i] += energy;
00192 entotHO += energy;
00193 } else if (subdet == static_cast<int>(HcalForward)) {
00194 meHFDepHit_->Fill(double(depth));
00195 meHFEtaHit_->Fill(double(eta));
00196 meHFPhiHit_->Fill(double(phi));
00197 meHFEneHit_->Fill(energy);
00198 meHFEneHit2_->Fill(energy);
00199 meHFTimHit_->Fill(time);
00200 meHFL10Ene_->Fill(log10en);
00201 if( log10i >=0 && log10i < 140 ) encontHF[log10i] += energy;
00202 entotHF += energy;
00203 }
00204 }
00205 }
00206 }
00207 if( entotHB != 0 ) for( int i=0; i<140; i++ ) meHBL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHB[i]/entotHB );
00208 if( entotHE != 0 ) for( int i=0; i<140; i++ ) meHEL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHE[i]/entotHE );
00209 if( entotHF != 0 ) for( int i=0; i<140; i++ ) meHFL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHF[i]/entotHF );
00210 if( entotHO != 0 ) for( int i=0; i<140; i++ ) meHOL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHO[i]/entotHO );
00211
00212 if (dbe_) {
00213 meAllNHit_->Fill(double(nHit));
00214 meBadDetHit_->Fill(double(nBad1));
00215 meBadSubHit_->Fill(double(nBad2));
00216 meBadIdHit_->Fill(double(nBad));
00217 meHBNHit_->Fill(double(nHB));
00218 meHENHit_->Fill(double(nHE));
00219 meHONHit_->Fill(double(nHO));
00220 meHFNHit_->Fill(double(nHF));
00221 }
00222 LogDebug("HcalSim") << "HcalSimHitStudy::analyzeHits: HB " << nHB
00223 << " HE " << nHE << " HO " << nHO << " HF " << nHF
00224 << " Bad " << nBad << " All " << nHit;
00225
00226 }
00227
1.7.3