//
//   Clustering; Feb 16, 2006; use rHits.C staff too
//               Jun 01, 2006 - tune for WSU Linux cluster
//               Nov 01, 2006 - new loop of pi0/gamma discriminations
//
#include "rHits.h"
#include "clust.h"
#include "/home/pavlinov/macros/macroIncludePai.h"
#include "/home/pavlinov/macros/TAliasPAI.h"
#include "/home/pavlinov/macros/macroIncludeAlice.h"
#include "AliReconstruction.h"

extern TString file;
extern AliRunLoader *RL;
extern AliEMCALLoader *eg;
extern  AliEMCALGeometry *emcalg;
// hists 
extern TList *lh, *lJetHists, *lGH, *lHitsHists, *lKineHists;

class TGeoNode; TGeoNode *XEN1 = 0;

class AliEMCALClusterizerv1; AliEMCALClusterizerv1 *cl=0;

TList *lhRP=0, *lhPC=0;
Double_t MOMENTUM=1.;
TString DivOpt; // 1X1, 2X2 or 3X3
char *optList[]={"1X1","2X2","3X3"};
int   nOptList = sizeof(optList)/sizeof(char*);

void clust(int nmax, int var,const char* opt)
{ 
  printf("\n <<<<<<<<<<<<<<< Cluster finder started  (clust(int nmax, int var,const char* opt)) <<<<<<<<<<<<<<< \n");
  gROOT->ProcessLine(".!rm EMCAL.RecPoints*.root"); // clean up old staff  
  if(nmax);
  defineInputDirectory(var,opt);
  if(!initAlice(file.Data())) return;
  double p = pai::getMomentumFromDirName(), eThSeed=0.1, minEofCell=0.016;
  if(p<=1.) eThSeed = 0.05;

  cl = new AliEMCALClusterizerv1(file.Data(),"Folders");
  cl->SetECAClusteringThreshold(eThSeed);
  cl->SetMinECut(minEofCell);
  if(nmax>0) cl->SetEventRange(0, nmax);

  gROOT->GetListOfBrowsables()->Add(cl);

  //  cl->ExecuteTask("tim all"); // "deb all" will fill histogramms
  cl->ExecuteTask("tim all"); // "deb all" will fill histogramms
  //cl->ExecuteTask("tim deb all"); // "deb all" will fill histogramms
  //cl->ExecuteTask(""); // "deb all" will fill histogramms
  printf("\n >>>>>>>>>>>>>>> Cluster finder ended  >>>>>>>>>>>>>> \n");

}

void digiTest(int nmax, int var,const char* opt)
{
  defineInputDirectory(var,opt);

  if(!initAlice(file.Data())) return;

  int maxEvent = RL->GetNumberOfEvents(), pri=0, fr=100;
  if(nmax>maxEvent || nmax<=0) nmax = maxEvent;
  printf("events in file %i -> %i will be read : fr %i \n", maxEvent, nmax, fr);
  TClonesArray  *digits=0;

  for (int nev=0; nev < nmax; nev++) {
    RL->GetEvent(nev);
    RL->LoadKinematics() ;
    RL->LoadDigits() ;

    pri = 0;
    if(nev==nmax-1 || nev%fr==0) pri = 1;

    digits = eg->Digits();

    for(int idigi=0; idigi<digits->GetEntries(); idigi++){
      cl->ExecuteTask();
    }

    if(pri) printf(" ******   Read event \t%i : N \t%i ********** \n", nev, digits->GetEntries());

  }
}

