[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSDD.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include <TFile.h>

#include "AliITSClusterFinderSDD.h"
#include "AliITSMapA1.h"
#include "AliITS.h"
#include "AliITSdigit.h"
#include "AliITSRawCluster.h"
#include "AliITSRecPoint.h"
#include "AliITSsegmentation.h"
#include "AliITSresponse.h"
#include "AliRun.h"


ClassImp(AliITSClusterFinderSDD)

//----------------------------------------------------------
AliITSClusterFinderSDD::AliITSClusterFinderSDD
(AliITSsegmentation *seg, AliITSresponse *response, TClonesArray *digits, TClonesArray *recp)   
{
  // constructor

    fSegmentation=seg;
    fResponse=response;
    fDigits=digits;
    fClusters=recp;
    fNclusters= fClusters->GetEntriesFast();
    SetCutAmplitude();
    SetDAnode();
    SetDTime();
    SetMinPeak();
    SetMinNCells();
    SetMaxNCells();
    SetTimeCorr();
    fMap=new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude);

}

//_____________________________________________________________________________
AliITSClusterFinderSDD::AliITSClusterFinderSDD()
{
  // constructor
    fSegmentation=0;
    fResponse=0;
    fDigits=0;
    fClusters=0;
    fNclusters=0;
    fMap=0;
    SetCutAmplitude();
    SetDAnode();
    SetDTime();
    SetMinPeak();
    SetMinNCells();
    SetMaxNCells();
    SetTimeCorr();

}

//_____________________________________________________________________________
AliITSClusterFinderSDD::~AliITSClusterFinderSDD()
{
    // destructor

    if(fMap) delete fMap;

}
//__________________________________________________________________________
AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){
  //     Copy Constructor 
  if(&source == this) return;
  this->fClusters = source.fClusters ;
  this->fNclusters = source.fNclusters ;
  this->fMap = source.fMap ;
  this->fCutAmplitude = source.fCutAmplitude ;
  this->fDAnode = source.fDAnode ;
  this->fDTime = source.fDTime ;
  this->fTimeCorr = source.fTimeCorr ;
  this->fMinPeak = source.fMinPeak ;
  this->fMinNCells = source.fMinNCells ;
  this->fMaxNCells = source.fMaxNCells ;
  return;
}

//_________________________________________________________________________
AliITSClusterFinderSDD& 
  AliITSClusterFinderSDD::operator=(const AliITSClusterFinderSDD &source) {
  //    Assignment operator
  if(&source == this) return *this;
  this->fClusters = source.fClusters ;
  this->fNclusters = source.fNclusters ;
  this->fMap = source.fMap ;
  this->fCutAmplitude = source.fCutAmplitude ;
  this->fDAnode = source.fDAnode ;
  this->fDTime = source.fDTime ;
  this->fTimeCorr = source.fTimeCorr ;
  this->fMinPeak = source.fMinPeak ;
  this->fMinNCells = source.fMinNCells ;
  this->fMaxNCells = source.fMaxNCells ;
  return *this;
}


//_____________________________________________________________________________

