epic-lfhcal-tbana/NewStructure/TileSpectra.cc

0001 #include "TileSpectra.h"
0002 #include "TFitResult.h"
0003 #include "TFitResultPtr.h"
0004
0005 ClassImp(TileSpectra);
0006
0007 int TileSpectra::GetCellID(){
0008   return cellID;
0009 }
0010
0011 bool TileSpectra::FillCAEN(double l, double h){
0012   hspectraLG.Fill(l);
0013   hspectraHG.Fill(h);
0014   hspectraLGHG.Fill(l,h);
0015   hspectraHGLG.Fill(h,l);
0016   return true;
0017 }
0018
0019 bool TileSpectra::FillHGCROC(double adc, double toa, double tot){
0020   hspectraHG.Fill(adc);
0021   hspectraTOT.Fill(toa);
0022   hspectraTOA.Fill(tot);
0023   hADCTOT.Fill(adc, tot);
0024   hTOAADC.Fill(toa, adc);
0025   return true;
0026 }
0027
0028 bool TileSpectra::FillSpectraCAEN(double l, double h){
0029   hspectraLG.Fill(l);
0030   hspectraHG.Fill(h);
0031   return true;
0032 }
0033
0034 bool TileSpectra::FillSpectraHGCROC(double adc, double toa, double tot){
0035   hspectraHG.Fill(adc);
0036   hspectraTOT.Fill(toa);
0037   hspectraTOA.Fill(tot);
0038   return true;
0039 }
0040
0041 bool TileSpectra::FillExtCAEN(double l, double h, double e, double lheq){
0042
0043   if (ROType != ReadOut::Type::Caen){
0044     std::cout << "\n\n ******************************************************** \n\n" << std::endl;
0045     std::cout << "ERROR:: You are using a CAEN filling function to fill HGCROC hists! Aborting!" << std::endl;
0046     std::cout << "\n\n ******************************************************** \n\n" << std::endl;
0047     return false;
0048   }
0049   hspectraLG.Fill(l);
0050   hspectraHG.Fill(h);
0051   if (extend ==  1){
0052     hcombined.Fill(e);
0053     if (h < 3500)
0054       hspectraLGHG.Fill(l,e);
0055     hspectraHGLG.Fill(l,lheq-h);
0056   } else if (extend ==  2){
0057     hspectraLGHG.Fill(l,h);
0058     hcorr.Fill(l,h);
0059   }
0060   return true;
0061 }
0062
0063 bool TileSpectra::FillExtHGCROC(double adc = 0., double toa= 0., double tot= 0., int sample = 0){
0064   if (ROType != ReadOut::Type::Hgcroc){
0065     std::cout << "\n\n ******************************************************** \n\n" << std::endl;
0066     std::cout << "ERROR:: You are using a HGCROC filling function to fill CAEN hists! Aborting!" << std::endl;
0067     std::cout << "\n\n ******************************************************** \n\n" << std::endl;
0068     return false;
0069   }
0070
0071   hspectraHG.Fill(adc);
0072   hspectraTOA.Fill(toa);
0073   hspectraTOT.Fill(tot);
0074   if (extend ==  2 ){
0075     hADCTOT.Fill(adc,tot);
0076   }
0077   if (extend == 3){
0078     hcorrTOAADC.Fill(toa,adc);
0079     hTOAADC.Fill(toa,adc);
0080     hcorrTOASample.Fill(toa,sample);
0081   }
0082   return true;
0083 }
0084
0085
0086 bool TileSpectra::FillExtHGCROCPed(std::vector<int> samples, double h){
0087   // fill pedestal only (1st sample of waveform)
0088   hspectraHG.Fill(h);
0089   // fill all samples of waveform
0090   for (int k = 0; k < (int)samples.size(); k++ ){
0091    hspectraLG.Fill(samples.at(k));
0092   }
0093   // fill correlation sample vs ADC
0094   for (int k = 0; k < (int)samples.size(); k++ ){
0095     hcorr.Fill(k,samples.at(k));
0096   }
0097   return true;
0098 }
0099
0100 bool TileSpectra::FillWaveform(std::vector<int> samples, double ped = 0){
0101  for (int k = 0; k < (int)samples.size(); k++ ){
0102    hcorr.Fill(k,samples.at(k)-ped);
0103    if (extend == 3) hWaveForm.Fill(k,samples.at(k)-ped);
0104  }
0105  return true;
0106 }
0107
0108 bool TileSpectra::FillWaveformVsTime(std::vector<int> samples, double toa = 0, double ped = 0, int offset = 0){
0109   // more infos
0110   for (int k = 0; k < (int)samples.size(); k++ ){
0111     hcorr.Fill((k-offset+1)*25000-toa*25,samples.at(k)-ped);
0112     if (extend == 3) hWaveForm.Fill((k-offset+1)*25000-toa*25,samples.at(k)-ped);
0113   }
0114   return true;
0115 }
0116
0117
0118 bool TileSpectra::FillTrigger(double t){
0119   if (!bTriggPrim) bTriggPrim =true;
0120   hTriggPrim.Fill(t);
0121   return true;
0122 }
0123
0124
0125 bool TileSpectra::FillCorrCAEN(double l, double h){
0126   hspectraLGHG.Fill(l,h);
0127   hspectraHGLG.Fill(h,l);
0128   return true;
0129 }
0130
0131 bool TileSpectra::FillCorrHGCROC(double adc, double toa, double tot){
0132   hADCTOT.Fill(adc, tot);
0133   hTOAADC.Fill(toa, adc);
0134   return true;
0135 }
0136
0137
0138 bool TileSpectra::FitNoise(double* out, int year = -1, bool isNoiseTrigg = false){        //[0] LG mean, [2] LG sigma, [4] HG mean, [6] HG sigma errors uneven numbers
0139   TFitResultPtr result;
0140   if (ROType == ReadOut::Type::Caen) {
0141     // estimate LG pedestal per channel
0142     BackgroundLG=TF1(Form("fped%sLGCellID%d",TileName.Data(),cellID),"gaus",0,400);
0143     BackgroundLG.SetNpx(400);
0144     if (year == 2023){
0145       BackgroundLG.SetParameter(1,50);
0146       BackgroundLG.SetParLimits(1,40,60);     // might need to make these values settable
0147       BackgroundLG.SetParameter(2,4);
0148       BackgroundLG.SetParLimits(2,0,10);     // might need to make these values settable
0149       BackgroundLG.SetRange(0,70);
0150     } else {
0151       BackgroundLG.SetParameter(1,hspectraLG.GetMean());
0152       BackgroundLG.SetParLimits(1,0,hspectraLG.GetMean()+100);     // might need to make these values settable
0153       BackgroundLG.SetParameter(2,10);
0154       BackgroundLG.SetParLimits(2,0,100);     // might need to make these values settable
0155     }
0156     double maxLG = 0;
0157     if (isNoiseTrigg){
0158       maxLG = GetMaxXInRangeLG(0,300);
0159       BackgroundLG.SetParameter(1,maxLG);
0160       BackgroundLG.SetParLimits(1,maxLG-5,maxLG+5);
0161       if (debug > 1) std::cout << "reset LG: " << maxLG << std::endl;
0162     }
0163
0164     if (!isNoiseTrigg){
0165       BackgroundLG.SetParLimits(0,0,hspectraLG.GetEntries());
0166       BackgroundLG.SetParameter(0,hspectraLG.GetEntries()/5);
0167       result=hspectraLG.Fit(&BackgroundLG,"QRMEN0S"); // initial fit
0168       double minLGFit = result->Parameter(1)-2*result->Parameter(2);
0169       double maxLGFit = result->Parameter(1)+1*result->Parameter(2);
0170       if (debug > 1) std::cout << "LG: " << minLGFit << "\t" << maxLGFit << "\t" << hspectraLG.GetEntries() << "\t" << hspectraLG.GetMean()<< std::endl;
0171       result=hspectraLG.Fit(&BackgroundLG,"QRMEN0S","", minLGFit, maxLGFit);  // limit to 2sigma
0172     } else {
0173       result=hspectraLG.Fit(&BackgroundLG,"QRMEN0S","", maxLG-20, maxLG+30);  // limit to 2sigma
0174       if (debug > 1) std::cout <<"LG: " << result->Parameter(0) << "\t"<< result->Parameter(1) << "\t" << result->Parameter(2) << std::endl;
0175     }
0176     bpedLG=true;
0177     calib->PedestalMeanL=result->Parameter(1);//Or maybe we do not want to do it automatically, only if =0?
0178     calib->PedestalSigL =result->Parameter(2);//Or maybe we do not want to do it automatically, only if =0?
0179     out[0]=result->Parameter(1);
0180     out[1]=result->Error(1);
0181     out[2]=result->Parameter(2);
0182     out[3]=result->Error(2);
0183
0184     // estimate HG pedestal per channel
0185     BackgroundHG=TF1(Form("fped%sHGCellID%d",TileName.Data(),cellID),"gaus",0,400);
0186     BackgroundHG.SetNpx(400);
0187     if (year == 2023){
0188       BackgroundHG.SetParameter(1,60);
0189       BackgroundHG.SetParLimits(1,0,100);     // might need to make these values settable
0190       BackgroundHG.SetRange(0,200);
0191     } else {
0192       BackgroundHG.SetParameter(1,hspectraHG.GetMean());
0193       BackgroundHG.SetParLimits(1,0,hspectraHG.GetMean()+100);     // might need to make these values settable
0194     }
0195     BackgroundHG.SetParameter(2,10);
0196     BackgroundHG.SetParLimits(2,0,100);     // might need to make these values settable
0197
0198     double maxHG = 0;
0199     if (isNoiseTrigg){
0200       maxHG = GetMaxXInRangeHG(0,300);
0201       BackgroundHG.SetParameter(2,20);
0202       BackgroundHG.SetParameter(1,maxHG);
0203       BackgroundHG.SetParLimits(1,maxHG-5,maxHG+5);
0204       if (debug > 1) std::cout << "reset HG: " << maxHG << std::endl;
0205     }
0206
0207     if (!isNoiseTrigg){
0208       BackgroundHG.SetParLimits(0,0,hspectraHG.GetEntries());
0209       BackgroundHG.SetParameter(0,hspectraHG.GetEntries()/10);
0210       result=hspectraHG.Fit(&BackgroundHG,"QRMEN0S");      // initial fit
0211       double minHGFit = result->Parameter(1)-2*result->Parameter(2);
0212       double maxHGFit = result->Parameter(1)+1*result->Parameter(2);
0213       if (debug > 1) std::cout <<"HG: " << minHGFit << "\t" << maxHGFit << "\t" << hspectraHG.GetEntries() << "\t" << hspectraHG.GetMean()<< std::endl;
0214       result=hspectraHG.Fit(&BackgroundHG,"QRMEN0S","",minHGFit, maxHGFit);  // limit to 2sigma range of previous fit
0215     } else {
0216       result=hspectraHG.Fit(&BackgroundHG,"QRMEN0S","",maxHG-40, maxHG+60);
0217       if (debug > 1) std::cout <<"HG: " << result->Parameter(0) << "\t"<< result->Parameter(1) << "\t" << result->Parameter(2) << std::endl;
0218     }
0219     bpedHG=true;
0220
0221     calib->PedestalMeanH=result->Parameter(1);//Or maybe we do not want to do it automatically, only if =0?
0222     calib->PedestalSigH =result->Parameter(2);//Or maybe we do not want to do it automatically, only if =0?
0223     out[4]=result->Parameter(1);
0224     out[5]=result->Error(1);
0225     out[6]=result->Parameter(2);
0226     out[7]=result->Error(2);
0227
0228     return true;
0229   } else if (ROType == ReadOut::Type::Hgcroc) {
0230     /*
0231     QRMEN0S
0232     Q: Quiet mode
0233     R: Use TF1::SetRange
0234     M: IMPROVE algorithm
0235     E: Minos error estimation
0236     N: Don't store graphics and don't draw
0237     0: Don't draw, but do store?
0238     S: Return full result in TFitResultPtr
0239     QNSWW
0240     */
0241     // hspectraHG.Rebin(2);
0242     BackgroundHG = TF1(Form("fped%sHGCellID%d", TileName.Data(), cellID), "gaus", 0, 400);
0243     if (debug > 2) {
0244       std::cout << "Histogram has " << hspectraHG.GetEntries() << " entries" << std::endl;
0245       std::cout << "Mean is " << hspectraHG.GetMean() << std::endl;
0246       std::cout << "Standard deviation is " << hspectraHG.GetStdDev() << std::endl;
0247     }
0248
0249     // First iteration
0250     // BackgroundHG.SetParLimits(0, 0, hspectraHG.GetEntries());
0251     // BackgroundHG.SetParameter(0, hspectraHG.GetEntries() / 5);
0252     // BackgroundHG.SetParameter(1, hspectraHG.GetMean());
0253     // BackgroundHG.SetParLimits(1, 0, hspectraHG.GetMean() + 100);
0254     // BackgroundHG.SetParameter(2, hspectraHG.GetStdDev());
0255     // BackgroundHG.SetParLimits(2, 0, 100);
0256
0257     result = hspectraHG.Fit(&BackgroundHG, "QNSWW");
0258     if ((int)result != 0 || result->IsValid() != true) {
0259       if (debug > 1) std::cout << "FIT FAILED FOR CELL " << cellID << ", FIRST ITERATION" << std::endl;
0260       return false;
0261     }
0262
0263     // Second iteration
0264     double minHGFit = result->Parameter(1) - 2 * result->Parameter(2);
0265     double maxHGFit = result->Parameter(1) + 1 * result->Parameter(2);
0266     if (debug > 1) std::cout << "HG: " << minHGFit << "\t" << maxHGFit << "\t" << hspectraHG.GetEntries() << "\t" << hspectraHG.GetMean()<< std::endl;
0267     result = hspectraHG.Fit(&BackgroundHG, "QNSWW", "", minHGFit, maxHGFit);  // limit to 2sigma
0268     if ((int)result != 0 || result->IsValid() != true){
0269     if (debug > 1) std::cout << "FIT FAILED FOR CELL " << cellID << ", SECOND ITERATION" << std::endl;
0270       return false;
0271     }
0272
0273     bpedHG=true; // We're (I'm) being lazy and just calling the HGCROC ADC the high gain.  maybe we use LG for TOT info?
0274     calib->PedestalMeanH=result->Parameter(1);
0275     calib->PedestalSigH =result->Parameter(2);
0276     out[4]=result->Parameter(1);
0277     out[5]=result->Error(1);
0278     out[6]=result->Parameter(2);
0279     out[7]=result->Error(2);
0280     return true;
0281   }
0282
0283   return false;
0284 }
0285
0286 void TileSpectra::FitFixedNoise(){
0287
0288   // estimate LG pedestal per channel
0289   BackgroundLG=TF1(Form("fpedLGCellID%d",cellID),"gaus",-10,10);
0290   BackgroundLG.SetNpx(400);
0291   BackgroundLG.SetParLimits(1,-2,2);
0292   BackgroundLG.SetParameter(2,calib->PedestalSigL);
0293
0294   hspectraLG.Fit(&BackgroundLG,"QIRN0"); // initial fit
0295   bpedLG=BackgroundLG.IsValid();
0296
0297   // estimate HG pedestal per channel
0298   BackgroundHG=TF1(Form("fpedHGCellID%d",cellID),"gaus",-20,20);
0299   BackgroundHG.SetNpx(400);
0300   // BackgroundHG.FixParameter(1,0);
0301   BackgroundLG.SetParLimits(1,-2,2);
0302   BackgroundHG.SetParameter(2,calib->PedestalSigH);
0303
0304   hspectraHG.Fit(&BackgroundHG,"QIRN0"); // initial fit
0305   bpedHG=BackgroundHG.IsValid();
0306
0307   return;
0308 }
0309
0310 void TileSpectra::InitializeNoiseFitsFromCalib(){
0311
0312   // estimate LG pedestal per channel
0313   BackgroundLG=TF1(Form("fpedLGCellID%d",cellID),"gaus",-20,20);
0314   BackgroundLG.SetNpx(400);
0315   BackgroundLG.FixParameter(1,0);
0316   BackgroundLG.FixParameter(2,calib->PedestalSigL);
0317   hspectraLG.Fit(&BackgroundLG,"QIRN0"); // initial fit
0318   BackgroundLG.FixParameter(1,0);
0319   BackgroundLG.FixParameter(2,calib->PedestalSigL);
0320   bpedLG=true;
0321
0322   // estimate HG pedestal per channel
0323   BackgroundHG=TF1(Form("fpedHGCellID%d",cellID),"gaus",-30,30);
0324   BackgroundHG.SetNpx(400);
0325   BackgroundHG.FixParameter(1,0);
0326   BackgroundHG.FixParameter(2,calib->PedestalSigH);
0327   hspectraHG.Fit(&BackgroundHG,"QIRN0"); // initial fit
0328   BackgroundHG.FixParameter(1,0);
0329   BackgroundHG.FixParameter(2,calib->PedestalSigH);
0330   bpedHG=true;
0331
0332   return;
0333 }
0334
0335
0336 bool TileSpectra::FitMipHG(double* out, double* outErr, int verbosity, int year, bool impE = false, double vov = -1000, double avmip = -1000){
0337
0338   // Once again, here are the Landau * Gaussian parameters:
0339   //   par[0]=Width (scale) parameter of Landau density
0340   //   par[1]=Most Probable (MP, location) parameter of Landau density
0341   //   par[2]=Total area (integral -inf to inf, normalization constant)
0342   //   par[3]=Width (sigma) of convoluted Gaussian function
0343   //
0344
0345   if (verbosity > 2) std::cout << "FitHG cell ID: " << cellID << std::endl;
0346
0347   TString funcName = Form("fmip%sHGCellID%d",TileName.Data(),cellID);
0348   bmipHG           = false;
0349
0350   if (calib->BadChannel != -64 && calib->BadChannel < 1 ){
0351     if (verbosity > 0) std::cout << "==========> Skipped HG cell " << cellID << " channel dead" << std::endl;
0352     return false;
0353   }
0354
0355   double fitrange[2]      = {50, 2000};
0356   if (ROType == ReadOut::Type::Hgcroc){
0357     // hspectraHG.Rebin(8);
0358     fitrange[0] = 16;
0359     fitrange[1] = 500;
0360   }
0361   if (impE){
0362     fitrange[0] = 0.6*avmip;
0363     fitrange[1] = 3*avmip;
0364   }
0365   if (year == 2023 && fitrange[0] < 100)
0366     fitrange[0] = 100;
0367
0368   double intArea    = hspectraHG.Integral(hspectraHG.FindBin(fitrange[0]),hspectraHG.FindBin(fitrange[1]));
0369   double intNoise   = hspectraHG.Integral(hspectraHG.FindBin(-2*calib->PedestalSigH),hspectraHG.FindBin(+2*calib->PedestalSigH));
0370   double intAN3s    = hspectraHG.Integral(hspectraHG.FindBin(+3*calib->PedestalSigH),hspectraHG.FindBin(fitrange[1]));
0371
0372   if (intArea/intNoise < 1e-5 && intAN3s > 200){
0373     if (verbosity > 0) std::cout << "==========> Skipped HG cell " << cellID << " S/B too small!" << std::endl;
0374     return false;
0375   }
0376   double startvalues[4]   = {50, 300, intArea, calib->PedestalSigH};
0377   double parlimitslo[4]   = {0.5, 50, 1.0, calib->PedestalSigH*0.01};
0378   double parlimitshi[4]   = {500, 1000, intArea*5, calib->PedestalSigH*40};
0379   if (year == 2023){
0380     startvalues[0]  = 200;
0381     startvalues[1]  = 500;
0382     parlimitslo[1]  = 100;
0383     parlimitshi[0]  = 1000;
0384     parlimitshi[1]  = 1500;
0385   } else if (ROType == ReadOut::Type::Hgcroc){
0386     parlimitslo[1]  = 20;
0387     fitrange[0]     = 15;
0388   }
0389
0390   if (impE && (avmip =! -1000)){
0391     startvalues[1]  = avmip;
0392     parlimitslo[1]  = 0.5*avmip;
0393     parlimitshi[1]  = 1.7*avmip;
0394   }
0395   if (ROType == ReadOut::Type::Caen){
0396     if (vov != -1000){
0397       if (verbosity > 1) std::cout << "adjusting according to V_ov: " << vov<< std::endl;
0398       if (vov < 2.5){
0399         parlimitslo[1]  = 20;
0400         fitrange[0]     = 15;
0401       } else if (vov > 5 ){
0402         fitrange[0]     = 150;
0403       }
0404     }
0405   }
0406
0407   if (verbosity > 1) {
0408     std::cout << "Fit range: " << fitrange[0] << "\t" << fitrange[1] << std::endl;
0409     for (int i=0; i<4; i++) {
0410       std::cout << "parameter " << i << ": " << startvalues[i] << "\t" << fitrange[0] << "\t" << fitrange[1] << std::endl;
0411     }
0412   }
0413
0414   SignalHG = TF1(funcName.Data(),langaufun,fitrange[0],fitrange[1],4);
0415   SignalHG.SetNpx(1000);
0416   SignalHG.SetParameters(startvalues);
0417   SignalHG.SetParNames("Width","MP","Area","GSigma");
0418
0419   for (int i=0; i<4; i++) {
0420     SignalHG.SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0421   }
0422
0423   TString fitOption = "";
0424   if (impE){
0425     fitOption = "QRLMNE0";
0426     if (verbosity > 2)
0427       fitOption = "RLMNE0";
0428   } else {
0429     fitOption = "QRLN0";
0430     if (verbosity > 2)
0431       fitOption = "RLN0";
0432   }
0433
0434   int fitStatus = hspectraHG.Fit(&SignalHG,fitOption);   // fit within specified range, use ParLimits, do not plot
0435   // Minuit status codes:
0436   // 4000 - Successful
0437   //    0 - Successful
0438   // 4070 - Problems
0439   //   70 - Problems
0440   //
0441   if (!SignalHG.IsValid())
0442     return false;
0443
0444   int limitStatus = 0;
0445   for (int i=0; i<4; i++) {
0446     if ( TMath::Abs(SignalHG.GetParameter(i) - parlimitslo[i]) < 1e-5 || TMath::Abs(SignalHG.GetParameter(i) - parlimitshi[i]) < 1e-5 ) {
0447       limitStatus++;
0448       if (verbosity > 0) std::cout << i << "\t" << SignalHG.GetParameter(i) << "\t : \t"<< parlimitslo[i] << "\t" << parlimitshi[i] << std::endl;
0449     }
0450   }
0451   if (verbosity > 1){
0452     std::cout << "Fit status HG " << cellID << " \t" << fitStatus << "\t limit reached: " << limitStatus  << std::endl;
0453   }
0454
0455   if (!(fitStatus == 4000 || fitStatus == 0 || fitStatus == 4070 || fitStatus == 70 )){ // only accept fits which succeeded in general
0456     if (verbosity > 0) std::cout << "==========> Skipped HG cell " << cellID << " fit failed" << std::endl;
0457     return false;
0458   }
0459   if (limitStatus > 0){                        // don't accept fits which reached the set limits
0460     if (verbosity > 0) std::cout << "==========> Skipped HG cell " << cellID << " too many limits reached" << std::endl;
0461     return false;
0462   }
0463   bmipHG = true;
0464
0465   if (bmipHG){
0466     SignalHG.GetParameters(out);    // obtain fit parameters
0467     for (int i=0; i<4; i++) {
0468       outErr[i] = SignalHG.GetParError(i);     // obtain fit parameter errors
0469     }
0470     outErr[4] = SignalHG.GetChisquare();  // obtain chi^2
0471     outErr[5] = SignalHG.GetNDF();           // obtain ndf
0472
0473     double SNRPeak, SNRFWHM;
0474     langaupro(out,SNRPeak,SNRFWHM);
0475
0476     calib->ScaleH       = SNRPeak;
0477     calib->ScaleWidthH  = SNRFWHM;
0478     out[4]    = SNRPeak;
0479     out[5]    = SNRFWHM;
0480   }
0481   return bmipHG;
0482 }
0483
0484
0485 bool TileSpectra::FitMipLG(double* out, double* outErr, int verbosity, int year, bool impE = false, double avmip = 1){
0486
0487   // Once again, here are the Landau * Gaussian parameters:
0488   //   par[0]=Width (scale) parameter of Landau density
0489   //   par[1]=Most Probable (MP, location) parameter of Landau density
0490   //   par[2]=Total area (integral -inf to inf, normalization constant)
0491   //   par[3]=Width (sigma) of convoluted Gaussian function
0492   //
0493
0494   TString funcName = Form("fmip%sLGCellID%d",TileName.Data(),cellID);
0495
0496   double fitrange[2]      = {0, 500};
0497   if (impE){
0498     fitrange[0] = 0.5*avmip;
0499     fitrange[1] = 4*avmip;
0500   }
0501
0502   if (calib->BadChannel != -64 && calib->BadChannel < 1 ){
0503     if (verbosity > 0) std::cout << "==========> Skipped LG cell " << cellID << " channel dead" << std::endl;
0504     return false;
0505   }
0506
0507   double intArea    = hspectraLG.Integral(hspectraLG.FindBin(fitrange[0]),hspectraLG.FindBin(fitrange[1]));
0508   double intNoise   = hspectraLG.Integral(hspectraLG.FindBin(-2*calib->PedestalSigL),hspectraLG.FindBin(+2*calib->PedestalSigL));
0509   double intAN3s    = hspectraLG.Integral(hspectraLG.FindBin(+3*calib->PedestalSigL),hspectraLG.FindBin(fitrange[1]));
0510
0511   if (intArea/intNoise < 1e-5 && intAN3s > 200){
0512     if (verbosity > 0) std::cout << "==========> Skipped LG cell " << cellID << " S/B too small!" << std::endl;
0513     return false;
0514   }
0515   double startvalues[4]   = {calib->PedestalSigL, 20, intArea, calib->PedestalSigL};
0516   double parlimitslo[4]   = {0.5, 0, 1.0, calib->PedestalSigL*0.1};
0517   double parlimitshi[4]   = {calib->PedestalSigL*10, 600, intArea*5, calib->PedestalSigL*10};
0518
0519   SignalLG = TF1(funcName.Data(),langaufun,fitrange[0],fitrange[1],4);
0520   SignalLG.SetNpx(1000);
0521   SignalLG.SetParameters(startvalues);
0522   SignalLG.SetParNames("Width","MP","Area","GSigma");
0523
0524   for (int i=0; i<4; i++) {
0525     SignalLG.SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0526   }
0527
0528   TString fitOption = "";
0529   if (impE){
0530     fitOption = "QRLMNE0";
0531     if (verbosity > 1)
0532       fitOption = "RLMNE0";
0533   } else {
0534     fitOption = "QRLN0";
0535     if (verbosity > 1)
0536       fitOption = "RLN0";
0537   }
0538
0539   int fitStatus = hspectraLG.Fit(&SignalLG,fitOption);   // fit within specified range, use ParLimits, do not plot
0540   // Minuit status codes:
0541   // 4000 - Successful
0542   //    0 - Successful
0543   // 4070 - Problems
0544   //   70 - Problems
0545   //
0546   if (!SignalLG.IsValid())
0547     return false;
0548
0549   int limitStatus = 0;
0550   for (int i=0; i<4; i++) {
0551     if ( TMath::Abs(SignalLG.GetParameter(i) - parlimitslo[i]) < 1e-5 || TMath::Abs(SignalLG.GetParameter(i) - parlimitshi[i]) < 1e-5 ) {
0552       limitStatus++;
0553       if (verbosity > 0) std::cout << i << "\t" << SignalLG.GetParameter(i) << "\t : \t"<< parlimitslo[i] << "\t" << parlimitshi[i] << std::endl;
0554     }
0555   }
0556   if (verbosity > 1){
0557     std::cout << "Fit status LG " << cellID << " \t" << fitStatus << "\t limit reached: " << limitStatus  << std::endl;
0558   }
0559
0560   if (!(fitStatus == 4000 || fitStatus == 0 || fitStatus == 4070 || fitStatus == 70)){ // only accept fits which succeeded or problems in general
0561     if (verbosity > 0) std::cout << "==========> Skipped LG cell " << cellID << " fit failed" << std::endl;
0562     return false;
0563   }
0564   if (limitStatus > 0){                        // don't accept fits which reached the set limits
0565     if (verbosity > 0) std::cout << "==========> Skipped LG cell " << cellID << " too many limits reached" << std::endl;
0566     return false;
0567   }
0568   bmipLG = true;
0569
0570   if (bmipLG){
0571     SignalLG.GetParameters(out);    // obtain fit parameters
0572     for (int i=0; i<4; i++) {
0573       outErr[i] = SignalLG.GetParError(i);     // obtain fit parameter errors
0574     }
0575     outErr[4] = SignalLG.GetChisquare();  // obtain chi^2
0576     outErr[5] = SignalLG.GetNDF();           // obtain ndf
0577
0578     double SNRPeak, SNRFWHM;
0579     langaupro(out,SNRPeak,SNRFWHM);
0580
0581     calib->ScaleL       = SNRPeak;
0582     calib->ScaleWidthL  = SNRFWHM;
0583     out[4]    = SNRPeak;
0584     out[5]    = SNRFWHM;
0585   }
0586   return bmipLG;
0587 }
0588
0589
0590 bool TileSpectra::FitCorrCAEN(int verbosity){
0591   if (verbosity > 2) std::cout << "FitCorr cell ID: " << cellID << std::endl;
0592   TString funcName = Form("fcorr%sLGHGCellID%d",TileName.Data(),cellID);
0593
0594   Double_t fitRangeLG[2]  = {20., 250.};
0595   Double_t fitRangeHG[2]  = {350., 3500.};
0596   int fitStatus   = 0;
0597   int limitStatus = 0;
0598
0599
0600   LGHGcorr =  TF1(funcName.Data(),"pol1", fitRangeLG[0], fitRangeLG[1]);
0601   LGHGcorr.SetParameter(0,0.);
0602   LGHGcorr.SetParameter(1,10.);
0603   LGHGcorr.SetParLimits(1,0,100.);
0604
0605   fitStatus = hspectraLGHG.Fit(&LGHGcorr,"QRMNE0");
0606
0607   if (!(LGHGcorr.IsValid())){
0608     if (verbosity > 0) std::cout << "Skipped LGHG cell " << cellID << " fit failed" << std::endl;
0609     bcorrLGHG = false;
0610   } else {
0611     if ( LGHGcorr.GetParameter(1) == 0. || LGHGcorr.GetParameter(1) == 100. )
0612       limitStatus++;
0613     if (!(fitStatus == 4000 || fitStatus == 0)){ // only accept fits which succeeded in general
0614       if (verbosity > 0) std::cout << "Skipped LGHG cell " << cellID << " fit failed" << std::endl;
0615       bcorrLGHG = false;
0616     } else if (limitStatus > 0){                        // don't accept fits which reached the set limits
0617       if (verbosity > 0) std::cout << "Skipped LGHG cell " << cellID << " too many limits reached" << std::endl;
0618       bcorrLGHG = false;
0619     } else {
0620       bcorrLGHG= true;
0621     }
0622   }
0623   if (bcorrLGHG){
0624     calib->LGHGCorr     = LGHGcorr.GetParameter(1);
0625     calib->LGHGCorrOff  = LGHGcorr.GetParameter(0);
0626   }
0627   funcName = Form("fcorr%sHGLGCellID%d",TileName.Data(),cellID);
0628   HGLGcorr =  TF1(funcName.Data(),"pol1",fitRangeHG[0],fitRangeHG[1]);
0629   HGLGcorr.SetParameter(0,0.);
0630   HGLGcorr.SetParameter(1,0.1);
0631   HGLGcorr.SetParLimits(1,0.,1.);
0632
0633   fitStatus = hspectraHGLG.Fit(&HGLGcorr,"QRMNE0");
0634   limitStatus = 0;
0635
0636   if (!(HGLGcorr.IsValid())){
0637     if (verbosity > 0) std::cout << "Skipped HGLG cell " << cellID << " fit failed" << std::endl;
0638     bcorrHGLG = false;
0639   } else {
0640     if ( HGLGcorr.GetParameter(1) == 0. || HGLGcorr.GetParameter(1) == 1. )
0641       limitStatus++;
0642
0643     if (!(fitStatus == 4000 || fitStatus == 0)){ // only accept fits which succeeded in general
0644       if (verbosity > 0) std::cout << "Skipped HGLG cell " << cellID << " fit failed" << std::endl;
0645       bcorrHGLG = false;
0646     } else if (limitStatus > 0){                        // don't accept fits which reached the set limits
0647       if (verbosity > 0) std::cout << "Skipped HGLG cell " << cellID << " too many limits reached" << std::endl;
0648       bcorrHGLG = false;
0649     } else {
0650       bcorrHGLG= true;
0651     }
0652   }
0653   if (bcorrHGLG){
0654     calib->HGLGCorr     = HGLGcorr.GetParameter(1);
0655     calib->HGLGCorrOff  = HGLGcorr.GetParameter(0);
0656   }
0657   return true;
0658 }
0659
0660 bool TileSpectra::FitLGHGCorr(int verbosity, bool resetCalib){
0661   if (verbosity > 2) std::cout << "FitCorr cell ID: " << cellID << std::endl;
0662   TString funcName = Form("fcorr%sLGHGCellID%d",TileName.Data(),cellID);
0663
0664   Double_t fitRangeLG[2]  = {calib->PedestalSigL*5, 250.};
0665   int fitStatus   = 0;
0666   int limitStatus = 0;
0667
0668
0669   LGHGcorr =  TF1(funcName.Data(),"pol1", fitRangeLG[0], fitRangeLG[1]);
0670   LGHGcorr.SetParameter(0,0.);
0671   LGHGcorr.SetParameter(1,10.);
0672   LGHGcorr.SetParLimits(1,0,1000.);
0673
0674   fitStatus = hspectraLGHG.Fit(&LGHGcorr,"QRMNE0");
0675
0676   if (!(LGHGcorr.IsValid())){
0677     if (verbosity > 0) std::cout << "Skipped LGHG cell " << cellID << " fit failed" << std::endl;
0678     bcorrLGHG = false;
0679   } else {
0680     if ( LGHGcorr.GetParameter(1) == 0. || LGHGcorr.GetParameter(1) == 1000. )
0681       limitStatus++;
0682     if (!(fitStatus == 4000 || fitStatus == 0)){ // only accept fits which succeeded in general
0683       if (verbosity > 0) std::cout << "Skipped LGHG cell " << cellID << " fit failed" << std::endl;
0684       bcorrLGHG = false;
0685     } else if (limitStatus > 0){                        // don't accept fits which reached the set limits
0686       if (verbosity > 0) std::cout << "Skipped LGHG cell " << cellID << " too many limits reached" << std::endl;
0687       bcorrLGHG = false;
0688     } else {
0689       bcorrLGHG= true;
0690     }
0691   }
0692   if (bcorrLGHG && resetCalib){
0693     calib->LGHGCorr     = LGHGcorr.GetParameter(1);
0694     calib->LGHGCorrOff  = LGHGcorr.GetParameter(0);
0695   }
0696
0697   return true;
0698 }
0699
0700 bool TileSpectra::FitPedConstWave(int verbosity){
0701   if (verbosity > 2) std::cout << "FitCorr cell ID: " << cellID << std::endl;
0702   TString funcName = Form("fcorr%sPedWaveCellID%d",TileName.Data(),cellID);
0703
0704   Double_t fitRange[2]  = {0, hcorr.GetXaxis()->GetBinCenter(hcorr.GetNbinsX())};
0705   int fitStatus   = 0;
0706   int limitStatus = 0;
0707
0708   pedWave =  TF1(funcName.Data(),"pol0", fitRange[0], fitRange[1]);
0709   pedWave.SetParameter(0,80.);
0710   pedWave.SetParLimits(0,40.,120.);
0711
0712   hcorr.RebinY(4);
0713   fitStatus = hcorr.Fit(&pedWave,"QRMNE0");
0714
0715   if (!(pedWave.IsValid())){
0716     if (verbosity > 0) std::cout << "Skipped ped wave cell " << cellID << " fit failed" << std::endl;
0717     bpedWave = false;
0718   } else {
0719     if ( pedWave.GetParameter(0) == 40. || pedWave.GetParameter(0) == 120. )
0720       limitStatus++;
0721     if (!(fitStatus == 4000 || fitStatus == 0)){ // only accept fits which succeeded in general
0722       if (verbosity > 0) std::cout << "Skipped ped wave cell " << cellID << " fit failed" << std::endl;
0723       bpedWave = false;
0724     } else if (limitStatus > 0){                        // don't accept fits which reached the set limits
0725       if (verbosity > 0) std::cout << "Skipped ped wave cell " << cellID << " too many limits reached" << std::endl;
0726       bpedWave = false;
0727     } else {
0728       bpedWave= true;
0729     }
0730   }
0731
0732   if (bpedWave){
0733     calib->PedestalMeanL    = pedWave.GetParameter(0);
0734     calib->PedestalSigL     = pedWave.GetParError(0);
0735   }
0736
0737   return bpedWave;
0738 }
0739
0740
0741 bool TileSpectra::FitNoiseWithBG(double* out){
0742   return true;
0743 }
0744
0745 TH1D* TileSpectra::GetHG(){
0746   return &hspectraHG;
0747 }
0748
0749 TH1D* TileSpectra::GetLG(){
0750   return &hspectraLG;
0751 }
0752
0753 TH1D* TileSpectra::GetTriggPrim(){
0754   return &hTriggPrim;
0755 }
0756
0757 TH1D* TileSpectra::GetTOT(){
0758   return &hspectraTOT;
0759 }
0760
0761 TH1D* TileSpectra::GetTOA(){
0762   return &hspectraTOA;
0763 }
0764
0765 TH1D* TileSpectra::GetComb(){
0766   return &hcombined;
0767 }
0768
0769 TProfile* TileSpectra::GetLGHGcorr(){
0770   return &hspectraLGHG;
0771 }
0772 TProfile* TileSpectra::GetHGLGcorr(){
0773   return &hspectraHGLG;
0774 }
0775
0776 TProfile* TileSpectra::GetWave1D(){
0777   return &hWaveForm;
0778 }
0779
0780
0781 TProfile* TileSpectra::GetTOAADC(){
0782   return &hTOAADC;
0783 }
0784 TProfile* TileSpectra::GetADCTOT(){
0785   return &hADCTOT;
0786 }
0787
0788 TH2D* TileSpectra::GetCorr(){
0789   return &hcorr;
0790 }
0791
0792 TH2D* TileSpectra::GetCorrTOAADC(){
0793   return &hcorrTOAADC;
0794 }
0795
0796 TH2D* TileSpectra::GetCorrTOASample(){
0797   return &hcorrTOASample;
0798 }
0799
0800 TF1* TileSpectra::GetBackModel(int lh){
0801   if(lh==0 && bpedLG){
0802     return &BackgroundLG;
0803   }
0804   else if (lh==1 && bpedHG){
0805     return &BackgroundHG;
0806   }
0807   else return nullptr;
0808 }
0809
0810
0811
0812 TF1* TileSpectra::GetSignalModel(int lh){
0813   if(lh==0 && bmipLG){
0814     return &SignalLG;
0815   }
0816   else if (lh==1 && bmipHG){
0817     return &SignalHG;
0818   }
0819   else return nullptr;
0820 }
0821
0822 TileCalib* TileSpectra::GetCalib(){
0823   return calib;
0824 }
0825
0826 TF1* TileSpectra::GetCorrModel( int opt ){
0827   if(opt==0 && bcorrLGHG){
0828     return &LGHGcorr;
0829   }
0830   else if (opt==1 && bcorrHGLG){
0831     return &HGLGcorr;
0832   }
0833   else if (opt==2 && bpedWave){
0834     return &pedWave;
0835   }
0836   else return nullptr;
0837 }
0838
0839 void TileSpectra::Write( bool wFits = true){
0840   // write correct histos for CAEN output
0841   if (ROType != ReadOut::Type::Caen){
0842     hspectraHG.Write(hspectraHG.GetName(), kOverwrite);
0843     hspectraLG.Write(hspectraLG.GetName(), kOverwrite);
0844     hspectraLGHG.Write(hspectraLGHG.GetName(), kOverwrite);
0845     hspectraHGLG.Write(hspectraHGLG.GetName(), kOverwrite);
0846     if (bTriggPrim) hTriggPrim.Write(hTriggPrim.GetName(), kOverwrite);
0847     if ( wFits ){
0848       if(bpedHG)BackgroundHG.Write(BackgroundHG.GetName(), kOverwrite);
0849       if(bmipHG)SignalHG.Write(SignalHG.GetName(), kOverwrite);
0850       if(bpedLG)BackgroundLG.Write(BackgroundLG.GetName(), kOverwrite);
0851       if(bmipLG)SignalLG.Write(SignalLG.GetName(), kOverwrite);
0852       if(bcorrLGHG)LGHGcorr.Write(LGHGcorr.GetName(), kOverwrite);
0853       if(bcorrHGLG)HGLGcorr.Write(HGLGcorr.GetName(), kOverwrite);
0854     }
0855   // write correct histos for HGCROC output
0856   } else if (ROType != ReadOut::Type::Hgcroc) {
0857     hspectraHG.Write(hspectraHG.GetName(), kOverwrite);
0858     hspectraTOT.Write(hspectraTOT.GetName(), kOverwrite);
0859     hspectraTOA.Write(hspectraTOA.GetName(), kOverwrite);
0860     if (bTriggPrim) hTriggPrim.Write(hTriggPrim.GetName(), kOverwrite);
0861     hADCTOT.Write(hspectraTOT.GetName(), kOverwrite);
0862     hTOAADC.Write(hspectraTOA.GetName(), kOverwrite);
0863     if ( wFits ){
0864       if(bpedHG)BackgroundHG.Write(BackgroundHG.GetName(), kOverwrite);
0865       if(bmipHG)SignalHG.Write(SignalHG.GetName(), kOverwrite);
0866       if(bpedWave)pedWave.Write(pedWave.GetName(), kOverwrite);
0867       if(bwave)wave.Write(wave.GetName(), kOverwrite);
0868     }
0869   }
0870 }
0871
0872 void TileSpectra::WriteExt( bool wFits = true){
0873   // write correct histos for CAEN output
0874   if (ROType != ReadOut::Type::Caen){
0875     hspectraHG.Write(hspectraHG.GetName(), kOverwrite);
0876     hspectraLG.Write(hspectraLG.GetName(), kOverwrite);
0877     hspectraLGHG.Write(hspectraLGHG.GetName(), kOverwrite);
0878     if (extend == 1){
0879       hcombined.Write(hcombined.GetName(), kOverwrite);
0880       hspectraHGLG.Write(hspectraHGLG.GetName(), kOverwrite);
0881     } else if (extend == 2 || extend == 3 ){
0882       hcorr.Write(hcorr.GetName(), kOverwrite);
0883     }
0884     if ( wFits ){
0885       if(bpedHG)BackgroundHG.Write(BackgroundHG.GetName(), kOverwrite);
0886       if(bmipHG)SignalHG.Write(SignalHG.GetName(), kOverwrite);
0887       if(bpedLG)BackgroundLG.Write(BackgroundLG.GetName(), kOverwrite);
0888       if(bmipLG)SignalLG.Write(SignalLG.GetName(), kOverwrite);
0889       if(bcorrLGHG)LGHGcorr.Write(LGHGcorr.GetName(), kOverwrite);
0890     }
0891   // write correct histos for HGCROC output
0892   } else if (ROType != ReadOut::Type::Hgcroc) {
0893     hspectraHG.Write(hspectraHG.GetName(), kOverwrite);
0894     if (extend > 1 ) hcorr.Write(hcorr.GetName(), kOverwrite);
0895
0896     if (extend < 4){
0897       hspectraTOT.Write(hspectraTOT.GetName(), kOverwrite);
0898       hspectraTOA.Write(hspectraTOA.GetName(), kOverwrite);
0899     }
0900     if (extend == 2){
0901       hADCTOT.Write(hADCTOT.GetName(), kOverwrite);
0902     } else if (extend == 3){
0903       hWaveForm.Write(hWaveForm.GetName(), kOverwrite);
0904     } else if (extend == 4){
0905       hspectraLG.Write(hspectraLG.GetName(), kOverwrite);
0906     }
0907
0908     if ( wFits ){
0909       if(bpedHG)BackgroundHG.Write(BackgroundHG.GetName(), kOverwrite);
0910       if(bmipHG)SignalHG.Write(SignalHG.GetName(), kOverwrite);
0911       if(bpedWave)pedWave.Write(pedWave.GetName(), kOverwrite);
0912       if(bwave)wave.Write(wave.GetName(), kOverwrite);
0913     }
0914   }
0915 }
0916
0917 double TileSpectra::langaufun(double *x, double *par) {
0918
0919   //Fit parameters:
0920   //par[0]=Width (scale) parameter of Landau density
0921   //par[1]=Most Probable (MP, location) parameter of Landau density
0922   //par[2]=Total area (integral -inf to inf, normalization constant)
0923   //par[3]=Width (sigma) of convoluted Gaussian function
0924   //
0925   //In the Landau distribution (represented by the CERNLIB approximation),
0926   //the maximum is located at x=-0.22278298 with the location parameter=0.
0927   //This shift is corrected within this function, so that the actual
0928   //maximum is identical to the MP parameter.
0929
0930   // Numeric constants
0931   static double invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
0932   static double mpshift  = -0.22278298;       // Landau maximum location
0933
0934   // Control constants
0935   static double np = 100.0;      // number of convolution steps
0936   static double sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
0937
0938   // Variables
0939   double xx;
0940   double mpc;
0941   double fland;
0942   double sum = 0.0;
0943   double xlow,xupp;
0944   double step;
0945   double i;
0946
0947
0948   // MP shift correction
0949   mpc = par[1] - mpshift * par[0];
0950
0951   // Range of convolution integral
0952   xlow = x[0] - sc * par[3];
0953   xupp = x[0] + sc * par[3];
0954
0955   step = (xupp-xlow) / np;
0956
0957   // Convolution integral of Landau and Gaussian by sum
0958   for(i=1.0; i<=np/2; i++) {
0959     xx = xlow + (i-.5) * step;
0960     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
0961     sum += fland * TMath::Gaus(x[0],xx,par[3]);
0962
0963     xx = xupp - (i-.5) * step;
0964     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
0965     sum += fland * TMath::Gaus(x[0],xx,par[3]);
0966   }
0967
0968   return (par[2] * step * sum * invsq2pi / par[3]);
0969 }
0970
0971
0972 int TileSpectra::langaupro(double *params, double &maxx, double &FWHM) {
0973
0974   // Searches for the location (x value) at the maximum of the
0975   // Landau-Gaussian convolute and its full width at half-maximum.
0976   //
0977   // The search is probably not very efficient, but it's a first try.
0978   double p,x,fy,fxr,fxl;
0979   double step;
0980   double l,lold;
0981   int i = 0;
0982   int MAXCALLS = 10000;
0983
0984   // Search for maximum
0985   p = params[1] - 0.1 * params[0];
0986   step = 0.05 * params[0];
0987   lold = -2.0;
0988   l    = -1.0;
0989
0990   while ( (l != lold) && (i < MAXCALLS) ) {
0991     i++;
0992     lold = l;
0993     x = p + step;
0994     l = langaufun(&x,params);
0995     if (l < lold)
0996       step = -step/10;
0997     p += step;
0998   }
0999
1000   if (i == MAXCALLS)
1001     return (-1);
1002   maxx = x;
1003   fy = l/2;
1004
1005   // Search for right x location of fy
1006   p = maxx + params[0];
1007   step = params[0];
1008   lold = -2.0;
1009   l    = -1e300;
1010   i    = 0;
1011
1012   while ( (l != lold) && (i < MAXCALLS) ) {
1013     i++;
1014     lold = l;
1015     x = p + step;
1016     l = TMath::Abs(langaufun(&x,params) - fy);
1017     if (l > lold)
1018       step = -step/10;
1019     p += step;
1020   }
1021
1022   if (i == MAXCALLS)
1023     return (-2);
1024   fxr = x;
1025
1026   // Search for left x location of fy
1027   p = maxx - 0.5 * params[0];
1028   step = -params[0];
1029   lold = -2.0;
1030   l    = -1e300;
1031   i    = 0;
1032
1033   while ( (l != lold) && (i < MAXCALLS) ) {
1034     i++;
1035     lold = l;
1036     x = p + step;
1037     l = TMath::Abs(langaufun(&x,params) - fy);
1038     if (l > lold)
1039       step = -step/10;
1040     p += step;
1041   }
1042
1043   if (i == MAXCALLS)
1044     return (-3);
1045
1046
1047   fxl = x;
1048   FWHM = fxr - fxl;
1049   return (0);
1050 }
1051
1052 //__________________________________________________________________________________________________________
1053 // find bin with largest content in given range and return X
1054 //__________________________________________________________________________________________________________
1055 double TileSpectra::GetMaxXInRangeHG(double minX = -10000, double maxX = -10000) {
1056   Double_t largestContent     = 0;
1057   Int_t minBin = 1;
1058   Int_t maxBin = hspectraHG.GetNbinsX()+1;
1059   if (minX != -10000) minBin = hspectraHG.GetXaxis()->FindBin(minX);
1060   if (maxX != -10000) maxBin = hspectraHG.GetXaxis()->FindBin(maxX)+0.0001;
1061   Int_t largestBin = minBin;
1062   for (Int_t i= minBin; i < maxBin; i++){
1063     if (largestContent < hspectraHG.GetBinContent(i)){
1064       largestContent = hspectraHG.GetBinContent(i);
1065       largestBin = i;
1066     }
1067   }
1068   return hspectraHG.GetBinCenter(largestBin);
1069 }
1070
1071
1072
1073 //__________________________________________________________________________________________________________
1074 // find bin with largest content in given range and return X
1075 //__________________________________________________________________________________________________________
1076 double TileSpectra::GetMaxXInRangeLG(double minX = -10000, double maxX = -10000) {
1077   Double_t largestContent     = 0;
1078   Int_t minBin = 1;
1079   Int_t maxBin = hspectraLG.GetNbinsX()+1;
1080   if (minX != -10000) minBin = hspectraLG.GetXaxis()->FindBin(minX);
1081   if (maxX != -10000) maxBin = hspectraLG.GetXaxis()->FindBin(maxX)+0.0001;
1082   Int_t largestBin = minBin;
1083   for (Int_t i= minBin; i < maxBin; i++){
1084     if (largestContent < hspectraLG.GetBinContent(i)){
1085       largestContent = hspectraLG.GetBinContent(i);
1086       largestBin = i;
1087     }
1088   }
1089   return hspectraLG.GetBinCenter(largestBin);
1090 }
1091
1092 //__________________________________________________________________________________________________________
1093 // find bad channels
1094 //__________________________________________________________________________________________________________
1095 short TileSpectra::DetermineBadChannel(){
1096   short bc = -64;
1097
1098   FitFixedNoise();
1099
1100   // double sigL   = calib->PedestalSigL;
1101   // double sigH   = calib->PedestalSigH;
1102
1103   double sigLN   = BackgroundLG.GetParameter(2);
1104   double sigHN   = BackgroundHG.GetParameter(2);
1105   // std::cout << "HG: " << sigH << "\t" << sigHN <<  "\t LG: " << sigL << "\t" << sigLN  << std::endl;
1106
1107 //   double bgH3    = hspectraHG.Integral(hspectraHG.FindBin(-sigH*3-0.01),hspectraHG.FindBin(sigH*3+0.01));
1108 //   double bgL3    = hspectraLG.Integral(hspectraLG.FindBin(-sigL*3-0.01),hspectraLG.FindBin(sigL*3+0.01));
1109 //   double intH3   = hspectraHG.Integral(hspectraHG.FindBin(sigH*3+0.01),hspectraHG.GetNbinsX()-1);
1110 //   double intL3   = hspectraLG.Integral(hspectraLG.FindBin(sigL*3+0.01),hspectraLG.GetNbinsX()-1);
1111 //   double intH5   = hspectraHG.Integral(hspectraHG.FindBin(sigH*5+0.01),hspectraHG.GetNbinsX()-1);
1112 //   double intL5   = hspectraLG.Integral(hspectraLG.FindBin(sigL*5+0.01),hspectraLG.GetNbinsX()-1);
1113 //
1114 //   if (intH5 > 100 && intL5 > 100){
1115 //     bc = 3;
1116 //   } else if (intH5 > 100 && intL3 > 100){
1117 //     bc = 2;
1118 //     std::cout << bc << "  \t HG: " << intH5  << "\t" << intH5/bgH3 << "\t" << " LG: "<< intL3  << "\t" << intL3/bgL3 << std::endl;
1119 //   } else if (intH3 > 100 && intL3 > 100){
1120 //     bc = 1;
1121 //     std::cout << bc << "  \t HG: " << intH3  << "\t" << intH3/bgH3 << "\t" << " LG: "<< intL3  << "\t"  << intL3/bgL3 << std::endl;
1122 //   } else {
1123 //     bc = 0;
1124 //     std::cout << bc << "  \t HG: " << intH3  << "\t" << intH3/bgH3 << "\t" << " LG: "<< intL3  << "\t"  << intL3/bgL3 << std::endl;
1125 //   }
1126   BackgroundHG.SetRange(-100,400);
1127   // double nH3  = BackgroundHG.Integral(-3*sigHN,5*sigHN);
1128   // double nH5  = BackgroundHG.Integral(-5*sigHN,7*sigHN);
1129   double nH3  = hspectraHG.Integral(hspectraHG.FindBin(-3*sigHN),hspectraHG.FindBin(3*sigHN));
1130   double nH5  = hspectraHG.Integral(hspectraHG.FindBin(-5*sigHN),hspectraHG.FindBin(5*sigHN));
1131   double fH   = hspectraHG.Integral(1,hspectraHG.GetNbinsX()-1);
1132   BackgroundLG.SetRange(-100,400);
1133   double nL3  = BackgroundLG.Integral(-3*sigLN,3*sigLN);
1134   double nL5  = BackgroundLG.Integral(-5*sigLN,5*sigLN);
1135   double fL   = hspectraLG.Integral(1,hspectraLG.GetNbinsX()-1);
1136
1137   double rH3  = (fH > 0) ? (fH-nH3)/fH : 0;
1138   double rH5  = (fH > 0) ? (fH-nH5)/fH : 0;
1139   double rL3  = (fL > 0) ? (fL-nL3)/fL : 0;
1140   double rL5  = (fL > 0) ? (fL-nL5)/fL : 0;
1141
1142   if (rH5 > 1e-6 && rL5 > 1e-6){
1143     bc = 3;
1144   } else if (rH5 > 1e-6 && rL3 > 0.005){
1145     bc = 2;
1146   } else if (rH3 > 0.005 && rL3 > 0.005){
1147     bc = 1;
1148   } else {
1149     bc = 0;
1150   }
1151   calib->BadChannel = bc;
1152
1153   return bc;
1154 }
1155
1156 void TileSpectra::SetBadChannelInCalib(short s ){
1157   calib->BadChannel = s;
1158 }