[u/mrichter/AliRoot.git] / PWG1 / AliIntSpotEstimator.cxx

#include "AliIntSpotEstimator.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliVertexerTracks.h"
#include "AliLog.h"
#include <TH1.h>
#include <TH1F.h>
#include <TNtuple.h>
#include <TObjArray.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TStyle.h>


ClassImp(AliIntSpotEstimator)

//______________________________________________________________________________________
AliIntSpotEstimator::AliIntSpotEstimator(const char* name,Double_t outcut,Int_t ntrIP,
					 Int_t nPhiBins,Int_t nestb,
					 Double_t estmin,Double_t estmax,
					 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
					 Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple) 
: TNamed(name,""),
  fEvProc(0),fIPCenterStat(0),fMinTracksForIP(ntrIP>2?ntrIP:2),fOutlierCut(outcut),fEstimIP(0),
  fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
{
  InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
}

//______________________________________________________________________________________                                    
AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef) 
  : TNamed("IPEstimator",""),
    fEvProc(0),fIPCenterStat(0),fMinTracksForIP(2),fOutlierCut(1e-4),
    fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
{
  if (initDef) InitEstimators();
  //
}

//______________________________________________________________________________________
AliIntSpotEstimator::~AliIntSpotEstimator()
{
  delete fEstimIP;
  delete fEstimVtx;
  delete fEstimTrc;
  delete fVertexer;
  delete fNtuple;
}

//______________________________________________________________________________________
Bool_t AliIntSpotEstimator::ProcessIPCenter(const AliESDVertex* vtx)
{
  // account the vertex in IP center estimation
  if (!vtx) return kFALSE;
  double xyz[3];
  vtx->GetXYZ(xyz);
  double r = xyz[0]*xyz[0] + xyz[1]*xyz[1];
  if (r>2.) return kFALSE;
  for (int i=3;i--;) {
    fIPCenter[i] += xyz[i];
    fIPCen2[i] += xyz[i]*xyz[i];
  }
  fIPCenterStat++;
  fHVtxXY->Fill(xyz[0],xyz[1]);
  return kTRUE;
  //
}

//______________________________________________________________________________________
Bool_t AliIntSpotEstimator::ProcessEvent(const AliESDEvent* esd, const AliESDVertex* vtx)
{
  // Note: this method will modify the indices of the provided vertex
  // Account the vertex in widths estimation
  // if vtx is provided its tracks will be used otherwise use the primary vertex of the event
  if (!IsValid()) {AliError("Not initialized yet"); return -999;}
  if (!fTracks)   fTracks   = new TObjArray(50);
  if (!fVertexer) fVertexer = new AliVertexerTracks();
  //
  fEvProc++;
  int nTracks;
  if (!vtx) vtx = (AliESDVertex*) esd->GetPrimaryVertex();
  if (!vtx || (nTracks=vtx->GetNIndices())<GetMinTracks()+1) return kFALSE;
  //
  UShort_t *trackID = (UShort_t*)vtx->GetIndices(); 
  for (int itr=0;itr<nTracks;itr++) fTracks->Add(esd->GetTrack(trackID[itr]));
  //
  fVertexer->SetFieldkG( esd->GetMagneticField() );
  return ProcessTracks();
  //
}

//______________________________________________________________________________________
Bool_t AliIntSpotEstimator::ProcessEvent(const TObjArray* tracks)
{
  // This method allows to process the vertex made of arbitrary tracks
  // Account the vertex in widths estimation with externally provided tracks 
  // The indices of the vertex must correspond to array of tracks
  //
  if (!IsValid()) {AliError("Not initialized yet"); return -999;}
  if (!fVertexer) fVertexer = new AliVertexerTracks();
  //
  fEvProc++;
  int nTracks = tracks->GetEntriesFast();
  if ( nTracks<GetMinTracks()+1 ) return kFALSE;
  for (int itr=0;itr<nTracks;itr++) fTracks->Add(tracks->At(itr));
  //
  return ProcessTracks();
}

//______________________________________________________________________________________
Bool_t AliIntSpotEstimator::ProcessTracks()
{
  // account the vertex made of provided tracks
  //
  int nTracks = fTracks->GetEntriesFast();
  UShort_t *trackID = new UShort_t[nTracks];
  for (int i=nTracks;i--;) trackID[i] = i;
  double xyzDCA[3],ddca[2],covdca[3];
  int nTracks1 = nTracks - 1;
  AliExternalTrackParam *selTrack,*movTrack=0;
  UShort_t selTrackID,movTrackID=0;
  //
  AliESDVertex* recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
  if (!recNewVtx || ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
    if (recNewVtx) delete recNewVtx;
    fTracks->Clear();
    delete[] trackID;
    return kFALSE;
  }
  if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
  //
  double pmn = GetTrackMinP();
  double pmx = GetTrackMaxP();
  double fieldVal = fVertexer->GetFieldkG();
  for (int itr=0;itr<nTracks;itr++) {
    selTrack   = (AliExternalTrackParam*) (*fTracks)[itr]; // track to probe
    double pTrack = selTrack->GetP();
    if (!IsZero(fieldVal) && (pTrack<pmn || pTrack>pmx)) continue;
    selTrackID = trackID[itr];
    //
    if (itr<nTracks1) {
      movTrack   = (AliExternalTrackParam*) (*fTracks)[nTracks1];   // save the track
      movTrackID = trackID[nTracks1];
      (*fTracks)[itr] = movTrack;  
      trackID[itr] = movTrackID;
    }
    fTracks->RemoveAt(nTracks1);     // move the last track in the place of the probed one
    //
    // refit the vertex w/o the probed track
    recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
    if (recNewVtx) {
      double told = selTrack->GetX();  // store the original track position
      selTrack->PropagateToDCA(recNewVtx,fieldVal,1e4,ddca,covdca); // in principle, done in the vertexer
      selTrack->GetXYZ(xyzDCA);
      //
      double phiTrack = selTrack->Phi();
      double cs = TMath::Cos(phiTrack);
      double sn = TMath::Sin(phiTrack);
      double trDCA = (xyzDCA[0]-fIPCenIni[0])         *sn - (xyzDCA[1]-fIPCenIni[1])         *cs;  // track signed DCA to origin
      double vtDCA = (recNewVtx->GetXv()-fIPCenIni[0])*sn - (recNewVtx->GetYv()-fIPCenIni[1])*cs;  // vertex signed DCA to origin
      UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
      selTrack->PropagateTo(told,fieldVal);    // restore the track
      if (fNtuple) {
	static float ntf[8];
	ntf[0] = float(nTracks1);
	ntf[1] = recNewVtx->GetXv();
	ntf[2] = recNewVtx->GetYv();
	ntf[3] = recNewVtx->GetZv();
	ntf[4] = xyzDCA[0];
	ntf[5] = xyzDCA[1];
	ntf[6] = phiTrack;
	ntf[7] = pTrack;
	fNtuple->Fill(ntf);
      }
    }
    delete recNewVtx;
    // restore the track indices
    (*fTracks)[itr] = selTrack; 
    trackID[itr] = selTrackID;
    if (itr<nTracks1) {
      (*fTracks)[nTracks1] = movTrack; 
      trackID[nTracks1] = movTrackID;
    }
  }
  //
  fTracks->Clear();
  delete[] trackID;
  return kTRUE;
  //
}

