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"
14 //////////////////////////////////////////////////////////////////////
15 // AliITSVertexerIons is a class for full 3D primary vertex //
16 // finding optimized for Ion-Ion interactions //
21 // Written by Giuseppe Lo Re and Francesco Riggi //
22 // Giuseppe.Lore@ct.infn.it //
23 // Franco.Riggi@ct.infn.it //
25 // Release date: Mar 2004 //
28 //////////////////////////////////////////////////////////////////////
31 ClassImp(AliITSVertexerIons)
33 //______________________________________________________________________
34 AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
38 // Default Constructor
46 //______________________________________________________________________
47 AliITSVertexerIons::~AliITSVertexerIons() {
52 //______________________________________________________________________
53 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(TTree *itsClusterTree){
54 // Defines the AliESDVertex for the current event
58 AliITSDetTypeRec detTypeRec;
60 detTypeRec.SetTreeAddressR(itsClusterTree);
62 TClonesArray *recpoints = detTypeRec.RecPoints();
67 Int_t nopoints1=40000;
68 Int_t nopoints2=40000;
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];
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);
97 l[0]=pnt->GetDetLocalX();
99 l[2]=pnt->GetDetLocalZ();
102 pnt->GetGlobalXYZ(p);
103 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
105 if(i<80 && TMath::Abs(p[2])<14.35) {
115 if(i>=80 && TMath::Abs(p[2])<14.35) {
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);
130 return fCurrentVertex;
134 Error("FindVertexForCurrentEvent","No points in the pixer layers");
135 return fCurrentVertex;
137 AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
139 Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
140 Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
142 Double_t x0 = 2.86*asparx;
143 Double_t y0 = 2.86*aspary;
145 AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));
147 for(Int_t i=0;i<np1;i++) {
150 PhiFunc(x1[i],y1[i],phi1[i]);
151 r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
153 for(Int_t i=0;i<np2;i++) {
156 PhiFunc(x2[i],y2[i],phi2[i]);
157 r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
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);
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));
181 max = hzv->GetMaximum();
182 binMax=hzv->GetMaximumBin();
183 maxcenter=hzv->GetBinCenter(binMax);
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");
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};
199 sprintf(name,"Vertex");
200 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
201 return fCurrentVertex;
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]));
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);
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");
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");
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};
255 sprintf(name,"Vertex");
256 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
258 return fCurrentVertex;
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;;
270 //________________________________________________________
271 void AliITSVertexerIons::PrintStatus() const {
272 // Print current status
273 cout <<"=======================================================\n";
274 if(fCurrentVertex)fCurrentVertex->PrintStatus();
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.
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;}
288 Double_t max=h->GetBinCenter(maxbin);