From 5845466564a3b3fa7e59cb0abaee13b5949753c7 Mon Sep 17 00:00:00 2001 From: Shervin Nourbakhsh Date: Mon, 17 Oct 2016 11:29:54 +0200 Subject: [PATCH 1/7] Removed nspills and updated documentation of TextSource --- RawToDigi/plugins/HGCalTBTextSource.cc | 1 - RawToDigi/plugins/HGCalTBTextSource.h | 64 +++++++++++++++++++++++--- doc/instructions.doc | 3 +- doc/main_page.doc | 6 ++- fulldoc | 4 +- test_cfg.py | 10 +--- 6 files changed, 69 insertions(+), 19 deletions(-) diff --git a/RawToDigi/plugins/HGCalTBTextSource.cc b/RawToDigi/plugins/HGCalTBTextSource.cc index bd5b395..5237c9a 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.cc +++ b/RawToDigi/plugins/HGCalTBTextSource.cc @@ -40,7 +40,6 @@ bool HGCalTBTextSource::setRunAndEventInfo(edm::EventID& id, edm::TimeValue_t& t id = edm::EventID(m_run, m_spill, m_event); time = (edm::TimeValue_t) m_time; - if(m_spill <= NSpills) return true; else return false; } diff --git a/RawToDigi/plugins/HGCalTBTextSource.h b/RawToDigi/plugins/HGCalTBTextSource.h index afc8fda..b97ff7a 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.h +++ b/RawToDigi/plugins/HGCalTBTextSource.h @@ -8,13 +8,68 @@ #include #include #include + +/** + \page TXTFORMAT_PAGE Data txt input format: + \tableofcontents + Return to the main page: \ref index + \section TXTFORMAT Original data txt format: + - Spill header containing spill time, run number and event number in the spill + \verbatim + STARTING SPILL READ AT TIME (1us): 0x55DC7602 RUN: 802 EVENT: 150 + \endverbatim + - Board header with the FMC-IO identification number + \verbatim + Board header: on FMC-IO 11, trig_count in mem= 150, sk_status = 1 + \endverbatim + - Event header with event number, timestamp of the event and the global trigger time + \verbatim + Event header for event 0 with (200ns) timestamp 0xB803BA11 global tts(us) 0x00000000 and CKOV= 0 + \endverbatim + - 68 lines with: eventNumber, channelID, 32bit word formed by two 16bit words for low gain values of two channels, 32bit word formed by two 16bit words for high gain vaues of two channels + the 68 lines have then the information from two skirocs + \verbatim + 0 0 0x11B411B8 0x11941188 + \endverbatim + + The original txt format is not suitable for being processed by CMSSW, events are not ordered. For each board all the events are dumped. + + \section NEWTXTFORMAT_ Rearranged data txt format: + In order to be able to process the data with CMSSW, we need to rearrange the txt in order to have for each event all the boards. + + This is done with the script \verbatim ./scripts/rearrangeTxtFile.sh \endverbatim + + This step is done centrally and files are available in ... + + The txt file is processed by HGCalTBTextSource + + \section EXAMPLE Usage example in your python cfg: + \code + process.source = cms.Source("HGCalTBTextSource", + fileNames=cms.untracked.vstring("file:myfile1.txt"), + ) + \endcode + + Return to the main page: \ref index + + \example test_cfg.py + +*/ + /** * \class HGCalTBTextSource HGCal/RawToDigi/plugins/HGCalTBTextSource.h * * \brief convert data from txt file to FEDRawData * - * \todo replace c-like scanf with c++ versions - * \todo change run and fed IDs (now are hardcoded) + * \details + * For info about the txt input format see: \ref TXTFORMAT_PAGE + * \author Shervin Nourbakhsh + * \author Rajdeep Mohan Chatterjee + * \author Jeremy Mans + * + * \todo fedID is taken from cfg, it should be taken from the txt file + * + * \example test_cfg.py source module */ class HGCalTBTextSource : public edm::ProducerSourceFromFiles { @@ -22,10 +77,8 @@ class HGCalTBTextSource : public edm::ProducerSourceFromFiles public: explicit HGCalTBTextSource(const edm::ParameterSet & pset, edm::InputSourceDescription const& desc) : edm::ProducerSourceFromFiles(pset, desc, true), m_file(0), - NSpills(pset.getUntrackedParameter("nSpills", 6)) { - - m_sourceId = pset.getUntrackedParameter("fed", 1000); /// \todo check and read from file? + m_sourceId = pset.getUntrackedParameter("fed", 1000); /// \todo read from file produces(); } @@ -47,6 +100,5 @@ class HGCalTBTextSource : public edm::ProducerSourceFromFiles FILE* m_file; unsigned int m_time; unsigned int m_event, m_run, m_spill, max_boards; - unsigned int NSpills;//Read while running how many spills we wish to run over int m_sourceId; }; diff --git a/doc/instructions.doc b/doc/instructions.doc index d3ffbe5..0329873 100644 --- a/doc/instructions.doc +++ b/doc/instructions.doc @@ -1,5 +1,6 @@ /** -\page instructions_ Instructions +\page INSTRUCTIONS_PAGE Download instructions + \par Prerequisites for using GIT with github Please make sure you have your github account, you have properly added your private key such that you are able to use the SSH protocol diff --git a/doc/main_page.doc b/doc/main_page.doc index 6a8b185..1171a1c 100644 --- a/doc/main_page.doc +++ b/doc/main_page.doc @@ -2,7 +2,11 @@ \mainpage \tableofcontents -Download instructions: \ref instructions_ +\section INSTRUCTIONS_ Download instructions +\subpage INSTRUCTIONS_PAGE + +\section TXTFORMAT_ Txt file format +\subpage TXTFORMAT_PAGE \section INPUTS_ Input files: - channel mapping: file that contains the conversion between electronicID (SKIROC chip number) to detector ID (DetId) diff --git a/fulldoc b/fulldoc index e0a5e12..0d3bc30 100644 --- a/fulldoc +++ b/fulldoc @@ -396,7 +396,7 @@ LOOKUP_CACHE_SIZE = 0 # normally produced when WARNINGS is set to YES. # The default value is: NO. -EXTRACT_ALL = YES +EXTRACT_ALL = NO # If the EXTRACT_PRIVATE tag is set to YES all private members of a class will # be included in the documentation. @@ -801,7 +801,7 @@ EXCLUDE_SYMBOLS = # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = +EXAMPLE_PATH = test_cfg.py # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and diff --git a/test_cfg.py b/test_cfg.py index 4f17fa3..2cd155e 100644 --- a/test_cfg.py +++ b/test_cfg.py @@ -41,11 +41,6 @@ VarParsing.VarParsing.varType.int, '0: if runType is PED then do Digi, if runType is HGC_Run then do Digi and Reco (not implemented yet); 1: do Digi; 2: only Reco (not implemented yet); 3: Digi + highgain_correlation_cm; 4: event display sequence; 5: highgain_correlation_cm + event display sequence') -options.register('nSpills', - 15, - VarParsing.VarParsing.multiplicity.singleton, - VarParsing.VarParsing.varType.int, - 'Number of spills in run') options.register('pedestalsHighGain', 'CondObjects/data/Ped_HighGain_L8.txt', @@ -91,14 +86,13 @@ process.load('HGCal.StandardSequences.LocalReco_cff') process.load('HGCal.StandardSequences.dqm_cff') - +## [source module] process.source = cms.Source("HGCalTBTextSource", run=cms.untracked.int32(options.runNumber), ### maybe this should be read from the file #fileNames=cms.untracked.vstring("file:Raw_data_New.txt") ### here a vector is provided, but in the .cc only the first one is used TO BE FIXED fileNames=cms.untracked.vstring("file:%s/%s_Output_%06d.txt"%(options.dataFolder,options.runType,options.runNumber)), ### here a vector is provided, but in the .cc only the first one is used TO BE FIXED - nSpills=cms.untracked.uint32(options.nSpills), ) - +## [source module] ###### From 45ee5d26829c4f9285940b558f2a1df211038316 Mon Sep 17 00:00:00 2001 From: Shervin Nourbakhsh Date: Mon, 17 Oct 2016 12:21:21 +0200 Subject: [PATCH 2/7] Cleaning and improvements in code clarity --- DataFormats/interface/SKIROCParameters.h | 5 +- .../interface/HGCalTBGeometryParameters.h | 6 +- RawToDigi/plugins/HGCalTBRawToDigi.cc | 61 +++++++------------ RawToDigi/plugins/HGCalTBRawToDigi.h | 3 + RawToDigi/plugins/HGCalTBTextSource.cc | 42 +++++++------ RawToDigi/plugins/HGCalTBTextSource.h | 25 +++++--- fulldoc | 2 +- 7 files changed, 70 insertions(+), 74 deletions(-) diff --git a/DataFormats/interface/SKIROCParameters.h b/DataFormats/interface/SKIROCParameters.h index 3d5e2a9..8d70fe7 100644 --- a/DataFormats/interface/SKIROCParameters.h +++ b/DataFormats/interface/SKIROCParameters.h @@ -3,9 +3,8 @@ namespace SKIROC { -static const unsigned char NCHANNELS = 64; ///< number of channels read by one SKIROC -static const unsigned char MAXSAMPLES = 2; ///< one sample for high gain and one for low gain // 15 is the real max of the skiroc -static const unsigned int NLAYERS = 1; ///< number of layers \todo this should not be hard coded here + static const unsigned char NCHANNELS = 64; ///< number of channels read by one SKIROC + static const unsigned char MAXSAMPLES = 2; ///< one sample for high gain and one for low gain // 15 is the real max of the skiroc } diff --git a/Geometry/interface/HGCalTBGeometryParameters.h b/Geometry/interface/HGCalTBGeometryParameters.h index 45c0db9..ca975fe 100644 --- a/Geometry/interface/HGCalTBGeometryParameters.h +++ b/Geometry/interface/HGCalTBGeometryParameters.h @@ -1,4 +1,6 @@ -/* now we have only 2 layers */ -#define MAXSKIROCS 16 /* MAXLAYERS = MAXSKIROCS/2 */ #define MAXLAYERS 8 +/* */ +#define MAXSKIROCS_PER_BOARD 2 +/* now we have only 2 layers */ +#define MAXSKIROCS MAXLAYERS*MAXSKIROCS_PER_BOARD diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.cc b/RawToDigi/plugins/HGCalTBRawToDigi.cc index c995dfb..6afe8f7 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.cc +++ b/RawToDigi/plugins/HGCalTBRawToDigi.cc @@ -5,6 +5,7 @@ using namespace std; unsigned int gray_to_binary (unsigned int gray); + HGCalTBRawToDigi::HGCalTBRawToDigi(edm::ParameterSet const& conf): dataTag_(conf.getParameter("InputLabel")), fedId_(conf.getUntrackedParameter("fedId")), @@ -36,52 +37,32 @@ void HGCalTBRawToDigi::produce(edm::Event& e, const edm::EventSetup& c) const FEDRawData& fed = rawraw->FEDData(fedId_); if (fed.size() != 0) { /// \todo Exception if 0???? // we can figure out the number of samples from the size of the raw data - int nsamples = fed.size() / (sizeof(uint16_t) * SKIROC::NCHANNELS * 2); // 2 is for ADC and TDC - digis = std::auto_ptr(new SKIROC2DigiCollection(nsamples)); + size_t nSkirocs = fed.size() / (sizeof(uint16_t) * SKIROC::NCHANNELS * SKIROC::MAXSAMPLES); + size_t nBoards = nSkirocs / MAXSKIROCS_PER_BOARD; + digis = std::auto_ptr(new SKIROC2DigiCollection(nSkirocs * SKIROC::NCHANNELS * SKIROC::MAXSAMPLES)); const uint16_t* pdata = (const uint16_t*)(fed.data()); // we start from the back... - int ptr = fed.size() / sizeof(uint16_t) - 1; - /* - printf("Starting on SKIROC %x\n", pdata[ptr]); - ptr--; // now we are pointing at a relatively-useless header word - ptr--; // now we are pointing at the first TDC word - */ - - for (int ichan = 0; ichan < SKIROC::NCHANNELS; ichan++) { - for (int ski = MAXSKIROCS; ski >= 0; ski--) { - HGCalTBElectronicsId eid(ski, ichan); - if (!essource_.emap_.existsEId(eid.rawId())) { -// std::cout << "We do not have a mapping for " << eid; - } else { - HGCalTBDetId did = essource_.emap_.eid2detId(eid); - digis->addDataFrame(did); - if((ski % 2) == 0){ - if((ichan % 2) == 0){ - ptradc1 = ptr - ichan*2 - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)); - ptradc2 = ptr - ichan*2 - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)) - (2*SKIROC::NCHANNELS); - } - else{ - ptradc1 = ptr - (ichan*2 - 1) - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)); - ptradc2 = ptr - (ichan*2 - 1) - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)) - (2*SKIROC::NCHANNELS); - } + size_t ptr = fed.size() / sizeof(uint16_t) - 1; + size_t ski = 0; + for (unsigned int i_board = 0 ; i_board < nBoards; ++i_board) { + for(size_t i_skiroc = 0; i_skiroc < MAXSKIROCS_PER_BOARD; ++i_skiroc) { + for (int ichan = 0; ichan < SKIROC::NCHANNELS; ichan++) { + HGCalTBElectronicsId eid(ski, ichan); + if (essource_.emap_.existsEId(eid.rawId())) { + HGCalTBDetId did = essource_.emap_.eid2detId(eid); + digis->addDataFrame(did); +#ifdef DEBUG + if(i_board == 0 && i_skiroc == 0) std::cout << (pdata & 0xFFF) << (pdata++ & 0xFFF) << (pdata & 0xFFF) << std::endl; +#endif + digis->backDataFrame().setSample(0, gray_to_binary(pdata++ & 0xFFF), gray_to_binary( pdata++ & 0xFFF), 0); + } - else{ - if((ichan % 2) == 0){ - ptradc1 = ptr - (ichan*2 + 2) - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)); - ptradc2 = ptr - (ichan*2 + 2) - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)) - (2*SKIROC::NCHANNELS); - } - else{ - ptradc1 = ptr - (ichan*2 + 1) - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)); - ptradc2 = ptr - (ichan*2 + 1) - (((MAXSKIROCS - ski)/2)*(4*SKIROC::NCHANNELS)) - (2*SKIROC::NCHANNELS); - } - } - digis->backDataFrame().setSample(0, gray_to_binary(pdata[ptradc1] & 0xFFF), gray_to_binary( pdata[ptradc2] & 0xFFF), 0); -// cout< 0 diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.h b/RawToDigi/plugins/HGCalTBRawToDigi.h index bd7ab05..bb2b65a 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.h +++ b/RawToDigi/plugins/HGCalTBRawToDigi.h @@ -19,6 +19,9 @@ /** * \class HGCal/RawToDigi/plugins/HGCalTBRawToDigi.h HGCalTBRawToDigi.h HGCalTBRawToDigi * \brief Produces a digi collection starting from FEDRawData + * \author Shervin Nourbakhsh (UMN) + * \author Rajdeep Mohan Chatterjee (UMN) + * \author Jeremy Mans (UMN) */ diff --git a/RawToDigi/plugins/HGCalTBTextSource.cc b/RawToDigi/plugins/HGCalTBTextSource.cc index 5237c9a..a889134 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.cc +++ b/RawToDigi/plugins/HGCalTBTextSource.cc @@ -6,6 +6,8 @@ #include "HGCal/Geometry/interface/HGCalTBGeometryParameters.h" #include "HGCal/Geometry/interface/HGCalTBSpillParameters.h" using namespace std; +#define TXTLINES 69 + //#define DEBUG /** RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=09 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=12 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=14 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=17 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=18 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=19 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=22 RUN=000880 SPILL=01 EVENT=000000 GLOBALTIME=0x01B6EAF2 BOARD=23 @@ -76,12 +78,12 @@ bool HGCalTBTextSource::readLines() board.clear(); } - char buff[1024], buff_SK0[1024], buff_SK1[1024]; + char buff[1024], buff_SK[1024]; buff[0] = 0; - unsigned int data_sk0, data_sk1; + unsigned int data_sk; // loop over all the lines - for(unsigned int i = 0; i < 69 && !feof(m_file); ++i) { + for(unsigned int i = 0; i < TXTLINES && !feof(m_file); ++i) { buff[0] = 0; fgets(buff, 1000, m_file); // loop over one line of the text file @@ -89,25 +91,26 @@ bool HGCalTBTextSource::readLines() std::istringstream buffer(b); unsigned int board_counter = 0; if((i < 1) || (i > 64)) continue;// Only these are data words, the rest dont interest us for now - while( buffer.peek() != EOF && board_counter < MAXLAYERS) { //buffer.good() gives compilation errors - // read the data of the two skirocs of one board - buffer >> buff_SK1; - buffer >> buff_SK0; - data_sk0 = strtoul(buff_SK0,NULL,0); - data_sk1 = strtoul(buff_SK1,NULL,0); + while( buffer.peek() != EOF && buffer.good()) { //continue for all the boards - // extra security - m_lines[board_counter].push_back(data_sk0); - m_lines[board_counter].push_back(data_sk1); - + for(size_t iSkiroc = 0; iSkiroc < MAXSKIROCS_PER_BOARD; ++iSkiroc) { + // read the data of the one skiroc of one board + buffer >> buff_SK; + data_sk = strtoul(buff_SK, NULL, 0); + m_lines[board_counter][iSkiroc].push_back(data_sk); + //m_lines[board_counter].push_back(data_sk1); + } ++board_counter; if(board_counter > max_boards) max_boards = board_counter; } } -// if(sscanf(buffer, "Event header for event %x with (200ns) timestamp %x", &tmp_event, &m_time) == 2) { return !m_lines.empty(); } +/** DATA in the FEDRAWDATA in the following order: + * BOARD, SKIROC, CHANNEL + * so [Board0 [Skiroc0 [Channel0][Channel1][Channel2]...][Skiroc1 [Channel0][Channel1][Channel2]...]] + */ void HGCalTBTextSource::produce(edm::Event & e) { std::auto_ptr bare_product(new FEDRawDataCollection()); @@ -116,10 +119,12 @@ void HGCalTBTextSource::produce(edm::Event & e) // make sure there are an even number of 32-bit-words (a round number of 64 bit words... for (unsigned int i_board = 0 ; i_board < max_boards; ++i_board) { - auto board = m_lines[i_board]; - for (auto skiword : board) { - skiwords.push_back(skiword >> 16); - skiwords.push_back(skiword); + for(size_t i_skiroc = 0; i_skiroc < MAXSKIROCS_PER_BOARD; ++i_skiroc) { + auto board = m_lines[i_board][i_skiroc]; + for (auto skiword : board) { + skiwords.push_back(skiword >> 16); + skiwords.push_back(skiword); + } } } @@ -138,7 +143,6 @@ void HGCalTBTextSource::fillDescriptions(edm::ConfigurationDescriptions& descrip desc.setComment("TEST"); desc.addUntracked("run", 101); desc.addUntracked >("fileNames"); - desc.addUntracked("nSpills", 6); descriptions.add("source", desc); } diff --git a/RawToDigi/plugins/HGCalTBTextSource.h b/RawToDigi/plugins/HGCalTBTextSource.h index b97ff7a..8c73055 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.h +++ b/RawToDigi/plugins/HGCalTBTextSource.h @@ -46,30 +46,32 @@ \section EXAMPLE Usage example in your python cfg: \code process.source = cms.Source("HGCalTBTextSource", - fileNames=cms.untracked.vstring("file:myfile1.txt"), + fileNames=cms.untracked.vstring("file:myfile1.txt"), ) \endcode Return to the main page: \ref index - \example test_cfg.py + \example test_cfg.py */ /** * \class HGCalTBTextSource HGCal/RawToDigi/plugins/HGCalTBTextSource.h * - * \brief convert data from txt file to FEDRawData + * \brief Convert data from txt file to FEDRawData * * \details * For info about the txt input format see: \ref TXTFORMAT_PAGE - * \author Shervin Nourbakhsh - * \author Rajdeep Mohan Chatterjee - * \author Jeremy Mans + * \author Shervin Nourbakhsh (UMN) + * \author Rajdeep Mohan Chatterjee (UMN) + * \author Jeremy Mans (UMN) * * \todo fedID is taken from cfg, it should be taken from the txt file * - * \example test_cfg.py source module + * DATA in the FEDRAWDATA in the following order: + * BOARD, SKIROC, CHANNEL + * so [Board0 [Skiroc0 [Channel0][Channel1][Channel2]...][Skiroc1 [Channel0][Channel1][Channel2]...]] */ class HGCalTBTextSource : public edm::ProducerSourceFromFiles { @@ -92,11 +94,16 @@ class HGCalTBTextSource : public edm::ProducerSourceFromFiles private: bool setRunAndEventInfo(edm::EventID& id, edm::TimeValue_t& time, edm::EventAuxiliary::ExperimentType&); - virtual void produce(edm::Event & e); + + /** DATA in the FEDRAWDATA in the following order: + * BOARD, SKIROC, CHANNEL + * so [Board0 [Skiroc0 [Channel0][Channel1][Channel2]...][Skiroc1 [Channel0][Channel1][Channel2]...]] + */ + virtual void produce(edm::Event & e); ///test bool readHeader(void); bool readLines(void); - std::array< std::vector < unsigned int> , MAXLAYERS> m_lines; + std::array< std::array< std::vector < unsigned int> , MAXLAYERS>, MAXSKIROCS_PER_BOARD> m_lines; FILE* m_file; unsigned int m_time; unsigned int m_event, m_run, m_spill, max_boards; diff --git a/fulldoc b/fulldoc index 0d3bc30..cce604f 100644 --- a/fulldoc +++ b/fulldoc @@ -447,7 +447,7 @@ EXTRACT_ANON_NSPACES = NO # section is generated. This option has no effect if EXTRACT_ALL is enabled. # The default value is: NO. -HIDE_UNDOC_MEMBERS = NO +HIDE_UNDOC_MEMBERS = YES # If the HIDE_UNDOC_CLASSES tag is set to YES, doxygen will hide all # undocumented classes that are normally visible in the class hierarchy. If set From c4ab7e8ceb99bc97da2e8b572f1e949a7bb37eba Mon Sep 17 00:00:00 2001 From: Shervin Nourbakhsh Date: Mon, 17 Oct 2016 14:24:31 +0200 Subject: [PATCH 3/7] Update with fixes and cleaning of DigiPlotter --- Geometry/interface/HGCalTBCellParameters.h | 2 + Geometry/interface/HGCalTBCellVertices.h | 11 ++ Geometry/interface/HGCalTBSpillParameters.h | 1 - Geometry/interface/HGCalTBTopology.h | 9 +- RawToDigi/plugins/DigiPlotter.cc | 195 +++++++++----------- RawToDigi/plugins/HGCalTBRawToDigi.cc | 25 ++- RawToDigi/plugins/HGCalTBRawToDigi.h | 1 - RawToDigi/plugins/HGCalTBTextSource.cc | 23 +-- RawToDigi/plugins/HGCalTBTextSource.h | 4 +- Reco/plugins/Layer_Sum_Analyzer.cc | 1 - Reco/plugins/RecHitPlotter_HighGain_New.cc | 3 +- test_cfg.py | 3 +- 12 files changed, 139 insertions(+), 139 deletions(-) delete mode 100644 Geometry/interface/HGCalTBSpillParameters.h diff --git a/Geometry/interface/HGCalTBCellParameters.h b/Geometry/interface/HGCalTBCellParameters.h index 494f5a1..9cc06b2 100644 --- a/Geometry/interface/HGCalTBCellParameters.h +++ b/Geometry/interface/HGCalTBCellParameters.h @@ -6,5 +6,7 @@ namespace HGCAL_TB_CELL static const double FULL_CELL_SIDE = 0.6496345; // in cm static const double CALIB_PAD_SIDE = 0.24; // in cm static const double MOUSE_BITE_SIDE = 0.39; // in cm +static const unsigned int HalfHexVertices = 4; +static const unsigned int FullHexVertices = 6; } #endif diff --git a/Geometry/interface/HGCalTBCellVertices.h b/Geometry/interface/HGCalTBCellVertices.h index 5c5cb2a..aa0934c 100644 --- a/Geometry/interface/HGCalTBCellVertices.h +++ b/Geometry/interface/HGCalTBCellVertices.h @@ -25,6 +25,12 @@ class HGCalTBCellVertices HGCalTBCellVertices(); ///< Constructor from cell \b iu & \b iv, valid sensorSizes are 128 and 256 std::vector> GetCellCoordinates(int layer, int sensor_iu, int sensor_iv, int iu, int iv, int sensorsize, bool flipX = false); ///< returns the coordinates of each vertex of cell in the lab frame \b (x,y) + + + inline std::vector> GetCellCoordinatesForPlots(const HGCalTBDetId& detId, int sensorSize) + { + return GetCellCoordinatesForPlots(detId.layer(), detId.sensorIU(), detId.sensorIV(), detId.iu(), detId.iv(), sensorSize); + } inline std::vector> GetCellCoordinatesForPlots(int layer, int sensor_iu, int sensor_iv, int iu, int iv, int sensorsize) { @@ -33,6 +39,10 @@ class HGCalTBCellVertices std::pair GetCellCentreCoordinates(int layer, int sensor_iu, int sensor_iv, int iu, int iv, int sensorsize, bool flipX = false); ///< returns the center of the cell in absolute coordinates: \b (x,y) + inline std::pair GetCellCentreCoordinatesForPlots(const HGCalTBDetId& detId, int sensorSize) + { + return GetCellCentreCoordinatesForPlots(detId.layer(), detId.sensorIU(), detId.sensorIV(), detId.iu(), detId.iv(), sensorSize); + } inline std::pair GetCellCentreCoordinatesForPlots(int layer, int sensor_iu, int sensor_iv, int iu, int iv, int sensorsize) { return GetCellCentreCoordinates(layer, sensor_iu, sensor_iv, iu, iv, sensorsize, true); @@ -41,6 +51,7 @@ class HGCalTBCellVertices // void CellType(int iu, int v, bool ValidFlag);// 1 for full hex, 2 for half hex and 3 for the pentagons(to be implemented later) private: + double a = HGCAL_TB_CELL::FULL_CELL_SIDE; // Size in terms of 1 unit of x/y co-ordinates of a cell side which is 0.064 cm double x_a = sqrt(3) / 2; // cosine pi/6 double y_a = 1 / 2.; // sine pi/6 diff --git a/Geometry/interface/HGCalTBSpillParameters.h b/Geometry/interface/HGCalTBSpillParameters.h deleted file mode 100644 index 07dfd03..0000000 --- a/Geometry/interface/HGCalTBSpillParameters.h +++ /dev/null @@ -1 +0,0 @@ -#define EVENTSPERSPILL 400 diff --git a/Geometry/interface/HGCalTBTopology.h b/Geometry/interface/HGCalTBTopology.h index d085f40..8742483 100644 --- a/Geometry/interface/HGCalTBTopology.h +++ b/Geometry/interface/HGCalTBTopology.h @@ -1,6 +1,6 @@ #ifndef HGCAL_GEOMETRY_HGCALTBTOPOLOGY_H #define HGCAL_GEOMETRY_HGCALTBTOPOLOGY_H 1 - +#include "HGCal/DataFormats/interface/HGCalTBDetId.h" /** \class HGCalTBTopology * * Reference: https://indico.cern.ch/event/456955/ @@ -14,7 +14,14 @@ class HGCalTBTopology public: // valid sensorSizes are 128 and 256 bool iu_iv_valid(int layer, int sensor_iu, int sensor_iv, int iu, int iv, int sensorSize) const; + + inline bool iu_iv_valid(const HGCalTBDetId& detId, int sensorSize) const + { + return iu_iv_valid(detId.layer(), detId.sensorIU(), detId.sensorIV(), detId.iu(), detId.iv(), sensorSize); + }; + double Cell_Area(int cell_type) const;//returns area in cm*cm + }; #endif diff --git a/RawToDigi/plugins/DigiPlotter.cc b/RawToDigi/plugins/DigiPlotter.cc index 17eb028..f6fcfb2 100644 --- a/RawToDigi/plugins/DigiPlotter.cc +++ b/RawToDigi/plugins/DigiPlotter.cc @@ -5,14 +5,16 @@ // /**\class DigiPlotter DigiPlotter.cc HGCal/DigiPlotter/plugins/DigiPlotter.cc - Description: [one line class summary] - - Implementation: - [Notes on implementation] + Description: Plugin to make 2D and 1D histograms of digis + + Implementation: + \author Rajdeep Mohan Chatterjee + \author Shervin Nourbakhsh */ // // Original Author: Rajdeep Mohan Chatterjee // Created: Mon, 15 Feb 2016 09:47:43 GMT +// Modified by Shervin Nourbakhsh // // @@ -31,10 +33,9 @@ #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ServiceRegistry/interface/Service.h" -#include "HGCal/DataFormats/interface/HGCalTBRecHitCollections.h" #include "HGCal/DataFormats/interface/HGCalTBDetId.h" -#include "HGCal/DataFormats/interface/HGCalTBRecHit.h" #include "HGCal/Geometry/interface/HGCalTBCellVertices.h" +#include "HGCal/Geometry/interface/HGCalTBCellParameters.h" #include "HGCal/Geometry/interface/HGCalTBTopology.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" #include "HGCal/CondObjects/interface/HGCalElectronicsMap.h" @@ -43,8 +44,6 @@ #include "HGCal/DataFormats/interface/HGCalTBDataFrameContainers.h" #include "HGCal/Geometry/interface/HGCalTBGeometryParameters.h" -using namespace std; - // // class declaration @@ -68,7 +67,6 @@ class DigiPlotter : public edm::one::EDAnalyzer void analyze(const edm::Event& , const edm::EventSetup&) override; virtual void endJob() override; // ----------member data --------------------------- - bool DEBUG = 0; HGCalTBTopology IsCellValid; HGCalTBCellVertices TheCell; std::string mapfile_ = "HGCal/CondObjects/data/map_CERN_8Layers_Sept2016.txt"; @@ -76,25 +74,21 @@ class DigiPlotter : public edm::one::EDAnalyzer HGCalElectronicsMap emap_; } essource_; int sensorsize = 128;// The geometry for a 256 cell sensor hasnt been implemted yet. Need a picture to do this. - std::vector> CellXY; - std::pair CellCentreXY; + std::vector>::const_iterator it; - const static int NSAMPLES = 2; - TH2Poly *h_digi_layer[NSAMPLES][MAXLAYERS]; - TH1F *h_digi_layer_summed[NSAMPLES][MAXLAYERS]; - TProfile *h_digi_layer_profile[NSAMPLES][MAXLAYERS]; - const static int cellx = 15; - const static int celly = 15; + TH2Poly *h_digi_layer[SKIROC::MAXSAMPLES][MAXLAYERS]; + TH1F *h_digi_layer_summed[SKIROC::MAXSAMPLES][MAXLAYERS]; + TProfile *h_digi_layer_profile[SKIROC::MAXSAMPLES][MAXLAYERS]; int Sensor_Iu = 0; int Sensor_Iv = 0; - TH2F* Noise_2D_Profile[NSAMPLES][MAXLAYERS]; - TH1F *h_digi_layer_channel[MAXSKIROCS][64][NSAMPLES]; -// TH1F *h_digi_layer_cell_event[NSAMPLES][MAXLAYERS][cellx][celly][512]; + TH2F* Noise_2D_Profile[SKIROC::MAXSAMPLES][MAXLAYERS]; + TH1F *h_digi_layer_channel[MAXSKIROCS][64][SKIROC::MAXSAMPLES]; +// TH1F *h_digi_layer_cell_event[SKIROC::MAXSAMPLES][MAXLAYERS][cellx][celly][512]; char name[50], title[50]; double ADC_Sum_SKI_Layer[2][MAXLAYERS][2]; // 2 SKIROCs per layer, High gain and low gain ADC HARD CODED int Cell_Count_SKI_Layer[2][4]; // 2 SKIROCs per layer, High gain and low gain ADC HARD CODED - string m_pedestalsHighGain; - string m_pedestalsLowGain; + std::string m_pedestalsHighGain; + std::string m_pedestalsLowGain; }; // @@ -114,77 +108,67 @@ DigiPlotter::DigiPlotter(const edm::ParameterSet& iConfig) usesResource("TFileService"); edm::Service fs; consumesMany(); - const int HalfHexVertices = 4; - double HalfHexX[HalfHexVertices] = {0.}; - double HalfHexY[HalfHexVertices] = {0.}; - const int FullHexVertices = 6; - double FullHexX[FullHexVertices] = {0.}; - double FullHexY[FullHexVertices] = {0.}; - for(int nsample = 0; nsample < NSAMPLES; nsample++) { + + double CellXs[HGCAL_TB_CELL::FullHexVertices] = {0.}; + double CellYs[HGCAL_TB_CELL::FullHexVertices] = {0.}; + + for(int i_sample = 0; i_sample < SKIROC::MAXSAMPLES; i_sample++) { for(int nlayers = 0; nlayers < MAXLAYERS; nlayers++) { - sprintf(name, "Noise_2D_Profile_ADC%i_Layer%i", nsample, nlayers); - sprintf(title, "Noise 2D Profile ADC%i Layer%i", nsample, nlayers); - Noise_2D_Profile[nsample][nlayers] = fs->make(name, title, 128, 0, 127, 2000, -1000, 1000); + sprintf(name, "Noise_2D_Profile_ADC%i_Layer%i", i_sample, nlayers); + sprintf(title, "Noise 2D Profile ADC%i Layer%i", i_sample, nlayers); + Noise_2D_Profile[i_sample][nlayers] = fs->make(name, title, 128, 0, 127, 2000, -1000, 1000); } } + for(int ISkiroc = 1; ISkiroc <= MAXSKIROCS; ISkiroc++) { - for(int Channel = 0; Channel < 64; Channel++) { - for(int iii = 0; iii < NSAMPLES; iii++) { + for(int Channel = 0; Channel < SKIROC::NCHANNELS; Channel++) { + for(int iii = 0; iii < SKIROC::MAXSAMPLES; iii++) { sprintf(name, "Ski_%i_Channel_%i_ADC%i", ISkiroc, Channel, iii); sprintf(title, "Ski %i Channel %i ADC%i", ISkiroc, Channel, iii); h_digi_layer_channel[ISkiroc - 1][Channel][iii] = fs->make(name, title, 4096, 0., 4095.); } } } - int iii = 0; - for(int nsample = 0; nsample < NSAMPLES; nsample++) { + + for(int i_sample = 0; i_sample < SKIROC::MAXSAMPLES; i_sample++) {// one histogram for each sample (high and low gain) for(int nlayers = 0; nlayers < MAXLAYERS; nlayers++) { -//Booking a "hexagonal" histograms to display the sum of Digis for NSAMPLES, in 1 SKIROC in 1 layer. To include all layers soon. Also the 1D Digis per cell in a sensor is booked here for NSAMPLES. - sprintf(name, "FullLayer_ADC%i_Layer%i", nsample, nlayers + 1); - sprintf(title, "Sum of adc counts per cell for ADC%i Layer%i", nsample, nlayers + 1); - h_digi_layer[nsample][nlayers] = fs->make(); - h_digi_layer[nsample][nlayers]->SetName(name); - h_digi_layer[nsample][nlayers]->SetTitle(title); - sprintf(name, "FullLayer_ADC%i_Layer%i_summed", nsample, nlayers + 1); - sprintf(title, "Sum of adc counts for all cells in ADC%i Layer%i", nsample, nlayers + 1); - h_digi_layer_summed[nsample][nlayers] = fs->make(name, title, 4096, 0., 4095.); - h_digi_layer_summed[nsample][nlayers]->GetXaxis()->SetTitle("Digis[adc counts]"); - sprintf(name, "FullLayer_ADC%i_Layer%i_profile", nsample, nlayers + 1); - sprintf(title, "profile of adc counts for all cells in ADC%i Layer%i", nsample, nlayers + 1); - h_digi_layer_profile[nsample][nlayers] = fs->make(name, title, 128, 0, 127, 0., 4095.); - h_digi_layer_profile[nsample][nlayers]->GetXaxis()->SetTitle("Channel #"); - h_digi_layer_profile[nsample][nlayers]->GetYaxis()->SetTitle("ADC counts"); +//Booking a "hexagonal" histograms to display the sum of Digis for SKIROC::MAXSAMPLES, in 1 SKIROC in 1 layer. To include all layers soon. Also the 1D Digis per cell in a sensor is booked here for SKIROC::MAXSAMPLES. + sprintf(name, "FullLayer_ADC%i_Layer%i", i_sample, nlayers + 1); + sprintf(title, "Sum of adc counts per cell for ADC%i Layer%i", i_sample, nlayers + 1); + h_digi_layer[i_sample][nlayers] = fs->make(); + h_digi_layer[i_sample][nlayers]->SetName(name); + h_digi_layer[i_sample][nlayers]->SetTitle(title); + sprintf(name, "FullLayer_ADC%i_Layer%i_summed", i_sample, nlayers + 1); + sprintf(title, "Sum of adc counts for all cells in ADC%i Layer%i", i_sample, nlayers + 1); + h_digi_layer_summed[i_sample][nlayers] = fs->make(name, title, 4096, 0., 4095.); + h_digi_layer_summed[i_sample][nlayers]->GetXaxis()->SetTitle("Digis[adc counts]"); + sprintf(name, "FullLayer_ADC%i_Layer%i_profile", i_sample, nlayers + 1); + sprintf(title, "profile of adc counts for all cells in ADC%i Layer%i", i_sample, nlayers + 1); + h_digi_layer_profile[i_sample][nlayers] = fs->make(name, title, 128, 0, 127, 0., 4095.); + h_digi_layer_profile[i_sample][nlayers]->GetXaxis()->SetTitle("Channel #"); + h_digi_layer_profile[i_sample][nlayers]->GetYaxis()->SetTitle("ADC counts"); for(int iv = -7; iv < 8; iv++) { for(int iu = -7; iu < 8; iu++) { if(!IsCellValid.iu_iv_valid(nlayers, Sensor_Iu, Sensor_Iv, iu, iv, sensorsize)) continue; - CellXY = TheCell.GetCellCoordinatesForPlots(nlayers, Sensor_Iu, Sensor_Iv, iu, iv, sensorsize); - int NumberOfCellVertices = CellXY.size(); - iii = 0; - if(NumberOfCellVertices == 4) { - for(it = CellXY.begin(); it != CellXY.end(); it++) { - HalfHexX[iii] = it->first; - HalfHexY[iii++] = it->second; - } -//Somehow cloning of the TH2Poly was not working. Need to look at it. Currently physically booked another one. - h_digi_layer[nsample][nlayers]->AddBin(NumberOfCellVertices, HalfHexX, HalfHexY); - } else if(NumberOfCellVertices == 6) { - iii = 0; - for(it = CellXY.begin(); it != CellXY.end(); it++) { - FullHexX[iii] = it->first; - FullHexY[iii++] = it->second; - } - h_digi_layer[nsample][nlayers]->AddBin(NumberOfCellVertices, FullHexX, FullHexY); + std::vector> CellXY = TheCell.GetCellCoordinatesForPlots(nlayers, Sensor_Iu, Sensor_Iv, iu, iv, sensorsize); + size_t vertex_idx = 0; + for(it = CellXY.begin(); it != CellXY.end(); it++) { + CellXs[vertex_idx] = it->first; + CellYs[vertex_idx++] = it->second; } +//Somehow cloning of the TH2Poly was not working. Need to look at it. Currently physically booked another one. + h_digi_layer[i_sample][nlayers]->AddBin(CellXY.size(), CellXs, CellYs); + }//loop over iu }//loop over iv }//loop over nlayers - }//loop over nsamples - - m_pedestalsHighGain = iConfig.getUntrackedParameter("pedestalsHighGain", ""); - m_pedestalsLowGain = iConfig.getUntrackedParameter("pedestalsLowGain", ""); + }//loop over i_samples + + m_pedestalsHighGain = iConfig.getUntrackedParameter("pedestalsHighGain", ""); + m_pedestalsLowGain = iConfig.getUntrackedParameter("pedestalsLowGain", ""); }//contructor ends here @@ -228,16 +212,12 @@ DigiPlotter::analyze(const edm::Event& event, const edm::EventSetup& setup) //////////////////////////////////Evaluate average pedestal per event to subtract out////////////////////////////////// for(SKIROC2DigiCollection::const_iterator k = Coll.begin(); k != Coll.end(); k++) { - const SKIROC2DataFrame& SKI_1 = *k ; - int n_layer = (SKI_1.detid()).layer(); - int n_sensor_IU = (SKI_1.detid()).sensorIU(); - int n_sensor_IV = (SKI_1.detid()).sensorIV(); - int n_cell_iu = (SKI_1.detid()).iu(); - int n_cell_iv = (SKI_1.detid()).iv(); - uint32_t EID = essource_.emap_.detId2eid(SKI_1.detid()); + const SKIROC2DataFrame& SKI = *k ; + const HGCalTBDetId& detId = SKI.detid(); + uint32_t EID = essource_.emap_.detId2eid(SKI.detid()); HGCalTBElectronicsId eid(EID); - if(DEBUG) cout << endl << " Layer = " << n_layer << " Sensor IU = " << n_sensor_IU << " Sensor IV = " << n_sensor_IV << " Cell iu = " << n_cell_iu << " Cell iu = " << n_cell_iv << endl; - if(!IsCellValid.iu_iv_valid(n_layer, n_sensor_IU, n_sensor_IV, n_cell_iu, n_cell_iv, sensorsize)) continue; + if(!IsCellValid.iu_iv_valid(detId, sensorsize)) continue; +; // ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1][1] += SKI_1[0].adcHigh(); // ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1][0] += SKI_1[0].adcLow(); // Cell_Count_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1] += 1; @@ -245,35 +225,32 @@ DigiPlotter::analyze(const edm::Event& event, const edm::EventSetup& setup) -// cout << "SKIROC2 Digis: " << i->provenance()->branchName() << endl; for(SKIROC2DigiCollection::const_iterator j = Coll.begin(); j != Coll.end(); j++) { const SKIROC2DataFrame& SKI = *j ; - int n_layer = (SKI.detid()).layer(); - int n_sensor_IU = (SKI.detid()).sensorIU(); - int n_sensor_IV = (SKI.detid()).sensorIV(); - int n_cell_iu = (SKI.detid()).iu(); - int n_cell_iv = (SKI.detid()).iv(); + const HGCalTBDetId& detId = SKI.detid(); uint32_t EID = essource_.emap_.detId2eid(SKI.detid()); HGCalTBElectronicsId eid(EID); - if(DEBUG) cout << endl << " Layer = " << n_layer << " Sensor IU = " << n_sensor_IU << " Sensor IV = " << n_sensor_IV << " Cell iu = " << n_cell_iu << " Cell iu = " << n_cell_iv << endl; - if(!IsCellValid.iu_iv_valid(n_layer, n_sensor_IU, n_sensor_IV, n_cell_iu, n_cell_iv, sensorsize)) continue; - CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots(n_layer, n_sensor_IU, n_sensor_IV, n_cell_iu, n_cell_iv, sensorsize); +#ifdef DEBUG + std::cout << std::endl << " Layer = " << detId.layer() << " Sensor IU = " << detId.sensorIU() << " Sensor IV = " << detId.sensorIV() << " Cell iu = " << detId.iu() << " Cell iu = " << detId.iv() << std::endl; +#endif + if(!IsCellValid.iu_iv_valid(detId, sensorsize)) continue; + std::pair CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots(detId, sensorsize); double iux = (CellCentreXY.first < 0 ) ? (CellCentreXY.first + 0.0001) : (CellCentreXY.first - 0.0001) ; double iyy = (CellCentreXY.second < 0 ) ? (CellCentreXY.second + 0.0001) : (CellCentreXY.second - 0.0001); int nsample = 0; - h_digi_layer[nsample][n_layer - 1]->Fill(iux , iyy, SKI[nsample].adcLow()); - h_digi_layer_profile[nsample][n_layer - 1]->Fill(counter1++, SKI[nsample].adcLow(), 1); -// h_digi_layer_summed[nsample][n_layer - 1]->Fill(ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1][0]); + h_digi_layer[nsample][detId.layer() - 1]->Fill(iux , iyy, SKI[nsample].adcLow()); + h_digi_layer_profile[nsample][detId.layer() - 1]->Fill(counter1++, SKI[nsample].adcLow(), 1); +// h_digi_layer_summed[nsample][detId.layer() - 1]->Fill(ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(detId.layer() - 1) - 1][detId.layer() - 1][0]); if(eid.iskiroc() > 0) h_digi_layer_channel[eid.iskiroc() - 1][eid.ichan()][nsample]->Fill(SKI[nsample].adcLow()); nsample = 1; - h_digi_layer[nsample][n_layer - 1]->Fill(iux , iyy, SKI[nsample - 1].adcHigh()); -// Noise_2D_Profile[nsample][n_layer - 1]->Fill(); - h_digi_layer_profile[nsample][n_layer - 1]->Fill(counter2++, SKI[nsample - 1].adcHigh(), 1); -// h_digi_layer_summed[nsample][n_layer - 1]->Fill(ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1][1]); -// if(((SKI.detid()).cellType() != 4) && (eid.ichan() == 0) ) cout<GetMean() << endl; - fs2 << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][0]->GetMean() << endl; + fs1 << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][1]->GetMean() << std::endl; + fs2 << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][0]->GetMean() << std::endl; } } } @@ -330,8 +307,8 @@ DigiPlotter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) edm::ParameterSetDescription desc; // desc.setUnknown(); // descriptions.addDefault(desc); - desc.addUntracked("pedestalsHighGain", ""); - desc.addUntracked("pedestalsLowGain", ""); + desc.addUntracked("pedestalsHighGain", ""); + desc.addUntracked("pedestalsLowGain", ""); descriptions.add("hgcaltbdigisplotter", desc); } diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.cc b/RawToDigi/plugins/HGCalTBRawToDigi.cc index 6afe8f7..28d44c8 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.cc +++ b/RawToDigi/plugins/HGCalTBRawToDigi.cc @@ -1,8 +1,15 @@ -#include #include "HGCal/RawToDigi/plugins/HGCalTBRawToDigi.h" + +// provides the maximum number of layers and number of skirocs per board #include "HGCal/Geometry/interface/HGCalTBGeometryParameters.h" -#include "HGCal/Geometry/interface/HGCalTBSpillParameters.h" -using namespace std; + +// provide the number of channels in a skiroc and number of samples +#include "HGCal/DataFormats/interface/SKIROCParameters.h" + +//#define DEBUG +#ifdef DEBUG +#include +#endif unsigned int gray_to_binary (unsigned int gray); @@ -42,9 +49,7 @@ void HGCalTBRawToDigi::produce(edm::Event& e, const edm::EventSetup& c) digis = std::auto_ptr(new SKIROC2DigiCollection(nSkirocs * SKIROC::NCHANNELS * SKIROC::MAXSAMPLES)); const uint16_t* pdata = (const uint16_t*)(fed.data()); - // we start from the back... - size_t ptr = fed.size() / sizeof(uint16_t) - 1; - size_t ski = 0; + size_t ski = 0; // the skirocs have an absolute numbering, start counting from the first board till the last for (unsigned int i_board = 0 ; i_board < nBoards; ++i_board) { for(size_t i_skiroc = 0; i_skiroc < MAXSKIROCS_PER_BOARD; ++i_skiroc) { for (int ichan = 0; ichan < SKIROC::NCHANNELS; ichan++) { @@ -53,13 +58,13 @@ void HGCalTBRawToDigi::produce(edm::Event& e, const edm::EventSetup& c) HGCalTBDetId did = essource_.emap_.eid2detId(eid); digis->addDataFrame(did); #ifdef DEBUG - if(i_board == 0 && i_skiroc == 0) std::cout << (pdata & 0xFFF) << (pdata++ & 0xFFF) << (pdata & 0xFFF) << std::endl; + if(i_board == 0) std::cout << (*pdata & 0xFFF) << "\t" << (*(pdata+1) & 0xFFF) << "\t" << (*(pdata+2) & 0xFFF) << std::endl; #endif - digis->backDataFrame().setSample(0, gray_to_binary(pdata++ & 0xFFF), gray_to_binary( pdata++ & 0xFFF), 0); - + digis->backDataFrame().setSample(0, gray_to_binary(*(pdata) & 0xFFF), gray_to_binary( *(pdata+1) & 0xFFF), 0); + pdata++;pdata++; } } - ++ski; + ++ski; //increment the absolute ID of the skiroc } } diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.h b/RawToDigi/plugins/HGCalTBRawToDigi.h index bb2b65a..6961b30 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.h +++ b/RawToDigi/plugins/HGCalTBRawToDigi.h @@ -12,7 +12,6 @@ #include "HGCal/DataFormats/interface/HGCalTBDataFrameContainers.h" -#include "HGCal/DataFormats/interface/SKIROCParameters.h" #include diff --git a/RawToDigi/plugins/HGCalTBTextSource.cc b/RawToDigi/plugins/HGCalTBTextSource.cc index a889134..356155d 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.cc +++ b/RawToDigi/plugins/HGCalTBTextSource.cc @@ -4,8 +4,6 @@ #include "stdlib.h" #include "HGCal/RawToDigi/plugins/HGCalTBTextSource.h" #include "HGCal/Geometry/interface/HGCalTBGeometryParameters.h" -#include "HGCal/Geometry/interface/HGCalTBSpillParameters.h" -using namespace std; #define TXTLINES 69 //#define DEBUG @@ -42,7 +40,7 @@ bool HGCalTBTextSource::setRunAndEventInfo(edm::EventID& id, edm::TimeValue_t& t id = edm::EventID(m_run, m_spill, m_event); time = (edm::TimeValue_t) m_time; - else return false; + return true; } // sets m_run, m_spill, m_even, m_time @@ -75,7 +73,9 @@ bool HGCalTBTextSource::readLines() { max_boards = 0; for( auto& board : m_lines) { - board.clear(); + for( auto& skiroc : board){ + skiroc.clear(); + } } char buff[1024], buff_SK[1024]; @@ -86,25 +86,26 @@ bool HGCalTBTextSource::readLines() for(unsigned int i = 0; i < TXTLINES && !feof(m_file); ++i) { buff[0] = 0; fgets(buff, 1000, m_file); + if((i < 1) || (i > 64)) continue;// Only these are data words, the rest dont interest us for now // loop over one line of the text file std::string b = buff; std::istringstream buffer(b); unsigned int board_counter = 0; - if((i < 1) || (i > 64)) continue;// Only these are data words, the rest dont interest us for now - while( buffer.peek() != EOF && buffer.good()) { + + while( buffer.peek() != '\n' && buffer.good()) { //continue for all the boards for(size_t iSkiroc = 0; iSkiroc < MAXSKIROCS_PER_BOARD; ++iSkiroc) { // read the data of the one skiroc of one board buffer >> buff_SK; data_sk = strtoul(buff_SK, NULL, 0); m_lines[board_counter][iSkiroc].push_back(data_sk); - //m_lines[board_counter].push_back(data_sk1); } ++board_counter; if(board_counter > max_boards) max_boards = board_counter; } + } - return !m_lines.empty(); + return true; } /** DATA in the FEDRAWDATA in the following order: @@ -118,10 +119,10 @@ void HGCalTBTextSource::produce(edm::Event & e) std::vector skiwords; // make sure there are an even number of 32-bit-words (a round number of 64 bit words... - for (unsigned int i_board = 0 ; i_board < max_boards; ++i_board) { + for (size_t i_board = 0 ; i_board < max_boards; ++i_board) { for(size_t i_skiroc = 0; i_skiroc < MAXSKIROCS_PER_BOARD; ++i_skiroc) { - auto board = m_lines[i_board][i_skiroc]; - for (auto skiword : board) { + auto& board = m_lines[i_board][i_skiroc]; + for (auto& skiword : board) { skiwords.push_back(skiword >> 16); skiwords.push_back(skiword); } diff --git a/RawToDigi/plugins/HGCalTBTextSource.h b/RawToDigi/plugins/HGCalTBTextSource.h index 8c73055..cf29903 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.h +++ b/RawToDigi/plugins/HGCalTBTextSource.h @@ -78,7 +78,7 @@ class HGCalTBTextSource : public edm::ProducerSourceFromFiles public: explicit HGCalTBTextSource(const edm::ParameterSet & pset, edm::InputSourceDescription const& desc) : edm::ProducerSourceFromFiles(pset, desc, true), - m_file(0), + m_file(0) { m_sourceId = pset.getUntrackedParameter("fed", 1000); /// \todo read from file produces(); @@ -103,7 +103,7 @@ class HGCalTBTextSource : public edm::ProducerSourceFromFiles bool readHeader(void); bool readLines(void); - std::array< std::array< std::vector < unsigned int> , MAXLAYERS>, MAXSKIROCS_PER_BOARD> m_lines; + std::array< std::array< std::vector < unsigned int> , MAXSKIROCS_PER_BOARD >, MAXLAYERS > m_lines; FILE* m_file; unsigned int m_time; unsigned int m_event, m_run, m_spill, max_boards; diff --git a/Reco/plugins/Layer_Sum_Analyzer.cc b/Reco/plugins/Layer_Sum_Analyzer.cc index 5ea551d..25645c3 100644 --- a/Reco/plugins/Layer_Sum_Analyzer.cc +++ b/Reco/plugins/Layer_Sum_Analyzer.cc @@ -40,7 +40,6 @@ #include "HGCal/DataFormats/interface/HGCalTBDataFrameContainers.h" #include "HGCal/Geometry/interface/HGCalTBCellVertices.h" #include "HGCal/Geometry/interface/HGCalTBCellParameters.h" -#include "HGCal/Geometry/interface/HGCalTBSpillParameters.h" // chooses which particle to look at. Inverts threshold filtering. // if nothing is selected, electrons are the default diff --git a/Reco/plugins/RecHitPlotter_HighGain_New.cc b/Reco/plugins/RecHitPlotter_HighGain_New.cc index 100a6f8..debde56 100644 --- a/Reco/plugins/RecHitPlotter_HighGain_New.cc +++ b/Reco/plugins/RecHitPlotter_HighGain_New.cc @@ -47,7 +47,6 @@ #include "HGCal/CondObjects/interface/HGCalCondObjectTextIO.h" #include "HGCal/DataFormats/interface/HGCalTBElectronicsId.h" #include "HGCal/Geometry/interface/HGCalTBGeometryParameters.h" -#include "HGCal/Geometry/interface/HGCalTBSpillParameters.h" #define MAXVERTICES 6 using namespace std; @@ -317,7 +316,7 @@ RecHitPlotter_HighGain_New::analyze(const edm::Event& event, const edm::EventSet } TH2Poly *h_RecHit_layer[MAXLAYERS]; - int evId = event.id().event() - 1; + int evId = event.id().event() - 1;// eventID starts from 1 int iLayer = (evId % (MAXLAYERS * EVENTSPERSPILL)) / (EVENTSPERSPILL); cout << endl << " iLayer= " << iLayer << endl; h_RecHit_layer[iLayer] = fs->make(); diff --git a/test_cfg.py b/test_cfg.py index 2cd155e..2dbd8b0 100644 --- a/test_cfg.py +++ b/test_cfg.py @@ -6,7 +6,7 @@ options = VarParsing.VarParsing('standard') # avoid the options: maxEvents, files, secondaryFiles, output, secondaryOutput because they are already defined in 'standard' options.register('dataFolder', - '/afs/cern.ch/work/r/rchatter/Final_Event_Builder/CMSSW_8_0_1/src/HGCal/tmpOut/', + './', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string, 'folder containing raw text input') @@ -67,6 +67,7 @@ os.system("mkdir -p " + options.outputFolder) if (options.runType != "PED" and options.runType != "HGCRun"): + print options sys.exit("Error: only runtypes PED and HGCRun supported for now; given runType was %s"%(options.runType)) if (options.runType == "PED"): From 86e47e3fb879d66a883e943b33f998090cd07857 Mon Sep 17 00:00:00 2001 From: Shervin Nourbakhsh Date: Mon, 17 Oct 2016 15:29:58 +0200 Subject: [PATCH 4/7] More fixes --- RawToDigi/plugins/DigiPlotter.cc | 27 +++++++++---------- test_cfg.py | 45 ++++++++++++++++++++------------ 2 files changed, 41 insertions(+), 31 deletions(-) diff --git a/RawToDigi/plugins/DigiPlotter.cc b/RawToDigi/plugins/DigiPlotter.cc index f6fcfb2..cd57004 100644 --- a/RawToDigi/plugins/DigiPlotter.cc +++ b/RawToDigi/plugins/DigiPlotter.cc @@ -234,20 +234,17 @@ DigiPlotter::analyze(const edm::Event& event, const edm::EventSetup& setup) std::cout << std::endl << " Layer = " << detId.layer() << " Sensor IU = " << detId.sensorIU() << " Sensor IV = " << detId.sensorIV() << " Cell iu = " << detId.iu() << " Cell iu = " << detId.iv() << std::endl; #endif if(!IsCellValid.iu_iv_valid(detId, sensorsize)) continue; - std::pair CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots(detId, sensorsize); + std::pair CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots(detId, sensorsize); double iux = (CellCentreXY.first < 0 ) ? (CellCentreXY.first + 0.0001) : (CellCentreXY.first - 0.0001) ; double iyy = (CellCentreXY.second < 0 ) ? (CellCentreXY.second + 0.0001) : (CellCentreXY.second - 0.0001); + int nsample = 0; h_digi_layer[nsample][detId.layer() - 1]->Fill(iux , iyy, SKI[nsample].adcLow()); h_digi_layer_profile[nsample][detId.layer() - 1]->Fill(counter1++, SKI[nsample].adcLow(), 1); -// h_digi_layer_summed[nsample][detId.layer() - 1]->Fill(ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(detId.layer() - 1) - 1][detId.layer() - 1][0]); if(eid.iskiroc() > 0) h_digi_layer_channel[eid.iskiroc() - 1][eid.ichan()][nsample]->Fill(SKI[nsample].adcLow()); nsample = 1; h_digi_layer[nsample][detId.layer() - 1]->Fill(iux , iyy, SKI[nsample - 1].adcHigh()); -// Noise_2D_Profile[nsample][detId.layer() - 1]->Fill(); h_digi_layer_profile[nsample][detId.layer() - 1]->Fill(counter2++, SKI[nsample - 1].adcHigh(), 1); -// h_digi_layer_summed[nsample][detId.layer() - 1]->Fill(ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(detId.layer() - 1) - 1][detId.layer() - 1][1]); -// if(((SKI.detid()).cellType() != 4) && (eid.ichan() == 0) ) std::cout<GetMean() << std::endl; - fs2 << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][0]->GetMean() << std::endl; + f_highGain << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][1]->GetMean() << std::endl; + f_lowGain << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][0]->GetMean() << std::endl; } } } diff --git a/test_cfg.py b/test_cfg.py index 2dbd8b0..0309f6f 100644 --- a/test_cfg.py +++ b/test_cfg.py @@ -36,9 +36,9 @@ 'Type of run: HGCRun for run with beam on, PED for pedestal run, Unknown otherwise') options.register('chainSequence', - 0, - VarParsing.VarParsing.multiplicity.singleton, - VarParsing.VarParsing.varType.int, + "", + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, '0: if runType is PED then do Digi, if runType is HGC_Run then do Digi and Reco (not implemented yet); 1: do Digi; 2: only Reco (not implemented yet); 3: Digi + highgain_correlation_cm; 4: event display sequence; 5: highgain_correlation_cm + event display sequence') @@ -95,7 +95,10 @@ ) ## [source module] - +process.options = cms.untracked.PSet( + wantSummary = cms.untracked.bool(True), +# SkipEvent = cms.untracked.vstring('ProductNotFound'), +) ###### process.hgcaltbdigisplotter.pedestalsHighGain = cms.untracked.string(options.pedestalsHighGain) process.hgcaltbdigisplotter.pedestalsLowGain = cms.untracked.string(options.pedestalsLowGain) @@ -115,11 +118,11 @@ fileName = cms.untracked.string(options.output) ) -# process.TFileService = cms.Service("TFileService", fileName = cms.string("HGC_Output_6_Reco_Display.root") ) -if (options.chainSequence == 1): - process.TFileService = cms.Service("TFileService", fileName = cms.string("%s/%s_Output_%06d_Digi.root"%(options.outputFolder,options.runType,options.runNumber))) -elif (options.chainSequence == 3 or options.chainSequence == 4 or options.chainSequence == 5): - process.TFileService = cms.Service("TFileService", fileName = cms.string("%s/%s_Output_%06d_Reco.root"%(options.outputFolder,options.runType,options.runNumber))) +process.TFileService = cms.Service("TFileService", fileName = cms.string("HGC_Output.root") ) +#if (options.chainSequence == 1): +# process.TFileService = cms.Service("TFileService", fileName = cms.string("%s/%s_Output_%06d_Digi.root"%(options.outputFolder,options.runType,options.runNumber))) +#elif (options.chainSequence == 3 or options.chainSequence == 4 or options.chainSequence == 5): +# process.TFileService = cms.Service("TFileService", fileName = cms.string("%s/%s_Output_%06d_Reco.root"%(options.outputFolder,options.runType,options.runNumber))) # process.TFileService = cms.Service("TFileService", fileName = cms.string("HGC_Output_6_Reco.root") ) #process.TFileService = cms.Service("TFileService", fileName = cms.string("HGC_Output_6_Reco_Layer.root") ) #process.TFileService = cms.Service("TFileService", fileName = cms.string("HGC_Output_6_Reco_Cluster.root") ) @@ -140,13 +143,21 @@ ################Miscellaneous############################################################################## #process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbrechits*process.FourLayerRecHitPlotterMax) -if (options.chainSequence == 1): - process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbdigisplotter) -elif (options.chainSequence == 3): - process.p =cms.Path(process.hgcaltbdigis) -elif (options.chainSequence == 4): - process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbrechits*process.hgcaltbrechitsplotter_highgain_new) -elif (options.chainSequence == 5): - process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbrechits*process.hgcaltbrechitsplotter_highgain_correlation_cm*process.hgcaltbrechitsplotter_highgain_new) +process.p = cms.Path() +for seq in options.chainSequence: + if(seq == "DIGI"): + process.p *= process.hgcaltbdigis + if(seq == "DIGIPLOT" or seq == "PED"): + process.p *= process.hgcaltbdigisplotter + if(seq == "RECO"): + process.p *= process.hgcaltbrechits +# if (options.chainSequence == 1): +# process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbdigisplotter) +# elif (options.chainSequence == 3): +# process.p =cms.Path(process.hgcaltbdigis) +# elif (options.chainSequence == 4): +# process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbrechits*process.hgcaltbrechitsplotter_highgain_new) +# elif (options.chainSequence == 5): +# process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbrechits*process.hgcaltbrechitsplotter_highgain_correlation_cm*process.hgcaltbrechitsplotter_highgain_new) process.end = cms.EndPath(process.output) From 6309a27c5433e7262ac822eee164711f9074cb23 Mon Sep 17 00:00:00 2001 From: Shervin Nourbakhsh Date: Mon, 17 Oct 2016 15:40:09 +0200 Subject: [PATCH 5/7] Adding option not to dump pedestals, otherwise it can be dangerous --- RawToDigi/plugins/DigiPlotter.cc | 8 +++++++- RawToDigi/python/hgcaltbdigis_cfi.py | 2 +- test_cfg.py | 2 ++ 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/RawToDigi/plugins/DigiPlotter.cc b/RawToDigi/plugins/DigiPlotter.cc index cd57004..395a39d 100644 --- a/RawToDigi/plugins/DigiPlotter.cc +++ b/RawToDigi/plugins/DigiPlotter.cc @@ -89,6 +89,7 @@ class DigiPlotter : public edm::one::EDAnalyzer int Cell_Count_SKI_Layer[2][4]; // 2 SKIROCs per layer, High gain and low gain ADC HARD CODED std::string m_pedestalsHighGain; std::string m_pedestalsLowGain; + bool _dumpNewPedestals; }; // @@ -102,7 +103,8 @@ class DigiPlotter : public edm::one::EDAnalyzer // // constructors and destructor // -DigiPlotter::DigiPlotter(const edm::ParameterSet& iConfig) +DigiPlotter::DigiPlotter(const edm::ParameterSet& iConfig): + _dumpNewPedestals(false) { //now do what ever initialization is needed usesResource("TFileService"); @@ -169,6 +171,7 @@ DigiPlotter::DigiPlotter(const edm::ParameterSet& iConfig) m_pedestalsHighGain = iConfig.getUntrackedParameter("pedestalsHighGain", ""); m_pedestalsLowGain = iConfig.getUntrackedParameter("pedestalsLowGain", ""); + _dumpNewPedestals = iConfig.getUntrackedParameter("dumpNewPedestals", false); }//contructor ends here @@ -271,6 +274,7 @@ DigiPlotter::beginJob() void DigiPlotter::endJob() { + if(_dumpNewPedestals){ int Code = 0; int SENSOR_IX = 0; int SENSOR_IV = 0; @@ -294,6 +298,7 @@ DigiPlotter::endJob() } } } + } } // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ @@ -308,6 +313,7 @@ DigiPlotter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) // descriptions.addDefault(desc); desc.addUntracked("pedestalsHighGain", ""); desc.addUntracked("pedestalsLowGain", ""); + desc.addUntracked("dumpNewPedestals", false); descriptions.add("hgcaltbdigisplotter", desc); } diff --git a/RawToDigi/python/hgcaltbdigis_cfi.py b/RawToDigi/python/hgcaltbdigis_cfi.py index 530a1d1..3a4bf70 100644 --- a/RawToDigi/python/hgcaltbdigis_cfi.py +++ b/RawToDigi/python/hgcaltbdigis_cfi.py @@ -5,6 +5,6 @@ fedId=cms.untracked.int32(1000), ## list of FEDs should be know from the setup electronicsMap=cms.untracked.string("HGCal/CondObjects/data/map_CERN_8Layers_Sept2016.txt"), ### electronicsMap written by hand - + dumpNewPedestals = cms.bool(False), ) diff --git a/test_cfg.py b/test_cfg.py index 0309f6f..b203c0f 100644 --- a/test_cfg.py +++ b/test_cfg.py @@ -149,6 +149,8 @@ process.p *= process.hgcaltbdigis if(seq == "DIGIPLOT" or seq == "PED"): process.p *= process.hgcaltbdigisplotter + if(seq == "PED"): + process.hgcaltbdigisplotter.dumpNewPedestals = cms.untracked.bool(True) if(seq == "RECO"): process.p *= process.hgcaltbrechits # if (options.chainSequence == 1): From 7cbae011a1ea8053ad573f9cb6224c518f85d6da Mon Sep 17 00:00:00 2001 From: Shervin Nourbakhsh Date: Mon, 17 Oct 2016 15:45:39 +0200 Subject: [PATCH 6/7] format code and minor fixes --- DataFormats/interface/SKIROCParameters.h | 4 +- Geometry/interface/HGCalTBCellVertices.h | 4 +- RawToDigi/plugins/DigiPlotter.cc | 54 +- RawToDigi/plugins/HGCalTBRawToDigi.cc | 7 +- RawToDigi/plugins/HGCalTBTextSource.cc | 2 +- Reco/plugins/HGCalTBRecHitProducer.cc | 4 +- Reco/plugins/HGCalTBRecHitProducer.h | 4 +- Reco/plugins/Layer_Sum_Analyzer.cc | 589 +++++++++--------- .../RecHitPlotter_HighGain_Correlation_CM.cc | 2 +- Reco/plugins/RecHitPlotter_HighGain_New.cc | 1 + test_cfg.py | 4 +- 11 files changed, 339 insertions(+), 336 deletions(-) diff --git a/DataFormats/interface/SKIROCParameters.h b/DataFormats/interface/SKIROCParameters.h index 8d70fe7..02d8fbe 100644 --- a/DataFormats/interface/SKIROCParameters.h +++ b/DataFormats/interface/SKIROCParameters.h @@ -3,8 +3,8 @@ namespace SKIROC { - static const unsigned char NCHANNELS = 64; ///< number of channels read by one SKIROC - static const unsigned char MAXSAMPLES = 2; ///< one sample for high gain and one for low gain // 15 is the real max of the skiroc +static const unsigned char NCHANNELS = 64; ///< number of channels read by one SKIROC +static const unsigned char MAXSAMPLES = 2; ///< one sample for high gain and one for low gain // 15 is the real max of the skiroc } diff --git a/Geometry/interface/HGCalTBCellVertices.h b/Geometry/interface/HGCalTBCellVertices.h index aa0934c..ad5c1bf 100644 --- a/Geometry/interface/HGCalTBCellVertices.h +++ b/Geometry/interface/HGCalTBCellVertices.h @@ -25,8 +25,8 @@ class HGCalTBCellVertices HGCalTBCellVertices(); ///< Constructor from cell \b iu & \b iv, valid sensorSizes are 128 and 256 std::vector> GetCellCoordinates(int layer, int sensor_iu, int sensor_iv, int iu, int iv, int sensorsize, bool flipX = false); ///< returns the coordinates of each vertex of cell in the lab frame \b (x,y) - - + + inline std::vector> GetCellCoordinatesForPlots(const HGCalTBDetId& detId, int sensorSize) { return GetCellCoordinatesForPlots(detId.layer(), detId.sensorIU(), detId.sensorIV(), detId.iu(), detId.iv(), sensorSize); diff --git a/RawToDigi/plugins/DigiPlotter.cc b/RawToDigi/plugins/DigiPlotter.cc index 395a39d..db5ddcf 100644 --- a/RawToDigi/plugins/DigiPlotter.cc +++ b/RawToDigi/plugins/DigiPlotter.cc @@ -6,7 +6,7 @@ /**\class DigiPlotter DigiPlotter.cc HGCal/DigiPlotter/plugins/DigiPlotter.cc Description: Plugin to make 2D and 1D histograms of digis - + Implementation: \author Rajdeep Mohan Chatterjee \author Shervin Nourbakhsh @@ -162,13 +162,13 @@ DigiPlotter::DigiPlotter(const edm::ParameterSet& iConfig): } //Somehow cloning of the TH2Poly was not working. Need to look at it. Currently physically booked another one. h_digi_layer[i_sample][nlayers]->AddBin(CellXY.size(), CellXs, CellYs); - + }//loop over iu }//loop over iv }//loop over nlayers }//loop over i_samples - + m_pedestalsHighGain = iConfig.getUntrackedParameter("pedestalsHighGain", ""); m_pedestalsLowGain = iConfig.getUntrackedParameter("pedestalsLowGain", ""); _dumpNewPedestals = iConfig.getUntrackedParameter("dumpNewPedestals", false); @@ -220,7 +220,7 @@ DigiPlotter::analyze(const edm::Event& event, const edm::EventSetup& setup) uint32_t EID = essource_.emap_.detId2eid(SKI.detid()); HGCalTBElectronicsId eid(EID); if(!IsCellValid.iu_iv_valid(detId, sensorsize)) continue; -; + ; // ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1][1] += SKI_1[0].adcHigh(); // ADC_Sum_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1][0] += SKI_1[0].adcLow(); // Cell_Count_SKI_Layer[eid.iskiroc() - 2*(n_layer - 1) - 1][n_layer - 1] += 1; @@ -250,7 +250,7 @@ DigiPlotter::analyze(const edm::Event& event, const edm::EventSetup& setup) h_digi_layer_profile[nsample][detId.layer() - 1]->Fill(counter2++, SKI[nsample - 1].adcHigh(), 1); if(eid.iskiroc() > 0) h_digi_layer_channel[eid.iskiroc() - 1][eid.ichan()][nsample]->Fill(SKI[nsample - 1].adcHigh()); } - + } } else { edm::LogWarning("DQM") << "No SKIROC2 Digis"; @@ -274,31 +274,31 @@ DigiPlotter::beginJob() void DigiPlotter::endJob() { - if(_dumpNewPedestals){ - int Code = 0; - int SENSOR_IX = 0; - int SENSOR_IV = 0; - std::ofstream f_highGain, f_lowGain; - - f_highGain.open(m_pedestalsHighGain.c_str()); - f_highGain << "SCHEME_CODE 0" << std::endl; - f_highGain << "# CODE LAYER SENSOR_IX SENSOR_IV IX IV TYPE VALUE" << std::endl; - - f_lowGain.open(m_pedestalsLowGain.c_str()); - f_lowGain << "SCHEME_CODE 0" << std::endl; - f_lowGain << "# CODE LAYER SENSOR_IX SENSOR_IV IX IV TYPE VALUE" << std::endl; - - for(int ISkiroc = 1; ISkiroc <= MAXSKIROCS; ISkiroc++) { - for(int Channel = 0; Channel < SKIROC::NCHANNELS; Channel++) { - HGCalTBElectronicsId ElId(ISkiroc, Channel); - HGCalTBDetId DetId = essource_.emap_.eid2detId(ElId); - if(DetId.layer() != 0) { - f_highGain << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][1]->GetMean() << std::endl; - f_lowGain << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][0]->GetMean() << std::endl; + if(_dumpNewPedestals) { + int Code = 0; + int SENSOR_IX = 0; + int SENSOR_IV = 0; + std::ofstream f_highGain, f_lowGain; + + f_highGain.open(m_pedestalsHighGain.c_str()); + f_highGain << "SCHEME_CODE 0" << std::endl; + f_highGain << "# CODE LAYER SENSOR_IX SENSOR_IV IX IV TYPE VALUE" << std::endl; + + f_lowGain.open(m_pedestalsLowGain.c_str()); + f_lowGain << "SCHEME_CODE 0" << std::endl; + f_lowGain << "# CODE LAYER SENSOR_IX SENSOR_IV IX IV TYPE VALUE" << std::endl; + + for(int ISkiroc = 1; ISkiroc <= MAXSKIROCS; ISkiroc++) { + for(int Channel = 0; Channel < SKIROC::NCHANNELS; Channel++) { + HGCalTBElectronicsId ElId(ISkiroc, Channel); + HGCalTBDetId DetId = essource_.emap_.eid2detId(ElId); + if(DetId.layer() != 0) { + f_highGain << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][1]->GetMean() << std::endl; + f_lowGain << " " << Code << " " << DetId.layer() << " " << SENSOR_IX << " " << SENSOR_IV << " " << DetId.iu() << " " << DetId.iv() << " " << " " << DetId.cellType() << " " << h_digi_layer_channel[ISkiroc - 1][Channel][0]->GetMean() << std::endl; + } } } } - } } // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.cc b/RawToDigi/plugins/HGCalTBRawToDigi.cc index 28d44c8..682503f 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.cc +++ b/RawToDigi/plugins/HGCalTBRawToDigi.cc @@ -58,10 +58,11 @@ void HGCalTBRawToDigi::produce(edm::Event& e, const edm::EventSetup& c) HGCalTBDetId did = essource_.emap_.eid2detId(eid); digis->addDataFrame(did); #ifdef DEBUG - if(i_board == 0) std::cout << (*pdata & 0xFFF) << "\t" << (*(pdata+1) & 0xFFF) << "\t" << (*(pdata+2) & 0xFFF) << std::endl; + if(i_board == 0) std::cout << (*pdata & 0xFFF) << "\t" << (*(pdata + 1) & 0xFFF) << "\t" << (*(pdata + 2) & 0xFFF) << std::endl; #endif - digis->backDataFrame().setSample(0, gray_to_binary(*(pdata) & 0xFFF), gray_to_binary( *(pdata+1) & 0xFFF), 0); - pdata++;pdata++; + digis->backDataFrame().setSample(0, gray_to_binary(*(pdata) & 0xFFF), gray_to_binary( *(pdata + 1) & 0xFFF), 0); + pdata++; + pdata++; } } ++ski; //increment the absolute ID of the skiroc diff --git a/RawToDigi/plugins/HGCalTBTextSource.cc b/RawToDigi/plugins/HGCalTBTextSource.cc index 356155d..660f5d2 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.cc +++ b/RawToDigi/plugins/HGCalTBTextSource.cc @@ -73,7 +73,7 @@ bool HGCalTBTextSource::readLines() { max_boards = 0; for( auto& board : m_lines) { - for( auto& skiroc : board){ + for( auto& skiroc : board) { skiroc.clear(); } } diff --git a/Reco/plugins/HGCalTBRecHitProducer.cc b/Reco/plugins/HGCalTBRecHitProducer.cc index 5f65c58..6921ec0 100644 --- a/Reco/plugins/HGCalTBRecHitProducer.cc +++ b/Reco/plugins/HGCalTBRecHitProducer.cc @@ -17,8 +17,8 @@ HGCalTBRecHitProducer::HGCalTBRecHitProducer(const edm::ParameterSet& cfg) HGCalCondObjectTextIO io(0); edm::FileInPath fip(_mapFile); - if (!io.load(fip.fullPath(), essource_.emap_)) { - throw cms::Exception("Unable to load electronics map"); + if (!io.load(fip.fullPath(), essource_.emap_)) { + throw cms::Exception("Unable to load electronics map"); }; } diff --git a/Reco/plugins/HGCalTBRecHitProducer.h b/Reco/plugins/HGCalTBRecHitProducer.h index 961df11..33e3771 100644 --- a/Reco/plugins/HGCalTBRecHitProducer.h +++ b/Reco/plugins/HGCalTBRecHitProducer.h @@ -46,8 +46,8 @@ class HGCalTBRecHitProducer : public edm::EDProducer std::vector _LG2HG_value; std::string _mapFile; struct { - HGCalElectronicsMap emap_; - } essource_; + HGCalElectronicsMap emap_; + } essource_; }; diff --git a/Reco/plugins/Layer_Sum_Analyzer.cc b/Reco/plugins/Layer_Sum_Analyzer.cc index 25645c3..998ed2d 100644 --- a/Reco/plugins/Layer_Sum_Analyzer.cc +++ b/Reco/plugins/Layer_Sum_Analyzer.cc @@ -60,15 +60,15 @@ double MIP2ParticleCalib = 1.06; // CERN to pion 125GeV //double LayerWeight[16] = {1.4091566745180932, 0.7020676448403224, 0.6054055986179145, 0.7020676448403224, 0.6054055986179145, 0.8415931435769973, 0.8061197656138868, 0.9811186423136724, 0.8061197656138868, 0.9811186423136724, 0.8061197656138868, 1.4462036381025891, 1.2075480996058319, 1.4462036381025891, 1.2075480996058319, 1.4462036381025891}; -double LayerWeight_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; -double X0depth_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; +double LayerWeight_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; +double X0depth_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; double LayerWeight_8L_conf1[16] = {33.074, 13.184, 14.17, 9.788, 9.766, 9.766, 16.339, 14.129, 0., 0., 0., 0., 0., 0., 0., 0.}; double X0depth_8L_conf1[16] = {6.268, 1.131, 1.131, 1.362, 0.574, 1.301, 0.574, 2.42, 0., 0., 0., 0., 0., 0., 0., 0.}; double LayerWeight_8L_conf2[16] = {35.866, 30.864, 28.803, 23.095, 20.657, 19.804, 36.322, 27.451, 0., 0., 0., 0., 0., 0., 0., 0.}; double X0depth_8L_conf2[16] = {5.048, 3.412, 3.412, 2.866, 2.512, 1.625, 2.368, 6.021, 0., 0., 0., 0., 0., 0., 0., 0.}; double weights2GeV = 1.e-03; -double weights2MIP = 52.8/63.6; // rescale weights from mean to MPV +double weights2MIP = 52.8 / 63.6; // rescale weights from mean to MPV //double LayerWeight[16] = {0.4847555727337982, 1.0214605968539232, 0.4847555727337982, 1.0214605968539232, 0.4847555727337982, 1.1420105918768606, 0.6423912113800805, 1.2625605868997982, 0.6423912113800805, 1.2625605868997982, 0.6423912113800805, 1.6643939036429232, 0.9576624886726451, 1.6643939036429232, 0.9576624886726451, 1.6643939036429232};// dE/dx weights @@ -104,7 +104,7 @@ class Layer_Sum_Analyzer : public edm::one::EDAnalyzer> CellXY; std::pair CellCentreXY; - HGCalTBTopology IsCellValid; + HGCalTBTopology IsCellValid; HGCalTBCellVertices TheCell; - double maxdist = (1 + sqrt (3) / 2) * HGCAL_TB_CELL::FULL_CELL_SIDE; // <<< FIXME maxdist > HGCAL_TB_CELL::FULL_CELL_SIDE !! + double maxdist = (1 + sqrt (3) / 2) * HGCAL_TB_CELL::FULL_CELL_SIDE; // <<< FIXME maxdist > HGCAL_TB_CELL::FULL_CELL_SIDE !! - double Weights_L[MAXLAYERS]; - double X0_L[MAXLAYERS]; - double ADCtoMIP[MAXLAYERS]; + double Weights_L[MAXLAYERS]; + double X0_L[MAXLAYERS]; + double ADCtoMIP[MAXLAYERS]; - TH1F *h_CM_layer[MAXSKIROCS]; + TH1F *h_CM_layer[MAXSKIROCS]; TH1F *h_sum_layer[MAXLAYERS], *h_layer_seven[MAXLAYERS], *h_layer_nineteen[MAXLAYERS], *h_Seed_layer[MAXLAYERS]; TH1F *h_sum_layer_AbsW[MAXLAYERS], *h_layer_seven_AbsW[MAXLAYERS], *h_layer_nineteen_AbsW[MAXLAYERS], *h_Seed_layer_AbsW[MAXLAYERS]; - TH1F *h_x_layer[MAXLAYERS], *h_y_layer[MAXLAYERS]; + TH1F *h_x_layer[MAXLAYERS], *h_y_layer[MAXLAYERS]; TH2F *h_x_y_layer[MAXLAYERS]; - TH1F* h_Radius[MAXLAYERS]; - TH1F *h_E1oE7_layer[MAXLAYERS], *h_E1oE19_layer[MAXLAYERS], *h_E7oE19_layer[MAXLAYERS]; + TH1F* h_Radius[MAXLAYERS]; + TH1F *h_E1oE7_layer[MAXLAYERS], *h_E1oE19_layer[MAXLAYERS], *h_E7oE19_layer[MAXLAYERS]; - TH1F *h_sum_all, *h_seven_all, *h_nineteen_all; - TH1F *h_sum_all_AbsW, *h_seven_all_AbsW, *h_nineteen_all_AbsW; - TProfile *tp_E1_vs_layer, *tp_E7_vs_layer, *tp_E19_vs_layer; - TProfile *tp_E1oSumL_vs_layer, *tp_E7oSumL_vs_layer, *tp_E19oSumL_vs_layer; - TProfile *tp_E1oE7_vs_layer, *tp_E1oE19_vs_layer, *tp_E7oE19_vs_layer; + TH1F *h_sum_all, *h_seven_all, *h_nineteen_all; + TH1F *h_sum_all_AbsW, *h_seven_all_AbsW, *h_nineteen_all_AbsW; + TProfile *tp_E1_vs_layer, *tp_E7_vs_layer, *tp_E19_vs_layer; + TProfile *tp_E1oSumL_vs_layer, *tp_E7oSumL_vs_layer, *tp_E19oSumL_vs_layer; + TProfile *tp_E1oE7_vs_layer, *tp_E1oE19_vs_layer, *tp_E7oE19_vs_layer; TH2F *HighGain_LowGain_2D; TH2F *Energy_LowGain_2D; - int EVENT = 0; + int EVENT = 0; map Time_Stamp; map Delta_Time_Stamp; @@ -150,7 +150,7 @@ class Layer_Sum_Analyzer : public edm::one::EDAnalyzer fs; @@ -159,28 +159,28 @@ Layer_Sum_Analyzer::Layer_Sum_Analyzer(const edm::ParameterSet& iConfig) //booking the histos for(int layer = 0; layer < MAXLAYERS; layer++) { - h_sum_layer[layer] = fs->make(Form("sumAll_Layer%d", layer+1), "", 40010, -10, 40000); - h_Seed_layer[layer] = fs->make(Form("h_Seed_layer%d",layer+1), Form("h_Seed_layer%d",layer+1), 40010, -10, 40000); - h_layer_seven[layer] = fs->make(Form("sum7_Layer%d", layer+1),"", 40010, -10, 40000); - h_layer_nineteen[layer] = fs->make(Form("sum19_Layer%d", layer+1), "", 40010, -10, 40000); - - h_sum_layer_AbsW[layer] = fs->make(Form("sumAll_Layer%d_AbsW", layer+1), "", 40010, -10, 40000); - h_Seed_layer_AbsW[layer] = fs->make(Form("h_Seed_layer%d_AbsW",layer+1), Form("h_Seed_layer%d",layer+1), 40010, -10, 40000); - h_layer_seven_AbsW[layer] = fs->make(Form("sum7_Layer%d_AbsW", layer+1),"", 40010, -10, 40000); - h_layer_nineteen_AbsW[layer] = fs->make(Form("sum19_Layer%d_AbsW", layer+1), "", 40010, -10, 40000); - - h_x_layer[layer] = fs->make(Form("X_Layer%d", layer+1), "", 2000, -10., 10. ); - h_y_layer[layer] = fs->make(Form("Y_Layer%d", layer+1), "", 2000, -10., 10. ); - h_x_y_layer[layer] = fs->make(Form("YvsX_Layer%d", layer+1), "", 2000, -10., 10., 2000, -10., 10. ); - - h_E1oE7_layer[layer] = fs->make(Form("h_E1oE7_layer%d",layer+1), Form("h_E1oE7_layer%d",layer+1), 5000, -5, 5); - h_E1oE19_layer[layer] = fs->make(Form("h_E1oE19_layer%d",layer+1), Form("h_E1oE19_layer%d",layer+1), 5000, -5, 5); - h_E7oE19_layer[layer] = fs->make(Form("h_E7oE19_layer%d",layer+1), Form("h_E7oE19_layer%d",layer+1), 5000, -5, 5); - h_Radius[layer] = fs->make(Form("h_Radius_layer%d",layer+1), Form("h_Radius_layer%d",layer+1), 5000, 0., 20.); + h_sum_layer[layer] = fs->make(Form("sumAll_Layer%d", layer + 1), "", 40010, -10, 40000); + h_Seed_layer[layer] = fs->make(Form("h_Seed_layer%d", layer + 1), Form("h_Seed_layer%d", layer + 1), 40010, -10, 40000); + h_layer_seven[layer] = fs->make(Form("sum7_Layer%d", layer + 1), "", 40010, -10, 40000); + h_layer_nineteen[layer] = fs->make(Form("sum19_Layer%d", layer + 1), "", 40010, -10, 40000); + + h_sum_layer_AbsW[layer] = fs->make(Form("sumAll_Layer%d_AbsW", layer + 1), "", 40010, -10, 40000); + h_Seed_layer_AbsW[layer] = fs->make(Form("h_Seed_layer%d_AbsW", layer + 1), Form("h_Seed_layer%d", layer + 1), 40010, -10, 40000); + h_layer_seven_AbsW[layer] = fs->make(Form("sum7_Layer%d_AbsW", layer + 1), "", 40010, -10, 40000); + h_layer_nineteen_AbsW[layer] = fs->make(Form("sum19_Layer%d_AbsW", layer + 1), "", 40010, -10, 40000); + + h_x_layer[layer] = fs->make(Form("X_Layer%d", layer + 1), "", 2000, -10., 10. ); + h_y_layer[layer] = fs->make(Form("Y_Layer%d", layer + 1), "", 2000, -10., 10. ); + h_x_y_layer[layer] = fs->make(Form("YvsX_Layer%d", layer + 1), "", 2000, -10., 10., 2000, -10., 10. ); + + h_E1oE7_layer[layer] = fs->make(Form("h_E1oE7_layer%d", layer + 1), Form("h_E1oE7_layer%d", layer + 1), 5000, -5, 5); + h_E1oE19_layer[layer] = fs->make(Form("h_E1oE19_layer%d", layer + 1), Form("h_E1oE19_layer%d", layer + 1), 5000, -5, 5); + h_E7oE19_layer[layer] = fs->make(Form("h_E7oE19_layer%d", layer + 1), Form("h_E7oE19_layer%d", layer + 1), 5000, -5, 5); + h_Radius[layer] = fs->make(Form("h_Radius_layer%d", layer + 1), Form("h_Radius_layer%d", layer + 1), 5000, 0., 20.); } - for(int ski=0; skimake(Form("h_CM_skiroc%d",ski+1), "", 5000, -500, 500); + for(int ski = 0; ski < MAXSKIROCS; ++ski) { + h_CM_layer[ski] = fs->make(Form("h_CM_skiroc%d", ski + 1), "", 5000, -500, 500); } h_sum_all = fs->make("h_sumAll_AllLayers", "", 40010, -10, 40000); @@ -199,44 +199,43 @@ Layer_Sum_Analyzer::Layer_Sum_Analyzer(const edm::ParameterSet& iConfig) HighGain_LowGain_2D = fs->make("h2_HGvsLG", "", 4000, 0, 4000, 4000, 0, 4000); Energy_LowGain_2D = fs->make("h2_EvsLG", "", 4000, 0, 4000, 4000, 0, 4000); - tp_E1_vs_layer = fs->make("tp_E1_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E7_vs_layer = fs->make("tp_E7_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E19_vs_layer = fs->make("tp_E19_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E1oSumL_vs_layer = fs->make("tp_E1oSumL_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E7oSumL_vs_layer = fs->make("tp_E7oSumL_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E19oSumL_vs_layer = fs->make("tp_E19oSumL_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E1oE7_vs_layer = fs->make("tp_E1oE7_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E1oE19_vs_layer = fs->make("tp_E1oE19_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - tp_E7oE19_vs_layer = fs->make("tp_E7oE19_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); - - - //loading the proper weights - if(MAXLAYERS != 8) { - std::cout << " update weights " << std::endl; - return; + tp_E1_vs_layer = fs->make("tp_E1_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E7_vs_layer = fs->make("tp_E7_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E19_vs_layer = fs->make("tp_E19_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E1oSumL_vs_layer = fs->make("tp_E1oSumL_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E7oSumL_vs_layer = fs->make("tp_E7oSumL_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E19oSumL_vs_layer = fs->make("tp_E19oSumL_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E1oE7_vs_layer = fs->make("tp_E1oE7_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E1oE19_vs_layer = fs->make("tp_E1oE19_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + tp_E7oE19_vs_layer = fs->make("tp_E7oE19_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); + + + //loading the proper weights + if(MAXLAYERS != 8) { + std::cout << " update weights " << std::endl; + return; } - for(int iL=0; iL>> Layer_Sum_Analyzer::analyze " << std::endl; - EVENT = (event.id()).event(); - - if((event.id()).event() == 1) { - Time_Temp = event.time().value(); - Time_Stamp[EVENT] = Time_Temp; - Delta_Time_Stamp[EVENT] = 0.; - } else { - Time_Stamp[EVENT] = event.time().value(); - Delta_Time_Stamp[EVENT] = event.time().value() - Time_Temp; - Time_Temp = event.time().value(); - } - - //opening Rechits - edm::Handle Rechits; - event.getByToken(HGCalTBRecHitCollection_, Rechits); - - // looping over each rechit to fill histogram - double commonmode[MAXSKIROCS]; - int cm_num[MAXSKIROCS]; - for(int iS=0; iS REPLACE WITH EXTERNAL CLASS => FIX - for(auto Rechit : *Rechits){ - - if(!IsCellValid.iu_iv_valid((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize)) continue; - - //getting electronics ID - uint32_t EID = essource_.emap_.detId2eid(Rechit.id()); - HGCalTBElectronicsId eid(EID); - - int n_layer = (Rechit.id()).layer() - 1; - int n_cell_type = (Rechit.id()).cellType(); - int n_skiroc = eid.iskiroc() - 1; - - - //getting X and Y coordinates - CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize); - - HighGain_LowGain_2D->Fill(Rechit.energyLow(), Rechit.energyHigh()); - Energy_LowGain_2D->Fill(Rechit.energy(), Rechit.energyHigh()); - - // needed? FIXME - if(n_cell_type != 0) continue; - - if(Rechit.energy() > max[n_layer]) { - max[n_layer] = Rechit.energy(); - max_x[n_layer] = CellCentreXY.first; - max_y[n_layer] = CellCentreXY.second; - } - - if((Rechit.energy()) / ADCtoMIP[n_skiroc] <= CMTHRESHOLD) { - commonmode[n_skiroc] += Rechit.energy(); - cm_num[n_skiroc]++; - } - - }//Rechit loop ends here - // std::cout << " >>> found commonmode = " << commonmode << std::endl; - - - for(int iS=0; iS>> Layer_Sum_Analyzer::analyze " << std::endl; + EVENT = (event.id()).event(); + + if((event.id()).event() == 1) { + Time_Temp = event.time().value(); + Time_Stamp[EVENT] = Time_Temp; + Delta_Time_Stamp[EVENT] = 0.; + } else { + Time_Stamp[EVENT] = event.time().value(); + Delta_Time_Stamp[EVENT] = event.time().value() - Time_Temp; + Time_Temp = event.time().value(); + } + + //opening Rechits + edm::Handle Rechits; + event.getByToken(HGCalTBRecHitCollection_, Rechits); + + // looping over each rechit to fill histogram + double commonmode[MAXSKIROCS]; + int cm_num[MAXSKIROCS]; + for(int iS = 0; iS < MAXSKIROCS; ++iS) { + commonmode[iS] = 0.; + cm_num[iS] = 0; + } + double max[MAXLAYERS], max_x[MAXLAYERS], max_y[MAXLAYERS]; + for(int iL = 0; iL < MAXLAYERS; ++iL) { + max[iL] = max_x[iL] = max_y[iL] = 0.; + } + + + //CM subtraction => REPLACE WITH EXTERNAL CLASS => FIX + for(auto Rechit : *Rechits) { + + if(!IsCellValid.iu_iv_valid((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize)) continue; + + //getting electronics ID + uint32_t EID = essource_.emap_.detId2eid(Rechit.id()); + HGCalTBElectronicsId eid(EID); + + int n_layer = (Rechit.id()).layer() - 1; + int n_cell_type = (Rechit.id()).cellType(); + int n_skiroc = eid.iskiroc() - 1; + + + //getting X and Y coordinates + CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize); + + HighGain_LowGain_2D->Fill(Rechit.energyLow(), Rechit.energyHigh()); + Energy_LowGain_2D->Fill(Rechit.energy(), Rechit.energyHigh()); + + // needed? FIXME + if(n_cell_type != 0) continue; + + if(Rechit.energy() > max[n_layer]) { + max[n_layer] = Rechit.energy(); + max_x[n_layer] = CellCentreXY.first; + max_y[n_layer] = CellCentreXY.second; + } + + if((Rechit.energy()) / ADCtoMIP[n_skiroc] <= CMTHRESHOLD) { + commonmode[n_skiroc] += Rechit.energy(); + cm_num[n_skiroc]++; + } + + }//Rechit loop ends here + // std::cout << " >>> found commonmode = " << commonmode << std::endl; + + + for(int iS = 0; iS < MAXSKIROCS; ++iS) { + commonmode[iS] = commonmode[iS] / cm_num[iS]; + // std::cout << " value = " << commonmode[iS] << std::endl; + } + + + edm::Handle Rechits1; + event.getByToken(HGCalTBRecHitCollection_, Rechits1); + + // looping over each rechit to fill histogram + double allcells_sum[MAXLAYERS], sevencells_sum[MAXLAYERS], nineteencells_sum[MAXLAYERS], radius[MAXLAYERS]; + double seedEnergy[MAXLAYERS]; + double x_tmp[MAXLAYERS], y_tmp[MAXLAYERS]; + int num[MAXLAYERS], sevennum[MAXLAYERS], nineteennum[MAXLAYERS]; + + for(int iL = 0; iL < MAXLAYERS; ++iL) { + allcells_sum[iL] = sevencells_sum[iL] = nineteencells_sum[iL] = radius[iL] = 0.; + x_tmp[iL] = y_tmp[iL] = seedEnergy[iL] = 0.; + num[iL] = sevennum[iL] = nineteennum[iL] = 0; + } + + for(auto Rechit1 : *Rechits1) { + + if(!IsCellValid.iu_iv_valid((Rechit1.id()).layer(), (Rechit1.id()).sensorIU(), (Rechit1.id()).sensorIV(), (Rechit1.id()).iu(), (Rechit1.id()).iv(), sensorsize)) continue; + + //getting electronics ID + uint32_t EID = essource_.emap_.detId2eid(Rechit1.id()); + HGCalTBElectronicsId eid(EID); + + int eLayer = (Rechit1.id()).layer() - 1; + int eCellType = (Rechit1.id()).cellType(); + int eSkiroc = eid.iskiroc() - 1; + + + //getting X and Y coordinates + CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots((Rechit1.id()).layer(), (Rechit1.id()).sensorIU(), (Rechit1.id()).sensorIV(), (Rechit1.id()).iu(), (Rechit1.id()).iv(), sensorsize); + + //FIXME >> needed? + if(eCellType != 0) continue; + + // std::cout << Rechit1.energy() << std::endl; + + radius[eLayer] = sqrt( pow(CellCentreXY.first - max_x[eLayer], 2) + pow(CellCentreXY.second - max_y[eLayer], 2) ); + + float energyCMsub = (Rechit1.energy() - commonmode[eSkiroc]) / ADCtoMIP[eSkiroc]; + // std::cout << "val = " << commonmode[eSkiroc] << std::endl; + if(energyCMsub > CMTHRESHOLD) { + allcells_sum[eLayer] += energyCMsub; + if(energyCMsub > seedEnergy[eLayer]) seedEnergy[eLayer] = energyCMsub; + } + + // countin all recHits even below CM threshold + num[eLayer]++; + + h_Radius[eLayer]->Fill(radius[eLayer]); + + //FIXME navigation to 7cells + if((radius[eLayer] < maxdist && sevennum[eLayer] < 7) && (energyCMsub > CMTHRESHOLD)) { + sevencells_sum[eLayer] += energyCMsub; + sevennum[eLayer]++; + } + + //FIXME navigation to 19cells + if((radius[eLayer] < 1.95 * maxdist && nineteennum[eLayer] < 19) && (energyCMsub > CMTHRESHOLD)) { // nineteencells_sum += (LayerWeight[LAYER]*(Rechit1.energyHigh() - commonmode))/ ADCtoMIP[LAYER]; - nineteencells_sum[eLayer] += energyCMsub; - x_tmp[eLayer] += CellCentreXY.first * energyCMsub; - y_tmp[eLayer] += CellCentreXY.second * energyCMsub; - nineteennum[eLayer]++; - } - } - - for(int iS=0; iS 32.) continue; + if(RecHit1.energyHigh() > 32.) continue; if((RecHit1.id()).cellType() == 0) { // Full_Cell[(RecHit1.id()).layer() - 1]->Fill(RecHit1.energyHigh()); Cell_counter[(RecHit1.id()).layer() - 1]++; diff --git a/Reco/plugins/RecHitPlotter_HighGain_New.cc b/Reco/plugins/RecHitPlotter_HighGain_New.cc index debde56..142df8d 100644 --- a/Reco/plugins/RecHitPlotter_HighGain_New.cc +++ b/Reco/plugins/RecHitPlotter_HighGain_New.cc @@ -317,6 +317,7 @@ RecHitPlotter_HighGain_New::analyze(const edm::Event& event, const edm::EventSet TH2Poly *h_RecHit_layer[MAXLAYERS]; int evId = event.id().event() - 1;// eventID starts from 1 + int EVENTSPERSPILL = 400; int iLayer = (evId % (MAXLAYERS * EVENTSPERSPILL)) / (EVENTSPERSPILL); cout << endl << " iLayer= " << iLayer << endl; h_RecHit_layer[iLayer] = fs->make(); diff --git a/test_cfg.py b/test_cfg.py index b203c0f..2fe6e83 100644 --- a/test_cfg.py +++ b/test_cfg.py @@ -147,12 +147,14 @@ for seq in options.chainSequence: if(seq == "DIGI"): process.p *= process.hgcaltbdigis - if(seq == "DIGIPLOT" or seq == "PED"): + if(seq == "DIGIDQM" or seq == "PED"): process.p *= process.hgcaltbdigisplotter if(seq == "PED"): process.hgcaltbdigisplotter.dumpNewPedestals = cms.untracked.bool(True) if(seq == "RECO"): process.p *= process.hgcaltbrechits + if(seq == "RECODQM"): + process.p *= process.hgcaltbrechitsplotter_highgain_new # if (options.chainSequence == 1): # process.p =cms.Path(process.hgcaltbdigis*process.hgcaltbdigisplotter) # elif (options.chainSequence == 3): From 49f7dc810fa7dbf78e230bcee6e13acbd39c4a26 Mon Sep 17 00:00:00 2001 From: Rajdeep Date: Wed, 19 Oct 2016 08:17:55 +0200 Subject: [PATCH 7/7] Working unpacker cross-checked with newTextInputFormat_Working + changes made by Arabella to seamlessly include HG->LG switching and different configurations: ###WARNING## DANGEROUSLY SLOW --- RawToDigi/plugins/HGCalTBRawToDigi.cc | 24 +- Reco/plugins/HGCalTBRecHitProducer.cc | 19 +- Reco/plugins/HGCalTBRecHitProducer.h | 6 +- Reco/plugins/Layer_Sum_Analyzer.cc | 648 ++++++++++++----------- Reco/python/hgcaltbrechitplotter_cfi.py | 4 +- Reco/python/hgcaltbrechitproducer_cfi.py | 4 +- test_cfg_newEB.py | 5 +- 7 files changed, 387 insertions(+), 323 deletions(-) diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.cc b/RawToDigi/plugins/HGCalTBRawToDigi.cc index 682503f..cb55b02 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.cc +++ b/RawToDigi/plugins/HGCalTBRawToDigi.cc @@ -11,6 +11,8 @@ #include #endif +using namespace std; + unsigned int gray_to_binary (unsigned int gray); HGCalTBRawToDigi::HGCalTBRawToDigi(edm::ParameterSet const& conf): @@ -49,23 +51,29 @@ void HGCalTBRawToDigi::produce(edm::Event& e, const edm::EventSetup& c) digis = std::auto_ptr(new SKIROC2DigiCollection(nSkirocs * SKIROC::NCHANNELS * SKIROC::MAXSAMPLES)); const uint16_t* pdata = (const uint16_t*)(fed.data()); - size_t ski = 0; // the skirocs have an absolute numbering, start counting from the first board till the last +// size_t ski = 0; // the skirocs have an absolute numbering, start counting from the first board till the last(Activate post consistent change in the EMAP) for (unsigned int i_board = 0 ; i_board < nBoards; ++i_board) { - for(size_t i_skiroc = 0; i_skiroc < MAXSKIROCS_PER_BOARD; ++i_skiroc) { - for (int ichan = 0; ichan < SKIROC::NCHANNELS; ichan++) { - HGCalTBElectronicsId eid(ski, ichan); +//The way the entries have been pushed they correspond to the Map a Ski 2,1,4,3,6,5...16,15 so for a simpler logic the EMAP will have to be modified that will be done later. +// for(size_t i_skiroc = 0; i_skiroc < MAXSKIROCS_PER_BOARD; ++i_skiroc) {//(Activate Post consistent change in EMAP) + for(size_t i_skiroc = MAXSKIROCS_PER_BOARD; i_skiroc >= 1; i_skiroc--) {//(De-Activate Post consistent change in EMAP) + for (int ichan = SKIROC::NCHANNELS - 1; ichan >= 0; ichan--) { +// HGCalTBElectronicsId eid(ski, ichan);//(Activate Post consistent change in EMAP) + HGCalTBElectronicsId eid(2*i_board + i_skiroc, ichan);//(De-Activate Post consistent change in EMAP) if (essource_.emap_.existsEId(eid.rawId())) { HGCalTBDetId did = essource_.emap_.eid2detId(eid); digis->addDataFrame(did); #ifdef DEBUG if(i_board == 0) std::cout << (*pdata & 0xFFF) << "\t" << (*(pdata + 1) & 0xFFF) << "\t" << (*(pdata + 2) & 0xFFF) << std::endl; +cout<backDataFrame().setSample(0, gray_to_binary(*(pdata) & 0xFFF), gray_to_binary( *(pdata + 1) & 0xFFF), 0); - pdata++; - pdata++; + + digis->backDataFrame().setSample(0, gray_to_binary(*(pdata + SKIROC::NCHANNELS) & 0xFFF), gray_to_binary( *(pdata) & 0xFFF), 0); } + pdata++;//Note this has to be outside the if condition, as for the test channel we wont enter the if condition but still have to increment the pointer address } - ++ski; //increment the absolute ID of the skiroc +//We have pushed back the high gain and low gain entries which are independent FED entries for all SKIROC::NCHANNELS so we need to jump by SKIROC::NCHANNELS + pdata = pdata + SKIROC::NCHANNELS; +// ++ski; //increment the absolute ID of the skiroc(Activate post consistent change in the EMAP) } } diff --git a/Reco/plugins/HGCalTBRecHitProducer.cc b/Reco/plugins/HGCalTBRecHitProducer.cc index 6921ec0..85ed06b 100644 --- a/Reco/plugins/HGCalTBRecHitProducer.cc +++ b/Reco/plugins/HGCalTBRecHitProducer.cc @@ -9,16 +9,26 @@ HGCalTBRecHitProducer::HGCalTBRecHitProducer(const edm::ParameterSet& cfg) _gainsLow_filename(cfg.getParameter("gainLow")), _gainsHigh_filename(cfg.getParameter("gainHigh")), _adcSaturation(cfg.getParameter("adcSaturation")), - _LG2HG_value(cfg.getParameter >("LG2HG_CERN")), - _mapFile(cfg.getParameter("mapFile")) + //_LG2HG_value(cfg.getParameter >("LG2HG_CERN")), + //_mapFile(cfg.getParameter("mapFile")), + _layers_config(cfg.getParameter("layers_config")) { produces (outputCollectionName); // std::cout << " >>> _LG2HG_value size = " << _LG2HG_value.size() << std::endl; + if(_layers_config == 0){ + _LG2HG_value = cfg.getParameter >("LG2HG_FNAL"); + _mapFile = cfg.getParameter("mapFile_FNAL"); + } + else{ + _LG2HG_value = cfg.getParameter >("LG2HG_CERN"); + _mapFile = cfg.getParameter("mapFile_CERN"); + } + HGCalCondObjectTextIO io(0); edm::FileInPath fip(_mapFile); - if (!io.load(fip.fullPath(), essource_.emap_)) { - throw cms::Exception("Unable to load electronics map"); + if (!io.load(fip.fullPath(), essource_.emap_)) { + throw cms::Exception("Unable to load electronics map"); }; } @@ -82,6 +92,7 @@ void HGCalTBRecHitProducer::produce(edm::Event& event, const edm::EventSetup& iS uint32_t EID = essource_.emap_.detId2eid(recHit.id()); HGCalTBElectronicsId eid(EID); + energy = ( energyHigh < _adcSaturation ) ? energyHigh : energyLow * _LG2HG_value.at(eid.iskiroc() - 1); recHit.setEnergy(energy); diff --git a/Reco/plugins/HGCalTBRecHitProducer.h b/Reco/plugins/HGCalTBRecHitProducer.h index 33e3771..0a6f672 100644 --- a/Reco/plugins/HGCalTBRecHitProducer.h +++ b/Reco/plugins/HGCalTBRecHitProducer.h @@ -2,7 +2,6 @@ #define HGCALTBRECHITPRODUCER_H /** \class Reco/plugins/HGCalTBRecHitProducer.h HGCalTBRecHitProducer HGCalTBRecHitProducer \brief - \author Shervin Nourbakhsh */ @@ -45,9 +44,10 @@ class HGCalTBRecHitProducer : public edm::EDProducer int _adcSaturation; std::vector _LG2HG_value; std::string _mapFile; + int _layers_config; struct { - HGCalElectronicsMap emap_; - } essource_; + HGCalElectronicsMap emap_; + } essource_; }; diff --git a/Reco/plugins/Layer_Sum_Analyzer.cc b/Reco/plugins/Layer_Sum_Analyzer.cc index 998ed2d..60eadca 100644 --- a/Reco/plugins/Layer_Sum_Analyzer.cc +++ b/Reco/plugins/Layer_Sum_Analyzer.cc @@ -60,21 +60,23 @@ double MIP2ParticleCalib = 1.06; // CERN to pion 125GeV //double LayerWeight[16] = {1.4091566745180932, 0.7020676448403224, 0.6054055986179145, 0.7020676448403224, 0.6054055986179145, 0.8415931435769973, 0.8061197656138868, 0.9811186423136724, 0.8061197656138868, 0.9811186423136724, 0.8061197656138868, 1.4462036381025891, 1.2075480996058319, 1.4462036381025891, 1.2075480996058319, 1.4462036381025891}; -double LayerWeight_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; -double X0depth_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; +double LayerWeight_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; +double X0depth_16L_FNAL[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; double LayerWeight_8L_conf1[16] = {33.074, 13.184, 14.17, 9.788, 9.766, 9.766, 16.339, 14.129, 0., 0., 0., 0., 0., 0., 0., 0.}; double X0depth_8L_conf1[16] = {6.268, 1.131, 1.131, 1.362, 0.574, 1.301, 0.574, 2.42, 0., 0., 0., 0., 0., 0., 0., 0.}; double LayerWeight_8L_conf2[16] = {35.866, 30.864, 28.803, 23.095, 20.657, 19.804, 36.322, 27.451, 0., 0., 0., 0., 0., 0., 0., 0.}; double X0depth_8L_conf2[16] = {5.048, 3.412, 3.412, 2.866, 2.512, 1.625, 2.368, 6.021, 0., 0., 0., 0., 0., 0., 0., 0.}; double weights2GeV = 1.e-03; -double weights2MIP = 52.8 / 63.6; // rescale weights from mean to MPV +double MIP2GeV_sim = 52.81e-06; +double weights2MIP = 52.8/63.6; // rescale weights from mean to MPV //double LayerWeight[16] = {0.4847555727337982, 1.0214605968539232, 0.4847555727337982, 1.0214605968539232, 0.4847555727337982, 1.1420105918768606, 0.6423912113800805, 1.2625605868997982, 0.6423912113800805, 1.2625605868997982, 0.6423912113800805, 1.6643939036429232, 0.9576624886726451, 1.6643939036429232, 0.9576624886726451, 1.6643939036429232};// dE/dx weights //double LayerSumWeight = 1.; -const int CMTHRESHOLD = 30;// anything less than this value is added to the commonmode sum +//const int CMTHRESHOLD = 30;// anything less than this value is added to the commonmode sum +const int CMTHRESHOLD = 2;// anything less than this value is added to the commonmode sum // applied to all layers sum after commonmode subtraction and the ADC to MIP conversion const double ALLCELLS_THRESHOLD = 50.; @@ -104,42 +106,44 @@ class Layer_Sum_Analyzer : public edm::one::EDAnalyzer> CellXY; std::pair CellCentreXY; - HGCalTBTopology IsCellValid; + HGCalTBTopology IsCellValid; HGCalTBCellVertices TheCell; - double maxdist = (1 + sqrt (3) / 2) * HGCAL_TB_CELL::FULL_CELL_SIDE; // <<< FIXME maxdist > HGCAL_TB_CELL::FULL_CELL_SIDE !! + double maxdist = (1 + sqrt (3) / 2) * HGCAL_TB_CELL::FULL_CELL_SIDE; // <<< FIXME maxdist > HGCAL_TB_CELL::FULL_CELL_SIDE !! - double Weights_L[MAXLAYERS]; - double X0_L[MAXLAYERS]; - double ADCtoMIP[MAXLAYERS]; + double Weights_L[MAXLAYERS]; + double X0_L[MAXLAYERS]; + double ADCtoMIP[MAXLAYERS]; - TH1F *h_CM_layer[MAXSKIROCS]; + TH1F *h_CM_layer[MAXSKIROCS]; TH1F *h_sum_layer[MAXLAYERS], *h_layer_seven[MAXLAYERS], *h_layer_nineteen[MAXLAYERS], *h_Seed_layer[MAXLAYERS]; - TH1F *h_sum_layer_AbsW[MAXLAYERS], *h_layer_seven_AbsW[MAXLAYERS], *h_layer_nineteen_AbsW[MAXLAYERS], *h_Seed_layer_AbsW[MAXLAYERS]; - TH1F *h_x_layer[MAXLAYERS], *h_y_layer[MAXLAYERS]; + TH1F *h_sum_layer_AbsW_Mip[MAXLAYERS], *h_layer_seven_AbsW_Mip[MAXLAYERS], *h_layer_nineteen_AbsW_Mip[MAXLAYERS], *h_Seed_layer_AbsW_Mip[MAXLAYERS]; + TH1F *h_sum_layer_AbsW_GeV[MAXLAYERS], *h_layer_seven_AbsW_GeV[MAXLAYERS], *h_layer_nineteen_AbsW_GeV[MAXLAYERS], *h_Seed_layer_AbsW_GeV[MAXLAYERS]; + TH1F *h_x_layer[MAXLAYERS], *h_y_layer[MAXLAYERS]; TH2F *h_x_y_layer[MAXLAYERS]; - TH1F* h_Radius[MAXLAYERS]; - TH1F *h_E1oE7_layer[MAXLAYERS], *h_E1oE19_layer[MAXLAYERS], *h_E7oE19_layer[MAXLAYERS]; + TH1F* h_Radius[MAXLAYERS]; + TH1F *h_E1oE7_layer[MAXLAYERS], *h_E1oE19_layer[MAXLAYERS], *h_E7oE19_layer[MAXLAYERS]; - TH1F *h_sum_all, *h_seven_all, *h_nineteen_all; - TH1F *h_sum_all_AbsW, *h_seven_all_AbsW, *h_nineteen_all_AbsW; - TProfile *tp_E1_vs_layer, *tp_E7_vs_layer, *tp_E19_vs_layer; - TProfile *tp_E1oSumL_vs_layer, *tp_E7oSumL_vs_layer, *tp_E19oSumL_vs_layer; - TProfile *tp_E1oE7_vs_layer, *tp_E1oE19_vs_layer, *tp_E7oE19_vs_layer; + TH1F *h_sum_all, *h_seven_all, *h_nineteen_all; + TH1F *h_sum_all_AbsW_Mip, *h_seven_all_AbsW_Mip, *h_nineteen_all_AbsW_Mip; + TH1F *h_sum_all_AbsW_GeV, *h_seven_all_AbsW_GeV, *h_nineteen_all_AbsW_GeV; + TProfile *tp_E1_vs_layer, *tp_E7_vs_layer, *tp_E19_vs_layer; + TProfile *tp_E1oSumL_vs_layer, *tp_E7oSumL_vs_layer, *tp_E19oSumL_vs_layer; + TProfile *tp_E1oE7_vs_layer, *tp_E1oE19_vs_layer, *tp_E7oE19_vs_layer; TH2F *HighGain_LowGain_2D; TH2F *Energy_LowGain_2D; - int EVENT = 0; + int EVENT = 0; map Time_Stamp; map Delta_Time_Stamp; @@ -150,7 +154,7 @@ class Layer_Sum_Analyzer : public edm::one::EDAnalyzer fs; @@ -159,28 +163,33 @@ Layer_Sum_Analyzer::Layer_Sum_Analyzer(const edm::ParameterSet& iConfig) //booking the histos for(int layer = 0; layer < MAXLAYERS; layer++) { - h_sum_layer[layer] = fs->make(Form("sumAll_Layer%d", layer + 1), "", 40010, -10, 40000); - h_Seed_layer[layer] = fs->make(Form("h_Seed_layer%d", layer + 1), Form("h_Seed_layer%d", layer + 1), 40010, -10, 40000); - h_layer_seven[layer] = fs->make(Form("sum7_Layer%d", layer + 1), "", 40010, -10, 40000); - h_layer_nineteen[layer] = fs->make(Form("sum19_Layer%d", layer + 1), "", 40010, -10, 40000); - - h_sum_layer_AbsW[layer] = fs->make(Form("sumAll_Layer%d_AbsW", layer + 1), "", 40010, -10, 40000); - h_Seed_layer_AbsW[layer] = fs->make(Form("h_Seed_layer%d_AbsW", layer + 1), Form("h_Seed_layer%d", layer + 1), 40010, -10, 40000); - h_layer_seven_AbsW[layer] = fs->make(Form("sum7_Layer%d_AbsW", layer + 1), "", 40010, -10, 40000); - h_layer_nineteen_AbsW[layer] = fs->make(Form("sum19_Layer%d_AbsW", layer + 1), "", 40010, -10, 40000); - - h_x_layer[layer] = fs->make(Form("X_Layer%d", layer + 1), "", 2000, -10., 10. ); - h_y_layer[layer] = fs->make(Form("Y_Layer%d", layer + 1), "", 2000, -10., 10. ); - h_x_y_layer[layer] = fs->make(Form("YvsX_Layer%d", layer + 1), "", 2000, -10., 10., 2000, -10., 10. ); - - h_E1oE7_layer[layer] = fs->make(Form("h_E1oE7_layer%d", layer + 1), Form("h_E1oE7_layer%d", layer + 1), 5000, -5, 5); - h_E1oE19_layer[layer] = fs->make(Form("h_E1oE19_layer%d", layer + 1), Form("h_E1oE19_layer%d", layer + 1), 5000, -5, 5); - h_E7oE19_layer[layer] = fs->make(Form("h_E7oE19_layer%d", layer + 1), Form("h_E7oE19_layer%d", layer + 1), 5000, -5, 5); - h_Radius[layer] = fs->make(Form("h_Radius_layer%d", layer + 1), Form("h_Radius_layer%d", layer + 1), 5000, 0., 20.); + h_sum_layer[layer] = fs->make(Form("sumAll_Layer%d", layer+1), "", 40010, -10, 40000); + h_Seed_layer[layer] = fs->make(Form("h_Seed_layer%d",layer+1), Form("h_Seed_layer%d",layer+1), 40010, -10, 40000); + h_layer_seven[layer] = fs->make(Form("sum7_Layer%d", layer+1),"", 40010, -10, 40000); + h_layer_nineteen[layer] = fs->make(Form("sum19_Layer%d", layer+1), "", 40010, -10, 40000); + + h_sum_layer_AbsW_Mip[layer] = fs->make(Form("sumAll_Layer%d_AbsW_Mip", layer+1), "", 40010, -10, 40000); + h_Seed_layer_AbsW_Mip[layer] = fs->make(Form("h_Seed_layer%d_AbsW_Mip",layer+1), Form("h_Seed_layer%d",layer+1), 40010, -10, 40000); + h_layer_seven_AbsW_Mip[layer] = fs->make(Form("sum7_Layer%d_AbsW_Mip", layer+1),"", 40010, -10, 40000); + h_layer_nineteen_AbsW_Mip[layer] = fs->make(Form("sum19_Layer%d_AbsW_Mip", layer+1), "", 40010, -10, 40000); + + h_sum_layer_AbsW_GeV[layer] = fs->make(Form("sumAll_Layer%d_AbsW_GeV", layer+1), "", 40010, -10, 40000); + h_Seed_layer_AbsW_GeV[layer] = fs->make(Form("h_Seed_layer%d_AbsW_GeV",layer+1), Form("h_Seed_layer%d",layer+1), 40010, -10, 40000); + h_layer_seven_AbsW_GeV[layer] = fs->make(Form("sum7_Layer%d_AbsW_GeV", layer+1),"", 40010, -10, 40000); + h_layer_nineteen_AbsW_GeV[layer] = fs->make(Form("sum19_Layer%d_AbsW_GeV", layer+1), "", 40010, -10, 40000); + + h_x_layer[layer] = fs->make(Form("X_Layer%d", layer+1), "", 2000, -10., 10. ); + h_y_layer[layer] = fs->make(Form("Y_Layer%d", layer+1), "", 2000, -10., 10. ); + h_x_y_layer[layer] = fs->make(Form("YvsX_Layer%d", layer+1), "", 2000, -10., 10., 2000, -10., 10. ); + + h_E1oE7_layer[layer] = fs->make(Form("h_E1oE7_layer%d",layer+1), Form("h_E1oE7_layer%d",layer+1), 5000, -5, 5); + h_E1oE19_layer[layer] = fs->make(Form("h_E1oE19_layer%d",layer+1), Form("h_E1oE19_layer%d",layer+1), 5000, -5, 5); + h_E7oE19_layer[layer] = fs->make(Form("h_E7oE19_layer%d",layer+1), Form("h_E7oE19_layer%d",layer+1), 5000, -5, 5); + h_Radius[layer] = fs->make(Form("h_Radius_layer%d",layer+1), Form("h_Radius_layer%d",layer+1), 5000, 0., 20.); } - for(int ski = 0; ski < MAXSKIROCS; ++ski) { - h_CM_layer[ski] = fs->make(Form("h_CM_skiroc%d", ski + 1), "", 5000, -500, 500); + for(int ski=0; skimake(Form("h_CM_skiroc%d",ski+1), "", 5000, -500, 500); } h_sum_all = fs->make("h_sumAll_AllLayers", "", 40010, -10, 40000); @@ -189,56 +198,68 @@ Layer_Sum_Analyzer::Layer_Sum_Analyzer(const edm::ParameterSet& iConfig) h_sum_all->Sumw2(); h_seven_all->Sumw2(); h_nineteen_all->Sumw2(); - h_sum_all_AbsW = fs->make("h_sumAll_AllLayers_AbsW", "", 40010, -10, 40000); - h_seven_all_AbsW = fs->make("h_sum7_AllLayers_AbsW", "", 40010, -10, 40000); - h_nineteen_all_AbsW = fs->make("h_sum19_AllLayers_AbsW", "", 40010, -10, 40000); - h_sum_all_AbsW->Sumw2(); - h_seven_all_AbsW->Sumw2(); - h_nineteen_all_AbsW->Sumw2(); + h_sum_all_AbsW_Mip = fs->make("h_sumAll_AllLayers_AbsW_Mip", "", 40010, -10, 40000); + h_seven_all_AbsW_Mip = fs->make("h_sum7_AllLayers_AbsW_Mip", "", 40010, -10, 40000); + h_nineteen_all_AbsW_Mip = fs->make("h_sum19_AllLayers_AbsW_Mip", "", 40010, -10, 40000); + h_sum_all_AbsW_Mip->Sumw2(); + h_seven_all_AbsW_Mip->Sumw2(); + h_nineteen_all_AbsW_Mip->Sumw2(); + + h_sum_all_AbsW_GeV = fs->make("h_sumAll_AllLayers_AbsW_GeV", "", 40010, -10, 40000); + h_seven_all_AbsW_GeV = fs->make("h_sum7_AllLayers_AbsW_GeV", "", 40010, -10, 40000); + h_nineteen_all_AbsW_GeV = fs->make("h_sum19_AllLayers_AbsW_GeV", "", 40010, -10, 40000); + h_sum_all_AbsW_GeV->Sumw2(); + h_seven_all_AbsW_GeV->Sumw2(); + h_nineteen_all_AbsW_GeV->Sumw2(); HighGain_LowGain_2D = fs->make("h2_HGvsLG", "", 4000, 0, 4000, 4000, 0, 4000); Energy_LowGain_2D = fs->make("h2_EvsLG", "", 4000, 0, 4000, 4000, 0, 4000); - tp_E1_vs_layer = fs->make("tp_E1_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E7_vs_layer = fs->make("tp_E7_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E19_vs_layer = fs->make("tp_E19_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E1oSumL_vs_layer = fs->make("tp_E1oSumL_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E7oSumL_vs_layer = fs->make("tp_E7oSumL_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E19oSumL_vs_layer = fs->make("tp_E19oSumL_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E1oE7_vs_layer = fs->make("tp_E1oE7_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E1oE19_vs_layer = fs->make("tp_E1oE19_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - tp_E7oE19_vs_layer = fs->make("tp_E7oE19_vs_layer", "", MAXLAYERS + 1, 0, MAXLAYERS + 1); - - - //loading the proper weights - if(MAXLAYERS != 8) { - std::cout << " update weights " << std::endl; - return; + tp_E1_vs_layer = fs->make("tp_E1_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E7_vs_layer = fs->make("tp_E7_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E19_vs_layer = fs->make("tp_E19_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E1oSumL_vs_layer = fs->make("tp_E1oSumL_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E7oSumL_vs_layer = fs->make("tp_E7oSumL_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E19oSumL_vs_layer = fs->make("tp_E19oSumL_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E1oE7_vs_layer = fs->make("tp_E1oE7_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E1oE19_vs_layer = fs->make("tp_E1oE19_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + tp_E7oE19_vs_layer = fs->make("tp_E7oE19_vs_layer", "", MAXLAYERS+1, 0, MAXLAYERS+1); + + + //loading the proper weights + if(MAXLAYERS != 8) { + std::cout << " update weights " << std::endl; + return; } - for(int iL = 0; iL < MAXLAYERS; ++iL) { - if(layers_config_ == 0) { - Weights_L[iL] = LayerWeight_16L_FNAL[iL]; - X0_L[iL] = X0depth_16L_FNAL[iL]; - ADCtoMIP[iL] = ADCtoMIP_FNAL[iL]; - } else { - ADCtoMIP[iL] = ADCtoMIP_CERN[iL]; - if(layers_config_ == 1) { - Weights_L[iL] = LayerWeight_8L_conf1[iL]; - X0_L[iL] = X0depth_8L_conf1[iL]; - } - if(layers_config_ == 2) { - Weights_L[iL] = LayerWeight_8L_conf2[iL]; - X0_L[iL] = X0depth_8L_conf2[iL]; - } - if(layers_config_ == -1) { - Weights_L[iL] = 1.; - X0_L[iL] = 0.; - } - } + for(int iL=0; iL("mapFile_FNAL"); + } + else{ + ADCtoMIP[iL] = ADCtoMIP_CERN[iL]; + mapfile_ = iConfig.getParameter("mapFile_CERN"); + if(layers_config_ == 1){ + Weights_L[iL] = LayerWeight_8L_conf1[iL]; + X0_L[iL] = X0depth_8L_conf1[iL]; + } + if(layers_config_ == 2){ + Weights_L[iL] = LayerWeight_8L_conf2[iL]; + X0_L[iL] = X0depth_8L_conf2[iL]; + } + if(layers_config_ == -1){ + Weights_L[iL] = 1.; + X0_L[iL] = 0.; + } + } } + + }//constructor ends here @@ -256,224 +277,245 @@ void Layer_Sum_Analyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) { - // std::cout << " >>> Layer_Sum_Analyzer::analyze " << std::endl; - EVENT = (event.id()).event(); - - if((event.id()).event() == 1) { - Time_Temp = event.time().value(); - Time_Stamp[EVENT] = Time_Temp; - Delta_Time_Stamp[EVENT] = 0.; - } else { - Time_Stamp[EVENT] = event.time().value(); - Delta_Time_Stamp[EVENT] = event.time().value() - Time_Temp; - Time_Temp = event.time().value(); - } - - //opening Rechits - edm::Handle Rechits; - event.getByToken(HGCalTBRecHitCollection_, Rechits); - - // looping over each rechit to fill histogram - double commonmode[MAXSKIROCS]; - int cm_num[MAXSKIROCS]; - for(int iS = 0; iS < MAXSKIROCS; ++iS) { - commonmode[iS] = 0.; - cm_num[iS] = 0; - } - double max[MAXLAYERS], max_x[MAXLAYERS], max_y[MAXLAYERS]; - for(int iL = 0; iL < MAXLAYERS; ++iL) { - max[iL] = max_x[iL] = max_y[iL] = 0.; - } - - - //CM subtraction => REPLACE WITH EXTERNAL CLASS => FIX - for(auto Rechit : *Rechits) { - - if(!IsCellValid.iu_iv_valid((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize)) continue; - - //getting electronics ID - uint32_t EID = essource_.emap_.detId2eid(Rechit.id()); - HGCalTBElectronicsId eid(EID); - - int n_layer = (Rechit.id()).layer() - 1; - int n_cell_type = (Rechit.id()).cellType(); - int n_skiroc = eid.iskiroc() - 1; - - - //getting X and Y coordinates - CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize); - - HighGain_LowGain_2D->Fill(Rechit.energyLow(), Rechit.energyHigh()); - Energy_LowGain_2D->Fill(Rechit.energy(), Rechit.energyHigh()); - - // needed? FIXME - if(n_cell_type != 0) continue; - - if(Rechit.energy() > max[n_layer]) { - max[n_layer] = Rechit.energy(); - max_x[n_layer] = CellCentreXY.first; - max_y[n_layer] = CellCentreXY.second; - } - - if((Rechit.energy()) / ADCtoMIP[n_skiroc] <= CMTHRESHOLD) { - commonmode[n_skiroc] += Rechit.energy(); - cm_num[n_skiroc]++; - } - - }//Rechit loop ends here - // std::cout << " >>> found commonmode = " << commonmode << std::endl; - - - for(int iS = 0; iS < MAXSKIROCS; ++iS) { - commonmode[iS] = commonmode[iS] / cm_num[iS]; - // std::cout << " value = " << commonmode[iS] << std::endl; - } - - - edm::Handle Rechits1; - event.getByToken(HGCalTBRecHitCollection_, Rechits1); - - // looping over each rechit to fill histogram - double allcells_sum[MAXLAYERS], sevencells_sum[MAXLAYERS], nineteencells_sum[MAXLAYERS], radius[MAXLAYERS]; - double seedEnergy[MAXLAYERS]; - double x_tmp[MAXLAYERS], y_tmp[MAXLAYERS]; - int num[MAXLAYERS], sevennum[MAXLAYERS], nineteennum[MAXLAYERS]; - - for(int iL = 0; iL < MAXLAYERS; ++iL) { - allcells_sum[iL] = sevencells_sum[iL] = nineteencells_sum[iL] = radius[iL] = 0.; - x_tmp[iL] = y_tmp[iL] = seedEnergy[iL] = 0.; - num[iL] = sevennum[iL] = nineteennum[iL] = 0; - } - - for(auto Rechit1 : *Rechits1) { - - if(!IsCellValid.iu_iv_valid((Rechit1.id()).layer(), (Rechit1.id()).sensorIU(), (Rechit1.id()).sensorIV(), (Rechit1.id()).iu(), (Rechit1.id()).iv(), sensorsize)) continue; - - //getting electronics ID - uint32_t EID = essource_.emap_.detId2eid(Rechit1.id()); - HGCalTBElectronicsId eid(EID); - - int eLayer = (Rechit1.id()).layer() - 1; - int eCellType = (Rechit1.id()).cellType(); - int eSkiroc = eid.iskiroc() - 1; - - - //getting X and Y coordinates - CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots((Rechit1.id()).layer(), (Rechit1.id()).sensorIU(), (Rechit1.id()).sensorIV(), (Rechit1.id()).iu(), (Rechit1.id()).iv(), sensorsize); - - //FIXME >> needed? - if(eCellType != 0) continue; - - // std::cout << Rechit1.energy() << std::endl; - - radius[eLayer] = sqrt( pow(CellCentreXY.first - max_x[eLayer], 2) + pow(CellCentreXY.second - max_y[eLayer], 2) ); - - float energyCMsub = (Rechit1.energy() - commonmode[eSkiroc]) / ADCtoMIP[eSkiroc]; - // std::cout << "val = " << commonmode[eSkiroc] << std::endl; - if(energyCMsub > CMTHRESHOLD) { - allcells_sum[eLayer] += energyCMsub; - if(energyCMsub > seedEnergy[eLayer]) seedEnergy[eLayer] = energyCMsub; - } - - // countin all recHits even below CM threshold - num[eLayer]++; - - h_Radius[eLayer]->Fill(radius[eLayer]); - - //FIXME navigation to 7cells - if((radius[eLayer] < maxdist && sevennum[eLayer] < 7) && (energyCMsub > CMTHRESHOLD)) { - sevencells_sum[eLayer] += energyCMsub; - sevennum[eLayer]++; - } - - //FIXME navigation to 19cells - if((radius[eLayer] < 1.95 * maxdist && nineteennum[eLayer] < 19) && (energyCMsub > CMTHRESHOLD)) { + // std::cout << " >>> Layer_Sum_Analyzer::analyze " << std::endl; + EVENT = (event.id()).event(); + + if((event.id()).event() == 1) { + Time_Temp = event.time().value(); + Time_Stamp[EVENT] = Time_Temp; + Delta_Time_Stamp[EVENT] = 0.; + } else { + Time_Stamp[EVENT] = event.time().value(); + Delta_Time_Stamp[EVENT] = event.time().value() - Time_Temp; + Time_Temp = event.time().value(); + } + + //opening Rechits + edm::Handle Rechits; + event.getByToken(HGCalTBRecHitCollection_, Rechits); + + // looping over each rechit to fill histogram + double commonmode[MAXSKIROCS]; + int cm_num[MAXSKIROCS]; + for(int iS=0; iS REPLACE WITH EXTERNAL CLASS => FIX + for(auto Rechit : *Rechits){ + + if(!IsCellValid.iu_iv_valid((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize)) continue; + + //getting electronics ID + uint32_t EID = essource_.emap_.detId2eid(Rechit.id()); + HGCalTBElectronicsId eid(EID); + + int n_layer = (Rechit.id()).layer() - 1; + int n_cell_type = (Rechit.id()).cellType(); + int n_skiroc = eid.iskiroc() - 1; + + + //getting X and Y coordinates + CellCentreXY = TheCell.GetCellCentreCoordinatesForPlots((Rechit.id()).layer(), (Rechit.id()).sensorIU(), (Rechit.id()).sensorIV(), (Rechit.id()).iu(), (Rechit.id()).iv(), sensorsize); + + HighGain_LowGain_2D->Fill(Rechit.energyLow(), Rechit.energyHigh()); + Energy_LowGain_2D->Fill(Rechit.energy(), Rechit.energyHigh()); + + // needed? FIXME + if(n_cell_type != 0 && n_cell_type != 4) continue; + + if(Rechit.energy() > max[n_layer]) { + max[n_layer] = Rechit.energy(); + max_x[n_layer] = CellCentreXY.first; + max_y[n_layer] = CellCentreXY.second; + } + + if((Rechit.energy()) / ADCtoMIP[n_layer] <= CMTHRESHOLD) { + commonmode[n_skiroc] += Rechit.energy(); + cm_num[n_skiroc]++; + } + + }//Rechit loop ends here + // std::cout << " >>> found commonmode = " << commonmode << std::endl; + + + for(int iS=0; iS