]>
Commit | Line | Data |
---|---|---|
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 | ||
16 | ClassImp(AliIntSpotEstimator) | |
17 | ||
18 | //______________________________________________________________________________________ | |
d631e5e7 | 19 | AliIntSpotEstimator::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 | //______________________________________________________________________________________ | |
36 | AliIntSpotEstimator::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 | //______________________________________________________________________________________ | |
49 | AliIntSpotEstimator::~AliIntSpotEstimator() | |
50 | { | |
51 | delete fEstimIP; | |
52 | delete fEstimVtx; | |
53 | delete fEstimTrc; | |
54 | delete fVertexer; | |
d631e5e7 | 55 | delete fNtuple; |
5b09c01f | 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]; | |
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 | //______________________________________________________________________________________ | |
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); | |
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 |
168 | double vtDCA = (recNewVtx->GetXv()-fIPCenIni[0])*sn - (recNewVtx->GetYv()-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); | |
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 | } | |
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 | //______________________________________________________________________________________ | |
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 | // | |
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 | //______________________________________________________________________________________ | |
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, | |
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 | 273 | Double_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 | 301 | Double_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 | 327 | Double_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 | 354 | Double_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 | //______________________________________________________________________________________ | |
382 | void 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 | //______________________________________________________________________________________ | |
419 | void 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 | //_____________________________________________________________________ | |
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); | |
d631e5e7 | 439 | if (fHVtxXY && src.fHVtxXY) fHVtxXY->Add(src.fHVtxXY); |
5b09c01f | 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 | // | |
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 | //_____________________________________________________________________ | |
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); | |
e93e44d0 | 557 | if (!entry) continue; |
5b09c01f | 558 | (*this) += *entry; |
559 | count++; | |
560 | } | |
561 | return count; | |
562 | // | |
563 | } | |
564 | ||
d631e5e7 | 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 | } |