]>
Commit | Line | Data |
---|---|---|
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 | } |