#include "AliPaiTestKine.h"
#include "TAliasPAI.h"

#include <assert.h>

#include <TROOT.h>
#include <TSystem.h>
#include <TMath.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TFile.h>
#include <TTree.h>
#include <TKey.h>
#include <TParticle.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TLatex.h>
#include <TDatime.h>
#include <TStyle.h>
#include <TString.h>
#include <TVector2.h>

#include "AliRun.h"
#include "AliStack.h"
#include "AliMC.h"
#include "AliHeader.h"
#include "AliGenerator.h"
#include "AliRunLoader.h"

// EMCAL staff
#include "AliEMCAL.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALHit.h"
#include "AliEMCALLoader.h"

#include "TGeant3.h"
#include "TVirtualMC.h"
// Work with output from MC - see rt.h  and rt.C; 17-mar-2003  

ClassImp(AliPaiTestKine)

  Int_t AliPaiTestKine::fDebug = 0;

Int_t nz=0, nphi=0, muMode=0;
Double_t etaMin=0., etaMax=0., phiMin=0.,phiMax=0., rEMCAL=0;
TParticle *mpart=0;

AliPaiTestKine::AliPaiTestKine(const char* name, const bool key, const int hist)
{  // you must create this object after definition of gAlice
  SetName(name);
  if(key) gROOT->GetListOfBrowsables()->Add(this);

  fGenType  = kNoGenerator; 

  fHistMode = hist;
  fL1 = fL2 = 0;
  fC1 = fC2 = fC3 = fC4 = fC5 = fC6 = 0;
  fP3 = fP4 = fP5 = fP6 = 0;

  switch (hist){
  case 1: Book();
          GetGeneratorType();
  break;
  case 2: ReadHists("RES/MU/mupGun.root");
  break;
  fHistMode = 0; // no hist
  default: printf("No histograms -> try ReadHists()\n");
  }
}

AliPaiTestKine::~AliPaiTestKine()
{ }

void AliPaiTestKine::GetGeneratorType()
{
  if(gAlice) {
     fGenName = gAlice->Generator()->GetName();
     if(fGenName.Contains("Pythia")) {
        fGenType = kPythia;
     } else if(fGenName.Contains("Hijing")) {
        fGenType = kHijing;
     } else if(fGenName.Contains("HIGINGpara")) {
        fGenType = kHijingPara;
     }
     TString name = GetName();
     name += " : ";
     name += fGenName;
     SetName(name.Data());
     printf(" Generator : %s -> %i \n", fGenName.Data(), fGenType);
  }
}

void AliPaiTestKine::Book()
{
  TH1::AddDirectory(kFALSE);
  TH2::AddDirectory(kFALSE);

  BookKine();
  BookHits();
}

void AliPaiTestKine::BookKine()
{
  //  TString smod(this->GetName());
  //  if(smod.Contains("mu"), TString::kIgnoreCase) muMode=1;
  printf("<I> AliPaiTestKine::BookKine() : muMode %i \n", muMode); 
  fL1 = new TList();
  fL1->SetName("primaryKinematics");
  Add(fL1);
  
  fL1->Add(new TH1F("hNPrim","Number of primary particles", 1000, 0., 10000.)); 
  fL1->Add(new TH2F("hPid","Geant Id of primary particles", 75,0.5,75.5, 2,0.5,2.5));
  fL1->Add(new TH1F("hPt","P_{t} of primary particles", 100, 0, 100));
  fL1->Add(new TH1F("hPz","P_{z} of primary particles", 2000, -1000., 1000.));
  fL1->Add(new TH1F("hE","Energy of primary particles", 100, 0., 100.));

  fL1->Add(new TH1F("hPhi","#phi of primary particles", 100, 0.0, 2.*TMath::Pi()));
  fL1->Add(new TH1F("hEta","#eta of primary particles", 700, -7., 7.));    // 6 histo

  fL1->Add(new TH2F("hEinEMCALKineTmp","Energy in EMCAL from kinematics(one event)", 
		    nz,etaMin,etaMax, nphi,phiMin,phiMax));                // 7 histo
  fL1->Add(new TH1F("hEinEMCALKine","Energy in EMCAL from kinematics", 400, 0., 400.));  // 8 histo
  // proy=1 - primary particles ( for checking); 2 - secondary particles
  fL1->Add(new TH2F("hRofVertexes"," Radius of vertexes ", 1000,0.,1000., 6,0.5,6.5)); // 9 
  // prox : 1 - vertexes of all primary   particles |eta|<0.7;
  //      : 2 - vertexes of all secondary particles;|eta|<0.7
  //      vertexes of particles which gave hist to EMCAL 
  //      : 3 - vertexes of tracks which deposit energy itself (mHit->GetTrack()) 
  //      : 4 - vertexes of tracks with weight eloss 
  //      : 5 - vertexes of tracks which primary look to the EMCAL with weight eloss 
  //      : 6 - vertexes of tracks which primary don't look to the EMCAL with weight eloss 
  fL1->Add(new TH2F("hTimeOfHitsEloss"," Times of hits (nanosec) with weight eloss", 
  1000, 0., 100., 3,0.5,3.5)); // 10
  //      : 1 - time distribution for all hits 
  //      : 2 - time distribution for "direct" hits
  //      : 3 - time distribution for "indirect" hits
  printf(" List %s has %i histograms : Capacity %i\n", fL1->GetName(), fL1->GetSize(), fL1->Capacity());
  fL1->ls();
}