//______________________________________________________________________________________
void AliIntSpotEstimator::UpdateEstimators(double rvD, double rtD, double nTracks, double pTrack, double phiTrack)
{
  // update the estimator values
  double estIP  = rvD*rtD;
  double estVtx = rvD*(rvD - rtD);
  double estTrc = rtD*(rtD - rvD);
  //
  if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
  fEstimVtx->Fill(nTracks, estVtx);
  if (pTrack<1e-6) pTrack = GetTrackMinP()+1e6;
  fEstimTrc->Fill(1./pTrack,estTrc);
  //
}
 
//______________________________________________________________________________________
void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t estmin,Double_t estmax,
					 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
					 Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
{
  Clear();
  // regularize binning/limits for DCA resolution
  nPBins = nPBins<1 ? 1: nPBins;
  pmn    = pmn>0.1  ? pmn : 0.1;
  pmx    = pmx>pmn  ? pmx : pmn+0.1;
  //
  // regularize binning/limits for vertex resolution
  ntMn = ntMn>2 ? ntMn : 2;
  ntMx = ntMx>ntMn ? ntMx : ntMn;
  ntrBins = ntrBins<1 ? 1:ntrBins;
  int step = (ntMx-ntMn)/ntrBins;
  if (step<1) step = 1;
  ntrBins = (ntMx-ntMn)/step;
  //
  // regularize binning for IP sigma vs angle
  nPhiBins = nPhiBins>1 ? nPhiBins : 1;
  //
  // regularize binning/limits for estimators
  nestb  = nestb>300 ? nestb:300;
  estmin = estmin<-2.e-2 ? estmin : -2.e-2;
  estmax = estmax> 4.e-2 ? estmax :  4.e-2;
  //
  TString nm;
  nm = GetName();  
  nm += "diamondEst";
  fEstimIP  = new TH2F(nm.Data(),nm.Data(),nPhiBins,0.,2.*TMath::Pi(),nestb,estmin,estmax);
  //  
  nm = GetName();  
  nm += "VResEst";  
  fEstimVtx = new TH2F(nm.Data(),nm.Data(),ntrBins,ntMn,ntMx,         nestb,estmin,estmax);
  // 
  nm = GetName();  
  nm += "dcaEst";   
  fEstimTrc = new TH2F(nm.Data(),nm.Data(),nPBins,1./pmx,1./pmn,      nestb,estmin,estmax);
  // 
  nm = GetName();  
  nm += "VtxXY";   
  fHVtxXY = new TH2F(nm.Data(),nm.Data(),200, -1,1, 200,-1,1);
  //
  if (ntuple) {
    nm = GetName();  
    nm += "ntuple";
    fNtuple = new TNtuple(nm.Data(),nm.Data(),"ntrack:xv:yv:zv:xt:yt:phi:p");
  }
  //
  fVertexer = new AliVertexerTracks(); // the field is set dynamically
  fVertexer->SetConstraintOff();
  fTracks = new TObjArray(50);
  //

}

//______________________________________________________________________________________                                                                                             
Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
{                                                                                                                                                                                    
  // get estimate for the IP sigma
  if (!IsValid()) {AliError("Not initialized yet"); return -999;}
  double cxe,cye;
  double cx  = GetIPCenter(0,&cxe) - GetIPCenIni(0);
  double cy  = GetIPCenter(1,&cye) - GetIPCenIni(1);
  TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
  double merr = 0;
  double est = CalcMean(proj, fOutlierCut, &merr) - (cx*cx + cy*cy)/2.;
  if (est>0) {
    est = TMath::Sqrt(est);
    if (err) {
      *err = 0;
      *err = merr*merr;
      *err += cx*cx*cxe*cxe + cy*cy*cye*cye;
      *err = TMath::Sqrt(*err)/est/2.;
    }
  }
  else {
    est = 0;
    if (err) *err = 0;
  }
  delete proj;
  return est;
}

//______________________________________________________________________________________                                                                                             
Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
{                                                                                                                                                                                    
  // get estimate for the IP sigma^2
  if (!IsValid()) {AliError("Not initialized yet"); return -999;}
  int bin = fEstimVtx->GetXaxis()->FindBin(ntr);
  if (bin<1 || bin>GetNTrackBins()) {
    AliError(Form("Requested vertex multiplicity %d out of defined %d-%d range",
		  ntr,GetMinTracks(),GetMaxTracks()));
    return -1;
  }
  TH1* proj = fEstimVtx->ProjectionY("vrProj",bin,bin,"e");
  double est = CalcMean(proj, fOutlierCut, err);
  delete proj;
  if (est>0) {
    est = TMath::Sqrt(est);
    if (err) *err /= 2*est;
  }
  else {
    est = 0;
    if (err) *err = 0;
  }
  return est;
  //
}

