Hello experts,
I am trying to fit over 100k TH1D with convoluted Landau distributions but my macro will stop working and exit root automatically after fitting about 380 histograms. I am running ROOT 6.08/00 on my own laptop (Memory 8 GB 1600 MHz DDR3). While my laptop is probably not a very powerful machine, I am wondering if there is any problem with my code (attached below).
I tried to be very careful about the memory management, for example, I reset the smart pointers when they are done even before it goes out of the scope. I also deleted all the raw pointers when they are done. Therefore I am not sure what was causing the memory issue. It will be nice if you could take a look at the following macro. Sorry for the long code, I tried to only keep the main part already.
void FitnMIPhisto(TH1D* fitThis);
#define nMipsMax 4 // what is the maximum number of MIPs you want to consider?
TF1* MipPeak[nMipsMax];
Double_t myfunc(Double_t* x, Double_t* param); // Fit Function used by Minuit
void AnalyzeThreeDhistoLite(TString infilename="Histos.root", TString mOutputFileName = "Fit_test"){
//Histos.root contains a lot of TH3D. This macro takes the TH3D, project it to different bins and get a TH1D, then the TH1D will be fitted by the FitnMIPhisto(TH1D* fitThis) method.
double pi=TMath::Pi();
//Create a title page for the PDF
unique_ptr<TCanvas> TitlePage(new TCanvas("Title","Title",1400,1000));
TPaveText tptitle(0.1,0.1,0.8,0.8);
tptitle.Draw();
TString OutputPDFName = "Analysis_"+mOutputFileName+".pdf";
TitlePage->SaveAs(Form("%s%s",OutputPDFName.Data(),"("));//cannot .reset() for now because it will be saved again at the end
double centpercent[]={80,70,60,50,40,30,20,10,5,0};
int vzrange[]={-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40};
for(int cent=0;cent<1;cent++){
for(int ew=0;ew<1;ew++){
for(int ivz=0;ivz<1;ivz++){
for(int ieta=0; ieta<16; ieta++){
//----------------Create a sub title page for each eta bin---------------------
TPaveText tptsec(0.1,0.3,0.9,0.7);
unique_ptr<TCanvas> SectionPage(new TCanvas("ringsec","ringsec",1400,1000));
//SectionPage->cd();
tptsec.AddText(Form("ring %d",ieta+1));
tptsec.Draw();
SectionPage->SaveAs(OutputPDFName);
SectionPage.reset();
//--------------Open the root file and get the histograms needed------------
unique_ptr<TFile> tf(new TFile(infilename.Data(),"READ"));
unique_ptr<TH1F> mVz((TH1F*)tf->Get(Form("VzCent%d",cent)));
unique_ptr<TH3F> ThreeDhist((TH3F*)tf->Get(Form("dNdphidnMIPdetaCent%dEW%dVz%d",cent,ew,ivz)));
mVz->SetDirectory(0);
ThreeDhist->SetDirectory(0);
tf.reset();//release the memory
double Nevents=mVz->GetBinContent(ivz+1);
mVz.reset();
ThreeDhist->Scale(1.0/Nevents);
ThreeDhist->Scale(1.0/(ThreeDhist->GetXaxis()->GetBinWidth(1)*ThreeDhist->GetYaxis()->GetBinWidth(1)));
double etaValue = ThreeDhist->GetZaxis()->GetBinCenter(ieta+1);
double etaBinWidth = ThreeDhist->GetZaxis()->GetBinWidth(ieta+1);
ThreeDhist->GetZaxis()->SetRange(ieta+1,ieta+1);
unique_ptr<TH2F> TwoDhist((TH2F*)ThreeDhist->Project3D("yxe"));
TwoDhist->SetDirectory(0);//TwoDhist belongs to gDirectory which is tf.
ThreeDhist.reset();// Now finished with ThreeDhist, free the memory
unique_ptr<TCanvas> c1(new TCanvas(Form("dNdnMIPfitting"),Form("dNdnMIPfitting"),1400,1000));
c1->Divide(4,3);//for 4*3=12 nMIP dis. fitting
//--------------loop over phi bins-------------
for(int iphi=0;iphi<12;iphi++){
double phiValue = TwoDhist->GetXaxis()->GetBinCenter(iphi+1);
c1->cd(iphi+1)->SetLogy();//interesting, has to cd() here
TH1D* dNdnMIPprojection=TwoDhist->ProjectionY(Form("Cent%dEW%dVz%dEta%dPhi%d",cent,ew,ivz,ieta,iphi),iphi+1,iphi+1,"e");//have to give it a unique name
FitnMIPhisto(dNdnMIPprojection);
dNdnMIPprojection->SetTitle(Form("Cent=%d EW=%d ivz=%d eta=%f phi=%f",cent,ew,ivz, etaValue,phiValue));
dNdnMIPprojection->DrawCopy();
dNdnMIPprojection->DrawCopy("same"); // yes, again, so it gets on top of the stats box
delete dNdnMIPprojection;
for(int i=0;i<nMipsMax;i++){
delete MipPeak[i];
}
}//end looping over phi bins, round one [0,11]
c1->SaveAs(OutputPDFName);
c1.reset();
}//end looping over eta bins
}//end looping over vz
}//end looping over ew
}//end looping over cent
cout<<"finish fitting"<<endl;
TitlePage->SaveAs(Form("%s%s",OutputPDFName.Data(),"]"));
}
// ---------------------------------- Multi-MIP fit ------------------------------------------
void FitnMIPhisto(TH1D* fitThis){
double SingleMipPeakStartingValue = 1.0; // where is the 1-mip peak?
double FitRangeLow = 0.75; // low edge of range along the x-axis
double FitRangeHigh = 8.0; // high edge of range along the x-axis
// (1) ====================== Set Up the functions ===============================
double xlo = fitThis->GetXaxis()->GetXmin();
double xhi = 4.0*fitThis->GetXaxis()->GetXmax();
TF1 func("MultiMipFit",myfunc,xlo,xhi,nMipsMax+2);
func.SetLineWidth(1);
// ---------------- here are the individual functions corresponding to 1,2,3,4... MIPs. func uses them...
MipPeak[0] = new TF1("1MIP","TMath::Landau(x,[0],[1],1)",xlo,xhi);
for (Int_t nMIP=2; nMIP<=nMipsMax; nMIP++){
TF1Convolution* c = new TF1Convolution(MipPeak[nMIP-2],MipPeak[0],xlo,xhi,true);
MipPeak[nMIP-1] = new TF1(Form("%dMIPs",nMIP),c,xlo,xhi,2*nMIP);
}
// (2) ======================= Set up the fit ======================================
for (Int_t nmip=0; nmip<nMipsMax; nmip++){
func.SetParName(nmip,Form("%dMIPweight",nmip+1));
}
func.SetParName(nMipsMax,"MPV");
func.SetParName(nMipsMax+1,"WIDbyMPV");
func.SetParameter(nMipsMax+1,0.15);
func.SetParameter(nMipsMax,SingleMipPeakStartingValue);
// (3) ======================= Do the fit ============================================
fitThis->Fit("MultiMipFit","","",FitRangeLow,FitRangeHigh);
}
// ------------------------------- here is the fitting function -----------------------------------
Double_t myfunc(Double_t* x, Double_t* param){
// parameters 0...(nMipsMax-1) are the weights of the N-MIP peaks
// and the last two parameters, index nMipsMax and nMipsMax+1, are single-MIP MPV and WID/MPV, respectively
Double_t ADC = x[0];
Double_t SingleMipMPV = param[nMipsMax];
Double_t WID = SingleMipMPV*param[nMipsMax+1];
Double_t fitval=0.0;
for (Int_t nMip=0; nMip<nMipsMax; nMip++){
Double_t Weight = abs(param[nMip]);
for (Int_t imip=0; imip<2*(nMip+1); imip+=2){
MipPeak[nMip]->SetParameter(imip,SingleMipMPV);
MipPeak[nMip]->SetParameter(imip+1,WID);
}
fitval += Weight*(*MipPeak[nMip])(ADC);
}
return fitval;
}
// -------------------------------------------------------------------------------------------------
I was running it without compilation.
Thanks,
Xiaoyu
_ROOT Version: ROOT 6.08/00
Platform: Not Provided
Compiler: Not Provided