1 #include <AliITSVertexerIons.h>
13 #include "AliITSgeom.h"
14 #include "AliITSLoader.h"
15 #include "AliITSRecPoint.h"
16 #include "AliGenerator.h"
22 #include <AliESDVertex.h>
23 #include <TObjArray.h>
25 #include <AliITSVertexerPPZ.h>
26 //////////////////////////////////////////////////////////////////////
27 // AliITSVertexerIons is a class for full 3D primary vertex //
28 // finding optimized for Ion-Ion interactions //
33 // Written by Giuseppe Lo Re and Francesco Riggi //
35 // Giuseppe.Lore@ct.infn.it //
36 // Franco.Riggi@ct.infn.it //
38 // Release date: May 2001 //
41 //////////////////////////////////////////////////////////////////////
43 ClassImp(AliITSVertexerIons)
47 //______________________________________________________________________
48 AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
49 // Default Constructor
55 //______________________________________________________________________
56 AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
57 // Standard constructor
64 //______________________________________________________________________
65 AliITSVertexerIons::~AliITSVertexerIons() {
70 //______________________________________________________________________
71 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
72 // Defines the AliESDVertex for the current event
75 Double_t resolution[3];
77 AliRunLoader *rl = AliRunLoader::GetRunLoader();
78 AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
80 fITS =(AliITS *)gAlice->GetDetector("ITS");
82 Error("FindVertexForCurrentEvent","AliITS object was not found");
83 return fCurrentVertex;
86 fITS->SetTreeAddress();
87 AliITSgeom *g2 = fITS->GetITSgeom();
88 TClonesArray *recpoints = fITS->RecPoints();
89 TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
92 branch = itsloader->TreeR()->GetBranch("Clusters");
93 branch->SetAddress(&clusters);
96 branch = itsloader->TreeR()->GetBranch("ITSRecPoints");
106 //------------ Rough Vertex evaluation ---------------------------------
108 Int_t i,npoints,ipoint,j,k,max,binmax;
117 TH1F *hITSz1 = new TH1F("hITSz1","",100,-14.35,14.35);
118 TTree *tr = itsloader->TreeR();
119 for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
121 fITS->ResetRecPoints();
124 Clusters2RecPoints(clusters,i,recpoints);
126 npoints = recpoints->GetEntries();
127 for (ipoint=0;ipoint<npoints;ipoint++) {
128 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
133 if(i<80 && TMath::Abs(p[2])<14.35) {
140 if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++;
144 nopoints1 = (Int_t)(hITSz1->GetEntries());
145 if (nopoints1 == 0) {
147 return fCurrentVertex;
149 aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
150 aspar[1] = (mypiu-mymeno)/(mypiu+mymeno);
152 vzero[0] = 5.24441*aspar[0];
153 vzero[1] = 5.24441*aspar[1];
155 zCentroid= TMath::Abs(hITSz1->GetMean());
156 vzero[2] = -5.31040e-02+1.42198*zCentroid+7.44718e-01*TMath::Power(zCentroid,2)
157 -5.73426e-01*TMath::Power(zCentroid,3)+2.01500e-01*TMath::Power(zCentroid,4)
158 -3.34118e-02*TMath::Power(zCentroid,5)+2.20816e-03*TMath::Power(zCentroid,6);
160 if(hITSz1->GetMean()<0) vzero[2] = -vzero[2];
162 //cout << "\nXvzero: " << vzero[0] << " cm" << "";
163 //cout << "\nYvzero: " << vzero[1] << " cm" << "";
164 //cout << "\nZvzero: " << vzero[2] << " cm" << "\n";
168 Double_t dphi,r,deltaZ,a,b;
173 Double_t deltaPhiZ = 0.08;
177 for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
178 gAlice->Field()->Field(origin,mag);
180 deltaPhiZ = deltaPhiZ*TMath::Abs(mag[2])/2;
181 Double_t deltaPhiXY = 1.0;
183 //cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n";
184 //cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n";
186 Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
187 z1=new Double_t[nopoints1];
188 z2=new Double_t[nopoints2];
189 y1=new Double_t[nopoints1];
190 y2=new Double_t[nopoints2];
191 x1=new Double_t[nopoints1];
192 x2=new Double_t[nopoints2];
193 phi1=new Double_t[nopoints1];
194 phi2=new Double_t[nopoints2];
195 r1=new Double_t[nopoints1];
196 r2=new Double_t[nopoints2];
198 deltaZ = 4.91617e-01+2.67567e-02*vzero[2]+1.49626e-02*TMath::Power(vzero[2],2);
199 Float_t multFactorZ = 28000./(Float_t)nopoints1;
200 Int_t nbin=(Int_t)((deltaZ/0.005)/multFactorZ);
202 Int_t *vectorBinZ,*vectorBinXY;
203 vectorBinZ=new Int_t[nbin];
204 vectorBinXY=new Int_t[nbinxy];
207 Double_t sigma,averagebg;
209 TH1D *hITSZv = new TH1D("hITSZv","",nbin,vzero[2]-deltaZ,vzero[2]+deltaZ);
210 TH1D *hITSXv = new TH1D("hITSXv","",nbinxy,-3,3);
211 TH1D *hITSYv = new TH1D("hITSYv","",nbinxy,-3,3);
213 //cout << "deltaZeta: " << deltaZ << " cm" << "\n";
218 hITSZv->Add(hITSZv,-1);
219 hITSXv->Add(hITSXv,-1);
220 hITSYv->Add(hITSYv,-1);
225 for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++)
227 fITS->ResetRecPoints();
228 itsloader->TreeR()->GetEvent(i);
230 Clusters2RecPoints(clusters,i,recpoints);
232 npoints = recpoints->GetEntries();
233 for (ipoint=0;ipoint<npoints;ipoint++) {
235 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
241 if(i<80 && TMath::Abs(p[2])<14.35) {
244 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
249 phi1[np1]=PhiFunc(p);
253 if(i>=80 && TMath::Abs(p[2])<14.35) {
256 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
261 phi2[np2]=PhiFunc(p);
268 //------------------ Correlation between rec points ----------------------
271 if(np1<fNpThreshold) {
272 /* cout << "Warning: AliITSVertexerIons finder is not reliable for low multiplicity events. May be you have to use AliITSVertexerPPZ.\n";
273 position[0]=vzero[0];
274 position[1]=vzero[1];
275 position[2]=vzero[2];
283 sprintf(name,"Vertex");
284 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
285 return fCurrentVertex;*/
287 Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
288 Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
289 AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("null");
290 fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
292 return fCurrentVertex;
296 for(j=0; j<(np2)-1; j++) {
297 for(k=0; k<(np1)-1; k++) {
298 dphi=TMath::Abs(phi1[k]-phi2[j]);
299 if(dphi>180) dphi = 360-dphi;
300 if(dphi<deltaPhiZ && TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-vzero[2])
301 <deltaZ) hITSZv->Fill(z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1));
305 //cout << "\nNumber of used pairs: \n";
309 max=(Int_t) hITSZv->GetMaximum();
310 binmax=hITSZv->GetMaximumBin();
312 for(i=0;i<nbin;i++) vectorBinZ[i]=(Int_t)hITSZv->GetBinContent(i);
313 for(i=0;i<10;i++) f1=f1+vectorBinZ[i]/10;
314 for(i=nbin-10;i<nbin;i++) f2=f2+vectorBinZ[i]/10;
316 for(i=0;i<nbin;i++) {
317 if(vectorBinZ[i]-averagebg>(max-averagebg)*0.4 &&
318 vectorBinZ[i]-averagebg<(max-averagebg)*0.7) {
319 sigma=hITSZv->GetBinCenter(binmax)-hITSZv->GetBinCenter(i);
320 sigma=TMath::Abs(sigma);
321 if(sigma==0) sigma=0.05;
326 TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);
328 fz->SetParameter(0,max);
329 if(niter==0) {Double_t temp = hITSZv->GetBinCenter(binmax); vzero[2]=temp;}
330 fz->SetParameter(1,vzero[2]);
331 fz->SetParameter(2,sigma);
332 fz->SetParameter(3,averagebg);
333 fz->SetParLimits(2,0,999);
334 fz->SetParLimits(3,0,999);
336 hITSZv->Fit("fz","RQ0");
338 snr[2] = fz->GetParameter(0)/fz->GetParameter(3);
340 Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n");
341 Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n");
343 return fCurrentVertex;
346 position[2] = fz->GetParameter(1);
349 position[2]=position[2]-TMath::Abs(position[2])*1.11/10000;
353 position[2]=position[2]+TMath::Abs(position[2])*1.11/10000;
356 resolution[2] = fz->GetParError(1);
358 for(j=0; j<(np2)-1; j++) {
359 for(k=0; k<(np1)-1; k++) {
360 dphi=TMath::Abs(phi1[k]-phi2[j]);
361 if(dphi>180) dphi = 360-dphi;
362 if(dphi>deltaPhiXY) continue;
363 if(TMath::Abs((z1[k]-(z2[j]-z1[k])/((r2[j]/r1[k])-1))-position[2])
365 hITSXv->Fill(vzero[0]+(x2[j]-((x2[j]-x1[k])/(y2[j]-y1[k]))*y2[j]));
366 hITSYv->Fill(vzero[1]+(y2[j]-((y2[j]-y1[k])/(x2[j]-x1[k]))*x2[j]));
371 TF1 *fx = new TF1("fx","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[0]-0.5,vzero[0]+0.5);
373 max=(Int_t) hITSXv->GetMaximum();
374 binmax=hITSXv->GetMaximumBin();
377 for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSXv->GetBinContent(i);
378 for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
379 for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
381 for(i=0;i<nbinxy;i++) {
382 if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 &&
383 vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
384 sigma=hITSXv->GetBinCenter(binmax)-hITSXv->GetBinCenter(i);
385 sigma=TMath::Abs(sigma);
386 if(sigma==0) sigma=0.05;
391 fx->SetParameter(0,max);
392 fx->SetParameter(1,vzero[0]);
393 fx->SetParameter(2,sigma);
394 fx->SetParameter(3,averagebg);
396 hITSXv->Fit("fx","RQ0");
398 snr[0] = fx->GetParameter(0)/fx->GetParameter(3);
401 Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n");
402 Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n");
404 return fCurrentVertex;
407 position[0]=fx->GetParameter(1);
408 resolution[0]=fx->GetParError(1);
411 TF1 *fy = new TF1("fy","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",vzero[1]-0.5,vzero[1]+0.5);
413 max=(Int_t) hITSYv->GetMaximum();
414 binmax=hITSYv->GetMaximumBin();
417 for(i=0;i<nbinxy;i++) vectorBinXY[i]=(Int_t)hITSYv->GetBinContent(i);
418 for(i=0;i<10;i++) f1=f1+vectorBinXY[i]/10;
419 for(i=nbinxy-10;i<nbinxy;i++) f2=f2+vectorBinXY[i]/10;
421 for(i=0;i<nbinxy;i++) {
422 if(vectorBinXY[i]-averagebg>(max-averagebg)*0.4 &&
423 vectorBinXY[i]-averagebg<(max-averagebg)*0.7) {
424 sigma=hITSYv->GetBinCenter(binmax)-hITSYv->GetBinCenter(i);
425 sigma=TMath::Abs(sigma);
426 if(sigma==0) sigma=0.05;
430 fy->SetParameter(0,max);
431 fy->SetParameter(1,vzero[1]);
432 fy->SetParameter(2,sigma);
433 fy->SetParameter(3,averagebg);
435 hITSYv->Fit("fy","RQ0");
437 snr[1] = fy->GetParameter(0)/fy->GetParameter(3);
439 Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n");
440 Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n");
442 return fCurrentVertex;
445 position[1]=fy->GetParameter(1);
446 resolution[1]=fy->GetParError(1);
450 /*cout << "iter = " << niter << " -> x = " << position[0] << endl;
451 cout << "iter = " << niter << " -> y = " << position[1] << endl;
452 cout << "iter = " << niter << " -> z = " << position[2] << endl;
457 vzero[0] = position[0];
458 vzero[1] = position[1];
459 vzero[2] = position[2];
461 if(niter<3) goto start; // number of iterations
464 cout<<"Number of iterations: "<<niter<<endl;
465 cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl;
477 delete [] vectorBinZ;
478 delete [] vectorBinXY;
484 // sprintf(name,"Vertex_%d",evnumber);
485 if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
486 sprintf(name,"Vertex");
487 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
488 return fCurrentVertex;
491 //______________________________________________________________________
492 Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
495 if(p[1]>0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578);
496 if(p[1]>0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
497 if(p[1]<0 && p[0]<0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+180;
498 if(p[1]<0 && p[0]>0) phi=(TMath::ATan((Double_t)(p[1]/p[0]))*57.29578)+360;
502 //______________________________________________________________________
503 void AliITSVertexerIons::FindVertices(){
504 // computes the vertices of the events in the range FirstEvent - LastEvent
505 AliRunLoader *rl = AliRunLoader::GetRunLoader();
506 AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
507 itsloader->LoadRecPoints("read");
508 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
510 FindVertexForCurrentEvent(i);
512 WriteCurrentVertex();
516 cout<<"Vertex not found for event "<<i<<endl;
522 //________________________________________________________
523 void AliITSVertexerIons::PrintStatus() const {
524 // Print current status
525 cout <<"=======================================================\n";
526 cout <<" Debug flag: "<<fDebug<<endl;
527 cout<<"First event to be processed "<<fFirstEvent;
528 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
529 if(fCurrentVertex)fCurrentVertex->PrintStatus();