// gROOT->Reset();

// Christian Klein-Boesing CKB 
// mailto:stevero@ikp.uni-muenster.de

// Macro to extract slopes and correction factors
// on root prompt load this macro 
// root[0] .L GetCF.C
// and type
// root[1] Draw_Fit();
// which will fit the slopes in the specified slopesFile and 
// write them out into a file xslopes.dat
// this file can be used to determine correction factors with 
// root[2] angular_all("xslopes.dat") 
// this will determine the correction factors realtive to a angular correctet
// mean of the slope for all PbGl sectors (angular_Re1 will do it 
// relative to E1) you will get an output file which can be used
// in the same way as PbSc mipcorrection files
// Another usefull method of this macro is
// CombineTwoCF(char *file1,char *file2,int invert1 = 0,float globalfactor = 1.) which allows you to cobine different correction factors (as the name says :-)






int MakeSwkey(int,int,int,int);
 
// PbGl geometry parameters

const int kArraySize = 9216;
const int kSec = 2;
const int kMaxY = 48;
const int kMaxZ = 96;
double SLOPES[2][kMaxY][kMaxZ];
const  Double_t x0=543.2;
const  Double_t yStep=4.1105;
const  Double_t zStep=4.1040;


const double PI = 3.1415927;
const float fMaxCorr = 2.; // allow only correction factors less than 2 otherwise one get's extra singers !!!






// char *slopesFile = "/phenix/data24/kleinb/run02/slopes/corr5/slopes_5.root";
char *slopesFile = "/phenix/data24/kleinb/run02/slopes/corr2_2/slopes2_2.root";

char *foutSlopes = "xslopes.dat";
char *fCFraw     = "xCFraw.dat";
char *fCFang_corr= "xCF_ang_corr.dat"; 

TH2F *angleE[2] = {2*0};



void Draw_Fit(){

  Int_t y_ind,z_ind;
  double slope;
  bool fail;


  gROOT->SetStyle("Plain");


 
  TCanvas* mycanvas = new TCanvas("mycanvas","",200,10,700,500); 
  mycanvas->SetLogy();

  TFile *f1 = new TFile(slopesFile);

  TH2F *hslopes = (TH2F*)f1->Get("hslope");

  TH1F *Ecal_Tower[4608];

 
  TH1D *hsE0 = new TH1D("hsE0","Slopes E0",400,-80,20);
  TH2D *secE0 = new TH2D("secE0","Slopes E0",kMaxZ,-0.5,(1.*kMaxZ)-0.5,
			 kMaxY,-0.5,(1.*kMaxY)-0.5);
  
  TH1D *hsE1 = new TH1D("hsE1","Slopes E1",400,-80,20);
  TH2D *secE1 = new TH2D("secE1","Slopes E1",kMaxZ,-0.5,(1.*kMaxZ)-0.5,
			 kMaxY,-0.5,(1.*kMaxY)-0.5);
  slope = -9999999;
  fail = true;

  Char_t next;
  
  for(int index=0;index<kArraySize;index++){
    // for(int index=0;index<150;index++){

    int nbin = (hslopes->GetXaxis())->FindBin(1.*index);
    TH1D *h1 = hslopes->ProjectionY("h1",nbin,nbin,"E");


    int sec = index/4608;
    int y_ind = index/kMaxZ-sec*kMaxY;
    int z_ind = index-kMaxZ*y_ind-4608*sec;
 
    //    cout << sec << " "<< y_ind << " " << z_ind << endl;
   
    // continue;
    h1->Rebin(2);
    slope = mfit(h1);

    
    // if(!index%400){
    /*
      h1->GetXaxis()->SetRange(0,100);
      h1->Draw();
      mycanvas->Update();
      char c = getchar();
      if(c=='q')return;
    */
    // }
    
      cout << "SLOPE  " << slope << " INDEX " << index << endl;
    SLOPES[sec][y_ind][z_ind] = slope;
    // testing of coordinates
    /*
    if(sec==0&&y_ind%4){
      SLOPES[sec][y_ind][z_ind] = 20;
    }
    if(sec==1&&z_ind%4){
      SLOPES[sec][y_ind][z_ind] = 20;
    }
    */
    if(slope!=-9999999){
      if(sec==1){
	hsE1->Fill(slope);
	secE1->Fill(z_ind,y_ind,slope);
      }
      if(sec==0){
	hsE0->Fill(slope);
	secE0->Fill(z_ind,y_ind,slope);
      }
	 
    }
    //   Ausgabe << y_ind << " " << z_ind << " " << slope << endl;
  }// for ... index < ArraySize
    

 
// secE0->Draw("colz");
   
  ofstream out1(foutSlopes);
  for(int i = 0;i<kArraySize;i++){
    int sec = i/4608;
    int y_ind = i/kMaxZ-sec*kMaxY;
    int z_ind = i-kMaxZ*y_ind-4608*sec;
    int arm = 1;
    out1 << MakeSwkey(arm,sec,y_ind,z_ind) << " " << SLOPES[sec][y_ind][z_ind] << endl;

  }
  out1.close();
  
  hsE0->Fit("gaus","E0");
  hsE0->Draw();
 
  double Ref[2];
  Ref[0] = gaus->GetParameter(1);

  hsE1->Fit("gaus","E0");
  hsE1->Draw("same");

  Ref[1] = gaus->GetParameter(1);


  ofstream out2(fCFraw);

  for(int i = 0;i<kArraySize;i++){
    int sec = i/4608;
    int y_ind = i/kMaxZ-sec*kMaxY;
    int z_ind = i-kMaxZ*y_ind-4608*sec;
    int arm = 1;

    int status = 1;
    float cf_err = 0;
    float cf = 1.;
    
    if(SLOPES[sec][y_ind][z_ind]!=0&&(SLOPES[sec][y_ind][z_ind]>fMaxCorr*Ref[sec])){
      cf = Ref[sec]/SLOPES[sec][y_ind][z_ind];
    }
    else{
      status = 0;
    }
    
    out2 << MakeSwkey(arm,sec,y_ind,z_ind) << " " << cf << " " << cf_err 
	 << " " << status << endl;
  }
  out2.close();
  

}

