Improved cuts for the reconstruction of V0s (M.Ivanov)
[u/mrichter/AliRoot.git] / ITS / AliITSclusterTable.cxx
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  **************************************************************************/
15
16 //////////////////////////////////////////////////////////////////////////
17 // Class used to simplify some operations with clusters.                 //
18 // -Function FillArray fills an array wich contains, for each            //
19 //  ITS module, an array with the indices of all the clusters detected   //
20 //  by the module. The indices correspond to the cluster indices in class// 
21 //  AliITSlayer of AliITStrackerV2.                                      //
22 //  This function is used in AliITStrackerSA::FindTracks.                // 
23 // -Function FillArrayLabel fills an array wich contains, for each       //
24 //  particle label, and for each layer, the information on clusters:     //
25 //  0 if there is no cluster, 1 if there is a cluster with this label.   //
26 //  This function is used to define trackable tracks.                    //
27 // -Function FillArrayCoorAngles fills arrays wich contains, for each    //
28 //  layer, the global coordinates, errors on x,y,z and angles lambda     //
29 //  and phi for each cluster                                             //
30 ///////////////////////////////////////////////////////////////////////////
31
32
33 #include <stdlib.h>
34 #include "AliITSclusterTable.h"
35 #include "AliITSclusterV2.h"
36 #include "AliITSgeom.h"
37 #include "AliITStrackerSA.h"
38 #include<TBranch.h>
39 #include<TClonesArray.h>
40 #include<TTree.h>
41 ClassImp(AliITSclusterTable)
42
43 //__________________________________________________________
44 AliITSclusterTable::AliITSclusterTable(){
45 // Default constructor
46   fDet=0;
47   fNCl = 0;
48   fLbl=0;
49   fGeom = 0;
50   fTracker = 0;
51   fPhiList = 0; 
52   fLambdaList = 0;
53   fXList = 0;
54   fYList =0;
55   fZList =0;
56   fSxList = 0;
57   fSyList =0;
58   fSzList =0;
59   for(Int_t i=0;i<3;i++){fPrimaryVertex[i]=0;}
60 }
61
62 //__________________________________________________________
63 AliITSclusterTable::AliITSclusterTable(AliITSgeom* geom, AliITStrackerSA* tracker,Double_t* primaryVertex) {
64 // Standard constructor
65   fDet=0;
66   Int_t nm = geom->GetIndexMax();
67   fNCl = new Int_t[nm];
68   for(Int_t i=0;i<nm;i++){fNCl[i]=0;}
69   fLbl=0;
70   fGeom = geom;
71   fTracker = tracker;
72   for(Int_t i=0;i<3;i++){fPrimaryVertex[i]=primaryVertex[i];}
73   fPhiList = 0;
74   fLambdaList = 0;
75   fXList = 0;
76   fYList = 0;
77   fZList =0;
78   fSxList = 0;
79   fSyList =0;
80   fSzList =0;
81 }
82
83 //______________________________________________________________________
84 AliITSclusterTable::AliITSclusterTable(const AliITSclusterTable &tab) : 
85                     TObject(tab) {
86   // Copy constructor
87   // Copies are not allowed. The method is protected to avoid misuse.
88   Error("AliITSclusterTable","Copy constructor not allowed\n");
89 }
90
91 //______________________________________________________________________
92 AliITSclusterTable& AliITSclusterTable::operator=(const 
93                     AliITSclusterTable& /* tab */){
94   // Assignment operator
95   // Assignment is not allowed. The method is protected to avoid misuse.
96   Error("= operator","Assignment operator not allowed\n");
97   return *this;
98 }
99
100 //__________________________________________________________
101 AliITSclusterTable::~AliITSclusterTable(){
102 // Destructor
103   Int_t nm = fGeom->GetIndexMax();
104   if(fDet){
105     for (Int_t i=0; i<nm; i++){
106       delete fDet[i];
107     }
108     delete[] fDet;
109   }
110   if (fLbl)delete [] fLbl; // memory leak!
111   if (fNCl)delete [] fNCl;
112   for (Int_t i = 0; i < fGeom->GetNlayers(); i++) {
113     if (fPhiList) delete fPhiList[i];
114     if (fLambdaList) delete fLambdaList[i];
115     if (fXList) delete fXList[i];
116     if (fYList) delete fYList[i];
117     if (fZList) delete fZList[i];
118     if (fSxList) delete fSxList[i];
119     if (fSyList) delete fSyList[i];
120     if (fSzList) delete fSzList[i];
121   }
122   if (fPhiList) delete[] fPhiList;
123   if (fLambdaList) delete[] fLambdaList;
124   if (fXList) delete[] fXList;
125   if (fYList) delete[] fYList;
126   if (fZList) delete[] fZList;
127   if (fSxList) delete[] fSxList;
128   if (fSyList) delete[] fSyList;
129   if (fSzList) delete[] fSzList;
130 }
131
132 //__________________________________________________________
133
134 void AliITSclusterTable::FillArray(TTree* clusterTree){
135   
136   //
137   Int_t nm = fGeom->GetIndexMax();
138   fDet = new TArrayI*[nm];
139
140   TArrayI** vect = new TArrayI*[fGeom->GetNlayers()];
141   
142   Int_t * firstmod = new Int_t[fGeom->GetNlayers()+1];
143   firstmod[fGeom->GetNlayers()]=fGeom->GetIndexMax();  // upper limit
144   for(Int_t nlayer=0;nlayer<fGeom->GetNlayers();nlayer++){
145     firstmod[nlayer] = fGeom->GetModuleIndex(nlayer+1,1,1);
146     Int_t ncl = fTracker->GetNumberOfClustersLayer(nlayer);
147     vect[nlayer]=new TArrayI(ncl); 
148     for(Int_t j=0;j<ncl;j++){
149       AliITSclusterV2* cl = fTracker->GetClusterLayer(nlayer,j);
150       vect[nlayer]->AddAt(cl->GetDetectorIndex()+firstmod[nlayer],j);
151     }
152   }
153     
154   TBranch *brancht=(TBranch*)clusterTree->GetBranch("Clusters");
155   if(!brancht) Warning("FillArray","No cluster branch");
156   TClonesArray* clus = new TClonesArray("AliITSclusterV2",10000);
157   brancht->SetAddress(&clus);
158  
159   
160   for(Int_t mod=0;mod<nm;mod++){
161     Int_t nc=0;
162     clusterTree->GetEvent(mod);
163     Int_t ncl = clus->GetEntries();
164     fDet[mod] = new TArrayI(ncl);
165     fNCl[mod]= ncl;
166     Int_t nlr = FindIndex(fGeom->GetNlayers(),firstmod,mod);
167     if(nlr<0){
168       Fatal("FillArray","Wrong module number %d . Limits: %d , %d",mod,firstmod[0],firstmod[fGeom->GetNlayers()+1]);
169       exit(1);
170     }
171     else {
172       for(Int_t n=0;n<vect[nlr]->GetSize();n++){
173         Int_t mm=vect[nlr]->At(n);
174         if (nc>=fDet[mod]->GetSize()) fDet[mod]->Set(nc*2+10);  
175         if(mm==mod) {(*fDet[mod])[nc]=n; nc+=1; }
176       }
177     }
178   }
179
180
181   clus->Delete();
182   delete clus;
183   for(Int_t n=0;n<fGeom->GetNlayers();n++)delete vect[n];
184   delete [] vect;
185   delete [] firstmod;
186 }
187
188 //_________________________________________________________________
189 void AliITSclusterTable::FillArrayLabel(Int_t numberofparticles){
190   //
191
192
193   fLbl = new TArrayI*[numberofparticles];
194   const Int_t knm =fGeom->GetNlayers();
195   for(Int_t nlab=0;nlab<numberofparticles;nlab++){
196     fLbl[nlab] = new TArrayI(knm);
197     Int_t * nn = new Int_t[knm]; 
198     for(Int_t i=0;i<knm;i++)nn[i]=0;
199     for(Int_t nlayer=0;nlayer<knm;nlayer++){
200       Int_t ncl = fTracker->GetNumberOfClustersLayer(nlayer);
201       while(ncl--){
202         AliITSclusterV2* cl = fTracker->GetClusterLayer(nlayer,ncl);
203         if(cl->IsUsed()==1) continue;
204         if(cl->GetLabel(0)==nlab || cl->GetLabel(1)==nlab || cl->GetLabel(2)==nlab){
205           nn[nlayer]+=1;
206           cl->Use();
207           break;
208         }
209       }     
210       fLbl[nlab]->AddAt(nn[nlayer],nlayer);
211     }
212     delete [] nn;
213   }
214
215 }
216
217 //_______________________________________________________________
218 Int_t AliITSclusterTable::ThisParticleIsTrackable(Int_t label,Int_t numberofpoints){
219
220   //Returns 1 if particle with label "label" is trackable. 
221    
222   Int_t nb=0;
223   for(Int_t i=0;i<fGeom->GetNlayers();i++){
224     Int_t ncl = fLbl[label]->At(i);
225     if(ncl>0) nb++;
226   }
227   if(nb>=numberofpoints) return 1;
228   else return 0;
229
230 }
231
232 //_______________________________________________________________
233 Int_t AliITSclusterTable::FindIndex(Int_t ndim, Int_t *ptr, Int_t value){
234 // ptr[ndim+1] is an array of integers
235 // ndim is its dimension minus 1
236 // value is a number such as:  ptr[0]<=value < ptr[ndim]
237 // if ptr[i]<=value<ptr[i+1] the index i is returned;  -1 if out of bounds
238   Int_t retval = -1;
239   for(Int_t k=0;k<ndim;k++){
240     if(value>=ptr[k] && value <ptr[k+1]){
241       retval = k;
242       break;
243     }
244   }
245   return retval;
246 }
247
248 void AliITSclusterTable::FillArrayCoorAngles(){
249   //Fill arrays with phi,lambda and indices of clusters for each layer
250   Info("FillArrayCoorAngles","Filling Array...");
251
252   fPhiList = new TArrayD*[fGeom->GetNlayers()];
253   fLambdaList = new TArrayD*[fGeom->GetNlayers()];
254   fXList = new TArrayD*[fGeom->GetNlayers()];
255   fYList = new TArrayD*[fGeom->GetNlayers()];
256   fZList = new TArrayD*[fGeom->GetNlayers()];
257   fSxList = new TArrayD*[fGeom->GetNlayers()];
258   fSyList = new TArrayD*[fGeom->GetNlayers()];
259   fSzList = new TArrayD*[fGeom->GetNlayers()];
260
261   Int_t * firstmod = new Int_t[fGeom->GetNlayers()+1];
262   firstmod[fGeom->GetNlayers()]=fGeom->GetIndexMax();  // upper limit
263
264   for(Int_t nlay=0;nlay<fGeom->GetNlayers();nlay++){
265     firstmod[nlay] = fGeom->GetModuleIndex(nlay+1,1,1);
266     Int_t ncl = fTracker->GetNumberOfClustersLayer(nlay);
267     fPhiList[nlay] = new TArrayD(ncl);
268     fLambdaList[nlay]=new TArrayD(ncl);
269     fXList[nlay]=new TArrayD(ncl);
270     fYList[nlay]=new TArrayD(ncl);
271     fZList[nlay]=new TArrayD(ncl);
272     fSxList[nlay]=new TArrayD(ncl);
273     fSyList[nlay]=new TArrayD(ncl);
274     fSzList[nlay]=new TArrayD(ncl);
275
276     for(Int_t j=0;j<ncl;j++){
277       AliITSclusterV2* cl = fTracker->GetClusterLayer(nlay,j);
278       Double_t phi=0;Double_t lambda=0;
279       Double_t x=0;Double_t y=0;Double_t z=0;
280       Double_t sx=0;Double_t sy=0;Double_t sz=0;
281       Int_t module = cl->GetDetectorIndex()+firstmod[nlay];
282       GetCoorAngles(cl,module,phi,lambda,x,y,z);
283       GetCoorErrors(cl,module,sx,sy,sz);
284       fPhiList[nlay]->AddAt(phi,j);
285       fLambdaList[nlay]->AddAt(lambda,j);
286       fXList[nlay]->AddAt(x,j);
287       fYList[nlay]->AddAt(y,j);
288       fZList[nlay]->AddAt(z,j);
289       fSxList[nlay]->AddAt(sx,j);
290       fSyList[nlay]->AddAt(sy,j);
291       fSzList[nlay]->AddAt(sz,j);
292       
293     }
294     
295   }
296   
297   delete [] firstmod;
298 }
299
300 void AliITSclusterTable::GetCoorAngles(AliITSclusterV2* cl,Int_t module,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z){
301   //Returns values of phi (azimuthal) and lambda angles for a given cluster
302   
303   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
304   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
305   Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);     
306
307   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
308   Double_t phi1=TMath::Pi()/2+alpha;
309   if (lay==1) phi1+=TMath::Pi();
310
311   Double_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
312   Double_t r=tx*cp+ty*sp;
313
314   x= r*cp - cl->GetY()*sp;
315   y= r*sp + cl->GetY()*cp;
316   z=cl->GetZ();
317   
318   phi=TMath::ATan2(y,x);
319   lambda=TMath::ATan2(z-fPrimaryVertex[2],TMath::Sqrt((x-fPrimaryVertex[0])*(x-fPrimaryVertex[0])+(y-fPrimaryVertex[1])*(y-fPrimaryVertex[1])));
320 }
321
322 void AliITSclusterTable::GetCoorErrors(AliITSclusterV2* cl, Int_t module,Double_t &sx,Double_t &sy, Double_t &sz){
323
324   //returns x,y,z of cluster in global coordinates
325
326   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
327   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
328  
329   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
330   Double_t phi=TMath::Pi()/2+alpha;
331   if (lay==1) phi+=TMath::Pi();
332
333   Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
334
335  
336   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
337   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
338   sz = TMath::Sqrt(cl->GetSigmaZ2());
339
340 }