// Z->mumu
// use with loadmus.C and goodrun_winter03_short.txt
// looks at Run II winter (02-03) data
// Not a fully debugged code.  Please replace this with a fully debugged code.
// August 22, 2003.


#include <iostream>
#include <sstream>
#include <math.h>
#include "TChain.h"
#include "TTree.h"
#include "TFile.h"
#include "TLeaf.h"
#include "TBranch.h"
#include "TObjString.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPostScript.h"
#include "TLegend.h"
#include "TROOT.h"
#include "TGraphErrors.h"
#include "TH1F.h"
#include "TH2F.h"

#include "UCNtuple/myuber_t.hh"

#include <vector>
#include <fstream>
#include <iomanip>
#include <algorithm>

bool GoodTrackQ(mytrack_t& mytrack, int i);

void declarestructs(TTree* tree, myuber_t* uber);
int goodrun(int runno);
bool StandardMuonQ(mymuon_t& mymuon, int i, mytrack_t& mytrack, int runno);
bool LooseMuonQ(mymuon_t& mymuon, int i, mytrack_t& mytrack, int runno);
float invMass(int, int, int, Float_t et1, Float_t theta1, Float_t phi1, Float_t et2, Float_t theta2, Float_t phi2);



TH1F samesign("samesign", "Z same sign", 100, 0, 200);
TH1F zmumuUK("zmumuUK", "Z to mumu UK with corrections", 100, 0, 200);

TFile *Zfile = new TFile("zmu.root", "RECREATE"); //saves the compiled program into this file

