// Read in data file produced by the fwdsct_common_mode.cpp 
// macro and plot correlation matrix
// SJMP 18/10/2000
//
// ChangeLog:
// 23 October 2000 SJMP Added more options, not 0,1 but -1,1
//

// 1536 * 1536 matrix : entire module
void fwdsct_common_mode_matrix(void) {
   printf("Usage:\n");
   return;
}

void fwdsct_common_mode_matrix(Char_t * filename) {
   
   Int_t chan_start  = 0;
   Int_t chan_end    = 1535;
   fwdsct_common_mode_matrix(filename,chan_start,chan_end,1000);
   return;
}

// 768 * 768 : one link
void fwdsct_common_mode_matrix(Char_t * filename, Int_t link) {
   
   if ( !(link == 0 || link == 1) ) {
      printf("Illegal link: %d\n",link);
      printf("Exiting..\n");
      return;
   }
   Int_t chan_start  =   0+link*768;
   Int_t chan_end    = 767+link*768;
   fwdsct_common_mode_matrix(filename,chan_start,chan_end,1000);
   return;
}

// 128 * 128 : one chip
void fwdsct_common_mode_matrix(Char_t * filename, Int_t link, Int_t chip) {

   if ( chip < 0 || chip > 5 ) {
      printf("Illegal chip number: %d (possible values: 0-5)\n",chip);
      printf("Exiting..\n");
      return;
   }
   if ( !(link == 0 || link == 1) ) {
      printf("Illegal link: %d (possible values: 0,1)\n",link);
      printf("Exiting..\n");
      return;
   }
   Int_t chan_start  =  chip   *128  +link*768;
   Int_t chan_end    = (chip+1)*128-1+link*768;
   fwdsct_common_mode_matrix(filename,chan_start,chan_end,1000);
   return;
}