void AliITSClusterFinderSDD::Find1DClusters()
{
  // find 1D clusters

    AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");

   // retrieve the parameters 
   Int_t fNofMaps = fSegmentation->Npz();
   Int_t fMaxNofSamples = fSegmentation->Npx();
   Int_t fNofAnodes = fNofMaps/2;
   Int_t dummy=0;
   Float_t fTimeStep = fSegmentation->Dpx(dummy);
   Float_t fSddLength = fSegmentation->Dx();
   Float_t fDriftSpeed = fResponse->DriftSpeed();
   
   Float_t anodePitch = fSegmentation->Dpz(dummy);
   // map the signal
   fMap->SetThreshold(fCutAmplitude);
   fMap->FillMap();

   Float_t maxadc = fResponse->MaxAdc();    
   Float_t topValue = fResponse->MagicValue();
   Float_t norm = maxadc/topValue;

   Int_t nofFoundClusters = 0;
   Int_t i;
   Float_t **dfadc = new Float_t*[fNofAnodes];
   for(i=0;i<fNofAnodes;i++) dfadc[i] = new Float_t[fMaxNofSamples];
   Float_t fadc = 0.;
   Float_t fadc1 = 0.;
   Float_t fadc2 = 0.;
   Int_t j,k,idx,l,m;
   for(j=0;j<2;j++) {
     for(k=0;k<fNofAnodes;k++) {
       idx = j*fNofAnodes+k;
       // signal (fadc) & derivative (dfadc)
       dfadc[k][255]=0.;
       for(l=0; l<fMaxNofSamples; l++) {
	 fadc2=(Float_t)fMap->GetSignal(idx,l);
	 if(l>0) fadc1=(Float_t)fMap->GetSignal(idx,l-1);
	 if(l>0) dfadc[k][l-1] = fadc2-fadc1;
       } // samples
     } // anodes
    
     for(k=0;k<fNofAnodes;k++) {
       //cout << "Anode: " << k+1 << ", Wing: " << j+1 << endl;
       idx = j*fNofAnodes+k;
 
       Int_t imax = 0;
       Int_t imaxd = 0;
       Int_t it=0;
       while(it <= fMaxNofSamples-3) {
	
	 imax = it;
	 imaxd = it;
	 // maximum of signal	  
	
	 Float_t fadcmax = 0.;
	 Float_t dfadcmax = 0.;
	 Int_t lthrmina = 1;
	 //	if(it >= 60) lthrmina = 2;
	 //      if(it >= 100) lthrmina = 3;
	 Int_t lthrmint = 2;
	 //if(it >= 60) lthrmint = 3;
	 //if(it >= 100) lthrmint = 4;

	 Int_t lthra = 1;
	 Int_t lthrt = 0;
	
	 for(m=0;m<20;m++) {
	   Int_t id = it+m;
	   if(id>=fMaxNofSamples) break;
	   fadc=(float)fMap->GetSignal(idx,id);
	   if(fadc > fadcmax) { fadcmax = fadc; imax = id;}
	   if(fadc > (float)fCutAmplitude) { 
	     lthrt++; 
	   }

	   if(dfadc[k][id] > dfadcmax) {
	     dfadcmax = dfadc[k][id];
	     imaxd = id;
	   }
	 }
	 it = imaxd;
	
	 if(fMap->TestHit(idx,imax) == kEmpty) {it++; continue;}

	 // cluster charge
	 Int_t tstart = it-1;
	 if( tstart<0 ) tstart = 0;
	
	 Bool_t ilcl = 0;
	 if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
	 //printf("ilcl %d\n",ilcl);

	 if(ilcl) {
	   nofFoundClusters++;
	   Int_t tstop = tstart;
	   Float_t dfadcmin = 10000.;
	   Int_t ij;
	   for(ij=0; ij<20; ij++) {
	     if(dfadc[k][it+ij] < dfadcmin) {
	       tstop = it+ij+1;
	       dfadcmin = dfadc[k][it+ij];
	     }
	   }

	   Float_t clusterCharge = 0.;
	   Float_t clusterAnode = k+0.5;
	   Float_t clusterTime = 0.;
	   Float_t clusterMult = 0.;
	   Float_t clusterPeakAmplitude = 0.;
	   Int_t its,peakpos=-1;
	   Float_t n, baseline;
	   fResponse->GetNoiseParam(n,baseline);
	   n *= norm;
	   baseline *= norm;
	   for(its=tstart; its<=tstop; its++) {
	     fadc=(float)fMap->GetSignal(idx,its);
	     if(fadc>baseline)
	       fadc-=baseline;
	     else
	       fadc=0.;
	     clusterCharge += fadc;
	     // as a matter of fact we should take the peak pos before FFT
	     // to get the list of tracks !!!
	     if(fadc > clusterPeakAmplitude) {
	       clusterPeakAmplitude = fadc;
	       //peakpos=fMap->GetHitIndex(idx,its);
	       Int_t shift=(int)(fTimeCorr/fTimeStep);
	       if(its>shift && its<(fMaxNofSamples-shift)) peakpos=fMap->GetHitIndex(idx,its+shift);
	       else peakpos=fMap->GetHitIndex(idx,its);
	       if(peakpos<0) peakpos=fMap->GetHitIndex(idx,its);
	     }
	     clusterTime += fadc*its;
	     clusterMult++;
	     if(its == tstop) {
	       // charge from ADC back to nA 
	       //clusterCharge /= norm;
	       if(clusterCharge <= 0.) printf("clusterCharge %f norm %f\n",clusterCharge,norm);
	       clusterTime /= (clusterCharge/fTimeStep);   // ns
	       clusterCharge *= (fTimeStep/160.);          // keV
	       if(clusterTime > fTimeCorr) clusterTime -= fTimeCorr;   // ns
	     }
	   }
	   // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;

	   Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)*anodePitch;
	   Float_t clusterDriftPath = clusterTime*fDriftSpeed;
	   if(TMath::Abs(clusterDriftPath) > fSddLength) {
	     Warning("AliITSClusterFinderSDD","Cluster drift path %f bigger then the detector size - please parametrise the time correction as a function of the drift time!",clusterDriftPath);
	   }
	   clusterDriftPath = fSddLength-clusterDriftPath;

	   if(clusterCharge <= 0.) break;

	   AliITSRawClusterSDD clust(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,peakpos,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
	   iTS->AddCluster(1,&clust);
	   it = tstop;
	} // ilcl
	
	it++;
	
      } // while (samples)
    } // anodes
  } // detectors (2)


   //fMap->ClearMap(); 
  
  for(i=0;i<fNofAnodes;i++) delete[] dfadc[i];
  delete [] dfadc;

  return;

}