int MakeSwkey(int arm, int sec, int iy, int iz)
{

  /* This is an encoded (arm, sector, iy, iz) of the central tower of the
     cluster. The encoding is done as follows:

     swkey = 100000*emc->arm    // 0 or 1
     +      10000*emc->sector   // 0 - 3
     +        100*emc->ind[1]   // iy 0-35 or 0-47
     +            emc->ind[0];  // iz 0-71 or 0-95

   */

  int swkey = 100000 * arm + 10000 * sec + 100 * iy + iz;
  return swkey;

}

void DecodeSwkey(int swkey,int &arm,int &sec,int &iy,int &iz){
  
  arm    = swkey/100000;
  swkey -= arm*100000;
  sec    =  swkey/10000;
  swkey -= sec*10000;
  iy     = swkey/100;
  swkey -= iy*100;
  iz     = swkey;
  return;
}


float mfit(TH1* hist,int drawflag = 0)
{
  
  //  double r1=3.16507354941686489e-03;
  //  double r2=1.95527179729018913e-04;

  // Maxims Values for Pbgl
  
  /*
  double r1=1.30561635349989678e-03;
  double r2=8.73425562667767868e-05;
  */  

    // PbSc

  
  // year02 due to Zero supression some Spectra have  noise peak
  // some don't fix lower fitlimit to 0.15 (below)


  float start = 0.15;
  int startbin = hist->FindBin(start);
  int startval = hist->GetBinContent(startbin);
  double rat = 1e-01; // fit over 2 decades

  double r1=1.e-02;
  double r2=1.e-04;
  

  double slope;

  int nbins=hist->GetNbinsX();
  double max=hist->GetMaximum();
  int maxbin=hist->GetMaximumBin();

  //  cout<<nbins<<" "<<max<<" "<<maxbin<<endl;
  cout<<"upper thresh "<<max*r1<<endl;

  int up3bins[15], low3bins[15];
  unsigned short upp=0;
  unsigned short lowp=0;

  //=====> search for upper crossing
  for(int ibin=maxbin+1; ibin<=nbins; ibin++){

    double binval=hist->GetBinContent(ibin);
    if(binval>0. && binval<max*r1){

      up3bins[upp]=ibin;
      upp++;

    } else if (binval>0.) upp=0;

    if(upp==5) break;

  }


 
  


  cout<<"upp= "<<upp<<endl;
  //  for(int ind=0; ind<15; ind++) cout<<hist->GetBinCenter(up3bins[ind])<<endl;



  //=====> search for lower crossing
  for(int ibin=maxbin+1; ibin<=nbins; ibin++){

    double binval=hist->GetBinContent(ibin);
    if(binval>0. && binval<max*r2){

      low3bins[lowp]=ibin;
      lowp++;

    } else if(binval>0.) lowp=0;

    if(lowp==5) break;

  }


  // CKB 

  upp = 5;
  up3bins[0] = startbin;

  lowp = 0;

  for(int ibin=maxbin+1; ibin<=nbins; ibin++){
    
    double binval=hist->GetBinContent(ibin);
    //  cout << binval << " " << startval*rat << endl;
    if(binval>0. && binval< (startval*rat)){

      // cout << ibin << endl;

      low3bins[lowp]=ibin;
      lowp++;

    } else if(binval>0.) lowp=0;

    if(lowp==5) break;

  }

  // BKC

  cout<<"lowp= "<<lowp<<endl;
  //  for(int ind=0; ind<15; ind++) cout<<hist->GetBinCenter(low3bins[ind])<<endl;

  double leftlimit, rightlimit;
  bool fail=false;

  if(upp==5 && lowp==5){ // both crossings found

    leftlimit=hist->GetBinCenter(up3bins[0]);
    rightlimit=hist->GetBinCenter(low3bins[0]);

  } else if(upp==5 && lowp<5){ // only upper crossing found

    leftlimit=hist->GetBinCenter(up3bins[0]);
    rightlimit=hist->GetBinCenter(nbins);
    // CKB
    fail = true;

    if(nbins-up3bins[0]<100) fail=true;

  } else fail=true;

  if(!fail){
    hist->Fit("expo", "EQ0", "", leftlimit, rightlimit);
    TF1 *fitfunc=hist->GetFunction("expo");
    slope=fitfunc->GetParameter(1);
    if(drawflag){
      hist->GetXaxis()->SetRange(0,up3bins[0]+30);
      mycanvas->Update();
      getchar();
    }
  }
  else {
    cout<<" MANUAL INTERVENTION NEEDED"<<endl;
    slope = -9999999;
  }
      return slope;
  //  return fail;
}

void DrawSlopes(){

  ifstream in0,in1;
  double value;

  TH1D *h0 = new TH1D("h0","Slopes",200,-60,20);
  TH1D *h1 = new TH1D("h1","Slopes",200,-60,20);

  TH2D *hsec0 = new TH2D("hsec0","Slopes",MaxZ,0.,(1.*MaxZ),MaxY,0.,(1.*MaxY));
  TH2D *hsec1 = new TH2D("hsec1","Slopes",MaxZ,0.,(1.*MaxZ),MaxY,0.,(1.*MaxY));
  


//    h1->SetXTitle("-k");
//    h1->SetLineColor(2);
//    h1->Draw();
  
  h0->SetXTitle("-k");
  h0->SetLineColor(6);
  h0->Draw("");
  h12->Draw("col");
}

