attenuation_calibration_comparisons.C

Go to the documentation of this file.

//Author: Prabhjot Singh (prabhjot@fnal.gov). 
//Dated April 19 2016 03:00 am
//This program is a template macro to plot and compare
//a. uncorrected PE/cm profiles,
//b. threshold and shielding corrections, 
//c. attenution fitting plots for any randomly chosen plane and cell,
//d. calibration constants, coeff_exp, attenlength, background and chisq, 
//e. list and number of calibrated channels for the following cases
//1. ND MC, 
//2. ND Data, 
//3. FD MC, 
//4. FD Data,
//5. Tricell vs. Trajectory (two cases in one time)
//6. Any two cases in one time e.g. 
//              a. different periods in same analysis, 
//              b. new calibration vs. old calibration, 
//              c. one analyis vs. another
//7. Many vs. one case, e.g. all 2nd analysis cases vs. 1st analysis e.g. period1, period2, epoch3b, epoch3c vs. 1st analysis.
//
//8. You have to set the booleans given at the start of the program to select one or more cases.
//

#include "HistogramAttr.h"
#include "Header.h"

bool uncorrectedPEplots         = false;     // set true if you want plots of uncorrected PE/cm vs. W
bool thresholdshadowingplots    = false;    // set true if you want plots of threshold and shadowing correction factors.
bool correctedPEplots           = true;    // set true if you want plots of corrected PE/cm vs. W. That is, attenuation fitting plots and thier ratios.
bool calibrationconstantsplots  = false;    // set it to true if you wans to make plots of calibration constants

bool nd                         = true;     // set true for the ND comparisons
bool fd                         = false;    // set true for the FD comparisons
bool mc                         = false;    // set true for the Monte Carlo based comparisons
bool data                       = true;     // set true for the Data based comparisons
bool tricellvstrajectory        = false;    // set true for the tricell vs. trajectory comparisons
bool onetoone                   = true;     // set true for one vs. one case e.g. point 6 above
bool manytoone                  = false;    // set true for many vs. one case e..g. point 7 above
bool nddatamuoncatcher          = true;     // set true for the ND data muon catcher plots

attenuation_calibration_comparisons(){
    gStyle->SetOptStat(" ");
    //-----------------------------------------------------------------------------------------------------
    //set of all .root files for uncorrected PE/cm vs. w plots
    char* c1            = "/nova/ana/users/prabhjot/gridjob/Work_with_Adam_and_Brian/nd/NDPCHitData_S150224/mainbody_and_muoncatcher/mainboby/ConsolidateViewsMode_false/Fromgrid_aftercosmiccalib/histofiles/calibNDPCHitData_hist_mainbody_ConsolidateViewsMode_false_S150224.root"; // base file
    char* c2            = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/calib_hist_2ndAna_fulldetector_nddata_mainbody_period2ddactivity.root"; // 2nd file

    //.root file for uncorrected PE/cm for tricell vs trajectory study
    char* c3            = ""; 

    //set of .root files for threshold and shadowing corrections
    char* cthreshold1   = "";   // base file
    char* cthreshold2   = "";   // 2nd file

    //set of .root files to plot attenuation fitting plots
    char* cattenuation1 = "/nova/app/users/prabhjot/nova_test_S150224/AttenuationCalibration/nd/pchitdata/fromgrid/mainbody_and_muoncatcher/mainboby/ConsolidateViewsMode_false/Atten_Fit_Results.root";  // base file
    char* cattenuation2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/FitAttenuation/Atten_Fit_Results.root";   // 2nd file

    //set of .root files to plot attenuation fitting plots for the ND data muon catcher
    char* cattenuation_muoncatcher1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root";   // base file
    char* cattenuation_muoncatcher2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3b_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root";   // 2nd file

    //set of attenuation constants .txt files for attenuation constants comparisons
    char* cconstants1   = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/comparison_plots/calib_atten_consts_uniq_sort_1st_ana.txt";   // base file
    char* cconstants2   = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/comparison_plots/calib_atten_consts_uniq_sort_period2.txt";   // 2nd file
    //---------------------------------------------------------------------------------------------------
    //---------------------------------------------------------------------------------------------------

    TString hist1;              // base histogram from the base file. You want to compare other histograms wrt base histogram.
    TString hist2;              // histogram from the second file that you want to compare with the base histogram.
    TString det;                // name of the detector.
    TString mcdata;             // labeling for the Monte Carlo or the data .
    TString analysis1;          // labeling for the base file.
    TString analysis2;          // labelling for the second file.
    TString var;                // variable labelling on the x axis. It is W (cm).
    TString plane;              // plane number.
    TString cell;               // cell number.

    analysis1       = "2nd_analysis_period2";
    analysis2       = "2nd_analysis_epoch3b";

    hist1           = "h_WPE_corr_xy";  
    hist2           = "h_WPE_corr_xy";  //_corr_xy_ for the tricell, _corr_traj_ for the trajectory

    var             = "W";
    plane           = "003";
    cell            = "023";

    //ND MC two cases comparison
    if(nd && onetoone){
        det             = "nd";
        if(mc) mcdata   = "mc";
        if(data) mcdata = "data";

        if(tricellvstrajectory && uncorrectedPEplots)   uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
        else if(uncorrectedPEplots)                     uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);

        if(mc && thresholdshadowingplots)               thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);

        if(data && nddatamuoncatcher && correctedPEplots)corrected_PE(cattenuation_muoncatcher1, cattenuation_muoncatcher2, det, mcdata, analysis1, analysis2, plane, cell);
        else if(correctedPEplots)                       corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);

        if(calibrationconstantsplots)                   calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2); 
    }//end of if(nd && mc && onetoone)

    //FD two case comparison
    if(fd && onetoone){
        det             = "fd";
        if(mc) mcdata   = "mc";
        if(data) mcdata = "data";

        if(tricellvstrajectory && uncorrectedPEplots)   uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
        else if(uncorrectedPEplots)                     uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);

        if(mc && thresholdshadowingplots)               thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);

        if(correctedPEplots)                            corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);

        if(calibrationconstantsplots)                   calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
    }//end of if(fd && mc && onetoone)
}//attenuation_calibration_comparisons()
/////////////////////////////////////////////////////////////////////////////////