//______________________________________________________________________________________                                                                                             
Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
{                                                                                                                                                                                    
  // get estimate for the IP sigma^2
  if (!IsValid()) {AliError("Not initialized yet"); return -999;}
  if (pt<1e-6) pt = GetTrackMinP()+1e6;
  pt = 1./pt;
  int bin = fEstimTrc->GetXaxis()->FindBin(pt);
  if (bin<1 || bin>GetNPBins()) {
    AliError(Form("Requested track P %.2f out of defined %.2f-%.2f range",1/pt,GetTrackMinP(),GetTrackMaxP()));
    return -1;
  }
  TH1* proj = fEstimTrc->ProjectionY("trProj",bin,bin,"e");
  double est = CalcMean(proj, fOutlierCut, err);
  delete proj;
   if (est>0) {
    est = TMath::Sqrt(est);
    if (err) *err /= 2*est;
  }
  else {
    est = 0;
    if (err) *err = 0;
  }
  return est;
  //
}

//______________________________________________________________________________________                                                                                             
Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
{
  // calculate mean applying the cut on the outliers
  //
  double max = histo->GetMaximum();
  double cut = (ctfact>0&&ctfact<1.) ? max*ctfact : 0;//TMath::Max(1.0,max*ctfact) : 0;
  int nb = histo->GetNbinsX();
  double mean = 0.,cumul = 0, rms = 0;
  for (int i=1;i<=nb;i++) {
    double vl = histo->GetBinContent(i) - cut;
    if (vl<1e-10) continue;
    double x = histo->GetBinCenter(i);
    mean  += vl*x;
    rms   += vl*x*x;
    cumul += vl;
  }
  //
  mean = cumul>0 ? mean/cumul : 0;
  rms -= mean*mean*cumul;
  if (err) {
    *err = cumul > 1 ? rms/(cumul-1) : 0;
    if (*err>0) *err = TMath::Sqrt(*err/cumul);
  }
  //
  return mean;
}

//______________________________________________________________________________________
void AliIntSpotEstimator::Print(Option_t *) const                                                                                                                              
{   
  if (!IsValid()) {printf("Not initialized yet\n"); return;}
  //
  double cx,cy,cz,cxe,cye,cze;
  cx = GetIPCenter(0,&cxe);
  cy = GetIPCenter(1,&cye);
  cz = GetIPCenter(2,&cze);
  printf("Processed %d events\n",fEvProc);
  printf("Estimator for IP center: %+.4e+-%.3e | %+.4e+-%.3e | %+.4e+-%.3e\n",
	 cx,cxe,cy,cye,cz,cze);
  printf("Initial IP center was  : %+.4e %+.4e %+.4e\n",
	 GetIPCenIni(0),GetIPCenIni(1),GetIPCenIni(2));
  double sgIP = GetIPSigma(0,&cxe);
  printf("Estimator for IP sigma : %.4e+-%.3e\n",sgIP,cxe);
  //
  printf("Estimators for vertex resolution vs Ntracks:\n");
  
  for (int i=1;i<=GetNTrackBins();i++) {
    double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
    if (IsZero(sig)) continue;
    int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
    int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
    printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
  }
  //
  printf("Estimators for track DCA resolution vs P:\n");
  for (int i=1;i<=GetNPBins();i++) {
    double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
    if (IsZero(sig)) continue;
    double pmax = 1./fEstimTrc->GetXaxis()->GetBinLowEdge(i);
    double pmin = 1./(fEstimTrc->GetXaxis()->GetBinLowEdge(i)+fEstimTrc->GetXaxis()->GetBinWidth(i));
    printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
  }
}

//______________________________________________________________________________________
void AliIntSpotEstimator::Clear(Option_t *)                                                                                   
{  
  // reset all
  fEvProc = 0;
  fIPCenterStat = 0;
  for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
  if (fEstimIP)  fEstimIP->Reset();
  if (fEstimVtx) fEstimVtx->Reset();
  if (fEstimTrc) fEstimTrc->Reset();
}

//_____________________________________________________________________
AliIntSpotEstimator &AliIntSpotEstimator::operator += (const AliIntSpotEstimator &src)
{
  fEvProc        += src.fEvProc;
  fIPCenterStat  += src.fIPCenterStat;
  for (int i=3;i--;) fIPCenter[i] += src.fIPCenter[i];
  if (fEstimIP  && src.fEstimIP ) fEstimIP->Add(src.fEstimIP);
  if (fEstimVtx && src.fEstimVtx) fEstimVtx->Add(src.fEstimVtx);
  if (fEstimTrc && src.fEstimTrc) fEstimTrc->Add(src.fEstimTrc);
  if (fHVtxXY   && src.fHVtxXY)   fHVtxXY->Add(src.fHVtxXY);
  //
  return *this;
}

