{
//////////////////////////////////////////////////////
//
// RESIDUAL_ANALYSIS.C:
//
// Root macro performing ACDC runs residuals analysis.
// For monitoring purposes.
//
// This macro uses the sample of 1000 tracks extracted 
// from each alignment run. These tracks are not used for 
// the alignment, and are fitted before and after it, so
// that residuals could be compared
//
// S.V. 25/07/06
//
//////////////////////////////////////////////////////


gROOT->Reset();

//
// Macro options
//

// Which side you want to check (1 == RIGHT / 0 == LEFT) ?
int side = 1; 

// Number of modules aligned in this side
const int nstations = 21; 

// Which modules you want to check ?
int n_modules = 6;
int modules[] = {10, 11, 12, 13, 14, 15};

// In how many R and Phi slices do you want to check your residuals ? 
const int    n_slices   = 6;  

// Min. number of hit per slice (in order to get relevant graphs)
int n_hitmin = 1; 

//
// Do you want to run on one file or on a bunch of file ?
// 
// --> If you choose one file you will have classic residual plots
// --> If you choose many files you will have sumary plots (mean res. before and after alignment) 
//

bool multifiles = false; 

int first_file = 1;
int last_file  = 1;



//
// Initialization
//

char    data_file[80];
char    canvas_title[100], canvas_name[3], module_num[30];

int  t_step, t_side; 

float   t_vx, t_vy, t_vz, t_station, t_resX, t_resY, t_X, t_Y;

int     N_polar[4*nstations][n_slices];

float  Coord_polar[4*nstations][n_slices];
float  Residual_polar[4*nstations][n_slices];
float  RMS_polar[4*nstations][n_slices];
float  RMS_polar2[4*nstations][n_slices];

TH1F *R_mean_before = new TH1F("mean_R_residual_before","",60,-0.1,0.1);
TH1F *R_mean_after  = new TH1F("mean_R_residual_after","",60,-0.03,0.03);
TH1F *P_mean_before = new TH1F("mean_Phi_residual_before","",60,-0.02,0.02);
TH1F *P_mean_after  = new TH1F("mean_Phi_residual_after","",60,-0.005,0.005);

float Pi = 4.*atan(1.);
float div_R = (45. - 8.)/n_slices; 
float div_P = Pi/n_slices; 

float z_right[] = {-158.85, -131.15, -98.85, -71.15, -38.85, -11.15,
		    21.15, 48.85, 81.15, 108.85, 141.15, 168.85,
		    201.15, 228.85, 261.15, 288.85, 451.15, 598.85,
		    651.15, 698.85, 751.15};

float z_left[] = {-176.15, -143.85, -116.15, -83.85, -56.15, -23.85,
		   3.85, 36.15, 63.85, 96.15, 123.85, 156.15,
		   183.85, 216.15, 243.85, 276.15, 433.85, 586.15,
		   633.85, 686.15, 733.85};

float z_module[21];

for (int i=0; i<21; ++i) z_module[i] = 0.;

for (int i=0; i<n_modules; ++i) 
{ 
  if (modules[i] < 0 || modules[i] >= nstations) continue;

  (side == 0)
    ? z_module[modules[i]] = z_left[modules[i]]
    : z_module[modules[i]] = z_right[modules[i]];
}

//
// Then get the result file (or loop on the requested set) 
//

for (int k=first_file; k<=last_file; ++k)
{
  cout << "Processing run " << k << endl; 

  if (k == 666) continue; // suppress the files you don't want

  sprintf(data_file,  "Align.root",k);  // Result file

  TFile *res = new TFile(data_file);
  TTree *check = (TTree*)res->Get("Align/Trackinfo");

  check->SetBranchAddress("step"    ,&t_step);
  check->SetBranchAddress("side"    ,&t_side);
  check->SetBranchAddress("vx"      ,&t_vx);
  check->SetBranchAddress("vy"      ,&t_vy);
  check->SetBranchAddress("vz"      ,&t_vz);
  check->SetBranchAddress("station" ,&t_station);
  check->SetBranchAddress("X_clus"  ,&t_X);
  check->SetBranchAddress("Y_clus"  ,&t_Y);
  check->SetBranchAddress("res_X"   ,&t_resX);
  check->SetBranchAddress("res_Y"   ,&t_resY);

  int nentries   = check->GetEntries();

  for (int i=0; i<4*nstations; ++i) 
  { 
    for (int j=0; j<n_slices; ++j) 
    { 
      N_polar[i][j]        = 0;
      Coord_polar[i][j]    = 0.;
      Residual_polar[i][j] = 0.;
      RMS_polar[i][j]      = 0.;
      RMS_polar2[i][j]     = 0.;
    }
  }


  // Nomenclature of the rows:
  //
  // i = 4*nstat + j with j = 0,1,2,3 
  //
  //  j = 0 : initial R   
  //  j = 1 : initial Phi 
  //  j = 2 : final R
  //  j = 3 : final Phi

  // First loop for means  
  
  for (Int_t i=0; i<nentries; i++) 
  {       
    check->GetEntry(i);

    if (t_side != side) continue;
    
    int nstat = -1;

    for (int j=0; j<nstations; ++j) // Get the station number
    {
      if (z_module[j] == t_station) nstat = j;
    }    

    if (nstat == -1) continue;

    int nrank = 4*nstat+2*t_step;  // Which row will we fill

    // Compute R/Phi, and the corresponding residuals

    float radius = sqrt((t_X-t_resX)*(t_X-t_resX)+(t_Y-t_resY)*(t_Y-t_resY));
    float phi    = atan2(t_X-t_resX,t_Y-t_resY);

    float resiR  = sqrt(t_X*t_X+t_Y*t_Y) - radius;
    float resiP  = atan2(t_X,t_Y) - phi;


    cout << radius << " / " << resiR << endl;


    for (int j=0; j<n_slices; ++j) // Fill the correct slice
    {  
      if (radius > 8.+div_R*j && radius < 8.+div_R*(j+1))
      {
	N_polar[nrank][j]++;
	Coord_polar[nrank][j] += radius;
	Residual_polar[nrank][j] += resiR;
      }

      if (phi > div_P*j && phi < div_P*(j+1))
      {
	N_polar[nrank+1][j]++;
	Coord_polar[nrank+1][j] += phi;
	Residual_polar[nrank+1][j] += resiP;
      }
    }
  }

  // Then we got the means

  for (int i=0; i<4*nstations; ++i) 
  { 
    for (int j=0; j<n_slices; ++j) 
    { 
      if (N_polar[i][j] >= n_hitmin)
      {
	Coord_polar[i][j]    /= N_polar[i][j];
	Residual_polar[i][j] /= N_polar[i][j];

	if (i%4 == 0) R_mean_before->Fill(Residual_polar[i][j]);
	if (i%4 == 1) P_mean_before->Fill(Residual_polar[i][j]);
	if (i%4 == 2) R_mean_after->Fill(Residual_polar[i][j]);
	if (i%4 == 3) P_mean_after->Fill(Residual_polar[i][j]);
      }
      else
      {
	N_polar[i][j]        = 0; 
	Coord_polar[i][j]    = 0.;
	Residual_polar[i][j] = 0.;
      }
    }
  }


  //
  // Then a second loop for the RMS
  //

  if (multifiles) continue;

  for (Int_t i=0; i<nentries; i++) 
  {       
    check->GetEntry(i);
    
    int nstat = -1;

    for (int j=0; j<nstations; ++j)
    {
      if (z_module[j] == t_station) nstat = j;
    }

    if (nstat == -1) continue;
    
    int nrank = 4*nstat+2*t_step;

    float radius = sqrt((t_X-t_resX)*(t_X-t_resX)+(t_Y-t_resY)*(t_Y-t_resY));
    float phi    = atan2(t_X-t_resX,t_Y-t_resY);

    float resiR  = sqrt(t_X*t_X+t_Y*t_Y) - radius;
    float resiP  = atan2(t_X,t_Y) - phi;

    for (int j=0; j<n_slices; ++j) 
    {  
      if (radius > 8.+div_R*j && radius < 8.+div_R*(j+1))
      {
	RMS_polar[nrank][j]  += (resiR-Residual_polar[nrank][j])**2;
	RMS_polar2[nrank][j] += (radius-Coord_polar[nrank][j])**2;
      }

      if (phi > div_P*j && phi < div_P*(j+1))
      {
	RMS_polar[nrank+1][j]  += (resiP-Residual_polar[nrank+1][j])**2;
	RMS_polar2[nrank+1][j] += (phi-Coord_polar[nrank+1][j])**2;
      }
    }
  }
  
  // Then we got the RMSs

  for (int i=0; i<4*nstations; ++i) 
  { 
    for (int j=0; j<n_slices; ++j) 
    { 
      if (N_polar[i][j] >= n_hitmin)
      {
	RMS_polar[i][j] /= N_polar[i][j];
	RMS_polar[i][j] = sqrt(RMS_polar[i][j]);
	RMS_polar2[i][j] /= N_polar[i][j];
	RMS_polar2[i][j] = sqrt(RMS_polar2[i][j]);
      }
    }
  }
}


//
///
//// Then finally make some cosmetics
///
//

/////////////STAT////////   
gStyle -> SetStatW(0.3);
gStyle -> SetStatH(0.2);
gStyle -> SetStatColor(41);   
gStyle -> SetStatX(0.8);
gStyle -> SetStatY(0.8);   
gStyle -> SetStatFont(42);
gStyle -> SetOptStat(0);
gStyle -> SetOptFit(111);
gStyle -> SetPalette(45);
/////////////////////////

if (!multifiles)
{

  TH2F *GraphRi = new TH2F("GraphRi","",300,5.,50.,100,-0.05,0.05);
  TH2F *GraphRf = new TH2F("GraphRf","",300,5.,50.,100,-0.05,0.05);
  TH2F *GraphPi = new TH2F("GraphPi","",300,0.,3.5,100,-0.005,0.005);
  TH2F *GraphPf = new TH2F("GraphPf","",300,0.,3.5,100,-0.005,0.005);
  
  TGraphErrors *r_before = new TGraphErrors(n_slices,Coord_polar[0],Residual_polar[0],RMS_polar2[0],RMS_polar[0]);
  TGraphErrors *p_before = new TGraphErrors(n_slices,Coord_polar[0],Residual_polar[0],RMS_polar2[0],RMS_polar[0]);
  TGraphErrors *r_after  = new TGraphErrors(n_slices,Coord_polar[0],Residual_polar[0],RMS_polar2[0],RMS_polar[0]);
  TGraphErrors *p_after  = new TGraphErrors(n_slices,Coord_polar[0],Residual_polar[0],RMS_polar2[0],RMS_polar[0]);
  TPaveText *pt = new TPaveText(35.,0.025,45.,0.04,"br");
  
  for (int i=0; i<nstations; ++i) 
  { 
    if (z_module[i] == 0.) continue;

    sprintf(canvas_title, "Track residual values for module %d",i); 
    sprintf(canvas_name,  "c%d",i); 
    sprintf(module_num,  "Module %d",i); 
 
    r_before = new TGraphErrors(n_slices,Coord_polar[4*i],Residual_polar[4*i],RMS_polar2[4*i],RMS_polar[4*i]);
    p_before = new TGraphErrors(n_slices,Coord_polar[4*i+1],Residual_polar[4*i+1],RMS_polar2[4*i+1],RMS_polar[4*i+1]);
    r_after  = new TGraphErrors(n_slices,Coord_polar[4*i+2],Residual_polar[4*i+2],RMS_polar2[4*i+2],RMS_polar[4*i+2]);
    p_after  = new TGraphErrors(n_slices,Coord_polar[4*i+3],Residual_polar[4*i+3],RMS_polar2[4*i+3],RMS_polar[4*i+3]);

    c1 = new TCanvas(canvas_name,canvas_title,200,500,700,800);
    c1->SetFillColor(0);
    c1->SetBorderMode(0); 
    c1->Divide(2,2);  

    c1->cd(1);
    r_before->SetMarkerStyle(21); 
    r_before->SetMarkerColor(1);
    GraphRi->GetXaxis()->SetTitle("Hit radius (in mm)");
    GraphRi->GetYaxis()->SetTitleSize(0.04);
    GraphRi->GetYaxis()->SetTitleOffset(1.5);
    GraphRi->GetYaxis()->SetTitle("Residual before alignment (in mm)");
    GraphRi->Draw();
    pt = new TPaveText(32.,0.028,48.,0.04,"br");
    pt->SetTextSize(0.05);
    TText *text = pt->AddText(module_num);
    pt->Draw();
    r_before->Draw("P");

    c1->cd(2);
    r_after->SetMarkerStyle(21); 
    r_after->SetMarkerColor(1);
    GraphRf->GetXaxis()->SetTitle("Hit radius (in mm)");
    GraphRf->GetYaxis()->SetTitleSize(0.04);
    GraphRf->GetYaxis()->SetTitleOffset(1.5);
    GraphRf->GetYaxis()->SetTitle("Residual after alignment (in mm)");
    GraphRf->Draw();
    r_after->Draw("P");
    
    c1->cd(3);
    p_before->SetMarkerStyle(21); 
    p_before->SetMarkerColor(1);
    GraphPi->GetXaxis()->SetTitle("Hit angle (in rad)");
    GraphPi->GetYaxis()->SetTitleSize(0.04);
    GraphPi->GetYaxis()->SetTitleOffset(1.8);
    GraphPi->GetYaxis()->SetTitle("Residual before alignment (in rad)");
    GraphPi->Draw();
    p_before->Draw("P");
    
    c1->cd(4);
    p_after->SetMarkerStyle(21); 
    p_after->SetMarkerColor(1);
    GraphPf->GetXaxis()->SetTitle("Hit angle (in rad)");
    GraphPf->GetYaxis()->SetTitleSize(0.04);
    GraphPf->GetYaxis()->SetTitleOffset(1.8);
    GraphPf->GetYaxis()->SetTitle("Residual after alignment (in rad)");
    GraphPf->Draw();
    p_after->Draw("P");
    
    c1->Update();
  }
}
else
{

c4 = new TCanvas("c4","",200,500,700,800);
c4->SetFillColor(0);
c4->SetBorderMode(0); 
c4->Divide(2,2);     

c4->cd(1);
R_mean_before->GetXaxis()->SetTitle("Mean X residual bef. alignment (in mm)");
R_mean_before->GetXaxis()->SetTitleSize(0.04);
R_mean_before->Draw();
R_mean_before->Fit("gaus");

c4->cd(2);
R_mean_after->GetXaxis()->SetTitle("Mean Y residual bef. alignment (in mm)");
R_mean_after->GetXaxis()->SetTitleSize(0.04);
R_mean_after->Draw();
R_mean_after->Fit("gaus");

c4->cd(3);
P_mean_before->GetXaxis()->SetTitle("Mean X residual aft. alignment (in mm)");
P_mean_before->GetXaxis()->SetTitleSize(0.04);
P_mean_before->Draw();
P_mean_before->Fit("gaus");

c4->cd(4);
P_mean_after->GetXaxis()->SetTitle("Mean Y residual aft. alignment (in mm)");
P_mean_after->GetXaxis()->SetTitleSize(0.04);
P_mean_after->Draw();
P_mean_after->Fit("gaus");

c4->Update();
}

}