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     /*
 247     cout<<"----------------------------------------------------------------"<<endl;
 248     cout<<"module="<<module<<" , # of Pdigits="<<fNDigitsP<<" , # of Ndigits="<<fNDigitsN<<endl;
 249
 250     cout<<"  Pside"<<endl;
 251     cout<<"       "<<((AliITSdigitSSD*)lDigits[lDigitsIndexP[0]])->GetStripNumber()<<" "<<
 252       ((AliITSdigitSSD*)lDigits[lDigitsIndexP[0]])->GetSignal()<<endl;
 253     */
 254
 255     for (i=1; i<fNDigitsP; i++) {
 256
 257       //reads new digit
 258       currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->GetStripNumber();
 259       //      cout<<"       "<<currentstripNo<<" "<<((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->GetSignal()<<endl;
 260
 261       if((((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->GetStripNumber()) ==  (currentstripNo-1) )
 262         dbuffer[dnumber++]=lDigitsIndexP[i];
 263
 264       else{
 265
 266         // check signal
 267         for(j=0;j<dnumber;j++) {
 268
 269           if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 270               > fgkNoiseThreshold* ((AliITSCalibrationSSD*)GetResp(module))->
 271               GetNoiseP( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 272             flag+=1;
 273
 274         }
 275
 276         //if(flag==dnumber) {
 277         if(flag>0) {
 278           //create a new one side cluster
 279
 280           //      cout<<"          new cluster with "<<dnumber<<" digits"<<endl;
 281
 282           new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
 283                                                         Digits(),
 284                                                         fgkSIDEP);
 285
 286           flag=0;
 287         }
 288
 289         flag=0;
 290         dbuffer[0]=lDigitsIndexP[i];
 291         dnumber = 1;
 292
 293       } // end if else
 294
 295     } // end loop over fNDigitsP
 296
 297     // check signal
 298     for(j=0;j<dnumber;j++) {
 299
 300       if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 301           > fgkNoiseThreshold*((AliITSCalibrationSSD*)GetResp(module))->
 302           GetNoiseP( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 303         flag+=1;
 304
 305     }
 306
 307     //if(flag==dnumber) {
 308     if(flag>0) {
 309       //create a new one side cluster
 310
 311       //      cout<<"          new cluster with "<<dnumber<<" digits"<<endl;
 312
 313       new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,
 314                                                     Digits(),
 315                                                     fgkSIDEP);
 316
 317       flag=0;
 318     }
 319
 320     flag=0;
 321
 322     //process N side
 323     //for comments, see above
 324     dnumber = 1;
 325     dbuffer[0]=lDigitsIndexN[0];
 326     //If next digit is a neigh. of previous, adds to last clust. this digit
 327
 328     //    cout<<"  Nside"<<endl;
 329     //    cout<<"       "<<((AliITSdigitSSD*)lDigits[lDigitsIndexN[0]])->GetStripNumber()<<" "<<
 330     //      ((AliITSdigitSSD*)lDigits[lDigitsIndexN[0]])->GetSignal()<<endl;
 331
 332     for (i=1; i<fNDigitsN; i++) {
 333       currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->GetStripNumber();
 334       //      cout<<"       "<<currentstripNo<<" "<<((AliITSdigitSSD*)lDigits[lDigitsIndexN[i]])->GetSignal()<<endl;
 335
 336       if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->GetStripNumber()) == (currentstripNo-1) )
 337         dbuffer[dnumber++]=lDigitsIndexN[i];
 338
 339       else {
 340
 341         // check signal
 342         for(j=0;j<dnumber;j++) {
 343
 344           if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 345               > fgkNoiseThreshold*((AliITSCalibrationSSD*)GetResp(module))->
 346               GetNoiseN( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 347             flag+=1;
 348
 349         }
 350
 351         //if(flag==dnumber) {
 352         if(flag>0) {
 353           //create a new one side cluster
 354
 355           //      cout<<"          new cluster with "<<dnumber<<" digits"<<endl;
 356
 357           new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
 358                                                         Digits(),
 359                                                         fgkSIDEN);
 360
 361           flag=0;
 362         }
 363
 364         flag=0;
 365         dbuffer[0]=lDigitsIndexN[i];
 366         dnumber = 1;
 367       } // end if else
 368     } // end loop over fNDigitsN
 369
 370     // check signal
 371     for(j=0;j<dnumber;j++) {
 372
 373       if( (((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetSignal())
 374           > fgkNoiseThreshold*((AliITSCalibrationSSD*)GetResp(module))->
 375           GetNoiseN( ((AliITSdigitSSD*)lDigits[dbuffer[j]])->GetStripNumber()) )
 376         flag+=1;
 377
 378     }
 379
 380     //if(flag==dnumber) {
 381     if(flag>0) {
 382       //create a new one side cluster
 383
 384       //      cout<<"          new cluster with "<<dnumber<<" digits"<<endl;
 385
 386       new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,
 387                                                     Digits(),
 388                                                     fgkSIDEN);
 389
 390       flag=0;
 391     }
 392
 393     flag=0;
 394     delete [] dbuffer;
 395
 396     //    if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP
 397     //             <<"  fNClusterN ="<<fNClusterN<<"\n";
 398 }
 399 //______________________________________________________________________
 400 void AliITSClusterFinderSSD::SeparateOverlappedClusters(){
 401     // overlapped clusters separation
 402     register Int_t i; //iterator
 403     Float_t  factor=0.75;            // How many percent must be lower signal
 404                                      // on the middle one digit
 405                                      // from its neighbours
 406     Int_t    signal0;              //signal on the strip before the current one
 407     Int_t    signal1;              //signal on the current one signal
 408     Int_t    signal2;              //signal on the strip after the current one
 409     TArrayI *splitlist;              //  List of splits
 410     Int_t    numerofsplits=0;        // number of splits
 411     Int_t    initPsize = fNClusterP; //initial size of the arrays
 412     Int_t    initNsize = fNClusterN; //we have to keep it because it will grow
 413                                      // in this function and it doasn't make
 414                                      // sense to pass through it again
 415     splitlist = new TArrayI(800);
 416
 417     for (i=0;i<initPsize;i++){
 418         if (( ((AliITSclusterSSD*)(*fClusterP)[i])->
 419               GetNumOfDigits())==1) continue;
 420         if (( ((AliITSclusterSSD*)(*fClusterP)[i])->
 421               GetNumOfDigits())==2) continue;
 422         Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
 423         for (Int_t j=1; j<nj; j++){
 424             signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
 425             signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
 426             signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
 427             //if signal is less then factor*signal of its neighbours
 428             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) ){
 429                 (*splitlist)[numerofsplits++]=j;
 430             } // end if
 431         } // end loop over number of digits
 432         //split this cluster if necessary
 433         if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
 434         numerofsplits=0;
 435         //in signed places (splitlist)
 436     } // end loop over clusters on Pside
 437
 438     for (i=0;i<initNsize;i++) {
 439         if (( ((AliITSclusterSSD*)(*fClusterN)[i])->
 440               GetNumOfDigits())==1) continue;
 441         if (( ((AliITSclusterSSD*)(*fClusterN)[i])->
 442               GetNumOfDigits())==2) continue;
 443         Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
 444         for (Int_t j=1; j<nj; j++){
 445             signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
 446             signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
 447             signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
 448             //if signal is less then factor*signal of its neighbours
 449             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
 450                 (*splitlist)[numerofsplits++]=j;
 451         } // end loop over number of digits
 452         //split this cluster into more clusters
 453         if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
 454         numerofsplits=0;
 455         //in signed places (splitlist)
 456     } // end loop over clusters on Nside
 457
 458     delete splitlist;
 459 }
 460 //______________________________________________________________________
 461 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits,
 462                                           Int_t index, Bool_t side){
 463     //This function splits one side cluster into more clusters
 464     //number of splits is defined by "nsplits"
 465     //Place of splits are defined in the TArray "list"
 466     // For further optimisation: Replace this function by two
 467     // specialised ones (each for one side)
 468     // save one "if"
 469     //For comlete comments see AliITSclusterSSD::SplitCluster
 470     register Int_t i; //iterator
 471     AliITSclusterSSD* curentcluster;
 472     Int_t   *tmpdigits = new Int_t[100];
 473     Int_t    nn;
 474
 475     // side true means P side
 476     if (side) {
 477         curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
 478         for (i = nsplits; i>0 ;i--) {
 479             nn=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
 480             new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(nn,tmpdigits,
 481                                                             Digits(),side);
 482             ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
 483                                                       SetLeftNeighbour(kTRUE);
 484             //if left cluster had neighbour on the right before split
 485             //new should have it too
 486             if ( curentcluster->GetRightNeighbour() )
 487                 ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
 488                                                      SetRightNeighbour(kTRUE);
 489             else curentcluster->SetRightNeighbour(kTRUE);
 490             fNClusterP++;
 491         } // end loop over nplits
 492     } else {
 493         curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
 494         for (i = nsplits; i>0 ;i--) {
 495             nn=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
 496             new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(nn,tmpdigits,
 497                                                             Digits(),side);
 498             ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
 499                                                     SetRightNeighbour(kTRUE);
 500             if (curentcluster->GetRightNeighbour())
 501                 ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
 502                                                      SetRightNeighbour(kTRUE);
 503             else curentcluster->SetRightNeighbour(kTRUE);
 504             fNClusterN++;
 505         } // end loop over nplits
 506     } // end if side
 507     delete []tmpdigits;
 508 }
 509 //______________________________________________________________________
 510 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end){
 511     // sort digits on the P side
 512     Int_t right;
 513     Int_t left;
 514
 515     if (start != (end - 1) ){
 516         left=this->SortDigitsP(start,(start+end)/2);
 517         right=this->SortDigitsP((start+end)/2,end);
 518         return (left || right);
 519     }else{
 520         left =  ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[start]]))->
 521                                                               GetStripNumber();
 522         right= ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[end]]))->
 523                                                               GetStripNumber();
 524         if( left > right ){
 525             Int_t tmp = (*fDigitsIndexP)[start];
 526             (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
 527             (*fDigitsIndexP)[end]=tmp;
 528             return 1;
 529         }else return 0;
 530     } // end if
 531 }
 532 //______________________________________________________________________
 533 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end){
 534     // sort digits on the N side
 535     Int_t right;
 536     Int_t left;
 537
 538     if (start != (end - 1)){
 539         left=this->SortDigitsN(start,(start+end)/2);
 540         right=this->SortDigitsN((start+end)/2,end);
 541         return (left || right);
 542     }else{
 543         left =((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[start]]))->
 544                                                               GetStripNumber();
 545         right=((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[end]]))->
 546                                                               GetStripNumber();
 547         if ( left > right ){
 548             Int_t tmp = (*fDigitsIndexN)[start];
 549             (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
 550             (*fDigitsIndexN)[end]=tmp;
 551             return 1;
 552         }else return 0;
 553     } // end if
 554 }
 555 //______________________________________________________________________
 556 void AliITSClusterFinderSSD::FillDigitsIndex(){
 557     //Fill the indexes of the clusters belonging to a given ITS module
 558
 559     Int_t noentries;
 560     Int_t i;
 561     Int_t gain=0;
 562     Int_t signal=0;
 563
 564     noentries = fDigits->GetEntriesFast();
 565
 566     if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(noentries);
 567     if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(noentries);
 568
 569     AliITSdigitSSD *dig;
 570
 571     for ( i = 0 ; i< noentries; i++ ) {
 572
 573       dig = (AliITSdigitSSD*)GetDigit(i);
 574
 575       gain=(Int_t) ((AliITSCalibrationSSD*)GetResp(fModule))->GetGainP(dig->GetStripNumber());
 576       signal=gain*dig->GetSignal();
 577       dig->SetSignal(signal);
 578
 579       if(dig->IsSideP()) fDigitsIndexP->AddAt(i,fNDigitsP++);
 580       else fDigitsIndexN->AddAt(i,fNDigitsN++);
 581
 582     } // end for i
 583
 584     //    delete [] psidx;
 585     //delete [] nsidx;
 586
 587     //    if (debug) cout<<"Digits :  P = "<<fNDigitsP<<"   N = "<<fNDigitsN<<endl;
 588 }
 589 //______________________________________________________________________
 590 void AliITSClusterFinderSSD::SortDigits(){
 591     // sort digits
 592     Int_t i;
 593
 594     if(fNDigitsP>1) for (i=0;i<fNDigitsP-1;i++)
 595         if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
 596     if(fNDigitsN>1) for (i=0;i<fNDigitsN-1;i++)
 597         if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
 598 }
 599 //______________________________________________________________________
 600 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN) const{
 601     // fill cluster index array
 602     register Int_t i;
 603
 604     for (i=0; i<fNClusterP;i++) arrayP[i]=i;
 605     for (i=0; i<fNClusterN;i++) arrayN[i]=i;
 606 }
 607 //______________________________________________________________________
 608 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN){
 609     // sort clusters
 610     Int_t i;
 611
 612     if(fNClusterP>1) for (i=0;i<fNClusterP-1;i++)
 613         if (SortClustersP(0,(fNClusterP-1),arrayP)==0)  break;
 614   if(fNClusterN>1) for (i=0;i<fNClusterN-1;i++)
 615       if (SortClustersN(0,(fNClusterN-1),arrayN)==0)  break;
 616 }
 617 //______________________________________________________________________
 618 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end,
 619                                             Int_t *array){
 620     //Sort P side clusters
 621     Int_t right;
 622     Int_t left;
 623
 624     if (start != (end - 1) ) {
 625         left=this->SortClustersP(start,(start+end)/2,array);
 626         right=this->SortClustersP((start+end)/2,end,array);
 627         return (left || right);
 628     } else {
 629         left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
 630                                                          GetDigitStripNo(0);
 631         right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
 632                                                          GetDigitStripNo(0);
 633         if(left>right) {
 634             Int_t tmp = array[start];
 635             array[start]=array[end];
 636             array[end]=tmp;
 637             return 1;
 638         } else return 0;
 639     } // end if
 640 }
 641 //______________________________________________________________________
 642 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end,
 643                                             Int_t *array){
 644     //Sort N side clusters
 645     Int_t right;
 646     Int_t left;
 647
 648     if (start != (end - 1) ) {
 649         left=this->SortClustersN(start,(start+end)/2,array);
 650         right=this->SortClustersN((start+end)/2,end,array);
 651         return (left || right);
 652     } else {
 653         left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
 654                                                          GetDigitStripNo(0);
 655         right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
 656                                                          GetDigitStripNo(0);
 657         if( left > right) {
 658             Int_t tmp = array[start];
 659             array[start]=array[end];
 660             array[end]=tmp;
 661             return 1;
 662         } else return 0;
 663     } // end if
 664 }
 665 //______________________________________________________________________
 666 void AliITSClusterFinderSSD::ClustersToPackages(){
 667     // fill packages
 668
 669     Int_t *oneSclP = new Int_t[fNClusterP];//I want to have sorted 1 S clusters
 670     Int_t *oneSclN = new Int_t[fNClusterN];//I can not sort it in TClonesArray
 671                                            //so, I create table of indexes and
 672                                            //sort it
 673                                            //I do not use TArrayI on purpose
 674                                            //MB: well, that's not true that one
 675                                            //cannot sort objs in TClonesArray
 676     AliITSclusterSSD *currentP;
 677     AliITSclusterSSD *currentN;
 678     Int_t j1, j2;
 679
 680     //Fills in One Side Clusters Index Array
 681     FillClIndexArrays(oneSclP,oneSclN);
 682     //Sorts filled Arrays
 683     //SortClusters(oneSclP,oneSclN);
 684
 685     fNPackages=1;
 686     new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
 687
 688     //This part was includede by Boris Batiounia in March 2001.
 689     // Take all recpoint pairs (x coordinates) in both P and N sides
 690     // to calculate z coordinates of the recpoints
 691
 692     for (j1=0;j1<fNClusterP;j1++) {
 693         currentP = GetPSideCluster(oneSclP[j1]);
 694         Double_t xP = currentP->GetPosition();
 695         Float_t signalP = currentP->GetTotalSignal();
 696         for (j2=0;j2<fNClusterN;j2++) {
 697             currentN = GetNSideCluster(oneSclN[j2]);
 698             Double_t xN = currentN->GetPosition();
 699             Float_t signalN = currentN->GetTotalSignal();
 700             CreateNewRecPoint(xP,1,xN,1,signalP,signalN,currentP,currentN,
 701                               0.75);
 702         } // end for j2
 703     } // end for j1
 704
 705     delete [] oneSclP;
 706     delete [] oneSclN;
 707 }
 708 //______________________________________________________________________
 709 Bool_t AliITSClusterFinderSSD::CreateNewRecPoint(Float_t P,Float_t dP,
 710                                                  Float_t N, Float_t dN,
 711                                                  Float_t SigP,Float_t SigN,
 712                                                  AliITSclusterSSD *clusterP,
 713                                                  AliITSclusterSSD *clusterN,
 714                                                  Stat_t prob){
 715     // create the recpoints
 716     const Float_t kADCtoKeV = 2.16;
 717     // 50 ADC units -> 30000 e-h pairs; 1e-h pair -> 3.6e-3 KeV;
 718     // 1 ADC unit -> (30000/50)*3.6e-3 = 2.16 KeV
 719     const Float_t kconv = 1.0e-4;
 720     const Float_t kRMSx = 20.0*kconv;
 721     const Float_t kRMSz = 800.0*kconv;
 722     Int_t n=0;
 723     Int_t *tr;
 724     Int_t ntracks;
 725
 726     Int_t lay,lad,det;
 727     fDetTypeRec->GetITSgeom()->GetModuleId(fModule,lay,lad,det);
 728     Int_t ind=(lad-1)*fDetTypeRec->GetITSgeom()->GetNdetectors(lay)+(det-1);
 729     Int_t lyr=(lay-1);
 730
 731     if (GetCrossing(P,N)) {
 732         //GetCrossingError(dP,dN);
 733         dP = dN = prob = 0.0; // to remove unused variable warning.
 734         AliITSRawClusterSSD cnew;
 735         Int_t nstripsP=clusterP->GetNumOfDigits();
 736         Int_t nstripsN=clusterN->GetNumOfDigits();
 737         Float_t signal = 0;
 738         Float_t dedx = 0;
 739         //      Float_t meannoiseP=clusterP->GetMeanNoise();
 740         //Float_t meannoiseN=clusterN->GetMeanNoise();
 741         //cout<<meannoiseP<<" "<<meannoiseN<<endl;
 742         if(SigP>SigN) {
 743             signal = SigP;
 744             dedx = SigP*kADCtoKeV;
 745         }else{
 746             signal = SigN;
 747             dedx = SigN*kADCtoKeV;
 748         } // end if SigP>SigN
 749      tr = (Int_t*) clusterP->GetTracks(n);
 750      ntracks = clusterP->GetNTracks();
 751
 752      //cnew.SetDigitsClusterP(clusterP);
 753      //cnew.SetDigitsClusterN(clusterN);
 754
 755      cnew.SetSignalP(SigP);
 756      cnew.SetSignalN(SigN);
 757      cnew.SetMultiplicity(nstripsP);
 758      cnew.SetMultN(nstripsN);
 759      cnew.SetQErr(TMath::Abs(SigP-SigN));
 760      cnew.SetNTrack(ntracks);
 761      //cnew.SetMeanNoiseP(meannoiseP);
 762      //cnew.SetMeanNoiseN(meannoiseN);
 763        fDetTypeRec->AddCluster(2,&cnew);
 764        //fITS->AddCluster(2,&cnew);
 765        //AliITSRecPoint rnew;
 766         AliITSRecPoint rnew(fDetTypeRec->GetITSgeom());
 767         rnew.SetXZ(fModule,P*kconv,N*kconv);
 768         //rnew.SetX(P*kconv);
 769         //rnew.SetZ(N*kconv);
 770      rnew.SetQ(signal);
 771      rnew.SetdEdX(dedx);
 772         rnew.SetSigmaDetLocX2( kRMSx* kRMSx);
 773         //     rnew.SetSigmaX2( kRMSx* kRMSx);
 774     rnew.SetSigmaZ2( kRMSz* kRMSz);
 775
 776          rnew.SetLabel(tr[0],0);
 777         rnew.SetLabel(tr[1],1);
 778         rnew.SetLabel(tr[2],2);
 779         rnew.SetDetectorIndex(ind);
 780         rnew.SetLayer(lyr);
 781
 782         rnew.SetNy(nstripsP);
 783         rnew.SetNz(nstripsN);
 784
 785     //rnew.fTracks[0]=tr[0];
 786     // rnew.fTracks[1]=tr[1];
 787     //rnew.fTracks[2]=tr[2];
 788      //rnew.SetMultP(nstripsP);
 789      //rnew.SetMultN(nstripsN);
 790         fDetTypeRec->AddRecPoint(rnew);
 791         //    fITS->AddRecPoint(rnew);
 792      return kTRUE;
 793     } // end if
 794     return kFALSE;
 795 }
 796 //______________________________________________________________________
 797 void  AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo){
 798     // calculate the step factor for matching clusters
 799     // 95 is the pitch, 4000 - dimension along z ?
 800     //Float_t dz=fSegmentation->Dz();
 801     Float_t dz=GetSeg()->Dz();
 802
 803     fSFF = ( (Int_t)  (Psteo*dz/fPitch ) );// +1;
 804     fSFB = ( (Int_t)  (Nsteo*dz/fPitch ) );// +1;
 805 }
 806 //______________________________________________________________________
 807 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx){
 808     // get P side clusters
 809
 810     if((idx<0)||(idx>=fNClusterP)){
 811         printf("AliITSClusterFinderSSD::GetPSideCluster: index out of range\n");
 812         return 0;
 813     }else{
 814         return (AliITSclusterSSD*)((*fClusterP)[idx]);
 815     } // end if
 816 }
 817 //______________________________________________________________________
 818 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx){
 819     // get N side clusters
 820
 821     if((idx<0)||(idx>=fNClusterN)){
 822         printf("AliITSClusterFinderSSD::GetNSideCluster: index out of range\n");
 823         return 0;
 824     }else{
 825         return (AliITSclusterSSD*)((*fClusterN)[idx]);
 826     } // end if
 827 }
 828 //______________________________________________________________________
 829 AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side){
 830     // Get cluster
 831
 832     return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
 833 }
 834 //______________________________________________________________________
 835 Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N){
 836     // get crossing
 837     // This function was rivised and changed by Boris Batiounia in March 2001
 838     Float_t dx = GetSeg()->Dx(); // detector size in x direction, microns
 839     Float_t dz = GetSeg()->Dz(); // detector size in z direction, microns
 840     //Float_t dx = fSegmentation->Dx(); // detector size in x direction, microns
 841     //Float_t dz = fSegmentation->Dz(); // detector size in z direction, microns
 842     //cout<<dx<<" "<<dz<<endl;
 843     Float_t xL; // x local coordinate
 844     Float_t zL; // z local coordinate
 845     Float_t x;  // x = xL + dx/2
 846     Float_t z;  // z = zL + dz/2
 847     Float_t xP; // x coordinate in the P side from the first P strip
 848     Float_t xN; // x coordinate in the N side from the first N strip
 849     Float_t stereoP,stereoN;
 850
 851     //fSegmentation->Angles(stereoP,stereoN);
 852     GetSeg()->Angles(stereoP,stereoN);
 853
 854     //cout<<stereoP<<" "<<stereoN<<" "<<P<<" "<<N<<endl;
 855
 856     //cout<<" P="<<P<<", N="<<N<<endl;
 857
 858     fTanP=TMath::Tan(stereoP);
 859     fTanN=TMath::Tan(stereoN);
 860     Float_t kP = fTanP; // Tangent of 0.0075 mrad
 861     Float_t kN = fTanN; // Tangent of 0.0275 mrad
 862     P *= fPitch;
 863     N *= fPitch;
 864
 865     xP = P;
 866     xN = N;
 867     //    xP = N;      // change the mistake for the P/N
 868     //xN = P;      // coordinates correspondence in this function
 869
 870     z=(xN-xP+dz*kN)/(kP+kN);
 871     x=xP+kP*z;
 872
 873     xL = x - dx/2;
 874     zL = z - dz/2;
 875     P = xL;
 876     N = zL;
 877     //cout<<"P= "<<P<<" , N= "<<N<<" , dx= "<<dx<<endl;
 878
 879     //cout<<"P="<<P<<", N="<<N<<endl;
 880
 881     if(TMath::Abs(xL) > dx/2 || TMath::Abs(zL) > dz/2) return kFALSE;
 882
 883     // Check that xL and zL are inside the detector for the
 884     // correspondent xP and xN coordinates
 885
 886     return kTRUE;
 887 }
 888 //______________________________________________________________________
 889 void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN){
 890     // get crossing error
 891     Float_t dz, dx;
 892
 893     dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
 894     dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
 895                  TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
 896     dN = dz;
 897     dP = dx;
 898 }