Implementation of single event reconstruction. Removal of the run loaders from the...
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerIons.cxx
CommitLineData
6cae184e 1#include "AliRunLoader.h"
7d62fb64 2#include "AliITSDetTypeRec.h"
b5d712e2 3#include "AliITSVertexerIons.h"
cd706e57 4#include "AliITSVertexerZ.h"
b5d712e2 5#include "AliESDVertex.h"
32e63e47 6#include "AliITSgeomTGeo.h"
88cb7938 7#include "AliITSLoader.h"
c5f0f3c1 8#include "AliITSRecPoint.h"
cd706e57 9#include "AliLog.h"
b5d712e2 10#include <TTree.h>
c5f0f3c1 11#include <TH1.h>
12#include <TF1.h>
b5d712e2 13
c5f0f3c1 14//////////////////////////////////////////////////////////////////////
15// AliITSVertexerIons is a class for full 3D primary vertex //
16// finding optimized for Ion-Ion interactions //
17// //
18// //
19// //
20// //
21// Written by Giuseppe Lo Re and Francesco Riggi //
22// Giuseppe.Lore@ct.infn.it //
c5f0f3c1 23// Franco.Riggi@ct.infn.it //
24// //
fe4280bf 25// Release date: Mar 2004 //
c5f0f3c1 26// //
27// //
28//////////////////////////////////////////////////////////////////////
29
c5f0f3c1 30
b5d712e2 31ClassImp(AliITSVertexerIons)
32
fe4280bf 33//______________________________________________________________________
8221b41b 34 AliITSVertexerIons::AliITSVertexerIons():AliITSVertexer(),
35fNpThreshold(0),
36fMaxDeltaPhi(0),
37fMaxDeltaZ(0){
c5f0f3c1 38 // Default Constructor
39
7d62fb64 40 //fITS = 0;
83028144 41 SetNpThreshold();
fe4280bf 42 SetMaxDeltaPhi();
43 SetMaxDeltaZ();
c5f0f3c1 44}
45
46//______________________________________________________________________
8221b41b 47AliITSVertexerIons::AliITSVertexerIons(TString fn):AliITSVertexer(fn),
48fNpThreshold(0),
49fMaxDeltaPhi(0),
50fMaxDeltaZ(0) {
c5f0f3c1 51 // Standard constructor
88cb7938 52
7d62fb64 53 //fITS = 0;
83028144 54 SetNpThreshold();
fe4280bf 55 SetMaxDeltaPhi();
b5d712e2 56 SetMaxDeltaZ();
57}
8221b41b 58/*
b5d712e2 59//______________________________________________________________________
60AliITSVertexerIons::AliITSVertexerIons(const AliITSVertexerIons &source):AliITSVertexer(source) {
61 // Copy constructor
62 // Copies are not allowed. The method is protected to avoid misuse.
63 Error("AliITSVertexerIons","Copy constructor not allowed\n");
64}
8221b41b 65*/
b5d712e2 66//_________________________________________________________________________
8221b41b 67//AliITSVertexerIons& AliITSVertexerIons::operator=(const AliITSVertexerIons &/*source*/) {
b5d712e2 68 // Assignment operator
69 // Assignment is not allowed. The method is protected to avoid misuse.
8221b41b 70 //Error("= operator","Assignment operator not allowed\n");
71 //return *this;
72//}
c5f0f3c1 73
74
75//______________________________________________________________________
76AliITSVertexerIons::~AliITSVertexerIons() {
77 // Default Destructor
7d62fb64 78 //fITS = 0;
c5f0f3c1 79}
80
81//______________________________________________________________________
b5d712e2 82AliESDVertex* AliITSVertexerIons::FindVertexForCurrentEvent(Int_t evnumber){
83// Defines the AliESDVertex for the current event
fe4280bf 84
c5f0f3c1 85 fCurrentVertex = 0;
fe4280bf 86
88cb7938 87 AliRunLoader *rl = AliRunLoader::GetRunLoader();
7d62fb64 88 AliITSLoader* itsloader = (AliITSLoader*)rl->GetLoader("ITSLoader");
32e63e47 89 /*
7d62fb64 90 TDirectory * olddir = gDirectory;
91 rl->CdGAFile();
92 AliITSgeom* g2 = (AliITSgeom*)gDirectory->Get("AliITSgeom");
93 olddir->cd();
32e63e47 94 */
7d62fb64 95
fe4280bf 96 TTree *tr = itsloader->TreeR();
7d62fb64 97 AliITSDetTypeRec detTypeRec;
98
99 detTypeRec.SetTreeAddressR(tr);
100
101 TClonesArray *recpoints = detTypeRec.RecPoints();
102 AliITSRecPoint *pnt;
103
88cb7938 104
fe4280bf 105 Int_t npoints=0;
106 Int_t nopoints1=40000;
107 Int_t nopoints2=40000;
32e63e47 108 Float_t p[3];
c5f0f3c1 109
83028144 110 Double_t *z1, *z2, *y1, *y2, *x1, *x2, *phi1, *phi2, *r1, *r2;
111 z1=new Double_t[nopoints1];
112 z2=new Double_t[nopoints2];
113 y1=new Double_t[nopoints1];
114 y2=new Double_t[nopoints2];
115 x1=new Double_t[nopoints1];
116 x2=new Double_t[nopoints2];
117 phi1=new Double_t[nopoints1];
118 phi2=new Double_t[nopoints2];
119 r1=new Double_t[nopoints1];
120 r2=new Double_t[nopoints2];
c5f0f3c1 121
fe4280bf 122 Double_t mxpiu = 0;
123 Double_t mxmeno = 0;
124 Double_t mypiu = 0;
125 Double_t mymeno = 0;
126 Double_t r=0;
127
128 Int_t np1=0, np2=0;
32e63e47 129 for(Int_t i=AliITSgeomTGeo::GetModuleIndex(1,1,1);i<=AliITSgeomTGeo::GetModuleIndex(2,1,1)-1;i++) {
7d62fb64 130 detTypeRec.ResetRecPoints();
fe4280bf 131 tr->GetEvent(i);
132 npoints = recpoints->GetEntries();
133 for (Int_t ipoint=0;ipoint<npoints;ipoint++) {
134 pnt = (AliITSRecPoint*)recpoints->UncheckedAt(ipoint);
32e63e47 135 /*
00a7cc50 136 l[0]=pnt->GetDetLocalX();
fe4280bf 137 l[1]=0;
00a7cc50 138 l[2]=pnt->GetDetLocalZ();
fe4280bf 139 g2->LtoG(i, l, p);
32e63e47 140 */
141 pnt->GetGlobalXYZ(p);
fe4280bf 142 r=TMath::Sqrt(TMath::Power(p[0],2)+TMath::Power(p[1],2));
143
144 if(i<80 && TMath::Abs(p[2])<14.35) {
145 y1[np1]=p[1];
146 x1[np1]=p[0];
147 z1[np1]=p[2];
148 if(p[0]>0) mxpiu++;
149 if(p[0]<0) mxmeno++;
150 if(p[1]>0) mypiu++;
151 if(p[1]<0) mymeno++;
152 np1++;
2257f27e 153 }
fe4280bf 154 if(i>=80 && TMath::Abs(p[2])<14.35) {
155 y2[np2]=p[1];
156 x2[np2]=p[0];
157 z2[np2]=p[2];
158 np2++;
88cb7938 159 }
160 }
fe4280bf 161 }
c5f0f3c1 162
83028144 163 if(np1<fNpThreshold) {
cd706e57 164 Warning("FindVertexForCurrentEvent","AliITSVertexerIons finder is not reliable for low multiplicity events. Switching to AliITSVertexerZ with default parameters...\n");
83028144 165 Warning("FindVertexForCurrentEvent","N rec points = %d - Threshold is %d",np1,fNpThreshold);
cd706e57 166 AliITSVertexerZ *dovert = new AliITSVertexerZ("default");
83028144 167 fCurrentVertex =dovert->FindVertexForCurrentEvent(rl->GetEventNumber());
168 delete dovert;
169 return fCurrentVertex;
170 }
171
fe4280bf 172 if(!np1 || !np2) {
173 Error("FindVertexForCurrentEvent","No points in the pixer layers");
83028144 174 return fCurrentVertex;
fe4280bf 175 }
cd706e57 176 AliDebug(1,Form("N. points layer 1 and 2 : %d %d",np1,np2));
fe4280bf 177
178 Double_t asparx = (mxpiu-mxmeno)/(mxpiu+mxmeno);
179 Double_t aspary = (mypiu-mymeno)/(mypiu+mymeno);
180
181 Double_t x0 = 2.86*asparx;
182 Double_t y0 = 2.86*aspary;
183
cd706e57 184 AliDebug(1,Form("Rough xy vertex = %f %f",x0,y0));
fe4280bf 185
186 for(Int_t i=0;i<np1;i++) {
187 x1[i]-=x0;
188 y1[i]-=y0;
189 PhiFunc(x1[i],y1[i],phi1[i]);
190 r1[i]=TMath::Sqrt(x1[i]*x1[i]+y1[i]*y1[i]);
88cb7938 191 }
fe4280bf 192 for(Int_t i=0;i<np2;i++) {
193 x2[i]-=x0;
194 y2[i]-=y0;
195 PhiFunc(x2[i],y2[i],phi2[i]);
196 r2[i]=TMath::Sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
83028144 197 }
fe4280bf 198
199 Int_t nbinxy=400;
200 Int_t nbinz=400;
201 TH1F *hxv=new TH1F("hxv","",nbinxy,-5,5);
202 TH1F *hyv=new TH1F("hyv","",nbinxy,-5,5);
203 TH1F *hzv=new TH1F("hzv","",nbinz,-15,15);
204
205 Double_t dphi;
206 for(Int_t j=0; j<np1; j++) {
207 for(Int_t k=0; k<np2; k++) {
208 dphi=TMath::Abs(phi2[k]-phi1[j]);
88cb7938 209 if(dphi>180) dphi = 360-dphi;
fe4280bf 210 if(dphi>fMaxDeltaPhi) continue;
211 hzv->Fill(z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1));
212 }
88cb7938 213 }
fe4280bf 214
215 // Fitting ...
216 Double_t max;
7d62fb64 217 Int_t binMax;
b5d712e2 218 Double_t maxcenter;
83028144 219
fe4280bf 220 max = hzv->GetMaximum();
7d62fb64 221 binMax=hzv->GetMaximumBin();
222 maxcenter=hzv->GetBinCenter(binMax);
fe4280bf 223 Double_t dxy=1.5;
83028144 224
b5d712e2 225 TF1 *fz = new TF1 ("fz","([0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2])))+[3]",maxcenter-dxy,maxcenter+dxy);
fe4280bf 226 fz->SetParameter(0,max);
b5d712e2 227 fz->SetParameter(1,maxcenter);
fe4280bf 228 fz->SetParameter(2,0.1);
229 fz->SetLineColor(kRed);
230 hzv->Fit("fz","RQ0");
231
232 if(TMath::Abs(x0)>0.4 || TMath::Abs(y0)>0.4) {
233 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");
234 Double_t position[3]={0,0,fz->GetParameter(1)};
235 Double_t resolution[3]={0,0,0};
236 Double_t snr[3]={0,0,0};
237 Char_t name[30];
cd706e57 238 AliDebug(1,Form("Vertex found for event %d",evnumber));
fe4280bf 239 sprintf(name,"Vertex");
240 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
83028144 241 return fCurrentVertex;
88cb7938 242 }
83028144 243
fe4280bf 244 for(Int_t j=0; j<np1; j++) {
245 for(Int_t k=0; k<np2; k++) {
246 if(TMath::Abs((z1[j]-(z2[k]-z1[j])/((r2[k]/r1[j])-1))-fz->GetParameter(1))>fMaxDeltaZ) continue;
247 if(y2[k]==y1[j]) continue;
248 hxv->Fill(x0+(x2[k]-((x2[k]-x1[j])/(y2[k]-y1[j]))*y2[k]));
249 if(x2[k]==x1[j]) continue;
250 hyv->Fill(y0+(y2[k]-((y2[k]-y1[j])/(x2[k]-x1[j]))*x2[k]));
88cb7938 251 }
252 }
fe4280bf 253
254 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);
255 fx->SetParameter(0,100);
256 Double_t dist=0.3;
b5d712e2 257 Double_t xapprox=FindMaxAround(x0,hxv,dist);
cd706e57 258 AliDebug(1,Form("xapprox = %f",xapprox));
b5d712e2 259 fx->SetParameter(1,xapprox);
260 Double_t difcentroid=0.07;
261 fx->SetParLimits(1,xapprox-difcentroid,xapprox+difcentroid);
fe4280bf 262 fx->SetParameter(2,0.1);
263 fx->SetLineColor(kRed);
264 hxv->Fit("fx","RQW0");
83028144 265
fe4280bf 266 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);
267 fy->SetParameter(0,100);
b5d712e2 268 Double_t yapprox=FindMaxAround(y0,hyv,dist);
cd706e57 269 AliDebug(1,Form("yapprox = %f",yapprox));
b5d712e2 270 fy->SetParameter(1,yapprox);
271 fy->SetParLimits(1,yapprox-difcentroid,yapprox+difcentroid);
fe4280bf 272 fy->SetParameter(2,0.1);
273 fy->SetLineColor(kRed);
274 hyv->Fit("fy","RQW0");
83028144 275
83028144 276 delete [] z1;
277 delete [] z2;
278 delete [] y1;
279 delete [] y2;
280 delete [] x1;
281 delete [] x2;
88cb7938 282 delete [] r1;
283 delete [] r2;
284 delete [] phi1;
285 delete [] phi2;
fe4280bf 286 delete hxv;
287 delete hyv;
288 delete hzv;
289
290 Double_t position[3]={fx->GetParameter(1),fy->GetParameter(1),fz->GetParameter(1)};
291 Double_t resolution[3]={fx->GetParameter(2),fy->GetParameter(2),fz->GetParameter(2)};
292 Double_t snr[3]={0,0,0};
293
88cb7938 294 Char_t name[30];
cd706e57 295 AliDebug(1,Form("Vertex found for event %d",evnumber));
88cb7938 296 sprintf(name,"Vertex");
d681bb2d 297 fCurrentVertex = new AliESDVertex(position,resolution,snr,name);
7d62fb64 298
88cb7938 299 return fCurrentVertex;
c5f0f3c1 300}
301
302//______________________________________________________________________
fe4280bf 303void AliITSVertexerIons::PhiFunc(Double_t &x,Double_t &y,Double_t &phi) {
b5d712e2 304 // Method for the computation of the phi angle given x and y. The result is in degrees.
fe4280bf 305 if(y>0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578);
306 if(y>0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
307 if(y<0 && x<0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+180;
308 if(y<0 && x>0) phi=(TMath::ATan((Double_t)(y/x))*57.29578)+360;;
c5f0f3c1 309}
310
311//______________________________________________________________________
312void AliITSVertexerIons::FindVertices(){
313 // computes the vertices of the events in the range FirstEvent - LastEvent
88cb7938 314 AliRunLoader *rl = AliRunLoader::GetRunLoader();
83028144 315 AliITSLoader* itsloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
316 itsloader->LoadRecPoints("read");
ab6a6f40 317 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
88cb7938 318 rl->GetEvent(i);
c5f0f3c1 319 FindVertexForCurrentEvent(i);
320 if(fCurrentVertex){
321 WriteCurrentVertex();
322 }
323 else {
cd706e57 324 AliDebug(1,Form("Vertex not found for event %d",i));
c5f0f3c1 325 }
326 }
c5f0f3c1 327}
328
329//________________________________________________________
330void AliITSVertexerIons::PrintStatus() const {
331 // Print current status
88cb7938 332 cout <<"=======================================================\n";
c5f0f3c1 333 cout<<"First event to be processed "<<fFirstEvent;
334 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
335 if(fCurrentVertex)fCurrentVertex->PrintStatus();
336}
fe4280bf 337
338//________________________________________________________
339Double_t AliITSVertexerIons::FindMaxAround(Double_t point, TH1F *h, Double_t distance) {
b5d712e2 340 // It finds the maximum of h within point-distance and distance+point.
341 Int_t maxcontent=0;
342 Int_t maxbin=0;
fe4280bf 343 for(Int_t i=0;i<h->GetNbinsX();i++) {
344 Int_t content=(Int_t)h->GetBinContent(i);
345 Double_t center=(Double_t)h->GetBinCenter(i);
1360bc6f 346 if(TMath::Abs(center-point)>distance) continue;
b5d712e2 347 if(content>maxcontent) {maxcontent=content;maxbin=i;}
fe4280bf 348 }
b5d712e2 349 Double_t max=h->GetBinCenter(maxbin);
fe4280bf 350 return max;
351}
352
353