//////////////////////////////////////////////////////////////////////////
//                                                                      //
// ATLFast PhotonMaker class.                                         //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include <TMCParticle.h>
#include <TRandom.h>

#include "ATLFast.h"
#include "ATLFMCMaker.h"
#include "ATLFClusterMaker.h"
#include "ATLFPhotonMaker.h"
#include "ATLFCluster.h"
#include "ATLFPhoton.h"

const Float_t kPI   = TMath::Pi();
const Float_t k2PI  = 2*kPI;
const Int_t kMAXPHOT = 100;

ClassImp(ATLFPhotonMaker)

//_____________________________________________________________________________
 ATLFPhotonMaker::ATLFPhotonMaker()
{
   m_Nphotons = 0;
}

//_____________________________________________________________________________
 ATLFPhotonMaker::ATLFPhotonMaker(const char *name, const char *title)
                 :ATLFMaker(name,title)
{
   SetMinPT();
   SetMaxEta();
   SetRconeMatch();
   SetRconeSep();
   SetRconeDep();
   SetMaxEdep();

   m_Fruits     = new TClonesArray("ATLFPhoton",2, kFALSE);
   m_BranchName = "Photons";
   m_Nphotons   = 0;
   Save();
}

//_____________________________________________________________________________
 ATLFPhotonMaker::~ATLFPhotonMaker()
{
   //dummy
//printf("ATLFPhotonMaker destructor calledn");
}

//_____________________________________________________________________________
 void  ATLFPhotonMaker::AddPhoton(Int_t code, Int_t mcparticle, Int_t mother, Float_t eta, Float_t phi, Float_t pt)
{
//            Add a new photon to the list of photons

 //Note the use of the "new with placement" to create a new photon object.
 //This complex "new" works in the following way:
 //   photons[i] is the value of the pointer for photon number i in the TClonesArray
 //   if it is zero, then a new photon must be generated. This typically
 //   will happen only at the first events
 //   If it is not zero, then the already existing object is overwritten
 //   by the new photon parameters.
 // This technique should save a huge amount of time otherwise spent
 // in the operators new and delete.

   TClonesArray &photons = *(TClonesArray*)m_Fruits;
   new(photons[m_Nphotons++]) ATLFPhoton(code,mcparticle,mother,eta,phi,pt);
}

//_____________________________________________________________________________
 void ATLFPhotonMaker::Clear(Option_t *option)
{
//    Reset Photon Maker

   m_Nphotons = 0;
   ATLFMaker::Clear(option);
}


//_____________________________________________________________________________
 void ATLFPhotonMaker::Init()
{
//  Create histograms
   m_Mult          = new TH1F("PhoMult","photon multiplicity isolated",10,0,10);
   m_MultHard      = new TH1F("PhoMultHard","photon multiplicity hard",10,0,10);
   m_MultHardIsol  = new TH1F("PhoMultHardIsol","photon multiplicity hard + isol",10,0,10);
   m_E             = new TH1F("PhoE","(E-Ecru)/Ecru  vers E photon ISOLATED",50,0,100);
   m_Theta         = new TH1F("PhoTheta","(theta-thetacru) vers E photon ISOLATED",50,0,100);
   m_Counter       = new TH1F("PhoCounter","counter vers E photon ISOLATED",50,0,100);
   m_Mass2ph       = new TH1F("PhoMass2ph","pho-pho mass",50,90,110);
}

//_____________________________________________________________________________
 void ATLFPhotonMaker::Finish()
{
// Function called by ATLFast::Finish at the end of the job

   m_E->Divide(m_Counter);
   m_Theta->Divide(m_Counter);
}

