#define Sample_cxx
#include "Sample.h"

/*************************************************
Default constructor
*************************************************/
Sample::Sample()
{
	SetType("");
	SetName("");
	SetXSection(1);
	SetKFactor(1);
	SetLuminosity(1);
	SetMETMin(0);
	SetMETMax(1000);
	SetHTMin(0);
	SetHTMax(1000);
	SetMETBins(5);
	SetHTBins(5);

	scale_factor = k_factor*luminosity*sigma;		//Still must be divided by n
	metbinsize = (metmax - metmin)/metbins;
	htbinsize = (htmax - htmin)/htbins;

	table[0].ResizeTo(1, metbins);						//Initialize tables -- special case for monojet table
	table[1].ResizeTo(htbins, metbins);
	table[2].ResizeTo(htbins, metbins);
	table[3].ResizeTo(htbins, metbins);

	monojet_hardLeptonsFailures = 0;		//Cut failure counters
	monojet_highJetEtaFailures = 0;
	monojet_weakJetFailures = 0;
	monojet_strongJetFailures = 0;
	monojet_metFailures = 0;
	monojet_metAndJetsAreCloseFailures = 0;

	dijet_hardLeptonsFailures = 0;
	dijet_highJetEtaFailures = 0;
	dijet_weakJetFailures = 0;
	dijet_strongJetFailures = 0;
	dijet_metFailures = 0;
	dijet_metAndJetsAreCloseFailures = 0;

	trijet_hardLeptonsFailures = 0;
	trijet_highJetEtaFailures = 0;
	trijet_weakJetFailures = 0;
	trijet_strongJetFailures = 0;
	trijet_metFailures = 0;
	trijet_metAndJetsAreCloseFailures = 0;

	tetrajet_hardLeptonsFailures = 0;
	tetrajet_highJetEtaFailures = 0;
	tetrajet_weakJetFailures = 0;
	tetrajet_sphericityFailures = 0;
	tetrajet_metFailures1 = 0;
	tetrajet_metFailures2 = 0;
	tetrajet_metAndJetsAreCloseFailures = 0;
}

Sample::Sample(TString type0, TString name0, Double_t sigma0, Double_t k_factor0, Double_t luminosity0, Double_t metmin0, Double_t metmax0, Double_t htmin0, Double_t htmax0, Int_t metbins0, Int_t htbins0)
{
	SetType(type0);
	SetName(name0);
	SetXSection(sigma0);
	SetKFactor(k_factor0);
	SetLuminosity(luminosity0);
	SetMETMin(metmin0);
	SetMETMax(metmax0);
	SetHTMin(htmin0);
	SetHTMax(htmax0);
	SetMETBins(metbins0);
	SetHTBins(htbins0);

	scale_factor = k_factor*luminosity*sigma;		//Still must be divided by n
	metbinsize = (metmax - metmin)/metbins;
	htbinsize = (htmax - htmin)/htbins;

	table[0].ResizeTo(1, metbins);						//Initialize tables -- special case for monojet table
	table[1].ResizeTo(htbins, metbins);
	table[2].ResizeTo(htbins, metbins);
	table[3].ResizeTo(htbins, metbins);

	monojet_hardLeptonsFailures = 0;		//Cut failure counters
	monojet_highJetEtaFailures = 0;
	monojet_weakJetFailures = 0;
	monojet_strongJetFailures = 0;
	monojet_metFailures = 0;
	monojet_metAndJetsAreCloseFailures = 0;

	dijet_hardLeptonsFailures = 0;
	dijet_highJetEtaFailures = 0;
	dijet_weakJetFailures = 0;
	dijet_strongJetFailures = 0;
	dijet_metFailures = 0;
	dijet_metAndJetsAreCloseFailures = 0;

	trijet_hardLeptonsFailures = 0;
	trijet_highJetEtaFailures = 0;
	trijet_weakJetFailures = 0;
	trijet_strongJetFailures = 0;
	trijet_metFailures = 0;
	trijet_metAndJetsAreCloseFailures = 0;

	tetrajet_hardLeptonsFailures = 0;
	tetrajet_highJetEtaFailures = 0;
	tetrajet_weakJetFailures = 0;
	tetrajet_sphericityFailures = 0;
	tetrajet_metFailures1 = 0;
	tetrajet_metFailures2 = 0;
	tetrajet_metAndJetsAreCloseFailures = 0;
}

/*************************************************
Default destructor
*************************************************/
Sample::~Sample()
{

}

void Sample::SetType(TString s)
{
	type = s;
}

TString Sample::GetType()
{
	return type;
}

void Sample::SetName(TString s)
{
	name = s;
}

TString Sample::GetName()
{
	return name;
}

void Sample::SetXSection(Double_t d)
{
	sigma = d;
}

Double_t Sample::GetXSection()
{
	return sigma;
}

void Sample::SetKFactor(Double_t d)
{
	k_factor = d;
}

Double_t Sample::GetKFactor()
{
	return k_factor;
}

void Sample::SetLuminosity(Double_t d)
{
	luminosity = d;
}

Double_t Sample::GetLuminosity()
{
	return luminosity;
}

Double_t Sample::GetNumOfEvents()
{
	return n;
}

Double_t Sample::GetScaleFactor()
{
	return scale_factor;
}

void Sample::SetMETMin(Double_t d)
{
	metmin = d;
}

Double_t Sample::GetMETMin()
{
	return metmin;
}

void Sample::SetMETMax(Double_t d)
{
	metmax = d;
}

Double_t Sample::GetMETMax()
{
	return metmax;
}

void Sample::SetHTMin(Double_t d)
{
	htmin = d;
}

Double_t Sample::GetHTMin()
{
	return htmin;
}

void Sample::SetHTMax(Double_t d)
{
	htmax = d;
}

Double_t Sample::GetHTMax()
{
	return htmax;
}

void Sample::SetMETBins(Int_t i)
{
	metbins = i;
}

Int_t Sample::GetMETBins()
{
	return metbins;
}

void Sample::SetHTBins(Int_t i)
{
	htbins = i;
}

Int_t Sample::GetHTBins()
{
	return htbins;
}

Double_t Sample::GetMETBinSize()
{
	return metbinsize;
}

Double_t Sample::GetHTBinSize()
{
	return htbinsize;
}

