/* ecalBasic.cxx

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

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

gROOT->Reset();

#define M_PI 3.1415926

#include <iostream.h>
#include <math.h>
#include <stdlib.h>

#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;
}

void ecalBasic(){
   double   m;
   Particle pi0(0.1349766), gamma1(0.0), gamma2(0.0);
   Detector det;

// The canvas
   TCanvas *c1  = new TCanvas("c1","ECAL Simulation",10,10,600,750);

// The histograms
   TH1D *h1 = new TH1D("h1","Two-photon incalriant mass without ECAL",50, 0.0, 0.2);
   TH1D *h2 = new TH1D("h2","Two-photon incalriant mass with    ECAL",50, 0.0, 0.2);
   pad1     = new TPad("pad1","Two-photon incalriant mass without ECAL",0.05,0.50,0.95,0.95,21);
   pad2     = new TPad("pad2","Two-photon incalriant mass with    ECAL",0.05,0.05,0.95,0.45,21);
   pad1->Draw();
   pad2->Draw();

       for(int iEvent = 1; iEvent<=5000; iEvent++)
       {
               // 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);
               h1->Fill(m);
               
               // 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);
               h2->Fill(m);
       }

  pad1->cd(); h1->Draw();
  pad2->cd(); h2->Draw();

}