//************************************************
// Author: Federica Lionetto
// Created on: 03/12/2014
//************************************************

/*
   CMNSubMean is one of the CMN subtraction algorithms.

   Pedestal subtracted ADC counts ---> Calculate mean and rms ---> Exclude signal (requiring |ADC-mean|<factor*rms) ---> Calculate mean and rms ---> Exclude signal (requiring |ADC-mean|<factor*rms) ---> Calculate mean ---> Subtract mean ---> Pedestal and CMN subtracted ADC counts.

   It loops over events, creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex_1evt), calculates the mean and the rms, ignores strips that may potentially contain a signal and creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig_1evt), calculates the mean and the rms, ignores strips that may potentially contain a signal and creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig2_1evt), calculates the average, subtracts it to every strip, and creates the event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex_1evt).
   If there is no common mode noise, the pedestal and the pedestal and CMN versions will look the same.

   It creates the distribution of the CMN (hCMNBeetlex).

   It creates the distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex) and the distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex), comparing the two.

   It creates the distribution of the RMS of the ADC count distribution for the following cases:
   - pedestal subtracted ADC counts;
   - pedestal subtracted ADC counts, after the first outlier removal;
   - pedestal subtracted ADC counts, after the second outlier removal.

   In beam off data the ADC count distribution after CMN subtraction is not symmetric.
   This might be explained by muons coming from the beam line and acting as "signal".
   To check this hypothesis, it is interesting to look at the distribution of the Beetle channels contributing to the tail on the right (TH1F *hchTailBeetlex) and the distribution of the event numbers where this happens (TH1F *hevTailBeetlex).

   The configuration of the CMN algorithm is defined by factor.
   The other arguments are:
   - the number of one event to display, to check that the algorithm is working properly;
   - the string corresponding to the input ROOT file, where to get the pedestal subtracted ADC counts from;
   - the pointer to the output ROOT file, where to save the output of this algorithm;
   - the string corresponding to the path where to save the pictures.

   CMNSubMean2 will create a directory in the output ROOT file called "CMNSubMean2", where the tree (with CMN1, CMN2, and ADCPedCMNSub) and the histograms, graphs, and canvases will be stored.
   */

// Header guard.
#ifndef __CMNSUBMEAN2_C_INCLUDED__
#define __CMNSUBMEAN2_C_INCLUDED__

void CMNSubMean2(float factor, int sample, string input_ROOT, TFile *output, TString path_to_figures); 

