CastorCalib/src/CastorLedAnalysis.cc

0001
0002 #include "CalibFormats/CastorObjects/interface/CastorCoderDb.h"
0003 #include "CalibFormats/CastorObjects/interface/CastorCalibrations.h"
0004 #include "CalibFormats/CastorObjects/interface/CastorDbService.h"
0005 #include "CalibFormats/CastorObjects/interface/CastorDbRecord.h"
0006 #include "CondFormats/CastorObjects/interface/CastorQIECoder.h"
0007 #include "CondFormats/CastorObjects/interface/CastorPedestals.h"
0008 #include "CondFormats/CastorObjects/interface/CastorPedestalWidths.h"
0009
0010 #include "CalibCalorimetry/CastorCalib/interface/CastorLedAnalysis.h"
0011 #include <TFile.h>
0012 #include <cmath>
0013
0014 using namespace std;
0015
0016 CastorLedAnalysis::CastorLedAnalysis(const edm::ParameterSet& ps) {
0017   // init
0018   evt = 0;
0019   sample = 0;
0020   m_file = nullptr;
0021   char output[100]{0};
0022   // output files
0023   for (int k = 0; k < 4; k++)
0024     state.push_back(true);  // 4 cap-ids (do we care?)
0025   m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
0026   m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML", "");
0027   if (!m_outputFileText.empty()) {
0028     cout << "Castor LED results will be saved to " << m_outputFileText.c_str() << endl;
0029     m_outFile.open(m_outputFileText.c_str());
0030   }
0031   m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
0032   if (!m_outputFileROOT.empty()) {
0033     cout << "Castor LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
0034   }
0035
0036   m_nevtsample = ps.getUntrackedParameter<int>("nevtsample", 9999999);
0037   if (m_nevtsample < 1)
0038     m_nevtsample = 9999999;
0039   m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag", 0);
0040   if (m_hiSaveflag < 0)
0041     m_hiSaveflag = 0;
0042   if (m_hiSaveflag > 0)
0043     m_hiSaveflag = 1;
0044   m_fitflag = ps.getUntrackedParameter<int>("analysisflag", 2);
0045   if (m_fitflag < 0)
0046     m_fitflag = 0;
0047   if (m_fitflag > 4)
0048     m_fitflag = 4;
0049   m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
0050   if (m_startTS < 0)
0051     m_startTS = 0;
0052   m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
0053   m_usecalib = ps.getUntrackedParameter<bool>("usecalib", false);
0054   m_logFile.open("CastorLedAnalysis.log");
0055
0056   int runNum = ps.getUntrackedParameter<int>("runNumber", 999999);
0057
0058   // histogram booking
0059   castorHists.ALLLEDS = new TH1F("Castor All LEDs", "HF All Leds", 10, 0, 9);
0060   castorHists.LEDRMS = new TH1F("Castor All LED RMS", "HF All LED RMS", 100, 0, 3);
0061   castorHists.LEDMEAN = new TH1F("Castor All LED Means", "HF All LED Means", 100, 0, 9);
0062   castorHists.CHI2 = new TH1F("Castor Chi2 by ndf for Landau fit", "HF Chi2/ndf Landau", 200, 0., 50.);
0063
0064   //XML file header
0065   m_outputFileXML.open(m_outputFileX.c_str());
0066   snprintf(output, sizeof output, "<?xml version='1.0' encoding='UTF-8'?>");
0067   m_outputFileXML << output << endl;
0068   snprintf(output, sizeof output, "<ROOT>");
0069   m_outputFileXML << output << endl << endl;
0070   snprintf(output, sizeof output, "  <HEADER>");
0071   m_outputFileXML << output << endl;
0072   snprintf(output, sizeof output, "    <TYPE>");
0073   m_outputFileXML << output << endl;
0074   snprintf(output, sizeof output, "      <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>");
0075   m_outputFileXML << output << endl;
0076   snprintf(output, sizeof output, "      <NAME>HCAL LED Timing</NAME>");
0077   m_outputFileXML << output << endl;
0078   snprintf(output, sizeof output, "    </TYPE>");
0079   m_outputFileXML << output << endl;
0080   snprintf(output, sizeof output, "    <RUN>");
0081   m_outputFileXML << output << endl;
0082   snprintf(output, sizeof output, "      <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>");
0083   m_outputFileXML << output << endl;
0084   snprintf(output, sizeof output, "      <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
0085   m_outputFileXML << output << endl;
0086   snprintf(output, sizeof output, "      <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>");
0087   m_outputFileXML << output << endl;
0088   snprintf(output, sizeof output, "      <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>");
0089   m_outputFileXML << output << endl;
0090   snprintf(output, sizeof output, "    </RUN>");
0091   m_outputFileXML << output << endl;
0092   snprintf(output, sizeof output, "  </HEADER>");
0093   m_outputFileXML << output << endl;
0094   snprintf(output, sizeof output, "<!-- Tags secton -->");
0095   m_outputFileXML << output << endl;
0096   snprintf(output, sizeof output, "  <ELEMENTS>");
0097   m_outputFileXML << output << endl;
0098   snprintf(output, sizeof output, "    <DATA_SET id='-1'/>");
0099   m_outputFileXML << output << endl;
0100   snprintf(output, sizeof output, "      <IOV id='1'>");
0101   m_outputFileXML << output << endl;
0102   snprintf(output, sizeof output, "        <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>");
0103   m_outputFileXML << output << endl;
0104   snprintf(output, sizeof output, "        <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>");
0105   m_outputFileXML << output << endl;
0106   snprintf(output, sizeof output, "      </IOV>");
0107   m_outputFileXML << output << endl;
0108   snprintf(output, sizeof output, "      <TAG id='2' mode='auto'>");
0109   m_outputFileXML << output << endl;
0110   snprintf(output, sizeof output, "        <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
0111   m_outputFileXML << output << endl;
0112   snprintf(output, sizeof output, "        <DETECTOR_NAME>HCAL</DETECTOR_NAME>");
0113   m_outputFileXML << output << endl;
0114   snprintf(output, sizeof output, "        <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>");
0115   m_outputFileXML << output << endl;
0116   snprintf(output, sizeof output, "      </TAG>");
0117   m_outputFileXML << output << endl;
0118   snprintf(output, sizeof output, "  </ELEMENTS>");
0119   m_outputFileXML << output << endl;
0120   snprintf(output, sizeof output, "  <MAPS>");
0121   m_outputFileXML << output << endl;
0122   snprintf(output, sizeof output, "      <TAG idref ='2'>");
0123   m_outputFileXML << output << endl;
0124   snprintf(output, sizeof output, "        <IOV idref='1'>");
0125   m_outputFileXML << output << endl;
0126   snprintf(output, sizeof output, "          <DATA_SET idref='-1' />");
0127   m_outputFileXML << output << endl;
0128   snprintf(output, sizeof output, "        </IOV>");
0129   m_outputFileXML << output << endl;
0130   snprintf(output, sizeof output, "      </TAG>");
0131   m_outputFileXML << output << endl;
0132   snprintf(output, sizeof output, "  </MAPS>");
0133   m_outputFileXML << output << endl;
0134 }
0135
0136 //-----------------------------------------------------------------------------
0137 CastorLedAnalysis::~CastorLedAnalysis() {
0138   ///All done, clean up!!
0139
0140   for (_meol = castorHists.LEDTRENDS.begin(); _meol != castorHists.LEDTRENDS.end(); _meol++) {
0141     for (int i = 0; i < 15; i++)
0142       _meol->second[i].first->Delete();
0143   }
0144
0145   castorHists.ALLLEDS->Delete();
0146   castorHists.LEDRMS->Delete();
0147   castorHists.LEDMEAN->Delete();
0148   castorHists.CHI2->Delete();
0149 }
0150
0151 //-----------------------------------------------------------------------------
0152 void CastorLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
0153   // open the histogram file, create directories within
0154   m_file = new TFile(m_outputFileROOT.c_str(), "RECREATE");
0155   m_file->mkdir("Castor");
0156   m_file->cd();
0157 }
0158
0159 //-----------------------------------------------------------------------------
0160 void CastorLedAnalysis::GetLedConst(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0161   double time2 = 0;
0162   double time1 = 0;
0163   double time3 = 0;
0164   double time4 = 0;
0165   double dtime2 = 0;
0166   double dtime1 = 0;
0167   double dtime3 = 0;
0168   double dtime4 = 0;
0169   char output[100]{0};
0170
0171   if (!m_outputFileText.empty()) {
0172     if (m_fitflag == 0 || m_fitflag == 2)
0173       m_outFile << "Det Eta,Phi,D   Mean    Error" << std::endl;
0174     else if (m_fitflag == 1 || m_fitflag == 3)
0175       m_outFile << "Det Eta,Phi,D   Peak    Error" << std::endl;
0176     else if (m_fitflag == 4)
0177       m_outFile << "Det Eta,Phi,D   Mean    Error      Peak    Error       MeanEv  Error       PeakEv  Error"
0178                 << std::endl;
0179   }
0180   for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
0181     // scale the LED pulse to 1 event
0182     _meol->second[10].first->Scale(1. / evt_curr);
0183     if (m_fitflag == 0 || m_fitflag == 4) {
0184       time1 = _meol->second[10].first->GetMean();
0185       dtime1 = _meol->second[10].first->GetRMS() / sqrt((float)evt_curr * (m_endTS - m_startTS + 1));
0186     }
0187     if (m_fitflag == 1 || m_fitflag == 4) {
0188       // put proper errors
0189       for (int j = 0; j < 10; j++)
0190         _meol->second[10].first->SetBinError(j + 1, _meol->second[j].first->GetRMS() / sqrt((float)evt_curr));
0191     }
0192     if (m_fitflag == 1 || m_fitflag == 3 || m_fitflag == 4) {
0193       _meol->second[10].first->Fit("landau", "Q");
0194       //      _meol->second[10].first->Fit("gaus","Q");
0195       TF1* fit = _meol->second[10].first->GetFunction("landau");
0196       //      TF1 *fit = _meol->second[10].first->GetFunction("gaus");
0197       time2 = fit->GetParameter(1);
0198       dtime2 = fit->GetParError(1);
0199     }
0200     if (m_fitflag == 2 || m_fitflag == 4) {
0201       time3 = _meol->second[12].first->GetMean();
0202       dtime3 = _meol->second[12].first->GetRMS() / sqrt((float)_meol->second[12].first->GetEntries());
0203     }
0204     if (m_fitflag == 3 || m_fitflag == 4) {
0205       time4 = _meol->second[13].first->GetMean();
0206       dtime4 = _meol->second[13].first->GetRMS() / sqrt((float)_meol->second[13].first->GetEntries());
0207     }
0208     for (int i = 0; i < 10; i++) {
0209       _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
0210       _meol->second[i].first->GetYaxis()->SetTitle("Counts");
0211       //      if(m_hiSaveflag>0)_meol->second[i].first->Write();
0212     }
0213     _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
0214     _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
0215     if (m_hiSaveflag > 0)
0216       _meol->second[10].first->Write();
0217     _meol->second[10].second.first[0].push_back(time1);
0218     _meol->second[10].second.first[1].push_back(dtime1);
0219     _meol->second[11].second.first[0].push_back(time2);
0220     _meol->second[11].second.first[1].push_back(dtime2);
0221     _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
0222     _meol->second[12].first->GetYaxis()->SetTitle("Counts");
0223     if (m_fitflag == 2 && m_hiSaveflag > 0)
0224       _meol->second[12].first->Write();
0225     _meol->second[12].second.first[0].push_back(time3);
0226     _meol->second[12].second.first[1].push_back(dtime3);
0227     _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
0228     _meol->second[13].first->GetYaxis()->SetTitle("Counts");
0229     if (m_fitflag > 2 && m_hiSaveflag > 0)
0230       _meol->second[13].first->Write();
0231     _meol->second[13].second.first[0].push_back(time4);
0232     _meol->second[13].second.first[1].push_back(dtime4);
0233     _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
0234     _meol->second[14].first->GetYaxis()->SetTitle("Counts");
0235     if (m_fitflag > 2 && m_hiSaveflag > 0)
0236       _meol->second[14].first->Write();
0237     _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
0238     _meol->second[15].first->GetYaxis()->SetTitle("Counts");
0239     if (m_fitflag > 2 && m_hiSaveflag > 0)
0240       _meol->second[15].first->Write();
0241     _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
0242     _meol->second[16].first->Write();
0243
0244     // Ascii printout (need to modify to include new info)
0245     HcalDetId detid = _meol->first;
0246
0247     if (!m_outputFileText.empty()) {
0248       if (m_fitflag == 0) {
0249         m_outFile << detid << "   " << time1 << " " << dtime1 << std::endl;
0250         snprintf(output, sizeof output, "  <DATA_SET>");
0251         m_outputFileXML << output << endl;
0252         snprintf(output, sizeof output, "    <VERSION>version:1</VERSION>");
0253         m_outputFileXML << output << endl;
0254         snprintf(output, sizeof output, "    <CHANNEL>");
0255         m_outputFileXML << output << endl;
0256         snprintf(output, sizeof output, "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
0257         m_outputFileXML << output << endl;
0258         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0259         m_outputFileXML << output << endl;
0260         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0261         m_outputFileXML << output << endl;
0262         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0263         m_outputFileXML << output << endl;
0264         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0265         m_outputFileXML << output << endl;
0266         if (detid.subdet() == 1)
0267           snprintf(output, sizeof output, "      <DETECTOR_NAME>HB</DETECTOR_NAME>");
0268         if (detid.subdet() == 2)
0269           snprintf(output, sizeof output, "      <DETECTOR_NAME>HE</DETECTOR_NAME>");
0270         if (detid.subdet() == 3)
0271           snprintf(output, sizeof output, "      <DETECTOR_NAME>HO</DETECTOR_NAME>");
0272         if (detid.subdet() == 4)
0273           snprintf(output, sizeof output, "      <DETECTOR_NAME>HF</DETECTOR_NAME>");
0274         m_outputFileXML << output << endl;
0275         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0276         m_outputFileXML << output << endl;
0277         snprintf(output, sizeof output, "    </CHANNEL>");
0278         m_outputFileXML << output << endl;
0279         snprintf(output, sizeof output, "    <DATA>");
0280         m_outputFileXML << output << endl;
0281         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time1);
0282         m_outputFileXML << output << endl;
0283         snprintf(output, sizeof output, "      <OFFSET_TIME> 0</OFFSET_TIME>");
0284         m_outputFileXML << output << endl;
0285         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
0286         m_outputFileXML << output << endl;
0287         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0288         m_outputFileXML << output << endl;
0289         snprintf(output, sizeof output, "      <STATUS_WORD>  0</STATUS_WORD>");
0290         m_outputFileXML << output << endl;
0291         snprintf(output, sizeof output, "    </DATA>");
0292         m_outputFileXML << output << endl;
0293         snprintf(output, sizeof output, "  </DATA_SET>");
0294         m_outputFileXML << output << endl;
0295
0296       } else if (m_fitflag == 1) {
0297         m_outFile << detid << "   " << time2 << " " << dtime2 << std::endl;
0298         snprintf(output, sizeof output, "  <DATA_SET>");
0299         m_outputFileXML << output << endl;
0300         snprintf(output, sizeof output, "    <VERSION>version:1</VERSION>");
0301         m_outputFileXML << output << endl;
0302         snprintf(output, sizeof output, "    <CHANNEL>");
0303         m_outputFileXML << output << endl;
0304         snprintf(output, sizeof output, "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
0305         m_outputFileXML << output << endl;
0306         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0307         m_outputFileXML << output << endl;
0308         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0309         m_outputFileXML << output << endl;
0310         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0311         m_outputFileXML << output << endl;
0312         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0313         m_outputFileXML << output << endl;
0314         if (detid.subdet() == 1)
0315           snprintf(output, sizeof output, "      <DETECTOR_NAME>HB</DETECTOR_NAME>");
0316         if (detid.subdet() == 2)
0317           snprintf(output, sizeof output, "      <DETECTOR_NAME>HE</DETECTOR_NAME>");
0318         if (detid.subdet() == 3)
0319           snprintf(output, sizeof output, "      <DETECTOR_NAME>HO</DETECTOR_NAME>");
0320         if (detid.subdet() == 4)
0321           snprintf(output, sizeof output, "      <DETECTOR_NAME>HF</DETECTOR_NAME>");
0322         m_outputFileXML << output << endl;
0323         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0324         m_outputFileXML << output << endl;
0325         snprintf(output, sizeof output, "    </CHANNEL>");
0326         m_outputFileXML << output << endl;
0327         snprintf(output, sizeof output, "    <DATA>");
0328         m_outputFileXML << output << endl;
0329         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time2);
0330         m_outputFileXML << output << endl;
0331         snprintf(output, sizeof output, "      <OFFSET_TIME> 0</OFFSET_TIME>");
0332         m_outputFileXML << output << endl;
0333         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime2);
0334         m_outputFileXML << output << endl;
0335         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0336         m_outputFileXML << output << endl;
0337         snprintf(output, sizeof output, "      <STATUS_WORD>  0</STATUS_WORD>");
0338         m_outputFileXML << output << endl;
0339         snprintf(output, sizeof output, "    </DATA>");
0340         m_outputFileXML << output << endl;
0341         snprintf(output, sizeof output, "  </DATA_SET>");
0342         m_outputFileXML << output << endl;
0343       }
0344
0345       else if (m_fitflag == 2) {
0346         m_outFile << detid << "   " << time3 << " " << dtime3 << std::endl;
0347         snprintf(output, sizeof output, "  <DATA_SET>");
0348         m_outputFileXML << output << endl;
0349         snprintf(output, sizeof output, "    <VERSION>version:1</VERSION>");
0350         m_outputFileXML << output << endl;
0351         snprintf(output, sizeof output, "    <CHANNEL>");
0352         m_outputFileXML << output << endl;
0353         snprintf(output, sizeof output, "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
0354         m_outputFileXML << output << endl;
0355         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0356         m_outputFileXML << output << endl;
0357         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0358         m_outputFileXML << output << endl;
0359         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0360         m_outputFileXML << output << endl;
0361         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0362         m_outputFileXML << output << endl;
0363         if (detid.subdet() == 1)
0364           snprintf(output, sizeof output, "      <DETECTOR_NAME>HB</DETECTOR_NAME>");
0365         if (detid.subdet() == 2)
0366           snprintf(output, sizeof output, "      <DETECTOR_NAME>HE</DETECTOR_NAME>");
0367         if (detid.subdet() == 3)
0368           snprintf(output, sizeof output, "      <DETECTOR_NAME>HO</DETECTOR_NAME>");
0369         if (detid.subdet() == 4)
0370           snprintf(output, sizeof output, "      <DETECTOR_NAME>HF</DETECTOR_NAME>");
0371         m_outputFileXML << output << endl;
0372         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0373         m_outputFileXML << output << endl;
0374         snprintf(output, sizeof output, "    </CHANNEL>");
0375         m_outputFileXML << output << endl;
0376         snprintf(output, sizeof output, "    <DATA>");
0377         m_outputFileXML << output << endl;
0378         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time3);
0379         m_outputFileXML << output << endl;
0380         snprintf(output, sizeof output, "      <OFFSET_TIME> 0</OFFSET_TIME>");
0381         m_outputFileXML << output << endl;
0382         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime3);
0383         m_outputFileXML << output << endl;
0384         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0385         m_outputFileXML << output << endl;
0386         snprintf(output, sizeof output, "      <STATUS_WORD>  0</STATUS_WORD>");
0387         m_outputFileXML << output << endl;
0388         snprintf(output, sizeof output, "    </DATA>");
0389         m_outputFileXML << output << endl;
0390         snprintf(output, sizeof output, "  </DATA_SET>");
0391         m_outputFileXML << output << endl;
0392       } else if (m_fitflag == 3) {
0393         m_outFile << detid << "   " << time4 << " " << dtime4 << std::endl;
0394         snprintf(output, sizeof output, "  <DATA_SET>");
0395         m_outputFileXML << output << endl;
0396         snprintf(output, sizeof output, "    <VERSION>version:1</VERSION>");
0397         m_outputFileXML << output << endl;
0398         snprintf(output, sizeof output, "    <CHANNEL>");
0399         m_outputFileXML << output << endl;
0400         snprintf(output, sizeof output, "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
0401         m_outputFileXML << output << endl;
0402         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0403         m_outputFileXML << output << endl;
0404         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0405         m_outputFileXML << output << endl;
0406         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0407         m_outputFileXML << output << endl;
0408         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0409         m_outputFileXML << output << endl;
0410         if (detid.subdet() == 1)
0411           snprintf(output, sizeof output, "      <DETECTOR_NAME>HB</DETECTOR_NAME>");
0412         if (detid.subdet() == 2)
0413           snprintf(output, sizeof output, "      <DETECTOR_NAME>HE</DETECTOR_NAME>");
0414         if (detid.subdet() == 3)
0415           snprintf(output, sizeof output, "      <DETECTOR_NAME>HO</DETECTOR_NAME>");
0416         if (detid.subdet() == 4)
0417           snprintf(output, sizeof output, "      <DETECTOR_NAME>HF</DETECTOR_NAME>");
0418         m_outputFileXML << output << endl;
0419         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0420         m_outputFileXML << output << endl;
0421         snprintf(output, sizeof output, "    </CHANNEL>");
0422         m_outputFileXML << output << endl;
0423         snprintf(output, sizeof output, "    <DATA>");
0424         m_outputFileXML << output << endl;
0425         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time4);
0426         m_outputFileXML << output << endl;
0427         snprintf(output, sizeof output, "      <OFFSET_TIME> 0</OFFSET_TIME>");
0428         m_outputFileXML << output << endl;
0429         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
0430         m_outputFileXML << output << endl;
0431         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0432         m_outputFileXML << output << endl;
0433         snprintf(output, sizeof output, "      <STATUS_WORD>  0</STATUS_WORD>");
0434         m_outputFileXML << output << endl;
0435         snprintf(output, sizeof output, "    </DATA>");
0436         m_outputFileXML << output << endl;
0437         snprintf(output, sizeof output, "  </DATA_SET>");
0438         m_outputFileXML << output << endl;
0439       }
0440
0441       else if (m_fitflag == 4) {
0442         m_outFile << detid << "   " << time1 << " " << dtime1 << "   " << time2 << " " << dtime2 << "   " << time3
0443                   << " " << dtime3 << "   " << time4 << " " << dtime4 << std::endl;
0444       }
0445     }
0446   }
0447 }
0448
0449 //-----------------------------------------------------------------------------
0450 void CastorLedAnalysis::LedSampleAnalysis() {
0451   // it is called every m_nevtsample events (a sample) and the end of run
0452   char LedSampleNum[20];
0453
0454   sprintf(LedSampleNum, "LedSample_%d", sample);
0455   m_file->cd();
0456   m_file->mkdir(LedSampleNum);
0457   m_file->cd(LedSampleNum);
0458
0459   // Compute LED constants
0460   GetLedConst(castorHists.LEDTRENDS);
0461 }
0462
0463 //-----------------------------------------------------------------------------
0464 void CastorLedAnalysis::LedTrendings(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0465   for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
0466     char name[1024];
0467     HcalDetId detid = _meol->first;
0468     sprintf(name, "LED timing trend, eta=%d phi=%d depth=%d", detid.ieta(), detid.iphi(), detid.depth());
0469     int bins = _meol->second[10 + m_fitflag].second.first[0].size();
0470     float lo = 0.5;
0471     float hi = (float)bins + 0.5;
0472     _meol->second[10 + m_fitflag].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0473
0474     std::vector<double>::iterator sample_it;
0475     // LED timing - put content and errors
0476     int j = 0;
0477     for (sample_it = _meol->second[10 + m_fitflag].second.first[0].begin();
0478          sample_it != _meol->second[10 + m_fitflag].second.first[0].end();
0479          ++sample_it) {
0480       _meol->second[10 + m_fitflag].second.second[0]->SetBinContent(++j, *sample_it);
0481     }
0482     j = 0;
0483     for (sample_it = _meol->second[10 + m_fitflag].second.first[1].begin();
0484          sample_it != _meol->second[10 + m_fitflag].second.first[1].end();
0485          ++sample_it) {
0486       _meol->second[10 + m_fitflag].second.second[0]->SetBinError(++j, *sample_it);
0487     }
0488     sprintf(name, "Sample (%d events)", m_nevtsample);
0489     _meol->second[10 + m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
0490     _meol->second[10 + m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
0491     _meol->second[10 + m_fitflag].second.second[0]->Write();
0492   }
0493 }
0494
0495 //-----------------------------------------------------------------------------
0496 void CastorLedAnalysis::LedDone() {
0497   // First process the last sample (remaining events).
0498   if (evt % m_nevtsample != 0)
0499     LedSampleAnalysis();
0500
0501   // Now do the end of run analysis: trending histos
0502   if (sample > 1 && m_fitflag != 4) {
0503     m_file->cd();
0504     m_file->cd("Castor");
0505     LedTrendings(castorHists.LEDTRENDS);
0506   }
0507
0508   // Write other histograms.
0509   m_file->cd();
0510   m_file->cd("Castor");
0511   castorHists.ALLLEDS->Write();
0512   castorHists.LEDRMS->Write();
0513   castorHists.LEDMEAN->Write();
0514
0515   // Write the histo file and close it
0516   //  m_file->Write();
0517   m_file->Close();
0518   cout << "Castor histograms written to " << m_outputFileROOT.c_str() << endl;
0519 }
0520
0521 //-----------------------------------------------------------------------------
0522 void CastorLedAnalysis::processLedEvent(const CastorDigiCollection& castor, const CastorDbService& cond) {
0523   evt++;
0524   sample = (evt - 1) / m_nevtsample + 1;
0525   evt_curr = evt % m_nevtsample;
0526   if (evt_curr == 0)
0527     evt_curr = m_nevtsample;
0528
0529   // HF/Castor
0530   if (castor.empty()) {
0531     edm::LogError("CastorLedAnalysis") << "Event with " << (int)castor.size() << "Castor Digis passed." << std::endl;
0532     return;
0533   }
0534   for (CastorDigiCollection::const_iterator j = castor.begin(); j != castor.end(); ++j) {
0535     const CastorDataFrame digi = (const CastorDataFrame)(*j);
0536     _meol = castorHists.LEDTRENDS.find(digi.id());
0537     if (_meol == castorHists.LEDTRENDS.end()) {
0538       SetupLEDHists(2, digi.id(), castorHists.LEDTRENDS);
0539     }
0540     LedCastorHists(digi.id(), digi, castorHists.LEDTRENDS, cond);
0541   }
0542
0543   // Call the function every m_nevtsample events
0544   if (evt % m_nevtsample == 0)
0545     LedSampleAnalysis();
0546 }
0547 //----------------------------------------------------------------------------
0548 void CastorLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0549   string type = "HBHE";
0550   if (id == 1)
0551     type = "HO";
0552   if (id == 2)
0553     type = "HF";
0554
0555   _meol = toolT.find(detid);
0556   if (_meol == toolT.end()) {
0557     // if histos for this channel do not exist, create them
0558     map<int, LEDBUNCH> insert;
0559     char name[1024];
0560     for (int i = 0; i < 10; i++) {
0561       sprintf(name,
0562               "%s Pulse height, eta=%d phi=%d depth=%d TS=%d",
0563               type.c_str(),
0564               detid.ieta(),
0565               detid.iphi(),
0566               detid.depth(),
0567               i);
0568       insert[i].first = new TH1F(name, name, 200, 0., 2000.);
0569     }
0570     sprintf(name, "%s LED Mean pulse, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0571     insert[10].first = new TH1F(name, name, 10, -0.5, 9.5);
0572     sprintf(name, "%s LED Pulse, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0573     insert[11].first = new TH1F(name, name, 10, -0.5, 9.5);
0574     sprintf(name, "%s Mean TS, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0575     insert[12].first = new TH1F(name, name, 200, 0., 10.);
0576     sprintf(name, "%s Peak TS, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0577     insert[13].first = new TH1F(name, name, 200, 0., 10.);
0578     sprintf(name, "%s Peak TS error, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0579     insert[14].first = new TH1F(name, name, 200, 0., 0.05);
0580     sprintf(name, "%s Fit chi2, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0581     insert[15].first = new TH1F(name, name, 100, 0., 50.);
0582     sprintf(
0583         name, "%s Integrated Signal, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0584     insert[16].first = new TH1F(name, name, 500, 0., 5000.);
0585
0586     toolT[detid] = insert;
0587   }
0588 }
0589 //-----------------------------------------------------------------------------
0590 void CastorLedAnalysis::LedCastorHists(const HcalDetId& detid,
0591                                        const CastorDataFrame& ledDigi,
0592                                        map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0593                                        const CastorDbService& cond) {
0594   map<int, LEDBUNCH> _mei;
0595   _meol = toolT.find(detid);
0596   _mei = _meol->second;
0597   // Rest the histos if we're at the end of a 'bunch'
0598   if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0599     for (int k = 0; k < (int)state.size(); k++)
0600       state[k] = false;
0601     for (int i = 0; i < 16; i++)
0602       _mei[i].first->Reset();
0603   }
0604
0605   // now we have the signal in fC, let's get rid of that darn pedestal
0606   // Most of this is borrowed from CastorSimpleReconstructor, so thanks Jeremy/Phil
0607
0608   float max_fC = 0;
0609   float ta = 0;
0610   m_coder = cond.getCastorCoder(detid);
0611   m_ped = cond.getPedestal(detid);
0612   m_shape = cond.getCastorShape();
0613   //cout << "New Digi!!!!!!!!!!!!!!!!!!!!!!" << endl;
0614   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0615     int capid = ledDigi[TS].capid();
0616     // BE CAREFUL: this is assuming peds are stored in ADCs
0617     int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
0618     if (adc < 0) {
0619       adc = 0;
0620     }  // to prevent negative adcs after ped subtraction, which should really only happen
0621        // if you're using the wrong peds.
0622     double fC = m_coder->charge(*m_shape, adc, capid);
0623     //ta = (fC - m_ped->getValue(capid));
0624     ta = fC;
0625     //cout << "DetID: " << detid << "  CapID: " << capid << "  ADC: " << adc << "  Ped: " << m_ped->getValue(capid) << "  fC: " << fC << endl;
0626     _mei[TS].first->Fill(ta);
0627     _mei[10].first->AddBinContent(TS + 1, ta);  // This is average pulse, could probably do better (Profile?)
0628     if (m_fitflag > 1) {
0629       if (TS == m_startTS)
0630         _mei[11].first->Reset();
0631       _mei[11].first->SetBinContent(TS + 1, ta);
0632     }
0633
0634     // keep track of max TS and max amplitude (in fC)
0635     if (ta > max_fC) {
0636       max_fC = ta;
0637     }
0638   }
0639
0640   // Now we have a sample with pedestals subtracted and in units of fC
0641   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
0642   // we now want to use Phil's timing correction.  This is not necessary
0643   // if we are performing a Landau fit (m_fitflag = 3)
0644
0645   float sum = 0.;
0646   for (int i = 0; i < 10; i++)
0647     sum = sum + _mei[11].first->GetBinContent(i + 1);
0648   if (sum > 100) {
0649     if (m_fitflag == 2 || m_fitflag == 4) {
0650       float timmean = _mei[11].first->GetMean();  // let's use Phil's way instead
0651       float timmeancorr = BinsizeCorr(timmean);
0652       _mei[12].first->Fill(timmeancorr);
0653     }
0654     _mei[16].first->Fill(
0655         _mei[11].first->Integral());  // Integrated charge (may be more usfull to convert to Energy first?)
0656     if (m_fitflag == 3 || m_fitflag == 4) {
0657       _mei[11].first->Fit("landau", "Q");
0658       TF1* fit = _mei[11].first->GetFunction("landau");
0659       _mei[13].first->Fill(fit->GetParameter(1));
0660       _mei[14].first->Fill(fit->GetParError(1));
0661       _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0662     }
0663   }
0664 }
0665
0666 //-----------------------------------------------------------------------------
0667 float CastorLedAnalysis::BinsizeCorr(float time) {
0668   // this is the bin size correction to be applied for laser data (from Andy),
0669   // it comes from a pulse shape measured from TB04 data (from Jordan)
0670   // This should eventually be replaced with the more thorough treatment from Phil
0671
0672   float corrtime = 0.;
0673   static const float tstrue[32] = {0.003, 0.03425, 0.06548, 0.09675, 0.128, 0.15925, 0.1905, 0.22175,
0674                                    0.253, 0.28425, 0.3155,  0.34675, 0.378, 0.40925, 0.4405, 0.47175,
0675                                    0.503, 0.53425, 0.5655,  0.59675, 0.628, 0.65925, 0.6905, 0.72175,
0676                                    0.753, 0.78425, 0.8155,  0.84675, 0.878, 0.90925, 0.9405, 0.97175};
0677   static const float tsreco[32] = {-0.00422, 0.01815, 0.04409, 0.07346, 0.09799, 0.12192, 0.15072, 0.18158,
0678                                    0.21397,  0.24865, 0.28448, 0.31973, 0.35449, 0.39208, 0.43282, 0.47244,
0679                                    0.5105,   0.55008, 0.58827, 0.62828, 0.6717,  0.70966, 0.74086, 0.77496,
0680                                    0.80843,  0.83472, 0.86044, 0.8843,  0.90674, 0.92982, 0.95072, 0.9726};
0681
0682   int inttime = (int)time;
0683   float restime = time - inttime;
0684   for (int i = 0; i <= 32; i++) {
0685     float lolim = 0.;
0686     float uplim = 1.;
0687     float tsdown;
0688     float tsup;
0689     if (i > 0) {
0690       lolim = tsreco[i - 1];
0691       tsdown = tstrue[i - 1];
0692     } else
0693       tsdown = tstrue[31] - 1.;
0694     if (i < 32) {
0695       uplim = tsreco[i];
0696       tsup = tstrue[i];
0697     } else
0698       tsup = tstrue[0] + 1.;
0699     if (restime >= lolim && restime < uplim) {
0700       corrtime = (tsdown * (uplim - restime) + tsup * (restime - lolim)) / (uplim - lolim);
0701     }
0702   }
0703   corrtime += inttime;
0704
0705   return corrtime;
0706 }
0707 //-----------------------------------------------------------------------------