// 14-jan-2003; 17-mar-2003 - return to business
#include "AliBgInfoForJet.h"
#include <TROOT.h>
#include <TBrowser.h>
#include <TList.h>
#include <TH1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TStyle.h>
#include <TArrayD.h>
#include <TVector3.h>
#include <AliEMCALJetMicroDst.h>

ClassImp(AliBgInfoForJet)

Double_t *ar=0;
Char_t   *text[3]={"EMCAL","ChTr","GRID"};
TArrayD  Ra;  // array of jet radius - 26-mar-2003

AliBgInfoForJet::AliBgInfoForJet(const Int_t mode, const bool putToBrowser)
{
   fNp           = 0;
   fMaxThreshold = 10.;
   fC1 = fC2 = fC3 = fC4 = 0;
   fP1 = fP2 = fP3 = fP4 = 0;

   fMode = mode;
   if(putToBrowser) gROOT->GetListOfBrowsables()->Add(this);
   if(mode>0) {
      fMdst = new  AliEMCALJetMicroDst;
      fMdst->Open(mode);
   }

   TString st(fMdst->GetTitle());
   st += " : ";
   st += mode;
   SetName(st.Data());
   Init();
}

AliBgInfoForJet::~AliBgInfoForJet() 
{
}

void AliBgInfoForJet::Init()
{
  Book1();
  Book2();
}

void AliBgInfoForJet::Book1()
{
  // Main characteristics of BG
  gROOT->cd();
  Double_t maxCell=14000., maxCellE=40., maxSumE = 6000.;
  if(fMdst->IsPythiaDst()) {
      maxCell  = 1000.;
      maxCellE = 200.;
      maxSumE  = 400.;
  }
  fHNcellEMCAL  = new TH1F("fHNcellEMCAL", " #cell with Et > 0", 1000, 0.0, maxCell);
  fHEtCellEMCAL = new TH1F("fHEtCellEMCAL", " Et of cell ", 400, 0.0, maxCellE);
  fHSumEMCAL    = new TH1F("fHSumEMCAL", " Sum Et in EMCAL ", 400, 0.0, maxSumE);

  fHNcellGRID   = new TH1F("fHNcellGRID", " #cell with Et > 0 for GRID", 1000, 0.0, maxCell);
  fHEtCellGRID  = new TH1F("fHEtCellGRID", " Et of cell in GRID", 400, 0.0, maxCellE);
  fHSumGRID     = new TH1F("fHSumGRID", " Sum Et in GRID ", 400, 0.0, maxSumE);

  fHNcellCHTR   = new TH1F("fHNcellCHTR", " #cell with Et > 0 for CHTR", 1000, 0.0, maxCell);
  fHEtCellCHTR  = new TH1F("fHEtCellCHTR", " Et of cell in CHTR", 400, 0.0, maxCellE);
  fHSumCHTR     = new TH1F("fHSumCHTR", " Sum Et in CHTR ", 400, 0.0, maxSumE);

  fListHist1 = AliEMCALJetMicroDst::MoveHistsToList("List 1", kFALSE); 
}

