[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 "AliITSMapA2.h"
#include "AliITSMapA1.h"
#include "AliITS.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();
    SetNCells();
    SetMap();
}

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

}

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

    if(fMap) delete fMap;
    if(fMapA2) delete fMapA2;

}
//__________________________________________________________________________
AliITSClusterFinderSDD::AliITSClusterFinderSDD(const AliITSClusterFinderSDD &source){
  //     Copy Constructor 
  if(&source == this) return;
  this->fClusters = source.fClusters ;
  this->fNclusters = source.fNclusters ;
  this->fMap = source.fMap ;
  this->fMapA2 = source.fMapA2 ;
  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->fMapA2 = source.fMapA2 ;
  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(!fMapA2) fMapA2=new AliITSMapA2(fSegmentation,fDigits,(double)fCutAmplitude);
    if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);

}
//_____________________________________________________________________________

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
   fMapA2->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;
  Int_t i;
  Float_t **dfadc = new Float_t*[fNofMaps];
  for(i=0;i<fNofMaps;i++) dfadc[i] = new Float_t[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)fMapA2->GetSignal(idx,l);
	fadc1=(Float_t)fMapA2->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=fMapA2->GetSignal(idx,id);
	  if(fadc > fadcmax) {
	    fadcmax = fadc;
	    if(fadc > 0) { lthra++; lthrt++; }
	    imax = id;
	  }
	  if(dfadc[k][id] > dfadcmax) {
	    dfadcmax = dfadc[k][id];
	    imaxd = id;
	  }
	}
	it = imaxd;
	// skip if no signal over threshold
	if(fMapA2->TestHit(idx,imax) == kEmpty) {it++; continue;}

	if(k>0) {
	  if(fMapA2->TestHit(idx-1,imax) != kEmpty) lthra++;
	}
	if(k<fNofAnodes-1)
	  if(fMapA2->TestHit(idx+1,imax) != kEmpty) lthra++;

	if(imax>0) {
	  if(fMapA2->TestHit(idx,imax-1) != kEmpty) lthrt++;
	}
	if(imax<fMaxNofSamples)
	  if(fMapA2->TestHit(idx,imax+1) != kEmpty) 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;
	  Float_t n, baseline;
	  fResponse->GetNoiseParam(n,baseline);
	  for(its=tstart; its<=tstop; its++) {
            fadc=fMapA2->GetSignal(idx,its)-baseline;
	    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)


  fMapA2->ClearMap();
  
  for(i=0;i<fNofMaps;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(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

  fMap->FillMap();

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