[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 "AliITSClusterFinderSDD.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();
    printf("SDD: fNclusters %d\n",fNclusters);
    SetCutAmplitude();
    SetDAnode();
    SetDTime();
    SetMap();
    SetMinPeak();
    SetNCells();
}

//_____________________________________________________________________________
AliITSClusterFinderSDD::AliITSClusterFinderSDD()
{
  // constructor
    fSegmentation=0;
    fResponse=0;
    fDigits=0;
    fClusters=0;
    fNclusters=0;
    SetCutAmplitude();
    SetDAnode();
    SetDTime();
    SetMap();
    SetMinPeak();
    SetNCells();

}

//__________________________________________________________________________
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->fMinPeak = source.fMinPeak ;
  this->fMinNCells = source.fMinNCells ;
  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->fMinPeak = source.fMinPeak ;
  this->fMinNCells = source.fMinNCells ;
  return *this;
}

//_____________________________________________________________________________
void AliITSClusterFinderSDD::SetMap()
{
  // set map
    if(!fMap) fMap=new AliITSMapA2(fSegmentation);

}
//_____________________________________________________________________________
void AliITSClusterFinderSDD::FillMap()
{
  // fCoord1 = anode #
  // fCoord2 = time sample

  if (!fDigits) return; 
 
  Int_t ndigits = fDigits->GetEntriesFast();
  //printf("FillMap: ndigits %d\n",ndigits);
  if (!ndigits) return;

    AliITSdigitSDD *dig;
    Int_t ndig;
    for(ndig=0; ndig<ndigits; ndig++) {
        dig = (AliITSdigitSDD*)fDigits->UncheckedAt(ndig);
        Double_t signal=dig->fSignal;
	//printf("FillMap: ndig fCoord1 fCoord2 signal %d %d %d %f\n",ndig,dig->fCoord1,dig->fCoord2,signal);
        fMap->SetHit(dig->fCoord1,dig->fCoord2,signal);
    }


}
//_____________________________________________________________________________

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
   FillMap();

   // Piergiorgio's code - do not subtract baseline since we start
   // from digits and do not duplicate arrays, i.e. use fMap instead
   // of Float_t fadc[2*fNofAnodes][fMaxNofSamples];
 

  Int_t nofFoundClusters = 0;
  Float_t dfadc[fNofMaps][fMaxNofSamples];
  Float_t fadc, fadc1, fadc2;
  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)
      for(l=0; l<fMaxNofSamples; l++) {
	fadc2=(Float_t)fMap->GetSignal(idx,l);
	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++) {
      idx = j*fNofAnodes+k;
 
      Int_t imax = 0;
      Int_t imaxd = 0;
      Int_t it=0;
      while(it <= fMaxNofSamples-3) {
	//	cout << "sample: " << it << endl;
	
	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 = 0;
	Int_t lthrt = 0;
	for(m=0;m<10;m++) {
	  Int_t id = it+m;
	  if(id>=fMaxNofSamples) break;
          fadc=fMap->GetSignal(idx,id);
	  if(fadc > fadcmax) {
	    fadcmax = fadc;
	    if(fadc > fCutAmplitude) { lthra++; lthrt++; }
	    imax = id;
	  }
	  if(dfadc[k][id] > dfadcmax) {
	    dfadcmax = dfadc[k][id];
	    imaxd = id;
	  }
	}
	it = imaxd;
	// skip if no signal over threshold
	if(fMap->GetSignal(idx,imax) < fCutAmplitude) {it++; continue;}

	if(k>0) {
	  if(fMap->GetSignal(idx-1,imax) > fCutAmplitude) lthra++;
	}
	if(k<fNofAnodes-1)
	  if(fMap->GetSignal(idx+1,imax) > fCutAmplitude) lthra++;

	if(imax>0) {
	  if(fMap->GetSignal(idx,imax-1) > fCutAmplitude) lthrt++;
	}
	if(imax<fMaxNofSamples)
	  if(fMap->GetSignal(idx,imax+1) > fCutAmplitude) lthrt++;
	
	// cluster charge
	Int_t tstart = it-1;
	
	Bool_t ilcl = 0;
	if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
	
	if(ilcl) {
	  nofFoundClusters++;
	  Int_t tstop = tstart;
	  Float_t dfadcmin = 10000.;
          Int_t ij;
	  for(ij=0; ij<10; 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;
	  for(its=tstart; its<=tstop; its++) {
            fadc=fMap->GetSignal(idx,its);
	    clusterCharge += fadc;
	    if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc;
	    clusterTime += fadc*its;
	    clusterMult++;
	    if(its == tstop) {
	       clusterTime /= (clusterCharge/fTimeStep);   // ns
	       clusterCharge *= (fTimeStep/160.);          // keV
	       if(clusterTime > 58.2) clusterTime -= 58.2;   // ns
	      /*
	      else {
		cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl;
		cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
	      }
	      */
	      }
	  }
	  // cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;

	  Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch;
	  Float_t clusterDriftPath = clusterTime*fDriftSpeed;
	  clusterDriftPath = fSddLength-clusterDriftPath;

          if(clusterCharge < 0.) break;

	  //printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime);

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

  Int_t nofClusters = fClusters->GetEntriesFast();
  nofClusters -= fNclusters;

  //printf("SDD- Find1Dclust: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);

  fMap->ClearMap();

  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;

  //printf("SDD- GroupClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);

  AliITSRawClusterSDD *clusterI;
  AliITSRawClusterSDD *clusterJ;

  Int_t label[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->Print();
      //      clusterJ->Print();
      clusterI->Add(clusterJ);
      label[j] = 1;
      fClusters->RemoveAt(j);
    } // J clusters  
    label[i] = 1;
  } // I clusters
  fClusters->Compress();
  return;

}

//_____________________________________________________________________________

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

  //printf("SDD- SelectClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);

  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(rmflg) fClusters->RemoveAt(i);
  } // I clusters
  fClusters->Compress();
  return;

}

