]>
Commit | Line | Data |
---|---|---|
6cae184e | 1 | #include "AliRunLoader.h" |
7d62fb64 | 2 | #include "AliITSDetTypeRec.h" |
b5d712e2 | 3 | #include "AliITSVertexerIons.h" |
cd706e57 | 4 | #include "AliITSVertexerZ.h" |
b5d712e2 | 5 | #include "AliESDVertex.h" |
32e63e47 | 6 | #include "AliITSgeomTGeo.h" |
88cb7938 | 7 | #include "AliITSLoader.h" |
c5f0f3c1 | 8 | #include "AliITSRecPoint.h" |
cd706e57 | 9 | #include "AliLog.h" |
b5d712e2 | 10 | #include <TTree.h> |
c5f0f3c1 | 11 | #include <TH1.h> |
12 | #include <TF1.h> | |
b5d712e2 | 13 | |
c5f0f3c1 | 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 // | |
c5f0f3c1 | 23 | // Franco.Riggi@ct.infn.it // |
24 | // // | |
fe4280bf | 25 | // Release date: Mar 2004 // |
c5f0f3c1 | 26 | // // |
27 | // // | |
28 | ////////////////////////////////////////////////////////////////////// | |
29 | ||
c5f0f3c1 | 30 | |
b5d712e2 | 31 | ClassImp(AliITSVertexerIons) |
32 | ||
fe4280bf | 33 | //______________________________________________________________________ |
8221b41b | 34 | AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(), |
35 | fNpThreshold(0), | |
36 | fMaxDeltaPhi(0), | |
37 | fMaxDeltaZ(0){ | |
c5f0f3c1 | 38 | // Default Constructor |
39 | ||
7d62fb64 | 40 | //fITS = 0; |
83028144 | 41 | SetNpThreshold(); |
fe4280bf | 42 | SetMaxDeltaPhi(); |
43 | SetMaxDeltaZ(); | |
c5f0f3c1 | 44 | } |
45 | ||
46 | //______________________________________________________________________ | |
8221b41b | 47 | AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn), |
48 | fNpThreshold(0), | |
49 | fMaxDeltaPhi(0), | |
50 | fMaxDeltaZ(0) { | |
c5f0f3c1 | 51 | // Standard constructor |
88cb7938 | 52 | |
7d62fb64 | 53 | //fITS = 0; |
83028144 | 54 | SetNpThreshold(); |
fe4280bf | 55 | SetMaxDeltaPhi(); |
b5d712e2 | 56 | SetMaxDeltaZ(); |
57 | } | |
8221b41b | 58 | /* |
b5d712e2 | 59 | //______________________________________________________________________ |
60 | AliITSVertexerIons::AliITSVertexerIons(const AliITSVertexerIons &source):AliITSVertexer(source) { | |
61 | // Copy constructor | |
62 | // Copies are not allowed. The method is protected to avoid misuse. | |
63 | Error("AliITSVertexerIons","Copy constructor not allowed\n"); | |
64 | } | |
8221b41b | 65 | */ |
b5d712e2 | 66 | //_________________________________________________________________________ |
8221b41b | 67 | //AliITSVertexerIons& AliITSVertexerIons::operator=(const AliITSVertexerIons &/*source*/) { |
b5d712e2 | 68 | // Assignment operator |
69 | // Assignment is not allowed. The method is protected to avoid misuse. | |
8221b41b | 70 | //Error("= operator","Assignment operator not allowed\n"); |
71 | //return *this; | |
72 | //} | |
c5f0f3c1 | 73 | |
74 | ||
75 | //______________________________________________________________________ | |
76 | AliITSVertexerIons::~AliITSVertexerIons() { | |
77 | // Default Destructor | |
7d62fb64 | 78 | //fITS = 0; |
c5f0f3c1 | 79 | } |
80 | ||
81 | //______________________________________________________________________ | |
b5d712e2 | 82 | AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){ |
83 | // Defines the AliESDVertex for the current event | |
fe4280bf | 84 | |
c5f0f3c1 | 85 | fCurrentVertex = 0; |
fe4280bf | 86 | |
88cb7938 | 87 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
7d62fb64 | 88 | AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader"); |
32e63e47 | 89 | /* |
7d62fb64 | 90 | TDirectory * olddir = gDirectory; |
91 | rl->CdGAFile(); | |
92 | AliITSgeom* g2 = (AliITSgeom*)gDirectory->Get("AliITSgeom"); | |
93 | olddir->cd(); | |
32e63e47 | 94 | */ |
7d62fb64 | 95 | |
fe4280bf | 96 | TTree *tr = itsloader->TreeR(); |
7d62fb64 | 97 | AliITSDetTypeRec detTypeRec; |
98 | ||
99 | detTypeRec.SetTreeAddressR(tr); | |
100 | ||
101 | TClonesArray *recpoints = detTypeRec.RecPoints(); | |
102 | AliITSRecPoint *pnt; | |
103 | ||
88cb7938 | 104 | |
fe4280bf | 105 | Int_t npoints=0; |
106 | Int_t nopoints1=40000; | |
107 | Int_t nopoints2=40000; | |
32e63e47 | 108 | Float_t p[3]; |
c5f0f3c1 | 109 | |
83028144 | 110 | Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2; |
111 | z1=new Double_t[nopoints1]; | |
112 | z2=new Double_t[nopoints2]; | |
113 | y1=new Double_t[nopoints1]; | |
114 | y2=new Double_t[nopoints2]; | |
115 | x1=new Double_t[nopoints1]; | |
116 | x2=new Double_t[nopoints2]; | |
117 | phi1=new Double_t[nopoints1]; | |
118 | phi2=new Double_t[nopoints2]; | |
119 | r1=new Double_t[nopoints1]; | |
120 | r2=new Double_t[nopoints2]; | |
c5f0f3c1 | 121 | |
fe4280bf | 122 | Double_t mxpiu = 0; |
123 | Double_t mxmeno = 0; | |
124 | Double_t mypiu = 0; | |
125 | Double_t mymeno = 0; | |
126 | Double_t r=0; | |
127 | ||
128 | Int_t np1=0, np2=0; | |
32e63e47 | 129 | for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) { |
7d62fb64 | 130 | detTypeRec.ResetRecPoints(); |
fe4280bf | 131 | tr->GetEvent(i); |
132 | npoints = recpoints->GetEntries(); | |
133 | for (Int_t ipoint=0;ipoint<npoints;ipoint++) { | |
134 | pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint); | |
32e63e47 | 135 | /* |
00a7cc50 | 136 | l[0]=pnt->GetDetLocalX(); |
fe4280bf | 137 | l[1]=0; |
00a7cc50 | 138 | l[2]=pnt->GetDetLocalZ(); |
fe4280bf | 139 | g2->LtoG(i, l, p); |
32e63e47 | 140 | */ |
141 | pnt->GetGlobalXYZ(p); | |
fe4280bf | 142 | r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2)); |
143 | ||
144 | if(i<80 && TMath::Abs(p[2])<14.35) { | |
145 | y1[np1]=p[1]; | |
146 | x1[np1]=p[0]; | |
147 | z1[np1]=p[2]; | |
148 | if(p[0]>0) mxpiu++; | |
149 | if(p[0]<0) mxmeno++; | |
150 | if(p[1]>0) mypiu++; | |
151 | if(p[1]<0) mymeno++; | |
152 | np1++; | |
2257f27e | 153 | } |
fe4280bf | 154 | if(i>=80 && TMath::Abs(p[2])<14.35) { |
155 | y2[np2]=p[1]; | |
156 | x2[np2]=p[0]; | |
157 | z2[np2]=p[2]; | |
158 | np2++; | |
88cb7938 | 159 | } |
160 | } | |
fe4280bf | 161 | } |
c5f0f3c1 | 162 | |
83028144 | 163 | if(np1<fNpThreshold) { |
cd706e57 | 164 | Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n"); |
83028144 | 165 | Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold); |
cd706e57 | 166 | AliITSVertexerZ *dovert = new AliITSVertexerZ("default"); |
83028144 | 167 | fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber()); |
168 | delete dovert; | |
169 | return fCurrentVertex; | |
170 | } | |
171 | ||
fe4280bf | 172 | if(!np1 || !np2) { |
173 | Error("FindVertexForCurrentEvent","No points in the pixer layers"); | |
83028144 | 174 | return fCurrentVertex; |
fe4280bf | 175 | } |
cd706e57 | 176 | AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2)); |
fe4280bf | 177 | |
178 | Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno); | |
179 | Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno); | |
180 | ||
181 | Double_t x0 = 2.86*asparx; | |
182 | Double_t y0 = 2.86*aspary; | |
183 | ||
cd706e57 | 184 | AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0)); |
fe4280bf | 185 | |
186 | for(Int_t i=0;i<np1;i++) { | |
187 | x1[i]-=x0; | |
188 | y1[i]-=y0; | |
189 | PhiFunc(x1[i],y1[i],phi1[i]); | |
190 | r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]); | |
88cb7938 | 191 | } |
fe4280bf | 192 | for(Int_t i=0;i<np2;i++) { |
193 | x2[i]-=x0; | |
194 | y2[i]-=y0; | |
195 | PhiFunc(x2[i],y2[i],phi2[i]); | |
196 | r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]); | |
83028144 | 197 | } |
fe4280bf | 198 | |
199 | Int_t nbinxy=400; | |
200 | Int_t nbinz=400; | |
201 | TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5); | |
202 | TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5); | |
203 | TH1F *hzv=new TH1F("hzv","",nbinz,-15,15); | |
204 | ||
205 | Double_t dphi; | |
206 | for(Int_t j=0; j<np1; j++) { | |
207 | for(Int_t k=0; k<np2; k++) { | |
208 | dphi=TMath::Abs(phi2[k]-phi1[j]); | |
88cb7938 | 209 | if(dphi>180) dphi = 360-dphi; |
fe4280bf | 210 | if(dphi>fMaxDeltaPhi) continue; |
211 | hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1)); | |
212 | } | |
88cb7938 | 213 | } |
fe4280bf | 214 | |
215 | // Fitting ... | |
216 | Double_t max; | |
7d62fb64 | 217 | Int_t binMax; |
b5d712e2 | 218 | Double_t maxcenter; |
83028144 | 219 | |
fe4280bf | 220 | max = hzv->GetMaximum(); |
7d62fb64 | 221 | binMax=hzv->GetMaximumBin(); |
222 | maxcenter=hzv->GetBinCenter(binMax); | |
fe4280bf | 223 | Double_t dxy=1.5; |
83028144 | 224 | |
b5d712e2 | 225 | TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy); |
fe4280bf | 226 | fz->SetParameter(0,max); |
b5d712e2 | 227 | fz->SetParameter(1,maxcenter); |
fe4280bf | 228 | fz->SetParameter(2,0.1); |
229 | fz->SetLineColor(kRed); | |
230 | hzv->Fit("fz","RQ0"); | |
231 | ||
232 | if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) { | |
233 | 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"); | |
234 | Double_t position[3]={0,0,fz->GetParameter(1)}; | |
235 | Double_t resolution[3]={0,0,0}; | |
236 | Double_t snr[3]={0,0,0}; | |
237 | Char_t name[30]; | |
cd706e57 | 238 | AliDebug(1,Form("Vertex found for event %d",evnumber)); |
fe4280bf | 239 | sprintf(name,"Vertex"); |
240 | fCurrentVertex = new AliESDVertex(position,resolution,snr,name); | |
83028144 | 241 | return fCurrentVertex; |
88cb7938 | 242 | } |
83028144 | 243 | |
fe4280bf | 244 | for(Int_t j=0; j<np1; j++) { |
245 | for(Int_t k=0; k<np2; k++) { | |
246 | if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue; | |
247 | if(y2[k]==y1[j]) continue; | |
248 | hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k])); | |
249 | if(x2[k]==x1[j]) continue; | |
250 | hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k])); | |
88cb7938 | 251 | } |
252 | } | |
fe4280bf | 253 | |
254 | 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); | |
255 | fx->SetParameter(0,100); | |
256 | Double_t dist=0.3; | |
b5d712e2 | 257 | Double_t xapprox=FindMaxAround(x0,hxv,dist); |
cd706e57 | 258 | AliDebug(1,Form("xapprox = %f",xapprox)); |
b5d712e2 | 259 | fx->SetParameter(1,xapprox); |
260 | Double_t difcentroid=0.07; | |
261 | fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid); | |
fe4280bf | 262 | fx->SetParameter(2,0.1); |
263 | fx->SetLineColor(kRed); | |
264 | hxv->Fit("fx","RQW0"); | |
83028144 | 265 | |
fe4280bf | 266 | 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); |
267 | fy->SetParameter(0,100); | |
b5d712e2 | 268 | Double_t yapprox=FindMaxAround(y0,hyv,dist); |
cd706e57 | 269 | AliDebug(1,Form("yapprox = %f",yapprox)); |
b5d712e2 | 270 | fy->SetParameter(1,yapprox); |
271 | fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid); | |
fe4280bf | 272 | fy->SetParameter(2,0.1); |
273 | fy->SetLineColor(kRed); | |
274 | hyv->Fit("fy","RQW0"); | |
83028144 | 275 | |
83028144 | 276 | delete [] z1; |
277 | delete [] z2; | |
278 | delete [] y1; | |
279 | delete [] y2; | |
280 | delete [] x1; | |
281 | delete [] x2; | |
88cb7938 | 282 | delete [] r1; |
283 | delete [] r2; | |
284 | delete [] phi1; | |
285 | delete [] phi2; | |
fe4280bf | 286 | delete hxv; |
287 | delete hyv; | |
288 | delete hzv; | |
289 | ||
290 | Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)}; | |
291 | Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)}; | |
292 | Double_t snr[3]={0,0,0}; | |
293 | ||
88cb7938 | 294 | Char_t name[30]; |
cd706e57 | 295 | AliDebug(1,Form("Vertex found for event %d",evnumber)); |
88cb7938 | 296 | sprintf(name,"Vertex"); |
d681bb2d | 297 | fCurrentVertex = new AliESDVertex(position,resolution,snr,name); |
7d62fb64 | 298 | |
88cb7938 | 299 | return fCurrentVertex; |
c5f0f3c1 | 300 | } |
301 | ||
302 | //______________________________________________________________________ | |
fe4280bf | 303 | void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) { |
b5d712e2 | 304 | // Method for the computation of the phi angle given x and y. The result is in degrees. |
fe4280bf | 305 | if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578); |
306 | if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180; | |
307 | if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180; | |
308 | if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;; | |
c5f0f3c1 | 309 | } |
310 | ||
311 | //______________________________________________________________________ | |
312 | void AliITSVertexerIons::FindVertices(){ | |
313 | // computes the vertices of the events in the range FirstEvent - LastEvent | |
88cb7938 | 314 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); |
83028144 | 315 | AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader"); |
316 | itsloader->LoadRecPoints("read"); | |
ab6a6f40 | 317 | for(Int_t i=fFirstEvent;i<=fLastEvent;i++){ |
88cb7938 | 318 | rl->GetEvent(i); |
c5f0f3c1 | 319 | FindVertexForCurrentEvent(i); |
320 | if(fCurrentVertex){ | |
321 | WriteCurrentVertex(); | |
322 | } | |
323 | else { | |
cd706e57 | 324 | AliDebug(1,Form("Vertex not found for event %d",i)); |
c5f0f3c1 | 325 | } |
326 | } | |
c5f0f3c1 | 327 | } |
328 | ||
329 | //________________________________________________________ | |
330 | void AliITSVertexerIons::PrintStatus() const { | |
331 | // Print current status | |
88cb7938 | 332 | cout <<"=======================================================\n"; |
c5f0f3c1 | 333 | cout<<"First event to be processed "<<fFirstEvent; |
334 | cout<<"\n Last event to be processed "<<fLastEvent<<endl; | |
335 | if(fCurrentVertex)fCurrentVertex->PrintStatus(); | |
336 | } | |
fe4280bf | 337 | |
338 | //________________________________________________________ | |
339 | Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) { | |
b5d712e2 | 340 | // It finds the maximum of h within point-distance and distance+point. |
341 | Int_t maxcontent=0; | |
342 | Int_t maxbin=0; | |
fe4280bf | 343 | for(Int_t i=0;i<h->GetNbinsX();i++) { |
344 | Int_t content=(Int_t)h->GetBinContent(i); | |
345 | Double_t center=(Double_t)h->GetBinCenter(i); | |
1360bc6f | 346 | if(TMath::Abs(center-point)>distance) continue; |
b5d712e2 | 347 | if(content>maxcontent) {maxcontent=content;maxbin=i;} |
fe4280bf | 348 | } |
b5d712e2 | 349 | Double_t max=h->GetBinCenter(maxbin); |
fe4280bf | 350 | return max; |
351 | } | |
352 | ||
353 |