CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CSCAFEBThrAnalysis.cc
Go to the documentation of this file.
6 #include "TMath.h"
7 
8 class CSCFitAFEBThr;
9 
11 
12 hist_file=0; // set to null
13 
14 nmbev=0;
16 
17 npulses=100;
18 vecDac.clear();
19 BegDac=1;
20 EndDac=29;
21 EvDac=1;
22 StepDac=1;
23 indDac=0;
24 
25 vecDacOccup.assign(150,0);
26 m_wire_dac.clear();
27 m_res_for_db.clear();
28 mh_ChanEff.clear();
29 mh_FirstTime.clear();
30 mh_AfebDac.clear();
31 
32 }
33 
36  hist_file=new TFile(histoFileName.c_str(),"RECREATE");
37  hist_file->cd();
38 }
39 
40 void CSCAFEBThrAnalysis::bookForId(int flag, const int& idint,
41  const std::string& idstring ) {
42  hist_file->cd();
43 
44  std::ostringstream ss;
45 
46  if(flag==100) {
47  ss <<idint<<"_Anode_First_Time";
48  mh_FirstTime[idint]=new TH2F(ss.str().c_str(),"",675,0.0,675.0,50,0.0,10.0);
49  mh_FirstTime[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
50  mh_FirstTime[idint]->GetYaxis()->SetTitle("Anode First Time Bin");
51  mh_FirstTime[idint]->SetOption("BOX");
52  ss.str(""); // clear
53  }
54 
55  if(flag==101) {
56  ss <<idint<<"_Anode_Chan_Eff";
57  mh_ChanEff[idint]=new TH1F(ss.str().c_str(),"",675,0.0,675.0);
58  mh_ChanEff[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
59  mh_ChanEff[idint]->GetYaxis()->SetTitle("Entries");
60  ss.str(""); // clear
61  }
62 
63  if(flag==200) {
64  ss <<idint<<"_Anode_AfebDac";
65  mh_AfebDac[idint]=new TH2F(ss.str().c_str(),"",75,0.0,75.0, 50,0.0,50.0);
66  mh_AfebDac[idint]->GetXaxis()->SetTitle("Threshold DAC");
67  mh_AfebDac[idint]->GetYaxis()->SetTitle("AFEB Channel Occupancy");
68  mh_AfebDac[idint]->SetOption("COL");
69  ss.str(""); // clear
70  }
71 
72  if(flag==300) {
73  ss <<idint<<"_Anode_AfebThrPar";
74  mh_AfebThrPar[idint]=new TH2F(ss.str().c_str(),"",700,0.0,700.0, 50,0.0,50.0);
75  mh_AfebThrPar[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
76  mh_AfebThrPar[idint]->GetYaxis()->SetTitle("AFEB Channel Threshold (DAC)");
77  mh_AfebThrPar[idint]->SetOption("BOX");
78  ss.str(""); // clear
79  }
80 
81  if(flag==400) {
82  ss <<idint<<"_Anode_AfebNoisePar";
83  mh_AfebNoisePar[idint]=new TH2F(ss.str().c_str(),"",700,0.0,700.0, 50,0.0,5.0);
84  mh_AfebNoisePar[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
85  mh_AfebNoisePar[idint]->GetYaxis()->SetTitle("AFEB Channel Noise (DAC)");
86  mh_AfebNoisePar[idint]->SetOption("BOX");
87  ss.str(""); // clear
88  }
89 
90  if(flag==500) {
91  ss <<idint<<"_Anode_AfebNDF";
92  mh_AfebNDF[idint]=new TH2F(ss.str().c_str(),"",700,0.0,700.0, 25,-5.0,20.0);
93  mh_AfebNDF[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
94  mh_AfebNDF[idint]->GetYaxis()->SetTitle("AFEB Channel Fit NDF");
95  mh_AfebNDF[idint]->SetOption("BOX");
96  ss.str(""); // clear
97  }
98 
99  if(flag==600) {
100  ss <<idint<<"_Anode_AfebChi2perNDF";
101  mh_AfebChi2perNDF[idint]=new TH2F(ss.str().c_str(),"",700,0.0,700.0, 50,0.0,10.0);
102  mh_AfebChi2perNDF[idint]->GetXaxis()->SetTitle("(AFEB-1)*16+ch");
103  mh_AfebChi2perNDF[idint]->GetYaxis()->SetTitle("AFEB Channel Fit Chi2/NDF");
104  mh_AfebChi2perNDF[idint]->SetOption("BOX");
105  ss.str(""); // clear
106  }
107 }
108 
109 void CSCAFEBThrAnalysis::hf1ForId(std::map<int, TH1*>& mp, int flag,
110 const int& id, float& x, float w) {
111 
112  std::map<int,TH1*>::iterator h;
113  h=mp.find(id);
114  if (h==mp.end()) {
115  bookForId(flag,id,"");
116  h=mp.find(id);
117  }
118  h->second->Fill(x,w);
119 }
120 
121 void CSCAFEBThrAnalysis::hf2ForId(std::map<int, TH2*>& mp, int flag,
122 const int& id, float& x, float& y, float w) {
123 
124  std::map<int,TH2*>::iterator h;
125  h=mp.find(id);
126  if (h==mp.end()) {
127  bookForId(flag,id,"");
128  h=mp.find(id);
129  }
130  h->second->Fill(x,y,w);
131 }
132 
133 
134 /* Analyze the hits */
136 
137 std::ostringstream ss;
138 std::map<int,std::vector<int> >::iterator intIt;
139 std::map<int, std::vector<std::vector<int> > >::iterator wiredacIt;
140 
141 std::vector<int> vec;
142 int afeb,ch;
143 float x,y;
144 
145 m_wire_ev.clear();
146 
148 
149 nmbev++;
150 //indDac=(nmbev-1)/EvDac;
151 float dac=BegDac+StepDac*indDac;
152 if(vecDac.size() <= indDac) vecDac.push_back(dac);
153 
154 //std::cout<<" Event "<<nmbev;
155 //std::cout<<" "<<indDac<<" "<<vecDac[indDac]<<std::endl;
156 
157 //Anode wires
158 
160  if(wirecltn.begin() == wirecltn.end()) nmbev_no_wire++;
161 
162  if(wirecltn.begin() != wirecltn.end()) {
163 
165 
166  for (wiredetUnitIt=wirecltn.begin();
167  wiredetUnitIt!=wirecltn.end();
168  ++wiredetUnitIt){
169 
170  const CSCDetId& id = (*wiredetUnitIt).first;
171 
172  const int idchamber=id.endcap()*10000 +id.station()*1000+
173  id.ring()*100 +id.chamber();
174  const int idlayer =id.endcap()*100000+id.station()*10000+
175  id.ring()*1000+id.chamber()*10+id.layer();
176 
177  // std::cout<<idchamber<<" "<<idlayer<<std::endl;
178 
179  const int maxwire=csctoafeb.getMaxWire(id.station(),id.ring());
180  std::vector<int> wireplane(maxwire,0);
181 
182  const CSCWireDigiCollection::Range& range = (*wiredetUnitIt).second;
184  range.first; digiIt!=range.second; ++digiIt){
185 
186  int iwire=(*digiIt).getWireGroup();
187  if(iwire<=maxwire) {
188  if(wireplane[iwire-1]==0) {
189  wireplane[iwire-1]=(*digiIt).getBeamCrossingTag()+1;
190  ch=csctoafeb.getAfebCh(id.layer(),(*digiIt).getWireGroup());
191  afeb=csctoafeb.getAfebPos(id.layer(),(*digiIt).getWireGroup());
192 
194 
195  x=(afeb-1)*16+ch;
196  y=wireplane[iwire-1];
197  hf2ForId(mh_FirstTime, 100, idchamber,x, y, 1.0);
198 
199  } // end if wireplane[iwire-1]==0
200  } // end if iwire<=csctoafeb.getMaxWire(id.station(),id.ring()
201  } // end of for digis in layer
202 
204 
205  if(m_wire_ev.find(idlayer) == m_wire_ev.end())
206  m_wire_ev[idlayer]=wireplane;
207 
208  } // end of cycle on detUnit
209 
210 
212 
213  for(intIt=m_wire_ev.begin(); intIt!=m_wire_ev.end(); ++intIt) {
214  const int idwirev=(*intIt).first;
215  const std::vector<int> wiretemp=(*intIt).second; // What for?
216  int nsize=EndDac-BegDac+1;
217  std::vector<int> zer(nsize,0);
218 
219  wiredacIt=m_wire_dac.find(idwirev);
220  if(wiredacIt==m_wire_dac.end()) {
221  for(unsigned int j=0;j<wiretemp.size();j++)
222  m_wire_dac[idwirev].push_back(zer);
223  wiredacIt=m_wire_dac.find(idwirev);
224 // std::cout<<idwirev<<" "<<wiredacIt->second.size()<<" "<<
225 // wiredacIt->second[0].size()<<std::endl;
226  }
227  for(unsigned int i=0;i<(*intIt).second.size();i++)
228  if((*intIt).second[i]>0) wiredacIt->second[i][indDac]=
229  wiredacIt->second[i][indDac]+1;
230  } // end of adding hits to the map of wire/DAC
231  } // end of if wire collection not empty
232  indDac++;
233  if(dac==(float)EndDac) indDac=0;
234 } // end of void CSCAFEBThrAnalysis
235 
236 
238 
239  float x,y;
240 
241  //This is for DB transfer
242  CSCobject *cn = new CSCobject();
243  // Unused variable dbon caused compiler warning.
244  //condbon *dbon = new condbon();
245  new condbon();
246  // Note: cn and dbon pointers are never freed, I think. This leaks memory!!!
247 
248  std::map<int, std::vector<std::vector<int> > >::iterator mwiredacIt;
249  std::map<int, std::vector<std::vector<float> > >::iterator mresfordbIt;
250  std::vector<int>::iterator vecIt;
251 
252  std::cout<<"Events analyzed "<<nmbev<<std::endl;
253  std::cout<<"Events no anodes "<<nmbev_no_wire<<std::endl<<std::endl;
254 
255  std::cout<<"DAC occupancy"<<std::endl;
256  size_t ndacsize=EndDac-BegDac+1;
257  for(size_t i=0;i<ndacsize;i++) std::cout
258  <<" "<<vecDacOccup[i];
259  std::cout<<"\n\n"<<std::endl;
260 
261  std::vector<float> inputx;
262  std::vector<float> inputy;
263 
264  std::vector<float> mypar(2, 0.0);
265  std::vector<float> ermypar(2, 0.0);
266  float ercorr, chisq, edm;
267  int ndf,niter;
268 
269  int ch, afeb, idchmb;
270 
271  CSCFitAFEBThr * fitAnodeThr;
272 
273  std::vector<float> fitres(4,0);
274 
275 // std::cout<<"Threshold curves:\n"<<std::endl;
276 
277  for(mwiredacIt=m_wire_dac.begin();mwiredacIt!=m_wire_dac.end();++mwiredacIt){
278  int idwiredac=(*mwiredacIt).first;
279 
280  int layer=idwiredac-(idwiredac/10)*10;
281  idchmb=idwiredac/10;
282 
283  for(int unsigned iwire=0; iwire<mwiredacIt->second.size();iwire++) {
284 
285  afeb=csctoafeb.getAfebPos(layer,iwire+1);
286  ch=csctoafeb.getAfebCh(layer,iwire+1);
287  int afebid=(idwiredac/10)*100+csctoafeb.getAfebPos(layer,iwire+1);
288 
289  indDac=0;
290  for(vecIt=mwiredacIt->second[iwire].begin();
291  vecIt!=mwiredacIt->second[iwire].end(); ++vecIt) {
292 
293  x=vecDac[indDac]; y=*vecIt;
294  hf2ForId(mh_AfebDac, 200, afebid,x, y,1.0);
295 
297  if(indDac==0) {
298  x=(afeb-1)*16+ch;
299  y=*vecIt;
300  hf1ForId(mh_ChanEff, 101, idchmb, x, y);
301  }
302 
303  inputx.push_back(x);
304  inputy.push_back(y);
305 
306  indDac++;
307  }
308  // end of DAC cycle to form vectors of input data (inputx,inputy)for fit
309 
310 // std::cout<<afebid<<" "<<ch<<std::endl;
311 // for(unsigned int i=0;i<inputx.size();i++)
312 // std::cout<<" "<<inputy[i];
313 // std::cout<<std::endl;
314 
315 
316  for(unsigned int i=0;i<2;i++) {mypar[i]=0.0; ermypar[i]=0.0;}
317  ercorr=0.0;
318  chisq=0.0;
319  ndf=0;
320  niter=0;
321  edm=0.0;
322 
324  fitAnodeThr=new CSCFitAFEBThr();
325 
326 fitAnodeThr->ThresholdNoise(inputx,inputy,npulses,vecDacOccup,mypar,ermypar,ercorr,chisq,ndf,niter,edm);
327  delete fitAnodeThr;
328 
329 // std::cout<<"Fit "<<mypar[0]<<" "<<mypar[1]<<" "<<ndf<<" "<<chisq
330 // <<std::endl<<std::endl;
331 
333 
334  x=(afeb-1)*16+ch;
335 
337  y=mypar[0];
338  hf2ForId(mh_AfebThrPar, 300, idchmb,x, y, 1.0);
339 
341  y=mypar[1];
342  hf2ForId(mh_AfebNoisePar, 400, idchmb,x, y, 1.0);
343 
345  y=ndf;
346  hf2ForId(mh_AfebNDF, 500, idchmb,x, y, 1.0);
347 
349  y=0.0;
350  if(ndf>0) y=chisq/(float)ndf;
351  hf2ForId(mh_AfebChi2perNDF, 600, idchmb,x, y, 1.0);
352 
354  fitres[0]=mypar[0];
355  fitres[1]=mypar[1];
356  fitres[2]=ndf;
357  fitres[3]=0.0;
358  if(ndf>0) fitres[3]=chisq/(float)ndf;
359 
361 
362  mresfordbIt=m_res_for_db.find(idwiredac);
363  if(mresfordbIt==m_res_for_db.end())
364  m_res_for_db[idwiredac].push_back(fitres);
365  else m_res_for_db[idwiredac].push_back(fitres);
366 
367  inputx.clear();
368  inputy.clear();
369 
370 } // end for(int iwire=0)
371 } // end of iteration thru m_wire_dac map
372 
373  std::cout<<"Size of map for DB "<<m_res_for_db.size()<<std::endl;
374 
375  std::cout<<"The following CSCs will go to DB"<<std::endl<<std::endl;
376 
377  for(mresfordbIt=m_res_for_db.begin(); mresfordbIt!=m_res_for_db.end();
378  ++mresfordbIt) {
379  int idlayer=(*mresfordbIt).first;
380  idchmb=idlayer/10;
381  int layer=idlayer-idchmb*10;
382  std::cout<<"CSC "<<idchmb<<" Layer "<<layer<<" "
383  <<(*mresfordbIt).second.size()<<std::endl;
384  }
385 
386  for(mresfordbIt=m_res_for_db.begin(); mresfordbIt!=m_res_for_db.end();
387  ++mresfordbIt) {
388  int idlayer=(*mresfordbIt).first;
389 
390  //This is for DB transfer
391  int size = (*mresfordbIt).second.size();
392  cn->obj[idlayer].resize(size);
393 
394  for (unsigned int i=0;i<(*mresfordbIt).second.size();i++) {
395  std::cout<<idlayer<<" "<<i+1<<" ";
396 
397  for(int j=0;j<4;j++)
398  std::cout<< (*mresfordbIt).second[i][j]<<" ";
399  std::cout<<std::endl;
400 
401  //This is for DB transfer
402  cn->obj[idlayer][i].resize(4);
403  cn->obj[idlayer][i][0] = (*mresfordbIt).second[i][0];
404  cn->obj[idlayer][i][1] = (*mresfordbIt).second[i][1];
405  cn->obj[idlayer][i][2] = (*mresfordbIt).second[i][2];
406  cn->obj[idlayer][i][3] = (*mresfordbIt).second[i][3];
407  }
408  }
409 
410  //send data to DB
411  //dbon->cdbon_last_run("afeb_thresholds",&run);
412  //std::cout<<"Last AFEB thresholds run "<<run<<" for run file "<<myname<<" saved "<<myTime<<std::endl;
413  //if(debug) dbon->cdbon_write(cn,"afeb_thresholds",run+1,myTime);
414 
415  if(hist_file!=0) { // if there was a histogram file...
416  hist_file->Write(); // write out the histrograms
417  delete hist_file; // close and delete the file
418  hist_file=0; // set to zero to clean up
419  std::cout << "Hist. file was closed\n";
420  }
421 
422  std::cout<<" End of CSCAFEBThrAnalysis"<<std::endl;
423 }
int i
Definition: DBlmapReader.cc:9
int getAfebCh(int layer, int wiregroup) const
return AFEB channel number
Definition: CSCToAFEB.cc:14
std::map< int, TH2 * > mh_FirstTime
const double w
Definition: UKUtility.cc:23
tuple histoFileName
Definition: diJetCalib.py:108
std::map< int, std::vector< std::vector< float > > > m_res_for_db
std::map< int, std::vector< int > > m_wire_ev
Maps - per event, threshold curve, fit results.
virtual bool ThresholdNoise(const std::vector< float > &inputx, const std::vector< float > &inputy, const int &npulses, std::vector< int > &dacoccup, std::vector< float > &mypar, std::vector< float > &ermypar, float &ercorr, float &chisq, int &ndf, int &niter, float &edm) const
int endcap() const
Definition: CSCDetId.h:106
std::map< int, TH1 * > mh_ChanEff
Histogram maps.
T x() const
Cartesian x coordinate.
std::vector< int > vecDacOccup
std::map< int, TH2 * > mh_AfebChi2perNDF
int j
Definition: DBlmapReader.cc:9
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
int getMaxWire(int station, int ring) const
return max. number of wiregroups per layer
Definition: CSCToAFEB.cc:50
std::map< int, TH2 * > mh_AfebDac
std::vector< float > vecDac
const CSCToAFEB csctoafeb
Layer, wire to AFEB, channel conversion.
std::vector< CSCWireDigi >::const_iterator const_iterator
std::map< int, TH2 * > mh_AfebNDF
int nmbev
Statistics.
void hf1ForId(std::map< int, TH1 * > &mp, int flag, const int &id, float &x, float w)
TFile * hist_file
ROOT hist file.
int getAfebPos(int layer, int wiregroup) const
return AFEB position number
Definition: CSCToAFEB.cc:20
std::map< int, TH2 * > mh_AfebNoisePar
std::map< int, std::vector< std::vector< int > > > m_wire_dac
tuple cout
Definition: gather_cfg.py:121
std::pair< const_iterator, const_iterator > Range
void setup(const std::string &histoFileName)
void analyze(const CSCWireDigiCollection &wirecltn)
std::map< int, std::vector< std::vector< float > > > obj
Definition: CSCobject.h:12
std::map< int, TH2 * > mh_AfebThrPar
void hf2ForId(std::map< int, TH2 * > &mp, int flag, const int &id, float &x, float &y, float w)
tuple size
Write out results.
void bookForId(int flag, const int &idint, const std::string &ids)