diff --git a/DataFormats/interface/SKIROCParameters.h b/DataFormats/interface/SKIROCParameters.h index 3d5e2a9..02d8fbe 100644 --- a/DataFormats/interface/SKIROCParameters.h +++ b/DataFormats/interface/SKIROCParameters.h @@ -5,7 +5,6 @@ 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 } 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..ad5c1bf 100644 --- a/Geometry/interface/HGCalTBCellVertices.h +++ b/Geometry/interface/HGCalTBCellVertices.h @@ -26,6 +26,12 @@ class HGCalTBCellVertices 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) { return GetCellCoordinates(layer, sensor_iu, sensor_iv, iu, iv, sensorsize, true); @@ -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/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/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..db5ddcf 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] + Description: Plugin to make 2D and 1D histograms of digis - Implementation: - [Notes on implementation] + 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,22 @@ 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; + bool _dumpNewPedestals; }; // @@ -108,83 +103,75 @@ 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"); 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 + }//loop over i_samples - m_pedestalsHighGain = iConfig.getUntrackedParameter("pedestalsHighGain", ""); - m_pedestalsLowGain = iConfig.getUntrackedParameter("pedestalsLowGain", ""); + m_pedestalsHighGain = iConfig.getUntrackedParameter("pedestalsHighGain", ""); + m_pedestalsLowGain = iConfig.getUntrackedParameter("pedestalsLowGain", ""); + _dumpNewPedestals = iConfig.getUntrackedParameter("dumpNewPedestals", false); }//contructor ends here @@ -228,16 +215,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,32 +228,26 @@ 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); 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; + 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; + } } } } @@ -330,8 +311,9 @@ 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", ""); + desc.addUntracked("dumpNewPedestals", false); descriptions.add("hgcaltbdigisplotter", desc); } diff --git a/RawToDigi/plugins/HGCalTBRawToDigi.cc b/RawToDigi/plugins/HGCalTBRawToDigi.cc index c995dfb..cb55b02 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.cc +++ b/RawToDigi/plugins/HGCalTBRawToDigi.cc @@ -1,10 +1,20 @@ -#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" + +// provide the number of channels in a skiroc and number of samples +#include "HGCal/DataFormats/interface/SKIROCParameters.h" + +//#define DEBUG +#ifdef DEBUG +#include +#endif + 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 +46,37 @@ 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 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) { +//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 + SKIROC::NCHANNELS) & 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..6961b30 100644 --- a/RawToDigi/plugins/HGCalTBRawToDigi.h +++ b/RawToDigi/plugins/HGCalTBRawToDigi.h @@ -12,13 +12,15 @@ #include "HGCal/DataFormats/interface/HGCalTBDataFrameContainers.h" -#include "HGCal/DataFormats/interface/SKIROCParameters.h" #include /** * \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 bd5b395..660f5d2 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.cc +++ b/RawToDigi/plugins/HGCalTBTextSource.cc @@ -4,8 +4,8 @@ #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 /** 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 @@ -40,8 +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; - if(m_spill <= NSpills) return true; - else return false; + return true; } // sets m_run, m_spill, m_even, m_time @@ -74,41 +73,45 @@ bool HGCalTBTextSource::readLines() { max_boards = 0; for( auto& board : m_lines) { - board.clear(); + for( auto& skiroc : board) { + skiroc.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); + 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 && 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() != '\n' && 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); + } ++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(); + return true; } +/** 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,11 +119,13 @@ 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) { - auto board = m_lines[i_board]; - for (auto skiword : board) { - skiwords.push_back(skiword >> 16); - skiwords.push_back(skiword); + 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) { + skiwords.push_back(skiword >> 16); + skiwords.push_back(skiword); + } } } @@ -139,7 +144,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 afc8fda..cf29903 100644 --- a/RawToDigi/plugins/HGCalTBTextSource.h +++ b/RawToDigi/plugins/HGCalTBTextSource.h @@ -8,24 +8,79 @@ #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 + * \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 (UMN) + * \author Rajdeep Mohan Chatterjee (UMN) + * \author Jeremy Mans (UMN) + * + * \todo fedID is taken from cfg, it should be taken from the txt file + * + * 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 { 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_file(0) { - - m_sourceId = pset.getUntrackedParameter("fed", 1000); /// \todo check and read from file? + m_sourceId = pset.getUntrackedParameter("fed", 1000); /// \todo read from file produces(); } @@ -39,14 +94,18 @@ 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> , MAXSKIROCS_PER_BOARD >, MAXLAYERS > m_lines; 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/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/Reco/plugins/HGCalTBRecHitProducer.cc b/Reco/plugins/HGCalTBRecHitProducer.cc index 5f65c58..85ed06b 100644 --- a/Reco/plugins/HGCalTBRecHitProducer.cc +++ b/Reco/plugins/HGCalTBRecHitProducer.cc @@ -9,12 +9,22 @@ 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_)) { @@ -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 961df11..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,6 +44,7 @@ class HGCalTBRecHitProducer : public edm::EDProducer int _adcSaturation; std::vector _LG2HG_value; std::string _mapFile; + int _layers_config; struct { HGCalElectronicsMap emap_; } essource_; diff --git a/Reco/plugins/Layer_Sum_Analyzer.cc b/Reco/plugins/Layer_Sum_Analyzer.cc index 5ea551d..60eadca 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 @@ -69,13 +68,15 @@ double X0depth_8L_conf1[16] = {6.268, 1.131, 1.131, 1.362, 0.574, 1.301, 0.574, 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 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.; @@ -110,7 +111,7 @@ class Layer_Sum_Analyzer : public edm::one::EDAnalyzer> CellXY; std::pair CellCentreXY; @@ -124,7 +125,8 @@ class Layer_Sum_Analyzer : public edm::one::EDAnalyzermake(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_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. ); @@ -190,12 +198,19 @@ 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); @@ -222,9 +237,11 @@ Layer_Sum_Analyzer::Layer_Sum_Analyzer(const edm::ParameterSet& iConfig) Weights_L[iL] = LayerWeight_16L_FNAL[iL]; X0_L[iL] = X0depth_16L_FNAL[iL]; ADCtoMIP[iL] = ADCtoMIP_FNAL[iL]; + mapfile_ = iConfig.getParameter("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]; @@ -241,6 +258,8 @@ Layer_Sum_Analyzer::Layer_Sum_Analyzer(const edm::ParameterSet& iConfig) } + + }//constructor ends here @@ -309,7 +328,7 @@ Layer_Sum_Analyzer::analyze(const edm::Event& event, const edm::EventSetup& setu Energy_LowGain_2D->Fill(Rechit.energy(), Rechit.energyHigh()); // needed? FIXME - if(n_cell_type != 0) continue; + if(n_cell_type != 0 && n_cell_type != 4) continue; if(Rechit.energy() > max[n_layer]) { max[n_layer] = Rechit.energy(); @@ -317,7 +336,7 @@ Layer_Sum_Analyzer::analyze(const edm::Event& event, const edm::EventSetup& setu max_y[n_layer] = CellCentreXY.second; } - if((Rechit.energy()) / ADCtoMIP[n_skiroc] <= CMTHRESHOLD) { + if((Rechit.energy()) / ADCtoMIP[n_layer] <= CMTHRESHOLD) { commonmode[n_skiroc] += Rechit.energy(); cm_num[n_skiroc]++; } @@ -326,9 +345,10 @@ Layer_Sum_Analyzer::analyze(const edm::Event& event, const edm::EventSetup& setu // std::cout << " >>> found commonmode = " << commonmode << std::endl; - for(int iS=0; iS