//_____________________________________________________________________________
 Int_t ATLFPhotonMaker::Make()
{
//.............................................
//  This function searches for isolated photons, by scanning through
//  the list of MC particles. If a photon is found, its energy is
//  smeared using function RESPHO  and its polar angle using RESTHE.
//  Different energy-smearing for low and high luminosity can be invoked.
//  Isolated photons are stored in The TClonesArray of photons
//  and the energy clusters associated with them are removed.
//  Non-isolated photons are not considered further.
//.............................................
   Int_t phoCode[kMAXPHOT],    phoPart[kMAXPHOT],  phoMother[kMAXPHOT];
   Float_t phoEta[kMAXPHOT],   phoPhi[kMAXPHOT],   phoPT[kMAXPHOT];

   Int_t sortindex[kMAXPHOT];
   Int_t nphot  = 0;
   Int_t i, k, ic, ij, KF;

//........................................................
//... Get pointers to MCparticles arrays and ClonesArray
//........................................................
   ATLFMCMaker *mcarlo = gATLFast->MCMaker();
   TClonesArray *particles = mcarlo->Fruits();
   TMCParticle *part, *parent;
   Float_t *pPT  = mcarlo->PT();
   Float_t *pEta = mcarlo->Eta();
   Float_t *pPhi = mcarlo->Phi();
   Int_t   *Index= mcarlo->Index();
   

//........................................................
//... Get pointers to clustermaker arrays and ClonesArray
//........................................................
   ATLFClusterMaker *clustermaker = gATLFast->ClusterMaker();
   TClonesArray *clusters = clustermaker->Fruits();
   ATLFCluster *cluster;
   Float_t *cellEta   = clustermaker->CellEta();
   Float_t *cellPhi   = clustermaker->CellPhi();
   Float_t *cellSumPT = clustermaker->CellSumPT();
   Int_t nclusters    = clustermaker->Nclusters();
   Int_t ncells       = clustermaker->Ncells();


//........................................................
//....look for isolated photons, sort clusters common
//........................................................
   Bool_t isolated;
   Float_t pxphot, pyphot, pzphot, eephot;
   Float_t pz, pt, eta, phi, atheta, abseta, deta, dphi;
   Float_t adphi, dr, ddr, edep;
   Float_t enecru, thetacru, sigma1;
   Float_t ene   = 0;
   Float_t theta = 0;
   Int_t lclu;
   Float_t RconeSep2   = m_RconeSep*m_RconeSep;
   Float_t RconeDep2   = m_RconeDep*m_RconeDep;
   Float_t RconeMatch2 = m_RconeMatch*m_RconeMatch;

   for (Int_t ind=0;ind<mcarlo->Nstables();ind++) {
      i = Index[ind];
      if (i < mcarlo->Nstart()) continue;
      part = (TMCParticle*)particles->UncheckedAt(i);
      KF   = part->GetKF();
      pt     = pPT[i];
      if (pt <= 0) continue;
      if (TMath::Abs(KF) != 22) continue;
                 //-- analyse  photons
      isolated = kTRUE;
      lclu   = -1;
      pz     = part->GetPz();
      eta    = pEta[i];
      phi    = pPhi[i];
      enecru = part->GetEnergy();
      thetacru = TMath::ACos(pz/enecru);
      
                 //.... apply smearing
      Float_t sigma2, theta1, etaphot, atheta1;
      if (gATLFast->Smearing() == 1) {
          ene    = enecru;
          if (ene <= 0) continue;
          pzphot = pz;
          eephot = enecru;
          if (TMath::Abs(pzphot/eephot) < 1) {
             theta = TMath::ACos(pzphot/eephot);
          } else {
             if (pzphot > 0) theta = 0;
             if (pzphot < 0) theta = kPI;
          }
          atheta = TMath::Abs(TMath::Tan(0.5*theta));
          if (atheta < 0.0001) atheta = 0.0001;
          eta = -TMath::Log(atheta);
          sigma2 = SmearTheta(ene, eta,theta);
          theta1 = theta + sigma2;
          pzphot = ene*TMath::Cos(theta1);
          atheta1= TMath::Abs(TMath::Tan(0.5*theta1));
          if (atheta1 < 0.0001) atheta1 = 0.0001;
          etaphot= -TMath::Log(atheta1);
          sigma1 = PhotonResolution(gATLFast->Luminosity(),ene,pt,etaphot);
          pxphot = part->GetPx()*(1+sigma1);
          pyphot = part->GetPy()*(1+sigma1);
          pzphot = pzphot*(1+sigma1);
          eephot = enecru*(1+sigma1);
          if (TMath::Abs(pzphot/eephot) < 1) {
             theta = TMath::ACos(pzphot/eephot);
          } else {
             if (pzphot > 0) theta = 0;
             if (pzphot < 0) theta = kPI;
          }
          atheta = TMath::Abs(TMath::Tan(0.5*theta));
          if (atheta < 0.0001) atheta = 0.0001;
          eta  = -TMath::Log(atheta);
          phi  = mcarlo->Angle(pxphot, pyphot);
          pt   = TMath::Sqrt(pxphot*pxphot + pyphot*pyphot);
          ene  = eephot;
      }
      if (pt < m_MinPT) continue;
      abseta = TMath::Abs(eta);
      if (abseta > m_MaxEta) continue;

           //........................................................
           //....mark photon cluster
           //........................................................
      dr = 100;
      for (k=0;k<nclusters;k++) {
         cluster = (ATLFCluster*)clusters->UncheckedAt(k);
         deta    = eta - cluster->Eta();
         dphi    = phi - cluster->Phi();
         adphi   = TMath::Abs(dphi) - k2PI;
         if (TMath::Abs(dphi) > kPI) ddr = deta*deta +adphi*adphi;
         else                        ddr = deta*deta+dphi*dphi;
         if (ddr < dr) {lclu = k; dr = ddr;}
      }
      if (dr > RconeMatch2) lclu = -1;
           //........................................................
           //....check if photon isolated from clusters
           //........................................................
      for (k=0;k<nclusters;k++) {
         if (k == lclu) {
            ddr = 100;
         } else {
            cluster = (ATLFCluster*)clusters->UncheckedAt(k);
            deta    = eta - cluster->Eta();
            dphi    = phi - cluster->Phi();
            adphi   = TMath::Abs(dphi) - k2PI;
            if (TMath::Abs(dphi) > kPI) ddr = deta*deta +adphi*adphi;
            else                        ddr = deta*deta+dphi*dphi;
         }
         if (ddr < RconeSep2) isolated = kFALSE;
      }

           //........................................................
           //....check on energy deposition of cells EDEP in cone m_MaxEdep
           //........................................................
      edep = 0;
      for (ic=0;ic<ncells;ic++) {
         deta  = eta - cellEta[ic];
         dphi  = phi - cellPhi[ic];
         adphi = TMath::Abs(dphi) - k2PI;
         if (TMath::Abs(dphi) > kPI)  ddr = deta*deta + adphi*adphi;
         else                         ddr = deta*deta + dphi*dphi;
         if (ddr < RconeDep2) edep += cellSumPT[ic];
      }
      if (edep-pt > m_MaxEdep) isolated = kFALSE;

      parent = (TMCParticle*)particles->UncheckedAt(part->GetParent()-1);
      if (isolated) { // fill isolated photon structure
           //........................................................
           //....remove pho-cluster from array of clustesr
           //........................................................
         if (lclu >= 0) {nclusters--; clustermaker->RemoveCluster(lclu);}

         phoCode[nphot]    = KF;
         phoPart[nphot]    = i; 
         phoMother[nphot]  = parent->GetKF();
         phoEta[nphot]     = eta;
         phoPhi[nphot]     = phi;
         phoPT[nphot]      = pt;
         nphot++;
           //........................................................
           //.....monitor energy smearing for isolated photons
           //........................................................
         m_E->Fill(part->GetEnergy(), TMath::Abs((ene-enecru)/enecru));
         m_Theta->Fill(part->GetEnergy(),TMath::Abs(theta-thetacru));
         m_Counter->Fill(part->GetEnergy());
     }
 }  // end of loop on primary particles


    m_Mult->Fill(nphot+0.01);

//........................................................
//.....arrange isolated photons in falling E_T sequence
//........................................................
   Int_t photCode[kMAXPHOT],    photPart[kMAXPHOT],  photMother[kMAXPHOT];
   Float_t photEta[kMAXPHOT],   photPhi[kMAXPHOT],   photPT[kMAXPHOT];

   gATLFast->SortDown(nphot, phoPT, sortindex);

   for (i=0;i<nphot;i++) {
      ij = sortindex[i];
      photCode[i]   = phoCode[ij];
      photPart[i]   = phoPart[ij];
      photMother[i] = phoMother[ij];
      photEta[i]    = phoEta[ij] ;
      photPhi[i]    = phoPhi[ij];
      photPT[i]     = phoPT[ij];
   }

   Float_t pxphot1, pyphot1, pzphot1, eephot1;
   Float_t pxphot2, pyphot2, pzphot2, eephot2;
   Float_t aphotphot;

//........................................................
//......check mass distribution
//........................................................
   if (nphot == 2) {
      pxphot1 = photPT[0]*TMath::Cos(photPhi[0]);
      pyphot1 = photPT[0]*TMath::Sin(photPhi[0]);
      pzphot1 = photPT[0]*TMath::SinH(photEta[0]);
      eephot1 = photPT[0]*TMath::CosH(photEta[0]);
      pxphot2 = photPT[1]*TMath::Cos(photPhi[1]);
      pyphot2 = photPT[1]*TMath::Sin(photPhi[1]);
      pzphot2 = photPT[1]*TMath::SinH(photEta[1]);
      eephot2 = photPT[1]*TMath::CosH(photEta[1]);
      aphotphot  = (eephot1+eephot2)*(eephot1+eephot2) - (pxphot1+pxphot2)*(pxphot1+pxphot2)
                  -(pyphot1+pyphot2)*(pyphot1+pyphot2) - (pzphot1+pzphot2)*(pzphot1+pzphot2);
      if (aphotphot > 0) aphotphot = TMath::Sqrt(aphotphot);
      m_Mass2ph->Fill(aphotphot);
   }

//........................................................
//....cross-check with partons
//........................................................
   Int_t iphot    = 0;
   Int_t iphotiso = 0;
   Float_t ener;
   Int_t KF2;
   TMCParticle *part2;
   Int_t nstop = mcarlo->Nstop();
   for (i=0;i<nstop;i++) {
      part = (TMCParticle*)particles->UncheckedAt(i);
      KF   = part->GetKF();
      if (TMath::Abs(KF) == 22) {
         abseta   = TMath::Abs(pEta[i]);
         ener     = 0;
         isolated = kTRUE;
         for (k=0;k<nstop;k++) {
            part2  = (TMCParticle*)particles->UncheckedAt(k);
            KF2    = TMath::Abs(part2->GetKF());
            if ((KF2 <= 21) && (i != k) && (KF2 != 12) && (KF2 != 14) && (KF2 !=16)) {
               deta  = pEta[i] - pEta[k];
               dphi  = pPhi[i] - pPhi[k];
               adphi = TMath::Abs(dphi) - k2PI;
               if (TMath::Abs(dphi) > kPI) ddr = deta*deta + adphi*adphi;
               else                        ddr = deta*deta + dphi*dphi;
               if (ddr < RconeSep2 && pPT[k] > clustermaker->MinETCalo()) isolated = kFALSE;
               if (ddr < RconeDep2 && pPT[k] < clustermaker->MinETCalo()) ener += pPT[k];
            }
         }
         if (ener > m_MaxEdep) isolated = kFALSE;
         if (abseta < m_MaxEta && pPT[i] > m_MinPT) {
            iphot++;
            if (isolated)  iphotiso++;
         }
      }
   }
   m_MultHard->Fill(iphot+0.01);
   m_MultHardIsol->Fill(iphotiso+0.01);

//........................................................
//....Store photons in Root ClonesArray
//........................................................
// AddPhoton(Int_t code, Int_t mcparticle, Int_t mother, Int_t use, Float_t eta, Float_t phi, Float_t pt)
   m_Nphotons = 0;
   for (k=0;k<nphot;k++) {
      AddPhoton(photCode[k],
              photPart[k],
              photMother[k],
              photEta[k],
              photPhi[k],
              photPT[k]);
   }
   return 0;
}

