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

/*
CMNSubMeanOddEven is one of the CMN subtraction algorithms.

For odd and for even Beetle channels, separately.
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 (hADCPedSubBeetlexOdd_1evt and hADCPedSubBeetlexEven_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 (hADCPedSubBeetlexOddNoSig_1evt and hADCPedSubBeetlexEvenNoSig_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 (hADCPedSubBeetlexOddNoSig2_1evt and hADCPedSubBeetlexEvenNoSig2_1evt), calculates the average, subtracts it to every strip, and creates the event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlexOdd_1evt and hADCPedCMNSubBeetlexEven_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 (hCMNBeetlexOdd and hCMNBeetlexEven).

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.

CMNSubMeanOddEven will create a directory in the output ROOT file called "CMNSubMeanOddEven", where the tree (with CMN1Odd, CMN1Even, CMN2Odd, CMN2Even, and ADCPedCMNSub) and the histograms, graphs, and canvases will be stored.
*/

// Header guard.
#ifndef __CMNSUBMEANODDEVEN_C_INCLUDED__
#define __CMNSUBMEANODDEVEN_C_INCLUDED__

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

void CMNSubMeanOddEven(float factor, int sample, string input_ROOT, TFile *output, TString path_to_figures)
{
  cout << "**************************************************" << endl;
  cout << "Computing CMN according to the mean odd/even 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 dirCMNSubMeanOddEven = "CMNSubMeanOddEven";
  output->cd();
  output->mkdir(dirCMNSubMeanOddEven);
  output->cd(dirCMNSubMeanOddEven);

  TString alg_path_to_figures = path_to_figures+dirCMNSubMeanOddEven;
  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 CMN1Odd = 0.;
  float CMN1Even = 0.;
  float CMN2Odd = 0.;
  float CMN2Even = 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("CMN1Odd",&CMN1Odd);
  EventInfoProcessed->Branch("CMN1Even",&CMN1Even);
  EventInfoProcessed->Branch("CMN2Odd",&CMN2Odd);
  EventInfoProcessed->Branch("CMN2Even",&CMN2Even);
  EventInfoProcessed->Branch("ADCPedCMNSub",&ADCPedCMNSub);

  // Event-by-event histograms.

  // Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlexOdd_1evt and hADCPedSubBeetlexEven_1evt).

  // Odd Beetle channels.
  TH1F *hADCPedSubBeetle1Odd_1evt = new TH1F("hADCPedSubBeetle1Odd_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedSubBeetle1Odd_1evt = new TCanvas("cADCPedSubBeetle1Odd_1evt","",400,300);
  InitHist(hADCPedSubBeetle1Odd_1evt,"Pedestal subtracted ADC counts on odd channels of Beetle 1","ADC counts","");

  TH1F *hADCPedSubBeetle2Odd_1evt = new TH1F("hADCPedSubBeetle2Odd_1evt","",400,-200,200); //Modify it.
  TCanvas *cADCPedSubBeetle2Odd_1evt = new TCanvas("cADCPedSubBeetle2Odd_1evt","",400,300);
  InitHist(hADCPedSubBeetle2Odd_1evt,"Pedestal subtracted ADC counts on odd channels of Beetle 2","ADC counts","");

  // Even Beetle channels.
  TH1F *hADCPedSubBeetle1Even_1evt = new TH1F("hADCPedSubBeetle1Even_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedSubBeetle1Even_1evt = new TCanvas("cADCPedSubBeetle1Even_1evt","",400,300);
  InitHist(hADCPedSubBeetle1Even_1evt,"Pedestal subtracted ADC counts on even channels of Beetle 1","ADC counts","");

  TH1F *hADCPedSubBeetle2Even_1evt = new TH1F("hADCPedSubBeetle2Even_1evt","",400,-200,200); //Modify it.
  TCanvas *cADCPedSubBeetle2Even_1evt = new TCanvas("cADCPedSubBeetle2Even_1evt","",400,300);
  InitHist(hADCPedSubBeetle2Even_1evt,"Pedestal subtracted ADC counts on even channels of Beetle 2","ADC counts","");

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

  // First iteration.

  // Odd Beetle channels.
  TH1F *hADCPedSubBeetle1OddNoSig_1evt = new TH1F("hADCPedSubBeetle1OddNoSig_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedSubBeetle1OddNoSig_1evt = new TCanvas("cADCPedSubBeetle1OddNoSig_1evt","",400,300);
  InitHist(hADCPedSubBeetle1OddNoSig_1evt,"Pedestal subtracted ADC counts on odd channels of Beetle 1 - excluding signal","ADC counts","");

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

  // Even Beetle channels.
  TH1F *hADCPedSubBeetle1EvenNoSig_1evt = new TH1F("hADCPedSubBeetle1EvenNoSig_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedSubBeetle1EvenNoSig_1evt = new TCanvas("cADCPedSubBeetle1EvenNoSig_1evt","",400,300);
  InitHist(hADCPedSubBeetle1EvenNoSig_1evt,"Pedestal subtracted ADC counts on even channels of Beetle 1 - excluding signal","ADC counts","");

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

  // Second iteration.

  // Odd Beetle channels.
  TH1F *hADCPedSubBeetle1OddNoSig2_1evt = new TH1F("hADCPedSubBeetle1OddNoSig2_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedSubBeetle1OddNoSig2_1evt = new TCanvas("cADCPedSubBeetle1OddNoSig2_1evt","",400,300);
  InitHist(hADCPedSubBeetle1OddNoSig2_1evt,"Pedestal subtracted ADC counts on odd channels of Beetle 1 - excluding signal","ADC counts","");

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

  // Even Beetle channels.
  TH1F *hADCPedSubBeetle1EvenNoSig2_1evt = new TH1F("hADCPedSubBeetle1EvenNoSig2_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedSubBeetle1EvenNoSig2_1evt = new TCanvas("cADCPedSubBeetle1EvenNoSig2_1evt","",400,300);
  InitHist(hADCPedSubBeetle1EvenNoSig2_1evt,"Pedestal subtracted ADC counts on even channels of Beetle 1 - excluding signal","ADC counts","");

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

  // Event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlexOdd_1evt and hADCPedCMNSubBeetlexEven_1evt).

  // Odd Beetle channels.
  TH1F *hADCPedCMNSubBeetle1Odd_1evt = new TH1F("hADCPedCMNSubBeetle1Odd_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedCMNSubBeetle1Odd_1evt = new TCanvas("cADCPedCMNSubBeetle1Odd_1evt","",400,300);
  InitHist(hADCPedCMNSubBeetle1Odd_1evt,"Pedestal and CMN subtracted ADC counts on odd channels of Beetle 1","ADC counts","");

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

  // Even Beetle channels.
  TH1F *hADCPedCMNSubBeetle1Even_1evt = new TH1F("hADCPedCMNSubBeetle1Even_1evt","",400,-200,200); // Modify it.
  TCanvas *cADCPedCMNSubBeetle1Even_1evt = new TCanvas("cADCPedCMNSubBeetle1Even_1evt","",400,300);
  InitHist(hADCPedCMNSubBeetle1Even_1evt,"Pedestal and CMN subtracted ADC counts on even channels of Beetle 1","ADC counts","");

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

  // Histograms containing the information on all events.

  // CMN1.

  // Odd Beetle channels.
  TH1F *hCMNBeetle1Odd = new TH1F("hCMNBeetle1Odd","",100,-50,50); // Modify it.
  TCanvas *cCMNBeetle1Odd = new TCanvas("cCMNBeetle1Odd","",400,300);
  InitHist(hCMNBeetle1Odd,"CMN on odd channels of Beetle 1","ADC counts","");

  // Even Beetle channels.
  TH1F *hCMNBeetle1Even = new TH1F("hCMNBeetle1Even","",100,-50,50); // Modify it.
  TCanvas *cCMNBeetle1Even = new TCanvas("cCMNBeetle1Even","",400,300);
  InitHist(hCMNBeetle1Even,"CMN on even channels of Beetle 1","ADC counts","");

  // CMN2.

  // Odd Beetle channels.
  TH1F *hCMNBeetle2Odd = new TH1F("hCMNBeetle2Odd","",100,-50,50); // Modify it.
  TCanvas *cCMNBeetle2Odd = new TCanvas("cCMNBeetle2Odd","",400,300);
  InitHist(hCMNBeetle2Odd,"CMN on odd channels of Beetle 2","ADC counts","");

  // Even Beetle channels.
  TH1F *hCMNBeetle2Even = new TH1F("hCMNBeetle2Even","",100,-50,50); // Modify it.
  TCanvas *cCMNBeetle2Even = new TCanvas("cCMNBeetle2Even","",400,300);
  InitHist(hCMNBeetle2Even,"CMN on even channels of 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 meanBeetle1Odd_temp;
  float meanBeetle1Even_temp;
  float rmsBeetle1Odd_temp;
  float rmsBeetle1Even_temp;

  // Beetle2.

  float meanBeetle2Odd_temp;
  float meanBeetle2Even_temp;
  float rmsBeetle2Odd_temp;
  float rmsBeetle2Even_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();

      hADCPedSubBeetle1Odd_1evt->Reset();
      hADCPedSubBeetle2Odd_1evt->Reset();
      hADCPedSubBeetle1OddNoSig_1evt->Reset();
      hADCPedSubBeetle2OddNoSig_1evt->Reset();
      hADCPedSubBeetle1OddNoSig2_1evt->Reset();
      hADCPedSubBeetle2OddNoSig2_1evt->Reset();
      hADCPedCMNSubBeetle1Odd_1evt->Reset();
      hADCPedCMNSubBeetle2Odd_1evt->Reset();

      hADCPedSubBeetle1Even_1evt->Reset();
      hADCPedSubBeetle2Even_1evt->Reset();
      hADCPedSubBeetle1EvenNoSig_1evt->Reset();
      hADCPedSubBeetle2EvenNoSig_1evt->Reset();
      hADCPedSubBeetle1EvenNoSig2_1evt->Reset();
      hADCPedSubBeetle2EvenNoSig2_1evt->Reset();
      hADCPedCMNSubBeetle1Even_1evt->Reset();
      hADCPedCMNSubBeetle2Even_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));
	if (iChannel%2!=0)
	  hADCPedSubBeetle1Odd_1evt->Fill((float)ADCProcessed->at(iChannel));
	else
	  hADCPedSubBeetle1Even_1evt->Fill((float)ADCProcessed->at(iChannel));
      }
    }

    // Exclude signal.
    meanBeetle1Odd_temp = hADCPedSubBeetle1Odd_1evt->GetMean();
    rmsBeetle1Odd_temp = hADCPedSubBeetle1Odd_1evt->GetRMS();

    meanBeetle1Even_temp = hADCPedSubBeetle1Even_1evt->GetMean();
    rmsBeetle1Even_temp = hADCPedSubBeetle1Even_1evt->GetRMS();

    hRMSPedSubBeetle1->Fill(rmsBeetle1Odd_temp);
    hRMSPedSubBeetle1->Fill(rmsBeetle1Even_temp);

    // Fill histograms again for Beetle 1.
    for (int iChannel=0;iChannel<NBeetle;++iChannel) {
      if (maskBeetle1[iChannel] == 0)
      {
	if (iChannel%2!=0)
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1Odd_temp)<factor*rmsBeetle1Odd_temp)
	    hADCPedSubBeetle1OddNoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
	else
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1Even_temp)<factor*rmsBeetle1Even_temp)
	    hADCPedSubBeetle1EvenNoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
      }
    }

    // Exclude signal.
    meanBeetle1Odd_temp = hADCPedSubBeetle1OddNoSig_1evt->GetMean();
    rmsBeetle1Odd_temp = hADCPedSubBeetle1OddNoSig_1evt->GetRMS();

    meanBeetle1Even_temp = hADCPedSubBeetle1EvenNoSig_1evt->GetMean();
    rmsBeetle1Even_temp = hADCPedSubBeetle1EvenNoSig_1evt->GetRMS();

    hRMSPedSubBeetle1NoSig->Fill(rmsBeetle1Odd_temp);
    hRMSPedSubBeetle1NoSig->Fill(rmsBeetle1Even_temp);

    // Fill histograms again for Beetle 1.
    for (int iChannel=0;iChannel<NBeetle;++iChannel) {
      if (maskBeetle1[iChannel] == 0)
      {
	if (iChannel%2!=0)
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1Odd_temp)<factor*rmsBeetle1Odd_temp)
	    hADCPedSubBeetle1OddNoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
	else
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle1Even_temp)<factor*rmsBeetle1Even_temp)
	    hADCPedSubBeetle1EvenNoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
      }
    }

    // Calculate CMN1.

    // Odd Beetle channels.
    CMN1Odd = hADCPedSubBeetle1OddNoSig2_1evt->GetMean();
    rmsBeetle1Odd_temp = hADCPedSubBeetle1OddNoSig2_1evt->GetRMS();

    hRMSPedSubBeetle1NoSig2->Fill(rmsBeetle1Odd_temp);

    hCMNBeetle1Odd->Fill(CMN1Odd);

    // Even Beetle channels.
    CMN1Even = hADCPedSubBeetle1EvenNoSig2_1evt->GetMean();
    rmsBeetle1Even_temp = hADCPedSubBeetle1EvenNoSig2_1evt->GetRMS();

    hRMSPedSubBeetle1NoSig2->Fill(rmsBeetle1Even_temp);

    hCMNBeetle1Even->Fill(CMN1Even);

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

    // Exclude signal.
    meanBeetle2Odd_temp = hADCPedSubBeetle2Odd_1evt->GetMean();
    rmsBeetle2Odd_temp = hADCPedSubBeetle2Odd_1evt->GetRMS();

    meanBeetle2Even_temp = hADCPedSubBeetle2Even_1evt->GetMean();
    rmsBeetle2Even_temp = hADCPedSubBeetle2Even_1evt->GetRMS();

    hRMSPedSubBeetle2->Fill(rmsBeetle2Odd_temp);
    hRMSPedSubBeetle2->Fill(rmsBeetle2Even_temp);

    // Fill histograms again for Beetle 2.
    for (int iChannel=NBeetle;iChannel<N;++iChannel) {
      if (maskBeetle2[iChannel-NBeetle] == 0)
      {
	if (iChannel%2!=0)
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2Odd_temp)<factor*rmsBeetle2Odd_temp)
	    hADCPedSubBeetle2OddNoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
	else
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2Even_temp)<factor*rmsBeetle2Even_temp)
	    hADCPedSubBeetle2EvenNoSig_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
      }
    }

    // Exclude signal.
    meanBeetle2Odd_temp = hADCPedSubBeetle2OddNoSig_1evt->GetMean();
    rmsBeetle2Odd_temp = hADCPedSubBeetle2OddNoSig_1evt->GetRMS();

    meanBeetle2Even_temp = hADCPedSubBeetle2EvenNoSig_1evt->GetMean();
    rmsBeetle2Even_temp = hADCPedSubBeetle2EvenNoSig_1evt->GetRMS();

    hRMSPedSubBeetle2NoSig->Fill(rmsBeetle2Odd_temp);
    hRMSPedSubBeetle2NoSig->Fill(rmsBeetle2Even_temp);

    // Fill histograms again for Beetle 2.
    for (int iChannel=NBeetle;iChannel<N;++iChannel) {
      if (maskBeetle2[iChannel-NBeetle] == 0)
      {
	if (iChannel%2!=0)
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2Odd_temp)<factor*rmsBeetle2Odd_temp)
	    hADCPedSubBeetle2OddNoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
	else
	{
	  if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2Even_temp)<factor*rmsBeetle2Even_temp)
	    hADCPedSubBeetle2EvenNoSig2_1evt->Fill((float)ADCProcessed->at(iChannel));
	}
      }
    }

    // Calculate CMN2.

    // Odd Beetle channels.
    CMN2Odd = hADCPedSubBeetle2OddNoSig2_1evt->GetMean();
    rmsBeetle2Odd_temp = hADCPedSubBeetle2OddNoSig2_1evt->GetRMS();

    hRMSPedSubBeetle2NoSig2->Fill(rmsBeetle2Odd_temp);

    hCMNBeetle2Odd->Fill(CMN2Odd);

    // Even Beetle channels.
    CMN2Even = hADCPedSubBeetle2EvenNoSig2_1evt->GetMean();
    rmsBeetle2Even_temp = hADCPedSubBeetle2EvenNoSig2_1evt->GetRMS();

    hRMSPedSubBeetle2NoSig2->Fill(rmsBeetle2Even_temp);

    hCMNBeetle2Even->Fill(CMN2Even);

    // Subtract CMN1.
    for (int iChannel=0;iChannel<NBeetle;++iChannel) {
      if (iChannel%2!=0)
      {
          ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN1Odd);
          hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN1Odd);
      }
      else
      {
          ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN1Even);
          hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN1Even);
      }
      if (maskBeetle1[iChannel] == 0)
      {
	if (iChannel%2!=0)
	{
	  hADCPedCMNSubBeetle1->Fill((float)ADCProcessed->at(iChannel)-CMN1Odd);
	  hADCPedCMNSubBeetle1Odd_1evt->Fill((float)ADCProcessed->at(iChannel)-CMN1Odd);
	}
	else
	{
	  hADCPedCMNSubBeetle1->Fill((float)ADCProcessed->at(iChannel)-CMN1Even);
	  hADCPedCMNSubBeetle1Even_1evt->Fill((float)ADCProcessed->at(iChannel)-CMN1Even);
	}
	if (((((float)ADCProcessed->at(iChannel)-CMN1Odd) > 30.) && (iChannel%2!=0)) || ((((float)ADCProcessed->at(iChannel)-CMN1Even) > 30.) && (iChannel%2==0)))
	{
	  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) {
      if (iChannel%2!=0)
      {
          ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN2Odd);
          hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN2Odd);
      }
          else
          {
	ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN2Even);
          hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN2Even);
          }
          if (maskBeetle2[iChannel-NBeetle] == 0)
      {
	if (iChannel%2!=0)
	{
	  hADCPedCMNSubBeetle2->Fill((float)ADCProcessed->at(iChannel)-CMN2Odd);
	  hADCPedCMNSubBeetle2Odd_1evt->Fill((float)ADCProcessed->at(iChannel)-CMN2Odd);
	}
	else
	{
	  hADCPedCMNSubBeetle2->Fill((float)ADCProcessed->at(iChannel)-CMN2Even);
	  hADCPedCMNSubBeetle2Even_1evt->Fill((float)ADCProcessed->at(iChannel)-CMN2Even);
	}

	if (((((float)ADCProcessed->at(iChannel)-CMN2Odd) > 30.) && (iChannel%2!=0)) || ((((float)ADCProcessed->at(iChannel)-CMN2Even) > 30.) && (iChannel%2==0)))
	{
	  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) {
      // Odd Beetle channels.
      DrawHist(cADCPedSubBeetle1Odd_1evt,hADCPedSubBeetle1Odd_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle2Odd_1evt,hADCPedSubBeetle2Odd_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle1OddNoSig_1evt,hADCPedSubBeetle1OddNoSig_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle2OddNoSig_1evt,hADCPedSubBeetle2OddNoSig_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle1OddNoSig2_1evt,hADCPedSubBeetle1OddNoSig2_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle2OddNoSig2_1evt,hADCPedSubBeetle2OddNoSig2_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedCMNSubBeetle1Odd_1evt,hADCPedCMNSubBeetle1Odd_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedCMNSubBeetle2Odd_1evt,hADCPedCMNSubBeetle2Odd_1evt,"",alg_path_to_figures);

      // Even Beetle channels.
      DrawHist(cADCPedSubBeetle1Even_1evt,hADCPedSubBeetle1Even_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle2Even_1evt,hADCPedSubBeetle2Even_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle1EvenNoSig_1evt,hADCPedSubBeetle1EvenNoSig_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle2EvenNoSig_1evt,hADCPedSubBeetle2EvenNoSig_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle1EvenNoSig2_1evt,hADCPedSubBeetle1EvenNoSig2_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedSubBeetle2EvenNoSig2_1evt,hADCPedSubBeetle2EvenNoSig2_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedCMNSubBeetle1Even_1evt,hADCPedCMNSubBeetle1Even_1evt,"",alg_path_to_figures);
      DrawHist(cADCPedCMNSubBeetle2Even_1evt,hADCPedCMNSubBeetle2Even_1evt,"",alg_path_to_figures);
    }
    
  }

  // Odd Beetle channels.
  DrawHist(cCMNBeetle1Odd,hCMNBeetle1Odd,"",alg_path_to_figures);
  DrawHist(cCMNBeetle2Odd,hCMNBeetle2Odd,"",alg_path_to_figures);

  // Even Beetle channels.
  DrawHist(cCMNBeetle1Even,hCMNBeetle1Even,"",alg_path_to_figures);
  DrawHist(cCMNBeetle2Even,hCMNBeetle2Even,"",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
