CMSSW Reference Manual

00001 
00002 #include "CalibFormats/CastorObjects/interface/CastorCoderDb.h"
00003 #include "CalibFormats/CastorObjects/interface/CastorCalibrations.h"
00004 #include "CalibFormats/CastorObjects/interface/CastorDbService.h"
00005 #include "CalibFormats/CastorObjects/interface/CastorDbRecord.h"
00006 #include "CondFormats/CastorObjects/interface/CastorQIECoder.h"
00007 #include "CondFormats/CastorObjects/interface/CastorPedestals.h"
00008 #include "CondFormats/CastorObjects/interface/CastorPedestalWidths.h"
00009 
00010 #include "CalibCalorimetry/CastorCalib/interface/CastorLedAnalysis.h"
00011 #include "CalibCalorimetry/CastorCalib/interface/CastorAlgoUtils.h"
00012 #include <TFile.h>
00013 #include <math.h>
00014 
00015 using namespace std;
00016 
00017 
00018 CastorLedAnalysis::CastorLedAnalysis(const edm::ParameterSet& ps)
00019 {
00020   // init
00021   evt=0;
00022   sample=0;
00023   m_file=0;
00024   // output files
00025   for(int k=0;k<4;k++) state.push_back(true); // 4 cap-ids (do we care?)
00026   m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
00027   m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML","");
00028   if ( m_outputFileText.size() != 0 ) {
00029     cout << "Castor LED results will be saved to " << m_outputFileText.c_str() << endl;
00030     m_outFile.open(m_outputFileText.c_str());
00031   } 
00032   m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
00033   if ( m_outputFileROOT.size() != 0 ) {
00034     cout << "Castor LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
00035   }
00036 
00037   m_nevtsample = ps.getUntrackedParameter<int>("nevtsample",9999999);
00038   if(m_nevtsample<1)m_nevtsample=9999999;
00039   m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag",0);
00040   if(m_hiSaveflag<0)m_hiSaveflag=0;
00041   if(m_hiSaveflag>0)m_hiSaveflag=1;
00042   m_fitflag = ps.getUntrackedParameter<int>("analysisflag",2);
00043   if(m_fitflag<0)m_fitflag=0;
00044   if(m_fitflag>4)m_fitflag=4;
00045   m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
00046   if(m_startTS<0) m_startTS=0;
00047   m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
00048   m_usecalib = ps.getUntrackedParameter<bool>("usecalib",false);
00049   m_logFile.open("CastorLedAnalysis.log");
00050 
00051   int runNum = ps.getUntrackedParameter<int>("runNumber",999999);
00052 
00053   // histogram booking
00054   castorHists.ALLLEDS = new TH1F("Castor All LEDs","HF All Leds",10,0,9);
00055   castorHists.LEDRMS= new TH1F("Castor All LED RMS","HF All LED RMS",100,0,3);
00056   castorHists.LEDMEAN= new TH1F("Castor All LED Means","HF All LED Means",100,0,9);
00057   castorHists.CHI2= new TH1F("Castor Chi2 by ndf for Landau fit","HF Chi2/ndf Landau",200,0.,50.);
00058 
00059   //XML file header
00060   m_outputFileXML.open(m_outputFileX.c_str());
00061   sprintf(output, "<?xml version='1.0' encoding='UTF-8'?>");
00062   m_outputFileXML << output << endl;
00063   sprintf(output, "<ROOT>");
00064   m_outputFileXML << output << endl << endl;
00065   sprintf(output, "  <HEADER>");
00066   m_outputFileXML << output << endl;
00067   sprintf(output, "    <TYPE>");
00068   m_outputFileXML << output << endl;
00069   sprintf(output, "      <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>");
00070   m_outputFileXML << output << endl;
00071   sprintf(output, "      <NAME>HCAL LED Timing</NAME>");
00072   m_outputFileXML << output << endl;
00073   sprintf(output, "    </TYPE>");
00074   m_outputFileXML << output << endl;
00075   sprintf(output, "    <RUN>");
00076   m_outputFileXML << output << endl;
00077   sprintf(output, "      <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>");
00078   m_outputFileXML << output << endl;
00079   sprintf(output, "      <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
00080   m_outputFileXML << output << endl;
00081   sprintf(output, "      <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>");
00082   m_outputFileXML << output << endl;
00083   sprintf(output, "      <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>");
00084   m_outputFileXML << output << endl;
00085   sprintf(output, "    </RUN>");
00086   m_outputFileXML << output << endl;
00087   sprintf(output, "  </HEADER>");
00088   m_outputFileXML << output << endl;
00089   sprintf(output, "<!-- Tags secton -->");
00090   m_outputFileXML << output << endl;
00091   sprintf(output, "  <ELEMENTS>");
00092   m_outputFileXML << output << endl;
00093   sprintf(output, "    <DATA_SET id='-1'/>");
00094   m_outputFileXML << output << endl;
00095   sprintf(output, "      <IOV id='1'>");
00096   m_outputFileXML << output << endl;
00097   sprintf(output, "        <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>");
00098   m_outputFileXML << output << endl;
00099   sprintf(output, "        <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>");
00100   m_outputFileXML << output << endl;
00101   sprintf(output, "      </IOV>");
00102   m_outputFileXML << output << endl;
00103   sprintf(output, "      <TAG id='2' mode='auto'>");
00104   m_outputFileXML << output << endl;
00105   sprintf(output, "        <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
00106   m_outputFileXML << output << endl;
00107   sprintf(output, "        <DETECTOR_NAME>HCAL</DETECTOR_NAME>");
00108   m_outputFileXML << output << endl;
00109   sprintf(output, "        <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>");
00110   m_outputFileXML << output << endl;
00111   sprintf(output, "      </TAG>");
00112   m_outputFileXML << output << endl;
00113   sprintf(output, "  </ELEMENTS>");
00114   m_outputFileXML << output << endl;
00115   sprintf(output, "  <MAPS>");
00116   m_outputFileXML << output << endl;
00117   sprintf(output, "      <TAG idref ='2'>");
00118   m_outputFileXML << output << endl;
00119   sprintf(output, "        <IOV idref='1'>");
00120   m_outputFileXML << output << endl;
00121   sprintf(output, "          <DATA_SET idref='-1' />");
00122   m_outputFileXML << output << endl;
00123   sprintf(output, "        </IOV>");
00124   m_outputFileXML << output << endl;
00125   sprintf(output, "      </TAG>");
00126   m_outputFileXML << output << endl;
00127   sprintf(output, "  </MAPS>");
00128   m_outputFileXML << output << endl;
00129 
00130 
00131 }
00132 
00133 //-----------------------------------------------------------------------------
00134 CastorLedAnalysis::~CastorLedAnalysis(){
00136 
00137   for(_meol=castorHists.LEDTRENDS.begin(); _meol!=castorHists.LEDTRENDS.end(); _meol++){
00138     for(int i=0; i<15; i++) _meol->second[i].first->Delete();
00139   }
00140 
00141   castorHists.ALLLEDS->Delete();
00142   castorHists.LEDRMS->Delete();
00143   castorHists.LEDMEAN->Delete();
00144   castorHists.CHI2->Delete();
00145 }
00146 
00147 //-----------------------------------------------------------------------------
00148 void CastorLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
00149   // open the histogram file, create directories within
00150   m_file=new TFile(m_outputFileROOT.c_str(),"RECREATE");
00151   m_file->mkdir("Castor");
00152   m_file->cd();
00153 }
00154 
00155 //-----------------------------------------------------------------------------
00156 void CastorLedAnalysis::GetLedConst(map<HcalDetId, map<int,LEDBUNCH> > &toolT){
00157   double time2=0; double time1=0; double time3=0; double time4=0;
00158   double dtime2=0; double dtime1=0; double dtime3=0; double dtime4=0;
00159 
00160   if (m_outputFileText!=""){
00161     if(m_fitflag==0 || m_fitflag==2) m_outFile<<"Det Eta,Phi,D   Mean    Error"<<std::endl;
00162     else if(m_fitflag==1 || m_fitflag==3) m_outFile<<"Det Eta,Phi,D   Peak    Error"<<std::endl;
00163     else if(m_fitflag==4) m_outFile<<"Det Eta,Phi,D   Mean    Error      Peak    Error       MeanEv  Error       PeakEv  Error"<<std::endl;
00164   }
00165   for(_meol=toolT.begin(); _meol!=toolT.end(); _meol++){
00166 // scale the LED pulse to 1 event
00167     _meol->second[10].first->Scale(1./evt_curr);
00168     if(m_fitflag==0 || m_fitflag==4){
00169       time1 = _meol->second[10].first->GetMean();
00170       dtime1 = _meol->second[10].first->GetRMS()/sqrt((float)evt_curr*(m_endTS-m_startTS+1));
00171     }
00172     if(m_fitflag==1 || m_fitflag==4){
00173 // put proper errors
00174       for(int j=0; j<10; j++) _meol->second[10].first->SetBinError(j+1,_meol->second[j].first->GetRMS()/sqrt((float)evt_curr));
00175     }
00176     if(m_fitflag==1 || m_fitflag==3 || m_fitflag==4){
00177       _meol->second[10].first->Fit("landau","Q");
00178 //      _meol->second[10].first->Fit("gaus","Q");
00179       TF1 *fit = _meol->second[10].first->GetFunction("landau");
00180 //      TF1 *fit = _meol->second[10].first->GetFunction("gaus");
00181       time2=fit->GetParameter(1);
00182       dtime2=fit->GetParError(1);
00183     }
00184     if(m_fitflag==2 || m_fitflag==4){
00185       time3 = _meol->second[12].first->GetMean();
00186       dtime3 = _meol->second[12].first->GetRMS()/sqrt((float)_meol->second[12].first->GetEntries());
00187     }
00188     if(m_fitflag==3 || m_fitflag==4){
00189       time4 = _meol->second[13].first->GetMean();
00190       dtime4 = _meol->second[13].first->GetRMS()/sqrt((float)_meol->second[13].first->GetEntries());
00191     }
00192     for (int i=0; i<10; i++){
00193       _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
00194       _meol->second[i].first->GetYaxis()->SetTitle("Counts");
00195 //      if(m_hiSaveflag>0)_meol->second[i].first->Write();
00196     }
00197     _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
00198     _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
00199     if(m_hiSaveflag>0)_meol->second[10].first->Write();
00200     _meol->second[10].second.first[0].push_back(time1);
00201     _meol->second[10].second.first[1].push_back(dtime1);
00202     _meol->second[11].second.first[0].push_back(time2);
00203     _meol->second[11].second.first[1].push_back(dtime2);
00204     _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
00205     _meol->second[12].first->GetYaxis()->SetTitle("Counts");
00206     if(m_fitflag==2 && m_hiSaveflag>0)_meol->second[12].first->Write();
00207     _meol->second[12].second.first[0].push_back(time3);
00208     _meol->second[12].second.first[1].push_back(dtime3);
00209     _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
00210     _meol->second[13].first->GetYaxis()->SetTitle("Counts");
00211     if(m_fitflag>2 && m_hiSaveflag>0)_meol->second[13].first->Write();
00212     _meol->second[13].second.first[0].push_back(time4);
00213     _meol->second[13].second.first[1].push_back(dtime4);
00214     _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
00215     _meol->second[14].first->GetYaxis()->SetTitle("Counts");
00216     if(m_fitflag>2 && m_hiSaveflag>0)_meol->second[14].first->Write();
00217     _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
00218     _meol->second[15].first->GetYaxis()->SetTitle("Counts");
00219     if(m_fitflag>2 && m_hiSaveflag>0)_meol->second[15].first->Write();
00220     _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
00221     _meol->second[16].first->Write();
00222 
00223 
00224 // Ascii printout (need to modify to include new info)
00225     HcalDetId detid = _meol->first;
00226 
00227     if (m_outputFileText!=""){
00228       if(m_fitflag==0) {
00229         m_outFile<<detid<<"   "<<time1<<" "<<dtime1<<std::endl;
00230         sprintf(output, "  <DATA_SET>");
00231         m_outputFileXML << output << endl;
00232         sprintf(output, "    <VERSION>version:1</VERSION>");
00233         m_outputFileXML << output << endl;
00234         sprintf(output, "    <CHANNEL>");
00235         m_outputFileXML << output << endl;
00236         sprintf(output, "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
00237         m_outputFileXML << output << endl;
00238         sprintf(output, "      <ETA>%2i</ETA>", detid.ietaAbs() );
00239         m_outputFileXML << output << endl;
00240         sprintf(output, "      <PHI>%2i</PHI>", detid.iphi() );
00241         m_outputFileXML << output << endl;
00242         sprintf(output, "      <DEPTH>%2i</DEPTH>", detid.depth() );
00243         m_outputFileXML << output << endl;
00244         sprintf(output, "      <Z>%2i</Z>", detid.zside() );
00245         m_outputFileXML << output << endl;
00246         if(detid.subdet() == 1) sprintf(output, "      <DETECTOR_NAME>HB</DETECTOR_NAME>");
00247         if(detid.subdet() == 2) sprintf(output, "      <DETECTOR_NAME>HE</DETECTOR_NAME>");
00248         if(detid.subdet() == 3) sprintf(output, "      <DETECTOR_NAME>HO</DETECTOR_NAME>");
00249         if(detid.subdet() == 4) sprintf(output, "      <DETECTOR_NAME>HF</DETECTOR_NAME>");
00250         m_outputFileXML << output << endl;
00251         sprintf(output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId() );
00252         m_outputFileXML << output << endl;
00253         sprintf(output, "    </CHANNEL>");
00254         m_outputFileXML << output << endl;
00255         sprintf(output, "    <DATA>");
00256         m_outputFileXML << output << endl;
00257         sprintf(output, "      <MEAN_TIME>%7f</MEAN_TIME>", time1);
00258         m_outputFileXML << output << endl;
00259         sprintf(output, "      <OFFSET_TIME> 0</OFFSET_TIME>");
00260         m_outputFileXML << output << endl;
00261         sprintf(output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
00262         m_outputFileXML << output << endl;
00263         sprintf(output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag+1);
00264         m_outputFileXML << output << endl;
00265         sprintf(output, "      <STATUS_WORD>  0</STATUS_WORD>");
00266         m_outputFileXML << output << endl;
00267         sprintf(output, "    </DATA>");
00268         m_outputFileXML << output << endl;
00269         sprintf(output, "  </DATA_SET>");
00270         m_outputFileXML << output << endl;
00271 
00272         }
00273       else if(m_fitflag==1){
00274         m_outFile<<detid<<"   "<<time2<<" "<<dtime2<<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>", time2);
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>", dtime2);
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==2){
00319         m_outFile<<detid<<"   "<<time3<<" "<<dtime3<<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>", time3);
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>", dtime3);
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       else if(m_fitflag==3){
00363         m_outFile<<detid<<"   "<<time4<<" "<<dtime4<<std::endl;
00364         sprintf(output, "  <DATA_SET>");
00365         m_outputFileXML << output << endl;
00366         sprintf(output, "    <VERSION>version:1</VERSION>");
00367         m_outputFileXML << output << endl;
00368         sprintf(output, "    <CHANNEL>");
00369         m_outputFileXML << output << endl;
00370         sprintf(output, "      <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
00371         m_outputFileXML << output << endl;
00372         sprintf(output, "      <ETA>%2i</ETA>", detid.ietaAbs() );
00373         m_outputFileXML << output << endl;
00374         sprintf(output, "      <PHI>%2i</PHI>", detid.iphi() );
00375         m_outputFileXML << output << endl;
00376         sprintf(output, "      <DEPTH>%2i</DEPTH>", detid.depth() );
00377         m_outputFileXML << output << endl;
00378         sprintf(output, "      <Z>%2i</Z>", detid.zside() );
00379         m_outputFileXML << output << endl;
00380         if(detid.subdet() == 1) sprintf(output, "      <DETECTOR_NAME>HB</DETECTOR_NAME>");
00381         if(detid.subdet() == 2) sprintf(output, "      <DETECTOR_NAME>HE</DETECTOR_NAME>");
00382         if(detid.subdet() == 3) sprintf(output, "      <DETECTOR_NAME>HO</DETECTOR_NAME>");
00383         if(detid.subdet() == 4) sprintf(output, "      <DETECTOR_NAME>HF</DETECTOR_NAME>");
00384         m_outputFileXML << output << endl;
00385         sprintf(output, "      <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId() );
00386         m_outputFileXML << output << endl;
00387         sprintf(output, "    </CHANNEL>");
00388         m_outputFileXML << output << endl;
00389         sprintf(output, "    <DATA>");
00390         m_outputFileXML << output << endl;
00391         sprintf(output, "      <MEAN_TIME>%7f</MEAN_TIME>", time4);
00392         m_outputFileXML << output << endl;
00393         sprintf(output, "      <OFFSET_TIME> 0</OFFSET_TIME>");
00394         m_outputFileXML << output << endl;
00395         sprintf(output, "      <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
00396         m_outputFileXML << output << endl;
00397         sprintf(output, "      <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag+1);
00398         m_outputFileXML << output << endl;
00399         sprintf(output, "      <STATUS_WORD>  0</STATUS_WORD>");
00400         m_outputFileXML << output << endl;
00401         sprintf(output, "    </DATA>");
00402         m_outputFileXML << output << endl;
00403         sprintf(output, "  </DATA_SET>");
00404         m_outputFileXML << output << endl;
00405         }
00406 
00407       else if(m_fitflag==4){
00408         m_outFile<<detid<<"   "<<time1<<" "<<dtime1<<"   "<<time2<<" "<<dtime2<<"   "<<time3<<" "<<dtime3<<"   "<<time4<<" "<<dtime4<<std::endl;
00409         }
00410     }
00411   }
00412 }
00413 
00414 //-----------------------------------------------------------------------------
00415 void CastorLedAnalysis::LedSampleAnalysis(){
00416   // it is called every m_nevtsample events (a sample) and the end of run
00417   char LedSampleNum[20];
00418 
00419   sprintf(LedSampleNum,"LedSample_%d",sample);
00420   m_file->cd();
00421   m_file->mkdir(LedSampleNum);
00422   m_file->cd(LedSampleNum);
00423 
00424 // Compute LED constants
00425   GetLedConst(castorHists.LEDTRENDS);
00426 }
00427 
00428 //-----------------------------------------------------------------------------
00429 void CastorLedAnalysis::LedTrendings(map<HcalDetId, map<int,LEDBUNCH> > &toolT)
00430 {
00431 
00432   for(_meol=toolT.begin(); _meol!=toolT.end(); _meol++){
00433     char name[1024];
00434     HcalDetId detid = _meol->first;
00435     sprintf(name,"LED timing trend, eta=%d phi=%d depth=%d",detid.ieta(),detid.iphi(),detid.depth());
00436     int bins = _meol->second[10+m_fitflag].second.first[0].size();
00437     float lo =0.5;
00438     float hi = (float)bins+0.5;
00439     _meol->second[10+m_fitflag].second.second.push_back(new TH1F(name,name,bins,lo,hi));
00440 
00441     std::vector<double>::iterator sample_it;
00442 // LED timing - put content and errors
00443     int j=0;
00444     for(sample_it=_meol->second[10+m_fitflag].second.first[0].begin();
00445         sample_it!=_meol->second[10+m_fitflag].second.first[0].end();sample_it++){
00446       _meol->second[10+m_fitflag].second.second[0]->SetBinContent(++j,*sample_it);
00447     }
00448     j=0;
00449     for(sample_it=_meol->second[10+m_fitflag].second.first[1].begin();
00450         sample_it!=_meol->second[10+m_fitflag].second.first[1].end();sample_it++){
00451       _meol->second[10+m_fitflag].second.second[0]->SetBinError(++j,*sample_it);
00452     }
00453     sprintf(name,"Sample (%d events)",m_nevtsample);
00454     _meol->second[10+m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
00455     _meol->second[10+m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
00456     _meol->second[10+m_fitflag].second.second[0]->Write();
00457   }
00458 }
00459 
00460 //-----------------------------------------------------------------------------
00461 void CastorLedAnalysis::LedDone() 
00462 {
00463 
00464 // First process the last sample (remaining events).
00465   if(evt%m_nevtsample!=0) LedSampleAnalysis();
00466 
00467 // Now do the end of run analysis: trending histos
00468   if(sample>1 && m_fitflag!=4){
00469     m_file->cd();
00470     m_file->cd("Castor");
00471     LedTrendings(castorHists.LEDTRENDS);
00472   }
00473 
00474   // Write other histograms.
00475   m_file->cd();
00476   m_file->cd("Castor");
00477   castorHists.ALLLEDS->Write();
00478   castorHists.LEDRMS->Write();
00479   castorHists.LEDMEAN->Write();
00480 
00481   // Write the histo file and close it
00482 //  m_file->Write();
00483   m_file->Close();
00484   cout << "Castor histograms written to " << m_outputFileROOT.c_str() << endl;
00485 }
00486 
00487 //-----------------------------------------------------------------------------
00488 void CastorLedAnalysis::processLedEvent(const CastorDigiCollection& castor,
00489                                         const CastorDbService& cond)
00490 {
00491   evt++;
00492   sample = (evt-1)/m_nevtsample +1;
00493   evt_curr = evt%m_nevtsample;
00494   if(evt_curr==0)evt_curr=m_nevtsample;
00495 
00496   // HF/Castor
00497   try{
00498     if(!castor.size()) throw (int)castor.size();
00499     for (CastorDigiCollection::const_iterator j=castor.begin(); j!=castor.end(); j++){
00500       const CastorDataFrame digi = (const CastorDataFrame)(*j);
00501       _meol = castorHists.LEDTRENDS.find(digi.id());
00502       if (_meol==castorHists.LEDTRENDS.end()){
00503         SetupLEDHists(2,digi.id(),castorHists.LEDTRENDS);
00504       }
00505       LedCastorHists(digi.id(),digi,castorHists.LEDTRENDS,cond);
00506     }
00507   } 
00508   catch (int i ) {
00509 //    m_logFile << "Event with " << i<<" Castor Digis passed." << std::endl;
00510   } 
00511 
00512   // Call the function every m_nevtsample events
00513   if(evt%m_nevtsample==0) LedSampleAnalysis();
00514 
00515 }
00516 //----------------------------------------------------------------------------
00517 void CastorLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int,LEDBUNCH> > &toolT) {
00518 
00519   string type = "HBHE";
00520   if(id==1) type = "HO";
00521   if(id==2) type = "HF";
00522 
00523   _meol = toolT.find(detid);
00524   if (_meol==toolT.end()){
00525 // if histos for this channel do not exist, create them
00526     map<int,LEDBUNCH> insert;
00527     char name[1024];
00528     for(int i=0; i<10; i++){
00529       sprintf(name,"%s Pulse height, eta=%d phi=%d depth=%d TS=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth(),i);
00530       insert[i].first =  new TH1F(name,name,200,0.,2000.);
00531     }
00532     sprintf(name,"%s LED Mean pulse, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00533     insert[10].first =  new TH1F(name,name,10,-0.5,9.5);
00534     sprintf(name,"%s LED Pulse, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00535     insert[11].first =  new TH1F(name,name,10,-0.5,9.5);
00536     sprintf(name,"%s Mean TS, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00537     insert[12].first =  new TH1F(name,name,200,0.,10.);
00538     sprintf(name,"%s Peak TS, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00539     insert[13].first =  new TH1F(name,name,200,0.,10.);
00540     sprintf(name,"%s Peak TS error, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00541     insert[14].first =  new TH1F(name,name,200,0.,0.05);
00542     sprintf(name,"%s Fit chi2, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00543     insert[15].first =  new TH1F(name,name,100,0.,50.);
00544     sprintf(name,"%s Integrated Signal, eta=%d phi=%d depth=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
00545     insert[16].first =  new TH1F(name,name,500,0.,5000.);
00546 
00547     toolT[detid] = insert;
00548   }
00549 }
00550 //-----------------------------------------------------------------------------
00551 void CastorLedAnalysis::LedCastorHists(const HcalDetId& detid, const CastorDataFrame& ledDigi, map<HcalDetId, map<int,LEDBUNCH> > &toolT, const CastorDbService& cond) {
00552 
00553   map<int,LEDBUNCH> _mei;
00554   _meol = toolT.find(detid);
00555   _mei = _meol->second;
00556   // Rest the histos if we're at the end of a 'bunch'
00557   if((evt-1)%m_nevtsample==0 && state[0]){
00558     for(int k=0; k<(int)state.size();k++) state[k]=false;
00559     for(int i=0; i<16; i++) _mei[i].first->Reset();
00560   }
00561 
00562   // now we have the signal in fC, let's get rid of that darn pedestal
00563   // Most of this is borrowed from CastorSimpleReconstructor, so thanks Jeremy/Phil
00564 
00565   int maxTS = -1;
00566   float max_fC = 0;
00567   float ta = 0;
00568   m_coder = cond.getCastorCoder(detid);
00569   m_ped = cond.getPedestal(detid);
00570   m_shape = cond.getCastorShape();
00571   //cout << "New Digi!!!!!!!!!!!!!!!!!!!!!!" << endl;
00572   for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++){
00573     int capid = ledDigi[TS].capid();
00574     // BE CAREFUL: this is assuming peds are stored in ADCs
00575     int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
00576     if (adc < 0){ adc = 0; }  // to prevent negative adcs after ped subtraction, which should really only happen
00577                               // if you're using the wrong peds.
00578     double fC = m_coder->charge(*m_shape,adc,capid);
00579     //ta = (fC - m_ped->getValue(capid));
00580     ta = fC;
00581     //cout << "DetID: " << detid << "  CapID: " << capid << "  ADC: " << adc << "  Ped: " << m_ped->getValue(capid) << "  fC: " << fC << endl;
00582     _mei[TS].first->Fill(ta);
00583     _mei[10].first->AddBinContent(TS+1,ta);  // This is average pulse, could probably do better (Profile?)
00584     if(m_fitflag>1){
00585       if(TS==m_startTS)_mei[11].first->Reset();
00586       _mei[11].first->SetBinContent(TS+1,ta);
00587     }
00588 
00589     // keep track of max TS and max amplitude (in fC)
00590     if (ta > max_fC){
00591       max_fC = ta;
00592       maxTS = TS;
00593     }
00594   }
00595 
00596   // Now we have a sample with pedestals subtracted and in units of fC
00597   // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
00598   // we now want to use Phil's timing correction.  This is not necessary
00599   // if we are performing a Landau fit (m_fitflag = 3)
00600 
00601   float sum=0.;
00602   for(int i=0; i<10; i++)sum=sum+_mei[11].first->GetBinContent(i+1);
00603   if(sum>100){
00604     if(m_fitflag==2 || m_fitflag==4){
00605       float timmean=_mei[11].first->GetMean();  // let's use Phil's way instead
00606       float timmeancorr=BinsizeCorr(timmean);
00607       _mei[12].first->Fill(timmeancorr);
00608     }
00609     _mei[16].first->Fill(_mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
00610     if(m_fitflag==3 || m_fitflag==4){
00611       _mei[11].first->Fit("landau","Q");
00612       TF1 *fit = _mei[11].first->GetFunction("landau");
00613       _mei[13].first->Fill(fit->GetParameter(1));
00614       _mei[14].first->Fill(fit->GetParError(1));
00615       _mei[15].first->Fill(fit->GetChisquare()/fit->GetNDF());
00616     }
00617   }
00618 
00619 
00620 
00621 }
00622 
00623 
00624 //-----------------------------------------------------------------------------
00625 float CastorLedAnalysis::BinsizeCorr(float time) {
00626 
00627 // this is the bin size correction to be applied for laser data (from Andy),
00628 // it comes from a pulse shape measured from TB04 data (from Jordan)
00629 // This should eventually be replaced with the more thorough treatment from Phil
00630 
00631   float corrtime=0.;
00632   static const float tstrue[32]={0.003, 0.03425, 0.06548, 0.09675, 0.128,
00633  0.15925, 0.1905, 0.22175, 0.253, 0.28425, 0.3155, 0.34675, 0.378, 0.40925,
00634  0.4405, 0.47175, 0.503, 0.53425, 0.5655, 0.59675, 0.628, 0.65925, 0.6905,
00635  0.72175, 0.753, 0.78425, 0.8155, 0.84675, 0.878, 0.90925, 0.9405, 0.97175};
00636   static const float tsreco[32]={-0.00422, 0.01815, 0.04409, 0.07346, 0.09799,
00637  0.12192, 0.15072, 0.18158, 0.21397, 0.24865, 0.28448, 0.31973, 0.35449,
00638  0.39208, 0.43282, 0.47244, 0.5105, 0.55008, 0.58827, 0.62828, 0.6717, 0.70966,
00639  0.74086, 0.77496, 0.80843, 0.83472, 0.86044, 0.8843, 0.90674, 0.92982,
00640  0.95072, 0.9726};
00641 
00642  int inttime=(int)time;
00643  float restime=time-inttime;
00644  for(int i=0; i<=32; i++) {
00645    float lolim=0.; float uplim=1.; float tsdown; float tsup;
00646    if(i>0){
00647      lolim=tsreco[i-1];
00648      tsdown=tstrue[i-1];
00649    }
00650    else tsdown=tstrue[31]-1.;
00651    if(i<32){
00652      uplim=tsreco[i];
00653      tsup=tstrue[i];
00654    }
00655    else tsup=tstrue[0]+1.;
00656    if(restime>=lolim && restime<uplim){
00657       corrtime=(tsdown*(uplim-restime)+tsup*(restime-lolim)) / (uplim-lolim);
00658     }
00659   }
00660   corrtime+=inttime;
00661 
00662  return corrtime;
00663 }
00664 //-----------------------------------------------------------------------------
00665
CastorLedAnalysis.cc
