using namespace std;

void ask()

{
  double pz_bins[] = {0.30, 0.50, 0.62, 0.75, 0.88, 1.00,1.20,1.50, 2.0, 4.00,6.00,8.00,10.00,12.00,14.00,17.00,20.00,24.00,28.00,34.00,40.00,48.00,56.00,68.00,80.00};
  unsigned long n_pz_bins = sizeof(pz_bins) / sizeof(double) - 1;
  double pT_bins[] = {0.00, 0.10, 0.20, 0.30, 0.40, 0.50, 2.0};
  unsigned long n_pT_bins = sizeof(pT_bins) / sizeof(double) - 1;
 auto c2 = new TCanvas();
 
  TH2D *he = new TH2D("he"," ",
                     n_pz_bins, pz_bins, n_pT_bins, pT_bins);  
  TH2D *he_a = new TH2D("he_a"," ",
                     n_pz_bins, pz_bins, n_pT_bins, pT_bins);  
  TH2D *he_b = new TH2D("he_b"," ",
                     n_pz_bins, pz_bins, n_pT_bins, pT_bins);  
  TH2D *he_c = new TH2D("he_c"," ",
                     n_pz_bins, pz_bins, n_pT_bins, pT_bins);                                                           
                     
     double stats_error[] = {0.97, 1.05, 1.46, 1.97, 0.0, 0.0,
1.22, 1.18, 1.57, 2.12, 5.39, 0.0,
1.41, 1.33, 1.59, 2.09, 3.12, 0.0,
1.74, 1.48, 1.67, 2.19, 2.84, 0.0,
1.87, 1.69, 1.85, 2.87, 2.90, 0.0,
2.04, 1.83, 2.06, 2.46, 4.18, 0.0,
2.35, 2.07, 2.20, 2.51, 3.11, 0.0,
2.90, 2.44, 2.41, 2.70, 3.14, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0,        //missing data [1.5,2.00)
11.19, 7.63, 8.72, 9.31, 0.0, 0.0,
4.78, 3.55, 3.82, 3.88, 7.89, 0.0,
3.55, 6.28, 1.70, 1.86, 2.64, 5.87,
4.01, 1.95, 1.72, 2.10, 2.79, 4.18,
4.98, 2.43, 2.00, 2.14, 2.40, 2.96,
4.62, 2.50, 1.88, 1.87, 2.12, 2.15,
5.78, 2.99, 2.17, 2.45, 2.55, 1.88,
5.79, 3.22, 2.98, 2.52, 2.72, 2.02,
6.65, 5.74, 3.92, 3.49, 3.35, 2.05,
11.39, 5.12, 4.06, 3.45, 3.36, 2.13,
8.56, 5.38, 4.20, 3.51, 3.44, 1.88,
9.82, 6.32, 4.63, 4.11, 3.83, 2.01,
11.91, 10.33, 5.87, 5.31, 5.21, 2.46,
11.86, 9.19, 6.28, 5.64, 6.22, 2.77, 
0.0, 0.0, 10.32, 10.39, 10.15, 4.69};

unsigned long n_stats_error = sizeof(stats_error)/sizeof(double);
 cout << "Number of errors data points :" <<n_stats_error<<endl;   
 
 c2->cd();   
 he->SetStats(kFALSE);
 
 for ( int k = 1; k<=n_pz_bins; k++)
    for ( int l = 1; l <= n_pT_bins; l++)
     he->SetBinContent(k, l, stats_error[((k - 1) * n_pT_bins + (l - 1))]); 
            
     gStyle->SetPalette(70);  
     he->GetXaxis()->SetRangeUser(0.0,68.0);
     he->GetYaxis()->SetRangeUser(0.0,0.8);  
     he->Draw("COLZ");
     
 for ( int ka = 1; ka<=n_pz_bins; ka++)
    for ( int la = 1; la <= n_pT_bins; la++)
     he_a->SetBinContent(ka ,la , stats_error[((ka - 1) * n_pT_bins + (la - 1))] > 0.0 && stats_error[((ka - 1) * n_pT_bins + (la - 1))] < 2.5 );     
     
    he_a->SetMarkerSize(0.8); 
    he_a->SetMarkerStyle(kFullDotLarge); 
    he_a->Draw("SAME");  
    
  for ( int kb = 1; kb<=n_pz_bins; kb++) 
    for ( int lb = 1; lb <= n_pT_bins; lb++)
     he_b->SetBinContent(kb, lb, stats_error[((kb - 1) * n_pT_bins + (lb - 1))]> 2.5 && stats_error[((kb - 1) * n_pT_bins + (lb - 1))] < 5.0 );     
     
    he_b->SetMarkerStyle(kCircle); 
    he_b->SetMarkerSize(0.6);
    he_b->Draw("SAME"); 
    
  for ( int kc = 1; kc<=n_pz_bins; kc++) 
    for ( int lc = 1; lc <= n_pT_bins; lc++)
    he_c->SetBinContent(kc, lc, stats_error[((kc - 1) * n_pT_bins + (lc - 1))] > 5.0 );     
     
    he_c->SetMarkerStyle(kPlus); 
    he_c->SetMarkerSize(0.7);
    he_c->Draw("SAME");   
   c2->Update();
                      
  
}