void zsearch(TTree* oldtree) {
  //begin zsearch

  myuber_t uber;
  TTree* markchange;
  if (! oldtree) return;
  cout << "Expect to process " << oldtree->GetEntries() << " entries" << endl;
  oldtree->SetBranchStatus("*", 1);
  oldtree->GetEntry(0);
  
  markchange = oldtree->GetTree();
  
  const char* fname = '\0';
  if (oldtree->IsA() == TChain::Class())
    fname = ((TChain*) oldtree)->GetFile()->GetPath();
  //cout<<"fname"<<endl;
  
#ifdef __CINT__
  gInterpreter->LoadMacro("declarestructs.C");
#endif
  
  declarestructs(oldtree, &uber);
  zmumuUK.Reset();
  samesign.Reset();
  int nmass = 0;
  int nmassUK = 0; 
  int nz = 0;

 TH1F *zmumu = new TH1F("zmumu", "Z to mumu with corrections", 100, 0, 200);

 TH1F *ptN1 = new TH1F("ptN1", "Pt of Z candiate muons (GeV)", 100, 20, 120);
 TH1F *z0N1 = new TH1F("z0N1", "Z0 of Z candiate muons (cm)", 100, -70, 70);
 TH1F *d0N1 = new TH1F("d0N1", "D0 of Z candiate muons (cm)", 50, 0, .2);
 TH1F *isoN1 = new TH1F("isoN1", "Iso/pt of Z candiate muons", 100, 0, .11);
 TH1F *trkN1 = new TH1F("trkN1", "Nhits of Z candiate muons", 200, 0, 200);
 TH1F *emN1 = new TH1F("emN1", "EM Energy of Z candiate muons (GeV)", 100, 0, 2.0);
 TH1F *hadeN1 = new TH1F("hadeN1", "HAD Energy of Z candiate muons (GeV)", 100, 0, 6);
 TH1F *deltaZ = new TH1F("deltaZ", "Delta Z of Z candiate muons (cm)", 50, -5,5);
 TH1F *cmudx = new TH1F("cmudx", "CMU delta x of Z candiate muons (cm)", 100,-4,4);
 TH1F *cmupdx = new TH1F("cmupdx", "CMP delta x of Z candiate muons (cm)", 100,-7,7);
 TH1F *cmudz = new TH1F("cmudz", "CMU delta z of Z candiate muons (cm)", 100,-12,12);
 TH1F *etaN1 = new TH1F("etaN1", "Eta of Z candiate muons", 100,-1.5,1.5);	
 TH1F *phi0N1 = new TH1F("phi0N1", "BC Phi0 of Z candiate muons", 100,0,6.3);
 TH1F *cmupOnlydx = new TH1F("cmupOnlydx", "CMU delta x of Z candiate CMUP muons (cm)", 100,-5,5);
 TH1F *cmxOnlydx = new TH1F("cmxOnlydx", "CMP delta x of Z candiate CMX muons (cm)", 100,-7,7);
 TH1F *cmuOnlydz = new TH1F("cmuOnlydz", "CMU delta z of Z candiate CMUP muons (cm)", 100,-12,12);
 TH1F *cmuOnlydx = new TH1F("cmuOnlydx", "CMU delta x of Z candiate CMUP muons (cm)", 100,-12,12);
 
// Number of Muon Hits, Run Number, and Event Number Histograms
 TH1F *numhits = new TH1F("numhits", "Number of Muon Hits", 97,0,97);
 TH1F *runnumber = new TH1F("runnumber", "The Run Number", 100,140000,160000);
 TH1F *eventnum = new TH1F("eventnum", "The Event Number", 100,0,10000000);

// to plot histograms of run number and event number with all the data:
 TH1F *runnumall = new TH1F("runnumall", "The Run Number of all the Data", 100,140000,160000);
 TH1F *eventnumall = new TH1F("eventnumall", "The Event Number of all the Data", 100,0,10000000);


// forward backward muon charge asymmetry
 TH1F *Afb = new TH1F("Afb", "Forward Backward Muon Asymmetry Versus Z Mass", 7, 60, 200);
 TH1F *PEta = new TH1F("PEta", "Number of Positive Muons Versus Eta", 30, -1.5, 1.5);
 TH1F *NEta = new TH1F("NEta", "Number of Negative Muons Versus Eta", 30, -1.5, 1.5);


// DEFINING VARIABLES
// finding the first and last run numbers
   int r = 200000;
   int r2 = 0; 

// finding the run and event numbers of the highest mass event
   
   int highmass=0;
   int highrun=0;
   int highevent=0;


// printing the values onto the screen
   int nn = 1;

// Counting the number of runs and events
   int runno_new = 0;
   int numruns = 0;
   int eventno_new = 0;
   int numevents = 0;

// forward backward muon asymmetry
   Float_t F1 = 0;
   Float_t B1 = 0;
   Float_t F2 = 0;
   Float_t B2 = 0;
   Float_t F3 = 0;
   Float_t B3 = 0;
   Float_t A1 = 0;
   Float_t A2 = 0;
   Float_t A3 = 0;
   Float_t A[100];
   Float_t F[100];
   Float_t B[100];


   for (int i=0; i<100; i++)
   {
     A[i] = 0.0;
     F[i] = 0.0;
     B[i] = 0.0;
   }
   


  for (int i = 0; i < oldtree->GetEntries(); i++) {
    // only way I can get this to be reliable
    oldtree->GetEntry(i);
    declarestructs(oldtree->GetTree(), &uber);
    oldtree->GetEntry(i);
    

// to plot histograms of run number and event number with all the data   
     runnumall->Fill(uber.header.run);
     eventnumall->Fill(uber.header.event);

// finding the first and last run numbers
   int runno = uber.header.run;

   if (runno < r)
   {
      r = runno;
   }
   
   if (runno > r2)
   {
      r2 = runno;
   }
   


// Counting the number of runs and events
   int eventno = uber.header.event;
   if (runno_new != runno)
   {
      runno_new = runno;
      numruns++;
   }   

   if (eventno_new != eventno)
   {
      eventno_new = eventno;
      numevents++;
   }





    if (i % 1000 == 0) {
      cout << "entry: " << i << " run " << uber.header.run 
	   << " event " << uber.header.event << endl;
    }
    
    
    // Remove bad_runs in data
    //    if (goodrun(uber.header.run)==0) continue;  // this generate segmentation fault
    //                                                   why? what is the difference with below method
    
    Int_t igood = goodrun(uber.header.run);
    if (igood==0) continue;
    
    // check CMUP18 trigger (Warning: there is "MUON_CMUP18_VOLUNTEER", we shoudn't take this)
    // String name, "MUON_CMUP18" will pick up volunteer trigger too.
    // So additional L1 trigger is imposed not to pick up volunteer trigger.
    bool passes_L3trigger = false;
    for (int k = 0; k < uber.trigger.nl3trigs; k++) {
      passes_L3trigger |= (strstr(uber.trigger.l3trigs[k], "MUON_CMUP18_v") != NULL);
    }
    if(!passes_L3trigger) continue;
    
    bool passes_L1trigger = false;
    for (int k = 0; k < uber.trigger.nl1trigs; k++) {
      passes_L1trigger |= (strstr(uber.trigger.l1trigs[k], "L1_CMUP6_PT4_v") != NULL);
    }
    if(!passes_L1trigger) continue;

    int nsecondary = 0;
    int muind[MAXMUO];
    int m[MAXMUO];
    int primInd[MAXMUO];
    int secInd[MAXMUO];
    
    Double_t pt[MAXMUO];
    Double_t p[MAXMUO];
    Double_t theta[MAXMUO];
    
    
    Int_t nprimary = 0;
 
    if (uber.muon.nmuon > 0) {
      for (int j = 0; j < uber.muon.nmuon; j++) {
	if((uber.muon.cosmic&0x1)==0x1) break;  //if cosmic, exit
       
	bool thismupasses = false;
        m[j] = uber.muon.trind[j];
	pt[j] = uber.muon.q[j]/(uber.muon.q[j]/uber.track.bcpt[m[j]]-.00042-.00116*sin(uber.track.bcphi0[m[j]]+0.3));
	Float_t angle = uber.track.bccottheta[m[j]];
	theta[j] = atan(1/angle);
	p[j] = pt[j]/sin(theta[j]);
 
	if(StandardMuonQ(uber.muon, j, uber.track, uber.header.run)){
	  primInd[nprimary] = j;
	  nprimary++;
	}
	else if(LooseMuonQ(uber.muon, j, uber.track, uber.header.run)){
	  secInd[nsecondary] = j;
	  nsecondary++;
	}
    
      }

	if (nprimary + nsecondary > 1) {  //nmu + nel >= 2;
	  if(nprimary > 1){
	    muind[0] = primInd[0];
	    muind[1] = primInd[1];
	  }
	  if(nprimary == 1){
	    muind[0] = primInd[0];
	    muind[1] = secInd[0];
	  }
	  if(nprimary == 0)continue;
	  
	  Float_t pt0 = pt[muind[0]];
	  Float_t pt1 = pt[muind[1]];
	  Float_t theta0 =  2.*atan2(exp(-(uber.muon.eveta[muind[0]])),1);
	  Float_t theta1 =  2.*atan2(exp(-(uber.muon.eveta[muind[1]])),1);
	  Float_t phi0 = uber.muon.phi0[muind[0]];
	  Float_t phi1 = uber.muon.phi0[muind[1]];
	  Int_t trk0 = uber.muon.trind[muind[0]];
	  Int_t trk1 = uber.muon.trind[muind[1]];
	  



	  if(uber.track.bcpt[trk0]!=0.) {
	    pt0 = uber.muon.q[muind[0]]/( uber.muon.q[muind[0]]/uber.track.bcpt[trk0] - .00042 - .00116*sin(uber.track.bcphi0[trk0]+0.3 ) );
	    //	    theta0 = atan(1/uber.track.bccottheta[trk0]);
	    theta0 =  2.*atan2(exp(-(uber.track.bccottheta[trk0])),1); 
	    phi0 = uber.track.bcphi0[trk0];
	  } 
	  if(uber.track.bcpt[trk1]!=0.) {
	    pt1 = uber.muon.q[muind[1]]/( uber.muon.q[muind[1]]/uber.track.bcpt[trk1] - .00042 - .00116*sin(uber.track.bcphi0[trk1]+0.3 ) );
	    //   theta1 = atan(1/uber.track.bccottheta[trk1]); 
	    theta1 =  2.*atan2(exp(-(uber.track.bccottheta[trk1])),1); 
	    phi1 = uber.track.bcphi0[trk1];
	  } 



	  Float_t invmass = sqrt(2*pt0*pt1*(1/(sin(theta0)*sin(theta1)) - cos(phi0)*cos(phi1)-sin(phi0)*sin(phi1)-1/tan(theta1)*1/tan(theta0)));
  



// Calculating the invariant mass of the muons using muind
	  Float_t px0 = uber.muon.px[muind[0]];
	  Float_t px1 = uber.muon.px[muind[1]];
	  Float_t py0 = uber.muon.py[muind[0]];
	  Float_t py1 = uber.muon.py[muind[1]];
	  Float_t pz0 = uber.muon.pz[muind[0]];
	  Float_t pz1 = uber.muon.pz[muind[1]];
	  Float_t e0 = uber.muon.e[muind[0]];
	  Float_t e1 = uber.muon.e[muind[1]];



/*  Calculating the invariant mass of the muons using beam constraint bcpt, bccottheta, and bcphi0 (using Lauren's way of calculating these) */

	

// finding the run and event numbers of the highest mass event
   
   if (highmass < invmass)
   {
      highmass = invmass;
      highrun = uber.header.run;
      highevent = uber.header.event;
   }


   if (invmass > 200)
	cout << "   Run " << uber.header.run << " Event " << uber.header.event << " has a mass of " << invmass << endl;


// looking at run number 152598 event number 298854 which has two positive muons
  if (uber.header.event == 298854)
    {
     cout << " run: " << uber.header.run << " event: " << uber.header.event << endl;
     cout << " The charge of muon 0 is " << uber.muon.q[muind[0]] << endl;
     cout << " The charge of muon 1 is " << uber.muon.q[muind[1]] << endl;

     cout << " The eta of muon 0 is " << uber.muon.eveta[muind[0]] << endl;
     cout << " The eta of muon 1 is " << uber.muon.eveta[muind[1]] << endl;

     cout << " The pt of muon 0 is " << uber.muon.pt[muind[0]] << endl;
     cout << " The pt of muon 1 is " << uber.muon.pt[muind[1]] << endl;
     cout << " The invariant mass is " << invmass << endl;
    }
     




if(uber.muon.q[muind[0]]*uber.muon.q[muind[1]]< 0  ){


// To plot asymmetry versus Z mass
      for(int i=1; i<=7; i++)
	{
	if((invmass>=Afb->GetBinLowEdge(i))&&(invmass<Afb->GetBinLowEdge(i+1)))
	{
	   if((uber.muon.q[muind[0]]==-1)&&(uber.muon.eveta[muind[0]]>uber.muon.eveta[muind[1]]))
	   F[i-1]++;
	   if((uber.muon.q[muind[1]]==-1)&&(uber.muon.eveta[muind[0]]>uber.muon.eveta[muind[1]]))
	   B[i-1]++;
	   if((uber.muon.q[muind[0]]==-1)&&(uber.muon.eveta[muind[0]]<uber.muon.eveta[muind[1]]))
	   B[i-1]++;
	   if((uber.muon.q[muind[1]]==-1)&&(uber.muon.eveta[muind[0]]<uber.muon.eveta[muind[1]]))
	   F[i-1]++;
	}
	if ((F[i-1]+B[i-1]) != 0)
          A[i-1]=((F[i-1]-B[i-1])/(F[i-1]+B[i-1]));
	Afb->SetBinContent(i,A[i-1]);
        }


// To plot the number of positive and negative muons versus Eta

	 if(uber.muon.q[muind[0]]==1)
	     PEta->Fill(uber.muon.eveta[muind[0]]);
	
	 if(uber.muon.q[muind[1]]==1)
	     PEta->Fill(uber.muon.eveta[muind[1]]);
	
	 if(uber.muon.q[muind[0]]==-1)
	     NEta->Fill(uber.muon.eveta[muind[0]]);
	
	 if(uber.muon.q[muind[1]]==-1)
	     NEta->Fill(uber.muon.eveta[muind[1]]);
	

	    ptN1->Fill(pt0);
	    ptN1->Fill(pt1);
	    emN1->Fill(uber.muon.em[muind[0]]);
	    hadeN1->Fill(uber.muon.had[muind[0]]);
	    isoN1->Fill(uber.muon.EtCone4[muind[0]]/pt0);
	    d0N1->Fill(fabs(uber.muon.d0[muind[0]]));
	    z0N1->Fill(uber.muon.ztrack[muind[0]]);
	    deltaZ->Fill(uber.muon.ztrack[muind[0]]-uber.muon.ztrack[muind[1]]);
	    emN1->Fill(uber.muon.em[muind[1]]);
	    hadeN1->Fill(uber.muon.had[muind[1]]);
	    isoN1->Fill(uber.muon.EtCone4[muind[1]]/pt1);
	    d0N1->Fill(fabs(uber.muon.d0[muind[1]]));
	    z0N1->Fill(uber.muon.ztrack[muind[1]]);
	    cmupdx->Fill(uber.muon.CMPDelX[muind[0]]);
	    cmudx->Fill(uber.muon.CMUDelX[muind[0]]);
	    cmudz->Fill(uber.muon.CMUDelZ[muind[0]]);
	    cmupdx->Fill(uber.muon.CMPDelX[muind[1]]);
	    cmudx->Fill(uber.muon.CMUDelX[muind[1]]);
	    cmudz->Fill(uber.muon.CMUDelZ[muind[1]]);
	    phi0N1->Fill(uber.track.bcphi0[m[muind[0]]]);
	    phi0N1->Fill(uber.track.bcphi0[m[muind[1]]]);
	    etaN1->Fill(uber.muon.eveta[muind[0]]);
	    etaN1->Fill(uber.muon.eveta[muind[1]]);
	    trkN1->Fill(uber.muon.naxhits[muind[1]]+uber.muon.nsthits[muind[1]]+uber.muon.naxhits[muind[0]]+uber.muon.nsthits[muind[0]]);	
            
// Number of Muon Hits, Run Number, and Event Number Histograms
	    numhits->Fill(uber.muon.naxhits[muind[0]]+uber.muon.nsthits[muind[0]]);
	    numhits->Fill(uber.muon.naxhits[muind[1]]+uber.muon.nsthits[muind[1]]);
            runnumber->Fill(uber.header.run);
	    eventnum->Fill(uber.header.event);

	    zmumu->Fill(invmass);
	    nz++;

	    bool CMUP0 = 
	      ( (uber.muon.det[muind[0]]&DET_CMU)&&(uber.muon.det[muind[0]]&DET_CMP) &&
		(fabs(uber.muon.CMUDelX[muind[0]]) < 3.0) && // units are cm
		(fabs(uber.muon.CMPDelX[muind[0]]) < 5.0) ); // units are cm
	    
	    bool CMU_or_CMNP0 = 
	      ( (uber.muon.det[muind[0]]&DET_CMU)&&(!(uber.muon.det[muind[0]]&DET_CMP)) &&
		(fabs(uber.muon.CMUDelX[muind[0]]) < 3.0) );
	    
	    bool CMP0 = 
	      ( (uber.muon.det[muind[0]]&DET_CMP)&&(!(uber.muon.det[muind[0]]&DET_CMU)) &&   
		(fabs(uber.muon.CMPDelX[muind[0]]) < 5.0) );
	    
	    bool CMX0 = 
	      ( (uber.muon.det[muind[0]]&DET_CMX) &&    
		(fabs(uber.muon.CMXDelX[muind[0]]) < 6.0) );

	    bool CMUP1 = 
	      ( (uber.muon.det[muind[1]]&DET_CMU)&&(uber.muon.det[muind[1]]&DET_CMP) &&
		(fabs(uber.muon.CMUDelX[muind[1]]) < 3.0) && // units are cm
		(fabs(uber.muon.CMPDelX[muind[1]]) < 5.0) ); // units are cm
	    
	    bool CMU_or_CMNP1 = 
	      ( (uber.muon.det[muind[1]]&DET_CMU)&&(!(uber.muon.det[muind[1]]&DET_CMP)) &&
		(fabs(uber.muon.CMUDelX[muind[1]]) < 3.0) );
	    
	    bool CMP1 = 
	      ( (uber.muon.det[muind[1]]&DET_CMP)&&(!(uber.muon.det[muind[1]]&DET_CMU)) &&   
		(fabs(uber.muon.CMPDelX[muind[1]]) < 5.0) );
	    
	    bool CMX1 = 
	      ( (uber.muon.det[muind[1]]&DET_CMX) &&    
		(fabs(uber.muon.CMXDelX[muind[1]]) < 6.0) );

	    if(CMUP0)cmupOnlydx->Fill(uber.muon.CMPDelX[muind[0]]);
	    if(CMUP1)cmupOnlydx->Fill(uber.muon.CMPDelX[muind[1]]);

	    if(CMX0)cmxOnlydx->Fill(uber.muon.CMXDelX[muind[0]]);
	    if(CMX1)cmxOnlydx->Fill(uber.muon.CMXDelX[muind[1]]);
	    

	    if(CMUP0)cmuOnlydz->Fill(uber.muon.CMUDelZ[muind[0]]);
	    if(CMUP1)cmuOnlydz->Fill(uber.muon.CMUDelZ[muind[1]]);

	    if(CMUP0)cmuOnlydx->Fill(uber.muon.CMUDelX[muind[0]]);
	    if(CMUP1)cmuOnlydx->Fill(uber.muon.CMUDelX[muind[1]]);



	  }
	  else{samesign.Fill(invmass);
	  }
	  if(66 < invmass && invmass < 116)nmass++;
	  Float_t massll = invMass(nn,runno,eventno,pt0,theta0,phi0,pt1,theta1,phi1);
	  zmumuUK.Fill(massll);
	  if(66 < massll && massll < 116)nmassUK++;
	}
	
      }
    }
  Zfile->Write();
  cout<<"Events: "<<nz<<endl;
  cout<<"events in mass range: "<<nmass<<endl;
  cout<<"UK events in mass range: "<<nmassUK<<endl;



// Counting the number of runs and events
  cout << "   The total number of runs: " << numruns << endl;
  cout << "   The total number of events: " << numevents << endl;


// finding the first and last run numbers
  cout << "   The first run number:  "<<r<<endl;
  cout << "   The last run number:  "<<r2<<endl;


// forward backward muon asymmetry
  cout << "   The number of forward events for Z's with a mass >60 and <80: " << F[0] << endl;
  cout << "   The number of backward events for Z's with a mass >60 and <80: " << B[0] << endl;
  cout << "   The number of forward events for Z's with a mass >80 and <100: " << F[1] << endl;
  cout << "   The number of backward events for Z's with a mass >80 and <100: " << B[1] << endl;
  cout << "   The number of forward events for Z's with a mass >100: " << F[2] << endl;
  cout << "   The number of backward events for Z's with a mass >100: " << B[2] << endl;

  cout << "   The asymmetry for Z's with a mass >60 and <80: " << A[0] << endl;
  cout << "   The asymmetry for Z's with a mass >80 and <100: " << A[1] << endl;
  cout << "   The asymmetry for Z's with a mass >100: " << A[2] << endl;




// finding the run and event numbers of the highest mass event

  cout << "   The event with the highest mass has a run number: " << highrun << ", eventnum: " << highevent << ", and a mass of: " << highmass << endl;
  
}
  