void AliBgInfoForJet::Book2()
{
  gROOT->cd();
  DefineArrays();

  fHOccupancyEMCAL = new TH1F("fHOccupancyEMCAL","Occupancy of EMCAL vs cell threshold", 
  fNp, ar);
  fHEtMeanEMCAL     = new TH1F("fHEtMeanEMCAL","#bar{E_{t}} in EMCAL vs cell threshold",
  fNp, ar);
  fHEtRMSEMCAL     = new TH1F("fHEtRMSEMCAL","RMS E_{t} in EMCAL vs cell threshold",
  fNp, ar);
  // 20-mar-2003
  fHncellmore0EMCAL = new TH1F("fHncellmore0EMCAL","#event with cell(s) EMCAL above threshold(%)", 
  fNp, ar);

  // Working hists.
  fHNcellTmp = new TH1F("fHNcellTmp", " #cell with Et > EtCut", 14001, -0.5, 14000.5);
  fHSumTmp   = new TH1F("fHSumTmp", " Sum Et in EMCAL for cell with Et > EtCut", 6000, 0.0, 6000.0);

  // hists for EMCAL, ChTr and grid 20-mar-2003
  Char_t nam[50], tit[100];
  for(Int_t i=0; i<3; i++){
     sprintf(nam,"Occu%s",text[i]);
     sprintf(tit,"Occupancy of %s vs threshold",text[i]);
     fhOccu[i] = new TH1F(nam,tit, fNp, ar);

     sprintf(nam,"EtMean%s",text[i]);
     sprintf(tit,"mean Et in %s vs threshold",text[i]);
     fhEtMean[i] = new TH1F(nam,tit, fNp, ar);

     sprintf(nam,"EtRMS%s",text[i]);
     sprintf(tit,"RMS Et in %s vs threshold",text[i]);
     fhEtRMS[i] = new TH1F(nam,tit, fNp, ar);

     sprintf(nam,"Cellmore0%s",text[i]);
     sprintf(tit,"#event with cell(s) %s above threshold(%)",text[i]);
     fhNcellmore0[i] = new TH1F(nam,tit, fNp, ar);
  }

  if(fMdst->IsPythiaDst()) {
     // Working hist. ; Ra - array for jet cone;
     fHSumTmpJet = new TH1F("fHSumTmpJet","Sum Et in Cone for cell Et>EtCut(in event)", 
     Ra.GetSize()-1, Ra.GetArray());
     fHSumTmp2Jet = new TH2F("fHSumTmp2Jet","Sum Et in Cone for cell Et > EtCut vs event", 
     3000, 0.0, 300., Ra.GetSize()-1, Ra.GetArray());
     fHSumTmp2Jet->Sumw2(); 
     // Final
     for(Int_t i=0; i<3; i++){
        sprintf(nam,"fHSum2Jet%s",text[i]);
        sprintf(tit,"Et in jet cone for %s above threshold",text[i]);
        fHSumJet[i] = new TH2F(nam, tit, 
        fNp,ar, Ra.GetSize()-1, Ra.GetArray());
	//        fHSumJet[i]->Sumw2();
        sprintf(nam,"fHRelResJet%s",text[i]);
        sprintf(tit,"Relative resolution ; %s above threshold",text[i]);
        fHRelResJet[i] = new TH2F(nam, tit, fNp,ar, Ra.GetSize()-1, Ra.GetArray());
     }
  }
  
  fListHist2 = AliEMCALJetMicroDst::MoveHistsToList("List 2", kFALSE); 
}

void AliBgInfoForJet::DefineArrays(Int_t pri)
{
  if(fMdst->IsPythiaDst()) {
     fNp = 20;
     ar  = new Double_t[fNp+1];
     for(Int_t i=0;  i<=5; i++)  ar[i] = 0.2*i;               // step 0.2 in [0,1] 
     for(Int_t i=6;  i<=10; i++) ar[i] = ar[5] + 0.4*(i-5);    // step 0.4 in [1,3] 
     for(Int_t i=11; i<=15; i++) ar[i] = ar[10] + 0.8*(i-10);  // step 0.8 in [3,7] 
     for(Int_t i=16; i<=20; i++) ar[i] = ar[15] + 2.0*(i-15);  // step 2.  in [7,17] 
     ar[20] = 20.;
  } else {
     fNp = 40;
     ar  = new Double_t[fNp+1];
     for(Int_t i=0;  i<=10; i++) ar[i] = 0.1*i;  // step 0.1 in [0,1] 
     for(Int_t i=11; i<=20; i++) ar[i] = 1. + 0.2*(i-10); // step 0.2 in [1,3] 
     for(Int_t i=21; i<=30; i++) ar[i] = 3. + 0.4*(i-20); // step 0.4 in [3,7] 
     for(Int_t i=31; i<=40; i++) ar[i] = 7. + 1.0*(i-30); // step 1.0 in [7,17]
     ar[40] = 20.;
     fNp = 38; // 18-mar-2003
  }

  Double_t r[4]={0.0, 0.3, 0.4, 0.5};
  Ra.Set(4,r);

  if(pri>0) {
     for(Int_t i=0;  i<=fNp; i++) printf(" ar[%3i] %5.1f \n", i,ar[i]);
     for(Int_t i=0;  i<Ra.GetSize(); i++) printf(" Ra[%1i] %5.1f \n",i,Ra[i]);
  }
}

