]>
Commit | Line | Data |
---|---|---|
1 | #include "AliRunLoader.h" | |
2 | #include "AliITSDetTypeRec.h" | |
3 | #include "AliITSVertexerIons.h" | |
4 | #include "AliITSVertexerZ.h" | |
5 | #include "AliESDVertex.h" | |
6 | #include "AliITSgeomTGeo.h" | |
7 | #include "AliITSLoader.h" | |
8 | #include "AliITSRecPoint.h" | |
9 | #include "AliLog.h" | |
10 | #include <TTree.h> | |
11 | #include <TH1.h> | |
12 | #include <TF1.h> | |
13 | ||
14 | ////////////////////////////////////////////////////////////////////// | |
15 | // AliITSVertexerIons is a class for full 3D primary vertex // | |
16 | // finding optimized for Ion-Ion interactions // | |
17 | // // | |
18 | // // | |
19 | // // | |
20 | // // | |
21 | // Written by Giuseppe Lo Re and Francesco Riggi // | |
22 | // Giuseppe.Lore@ct.infn.it // | |
23 | // Franco.Riggi@ct.infn.it // | |
24 | // // | |
25 | // Release date: Mar 2004 // | |
26 | // // | |
27 | // // | |
28 | ////////////////////////////////////////////////////////////////////// | |
29 | ||
30 | ||
31 | ClassImp(AliITSVertexerIons) | |
32 | ||
33 | //______________________________________________________________________ | |
34 | AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(), | |
35 | fNpThreshold(0), | |
36 | fMaxDeltaPhi(0), | |
37 | fMaxDeltaZ(0){ | |
38 | // Default Constructor | |
39 | ||
40 | //fITS = 0; | |
41 | SetNpThreshold(); | |
42 | SetMaxDeltaPhi(); | |
43 | SetMaxDeltaZ(); | |
44 | } | |
45 | ||
46 | //______________________________________________________________________ | |
47 | AliITSVertexerIons::~AliITSVertexerIons() { | |
48 | // Default Destructor | |
49 | //fITS = 0; | |
50 | } | |
51 | ||
52 | //______________________________________________________________________ | |
53 | AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(TTree *itsClusterTree){ | |
54 | // Defines the AliESDVertex for the current event | |
55 | ||
56 | fCurrentVertex = 0; | |
57 | ||
58 | AliITSDetTypeRec detTypeRec; | |
59 | ||
60 | detTypeRec.SetTreeAddressR(itsClusterTree); | |
61 | ||
62 | TClonesArray *recpoints = detTypeRec.RecPoints(); | |
63 | AliITSRecPoint *pnt; | |
64 | ||
65 | ||
66 | Int_t npoints=0; | |
67 | Int_t nopoints1=40000; | |
68 | Int_t nopoints2=40000; | |
69 | Float_t p[3]; | |
70 | ||
71 | Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2; | |
72 | z1=new Double_t[nopoints1]; | |
73 | z2=new Double_t[nopoints2]; | |
74 | y1=new Double_t[nopoints1]; | |
75 | y2=new Double_t[nopoints2]; | |
76 | x1=new Double_t[nopoints1]; | |
77 | x2=new Double_t[nopoints2]; | |
78 | phi1=new Double_t[nopoints1]; | |
79 | phi2=new Double_t[nopoints2]; | |
80 | r1=new Double_t[nopoints1]; | |
81 | r2=new Double_t[nopoints2]; | |
82 | ||
83 | Double_t mxpiu = 0; | |
84 | Double_t mxmeno = 0; | |
85 | Double_t mypiu = 0; | |
86 | Double_t mymeno = 0; | |
87 | Double_t r=0; | |
88 | ||
89 | Int_t np1=0, np2=0; | |
90 | for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) { | |
91 | detTypeRec.ResetRecPoints(); | |
92 | itsClusterTree->GetEvent(i); | |
93 | npoints = recpoints->GetEntries(); | |
94 | for (Int_t ipoint=0;ipoint<npoints;ipoint++) { | |
95 | pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint); | |
96 | /* | |
97 | l[0]=pnt->GetDetLocalX(); | |
98 | l[1]=0; | |
99 | l[2]=pnt->GetDetLocalZ(); | |
100 | g2->LtoG(i, l, p); | |
101 | */ | |
102 | pnt->GetGlobalXYZ(p); | |
103 | r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); | |
104 | ||
105 | if(i<80 && TMath::Abs(p[2])<14.35) { | |
106 | y1[np1]=p[1]; | |
107 | x1[np1]=p[0]; | |
108 | z1[np1]=p[2]; | |
109 | if(p[0]>0) mxpiu++; | |
110 | if(p[0]<0) mxmeno++; | |
111 | if(p[1]>0) mypiu++; | |
112 | if(p[1]<0) mymeno++; | |
113 | np1++; | |
114 | } | |
115 | if(i>=80 && TMath::Abs(p[2])<14.35) { | |
116 | y2[np2]=p[1]; | |
117 | x2[np2]=p[0]; | |
118 | z2[np2]=p[2]; | |
119 | np2++; | |
120 | } | |
121 | } | |
122 | } | |
123 | ||
124 | if(np1<fNpThreshold) { | |
125 | Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n"); | |
126 | Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold); | |
127 | AliITSVertexerZ *dovert = new AliITSVertexerZ(); | |
128 | fCurrentVertex =dovert->FindVertexForCurrentEvent(itsClusterTree); | |
129 | delete dovert; | |
130 | return fCurrentVertex; | |
131 | } | |
132 | ||
133 | if(!np1 || !np2) { | |
134 | Error("FindVertexForCurrentEvent","No points in the pixer layers"); | |
135 | return fCurrentVertex; | |
136 | } | |
137 | AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2)); | |
138 | ||
139 | Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno); | |
140 | Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno); | |
141 | ||
142 | Double_t x0 = 2.86*asparx; | |
143 | Double_t y0 = 2.86*aspary; | |
144 | ||
145 | AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0)); | |
146 | ||
147 | for(Int_t i=0;i<np1;i++) { | |
148 | x1[i]-=x0; | |
149 | y1[i]-=y0; | |
150 | PhiFunc(x1[i],y1[i],phi1[i]); | |
151 | r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]); | |
152 | } | |
153 | for(Int_t i=0;i<np2;i++) { | |
154 | x2[i]-=x0; | |
155 | y2[i]-=y0; | |
156 | PhiFunc(x2[i],y2[i],phi2[i]); | |
157 | r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]); | |
158 | } | |
159 | ||
160 | Int_t nbinxy=400; | |
161 | Int_t nbinz=400; | |
162 | TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5); | |
163 | TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5); | |
164 | TH1F *hzv=new TH1F("hzv","",nbinz,-15,15); | |
165 | ||
166 | Double_t dphi; | |
167 | for(Int_t j=0; j<np1; j++) { | |
168 | for(Int_t k=0; k<np2; k++) { | |
169 | dphi=TMath::Abs(phi2[k]-phi1[j]); | |
170 | if(dphi>180) dphi = 360-dphi; | |
171 | if(dphi>fMaxDeltaPhi) continue; | |
172 | hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1)); | |
173 | } | |
174 | } | |
175 | ||
176 | // Fitting ... | |
177 | Double_t max; | |
178 | Int_t binMax; | |
179 | Double_t maxcenter; | |
180 | ||
181 | max = hzv->GetMaximum(); | |
182 | binMax=hzv->GetMaximumBin(); | |
183 | maxcenter=hzv->GetBinCenter(binMax); | |
184 | Double_t dxy=1.5; | |
185 | ||
186 | TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy); | |
187 | fz->SetParameter(0,max); | |
188 | fz->SetParameter(1,maxcenter); | |
189 | fz->SetParameter(2,0.1); | |
190 | fz->SetLineColor(kRed); | |
191 | hzv->Fit("fz","RQ0"); | |
192 | ||
193 | if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) { | |
194 | Warning("FindVertexForCurrentEvent","AliITSVertexerIonsGeneral finder is not reliable for events with x and y vertex coordinates too far from the origin of the reference system. Their values will be set to 0.\n"); | |
195 | Double_t position[3]={0,0,fz->GetParameter(1)}; | |
196 | Double_t resolution[3]={0,0,0}; | |
197 | Double_t snr[3]={0,0,0}; | |
198 | Char_t name[30]; | |
199 | sprintf(name,"Vertex"); | |
200 | fCurrentVertex = new AliESDVertex(position,resolution,snr,name); | |
201 | return fCurrentVertex; | |
202 | } | |
203 | ||
204 | for(Int_t j=0; j<np1; j++) { | |
205 | for(Int_t k=0; k<np2; k++) { | |
206 | if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue; | |
207 | if(y2[k]==y1[j]) continue; | |
208 | hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k])); | |
209 | if(x2[k]==x1[j]) continue; | |
210 | hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k])); | |
211 | } | |
212 | } | |
213 | ||
214 | TF1 *fx = new TF1 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",x0-dxy,x0+dxy); | |
215 | fx->SetParameter(0,100); | |
216 | Double_t dist=0.3; | |
217 | Double_t xapprox=FindMaxAround(x0,hxv,dist); | |
218 | AliDebug(1,Form("xapprox = %f",xapprox)); | |
219 | fx->SetParameter(1,xapprox); | |
220 | Double_t difcentroid=0.07; | |
221 | fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid); | |
222 | fx->SetParameter(2,0.1); | |
223 | fx->SetLineColor(kRed); | |
224 | hxv->Fit("fx","RQW0"); | |
225 | ||
226 | TF1 *fy = new TF1 ("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]+[4]*x+[5]*x*x",y0-dxy,y0+dxy); | |
227 | fy->SetParameter(0,100); | |
228 | Double_t yapprox=FindMaxAround(y0,hyv,dist); | |
229 | AliDebug(1,Form("yapprox = %f",yapprox)); | |
230 | fy->SetParameter(1,yapprox); | |
231 | fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid); | |
232 | fy->SetParameter(2,0.1); | |
233 | fy->SetLineColor(kRed); | |
234 | hyv->Fit("fy","RQW0"); | |
235 | ||
236 | delete [] z1; | |
237 | delete [] z2; | |
238 | delete [] y1; | |
239 | delete [] y2; | |
240 | delete [] x1; | |
241 | delete [] x2; | |
242 | delete [] r1; | |
243 | delete [] r2; | |
244 | delete [] phi1; | |
245 | delete [] phi2; | |
246 | delete hxv; | |
247 | delete hyv; | |
248 | delete hzv; | |
249 | ||
250 | Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)}; | |
251 | Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)}; | |
252 | Double_t snr[3]={0,0,0}; | |
253 | ||
254 | Char_t name[30]; | |
255 | sprintf(name,"Vertex"); | |
256 | fCurrentVertex = new AliESDVertex(position,resolution,snr,name); | |
257 | ||
258 | return fCurrentVertex; | |
259 | } | |
260 | ||
261 | //______________________________________________________________________ | |
262 | void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) { | |
263 | // Method for the computation of the phi angle given x and y. The result is in degrees. | |
264 | if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578); | |
265 | if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180; | |
266 | if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180; | |
267 | if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;; | |
268 | } | |
269 | ||
270 | //________________________________________________________ | |
271 | void AliITSVertexerIons::PrintStatus() const { | |
272 | // Print current status | |
273 | cout <<"=======================================================\n"; | |
274 | if(fCurrentVertex)fCurrentVertex->PrintStatus(); | |
275 | } | |
276 | ||
277 | //________________________________________________________ | |
278 | Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) { | |
279 | // It finds the maximum of h within point-distance and distance+point. | |
280 | Int_t maxcontent=0; | |
281 | Int_t maxbin=0; | |
282 | for(Int_t i=0;i<h->GetNbinsX();i++) { | |
283 | Int_t content=(Int_t)h->GetBinContent(i); | |
284 | Double_t center=(Double_t)h->GetBinCenter(i); | |
285 | if(TMath::Abs(center-point)>distance) continue; | |
286 | if(content>maxcontent) {maxcontent=content;maxbin=i;} | |
287 | } | |
288 | Double_t max=h->GetBinCenter(maxbin); | |
289 | return max; | |
290 | } | |
291 | ||
292 |