Reference Manual

Go to the documentation of this file.
00001 
00002 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
00003 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
00004 #include "CondFormats/HcalObjects/interface/HcalPedestals.h"
00005 
00006 #include "CalibCalorimetry/HcalAlgos/interface/HcalLedAnalysis.h"
00007 #include "TFile.h"
00008 #include <math.h>
00009 using namespace std;
00010 
00011 
00012 HcalLedAnalysis::HcalLedAnalysis(const edm::ParameterSet& ps)
00013 {
00014   // init
00015 
00016   m_coder = 0;
00017   m_ped   = 0;
00018   m_shape = 0;
00019   evt=0;
00020   sample=0;
00021   m_file=0;
00022   // output files
00023   for(int k=0;k<4;k++) state.push_back(true); // 4 cap-ids (do we care?)
00024   m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
00025   m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML","");
00026   if ( m_outputFileText.size() != 0 ) {
00027     cout << "Hcal LED results will be saved to " << m_outputFileText.c_str() << endl;
00028     m_outFile.open(m_outputFileText.c_str());
00029   } 
00030   m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
00031   if ( m_outputFileROOT.size() != 0 ) {
00032     cout << "Hcal LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
00033   }
00034 
00035   m_nevtsample = ps.getUntrackedParameter<int>("nevtsample",9999999);
00036   if(m_nevtsample<1)m_nevtsample=9999999;
00037   m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag",0);
00038   if(m_hiSaveflag<0)m_hiSaveflag=0;
00039   if(m_hiSaveflag>0)m_hiSaveflag=1;
00040   m_fitflag = ps.getUntrackedParameter<int>("analysisflag",2);
00041   if(m_fitflag<0)m_fitflag=0;
00042   if(m_fitflag>4)m_fitflag=4;
00043   m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
00044   if(m_startTS<0) m_startTS=0;
00045   m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
00046   m_usecalib = ps.getUntrackedParameter<bool>("usecalib",false);
00047   m_logFile.open("HcalLedAnalysis.log");
00048 
00049   int runNum = ps.getUntrackedParameter<int>("runNumber",999999);
00050 
00051   // histogram booking
00052   hbHists.ALLLEDS = new TH1F("HBHE All LEDs","HB/HE All Leds",10,0,9);
00053   hbHists.LEDRMS= new TH1F("HBHE All LED RMS","HB/HE All LED RMS",100,0,3);
00054   hbHists.LEDMEAN= new TH1F("HBHE All LED Means","HB/HE All LED Means",100,0,9);
00055   hbHists.CHI2= new TH1F("HBHE Chi2 by ndf for Landau fit","HB/HE Chi2/ndf Landau",200,0.,50.);
00056 
00057   hoHists.ALLLEDS = new TH1F("HO All LEDs","HO All Leds",10,0,9);
00058   hoHists.LEDRMS= new TH1F("HO All LED RMS","HO All LED RMS",100,0,3);
00059   hoHists.LEDMEAN= new TH1F("HO All LED Means","HO All LED Means",100,0,9);
00060   hoHists.CHI2= new TH1F("HO Chi2 by ndf for Landau fit","HO Chi2/ndf Landau",200,0.,50.);
00061 
00062   hfHists.ALLLEDS = new TH1F("HF All LEDs","HF All Leds",10,0,9);
00063   hfHists.LEDRMS= new TH1F("HF All LED RMS","HF All LED RMS",100,0,3);
00064   hfHists.LEDMEAN= new TH1F("HF All LED Means","HF All LED Means",100,0,9);
00065   hfHists.CHI2= new TH1F("HF Chi2 by ndf for Landau fit","HF Chi2/ndf Landau",200,0.,50.);
00066 
00067 
00068   //XML file header
00069   m_outputFileXML.open(m_outputFileX.c_str());
00070 
00071   m_outputFileXML << "<?xml version='1.0' encoding='UTF-8'?>" << endl;
00072 
00073   m_outputFileXML << "<ROOT>" << endl;
00074 
00075   m_outputFileXML << "  <HEADER>" << endl;
00076 
00077   m_outputFileXML << "    <TYPE>" << endl;
00078 
00079   m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>" << endl;
00080 
00081   m_outputFileXML << "      <NAME>HCAL LED Timing</NAME>" << endl;
00082 
00083   m_outputFileXML << "    </TYPE>" << endl;
00084 
00085   m_outputFileXML << "    <RUN>" << endl;
00086 
00087   m_outputFileXML << "      <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>" << endl;
00088 
00089   sprintf(output, "      <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
00090   m_outputFileXML << output << endl;
00091 
00092   m_outputFileXML << "      <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>" << endl;
00093 
00094   m_outputFileXML << "      <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
00095 
00096   m_outputFileXML << "    </RUN>" << endl;
00097 
00098   m_outputFileXML << "  </HEADER>" << endl;
00099 
00100   m_outputFileXML << "<!-- Tags secton -->" << endl;
00101 
00102   m_outputFileXML <<  "  <ELEMENTS>" << endl;
00103 
00104   m_outputFileXML <<  "    <DATA_SET id='-1'/>" << endl;
00105 
00106   m_outputFileXML << "      <IOV id='1'>" << endl;
00107 
00108   m_outputFileXML << "        <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>" << endl;
00109 
00110   m_outputFileXML <<  "        <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>" << endl;
00111 
00112   m_outputFileXML << "      </IOV>" << endl;
00113 
00114   m_outputFileXML << "      <TAG id='2' mode='auto'>" << endl;
00115 
00116   sprintf(output, "        <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
00117   m_outputFileXML << output << endl;
00118 
00119   m_outputFileXML << "        <DETECTOR_NAME>HCAL</DETECTOR_NAME>" << endl;
00120 
00121   m_outputFileXML << "        <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
00122 
00123   m_outputFileXML << "      </TAG>" << endl;
00124 
00125   m_outputFileXML << "  </ELEMENTS>" << endl;
00126 
00127   m_outputFileXML << "  <MAPS>" << endl;
00128 
00129   m_outputFileXML << "      <TAG idref ='2'>" << endl;
00130 
00131   m_outputFileXML << "        <IOV idref='1'>" << endl;
00132 
00133   m_outputFileXML << "          <DATA_SET idref='-1' />" << endl;
00134 
00135   m_outputFileXML << "        </IOV>" << endl;
00136 
00137   m_outputFileXML << "      </TAG>" << endl;
00138 
00139   m_outputFileXML <<   "  </MAPS>" << endl;
00140 
00141 }
00142 
00143 //-----------------------------------------------------------------------------
00144 HcalLedAnalysis::~HcalLedAnalysis(){
00146   for(_meol=hbHists.LEDTRENDS.begin(); _meol!=hbHists.LEDTRENDS.end(); _meol++){
00147     for(int i=0; i<15; i++) _meol->second[i].first->Delete();
00148   }
00149   for(_meol=hoHists.LEDTRENDS.begin(); _meol!=hoHists.LEDTRENDS.end(); _meol++){
00150     for(int i=0; i<15; i++) _meol->second[i].first->Delete();
00151   }
00152   for(_meol=hfHists.LEDTRENDS.begin(); _meol!=hfHists.LEDTRENDS.end(); _meol++){
00153     for(int i=0; i<15; i++) _meol->second[i].first->Delete();
00154   }
00155   hbHists.ALLLEDS->Delete();
00156   hbHists.LEDRMS->Delete();
00157   hbHists.LEDMEAN->Delete();
00158   hbHists.CHI2->Delete();
00159 
00160   hoHists.ALLLEDS->Delete();
00161   hoHists.LEDRMS->Delete();
00162   hoHists.LEDMEAN->Delete();
00163   hoHists.CHI2->Delete();
00164 
00165   hfHists.ALLLEDS->Delete();
00166   hfHists.LEDRMS->Delete();
00167   hfHists.LEDMEAN->Delete();
00168   hfHists.CHI2->Delete();
00169 }
00170 
00171 //-----------------------------------------------------------------------------
00172 void HcalLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
00173   // open the histogram file, create directories within
00174   m_file=new TFile(m_outputFileROOT.c_str(),"RECREATE");
00175   m_file->mkdir("HBHE");
00176   m_file->cd();
00177   m_file->mkdir("HO");
00178   m_file->cd();
00179   m_file->mkdir("HF");
00180   m_file->cd();
00181   m_file->mkdir("Calib");
00182   m_file->cd();
00183 }
00184 
00185 //-----------------------------------------------------------------------------
00186 /*
00187 void HcalLedAnalysis::doPeds(const HcalPedestal* fInputPedestals){
00188 // put all pedestals in a map m_AllPedVals, to be used in processLedEvent -
00189 // sorry, this is the only way I was able to implement pedestal subtraction
00190 
00191 // DEPRECATED
00192 // This is no longer useful, better ways of doing it -A
00193   HcalDetId detid;
00194   map<int,float> PedVals;
00195   pedCan = fInputPedestals;
00196   if(pedCan){
00197     std::vector<DetId> Channs=pedCan->getAllChannels();
00198     for (int i=0; i<(int)Channs.size(); i++){
00199       detid=HcalDetId (Channs[i]);
00200       for (int icap=0; icap<4; icap++) PedVals[icap]=pedCan->getValue(detid,icap);
00201       m_AllPedVals[detid]=PedVals;
00202     }
00203   }
00204 }
00205 */
00206 //-----------------------------------------------------------------------------
00207 void HcalLedAnalysis::GetLedConst(map<HcalDetId, map<int,LEDBUNCH> > &toolT){
00208   double time2=0; double time1=0; double time3=0; double time4=0;
00209   double dtime2=0; double dtime1=0; double dtime3=0; double dtime4=0;
00210 
00211   if (m_outputFileText!=""){
00212     if(m_fitflag==0 || m_fitflag==2) m_outFile<<"Det Eta,Phi,D   Mean    Error"<<std::endl;
00213     else if(m_fitflag==1 || m_fitflag==3) m_outFile<<"Det Eta,Phi,D   Peak    Error"<<std::endl;
00214     else if(m_fitflag==4) m_outFile<<"Det Eta,Phi,D   Mean    Error      Peak    Error       MeanEv  Error       PeakEv  Error"<<std::endl;
00215   }
00216   for(_meol=toolT.begin(); _meol!=toolT.end(); _meol++){
00217 // scale the LED pulse to 1 event
00218     _meol->second[10].first->Scale(1./evt_curr);
00219     if(m_fitflag==0 || m_fitflag==4){
00220       time1 = _meol->second[10].first->GetMean();
00221       dtime1 = _meol->second[10].first->GetRMS()/sqrt((float)evt_curr*(m_endTS-m_startTS+1));
00222     }
00223     if(m_fitflag==1 || m_fitflag==4){
00224 // put proper errors
00225       for(int j=0; j<10; j++) _meol->second[10].first->SetBinError(j+1,_meol->second[j].first->GetRMS()/sqrt((float)evt_curr));
00226     }
00227     if(m_fitflag==1 || m_fitflag==3 || m_fitflag==4){
00228       _meol->second[10].first->Fit("landau","Q");
00229 //      _meol->second[10].first->Fit("gaus","Q");
00230       TF1 *fit = _meol->second[10].first->GetFunction("landau");
00231 //      TF1 *fit = _meol->second[10].first->GetFunction("gaus");
00232       time2=fit->GetParameter(1);
00233       dtime2=fit->GetParError(1);
00234     }
00235     if(m_fitflag==2 || m_fitflag==4){
00236       time3 = _meol->second[12].first->GetMean();
00237       dtime3 = _meol->second[12].first->GetRMS()/sqrt((float)_meol->second[12].first->GetEntries());
00238     }
00239     if(m_fitflag==3 || m_fitflag==4){
00240       time4 = _meol->second[13].first->GetMean();
00241       dtime4 = _meol->second[13].first->GetRMS()/sqrt((float)_meol->second[13].first->GetEntries());
00242     }
00243     for (int i=0; i<10; i++){
00244       _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
00245       _meol->second[i].first->GetYaxis()->SetTitle("Counts");
00246 //      if(m_hiSaveflag>0)_meol->second[i].first->Write();
00247     }
00248     _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
00249     _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
00250     if(m_hiSaveflag>0)_meol->second[10].first->Write();
00251     _meol->second[10].second.first[0].push_back(time1);
00252     _meol->second[10].second.first[1].push_back(dtime1);
00253     _meol->second[11].second.first[0].push_back(time2);
00254     _meol->second[11].second.first[1].push_back(dtime2);
00255     _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
00256     _meol->second[12].first->GetYaxis()->SetTitle("Counts");
00257     if(m_fitflag==2 && m_hiSaveflag>0)_meol->second[12].first->Write();
00258     _meol->second[12].second.first[0].push_back(time3);
00259     _meol->second[12].second.first[1].push_back(dtime3);
00260     _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
00261     _meol->second[13].first->GetYaxis()->SetTitle("Counts");
00262     if(m_fitflag>2 && m_hiSaveflag>0)_meol->second[13].first->Write();
00263     _meol->second[13].second.first[0].push_back(time4);
00264     _meol->second[13].second.first[1].push_back(dtime4);
00265     _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
00266     _meol->second[14].first->GetYaxis()->SetTitle("Counts");
00267     if(m_fitflag>2 && m_hiSaveflag>0)_meol->second[14].first->Write();
00268     _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
00269     _meol->second[15].first->GetYaxis()->SetTitle("Counts");
00270     if(m_fitflag>2 && m_hiSaveflag>0)_meol->second[15].first->Write();
00271     _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
00272     _meol->second[16].first->Write();
00273 
00274 
00275 // Ascii printout (need to modify to include new info)
00276     HcalDetId detid = _meol->first;
00277 
00278     if (m_outputFileText!=""){
00279       if(m_fitflag==0) {
00280         m_outFile<<detid<<"   "<<time1<<" "<<dtime1<<std::endl;
00281         m_outputFileXML << "  <DATA_SET>" << endl;
00282         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
00283         m_outputFileXML << "    <CHANNEL>" << endl;
00284         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
00285         sprintf(output, "      <ETA>%2i</ETA>", detid.ietaAbs() );
00286         m_outputFileXML << output << endl;
00287         sprintf(output, "      <PHI>%2i</PHI>", detid.iphi() );
00288         m_outputFileXML << output << endl;
00289         sprintf(output, "      <DEPTH>%2i</DEPTH>", detid.depth() );
00290         m_outputFileXML << output << endl;
00291         sprintf(output, "      <Z>%2i</Z>", detid.zside() );
00292         m_outputFileXML << output << endl;
00293 
00294         if(detid.subdet() == 1)  m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
00295         if(detid.subdet() == 2)  m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
00296         if(detid.subdet() == 3)  m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
00297         if(detid.subdet() == 4)  m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
00298         sprintf(output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId() );
00299         m_outputFileXML << output << endl;
00300         m_outputFileXML << "    </CHANNEL>" << endl;
00301         m_outputFileXML << "    <DATA>" << endl;
00302         sprintf(output, "      <MEAN_TIME>%7f</MEAN_TIME>", time1);
00303         m_outputFileXML << output << endl;
00304         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
00305         sprintf(output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
00306         m_outputFileXML << output << endl;
00307         sprintf(output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag+1);
00308         m_outputFileXML << output << endl;
00309         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
00310         m_outputFileXML << "    </DATA>" << endl;
00311         m_outputFileXML << "  </DATA_SET>" << endl;
00312 
00313         }
00314       else if(m_fitflag==1){
00315         m_outFile<<detid<<"   "<<time2<<" "<<dtime2<<std::endl;
00316         m_outputFileXML << "  <DATA_SET>" << endl;
00317         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
00318         m_outputFileXML << "    <CHANNEL>" << endl;
00319         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
00320         sprintf(output, "      <ETA>%2i</ETA>", detid.ietaAbs() );
00321         m_outputFileXML << output << endl;
00322         sprintf(output, "      <PHI>%2i</PHI>", detid.iphi() );
00323         m_outputFileXML << output << endl;
00324         sprintf(output, "      <DEPTH>%2i</DEPTH>", detid.depth() );
00325         m_outputFileXML << output << endl;
00326         sprintf(output, "      <Z>%2i</Z>", detid.zside() );
00327         m_outputFileXML << output << endl;
00328         if(detid.subdet() == 1)  m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>"<< endl;
00329         if(detid.subdet() == 2)  m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>"<< endl;
00330         if(detid.subdet() == 3)  m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>"<< endl;
00331         if(detid.subdet() == 4)  m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>"<< endl;
00332         sprintf(output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId() );
00333         m_outputFileXML << output << endl;
00334         m_outputFileXML << "    </CHANNEL>" << endl;
00335         m_outputFileXML << "    <DATA>" << endl;
00336         sprintf(output, "      <MEAN_TIME>%7f</MEAN_TIME>", time2);
00337         m_outputFileXML << output << endl;
00338         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
00339         sprintf(output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime2);
00340         m_outputFileXML << output << endl;
00341         sprintf(output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag+1);
00342         m_outputFileXML << output << endl;
00343         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
00344         m_outputFileXML << "    </DATA>" << endl;
00345         m_outputFileXML << "  </DATA_SET>" << endl;
00346         }
00347 
00348       else if(m_fitflag==2){
00349         m_outFile<<detid<<"   "<<time3<<" "<<dtime3<<std::endl;
00350         m_outputFileXML << "  <DATA_SET>" << endl;
00351         m_outputFileXML << "    <VERSION>version:1</VERSION>" << endl;
00352         m_outputFileXML << "    <CHANNEL>" << endl;
00353         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
00354         sprintf(output, "      <ETA>%2i</ETA>", detid.ietaAbs() );
00355         m_outputFileXML << output << endl;
00356         sprintf(output, "      <PHI>%2i</PHI>", detid.iphi() );
00357         m_outputFileXML << output << endl;
00358         sprintf(output, "      <DEPTH>%2i</DEPTH>", detid.depth() );
00359         m_outputFileXML << output << endl;
00360         sprintf(output, "      <Z>%2i</Z>", detid.zside() );
00361         m_outputFileXML << output << endl;
00362         if(detid.subdet() == 1) m_outputFileXML << "      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
00363         if(detid.subdet() == 2) m_outputFileXML << "      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
00364         if(detid.subdet() == 3) m_outputFileXML << "      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
00365         if(detid.subdet() == 4) m_outputFileXML << "      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
00366         sprintf(output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId() );
00367         m_outputFileXML << output << endl;
00368         m_outputFileXML << "    </CHANNEL>" << endl;
00369         m_outputFileXML << "    <DATA>" << endl;
00370         sprintf(output, "      <MEAN_TIME>%7f</MEAN_TIME>", time3);
00371         m_outputFileXML << output << endl;
00372         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
00373         sprintf(output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime3);
00374         m_outputFileXML << output << endl;
00375         sprintf(output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag+1);
00376         m_outputFileXML << output << endl;
00377         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
00378         m_outputFileXML << "    </DATA>" << endl;
00379         m_outputFileXML << "  </DATA_SET>" << endl;
00380         }
00381       else if(m_fitflag==3){
00382         m_outFile<<detid<<"   "<<time4<<" "<<dtime4<<std::endl;
00383         m_outputFileXML << "  <DATA_SET>" << endl;
00384         m_outputFileXML <<"    <VERSION>version:1</VERSION>" << endl;
00385         m_outputFileXML << "    <CHANNEL>" << endl;
00386         m_outputFileXML << "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
00387         sprintf(output, "      <ETA>%2i</ETA>", detid.ietaAbs() );
00388         m_outputFileXML << output << endl;
00389         sprintf(output, "      <PHI>%2i</PHI>", detid.iphi() );
00390         m_outputFileXML << output << endl;
00391         sprintf(output, "      <DEPTH>%2i</DEPTH>", detid.depth() );
00392         m_outputFileXML << output << endl;
00393         sprintf(output, "      <Z>%2i</Z>", detid.zside() );
00394         m_outputFileXML << output << endl;
00395         if(detid.subdet() == 1) m_outputFileXML <<"      <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
00396         if(detid.subdet() == 2) m_outputFileXML <<"      <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
00397         if(detid.subdet() == 3) m_outputFileXML <<"      <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
00398         if(detid.subdet() == 4) m_outputFileXML <<"      <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
00399         sprintf(output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId() );
00400         m_outputFileXML << output << endl;
00401         m_outputFileXML << "    </CHANNEL>" << endl;
00402         m_outputFileXML << "    <DATA>" << endl;
00403         sprintf(output, "      <MEAN_TIME>%7f</MEAN_TIME>", time4);
00404         m_outputFileXML << output << endl;
00405         m_outputFileXML << "      <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
00406         sprintf(output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
00407         m_outputFileXML << output << endl;
00408         sprintf(output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag+1);
00409         m_outputFileXML << output << endl;
00410         m_outputFileXML << "      <STATUS_WORD>  0</STATUS_WORD>" << endl;
00411         m_outputFileXML << "    </DATA>" << endl;
00412         m_outputFileXML << "  </DATA_SET>" << endl;
00413         }
00414 
00415       else if(m_fitflag==4){
00416         m_outFile<<detid<<"   "<<time1<<" "<<dtime1<<"   "<<time2<<" "<<dtime2<<"   "<<time3<<" "<<dtime3<<"   "<<time4<<" "<<dtime4<<std::endl;
00417         }
00418     }
00419   }
00420 }
00421 
00422 //-----------------------------------------------------------------------------
00423 void HcalLedAnalysis::LedSampleAnalysis(){
00424   // it is called every m_nevtsample events (a sample) and the end of run
00425   char LedSampleNum[20];
00426 
00427   sprintf(LedSampleNum,"LedSample_%d",sample);
00428   m_file->cd();
00429   m_file->mkdir(LedSampleNum);
00430   m_file->cd(LedSampleNum);
00431 
00432 // Compute LED constants for each HB/HE, HO, HF
00433   GetLedConst(hbHists.LEDTRENDS);
00434   GetLedConst(hoHists.LEDTRENDS);
00435   GetLedConst(hfHists.LEDTRENDS);
00436 }
00437 
00438 //-----------------------------------------------------------------------------
00439 void HcalLedAnalysis::LedTrendings(map<HcalDetId, map<int,LEDBUNCH> > &toolT)
00440 {
00441 
00442   for(_meol=toolT.begin(); _meol!=toolT.end(); _meol++){
00443     char name[1024];
00444     HcalDetId detid = _meol->first;
00445     sprintf(name,"LED timing trend, eta=%d phi=%d depth=%d",detid.ieta(),detid.iphi(),detid.depth());
00446     int bins = _meol->second[10+m_fitflag].second.first[0].size();
00447     float lo =0.5;
00448     float hi = (float)bins+0.5;
00449     _meol->second[10+m_fitflag].second.second.push_back(new TH1F(name,name,bins,lo,hi));
00450 
00451     std::vector<double>::iterator sample_it;
00452 // LED timing - put content and errors
00453     int j=0;
00454     for(sample_it=_meol->second[10+m_fitflag].second.first[0].begin();
00455         sample_it!=_meol->second[10+m_fitflag].second.first[0].end();sample_it++){
00456       _meol->second[10+m_fitflag].second.second[0]->SetBinContent(++j,*sample_it);
00457     }
00458     j=0;
00459     for(sample_it=_meol->second[10+m_fitflag].second.first[1].begin();
00460         sample_it!=_meol->second[10+m_fitflag].second.first[1].end();sample_it++){
00461       _meol->second[10+m_fitflag].second.second[0]->SetBinError(++j,*sample_it);
00462     }
00463     sprintf(name,"Sample (%d events)",m_nevtsample);
00464     _meol->second[10+m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
00465     _meol->second[10+m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
00466     _meol->second[10+m_fitflag].second.second[0]->Write();
00467   }
00468 }
00469 
00470 //-----------------------------------------------------------------------------
00471 void HcalLedAnalysis::LedDone() 
00472 {
00473 
00474 // First process the last sample (remaining events).
00475   if(evt%m_nevtsample!=0) LedSampleAnalysis();
00476 
00477 // Now do the end of run analysis: trending histos
00478   if(sample>1 && m_fitflag!=4){
00479     m_file->cd();
00480     m_file->cd("HBHE");
00481     LedTrendings(hbHists.LEDTRENDS);
00482     m_file->cd();
00483     m_file->cd("HO");
00484     LedTrendings(hoHists.LEDTRENDS);
00485     m_file->cd();
00486     m_file->cd("HF");
00487     LedTrendings(hfHists.LEDTRENDS);
00488   }
00489 
00490   // Write other histograms.
00491   // HB
00492   m_file->cd();
00493   m_file->cd("HBHE");
00494   hbHists.ALLLEDS->Write();
00495   hbHists.LEDRMS->Write();
00496   hbHists.LEDMEAN->Write();
00497   // HO
00498   m_file->cd();
00499   m_file->cd("HO");
00500   hoHists.ALLLEDS->Write();
00501   hoHists.LEDRMS->Write();
00502   hoHists.LEDMEAN->Write();
00503   // HF
00504   m_file->cd();
00505   m_file->cd("HF");
00506   hfHists.ALLLEDS->Write();
00507   hfHists.LEDRMS->Write();
00508   hfHists.LEDMEAN->Write();
00509   // Calib
00510   m_file->cd();
00511   m_file->cd("Calib");
00512   for(_meca=calibHists.begin(); _meca!=calibHists.end(); _meca++){
00513     _meca->second.avePulse->Write();
00514     _meca->second.integPulse->Write();
00515   }
00516 
00517   // Write the histo file and close it
00518 //  m_file->Write();
00519   m_file->Close();
00520   cout << "Hcal histograms written to " << m_outputFileROOT.c_str() << endl;
00521 }
00522 
00523 //-----------------------------------------------------------------------------
00524 void HcalLedAnalysis::processLedEvent(const HBHEDigiCollection& hbhe,
00525                                         const HODigiCollection& ho,
00526                                         const HFDigiCollection& hf,
00527                                         const HcalCalibDigiCollection calib,
00528                                         const HcalDbService& cond)
00529 {
00530   evt++;
00531   sample = (evt-1)/m_nevtsample +1;
00532   evt_curr = evt%m_nevtsample;
00533   if(evt_curr==0)evt_curr=m_nevtsample;
00534 
00535   // Calib
00536 
00537   if (m_usecalib){
00538     try{
00539       if(!calib.size()) throw (int)calib.size();
00540       // this is effectively a loop over electronic channels
00541       for (HcalCalibDigiCollection::const_iterator j=calib.begin(); j!=calib.end(); j++){
00542         const HcalCalibDataFrame digi = (const HcalCalibDataFrame)(*j);   
00543         HcalElectronicsId elecId = digi.elecId();
00544         HcalCalibDetId calibId = digi.id();
00545         ProcessCalibEvent(elecId.fiberChanId(),calibId,digi);  //Shouldn't depend on anything in elecId but not sure how else to do it 
00546       }
00547     }
00548     catch (int i ) {
00549     //  m_logFile<< "Event with " << i<<" Calib Digis passed." << std::endl;
00550     }
00551   }
00552 
00553 
00554   // HB + HE
00555   try{
00556     if(!hbhe.size()) throw (int)hbhe.size();
00557 // this is effectively a loop over electronic channels
00558     for (HBHEDigiCollection::const_iterator j=hbhe.begin(); j!=hbhe.end(); j++){
00559       const HBHEDataFrame digi = (const HBHEDataFrame)(*j);
00560       for(int k=0; k<(int)state.size();k++) state[k]=true;
00561       // See if histos exist for this channel, and if not, create them
00562       _meol = hbHists.LEDTRENDS.find(digi.id());
00563       if (_meol==hbHists.LEDTRENDS.end()){
00564         SetupLEDHists(0,digi.id(),hbHists.LEDTRENDS);
00565       }
00566       LedHBHEHists(digi.id(),digi,hbHists.LEDTRENDS,cond);
00567     }
00568   }
00569   catch (int i ) {
00570 //    m_logFile<< "Event with " << i<<" HBHE Digis passed." << std::endl;
00571   } 
00572 
00573   // HO
00574   try{
00575     if(!ho.size()) throw (int)ho.size();
00576     for (HODigiCollection::const_iterator j=ho.begin(); j!=ho.end(); j++){
00577       const HODataFrame digi = (const HODataFrame)(*j);
00578       _meol = hoHists.LEDTRENDS.find(digi.id());
00579       if (_meol==hoHists.LEDTRENDS.end()){
00580         SetupLEDHists(1,digi.id(),hoHists.LEDTRENDS);
00581       }
00582       LedHOHists(digi.id(),digi,hoHists.LEDTRENDS,cond);
00583     }        
00584   } 
00585   catch (int i ) {
00586 //    m_logFile << "Event with " << i<<" HO Digis passed." << std::endl;
00587   } 
00588 
00589   // HF
00590   try{
00591     if(!hf.size()) throw (int)hf.size();
00592     for (HFDigiCollection::const_iterator j=hf.begin(); j!=hf.end(); j++){
00593       const HFDataFrame digi = (const HFDataFrame)(*j);
00594       _meol = hfHists.LEDTRENDS.find(digi.id());
00595       if (_meol==hfHists.LEDTRENDS.end()){
00596         SetupLEDHists(2,digi.id(),hfHists.LEDTRENDS);
00597       }
00598       LedHFHists(digi.id(),digi,hfHists.LEDTRENDS,cond);
00599     }
00600   } 
00601   catch (int i ) {
00602 //    m_logFile << "Event with " << i<<" HF Digis passed." << std::endl;
00603   } 
00604 
00605   // Call the function every m_nevtsample events
00606   if(evt%m_nevtsample==0) LedSampleAnalysis();
00607 
00608 }
00609 //----------------------------------------------------------------------------
00610 void HcalLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int,LEDBUNCH> > &toolT) {
00611 
00612   string type = "HBHE";
00613   if(id==1) type = "HO";
00614   if(id==2) type = "HF";
00615 
00616   _meol = toolT.find(detid);
00617   if (_meol==toolT.end()){
00618 // if histos for this channel do not exist, create them
00619     map<int,LEDBUNCH> insert;
00620     char name[1024];
00621     for(int i=0; i<10; i++){
00622       sprintf(name,"%s Pulse height, eta=%d phi=%d depth=%d TS=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth(),i);
00623       insert[i].first =  new TH1F(name,name,200,0.,2000.);
00624     }
00625     sprintf(name,"%s LED Mean pulse, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00626     insert[10].first =  new TH1F(name,name,10,-0.5,9.5);
00627     sprintf(name,"%s LED Pulse, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00628     insert[11].first =  new TH1F(name,name,10,-0.5,9.5);
00629     sprintf(name,"%s Mean TS, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00630     insert[12].first =  new TH1F(name,name,200,0.,10.);
00631     sprintf(name,"%s Peak TS, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00632     insert[13].first =  new TH1F(name,name,200,0.,10.);
00633     sprintf(name,"%s Peak TS error, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00634     insert[14].first =  new TH1F(name,name,200,0.,0.05);
00635     sprintf(name,"%s Fit chi2, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00636     insert[15].first =  new TH1F(name,name,100,0.,50.);
00637     sprintf(name,"%s Integrated Signal, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00638     insert[16].first =  new TH1F(name,name,500,0.,5000.);
00639 
00640     toolT[detid] = insert;
00641   }
00642 }
00643 //-----------------------------------------------------------------------------
00644 void HcalLedAnalysis::LedHBHEHists(const HcalDetId& detid, const HBHEDataFrame& ledDigi, map<HcalDetId, map<int,LEDBUNCH> > &toolT, const HcalDbService& cond){
00645 
00646   map<int,LEDBUNCH> _mei;
00647   _meol = toolT.find(detid);
00648   _mei = _meol->second;
00649 
00650   // Reset the histos if we're at the end of a 'bunch'
00651   if((evt-1)%m_nevtsample==0 && state[0]){
00652     for(int k=0; k<(int)state.size();k++) state[k]=false;
00653     for(int i=0; i<16; i++) _mei[i].first->Reset();
00654   }
00655 
00656   // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
00657 
00658 
00659   //  int maxTS = -1;
00660   float max_fC = 0;
00661   float ta = 0;
00662   m_coder = cond.getHcalCoder(detid);
00663   m_ped = cond.getPedestal(detid);
00664   m_shape = cond.getHcalShape(m_coder);
00665   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++){
00666     int capid = ledDigi[TS].capid();
00667     int adc = ledDigi[TS].adc();
00668     double fC = m_coder->charge(*m_shape,adc,capid);
00669     ta = (fC - m_ped->getValue(capid));
00670     //cout << "DetID: " << detid << "  CapID: " << capid << "  ADC: " << adc << "  fC: " << fC << endl;
00671     _mei[TS].first->Fill(ta);
00672     _mei[10].first->AddBinContent(TS+1,ta);  // This is average pulse, could probably do better (Profile?)
00673     if(m_fitflag>1){
00674       if(TS==m_startTS)_mei[11].first->Reset();
00675       _mei[11].first->SetBinContent(TS+1,ta);
00676     }
00677     // keep track of max TS and max amplitude (in fC)
00678     if (ta > max_fC){
00679       max_fC = ta;
00680       //      maxTS = TS;
00681     }
00682   }
00683 
00684   // Now we have a sample with pedestals subtracted and in units of fC
00685   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
00686   // we now want to use Phil's timing correction.  This is not necessary
00687   // if we are performing a Landau fit (m_fitflag = 3)
00688 
00689   float sum=0.;
00690   for(int i=0; i<10; i++)sum=sum+_mei[11].first->GetBinContent(i+1);
00691   if(sum>100){
00692     if(m_fitflag==2 || m_fitflag==4){
00693       float timmean=_mei[11].first->GetMean();  // let's use Phil's way instead
00694       float timmeancorr=BinsizeCorr(timmean);
00695       _mei[12].first->Fill(timmeancorr);
00696     }
00697     _mei[16].first->Fill(_mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
00698     if(m_fitflag==3 || m_fitflag==4){
00699       _mei[11].first->Fit("landau","Q");
00700       TF1 *fit = _mei[11].first->GetFunction("landau");
00701       _mei[13].first->Fill(fit->GetParameter(1));
00702       _mei[14].first->Fill(fit->GetParError(1));
00703       _mei[15].first->Fill(fit->GetChisquare()/fit->GetNDF());
00704     }
00705   }
00706 
00707 }
00708 
00709 //-----------------------------------------------------------------------------
00710 void HcalLedAnalysis::LedHOHists(const HcalDetId& detid, const HODataFrame& ledDigi, map<HcalDetId, map<int,LEDBUNCH> > &toolT, const HcalDbService& cond) {
00711 
00712   map<int,LEDBUNCH> _mei;
00713   _meol = toolT.find(detid);
00714   _mei = _meol->second;
00715   // Rest the histos if we're at the end of a 'bunch'
00716   if((evt-1)%m_nevtsample==0 && state[0]){
00717     for(int k=0; k<(int)state.size();k++) state[k]=false;
00718     for(int i=0; i<16; i++) _mei[i].first->Reset();
00719   }
00720 
00721   // now we have the signal in fC, let's get rid of that darn pedestal
00722   // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
00723 
00724   //  int maxTS = -1;
00725   float max_fC = 0;
00726   float ta = 0;
00727   m_coder = cond.getHcalCoder(detid);
00728   m_ped = cond.getPedestal(detid);
00729   m_shape = cond.getHcalShape(m_coder);
00730   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++){
00731     int capid = ledDigi[TS].capid();
00732     int adc = ledDigi[TS].adc();
00733     double fC = m_coder->charge(*m_shape,adc,capid);
00734     ta = (fC - m_ped->getValue(capid));
00735     _mei[TS].first->Fill(ta);
00736     _mei[10].first->AddBinContent(TS+1,ta);  // This is average pulse, could probably do better (Profile?)
00737     if(m_fitflag>1){
00738       if(TS==m_startTS)_mei[11].first->Reset();
00739       _mei[11].first->SetBinContent(TS+1,ta);
00740     }
00741     // keep track of max TS and max amplitude (in fC)
00742     if (ta > max_fC){
00743       max_fC = ta;
00744       //      maxTS = TS;
00745     }
00746   }
00747 
00748   // Now we have a sample with pedestals subtracted and in units of fC
00749   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
00750   // we now want to use Phil's timing correction.  This is not necessary
00751   // if we are performing a Landau fit (m_fitflag = 3)
00752 
00753   float sum=0.;
00754   for(int i=0; i<10; i++)sum=sum+_mei[11].first->GetBinContent(i+1);
00755   if(sum>100){
00756     if(m_fitflag==2 || m_fitflag==4){
00757       float timmean=_mei[11].first->GetMean();  // let's use Phil's way instead
00758       float timmeancorr=BinsizeCorr(timmean);
00759       _mei[12].first->Fill(timmeancorr);
00760     }
00761     _mei[16].first->Fill(_mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
00762     if(m_fitflag==3 || m_fitflag==4){
00763       _mei[11].first->Fit("landau","Q");
00764       TF1 *fit = _mei[11].first->GetFunction("landau");
00765       _mei[13].first->Fill(fit->GetParameter(1));
00766       _mei[14].first->Fill(fit->GetParError(1));
00767       _mei[15].first->Fill(fit->GetChisquare()/fit->GetNDF());
00768     }
00769   }
00770 
00771 }
00772 
00773 //-----------------------------------------------------------------------------
00774 void HcalLedAnalysis::LedHFHists(const HcalDetId& detid, const HFDataFrame& ledDigi, map<HcalDetId, map<int,LEDBUNCH> > &toolT, const HcalDbService& cond) {
00775 
00776   map<int,LEDBUNCH> _mei;
00777   _meol = toolT.find(detid);
00778   _mei = _meol->second;
00779   // Rest the histos if we're at the end of a 'bunch'
00780   if((evt-1)%m_nevtsample==0 && state[0]){
00781     for(int k=0; k<(int)state.size();k++) state[k]=false;
00782     for(int i=0; i<16; i++) _mei[i].first->Reset();
00783   }
00784 
00785   // now we have the signal in fC, let's get rid of that darn pedestal
00786   // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
00787 
00788   //  int maxTS = -1;
00789   float max_fC = 0;
00790   float ta = 0;
00791   m_coder = cond.getHcalCoder(detid);
00792   m_ped = cond.getPedestal(detid);
00793   m_shape = cond.getHcalShape(m_coder);
00794   //cout << "New Digi!!!!!!!!!!!!!!!!!!!!!!" << endl;
00795   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++){
00796     int capid = ledDigi[TS].capid();
00797     // BE CAREFUL: this is assuming peds are stored in ADCs
00798     int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
00799     if (adc < 0){ adc = 0; }  // to prevent negative adcs after ped subtraction, which should really only happen
00800                               // if you're using the wrong peds.
00801     double fC = m_coder->charge(*m_shape,adc,capid);
00802     //ta = (fC - m_ped->getValue(capid));
00803     ta = fC;
00804     //cout << "DetID: " << detid << "  CapID: " << capid << "  ADC: " << adc << "  Ped: " << m_ped->getValue(capid) << "  fC: " << fC << endl;
00805     _mei[TS].first->Fill(ta);
00806     _mei[10].first->AddBinContent(TS+1,ta);  // This is average pulse, could probably do better (Profile?)
00807     if(m_fitflag>1){
00808       if(TS==m_startTS)_mei[11].first->Reset();
00809       _mei[11].first->SetBinContent(TS+1,ta);
00810     }
00811 
00812     // keep track of max TS and max amplitude (in fC)
00813     if (ta > max_fC){
00814       max_fC = ta;
00815       //      maxTS = TS;
00816     }
00817   }
00818 
00819   // Now we have a sample with pedestals subtracted and in units of fC
00820   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
00821   // we now want to use Phil's timing correction.  This is not necessary
00822   // if we are performing a Landau fit (m_fitflag = 3)
00823 
00824   float sum=0.;
00825   for(int i=0; i<10; i++)sum=sum+_mei[11].first->GetBinContent(i+1);
00826   if(sum>100){
00827     if(m_fitflag==2 || m_fitflag==4){
00828       float timmean=_mei[11].first->GetMean();  // let's use Phil's way instead
00829       float timmeancorr=BinsizeCorr(timmean);
00830       _mei[12].first->Fill(timmeancorr);
00831     }
00832     _mei[16].first->Fill(_mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
00833     if(m_fitflag==3 || m_fitflag==4){
00834       _mei[11].first->Fit("landau","Q");
00835       TF1 *fit = _mei[11].first->GetFunction("landau");
00836       _mei[13].first->Fill(fit->GetParameter(1));
00837       _mei[14].first->Fill(fit->GetParError(1));
00838       _mei[15].first->Fill(fit->GetChisquare()/fit->GetNDF());
00839     }
00840   }
00841 
00842 
00843 
00844 }
00845 
00846 
00847 //-----------------------------------------------------------------------------
00848 float HcalLedAnalysis::BinsizeCorr(float time) {
00849 
00850 // this is the bin size correction to be applied for laser data (from Andy),
00851 // it comes from a pulse shape measured from TB04 data (from Jordan)
00852 // This should eventually be replaced with the more thorough treatment from Phil
00853 
00854   float corrtime=0.;
00855   static const float tstrue[32]={0.003, 0.03425, 0.06548, 0.09675, 0.128,
00856  0.15925, 0.1905, 0.22175, 0.253, 0.28425, 0.3155, 0.34675, 0.378, 0.40925,
00857  0.4405, 0.47175, 0.503, 0.53425, 0.5655, 0.59675, 0.628, 0.65925, 0.6905,
00858  0.72175, 0.753, 0.78425, 0.8155, 0.84675, 0.878, 0.90925, 0.9405, 0.97175};
00859   static const float tsreco[32]={-0.00422, 0.01815, 0.04409, 0.07346, 0.09799,
00860  0.12192, 0.15072, 0.18158, 0.21397, 0.24865, 0.28448, 0.31973, 0.35449,
00861  0.39208, 0.43282, 0.47244, 0.5105, 0.55008, 0.58827, 0.62828, 0.6717, 0.70966,
00862  0.74086, 0.77496, 0.80843, 0.83472, 0.86044, 0.8843, 0.90674, 0.92982,
00863  0.95072, 0.9726};
00864 
00865  int inttime=(int)time;
00866  float restime=time-inttime;
00867  for(int i=0; i<=32; i++) {
00868    float lolim=0.; float uplim=1.; float tsdown; float tsup;
00869    if(i>0){
00870      lolim=tsreco[i-1];
00871      tsdown=tstrue[i-1];
00872    }
00873    else tsdown=tstrue[31]-1.;
00874    if(i<32){
00875      uplim=tsreco[i];
00876      tsup=tstrue[i];
00877    }
00878    else tsup=tstrue[0]+1.;
00879    if(restime>=lolim && restime<uplim){
00880       corrtime=(tsdown*(uplim-restime)+tsup*(restime-lolim)) / (uplim-lolim);
00881     }
00882   }
00883   corrtime+=inttime;
00884 
00885  return corrtime;
00886 }
00887 //-----------------------------------------------------------------------------
00888 
00889 // Will try to implement Phil's time slew correction here at some point
00890 
00891 
00892 //-----------------------------------------------------------------------------
00893 void HcalLedAnalysis::ProcessCalibEvent(int fiberChan, HcalCalibDetId calibId, const HcalCalibDataFrame digi){
00894 
00895   _meca = calibHists.find(calibId);
00896   if (_meca==calibHists.end()){
00897   // if histos for this channel do not exist, first create them
00898     char name[1024];
00899     std::string prefix;
00900     if (calibId.calibFlavor()==HcalCalibDetId::CalibrationBox) {
00901       std::string sector=(calibId.hcalSubdet()==HcalBarrel)?("HB"):
00902         (calibId.hcalSubdet()==HcalEndcap)?("HE"):
00903         (calibId.hcalSubdet()==HcalOuter)?("HO"):
00904         (calibId.hcalSubdet()==HcalForward)?("HF"):"";
00905       sprintf(name,"%s %+d iphi=%d %s",sector.c_str(),calibId.ieta(),calibId.iphi(),calibId.cboxChannelString().c_str());
00906       prefix=name;
00907     }
00908     
00909     sprintf(name,"%s Pin Diode Mean",prefix.c_str());
00910     calibHists[calibId].avePulse = new TProfile(name,name,10,-0.5,9.5,0,1000);
00911     sprintf(name,"%s Pin Diode Current Pulse",prefix.c_str());
00912     calibHists[calibId].thisPulse = new TH1F(name,name,10,-0.5,9.5);
00913     sprintf(name,"%s Pin Diode Integrated Pulse",prefix.c_str());
00914     calibHists[calibId].integPulse = new TH1F(name,name,200,0,500);    
00915   }
00916   else {
00917     for (int i=m_startTS; i<digi.size() && i<=m_endTS; i++) {
00918       calibHists[calibId].avePulse->Fill(i,digi.sample(i).adc());
00919       calibHists[calibId].thisPulse->SetBinContent(i+1,digi.sample(i).adc());
00920     }
00921     calibHists[calibId].integPulse->Fill(calibHists[calibId].thisPulse->Integral());
00922   }
00923 }
00924
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/CalibCalorimetry/HcalAlgos/src/HcalLedAnalysis.cc
