Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
towerid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file towerid.cc
1 #include "towerid.h"
2 #include <fun4all/Fun4AllReturnCodes.h>
3 #include <phool/PHCompositeNode.h>
4 
5 //Fun4All
6 #include <fun4all/Fun4AllReturnCodes.h>
7 #include <fun4all/Fun4AllServer.h>
8 #include <fun4all/Fun4AllHistoManager.h>
9 #include <phool/PHCompositeNode.h>
10 #include <phool/getClass.h>
11 #include <phool/phool.h>
12 #include <ffaobjects/EventHeader.h>
13 
14 //ROOT
15 #include <TFile.h>
16 #include <TTree.h>
17 #include <TH1F.h>
18 #include <TH2F.h>
19 #include <TSystem.h>
20 #include <TPad.h>
21 #include <TLine.h>
22 #include <TStyle.h>
23 //Tower stuff
24 #include <calobase/TowerInfoContainer.h>
25 #include <calobase/TowerInfo.h>
26 #include <calobase/TowerInfoDefs.h>
27 
28 #include <cdbobjects/CDBTTree.h>
29 
30 //________________________________
31 towerid::towerid(const std::string &name, const std::string &cdbtreename, float adccut_sg, float adccut_k, float sigmas_lo, float sigmas_hi, float SG_f, float Kur_f, float region_f):
32 SubsysReco(name)
33 ,T(nullptr)
34  ,Outfile(name)
35  ,cdbtreename(cdbtreename)
36 ,adccut_sg(adccut_sg) // The minimum ADC required for a tower with Saint Gobain fibers to register a hit
37 ,adccut_k(adccut_k) // Minimum ADC for Kurary fibers to register a hit
38  ,sigmas_lo(sigmas_lo) // # of standard deviations from the mode for a cold tower
39  ,sigmas_hi(sigmas_hi) // # of standard deviations from the mode for a hot tower
40 ,SG_f(SG_f) // Fiducial cut (artificial maximum) for Saint Gobain towers
41  ,Kur_f(Kur_f) // Fiducial cut for Kurary
42  ,region_f(region_f) // Fiducial cut for Sectors and IBs
43 {
44  std::cout << "towerid::towerid(const std::string &name) Calling ctor" << std::endl;
45 }
46 //__________________________________
47 towerid::~towerid()
48 {
49  std::cout << "towerid::~towerid() Calling dtor" << std::endl;
50 }
51 //_____________________________
52 int towerid::Init(PHCompositeNode *topNode)
53 {
54  //initialize output file with hot channels
55 
56  //initialize tree with SG vs Kurary fiber information
57  fchannels = new TFile("channels.root");
58  channels = (TTree*)fchannels->Get("myTree");
59  channels->SetBranchAddress("fiber_type",&fiber_type);
60 
61  std::cout << "towerid::Init(PHCompositeNode *topNode) Initializing" << std::endl;
62  return Fun4AllReturnCodes::EVENT_OK;
63 
64 }
65 //_____________________________
66 int towerid::InitRun(PHCompositeNode *topNode)
67 {
68  std::cout << "towerid::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
69  return Fun4AllReturnCodes::EVENT_OK;
70 }
71 //_____________________________
72 int towerid::process_event(PHCompositeNode *topNode)
73 {
74 
75  //Get TowerInfoContainer
76  TowerInfoContainer *emcTowerContainer;
77  emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode,"TOWERS_CEMC");
78  if(!emcTowerContainer)
79  {
80  std::cout << PHWHERE << "towerid::process_event Could not find node TOWERS_CEMC" << std::endl;
81  return Fun4AllReturnCodes::ABORTEVENT;
82  }
83  //iterate through all towers, incrementing their Frequency arrays if they record a hit
84 
85  bool goodevent = false;
86  int tower_range = emcTowerContainer->size();
87  for(int j = 0; j < tower_range; j++){
88 
89  double energy = emcTowerContainer -> get_tower_at_channel(j) -> get_energy();
90  channels->GetEntry(j);
91  if((fiber_type ==0) && (energy > adccut_sg)){
92  towerF[j]++;
93  goodevent = true; //counter of events with nonzero EMCal energy
94  sectorF[j/nTowersSec]++;
95  ibF[j/nTowersIB]++;
96  }
97  else if ((fiber_type ==1) && (energy > adccut_k)){
98  towerF[j]++;
99  goodevent = true;
100  sectorF[j/nTowersSec]++;
101  ibF[j/nTowersIB]++;
102  }
103 
104  towerE[j]+=energy;
105  sectorE[j/nSectors]+=energy;
106  ibE[j/nTowersIB]+=energy;
107 
108 
109 
110  }
111  if(goodevent == true){goodevents++;}
112  return Fun4AllReturnCodes::EVENT_OK;
113 }
114 //__________________________
115 int towerid::ResetEvent(PHCompositeNode *topNode)
116 {
117  return Fun4AllReturnCodes::EVENT_OK;
118 }
119 //__________________________
120 void towerid::FillHistograms(const int runnumber, const int segment)
121 {
122 
123  const int nBins = 101;
124 
125  out = new TFile(Outfile.c_str(),"RECREATE");
126 
127  Fspec = new TH2F("Fspec","Fspec",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1+1/(2*nBins));
128  Fspec_SG = new TH2F("Fspec_SG","Fspec_SG",nTowers+1,-0.5,nTowers+0.5,nBins,0.-SG_f/(2*nBins),SG_f+SG_f/(2*nBins));
129  Fspec_K = new TH2F("Fspec_K","Fspec_K",nTowers+1,-0.5,nTowers+0.5,nBins,0.-Kur_f/(2*nBins),Kur_f+Kur_f/(2*nBins));
130  Fspec_sector = new TH2F("Fspec_sector","Fspec_sector",nSectors+1,-0.5,nSectors+0.5,nBins,0.-region_f/(2*nBins),region_f+region_f/(2*nBins));
131  Fspec_IB = new TH2F("Fspec_IB","Fspec_IB",nIB+1,-0.5,nIB+0.5,nBins,0.-region_f/(2*nBins),region_f+region_f/(2*nBins));
132 
133  Fspeci = new TH2F("Fspeci","Fspec_init",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
134  Fspeci_SG = new TH2F("Fspeci_SG","Fspeci_SG",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
135  Fspeci_K = new TH2F("Fspeci_K","Fspeci_K",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
136  Fspeci_sector = new TH2F("Fspeci_sector","Fspeci_sector",nSectors+1,-0.5,nSectors+0.5,nBins,0.-1/(2*nBins),1);
137  Fspeci_IB = new TH2F("Fspeci_IB","Fspeci_IB",nIB+1,-0.5,nIB+0.5,nBins,0.-1/(2*nBins),1);
138 
139  Espec = new TH2F("Espec","Espec",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
140  Espec_SG = new TH2F("Espec_SG","Espec_SG",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
141  Espec_K = new TH2F("Espec_K","Espec_K",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
142  Espec_sector = new TH2F("Espec_sector","Espec_sector",nSectors,-0.5,nSectors-0.5,nBins,0.-1/(2*nBins),1);
143  Espec_IB = new TH2F("Espec_IB","Espec_IB",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
144 
145  // hEventCounter = new TH1F("goodEventCounter","goodEventCounter",1,-0.5,0.5);
146  // hEventCounter -> SetBinContent(1,goodevents);
147 
148  //fill histograms
149 
150  for(int i = 0; i < nTowers; i++){
151 
152  Fspec->Fill(i,towerF[i]/goodevents);
153  Espec->Fill(i,towerE[i]/goodevents);
154  Fspeci->Fill(i,towerF[i]/goodevents);
155 
156  channels->GetEntry(i);
157  if(fiber_type == 0){
158  Fspec_SG->Fill(i,towerF[i]/goodevents);
159  Espec_SG->Fill(i,towerE[i]/goodevents);
160  Fspeci_SG->Fill(i,towerF[i]/goodevents);
161  }
162  else{
163  Fspec_K->Fill(i,towerF[i]/goodevents);
164  Espec_K->Fill(i,towerE[i]/goodevents);
165  Fspeci_K->Fill(i,towerF[i]/goodevents);
166  }
167  }
168  for(int j = 0; j<nIB; j++){
169  Fspec_IB->Fill(j,1.0*ibF[j]/nTowersIB/goodevents);
170  Espec_IB->Fill(j,1.0*ibE[j]/nTowersIB/goodevents);
171  Fspeci_IB->Fill(j,1.0*ibF[j]/nTowersIB/goodevents);
172  }
173  for(int k = 0; k<nSectors; k++){
174  Fspec_sector->Fill(k,1.0*sectorF[k]/nTowersSec/goodevents);
175  Espec_sector->Fill(k,1.0*sectorE[k]/nTowersSec/goodevents);
176  Fspeci_sector->Fill(k,1.0*sectorF[k]/nTowersSec/goodevents);
177  }
178 
179  //kill zeros:
180 
181  // // Fspec -> SetBinContent(1,0);
182  // // Fspec_SG -> SetBinContent(1,0);
183  // // Fspec_K -> SetBinContent(1,0);
184  // // Fspec_IB -> SetBinContent(1,0);
185  // // Fspec_sector -> SetBinContent(1,0);
186 
187  out -> cd();
188 
189  // hEventCounter -> Write();
190  // delete hEventCounter;
191 
192  Fspec->Write();
193  delete Fspec;
194 
195  Fspec_K->Write();
196  delete Fspec_K;
197 
198  Fspec_SG->Write();
199  delete Fspec_SG;
200 
201  Fspec_sector->Write();
202  delete Fspec_sector;
203 
204  Fspec_IB->Write();
205  delete Fspec_IB;
206 
207  Espec->Write();
208  delete Espec;
209 
210  Espec_SG->Write();
211  delete Espec_SG;
212 
213  Espec_K->Write();
214  delete Espec_K;
215 
216  Espec_sector->Write();
217  delete Espec_sector;
218 
219  Espec_IB->Write();
220  delete Espec_IB;
221 
222  Fspeci->Write();
223  delete Fspeci;
224 
225  Fspeci_K->Write();
226  delete Fspeci_K;
227 
228  Fspeci_SG->Write();
229  delete Fspeci_SG;
230 
231  Fspeci_sector->Write();
232  delete Fspeci_sector;
233 
234  Fspeci_IB->Write();
235  delete Fspeci_IB;
236 
237  out -> Close();
238  delete out;
239 
240 }
241 //__________________________
242 void towerid::CalculateCutOffs(const int runnumber)
243 {
244  TFile *histsIn = new TFile(Form("output/%d/DST_CALOR-%08d_badTowerMapTree.root",runnumber,runnumber));
245 
246  float cutoffFreq_SG;
247  float cutoffFreq_K;
248  float cutoffFreq;
249 
250  // float cutoffFreq_sector;
251  float cutoffFreq_IB;
252 
253  float cutoffFreq_SG_lo;
254  float cutoffFreq_K_lo;
255 
256  // float cutoffFreq_sector_lo;
257  float cutoffFreq_IB_lo;
258 
259  // hEventCounter = (TH1F*)histsIn -> Get("goodEventCounter");
260  // goodevents = hEventCounter -> GetBinContent(1);
261 
262 
263  //calculate initial cut off values
264  TH1F *dummyProj = NULL;
265  gStyle -> SetOptStat(0);
266  Fspec_SG = (TH2F*)histsIn -> Get("Fspec_SG");
267  dummyProj = (TH1F*)Fspec_SG -> ProjectionY("dummy");
268  dummyProj -> SetBinContent(1,0);
269  cutoffFreq_SG = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
270  cutoffFreq_SG_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
271 
272  Fspec_K = (TH2F*)histsIn -> Get("Fspec_K");
273  dummyProj = (TH1F*)Fspec_K -> ProjectionY("dummy");
274  dummyProj -> SetBinContent(1,0);
275 
276  cutoffFreq_K = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
277  cutoffFreq_K_lo = -1*Fspec_K->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
278 
279 
280  Fspec = (TH2F*)histsIn -> Get("Fspec");
281  dummyProj = (TH1F*)Fspec -> ProjectionY("dummy");
282  dummyProj -> SetBinContent(1,0);
283  cutoffFreq = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
284 
285 
286  Fspec_IB = (TH2F*)histsIn -> Get("Fspec_IB");
287  dummyProj = (TH1F*)Fspec_IB -> ProjectionY("dummy");
288  dummyProj -> SetBinContent(1,0);
289 
290  cutoffFreq_IB = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
291  cutoffFreq_IB_lo = -1*Fspec_IB->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
292 
293  // Fspec_sector = (TH2F*)histsIn -> Get("Fspec_sector");
294  // dummyProj = (TH1F*)Fspec_sector -> ProjectionY("dummy");
295  // dummyProj -> SetBinContent(1,0);
296 
297  // cutoffFreq_sector = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
298  // cutoffFreq_sector_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
299 
300  //get the rest of the histograms we nede
301  Fspeci_K = (TH2F*)histsIn -> Get("Fspeci_K");
302 
303  Fspeci_SG = (TH2F*)histsIn -> Get("Fspeci_SG");
304  Fspeci_SG -> SetTitle(";Tower ID; Hit Frequency");
305 
306 
307  Fspeci = (TH2F*)histsIn -> Get("Fspeci");
308  Fspeci_IB = (TH2F*)histsIn -> Get("Fspeci_IB");
309  Fspeci_IB -> SetTitle(";IB ID; Hit Frequency");
310 
311  Fspeci_sector = (TH2F*)histsIn -> Get("Fspeci_sector");
312 
313 
314 
315  //initial hot tower calculation
316  for(int i = 0; i < nTowers; i++){
317 
318  channels->GetEntry(i);
319 
320  TH1F *fSpecProj_SG = (TH1F*)Fspeci_SG->ProjectionY("dummySG",i+1,i+1);
321  TH1F *fSpecProj_K = (TH1F*)Fspeci_K->ProjectionY("dummyK",i+1,i+1);
322  TH1F *fSpecProj = (TH1F*)Fspeci->ProjectionY("dummySpec",i+1,i+1);
323 
324  if(fiber_type == 0 && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) > cutoffFreq_SG){
325  hottowers[i]++;
326  }
327  else if(fiber_type == 1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > cutoffFreq_K){
328  hottowers[i]++;
329  }
330  else if(fSpecProj -> GetBinCenter(fSpecProj->GetMaximumBin()) == 0){
331  deadtowers[i]++;
332  }
333  else if (fiber_type == 0 && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) < cutoffFreq_SG_lo && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) >0 ){
334  coldtowers[i]++;
335  }
336  else if (fiber_type ==1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) < cutoffFreq_K_lo && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > 0){
337  coldtowers[i]++;
338  }
339 
340  }
341  //go through and look for hot/cold "regions", e.g. sectors and interface boards
342  hot_regions = 0;
343  cold_regions = 0;
344 
345  //Interface board QA, create projection
346 
347  //hotIB
348  for(int j = 0; j < nIB; j++){
349  TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
350 
351  if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())>cutoffFreq_IB){
352  hot_regions = 1;
353  std::cout << "IB " << j << " is hot with ADC rate" << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())<< " > "<< cutoffFreq_IB << std::endl;
354 
355  hotIB[j]++;
356  for(int j1 = 0; j1 < nTowersIB; j1++ ){
357  hottowers[j*nTowersIB+j1]++;
358  for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1, biny+1,-10);
359  //towerF[j*nTowersIB+j1]=0;//mask tower's contribution so it's no longer used to calculate cutoffs
360  }
361  for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB->SetBinContent(j+1,biny+1,-10);
362  //ibF[j]=-1;
363  }
364  }
365 
366  //cold IB
367  for(int j = 0; j < nIB; j++){
368  TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
369 
370  if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())<cutoffFreq_IB_lo){
371  hot_regions = 1;
372  std::cout << "IB " << j << " is cold with ADC rate" << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin()) << " < " << cutoffFreq_IB_lo << std::endl;
373  hotIB[j]++;
374 
375  for(int j1 = 0; j1< nTowersIB; j1++ ){
376  hottowers[j*nTowersIB+j1]++;
377  for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1,biny+1,-10);
378  //towerF[j*nTowersIB+j1]=0;
379  }
380  //ibF[j] = -1;
381  for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB -> SetBinContent(j+1,biny+1,-10);
382  }
383  }
384 
385  /*Removing sector level qa
386  Essentially the issue is that removing IB's biases the determination of hot and cold sectors.
387  Like removing one IB's worth of towers is just going to make a sector cold, or having 1 hot IB can make a sector hot.
388  One could, in principle, re-scale the response by the fraction of missing IB's, but it isn't clear that that's necessary.
389  Towers are necessarily effected by the behavior of their IB, but it isn't clear that a single sector can be driven
390  bad by the behavior of a single IB. - AH*/
391 
392 
393  //QA the sectors
394  //hot sector
395  // for(int k = 0; k < nSectors; k++){
396  // TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);
397 
398  // if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) > cutoffFreq_sector){
399  // hot_regions = 1;
400  // std::cout << "sector " << k << "is hot with ADC rate " << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " > " << cutoffFreq_sector << std::endl;
401  // hotsectors[k]++;
402 
403  // for(int k1 = 0; k1 < nTowersSec; k1++){
404  // hottowers[k*nTowersSec+k1]++;
405  // towerF[k*nTowersSec+k1] = 0;
406  // //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
407  // }
408  // sectorF[k] = -1;
409  // //fSpecSecProj->SetBinContent(k,0);
410  // //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
411  // }
412  // }
413 
414 
415 
416  // //cold sector
417  // for(int k = 0; k < nSectors; k++){
418  // TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);
419 
420  // if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) < cutoffFreq_sector_lo){
421  // cold_regions = 1;
422  // std::cout << "sector " << k << "is cold with ADC rate" << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " < " << cutoffFreq_sector_lo << std::endl;
423  // coldsectors[k]++;
424 
425  // for(int k1 = 0; k1 < nIB; k1++){
426  // coldtowers[k*nIB+k1]++;
427  // //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
428  // towerF[k*nTowersSec+k1] = 0;
429  // }
430  // //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
431  // sectorF[k]=-1;
432  // }
433  // }
434 
435  std::cout << "removing hot/cold regions" << std::endl;
436 
437  //go through again and calculate cut offs with bad regions masked until the calculation is made without hot regions
438  while(hot_regions == 1 || cold_regions ==1){
439  std::cout << "hot/cold IB or sector detected. Running another pass for hot towers" << std::endl;
440  // Fspec->Reset();
441  // Fspec_SG->Reset();
442  // Fspec_K->Reset();
443  // Fspec_IB->Reset();
444  //Fspec_sector->Reset();
445 
446  //not used for now
447  // Espec->Reset();
448  // Espec_SG->Reset();
449  // Espec_K->Reset();
450  // Espec_IB->Reset();
451  // Espec_sector->Reset();
452  //for(int i = 0; i < nTowers; i++){
453 
454  //Fspec->Fill(i,towerF[i]/goodevents);
455  // Espec->Fill(i,towerE[i]/goodevents);
456  //channels->GetEntry(i);
457  //if(fiber_type == 0){
458  //Fspec_SG->Fill(i,towerF[i]/goodevents);
459  //Espec_SG->Fill(i,towerE[i]/goodevents);
460  //}
461  //else{
462  //Fspec_K->Fill(i,towerF[i]/goodevents);
463  //Espec_K->Fill(i,towerE[i]/goodevents);
464  //}
465  //}
466 
467  //for(int j = 0; j<nIB; j++){
468  //Fspec_IB->Fill(j,1.0*ibF[j]/nTowersIB/goodevents);
469  //Espec_IB->Fill(1,1.0*ibE[j]/64/goodevents);
470  //}
471  //for(int k = 0; k<nSectors; k++){
472  //Fspec_sector->Fill(k,1.0*sectorF[k]/nTowersSec/goodevents);
473  //Espec_sector->Fill(1,1.0*sectorE[k]/384/goodevents);
474  //}
475 
476  //kill zeros
477 
478  // Fspec->SetBinContent(1,0);
479  // Fspec_SG->SetBinContent(1,0);
480  // Fspec_K->SetBinContent(1,0);
481  // Fspec_IB->SetBinContent(1,0);
482  // Fspec_sector->SetBinContent(1,0);
483 
484  dummyProj = (TH1F*)Fspec_SG -> ProjectionY("dummy");
485  dummyProj -> SetBinContent(1,0);
486 
487 
488  cutoffFreq_SG = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
489  cutoffFreq_SG_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
490 
491  dummyProj = (TH1F*)Fspec_K -> ProjectionY("dummy");
492  dummyProj -> SetBinContent(1,0);
493 
494  cutoffFreq_K = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
495  cutoffFreq_K_lo = -1*Fspec_K->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
496 
497  dummyProj = (TH1F*)Fspec -> ProjectionY("dummy");
498  dummyProj -> SetBinContent(1,0);
499  cutoffFreq = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
500 
501  dummyProj = (TH1F*)Fspec_IB -> ProjectionY("dummy");
502  dummyProj -> SetBinContent(1,0);
503  cutoffFreq_IB = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
504  cutoffFreq_IB_lo = -1*Fspec_IB->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
505 
506  // dummyProj = (TH1F*)Fspec_sector -> ProjectionY("dummy");
507  // dummyProj -> SetBinContent(1,0);
508  // cutoffFreq_sector = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
509  // cutoffFreq_sector_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
510 
511 
512  for(int i = 0; i < nTowers; i++){
513 
514  channels->GetEntry(i);
515  TH1F *fSpecProj_SG = (TH1F*)Fspeci_SG->ProjectionY("dummySG",i+1,i+1);
516  TH1F *fSpecProj_K = (TH1F*)Fspeci_K->ProjectionY("dummyK",i+1,i+1);
517  TH1F *fSpecProj = (TH1F*)Fspeci->ProjectionY("dummySpec",i+1,i+1);
518 
519  if(fiber_type == 0 && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) > cutoffFreq_SG){
520  hottowers[i]++;
521  }
522  else if(fiber_type == 1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > cutoffFreq_K){
523  hottowers[i]++;
524  }
525  else if(fSpecProj -> GetBinCenter(fSpecProj->GetMaximumBin()) ==0){
526  deadtowers[i]++;
527  }
528  else if (fiber_type == 0 && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) < cutoffFreq_SG_lo && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) >0 ){
529  coldtowers[i]++;
530  }
531  else if (fiber_type ==1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) < cutoffFreq_K_lo && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > 0){
532  coldtowers[i]++;
533  }
534  }
535 
536  hot_regions = 0;
537  cold_regions = 0;
538 
539  //Re-QA IB's
540  //hotIB
541  for(int j = 0; j< nIB; j++){
542  if(ibF[j]==-1)continue; //IB has been marked hot already, don't bother with it
543  TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
544 
545  if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())>cutoffFreq_IB){
546  hot_regions = 1;
547  std::cout << "IB " << j << "is hot with ADC rate " << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin()) << " > " << cutoffFreq_IB << std::endl;
548  hotIB[j]++;
549 
550  for(int j1 = 0; j1 < nTowersIB; j1++ ){
551  hottowers[j*nTowersIB+j1]++;
552  //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1,biny+1, 0);
553  towerF[j*nTowersIB+j1] = 0;//mask tower's contribution so it's no longer used to calculate cutoffs
554  }
555  //for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB->SetBinContent(j+1,biny+1,0);
556  ibF[j] = -1;
557  }
558  }
559 
560  //cold IB
561  for(int j = 0; j<nIB; j++){
562  if(ibF[j]==-1)continue; //IB has been marked hot already, don't bother with it
563 
564  TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
565 
566  if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())<cutoffFreq_IB_lo){
567  hot_regions = 1;
568  std::cout << "IB " << j << " is cold with ADC rate " << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin()) << " < " << cutoffFreq_IB_lo << std::endl;
569  hotIB[j]++;
570 
571  for(int j1 = 0; j1< nTowersIB; j1++ ){
572  hottowers[j*nTowersIB+j1]++;
573  //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1,biny+1,0);
574  towerF[j*nTowersIB+j1] = 0;
575  }
576  ibF[j] = -1;
577  //for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB->SetBinContent(j+1, biny+1,0);
578  }
579  }
580 
581  //Re-QA sectors
582  //hot sectors
583  // for(int k = 0; k<nSectors; k++){
584  // if(sectorF[k]==-1) continue;//sector has been marked hot already
585  // TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);
586 
587  // if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) > cutoffFreq_sector){
588  // hot_regions = 1;
589  // std::cout << "sector " << k << " is hot with ADC rate " << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " > " <<cutoffFreq_sector << std::endl;
590  // hotsectors[k]++;
591 
592  // for(int k1 = 0; k1 < nTowersSec; k1++){
593  // hottowers[k*nTowersSec+k1]++;
594 
595  // //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
596  // towerF[k*nTowersSec+k1] = 0;
597  // }
598  // sectorF[k] = -1;
599  // //fSpecSecProj->SetBinContent(k,0);
600  // //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
601  // }
602  // }
603 
604  //cold sector
605  // for(int k = 0; k<nSectors; k++){
606  // if(sectorF[k]==-1) continue;
607  // TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);
608 
609  // if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) <cutoffFreq_sector_lo){
610  // cold_regions = 1;
611  // std::cout << "sector " << k << " is cold with ADC rate " << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " < " <<cutoffFreq_sector_lo << std::endl;
612  // coldsectors[k]++;
613 
614  // for(int k1 = 0; k1 < nIB; k1++){
615  // coldtowers[k*nIB+k1]++;
616  // //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
617  // towerF[k*nTowersSec+k1] = 0;
618  // }
619  // //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
620  // sectorF[k] = -1;
621  // }
622  // }
623 
624  }
625 
626  std::cout << "towerid::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
627  std::cout << "Saint Gobain Cutoff Frequency: " << cutoffFreq_SG << std::endl;
628  std::cout << "Kurary Cutoff Frequency: " << cutoffFreq_K << std::endl;
629  std::cout << "Overall Tower Cutoff Frequency: " << cutoffFreq << std::endl;
630  std::cout << "IB Cutoff hit Frequency: " << cutoffFreq_IB << std::endl;
631  //std::cout << "Sector Cutoff hit Frequency: " << cutoffFreq_sector << std::endl;
632 }
633 
634 //____________________________________________________________________________..
635 void towerid::WriteCDBTree(const int runnumber)
636 {
637 
638 
639  //TFile *cdbOut = new TFile(Form("cdbMaps/%d/DST_CALOR-%08d_badTowerMapCDBTTree.root",runnumber,runnumber),"RECREATE");
640  TFile *cdbOut = new TFile(Form("cdbMaps/%d/CEMC-%08d_badTowerMapCDBTTree.root",runnumber,runnumber),"RECREATE");
641 
642  T = new TTree("T","T");
643  T -> Branch("hot_channels",&m_hot_channels);
644 
645  CDBTTree*cdbttree = new CDBTTree(cdbtreename);
646  std::string fieldname = "status";
647  for(int i = 0; i<nTowers; i++){
648 
649  unsigned int key = TowerInfoDefs::encode_emcal(i);
650 
651  if(hottowers[i] >= 0.5){
652  m_hot_channels = 2;
653  T->Fill();
654  cdbttree->SetIntValue(key,fieldname,2);
655  }
656  else if(deadtowers[i] >= 0.5){
657  m_hot_channels = 1;
658  T->Fill();
659  cdbttree->SetIntValue(key,fieldname,1);
660  }
661  else if(coldtowers[i] > 0.5){
662  m_hot_channels = 3;
663  T->Fill();
664  cdbttree->SetIntValue(key,fieldname,3);
665  }
666  else{
667  m_hot_channels = 0;
668  T->Fill();
669  cdbttree->SetIntValue(key,fieldname,0);
670  }
671  }
672 
673  fchannels->Close();
674  delete fchannels;
675 
676  cdbOut->cd();
677  T -> Write();
678  cdbttree->Commit();
679  //cdbttree->Print();
680  cdbttree->WriteCDBTTree();
681  delete cdbttree;
682  cdbOut->Close();
683 
684 }
685 
686 //____________________________________________________________________________..
687 int towerid::Reset(PHCompositeNode *topNode)
688 {
689  std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
690  return Fun4AllReturnCodes::EVENT_OK;
691 }
692 //____________________________________________________
693 void towerid::Print(const std::string &what) const
694 {
695  std::cout << "towerid::Print(const std::string &what) const Printing info for " << what << std::endl;
696 }
697 //____________________________________________________
698 int towerid::End(PHCompositeNode *topNode)
699 {
700  std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
701  return Fun4AllReturnCodes::EVENT_OK;
702 }
703 //____________________________________________________
704 int towerid::EndRun(int runnumber)
705 {
706  std::cout << "towerid::EndRun: this is the end of run: "<< runnumber << std::endl;
707  return Fun4AllReturnCodes::EVENT_OK;
708 }
709 //____________________________________________________
710 
711 
712 
713 
714 
715 
716