CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalGainsCheck.cc
Go to the documentation of this file.
2 
4 {
5  rootfile = ps.getUntrackedParameter<std::string>("rootfile","null");
6  outfile = ps.getUntrackedParameter<std::string>("outFile","null");
7  dumpupdate = ps.getUntrackedParameter<std::string>("dumpUpdateGainsTo","null");
8  dumprefs = ps.getUntrackedParameter<std::string>("dumpRefGainsTo","null");
9  emapflag = ps.getUntrackedParameter<bool>("checkEmap",false);
10  validategainsflag = ps.getUntrackedParameter<bool>("validateGains",false);
11  epsilon = ps.getUntrackedParameter<double>("deltaG",1000000);
12 }
13 
15 {
16  f = new TFile(rootfile.c_str(),"RECREATE");
17 
18  //book histos:
19  ocMapUp = new TH2F("ocMapUp","occupancy_map_updated_gains",83,-41.5,41.5,72,0.5,72.5);
20  ocMapRef = new TH2F("ocMapRef","occupancy_map_updated_gains",83,-41.5,41.5,72,0.5,72.5);
21 // valMapUp;
22 // valMapRef;
23 
24  diffUpRefCap0 = new TH1F("diffUpRefCap0","difference_update_reference_Cap0",100,-0.5,0.5);
25  ratioUpRefCap0 = new TH1F("ratioUpRefCap0", "ration_update_reference_Cap0",100,0.5,1.5);
26  gainsUpCap0 = new TH1F("gainsUpCap0","gains_update_Cap0",100,0.0,0.6);
27  gainsRefCap0 = new TH1F("gainsRefCap0","gains_reference_Cap0",100,0.0,0.6);
28 
29  diffUpRefCap1 = new TH1F("diffUpRefCap1","difference_update_reference_Cap1",100,-0.5,0.5);
30  ratioUpRefCap1 = new TH1F("ratioUpRefCap1", "ration_update_reference_Cap1",100,0.5,1.5);
31  gainsUpCap1 = new TH1F("gainsUpCap1","gains_update_Cap1",100,0.0,0.6);
32  gainsRefCap1 = new TH1F("gainsRefCap1","gains_reference_Cap1",100,0.0,0.6);
33 
34  diffUpRefCap2 = new TH1F("diffUpRefCap2","difference_update_reference_Cap2",100,-0.5,0.5);
35  ratioUpRefCap2 = new TH1F("ratioUpRefCap2", "ration_update_reference_Cap2",100,0.5,1.5);
36  gainsUpCap2 = new TH1F("gainsUpCap2","gains_update_Cap2",100,0.0,0.6);
37  gainsRefCap2 = new TH1F("gainsRefCap2","gains_reference_Cap2",100,0.0,0.6);
38 
39  diffUpRefCap3 = new TH1F("diffUpRefCap3","difference_update_reference_Cap3",100,-0.5,0.5);
40  ratioUpRefCap3 = new TH1F("ratioUpRefCap3", "ration_update_reference_Cap3",100,0.5,1.5);
41  gainsUpCap3 = new TH1F("gainsUpCap3","gains_update_Cap3",100,0.0,0.6);
42  gainsRefCap3 = new TH1F("gainsRefCap3","gains_reference_Cap3",100,0.0,0.6);
43 
44  // gainsUpCap0vsEta = new TGraph("gainsUpCap0vsEta","gains_update_Cap0_vsEta",100,-41,0.6);
45  // gainsRefCap0vsEta = new TGraph("gainsRefCap0vsEta","gains_reference_Cap0_vsEta",100,0.0,0.6);
46 }
47 
48 
50 {
51  using namespace edm::eventsetup;
52  bool epsilonflag = false;
53  bool notequalsflag = false;
54  // get new gains
55  edm::ESHandle<HcalGains> newGains;
56  es.get<HcalGainsRcd>().get("update",newGains);
57  const HcalGains* myNewGains = newGains.product();
58 
59  // get reference gains
60  edm::ESHandle<HcalGains> refGains;
61  es.get<HcalGainsRcd>().get("reference",refGains);
62  const HcalGains* myRefGains = refGains.product();
63 
64  // get e-map from reference
66  es.get<HcalElectronicsMapRcd>().get("reference",refEMap);
67  const HcalElectronicsMap* myRefEMap = refEMap.product();
68 
69 
70  // dump gains:
71  if(dumpupdate.compare("null")!=0){
72  std::ofstream outStream(dumpupdate.c_str());
73  std::cout << "--- Dumping Gains - update ---" << std::endl;
74  HcalDbASCIIIO::dumpObject (outStream, (*myNewGains) );
75  }
76  if(dumprefs.compare("null")!=0){
77  std::ofstream outStream2(dumprefs.c_str());
78  std::cout << "--- Dumping Gains - reference ---" << std::endl;
79  HcalDbASCIIIO::dumpObject (outStream2, (*myRefGains) );
80  }
81  // get the list of all channels
82  std::vector<DetId> listNewChan = myNewGains->getAllChannels();
83  std::vector<DetId> listRefChan = myRefGains->getAllChannels();
84 
85  std::vector<DetId>::const_iterator cell;
86 
87  //plots: occupancy map, value map, difference, ratio, gains:
88  for (std::vector<DetId>::const_iterator it = listRefChan.begin(); it!=listRefChan.end(); it++)
89  {
90  HcalGenericDetId myId(*it);
91  // ocMapRef->Fill(myId->);
92 
93  float valCap0 = myRefGains->getValues(*it)->getValue(0);
94  float valCap1 = myRefGains->getValues(*it)->getValue(1);
95  float valCap2 = myRefGains->getValues(*it)->getValue(2);
96  float valCap3 = myRefGains->getValues(*it)->getValue(3);
97 
98  gainsRefCap0->Fill(valCap0);
99  gainsRefCap1->Fill(valCap1);
100  gainsRefCap2->Fill(valCap2);
101  gainsRefCap3->Fill(valCap3);
102 
103  cell = std::find(listNewChan.begin(), listNewChan.end(), (*it));
104  if (cell != listNewChan.end() ) //found
105  {
106  float valCap0up = myNewGains->getValues(*it)->getValue(0);
107  float valCap1up = myNewGains->getValues(*it)->getValue(1);
108  float valCap2up = myNewGains->getValues(*it)->getValue(2);
109  float valCap3up = myNewGains->getValues(*it)->getValue(3);
110 
111  diffUpRefCap0->Fill(valCap0up - valCap0);
112  diffUpRefCap1->Fill(valCap1up - valCap1);
113  diffUpRefCap2->Fill(valCap2up - valCap2);
114  diffUpRefCap3->Fill(valCap3up - valCap3);
115 
116  if(fabs(valCap0up - valCap0) > epsilon) epsilonflag = true;
117  if(fabs(valCap1up - valCap1) > epsilon) epsilonflag = true;
118  if(fabs(valCap2up - valCap2) > epsilon) epsilonflag = true;
119  if(fabs(valCap3up - valCap3) > epsilon) epsilonflag = true;
120 
121  if(valCap0up != valCap0) notequalsflag = true;
122  if(valCap1up != valCap1) notequalsflag = true;
123  if(valCap2up != valCap2) notequalsflag = true;
124  if(valCap3up != valCap3) notequalsflag = true;
125 
126  ratioUpRefCap0->Fill(valCap0up / valCap0);
127  ratioUpRefCap1->Fill(valCap1up / valCap1);
128  ratioUpRefCap2->Fill(valCap2up / valCap2);
129  ratioUpRefCap3->Fill(valCap3up / valCap3);
130  }
131  }
132  for (std::vector<DetId>::const_iterator it = listNewChan.begin(); it!=listNewChan.end(); it++)
133  {
134  float valCap0 = myNewGains->getValues(*it)->getValue(0);
135  float valCap1 = myNewGains->getValues(*it)->getValue(1);
136  float valCap2 = myNewGains->getValues(*it)->getValue(2);
137  float valCap3 = myNewGains->getValues(*it)->getValue(3);
138 
139  gainsUpCap0->Fill(valCap0);
140  gainsUpCap1->Fill(valCap1);
141  gainsUpCap2->Fill(valCap2);
142  gainsUpCap3->Fill(valCap3);
143  }
144 
145  if(epsilon != 1000000){
146  if(epsilonflag) throw cms::Exception("DataDoesNotMatch") << "Values differ by more than deltaG" << std::endl;
147  }else{
148  std::cout << "These gains do not differ by more than deltaG" << std::endl;
149  }
150 
151  if(validategainsflag){
152  if(notequalsflag) throw cms::Exception("DataDoesNotMatch") << "Values do not match" << std::endl;
153  }else{
154  std::cout << "These gains are identical" << std::endl;
155  }
156 
157  // go through list of valid channels from reference, look up if conditions exist for update
158  // push back into new vector the corresponding updated conditions,
159  // or if it doesn't exist, the reference
160 
161  if(outfile.compare("null")!=0){
162  HcalGains *resultGains = new HcalGains(refGains->topo());
163  for (std::vector<DetId>::const_iterator it = listRefChan.begin(); it != listRefChan.end(); it++)
164  {
165  DetId mydetid = *it;
166  HcalGenericDetId myId(*it);
167  cell = std::find(listNewChan.begin(), listNewChan.end(), mydetid);
168  if (cell == listNewChan.end()) // not present in new list, take old conditions
169  {
170  const HcalGain* item = myRefGains->getValues(mydetid);
171  std::cout << "o";
172  resultGains->addValues(*item);
173  }
174  else // present in new list, take new conditions
175  {
176  const HcalGain* item = myNewGains->getValues(mydetid);
177  std::cout << "n";
178  resultGains->addValues(*item);
179  }
180  }
181  std::cout << std::endl;
182 
183  std::vector<DetId> listResult = resultGains->getAllChannels();
184  // get the e-map list of channels
185  if(emapflag){
186  std::vector<HcalGenericDetId> listEMap = myRefEMap->allPrecisionId();
187  // look up if emap channels are all present in pedestals, if not then cerr
188  for (std::vector<HcalGenericDetId>::const_iterator it = listEMap.begin(); it != listEMap.end(); it++)
189  {
190  DetId mydetid = DetId(it->rawId());
191  if (std::find(listResult.begin(), listResult.end(), mydetid ) == listResult.end() )
192  {
193  std::cout << "Conditions not found for DetId = " << HcalGenericDetId(it->rawId()) << std::endl;
194  }
195  }
196  }
197 
198  // dump the resulting list of pedestals into a file
199  // std::ostringstream filename3;
200  // filename3 << "test_combined.txt";
201  std::ofstream outStream3(outfile.c_str());
202  std::cout << "--- Dumping Gains - the combined ones ---" << std::endl;
203  HcalDbASCIIIO::dumpObject (outStream3, (*resultGains) );
204 
205  }
206  // const float* values = myped->getValues (channelID);
207  // if (values) std::cout << "pedestals for channel " << channelID << ": "
208  // << values [0] << '/' << values [1] << '/' << values [2] << '/' << values [3] << std::endl;
209 
210 }
211 
212 
213 // ------------ method called once each job just after ending the event loop ------------
214 void
216 {
217  if(rootfile.compare("null")!=0){
218  f->Write();
219 }
220  f->Close();
221 
222 }
223 
TH1F * diffUpRefCap1
T getUntrackedParameter(std::string const &, T const &) const
virtual void beginJob()
TH1F * diffUpRefCap3
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
TH1F * ratioUpRefCap2
virtual void endJob()
bool ev
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
TH1F * diffUpRefCap2
TH1F * diffUpRefCap0
void analyze(const edm::Event &ev, const edm::EventSetup &es)
std::vector< DetId > getAllChannels() const
HcalGainsCheck(edm::ParameterSet const &ps)
const float * getValues() const
get value for all capId = 0..3
Definition: HcalGain.h:20
std::string dumpupdate
TH1F * ratioUpRefCap3
TH1F * ratioUpRefCap0
Definition: DetId.h:18
std::string outfile
const T & get() const
Definition: EventSetup.h:56
bool dumpObject(std::ostream &fOutput, const HcalPedestals &fObject)
TH1F * ratioUpRefCap1
tuple cout
Definition: gather_cfg.py:121
std::string dumprefs
bool addValues(const Item &myItem)
std::string rootfile