Reference Manual

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