]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliIntSpotEstimator.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adding the estimator for the interaction spot size (R. Shahoyan)
[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>
9#include <TObjArray.h>
10#include <TCanvas.h>
11#include <TPaveText.h>
12#include <TStyle.h>
13
14
15ClassImp(AliIntSpotEstimator)
16
17//______________________________________________________________________________________
18AliIntSpotEstimator::AliIntSpotEstimator(const char* name,Double_t outcut,
19 Int_t nPhiBins,Int_t nestb,
20 Double_t estmin,Double_t estmax,
21 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
22 Int_t nPBins,Double_t pmn,Double_t pmx)
23: TNamed(name,""),
24 fEvProc(0),fIPCenterStat(0),fOutlierCut(outcut),fEstimIP(0),fEstimVtx(0),fEstimTrc(0),
25 fVertexer(0),fTracks(0)
26{
27 InitEstimators(nPhiBins,nestb,estmin,estmax,ntrBins,ntMn,ntMx,nPBins,pmn,pmx);
28}
29
30//______________________________________________________________________________________
31AliIntSpotEstimator::AliIntSpotEstimator(Bool_t initDef)
32 : TNamed("IPEstimator",""),
33 fEvProc(0),fIPCenterStat(0),fOutlierCut(1e-4),
34 fEstimIP(0),fEstimVtx(0),fEstimTrc(0),fVertexer(0),fTracks(0)
35{
36 if (initDef) InitEstimators();
37 //
38}
39
40//______________________________________________________________________________________
41AliIntSpotEstimator::~AliIntSpotEstimator()
42{
43 delete fEstimIP;
44 delete fEstimVtx;
45 delete fEstimTrc;
46 delete fVertexer;
47}
48
49//______________________________________________________________________________________
50Bool_t AliIntSpotEstimator::ProcessIPCenter(const AliESDVertex* vtx)
51{
52 // account the vertex in IP center estimation
53 if (!vtx) return kFALSE;
54 double xyz[3];
55 vtx->GetXYZ(xyz);
56 double r = xyz[0]*xyz[0] + xyz[1]*xyz[1];
57 if (r>1.) return kFALSE;
58 for (int i=3;i--;) fIPCenter[i] += xyz[i];
59 fIPCenterStat++;
60 return kTRUE;
61 //
62}
63
64//______________________________________________________________________________________
65Bool_t AliIntSpotEstimator::ProcessEvent(const AliESDEvent* esd, const AliESDVertex* vtx)
66{
67 // Note: this method will modify the indices of the provided vertex
68 // Account the vertex in widths estimation
69 // if vtx is provided its tracks will be used otherwise use the primary vertex of the event
70 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
71 if (!fTracks) fTracks = new TObjArray(50);
72 if (!fVertexer) fVertexer = new AliVertexerTracks();
73 //
74 fEvProc++;
75 int nTracks;
76 if (!vtx) vtx = (AliESDVertex*) esd->GetPrimaryVertex();
77 if (!vtx || (nTracks=vtx->GetNIndices())<GetMinTracks()+1) return kFALSE;
78 //
79 UShort_t *trackID = (UShort_t*)vtx->GetIndices();
80 for (int itr=0;itr<nTracks;itr++) fTracks->Add(esd->GetTrack(trackID[itr]));
81 //
82 fVertexer->SetFieldkG( esd->GetMagneticField() );
83 return ProcessTracks();
84 //
85}
86
87//______________________________________________________________________________________
88Bool_t AliIntSpotEstimator::ProcessEvent(const TObjArray* tracks)
89{
90 // This method allows to process the vertex made of arbitrary tracks
91 // Account the vertex in widths estimation with externally provided tracks
92 // The indices of the vertex must correspond to array of tracks
93 //
94 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
95 if (!fVertexer) fVertexer = new AliVertexerTracks();
96 //
97 fEvProc++;
98 int nTracks = tracks->GetEntriesFast();
99 if ( nTracks<GetMinTracks()+1 ) return kFALSE;
100 for (int itr=0;itr<nTracks;itr++) fTracks->Add(tracks->At(itr));
101 //
102 return ProcessTracks();
103}
104
105//______________________________________________________________________________________
106Bool_t AliIntSpotEstimator::ProcessTracks()
107{
108 // account the vertex made of provided tracks
109 //
110 int nTracks = fTracks->GetEntriesFast();
111 UShort_t *trackID = new UShort_t[nTracks];
112 for (int i=nTracks;i--;) trackID[i] = i;
113 double xyzDCA[3],ddca[2],covdca[3];
114 int nTracks1 = nTracks - 1;
115 AliExternalTrackParam *selTrack,*movTrack=0;
116 UShort_t selTrackID,movTrackID=0;
117 //
118 AliESDVertex* recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
119 if (!recNewVtx ||
120 ((nTracks=recNewVtx->GetNIndices())<GetMinTracks()) || // in case some tracks were rejected
121 !ProcessIPCenter(recNewVtx)) {
122 if (recNewVtx) delete recNewVtx;
123 fTracks->Clear();
124 delete[] trackID;
125 return kFALSE;
126 }
127 //
128 double pmn = GetTrackMinP();
129 double pmx = GetTrackMaxP();
130 double fieldVal = fVertexer->GetFieldkG();
131 for (int itr=0;itr<nTracks;itr++) {
132 selTrack = (AliExternalTrackParam*) (*fTracks)[itr]; // track to probe
133 double pTrack = selTrack->GetP();
134 if (!IsZero(fieldVal) && (pTrack<pmn || pTrack>pmx)) continue;
135 selTrackID = trackID[itr];
136 //
137 if (itr<nTracks1) {
138 movTrack = (AliExternalTrackParam*) (*fTracks)[nTracks1]; // save the track
139 movTrackID = trackID[nTracks1];
140 (*fTracks)[itr] = movTrack;
141 trackID[itr] = movTrackID;
142 }
143 fTracks->RemoveAt(nTracks1); // move the last track in the place of the probed one
144 //
145 // refit the vertex w/o the probed track
146 recNewVtx = fVertexer->VertexForSelectedTracks(fTracks,trackID,kTRUE,kFALSE,kFALSE);
147 if (recNewVtx) {
148 double told = selTrack->GetX(); // store the original track position
149 selTrack->PropagateToDCA(recNewVtx,fieldVal,1e4,ddca,covdca); // in principle, done in the vertexer
150 selTrack->GetXYZ(xyzDCA);
151 //
152 double phiTrack = selTrack->Phi();
153 double cs = TMath::Cos(phiTrack);
154 double sn = TMath::Sin(phiTrack);
155 double trDCA = xyzDCA[0] *sn - xyzDCA[1] *cs; // track signed DCA to origin
156 double vtDCA = recNewVtx->GetXv()*sn - recNewVtx->GetYv()*cs; // vertex signed DCA to origin
157 UpdateEstimators(vtDCA,trDCA, nTracks1, pTrack, phiTrack);
158 selTrack->PropagateTo(told,fieldVal); // restore the track
159 }
160 delete recNewVtx;
161 // restore the track indices
162 (*fTracks)[itr] = selTrack;
163 trackID[itr] = selTrackID;
164 if (itr<nTracks1) {
165 (*fTracks)[nTracks1] = movTrack;
166 trackID[nTracks1] = movTrackID;
167 }
168 }
169 //
170 fTracks->Clear();
171 delete[] trackID;
172 return kTRUE;
173 //
174}
175
176//______________________________________________________________________________________
177void AliIntSpotEstimator::UpdateEstimators(double rvD, double rtD, double nTracks, double pTrack, double phiTrack)
178{
179 // update the estimator values
180 double estIP = rvD*rtD;
181 double estVtx = rvD*(rvD - rtD);
182 double estTrc = rtD*(rtD - rvD);
183 //
184 fEstimIP->Fill(phiTrack, estIP);
185 fEstimVtx->Fill(nTracks, estVtx);
186 if (pTrack<1e-6) pTrack = GetTrackMinP()+1e6;
187 fEstimTrc->Fill(1./pTrack,estTrc);
188 //
189}
190
191//______________________________________________________________________________________
192void AliIntSpotEstimator::InitEstimators(Int_t nPhiBins,Int_t nestb,Double_t estmin,Double_t estmax,
193 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
194 Int_t nPBins,Double_t pmn,Double_t pmx)
195{
196 Clear();
197 // regularize binning/limits for DCA resolution
198 nPBins = nPBins<1 ? 1: nPBins;
199 pmn = pmn>0.1 ? pmn : 0.1;
200 pmx = pmx>pmn ? pmx : pmn+0.1;
201 //
202 // regularize binning/limits for vertex resolution
203 ntMn = ntMn>2 ? ntMn : 2;
204 ntMx = ntMx>ntMn ? ntMx : ntMn;
205 ntrBins = ntrBins<1 ? 1:ntrBins;
206 int step = (ntMx-ntMn)/ntrBins;
207 if (step<1) step = 1;
208 ntrBins = (ntMx-ntMn)/step;
209 //
210 // regularize binning for IP sigma vs angle
211 nPhiBins = nPhiBins>1 ? nPhiBins : 1;
212 //
213 // regularize binning/limits for estimators
214 nestb = nestb>300 ? nestb:300;
215 estmin = estmin<-2.e-2 ? estmin : -2.e-2;
216 estmax = estmax> 4.e-2 ? estmax : 4.e-2;
217 //
218 TString nm;
219 nm = GetName();
220 nm += "diamondEst";
221 fEstimIP = new TH2F(nm.Data(),nm.Data(),nPhiBins,0.,2.*TMath::Pi(),nestb,estmin,estmax);
222 //
223 nm = GetName();
224 nm += "VResEst";
225 fEstimVtx = new TH2F(nm.Data(),nm.Data(),ntrBins,ntMn,ntMx, nestb,estmin,estmax);
226 //
227 nm = GetName();
228 nm += "dcaEst";
229 fEstimTrc = new TH2F(nm.Data(),nm.Data(),nPBins,1./pmx,1./pmn, nestb,estmin,estmax);
230 //
231 fVertexer = new AliVertexerTracks(); // the field is set dynamically
232 fVertexer->SetConstraintOff();
233 fTracks = new TObjArray(50);
234 //
235
236}
237
238//______________________________________________________________________________________
239Double_t AliIntSpotEstimator::GetIPSigma(Int_t bin) const
240{
241 // get estimate for the IP sigma
242 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
243 double cx = GetIPCenter(0);
244 double cy = GetIPCenter(1);
245 TH1* proj = fEstimIP->ProjectionY("ipProj",bin<1 ? 1:bin, bin<1 ? GetNPhiBins():bin,"e");
246 double est = CalcMean(proj, fOutlierCut) - (cx*cx + cy*cy)/2.;
247 delete proj;
248 return est>0 ? TMath::Sqrt(est) : 0;
249}
250
251//______________________________________________________________________________________
252Double_t AliIntSpotEstimator::GetVtxSigma(int ntr) const
253{
254 // get estimate for the IP sigma^2
255 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
256 int bin = fEstimVtx->GetXaxis()->FindBin(ntr);
257 if (bin<1 || bin>GetNTrackBins()) {
258 AliError(Form("Requested vertex multiplicity %d out of defined %d-%d range",
259 ntr,GetMinTracks(),GetMaxTracks()));
260 return -1;
261 }
262 TH1* proj = fEstimVtx->ProjectionY("vrProj",bin,bin,"e");
263 double est = CalcMean(proj, fOutlierCut);
264 delete proj;
265 return est>0 ? TMath::Sqrt(est) : 0;
266 //
267}
268
269//______________________________________________________________________________________
270Double_t AliIntSpotEstimator::GetDCASigma(double pt) const
271{
272 // get estimate for the IP sigma^2
273 if (!IsValid()) {AliError("Not initialized yet"); return -999;}
274 if (pt<1e-6) pt = GetTrackMinP()+1e6;
275 pt = 1./pt;
276 int bin = fEstimTrc->GetXaxis()->FindBin(pt);
277 if (bin<1 || bin>GetNPBins()) {
278 AliError(Form("Requested track P %.2f out of defined %.2f-%.2f range",1/pt,GetTrackMinP(),GetTrackMaxP()));
279 return -1;
280 }
281 TH1* proj = fEstimTrc->ProjectionY("trProj",bin,bin,"e");
282 double est = CalcMean(proj, fOutlierCut);
283 delete proj;
284 return est>0 ? TMath::Sqrt(est) : 0;
285 //
286}
287
288//______________________________________________________________________________________
289Double_t AliIntSpotEstimator::CalcMean(TH1* histo, Double_t ctfact)
290{
291 // calculate mean applying the cut on the outliers
292 //
293 double max = histo->GetMaximum();
294 double cut = (ctfact>0&&ctfact<1.) ? TMath::Max(1.0,max*ctfact) : 0;
295 int nb = histo->GetNbinsX();
296 double mean = 0.,cumul = 0;
297 for (int i=1;i<=nb;i++) {
298 double vl = histo->GetBinContent(i) - cut;
299 if (vl<1e-10) continue;
300 mean += vl*histo->GetBinCenter(i);
301 cumul += vl;
302 }
303 //
304 return cumul>0 ? mean/cumul : 0;
305 //
306}
307
308//______________________________________________________________________________________
309void AliIntSpotEstimator::Print(Option_t *) const
310{
311 if (!IsValid()) {printf("Not initialized yet\n"); return;}
312 //
313 printf("Processed %d events\n",fEvProc);
314 printf("Estimator for IP center: %+.4e %+.4e %+.4e\n",
315 GetIPCenter(0),GetIPCenter(1),GetIPCenter(2));
316 double sgIP = GetIPSigma();
317 printf("Estimator for IP sigma : %.4e\n",sgIP);
318 //
319 printf("Estimators for vertex resolution vs Ntracks:\n");
320
321 for (int i=1;i<=GetNTrackBins();i++) {
322 double sig = GetVtxSigma( (int)fEstimVtx->GetXaxis()->GetBinCenter(i) );
323 if (IsZero(sig)) continue;
324 int tmin = TMath::Nint(fEstimVtx->GetXaxis()->GetBinLowEdge(i));
325 int tmax = tmin + int(fEstimVtx->GetXaxis()->GetBinWidth(i));
326 printf("%3d-%3d : %.4e\n",tmin,tmax,sig);
327 }
328 //
329 printf("Estimators for track DCA resolution vs P:\n");
330 for (int i=1;i<=GetNPBins();i++) {
331 double sig = GetDCASigma( 1./fEstimTrc->GetXaxis()->GetBinCenter(i) );
332 if (IsZero(sig)) continue;
333 double pmax = 1./fEstimTrc->GetXaxis()->GetBinLowEdge(i);
334 double pmin = 1./(fEstimTrc->GetXaxis()->GetBinLowEdge(i)+fEstimTrc->GetXaxis()->GetBinWidth(i));
335 printf("%.2f-%.2f : %.4e\n",pmin,pmax,sig);
336 }
337}
338
339//______________________________________________________________________________________
340void AliIntSpotEstimator::Clear(Option_t *)
341{
342 // reset all
343 fEvProc = 0;
344 fIPCenterStat = 0;
345 fIPCenter[0] = fIPCenter[1] = fIPCenter[2] = 0;
346 if (fEstimIP) fEstimIP->Reset();
347 if (fEstimVtx) fEstimVtx->Reset();
348 if (fEstimTrc) fEstimTrc->Reset();
349}
350
351//_____________________________________________________________________
352AliIntSpotEstimator &AliIntSpotEstimator::operator += (const AliIntSpotEstimator &src)
353{
354 fEvProc += src.fEvProc;
355 fIPCenterStat += src.fIPCenterStat;
356 for (int i=3;i--;) fIPCenter[i] += src.fIPCenter[i];
357 if (fEstimIP && src.fEstimIP ) fEstimIP->Add(src.fEstimIP);
358 if (fEstimVtx && src.fEstimVtx) fEstimVtx->Add(src.fEstimVtx);
359 if (fEstimTrc && src.fEstimTrc) fEstimTrc->Add(src.fEstimTrc);
360 //
361 return *this;
362}
363
364//_____________________________________________________________________
365TCanvas* AliIntSpotEstimator::CreateReport(const char* outname)
366{
367 // prepare canvas with report
368 TCanvas *cnv = new TCanvas(GetName(), GetName(),5,5,700,1000);
369 gStyle->SetOptStat(0);
370 gStyle->SetTitleH(0.07);
371 gStyle->SetTitleW(0.7);
372 gStyle->SetTitleY(1);
373 gStyle->SetTitleX(0.2);
374 //
375 char buff[200];
376 //
377 double sigIPtot = GetIPSigma()*1e4;
378 //
379 cnv->Divide(2,2);
380 cnv->cd(1);
381 //
382 TPaveText *pt = new TPaveText(0.05,0.05,0.95,0.95,"blNDC");
383 sprintf(buff,"%s | Outliers Cut : %.2e",GetName(),fOutlierCut);
384 pt->AddText(buff);
385 //
386 sprintf(buff,"Processed Events:\n%d",fEvProc);
387 pt->AddText(buff);
388 //
389 sprintf(buff,"Accepted Events\n%d",fIPCenterStat);
390 pt->AddText(buff);
391 //
392 sprintf(buff,"Int.Spot X:Y:Z (cm)\n%+.4e : %+.4e : %+.4e",GetIPCenter(0),GetIPCenter(1),GetIPCenter(2));
393 pt->AddText(buff);
394 //
395 sprintf(buff,"Int.Spot #sigma (#mum):\n%d",int(sigIPtot));
396 pt->AddText(buff);
397 pt->Draw();
398 gPad->Modified();
399 //
400 cnv->cd(2);
401 gPad->SetLeftMargin(0.2);
402 TH1* iph = fEstimIP->ProjectionX();
403 iph->Reset();
404 iph->SetTitle("Int.Spot size vs #phi");
405 for (int i=1;i<=iph->GetNbinsX();i++) iph->SetBinContent(i,GetIPSigma(i)*1e4);
406 iph->GetXaxis()->SetTitle("#phi");
407 iph->GetYaxis()->SetTitle("#sigma_{IP} [#mum]");
408 iph->SetMarkerStyle(20);
409 iph->Draw("p");
410 gPad->Modified();
411 iph->SetTitleOffset(2.5,"Y");
412 //
413 cnv->cd(3);
414 gPad->SetLeftMargin(0.2);
415 TH1* vrh = fEstimVtx->ProjectionX();
416 vrh->Reset();
417 vrh->SetTitle("Vertex resolution vs N tracks");
418 for (int i=1;i<=vrh->GetNbinsX();i++) vrh->SetBinContent(i,GetVtxSigma(TMath::Nint(vrh->GetBinCenter(i)))*1e4);
419 vrh->GetXaxis()->SetTitle("n tracks");
420 vrh->GetYaxis()->SetTitle("#sigma_{VTX} [#mum]");
421 vrh->SetMarkerStyle(20);
422 vrh->Draw("p");
423 gPad->Modified();
424 vrh->SetTitleOffset(2.5,"Y");
425 //
426 cnv->cd(4);
427 gPad->SetLeftMargin(0.2);
428 TH1* trh = fEstimTrc->ProjectionX();
429 trh->Reset();
430 trh->SetTitle("Track DCA resolution vs 1/P");
431 for (int i=1;i<=trh->GetNbinsX();i++) trh->SetBinContent(i,GetDCASigma(1./trh->GetBinCenter(i))*1e4);
432 trh->GetXaxis()->SetTitle("1/p [GeV]");
433 trh->GetYaxis()->SetTitle("#sigma_{DCA} [#mum]");
434 trh->SetMarkerStyle(20);
435 trh->Draw("p");
436 gPad->Modified();
437 trh->SetTitleOffset(2.5,"Y");
438 //
439 cnv->cd();
440 gPad->Modified();
441 //
442 if (outname) cnv->Print(outname);
443 //
444 return cnv;
445}
446
447
448//_____________________________________________________________________
449Long64_t AliIntSpotEstimator::Merge(TCollection *coll)
450{
451 // Merge estimators: used to put together the results of parallel processing
452 if(!coll) return 0;
453 if (coll->IsEmpty()) return 1;
454 TIterator* iter = coll->MakeIterator();
455 TObject* obj;
456 int count = 0;
457 while ((obj = iter->Next())) {
458 AliIntSpotEstimator* entry = dynamic_cast<AliIntSpotEstimator*>(obj);
459 (*this) += *entry;
460 count++;
461 }
462 return count;
463 //
464}
465