void recPoint(int nmax, int var,const char* opt)
{ // Corr. Nov 02,2006
  printf("\n <<<<<<<<<<<<<<< Rec.points analysis is started  (recPoint(int nmax, int var,const char* opt)) <<<<<<<<<<<<<<< \n");

  TString OPT(opt);
  defineInputDirectory(var,opt);
  if(!initAlice(file.Data())) return;
  double momSig = 0.1*TMath::Sqrt(MOMENTUM), eRP=0.0, eSumRP=0., eSumPC=0.;
  if(MOMENTUM >= 20.) eRP = MOMENTUM - 0.5; // nonlinearity 

  int maxEvent = RL->GetNumberOfEvents(), pri=0, fr=100;
  if(nmax>maxEvent || nmax<=0) nmax = maxEvent;
  printf("events in file %i -> %i will be read : fr %i \n", maxEvent, nmax, fr);

  lhRP = bookHistsForRecPoints("RP");
  gROOT->GetListOfBrowsables()->Add(lhRP);
  lhPC = bookHistsForRecPoints("PC");
  gROOT->GetListOfBrowsables()->Add(lhPC);
  lKineHists = bookKineHists();
  gROOT->GetListOfBrowsables()->Add(lKineHists);
  lHitsHists = bookHitsHists(emcalg);
  lHitsHists->SetName("indexesForPC");
  gROOT->GetListOfBrowsables()->Add(lHitsHists);

  TObjArray *recPoints=0;
  const AliEMCALRecPoint *rp=0, *rp0=0, *rp1=0;
  double erp=0.0, effMass=0.0; // have to calculate
  float *eCells=0, lambda[2];
  int absIdCellMax=0, absId=0;
  TVector3 lpos, gpos, vg1, vg2;
  TParticle *p=0, *pg1=0, *pg2=0;
  double dphi=0., sphericity = 0;
  TLorentzVector pi0;

  int keyPi0 = 0, nRP=0, nPC=0; // 
  for (int nev=0; nev < nmax; nev++) {
    // As in AliEMCALReconstructor::FillESD - Feb 27,2007
    // I had problem with first event before !
    RL->LoadRecPoints();
    RL->LoadKinematics();
    RL->LoadDigits();
    RL->LoadHits();
    RL->GetEvent(nev);
 
    pri = 0;
    if(nev==nmax-1 || nev%fr==0) pri = 1;

    recPoints = eg->RecPoints();
    int nrp = recPoints->GetEntries();
    erp = 0.0;
    // Kinematics is o'k
    p = RL->Stack()->Particle(0);
    if(p->GetPdgCode() == 111) { // pi0; Nov 3,2006
      keyPi0 = 1; 
      if(p->GetNDaughters() == 2) { // gamma, gamma case
        pg1 = RL->Stack()->Particle(1);
        pg2 = RL->Stack()->Particle(2);
        vg1.SetXYZ(pg1->Px(), pg1->Py(), pg1->Pz());
        vg2.SetXYZ(pg2->Px(), pg2->Py(), pg2->Pz());
        pai::fillH1F(lhRP, 25, vg1.Angle(vg2)*TMath::RadToDeg());
      }
    }
    fillKineHists();

    if(pri) printf(" ******   Read event %i : %i rec.point  ********** \n", nev, nrp);
    nRP=nPC=0;
    TList *lh = 0;
    eSumRP = eSumPC = 0.;
    for(int irp=0; irp<nrp; irp++) {
      rp    = eg->RecPoint(irp);
      if(rp == 0) continue;
      if(rp->GetClusterType() == AliESDCaloCluster::kPseudoCluster) {
        lh = lhPC; nPC++;
        eSumPC += rp->GetPointEnergy();
      } else {
        if     (nRP==0) rp0 = rp;
        else if(nRP==1) rp1 = rp;
        lh = lhRP; nRP++;
        eSumRP += rp->GetPointEnergy();
      }
      pai::fillH1F(lh, 33,  double(rp->GetClusterType()));

      pai::fillH1F(lh, 1,  double(rp->GetMultiplicity()) );
      pai::fillH1F(lh, 2,  double(rp->GetCoreEnergy()) ); // have to be calculated
      pai::fillH1F(lh, 3,  double(rp->GetMaximalEnergy()) );

      eCells = rp->GetEnergiesList();
      erp = 0.0;
      int nmult = rp->GetMultiplicity(); 
      Float_t *ce = rp->GetEnergiesList();
      Float_t *timeFloat = rp->GetTimeList();
      Int_t   *digitInts = rp->GetAbsId();
      for(int ic=0; ic<rp->GetMultiplicity(); ic++) {
        erp  += double(eCells[ic]);
        absId = rp->GetAbsId(ic);
        if(rp->GetClusterType() == AliESDCaloCluster::kPseudoCluster)
        fillIndexHists(absId);       
        pai::fillH1F(lh, 34,  double(ce[ic]));
        pai::fillH1F(lh, 35,  double(timeFloat[ic]*1.e+11)); // 10 ps units
        pai::fillH1F(lh, 36,  double(digitInts[ic])); 
      }
      pai::fillH1F(lh, 4, erp);
      absIdCellMax =  rp->GetAbsId(0);
      if(absIdCellMax>=0) pai::fillH1F(lh, 12, double(absIdCellMax));
      
      // Local coordinates 
      rp->GetLocalPosition(lpos);
      for(int i=0; i<3; i++) pai::fillH1F(lh, 5+i, lpos[i]); 
      pai::fillH1F(lh, 8, rp->GetDispersion());

      rp->GetGlobalPosition(gpos);
      // Global coordinates
      for(int i=0; i<3; i++) pai::fillH1F(lh, 13+i, gpos[i]); 
      pai::fillH2F(lh, 16, gpos[0], gpos[1]);
      pai::fillH2F(lh, 17,  TVector2::Phi_0_2pi(gpos.Phi())*TMath::RadToDeg(), gpos[2]);
      // Compare reco vs initial for RP
      if(lh == lhPC) continue;

      dphi =  TVector2::Phi_0_2pi(p->Phi()) - TVector2::Phi_0_2pi(gpos.Phi());   
      pai::fillH1F(lh, 9,  p->Energy() - erp);
      pai::fillH1F(lh, 10, p->Theta() - gpos.Theta());            
      pai::fillH1F(lh, 11, dphi);
      if(dphi>0.1){
	//        printf(" p->Phi() %f : gpos.Phi() %f : dphi %f \n", p->Phi(), gpos.Phi(), dphi);
      }
      rp->GetElipsAxis(lambda);
      if(lambda[0]+lambda[1]) {
        sphericity = (lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
      }
      if(nRP == 1) {
        pai::fillH1F(lh, 24, rp->GetPointEnergy());
          for(int i=0; i<2; i++) pai::fillH1F(lh, 22+i, double(lambda[i]));
	  if(rp->GetPointEnergy() > (eRP-3.*momSig)){
            for(int i=0; i<2; i++) pai::fillH1F(lh, 26+i, double(lambda[i]));
            pai::fillH1F(lh, 28, sphericity);
	  }
	  // Mip business
          pai::fillH1F(lh, 29, erp);
          if     (nmult==1) pai::fillH1F(lh, 30, erp);
          else if(nmult==2) pai::fillH1F(lh, 31, erp);
          else if(nmult >2) pai::fillH1F(lh, 32, erp);
      }
    }
    pai::fillH1F(lhPC, 0, double(nPC));
    pai::fillH1F(lhPC, 37, eSumPC);
    pai::fillH1F(lhRP, 0, double(nRP));
    pai::fillH1F(lhRP, 37, eSumRP);

    if(nRP==2) {
      effMass = effectiveMass(rp0, rp1, pi0);
      pai::fillH1F(lhRP, 18, effMass);
      pai::fillH1F(lhRP, 20, pi0.Vect().Mag());
    }

  }

  TH1F *heffMass = (TH1F*)lhRP->At(18);
  if(heffMass->GetEntries() > 100.) drawPi0(lhRP);

  //  drawLocalCoordHists(lhRP);
  // drawInitialKine(lKineHists);
  // if(lHitsHists) drawIndexes(lHitsHists);

  pai::saveListOfHists(lhRP, "rpHists.root", kTRUE);
  pai::saveListOfHists(lhPC, "rpHists.root", kTRUE,"UPDATE");
  pai::saveListOfHists(lKineHists, "rpHists.root", kTRUE,"UPDATE");
  pai::saveListOfHists(lHitsHists, "rpHists.root", kTRUE,"UPDATE");
  
  TH1F *h = (TH1F*)lhRP->At(33);
  printf(" #number of pseudo cluster %i \n", int(h->GetBinContent(h->GetXaxis()->FindBin(0.))));
  printf(" #number of rec.points     %i \n", int(h->GetBinContent(h->GetXaxis()->FindBin(1.))));

  printf("\n >>>>>>>>>>>>>>> Rec.points analysis is ended  (recPoint(int nmax, int var,const char* opt)) >>>>>>>>>>>>>>> \n");
}

void rpWithEsd(int n, int deb)
{ // Mar 6, 2007
  /*
.L /home/pavlinov/ALICE/SHISHKEBAB/rHits.C+
.L /home/pavlinov/ALICE/SHISHKEBAB/clust.C+
rpWithEsd();
   */
  if(deb>0) AliLog::SetGlobalDebugLevel(deb);
  AliReconstruction *rec = new AliReconstruction;
  gROOT->GetListOfBrowsables()->Add(rec);

  rec->SetRunTracking("");
  rec->SetRunVertexFinder(kFALSE);

  rec->SetRunLocalReconstruction("EMCAL");  //only do emcal
  rec->SetFillESD("EMCAL");

  //  if(n>0) rec->SetEventRange(0, n); // Mar 12, 2007 - for testing
  rec->SetEventRange(3510, 10); // Mar 12, 2007

  rec->Run("");

  printf(" # events %i : deb level %i \n", n, deb);
}

double effectiveMass(const AliEMCALRecPoint *rp1, const AliEMCALRecPoint *rp2, TLorentzVector &pi0)
{ // Nov 2, 2006 - just for 2 RP
  static TLorentzVector p1, p2; 
  static double em=0.0, angle=0.0;
  static TVector3 diff, v2;

  getLorentzVector(rp1, p1); 
  getLorentzVector(rp2, p2);

  pi0 = p1 + p2;
  angle = p1.Angle(p2.Vect());
  pai::fillH1F(lhRP, 19, angle*TMath::RadToDeg());

  rp1->GetGlobalPosition(diff);
  rp2->GetGlobalPosition(v2);
  diff -= v2;
  pai::fillH1F(lhRP, 21, diff.Mag());

  //em = TMath::Sqrt(2.*p1.E()*p2.E()*(1. - TMath::Cos(angle)));
  em = pi0.M();
  return em;
}

void getLorentzVector(const AliEMCALRecPoint *rp, TLorentzVector &p)
{
  static double e;
  static TVector3 gpos;

  e = rp->GetPointEnergy();
  rp->GetGlobalPosition(gpos);
  gpos.SetMag(e);
  p.SetVectM(gpos, 0.0);
}

TList *bookHistsForRecPoints(const char *opt)
{
  double ADCchannelEC = 0.00305; // ~3mev per adc count  
  if(opt);
  gROOT->cd();

  int nmax=100;
  double p = pai::getMomentumFromDirName(), emax=0.0, maxAngle=0.0;
  if(p<0) emax = 10.;
  else    emax = 2.*p;
  MOMENTUM = p;
  TString tv(Form("%iGEV", int(p)));
  maxAngle = 0.135/p*10.;
  if(p>10) maxAngle *= 1.6; 
  printf(" bookHistsForRecPoints : %s : maxAngle = %7.3f\n", tv.Data(), maxAngle*TMath::RadToDeg());

  new TH1F(Form("00_hNum%s%s",opt,tv.Data()), Form("#rec.points : %s", tv.Data()), 101, -0.5,100.5);
  new TH1F(Form("01_hNumCellsIn%s%s",opt,tv.Data()), Form("#cells in rec.points : %s", tv.Data()), 41, -0.5,40.5);  
  new TH1F(Form("02_hECore%s%s",opt,tv.Data()), Form("core energy of rec.points : %s", tv.Data()), nmax, 0.0,emax);
  new TH1F(Form("03_hEMaximalEenergy%s%s",opt,tv.Data()), Form("highest energy in the rec.points : %s", 
  tv.Data()), nmax, 0.0,emax);
  new TH1F(Form("04_hEof%s%s",opt,tv.Data()), Form("energy of rec.points : %s", tv.Data()), nmax, 0.0,emax); //4 

  new TH1F(Form("05_hXL%s%s",opt,tv.Data()), Form("local X : %s", tv.Data()), 120, -3.0,3.0);    // 5
  new TH1F(Form("06_hYL%s%s",opt,tv.Data()), Form("local Y : %s", tv.Data()), 140,  -70.0,70.0); // 6
  new TH1F(Form("07_hZL%s%s",opt,tv.Data()), Form("local Z : %s", tv.Data()), 180, -180.,180.0); // 7
  new TH1F(Form("08_hDisp%s%s",opt,tv.Data()), Form("shower dispersion: %s", tv.Data()), 100,0.0,10.0); // 8
  // Kinematics vs clusters  
  double sig = 0.1*sqrt(p);
  new TH1F(Form("09_hEkine-Eclust%s%s",opt,tv.Data()),Form("E(kine) - E(RP): %s", tv.Data()),100,-3*sig, 3*sig); // 9
  new TH1F(Form("10_hThetakine-Thetaclust%s%s",opt,tv.Data()), Form("#theta(kine)-#theta(RP): %s", tv.Data()), 100,-0.01,+0.01); // 10
  new TH1F(Form("11_hPhikine-Phiclust%s%s",opt,tv.Data()), Form("#phi(kine)-#phi(RP): %s", tv.Data()),100, -0.01,+0.01); // 11
  // Id with max energy
  int ncells = emcalg->GetNCells();
  new TH1F("12_hidOfMaxCellOf%s%s"," id of max cell from RP ", ncells, 0.5, double(ncells)-0.5); // 12
  // Global coordinates
  new TH1F(Form("13_hXGlob%s%s",opt,tv.Data()), Form("global X : %s", tv.Data()), 1000, -500., 500.);  // 13
  new TH1F(Form("14_hYGlob%s%s",opt,tv.Data()), Form("global Y : %s", tv.Data()), 1000, -500., 500.);  // 14
  new TH1F(Form("15_hZGlob%s%s",opt,tv.Data()), Form("global Z : %s", tv.Data()), 700,  -350., 350.);  // 15
  new TH2F(Form("16_hXYGlob%s%s",opt,tv.Data()), Form("global X:Y : %s", tv.Data()), 200,-500.,500., 200,-500.,500.);  // 16
  new TH2F(Form("17_hPhiZGlob%s%s",opt,tv.Data()), Form("global #phi:z : %s", tv.Data()), 110,80.,190., 60,-300.,+300.);  // 17
  // Nov 01, 2006
  new TH1F(Form("18_hGGEfMass%s%s",opt,tv.Data()), Form("eff.mass of RP,RP : %s ", tv.Data()), 
  100,  0.0, 0.5);  // 18
  new TH1F(Form("19_hAlphaOfRpRp%s",tv.Data()), Form("#alpha(degree) of RP,RP : %s ", tv.Data()), 
	   60,  0.0, maxAngle*TMath::RadToDeg());  // 19
  new TH1F(Form("20_hMomOfRpRp%s",tv.Data()), Form("#momentum of RP,RP : %s ", tv.Data()), 
  200,  0.0, emax);  // 20
  new TH1F(Form("21_hDistOfRpRp%s",tv.Data()), Form("distance between of RP,RP : %s ", tv.Data()), 
	   100,  0.0, maxAngle*500.);  // 21
  // One cluster 
  new TH1F(Form("22_hMaxL1_1cl%s",tv.Data()), Form("max #lambda_{1} - 1 cluster", tv.Data()), 300,  0.0, 3.);  // 22
  new TH1F(Form("23_hMinL2_1cl%s",tv.Data()), Form("min #lambda_{2} - 1 cluster", tv.Data()), 300,  0.0, 3.);  // 23
  new TH1F(Form("24_hMom_1cl%s",tv.Data()), Form("Momentum with 1 cluster", tv.Data()), 200,  0.0, emax);  // 24
  new TH1F(Form("25_hAlphaOfGG%s",tv.Data()), Form("#alpha(degree) of #gamma,#gamma(initial kine) : %s ", tv.Data()), 
	   60,  0.0, maxAngle*TMath::RadToDeg());  // 25
  // With some cuts 
  new TH1F(Form("26_hMaxL1_1cl_cuts%s",tv.Data()), Form("max #lambda_{1} with energy cut", tv.Data()), 300,  0.0, 3.);  // 26
  new TH1F(Form("27_hMinL2_1cl_cuts%s",tv.Data()), Form("min #lambda_{2} with energy cut", tv.Data()), 300,  0.0, 3.);  // 27
  new TH1F(Form("28_hSphericity%s",tv.Data()), Form("Sphericity (#lambda_{1}-#lambda_{2})/(#lambda_{1}+#lambda_{2})", tv.Data()), 100,  0.0, 1.);  // 28
  /* Mip staf - see AliEMCALDigitizer::InitParameters()
  fADCchannelEC    = 0.00305; // 200./65536 - width of one ADC channel in GeV
  fADCpedestalEC   = 0.009 ;  // GeV
  */
  double fADCchannelEC = 0.00305;
  int nmip = 300;
  double xmi = fADCchannelEC/2., xma = xmi + fADCchannelEC*nmip;
  new TH1F(Form("29_hMip_%s",tv.Data()), Form("MIP(energy of rec.points) : %s", tv.Data()), nmip,xmi,xma); // 29
  new TH1F(Form("30_hMipNcell=1_%s",tv.Data()), Form("MIP(energy of rec.points) ncell=1 : %s", tv.Data()), nmip,xmi,xma); // 30
  new TH1F(Form("31_hMipNcell=2_%s",tv.Data()), Form("MIP(energy of rec.points) ncell=2 : %s", tv.Data()), nmip,xmi,xma); // 31
  new TH1F(Form("32_hMipNcell>2_%s",tv.Data()), Form("MIP(energy of rec.points) ncell>2 : %s", tv.Data()), nmip,xmi,xma); // 32
  // --
  new TH1F(Form("33_hClusteType%s",tv.Data()), Form("# number of clusters (0-PC, 1-RP) : %s", 
						    tv.Data()), 4,-1.5,2.5);
  int nmax1=10000;
  xmi = ADCchannelEC/2.;
  xma = xmi + ADCchannelEC*nmax1;
  new TH1F(Form("34_hCellEnergy%s%s",opt,tv.Data()), Form("cell(digit) energy in the %s: %s", 
							  opt, tv.Data()), nmax1,xmi,xma);
  new TH1F(Form("35_hCellTime%s%s",opt,tv.Data()), Form("cell(digit) time in the %s: %s", opt, tv.Data()),
							  1000, 0.0, 3.e+3); // ns/100 = 10 ps
  int nmaxCell = emcalg->GetNCells();
  new TH1F(Form("36_hCellAbsId%s%s",opt,tv.Data()), Form("cell(digit) Abs Id in the %s: %s", opt, tv.Data()),
  nmaxCell, -0.5, double(nmaxCell)-0.5);

  new TH1F(Form("37_hEofSumOf%s%s",opt,tv.Data()), Form("sum energy of %s  : %s", opt, tv.Data()), nmax, 0.0,emax); // Mar 5, 2007
  return pai::moveHistsToList(Form("lh%s",opt),kTRUE);
}

// Read all lists - Nov 10,2006
TList *readAllLists()
{ 
  // Corr Feb 08, 2007
  MOMENTUM = pai::getMomentumFromDirName();
  TString cwd(gSystem->DirName(""));
  TString stRemove("/data/r22b/ALICE/PROD/");
  if(cwd.Contains("/eliza5")) stRemove = "/eliza5/alice/pavlinov/PROD/"; // PDSF
  cwd.ReplaceAll(stRemove.Data(), "");
  cwd.ReplaceAll("/","_");
  gROOT->ProcessLine(".!ls hitsP*.root >  tmp.lst");
  gROOT->ProcessLine(".!ls sdigi*.root >> tmp.lst");
  gROOT->ProcessLine(".!ls digi*.root  >> tmp.lst");
  gROOT->ProcessLine(".!ls rp*.root    >> tmp.lst");

  ifstream fin;
  fin.open("tmp.lst");
  char comment[200];
  int  nf=0;

  TList *lall = new TList;
  lall->SetName("readAllLists");
  if(lall) gROOT->GetListOfBrowsables()->Add(lall);

  while(!fin.getline(comment,200).eof()) {
    TString nls(Form("%s_%i",cwd.Data(),nf++));
    TList *l = pai::readAll(comment, nls.Data());
    lall->Add(l);
    if(strstr(comment,"rpHists")) {
      lhRP = (TList*)l->At(0); 
    } else if(strstr(comment,"hits")) {
      lKineHists = (TList*)l->FindObject("kineHists"); 
    }
  }
  for(int i=0; i<nOptList; i++) if(cwd.Contains(optList[i])) {DivOpt = optList[i]; break;}
  printf(" readAllLists() : read %i lists : DivOpt -> %s\n", nf, DivOpt.Data());
  // Get Geometry from galice.root
  getGeometry("galice.root"); // see rHits.C
  gROOT->GetListOfBrowsables()->Add(emcalg);
  lall->Add(emcalg);
  return lall;
}

// Drawing - Jun 08, 2006
void drawLocalCoordHists(TList *l)
{
  if(l==0) l = lhRP;
  if(l==0) return;

  TCanvas *c = pai::canvas("RPLocal","RPL");  
  TPad    *pad = pai::pad("EMCAL local coordinates",c);
  pad->Divide(2,2);

  gStyle->SetOptStat(110001);
  pad->cd(1);
  TH1F* hXl = (TH1F*)l->At(5);
  pai::draw(hXl,2);

  pad->cd(2);
  TH1F* hYl = (TH1F*)l->At(6);
  pai::draw(hYl,2);

  pad->cd(3);
  TH1F* hZl = (TH1F*)l->At(7);
  pai::draw(hZl,2);

  c->Update();
}

void drawGlobalCoordHists(TList *l)
{
  if(l==0) l = lhRP;
  if(l==0) return;

  TCanvas *c = pai::canvas("RPglobal","RPL");  
  TPad    *pad = pai::pad("EMCAL global coordinates of RP",c);
  pad->Divide(2,2);

  gStyle->SetOptStat(110001);
  pad->cd(1);
  TH1F* h = (TH1F*)l->At(13);
  pai::draw(h,2);

  pad->cd(2);
  h = (TH1F*)l->At(14);
  pai::draw(h,2);

  pad->cd(3);
  h = (TH1F*)l->At(15);
  pai::draw(h,2);

  c->Update();
}

void drawRecoQA(TList *l)
{
  if(l==0) l = lhRP;
  if(l==0) return;

  TCanvas *c = pai::canvas("recoQA","recoQA");  
  TPad    *pad = pai::pad("QA for EMCAL reconstruction",c);
  pad->Divide(2,2);

  gStyle->SetOptStat(111111);

  pad->cd(1);
  TH1F *h = (TH1F*)l->At(9);
  pai::draw(h,2);

  pad->cd(2);
  h = (TH1F*)l->At(10);
  pai::draw(h,2);

  pad->cd(3);
  h = (TH1F*)l->At(11);
  pai::draw(h,3);

  c->Update();
}

void drawInitialKine(TList *l)
{
  if(l==0) {
    readAllLists();
    l = lKineHists;
  }
  if(l==0) return;

  TCanvas *c = pai::canvas("kine","kine");  
  TPad    *pad = pai::pad("Initial kinematics",c);
  pad->Divide(2,2);

  gStyle->SetOptStat(111011);
  pad->cd(1);
  TH1F* hPhi = (TH1F*)l->At(11);
  pai::draw(hPhi,2);

  pad->cd(2);
  TH1F* hTheta = (TH1F*)l->At(12);
  pai::draw(hTheta,2);

  pad->cd(3);
  TH1F* hVz = (TH1F*)l->At(15);
  pai::draw(hVz,2);

  if(0) {
    TH1F* hVx = (TH1F*)l->At(13);
    pai::draw(hVx,2,2,"same");

    TH1F* hVy = (TH1F*)l->At(14);
    pai::draw(hVy,2,2,"same");
  }
  c->Update();
}

void drawPi0(const char *opt)
{ 
  // Dec 12,2006
  TList *l[2];
  getRpLists(opt, l);
  TList *lpi0 = l[1];
  drawPi0(lpi0);
}

void drawPi0(TList *l)
{ // Nov 2,2006
  // Should stay at directory with output file
  if(l==0) {
    readAllLists();
    l = lhRP;
  }
  if(l==0) return;

  TH1F* h0 = (TH1F*)l->At(0);
  TCanvas *c = pai::canvas(Form("%s,pi0",DivOpt.Data()),"pi0");  
  TPad    *pad = pai::pad(Form("%s, %iGeV, %i evs, #pi^{0}->#gamma, #gamma", 
			       DivOpt.Data(), int(MOMENTUM), int(h0->GetEntries())),c);
  pad->Divide(2,3);

  pad->cd(1);
  TH1F* h1 = (TH1F*)l->At(18);
  gStyle->SetOptStat(1000000);
  gStyle->SetOptFit(111);
  if(h1->GetEntries() <= 0.0) return;
  TF1 *gausi = pai::gausi("Pi0", 0.1, 0.2, h1);
  gausi->SetParName(0,"N_{#pi^{0}}");
  gausi->SetParName(1,"m_{#pi^{0}}");
  gausi->SetParName(2,"#sigma");
  h1->Fit(gausi,"IL","",0.05,0.2);
  h1->Fit(gausi,"EM","",0.05,0.2);
  pai::draw(h1,2);
  double pi0Mass    = gausi->GetParameter(1);
  double errPi0Mass = gausi->GetParError(1);

  pad->cd(2);
  TH1F* h2 = (TH1F*)l->At(19);
  TH1F* h2_1 = (TH1F*)l->At(25);
  gStyle->SetOptFit(0);
  pai::draw(h2,2);
  pai::draw(h2_1,2,2,"same");
  h2->SetMaximum(h2_1->GetMaximum()*1.1); // initial kinematics

  pad->cd(3);
  TH1F* h3 = (TH1F*)l->At(20);
  gStyle->SetOptStat(1000000);
  gStyle->SetOptFit(111);
  h3->Fit("gausn");
  pai::draw(h3,2);
  double pi0Mom    = h3->GetFunction("gausn")->GetParameter(1);
  double errPi0Mom = h3->GetFunction("gausn")->GetParError(1);

  pad->cd(4);
  TH1F* h4 = (TH1F*)l->At(21);
  gStyle->SetOptFit(0);
  pai::draw(h4,2);

  pad->cd(5);
  TH1F* h5 = (TH1F*)l->At(22);
  TH1F* h5_1 = (TH1F*)l->At(23);
  gStyle->SetOptStat(000000);
  pai::draw(h5,2);
  pai::draw(h5_1,2,2,"same");
  if(h5_1->GetMaximum() > h5->GetMaximum()) h5->SetMaximum(h5_1->GetMaximum()*1.1);
  TLegend *leg = new TLegend(0.5,0.5, 0.8,0.8);
  leg->AddEntry(h5, "#lambda_{1} (largest)", "L");
  leg->AddEntry(h5_1,"#lambda_{2} (smallest)", "L");
  leg->Draw();

  pad->cd(6);
  TH1F* h6 = (TH1F*)l->At(24);
  gStyle->SetOptStat(1110);
  pai::draw(h6,2);

  c->Update();

  printf(" di gamma mass %5.4f +/- %5.4f (0.135 GeV) <- PDG value \n", pi0Mass ,errPi0Mass);
  printf("          mom  %5.2f +/- %5.2f (%f) <- initial value \n", pi0Mom ,errPi0Mom, MOMENTUM);
  // rp multiplicity
  TH1F* hrpMult = (TH1F*)l->At(0);
  printf(" ## Multiplicity ## \n");
  for(int i=1; i<=10; i++) {
    int nbin = int((*hrpMult)[i]);
    int mult = int(hrpMult->GetBinCenter(i));
    if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/hrpMult->GetEntries()); 
  }
}