//_____________________________________________________________________________

void AliITSClusterFinderSDD::GetRecPoints()
{
  // get rec points
  //static Int_t counter=0;

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

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

  //printf("SDD- GetRecPoints: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);

//  const Float_t kdEdXtoQ = 2.778e+2; // KeV -> number of e-hole pairs in Si
  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;
  for(i=0; i<nofClusters; i++) { 
    AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
    AliITSRecPoint rnew;
    rnew.SetX(clusterI->X());
    rnew.SetZ(clusterI->Z());
    rnew.SetQ(clusterI->Q());   // in KeV - should be ADC
    //rnew.SetdEdX((clusterI->Q())/kdEdXtoQ);
    rnew.SetdEdX(kconvGeV*clusterI->Q());
    rnew.SetSigmaX2(kRMSx*kRMSx);
    rnew.SetSigmaZ2(kRMSz*kRMSz);
    rnew.SetProbability(1.);
    iTS->AddRecPoint(rnew);
   //counter++;
  } // I clusters
  //printf("counter %d\n",counter);

}

//_____________________________________________________________________________

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
	16
	17	#include "AliITSClusterFinderSDD.h"
	18	#include "AliRun.h"
	19
	20
	21
	22	ClassImp(AliITSClusterFinderSDD)
	23
	24	//----------------------------------------------------------
	25	AliITSClusterFinderSDD::AliITSClusterFinderSDD
	26	(AliITSsegmentation seg, AliITSresponse response, TClonesArray digits, TClonesArray recp)
	27	{
	28	// constructor
	29	fSegmentation=seg;
	30	fResponse=response;
	31	fDigits=digits;
	32	fClusters=recp;
	33	fNclusters= fClusters->GetEntriesFast();
	34	printf("SDD: fNclusters %d\n",fNclusters);
	35	SetCutAmplitude();
	36	SetDAnode();
	37	SetDTime();
	38	SetMap();
	39	SetMinPeak();
	40	SetNCells();
	41	}
	42
	43	//_____________________________________________________________________________
	44	AliITSClusterFinderSDD::AliITSClusterFinderSDD()
	45	{
	46	// constructor
	47	fSegmentation=0;
	48	fResponse=0;
	49	fDigits=0;
	50	fClusters=0;
	51	fNclusters=0;
	52	SetCutAmplitude();
	53	SetDAnode();
	54	SetDTime();
	55	SetMap();
	56	SetMinPeak();
	57	SetNCells();
	58
	59	}
	60
	61	//__________________________________________________________________________
	62	AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){
	63	// Copy Constructor
	64	if(&source == this) return;
65	this->fClusters = source.fClusters ;
66	this->fNclusters = source.fNclusters ;
67	this->fMap = source.fMap ;
68	this->fCutAmplitude = source.fCutAmplitude ;
69	this->fDAnode = source.fDAnode ;
70	this->fDTime = source.fDTime ;
71	this->fMinPeak = source.fMinPeak ;
72	this->fMinNCells = source.fMinNCells ;
73	return;
74	}
75
76	//_________________________________________________________________________
77	AliITSClusterFinderSDD&
78	AliITSClusterFinderSDD::operator=(const AliITSClusterFinderSDD &source) {
79	// Assignment operator
80	if(&source == this) return *this;
81	this->fClusters = source.fClusters ;
82	this->fNclusters = source.fNclusters ;
83	this->fMap = source.fMap ;
84	this->fCutAmplitude = source.fCutAmplitude ;
85	this->fDAnode = source.fDAnode ;
86	this->fDTime = source.fDTime ;
87	this->fMinPeak = source.fMinPeak ;
88	this->fMinNCells = source.fMinNCells ;
89	return *this;
90	}
91
92	//_____________________________________________________________________________
93	void AliITSClusterFinderSDD::SetMap()
94	{
95	// set map
96	if(!fMap) fMap=new AliITSMapA2(fSegmentation);
97
98	}
99	//_____________________________________________________________________________
100	void AliITSClusterFinderSDD::FillMap()
101	{
102	// fCoord1 = anode #
103	// fCoord2 = time sample
104
105	if (!fDigits) return;
106
107	Int_t ndigits = fDigits->GetEntriesFast();
108	//printf("FillMap: ndigits %d\n",ndigits);
109	if (!ndigits) return;
110
111	AliITSdigitSDD *dig;
112	Int_t ndig;
113	for(ndig=0; ndig<ndigits; ndig++) {
114	dig = (AliITSdigitSDD*)fDigits->UncheckedAt(ndig);
115	Double_t signal=dig->fSignal;
116	//printf("FillMap: ndig fCoord1 fCoord2 signal %d %d %d %f\n",ndig,dig->fCoord1,dig->fCoord2,signal);
117	fMap->SetHit(dig->fCoord1,dig->fCoord2,signal);
118	}
119
120
121	}
122	//_____________________________________________________________________________
123
124	void AliITSClusterFinderSDD::Find1DClusters()
125	{
126	// find 1D clusters
127
128	AliITS iTS=(AliITS)gAlice->GetModule("ITS");
129
130	// retrieve the parameters
131	Int_t fNofMaps = fSegmentation->Npz();
132	Int_t fMaxNofSamples = fSegmentation->Npx();
133	Int_t fNofAnodes = fNofMaps/2;
134	Int_t dummy=0;
135	Float_t fTimeStep = fSegmentation->Dpx(dummy);
136	Float_t fSddLength = fSegmentation->Dx();
137	Float_t fDriftSpeed = fResponse->DriftSpeed();
138
139	Float_t anodePitch = fSegmentation->Dpz(dummy);
140	// map the signal
141	FillMap();
142
143	// Piergiorgio's code - do not subtract baseline since we start
144	// from digits and do not duplicate arrays, i.e. use fMap instead
145	// of Float_t fadc[2*fNofAnodes][fMaxNofSamples];
146
147
148	Int_t nofFoundClusters = 0;
149	Float_t dfadc[fNofMaps][fMaxNofSamples];
150	Float_t fadc, fadc1, fadc2;
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)
156	for(l=0; l<fMaxNofSamples; l++) {
157	fadc2=(Float_t)fMap->GetSignal(idx,l);
158	fadc1=(Float_t)fMap->GetSignal(idx,l-1);
159	if(l>0) dfadc[k][l-1] = fadc2-fadc1;
160	} // samples
161	} // anodes
162
163	for(k=0;k<fNofAnodes;k++) {
164	idx = j*fNofAnodes+k;
165
166	Int_t imax = 0;
167	Int_t imaxd = 0;
168	Int_t it=0;
169	while(it <= fMaxNofSamples-3) {
170	// cout << "sample: " << it << endl;
171
172	imax = it;
173	imaxd = it;
174	// maximum of signal
175
176	Float_t fadcmax = 0.;
177	Float_t dfadcmax = 0.;
178	Int_t lthrmina = 1;
179	// if(it >= 60) lthrmina = 2;
180	// if(it >= 100) lthrmina = 3;
181	Int_t lthrmint = 2;
182	//if(it >= 60) lthrmint = 3;
183	//if(it >= 100) lthrmint = 4;
184
185	Int_t lthra = 0;
186	Int_t lthrt = 0;
187	for(m=0;m<10;m++) {
188	Int_t id = it+m;
189	if(id>=fMaxNofSamples) break;
190	fadc=fMap->GetSignal(idx,id);
191	if(fadc > fadcmax) {
192	fadcmax = fadc;
193	if(fadc > fCutAmplitude) { lthra++; lthrt++; }
194	imax = id;
195	}
196	if(dfadc[k][id] > dfadcmax) {
197	dfadcmax = dfadc[k][id];
198	imaxd = id;
199	}
200	}
201	it = imaxd;
202	// skip if no signal over threshold
203	if(fMap->GetSignal(idx,imax) < fCutAmplitude) {it++; continue;}
204
205	if(k>0) {
206	if(fMap->GetSignal(idx-1,imax) > fCutAmplitude) lthra++;
207	}
208	if(k<fNofAnodes-1)
209	if(fMap->GetSignal(idx+1,imax) > fCutAmplitude) lthra++;
210
211	if(imax>0) {
212	if(fMap->GetSignal(idx,imax-1) > fCutAmplitude) lthrt++;
213	}
214	if(imax<fMaxNofSamples)
215	if(fMap->GetSignal(idx,imax+1) > fCutAmplitude) lthrt++;
216
217	// cluster charge
218	Int_t tstart = it-1;
219
220	Bool_t ilcl = 0;
221	if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
222
223	if(ilcl) {
224	nofFoundClusters++;
225	Int_t tstop = tstart;
226	Float_t dfadcmin = 10000.;
227	Int_t ij;
228	for(ij=0; ij<10; ij++) {
229	if(dfadc[k][it+ij] < dfadcmin) {
230	tstop = it+ij+1;
231	dfadcmin = dfadc[k][it+ij];
232	}
233	}
234
235	Float_t clusterCharge = 0.;
236	Float_t clusterAnode = k+0.5;
237	Float_t clusterTime = 0.;
238	Float_t clusterMult = 0.;
239	Float_t clusterPeakAmplitude = 0.;
240	Int_t its;
241	for(its=tstart; its<=tstop; its++) {
242	fadc=fMap->GetSignal(idx,its);
243	clusterCharge += fadc;
244	if(fadc > clusterPeakAmplitude) clusterPeakAmplitude = fadc;
245	clusterTime += fadc*its;
246	clusterMult++;
247	if(its == tstop) {
248	clusterTime /= (clusterCharge/fTimeStep); // ns
249	clusterCharge *= (fTimeStep/160.); // keV
250	if(clusterTime > 58.2) clusterTime -= 58.2; // ns
251	/*
252	else {
253	cout << "Warning: cluster with negative time " << clusterTime << ", peak ampl.: " << clusterPeakAmplitude << ", mult: " << clusterMult << ", charge: " << clusterCharge << endl;
254	cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
255	}
256	*/
257	}
258	}
259	// cout << "Anode: " << k << ", tstart: " << tstart << ", tstop: " << tstop << ", Charge: " << clusterCharge << endl;
260
261	Float_t clusteranodePath = (0.06 + clusterAnode - fNofAnodes/2)*anodePitch;
262	Float_t clusterDriftPath = clusterTime*fDriftSpeed;
263	clusterDriftPath = fSddLength-clusterDriftPath;
264
265	if(clusterCharge < 0.) break;
266
267	//printf("wing clusterMult clusterAnode clusterTime %d %f %f %f \n",j+1,clusterMult, clusterAnode, clusterTime);
268
269	AliITSRawClusterSDD *clust = new AliITSRawClusterSDD(j+1,clusterAnode,clusterTime,clusterCharge,clusterPeakAmplitude,0.,0.,clusterDriftPath,clusteranodePath,clusterMult);
270	//fClusters->Add(point);
271	iTS->AddCluster(1,clust);
272	it = tstop;
273	} // ilcl
274
275	it++;
276
277	} // while (samples)
278	} // anodes
279	} // detectors (2)
280
281	Int_t nofClusters = fClusters->GetEntriesFast();
282	nofClusters -= fNclusters;
283
284	//printf("SDD- Find1Dclust: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
285
286	fMap->ClearMap();
287
288	return;
289
290	}
291
292	//_____________________________________________________________________________
293	void AliITSClusterFinderSDD::GroupClusters()
294	{
295	// group clusters
296	Int_t dummy=0;
297	Float_t fTimeStep = fSegmentation->Dpx(dummy);
298
299
300	// get number of clusters for this module
301	Int_t nofClusters = fClusters->GetEntriesFast();
302	nofClusters -= fNclusters;
303
304	//printf("SDD- GroupClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
305
306	AliITSRawClusterSDD *clusterI;
307	AliITSRawClusterSDD *clusterJ;
308
309	Int_t label[nofClusters];
310	Int_t i,j;
311	for(i=0; i<nofClusters; i++) label[i] = 0;
312	for(i=0; i<nofClusters; i++) {
313	if(label[i] != 0) continue;
314	for(j=i+1; j<nofClusters; j++) {
315	if(label[j] != 0) continue;
316	clusterI = (AliITSRawClusterSDD*) fClusters->At(i);
317	clusterJ = (AliITSRawClusterSDD*) fClusters->At(j);
318	// 1.3 good
319	if(clusterI->T() < fTimeStep*60) fDAnode = 3.2;
320	if(clusterI->T() < fTimeStep*10) fDAnode = 1.2;
321	Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
322	if(!pair) continue;
323	// clusterI->Print();
324	// clusterJ->Print();
325	clusterI->Add(clusterJ);
326	label[j] = 1;
327	fClusters->RemoveAt(j);
328	} // J clusters
329	label[i] = 1;
330	} // I clusters
331	fClusters->Compress();
332	return;
333
334	}
335
336	//_____________________________________________________________________________
337
338	void AliITSClusterFinderSDD::SelectClusters()
339	{
340	// get number of clusters for this module
341	Int_t nofClusters = fClusters->GetEntriesFast();
342	nofClusters -= fNclusters;
343
344	//printf("SDD- SelectClusters: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
345
346	Int_t i;
347	for(i=0; i<nofClusters; i++) {
348	AliITSRawClusterSDD clusterI = (AliITSRawClusterSDD) fClusters->At(i);
349	Int_t rmflg = 0;
350	Float_t wy = 0.;
351	if(clusterI->Anodes() != 0.) {
352	wy = ((Float_t) clusterI->Samples())/clusterI->Anodes();
353	}
354	Float_t amp = clusterI->PeakAmpl();
355	if(amp < fMinPeak) rmflg = 1;
356	if(wy < fMinNCells) rmflg = 1;
357	if(rmflg) fClusters->RemoveAt(i);
358	} // I clusters
359	fClusters->Compress();
360	return;
361
362	}
363
364	//_____________________________________________________________________________
365
366	void AliITSClusterFinderSDD::GetRecPoints()
367	{
368	// get rec points
369	//static Int_t counter=0;
370
371	AliITS iTS=(AliITS)gAlice->GetModule("ITS");
372
373	// get number of clusters for this module
374	Int_t nofClusters = fClusters->GetEntriesFast();
375	nofClusters -= fNclusters;
376
377	//printf("SDD- GetRecPoints: fNclusters nofClusters %d %d \n",fNclusters, nofClusters);
378
379	// const Float_t kdEdXtoQ = 2.778e+2; // KeV -> number of e-hole pairs in Si
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
385	Int_t i;
386	for(i=0; i<nofClusters; i++) {
387	AliITSRawClusterSDD clusterI = (AliITSRawClusterSDD) fClusters->At(i);
388	AliITSRecPoint rnew;
389	rnew.SetX(clusterI->X());
390	rnew.SetZ(clusterI->Z());
391	rnew.SetQ(clusterI->Q()); // in KeV - should be ADC
392	//rnew.SetdEdX((clusterI->Q())/kdEdXtoQ);
393	rnew.SetdEdX(kconvGeV*clusterI->Q());
394	rnew.SetSigmaX2(kRMSx*kRMSx);
395	rnew.SetSigmaZ2(kRMSz*kRMSz);
396	rnew.SetProbability(1.);
397	iTS->AddRecPoint(rnew);
398	//counter++;
399	} // I clusters
400	//printf("counter %d\n",counter);
401
402	}
403
404	//_____________________________________________________________________________
405
406	void AliITSClusterFinderSDD::FindRawClusters()
407	{
408	// find raw clusters
409	Find1DClusters();
410	GroupClusters();
411	SelectClusters();
412	GetRecPoints();
413	}
414
415
416
417