void CMNSubMean2(float factor, int sample, string input_ROOT, TFile *output, TString path_to_figures)
{
    cout << "**************************************************" << endl;
    cout << "Computing CMN according to the mean 2 algorithm..." << endl;
    cout << "**************************************************" << endl;

    cout << "**************************************************" << endl;
    cout << "Configuration" << endl;
    cout << "factor: " << factor << endl;
    cout << "sample: " << sample << endl;
    cout << "**************************************************" << endl;

    // Open input ROOT file.
    cout << "Open input ROOT file: " << input_ROOT << endl;
    TFile *input = TFile::Open(TString(input_ROOT));

    input->cd();
    TTree *EventInfo = (TTree *)input->Get("EventInfo");
    int NEvents = EventInfo->GetEntries();

    TString dirCMNSubMean2 = "CMNSubMean2";
    output->cd();
    output->mkdir(dirCMNSubMean2);
    output->cd(dirCMNSubMean2);

    TString alg_path_to_figures = path_to_figures+dirCMNSubMean2;
    TString alg_path_to_make_figures = "mkdir -p "+alg_path_to_figures;
    system(alg_path_to_make_figures.Data());

    // Copy the EventInfo tree from the input to the output file and add some more information.
    float CMN1 = 0.;
    float CMN2 = 0.;
    std::vector<float> ADCPedCMNSub; // Raw ADC - pedestal - CMN.

    TTree *EventInfoProcessed = EventInfo->CloneTree(0);
    std::vector<float> *ADCProcessed = 0;
    TBranch *b_ADCProcessed = 0;
    EventInfo->SetBranchAddress("ADCProcessed",&ADCProcessed,&b_ADCProcessed);
    EventInfoProcessed->Branch("CMN1",&CMN1);
    EventInfoProcessed->Branch("CMN2",&CMN2);
    EventInfoProcessed->Branch("ADCPedCMNSub",&ADCPedCMNSub);

    // Event-by-event histograms.

    // Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex_1evt).
    TH1F *hADCPedSubBeetle1_1evt = new TH1F("hADCPedSubBeetle1_1evt","",400,-200,200); // Modify it.
    TCanvas *cADCPedSubBeetle1_1evt = new TCanvas("cADCPedSubBeetle1_1evt","",400,300);
    InitHist(hADCPedSubBeetle1_1evt,"Pedestal subtracted ADC counts on Beetle 1","ADC counts","");

    TH1F *hADCPedSubBeetle2_1evt = new TH1F("hADCPedSubBeetle2_1evt","",400,-200,200); //Modify it.
    TCanvas *cADCPedSubBeetle2_1evt = new TCanvas("cADCPedSubBeetle2_1evt","",400,300);
    InitHist(hADCPedSubBeetle2_1evt,"Pedestal subtracted ADC counts on Beetle 2","ADC counts","");

    // Ignore strips that may potentially contain a signal.
    // Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig_1evt).

    // First iteration.
    TH1F *hADCPedSubBeetle1NoSig_1evt = new TH1F("hADCPedSubBeetle1NoSig_1evt","",400,-200,200); // Modify it.
    TCanvas *cADCPedSubBeetle1NoSig_1evt = new TCanvas("cADCPedSubBeetle1NoSig_1evt","",400,300);
    InitHist(hADCPedSubBeetle1NoSig_1evt,"Pedestal subtracted ADC counts on Beetle 1 - excluding signal","ADC counts","");

    TH1F *hADCPedSubBeetle2NoSig_1evt = new TH1F("hADCPedSubBeetle2NoSig_1evt","",400,-200,200); //Modify it.
    TCanvas *cADCPedSubBeetle2NoSig_1evt = new TCanvas("cADCPedSubBeetle2NoSig_1evt","",400,300);
    InitHist(hADCPedSubBeetle2NoSig_1evt,"Pedestal subtracted ADC counts on Beetle 2 - excluding signal","ADC counts","");

    // Second iteration.
    TH1F *hADCPedSubBeetle1NoSig2_1evt = new TH1F("hADCPedSubBeetle1NoSig2_1evt","",400,-200,200); // Modify it.
    TCanvas *cADCPedSubBeetle1NoSig2_1evt = new TCanvas("cADCPedSubBeetle1NoSig2_1evt","",400,300);
    InitHist(hADCPedSubBeetle1NoSig2_1evt,"Pedestal subtracted ADC counts on Beetle 1 - excluding signal","ADC counts","");

    TH1F *hADCPedSubBeetle2NoSig2_1evt = new TH1F("hADCPedSubBeetle2NoSig2_1evt","",400,-200,200); //Modify it.
    TCanvas *cADCPedSubBeetle2NoSig2_1evt = new TCanvas("cADCPedSubBeetle2NoSig2_1evt","",400,300);
    InitHist(hADCPedSubBeetle2NoSig2_1evt,"Pedestal subtracted ADC counts on Beetle 2 - excluding signal","ADC counts","");

    // Event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex_1evt).
    TH1F *hADCPedCMNSubBeetle1_1evt = new TH1F("hADCPedCMNSubBeetle1_1evt","",400,-200,200); // Modify it.
    TCanvas *cADCPedCMNSubBeetle1_1evt = new TCanvas("cADCPedCMNSubBeetle1_1evt","",400,300);
    InitHist(hADCPedCMNSubBeetle1_1evt,"Pedestal and CMN subtracted ADC counts on Beetle 1","ADC counts","");

    TH1F *hADCPedCMNSubBeetle2_1evt = new TH1F("hADCPedCMNSubBeetle2_1evt","",400,-200,200); //Modify it.
    TCanvas *cADCPedCMNSubBeetle2_1evt = new TCanvas("cADCPedCMNSubBeetle2_1evt","",400,300);
    InitHist(hADCPedCMNSubBeetle2_1evt,"Pedestal and CMN subtracted ADC counts on Beetle 2","ADC counts","");

    // Histograms containing the information on all events.

    // CMN1.
    TH1F *hCMNBeetle1 = new TH1F("hCMNBeetle1","",100,-50,50); // Modify it.
    TCanvas *cCMNBeetle1 = new TCanvas("cCMNBeetle1","",400,300);
    InitHist(hCMNBeetle1,"CMN on Beetle 1","ADC counts","");

    // CMN2.
    TH1F *hCMNBeetle2 = new TH1F("hCMNBeetle2","",100,-50,50); // Modify it.
    TCanvas *cCMNBeetle2 = new TCanvas("cCMNBeetle2","",400,300);
    InitHist(hCMNBeetle2,"CMN on Beetle 2","ADC counts","");

    // Distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex).
    TH1F *hADCPedSubBeetle1 = new TH1F("hADCPedSubBeetle1","",400,-200,200);
    InitHist(hADCPedSubBeetle1,"Pedestal subtracted ADC counts on Beetle 1","ADC counts","");
    TH1F *hADCPedSubBeetle2 = new TH1F("hADCPedSubBeetle2","",400,-200,200);
    InitHist(hADCPedSubBeetle2,"Pedestal subtracted ADC counts on Beetle 2","ADC counts","");

    // Distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex).
    TH1F *hADCPedCMNSubBeetle1 = new TH1F("hADCPedCMNSubBeetle1","",400,-200,200);
    InitHist(hADCPedCMNSubBeetle1,"Pedestal and CMN subtracted ADC counts on Beetle 1","ADC counts","");
    TH1F *hADCPedCMNSubBeetle2 = new TH1F("hADCPedCMNSubBeetle2","",400,-200,200);
    InitHist(hADCPedCMNSubBeetle2,"Pedestal and CMN subtracted ADC counts on Beetle 2","ADC counts","");

    TCanvas *cADCCompBeetle1 = new TCanvas("cADCCompBeetle1","",400,300);
    TCanvas *cADCCompBeetle2 = new TCanvas("cADCCompBeetle2","",400,300);

    // Distribution of the RMS of the ADC count distribution for the following cases:
    // - pedestal subtracted ADC counts;
    // - pedestal subtracted ADC counts, after the first outlier removal;
    // - pedestal subtracted ADC counts, after the second outlier removal.
    TH1F *hRMSPedSubBeetle1 = new TH1F("hRMSPedSubBeetle1","",90,0,30);
    TH1F *hRMSPedSubBeetle1NoSig = new TH1F("hRMSPedSubBeetle1NoSig","",90,0,30);
    TH1F *hRMSPedSubBeetle1NoSig2 = new TH1F("hRMSPedSubBeetle1NoSig2","",90,0,30);
    TH1F *hRMSPedSubBeetle2 = new TH1F("hRMSPedSubBeetle2","",90,0,30);
    TH1F *hRMSPedSubBeetle2NoSig = new TH1F("hRMSPedSubBeetle2NoSig","",90,0,30);
    TH1F *hRMSPedSubBeetle2NoSig2 = new TH1F("hRMSPedSubBeetle2NoSig2","",90,0,30);

    InitHist(hRMSPedSubBeetle1,"RMS of pedestal subtracted ADC counts on Beetle 1","RMS (ADC counts)","");
    InitHist(hRMSPedSubBeetle1NoSig,"RMS of pedestal subtracted ADC counts on Beetle 1","RMS (ADC counts)","");
    InitHist(hRMSPedSubBeetle1NoSig2,"RMS of pedestal subtracted ADC counts on Beetle 1","RMS (ADC counts)","");
    InitHist(hRMSPedSubBeetle2,"RMS of pedestal subtracted ADC counts on Beetle 2","RMS (ADC counts)","");
    InitHist(hRMSPedSubBeetle2NoSig,"RMS of pedestal subtracted ADC counts on Beetle 2","RMS (ADC counts)","");
    InitHist(hRMSPedSubBeetle2NoSig2,"RMS of pedestal subtracted ADC counts on Beetle 2","RMS (ADC counts)","");

    TCanvas *cRMSBeetle1 = new TCanvas("cRMSBeetle1","",400,300);
    TCanvas *cRMSBeetle2 = new TCanvas("cRMSBeetle2","",400,300);

    // In beam off data the ADC count distribution after CMN subtraction is not symmetric.
    // This might be explained by muons coming from the beam line and acting as "signal".
    // To check this hypothesis, it is interesting to look at the distribution of the Beetle channels contributing to the tail on the right (TH1F *hchTailBeetlex) and the distribution of the event numbers where this happens (TH1F *hevTailBeetlex).
    TH1F *hchTailBeetle1 = new TH1F("hcdTailBeetle1","",NBeetle,0,NBeetle);
    TH1F *hchTailBeetle2 = new TH1F("hcdTailBeetle2","",NBeetle,NBeetle,N);

    InitHist(hchTailBeetle1,"Asymmetric tail on Beetle 1","Beetle channel","");
    InitHist(hchTailBeetle2,"Asymmetric tail on Beetle 2","Beetle channel","");

    TCanvas *cchTailBeetle1 = new TCanvas("cchTailBeetle1","",400,300);
    TCanvas *cchTailBeetle2 = new TCanvas("cchTailBeetle2","",400,300);

    TH1F *hevTailBeetle1 = new TH1F("hevTailBeetle1","",50,0,NEvents+1);
    TH1F *hevTailBeetle2 = new TH1F("hevTailBeetle2","",50,0,NEvents+1);

    InitHist(hevTailBeetle1,"Asymmetric tail on Beetle 1","Event number","");
    InitHist(hevTailBeetle2,"Asymmetric tail on Beetle 2","Event number","");

    TCanvas *cevTailBeetle1 = new TCanvas("cevTailBeetle1","",400,300);
    TCanvas *cevTailBeetle2 = new TCanvas("cevTailBeetle2","",400,300);

    // Add the histograms with the pedestal and CMN subtracted distribution of the ADC counts of each Beetle channel (TH1F *hADCPedCMNSub[N]) to the output file.
    TH1F *hADCPedCMNSub[N];
    for (int iChannel=0;iChannel<N;++iChannel) {
        hADCPedCMNSub[iChannel] = new TH1F(Form("hADCPedCMNSub%d",iChannel),"",1024,0,1024); 
        InitHist(hADCPedCMNSub[iChannel],Form("Pedestal and CMN subtracted ADC counts on Beetle channel %d",iChannel),"ADC counts","");
    }
    
    // Beetle 1.

    float meanBeetle1_temp;
    float rmsBeetle1_temp;

    // Beetle2.

    float meanBeetle2_temp;
    float rmsBeetle2_temp;

    // Fill histograms.
    for (int iEvent=0;iEvent<NEvents;++iEvent) {
        if (iEvent%1000==0)
            cout << "Processing event " << iEvent << endl;

        // Reset.
        if (iEvent>0)
        {
            ADCPedCMNSub.clear();
            hADCPedSubBeetle1_1evt->Reset();
            hADCPedSubBeetle2_1evt->Reset();
            hADCPedSubBeetle1NoSig_1evt->Reset();
            hADCPedSubBeetle2NoSig_1evt->Reset();
            hADCPedSubBeetle1NoSig2_1evt->Reset();
            hADCPedSubBeetle2NoSig2_1evt->Reset();
            hADCPedCMNSubBeetle1_1evt->Reset();
            hADCPedCMNSubBeetle2_1evt->Reset();
        }

        EventInfo->GetEntry(iEvent);

        // Fill histograms for Beetle 1.
        for (int iChannel=0;iChannel<NBeetle;++iChannel) {
            if (maskBeetle1[iChannel] == 0)
            {
                hADCPedSubBeetle1->Fill((float)ADCProcessed->at(iChannel));
                hADCPedSubBeetle1_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Exclude signal.
        meanBeetle1_temp = hADCPedSubBeetle1_1evt->GetMean();
        rmsBeetle1_temp = hADCPedSubBeetle1_1evt->GetRMS();

        hRMSPedSubBeetle1->Fill(rmsBeetle1_temp);

        // Fill histograms again for Beetle 1.
        for (int iChannel=0;iChannel<NBeetle;++iChannel) {
            if (maskBeetle1[iChannel] == 0)
            {
                if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1_temp)<factor*rmsBeetle1_temp)
                    hADCPedSubBeetle1NoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Exclude signal.
        meanBeetle1_temp = hADCPedSubBeetle1NoSig_1evt->GetMean();
        rmsBeetle1_temp = hADCPedSubBeetle1NoSig_1evt->GetRMS();

        hRMSPedSubBeetle1NoSig->Fill(rmsBeetle1_temp);

        // Fill histograms again for Beetle 1.
        for (int iChannel=0;iChannel<NBeetle;++iChannel) {
            if (maskBeetle1[iChannel] == 0)
            {
                if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1_temp)<factor*rmsBeetle1_temp)
                    hADCPedSubBeetle1NoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Calculate CMN1.
        CMN1 = hADCPedSubBeetle1NoSig2_1evt->GetMean();
        rmsBeetle1_temp = hADCPedSubBeetle1NoSig2_1evt->GetRMS();

        hRMSPedSubBeetle1NoSig2->Fill(rmsBeetle1_temp);

        hCMNBeetle1->Fill(CMN1);

        // Fill histograms for Beetle 2.
        for (int iChannel=NBeetle;iChannel<N;++iChannel) {
            if (maskBeetle2[iChannel-NBeetle] == 0)
            {
                hADCPedSubBeetle2->Fill((float)ADCProcessed->at(iChannel));
                hADCPedSubBeetle2_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Exclude signal.
        meanBeetle2_temp = hADCPedSubBeetle2_1evt->GetMean();
        rmsBeetle2_temp = hADCPedSubBeetle2_1evt->GetRMS();

        hRMSPedSubBeetle2->Fill(rmsBeetle2_temp);

        // Fill histograms again for Beetle 2.
        for (int iChannel=NBeetle;iChannel<N;++iChannel) {
            if (maskBeetle2[iChannel-NBeetle] == 0)
            {
                if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2_temp)<factor*rmsBeetle2_temp)
                    hADCPedSubBeetle2NoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Exclude signal.
        meanBeetle2_temp = hADCPedSubBeetle2NoSig_1evt->GetMean();
        rmsBeetle2_temp = hADCPedSubBeetle2NoSig_1evt->GetRMS();

        hRMSPedSubBeetle2NoSig->Fill(rmsBeetle2_temp);

        // Fill histograms again for Beetle 2.
        for (int iChannel=NBeetle;iChannel<N;++iChannel) {
            if (maskBeetle2[iChannel-NBeetle] == 0)
            {
                if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2_temp)<factor*rmsBeetle2_temp)
                    hADCPedSubBeetle2NoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Calculate CMN2.
        CMN2 = hADCPedSubBeetle2NoSig2_1evt->GetMean();
        rmsBeetle2_temp = hADCPedSubBeetle2NoSig2_1evt->GetRMS();

        hRMSPedSubBeetle2NoSig2->Fill(rmsBeetle2_temp);

        hCMNBeetle2->Fill(CMN2);

        // Subtract CMN1.
        for (int iChannel=0;iChannel<NBeetle;++iChannel) {
            ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN1);
            hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN1);
            if (maskBeetle1[iChannel] == 0)
            {
                hADCPedCMNSubBeetle1->Fill((float)ADCProcessed->at(iChannel)-CMN1);
                hADCPedCMNSubBeetle1_1evt->Fill((float)ADCProcessed->at(iChannel)-CMN1);
                if (((float)ADCProcessed->at(iChannel)-CMN1) > 30.)
                {
                    cout << "The event #" << iEvent << " has at least one Beetle channels (" << iChannel << ") contributing to the tail" << endl;
                    hchTailBeetle1->Fill(iChannel);
                    hevTailBeetle1->Fill(iEvent);
                }
            }
        }

        // Subtract CMN2.
        for (int iChannel=NBeetle;iChannel<N;++iChannel) {
            ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN2);
            hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN2);
            if (maskBeetle2[iChannel-NBeetle] == 0)
            {
                hADCPedCMNSubBeetle2->Fill((float)ADCProcessed->at(iChannel)-CMN2);
                hADCPedCMNSubBeetle2_1evt->Fill((float)ADCProcessed->at(iChannel)-CMN2);
                if (((float)ADCProcessed->at(iChannel)-CMN2) > 30.)
                {
                    cout << "The event #" << iEvent << " has at least one Beetle channels (" << iChannel << ") contributing to the tail" << endl;
                    hchTailBeetle2->Fill(iChannel);
                    hevTailBeetle2->Fill(iEvent);
                }
            }
        }

        EventInfoProcessed->Fill();

        if (iEvent==sample) {
            DrawHist(cADCPedSubBeetle1_1evt,hADCPedSubBeetle1_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedSubBeetle2_1evt,hADCPedSubBeetle2_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedSubBeetle1NoSig_1evt,hADCPedSubBeetle1NoSig_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedSubBeetle2NoSig_1evt,hADCPedSubBeetle2NoSig_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedSubBeetle1NoSig2_1evt,hADCPedSubBeetle1NoSig2_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedSubBeetle2NoSig2_1evt,hADCPedSubBeetle2NoSig2_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedCMNSubBeetle1_1evt,hADCPedCMNSubBeetle1_1evt,"",alg_path_to_figures);
            DrawHist(cADCPedCMNSubBeetle2_1evt,hADCPedCMNSubBeetle2_1evt,"",alg_path_to_figures);
        }

    }

    DrawHist(cCMNBeetle1,hCMNBeetle1,"",alg_path_to_figures);
    DrawHist(cCMNBeetle2,hCMNBeetle2,"",alg_path_to_figures);

    TLegend *legADCCompBeetle1 = CreateLegend2(hADCPedSubBeetle1,"ped sub",hADCPedCMNSubBeetle1,"ped + CMN sub");
    TLegend *legADCCompBeetle2 = CreateLegend2(hADCPedSubBeetle2,"ped sub",hADCPedCMNSubBeetle2,"ped + CMN sub");

    cADCCompBeetle1->SetLogy();
    cADCCompBeetle2->SetLogy();

    DrawHistCompare(cADCCompBeetle1,hADCPedSubBeetle1,hADCPedCMNSubBeetle1,legADCCompBeetle1,alg_path_to_figures);
    DrawHistCompare(cADCCompBeetle2,hADCPedSubBeetle2,hADCPedCMNSubBeetle2,legADCCompBeetle2,alg_path_to_figures);

    TLegend *legRMSBeetle1 = CreateLegend3(hRMSPedSubBeetle1,"ped sub",hRMSPedSubBeetle1NoSig,"after 1st outlier removal",hRMSPedSubBeetle1NoSig2,"after 2nd outlier removal","lpfw",0.42,0.59,0.94,0.89);
    TLegend *legRMSBeetle2 = CreateLegend3(hRMSPedSubBeetle2,"ped sub",hRMSPedSubBeetle2NoSig,"after 1st outlier removal",hRMSPedSubBeetle2NoSig2,"after 2nd outlier removal","lpfw",0.42,0.59,0.94,0.89);

    DrawHistCompare3(cRMSBeetle1,hRMSPedSubBeetle1,hRMSPedSubBeetle1NoSig,hRMSPedSubBeetle1NoSig2,legRMSBeetle1,alg_path_to_figures);
    DrawHistCompare3(cRMSBeetle2,hRMSPedSubBeetle2,hRMSPedSubBeetle2NoSig,hRMSPedSubBeetle2NoSig2,legRMSBeetle2,alg_path_to_figures);

    DrawHist(cchTailBeetle1,hchTailBeetle1,"E",alg_path_to_figures);
    DrawHist(cchTailBeetle2,hchTailBeetle2,"E",alg_path_to_figures);

    DrawHist(cevTailBeetle1,hevTailBeetle1,"E",alg_path_to_figures);
    DrawHist(cevTailBeetle2,hevTailBeetle2,"E",alg_path_to_figures);

    output->cd();

    // Close input ROOT file.
    input->Close();

    return;
}

#endif