void drawLambda(const char *opt)
{  // Nov 11, 2006
  /* const char *opt="2X2 50GEV two"
     drawLambda("2X2 50GEV two")
     drawLambda("3X3 1GEV")
     drawLambda("1X1 1GEV")
     drawLambda("2X2 5GEV")
   */ 
  TString strOpt(opt),prodDir("/data/r22b/ALICE/PROD/Oct2006/"); 
  TObjArray opts;
  pai::ParseString(strOpt, opts);
  printf(" drawLambda(\"%s\") : nopt = %i\n", opt, opts.GetEntries());
  TObjString *os = (TObjString*)opts.At(0);
  TString s0(os->GetString());
  os = (TObjString*)opts.At(1);
  TString s1(os->GetString()), s2;
  TList *lg = pai::readAll(Form("%s/%s/GAMMA/%s/rpHists.root", prodDir.Data(),s0.Data(),s1.Data()),
			   Form("lg_%s",s1.Data()));
  lg = (TList*)lg->At(0);
  TList *lpi0 = pai::readAll(Form("%s/%s/PI0/%s/rpHists.root", prodDir.Data(),s0.Data(),s1.Data()),
			   Form("lpi0_%s",s1.Data()));
  lpi0 = (TList*)lpi0->At(0);

  int ind0=25, ind = 1; // 25 - lambda with cut; 21 - no cut
  double xmax = 2.5;
  if(opts.GetEntries()>=3) {
    os = (TObjString*)opts.At(2);
    s2 = os->GetString();
    if(s2.Contains("two",TString::kIgnoreCase)) {
      ind++;
      xmax = 0.7;
    }
    printf(" %s : ind %i : xmax %5.2f\n", s2.Data(), ind, xmax);
  }

  TCanvas *c = pai::canvas("lambda","lambda", 10,20, 700, 500);  
  TPad    *pad = pai::pad(Form("%s, %s , #lambda_{%i} , #gamma vs #pi^{0}",s0.Data(), s1.Data(),ind),c);
  //  pad->Divide(2,3);
  pad->cd();

  TH1F* hgLam = (TH1F*)lg->At(ind0+ind);
  TH1F* hpi0Lam = (TH1F*)lpi0->At(ind0+ind);
  hpi0Lam->Scale(hgLam->Integral()/hpi0Lam->Integral());
  hgLam->SetAxisRange(0.0, xmax);

  gStyle->SetOptStat(0);

  pai::draw(hgLam,2);
  pai::draw(hpi0Lam,2,2,"same");
  if(hpi0Lam->GetMaximum() > hgLam->GetMaximum()) hgLam->SetMaximum(hpi0Lam->GetMaximum()*1.1);

  TLegend *leg = new TLegend(0.65,0.65, 0.85,0.85);
  leg->AddEntry(hgLam, Form("#gamma, %5.3f #pm %5.3f",hgLam->GetMean(),hgLam->GetRMS()), "L");
  leg->AddEntry(hpi0Lam, Form("#pi^{0}, %5.3f #pm %5.3f",hpi0Lam->GetMean(),hpi0Lam->GetRMS()), "L");
  leg->Draw();

  c->Update();
}