void AliPaiTestKine::BookHits()
{
  AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
  if (!pEMCAL) {
     printf("<E> No EMCAL in file !!!\n");
     return;
  }
  printf("<I> Name of EMCAL geometry is %s \n", pEMCAL->GetTitle());
  AliEMCALGeometry *geom  = pEMCAL->GetGeometry();
  //  AliEMCALGeometry *geom  = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");

  fL2 = new TList();
  fL2->SetName("EMCAL hits");
  Add(fL2);
  if(geom) {
     printf("Constant get from GEANT geometry ");
     nz     = geom->GetNZ();
     etaMin = geom->GetArm1EtaMin();
     etaMax = geom->GetArm1EtaMax();
     nphi   = geom->GetNPhi();
     phiMin = geom->GetArm1PhiMin()*TMath::Pi()/180.;
     phiMax = geom->GetArm1PhiMax()*TMath::Pi()/180.;
     rEMCAL = geom->GetIPDistance();
     printf(" eta %4.2f->%4.2f : phi %5.3f->%5.3f : rEMCAL %6.1f\n", 
     etaMin,etaMax, phiMin,phiMax,rEMCAL);  
  }

  fL2->Add(new TH2F("hEinEMCALDeTmp","Deposit Energy in EMCAL(one event)", 
		    nz,etaMin,etaMax, nphi,phiMin,phiMax));    // 0

  Int_t nMAX = 2000;
  Double_t deMAX = 200.;
  if (muMode == 1) {
     nMAX  = 1000;
     deMAX = 0.2;
  }
  fL2->Add(new TH1F("hEinEMCALDe","Deposit Energy in EMCAL", nMAX, 0., deMAX));  // 1

  fL2->Add(new TH1F("hSF","E(kine)/DE - samlpling fraction ", 100, 0.0, 20.));   // 2
  // for mu
  fL2->Add(new TH1F("hNormDeForMu","De * sin(#theta) for #mu^{+,-}", nMAX, 0.0, deMAX)); // 3
  // for checking EMCAL geometry
  nMAX = 250;
  //  Double_t xmi = geom->GetIPDistance();
  Double_t xmi = 457;
  Double_t xma = xmi + 0.2*nMAX;  
  fL2->Add(new TH1F("hRofEmcalHits","R of EMCAL hits", nMAX,xmi,xma)); // 4
  fL2->Add(new TH2F("hRvsZofEmcalHits","R vs abs(z) of EMCAL hits",    //5 
  nMAX,xmi,xma, 400, 0., 400.)); 
  // y=1 - pi0,e+,e-; 3 - n,nbar or K0L; 3 - other (must be charge particles)
  fL2->Add(new TH2F("hDeEMCALVsGidOneEv","De in EMCAL vs geant id in event", 75,0.5,75.5, 3,0.5,3.5)) ; //6
  fL2->Add(new TH2F("hDeEMCALVsGid","De in EMCAL vs geant id in event", 75,0.5,75.5, 3,0.5,3.5)) ;      //7
  fL2->Add(new TH2F("hRelContribution","%% in EMCAL vs for el.particle, neutr. hadr, and charge", 
  100,0.0,100., 3,0.5,3.5)) ; //8
  fL2->Add(new TH2F("hEtaPhiPartHitsInEmcalYgid1","#eta:#phi for #pi0 with hits in EMCAL(de)",
		    64,-8.,8., 60, -TMath::Pi(),+TMath::Pi())); // 9
  fL2->Add(new TH2F("hEtaPhiPartHitsInEmcalYgid2","#eta:#phi for #bar{n} or K^{0}_{L} with hits in EMCAL(de)",
		    64,-8.,8., 60, -TMath::Pi(),+TMath::Pi())); // 10
  fL2->Add(new TH2F("hEtaPhiPartHitsInEmcalYgid3","#eta:#phi other with hits in EMCAL(de)",
		    64,-8.,8., 60, -TMath::Pi(),+TMath::Pi())); // 11

  printf(" List %s has %i histograms\n", fL2->GetName(), fL2->GetSize());
  fL2->ls();
}

//Bool_t LookToEMCAL(TParticle *p, Double_t x0, Double_t y0, Double_t z0)
Bool_t AliPaiTestKine::LookToEMCAL(TParticle *p)
{ // 7-feb-2003
  // x = ax*t + x0; - ax,ay,az - cos of direction ax=px/p, ay=py/p, az=pz/p
  // y = ay*t + y0;
  // z = az*t + z0;
  static Double_t x0,y0,z0;
  static Double_t ax=0., ay=0., az=0., a2=0., b=0., c=0., t=0., tt=0., x=0.,y=0.,z=0., ctmp;
  static TVector3 v;

  if(!p) return kFALSE;
  if(p->Pt()<0.2) return kFALSE; // cut on field
  x0 = p->Vx();
  y0 = p->Vy();
  z0 = p->Vz();
  if(TMath::Sqrt(x0*x0 + y0*y0) >= rEMCAL) return kFALSE; // out of radius

  ax   = p->Px()/p->P();
  ay   = p->Py()/p->P();
  az   = p->Pz()/p->P();
  //  sign = TMath::Sign(1., az); 
  a2   = ax*ax+ay*ay;

  // coefficient square equation - t**2 + b*t = c
  b  = (2.*x0*ax+2.*y0*ay)/a2;
  c  = (rEMCAL*rEMCAL - x0*x0 - y0*y0)/a2;
  // t - must be positive - 7-feb-2003 
  ctmp = (c + b*b/4.);
  if(ctmp<0.) { // just for sure
     printf(" x0 %f y0 %f ax %f ay %f :: ",x0,y0,ax,ay);
     printf(" b %f c %f ctmp %f \n", b,c, ctmp); 
     ctmp = 0.;
  }
  ctmp = TMath::Sqrt(ctmp);
  t    = ctmp - b/2.;
  
  // x,y,z to vector
  v.SetXYZ(x0+ax*t, y0+ay*t, z0+az*t);
  if(v.Phi()<phiMin || v.Phi()>phiMax) return kFALSE;
  if(v.Eta()<etaMin || v.Eta()>etaMax) return kFALSE;
  //  printf(" v.Mag() %8.3f -> rEMCAL %8.3f t %6.1f eta %5.3f phi %5.3f\n", v.Pt(),rEMCAL, t,v.Eta(),v.Phi());
  //  printf(" ax %f %f \n", ax, v.Px()/v.Mag());
  //  printf(" ay %f %f \n", ay, v.Py()/v.Mag());
  //  printf(" az %f %f \n", az, v.Pz()/v.Mag());
  return kTRUE;
}

