Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-06-11 04:36:46 

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