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