void AliBgInfoForJet::Fill1(AliEMCALJetMicroDst *miDst)
{
  if(miDst == 0) miDst = fMdst;
  if(miDst == 0) {
    printf("<E> AliBgInfoForJet::Fill1 : Defined miDst !!! \n");
    return;
  }

  Double_t sumEt = 0, sumPt = 0, etCut = 0.0;
  TTree *tr = fMdst->fTree;
  
  for(Int_t ev=0; ev<tr->GetEntries(); ev++){       // cycle on events         
     fMdst->GetEntry(ev);
     fHNcellEMCAL->Fill(fMdst->ncell);
     fHNcellCHTR->Fill(fMdst->nchp);
     fHNcellGRID->Fill(fMdst->ngrid);


     sumEt = 0;
     for(Int_t ic=0; ic<fMdst->ncell; ic++){        // cycle on cells(EMCAL)
        fHEtCellEMCAL->Fill(fMdst->etcell[ic]);
        sumEt += fMdst->etcell[ic];
     }
     fHSumEMCAL->Fill(sumEt);

     sumPt = 0;
     for(Int_t ip=0; ip<fMdst->nchp; ip++){        // cycle on charge particles
        fHEtCellCHTR->Fill(fMdst->ppt[ip]);
        sumPt += fMdst->ppt[ip];
     }
     fHSumCHTR->Fill(sumPt);

     sumEt = 0;
     for(Int_t ic=0; ic<fMdst->ngrid; ic++){        // cycle on cells(GRID)
        fHEtCellGRID->Fill(fMdst->etgrid[ic]);
        sumEt += fMdst->etgrid[ic];
     }
     fHSumGRID->Fill(sumEt);
  }
}

void AliBgInfoForJet::Fill2(AliEMCALJetMicroDst *miDst)
{
  printf("<I> AliBgInfoForJet::Fill2 \n");
  if(miDst == 0) miDst = fMdst;
  if(miDst == 0) {
    printf("<E> : Defined miDst !!! \n");
    return;
  }

  Int_t ncell=0;
  Double_t sumEt = 0, etCut = 0.0;
  TTree *tr = fMdst->fTree;
  for(Int_t ip=1; ip<=fNp; ip++){                     // cycle on thresholds
     etCut = ar[ip-1]; 
     fHNcellTmp->Reset();
     fHSumTmp->Reset();
     for(Int_t ev=0; ev<tr->GetEntries(); ev++){       // cycle on events         
        fMdst->GetEntry(ev);
        ncell = 0;
        sumEt = 0;
        for(Int_t ic=0; ic<fMdst->ncell; ic++){        // cycle on cells
	  if(fMdst->etcell[ic]>etCut) {
	     ncell++;
             sumEt += fMdst->etcell[ic];
          }
        }
        fHNcellTmp->Fill(Double_t(ncell));
        fHSumTmp->Fill(sumEt);
	//        printf("%3i ev : ncell %4i sumEt %f\n", ev, ncell, sumEt); 
     }
     //     Int_t N = Int_t(fHSumTmp->Integral()), Nmore0 = Int_t(fHSumTmp->Integral(2, fHSumTmp->GetNbinsX()));
     Int_t N = Int_t(fHNcellTmp->Integral()), Nmore0 = Int_t(fHNcellTmp->Integral(2, fHNcellTmp->GetNbinsX()));
     if(Nmore0 == 0) break;  // Threshold is high already 

     fHOccupancyEMCAL->SetBinContent(ip, fHNcellTmp->GetMean());
     fHOccupancyEMCAL->SetBinError(ip, fHNcellTmp->GetRMS()); 
     fHEtMeanEMCAL->SetBinContent(ip, fHSumTmp->GetMean()); 
     fHEtMeanEMCAL->SetBinError(ip, fHSumTmp->GetRMS());
     fHEtRMSEMCAL->SetBinContent(ip, fHSumTmp->GetRMS());
     // 20-mar-2003
     Double_t prob, errProb;
     TAliasPAI::probAndErr(Nmore0, N, prob, errProb, 0);     
     fHncellmore0EMCAL->SetBinContent(ip, 100.*prob);
     fHncellmore0EMCAL->SetBinError(ip, 100.*errProb);
     printf(" %4i threshold %6.2f : N %3i : N(ncell>0) %3i : prob %4.1f+/-%4.2f%%\n", 
     ip-1, etCut, N, Nmore0, 100.*prob, 100.*errProb);
  }
}