//----------------------------------------------------------------------------------------------------------------------
uncorrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
    TFile* f1   = new TFile(c1, "read");
    TFile* f2   = new TFile(c2, "read");
    TString viewstr;
    TString viewStr;
    for(int view = 0; view < 2; ++view){
        if(nddatamuoncatcher) viewstr = view ? "muy" : "mux";
        else if(!nddatamuoncatcher) viewstr = view ? "y" : "x";
        if(nddatamuoncatcher) viewStr = view ? "muY" : "muX";
        else if(!nddatamuoncatcher) viewStr = view ? "Y" : "X";

        f1->cd();
        if(f1->GetListOfKeys()->Contains("calib")) 
            TH2F* h2D_1         = (TH2F*)f1->Get("calib/"+hist1+"_in"+viewStr);
        if(f1->GetListOfKeys()->Contains("make"))
            TH2F* h2D_1         = (TH2F*)f1->Get("make/"+hist1+"_in"+viewStr);
        TProfile* ProfileX_1    = h2D_1->ProfileX();
        if(det=="nd" && mcdata=="data") ProfileX_1->Rebin();

        f2->cd();
        if(f2->GetListOfKeys()->Contains("calib"))
            TH2F* h2D_2         = (TH2F*)f2->Get("calib/"+hist2+"_in"+viewStr);
        else if(f2->GetListOfKeys()->Contains("make"))
            TH2F* h2D_2         = (TH2F*)f2->Get("make/"+hist2+"_in"+viewStr);
        TProfile* ProfileX_2   = h2D_2->ProfileX();
        if(det=="nd" && mcdata=="data") ProfileX_2->Rebin();

        new TCanvas;
        gPad->SetGrid();
        HistogramAttr1D(ProfileX_1, "W (cm)", "Uncorrected PE/cm", ProfileX_1->GetXaxis()->GetBinCenter(1), ProfileX_1->GetXaxis()->GetBinCenter(ProfileX_1->GetXaxis()->GetNbins()), 15, 80, kBlue);
        ProfileX_1->Draw();
        ProfileX_1->SetLineColor(kBlue);
        ProfileX_2->SetLineColor(kRed);
        ProfileX_2->Draw("same");

        TLegend *legend = new TLegend(.12,.62,.64,.90);
        legend->SetBorderSize(0);
        legend->SetFillStyle(0);
        legend->AddEntry((TObject*)0, det+" "+mcdata," ");
        legend->AddEntry(ProfileX_1, analysis1+" "+viewStr+" view");
        legend->AddEntry(ProfileX_2, analysis2+" "+viewStr+" view");
        legend->Draw();

        gPad->Print("Images/uncorrected_PEcm_vs_w_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");

        //ratios
        new TCanvas;
        gPad->SetGrid();
        const TString hrationame = Form("hratio%s_%s_%s", det.Data(), mcdata.Data(), viewstr.Data());
        TH1D* hratio  = new TH1D(hrationame, "",ProfileX_1->GetXaxis()->GetNbins(), ProfileX_1->GetBinCenter(0), ProfileX_1->GetBinCenter(ProfileX_1->GetXaxis()->GetNbins()));
        Ratiohist(ProfileX_1, ProfileX_2, hratio, kRed, 24, "W (cm)", "Ratio", 0.85, 1.10);

        TLegend *legend = new TLegend(.12,.62,.64,.90);
        legend->SetBorderSize(0);
        legend->SetFillStyle(0);
        legend->SetTextSize(0.04);
        legend->AddEntry((TObject*)0, det+" "+mcdata," ");
        legend->AddEntry(hratio,analysis2+" to "+analysis1+" "+viewStr+" view");
        legend->Draw();
        gPad->Print("Images/uncorrected_PEcm_vs_w_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
    }//end of loop over views   
} // end of uncorrected_PE()
//-------------------------------------------------------------------------------------------------------------------------------
void  uncorrected_PE_tricellvstrajectory(char* c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
    uncorrected_PE(c1, c1, det, mcdata, analysis1, analysis2, hist1, hist2);
}//end of uncorrected_PE_tricellvstrajectory()
//-------------------------------------------------------------------------------------------------------------------------------
void thresholdshadowing(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var){
    TFile* f1   = new TFile(c1, "read");
    TFile* f2   = new TFile(c2, "read");
    for(int view = 0; view < 2; ++view){
        const TString viewstr = view ? "y" : "x";
        const TString viewStr = view ? "Y" : "X";

        f1->cd();
        TH1D* correction_1     = (TH1D*)f1->Get("a"+viewstr+"view");
        correction_1->SetLineColor(kBlack);

        f2->cd();
        TH1D* correction_2     = (TH1D*)f2->Get("a"+viewstr+"view");
        correction_2->SetLineColor(kRed);

        new TCanvas;
        gPad->SetGrid();
        if(correction_1->FindObject("stats")){
            TPaveStats* pave1 = (TPaveStats*)correction_1->FindObject("stats");
            pave1->Delete();
        }
        if(correction_2->FindObject("stats")){
            TPaveStats* pave2 = (TPaveStats*)correction_2->FindObject("stats");
            pave2->Delete();
        }
        correction_1->Draw();
        correction_2->Draw("same");
        if(var=="W") TString labelx = "W (cm)";
        HistogramAttr1D(correction_1, labelx.Data(), "correction factor", correction_1->GetXaxis()->GetBinCenter(1), correction_1->GetXaxis()->GetBinCenter(correction_1->GetXaxis()->GetNbins()), 0.5, 2.0, kBlack);

        TLegend *legend = new TLegend(.12,.62,.64,.90);
        legend->SetBorderSize(0);
        legend->SetFillStyle(0);
        legend->AddEntry((TObject*)0, det+" "+mcdata," ");
        legend->AddEntry(correction_1, analysis1+" "+viewStr+" view");
        legend->AddEntry(correction_2, analysis2+" "+viewStr+" view");
        legend->Draw();
        gPad->Print("Images/thresholdshadowing_vs_"+var+"_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");

        //ratios
        new TCanvas;
        gPad->SetGrid();
        const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), viewstr.Data(), var.Data());
        TH1D* hratio  = new TH1D(hrationame, "",correction_1->GetXaxis()->GetNbins(), correction_1->GetBinCenter(0), correction_1->GetBinCenter(correction_1->GetXaxis()->GetNbins()));
        Ratiohist(correction_1, correction_2, hratio, kRed, 24, labelx.Data(), "Correction factor ratio", 1.8, 1.2);

        TLegend *legend = new TLegend(.12,.62,.64,.90);
        legend->SetBorderSize(0);
        legend->SetFillStyle(0);
        legend->SetTextSize(0.04);
        legend->AddEntry((TObject*)0, det+" "+mcdata," ");
        legend->AddEntry(hratio,analysis2+" to "+analysis1+" "+viewStr+" view");
        legend->Draw();
        gPad->Print("Images/thresholdshadowinge_vs_"+var+"_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
    }//end of loop over views
}//end of thresholdshadowing()
//-------------------------------------------------------------------------------------------------------------------------------
void corrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1,  TString analysis2, TString plane, TString cell){

    if(nddatamuoncatcher){
        if(plane.Atoi() > 3 || plane.Atoi() < 2){
            std::cout << "It is a ND muon catcher data file and has only 002, 003 planes. Plane " << plane << " is out of range. " << std::endl;
            abort();
        }
        if(cell.Atoi() > 95 || cell.Atoi() < 0){
            std::cout << "It is a Data file. Cell " << cell << " is out of range. " << std::endl;
            abort();
        }
    }
    if(det.Data()=="nd"){
        if(mcdata.Data()=="mc"){
            if(plane.Atoi() > 3 || plane.Atoi() < 0){
                std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
                abort();
            }
            if(cell.Atoi() > 95 || cell.Atoi() < 0){
                std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
                abort();
            }
        }//end of MC condition
        if(mcdata.Data()=="data"){
            if(plane.Atoi() > 191 || plane.Atoi() < 0){
                std::cout << "It is a ND data file and has only 191 main body planes. Plane " << plane << " out of range. " << std::endl;
                abort();
            }
            if(cell.Atoi() > 95 || cell.Atoi() < 0){
                std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
                abort();
            }
        }//end of data condition
    }//end of the ND condition

    if(det.Data()=="fd"){
        if(mcdata.Data()=="mc"){
            if(plane.Atoi() > 1 || plane.Atoi() < 0){
                std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
                abort();
            }
            if(cell.Atoi() > 383 || cell.Atoi() < 0){
                std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
                abort();
            }
        }//end of MC condition
        else if(mcdata.Data()=="data"){
            if(plane.Atoi() > 895 || plane.Atoi() < 0){
                std::cout << "It is a Data file. Plane " << plane << " out of range. " << std::endl;
                abort();
            }
            if(cell.Atoi() > 383 || cell.Atoi() < 0){
                std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
                abort();
            }
        }//end of data condition
    }//end of FD condition

    TFile* f1   = new TFile(c1, "read");
    TFile* f2   = new TFile(c2, "read");

    TF1*   fit1;
    TF1*   fit2;

    f1->cd();
    TH1F* h1        = (TH1F*)f1->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
    Color_t color1  = kBlack;
    if(mcdata=="mc"){
        if(plane.Atoi()==0)      fit1  = (TF1*)h1->FindObject("fitX");
        else if(plane.Atoi()==1) fit1  = (TF1*)h1->FindObject("fitY");
        else if(plane.Atoi()==2) fit1  = (TF1*)h1->FindObject("fitMuX");
        else if(plane.Atoi()==3) fit1  = (TF1*)h1->FindObject("fitMuY");
        if(fit1) fit1->SetLineColor(color1);
    }//end of MC condition
    else if(mcdata=="data"){
        if(nddatamuoncatcher){
            if(plane.Atoi()==2) fit1        = (TF1*)h1->FindObject("fitMuX");
            else if(plane.Atoi()==3) fit1   = (TF1*)h1->FindObject("fitMuY");
            if(fit1) fit1->SetLineColor(color1);
        }//end of muon catcher condition
        else if(plane.Atoi()%2==0) fit1     = (TF1*)h1->FindObject("fitY");
        else if(plane.Atoi()%2!=0) fit1     = (TF1*)h1->FindObject("fitX");
        if(fit1) fit1->SetLineColor(color1);
    }//end of data condition

    f2->cd();
    TH1F* h2        = (TH1F*)f2->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
    Color_t color2  = kRed;
    if(mcdata=="mc"){
        if(plane.Atoi()==0)      fit2  = (TF1*)h2->FindObject("fitX");
        else if(plane.Atoi()==1) fit2  = (TF1*)h2->FindObject("fitY");
        else if(plane.Atoi()==2) fit2  = (TF1*)h2->FindObject("fitMuX");
        else if(plane.Atoi()==3) fit2  = (TF1*)h2->FindObject("fitMuY");
        if(fit2) fit2->SetLineColor(color2);
    }//end of MC condition
    else if(mcdata=="data"){
        if(nddatamuoncatcher){
            if(plane.Atoi()==2) fit2  = (TF1*)h2->FindObject("fitMuX");
            else if(plane.Atoi()==3) fit2  = (TF1*)h2->FindObject("fitMuY");
            if(fit2) fit2->SetLineColor(color2);
        }//end of muon catcher condition
        else if(plane.Atoi()%2==0) fit2       = (TF1*)h2->FindObject("fitY");
        else if(plane.Atoi()%2!=0) fit2  = (TF1*)h2->FindObject("fitX");
        if(fit2) fit2->SetLineColor(color2);
    }//end of data condition

    new TCanvas;
    gPad->SetGrid();
    h1->Draw();
    HistogramAttr1D(h1, "W (cm)", "Mean PE/cm", h1->GetXaxis()->GetBinCenter(1), h1->GetXaxis()->GetBinCenter(h1->GetXaxis()->GetNbins()), 0., 150., color1);
    h1->SetLineColor(color1);
    h1->SetMarkerColor(color1);
    h1->SetMarkerStyle(24);
    h2->Draw("same");
    h2->SetLineColor(color2);
    h2->SetMarkerColor(color2);
    h2->SetMarkerStyle(24);

    TLegend *legend = new TLegend(.12,.62,.64,.90);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.04);
    legend->AddEntry((TObject*)0, det+" "+mcdata," ");
    legend->AddEntry(h1, analysis1);
    legend->AddEntry(h2, analysis2);
    legend->Draw();
    gPad->Print("Images/Atten_Fit_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+plane+"_"+cell+".pdf");

    //ratio of data points
    new TCanvas;
    gPad->SetGrid();
    const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
    TH1D* hratio  = new TH1D(hrationame, "",h1->GetXaxis()->GetNbins(), h1->GetBinCenter(0), h1->GetBinCenter(h1->GetXaxis()->GetNbins()));
    Ratiohist(h1, h2, hratio, kRed, 24, "W (cm)", "Ratio of Mean PE/cm data points", 0.3, 1.0);

    TLegend *legend = new TLegend(.12,.62,.64,.90);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.04);
    legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
    legend->AddEntry(hratio,analysis2+" to "+analysis1);
    legend->Draw();
    gPad->Print("Images/Atten_Fit_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_plane"+plane+"_cell"+cell+".pdf");

    //ratio of fitting curve
    new TCanvas;
    gPad->SetGrid();
    const TString hratiofitname = Form("hratiofit%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
    if(det=="nd"){
        float x_limit   = 250.;
        int x_bins      = 2*x_limit;
    }
    else if(det=="fd"){
        float x_limit   = 800.;
        int x_bins      = 2*x_limit;
    }
    TH1D* hfit  = new TH1D("hfit", "", x_bins, -x_limit, +x_limit);
    HistogramAttr1D(hfit, "W (cm)", "Ratio of Mean PE/cm fitting curves", -x_limit, +x_limit, 0.55, 0.68, kWhite);
    hfit->Draw();
    TH1D* hratiofit  = new TH1D(hratiofitname, "", x_bins, -x_limit, +x_limit);
    Ratiofit(fit1, fit2, hratiofit, plane.Data(), cell.Data(), color2, 24, "W (cm)", "Ratio of Mean PE/cm fitting curves", 0.55, 0.68);
    TLegend *legend = new TLegend(.12,.62,.64,.90);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.04);
    legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
    legend->AddEntry(hratiofit,analysis2+" to "+analysis1);
    legend->Draw();
    gPad->Print("Images/Atten_Fit_fitratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_plane"+plane+"_cell"+cell+".pdf");
}//end of corrected_PE()
//---------------------------------------------------------------------------------------------------
void calibrationconstants(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2){
    ifstream in_1; 
    ifstream in_2;

    ofstream out;

    FILE* fout = fopen("uncalibrated_channels_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+".csv","w");
    fprintf(fout, "%s%50s%s\n", analysis1.Data(), " ", analysis2.Data());
    fprintf(fout, "channel,coeff_exp,atten_length,background,chisq, %9s, channel,coeff_exp,atten_length,background,chisq\n", " ");

    in_1.open(c1);
    in_2.open(c2);


    //x=coeff_exp, y=atten_length, z=background, a=edgx_low, b=edgx_high, c=coeff_low
    Float_t x_1, y_1, z_1, a_1, b_1, c_1, d_1, e_1, g_1;
    Float_t x_2, y_2, z_2, a_2, b_2, c_2, d_2, e_2, g_2;
    Float_t     chisq               = 0.2;

    Int_t       channels;
    Int_t       totalchannels1      = 0;
    Int_t       totalchannels2      = 0;
    Int_t       calibratedchannels1 = 0;
    Int_t       calibratedchannels2 = 0;
    Int_t       f_1;
    Int_t       f_2;
    Int_t       bins1;

    Double_t    limit1D     = 3.;
    Double_t    limit2D     = 1000.;

    if(det=="nd" && mcdata=="mc"    )   { channels  = 3063;     }
    if(det=="nd" && mcdata=="data"  )   { channels  = 1003063;  }
    if(det=="fd" && mcdata=="mc"    )   { channels  = 1383;     }
    if(det=="fd" && mcdata=="data"  )   { channels  = 895383;   }


    TH1F* hcoeffexp     = new TH1F("hcoeffexp",     "coeffexp difference distribution",     120, -limit1D, +limit1D    );
    TH1F* hattenlength  = new TH1F("hattenlength",  "attenlength difference distribution",  600, -limit1D, +limit1D    );
    TH1F* hbackground   = new TH1F("hbackground",   "background difference distribution",   120, -limit1D, +limit1D    );
    TH1F* hchisq        = new TH1F("hchisq",        "chisq difference distribution",        600, -limit1D, +limit1D    );

    TH2F* h2coeffexp    = new TH2F("h2coeffexp",   " " , 2*(Int_t)limit2D, -limit2D, limit2D,  2*(Int_t)limit2D, -limit2D, limit2D );
    TH2F* h2attenlength = new TH2F("h2attenlength"," " , 2*(Int_t)limit2D, -limit2D, limit2D,  2*(Int_t)limit2D, -limit2D, limit2D );
    TH2F* h2background  = new TH2F("h2background", " " , 2*(Int_t)limit2D, -limit2D, limit2D,  2*(Int_t)limit2D, -limit2D, limit2D );
    TH2F* h2chisq       = new TH2F("h2chisq"," ",                   200,0.,2.,    200,0.,2. );

    while(1){
        in_1 >> x_1 >> y_1 >> z_1 >> a_1 >> b_1 >> c_1 >> d_1 >> e_1 >> f_1 >> g_1;
        in_2 >> x_2 >> y_2 >> z_2 >> a_2 >> b_2 >> c_2 >> d_2 >> e_2 >> f_2 >> g_2;
        if(!in_1.good()) break;
        if(!in_2.good()) break;

        if(f_1<=channels && e_1 < chisq && e_2 < chisq){
            hcoeffexp->Fill((x_1 - x_2)/x_1);
            hattenlength->Fill((y_1 - y_2)/y_1);
            hbackground->Fill((z_1 - z_2)/z_1);

            h2coeffexp->Fill(x_2, x_1);
            h2attenlength->Fill(y_2, y_1);
            h2background->Fill(z_2, z_1);
        }//end of loop over channels and chisq conditions.
        if(f_1<=channels){
            totalchannels1++;
            totalchannels2++;
            if(e_1<chisq) calibratedchannels1++;
            if(e_2<chisq) calibratedchannels2++;

            hchisq->Fill((e_1 - e_2)/e_1);
            h2chisq->Fill(e_2, e_1);

            if(e_1>=chisq || e_2>=chisq){ 
                fprintf(fout, "%i,%g,%g,%g,%g,%20s,%i,%g,%g,%g,%g\n",f_1,x_1,y_1,z_1,e_1, " ", f_2,x_2,y_2,z_2,e_2);
                out << f_1 << setw(22) << x_1 << setw(22) << y_1 << setw(22) << z_1 << setw(22) << e_1 << f_2 << setw(22) << x_2 << setw(22) << y_2 << setw(22) << z_2 << setw(22) << e_2 << endl;
            }
        }//end of only channels condition for chisq plots. For chisq plots we have to consider all channels.
    }//end of while condition.

    cout << " " << endl;
    cout << "Total channels in "+analysis1+" = " << totalchannels1 << " and calibrated channels = " << calibratedchannels1 <<  " = " << ((float)calibratedchannels1/(float)totalchannels1)*100. << " %." << endl;
    cout << "Total channels in "+analysis2+" = " << totalchannels2 << " and calibrated channels = " << calibratedchannels2 << " = " << ((float)calibratedchannels2/(float)totalchannels2)*100. << " %." << endl;
    cout << " " << endl;

    //coeffexp plots
    new TCanvas;
    gPad->SetGrid();
    HistogramAttr1D(hcoeffexp, "#Deltacoeffexp/coeffexp", "# of channels", 0.5*hcoeffexp->GetXaxis()->GetBinCenter(1), 0.5*hcoeffexp->GetXaxis()->GetBinCenter(hcoeffexp->GetXaxis()->GetNbins()), 0., 7e5, kBlack);
    hcoeffexp->Draw();
    TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.05);
    legend->AddEntry((TObject*)0, det+" "+mcdata," ");
    legend->AddEntry(hcoeffexp,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
    legend->Draw();
    gPad->Print("Images/Delta_coeffexp_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    new TCanvas;
    gPad->SetGrid();
    gPad->SetLogz();
    HistogramAttr2D(h2coeffexp, "coeffexp_{"+analysis2+"}", "coeffexp_{"+analysis1+"}", " ", 0, 100, 0, 100);
    h2coeffexp->Draw("colz");
    gPad->Print("Images/coeffexp2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    //attenlength plots
    new TCanvas;
    gPad->SetGrid();
    HistogramAttr1D(hattenlength, "#Deltaattenlength/attenlength", "# of channels", 0.5*hattenlength->GetXaxis()->GetBinCenter(1), 0.5*hattenlength->GetXaxis()->GetBinCenter(hattenlength->GetXaxis()->GetNbins()), 0., 166500, kBlack);
    hattenlength->Draw();
    TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.05);
    legend->AddEntry((TObject*)0, det+" "+mcdata," ");
    legend->AddEntry(hattenlength,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
    legend->Draw();
    gPad->Print("Images/Delta_attenlength_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    new TCanvas;
    gPad->SetGrid();
    gPad->SetLogz();
    HistogramAttr2D(h2attenlength, "attenlength_{"+analysis2+"} (cm)", "attenlength_{"+analysis1+"} (cm)", " ", 0, 1000, 0, 1000);
    h2attenlength->Draw("colz");
    gPad->Print("Images/attenlength2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    //backgrounds plots
    new TCanvas;
    gPad->SetGrid();
    HistogramAttr1D(hbackground, "#Deltabackground/background", "# of channels", 0.5*hbackground->GetXaxis()->GetBinCenter(1), 0.5*hbackground->GetXaxis()->GetBinCenter(hbackground->GetXaxis()->GetNbins()), 0., 3e4, kBlack);
    hbackground->Draw();
    TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.05);
    legend->AddEntry((TObject*)0, det+" "+mcdata," ");
    legend->AddEntry(hbackground,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
    legend->Draw();
    gPad->Print("Images/Delta_background_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    new TCanvas;
    gPad->SetGrid();
    gPad->SetLogz();
    HistogramAttr2D(h2background, "background_{"+analysis2+"}", "background_{"+analysis1+"}", " ", -50, 50, -50, 50);
    h2background->Draw("colz");
    gPad->Print("Images/background2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    //chisq plots
    new TCanvas;
    gPad->SetGrid();
    HistogramAttr1D(hchisq, "#Deltachisq/chisq", "# of channels", 0.5*hchisq->GetXaxis()->GetBinCenter(1), 0.5*hchisq->GetXaxis()->GetBinCenter(hchisq->GetXaxis()->GetNbins()), 0., 4200., kBlack);
    hchisq->Draw();
    TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
    legend->SetBorderSize(0);
    legend->SetFillStyle(0);
    legend->SetTextSize(0.05);
    legend->AddEntry((TObject*)0, det+" "+mcdata," ");
    legend->AddEntry(hchisq,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
    legend->Draw();
    gPad->Print("Images/Delta_chisq_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    new TCanvas;
    gPad->SetGrid();
    gPad->SetLogz();
    HistogramAttr2D(h2chisq, "#chi^{2}_{"+analysis2+"}", "#chi^{2}_{"+analysis1+"}", " ", 0, 0.5, 0, 0.5);
    h2chisq->Draw("colz");
    gPad->Print("Images/chisq2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");

    in_1.close();
    in_2.close();
    out.close();
    fclose(fout);
}//end of calibrationconstants()
//---------------------------------------------------------------------------------------------------