#ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include #include "TFile.h" #include "TH1.h" #include "TH1D.h" #include "TString.h" #include "TGraph.h" #include "TRandom3.h" #include "TCanvas.h" #include "TDirectory.h" #include "RooRealVar.h" #include "RooDataSet.h" #include "RooDataHist.h" #include "RooGaussian.h" #include "RooPlot.h" using namespace RooFit ; using namespace std; struct Vavilov_Func { Vavilov_Func() {} double operator() (const double *x, const double *p) { double kappa = p[0]; double beta2 = p[1]; return p[4]*( pdf.Pdf( (x[0]-p[2])/p[3], kappa,beta2) ); } ROOT::Math::VavilovAccurate pdf; }; void Test(){ const Int_t nbfiles = 10;//100; //240 //number of files const Int_t nLayers = 16; //number of layers const Int_t nMicrons= 300; //number of slices const Int_t nParticles=50000; // nb of particles const Int_t Emaxdep= 50; Double_t edep; Double_t eventID; Double_t Detector; Double_t step; TH1::AddDirectory(kFALSE); vector ion= {"alpha", "ion_6_12_6","ion_14_28_14","ion_26_56_26"}; //{"alpha", "ion_6_12_6","ion_14_28_14","ion_26_56_26"} //{"proton","alpha", "He3", "ion_3_7_3", "ion_4_9_4", "ion_5_11_5", "ion_6_12_6", "ion_7_14_7", "ion_8_16_8", "ion_9_19_9", "ion_10_20_10", "ion_11_23_11", "ion_12_24_12" , "ion_13_27_13", "ion_14_28_14", "ion_26_56_26"}; Int_t nspecies = ion.size(); // number of particles to plot vector ener= {1, 2, 3, 5, 8, 10, 12, 15, 18, 20, 25, 30, 40, 50, 60, 80, 100}; //{1, 2, 3, 5, 8, 10, 12, 15, 18, 20, 25, 30, 40, 50, 60, 80, 100} Int_t nener = ener.size(); // number of particles to plot vector angle= {0, 1, 2, 5, 8, 10, 15, 20, 30, 40, 50, 70}; Int_t nangle = angle.size(); // number of particles to plot Int_t nI=0 ;Int_t nE=16; Int_t nA=0; //4 11 16 TString current_read_file_name = TString("processed/config300_")+ ion[nI] + TString("_")+ Form("%d", ener[nE])+ TString("MeV_nucl_") + Form("%d", angle[nA]) +TString("_deg.root"); cout<< "File "<< current_read_file_name << " opened" << endl; TFile *myfile= new TFile(current_read_file_name,"READ"); gDirectory->pwd(); myfile->cd("Microns"); TString title = TString("Edep Distribution ")+ TString(" Alpha 3MeV/nucl, ")+ Form("%d", nParticles); TIter next(gDirectory->GetListOfKeys()); TKey *key; Int_t count=0; TH1F *h1, *h5, *hmax, *htot; while ((key = (TKey*)next())) { TClass *cl = gROOT->GetClass(key->GetClassName()); if (!cl->InheritsFrom("TH1")) continue; TH1F *myh1 = (TH1F*)key->ReadObj(); if (count==0) { h1= (TH1F*)myh1->Clone(); h1->SetDirectory(0); htot= (TH1F*)myh1->Clone(); htot->SetDirectory(0); }else if (count==4) { htot->Add((TH1D*)myh1); h5= (TH1F*)htot->Clone(); h5->SetDirectory(0); }else if (count==61) { htot->Add((TH1D*)myh1); hmax= (TH1F*)myh1->Clone(); hmax->SetDirectory(0); }else { htot->Add(myh1); } ++count; } myfile->Close(); delete myfile; cout<<"Files closed"<SetTitle(ion[nI]+ TString(" ")+Form("%d", ener[nE])+ TString("MeV/nucl: Distribution of Edep in first micron")); h5->SetTitle(ion[nI]+ TString(" ")+Form("%d", ener[nE])+ TString("MeV/nucl: Distribution of Edep in 5 first microns")); htot->SetTitle(ion[nI]+ TString(" ")+Form("%d", ener[nE])+ TString("MeV/nucl: Distribution of Edep in 300um")); TH1F *h1bis=(TH1F*)h1->Clone(); TH1F *h1bisbis=(TH1F*)h1->Clone(); /////%%%%%% Roofit %%%%%/// Int_t nbins = h1->GetSize()-2; Double_t mean = h1->GetMean(); Double_t sigma = h1->GetStdDev() ; Double_t minX= h1->GetXaxis()->GetBinLowEdge(1); Double_t maxX = h1->GetXaxis()->GetBinUpEdge(nbins); RooRealVar x("E","E", minX, Emaxdep,"keV"); RooDataHist dh("dh","dh",x,Import(*h1)); // Gauss /*RooPlot* frame = x.frame(Title("Imported TH1 with Gaussian")); dh.plotOn(frame); RooRealVar mg("mean","mean",mean, mean*0.1,mean*5) ; RooRealVar sg("sigma","sigma",sigma,0.1*sigma,5*sigma) ; RooGaussian gauss("gauss","gauss",x,mg,sg) ; gauss.plotOn(frame); gauss.paramOn(frame,Format(("NE"),AutoPrecision(1)),Layout(.55,.9,.9)) ; TCanvas * c1 = new TCanvas("c1","c1", 1200, 800); frame->Draw() ; Double_t chi1 = frame->chiSquare(2); TPaveLabel *t1 = new TPaveLabel(0.55,0.6,0.9,0.72, Form("#chi^{2} = %f", chi1),"brNDC"); t1->Draw(); frame->addObject(t1) ; frame->Draw();*/ // Vavilov X Gauss RooRealVar mg2("mean","mean",mean,0.1*mean,5*mean) ; RooRealVar sg2("sigma","sigma",sigma,0.1*sigma,5*sigma) ; RooGaussian gauss2("gauss","gauss",x,mg2,sg2) ; RooRealVar kappa("kappa","mean landau",0.1, 0.01,10); RooRealVar beta2("beta2","sigma landau",0.01,0,0.5); RooAbsPdf *vavilov= bindPdf("vavilov",ROOT::Math::vavilov_accurate_pdf,x,kappa,beta2); //RooAbsPdf *vavilov= RooFunctorPdfBinding("vavilov",Vavilov_Func,x,kappa,beta2,mv,sg,amp); vavilov->Print() ; RooFFTConvPdf vxg("vxg", "Vavilov (X) gauss", x, *vavilov, gauss2); vxg.fitTo(dh); RooPlot* frame2 = x.frame(Title("Imported TH1 with Vavilov X Gaussian")); dh.plotOn(frame2); vxg.plotOn(frame2); vxg.paramOn(frame2,Format(("NE"),AutoPrecision(1)),Layout(.45,.9,.9)) ; TCanvas * c2 = new TCanvas("c2","c2", 1200, 800); frame2->Draw(); Double_t chi2 = frame2->chiSquare(3); TPaveLabel *t2 = new TPaveLabel(0.55,0.55,0.9,0.66, Form("#chi^{2} = %f", chi2),"brNDC"); t2->Draw(); frame2->addObject(t2) ; frame2->Draw(); /* // Landau X Gauss RooRealVar mg3("mean","mean",0) ; RooRealVar sg3("sigma","sigma",sigma,0.1*sigma,5*sigma) ; RooGaussian gauss3("gauss","gauss",x,mg3,sg3) ; RooRealVar ml("mu","mean landau",mean, mean-sigma, mean+sigma) ; RooRealVar sl("c","sigma landau",1,0.001,Emaxdep) ; RooLandau landau("lx","lx",x,ml,sl) ; RooFFTConvPdf lxg("lxg", "landau (X) gauss", x, landau, gauss3); lxg.fitTo(dh); RooPlot* frame3 = x.frame(Title("Imported TH1 with Landau X Gaussian")); dh.plotOn(frame3); lxg.plotOn(frame3); lxg.paramOn(frame3,Format(("NE"),AutoPrecision(1)),Layout(.55,.9,.9)) ; TCanvas * c3 = new TCanvas("c2","c2", 1200, 800); frame3->Draw(); Double_t chi3 = frame3->chiSquare(3); TPaveLabel *t3 = new TPaveLabel(0.55,0.6,0.9,0.72, Form("#chi^{2} = %f", chi3),"brNDC"); t3->Draw(); frame3->addObject(t3) ; frame3->Draw(); */ }