/* ecalBasic.cpp

   A basic simulation of a Electromagnetic CALorimeter to see
   detector effects on the photons.

   Prints four-momenta of photons originating from a pi0 decay and 
   two-photon invariant mass of two photons with and without ECAL simulation. 
   Note that, the decay channel of pi0 is: 
   pi0 -> gamma1 + gamma2

   An example output is:
   Four-momenta of each gammas and their invariant mass:
   Before entering ECAL:
   P4( gamma1 ) = (0.472459,    0.0657655,     0.0146966,    0.477240)
   P4( gamma2 ) = (0.527541,   -0.0657655,    -0.0146966,    0.531828)
   m = 0.134977
   After  entering ECAL:
   P4( gamma1 ) = (0.445518,    0.0620154,     0.0138586,    0.450027)
   P4( gamma2 ) = (0.495100,   -0.0617213,    -0.0137929,    0.499123)
   m = 0.126978
*/

#include <iostream>
#include <cmath>
#include <cstdlib>

using namespace std;

#include "Random.h"
#include "Particle.h"
#include "Detector.h"

// returns invariant mass of particles p1 and p2
double invm2(double* p1, double* p2){

     double m2 = pow(p1[3]+p2[3],2.) -
                 pow(p1[0]+p2[0],2.) -
                 pow(p1[1]+p2[1],2.) -
                 pow(p1[2]+p2[2],2.);

     double mass = (m2>0.0) ? sqrt(m2) : 0.0;

     return mass;
}

int main(){
       double   m;
       Particle pi0(0.1349766), gamma1(0.0), gamma2(0.0);
       Detector det;
       
       for(int iEvent = 1; iEvent<50; iEvent++)
       {
               cout << "Event # " << iEvent << endl;

               // Set pi0 momentum in the Lab frame
               double px = 1.0;
               double py = 0.0;
               double pz = 0.0;
               double p  = sqrt(px*px + py*py + pz*pz);
               double e  = sqrt(p*p+pi0.m*pi0.m);
               pi0.p4(px, py, pz, e);
               
               // do a two-body decay
               pi0.twoBodyDecay(gamma1, gamma2);
               
               // boost gammas
               gamma1.boost(pi0);
               gamma2.boost(pi0);
               
               // calculate invariant mass of two gammas
               m = invm2(gamma1.p, gamma2.p);

               cout << "Four-momenta of each gammas and their invariant mass:" << endl;
               cout << "Before entering ECAL:" << endl;
               cout << "P4( gamma1 ) = "; gamma1.print();
               cout << "P4( gamma2 ) = "; gamma2.print();
               cout << "m = " << m << endl << endl;
               
               // smear four momenta of each gamma tracks
               det.ecalFast(gamma1, 0.05,0.0,2.5);
               det.ecalFast(gamma2, 0.05,0.0,2.5);

               // re-calculate invariant mass
               m = invm2(gamma1.p, gamma2.p);

               cout << "After  entering ECAL:" << endl;
               cout << "P4( gamma1 ) = "; gamma1.print();
               cout << "P4( gamma2 ) = "; gamma2.print();
               cout << "m = " << m << endl;   

               cout << "_________________________________________________________________" << endl;
       }
       
  return 0;
}