CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/CalibCalorimetry/HcalStandardModules/src/HcalPedestalsAnalysis.cc

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