//_____________________________________________________________________________
void  AliITSClusterFinderSDD::GroupClusters()
{
  // group clusters
  Int_t dummy=0;
  Float_t fTimeStep = fSegmentation->Dpx(dummy);


  // get number of clusters for this module
  Int_t nofClusters = fClusters->GetEntriesFast();
  nofClusters -= fNclusters;

  AliITSRawClusterSDD *clusterI;
  AliITSRawClusterSDD *clusterJ;

  Int_t *label = new Int_t [nofClusters];
  Int_t i,j;
  for(i=0; i<nofClusters; i++) label[i] = 0;
  for(i=0; i<nofClusters; i++) { 
    if(label[i] != 0) continue;
    for(j=i+1; j<nofClusters; j++) { 
      if(label[j] != 0) continue;
      clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
      clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
      // 1.3 good
      if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
      if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
      Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
      if(!pair) continue;
      //      clusterI->PrintInfo();
      //      clusterJ->PrintInfo();
      clusterI->Add(clusterJ);
      label[j] = 1;
      fClusters->RemoveAt(j);
    } // J clusters  
    label[i] = 1;
  } // I clusters
  fClusters->Compress();
  
  delete [] label;
  return;

}

//_____________________________________________________________________________

void AliITSClusterFinderSDD::SelectClusters()
{
  // get number of clusters for this module
  Int_t nofClusters = fClusters->GetEntriesFast();
  nofClusters -= fNclusters;

  Int_t i;
  for(i=0; i<nofClusters; i++) { 
    AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
    Int_t rmflg = 0;
    Float_t wy = 0.;
    if(clusterI->Anodes() != 0.) {
      wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
    }
    Float_t amp = clusterI->PeakAmpl();
    if(amp < fMinPeak) rmflg = 1;  
    if(wy < fMinNCells) rmflg = 1;
    //if(wy > fMaxNCells) rmflg = 1;
    if(rmflg) fClusters->RemoveAt(i);
  } // I clusters
  fClusters->Compress();
  return;

}

//_____________________________________________________________________________

