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