bool GoodTrackQ(mytrack_t& mytrack, int i)
{
  if(i<0) return(false);
 
  bool ans = false;
 
  Int_t num_ax=0, num_axhits=0;
  Int_t num_st=0, num_sthits=0;
  for (int j = 0; j < 8; j++) {
    Int_t num_hit_sl = (mytrack.slhitcount[i]>>j*4)&0xf;
    
    if (fmod(j*1.0,2)>0) {
      num_axhits +=num_hit_sl; 
    } else {
      num_sthits +=num_hit_sl; 
    }
    
    if(num_hit_sl<7) continue;
    
    if (fmod(j*1.0,2)>0) {
      ++num_ax;}
    else {
      ++num_st;
    }
  }  
  ans = (num_ax>=3&&num_st>=3);
  ans &= (fabs(mytrack.z0[i])<60.);
 
  
  return(ans);
}

int goodrun(int runno) {

  static std::vector<int> goodrunlist;

  if (goodrunlist.empty()) {
    //    std::ifstream fin("goodruns");
    std::ifstream fin("goodrun_winter03_short.txt");
    int run;
    while (fin >> run) {
      goodrunlist.push_back(run);
      cout << " good " << run << endl;
    }

    fin.close();
    std::sort(goodrunlist.begin(), goodrunlist.end());
  }

  int good = std::binary_search(goodrunlist.begin(), goodrunlist.end(), runno);

  return good;
}


