new SDD preprocessor + removal of eff C++ warning (base) - E. Crescio
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinder.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 // Base Class used to find                                                //
18 // the reconstructed points for ITS                                       //
19 // See also AliITSClusterFinderSPD, AliITSClusterFinderSDD,               //
20 // AliITSClusterFinderSDD  AliITSClusterFinderV2                          //
21 ////////////////////////////////////////////////////////////////////////////
22
23 #include "AliITSClusterFinder.h"
24 #include "AliITSRecPoint.h"
25 #include "AliITSdigit.h"
26 #include "AliITSDetTypeRec.h"
27 #include "AliITSMap.h"
28
29 ClassImp(AliITSClusterFinder)
30
31 //----------------------------------------------------------------------
32 AliITSClusterFinder::AliITSClusterFinder():
33 TObject(),
34 fDebug(0),
35 fModule(0),
36 fDigits(0),
37 fNdigits(0),
38 fDetTypeRec(0),
39 fClusters(0),
40 fNRawClusters(0),
41 fMap(0),
42 fNperMax(0),
43 fDeclusterFlag(0),
44 fClusterSize(0),
45 fNPeaks(-1){
46     // default cluster finder
47     // Input:
48     //   none.
49     // Output:
50     //   none.
51     // Return:
52     //   A default constructed AliITSCulsterFinder
53 }
54 //----------------------------------------------------------------------
55 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp):
56 TObject(),
57 fDebug(0),
58 fModule(0),
59 fDigits(0),
60 fNdigits(0),
61 fDetTypeRec(dettyp),
62 fClusters(0),
63 fNRawClusters(0),
64 fMap(0),
65 fNperMax(0),
66 fDeclusterFlag(0),
67 fClusterSize(0),
68 fNPeaks(-1){
69     // Standard constructor for cluster finder
70     // Input:
71     //   AliITSsegmentation *seg  The segmentation class to be used
72     //   AliITSresponse     *res  The response class to be used
73     // Output:
74     //   none.
75     // Return:
76     //   A Standard constructed AliITSCulsterFinder
77
78     SetNperMax();
79     SetClusterSize();
80     SetDeclusterFlag();
81 }
82 //----------------------------------------------------------------------
83 AliITSClusterFinder::AliITSClusterFinder(AliITSDetTypeRec* dettyp,
84                                          TClonesArray *digits):
85 TObject(),
86 fDebug(0),
87 fModule(0),
88 fDigits(digits),
89 fNdigits(0),
90 fDetTypeRec(dettyp),
91 fClusters(0),
92 fNRawClusters(0),
93 fMap(0),
94 fNperMax(0),
95 fDeclusterFlag(0),
96 fClusterSize(0),
97 fNPeaks(-1){
98     // Standard + cluster finder constructor
99     // Input:
100     //   AliITSsegmentation *seg  The segmentation class to be used
101     //   AliITSresponse     *res  The response class to be used
102     //   TClonesArray    *digits  Array of digits to be used
103     // Output:
104     //   none.
105     // Return:
106     //   A Standard constructed AliITSCulsterFinder
107
108     fNdigits = fDigits->GetEntriesFast();
109     SetNperMax();
110     SetClusterSize();
111     SetDeclusterFlag();
112 }
113 /*
114 //______________________________________________________________________
115 AliITSClusterFinder::AliITSClusterFinder(const AliITSClusterFinder &source) : TObject(source) {
116   // Copy constructor
117   // Copies are not allowed. The method is protected to avoid misuse.
118   Fatal("AliITSClusterFinder","Copy constructor not allowed\n");
119 }
120 */
121
122 //______________________________________________________________________
123 //AliITSClusterFinder& AliITSClusterFinder::operator=(const AliITSClusterFinder& /* source */){
124   // Assignment operator
125   // Assignment is not allowed. The method is protected to avoid misuse.
126 //  Fatal("= operator","Assignment operator not allowed\n");
127 //  return *this;
128 //}
129
130 //----------------------------------------------------------------------
131 AliITSClusterFinder::~AliITSClusterFinder(){
132     // destructor cluster finder
133     // Input:
134     //   none.
135     // Output:
136     //   none.
137     // Return:
138     //   none.
139
140     if(fMap) {delete fMap;}
141     // Zero local pointers. Other classes own these pointers.
142     fMap          = 0;
143     fDigits       = 0;
144     fNdigits      = 0;
145     fNRawClusters = 0;
146     fNperMax      = 0;
147     fDeclusterFlag= 0;
148     fClusterSize  = 0;
149     fNPeaks       = 0;
150     fDetTypeRec   = 0;
151
152 }
153 //__________________________________________________________________________
154 void AliITSClusterFinder::InitGeometry(){
155
156
157  //Initialisation of ITS geometry
158   if(!fDetTypeRec->GetITSgeom()) {
159     Error("InitGeometry","ITS geom is null!");
160     return;
161   }
162   Int_t mmax=fDetTypeRec->GetITSgeom()->GetIndexMax();
163   if (mmax>2200) {
164     Fatal("AliITSClusterFinder","Too many ITS subdetectors !"); 
165   }
166   Int_t m;
167   for (m=0; m<mmax; m++) {
168     Int_t lay,lad,det; fDetTypeRec->GetITSgeom()->GetModuleId(m,lay,lad,det);
169     Float_t x,y,z;     fDetTypeRec->GetITSgeom()->GetTrans(lay,lad,det,x,y,z); 
170     Double_t rot[9];   fDetTypeRec->GetITSgeom()->GetRotMatrix(lay,lad,det,rot);
171     Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
172     Double_t ca=TMath::Cos(alpha), sa=TMath::Sin(alpha);
173     fYshift[m] = x*ca + y*sa;
174     fZshift[m] = (Double_t)z;
175     fNdet[m] = (lad-1)*fDetTypeRec->GetITSgeom()->GetNdetectors(lay) + (det-1);
176     fNlayer[m] = lay-1;
177   }
178 }
179
180
181
182 //----------------------------------------------------------------------
183 void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c){
184     // Add a raw cluster copy to the list
185     // Input:
186     //   Int_t       branch  The branch to which the cluster is to be added to
187     //   AliITSRawCluster *c The cluster to be added to the array of clusters
188     // Output:
189     //   none.
190     // Return:
191     //   none.
192
193    if(!fDetTypeRec) {
194     Error("AddCluster","fDetTypeRec is null!");
195     return;
196   }
197   fDetTypeRec->AddCluster(branch,c); 
198   fNRawClusters++;
199 }
200 //----------------------------------------------------------------------
201 void AliITSClusterFinder::AddCluster(Int_t branch, AliITSRawCluster *c, 
202                                      AliITSRecPoint &rp){
203     // Add a raw cluster copy to the list and the RecPoint
204     // Input:
205     //   Int_t       branch  The branch to which the cluster is to be added to
206     //   AliITSRawCluster *c The cluster to be added to the array of clusters
207     //   AliITSRecPoint  &rp The RecPoint to be added to the array of RecPoints
208     // Output:
209     //   none.
210     // Return:
211     //   none.
212   if(!fDetTypeRec) {
213     Error("AddCluster","fDetTypeRec is null!");
214     return;
215   }
216
217   fDetTypeRec->AddCluster(branch,c); 
218   fNRawClusters++;
219   fDetTypeRec->AddRecPoint(rp); 
220
221 }
222 /*
223 //______________________________________________________________________
224 void AliITSClusterFinder::CheckLabels(Int_t lab[3]) {
225   //------------------------------------------------------------
226   // Tries to find mother's labels
227   //------------------------------------------------------------
228
229   if(lab[0]<0 && lab[1]<0 && lab[2]<0) return; // In case of no labels just exit
230   // Check if simulation
231   AliMC* mc = gAlice->GetMCApp();
232   if(!mc)return;
233
234   Int_t ntracks = mc->GetNtrack();
235   for (Int_t i=0;i<3;i++){
236     Int_t label = lab[i];
237     if (label>=0 && label<ntracks) {
238       TParticle *part=(TParticle*)mc->Particle(label);
239       if (part->P() < 0.005) {
240         Int_t m=part->GetFirstMother();
241         if (m<0) {      
242           continue;
243         }
244         if (part->GetStatusCode()>0) {
245           continue;
246         }
247         lab[i]=m;       
248       }
249     }    
250   }
251   
252 }
253 */
254 //______________________________________________________________________
255 void AliITSClusterFinder::FindRawClusters(Int_t module){
256     // Default Cluster finder.
257     // Input:
258     //   Int_t module   Module number for which culster are to be found.
259     // Output:
260     //   none.
261     // Return:
262     //   none.
263     const Int_t kelms = 10;
264     Int_t ndigits = fDigits->GetEntriesFast();
265     TObjArray *digs = new TObjArray(ndigits);
266     TObjArray *clusts = new TObjArray(ndigits); // max # cluster
267     TObjArray *clust0=0; // A spacific cluster of digits
268     TObjArray *clust1=0; // A spacific cluster of digits
269     AliITSdigit *dig=0; // locat pointer to a digit
270     Int_t i=0,nc=0,j[4],k,k2=0;
271
272     // Copy all digits for this module into a local TObjArray.
273     for(i=0;i<ndigits;i++) digs->AddAt(new AliITSdigit(*(GetDigit(i))),i);
274     digs->Sort();
275     // First digit is a cluster.
276     i  = 0;
277     nc = 0;
278     clusts->AddAt(new TObjArray(kelms),nc);
279     clust0 = (TObjArray*)(clusts->At(nc));
280     clust0->AddAtFree(digs->At(i)); // move owner ship from digs to clusts
281     nc++;
282     for(i=1;i<ndigits;i++){
283         if(IsNeighbor(digs,i,j)){
284             dig = (AliITSdigit*)(digs->At(j[0]));
285             // Add to existing cluster. Find which cluster this digis 
286             for(k=0;k<nc;k++){
287                 clust0 = ((TObjArray*)(clusts->At(k)));
288                 if(clust0->IndexOf(dig)>=0) break;
289             } // end for k
290             if(k>=nc){
291                 Fatal("FindRawClusters","Digit not found as expected");
292             } // end if
293             if(j[1]>=0){
294                 dig = (AliITSdigit*)(digs->At(j[1]));
295                 // Add to existing cluster. Find which cluster this digis 
296                 for(k2=0;k2<nc;k2++){
297                     clust1 = ((TObjArray*)(clusts->At(k2)));
298                     if(clust1->IndexOf(dig)>=0) break;
299                 } // end for k2
300                 if(k2>=nc){
301                     Fatal("FindRawClusters","Digit not found as expected");
302                 } // end if
303             } // end if j[1]>=0
304             // Found cluster with neighboring digits add this one to it.
305             if(clust0==clust1){ // same cluster
306                 clust0->AddAtFree(digs->At(i));
307                 clust0 = 0; // finished with cluster. zero for safty
308                 clust1 = 0; // finished wit hcluster. zero for safty
309             }else{ // two different clusters which need to be merged.
310                 clust0->AddAtFree(digs->At(i)); // Add digit to this cluster.
311                 for(k=0;k<clust1->GetEntriesFast();k++){
312                     // move clust1 into clust0
313                     //move digit to this cluster
314                     clust0->AddAtFree(clust1->At(k));
315                     clust1->AddAt(0,k); // zero this one
316                 } // end for k
317                 delete clust1;
318                 clusts->AddAt(0,k2); // zero array of clusters element clust1
319                 clust0 = 0; // finished with cluster. zero for safty
320                 clust1 = 0; // finished wit hcluster. zero for safty
321             } // end if clust0==clust1
322         }else{// New cluster
323             clusts->AddAt(new TObjArray(kelms),nc);
324             clust0 = ((TObjArray*)(clusts->At(nc)));
325             // move owner ship from digs to clusts
326             clust0->AddAtFree(digs->At(i));
327             clust0 = 0; // finished with cluster. zero for safty
328             nc++;
329         } // End if IsNeighbor
330     } // end for i
331     // There are now nc clusters in clusts. Each element of clust is an
332     // array of digits which are clustered together.
333
334     // For each cluster call detector specific CreateRecPoints
335     for(i=0;i<nc;i++) CreateRecPoints((TObjArray*)(clusts->At(i)),module);
336
337     // clean up at the end.
338     for(i=0;i<nc;i++){ 
339         clust0 =(TObjArray*)(clusts->At(i));
340         // Digits deleted below, so zero this TObjArray
341         for(k=0;k<clust0->GetEntriesFast();k++) clust0->AddAt(0,k);
342         delete clust0; // Delete this TObjArray
343         clusts->AddAt(0,i); // Contents deleted above, so zero it.
344     } // end for i
345     delete clusts; // Delete this TObjArray/
346     // Delete the digits then the TObjArray which containted them.
347     for(i=0;i<ndigits;i++) delete ((AliITSdigit*)(digs->At(i)));
348     delete digs;
349 }
350 //______________________________________________________________________
351 Bool_t AliITSClusterFinder::IsNeighbor(TObjArray *digs,Int_t i,Int_t n[])const{
352     // Locagical function which checks to see if digit i has a neighbor.
353     // If so, then it returns kTRUE and its neighbor index j.
354     // This routine checks if the digits are side by side or one before the
355     // other. Requires that the array of digits be in proper order.
356     // Returns kTRUE in the following cases.
357     //                 ji   0j   if kdiagonal  j0    0i
358     //                 00   0i   if kdiagonal  0i    j0
359     // Inputs:
360     //    TObjArray *digs   Array to search for neighbors in
361     //    Int_t      i      Index of digit for which we are searching for
362     //                      a neighbor of.
363     // Output:
364     //    Int_t      j[4]   Index of one or more of the digits which is a
365     //                      neighbor of digit a index i.
366     // Return:
367     //    Bool_t            kTRUE if a neighbor was found kFALSE otherwise.
368     Int_t ix,jx,iz,jz,j;
369     const Bool_t kdiagonal=kFALSE;
370     Bool_t nei[4];
371
372     // No neighbors found if array empty.
373     if(digs->GetEntriesFast()<=0) return kFALSE;
374     // can not be a digit with first element or elements out or range
375     if(i<=0 || i>=digs->GetEntriesFast()) return kFALSE;
376
377     for(j=0;j<4;j++){n[j] = -1;nei[j]=kFALSE;}
378     ix = ((AliITSdigit*)(digs->At(i)))->GetCoord1();
379     iz = ((AliITSdigit*)(digs->At(i)))->GetCoord2();
380     for(j=0;j<i;j++){
381         jx = ((AliITSdigit*)(digs->At(j)))->GetCoord1();
382         jz = ((AliITSdigit*)(digs->At(j)))->GetCoord2();
383         if(jx+1==ix && jz  ==iz){n[0] = j;nei[0] = kTRUE;}
384         if(jx  ==ix && jz+1==iz){n[1] = j;nei[1] = kTRUE;}
385         if(jx+1==ix && jz+1==iz){n[2] = j;nei[2] = kTRUE;}
386         if(jx+1==ix && jz-1==iz){n[3] = j;nei[3] = kTRUE;}
387     } // end for k
388     if(nei[0]||nei[1]) return kTRUE;
389     if(kdiagonal&&(nei[2]||nei[3])) return kTRUE;
390     // no Neighbors found.
391     return kFALSE;
392 }
393
394 //______________________________________________________________________
395 void AliITSClusterFinder::Print(ostream *os) const{
396     //Standard output format for this class
397     // Inputs:
398     //    ostream *os   Output stream
399     // Output:
400     //    ostream *os   Output stream
401     // Return:
402     //    none.
403
404     *os << fDebug<<",";
405     *os << fModule<<",";
406     *os << fNdigits<<",";
407     *os << fNRawClusters<<",";
408     *os << fNperMax<<",";
409     *os << fDeclusterFlag<<",";
410     *os << fClusterSize<<",";
411     *os << fNPeaks<<endl;
412 }
413 /*
414 //______________________________________________________________________
415 void AliITSClusterFinder::RecPoints2Clusters
416 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
417   //------------------------------------------------------------
418   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS 
419   // subdetector indexed by idx 
420   //------------------------------------------------------------
421   if(!fDetTypeRec->GetITSgeom()) {
422     Error("RecPoints2Clusters","ITS geom is null!");
423     return;
424   }
425   Int_t lastSPD1=fDetTypeRec->GetITSgeom()->GetModuleIndex(2,1,1)-1;
426   TClonesArray &cl=*clusters;
427   Int_t ncl=points->GetEntriesFast();
428   for (Int_t i=0; i<ncl; i++) {
429     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
430     Float_t lp[5];
431     lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=lastSPD1) lp[0]*=-1; //SPD1
432     lp[1]=  -p->GetZ()+fZshift[idx];
433     lp[2]=p->GetSigmaX2();
434     lp[3]=p->GetSigmaZ2();
435     lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
436     Int_t lab[4]; 
437     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
438     lab[3]=fNdet[idx];
439     CheckLabels(lab);
440     Int_t dummy[3]={0,0,0};
441     new (cl[i]) AliITSclusterV2(lab,lp, dummy);
442   }  
443
444 */
445 //______________________________________________________________________
446 void AliITSClusterFinder::Read(istream *is)  {
447     //Standard input for this class
448     // Inputs:
449     //    istream *is   Input stream
450     // Output:
451     //    istream *is   Input stream
452     // Return:
453     //    none.
454
455     *is >> fDebug;
456     *is >> fModule;
457     *is >> fNdigits;
458     *is >> fNRawClusters;
459     *is >> fNperMax;
460     *is >> fDeclusterFlag;
461     *is >> fClusterSize;
462     *is >> fNPeaks;
463 }
464 //______________________________________________________________________
465 ostream &operator<<(ostream &os,AliITSClusterFinder &source){
466     // Standard output streaming function.
467     // Inputs:
468     //    ostream             *os     Output stream
469     //    AliITSClusterFinder &source Class to be printed
470     // Output:
471     //    ostream             *os     Output stream
472     // Return:
473     //    none.
474
475     source.Print(&os);
476     return os;
477 }
478 //______________________________________________________________________
479 istream &operator>>(istream &is,AliITSClusterFinder &source){
480     // Standard output streaming function.
481     // Inputs:
482     //    istream              *is      Input stream
483     //     AliITSClusterFinder &source  Class to be read in.
484     // Output:
485     //    istream              *is      Input stream
486     // Return:
487     //    none.
488
489     source.Read(&is);
490     return is;
491 }
492 /*
493 void AliITSClusterFinder::RecPoints2Clusters
494 (const TClonesArray *points, Int_t idx, TClonesArray *clusters) {
495   //------------------------------------------------------------
496   // Conversion AliITSRecPoint -> AliITSclusterV2 for the ITS 
497   // subdetector indexed by idx 
498   //------------------------------------------------------------
499   TClonesArray &cl=*clusters;
500   Int_t ncl=points->GetEntriesFast();
501   for (Int_t i=0; i<ncl; i++) {
502     AliITSRecPoint *p = (AliITSRecPoint *)points->UncheckedAt(i);
503     Float_t lp[5];
504     lp[0]=-(-p->GetX()+fYshift[idx]); if (idx<=fLastSPD1) lp[0]*=-1; //SPD1
505     lp[1]=  -p->GetZ()+fZshift[idx];
506     lp[2]=p->GetSigmaX2();
507     lp[3]=p->GetSigmaZ2();
508     lp[4]=p->GetQ()*36./23333.;  //electrons -> ADC
509     Int_t lab[4]; 
510     lab[0]=p->GetLabel(0); lab[1]=p->GetLabel(1); lab[2]=p->GetLabel(2);
511     lab[3]=fNdet[idx];
512     CheckLabels(lab);
513     Int_t dummy[3]={0,0,0};
514     new (cl[i]) AliITSclusterV2(lab,lp, dummy);
515   }  
516
517 */