void AliBgInfoForJet::Fill3(AliEMCALJetMicroDst *miDst)
{
  printf("<I> AliBgInfoForJet::Fill3 \n");
  if(miDst == 0) miDst = fMdst;
  if(miDst == 0) {
    printf("<E> : Defined miDst !!! \n");
    return;
  }

  Int_t ncell;
  Double_t sumEt=0.0, etCut = 0.0, prob=0.0, errProb=0.0, rc=0;
  TTree *tr = fMdst->fTree;
  TVector3 jet, cell;
  //  fNp = 1;
  for(Int_t ip=1; ip<=fNp; ip++){                      // cycle on thresholds
     etCut = ar[ip-1];

     for(Int_t i=0; i<3; i++){                         // cycle on emcal, ch.tr or grid 
        fHNcellTmp->Reset();
        fHSumTmp->Reset();
        if(fMdst->IsPythiaDst()) fHSumTmp2Jet->Reset();
        for(Int_t ev=0; ev<tr->GetEntries(); ev++){       // cycle on events         
           fMdst->GetEntry(ev);
           if(fMdst->IsPythiaDst() && fMdst->njet!=1) continue; // 1-apr-2003
           ncell = 0;
           sumEt = 0;
           if(fMdst->IsPythiaDst() && fMdst->njet==1) {
              jet.SetPtEtaPhi(fMdst->jet[0],fMdst->jetaw[0],fMdst->jphiw[0]);
              fHSumTmpJet->Reset();
           }

           Int_t nc;
           Float_t *etc=0, *etac=0, *phic=0;
           if(i==0) {
              nc   = fMdst->ncell;
              etc  = fMdst->etcell; 
           } else if(i==1) {
              nc   = fMdst->nchp;
              etc  = fMdst->ppt;
              etac = fMdst->peta; 
              phic = fMdst->pphi; 
           } else if(i==2) {
              nc  = fMdst->ngrid;
              etc = fMdst->etgrid; 
	   }

           for(Int_t ic=0; ic<nc; ic++){        // cycle on cells
	      if(etc[ic]>etCut) {
	         ncell++;
                 sumEt += etc[ic];
                 if(fMdst->IsPythiaDst() && fMdst->njet==1) {
                    if(i==0) {
		      cell = fMdst->GetCellVector(ic);
                    } else if(i==1) {
                       cell.SetPtEtaPhi(etc[ic],etac[ic],phic[ic]);
                    } else if(i==2) {
		       cell = fMdst->GetGridVector(ic);
                    }
                    rc = jet.DeltaR(cell);
                    fHSumTmpJet->Fill(rc, etc[ic]);
	         }
              }
           }                                    // end cycle on cells

           fHNcellTmp->Fill(Double_t(ncell));
           fHSumTmp->Fill(sumEt);
           if(fMdst->IsPythiaDst() && fMdst->njet==1) {
	     for(Int_t ix=1; ix<=fHSumTmpJet->GetNbinsX();ix++){
	        fHSumTmp2Jet->Fill(fHSumTmpJet->Integral(1,ix),fHSumTmpJet->GetBinCenter(ix));
             }
           }

        }  // end cycle on events

        Int_t N = Int_t(fHNcellTmp->Integral()), Nmore0 = Int_t(fHNcellTmp->Integral(2, fHNcellTmp->GetNbinsX()));

        fhOccu[i]->SetBinContent(ip, fHNcellTmp->GetMean());
        fhOccu[i]->SetBinError(ip, fHNcellTmp->GetRMS()); 
        fhEtMean[i]->SetBinContent(ip, fHSumTmp->GetMean()); 
        fhEtMean[i]->SetBinError(ip, fHSumTmp->GetRMS());
        fhEtRMS[i]->SetBinContent(ip, fHSumTmp->GetRMS());

        TAliasPAI::probAndErr(Nmore0, N, prob, errProb, 0);     
        printf("Case %1i %4i threshold %6.2f : N %3i : N(ncell>0) %3i : prob %4.1f+/-%4.2f%%\n", 
        i, ip-1, etCut, N, Nmore0, 100.*prob, 100.*errProb);
        printf(" Mean %f : RMS %f \n", fHSumTmp->GetMean(), fHSumTmp->GetRMS());
        if(Nmore0 > 0) {  // Threshold is high already 
           fhNcellmore0[i]->SetBinContent(ip, 100.*prob);
           fhNcellmore0[i]->SetBinError(ip, 100.*errProb);
        }
        if(fMdst->IsPythiaDst()) {
	   Double_t mean=0., rms=0.;
           TH1D* hid=0;
	   for(Int_t iy=1; iy<=fHSumTmp2Jet->GetNbinsY(); iy++){
	      hid =  fHSumTmp2Jet->ProjectionX("hidTmp",iy,iy);
              mean = hid->GetMean();
              rms  = hid->GetRMS();
	      fHSumJet[i]->SetBinContent(ip, iy, mean);
	      fHSumJet[i]->SetBinError(ip, iy, rms);
              printf(" jet %i : mean %f : rms %f \n", iy, mean, rms);
	      delete hid;
           }
        }
     }
  }
}

