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

/*
CMNSubMeanxStrips is one of the CMN subtraction algorithms.

For each group of 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.
Instead of considering all the Beetle channels together, split them in groups of x Beetle channels each:
- number of groups = NBeetle/x;
- first group = {0,...,x-1};
- second group = {x,...,2x-1};
- ...
- in general, i-group = {i*x,...,(i+1)*x-1}.

First part:
- for each event, calculate CMNx_y, with x={1,2} and y={0,...,number of groups} and put this information in the tree;
- for each event, calculate ADCPedCMNSub and put this information in the tree.

Second part:  
- use the information in the tree to create the histograms and graphs and save them.

CMNSubMeanxStrips will create a directory in the output ROOT file called "CMNSubMeanxStrips", where the tree and the histograms, graphs, and canvases will be stored.
*/

// Header guard.
#ifndef __CMNSUBMEANXSTRIPS_C_INCLUDED__
#define __CMNSUBMEANXSTRIPS_C_INCLUDED__


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

void CMNSubMeanxStrips(float x, float factor, int sample, string input_ROOT, TFile *output, TString path_to_figures)
{
  if (x==0)
  {
    cout << "Error! The number of Beetle channels per group is zero. The CMNSubMeanxStrips algorithm will not be executed." << endl;
    return;
  }
  if (NBeetle%(int)x!=0)
  {
    cout << "Error! The number of Beetle channels per group is not a divisor of the total number of Beetle channels. The CMNSubMeanxStrips algorithm will not be executed." << endl;
    return;
  } 
   
  cout << "**************************************************" << endl;
  cout << "Computing CMN according to the x strips algorithm..." << endl;
  cout << "**************************************************" << endl;

  cout << "**************************************************" << endl;
  cout << "Configuration" << endl;
  cout << "x: " << x << 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();

  char x_char[10];
  sprintf(x_char,"%d",(int)x);
  // cout << x_char << endl;
  TString dirCMNSubMeanxStrips = "CMNSubMean"+TString(x_char)+"Strips";
  // cout << dirCMNSubMeanxStrips << endl;
  output->cd();
  output->mkdir(dirCMNSubMeanxStrips);
  output->cd(dirCMNSubMeanxStrips);
  
  TString alg_path_to_figures = path_to_figures+dirCMNSubMeanxStrips;
  TString alg_path_to_make_figures = "mkdir -p "+alg_path_to_figures;
  system(alg_path_to_make_figures.Data());

  int groups = NBeetle/x;
  cout << "The CMN will be determined on " << groups << " groups of Beetle channels." << endl;

  // Copy the EventInfo tree from the input to the output file and add some more information.
  float CMN1[groups];
  float CMN2[groups];
  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    CMN1[iGroup] = 0.;
    CMN2[iGroup] = 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);
  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    EventInfoProcessed->Branch(Form("CMN1_%d",iGroup),&CMN1[iGroup]);
    EventInfoProcessed->Branch(Form("CMN2_%d",iGroup),&CMN2[iGroup]);
  }
  EventInfoProcessed->Branch("ADCPedCMNSub",&ADCPedCMNSub);

  // Event-by-event histograms.

  // Event-by-event distribution of the pedestal subtracted ADC counts (hADCPedSubBeetlex_1evt[groups]).
  TH1F *hADCPedSubBeetle1_1evt[groups];
  TCanvas *cADCPedSubBeetle1_1evt[groups];

  TH1F *hADCPedSubBeetle2_1evt[groups];
  TCanvas *cADCPedSubBeetle2_1evt[groups];

  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    hADCPedSubBeetle1_1evt[iGroup] = new TH1F(Form("hADCPedSubBeetle1Group%d_1evt",iGroup),"",400,-200,200); // Modify it.
    cADCPedSubBeetle1_1evt[iGroup] = new TCanvas(Form("cADCPedSubBeetle1Group%d_1evt",iGroup),"",400,300);
    InitHist(hADCPedSubBeetle1_1evt[iGroup],Form("Pedestal subtracted ADC counts on group %d of Beetle 1",iGroup),"ADC counts","");

    hADCPedSubBeetle2_1evt[iGroup] = new TH1F(Form("hADCPedSubBeetle2Group%d_1evt",iGroup),"",400,-200,200); //Modify it.
    cADCPedSubBeetle2_1evt[iGroup] = new TCanvas(Form("cADCPedSubBeetle2Group%d_1evt",iGroup),"",400,300);
    InitHist(hADCPedSubBeetle2_1evt[iGroup],Form("Pedestal subtracted ADC counts on group %d of Beetle 2",iGroup),"ADC counts","");
  }

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

  // First iteration.
  TH1F *hADCPedSubBeetle1NoSig_1evt[groups];
  TCanvas *cADCPedSubBeetle1NoSig_1evt[groups];

  TH1F *hADCPedSubBeetle2NoSig_1evt[groups];
  TCanvas *cADCPedSubBeetle2NoSig_1evt[groups];

  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    hADCPedSubBeetle1NoSig_1evt[iGroup] = new TH1F(Form("hADCPedSubBeetle1Group%dNoSig_1evt",iGroup),"",400,-200,200); // Modify it.
    cADCPedSubBeetle1NoSig_1evt[iGroup] = new TCanvas(Form("cADCPedSubBeetle1Group%dNoSig_1evt",iGroup),"",400,300);
    InitHist(hADCPedSubBeetle1NoSig_1evt[iGroup],Form("Pedestal subtracted ADC counts on group %d of Beetle 1 - excluding signal",iGroup),"ADC counts","");

    hADCPedSubBeetle2NoSig_1evt[iGroup] = new TH1F(Form("hADCPedSubBeetle2Group%dNoSig_1evt",iGroup),"",400,-200,200); //Modify it.
    cADCPedSubBeetle2NoSig_1evt[iGroup] = new TCanvas(Form("cADCPedSubBeetle2Group%dNoSig_1evt",iGroup),"",400,300);
    InitHist(hADCPedSubBeetle2NoSig_1evt[iGroup],Form("Pedestal subtracted ADC counts on group %d of Beetle 2 - excluding signal",iGroup),"ADC counts","");
  }

  // Second iteration.
  TH1F *hADCPedSubBeetle1NoSig2_1evt[groups];
  TCanvas *cADCPedSubBeetle1NoSig2_1evt[groups];

  TH1F *hADCPedSubBeetle2NoSig2_1evt[groups];
  TCanvas *cADCPedSubBeetle2NoSig2_1evt[groups];

  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    hADCPedSubBeetle1NoSig2_1evt[iGroup] = new TH1F(Form("hADCPedSubBeetle1Group%dNoSig2_1evt",iGroup),"",400,-200,200); // Modify it.
    cADCPedSubBeetle1NoSig2_1evt[iGroup] = new TCanvas(Form("cADCPedSubBeetle1Group%dNoSig2_1evt",iGroup),"",400,300);
    InitHist(hADCPedSubBeetle1NoSig2_1evt[iGroup],Form("Pedestal subtracted ADC counts on group %d of Beetle 1 - excluding signal",iGroup),"ADC counts","");

    hADCPedSubBeetle2NoSig2_1evt[iGroup] = new TH1F(Form("hADCPedSubBeetle2Group%dNoSig2_1evt",iGroup),"",400,-200,200); //Modify it.
    cADCPedSubBeetle2NoSig2_1evt[iGroup] = new TCanvas(Form("cADCPedSubBeetle2Group%dNoSig2_1evt",iGroup),"",400,300);
    InitHist(hADCPedSubBeetle2NoSig2_1evt[iGroup],Form("Pedestal subtracted ADC counts on group %d of Beetle 2 - excluding signal",iGroup),"ADC counts","");
  }

  // Event-by-event distribution of the pedestal and CMN subtracted ADC counts (hADCPedCMNSubBeetlex_1evt[groups]).
  TH1F *hADCPedCMNSubBeetle1_1evt[groups];
  TCanvas *cADCPedCMNSubBeetle1_1evt[groups];

  TH1F *hADCPedCMNSubBeetle2_1evt[groups];
  TCanvas *cADCPedCMNSubBeetle2_1evt[groups];

  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    hADCPedCMNSubBeetle1_1evt[iGroup] = new TH1F(Form("hADCPedCMNSubBeetle1Group%d_1evt",iGroup),"",400,-200,200); // Modify it.
    cADCPedCMNSubBeetle1_1evt[iGroup] = new TCanvas(Form("cADCPedCMNSubBeetle1Group%d_1evt",iGroup),"",400,300);
    InitHist(hADCPedCMNSubBeetle1_1evt[iGroup],Form("Pedestal and CMN subtracted ADC counts on group %d of Beetle 1",iGroup),"ADC counts","");

    hADCPedCMNSubBeetle2_1evt[iGroup] = new TH1F(Form("hADCPedCMNSubBeetle2Group%d_1evt",iGroup),"",400,-200,200); //Modify it.
    cADCPedCMNSubBeetle2_1evt[iGroup] = new TCanvas(Form("cADCPedCMNSubBeetle2Group%d_1evt",iGroup),"",400,300);
    InitHist(hADCPedCMNSubBeetle2_1evt[iGroup],Form("Pedestal and CMN subtracted ADC counts on group %d of Beetle 2",iGroup),"ADC counts","");
  }

  // Histograms containing the information on all events.

  // CMN1.
  TH1F *hCMNBeetle1[groups];
  TCanvas *cCMNBeetle1[groups];

  // CMN2.
  TH1F *hCMNBeetle2[groups];
  TCanvas *cCMNBeetle2[groups];

  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    hCMNBeetle1[iGroup] = new TH1F(Form("hCMNBeetle1Group%d",iGroup),"",100,-50,50); // Modify it.
    cCMNBeetle1[iGroup] = new TCanvas(Form("cCMNBeetle1Group%d",iGroup),"",400,300);
    InitHist(hCMNBeetle1[iGroup],Form("CMN on group %d of Beetle 1",iGroup),"ADC counts","");
    
    hCMNBeetle2[iGroup] = new TH1F(Form("hCMNBeetle2Group%d",iGroup),"",100,-50,50); // Modify it.
    cCMNBeetle2[iGroup] = new TCanvas(Form("cCMNBeetle2Group%d",iGroup),"",400,300);
    InitHist(hCMNBeetle2[iGroup],Form("CMN on group %d of Beetle 2",iGroup),"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[groups];
  float rmsBeetle1_temp[groups];

  // Beetle 2.
  float meanBeetle2_temp[groups];
  float rmsBeetle2_temp[groups];

  // Given the Beetle channel, find the appropriate group.
  int group = -1;

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

    // Reset.
    if (iEvent>0)
    {
      ADCPedCMNSub.clear();

      for (int iGroup=0;iGroup<groups;++iGroup)
      {
	hADCPedSubBeetle1_1evt[iGroup]->Reset();
	hADCPedSubBeetle2_1evt[iGroup]->Reset();
	hADCPedSubBeetle1NoSig_1evt[iGroup]->Reset();
	hADCPedSubBeetle2NoSig_1evt[iGroup]->Reset();
	hADCPedSubBeetle1NoSig2_1evt[iGroup]->Reset();
	hADCPedSubBeetle2NoSig2_1evt[iGroup]->Reset();
	hADCPedCMNSubBeetle1_1evt[iGroup]->Reset();
	hADCPedCMNSubBeetle2_1evt[iGroup]->Reset();
      }
    }

    EventInfo->GetEntry(iEvent);

    // Fill histograms for Beetle 1.
    // Beetle channels in i-group = {i*x,...,(i+1)*x-1}.
    for (int iChannel=0;iChannel<NBeetle;++iChannel) {
      if (maskBeetle1[iChannel] == 0)
      {
	hADCPedSubBeetle1->Fill((float)ADCProcessed->at(iChannel));
	group = iChannel/(int)x;
	// cout << iChannel << ", " << group << endl;
	hADCPedSubBeetle1_1evt[group]->Fill((float)ADCProcessed->at(iChannel));
      }
    }

    // Exclude signal.
    for (int iGroup=0;iGroup<groups;++iGroup)
    {
      meanBeetle1_temp[iGroup] = hADCPedSubBeetle1_1evt[iGroup]->GetMean();
      rmsBeetle1_temp[iGroup] = hADCPedSubBeetle1_1evt[iGroup]->GetRMS();
      hRMSPedSubBeetle1->Fill(rmsBeetle1_temp[iGroup]);
    }

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

    // Exclude signal.
    for (int iGroup=0;iGroup<groups;++iGroup)
    {
      meanBeetle1_temp[iGroup] = hADCPedSubBeetle1NoSig_1evt[iGroup]->GetMean();
      rmsBeetle1_temp[iGroup] = hADCPedSubBeetle1NoSig_1evt[iGroup]->GetRMS();
      hRMSPedSubBeetle1NoSig->Fill(rmsBeetle1_temp[iGroup]);
    }

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

    // Calculate CMN1.
    for (int iGroup=0;iGroup<groups;++iGroup)
    {
      CMN1[iGroup] = hADCPedSubBeetle1NoSig2_1evt[iGroup]->GetMean();
      rmsBeetle1_temp[iGroup] = hADCPedSubBeetle1NoSig2_1evt[iGroup]->GetRMS();
      
      hRMSPedSubBeetle1NoSig2->Fill(rmsBeetle1_temp[iGroup]);

      hCMNBeetle1[iGroup]->Fill(CMN1[iGroup]);
    }

    // Fill histograms for Beetle 2.
    // Beetle channels in i-group = {NBeetle+i*x,...,NBeetle+(i+1)*x-1}.
    for (int iChannel=NBeetle;iChannel<N;++iChannel) {
      if (maskBeetle2[iChannel-NBeetle] == 0)
      {
	hADCPedSubBeetle2->Fill((float)ADCProcessed->at(iChannel));
	group = (iChannel-NBeetle)/(int)x;
	hADCPedSubBeetle2_1evt[group]->Fill((float)ADCProcessed->at(iChannel));
      }
    }

    // Exclude signal.
    for (int iGroup=0;iGroup<groups;++iGroup)
    {
      meanBeetle2_temp[iGroup] = hADCPedSubBeetle2_1evt[iGroup]->GetMean();
      rmsBeetle2_temp[iGroup] = hADCPedSubBeetle2_1evt[iGroup]->GetRMS();
      hRMSPedSubBeetle2->Fill(rmsBeetle2_temp[iGroup]);
    }

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

    // Exclude signal.
    for (int iGroup=0;iGroup<groups;++iGroup)
    {
      meanBeetle2_temp[iGroup] = hADCPedSubBeetle2NoSig_1evt[iGroup]->GetMean();
      rmsBeetle2_temp[iGroup] = hADCPedSubBeetle2NoSig_1evt[iGroup]->GetRMS();
      hRMSPedSubBeetle2NoSig->Fill(rmsBeetle2_temp[iGroup]);
    }

    // Fill histograms again for Beetle 2.

    for (int iChannel=NBeetle;iChannel<N;++iChannel) {
      if (maskBeetle2[iChannel-NBeetle] == 0)
      {
	group = (iChannel-NBeetle)/(int)x;
    
    // For testing only!
    /*
    if (iEvent == sample) 
    {
    cout << "Group: " << group << endl;
    cout << "Factor for Beetle 2: " << factor << endl;
    cout << "RMS for Beetle 2: " << rmsBeetle2_temp[group] << endl;
    cout << "Factor*RMS for Beetle 2: " << factor*rmsBeetle2_temp[group] << endl;
    cout << abs((float)ADCProcessed->at(iChannel)-meanBeetle2_temp[group]) << endl;
    }
    */

    if (abs((float)ADCProcessed->at(iChannel)-meanBeetle2_temp[group])<factor*rmsBeetle2_temp[group])
	  hADCPedSubBeetle2NoSig2_1evt[group]->Fill((float)ADCProcessed->at(iChannel));
      }
    }


    // Calculate CMN2.
    for (int iGroup=0;iGroup<groups;++iGroup)
    {
      CMN2[iGroup] = hADCPedSubBeetle2NoSig2_1evt[iGroup]->GetMean();
      rmsBeetle2_temp[iGroup] = hADCPedSubBeetle2NoSig2_1evt[iGroup]->GetRMS();
      
      hRMSPedSubBeetle2NoSig2->Fill(rmsBeetle2_temp[iGroup]);
      
      hCMNBeetle2[iGroup]->Fill(CMN2[iGroup]);
    }

    // Subtract CMN1.
    for (int iChannel=0;iChannel<NBeetle;++iChannel) {
      group = iChannel/(int)x;
      ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN1[group]);
    hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN1[group]);

      if (maskBeetle1[iChannel] == 0)
      {
	hADCPedCMNSubBeetle1->Fill((float)ADCProcessed->at(iChannel)-CMN1[group]);
	hADCPedCMNSubBeetle1_1evt[group]->Fill((float)ADCProcessed->at(iChannel)-CMN1[group]);
	if (((float)ADCProcessed->at(iChannel)-CMN1[group]) > 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) {
      group = (iChannel-NBeetle)/(int)x;
      ADCPedCMNSub.push_back((float)ADCProcessed->at(iChannel)-CMN2[group]);
      hADCPedCMNSub[iChannel]->Fill((float)ADCProcessed->at(iChannel)-CMN2[group]);

      if (maskBeetle2[iChannel-NBeetle] == 0)
      {
	hADCPedCMNSubBeetle2->Fill((float)ADCProcessed->at(iChannel)-CMN2[group]);
	hADCPedCMNSubBeetle2_1evt[group]->Fill((float)ADCProcessed->at(iChannel)-CMN2[group]);
	if (((float)ADCProcessed->at(iChannel)-CMN2[group]) > 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) {
      for (int iGroup=0;iGroup<groups;++iGroup)
      {
	DrawHist(cADCPedSubBeetle1_1evt[iGroup],hADCPedSubBeetle1_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedSubBeetle2_1evt[iGroup],hADCPedSubBeetle2_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedSubBeetle1NoSig_1evt[iGroup],hADCPedSubBeetle1NoSig_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedSubBeetle2NoSig_1evt[iGroup],hADCPedSubBeetle2NoSig_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedSubBeetle1NoSig2_1evt[iGroup],hADCPedSubBeetle1NoSig2_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedSubBeetle2NoSig2_1evt[iGroup],hADCPedSubBeetle2NoSig2_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedCMNSubBeetle1_1evt[iGroup],hADCPedCMNSubBeetle1_1evt[iGroup],"",alg_path_to_figures);
	DrawHist(cADCPedCMNSubBeetle2_1evt[iGroup],hADCPedCMNSubBeetle2_1evt[iGroup],"",alg_path_to_figures);
      }
    }

  }

  for (int iGroup=0;iGroup<groups;++iGroup)
  {
    DrawHist(cCMNBeetle1[iGroup],hCMNBeetle1[iGroup],"",alg_path_to_figures);
    DrawHist(cCMNBeetle2[iGroup],hCMNBeetle2[iGroup],"",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