void sumLambda()
{
  // Nov 27, 2006
  // lambda = sqrt(variance); see - AliEMCALRecPoint::EvalElipsAxis
  // Dimension is length.

  // 1X1 
  double x1x1[]={1.,5.,10.,20.,30., 50.};
  double lg1x1[]={0.401,0.381,0.369,0.375,0.354,0.347};
  //double elg1x1[]={0.155,0.137,0.148,0.145,0.153, 0.150};
  double lpi1x1[]={0.387,1.035,0.675,0.494,0.439, 0.390};
  //double elpi1x1[]={0.172,0.302,0.161,0.118,0.137, 0.159};
  int    n1x1=sizeof(x1x1)/sizeof(double);
  // 2X2 
  double x2x2[]={1.,5., 10.,20.,30.,50.};
  double lg2x2[]={0.527,0.469, 0.460,0.457,0.457,0.448};
  double elg2x2[]={0.099,0.50, 0.063,0.090,0.095,0.056};
  double lpi2x2[]={0.518,0.649, 1.179,0.726,0.588,0.511};
  double elpi2x2[]={0.112,0.504,0.215,0.156,0.152,0.139};
  int    n2x2=sizeof(x2x2)/sizeof(double);
  double c2x2=1./2.;
  pai::multiply(n2x2, lg2x2,  c2x2);
  pai::multiply(n2x2,elg2x2,  c2x2);
  pai::multiply(n2x2, lpi2x2, c2x2);
  pai::multiply(n2x2,elpi2x2, c2x2);
  // 3X3 
  double x3x3[]={1.,5.,10.,20.,30.,50.};
  double lg3x3[]={0.603,0.519,0.501,0.490,0.490,0.487};
  double elg3x3[]={0.110,0.048,0.048,0.067,0.082,0.080};
  double lpi3x3[]={0.597,0.521,1.646,0.980,0.755,0.602};
  double elpi3x3[]={0.114,0.052,0.293,0.171,0.174,0.196};
  int    n3x3=sizeof(x3x3)/sizeof(double);
  double c3x3=1./3.;
  pai::multiply(n3x3, lg3x3,  c3x3);
  pai::multiply(n3x3,elg3x3,  c3x3);
  pai::multiply(n3x3, lpi3x3, c3x3);
  pai::multiply(n3x3,elpi3x3, c3x3);

  TCanvas *c = pai::canvas("sumLambda","sumLambda", 10,20, 700, 500);  
  TPad    *pad = pai::pad("Summary of #lambda_{1} , #gamma vs #pi^{0}, 1X1, 2X2, 3X3",c);
  pad->cd();

  TGraph *grPI3X3 = pai::drawGraph(n3x3, x3x3, lpi3x3, kBlue, 24, "AC*", 
		    " #lambda_{1} ", " momentum, GeV/C ", "#lambda_{1}", -1, "no");
  TH1* h = grPI3X3->GetHistogram();
  h->SetMinimum(0.0);
  h->SetMaximum(1.2);
  h->SetAxisRange(0.0, 60.);
  TGraph *grG3X3 = pai::drawGraph(n3x3, x3x3, lg3x3, kBlue, 24, "C*", 
			      "test", " momentum, GeV/C ", "#lambda_{1}", -1, "no");
  grG3X3->SetLineStyle(2);

  TGraph *grPI1X1 = pai::drawGraph(n1x1, x1x1, lpi1x1, kBlack, 24, "C*", 
			      "test", " momentum, GeV/C ", "#lambda_{1}", -1, "no");
  TGraph *grG1X1 = pai::drawGraph(n1x1, x1x1, lg1x1, kBlack, 24, "C*", 
			      "test", " momentum, GeV/C ", "#lambda_{1}", -1, "no");
  grG1X1->SetLineStyle(2);

  TGraph *grPI2X2 = pai::drawGraph(n2x2, x2x2, lpi2x2, kGreen, 24, "C*", 
			      "test", " momentum, GeV/C ", "#lambda_{1}", -1, "no");
  TGraph *grG2X2 = pai::drawGraph(n2x2, x2x2, lg2x2, kGreen, 24, "C*", 
			      "test", " momentum, GeV/C ", "#lambda_{1}", -1, "no");
  grG2X2->SetLineStyle(2);

  TLegend *leg = new TLegend(0.62,0.68, 0.95,0.98);
  leg->AddEntry(grPI1X1, "1X1", "lp");
  leg->AddEntry(grPI2X2, "2X2", "lp");
  leg->AddEntry(grPI3X3, "3X3", "lp");

  leg->Draw();

  c->Update();
}