void Sample::Process(TString s)
{
	FILE *oFile;
	Int_t i, j;
	PEvent* E;
	Int_t Tags[5] = {0,0,0,0,0};

	input_file_name = s;
	s.Remove(s.Length()-5, 5);		//Cut of ".lhco"
	output_file_name = s + "_cut";

	cout<<"Input file is "<<input_file_name<<", output file is "<<output_file_name<<endl;
	
	cout<<"Opening input file....."<<endl;
	iFile = fopen(input_file_name.Data(), "r");		//Open input file for reading
	TString summary_file_name = output_file_name+".txt";

	cout<<"Making output summary file....."<<endl;
	oFile = fopen(summary_file_name.Data(), "w");		//Open/make output summary file for writing

	TString output_file_names[4] = {"monojet.", "dijet.", "trijet.", "tetrajet."};
	FILE* samplePGSFile[4];
	for (Int_t i = 0; i < 4; i++)
	{
		output_file_names[i] += (output_file_name+".lhco");

		cout<<"Making output data file #"<<i<<"....."<<endl;
		samplePGSFile[i] = fopen(output_file_names[i].Data(), "w");	//Open/make output data files for writing
	}
	
	Char_t c = ' ';
	while (c != '\n')
	{
		fscanf(iFile, "%c", &c);
		cout<<c;
	}
	
	fscanf(iFile, "%d", &i);

	/*
	for (Int_t i = 0; i < 4; i++)
	{
		TString str = name;
		str.Append("h_mpt");
		str.Append(48+i);
		h_mpt[i] = new TH1D(str.Data(), "Missing p_T", 8, 100, 500);

		str = name;
		str.Append("h_meff");
		str.Append(48+i);
		h_meff[i] = new TH1D(str.Data(), "Effective mass", 8, 400, 2000);

		str = name;
		str.Append("h_pt1");
		str.Append(48+i);
		h_pt1[i] = new TH1D(str.Data(), "Hardest jet p_T", 10, 100, 250);

		str = name;
		str.Append("h_pt2");
		str.Append(48+i);
		h_pt2[i] = new TH1D(str.Data(), "Second hardest jet p_T", 10, 50, 250);
	}
	*/

	Int_t iEvent = 0;
	do				//Loop over all events in LHCO file
	{
		iEvent++;
		E = CreateEvent();
		j = 0;
		//fprintf(stdout,"sample is (%d,%d,%d,%d)\n",E->sample[0],E->sample[1],E->sample[2],E->sample[3]);
		//if (iEvent%10000 == 0)
			cout<<"Event #"<<iEvent<<"....."<<endl;
		for (i = 0; i < 4; i++)
		{
			if (E->sample[i] == 1)
			{
				OutputPGSEvent(samplePGSFile[i], E);
				Tags[i+1] = Tags[i+1] + 1;

				if (!OutOfBounds(E))
				{
					if (i == 0)											//MET vs. H_T tables
						table[i][0][int((E->met-metmin)/metbinsize)]++;
					else
						table[i][int((E->h_t-htmin)/htbinsize)][int((E->met-metmin)/metbinsize)]++;
				}

				//h_mpt[i]->Fill(E->met);				//Missing p_T histogram
				//h_meff[i]->Fill(E->met + E->h_t);	//Effective mass histogram
				//h_pt1[i]->Fill(E->head->pt);		//Hardest jet p_T histogram
				//if (E->head->next->pType == 4)		//If there are at least 2 jets
					//h_pt2[i]->Fill(E->head->next->pt);	//Also fill second hardest jet p_T histogram
			}
			j += E->sample[i];
		}

		if (!j)
			Tags[0]++;
		
		DeleteEvent(E);
	} while(!feof(iFile));	//end of do loop

	n = iEvent;

	scale_factor = luminosity*k_factor*sigma/n;

	for (Int_t i = 0; i < 4; i++)				//Rescale bins
	{
		table[i] = scale_factor*table[i];
		//(*h_mpt[i]) = scale_factor*(*h_mpt[i]);
		//(*h_meff[i]) = scale_factor*(*h_meff[i]);
		//(*h_pt1[i]) = scale_factor*(*h_pt1[i]);
		//(*h_pt2[i]) = scale_factor*(*h_pt2[i]);
	}

	fclose(iFile);
	for (Int_t i = 0; i < 4; i++)
	{
		fclose(samplePGSFile[i]);
		//delete outFileR[i];
	}
	
	fprintf(oFile, "End of reading the file.  Total number of events = %d \n", iEvent);
	fprintf(oFile, "Number of events thrown away = %d \n", Tags[0]);

	fprintf(oFile, "Monojet passes = %d \n", Tags[1]);
	fprintf(oFile, "Failures for lepton cut: %d \n", monojet_hardLeptonsFailures);
	fprintf(oFile, "Failures for eta cut: %d \n", monojet_highJetEtaFailures);
	fprintf(oFile, "Failures for weak jet cut: %d \n", monojet_weakJetFailures);
	fprintf(oFile, "Failures for strong jet cut: %d \n", monojet_strongJetFailures);
	fprintf(oFile, "Failures for MET<=100 cut: %d \n", monojet_metFailures);
	fprintf(oFile, "Failures for close MET and jets cut: %d \n\n", monojet_metAndJetsAreCloseFailures);

	fprintf(oFile, "Dijet passes = %d \n", Tags[2]);
	fprintf(oFile, "Failures for lepton cut: %d \n", dijet_hardLeptonsFailures);
	fprintf(oFile, "Failures for eta cut: %d \n", dijet_highJetEtaFailures);
	fprintf(oFile, "Failures for weak jet cut: %d \n", dijet_weakJetFailures);
	fprintf(oFile, "Failures for strong jet cut: %d \n", dijet_strongJetFailures);
	fprintf(oFile, "Failures for MET<=100 cut: %d \n", dijet_metFailures);
	fprintf(oFile, "Failures for close MET and jets cut: %d \n\n", dijet_metAndJetsAreCloseFailures);

	fprintf(oFile, "Trijet passes = %d \n", Tags[3]);
	fprintf(oFile, "Failures for lepton cut: %d \n", trijet_hardLeptonsFailures);
	fprintf(oFile, "Failures for eta cut: %d \n", trijet_highJetEtaFailures);
	fprintf(oFile, "Failures for weak jet cut: %d \n", trijet_weakJetFailures);
	fprintf(oFile, "Failures for strong jet cut: %d \n", trijet_strongJetFailures);
	fprintf(oFile, "Failures for MET<=100 cut: %d \n", trijet_metFailures);
	fprintf(oFile, "Failures for close MET and jets cut: %d \n\n", trijet_metAndJetsAreCloseFailures);

	fprintf(oFile, "Tetrajet passes = %d \n", Tags[4]);
	fprintf(oFile, "Failures for lepton cut: %d \n", tetrajet_hardLeptonsFailures);
	fprintf(oFile, "Failures for eta cut: %d \n", tetrajet_highJetEtaFailures);
	fprintf(oFile, "Failures for weak jet cut: %d \n", tetrajet_weakJetFailures);
	fprintf(oFile, "Failures for sphericity cut: %d \n", tetrajet_sphericityFailures);
	fprintf(oFile, "Failures for MET<=100 cut: %d \n", tetrajet_metFailures1);
	fprintf(oFile, "Failures for MET<=0.2*M_eff cut: %d \n", tetrajet_metFailures2);
	fprintf(oFile, "Failures for close MET and jets cut: %d \n\n", tetrajet_metAndJetsAreCloseFailures);

	fclose(oFile);
}