Double_t AliPaiTestKine::RofVertex(TParticle *p)
{ 
   if(!p) return -1.;
   return TMath::Sqrt(p->Vx()*p->Vx()+p->Vy()*p->Vy()+p->Vz()*p->Vz());
}

void AliPaiTestKine::Fill()
{ 
   FillKine(fL1);
   FillHits(fL2);
}

void AliPaiTestKine::FillKine(TList* l)
{
  if(l == 0 ) l = fL1;
  if(l == 0 ) return; // 23-jun-2003

  Int_t nprim =  gAlice->Stack()->GetNprimary();
  FillH1F(l,0,double(nprim));

  TH2F* hidK  = (TH2F*)GetHist(l,7);
  if(!hidK) return;

  hidK->Reset();
  Int_t gid=0, child1=0, status;
  double phi=0., eta=0., pt=0.0, energy=0.;
  double energySum=0., rV = 0.0;
  for(Int_t i=0; i<nprim; i++){
    if(fDebug>2) {
      printf(" %i(%i) parts. : %s has %i hists (Capacity %i) \n", 
      i, nprim, l->GetName(), l->GetSize(), l->Capacity());
      l->ls();
    }
    mpart  = gAlice->Stack()->Particle(i);
     child1 = mpart->GetFirstDaughter();

     if (child1 >= 0 && child1 < nprim) continue;

     gid    = gMC->IdFromPDG(mpart->GetPdgCode()); // gMC must be defined
     if(gid <= 0) { // discard g, q qbar and .. - it is simpler
       //        printf("%s PdgCode %i : gid = %i\n", mpart->GetName(), mpart->GetPdgCode(), gid); 
        continue;
     }

     status = mpart->GetStatusCode(); 
     if(fGenType==kPythia && status != 1) {
       // May be this is just for sure - 7-feb-2003
       //        printf(" <W> somethinf wrong \n"); 
       //        printf("gid = %2i Status %i : %s \n", gid, status, mpart->GetName()); 
        continue;
     }

     pt     = mpart->Pt();
     if(pt < 0.0001) continue;

     energy = mpart->Energy();
     phi    = mpart->Phi(); // from 0 to 2pi !!! 
     eta    = mpart->Eta();

     FillH2F(l,1,double(gid), 1.); // all particles
     if(TMath::Abs(eta) < 0.7) {
        FillH2F(l,1,double(gid),2.);
        rV = RofVertex(mpart);
	//        FillH2F(l, 9,rV,1.); // 23-jun-2003 - time solution
     }

     FillH1F(l,2,pt);
     FillH1F(l,3,mpart->Pz());
     FillH1F(l,4,energy);
     FillH1F(l,5,phi);
     FillH1F(l,6,eta);

     hidK->Fill(eta,phi,energy);

     //     Print(mpart, 1);
  }
  energySum = hidK->Integral();
  FillH1F(l,8, energySum);

  Int_t ntrack = (gAlice->GetHeader())->GetNtrack();
  for(Int_t i=nprim; i<ntrack; i++) {
    mpart  = gAlice->Stack()->Particle(i);
     //     status = mpart->GetStatusCode(); 
     if(TMath::Abs(eta) < 0.7) {
        rV= RofVertex(mpart);
	//        FillH2F(l, 9,rV,2.);  // 23-jun-2003 - time solution
     }
  }
}

