1 #include <AliITSVertexerIons.h>
13 #include "AliITSgeom.h"
14 #include "AliITSLoader.h"
15 #include "AliITSRecPoint.h"
16 #include "AliGenerator.h"
22 #include <AliITSVertex.h>
23 #include <TObjArray.h>
25 //////////////////////////////////////////////////////////////////////
26 // AliITSVertexerIons is a class for full 3D primary vertex //
27 // finding optimized for Ion-Ion interactions //
32 // Written by Giuseppe Lo Re and Francesco Riggi //
33 // Giuseppe.Lore@ct.infn.it //
35 // Franco.Riggi@ct.infn.it //
37 // Release date: May 2001 //
40 //////////////////////////////////////////////////////////////////////
42 ClassImp(AliITSVertexerIons)
46 //______________________________________________________________________
47 AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
48 // Default Constructor
53 //______________________________________________________________________
54 AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
55 // Standard constructor
61 //______________________________________________________________________
62 AliITSVertexerIons::~AliITSVertexerIons() {
67 //______________________________________________________________________
68 AliITSVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
69 // Defines the AliITSVertex for the current event
72 Double_t Resolution[3];
74 AliRunLoader *rl = AliRunLoader::GetRunLoader();
75 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
77 fITS =(AliITS *)gAlice->GetDetector("ITS");
79 Error("FindVertexForCurrentEvent","AliITS object was not found");
80 return fCurrentVertex;
83 fITS->SetTreeAddress();
84 AliITSgeom *g2 = fITS->GetITSgeom();
85 TClonesArray *recpoints = fITS->RecPoints();
93 //------------ Rough Vertex evaluation ---------------------------------
95 Int_t i,npoints,ipoint,j,k,max,BinMax;
104 TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35);
105 TTree *TR = ITSloader->TreeR();
106 for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
108 fITS->ResetRecPoints();
110 npoints = recpoints->GetEntries();
111 for (ipoint=0;ipoint<npoints;ipoint++) {
112 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
117 if(i<80 && TMath::Abs(p[2])<14.35) {
124 if(i>=80 && TMath::Abs(p[2])<14.35) NoPoints2++;
128 NoPoints1 = (Int_t)(hITSz1->GetEntries());
130 AsPar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
131 AsPar[1] = (mypiu-mymeno)/(mypiu+mymeno);
133 Vzero[0] = 5.24441*AsPar[0];
134 Vzero[1] = 5.24441*AsPar[1];
136 ZCentroid= TMath::Abs(hITSz1->GetMean());
137 Vzero[2] = -5.31040e-02+1.42198*ZCentroid+7.44718e-01*TMath::Power(ZCentroid,2)
138 -5.73426e-01*TMath::Power(ZCentroid,3)+2.01500e-01*TMath::Power(ZCentroid,4)
139 -3.34118e-02*TMath::Power(ZCentroid,5)+2.20816e-03*TMath::Power(ZCentroid,6);
141 if(hITSz1->GetMean()<0) Vzero[2] = -Vzero[2];
143 /*cout << "\nXvzero: " << Vzero[0] << " cm" << "";
144 cout << "\nYvzero: " << Vzero[1] << " cm" << "";
145 cout << "\nZvzero: " << Vzero[2] << " cm" << "\n";*/
149 Double_t dphi,r,DeltaZ,a,b;
154 Double_t DeltaPhiZ = 0.08;
158 for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
159 gAlice->Field()->Field(origin,B);
161 DeltaPhiZ = DeltaPhiZ*B[2]/2;
162 Double_t DeltaPhiXY = 1.0;
164 // cout << "\nDeltaPhiZ: " << DeltaPhiZ << " deg" << "\n";
165 // cout << "DeltaPhiXY: " << DeltaPhiXY << " deg" << "\n";
167 Double_t *Z1, *Z2, *Y1, *Y2, *X1, *X2, *phi1, *phi2, *r1, *r2;
168 Z1=new Double_t[NoPoints1];
169 Z2=new Double_t[NoPoints2];
170 Y1=new Double_t[NoPoints1];
171 Y2=new Double_t[NoPoints2];
172 X1=new Double_t[NoPoints1];
173 X2=new Double_t[NoPoints2];
174 phi1=new Double_t[NoPoints1];
175 phi2=new Double_t[NoPoints2];
176 r1=new Double_t[NoPoints1];
177 r2=new Double_t[NoPoints2];
179 DeltaZ = 4.91617e-01+2.67567e-02*Vzero[2]+1.49626e-02*TMath::Power(Vzero[2],2);
180 Float_t MulFactorZ = 28000./(Float_t)NoPoints1;
181 Int_t nbin=(Int_t)((DeltaZ/0.005)/MulFactorZ);
183 Int_t *VectorBinZ,*VectorBinXY;
184 VectorBinZ=new Int_t[nbin];
185 VectorBinXY=new Int_t[nbinxy];
188 Double_t sigma,MediaFondo;
190 TH1D *hITSZv = new TH1D("hITSZv","",nbin,Vzero[2]-DeltaZ,Vzero[2]+DeltaZ);
191 TH1D *hITSXv = new TH1D("hITSXv","",nbinxy,-3,3);
192 TH1D *hITSYv = new TH1D("hITSYv","",nbinxy,-3,3);
194 // cout << "DeltaZeta: " << DeltaZ << " cm" << "\n";
199 hITSZv->Add(hITSZv,-1.);
200 hITSXv->Add(hITSXv,-1.);
201 hITSYv->Add(hITSYv,-1.);
207 for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
209 fITS->ResetRecPoints();
210 ITSloader->TreeR()->GetEvent(i);
211 npoints = recpoints->GetEntries();
212 for (ipoint=0;ipoint<npoints;ipoint++) {
214 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
220 if(i<80 && TMath::Abs(p[2])<14.35) {
223 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
228 phi1[np1]=PhiFunc(p);
232 if(i>=80 && TMath::Abs(p[2])<14.35) {
235 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
240 phi2[np2]=PhiFunc(p);
247 //------------------ Correlation between rec points ----------------------
249 //cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl;
251 for(j=0; j<(np2)-1; j++) {
252 for(k=0; k<(np1)-1; k++) {
253 dphi=TMath::Abs(phi1[k]-phi2[j]);
254 if(dphi>180) dphi = 360-dphi;
255 if(dphi<DeltaPhiZ && TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Vzero[2])
256 <DeltaZ) hITSZv->Fill(Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1));
260 //cout << "\nNumber of used pairs: \n";
261 //cout << hITSZv->GetEntries() << '\n' << '\n';
264 max=(Int_t) hITSZv->GetMaximum();
265 BinMax=hITSZv->GetMaximumBin();
267 for(i=0;i<nbin;i++) VectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
268 for(i=0;i<10;i++) f1=f1+VectorBinZ[i]/10;
269 for(i=nbin-10;i<nbin;i++) f2=f2+VectorBinZ[i]/10;
270 MediaFondo=(f1+f2)/2;
271 for(i=0;i<nbin;i++) {
272 if(VectorBinZ[i]-MediaFondo>(max-MediaFondo)*0.4 &&
273 VectorBinZ[i]-MediaFondo<(max-MediaFondo)*0.7) {
274 sigma=hITSZv->GetBinCenter(BinMax)-hITSZv->GetBinCenter(i);
275 sigma=TMath::Abs(sigma);
276 if(sigma==0) sigma=0.05;
280 /*cout << "f1 " <<f1 <<endl;
281 cout << "f2 " <<f2 <<endl;
282 cout << "GetMaximumBin " <<hITSZv->GetMaximumBin() <<endl;
283 cout << "nbin " <<nbin <<endl;
284 cout << "max " << hITSZv->GetBinContent(BinMax)<<endl;
285 cout << "sigma " <<sigma <<endl;
286 cout << "Fondo " << MediaFondo<< endl;*/
288 TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);
290 fz->SetParameter(0,max);
291 if(niter==0) {Double_t temp = hITSZv->GetBinCenter(BinMax); Vzero[2]=temp;}
292 fz->SetParameter(1,Vzero[2]);
293 fz->SetParameter(2,sigma);
294 fz->SetParameter(3,MediaFondo);
295 fz->SetParLimits(2,0,999);
296 fz->SetParLimits(3,0,999);
298 hITSZv->Fit("fz","RMEQ0");
300 SNR[2] = fz->GetParameter(0)/fz->GetParameter(3);
302 cout << "\nNegative Signal to noise ratio for z!!!" << endl;
303 cout << "The algorithm cannot find the z vertex position." << endl;
308 Position[2] = fz->GetParameter(1);
311 Position[2]=Position[2]-TMath::Abs(Position[2])*1.11/10000;
315 Position[2]=Position[2]+TMath::Abs(Position[2])*1.11/10000;
318 Resolution[2] = fz->GetParError(1);
322 for(j=0; j<(np2)-1; j++) {
323 for(k=0; k<(np1)-1; k++) {
324 dphi=TMath::Abs(phi1[k]-phi2[j]);
325 if(dphi>180) dphi = 360-dphi;
326 if(dphi>DeltaPhiXY) continue;
327 if(TMath::Abs((Z1[k]-(Z2[j]-Z1[k])/((r2[j]/r1[k])-1))-Position[2])
329 hITSXv->Fill(Vzero[0]+(X2[j]-((X2[j]-X1[k])/(Y2[j]-Y1[k]))*Y2[j]));
330 hITSYv->Fill(Vzero[1]+(Y2[j]-((Y2[j]-Y1[k])/(X2[j]-X1[k]))*X2[j]));
336 ("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[0]-0.5,Vzero[0]+0.5);
338 max=(Int_t) hITSXv->GetMaximum();
339 BinMax=hITSXv->GetMaximumBin();
342 for(i=0;i<nbinxy;i++) VectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
343 for(i=0;i<10;i++) f1=f1+VectorBinXY[i]/10;
344 for(i=nbinxy-10;i<nbinxy;i++) f2=f2+VectorBinXY[i]/10;
345 MediaFondo=(f1+f2)/2;
346 for(i=0;i<nbinxy;i++) {
347 if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 &&
348 VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
349 sigma=hITSXv->GetBinCenter(BinMax)-hITSXv->GetBinCenter(i);
350 sigma=TMath::Abs(sigma);
351 if(sigma==0) sigma=0.05;
355 /*cout << "f1 " <<f1 <<endl;
356 cout << "f2 " <<f2 <<endl;
357 cout << "GetMaximumBin " <<hITSXv->GetMaximumBin() <<endl;
358 cout << "max " << hITSXv->GetBinContent(BinMax)<<endl;
359 cout << "sigma " <<sigma <<endl;
360 cout << "Fondo " << MediaFondo<< endl;
361 cout << "nbinxy " <<nbinxy <<endl;*/
363 fx->SetParameter(0,max);
364 fx->SetParameter(1,Vzero[0]);
365 fx->SetParameter(2,sigma);
366 fx->SetParameter(3,MediaFondo);
368 hITSXv->Fit("fx","RMEQ0");
370 SNR[0] = fx->GetParameter(0)/fx->GetParameter(3);
372 cout << "\nNegative Signal to noise ratio for x!!!" << endl;
373 cout << "The algorithm cannot find the x vertex position." << endl;
378 Position[0]=fx->GetParameter(1);
379 Resolution[0]=fx->GetParError(1);
382 TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",Vzero[1]-0.5,Vzero[1]+0.5);
384 max=(Int_t) hITSYv->GetMaximum();
385 BinMax=hITSYv->GetMaximumBin();
388 for(i=0;i<nbinxy;i++) VectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
389 for(i=0;i<10;i++) f1=f1+VectorBinXY[i]/10;
390 for(i=nbinxy-10;i<nbinxy;i++) f2=f2+VectorBinXY[i]/10;
391 MediaFondo=(f1+f2)/2;
392 for(i=0;i<nbinxy;i++) {
393 if(VectorBinXY[i]-MediaFondo>(max-MediaFondo)*0.4 &&
394 VectorBinXY[i]-MediaFondo<(max-MediaFondo)*0.7) {
395 sigma=hITSYv->GetBinCenter(BinMax)-hITSYv->GetBinCenter(i);
396 sigma=TMath::Abs(sigma);
397 if(sigma==0) sigma=0.05;
401 /*cout << "f1 " <<f1 <<endl;
402 cout << "f2 " <<f2 <<endl;
403 cout << "GetMaximumBin " <<hITSYv->GetMaximumBin() <<endl;
404 cout << "max " << hITSYv->GetBinContent(BinMax)<<endl;
405 cout << "sigma " <<sigma <<endl;
406 cout << "Fondo " << MediaFondo<< endl;
407 cout << "nbinxy " <<nbinxy <<endl;*/
409 fy->SetParameter(0,max);
410 fy->SetParameter(1,Vzero[1]);
411 fy->SetParameter(2,sigma);
412 fy->SetParameter(3,MediaFondo);
414 hITSYv->Fit("fy","RMEQ0");
416 SNR[1] = fy->GetParameter(0)/fy->GetParameter(3);
418 cout << "\nNegative Signal to noise ratio for y!!!" << endl;
419 cout << "The algorithm cannot find the y vertex position." << endl;
424 Position[1]=fy->GetParameter(1);
425 Resolution[1]=fy->GetParError(1);
430 Vzero[0] = Position[0];
431 Vzero[1] = Position[1];
432 Vzero[2] = Position[2];
434 if(niter<3) goto start; // number of iterations
437 cout<<"Number of iterations: "<<niter<<endl;
449 delete [] VectorBinZ;
450 delete [] VectorBinXY;
456 // sprintf(name,"Vertex_%d",evnumber);
457 sprintf(name,"Vertex");
458 fCurrentVertex = new AliITSVertex(Position,Resolution,SNR,name);
459 return fCurrentVertex;
462 //______________________________________________________________________
463 Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
466 if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
467 if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
468 if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
469 if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
473 //______________________________________________________________________
474 void AliITSVertexerIons::FindVertices(){
475 // computes the vertices of the events in the range FirstEvent - LastEvent
476 AliRunLoader *rl = AliRunLoader::GetRunLoader();
477 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
478 ITSloader->LoadRecPoints("read");
479 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
481 FindVertexForCurrentEvent(i);
483 WriteCurrentVertex();
487 cout<<"Vertex not found for event "<<i<<endl;
493 //________________________________________________________
494 void AliITSVertexerIons::PrintStatus() const {
495 // Print current status
496 cout <<"=======================================================\n";
497 cout <<" Debug flag: "<<fDebug<<endl;
498 cout<<"First event to be processed "<<fFirstEvent;
499 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
500 if(fCurrentVertex)fCurrentVertex->PrintStatus();