TH1F* effGammaVsLam1(const char *opt)
{
  // Nov 29, 2006 - request of Tom; unfinished
  // Result is not so clear
  /*
    effGammaVsLam1("2X2 10GEV draw");
   */

  TList *l[2];
  getRpLists(opt, l);
  TList *lg = l[0];

  TH1F *h26 = (TH1F*)lg->At(26), *hOut=0;

  //  TString optEff(""); // if "r" 
  hOut = pai::efficiencyFactor(h26, "", 0);

  pai::draw(hOut,2);

  return hOut;
}

TGraphErrors* suppPi0VsEffGamma(const char *opt)
{
  // Dec 1, 2006; see drawHSvsEf(TH1F *h1,TH1F *h2, TH1F *hn1,TH1F *hn2) in rHits.C
  /*
    suppPi0VsEffGamma("2X2 1GEV");
    suppPi0VsEffGamma("2X2 10GEV");
   */
  gROOT->SetStyle("Plain");

  TList *l[2];
  getRpLists(opt, l);
  TList *lg = l[0], *lpi0 = l[1];

  if(lg==0 || lpi0==0) return 0;

  gROOT->cd();
  TH1::AddDirectory(kTRUE);

  TString optEff(""); //  ef for x<x0
  Int_t normGam = 0, normPi = 0; 
  TH1F *hg26 = (TH1F*)lg->At(26);
  //  TH1F *hgamNorm = (TH1F*)lg->At(0);
  TH1F *hgamNorm = hg26;
  if(hgamNorm) normGam = (Int_t)hgamNorm->Integral();
  TH1  *hef  = pai::efficiencyFactor(hg26,optEff, normGam);

  TH1F *hpi26 = (TH1F*)lpi0->At(26);
  // TH1F *hpiNorm = (TH1F*)lpi0->At(0);
  TH1F *hpiNorm = hpi26;
  if(hpiNorm) normPi = (Int_t)hpiNorm->Integral();
  TH1  *hsup  = pai::suppressionFactor(hpi26,optEff, normPi);

  Int_t n=hef->GetNbinsX();
  TArrayD *x  = new TArrayD(n);
  TArrayD *ex = new TArrayD(n);
  TArrayD *y  = new TArrayD(n);
  TArrayD *ey = new TArrayD(n);
 
  Int_t ind=0, markerColor=1;
  Double_t step=0.001, xminEff=0.1;
  for(Int_t i=n-1; i>=0; i--){      
     (*x)[ind] = hef->GetBinContent(i+1);
    (*ex)[ind] = hef->GetBinError(i+1);
     (*y)[ind] = hsup->GetBinContent(i+1);
    (*ey)[ind] = hsup->GetBinError(i+1);

    if( hef->GetBinContent(i+1) ||  hsup->GetBinContent(i+1))
      printf(" ef %f : sup %f : %f \n", hef->GetBinContent(i+1), hsup->GetBinContent(i+1),
      hef->GetBinContent(i+1) * hsup->GetBinContent(i+1));
    if((*x)[ind]  < xminEff) continue; 
    //    if(ind>0 && (*x)[ind] == 1.) {ind++;break;}
    if(ind>0 && TMath::Abs((*x)[ind]-(*x)[ind-1])<step) continue;
       //       || TMath::Abs((*y)[ind]-(*y)[ind-1])) continue;
    //    if((*y)[ind] < 1.0 || (*ey)[ind]/(*y)[ind]>0.7)     continue;
    ind++;
  }
  printf(" #points %i \n", ind);

  TGraphErrors *gr = 0;
  gr = pai::drawGraphErrors(ind, x->GetArray(),y->GetArray(),ex->GetArray(),ey->GetArray(),
  markerColor,21+markerColor,"AP", opt, "#varepsilon_{#gamma}"," #pi^{0} supp.  ",999);
  gr->GetHistogram()->SetAxisRange(xminEff,1.);
  gr->GetHistogram()->SetMinimum(1.);
  gr->GetHistogram()->SetTitleOffset(0.02,"y");
  gPad->SetLogy(1);

  return gr;
}

