better log validation
[u/mrichter/AliRoot.git] / PWGPP / 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,""),
bc676576 25 fEvProc(0),fIPCenterStat(0),fMinTracksForIP(ntrIP>2?ntrIP:2),fOutlierCut(outcut),
26 fIPCenIni(),
27 fIPCenter(),
28 fIPCen2(),
29 fEstimIP(0),
d631e5e7 30 fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
5b09c01f 31{
d631e5e7 32 InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx,ntuple);
5b09c01f 33}
34
35//______________________________________________________________________________________
36AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef)
37 : TNamed("IPEstimator",""),
d631e5e7 38 fEvProc(0),fIPCenterStat(0),fMinTracksForIP(2),fOutlierCut(1e-4),
bc676576 39 fIPCenIni(),
40 fIPCenter(),
41 fIPCen2(),
d631e5e7 42 fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fHVtxXY(0),fNtuple(0),fVertexer(0),fTracks(0)
5b09c01f 43{
44 if (initDef) InitEstimators();
45 //
46}
47
48//______________________________________________________________________________________
49AliIntSpotEstimator::~AliIntSpotEstimator()
50{
51 delete fEstimIP;
52 delete fEstimVtx;
53 delete fEstimTrc;
54 delete fVertexer;
d631e5e7 55 delete fNtuple;
5b09c01f 56}
57
58//______________________________________________________________________________________
59Bool_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];
d631e5e7 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 }
5b09c01f 71 fIPCenterStat++;
d631e5e7 72 fHVtxXY->Fill(xyz[0],xyz[1]);
5b09c01f 73 return kTRUE;
74 //
75}
76
77//______________________________________________________________________________________
78Bool_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//______________________________________________________________________________________
101Bool_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//______________________________________________________________________________________
119Bool_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);
d631e5e7 132 if (!recNewVtx || ((nTracks=recNewVtx->GetNIndices())<GetMinTracks())) {
5b09c01f 133 if (recNewVtx) delete recNewVtx;
134 fTracks->Clear();
135 delete[] trackID;
136 return kFALSE;
137 }
d631e5e7 138 if (nTracks>=fMinTracksForIP) ProcessIPCenter(recNewVtx);
5b09c01f 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);
d631e5e7 167 double trDCA = (xyzDCA[0]-fIPCenIni[0]) *sn - (xyzDCA[1]-fIPCenIni[1]) *cs; // track signed DCA to origin
e690d4d0 168 double vtDCA = (recNewVtx->GetX()-fIPCenIni[0])*sn - (recNewVtx->GetY()-fIPCenIni[1])*cs; // vertex signed DCA to origin
5b09c01f 169 UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
170 selTrack->PropagateTo(told,fieldVal); // restore the track
d631e5e7 171 if (fNtuple) {
172 static float ntf[8];
173 ntf[0] = float(nTracks1);
e690d4d0 174 ntf[1] = recNewVtx->GetX();
175 ntf[2] = recNewVtx->GetY();
176 ntf[3] = recNewVtx->GetZ();
d631e5e7 177 ntf[4] = xyzDCA[0];
178 ntf[5] = xyzDCA[1];
179 ntf[6] = phiTrack;
180 ntf[7] = pTrack;
181 fNtuple->Fill(ntf);
182 }
5b09c01f 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//______________________________________________________________________________________
201void 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 //
d631e5e7 208 if (nTracks >= fMinTracksForIP) fEstimIP->Fill(phiTrack, estIP);
5b09c01f 209 fEstimVtx->Fill(nTracks, estVtx);
210 if (pTrack<1e-6) pTrack = GetTrackMinP()+1e6;
211 fEstimTrc->Fill(1./pTrack,estTrc);
212 //
213}
214
215//______________________________________________________________________________________
216void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t estmin,Double_t estmax,
217 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
d631e5e7 218 Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t ntuple)
5b09c01f 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);
d631e5e7 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 }
5b09c01f 264 //
265 fVertexer = new AliVertexerTracks(); // the field is set dynamically
266 fVertexer->SetConstraintOff();
267 fTracks = new TObjArray(50);
268 //
269
270}
271
272//______________________________________________________________________________________
d631e5e7 273Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin, Double_t *err) const
5b09c01f 274{
275 // get estimate for the IP sigma
276 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
d631e5e7 277 double cxe,cye;
278 double cx = GetIPCenter(0,&cxe) - GetIPCenIni(0);
279 double cy = GetIPCenter(1,&cye) - GetIPCenIni(1);
5b09c01f 280 TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
d631e5e7 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 }
5b09c01f 296 delete proj;
d631e5e7 297 return est;
5b09c01f 298}
299
300//______________________________________________________________________________________
d631e5e7 301Double_t AliIntSpotEstimator::GetVtxSigma(int ntr, double* err) const
5b09c01f 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");
d631e5e7 312 double est = CalcMean(proj, fOutlierCut, err);
5b09c01f 313 delete proj;
d631e5e7 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;
5b09c01f 323 //
324}
325
326//______________________________________________________________________________________
d631e5e7 327Double_t AliIntSpotEstimator::GetDCASigma(double pt, double *err) const
5b09c01f 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");
d631e5e7 339 double est = CalcMean(proj, fOutlierCut, err);
5b09c01f 340 delete proj;
d631e5e7 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;
5b09c01f 350 //
351}
352
353//______________________________________________________________________________________
d631e5e7 354Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact, Double_t *err)
5b09c01f 355{
356 // calculate mean applying the cut on the outliers
357 //
358 double max = histo->GetMaximum();
d631e5e7 359 double cut = (ctfact>0&&ctfact<1.) ? max*ctfact : 0;//TMath::Max(1.0,max*ctfact) : 0;
5b09c01f 360 int nb = histo->GetNbinsX();
d631e5e7 361 double mean = 0.,cumul = 0, rms = 0;
5b09c01f 362 for (int i=1;i<=nb;i++) {
363 double vl = histo->GetBinContent(i) - cut;
364 if (vl<1e-10) continue;
d631e5e7 365 double x = histo->GetBinCenter(i);
366 mean += vl*x;
367 rms += vl*x*x;
5b09c01f 368 cumul += vl;
369 }
370 //
d631e5e7 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 }
5b09c01f 377 //
d631e5e7 378 return mean;
5b09c01f 379}
380
381//______________________________________________________________________________________
382void AliIntSpotEstimator::Print(Option_t *) const
383{
384 if (!IsValid()) {printf("Not initialized yet\n"); return;}
385 //
d631e5e7 386 double cx,cy,cz,cxe,cye,cze;
387 cx = GetIPCenter(0,&cxe);
388 cy = GetIPCenter(1,&cye);
389 cz = GetIPCenter(2,&cze);
5b09c01f 390 printf("Processed %d events\n",fEvProc);
d631e5e7 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);
5b09c01f 397 //
398 printf("Estimators for vertex resolution vs Ntracks:\n");
399
400 for (int i=1;i<=GetNTrackBins();i++) {
d631e5e7 401 double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i), &cxe );
5b09c01f 402 if (IsZero(sig)) continue;
403 int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
404 int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
d631e5e7 405 printf("%3d-%3d : %.4e+-%.3e\n",tmin,tmax,sig,cxe);
5b09c01f 406 }
407 //
408 printf("Estimators for track DCA resolution vs P:\n");
409 for (int i=1;i<=GetNPBins();i++) {
d631e5e7 410 double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i), &cxe );
5b09c01f 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));
d631e5e7 414 printf("%.2f-%.2f : %.4e+-%.3e\n",pmin,pmax,sig, cxe);
5b09c01f 415 }
416}
417
418//______________________________________________________________________________________
419void AliIntSpotEstimator::Clear(Option_t *)
420{
421 // reset all
422 fEvProc = 0;
423 fIPCenterStat = 0;
d631e5e7 424 for (int i=3;i--;) fIPCenter[i] = fIPCenIni[i] = 0.;
5b09c01f 425 if (fEstimIP) fEstimIP->Reset();
426 if (fEstimVtx) fEstimVtx->Reset();
427 if (fEstimTrc) fEstimTrc->Reset();
428}
429
430//_____________________________________________________________________
431AliIntSpotEstimator &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);
d631e5e7 439 if (fHVtxXY && src.fHVtxXY) fHVtxXY->Add(src.fHVtxXY);
5b09c01f 440 //
441 return *this;
442}
443
444//_____________________________________________________________________
445TCanvas* 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 //
e93e44d0 455 const Int_t nc=200;
456 char buff[nc];
5b09c01f 457 //
5b09c01f 458 cnv->Divide(2,2);
459 cnv->cd(1);
460 //
461 TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
e93e44d0 462 snprintf(buff,nc,"%s | Outliers Cut : %.2e",GetName(),fOutlierCut);
5b09c01f 463 pt->AddText(buff);
464 //
e93e44d0 465 snprintf(buff,nc,"Processed Events:\n%d",fEvProc);
5b09c01f 466 pt->AddText(buff);
467 //
e93e44d0 468 snprintf(buff,nc,"Accepted Events\n%d",fIPCenterStat);
5b09c01f 469 pt->AddText(buff);
470 //
d631e5e7 471 double cx,cy,cz,cxe,cye,cze;
472 cx = GetIPCenter(0,&cxe);
473 cy = GetIPCenter(1,&cye);
474 cz = GetIPCenter(2,&cze);
475 //
e93e44d0 476 snprintf(buff,nc,"Int.Spot (#mum)\n%+d#pm%d\n%+d#pm%d\n%+d#pm%d",
d631e5e7 477 int(cx*1e4),int(cxe*1e4),int(cy*1e4),int(cye*1e4),int(cz*1e4),int(cze*1e4));
5b09c01f 478 pt->AddText(buff);
479 //
d631e5e7 480 cx = GetIPSigma(0,&cxe);
e93e44d0 481 snprintf(buff,nc,"Int.Spot #sigma (#mum):\n%d#pm%d",int(cx*1e4),int(cxe*1e4));
5b09c01f 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");
d631e5e7 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 }
5b09c01f 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");
d631e5e7 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 }
5b09c01f 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");
d631e5e7 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 }
5b09c01f 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//_____________________________________________________________________
547Long64_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);
e93e44d0 557 if (!entry) continue;
5b09c01f 558 (*this) += *entry;
559 count++;
560 }
561 return count;
562 //
563}
564
d631e5e7 565//_____________________________________________________________________
566Double_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}