/*************************************************
CreateEvent takes the next PGS from iFile, reads it in, and then tags the event as a monojet, dijet, trijet, or tetrajet
*************************************************/
PEvent* Sample::CreateEvent()
{
	PObject* C;
	Int_t Tag;

	//fprintf(stdout,"Start reading the event\n");
	PEvent* E = ReadEvent();
	//fprintf(stdout, "Done reading the event\n");
	//PrintEvent(stdout,E->head);
	
	E->head=SortJets(E->head);
	//PrintEvent(stdout,E->head);

	MergeMuons(E);
	//PrintEvent(stdout,E->head);
	
	ConvertTaus(E);
	//PrintEvent(stdout,E->head);

	FiducialJets(E);
	//PrintEvent(stdout,E->head);

	E->head = SortJets(E->head);
	//PrintEvent(stdout,E->head);
	//Check everything working well up to this point

	E->njet = 0;
	E->nbjet = 0;
	E->nelectronp = 0;
	E->nelectronm = 0;
	E->nmuonp = 0;
	E->nmuonm = 0;
	E->ntaup = 0;
	E->ntaup = 0;
	E->nphoton = 0;
	E->met = 0;
	E->h_t = 0;

	for (C = E->head; C != NULL; C = C->next)
	{
		if (C->pType == 0)
			E->nphoton++;
		if (C->pType == 1)
		{
			if (C->ntrk < 0)
				E->nelectronm++;
			if (C->ntrk > 0)
				E->nelectronp++;
		}
		if (C->pType == 2)
		{
			if (C->ntrk < 0)
				E->nmuonm++;
			if (C->ntrk > 0)
				E->nmuonp++;
		}
		if (C->pType == 3)
		{
			if (C->ntrk < 0)
				E->ntaum++;
			if (C->ntrk > 0)
				E->ntaup++;
		}
		if (C->pType == 4)
		{
			E->njet++;
			E->h_t+=C->pt;
		}
		if (C->btag > 0)
			E->nbjet++;
		if (C->pType == 6)
			E->met+=C->pt;
		
		//if (C->pType != 6)
			//E->h_t+=C->pt;
	}
	
	//At this point, we know some accounting for this event
	Tag = JetTag(E);	//Run the cuts on the event
	if (Tag != 0 && Tag != 1 && Tag != 2 && Tag != 4 && Tag != 8)
	{
		fprintf(stdout, "Tag: %d\n", Tag);
		PrintEvent(stdout, E->head);
	}
    
	return(E);
}

/*************************************************
ReadEvent reads in a PGS event, builds the linked list of objects in the event, and returns the head of the linked list 
*************************************************/
PEvent* Sample::ReadEvent()
{
	PObject *head = NULL, *Current = NULL;
	PEvent* evt = new PEvent;
	
	Int_t pID, pType;
	Double_t eta, phi, pt, jmas, ntrk, btag, hem, dum;
	if (iFile == NULL)
		return NULL;

	fscanf(iFile, "%d", &pID);
	//fprintf(stdout,"PID1: %d\t",pID);
	evt->trigger1 = pID;
	fscanf(iFile, "%d", &pID);
	//fprintf(stdout,"PID2: %d\n",pID);
	evt->trigger2 = pID;
	fscanf(iFile, "%d", &pID);
	//fprintf(stdout,"PID3: %d\n",pID);

	do
	{
		if (Current == NULL)
		{
			head = new PObject;
			Current = head;
			Current->prev = NULL;
			Current->next = NULL;
		}
		else
		{
			Current->next = new PObject;
			Current->next->prev = Current;
			Current = Current->next;
			Current->next = NULL;
		}
		
		Current->pID = pID;
		fscanf(iFile, "%d", &pType);
		Current->pType = pType;
		fscanf(iFile, "%lf", &eta);
		Current->eta = eta;
		fscanf(iFile, "%lf", &phi);
		Current->phi = phi;
		fscanf(iFile, "%lf", &pt);
		Current->pt = pt;
		fscanf(iFile, "%lf", &jmas);
		Current->jmas = jmas;
		fscanf(iFile, "%lf", &ntrk);
		Current->ntrk = ntrk;
		fscanf(iFile, "%lf", &btag);
		Current->btag = btag;
		fscanf(iFile, "%lf", &hem);
		Current->hem = hem;
		fscanf(iFile, "%lf", &dum);
		fscanf(iFile, "%lf", &dum);
		fscanf(iFile, "%d", &pID); //the next number, if it's zero, it's signal the next event
		//fprintf(stdout,"pID: %d\n",pID);
	} while (pID != 0 && !feof(iFile));

	evt->head = head;
	return evt;
}

