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
  18 Adding rekonstruction facilities
  19 Piotr Krzysztof Skowronski
  20 December 1999.
  21
  22 //Piotr Krzysztof Skowronski
  23 //Warsaw University of Technology
  24 //skowron@if.pw.edu.pl
  25
  26 */
  27
  28 /*
  29 15 -18 V 2000
  30 Eroor counting routines
  31 Automatic combination routines improved (traps)
  32
  33 */
  34
  35 #include <iostream.h>
  36 #include <TArrayI.h>
  37
  38
  39 #include "AliRun.h"
  40 #include "AliITS.h"
  41 #include "AliITSdigit.h"
  42 #include "AliITSRawCluster.h"
  43 #include "AliITSRecPoint.h"
  44 #include "AliITSMapA1.h"
  45 #include "AliITSClusterFinderSSD.h"
  46 #include "AliITSclusterSSD.h"
  47 #include "AliITSpackageSSD.h"
  48 #include "AliITSsegmentation.h"
  49
  50 const Bool_t AliITSClusterFinderSSD::fgkSIDEP=kTRUE;
  51 const Bool_t AliITSClusterFinderSSD::fgkSIDEN=kFALSE;
  52
  53 const Int_t debug=0;
  54
  55 ClassImp(AliITSClusterFinderSSD)
  56
  57 //____________________________________________________________________
  58 //
  59 //  Constructor
  60 //____________________________________________________________________
  61 //
  62
  63
  64 AliITSClusterFinderSSD::AliITSClusterFinderSSD(AliITSsegmentation *seg, TClonesArray *digits)
  65 {
  66 //Standard constructor
  67
  68     fSegmentation=seg;
  69     fDigits=digits;
  70
  71     fMap = new AliITSMapA1(fSegmentation,fDigits);
  72
  73     fITS=(AliITS*)gAlice->GetModule("ITS");
  74
  75     fClusterP =  new TClonesArray ("AliITSclusterSSD",200);
  76     fNClusterP =0;
  77
  78     fClusterN=  new TClonesArray ("AliITSclusterSSD",200);
  79     fNClusterN =0;
  80
  81     fPackages =  new TClonesArray ("AliITSpackageSSD",200);    //packages
  82     fNPackages = 0;
  83
  84
  85     fDigitsIndexP      =  new TArrayI(300);
  86     fNDigitsP          =  0;
  87
  88     fDigitsIndexN      =  new TArrayI(300);
  89     fNDigitsN          =  0;
  90
  91     SetAlpha1(1000);
  92     SetAlpha2(1000);
  93     SetAlpha3(1000);
  94
  95
  96     fPitch = fSegmentation->Dpx(0);
  97     Float_t StereoP,StereoN;
  98     fSegmentation->Angles(StereoP,StereoN);
  99     fTanP=TMath::Tan(StereoP);
 100     fTanN=TMath::Tan(StereoN);
 101
 102     fPNsignalRatio=7./8.;         // warning: hard-wired number
 103
 104 }
 105
 106 //-------------------------------------------------------
 107 AliITSClusterFinderSSD::~AliITSClusterFinderSSD() {
 108 // Default destructor
 109
 110     delete fClusterP;
 111     delete fClusterN;
 112     delete fPackages;
 113     delete fDigitsIndexP;
 114     delete fDigitsIndexN;
 115
 116     delete fMap;
 117
 118 }
 119
 120 //-------------------------------------------------------
 121 void AliITSClusterFinderSSD::InitReconstruction()
 122 {
 123 // initialization of the cluster finder
 124
 125   register Int_t i; //iterator
 126
 127   for (i=0;i<fNClusterP;i++)
 128     {
 129       fClusterP->RemoveAt(i);
 130     }
 131   fNClusterP  =0;
 132   for (i=0;i<fNClusterN;i++)
 133     {
 134       fClusterN->RemoveAt(i);
 135     }
 136   fNClusterN=0;
 137
 138   for (i=0;i<fNPackages;i++)
 139     {
 140       fPackages->RemoveAt(i);
 141     }
 142
 143   fNPackages = 0;
 144   fNDigitsP=0;
 145   fNDigitsN=0;
 146
 147   Float_t StereoP,StereoN;
 148   fSegmentation->Angles(StereoP,StereoN);
 149
 150   CalcStepFactor (StereoP,StereoN);
 151
 152   if (debug) cout<<"fSFF = "<<fSFF<<"  fSFB = "<<fSFB<<"\n";
 153 }
 154
 155
 156 //---------------------------------------------
 157 void AliITSClusterFinderSSD::FindRawClusters()
 158 {
 159 // This function findes out all clusters belonging to one module
 160 // 1. Zeroes all space after previous module reconstruction
 161 // 2. Finds all neighbouring digits
 162 // 3. If necesery, resolves for each group of neighbouring digits
 163 //    how many clusters creates it.
 164 // 4. Creates packages
 165 // 5. Creates clusters
 166
 167   InitReconstruction();  //ad. 1
 168   fMap->FillMap();
 169   FillDigitsIndex();
 170   SortDigits();
 171   FindNeighbouringDigits(); //ad. 2
 172   SeparateOverlappedClusters();  //ad. 3
 173   ClustersToPackages();  //ad. 4
 174   ConsumeClusters();
 175   PackagesToPoints();   //ad. 5
 176   ReconstructNotConsumedClusters();
 177
 178   fMap->ClearMap();
 179 }
 180
 181
 182 //-------------------------------------------------
 183 void AliITSClusterFinderSSD::FindNeighbouringDigits()
 184 {
 185 //If there are any digits on this side, create 1st Cluster,
 186 // add to it this digit, and increment number of clusters
 187
 188
 189  register Int_t i;
 190
 191
 192  if ((fNDigitsP>0 )  && (fNDigitsN > 0 )) {
 193
 194    Int_t currentstripNo;
 195    Int_t *dbuffer = new Int_t [300];   //buffer for strip numbers
 196    Int_t dnumber;    //curent number of digits in buffer
 197    TArrayI      &lDigitsIndexP = *fDigitsIndexP;
 198    TArrayI      &lDigitsIndexN = *fDigitsIndexN;
 199    TObjArray    &lDigits       = *(Digits());
 200    TClonesArray &lClusterP     = *fClusterP;
 201    TClonesArray &lClusterN     = *fClusterN;
 202
 203    //process P side
 204    dnumber = 1;
 205    dbuffer[0]=lDigitsIndexP[0];
 206    //If next digit is a neighbour of previous, adds to last cluster this digit
 207    for (i=1; i<fNDigitsP; i++) {
 208      //reads new digit
 209      currentstripNo = ((AliITSdigitSSD*)lDigits[lDigitsIndexP[i]])->
 210                                                             GetStripNumber();
 211      //check if it is a neighbour of a previous one
 212      if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexP[i-1]])->GetStripNumber())
 213            ==  (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexP[i];
 214      else  {
 215        //create a new one side cluster
 216        new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEP);
 217        dbuffer[0]=lDigitsIndexP[i];
 218        dnumber = 1;
 219      }
 220    } // end loop over fNDigitsP
 221    new(lClusterP[fNClusterP++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEP);
 222
 223
 224    //process N side
 225    //for comments, see above
 226    dnumber = 1;
 227    dbuffer[0]=lDigitsIndexN[0];
 228    //If next digit is a neighbour of previous, adds to last cluster this digit
 229    for (i=1; i<fNDigitsN; i++) {
 230      currentstripNo = ((AliITSdigitSSD*)(lDigits[lDigitsIndexN[i]]))->
 231                                                             GetStripNumber();
 232      if ( (((AliITSdigitSSD*)lDigits[lDigitsIndexN[i-1]])->GetStripNumber())
 233            == (currentstripNo-1) ) dbuffer[dnumber++]=lDigitsIndexN[i];
 234      else {
 235        new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEN);
 236        dbuffer[0]=lDigitsIndexN[i];
 237        dnumber = 1;
 238      }
 239    } // end loop over fNDigitsN
 240    new(lClusterN[fNClusterN++]) AliITSclusterSSD(dnumber,dbuffer,Digits(),fgkSIDEN);
 241    delete [] dbuffer;
 242
 243  } // end condition on  NDigits
 244
 245  if (debug) cout<<"\n Found clusters: fNClusterP = "<<fNClusterP<<"  fNClusterN ="<<fNClusterN<<"\n";
 246
 247 }
 248 /**********************************************************************/
 249
 250
 251 void AliITSClusterFinderSSD::SeparateOverlappedClusters()
 252 {
 253 // overlapped clusters separation
 254
 255   register Int_t i; //iterator
 256
 257   Float_t  factor=0.75;            // How many percent must be lower signal
 258                                    // on the middle one digit
 259                                    // from its neighbours
 260   Int_t    signal0;                //signal on the strip before the current one
 261   Int_t    signal1;                //signal on the current one signal
 262   Int_t    signal2;                //signal on the strip after the current one
 263   TArrayI *splitlist;              //  List of splits
 264   Int_t    numerofsplits=0;        // number of splits
 265   Int_t    initPsize = fNClusterP; //initial size of the arrays
 266   Int_t    initNsize = fNClusterN; //we have to keep it because it will grow
 267                                    // in this function and it doasn't make
 268                                    // sense to pass through it again
 269
 270   splitlist = new TArrayI(300);
 271
 272   for (i=0;i<initPsize;i++)
 273   {
 274     if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==1) continue;
 275     if (( ((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits())==2) continue;
 276         Int_t nj=(((AliITSclusterSSD*)(*fClusterP)[i])->GetNumOfDigits()-1);
 277         for (Int_t j=1; j<nj; j++)
 278           {
 279             signal1=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j);
 280             signal0=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j-1);
 281             signal2=((AliITSclusterSSD*)(*fClusterP)[i])->GetDigitSignal(j+1);
 282             //if signal is less then factor*signal of its neighbours
 283             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
 284              {
 285                (*splitlist)[numerofsplits++]=j;
 286              }
 287           } // end loop over number of digits
 288           //split this cluster if necessary
 289           if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEP);
 290           numerofsplits=0;
 291
 292           //in signed places (splitlist)
 293   } // end loop over clusters on Pside
 294
 295   for (i=0;i<initNsize;i++) {
 296     if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==1) continue;
 297     if (( ((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits())==2) continue;
 298        Int_t nj=(((AliITSclusterSSD*)(*fClusterN)[i])->GetNumOfDigits()-1);
 299        for (Int_t j=1; j<nj; j++)
 300           {
 301             signal1=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j);
 302             signal0=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j-1);
 303             signal2=((AliITSclusterSSD*)(*fClusterN)[i])->GetDigitSignal(j+1);
 304             //if signal is less then factor*signal of its neighbours
 305             if (  (signal1<(factor*signal0)) && (signal1<(factor*signal2)) )
 306                (*splitlist)[numerofsplits++]=j;
 307           } // end loop over number of digits
 308           //split this cluster into more clusters
 309           if(numerofsplits>0) SplitCluster(splitlist,numerofsplits,i,fgkSIDEN);
 310           numerofsplits=0;                                                              //in signed places (splitlist)
 311   } // end loop over clusters on Nside
 312
 313   delete splitlist;
 314 }
 315
 316 //-------------------------------------------------------
 317 void AliITSClusterFinderSSD::SplitCluster(TArrayI *list, Int_t nsplits, Int_t index, Bool_t side)
 318 {
 319   //This function splits one side cluster into more clusters
 320   //number of splits is defined by "nsplits"
 321   //Place of splits are defined in the TArray "list"
 322
 323   // For further optimisation: Replace this function by two
 324   // specialised ones (each for one side)
 325   // save one "if"
 326
 327   //For comlete comments see AliITSclusterSSD::SplitCluster
 328
 329
 330   register Int_t i; //iterator
 331
 332   AliITSclusterSSD* curentcluster;
 333   Int_t   *tmpdigits = new Int_t[100];
 334   Int_t    NN;
 335   // side true means P side
 336   if (side) {
 337      curentcluster =((AliITSclusterSSD*)((*fClusterP)[index])) ;
 338      for (i = nsplits; i>0 ;i--) {
 339          NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
 340          new ((*fClusterP)[fNClusterP]) AliITSclusterSSD(NN,tmpdigits,Digits(),side);
 341          ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
 342                                                       SetLeftNeighbour(kTRUE);
 343          //if left cluster had neighbour on the right before split
 344          //new should have it too
 345          if ( curentcluster->GetRightNeighbour() )
 346                      ( (AliITSclusterSSD*)((*fClusterP)[fNClusterP]) )->
 347                                                      SetRightNeighbour(kTRUE);
 348          else curentcluster->SetRightNeighbour(kTRUE);
 349          fNClusterP++;
 350      } // end loop over nplits
 351   } else {
 352      curentcluster =((AliITSclusterSSD*)((*fClusterN)[index]));
 353      for (i = nsplits; i>0 ;i--) {
 354          NN=curentcluster->SplitCluster((*list)[(i-1)],tmpdigits);
 355          new ((*fClusterN)[fNClusterN]) AliITSclusterSSD(NN,tmpdigits,Digits(),side);
 356          ((AliITSclusterSSD*)((*fClusterN)[fNClusterN]))->
 357                                                     SetRightNeighbour(kTRUE);
 358          if (curentcluster->GetRightNeighbour())
 359                       ( (AliITSclusterSSD*)( (*fClusterN)[fNClusterN]) )->
 360                                                      SetRightNeighbour(kTRUE);
 361          else curentcluster->SetRightNeighbour(kTRUE);
 362          fNClusterN++;
 363      } // end loop over nplits
 364   } // end if side
 365   delete []tmpdigits;
 366
 367 }
 368
 369
 370 //-------------------------------------------------
 371 Int_t AliITSClusterFinderSSD::SortDigitsP(Int_t start, Int_t end)
 372 {
 373 // sort digits on the P side
 374
 375
 376   Int_t right;
 377   Int_t left;
 378   if (start != (end - 1) )
 379     {
 380       left=this->SortDigitsP(start,(start+end)/2);
 381       right=this->SortDigitsP((start+end)/2,end);
 382       return (left || right);
 383     }
 384   else
 385    {
 386     left =  ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[start]]))->GetStripNumber();
 387     right= ((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexP)[end]]))->GetStripNumber();
 388     if( left > right )
 389      {
 390        Int_t tmp = (*fDigitsIndexP)[start];
 391        (*fDigitsIndexP)[start]=(*fDigitsIndexP)[end];
 392        (*fDigitsIndexP)[end]=tmp;
 393        return 1;
 394      }
 395     else return 0;
 396    }
 397 }
 398
 399
 400 //--------------------------------------------------
 401
 402 Int_t AliITSClusterFinderSSD::SortDigitsN(Int_t start, Int_t end)
 403 {
 404 // sort digits on the N side
 405
 406   Int_t right;
 407   Int_t left;
 408   if (start != (end - 1))
 409     {
 410       left=this->SortDigitsN(start,(start+end)/2);
 411       right=this->SortDigitsN((start+end)/2,end);
 412       return (left || right);
 413     }
 414   else
 415    {
 416     left =((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[start]]))->GetStripNumber();
 417     right=((AliITSdigitSSD*)((*(Digits()))[(*fDigitsIndexN)[end]]))->GetStripNumber();
 418     if ( left > right )
 419       {
 420         Int_t tmp = (*fDigitsIndexN)[start];
 421         (*fDigitsIndexN)[start]=(*fDigitsIndexN)[end];
 422         (*fDigitsIndexN)[end]=tmp;
 423         return 1;
 424       }else return 0;
 425    }
 426 }
 427
 428
 429 //------------------------------------------------
 430 void AliITSClusterFinderSSD::FillDigitsIndex()
 431 {
 432  //Fill the indexes of the clusters belonging to a given ITS module
 433
 434  Int_t PNs=0, NNs=0;
 435  Int_t tmp,bit,k;
 436  Int_t N;
 437  Int_t i;
 438
 439  N = fDigits->GetEntriesFast();
 440
 441  Int_t* PSidx = new Int_t [N*sizeof(Int_t)];
 442  Int_t* NSidx = new Int_t [N*sizeof(Int_t)];
 443  if (fDigitsIndexP==NULL) fDigitsIndexP = new TArrayI(N);
 444  if (fDigitsIndexN==NULL) fDigitsIndexN = new TArrayI(N);
 445
 446  AliITSdigitSSD *dig;
 447
 448  for ( i = 0 ; i< N; i++ ) {
 449       dig = (AliITSdigitSSD*)GetDigit(i);
 450       if(dig->IsSideP()) {
 451            bit=1;
 452            tmp=dig->GetStripNumber();
 453            // I find this totally unnecessary - it's just a
 454            // CPU consuming double check
 455            for( k=0;k<PNs;k++)
 456             {
 457              if (tmp==PSidx[k])
 458               {
 459                  if (debug) cout<<"Such a digit exists \n";
 460                  bit=0;
 461               }
 462            }
 463            // end comment
 464         if(bit) {
 465             fDigitsIndexP->AddAt(i,fNDigitsP++);
 466             PSidx[PNs++]=tmp;
 467         }
 468       } else {
 469          bit=1;
 470          tmp=dig->GetStripNumber();
 471          // same as above
 472          for( k=0;k<NNs;k++)
 473           {
 474            if (tmp==NSidx[k])
 475             {
 476              if (debug) cout<<"Such a digit exists \n";
 477              bit=0;
 478             }
 479           }
 480          // end comment
 481          if (bit) {
 482           fDigitsIndexN->AddAt(i,fNDigitsN++);
 483           NSidx[NNs++] =tmp;
 484          }
 485       }
 486  }
 487
 488  delete [] PSidx;
 489  delete [] NSidx;
 490
 491  if (debug) cout<<"Digits :  P = "<<fNDigitsP<<"   N = "<<fNDigitsN<<endl;
 492
 493 }
 494
 495
 496 //-------------------------------------------
 497
 498 void AliITSClusterFinderSSD::SortDigits()
 499 {
 500 // sort digits
 501
 502   Int_t i;
 503   if(fNDigitsP>1)
 504   for (i=0;i<fNDigitsP-1;i++)
 505   if (SortDigitsP(0,(fNDigitsP-1-i))==0) break;
 506
 507   if(fNDigitsN>1)
 508     for (i=0;i<fNDigitsN-1;i++)
 509   if(SortDigitsN(0,(fNDigitsN-1-i))==0) break;
 510 }
 511
 512
 513
 514 //----------------------------------------------
 515 void AliITSClusterFinderSSD::FillClIndexArrays(Int_t* arrayP, Int_t *arrayN)
 516 {
 517 // fill cluster index array
 518
 519   register Int_t i;
 520   for (i=0; i<fNClusterP;i++)
 521     {
 522       arrayP[i]=i;
 523     }
 524   for (i=0; i<fNClusterN;i++)
 525     {
 526       arrayN[i]=i;
 527     }
 528 }
 529
 530
 531 //------------------------------------------------------
 532 void AliITSClusterFinderSSD::SortClusters(Int_t* arrayP, Int_t *arrayN)
 533 {
 534 // sort clusters
 535
 536   Int_t i;
 537   if(fNClusterP>1)
 538     for (i=0;i<fNClusterP-1;i++)
 539       if (SortClustersP(0,(fNClusterP-1),arrayP)==0)  break;
 540
 541
 542   if(fNClusterN>1)
 543     for (i=0;i<fNClusterN-1;i++)
 544       if (SortClustersN(0,(fNClusterN-1),arrayN)==0)  break;
 545
 546 }
 547
 548
 549
 550 //---------------------------------------------------
 551 Int_t AliITSClusterFinderSSD::SortClustersP(Int_t start, Int_t end, Int_t *array)
 552 {
 553 //Sort P side clusters
 554
 555
 556   Int_t right;
 557   Int_t left;
 558   if (start != (end - 1) ) {
 559       left=this->SortClustersP(start,(start+end)/2,array);
 560       right=this->SortClustersP((start+end)/2,end,array);
 561       return (left || right);
 562   } else {
 563       left =((AliITSclusterSSD*)((*fClusterP)[array[start]]))->
 564                                                          GetDigitStripNo(0);
 565       right=((AliITSclusterSSD*)((*fClusterP)[array[ end ]]))->
 566                                                          GetDigitStripNo(0);
 567       if(left>right) {
 568          Int_t tmp = array[start];
 569          array[start]=array[end];
 570          array[end]=tmp;
 571          return 1;
 572       } else return 0;
 573   }
 574
 575
 576 }
 577
 578
 579
 580 //-------------------------------------------------------
 581 Int_t AliITSClusterFinderSSD::SortClustersN(Int_t start, Int_t end, Int_t *array)
 582 {
 583 //Sort N side clusters
 584
 585
 586   Int_t right;
 587   Int_t left;
 588
 589   if (start != (end - 1) ) {
 590       left=this->SortClustersN(start,(start+end)/2,array);
 591       right=this->SortClustersN((start+end)/2,end,array);
 592       return (left || right);
 593   } else {
 594       left =((AliITSclusterSSD*)((*fClusterN)[array[start]]))->
 595                                                          GetDigitStripNo(0);
 596       right=((AliITSclusterSSD*)((*fClusterN)[array[ end ]]))->
 597                                                          GetDigitStripNo(0);
 598       if( left > right) {
 599          Int_t tmp = array[start];
 600          array[start]=array[end];
 601          array[end]=tmp;
 602          return 1;
 603       } else return 0;
 604     }
 605
 606 }
 607
 608
 609
 610 //-------------------------------------------------------
 611 void AliITSClusterFinderSSD::ClustersToPackages()
 612 {
 613 // fill packages
 614
 615   Int_t *oneSclP = new Int_t[fNClusterP]; //I want to have sorted 1 S clusters
 616   Int_t *oneSclN = new Int_t[fNClusterN]; //I can not sort it in TClonesArray
 617                                           //so, I create table of indexes and
 618                                           //sort it
 619                                           //I do not use TArrayI on purpose
 620                                           // MB: well, that's not true that one
 621                                           //cannot sort objects in TClonesArray
 622
 623   register Int_t i;    //iterator
 624   AliITSpackageSSD *currentpkg;
 625   AliITSclusterSSD *currentP;
 626   AliITSclusterSSD *currentN;
 627   AliITSclusterSSD *fcurrN;
 628
 629   Int_t     lastNclIndex=0;          //index of last N sidecluster
 630   Int_t     Nclpack=0;               //number of P clusters in current package
 631   Int_t     Pclpack=0;               //number of N clusters in current package
 632   Int_t     firstStripP;
 633   Int_t     lastStripP;
 634   Int_t     firstStripN;
 635   Int_t     lastStripN;
 636   Int_t     tmplastNclIndex;
 637
 638 //Fills in One Side Clusters Index Array
 639   FillClIndexArrays(oneSclP,oneSclN);
 640 //Sorts filled Arrays
 641   SortClusters(oneSclP,oneSclN);
 642
 643
 644   fNPackages=1;
 645   new ((*fPackages)[0]) AliITSpackageSSD(fClusterP,fClusterN);
 646   currentpkg = (AliITSpackageSSD*)((*fPackages)[0]);
 647
 648   //main loop on sorted P side clusters
 649   Int_t k;
 650   for (i=0;i<fNClusterP;i++) {
 651       //Take new P side cluster
 652       currentP = GetPSideCluster(oneSclP[i]);
 653       currentN = GetNSideCluster(oneSclN[lastNclIndex]);
 654       // take a new P cluster or not ?
 655       if(IsCrossing(currentP,currentN)) {
 656           // don't take a new P cluster
 657           Pclpack++;
 658           currentP->AddCross(oneSclN[lastNclIndex]);
 659           currentN->AddCross(oneSclP[i]);
 660           currentpkg->AddPSideCluster(oneSclP[i]);
 661           if (Nclpack==0) {
 662               currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
 663               Nclpack++;
 664           }
 665           //check if previous N side clusters crosses with it too
 666           for(k=1;k<fSFB+1;k++) {
 667               //check if we are not out of array
 668               if ((lastNclIndex-k)>-1) {
 669                   fcurrN = GetNSideCluster( oneSclN[lastNclIndex-k] );
 670                   if( IsCrossing(currentP,fcurrN) ) {
 671                       currentP->AddCross(oneSclN[lastNclIndex-k]);
 672                       fcurrN->AddCross(oneSclP[i]);
 673                  } else break; //There is no sense to check rest of clusters
 674               } else break;
 675           }
 676           tmplastNclIndex=lastNclIndex;
 677           //Check if next N side clusters crosses with it too
 678           for (Int_t k=1;k<fSFF+1;k++) {
 679               //check if we are not out of array
 680               if ((tmplastNclIndex+k)<fNClusterN) {
 681                   fcurrN = GetNSideCluster( oneSclN[tmplastNclIndex+k] );
 682                   if(IsCrossing(currentP,fcurrN) ) {
 683                      lastNclIndex++;
 684                      fcurrN->AddCross(oneSclP[i]);
 685                      currentP->AddCross(oneSclN[tmplastNclIndex+k]);
 686                      currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
 687                      Nclpack++;
 688                   }else break;
 689               } else break;
 690           }
 691
 692           // end of package
 693           //if( IsCrossing )
 694       } else {
 695          lastStripP  = currentP->GetLastDigitStripNo();
 696          lastStripN  = currentN->GetLastDigitStripNo();
 697
 698          if((lastStripN-fSFF) < lastStripP ) {
 699             //Take new PCluster
 700             if((Nclpack>0)&& (Pclpack>0)) {
 701               new ((*fPackages)[fNPackages]) AliITSpackageSSD(fClusterP,fClusterN);
 702               currentpkg = (AliITSpackageSSD*)((*fPackages)[fNPackages++]);
 703             }
 704
 705             //so we have to take another cluster on side N and check it
 706             //we have to check until taken cluster's last strip will
 707             //be > last of this cluster(P)
 708             //so it means that there is no corresponding cluster on side N
 709             //to this cluster (P)
 710             //so we have to continue main loop with new "lastNclIndex"
 711
 712             Nclpack=0;
 713             Pclpack=0;
 714             //We are not sure that next N cluter will cross with this P cluster
 715             //There might one, or more, clusters that do not have any
 716             //corresponding cluster on the other side
 717             for(;;) {
 718                lastNclIndex++;
 719                //Check if we are not out of array
 720                if (lastNclIndex<fNClusterN) {
 721                   currentN = GetNSideCluster(oneSclN[lastNclIndex]);
 722                   if ( IsCrossing(currentP, currentN) ){
 723                      Nclpack++;
 724                      Pclpack++;
 725                      currentP->AddCross(oneSclN[lastNclIndex]);
 726                      currentN->AddCross(oneSclP[i]);
 727                      currentpkg->AddNSideCluster(oneSclN[lastNclIndex]);
 728                      currentpkg->AddPSideCluster(oneSclP[i]);
 729                      //Check, if next N side clusters crosses with it too
 730                      tmplastNclIndex=lastNclIndex;
 731                      for (Int_t k=1;k<fSFF+1;k++) {
 732                        //Check if we are not out of array
 733                        if ((tmplastNclIndex+k)<fNClusterN) {
 734                           fcurrN = GetNSideCluster(oneSclN[tmplastNclIndex+k]);
 735                           if( IsCrossing(currentP, fcurrN) ) {
 736                              Nclpack++;
 737                              lastNclIndex++;
 738                              currentP->AddCross(oneSclN[tmplastNclIndex+k]);
 739                              fcurrN->AddCross(oneSclP[i]);
 740                              currentpkg->AddNSideCluster(oneSclN[tmplastNclIndex+k]);
 741                           }else break;
 742                        } else break; //we are out of array
 743                      } // end loop
 744                      break;
 745                   } else {
 746                      firstStripP = currentP->GetFirstDigitStripNo();
 747                      firstStripN = currentN->GetFirstDigitStripNo();
 748                      if((firstStripN+fSFB)>=firstStripP) break;
 749                      else continue;
 750                   }
 751                } else goto EndOfFunction;
 752             } // end for(;;)
 753          }  else  //EndOfPackage
 754            {
 755              continue;
 756            } //else EndOfPackage
 757
 758       }//else IsCrossing
 759
 760   }//main for
 761
 762  EndOfFunction:
 763   if ((Nclpack<1) || (Pclpack<1)) fNPackages--;
 764
 765    delete oneSclP;
 766    delete oneSclN;
 767
 768 }
 769
 770
 771
 772 //-----------------------------------------------
 773 void AliITSClusterFinderSSD::PackagesToPoints()
 774 {
 775 // find recpoints from packages
 776
 777  register Int_t i;
 778  AliITSpackageSSD *currentpkg;
 779  Int_t NumNcl;
 780  Int_t NumPcl;
 781  Int_t clusterIndex;
 782  Bool_t clusterSide;
 783
 784  for (i=0;i<fNPackages;i++) {
 785     //get pointer to next package
 786     currentpkg = (AliITSpackageSSD*)((*fPackages)[i]);
 787     NumNcl = currentpkg->GetNumOfClustersN();
 788     NumPcl = currentpkg->GetNumOfClustersP();
 789     if(debug) cout<<"New Package\nNumPcl ="<<NumPcl<<" NumNcl ="<<NumNcl<<"\n";
 790
 791     while(((NumPcl>=2)&&(NumNcl>2))||((NumPcl>2)&&(NumNcl>=2))) {
 792         //This is case of "big" pacakage
 793         //conditions see 3 lines above
 794         //if in the big package exists cluster with one cross
 795         //we can reconstruct this point without any geometrical ambiguities
 796         if(currentpkg->GetClusterWithOneCross(clusterIndex, clusterSide) )
 797         {
 798             ResolveClusterWithOneCross(currentpkg, clusterIndex, clusterSide);
 799         } else if (clusterIndex == -2) {
 800             NumPcl = 0;
 801             NumNcl = 0;
 802             break;
 803         } else {
 804             if ( (NumNcl==NumPcl) && (NumPcl<10)) {
 805                //if there is no cluster with one cross
 806                //we can resolve whole package at once
 807                //by finding best combination
 808                //we can make combinations when NumNcl==NumPcl
 809                if (ResolvePackageBestCombin(currentpkg)) {
 810                    //sometimes creating combinations fail,
 811                    //but it happens very rarely
 812                    NumPcl = 0;
 813                    NumNcl = 0;
 814                    break;
 815                } else ResolveOneBestMatchingPoint(currentpkg);
 816             } else {
 817                ResolveOneBestMatchingPoint(currentpkg);
 818             }
 819         }
 820         NumNcl = currentpkg->GetNumOfClustersN();
 821         NumPcl = currentpkg->GetNumOfClustersP();
 822     } // end while
 823     if ((NumPcl==1)&&(NumNcl==1)) {
 824       ResolveSimplePackage(currentpkg);
 825        continue;
 826     }
 827     if (NumPcl==1) {
 828       ResolvePackageWithOnePSideCluster(currentpkg);
 829       continue;
 830     }
 831     if (NumNcl==1) {
 832       ResolvePackageWithOneNSideCluster(currentpkg);
 833       continue;
 834     }
 835     if ((NumPcl==2)&&(NumNcl==2)) {
 836       ResolveTwoForTwoPackage(currentpkg);
 837       continue;
 838     }
 839     if ((NumPcl==0)&&(NumNcl>0)) { }
 840     if ((NumNcl==0)&&(NumPcl>0)) { }
 841
 842   } // end loop over fNPackages
 843
 844
 845 }
 846
 847
 848 //----------------------------------------------------------
 849
 850 void AliITSClusterFinderSSD::
 851 ResolveClusterWithOneCross(AliITSpackageSSD *currentpkg, Int_t clusterIndex, Bool_t clSide)
 852 {
 853 // resolve clusters with one cross
 854 /*
 855
 856 ie:
 857     \    /
 858      \  /
 859       \/
 860       /\
 861      /  \
 862     /    \
 863    /      \
 864   /        \
 865  /          \
 866 */
 867
 868
 869
 870    if (clSide == fgkSIDEP) ResolvePClusterWithOneCross(currentpkg,clusterIndex);
 871    else  ResolveNClusterWithOneCross(currentpkg,clusterIndex);
 872
 873 }
 874
 875
 876 //---------------------------------------------------------
 877 void AliITSClusterFinderSSD::
 878 ResolvePClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
 879 {
 880 // resolve clusters with one cross on P side
 881 /*
 882 There is cluster (side P) which crosses with only one cluster on the other
 883 side (N)
 884
 885 ie:
 886     \    /  \/   /
 887      \  /   /\  /
 888       \/   /  \/
 889       /\  /   /\
 890      /  \/   /  \    .....
 891     /   /\  /    \
 892    /   /  \/      \
 893   /   /   /\       \
 894  /   /   /  \       \
 895 */
 896
 897
 898   AliITSclusterSSD * clusterP;
 899   AliITSclusterSSD * clusterN;
 900
 901   Float_t posClusterP;           //Cluster P position in strip coordinates
 902   Float_t posClusterN;           //Cluster N position in strip coordinates
 903
 904   Float_t posErrorClusterP;
 905   Float_t posErrorClusterN;
 906
 907   Float_t sigClusterP;
 908   Float_t sigClusterN;
 909
 910   Float_t sigErrorClusterP;
 911   Float_t sigErrorClusterN;
 912
 913   Float_t ps;
 914   Float_t ns;
 915
 916   Float_t Chicomb;
 917   Int_t clusterIdx;
 918
 919   if (debug) cout<<"ResolvePClusterWithOneCross\n";
 920
 921   clusterP=pkg->GetPSideCluster(clusterIndex);
 922   posClusterP = GetClusterZ(clusterP);
 923   posErrorClusterP = clusterP->GetPositionError();
 924   sigClusterP = ps = clusterP->GetTotalSignal();
 925   sigErrorClusterP = clusterP->GetTotalSignalError();
 926   //carefully, it returns index in TClonesArray
 927
 928   clusterIdx = clusterP->GetCross(0);
 929   clusterN=this->GetNSideCluster(clusterIdx);
 930   ns = clusterN->GetTotalSignal();
 931   posClusterN = GetClusterZ(clusterN);
 932   posErrorClusterN = clusterN->GetPositionError();
 933   pkg->DelCluster(clusterIndex,fgkSIDEP);
 934   sigClusterN = ps/fPNsignalRatio;
 935   // there is no sonse to check how signal ratio is far from perfect
 936   // matching line if the if below it is true
 937   if (ns < sigClusterN) {
 938       sigClusterN=ns;
 939       if (debug) cout<<"n1 < p1/fPNsignalRatio";
 940       if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<" ... \n";
 941       pkg->DelClusterOI(clusterIdx,fgkSIDEN);
 942   } else {
 943       //Let's see how signal ratio is far from perfect matching line
 944       Chicomb = DistToPML(ps,ns);
 945       if (debug) cout<<"Chic "<<Chicomb<<"\n";
 946       if (Chicomb > fAlpha2) {
 947          //it is near, so we can risk throwing this cluster away too
 948          if (debug) cout<<"Attempting to del cluster N "<<clusterIdx<<"...\n";
 949          pkg->DelClusterOI(clusterIdx,fgkSIDEN);
 950       } else {
 951          clusterN->CutTotalSignal(sigClusterN);
 952          if (debug) cout <<"Signal cut  |||||||||||||\n";
 953       }
 954   }
 955   sigErrorClusterN= clusterN->GetTotalSignalError();
 956   CreateNewRecPoint(posClusterP,posErrorClusterP,posClusterN,posErrorClusterN,
 957                     sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
 958                     clusterP, clusterN, 0.75);
 959
 960 }
 961
 962
 963 //------------------------------------------------------
 964 void AliITSClusterFinderSSD::
 965 ResolveNClusterWithOneCross(AliITSpackageSSD *pkg, Int_t clusterIndex)
 966 {
 967 // resolve clusters with one cross on N side
 968
 969
 970   AliITSclusterSSD * clusterP;
 971   AliITSclusterSSD * clusterN;
 972
 973   Float_t posClusterP;              //Cluster P position in strip coordinates
 974   Float_t posClusterN;              //Cluster N position in strip coordinates
 975
 976   Float_t posErrorClusterP;
 977   Float_t posErrorClusterN;
 978
 979   Float_t sigClusterP;
 980   Float_t sigClusterN;
 981
 982   Float_t sigErrorClusterP;
 983   Float_t sigErrorClusterN;
 984
 985   Float_t ps;
 986   Float_t ns;
 987
 988   Float_t Chicomb;
 989   Int_t clusterIdx;
 990
 991   if (debug) cout<<"ResolveNClusterWithOneCross\n";
 992
 993   clusterN=pkg->GetNSideCluster(clusterIndex);
 994   posClusterN = GetClusterZ(clusterN);
 995   posErrorClusterN = clusterN->GetPositionError();
 996   sigClusterN = ns = clusterN->GetTotalSignal();
 997   sigErrorClusterN = clusterN->GetTotalSignalError();
 998   //carefully, it returns index in TClonesArray
 999   clusterIdx = clusterN->GetCross(0);
1000   clusterP=this->GetPSideCluster(clusterIdx);
1001   ps = clusterP->GetTotalSignal();
1002   posClusterP = GetClusterZ(clusterP);
1003   posErrorClusterP = clusterP->GetPositionError();
1004   pkg->DelCluster(clusterIndex,fgkSIDEN);
1005   sigClusterP=ns*fPNsignalRatio;
1006   // there is no sonse to check how signal ratio is far from perfect
1007   // matching line if the if below it is true
1008   if (ps < sigClusterP) {
1009       sigClusterP = ps;
1010       if (debug) cout<<"ps < ns*fPNsignalRatio";
1011       if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<" ... \n";
1012       pkg->DelClusterOI(clusterIdx,fgkSIDEP);
1013   } else {
1014       //Let's see how signal ratio is far from perfect matching line
1015        Chicomb = DistToPML(ps,ns);
1016        if (debug) cout<<"Chic "<<Chicomb<<"\n";
1017        if (Chicomb > fAlpha2) {
1018           //it is near, so we can risk frowing this cluster away too
1019           if (debug) cout<<"Attempting to del cluster P "<<clusterIdx<<"...\n";
1020           pkg->DelClusterOI(clusterIdx,fgkSIDEP);
1021        } else {
1022           clusterN->CutTotalSignal(sigClusterP);
1023           if (debug) cout <<"Signal cut  ------------\n";
1024        }
1025   }
1026   sigErrorClusterP= clusterP->GetTotalSignalError();
1027   CreateNewRecPoint( posClusterP,posErrorClusterP,posClusterN,posErrorClusterN,
1028                      sigClusterP+sigClusterN,sigErrorClusterN+sigErrorClusterP,
1029                      clusterP, clusterN, 0.75);
1030
1031
1032 }
1033
1034
1035
1036 //-------------------------------------------------------
1037
1038 Bool_t AliITSClusterFinderSSD::
1039 ResolvePackageBestCombin(AliITSpackageSSD *pkg)
1040 {
1041 //find best combination
1042
1043   if (debug) cout<<"NumNcl==NumPcl ("<<pkg->GetNumOfClustersN()
1044       <<"=="<<pkg->GetNumOfClustersP()<<"); Generating combinations ... \n";
1045
1046
1047   AliITSclusterSSD * clusterP;
1048   AliITSclusterSSD * clusterN;
1049
1050   Int_t Ncombin; //Number of combinations
1051   Int_t itera;   //iterator
1052   Int_t sizet=1;   //size of array to allocate
1053
1054   Int_t NP = pkg->GetNumOfClustersP();
1055   for (itera =2; itera <= NP ;itera ++) {
1056      sizet=sizet*itera;
1057      if (sizet > 10000) {
1058        sizet=10000;
1059        break;
1060      }
1061   }
1062
1063   Int_t** combin = new Int_t*[sizet]; //2D array to keep combinations in
1064
1065   for (itera =0; itera <sizet;itera++) {
1066    combin[itera] = new Int_t[NP+1];
1067   }
1068
1069   pkg->GetAllCombinations(combin,Ncombin,sizet);
1070
1071   if (Ncombin==0) {
1072     if (debug) cout<<"No combin Error";
1073     return kFALSE;
1074   }
1075 //calculate which combination fits the best to perfect matching line
1076   Int_t bc = GetBestComb(combin,Ncombin,NP,pkg);
1077   if (debug) cout<<"  bc "<<bc <<"\n";
1078
1079   for (itera =0; itera < NP; itera ++) {
1080       clusterP = pkg->GetPSideCluster(itera);
1081       //carefully here
1082       //becase AliITSclusterSSD::GetCross returns index in
1083       //AliITSmoduleSSD.fClusterP, which is put to "combin"
1084       clusterN = GetNSideCluster(combin[bc][itera]);
1085       CreateNewRecPoint(clusterP, clusterN, 0.75);
1086   }
1087
1088   for (itera =0; itera <sizet;itera++) {
1089      delete [](combin[itera]);
1090   }
1091   delete [] combin;
1092   return kTRUE;
1093
1094 }
1095
1096
1097 //----------------------------------------------------
1098 void AliITSClusterFinderSSD::
1099 ResolveOneBestMatchingPoint(AliITSpackageSSD *pkg)
1100 {
1101 // find best matching point
1102
1103  Int_t ni, pi;
1104
1105  Int_t prvP, prvN;
1106  Int_t nextP, nextN=0;
1107
1108
1109  AliITSclusterSSD * clusterP;
1110  AliITSclusterSSD * clusterN;
1111
1112  Bool_t split = kFALSE;
1113
1114  if (debug) cout<<"ResolveOneBestMatchingPoint \n";
1115
1116  GetBestMatchingPoint(pi, ni, pkg);
1117
1118  clusterP = GetPSideCluster(pi);
1119  clusterN = GetNSideCluster(ni);
1120
1121  CreateNewRecPoint(clusterP, clusterN, 0.75);
1122
1123  if ((nextP=pkg->GetNextPIdx(pi))!=-1)
1124    if ((nextN=pkg->GetNextNIdx(ni))!=-1)
1125      if ((prvP=pkg->GetPrvPIdx(pi))!=-1)
1126        if ((prvN=pkg->GetPrvNIdx(ni))!=-1)
1127           if( !(GetPSideCluster(prvP)->IsCrossingWith(nextN)))
1128             if( !(GetPSideCluster(nextP)->IsCrossingWith(prvN)))
1129               {
1130                  split = kTRUE;
1131               }
1132
1133
1134  pkg->DelClusterOI(pi, fgkSIDEP);
1135  pkg->DelClusterOI(ni, fgkSIDEN);
1136
1137  if (split) {
1138    if (debug) cout<<"spltting package ...\n";
1139    new ((*fPackages)[fNPackages]) AliITSpackageSSD(fClusterP,fClusterN);
1140    pkg->
1141      SplitPackage(nextP,nextN,(AliITSpackageSSD*)((*fPackages)[fNPackages++]));
1142  }
1143
1144 }
1145
1146
1147 //--------------------------------------------------
1148 void  AliITSClusterFinderSSD::ResolveSimplePackage(AliITSpackageSSD *pkg)
1149 {
1150 // resolve simple package
1151
1152   AliITSclusterSSD * clusterP;
1153   AliITSclusterSSD * clusterN;
1154
1155   clusterP         = pkg->GetPSideCluster(0);
1156
1157   clusterN         = pkg->GetNSideCluster(0);
1158
1159   CreateNewRecPoint(clusterP, clusterN, 1.0);
1160
1161
1162 }
1163
1164
1165 //--------------------------------------------------
1166 void AliITSClusterFinderSSD:: ResolvePackageWithOnePSideCluster(AliITSpackageSSD *pkg) {
1167 // resolve P side clusters from packages
1168
1169 /*
1170  \   \   \  /
1171   \   \   \/
1172    \   \  /\
1173     \   \/  \
1174      \  /\   \
1175       \/  \   \
1176       /\   \   \
1177      /  \   \   \
1178
1179 */
1180
1181
1182 /************************/
1183 /**************************/
1184 /*XP SP itg jest tylko jeden nie musi byc tablica i w peetli nie trzeba po niej iterowcac*/
1185 /***************************/
1186
1187
1188  Int_t k;
1189  AliITSclusterSSD * clusterP;
1190  AliITSclusterSSD * clusterN;
1191
1192  Int_t NN=pkg->GetNumOfClustersN();
1193  Float_t sumsig = 0;
1194
1195  Float_t   XP;
1196  Float_t   XPerr;
1197
1198  Float_t * XN = new Float_t[NN];
1199  Float_t * XNerr = new Float_t[NN];
1200
1201  Float_t * SP = new Float_t[NN];
1202  Float_t * SN = new Float_t[NN];
1203
1204  Float_t * SPerr = new Float_t[NN];
1205  Float_t * SNerr = new Float_t[NN];
1206
1207  Float_t p1;
1208  Float_t p1err;
1209
1210  clusterP = pkg->GetPSideCluster(0);
1211  p1       = clusterP->GetTotalSignal();
1212
1213  XP    = GetClusterZ(clusterP);
1214  XPerr = clusterP->GetPositionError();
1215  p1err = clusterP->GetTotalSignalError();
1216
1217  for (k=0;k<NN;k++) {
1218     SN[k]    = pkg->GetNSideCluster(k)->GetTotalSignal();
1219     SNerr[k] = pkg->GetNSideCluster(k)->GetTotalSignalError();
1220     sumsig   += SN[k];
1221  }
1222  for (k=0;k<NN;k++) {
1223     clusterN = pkg->GetNSideCluster(k);
1224     SP[k]= p1*SN[k]/sumsig;
1225     SPerr[k] = p1err*SN[k]/sumsig;
1226     XN[k]=GetClusterZ(clusterN);
1227     XNerr[k]= clusterN->GetPositionError();
1228     CreateNewRecPoint(XP,XPerr, XN[k],XNerr[k], SP[k]+SN[k], SPerr[k]+SNerr[k], clusterP, clusterN, 1.0);
1229
1230  }
1231
1232  delete [] XN;
1233  delete [] XNerr;
1234  delete [] SP;
1235  delete [] SN;
1236  delete [] SPerr;
1237  delete [] SNerr;
1238
1239
1240 }
1241
1242
1243 //---------------------------------------------------------
1244 void AliITSClusterFinderSSD::ResolvePackageWithOneNSideCluster(AliITSpackageSSD *pkg) {
1245 // resolve N side clusters from packages
1246
1247 /*
1248     \    /   /   /
1249      \  /   /   /
1250       \/   /   /
1251       /\  /   /
1252      /  \/   /
1253     /   /\  /
1254    /   /  \/
1255   /   /   /\
1256  /   /   /  \
1257 */
1258
1259
1260  Int_t k;
1261  AliITSclusterSSD * clusterP;
1262  AliITSclusterSSD * clusterN;
1263
1264  Int_t NP=pkg->GetNumOfClustersP();
1265  Float_t sumsig = 0;
1266
1267  Float_t * XP = new Float_t[NP];
1268  Float_t * XPerr = new Float_t[NP];
1269
1270  Float_t   XN;
1271  Float_t   XNerr;
1272
1273  Float_t * SP = new Float_t[NP];
1274  Float_t * SN = new Float_t[NP];
1275
1276  Float_t * SPerr = new Float_t[NP];
1277  Float_t * SNerr = new Float_t[NP];
1278
1279  Float_t n1;
1280  Float_t n1err;
1281
1282  clusterN = pkg->GetNSideCluster(0);
1283  n1       = clusterN->GetTotalSignal();
1284
1285  XN = GetClusterZ(clusterN);
1286  XNerr = clusterN->GetPositionError();
1287
1288  n1err=clusterN->GetTotalSignalError();
1289
1290  for (k=0;k<NP;k++) {
1291     SP[k] = pkg->GetPSideCluster(k)->GetTotalSignal();
1292     sumsig += SP[k];
1293     SPerr[k] = pkg->GetPSideCluster(k)->GetTotalSignalError();
1294  }
1295
1296  for (k=0;k<NP;k++) {
1297     clusterP = pkg->GetPSideCluster(k);
1298     SN[k]= n1*SP[k]/sumsig;
1299     XP[k]=GetClusterZ(clusterP);
1300     XPerr[k]= clusterP->GetPositionError();
1301     SNerr[k] = n1err*SP[k]/sumsig;
1302     CreateNewRecPoint(XP[k],XPerr[k], XN,XNerr, SP[k]+SN[k], SPerr[k]+SNerr[k],clusterP, clusterN, 1.0);
1303   }
1304
1305  delete [] XP;
1306  delete [] XPerr;
1307  delete [] SP;
1308  delete [] SN;
1309  delete [] SPerr;
1310  delete [] SNerr;
1311
1312 }
1313
1314 /**********************************************************************/
1315 /*********      2  X   2  *********************************************/
1316 /**********************************************************************/
1317
1318 void AliITSClusterFinderSSD::
1319 ResolveTwoForTwoPackage(AliITSpackageSSD *pkg)
1320 {
1321 // resolve 2x2 packages
1322
1323   AliITSclusterSSD *clusterP1 = pkg->GetPSideCluster(0);
1324   AliITSclusterSSD *clusterP2 = pkg->GetPSideCluster(1);
1325   AliITSclusterSSD *clusterN1 = pkg->GetNSideCluster(0);
1326   AliITSclusterSSD *clusterN2 = pkg->GetNSideCluster(1);
1327
1328   AliITSclusterSSD *tmp;
1329
1330   Float_t p1sig;
1331   Float_t p2sig;
1332   Float_t n1sig;
1333   Float_t n2sig;
1334
1335   Float_t p1sigErr;
1336   Float_t p2sigErr;
1337   Float_t n1sigErr;
1338   Float_t n2sigErr;
1339
1340   Float_t ZposP1;
1341   Float_t ZposP2;
1342   Float_t ZposN1;
1343   Float_t ZposN2;
1344
1345   Float_t ZposErrP1;
1346   Float_t ZposErrP2;
1347   Float_t ZposErrN1;
1348   Float_t ZposErrN2;
1349
1350   Float_t D12;
1351
1352   Float_t Chicomb1;
1353   Float_t Chicomb2;
1354
1355   if(clusterP1->GetDigitStripNo(0) > clusterP2->GetDigitStripNo(0)) {
1356      if (debug) cout<<"P strips flip\n";
1357      tmp       = clusterP1;
1358      clusterP1 = clusterP2;
1359      clusterP2 = tmp;
1360   }
1361   if(clusterN1->GetDigitStripNo(0) > clusterN2->GetDigitStripNo(0)) {
1362      if (debug) cout<<"N strips flip\n";
1363      tmp       = clusterN1;
1364      clusterN1 = clusterN2;
1365      clusterN2 = tmp;
1366   }
1367   p1sig = clusterP1->GetTotalSignal();
1368   p2sig = clusterP2->GetTotalSignal();
1369   n1sig = clusterN1->GetTotalSignal();
1370   n2sig = clusterN2->GetTotalSignal();
1371
1372
1373   p1sigErr = clusterP1->GetTotalSignalError();
1374   n1sigErr = clusterN1->GetTotalSignalError();
1375   p2sigErr = clusterP2->GetTotalSignalError();
1376   n2sigErr = clusterN2->GetTotalSignalError();
1377
1378   ZposP1 = clusterP1->GetPosition();
1379   ZposN1 = clusterN1->GetPosition();
1380   ZposP2 = clusterP2->GetPosition();
1381   ZposN2 = clusterN2->GetPosition();
1382
1383   ZposErrP1 = clusterP1->GetPositionError();
1384   ZposErrN1 = clusterN1->GetPositionError();
1385   ZposErrP2 = clusterP2->GetPositionError();
1386   ZposErrN2 = clusterN2->GetPositionError();
1387
1388   //in this may be two types:
1389   //    1.When each cluster crosses with 2 clusters on the other side
1390   //            gives 2 ghosts and 2 real points
1391   //
1392   //    2.When two clusters (one an each side) crosses with only one on
1393   //           the other side and two crosses (one on the each side) with
1394   //           two on the other gives 2 or 3 points,
1395
1396   if (debug) cout<<"2 for 2 ambiguity ...";
1397   /***************************/
1398   /*First sort of combination*/
1399   /***************************/
1400
1401   if((clusterP1->GetCrossNo()==2) && (clusterN1->GetCrossNo()==2)) {
1402      if (debug) cout<<"All clusters has 2 crosses\n";
1403      D12 = TMath::Sqrt((Float_t)((p1sig -p2sig)*(p1sig -p2sig) + (n1sig -n2sig)*(n1sig -n2sig)));
1404
1405      if(debug) cout<<"Distance between points in P(N) plane D12 = "<<D12<<"\n";
1406
1407      /*********************************************/
1408      /*resolving ambiguities:                     */
1409      /*we count Chicomb's and we take combination */
1410      /*giving lower Chicomb as a real points      */
1411      /*Keep only better combinantion              */
1412      /*********************************************/
1413
1414      if (D12 > (fAlpha3*17.8768)) {
1415        if (debug) cout<<"decided to take only one pair \n";
1416        Chicomb1 = DistToPML(p1sig,n1sig) + DistToPML(p2sig,n2sig);
1417        Chicomb2 = DistToPML(p2sig,n1sig) + DistToPML(p1sig,n2sig);
1418        if (debug) {
1419          cout<<" 00 11 combination : "<<Chicomb1<<" 01 10 combination : "<<Chicomb2<<" \n";
1420          cout<<"p1 = "<<p1sig<<"  n1 = "<<n1sig<<"  p2 = "<<p2sig<<"  n2 = "<<n2sig<<"\n";
1421        }
1422        if (Chicomb1 < Chicomb2) {
1423           if (debug) cout<<"00 11";
1424           CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 0.75);
1425           CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 0.75);
1426
1427
1428           //second cominantion
1429        } else {
1430           if (debug) cout<<"01 10";
1431           CreateNewRecPoint(ZposP1,0, ZposN2,0, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.75);
1432           CreateNewRecPoint(ZposP2,0, ZposN1,0, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.75);
1433        } //end second combinantion
1434        //if (D12 > fAlpha3*17.8768)
1435        //keep all combinations
1436      } else {
1437         if (debug) cout<<"We decide to take all points\n";
1438         CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 0.5);
1439         CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 0.5);
1440         CreateNewRecPoint(ZposP1,ZposErrP1, ZposN2,ZposErrN2, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.5);
1441         CreateNewRecPoint(ZposP2,ZposErrP2, ZposN1,ZposErrN1, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.5);
1442      }
1443   }
1444   // ad.2 Second type of combination
1445   else {
1446     Chicomb1 = DistToPML(p1sig,n1sig) + DistToPML(p2sig,n2sig);
1447     if (debug) cout<<"\nhere can be reconstructed 3 points: chicomb = "<<Chicomb1<<"\n";
1448
1449     if (Chicomb1<fAlpha1) {
1450        if (debug) cout<<"\nWe decided to take 3rd point";
1451        if (clusterP1->GetCrossNo()==1) {
1452           if (debug) cout<<"...  P1 has one cross\n";
1453           n1sig = p1sig/fPNsignalRatio;
1454           p2sig = n2sig*fPNsignalRatio;
1455
1456           clusterN1->CutTotalSignal(n1sig);
1457           clusterP2->CutTotalSignal(p2sig);
1458
1459           CreateNewRecPoint(ZposP2,ZposErrP2, ZposN1,ZposErrN1, p2sig+n1sig, p2sigErr+n1sigErr, clusterP2, clusterN1, 0.5);
1460
1461           n1sig = clusterN1->GetTotalSignal();
1462           p2sig = clusterP2->GetTotalSignal();
1463
1464        } else {
1465           if (debug) cout<<"...  N1 has one cross\n";
1466
1467           n2sig=p2sig/fPNsignalRatio;
1468           p1sig=n1sig*fPNsignalRatio;
1469
1470           clusterN2->CutTotalSignal(n2sig);
1471           clusterP1->CutTotalSignal(p1sig);
1472
1473           CreateNewRecPoint(ZposP1,ZposErrP1, ZposN2,ZposErrN2, p1sig+n2sig, p1sigErr+n2sigErr, clusterP1, clusterN2, 0.5);
1474
1475           n2sig=clusterN2->GetTotalSignal();
1476           p1sig=clusterP1->GetTotalSignal();
1477        }
1478     } else {
1479        if (debug) cout<<"\nWe decided  NOT to take 3rd point\n";
1480     }
1481
1482     CreateNewRecPoint(ZposP1,ZposErrP1, ZposN1,ZposErrN1, p1sig+n1sig, p1sigErr+n1sigErr, clusterP1, clusterN1, 1.0);
1483     CreateNewRecPoint(ZposP2,ZposErrP2, ZposN2,ZposErrN2, p2sig+n2sig, p2sigErr+n2sigErr, clusterP2, clusterN2, 1.0);
1484
1485   }
1486
1487 }
1488
1489
1490
1491 //------------------------------------------------------
1492 Bool_t AliITSClusterFinderSSD::
1493 CreateNewRecPoint(Float_t P, Float_t dP, Float_t N, Float_t dN,
1494                   Float_t Sig,Float_t dSig,
1495                   AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN,
1496                   Stat_t prob)
1497 {
1498 // create the recpoints
1499
1500
1501   const Float_t kdEdXtoQ = 2.778e+8;
1502   const Float_t kconv = 1.0e-4;
1503   const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1504   const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1505
1506
1507   Float_t p=P;
1508   Float_t n=N;
1509   Int_t stripP, stripN;
1510   Int_t sigP, sigN;
1511   AliITSdigitSSD *digP, *digN;
1512   if (GetCrossing(P,N)) {
1513      GetCrossingError(dP,dN);
1514      AliITSRawClusterSSD cnew;
1515      Int_t nstripsP=clusterP->GetNumOfDigits();
1516      Int_t nstripsN=clusterN->GetNumOfDigits();
1517      cnew.fMultiplicity=nstripsP;
1518      cnew.fMultiplicityN=nstripsN;
1519      /*
1520      if (nstripsP > 100) {
1521         printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1522         nstripsP=100;
1523      }
1524      for(int i=0;i<nstripsP;i++) {
1525        // check if 'clusterP->GetDigitStripNo(i)' returns the digit index
1526        cnew.fIndexMap[i] = clusterP->GetDigitStripNo(i);
1527      }
1528      if (nstripsN > 100) {
1529         printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1530         nstripsN=100;
1531      }
1532      for(int i=0;i<nstripsN;i++) {
1533        // check if 'clusterN->GetDigitStripNo(i)' returns the digit index
1534        cnew.fIndexMapN[i] = clusterN->GetDigitStripNo(i);
1535      }
1536      */
1537      cnew.fQErr=dSig;
1538      //cnew.fProbability=(float)prob;
1539      fITS->AddCluster(2,&cnew);
1540      fSegmentation->GetPadIxz(P,N,stripP,stripN);
1541      digP = (AliITSdigitSSD*)fMap->GetHit(1,stripP);
1542      digN = (AliITSdigitSSD*)fMap->GetHit(0,stripN);
1543      //printf("SSD: digP digN %p %p\n",digP,digN);
1544      if(digP) sigP = digP->fSignal;
1545      else sigP=0;
1546      if(digN) sigN = digN->fSignal;
1547      else sigN=0;
1548      if (!digP && !digN) {
1549        Error("CreateNewRecPoint","cannot find the digit!");
1550        return kFALSE;
1551      }
1552      // add the rec point info
1553      AliITSRecPoint rnew;
1554      rnew.SetX(P*kconv);
1555      rnew.SetZ(N*kconv);
1556      rnew.SetQ(Sig);
1557      rnew.SetdEdX(Sig/kdEdXtoQ);
1558      rnew.SetSigmaX2( kRMSx* kRMSx);
1559      rnew.SetSigmaZ2( kRMSz* kRMSz);
1560      //rnew.SetProbability((float)prob);
1561      if(sigP > sigN) {
1562         rnew.fTracks[0]=digP->fTracks[0];
1563         rnew.fTracks[1]=digP->fTracks[1];
1564         rnew.fTracks[2]=digP->fTracks[2];
1565         //printf("sigP>sigN: %d %d %d\n",digP->fTracks[0],digP->fTracks[1],digP->fTracks[2]);
1566      } else {
1567         rnew.fTracks[0]=digN->fTracks[0];
1568         rnew.fTracks[1]=digN->fTracks[1];
1569         rnew.fTracks[2]=digN->fTracks[2];
1570         //printf("sigN>sigP: %d %d %d\n",digN->fTracks[0],digN->fTracks[1],digN->fTracks[2]);
1571      }
1572      //printf("SSD: track1 track2 track3 X Z %d %d %d %f %f\n",rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],P*kconv,N*kconv);
1573      fITS->AddRecPoint(rnew);
1574      /*
1575      // it was
1576      fPointsM->AddLast( (TObject*)
1577                    new AliITSRecPointSSD(P,dP,0,0,N,dN,Sig,dSig,fLayer,fLad,fDet,clusterP,clusterN,prob) );
1578      */
1579
1580      if(debug) cout<<"\n"<<": ImpactPoint Created: X = "<<P<<" Z = "<<N<<";   P = "<<p<<"  N = "<<n<<"\n";
1581      return kTRUE;
1582   } else {
1583      if (debug) {
1584        cout<<"\n"<<": ATENTION : stupid ImpactPoit Point X = "<<P<<" Z = "<<N<<";  P = "<<p<<"  N = "<<n<<"\n";
1585      }
1586      return kFALSE;
1587   }
1588 }
1589
1590
1591
1592 /**********************************************************************/
1593 Bool_t  AliITSClusterFinderSSD::
1594 CreateNewRecPoint(AliITSclusterSSD *clusterP, AliITSclusterSSD *clusterN, Stat_t prob)
1595 {
1596 // create recpoints
1597
1598   Float_t posClusterP;  //Cluster P position in strip coordinates
1599   Float_t posClusterN;  //Cluster N position in strip coordinates
1600
1601   Float_t posErrorClusterP;
1602   Float_t posErrorClusterN;
1603
1604   Float_t sigClusterP;
1605   Float_t sigClusterN;
1606
1607   Float_t sigErrorClusterP;
1608   Float_t sigErrorClusterN;
1609
1610   posClusterP      = clusterP->GetPosition();
1611   posErrorClusterP = clusterP->GetPositionError();
1612   sigClusterP      = clusterP->GetTotalSignal();
1613   sigErrorClusterP = clusterP->GetTotalSignalError();
1614
1615   posClusterN      = clusterN->GetPosition();
1616   posErrorClusterN = clusterN->GetPositionError();
1617   sigClusterN      = clusterN->GetTotalSignal();
1618   sigErrorClusterN = clusterN->GetTotalSignalError();
1619
1620   return CreateNewRecPoint( posClusterP,posErrorClusterP,  posClusterN,posErrorClusterN,
1621                      sigClusterP+sigClusterN, sigErrorClusterN+sigErrorClusterP,
1622                      clusterP, clusterN, prob);
1623
1624 }
1625
1626 //--------------------------------------------------
1627 Bool_t AliITSClusterFinderSSD::IsCrossing(AliITSclusterSSD* p, AliITSclusterSSD* n)
1628 {
1629 // check for crossings
1630
1631   Float_t x = p->GetPosition();
1632   Float_t y = n->GetPosition();
1633   return GetCrossing(x,y);
1634 }
1635
1636
1637 //----------------------------------------------
1638 Bool_t AliITSClusterFinderSSD::Strip2Local( Float_t stripP, Float_t stripN, Float_t &Z,Float_t &X)
1639 {
1640 // convert between strip numbers and local coordinates
1641
1642 /*
1643  Z = (stripN-stripP)*fFactorOne;
1644  X = (stripN + stripP - fStripZeroOffset)*fFactorTwo;
1645 */
1646
1647  Float_t P =stripP;
1648  Float_t N =stripN;
1649
1650  GetCrossing(P,N);
1651  X=P;
1652  Z=N;
1653  if (debug) cout<<"P="<<stripP<<" N="<<stripN<<"   X = "<<X<<" Z = "<<Z<<"\n";
1654  if ((Z<2.1)&&(Z>-2.1)&&(X<3.65)&&(X>-3.65)) return kTRUE;
1655     else return kFALSE;
1656
1657 }
1658
1659
1660 //-----------------------------------------------------------
1661 Float_t  AliITSClusterFinderSSD::GetClusterZ(AliITSclusterSSD* clust)
1662 {
1663 // get Z coordinate of the cluster
1664
1665   return clust->GetPosition();
1666
1667 }
1668
1669 //-------------------------------------------------------------
1670 Int_t AliITSClusterFinderSSD::GetBestComb
1671 (Int_t** comb,Int_t Ncomb, Int_t Ncl, AliITSpackageSSD * pkg)
1672 {
1673 //returns index of best combination in "comb"
1674 //comb : sets of combinations
1675 //       in given combination on place "n" is index of
1676 //       Nside cluster coresponding to "n" P side cluster
1677 //
1678 //Ncomb : number of combinations == number of rows in "comb"
1679 //
1680 //Ncl   : number of clusters in each row  == number of columns in "comb"
1681 //
1682 //pkg   : package
1683
1684
1685   Float_t currbval=-1;  //current best value,
1686                         //starting value is set to number negative,
1687                         //to be changed in first loop turn automatically
1688
1689   Int_t  curridx=0;     //index of combination giving currently the best chi^2
1690   Float_t chi;
1691   Float_t ps, ns;       //signal of P cluster and N cluster
1692
1693   for (Int_t i=0;i<Ncomb;i++)
1694     {
1695       chi=0;
1696
1697       for (Int_t j=0;j<Ncl;j++)
1698       {
1699         ps = pkg->GetPSideCluster(j)->GetTotalSignal();  //carrefully here, different functions
1700         ns = GetNSideCluster(comb[i][j])->GetTotalSignal();
1701         chi+=DistToPML(ps,ns);
1702       }
1703
1704       if (currbval==-1) currbval = chi;
1705       if (chi<currbval)
1706         {
1707           curridx = i;
1708           currbval  =chi;
1709         }
1710     }
1711
1712
1713   return curridx;
1714
1715 }
1716
1717 //--------------------------------------------------
1718 void AliITSClusterFinderSSD::GetBestMatchingPoint
1719 (Int_t & ip, Int_t & in, AliITSpackageSSD* pkg )
1720 {
1721 // find the best matching point
1722
1723   if (debug) pkg->PrintClusters();
1724
1725   Float_t ps, ns;  //signals on side p & n
1726
1727   Int_t nc;       // number of crosses in given p side cluster
1728   Int_t n,p;
1729   AliITSclusterSSD *curPcl, *curNcl;  //current p side cluster
1730   Float_t  bestchi, chi;
1731   ip=p=pkg->GetPSideClusterIdx(0);
1732   in=n=pkg->GetNSideClusterIdx(0);
1733
1734   bestchi=DistToPML(  pkg->GetPSideCluster(0)->GetTotalSignal(),
1735                   pkg->GetNSideCluster(0)->GetTotalSignal()  );
1736   for (Int_t i = 0; i< pkg->GetNumOfClustersP(); i++)
1737     {
1738
1739       p =  pkg->GetPSideClusterIdx(i);
1740       curPcl= GetPSideCluster(p);
1741       nc=curPcl->GetCrossNo();
1742       ps=curPcl->GetTotalSignal();
1743
1744       for (Int_t j = 0; j< nc; j++)
1745         {
1746          n=curPcl->GetCross(j);
1747          curNcl= GetNSideCluster(n);
1748          ns=curNcl->GetTotalSignal();
1749          chi = DistToPML(ps, ns);
1750
1751          if (chi <= bestchi)
1752           {
1753             bestchi = chi;
1754             in = n;
1755             ip = p;
1756           }
1757         }
1758     }
1759
1760 }
1761
1762
1763 //------------------------------------------------------
1764 void  AliITSClusterFinderSSD::CalcStepFactor(Float_t Psteo, Float_t Nsteo )
1765 {
1766 // calculate the step factor for matching clusters
1767
1768
1769   // 95 is the pitch, 4000 - dimension along z ?
1770   /*
1771   fSFF = ( (Int_t)  (Psteo*40000/95 ) );// +1;
1772   fSFB = ( (Int_t)  (Nsteo*40000/95 ) );// +1;
1773   */
1774
1775
1776   Float_t dz=fSegmentation->Dz();
1777
1778   fSFF = ( (Int_t)  (Psteo*dz/fPitch ) );// +1;
1779   fSFB = ( (Int_t)  (Nsteo*dz/fPitch ) );// +1;
1780
1781 }
1782
1783
1784 //-----------------------------------------------------------
1785 AliITSclusterSSD* AliITSClusterFinderSSD::GetPSideCluster(Int_t idx)
1786 {
1787 // get P side clusters
1788
1789
1790
1791
1792   if((idx<0)||(idx>=fNClusterP))
1793     {
1794       printf("AliITSClusterFinderSSD::GetPSideCluster  : index out of range\n");
1795       return 0;
1796     }
1797   else
1798     {
1799       return (AliITSclusterSSD*)((*fClusterP)[idx]);
1800     }
1801 }
1802
1803 //-------------------------------------------------------
1804 AliITSclusterSSD* AliITSClusterFinderSSD::GetNSideCluster(Int_t idx)
1805 {
1806 // get N side clusters
1807
1808   if((idx<0)||(idx>=fNClusterN))
1809     {
1810       printf("AliITSClusterFinderSSD::GetNSideCluster  : index out of range\n");
1811       return 0;
1812     }
1813   else
1814     {
1815       return (AliITSclusterSSD*)((*fClusterN)[idx]);
1816     }
1817 }
1818
1819 //--------------------------------------------------------
1820 AliITSclusterSSD* AliITSClusterFinderSSD::GetCluster(Int_t idx, Bool_t side)
1821 {
1822 // Get cluster
1823
1824   return (side) ? GetPSideCluster(idx) : GetNSideCluster(idx);
1825 }
1826
1827
1828 //--------------------------------------------------------
1829 void AliITSClusterFinderSSD::ConsumeClusters()
1830 {
1831 // remove clusters from the list of unmatched clusters
1832
1833  for ( Int_t i=0;i<fNPackages;i++)
1834   {
1835     ((AliITSpackageSSD*)((*fPackages)[i]))->ConsumeClusters();
1836   }
1837
1838 }
1839
1840
1841 //--------------------------------------------------------
1842 void AliITSClusterFinderSSD::ReconstructNotConsumedClusters()
1843 {
1844 // reconstruct remaining non-crossing clusters
1845
1846   //printf("ReconstructNotConsumedClusters !!!\n");
1847
1848   Int_t i;
1849   AliITSclusterSSD *cluster;
1850   Float_t pos;
1851   Float_t sig;
1852   Float_t sigerr;
1853   Float_t x1,x2,z1,z2;
1854
1855   Float_t Dx = fSegmentation->Dx();
1856   Float_t Dz = fSegmentation->Dz();
1857
1858   const Float_t kdEdXtoQ = 2.778e+8; // GeV -> number of e-hole pairs
1859   const Float_t kconv = 1.0e-4;
1860   const Float_t kRMSx = 20.0*kconv; // microns->cm ITS TDR Table 1.3
1861   const Float_t kRMSz = 830.0*kconv; // microns->cm ITS TDR Table 1.3
1862
1863   Int_t stripP,stripN;
1864   AliITSdigitSSD *dig;
1865   for (i=0;i<fNClusterP;i++)
1866     {
1867       //printf("P side cluster: cluster number %d\n",i);
1868       cluster = GetPSideCluster(i);
1869       if (!cluster->IsConsumed())
1870         {
1871           if( (pos = cluster->GetPosition()) < fSFB)
1872             {
1873               sig = cluster->GetTotalSignal();
1874
1875               sig += sig/fPNsignalRatio;
1876
1877               sigerr = cluster->GetTotalSignalError();
1878               x1 = -Dx/2 + pos *fPitch;
1879               z1 = -Dz/2;
1880
1881               x2 = pos;
1882               z2 = 1;
1883               GetCrossing (x2,z2);
1884               fSegmentation->GetPadIxz((x1+x2)/2,(z1+z2)/2,stripP,stripN);
1885               dig = (AliITSdigitSSD*)fMap->GetHit(1,stripP);
1886               //printf("dig %p\n",dig);
1887               if (!dig) printf("SSD: check the digit!  dig %p\n",dig);
1888               if(!dig) continue;
1889               AliITSRawClusterSSD cnew;
1890               Int_t nstripsP=cluster->GetNumOfDigits();
1891               cnew.fMultiplicity=nstripsP;
1892               cnew.fMultiplicityN=0;
1893               /*
1894               if (nstripsP > 100) {
1895                 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsP);
1896                 nstripsP=100;
1897               }
1898               for(int k=0;k<nstripsP;k++) {
1899                 // check if 'clusterP->GetDigitStripNo(i)' returns
1900                 // the digit index and not smth else
1901                 cnew.fIndexMap[k] = cluster->GetDigitStripNo(k);
1902               }
1903               */
1904               cnew.fQErr=sigerr;
1905               //cnew.fProbability=0.75;
1906               //printf("AddCluster: %p nstripsP %d\n",&cnew,nstripsP);
1907               fITS->AddCluster(2,&cnew);
1908               // add the rec point info
1909               AliITSRecPoint rnew;
1910               rnew.SetX(kconv*(x1+x2)/2);
1911               rnew.SetZ(kconv*(z1+z2)/2);
1912               rnew.SetQ(sig);
1913               rnew.SetdEdX(sig/kdEdXtoQ);
1914               rnew.SetSigmaX2( kRMSx* kRMSx);
1915               rnew.SetSigmaZ2( kRMSz* kRMSz);
1916               //rnew.SetProbability(0.75);
1917               rnew.fTracks[0]=dig->fTracks[0];
1918               rnew.fTracks[1]=dig->fTracks[1];
1919               rnew.fTracks[2]=dig->fTracks[2];
1920               //printf("digP: %d %d %d\n",dig->fTracks[0],dig->fTracks[1],dig->fTracks[2]);
1921               //printf("SSD - P : track1 track2 track3  X Z %d %d %d %f %f\n",rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],kconv*(x1+x2)/2,kconv*(z1+z2)/2);
1922               fITS->AddRecPoint(rnew);
1923               /*
1924               fPointsM->AddLast( (TObject*)
1925                    new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1926               */
1927
1928               if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1929                   <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<";   Pos = "<<pos;
1930             }
1931         }
1932     }
1933
1934   for (i=0;i<fNClusterN;i++)
1935     {
1936       //printf("N side cluster: cluster number %d\n",i);
1937       cluster = GetNSideCluster(i);
1938       if (!cluster->IsConsumed())
1939         {
1940           if( (pos = cluster->GetPosition()) < fSFF)
1941             {
1942               sig = cluster->GetTotalSignal();
1943
1944               sig += sig*fPNsignalRatio;
1945
1946               sigerr = cluster->GetTotalSignalError();
1947
1948               x1 = -Dx/2 + pos *fPitch;
1949               z1 = Dz/2;
1950
1951               x2 = 1;
1952               z2 = pos;
1953
1954               GetCrossing (x2,z2);
1955               fSegmentation->GetPadIxz((x1+x2)/2,(z1+z2)/2,stripP,stripN);
1956               dig = (AliITSdigitSSD*)fMap->GetHit(0,stripN);
1957               if (!dig) printf("SSD: check the digit! dig %p\n",dig);
1958               if(!dig) continue;
1959               AliITSRawClusterSSD cnew;
1960               Int_t nstripsN=cluster->GetNumOfDigits();
1961               cnew.fMultiplicity=0;
1962               cnew.fMultiplicityN=nstripsN;
1963               /*
1964               if (nstripsN > 100) {
1965                 printf("multiplicity > 100 - increase the dimension of the arrays %d\n",nstripsN);
1966                 nstripsN=100;
1967               }
1968               for(int k=0;k<nstripsN;k++) {
1969                 // check if 'clusterP->GetDigitStripNo(i)' returns
1970                 // the digit index and not smth else
1971                 cnew.fIndexMapN[k] = cluster->GetDigitStripNo(k);
1972               }
1973               */
1974               cnew.fQErr=sigerr;
1975               //cnew.fProbability=0.75;
1976               fITS->AddCluster(2,&cnew);
1977               AliITSRecPoint rnew;
1978               rnew.SetX(kconv*(x1+x2)/2);
1979               rnew.SetZ(kconv*(z1+z2)/2);
1980               rnew.SetQ(sig);
1981               rnew.SetdEdX(sig/kdEdXtoQ);
1982               rnew.SetSigmaX2( kRMSx* kRMSx);
1983               rnew.SetSigmaZ2( kRMSz* kRMSz);
1984               //rnew.SetProbability(0.75);
1985               rnew.fTracks[0]=dig->fTracks[0];
1986               rnew.fTracks[1]=dig->fTracks[1];
1987               rnew.fTracks[2]=dig->fTracks[2];
1988               //printf("digN: %d %d %d\n",dig->fTracks[0],dig->fTracks[1],dig->fTracks[2]);
1989               //printf("SSD -N : track1 track2 track3  X Z %d %d %d %f %f\n",rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],kconv*(x1+x2)/2,kconv*(z1+z2)/2);
1990               fITS->AddRecPoint(rnew);
1991               /*
1992               fPointsM->AddLast( (TObject*)
1993                    new AliITSRecPointSSD((x1+x2)/2 ,0,0,0,(z1+z2)/2,0,sig,sigerr,fLayer,fLad,fDet,cluster, 0.75) );
1994               */
1995
1996               if(debug) cout<<"\n"<<": SINGLE SIDE ImpactPoint Created: X = "
1997                   <<(x1+x2)/2<<" Z = "<<(z1+z2)/2<<";   Pos = "<<pos;
1998
1999             }
2000         }
2001     }
2002
2003
2004 }
2005
2006 //_______________________________________________________________________
2007
2008 Bool_t AliITSClusterFinderSSD::GetCrossing (Float_t &P, Float_t &N)
2009 {
2010 // get crossing
2011
2012    Float_t Dx = fSegmentation->Dx();
2013    Float_t Dz = fSegmentation->Dz();
2014
2015    //P==X
2016    //N==Z
2017
2018     Float_t dx = 0.1;
2019     P *= fPitch;
2020     N *= fPitch;
2021
2022     //P = (kZ * fTan + N + P)/2.0;         // x coordinate
2023     //N = kZ - (P-N)/fTan;                 // z coordinate
2024
2025     if(fTanP + fTanN == 0) return kFALSE;
2026
2027     N = (N - P + fTanN * Dz) / (fTanP + fTanN);  // X coordinate
2028     P = P + fTanP * N;                           // Y coordinate
2029
2030     P -= Dx/2;
2031     N -= Dz/2;
2032
2033     //   N = -(N - P + kZ/2*(fTanP + fTanN))/(fTanP + fTanN);
2034     //  P = (fTanP*(N-kX/2-kZ*fTanN/2) + fTanN*(P-kX/2-kZ*fTanP/2) )/(fTanP + fTanN);
2035
2036     if ( ( N < -(Dz/2+dx) ) || ( N > (Dz/2+dx) ) ) return kFALSE;
2037     if ( ( P < -(Dx/2+dx) ) || ( P > (Dx/2+dx) ) ) return kFALSE;
2038
2039     return kTRUE;
2040 }
2041
2042
2043 //_________________________________________________________________________
2044
2045 void AliITSClusterFinderSSD::GetCrossingError(Float_t& dP, Float_t& dN)
2046 {
2047 // get crossing error
2048
2049   Float_t dz, dx;
2050
2051   dz = TMath::Abs(( dP + dN )*fPitch/(fTanP + fTanN) );
2052   dx = fPitch*(TMath::Abs(dP*(1 - fTanP/(fTanP + fTanN))) +
2053                TMath::Abs(dN *fTanP/(fTanP + fTanN) ));
2054
2055   dN = dz;
2056   dP = dx;
2057 }