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