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

/*
   CMNSubMean is one of the CMN subtraction algorithms.

   Pedestal subtracted ADC counts ---> Exclude signal (requiring |ADC|<th) ---> 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), ignores strips that may potentially contain a signal (by setting a threshold) and creates the event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig_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.

   The configuration of the CMN algorithm is defined by th1 and th2.
   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.

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

// Header guard.
#ifndef __CMNSubMean_C_INCLUDED__
#define __CMNSubMean_C_INCLUDED__

void CMNSubMean(float th1, float th2, int sample, string input_ROOT, TFile *output, TString path_to_figures);

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

    cout << "**************************************************" << endl;
    cout << "Configuration" << endl;
    cout << "th1: " << th1 << endl;
    cout << "th2: " << th2 << 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 dirCMNSubMean = "CMNSubMean";
    output->cd();
    output->mkdir(dirCMNSubMean);
    output->cd(dirCMNSubMean);

    TString alg_path_to_figures = path_to_figures+dirCMNSubMean;
    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 (by setting a threshold).
    // Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexNoSig_1evt).
    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","");

    // 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","",2*th1,-th1,th1); // Modify it.
    TCanvas *cCMNBeetle1 = new TCanvas("cCMNBeetle1","",400,300);
    InitHist(hCMNBeetle1,"CMN on Beetle 1","ADC counts","");

    // CMN2.
    TH1F *hCMNBeetle2 = new TH1F("hCMNBeetle2","",2*th2,-th2,th2); // 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);

    // 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","");
    }
    
    // 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();
            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));
                if (abs((float)ADCProcessed->at(iChannel))<th1)
                    hADCPedSubBeetle1NoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Calculate CMN1.
        CMN1 = hADCPedSubBeetle1NoSig_1evt->GetMean();
        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));
                if (abs((float)ADCProcessed->at(iChannel))<th2)
                    hADCPedSubBeetle2NoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
            }
        }

        // Calculate CMN2.
        CMN2 = hADCPedSubBeetle2NoSig_1evt->GetMean();
        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);
            }
        }

        // 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);
            }
        }

        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(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);  

    output->cd();

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

    return;
}

#endif