void AliPaiTestKine::FillHits(TList* l)
{
  if(l == 0 ) l = fL2;
  if(l == 0 ) return; // 23-jun-2003
  // see AliEMCALJetFinder::FillFromHits
  Int_t nprim = (gAlice->GetHeader())->GetNprimary();
  AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");  
  // Mar 14, 2007
  //  AliRunLoader*   rl = AliRunLoader::GetRunLoader();
  //AliEMCALLoader* emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
  // --
  TTree *treeH = 0;
  //treeH = gAlice->TreeH(); // should be change

  Int_t ntracks = (Int_t) treeH->GetEntries();
  
  TH2F* hidDe = (TH2F*)GetHist(l,0);
  hidDe->Reset();
  TH2* hidDeGid    = (TH2F*)GetHist(fL2,6);
  hidDeGid->Reset();
  TH2* hidDeGidAll = (TH2F*)GetHist(fL2,7);

  Int_t nbytes = 0, gid = 0, pdgCode=0;
  double phi=0., eta=0., yGid=0.0, peta=0., pphi=0.;
  double energySum=0.,  deSum=0., sf=0., timeH = 0.0;
  printf(" # of primaries %i : # of tracks in TreeH %i : %s \n", nprim, ntracks, treeH->GetName());
  Double_t x=0.,y=0.,z=0.,eloss=0.,r=0.,theta=0., rV=0.0;
  TParticle *mpartK    = 0; // particle from gAlice->Stack()->Particle()
  TParticle *mpDeposit = 0; // particle from gAlice->Particle(hit->GetTrack())
  Bool_t lookToEMCAL = kFALSE;
  for (Int_t track=0; track<ntracks;track++) {
     mpartK = gAlice->Stack()->Particle(track); // 31-jan-2002
     gAlice->ResetHits();              // What is this - 14-oct-2002
     nbytes += treeH->GetEvent(track);
//
//  Loop over hits
//
     AliEMCALHit* mHit=0;
     Int_t prim = 0, nh = 0;
     for (mHit=(AliEMCALHit*)pEMCAL->FirstHit(-1); mHit; mHit=(AliEMCALHit*) pEMCAL->NextHit()){
         x      =   (Double_t)mHit->X();             // x-pos of hit
	 y      =   (Double_t)mHit->Y();             // y-pos
	 z      =   (Double_t)mHit->Z();             // z-pos
	 eloss  =   (Double_t)mHit->GetEnergy();     // deposited energy
         timeH  =   (Double_t)mHit->GetTime()*1.e+9; // time in nanosec
	 //y         printf(" time %f \n", timeH);
	 r      =    TMath::Sqrt(x*x+y*y);
	 theta  =    TMath::ATan2(r,z);
	 eta    =   -TMath::Log(TMath::Tan(theta/2.));
	 phi    =    TMath::ATan2(y,x);
         hidDe->Fill(eta,phi,eloss);
         FillH1F(l,4, r);
         FillH2F(l,5, r, TMath::Abs(z));
	 if (nh==0) {
	 // Type of particle
	    mpDeposit = gAlice->Stack()->Particle(mHit->GetTrack());
            rV = RofVertex(mpDeposit);
            FillH2F(fL1, 9,rV, 3.);

            mpart     = GetPrimaryParent(mHit);
            lookToEMCAL = LookToEMCAL(mpart);
            pdgCode = mpart->GetPdgCode();
            gid     = gMC->IdFromPDG(pdgCode);
            peta    = mpart->Eta();
            pphi    = TVector2::Phi_mpi_pi(mpart->Phi());

            if(gid <0){ 
               yGid = 0.; 
            }else if(IsGammaPositronElectron(gid)){ 
               yGid=1;
            } else if(IsNeutralHadron(gid)) {
               yGid=2;
	    } else {
               yGid=3; 
            } 
         }
         FillH2F(fL1, 9,rV, 4., eloss);
         FillH2F(fL1, 10,timeH, 1.,eloss);
         if(lookToEMCAL) { // primary look to the EMCAL
            FillH2F(fL1, 9,rV, 5., eloss);
            FillH2F(fL1, 10,timeH, 2.,eloss);
         } else {
            FillH2F(fL1, 9,rV, 6., eloss);
            FillH2F(fL1, 10,timeH, 3.,eloss);
         }
	 //	 printf(" mpart %x pdg %i gid %2i : P %7.2f ETA %6.3f PHI %6.3f yGid %i\n", 
	 //         mpart, pdgCode, gid, mpart->P(), peta, pphi, Int_t(yGid));
	 //         Print(mpDeposit);

         FillH2F(l, 8+Int_t(yGid), peta, pphi, eloss); 
         hidDeGid->Fill(Double_t(gid), yGid, eloss);
         
	 //	 printf(" nh %3i eloss %f r %7.1f eta %5.3f phi %5.3f\n", nh, eloss, r, eta, phi);
         nh++;
     }
     //     if(nh>0) printf(" Track %i has %i hit(s) mpartK %x P %7.2f\n", track, nh, mpartK, mpartK->P());
     //     if(track>20) assert(0); //for testing
  }
  deSum = hidDe->Integral();
  FillH1F(l,1, deSum);

  //printf("Energy %9.3f -> de %9.3f\n",energySum, deSum); 
  if(deSum && energySum && nprim==1 && gid==1) {
     sf = energySum/deSum;
     FillH1F(l,2, sf);
  }
  if(nprim==1 && (gid==5||gid==6)) {
     theta = mpart->Theta();
     FillH1F(l,3, deSum*TMath::Sin(theta));
  }

  Double_t sum = hidDeGid->Integral(), per = 0;
  if(sum>0.0) {
     hidDeGidAll->Add(hidDeGid);
     for(Int_t ip=1; ip<=3; ip++){
        per =  100.*hidDeGid->Integral(1,hidDeGid->GetNbinsX(), ip,ip)/sum;
        FillH2F(l,8, per, Double_t(ip));
     }
  }
}

Bool_t AliPaiTestKine::IsGammaPositronElectron(const Int_t gid)
{
  if(gid==1 || gid==2 || gid==3 || gid==7) return kTRUE; // sometimes pi0 is stable
  else                 return kFALSE;
}

Bool_t AliPaiTestKine::IsNeutralHadron(const Int_t gid)
{
  //      K0L         n         nbar
  if   (gid==10 || gid==13 || gid==25) return kTRUE;
  else                                 return kFALSE;
}

Bool_t AliPaiTestKine::IsStableParticleForPythia(const Int_t gid)
{ // 4-feb-2003
  const  Int_t N=11; 
  static Int_t gidStable[N] = {1, 7,8,9, 10,11,12, 13,14,15,25};
  for(Int_t i=0; i<N; i++){
    if(gid == gidStable[i]) return kTRUE;
  }
  return kFALSE;
}

Bool_t AliPaiTestKine::IsStableParticleForHijing(const Int_t gid)
{ // 6-feb-2003 : g3->Gprint("part")
  // may be not so correct
  if(gid>=1  && gid <=16) return kTRUE; // no eta(17)
  if(gid>=18 && gid <=31) return kTRUE; // no omega(24)
  if(gid>=35 && gid <=38) return kTRUE; // D-mezons
  if(gid>=52 && gid <=55) return kTRUE; // Ds-mezons
  if(gid>=58 && gid <=71) return kTRUE; // Ds-mezons
  return kFALSE;
}

