#ifndef HPARTICLEBEAMTILT_H
#define HPARTICLEBEAMTILT_H

#include "hgeomrotation.h"
#include "heventheader.h"
#include "hgeomvector.h"

#include "TMath.h"
#include "TLorentzVector.h"


using namespace std;

class HParticleBeamTilt : public TObject {
protected:
   HGeomRotation r;
   Int_t beamR;
   HEventHeader *evHeader;
   UInt_t rumId;
   Bool_t checkRunNum;
   
   Double_t Ek;
   Double_t gm2; 
   Double_t p0;
   TLorentzVector beam;
   Double_t T0;
public:
   HParticleBeamTilt(void) { clear(); }
   ~HParticleBeamTilt(void) {}
   void setEventHeader(HEventHeader *eh) { evHeader = eh; }

   void calcBeamTiltFind(Double_t alpha,Double_t beta) {
     // For titl finder only !
     checkRunNum = kFALSE;
     calcBeamTilt(alpha,beta);
   }
   
   Double_t beamMom(void) {return testBeamNum() ? p0  : 0.;}
   Double_t gamma2(void)  {return testBeamNum() ? gm2 : 0.;}
   Double_t beamEk(void)  {return testBeamNum() ? Ek  : 0.;}
   
   void setT0(Double_t t0) { T0 = t0; }
   
   const TLorentzVector *beamLVector(void) {return testBeamNum() ? &beam : NULL;}
   
   Bool_t correctAngles(Double_t theta,Double_t phi,Double_t& thetaC,Double_t& phiC) {
     // Correction of theta and phi angles for beam tilt:
     
     if( checkRunNum ) {
        if(evHeader == NULL) {
           Error("correctAngles()","Pointer to the HEventHeader == NULL!");
           exit(1);
        }
        if( !testBeamNum() ) {
           Error("correctAngles()","runId not from Feb22 beamtime");
           exit(1);
        }
     }
     
     Double_t dZ  = TMath::Cos(theta*TMath::DegToRad());
     Double_t dxy = TMath::Sqrt(1.-dZ*dZ);
     Double_t dX  = dxy*TMath::Cos(phi*TMath::DegToRad());
     Double_t dY  = dxy*TMath::Sin(phi*TMath::DegToRad());
     HGeomVector v;
     v.setXYZ(dX,dY,dZ);

     HGeomVector pr = r*v;

     Double_t pt =   TMath::Sqrt(pr(0)*pr(0)+pr(1)*pr(1));
     phiC   = TMath::ATan2(pr(1),pr(0))*TMath::RadToDeg();
     thetaC = TMath::ATan2(pt,pr(2))*TMath::RadToDeg();
     if(phiC < 0.) phiC += 360.;
     
     return kTRUE;
   }
   
private: 
   Bool_t testBeamNum(void) {
     UInt_t runN = evHeader->getEventRunNumber();
     if(runN == rumId) return kTRUE;
     rumId = runN;
     Int_t bn = 0;
     if     (rumId < 443848018) bn = 0;
     else if(rumId < 444942513) bn = 1;
     else if(rumId < 445121924) bn = 2;
     else if(rumId < 446061982) bn = 3;
     else if(rumId < 446741216) bn = 4;
     else if(rumId < 446828471) bn = 5;
     if(beamR == bn) return kTRUE;
     beamR = bn;
     if     (beamR == 1) calcBeamTilt( 0.085,0.074);  // beam I  , deg 0.085,0.074
     else if(beamR == 2) calcBeamTilt( 0.053,0.049);  // beam IIa, deg 0.053,0.049
     else if(beamR == 3) calcBeamTilt( 0.024,0.051);  // beam IIb, deg 0.024,0.051
     else if(beamR == 4) calcBeamTilt( 0.044,0.063);  // beam III, deg 0.044,0.063
     else if(beamR == 5) calcBeamTilt(-0.016,0.087);  // beam IV , deg
     else {
       clear();
       return kFALSE;
     }
      
     Double_t mp = 938.27231;
     if(T0 > 0.) Ek = T0;
     else        Ek = beamR<5 ? 4528. : 1602.;  //4500 1580;
     gm2 = (Ek+2*mp)/(2*mp);
     p0  = TMath::Sqrt((Ek+mp)*(Ek+mp) - mp*mp);
     
     HGeomVector bm(0.,0.,p0);
     HGeomRotation ri = r.inverse();
     bm = ri*bm;
     beam.SetXYZM(bm.getX(),bm.getY(),bm.getZ(),mp);
     return kTRUE;
   }
      
   void calcBeamTilt(Double_t alpha,Double_t beta) {
     alpha  *= TMath::DegToRad();
     beta   *= TMath::DegToRad();

     Double_t cosAlpha = TMath::Cos(alpha);
     Double_t sinAlpha = TMath::Sin(alpha);
     Double_t cosBeta  = TMath::Cos(beta);
     Double_t sinBeta  = TMath::Sin(beta);

     Double_t rot[9];
     rot[0] = cosBeta;
     rot[1] = 0.;
     rot[2] = sinBeta;
     rot[3] = sinAlpha*sinBeta;
     rot[4] = cosAlpha;
     rot[5] = -sinAlpha*cosBeta;
     rot[6] = -cosAlpha*sinBeta;
     rot[7] = sinAlpha;
     rot[8] = cosAlpha*cosBeta;

     r.setMatrix(rot);
     r.invert();
     return;
   }
   
   void clear(void) {
     beamR       = -1;
     evHeader    = NULL;
     checkRunNum = kTRUE;
     rumId       = 0;
     Ek          = 0.;
     gm2         = 0.;
     p0          = 0.;
     T0          = 0.;
   }
   
   ClassDef(HParticleBeamTilt,0)
};

#endif /*!HPARTICLEBEAMTILT_H*/