///////////////////////////////////////////////////////////////////////////
// Here comes the Stuff for the angular Dependence
///////////////////////////////////////////////////////////////////////////


Double_t fitf(Double_t *x, Double_t *par){
  // Fit function for angular dependence of correction
  Double_t fitval;
  Double_t arg1, arg2;
  // arg1 = par[1]*x[0];
  //arg2 = par[3]*x[0];
  //fitval = par[0]*(TMath::Exp(arg1)) +1.-par[2] + par[4]*sin(arg2);
  //fitval par[0]*(TMath::Exp(arg1)) +1.-par[2];
  //  fitval = par[0]*(TMath::Exp(arg1))+par[2];
  fitval = par[0] + par[1] * x[0] + par[2] *x[0]*x[0];
  return fitval;
}

void angular(const char* inSlopes){
  // Start again with angiular correction
  // 1. Get the Angular dependence of inverse Slopes
  
  GetAngles();
  
  const Int_t nba = 25;
  Double_t xla = 0.;
  Double_t xua = 25.;
  Double_t wba;
  wba = (xua - xla)/static_cast<Int_t>(nba);
 
 
  // Fit histo 1
  Double_t min=0.;
  Double_t max=22.;
  // TF1 *fitf[2] = {2*0}; 
  // fitf[0] = new TF1("f1",fitf,min,max,3);

  TF1 *f1 = new TF1("f1",fitf,min,max,3);
  TF1 *f0 = new TF1("f0",fitf,min,max,3);

  TF1 *fitf[2] = {f0,f1};



  // Book all that's needed for the inverse Slopes
  
  TH1D *inv1[2] = {2*0};
  TH2D *inv2[2] = {2*0};
  TProfile *inv3[2] = {2*0};

  TH1D *inv11[2] = {2*0.};
  TH2D *inv12[2] = {2*0};
  TProfile *inv13[2] = {2*0};

  TH1D *inv110[2] = {2*0.};




  inv1[0] = new TH1D("inv1E0","Inverse slopes",500,-1.,0.);
  inv2[0] = new TH2D("inv2E0","Inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv3[0] = new TProfile("inv3E0","Inv Slope vs angle",nba,xla,xua,""); 
  inv11[0] = new TH1D("inv11E0","corr inverse slopes",500,-1.,0.);
  inv12[0] = new TH2D("inv12E0"," corr inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv13[0] = new TProfile("inv13E0","corr inv Slope vs angle",nba,xla,xua,""); 
  inv110[0] = new TH1D("inv110E0","CF from ang Corr inverse Slopes",400,0.,4.);
 

  inv1[1] = new TH1D("inv1E1","Inverse slopes",500,-1.,0.);
  inv2[1] = new TH2D("inv2E1","Inverse slopes vs angle",nba,xla,xua,500,-1.,0.);

  inv3[1] = new TProfile("inv3E1","Inv Slope vs angle",nba,xla,xua,""); 
  TH1D *inv11[1] = new TH1D("inv11E1","corr inverse slopes",500,-1.,0.);
  TH2D *inv12[1] = new TH2D("inv12E1"," corr inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv13[1] = new TProfile("inv13E1","corr inv Slope vs angle",nba,xla,xua,""); 
  inv110[1] = new TH1D("inv110E1","CF from ang Corr inverse Slopes",400,0.,4.);

  TH1F *htemp = new TH1F("htemp","",200,-50,0);

  // Open the Slopes File 




  
  ifstream in1;  
  in1.open(inSlopes,ios::in);
  
  
  int swkey;
  float value;

  if(in1.good()){
    cout << "Opening " << inSlopes << endl;
    while(in1 >> swkey >> value){
     
      int iy,iz,arm,sec;

      DecodeSwkey(swkey,arm,sec,iy,iz);

      SLOPES[sec][iy][iz] = value;
            
      if(value!=0&&value>-1000){
	
	int nbin = angleE[sec]->FindBin(iz,iy);
	float theta = angleE[sec]->GetBinContent(nbin);
	if(theta>0&&theta<2){
	  htemp->Fill(value);
	}

	inv1[sec]->Fill(1./value);
	inv2[sec]->Fill(theta,1./value);
	inv3[sec]->Fill(theta,1./value);
       
      }
    } // kArraySize

    fitf[0]->SetParameters(-0.1,1,1);
    fitf[1]->SetParameters(-0.1,1,1);
    

    inv3[0]->Draw();
   
    inv3[0]->Fit("f0","R");
    inv3[1]->Fit("f1","R");
   

    // Correct the inv Slopes for angular Dependence
  
    for(int sec = 0;sec < 2;sec++){
      for(int iy=0;iy<kMaxY;iy++){
	for(int iz=0;iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin);
	  
	  // cout << sec << " " << iy << " " << iz << endl;
	  if(SLOPES[sec][iy][iz]!=0&&SLOPES[sec][iy][iz]>-1000){
	    
	  
	    float invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf[sec]->Eval(theta)*fitf[sec]->Eval(0);
	    inv11[sec]->Fill(invSlopeNew);
	    inv12[sec]->Fill(theta,invSlopeNew);
	    inv13[sec]->Fill(theta,invSlopeNew);
	  }	
	} // iy < Maxy
      } // iz < Maxz
    }// sec < 2
  
    inv11[0]->Fit("gaus");
    double MeanInv11[2];

    MeanInv11[0] = gaus->GetParameter(1);
    inv11[1]->Fit("gaus");
    MeanInv11[1] = gaus->GetParameter(1);

    // And loop again for CF
    ofstream   outCF1;
    outCF1.open(fCFang_corr,ios::out);
    
    for(int sec = 0;sec<2;sec++){
      for(int iy=0; iy<kMaxY; iy++){
	for(int iz=0; iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin); 
	  
	  float invSlopeNew;
	  float CFInvAngCorr = 1;
	  float cf_err = 0;
	  int status = 1;
	  int arm = 1;
	  
	  if(SLOPES[sec][iy][iz]!=0){
	    invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf[sec]->Eval(theta)*fitf[sec]->Eval(0);
	    CFInvAngCorr = MeanInv11[sec]/invSlopeNew;
	    if(CFInvAngCorr<=fMaxCorr){
	      inv110[sec]->Fill(CFInvAngCorr);
	    }	
	    else{
	      status = 0;
	      CFInvAngCorr = 1;
	    }
	  }
	  else{
	    status = 0;
	  }	  
	  outCF1 <<MakeSwkey(arm,sec,iy,iz) << " " << CFInvAngCorr 
		 << " " << cf_err << " " << status <<  endl;      
	} // iy < Maxy
      } // iz < Maxz     
    }// sec < 2
  } // in1.good
  
  in1.close();
  outCF1.close();

}


void angular_Ref1(const char* inSlopes){
  // Start again with angiular correction
  // 1. Get the Angular dependence of inverse Slopes
  
  GetAngles();
  
  const Int_t nba = 25;
  Double_t xla = 0.;
  Double_t xua = 25.;
  Double_t wba;
  wba = (xua - xla)/static_cast<Int_t>(nba);
 
 
  // Fit histo 1
  Double_t min=0.;
  Double_t max=22.;
  // TF1 *fitf[2] = {2*0}; 
  // fitf[0] = new TF1("f1",fitf,min,max,3);

  TF1 *f1 = new TF1("f1",fitf,min,max,3);
  TF1 *f0 = new TF1("f0",fitf,min,max,3);

  TF1 *fitf[2] = {f0,f1};



  // Book all that's needed for the inverse Slopes
  
  TH1D *inv1[2] = {2*0};
  TH2D *inv2[2] = {2*0};
  TProfile *inv3[2] = {2*0};

  TH1D *inv11[2] = {2*0.};
  TH2D *inv12[2] = {2*0};
  TProfile *inv13[2] = {2*0};

  TH1D *inv110[2] = {2*0.};




  inv1[0] = new TH1D("inv1E0","Inverse slopes",500,-1.,0.);
  inv2[0] = new TH2D("inv2E0","Inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv3[0] = new TProfile("inv3E0","Inv Slope vs angle",nba,xla,xua,""); 
  inv11[0] = new TH1D("inv11E0","corr inverse slopes",500,-1.,0.);
  inv12[0] = new TH2D("inv12E0"," corr inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv13[0] = new TProfile("inv13E0","corr inv Slope vs angle",nba,xla,xua,""); 
  inv110[0] = new TH1D("inv110E0","CF from ang Corr inverse Slopes",400,0.,4.);
 

  inv1[1] = new TH1D("inv1E1","Inverse slopes",500,-1.,0.);
  inv2[1] = new TH2D("inv2E1","Inverse slopes vs angle",nba,xla,xua,500,-1.,0.);

  inv3[1] = new TProfile("inv3E1","Inv Slope vs angle",nba,xla,xua,""); 
  TH1D *inv11[1] = new TH1D("inv11E1","corr inverse slopes",500,-1.,0.);
  TH2D *inv12[1] = new TH2D("inv12E1"," corr inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv13[1] = new TProfile("inv13E1","corr inv Slope vs angle",nba,xla,xua,""); 
  inv110[1] = new TH1D("inv110E1","CF from ang Corr inverse Slopes",400,0.,4.);

  TH1F *htemp = new TH1F("htemp","",200,-50,0);

  // Open the Slopes File 




  
  ifstream in1;  
  in1.open(inSlopes,ios::in);
  
  
  int swkey;
  float value;

  if(in1.good()){
    cout << "Opening " << inSlopes << endl;
    while(in1 >> swkey >> value){
     
      int iy,iz,arm,sec;

      DecodeSwkey(swkey,arm,sec,iy,iz);

      SLOPES[sec][iy][iz] = value;
            
      if(value!=0&&value>-1000){
	
	int nbin = angleE[sec]->FindBin(iz,iy);
	float theta = angleE[sec]->GetBinContent(nbin);
	if(theta>0&&theta<2){
	  htemp->Fill(value);
	}

	inv1[sec]->Fill(1./value);
	inv2[sec]->Fill(theta,1./value);
	inv3[sec]->Fill(theta,1./value);
       
      }
    } // kArraySize

    fitf[0]->SetParameters(-0.1,1,1);
    fitf[1]->SetParameters(-0.1,1,1);
    

    inv3[0]->Draw();
   
    inv3[0]->Fit("f0","R");
    inv3[1]->Fit("f1","R");
   

    // Correct the inv Slopes for angular Dependence
  
    for(int sec = 0;sec < 2;sec++){
      for(int iy=0;iy<kMaxY;iy++){
	for(int iz=0;iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin);
	  
	  // cout << sec << " " << iy << " " << iz << endl;
	  if(SLOPES[sec][iy][iz]!=0&&SLOPES[sec][iy][iz]>-1000){
	    
	  
	    float invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf[1]->Eval(theta)*fitf[1]->Eval(0);
	    inv11[sec]->Fill(invSlopeNew);
	    inv12[sec]->Fill(theta,invSlopeNew);
	    inv13[sec]->Fill(theta,invSlopeNew);
	  }	
	} // iy < Maxy
      } // iz < Maxz
    }// sec < 2
  
    inv11[0]->Fit("gaus");
    double MeanInv11[2];

    MeanInv11[0] = gaus->GetParameter(1);
    inv11[1]->Fit("gaus");
    MeanInv11[1] = gaus->GetParameter(1);

    // And loop again for CF
    ofstream   outCF1;
    outCF1.open(fCFang_corr,ios::out);
    
    for(int sec = 0;sec<2;sec++){
      for(int iy=0; iy<kMaxY; iy++){
	for(int iz=0; iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin); 
	  
	  float invSlopeNew;
	  float CFInvAngCorr = 1;
	  float cf_err = 0;
	  int status = 1;
	  int arm = 1;
	  
	  if(SLOPES[sec][iy][iz]!=0){
	    invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf[1]->Eval(theta)*fitf[1]->Eval(0);
	    CFInvAngCorr = MeanInv11[1]/invSlopeNew;
	    if(CFInvAngCorr<=fMaxCorr){
	      inv110[sec]->Fill(CFInvAngCorr);
	    }	
	    else{
	      status = 0;
	      CFInvAngCorr = 1;
	    }
	  }
	  else{
	    status = 0;
	  }	  
	  outCF1 <<MakeSwkey(arm,sec,iy,iz) << " " << CFInvAngCorr 
		 << " " << cf_err << " " << status <<  endl;      
	} // iy < Maxy
      } // iz < Maxz     
    }// sec < 2
  } // in1.good
  
  in1.close();
  outCF1.close();

}


void angular_all(const char* inSlopes){
  // Start again with angiular correction
  // 1. Get the Angular dependence of inverse Slopes
  

  cout << "fMaxCorr " << fMaxCorr << endl;

  GetAngles();
  
  const Int_t nba = 25;
  Double_t xla = 0.;
  Double_t xua = 25.;
  Double_t wba;
  wba = (xua - xla)/static_cast<Int_t>(nba);
 
 
  // Fit histo 1
  Double_t min=0.;
  Double_t max=22.;
  // TF1 *fitf[2] = {2*0}; 
  TF1 *fitf = new TF1("f1",fitf,min,max,3);





  // Book all that's needed for the inverse Slopes
  
  TH1D *inv1 = 0;
  TH2D *inv2 = 0;
  TProfile *inv3 = 0;

  TH1D *inv11 = 0;
  TH2D *inv12 = 0;
  TProfile *inv13 = 0;

  TH1D *inv110 = 0;




  inv1 = new TH1D("inv1E0E1","Inverse slopes",500,-1.,0.);
  inv2 = new TH2D("inv2E0E1","Inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv3 = new TProfile("inv3E0E1","Inv Slope vs angle",nba,xla,xua,""); 
  inv11 = new TH1D("inv11E0E1","corr inverse slopes",500,-1.,0.);
  inv12 = new TH2D("inv12E0E1"," corr inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv13 = new TProfile("inv13E0E1","corr inv Slope vs angle",nba,xla,xua,""); 
  inv110 = new TH1D("inv110E0E1","CF from ang Corr inverse Slopes",400,0.,4.);
 
  TH1F *htemp = new TH1F("htemp","",200,-50,0);

  // Open the Slopes File 




  
  ifstream in1;  
  in1.open(inSlopes,ios::in);
  
  
  int swkey;
  float value;

  if(in1.good()){
    cout << "Opening " << inSlopes << endl;
    while(in1 >> swkey >> value){
     
      int iy,iz,arm,sec;

      DecodeSwkey(swkey,arm,sec,iy,iz);

      SLOPES[sec][iy][iz] = value;
            
      if(value!=0&&value>-1000){
	
	int nbin = angleE[sec]->FindBin(iz,iy);
	float theta = angleE[sec]->GetBinContent(nbin);
	if(theta>0&&theta<2){
	  htemp->Fill(value);
	}

	inv1->Fill(1./value);
	inv2->Fill(theta,1./value);
	inv3->Fill(theta,1./value);
       
      }
    } // kArraySize

    fitf->SetParameters(-0.1,1,1);
   
    

    inv3->Draw();
   
    inv3->Fit("f1","R");
   
   

    // Correct the inv Slopes for angular Dependence
  
    for(int sec = 0;sec < 2;sec++){
      for(int iy=0;iy<kMaxY;iy++){
	for(int iz=0;iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin);
	  
	  // cout << sec << " " << iy << " " << iz << endl;
	  if(SLOPES[sec][iy][iz]!=0&&SLOPES[sec][iy][iz]>-1000){
	    
	  
	    float invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf->Eval(theta)*fitf->Eval(0);
	    inv11->Fill(invSlopeNew);
	    inv12->Fill(theta,invSlopeNew);
	    inv13->Fill(theta,invSlopeNew);
	  }	
	} // iy < Maxy
      } // iz < Maxz
    }// sec < 2
  
    inv11->Fit("gaus");
    double MeanInv11;
    double SigInv11;

    MeanInv11 = gaus->GetParameter(1);
    SigInv11 = gaus->GetParameter(2);


    // And loop again for CF
    ofstream   outCF1;
    outCF1.open(fCFang_corr,ios::out);
    
    for(int sec = 0;sec<2;sec++){
      for(int iy=0; iy<kMaxY; iy++){
	for(int iz=0; iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin); 
	  
	  float invSlopeNew;
	  float CFInvAngCorr = 1;
	  float cf_err = 0;
	  int status = 1;
	  int arm = 1;
	  
	  if(SLOPES[sec][iy][iz]!=0){
	    invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf->Eval(theta)*fitf->Eval(0);
	    CFInvAngCorr = MeanInv11/invSlopeNew;
	    if(CFInvAngCorr<=fMaxCorr){
	  
	      inv110->Fill(CFInvAngCorr);
	    }	
	    else{
	      status = 0;
	      CFInvAngCorr = 1;
	    }
	  }
	  else{
	    status = 0;
	  }	  
	  outCF1 <<MakeSwkey(arm,sec,iy,iz) << " " << CFInvAngCorr 
		 << " " << cf_err << " " << status <<  endl;      
	} // iy < Maxy
      } // iz < Maxz     
    }// sec < 2
  } // in1.good
  
  in1.close();
  outCF1.close();

}

void angular_all_special(const char* inSlopes){
  // Start again with angiular correction
  // 1. Get the Angular dependence of inverse Slopes  


  GetAngles();

  int statusfail = 0;
  
  const Int_t nba = 25;
  Double_t xla = 0.;
  Double_t xua = 25.;
  Double_t wba;
  wba = (xua - xla)/static_cast<Int_t>(nba);
 
 
  // Fit histo 1
  Double_t min=0.;
  Double_t max=22.;
  // TF1 *fitf[2] = {2*0}; 
  TF1 *fitf = new TF1("f1",fitf,min,max,3);





  // Book all that's needed for the inverse Slopes
  
  TH1D *inv1 = 0;
  TH2D *inv2 = 0;
  TProfile *inv3 = 0;

  TH1D *inv11 = 0;
  TH2D *inv12 = 0;
  TProfile *inv13 = 0;

  TH1D *inv110 = 0;




  inv1 = new TH1D("inv1E0E1","Inverse slopes",500,-1.,0.);
  inv2 = new TH2D("inv2E0E1","Inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv3 = new TProfile("inv3E0E1","Inv Slope vs angle",nba,xla,xua,""); 
  inv11 = new TH1D("inv11E0E1","corr inverse slopes",500,-1.,0.);
  inv12 = new TH2D("inv12E0E1"," corr inverse slopes vs angle",nba,xla,xua,500,-1.,0.);
  inv13 = new TProfile("inv13E0E1","corr inv Slope vs angle",nba,xla,xua,""); 
  inv110 = new TH1D("inv110E0E1","CF from ang Corr inverse Slopes",400,0.,4.);
 
  TH1F *htemp = new TH1F("htemp","",200,-50,0);

  // Open the Slopes File 




  
  ifstream in1;  
  in1.open(inSlopes,ios::in);
  
  
  int swkey;
  float value;

  if(in1.good()){
    cout << "Opening " << inSlopes << endl;
    while(in1 >> swkey >> value){
     
      int iy,iz,arm,sec;

      DecodeSwkey(swkey,arm,sec,iy,iz);

      SLOPES[sec][iy][iz] = value;
            
      if(value!=0&&value>-1000){
	
	int nbin = angleE[sec]->FindBin(iz,iy);
	float theta = angleE[sec]->GetBinContent(nbin);
	if(theta>0&&theta<2){
	  htemp->Fill(value);
	}

	inv1->Fill(1./value);
	inv2->Fill(theta,1./value);
	inv3->Fill(theta,1./value);
       
      }
    } // kArraySize

    fitf->SetParameters(-0.1,1,1);
   
    

    inv3->Draw();
   
    inv3->Fit("f1","R");
   
   

    // Correct the inv Slopes for angular Dependence
  
    for(int sec = 0;sec < 2;sec++){
      for(int iy=0;iy<kMaxY;iy++){
	for(int iz=0;iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin);
	  
	  // cout << sec << " " << iy << " " << iz << endl;
	  if(SLOPES[sec][iy][iz]!=0&&SLOPES[sec][iy][iz]>-1000){
	    
	  
	    float invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf->Eval(theta)*fitf->Eval(0);
	    inv11->Fill(invSlopeNew);
	    inv12->Fill(theta,invSlopeNew);
	    inv13->Fill(theta,invSlopeNew);
	  }	
	} // iy < Maxy
      } // iz < Maxz
    }// sec < 2
  
    inv11->Fit("gaus");
    double MeanInv11;
    double SigInv11;

    MeanInv11 = gaus->GetParameter(1);
    SigInv11 = gaus->GetParameter(2);


    // And loop again for CF
    ofstream   outCF1;
    outCF1.open(fCFang_corr,ios::out);
    
    for(int sec = 0;sec<2;sec++){
      for(int iy=0; iy<kMaxY; iy++){
	for(int iz=0; iz<kMaxZ;iz++){ 
	  
	  int nbin = angleE[sec]->FindBin(iz,iy);
	  float theta = angleE[sec]->GetBinContent(nbin); 
	  
	  float invSlopeNew;
	  float CFInvAngCorr = 1;
	  float cf_err = 0;
	  int status = 1;
	  int arm = 1;
	  
	  if(SLOPES[sec][iy][iz]!=0){
	    invSlopeNew = 1./SLOPES[sec][iy][iz]/fitf->Eval(theta)*fitf->Eval(0);
	    CFInvAngCorr = MeanInv11/invSlopeNew;

	    if(CFInvAngCorr<=fMaxCorr&&TMath::fabs(MeanInv11-invSlopeNew)<5.*SigInv11){
	     
	      inv110->Fill(CFInvAngCorr);
	    }	
	    else{
	      status = 0;	      
	      CFInvAngCorr = 1;
	      statusfail++;
	    }
	  }
	  else{
	    status = 0;
	  }	  
	  outCF1 <<MakeSwkey(arm,sec,iy,iz) << " " << CFInvAngCorr 
		 << " " << cf_err << " " << status <<  endl;      
	} // iy < Maxy
      } // iz < Maxz     
    }// sec < 2
  } // in1.good
  
  cout << statusfail << " with failed status " << endl;

  in1.close();
  outCF1.close();

}


void GetAngles(){

  gSystem->Load("libemcOM.so"); // database access
  gSystem->Load("libemc.so"); // needed for geometry module


  angleE[0] = new TH2F("angleE0","theta from vtx(0,0,0) to module y,z",
		       kMaxZ,-0.5,(1.*kMaxZ)-0.5,kMaxY,-0.5,(1.*kMaxY)-0.5);
  angleE[1] = new TH2F("angleE1","theta from vtx(0,0,0) to module y,z",
		       kMaxZ,-0.5,(1.*kMaxZ)-0.5,kMaxY,-0.5,(1.*kMaxY)-0.5);

  mEmcGeometryModule *mEmcGeometry = new mEmcGeometryModule();
  

 


  for(int is = 6;is<8;is++){
    for(int iy = 0;iy<kMaxY;iy++){
      for(int iz = 0;iz<kMaxZ;iz++){
	int index = iz + iy * kMaxZ;
	
	// local tower pos
	float pos[3];
	mEmcGeometry->GetTowerPosLocal(is,index,pos[0],pos[1],pos[2]);
	
	// 
	float Vert[3];
	mEmcGeometry->GlobalToLocal(0.,0.,0.,is,Vert[0],Vert[1],Vert[2]);
	
	float vx = pos[0] - Vert[0];
	float vy = pos[1] - Vert[1];
	float vz = -Vert[2];
	
	// From this point X coord in sector frame is Z coord in Global Coord System !!!
	double sinT = sqrt((vx*vx+vy*vy)/(vx*vx+vy*vy+vz*vz));
	float theta = TMath::ASin(sinT)*180/TMath::Pi();
	if(is==6){ // E1
	  angleE[1]->Fill(iz,iy,theta);
	}
	if(is==7){ // E0
	  angleE[0]->Fill(iz,iy,theta);
	}

      }
    }
  }
  return;

} 




void CompareSlopes(char *file1,char *file2){



  float cf_1[kSec][kMaxY][kMaxZ] = {kSec*kMaxZ*kMaxY*0.};
  float cf_2[kSec][kMaxY][kMaxZ] = {kSec*kMaxZ*kMaxY*0.};


  ifstream in1;
  in1.open(file1);
  int iy,iz,isec,iarm;
  float value;
  int swkey;
  int status = 0;
  float error;


  if(!in1.good())cout << "Can't open " << file1 << endl;
  while(in1>>swkey>>value){

    DecodeSwkey(swkey,iarm,isec,iy,iz);    
    cf_1[isec][iy][iz] = value;
    if(value<-100)cf_1[isec][iy][iz]=0;

  }
  in1.close();
  
  in1.open(file2);

  if(!in1.good())cout << "Can't open " << file2 << endl;

  while(in1>>swkey>>value){
    DecodeSwkey(swkey,iarm,isec,iy,iz);  
    cf_2[isec][iy][iz] = value;
    if(value<-100)cf_2[isec][iy][iz]=0;    
  }
  in1.close();


  TH1F *hDiff = new TH1F("hDiff","CF1 - CF2",100,-2.,2);
  TH1F *hRatio = new TH1F("hRatio","CF1/CF2",500,0.,2.5);

  TH2F *hDiffSec[2];
  hDiffSec[0] = new TH2F("hDiffE0","CF1 - CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);
  hDiffSec[1] = new TH2F("hDiffE1"," CF1 - CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);

  TH2F *hRatSec[2];
  hRatSec[0] = new TH2F("hRatE0","CF1 /CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);
  hRatSec[1] = new TH2F("hRatE1"," CF1 /CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);


  TH2F *hcorr = new TH2F("hcorr","1 vs 2",400,-30,-10.,400,-30.,-10.);

  for(isec = 0;isec < kSec;isec++){
    for(iy = 0;iy < kMaxY;iy++){
      for(iz = 0;iz < kMaxZ;iz++){

	/*
	int iflipsec = (isec+1)%2;
	int iflipy = (iy+10)%kMaxY;
	int iflipz = (iz+20)%kMaxY;
	*/
	
	int iflipsec = isec;
	int iflipy = iy;
	int iflipz = iz;
	

	if( cf_1[isec][iy][iz]!= 0 && cf_2[iflipsec][iflipy][iflipz]!=0){
	  
      

	  int ymax = kMaxY - 15;
	  int ymin = 15;
	  int zmax = kMaxZ - 30;
	  int zmin = 30;

	  //  if(iy<ymax&&iy>ymin&&iz<zmax&&iz>zmin){// consider only the centers
	  float diff = cf_1[isec][iy][iz] - cf_2[iflipsec][iflipy][iflipz];
	  float ratio = cf_1[isec][iy][iz]/cf_2[iflipsec][iflipy][iflipz];
	  hDiff->Fill(diff);
	  hRatio->Fill(ratio);
	  hDiffSec[isec]->Fill(iz,iy,diff);
	  hRatSec[isec]->Fill(iz,iy,ratio);
	  hcorr->Fill(cf_1[isec][iy][iz],cf_2[iflipsec][iflipy][iflipz]);
	  // }
	    
	}
      } 
    } 
  }

}







void CompareCF(char *file1,char *file2){



  float cf_1[kSec][kMaxY][kMaxZ] = {kSec*kMaxZ*kMaxY*0.};
  float cf_2[kSec][kMaxY][kMaxZ] = {kSec*kMaxZ*kMaxY*0.};


  ifstream in1;
  in1.open(file1);
  int iy,iz,isec,iarm;
  float value;
  int swkey;
  int status = 0;
  float error;


  if(!in1.good())cout << "Can't open " << file1 << endl;
  while(in1>>swkey>>value>>error>>status){

    DecodeSwkey(swkey,iarm,isec,iy,iz);    
    cf_1[isec][iy][iz] = value;
    if(status==0||value>2)cf_1[isec][iy][iz]=0;

  }
  in1.close();
  
  in1.open(file2);

  if(!in1.good())cout << "Can't open " << file2 << endl;

  while(in1>>swkey>>value>>error>>status){
    DecodeSwkey(swkey,iarm,isec,iy,iz);  
    cf_2[isec][iy][iz] = value;
    if(status==0||value>2)cf_2[isec][iy][iz]=0;
  }
  in1.close();


  TH1F *hDiff = new TH1F("hDiff","CF1 - CF2",100,-2.,2);
  TH1F *hRatio = new TH1F("hRatio","CF1/CF2",500,0.,2.5);

  TH2F *hDiffSec[2];
  hDiffSec[0] = new TH2F("hDiffE0","CF1 - CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);
  hDiffSec[1] = new TH2F("hDiffE1"," CF1 - CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);

  TH2F *hRatSec[2];
  hRatSec[0] = new TH2F("hRatE0","CF1 /CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);
  hRatSec[1] = new TH2F("hRatE1"," CF1 /CF2",kMaxZ,-0.5,kMaxZ-0.5,
		 kMaxY,-0.5,kMaxY-0.5);


  TH2F *hcorr = new TH2F("hcorr","1 vs 2",200,0,2,200,0,2);
  TH2F *hacorr = new TH2F("hacorr","1 vs 2",200,0,2,200,0,2);

  for(isec = 0;isec < kSec;isec++){
    for(iy = 0;iy < kMaxY;iy++){
      for(iz = 0;iz < kMaxZ;iz++){

	/*
	int iflipsec = (isec+1)%2;
	int iflipy = (iy+10)%kMaxY;
	int iflipz = (iz+20)%kMaxY;
	*/
	
	int iflipsec = isec;
	int iflipy = iy;
	int iflipz = iz;
	
	if(cf_1[isec][iy][iz]!= 0 && cf_1[iflipsec][iflipy][iflipz]!=0){
	  hacorr->Fill(cf_1[isec][iy][iz],1/cf_1[iflipsec][iflipy][iflipz]);
	}

	

	if( cf_1[isec][iy][iz]!= 0 && cf_2[iflipsec][iflipy][iflipz]!=0){
	  
      

	  int ymax = kMaxY - 15;
	  int ymin = 15;
	  int zmax = kMaxZ - 30;
	  int zmin = 30;

	  //  if(iy<ymax&&iy>ymin&&iz<zmax&&iz>zmin){// consider only the centers
	  float diff = cf_1[isec][iy][iz] - cf_2[iflipsec][iflipy][iflipz];
	  float ratio = cf_1[isec][iy][iz]/cf_2[iflipsec][iflipy][iflipz];
	  hDiff->Fill(diff);
	  hRatio->Fill(ratio);
	  hDiffSec[isec]->Fill(iz,iy,diff);
	  hRatSec[isec]->Fill(iz,iy,ratio);
	  hcorr->Fill(cf_1[isec][iy][iz],cf_2[iflipsec][iflipy][iflipz]);
	  
	  // }
	    
	}
      } 
    } 
  }

}





void CombineTwoCF(char *file1,char *file2,int invert1 = 0,float globalfactor = 1.){



  float cf_1[kSec][kMaxY][kMaxZ] = {kSec*kMaxZ*kMaxY*0.};
  float cf_2[kSec][kMaxY][kMaxZ] = {kSec*kMaxZ*kMaxY*0.};


  ifstream in1;
  in1.open(file1);
  int iy,iz,isec,iarm;
  float value;
  int swkey;
  int status = 0;
  float error;


  TH1F *hNew = new TH1F("hNew"," new cf ",200,0.,4.);

  if(!in1.good())cout << "Can't open " << file1 << endl;
  while(in1>>swkey>>value>>error>>status){

    DecodeSwkey(swkey,iarm,isec,iy,iz);    
    cf_1[isec][iy][iz] = value;
    if(status==0)cf_1[isec][iy][iz]=1;

  }
  in1.close();
  in1.open(file2);

  ofstream out1("xCF_comb.dat");

  if(!in1.good())cout << "Can't open " << file2 << endl;

  while(in1>>swkey>>value>>error>>status){
    DecodeSwkey(swkey,iarm,isec,iy,iz);  
    cf_2[isec][iy][iz] = value;
    if(status==0)cf_2[isec][iy][iz]=1;

    if(invert1){
      if(cf_1[isec][iy][iz]!=0){

	float cf = cf_2[isec][iy][iz]/cf_1[isec][iy][iz]*globalfactor;
	out1 << swkey << " "
	     << cf << " " 
	     << error << " " 
	     << status << endl;
	hNew->Fill(cf);
      }
      else{
	float cf = cf_2[isec][iy][iz]*globalfactor;
	out1 << swkey << " "
	     << cf << " " 
	     << error << " " 
	     << status << endl;
	hNew->Fill(cf);
      }
    }
    else{ 
      float cf = cf_2[isec][iy][iz]*cf_1[isec][iy][iz]*globalfactor;
	out1 << swkey << " "
	     << cf << " " 
	     << error << " " 
	     << status << endl;
	hNew->Fill(cf);
    }

  }
  in1.close();
  out1.close();

}