/*************************************************
MergeMuons merges any muons into jets
*************************************************/
Int_t Sample::MergeMuons(PEvent *E)
{
  PEvent *muList = new PEvent;
  muList->head= NULL;

  Int_t nMergeable = LabelMuons(E,muList); //get the number of mergeable muons
  
  //fprintf(stdout, "Mergemuon\n");
  //PrintEvent(stdout,muList->head);
  if(nMergeable >=1) {
    muList->head = PtSort(muList->head, 2,-1);  
    //fprintf(stdout, "Mergeable muon list 2 %d\n",nMergeable);
    //PrintEvent(stdout,muList->head);
  }  
  PObject *Ce,*Cm, *nearestObject,*newHead; 
  Double_t deltaR = 0;

  //Now delete all the mergeable Muon from the list
  newHead = E->head;  
  for (Ce=E->head; Ce!=NULL;) {
    //fprintf(stdout,"C=%x,\tC->prev=%x,\tC->next=%x\n",Ce,Ce->prev,Ce->next);
    Cm = Ce->next;
    if(Ce->pType == 2 && Ce->mergeable){
      //fprintf(stdout,"delete this object %x\n",Ce);
      if(Ce->prev == NULL){  //deleting the first object
	newHead = Ce->next;
      }else{
	Ce->prev->next = Ce->next;
      }
      if(Ce->next!=NULL)
	Ce->next->prev = Ce->prev;

      delete Ce;
    }
    Ce = Cm;
  }
  E->head = newHead;
  //fprintf(stdout, "Leftover list\n");
  //PrintEvent(stdout,E->head);  
  
  for (Cm=muList->head; Cm!=NULL; Cm=Cm->next)  // go thru the list the and process each mergeable muon
  {
    //Now for each of the mergeable muon on the list, to thru each of the object in the main event, compute the object that has smallest DeltaR to the it 
    //(least negative if we take -DeltaR and put -99 if the object is other than a jet) keep track of its index in the link list and 
    //(1) if the closest object is a jet => replace the object with the new jet that has the 4 mom as sum of old jet + muon and delete the muon
    //(2) if the closest object is NOT a jet => just delete the object
    Double_t minDeltaR = -99;
    for (Ce=E->head; Ce!= NULL; Ce=Ce->next) {
      //Double_t dR(Double_t eta1,Double_t phi1, Double_t eta2,Double_t phi2);
      
      if (Ce->pType==4) {
	deltaR = -dR(Ce->eta,Ce->phi,Cm->eta,Cm->phi);
      }

      if(deltaR > minDeltaR){
	minDeltaR = deltaR;
        nearestObject = Ce;
      }
    } 
    //fprintf(stdout, "the nearest object of Cm = %3.2f,\t is %x\n",Cm->pt,nearestObject);
    // now check if the nearestObject is a Jet, then merge and delete the muon
    if(nearestObject->pType==4){
      FourVect *temp= add2FourVect(convertTo4Vect(Cm),convertTo4Vect(nearestObject)); 
      convertFrom4Vect(nearestObject,temp);
    }
    //PrintEvent(stdout,E->head);   
  }
  //fprintf(stdout, "Resulted list\n");
  //PrintEvent(stdout,E->head);  
  delete muList;
  return(1);
}


FourVect* Sample::convertTo4Vect(PObject *obj){

  FourVect *res = new FourVect;
  res->e = sqrt(pow(obj->pt*cosh(obj->eta),2)+pow(obj->jmas,2));
  res->pz = obj->pt*sinh(obj->eta);
  res->px = obj->pt*cos(obj->phi);
  res->py = obj->pt*sin(obj->phi);
  return res;
}

void Sample::convertFrom4Vect(PObject *obj, FourVect *vect4)
{
  obj->pt =  sqrt(pow(vect4->px,2) + pow(vect4->py,2));
  obj->eta = TMath::ASinH( (vect4->pz)/(obj->pt) );
  obj->phi = atan2(vect4->py,vect4->px);
  if(obj->phi < 0) obj->phi += 2 *Pi;
  obj->jmas = sqrt(pow(vect4->e,2) - pow(vect4->px,2)- pow(vect4->py,2) - pow(vect4->pz,2));
}
  
//
// Int_t LabelMuons(PEvent * E)
// Labels muons as mergeable or Isolated or Neither, and return a Events (a list of mergeable muon)
//
FourVect* Sample::add2FourVect(FourVect *v1,FourVect *v2){

  FourVect *vres = new FourVect;
  vres->e = v1->e + v2->e;
  vres->px = v1->px + v2->px;
  vres->py = v1->py + v2->py;
  vres->pz = v1->pz + v2->pz;
  return vres;
}


PObject* Sample::CopyObject(PObject *Obj)
{
  PObject *nObj = new PObject; //create the new object
  nObj->prev = NULL;
  nObj->next = NULL;
  nObj->pID = Obj->pID;
  nObj->pType= Obj->pType;
  nObj->phi=Obj->phi;
  nObj->eta=Obj->eta;
  nObj->pt= Obj->pt;
  nObj->jmas= Obj->jmas;
  nObj->ntrk=Obj->ntrk;
  nObj->btag= Obj->btag;
  nObj->hem= Obj->hem;
  nObj->mergeable=Obj->mergeable;
  nObj->isolated =Obj->isolated;  
  return nObj;
}


Int_t Sample::LabelMuons(PEvent * E, PEvent *MergeableMuon)
{
  PObject * C,*mu;
  Double_t PTIso,ETRatio;
  Int_t M=0;
  //fprintf(stdout, "LabelMuons\t");

  for(C=E->head;C!=NULL;C=C->next)
  {
    if(C->pType==2)
    {
      PTIso=floor(C->hem); 
      ETRatio=C->hem-PTIso;
      if(ETRatio>= 0.1*9/8 || PTIso >= 5.0)
      {
       C->mergeable=kTRUE;
       //now make that list
       if(M==0) {
	 MergeableMuon->head = CopyObject(C);
         mu =MergeableMuon->head;
	 mu->prev = NULL;
	 mu->next = NULL;
         //fprintf(stdout,"First mergeable muon head=%x\n",mu);
       } else {
         //fprintf(stdout,"First mergeable muon current =%x of M=%d\n",mu,M);
	 mu->next= CopyObject(C);
         mu->next->prev = mu;
         mu = mu->next;
         mu->next = NULL;
       }
       M++;
      }
      else
	{C->mergeable=kFALSE; /*fprintf(stdout,"immergeable!!!\n");*/ }
      if(PTIso<5.0&& ETRatio<0.1*9/8)
	{C->isolated=kTRUE; /*fprintf(stdout,"isolatedMuon!!!\n");*/ }
      else
       C->isolated=kFALSE;
    }
  }
  return(M);
}


// the event appears like this:
// (1) Hardest Jet
// (2) Second Hardest jet, etc...
// (3) all other objects
// (4) MET

Int_t Sample::Construct4JLeptonVeto(PEvent *E, PObject* METJetList[]){
  // return 0 => it's vetoed by the lepton, number of jet and MET cut, return 1 => passes the these vetoes 
  // the [0]=>MET, [1]=>Hardest jet, [4]=> 4th hardest jet
  PObject *C;
  Int_t jetCount = 0;

  if(E->njet < 1)  return 0;  // get out right away, no jet at all
  if(E->met < 40) return 0; // get out right away, too soft a MET

  for (C=E->head; C!=NULL; C=C->next)
  {
    if(C->pType == 4 && jetCount <4){
      jetCount ++;
      METJetList[jetCount] = C;
    }
    if(C->pType==6)   METJetList[0]= C;

    //now do the isolated muon and electron veto
    if(C->pType == 2) {// if it's muon
      if(C->isolated && C->pt >= 10) return 0; // break it out
    }
    if(C->pType == 1) 
      if(C->pt >=10)  return 0; 

  }
  return (1); // OK to proceed
}

//
// Int_t MonoJetTag(PObject * head)
// Tags monojets, returns 0 for an event that doesn't pass the cuts
// and 1 for events that do.
//