void AliBgInfoForJet::Draw1(Int_t mode)
{
  const Int_t ntit=4;
  Char_t *tit[ntit]={"BG2000","BG4000","BG8000","HijingCentral"};

  fC1 = TAliasPAI::canvas("Mean BG ch"); 
  TString st("Mean BG characteristic : ");
  st += fMdst->GetTitle();
  fP1 = TAliasPAI::pad(st.Data(), fC1);
  Int_t ny=2, ip=1;
  if(mode>0) {
     ny=3;
     fP1->Divide(1,ny);
  } else {
     fP1->Divide(ny,1);
  }

  Double_t y=0., x=0.;
  if(mode>0) {
     fP1->cd(ip++);
     TAliasPAI::drawHist(fHNcellEMCAL, 2);
     TAliasPAI::drawHist(fHNcellCHTR,  3, 2, "same");
     TAliasPAI::drawHist(fHNcellGRID,  2, 3, "same");
     fHNcellEMCAL->SetMaximum(fHNcellCHTR->GetMaximum()*1.1);

     y = fHNcellEMCAL->GetMaximum();
     x = 4000.;
     if(fMdst->IsPythiaDst()) x=200;

     for(Int_t i=0; i<3; i++) {
        TAliasPAI::lat(text[i], x, y*(0.8-0.2*i), 12, 0.06, i+1);
     }
  }

  fP1->cd(ip++);
  TAliasPAI::drawHist(fHEtCellEMCAL,2);
  TAliasPAI::drawHist(fHEtCellCHTR,3,2,"same");
  TAliasPAI::drawHist(fHEtCellGRID,2,3,"same");
  gPad->SetLogy(1);
  if(mode==0) {
     y = fHEtCellEMCAL->GetMaximum();
     if(fMdst->IsPythiaDst()) y = fHEtCellCHTR->GetMaximum();
     if(gPad->GetLogy()) y = TMath::Log10(y);
     x = 10.;
     TAliasPAI::dispLine(x, y*0.8, 1,"EMCAL", 200., 1, 2, 400., 0.0, 0.06);
     TAliasPAI::dispLine(x, y/0.6, 2,"Ch.P", 200., 2, 2, 400., 0.0, 0.06);
     TAliasPAI::dispLine(x, y/0.4, 2,"GRID", 200., 3, 2, 400., 0.0, 0.06);
  }

  fP1->cd(ip++);
  TAliasPAI::drawHist(fHSumEMCAL, 2); 
  TAliasPAI::drawHist(fHSumCHTR, 3, 2, "same");
  TAliasPAI::drawHist(fHSumGRID, 2, 3, "same");
  fC1->Update();
}

