{
#include <fstream.h>
gROOT->LoadMacro("AtlasUtils.C");
ATLAS->SetPadTickX(1);
ATLAS->SetPadRightMargin(0.12);
ATLAS->SetPadLeftMargin(0.12);
gStyle->SetOptStat(0); 
c1 = new TCanvas("c1","A Simple Graph",500,10,700,500);
c1->Divide(2,3);
c1->cd(1);
//
double Eg = 1.12, EI=1.09, kB=8.617343E-5, TR=293.15, T0=273.15;
double DT_barrel = 3.6, DT_ECouter=15.5, DT_ECmiddle=15.5, DT_ECinner=13.1;  
//from Paul (2011.7.14) For the endcap outer and middle modules I use 15.5 degrees C 
//(silicon cooler than hybrid temp) and for the inners I use 13.1 degrees C.
 
double Temp[3][10][13][56];  //Temp[bec][layer][eta][phi]
double Temp1[3][10][13][56];  //Temp[bec][layer][eta][phi]
double HV[3][10][13][56];  //HV[bec][layer][eta][phi]
double HVLC[3][10][13][56];
double HVLC_norm[3][10][13][56];  
double Phi[56];
double HVcut=140., TEMPcut_high=14.0, TEMPcut_low= -4.;
int barrel[10]={0,0,0,32,40,48,56,0,0,0};
for (int i=0;i<56;i++) Phi[i]=i;
//
for(int i=0;i<3;i++){
  for(int j=0;j<10;j++){
    for(int k=0;k<13;k++){
      for(int l=0;l<56;l++){ 
         HVLC[i][j][k][l] = -999.; 
         HVLC_norm[i][j][k][l] = -999.; 
} } } }
//---------------------------------
int bec,layer,phi,eta,date,time;
//------------------read in TM0---------------------------
int bec0=9,layer0=9,phi0=9,eta0=9,date0=9,time0=9;
double data, data0=-999.;
int idata=0;
double sum=0.;
ifstream fin("sct_TM0_000016.txt");
while(fin>>bec>>layer>>phi>>eta>>data>>date>>time)
{
  if(bec!=bec0 || layer!=layer0 || phi!=phi0 || eta!=eta0){
    if(idata != 0) {
      if(bec0==1) eta0 += 6;
      Temp[bec0][layer0][eta0][phi0] = sum/idata;
//    if(bec0==1 && layer0 ==3 && eta0 == 0) 
//    cout<<bec0<<" "<<layer0<<" "<<eta0-6<<" "<<phi0<<" " 
//     <<Temp[bec0][layer0][eta0][phi0]<<endl;
    }
    bec0 = bec ; layer0 = layer ; phi0 = phi ; eta0 = eta; 
    sum = data;
    idata = 1;
  } else {
    sum += data;
    idata++;
  }
}
fin.close();
//------------------read in TM1---------------------------
ifstream fin("sct_TM1_000016.txt");
while(fin>>bec>>layer>>phi>>eta>>data>>date>>time)
{
  if(bec!=bec0 || layer!=layer0 || phi!=phi0 || eta!=eta0){
    if(idata != 0) {
      if(bec0==1) eta0 += 6;
      Temp1[bec0][layer0][eta0][phi0] = sum/idata;
//    if(bec0==1 && layer0 ==3 && eta0 == 0) 
//   cout<<bec0<<" "<<layer0<<" "<<eta0-6<<" "<<phi0<<" "
//    <<Temp1[bec0][layer0][eta0][phi0]<<endl;
    }
    bec0 = bec ; layer0 = layer ; phi0 = phi ; eta0 = eta; 
    sum = data;
    idata = 1;
  } else {
    sum += data;
    idata++;
  }
}
fin.close();
//-----------take average of two hybrid temperature of barrel modules ----------- 
int BEC=1; 
for(int Layer=3; Layer<7; Layer++) {
  for(int Eta=0; Eta<13 ;Eta++) {
    if( Eta == 6 ) continue;
    int nphi = barrel[Layer];
    for (int iphi=0; iphi<nphi; iphi++) {
//       cout<<BEC<<" "<<Layer<<" "<<Eta<<" "<<iphi<<" : "<<
//       Temp[BEC][Layer][Eta][iphi] <<" "<<Temp1[BEC][Layer][Eta][iphi]<<endl;
       if( Temp[BEC][Layer][Eta][iphi] > TEMPcut_low) {
           if (Temp1[BEC][Layer][Eta][iphi] > TEMPcut_low) {
             Temp[BEC][Layer][Eta][iphi]+= Temp1[BEC][Layer][Eta][iphi];
             Temp[BEC][Layer][Eta][iphi] /= 2.;
           }  
        } else {
           Temp[BEC][Layer][Eta][iphi] = Temp1[BEC][Layer][Eta][iphi];
        }
//      cout<<BEC<<" "<<Layer<<" "<<Eta<<" "<<iphi<<" ---->" <<
//       Temp[BEC][Layer][Eta][iphi] <<endl;
    }
  }
}
//--------------------read in HVVoltage-------------------------
int bec0=9,layer0=9,phi0=9,eta0=9,date0=9,time0=9;
double data, data0=-999.;
int idata=0;
double sum=0.;
ifstream fin("sct_HVVoltage_000016.txt");
while(fin>>bec>>layer>>phi>>eta>>data>>date>>time)
{
  if(bec!=bec0 || layer!=layer0 || phi!=phi0 || eta!=eta0){
    if(idata != 0) {
      if(bec0==1) eta0 += 6;
      HV[bec0][layer0][eta0][phi0] = sum/idata;
//    if(bec0==1 && layer0 ==3 && eta0 == 0) 
//    cout<<bec0<<" "<<layer0<<" "<<eta0-6<<" "<<phi0
//     <<" "<<HV[bec0][layer0][eta0][phi0]<<endl;
    }
    bec0 = bec ; layer0 = layer ; phi0 = phi ; eta0 = eta; 
    sum = data;
    idata = 1;
  } else {
    sum += data;
    idata++;
  }
}
fin.close();
//------------------read in HVCurrent---------------------------
int bec0=9,layer0=9,phi0=9,eta0=9,date0=9,time0=9;
double data, data0=-999.;
int idata=0;
double sum=0.;
ifstream fin("sct_HVCurrent_000016.txt");
while(fin>>bec>>layer>>phi>>eta>>data>>date>>time)
{
  if(bec!=bec0 || layer!=layer0 || phi!=phi0 || eta!=eta0){
    if(idata != 0) {
      if(bec0==1) eta0 += 6;
      HVLC[bec0][layer0][eta0][phi0] = sum/idata;
//    if(bec0==1 && layer0 ==3 && eta0 == 0) 
//    cout<<bec0<<" "<<layer0<<" "<<eta0-6<<" "<<phi0<<" "
//    <<HVLC[bec0][layer0][eta0][phi0]<<endl;
    }
    bec0 = bec ; layer0 = layer ; phi0 = phi ; eta0 = eta; 
    sum = data/1000.;
    idata = 1;
  } else {
    sum += data/1000.;
    idata++;
  }
}
fin.close();
//------------------Read in CiS module numbers---------------------------
int nCiS = 0.;
double dummy1,dummy2;
int becCiS[1000],layerCiS[1000],phiCiS[1000],etaCiS[1000];
ifstream fin1("CiS.txt");
while(fin1>>dummy1>>dummy2>>becCiS[nCiS]>>layerCiS[nCiS]>>phiCiS[nCiS]>>etaCiS[nCiS]){
  etaCiS[nCiS] -= 6; // change eta numbering from 6-8 to 0-2.
  layerCiS[nCiS] += 1; // change Disk numbering from 0-8 to 1-9. 
  if(nCiS<10) { cout<<"CiS="<<nCiS<<" "<<becCiS[nCiS]<<" "<<layerCiS[nCiS]<<" "<<  
        phiCiS[nCiS]<<" "<<etaCiS[nCiS]<<endl;}
  nCiS++;
}
cout << "number of CiS modules ="<< nCiS +1 <<endl;
fin1.close();
//-----plot individual phi distribution----------------------
int icd;
//-----------plot barrel----------------
int BEC=1; 
for(int Layer=0; Layer<10; Layer++) {
  if(Layer < 3 || Layer > 6) continue;
  for(int Eta=0; Eta<13 ;Eta++) {
    int etavalue = Eta -6;
    if( Eta == 6 ) continue;
    icd = Eta + 1;
    if( Eta > 6 ) icd = Eta - 6 ;
    c1 -> cd(icd);
    int nphi = barrel[Layer];
    for (int iphi=0; iphi<nphi; iphi++) {
//------- cut out bad modules -------------------------
       if( HV[BEC][Layer][Eta][iphi] <= HVcut) continue;
       double Tmodule = Temp[BEC][Layer][Eta][iphi];       
       if( Tmodule <= TEMPcut_low) continue;
       if( Tmodule >= TEMPcut_high) continue;
//--------------------- calculate to -10C ----------------------- 
       Tmodule -= DT_barrel; // T difference between hybrid and module (FEM analysis)
       double W= pow((T0-10.)/(T0+Tmodule),2)
                 * exp(Eg/2./kB*(1./(T0+Tmodule)-1./(T0-10.)));
       HVLC_norm[BEC][Layer][Eta][iphi] = W * HVLC[BEC][Layer][Eta][iphi];
    }
    h1 = new TGraph(nphi,Phi,HVLC_norm[BEC][Layer][Eta]);
    h1->SetMarkerSize(0.5);
    h1->SetMaximum(10.);
    h1->SetMinimum(0.);
    h1->GetXaxis()->SetTitle("phi"); h1->GetXaxis()->SetTitleOffset(0.8);
    h1->GetXaxis()->SetTitleSize(0.10);
    h1->GetXaxis()->SetLabelSize(0.08);    
    h1->GetYaxis()->SetTitle("I_{leak}@-10^{o}C (#muA)"); 
    h1->GetYaxis()->SetTitleOffset(0.5);
    h1->GetYaxis()->SetTitleSize(0.10);
    h1->GetYaxis()->SetLabelSize(0.08);
    h1->SetMarkerColor(2);
    h1->GetXaxis()->SetLimits(0.,nphi);
    h1->Draw("ALP");
    cid = new char[40];
    sprintf(cid,"Barrel B%d, #eta %d",Layer,etavalue);
    cid[39]=0;
    yleg = 0.78;
    leg = new TLegend(0.55 ,yleg ,0.85,yleg+0.12,cid);
    leg -> SetFillColor(10);
    leg->SetBorderSize(0); leg->SetTextSize(0.09); leg->Draw();
    cfile = new char[25];sprintf(cfile,"HVLC_Barrel_norm_%d_1.pdf",Layer);cfile[24]=0;
    if( Eta == 5 ) c1->Print(cfile);
    cfile = new char[25];sprintf(cfile,"HVLC_Barrel_norm_%d_2.pdf",Layer);cfile[24]=0;
    if( Eta == 12 )c1->Print(cfile);
  }
}
//-------------plot endcaps--------------------
double Phi_CiS[60], HVLC_CiS[60];
for(int BEC=0; BEC<3; BEC+=2){ 
  for(int Eta=0; Eta<3 ;Eta++) {
  c1->Clear();
  c1->Divide(3,3);
    for(int Layer=1; Layer<10; Layer++) {
      c1 -> cd(Layer);
      if(HVLC[BEC][Layer][Eta][0]<-100.) continue;
      int nphi = 52;
      if(Eta !=0) nphi=39;
      int kCiS = 0;
      for (int iphi=0; iphi<nphi; iphi++) {
//------- cut out bad modules -------------------------
        if( HV[BEC][Layer][Eta][iphi] <= HVcut) continue;
        double Tmodule = Temp[BEC][Layer][Eta][iphi];       
        if( Tmodule <= TEMPcut_low) continue;
        if( Tmodule >= TEMPcut_high) continue;
//--------------------- calculate to -10C ----------------------- 
//-------------T difference between hybrid and module (FEM analysis)
        if( Eta==0) Tmodule -= DT_ECouter;
        if( Eta==1) Tmodule -= DT_ECmiddle;
        if( Eta==2) Tmodule -= DT_ECinner;
        double W= pow((T0-10.)/(T0+Tmodule),2)
                  * exp(Eg/2./kB*(1./(T0+Tmodule)-1./(T0-10.)));
        HVLC_norm[BEC][Layer][Eta][iphi] = W * HVLC[BEC][Layer][Eta][iphi];
      }
//------ check if CiS modules are contained---------------------
       if(Eta != 0 ){
         for (int iphi=0; iphi<nphi; iphi++) {
           for (int iCiS=0; iCiS<nCiS ; iCiS++) {
             if( BEC != becCiS[iCiS]) continue;
             if( Layer != layerCiS[iCiS]) continue;
             if( Eta != etaCiS[iCiS]) continue;
             if( iphi != phiCiS[iCiS]) continue;
//-------------- CiS module found -------------------
              Phi_CiS[kCiS] = iphi;
              HVLC_CiS[kCiS] = HVLC_norm[BEC][Layer][Eta][iphi]; 
              kCiS++;
           }
         }         
       }
       cout<<"CiS module number = "<<kCiS<<endl;
//-----------------now plot graphs--------------------------------
      h1 = new TGraph(nphi,Phi,HVLC_norm[BEC][Layer][Eta]);
      if( kCiS > 0 ) {
        h2 = new TGraph(kCiS,Phi_CiS,HVLC_CiS);
        h2->SetMarkerSize(0.5);h2->SetMarkerColor(9);
      }
      h1->SetMarkerSize(0.5);
      double ymax=10.;
      if(Eta==2) ymax=6.;
      h1->SetMaximum(ymax);
      h1->SetMinimum(0.);
      h1->GetXaxis()->SetTitle("phi"); h1->GetXaxis()->SetTitleOffset(0.8);
      h1->GetXaxis()->SetTitleSize(0.10);
      h1->GetXaxis()->SetLabelSize(0.08);    
      h1->GetYaxis()->SetTitle("I_{leak}@-10^{o}C (#muA)"); 
      h1->GetYaxis()->SetTitleOffset(0.7);
      h1->GetYaxis()->SetTitleSize(0.07);
      h1->GetYaxis()->SetLabelSize(0.08);
      h1->SetMarkerColor(2);
      h1->GetXaxis()->SetLimits(0.,nphi);
      h1->Draw("ALP");
      if( kCiS !=0 ) h2-> Draw("SAME P");
//--------------Hamamatsu vs CiS--------------------
    t1 = new TLatex(2,ymax*0.85,"#bullet Hamamatsu"); t1->SetTextColor(2); 
    t2 = new TLatex(2,ymax*0.75,"#bullet CiS"); t2->SetTextColor(9);
    t1->SetTextSize(0.06); t2->SetTextSize(0.06);
    t1->Draw(); t2->Draw();
//---------------Layer and outer/middle/inner------------
      cid = new char[40];    cid[39]=0;
      if(BEC==0){
        sprintf(cid,"ECC Layer %d, outer",Layer);
        if(Eta==1) sprintf(cid,"ECC Layer %d middle",Layer);
        if(Eta==2) sprintf(cid,"ECC Layer %d inner",Layer);
      } else {
        sprintf(cid,"ECA Layer %d, outer",Layer);
        if(Eta==1) sprintf(cid,"ECA Layer %d middle",Layer);
        if(Eta==2) sprintf(cid,"ECA Layer %d inner",Layer);
      }
      leg = new TLegend(0.42 ,0.78 ,0.85,0.90,cid);leg -> SetFillColor(10);
      leg->SetBorderSize(0); leg->SetTextSize(0.07); leg->Draw();   
    }
  cfile = new char[25]; cfile[24]=0;
  if(BEC==0) sprintf(cfile,"HVLC_ECC_norm_%d.pdf",Eta);
  if(BEC==2) sprintf(cfile,"HVLC_ECA_norm_%d.pdf",Eta);
  c1->Print(cfile);
  }
}
//------------------------------------------   
}