Int_t Sample::MonoJetTag(PEvent *E)
{
	/*
  PObject *metJetList[5];  
  if(Construct4JLeptonVeto(E,metJetList)== 0) 
    return 0;// get out now

  // Now go thru 1,2,3,4 jets and MET to make the cuts

  if( (metJetList[1]->pt < 150) || fabs(metJetList[1]->eta) > 0.8) return 0;  //hardest jet has pt=>150 and eta<=0.8
  
  if( dPhi(metJetList[0]->phi,metJetList[1]->phi) < 90.0/180*Pi ) return 0;   //deltaPhi (Jet1,MET) >= 90*PI/180

  if(E->njet > 1){
    if( metJetList[2]->pt >=35 && fabs(metJetList[2]->eta) <= 0.8 ) return 0;  
    if( dPhi(metJetList[0]->phi,metJetList[2]->phi) < 50.0/180.0*Pi ) return 0;
    
    if( dPhi(metJetList[1]->phi,metJetList[2]->phi) > 165.0/180.0*Pi ) return 0;
  }
  if(E->njet >2){
    if(metJetList[3]->pt >= 35.0) return 0;
  }
  if(E->njet >3){
    if(metJetList[4]->pt >= 20.0) return 0;
  }

  return 1;
  */

	Int_t good = 1;

	if (HardLeptons(E))		//Make sure there are no hard leptons
	{
		good = 0;
		monojet_hardLeptonsFailures++;
	}
	if (HighJetEta(E,1))		//Make sure the hardest jet doesn't have too high eta value
	{
		good = 0;
		monojet_highJetEtaFailures++;
	}
	if (WeakJet(E,1))			//Make sure the hardest jet is sufficiently hard
	{
		good = 0;
		monojet_weakJetFailures++;
	}
	if (StrongJet(E,1))	//Make sure the other jets are sufficiently soft
	{
		good = 0;
		monojet_strongJetFailures++;
	}
	if (E->met <= 100)		//Make sure the MET is sufficiently high
	{
		good = 0;
		monojet_metFailures++;
	}
	if (METAndJetsAreClose(E,1))	//Make sure the MET and the hardest jet aren't too close
	{
		good = 0;
		monojet_metAndJetsAreCloseFailures++;
	}

	return good;
}


Int_t Sample::DiJetTag(PEvent *E)
{
	/*
  PObject *metJetList[5];
  PObject *C;  
  if(Construct4JLeptonVeto(E,metJetList)== 0) return 0;// get out now

  if(E->njet < 2) return 0;
  // Now go thru 1,2,3,4 jets and MET to make the cuts

  if( (metJetList[1]->pt < 35) || fabs(metJetList[1]->eta) > 0.8) return 0;  //hardest jet has pt>150 and eta<0.8
  if( (metJetList[2]->pt < 35) || fabs(metJetList[2]->eta) > 0.8) return 0;  //hardest jet has pt>150 and eta<0.8
  
  if( dPhi(metJetList[0]->phi,metJetList[1]->phi) < 90.0/180.0*Pi ) return 0;   //deltaPhi (Jet1,MET) >= 90*PI/180
  if( dPhi(metJetList[0]->phi,metJetList[2]->phi) < 50.0/180.0*Pi ) return 0;
  if( dPhi(metJetList[1]->phi,metJetList[2]->phi) > 165.0/180.0*Pi ) return 0;

  if(E->njet > 2 ){
    if(metJetList[3]->pt >= 35.0) return 0;
  }
  if(E->njet > 3){
    if(metJetList[4]->pt >= 20.0) return 0;
  }
  //Additional Phi cut between the Jets and MET direction, it must be greater than 40 degree
  for (C=E->head; C!=NULL; C=C->next)
    if(C->pType == 4 && dPhi(C->phi,metJetList[0]->phi) <= 40.0*Pi/180.0) 
      return 0;
  

  return 1;
  */

	if(E->njet < 2)
		return 0;

	Int_t good = 1;

	if (HardLeptons(E))		//Make sure there are no hard leptons
	{
		good = 0;
		dijet_hardLeptonsFailures++;
	}
	if (HighJetEta(E,2))		//Make sure the hardest jet doesn't have too high eta value
	{
		good = 0;
		dijet_highJetEtaFailures++;
	}
	if (WeakJet(E,2))			//Make sure the hardest jet is sufficiently hard
	{
		good = 0;
		dijet_weakJetFailures++;
	}
	if (StrongJet(E,2))	//Make sure the other jets are sufficiently soft
	{
		good = 0;
		dijet_strongJetFailures++;
	}
	if (E->met <= 100)		//Make sure the MET is sufficiently high
	{
		good = 0;
		dijet_metFailures++;
	}
	if (METAndJetsAreClose(E,2))	//Make sure the MET and the hardest jet aren't too close
	{
		good = 0;
		dijet_metAndJetsAreCloseFailures++;
	}

	return good;
}

Int_t Sample::TriJetTag(PEvent *E)
{
	/*
  PObject *metJetList[5];
  PObject *C;  
  if(Construct4JLeptonVeto(E,metJetList)== 0) return 0;// get out now

  if(E->njet < 3) return 0;
  // Now go thru 1,2,3,4 jets and MET to make the cuts

  if( (metJetList[1]->pt < 35) || fabs(metJetList[1]->eta) > 0.8) return 0;  //hardest jet has pt>150 and eta<0.8
  if( (metJetList[2]->pt < 35) || fabs(metJetList[2]->eta) > 0.8) return 0;  //hardest jet has pt>150 and eta<0.8
  if( metJetList[3]->pt < 35.0) return 0;
  
  if( dPhi(metJetList[0]->phi,metJetList[1]->phi) < 90.0/180*Pi ) return 0;   //deltaPhi (Jet1,MET) >= 90*PI/180
  if( dPhi(metJetList[0]->phi,metJetList[2]->phi) < 50.0/180.0*Pi ) return 0;
  
  if( dPhi(metJetList[1]->phi,metJetList[2]->phi) > 165.0/180.0*Pi ) return 0;

  if(E->njet > 3){
    if(metJetList[4]->pt >= 20.0) return 0;
    }

  return 1;
  */

	if(E->njet < 3)
		return 0;

	Int_t good = 1;

	if (HardLeptons(E))		//Make sure there are no hard leptons
	{
		good = 0;
		trijet_hardLeptonsFailures++;
	}
	if (HighJetEta(E,3))		//Make sure the hardest jet doesn't have too high eta value
	{
		good = 0;
		trijet_highJetEtaFailures++;
	}
	if (WeakJet(E,3))			//Make sure the hardest jet is sufficiently hard
	{
		good = 0;
		trijet_weakJetFailures++;
	}
	if (StrongJet(E,3))	//Make sure the other jets are sufficiently soft
	{
		good = 0;
		trijet_strongJetFailures++;
	}
	if (E->met <= 100)		//Make sure the MET is sufficiently high
	{
		good = 0;
		trijet_metFailures++;
	}
	if (METAndJetsAreClose(E,3))	//Make sure the MET and the hardest jet aren't too close
	{
		good = 0;
		trijet_metAndJetsAreCloseFailures++;
	}

	return good;
}

