Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:13:10 

0001 #include "towerid.h"
0002 #include <fun4all/Fun4AllReturnCodes.h>
0003 #include <phool/PHCompositeNode.h>
0004 
0005 //Fun4All
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/Fun4AllServer.h>
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>
0012 #include <ffaobjects/EventHeader.h>
0013 
0014 //ROOT
0015 #include <TFile.h>
0016 #include <TTree.h>
0017 #include <TH1F.h>
0018 #include <TH2F.h>
0019 #include <TSystem.h>
0020 #include <TPad.h>
0021 #include <TLine.h>
0022 #include <TStyle.h>
0023 //Tower stuff
0024 #include <calobase/TowerInfoContainer.h>
0025 #include <calobase/TowerInfo.h>
0026 #include <calobase/TowerInfoDefs.h>
0027 
0028 #include <cdbobjects/CDBTTree.h>
0029 
0030 //________________________________
0031 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): 
0032 SubsysReco(name)
0033 ,T(nullptr)
0034   ,Outfile(name)    
0035   ,cdbtreename(cdbtreename)
0036 ,adccut_sg(adccut_sg)       // The minimum ADC required for a tower with Saint Gobain fibers to register a hit
0037 ,adccut_k(adccut_k)     // Minimum ADC for Kurary fibers to register a hit
0038   ,sigmas_lo(sigmas_lo)     // # of standard deviations from the mode for a cold tower
0039   ,sigmas_hi(sigmas_hi)     // # of standard deviations from the mode for a hot tower
0040 ,SG_f(SG_f)         // Fiducial cut (artificial maximum) for Saint Gobain towers
0041   ,Kur_f(Kur_f)         // Fiducial cut for Kurary
0042   ,region_f(region_f)       // Fiducial cut for Sectors and IBs 
0043 {
0044   std::cout << "towerid::towerid(const std::string &name) Calling ctor" << std::endl;
0045 }
0046 //__________________________________
0047 towerid::~towerid()
0048 {
0049   std::cout << "towerid::~towerid() Calling dtor" << std::endl;
0050 }
0051 //_____________________________
0052 int towerid::Init(PHCompositeNode *topNode)
0053 {
0054   //initialize output file with hot channels
0055   
0056   //initialize tree with SG vs Kurary fiber information
0057   fchannels = new TFile("channels.root");
0058   channels = (TTree*)fchannels->Get("myTree");
0059   channels->SetBranchAddress("fiber_type",&fiber_type);
0060   
0061   std::cout << "towerid::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0062   return Fun4AllReturnCodes::EVENT_OK;
0063 
0064 }
0065 //_____________________________
0066 int towerid::InitRun(PHCompositeNode *topNode)
0067 {
0068   std::cout << "towerid::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0069   return Fun4AllReturnCodes::EVENT_OK;
0070 }
0071 //_____________________________
0072 int towerid::process_event(PHCompositeNode *topNode)
0073 {
0074 
0075   //Get TowerInfoContainer
0076   TowerInfoContainer *emcTowerContainer;
0077   emcTowerContainer =  findNode::getClass<TowerInfoContainer>(topNode,"TOWERS_CEMC");
0078   if(!emcTowerContainer)
0079     {
0080       std::cout << PHWHERE << "towerid::process_event Could not find node TOWERS_CEMC"  << std::endl;
0081       return Fun4AllReturnCodes::ABORTEVENT;
0082     }
0083   //iterate through all towers, incrementing their Frequency arrays if they record a hit
0084 
0085   bool goodevent = false;
0086   int tower_range = emcTowerContainer->size();
0087   for(int j = 0; j < tower_range; j++){
0088 
0089     double energy = emcTowerContainer -> get_tower_at_channel(j) -> get_energy();
0090     channels->GetEntry(j);
0091     if((fiber_type ==0) && (energy > adccut_sg)){
0092       towerF[j]++;
0093       goodevent = true; //counter of events with nonzero EMCal energy
0094       sectorF[j/nTowersSec]++;
0095       ibF[j/nTowersIB]++;
0096     }
0097     else if ((fiber_type ==1) && (energy > adccut_k)){
0098       towerF[j]++;
0099       goodevent = true;
0100       sectorF[j/nTowersSec]++;
0101       ibF[j/nTowersIB]++;
0102     }
0103                 
0104     towerE[j]+=energy;
0105     sectorE[j/nSectors]+=energy;
0106     ibE[j/nTowersIB]+=energy;
0107         
0108             
0109    
0110   }
0111   if(goodevent == true){goodevents++;}
0112   return Fun4AllReturnCodes::EVENT_OK;
0113 }
0114 //__________________________
0115 int towerid::ResetEvent(PHCompositeNode *topNode)
0116 {
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119 //__________________________
0120 void towerid::FillHistograms(const int runnumber, const int segment)
0121 {   
0122   
0123   const int nBins = 101;
0124  
0125   out = new TFile(Outfile.c_str(),"RECREATE");
0126 
0127   Fspec = new TH2F("Fspec","Fspec",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1+1/(2*nBins));
0128   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));
0129   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));  
0130   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));
0131   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));
0132   
0133   Fspeci = new TH2F("Fspeci","Fspec_init",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
0134   Fspeci_SG = new TH2F("Fspeci_SG","Fspeci_SG",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
0135   Fspeci_K = new TH2F("Fspeci_K","Fspeci_K",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
0136   Fspeci_sector = new TH2F("Fspeci_sector","Fspeci_sector",nSectors+1,-0.5,nSectors+0.5,nBins,0.-1/(2*nBins),1);
0137   Fspeci_IB = new TH2F("Fspeci_IB","Fspeci_IB",nIB+1,-0.5,nIB+0.5,nBins,0.-1/(2*nBins),1);
0138 
0139   Espec = new TH2F("Espec","Espec",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
0140   Espec_SG = new TH2F("Espec_SG","Espec_SG",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
0141   Espec_K = new TH2F("Espec_K","Espec_K",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
0142   Espec_sector = new TH2F("Espec_sector","Espec_sector",nSectors,-0.5,nSectors-0.5,nBins,0.-1/(2*nBins),1);
0143   Espec_IB = new TH2F("Espec_IB","Espec_IB",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
0144   
0145   // hEventCounter = new TH1F("goodEventCounter","goodEventCounter",1,-0.5,0.5);
0146   // hEventCounter -> SetBinContent(1,goodevents);
0147   
0148   //fill histograms
0149  
0150   for(int i = 0; i < nTowers; i++){
0151 
0152     Fspec->Fill(i,towerF[i]/goodevents);
0153     Espec->Fill(i,towerE[i]/goodevents);
0154     Fspeci->Fill(i,towerF[i]/goodevents);
0155 
0156     channels->GetEntry(i);
0157     if(fiber_type == 0){
0158       Fspec_SG->Fill(i,towerF[i]/goodevents);
0159       Espec_SG->Fill(i,towerE[i]/goodevents);
0160       Fspeci_SG->Fill(i,towerF[i]/goodevents);
0161     }
0162     else{
0163       Fspec_K->Fill(i,towerF[i]/goodevents);
0164       Espec_K->Fill(i,towerE[i]/goodevents);
0165       Fspeci_K->Fill(i,towerF[i]/goodevents);
0166     }
0167   }
0168   for(int j = 0; j<nIB; j++){
0169     Fspec_IB->Fill(j,1.0*ibF[j]/nTowersIB/goodevents);
0170     Espec_IB->Fill(j,1.0*ibE[j]/nTowersIB/goodevents);
0171     Fspeci_IB->Fill(j,1.0*ibF[j]/nTowersIB/goodevents);
0172   }
0173   for(int k = 0; k<nSectors; k++){
0174     Fspec_sector->Fill(k,1.0*sectorF[k]/nTowersSec/goodevents);
0175     Espec_sector->Fill(k,1.0*sectorE[k]/nTowersSec/goodevents);
0176     Fspeci_sector->Fill(k,1.0*sectorF[k]/nTowersSec/goodevents);
0177   }
0178 
0179   //kill zeros: 
0180 
0181   // // Fspec -> SetBinContent(1,0);
0182   // // Fspec_SG -> SetBinContent(1,0);
0183   // // Fspec_K -> SetBinContent(1,0);
0184   // // Fspec_IB -> SetBinContent(1,0);
0185   // // Fspec_sector -> SetBinContent(1,0);
0186   
0187   out -> cd();
0188   
0189   // hEventCounter -> Write();
0190   // delete hEventCounter;
0191 
0192   Fspec->Write();
0193   delete Fspec;
0194   
0195   Fspec_K->Write();
0196   delete Fspec_K;
0197   
0198   Fspec_SG->Write();
0199   delete Fspec_SG;
0200   
0201   Fspec_sector->Write();
0202   delete Fspec_sector;
0203 
0204   Fspec_IB->Write();
0205   delete Fspec_IB;
0206   
0207   Espec->Write();
0208   delete Espec;
0209 
0210   Espec_SG->Write();
0211   delete Espec_SG;
0212 
0213   Espec_K->Write();
0214   delete Espec_K;
0215   
0216   Espec_sector->Write();
0217   delete Espec_sector;
0218 
0219   Espec_IB->Write();
0220   delete Espec_IB;
0221 
0222   Fspeci->Write();
0223   delete Fspeci;
0224 
0225   Fspeci_K->Write();
0226   delete Fspeci_K;
0227 
0228   Fspeci_SG->Write();
0229   delete Fspeci_SG;
0230 
0231   Fspeci_sector->Write();
0232   delete Fspeci_sector;
0233   
0234   Fspeci_IB->Write();
0235   delete Fspeci_IB; 
0236   
0237   out -> Close();
0238   delete out;
0239   
0240 }
0241 //__________________________
0242 void towerid::CalculateCutOffs(const int runnumber)
0243 {   
0244   TFile *histsIn = new TFile(Form("output/%d/DST_CALOR-%08d_badTowerMapTree.root",runnumber,runnumber));
0245     
0246   float cutoffFreq_SG;
0247   float cutoffFreq_K;
0248   float cutoffFreq;
0249     
0250   // float cutoffFreq_sector;
0251   float cutoffFreq_IB;
0252 
0253   float cutoffFreq_SG_lo;
0254   float cutoffFreq_K_lo;
0255 
0256   // float cutoffFreq_sector_lo;
0257   float cutoffFreq_IB_lo;
0258 
0259   // hEventCounter = (TH1F*)histsIn -> Get("goodEventCounter");
0260   // goodevents = hEventCounter -> GetBinContent(1);
0261    
0262  
0263   //calculate initial cut off values
0264   TH1F *dummyProj = NULL;
0265   gStyle -> SetOptStat(0);
0266   Fspec_SG = (TH2F*)histsIn -> Get("Fspec_SG");
0267   dummyProj = (TH1F*)Fspec_SG -> ProjectionY("dummy");
0268   dummyProj -> SetBinContent(1,0);
0269   cutoffFreq_SG = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0270   cutoffFreq_SG_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0271 
0272   Fspec_K = (TH2F*)histsIn -> Get("Fspec_K");
0273   dummyProj = (TH1F*)Fspec_K -> ProjectionY("dummy");
0274   dummyProj -> SetBinContent(1,0);
0275 
0276   cutoffFreq_K = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0277   cutoffFreq_K_lo = -1*Fspec_K->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0278   
0279 
0280   Fspec = (TH2F*)histsIn -> Get("Fspec");
0281   dummyProj = (TH1F*)Fspec -> ProjectionY("dummy");
0282   dummyProj -> SetBinContent(1,0);
0283   cutoffFreq = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin()); 
0284   
0285 
0286   Fspec_IB = (TH2F*)histsIn -> Get("Fspec_IB");
0287   dummyProj = (TH1F*)Fspec_IB -> ProjectionY("dummy");
0288   dummyProj -> SetBinContent(1,0);
0289 
0290   cutoffFreq_IB = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0291   cutoffFreq_IB_lo = -1*Fspec_IB->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0292                    
0293                    // Fspec_sector = (TH2F*)histsIn -> Get("Fspec_sector");
0294   // dummyProj = (TH1F*)Fspec_sector -> ProjectionY("dummy");
0295   // dummyProj -> SetBinContent(1,0);
0296                    
0297   // cutoffFreq_sector = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0298   // cutoffFreq_sector_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0299 
0300   //get the rest of the histograms we nede  
0301   Fspeci_K = (TH2F*)histsIn -> Get("Fspeci_K");
0302   
0303   Fspeci_SG = (TH2F*)histsIn -> Get("Fspeci_SG");
0304   Fspeci_SG -> SetTitle(";Tower ID; Hit Frequency");
0305   
0306   
0307   Fspeci = (TH2F*)histsIn -> Get("Fspeci");  
0308   Fspeci_IB = (TH2F*)histsIn -> Get("Fspeci_IB");
0309   Fspeci_IB -> SetTitle(";IB ID; Hit Frequency");
0310 
0311   Fspeci_sector = (TH2F*)histsIn -> Get("Fspeci_sector");
0312  
0313  
0314 
0315   //initial hot tower calculation
0316   for(int i = 0; i < nTowers; i++){
0317     
0318     channels->GetEntry(i);
0319 
0320     TH1F *fSpecProj_SG = (TH1F*)Fspeci_SG->ProjectionY("dummySG",i+1,i+1);
0321     TH1F *fSpecProj_K = (TH1F*)Fspeci_K->ProjectionY("dummyK",i+1,i+1);
0322     TH1F *fSpecProj = (TH1F*)Fspeci->ProjectionY("dummySpec",i+1,i+1);
0323     
0324     if(fiber_type == 0 && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) > cutoffFreq_SG){
0325       hottowers[i]++;
0326     }
0327     else if(fiber_type == 1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > cutoffFreq_K){
0328       hottowers[i]++;
0329     }
0330     else if(fSpecProj -> GetBinCenter(fSpecProj->GetMaximumBin()) == 0){
0331       deadtowers[i]++;
0332     }
0333     else if (fiber_type == 0 &&  fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) < cutoffFreq_SG_lo &&  fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) >0 ){
0334       coldtowers[i]++;
0335     }
0336     else if (fiber_type ==1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) < cutoffFreq_K_lo &&  fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > 0){
0337       coldtowers[i]++;  
0338     }
0339     
0340   }
0341   //go through and look for hot/cold "regions", e.g. sectors and interface boards
0342   hot_regions = 0;
0343   cold_regions = 0; 
0344   
0345   //Interface board QA, create projection
0346 
0347   //hotIB
0348   for(int j = 0; j < nIB; j++){
0349     TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
0350 
0351     if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())>cutoffFreq_IB){
0352       hot_regions = 1;
0353       std::cout << "IB " << j << " is hot with ADC rate" << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())<< " > "<< cutoffFreq_IB  << std::endl;
0354       
0355       hotIB[j]++;
0356       for(int j1 = 0; j1 < nTowersIB; j1++ ){
0357     hottowers[j*nTowersIB+j1]++;
0358     for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1, biny+1,-10);
0359     //towerF[j*nTowersIB+j1]=0;//mask tower's contribution so it's no longer used to calculate cutoffs
0360       }
0361       for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB->SetBinContent(j+1,biny+1,-10);
0362       //ibF[j]=-1;
0363     }
0364   }
0365   
0366   //cold IB
0367   for(int j = 0; j < nIB; j++){
0368     TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
0369 
0370     if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())<cutoffFreq_IB_lo){
0371       hot_regions = 1;
0372       std::cout << "IB " << j << " is cold with ADC rate" << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin()) << " < " << cutoffFreq_IB_lo << std::endl;
0373       hotIB[j]++;
0374       
0375       for(int j1 = 0; j1< nTowersIB; j1++ ){
0376     hottowers[j*nTowersIB+j1]++;
0377     for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1,biny+1,-10);
0378     //towerF[j*nTowersIB+j1]=0; 
0379       }
0380       //ibF[j] = -1;
0381       for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB -> SetBinContent(j+1,biny+1,-10);
0382     }
0383   }
0384   
0385   /*Removing sector level qa
0386   Essentially the issue is that removing IB's biases the determination of hot and cold sectors. 
0387   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.
0388   One could, in principle, re-scale the response by the fraction of missing IB's, but it isn't clear that that's necessary. 
0389   Towers are necessarily effected by the behavior of their IB, but it isn't clear that a single sector can be driven 
0390   bad by the behavior of a single IB. - AH*/
0391 
0392 
0393   //QA the sectors
0394   //hot sector
0395   // for(int k = 0; k < nSectors; k++){
0396   //   TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);    
0397     
0398   //   if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) > cutoffFreq_sector){
0399   //     hot_regions = 1;
0400   //     std::cout << "sector " << k << "is hot with ADC rate " << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " > " << cutoffFreq_sector  << std::endl;
0401   //     hotsectors[k]++;
0402     
0403   //     for(int k1 = 0; k1 < nTowersSec; k1++){
0404   //    hottowers[k*nTowersSec+k1]++;
0405   //    towerF[k*nTowersSec+k1] = 0;
0406   //    //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
0407   //     }
0408   //     sectorF[k] = -1;
0409   //     //fSpecSecProj->SetBinContent(k,0);
0410   //     //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
0411   //   }
0412   // }
0413 
0414   
0415 
0416   // //cold sector
0417   // for(int k = 0; k < nSectors; k++){
0418   //   TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);    
0419       
0420   //   if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) < cutoffFreq_sector_lo){
0421   //     cold_regions = 1;
0422   //     std::cout << "sector " << k << "is cold  with ADC rate" << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " < " << cutoffFreq_sector_lo << std::endl;
0423   //     coldsectors[k]++;
0424      
0425   //     for(int k1 = 0; k1 < nIB; k1++){
0426   //    coldtowers[k*nIB+k1]++;
0427   //    //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
0428   //    towerF[k*nTowersSec+k1] = 0;
0429   //     }
0430   //     //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
0431   //     sectorF[k]=-1;
0432   //   }
0433   // }
0434   
0435   std::cout << "removing hot/cold regions" << std::endl;
0436 
0437   //go through again and calculate cut offs with bad regions masked until the calculation is made without hot regions
0438   while(hot_regions == 1 || cold_regions ==1){
0439     std::cout << "hot/cold IB or sector detected. Running another pass for hot towers" << std::endl;
0440     // Fspec->Reset();
0441     // Fspec_SG->Reset();
0442     // Fspec_K->Reset();
0443     // Fspec_IB->Reset();
0444     //Fspec_sector->Reset();
0445 
0446     //not used for now
0447     // Espec->Reset();
0448     // Espec_SG->Reset();
0449     // Espec_K->Reset();
0450     // Espec_IB->Reset();
0451     // Espec_sector->Reset();
0452     //for(int i = 0; i < nTowers; i++){
0453 
0454     //Fspec->Fill(i,towerF[i]/goodevents);
0455       // Espec->Fill(i,towerE[i]/goodevents);
0456       //channels->GetEntry(i);
0457       //if(fiber_type == 0){
0458     //Fspec_SG->Fill(i,towerF[i]/goodevents);
0459     //Espec_SG->Fill(i,towerE[i]/goodevents);
0460     //}
0461     //else{
0462     //Fspec_K->Fill(i,towerF[i]/goodevents);
0463     //Espec_K->Fill(i,towerE[i]/goodevents);
0464     //}
0465     //}
0466     
0467     //for(int j = 0; j<nIB; j++){
0468     //Fspec_IB->Fill(j,1.0*ibF[j]/nTowersIB/goodevents);
0469       //Espec_IB->Fill(1,1.0*ibE[j]/64/goodevents);
0470     //}
0471     //for(int k = 0; k<nSectors; k++){
0472       //Fspec_sector->Fill(k,1.0*sectorF[k]/nTowersSec/goodevents);
0473       //Espec_sector->Fill(1,1.0*sectorE[k]/384/goodevents);
0474     //}
0475 
0476     //kill zeros
0477 
0478     // Fspec->SetBinContent(1,0);
0479     // Fspec_SG->SetBinContent(1,0);
0480     // Fspec_K->SetBinContent(1,0);
0481     // Fspec_IB->SetBinContent(1,0);
0482     // Fspec_sector->SetBinContent(1,0);
0483     
0484     dummyProj = (TH1F*)Fspec_SG -> ProjectionY("dummy");
0485     dummyProj -> SetBinContent(1,0);
0486    
0487 
0488     cutoffFreq_SG = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0489     cutoffFreq_SG_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0490 
0491     dummyProj = (TH1F*)Fspec_K -> ProjectionY("dummy");
0492     dummyProj -> SetBinContent(1,0);
0493 
0494     cutoffFreq_K = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0495     cutoffFreq_K_lo = -1*Fspec_K->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0496 
0497     dummyProj = (TH1F*)Fspec -> ProjectionY("dummy");
0498     dummyProj -> SetBinContent(1,0);
0499     cutoffFreq = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin()); 
0500 
0501     dummyProj = (TH1F*)Fspec_IB -> ProjectionY("dummy");
0502     dummyProj -> SetBinContent(1,0);
0503     cutoffFreq_IB = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0504     cutoffFreq_IB_lo = -1*Fspec_IB->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0505 
0506     // dummyProj = (TH1F*)Fspec_sector -> ProjectionY("dummy");
0507     // dummyProj -> SetBinContent(1,0);
0508     // cutoffFreq_sector = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0509     // cutoffFreq_sector_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0510     
0511        
0512     for(int i = 0; i < nTowers; i++){
0513       
0514       channels->GetEntry(i);
0515       TH1F *fSpecProj_SG = (TH1F*)Fspeci_SG->ProjectionY("dummySG",i+1,i+1);
0516       TH1F *fSpecProj_K = (TH1F*)Fspeci_K->ProjectionY("dummyK",i+1,i+1);
0517       TH1F *fSpecProj = (TH1F*)Fspeci->ProjectionY("dummySpec",i+1,i+1);
0518 
0519       if(fiber_type == 0 && fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) > cutoffFreq_SG){
0520     hottowers[i]++;
0521       }
0522       else if(fiber_type == 1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > cutoffFreq_K){
0523     hottowers[i]++;
0524       }
0525       else if(fSpecProj -> GetBinCenter(fSpecProj->GetMaximumBin()) ==0){
0526     deadtowers[i]++;
0527       }
0528       else if (fiber_type == 0 &&  fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) < cutoffFreq_SG_lo &&  fSpecProj_SG->GetBinCenter(fSpecProj_SG->GetMaximumBin()) >0 ){
0529     coldtowers[i]++;
0530       }
0531       else if (fiber_type ==1 && fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) < cutoffFreq_K_lo &&  fSpecProj_K->GetBinCenter(fSpecProj_K->GetMaximumBin()) > 0){
0532     coldtowers[i]++;    
0533       }
0534     }
0535 
0536     hot_regions = 0;
0537     cold_regions = 0;
0538    
0539     //Re-QA IB's
0540     //hotIB
0541     for(int j = 0; j< nIB; j++){
0542       if(ibF[j]==-1)continue; //IB has been marked hot already, don't bother with it
0543       TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
0544       
0545       if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())>cutoffFreq_IB){
0546     hot_regions = 1;
0547     std::cout << "IB " << j << "is hot with ADC rate " << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin()) << " > " << cutoffFreq_IB << std::endl;
0548     hotIB[j]++;
0549 
0550     for(int j1 = 0; j1 < nTowersIB; j1++ ){
0551       hottowers[j*nTowersIB+j1]++;
0552       //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1,biny+1, 0);
0553       towerF[j*nTowersIB+j1] = 0;//mask tower's contribution so it's no longer used to calculate cutoffs
0554     }
0555     //for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB->SetBinContent(j+1,biny+1,0);
0556     ibF[j] = -1;
0557       }
0558     }
0559   
0560     //cold IB
0561     for(int j = 0; j<nIB; j++){
0562       if(ibF[j]==-1)continue; //IB has been marked hot already, don't bother with it
0563 
0564       TH1F *fSpecIBProj = (TH1F*)Fspeci_IB->ProjectionY("dummyIB",j+1,j+1);
0565       
0566       if(fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin())<cutoffFreq_IB_lo){
0567     hot_regions = 1;
0568     std::cout << "IB " << j << " is cold with ADC rate " << fSpecIBProj->GetBinCenter(fSpecIBProj->GetMaximumBin()) << " < " << cutoffFreq_IB_lo << std::endl;
0569     hotIB[j]++;
0570 
0571     for(int j1 = 0; j1< nTowersIB; j1++ ){
0572       hottowers[j*nTowersIB+j1]++;
0573       //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci -> SetBinContent(j*nTowersIB+j1+1,biny+1,0);
0574       towerF[j*nTowersIB+j1] = 0;   
0575     }
0576     ibF[j] = -1;
0577     //for(int biny = 0; biny < Fspeci_IB->GetNbinsY(); biny++)Fspeci_IB->SetBinContent(j+1, biny+1,0);
0578       }
0579     }
0580     
0581     //Re-QA sectors
0582     //hot sectors
0583     // for(int k = 0; k<nSectors; k++){
0584     //   if(sectorF[k]==-1) continue;//sector has been marked hot already
0585     //   TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);    
0586       
0587     //   if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) > cutoffFreq_sector){
0588     //  hot_regions = 1;
0589     //  std::cout << "sector " << k << " is hot with ADC rate " << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) << " > " <<cutoffFreq_sector  << std::endl;
0590     //  hotsectors[k]++;
0591     
0592     //  for(int k1 = 0; k1 < nTowersSec; k1++){
0593     //    hottowers[k*nTowersSec+k1]++;
0594      
0595     //    //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
0596     //    towerF[k*nTowersSec+k1] = 0;
0597     //  }
0598     //  sectorF[k] = -1;
0599     //  //fSpecSecProj->SetBinContent(k,0);
0600     //  //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
0601     //   }
0602     // }
0603 
0604     //cold sector
0605     // for(int k = 0; k<nSectors; k++){
0606     //   if(sectorF[k]==-1) continue;
0607     //   TH1F *fSpecSecProj = (TH1F*)Fspeci_sector->ProjectionY("dummySec",k+1,k+1);    
0608 
0609     //   if(fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin()) <cutoffFreq_sector_lo){
0610     //  cold_regions = 1;
0611     //  std::cout << "sector " << k << " is cold  with ADC rate " << fSpecSecProj->GetBinCenter(fSpecSecProj->GetMaximumBin())  << " < " <<cutoffFreq_sector_lo << std::endl;
0612     //  coldsectors[k]++;
0613     
0614     //  for(int k1 = 0; k1 < nIB; k1++){
0615     //    coldtowers[k*nIB+k1]++;
0616     //    //for(int biny = 0; biny < Fspeci->GetNbinsY(); biny++) Fspeci->SetBinContent(k*nTowersSec+k1, biny+1, 0);
0617     //     towerF[k*nTowersSec+k1] = 0;
0618     //  }
0619     //  //for(int biny = 0; biny < Fspeci_sector->GetNbinsY(); biny++) Fspeci_sector->SetBinContent(k+1, biny+1, 0);
0620     //  sectorF[k] = -1;
0621     //   }
0622     // } 
0623     
0624   }
0625   
0626   std::cout << "towerid::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0627   std::cout << "Saint Gobain Cutoff Frequency: " << cutoffFreq_SG << std::endl;
0628   std::cout <<  "Kurary Cutoff Frequency: " << cutoffFreq_K << std::endl;
0629   std::cout << "Overall Tower Cutoff Frequency: " << cutoffFreq << std::endl;
0630   std::cout << "IB Cutoff hit Frequency: " << cutoffFreq_IB << std::endl;
0631   //std::cout <<  "Sector Cutoff hit Frequency: " << cutoffFreq_sector << std::endl;
0632 }
0633 
0634 //____________________________________________________________________________..
0635 void towerid::WriteCDBTree(const int runnumber)
0636 {
0637   
0638   
0639   //TFile *cdbOut = new TFile(Form("cdbMaps/%d/DST_CALOR-%08d_badTowerMapCDBTTree.root",runnumber,runnumber),"RECREATE");
0640   TFile *cdbOut = new TFile(Form("cdbMaps/%d/CEMC-%08d_badTowerMapCDBTTree.root",runnumber,runnumber),"RECREATE");
0641   
0642   T = new TTree("T","T");
0643   T -> Branch("hot_channels",&m_hot_channels);
0644   
0645   CDBTTree*cdbttree = new CDBTTree(cdbtreename);
0646   std::string fieldname = "status";
0647   for(int i = 0; i<nTowers; i++){
0648 
0649     unsigned int key = TowerInfoDefs::encode_emcal(i);
0650     
0651     if(hottowers[i] >= 0.5){
0652       m_hot_channels = 2;
0653       T->Fill();
0654       cdbttree->SetIntValue(key,fieldname,2);
0655     }
0656     else if(deadtowers[i] >= 0.5){
0657       m_hot_channels = 1;
0658       T->Fill();
0659       cdbttree->SetIntValue(key,fieldname,1);
0660     }
0661     else if(coldtowers[i] > 0.5){
0662       m_hot_channels = 3;
0663       T->Fill();
0664       cdbttree->SetIntValue(key,fieldname,3);
0665     }
0666     else{
0667       m_hot_channels = 0;
0668       T->Fill();
0669       cdbttree->SetIntValue(key,fieldname,0);
0670     }
0671   }
0672 
0673   fchannels->Close();
0674   delete fchannels;
0675   
0676   cdbOut->cd();
0677   T -> Write();
0678   cdbttree->Commit();
0679   //cdbttree->Print();
0680   cdbttree->WriteCDBTTree();
0681   delete cdbttree;
0682   cdbOut->Close();
0683 
0684 }
0685 
0686 //____________________________________________________________________________..
0687 int towerid::Reset(PHCompositeNode *topNode)
0688 {
0689   std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0690   return Fun4AllReturnCodes::EVENT_OK;
0691 }
0692 //____________________________________________________
0693 void towerid::Print(const std::string &what) const
0694 {
0695   std::cout << "towerid::Print(const std::string &what) const Printing info for " << what << std::endl;
0696 }
0697 //____________________________________________________
0698 int towerid::End(PHCompositeNode *topNode)
0699 {
0700   std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0701   return Fun4AllReturnCodes::EVENT_OK;
0702 }
0703 //____________________________________________________
0704 int towerid::EndRun(int runnumber)
0705 {
0706   std::cout << "towerid::EndRun: this is the end of run: "<< runnumber << std::endl;
0707   return Fun4AllReturnCodes::EVENT_OK;
0708 }
0709 //____________________________________________________
0710 
0711 
0712 
0713 
0714 
0715 
0716
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!