//////////////////////////////////////////////////////////////////////////
//
// Function: getPlugGainCorrection
// Purpose:  Apply corrections determined by h. Budd to the PEM and PHA
//           towers to account for the gain changes
//          
//////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <iostream>
#include <string>

void getPlugGainCorrection(int irun, float corrpem[52][48],
                                     float corrpha[52][48])
{

  enum {nrun = 215};
  int runs[nrun] = 
  {
     134898,139524,139624,140024,140282,140428,141403,141481,141583,141644,
     141914,141963,142122,142188,142516,142655,142926,143062,143233,143277,
     143612,143659,143766,143878,143971,143987,144246,144436,144848,144895,
     145020,145090,145165,145238,145282,145362,145403,145440,145631,145657,
     146258,146313,146388,146406,146504,146692,146791,146883,146931,146987,
     147054,147104,147182,147265,147266,147268,147466,147750,147814,147881,
     148121,148253,148276,148486,148574,148661,148697,148786,148792,148959,
     148965,149103,149360,149405,149418,149521,149692,149730,149786,150295,
     150698,150743,151181,151187,151217,151403,151783,151857,151943,151945,
     151988,151989,152107,152173,152230,152253,152286,152342,152429,152475,
     152545,152702,152787,152818,152905,152975,153248,153288,153333,153400,
     153428,153638,153688,153708,153710,153723,153899,153961,153992,154191,
     154228,154410,154493,154666,154685,154718,154746,154784,154927,155104,
     155167,155349,155690,155733,155789,155808,155870,155884,155931,155985,
     156161,156353,156410,156501,156502,158949,159129,159239,159272,159301,
     159719,159742,159764,160241,160277,160370,160456,160498,160620,160674,
     160720,160770,160908,160953,161001,161060,162136,162415,162506,162596,
     162911,162950,163055,163077,163416,163482,163535,164033,164212,164323,
     164380,164647,164974,165169,165237,165449,165494,165737,165778,165850,
     165884,165978,166081,166340,166443,166554,166676,166759,166949,166998,
     167208,167344,167791,167877,168136,168455,168521,168577,168624,168658,
     168717,168798,168834,168905,200000
  };

// Normalization factor to get the Z peak for each PEM quadrant in the 
// same place. The hypothesis is that there is problem with the absolute
// energy calibration of the Co60 calibration sources in each quadrant.
// East PEM quadrant 1 is the reference as it is the only quadrant
// cross checked with a Cs137 source. The scorr1 table is applied to both
// the PEM and PHA, and is only for the 2002 data (run < 156600). The
// scorr2 table is for run > 156600.

  float scorr1[2][4] = {
     { 1.0503, 1.0229, 1.0398, 1.0264 },	// West quadrants
     { 1.0000, 1.0118, 1.0203, 1.0303 }		// East quadrants
  };
  float scorr2[2][4] = {
     { 1.0445, 1.0051, 1.0335, 1.0400 },
     { 1.0109, 1.0053, 1.0229, 1.0568 }
  };
   
  float cpemL[52][48],cpemH[52][48];
  float cphaL[52][48],cphaH[52][48];

  for (int i=0; i<52; i++) {
  for (int j=0; j<48; j++) { 
    corrpem[i][j] = 1.0;
    corrpha[i][j] = 1.0;
    cpemL[i][j] = cpemH[i][j] = 1.0;
    cphaL[i][j] = cphaH[i][j] = 1.0;
  }
  }

  if (irun < runs[0]) { 
    std::cout << "getPlugGainCorrection: " << irun << " "
              << "corr = 1" << std::endl;
    return;
  }

  int irtbl = nrun-2;
  if (irun < runs[nrun-1]) {
    for (int i = 0; i < nrun-1; i++) {
      if ( irun >= runs[i] && irun < runs[i+1] ) {
        irtbl = i;
        break;
      }
    }
  }

  std::string lrun, hrun;
  char crunL[20], crunH[20];
  sprintf(crunL, "%d", runs[irtbl  ]);
  sprintf(crunH, "%d", runs[irtbl+1]);


  lrun  = crunL;
  hrun  = crunH;

  std::string srunL, srunH;
  srunL = "/cdf/data38a/s1/calorimeter/pcgain/plug_cor_"
          + lrun + ".dat";
  srunH = "/cdf/data38a/s1/calorimeter/pcgain/plug_cor_"
          + hrun + ".dat";

  FILE *pcfileL, *pcfileH;

  pcfileL=fopen(srunL.c_str(),"r");
  pcfileH=fopen(srunH.c_str(),"r");

  if ( !pcfileL ) {
    std::cout << "getPlugGainCorrection: file not found - " 
              << srunL << std::endl;
    return;
  }
  if ( !pcfileH ) {
    std::cout << "getPlugGainCorrection: file not found - " 
              << srunH << std::endl;
    return;
  }

  int idet, iside, iwedge, itower;
  float pcorr;

//
// The iphi index for itower=0,1,2,3 goes from 0-47 because it includes the
// itile index. To access these towers in CalData, this 0-47 phi index must
// decoded into a wedge (0-23) and itile. This is ieta=14,15,36,37.
// 

  while ( fscanf(pcfileL, "%i %i %i %i %f",
                 &idet, &iside, &iwedge, &itower, &pcorr) != EOF ) {

    int ieta  = (iside == 0) ? 25 - (10 + itower/2) :
                               26 + (10 + itower/2);
    int itile = itower%2;
    int iphi  = 2*iwedge + itile;
    if (pcorr > 0.01) {
      float fscorr( irun > 156600 ? scorr2[iside][iwedge/6] :
				    scorr1[iside][iwedge/6] );
      if (idet==0) cpemL[ieta][iphi]=fscorr/pcorr;
      if (idet==1) cphaL[ieta][iphi]=fscorr/pcorr;
    }
  }
  fclose( pcfileL );

  while ( fscanf(pcfileH, "%i %i %i %i %f",
                 &idet, &iside, &iwedge, &itower, &pcorr) != EOF ) {

    int ieta  = (iside == 0) ? 25 - (10 + itower/2) :
                               26 + (10 + itower/2);
    int itile = itower%2;
    int iphi  = 2*iwedge + itile;
    if (pcorr > 0.01) {
      float fscorr( irun > 156600 ? scorr2[iside][iwedge/6] :
				    scorr1[iside][iwedge/6] );
      if (idet==0) cpemH[ieta][iphi]=fscorr/pcorr;
      if (idet==1) cphaH[ieta][iphi]=fscorr/pcorr;
    }
  }
  fclose( pcfileH );

  float irunL = runs[irtbl  ];
  float irunH = runs[irtbl+1];
  float frun  = irun;

  float f = (frun - irunL)/(irunH - irunL);

  for (int i = 0; i < 52; i++) {
  for (int j = 0; j < 48; j++) {
    corrpem[i][j] = (1.0 - f)*cpemL[i][j] + f*cpemH[i][j];
    corrpha[i][j] = (1.0 - f)*cphaL[i][j] + f*cphaH[i][j];
  }
  }

  std::cout << "getPlugGainCorrection: " << irun
            << " (" << runs[irtbl] << ":" << runs[irtbl+1]
            << ")" << std::endl;

}