FlowMC/Flow/calcRes.c

#include <TMath.h>
#include "../Utils/Funcs.cc"

void calcRes(TString inName, TString outName)
{
    TFile *fi = new TFile(inName.Data(),"read");

    // Input histograms
    TProfile *p_res2RXNEvsRXNW[3], *p_res3RXNEvsRXNW[3];
    TProfile *p_res2TPCEvsTPCW[4], *p_res3TPCEvsTPCW[4];
    TProfile *p_res2RXNEvsTPCM[3], *p_res3RXNEvsTPCM[3];
    TProfile *p_res2RXNWvsTPCM[3], *p_res3RXNWvsTPCM[3];

    for (int iGap=0; iGap<4; iGap++)
    {
        p_res2TPCEvsTPCW[iGap] = (TProfile*) fi->Get(Form("p_res2TPCEvsTPCW_gap%i",iGap));
        p_res3TPCEvsTPCW[iGap] = (TProfile*) fi->Get(Form("p_res3TPCEvsTPCW_gap%i",iGap));
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        p_res2RXNEvsRXNW[irxn] = (TProfile*) fi->Get(Form("p_res2RXNEvsRXNW_%i",irxn));
        p_res3RXNEvsRXNW[irxn] = (TProfile*) fi->Get(Form("p_res3RXNEvsRXNW_%i",irxn));

        // p_res2RXNEvsTPCM[irxn] = (TProfile*) fi->Get(Form("p_res2RXNEvsTPCM_%i",irxn));
        // p_res3RXNEvsTPCM[irxn] = (TProfile*) fi->Get(Form("p_res3RXNEvsTPCM_%i",irxn));
        // p_res2RXNWvsTPCM[irxn] = (TProfile*) fi->Get(Form("p_res2RXNWvsTPCM_%i",irxn));
        // p_res3RXNWvsTPCM[irxn] = (TProfile*) fi->Get(Form("p_res3RXNWvsTPCM_%i",irxn));
    }

    // Output histograms
    TH1F *hres2TPC[4], *hres3TPC[4];
    TH1F *hres2RXN[3], *hres3RXN[3];
    TH1F *hres2RXNfull[3], *hres3RXNfull[3];
    for (int iGap=0; iGap<4; iGap++)
    {
        hres2TPC[iGap] = new TH1F(Form("hres2TPC_gap%i",iGap),Form("hres2TPC_gap%i",iGap),10,0,100);
        hres3TPC[iGap] = new TH1F(Form("hres3TPC_gap%i",iGap),Form("hres3TPC_gap%i",iGap),10,0,100);
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        hres2RXN[irxn] = new TH1F(Form("hres2RXN_part%i",irxn),Form("hres2RXN_part%i",irxn),10,0,100);
        hres3RXN[irxn] = new TH1F(Form("hres3RXN_part%i",irxn),Form("hres3RXN_part%i",irxn),10,0,100);

        hres2RXNfull[irxn] = new TH1F(Form("hres2RXNfull_part%i",irxn),Form("hres2RXNfull_part%i",irxn),10,0,100);
        hres3RXNfull[irxn] = new TH1F(Form("hres3RXNfull_part%i",irxn),Form("hres3RXNfull_part%i",irxn),10,0,100);
    }

    // 2-sub TPC method
    Double_t res, res2;
    for (int iGap=0; iGap<4; iGap++)
    {
        for (int ibin = 0; ibin<hres2TPC[iGap]->GetNbinsX(); ibin++)
        {
            res2 = p_res2TPCEvsTPCW[iGap]->GetBinContent(ibin+1);
            if (res2 > 0) res = TMath::Sqrt(res2);
            hres2TPC[iGap]->SetBinContent(ibin+1,res);
            hres2TPC[iGap]->SetBinError(ibin+1,p_res2TPCEvsTPCW[iGap]->GetBinError(ibin+1));
        }

        for (int ibin = 0; ibin<hres3TPC[iGap]->GetNbinsX(); ibin++)
        {
            res2 = p_res3TPCEvsTPCW[iGap]->GetBinContent(ibin+1);
            if (res2 > 0) res = TMath::Sqrt(res2);
            hres3TPC[iGap]->SetBinContent(ibin+1,res);
            hres3TPC[iGap]->SetBinError(ibin+1,p_res3TPCEvsTPCW[iGap]->GetBinError(ibin+1));
        }
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        for (int ibin = 0; ibin<hres2RXN[irxn]->GetNbinsX(); ibin++)
        {
            res2 = p_res2RXNEvsRXNW[irxn]->GetBinContent(ibin+1);
            if (res2 > 0) res = TMath::Sqrt(res2);
            hres2RXN[irxn]->SetBinContent(ibin+1,res);
            hres2RXN[irxn]->SetBinError(ibin+1,p_res2RXNEvsRXNW[irxn]->GetBinError(ibin+1));
        }

        for (int ibin = 0; ibin<hres3RXN[irxn]->GetNbinsX(); ibin++)
        {
            res2 = p_res3RXNEvsRXNW[irxn]->GetBinContent(ibin+1);
            if (res2 > 0) res = TMath::Sqrt(res2);
            hres3RXN[irxn]->SetBinContent(ibin+1,res);
            hres3RXN[irxn]->SetBinError(ibin+1,p_res3RXNEvsRXNW[irxn]->GetBinError(ibin+1));
        }
    }
    Double_t chi;
    for (int irxn=0; irxn<3; irxn++)
    {
        for (int ibin = 0; ibin<hres2RXNfull[irxn]->GetNbinsX(); ibin++)
        {
            res2 = p_res2RXNEvsRXNW[irxn]->GetBinContent(ibin+1);
            if (res2 > 0)
            {
                res = TMath::Sqrt(res2);
                chi = GetChi(res,1,50);
                chi *= TMath::Sqrt(2.);
                res = GetRes(chi,1);
                hres2RXNfull[irxn]->SetBinContent(ibin+1,res);
                hres2RXNfull[irxn]->SetBinError(ibin+1,p_res2RXNEvsRXNW[irxn]->GetBinError(ibin+1));
            }
        }

        for (int ibin = 0; ibin<hres3RXNfull[irxn]->GetNbinsX(); ibin++)
        {
            res2 = p_res3RXNEvsRXNW[irxn]->GetBinContent(ibin+1);
            if (res2 > 0)
            {
                res = TMath::Sqrt(res2);
                chi = GetChi(res,1,50);
                chi *= TMath::Sqrt(2.);
                res = GetRes(chi,1);
                hres3RXNfull[irxn]->SetBinContent(ibin+1,res);
                hres3RXNfull[irxn]->SetBinError(ibin+1,TMath::Sqrt(2)*p_res3RXNEvsRXNW[irxn]->GetBinError(ibin+1));
            }
        }
    }

    std::cout << " ResV2_Full <cosLR>" << std::endl;
    for (int ibin=0; ibin<hres2RXNfull[0]->GetNbinsX();ibin++)
    {
        std::cout << hres2RXNfull[0]->GetBinContent(ibin+1) << " " << p_res2RXNEvsRXNW[0]->GetBinContent(ibin+1) << std::endl;
    }


    TFile *fo = new TFile(outName.Data(),"recreate");
    fo->cd();
    for (int iGap=0; iGap<4; iGap++)
    {
        hres2TPC[iGap]->Write();
        hres3TPC[iGap]->Write();
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        hres2RXN[irxn]->Write();
        hres3RXN[irxn]->Write();
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        hres2RXNfull[irxn]->Write();
        hres3RXNfull[irxn]->Write();
    }
    fo->Close();
}