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