Fixes for Coverity warnings (M. van Leeuwen)
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerCosmics.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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  **************************************************************************/
15
16 #include <TClonesArray.h>
17 #include <TTree.h>
18 #include "AliLog.h"
19 #include "AliESDVertex.h"
20 #include "AliITSgeomTGeo.h"
21 #include "AliITSRecPoint.h"
22 #include "AliITSReconstructor.h"
23 #include "AliITSVertexerCosmics.h"
24 #include "AliStrLine.h"
25
26 //------------------------------------------------------------------------
27 // This class implements a method to construct a "fake" primary
28 // vertex for cosmic events in which the muon crosses one of 5 inner
29 // ITS layers. A fake primary vertex is needed for the reconstruction,
30 // with e.g. AliITStrackerSA, of the two tracks produced by the muon 
31 // in the ITS.
32 //   We build pairs of clusters on a given layer and define the fake vertex as
33 // the mid-point of the straight line joining the two clusters.
34 //   We use the innermost layer that has at least two clusters.
35 //   We reject the background by requiring at least one cluster on the outer
36 // layer, closer than fMaxDistOnOuterLayer to the tracklet prolongation.
37 //   We can reject (potentially pathological) events with the muon track
38 // tangential to the layer by the requiring the radial position of
39 // the vertex to be smaller than fMaxVtxRadius.
40 //   Due to background clusters, more than one vertex per event can 
41 // be found. We consider the first found.
42 //   The errors on x,y,z of the vertex are calculated as errors on the mean
43 // of clusters coordinates. Non-diag elements of vertex cov. mat. are set to 0.
44 //   The number of contributors set in the AliESDVertex object is the
45 // number of the layer on which the tracklet was built; if this number is -1, 
46 // the procedure could not find a vertex position and by default 
47 // the vertex coordinates are set to (0,0,0) with large errors (100,100,100)
48 //
49 // Origin: A.Dainese, andrea.dainese@lnl.infn.it
50 //-------------------------------------------------------------------------
51
52 ClassImp(AliITSVertexerCosmics)
53
54 //-------------------------------------------------------------------------
55 AliITSVertexerCosmics::AliITSVertexerCosmics():AliITSVertexer(),
56 fMaxDistOnOuterLayer(0),
57 fMinDist2Vtxs(0)
58 {
59   // Default constructor
60   SetFirstLastModules(0,0,79);
61   SetFirstLastModules(1,80,239);
62   SetFirstLastModules(2,240,323);
63   SetFirstLastModules(3,324,499);
64   SetFirstLastModules(4,500,1247);
65   SetFirstLastModules(5,1248,2197);
66   /*
67   SetMaxVtxRadius(0,3.5);
68   SetMaxVtxRadius(1,6.5);
69   SetMaxVtxRadius(2,14.5);
70   SetMaxVtxRadius(3,23.5);
71   SetMaxVtxRadius(4,37.5);
72   SetMaxVtxRadius(5,42.5);
73   */  
74   SetMaxVtxRadius(0,5.5);
75   SetMaxVtxRadius(1,8.5);
76   SetMaxVtxRadius(2,18.5);
77   SetMaxVtxRadius(3,28.5);
78   SetMaxVtxRadius(4,39.5);
79   SetMaxVtxRadius(5,48.5);
80   
81   SetMaxDistOnOuterLayer();
82   SetMinDist2Vtxs();
83 }
84 //--------------------------------------------------------------------------
85 AliESDVertex* AliITSVertexerCosmics::FindVertexForCurrentEvent(TTree *itsClusterTree) 
86 {
87   // Defines the AliESDVertex for the current event
88
89   fCurrentVertex = 0;
90
91   TClonesArray *recpoints=new TClonesArray("AliITSRecPoint",10000);
92   itsClusterTree->SetBranchAddress("ITSRecPoints",&recpoints);
93
94   Int_t lay,lad,det; 
95
96   Double_t pos[3]={0.,0.,0.};
97   Double_t err[3]={100.,100.,100.};
98   Int_t ncontributors = -1;
99
100   // Search for innermost layer with at least two clusters 
101   // on two different modules
102   Int_t ilayer=0,ilayer2=0;
103   Int_t nHitModulesSPDinner=0;
104   while(ilayer<AliITSgeomTGeo::GetNLayers()) {
105     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) {
106       ilayer++;
107       continue;
108     }
109     Int_t nHitModules=0;
110     for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer]; imodule++) {
111       itsClusterTree->GetEvent(imodule);
112       AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
113       lay -= 1;  // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
114       if(lay!=ilayer) AliFatal("Layer mismatch!");
115       if(recpoints->GetEntriesFast()>0) nHitModules++;
116     }
117     if(ilayer==0) nHitModulesSPDinner=nHitModules;
118     if(nHitModules>=2) break;
119     ilayer++;
120   }
121
122   ilayer2=ilayer+1;
123   while(ilayer2<6) {
124     if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer2)) break;
125     ilayer2++;
126   }
127
128   // try tracklet on SPD2 and point on SPD1
129   if(ilayer==1 && !AliITSReconstructor::GetRecoParam()->GetLayersToSkip(0) &&
130      nHitModulesSPDinner>0) { ilayer=0; ilayer2=1; }
131
132   if(ilayer>4 || ilayer2>5) {
133     AliWarning("Not enough clusters");
134     delete recpoints;
135     fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
136     fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
137     fCurrentVertex->SetNContributors(ncontributors);
138     return fCurrentVertex;
139   }
140
141
142   const Int_t arrSize = 1000;
143   Float_t xclInnLay[arrSize],yclInnLay[arrSize],zclInnLay[arrSize],modclInnLay[arrSize];
144   Float_t e2xclInnLay[arrSize],e2yclInnLay[arrSize],e2zclInnLay[arrSize];
145   Float_t e2xclOutLay[arrSize],e2yclOutLay[arrSize],e2zclOutLay[arrSize];
146   Int_t nclInnLayStored=0;
147   Float_t xclOutLay[arrSize],yclOutLay[arrSize],zclOutLay[arrSize],modclOutLay[arrSize];
148   Int_t nclOutLayStored=0;
149   Int_t nRecPoints,nRecPointsInnLay=0;
150
151   Float_t gc[3],gcov[6];
152
153   Float_t matchOutLayValue;
154   Float_t distxyInnLay,distxyInnLayBest=0.;
155   Double_t p1[3],p2[3],p3[3];
156
157   Float_t xvtx,yvtx,zvtx,rvtx;
158   Int_t i1InnLayBest=-1,i2InnLayBest=-1;
159
160
161   // Collect clusters in the selected layer and the outer one
162   for(Int_t imodule=fFirst[ilayer]; imodule<=fLast[ilayer2]; imodule++) {
163     itsClusterTree->GetEvent(imodule);
164     AliITSgeomTGeo::GetModuleId(imodule,lay,lad,det);
165     lay -= 1; // AliITSgeomTGeo gives layer from 1 to 6, we want 0 to 5
166     nRecPoints=recpoints->GetEntriesFast();
167     if(imodule<=fLast[ilayer]) nRecPointsInnLay += nRecPoints;
168     //printf("cosmics: module %d clusters %d\n",imodule,nRecPoints);
169     for(Int_t irp=0; irp<nRecPoints; irp++) {
170       AliITSRecPoint *rp=(AliITSRecPoint*)recpoints->UncheckedAt(irp);
171       // Local coordinates of this recpoint
172       rp->GetGlobalXYZ(gc);
173       if(lay==ilayer) { // store InnLay clusters
174         xclInnLay[nclInnLayStored]=gc[0];
175         yclInnLay[nclInnLayStored]=gc[1];
176         zclInnLay[nclInnLayStored]=gc[2];
177         rp->GetGlobalCov(gcov);
178         e2xclInnLay[nclInnLayStored]=gcov[0];
179         e2yclInnLay[nclInnLayStored]=gcov[3];
180         e2zclInnLay[nclInnLayStored]=gcov[5];
181         modclInnLay[nclInnLayStored]=imodule;
182         nclInnLayStored++;
183       }
184       if(lay==ilayer2) { // store OutLay clusters
185         xclOutLay[nclOutLayStored]=gc[0];
186         yclOutLay[nclOutLayStored]=gc[1];
187         zclOutLay[nclOutLayStored]=gc[2];
188         rp->GetGlobalCov(gcov);
189         e2xclOutLay[nclOutLayStored]=gcov[0];
190         e2yclOutLay[nclOutLayStored]=gcov[3];
191         e2zclOutLay[nclOutLayStored]=gcov[5];
192         modclOutLay[nclOutLayStored]=imodule;
193         nclOutLayStored++;
194       }
195       if(nclInnLayStored>arrSize || nclOutLayStored>arrSize) {
196         //AliFatal("More than arrSize clusters per layer");
197         AliWarning("Too many clusters per layer");
198         delete recpoints;
199         fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
200         fCurrentVertex->SetTitle("cosmics fake vertex (failed)");
201         fCurrentVertex->SetNContributors(ncontributors);
202         return fCurrentVertex;
203       }
204     }// end clusters in a module
205   }// end modules
206
207
208   Int_t i1InnLay,i2InnLay,iOutLay;
209
210   // build fake vertices
211   //printf("Building tracklets on layer %d\n",ilayer);
212
213   // InnLay - first cluster
214   for(i1InnLay=0; i1InnLay<nclInnLayStored; i1InnLay++) { 
215     p1[0]=xclInnLay[i1InnLay]; 
216     p1[1]=yclInnLay[i1InnLay]; 
217     p1[2]=zclInnLay[i1InnLay];
218     // InnLay - second cluster
219     for(i2InnLay=i1InnLay+1; i2InnLay<nclInnLayStored; i2InnLay++) { 
220       if(modclInnLay[i1InnLay]==modclInnLay[i2InnLay]) continue;
221       p2[0]=xclInnLay[i2InnLay]; 
222       p2[1]=yclInnLay[i2InnLay]; 
223       p2[2]=zclInnLay[i2InnLay];
224       // look for point on OutLay
225       AliStrLine innLayline(p1,p2,kTRUE);
226       for(iOutLay=0; iOutLay<nclOutLayStored; iOutLay++) {
227         p3[0]=xclOutLay[iOutLay]; 
228         p3[1]=yclOutLay[iOutLay]; 
229         p3[2]=zclOutLay[iOutLay];
230         //printf("(%f,%f) (%f,%f)     (%f,%f) %f\n",p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],innLayline.GetDistFromPoint(p3));
231         matchOutLayValue=innLayline.GetDistFromPoint(p3);
232         distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
233         if(matchOutLayValue<fMaxDistOnOuterLayer &&
234            distxyInnLay>distxyInnLayBest) { 
235           //printf("found\n");
236           distxyInnLayBest=distxyInnLay;
237           i1InnLayBest=i1InnLay;
238           i2InnLayBest=i2InnLay;
239         }
240       }
241     } // InnLay - second cluster
242   } // InnLay - first cluster
243
244   if(i1InnLayBest>-1 && i2InnLayBest>-1) { 
245     xvtx = 0.5*(xclInnLay[i1InnLayBest]+xclInnLay[i2InnLayBest]);
246     yvtx = 0.5*(yclInnLay[i1InnLayBest]+yclInnLay[i2InnLayBest]);
247     zvtx = 0.5*(zclInnLay[i1InnLayBest]+zclInnLay[i2InnLayBest]);
248     rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
249     if(rvtx<fMaxVtxRadius[ilayer]) {
250       ncontributors = ilayer;
251       pos[0] = xvtx;
252       pos[1] = yvtx;
253       pos[2] = zvtx;
254       err[0]=TMath::Sqrt(0.25*(e2xclInnLay[i1InnLayBest]+e2xclInnLay[i2InnLayBest])); 
255       err[1]=TMath::Sqrt(0.25*(e2yclInnLay[i1InnLayBest]+e2yclInnLay[i2InnLayBest])); 
256       err[2]=TMath::Sqrt(0.25*(e2zclInnLay[i1InnLayBest]+e2zclInnLay[i2InnLayBest]));
257     }
258
259   } else { // give it a try exchanging InnLay and OutLay
260
261     // OutLay - first cluster
262     for(i1InnLay=0; i1InnLay<nclOutLayStored; i1InnLay++) { 
263       p1[0]=xclOutLay[i1InnLay]; 
264       p1[1]=yclOutLay[i1InnLay]; 
265       p1[2]=zclOutLay[i1InnLay];
266       // OutLay - second cluster
267       for(i2InnLay=i1InnLay+1; i2InnLay<nclOutLayStored; i2InnLay++) { 
268         if(modclOutLay[i1InnLay]==modclOutLay[i2InnLay]) continue;
269         p2[0]=xclOutLay[i2InnLay]; 
270         p2[1]=yclOutLay[i2InnLay]; 
271         p2[2]=zclOutLay[i2InnLay];
272         // look for point on InnLay
273         AliStrLine outLayline(p1,p2,kTRUE);
274         for(iOutLay=0; iOutLay<nclInnLayStored; iOutLay++) {
275           p3[0]=xclInnLay[iOutLay]; 
276           p3[1]=yclInnLay[iOutLay]; 
277           p3[2]=zclInnLay[iOutLay];
278           //printf(" %f\n",InnLayline.GetDistFromPoint(p3));
279           matchOutLayValue=outLayline.GetDistFromPoint(p3);
280           distxyInnLay = (p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1]);
281           if(matchOutLayValue<fMaxDistOnOuterLayer &&
282              distxyInnLay>distxyInnLayBest) { 
283             distxyInnLayBest=distxyInnLay;
284             i1InnLayBest=i1InnLay;
285             i2InnLayBest=i2InnLay;
286           }
287         }
288       } // OutLay - second cluster
289     } // OutLay - first cluster
290
291     if(i1InnLayBest>-1 && i2InnLayBest>-1) { 
292       xvtx = 0.5*(xclOutLay[i1InnLayBest]+xclOutLay[i2InnLayBest]);
293       yvtx = 0.5*(yclOutLay[i1InnLayBest]+yclOutLay[i2InnLayBest]);
294       zvtx = 0.5*(zclOutLay[i1InnLayBest]+zclOutLay[i2InnLayBest]);
295       rvtx = TMath::Sqrt(xvtx*xvtx+yvtx*yvtx);
296       if(rvtx<fMaxVtxRadius[ilayer]) {
297         ncontributors = ilayer2;
298         pos[0] = xvtx;
299         pos[1] = yvtx;
300         pos[2] = zvtx;
301         err[0]=TMath::Sqrt(0.25*(e2xclOutLay[i1InnLayBest]+e2xclOutLay[i2InnLayBest])); 
302         err[1]=TMath::Sqrt(0.25*(e2yclOutLay[i1InnLayBest]+e2yclOutLay[i2InnLayBest])); 
303         err[2]=TMath::Sqrt(0.25*(e2zclOutLay[i1InnLayBest]+e2zclOutLay[i2InnLayBest]));
304       }
305     }
306   } // give it a try exchanging InnLay and OutLay
307
308   fCurrentVertex = new AliESDVertex(pos,err,"cosmics");
309   fCurrentVertex->SetTitle("cosmics fake vertex");
310   fCurrentVertex->SetNContributors(ncontributors);
311   //fCurrentVertex->Print();
312   if(fComputeMultiplicity) FindMultiplicity(itsClusterTree);
313
314   delete recpoints;
315
316   return fCurrentVertex;
317 }  
318
319 //-------------------------------------------------------------------------
320 void AliITSVertexerCosmics::PrintStatus() const 
321 {
322   // Print current status
323   printf("=======================================================\n");
324   printf(" fMaxDistOnOuterLayer: %f\n",fMaxDistOnOuterLayer);
325   printf(" fMaxVtxRadius[0]:  %f\n",fMaxVtxRadius[0]);
326   printf(" fMinDist2Vtxs:  %f\n",fMinDist2Vtxs);
327   printf("=======================================================\n");
328 }
329 //-------------------------------------------------------------------------