void AliITSClusterFinderSDD::GetRecPoints()
{
  // get rec points

  AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");

  // get number of clusters for this module
  Int_t nofClusters = fClusters->GetEntriesFast();
  nofClusters -= fNclusters;

  const Float_t kconvGeV = 1.e-6; // GeV -> KeV
  const Float_t kconv = 1.0e-4; 
  const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
  const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3


  Int_t i;
  Int_t ix, iz, idx=-1;
  AliITSdigitSDD *dig=0;
  //Int_t maxt=fSegmentation->Npx();
  Int_t ndigits=fDigits->GetEntriesFast();
  for(i=0; i<nofClusters; i++) { 
    AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)fClusters->At(i);
    if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI);
    if(clusterI) idx=clusterI->PeakPos();
    if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits);
    // try peak neighbours - to be done 
    if(idx && idx <= ndigits) dig = (AliITSdigitSDD*)fDigits->UncheckedAt(idx);
    if(!dig) {
	// try cog
	fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
	dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
	// if null try neighbours
	if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix); 
	if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix+1); 
        if (!dig) printf("SDD: cannot assign the track number!\n");
    }

    AliITSRecPoint rnew;
    rnew.SetX(clusterI->X());
    rnew.SetZ(clusterI->Z());
    rnew.SetQ(clusterI->Q());   // in KeV - should be ADC
    rnew.SetdEdX(kconvGeV*clusterI->Q());
    rnew.SetSigmaX2(kRMSx*kRMSx);
    rnew.SetSigmaZ2(kRMSz*kRMSz);
    if(dig) rnew.fTracks[0]=dig->fTracks[0];
    if(dig) rnew.fTracks[1]=dig->fTracks[1];
    if(dig) rnew.fTracks[2]=dig->fTracks[2];
    //printf("SDD: i %d track1 track2 track3 %d %d %d x y %f %f\n",i,rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],clusterI->X(),clusterI->Z());
    iTS->AddRecPoint(rnew);
  } // I clusters

  fMap->ClearMap();
}

//_____________________________________________________________________________

void AliITSClusterFinderSDD::FindRawClusters()
{
  // find raw clusters
    Find1DClusters();
    GroupClusters();
    SelectClusters();
    GetRecPoints();
}
Commit	Line	Data
b0f5e3fc	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
78a228db	16	#include <TFile.h>
b0f5e3fc	17
b0f5e3fc	18	#include "AliITSClusterFinderSDD.h"
e8189707	19	#include "AliITSMapA1.h"
e8189707	20	#include "AliITS.h"
78a228db	21	#include "AliITSdigit.h"
	22	#include "AliITSRawCluster.h"
	23	#include "AliITSRecPoint.h"
	24	#include "AliITSsegmentation.h"
	25	#include "AliITSresponse.h"
b0f5e3fc	26	#include "AliRun.h"
	27
	28
	29
	30	ClassImp(AliITSClusterFinderSDD)
	31
	32	//----------------------------------------------------------
	33	AliITSClusterFinderSDD::AliITSClusterFinderSDD
	34	(AliITSsegmentation seg, AliITSresponse response, TClonesArray digits, TClonesArray recp)
	35	{
	36	// constructor
78a228db	37
b0f5e3fc	38	fSegmentation=seg;
	39	fResponse=response;
	40	fDigits=digits;
	41	fClusters=recp;
	42	fNclusters= fClusters->GetEntriesFast();
b0f5e3fc	43	SetCutAmplitude();
	44	SetDAnode();
	45	SetDTime();
b0f5e3fc	46	SetMinPeak();
78a228db	47	SetMinNCells();
	48	SetMaxNCells();
	49	SetTimeCorr();
	50	fMap=new AliITSMapA1(fSegmentation,fDigits,fCutAmplitude);
	51
b0f5e3fc	52	}
	53
	54	//_____________________________________________________________________________
	55	AliITSClusterFinderSDD::AliITSClusterFinderSDD()
	56	{
	57	// constructor
	58	fSegmentation=0;
	59	fResponse=0;
	60	fDigits=0;
	61	fClusters=0;
	62	fNclusters=0;
e8189707	63	fMap=0;
b0f5e3fc	64	SetCutAmplitude();
	65	SetDAnode();
	66	SetDTime();
b0f5e3fc	67	SetMinPeak();
78a228db	68	SetMinNCells();
	69	SetMaxNCells();
	70	SetTimeCorr();
b0f5e3fc	71
	72	}
	73
