ITS/AliITSClusterFinderSDD.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   $Id$
  17   $Log$
  18   Revision 1.38  2005/02/15 13:39:35  masera
  19   V2 clusterer moved to the standard framework. V2 clusters and recpoints are still different objects for the moment
  20
  21   Revision 1.37  2004/06/10 21:00:24  nilsen
  22   Modifications associated with remerging the Ba/Sa and Dubna pixel simulations,
  23   some cleaning of general code (including coding convensions), and adding some
  24   protections associated with SetDefaults/SetDefaultSimulations which should help
  25   with the Test beam simulations. Details below. The default SPD simulation for
  26   the general ITS runs/geometry is still the Ba/Sa, but for the Test beam
  27   geometries this has been changed to the merged versions.
  28   File: AliITS.cxx                         Modified
  29   File: AliITS.h                           Modified
  30         In lined many one-two line functions. Added some protection to
  31         SetDefaults(), SetDefaultSimulation(), and SetDefaultClusterFinders(),
  32         such that they should now even work when only one detector type has
  33         been defined (as it should be for the test beams...). Some mostly
  34         cosmetic issues associated with getting branch names for digits. And
  35         Generally some cleaning up of the code.
  36   File: AliITSClusterFinder.cxx            Modified
  37   File: AliITSClusterFinder.h              Modified
  38         Did some additional consolidation of data into the base class, added
  39         TClonesArray *fClusters, a fDebug, and fModule variables. Otherwise
  40         some cosmetic and coding conversion changes.
  41   File: AliITSClusterFinderSDD.cxx         Modified
  42   File: AliITSClusterFinderSDD.h           Modified
  43         Changes to be consistent with the modified base class, and cosmetic
  44         and coding conversion changes.
  45   File: AliITSClusterFinderSPD.cxx         Modified
  46   File: AliITSClusterFinderSPD.h           Modified
  47         Changes to be consistent with the modified base class, and cosmetic
  48         and coding conversion changes.
  49   File: AliITSClusterFinderSPDdubna.h       Removed
  50   File: AliITSClusterFinderSPDdubna.cxx     Removed
  51         Since we have ClusterFinderSPD and V2 and this version isn't being
  52         maintained, it is being retired.
  53   File: AliITSClusterFinderSSD.cxx         Modified
  54   File: AliITSClusterFinderSSD.h           Modified
  55         Changes to be consistent with the modified base class, and cosmetic
  56         and coding conversion changes.
  57   File: AliITSDetType.cxx                  Modified
  58   File: AliITSDetType.h                    Modified
  59         Added a new class variable to indicate what the detector type is
  60         AliITSDetector fDetType;  values of kSPD, kSDD, kSSD, .... Otherwise
  61         cosmetic and Coding convention changes.
  62   File: AliITSLoader.cxx                   Modified
  63   File: AliITSLoader.h                     Modified
  64         Some changes which are not complete. The idea is to be able to get,
  65         simply via one call, a specific hit, Sdigit, digit, RecPoint,...
  66         without all of the usual over head of initializing TClonesArrays setting
  67         branch addresses and the like. Work is far form ready.
  68   File: AliITSdcsSSD.cxx                   Modified
  69         Some nearly cosmetic changes necessary due to changes to response and
  70         segmentation class'.
  71   File: AliITSgeom.h                       Modified
  72         In the definition of AliITSDetector type, added kND=-1, no detector
  73         defined. Expect to use it later(?).
  74   File: AliITSresponse.h                   Modified
  75         Basically cosmetic. Mostly changing Float_t to Double_t.
  76   File: AliITSresponseSDD.cxx              Modified
  77   File: AliITSresponseSDD.h                Modified
  78         Basically the cosmetic and Float_t to Double_t
  79   File: AliITSresponseSPD.cxx              Modified
  80   File: AliITSresponseSPD.h                Modified
  81         Mostly Float_t to Double_t and added in the IsPixelDead function for
  82         the dubna version (otherwise the merging had been done).
  83   File: AliITSresponseSPDdubna.h           Removed
  84   File: AliITSresponseSPDdubna.cxx         Removed
  85         We should be able to remove this class now. AliITSresponseSPD is now
  86         used for both the Bari-Salerno and the dubna models.
  87   File: AliITSresponseSSD.cxx              Modified
  88   File: AliITSresponseSSD.h                Modified
  89         Float_t to Double_t changes.
  90   File: AliITSsegmentation.h               Modified
  91         Made LocaltoDet return a Bool_t. Now if the x,z location is outside
  92         of the volume, it returns kFALSE. see below.
  93   File: AliITSsegmentationSDD.cxx          Modified
  94   File: AliITSsegmentationSDD.h            Modified
  95         Made LocaltoDet return a Bool_t. Now if the x,z location is outside
  96         of the volume, it returns kFALSE.
  97   File: AliITSsegmentationSPD.cxx          Modified
  98   File: AliITSsegmentationSPD.h            Modified
  99         Made LocaltoDet return a Bool_t. Now if the x,z location is outside
 100         of the volume, it returns kFALSE.
 101   File: AliITSsegmentationSSD.cxx          Modified
 102   File: AliITSsegmentationSSD.h            Modified
 103         Made LocaltoDet return a Bool_t. Now if the x,z location is outside
 104         of the volume, it returns kFALSE. see below.
 105   File: AliITSsimulation.cxx               Modified
 106   File: AliITSsimulation.h                 Modified
 107         Added fDebug variable, new Constructor for use below. Cosmetic and
 108         coding convention changes
 109   File: AliITSsimulationSDD.cxx            Modified
 110   File: AliITSsimulationSDD.h              Modified
 111         Added new Constructor, removed redundant variables and Cosmetic and
 112         coding convention changes.
 113   File: AliITSsimulationSPD.cxx            Modified
 114   File: AliITSsimulationSPD.h              Modified
 115         Removed some dead code, made changes as needed by the changes above
 116         (response and segmentation classes...). a few cosmetic and coding
 117         convention changes.
 118   File: AliITSsimulationSPDdubna.cxx       Modified
 119   File: AliITSsimulationSPDdubna.h         Modified
 120         New merged version, implemented new and old coupling with switch,
 121         coding convention and similar changes. (found 1 bugs, missing
 122         ! in front of if(mod-LineSegmentL(....,).
 123   File: AliITSsimulationSSD.cxx            Modified
 124   File: AliITSsimulationSSD.h              Modified
 125         removed redundant variables with base class. Fixed for coding
 126         convention and other cosmetic changes.
 127   File: AliITSvSDD03.cxx                   Modified
 128   File: AliITSvSPD02.cxx                   Modified
 129   File: AliITSvSSD03.cxx                   Modified
 130         These two have their private versions of SetDefaults and
 131         SetDefaultSimulation which have been similarly protected as in AliITS.cxx
 132   File: ITSLinkDef.h                       Modified
 133   File: libITS.pkg                         Modified
 134         Versions which include v11 geometry and other private changes
 135
 136   Revision 1.36  2004/01/27 16:12:03  masera
 137   Coding conventions for AliITSdigitXXX classes and AliITSTrackerV1
 138
 139   Revision 1.35  2003/11/10 16:33:50  masera
 140   Changes to obey our coding conventions
 141
 142   Revision 1.34  2003/09/11 13:48:52  masera
 143   Data members of AliITSdigit classes defined as protected (They were public)
 144
 145   Revision 1.33  2003/07/21 14:20:51  masera
 146   Fix to track labes in SDD Rec-points
 147
 148   Revision 1.31.2.1  2003/07/16 13:18:04  masera
 149   Proper fix to track labels associated to SDD rec-points
 150
 151   Revision 1.31  2003/05/19 14:44:41  masera
 152   Fix to track labels associated to SDD rec-points
 153
 154   Revision 1.30  2003/03/03 16:34:35  masera
 155   Corrections to comply with coding conventions
 156
 157   Revision 1.29  2002/10/25 18:54:22  barbera
 158   Various improvements and updates from B.S.Nilsen and T. Virgili
 159
 160   Revision 1.28  2002/10/22 14:45:29  alibrary
 161   Introducing Riostream.h
 162
 163   Revision 1.27  2002/10/14 14:57:00  hristov
 164   Merging the VirtualMC branch to the main development branch (HEAD)
 165
 166   Revision 1.23.4.2  2002/10/14 13:14:07  hristov
 167   Updating VirtualMC to v3-09-02
 168
 169   Revision 1.26  2002/09/09 17:23:28  nilsen
 170   Minor changes in support of changes to AliITSdigitS?D class'.
 171
 172   Revision 1.25  2002/05/10 22:29:40  nilsen
 173   Change my Massimo Masera in the default constructor to bring things into
 174   compliance.
 175
 176   Revision 1.24  2002/04/24 22:02:31  nilsen
 177   New SDigits and Digits routines, and related changes,  (including new
 178   noise values).
 179
 180  */
 181 ///////////////////////////////////////////////////////////////////////////
 182 //  Cluster finder                                                       //
 183 //  for Silicon                                                          //
 184 //  Drift Detector                                                       //
 185 //////////////////////////////////////////////////////////////////////////
 186 #include <Riostream.h>
 187
 188 #include <TMath.h>
 189 #include <math.h>
 190
 191 #include "AliITSClusterFinderSDD.h"
 192 #include "AliITSMapA1.h"
 193 #include "AliITS.h"
 194 #include "AliITSdigitSDD.h"
 195 #include "AliITSRawClusterSDD.h"
 196 #include "AliITSRecPoint.h"
 197 #include "AliITSsegmentationSDD.h"
 198 #include "AliITSresponseSDD.h"
 199 #include "AliRun.h"
 200
 201 ClassImp(AliITSClusterFinderSDD)
 202
 203 //______________________________________________________________________
 204 AliITSClusterFinderSDD::AliITSClusterFinderSDD():
 205 AliITSClusterFinder(),
 206 fNclusters(0),
 207 fDAnode(0.0),
 208 fDTime(0.0),
 209 fTimeCorr(0.0),
 210 fCutAmplitude(0),
 211 fMinPeak(0),
 212 fMinCharge(0),
 213 fMinNCells(0),
 214 fMaxNCells(0){
 215     // default constructor
 216 }
 217 //______________________________________________________________________
 218 AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSsegmentation *seg,
 219                                                AliITSresponse *response,
 220                                                TClonesArray *digits,
 221                                                TClonesArray *recp):
 222 AliITSClusterFinder(seg,response),
 223 fNclusters(0),
 224 fDAnode(0.0),
 225 fDTime(0.0),
 226 fTimeCorr(0.0),
 227 fCutAmplitude(0),
 228 fMinPeak(0),
 229 fMinCharge(0),
 230 fMinNCells(0),
 231 fMaxNCells(0){
 232     // standard constructor
 233
 234     SetDigits(digits);
 235     SetClusters(recp);
 236     SetCutAmplitude();
 237     SetDAnode();
 238     SetDTime();
 239     SetMinPeak((Int_t)(((AliITSresponseSDD*)GetResp())->
 240                        GetNoiseAfterElectronics()*5));
 241     //    SetMinPeak();
 242     SetMinNCells();
 243     SetMaxNCells();
 244     SetTimeCorr();
 245     SetMinCharge();
 246     SetMap(new AliITSMapA1(GetSeg(),Digits(),fCutAmplitude));
 247 }
 248 //______________________________________________________________________
 249 void AliITSClusterFinderSDD::SetCutAmplitude(Double_t nsigma){
 250     // set the signal threshold for cluster finder
 251     Double_t baseline,noise,noiseAfterEl;
 252
 253     GetResp()->GetNoiseParam(noise,baseline);
 254     noiseAfterEl = ((AliITSresponseSDD*)GetResp())->GetNoiseAfterElectronics();
 255     fCutAmplitude = (Int_t)((baseline + nsigma*noiseAfterEl));
 256 }
 257 //______________________________________________________________________
 258 void AliITSClusterFinderSDD::Find1DClusters(){
 259     // find 1D clusters
 260
 261     // retrieve the parameters
 262     Int_t fNofMaps       = GetSeg()->Npz();
 263     Int_t fMaxNofSamples = GetSeg()->Npx();
 264     Int_t fNofAnodes     = fNofMaps/2;
 265     Int_t dummy          = 0;
 266     Double_t fTimeStep    = GetSeg()->Dpx(dummy);
 267     Double_t fSddLength   = GetSeg()->Dx();
 268     Double_t fDriftSpeed  = GetResp()->DriftSpeed();
 269     Double_t anodePitch   = GetSeg()->Dpz(dummy);
 270
 271     // map the signal
 272     Map()->ClearMap();
 273     Map()->SetThreshold(fCutAmplitude);
 274     Map()->FillMap();
 275
 276     Double_t noise;
 277     Double_t baseline;
 278     GetResp()->GetNoiseParam(noise,baseline);
 279
 280     Int_t nofFoundClusters = 0;
 281     Int_t i;
 282     Double_t **dfadc = new Double_t*[fNofAnodes];
 283     for(i=0;i<fNofAnodes;i++) dfadc[i] = new Double_t[fMaxNofSamples];
 284     Double_t fadc  = 0.;
 285     Double_t fadc1 = 0.;
 286     Double_t fadc2 = 0.;
 287     Int_t j,k,idx,l,m;
 288     for(j=0;j<2;j++) {
 289         for(k=0;k<fNofAnodes;k++) {
 290             idx = j*fNofAnodes+k;
 291             // signal (fadc) & derivative (dfadc)
 292             dfadc[k][255]=0.;
 293             for(l=0; l<fMaxNofSamples; l++) {
 294                 fadc2=(Double_t)Map()->GetSignal(idx,l);
 295                 if(l>0) fadc1=(Double_t)Map()->GetSignal(idx,l-1);
 296                 if(l>0) dfadc[k][l-1] = fadc2-fadc1;
 297             } // samples
 298         } // anodes
 299
 300         for(k=0;k<fNofAnodes;k++) {
 301             if(GetDebug(5)) cout<<"Anode: "<<k+1<<", Wing: "<<j+1<< endl;
 302             idx = j*fNofAnodes+k;
 303             Int_t imax  = 0;
 304             Int_t imaxd = 0;
 305             Int_t it    = 0;
 306             while(it <= fMaxNofSamples-3) {
 307                 imax  = it;
 308                 imaxd = it;
 309                 // maximum of signal
 310                 Double_t fadcmax  = 0.;
 311                 Double_t dfadcmax = 0.;
 312                 Int_t lthrmina   = 1;
 313                 Int_t lthrmint   = 3;
 314                 Int_t lthra      = 1;
 315                 Int_t lthrt      = 0;
 316                 for(m=0;m<20;m++) {
 317                     Int_t id = it+m;
 318                     if(id>=fMaxNofSamples) break;
 319                     fadc=(float)Map()->GetSignal(idx,id);
 320                     if(fadc > fadcmax) { fadcmax = fadc; imax = id;}
 321                     if(fadc > (float)fCutAmplitude)lthrt++;
 322                     if(dfadc[k][id] > dfadcmax) {
 323                         dfadcmax = dfadc[k][id];
 324                         imaxd    = id;
 325                     } // end if
 326                 } // end for m
 327                 it = imaxd;
 328                 if(Map()->TestHit(idx,imax) == kEmpty) {it++; continue;}
 329                 // cluster charge
 330                 Int_t tstart = it-2;
 331                 if(tstart < 0) tstart = 0;
 332                 Bool_t ilcl = 0;
 333                 if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
 334                 if(ilcl) {
 335                     nofFoundClusters++;
 336                     Int_t tstop      = tstart;
 337                     Double_t dfadcmin = 10000.;
 338                     Int_t ij;
 339                     for(ij=0; ij<20; ij++) {
 340                         if(tstart+ij > 255) { tstop = 255; break; }
 341                         fadc=(float)Map()->GetSignal(idx,tstart+ij);
 342                         if((dfadc[k][tstart+ij] < dfadcmin) &&
 343                            (fadc > fCutAmplitude)) {
 344                             tstop = tstart+ij+5;
 345                             if(tstop > 255) tstop = 255;
 346                             dfadcmin = dfadc[k][it+ij];
 347                         } // end if
 348                     } // end for ij
 349
 350                     Double_t clusterCharge = 0.;
 351                     Double_t clusterAnode  = k+0.5;
 352                     Double_t clusterTime   = 0.;
 353                     Int_t   clusterMult   = 0;
 354                     Double_t clusterPeakAmplitude = 0.;
 355                     Int_t its,peakpos     = -1;
 356                     Double_t n, baseline;
 357                     GetResp()->GetNoiseParam(n,baseline);
 358                     for(its=tstart; its<=tstop; its++) {
 359                         fadc=(float)Map()->GetSignal(idx,its);
 360                         if(fadc>baseline) fadc -= baseline;
 361                         else fadc = 0.;
 362                         clusterCharge += fadc;
 363                         // as a matter of fact we should take the peak
 364                         // pos before FFT
 365                         // to get the list of tracks !!!
 366                         if(fadc > clusterPeakAmplitude) {
 367                             clusterPeakAmplitude = fadc;
 368                             //peakpos=Map()->GetHitIndex(idx,its);
 369                             Int_t shift = (int)(fTimeCorr/fTimeStep);
 370                             if(its>shift && its<(fMaxNofSamples-shift))
 371                                 peakpos  = Map()->GetHitIndex(idx,its+shift);
 372                             else peakpos = Map()->GetHitIndex(idx,its);
 373                             if(peakpos<0) peakpos =Map()->GetHitIndex(idx,its);
 374                         } // end if
 375                         clusterTime += fadc*its;
 376                         if(fadc > 0) clusterMult++;
 377                         if(its == tstop) {
 378                             clusterTime /= (clusterCharge/fTimeStep);   // ns
 379                             if(clusterTime>fTimeCorr) clusterTime -=fTimeCorr;
 380                             //ns
 381                         } // end if
 382                     } // end for its
 383
 384                     Double_t clusteranodePath = (clusterAnode - fNofAnodes/2)*
 385                                                  anodePitch;
 386                     Double_t clusterDriftPath = clusterTime*fDriftSpeed;
 387                     clusterDriftPath = fSddLength-clusterDriftPath;
 388                     if(clusterCharge <= 0.) break;
 389                     AliITSRawClusterSDD clust(j+1,//i
 390                                               clusterAnode,clusterTime,//ff
 391                                               clusterCharge, //f
 392                                               clusterPeakAmplitude, //f
 393                                               peakpos, //i
 394                                               0.,0.,clusterDriftPath,//fff
 395                                               clusteranodePath, //f
 396                                               clusterMult, //i
 397                                               0,0,0,0,0,0,0);//7*i
 398                     fITS->AddCluster(1,&clust);
 399                     it = tstop;
 400                 } // ilcl
 401                 it++;
 402             } // while (samples)
 403         } // anodes
 404     } // detectors (2)
 405
 406     for(i=0;i<fNofAnodes;i++) delete[] dfadc[i];
 407     delete [] dfadc;
 408
 409     return;
 410 }
 411 //______________________________________________________________________
 412 void AliITSClusterFinderSDD::Find1DClustersE(){
 413     // find 1D clusters
 414     // retrieve the parameters
 415     Int_t fNofMaps = GetSeg()->Npz();
 416     Int_t fMaxNofSamples = GetSeg()->Npx();
 417     Int_t fNofAnodes = fNofMaps/2;
 418     Int_t dummy=0;
 419     Double_t fTimeStep = GetSeg()->Dpx( dummy );
 420     Double_t fSddLength = GetSeg()->Dx();
 421     Double_t fDriftSpeed = GetResp()->DriftSpeed();
 422     Double_t anodePitch = GetSeg()->Dpz( dummy );
 423     Double_t n, baseline;
 424     GetResp()->GetNoiseParam( n, baseline );
 425     // map the signal
 426     Map()->ClearMap();
 427     Map()->SetThreshold( fCutAmplitude );
 428     Map()->FillMap();
 429
 430     Int_t nClu = 0;
 431     //        cout << "Search  cluster... "<< endl;
 432     for( Int_t j=0; j<2; j++ ){
 433         for( Int_t k=0; k<fNofAnodes; k++ ){
 434             Int_t idx = j*fNofAnodes+k;
 435             Bool_t on = kFALSE;
 436             Int_t start = 0;
 437             Int_t nTsteps = 0;
 438             Double_t fmax = 0.;
 439             Int_t lmax = 0;
 440             Double_t charge = 0.;
 441             Double_t time = 0.;
 442             Double_t anode = k+0.5;
 443             Int_t peakpos = -1;
 444             for( Int_t l=0; l<fMaxNofSamples; l++ ){
 445                 Double_t fadc = (Double_t)Map()->GetSignal( idx, l );
 446                 if( fadc > 0.0 ){
 447                     if( on == kFALSE && l<fMaxNofSamples-4){
 448                         // star RawCluster (reset var.)
 449                         Double_t fadc1 = (Double_t)Map()->GetSignal( idx, l+1 );
 450                         if( fadc1 < fadc ) continue;
 451                         start = l;
 452                         fmax = 0.;
 453                         lmax = 0;
 454                         time = 0.;
 455                         charge = 0.;
 456                         on = kTRUE;
 457                         nTsteps = 0;
 458                     } // end if on...
 459                     nTsteps++ ;
 460                     if( fadc > baseline ) fadc -= baseline;
 461                     else fadc=0.;
 462                     charge += fadc;
 463                     time += fadc*l;
 464                     if( fadc > fmax ){
 465                         fmax = fadc;
 466                         lmax = l;
 467                         Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5);
 468                         if( l > shift && l < (fMaxNofSamples-shift) )
 469                             peakpos = Map()->GetHitIndex( idx, l+shift );
 470                         else
 471                             peakpos = Map()->GetHitIndex( idx, l );
 472                         if( peakpos < 0) peakpos = Map()->GetHitIndex(idx,l);
 473                     } // end if fadc
 474                 }else{ // end fadc>0
 475                     if( on == kTRUE ){
 476                         if( nTsteps > 2 ){
 477                             //  min # of timesteps for a RawCluster
 478                             // Found a RawCluster...
 479                             Int_t stop = l-1;
 480                             time /= (charge/fTimeStep);   // ns
 481                                 // time = lmax*fTimeStep;   // ns
 482                             if( time > fTimeCorr ) time -= fTimeCorr;   // ns
 483                             Double_t anodePath =(anode-fNofAnodes/2)*anodePitch;
 484                             Double_t driftPath = time*fDriftSpeed;
 485                             driftPath = fSddLength-driftPath;
 486                             AliITSRawClusterSDD clust(j+1,anode,time,charge,
 487                                                       fmax, peakpos,0.,0.,
 488                                                       driftPath,anodePath,
 489                                                       nTsteps,start,stop,
 490                                                       start, stop, 1, k, k );
 491                             fITS->AddCluster( 1, &clust );
 492                             if(GetDebug(5)) clust.PrintInfo();
 493                             nClu++;
 494                         } // end if nTsteps
 495                         on = kFALSE;
 496                     } // end if on==kTRUE
 497                 } // end if fadc>0
 498             } // samples
 499         } // anodes
 500     } // wings
 501     if(GetDebug(3)) cout << "# Rawclusters " << nClu << endl;
 502     return;
 503 }
 504 //_______________________________________________________________________
 505 Int_t AliITSClusterFinderSDD::SearchPeak(Double_t *spect,Int_t xdim,Int_t zdim,
 506                                          Int_t *peakX, Int_t *peakZ,
 507                                          Double_t *peakAmp, Double_t minpeak ){
 508     // search peaks on a 2D cluster
 509     Int_t npeak = 0;    // # peaks
 510     Int_t i,j;
 511     // search peaks
 512     for( Int_t z=1; z<zdim-1; z++ ){
 513         for( Int_t x=1; x<xdim-2; x++ ){
 514             Double_t sxz = spect[x*zdim+z];
 515             Double_t sxz1 = spect[(x+1)*zdim+z];
 516             Double_t sxz2 = spect[(x-1)*zdim+z];
 517             // search a local max. in s[x,z]
 518             if( sxz < minpeak || sxz1 <= 0 || sxz2 <= 0 ) continue;
 519             if( sxz >= spect[(x+1)*zdim+z  ] && sxz >= spect[(x-1)*zdim+z  ] &&
 520                 sxz >= spect[x*zdim    +z+1] && sxz >= spect[x*zdim    +z-1] &&
 521                 sxz >= spect[(x+1)*zdim+z+1] && sxz >= spect[(x+1)*zdim+z-1] &&
 522                 sxz >= spect[(x-1)*zdim+z+1] && sxz >= spect[(x-1)*zdim+z-1] ){
 523                 // peak found
 524                 peakX[npeak] = x;
 525                 peakZ[npeak] = z;
 526                 peakAmp[npeak] = sxz;
 527                 npeak++;
 528             } // end if ....
 529         } // end for x
 530     } // end for z
 531     // search groups of peaks with same amplitude.
 532     Int_t *flag = new Int_t[npeak];
 533     for( i=0; i<npeak; i++ ) flag[i] = 0;
 534     for( i=0; i<npeak; i++ ){
 535         for( j=0; j<npeak; j++ ){
 536             if( i==j) continue;
 537             if( flag[j] > 0 ) continue;
 538             if( peakAmp[i] == peakAmp[j] &&
 539                 TMath::Abs(peakX[i]-peakX[j])<=1 &&
 540                 TMath::Abs(peakZ[i]-peakZ[j])<=1 ){
 541                 if( flag[i] == 0) flag[i] = i+1;
 542                 flag[j] = flag[i];
 543             } // end if ...
 544         } // end for j
 545     } // end for i
 546     // make average of peak groups
 547     for( i=0; i<npeak; i++ ){
 548         Int_t nFlag = 1;
 549         if( flag[i] <= 0 ) continue;
 550         for( j=0; j<npeak; j++ ){
 551             if( i==j ) continue;
 552             if( flag[j] != flag[i] ) continue;
 553             peakX[i] += peakX[j];
 554             peakZ[i] += peakZ[j];
 555             nFlag++;
 556             npeak--;
 557             for( Int_t k=j; k<npeak; k++ ){
 558                 peakX[k] = peakX[k+1];
 559                 peakZ[k] = peakZ[k+1];
 560                 peakAmp[k] = peakAmp[k+1];
 561                 flag[k] = flag[k+1];
 562             } // end for k
 563             j--;
 564         } // end for j
 565         if( nFlag > 1 ){
 566             peakX[i] /= nFlag;
 567             peakZ[i] /= nFlag;
 568         } // end fi nFlag
 569     } // end for i
 570     delete [] flag;
 571     return( npeak );
 572 }
 573 //______________________________________________________________________
 574 void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Double_t *par,
 575                                        Double_t *spe, Double_t *integral){
 576     // function used to fit the clusters
 577     // par -> parameters..
 578     // par[0]  number of peaks.
 579     // for each peak i=1, ..., par[0]
 580     //                 par[i] = Ampl.
 581     //                 par[i+1] = xpos
 582     //                 par[i+2] = zpos
 583     //                 par[i+3] = tau
 584     //                 par[i+4] = sigma.
 585     Int_t electronics = GetResp()->Electronics(); // 1 = PASCAL, 2 = OLA
 586     const Int_t knParam = 5;
 587     Int_t npeak = (Int_t)par[0];
 588
 589     memset( spe, 0, sizeof( Double_t )*zdim*xdim );
 590
 591     Int_t k = 1;
 592     for( Int_t i=0; i<npeak; i++ ){
 593         if( integral != 0 ) integral[i] = 0.;
 594         Double_t sigmaA2 = par[k+4]*par[k+4]*2.;
 595         Double_t t2 = par[k+3];   // PASCAL
 596         if( electronics == 2 ) { t2 *= t2; t2 *= 2; } // OLA
 597         for( Int_t z=0; z<zdim; z++ ){
 598             for( Int_t x=0; x<xdim; x++ ){
 599                 Double_t z2 = (z-par[k+2])*(z-par[k+2])/sigmaA2;
 600                 Double_t x2 = 0.;
 601                 Double_t signal = 0.;
 602                 if( electronics == 1 ){ // PASCAL
 603                     x2 = (x-par[k+1]+t2)/t2;
 604                     signal = (x2>0.) ? par[k]*x2*exp(-x2+1.-z2) :0.0; // RCCR2
 605                 //  signal =(x2>0.) ? par[k]*x2*x2*exp(-2*x2+2.-z2 ):0.0;//RCCR
 606                 }else if( electronics == 2 ) { // OLA
 607                     x2 = (x-par[k+1])*(x-par[k+1])/t2;
 608                     signal = par[k]  * exp( -x2 - z2 );
 609                 } else {
 610                     Warning("PeakFunc","Wrong SDD Electronics = %d",
 611                             electronics);
 612                     // exit( 1 );
 613                 } // end if electronicx
 614                 spe[x*zdim+z] += signal;
 615                 if( integral != 0 ) integral[i] += signal;
 616             } // end for x
 617         } // end for z
 618         k += knParam;
 619     } // end for i
 620     return;
 621 }
 622 //__________________________________________________________________________
 623 Double_t AliITSClusterFinderSDD::ChiSqr( Int_t xdim, Int_t zdim, Double_t *spe,
 624                                         Double_t *speFit ) const{
 625     // EVALUATES UNNORMALIZED CHI-SQUARED
 626     Double_t chi2 = 0.;
 627     for( Int_t z=0; z<zdim; z++ ){
 628         for( Int_t x=1; x<xdim-1; x++ ){
 629             Int_t index = x*zdim+z;
 630             Double_t tmp = spe[index] - speFit[index];
 631             chi2 += tmp*tmp;
 632         } // end for x
 633     } // end for z
 634     return( chi2 );
 635 }
 636 //_______________________________________________________________________
 637 void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Double_t *param,
 638                                     Double_t *prm0,Double_t *steprm,
 639                                     Double_t *chisqr,Double_t *spe,
 640                                     Double_t *speFit ){
 641     //
 642     Int_t   k, nnn, mmm, i;
 643     Double_t p1, delta, d1, chisq1, p2, chisq2, t, p3, chisq3, a, b, p0, chisqt;
 644     const Int_t knParam = 5;
 645     Int_t npeak = (Int_t)param[0];
 646     for( k=1; k<(npeak*knParam+1); k++ ) prm0[k] = param[k];
 647     for( k=1; k<(npeak*knParam+1); k++ ){
 648         p1 = param[k];
 649         delta = steprm[k];
 650         d1 = delta;
 651         // ENSURE THAT STEP SIZE IS SENSIBLY LARGER THAN MACHINE ROUND OFF
 652         if( fabs( p1 ) > 1.0E-6 )
 653             if ( fabs( delta/p1 ) < 1.0E-4 ) delta = p1/1000;
 654             else  delta = (Double_t)1.0E-4;
 655         //  EVALUATE CHI-SQUARED AT FIRST TWO SEARCH POINTS
 656         PeakFunc( xdim, zdim, param, speFit );
 657         chisq1 = ChiSqr( xdim, zdim, spe, speFit );
 658         p2 = p1+delta;
 659         param[k] = p2;
 660         PeakFunc( xdim, zdim, param, speFit );
 661         chisq2 = ChiSqr( xdim, zdim, spe, speFit );
 662         if( chisq1 < chisq2 ){
 663             // REVERSE DIRECTION OF SEARCH IF CHI-SQUARED IS INCREASING
 664             delta = -delta;
 665             t = p1;
 666             p1 = p2;
 667             p2 = t;
 668             t = chisq1;
 669             chisq1 = chisq2;
 670             chisq2 = t;
 671         } // end if
 672         i = 1; nnn = 0;
 673         do {   // INCREMENT param(K) UNTIL CHI-SQUARED STARTS TO INCREASE
 674             nnn++;
 675             p3 = p2 + delta;
 676             mmm = nnn - (nnn/5)*5;  // multiplo de 5
 677             if( mmm == 0 ){
 678                 d1 = delta;
 679                 // INCREASE STEP SIZE IF STEPPING TOWARDS MINIMUM IS TOO SLOW
 680                 delta *= 5;
 681             } // end if
 682             param[k] = p3;
 683             // Constrain paramiters
 684             Int_t kpos = (k-1) % knParam;
 685             switch( kpos ){
 686             case 0 :
 687                 if( param[k] <= 20 ) param[k] = fMinPeak;
 688                 break;
 689             case 1 :
 690                 if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k];
 691                 break;
 692             case 2 :
 693                 if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k];
 694                 break;
 695             case 3 :
 696                 if( param[k] < .5 ) param[k] = .5;
 697                 break;
 698             case 4 :
 699                 if( param[k] < .288 ) param[k] = .288;// 1/sqrt(12) = 0.288
 700                 if( param[k] > zdim*.5 ) param[k] = zdim*.5;
 701                 break;
 702             }; // end switch
 703             PeakFunc( xdim, zdim, param, speFit );
 704             chisq3 = ChiSqr( xdim, zdim, spe, speFit );
 705             if( chisq3 < chisq2 && nnn < 50 ){
 706                 p1 = p2;
 707                 p2 = p3;
 708                 chisq1 = chisq2;
 709                 chisq2 = chisq3;
 710             }else i=0;
 711         } while( i );
 712         // FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS
 713         a = chisq1*(p2-p3)+chisq2*(p3-p1)+chisq3*(p1-p2);
 714         b = chisq1*(p2*p2-p3*p3)+chisq2*(p3*p3-p1*p1)+chisq3*(p1*p1-p2*p2);
 715         if( a!=0 ) p0 = (Double_t)(0.5*b/a);
 716         else p0 = 10000;
 717         //--IN CASE OF NEARLY EQUAL CHI-SQUARED AND TOO SMALL STEP SIZE PREVENT
 718         //   ERRONEOUS EVALUATION OF PARABOLA MINIMUM
 719         //---NEXT TWO LINES CAN BE OMITTED FOR HIGHER PRECISION MACHINES
 720         //dp = (Double_t) max (fabs(p3-p2), fabs(p2-p1));
 721         //if( fabs( p2-p0 ) > dp ) p0 = p2;
 722         param[k] = p0;
 723         // Constrain paramiters
 724         Int_t kpos = (k-1) % knParam;
 725         switch( kpos ){
 726         case 0 :
 727             if( param[k] <= 20 ) param[k] = fMinPeak;
 728             break;
 729         case 1 :
 730             if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k];
 731             break;
 732         case 2 :
 733             if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k];
 734             break;
 735         case 3 :
 736             if( param[k] < .5 ) param[k] = .5;
 737             break;
 738         case 4 :
 739             if( param[k] < .288 ) param[k] = .288;  // 1/sqrt(12) = 0.288
 740             if( param[k] > zdim*.5 ) param[k] = zdim*.5;
 741             break;
 742         }; // end switch
 743         PeakFunc( xdim, zdim, param, speFit );
 744         chisqt = ChiSqr( xdim, zdim, spe, speFit );
 745         // DO NOT ALLOW ERRONEOUS INTERPOLATION
 746         if( chisqt <= *chisqr ) *chisqr = chisqt;
 747         else param[k] = prm0[k];
 748         // OPTIMIZE SEARCH STEP FOR EVENTUAL NEXT CALL OF MINIM
 749         steprm[k] = (param[k]-prm0[k])/5;
 750         if( steprm[k] >= d1 ) steprm[k] = d1/5;
 751     } // end for k
 752     // EVALUATE FIT AND CHI-SQUARED FOR OPTIMIZED PARAMETERS
 753     PeakFunc( xdim, zdim, param, speFit );
 754     *chisqr = ChiSqr( xdim, zdim, spe, speFit );
 755     return;
 756 }
 757 //_________________________________________________________________________
 758 Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim,
 759                                            Double_t *param, Double_t *spe,
 760                                            Int_t *niter, Double_t *chir ){
 761     // fit method from Comput. Phys. Commun 46(1987) 149
 762     const Double_t kchilmt = 0.01;  //        relative accuracy
 763     const Int_t   knel = 3;        //        for parabolic minimization
 764     const Int_t   knstop = 50;     //        Max. iteration number
 765     const Int_t   knParam = 5;
 766     Int_t npeak = (Int_t)param[0];
 767     // RETURN IF NUMBER OF DEGREES OF FREEDOM IS NOT POSITIVE
 768     if( (xdim*zdim - npeak*knParam) <= 0 ) return( -1 );
 769     Double_t degFree = (xdim*zdim - npeak*knParam)-1;
 770     Int_t   n, k, iterNum = 0;
 771     Double_t *prm0 = new Double_t[npeak*knParam+1];
 772     Double_t *step = new Double_t[npeak*knParam+1];
 773     Double_t *schi = new Double_t[npeak*knParam+1];
 774     Double_t *sprm[3];
 775     sprm[0] = new Double_t[npeak*knParam+1];
 776     sprm[1] = new Double_t[npeak*knParam+1];
 777     sprm[2] = new Double_t[npeak*knParam+1];
 778     Double_t  chi0, chi1, reldif, a, b, prmin, dp;
 779     Double_t *speFit = new Double_t[ xdim*zdim ];
 780     PeakFunc( xdim, zdim, param, speFit );
 781     chi0 = ChiSqr( xdim, zdim, spe, speFit );
 782     chi1 = chi0;
 783     for( k=1; k<(npeak*knParam+1); k++) prm0[k] = param[k];
 784         for( k=1 ; k<(npeak*knParam+1); k+=knParam ){
 785             step[k] = param[k] / 20.0 ;
 786             step[k+1] = param[k+1] / 50.0;
 787             step[k+2] = param[k+2] / 50.0;
 788             step[k+3] = param[k+3] / 20.0;
 789             step[k+4] = param[k+4] / 20.0;
 790         } // end for k
 791     Int_t out = 0;
 792     do{
 793         iterNum++;
 794             chi0 = chi1;
 795             Minim( xdim, zdim, param, prm0, step, &chi1, spe, speFit );
 796             reldif = ( chi1 > 0 ) ? ((Double_t) fabs( chi1-chi0)/chi1 ) : 0;
 797         // EXIT conditions
 798         if( reldif < (float) kchilmt ){
 799             *chir  = (chi1>0) ? (float) TMath::Sqrt (chi1/degFree) :0;
 800             *niter = iterNum;
 801             out = 0;
 802             break;
 803         } // end if
 804         if( (reldif < (float)(5*kchilmt)) && (iterNum > knstop) ){
 805             *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/degFree):0;
 806             *niter = iterNum;
 807             out = 0;
 808             break;
 809         } // end if
 810         if( iterNum > 5*knstop ){
 811             *chir  = (chi1>0) ?(float) TMath::Sqrt (chi1/degFree):0;
 812             *niter = iterNum;
 813             out = 1;
 814             break;
 815         } // end if
 816         if( iterNum <= knel ) continue;
 817         n = iterNum - (iterNum/knel)*knel; // EXTRAPOLATION LIMIT COUNTER N
 818         if( n > 3 || n == 0 ) continue;
 819         schi[n-1] = chi1;
 820         for( k=1; k<(npeak*knParam+1); k++ ) sprm[n-1][k] = param[k];
 821         if( n != 3 ) continue;
 822         // -EVALUATE EXTRAPOLATED VALUE OF EACH PARAMETER BY FINDING MINIMUM OF
 823         //    PARABOLA DEFINED BY LAST THREE CALLS OF MINIM
 824         for( k=1; k<(npeak*knParam+1); k++ ){
 825             Double_t tmp0 = sprm[0][k];
 826             Double_t tmp1 = sprm[1][k];
 827             Double_t tmp2 = sprm[2][k];
 828             a  = schi[0]*(tmp1-tmp2) + schi[1]*(tmp2-tmp0);
 829             a += (schi[2]*(tmp0-tmp1));
 830             b  = schi[0]*(tmp1*tmp1-tmp2*tmp2);
 831             b += (schi[1]*(tmp2*tmp2-tmp0*tmp0)+(schi[2]*
 832                                              (tmp0*tmp0-tmp1*tmp1)));
 833             if ((double)a < 1.0E-6) prmin = 0;
 834             else prmin = (float) (0.5*b/a);
 835             dp = 5*(tmp2-tmp0);
 836             if( fabs(prmin-tmp2) > fabs(dp) ) prmin = tmp2+dp;
 837             param[k] = prmin;
 838             step[k]  = dp/10; // OPTIMIZE SEARCH STEP
 839         } // end for k
 840     } while( kTRUE );
 841     delete [] prm0;
 842     delete [] step;
 843     delete [] schi;
 844     delete [] sprm[0];
 845     delete [] sprm[1];
 846     delete [] sprm[2];
 847     delete [] speFit;
 848     return( out );
 849 }
 850
 851 //______________________________________________________________________
 852 void AliITSClusterFinderSDD::ResolveClusters(){
 853     // The function to resolve clusters if the clusters overlapping exists
 854     Int_t i;
 855     // get number of clusters for this module
 856     Int_t nofClusters = NClusters();
 857     nofClusters -= fNclusters;
 858     Int_t fNofMaps = GetSeg()->Npz();
 859     Int_t fNofAnodes = fNofMaps/2;
 860     //Int_t fMaxNofSamples = GetSeg()->Npx();
 861     Int_t dummy=0;
 862     Double_t fTimeStep = GetSeg()->Dpx( dummy );
 863     Double_t fSddLength = GetSeg()->Dx();
 864     Double_t fDriftSpeed = GetResp()->DriftSpeed();
 865     Double_t anodePitch = GetSeg()->Dpz( dummy );
 866     Double_t n, baseline;
 867     GetResp()->GetNoiseParam( n, baseline );
 868     Int_t electronics = GetResp()->Electronics(); // 1 = PASCAL, 2 = OLA
 869
 870     for( Int_t j=0; j<nofClusters; j++ ){
 871         // get cluster information
 872         AliITSRawClusterSDD *clusterJ=(AliITSRawClusterSDD*) Cluster(j);
 873         Int_t astart = clusterJ->Astart();
 874         Int_t astop = clusterJ->Astop();
 875         Int_t tstart = clusterJ->Tstartf();
 876         Int_t tstop = clusterJ->Tstopf();
 877         Int_t wing = (Int_t)clusterJ->W();
 878         if( wing == 2 ){
 879             astart += fNofAnodes;
 880             astop  += fNofAnodes;
 881         } // end if
 882         Int_t xdim = tstop-tstart+3;
 883         Int_t zdim = astop-astart+3;
 884         if( xdim > 50 || zdim > 30 ) {
 885             Warning("ResolveClusters","xdim: %d , zdim: %d ",xdim,zdim);
 886             continue;
 887         }
 888         Double_t *sp = new Double_t[ xdim*zdim+1 ];
 889         memset( sp, 0, sizeof(Double_t)*(xdim*zdim+1) );
 890
 891         // make a local map from cluster region
 892         for( Int_t ianode=astart; ianode<=astop; ianode++ ){
 893             for( Int_t itime=tstart; itime<=tstop; itime++ ){
 894                 Double_t fadc = Map()->GetSignal( ianode, itime );
 895                 if( fadc > baseline ) fadc -= (Double_t)baseline;
 896                 else fadc = 0.;
 897                 Int_t index = (itime-tstart+1)*zdim+(ianode-astart+1);
 898                 sp[index] = fadc;
 899             } // time loop
 900         } // anode loop
 901
 902         // search peaks on cluster
 903         const Int_t kNp = 150;
 904         Int_t peakX1[kNp];
 905         Int_t peakZ1[kNp];
 906         Double_t peakAmp1[kNp];
 907         Int_t npeak = SearchPeak(sp,xdim,zdim,peakX1,peakZ1,peakAmp1,fMinPeak);
 908
 909         // if multiple peaks, split cluster
 910         if( npeak >= 1 ){
 911             //        cout << "npeak " << npeak << endl;
 912             //        clusterJ->PrintInfo();
 913             Double_t *par = new Double_t[npeak*5+1];
 914             par[0] = (Double_t)npeak;
 915             // Initial parameters in cell dimentions
 916             Int_t k1 = 1;
 917             for( i=0; i<npeak; i++ ){
 918                 par[k1] = peakAmp1[i];
 919                 par[k1+1] = peakX1[i]; // local time pos. [timebin]
 920                 par[k1+2] = peakZ1[i]; // local anode pos. [anodepitch]
 921                 if( electronics == 1 ) par[k1+3] = 2.; // PASCAL
 922                 else if(electronics==2) par[k1+3] = 0.7;//tau [timebin] OLA
 923                 par[k1+4] = .4;    // sigma        [anodepich]
 924                 k1 += 5;
 925             } // end for i
 926             Int_t niter;
 927             Double_t chir;
 928             NoLinearFit( xdim, zdim, par, sp, &niter, &chir );
 929             Double_t peakX[kNp];
 930             Double_t peakZ[kNp];
 931             Double_t sigma[kNp];
 932             Double_t tau[kNp];
 933             Double_t peakAmp[kNp];
 934             Double_t integral[kNp];
 935             //get integrals => charge for each peak
 936             PeakFunc( xdim, zdim, par, sp, integral );
 937             k1 = 1;
 938             for( i=0; i<npeak; i++ ){
 939                 peakAmp[i] = par[k1];
 940                 peakX[i]   = par[k1+1];
 941                 peakZ[i]   = par[k1+2];
 942                 tau[i]     = par[k1+3];
 943                 sigma[i]   = par[k1+4];
 944                 k1+=5;
 945             } // end for i
 946             // calculate parameter for new clusters
 947             for( i=0; i<npeak; i++ ){
 948                 AliITSRawClusterSDD clusterI( *clusterJ );
 949
 950                 Int_t newAnode = peakZ1[i]-1 + astart;
 951
 952             //    Int_t newiTime = peakX1[i]-1 + tstart;
 953             //    Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5);
 954             //    if( newiTime > shift && newiTime < (fMaxNofSamples-shift) )
 955             //        shift = 0;
 956             //    Int_t peakpos = Map()->GetHitIndex(newAnode,newiTime+shift );
 957             //    clusterI.SetPeakPos( peakpos );
 958
 959                 clusterI.SetPeakAmpl( peakAmp1[i] );
 960                 Double_t newAnodef = peakZ[i] - 0.5 + astart;
 961                 Double_t newiTimef = peakX[i] - 1 + tstart;
 962                 if( wing == 2 ) newAnodef -= fNofAnodes;
 963                 Double_t anodePath = (newAnodef - fNofAnodes/2)*anodePitch;
 964                 newiTimef *= fTimeStep;
 965                 if( newiTimef > fTimeCorr ) newiTimef -= fTimeCorr;
 966                 if( electronics == 1 ){
 967                 //    newiTimef *= 0.999438;    // PASCAL
 968                 //    newiTimef += (6./fDriftSpeed - newiTimef/3000.);
 969                 }else if( electronics == 2 )
 970                     newiTimef *= 0.99714;    // OLA
 971
 972                 Int_t timeBin = (Int_t)(newiTimef/fTimeStep+0.5);
 973                 Int_t peakpos = Map()->GetHitIndex( newAnode, timeBin );
 974                 if( peakpos < 0 ) {
 975                     for( Int_t ii=0; ii<3; ii++ ) {
 976                         peakpos = Map()->GetHitIndex( newAnode, timeBin+ii );
 977                         if( peakpos > 0 ) break;
 978                         peakpos = Map()->GetHitIndex( newAnode, timeBin-ii );
 979                         if( peakpos > 0 ) break;
 980                     }
 981                 }
 982
 983                 if( peakpos < 0 ) {
 984                     //Warning("ResolveClusters",
 985                     //        "Digit not found for cluster");
 986                     //if(GetDebug(3)) clusterI.PrintInfo();
 987                    continue;
 988                 }
 989                 clusterI.SetPeakPos( peakpos );
 990                 Double_t driftPath = fSddLength - newiTimef * fDriftSpeed;
 991                 Double_t sign = ( wing == 1 ) ? -1. : 1.;
 992                 clusterI.SetX( driftPath*sign * 0.0001 );
 993                 clusterI.SetZ( anodePath * 0.0001 );
 994                 clusterI.SetAnode( newAnodef );
 995                 clusterI.SetTime( newiTimef );
 996                 clusterI.SetAsigma( sigma[i]*anodePitch );
 997                 clusterI.SetTsigma( tau[i]*fTimeStep );
 998                 clusterI.SetQ( integral[i] );
 999
1000                 fITS->AddCluster( 1, &clusterI );
1001             } // end for i
1002             Clusters()->RemoveAt( j );
1003             delete [] par;
1004         } else {  // something odd
1005             Warning( "ResolveClusters",
1006                      "--- Peak not found!!!!  minpeak=%d ,cluster peak= %f"
1007                      " , module= %d",
1008                      fMinPeak, clusterJ->PeakAmpl(),GetModule());
1009             clusterJ->PrintInfo();
1010             Warning( "ResolveClusters"," xdim= %d zdim= %d", xdim-2, zdim-2 );
1011         }
1012         delete [] sp;
1013     } // cluster loop
1014     Clusters()->Compress();
1015 //    Map()->ClearMap();
1016 }
1017 //________________________________________________________________________
1018 void  AliITSClusterFinderSDD::GroupClusters(){
1019     // group clusters
1020     Int_t dummy=0;
1021     Double_t fTimeStep = GetSeg()->Dpx(dummy);
1022     // get number of clusters for this module
1023     Int_t nofClusters = NClusters();
1024     nofClusters -= fNclusters;
1025     AliITSRawClusterSDD *clusterI;
1026     AliITSRawClusterSDD *clusterJ;
1027     Int_t *label = new Int_t [nofClusters];
1028     Int_t i,j;
1029     for(i=0; i<nofClusters; i++) label[i] = 0;
1030     for(i=0; i<nofClusters; i++) {
1031         if(label[i] != 0) continue;
1032         for(j=i+1; j<nofClusters; j++) {
1033             if(label[j] != 0) continue;
1034             clusterI = (AliITSRawClusterSDD*) Cluster(i);
1035             clusterJ = (AliITSRawClusterSDD*) Cluster(j);
1036             // 1.3 good
1037             if(clusterI->T() < fTimeStep*60) fDAnode = 4.2;  // TB 3.2
1038             if(clusterI->T() < fTimeStep*10) fDAnode = 1.5;  // TB 1.
1039             Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
1040             if(!pair) continue;
1041             if(GetDebug(4)){
1042                 clusterI->PrintInfo();
1043                 clusterJ->PrintInfo();
1044             } // end if GetDebug
1045             clusterI->Add(clusterJ);
1046             label[j] = 1;
1047             Clusters()->RemoveAt(j);
1048             j=i; // <- Ernesto
1049         } // J clusters
1050         label[i] = 1;
1051     } // I clusters
1052     Clusters()->Compress();
1053
1054     delete [] label;
1055     return;
1056 }
1057 //________________________________________________________________________
1058 void AliITSClusterFinderSDD::SelectClusters(){
1059     // get number of clusters for this module
1060     Int_t nofClusters = NClusters();
1061
1062     nofClusters -= fNclusters;
1063     Int_t i;
1064     for(i=0; i<nofClusters; i++) {
1065         AliITSRawClusterSDD *clusterI =(AliITSRawClusterSDD*) Cluster(i);
1066         Int_t rmflg = 0;
1067         Double_t wy = 0.;
1068         if(clusterI->Anodes() != 0.) {
1069             wy = ((Double_t) clusterI->Samples())/clusterI->Anodes();
1070         } // end if
1071         Int_t amp = (Int_t) clusterI->PeakAmpl();
1072         Int_t cha = (Int_t) clusterI->Q();
1073         if(amp < fMinPeak) rmflg = 1;
1074         if(cha < fMinCharge) rmflg = 1;
1075         if(wy < fMinNCells) rmflg = 1;
1076         //if(wy > fMaxNCells) rmflg = 1;
1077         if(rmflg) Clusters()->RemoveAt(i);
1078     } // I clusters
1079     Clusters()->Compress();
1080     return;
1081 }
1082
1083 //______________________________________________________________________
1084 void AliITSClusterFinderSDD::GetRecPoints(){
1085     // get rec points
1086
1087     // get number of clusters for this module
1088     Int_t nofClusters = NClusters();
1089     nofClusters -= fNclusters;
1090     const Double_t kconvGeV = 1.e-6; // GeV -> KeV
1091     const Double_t kconv = 1.0e-4;
1092     const Double_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
1093     const Double_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
1094     Int_t i;
1095     Int_t ix, iz, idx=-1;
1096     AliITSdigitSDD *dig=0;
1097     Int_t ndigits=NDigits();
1098     for(i=0; i<nofClusters; i++) {
1099         AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)Cluster(i);
1100         if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI);
1101         if(clusterI) idx=clusterI->PeakPos();
1102         if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits);
1103         // try peak neighbours - to be done
1104         if(idx&&idx<= ndigits) dig =(AliITSdigitSDD*)GetDigit(idx);
1105         if(!dig) {
1106             // try cog
1107             GetSeg()->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
1108             dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix-1);
1109             // if null try neighbours
1110             if (!dig) dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix);
1111             if (!dig) dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix+1);
1112             if (!dig) printf("SDD: cannot assign the track number!\n");
1113         } //  end if !dig
1114         AliITSRecPoint rnew;
1115         rnew.SetX(clusterI->X());
1116         rnew.SetZ(clusterI->Z());
1117         rnew.SetQ(clusterI->Q());   // in KeV - should be ADC
1118         rnew.SetdEdX(kconvGeV*clusterI->Q());
1119         rnew.SetSigmaX2(kRMSx*kRMSx);
1120         rnew.SetSigmaZ2(kRMSz*kRMSz);
1121
1122         if(dig) rnew.fTracks[0]=dig->GetTrack(0);
1123         if(dig) rnew.fTracks[1]=dig->GetTrack(1);
1124         if(dig) rnew.fTracks[2]=dig->GetTrack(2);
1125
1126         fITS->AddRecPoint(rnew);
1127     } // I clusters
1128 //    Map()->ClearMap();
1129 }
1130 //______________________________________________________________________
1131 void AliITSClusterFinderSDD::FindRawClusters(Int_t mod){
1132     // find raw clusters
1133
1134     SetModule(mod);
1135     Find1DClustersE();
1136     GroupClusters();
1137     SelectClusters();
1138     ResolveClusters();
1139     GetRecPoints();
1140 }
1141 //_______________________________________________________________________
1142 void AliITSClusterFinderSDD::PrintStatus() const{
1143     // Print SDD cluster finder Parameters
1144
1145     cout << "**************************************************" << endl;
1146     cout << " Silicon Drift Detector Cluster Finder Parameters " << endl;
1147     cout << "**************************************************" << endl;
1148     cout << "Number of Clusters: " << fNclusters << endl;
1149     cout << "Anode Tolerance: " << fDAnode << endl;
1150     cout << "Time  Tolerance: " << fDTime << endl;
1151     cout << "Time  correction (electronics): " << fTimeCorr << endl;
1152     cout << "Cut Amplitude (threshold): " << fCutAmplitude << endl;
1153     cout << "Minimum Amplitude: " << fMinPeak << endl;
1154     cout << "Minimum Charge: " << fMinCharge << endl;
1155     cout << "Minimum number of cells/clusters: " << fMinNCells << endl;
1156     cout << "Maximum number of cells/clusters: " << fMaxNCells << endl;
1157     cout << "**************************************************" << endl;
1158 }