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