void AliBgInfoForJet::Draw2(Int_t mode)
{
  fC2 = TAliasPAI::canvas("BGvsEtCut"); 
  fC2->SetTitle("Canvas2");
  TString st("BG characteristic vs Et cut : ");
  st += text[0]; st += " : ";
  st += fMdst->GetTitle();
  fP2 = TAliasPAI::pad(st.Data(), fC2);
  fP2->Divide(2,2);

  Int_t ip=1;
  fP2->cd(ip++);
  gStyle->SetOptStat(0);
  TAliasPAI::drawHist(fHOccupancyEMCAL, 2);
  TAliasPAI::titles(fHOccupancyEMCAL, " E_{T}(cell) "," # cell above th. ");
  gPad->SetLogy(1);

  fP2->cd(ip++);
  gStyle->SetOptStat(0);
  TAliasPAI::drawHist(fHncellmore0EMCAL, 2);
  TAliasPAI::titles(fHncellmore0EMCAL, " E_{T}(cell) "," % of events ");

  fP2->cd(ip++);
  gStyle->SetOptStat(0);
  TAliasPAI::drawHist(fHEtMeanEMCAL, 2);
  TAliasPAI::titles(fHEtMeanEMCAL, " E_{T}(cell) "," mean E_{T}(EMCAL)");
  gPad->SetLogy(1);

  fP2->cd(ip++);
  gStyle->SetOptStat(0);
  TAliasPAI::drawHist(fHEtRMSEMCAL, 2);
  TAliasPAI::titles(fHEtRMSEMCAL, " E_{T}(cell) "," RMS E_{T}(EMCAL)");
  gPad->SetLogy(1);
  // fluctuation in jet cone
  char ctmp[200];
  Double_t rjet=0., c=0, ym=0.;
  TH1F *htmp=0;
  ym = fHEtRMSEMCAL->GetMaximum();
  printf(" ym %f \n", ym);
  for(Int_t i=0; i<3; i++){
    rjet = 0.5 - 0.1*i;
 // square of jet cone to square of EMCAL : (pi*r**2)/(2.*pi/3.*1.4)
    c = (rjet*rjet)*(2./3.*1.4);
    htmp = (TH1F*)fHEtRMSEMCAL->Clone();
    sprintf(ctmp,"hEtRmsR%3.1f",rjet);
    htmp->SetName(ctmp);
    sprintf(ctmp,"Et RMS for jet with R %3.1f",rjet);
    htmp->SetTitle(ctmp);
    htmp->Scale(TMath::Sqrt(c)); // 26-mar-2003 Dy = C*Dx -> rms(y) = sqrt(c)*rms(x)
    TAliasPAI::drawHist(htmp, 2, i+2, "same");
    sprintf(ctmp,"R %3.1f, c %5.2f",rjet, c);
    TAliasPAI::lat(ctmp, 4.0, TMath::Log10(ym)*(1.-0.2*i), 12, 0.05, i+2);
  }

  fC2->Update();
}