//_____________________________________________________________________
TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
{
  // prepare canvas with report
  TCanvas *cnv = new TCanvas(GetName(), GetName(),5,5,700,1000);
  gStyle->SetOptStat(0);
  gStyle->SetTitleH(0.07);
  gStyle->SetTitleW(0.7);
  gStyle->SetTitleY(1);
  gStyle->SetTitleX(0.2);
  //
  const Int_t nc=200;
  char buff[nc];
  //
  cnv->Divide(2,2);
  cnv->cd(1);
  //
  TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
  snprintf(buff,nc,"%s | Outliers Cut : %.2e",GetName(),fOutlierCut);
  pt->AddText(buff);
  //
  snprintf(buff,nc,"Processed Events:\n%d",fEvProc);
  pt->AddText(buff);
  //
  snprintf(buff,nc,"Accepted  Events\n%d",fIPCenterStat);
  pt->AddText(buff);
  //
  double cx,cy,cz,cxe,cye,cze;
  cx = GetIPCenter(0,&cxe);
  cy = GetIPCenter(1,&cye);
  cz = GetIPCenter(2,&cze);
  //  
  snprintf(buff,nc,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
	  int(cx*1e4),int(cxe*1e4),int(cy*1e4),int(cye*1e4),int(cz*1e4),int(cze*1e4));
  pt->AddText(buff);
  //
  cx = GetIPSigma(0,&cxe);
  snprintf(buff,nc,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx*1e4),int(cxe*1e4));
  pt->AddText(buff);
  pt->Draw();
  gPad->Modified();
  //
  cnv->cd(2);
  gPad->SetLeftMargin(0.2);
  TH1* iph = fEstimIP->ProjectionX();
  iph->Reset();
  iph->SetTitle("Int.Spot size vs #phi");
  for (int i=1;i<=iph->GetNbinsX();i++) {
    cx = GetIPSigma(i,&cxe);
    iph->SetBinContent(i,cx*1e4);
    iph->SetBinError(i,cxe*1e4);
  }
  iph->GetXaxis()->SetTitle("#phi");
  iph->GetYaxis()->SetTitle("#sigma_{IP} [#mum]");
  iph->SetMarkerStyle(20);
  iph->Draw("p");
  gPad->Modified();
  iph->SetTitleOffset(2.5,"Y");
  //
  cnv->cd(3);
  gPad->SetLeftMargin(0.2);
  TH1* vrh = fEstimVtx->ProjectionX();
  vrh->Reset();
  vrh->SetTitle("Vertex resolution vs N tracks");
  for (int i=1;i<=vrh->GetNbinsX();i++) {
    cx = GetVtxSigma( TMath::Nint(vrh->GetBinCenter(i)), &cxe);
    vrh->SetBinContent(i,cx*1e4);
    vrh->SetBinError(i,cxe*1e4);
  }
  vrh->GetXaxis()->SetTitle("n tracks");
  vrh->GetYaxis()->SetTitle("#sigma_{VTX} [#mum]");
  vrh->SetMarkerStyle(20);
  vrh->Draw("p");
  gPad->Modified();
  vrh->SetTitleOffset(2.5,"Y");
  //
  cnv->cd(4);
  gPad->SetLeftMargin(0.2);
  TH1* trh = fEstimTrc->ProjectionX();
  trh->Reset();
  trh->SetTitle("Track DCA resolution vs 1/P");
  for (int i=1;i<=trh->GetNbinsX();i++) {
    cx = GetDCASigma(1./trh->GetBinCenter(i), &cxe);
    trh->SetBinContent(i,cx*1e4);
    trh->SetBinError(i,cxe*1e4);
  }
  trh->GetXaxis()->SetTitle("1/p [GeV]");
  trh->GetYaxis()->SetTitle("#sigma_{DCA} [#mum]");
  trh->SetMarkerStyle(20);
  trh->Draw("p");
  gPad->Modified();
  trh->SetTitleOffset(2.5,"Y");
  //
  cnv->cd();
  gPad->Modified();
  //
  if (outname) cnv->Print(outname);
  //
  return cnv;
}


//_____________________________________________________________________
Long64_t AliIntSpotEstimator::Merge(TCollection *coll)
{
  // Merge estimators: used to put together the results of parallel processing
  if(!coll) return 0;
  if (coll->IsEmpty()) return 1;
  TIterator* iter = coll->MakeIterator();
  TObject* obj;
  int count = 0;
  while ((obj = iter->Next())) {
    AliIntSpotEstimator* entry = dynamic_cast<AliIntSpotEstimator*>(obj);
    if (!entry) continue;
    (*this) += *entry;
    count++;
  }
  return count;
  //
}

//_____________________________________________________________________
Double_t AliIntSpotEstimator::GetIPCenter(Int_t id,Double_t *err) const
{
  // calculate IP center in axis id and error
  double cen = fIPCenterStat>0 ? fIPCenter[id]/fIPCenterStat:0;
  if (err) {
    *err = fIPCenterStat>1 ? (fIPCen2[id] - cen*cen*fIPCenterStat)/(fIPCenterStat-1) : 0;
    *err = *err > 0 ? TMath::Sqrt(*err/fIPCenterStat) : 0;
  }
  return cen;
}
Commit	Line	Data
5b09c01f	1	#include "AliIntSpotEstimator.h"
	2	#include "AliESDEvent.h"
	3	#include "AliESDtrack.h"
	4	#include "AliESDVertex.h"
	5	#include "AliVertexerTracks.h"
	6	#include "AliLog.h"
	7	#include <TH1.h>
	8	#include <TH1F.h>
d631e5e7	9	#include <TNtuple.h>
5b09c01f	10	#include <TObjArray.h>
	11	#include <TCanvas.h>
	12	#include <TPaveText.h>
	13	#include <TStyle.h>
	14
	15
	16	ClassImp(AliIntSpotEstimator)
	17
	18	//______________________________________________________________________________________
d631e5e7	19	AliIntSpotEstimator::AliIntSpotEstimator(const char* name,Double_t outcut,Int_t ntrIP,
5b09c01f	20	Int_t nPhiBins,Int_t nestb,
	21	Double_t estmin,Double_t estmax,
	22	Int_t ntrBins,Int_t ntMn,Int_t ntMx,
d631e5e7	23	Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
5b09c01f	24	: TNamed(name,""),
d631e5e7	25	fEvProc(0),fIPCenterStat(0),fMinTracksForIP(ntrIP>2?ntrIP:2),fOutlierCut(outcut),fEstimIP(0),
d631e5e7	26	fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
5b09c01f	27	{
d631e5e7	28	InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
5b09c01f	29	}
	30
	31	//______________________________________________________________________________________
	32	AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef)
	33	: TNamed("IPEstimator",""),
