// 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(Int_t run, Int_t scan) {
   
   Int_t chan_start  = 0;
   Int_t chan_end    = 1535;
   fwdsct_common_mode_matrix(run,scan,chan_start,chan_end,1000);
   return;
}

// 768 * 768 : one link
void fwdsct_common_mode_matrix(Int_t run, Int_t scan, 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(run,scan,chan_start,chan_end,1000);
   return;
}

// 128 * 128 : one chip
void fwdsct_common_mode_matrix(Int_t run, Int_t scan, 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(run,scan,chan_start,chan_end,1000);
   return;
}

void fwdsct_common_mode_matrix(Int_t run, Int_t scan,
			       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, j;
   Char_t s[250];
   float occ_matrix[1536][1536];
   float occ_array[1536];
   float ope, mean_ope, sigma, mean_sigma, rms;
   float ope_chip[12];
   int ope_bin;
   for ( int x = 0 ; x <= 1535 ; ++x ) {
     occ_array[x] = 0. ;
     for ( int y = 0 ; y <= 1535 ; ++y ) {
       occ_matrix[x][y] = 0.;
     }
   }
   for ( i = 0 ; i <=11 ; ++i ) {
     ope_chip[i] = 0. ;
   }

   // 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,"F:\\sctvar\\data\\%d_%d.root",run,scan);
   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;
   }
   FwdSctSingleEvent * event = new FwdSctSingleEvent();
   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");
   TH2F * ope_his   = new TH2F("ope_his",  "occupancy per event"   ,
			   100.,0.,1.0,12,0,12);
   ope_his->SetXTitle("OPE");
   ope_his->SetYTitle("Chips");

   /*   TH1F * ope_his   = new TH1F("ope_his",  "occupancy per event"   ,
			   100.,0.,1.0);
   ope_his->SetXTitle("OPE");
   ope_his->SetYTitle("Events");*/


   mean_ope = 0.0;
   for ( i = 0 ; i < n ; ++i ) {
      ope = 0.0;
      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);
	 ope += (float)xbit;
	 j = (int)(x/128);
	 ope_chip[j] += (float)xbit;
     	 // opc->SetBinContent(x,opc->GetBinContent(x)+(float)xbit/(float)n);
	 occ_array[x] += (float)xbit;
	 for ( int y = x ; 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); */
	    if      ( xbit == 1  && ybit == 1 ) occ_matrix[x][y]++ ;
	    //  else if ( (xbit == 0 && ybit == 1) ||
	    //      (xbit == 1 && ybit == 0) ) occ_matrix[x][y]--;
	 }

      }//channel loop

      ope = ope / (chan_end-chan_start+1);
      for (j = 0; j <= 11; ++j){
	ope_chip[j] = ope_chip[j] /128;
      ope_bin = (int)(ope_chip[j] * 100);

      ope_his->SetBinContent(ope_his->GetBin(ope_bin,j+1),
			     ope_his->GetBinContent(ope_his->GetBin(ope_bin,j+1))+1.);

      }

      mean_ope += ope;

   }//event loop

   /*   mean_ope = mean_ope / n;
   sigma = TMath::Sqrt( (mean_ope * (1. - mean_ope) ) / (chan_end-chan_start+1));
   printf("sigma= %f\n",sigma);
   */ 
 
    // fill histograms
      for ( int x = chan_start ; x <= chan_end ; ++x ) {
	 opc->SetBinContent(x-chan_start+1,occ_array[x]/(float)n);
      for ( int y = x ; y <= chan_end ; ++y ) {
	opcxy->SetBinContent(opcxy->GetBin(x-chan_start+1,y-chan_start+1),occ_matrix[x][y]/(float)n);
       	opcxy->SetBinContent(opcxy->GetBin(y-chan_start+1,x-chan_start+1),occ_matrix[x][y]/(float)n);
      }
      }

   // open output file
   sprintf(s,"F:\\sctvar\\data\\%d_%d_chan%d_%d.root",
	   run,scan,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(3,1);
   occupancy_canvas->cd(1);
   opc->SetFillColor(46);
   opc->Draw();
   occupancy_canvas->cd(2);
   opcxy->Draw("colz");
   occupancy_canvas->cd(3);
   ope_his->Draw("box");
   sprintf(s,"F:\\sctvar\\data\\%d_%d_%d_%d_occupancy.ps",
	   run,scan,chan_start,chan_end);
   occupancy_canvas->SaveAs(s);
   opc->Write();
   opcxy->Write();
   ope_his->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-chan_start+1,y-chan_start+1),1);
	 }
	 else {
	    double av_x  = (double)opc->GetBinContent(x-chan_start+1);
	    double av_y  = (double)opc->GetBinContent(y-chan_start+1);
	    double av_xy = (double)opcxy->GetBinContent(opcxy->GetBin(x-chan_start+1,y-chan_start+1));
	    // 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-chan_start+1,y-chan_start+1),0);
	    }
	    else if ( fabs(av_y)   < diff ||
		      fabs(av_y-1) < diff ||
		      fabs(av_y+1) < diff ) {
	       matrix->SetBinContent(matrix->GetBin(x-chan_start+1,y-chan_start+1),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_y)-av_y*av_y);
	       float n_corr   = corr/(s_x*s_y);
	       matrix->SetBinContent(matrix->GetBin(x-chan_start+1,y-chan_start+1),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,"F:\\sctvar\\data\\%d_%d_%d_%d_matrix.ps",
	   run,scan,chan_start,chan_end);
   matrix_canvas->SaveAs(s);
   matrix->Write();
   hfile2->Write();
   //   hfile1->Close();
   //   hfile2->Close();

   return;
}