void suppPi0VsEffGammaNXN(const char *opt)
{
  // Dec 1, 2005
  /*
    suppPi0VsEffGammaNXN("1X1");
    suppPi0VsEffGammaNXN("2X2");
    suppPi0VsEffGammaNXN("3X3");
   */
  int P[]={1,5,10,20,30,50};
  int markerStyle[]={21,22,23,24,25,26};
  int markerColor[]={1,2,3,4,7,6};

  const int N = sizeof(P) / sizeof(int);
  TGraphErrors *grArr[N], *gr=0;;

  gROOT->SetStyle("Plain");

  double ymax=0;
  for(int i=0; i<N; i++) {
    grArr[i] = suppPi0VsEffGamma(Form("%s %iGEV",opt,P[i]));
    if(grArr[i] && ymax < grArr[i]->GetHistogram()->GetMaximum()) 
    ymax =  grArr[i]->GetHistogram()->GetMaximum();
  }

  TCanvas *c = pai::canvas("#gamma vs #pi^{0}", "GvsPi0", 15,50, 550, 600);
  TPad *p = pai::pad("#gamma efficiency vs #pi^{0} suppression", c);
  p->cd();

  TString grOpt("APC");
  TLegend *leg=0;
  for(int i=0; i<N; i++) {
    gr =  grArr[i];
    if(gr) {
      gr->SetMarkerStyle(markerStyle[i]);
      gr->SetMarkerColor(markerColor[i]);
      gr->SetLineColor(markerColor[i]);
      gr->Draw(grOpt.Data());
      if(grOpt.Contains("AP")) {
	//        gr->GetHistogram()->SetMaximum(2.*ymax);
        gr->GetHistogram()->SetMaximum(2.0e+4);
        gr->GetHistogram()->SetAxisRange(0.1,1.1);
        gr->SetTitle(opt);
        leg = new TLegend(0.62,0.68, 0.95,0.98);
        
        leg->Draw();
      }
      leg->AddEntry(gr, Form("%iGeV",P[i]), "lp");
      grOpt = "PC";
    }
  }
  gPad->SetLogy(1);
  

  c->Update();

}

void pi0ResAndEff(const char *opt, TGraphErrors **grRes, TGraphErrors **grEff)
{
  // Dec 4, 2006
  /*
    pi0ResAndEff("1X1");
    pi0ResAndEff("2X2");
    pi0ResAndEff("3X3");
   */
  double P[]={1,5,10,20,30}; // discard 50GEV
  double eP[]={0.,0.,0.,0.,0.};
  double x[6];

  TString OPT(opt);
  int N = sizeof(P) / sizeof(double);
  if(OPT.Contains("1X1")) N = 4;
  // N=6;
  TList *l[2], *lpi=0;

  TArrayD *res  = new TArrayD(N);
  TArrayD *eRes = new TArrayD(N);
  TArrayD *eff  = new TArrayD(N);
  TArrayD *eEff = new TArrayD(N);

  int np=0;
  double r,er, ef,eef;
  for(int i=0; i<N; i++) {
    getRpLists(Form("%s %iGEV",opt,int(P[i])), l);
    lpi = l[1];
    if(lpi) {
      TH1F *hpi18 = (TH1F*)lpi->At(18);
      TH1F *hpiNorm = (TH1F*)lpi->At(0);
      if(hpi18) {
        TF1 *fgi = pai::gausi("Pi0", 0.05, 0.2, hpi18);
        fgi->SetParameter(1,0.135);
        fgi->SetParameter(2,0.01);
        hpi18->Fit(fgi,"IL","", 0.05, 0.2);
        hpi18->Fit(fgi,"EM","", 0.05, 0.2);
        if(fgi->GetParameter(1)<0.16 && fgi->GetParameter(1)>0.10){
	  r  = fgi->GetParameter(2)/fgi->GetParameter(1)*100.; // in %
          er = r*(fgi->GetParError(2)/fgi->GetParameter(2)); 
	  (*res)[np]  = r;
          (*eRes)[np] = er;
          x[np] = P[i];

          int m = int(fgi->GetParameter(0));
          int n = int(hpiNorm->GetEntries());
	  pai::probAndErr(m,n,ef,eef,1);
          (*eff)[np]  = ef*100.; // to %
          (*eEff)[np] = eef*100.;// to %
          printf(" %i | P %6.1f res = %6.3f +/- %6.4f : eff %6.3f +/- %6.4f\n", 
		 np, x[np], r,er, ef, eef);
          np++;
	}
      }
    } else { // ??
    }
  }

  TCanvas *c = pai::canvas("#pi^{0} resolution and efficiency", "Pi0ResAndEff", 15,50, 550, 400);
  TPad *p = pai::pad(Form("%s, #pi^{0} resolution and efficiency", opt), c);
  p->Divide(2,1);

  //  TLegend *leg=0;
  p->cd(1);
  TGraphErrors *gr1 = pai::drawGraphErrors(np, P,res->GetArray(), eP,eRes->GetArray(),
  kBlack,24,"AP", Form("%s, #pi^{0} resolution", opt), " P(GEV)"," #pi^{0} resolution(%) ",999);
  gr1->SetLineWidth(2);
  gr1->GetHistogram()->SetMinimum(0.0);
  gr1->GetHistogram()->SetAxisRange(0.0, 50.);

  p->cd(2);
  TGraphErrors *gr2 = pai::drawGraphErrors(np, P,eff->GetArray(), eP,eEff->GetArray(),
  kBlack,24,"AP", Form("%s, #pi^{0} efficiency ", opt), " P(GEV)"," #pi^{0} efficiency (%) ",999);
  gr2->SetLineWidth(2);
  TH1 *h2 = gr2->GetHistogram();
  h2->SetMinimum(0.0);
  h2->SetMaximum(100.);
  gPad->SetLogy(1);
  //  gr2->GetHistogram()->SetAxisRange(0.0, 50.);
  
  c->Update();
  if(grRes) {
    grRes[0] = gr1;
    grEff[0] = gr2;
  }
}

