Coverity fixes
[u/mrichter/AliRoot.git] / PWG1 / AliIntSpotEstimator.cxx
CommitLineData
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
16ClassImp(AliIntSpotEstimator)
17
18//______________________________________________________________________________________
d631e5e7 19AliIntSpotEstimator::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,""),
d631e5e7 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)
5b09c01f 27{
d631e5e7 28 InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
5b09c01f 29}
30
31//______________________________________________________________________________________
32AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef)
33 : TNamed("IPEstimator",""),
d631e5e7 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)
5b09c01f 36{
37 if (initDef) InitEstimators();
38 //
39}
40
41//______________________________________________________________________________________
42AliIntSpotEstimator::~AliIntSpotEstimator()
43{
44 delete fEstimIP;
45 delete fEstimVtx;
46 delete fEstimTrc;
47 delete fVertexer;
d631e5e7 48 delete fNtuple;
5b09c01f 49}
50
51//______________________________________________________________________________________
52Bool_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];
d631e5e7 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 }
5b09c01f 64 fIPCenterStat++;
d631e5e7 65 fHVtxXY->Fill(xyz[0],xyz[1]);
5b09c01f 66 return kTRUE;
67 //
68}
69
70//______________________________________________________________________________________
71Bool_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//______________________________________________________________________________________
94Bool_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//______________________________________________________________________________________
112Bool_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);
d631e5e7 125 if (!recNewVtx || ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
5b09c01f 126 if (recNewVtx) delete recNewVtx;
127 fTracks->Clear();
128 delete[] trackID;
129 return kFALSE;
130 }
d631e5e7 131 if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
5b09c01f 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);
d631e5e7 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
5b09c01f 162 UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
163 selTrack->PropagateTo(told,fieldVal); // restore the track
d631e5e7 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 }
5b09c01f 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//______________________________________________________________________________________
194void 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 //
d631e5e7 201 if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
5b09c01f 202 fEstimVtx->Fill(nTracks, estVtx);
203 if (pTrack<1e-6) pTrack = GetTrackMinP()+1e6;
204 fEstimTrc->Fill(1./pTrack,estTrc);
205 //
206}
207
208//______________________________________________________________________________________
209void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t estmin,Double_t estmax,
210 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
d631e5e7 211 Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
5b09c01f 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);
d631e5e7 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 }
5b09c01f 257 //
258 fVertexer = new AliVertexerTracks(); // the field is set dynamically
259 fVertexer->SetConstraintOff();
260 fTracks = new TObjArray(50);
261 //
262
263}
264
265//______________________________________________________________________________________
d631e5e7 266Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
5b09c01f 267{
268 // get estimate for the IP sigma
269 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
d631e5e7 270 double cxe,cye;
271 double cx = GetIPCenter(0,&cxe) - GetIPCenIni(0);
272 double cy = GetIPCenter(1,&cye) - GetIPCenIni(1);
5b09c01f 273 TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
d631e5e7 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 }
5b09c01f 289 delete proj;
d631e5e7 290 return est;
5b09c01f 291}
292
293//______________________________________________________________________________________
d631e5e7 294Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
5b09c01f 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");
d631e5e7 305 double est = CalcMean(proj, fOutlierCut, err);
5b09c01f 306 delete proj;
d631e5e7 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;
5b09c01f 316 //
317}
318
319//______________________________________________________________________________________
d631e5e7 320Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
5b09c01f 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");
d631e5e7 332 double est = CalcMean(proj, fOutlierCut, err);
5b09c01f 333 delete proj;
d631e5e7 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;
5b09c01f 343 //
344}
345
346//______________________________________________________________________________________
d631e5e7 347Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
5b09c01f 348{
349 // calculate mean applying the cut on the outliers
350 //
351 double max = histo->GetMaximum();
d631e5e7 352 double cut = (ctfact>0&&ctfact<1.) ? max*ctfact : 0;//TMath::Max(1.0,max*ctfact) : 0;
5b09c01f 353 int nb = histo->GetNbinsX();
d631e5e7 354 double mean = 0.,cumul = 0, rms = 0;
5b09c01f 355 for (int i=1;i<=nb;i++) {
356 double vl = histo->GetBinContent(i) - cut;
357 if (vl<1e-10) continue;
d631e5e7 358 double x = histo->GetBinCenter(i);
359 mean += vl*x;
360 rms += vl*x*x;
5b09c01f 361 cumul += vl;
362 }
363 //
d631e5e7 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 }
5b09c01f 370 //
d631e5e7 371 return mean;
5b09c01f 372}
373
374//______________________________________________________________________________________
375void AliIntSpotEstimator::Print(Option_t *) const
376{
377 if (!IsValid()) {printf("Not initialized yet\n"); return;}
378 //
d631e5e7 379 double cx,cy,cz,cxe,cye,cze;
380 cx = GetIPCenter(0,&cxe);
381 cy = GetIPCenter(1,&cye);
382 cz = GetIPCenter(2,&cze);
5b09c01f 383 printf("Processed %d events\n",fEvProc);
d631e5e7 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);
5b09c01f 390 //
391 printf("Estimators for vertex resolution vs Ntracks:\n");
392
393 for (int i=1;i<=GetNTrackBins();i++) {
d631e5e7 394 double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
5b09c01f 395 if (IsZero(sig)) continue;
396 int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
397 int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
d631e5e7 398 printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
5b09c01f 399 }
400 //
401 printf("Estimators for track DCA resolution vs P:\n");
402 for (int i=1;i<=GetNPBins();i++) {
d631e5e7 403 double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
5b09c01f 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));
d631e5e7 407 printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
5b09c01f 408 }
409}
410
411//______________________________________________________________________________________
412void AliIntSpotEstimator::Clear(Option_t *)
413{
414 // reset all
415 fEvProc = 0;
416 fIPCenterStat = 0;
d631e5e7 417 for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
5b09c01f 418 if (fEstimIP) fEstimIP->Reset();
419 if (fEstimVtx) fEstimVtx->Reset();
420 if (fEstimTrc) fEstimTrc->Reset();
421}
422
423//_____________________________________________________________________
424AliIntSpotEstimator &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);
d631e5e7 432 if (fHVtxXY && src.fHVtxXY) fHVtxXY->Add(src.fHVtxXY);
5b09c01f 433 //
434 return *this;
435}
436
437//_____________________________________________________________________
438TCanvas* 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 //
e93e44d0 448 const Int_t nc=200;
449 char buff[nc];
5b09c01f 450 //
5b09c01f 451 cnv->Divide(2,2);
452 cnv->cd(1);
453 //
454 TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
e93e44d0 455 snprintf(buff,nc,"%s | Outliers Cut : %.2e",GetName(),fOutlierCut);
5b09c01f 456 pt->AddText(buff);
457 //
e93e44d0 458 snprintf(buff,nc,"Processed Events:\n%d",fEvProc);
5b09c01f 459 pt->AddText(buff);
460 //
e93e44d0 461 snprintf(buff,nc,"Accepted Events\n%d",fIPCenterStat);
5b09c01f 462 pt->AddText(buff);
463 //
d631e5e7 464 double cx,cy,cz,cxe,cye,cze;
465 cx = GetIPCenter(0,&cxe);
466 cy = GetIPCenter(1,&cye);
467 cz = GetIPCenter(2,&cze);
468 //
e93e44d0 469 snprintf(buff,nc,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
d631e5e7 470 int(cx*1e4),int(cxe*1e4),int(cy*1e4),int(cye*1e4),int(cz*1e4),int(cze*1e4));
5b09c01f 471 pt->AddText(buff);
472 //
d631e5e7 473 cx = GetIPSigma(0,&cxe);
e93e44d0 474 snprintf(buff,nc,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx*1e4),int(cxe*1e4));
5b09c01f 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");
d631e5e7 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 }
5b09c01f 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");
d631e5e7 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 }
5b09c01f 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");
d631e5e7 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 }
5b09c01f 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//_____________________________________________________________________
540Long64_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);
e93e44d0 550 if (!entry) continue;
5b09c01f 551 (*this) += *entry;
552 count++;
553 }
554 return count;
555 //
556}
557
d631e5e7 558//_____________________________________________________________________
559Double_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}