void AliBgInfoForJet::Draw3(Int_t iv, Int_t mode)
{ 
  printf("<I> AliBgInfoForJet::Draw3 : iv %i mode %i\n", iv, mode);
  if(iv<0 || iv>2) {
     iv = 0;
     printf("Wrong iv : iv = %i now : %s \n", iv, text[iv]);
  }
  Double_t xb=0.0, yb=0.0, mean=0.0, rms=0.0, pRes=0.0; // relative resolution - 7-apr-2003
  TString st("BGvsEtCut"), titx(" E_{T}(cell) ");
  if(iv==1) titx = "P_{t}(track)";
  st += text[iv];
  fC3 = TAliasPAI::canvas(st.Data());
  st  = "Canvas3"; st += iv;
  fC3->SetTitle(st.Data());

  st  = "BG characteristic vs Et cut : ";
  st += text[iv]; st += " : ";
  st += fMdst->GetTitle();
  fP3 = TAliasPAI::pad(st.Data(), fC3);
  fP3->Divide(2,2);

  Int_t ip=1;
  fP3->cd(ip++);
  gStyle->SetOptStat(0);
  fhOccu[iv]->SetMinimum(0.01);
  TAliasPAI::drawHist(fhOccu[iv], 2);
  TAliasPAI::titles(fhOccu[iv], titx.Data()," # cell above th. ");
  gPad->SetLogy(1);

  fP3->cd(ip++);
  gStyle->SetOptStat(0);
  TAliasPAI::drawHist(fhNcellmore0[iv], 2);
  TAliasPAI::titles(fhNcellmore0[iv], titx.Data()," % of events ");

  TString opt;
  fP3->cd(ip++);
  gStyle->SetOptStat(0);
  TH1F *hidNoErr = 0;
  TString nam=fhEtMean[iv]->GetName(); nam += "c";
  hidNoErr = TAliasPAI::copy(fhEtMean[iv], nam.Data(), "", 0);
  if(fMdst->IsPythiaDst()) opt="c";
  TAliasPAI::drawHist(hidNoErr, 4, 1, opt.Data());
  TAliasPAI::titles(hidNoErr, titx.Data()," mean E_{T}(EMCAL)");
  // In cone 
  TH1D *hid=0, *hidc=0, *hidRMS[3];
  if(fMdst->IsPythiaDst()) {
     TH2F *hid2 = fHSumJet[iv];
     Char_t nam[100];
     for(Int_t ix=1; ix<=3; ix++){
        sprintf(nam,"hMeanEtCone0%i%s\n", 2+ix, text[iv-1]); 
        hid = hid2->ProjectionX(nam,ix,ix,"e"); // "e" - is very important
        sprintf(nam,"hcMeanEtCone0%i%s\n", 2+ix, text[iv-1]); 
        hidc = TAliasPAI::copy(hid, nam);
        TAliasPAI::drawHist(hidc, 2, 5-ix, "SAME C");

        sprintf(nam,"hRmsEtCone0%i%s\n", 2+ix, text[iv-1]);
        hidRMS[ix-1] = TAliasPAI::create(hid, nam, "rms");
        yb = Ra[ix-1] + 0.01;
        printf(" Projection %i : yb %5.2f \n", ix, yb); 
        for(Int_t i=1; i<=hid->GetNbinsX(); i++){
	  xb   = hid->GetBinCenter(i);
          mean = hid->GetBinContent(i);  
          rms  = hid->GetBinError(i);  
	  hidRMS[ix-1]->Fill(xb, rms);
          if(mean>0.0) {
	     pRes = 100.*rms/mean; 
	     fHRelResJet[iv]->Fill(xb, yb, pRes);
             printf(" xb %f mean %f rms %f pRes %5.1f%%\n", xb, mean, rms, pRes); 
          }
        }
     }
  } else {
     gPad->SetLogy(1);
  }

  fP3->cd(ip++);
  gStyle->SetOptStat(0);
  fhEtRMS[iv]->SetMinimum(0.1);
  TAliasPAI::drawHist(fhEtRMS[iv], 4, 1, opt.Data());
  TAliasPAI::titles(fhEtRMS[iv], titx.Data()," RMS E_{T}(EMCAL)");
  // fluctuation in jet cone
  char ctmp[200];
  Double_t rjet=0., c=0, ym=0.;
  TH1F *htmp=0;
  ym = fhEtRMS[iv]->GetMaximum();
  printf(" (%s) ym %f \n", text[iv], ym);

  if(fMdst->IsPythiaDst()) {
     for(Int_t i=0; i<3; i++){
        TAliasPAI::drawHist(hidRMS[i], 4-i, i+2, "SAME C");
        rjet = 0.5 - 0.1*i;        
        sprintf(ctmp,"R %3.1f ",rjet);
        TAliasPAI::lat(ctmp, 6.0, ym * (0.8-0.1*i), 12, 0.05, i+2);
     }
  } else {
     gPad->SetLogy(1);
     for(Int_t i=0; i<3; i++){
        rjet = 0.5 - 0.1*i;        
    // square of jet cone to square of EMCAL : (pi*r**2)/(2.*pi/3.*1.4)
        c = (rjet*rjet)*(2./3.*1.4);
        htmp = (TH1F*)fhEtRMS[iv]->Clone();
        sprintf(ctmp,"hEtRmsR0%1i",Int_t(10.*rjet));
        htmp->SetName(ctmp);
        sprintf(ctmp,"Et RMS for jet with R %f3.1",rjet);
        htmp->SetTitle(ctmp);
        htmp->Scale(TMath::Sqrt(c)); // 26-mar-2003 Dy = C*Dx -> rms(y) = sqrt(c)*rms(x)
        TAliasPAI::drawHist(htmp, 2, i+2, "same");
        sprintf(ctmp,"R %3.1f, #sqrt{c} %5.2f",rjet, TMath::Sqrt(c));
        TAliasPAI::lat(ctmp, 6.0, TMath::Log10(ym)*(1.-0.2*i), 12, 0.05, i+2);
        printf(" %s -> \n", ctmp, htmp->GetMaximum());
     }
  }

  fC3->Update();
}