Bool_t AliPaiTestKine::IsStableParticle(const TParticle *p)
{ // 5-feb-2003 - for convenience
   if(!p || !gMC) return kFALSE;
   Int_t gid = gMC->IdFromPDG(p->GetPdgCode()); // gMC must be defined
   if(fGenType == kPythia) {
      return IsStableParticleForPythia(gid);
   } else if(fGenType == kHijing) {
      return IsStableParticleForHijing(gid);
   }
}

TParticle  *AliPaiTestKine::GetPrimaryParent(const AliEMCALHit *hit)
{ // 6-feb-2003  
   static TParticle *part=0;
   if(!hit) {
      printf("<E> AliPaiTestKine::GetPrimaryParent : input hit undefined!\n");
      return 0;
   }
   if        (fGenType == kPythia) {
      part   = gAlice->Stack()->Particle(GetPrimaryParentForPythia(hit->GetTrack()));
   } else if (fGenType == kHijing) {
      part   = gAlice->Stack()->Particle(GetPrimaryParentForHijing(hit->GetTrack())); // must be replace
   } else if (fGenType == kHijingPara) {
     // parametrized HIJING: charged and neutral pions and kaons, no decays
      part   = gAlice->Stack()->Particle(hit->GetPrimary());
   } else {
      printf("<E> AliPaiTestKine::GetPrimaryParent : wrong type of generator %i \n", fGenType);
   }
   return part;
}

Int_t  AliPaiTestKine::GetPrimaryParentForPythia(const Int_t in=0)
{ // 5-feb-2003  
  // PYTHIA: all particles except long lived (K0L) are decayed. History 
  // back to parton evolution.
  // pi0 also stable on generator level 
   static TParticle *inP = 0, *outP = 0;
   static Int_t nprimary=0, ntrack = 0, out = 0, status = 0;

   ntrack = (gAlice->GetHeader())->GetNtrack();
   if(in<=0 || in >= ntrack) return -1;   // no parent particle

   nprimary = (gAlice->GetHeader())->GetNprimary();
   inP      = gAlice->Stack()->Particle(in);
   out      = in;
   outP     = inP;
 CYCLE:
   status   = outP->GetStatusCode();
   if(status == 0) {
      out  = outP->GetFirstMother();
      outP = gAlice->Stack()->Particle(out);
      goto CYCLE;
   }
   // just for checking 
   if (out>=nprimary || IsStableParticle(outP) == kFALSE) {
      printf("<E> GetPrimaryParentForPythia() out %i > nprimary %i OR is not stable \n", out, nprimary);
      Print(outP);
   }
   return out;
}

Int_t  AliPaiTestKine::GetPrimaryParentForHijing(const Int_t in=0)
{ // 6-feb-2003  
  // HIJING:  particles with sizeable life-time are not decayes (D, B, ...). 
  // pi0 are also not decayed to save space. No parton evolution history.
   static TParticle *inP = 0, *outP = 0;
   static Int_t nprimary=0, ntrack = 0, out = 0, status = 0;

   ntrack = (gAlice->GetHeader())->GetNtrack();
   if(in<=0 || in >= ntrack) return -1;   // no parent particle

   nprimary = (gAlice->GetHeader())->GetNprimary();
   inP      = gAlice->Stack()->Particle(in);
   out      = in;
   outP     = inP;
 CYCLE:
   status   = outP->GetStatusCode();
   if(status==0 && out>=nprimary) {
      out  = outP->GetFirstMother();
      outP = gAlice->Stack()->Particle(out);
      goto CYCLE;
   }
   // just for checking 
   if (out>=nprimary || IsStableParticle(outP) == kFALSE) {
      printf("<E> GetPrimaryParentForPythia() out %i > nprimary %i OR is not stable \n", out, nprimary);
      Print(outP);
   }
   return out;
}
//Bool_t AliPsaiTestKine::

void AliPaiTestKine::SaveHists(const char* name)
{
  TString ntmp(name), dir, nf;
  if (ntmp.Contains("/")) {
     nf = ntmp; 
  }  else {
     dir = "RES/SF/";
     if (ntmp.Length()) nf  = dir + ntmp;
     else               nf  = dir + "sfTest";
  }

  if(nf.Contains(".root")==0) nf += ".root";
  printf("AliPaiTestKine::SaveHists => %s ", nf.Data());
  TFile f(nf.Data(),"RECREATE");
  if (f.IsOpen()){
     fL1->Write(fL1->GetName(), TObject::kSingleKey);
     //     fL2->Write(fL1->GetName(), TObject::kSingleKey); // I dont'k know how now 
     f.ls();
     f.Close();
  } else {
     printf(" !! <E> Problem with opening file");
  }
  printf("\n");
}

Int_t AliPaiTestKine::ReadHists(const char* name)
{
  if(strlen(name) == 0) return -1;

  printf("<I> AliPaiTestKine::ReadHists : File %s with hists will be opened\n", name); 
  TFile histFile(name, "READ");
  if(histFile.IsOpen()   == 0) return -2;
  histFile.ls();

  TKey *key=0;
  TObject *obj=0;
  TH1 *hid = 0;
  TIter next(histFile.GetListOfKeys());

  while((key=(TKey*)next())) {
     obj = key->ReadObj();
     TString tmp(obj->GetName());
     if(tmp.Contains("primaryKinematics", TString::kIgnoreCase)) {
        if (strcmp(obj->ClassName(),"TList")==0) {
	   fL1 = (TList*)obj;
           if(fL1->GetSize() == 0) return -5;
           Add(fL1);          // 4-nov-2001 
           for(Int_t i=0; i<fL1->GetSize(); i++){
	      TObject *oTmp = fL1->At(i);
              if(oTmp->InheritsFrom("TH1") == 0) continue;
              hid = (TH1*)oTmp;
              hid->SetDirectory(0);
              printf(" Hist %s : %s in list \n", hid->GetName(),hid->GetTitle());
           }

        } return -4; // this not TList
     } 
     return 0; 
  }
  return -3; // no object with this name
}

