Reference Manual

Go to the documentation of this file.
00001 #include "Validation/HcalHits/interface/SimHitsValidationHcal.h"
00002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00003 
00004 SimHitsValidationHcal::SimHitsValidationHcal(const edm::ParameterSet& ps) {
00005 
00006   g4Label  = ps.getUntrackedParameter<std::string>("moduleLabel","g4SimHits");
00007   hcalHits = ps.getUntrackedParameter<std::string>("HitCollection","HcalHits");
00008   verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
00009 
00010   edm::LogInfo("HitsValidationHcal") << "Module Label: " << g4Label << "   Hits: "
00011                                      << hcalHits;
00012 
00013   dbe_ = edm::Service<DQMStore>().operator->();
00014   if (dbe_) {
00015     if (verbose_) {
00016       dbe_->setVerbose(1);
00017       sleep (3);
00018       dbe_->showDirStructure();
00019     } else {
00020       dbe_->setVerbose(0);
00021     }
00022   }
00023 }
00024 
00025 SimHitsValidationHcal::~SimHitsValidationHcal() {}
00026 
00027 void SimHitsValidationHcal::beginJob() {
00028   
00029   if (dbe_) {
00030     edm::LogInfo("HitsValidationHcal") << "Histograms booked";
00031     dbe_->setCurrentFolder("HcalHitsV/SimHitsValidationHcal");
00032 
00033     //Histograms for Hits
00034     
00035     std::string divisions[nType]={"HB0","HB1","HE0+z","HE1+z","HE2+z","HE0-z","HE1-z",
00036                                   "HE2-z","HO0","HFL0+z","HFS0+z","HFL1+z","HFS1+z",
00037                                   "HFL2+z","HFS2+z","HFL3+z","HFS3+z","HFL0-z","HFS0-z",
00038                                   "HFL1-z","HFS1-z","HFL2-z","HFS2-z","HFL3-z","HFS3-z"};
00039     
00040     std::string divisions1[nType]={"HB depth1","HB depth2 ","HE+z depth1","HE+z depth2","HE+z depth3","HE-z depth1","HE-z depth2",
00041                                    "HE-z depth3","HO depth1","HFL+z depth1","HFS+z depth1","HFL+z depth2","HFS+z depth2",
00042                                    "HFL+z depth3","HFS+z depth3","HFL+z depth4","HFS+z depth4","HFL-z depth1","HFS-z depth1",
00043                                    "HFL-z depth2","HFS-z depth2","HFL-z depth3","HFS-z depth3 ","HFL-z depth4","HFS-z depth4"};
00044     
00045     double etaLow[nType]={-16,-16,16,16,16,-30,-30,-30,-15,29,29,29,29,29,29,29,29,
00046                           -41,-41,-41,-41,-41,-41,-41,-41};
00047     double etaHigh[nType]={16,16,30,30,30,-16,-16,-16,15,41,41,41,41,41,41,41,41,
00048                            -29,-29,-29,-29,-29,-29,-29,-29};
00049     int etaBins[nType]={32,32,14,14,14,14,14,14,30,12,12,12,12,12,12,12,12,
00050                         12,12,12,12,12,12,12,12};
00051     char name[40], title[100];
00052     
00053     for (int i=0; i<nType; ++i) {
00054       sprintf (name, "HcalHitEta%s", divisions[i].c_str());
00055       sprintf (title, "Hit energy as a function of eta tower index in %s", divisions1[i].c_str());
00056       meHcalHitEta_[i] = dbe_->book1D(name, title, etaBins[i], etaLow[i], etaHigh[i]);
00057       
00058       sprintf (name, "HcalHitTimeAEta%s", divisions[i].c_str());
00059       sprintf (title, "Hit time as a function of eta tower index in %s", divisions1[i].c_str());
00060       meHcalHitTimeEta_[i] = dbe_->book1D(name, title, etaBins[i], etaLow[i], etaHigh[i]);
00061       
00062       sprintf (name, "HcalHitE25%s", divisions[i].c_str());
00063       sprintf (title, "Energy in time window 0 to 25 for a tower in %s", divisions1[i].c_str());
00064       meHcalEnergyl25_[i] = dbe_->book2D(name, title, etaBins[i], etaLow[i], etaHigh[i], 72, 0., 72.);
00065       
00066       sprintf (name, "HcalHitE50%s", divisions[i].c_str());
00067       sprintf (title, "Energy in time window 0 to 50 for a tower in %s", divisions1[i].c_str());
00068       meHcalEnergyl50_[i] = dbe_->book2D(name, title, etaBins[i], etaLow[i], etaHigh[i], 72, 0., 72.);
00069       
00070       sprintf (name, "HalHitE100%s", divisions[i].c_str());
00071       sprintf (title, "Energy in time window 0 to 100 for a tower in %s", divisions1[i].c_str());
00072       meHcalEnergyl100_[i] = dbe_->book2D(name, title, etaBins[i], etaLow[i], etaHigh[i], 72, 0., 72.);
00073       
00074       sprintf (name, "HcalHitE250%s", divisions[i].c_str());
00075       sprintf (title, "Energy in time window 0 to 250 for a tower in %s", divisions1[i].c_str());
00076       meHcalEnergyl250_[i] = dbe_->book2D(name, title, etaBins[i], etaLow[i], etaHigh[i], 72, 0., 72.);
00077     }
00078 
00079     sprintf (name, "Energy_HB");
00080     meEnergy_HB = dbe_->book1D(name, name, 100,0,1);
00081     sprintf (name, "Energy_HE");
00082     meEnergy_HE = dbe_->book1D(name, name, 100,0,1);
00083     sprintf (name, "Energy_HO");
00084     meEnergy_HO = dbe_->book1D(name, name, 100,0,1);
00085     sprintf (name, "Energy_HF");
00086     meEnergy_HF = dbe_->book1D(name, name, 100,0,50);
00087     
00088     sprintf (name, "Time_HB");
00089     metime_HB = dbe_->book1D(name, name, 300,-150,150);
00090     sprintf (name, "Time_HE");
00091     metime_HE = dbe_->book1D(name, name, 300,-150,150);
00092     sprintf (name, "Time_HO");
00093     metime_HO = dbe_->book1D(name, name, 300,-150, 150);
00094     sprintf (name, "Time_HF");
00095     metime_HF = dbe_->book1D(name, name, 300,-150,150);
00096 
00097     sprintf (name, "Time_Enweighted_HB");
00098     metime_enweighted_HB = dbe_->book1D(name, name, 300,-150,150);
00099     sprintf (name, "Time_Enweighted_HE");
00100     metime_enweighted_HE = dbe_->book1D(name, name, 300,-150,150);
00101     sprintf (name, "Time_Enweighted_HO");
00102     metime_enweighted_HO = dbe_->book1D(name, name, 300,-150, 150);
00103     sprintf (name, "Time_Enweighted_HF");
00104     metime_enweighted_HF = dbe_->book1D(name, name, 300,-150,150);
00105   }
00106 }
00107 
00108 void SimHitsValidationHcal::endJob() {}
00109 
00110 void SimHitsValidationHcal::analyze(const edm::Event& e, const edm::EventSetup& ) {
00111   
00112   
00113   LogDebug("SimHitsValidationHcal") << "Run = " << e.id().run() << " Event = " 
00114                                     << e.id().event();
00115   
00116   std::vector<PCaloHit>               caloHits;
00117   edm::Handle<edm::PCaloHitContainer> hitsHcal;
00118   
00119   bool getHits = false;
00120   e.getByLabel(g4Label,hcalHits,hitsHcal); 
00121   if (hitsHcal.isValid()) getHits = true;
00122   
00123   LogDebug("SimHitsValidationHcal") << "SimHitsValidationHcal.: Input flags Hits " << getHits;
00124 
00125   if (getHits) {
00126     caloHits.insert(caloHits.end(),hitsHcal->begin(),hitsHcal->end());
00127     LogDebug("SimHitsValidationHcal") << "SimHitsValidationHcal: Hit buffer " 
00128                                       << caloHits.size(); 
00129     analyzeHits (caloHits);
00130   }
00131 }
00132 
00133 void SimHitsValidationHcal::analyzeHits (std::vector<PCaloHit>& hits) {
00134 
00135   int nHit = hits.size();
00136   double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0; 
00137   double timetotHB = 0, timetotHE = 0, timetotHF = 0, timetotHO = 0; 
00138   int    nHB=0, nHE=0, nHO=0, nHF=0;
00139   
00140   std::map<std::pair<HcalDetId,int>,energysum> map_try;
00141   map_try.clear();
00142   
00143   std::map<std::pair<HcalDetId,int>,energysum>::iterator itr;
00144   
00145   for (int i=0; i<nHit; i++) {
00146     double energy    = hits[i].energy();
00147     double time      = hits[i].time();
00148     unsigned int id_ = hits[i].id();
00149     HcalDetId id     = HcalDetId(id_);
00150     int itime        = (int)(time);
00151     int det          = id.det();
00152     int subdet       = id.subdet();
00153     int depth        = id.depth();
00154     int eta          = id.ieta();
00155     int phi          = id.iphi();
00156     unsigned int dep = hits[i].depth();
00157         
00158     int type         =-1;
00159     if (subdet == static_cast<int>(HcalBarrel)) {
00160       entotHB += energy;
00161       timetotHB += time;
00162       nHB++;
00163       type     = depth-1;
00164     } else if (subdet == static_cast<int>(HcalEndcap)) {
00165       entotHE += energy;
00166       timetotHE += time;
00167       nHE++;
00168       type     = depth+1;
00169       if (eta < 0) type += 3;
00170     } else if (subdet == static_cast<int>(HcalOuter)) {
00171       entotHO += energy;
00172       timetotHO += time;
00173       nHO++;
00174       type = 8;
00175     } else if (subdet == static_cast<int>(HcalForward)) {
00176       entotHF += energy;
00177       timetotHF += time;
00178       nHF++;
00179       type     = depth+8+2*dep;
00180       if (eta < 0) type += 8;
00181     }
00182 
00183     std::pair<HcalDetId,int> id0(id,type);
00184     //   LogDebug("HitsValidationHcal") << "Id and type " << id << ":" << type;
00185     energysum ensum;
00186     if (map_try.count(id0) != 0) ensum =  map_try[id0];
00187     if (itime<250) {
00188       ensum.e250 += energy;
00189       if (itime<100) {
00190         ensum.e100 += energy;
00191         if (itime<50) {
00192           ensum.e50 += energy;
00193           if (itime<25) ensum.e25 += energy;
00194         }
00195       }
00196     }
00197     map_try[id0] = ensum;
00198 
00199     LogDebug("SimHitsValidationHcal") << "Hit[" << i << "] ID " << std::hex << id_ 
00200                                       << std::dec << " " << id << std::dec<< " Det " << det << " Sub " 
00201                                       << subdet << " depth " << depth << " depthX "
00202                                       << dep << " Eta " << eta << " Phi " << phi 
00203                                       << " E " << energy << " time " << time
00204                                       << " type " << type;
00205 
00206     //    LogDebug("HitsValidationHcal") << "Hit[" << i << "] ID " << std::hex << id_ << "ID="<<std::dec << id << " Det " << det << " Sub " << subdet << " depth " << depth << " depthX " << dep << " Eta " << eta << " Phi " << phi << " E  " << energy << " time  " << time <<" itime  " << itime << " type " << type;
00207     
00208     double etax = eta - 0.5;
00209     if (eta < 0) etax += 1;
00210     if (dbe_ && type >= 0) {
00211       meHcalHitEta_[type]->Fill(etax,energy);
00212       meHcalHitTimeEta_[type]->Fill(etax,time);
00213     }
00214   }
00215   
00216   if (dbe_) {
00217     meEnergy_HB->Fill(entotHB);
00218     meEnergy_HE->Fill(entotHE);
00219     meEnergy_HF->Fill(entotHF);
00220     meEnergy_HO->Fill(entotHO);
00221 
00222     metime_HB->Fill(timetotHB);
00223     metime_HE->Fill(timetotHE);
00224     metime_HF->Fill(timetotHF);
00225     metime_HO->Fill(timetotHO);
00226     
00227     metime_enweighted_HB->Fill(timetotHB,entotHB);
00228     metime_enweighted_HE->Fill(timetotHE,entotHE);
00229     metime_enweighted_HF->Fill(timetotHF,entotHF);
00230     metime_enweighted_HO->Fill(timetotHO,entotHO);
00231   }
00232   
00233   for ( itr = map_try.begin() ; itr != map_try.end(); ++itr)   {
00234     if (dbe_ && (*itr).first.second >= 0) {
00235       HcalDetId id= (*itr).first.first;
00236       energysum ensum= (*itr).second;
00237       int eta = id.ieta();
00238       int phi = id.iphi();
00239       double etax= eta-0.5;
00240       double phix= phi-0.5;
00241       meHcalEnergyl25_[(*itr).first.second]->Fill(etax,phix,ensum.e25);
00242       meHcalEnergyl50_[(*itr).first.second]->Fill(etax,phix,ensum.e50);
00243       meHcalEnergyl100_[(*itr).first.second]->Fill(etax,phix,ensum.e100);
00244       meHcalEnergyl250_[(*itr).first.second]->Fill(etax,phix,ensum.e250);
00245 
00246       //LogDebug("HitsValidationHcal") <<" energy of tower ="<<(*itr).first<<"in time 25ns is== "<<(*itr).second.e25<<"in time 25-50ns=="<<(*itr).second.e50<<"in time 50-100ns=="<<(*itr).second.e100<<"in time 100-250 ns=="<<(*itr).second.e250;
00247     }
00248   }
00249   
00250 }
00251 
00252 DEFINE_FWK_MODULE(SimHitsValidationHcal);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Validation/HcalHits/src/SimHitsValidationHcal.cc
