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