Float_t invMass(int nn, int runno, int eventno, Float_t et1, Float_t theta1, Float_t phi1, Float_t et2, Float_t theta2, Float_t phi2)
{

  Float_t E_1 = et1/sin(theta1);
  Float_t Ex_1 = E_1*sin(theta1)*cos(phi1);
  Float_t Ey_1 = E_1*sin(theta1)*sin(phi1);
  Float_t Ez_1 = E_1*cos(theta1);

  Float_t E_2 = et2/sin(theta2);
  Float_t Ex_2 = E_2*sin(theta2)*cos(phi2);
  Float_t Ey_2 = E_2*sin(theta2)*sin(phi2);
  Float_t Ez_2 = E_2*cos(theta2);

  /*
  cout << " mass " << E_1 
       << "  "     << E_2 << endl;
  */

    if(theta1==0.0||theta2==0.0) {
      cout << "Warning in theta value " << endl;
    }

 Float_t check = pow((E_1+E_2),2) - pow((Ex_1+Ex_2),2) 
                 - pow((Ey_1+Ey_2),2) - pow((Ez_1+Ez_2),2);    
 /*
 cout << " mass_check " << Ex_1
      << setw(10) << Ez_1
      << setw(10) << Ex_2
      << setw(10) << Ez_2
      << setw(10) << check << endl;
  */
  Float_t mass  = sqrt(pow((E_1+E_2),2) - pow((Ex_1+Ex_2),2) 
		     - pow((Ey_1+Ey_2),2) - pow((Ez_1+Ez_2),2));
  
  return mass;

}


