[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx

/**************************************************************************
 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include <TMath.h>
#include "AliAODEvent.h"
#include "AliAODMCParticle.h"
#include "AliVertexingHFUtils.h"

/* $Id$ */

///////////////////////////////////////////////////////////////////
//                                                               //
// Class with functions useful for different D2H analyses        //
// - event plane resolution                                      //
// - <pt> calculation with side band subtraction                 //
// - tracklet multiplicity calculation                            //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

ClassImp(AliVertexingHFUtils)

//______________________________________________________________________
AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
  fK(1),
  fSubRes(1.),
  fMinEtaForTracklets(-1.),
  fMaxEtaForTracklets(1.)
{
  // Default contructor
}


//______________________________________________________________________
AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
  TObject(),
  fK(k),
  fSubRes(1.),
  fMinEtaForTracklets(-1.),
  fMaxEtaForTracklets(1.)
{
  // Standard constructor
}


//______________________________________________________________________
Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
  // compute chi from polynomial approximation
  Double_t c[5];
  if(k==1){ 
    c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  }
  else if(k==2){
    c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  }
  else return -1;
  return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
}

//______________________________________________________________________
Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
}


//______________________________________________________________________
Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
  // compute chi variable (=v2*sqrt(N)) from external values

  Double_t x1=0;
  Double_t x2=15;
  Double_t y1,y2;
  if(k==1){
    y1=ResolK1(x1)-res;
    y2=ResolK1(x2)-res;
  }
  else if(k==2){
    y1=Pol(x1,2)-res;
    y2=Pol(x2,2)-res;
  }
  else return -1;

  if(y1*y2>0) return -1;
  if(y1==0) return y1;
  if(y2==0) return y2;
  Double_t xmed,ymed;
  Int_t jiter=0;
  while((x2-x1)>0.0001){
    xmed=0.5*(x1+x2);
    if(k==1){
      y1=ResolK1(x1)-res;
      y2=ResolK1(x2)-res;
      ymed=ResolK1(xmed)-res;
    }
    else if(k==2){
      y1=Pol(x1,2)-res;
      y2=Pol(x2,2)-res;
      ymed=Pol(xmed,2)-res;
    }
    else return -1;
    if((y1*ymed)<0) x2=xmed;
    if((y2*ymed)<0)x1=xmed;
    if(ymed==0) return xmed;
    jiter++;
  }
  return 0.5*(x1+x2);
}

//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
  // computes event plane resolution starting from sub event resolution
  Double_t chisub=FindChi(resSub,k);
  Double_t chifull=chisub*TMath::Sqrt(2);
  if(k==1) return ResolK1(chifull);
  else if(k==2) return Pol(chifull,2);
  else{
    printf("k should be <=2\n");
    return 1.;
  }
}

