1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 #include "AliITSdigit.h"
21 #include "AliITSRawCluster.h"
22 #include "AliITSRecPoint.h"
23 #include "AliITSMapA1.h"
24 #include "AliITSClusterFinderSSD.h"
25 #include "AliITSclusterSSD.h"
26 #include "AliITSpackageSSD.h"
27 #include "AliITSsegmentation.h"
28 #include "AliITSgeom.h"
29 #include "AliITSsegmentationSSD.h"
31 const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
32 const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
35 ClassImp(AliITSClusterFinderSSD)
37 //______________________________________________________________________
40 //______________________________________________________________________
44 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg,
45 TClonesArray *digits){
46 //Standard constructor
51 fMap = new AliITSMapA1(fSegmentation,fDigits);
53 fITS=(AliITS*)gAlice->GetModule("ITS");
55 fClusterP = new TClonesArray ("AliITSclusterSSD",200);
58 fClusterN= new TClonesArray ("AliITSclusterSSD",200);
61 fPackages = new TClonesArray ("AliITSpackageSSD",200); //packages
64 fDigitsIndexP = new TArrayI(300);
67 fDigitsIndexN = new TArrayI(300);
70 fPitch = fSegmentation->Dpx(0);
71 fPNsignalRatio=7./8.; // warning: hard-wired number
74 //----------------------------------------------------------------------
75 AliITSClusterFinderSSD::~AliITSClusterFinderSSD() {
85 //----------------------------------------------------------------------
86 void AliITSClusterFinderSSD::InitReconstruction(){
87 // initialization of the cluster finder
88 register Int_t i; //iterator
90 for (i=0;i<fNClusterP;i++){
91 fClusterP->RemoveAt(i);
94 for (i=0;i<fNClusterN;i++){
95 fClusterN->RemoveAt(i);
99 for (i=0;i<fNPackages;i++){
100 fPackages->RemoveAt(i);
107 Float_t StereoP,StereoN;
108 fSegmentation->Angles(StereoP,StereoN);
110 CalcStepFactor (StereoP,StereoN);
112 if (debug) cout<<"fSFF = "<<fSFF<<" fSFB = "<<fSFB<<"\n";
114 //----------------------------------------------------------------------
115 void AliITSClusterFinderSSD::FindRawClusters(Int_t module){
116 // This function findes out all clusters belonging to one module
117 // 1. Zeroes all space after previous module reconstruction
118 // 2. Finds all neighbouring digits, create clusters
119 // 3. If necesery, resolves for each group of neighbouring digits
120 // how many clusters creates it.
121 // 4. Colculate the x and z coordinate
122 Int_t lay, lad, detect;
123 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
124 AliITSgeom *geom = aliITS->GetITSgeom();
125 geom->GetModuleId(module,lay, lad, detect);
126 if ( lay == 6 )((AliITSsegmentationSSD*)fSegmentation)->SetLayer(6);
127 if ( lay == 5 )((AliITSsegmentationSSD*)fSegmentation)->SetLayer(5);
129 InitReconstruction(); //ad. 1
133 FindNeighbouringDigits(); //ad. 2
134 //SeparateOverlappedClusters(); //ad. 3
135 ClustersToPackages(); //ad. 4
139 //----------------------------------------------------------------------
140 void AliITSClusterFinderSSD::FindNeighbouringDigits(){
141 //If there are any digits on this side, create 1st Cluster,
142 // add to it this digit, and increment number of clusters
145 if ((fNDigitsP>0 ) && (fNDigitsN > 0 )) {
147 Int_t currentstripNo;
148 Int_t *dbuffer = new Int_t [300]; //buffer for strip numbers
149 Int_t dnumber; //curent number of digits in buffer
150 TArrayI &lDigitsIndexP = *fDigitsIndexP;
151 TArrayI &lDigitsIndexN = *fDigitsIndexN;
152 TObjArray &lDigits = *(Digits());
153 TClonesArray &lClusterP = *fClusterP;
154 TClonesArray &lClusterN = *fClusterN;
158 dbuffer[0]=lDigitsIndexP[0];
159 // If next digit is a neighbour of previous, adds to last cluster
161 for (i=1; i<fNDigitsP; i++) {
163 currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->
165 //check if it is a neighbour of a previous one
166 if((((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->
167 GetStripNumber()) == (currentstripNo-1) ){
168 dbuffer[dnumber++]=lDigitsIndexP[i];
170 //create a new one side cluster
171 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,
172 dbuffer,Digits(),fgkSIDEP);
173 dbuffer[0]=lDigitsIndexP[i];
176 } // end loop over fNDigitsP
177 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
180 //for comments, see above
182 dbuffer[0]=lDigitsIndexN[0];
183 // If next digit is a neighbour of previous, adds to last cluster
185 for (i=1; i<fNDigitsN; i++) {
186 currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->
188 if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->
189 GetStripNumber()) == (currentstripNo-1) ){
190 dbuffer[dnumber++]=lDigitsIndexN[i];
192 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
194 dbuffer[0]=lDigitsIndexN[i];
197 } // end loop over fNDigitsN
198 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
201 } // end condition on NDigits
203 if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP<<
204 " fNClusterN ="<<fNClusterN<<"\n";
206 //----------------------------------------------------------------------
207 void AliITSClusterFinderSSD::SeparateOverlappedClusters(){
208 // overlapped clusters separation
209 register Int_t i; //iterator
210 Float_t factor=0.75; // How many percent must be lower signal
211 // on the middle one digit
212 // from its neighbours
213 Int_t signal0; //signal on strip before the current one
214 Int_t signal1; //signal on the current one signal
215 Int_t signal2; //signal on the strip after the current one
216 TArrayI *splitlist; // List of splits
217 Int_t numerofsplits=0; // number of splits
218 Int_t initPsize = fNClusterP;//initial size of the arrays
219 Int_t initNsize = fNClusterN;//we have to keep it because it will grow
220 // in this function and it doasn't make
221 // sense to pass through it again
223 splitlist = new TArrayI(300);
225 for (i=0;i<initPsize;i++){
226 if((((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==1) continue;
227 if((((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==2) continue;
228 Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
229 for (Int_t j=1; j<nj; j++){
230 signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
231 signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
232 signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
233 //if signal is less then factor*signal of its neighbours
234 if ( (signal1<(factor*signal0)) && (signal1<(factor*signal2)) ){
235 (*splitlist)[numerofsplits++]=j;
237 } // end loop over number of digits
238 //split this cluster if necessary
239 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
242 //in signed places (splitlist)
243 } // end loop over clusters on Pside
245 for (i=0;i<initNsize;i++) {
246 if((((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==1) continue;
247 if((((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==2) continue;
248 Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
249 for (Int_t j=1; j<nj; j++){
250 signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
251 signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
252 signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
253 //if signal is less then factor*signal of its neighbours
254 if ( (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
255 (*splitlist)[numerofsplits++]=j;
256 } // end loop over number of digits
257 //split this cluster into more clusters
258 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
259 numerofsplits=0; //in signed places (splitlist)
260 } // end loop over clusters on Nside
264 //----------------------------------------------------------------------
265 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits,
266 Int_t index, Bool_t side){
267 //This function splits one side cluster into more clusters
268 //number of splits is defined by "nsplits"
269 //Place of splits are defined in the TArray "list"
270 // For further optimisation: Replace this function by two
271 // specialised ones (each for one side)
273 //For comlete comments see AliITSclusterSSD::SplitCluster
274 register Int_t i; //iterator
275 AliITSclusterSSD* curentcluster;
276 Int_t *tmpdigits = new Int_t[100];
278 // side true means P side
281 curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
282 for (i = nsplits; i>0 ;i--) {
283 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
284 new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(NN,
285 tmpdigits,Digits(),side);
286 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
287 SetLeftNeighbour(kTRUE);
288 //if left cluster had neighbour on the right before split
289 //new should have it too
290 if ( curentcluster->GetRightNeighbour() )
291 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
292 SetRightNeighbour(kTRUE);
293 else curentcluster->SetRightNeighbour(kTRUE);
295 } // end loop over nplits
297 curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
298 for (i = nsplits; i>0 ;i--) {
299 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
300 new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(NN,
301 tmpdigits,Digits(),side);
302 ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
303 SetRightNeighbour(kTRUE);
304 if (curentcluster->GetRightNeighbour())
305 ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
306 SetRightNeighbour(kTRUE);
307 else curentcluster->SetRightNeighbour(kTRUE);
309 } // end loop over nplits
313 //----------------------------------------------------------------------
314 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end){
315 // sort digits on the P side
318 if (start != (end - 1) ){
319 left=this->SortDigitsP(start,(start+end)/2);
320 right=this->SortDigitsP((start+end)/2,end);
321 return (left || right);
323 left = ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[start]]))->
325 right= ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[end]]))->
328 Int_t tmp = (*fDigitsIndexP)[start];
329 (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
330 (*fDigitsIndexP)[end]=tmp;
336 //----------------------------------------------------------------------
337 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end){
338 // sort digits on the N side
342 if (start != (end - 1)){
343 left=this->SortDigitsN(start,(start+end)/2);
344 right=this->SortDigitsN((start+end)/2,end);
345 return (left || right);
347 left =((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[start]]))->
349 right=((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[end]]))->
352 Int_t tmp = (*fDigitsIndexN)[start];
353 (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
354 (*fDigitsIndexN)[end]=tmp;
359 //----------------------------------------------------------------------
360 void AliITSClusterFinderSSD::FillDigitsIndex(){
361 //Fill the indexes of the clusters belonging to a given ITS module
367 N = fDigits->GetEntriesFast();
369 Int_t* PSidx = new Int_t [N*sizeof(Int_t)];
370 Int_t* NSidx = new Int_t [N*sizeof(Int_t)];
371 if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(N);
372 if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(N);
376 for ( i = 0 ; i< N; i++ ) {
377 dig = (AliITSdigitSSD*)GetDigit(i);
380 tmp=dig->GetStripNumber();
381 // I find this totally unnecessary - it's just a
382 // CPU consuming double check
385 if (debug) cout<<"Such a digit exists \n";
391 fDigitsIndexP->AddAt(i,fNDigitsP++);
396 tmp=dig->GetStripNumber();
400 if (debug) cout<<"Such a digit exists \n";
406 fDigitsIndexN->AddAt(i,fNDigitsN++);
409 } // end for dig->IsSideP()
415 if (debug) cout<<"Digits : P = "<<fNDigitsP<<" N = "<<fNDigitsN<<endl;
418 //----------------------------------------------------------------------
419 void AliITSClusterFinderSSD::SortDigits(){
424 for (i=0;i<fNDigitsP-1;i++)
425 if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
428 for (i=0;i<fNDigitsN-1;i++)
429 if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
431 //----------------------------------------------------------------------
432 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN){
433 // fill cluster index array
436 for (i=0; i<fNClusterP;i++){
439 for (i=0; i<fNClusterN;i++){
443 //----------------------------------------------------------------------
444 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN){
449 for (i=0;i<fNClusterP-1;i++)
450 if (SortClustersP(0,(fNClusterP-1),arrayP)==0) break;
452 for (i=0;i<fNClusterN-1;i++)
453 if (SortClustersN(0,(fNClusterN-1),arrayN)==0) break;
455 //----------------------------------------------------------------------
456 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end,
458 //Sort P side clusters
462 if (start != (end - 1) ) {
463 left=this->SortClustersP(start,(start+end)/2,array);
464 right=this->SortClustersP((start+end)/2,end,array);
465 return (left || right);
467 left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
469 right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
472 Int_t tmp = array[start];
473 array[start]=array[end];
479 //-------------------------------------------------------
480 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end,
482 //Sort N side clusters
486 if (start != (end - 1) ) {
487 left=this->SortClustersN(start,(start+end)/2,array);
488 right=this->SortClustersN((start+end)/2,end,array);
489 return (left || right);
491 left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
493 right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
496 Int_t tmp = array[start];
497 array[start]=array[end];
503 //----------------------------------------------------------------------
504 void AliITSClusterFinderSSD::ClustersToPackages(){
506 Int_t *oneSclP = new Int_t[fNClusterP]; //I want to have sorted 1 S clusters
507 Int_t *oneSclN = new Int_t[fNClusterN]; //I can not sort it in TClonesArray
508 //so, I create table of indexes and
510 //I do not use TArrayI on purpose
511 // MB: well, that's not true that one
512 //cannot sort objects in TClonesArray
514 //AliITSpackageSSD *currentpkg;
515 AliITSclusterSSD *currentP;
516 AliITSclusterSSD *currentN;
519 //Fills in One Side Clusters Index Array
520 FillClIndexArrays(oneSclP,oneSclN);
521 //Sorts filled Arrays
522 SortClusters(oneSclP,oneSclN);
525 new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
526 //currentpkg = (AliITSpackageSSD*)((*fPackages)[0]);
528 // Take all pairs of recpoints x coordinates in both sides
529 // to calculate z coordinates of the recpoints
530 for (j1=0;j1<fNClusterP;j1++) {
531 currentP = GetPSideCluster(oneSclP[j1]);
532 Double_t xP = currentP->GetPosition();
533 //Int_t NxP = currentP->GetNumOfDigits();
534 Float_t signalP = currentP->GetTotalSignal();
535 for (j2=0;j2<fNClusterN;j2++) {
536 currentN = GetNSideCluster(oneSclN[j2]);
537 Double_t xN = currentN->GetPosition();
538 //Int_t NxN = currentN->GetNumOfDigits();
539 Float_t signalN = currentN->GetTotalSignal();
540 CreateNewRecPoint(xP,1,xN,1,signalP,signalN,currentP, currentN, 0.75);
551 //------------------------------------------------------
552 Bool_t AliITSClusterFinderSSD::
553 CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
554 Float_t SigP,Float_t SigN,
555 AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN,
558 // create the recpoints
559 const Float_t kADCtoKeV = 2.16;
560 // 50 ADC units -> 30000 e-h pairs; 1e-h pair -> 3.6e-3 KeV;
561 // 1 ADC unit -> (30000/50)*3.6e-3 = 2.16 KeV
562 const Float_t kconv = 1.0e-4;
564 const Float_t kRMSx = 20.0*kconv;
565 const Float_t kRMSz = 800.0*kconv;
569 if (GetCrossing(P,N)) {
571 //GetCrossingError(dP,dN);
572 AliITSRawClusterSSD cnew;
573 Int_t nstripsP=clusterP->GetNumOfDigits();
574 Int_t nstripsN=clusterN->GetNumOfDigits();
580 dedx = SigP*kADCtoKeV;
583 dedx = SigN*kADCtoKeV;
586 Float_t signalCut = TMath::Abs((SigP-SigN)/signal);
590 cnew.fMultiplicity=nstripsP;
591 cnew.fMultiplicityN=nstripsN;
592 cnew.fQErr=signalCut;
594 //fITS->AddCluster(2,&cnew);
596 if(signalCut<0.18) fITS->AddCluster(2,&cnew);
597 // the cut of the signals difference in P and N sides to
598 // remove the "ghosts"
600 tr = (Int_t*) clusterP->GetTracks(n);
607 rnew.SetSigmaX2( kRMSx* kRMSx);
608 rnew.SetSigmaZ2( kRMSz* kRMSz);
609 rnew.fTracks[0]=tr[0];
610 rnew.fTracks[1]=tr[1];
611 rnew.fTracks[2]=tr[2];
613 //fITS->AddRecPoint(rnew);
615 if(signalCut<0.18) fITS->AddRecPoint(rnew);
616 // the cut of the signals difference in P and N sides to
617 // remove the "ghosts"
625 //------------------------------------------------------
626 void AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo )
628 // calculate the step factor for matching clusters
631 // 95 is the pitch, 4000 - dimension along z ?
634 Float_t dz=fSegmentation->Dz();
636 fSFF = ( (Int_t) (Psteo*dz/fPitch ) );// +1;
637 fSFB = ( (Int_t) (Nsteo*dz/fPitch ) );// +1;
642 //-----------------------------------------------------------
643 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx)
645 // get P side clusters
650 if((idx<0)||(idx>=fNClusterP))
652 printf("AliITSClusterFinderSSD::GetPSideCluster : index out of range\n");
657 return (AliITSclusterSSD*)((*fClusterP)[idx]);
661 //-------------------------------------------------------
662 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx)
664 // get N side clusters
666 if((idx<0)||(idx>=fNClusterN))
668 printf("AliITSClusterFinderSSD::GetNSideCluster : index out of range\n");
673 return (AliITSclusterSSD*)((*fClusterN)[idx]);
677 //--------------------------------------------------------
678 AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side)
682 return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
685 //_______________________________________________________________________
687 Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
691 Float_t Dx = fSegmentation->Dx(); // detector size in x direction, microns
692 Float_t Dz = fSegmentation->Dz(); // detector size in z direction, microns
694 Float_t xL; // x local coordinate
695 Float_t zL; // z local coordinate
696 Float_t x; // x = xL + Dx/2
697 Float_t z; // z = zL + Dz/2
698 Float_t xP; // x coordinate in the P side from the first P strip
699 Float_t xN; // x coordinate in the N side from the first N strip
700 Float_t StereoP,StereoN;
702 fSegmentation->Angles(StereoP,StereoN);
703 fTanP=TMath::Tan(StereoP);
704 fTanN=TMath::Tan(StereoN);
705 Float_t kP = fTanP; // Tangent of 0.0075 mrad
706 Float_t kN = fTanN; // Tangent of 0.0275 mrad
710 xP = N; // change the mistake for the P/N
711 xN = P; // coordinates correspondence in this function
712 x = xP + kP*(Dz*kN-xP+xN)/(kP+kN);
713 z = (Dz*kN-xP+xN)/(kP+kN);
719 if(TMath::Abs(xL) > Dx/2 || TMath::Abs(zL) > Dz/2) return kFALSE;
720 // Check that xL and zL are inside the detector for the
721 // correspondent xP and xN coordinates
727 //_________________________________________________________________________
729 void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN)
731 // get crossing error
735 dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
736 dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
737 TMath::Abs(dN *fTanP/(fTanP + fTanN) ));