// Draw methods
void  AliPaiTestKine::DrawDeForMu()
{
  static TPad *p1=0;
  if (!fC1) {
      fC1 = new TCanvas("deForMu","Canvas , de for #mu meson",0,25,600,800);
      fC1->cd();
      TPad *padTit = new TPad("De for #mu meson ","Title Pad",0.02,0.96,0.98,1.0);
      padTit->SetFillColor(18); // See class TAttFill 
      padTit->Draw();
      padTit->cd();

      TString tit("De  for # mu meson"); 
      TDatime HistTime;
      tit += "    "; 
      tit += HistTime.AsString();
      TLatex *latex = new TLatex;
      latex->SetTextAlign(22);
      latex->SetTextSize(0.4);
      latex->DrawLatex(0.5,0.5,tit.Data());

      fC1->cd();
      p1 = new TPad("p1","Main Pad",0.0,0.00,1.00,0.96);
      p1->Draw();
  }
  p1->cd();
  p1->Clear();
  p1->Divide(1,2);
  p1->cd(1);
  DrawHist(GetHist(fL2,2), 2); // ??
  p1->cd(2);
  DrawHist(GetHist(fL1,4),2); //

  fC1->Update();
  printf("AliPaiTestKine::DrawDeForMu() finished !!\n");
}

void  AliPaiTestKine::DrawPictForHits()
{
  static TPad *p2=0;
  if (!fC2) {
      fC2 = new TCanvas("C2ForHits","C2, space hit's distribution",0,25,600,800);
      fC2->cd();
      TPad *padTit = new TPad("P2ForHits","P2",0.02,0.96,0.98,1.0);
      padTit->SetFillColor(18); // See class TAttFill 
      padTit->Draw();
      padTit->cd();

      TString tit("Space hit's distribution"); 
      TDatime HistTime;
      tit += "    "; 
      tit += HistTime.AsString();
      TLatex *latex = new TLatex;
      latex->SetTextAlign(22);
      latex->SetTextSize(0.4);
      latex->DrawLatex(0.5,0.5,tit.Data());

      fC2->cd();
      p2 = new TPad("p2","Main Pad",0.0,0.00,1.00,0.96);
      p2->Draw();
  }
  p2->cd();
  p2->Clear();
  p2->Divide(1,2);
  gStyle->SetOptStat(111111);

  p2->cd(1);
  DrawHist(GetHist(fL2,4), 2);
  GetHist(fL2,4)->SetXTitle(" R ");

  p2->cd(2);
  DrawHist(GetHist(fL2,5),2);
  GetHist(fL2,5)->SetXTitle(" R ");
  GetHist(fL2,5)->SetYTitle(" Z ");

  fC2->Update();
}

void  AliPaiTestKine::DrawDeVsGeantID()
{// 31-jan-2003
  TH1D *hidprox=0;
  TH2F *hid = 0;
  char ctmp[100];

  if(!fC3) {
     fC3 = TAliasPAI::canvas("De vs Geant ID");
     TString stmp("De vs type of hadrons : ");
     stmp += fGenName;
     stmp += " N ";
     //     stmp += (gAlice->GetHeader()->GetEvent()+1); // must be correct
     fP3 = TAliasPAI::pad(stmp.Data(), fC3);
  }
  fP3->Clear();
  fP3->Divide(1,2);

  fP3->cd(1);
  hid = (TH2F*)GetHist(fL2,7);
  Double_t scale = 100./hid->Integral();
  hid->Scale(scale);      // normalized on 100%
  TString opt("");
  Double_t max = 50.; // tune by hand now
  char* text[3] = {"#gamma, #pi^{0}, e^{+} or e^{-}"," n, #bar{n} or K^{0}_{L} "," other "};
  for(Int_t prox=1; prox<=3; prox++){
     sprintf(ctmp,"hdeVsGid%i1",prox);
     hidprox = hid->ProjectionX(ctmp, prox, prox);
     if(prox == 1) {
        hidprox->SetMaximum(max);
        hidprox->SetTitle("% of De in EMCAL vs Geant ID");
        TAliasPAI::titles(hidprox," Particle geant ID ", " %");
     }
     TAliasPAI::drawHist(hidprox, 5-prox, prox, opt.Data(), 1);

     TAliasPAI::dispLine(10., max-10-5*prox, prox, text[prox-1], 2., 1, 5-prox, 5., 0.0, 0.06);

     //     if(hidProx->GetMaximum() > max) max = hidProx->GetMaximum();
     opt = "same";
  }
  

  fP3->cd(2);
  hid = (TH2F*)GetHist(fL2,8);
  opt = "";
  for(Int_t prox=1; prox<=3; prox++){
     if(prox == 1) {
        hidprox->SetMaximum(35.);
        hidprox->SetTitle("% of De in EMCAL vs Geant ID vs event");
        TAliasPAI::titles(hidprox, " N "," % ");
     }
     sprintf(ctmp,"hPrcentVsGid%i1",prox);
     hidprox = hid->ProjectionX(ctmp, prox, prox);
     TAliasPAI::drawHist(hidprox, 5-prox, prox, opt.Data(), 1);
     opt = "same";
  }
}

