Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:34:53 

0001 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0002 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
0003 #include "CondFormats/HcalObjects/interface/HcalPedestals.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "CalibCalorimetry/HcalAlgos/interface/HcalLedAnalysis.h"
0006 #include "TFile.h"
0007 #include <cmath>
0008 using namespace std;
0009 
0010 HcalLedAnalysis::HcalLedAnalysis(const edm::ParameterSet& ps) {
0011   // init
0012 
0013   m_coder = nullptr;
0014   m_ped = nullptr;
0015   m_shape = nullptr;
0016   evt = 0;
0017   sample = 0;
0018   m_file = nullptr;
0019   // output files
0020   for (int k = 0; k < 4; k++)
0021     state.push_back(true);  // 4 cap-ids (do we care?)
0022   m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
0023   m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML", "");
0024   if (!m_outputFileText.empty()) {
0025     edm::LogInfo("HcalLedAnalysis") << "Hcal LED results will be saved to " << m_outputFileText.c_str() << endl;
0026     m_outFile.open(m_outputFileText.c_str());
0027   }
0028   m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
0029   if (!m_outputFileROOT.empty()) {
0030     edm::LogInfo("HcalLedAnalysis") << "Hcal LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
0031   }
0032 
0033   m_nevtsample = ps.getUntrackedParameter<int>("nevtsample", 9999999);
0034   if (m_nevtsample < 1)
0035     m_nevtsample = 9999999;
0036   m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag", 0);
0037   if (m_hiSaveflag < 0)
0038     m_hiSaveflag = 0;
0039   if (m_hiSaveflag > 0)
0040     m_hiSaveflag = 1;
0041   m_fitflag = ps.getUntrackedParameter<int>("analysisflag", 2);
0042   if (m_fitflag < 0)
0043     m_fitflag = 0;
0044   if (m_fitflag > 4)
0045     m_fitflag = 4;
0046   m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
0047   if (m_startTS < 0)
0048     m_startTS = 0;
0049   m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
0050   m_usecalib = ps.getUntrackedParameter<bool>("usecalib", false);
0051 
0052   int runNum = ps.getUntrackedParameter<int>("runNumber", 999999);
0053 
0054   // histogram booking
0055   hbHists.ALLLEDS = new TH1F("HBHE All LEDs", "HB/HE All Leds", 10, 0, 9);
0056   hbHists.LEDRMS = new TH1F("HBHE All LED RMS", "HB/HE All LED RMS", 100, 0, 3);
0057   hbHists.LEDMEAN = new TH1F("HBHE All LED Means", "HB/HE All LED Means", 100, 0, 9);
0058   hbHists.CHI2 = new TH1F("HBHE Chi2 by ndf for Landau fit", "HB/HE Chi2/ndf Landau", 200, 0., 50.);
0059 
0060   hoHists.ALLLEDS = new TH1F("HO All LEDs", "HO All Leds", 10, 0, 9);
0061   hoHists.LEDRMS = new TH1F("HO All LED RMS", "HO All LED RMS", 100, 0, 3);
0062   hoHists.LEDMEAN = new TH1F("HO All LED Means", "HO All LED Means", 100, 0, 9);
0063   hoHists.CHI2 = new TH1F("HO Chi2 by ndf for Landau fit", "HO Chi2/ndf Landau", 200, 0., 50.);
0064 
0065   hfHists.ALLLEDS = new TH1F("HF All LEDs", "HF All Leds", 10, 0, 9);
0066   hfHists.LEDRMS = new TH1F("HF All LED RMS", "HF All LED RMS", 100, 0, 3);
0067   hfHists.LEDMEAN = new TH1F("HF All LED Means", "HF All LED Means", 100, 0, 9);
0068   hfHists.CHI2 = new TH1F("HF Chi2 by ndf for Landau fit", "HF Chi2/ndf Landau", 200, 0., 50.);
0069 
0070   char output[100]{0};
0071 
0072   //XML file header
0073   m_outputFileXML.open(m_outputFileX.c_str());
0074 
0075   m_outputFileXML << "<?xml version='1.0' encoding='UTF-8'?>" << endl;
0076 
0077   m_outputFileXML << "<ROOT>" << endl;
0078 
0079   m_outputFileXML << "  <HEADER>" << endl;
0080 
0081   m_outputFileXML << "    <TYPE>" << endl;
0082 
0083   m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>" << endl;
0084 
0085   m_outputFileXML << "      <NAME>HCAL LED Timing</NAME>" << endl;
0086 
0087   m_outputFileXML << "    </TYPE>" << endl;
0088 
0089   m_outputFileXML << "    <RUN>" << endl;
0090 
0091   m_outputFileXML << "      <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>" << endl;
0092 
0093   snprintf(output, sizeof output, "      <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
0094   m_outputFileXML << output << endl;
0095 
0096   m_outputFileXML << "      <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>" << endl;
0097 
0098   m_outputFileXML << "      <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
0099 
0100   m_outputFileXML << "    </RUN>" << endl;
0101 
0102   m_outputFileXML << "  </HEADER>" << endl;
0103 
0104   m_outputFileXML << "<!-- Tags secton -->" << endl;
0105 
0106   m_outputFileXML << "  <ELEMENTS>" << endl;
0107 
0108   m_outputFileXML << "    <DATA_SET id='-1'/>" << endl;
0109 
0110   m_outputFileXML << "      <IOV id='1'>" << endl;
0111 
0112   m_outputFileXML << "        <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>" << endl;
0113 
0114   m_outputFileXML << "        <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>" << endl;
0115 
0116   m_outputFileXML << "      </IOV>" << endl;
0117 
0118   m_outputFileXML << "      <TAG id='2' mode='auto'>" << endl;
0119 
0120   snprintf(output, sizeof output, "        <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
0121   m_outputFileXML << output << endl;
0122 
0123   m_outputFileXML << "        <DETECTOR_NAME>HCAL</DETECTOR_NAME>" << endl;
0124 
0125   m_outputFileXML << "        <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
0126 
0127   m_outputFileXML << "      </TAG>" << endl;
0128 
0129   m_outputFileXML << "  </ELEMENTS>" << endl;
0130 
0131   m_outputFileXML << "  <MAPS>" << endl;
0132 
0133   m_outputFileXML << "      <TAG idref ='2'>" << endl;
0134 
0135   m_outputFileXML << "        <IOV idref='1'>" << endl;
0136 
0137   m_outputFileXML << "          <DATA_SET idref='-1' />" << endl;
0138 
0139   m_outputFileXML << "        </IOV>" << endl;
0140 
0141   m_outputFileXML << "      </TAG>" << endl;
0142 
0143   m_outputFileXML << "  </MAPS>" << endl;
0144 }
0145 
0146 //-----------------------------------------------------------------------------
0147 HcalLedAnalysis::~HcalLedAnalysis() {
0148   ///All done, clean up!!
0149   for (_meol = hbHists.LEDTRENDS.begin(); _meol != hbHists.LEDTRENDS.end(); _meol++) {
0150     for (int i = 0; i < 15; i++)
0151       _meol->second[i].first->Delete();
0152   }
0153   for (_meol = hoHists.LEDTRENDS.begin(); _meol != hoHists.LEDTRENDS.end(); _meol++) {
0154     for (int i = 0; i < 15; i++)
0155       _meol->second[i].first->Delete();
0156   }
0157   for (_meol = hfHists.LEDTRENDS.begin(); _meol != hfHists.LEDTRENDS.end(); _meol++) {
0158     for (int i = 0; i < 15; i++)
0159       _meol->second[i].first->Delete();
0160   }
0161   hbHists.ALLLEDS->Delete();
0162   hbHists.LEDRMS->Delete();
0163   hbHists.LEDMEAN->Delete();
0164   hbHists.CHI2->Delete();
0165 
0166   hoHists.ALLLEDS->Delete();
0167   hoHists.LEDRMS->Delete();
0168   hoHists.LEDMEAN->Delete();
0169   hoHists.CHI2->Delete();
0170 
0171   hfHists.ALLLEDS->Delete();
0172   hfHists.LEDRMS->Delete();
0173   hfHists.LEDMEAN->Delete();
0174   hfHists.CHI2->Delete();
0175 }
0176 
0177 //-----------------------------------------------------------------------------
0178 void HcalLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
0179   // open the histogram file, create directories within
0180   m_file = new TFile(m_outputFileROOT.c_str(), "RECREATE");
0181   m_file->mkdir("HBHE");
0182   m_file->cd();
0183   m_file->mkdir("HO");
0184   m_file->cd();
0185   m_file->mkdir("HF");
0186   m_file->cd();
0187   m_file->mkdir("Calib");
0188   m_file->cd();
0189 }
0190 
0191 //-----------------------------------------------------------------------------
0192 /*
0193 void HcalLedAnalysis::doPeds(const HcalPedestal* fInputPedestals){
0194 // put all pedestals in a map m_AllPedVals, to be used in processLedEvent -
0195 // sorry, this is the only way I was able to implement pedestal subtraction
0196 
0197 // DEPRECATED
0198 // This is no longer useful, better ways of doing it -A
0199   HcalDetId detid;
0200   map<int,float> PedVals;
0201   pedCan = fInputPedestals;
0202   if(pedCan){
0203     std::vector<DetId> Channs=pedCan->getAllChannels();
0204     for (int i=0; i<(int)Channs.size(); i++){
0205       detid=HcalDetId (Channs[i]);
0206       for (int icap=0; icap<4; icap++) PedVals[icap]=pedCan->getValue(detid,icap);
0207       m_AllPedVals[detid]=PedVals;
0208     }
0209   }
0210 }
0211 */
0212 //-----------------------------------------------------------------------------
0213 void HcalLedAnalysis::GetLedConst(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0214   double time2 = 0;
0215   double time1 = 0;
0216   double time3 = 0;
0217   double time4 = 0;
0218   double dtime2 = 0;
0219   double dtime1 = 0;
0220   double dtime3 = 0;
0221   double dtime4 = 0;
0222   char output[256]{0};
0223 
0224   if (!m_outputFileText.empty()) {
0225     if (m_fitflag == 0 || m_fitflag == 2)
0226       m_outFile << "Det Eta,Phi,D   Mean    Error" << std::endl;
0227     else if (m_fitflag == 1 || m_fitflag == 3)
0228       m_outFile << "Det Eta,Phi,D   Peak    Error" << std::endl;
0229     else if (m_fitflag == 4)
0230       m_outFile << "Det Eta,Phi,D   Mean    Error      Peak    Error       MeanEv  Error       PeakEv  Error"
0231                 << std::endl;
0232   }
0233   for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
0234     // scale the LED pulse to 1 event
0235     _meol->second[10].first->Scale(1. / evt_curr);
0236     if (m_fitflag == 0 || m_fitflag == 4) {
0237       time1 = _meol->second[10].first->GetMean();
0238       dtime1 = _meol->second[10].first->GetRMS() / sqrt((float)evt_curr * (m_endTS - m_startTS + 1));
0239     }
0240     if (m_fitflag == 1 || m_fitflag == 4) {
0241       // put proper errors
0242       for (int j = 0; j < 10; j++)
0243         _meol->second[10].first->SetBinError(j + 1, _meol->second[j].first->GetRMS() / sqrt((float)evt_curr));
0244     }
0245     if (m_fitflag == 1 || m_fitflag == 3 || m_fitflag == 4) {
0246       _meol->second[10].first->Fit("landau", "Q");
0247       //      _meol->second[10].first->Fit("gaus","Q");
0248       TF1* fit = _meol->second[10].first->GetFunction("landau");
0249       //      TF1 *fit = _meol->second[10].first->GetFunction("gaus");
0250       time2 = fit->GetParameter(1);
0251       dtime2 = fit->GetParError(1);
0252     }
0253     if (m_fitflag == 2 || m_fitflag == 4) {
0254       time3 = _meol->second[12].first->GetMean();
0255       dtime3 = _meol->second[12].first->GetRMS() / sqrt((float)_meol->second[12].first->GetEntries());
0256     }
0257     if (m_fitflag == 3 || m_fitflag == 4) {
0258       time4 = _meol->second[13].first->GetMean();
0259       dtime4 = _meol->second[13].first->GetRMS() / sqrt((float)_meol->second[13].first->GetEntries());
0260     }
0261     for (int i = 0; i < 10; i++) {
0262       _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
0263       _meol->second[i].first->GetYaxis()->SetTitle("Counts");
0264       //      if(m_hiSaveflag>0)_meol->second[i].first->Write();
0265     }
0266     _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
0267     _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
0268     if (m_hiSaveflag > 0)
0269       _meol->second[10].first->Write();
0270     _meol->second[10].second.first[0].push_back(time1);
0271     _meol->second[10].second.first[1].push_back(dtime1);
0272     _meol->second[11].second.first[0].push_back(time2);
0273     _meol->second[11].second.first[1].push_back(dtime2);
0274     _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
0275     _meol->second[12].first->GetYaxis()->SetTitle("Counts");
0276     if (m_fitflag == 2 && m_hiSaveflag > 0)
0277       _meol->second[12].first->Write();
0278     _meol->second[12].second.first[0].push_back(time3);
0279     _meol->second[12].second.first[1].push_back(dtime3);
0280     _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
0281     _meol->second[13].first->GetYaxis()->SetTitle("Counts");
0282     if (m_fitflag > 2 && m_hiSaveflag > 0)
0283       _meol->second[13].first->Write();
0284     _meol->second[13].second.first[0].push_back(time4);
0285     _meol->second[13].second.first[1].push_back(dtime4);
0286     _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
0287     _meol->second[14].first->GetYaxis()->SetTitle("Counts");
0288     if (m_fitflag > 2 && m_hiSaveflag > 0)
0289       _meol->second[14].first->Write();
0290     _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
0291     _meol->second[15].first->GetYaxis()->SetTitle("Counts");
0292     if (m_fitflag > 2 && m_hiSaveflag > 0)
0293       _meol->second[15].first->Write();
0294     _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
0295     _meol->second[16].first->Write();
0296 
0297     // Ascii printout (need to modify to include new info)
0298     HcalDetId detid = _meol->first;
0299 
0300     if (!m_outputFileText.empty()) {
0301       if (m_fitflag == 0) {
0302         m_outFile << detid << "   " << time1 << " " << dtime1 << std::endl;
0303         m_outputFileXML << "  <DATA_SET>" << endl;
0304         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
0305         m_outputFileXML << "    <CHANNEL>" << endl;
0306         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0307         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0308         m_outputFileXML << output << endl;
0309         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0310         m_outputFileXML << output << endl;
0311         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0312         m_outputFileXML << output << endl;
0313         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0314         m_outputFileXML << output << endl;
0315 
0316         if (detid.subdet() == 1)
0317           m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0318         if (detid.subdet() == 2)
0319           m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0320         if (detid.subdet() == 3)
0321           m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0322         if (detid.subdet() == 4)
0323           m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0324         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0325         m_outputFileXML << output << endl;
0326         m_outputFileXML << "    </CHANNEL>" << endl;
0327         m_outputFileXML << "    <DATA>" << endl;
0328         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time1);
0329         m_outputFileXML << output << endl;
0330         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0331         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
0332         m_outputFileXML << output << endl;
0333         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0334         m_outputFileXML << output << endl;
0335         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
0336         m_outputFileXML << "    </DATA>" << endl;
0337         m_outputFileXML << "  </DATA_SET>" << endl;
0338 
0339       } else if (m_fitflag == 1) {
0340         m_outFile << detid << "   " << time2 << " " << dtime2 << std::endl;
0341         m_outputFileXML << "  <DATA_SET>" << endl;
0342         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
0343         m_outputFileXML << "    <CHANNEL>" << endl;
0344         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0345         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0346         m_outputFileXML << output << endl;
0347         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0348         m_outputFileXML << output << endl;
0349         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0350         m_outputFileXML << output << endl;
0351         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0352         m_outputFileXML << output << endl;
0353         if (detid.subdet() == 1)
0354           m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0355         if (detid.subdet() == 2)
0356           m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0357         if (detid.subdet() == 3)
0358           m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0359         if (detid.subdet() == 4)
0360           m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0361         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0362         m_outputFileXML << output << endl;
0363         m_outputFileXML << "    </CHANNEL>" << endl;
0364         m_outputFileXML << "    <DATA>" << endl;
0365         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time2);
0366         m_outputFileXML << output << endl;
0367         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0368         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime2);
0369         m_outputFileXML << output << endl;
0370         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0371         m_outputFileXML << output << endl;
0372         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
0373         m_outputFileXML << "    </DATA>" << endl;
0374         m_outputFileXML << "  </DATA_SET>" << endl;
0375       }
0376 
0377       else if (m_fitflag == 2) {
0378         m_outFile << detid << "   " << time3 << " " << dtime3 << std::endl;
0379         m_outputFileXML << "  <DATA_SET>" << endl;
0380         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
0381         m_outputFileXML << "    <CHANNEL>" << endl;
0382         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0383         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0384         m_outputFileXML << output << endl;
0385         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0386         m_outputFileXML << output << endl;
0387         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0388         m_outputFileXML << output << endl;
0389         snprintf(output, sizeof output, "      <Z>%2i</Z>", detid.zside());
0390         m_outputFileXML << output << endl;
0391         if (detid.subdet() == 1)
0392           m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0393         if (detid.subdet() == 2)
0394           m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0395         if (detid.subdet() == 3)
0396           m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0397         if (detid.subdet() == 4)
0398           m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0399         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0400         m_outputFileXML << output << endl;
0401         m_outputFileXML << "    </CHANNEL>" << endl;
0402         m_outputFileXML << "    <DATA>" << endl;
0403         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time3);
0404         m_outputFileXML << output << endl;
0405         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0406         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime3);
0407         m_outputFileXML << output << endl;
0408         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0409         m_outputFileXML << output << endl;
0410         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
0411         m_outputFileXML << "    </DATA>" << endl;
0412         m_outputFileXML << "  </DATA_SET>" << endl;
0413       } else if (m_fitflag == 3) {
0414         m_outFile << detid << "   " << time4 << " " << dtime4 << std::endl;
0415         m_outputFileXML << "  <DATA_SET>" << endl;
0416         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
0417         m_outputFileXML << "    <CHANNEL>" << endl;
0418         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0419         snprintf(output, sizeof output, "      <ETA>%2i</ETA>", detid.ietaAbs());
0420         m_outputFileXML << output << endl;
0421         snprintf(output, sizeof output, "      <PHI>%2i</PHI>", detid.iphi());
0422         m_outputFileXML << output << endl;
0423         snprintf(output, sizeof output, "      <DEPTH>%2i</DEPTH>", detid.depth());
0424         m_outputFileXML << output << endl;
0425         sprintf(output, "      <Z>%2i</Z>", detid.zside());
0426         m_outputFileXML << output << endl;
0427         if (detid.subdet() == 1)
0428           m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0429         if (detid.subdet() == 2)
0430           m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0431         if (detid.subdet() == 3)
0432           m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0433         if (detid.subdet() == 4)
0434           m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0435         snprintf(output, sizeof output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0436         m_outputFileXML << output << endl;
0437         m_outputFileXML << "    </CHANNEL>" << endl;
0438         m_outputFileXML << "    <DATA>" << endl;
0439         snprintf(output, sizeof output, "      <MEAN_TIME>%7f</MEAN_TIME>", time4);
0440         m_outputFileXML << output << endl;
0441         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0442         snprintf(output, sizeof output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
0443         m_outputFileXML << output << endl;
0444         snprintf(output, sizeof output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0445         m_outputFileXML << output << endl;
0446         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
0447         m_outputFileXML << "    </DATA>" << endl;
0448         m_outputFileXML << "  </DATA_SET>" << endl;
0449       }
0450 
0451       else if (m_fitflag == 4) {
0452         m_outFile << detid << "   " << time1 << " " << dtime1 << "   " << time2 << " " << dtime2 << "   " << time3
0453                   << " " << dtime3 << "   " << time4 << " " << dtime4 << std::endl;
0454       }
0455     }
0456   }
0457 }
0458 
0459 //-----------------------------------------------------------------------------
0460 void HcalLedAnalysis::LedSampleAnalysis() {
0461   // it is called every m_nevtsample events (a sample) and the end of run
0462   char LedSampleNum[20];
0463 
0464   snprintf(LedSampleNum, sizeof LedSampleNum, "LedSample_%d", sample);
0465   m_file->cd();
0466   m_file->mkdir(LedSampleNum);
0467   m_file->cd(LedSampleNum);
0468 
0469   // Compute LED constants for each HB/HE, HO, HF
0470   GetLedConst(hbHists.LEDTRENDS);
0471   GetLedConst(hoHists.LEDTRENDS);
0472   GetLedConst(hfHists.LEDTRENDS);
0473 }
0474 
0475 //-----------------------------------------------------------------------------
0476 void HcalLedAnalysis::LedTrendings(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0477   for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
0478     char name[1024];
0479     HcalDetId detid = _meol->first;
0480     snprintf(name, sizeof name, "LED timing trend, eta=%d phi=%d depth=%d", detid.ieta(), detid.iphi(), detid.depth());
0481     int bins = _meol->second[10 + m_fitflag].second.first[0].size();
0482     float lo = 0.5;
0483     float hi = (float)bins + 0.5;
0484     _meol->second[10 + m_fitflag].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0485 
0486     std::vector<double>::iterator sample_it;
0487     // LED timing - put content and errors
0488     int j = 0;
0489     for (sample_it = _meol->second[10 + m_fitflag].second.first[0].begin();
0490          sample_it != _meol->second[10 + m_fitflag].second.first[0].end();
0491          ++sample_it) {
0492       _meol->second[10 + m_fitflag].second.second[0]->SetBinContent(++j, *sample_it);
0493     }
0494     j = 0;
0495     for (sample_it = _meol->second[10 + m_fitflag].second.first[1].begin();
0496          sample_it != _meol->second[10 + m_fitflag].second.first[1].end();
0497          ++sample_it) {
0498       _meol->second[10 + m_fitflag].second.second[0]->SetBinError(++j, *sample_it);
0499     }
0500     snprintf(name, sizeof name, "Sample (%d events)", m_nevtsample);
0501     _meol->second[10 + m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
0502     _meol->second[10 + m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
0503     _meol->second[10 + m_fitflag].second.second[0]->Write();
0504   }
0505 }
0506 
0507 //-----------------------------------------------------------------------------
0508 void HcalLedAnalysis::LedDone() {
0509   // First process the last sample (remaining events).
0510   if (evt % m_nevtsample != 0)
0511     LedSampleAnalysis();
0512 
0513   // Now do the end of run analysis: trending histos
0514   if (sample > 1 && m_fitflag != 4) {
0515     m_file->cd();
0516     m_file->cd("HBHE");
0517     LedTrendings(hbHists.LEDTRENDS);
0518     m_file->cd();
0519     m_file->cd("HO");
0520     LedTrendings(hoHists.LEDTRENDS);
0521     m_file->cd();
0522     m_file->cd("HF");
0523     LedTrendings(hfHists.LEDTRENDS);
0524   }
0525 
0526   // Write other histograms.
0527   // HB
0528   m_file->cd();
0529   m_file->cd("HBHE");
0530   hbHists.ALLLEDS->Write();
0531   hbHists.LEDRMS->Write();
0532   hbHists.LEDMEAN->Write();
0533   // HO
0534   m_file->cd();
0535   m_file->cd("HO");
0536   hoHists.ALLLEDS->Write();
0537   hoHists.LEDRMS->Write();
0538   hoHists.LEDMEAN->Write();
0539   // HF
0540   m_file->cd();
0541   m_file->cd("HF");
0542   hfHists.ALLLEDS->Write();
0543   hfHists.LEDRMS->Write();
0544   hfHists.LEDMEAN->Write();
0545   // Calib
0546   m_file->cd();
0547   m_file->cd("Calib");
0548   for (_meca = calibHists.begin(); _meca != calibHists.end(); _meca++) {
0549     _meca->second.avePulse->Write();
0550     _meca->second.integPulse->Write();
0551   }
0552 
0553   // Write the histo file and close it
0554   //  m_file->Write();
0555   m_file->Close();
0556   edm::LogInfo("HcalLedAnalysis") << "Hcal histograms written to " << m_outputFileROOT.c_str() << endl;
0557 }
0558 
0559 //-----------------------------------------------------------------------------
0560 void HcalLedAnalysis::processLedEvent(const HBHEDigiCollection& hbhe,
0561                                       const HODigiCollection& ho,
0562                                       const HFDigiCollection& hf,
0563                                       const HcalCalibDigiCollection& calib,
0564                                       const HcalDbService& cond) {
0565   evt++;
0566   sample = (evt - 1) / m_nevtsample + 1;
0567   evt_curr = evt % m_nevtsample;
0568   if (evt_curr == 0)
0569     evt_curr = m_nevtsample;
0570 
0571   // Calib
0572 
0573   if (m_usecalib) {
0574     if (calib.empty()) {
0575       edm::LogError("HcalLedAnalysis") << "Event with " << (int)calib.size() << " Calib Digis passed.";
0576       return;
0577     }
0578     // this is effectively a loop over electronic channels
0579     for (HcalCalibDigiCollection::const_iterator j = calib.begin(); j != calib.end(); ++j) {
0580       const HcalCalibDataFrame digi = (const HcalCalibDataFrame)(*j);
0581       HcalElectronicsId elecId = digi.elecId();
0582       HcalCalibDetId calibId = digi.id();
0583       ProcessCalibEvent(elecId.fiberChanId(),
0584                         calibId,
0585                         digi);  //Shouldn't depend on anything in elecId but not sure how else to do it
0586     }
0587   }
0588 
0589   // HB + HE
0590   if (hbhe.empty()) {
0591     edm::LogError("HcalLedAnalysis") << "Event with " << (int)hbhe.size() << " HBHE Digis passed.";
0592     return;
0593   }
0594   // this is effectively a loop over electronic channels
0595   for (HBHEDigiCollection::const_iterator j = hbhe.begin(); j != hbhe.end(); ++j) {
0596     const HBHEDataFrame digi = (const HBHEDataFrame)(*j);
0597     for (int k = 0; k < (int)state.size(); k++)
0598       state[k] = true;
0599     // See if histos exist for this channel, and if not, create them
0600     _meol = hbHists.LEDTRENDS.find(digi.id());
0601     if (_meol == hbHists.LEDTRENDS.end()) {
0602       SetupLEDHists(0, digi.id(), hbHists.LEDTRENDS);
0603     }
0604     LedHBHEHists(digi.id(), digi, hbHists.LEDTRENDS, cond);
0605   }
0606 
0607   // HO
0608   if (ho.empty()) {
0609     edm::LogError("HcalLedAnalysis") << "Event with " << (int)ho.size() << " HO Digis passed.";
0610     return;
0611   }
0612   for (HODigiCollection::const_iterator j = ho.begin(); j != ho.end(); ++j) {
0613     const HODataFrame digi = (const HODataFrame)(*j);
0614     _meol = hoHists.LEDTRENDS.find(digi.id());
0615     if (_meol == hoHists.LEDTRENDS.end()) {
0616       SetupLEDHists(1, digi.id(), hoHists.LEDTRENDS);
0617     }
0618     LedHOHists(digi.id(), digi, hoHists.LEDTRENDS, cond);
0619   }
0620 
0621   // HF
0622   if (hf.empty()) {
0623     edm::LogError("HcalLedAnalysis") << "Event with " << (int)hf.size() << " HF Digis passed.";
0624     return;
0625   }
0626   for (HFDigiCollection::const_iterator j = hf.begin(); j != hf.end(); ++j) {
0627     const HFDataFrame digi = (const HFDataFrame)(*j);
0628     _meol = hfHists.LEDTRENDS.find(digi.id());
0629     if (_meol == hfHists.LEDTRENDS.end()) {
0630       SetupLEDHists(2, digi.id(), hfHists.LEDTRENDS);
0631     }
0632     LedHFHists(digi.id(), digi, hfHists.LEDTRENDS, cond);
0633   }
0634 
0635   // Call the function every m_nevtsample events
0636   if (evt % m_nevtsample == 0)
0637     LedSampleAnalysis();
0638 }
0639 //----------------------------------------------------------------------------
0640 void HcalLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0641   string type = "HBHE";
0642   if (id == 1)
0643     type = "HO";
0644   if (id == 2)
0645     type = "HF";
0646 
0647   _meol = toolT.find(detid);
0648   if (_meol == toolT.end()) {
0649     // if histos for this channel do not exist, create them
0650     map<int, LEDBUNCH> insert;
0651     char name[1024];
0652     for (int i = 0; i < 10; i++) {
0653       snprintf(name,
0654                sizeof name,
0655                "%s Pulse height, eta=%d phi=%d depth=%d TS=%d",
0656                type.c_str(),
0657                detid.ieta(),
0658                detid.iphi(),
0659                detid.depth(),
0660                i);
0661       insert[i].first = new TH1F(name, name, 200, 0., 2000.);
0662     }
0663     snprintf(name,
0664              sizeof name,
0665              "%s LED Mean pulse, eta=%d phi=%d depth=%d",
0666              type.c_str(),
0667              detid.ieta(),
0668              detid.iphi(),
0669              detid.depth());
0670     insert[10].first = new TH1F(name, name, 10, -0.5, 9.5);
0671     snprintf(name,
0672              sizeof name,
0673              "%s LED Pulse, eta=%d phi=%d depth=%d",
0674              type.c_str(),
0675              detid.ieta(),
0676              detid.iphi(),
0677              detid.depth());
0678     insert[11].first = new TH1F(name, name, 10, -0.5, 9.5);
0679     snprintf(name,
0680              sizeof name,
0681              "%s Mean TS, eta=%d phi=%d depth=%d",
0682              type.c_str(),
0683              detid.ieta(),
0684              detid.iphi(),
0685              detid.depth());
0686     insert[12].first = new TH1F(name, name, 200, 0., 10.);
0687     snprintf(name,
0688              sizeof name,
0689              "%s Peak TS, eta=%d phi=%d depth=%d",
0690              type.c_str(),
0691              detid.ieta(),
0692              detid.iphi(),
0693              detid.depth());
0694     insert[13].first = new TH1F(name, name, 200, 0., 10.);
0695     snprintf(name,
0696              sizeof name,
0697              "%s Peak TS error, eta=%d phi=%d depth=%d",
0698              type.c_str(),
0699              detid.ieta(),
0700              detid.iphi(),
0701              detid.depth());
0702     insert[14].first = new TH1F(name, name, 200, 0., 0.05);
0703     snprintf(name,
0704              sizeof name,
0705              "%s Fit chi2, eta=%d phi=%d depth=%d",
0706              type.c_str(),
0707              detid.ieta(),
0708              detid.iphi(),
0709              detid.depth());
0710     insert[15].first = new TH1F(name, name, 100, 0., 50.);
0711     snprintf(name,
0712              sizeof name,
0713              "%s Integrated Signal, eta=%d phi=%d depth=%d",
0714              type.c_str(),
0715              detid.ieta(),
0716              detid.iphi(),
0717              detid.depth());
0718     insert[16].first = new TH1F(name, name, 500, 0., 5000.);
0719 
0720     toolT[detid] = insert;
0721   }
0722 }
0723 //-----------------------------------------------------------------------------
0724 void HcalLedAnalysis::LedHBHEHists(const HcalDetId& detid,
0725                                    const HBHEDataFrame& ledDigi,
0726                                    map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0727                                    const HcalDbService& cond) {
0728   map<int, LEDBUNCH> _mei;
0729   _meol = toolT.find(detid);
0730   _mei = _meol->second;
0731 
0732   // Reset the histos if we're at the end of a 'bunch'
0733   if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0734     for (int k = 0; k < (int)state.size(); k++)
0735       state[k] = false;
0736     for (int i = 0; i < 16; i++)
0737       _mei[i].first->Reset();
0738   }
0739 
0740   // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
0741 
0742   //  int maxTS = -1;
0743   float max_fC = 0;
0744   float ta = 0;
0745   m_coder = cond.getHcalCoder(detid);
0746   m_ped = cond.getPedestal(detid);
0747   m_shape = cond.getHcalShape(m_coder);
0748   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0749     int capid = ledDigi[TS].capid();
0750     int adc = ledDigi[TS].adc();
0751     double fC = m_coder->charge(*m_shape, adc, capid);
0752     ta = (fC - m_ped->getValue(capid));
0753     //cout << "DetID: " << detid << "  CapID: " << capid << "  ADC: " << adc << "  fC: " << fC << endl;
0754     _mei[TS].first->Fill(ta);
0755     _mei[10].first->AddBinContent(TS + 1, ta);  // This is average pulse, could probably do better (Profile?)
0756     if (m_fitflag > 1) {
0757       if (TS == m_startTS)
0758         _mei[11].first->Reset();
0759       _mei[11].first->SetBinContent(TS + 1, ta);
0760     }
0761     // keep track of max TS and max amplitude (in fC)
0762     if (ta > max_fC) {
0763       max_fC = ta;
0764       //      maxTS = TS;
0765     }
0766   }
0767 
0768   // Now we have a sample with pedestals subtracted and in units of fC
0769   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
0770   // we now want to use Phil's timing correction.  This is not necessary
0771   // if we are performing a Landau fit (m_fitflag = 3)
0772 
0773   float sum = 0.;
0774   for (int i = 0; i < 10; i++)
0775     sum = sum + _mei[11].first->GetBinContent(i + 1);
0776   if (sum > 100) {
0777     if (m_fitflag == 2 || m_fitflag == 4) {
0778       float timmean = _mei[11].first->GetMean();  // let's use Phil's way instead
0779       float timmeancorr = BinsizeCorr(timmean);
0780       _mei[12].first->Fill(timmeancorr);
0781     }
0782     _mei[16].first->Fill(
0783         _mei[11].first->Integral());  // Integrated charge (may be more usfull to convert to Energy first?)
0784     if (m_fitflag == 3 || m_fitflag == 4) {
0785       _mei[11].first->Fit("landau", "Q");
0786       TF1* fit = _mei[11].first->GetFunction("landau");
0787       _mei[13].first->Fill(fit->GetParameter(1));
0788       _mei[14].first->Fill(fit->GetParError(1));
0789       _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0790     }
0791   }
0792 }
0793 
0794 //-----------------------------------------------------------------------------
0795 void HcalLedAnalysis::LedHOHists(const HcalDetId& detid,
0796                                  const HODataFrame& ledDigi,
0797                                  map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0798                                  const HcalDbService& cond) {
0799   map<int, LEDBUNCH> _mei;
0800   _meol = toolT.find(detid);
0801   _mei = _meol->second;
0802   // Rest the histos if we're at the end of a 'bunch'
0803   if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0804     for (int k = 0; k < (int)state.size(); k++)
0805       state[k] = false;
0806     for (int i = 0; i < 16; i++)
0807       _mei[i].first->Reset();
0808   }
0809 
0810   // now we have the signal in fC, let's get rid of that darn pedestal
0811   // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
0812 
0813   //  int maxTS = -1;
0814   float max_fC = 0;
0815   float ta = 0;
0816   m_coder = cond.getHcalCoder(detid);
0817   m_ped = cond.getPedestal(detid);
0818   m_shape = cond.getHcalShape(m_coder);
0819   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0820     int capid = ledDigi[TS].capid();
0821     int adc = ledDigi[TS].adc();
0822     double fC = m_coder->charge(*m_shape, adc, capid);
0823     ta = (fC - m_ped->getValue(capid));
0824     _mei[TS].first->Fill(ta);
0825     _mei[10].first->AddBinContent(TS + 1, ta);  // This is average pulse, could probably do better (Profile?)
0826     if (m_fitflag > 1) {
0827       if (TS == m_startTS)
0828         _mei[11].first->Reset();
0829       _mei[11].first->SetBinContent(TS + 1, ta);
0830     }
0831     // keep track of max TS and max amplitude (in fC)
0832     if (ta > max_fC) {
0833       max_fC = ta;
0834       //      maxTS = TS;
0835     }
0836   }
0837 
0838   // Now we have a sample with pedestals subtracted and in units of fC
0839   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
0840   // we now want to use Phil's timing correction.  This is not necessary
0841   // if we are performing a Landau fit (m_fitflag = 3)
0842 
0843   float sum = 0.;
0844   for (int i = 0; i < 10; i++)
0845     sum = sum + _mei[11].first->GetBinContent(i + 1);
0846   if (sum > 100) {
0847     if (m_fitflag == 2 || m_fitflag == 4) {
0848       float timmean = _mei[11].first->GetMean();  // let's use Phil's way instead
0849       float timmeancorr = BinsizeCorr(timmean);
0850       _mei[12].first->Fill(timmeancorr);
0851     }
0852     _mei[16].first->Fill(
0853         _mei[11].first->Integral());  // Integrated charge (may be more usfull to convert to Energy first?)
0854     if (m_fitflag == 3 || m_fitflag == 4) {
0855       _mei[11].first->Fit("landau", "Q");
0856       TF1* fit = _mei[11].first->GetFunction("landau");
0857       _mei[13].first->Fill(fit->GetParameter(1));
0858       _mei[14].first->Fill(fit->GetParError(1));
0859       _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0860     }
0861   }
0862 }
0863 
0864 //-----------------------------------------------------------------------------
0865 void HcalLedAnalysis::LedHFHists(const HcalDetId& detid,
0866                                  const HFDataFrame& ledDigi,
0867                                  map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0868                                  const HcalDbService& cond) {
0869   map<int, LEDBUNCH> _mei;
0870   _meol = toolT.find(detid);
0871   _mei = _meol->second;
0872   // Rest the histos if we're at the end of a 'bunch'
0873   if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0874     for (int k = 0; k < (int)state.size(); k++)
0875       state[k] = false;
0876     for (int i = 0; i < 16; i++)
0877       _mei[i].first->Reset();
0878   }
0879 
0880   // now we have the signal in fC, let's get rid of that darn pedestal
0881   // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
0882 
0883   //  int maxTS = -1;
0884   float max_fC = 0;
0885   float ta = 0;
0886   m_coder = cond.getHcalCoder(detid);
0887   m_ped = cond.getPedestal(detid);
0888   m_shape = cond.getHcalShape(m_coder);
0889   //cout << "New Digi!!!!!!!!!!!!!!!!!!!!!!" << endl;
0890   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0891     int capid = ledDigi[TS].capid();
0892     // BE CAREFUL: this is assuming peds are stored in ADCs
0893     int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
0894     if (adc < 0) {
0895       adc = 0;
0896     }  // to prevent negative adcs after ped subtraction, which should really only happen
0897     // if you're using the wrong peds.
0898     double fC = m_coder->charge(*m_shape, adc, capid);
0899     //ta = (fC - m_ped->getValue(capid));
0900     ta = fC;
0901     //cout << "DetID: " << detid << "  CapID: " << capid << "  ADC: " << adc << "  Ped: " << m_ped->getValue(capid) << "  fC: " << fC << endl;
0902     _mei[TS].first->Fill(ta);
0903     _mei[10].first->AddBinContent(TS + 1, ta);  // This is average pulse, could probably do better (Profile?)
0904     if (m_fitflag > 1) {
0905       if (TS == m_startTS)
0906         _mei[11].first->Reset();
0907       _mei[11].first->SetBinContent(TS + 1, ta);
0908     }
0909 
0910     // keep track of max TS and max amplitude (in fC)
0911     if (ta > max_fC) {
0912       max_fC = ta;
0913       //      maxTS = TS;
0914     }
0915   }
0916 
0917   // Now we have a sample with pedestals subtracted and in units of fC
0918   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
0919   // we now want to use Phil's timing correction.  This is not necessary
0920   // if we are performing a Landau fit (m_fitflag = 3)
0921 
0922   float sum = 0.;
0923   for (int i = 0; i < 10; i++)
0924     sum = sum + _mei[11].first->GetBinContent(i + 1);
0925   if (sum > 100) {
0926     if (m_fitflag == 2 || m_fitflag == 4) {
0927       float timmean = _mei[11].first->GetMean();  // let's use Phil's way instead
0928       float timmeancorr = BinsizeCorr(timmean);
0929       _mei[12].first->Fill(timmeancorr);
0930     }
0931     _mei[16].first->Fill(
0932         _mei[11].first->Integral());  // Integrated charge (may be more usfull to convert to Energy first?)
0933     if (m_fitflag == 3 || m_fitflag == 4) {
0934       _mei[11].first->Fit("landau", "Q");
0935       TF1* fit = _mei[11].first->GetFunction("landau");
0936       _mei[13].first->Fill(fit->GetParameter(1));
0937       _mei[14].first->Fill(fit->GetParError(1));
0938       _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0939     }
0940   }
0941 }
0942 
0943 //-----------------------------------------------------------------------------
0944 float HcalLedAnalysis::BinsizeCorr(float time) {
0945   // this is the bin size correction to be applied for laser data (from Andy),
0946   // it comes from a pulse shape measured from TB04 data (from Jordan)
0947   // This should eventually be replaced with the more thorough treatment from Phil
0948 
0949   float corrtime = 0.;
0950   static const float tstrue[32] = {0.003, 0.03425, 0.06548, 0.09675, 0.128, 0.15925, 0.1905, 0.22175,
0951                                    0.253, 0.28425, 0.3155,  0.34675, 0.378, 0.40925, 0.4405, 0.47175,
0952                                    0.503, 0.53425, 0.5655,  0.59675, 0.628, 0.65925, 0.6905, 0.72175,
0953                                    0.753, 0.78425, 0.8155,  0.84675, 0.878, 0.90925, 0.9405, 0.97175};
0954   static const float tsreco[32] = {-0.00422, 0.01815, 0.04409, 0.07346, 0.09799, 0.12192, 0.15072, 0.18158,
0955                                    0.21397,  0.24865, 0.28448, 0.31973, 0.35449, 0.39208, 0.43282, 0.47244,
0956                                    0.5105,   0.55008, 0.58827, 0.62828, 0.6717,  0.70966, 0.74086, 0.77496,
0957                                    0.80843,  0.83472, 0.86044, 0.8843,  0.90674, 0.92982, 0.95072, 0.9726};
0958 
0959   int inttime = (int)time;
0960   float restime = time - inttime;
0961   for (int i = 0; i <= 32; i++) {
0962     float lolim = 0.;
0963     float uplim = 1.;
0964     float tsdown;
0965     float tsup;
0966     if (i > 0) {
0967       lolim = tsreco[i - 1];
0968       tsdown = tstrue[i - 1];
0969     } else
0970       tsdown = tstrue[31] - 1.;
0971     if (i < 32) {
0972       uplim = tsreco[i];
0973       tsup = tstrue[i];
0974     } else
0975       tsup = tstrue[0] + 1.;
0976     if (restime >= lolim && restime < uplim) {
0977       corrtime = (tsdown * (uplim - restime) + tsup * (restime - lolim)) / (uplim - lolim);
0978     }
0979   }
0980   corrtime += inttime;
0981 
0982   return corrtime;
0983 }
0984 //-----------------------------------------------------------------------------
0985 
0986 // Will try to implement Phil's time slew correction here at some point
0987 
0988 //-----------------------------------------------------------------------------
0989 void HcalLedAnalysis::ProcessCalibEvent(int fiberChan, HcalCalibDetId calibId, const HcalCalibDataFrame& digi) {
0990   _meca = calibHists.find(calibId);
0991   if (_meca == calibHists.end()) {
0992     // if histos for this channel do not exist, first create them
0993     char name[1024];
0994     std::string prefix;
0995     if (calibId.calibFlavor() == HcalCalibDetId::CalibrationBox) {
0996       std::string sector = (calibId.hcalSubdet() == HcalBarrel)    ? ("HB")
0997                            : (calibId.hcalSubdet() == HcalEndcap)  ? ("HE")
0998                            : (calibId.hcalSubdet() == HcalOuter)   ? ("HO")
0999                            : (calibId.hcalSubdet() == HcalForward) ? ("HF")
1000                                                                    : "";
1001       snprintf(name,
1002                sizeof name,
1003                "%s %+d iphi=%d %s",
1004                sector.c_str(),
1005                calibId.ieta(),
1006                calibId.iphi(),
1007                calibId.cboxChannelString().c_str());
1008       prefix = name;
1009     }
1010 
1011     snprintf(name, sizeof name, "%s Pin Diode Mean", prefix.c_str());
1012     calibHists[calibId].avePulse = new TProfile(name, name, 10, -0.5, 9.5, 0, 1000);
1013     snprintf(name, sizeof name, "%s Pin Diode Current Pulse", prefix.c_str());
1014     calibHists[calibId].thisPulse = new TH1F(name, name, 10, -0.5, 9.5);
1015     snprintf(name, sizeof name, "%s Pin Diode Integrated Pulse", prefix.c_str());
1016     calibHists[calibId].integPulse = new TH1F(name, name, 200, 0, 500);
1017   } else {
1018     for (int i = m_startTS; i < digi.size() && i <= m_endTS; i++) {
1019       calibHists[calibId].avePulse->Fill(i, digi.sample(i).adc());
1020       calibHists[calibId].thisPulse->SetBinContent(i + 1, digi.sample(i).adc());
1021     }
1022     calibHists[calibId].integPulse->Fill(calibHists[calibId].thisPulse->Integral());
1023   }
1024 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!