ITS/AliITSClusterFinderSSD.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 //  * The package was revised and changed by Boris Batiounia in the time     *
  18 //  * period of March - June 2001                                            *
  19 // **************************************************************************/
  20 //
  21 //#include <Riostream.h>
  22 #include <TArrayI.h>
  23
  24 //#include "AliRun.h"
  25 //#include "AliITS.h"
  26 #include "AliITSClusterFinderSSD.h"
  27 #include "AliITSDetTypeRec.h"
  28 #include "AliITSMapA1.h"
  29 #include "AliITSRawClusterSSD.h"
  30 #include "AliITSRecPoint.h"
  31 #include "AliITSdigitSSD.h"
  32 #include "AliITSclusterSSD.h"
  33 #include "AliITSpackageSSD.h"
  34 #include "AliITSsegmentationSSD.h"
  35 //#include "AliITSgeom.h"
  36 #include "AliITSCalibrationSSD.h"
  37 #include "AliLog.h"
  38
  39 const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
  40 const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
  41 const Int_t AliITSClusterFinderSSD::fgkNoiseThreshold=5;
  42 const Int_t debug=0;
  43
  44 ClassImp(AliITSClusterFinderSSD)
  45
  46   //____________________________________________________________________
  47   //
  48   //  Constructor
  49   //______________________________________________________________________
  50   AliITSClusterFinderSSD::AliITSClusterFinderSSD():
  51 AliITSClusterFinder(),
  52 fClusterP(0),
  53 fNClusterP(0),
  54 fClusterN(0),
  55 fNClusterN(0),
  56 fPackages(0),
  57 fNPackages(0),
  58 fDigitsIndexP(0),
  59 fNDigitsP(0),
  60 fDigitsIndexN(0),
  61 fNDigitsN(0),
  62 fPitch(0.0),
  63 fTanP(0.0),
  64 fTanN(0.0),
  65 fPNsignalRatio(0.0),
  66 fSFF(0),
  67 fSFB(0){
  68     //Default constructor
  69 }
  70
  71 //______________________________________________________________________
  72 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSDetTypeRec* dettyp,
  73                                                TClonesArray *digits):
  74 AliITSClusterFinder(dettyp,digits),
  75 fClusterP(0),
  76 fNClusterP(0),
  77 fClusterN(0),
  78 fNClusterN(0),
  79 fPackages(0),
  80 fNPackages(0),
  81 fDigitsIndexP(0),
  82 fNDigitsP(0),
  83 fDigitsIndexN(0),
  84 fNDigitsN(0),
  85 fPitch(0.0),
  86 fTanP(0.0),
  87 fTanN(0.0),
  88 fPNsignalRatio(0.0),
  89 fSFF(0),
  90 fSFB(0){
  91     //Standard constructor
  92
  93     SetDigits(digits);
  94     SetMap(new AliITSMapA1(GetSeg(),Digits()));
  95     fClusterP     = new TClonesArray ("AliITSclusterSSD",200);
  96     fNClusterP    = 0;
  97     fClusterN     = new TClonesArray ("AliITSclusterSSD",200);
  98     fNClusterN    = 0;
  99     fPackages     = new TClonesArray ("AliITSpackageSSD",200);    //packages
 100     fNPackages    = 0;
 101     fDigitsIndexP = new TArrayI(800);
 102     fNDigitsP     = 0;
 103     fDigitsIndexN = new TArrayI(800);
 104     fNDigitsN     = 0;
 105     fPitch        = GetSeg()->Dpx(0);
 106     fPNsignalRatio= 7./8.;    // warning: hard-wired number
 107 }
 108
 109 //______________________________________________________________________
 110 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSDetTypeRec* dettyp):
 111 AliITSClusterFinder(dettyp),
 112 fClusterP(0),
 113 fNClusterP(0),
 114 fClusterN(0),
 115 fNClusterN(0),
 116 fPackages(0),
 117 fNPackages(0),
 118 fDigitsIndexP(0),
 119 fNDigitsP(0),
 120 fDigitsIndexN(0),
 121 fNDigitsN(0),
 122 fPitch(0.0),
 123 fTanP(0.0),
 124 fTanN(0.0),
 125 fPNsignalRatio(0.0),
 126 fSFF(0),
 127 fSFB(0){
 128     //Standard constructor
 129
 130     fClusterP     = new TClonesArray ("AliITSclusterSSD",200);
 131     fNClusterP    = 0;
 132     fClusterN     = new TClonesArray ("AliITSclusterSSD",200);
 133     fNClusterN    = 0;
 134     fPackages     = new TClonesArray ("AliITSpackageSSD",200);    //packages
 135     fNPackages    = 0;
 136     fDigitsIndexP = new TArrayI(800);
 137     fNDigitsP     = 0;
 138     fDigitsIndexN = new TArrayI(800);
 139     fNDigitsN     = 0;
 140     fPitch        = GetSeg()->Dpx(0);
 141     fPNsignalRatio= 7./8.;    // warning: hard-wired number
 142 }
 143
 144 //______________________________________________________________________
 145 AliITSClusterFinderSSD::~AliITSClusterFinderSSD(){
 146   // Default destructor
 147
 148   delete fClusterP;
 149   delete fClusterN;
 150   delete fPackages;
 151   delete fDigitsIndexP;
 152   delete fDigitsIndexN;
 153   delete fMap;
 154 }
 155 //______________________________________________________________________
 156 void AliITSClusterFinderSSD::InitReconstruction(){
 157   // initialization of the cluster finder
 158
 159   register Int_t i; //iterator
 160
 161     for (i=0;i<fNClusterP;i++) fClusterP->RemoveAt(i);
 162     fNClusterP  =0;
 163     for (i=0;i<fNClusterN;i++) fClusterN->RemoveAt(i);
 164     fNClusterN=0;
 165     for (i=0;i<fNPackages;i++) fPackages->RemoveAt(i);
 166     fNPackages = 0;
 167     fNDigitsP  = 0;
 168     fNDigitsN  = 0;
 169     Float_t stereoP,stereoN;
 170     //fSegmentation->Angles(stereoP,stereoN);
 171      GetSeg()->Angles(stereoP,stereoN);
 172    CalcStepFactor(stereoP,stereoN);
 173     if (debug) cout<<"fSFF = "<<fSFF<<"  fSFB = "<<fSFB<<"\n";
 174 }
 175 //______________________________________________________________________
 176 void AliITSClusterFinderSSD::FindRawClusters(Int_t module){
 177   // This function findes out all clusters belonging to one module
 178   // 1. Zeroes all space after previous module reconstruction
 179   // 2. Finds all neighbouring digits, create clusters
 180   // 3. If necesery, resolves for each group of neighbouring digits
 181   //    how many clusters creates it.
 182   // 4. Colculate the x and z coordinate
 183   Int_t lay, lad, detect;
 184
 185   //cout<<"clusterfinder: this is module "<<module<<endl;
 186
 187   //    AliITSgeom *geom = fITS->GetITSgeom();
 188
 189   if(!fDetTypeRec->GetITSgeom()) {
 190     Error("FindRawClusters","ITS geom is null!");
 191     return;
 192   }
 193
 194   SetModule(module);
 195   fDetTypeRec->GetITSgeom()->GetModuleId(GetModule(),lay, lad, detect);
 196   //   geom->GetModuleId(module,lay, lad, detect);
 197
 198   //cout<<"layer = "<<lay<<endl;
 199
 200   if ( lay == 6 ) ((AliITSsegmentationSSD*)GetSeg())->SetLayer(6);
 201   if ( lay == 5 ) ((AliITSsegmentationSSD*)GetSeg())->SetLayer(5);
 202
 203
 204   InitReconstruction();  //ad. 1
 205   fMap->FillMap();
 206
 207   FillDigitsIndex();
 208   if ( (fNDigitsP==0 )  || (fNDigitsN == 0 ))  return;
 209
 210   SortDigits();
 211   FindNeighbouringDigits(module); //ad. 2
 212
 213   //SeparateOverlappedClusters();  //ad. 3
 214   ClustersToPackages();  //ad. 4
 215
 216   fMap->ClearMap();
 217
 218 }
 219 //______________________________________________________________________
 220 void AliITSClusterFinderSSD::FindNeighbouringDigits(Int_t module){
 221     //If there are any digits on this side, create 1st Cluster,
 222     // add to it this digit, and increment number of clusters
 223
 224     register Int_t i,j;
 225     Int_t flag=0;
 226
 227     //if ( (fNDigitsP==0 )  || (fNDigitsN == 0 ))  return;
 228
 229     Int_t currentstripNo;
 230     Int_t *dbuffer = new Int_t [800];   //buffer for strip numbers
 231     Int_t dnumber;    //curent number of digits in buffer
 232
 233     TArrayI      &lDigitsIndexP = *fDigitsIndexP;
 234     TArrayI      &lDigitsIndexN = *fDigitsIndexN;
 235
 236     TObjArray    &lDigits       = *(Digits());
 237
 238     TClonesArray &lClusterP     = *fClusterP;
 239     TClonesArray &lClusterN     = *fClusterN;
 240
 241     //process P side
 242     dnumber = 1;
 243     dbuffer[0]=lDigitsIndexP[0];
 244
 245     //If next digit is a neigh. of previous, adds to last clust. this digit
 246     for (i=1; i<fNDigitsP; i++) {
 247
 248       //reads new digit
 249       currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->GetStripNumber();
 250
 251       if((((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->GetStripNumber()) ==  (currentstripNo-1) )
 252         dbuffer[dnumber++]=lDigitsIndexP[i];
 253
 254       else{
 255
 256         // check signal
 257         for(j=0;j<dnumber;j++) {
 258
 259           if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 260               > fgkNoiseThreshold* ((AliITSCalibrationSSD*)GetResp(module))->
 261               GetNoiseP( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 262             flag+=1;
 263
 264         }
 265
 266         //if(flag==dnumber) {
 267         if(flag>0) {
 268           //create a new one side cluster
 269           new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
 270                                                         Digits(),
 271                                                         fgkSIDEP);
 272
 273           flag=0;
 274         }
 275
 276         flag=0;
 277         dbuffer[0]=lDigitsIndexP[i];
 278         dnumber = 1;
 279
 280       } // end if else
 281
 282     } // end loop over fNDigitsP
 283
 284     // check signal
 285     for(j=0;j<dnumber;j++) {
 286
 287       if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 288           > fgkNoiseThreshold*((AliITSCalibrationSSD*)GetResp(module))->
 289           GetNoiseP( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 290         flag+=1;
 291
 292     }
 293
 294     //if(flag==dnumber) {
 295     if(flag>0) {
 296       //create a new one side cluster
 297       new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
 298                                                     Digits(),
 299                                                     fgkSIDEP);
 300
 301       flag=0;
 302     }
 303
 304     flag=0;
 305
 306     //process N side
 307     //for comments, see above
 308     dnumber = 1;
 309     dbuffer[0]=lDigitsIndexN[0];
 310     //If next digit is a neigh. of previous, adds to last clust. this digit
 311
 312     for (i=1; i<fNDigitsN; i++) {
 313       currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->GetStripNumber();
 314
 315       if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->GetStripNumber()) == (currentstripNo-1) )
 316         dbuffer[dnumber++]=lDigitsIndexN[i];
 317
 318       else {
 319
 320         // check signal
 321         for(j=0;j<dnumber;j++) {
 322
 323           if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 324               > fgkNoiseThreshold*((AliITSCalibrationSSD*)GetResp(module))->
 325               GetNoiseN( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 326             flag+=1;
 327
 328         }
 329
 330         //if(flag==dnumber) {
 331         if(flag>0) {
 332           //create a new one side cluster
 333           new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
 334                                                         Digits(),
 335                                                         fgkSIDEN);
 336
 337           flag=0;
 338         }
 339
 340         flag=0;
 341         dbuffer[0]=lDigitsIndexN[i];
 342         dnumber = 1;
 343       } // end if else
 344     } // end loop over fNDigitsN
 345
 346     // check signal
 347     for(j=0;j<dnumber;j++) {
 348
 349       if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 350           > fgkNoiseThreshold*((AliITSCalibrationSSD*)GetResp(module))->
 351           GetNoiseN( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 352         flag+=1;
 353
 354     }
 355
 356     //if(flag==dnumber) {
 357     if(flag>0) {
 358       //create a new one side cluster
 359       new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
 360                                                     Digits(),
 361                                                     fgkSIDEN);
 362
 363       flag=0;
 364     }
 365
 366     flag=0;
 367     delete [] dbuffer;
 368
 369     if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP
 370                    <<"  fNClusterN ="<<fNClusterN<<"\n";
 371 }
 372 //______________________________________________________________________
 373 void AliITSClusterFinderSSD::SeparateOverlappedClusters(){
 374     // overlapped clusters separation
 375     register Int_t i; //iterator
 376     Float_t  factor=0.75;            // How many percent must be lower signal
 377                                      // on the middle one digit
 378                                      // from its neighbours
 379     Int_t    signal0;              //signal on the strip before the current one
 380     Int_t    signal1;              //signal on the current one signal
 381     Int_t    signal2;              //signal on the strip after the current one
 382     TArrayI *splitlist;              //  List of splits
 383     Int_t    numerofsplits=0;        // number of splits
 384     Int_t    initPsize = fNClusterP; //initial size of the arrays
 385     Int_t    initNsize = fNClusterN; //we have to keep it because it will grow
 386                                      // in this function and it doasn't make
 387                                      // sense to pass through it again
 388     splitlist = new TArrayI(800);
 389
 390     for (i=0;i<initPsize;i++){
 391         if (( ((AliITSclusterSSD*)(*fClusterP)[i])->
 392               GetNumOfDigits())==1) continue;
 393         if (( ((AliITSclusterSSD*)(*fClusterP)[i])->
 394               GetNumOfDigits())==2) continue;
 395         Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
 396         for (Int_t j=1; j<nj; j++){
 397             signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
 398             signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
 399             signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
 400             //if signal is less then factor*signal of its neighbours
 401             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) ){
 402                 (*splitlist)[numerofsplits++]=j;
 403             } // end if
 404         } // end loop over number of digits
 405         //split this cluster if necessary
 406         if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
 407         numerofsplits=0;
 408         //in signed places (splitlist)
 409     } // end loop over clusters on Pside
 410
 411     for (i=0;i<initNsize;i++) {
 412         if (( ((AliITSclusterSSD*)(*fClusterN)[i])->
 413               GetNumOfDigits())==1) continue;
 414         if (( ((AliITSclusterSSD*)(*fClusterN)[i])->
 415               GetNumOfDigits())==2) continue;
 416         Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
 417         for (Int_t j=1; j<nj; j++){
 418             signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
 419             signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
 420             signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
 421             //if signal is less then factor*signal of its neighbours
 422             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
 423                 (*splitlist)[numerofsplits++]=j;
 424         } // end loop over number of digits
 425         //split this cluster into more clusters
 426         if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
 427         numerofsplits=0;
 428         //in signed places (splitlist)
 429     } // end loop over clusters on Nside
 430
 431     delete splitlist;
 432 }
 433 //______________________________________________________________________
 434 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits,
 435                                           Int_t index, Bool_t side){
 436     //This function splits one side cluster into more clusters
 437     //number of splits is defined by "nsplits"
 438     //Place of splits are defined in the TArray "list"
 439     // For further optimisation: Replace this function by two
 440     // specialised ones (each for one side)
 441     // save one "if"
 442     //For comlete comments see AliITSclusterSSD::SplitCluster
 443     register Int_t i; //iterator
 444     AliITSclusterSSD* curentcluster;
 445     Int_t   *tmpdigits = new Int_t[100];
 446     Int_t    nn;
 447
 448     // side true means P side
 449     if (side) {
 450         curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
 451         for (i = nsplits; i>0 ;i--) {
 452             nn=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
 453             new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(nn,tmpdigits,
 454                                                             Digits(),side);
 455             ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
 456                                                       SetLeftNeighbour(kTRUE);
 457             //if left cluster had neighbour on the right before split
 458             //new should have it too
 459             if ( curentcluster->GetRightNeighbour() )
 460                 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
 461                                                      SetRightNeighbour(kTRUE);
 462             else curentcluster->SetRightNeighbour(kTRUE);
 463             fNClusterP++;
 464         } // end loop over nplits
 465     } else {
 466         curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
 467         for (i = nsplits; i>0 ;i--) {
 468             nn=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
 469             new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(nn,tmpdigits,
 470                                                             Digits(),side);
 471             ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
 472                                                     SetRightNeighbour(kTRUE);
 473             if (curentcluster->GetRightNeighbour())
 474                 ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
 475                                                      SetRightNeighbour(kTRUE);
 476             else curentcluster->SetRightNeighbour(kTRUE);
 477             fNClusterN++;
 478         } // end loop over nplits
 479     } // end if side
 480     delete []tmpdigits;
 481 }
 482 //______________________________________________________________________
 483 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end){
 484     // sort digits on the P side
 485     Int_t right;
 486     Int_t left;
 487
 488     if (start != (end - 1) ){
 489         left=this->SortDigitsP(start,(start+end)/2);
 490         right=this->SortDigitsP((start+end)/2,end);
 491         return (left || right);
 492     }else{
 493         left =  ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[start]]))->
 494                                                               GetStripNumber();
 495         right= ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[end]]))->
 496                                                               GetStripNumber();
 497         if( left > right ){
 498             Int_t tmp = (*fDigitsIndexP)[start];
 499             (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
 500             (*fDigitsIndexP)[end]=tmp;
 501             return 1;
 502         }else return 0;
 503     } // end if
 504 }
 505 //______________________________________________________________________
 506 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end){
 507     // sort digits on the N side
 508     Int_t right;
 509     Int_t left;
 510
 511     if (start != (end - 1)){
 512         left=this->SortDigitsN(start,(start+end)/2);
 513         right=this->SortDigitsN((start+end)/2,end);
 514         return (left || right);
 515     }else{
 516         left =((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[start]]))->
 517                                                               GetStripNumber();
 518         right=((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[end]]))->
 519                                                               GetStripNumber();
 520         if ( left > right ){
 521             Int_t tmp = (*fDigitsIndexN)[start];
 522             (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
 523             (*fDigitsIndexN)[end]=tmp;
 524             return 1;
 525         }else return 0;
 526     } // end if
 527 }
 528 //______________________________________________________________________
 529 void AliITSClusterFinderSSD::FillDigitsIndex(){
 530     //Fill the indexes of the clusters belonging to a given ITS module
 531
 532     Int_t noentries;
 533     Int_t i;
 534     Int_t gain=0;
 535     Int_t signal=0;
 536
 537     noentries = fDigits->GetEntriesFast();
 538
 539     if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(noentries);
 540     if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(noentries);
 541
 542     AliITSdigitSSD *dig;
 543
 544     for ( i = 0 ; i< noentries; i++ ) {
 545
 546       dig = (AliITSdigitSSD*)GetDigit(i);
 547
 548       gain=(Int_t) ((AliITSCalibrationSSD*)GetResp(fModule))->GetGainP(dig->GetStripNumber());
 549       signal=gain*dig->GetSignal();
 550       dig->SetSignal(signal);
 551
 552       if(dig->IsSideP()) fDigitsIndexP->AddAt(i,fNDigitsP++);
 553       else fDigitsIndexN->AddAt(i,fNDigitsN++);
 554
 555     } // end for i
 556
 557     //    delete [] psidx;
 558     //delete [] nsidx;
 559
 560     if (debug) cout<<"Digits :  P = "<<fNDigitsP<<"   N = "<<fNDigitsN<<endl;
 561 }
 562 //______________________________________________________________________
 563 void AliITSClusterFinderSSD::SortDigits(){
 564     // sort digits
 565     Int_t i;
 566
 567     if(fNDigitsP>1) for (i=0;i<fNDigitsP-1;i++)
 568         if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
 569     if(fNDigitsN>1) for (i=0;i<fNDigitsN-1;i++)
 570         if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
 571 }
 572 //______________________________________________________________________
 573 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN) const{
 574     // fill cluster index array
 575     register Int_t i;
 576
 577     for (i=0; i<fNClusterP;i++) arrayP[i]=i;
 578     for (i=0; i<fNClusterN;i++) arrayN[i]=i;
 579 }
 580 //______________________________________________________________________
 581 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN){
 582     // sort clusters
 583     Int_t i;
 584
 585     if(fNClusterP>1) for (i=0;i<fNClusterP-1;i++)
 586         if (SortClustersP(0,(fNClusterP-1),arrayP)==0)  break;
 587   if(fNClusterN>1) for (i=0;i<fNClusterN-1;i++)
 588       if (SortClustersN(0,(fNClusterN-1),arrayN)==0)  break;
 589 }
 590 //______________________________________________________________________
 591 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end,
 592                                             Int_t *array){
 593     //Sort P side clusters
 594     Int_t right;
 595     Int_t left;
 596
 597     if (start != (end - 1) ) {
 598         left=this->SortClustersP(start,(start+end)/2,array);
 599         right=this->SortClustersP((start+end)/2,end,array);
 600         return (left || right);
 601     } else {
 602         left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
 603                                                          GetDigitStripNo(0);
 604         right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
 605                                                          GetDigitStripNo(0);
 606         if(left>right) {
 607             Int_t tmp = array[start];
 608             array[start]=array[end];
 609             array[end]=tmp;
 610             return 1;
 611         } else return 0;
 612     } // end if
 613 }
 614 //______________________________________________________________________
 615 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end,
 616                                             Int_t *array){
 617     //Sort N side clusters
 618     Int_t right;
 619     Int_t left;
 620
 621     if (start != (end - 1) ) {
 622         left=this->SortClustersN(start,(start+end)/2,array);
 623         right=this->SortClustersN((start+end)/2,end,array);
 624         return (left || right);
 625     } else {
 626         left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
 627                                                          GetDigitStripNo(0);
 628         right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
 629                                                          GetDigitStripNo(0);
 630         if( left > right) {
 631             Int_t tmp = array[start];
 632             array[start]=array[end];
 633             array[end]=tmp;
 634             return 1;
 635         } else return 0;
 636     } // end if
 637 }
 638 //______________________________________________________________________
 639 void AliITSClusterFinderSSD::ClustersToPackages(){
 640     // fill packages
 641
 642     Int_t *oneSclP = new Int_t[fNClusterP];//I want to have sorted 1 S clusters
 643     Int_t *oneSclN = new Int_t[fNClusterN];//I can not sort it in TClonesArray
 644                                            //so, I create table of indexes and
 645                                            //sort it
 646                                            //I do not use TArrayI on purpose
 647                                            //MB: well, that's not true that one
 648                                            //cannot sort objs in TClonesArray
 649     AliITSclusterSSD *currentP;
 650     AliITSclusterSSD *currentN;
 651     Int_t j1, j2;
 652
 653     //Fills in One Side Clusters Index Array
 654     FillClIndexArrays(oneSclP,oneSclN);
 655     //Sorts filled Arrays
 656     //SortClusters(oneSclP,oneSclN);
 657
 658     fNPackages=1;
 659     new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
 660
 661     //This part was includede by Boris Batiounia in March 2001.
 662     // Take all recpoint pairs (x coordinates) in both P and N sides
 663     // to calculate z coordinates of the recpoints
 664
 665     for (j1=0;j1<fNClusterP;j1++) {
 666         currentP = GetPSideCluster(oneSclP[j1]);
 667         Double_t xP = currentP->GetPosition();
 668         Float_t signalP = currentP->GetTotalSignal();
 669         for (j2=0;j2<fNClusterN;j2++) {
 670             currentN = GetNSideCluster(oneSclN[j2]);
 671             Double_t xN = currentN->GetPosition();
 672             Float_t signalN = currentN->GetTotalSignal();
 673             CreateNewRecPoint(xP,1,xN,1,signalP,signalN,currentP,currentN,
 674                               0.75);
 675         } // end for j2
 676     } // end for j1
 677
 678     delete [] oneSclP;
 679     delete [] oneSclN;
 680 }
 681 //______________________________________________________________________
 682 Bool_t AliITSClusterFinderSSD::CreateNewRecPoint(Float_t P,Float_t dP,
 683                                                  Float_t N, Float_t dN,
 684                                                  Float_t SigP,Float_t SigN,
 685                                                  AliITSclusterSSD *clusterP,
 686                                                  AliITSclusterSSD *clusterN,
 687                                                  Stat_t prob){
 688     // create the recpoints
 689     const Float_t kADCtoKeV = 2.16;
 690     // 50 ADC units -> 30000 e-h pairs; 1e-h pair -> 3.6e-3 KeV;
 691     // 1 ADC unit -> (30000/50)*3.6e-3 = 2.16 KeV
 692     const Float_t kconv = 1.0e-4;
 693     const Float_t kRMSx = 20.0*kconv;
 694     const Float_t kRMSz = 800.0*kconv;
 695     Int_t n=0;
 696     Int_t *tr;
 697     Int_t ntracks;
 698
 699     Int_t lay,lad,det;
 700     fDetTypeRec->GetITSgeom()->GetModuleId(fModule,lay,lad,det);
 701     Int_t ind=(lad-1)*fDetTypeRec->GetITSgeom()->GetNdetectors(lay)+(det-1);
 702     Int_t lyr=(lay-1);
 703
 704     if (GetCrossing(P,N)) {
 705         //GetCrossingError(dP,dN);
 706         dP = dN = prob = 0.0; // to remove unused variable warning.
 707         AliITSRawClusterSSD cnew;
 708         Int_t nstripsP=clusterP->GetNumOfDigits();
 709         Int_t nstripsN=clusterN->GetNumOfDigits();
 710         Float_t signal = 0;
 711         Float_t dedx = 0;
 712         //      Float_t meannoiseP=clusterP->GetMeanNoise();
 713         //Float_t meannoiseN=clusterN->GetMeanNoise();
 714         //cout<<meannoiseP<<" "<<meannoiseN<<endl;
 715         if(SigP>SigN) {
 716             signal = SigP;
 717             dedx = SigP*kADCtoKeV;
 718         }else{
 719             signal = SigN;
 720             dedx = SigN*kADCtoKeV;
 721         } // end if SigP>SigN
 722      tr = (Int_t*) clusterP->GetTracks(n);
 723      ntracks = clusterP->GetNTracks();
 724
 725      //cnew.SetDigitsClusterP(clusterP);
 726      //cnew.SetDigitsClusterN(clusterN);
 727
 728      cnew.SetSignalP(SigP);
 729      cnew.SetSignalN(SigN);
 730      cnew.SetMultiplicity(nstripsP);
 731      cnew.SetMultN(nstripsN);
 732      cnew.SetQErr(TMath::Abs(SigP-SigN));
 733      cnew.SetNTrack(ntracks);
 734      //cnew.SetMeanNoiseP(meannoiseP);
 735      //cnew.SetMeanNoiseN(meannoiseN);
 736        fDetTypeRec->AddCluster(2,&cnew);
 737        //fITS->AddCluster(2,&cnew);
 738        //AliITSRecPoint rnew;
 739         AliITSRecPoint rnew(fDetTypeRec->GetITSgeom());
 740         rnew.SetXZ(fModule,P*kconv,N*kconv);
 741         //rnew.SetX(P*kconv);
 742         //rnew.SetZ(N*kconv);
 743      rnew.SetQ(signal);
 744      rnew.SetdEdX(dedx);
 745         rnew.SetSigmaDetLocX2( kRMSx* kRMSx);
 746         //     rnew.SetSigmaX2( kRMSx* kRMSx);
 747     rnew.SetSigmaZ2( kRMSz* kRMSz);
 748
 749          rnew.SetLabel(tr[0],0);
 750         rnew.SetLabel(tr[1],1);
 751         rnew.SetLabel(tr[2],2);
 752         rnew.SetDetectorIndex(ind);
 753         rnew.SetLayer(lyr);
 754
 755     //rnew.fTracks[0]=tr[0];
 756     // rnew.fTracks[1]=tr[1];
 757     //rnew.fTracks[2]=tr[2];
 758      //rnew.SetMultP(nstripsP);
 759      //rnew.SetMultN(nstripsN);
 760         fDetTypeRec->AddRecPoint(rnew);
 761         //    fITS->AddRecPoint(rnew);
 762      return kTRUE;
 763     } // end if
 764     return kFALSE;
 765 }
 766 //______________________________________________________________________
 767 void  AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo){
 768     // calculate the step factor for matching clusters
 769     // 95 is the pitch, 4000 - dimension along z ?
 770     //Float_t dz=fSegmentation->Dz();
 771     Float_t dz=GetSeg()->Dz();
 772
 773     fSFF = ( (Int_t)  (Psteo*dz/fPitch ) );// +1;
 774     fSFB = ( (Int_t)  (Nsteo*dz/fPitch ) );// +1;
 775 }
 776 //______________________________________________________________________
 777 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx){
 778     // get P side clusters
 779
 780     if((idx<0)||(idx>=fNClusterP)){
 781         printf("AliITSClusterFinderSSD::GetPSideCluster: index out of range\n");
 782         return 0;
 783     }else{
 784         return (AliITSclusterSSD*)((*fClusterP)[idx]);
 785     } // end if
 786 }
 787 //______________________________________________________________________
 788 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx){
 789     // get N side clusters
 790
 791     if((idx<0)||(idx>=fNClusterN)){
 792         printf("AliITSClusterFinderSSD::GetNSideCluster: index out of range\n");
 793         return 0;
 794     }else{
 795         return (AliITSclusterSSD*)((*fClusterN)[idx]);
 796     } // end if
 797 }
 798 //______________________________________________________________________
 799 AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side){
 800     // Get cluster
 801
 802     return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
 803 }
 804 //______________________________________________________________________
 805 Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N){
 806     // get crossing
 807     // This function was rivised and changed by Boris Batiounia in March 2001
 808     Float_t dx = GetSeg()->Dx(); // detector size in x direction, microns
 809     Float_t dz = GetSeg()->Dz(); // detector size in z direction, microns
 810     //Float_t dx = fSegmentation->Dx(); // detector size in x direction, microns
 811     //Float_t dz = fSegmentation->Dz(); // detector size in z direction, microns
 812     //cout<<dx<<" "<<dz<<endl;
 813     Float_t xL; // x local coordinate
 814     Float_t zL; // z local coordinate
 815     Float_t x;  // x = xL + dx/2
 816     Float_t z;  // z = zL + dz/2
 817     Float_t xP; // x coordinate in the P side from the first P strip
 818     Float_t xN; // x coordinate in the N side from the first N strip
 819     Float_t stereoP,stereoN;
 820
 821     //fSegmentation->Angles(stereoP,stereoN);
 822     GetSeg()->Angles(stereoP,stereoN);
 823
 824     //cout<<stereoP<<" "<<stereoN<<" "<<P<<" "<<N<<endl;
 825
 826     //cout<<" P="<<P<<", N="<<N<<endl;
 827
 828     fTanP=TMath::Tan(stereoP);
 829     fTanN=TMath::Tan(stereoN);
 830     Float_t kP = fTanP; // Tangent of 0.0075 mrad
 831     Float_t kN = fTanN; // Tangent of 0.0275 mrad
 832     P *= fPitch;
 833     N *= fPitch;
 834
 835     xP = P;
 836     xN = N;
 837     //    xP = N;      // change the mistake for the P/N
 838     //xN = P;      // coordinates correspondence in this function
 839
 840     z=(xN-xP+dz*kN)/(kP+kN);
 841     x=xP+kP*z;
 842
 843     xL = x - dx/2;
 844     zL = z - dz/2;
 845     P = xL;
 846     N = zL;
 847     //cout<<"P= "<<P<<" , N= "<<N<<" , dx= "<<dx<<endl;
 848
 849     //cout<<"P="<<P<<", N="<<N<<endl;
 850
 851     if(TMath::Abs(xL) > dx/2 || TMath::Abs(zL) > dz/2) return kFALSE;
 852
 853     // Check that xL and zL are inside the detector for the
 854     // correspondent xP and xN coordinates
 855
 856     return kTRUE;
 857 }
 858 //______________________________________________________________________
 859 void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN){
 860     // get crossing error
 861     Float_t dz, dx;
 862
 863     dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
 864     dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
 865                  TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
 866     dN = dz;
 867     dP = dx;
 868 }