X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSVertexerIons.cxx;h=65f3ddb0e967b4f84ca5c3fc7a465fa6dd7ee153;hb=b5a9f753d92077cec01150463deccec04904304a;hp=3e557ac97235fd6456dcf54735ea70feb38d1014;hpb=83028144a8b87d2decbc8bfa5ab1f99907442cbb;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSVertexerIons.cxx b/ITS/AliITSVertexerIons.cxx index 3e557ac9723..65f3ddb0e96 100644 --- a/ITS/AliITSVertexerIons.cxx +++ b/ITS/AliITSVertexerIons.cxx @@ -1,28 +1,16 @@ -#include -#include "stdlib.h" -#include -#include -#include -#include -#include -#include -#include - -#include "AliRun.h" -#include "AliITS.h" -#include "AliITSgeom.h" +#include "AliRunLoader.h" +#include "AliITSDetTypeRec.h" +#include "AliITSVertexerIons.h" +#include "AliITSVertexerZ.h" +#include "AliESDVertex.h" +#include "AliITSgeomTGeo.h" #include "AliITSLoader.h" #include "AliITSRecPoint.h" -#include "AliGenerator.h" -#include "AliMagF.h" - +#include "AliLog.h" +#include #include #include -#include -#include -#include -#include -#include + ////////////////////////////////////////////////////////////////////// // AliITSVertexerIons is a class for full 3D primary vertex // // finding optimized for Ion-Ion interactions // @@ -31,142 +19,55 @@ // // // // // Written by Giuseppe Lo Re and Francesco Riggi // -// // Giuseppe.Lore@ct.infn.it // // Franco.Riggi@ct.infn.it // // // -// Release date: May 2001 // +// Release date: Mar 2004 // // // // // ////////////////////////////////////////////////////////////////////// -ClassImp(AliITSVertexerIons) - - - //______________________________________________________________________ - AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() { +ClassImp(AliITSVertexerIons) + +//______________________________________________________________________ + AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(), +fNpThreshold(0), +fMaxDeltaPhi(0), +fMaxDeltaZ(0){ // Default Constructor - fITS = 0; + //fITS = 0; SetNpThreshold(); + SetMaxDeltaPhi(); + SetMaxDeltaZ(); } -//______________________________________________________________________ -AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) { - // Standard constructor - - fITS = 0; - SetNpThreshold(); -} - - //______________________________________________________________________ AliITSVertexerIons::~AliITSVertexerIons() { // Default Destructor - fITS = 0; + //fITS = 0; } //______________________________________________________________________ -AliITSVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){ - // Defines the AliITSVertex for the current event - fCurrentVertex = 0; - Double_t position[3]; - Double_t resolution[3]; - Double_t snr[3]; - AliRunLoader *rl = AliRunLoader::GetRunLoader(); - AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader"); - if(!fITS) { - fITS =(AliITS *)gAlice->GetDetector("ITS"); - if(!fITS) { - Error("FindVertexForCurrentEvent","AliITS object was not found"); - return fCurrentVertex; - } - } - fITS->SetTreeAddress(); - AliITSgeom *g2 = fITS->GetITSgeom(); - TClonesArray *recpoints = fITS->RecPoints(); - AliITSRecPoint *pnt; +AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(TTree *itsClusterTree){ +// Defines the AliESDVertex for the current event - Int_t nopoints1 = 0; - Int_t nopoints2 = 0; - Double_t vzero[3]; - Double_t aspar[2]; - - //------------ Rough Vertex evaluation --------------------------------- + fCurrentVertex = 0; - Int_t i,npoints,ipoint,j,k,max,binmax; - Double_t zCentroid; - Float_t l[3], p[3]; + AliITSDetTypeRec detTypeRec; - Double_t mxpiu = 0; - Double_t mxmeno = 0; - Double_t mypiu = 0; - Double_t mymeno = 0; - - TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35); - TTree *tr = itsloader->TreeR(); - for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) - { - fITS->ResetRecPoints(); - tr->GetEvent(i); - npoints = recpoints->GetEntries(); - for (ipoint=0;ipointUncheckedAt(ipoint); - l[0]=pnt->GetX(); - l[1]=0; - l[2]=pnt->GetZ(); - g2->LtoG(i, l, p); - if(i<80 && TMath::Abs(p[2])<14.35) { - if(p[0]>0) mxpiu++; - if(p[0]<0) mxmeno++; - if(p[1]>0) mypiu++; - if(p[1]<0) mymeno++; - hITSz1->Fill(p[2]); - } - if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++; - } - } - - nopoints1 = (Int_t)(hITSz1->GetEntries()); - - aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno); - aspar[1] = (mypiu-mymeno)/(mypiu+mymeno); - - vzero[0] = 5.24441*aspar[0]; - vzero[1] = 5.24441*aspar[1]; + detTypeRec.SetTreeAddressR(itsClusterTree); - zCentroid= TMath::Abs(hITSz1->GetMean()); - vzero[2] = -5.31040e-02+1.42198*zCentroid+7.44718e-01*TMath::Power(zCentroid,2) - -5.73426e-01*TMath::Power(zCentroid,3)+2.01500e-01*TMath::Power(zCentroid,4) - -3.34118e-02*TMath::Power(zCentroid,5)+2.20816e-03*TMath::Power(zCentroid,6); - - if(hITSz1->GetMean()<0) vzero[2] = -vzero[2]; - - /*cout << "\nXvzero: " << vzero[0] << " cm" << ""; - cout << "\nYvzero: " << vzero[1] << " cm" << ""; - cout << "\nZvzero: " << vzero[2] << " cm" << "\n";*/ + TClonesArray *recpoints = detTypeRec.RecPoints(); + AliITSRecPoint *pnt; - delete hITSz1; - Double_t dphi,r,deltaZ,a,b; - Int_t np1=0; - Int_t np2=0; - Int_t niter=0; - - Double_t deltaPhiZ = 0.08; - - Float_t mag[3]; - Float_t origin[3]; - for(Int_t okm=0;okm<3;okm++) origin[okm]=0; - gAlice->Field()->Field(origin,mag); - - deltaPhiZ = deltaPhiZ*mag[2]/2; - Double_t deltaPhiXY = 1.0; + Int_t npoints=0; + Int_t nopoints1=40000; + Int_t nopoints2=40000; + Float_t p[3]; - // cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n"; - // cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n"; - Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2; z1=new Double_t[nopoints1]; z2=new Double_t[nopoints2]; @@ -178,273 +79,160 @@ AliITSVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){ phi2=new Double_t[nopoints2]; r1=new Double_t[nopoints1]; r2=new Double_t[nopoints2]; - - deltaZ = 4.91617e-01+2.67567e-02*vzero[2]+1.49626e-02*TMath::Power(vzero[2],2); - Float_t multFactorZ = 28000./(Float_t)nopoints1; - Int_t nbin=(Int_t)((deltaZ/0.005)/multFactorZ); - Int_t nbinxy=250; - Int_t *vectorBinZ,*vectorBinXY; - vectorBinZ=new Int_t[nbin]; - vectorBinXY=new Int_t[nbinxy]; - Float_t f1= 0; - Float_t f2= 0; - Double_t sigma,averagebg; - - TH1D *hITSZv = new TH1D("hITSZv","",nbin,vzero[2]-deltaZ,vzero[2]+deltaZ); - TH1D *hITSXv = new TH1D("hITSXv","",nbinxy,-3,3); - TH1D *hITSYv = new TH1D("hITSYv","",nbinxy,-3,3); - - // cout << "deltaZeta: " << deltaZ << " cm" << "\n"; - - - start: - - hITSZv->Add(hITSZv,-1); - hITSXv->Add(hITSXv,-1); - hITSYv->Add(hITSYv,-1); - - np1=np2=0; - - - for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) - { - fITS->ResetRecPoints(); - itsloader->TreeR()->GetEvent(i); - npoints = recpoints->GetEntries(); - for (ipoint=0;ipointUncheckedAt(ipoint); - l[0]=pnt->GetX(); - l[1]=0; - l[2]=pnt->GetZ(); - g2->LtoG(i, l, p); - - if(i<80 && TMath::Abs(p[2])<14.35) { - p[0]=p[0]-vzero[0]; - p[1]=p[1]-vzero[1]; - r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); - y1[np1]=p[1]; - x1[np1]=p[0]; - z1[np1]=p[2]; - r1[np1]=r; - phi1[np1]=PhiFunc(p); - np1++; - } - if(i>=80 && TMath::Abs(p[2])<14.35) { - p[0]=p[0]-vzero[0]; - p[1]=p[1]-vzero[1]; - r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); - y2[np2]=p[1]; - x2[np2]=p[0]; - z2[np2]=p[2]; - r2[np2]=r; - phi2[np2]=PhiFunc(p); - np2++; - } - + Double_t mxpiu = 0; + Double_t mxmeno = 0; + Double_t mypiu = 0; + Double_t mymeno = 0; + Double_t r=0; + + Int_t np1=0, np2=0; + for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) { + detTypeRec.ResetRecPoints(); + itsClusterTree->GetEvent(i); + npoints = recpoints->GetEntries(); + for (Int_t ipoint=0;ipointUncheckedAt(ipoint); + /* + l[0]=pnt->GetDetLocalX(); + l[1]=0; + l[2]=pnt->GetDetLocalZ(); + g2->LtoG(i, l, p); + */ + pnt->GetGlobalXYZ(p); + r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); + + if(i<80 && TMath::Abs(p[2])<14.35) { + y1[np1]=p[1]; + x1[np1]=p[0]; + z1[np1]=p[2]; + if(p[0]>0) mxpiu++; + if(p[0]<0) mxmeno++; + if(p[1]>0) mypiu++; + if(p[1]<0) mymeno++; + np1++; + } + if(i>=80 && TMath::Abs(p[2])<14.35) { + y2[np2]=p[1]; + x2[np2]=p[0]; + z2[np2]=p[2]; + np2++; } } + } - //------------------ Correlation between rec points ---------------------- - - if(np1FindVertexForCurrentEvent(rl->GetEventNumber()); + AliITSVertexerZ *dovert = new AliITSVertexerZ(); + fCurrentVertex =dovert->FindVertexForCurrentEvent(itsClusterTree); delete dovert; return fCurrentVertex; } + if(!np1 || !np2) { + Error("FindVertexForCurrentEvent","No points in the pixer layers"); + return fCurrentVertex; + } + AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2)); - for(j=0; j<(np2)-1; j++) { - for(k=0; k<(np1)-1; k++) { - dphi=TMath::Abs(phi1[k]-phi2[j]); - if(dphi>180) dphi = 360-dphi; - if(dphiFill(z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1)); - } - } - - //cout << "\nNumber of used pairs: \n"; + Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno); + Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno); - a = vzero[2]-deltaZ; - b = vzero[2]+deltaZ; - max=(Int_t) hITSZv->GetMaximum(); - binmax=hITSZv->GetMaximumBin(); - sigma=0; - for(i=0;iGetBinContent(i); - for(i=0;i<10;i++) f1=f1+vectorBinZ[i]/10; - for(i=nbin-10;i(max-averagebg)*0.4 && - vectorBinZ[i]-averagebg<(max-averagebg)*0.7) { - sigma=hITSZv->GetBinCenter(binmax)-hITSZv->GetBinCenter(i); - sigma=TMath::Abs(sigma); - if(sigma==0) sigma=0.05; - } - } - - - TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b); - - fz->SetParameter(0,max); - if(niter==0) {Double_t temp = hITSZv->GetBinCenter(binmax); vzero[2]=temp;} - fz->SetParameter(1,vzero[2]); - fz->SetParameter(2,sigma); - fz->SetParameter(3,averagebg); - fz->SetParLimits(2,0,999); - fz->SetParLimits(3,0,999); - - hITSZv->Fit("fz","RQ0"); - - snr[2] = fz->GetParameter(0)/fz->GetParameter(3); - if(snr[2]<0.) { - Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n"); - Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n"); - // exit(123456789); - return fCurrentVertex; + Double_t x0 = 2.86*asparx; + Double_t y0 = 2.86*aspary; + + AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0)); + + for(Int_t i=0;iGetParameter(1); - if(position[2]<0) - { - position[2]=position[2]-TMath::Abs(position[2])*1.11/10000; - } - else - { - position[2]=position[2]+TMath::Abs(position[2])*1.11/10000; - } + for(Int_t i=0;iGetParError(1); - - for(j=0; j<(np2)-1; j++) { - for(k=0; k<(np1)-1; k++) { - dphi=TMath::Abs(phi1[k]-phi2[j]); + + Int_t nbinxy=400; + Int_t nbinz=400; + TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5); + TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5); + TH1F *hzv=new TH1F("hzv","",nbinz,-15,15); + + Double_t dphi; + for(Int_t j=0; j180) dphi = 360-dphi; - if(dphi>deltaPhiXY) continue; - if(TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-position[2]) - <4*resolution[2]) { - hITSXv->Fill(vzero[0]+(x2[j]-((x2[j]-x1[k])/(y2[j]-y1[k]))*y2[j])); - hITSYv->Fill(vzero[1]+(y2[j]-((y2[j]-y1[k])/(x2[j]-x1[k]))*x2[j])); - } - } + if(dphi>fMaxDeltaPhi) continue; + hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1)); + } } + + // Fitting ... + Double_t max; + Int_t binMax; + Double_t maxcenter; - TF1 *fx = new TF1("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[0]-0.5,vzero[0]+0.5); - - max=(Int_t) hITSXv->GetMaximum(); - binmax=hITSXv->GetMaximumBin(); - sigma=0; - f1=f2=0; - for(i=0;iGetBinContent(i); - for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10; - for(i=nbinxy-10;i(max-averagebg)*0.4 && - vectorBinXY[i]-averagebg<(max-averagebg)*0.7) { - sigma=hITSXv->GetBinCenter(binmax)-hITSXv->GetBinCenter(i); - sigma=TMath::Abs(sigma); - if(sigma==0) sigma=0.05; - } - } - - - fx->SetParameter(0,max); - fx->SetParameter(1,vzero[0]); - fx->SetParameter(2,sigma); - fx->SetParameter(3,averagebg); - - hITSXv->Fit("fx","RQ0"); - - snr[0] = fx->GetParameter(0)/fx->GetParameter(3); + max = hzv->GetMaximum(); + binMax=hzv->GetMaximumBin(); + maxcenter=hzv->GetBinCenter(binMax); + Double_t dxy=1.5; - if(snr[0]<0.) { - Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n"); - Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n"); - // exit(123456789); + TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy); + fz->SetParameter(0,max); + fz->SetParameter(1,maxcenter); + fz->SetParameter(2,0.1); + fz->SetLineColor(kRed); + hzv->Fit("fz","RQ0"); + + if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) { + Warning("FindVertexForCurrentEvent","AliITSVertexerIonsGeneral finder is not reliable for events with x and y vertex coordinates too far from the origin of the reference system. Their values will be set to 0.\n"); + Double_t position[3]={0,0,fz->GetParameter(1)}; + Double_t resolution[3]={0,0,0}; + Double_t snr[3]={0,0,0}; + Char_t name[30]; + sprintf(name,"Vertex"); + fCurrentVertex = new AliESDVertex(position,resolution,snr,name); return fCurrentVertex; } - else { - position[0]=fx->GetParameter(1); - resolution[0]=fx->GetParError(1); - } - TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[1]-0.5,vzero[1]+0.5); - - max=(Int_t) hITSYv->GetMaximum(); - binmax=hITSYv->GetMaximumBin(); - sigma=0; - f1=f2=0; - for(i=0;iGetBinContent(i); - for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10; - for(i=nbinxy-10;i(max-averagebg)*0.4 && - vectorBinXY[i]-averagebg<(max-averagebg)*0.7) { - sigma=hITSYv->GetBinCenter(binmax)-hITSYv->GetBinCenter(i); - sigma=TMath::Abs(sigma); - if(sigma==0) sigma=0.05; + for(Int_t j=0; jGetParameter(1))>fMaxDeltaZ) continue; + if(y2[k]==y1[j]) continue; + hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k])); + if(x2[k]==x1[j]) continue; + hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k])); } } + + TF1 *fx = new TF1 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",x0-dxy,x0+dxy); + fx->SetParameter(0,100); + Double_t dist=0.3; + Double_t xapprox=FindMaxAround(x0,hxv,dist); + AliDebug(1,Form("xapprox = %f",xapprox)); + fx->SetParameter(1,xapprox); + Double_t difcentroid=0.07; + fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid); + fx->SetParameter(2,0.1); + fx->SetLineColor(kRed); + hxv->Fit("fx","RQW0"); - fy->SetParameter(0,max); - fy->SetParameter(1,vzero[1]); - fy->SetParameter(2,sigma); - fy->SetParameter(3,averagebg); - - hITSYv->Fit("fy","RQ0"); - - snr[1] = fy->GetParameter(0)/fy->GetParameter(3); - if(snr[1]<0.) { - Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n"); - Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n"); - // exit(123456789); - return fCurrentVertex; - } - else { - position[1]=fy->GetParameter(1); - resolution[1]=fy->GetParError(1); - } + TF1 *fy = new TF1 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",y0-dxy,y0+dxy); + fy->SetParameter(0,100); + Double_t yapprox=FindMaxAround(y0,hyv,dist); + AliDebug(1,Form("yapprox = %f",yapprox)); + fy->SetParameter(1,yapprox); + fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid); + fy->SetParameter(2,0.1); + fy->SetLineColor(kRed); + hyv->Fit("fy","RQW0"); - - /*cout << "iter = " << niter << " -> x = " << position[0] << endl; - cout << "iter = " << niter << " -> y = " << position[1] << endl; - cout << "iter = " << niter << " -> z = " << position[2] << endl; - cout << "---\n";*/ - - niter++; - - vzero[0] = position[0]; - vzero[1] = position[1]; - vzero[2] = position[2]; - - if(niter<3) goto start; // number of iterations - - - cout<<"Number of iterations: "<GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)}; + Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)}; + Double_t snr[3]={0,0,0}; + Char_t name[30]; - // sprintf(name,"Vertex_%d",evnumber); - if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber); sprintf(name,"Vertex"); - fCurrentVertex = new AliITSVertex(position,resolution,snr,name); - return fCurrentVertex; -} - -//______________________________________________________________________ -Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) { - Double_t phi=0; + fCurrentVertex = new AliESDVertex(position,resolution,snr,name); - if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578); - if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180; - if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180; - if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360; - return phi; + return fCurrentVertex; } //______________________________________________________________________ -void AliITSVertexerIons::FindVertices(){ - // computes the vertices of the events in the range FirstEvent - LastEvent - AliRunLoader *rl = AliRunLoader::GetRunLoader(); - AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader"); - itsloader->LoadRecPoints("read"); - for(Int_t i=fFirstEvent;i<=fLastEvent;i++){ - rl->GetEvent(i); - FindVertexForCurrentEvent(i); - if(fCurrentVertex){ - WriteCurrentVertex(); - } - else { - if(fDebug>0){ - cout<<"Vertex not found for event "<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578); + if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180; + if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180; + if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;; } //________________________________________________________ void AliITSVertexerIons::PrintStatus() const { // Print current status cout <<"=======================================================\n"; - cout <<" Debug flag: "<PrintStatus(); } + +//________________________________________________________ +Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) { + // It finds the maximum of h within point-distance and distance+point. + Int_t maxcontent=0; + Int_t maxbin=0; + for(Int_t i=0;iGetNbinsX();i++) { + Int_t content=(Int_t)h->GetBinContent(i); + Double_t center=(Double_t)h->GetBinCenter(i); + if(TMath::Abs(center-point)>distance) continue; + if(content>maxcontent) {maxcontent=content;maxbin=i;} + } + Double_t max=h->GetBinCenter(maxbin); + return max; +} + +