Reference Manual

Go to the documentation of this file.
00001 // Original Author:  Steven Won
00002 // Original Author:  Alan Campbell for castor
00003 //         Created:  Fri May  2 15:34:43 CEST 2008
00004 // Written to replace the combination of HcalPedestalAnalyzer and HcalPedestalAnalysis 
00005 // This code runs 1000x faster and produces all outputs from a single run
00006 // (ADC, fC in .txt plus an .xml file)
00007 //
00008 #include <memory>
00009 #include "CalibCalorimetry/CastorCalib/interface/CastorPedestalsAnalysis.h"
00010 
00011 CastorPedestalsAnalysis::CastorPedestalsAnalysis(const edm::ParameterSet& ps) :
00012    castorDigiCollectionTag(ps.getParameter<edm::InputTag>("castorDigiCollectionTag"))
00013 {
00014    hiSaveFlag = ps.getUntrackedParameter<bool>("hiSaveFlag", false);
00015    dumpXML = ps.getUntrackedParameter<bool>("dumpXML", false);
00016    verboseflag = ps.getUntrackedParameter<bool>("verbose", false);
00017    firstTS = ps.getUntrackedParameter<int>("firstTS", 0);
00018    lastTS = ps.getUntrackedParameter<int>("lastTS", 9);   
00019    firsttime = true;
00020 }
00021 
00022 
00023 CastorPedestalsAnalysis::~CastorPedestalsAnalysis()
00024 {
00025    CastorPedestals* rawPedsItem = new CastorPedestals(true);
00026    CastorPedestalWidths* rawWidthsItem = new CastorPedestalWidths(true);
00027    CastorPedestals* rawPedsItemfc = new CastorPedestals(false);
00028    CastorPedestalWidths* rawWidthsItemfc = new CastorPedestalWidths(false);
00029 
00030    //Calculate pedestal constants
00031    std::cout << "Calculating Pedestal constants...\n";
00032    std::vector<NewPedBunch>::iterator bunch_it;
00033    for(bunch_it=Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00034    {
00035       if(bunch_it->usedflag){
00036 
00037       if(verboseflag) std::cout << "Analyzing channel sector= " << bunch_it->detid.sector() 
00038         << " module = " << bunch_it->detid.module() 
00039         << std::endl;
00040       //pedestal constant is the mean
00041       bunch_it->cap[0] /= bunch_it->num[0][0];
00042       bunch_it->cap[1] /= bunch_it->num[1][1];
00043       bunch_it->cap[2] /= bunch_it->num[2][2];
00044       bunch_it->cap[3] /= bunch_it->num[3][3];
00045       bunch_it->capfc[0] /= bunch_it->num[0][0];
00046       bunch_it->capfc[1] /= bunch_it->num[1][1];
00047       bunch_it->capfc[2] /= bunch_it->num[2][2];
00048       bunch_it->capfc[3] /= bunch_it->num[3][3];
00049       //widths are the covariance matrix--assumed symmetric
00050       bunch_it->sig[0][0] = (bunch_it->prod[0][0]/bunch_it->num[0][0])-(bunch_it->cap[0]*bunch_it->cap[0]);
00051       bunch_it->sig[0][1] = (bunch_it->prod[0][1]/bunch_it->num[0][1])-(bunch_it->cap[0]*bunch_it->cap[1]);
00052       bunch_it->sig[0][2] = (bunch_it->prod[0][2]/bunch_it->num[0][2])-(bunch_it->cap[0]*bunch_it->cap[2]);
00053       bunch_it->sig[0][3] = (bunch_it->prod[0][3]/bunch_it->num[0][3])-(bunch_it->cap[0]*bunch_it->cap[3]);
00054       bunch_it->sig[1][0] = (bunch_it->prod[1][0]/bunch_it->num[1][0])-(bunch_it->cap[1]*bunch_it->cap[0]);
00055       bunch_it->sig[1][1] = (bunch_it->prod[1][1]/bunch_it->num[1][1])-(bunch_it->cap[1]*bunch_it->cap[1]);
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][0] = (bunch_it->prod[2][0]/bunch_it->num[2][0])-(bunch_it->cap[2]*bunch_it->cap[0]);
00059       bunch_it->sig[2][1] = (bunch_it->prod[2][1]/bunch_it->num[2][1])-(bunch_it->cap[2]*bunch_it->cap[1]);
00060       bunch_it->sig[2][2] = (bunch_it->prod[2][2]/bunch_it->num[2][2])-(bunch_it->cap[2]*bunch_it->cap[2]);
00061       bunch_it->sig[2][3] = (bunch_it->prod[2][3]/bunch_it->num[2][3])-(bunch_it->cap[2]*bunch_it->cap[3]);
00062       bunch_it->sig[3][0] = (bunch_it->prod[3][0]/bunch_it->num[3][0])-(bunch_it->cap[3]*bunch_it->cap[0]);
00063       bunch_it->sig[3][1] = (bunch_it->prod[3][1]/bunch_it->num[3][1])-(bunch_it->cap[3]*bunch_it->cap[1]);
00064       bunch_it->sig[3][2] = (bunch_it->prod[3][2]/bunch_it->num[3][2])-(bunch_it->cap[3]*bunch_it->cap[2]);
00065       bunch_it->sig[3][3] = (bunch_it->prod[3][3]/bunch_it->num[3][3])-(bunch_it->cap[3]*bunch_it->cap[3]);
00066 
00067       bunch_it->sigfc[0][0] = (bunch_it->prodfc[0][0]/bunch_it->num[0][0])-(bunch_it->capfc[0]*bunch_it->capfc[0]);
00068       bunch_it->sigfc[0][1] = (bunch_it->prodfc[0][1]/bunch_it->num[0][1])-(bunch_it->capfc[0]*bunch_it->capfc[1]);
00069       bunch_it->sigfc[0][2] = (bunch_it->prodfc[0][2]/bunch_it->num[0][2])-(bunch_it->capfc[0]*bunch_it->capfc[2]);
00070       bunch_it->sigfc[0][3] = (bunch_it->prodfc[0][3]/bunch_it->num[0][3])-(bunch_it->capfc[0]*bunch_it->capfc[3]);
00071       bunch_it->sigfc[1][0] = (bunch_it->prodfc[1][0]/bunch_it->num[1][0])-(bunch_it->capfc[1]*bunch_it->capfc[0]);
00072       bunch_it->sigfc[1][1] = (bunch_it->prodfc[1][1]/bunch_it->num[1][1])-(bunch_it->capfc[1]*bunch_it->capfc[1]);
00073       bunch_it->sigfc[1][2] = (bunch_it->prodfc[1][2]/bunch_it->num[1][2])-(bunch_it->capfc[1]*bunch_it->capfc[2]);
00074       bunch_it->sigfc[1][3] = (bunch_it->prodfc[1][3]/bunch_it->num[1][3])-(bunch_it->capfc[1]*bunch_it->capfc[3]);
00075       bunch_it->sigfc[2][0] = (bunch_it->prodfc[2][0]/bunch_it->num[2][0])-(bunch_it->capfc[2]*bunch_it->capfc[0]);
00076       bunch_it->sigfc[2][1] = (bunch_it->prodfc[2][1]/bunch_it->num[2][1])-(bunch_it->capfc[2]*bunch_it->capfc[1]);
00077       bunch_it->sigfc[2][2] = (bunch_it->prodfc[2][2]/bunch_it->num[2][2])-(bunch_it->capfc[2]*bunch_it->capfc[2]);
00078       bunch_it->sigfc[2][3] = (bunch_it->prodfc[2][3]/bunch_it->num[2][3])-(bunch_it->capfc[2]*bunch_it->capfc[3]);
00079       bunch_it->sigfc[3][0] = (bunch_it->prodfc[3][0]/bunch_it->num[3][0])-(bunch_it->capfc[3]*bunch_it->capfc[0]);
00080       bunch_it->sigfc[3][1] = (bunch_it->prodfc[3][1]/bunch_it->num[3][1])-(bunch_it->capfc[3]*bunch_it->capfc[1]);
00081       bunch_it->sigfc[3][2] = (bunch_it->prodfc[3][2]/bunch_it->num[3][2])-(bunch_it->capfc[3]*bunch_it->capfc[2]);
00082       bunch_it->sigfc[3][3] = (bunch_it->prodfc[3][3]/bunch_it->num[3][3])-(bunch_it->capfc[3]*bunch_it->capfc[3]);
00083 
00084         for(int i = 0; i != 3; i++){
00085             CASTORMeans->Fill(bunch_it->cap[i]);
00086             CASTORWidths->Fill(bunch_it->sig[i][i]);
00087          }
00088 
00089       //if(bunch_it->detid.subdet() == 1){
00090  
00091 
00092 
00093       int fillphi = bunch_it->detid.sector();
00094       //if (bunch_it->detid.depth()==4) fillphi++;
00095 
00096   //    dephist[bunch_it->detid.module()-1]->Fill(bunch_it->detid.ieta(),fillphi,
00097    //             (bunch_it->cap[0]+bunch_it->cap[1]+bunch_it->cap[2]+bunch_it->cap[3])/4);
00098       dephist->Fill( bunch_it->detid.module(),fillphi,
00099                 (bunch_it->cap[0]+bunch_it->cap[1]+bunch_it->cap[2]+bunch_it->cap[3])/4);
00100 
00101       const CastorPedestal item(bunch_it->detid, bunch_it->cap[0], bunch_it->cap[1], bunch_it->cap[2], bunch_it->cap[3],
00102                               bunch_it->sig[0][0], bunch_it->sig[1][1], bunch_it->sig[2][2], bunch_it->sig[3][3]);
00103       rawPedsItem->addValues(item);
00104       CastorPedestalWidth widthsp(bunch_it->detid);
00105       widthsp.setSigma(0,0,bunch_it->sig[0][0]);
00106       widthsp.setSigma(0,1,bunch_it->sig[0][1]);
00107       widthsp.setSigma(0,2,bunch_it->sig[0][2]);
00108       widthsp.setSigma(0,3,bunch_it->sig[0][3]);
00109       widthsp.setSigma(1,0,bunch_it->sig[1][0]);
00110       widthsp.setSigma(1,1,bunch_it->sig[1][1]);
00111       widthsp.setSigma(1,2,bunch_it->sig[1][2]);
00112       widthsp.setSigma(1,3,bunch_it->sig[1][3]);
00113       widthsp.setSigma(2,0,bunch_it->sig[2][0]);
00114       widthsp.setSigma(2,1,bunch_it->sig[2][1]);
00115       widthsp.setSigma(2,2,bunch_it->sig[2][2]);
00116       widthsp.setSigma(2,3,bunch_it->sig[2][3]);
00117       widthsp.setSigma(3,0,bunch_it->sig[3][0]);
00118       widthsp.setSigma(3,1,bunch_it->sig[3][1]);
00119       widthsp.setSigma(3,2,bunch_it->sig[3][2]);
00120       widthsp.setSigma(3,3,bunch_it->sig[3][3]);
00121       rawWidthsItem->addValues(widthsp);
00122 
00123       const CastorPedestal itemfc(bunch_it->detid, bunch_it->capfc[0], bunch_it->capfc[1], bunch_it->capfc[2], bunch_it->capfc[3],
00124                               bunch_it->sigfc[0][0], bunch_it->sigfc[1][1], bunch_it->sigfc[2][2], bunch_it->sigfc[3][3]);
00125       rawPedsItemfc->addValues(itemfc);
00126       CastorPedestalWidth widthspfc(bunch_it->detid);
00127       widthspfc.setSigma(0,0,bunch_it->sigfc[0][0]);
00128       widthspfc.setSigma(0,1,bunch_it->sigfc[0][1]);
00129       widthspfc.setSigma(0,2,bunch_it->sigfc[0][2]);
00130       widthspfc.setSigma(0,3,bunch_it->sigfc[0][3]);
00131       widthspfc.setSigma(1,0,bunch_it->sigfc[1][0]);      
00132       widthspfc.setSigma(1,1,bunch_it->sigfc[1][1]);
00133       widthspfc.setSigma(1,2,bunch_it->sigfc[1][2]);
00134       widthspfc.setSigma(1,3,bunch_it->sigfc[1][3]);
00135       widthspfc.setSigma(2,0,bunch_it->sigfc[2][0]);
00136       widthspfc.setSigma(2,1,bunch_it->sigfc[2][1]);
00137       widthspfc.setSigma(2,2,bunch_it->sigfc[2][2]);
00138       widthspfc.setSigma(2,3,bunch_it->sigfc[2][3]);
00139       widthspfc.setSigma(3,0,bunch_it->sigfc[3][0]);
00140       widthspfc.setSigma(3,1,bunch_it->sigfc[3][1]);
00141       widthspfc.setSigma(3,2,bunch_it->sigfc[3][2]);
00142       widthspfc.setSigma(3,3,bunch_it->sigfc[3][3]);
00143       rawWidthsItemfc->addValues(widthspfc);
00144 
00145       }
00146    }
00147 
00148     // dump the resulting list of pedestals into a file
00149     std::ofstream outStream1(pedsADCfilename.c_str());
00150     CastorDbASCIIIO::dumpObject (outStream1, (*rawPedsItem) );
00151     std::ofstream outStream2(widthsADCfilename.c_str());
00152     CastorDbASCIIIO::dumpObject (outStream2, (*rawWidthsItem) );
00153 
00154     std::ofstream outStream3(pedsfCfilename.c_str());
00155     CastorDbASCIIIO::dumpObject (outStream3, (*rawPedsItemfc) );
00156     std::ofstream outStream4(widthsfCfilename.c_str());
00157     CastorDbASCIIIO::dumpObject (outStream4, (*rawWidthsItemfc) );
00158 
00159     if(dumpXML){
00160        std::ofstream outStream5(XMLfilename.c_str());
00161      //  CastorCondXML::dumpObject (outStream5, runnum, runnum, runnum, XMLtag, 1, (*rawPedsItem), (*rawWidthsItem)); 
00162     }
00163 
00164     if(hiSaveFlag){
00165        theFile->Write();
00166     }else{
00167        theFile->cd();
00168        theFile->cd("CASTOR");
00169        CASTORMeans->Write();
00170        CASTORWidths->Write();
00171  
00172     }
00173     theFile->cd();
00174         dephist->Write();
00175         dephist->SetDrawOption("colz");
00176         dephist->GetXaxis()->SetTitle("module");
00177         dephist->GetYaxis()->SetTitle("sector");
00178     
00179     //for (int n=0; n!= 4; n++) 
00180     //{
00181          //dephist[n]->Write();
00182          //dephist[n]->SetDrawOption("colz");
00183          //dephist[n]->GetXaxis()->SetTitle("i#eta");
00184          //dephist[n]->GetYaxis()->SetTitle("i#phi");
00185     //}
00186 
00187     std::stringstream tempstringout;
00188     tempstringout << runnum;
00189     std::string name1 = tempstringout.str() + "_pedplots_1d.png";
00190     std::string name2 = tempstringout.str() + "_pedplots_2d.png";
00191 
00192     TStyle *theStyle = new TStyle("style","null");
00193     theStyle->SetPalette(1,0);
00194     theStyle->SetCanvasDefH(1200); //Height of canvas
00195     theStyle->SetCanvasDefW(1600); //Width of canvas
00196 
00197     gStyle = theStyle;
00198 /*
00199     TCanvas * c1 = new TCanvas("c1","graph",1);
00200     c1->Divide(2,2);
00201     c1->cd(1);
00202     CASTORMeans->Draw();
00203     c1->SaveAs(name1.c_str());   
00204 
00205     theStyle->SetOptStat("n");
00206     gStyle = theStyle;
00207 
00208     TCanvas * c2 = new TCanvas("c2","graph",1);
00209  //   c2->Divide(2,2);
00210     c2->cd(1);
00211     dephist->Draw();
00212     dephist->SetDrawOption("colz");
00213     //c2->cd(2);
00214     //dephist[1]->Draw();
00215     //dephist[1]->SetDrawOption("colz");
00216     //c2->cd(3);
00217     //dephist[2]->Draw();
00218     //dephist[2]->SetDrawOption("colz");
00219     //c2->cd(4);
00220     //dephist[3]->Draw();
00221     //dephist[3]->SetDrawOption("colz");
00222     c2->SaveAs(name2.c_str());
00223 */
00224     std::cout << "Writing ROOT file... ";
00225     theFile->Close();
00226     std::cout << "ROOT file closed.\n";
00227 }
00228 
00229 // ------------ method called to for each event  ------------
00230 void
00231 CastorPedestalsAnalysis::analyze(const edm::Event& e, const edm::EventSetup& iSetup)
00232 {
00233    using namespace edm;
00234    using namespace std;
00235 
00236    edm::Handle<CastorDigiCollection> castor;
00237    e.getByLabel(castorDigiCollectionTag, castor);
00238 
00239    edm::ESHandle<CastorDbService> conditions;
00240    iSetup.get<CastorDbRecord>().get(conditions);
00241 
00242    const CastorQIEShape* shape = conditions->getCastorShape();
00243 
00244    if(firsttime)
00245    {
00246       runnum = e.id().run();
00247       std::string runnum_string;
00248       std::stringstream tempstringout;
00249       tempstringout << runnum;
00250       runnum_string = tempstringout.str();
00251       ROOTfilename = runnum_string + "-peds_ADC.root";
00252       pedsADCfilename = runnum_string + "-peds_ADC.txt";
00253       pedsfCfilename = runnum_string + "-peds_fC.txt";
00254       widthsADCfilename = runnum_string + "-widths_ADC.txt";
00255       widthsfCfilename = runnum_string + "-widths_fC.txt";
00256       XMLfilename = runnum_string + "-peds_ADC_complete.xml"; 
00257       XMLtag = "Castor_pedestals_" + runnum_string;
00258 
00259       theFile = new TFile(ROOTfilename.c_str(), "RECREATE");
00260       theFile->cd();
00261       // Create sub-directories
00262       theFile->mkdir("CASTOR"); 
00263       theFile->cd();
00264 
00265       CASTORMeans = new TH1F("All Ped Means CASTOR","All Ped Means CASTOR", 100, 0, 9);
00266       CASTORWidths = new TH1F("All Ped Widths CASTOR","All Ped Widths CASTOR", 100, 0, 3);
00267 
00268       dephist = new TH2F("Pedestals (ADC)","All Castor",14, 0., 14.5, 16, .5, 16.5);
00269      // dephist[0] = new TH2F("Pedestals (ADC)","Depth 1",89, -44, 44, 72, .5, 72.5);
00270      // dephist[1] = new TH2F("Pedestals (ADC)","Depth 2",89, -44, 44, 72, .5, 72.5);
00271      // dephist[2] = new TH2F("Pedestals (ADC)","Depth 3",89, -44, 44, 72, .5, 72.5);
00272      // dephist[3] = new TH2F("Pedestals (ADC)","Depth 4",89, -44, 44, 72, .5, 72.5);
00273 
00274       edm::ESHandle<CastorElectronicsMap> refEMap;
00275       iSetup.get<CastorElectronicsMapRcd>().get(refEMap);
00276       const CastorElectronicsMap* myRefEMap = refEMap.product();
00277       std::vector<HcalGenericDetId> listEMap = myRefEMap->allPrecisionId();
00278       for (std::vector<HcalGenericDetId>::const_iterator it = listEMap.begin(); it != listEMap.end(); it++)
00279       {     
00280          HcalGenericDetId mygenid(it->rawId());
00281          if(mygenid.isHcalCastorDetId())
00282          {
00283             NewPedBunch a;
00284             HcalCastorDetId chanid(mygenid.rawId());
00285             a.detid = chanid;
00286             a.usedflag = false;
00287             string type;
00288                         type = "CASTOR";
00289             for(int i = 0; i != 4; i++)
00290             {
00291                a.cap[i] = 0;
00292                a.capfc[i] = 0;
00293                for(int j = 0; j != 4; j++)
00294                {
00295                   a.sig[i][j] = 0;
00296                   a.sigfc[i][j] = 0;
00297                   a.prod[i][j] = 0;
00298                   a.prodfc[i][j] = 0;
00299                   a.num[i][j] = 0;
00300                }
00301             }
00302             Bunches.push_back(a);
00303          }
00304       }
00305       firsttime = false;
00306    }
00307 
00308    std::vector<NewPedBunch>::iterator bunch_it;
00309 
00310    for(CastorDigiCollection::const_iterator j = castor->begin(); j != castor->end(); j++)
00311    {
00312       const CastorDataFrame digi = (const CastorDataFrame)(*j);
00313       for(bunch_it = Bunches.begin(); bunch_it != Bunches.end(); bunch_it++)
00314          if(bunch_it->detid.rawId() == digi.id().rawId()) break;
00315       bunch_it->usedflag = true;
00316       for(int ts = firstTS; ts != lastTS+1; ts++)
00317       {
00318          const CastorQIECoder* coder = conditions->getCastorCoder(digi.id().rawId());
00319          bunch_it->num[digi.sample(ts).capid()][digi.sample(ts).capid()] += 1;
00320          bunch_it->cap[digi.sample(ts).capid()] += digi.sample(ts).adc();
00321          double charge1 = coder->charge(*shape, digi.sample(ts).adc(), digi.sample(ts).capid());
00322          bunch_it->capfc[digi.sample(ts).capid()] += charge1;
00323          bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts).capid()] += (digi.sample(ts).adc() * digi.sample(ts).adc());
00324          bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts).capid()] += charge1 * charge1;
00325          if((ts+1 < digi.size()) && (ts+1 < lastTS)){
00326             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += digi.sample(ts).adc()*digi.sample(ts+1).adc();
00327             double charge2 = coder->charge(*shape, digi.sample(ts+1).adc(), digi.sample(ts+1).capid());
00328             bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += charge1*charge2;
00329             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+1).capid()] += 1;
00330          }
00331          if((ts+2 < digi.size()) && (ts+2 < lastTS)){
00332             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += digi.sample(ts).adc()*digi.sample(ts+2).adc();
00333             double charge2 = coder->charge(*shape, digi.sample(ts+2).adc(), digi.sample(ts+2).capid());
00334             bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += charge1*charge2;
00335             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+2).capid()] += 1;
00336          }
00337          if((ts+3 < digi.size()) && (ts+3 < lastTS)){
00338             bunch_it->prod[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += digi.sample(ts).adc()*digi.sample(ts+3).adc();
00339             double charge2 = coder->charge(*shape, digi.sample(ts+3).adc(), digi.sample(ts+3).capid());
00340             bunch_it->prodfc[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += charge1*charge2;
00341             bunch_it->num[digi.sample(ts).capid()][digi.sample(ts+3).capid()] += 1;
00342          }
00343       }
00344    }
00345 
00346 
00347 //this is the last brace
00348 }
00349 
00350 //define this as a plug-in
00351 DEFINE_FWK_MODULE(CastorPedestalsAnalysis);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/CalibCalorimetry/CastorCalib/plugins/CastorPedestalsAnalysis.cc
