Add INT7 for pA.# Please enter the commit message for your changes. s starting
[u/mrichter/AliRoot.git] / 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->GetXv()-fIPCenIni[0])*sn - (recNewVtx->GetYv()-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->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       }
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 }