//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
  // computes event plane resolution starting from sub event correlation histogram
  if(!hSubEvCorr) return 1.;
  Double_t resSub=GetSubEvResol(hSubEvCorr);
  return GetFullEvResol(resSub,k);
}
//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
  // computes low limit event plane resolution starting from sub event correlation histogram
  if(!hSubEvCorr) return 1.;
  Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
  printf("%f\n",resSub);
  return GetFullEvResol(resSub,k);  
}
//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
  // computes high limit event plane resolution starting from sub event correlation histogram
  if(!hSubEvCorr) return 1.;
  Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
  printf("%f\n",resSub);
  return GetFullEvResol(resSub,k);  
}
//______________________________________________________________________
Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
  // counts tracklets in given eta range
  AliAODTracklets* tracklets=ev->GetTracklets();
  Int_t nTr=tracklets->GetNumberOfTracklets();
  Int_t count=0;
  for(Int_t iTr=0; iTr<nTr; iTr++){
    Double_t theta=tracklets->GetTheta(iTr);
    Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
    if(eta>mineta && eta<maxeta) count++;
  }
  return count;
}
//______________________________________________________________________
Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
  // counts generated particles in fgiven eta range

  Int_t nChargedMC=0;
  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
    AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
    Int_t charge = part->Charge();
    Double_t eta = part->Eta();
    if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
  } 
  return nChargedMC;
}
//______________________________________________________________________
Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
  // counts generated primary particles in given eta range

  Int_t nChargedMC=0;
  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
    AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
    Int_t charge = part->Charge();
    Double_t eta = part->Eta();
    if(charge!=0 && eta>mineta && eta<maxeta){
      if(part->IsPrimary())nChargedMC++;
    } 
  }  
  return nChargedMC;
}
//______________________________________________________________________
Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
  // counts generated primary particles in given eta range

  Int_t nChargedMC=0;
  for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
    AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
    Int_t charge = part->Charge();
    Double_t eta = part->Eta();
    if(charge!=0 && eta>mineta && eta<maxeta){
      if(part->IsPhysicalPrimary())nChargedMC++;
    } 
  }
  return nChargedMC;
}
//______________________________________________________________________
void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Int_t rebin){

  // Compute <pt> from 2D histogram M vs pt

  //Make 2D histos in the desired pt range
  Int_t start=hMassD->FindBin(ptmin);
  Int_t end=hMassD->FindBin(ptmax)-1;
  const Int_t nx=end-start;
  TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY()));
  for(Int_t ix=start;ix<end;ix++){
    for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
      hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
      hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
    }
  }

  Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
  Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
  Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
  Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
  Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
  Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
  //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);

  Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
  Int_t minBinBkgLow=2;
  Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
  Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
  Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
  Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
  Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
  Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
  Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
  Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
  //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);

  Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
  Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
  Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
  //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);

  TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
  TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
  TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);

  hMptBkgLo->Rebin(rebin);
  hMptBkgHi->Rebin(rebin);
  hMptSigReg->Rebin(rebin);

  hMptBkgLo->Sumw2();
  hMptBkgHi->Sumw2();
  TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
  hMptBkgLoScal->Scale(bkgSig/bkgLow);
  TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
  hMptBkgHiScal->Scale(bkgSig/bkgHi);

  TH1F* hMptBkgAver=0x0;
  hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
  hMptBkgAver->Add(hMptBkgHiScal);
  hMptBkgAver->Scale(0.5);
  TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
  hMptSig->Add(hMptBkgAver,-1.);   
 
  averagePt = hMptSig->GetMean();
  errorPt = hMptSig->GetMeanError();

  delete hMptBkgLo;
  delete hMptBkgHi;
  delete hMptSigReg;
  delete hMptBkgLoScal;
  delete hMptBkgHiScal;
  delete hMptBkgAver;
  delete hMassDpt;
  delete hMptSig;

}
Commit	Line	Data
a6c5a2e9	1	/**************************************************************************
	2	* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	#include <TMath.h>
	17	#include "AliAODEvent.h"
25504150	18	#include "AliAODMCParticle.h"
a6c5a2e9	19	#include "AliVertexingHFUtils.h"
a6c5a2e9	20
d947ea9c	21	/* $Id$ */
a6c5a2e9	22
	23	///////////////////////////////////////////////////////////////////
	24	// //
	25	// Class with functions useful for different D2H analyses //
	26	// - event plane resolution //
	27	// - <pt> calculation with side band subtraction //
