Reference Manual

Go to the documentation of this file.
00001 
00002 // Package:    ZDCDigiStudy
00003 // Class:      ZDCDigiStudy
00004 //
00005 /*
00006  Description: 
00007               This code has been developed to be a check for the ZDC sim. In 2009, it was found that the ZDC Simulation was unrealistic and needed repair. The aim of this code is to show the user the input and output of a ZDC MinBias simulation.
00008 
00009  Implementation:
00010       First a MinBias simulation should be run, it could be pythia,hijin,or hydjet. This will output a .root file which should have information about recoGenParticles, hcalunsuppresseddigis. Use this .root file as the input into the cfg.py which is found in the main directory of this package. This output will be another .root file which is meant to be viewed in a TBrowser.
00011 
00012 */
00013 //
00014 // Original Author: Jaime Gomez (U. of Maryland) with SIGNIFICANT assistance of Dr. Jefferey Temple (U. of Maryland) 
00015 //   
00016 //
00017 //         Created:  Summer 2012
00019 
00020 
00021 
00022 
00023 
00024 #include "Validation/HcalDigis/interface/ZDCDigiStudy.h"
00025 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
00026 
00027 #include "FWCore/Utilities/interface/Exception.h"
00028 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00029 
00030 ZDCDigiStudy::ZDCDigiStudy(const edm::ParameterSet& ps) {
00031 
00032   zdcHits = ps.getUntrackedParameter<std::string>("HitCollection","ZdcHits");
00033   outFile_ = ps.getUntrackedParameter<std::string>("outputFile", "zdcHitStudy.root");
00034   verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
00035   checkHit_= true;
00036 
00037   edm::LogInfo("ZDCDigiStudy") 
00038     //std::cout
00039     << "   Hits: "
00040     << zdcHits << " / "<< checkHit_ 
00041     << "   Output: " << outFile_;
00042 
00043   dbe_ = edm::Service<DQMStore>().operator->();
00044   if (dbe_) {
00045     if (verbose_) {
00046       dbe_->setVerbose(1);
00047       sleep (3);
00048       dbe_->showDirStructure();
00049     } else {
00050       dbe_->setVerbose(0);
00051     }
00052   }
00053 }
00054 
00055 ZDCDigiStudy::~ZDCDigiStudy() {}
00056 
00057 void ZDCDigiStudy::beginJob() {
00058   if (dbe_) {
00059     dbe_->setCurrentFolder("ZDCDigiValidation");
00060     //Histograms for Hits
00062 //# Below we are filling the histograms made in the .h file. The syntax is as follows:                                      #
00063 //# plot_code_name = dbe_->TypeofPlot[(1,2,3)-D,(F,I,D)]("Name as it will appear","Title",axis options);                    #
00064 //# They will be stored in the TFile subdirectory set by :    dbe_->setCurrentFolder("FolderIwant")                         #
00065 //# axis options are like (#ofbins,min,max)                                                                                 #
00067 
00068 
00069     if (checkHit_) {
00070 
00071 
00073 
00075       dbe_->setCurrentFolder("ZDCDigiValidation/ZDC_Digis/1D_fC");
00076       meZdcfCPHAD = dbe_->book1D("PHAD_TotalfC","PZDC_HAD_TotalfC",1000,-50,950);
00077       meZdcfCPHAD->setAxisTitle("Counts",2);
00078       meZdcfCPHAD->setAxisTitle("fC",1);
00080       meZdcfCPTOT = dbe_->book1D("PZDC_TotalfC","PZDC_TotalfC",1000,-50,950);
00081       meZdcfCPTOT->setAxisTitle("Counts",2);
00082       meZdcfCPTOT->setAxisTitle("fC",1);
00084       meZdcfCNHAD = dbe_->book1D("NHAD_TotalfC","NZDC_HAD_TotalfC",1000,-50,950);
00085       meZdcfCNHAD->setAxisTitle("Counts",2);
00086       meZdcfCNHAD->setAxisTitle("fC",1);
00088       meZdcfCNTOT = dbe_->book1D("NZDC_TotalfC","NZDC_TotalfC",1000,-50,950);
00089       meZdcfCNTOT->setAxisTitle("Counts",2);
00090       meZdcfCNTOT->setAxisTitle("fC",1);
00092 
00094      dbe_->setCurrentFolder("ZDCDigiValidation/ZDC_Digis/fCvsTS/PZDC");
00095      
00097      meZdcPEM1fCvsTS = dbe_->book1D("PEM1_fCvsTS","P-EM1_AveragefC_vsTS",10,0,9);
00098      meZdcPEM1fCvsTS->setAxisTitle("fC",2);
00099      meZdcPEM1fCvsTS->setAxisTitle("TS",1);
00102      meZdcPEM2fCvsTS = dbe_->book1D("PEM2_fCvsTS","P-EM2_AveragefC_vsTS",10,0,9);
00103      meZdcPEM2fCvsTS->setAxisTitle("fC",2);
00104      meZdcPEM2fCvsTS->setAxisTitle("TS",1);
00107      meZdcPEM3fCvsTS = dbe_->book1D("PEM3_fCvsTS","P-EM3_AveragefC_vsTS",10,0,9);
00108      meZdcPEM3fCvsTS->setAxisTitle("fC",2);
00109      meZdcPEM3fCvsTS->setAxisTitle("TS",1);
00112      meZdcPEM4fCvsTS = dbe_->book1D("PEM4_fCvsTS","P-EM4_AveragefC_vsTS",10,0,9);
00113      meZdcPEM4fCvsTS->setAxisTitle("fC",2);
00114      meZdcPEM4fCvsTS->setAxisTitle("TS",1);
00117      meZdcPEM5fCvsTS = dbe_->book1D("PEM5_fCvsTS","P-EM5_AveragefC_vsTS",10,0,9);
00118      meZdcPEM5fCvsTS->setAxisTitle("fC",2);
00119      meZdcPEM5fCvsTS->setAxisTitle("TS",1);
00122      meZdcPHAD1fCvsTS = dbe_->book1D("PHAD1_fCvsTS","P-HAD1_AveragefC_vsTS",10,0,9);
00123      meZdcPHAD1fCvsTS->setAxisTitle("fC",2);
00124      meZdcPHAD1fCvsTS->setAxisTitle("TS",1);
00127      meZdcPHAD2fCvsTS = dbe_->book1D("PHAD2_fCvsTS","P-HAD2_AveragefC_vsTS",10,0,9);
00128      meZdcPHAD2fCvsTS->setAxisTitle("fC",2);
00129      meZdcPHAD2fCvsTS->setAxisTitle("TS",1);
00132      meZdcPHAD3fCvsTS = dbe_->book1D("PHAD3_fCvsTS","P-HAD3_AveragefC_vsTS",10,0,9);
00133      meZdcPHAD3fCvsTS->setAxisTitle("fC",2);
00134      meZdcPHAD3fCvsTS->setAxisTitle("TS",1);
00137      meZdcPHAD4fCvsTS = dbe_->book1D("PHAD4_fCvsTS","P-HAD4_AveragefC_vsTS",10,0,9);
00138      meZdcPHAD4fCvsTS->setAxisTitle("fC",2);
00139      meZdcPHAD4fCvsTS->setAxisTitle("TS",1);
00141      dbe_->setCurrentFolder("ZDCDigiValidation/ZDC_Digis/fCvsTS/NZDC");
00142      
00144      meZdcNEM1fCvsTS = dbe_->book1D("NEM1_fCvsTS","N-EM1_AveragefC_vsTS",10,0,9);
00145      meZdcNEM1fCvsTS->setAxisTitle("fC",2);
00146      meZdcNEM1fCvsTS->setAxisTitle("TS",1);
00149      meZdcNEM2fCvsTS = dbe_->book1D("NEM2_fCvsTS","N-EM2_AveragefC_vsTS",10,0,9);
00150      meZdcNEM2fCvsTS->setAxisTitle("fC",2);
00151      meZdcNEM2fCvsTS->setAxisTitle("TS",1);
00154      meZdcNEM3fCvsTS = dbe_->book1D("NEM3_fCvsTS","N-EM3_AveragefC_vsTS",10,0,9);
00155      meZdcNEM3fCvsTS->setAxisTitle("fC",2);
00156      meZdcNEM3fCvsTS->setAxisTitle("TS",1);
00159      meZdcNEM4fCvsTS = dbe_->book1D("NEM4_fCvsTS","N-EM4_AveragefC_vsTS",10,0,9);
00160      meZdcNEM4fCvsTS->setAxisTitle("fC",2);
00161      meZdcNEM4fCvsTS->setAxisTitle("TS",1);
00164      meZdcNEM5fCvsTS = dbe_->book1D("NEM5_fCvsTS","N-EM5_AveragefC_vsTS",10,0,9);
00165      meZdcNEM5fCvsTS->setAxisTitle("fC",2);
00166      meZdcNEM5fCvsTS->setAxisTitle("TS",1);
00169      meZdcNHAD1fCvsTS = dbe_->book1D("NHAD1_fCvsTS","N-HAD1_AveragefC_vsTS",10,0,9);
00170      meZdcNHAD1fCvsTS->setAxisTitle("fC",2);
00171      meZdcNHAD1fCvsTS->setAxisTitle("TS",1);
00174      meZdcNHAD2fCvsTS = dbe_->book1D("NHAD2_fCvsTS","N-HAD2_AveragefC_vsTS",10,0,9);
00175      meZdcNHAD2fCvsTS->setAxisTitle("fC",2);
00176      meZdcNHAD2fCvsTS->setAxisTitle("TS",1);
00179      meZdcNHAD3fCvsTS = dbe_->book1D("NHAD3_fCvsTS","N-HAD3_AveragefC_vsTS",10,0,9);
00180      meZdcNHAD3fCvsTS->setAxisTitle("fC",2);
00181      meZdcNHAD3fCvsTS->setAxisTitle("TS",1);
00184      meZdcNHAD4fCvsTS = dbe_->book1D("NHAD4_fCvsTS","N-HAD4_AveragefC_vsTS",10,0,9);
00185      meZdcNHAD4fCvsTS->setAxisTitle("fC",2);
00186      meZdcNHAD4fCvsTS->setAxisTitle("TS",1);
00188 
00190     dbe_->setCurrentFolder("ZDCDigiValidation/ZDC_Digis/2D_EMvHAD");
00193     meZdcfCPEMvHAD = dbe_->book2D("PEMvPHAD","PZDC_EMvHAD",1000,-25,1000,1000,-25,1000);
00194     meZdcfCPEMvHAD->setAxisTitle("SumEM_fC",2);
00195     meZdcfCPEMvHAD->setAxisTitle("SumHAD_fC",1);
00196     meZdcfCPEMvHAD->getTH2F()->SetOption("colz");
00198     meZdcfCNEMvHAD = dbe_->book2D("NEMvNHAD","NZDC_EMvHAD",1000,-25,1000,1000,-25,1000);
00199     meZdcfCNEMvHAD->setAxisTitle("SumEM_fC",2);
00200     meZdcfCNEMvHAD->setAxisTitle("SumHAD_fC",1);
00201     meZdcfCNEMvHAD->getTH2F()->SetOption("colz");
00203 
00204     }
00205   }
00206 }
00207 
00208 void ZDCDigiStudy::endJob() {
00209   if (dbe_ && outFile_.size() > 0) dbe_->save(outFile_);
00210 }
00211 
00212 //void ZDCDigiStudy::analyze(const edm::Event& e, const edm::EventSetup& ) {
00213 void ZDCDigiStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup ) {
00215 
00216    using namespace edm;
00217    bool gotZDCDigis=true;
00218    
00219 
00220    Handle<ZDCDigiCollection> zdchandle; 
00221    if (!(iEvent.getByLabel("simHcalUnsuppressedDigis",zdchandle)))
00222    {
00223     gotZDCDigis=false; //this is a boolean set up to check if there are ZDCDigis in the input root file
00224    }
00225    if (!(zdchandle.isValid()))
00226    {
00227    gotZDCDigis=false; //if it is not there, leave it false
00228    }
00229   
00230    double totalPHADCharge=0;
00231    double totalNHADCharge=0;
00232    double totalPEMCharge=0;
00233    double totalNEMCharge=0;
00234    double totalPCharge=0;
00235    double totalNCharge=0;
00236    
00237 
00239       if (gotZDCDigis==true){
00240     for (ZDCDigiCollection::const_iterator zdc = zdchandle->begin();
00241         zdc!=zdchandle->end();
00242         ++zdc)
00243      {
00244        const ZDCDataFrame digi = (const ZDCDataFrame)(*zdc);
00245        //std::cout <<"CHANNEL = "<<zdc->id().channel()<<std::endl;
00246        
00247 
00248 
00250 
00251        if (digi.id().section()==2){  // require HAD
00252        if (digi.id().zside()==1)
00253            { // require POS
00254          for (int i=0;i<digi.size();++i) // loop over all 10 TS because each digi has 10 entries
00255          {
00256            if (digi.id().channel()==1){ //here i specify PHAD1
00257              meZdcPHAD1fCvsTS->Fill(i,digi.sample(i).nominal_fC()); //filling the plot name with the nominal fC value for each TS
00258              if (i==0) meZdcPHAD1fCvsTS->Fill(-1,1);  // on first iteration of loop, increment underflow bin
00259            }//NEW AVERAGE Thingy
00260            if (digi.id().channel()==2){
00261              meZdcPHAD2fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00262              if (i==0) meZdcPHAD2fCvsTS->Fill(-1,1);
00263                                       }
00264            if (digi.id().channel()==3){
00265              meZdcPHAD3fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00266              if (i==0) meZdcPHAD3fCvsTS->Fill(-1,1);
00267                                        }
00268            if (digi.id().channel()==4){
00269              meZdcPHAD4fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00270              if (i==0) meZdcPHAD4fCvsTS->Fill(-1,1);
00271                                       }
00272             totalPHADCharge+=digi.sample(i).nominal_fC(); //here i am looking for the total charge in PHAD so i sum over every TS and channel and add it up
00273          } // loop over all (10) TS for the given digi
00274            }
00275        else {
00276          for (int i=0; i<digi.size();++i)
00277            {
00278              if (digi.id().channel()==1){
00279                meZdcNHAD1fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00280                if (i==0) meZdcNHAD1fCvsTS->Fill(-1,1);
00281                                         }
00282              if (digi.id().channel()==2){
00283                meZdcNHAD2fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00284                if (i==0) meZdcNHAD2fCvsTS->Fill(-1,1);
00285                                          }
00286              if (digi.id().channel()==3){
00287                meZdcNHAD3fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00288                if (i==0) meZdcNHAD3fCvsTS->Fill(-1,1);
00289                                         }
00290              if (digi.id().channel()==4){
00291                meZdcNHAD4fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00292                if (i==0) meZdcNHAD4fCvsTS->Fill(-1,1);
00293                                          }
00294             totalNHADCharge+=digi.sample(i).nominal_fC();
00295            }
00296        }
00297        }
00299        if (digi.id().section()==1){//require EM....here i do the smae thing that i did above but now for P/N EM sections
00300          if (digi.id().zside()==1)
00301            {//require pos
00302              for (int i=0;i<digi.size();++i)
00303                {
00304                  if (digi.id().channel()==1){
00305                    meZdcPEM1fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00306                    if (i==0) meZdcPEM1fCvsTS->Fill(-1,1);
00307                                              }
00308                  if (digi.id().channel()==2){
00309                    meZdcPEM2fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00310                    if (i==0) meZdcPEM2fCvsTS->Fill(-1,1);
00311                                              }
00312                  if (digi.id().channel()==3){
00313                    meZdcPEM3fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00314                    if (i==0) meZdcPEM3fCvsTS->Fill(-1,1);
00315                                              }
00316                  if (digi.id().channel()==4){
00317                    meZdcPEM4fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00318                    if (i==0) meZdcPEM4fCvsTS->Fill(-1,1);
00319                                             }
00320                  if (digi.id().channel()==5){
00321                    meZdcPEM5fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00322                    if (i==0) meZdcPEM5fCvsTS->Fill(-1,1);
00323                                              }
00324                 totalPEMCharge+=digi.sample(i).nominal_fC();
00325                }
00326            }
00327          else {
00328            for (int i=0;i<digi.size();++i)
00329              {
00330                if (digi.id().channel()==1){
00331                  meZdcNEM1fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00332                  if (i==0) meZdcNEM1fCvsTS->Fill(-1,1);
00333                                           }
00334                if (digi.id().channel()==2){
00335                  meZdcNEM2fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00336                  if (i==0) meZdcNEM2fCvsTS->Fill(-1,1);
00337                                           }
00338                if (digi.id().channel()==3){
00339                  meZdcNEM3fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00340                  if (i==0) meZdcNEM3fCvsTS->Fill(-1,1);
00341                                           }
00342                if (digi.id().channel()==4){
00343                  meZdcNEM4fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00344                  if (i==0) meZdcNEM4fCvsTS->Fill(-1,1);
00345                                           }
00346                if (digi.id().channel()==5){
00347                  meZdcNEM5fCvsTS->Fill(i,digi.sample(i).nominal_fC());
00348                  if (i==0) meZdcNEM5fCvsTS->Fill(-1,1);
00349                                           }
00350                totalNEMCharge+=digi.sample(i).nominal_fC();
00351              }
00352          }
00353        }     
00354 
00355        totalPCharge=totalPHADCharge+totalPEMCharge;
00356        totalNCharge=totalNHADCharge+totalNEMCharge;
00357 
00358        /*       std::cout <<"CHANNEL = "<<digi.id().channel()<<std::endl;
00359        for (int i=0;i<digi.size();++i)
00360          std::cout <<"SAMPLE = "<<i<<"  ADC = "<<digi.sample(i).adc()<<" fC =  "<<digi.sample(i).nominal_fC()<<std::endl;
00361        */
00362        //  digi[i] should be the sample as digi.sample(i), I think
00363      } // loop on all (22) ZDC digis
00364      }
00366 
00367   
00368 
00369  
00370    // Now fill total charge histogram
00371    meZdcfCPEMvHAD->Fill(totalPHADCharge,totalPEMCharge);
00372    meZdcfCNEMvHAD->Fill(totalNHADCharge,totalNEMCharge);
00373    meZdcfCPHAD->Fill(totalPHADCharge);
00374    meZdcfCNHAD->Fill(totalNHADCharge);
00375    meZdcfCNTOT->Fill(totalNCharge);
00376    meZdcfCPTOT->Fill(totalPCharge);
00377 
00378 }
00379 
00381 
00382 void ZDCDigiStudy::endRun(const edm::Run& run, const edm::EventSetup& c)
00383 {
00384   int nevents=(meZdcPHAD1fCvsTS->getTH1F())->GetBinContent(0); //grab the number of digis that were read in and stored in the underflow bin, and call them Nevents
00385   (meZdcPHAD1fCvsTS->getTH1F())->Scale(1./nevents);  // divide histogram by nevents thereby creating an average..it was done this way so that in DQM when everything is done in parallel and added at the end then the average will add appropriately
00386 
00387   int nevents1=(meZdcPHAD2fCvsTS->getTH1F())->GetBinContent(0);
00388   (meZdcPHAD2fCvsTS->getTH1F())->Scale(1./nevents1);
00389 
00390   int nevents2=(meZdcPHAD3fCvsTS->getTH1F())->GetBinContent(0);
00391   (meZdcPHAD3fCvsTS->getTH1F())->Scale(1./nevents2);
00392 
00393   int nevents3=(meZdcPHAD4fCvsTS->getTH1F())->GetBinContent(0);
00394   (meZdcPHAD4fCvsTS->getTH1F())->Scale(1./nevents3);
00395 
00396   int nevents4=(meZdcNHAD1fCvsTS->getTH1F())->GetBinContent(0);
00397   (meZdcNHAD1fCvsTS->getTH1F())->Scale(1./nevents4);
00398 
00399   int nevents5=(meZdcNHAD2fCvsTS->getTH1F())->GetBinContent(0);
00400   (meZdcNHAD2fCvsTS->getTH1F())->Scale(1./nevents5);
00401 
00402   int nevents6=(meZdcNHAD3fCvsTS->getTH1F())->GetBinContent(0);
00403   (meZdcNHAD3fCvsTS->getTH1F())->Scale(1./nevents6);
00404 
00405   int nevents7=(meZdcNHAD4fCvsTS->getTH1F())->GetBinContent(0);
00406   (meZdcNHAD4fCvsTS->getTH1F())->Scale(1./nevents7);
00407 
00408   int nevents8=(meZdcPEM1fCvsTS->getTH1F())->GetBinContent(0);
00409   (meZdcPEM1fCvsTS->getTH1F())->Scale(1./nevents8);
00410 
00411   int nevents9=(meZdcPEM2fCvsTS->getTH1F())->GetBinContent(0);
00412   (meZdcPEM2fCvsTS->getTH1F())->Scale(1./nevents9);
00413 
00414   int nevents10=(meZdcPEM3fCvsTS->getTH1F())->GetBinContent(0);
00415   (meZdcPEM3fCvsTS->getTH1F())->Scale(1./nevents10);
00416 
00417   int nevents11=(meZdcPEM4fCvsTS->getTH1F())->GetBinContent(0);
00418   (meZdcPEM4fCvsTS->getTH1F())->Scale(1./nevents11);
00419 
00420   int nevents12=(meZdcPEM5fCvsTS->getTH1F())->GetBinContent(0);
00421   (meZdcPEM5fCvsTS->getTH1F())->Scale(1./nevents12);
00422 
00423   int nevents13=(meZdcNEM1fCvsTS->getTH1F())->GetBinContent(0);
00424   (meZdcNEM1fCvsTS->getTH1F())->Scale(1./nevents13);
00425 
00426   int nevents14=(meZdcNEM2fCvsTS->getTH1F())->GetBinContent(0);
00427   (meZdcNEM2fCvsTS->getTH1F())->Scale(1./nevents14);
00428 
00429   int nevents15=(meZdcNEM3fCvsTS->getTH1F())->GetBinContent(0);
00430   (meZdcNEM3fCvsTS->getTH1F())->Scale(1./nevents15);
00431 
00432   int nevents16=(meZdcNEM4fCvsTS->getTH1F())->GetBinContent(0);
00433   (meZdcNEM4fCvsTS->getTH1F())->Scale(1./nevents16);
00434 
00435   int nevents17=(meZdcNEM5fCvsTS->getTH1F())->GetBinContent(0);
00436   (meZdcNEM5fCvsTS->getTH1F())->Scale(1./nevents17);
00437 
00438 
00439 
00440 }
00441 
00442 //define this as a plug-in
00443 DEFINE_FWK_MODULE(ZDCDigiStudy);
00444 
00445 
00446 
00447 
00448 
00449 
00450 
00451 
00452 
00453
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Validation/HcalDigis/src/ZDCDigiStudy.cc
