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