e8189707	74	//_____________________________________________________________________________
	75	AliITSClusterFinderSDD::~AliITSClusterFinderSDD()
	76	{
	77	// destructor
	78
	79	if(fMap) delete fMap;
e8189707	80
e8189707	81	}
b0f5e3fc	82	//__________________________________________________________________________
	83	AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){
	84	// Copy Constructor
	85	if(&source == this) return;
	86	this->fClusters = source.fClusters ;
	87	this->fNclusters = source.fNclusters ;
	88	this->fMap = source.fMap ;
	89	this->fCutAmplitude = source.fCutAmplitude ;
	90	this->fDAnode = source.fDAnode ;
	91	this->fDTime = source.fDTime ;
78a228db	92	this->fTimeCorr = source.fTimeCorr ;
b0f5e3fc	93	this->fMinPeak = source.fMinPeak ;
b0f5e3fc	94	this->fMinNCells = source.fMinNCells ;
78a228db	95	this->fMaxNCells = source.fMaxNCells ;
b0f5e3fc	96	return;
	97	}
	98
	99	//_________________________________________________________________________
	100	AliITSClusterFinderSDD&
	101	AliITSClusterFinderSDD::operator=(const AliITSClusterFinderSDD &source) {
	102	// Assignment operator
	103	if(&source == this) return *this;
	104	this->fClusters = source.fClusters ;
	105	this->fNclusters = source.fNclusters ;
	106	this->fMap = source.fMap ;
	107	this->fCutAmplitude = source.fCutAmplitude ;
	108	this->fDAnode = source.fDAnode ;
	109	this->fDTime = source.fDTime ;
78a228db	110	this->fTimeCorr = source.fTimeCorr ;
b0f5e3fc	111	this->fMinPeak = source.fMinPeak ;
b0f5e3fc	112	this->fMinNCells = source.fMinNCells ;
78a228db	113	this->fMaxNCells = source.fMaxNCells ;
b0f5e3fc	114	return *this;
	115	}
	116
