Hello Experts.
I’m using TF1 Convolution to fit an array of histograms of Analog-to-Digital Converter Signal with a landau convoluted by a gaussian. Currently, not all histograms are being fitted with the convolution function. For the one’s that don’t get fit, I get the error: Error in TF1::Moment: Integral zero over range. I’m trying to fix it by fixing the number of points that is used to plot the fits with SetNofPointsFFT(). It’s currently at 75,000, but it doesn’t accommodate for each histograms even though each histogram is suppose to contain the same information (e.g. the same particle). I wish I can provide the root file that I’m using for this macro, but it’s larger than 3 MB. I’m currently producing another one that might be smaller, but we’ll see. Hope to see a response soon. If you need me to provide any more information, please let me know. Bye for now.
mip_histogram_fitter.C (17.1 KB)
-Charlie Clark
#include <iostream>
#include <fstream>
#include <vector>
#include <sstream>
using namespace std;
#include "TFile.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TPDF.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLatex.h"
#include "TString.h"
#include "TLine.h"
#include "TF1Convolution.h"
#include "TROOT.h"
//void drawTower(TH1D*, Int_t, Int_t, Float_t, Float_t);
//Bool_t isBadTower(Int_t);
stringstream ss;
void mip_histogram_fitter(const Char_t* mipRootfile = "./0506/mipcombined22_0426.root", const Char_t* psName = "./0506/mips_conv.pdf", const Char_t* rootFilename = "./0506/towerMipFits_conv_f10_1.root", const Char_t* mipOutputName = "./0506/mips_conv_f10_1.gains")
{
gROOT->Macro("loadMuDst.C");
//StEmcGeom *mEmcGeom = StEmcGeom::instance("bemc");
const Int_t nTowers = 4800;
Int_t towerStatus[nTowers];
Double_t fitMean[nTowers], fitMeanError[nTowers];
Double_t fitParams[5];
cout << "Input Filename: " << mipRootfile << endl;
cout << "Plot Filename: " << psName << endl;
cout << "Gain Filename: " << mipOutputName << endl;
gStyle->SetCanvasColor(10);
gStyle->SetCanvasBorderMode(0);
gStyle->SetStatColor(10);
gStyle->SetOptTitle(0);
gStyle->SetOptDate(0);
gStyle->SetOptFit(111);
gStyle->SetOptStat("e");
gStyle->SetPalette(1);
/* Declare pointers for histograms before output file. Instantiate after,
this allows us to use outfile->Write() to write everything at once */
TH1D *towerHisto[nTowers];
//TH1F *htowerTrackDelta[nTowers];
//TH2F *etaPhiMean, *etaPhiStatus, *etaPhiRatio;
TFile* outfile = new TFile(rootFilename,"RECREATE");
Float_t pi = TMath::Pi();
// Instantiate histograms
//etaPhiMean = new TH2F("etaPhiMean","Fit Means", 40, -1., 1., 120, -pi, pi);
//etaPhiStatus = new TH2F("etaPhiStatus","Status Codes", 40, -1., 1., 120, -pi, pi);
//etaPhiRatio = new TH2F("etaPhiMeanRatio","Mean Ratios", 40, -1., 1., 120, -pi, pi);
for (Int_t i = 0; i < nTowers; ++i)
towerHisto[i] = new TH1D(Form("towerHisto_%i",i+1),Form("Tower %i",i+1),250,-50.5,199.5);
TFile *rootfile = new TFile(mipRootfile);
for(Int_t y = 0; y < nTowers; ++y)
towerHisto[y]->Add((TH1D*)rootfile->Get(Form("tower_histo_%i",y+1)));
rootfile->Close();
if(rootfile) delete rootfile;
/****************** Begin MIP Fitting for each Tower ******************/
TF1 *gaussian[nTowers];
TF1 *landau[nTowers];
TF1Convolution *towerFit_Initial[nTowers];
TF1 *towerFit[nTowers];
TString fitName;
TString gausName;
TString landauName;
//TH1F *hRMS = new TH1F("hRMS","RMS",)
for(Int_t iTow = 0; iTow < nTowers; ++iTow){
Int_t softId = iTow+1;
//Float_t mTowerEta = 0., mTowerPhi = 0.;
//mEmcGeom->getEta(softId,mTowerEta);
//mEmcGeom->getPhi(softId,mTowerPhi);
//if(softId >= 3381 && softId <= 3540) continue;
//
//if (mTowerEta >= 0.1 && mTowerEta <= 0.15 && mTowerPhi <= -2.0 && mTowerPhi >= -2.2) cout << "TowerID: " << softId << endl;
towerStatus[iTow] = 0;
fitMean[iTow] = 0.;
fitMeanError[iTow] = 0.;
Double_t fitLow = towerHisto[iTow]->GetMean() - towerHisto[iTow]->GetRMS();
Double_t fitHigh = towerHisto[iTow]->GetMean() + towerHisto[iTow]->GetRMS()*2.;
towerHisto[iTow]->GetXaxis()->SetRangeUser(6.,100.);
cout << "Now fitting softId: " << softId << endl;
ss.str("");
ss << softId;
fitName = (TString)("towerFit_" + ss.str());
gausName = (TString)("gausFit_" + ss.str());
landauName = (TString)("landauFit_" + ss.str());
//Double_t fitLow = 0;
//Double_t fitHigh = 100;
landau[iTow] = new TF1(landauName,"landau",fitLow,fitHigh);
landau[iTow]->SetParameters(towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()),towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin()),towerHisto[iTow]->GetRMS());
//towerHisto[iTow]->Fit(landau[iTow],"RQ");
//Double_t mpv = landau[iTow]->GetParameter(1);
//Double_t rms_landau = landau[iTow]->GetParameter(2);
//Double_t tail_range_low = 10;
//Double_t tail_range_high = 40;
gaussian[iTow] = new TF1(gausName,"gaus",fitLow,fitHigh);
gaussian[iTow]->SetParameters(towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()),towerHisto[iTow]->GetMean(),towerHisto[iTow]->GetRMS());
//towerHisto[iTow]->Fit(gaussian[iTow],"RQ");
Double_t area_gaus_tails = gaussian[iTow]->Integral(fitLow,fitHigh);
//cout << "area_gaus_tails " << softId << " " << area_gaus_tails << endl;
Double_t area_tower_histograms = towerHisto[iTow]->Integral();
//cout << "area_tower_histograms " << softId << " " << area_gaus_tails << endl;
Double_t gaussian_scale_parameter = area_tower_histograms / area_gaus_tails;
cout << "Tower: " << softId << " " << "Gaussian Integral: " << area_gaus_tails << " " << "Histogram Integral: " << area_tower_histograms << " " << "Scale Parameter: " << gaussian_scale_parameter << endl;
towerFit_Initial[iTow] = new TF1Convolution(landau[iTow],gaussian[iTow],fitLow,fitHigh,true);
//towerFit_Initial->SetRange(fitLow,fitHigh);
towerFit_Initial[iTow]->SetNofPointsFFT(75000);
towerFit[iTow] = new TF1(fitName, *towerFit_Initial ,fitLow,fitHigh,5);
//towerFit[iTow]->SetNpx(90000);
cout << "Name of Fit " << towerFit[iTow]->GetName() << endl;
//towerFit[iTow] = new TF1(fitName,"[0]*Gaus(x,[1],[2])*Landau(x,[3],[4])",fitLow,fitHigh);
// Fit Parameters
towerFit[iTow]->SetParameter(0,towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin())); //Amplitude
towerFit[iTow]->SetParameter(1,towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin())); //Landau MPV
towerFit[iTow]->SetParameter(2,towerHisto[iTow]->GetRMS()); //Landau RMS
towerFit[iTow]->SetParameter(3,gaussian_scale_parameter); //Gaussian Scale Parameter
towerFit[iTow]->SetParameter(4,towerHisto[iTow]->GetRMS()); //Gaussian RMS
//towerFit[iTow]->SetParLimits(2,0,1000);
towerFit[iTow]->SetLineColor(kBlue);
towerFit[iTow]->SetLineWidth(1.5);
// If entries less than 25 we can't fit the tower, so we skip it
if(towerHisto[iTow]->Integral(56,156) < 25){
towerStatus[iTow] = 13;
//int isGood = 0;
//if (towerStatus[iTow] == 1) isGood = 1;
//else isGood = 0;
//etaPhiStatus->Fill(mTowerEta,mTowerPhi,towerStatus[iTow]);
//etaPhiStatus->Fill(mTowerEta,mTowerPhi,isGood);
continue;
}
towerHisto[iTow]->Fit(towerFit[iTow],"RQ");
towerStatus[iTow]+=1;
fitMean[iTow] = towerFit[iTow]->Mean(fitLow,fitHigh);
// Get the fit parameters in an array
Double_t fitParams[5] = {0.,0.,0.,0.,0.};
fitParams[0] = towerFit[iTow]->GetParameter(0);
fitParams[1] = towerFit[iTow]->GetParameter(1);
fitParams[2] = towerFit[iTow]->GetParameter(2);
fitParams[3] = towerFit[iTow]->GetParameter(3);
fitParams[4] = towerFit[iTow]->GetParameter(4);
//Calculate the variance of the fit
Double_t fitVariance = towerFit[iTow]->Variance(0., 100., fitParams);
Double_t fitSigma = sqrt(fitVariance);
Double_t entriesInFit = towerHisto[iTow]->Integral(51,151); // Bin 51 centered over zero
cout << "entriesInFit " << softId << " " << entriesInFit << endl;
//Set the status codes of the towers
// if(fitMean[iTow] < 6) towerStatus[iTow] += 10; //Bad mean
/*
if((softId)%20==0){
if (abs(fitMean[iTow]-towerHisto[iTow]->GetMean()) > 15)
towerStatus[iTow]+=8;
}
else{
if(abs(fitMean[iTow]-towerHisto[iTow]->GetMean()) > 10)
towerStatus[iTow]+=8;
}
*/
//Check for hard coded bad towers
//if(isBadTower(softId)) towerStatus[iTow] = 7;
//if(softId == 4) towerStatus[iTow] = 1;
//if(softId >= 3381 && softId <= 3540) towerStatus[iTow] = 9;
//Calculate error on mean for good statuses
if (towerStatus[iTow] == 1)fitMeanError[iTow] = fitSigma/sqrt(entriesInFit);
// Fill eta/phi for good statuses
/*
if (towerStatus[iTow] == 1){
etaPhiMean->Fill(mTowerEta,mTowerPhi,fitMean[iTow]);
etaPhiRatio->Fill(mTowerEta,mTowerPhi,fitMean[iTow]/towerHisto[iTow]->GetMean());
}*/
//etaPhiStatus->Fill(mTowerEta,mTowerPhi,towerStatus[iTow]);
}// Towers loop
/******************** End MIP Fitting for each Tower ********************/
/******************* Begin Drawing the Fitted Spectra *******************/
TPDF *ps = new TPDF(psName);
TCanvas *canvas = new TCanvas("can","can1",100,100,600.,800.);
Int_t pad;
cout << endl << "Begin Tower Drawing Loop" << endl;
for (Int_t iTow = 0; iTow < nTowers; ++iTow){
gStyle->SetOptFit(1);
Int_t softId = iTow + 1;
//Float_t mTowerEta = 0., mTowerPhi = 0.;
//mEmcGeom->getEta(softId,mTowerEta);
//mEmcGeom->getPhi(softId,mTowerPhi);
//if(softId >= 3381 && softId <= 3540) continue;
cout << "Drawing tower " << softId << endl;
if(iTow%20 == 0){
canvas->Update();
ps->NewPage();
canvas->Clear();
canvas->Divide(4,5);
pad = 1;
}
canvas->cd(pad);
if(towerStatus[iTow] > 1)
towerHisto[iTow]->SetLineColor(kRed);
// Draw the tower onto the pad
//drawTower(towerHisto[iTow],softId,towerStatus[iTow],mTowerEta,mTowerPhi);
towerHisto[iTow]->GetXaxis()->SetTitle("ADC");
towerHisto[iTow]->GetXaxis()->CenterTitle();
towerHisto[iTow]->Draw();
Double_t lineMax = towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()) + 30;
if(towerStatus[iTow] == 1){
TLine *avgValueLine = new TLine(fitMean[iTow],0.,fitMean[iTow],lineMax);
avgValueLine->SetLineColor(kRed);
avgValueLine->SetLineWidth(1.25);
avgValueLine->Draw("same");
}
pad++;
}//End tower draw loop
/*
canvas->Update();
ps->NewPage();
canvas->Clear();
canvas->Divide(1,2);
canvas->cd(1);
etaPhiMean->GetXaxis()->SetTitle("#eta");
etaPhiMean->GetXaxis()->CenterTitle();
etaPhiMean->GetYaxis()->SetTitle("#phi");
etaPhiMean->GetYaxis()->CenterTitle();
etaPhiMean->GetZaxis()->SetRangeUser(0,45);
etaPhiMean->Draw("colz");
canvas->cd(2);
//etaPhiStatus->GetZaxis()->SetRangeUser(0,20);
cout << endl << etaPhiStatus->GetNbinsX();
for (int ix = 1; ix <= etaPhiStatus->GetNbinsX(); ix++) {
for (int iy = 1; iy <= etaPhiStatus->GetNbinsY(); iy++) {
if (etaPhiStatus->GetBinContent(ix, iy) != 1) etaPhiStatus->SetBinContent(ix, iy, 0);
}
}
etaPhiStatus->GetZaxis()->SetRangeUser(0,2);
etaPhiStatus->GetXaxis()->SetTitle("#eta");
etaPhiStatus->GetXaxis()->CenterTitle();
etaPhiStatus->GetYaxis()->SetTitle("#phi");
etaPhiStatus->GetYaxis()->CenterTitle();
etaPhiStatus->Draw("colz");
canvas->Update();
ps->NewPage();
canvas->Clear();
canvas->Divide(1,2);
canvas->cd(1);
etaPhiRatio->GetXaxis()->SetTitle("#eta");
etaPhiRatio->GetXaxis()->CenterTitle();
etaPhiRatio->GetYaxis()->SetTitle("#phi");
etaPhiRatio->GetYaxis()->CenterTitle();
etaPhiRatio->Draw("colz");
canvas->cd(2);
TH2F *ratioClone = (TH2F*)etaPhiRatio->Clone();
ratioClone->GetZaxis()->SetRangeUser(0.7,1.);
ratioClone->GetXaxis()->SetTitle("#eta");
ratioClone->GetXaxis()->CenterTitle();
ratioClone->GetYaxis()->SetTitle("#phi");
ratioClone->GetYaxis()->CenterTitle();
ratioClone->Draw("colz");
ps->Close();
*/
cout << "Tower Drawing Loop Complete" << endl;
/******************** End Drawing the Fitted Spectra ********************/
/******************** Begin Output of MIP Fit Data ********************/
ofstream mipOutput(mipOutputName);
cout << "Now writing to Gains file: " << mipOutputName << endl;
Int_t nGood = 0;
Int_t nZero = 0;
Int_t nBad = 0;
/*
for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerFit[iTow]->GetName() << endl;
cout << " " << endl;
for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetMean() << endl;
cout << " " << endl;
for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetRMS() << endl;
cout << " " << endl;
for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin()) << endl;
cout << " " << endl;
for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin()) << endl;
cout << " " << endl;
for(Int_t iTow = 0; iTow < nTowers; ++iTow) cout << fitMean[iTow]/towerHisto[iTow]->GetMean() << endl;
*/
for(Int_t iTow = 0; iTow < nTowers; ++iTow){
//cout << fitMeanError[iTow] << endl;
//cout << "Inside of Tower loop to write to Gains file" << endl;
Int_t softId = iTow + 1;
Double_t mipAdc = fitMean[iTow];
Double_t mipError = fitMeanError[iTow];
// Double_t mipRMS = towerFit[iTow]->GetParameter(2);
Double_t histMean = towerHisto[iTow]->GetMean();
Double_t histRMS = towerHisto[iTow]->GetRMS();
Double_t histIP1 = towerHisto[iTow]->GetBinContent(towerHisto[iTow]->GetMaximumBin());
Double_t histIP3 = towerHisto[iTow]->GetBinCenter(towerHisto[iTow]->GetMaximumBin());
Double_t Ratio = fitMean[iTow]/towerHisto[iTow]->GetMean();
//cout << "Writing to Gains file" << endl;
mipOutput << softId << " " << mipAdc << " " << mipError << " " << towerStatus[iTow] << endl;
if(towerStatus[iTow]==1)nGood++;
if(towerStatus[iTow]==13)nZero++;
if(towerStatus[iTow]==7)nBad++;
}
cout << "Finished writing to Gains file: " << mipOutputName << endl;
mipOutput.close();
/********************* End Output of MIP Fit Data *********************/
outfile->Write();
outfile->Close();
cout << "Finished tower fits" << endl << endl;
cout << "Number of Good Statuses = " << nGood << endl;
cout << "Number of Bad Statuses = " << nBad << endl;
cout << "Number of Empty Towers = " << nZero << endl;
// Delete geometry pointer
//if (mEmcGeom) delete mEmcGeom;
}
// Load StEmcGeom for tower eta & phi
// End of main routine
/*
void drawTower(TH1D* histo, Int_t id, Int_t status, Float_t towerEta, Float_t towerPhi)
{
// Set histogram options
histo->GetXaxis()->SetTitle("ADC");
histo->GetXaxis()->CenterTitle();
histo->Draw(); // Used Sumw2() in histogram maker to calculate errors
// Write softId and status code on plot
Char_t towerTitle[100];
Char_t etaTitle[100];
Char_t phiTitle[100];
sprintf(etaTitle,"#eta:%1.2f",towerEta);
sprintf(phiTitle,"#phi:%1.2f",towerPhi);
TLatex etaLatex;
TLatex phiLatex;
etaLatex.SetTextSize(0.1);
phiLatex.SetTextSize(0.1);
sprintf(towerTitle,"%i",id);
TLatex titleLatex;
titleLatex.SetTextSize(0.15);
//if(status!=1) titleLatex.SetTextColor(kRed); // Non - colorblind friendly color
if(status!=1) titleLatex.SetTextColor(kOrange);
phiLatex.DrawTextNDC(0.55,0.2,phiTitle);
etaLatex.DrawTextNDC(0.55,0.33,etaTitle);
titleLatex.DrawTextNDC(0.55,0.45,towerTitle);
if(status!=1){
Char_t towerCode[100];
sprintf(towerCode,"%i",status);
TLatex statusLatex;
statusLatex.SetTextSize(0.15);
//statusLatex.SetTextColor(kRed); // Non - colorblind friendly color
statusLatex.SetTextColor(kOrange);
statusLatex.DrawTextNDC(0.47,0.78,towerCode);
}
}
Bool_t isBadTower(Int_t id) {
//Empty towers = {31,113,134,166,257,426,509,541,793,817,832,844,846,956,989,1005,1012,1027,1057,1081,1122,1158,1159,1180,1189,1190,1202,1203,1204,1207,1214,1222,1224,1242
// 1243,1244,1329,1341,1354,1369,1388,1434,1438,1443,1448,1533,1534,1575,1593,1594,1595,1597,1654,1668,1674,1679,1705,1807,1866,2019,2104,2111,2168,2214,2281
// 2282,2283,2386,2392}
//Hot towers = {365,375,533,939,945,1211,1223,1826,2032,2234,2300}
Int_t badTowers[] = {34,106,139,160,245,266,267,286,287,365,368,375,504,514,533,561,562,615,616,633,637,638,649,650,653,657,671,673,674,740,760,
813,814,837,857,875,899,901,939,945,953,954,993,1026,1044,1045,1046,1048,1080,1100,1198,1199,1200,1207,1211,1219,1220,1223,
1237,1280,1303,1320,1337,1341,1348,1353,1354,1377,1380,1407,1447,1448,1530,1574,1575,1597,1612,1620,1654,1753,
1765,1766,1826,1856,1877,1878,2032,2040,2073,2077,2092,2093,2097,2219,2234,2300,2305,2529,2569,
2589,2590,2613,2700,2834,2890,2969,3005,3061,3070,3071,3381,3384,3387,3391,3392,3393,3394,3395,
3396,3397,3398,3403,3407,3408,3411,3413,3416,3419,3426,3430,3432,3436,3437,3442,3445,3447,3448,
3451,3452,3455,3456,3458,3460,3464,3468,3474,3476,3478,3481,3483,3484,3485,3486,3488,3489,3496,
3497,3498,3499,3500,3501,3502,3503,3504,3507,3509,3510,3515,3516,3519,3523,3526,3528,3531,3532,3534,3535,3538,3539,3678,3679,3718,3738,3794,4017,4018,4019,4057,4177,4223,
4339,4560,4657,4677,4678,4684};
for (int i = 0; i < sizeof(badTowers); i++)
{
if (badTowers[i] == id) return true;
}
return false;
}*/