//LHC Modification: New cuts added, old ones commented out
//Note that the current cut on jets is a somewhat simplified version of the ATLAS zero-lepton cuts; eventually, we need
//to organize jets by eta first (central vs. non-central jets) and then make sure the central jets (|eta|<2.5) have 
//sufficiently high energies
Int_t Sample::TetraJetTag(PEvent *E)
{
  //PObject *metJetList[5];
  //PObject *C;  
  Int_t good = 1;
  
  //if(Construct4JLeptonVeto(E,metJetList)== 0) return 0;// get out now

  if(E->njet < 4) return 0;
  // Now go thru 1,2,3,4 jets and MET to make the cuts

  /*
  if( (metJetList[1]->pt < 35) || fabs(metJetList[1]->eta) > 0.8) return 0;  //hardest jet has pt>150 and eta<0.8
  if( (metJetList[2]->pt < 35) || fabs(metJetList[2]->eta) > 0.8) return 0;  //2nd hardest jet has pt>150 and eta<0.8
  if( metJetList[3]->pt < 35.0) return 0;
  if( metJetList[4]->pt < 20.0) return 0;
  
  if( dPhi(metJetList[0]->phi,metJetList[1]->phi) < 90.0/180*Pi ) return 0;   //deltaPhi (Jet1,MET) >= 90*PI/180
  if( dPhi(metJetList[0]->phi,metJetList[2]->phi) < 50.0/180.0*Pi ) return 0;
  if( dPhi(metJetList[1]->phi,metJetList[2]->phi) > 165.0/180.0*Pi ) return 0;
  */

	if (HardLeptons(E))		//Make sure there are no hard leptons
	{
		good = 0;
		tetrajet_hardLeptonsFailures++;
	}
	if (HighJetEta(E,4))		//Make sure the jets don't have too high eta values
	{
		good = 0;
		tetrajet_highJetEtaFailures++;
	}
	if (WeakJet(E,4))			//Make sure the jets are sufficiently hard
	{
		good = 0;
		tetrajet_weakJetFailures++;
	}
	if (CalcSph_T(E) <= 0.2)	//Make sure the sphericity is good
	{
		good = 0;
		tetrajet_sphericityFailures++;
	}
	if (E->met <= 100)		//Make sure the MET is sufficiently high
	{
		good = 0;
		tetrajet_metFailures1++;
	}
	if (E->met <= 0.2*CalcM_eff(E))
	{
		good = 0;
		tetrajet_metFailures2++;
	}
	if (METAndJetsAreClose(E,4))	//Make sure the MET and the jets aren't too close
	{
		good = 0;
		tetrajet_metAndJetsAreCloseFailures++;
	}

	/*
	if (HardLeptons(E))		//Make sure there are no hard leptons
		return 0;
	if (HighJetEta(E))		//Make sure the jets don't have too high eta values
		return 0;
	if (WeakJet(E))			//Make sure the jets are sufficiently hard
		return 0;
	if (CalcSph_T(E) <= 0.2)	//Make sure the sphericity is good
		return 0;
	if (E->met <= 100)		//Make sure the MET is sufficiently high
		return 0;
	if (E->met <= 0.2*CalcM_eff(E))
		return 0;
	if (METAndJetsAreClose(E))	//Make sure the MET and the jets aren't too close
		return 0;
	*/

	return good;
	//return 1;
}

//
// Int_t JetTag(PEvent *E)
// The global tagging functions
// It returns what the events were tagged with encoded in binary.
//
Int_t Sample::JetTag(PEvent *E)
{
  
  E->sample[1]=DiJetTag(E);
  E->sample[2]=TriJetTag(E);
  E->sample[3]=TetraJetTag(E);
  E->sample[0]=MonoJetTag(E);
  
  return(E->sample[0]+2*E->sample[1]+4*E->sample[2]+8*E->sample[3]);
}

//
// Int_t FiducialJets(PEvent *E)
// Removes jets that have pt<15 GeV or |eta|>2.5
//

Int_t Sample::FiducialJets(PEvent *E)
{
  //my code
  PObject * C, *T,*newHead = E->head;
  Int_t count=0;
  //fprintf(stdout, "FiducialJet\t"); 
   
  for(C=E->head; C!= NULL;)
  {
    //fprintf(stdout,"this object %d,\t%3.2f,\t%3.2f, C=%x,\tC->next=%x\tC->prev=%x\n", C->pType,C->pt,C->eta,C,C->next,C->prev);
    if(C->pType == 4 && (C->pt < PTMin || fabs(C->eta) >EtaMax)) //remove this jet
      { 
        T = C->next;
	if(C== newHead)
	  {
	    newHead = C->next;
            C->next->prev = NULL;
            
          }else
	  {
            C->prev->next = C->next;
            C->next->prev = C->prev;
	  }
	delete C;
        C= T; 
      }else
      {   
	if(C->pType==4) count++;
	  C=C->next;
      }	 
  }
  E->head = newHead;
  //fprintf(stdout,"Number of Fiducial Jets = %d\n",count); 
  return (count); 
}

//
// Int_t ConvertTaus(PEvent* E)
// Converts taus into jets
//

Int_t Sample::ConvertTaus(PEvent* E)
{
  PObject * C;
  Int_t count=0;
  //fprintf(stdout, "ConvertTaus\t");
  for(C=E->head;C!=NULL;C=C->next)
  {
   if(C->pType==3)
   {
     C->pType=4;
     count++;
   }
  } 
  return(count);
}

//
// Double_t dR(Double_t eta1,Double_t phi1, Double_t eta2,Double_t phi2)
// Computes the angular distance between two objects
//

Double_t Sample::dR(Double_t eta1,Double_t phi1, Double_t eta2,Double_t phi2)
{
  Double_t R, phi;
  phi= dPhi(phi1,phi2);
  R = sqrt( (eta1-eta2)*(eta1-eta2) + phi*phi);
  return(R);
}

//
// Double_t dPhi( Double_t phi1, Double_t phi2)
// Computes the azimuthal angular distance between two objects
//

Double_t Sample::dPhi( Double_t phi1, Double_t phi2)
{
Double_t phi;
phi = fabs( phi1-phi2);
if(phi>Pi)
  phi= fabs(2*Pi-phi);
return(phi);
}

//
// Int_t DeleteEvent(PEvent *E)
// Deletes an event
//

