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"
30 const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
31 const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
34 ClassImp(AliITSClusterFinderSSD)
36 //____________________________________________________________________
39 //____________________________________________________________________
43 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg, TClonesArray *digits)
45 //Standard constructor
50 fMap = new AliITSMapA1(fSegmentation,fDigits);
52 fITS=(AliITS*)gAlice->GetModule("ITS");
54 fClusterP = new TClonesArray ("AliITSclusterSSD",200);
57 fClusterN= new TClonesArray ("AliITSclusterSSD",200);
60 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
75 //-------------------------------------------------------
76 AliITSClusterFinderSSD::~AliITSClusterFinderSSD() {
88 //-------------------------------------------------------
89 void AliITSClusterFinderSSD::InitReconstruction()
91 // initialization of the cluster finder
93 register Int_t i; //iterator
95 for (i=0;i<fNClusterP;i++)
97 fClusterP->RemoveAt(i);
100 for (i=0;i<fNClusterN;i++)
102 fClusterN->RemoveAt(i);
106 for (i=0;i<fNPackages;i++)
108 fPackages->RemoveAt(i);
115 Float_t StereoP,StereoN;
116 fSegmentation->Angles(StereoP,StereoN);
118 CalcStepFactor (StereoP,StereoN);
120 if (debug) cout<<"fSFF = "<<fSFF<<" fSFB = "<<fSFB<<"\n";
124 //---------------------------------------------
125 void AliITSClusterFinderSSD::FindRawClusters(Int_t module)
127 // This function findes out all clusters belonging to one module
128 // 1. Zeroes all space after previous module reconstruction
129 // 2. Finds all neighbouring digits, create clusters
130 // 3. If necesery, resolves for each group of neighbouring digits
131 // how many clusters creates it.
132 // 4. Colculate the x and z coordinate
134 Int_t lay, lad, detect;
135 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
136 AliITSgeom *geom = aliITS->GetITSgeom();
137 geom->GetModuleId(module,lay, lad, detect);
139 ((AliITSsegmentationSSD*)fSegmentation)->SetLayer(6);
141 ((AliITSsegmentationSSD*)fSegmentation)->SetLayer(5);
143 InitReconstruction(); //ad. 1
147 FindNeighbouringDigits(); //ad. 2
148 //SeparateOverlappedClusters(); //ad. 3
149 ClustersToPackages(); //ad. 4
155 //-------------------------------------------------
156 void AliITSClusterFinderSSD::FindNeighbouringDigits()
158 //If there are any digits on this side, create 1st Cluster,
159 // add to it this digit, and increment number of clusters
165 if ((fNDigitsP>0 ) && (fNDigitsN > 0 )) {
167 Int_t currentstripNo;
168 Int_t *dbuffer = new Int_t [300]; //buffer for strip numbers
169 Int_t dnumber; //curent number of digits in buffer
170 TArrayI &lDigitsIndexP = *fDigitsIndexP;
171 TArrayI &lDigitsIndexN = *fDigitsIndexN;
172 TObjArray &lDigits = *(Digits());
173 TClonesArray &lClusterP = *fClusterP;
174 TClonesArray &lClusterN = *fClusterN;
178 dbuffer[0]=lDigitsIndexP[0];
179 //If next digit is a neighbour of previous, adds to last cluster this digit
180 for (i=1; i<fNDigitsP; i++) {
182 currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->
184 //check if it is a neighbour of a previous one
185 if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->GetStripNumber())
186 == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexP[i];
188 //create a new one side cluster
189 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEP);
190 dbuffer[0]=lDigitsIndexP[i];
193 } // end loop over fNDigitsP
194 new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEP);
198 //for comments, see above
200 dbuffer[0]=lDigitsIndexN[0];
201 //If next digit is a neighbour of previous, adds to last cluster this digit
202 for (i=1; i<fNDigitsN; i++) {
203 currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->
205 if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->GetStripNumber())
206 == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexN[i];
208 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEN);
209 dbuffer[0]=lDigitsIndexN[i];
212 } // end loop over fNDigitsN
213 new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEN);
216 } // end condition on NDigits
218 if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP<<" fNClusterN ="<<fNClusterN<<"\n";
222 //--------------------------------------------------------------
224 void AliITSClusterFinderSSD::SeparateOverlappedClusters()
226 // overlapped clusters separation
228 register Int_t i; //iterator
230 Float_t factor=0.75; // How many percent must be lower signal
231 // on the middle one digit
232 // from its neighbours
233 Int_t signal0; //signal on the strip before the current one
234 Int_t signal1; //signal on the current one signal
235 Int_t signal2; //signal on the strip after the current one
236 TArrayI *splitlist; // List of splits
237 Int_t numerofsplits=0; // number of splits
238 Int_t initPsize = fNClusterP; //initial size of the arrays
239 Int_t initNsize = fNClusterN; //we have to keep it because it will grow
240 // in this function and it doasn't make
241 // sense to pass through it again
243 splitlist = new TArrayI(300);
245 for (i=0;i<initPsize;i++)
247 if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==1) continue;
248 if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==2) continue;
249 Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
250 for (Int_t j=1; j<nj; j++)
252 signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
253 signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
254 signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
255 //if signal is less then factor*signal of its neighbours
256 if ( (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
258 (*splitlist)[numerofsplits++]=j;
260 } // end loop over number of digits
261 //split this cluster if necessary
262 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
265 //in signed places (splitlist)
266 } // end loop over clusters on Pside
268 for (i=0;i<initNsize;i++) {
269 if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==1) continue;
270 if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==2) continue;
271 Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
272 for (Int_t j=1; j<nj; j++)
274 signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
275 signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
276 signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
277 //if signal is less then factor*signal of its neighbours
278 if ( (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
279 (*splitlist)[numerofsplits++]=j;
280 } // end loop over number of digits
281 //split this cluster into more clusters
282 if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
283 numerofsplits=0; //in signed places (splitlist)
284 } // end loop over clusters on Nside
289 //-------------------------------------------------------
290 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits, Int_t index, Bool_t side)
292 //This function splits one side cluster into more clusters
293 //number of splits is defined by "nsplits"
294 //Place of splits are defined in the TArray "list"
296 // For further optimisation: Replace this function by two
297 // specialised ones (each for one side)
300 //For comlete comments see AliITSclusterSSD::SplitCluster
303 register Int_t i; //iterator
305 AliITSclusterSSD* curentcluster;
306 Int_t *tmpdigits = new Int_t[100];
308 // side true means P side
310 curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
311 for (i = nsplits; i>0 ;i--) {
312 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
313 new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(NN,tmpdigits,Digits(),side);
314 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
315 SetLeftNeighbour(kTRUE);
316 //if left cluster had neighbour on the right before split
317 //new should have it too
318 if ( curentcluster->GetRightNeighbour() )
319 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
320 SetRightNeighbour(kTRUE);
321 else curentcluster->SetRightNeighbour(kTRUE);
323 } // end loop over nplits
325 curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
326 for (i = nsplits; i>0 ;i--) {
327 NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
328 new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(NN,tmpdigits,Digits(),side);
329 ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
330 SetRightNeighbour(kTRUE);
331 if (curentcluster->GetRightNeighbour())
332 ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
333 SetRightNeighbour(kTRUE);
334 else curentcluster->SetRightNeighbour(kTRUE);
336 } // end loop over nplits
343 //-------------------------------------------------
344 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end)
346 // sort digits on the P side
351 if (start != (end - 1) )
353 left=this->SortDigitsP(start,(start+end)/2);
354 right=this->SortDigitsP((start+end)/2,end);
355 return (left || right);
359 left = ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[start]]))->GetStripNumber();
360 right= ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[end]]))->GetStripNumber();
363 Int_t tmp = (*fDigitsIndexP)[start];
364 (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
365 (*fDigitsIndexP)[end]=tmp;
373 //--------------------------------------------------
375 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end)
377 // sort digits on the N side
381 if (start != (end - 1))
383 left=this->SortDigitsN(start,(start+end)/2);
384 right=this->SortDigitsN((start+end)/2,end);
385 return (left || right);
389 left =((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[start]]))->GetStripNumber();
390 right=((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[end]]))->GetStripNumber();
393 Int_t tmp = (*fDigitsIndexN)[start];
394 (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
395 (*fDigitsIndexN)[end]=tmp;
402 //------------------------------------------------
403 void AliITSClusterFinderSSD::FillDigitsIndex()
405 //Fill the indexes of the clusters belonging to a given ITS module
412 N = fDigits->GetEntriesFast();
414 Int_t* PSidx = new Int_t [N*sizeof(Int_t)];
415 Int_t* NSidx = new Int_t [N*sizeof(Int_t)];
416 if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(N);
417 if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(N);
421 for ( i = 0 ; i< N; i++ ) {
422 dig = (AliITSdigitSSD*)GetDigit(i);
425 tmp=dig->GetStripNumber();
426 // I find this totally unnecessary - it's just a
427 // CPU consuming double check
432 if (debug) cout<<"Such a digit exists \n";
438 fDigitsIndexP->AddAt(i,fNDigitsP++);
443 tmp=dig->GetStripNumber();
449 if (debug) cout<<"Such a digit exists \n";
455 fDigitsIndexN->AddAt(i,fNDigitsN++);
464 if (debug) cout<<"Digits : P = "<<fNDigitsP<<" N = "<<fNDigitsN<<endl;
469 //-------------------------------------------
471 void AliITSClusterFinderSSD::SortDigits()
477 for (i=0;i<fNDigitsP-1;i++)
478 if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
481 for (i=0;i<fNDigitsN-1;i++)
482 if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
487 //----------------------------------------------
488 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN)
490 // fill cluster index array
493 for (i=0; i<fNClusterP;i++)
497 for (i=0; i<fNClusterN;i++)
504 //------------------------------------------------------
505 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN)
511 for (i=0;i<fNClusterP-1;i++)
512 if (SortClustersP(0,(fNClusterP-1),arrayP)==0) break;
516 for (i=0;i<fNClusterN-1;i++)
517 if (SortClustersN(0,(fNClusterN-1),arrayN)==0) break;
523 //---------------------------------------------------
524 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end, Int_t *array)
526 //Sort P side clusters
531 if (start != (end - 1) ) {
532 left=this->SortClustersP(start,(start+end)/2,array);
533 right=this->SortClustersP((start+end)/2,end,array);
534 return (left || right);
536 left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
538 right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
541 Int_t tmp = array[start];
542 array[start]=array[end];
553 //-------------------------------------------------------
554 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end, Int_t *array)
556 //Sort N side clusters
562 if (start != (end - 1) ) {
563 left=this->SortClustersN(start,(start+end)/2,array);
564 right=this->SortClustersN((start+end)/2,end,array);
565 return (left || right);
567 left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
569 right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
572 Int_t tmp = array[start];
573 array[start]=array[end];
583 //-------------------------------------------------------
584 void AliITSClusterFinderSSD::ClustersToPackages()
588 Int_t *oneSclP = new Int_t[fNClusterP]; //I want to have sorted 1 S clusters
589 Int_t *oneSclN = new Int_t[fNClusterN]; //I can not sort it in TClonesArray
590 //so, I create table of indexes and
592 //I do not use TArrayI on purpose
593 // MB: well, that's not true that one
594 //cannot sort objects in TClonesArray
596 //AliITSpackageSSD *currentpkg;
597 AliITSclusterSSD *currentP;
598 AliITSclusterSSD *currentN;
601 //Fills in One Side Clusters Index Array
602 FillClIndexArrays(oneSclP,oneSclN);
603 //Sorts filled Arrays
604 SortClusters(oneSclP,oneSclN);
607 new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
608 //currentpkg = (AliITSpackageSSD*)((*fPackages)[0]);
610 // Take all pairs of recpoints x coordinates in both sides
611 // to calculate z coordinates of the recpoints
612 for (j1=0;j1<fNClusterP;j1++) {
613 currentP = GetPSideCluster(oneSclP[j1]);
614 Double_t xP = currentP->GetPosition();
615 //Int_t NxP = currentP->GetNumOfDigits();
616 Float_t signalP = currentP->GetTotalSignal();
617 for (j2=0;j2<fNClusterN;j2++) {
618 currentN = GetNSideCluster(oneSclN[j2]);
619 Double_t xN = currentN->GetPosition();
620 //Int_t NxN = currentN->GetNumOfDigits();
621 Float_t signalN = currentN->GetTotalSignal();
622 CreateNewRecPoint(xP,1,xN,1,signalP,signalN,currentP, currentN, 0.75);
633 //------------------------------------------------------
634 Bool_t AliITSClusterFinderSSD::
635 CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
636 Float_t SigP,Float_t SigN,
637 AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN,
640 // create the recpoints
641 const Float_t kADCtoKeV = 2.16;
642 // 50 ADC units -> 30000 e-h pairs; 1e-h pair -> 3.6e-3 KeV;
643 // 1 ADC unit -> (30000/50)*3.6e-3 = 2.16 KeV
644 const Float_t kconv = 1.0e-4;
646 const Float_t kRMSx = 20.0*kconv;
647 const Float_t kRMSz = 800.0*kconv;
653 if (GetCrossing(P,N)) {
655 //GetCrossingError(dP,dN);
656 AliITSRawClusterSSD cnew;
657 Int_t nstripsP=clusterP->GetNumOfDigits();
658 Int_t nstripsN=clusterN->GetNumOfDigits();
664 dedx = SigP*kADCtoKeV;
667 dedx = SigN*kADCtoKeV;
670 Float_t signalCut = TMath::Abs((SigP-SigN)/signal);
673 tr = (Int_t*) clusterP->GetTracks(n);
674 NTracks = clusterP->GetNTracks();
675 //cout<<"NewRec: NTracks,tr0,tr1,tr2 ="<<NTracks<<","<<tr[0]<<","<<tr[1]<<","<<tr[2]<<endl;
676 if(NTracks>0 && signalCut<0.18) {
677 // the cut of the signals difference in P and N sides to
678 // remove the "ghosts"
681 cnew.fMultiplicity=nstripsP;
682 cnew.fMultiplicityN=nstripsN;
683 cnew.fQErr=signalCut;
684 cnew.fNtracks=NTracks;
686 fITS->AddCluster(2,&cnew);
693 rnew.SetSigmaX2( kRMSx* kRMSx);
694 rnew.SetSigmaZ2( kRMSz* kRMSz);
695 rnew.fTracks[0]=tr[0];
696 rnew.fTracks[1]=tr[1];
697 rnew.fTracks[2]=tr[2];
699 fITS->AddRecPoint(rnew);
708 //------------------------------------------------------
709 void AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo )
711 // calculate the step factor for matching clusters
714 // 95 is the pitch, 4000 - dimension along z ?
717 Float_t dz=fSegmentation->Dz();
719 fSFF = ( (Int_t) (Psteo*dz/fPitch ) );// +1;
720 fSFB = ( (Int_t) (Nsteo*dz/fPitch ) );// +1;
725 //-----------------------------------------------------------
726 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx)
728 // get P side clusters
733 if((idx<0)||(idx>=fNClusterP))
735 printf("AliITSClusterFinderSSD::GetPSideCluster : index out of range\n");
740 return (AliITSclusterSSD*)((*fClusterP)[idx]);
744 //-------------------------------------------------------
745 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx)
747 // get N side clusters
749 if((idx<0)||(idx>=fNClusterN))
751 printf("AliITSClusterFinderSSD::GetNSideCluster : index out of range\n");
756 return (AliITSclusterSSD*)((*fClusterN)[idx]);
760 //--------------------------------------------------------
761 AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side)
765 return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
768 //_______________________________________________________________________
770 Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
774 Float_t Dx = fSegmentation->Dx(); // detector size in x direction, microns
775 Float_t Dz = fSegmentation->Dz(); // detector size in z direction, microns
777 Float_t xL; // x local coordinate
778 Float_t zL; // z local coordinate
779 Float_t x; // x = xL + Dx/2
780 Float_t z; // z = zL + Dz/2
781 Float_t xP; // x coordinate in the P side from the first P strip
782 Float_t xN; // x coordinate in the N side from the first N strip
783 Float_t StereoP,StereoN;
785 fSegmentation->Angles(StereoP,StereoN);
786 fTanP=TMath::Tan(StereoP);
787 fTanN=TMath::Tan(StereoN);
788 Float_t kP = fTanP; // Tangent of 0.0075 mrad
789 Float_t kN = fTanN; // Tangent of 0.0275 mrad
793 xP = N; // change the mistake for the P/N
794 xN = P; // coordinates correspondence in this function
795 x = xP + kP*(Dz*kN-xP+xN)/(kP+kN);
796 z = (Dz*kN-xP+xN)/(kP+kN);
802 if(TMath::Abs(xL) > Dx/2 || TMath::Abs(zL) > Dz/2) return kFALSE;
803 // Check that xL and zL are inside the detector for the
804 // correspondent xP and xN coordinates
810 //_________________________________________________________________________
812 void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN)
814 // get crossing error
818 dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
819 dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
820 TMath::Abs(dN *fTanP/(fTanP + fTanN) ));