1 #include "AliITSDetTypeRec.h"
2 #include "AliITSVertexerIons.h"
3 #include "AliITSVertexerPPZ.h"
4 #include "AliESDVertex.h"
5 #include "AliITSgeom.h"
6 #include "AliITSLoader.h"
7 #include "AliITSRecPoint.h"
12 //////////////////////////////////////////////////////////////////////
13 // AliITSVertexerIons is a class for full 3D primary vertex //
14 // finding optimized for Ion-Ion interactions //
19 // Written by Giuseppe Lo Re and Francesco Riggi //
20 // Giuseppe.Lore@ct.infn.it //
21 // Franco.Riggi@ct.infn.it //
23 // Release date: Mar 2004 //
26 //////////////////////////////////////////////////////////////////////
29 ClassImp(AliITSVertexerIons)
31 //______________________________________________________________________
32 AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
36 // Default Constructor
44 //______________________________________________________________________
45 AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn),
49 // Standard constructor
57 //______________________________________________________________________
58 AliITSVertexerIons::AliITSVertexerIons(const AliITSVertexerIons &source):AliITSVertexer(source) {
60 // Copies are not allowed. The method is protected to avoid misuse.
61 Error("AliITSVertexerIons","Copy constructor not allowed\n");
64 //_________________________________________________________________________
65 //AliITSVertexerIons& AliITSVertexerIons::operator=(const AliITSVertexerIons &/*source*/) {
66 // Assignment operator
67 // Assignment is not allowed. The method is protected to avoid misuse.
68 //Error("= operator","Assignment operator not allowed\n");
73 //______________________________________________________________________
74 AliITSVertexerIons::~AliITSVertexerIons() {
79 //______________________________________________________________________
80 AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
81 // Defines the AliESDVertex for the current event
85 AliRunLoader *rl = AliRunLoader::GetRunLoader();
86 AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
87 TDirectory * olddir = gDirectory;
89 AliITSgeom* g2 = (AliITSgeom*)gDirectory->Get("AliITSgeom");
92 TTree *tr = itsloader->TreeR();
93 AliITSDetTypeRec detTypeRec;
95 detTypeRec.SetTreeAddressR(tr);
97 TClonesArray *recpoints = detTypeRec.RecPoints();
102 Int_t nopoints1=40000;
103 Int_t nopoints2=40000;
106 Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
107 z1=new Double_t[nopoints1];
108 z2=new Double_t[nopoints2];
109 y1=new Double_t[nopoints1];
110 y2=new Double_t[nopoints2];
111 x1=new Double_t[nopoints1];
112 x2=new Double_t[nopoints2];
113 phi1=new Double_t[nopoints1];
114 phi2=new Double_t[nopoints2];
115 r1=new Double_t[nopoints1];
116 r2=new Double_t[nopoints2];
125 for(Int_t i=g2->GetStartSPD();i<=g2->GetLastSPD();i++) {
126 detTypeRec.ResetRecPoints();
128 npoints = recpoints->GetEntries();
129 for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
130 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
131 l[0]=pnt->GetDetLocalX();
133 l[2]=pnt->GetDetLocalZ();
135 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
137 if(i<80 && TMath::Abs(p[2])<14.35) {
147 if(i>=80 && TMath::Abs(p[2])<14.35) {
156 if(np1<fNpThreshold) {
157 Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerPPZ with default parameters...\n");
158 Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
159 AliITSVertexerPPZ *dovert = new AliITSVertexerPPZ("default");
160 fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
162 return fCurrentVertex;
166 Error("FindVertexForCurrentEvent","No points in the pixer layers");
167 return fCurrentVertex;
169 if(fDebug>0) cout << "N. points layer 1 and 2 : " << np1 << " " << np2 << endl;
171 Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
172 Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
174 Double_t x0 = 2.86*asparx;
175 Double_t y0 = 2.86*aspary;
177 if(fDebug>0) cout << "Rough xy vertex = " << x0 << " " << y0 << endl;
179 for(Int_t i=0;i<np1;i++) {
182 PhiFunc(x1[i],y1[i],phi1[i]);
183 r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
185 for(Int_t i=0;i<np2;i++) {
188 PhiFunc(x2[i],y2[i],phi2[i]);
189 r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
194 TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
195 TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
196 TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
199 for(Int_t j=0; j<np1; j++) {
200 for(Int_t k=0; k<np2; k++) {
201 dphi=TMath::Abs(phi2[k]-phi1[j]);
202 if(dphi>180) dphi = 360-dphi;
203 if(dphi>fMaxDeltaPhi) continue;
204 hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
213 max = hzv->GetMaximum();
214 binMax=hzv->GetMaximumBin();
215 maxcenter=hzv->GetBinCenter(binMax);
218 TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);
219 fz->SetParameter(0,max);
220 fz->SetParameter(1,maxcenter);
221 fz->SetParameter(2,0.1);
222 fz->SetLineColor(kRed);
223 hzv->Fit("fz","RQ0");
225 if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
226 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");
227 Double_t position[3]={0,0,fz->GetParameter(1)};
228 Double_t resolution[3]={0,0,0};
229 Double_t snr[3]={0,0,0};
231 if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
232 sprintf(name,"Vertex");
233 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
234 return fCurrentVertex;
237 for(Int_t j=0; j<np1; j++) {
238 for(Int_t k=0; k<np2; k++) {
239 if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
240 if(y2[k]==y1[j]) continue;
241 hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
242 if(x2[k]==x1[j]) continue;
243 hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
247 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);
248 fx->SetParameter(0,100);
250 Double_t xapprox=FindMaxAround(x0,hxv,dist);
251 if(fDebug>0) cout << "xapprox = " << xapprox << endl;
252 fx->SetParameter(1,xapprox);
253 Double_t difcentroid=0.07;
254 fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
255 fx->SetParameter(2,0.1);
256 fx->SetLineColor(kRed);
257 hxv->Fit("fx","RQW0");
259 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);
260 fy->SetParameter(0,100);
261 Double_t yapprox=FindMaxAround(y0,hyv,dist);
262 if(fDebug>0) cout << "yapprox = " << yapprox << endl;
263 fy->SetParameter(1,yapprox);
264 fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
265 fy->SetParameter(2,0.1);
266 fy->SetLineColor(kRed);
267 hyv->Fit("fy","RQW0");
283 Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
284 Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
285 Double_t snr[3]={0,0,0};
288 if(fDebug>0)Info("FindVertexForCurrentEvent","Vertex found for event %d",evnumber);
289 sprintf(name,"Vertex");
290 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
292 return fCurrentVertex;
295 //______________________________________________________________________
296 void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
297 // Method for the computation of the phi angle given x and y. The result is in degrees.
298 if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
299 if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
300 if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
301 if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
304 //______________________________________________________________________
305 void AliITSVertexerIons::FindVertices(){
306 // computes the vertices of the events in the range FirstEvent - LastEvent
307 AliRunLoader *rl = AliRunLoader::GetRunLoader();
308 AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
309 itsloader->LoadRecPoints("read");
310 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
312 FindVertexForCurrentEvent(i);
314 WriteCurrentVertex();
318 cout<<"Vertex not found for event "<<i<<endl;
324 //________________________________________________________
325 void AliITSVertexerIons::PrintStatus() const {
326 // Print current status
327 cout <<"=======================================================\n";
328 cout <<" Debug flag: "<<fDebug<<endl;
329 cout<<"First event to be processed "<<fFirstEvent;
330 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
331 if(fCurrentVertex)fCurrentVertex->PrintStatus();
334 //________________________________________________________
335 Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
336 // It finds the maximum of h within point-distance and distance+point.
339 for(Int_t i=0;i<h->GetNbinsX();i++) {
340 Int_t content=(Int_t)h->GetBinContent(i);
341 Double_t center=(Double_t)h->GetBinCenter(i);
342 if(TMath::Abs(center-point)>distance) continue;
343 if(content>maxcontent) {maxcontent=content;maxbin=i;}
345 Double_t max=h->GetBinCenter(maxbin);