#include "Couples.h"

void compcol(Couples *cou1, Couples *cou2, int p, int col)
{
   Map2D *map1 = cou1->getmap(p);
   Map2D *map2 = cou2->getmap(p);
   TH1D *py1 = map1->ProjectionY("py1",col,col);
   TH1D *py2 = map2->ProjectionY("py2",col,col);
   py2->SetLineColor(kRed);
   if (py1->GetMaximum() > py2->GetMaximum()) {
      py1->Draw();
      py2->Draw("same");
   }
   else {
      py2->Draw();
      py1->Draw("same");
   }
}

void getmean(Couples *cou1, Couples *cou2, int p)
{
   Map2D *map1 = cou1->getmap(p);
   Map2D *map2 = cou2->getmap(p);
   double sum10=0, sum11=0;
   for (int i=1; i<=map1->GetNbinsX(); ++i) {
      for (int j=1; j<=map1->GetNbinsY(); ++j) {
         double h = map1->GetBinContent(i,j);
         if (h != 0) {
            sum11 += h;
            sum10 += 1;
         }
      }
   }
   double sum20=0, sum21=0;
   for (int i=1; i<=map2->GetNbinsX(); ++i) {
      for (int j=1; j<=map2->GetNbinsY(); ++j) {
         double h = map2->GetBinContent(i,j);
         if (h != 0) {
            sum21 += h;
            sum20 += 1;
         }
      }
   }
   printf("map means are %lf,%lf\n", sum11/sum10, sum21/sum20);
}

TH1D &combine(const Map2D *m1, const Map2D *m2)
{
   Map2D *combmean = new Map2D(m1);
   Map2D *combdiff = new Map2D(m1);
   combmean->Add(m2,1);
   combmean->Rescale(1,1,0.5);
   combdiff->Add(m2,-1);
   combmean->Mask(m1);
   combmean->Mask(m2);
   combdiff->Mask(m1);
   combdiff->Mask(m2);
   Map2D m1m(m1);
   Map2D m2m(m2);
   m1m.Mask(combmean,1);
   m2m.Mask(combmean,1);
   combmean->Fill(&m1m);
   combmean->Fill(&m2m);
   combmean->SetName("combmean");
   combdiff->SetName("combdiff");
   combmean->SetContour(100);
   combdiff->SetContour(100);
   return combdiff->Profile("combprof","",800,-1600,1600);
}