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