void fwdsct_common_mode_matrix(Char_t * filename,
			       Int_t chan_start, Int_t chan_end,
			       Int_t n_events) {

   if ( !(chan_start >=   0  && chan_end >=    0 &&
	  chan_start <= 1535 && chan_end <= 1535 &&
          chan_start <= chan_end              ) ) {
      printf("Illegal channel range %d - %d\n",chan_start,chan_end);
      printf("Exiting...\n");
      return;
   }
   printf("Generating matrix for channel %d - %d\n",chan_start,chan_end);
   int i;
   Char_t s[250];

   // set-up drawing style
   gROOT->SetStyle("Plain");
   gROOT->GetStyle("Plain")->SetPalette(1);
   gROOT->GetStyle("Plain")->SetOptStat(0);
   gROOT->GetStyle("Plain")->SetTitleColor(10);
   gROOT->GetStyle("Plain")->SetTitleW(1);
   gROOT->GetStyle("Plain")->SetTitleBorderSize(0);
   gROOT->GetStyle("Plain")->SetPadBottomMargin((float)0.15);
   gROOT->GetStyle("Plain")->SetPadTopMargin((float)0.15);
   gROOT->GetStyle("Plain")->SetPadLeftMargin((float)0.15);
   gROOT->GetStyle("Plain")->SetPadRightMargin((float)0.15);
   gROOT->GetStyle("Plain")->SetTitleOffset(1.5,"X");
   gROOT->GetStyle("Plain")->SetTitleOffset(1.5,"Y");

   // get data
   sprintf(s,"%s.root",filename);
   TFile * hfile1 = new TFile(s,"READ");
   if ( hfile1 == 0 ) {
      printf("Exiting...\n");
      return 1;
   }
   TTree * data = (TTree*)hfile1->Get("data");
   data->Print();
   if ( data == 0 ) {
      printf("Empty or incorrect data file, no tree data present\n");
      printf("Exiting...\n");
      return 1;
   }
   TBranch * branch = data->GetBranch("event");
   if ( branch == 0 ) {
      printf("Empty or incorrect data file, no branch event present\n");
      printf("Exiting...\n");
      return 1;
   }
   FwdSctEvent * event = new FwdSctEvent();
   branch->SetAddress(&event);
   int n = (int)data->GetEntries();
   if ( n_events > 0 ) n = n_events;
   
   // calcalute occupancies
   printf("Calculating occupancies...\n");
   TH1F * opc   = new TH1F("opc",  "average occupancy"   ,
			   chan_end-chan_start+1,chan_start,chan_end);
   opc->SetXTitle("channel");
   TH2F * opcxy = new TH2F("opcxy","average occupancy xy",
			   chan_end-chan_start+1,chan_start,chan_end,
			   chan_end-chan_start+1,chan_start,chan_end);
   opcxy->SetXTitle("channel");
   opcxy->SetYTitle("channel");
   for ( i = 0 ; i < n ; ++i ) {
      data->GetEvent(i);
      if ( i%1 == 0 ) printf("\tevent %d\n",i);
      for ( int x = chan_start ; x <= chan_end ; ++x ) {
	 char xbit = event->GetBit(x);
	 opc->SetBinContent(x,opc->GetBinContent(x)+(float)xbit/(float)n);
	 for ( int y = chan_start ; y <= chan_end ; ++y ) {
	    char ybit = event->GetBit(y);
	    int x_and_y = 0;
	    if      ( xbit == 1  && ybit == 1 ) x_and_y = +1;
	    else if ( (xbit == -1 || xbit == 1) &&
		      (ybit == -1 || ybit == 1) ) x_and_y = -1;
	    opcxy->SetBinContent(opcxy->GetBin(x,y),
				 opcxy->GetBinContent(opcxy->GetBin(x,y))+
				 (float)x_and_y/(float)n);
	 }
      }
   }

   // open output file
   sprintf(s,"%s_chan%d_%d.root",
	   filename,chan_start,chan_end);
   hfile2 = new TFile(s,"RECREATE");
   hfile2->cd();

   // plot and save
   TCanvas * occupancy_canvas = new TCanvas("occupancy_canvas",
					    "occupancy",1000,500);
   occupancy_canvas->Divide(2,1);
   occupancy_canvas->cd(1);
   opc->SetFillColor(46);
   opc->Draw();
   occupancy_canvas->cd(2);
   opcxy->Draw("colz");
   sprintf(s,"%s_%d_%d_occupancy.ps",
	   filename,chan_start,chan_end);
   occupancy_canvas->SaveAs(s);
   opc->Write();
   opcxy->Write();
   
   // calcalute matrix
   printf("Calculating matrix...\n");
   sprintf(s,"correlation matrix channel %d - %d",chan_start,chan_end);
   TH2F *matrix = new TH2F("matrix",s,
			   chan_end-chan_start+1,chan_start,chan_end,
			   chan_end-chan_start+1,chan_start,chan_end);
   matrix->SetXTitle("channel");
   matrix->SetYTitle("channel");
   for ( int x = chan_start ; x <= chan_end ; ++x ) {
      printf("\tchannel %d\n",x);
      for ( int y = chan_start ; y <= chan_end ; ++y ) {
	 // correlation = 1 for the same channels
	 if ( x == y ) {
	    matrix->SetBinContent(matrix->GetBin(x,y),1);
	 }
	 else {
	    double av_x  = (double)opc->GetBinContent(x);
	    double av_y  = (double)opc->GetBinContent(y);
	    double av_xy = (double)opcxy->GetBinContent(opcxy->GetBin(x,y));
	    // if channel is stuck or dead correlation = 0
	    float  diff  = 0.0001;
	    if ( fabs(av_x)   < diff ||
		 fabs(av_x-1) < diff ||
		 fabs(av_x+1) < diff ) {
	       matrix->SetBinContent(matrix->GetBin(x,y),0);
	    }
	    else if ( fabs(av_y)   < diff ||
		      fabs(av_y-1) < diff ||
		      fabs(av_y+1) < diff ) {
	       matrix->SetBinContent(matrix->GetBin(x,y),0);
	    }
	    else  {
	       float corr     = (av_xy-av_x*av_y);
	       float s_x      = TMath::Sqrt(fabs(av_x)-av_x*av_x);
	       float s_y      = TMath::Sqrt(fabs(av_x)-av_x*av_x);
	       float n_corr   = corr/(s_x*s_y);
	       matrix->SetBinContent(matrix->GetBin(x,y),n_corr);
	    }
	 }	 
      }
   }

   // draw matrices
   printf("Drawing matrix...\n");
   TCanvas * matrix_canvas = new TCanvas ("matrix_canvas",
					  "fwdsct_common_mode",750,750);
   matrix_canvas->SetFillColor(0);
   matrix->SetMaximum(1);
   matrix->SetMinimum(-1);
   matrix->Draw("colz");
   matrix_canvas->Update();
   sprintf(s,"%s_%d_%d_matrix.ps",
	   filename,chan_start,chan_end);
   matrix_canvas->SaveAs(s);
   matrix->Write();
   hfile2->Write();
   //   hfile1->Close();
   //   hfile2->Close();

   return;
}