Reference Manual

Go to the documentation of this file.
00001 // system include files
00002 #include <memory>
00003 #include <string>
00004 #include <iostream>
00005 
00006 // user include files
00007 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00008 #include "Calibration/HcalCalibAlgos/interface/Analyzer_minbias.h"
00009 #include "DataFormats/Provenance/interface/Provenance.h"
00010 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00011 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00012 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00013 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00014 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00015 #include "DataFormats/L1GlobalTrigger/interface/L1GtfeWord.h"
00016 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
00017 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerRecord.h"
00018 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
00019 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
00020 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
00021 
00022 #include "TFile.h"
00023 #include "TH1.h"
00024 #include "TH2.h"
00025 #include <fstream>
00026 #include <sstream>
00027 
00028 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
00029 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
00030 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
00031 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
00032 using namespace std;
00033 using namespace reco;
00034 //
00035 // constructors and destructor
00036 //
00037 namespace cms{
00038 Analyzer_minbias::Analyzer_minbias(const edm::ParameterSet& iConfig)
00039 {
00040   // get name of output file with histogramms
00041   fOutputFileName = iConfig.getUntrackedParameter<string>("HistOutFile"); 
00042   // get names of modules, producing object collections
00043   
00044   hbherecoMB = iConfig.getParameter<edm::InputTag>("hbheInputMB");
00045   horecoMB   = iConfig.getParameter<edm::InputTag>("hoInputMB");
00046   hfrecoMB   = iConfig.getParameter<edm::InputTag>("hfInputMB");
00047   
00048   hbherecoNoise = iConfig.getParameter<edm::InputTag>("hbheInputNoise");
00049   horecoNoise   = iConfig.getParameter<edm::InputTag>("hoInputNoise");
00050   hfrecoNoise   = iConfig.getParameter<edm::InputTag>("hfInputNoise");
00051   
00052   theRecalib = iConfig.getParameter<bool>("Recalib"); 
00053    
00054 //
00055 //
00056   for(int i=0; i<73; i++)
00057   {
00058      for(int j=0; j<43; j++)
00059      {
00060         noise_min[i][j] = 0.;
00061         noise_pl[i][j] = 0.;
00062      } 
00063   }
00064 //
00065 //
00066      
00067 }
00068 
00069 Analyzer_minbias::~Analyzer_minbias()
00070 {
00071  
00072    // do anything here that needs to be done at desctruction time
00073    // (e.g. close files, deallocate resources etc.)
00074 
00075 }
00076 
00077 void Analyzer_minbias::beginRun( const edm::Run& r, const edm::EventSetup& iSetup)
00078 {
00079   nevent_run = 0;
00080 }
00081 void Analyzer_minbias::endRun( const edm::Run& r, const edm::EventSetup& iSetup)
00082 {
00083  std::cout<<" Runnumber "<<r.run()<<" Nevents  "<<nevent_run<<std::endl;
00084 }
00085 
00086 void Analyzer_minbias::beginJob()
00087 {
00088    
00089    hOutputFile   = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
00090    
00091    myTree = new TTree("RecJet","RecJet Tree");
00092    myTree->Branch("mydet",  &mydet, "mydet/I");
00093    myTree->Branch("mysubd",  &mysubd, "mysubd/I");
00094    myTree->Branch("depth",  &depth, "depth/I");
00095    myTree->Branch("ieta",  &ieta, "ieta/I");
00096    myTree->Branch("iphi",  &iphi, "iphi/I");
00097    myTree->Branch("eta",  &eta, "eta/F");
00098    myTree->Branch("phi",  &phi, "phi/F");
00099    
00100    myTree->Branch("mom0_MB",  &mom0_MB, "mom0_MB/F");
00101    myTree->Branch("mom1_MB",  &mom1_MB, "mom1_MB/F");
00102    myTree->Branch("mom2_MB",  &mom2_MB, "mom2_MB/F");
00103    myTree->Branch("mom4_MB",  &mom4_MB, "mom4_MB/F");
00104    
00105    myTree->Branch("mom0_Noise",  &mom0_Noise, "mom0_Noise/F");
00106    myTree->Branch("mom1_Noise",  &mom1_Noise, "mom1_Noise/F");
00107    myTree->Branch("mom2_Noise",  &mom2_Noise, "mom2_Noise/F");
00108    myTree->Branch("mom4_Noise",  &mom4_Noise, "mom4_Noise/F");
00109    
00110    myTree->Branch("mom0_Diff",  &mom0_Diff, "mom0_Diff/F");
00111    myTree->Branch("mom1_Diff",  &mom1_Diff, "mom1_Diff/F");
00112    myTree->Branch("mom2_Diff",  &mom2_Diff, "mom2_Diff/F");
00113 
00114    myTree->Branch("occup",  &occup, "occup/F");
00115    
00116    std::cout<<" Before ordering Histos "<<std::endl;
00117   
00118    char str0[15];
00119    char str1[15];
00120 
00121    char str10[15];
00122    char str11[15];
00123 
00124    int k=0;
00125    nevent = 0;
00126 // Size of collections
00127 
00128    hHBHEsize_vs_run = new TH2F("hHBHEsize_vs_run","hHBHEsize_vs_run",500,111500.,112000.,6101,-100.5,6000.5);
00129    hHFsize_vs_run = new TH2F("hHFsize_vs_run","hHFsize_vs_run",500,111500.,112000.,6101,-100.5,6000.5);
00130    
00131    for(int i=1;i<73;i++){
00132     for(int j=1;j<43;j++){
00133 
00134        meannoise_pl[i][j] = 0.;
00135        meannoise_min[i][j] = 0.;
00136 
00137 //     for(int l=1;l<5;l++){
00138         k = i*1000+j;
00139         sprintf(str0,"mpl%d",k);
00140         sprintf(str1,"mmin%d",k);
00141 
00142         sprintf(str10,"vpl%d",k);
00143         sprintf(str11,"vmin%d",k);
00144 //      cout<<" "<<i<<" "<<j<<endl;
00145     if( j < 30 )
00146     {
00147 // first order moment
00148     hCalo1[i][j] = new TH1F(str0, "h0", 320, -10., 10.);
00149     hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
00150 
00151 // second order moment
00152     hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 20.);
00153     hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 20.);
00154     }
00155       else
00156       {
00157 // HF
00158 // first order moment
00159 //   cout<<" "<<i<<" "<<j<<" "<<k<<endl;
00160    if(j < 40)
00161    {
00162     hCalo1[i][j] = new TH1F(str0, "h0", 320, -10., 10.);
00163     hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
00164 //
00165 // second order moment
00166     hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 40.);
00167     hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 40.);
00168    }
00169      else
00170      {
00171     hCalo1[i][j] = new TH1F(str0,"h0" , 320, -10., 10.);
00172     hCalo2[i][j] = new TH1F(str1, "h1", 320, -10., 10.);
00173 
00174 // second order moment
00175     hCalo1mom2[i][j] = new TH1F(str10, "h10", 320, 0., 120.);
00176     hCalo2mom2[i][j] = new TH1F(str11, "h11", 320, 0., 120.);
00177 
00178      }
00179     } // HE/HF boundary
00180 //     } // l
00181     } // j
00182    } // i
00183 
00184 
00185     hbheNoiseE = new TH1F("hbheNoiseE","hbheNoiseE", 320, -10., 10.);
00186     hfNoiseE = new TH1F("hfNoiseE","hfNoiseE", 320, -10., 10.);
00187     hbheSignalE = new TH1F("hbheSignalE","hbheSignalE", 320, -10., 10.);
00188     hfSignalE = new TH1F("hfSignalE","hfSignalE", 320, -10., 10.);
00189     
00190 
00191    std::cout<<" After ordering Histos "<<std::endl;
00192 
00193    std::string ccc = "noise_0.dat";
00194 
00195    myout_hcal = new ofstream(ccc.c_str());
00196    if(!myout_hcal) cout << " Output file not open!!! "<<endl;
00197 
00198 //
00199    for (int i=0; i<5;i++)
00200    {
00201     for (int j=0; j<5;j++)
00202     {
00203      for (int k=0; k<73;k++)
00204      {
00205        for (int l=0; l<43;l++)
00206        {
00207         theMBFillDetMapPl0[i][j][k][l] = 0.;
00208         theMBFillDetMapPl1[i][j][k][l] = 0.;
00209         theMBFillDetMapPl2[i][j][k][l] = 0.;
00210         theMBFillDetMapPl4[i][j][k][l] = 0.;
00211 
00212         theMBFillDetMapMin0[i][j][k][l] = 0.;
00213         theMBFillDetMapMin1[i][j][k][l] = 0.;
00214         theMBFillDetMapMin2[i][j][k][l] = 0.;
00215         theMBFillDetMapMin4[i][j][k][l] = 0.;
00216 
00217 
00218         theNSFillDetMapPl0[i][j][k][l] = 0.;
00219         theNSFillDetMapPl1[i][j][k][l] = 0.;
00220         theNSFillDetMapPl2[i][j][k][l] = 0.;
00221         theNSFillDetMapPl4[i][j][k][l] = 0.;
00222 
00223         theNSFillDetMapMin0[i][j][k][l] = 0.;
00224         theNSFillDetMapMin1[i][j][k][l] = 0.;
00225         theNSFillDetMapMin2[i][j][k][l] = 0.;
00226         theNSFillDetMapMin4[i][j][k][l] = 0.;
00227 
00228         theDFFillDetMapPl0[i][j][k][l] = 0.;
00229         theDFFillDetMapPl1[i][j][k][l] = 0.;
00230         theDFFillDetMapPl2[i][j][k][l] = 0.;
00231         theDFFillDetMapMin0[i][j][k][l] = 0.;
00232         theDFFillDetMapMin1[i][j][k][l] = 0.;
00233         theDFFillDetMapMin2[i][j][k][l] = 0.;
00234        }
00235      }  
00236     }
00237    }    
00238     
00239    return ;
00240 }
00241 //
00242 //  EndJob
00243 //
00244 void Analyzer_minbias::endJob()
00245 {
00246    int ii=0;
00247       
00248    for (int i=1; i<5;i++)
00249    {
00250     for (int j=1; j<5;j++)
00251     {
00252      for (int k=1; k<73;k++)
00253      {
00254        for (int l=1; l<43;l++)
00255        {
00256           if(theMBFillDetMapPl0[i][j][k][l] > 0)
00257           { 
00258             mom0_MB = theMBFillDetMapPl0[i][j][k][l];
00259             mom1_MB = theMBFillDetMapPl1[i][j][k][l];
00260             mom2_MB = theMBFillDetMapPl2[i][j][k][l];
00261             mom4_MB = theMBFillDetMapPl4[i][j][k][l];
00262             mom0_Noise = theNSFillDetMapPl0[i][j][k][l];
00263             mom1_Noise = theNSFillDetMapPl1[i][j][k][l];
00264             mom2_Noise = theNSFillDetMapPl2[i][j][k][l];
00265             mom4_Noise = theNSFillDetMapPl4[i][j][k][l];
00266             mom0_Diff = theDFFillDetMapPl0[i][j][k][l];
00267             mom1_Diff = theDFFillDetMapPl1[i][j][k][l];
00268             mom2_Diff = theDFFillDetMapPl2[i][j][k][l];
00269             
00270             mysubd = i;
00271             depth = j;
00272             ieta = l;
00273             iphi = k;
00274             cout<<" Result Plus= "<<mysubd<<" "<<ieta<<" "<<iphi<<" mom0  "<<mom0_MB<<" mom1 "<<mom1_MB<<" mom2 "<<mom2_MB<<endl;
00275             myTree->Fill();
00276             ii++;
00277            } // Pl > 0
00278         
00279          
00280           if(theMBFillDetMapMin0[i][j][k][l] > 0)
00281           { 
00282             mom0_MB = theMBFillDetMapMin0[i][j][k][l];
00283             mom1_MB = theMBFillDetMapMin1[i][j][k][l];
00284             mom2_MB = theMBFillDetMapMin2[i][j][k][l];
00285             mom4_MB = theMBFillDetMapMin4[i][j][k][l];
00286             mom0_Noise = theNSFillDetMapMin0[i][j][k][l];
00287             mom1_Noise = theNSFillDetMapMin1[i][j][k][l];
00288             mom2_Noise = theNSFillDetMapMin2[i][j][k][l];
00289             mom4_Noise = theNSFillDetMapMin4[i][j][k][l];
00290             mom0_Diff = theDFFillDetMapMin0[i][j][k][l];
00291             mom1_Diff = theDFFillDetMapMin1[i][j][k][l];
00292             mom2_Diff = theDFFillDetMapMin2[i][j][k][l];
00293             
00294             
00295             mysubd = i;
00296             depth = j;
00297             ieta = -1*l;
00298             iphi = k;
00299             cout<<" Result Minus= "<<mysubd<<" "<<ieta<<" "<<iphi<<" mom0  "<<mom0_MB<<" mom1 "<<mom1_MB<<" mom2 "<<mom2_MB<<endl;
00300             myTree->Fill();
00301             ii++;
00302             
00303           } // Min>0  
00304       } // ieta
00305      } // iphi  
00306     } // depth
00307    } //subd    
00308    
00309    
00310    
00311    cout<<" Number of cells "<<ii<<endl; 
00312       
00313    hOutputFile->Write();   
00314 
00315    hOutputFile->cd();
00316    
00317    myTree->Write();
00318    
00319    hHBHEsize_vs_run->Write() ;
00320    hHFsize_vs_run->Write() ;
00321 
00322    for(int i=1;i<73;i++){
00323     for(int j=1;j<43;j++){
00324     hCalo1[i][j]->Write();
00325     hCalo2[i][j]->Write();
00326     hCalo1mom2[i][j]->Write();
00327     hCalo2mom2[i][j]->Write();
00328    }
00329   }
00330 
00331     hbheNoiseE->Write() ;
00332     hfNoiseE->Write() ;
00333     hbheSignalE->Write() ;
00334     hfSignalE->Write() ;
00335    
00336 
00337    hOutputFile->Close() ;
00338    
00339    cout<<" File is closed "<<endl;
00340    
00341    return ;
00342 }
00343 
00344 
00345 //
00346 // member functions
00347 //
00348 
00349 // ------------ method called to produce the data  ------------
00350 void
00351 Analyzer_minbias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00352 {
00353 
00354   std::cout<<" Start Analyzer_minbias::analyze "<<nevent<<std::endl;
00355   nevent++; 
00356   nevent_run++;
00357   using namespace edm;
00358 
00359   float rnnum = (float)iEvent.run(); 
00360 /*
00361   std::vector<Provenance const*> theProvenance;
00362   iEvent.getAllProvenance(theProvenance);
00363 
00364   for( std::vector<Provenance const*>::const_iterator ip = theProvenance.begin();
00365                                                       ip != theProvenance.end(); ip++)
00366   {
00367      cout<<" Print all module/label names "<<(**ip).moduleName()<<" "<<(**ip).moduleLabel()<<
00368      " "<<(**ip).productInstanceName()<<endl;
00369   }
00370 */
00371 
00372   edm::ESHandle<L1GtTriggerMenu> menuRcd;
00373   iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
00374   const L1GtTriggerMenu* menu = menuRcd.product();
00375   const AlgorithmMap& bitMap = menu->gtAlgorithmMap();
00376 
00377     edm::Handle<L1GlobalTriggerReadoutRecord> gtRecord;
00378     iEvent.getByLabel("gtDigisAlCaMB", gtRecord);
00379  
00380    if (!gtRecord.isValid()) {
00381 
00382 //     LogDebug("L1GlobalTriggerRecordProducer")
00383 //       << "\n\n Error: no L1GlobalTriggerReadoutRecord found with input tag "
00384 //       << m_l1GtReadoutRecord
00385 //       << "\n Returning empty L1GlobalTriggerRecord.\n\n"
00386 //       << std::endl;
00387       cout<<" No L1 trigger record "<<endl;
00388    } else {
00389 
00390   const DecisionWord dWord = gtRecord->decisionWord();
00391 
00392   for (CItAlgo itAlgo = bitMap.begin(); itAlgo != bitMap.end(); itAlgo++)
00393     {
00394       bool decision=menu->gtAlgorithmResult(itAlgo->first,dWord);
00395       if(decision == 1) std::cout<<" Trigger "<<itAlgo->first<<" "<<decision<<std::endl;
00396     }
00397 
00398   }
00399 
00400   const HcalRespCorrs* myRecalib=0;
00401   if( theRecalib ) {
00402 // Radek:   
00403   edm::ESHandle <HcalRespCorrs> recalibCorrs;
00404   iSetup.get<HcalRespCorrsRcd>().get("recalibrate",recalibCorrs);
00405   myRecalib = recalibCorrs.product();
00406 // end
00407   } // theRecalib
00408 
00409 // Noise part for HB HE
00410 
00411      double tmpNSFillDetMapPl1[5][5][73][43]; 
00412      double tmpNSFillDetMapPl2[5][5][73][43];
00413      double tmpNSFillDetMapMin1[5][5][73][43]; 
00414      double tmpNSFillDetMapMin2[5][5][73][43];
00415 
00416    for (int i=0; i<5;i++)
00417    {
00418     for (int j=0; j<5;j++)
00419     {
00420      for (int k=0; k<73;k++)
00421      {
00422        for (int l=0; l<43;l++)
00423        {
00424         tmpNSFillDetMapPl1[i][j][k][l] = 0.;
00425         tmpNSFillDetMapPl2[i][j][k][l] = 0.;
00426         tmpNSFillDetMapMin1[i][j][k][l] = 0.;
00427         tmpNSFillDetMapMin2[i][j][k][l] = 0.;
00428        }
00429      }  
00430     }
00431    }    
00432    edm::Handle<HBHERecHitCollection> hbheNormal;
00433    iEvent.getByLabel("hbhereco", hbheNormal);
00434    if(!hbheNormal.isValid()){  
00435      cout<<" hbheNormal failed "<<endl;
00436    } else {
00437        cout<<" The size of the normal collection "<<hbheNormal->size()<<endl;
00438    } 
00439 
00440 
00441    edm::Handle<HBHERecHitCollection> hbheNS;
00442    iEvent.getByLabel(hbherecoNoise, hbheNS);
00443 
00444 
00445    if(!hbheNS.isValid()){
00446      LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
00447      cout<<" No HBHE MS "<<endl;
00448      return ;
00449    }
00450 
00451   
00452    const HBHERecHitCollection HithbheNS = *(hbheNS.product());
00453    cout<<" HBHE NS size of collection "<<HithbheNS.size()<<endl;
00454    hHBHEsize_vs_run->Fill(rnnum,(float)HithbheNS.size());
00455 
00456    if(HithbheNS.size()!= 5184) {
00457           cout<<" HBHE problem "<<rnnum<<" "<<HithbheNS.size()<<endl;
00458           return;
00459    }
00460    edm::Handle<HBHERecHitCollection> hbheMB;
00461    iEvent.getByLabel(hbherecoMB, hbheMB);
00462 
00463    if(!hbheMB.isValid()){
00464      LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
00465      cout<<" No HBHE MB"<<endl;
00466      return ;
00467    }
00468 
00469   const HBHERecHitCollection HithbheMB = *(hbheMB.product());
00470   cout<<" HBHE MB size of collection "<<HithbheMB.size()<<endl;
00471   if(HithbheMB.size()!= 5184) {
00472      cout<<" HBHE problem "<<rnnum<<" "<<HithbheMB.size()<<endl;
00473      return;
00474   }
00475 
00476    edm::Handle<HFRecHitCollection> hfNS;
00477    iEvent.getByLabel(hfrecoNoise, hfNS);
00478 
00479    if(!hfNS.isValid()){
00480      LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
00481      cout<<" No HF NS "<<endl;
00482      return ;
00483    }
00484 
00485    const HFRecHitCollection HithfNS = *(hfNS.product());
00486    cout<<" HFE NS size of collection "<<HithfNS.size()<<endl;
00487    hHFsize_vs_run->Fill(rnnum,(float)HithfNS.size());
00488    if(HithfNS.size()!= 1728) {
00489           cout<<" HF problem "<<rnnum<<" "<<HithfNS.size()<<endl;
00490           return;
00491    }
00492 
00493    edm::Handle<HFRecHitCollection> hfMB;
00494    iEvent.getByLabel(hfrecoMB, hfMB);
00495 
00496    if(!hfMB.isValid()){
00497      LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
00498      cout<<" No HBHE MB"<<endl;
00499      return ;
00500    }
00501 
00502   const HFRecHitCollection HithfMB = *(hfMB.product());
00503   cout<<" HF MB size of collection "<<HithfMB.size()<<endl;
00504    if(HithfMB.size()!= 1728) {
00505           cout<<" HF problem "<<rnnum<<" "<<HithfMB.size()<<endl;
00506           return;
00507    }
00508 
00509 
00510 
00511   for(HBHERecHitCollection::const_iterator hbheItr=HithbheNS.begin(); hbheItr!=HithbheNS.end(); hbheItr++)
00512         {
00513 // Recalibration of energy
00514          float icalconst=1.;     
00515          DetId mydetid = hbheItr->id().rawId();
00516          if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
00517     
00518          HBHERecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
00519          
00520          double energyhit = aHit.energy();
00521          
00522          DetId id = (*hbheItr).detid(); 
00523          HcalDetId hid=HcalDetId(id);
00524  
00525  
00526  
00527          int mysu = ((hid).rawId()>>25)&0x7;
00528          if( hid.ieta() > 0 ) {
00529          theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
00530          theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
00531          theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
00532          theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
00533          
00534          tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
00535          tmpNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
00536          
00537          
00538          } else {
00539          theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
00540          theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
00541          theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
00542          theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
00543 
00544 
00545          tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
00546          tmpNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
00547          
00548          }  
00549          
00550          if(hid.depth() == 1) {
00551            hbheNoiseE->Fill(energyhit);
00552         
00553            if(energyhit<-2.) std::cout<<" Run "<<rnnum<<" ieta,iphi "<<hid.ieta()<<" "<<hid.iphi()<<energyhit<<std::endl;
00554  
00555         // if( hid.ieta() > 0 ) {
00556         //  hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit-noise_pl[hid.iphi()][hid.ieta()]);
00557         //  hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
00558         // } else {
00559         //  hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit-noise_min[hid.iphi()][abs(hid.ieta())]);
00560         //  hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
00561         // } // eta><0
00562         
00563          } // depth=1
00564          
00565          
00566         } // HBHE_NS
00567 
00568 
00569 // Signal part for HB HE
00570      
00571 
00572   for(HBHERecHitCollection::const_iterator hbheItr=HithbheMB.begin(); hbheItr!=HithbheMB.end(); hbheItr++)
00573         {
00574 // Recalibration of energy
00575          float icalconst=1.;     
00576          DetId mydetid = hbheItr->id().rawId();
00577          if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
00578     
00579          HBHERecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
00580          
00581          double energyhit = aHit.energy();
00582          
00583          DetId id = (*hbheItr).detid(); 
00584          HcalDetId hid=HcalDetId(id);
00585  
00586          int mysu = ((hid).rawId()>>25)&0x7;
00587          if( hid.ieta() > 0 ) {
00588          theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
00589          theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
00590          theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
00591          theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
00592          float mydiff = energyhit - tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
00593          
00594          
00595          theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+mydiff;
00596          theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(mydiff,2);
00597          } else {
00598          theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
00599          theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
00600          theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
00601          theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
00602 
00603 
00604          float mydiff = energyhit - tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
00605          theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+mydiff;
00606          theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(mydiff,2);
00607          }  
00608 
00609          
00610          if(hid.depth() == 1) {
00611          
00612          hbheSignalE->Fill(energyhit);
00613          
00614          if( hid.ieta() > 0 ) {
00615           hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit);
00616           hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
00617          } else {
00618           hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit);
00619           hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
00620          } // eta><0
00621          
00622          } // depth=1
00623          
00624          
00625         } // HBHE_MB
00626         
00627 // HF
00628  
00629   for(HFRecHitCollection::const_iterator hbheItr=HithfNS.begin(); hbheItr!=HithfNS.end(); hbheItr++)
00630         {
00631 // Recalibration of energy
00632          float icalconst=1.;     
00633          DetId mydetid = hbheItr->id().rawId();
00634          if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
00635     
00636          HFRecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
00637          
00638          double energyhit = aHit.energy();
00639 //
00640 // Remove PMT hits
00641 //               
00642          if(fabs(energyhit) > 40. ) continue;
00643                  
00644          DetId id = (*hbheItr).detid(); 
00645          HcalDetId hid=HcalDetId(id);
00646  
00647          int mysu = ((hid).rawId()>>25)&0x7;
00648          if( hid.ieta() > 0 ) {
00649          theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
00650          theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
00651          theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
00652          theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theNSFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
00653          
00654          tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
00655          tmpNSFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
00656          
00657          
00658          } else {
00659          theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
00660          theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
00661          theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
00662          theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theNSFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
00663 
00664  
00665          tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = energyhit;
00666          tmpNSFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = pow(energyhit,2);
00667          
00668          }  
00669          
00670          if(hid.depth() == 1) {
00671          hfNoiseE->Fill(energyhit);
00672             
00673          //if( hid.ieta() > 0 ) {
00674          // hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit-noise_pl[hid.iphi()][hid.ieta()]);
00675          // hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
00676          //} else {
00677          // hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit-noise_min[hid.iphi()][abs(hid.ieta())]);
00678          // hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
00679          //} // eta><0
00680          
00681          } // depth=1
00682          
00683         } // HBHE_NS
00684 
00685 
00686 // Signal part for HB HE
00687 
00688   for(HFRecHitCollection::const_iterator hbheItr=HithfMB.begin(); hbheItr!=HithfMB.end(); hbheItr++)
00689         {
00690 // Recalibration of energy
00691          float icalconst=1.;     
00692          DetId mydetid = hbheItr->id().rawId();
00693          if( theRecalib ) icalconst=myRecalib->getValues(mydetid)->getValue();
00694     
00695          HFRecHit aHit(hbheItr->id(),hbheItr->energy()*icalconst,hbheItr->time());
00696          
00697          double energyhit = aHit.energy();
00698 //
00699 // Remove PMT hits
00700 //       
00701          if(fabs(energyhit) > 40. ) continue;
00702          
00703          DetId id = (*hbheItr).detid(); 
00704          HcalDetId hid=HcalDetId(id);
00705  
00706          int mysu = ((hid).rawId()>>25)&0x7;
00707          if( hid.ieta() > 0 ) {
00708          theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl0[mysu][hid.depth()][hid.iphi()][hid.ieta()]+ 1.;
00709          theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit;
00710          theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,2);
00711          theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theMBFillDetMapPl4[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow(energyhit,4);
00712 
00713 
00714          theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit-tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
00715          theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()] =
00716          theDFFillDetMapPl2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow((energyhit-tmpNSFillDetMapPl1[mysu][hid.depth()][hid.iphi()][hid.ieta()]),2);
00717          } else {
00718          theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin0[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+ 1.;
00719          theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin1[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+energyhit;
00720          theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin2[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,2);
00721          theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())] = theMBFillDetMapMin4[mysu][hid.depth()][hid.iphi()][abs(hid.ieta())]+pow(energyhit,4);
00722  
00723          theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()] = theDFFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]+energyhit-tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()];
00724          theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()] =
00725          theDFFillDetMapMin2[mysu][hid.depth()][hid.iphi()][hid.ieta()]+pow((energyhit-tmpNSFillDetMapMin1[mysu][hid.depth()][hid.iphi()][hid.ieta()]),2);
00726          }  
00727 
00728          
00729          if(hid.depth() == 1) {
00730             hfSignalE->Fill(energyhit);
00731             
00732          if( hid.ieta() > 0 ) 
00733          {
00734           hCalo1[hid.iphi()][hid.ieta()]->Fill(energyhit);
00735           hCalo1mom2[hid.iphi()][hid.ieta()]->Fill(pow(energyhit,2));
00736          } else 
00737          {
00738           hCalo2[hid.iphi()][abs(hid.ieta())]->Fill(energyhit);
00739           hCalo2mom2[hid.iphi()][abs(hid.ieta())]->Fill(pow(energyhit,2));
00740          } // eta><0
00741          
00742          } // depth=1
00743          
00744         } // HF_MB
00745 
00746    std::cout<<" Event is finished "<<std::endl;
00747 }
00748 }
00749 //define this as a plug-in
00750 //DEFINE_ANOTHER_FWK_MODULE(Analyzer_minbias)
00751
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/Calibration/HcalCalibAlgos/src/Analyzer_minbias.cc
