PWGPP/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),
  26   fIPCenIni(),
  27   fIPCenter(),
  28   fIPCen2(),
  29   fEstimIP(0),
  30   fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
  31 {
  32   InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
  33 }
  34
  35 //______________________________________________________________________________________
  36 AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef)
  37   : TNamed("IPEstimator",""),
  38     fEvProc(0),fIPCenterStat(0),fMinTracksForIP(2),fOutlierCut(1e-4),
  39     fIPCenIni(),
  40     fIPCenter(),
  41     fIPCen2(),
  42     fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
  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;
  55   delete fNtuple;
  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];
  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   }
  71   fIPCenterStat++;
  72   fHVtxXY->Fill(xyz[0],xyz[1]);
  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);
 132   if (!recNewVtx || ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
 133     if (recNewVtx) delete recNewVtx;
 134     fTracks->Clear();
 135     delete[] trackID;
 136     return kFALSE;
 137   }
 138   if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
 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);
 167       double trDCA = (xyzDCA[0]-fIPCenIni[0])         *sn - (xyzDCA[1]-fIPCenIni[1])         *cs;  // track signed DCA to origin
 168       double vtDCA = (recNewVtx->GetX()-fIPCenIni[0])*sn - (recNewVtx->GetY()-fIPCenIni[1])*cs;  // vertex signed DCA to origin
 169       UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
 170       selTrack->PropagateTo(told,fieldVal);    // restore the track
 171       if (fNtuple) {
 172         static float ntf[8];
 173         ntf[0] = float(nTracks1);
 174         ntf[1] = recNewVtx->GetX();
 175         ntf[2] = recNewVtx->GetY();
 176         ntf[3] = recNewVtx->GetZ();
 177         ntf[4] = xyzDCA[0];
 178         ntf[5] = xyzDCA[1];
 179         ntf[6] = phiTrack;
 180         ntf[7] = pTrack;
 181         fNtuple->Fill(ntf);
 182       }
 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   //
 208   if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
 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,
 218                                          Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
 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);
 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   }
 264   //
 265   fVertexer = new AliVertexerTracks(); // the field is set dynamically
 266   fVertexer->SetConstraintOff();
 267   fTracks = new TObjArray(50);
 268   //
 269
 270 }
 271
 272 //______________________________________________________________________________________
 273 Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
 274 {
 275   // get estimate for the IP sigma
 276   if (!IsValid()) {AliError("Not initialized yet"); return -999;}
 277   double cxe,cye;
 278   double cx  = GetIPCenter(0,&cxe) - GetIPCenIni(0);
 279   double cy  = GetIPCenter(1,&cye) - GetIPCenIni(1);
 280   TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
 281   double merr = 0;
 282   double est = CalcMean(proj, fOutlierCut, &merr) - (cx*cx + cy*cy)/2.;
 283   if (est>0) {
 284     est = TMath::Sqrt(est);
 285     if (err) {
 286       *err = 0;
 287       *err = merr*merr;
 288       *err += cx*cx*cxe*cxe + cy*cy*cye*cye;
 289       *err = TMath::Sqrt(*err)/est/2.;
 290     }
 291   }
 292   else {
 293     est = 0;
 294     if (err) *err = 0;
 295   }
 296   delete proj;
 297   return est;
 298 }
 299
 300 //______________________________________________________________________________________
 301 Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
 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");
 312   double est = CalcMean(proj, fOutlierCut, err);
 313   delete proj;
 314   if (est>0) {
 315     est = TMath::Sqrt(est);
 316     if (err) *err /= 2*est;
 317   }
 318   else {
 319     est = 0;
 320     if (err) *err = 0;
 321   }
 322   return est;
 323   //
 324 }
 325
 326 //______________________________________________________________________________________
 327 Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
 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");
 339   double est = CalcMean(proj, fOutlierCut, err);
 340   delete proj;
 341    if (est>0) {
 342     est = TMath::Sqrt(est);
 343     if (err) *err /= 2*est;
 344   }
 345   else {
 346     est = 0;
 347     if (err) *err = 0;
 348   }
 349   return est;
 350   //
 351 }
 352
 353 //______________________________________________________________________________________
 354 Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
 355 {
 356   // calculate mean applying the cut on the outliers
 357   //
 358   double max = histo->GetMaximum();
 359   double cut = (ctfact>0&&ctfact<1.) ? max*ctfact : 0;//TMath::Max(1.0,max*ctfact) : 0;
 360   int nb = histo->GetNbinsX();
 361   double mean = 0.,cumul = 0, rms = 0;
 362   for (int i=1;i<=nb;i++) {
 363     double vl = histo->GetBinContent(i) - cut;
 364     if (vl<1e-10) continue;
 365     double x = histo->GetBinCenter(i);
 366     mean  += vl*x;
 367     rms   += vl*x*x;
 368     cumul += vl;
 369   }
 370   //
 371   mean = cumul>0 ? mean/cumul : 0;
 372   rms -= mean*mean*cumul;
 373   if (err) {
 374     *err = cumul > 1 ? rms/(cumul-1) : 0;
 375     if (*err>0) *err = TMath::Sqrt(*err/cumul);
 376   }
 377   //
 378   return mean;
 379 }
 380
 381 //______________________________________________________________________________________
 382 void AliIntSpotEstimator::Print(Option_t *) const
 383 {
 384   if (!IsValid()) {printf("Not initialized yet\n"); return;}
 385   //
 386   double cx,cy,cz,cxe,cye,cze;
 387   cx = GetIPCenter(0,&cxe);
 388   cy = GetIPCenter(1,&cye);
 389   cz = GetIPCenter(2,&cze);
 390   printf("Processed %d events\n",fEvProc);
 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);
 397   //
 398   printf("Estimators for vertex resolution vs Ntracks:\n");
 399
 400   for (int i=1;i<=GetNTrackBins();i++) {
 401     double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
 402     if (IsZero(sig)) continue;
 403     int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
 404     int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
 405     printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
 406   }
 407   //
 408   printf("Estimators for track DCA resolution vs P:\n");
 409   for (int i=1;i<=GetNPBins();i++) {
 410     double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
 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));
 414     printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
 415   }
 416 }
 417
 418 //______________________________________________________________________________________
 419 void AliIntSpotEstimator::Clear(Option_t *)
 420 {
 421   // reset all
 422   fEvProc = 0;
 423   fIPCenterStat = 0;
 424   for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
 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);
 439   if (fHVtxXY   && src.fHVtxXY)   fHVtxXY->Add(src.fHVtxXY);
 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   //
 455   const Int_t nc=200;
 456   char buff[nc];
 457   //
 458   cnv->Divide(2,2);
 459   cnv->cd(1);
 460   //
 461   TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
 462   snprintf(buff,nc,"%s | Outliers Cut : %.2e",GetName(),fOutlierCut);
 463   pt->AddText(buff);
 464   //
 465   snprintf(buff,nc,"Processed Events:\n%d",fEvProc);
 466   pt->AddText(buff);
 467   //
 468   snprintf(buff,nc,"Accepted  Events\n%d",fIPCenterStat);
 469   pt->AddText(buff);
 470   //
 471   double cx,cy,cz,cxe,cye,cze;
 472   cx = GetIPCenter(0,&cxe);
 473   cy = GetIPCenter(1,&cye);
 474   cz = GetIPCenter(2,&cze);
 475   //
 476   snprintf(buff,nc,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
 477           int(cx*1e4),int(cxe*1e4),int(cy*1e4),int(cye*1e4),int(cz*1e4),int(cze*1e4));
 478   pt->AddText(buff);
 479   //
 480   cx = GetIPSigma(0,&cxe);
 481   snprintf(buff,nc,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx*1e4),int(cxe*1e4));
 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");
 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   }
 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");
 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   }
 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");
 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   }
 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);
 557     if (!entry) continue;
 558     (*this) += *entry;
 559     count++;
 560   }
 561   return count;
 562   //
 563 }
 564
 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] - cen*cen*fIPCenterStat)/(fIPCenterStat-1) : 0;
 572     *err = *err > 0 ? TMath::Sqrt(*err/fIPCenterStat) : 0;
 573   }
 574   return cen;
 575 }