Well I did spend 12 weeks working on it, tbh it's quite easy to understand, you point it at data and it does the rest.
#define Class1_cxx
#include "Class1.h"
#include <TH2.h>
#include <TStyle.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TBranch.h>
#include <TH2F.h>
#include <cmath>
#include <math.h>
TH1F *h[100];
TH2F *h2[100];
void Class1::Loop()
{
cout<< "in Loop " <<endl;
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
double pi;
pi = 4 * atan(1.0);
Int_t nbytes = 0, nb = 0;
for (Long64_t jentry=0; jentry<nentries;jentry++) {
Long64_t ientry = LoadTree(jentry);
if (ientry < 0) break;
nb = fChain->GetEntry(jentry); nbytes += nb;
if(jentry%50000 == 0) cout << jentry << endl;
if(nTRK == 3) continue;
if(thrustMagAll > 0.995) continue;
if(eTotal > 8) continue;
if(nTracks > 2) continue;
if(eTotal < 5) continue;
float_t dth;
float_t theta1, theta2;
theta1 = acos(TRKcosthCM[0]);
theta2 = acos(TRKcosthCM[1]);
dth = theta1 + theta2;
dth = abs(dth);
if(dth > pi){
dth = 2*pi-dth;
}
if(dth < 2.5) continue;
float_t dphi = TRKphiCM[0] - TRKphiCM[1];
dphi = abs(dphi);
if(dphi > pi){
dphi = 2*pi-dphi;
}
if(dphi < 1.57 || dphi > 3) continue;
if(dphi > 3 && dth > 3) continue;
if(dphi > 2.9 && dth < 2.7) continue;
if(dphi > 2.9) continue;
double momentum1[3];
double momentum2[3];
double cosphi1;
cosphi1 = cos(TRKphiCM[0]);
double sinphi1;
sinphi1 = sin(TRKphiCM[0]);
double sinth1;
sinth1 = acos(TRKcosthCM[0]);
sinth1 = sin(sinth1);
double cosphi2;
cosphi2 = cos(TRKphiCM[1]);
double sinphi2;
sinphi2 = sin(TRKphiCM[1]);
double sinth2;
sinth2 = acos(TRKcosthCM[1]);
sinth2 = sin(sinth2);
momentum1[0] = TRKp3CM[0] * cosphi1 * sinth1;//x comp
momentum1[1] = TRKp3CM[0] * sinphi1 * sinth1;// y comp
momentum1[2] = TRKp3CM[0] * TRKcosthCM[0];// z comp
momentum2[0] = TRKp3CM[1] * cosphi2 * sinth2;//x comp
momentum2[1] = TRKp3CM[1] * sinphi2 * sinth2;// y comp
momentum2[2] = TRKp3CM[1] * TRKcosthCM[1];// z comp
double dmom[3];
for (int l = 0; l < 3; l++){
dmom[l] = momentum1[l] + momentum2[l];
dmom[l] = -1*dmom[l];
}
double dmom2;
dmom2 = (dmom[0]*dmom[0]) + (dmom[1]*dmom[1]) + (dmom[2]*dmom[2]);
double thetamissing;
thetamissing = dmom[2] / (pow(dmom2, 0.5));
thetamissing = acos(thetamissing);
double c, y, cprime;
c = 22;
y = thrustMagAll * 30;
cprime = y - eTotal;
if( cprime < c ) continue;
double energy1;
double energy2;
double mass;
bool echeck[2][6];
bool mucheck[2][16];
for(int j=0; j < nTRK; j++){ //loop for eSM
int eSM = eSelectorsMap[j];
for(int i=0; i < 6; i++){ //Loop for echeck
echeck[j][i] = eSM%2;
eSM = eSM/2;
}
}
for(int j=0; j < nTRK; j++){ //loop for eSM
int mSM = muSelectorsMap[j];
for(int i=0; i < 16; i++){ //Loop for echeck
mucheck[j][i] = mSM%2;
mSM = mSM/2;
}
}
//mu mu begin
if(mucheck[0][8] == 1 && mucheck[1][8] == 1){
energy1 = (TRKp3CM[0]*TRKp3CM[0]) + pow(0.105658, 2);
energy1 = pow(energy1, 0.5);
energy2 = (TRKp3CM[1]*TRKp3CM[1]) + pow(0.105658, 2);
energy2 = pow(energy2, 0.5);
mass = ((10.53 - energy1 - energy2)*(10.53 - energy1 - energy2)) - dmom2;
mass = pow(mass,0.5);
if(mass < 2) continue;
double dmomn = pow(dmom2,0.5);
h[1]->Fill(TRKp3CM[0]);
h[2]->Fill(TRKp3CM[1]);
h[3]->Fill(dmomn);
h[4]->Fill(mass);
h[5]->Fill(eTotal);
h[6]->Fill(thetamissing);
h2[2]->Fill(eTotal,thrustMagAll);
h[11]->Fill(dphi);
h[14]->Fill(dth);
h2[5]->Fill(dphi,dth);
h2[8]->Fill(dphi, eTotal);
}
//emu begin
if((mucheck[0][8] == 1 && echeck[1][4] == 1) || (echeck[0][4] == 1 && mucheck[1][8] == 1)){
if(echeck[0][4] ==1){
energy1 = (TRKp3CM[0]*TRKp3CM[0]) + pow(0.000511, 2);
energy1 = pow(energy1, 0.5);
energy2 = (TRKp3CM[1]*TRKp3CM[1]) + pow(0.105658, 2);
energy2 = pow(energy2, 0.5);
mass = ((10.53 - energy1 - energy2)*(10.53 - energy1 - energy2)) - dmom2;
mass = pow(mass,0.5);
if(mass < 2) continue;
}
if(echeck[1][4] ==1){
energy1 = (TRKp3CM[0]*TRKp3CM[0]) + pow(0.105658, 2);
energy1 = pow(energy1, 0.5);
energy2 = (TRKp3CM[1]*TRKp3CM[1]) + pow(0.000511, 2);
energy2 = pow(energy2, 0.5);
mass = ((10.53 - energy1 - energy2)*(10.53 - energy1 - energy2)) - dmom2;
mass = pow(mass,0.5);
if(mass < 2) continue;
}
h[7]->Fill(mass);
h[8]->Fill(eTotal);
h[9]->Fill(thetamissing);
h2[3]->Fill(eTotal,thrustMagAll);
h[12]->Fill(dphi);
h[15]->Fill(dth);
h2[6]->Fill(dphi,dth);
h2[9]->Fill(dphi, eTotal);
}
}
}
Theres a Script file containing all my graph plotting and the loaction of the data.
This was designed to run in ROOT, and when I did a full run it took 4 days to complete, using the most powerful computers in my department (we can access higher power ones but they are off site).