Print removed
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerPPZ.cxx
CommitLineData
c95d6417 1/**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
c5f0f3c1 15#include <AliITSVertexerPPZ.h>
16#include <Riostream.h>
17#include <TArrayF.h>
18#include <TDirectory.h>
19#include <TH1F.h>
20#include <TMath.h>
21#include <TObjArray.h>
22#include <TTree.h>
23#include <TClonesArray.h>
24#include "AliRun.h"
25#include "AliITS.h"
26#include "AliITSgeom.h"
88cb7938 27#include "AliITSLoader.h"
c5f0f3c1 28#include "AliITSRecPoint.h"
29/////////////////////////////////////////////////////////////////////////
30// //
31// This class is intended to compute the Z coordinate of the //
32// primary vertex for p-p interactions, or in general when the //
33// number of secondaries is too slow to use the class //
34// AliITSVertexerIons, which in turn is optimized for A-A collsions //
35// Origin: masera@to.infn.it 9/12/2002 //
36// //
37/////////////////////////////////////////////////////////////////////////
38ClassImp(AliITSVertexerPPZ)
39
40
41
42//______________________________________________________________________
43AliITSVertexerPPZ::AliITSVertexerPPZ():AliITSVertexer() {
44 // Default Constructor
45
46 SetDiffPhiMax(0);
47 fX0 = 0.;
48 fY0 = 0.;
49 SetFirstLayerModules(0);
50 SetSecondLayerModules(0);
51 fZFound = 0;
52 fZsig = 0.;
53 fITS = 0;
c95d6417 54 SetWindow(0);
c5f0f3c1 55}
56
88cb7938 57AliITSVertexerPPZ::AliITSVertexerPPZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn) {
c5f0f3c1 58 // Standard constructor
c5f0f3c1 59 SetDiffPhiMax();
60 fX0 = x0;
61 fY0 = y0;
62 SetFirstLayerModules();
63 SetSecondLayerModules();
64 fZFound = 0;
65 fZsig = 0.;
66 fITS = 0;
c95d6417 67 SetWindow();
c5f0f3c1 68}
69
70//______________________________________________________________________
71AliITSVertexerPPZ::~AliITSVertexerPPZ() {
72 // Default Destructor
73 fITS = 0;
74}
75
76//________________________________________________________
77void AliITSVertexerPPZ::EvalZ(TH1F *hist,Int_t sepa, Int_t ncoinc, TArrayF *zval) {
c95d6417 78 Float_t DeltaVal = hist->GetBinWidth(1)*fWindow; // max window in Z for searching
c5f0f3c1 79 fZFound=0;
80 fZsig=0;
81 Int_t N=0;
82 Int_t NbinNotZero=0;
83 Float_t totst = 0.;
84 Float_t totst2 = 0.;
85 Float_t curz = 0.;
86 for(Int_t i=1;i<=sepa;i++){
87 Float_t cont=hist->GetBinContent(i);
c95d6417 88 if(cont!=0)curz = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i);
c5f0f3c1 89 totst+=cont;
90 totst2+=cont*cont;
91 N++;
92 if(cont!=0)NbinNotZero++;
93 }
94 if(NbinNotZero==0){fZFound=-100; fZsig=-100; return;}
c95d6417 95 if(NbinNotZero==1){
96 fZFound = curz;
97 fZsig=0;
d681bb2d 98 fCurrentVertex = new AliESDVertex(fZFound,fZsig,NbinNotZero);
c95d6417 99 return;
100 }
101 Float_t errsq = totst2/(N-1)-totst*totst/N/(N-1);
102 if(errsq>=0){
103 totst2=TMath::Sqrt(errsq);
104 }
105 else {
106 Error("EvalZ","Negative variance: %d - %12.7f - %12.7f",N,totst2,totst);
107 fZFound=-100;
108 fZsig=-100;
109 return;
110 }
c5f0f3c1 111 totst /= N;
c95d6417 112 Float_t cut = totst+totst2*2.;
113 if(fDebug>1)cout<<"totst, totst2, cut: "<<totst<<", "<<totst2<<", "<<cut<<endl;
c5f0f3c1 114 Float_t val1=hist->GetBinLowEdge(sepa);
115 Float_t val2=hist->GetBinLowEdge(1);
c95d6417 116 Float_t valm = 0.;
117 Float_t zmax = 0.;
c5f0f3c1 118 for(Int_t i=1;i<=sepa;i++){
119 Float_t cont=hist->GetBinContent(i);
c95d6417 120 if(cont>valm){
121 valm = cont;
122 zmax = hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(1);
123 }
124 if(cont>cut){
c5f0f3c1 125 curz=hist->GetBinLowEdge(i);
126 if(curz<val1)val1=curz;
127 if(curz>val2)val2=curz;
128 }
129 }
130 val2+=hist->GetBinWidth(1);
c95d6417 131 if((val2-val1)>DeltaVal){
132 val1 = zmax-DeltaVal/2.;
133 val2 = zmax+DeltaVal/2.;
134 if(fDebug>0)cout<<"val1 and val2 recomputed\n";
135 }
136 if(fDebug>0)cout<<"Values for Z finding: "<<val1<<" "<<val2<<" "<<val2-val1<<endl;
c5f0f3c1 137 fZFound=0;
138 fZsig=0;
139 N=0;
140 for(Int_t i=0; i<ncoinc; i++){
141 Float_t z=zval->At(i);
142 if(z<val1)continue;
143 if(z>val2)continue;
c95d6417 144
c5f0f3c1 145 fZFound+=z;
146 fZsig+=z*z;
c95d6417 147 /* weights defined by the curvature
148 Float_t wei = 1./curv->At(i);
149 fZFound+=z*wei;
150 fZsig+=wei;
151 */
c5f0f3c1 152 N++;
153 }
c95d6417 154 if(N<1){fZFound=-110; fZsig=-110; return;}
155 if(N==1){
156 fZsig = 0;
d681bb2d 157 fCurrentVertex = new AliESDVertex(fZFound,fZsig,N);
c95d6417 158 return;
159 }
160 errsq = (fZsig/(N-1)-fZFound*fZFound/N/(N-1))/N;
161 if(errsq>=0.){
162 fZsig=TMath::Sqrt(errsq);
163 }
164 else {
165 Error("evalZ","Negative variance: %d - %12.7f %12.7f",N,fZsig,fZFound);
166 fZsig=0;
167 }
c5f0f3c1 168 fZFound=fZFound/N;
c95d6417 169 /* weights defined by the curvature
170 fZsig=1./fZsig;
171 fZFound*=fZsig;
172 fZsig = TMath::Sqrt(fZsig);
173 */
d681bb2d 174 fCurrentVertex = new AliESDVertex(fZFound,fZsig,N);
c5f0f3c1 175}
176
177//______________________________________________________________________
d681bb2d 178AliESDVertex* AliITSVertexerPPZ::FindVertexForCurrentEvent(Int_t evnumber){
179 // Defines the AliESDVertex for the current event
c5f0f3c1 180 fCurrentVertex = 0;
181 fZFound = -999;
182 fZsig = -999;
88cb7938 183 AliRunLoader *rl =AliRunLoader::GetRunLoader();
184 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
c5f0f3c1 185 if(!fITS) {
186 fITS = (AliITS*)gAlice->GetModule("ITS");
187 if(!fITS) {
188 Error("FindVertexForCurrentEvent","AliITS object was not found");
189 return fCurrentVertex;
190 }
191 }
88cb7938 192 fITS->SetTreeAddress();
c5f0f3c1 193 AliITSgeom *geom = fITS->GetITSgeom();
194 if(!geom) {
195 Error("FindVertexForCurrentEvent","ITS geometry is not defined");
196 return fCurrentVertex;
197 }
2257f27e 198 TTree *tR=0;
199 TClonesArray *itsRec = 0;
c5f0f3c1 200 Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
201 Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
202 Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
203 Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
204
2257f27e 205 tR = ITSloader->TreeR();
206 if(!tR){
c5f0f3c1 207 Error("FindVertexForCurrentEvent","TreeR not found");
208 return fCurrentVertex;
209 }
2257f27e 210 itsRec = fITS->RecPoints(); // uses slow points or fast points if slow are
c5f0f3c1 211 // missing
2257f27e 212 TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
213 TBranch *branch;
214 if(fUseV2Clusters){
215 branch = tR->GetBranch("Clusters");
216 branch->SetAddress(&clusters);
217 }
218 else {
219 branch = tR->GetBranch("ITSRecPoints");
220 if(!branch){
221 branch = tR->GetBranch("ITSRecPointsF");
222 }
c5f0f3c1 223 }
224 if(!branch){
225 Error("FindVertexForCurrentEvent","branch for ITS rec points not found");
226 return fCurrentVertex;
227 }
228 Float_t zave=0;
229 Float_t rmszav=0;
230 Float_t zave2=0;
231 Int_t firipixe=0;
232 for(Int_t module= fFirstL1; module<=fLastL1;module++){
233 branch->GetEvent(module);
2257f27e 234 if(fUseV2Clusters){
235 Clusters2RecPoints(clusters,module,itsRec);
236 }
237 Int_t nrecp1 = itsRec->GetEntries();
c5f0f3c1 238 for(Int_t i=0; i<nrecp1;i++){
2257f27e 239 AliITSRecPoint *current = (AliITSRecPoint*)itsRec->At(i);
c5f0f3c1 240 lc[0]=current->GetX();
241 lc[2]=current->GetZ();
242 geom->LtoG(module,lc,gc);
243 zave+=gc[2];
244 zave2+=gc[2]*gc[2];
245 firipixe++;
246 }
247 fITS->ResetRecPoints();
248 }
249 if(firipixe>1){
250 rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
251 zave=zave/firipixe;
c95d6417 252 if(fDebug>1)cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++\n Number of firing pixels: "<<firipixe<<endl;
c5f0f3c1 253 }
254 else {
255 fZFound = -200;
256 fZsig = -200;
257 Warning("FindVertexForCurrentEvent","No rec points on first layer for this event");
258 return fCurrentVertex;
259 }
260 Float_t zlim1=zave-rmszav;
261 Float_t zlim2=zave+rmszav;
262 Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
263 zlim2=zlim1 + sepa/10.;
264 TH1F *zvdis = new TH1F("z_ev","zv distr",sepa,zlim1,zlim2);
c95d6417 265 if(fDebug>0){
266 cout<<"Z limits: "<<zlim1<<" "<<zlim2<<"; Bins= "<<sepa<<endl;
267 cout<<"Bin width: "<<zvdis->GetBinWidth(1)*10000.<<" micron\n";
268 }
c5f0f3c1 269 Int_t sizarr=100;
270 TArrayF *zval = new TArrayF(sizarr);
c95d6417 271 // TArrayF *curv = new TArrayF(sizarr);
c5f0f3c1 272 Int_t ncoinc=0;
273 for(Int_t module= fFirstL1; module<=fLastL1;module++){
274 if(fDebug>0)cout<<"processing module "<<module<<" \r";
275 branch->GetEvent(module);
2257f27e 276 if(fUseV2Clusters){
277 Clusters2RecPoints(clusters,module,itsRec);
278 }
279 Int_t nrecp1 = itsRec->GetEntries();
c5f0f3c1 280 TObjArray *poiL1 = new TObjArray(nrecp1);
2257f27e 281 for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(itsRec->At(i),i);
c5f0f3c1 282 fITS->ResetRecPoints();
283 for(Int_t i=0; i<nrecp1;i++){
284 AliITSRecPoint *current = (AliITSRecPoint*)poiL1->At(i);
285 lc[0]=current->GetX();
286 lc[2]=current->GetZ();
287 geom->LtoG(module,lc,gc);
288 gc[0]-=fX0;
289 gc[1]-=fY0;
290 Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
291 Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
292 if(phi1<0)phi1=2*TMath::Pi()+phi1;
293 if(fDebug>1)cout<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1<<" \n";
294 for(Int_t modul2=fFirstL2; modul2<=fLastL2; modul2++){
295 branch->GetEvent(modul2);
2257f27e 296 if(fUseV2Clusters){
297 Clusters2RecPoints(clusters,modul2,itsRec);
298 }
299 Int_t nrecp2 = itsRec->GetEntries();
c5f0f3c1 300 for(Int_t j=0; j<nrecp2;j++){
2257f27e 301 AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
c5f0f3c1 302 lc2[0]=recp->GetX();
303 lc2[2]=recp->GetZ();
304 geom->LtoG(modul2,lc2,gc2);
305 gc2[0]-=fX0;
306 gc2[1]-=fY0;
307 Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
308 Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
309 Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
310 if(phi2<0)phi2=2.*TMath::Pi()+phi2;
311 Float_t diff = TMath::Abs(phi2-phi1);
312 if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff;
313 if(zr0>zlim1 && zr0<zlim2){
314 if(diff<fDiffPhiMax ){
315 zvdis->Fill(zr0);
316 zval->AddAt(zr0,ncoinc);
c95d6417 317 /* uncomment these lines to use curvature as a weight
318 Float_t cu = Curv(0.,0.,gc[0],gc[1],gc2[0],gc[1]);
319 curv->AddAt(cu,ncoinc);
320 */
c5f0f3c1 321 ncoinc++;
322 if(ncoinc==(sizarr-1)){
323 sizarr+=100;
324 zval->Set(sizarr);
c95d6417 325 //uncomment next line to use curvature as weight
326 // curv->Set(sizarr);
c5f0f3c1 327 }
328 }
329 }
330 }
331 fITS->ResetRecPoints();
332 }
333 }
334 delete poiL1;
335 } // loop on modules
336 if(fDebug>0){
337 cout<<endl<<"Number of coincidences = "<<ncoinc<<endl;
338 }
c95d6417 339 // EvalZ(zvdis,sepa,ncoinc,zval,curv);
c5f0f3c1 340 EvalZ(zvdis,sepa,ncoinc,zval);
341 delete zvdis;
342 delete zval;
c95d6417 343 // delete curv;
c5f0f3c1 344 if(fCurrentVertex){
345 char name[30];
346 sprintf(name,"Vertex_%d",evnumber);
88cb7938 347 // fCurrentVertex->SetName(name);
c5f0f3c1 348 fCurrentVertex->SetTitle("vertexer: PPZ");
349 }
350 return fCurrentVertex;
351}
352
353//______________________________________________________________________
354void AliITSVertexerPPZ::FindVertices(){
355 // computes the vertices of the events in the range FirstEvent - LastEvent
88cb7938 356 AliRunLoader *rl = AliRunLoader::GetRunLoader();
357 AliITSLoader* ITSloader = (AliITSLoader*) rl->GetLoader("ITSLoader");
358 ITSloader->ReloadRecPoints();
ab6a6f40 359 for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
88cb7938 360 cout<<"Processing event "<<i<<endl;
361 rl->GetEvent(i);
c5f0f3c1 362 FindVertexForCurrentEvent(i);
363 if(fCurrentVertex){
364 WriteCurrentVertex();
365 }
366 else {
367 if(fDebug>0){
368 cout<<"Vertex not found for event "<<i<<endl;
369 cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
370 }
371 }
372 }
c5f0f3c1 373}
374
375//________________________________________________________
376void AliITSVertexerPPZ::PrintStatus() const {
377 // Print current status
378 cout <<"=======================================================\n";
379 cout <<" First layer first and last modules: "<<fFirstL1<<", ";
380 cout <<fLastL1<<endl;
381 cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
382 cout <<fLastL2<<endl;
383 cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
c95d6417 384 cout <<" Window for Z search: "<<fWindow<<endl;
c5f0f3c1 385 cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
386 cout <<" Debug flag: "<<fDebug<<endl;
c5f0f3c1 387 cout<<"First event to be processed "<<fFirstEvent;
388 cout<<"\n Last event to be processed "<<fLastEvent<<endl;
389}
c95d6417 390
391//________________________________________________________
392Float_t AliITSVertexerPPZ::Curv(Double_t x1,Double_t y1,
393 Double_t x2,Double_t y2,
394 Double_t x3,Double_t y3)
395{
396 //-----------------------------------------------------------------
397 // Initial approximation of the track curvature (Y. Belikov) squared
398 //-----------------------------------------------------------------
399 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
400 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
401 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
402 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
403 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
404
405 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
406
407 Float_t val = (Float_t)(( -xr*yr/TMath::Sqrt(xr*xr+yr*yr))
408 *( -xr*yr/TMath::Sqrt(xr*xr+yr*yr)));
409 return val;
410}