b0f5e3fc	117
b0f5e3fc	118	//_____________________________________________________________________________
	119
	120	void AliITSClusterFinderSDD::Find1DClusters()
	121	{
	122	// find 1D clusters
	123
	124	AliITS iTS=(AliITS)gAlice->GetModule("ITS");
	125
	126	// retrieve the parameters
	127	Int_t fNofMaps = fSegmentation->Npz();
	128	Int_t fMaxNofSamples = fSegmentation->Npx();
	129	Int_t fNofAnodes = fNofMaps/2;
	130	Int_t dummy=0;
	131	Float_t fTimeStep = fSegmentation->Dpx(dummy);
	132	Float_t fSddLength = fSegmentation->Dx();
	133	Float_t fDriftSpeed = fResponse->DriftSpeed();
	134
	135	Float_t anodePitch = fSegmentation->Dpz(dummy);
	136	// map the signal
78a228db	137	fMap->SetThreshold(fCutAmplitude);
	138	fMap->FillMap();
	139
78a228db	140	Float_t maxadc = fResponse->MaxAdc();
	141	Float_t topValue = fResponse->MagicValue();
	142	Float_t norm = maxadc/topValue;
	143
	144	Int_t nofFoundClusters = 0;
	145	Int_t i;
	146	Float_t *dfadc = new Float_t[fNofAnodes];
	147	for(i=0;i<fNofAnodes;i++) dfadc[i] = new Float_t[fMaxNofSamples];
063150f2	148	Float_t fadc = 0.;
	149	Float_t fadc1 = 0.;
	150	Float_t fadc2 = 0.;
78a228db	151	Int_t j,k,idx,l,m;
	152	for(j=0;j<2;j++) {
	153	for(k=0;k<fNofAnodes;k++) {
	154	idx = j*fNofAnodes+k;
	155	// signal (fadc) & derivative (dfadc)
65d4384f	156	dfadc[k][255]=0.;
78a228db	157	for(l=0; l<fMaxNofSamples; l++) {
	158	fadc2=(Float_t)fMap->GetSignal(idx,l);
	159	if(l>0) fadc1=(Float_t)fMap->GetSignal(idx,l-1);
	160	if(l>0) dfadc[k][l-1] = fadc2-fadc1;
	161	} // samples
	162	} // anodes
b0f5e3fc	163
78a228db	164	for(k=0;k<fNofAnodes;k++) {
	165	//cout << "Anode: " << k+1 << ", Wing: " << j+1 << endl;
	166	idx = j*fNofAnodes+k;
b0f5e3fc	167
78a228db	168	Int_t imax = 0;
	169	Int_t imaxd = 0;
	170	Int_t it=0;
	171	while(it <= fMaxNofSamples-3) {
b0f5e3fc	172
78a228db	173	imax = it;
	174	imaxd = it;
	175	// maximum of signal
b0f5e3fc	176
78a228db	177	Float_t fadcmax = 0.;
	178	Float_t dfadcmax = 0.;
	179	Int_t lthrmina = 1;
	180	// if(it >= 60) lthrmina = 2;
	181	// if(it >= 100) lthrmina = 3;
	182	Int_t lthrmint = 2;
	183	//if(it >= 60) lthrmint = 3;
	184	//if(it >= 100) lthrmint = 4;
	185
	186	Int_t lthra = 1;
	187	Int_t lthrt = 0;
b0f5e3fc	188
78a228db	189	for(m=0;m<20;m++) {
	190	Int_t id = it+m;
	191	if(id>=fMaxNofSamples) break;
	192	fadc=(float)fMap->GetSignal(idx,id);
	193	if(fadc > fadcmax) { fadcmax = fadc; imax = id;}
	194	if(fadc > (float)fCutAmplitude) {
	195	lthrt++;
	196	}
	197
	198	if(dfadc[k][id] > dfadcmax) {
	199	dfadcmax = dfadc[k][id];
	200	imaxd = id;
	201	}
	202	}
	203	it = imaxd;
b0f5e3fc	204
78a228db	205	if(fMap->TestHit(idx,imax) == kEmpty) {it++; continue;}
	206
	207	// cluster charge
	208	Int_t tstart = it-1;
e3d05a9c	209	if( tstart<0 ) tstart = 0;
b0f5e3fc	210
78a228db	211	Bool_t ilcl = 0;
	212	if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
	213	//printf("ilcl %d\n",ilcl);
	214
78a228db	215	if(ilcl) {
	216	nofFoundClusters++;
	217	Int_t tstop = tstart;
	218	Float_t dfadcmin = 10000.;
	219	Int_t ij;
	220	for(ij=0; ij<20; ij++) {
	221	if(dfadc[k][it+ij] < dfadcmin) {
	222	tstop = it+ij+1;
	223	dfadcmin = dfadc[k][it+ij];
	224	}
	225	}
	226
	227	Float_t clusterCharge = 0.;
	228	Float_t clusterAnode = k+0.5;
	229	Float_t clusterTime = 0.;
	230	Float_t clusterMult = 0.;
	231	Float_t clusterPeakAmplitude = 0.;
	232	Int_t its,peakpos=-1;
	233	Float_t n, baseline;
	234	fResponse->GetNoiseParam(n,baseline);
210bc185	235	n *= norm;
210bc185	236	baseline *= norm;
78a228db	237	for(its=tstart; its<=tstop; its++) {
	238	fadc=(float)fMap->GetSignal(idx,its);
	239	if(fadc>baseline)
	240	fadc-=baseline;
	241	else
	242	fadc=0.;
	243	clusterCharge += fadc;
	244	// as a matter of fact we should take the peak pos before FFT
	245	// to get the list of tracks !!!
	246	if(fadc > clusterPeakAmplitude) {
	247	clusterPeakAmplitude = fadc;
	248	//peakpos=fMap->GetHitIndex(idx,its);
	249	Int_t shift=(int)(fTimeCorr/fTimeStep);
	250	if(its>shift && its<(fMaxNofSamples-shift)) peakpos=fMap->GetHitIndex(idx,its+shift);
	251	else peakpos=fMap->GetHitIndex(idx,its);
	252	if(peakpos<0) peakpos=fMap->GetHitIndex(idx,its);
	253	}
	254	clusterTime += fadc*its;
	255	clusterMult++;
	256	if(its == tstop) {
	257	// charge from ADC back to nA
063150f2	258	//clusterCharge /= norm;
78a228db	259	if(clusterCharge <= 0.) printf("clusterCharge %f norm %f\n",clusterCharge,norm);
b0f5e3fc	260	clusterTime /= (clusterCharge/fTimeStep); // ns
b0f5e3fc	261	clusterCharge *= (fTimeStep/160.); // keV
78a228db	262	if(clusterTime > fTimeCorr) clusterTime -= fTimeCorr; // ns
	263	}
	264	}
	265	// cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
	266
	267	Float_t clusteranodePath = (clusterAnode - fNofAnodes/2)*anodePitch;
	268	Float_t clusterDriftPath = clusterTime*fDriftSpeed;
063150f2	269	if(TMath::Abs(clusterDriftPath) > fSddLength) {
	270	Warning("AliITSClusterFinderSDD","Cluster drift path %f bigger then the detector size - please parametrise the time correction as a function of the drift time!",clusterDriftPath);
	271	}
78a228db	272	clusterDriftPath = fSddLength-clusterDriftPath;
	273
	274	if(clusterCharge <= 0.) break;
	275
210bc185	276	AliITSRawClusterSDD clust(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,peakpos,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
210bc185	277	iTS->AddCluster(1,&clust);
78a228db	278	it = tstop;
b0f5e3fc	279	} // ilcl
	280
	281	it++;
	282
	283	} // while (samples)
	284	} // anodes
	285	} // detectors (2)
	286
b0f5e3fc	287
78a228db	288	//fMap->ClearMap();
44b3710f	289
78a228db	290	for(i=0;i<fNofAnodes;i++) delete[] dfadc[i];
44b3710f	291	delete [] dfadc;
b0f5e3fc	292
	293	return;
	294
	295	}
	296
	297	//_____________________________________________________________________________
	298	void AliITSClusterFinderSDD::GroupClusters()
	299	{
	300	// group clusters
	301	Int_t dummy=0;
	302	Float_t fTimeStep = fSegmentation->Dpx(dummy);
	303
	304
	305	// get number of clusters for this module
	306	Int_t nofClusters = fClusters->GetEntriesFast();
	307	nofClusters -= fNclusters;
	308
b0f5e3fc	309	AliITSRawClusterSDD *clusterI;
	310	AliITSRawClusterSDD *clusterJ;
	311
44b3710f	312	Int_t *label = new Int_t [nofClusters];
b0f5e3fc	313	Int_t i,j;
	314	for(i=0; i<nofClusters; i++) label[i] = 0;
	315	for(i=0; i<nofClusters; i++) {
	316	if(label[i] != 0) continue;
	317	for(j=i+1; j<nofClusters; j++) {
	318	if(label[j] != 0) continue;
	319	clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
	320	clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
	321	// 1.3 good
	322	if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
	323	if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
	324	Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
	325	if(!pair) continue;
e8189707	326	// clusterI->PrintInfo();
e8189707	327	// clusterJ->PrintInfo();
b0f5e3fc	328	clusterI->Add(clusterJ);
	329	label[j] = 1;
	330	fClusters->RemoveAt(j);
	331	} // J clusters
	332	label[i] = 1;
	333	} // I clusters
	334	fClusters->Compress();
44b3710f	335
44b3710f	336	delete [] label;
b0f5e3fc	337	return;
	338
	339	}
	340
	341	//_____________________________________________________________________________
	342
	343	void AliITSClusterFinderSDD::SelectClusters()
	344	{
	345	// get number of clusters for this module
	346	Int_t nofClusters = fClusters->GetEntriesFast();
	347	nofClusters -= fNclusters;
	348
b0f5e3fc	349	Int_t i;
	350	for(i=0; i<nofClusters; i++) {
	351	AliITSRawClusterSDD clusterI = (AliITSRawClusterSDD) fClusters->At(i);
	352	Int_t rmflg = 0;
	353	Float_t wy = 0.;
	354	if(clusterI->Anodes() != 0.) {
	355	wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
	356	}
	357	Float_t amp = clusterI->PeakAmpl();
	358	if(amp < fMinPeak) rmflg = 1;
	359	if(wy < fMinNCells) rmflg = 1;
78a228db	360	//if(wy > fMaxNCells) rmflg = 1;
b0f5e3fc	361	if(rmflg) fClusters->RemoveAt(i);
	362	} // I clusters
	363	fClusters->Compress();
	364	return;
	365
	366	}
	367
	368	//_____________________________________________________________________________
	369
	370	void AliITSClusterFinderSDD::GetRecPoints()
	371	{
	372	// get rec points
b0f5e3fc	373
	374	AliITS iTS=(AliITS)gAlice->GetModule("ITS");
	375
	376	// get number of clusters for this module
	377	Int_t nofClusters = fClusters->GetEntriesFast();
	378	nofClusters -= fNclusters;
	379
b0f5e3fc	380	const Float_t kconvGeV = 1.e-6; // GeV -> KeV
	381	const Float_t kconv = 1.0e-4;
	382	const Float_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
	383	const Float_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
	384
e8189707	385
78a228db	386	Int_t i;
	387	Int_t ix, iz, idx=-1;
	388	AliITSdigitSDD *dig=0;
063150f2	389	//Int_t maxt=fSegmentation->Npx();
78a228db	390	Int_t ndigits=fDigits->GetEntriesFast();
b0f5e3fc	391	for(i=0; i<nofClusters; i++) {
78a228db	392	AliITSRawClusterSDD clusterI = (AliITSRawClusterSDD)fClusters->At(i);
	393	if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI);
	394	if(clusterI) idx=clusterI->PeakPos();
	395	if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits);
	396	// try peak neighbours - to be done
	397	if(idx && idx <= ndigits) dig = (AliITSdigitSDD*)fDigits->UncheckedAt(idx);
	398	if(!dig) {
	399	// try cog
	400	fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
	401	dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix-1);
	402	// if null try neighbours
	403	if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix);
	404	if (!dig) dig = (AliITSdigitSDD*)fMap->GetHit(iz-1,ix+1);
	405	if (!dig) printf("SDD: cannot assign the track number!\n");
	406	}
	407
