CMSSW Reference Manual

00001 // Original Author:  Steven Won
00002 //         Created:  Fri May  2 15:34:43 CEST 2008
00003 // Written to replace the combination of HcalPedestalAnalyzer and HcalPedestalAnalysis 
00004 // This code runs 1000x faster and produces all outputs from a single run
00005 // (ADC, fC in .txt plus an .xml file)
00006 //
00007 #include <memory>
00008 #include "CalibCalorimetry/HcalStandardModules/interface/HcalPedestalsAnalysis.h"
00009 
00010 HcalPedestalsAnalysis::HcalPedestalsAnalysis(const edm::ParameterSet& ps)
00011 {
00012    hiSaveFlag = ps.getUntrackedParameter<bool>("hiSaveFlag", false);
00013    dumpXML = ps.getUntrackedParameter<bool>("dumpXML", true);
00014    verboseflag = ps.getUntrackedParameter<bool>("verbose", false);
00015    firstTS = ps.getUntrackedParameter<int>("firstTS", 0);
00016    lastTS = ps.getUntrackedParameter<int>("lastTS", 9);   
00017    firsttime = true;
00018 }
00019 
00020 
00021 HcalPedestalsAnalysis::~HcalPedestalsAnalysis()
00022 {
00023    HcalPedestals* rawPedsItem = new HcalPedestals();
00024    HcalPedestalWidths* rawWidthsItem = new HcalPedestalWidths();
00025 //   HcalPedestals* rawPedsItemfc = new HcalPedestals();
00026 //   HcalPedestalWidths* rawWidthsItemfc = new HcalPedestalWidths();
00027    HcalZSThresholds* ZSItem = new HcalZSThresholds();
00028 
00029    //Calculate pedestal constants
00030    std::cout << "Calculating Pedestal constants...\n";
00031    std::vector<NewPedBunch>::iterator bunch_it;
00032    for(bunch_it=Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00033    {
00034       if(bunch_it->usedflag){
00035 
00036       if(verboseflag) std::cout << "Analyzing channel phi= " << bunch_it->detid.iphi() 
00037         << " eta = " << bunch_it->detid.ieta() << " depth = " << bunch_it->detid.depth()
00038         << std::endl;
00039       //pedestal constant is the mean
00040       bunch_it->cap[0] /= bunch_it->num[0][0];
00041       bunch_it->cap[1] /= bunch_it->num[1][1];
00042       bunch_it->cap[2] /= bunch_it->num[2][2];
00043       bunch_it->cap[3] /= bunch_it->num[3][3];
00044 //      bunch_it->capfc[0] /= bunch_it->num[0][0];
00045 //      bunch_it->capfc[1] /= bunch_it->num[1][1];
00046 //      bunch_it->capfc[2] /= bunch_it->num[2][2];
00047 //      bunch_it->capfc[3] /= bunch_it->num[3][3];
00048       //widths are the covariance matrix--assumed symmetric
00049       bunch_it->sig[0][0] = (bunch_it->prod[0][0]/bunch_it->num[0][0])-(bunch_it->cap[0])*(bunch_it->cap[0]);
00050       bunch_it->sig[1][1] = (bunch_it->prod[1][1]/bunch_it->num[1][1])-(bunch_it->cap[1])*(bunch_it->cap[1]);
00051       bunch_it->sig[2][2] = (bunch_it->prod[2][2]/bunch_it->num[2][2])-(bunch_it->cap[2])*(bunch_it->cap[2]);
00052       bunch_it->sig[3][3] = (bunch_it->prod[3][3]/bunch_it->num[3][3])-(bunch_it->cap[3])*(bunch_it->cap[3]);
00053       bunch_it->sig[0][1] = (bunch_it->prod[0][1])/(bunch_it->num[0][1])-(bunch_it->cap[0]*bunch_it->cap[1]);
00054       bunch_it->sig[0][2] = (bunch_it->prod[0][2])/(bunch_it->num[0][2])-(bunch_it->cap[0]*bunch_it->cap[2]);
00055       bunch_it->sig[0][3] = (bunch_it->prod[0][3])/(bunch_it->num[0][3])-(bunch_it->cap[0]*bunch_it->cap[3]);
00056       bunch_it->sig[1][2] = (bunch_it->prod[1][2])/(bunch_it->num[1][2])-(bunch_it->cap[1]*bunch_it->cap[2]);
00057       bunch_it->sig[1][3] = (bunch_it->prod[1][3])/(bunch_it->num[1][3])-(bunch_it->cap[1]*bunch_it->cap[3]);
00058       bunch_it->sig[2][3] = (bunch_it->prod[2][3])/(bunch_it->num[2][3])-(bunch_it->cap[2]*bunch_it->cap[3]);
00059 
00060 //      bunch_it->sigfc[0][0] = (bunch_it->prodfc[0][0]/bunch_it->num[0][0])-(bunch_it->capfc[0])*(bunch_it->capfc[0]);
00061 //      bunch_it->sigfc[1][1] = (bunch_it->prodfc[1][1]/bunch_it->num[1][1])-(bunch_it->capfc[1])*(bunch_it->capfc[1]);
00062 //      bunch_it->sigfc[2][2] = (bunch_it->prodfc[2][2]/bunch_it->num[2][2])-(bunch_it->capfc[2])*(bunch_it->capfc[2]);
00063 //      bunch_it->sigfc[3][3] = (bunch_it->prodfc[3][3]/bunch_it->num[3][3])-(bunch_it->capfc[3])*(bunch_it->capfc[3]);
00064 //      bunch_it->sigfc[0][1] = (bunch_it->prodfc[0][1]/(bunch_it->num[0][1]))-(bunch_it->capfc[0]*bunch_it->capfc[1]);
00065 //      bunch_it->sigfc[0][2] = (bunch_it->prodfc[0][2]/(bunch_it->num[0][2]))-(bunch_it->capfc[0]*bunch_it->capfc[2]);
00066 //      bunch_it->sigfc[0][3] = (bunch_it->prodfc[0][3]/(bunch_it->num[0][3]))-(bunch_it->capfc[0]*bunch_it->capfc[3]);
00067 //      bunch_it->sigfc[1][2] = (bunch_it->prodfc[1][2]/(bunch_it->num[1][2]))-(bunch_it->capfc[1]*bunch_it->capfc[2]);
00068 //      bunch_it->sigfc[1][3] = (bunch_it->prodfc[1][3]/(bunch_it->num[1][3]))-(bunch_it->capfc[1]*bunch_it->capfc[3]);
00069 //      bunch_it->sigfc[2][3] = (bunch_it->prodfc[2][3]/(bunch_it->num[2][3]))-(bunch_it->capfc[2]*bunch_it->capfc[3]);
00070 
00071       if(bunch_it->detid.subdet() == 1){
00072          for(int i = 0; i != 3; i++){
00073             HBMeans->Fill(bunch_it->cap[i]);
00074             HBWidths->Fill(bunch_it->sig[i][i]);
00075          }
00076       }
00077       if(bunch_it->detid.subdet() == 2){
00078          for(int i = 0; i != 3; i++){
00079             HEMeans->Fill(bunch_it->cap[i]);
00080             HEWidths->Fill(bunch_it->sig[i][i]);
00081          }
00082       }
00083       if(bunch_it->detid.subdet() == 3){
00084          for(int i = 0; i != 3; i++){
00085             HOMeans->Fill(bunch_it->cap[i]);
00086             HOWidths->Fill(bunch_it->sig[i][i]);
00087          }
00088       }
00089       if(bunch_it->detid.subdet() == 4){
00090          for(int i = 0; i != 3; i++){
00091             HFMeans->Fill(bunch_it->cap[i]);
00092             HFWidths->Fill(bunch_it->sig[i][i]);
00093          }
00094       }
00095 
00096       const HcalPedestal item(bunch_it->detid, bunch_it->cap[0], bunch_it->cap[1], bunch_it->cap[2], bunch_it->cap[3]);
00097       rawPedsItem->addValues(item);
00098       HcalPedestalWidth widthsp(bunch_it->detid);
00099       widthsp.setSigma(0,0,bunch_it->sig[0][0]);
00100       widthsp.setSigma(0,1,bunch_it->sig[0][1]);
00101       widthsp.setSigma(0,2,bunch_it->sig[0][2]);
00102       widthsp.setSigma(0,3,bunch_it->sig[0][3]);
00103       widthsp.setSigma(1,1,bunch_it->sig[1][1]);
00104       widthsp.setSigma(1,2,bunch_it->sig[1][2]);
00105       widthsp.setSigma(1,3,bunch_it->sig[1][3]);
00106       widthsp.setSigma(2,2,bunch_it->sig[2][2]);
00107       widthsp.setSigma(2,3,bunch_it->sig[2][3]);
00108       widthsp.setSigma(3,3,bunch_it->sig[3][3]);
00109       rawWidthsItem->addValues(widthsp);
00110 
00111 /*      const HcalPedestal itemfc(bunch_it->detid, bunch_it->capfc[0], bunch_it->capfc[1],
00112                                    bunch_it->capfc[2], bunch_it->capfc[3]);
00113       rawPedsItemfc->addValues(itemfc);
00114       HcalPedestalWidth widthspfc(bunch_it->detid);
00115       widthspfc.setSigma(0,0,bunch_it->sigfc[0][0]);
00116       widthspfc.setSigma(0,1,bunch_it->sigfc[0][1]);
00117       widthspfc.setSigma(0,2,bunch_it->sigfc[0][2]);
00118       widthspfc.setSigma(0,3,bunch_it->sigfc[0][3]);
00119       widthspfc.setSigma(1,1,bunch_it->sigfc[1][1]);
00120       widthspfc.setSigma(1,2,bunch_it->sigfc[1][2]);
00121       widthspfc.setSigma(1,3,bunch_it->sigfc[1][3]);
00122       widthspfc.setSigma(2,2,bunch_it->sigfc[2][2]);
00123       widthspfc.setSigma(2,3,bunch_it->sigfc[2][3]);
00124       widthspfc.setSigma(3,3,bunch_it->sigfc[3][3]);
00125       rawWidthsItemfc->addValues(widthspfc);*/
00126 
00127       int threshold = (int)floor(
00128                                 (3 * sqrt(bunch_it->sig[0][0]) + bunch_it->cap[0] +
00129                                  3 * sqrt(bunch_it->sig[1][1]) + bunch_it->cap[1] +
00130                                  3 * sqrt(bunch_it->sig[2][2]) + bunch_it->cap[2] +
00131                                  3 * sqrt(bunch_it->sig[3][3]) + bunch_it->cap[3]) / 4);
00132       const HcalZSThreshold zerosup(bunch_it->detid, threshold);
00133       ZSItem->addValues(zerosup);
00134       }
00135    }
00136 
00137     // dump the resulting list of pedestals into a file
00138     std::ofstream outStream1(pedsADCfilename.c_str());
00139     HcalDbASCIIIO::dumpObject (outStream1, (*rawPedsItem) );
00140     std::ofstream outStream2(widthsADCfilename.c_str());
00141     HcalDbASCIIIO::dumpObject (outStream2, (*rawWidthsItem) );
00142 
00143 /*    std::ofstream outStream3(pedsfCfilename.c_str());
00144     HcalDbASCIIIO::dumpObject (outStream3, (*rawPedsItemfc) );
00145     std::ofstream outStream4(widthsfCfilename.c_str());
00146     HcalDbASCIIIO::dumpObject (outStream4, (*rawWidthsItemfc) );
00147 */
00148 
00149     std::ofstream outStream8(ZSfilename.c_str());
00150     HcalDbASCIIIO::dumpObject (outStream8, (*ZSItem) );
00151 
00152     if(dumpXML){
00153        std::ofstream outStream5(XMLfilename.c_str());
00154        HcalDbXml::dumpObject (outStream5, runnum, 0, 2147483647, XMLtag, 1, (*rawPedsItem), (*rawWidthsItem)); 
00155     }
00156 
00157     if(hiSaveFlag){
00158        theFile->Write();
00159     }else{
00160        theFile->cd();
00161        theFile->cd("HB");
00162        HBMeans->Write();
00163        HBWidths->Write();
00164        theFile->cd();
00165        theFile->cd("HF");
00166        HFMeans->Write();
00167        HFWidths->Write();
00168        theFile->cd();
00169        theFile->cd("HE");
00170        HEMeans->Write();
00171        HEWidths->Write();
00172        theFile->cd();
00173        theFile->cd("HO");
00174        HOMeans->Write();
00175        HOWidths->Write();
00176     }
00177 
00178     std::cout << "Writing ROOT file... ";
00179     theFile->Close();
00180     std::cout << "ROOT file closed.\n";
00181 }
00182 
00183 // ------------ method called to for each event  ------------
00184 void
00185 HcalPedestalsAnalysis::analyze(const edm::Event& e, const edm::EventSetup& iSetup)
00186 {
00187    using namespace edm;
00188    using namespace std;
00189 
00190    edm::Handle<HBHEDigiCollection> hbhe;              e.getByType(hbhe);
00191    edm::Handle<HODigiCollection> ho;                  e.getByType(ho);
00192    edm::Handle<HFDigiCollection> hf;                  e.getByType(hf);
00193    edm::ESHandle<HcalDbService> conditions;
00194    iSetup.get<HcalDbRecord>().get(conditions);
00195 
00196 //   const HcalQIEShape* shape = conditions->getHcalShape();
00197 
00198    if(firsttime)
00199    {
00200       runnum = e.id().run();
00201       std::string runnum_string;
00202       std::stringstream tempstringout;
00203       tempstringout << runnum;
00204       runnum_string = tempstringout.str();
00205       ROOTfilename = runnum_string + "-peds_ADC.root";
00206       pedsADCfilename = runnum_string + "-peds_ADC.txt";
00207 //      pedsfCfilename = runnum_string + "-peds_fC.txt";
00208       widthsADCfilename = runnum_string + "-widths_ADC.txt";
00209 //      widthsfCfilename = runnum_string + "-widths_fC.txt";
00210       XMLfilename = runnum_string + "-peds_ADC_complete.xml"; 
00211       XMLtag = "Hcal_pedestals_" + runnum_string;
00212       ZSfilename = runnum_string + "-ZSThresholds.txt";
00213 
00214       theFile = new TFile(ROOTfilename.c_str(), "RECREATE");
00215       theFile->cd();
00216       // Create sub-directories
00217       theFile->mkdir("HB"); 
00218       theFile->mkdir("HE");
00219       theFile->mkdir("HF");
00220       theFile->mkdir("HO");
00221       theFile->cd();
00222 
00223       HBMeans = new TH1F("All Ped Means HB","All Ped Means HB", 100, 0, 9);
00224       HBWidths = new TH1F("All Ped Widths HB","All Ped Widths HB", 100, 0, 3);
00225       HEMeans = new TH1F("All Ped Means HE","All Ped Means HE", 100, 0, 9);
00226       HEWidths = new TH1F("All Ped Widths HE","All Ped Widths HE", 100, 0, 3);
00227       HFMeans = new TH1F("All Ped Means HF","All Ped Means HF", 100, 0, 9);
00228       HFWidths = new TH1F("All Ped Widths HF","All Ped Widths HF", 100, 0, 3);
00229       HOMeans = new TH1F("All Ped Means HO","All Ped Means HO", 100, 0, 9);
00230       HOWidths = new TH1F("All Ped Widths HO","All Ped Widths HO", 100, 0, 3);
00231 
00232       edm::ESHandle<HcalElectronicsMap> refEMap;
00233       iSetup.get<HcalElectronicsMapRcd>().get(refEMap);
00234       const HcalElectronicsMap* myRefEMap = refEMap.product();
00235       std::vector<HcalGenericDetId> listEMap = myRefEMap->allPrecisionId();
00236       for (std::vector<HcalGenericDetId>::const_iterator it = listEMap.begin(); it != listEMap.end(); it++)
00237       {     
00238          HcalGenericDetId mygenid(it->rawId());
00239          if(mygenid.isHcalDetId())
00240          {
00241             NewPedBunch a;
00242             HcalDetId chanid(mygenid.rawId());
00243             a.detid = chanid;
00244             a.usedflag = false;
00245             string type;
00246             if(chanid.subdet() == 1) type = "HB";
00247             if(chanid.subdet() == 2) type = "HE";
00248             if(chanid.subdet() == 3) type = "HO";
00249             if(chanid.subdet() == 4) type = "HF";
00250             for(int i = 0; i != 4; i++)
00251             {
00252                a.cap[i] = 0;
00253 //               a.capfc[i] = 0;
00254                for(int j = 0; j != 4; j++)
00255                {
00256                   a.sig[i][j] = 0;
00257 //                  a.sigfc[i][j] = 0;
00258                   a.prod[i][j] = 0;
00259 //                  a.prodfc[i][j] = 0;
00260                   a.num[i][j] = 0;
00261                }
00262             }
00263             Bunches.push_back(a);
00264          }
00265       }
00266       firsttime = false;
00267    }
00268 
00269    std::vector<NewPedBunch>::iterator bunch_it;
00270 
00271    for(HBHEDigiCollection::const_iterator j = hbhe->begin(); j != hbhe->end(); j++)
00272    {
00273       const HBHEDataFrame digi = (const HBHEDataFrame)(*j);
00274       for(bunch_it = Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00275          if(bunch_it->detid.rawId() == digi.id().rawId()) break;
00276       bunch_it->usedflag = true;
00277       for(int ts = firstTS; ts != lastTS+1; ts++)
00278       {
00279 //         const HcalQIECoder* coder = conditions->getHcalCoder(digi.id().rawId());
00280          bunch_it->num[digi.sample(ts).capid()][digi.sample(ts).capid()] += 1;
00281          bunch_it->cap[digi.sample(ts).capid()] += digi.sample(ts).adc();
00282 //         double charge1 = coder->charge(*shape, digi.sample(ts).adc(), digi.sample(ts).capid());
00283 //         bunch_it->capfc[digi.sample(ts).capid()] += charge1;
00284          bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts).capid()] += (digi.sample(ts).adc() * digi.sample(ts).adc());
00285 //         bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts).capid()] += charge1 * charge1;
00286          if((ts+1 < digi.size()) && (ts+1 < lastTS)){
00287             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += digi.sample(ts).adc()*digi.sample(ts+1).adc();
00288 //            double charge2 = coder->charge(*shape, digi.sample(ts+1).adc(), digi.sample(ts+1).capid());
00289 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += charge1*charge2;
00290             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += 1;
00291          }
00292          if((ts+2 < digi.size()) && (ts+2 < lastTS)){
00293             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += digi.sample(ts).adc()*digi.sample(ts+2).adc();
00294 //            double charge2 = coder->charge(*shape, digi.sample(ts+2).adc(), digi.sample(ts+2).capid());
00295 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += charge1*charge2;
00296             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += 1;
00297          }
00298          if((ts+3 < digi.size()) && (ts+3 < lastTS)){
00299             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += digi.sample(ts).adc()*digi.sample(ts+3).adc();
00300 //            double charge2 = coder->charge(*shape, digi.sample(ts+3).adc(), digi.sample(ts+3).capid());
00301 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += charge1*charge2;
00302             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += 1;
00303          }
00304       }
00305    }
00306 
00307    for(HODigiCollection::const_iterator j = ho->begin(); j != ho->end(); j++)
00308    {
00309       const HODataFrame digi = (const HODataFrame)(*j);
00310       for(bunch_it = Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00311          if(bunch_it->detid.rawId() == digi.id().rawId()) break;
00312       bunch_it->usedflag = true;
00313       for(int ts = firstTS; ts <= lastTS; ts++)
00314       {
00315 //         const HcalQIECoder* coder = conditions->getHcalCoder(digi.id().rawId());
00316          bunch_it->num[digi.sample(ts).capid()][digi.sample(ts).capid()] += 1;
00317          bunch_it->cap[digi.sample(ts).capid()] += digi.sample(ts).adc();
00318 //         double charge1 = coder->charge(*shape, digi.sample(ts).adc(), digi.sample(ts).capid());
00319 //         bunch_it->capfc[digi.sample(ts).capid()] += charge1;
00320          bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts).capid()] += (digi.sample(ts).adc() * digi.sample(ts).adc());
00321 //         bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts).capid()] += charge1 * charge1;
00322          if((ts+1 < digi.size()) && (ts+1 < lastTS)){
00323             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += digi.sample(ts).adc()*digi.sample(ts+1).adc();
00324 //            double charge2 = coder->charge(*shape, digi.sample(ts+1).adc(), digi.sample(ts+1).capid());
00325 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += charge1*charge2;
00326             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += 1;
00327          }
00328          if((ts+2 < digi.size()) && (ts+2 < lastTS)){
00329             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += digi.sample(ts).adc()*digi.sample(ts+2).adc();
00330 //            double charge2 = coder->charge(*shape, digi.sample(ts+2).adc(), digi.sample(ts+2).capid());
00331 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += charge1*charge2;
00332             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += 1;
00333          }
00334          if((ts+3 < digi.size()) && (ts+3 < lastTS)){
00335             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += digi.sample(ts).adc()*digi.sample(ts+3).adc();
00336 //            double charge2 = coder->charge(*shape, digi.sample(ts+3).adc(), digi.sample(ts+3).capid());
00337 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += charge1*charge2;
00338             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += 1;
00339          }
00340       }
00341    }
00342 
00343    for(HFDigiCollection::const_iterator j = hf->begin(); j != hf->end(); j++)
00344    {
00345       const HFDataFrame digi = (const HFDataFrame)(*j);
00346       for(bunch_it = Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00347          if(bunch_it->detid.rawId() == digi.id().rawId()) break;
00348       bunch_it->usedflag = true;
00349       for(int ts = firstTS; ts <= lastTS; ts++)
00350       {
00351 //         const HcalQIECoder* coder = conditions->getHcalCoder(digi.id().rawId());
00352          bunch_it->num[digi.sample(ts).capid()][digi.sample(ts).capid()] += 1;
00353          bunch_it->cap[digi.sample(ts).capid()] += digi.sample(ts).adc();
00354 //         double charge1 = coder->charge(*shape, digi.sample(ts).adc(), digi.sample(ts).capid());
00355 //         bunch_it->capfc[digi.sample(ts).capid()] += charge1;
00356          bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts).capid()] += (digi.sample(ts).adc() * digi.sample(ts).adc());
00357 //         bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts).capid()] += charge1 * charge1;
00358          if((ts+1 < digi.size()) && (ts+1 < lastTS)){
00359             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += digi.sample(ts).adc()*digi.sample(ts+1).adc();
00360 //            double charge2 = coder->charge(*shape, digi.sample(ts+1).adc(), digi.sample(ts+1).capid());
00361 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += charge1*charge2;
00362             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += 1;
00363          }
00364          if((ts+2 < digi.size()) && (ts+2 < lastTS)){
00365             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += digi.sample(ts).adc()*digi.sample(ts+2).adc();
00366 //            double charge2 = coder->charge(*shape, digi.sample(ts+2).adc(), digi.sample(ts+2).capid());
00367 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += charge1*charge2;
00368             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += 1;
00369          }
00370          if((ts+3 < digi.size()) && (ts+3 < lastTS)){
00371             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += digi.sample(ts).adc()*digi.sample(ts+3).adc();
00372 //            double charge2 = coder->charge(*shape, digi.sample(ts+3).adc(), digi.sample(ts+3).capid());
00373 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += charge1*charge2;
00374             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += 1;
00375          }
00376       }
00377 
00378 /* Once I figure out how to unpack Calib digis they go here
00379       const HFDataFrame digi = (const HFDataFrame)(*j);
00380       for(bunch_it = Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00381          if(bunch_it->detid.rawId() == digi.id().rawId()) break;
00382       bunch_it->usedflag = true;
00383       for(int ts = firstTS; ts <= lastTS; ts++)
00384       {
00385 //         const HcalQIECoder* coder = conditions->getHcalCoder(digi.id().rawId());
00386          bunch_it->num[digi.sample(ts).capid()][digi.sample(ts).capid()] += 1;
00387          bunch_it->cap[digi.sample(ts).capid()] += digi.sample(ts).adc();
00388 //         double charge1 = coder->charge(*shape, digi.sample(ts).adc(), digi.sample(ts).capid());
00389 //         bunch_it->capfc[digi.sample(ts).capid()] += charge1;
00390          bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts).capid()] += (digi.sample(ts).adc() * digi.sample(ts).adc());
00391 //         bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts).capid()] += charge1 * charge1;
00392          if((ts+1 < digi.size()) && (ts+1 < lastTS)){
00393             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += digi.sample(ts).adc()*digi.sample(ts+1).adc();
00394 //            double charge2 = coder->charge(*shape, digi.sample(ts+1).adc(), digi.sample(ts+1).capid());
00395 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += charge1*charge2;
00396             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += 1;
00397          }
00398          if((ts+2 < digi.size()) && (ts+2 < lastTS)){
00399             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += digi.sample(ts).adc()*digi.sample(ts+2).adc();
00400 //            double charge2 = coder->charge(*shape, digi.sample(ts+2).adc(), digi.sample(ts+2).capid());
00401 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += charge1*charge2;
00402             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += 1;
00403          }
00404          if((ts+3 < digi.size()) && (ts+3 < lastTS)){
00405             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += digi.sample(ts).adc()*digi.sample(ts+3).adc();
00406 //            double charge2 = coder->charge(*shape, digi.sample(ts+3).adc(), digi.sample(ts+3).capid());
00407 //            bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += charge1*charge2;
00408             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += 1;
00409          }
00410 */
00411 
00412    }
00413 //this is the last brace
00414 }
00415 
00416 //define this as a plug-in
00417 DEFINE_FWK_MODULE(HcalPedestalsAnalysis);
HcalPedestalsAnalysis.cc