d631e5e7	34	fEvProc(0),fIPCenterStat(0),fMinTracksForIP(2),fOutlierCut(1e-4),
d631e5e7	35	fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
5b09c01f	36	{
	37	if (initDef) InitEstimators();
	38	//
	39	}
	40
	41	//______________________________________________________________________________________
	42	AliIntSpotEstimator::~AliIntSpotEstimator()
	43	{
	44	delete fEstimIP;
	45	delete fEstimVtx;
	46	delete fEstimTrc;
	47	delete fVertexer;
d631e5e7	48	delete fNtuple;
5b09c01f	49	}
	50
	51	//______________________________________________________________________________________
	52	Bool_t AliIntSpotEstimator::ProcessIPCenter(const AliESDVertex* vtx)
	53	{
	54	// account the vertex in IP center estimation
	55	if (!vtx) return kFALSE;
	56	double xyz[3];
	57	vtx->GetXYZ(xyz);
	58	double r = xyz[0]xyz[0] + xyz[1]xyz[1];
d631e5e7	59	if (r>2.) return kFALSE;
	60	for (int i=3;i--;) {
	61	fIPCenter[i] += xyz[i];
	62	fIPCen2[i] += xyz[i]*xyz[i];
	63	}
5b09c01f	64	fIPCenterStat++;
d631e5e7	65	fHVtxXY->Fill(xyz[0],xyz[1]);
5b09c01f	66	return kTRUE;
	67	//
	68	}
	69
	70	//______________________________________________________________________________________
	71	Bool_t AliIntSpotEstimator::ProcessEvent(const AliESDEvent* esd, const AliESDVertex* vtx)
	72	{
	73	// Note: this method will modify the indices of the provided vertex
	74	// Account the vertex in widths estimation
	75	// if vtx is provided its tracks will be used otherwise use the primary vertex of the event
	76	if (!IsValid()) {AliError("Not initialized yet"); return -999;}
	77	if (!fTracks) fTracks = new TObjArray(50);
	78	if (!fVertexer) fVertexer = new AliVertexerTracks();
	79	//
	80	fEvProc++;
	81	int nTracks;
	82	if (!vtx) vtx = (AliESDVertex*) esd->GetPrimaryVertex();
	83	if (!vtx \|\| (nTracks=vtx->GetNIndices())<GetMinTracks()+1) return kFALSE;
	84	//
	85	UShort_t trackID = (UShort_t)vtx->GetIndices();
	86	for (int itr=0;itr<nTracks;itr++) fTracks->Add(esd->GetTrack(trackID[itr]));
	87	//
	88	fVertexer->SetFieldkG( esd->GetMagneticField() );
	89	return ProcessTracks();
	90	//
	91	}
	92
	93	//______________________________________________________________________________________
	94	Bool_t AliIntSpotEstimator::ProcessEvent(const TObjArray* tracks)
	95	{
	96	// This method allows to process the vertex made of arbitrary tracks
	97	// Account the vertex in widths estimation with externally provided tracks
	98	// The indices of the vertex must correspond to array of tracks
	99	//
	100	if (!IsValid()) {AliError("Not initialized yet"); return -999;}
	101	if (!fVertexer) fVertexer = new AliVertexerTracks();
	102	//
	103	fEvProc++;
	104	int nTracks = tracks->GetEntriesFast();
	105	if ( nTracks<GetMinTracks()+1 ) return kFALSE;
	106	for (int itr=0;itr<nTracks;itr++) fTracks->Add(tracks->At(itr));
	107	//
	108	return ProcessTracks();
	109	}
	110
	111	//______________________________________________________________________________________
	112	Bool_t AliIntSpotEstimator::ProcessTracks()
	113	{
	114	// account the vertex made of provided tracks
	115	//
	116	int nTracks = fTracks->GetEntriesFast();
	117	UShort_t *trackID = new UShort_t[nTracks];
	118	for (int i=nTracks;i--;) trackID[i] = i;
	119	double xyzDCA[3],ddca[2],covdca[3];
	120	int nTracks1 = nTracks - 1;
	121	AliExternalTrackParam selTrack,movTrack=0;
	122	UShort_t selTrackID,movTrackID=0;
	123	//
	124	AliESDVertex* recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
d631e5e7	125	if (!recNewVtx \|\| ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
5b09c01f	126	if (recNewVtx) delete recNewVtx;
	127	fTracks->Clear();
	128	delete[] trackID;
	129	return kFALSE;
	130	}
d631e5e7	131	if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
5b09c01f	132	//
	133	double pmn = GetTrackMinP();
	134	double pmx = GetTrackMaxP();
	135	double fieldVal = fVertexer->GetFieldkG();
	136	for (int itr=0;itr<nTracks;itr++) {
	137	selTrack = (AliExternalTrackParam) (fTracks)[itr]; // track to probe
	138	double pTrack = selTrack->GetP();
	139	if (!IsZero(fieldVal) && (pTrack<pmn \|\| pTrack>pmx)) continue;
	140	selTrackID = trackID[itr];
	141	//
	142	if (itr<nTracks1) {
	143	movTrack = (AliExternalTrackParam) (fTracks)[nTracks1]; // save the track
	144	movTrackID = trackID[nTracks1];
	145	(*fTracks)[itr] = movTrack;
	146	trackID[itr] = movTrackID;
	147	}
	148	fTracks->RemoveAt(nTracks1); // move the last track in the place of the probed one
	149	//
	150	// refit the vertex w/o the probed track
	151	recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
	152	if (recNewVtx) {
	153	double told = selTrack->GetX(); // store the original track position
	154	selTrack->PropagateToDCA(recNewVtx,fieldVal,1e4,ddca,covdca); // in principle, done in the vertexer
	155	selTrack->GetXYZ(xyzDCA);
	156	//
	157	double phiTrack = selTrack->Phi();
	158	double cs = TMath::Cos(phiTrack);
	159	double sn = TMath::Sin(phiTrack);
d631e5e7	160	double trDCA = (xyzDCA[0]-fIPCenIni[0]) sn - (xyzDCA[1]-fIPCenIni[1]) cs; // track signed DCA to origin
d631e5e7	161	double vtDCA = (recNewVtx->GetXv()-fIPCenIni[0])sn - (recNewVtx->GetYv()-fIPCenIni[1])cs; // vertex signed DCA to origin
5b09c01f	162	UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
5b09c01f	163	selTrack->PropagateTo(told,fieldVal); // restore the track
d631e5e7	164	if (fNtuple) {
	165	static float ntf[8];
	166	ntf[0] = float(nTracks1);
	167	ntf[1] = recNewVtx->GetXv();
	168	ntf[2] = recNewVtx->GetYv();
	169	ntf[3] = recNewVtx->GetZv();
	170	ntf[4] = xyzDCA[0];
	171	ntf[5] = xyzDCA[1];
	172	ntf[6] = phiTrack;
	173	ntf[7] = pTrack;
	174	fNtuple->Fill(ntf);
	175	}
5b09c01f	176	}
	177	delete recNewVtx;
	178	// restore the track indices
	179	(*fTracks)[itr] = selTrack;
	180	trackID[itr] = selTrackID;
	181	if (itr<nTracks1) {
	182	(*fTracks)[nTracks1] = movTrack;
	183	trackID[nTracks1] = movTrackID;
	184	}
	185	}
	186	//
	187	fTracks->Clear();
	188	delete[] trackID;
	189	return kTRUE;
	190	//
	191	}
	192
	193	//______________________________________________________________________________________
	194	void AliIntSpotEstimator::UpdateEstimators(double rvD, double rtD, double nTracks, double pTrack, double phiTrack)
	195	{
	196	// update the estimator values
	197	double estIP = rvD*rtD;
	198	double estVtx = rvD*(rvD - rtD);
	199	double estTrc = rtD*(rtD - rvD);
	200	//
d631e5e7	201	if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
5b09c01f	202	fEstimVtx->Fill(nTracks, estVtx);
	203	if (pTrack<1e-6) pTrack = GetTrackMinP()+1e6;
	204	fEstimTrc->Fill(1./pTrack,estTrc);
	205	//
	206	}
	207
	208	//______________________________________________________________________________________
	209	void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t estmin,Double_t estmax,
	210	Int_t ntrBins,Int_t ntMn,Int_t ntMx,