25504150	28	// - tracklet multiplicity calculation //
a6c5a2e9	29	// Origin: F.Prino, Torino, prino@to.infn.it //
	30	// //
	31	///////////////////////////////////////////////////////////////////
	32
	33	ClassImp(AliVertexingHFUtils)
	34
	35	//______________________________________________________________________
	36	AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
	37	fK(1),
	38	fSubRes(1.),
	39	fMinEtaForTracklets(-1.),
	40	fMaxEtaForTracklets(1.)
	41	{
	42	// Default contructor
	43	}
	44
	45
	46	//______________________________________________________________________
	47	AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
	48	TObject(),
	49	fK(k),
	50	fSubRes(1.),
	51	fMinEtaForTracklets(-1.),
	52	fMaxEtaForTracklets(1.)
	53	{
	54	// Standard constructor
	55	}
	56
	57
	58	//______________________________________________________________________
	59	Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
	60	// compute chi from polynomial approximation
	61	Double_t c[5];
	62	if(k==1){
	63	c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
	64	}
	65	else if(k==2){
	66	c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
	67	}
	68	else return -1;
	69	return c[0]x+c[1]xx+c[2]xxx+c[3]xxxx+c[4]xxxx*x;
	70	}
	71
	72	//______________________________________________________________________
	73	Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
	74	return TMath::Sqrt(TMath::Pi()/8)xTMath::Exp(-xx/4)(TMath::BesselI0(xx/4)+TMath::BesselI1(xx/4));
	75	}
	76
	77
	78	//______________________________________________________________________
	79	Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){
	80	// compute chi variable (=v2*sqrt(N)) from external values
	81
	82	Double_t x1=0;
	83	Double_t x2=15;
	84	Double_t y1,y2;
	85	if(k==1){
	86	y1=ResolK1(x1)-res;
	87	y2=ResolK1(x2)-res;
	88	}
	89	else if(k==2){
	90	y1=Pol(x1,2)-res;
	91	y2=Pol(x2,2)-res;
	92	}
93	else return -1;
94
95	if(y1*y2>0) return -1;
96	if(y1==0) return y1;
97	if(y2==0) return y2;
98	Double_t xmed,ymed;
99	Int_t jiter=0;
100	while((x2-x1)>0.0001){
101	xmed=0.5*(x1+x2);
102	if(k==1){
103	y1=ResolK1(x1)-res;
104	y2=ResolK1(x2)-res;
105	ymed=ResolK1(xmed)-res;
106	}
107	else if(k==2){
108	y1=Pol(x1,2)-res;
109	y2=Pol(x2,2)-res;
110	ymed=Pol(xmed,2)-res;
111	}
112	else return -1;
113	if((y1*ymed)<0) x2=xmed;
114	if((y2*ymed)<0)x1=xmed;
115	if(ymed==0) return xmed;
116	jiter++;
117	}
118	return 0.5*(x1+x2);
119	}
120
121	//______________________________________________________________________
122	Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
123	// computes event plane resolution starting from sub event resolution
124	Double_t chisub=FindChi(resSub,k);
125	Double_t chifull=chisub*TMath::Sqrt(2);
126	if(k==1) return ResolK1(chifull);
127	else if(k==2) return Pol(chifull,2);
128	else{
129	printf("k should be <=2\n");
130	return 1.;
131	}
132	}
133
134	//______________________________________________________________________
135	Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
136	// computes event plane resolution starting from sub event correlation histogram
137	if(!hSubEvCorr) return 1.;
138	Double_t resSub=GetSubEvResol(hSubEvCorr);
139	return GetFullEvResol(resSub,k);
140	}
141	//______________________________________________________________________
142	Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
143	// computes low limit event plane resolution starting from sub event correlation histogram
144	if(!hSubEvCorr) return 1.;
145	Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
146	printf("%f\n",resSub);
147	return GetFullEvResol(resSub,k);
148	}
149	//______________________________________________________________________
150	Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
151	// computes high limit event plane resolution starting from sub event correlation histogram
152	if(!hSubEvCorr) return 1.;
153	Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
154	printf("%f\n",resSub);
155	return GetFullEvResol(resSub,k);
156	}
157	//______________________________________________________________________
158	Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
159	// counts tracklets in given eta range
160	AliAODTracklets* tracklets=ev->GetTracklets();
161	Int_t nTr=tracklets->GetNumberOfTracklets();
162	Int_t count=0;
163	for(Int_t iTr=0; iTr<nTr; iTr++){
164	Double_t theta=tracklets->GetTheta(iTr);
165	Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
166	if(eta>mineta && eta<maxeta) count++;
167	}
168	return count;
169	}
86d84fa4	170	//______________________________________________________________________
25504150	171	Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
	172	// counts generated particles in fgiven eta range
	173
	174	Int_t nChargedMC=0;
	175	for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
	176	AliAODMCParticle part=(AliAODMCParticle)arrayMC->UncheckedAt(i);
	177	Int_t charge = part->Charge();
	178	Double_t eta = part->Eta();
	179	if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
	180	}
	181	return nChargedMC;
	182	}
	183	//______________________________________________________________________
	184	Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
	185	// counts generated primary particles in given eta range
	186
	187	Int_t nChargedMC=0;
	188	for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
	189	AliAODMCParticle part=(AliAODMCParticle)arrayMC->UncheckedAt(i);
	190	Int_t charge = part->Charge();
	191	Double_t eta = part->Eta();
	192	if(charge!=0 && eta>mineta && eta<maxeta){
	193	if(part->IsPrimary())nChargedMC++;
	194	}
	195	}
	196	return nChargedMC;
	197	}
	198	//______________________________________________________________________
	199	Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
	200	// counts generated primary particles in given eta range
	201
	202	Int_t nChargedMC=0;
	203	for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
	204	AliAODMCParticle part=(AliAODMCParticle)arrayMC->UncheckedAt(i);
	205	Int_t charge = part->Charge();
	206	Double_t eta = part->Eta();
	207	if(charge!=0 && eta>mineta && eta<maxeta){
	208	if(part->IsPhysicalPrimary())nChargedMC++;
	209	}
	210	}
	211	return nChargedMC;
	212	}
	213	//______________________________________________________________________
