PWG1/AliIntSpotEstimator.cxx

   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>
   9 #include <TNtuple.h>
  10 #include <TObjArray.h>
  11 #include <TCanvas.h>
  12 #include <TPaveText.h>
  13 #include <TStyle.h>
  14
  15
  16 ClassImp(AliIntSpotEstimator)
  17
  18 //______________________________________________________________________________________
  19 AliIntSpotEstimator::AliIntSpotEstimator(const char* name,Double_t outcut,Int_t ntrIP,
  20                                          Int_t nPhiBins,Int_t nestb,
  21                                          Double_t estmin,Double_t estmax,
  22                                          Int_t ntrBins,Int_t ntMn,Int_t ntMx,
  23                                          Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
  24 : TNamed(name,""),
  25   fEvProc(0),fIPCenterStat(0),fMinTracksForIP(ntrIP>2?ntrIP:2),fOutlierCut(outcut),fEstimIP(0),
  26   fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
  27 {
  28   InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
  29 }
  30
  31 //______________________________________________________________________________________
  32 AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef)
  33   : TNamed("IPEstimator",""),
  34     fEvProc(0),fIPCenterStat(0),fMinTracksForIP(2),fOutlierCut(1e-4),
  35     fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
  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;
  48   delete fNtuple;
  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];
  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   }
  64   fIPCenterStat++;
  65   fHVtxXY->Fill(xyz[0],xyz[1]);
  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);
 125   if (!recNewVtx || ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
 126     if (recNewVtx) delete recNewVtx;
 127     fTracks->Clear();
 128     delete[] trackID;
 129     return kFALSE;
 130   }
 131   if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
 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);
 160       double trDCA = (xyzDCA[0]-fIPCenIni[0])         *sn - (xyzDCA[1]-fIPCenIni[1])         *cs;  // track signed DCA to origin
 161       double vtDCA = (recNewVtx->GetXv()-fIPCenIni[0])*sn - (recNewVtx->GetYv()-fIPCenIni[1])*cs;  // vertex signed DCA to origin
 162       UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
 163       selTrack->PropagateTo(told,fieldVal);    // restore the track
 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       }
 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   //
 201   if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
 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,
 211                                          Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
 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);
 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   }
 257   //
 258   fVertexer = new AliVertexerTracks(); // the field is set dynamically
 259   fVertexer->SetConstraintOff();
 260   fTracks = new TObjArray(50);
 261   //
 262
 263 }
 264
 265 //______________________________________________________________________________________
 266 Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
 267 {
 268   // get estimate for the IP sigma
 269   if (!IsValid()) {AliError("Not initialized yet"); return -999;}
 270   double cxe,cye;
 271   double cx  = GetIPCenter(0,&cxe) - GetIPCenIni(0);
 272   double cy  = GetIPCenter(1,&cye) - GetIPCenIni(1);
 273   TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
 274   double merr = 0;
 275   double est = CalcMean(proj, fOutlierCut, &merr) - (cx*cx + cy*cy)/2.;
 276   if (est>0) {
 277     est = TMath::Sqrt(est);
 278     if (err) {
 279       *err = 0;
 280       *err = merr*merr;
 281       *err += cx*cx*cxe*cxe + cy*cy*cye*cye;
 282       *err = TMath::Sqrt(*err)/est/2.;
 283     }
 284   }
 285   else {
 286     est = 0;
 287     if (err) *err = 0;
 288   }
 289   delete proj;
 290   return est;
 291 }
 292
 293 //______________________________________________________________________________________
 294 Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
 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");
 305   double est = CalcMean(proj, fOutlierCut, err);
 306   delete proj;
 307   if (est>0) {
 308     est = TMath::Sqrt(est);
 309     if (err) *err /= 2*est;
 310   }
 311   else {
 312     est = 0;
 313     if (err) *err = 0;
 314   }
 315   return est;
 316   //
 317 }
 318
 319 //______________________________________________________________________________________
 320 Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
 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");
 332   double est = CalcMean(proj, fOutlierCut, err);
 333   delete proj;
 334    if (est>0) {
 335     est = TMath::Sqrt(est);
 336     if (err) *err /= 2*est;
 337   }
 338   else {
 339     est = 0;
 340     if (err) *err = 0;
 341   }
 342   return est;
 343   //
 344 }
 345
 346 //______________________________________________________________________________________
 347 Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
 348 {
 349   // calculate mean applying the cut on the outliers
 350   //
 351   double max = histo->GetMaximum();
 352   double cut = (ctfact>0&&ctfact<1.) ? max*ctfact : 0;//TMath::Max(1.0,max*ctfact) : 0;
 353   int nb = histo->GetNbinsX();
 354   double mean = 0.,cumul = 0, rms = 0;
 355   for (int i=1;i<=nb;i++) {
 356     double vl = histo->GetBinContent(i) - cut;
 357     if (vl<1e-10) continue;
 358     double x = histo->GetBinCenter(i);
 359     mean  += vl*x;
 360     rms   += vl*x*x;
 361     cumul += vl;
 362   }
 363   //
 364   mean = cumul>0 ? mean/cumul : 0;
 365   rms -= mean*mean*cumul;
 366   if (err) {
 367     *err = cumul > 1 ? rms/(cumul-1) : 0;
 368     if (*err>0) *err = TMath::Sqrt(*err/cumul);
 369   }
 370   //
 371   return mean;
 372 }
 373
 374 //______________________________________________________________________________________
 375 void AliIntSpotEstimator::Print(Option_t *) const
 376 {
 377   if (!IsValid()) {printf("Not initialized yet\n"); return;}
 378   //
 379   double cx,cy,cz,cxe,cye,cze;
 380   cx = GetIPCenter(0,&cxe);
 381   cy = GetIPCenter(1,&cye);
 382   cz = GetIPCenter(2,&cze);
 383   printf("Processed %d events\n",fEvProc);
 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);
 390   //
 391   printf("Estimators for vertex resolution vs Ntracks:\n");
 392
 393   for (int i=1;i<=GetNTrackBins();i++) {
 394     double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
 395     if (IsZero(sig)) continue;
 396     int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
 397     int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
 398     printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
 399   }
 400   //
 401   printf("Estimators for track DCA resolution vs P:\n");
 402   for (int i=1;i<=GetNPBins();i++) {
 403     double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
 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));
 407     printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
 408   }
 409 }
 410
 411 //______________________________________________________________________________________
 412 void AliIntSpotEstimator::Clear(Option_t *)
 413 {
 414   // reset all
 415   fEvProc = 0;
 416   fIPCenterStat = 0;
 417   for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
 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);
 432   if (fHVtxXY   && src.fHVtxXY)   fHVtxXY->Add(src.fHVtxXY);
 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   //
 448   const Int_t nc=200;
 449   char buff[nc];
 450   //
 451   cnv->Divide(2,2);
 452   cnv->cd(1);
 453   //
 454   TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
 455   snprintf(buff,nc,"%s | Outliers Cut : %.2e",GetName(),fOutlierCut);
 456   pt->AddText(buff);
 457   //
 458   snprintf(buff,nc,"Processed Events:\n%d",fEvProc);
 459   pt->AddText(buff);
 460   //
 461   snprintf(buff,nc,"Accepted  Events\n%d",fIPCenterStat);
 462   pt->AddText(buff);
 463   //
 464   double cx,cy,cz,cxe,cye,cze;
 465   cx = GetIPCenter(0,&cxe);
 466   cy = GetIPCenter(1,&cye);
 467   cz = GetIPCenter(2,&cze);
 468   //
 469   snprintf(buff,nc,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
 470           int(cx*1e4),int(cxe*1e4),int(cy*1e4),int(cye*1e4),int(cz*1e4),int(cze*1e4));
 471   pt->AddText(buff);
 472   //
 473   cx = GetIPSigma(0,&cxe);
 474   snprintf(buff,nc,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx*1e4),int(cxe*1e4));
 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");
 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   }
 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");
 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   }
 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");
 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   }
 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);
 550     if (!entry) continue;
 551     (*this) += *entry;
 552     count++;
 553   }
 554   return count;
 555   //
 556 }
 557
 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] - cen*cen*fIPCenterStat)/(fIPCenterStat-1) : 0;
 565     *err = *err > 0 ? TMath::Sqrt(*err/fIPCenterStat) : 0;
 566   }
 567   return cen;
 568 }