void sumPi0Res()
{
  // Dec 6, 2006
  TGraphErrors *grResAll[3], *grEffAll[3], *grRes[1], *grEff[1], *gr=0;
  for(int i=1; i<=3; i++) {
    TString OPT(Form("%i%s%i",i,"X", i)); // Be careful with Form() !!
    pi0ResAndEff(OPT.Data(), grRes, grEff);
    grResAll[i-1] = grRes[0];
    grEffAll[i-1] = grEff[0];
  }

  TCanvas *c = pai::canvas("#pi^{0} resolution and efficiency vs division", "Pi0ResAndEffvsDivision", 15,50, 515,650);
  TPad *p = pai::pad("#pi^{0} resolution and efficiency vs division", c);
  p->Divide(1,2);

  p->cd(1);
  TString grOpt("APC");
  TLegend *leg=0;
  for(int i=1; i<=3; i++) {
    int ind = 3 - i;
    gr = grResAll[ind];
    gr->Draw(grOpt.Data());
    gr->SetLineWidth(2);
    gr->SetLineColor(i);
    gr->SetMarkerColor(i);
    if(grOpt=="APC"){
      gr->GetHistogram()->SetTitle("#pi^{0} resolution");
      gr->GetHistogram()->SetMaximum(20.);
      leg = new TLegend(0.6,0.1, 0.90,0.42);
      //      leg->Draw();
      grOpt = "PC";
    }
    TString OPT(Form("%i%s%i",ind+1,"X", ind+1));
    leg->AddEntry(gr, OPT.Data(), "lp");
  }

  p->cd(2);
  grOpt = "APC";
  gPad->SetLogy(1);
  for(int i=1; i<=3; i++) {
    int ind = 3 - i;
    gr = grEffAll[ind];
    gr->Draw(grOpt.Data());
    gr->SetLineWidth(2);
    gr->SetLineColor(i);
    gr->SetMarkerColor(i);
    if(grOpt=="APC"){
      gr->GetHistogram()->SetTitle("#pi^{0} efficiency");
      gr->GetHistogram()->SetMaximum(100.);
      leg = new TLegend(0.12,0.12, 0.46,0.46);
      leg->Draw();
      grOpt = "PC";
    }
    TString OPT(Form("%i%s%i",ind+1,"X", ind+1));
    leg->AddEntry(gr, OPT.Data(), "lp");
  }

  c->Update();
}

void fitHistByFunction(const char* nh,  int iloss, int p)
{
  /* 
L ~/ALICE/SHISHKEBAB/clust.C;
    Jan 10,2007 
   fitHistByFunction("hDE");
   fitHistByFunction("hDEIn");
   fitHistByFunction("29_hMip_8GEV");
   fitHistByFunction("30_hMipNcell=1_8GEV");
   fitHistByFunction("31_hMipNcell=2_8GEV");
   fitHistByFunction("32_hMipNcell>2_8GEV");
      Jan 17,2007 - 16GeV; 16GEV_ILOSS3_IHADR0
   fitHistByFunction("hDE", 3, 16);
   fitHistByFunction("29_hMip_16GEV", 3, 16);
   fitHistByFunction("30_hMipNcell=1_16GEV", 3, 16);
   fitHistByFunction("31_hMipNcell=2_16GEV", 3, 16);
      Jan 17,2007 - 32GeV; 32GEV_ILOSS3_IHADR0
   fitHistByFunction("7_hDEMip", 3, 32);
   fitHistByFunction("29_hMip_32GEV", 3, 32);
   fitHistByFunction("30_hMipNcell=1_32GEV", 3, 32);
   fitHistByFunction("31_hMipNcell=2_32GEV", 3, 32);
      Jan 18,2007 - 4GeV; 4GEV_ILOSS3_IHADR0
   fitHistByFunction("7_hDEMip", 3, 4);
   fitHistByFunction("29_hMip_4GEV", 3, 4);
   fitHistByFunction("30_hMipNcell=1_4GEV", 3, 4);
   fitHistByFunction("31_hMipNcell=2_4GEV", 3, 4);
      Jan 21,2007 - 2GeV; 2GEV_ILOSS3_IHADR0
   fitHistByFunction("7_hDEMip", 3, 2);
   fitHistByFunction("29_hMip_2GEV", 3, 2);
   fitHistByFunction("30_hMipNcell=1_2GEV", 3, 2);
   fitHistByFunction("31_hMipNcell=2_2GEV", 3, 2);
      Jan 21,2007 - 4GeV; 4GEV_ILOSS3_IHADR0_10KEV
   fitHistByFunction("7_hDEMip", 3, 4);
   fitHistByFunction("29_hMip_4GEV", 3, 4);
   fitHistByFunction("30_hMipNcell=1_4GEV", 3, 4);
   fitHistByFunction("31_hMipNcell=2_4GEV", 3, 4);
   */
  static TList *l = 0;
  if(l==0) l = readAllLists();
  TH1   *h = pai::findHistByName(l,nh);
  if(h==0) {
    printf(" <W> no histogram with name %s \n", nh);
    return;
  }

  TString snh(nh);
  // if(nf); // now only Landau convoluted with gauss
  double xmi = 0.020, xma = 0.033, width = 2.880e-4, MP = 0.02325;  
  if(snh.Contains("hDE")) {
  // Landau convoluted with gauss - Mar 09,2006; 
  //static TF1 *fitLandauWithGaus(TH1 *h, const char *addName, double xmi, double xma, double width, 
  //double MP,  double area, double gSigma, double binWidth, int setParLimits=0);
    //    xmi = 0.020; 
    xmi = 0.0210; 
    xma = 0.033; width = 2.880e-4; MP = 0.02325;
    if     (p>=16) {xmi=0.02; xma=0.04;}
    else if(p>=4)  {xmi=0.02; xma=0.03;}
    else if(p>=2)  {xmi=0.022; xma=0.032;}
    double area = h->Integral(10, h->GetNbinsX()), gSigma=1.3e-3; 
    double binWidth = h->GetBinWidth(1);
    // binWidth = 0.0;
    int setParLimits = 0;
    TF1 *f = pai::fitLandauWithGaus(h, nh, xmi, xma, width, MP,
             area, gSigma, binWidth, setParLimits);
    //    pai::addToTitle(h, Form(", ILOSS=%i, %iGeV/C",iloss,p));
    pai::addToTitle(h, Form(", ILOSS=%i, %iGeV/C, 10KeV",iloss,p));
  } else {
    xmi = 0.25; xma = 0.42; width = 3.69676e-03; MP = 2.91373e-01;
    double area = h->Integral(10, h->GetNbinsX()), gSigma=1.8e-2; 
    double binWidth = h->GetBinWidth(1);
    int setParLimits = 0;

    if     (p==4)  {xma=0.38;}
    else if(p==2)  {xmi=0.28; xma=0.44;}

    TF1 *f = pai::fitLandauWithGaus(h, nh, xmi, xma, width, MP,
             area, gSigma, binWidth, setParLimits);
    //    pai::addToTitle(h, Form(", ILOSS=%i",iloss));
    pai::addToTitle(h, Form(", ILOSS=%i, 10KeV",iloss));
  }
  gStyle->SetOptFit(111);
  h->SetAxisRange(xmi*0.9, xma*1.1);
  pai::draw(h,2);
}