86d84fa4	214	void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Int_t rebin){
	215
	216	// Compute <pt> from 2D histogram M vs pt
	217
	218	//Make 2D histos in the desired pt range
	219	Int_t start=hMassD->FindBin(ptmin);
	220	Int_t end=hMassD->FindBin(ptmax)-1;
	221	const Int_t nx=end-start;
	222	TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY()));
	223	for(Int_t ix=start;ix<end;ix++){
	224	for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
	225	hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
	226	hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
	227	}
	228	}
	229
	230	Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
	231	Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
	232	Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
	233	Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
	234	Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
	235	Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
	236	// printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
	237
	238	Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
	239	Int_t minBinBkgLow=2;
	240	Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
	241	Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
	242	Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
	243	Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
	244	Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
	245	Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
	246	Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
	247	Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
	248	// printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
	249
	250	Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
	251	Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
	252	Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
	253	// printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
	254
	255	TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
	256	TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
	257	TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
	258
	259	hMptBkgLo->Rebin(rebin);
	260	hMptBkgHi->Rebin(rebin);
	261	hMptSigReg->Rebin(rebin);
	262
	263	hMptBkgLo->Sumw2();
	264	hMptBkgHi->Sumw2();
	265	TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
	266	hMptBkgLoScal->Scale(bkgSig/bkgLow);
	267	TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
	268	hMptBkgHiScal->Scale(bkgSig/bkgHi);
	269
	270	TH1F* hMptBkgAver=0x0;
	271	hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
	272	hMptBkgAver->Add(hMptBkgHiScal);
	273	hMptBkgAver->Scale(0.5);
	274	TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
	275	hMptSig->Add(hMptBkgAver,-1.);
	276
	277	averagePt = hMptSig->GetMean();
278	errorPt = hMptSig->GetMeanError();
279
280	delete hMptBkgLo;
281	delete hMptBkgHi;
282	delete hMptSigReg;
283	delete hMptBkgLoScal;
284	delete hMptBkgHiScal;
285	delete hMptBkgAver;
286	delete hMassDpt;
287	delete hMptSig;
288
289	}