void  AliPaiTestKine::DrawKine()
{// 04-feb-2003
  TH1D *hidprox=0;
  TH2F *hid  = 0;
  TH1  *hid1 = 0;
  char ctmp[100];

  if(!fC4) {
     fC4 = TAliasPAI::canvas("Initial Geant Id");
     TString stmp("Initial kinematics : ");
     stmp += fGenName;
     stmp += " N ";
     //     stmp += (gAlice->GetHeader()->GetEvent()+1); must be correct
     fP4 = TAliasPAI::pad(stmp.Data(), fC4);
  }
  fP4->Clear();
  fP4->Divide(2,2);

  fP4->cd(1);
  hid = (TH2F*)GetHist(fL1, 1);
  TString opt("");
  Int_t nev = gAlice->GetEvNumber()+1;
  char* text[2]={"All particles ","|#eta|<0.7"};
  Double_t max = 0., c[2]={0.8, 0.6};
  for(Int_t prox=1; prox<=2; prox++){
     sprintf(ctmp,"gid%i",prox);
     hidprox = hid->ProjectionX(ctmp, prox, prox);
     if(nev > 0) hidprox->Scale(1./Double_t(nev));
     fL1->Add(hidprox);

     if(prox==1) {
        max = hidprox->GetMaximum();
        TAliasPAI::titles((TH1*)hidprox, "Particle geant ID  "," % ");
        hidprox->SetMinimum(0.1);
        //hidprox->SetXTitle("Particle geant ID  ");
	//        hidprox->SetXTitle(" % ");
     }
     TAliasPAI::drawHist(hidprox, 4-prox, prox, opt.Data(), prox);
     TAliasPAI::dispLine(12., max*c[prox-1], prox, text[prox-1], 2., prox, 5-prox, 5., 0.0, 0.06);

     opt = "same";
  }
  gPad->SetLogy(1);

  fP4->cd(2);
  hid1 = GetHist(fL1,2);
  TAliasPAI::drawHist(hid1, 2);
  TAliasPAI::titles(hid1, " P_{t}", " N ");
  gPad->SetLogy(1);

  fP4->cd(3);
  hid1 = GetHist(fL1,3);
  TAliasPAI::drawHist(hid1, 2);
  TAliasPAI::titles(hid1, " P_{z}" , " N ");
  gPad->SetLogy(1);

  fP4->cd(4);
  hid1 = GetHist(fL1,6);
  TAliasPAI::drawHist(hid1, 2);
  TAliasPAI::titles(hid1 ," #eta ", " N ");
}

void  AliPaiTestKine::DrawRofVertexes(Int_t gr)
{// 07-feb-2003
  TH1D *hidprox=0;
  TH2F *hid  = 0;
  TH1  *hid1 = 0;
  char ctmp[100];
  TLatex *lat = 0;

  TString stmp(" R of vertexes ");
  if(!fC5) {
     fC5 = TAliasPAI::canvas("R of vertexes");
     if(gr == 2) {
        stmp = " Vertexes of particles with EMCAL hits ";
     }
     stmp += fGenName;
     //     stmp += " N ";
     //     stmp += (gAlice->GetHeader()->GetEvent()+1); must be correct
     fP5 = TAliasPAI::pad(stmp.Data(), fC5);
  }
  fP5->Clear();
  fP5->Divide(1, 3);

  hid = (TH2F*)GetHist(fL1, 9);
  TString opt("");
  Double_t max = 0.0, sumall = 0.0, sum = 0.0, per = 0.0;
  if(gr == 1) {
     Char_t *tit[3]={"Vertexes for primary particles","Vertexes for secondary particles",
		 "Vertexes for tracks with EMCAL hits"};
     for(Int_t prox=1; prox<=3; prox++){
        fP5->cd(prox);
        sprintf(ctmp,"hR%i",prox);
        hidprox = hid->ProjectionX(ctmp, prox, prox);
        hidprox->SetTitle(tit[prox-1]);
        fL1->Add(hidprox);

        if(hidprox->Integral()) {
           TAliasPAI::drawHist(hidprox, 2, 1, opt.Data());
           TAliasPAI::titles((TH1*)hidprox, " R ", " N ");
        }
        if(prox==1) {
           max    = hidprox->GetMaximum();
        } else {
           hidprox->SetMaximum(max);
        }
        gPad->SetLogy(1);
        hidprox->SetMinimum(0.1);
        //     opt = "same";
      }
  } else {
      fC5->SetName("Vertexes of particles with EMCAL hits ");
      opt = "";
      Char_t *tit[3]={"All particles ","Primary particles which look to EMCAL",
		 "All others"};
      for(Int_t prox=4; prox<=hid->GetNbinsY(); prox++){
        fP5->cd(prox-3);
        sprintf(ctmp,"hR%i",prox);
        hidprox = hid->ProjectionX(ctmp, prox, prox);
        hidprox->SetTitle(tit[prox-4]);
        fL1->Add(hidprox);

        if(prox==4) {
           max = hidprox->GetMaximum();
           sumall = hidprox->Integral();
        } else {
           sum = hidprox->Integral();
           per = 100.*sum/sumall;
           hidprox->SetMaximum(1.1 * max);
           if(prox == 5) sprintf(ctmp,"Direct   contribution %5.1f%% ",per);
           else          sprintf(ctmp,"Indirect contribution %5.1f%% ",per);
        }
        gPad->SetLogy(1);
        hidprox->SetMinimum(0.1);
        gStyle->SetOptStat(1000000);

        if(hidprox->Integral()) {
           TAliasPAI::drawHist(hidprox, 3, prox-3, opt.Data());
           TAliasPAI::titles((TH1*)hidprox, " R ", " De ");
        }
        if(prox!=4) {
           lat = new TLatex(400., TMath::Log10(max/3.), ctmp);
           lat->SetLineColor(prox-4);
	   lat->Draw();
        }
	//        opt = "same";
      }
  }
}