d631e5e7	211	Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
5b09c01f	212	{
	213	Clear();
	214	// regularize binning/limits for DCA resolution
	215	nPBins = nPBins<1 ? 1: nPBins;
	216	pmn = pmn>0.1 ? pmn : 0.1;
	217	pmx = pmx>pmn ? pmx : pmn+0.1;
	218	//
	219	// regularize binning/limits for vertex resolution
	220	ntMn = ntMn>2 ? ntMn : 2;
	221	ntMx = ntMx>ntMn ? ntMx : ntMn;
	222	ntrBins = ntrBins<1 ? 1:ntrBins;
	223	int step = (ntMx-ntMn)/ntrBins;
	224	if (step<1) step = 1;
	225	ntrBins = (ntMx-ntMn)/step;
	226	//
	227	// regularize binning for IP sigma vs angle
	228	nPhiBins = nPhiBins>1 ? nPhiBins : 1;
	229	//
	230	// regularize binning/limits for estimators
	231	nestb = nestb>300 ? nestb:300;
	232	estmin = estmin<-2.e-2 ? estmin : -2.e-2;
	233	estmax = estmax> 4.e-2 ? estmax : 4.e-2;
	234	//
	235	TString nm;
	236	nm = GetName();
	237	nm += "diamondEst";
	238	fEstimIP = new TH2F(nm.Data(),nm.Data(),nPhiBins,0.,2.*TMath::Pi(),nestb,estmin,estmax);
	239	//
	240	nm = GetName();
	241	nm += "VResEst";
	242	fEstimVtx = new TH2F(nm.Data(),nm.Data(),ntrBins,ntMn,ntMx, nestb,estmin,estmax);
	243	//
	244	nm = GetName();
	245	nm += "dcaEst";
	246	fEstimTrc = new TH2F(nm.Data(),nm.Data(),nPBins,1./pmx,1./pmn, nestb,estmin,estmax);
d631e5e7	247	//
	248	nm = GetName();
	249	nm += "VtxXY";
	250	fHVtxXY = new TH2F(nm.Data(),nm.Data(),200, -1,1, 200,-1,1);
	251	//
	252	if (ntuple) {
	253	nm = GetName();
	254	nm += "ntuple";
	255	fNtuple = new TNtuple(nm.Data(),nm.Data(),"ntrack:xv:yv:zv:xt:yt:phi:p");
	256	}
5b09c01f	257	//
	258	fVertexer = new AliVertexerTracks(); // the field is set dynamically
	259	fVertexer->SetConstraintOff();
	260	fTracks = new TObjArray(50);
	261	//
	262
	263	}
	264
	265	//______________________________________________________________________________________
d631e5e7	266	Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
5b09c01f	267	{
	268	// get estimate for the IP sigma
	269	if (!IsValid()) {AliError("Not initialized yet"); return -999;}
d631e5e7	270	double cxe,cye;
	271	double cx = GetIPCenter(0,&cxe) - GetIPCenIni(0);
	272	double cy = GetIPCenter(1,&cye) - GetIPCenIni(1);
5b09c01f	273	TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
d631e5e7	274	double merr = 0;
	275	double est = CalcMean(proj, fOutlierCut, &merr) - (cxcx + cycy)/2.;
	276	if (est>0) {
	277	est = TMath::Sqrt(est);
	278	if (err) {
	279	*err = 0;
	280	err = merrmerr;
	281	err += cxcxcxecxe + cycycye*cye;
	282	err = TMath::Sqrt(err)/est/2.;
	283	}
	284	}
	285	else {
	286	est = 0;
	287	if (err) *err = 0;
	288	}
5b09c01f	289	delete proj;
d631e5e7	290	return est;
5b09c01f	291	}
	292
	293	//______________________________________________________________________________________
d631e5e7	294	Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
5b09c01f	295	{
	296	// get estimate for the IP sigma^2
	297	if (!IsValid()) {AliError("Not initialized yet"); return -999;}
	298	int bin = fEstimVtx->GetXaxis()->FindBin(ntr);
	299	if (bin<1 \|\| bin>GetNTrackBins()) {
	300	AliError(Form("Requested vertex multiplicity %d out of defined %d-%d range",
	301	ntr,GetMinTracks(),GetMaxTracks()));
	302	return -1;
	303	}
	304	TH1* proj = fEstimVtx->ProjectionY("vrProj",bin,bin,"e");
d631e5e7	305	double est = CalcMean(proj, fOutlierCut, err);
5b09c01f	306	delete proj;
d631e5e7	307	if (est>0) {
	308	est = TMath::Sqrt(est);
	309	if (err) err /= 2est;
	310	}
	311	else {
	312	est = 0;
	313	if (err) *err = 0;
	314	}
	315	return est;
5b09c01f	316	//
	317	}
	318
	319	//______________________________________________________________________________________
