LumiProducer/test/standaloneScanLumiRaw.cpp

0001 #include "TFile.h"
0002 #include "TTree.h"
0003 #include "TBranch.h"
0004 #include "RecoLuminosity/LumiProducer/interface/LumiRawDataStructures.h"
0005 #include "RecoLuminosity/LumiProducer/interface/LumiCorrector.h"
0006 #include <iostream>
0007 #include <iomanip>
0008 #include <memory>
0009 #include <vector>
0010 #include <map>
0011 //#include <cmath>
0012 #include <algorithm>
0013 #include <fstream>
0014
0015 /** This programm scans a given lumi raw data file and print out the content
0016  **/
0017 struct beaminfo {
0018   unsigned int bxidx;
0019   unsigned int timestamp;
0020   float lumival;
0021   float beam1_intensity;
0022   float beam2_intensity;
0023 };
0024
0025 int main(int argc, char **argv) {
0026   const char *filename = "file:test.root";
0027   //default file to read. file name is taken from command argument
0028   if (argc > 1) {
0029     filename = argv[1];
0030   }
0031   //TFile *myfile=new TFile(filename,"READ");
0032   TFile *myfile = TFile::Open(filename);
0033
0034   HCAL_HLX::LUMI_SECTION *myLumiSection = new HCAL_HLX::LUMI_SECTION;
0035   HCAL_HLX::LUMI_SECTION_HEADER *myLumiHeader = &(myLumiSection->hdr);
0036   HCAL_HLX::LUMI_SUMMARY *myLumiSummary = &(myLumiSection->lumiSummary);
0037   HCAL_HLX::LUMI_DETAIL *myLumiDetail = &(myLumiSection->lumiDetail);
0038
0039   TTree *hlxTree = (TTree *)myfile->Get("HLXData");
0040   if (!hlxTree)
0041     std::cout << "no hlx data" << std::endl;
0042   hlxTree->SetBranchAddress("Header.", &myLumiHeader);
0043   hlxTree->SetBranchAddress("Summary.", &myLumiSummary);
0044   hlxTree->SetBranchAddress("Detail.", &myLumiDetail);
0045   size_t hlxentries = hlxTree->GetEntries();
0046   //std::cout<<"hlxentries "<<hlxentries<<std::endl;
0047   std::map<unsigned int, std::vector<beaminfo> > bxlumis;
0048   std::vector<beaminfo> tmpbx;
0049   unsigned int ncollidingbx = 0;
0050   for (size_t i = 0; i < hlxentries; ++i) {
0051     hlxTree->GetEntry(i);
0052     ncollidingbx = myLumiHeader->numBunches;
0053     //std::cout<<"Lumi summary for run : "<<myLumiHeader->runNumber<<" : LS : "<<myLumiHeader->sectionNumber<<" "<<myLumiHeader->timestamp<<" "<<myLumiHeader->numBunches<<std::endl;
0054
0055     //std::cout<<std::setw(20)<<"lumi details : "<<std::endl;
0056     unsigned int hlxls = myLumiHeader->sectionNumber;
0057     unsigned int ts = myLumiHeader->timestamp;
0058     tmpbx.clear();
0059     for (size_t j = 0; j < 3564; ++j) {
0060       //std::cout<<std::setw(20)<<"    BX : "<<j<<" : OccLumi : "<<myLumiDetail->OccLumi[0][j]<<std::endl;
0061       beaminfo b;
0062       b.timestamp = ts;
0063       b.bxidx = j;
0064       b.lumival = myLumiDetail->OccLumi[0][j];
0065       b.beam1_intensity = 0.;
0066       b.beam2_intensity = 0.;
0067       tmpbx.push_back(b);
0068     }
0069     bxlumis.insert(std::make_pair(hlxls, tmpbx));
0070   }
0071
0072   LumiCorrector corr;
0073   std::map<unsigned int, std::vector<beaminfo> >::iterator mapIt;
0074   std::map<unsigned int, std::vector<beaminfo> >::iterator itBeg = bxlumis.begin();
0075   std::map<unsigned int, std::vector<beaminfo> >::iterator itEnd = bxlumis.end();
0076
0077   for (mapIt = itBeg; mapIt != itEnd; ++mapIt) {
0078     float totlumi = 0.;
0079     std::vector<beaminfo>::iterator thislslumisBeg = mapIt->second.begin();
0080     std::vector<beaminfo>::iterator thislslumisEnd = mapIt->second.end();
0081     for (std::vector<beaminfo>::iterator it = thislslumisBeg; it != thislslumisEnd; ++it) {
0082       totlumi += it->lumival;
0083     }
0084     for (std::vector<beaminfo>::iterator it = thislslumisBeg; it != thislslumisEnd; ++it) {
0085       float thecorrector = corr.TotalNormOcc1(totlumi * 1.0e-3, ncollidingbx);
0086       float correctedbxlumi = thecorrector * (it->lumival);
0087       it->lumival = correctedbxlumi;
0088     }
0089   }
0090   //lsnum++;
0091   //std::cout<<std::setw(20)<<"#LS "<<lsnum<<" timestamp "<<thislstimestamp<<std::endl;
0092
0093   TTree *diptree = (TTree *)myfile->Get("DIPCombined");
0094
0095   if (diptree) {
0096     std::unique_ptr<HCAL_HLX::DIP_COMBINED_DATA> dipdata(new HCAL_HLX::DIP_COMBINED_DATA);
0097     diptree->SetBranchAddress("DIPCombined.", &dipdata);
0098     size_t ndipentries = diptree->GetEntries();
0099     unsigned int dipls = 0;
0100     if (ndipentries > 0) {
0101       for (size_t i = 0; i < 1; ++i) {
0102         diptree->GetEntry(i);
0103         //unsigned int fillnumber=dipdata->FillNumber;
0104         dipls = dipdata->sectionNumber;
0105         std::map<unsigned int, std::vector<beaminfo> >::iterator dipIt = bxlumis.end();
0106         if (bxlumis.find(dipls) != dipIt) {
0107           dipIt = bxlumis.find(dipls);
0108         }
0109         for (unsigned int i = 0; i < 3564; ++i) {
0110           float beam1in = dipdata->Beam[0].averageBunchIntensities[i];
0111           float beam2in = dipdata->Beam[1].averageBunchIntensities[i];
0112           if (dipIt != bxlumis.end()) {
0113             dipIt->second[i].beam1_intensity = beam1in;
0114             dipIt->second[i].beam2_intensity = beam2in;
0115           }
0116         }
0117       }
0118     }
0119   }
0120
0121   std::ofstream outfile;
0122   outfile.open("out.txt");
0123   for (mapIt = itBeg; mapIt != itEnd; ++mapIt) {
0124     outfile << "# " << mapIt->first << std::endl;
0125     std::vector<beaminfo>::iterator thislslumisBeg = mapIt->second.begin();
0126     std::vector<beaminfo>::iterator thislslumisEnd = mapIt->second.end();
0127     for (std::vector<beaminfo>::iterator it = thislslumisBeg; it != thislslumisEnd; ++it) {
0128       outfile << it->timestamp << "," << it->bxidx << "," << it->lumival << "," << it->beam1_intensity << ","
0129               << it->beam2_intensity << std::endl;
0130     }
0131   }
0132   outfile.close();
0133 }
0134 //TFile * myfile=TFile::Open("rfio:/castor/cern.ch/cms/store/lumi/200912/CMS_LUMI_RAW_20091212_000124025_0001_1.root");
0135 //HCAL_HLX::RUN_SUMMARY *myRunSummary = new HCAL_HLX::RUN_SUMMARY;
0136 //TTree *runsummaryTree = (TTree *) myfile->Get("RunSummary");
0137 //if(!runsummaryTree) std::cout<<"no run summary data"<<std::endl;
0138 //runsummaryTree->SetBranchAddress("RunSummary.",&myRunSummary);
0139 //size_t runsummaryentries=runsummaryTree->GetEntries();
0140 //std::cout<<"n run summary entries "<<runsummaryentries<<std::endl;
0141 //for(size_t i=0;i<runsummaryentries;++i){
0142 //runsummaryTree->GetEntry(i);
0143 // std::cout<<"Summary for run : "<<myRunSummary->runNumber<<std::endl;
0144 // std::cout<<std::setw(20)<<"timestamp : "<<myRunSummary->timestamp<<" : timestamp micros : "<<myRunSummary->timestamp_micros<<" : start orbit : "<<myRunSummary->startOrbitNumber<<" : end orbit : "<<myRunSummary->endOrbitnumber<<" : fill number : "<<myRunSummary->fillNumber<<" : number CMS LS : "<<myRunSummary->numberCMSLumiSections<<" : number DAQ LS : "<<myRunSummary->numberLumiDAQLumiSections<<std::endl;
0145 //}
0146 /**
0147   HCAL_HLX::LEVEL1_TRIGGER *myTRG = new HCAL_HLX::LEVEL1_TRIGGER;
0148   TTree *trgTree = (TTree *) myfile->Get("L1Trigger");
0149   if(!trgTree) std::cout<<"no trg data"<<std::endl;
0150   trgTree->SetBranchAddress("L1Trigger.",&myTRG);
0151   size_t trgentries=trgTree->GetEntries();
0152   for(size_t i=0;i<trgentries;++i){
0153     trgTree->GetEntry(i);
0154     //std::cout<<"trg runnumber "<<myTRG->runNumber<<std::endl;
0155     std::cout<<"TRG for run : "<< myTRG->runNumber<<" : LS : "<<myTRG->sectionNumber<<" : deadtime : "<< myTRG->deadtiecount<<std::endl;
0156     for( unsigned int j=0; j<128; ++j){
0157       std::cout<<std::setw(20)<<"GT Algo  "<<j;
0158       std::cout<<" : path : "<< myTRG->GTAlgo[j].pathName<<" : counts : "<< myTRG->GTAlgo[j].counts<<" : prescale : "<< myTRG->GTAlgo[j].prescale<<std::endl;
0159     }
0160     for( unsigned int k=0; k<64; ++k){
0161       std::cout<<std::setw(20)<<"GT Tech : "<<k;
0162       std::cout<<" : path : "<< myTRG->GTTech[k].pathName<<" : counts : "<< myTRG->GTTech[k].counts<<" : prescale : "<< myTRG->GTTech[k].prescale<<std::endl;
0163     }
0164   }
0165   **/
0166 /**
0167   HCAL_HLX::HLTRIGGER *myHLT = new HCAL_HLX::HLTRIGGER;
0168   TTree *hltTree = (TTree *) myfile->Get("HLTrigger");
0169   if(!hltTree) std::cout<<"no hlt data"<<std::endl;
0170   hltTree->SetBranchAddress("HLTrigger.",&myHLT);
0171   size_t hltentries=hltTree->GetEntries();
0172   for(size_t i=0;i<hltentries;++i){
0173     hltTree->GetEntry(i);
0174     std::cout<<"HLT for run : "<< myHLT->runNumber<<":  LS  : "<<myHLT->sectionNumber<<" : total hlt paths : "<<myHLT->numPaths<<std::endl;
0175     for( unsigned int j=0; j<myHLT->numPaths; ++j){
0176       std::cout<<std::setw(20)<<"HLTConfigId : "<<myHLT->HLTPaths[j].HLTConfigId<<"path : "<<myHLT->HLTPaths[j].PathName<<" : L1Pass : "<<myHLT->HLTPaths[j].L1Pass<<" : PSPass : "<<myHLT->HLTPaths[j].PSPass<<" : PAccept : "<<myHLT->HLTPaths[j].PAccept<<" : PExcept : "<<myHLT->HLTPaths[j].PExcept<<" : PReject : "<<myHLT->HLTPaths[j].PReject<<" : PSIndex : "<<myHLT->HLTPaths[j].PSIndex<<" : Prescale : "<<myHLT->HLTPaths[j].Prescale<<std::endl;
0177     }
0178   }
0179   **/
0180 /**
0181   HCAL_HLX::LUMI_SECTION *myLumiSection=new HCAL_HLX::LUMI_SECTION;
0182   HCAL_HLX::LUMI_SECTION_HEADER *myLumiHeader = &(myLumiSection->hdr);
0183   HCAL_HLX::LUMI_SUMMARY *myLumiSummary = &(myLumiSection->lumiSummary);
0184   HCAL_HLX::LUMI_DETAIL *myLumiDetail = &(myLumiSection->lumiDetail);
0185
0186   TTree *hlxTree = (TTree *) myfile->Get("HLXData");
0187   if(!hlxTree) std::cout<<"no hlx data"<<std::endl;
0188   hlxTree->SetBranchAddress("Header.",&myLumiHeader);
0189   hlxTree->SetBranchAddress("Summary.",&myLumiSummary);
0190   hlxTree->SetBranchAddress("Detail.",&myLumiDetail);
0191   size_t hlxentries=hlxTree->GetEntries();
0192   std::cout<<"hlxentries "<<hlxentries<<std::endl;
0193   for(size_t i=0;i<hlxentries;++i){
0194     hlxTree->GetEntry(i);
0195     std::cout<<"Lumi summary for run : "<<myLumiHeader->runNumber<<" : LS : "<<myLumiHeader->sectionNumber<<" cmsalive value: "<<myLumiHeader->bCMSLive<<std::endl;
0196     bool a=true;
0197     if(typeid(myLumiHeader->bCMSLive)==typeid(a)){
0198       std::cout<<"is bool type"<<std::endl;
0199       std::cout<<"normal bool "<<a<<std::endl;
0200       std::cout<<"cms alive bool "<< myLumiHeader->bCMSLive<< std::endl;
0201     }else{
0202       std::cout<<"not bool type"<<std::endl;
0203       std::cout<<"cms alive"<< myLumiHeader->bCMSLive<< std::endl;
0204     }
0205     //std::cout<<std::setw(20)<<"deadtime norm : "<<myLumiSummary->DeadTimeNormalization<<" : LHC norm : "<<myLumiSummary->LHCNormalization<<" : instantlumi : "<<myLumiSummary->InstantLumi<<" : instantlumiErr : "<<myLumiSummary->InstantLumiErr<<" : instantlumiQlty : "<<myLumiSummary->InstantLumiQlty<<std::endl;
0206     //std::cout<<std::setw(20)<<"lumi details : "<<std::endl;
0207     //for(size_t j=0;j<HCAL_HLX_MAX_BUNCHES;++j){
0208     //  std::cout<<std::setw(20)<<"    LHCLumi : "<<myLumiDetail->LHCLumi[j]<<" : ETLumi : "<<myLumiDetail->ETLumi[j]<<" : ETLumiErr : "<<myLumiDetail->ETLumiErr[j]<<" : ETLumiQlty : "<<myLumiDetail->ETLumiQlty[j]<<" : ETBXNormalization : "<<myLumiDetail->ETBXNormalization[j]<<std::endl;
0209     //}
0210   }
0211   **/
0212 //}