Int_t Sample::DeleteEvent(PEvent *E)
{
  PObject *C;
  for(C=E->head;C->next!=NULL;)
  {
    C=C->next;
    delete C->prev;
  }
  delete C;
  delete E;
  return(1);
}

//
// Int_t PrintEvent(FILE* OFile, PObject * Head)
// Prints out a barebones list of objects in an event
//

Int_t Sample::PrintEvent(FILE* OFile, PObject * head)
{
  PObject *C;
  
  for(C=head;C!=NULL;C=C->next)
  {
    //fprintf(OFile,"%x \n",C);  
    fprintf(OFile,"%d\t%3.2f\t%3.2f\t%3.2f\t%3.2f\t%3.2f\t%x\t%x\t%x\n", C->pType, C->pt, C->eta,C->phi, C->jmas,C->hem,C,C->prev, C->next);
  }
  //fprintf(OFile,"get out = %x \n",C);    
  return 1;
}

//
// Int_t OutputPGSEvent(FILE * OFile, PObject *head)
// Outputs PGS style event
//

Int_t Sample::OutputPGSEvent(FILE * OFile, PEvent *evt)
{
  Int_t i=0;
  PObject* C;
  fprintf(OFile,"%d   %d   %d\n", i, evt->trigger1, evt->trigger2);
  for(C=evt->head;C!=NULL;C=C->next)
  {
    i++;
    fprintf(OFile,"%d  %d ",i, C->pType);
    fprintf(OFile,"%3.3f  %3.3f ",C->eta, C->phi);
    fprintf(OFile,"%3.2f  %3.2f ",C->pt, C->jmas);
    fprintf(OFile,"%3.1f  %3.1f ",C->ntrk, C->btag);
    fprintf(OFile,"%3.2f 0.0 0.0\n",C->hem);
  }
  
  return(1);
}


//
// Int_t OutputEvent(FILE * OFile, PEvent *E)
// Prints a brief su;mmary of an event including 
// what sample it was tagged under
//

Int_t Sample::OutputEvent(FILE * OFile, PEvent *E)
{
  fprintf(OFile,"njet=%d\tMET=%3.2f\tHT=%3.2f\tnJ=(%d,%d,%d,%d)\n",E->njet, E->met, E->h_t,E->sample[0], E->sample[1], E->sample[2], E->sample[3]);
  return(1);
}

//
// PObject * SortJets(PObject * head)
// This function sorts the link list so that
// jets occur before any other objects and are pt ordered
//
PObject* Sample::SortJets(PObject * head) 
{  return PtSort(head,4,1);
}
//typeSort 1 -> higher to lower
//typeSort -1 -> lower to higher
PObject* Sample::PtSort(PObject * head, Int_t particleType,Int_t typeSort)
{
  PObject * NewHead=head, *C,*T;
  Double_t pt1,pt2;
  //fprintf(stdout, "SortJets, %x, %x \n", head, head->next);
 
  for(C=head->next;C!=NULL;C=C->next)
  {
    if(C->pType==particleType)
      pt1=C->pt; 
    else{
      if (typeSort ==1)
	pt1=0;
      else pt1 = OVERMAXPT;
    }

    if(C->prev->pType==particleType)
      pt2=C->prev->pt; 
    else {
      if (typeSort ==1)
	pt2=0;
      else pt2 = OVERMAXPT;
    }
    
    // for the case of two jet having the same pt, the one more central would be considered higher pt
    if(pt1 == pt2 && pt1!=0) {
      if(fabs(C->eta) < fabs(C->prev->eta))  // the pt1 is more central
	pt1 = pt1+  0.00001; // make sure pt1 is higher
      else  //pt2 is more central
        pt2 = pt2 + 0.00001;
    }

    if( (typeSort == 1 && pt2<pt1) || (typeSort==-1 && pt2 >pt1) )  
    {  //need to swap??? 
      if(C->prev==NewHead)
      {
        T=C->prev;  
        T->next=C->next; 
        if(T->next!=NULL)
	  T->next->prev=T;
        T->prev=C;  
        C->prev=NULL;
        C->next=T;
        NewHead=C;
      }
      else
      {
        if(C->next!=NULL)
        {
          T = C-> prev;
          T->next = C->next;
          C->next = T;
          C->prev= T->prev;
          T->prev = C;
          C->prev->next=C;
          T->next->prev = T;
          C=C->prev; 	  //going back
        }
        else
        {
          T=C->prev;
          C->prev=T->prev;
          C->next=T;
          C->prev->next=C;
          T->prev=C;
          T->next=NULL;
          C=C->prev;  // going back
        }
      }
      
    }
  } 

  //fprintf(stdout,"return New head = %x, next %x\n", NewHead, NewHead->next);
  return(NewHead);
}

//LHC Modification: Check for close MET and jets
Bool_t Sample::METAndJetsAreClose(PEvent* E, Int_t jets)
{
	PObject* met;
	//First we must find the MET
	for (met = E->head; met != NULL; met = met->next)
	{
		if (met->pType == 6)
			break;
	}

	//if ( (jets == 1) || (jets == 2) )
	//{
		//Now compare MET to the 3 hardest jets, or however many jets there are
		PObject* cur = E->head;
		for (Int_t j = 0; j < 3; j++)
		{
			if (cur == NULL)
				break;
			if (dPhi(cur->phi, met->phi) <= 0.2)
				return kTRUE;

			cur = cur->next;
		}
	//}

	/*
	if ( (jets == 3) || (jets == 4) )
	{
		//Now compare MET to the 3 hardest jets
		PObject* cur = E->head;
		for (Int_t j = 0; j < 3; j++)
		{
			if (dPhi(cur->phi, met->phi) <= 0.2)
				return kTRUE;

			cur = cur->next;
		}
	}
	*/

	return kFALSE;
}

//LHC Modification: Check for high jet eta values on hardest prime jets
Bool_t Sample::HighJetEta(PEvent* E, Int_t jets)
{
	PObject* cur = E->head;

	for (Int_t j = 0; j < jets; j++)	//Check only for hardest prime jets
	{
		if (cur == NULL)
			break;
		if (fabs(cur->eta) >= 2.5)	//Cut if it has high eta
			return kTRUE;
		cur = cur->next;
	}

	return kFALSE;
}

//LHC Modification: Check for weak jets
Bool_t Sample::WeakJet(PEvent* E, Int_t jets)
{
	PObject* cur = E->head;
	if (cur->pt <= 100)		//Make sure hardest jet (which should be first) has at least 100 GeV
		return kTRUE;

	//Now make sure hardest prime jets have at least 50 GeV
	for (Int_t j = 0; j < jets; j++)
	{
		if (cur->pt <= 50)
			return kTRUE;
		
		cur = cur->next;
	}

	return kFALSE;
}