d631e5e7	320	Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
5b09c01f	321	{
	322	// get estimate for the IP sigma^2
	323	if (!IsValid()) {AliError("Not initialized yet"); return -999;}
	324	if (pt<1e-6) pt = GetTrackMinP()+1e6;
	325	pt = 1./pt;
	326	int bin = fEstimTrc->GetXaxis()->FindBin(pt);
	327	if (bin<1 \|\| bin>GetNPBins()) {
	328	AliError(Form("Requested track P %.2f out of defined %.2f-%.2f range",1/pt,GetTrackMinP(),GetTrackMaxP()));
	329	return -1;
	330	}
	331	TH1* proj = fEstimTrc->ProjectionY("trProj",bin,bin,"e");
d631e5e7	332	double est = CalcMean(proj, fOutlierCut, err);
5b09c01f	333	delete proj;
d631e5e7	334	if (est>0) {
	335	est = TMath::Sqrt(est);
	336	if (err) err /= 2est;
	337	}
	338	else {
	339	est = 0;
	340	if (err) *err = 0;
	341	}
	342	return est;
5b09c01f	343	//
	344	}
	345
	346	//______________________________________________________________________________________
d631e5e7	347	Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
5b09c01f	348	{
	349	// calculate mean applying the cut on the outliers
	350	//
	351	double max = histo->GetMaximum();
d631e5e7	352	double cut = (ctfact>0&&ctfact<1.) ? maxctfact : 0;//TMath::Max(1.0,maxctfact) : 0;
5b09c01f	353	int nb = histo->GetNbinsX();
d631e5e7	354	double mean = 0.,cumul = 0, rms = 0;
5b09c01f	355	for (int i=1;i<=nb;i++) {
	356	double vl = histo->GetBinContent(i) - cut;
	357	if (vl<1e-10) continue;
d631e5e7	358	double x = histo->GetBinCenter(i);
	359	mean += vl*x;
	360	rms += vlxx;
5b09c01f	361	cumul += vl;
	362	}
	363	//
d631e5e7	364	mean = cumul>0 ? mean/cumul : 0;
	365	rms -= meanmeancumul;
	366	if (err) {
	367	*err = cumul > 1 ? rms/(cumul-1) : 0;
	368	if (err>0) err = TMath::Sqrt(*err/cumul);
	369	}
5b09c01f	370	//
d631e5e7	371	return mean;
5b09c01f	372	}
	373
	374	//______________________________________________________________________________________
	375	void AliIntSpotEstimator::Print(Option_t *) const
	376	{
	377	if (!IsValid()) {printf("Not initialized yet\n"); return;}
	378	//
d631e5e7	379	double cx,cy,cz,cxe,cye,cze;
	380	cx = GetIPCenter(0,&cxe);
	381	cy = GetIPCenter(1,&cye);
	382	cz = GetIPCenter(2,&cze);
5b09c01f	383	printf("Processed %d events\n",fEvProc);
d631e5e7	384	printf("Estimator for IP center: %+.4e+-%.3e \| %+.4e+-%.3e \| %+.4e+-%.3e\n",
	385	cx,cxe,cy,cye,cz,cze);
	386	printf("Initial IP center was : %+.4e %+.4e %+.4e\n",
	387	GetIPCenIni(0),GetIPCenIni(1),GetIPCenIni(2));
	388	double sgIP = GetIPSigma(0,&cxe);
	389	printf("Estimator for IP sigma : %.4e+-%.3e\n",sgIP,cxe);
5b09c01f	390	//
	391	printf("Estimators for vertex resolution vs Ntracks:\n");
	392
	393	for (int i=1;i<=GetNTrackBins();i++) {
d631e5e7	394	double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
5b09c01f	395	if (IsZero(sig)) continue;
	396	int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
	397	int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
d631e5e7	398	printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
5b09c01f	399	}
	400	//
	401	printf("Estimators for track DCA resolution vs P:\n");
	402	for (int i=1;i<=GetNPBins();i++) {
d631e5e7	403	double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
5b09c01f	404	if (IsZero(sig)) continue;
	405	double pmax = 1./fEstimTrc->GetXaxis()->GetBinLowEdge(i);
	406	double pmin = 1./(fEstimTrc->GetXaxis()->GetBinLowEdge(i)+fEstimTrc->GetXaxis()->GetBinWidth(i));
d631e5e7	407	printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
5b09c01f	408	}
	409	}
	410
	411	//______________________________________________________________________________________
	412	void AliIntSpotEstimator::Clear(Option_t *)
	413	{
	414	// reset all
	415	fEvProc = 0;
	416	fIPCenterStat = 0;
d631e5e7	417	for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
5b09c01f	418	if (fEstimIP) fEstimIP->Reset();
	419	if (fEstimVtx) fEstimVtx->Reset();
	420	if (fEstimTrc) fEstimTrc->Reset();
	421	}
	422
	423	//_____________________________________________________________________
	424	AliIntSpotEstimator &AliIntSpotEstimator::operator += (const AliIntSpotEstimator &src)
	425	{
	426	fEvProc += src.fEvProc;
	427	fIPCenterStat += src.fIPCenterStat;
	428	for (int i=3;i--;) fIPCenter[i] += src.fIPCenter[i];
	429	if (fEstimIP && src.fEstimIP ) fEstimIP->Add(src.fEstimIP);
	430	if (fEstimVtx && src.fEstimVtx) fEstimVtx->Add(src.fEstimVtx);
	431	if (fEstimTrc && src.fEstimTrc) fEstimTrc->Add(src.fEstimTrc);
d631e5e7	432	if (fHVtxXY && src.fHVtxXY) fHVtxXY->Add(src.fHVtxXY);
5b09c01f	433	//
	434	return *this;
	435	}
	436
	437	//_____________________________________________________________________
	438	TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
	439	{
	440	// prepare canvas with report
	441	TCanvas *cnv = new TCanvas(GetName(), GetName(),5,5,700,1000);
	442	gStyle->SetOptStat(0);
	443	gStyle->SetTitleH(0.07);
	444	gStyle->SetTitleW(0.7);
	445	gStyle->SetTitleY(1);
	446	gStyle->SetTitleX(0.2);
	447	//
e93e44d0	448	const Int_t nc=200;
e93e44d0	449	char buff[nc];
5b09c01f	450	//
5b09c01f	451	cnv->Divide(2,2);
	452	cnv->cd(1);
	453	//
	454	TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
e93e44d0	455	snprintf(buff,nc,"%s \| Outliers Cut : %.2e",GetName(),fOutlierCut);
5b09c01f	456	pt->AddText(buff);
5b09c01f	457	//
e93e44d0	458	snprintf(buff,nc,"Processed Events:\n%d",fEvProc);
5b09c01f	459	pt->AddText(buff);
5b09c01f	460	//
e93e44d0	461	snprintf(buff,nc,"Accepted Events\n%d",fIPCenterStat);
5b09c01f	462	pt->AddText(buff);
5b09c01f	463	//
d631e5e7	464	double cx,cy,cz,cxe,cye,cze;
	465	cx = GetIPCenter(0,&cxe);
	466	cy = GetIPCenter(1,&cye);
	467	cz = GetIPCenter(2,&cze);
	468	//
e93e44d0	469	snprintf(buff,nc,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
d631e5e7	470	int(cx1e4),int(cxe1e4),int(cy1e4),int(cye1e4),int(cz1e4),int(cze1e4));
5b09c01f	471	pt->AddText(buff);
5b09c01f	472	//
d631e5e7	473	cx = GetIPSigma(0,&cxe);
e93e44d0	474	snprintf(buff,nc,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx1e4),int(cxe1e4));
5b09c01f	475	pt->AddText(buff);
	476	pt->Draw();
	477	gPad->Modified();
	478	//
	479	cnv->cd(2);
	480	gPad->SetLeftMargin(0.2);
	481	TH1* iph = fEstimIP->ProjectionX();
	482	iph->Reset();
	483	iph->SetTitle("Int.Spot size vs #phi");