//_____________________________________________________________________________
bool StandardMuonQ(mymuon_t& mymuon, int i, mytrack_t& mytrack, int runno)
{
  //need to check on id for mytrack


  bool ans = false;

  bool energyFlow = (mymuon.e[i]>100.) ?
       	 ( (mymuon.em[i] < 2.0 + 0.0115*(mymuon.e[i]-100)) 
       && (mymuon.had[i] < 6.0 + 0.0280*(mymuon.e[i]-100)) ):
       	  (mymuon.em[i] < 2.0) && (mymuon.had[i] < 6.0);

  Float_t pt;

  if(mytrack.bcpt[mymuon.trind[i]]!=0.) {
    pt = mymuon.q[i]/
         ( mymuon.q[i]/mytrack.bcpt[mymuon.trind[i]]
          -.00042-.00116*sin(mytrack.bcphi0[mymuon.trind[i]]+0.3) );
  } else {
    pt=mymuon.pt[i];
  }
 
  bool highpt = ( pt >  20.0);

  bool iso = mymuon.EtCone4[i]/pt < 0.1;
 
  Float_t d0cor = mymuon.d0[i]; 
  bool goodTrack = GoodTrackQ(mytrack, mymuon.trind[i]);

  if(mytrack.svxhits[mymuon.trind[i]]>0) {
    goodTrack = goodTrack && (fabs(d0cor)<0.02);
  } else {
    goodTrack = goodTrack && (fabs(d0cor)<0.2);
  }


  bool CMUP = 
    ( (mymuon.det[i]&DET_CMU)&&(mymuon.det[i]&DET_CMP) &&
      (fabs(mymuon.CMUDelX[i]) < 3.0) && // units are cm
      (fabs(mymuon.CMPDelX[i]) < 5.0) ); // units are cm

  bool CMU_or_CMNP = 
    ( (mymuon.det[i]&DET_CMU)&&(!(mymuon.det[i]&DET_CMP)) &&
      (fabs(mymuon.CMUDelX[i]) < 3.0) );

  bool CMP = 
    ( (mymuon.det[i]&DET_CMP)&&(!(mymuon.det[i]&DET_CMU)) &&   
      (fabs(mymuon.CMPDelX[i]) < 5.0) );

  bool CMX = 
    ( (mymuon.det[i]&DET_CMX) &&    
     (fabs(mymuon.CMXDelX[i]) < 6.0) );

  // top analysis
  //  bool good_det = (runno >=150799) ?
  //    (CMUP || CMX ) : CMUP;

  //Z cross section analysis

  bool good_det = CMUP;

  bool answer1 =  ( good_det
	       && energyFlow 
	       && highpt && iso
  	       && goodTrack );

  return(answer1);
}





