fix coverity
[u/mrichter/AliRoot.git] / 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 }