ITS/AliITSClusterFinderSSD.cxx

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