void  AliPaiTestKine::DrawTimesOfHits()
{// 07-feb-2003
  TH1D *hidprox=0;
  TH2F *hid  = 0;
  char ctmp[100];
  TLatex *lat;

  TString stmp(" Times of hits ");
  if(!fC6) {
     fC6 = TAliasPAI::canvas("cTime");
     fC6->SetTitle("Time");
     stmp += fGenName;
     fP6 = TAliasPAI::pad(stmp.Data(), fC6);
  }

  hid = (TH2F*)GetHist(fL1, 10);
  if(hid == 0) return;
  if(hid->Integral()<2) return;
 
  TString opt("");
  fP6->cd(1);
  Char_t *tit[3]={" All "," Direct hits"," Indirect hits "};
  Char_t *smod[2]={"#tau<100nsc", "#tau<50nsc"};
  Double_t sum[2], per[2], max, ymax;
  Int_t bin50;
  for(Int_t prox=1; prox<=3; prox++){
     sprintf(ctmp,"hTime%i",prox);
     hidprox = hid->ProjectionX(ctmp, prox, prox);
     hidprox->SetTitle(tit[prox-1]);
     fL1->Add(hidprox);
     bin50 = hidprox->GetXaxis()->FindBin(50.);
     if(hidprox->Integral()) {
        TAliasPAI::drawHist(hidprox, 5-prox, prox, opt.Data());
        TAliasPAI::titles((TH1*)hidprox, " time(nsec) ", " De ");
     }
     if(prox==1) {
        gPad->SetLogy(1);
        sum[0] = hidprox->Integral();
        sum[1] = hidprox->Integral(1, bin50);
        max    = hidprox->GetMaximum();
     } else {
        per[0] = 100.*hidprox->Integral()/sum[0];
        per[1] = 100.*hidprox->Integral(1,bin50)/sum[1];
        ymax   = TMath::Log10(max/3.);
        if(prox==3) ymax   = TMath::Log10(max/30.);
        for(Int_t i=0; i<2; i++){
           if     (prox == 2) sprintf(ctmp,"Direct   contribution ; %s ; %5.1f%% ", smod[i], per[i]);
           else if(prox == 3) sprintf(ctmp,"Indirect contribution ; %s ; %5.1f%% ", smod[i], per[i]);
           printf("%s\n", ctmp);
           if(i==1) ymax /= 2.; 
           lat = new TLatex(30., ymax, ctmp);
           lat->SetTextColor(prox);
	   lat->Draw();
        }
     }
     opt = "same";
  }
}

//
//// Service methods
//
void AliPaiTestKine::Print(TParticle *p,  Int_t mode)
{
  if(!p || !gMC) return;

  if(mode >=0 ) {
     Int_t pdg    = p->GetPdgCode();
     Int_t gid    = gMC->IdFromPDG(pdg); 
     Int_t status = p->GetStatusCode();
     Int_t m1     = p->GetFirstMother();
     Int_t df     = p->GetFirstDaughter();
     Int_t dl     = p->GetLastDaughter();
     printf("%s pdg %5i gid %2i status %i df %i dl %i ", p->GetName(), pdg, gid, status, df,dl);
  }

  if(mode == 1) {
     printf("Vx %f Vy %f Vz %f ", p->Vx(), p->Vy(), p->Vz());
  }
  printf("\n");
}

TH1* AliPaiTestKine::GetHist(const TList *l, const int i)
{
  static TObject *obj=0;
  if(!l) {
    printf("AliPaiTestKine::GetHist : zero TList of hists\n");
    return 0;
  }
  //  if(fDebug>3) printf("AliPaiTestKine::GetHist ## TList %s : index %i\n", l->GetName(), i);
  obj = l->At(i);
  if(obj) {
    if(obj->InheritsFrom("TH1")) return (TH1*)obj;
    else return 0;
  } else return 0; 
}

void AliPaiTestKine::FillH1F(const TList *l, const int i, const double x, const double w)
{
  static TH1  *hid=0;
  static TH1F *h1F=0;

  if(!l) return;
  hid = GetHist(l,i);

  if(!hid) return;
  h1F = (TH1F*)hid;

  if(!h1F) return;
  h1F->Fill(x,w);    
}

void AliPaiTestKine::FillH2F(const TList *l, const int i, 
const double x, const double y, const double w)
{
  static TH1  *hid=0;
  static TH2F *h2f=0;

  if(!l) return;
  hid = GetHist(l,i);

  if(!hid) return;
  h2f = (TH2F*)hid;

  if(!h2f) return;
  h2f->Fill(x,y,w);    
}

void AliPaiTestKine::DrawHist(TH1* hid,Int_t lineWidth, Int_t lineColor, const char* opt, Int_t lineStyle)
{
  if(hid) {
    if(lineWidth) hid->SetLineWidth(lineWidth);
    if(lineColor) hid->SetLineColor(lineColor);
    if(lineStyle) hid->SetLineStyle(lineStyle);
    hid->Draw(opt);
  } else {
    printf("AliPaiTestKine::DrawHist : NO hist !!\n");
  }
}

void AliPaiTestKine::CheckEMCALGeometry()
{
  // 20-dec-2002 - correct for new geometry
  /* Have to be correct - Mar 14, 2007
  AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
  if(!pEMCAL) return;
  AliEMCALGeometry *geom  = AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
  if(!geom) return;

  printf(" EMCAL geometry : NAME %s TITLE  %s \n", geom->GetName(), geom->GetTitle());
  printf(" Nlayers : EMCAL %i  HC %i \n", geom->GetNECLayers(), geom->GetNHCLayers());
  */
}

void AliPaiTestKine::Define_gMC()
{// 20-jun-2003 - for convenience
  /* Have to be correct - Mar 14, 2007
  if(gMC==0) {
    new   TGeant3("C++ Interface to Geant3");
    ((TGeant3*)gMC)->DefineParticles(); // gMC->IdFromPDG(22)
    printf("gMC is defined now \n");
  } else {
    printf("gMC is defined already \n");
  }
  */
}