void AliBgInfoForJet::Draw4()
{
  if(!fMdst->IsPythiaDst()) return;
 
  fC4 = TAliasPAI::canvas("RelRes"); 
  fC4->SetTitle("Canvas4");
  TString st("Relative resolution : "), opt;
  //  st += text[0]; st += " : ";
  st += fMdst->GetTitle();
  fP4 = TAliasPAI::pad(st.Data(), fC4);
  fP4->cd();
  fP4->Divide(2,2);
  TH1D *hid = 0;
  Char_t nam[100];
  Int_t col[3]={4,3,2};
  for(Int_t i=0; i<3; i++){
     TH2F *hid2 = fHRelResJet[i];
     fP4->cd(i+1);
     opt = "";
     for(Int_t ix=1; ix<=3; ix++){
        sprintf(nam, "hRelRes%sJR%i", text[i], Int_t(10.*Ra[ix]));
        hid = hid2->ProjectionX(nam,ix,ix);
        hid->SetTitle(text[i]);
        TAliasPAI::drawHist(hid, 5-ix, col[ix-1], opt.Data(), 1);
        opt = "same";
     }
  }

  fP4->cd(4);
  for(Int_t ix=1; ix<=3; ix++){
     sprintf(nam, "Jet Radius %3.1f", Ra[ix]);       
     TAliasPAI::lat(nam, 0.1, 1.-0.2*ix, 12, 0.08, col[ix-1]);
  }  

  fC4->Update();
}

void AliBgInfoForJet::Print1(Int_t mode)
{
}

void AliBgInfoForJet::Browse(TBrowser* b)
{
   if(fListHist1)  b->Add((TObject*)fListHist1);
   if(fListHist2)  b->Add((TObject*)fListHist2);
}

Bool_t AliBgInfoForJet::IsFolder() const  
{
  if(fListHist1 || fListHist2) return kTRUE;
  else          return kFALSE;
}