Bool_t Sample::StrongJet(PEvent* E, Int_t jets)
{
	PObject* cur = E->head;

	for (Int_t j = 0; j < jets; j++)	//Skip the "prime" jets
		cur = cur->next;

	if (jets < 4)		//Only check this for monojet, dijet, and trijet cases
	{
		//Make sure other jets (after hardest one) are either sufficiently soft or sufficiently non-central
		for (; cur != NULL; cur = cur->next)
		{
			if (cur->pType == 4)			//If it's a jet
			{
				if ( (cur->pt > 50) || (fabs(cur->eta) < 2.5) )	//Cut if it has high pt or is too central
					return kTRUE;
			}
		}
	}

	return kFALSE;
}

//LHC Modification: Check for hard leptons
Bool_t Sample::HardLeptons(PEvent* E)
{
	for (PObject* cur = E->head; cur != NULL; cur = cur->next)
	{
		if (cur->pType == 1 || cur->pType == 2)	//For electrons and muons
		{
			if (cur->pt > 20)
				return kTRUE;
		}
	}

	return kFALSE;
}

//LHC Modification: Computes effective mass by summing pt of 4 hardest jets, all leptons, and MET
//Note: This only works for the tetrajet case right now!
Double_t Sample::CalcM_eff(PEvent* E)
{
	Double_t m_eff = 0;
	Double_t jetcount = 0;

	for (PObject* cur = E->head; cur != NULL; cur = cur->next)
	{
		if (cur->pType == 4 && jetcount <= 4)	//Sum if object is among 4 hardest jets
		{
			m_eff += cur->pt;
			jetcount++;
		}
		if (cur->pType == 1 || cur->pType == 2 /*|| cur->pType == 3*/)	//Sum if object is a lepton
			m_eff += cur->pt;
	}

	return m_eff + E->met;	//Finally, add MET
}

//LHC Modification: Computes the transverse sphericity
Double_t Sample::CalcSph_T(PEvent* E)
{
	Double_t lambda1, lambda2;
	Double_t a, b, c;
	a = b = c = 0;
	PObject* cur;

	cur = E->head;
	while (kTRUE)
	{
		if (cur == NULL || cur->pType == 6)
			break;
		a += pow((cur->pt)*cos(cur->phi), 2);
		b += pow(cur->pt, 2)*cos(cur->phi)*sin(cur->phi);
		c += pow((cur->pt)*sin(cur->phi), 2);

		cur = cur->next;
	}
	lambda1 = 0.5*( (a + c) - sqrt(pow(a, 2) + pow(c, 2) + 4*pow(b, 2) - 2*a*c) );
	lambda2 = 0.5*( (a + c) + sqrt(pow(a, 2) + pow(c, 2) + 4*pow(b, 2) - 2*a*c) );

	return 2*lambda1/(lambda1+lambda2);
}

//OutOfBounds functions: return -1 for values too low, 0 for in bounds, and 1 for too high
Bool_t Sample::OutOfBounds(PEvent* ev)
{
	if (ev->met > metmax || ev->met < metmin || ev->h_t > htmax || ev->h_t < htmin)
		return kTRUE;

	return kFALSE;
}

//
//  Int_t CountObjects(PObject * head)
//  Counts the number of objects in a linked list
//  This is mostly a debugging function
//

Int_t Sample::CountObjects(PObject * head)
{
  Int_t c=0;
  PObject * C;
  for(C=head;C!=NULL;C=C->next)
   c++;
  return(c); 
}

//
//  Int_t FixList(PObject * head)
//  Fixes the linked list to make sure it's properly linked
//  This is a debugging function that shouldn't regularly be called.
//

Int_t Sample::FixList(PObject * head)
{
  Int_t fix=0;
  PObject *C;
  for(C=head;C!=NULL;C=C->next)
  {
    if(C->next!=NULL)
    {
       if(C->next->prev!=C)
       {
         fix++;
         C->next->prev=C;
       }
    }
  }
  if(fix!=0) 
     fprintf(stdout,"Fixes: %d\n", fix); 
  return(fix);
}

/*
Sample::Sample(Double_t xsection, Double_t l, Double_t k)
{
	//Default initialization
	sigma = xsection;
	luminosity = l;
	k_factor = k;
	xbinsize = (METMAX - METMIN)/METBINS;
	ybinsize = (HTMAX - HTMIN)/HTBINS;

	SL_monojet.ResizeTo(METBINS, HTBINS);			//Initialize matrices' dimensions
	SL_dijet.ResizeTo(METBINS, HTBINS);
	SL_trijet.ResizeTo(METBINS, HTBINS);
	SL_tetrajet.ResizeTo(METBINS, HTBINS);
	for (Int_t i = 0; i < 4; i++)
	{
		background_monojet[i].ResizeTo(METBINS, HTBINS);
		background_dijet[i].ResizeTo(METBINS, HTBINS);
		background_trijet[i].ResizeTo(METBINS, HTBINS);
		background_tetrajet[i].ResizeTo(METBINS, HTBINS);
	}
	signal_monojet.ResizeTo(METBINS, HTBINS);
	signal_dijet.ResizeTo(METBINS, HTBINS);
	signal_trijet.ResizeTo(METBINS, HTBINS);
	signal_tetrajet.ResizeTo(METBINS, HTBINS);

	//monojet_table = new TH2D("monojet_table", "Monojet", METBINS, METMIN, METMAX, HTBINS, HTMIN, HTMAX);
	//dijet_table = new TH2D("dijet_table", "Dijet", METBINS, METMIN, METMAX, HTBINS, HTMIN, HTMAX);
	//trijet_table = new TH2D("trijet_table", "Trijet", METBINS, METMIN, METMAX, HTBINS, HTMIN, HTMAX);
	//tetrajet_table = new TH2D("tetrajet_table", "Tetrajet", METBINS, METMIN, METMAX, HTBINS, HTMIN, HTMAX);

	//tetrajet_table = new TH2D("tetrajet_table", "Tetrajet", 10, 0, 2000, 10, 0, 3000);
	//tetrajet_table->Draw();

	TCanvas* c1 = new TCanvas;
	TCanvas* c2 = new TCanvas;
	TCanvas* c3 = new TCanvas;
	TCanvas* c4 = new TCanvas;
	c1->cd();
	monojet_table->Draw("TEXT");
	c2->cd();
	dijet_table->Draw("TEXT");
	c3->cd();
	trijet_table->Draw("TEXT");
	c4->cd();
	tetrajet_table->Draw("TEXT");

	//Int_t dummy;
	//cin>>dummy;
}

*/