b0f5e3fc	408	AliITSRecPoint rnew;
	409	rnew.SetX(clusterI->X());
	410	rnew.SetZ(clusterI->Z());
	411	rnew.SetQ(clusterI->Q()); // in KeV - should be ADC
b0f5e3fc	412	rnew.SetdEdX(kconvGeV*clusterI->Q());
	413	rnew.SetSigmaX2(kRMSx*kRMSx);
	414	rnew.SetSigmaZ2(kRMSz*kRMSz);
78a228db	415	if(dig) rnew.fTracks[0]=dig->fTracks[0];
	416	if(dig) rnew.fTracks[1]=dig->fTracks[1];
	417	if(dig) rnew.fTracks[2]=dig->fTracks[2];
	418	//printf("SDD: i %d track1 track2 track3 %d %d %d x y %f %f\n",i,rnew.fTracks[0],rnew.fTracks[1],rnew.fTracks[2],clusterI->X(),clusterI->Z());
b0f5e3fc	419	iTS->AddRecPoint(rnew);
b0f5e3fc	420	} // I clusters
b0f5e3fc	421
e8189707	422	fMap->ClearMap();
b0f5e3fc	423	}
	424
	425	//_____________________________________________________________________________
	426
	427	void AliITSClusterFinderSDD::FindRawClusters()
	428	{
	429	// find raw clusters
	430	Find1DClusters();
	431	GroupClusters();
	432	SelectClusters();
	433	GetRecPoints();
	434	}