Primary vertex reconstruction and standalone ITS tracking in the reconstruction chain
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx
1 #include <AliITSVertexerIons.h>
2 #include "stdlib.h"
3 #include <TMath.h>
4 #include <TRandom.h>
5 #include <TObjArray.h>
6 #include <TROOT.h>
7 #include <TFile.h>
8 #include <TTree.h>
9 #include <Riostream.h>
10
11 #include "AliRun.h"
12 #include "AliITS.h"
13 #include "AliITSgeom.h"
14 #include "AliITSLoader.h"
15 #include "AliITSRecPoint.h"
16 #include "AliGenerator.h"
17 #include "AliMagF.h"
18
19 #include <TH1.h>
20 #include <TF1.h>
21 #include <TCanvas.h>
22 #include <AliESDVertex.h>
23 #include <TObjArray.h>
24 #include <TObject.h>
25 #include <AliITSVertexerPPZ.h>
26 //////////////////////////////////////////////////////////////////////
27 // AliITSVertexerIons is a class for full 3D primary vertex         //
28 // finding optimized for Ion-Ion interactions                       //
29 //                                                                  // 
30 //                                                                  //
31 //                                                                  //
32 //                                                                  //
33 // Written by Giuseppe Lo Re and Francesco Riggi                    //
34 //
35 // Giuseppe.Lore@ct.infn.it                                         //
36 // Franco.Riggi@ct.infn.it                                          //
37 //                                                                  //
38 // Release date: May 2001                                           //
39 //                                                                  //
40 //                                                                  //       
41 //////////////////////////////////////////////////////////////////////
42
43 ClassImp(AliITSVertexerIons)
44
45
46
47   //______________________________________________________________________
48   AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer() {
49   // Default Constructor
50
51   fITS = 0;
52   SetNpThreshold();
53 }
54
55 //______________________________________________________________________
56 AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn) {
57   // Standard constructor
58   
59   fITS = 0;
60   SetNpThreshold();
61 }
62
63
64 //______________________________________________________________________
65 AliITSVertexerIons::~AliITSVertexerIons() {
66   // Default Destructor
67   fITS = 0;
68 }
69
70 //______________________________________________________________________
71 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
72   // Defines the AliESDVertex for the current event
73   fCurrentVertex = 0;
74   Double_t position[3];
75   Double_t resolution[3];
76   Double_t snr[3];
77   AliRunLoader *rl = AliRunLoader::GetRunLoader();
78   AliITSLoader* itsloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
79   if(!fITS) { 
80     fITS =(AliITS *)gAlice->GetDetector("ITS");
81     if(!fITS) {
82       Error("FindVertexForCurrentEvent","AliITS object was not found");
83       return fCurrentVertex;
84     }
85   }
86   fITS->SetTreeAddress();
87   AliITSgeom *g2 = fITS->GetITSgeom(); 
88   TClonesArray  *recpoints = fITS->RecPoints();
89   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
90   TBranch *branch;
91   if(fUseV2Clusters){
92     branch = itsloader->TreeR()->GetBranch("Clusters");
93     branch->SetAddress(&clusters);
94   }
95   else {
96     branch = itsloader->TreeR()->GetBranch("ITSRecPoints");
97   }
98
99   AliITSRecPoint *pnt;
100
101   Int_t nopoints1 = 0;                                         
102   Int_t nopoints2 = 0;
103   Double_t vzero[3];
104   Double_t aspar[2];
105      
106   //------------ Rough Vertex evaluation ---------------------------------   
107
108   Int_t i,npoints,ipoint,j,k,max,binmax;
109   Double_t zCentroid;
110   Float_t l[3], p[3];
111
112   Double_t mxpiu = 0;
113   Double_t mxmeno = 0;
114   Double_t mypiu = 0;
115   Double_t mymeno = 0;
116    
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++) 
120     {
121       fITS->ResetRecPoints();
122       tr->GetEvent(i); 
123       if(fUseV2Clusters){
124         Clusters2RecPoints(clusters,i,recpoints);
125       }
126       npoints = recpoints->GetEntries();
127       for (ipoint=0;ipoint<npoints;ipoint++) {
128         pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
129         l[0]=pnt->GetX();
130         l[1]=0;
131         l[2]=pnt->GetZ();
132         g2->LtoG(i, l, p);
133         if(i<80 && TMath::Abs(p[2])<14.35) {         
134           if(p[0]>0) mxpiu++; 
135           if(p[0]<0) mxmeno++; 
136           if(p[1]>0) mypiu++; 
137           if(p[1]<0) mymeno++; 
138           hITSz1->Fill(p[2]);       
139         }
140         if(i>=80 && TMath::Abs(p[2])<14.35) nopoints2++;
141       }       
142     }  
143   
144   nopoints1 = (Int_t)(hITSz1->GetEntries()); 
145   if (nopoints1 == 0) {
146     delete hITSz1;
147     return fCurrentVertex;
148   }         
149   aspar[0] = (mxpiu-mxmeno)/(mxpiu+mxmeno);
150   aspar[1] = (mypiu-mymeno)/(mypiu+mymeno); 
151    
152   vzero[0] = 5.24441*aspar[0]; 
153   vzero[1] = 5.24441*aspar[1]; 
154
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);
159    
160   if(hITSz1->GetMean()<0) vzero[2] = -vzero[2];
161    
162   //cout << "\nXvzero: " << vzero[0] << " cm" << "";
163   //cout << "\nYvzero: " << vzero[1] << " cm" << "";
164   //cout << "\nZvzero: " << vzero[2] << " cm" << "\n";
165
166   delete hITSz1;
167
168   Double_t dphi,r,deltaZ,a,b;
169   Int_t np1=0;
170   Int_t np2=0;
171   Int_t niter=0;
172    
173   Double_t deltaPhiZ = 0.08;
174    
175   Float_t mag[3];
176   Float_t origin[3];
177   for(Int_t okm=0;okm<3;okm++) origin[okm]=0;
178   gAlice->Field()->Field(origin,mag);
179
180   deltaPhiZ = deltaPhiZ*TMath::Abs(mag[2])/2;
181   Double_t deltaPhiXY = 1.0;   
182
183   //cout << "\ndeltaPhiZ: " << deltaPhiZ << " deg" << "\n";   
184   //cout << "deltaPhiXY: " << deltaPhiXY << " deg" << "\n";   
185    
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];
197         
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);
201   Int_t nbinxy=250;
202   Int_t *vectorBinZ,*vectorBinXY;
203   vectorBinZ=new Int_t[nbin];
204   vectorBinXY=new Int_t[nbinxy];
205   Float_t f1= 0;
206   Float_t f2= 0;
207   Double_t sigma,averagebg;
208      
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);
212
213   //cout << "deltaZeta: " << deltaZ << " cm" << "\n";   
214    
215    
216  start:   
217
218   hITSZv->Add(hITSZv,-1);
219   hITSXv->Add(hITSXv,-1);
220   hITSYv->Add(hITSYv,-1);
221
222   np1=np2=0;
223
224
225   for(i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) 
226     {
227       fITS->ResetRecPoints(); 
228       itsloader->TreeR()->GetEvent(i);
229      if(fUseV2Clusters){
230         Clusters2RecPoints(clusters,i,recpoints);
231       }
232       npoints = recpoints->GetEntries();
233       for (ipoint=0;ipoint<npoints;ipoint++) {
234                
235         pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
236         l[0]=pnt->GetX();
237         l[1]=0;
238         l[2]=pnt->GetZ();
239         g2->LtoG(i, l, p);
240               
241         if(i<80 && TMath::Abs(p[2])<14.35) {
242           p[0]=p[0]-vzero[0];
243           p[1]=p[1]-vzero[1];
244           r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
245           y1[np1]=p[1];
246           x1[np1]=p[0];
247           z1[np1]=p[2]; 
248           r1[np1]=r;
249           phi1[np1]=PhiFunc(p);
250           np1++;
251         }
252
253         if(i>=80 &&  TMath::Abs(p[2])<14.35) { 
254           p[0]=p[0]-vzero[0];
255           p[1]=p[1]-vzero[1];
256           r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
257           y2[np2]=p[1];
258           x2[np2]=p[0];
259           z2[np2]=p[2];
260           r2[np2]=r;
261           phi2[np2]=PhiFunc(p);
262           np2++;
263         }
264                         
265       }
266     }
267
268   //------------------ Correlation between rec points ----------------------
269     
270   
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];
276           resolution[0]=-123;
277           resolution[1]=-123;
278           resolution[2]=-123;
279           snr[0]=-123;
280           snr[1]=-123;
281           snr[2]=-123;
282           Char_t name[30];
283           sprintf(name,"Vertex");
284           fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
285           return fCurrentVertex;*/
286
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());
291     delete dovert;
292     return fCurrentVertex;
293   }
294
295
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));        
302     }
303   }
304       
305   //cout << "\nNumber of used pairs: \n";
306
307   a = vzero[2]-deltaZ; 
308   b = vzero[2]+deltaZ; 
309   max=(Int_t) hITSZv->GetMaximum();
310   binmax=hITSZv->GetMaximumBin();
311   sigma=0;
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;
315   averagebg=(f1+f2)/2;
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;
322     }
323   }
324   
325   
326   TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",a,b);  
327   
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);
335   
336   hITSZv->Fit("fz","RQ0");
337   
338   snr[2] = fz->GetParameter(0)/fz->GetParameter(3);
339   if(snr[2]<0.) { 
340     Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for z!!!\n");
341     Error("FindVertexForCurrentEvent","The algorithm cannot find the z vertex position.\n");
342     //    exit(123456789);
343     return fCurrentVertex;
344   }
345   else {
346     position[2] = fz->GetParameter(1);
347     if(position[2]<0) 
348       {
349         position[2]=position[2]-TMath::Abs(position[2])*1.11/10000;
350       }
351     else
352       {
353         position[2]=position[2]+TMath::Abs(position[2])*1.11/10000;
354       }
355   }
356   resolution[2] = fz->GetParError(1);
357   
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])
364          <4*resolution[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]));
367       }
368     }
369   }
370   
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);  
372   
373   max=(Int_t) hITSXv->GetMaximum();
374   binmax=hITSXv->GetMaximumBin();
375   sigma=0;
376   f1=f2=0;
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;
380   averagebg=(f1+f2)/2;
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;
387     }
388   }
389   
390   
391   fx->SetParameter(0,max);
392   fx->SetParameter(1,vzero[0]);
393   fx->SetParameter(2,sigma); 
394   fx->SetParameter(3,averagebg);
395   
396   hITSXv->Fit("fx","RQ0"); 
397   
398   snr[0] = fx->GetParameter(0)/fx->GetParameter(3);
399   
400   if(snr[0]<0.) { 
401     Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for x!\n");
402     Error("FindVertexForCurrentEvent","The algorithm cannot find the x vertex position\n");
403     //   exit(123456789);  
404     return fCurrentVertex;
405   }
406   else {
407     position[0]=fx->GetParameter(1);
408     resolution[0]=fx->GetParError(1);
409   }
410   
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);  
412   
413   max=(Int_t) hITSYv->GetMaximum();
414   binmax=hITSYv->GetMaximumBin();
415   sigma=0;
416   f1=f2=0;
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;
420   averagebg=(f1+f2)/2;
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;
427     }
428   }
429   
430   fy->SetParameter(0,max);
431   fy->SetParameter(1,vzero[1]);
432   fy->SetParameter(2,sigma); 
433   fy->SetParameter(3,averagebg);
434   
435   hITSYv->Fit("fy","RQ0"); 
436   
437   snr[1] = fy->GetParameter(0)/fy->GetParameter(3);
438   if(snr[1]<0.) {   
439     Error("FindVertexForCurrentEvent","\nNegative Signal to noise ratio for y!\n");
440     Error("FindVertexForCurrentEvent","The algorithm cannot find the y vertex position.\n");
441     //    exit(123456789); 
442     return fCurrentVertex;
443   }
444   else {
445     position[1]=fy->GetParameter(1);
446     resolution[1]=fy->GetParError(1);
447   }
448   
449     
450   /*cout << "iter = " << niter << " -> x = " << position[0] << endl;
451     cout << "iter = " << niter << " -> y = " << position[1] << endl;
452     cout << "iter = " << niter << " -> z = " << position[2] << endl;
453     cout << "---\n";*/
454
455   niter++;
456    
457   vzero[0] = position[0];
458   vzero[1] = position[1];
459   vzero[2] = position[2];
460
461   if(niter<3) goto start;  // number of iterations
462    
463
464   cout<<"Number of iterations: "<<niter<<endl;
465   cout << "\nNo. of Points on the two pixel layers: "<<np1<<" "<<np2<<endl; 
466  
467   delete [] z1;
468   delete [] z2;
469   delete [] y1;
470   delete [] y2;
471   delete [] x1;
472   delete [] x2;
473   delete [] r1;
474   delete [] r2;
475   delete [] phi1;
476   delete [] phi2;
477   delete [] vectorBinZ;
478   delete [] vectorBinXY;
479    
480   delete hITSZv;
481   delete hITSXv;
482   delete hITSYv;
483   Char_t name[30];
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;
489 }
490
491 //______________________________________________________________________
492 Double_t AliITSVertexerIons::PhiFunc(Float_t p[]) {
493   Double_t phi=0;
494
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;
499   return phi;
500 }
501
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++){
509     rl->GetEvent(i);
510     FindVertexForCurrentEvent(i);
511     if(fCurrentVertex){
512       WriteCurrentVertex();
513     }
514     else {
515       if(fDebug>0){
516         cout<<"Vertex not found for event "<<i<<endl;
517       }
518     }
519   }
520 }
521
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();
530 }