OmegaQuadrupoleFourEnergiesCombination.C

Go to the documentation of this file.

TString CompleteFilename(TString mu,
                TString Q,
                TString Psf
                        )
{
        TString FileEnd = "mu"+mu+",Q"+Q+"_OQP_Psf"+Psf+".root";
        return FileEnd;

}

Double_t GetWidthkeV(TF1 *func, Double_t Whichwidth=0.5)
{
        Double_t Max = func->GetMaximum();
        //cout << "max"<<Max<<endl;
        //cout << "Xmax"<<func->GetMaximumX()<<endl;

        Double_t XMin = func->GetX(Max*Whichwidth,func->GetMaximumX()-0.20,func->GetMaximumX());
        Double_t XMax = func->GetX(Max*Whichwidth,func->GetMaximumX(),func->GetMaximumX()+0.20);

        return (XMax-XMin)*1000;
}




Double_t QuadGausPlusStep (Double_t *x,Double_t *par)
{

                        Double_t Gausian1 = par[0] * TMath::Exp(- 0.5*pow((x[0] - par[1])/par[2],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
                        Double_t Gausian2 = par[3] * TMath::Exp(- 0.5*pow((x[0] - par[4])/par[5],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
                        Double_t Gausian3 = par[6] * TMath::Exp(- 0.5*pow((x[0] - par[7])/par[8],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
                        Double_t Gausian4 = par[9] * TMath::Exp(- 0.5*pow((x[0] - par[10])/par[11],2)); // par[0]: amplitude, par[1]: x[0]0; par[2]: sigma
                        Double_t SmoothedStep = par[12] * TMath::Erfc( (x[0] - par[1]) /(TMath::Sqrt(2)*par[2])); // par[3]: amplitude

        return Gausian1+Gausian2+Gausian3+Gausian4+SmoothedStep ;
}


int OmegaQuadrupoleFourEnergiesCombination(
                        TString mu = "-2.5",
                        TString Q = "-2.8",
                        TString Psf = "1",
                        TString Filename1 ="/data/work/kpha1/steinen/Gamma/Ana/CombinedData/Combined_Ana_Geo43_E0.5189MeV_Evts2860000_FileEvts11440_Gen1_ST12_mu-2.5,Q-2.8_OQP_Psf1.root",
                        TString Filename2 ="/data/work/kpha1/steinen/Gamma/Ana/CombinedData/Combined_Ana_Geo43_E0.5198MeV_Evts2580000_FileEvts10320_Gen1_ST12_mu-2.5,Q-2.8_OQP_Psf1.root",
                        TString Filename3 ="/data/work/kpha1/steinen/Gamma/Ana/CombinedData/Combined_Ana_Geo43_E0.5208MeV_Evts2330000_FileEvts9320_Gen1_ST12_mu-2.5,Q-2.8_OQP_Psf1.root",
                        TString Filename4 ="/data/work/kpha1/steinen/Gamma/Ana/CombinedData/Combined_Ana_Geo43_E0.5215MeV_Evts2110000_FileEvts8440_Gen1_ST12_mu-2.5,Q-2.8_OQP_Psf1.root"
        )
{
//      TString Filename1 ="Combined_Ana_Geo43_E0.5*MeV_Evts2860000_FileEvts11440_Gen1_ST12_"+CompleteFilename(mu,Q,Psf);
//      TString Filename2 ="Combined_Ana_Geo43_E0.5*MeV_Evts2580000_FileEvts10320_Gen1_ST12_"+CompleteFilename(mu,Q,Psf);
//      TString Filename3 ="Combined_Ana_Geo43_E0.5*MeV_Evts2330000_FileEvts9320_Gen1_ST12_"+CompleteFilename(mu,Q,Psf);
//      TString Filename4 ="Combined_Ana_Geo43_E0.5*MeV_Evts2110000_FileEvts8440_Gen1_ST12_"+CompleteFilename(mu,Q,Psf);


        andi::setCustomStyle();

        TString Path = "/data/work/kpha1/steinen/Gamma/Ana/CombinedData/";
        TString Filenames[4];

                Filenames[0]=Filename1;
                Filenames[1]=Filename2;
                Filenames[2]=Filename3;
                Filenames[3]=Filename4;

                cout<< Filenames[0].Data()<<endl;
        TFile *InputFiles[4];
        TFile *OutputFile;
        TH1D *hAddUp;// = new TH1D("hGammaSpectrumCombined","Combined Gamma spectrum; Energy [MeV];Counts [a.u.]",7000,1,0.7)
        TH1D *hTemporary[4];
        TF1 *fGausFit[4];
        //TF1 *fGausFitCombined = new TF1("fGausFitCombined","gaus(0)+gaus(3)+gaus(6)+gaus(9)+pol2(12)",0.515,0.525);
        TF1 *fGausFitCombined = new TF1("fGausFitCombined",QuadGausPlusStep,0.515,0.525,13);
        TString Fitname;
        TCanvas *canCombo= new TCanvas("cCombo","cCombo",800,600);
        TCanvas *canAdd= new TCanvas("cAdd","cAdd",800,600);
        TCanvas *can[4];
        can[0]=new TCanvas("c0","c0",800,600);
        can[1]=new TCanvas("c1","c1",800,600);
        can[2]=new TCanvas("c2","c2",800,600);
        can[3]=new TCanvas("c3","c3",800,600);
        for (int i = 0; i<4; i++)
        {
                InputFiles[i]= new TFile(Filenames[i].Data());
                InputFiles[i]->GetObject("GamEnergy;1",hTemporary[i]);
                hTemporary[i]->SetDirectory(0);
                InputFiles[i]->Close();
                canCombo->cd();

                hTemporary[i]->Draw();
                can[i]->cd();
                Fitname.Form("fGausFit%d",i);
                fGausFit[i]= new TF1(Fitname.Data(),"gaus",0.515,0.525);
                hTemporary[i]->Fit(Fitname.Data(),"r");
                for (int j=0;j<3;j++)
                        fGausFitCombined->SetParameter(3*i+j,fGausFit[i]->GetParameter(j));
                fGausFitCombined->FixParameter(1+3*i,fGausFitCombined->GetParameter(1+3*i));
                hTemporary[i]->Draw("");
        }
        hAddUp = new TH1D("hAddUp","4 lines added up;Energy [MeV];Counts [a.u.]",6000,0,0.6);
        for(int i = 0; i<4;i++)
        {
                hAddUp->Add(hTemporary[i],1);
        }
        canAdd->cd();
        fGausFitCombined->SetLineColor(kRed);
        hAddUp->Fit("fGausFitCombined","r");
        //hAddUp->Draw();

        //fGausFitCombined->Draw("same");

        for(int i = 0; i<4;i++)
                {
                fGausFit[i]->Draw("same");
                cout <<"gaus " << i << endl;
                //cout << fGausFit[i]->GetParameter(0) << "\t" << fGausFitCombined->GetParameter(0+i*3) << endl;
                //cout << fGausFit[i]->GetParameter(1) << "\t" << fGausFitCombined->GetParameter(1+i*3) << endl;
                //cout << fGausFit[i]->GetParameter(2) << "\t" << fGausFitCombined->GetParameter(2+i*3) << endl;
                }

        Double_t Mean =fGausFitCombined->Mean(0.516,0.524);
        Double_t StandardDev =sqrt(fGausFitCombined->Variance(0.516,0.524));
        //Double_t Dip = fGausFitCombined->GetMinimumX(fGausFit[1]->GetParameter(1),fGausFit[2]->GetParameter(1));                      // the dip does not work for mu = -2.0, there is simply no dip!!!
        Double_t Skewness = fGausFitCombined->CentralMoment(3,0.516,0.524);
        Double_t Kurtosis = fGausFitCombined->CentralMoment(4,0.516,0.524);
        Double_t MaxX = fGausFitCombined->GetMaximumX(0.516,0.524);
        cout <<GetWidthkeV(fGausFitCombined,0.5)<<endl;
        cout <<GetWidthkeV(fGausFitCombined,0.3)<<endl;
        cout <<GetWidthkeV(fGausFitCombined,0.1)<<endl;
        cout <<"Variance"<<StandardDev<<endl;
        cout <<"Mean"<<Mean<<endl;
        //cout << "Dip"<< Dip<<endl;
        cout <<"Skewness"<<Skewness<<endl;
        cout <<"Kurtosis"<<Kurtosis<<endl;

        //cout << "mu" << "\t"<<"Q"<<"\t"<<"Psf"<< "\t"<< "Mean" << "\t" << "StandardDev" << "\t" << "Dip" << endl;
        //cout << mu.Data() << "\t"<<Q.Data()<<"\t"<<Psf.Data()<< "\t"<< Mean << "\t" << StandardDev << "\t" << Dip << endl;

        TString Outfile = Path + "OmegaResults/OmegaResults.txt";
        ifstream outfile (Outfile.Data());                      // check if file exist, if not, header line will be written later on
        bool FileExists =outfile.good();
        outfile.close();
        ofstream outfile2 (Outfile.Data(),std::ofstream::app|std::ofstream::out);
        if (!FileExists)
                        outfile2 << "mu" << "\t"<<"Q"<<"\t"<<"Psf"<< "\t"<< "MeanPos" << "\t" << "StandardDev" << "\t" << "Skewness"<< "\t" << "Kurtosis" << "\t" << "MaximumX" << endl;
        outfile2<< mu.Data() << "\t"<<Q.Data()<<"\t"<<Psf.Data()<< "\t"<< Mean << "\t" << StandardDev << "\t" << Skewness <<"\t" << Kurtosis<< "\t" << MaxX  << endl;
        outfile2.close();

        TString RootOutName = Path + "/OmegaResults/mu"+mu+"_Q"+Q +"__Psf"+Psf+".root"  ;
        TFile *RootOutfile = new TFile(RootOutName.Data(),"RECREATE");
                hAddUp->Write();
                fGausFitCombined->Write();
                canAdd->Write();
                for(int i = 0;i<4; i++)
                {
                        fGausFit[i]->Write();
                        hTemporary[i]->Write();
                        can[i]->Write();
                }
                RootOutfile->Close();
  return 0;
}