d631e5e7	484	for (int i=1;i<=iph->GetNbinsX();i++) {
	485	cx = GetIPSigma(i,&cxe);
	486	iph->SetBinContent(i,cx*1e4);
	487	iph->SetBinError(i,cxe*1e4);
	488	}
5b09c01f	489	iph->GetXaxis()->SetTitle("#phi");
	490	iph->GetYaxis()->SetTitle("#sigma_{IP} [#mum]");
	491	iph->SetMarkerStyle(20);
	492	iph->Draw("p");
	493	gPad->Modified();
	494	iph->SetTitleOffset(2.5,"Y");
	495	//
	496	cnv->cd(3);
	497	gPad->SetLeftMargin(0.2);
	498	TH1* vrh = fEstimVtx->ProjectionX();
	499	vrh->Reset();
	500	vrh->SetTitle("Vertex resolution vs N tracks");
d631e5e7	501	for (int i=1;i<=vrh->GetNbinsX();i++) {
	502	cx = GetVtxSigma( TMath::Nint(vrh->GetBinCenter(i)), &cxe);
	503	vrh->SetBinContent(i,cx*1e4);
	504	vrh->SetBinError(i,cxe*1e4);
	505	}
5b09c01f	506	vrh->GetXaxis()->SetTitle("n tracks");
	507	vrh->GetYaxis()->SetTitle("#sigma_{VTX} [#mum]");
	508	vrh->SetMarkerStyle(20);
	509	vrh->Draw("p");
	510	gPad->Modified();
	511	vrh->SetTitleOffset(2.5,"Y");
	512	//
	513	cnv->cd(4);
	514	gPad->SetLeftMargin(0.2);
	515	TH1* trh = fEstimTrc->ProjectionX();
	516	trh->Reset();
	517	trh->SetTitle("Track DCA resolution vs 1/P");
d631e5e7	518	for (int i=1;i<=trh->GetNbinsX();i++) {
	519	cx = GetDCASigma(1./trh->GetBinCenter(i), &cxe);
	520	trh->SetBinContent(i,cx*1e4);
	521	trh->SetBinError(i,cxe*1e4);
	522	}
5b09c01f	523	trh->GetXaxis()->SetTitle("1/p [GeV]");
	524	trh->GetYaxis()->SetTitle("#sigma_{DCA} [#mum]");
	525	trh->SetMarkerStyle(20);
	526	trh->Draw("p");
	527	gPad->Modified();
	528	trh->SetTitleOffset(2.5,"Y");
	529	//
	530	cnv->cd();
	531	gPad->Modified();
	532	//
	533	if (outname) cnv->Print(outname);
	534	//
	535	return cnv;
	536	}
	537
	538
	539	//_____________________________________________________________________
	540	Long64_t AliIntSpotEstimator::Merge(TCollection *coll)
	541	{
	542	// Merge estimators: used to put together the results of parallel processing
	543	if(!coll) return 0;
	544	if (coll->IsEmpty()) return 1;
	545	TIterator* iter = coll->MakeIterator();
	546	TObject* obj;
	547	int count = 0;
	548	while ((obj = iter->Next())) {
	549	AliIntSpotEstimator* entry = dynamic_cast<AliIntSpotEstimator*>(obj);
e93e44d0	550	if (!entry) continue;
5b09c01f	551	(this) += entry;
	552	count++;
	553	}
	554	return count;
	555	//
	556	}
	557
d631e5e7	558	//_____________________________________________________________________
	559	Double_t AliIntSpotEstimator::GetIPCenter(Int_t id,Double_t *err) const
	560	{
	561	// calculate IP center in axis id and error
	562	double cen = fIPCenterStat>0 ? fIPCenter[id]/fIPCenterStat:0;
	563	if (err) {
	564	err = fIPCenterStat>1 ? (fIPCen2[id] - cencen*fIPCenterStat)/(fIPCenterStat-1) : 0;
	565	err = err > 0 ? TMath::Sqrt(*err/fIPCenterStat) : 0;
	566	}
	567	return cen;
	568	}