void sumProtonNoHadrInteraction()
{ // Jan 12, 2007
  double sf  = 12.46;      // Sampling factor = 12.46 +/- 0.004 (from linear fit)
  double psi = 0.996*sf;   // psi in mev (8GeV)
  // double mop = 22.3183*sf; // mop in Mev from table - to low
  double mop = 23.1392*sf; // mop in Mev from Geant
  /*
     2.064 *1.032*0.16*77=26.2422 - mean energy lost from Geant
     shift   -3.250 from SIS
     mop = 26.2422 - 3.250 = 22.9922
     shift   -3.103 from langau
     mop = 26.2422 - 3.103 = 23.1392
  */
  // First point from hDE, second from rec. point with ncell=1
  /* ILOSS = 2,(D) Continuous energy loss without generation of delta-rays and full 
     Landau-Vavilov-Gauss fluctuations. In this case the variable IDRAY is forced 
     to 0 to avoid double counting of fluctuations.
  */
  static Double_t mpshift  = -0.22278298;       // Landau maximum location
  double  mopLoss2[2]   = {(22.85-mpshift*0.9837)*sf, 282.8};
  double emopLoss2[2]   = {0.02*sf,    0.4};
  double  psiLoss2[2]   = {0.9837*sf, 10.19};
  double epsiLoss2[2]   = {0.0123*sf, 0.33};
  double  gsigmaLoss2[2]= {0.0,      9.268};
  double egsigmaLoss2[2]= {0.0,      0.725};
  double  x[2]={1.,2.};
  double ex[2]={0.,0.};
  /* ILOSS = 3(1),Continuous energy loss with generation of delta-rays above 
      DCUTE (common /GCUTS/) and restricted Landau fluctuations below DCUTE.
  */
  double  mopLoss3[2]   = {23.14*sf, 290.4};
  double emopLoss3[2]   = {0.016*sf,   0.5};
  double  psiLoss3[2]   = {0.3714*sf, 2.334};
  double epsiLoss3[2]   = {0.0157*sf, 0.176};
  double  gsigmaLoss3[2]= {1.233*sf,  17.53};
  double egsigmaLoss3[2]= {0.029*sf,  0.350};

  /*  ILOSS = 4,Energy loss without fluctuation. The value obtained from 
      the tables is used directly.
   */
  double  mopLoss4[2]   = {23.11*sf, 291.6};
  double emopLoss4[2]   = {0.016*sf,   0.5};
  double  psiLoss4[2]   = {0.4723*sf, 2.268};
  double epsiLoss4[2]   = {0.0160*sf, 0.167};
  double  gsigmaLoss4[2]= {0.9879*sf, 17.00};
  double egsigmaLoss4[2]= {0.0291*sf,  0.350};

  TCanvas *c = pai::canvas("pIloss2","pIloss2", 10,20, 700, 500);  
  TPad    *pad = pai::pad("Summary of proton, no hadron interaction",c);
  pad->Divide(3,1);

  Int_t markerColor=1;
  TGraphErrors *gr1 = 0, *gr2 = 0, *gr3 = 0;
  TGraphErrors *grl3_1 = 0, *grl3_2 = 0, *grl3_3 = 0;
  TGraphErrors *grl4_1 = 0, *grl4_2 = 0, *grl4_3 = 0;
  int keyLoss4=0;

  // ILOSS = 2
  pad->cd(1);
  gr1 = pai::drawGraphErrors(2, x,mopLoss2, ex,emopLoss2,
  markerColor,21+markerColor,"AP", " MOP ", " de or e"," MeV ",999);
  gr1->GetHistogram()->SetMinimum(280.);
  gr1->GetHistogram()->SetMaximum(295.);
  gr1->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
  // ILOSS = 3
  markerColor=2;
  grl3_1 = pai::drawGraphErrors(2, x,mopLoss3, ex,emopLoss3,
  markerColor,21+markerColor,"P", "", " de or e"," MeV ",999);
  if(keyLoss4) {
  // ILOSS = 4
    markerColor=3;
    grl4_1 = pai::drawGraphErrors(2, x,mopLoss4, ex,emopLoss4,
    markerColor,21+markerColor,"P", "", " de or e"," MeV ",999);
  }
  TLegend *leg = new TLegend(0.2,0.72, 0.9,0.88);
  TLegendEntry *le1 = leg->AddEntry(gr1, "ILOSS=2", "L");
  //  printf(" text size %f  \n", le->GetTextSize());
  TLegendEntry *le2 = leg->AddEntry(grl3_1, "ILOSS=3(1)", "L");
  leg->Draw();
  le1->SetTextSize(0.1);
  le2->SetTextSize(0.1);
  le2->SetTextColor(2);
  //    TLatex *lat = 0;
  //  lat = pai::lat(funExp.Data(),20., 1., 12, 0.06, 1); 
  // mop value from BB formula - Jan 16, 2006
  TLine *mopLine = new TLine(0.93, mop, 2.07, mop);
  mopLine->Draw("same");

  // ----------------------------------------------------------------

  // ILOSS = 2
  pad->cd(2);
  markerColor=1;
  gr2 = pai::drawGraphErrors(2, x,psiLoss2, ex,epsiLoss2,
  markerColor,21+markerColor,"AP", " #xi (width) ", " de or e"," MeV ",999);
  gr2->GetHistogram()->SetMinimum(0.);
  gr2->GetHistogram()->SetMaximum(20.);
  gr2->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
  // ILOSS = 3
  markerColor=2;
  grl3_2 = pai::drawGraphErrors(2, x,psiLoss3, ex,epsiLoss3,
  markerColor,21+markerColor,"P", "", " de or e"," MeV ",999);
  // ILOSS = 4
  if(keyLoss4) {
    markerColor=3;
    grl4_3 = pai::drawGraphErrors(2, x,psiLoss4, ex,epsiLoss4,
    markerColor,21+markerColor,"P", "", " de or e"," MeV ",999);
  }
  // psi value from BB formula - Jan 16, 2006
  TLine *psiLine = new TLine(0.93, psi, 2.07, psi);
  psiLine->Draw("same");
  // ----------------------------------------------------------------

  // ILOSS = 2
  pad->cd(3);
  markerColor=1;
  gr3 = pai::drawGraphErrors(2, x,gsigmaLoss2, ex,egsigmaLoss2,
  markerColor,21+markerColor,"AP", " #sigma (resolution) ", " de or e"," MeV ",999);
  gr3->GetHistogram()->SetMinimum(0.);
  gr3->GetHistogram()->SetMaximum(20.);
  gr3->GetHistogram()->GetYaxis()->SetLabelSize(0.06);
  // ILOSS = 3
  markerColor=2;
  grl3_3 = pai::drawGraphErrors(2, x,gsigmaLoss3, ex,egsigmaLoss3,
  markerColor,21+markerColor,"P", "", " de or e"," MeV ",999);
  if(keyLoss4) {
  // ILOSS = 4
    markerColor=3;
    grl4_3 = pai::drawGraphErrors(2, x,gsigmaLoss4, ex,egsigmaLoss4,
    markerColor,21+markerColor,"P", "", " de or e"," MeV ",999);
  }
  // 
  // ----------------------------------------------------------------

  c->Update();
}

void getRpLists(const char *opt, TList *l[2])
{ 
  // Service method - Dec 1, 2006
  TString strOpt(opt), prodDir("/data/r22b/ALICE/PROD/Oct2006/"); 
  TObjArray opts;
  pai::ParseString(strOpt, opts);
  printf(" getRpLists((\"%s\") : nopt = %i\n", opt, opts.GetEntries());
  TObjString *os = (TObjString*)opts.At(0);
  TString s0(os->GetString());
  os = (TObjString*)opts.At(1);
  TString s1(os->GetString());

  l[0] = pai::readAll(Form("%s/%s/GAMMA/%s/rpHists.root", prodDir.Data(),s0.Data(),s1.Data()),
			   Form("lg_%s",s1.Data()));
  if(l[0]) l[0] = (TList*)l[0]->At(0);
  l[1] = pai::readAll(Form("%s/%s/PI0/%s/rpHists.root", prodDir.Data(),s0.Data(),s1.Data()),
			   Form("lpi0_%s",s1.Data()));
  if(l[1]) l[1] = (TList*)l[1]->At(0);
}

void testGeo()
{ // work perfect !!
  gSystem->Load("libGeom.dll");
  TGeoManager::Import("geometry.root");
  //  TGeoVolume *xen1 = gGeoManager->GetMasterVolume();
  TGeoNode *tn = gGeoManager->GetTopNode();
  int i=0;
  for(i=0; i<tn->GetNdaughters(); i++) {
    XEN1 = tn->GetDaughter(i);
    TString ns(XEN1->GetName());
    if(ns.Contains("XEN1")) break;
  }
  printf(" i %i : XEN1 %s \n", i, XEN1->GetName());
  //  TGeoNode *sm = xen1->GetNode(1);
  TGeoNodeMatrix *sm = 0;
  TGeoMatrix *mat=0;
  for(i=0; i<XEN1->GetNdaughters(); i++) {
    sm = (TGeoNodeMatrix*)XEN1->GetDaughter(i);
    printf(" i %i : SM %s : %s \n", i, sm->GetName(), sm->ClassName());
    mat = sm->GetMatrix();
    double loc1[3]={0.,0.,-171.}, loc2[3]={999.,999.,999.}, glob[3];
    mat->LocalToMaster(loc1, glob);
    mat->MasterToLocal(glob, loc2);
    for(int j=0;j<3;j++) {
      printf(" %i loc : %f %f : global %f \n", j, loc1[j], loc2[j], glob[j]);  
    }
    //    mat->Dump();
  }
}

void testAliEMCALGeometry(char *opt)
{ // Nov 7,2006
  AliEMCALGeometry *g = AliEMCALGeometry::GetInstance(opt,"");
  g->PrintGeometry();

  gROOT->GetListOfBrowsables()->Add(g);

}

/* Jun 05,2006

.L ${HOME}/ALICE/SHISHKEBAB/rHits.C+
.L ${HOME}/ALICE/SHISHKEBAB/clust.C+
//clust(0,0,"2x2")
recPoint(0,0,"2x2")
//testGeo();

*/