//_____________________________________________________________________________
 Float_t ATLFPhotonMaker::PhotonResolution(Int_t keylum, Float_t ene, Float_t pt, Float_t eta)
{
//     parametrizes photon  resolution for low and high luminosity
//     parametrization from L. Poggioli and F. Gianotti
//         keylum=1   -- low  luminosity
//         keylum=2   -- high luminosity

   Float_t aa, bb, sigpu, sigph, sigph1;
   Float_t rpilup = 0;
   Float_t abseta = TMath::Abs(eta);

   while(1) {
      while(1) {
         gRandom->Rannor(aa,bb);
         sigph1 = aa*0.10/sqrt(ene); 
         if(1.+sigph1 > 0) break;
      }
      sigph=sigph1;
      if (abseta < 1.4) {
         while(1) {
            gRandom->Rannor(aa,bb);
            sigph1 =  aa*0.245/pt + bb*0.007;
            if(1.+sigph1 > 0) break;
         }
      } else {
         while(1) {
            gRandom->Rannor(aa,bb);
            sigph1 = aa* 0.306*((2.4-abseta)+0.228)/ene + bb*0.007;
            if(1.+sigph1 > 0) break;
         }
      }
      sigph = sigph + sigph1;
      if (keylum == 2) {
         if(abseta < 0.6)                 rpilup = 0.32;
         if(abseta > 0.6 && abseta < 1.4) rpilup = 0.295;
         if(abseta > 1.4)                 rpilup = 0.27;
         while(1) {
            gRandom->Rannor(aa,bb);
            sigpu = (aa*rpilup/pt);
            if(1.+sigpu > 0) break;
         }
         sigph = sigpu + sigph;
      }
      if(1.+sigph > 0) break;
   } 
   return sigph;
}

//_____________________________________________________________________________
 Float_t ATLFPhotonMaker::SmearTheta(Float_t ene, Float_t eta, Float_t theta)
{
//     parametrizes smearing in photon position  theta
//     parametrization from F. Gianotti

   const Float_t kPI = 3.141592653;
   Float_t aa, bb;
   Float_t sigph2 = 0;
   Float_t abseta = TMath::Abs(eta);
   Float_t sqrene = TMath::Sqrt(ene);

   while(1) {
      gRandom->Rannor(aa,bb);
      if(abseta < 0.8)                  sigph2 = aa*0.065/sqrene;
      if(abseta >= 0.8 && abseta < 1.4) sigph2 = aa*0.050/sqrene;
      if(abseta >= 1.4 && abseta < 2.5) sigph2 = aa*0.040/sqrene;
      if(abseta >= 2.5)                 sigph2 = 0;
      if(TMath::Abs(theta+sigph2) <= kPI) break;
   }
   return sigph2;
}

//_____________________________________________________________________________
 void ATLFPhotonMaker::PrintInfo()
{
   ATLFMaker::PrintInfo();
}