//_____________________________________________________________________________
bool LooseMuonQ(mymuon_t& mymuon, int i, mytrack_t& mytrack, int runno)
{
  //CMIO muons


  bool ans = false;

  bool energyFlow = (mymuon.e[i]>100.) ?
       	 ( (mymuon.em[i] < 2.0 + 0.0115*(mymuon.e[i]-100)) 
       && (mymuon.had[i] < 6.0 + 0.0280*(mymuon.e[i]-100)) ):
       	  (mymuon.em[i] < 2.0) && (mymuon.had[i] < 6.0);

  Float_t pt;

  if(mytrack.bcpt[mymuon.trind[i]]!=0.) {
    pt = mymuon.q[i]/
         ( mymuon.q[i]/mytrack.bcpt[mymuon.trind[i]]
          -.00042-.00116*sin(mytrack.bcphi0[mymuon.trind[i]]+0.3) );
  } else {
    pt=mymuon.pt[i];
  }
 
  bool highpt = ( pt >  20.0);

  bool iso = mymuon.EtCone4[i]/pt < 0.1;
 
  Float_t d0cor = mymuon.d0[i]; 
  bool goodTrack = GoodTrackQ(mytrack, mymuon.trind[i]);

  if(mytrack.svxhits[mymuon.trind[i]]>0) {
    goodTrack = goodTrack && (fabs(d0cor)<0.02);
  } else {
    goodTrack = goodTrack && (fabs(d0cor)<0.2);
  }

  bool answer2 =  ( energyFlow 
	        && highpt && iso
		&& goodTrack );
  return(answer2);
}