#include <iostream>
#include <TClonesArray.h>
#include <TTree.h>
#ifdef __CLING__
R__LOAD_LIBRARY(libDelphes)
#include "classes/DelphesClasses.h"
#include "external/ExRootAnalysis/ExRootTreeReader.h"
#endif
 #define PI 3.1415927

using namespace std;

void NewDelphes( const char *inputfile, const char *outfile){
  gSystem->Load("libDelphes");

  TChain chain("Delphes");
  chain.Add(inputfile);

  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();

  TClonesArray *branchParticle = treeReader->UseBranch("Particle");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");

  //new tree file
  double MET[numberOfEntries];
  double axialMET[numberOfEntries];
  double fracPT[numberOfEntries];
  double ETAl[numberOfEntries];
  double dPHI[numberOfEntries];
  double lepM[numberOfEntries];
  double lepPt[numberOfEntries];
  double ETAmiss[numberOfEntries];
  double PHImiss[numberOfEntries];
  double ETAlep[numberOfEntries];
  double PHIlep[numberOfEntries];
  double alphaT[numberOfEntries];

  TFile *newfile = new TFile(outfile, "recreate");
  TTree *tree = new TTree("T","Delphes");

  tree->Branch("lepPt",&lepPt[numberOfEntries]);
  tree->Branch("ETAlep",&ETAlep[numberOfEntries]);
  tree->Branch("PHIlep",&PHIlep[numberOfEntries]);
  tree->Branch("lepM",&lepM[numberOfEntries]);
  tree->Branch("dPHI",&dPHI[numberOfEntries]);
  tree->Branch("MissingET",&MET[numberOfEntries]);
  tree->Branch("ETAmiss",&ETAmiss[numberOfEntries]);
  tree->Branch("PHImiss",&PHImiss[numberOfEntries]);
  tree->Branch("axialMET",&axialMET[numberOfEntries]);
  tree->Branch("fracPT",&fracPT[numberOfEntries]);
  tree->Branch("alphaT",&alphaT[numberOfEntries]);


   Electron *elec1, *elec2;
   Muon *mu1, *mu2;
   MissingET *MET1, *MET2, *misset;
   GenParticle *par,*par1;


   TLorentzVector Jet1P4,Jet2P4,pho1P4,pho2P4,lepP4,lep1P4,lep2P4;

   //Loop over all events
   for(Int_t entry = 0; entry < numberOfEntries; entry++){
     //Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
   if(entry<50){   printf("Event-%d\n",entry);}


     for(Int_t i=0;i<branchParticle->GetEntries();i++){
   //Set par to be the Z-boson
         par1=(GenParticle*) branchParticle->At(i);

         for(Int_t j=0;j<branchMissingET->GetEntries();j++){
          misset  = (MissingET *) branchMissingET->At(j);

 	  double deltaPhi;
           if(fabs(misset->Phi-par1->Phi)<=PI) deltaPhi=misset->Phi-par1->Phi;
           else if(misset->Phi-par1->Phi>PI) deltaPhi=misset->Phi-par1->Phi-2*PI;
           else deltaPhi=misset->Phi-par1->Phi+2*PI;
         MET[entry] = misset->MET;
	 ETAmiss[entry] = misset->Eta;
	 PHImiss[entry] = misset->Phi;
         if ((par1->PID ==23||par1->PID==-23)) //continue;
{
         fracPT[entry] = (misset->MET - par1->PT)/(par1->PT);
         axialMET[entry] = (misset->MET) * (TMath::Cos(deltaPhi));

}
       }
   }
    if(branchElectron->GetEntries() > 1)
    {
      // Take first two electrons
      elec1 = (Electron *) branchElectron->At(0);
      elec2 = (Electron *) branchElectron->At(1);

      lepPt[entry] = (elec1->PT) + (elec2->PT);
      ETAlep[entry] = (elec1->Eta) + (elec2->Eta);
      PHIlep[entry] = (elec1->Phi) + (elec2->Phi);
      lepM[entry] = ((elec1->P4()) + (elec2->P4())).M(); //Electron Invariant Mass
      dPHI[entry] = (elec1->Phi)-(elec2->Phi);
}
    if(branchMuon->GetEntries() > 1)
    {
      // Take first two Muons
      mu1 = (Muon *) branchMuon->At(0);
      mu2 = (Muon *) branchMuon->At(1);

      lepPt[entry] = (mu1->PT) + (mu2->PT);
      ETAlep[entry] = (mu1->Eta) + (mu2->Eta);
      PHIlep[entry] = (mu1->Phi) + (mu2->Phi);
      lepM[entry] = ((mu1->P4()) + (mu2->P4())).M(); //Muon Invariant Mass
      dPHI[entry] = (mu1->Phi)-(mu2->Phi);

}

         
      if(entry<50){cout<<"AxialMET = "<<fracPT[entry]<<endl;}
  //write a new tree~~~~
  tree->Fill();
if(entry<100){cout<<"Lepton PT = "<<lepPt[entry]<<endl;}

}

  tree->Write();
  newfile->Close();
}