{


  TFile * f1 = new TFile("PIDFile.root");
  TH2F * hPID = (TH2F*) f1->Get("ipid");

  gStyle->SetOptStat("");  /// no state
  gStyle->SetOptTitle(0);  /// no title

  TCanvas * ccc = new TCanvas("ccc", "ccc", 1000, 1000);

  const int nBinX = hPID->GetNbinsX();
  const int nBinY = hPID->GetNbinsY();
  const double xBinWidth = hPID->GetXaxis()->GetBinWidth(1);
  const double yBinWidth = hPID->GetYaxis()->GetBinWidth(1);

  
  double skewPar[6] = {10.9179, 
                       0.00416034,
                       -233.306,
                       3607.16,
                       1538.37,
                       0.000430609};
                       
  const double c = 299.792458; /// mm/ns
  const double FL = 56423 ;  /// mm 
  const double Brho = 9.04; /// T.mm   
  const double ma = 931.5;
  const double beta0 = Brho * c / TMath::Sqrt( 9*ma*ma + Brho*Brho*c*c); // for AoQ = 3
  const double tofOffset = 504.69; /// ns
  
  
  TH2F * hPID2 = new TH2F("hPID2", "PID; ToF; dE", 500, -280, -210, 500, 0, 5500); 
  TH2F * hZAoQ = new TH2F("hZAoQ", "PID; A/Q; Z", 400, 2.0, 4.0, 400, 1, 17); 
  TH1F * hZ    = new TH1F("hZ", "Z; Z; count", 500, 1, 17); 


  for( int i = 1; i < nBinX ; i++){
    for( int j = 1; j < nBinY ; j++){
      double xpos = hPID->GetXaxis()->GetBinCenter(i) + (gRandom->Rndm()-0.5) * xBinWidth;
      double ypos = hPID->GetYaxis()->GetBinCenter(j) + (gRandom->Rndm()-0.5) * yBinWidth;
      
      double zpos = hPID->GetBinContent(i,j);
      
      
      ///===== Skew correction
      double ToF = xpos - (skewPar[0] * exp( - skewPar[1] * ypos ) + skewPar[2] - skewPar[5] * (skewPar[3] - ypos) * exp( - pow( (skewPar[3]-ypos)/skewPar[4],2) ) ) - 234 + 12; 
      
      
      hPID2->Fill( ToF, ypos, zpos);
      
      double beta = FL/c/(ToF + tofOffset);
      double gamma = 1./TMath::Sqrt(1.-beta*beta);
      
      double Z = sqrt((ypos + 11.9473) / 23.097) * TMath::Power(beta / beta0, 1.3)  ;
      double AoQ = c * Brho / gamma/ beta / ma;

      Z = -0.732942 + 1.47488 * Z  - 0.100088 * Z*Z  + 0.009867 * Z*Z*Z - 0.000476 * Z*Z*Z*Z + 8.90519e-6 * Z*Z*Z*Z*Z;
    
      hZAoQ->Fill(AoQ, Z, zpos);
      hZ->Fill(Z, zpos);
      
    }
  }

  
  double gamma = 4.0;
  int start = 0;
  int color = kColorPrintableOnGrey;  
  gStyle->SetPalette(color);TColor::InvertPalette();
  
  
  double zMax = hZAoQ->GetMaximum();
  
  const int nCon = 200;  /// number of contours
  double contours[nCon];
  
  for( int i = 0; i < nCon; i++){
    double xpos = i*1.0/nCon;
    double ypos = TMath::Power(xpos, gamma) *0.95;
    contours[i] = start + ypos*(zMax - start);
  }

  ///ccc->Divide(2,1);
  
  ccc->cd(1); ///ccc->cd(1)->SetLogz();
  hZAoQ->SetContour(nCon, contours);
  hZAoQ->Draw("col");
  
  TEllipse * si42 = new TEllipse(3.01, 14, 0.05, 0.55);
  si42->SetFillStyle(0);
  si42->SetLineColor(2);
  si42->SetLineWidth(3);
  
  si42->Draw("same");
  
  TLatex text;
  text.SetTextColor(2);
  text.SetTextSize(0.05);
  text.DrawLatex(3.05, 14.2, "^{42}Si");

  //ccc->cd(2);
  //hZ->Draw("hist");
  
  ccc->SaveAs(Form("iPID_%3.1f_%d_%d.pdf", gamma, start, color));



}