[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSPD.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 "AliITSClusterFinderSPD.h"
#include "AliRun.h"


ClassImp(AliITSClusterFinderSPD)

//----------------------------------------------------------
AliITSClusterFinderSPD::AliITSClusterFinderSPD
(AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
{
  // constructor
    fSegmentation=seg;
    fDigits=digits;
    fClusters=recp;
    fNclusters= fClusters->GetEntriesFast();
    SetDx();
    SetDz();
    SetMap();
    SetNCells();
}

//_____________________________________________________________________________
AliITSClusterFinderSPD::AliITSClusterFinderSPD()
{
  // constructor
  fSegmentation=0;
  fDigits=0;
  fClusters=0;
  fNclusters=0;
  SetDx();
  SetDz();
  SetMap();
  SetNCells();
  
}

//__________________________________________________________________________
AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source){
  //     Copy Constructor 
  if(&source == this) return;
  this->fClusters = source.fClusters ;
  this->fNclusters = source.fNclusters ;
  this->fMap = source.fMap ;
  this->fDz = source.fDz ;
  this->fDx = source.fDx ;
  this->fMinNCells = source.fMinNCells ;
  return;
}

//_________________________________________________________________________
AliITSClusterFinderSPD& 
  AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD &source) {
  //    Assignment operator
  if(&source == this) return *this;
  this->fClusters = source.fClusters ;
  this->fNclusters = source.fNclusters ;
  this->fMap = source.fMap ;
  this->fDz = source.fDz ;
  this->fDx = source.fDx ;
  this->fMinNCells = source.fMinNCells ;
  return *this;
}

//_____________________________________________________________________________
void AliITSClusterFinderSPD::SetMap()
{
  // set map
  if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
  
}

//_____________________________________________________________________________

void AliITSClusterFinderSPD::Find1DClusters()
{
  // Find one dimensional clusters, i.e.
  // in r*phi(x) direction for each colunm in z direction
  
  AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
  
  // retrieve the parameters 
  Int_t fNofPixels = fSegmentation->Npz(); 
  Int_t fMaxNofSamples = fSegmentation->Npx();
  
  // read in digits -> do not apply threshold 
  // signal in fired pixels is always 1
  fMap->FillMap();
  
  Int_t nofFoundClusters = 0;
  
  Int_t k,it,m;
  for(k=0;k<fNofPixels;k++) {
    
    Int_t mmax = 10;  // a size of the window for the cluster finding
    
    for(it=0;it<fMaxNofSamples;it++) {
      
      Int_t lclx = 0;
      Int_t xstart = 0;
      Int_t xstop = 0;
      Int_t id = 0;
      Int_t ilcl =0;
      
      for(m=0;m<mmax;m++) {  // find the cluster inside the window
	id = it+m;
	if(id >= fMaxNofSamples) break;    // ! no possible for the fadc 
	
	if(fMap->GetHitIndex(k,id)) {   // start of the cluster
	  lclx += 1;
	  if(lclx == 1) xstart = id;
	  
	}
	
	if(lclx > 0 && !fMap->GetHitIndex(k,id)) {  
	  // end of cluster if a gap exists
	  xstop = id-1;
	  ilcl = 1;
	  break;
	}            
	
      }   //  end of m-loop
      
      if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
      // continue the "it" loop
      
      if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
      
      if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {  
	// cluster end is outside of the window,
	mmax += 5;                 // increase mmax and repeat the cluster
	// finding
	it -= 1;
      }
      
      if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
	xstop = fMaxNofSamples - 1;          // the end cluster exists
	ilcl = 1;
      } 
      
      // ---  Calculate z and x coordinates for one dimensional clusters
      
      if(ilcl == 1) {         // new cluster exists
	it = id;
	mmax = 10;
	    nofFoundClusters++;
	    Float_t clusterCharge = 0.;
            // get this from segmentation when this will be implemented
            Float_t zpitch = fSegmentation->Dpz(k+1); 
	    Float_t clusterZ, dummyX; 
            Int_t dummy=0;
            fSegmentation->GetCellCxz(dummy,k,dummyX,clusterZ);
            Float_t zstart = clusterZ - 0.5*zpitch;
            Float_t zstop = clusterZ + 0.5*zpitch;
	    Float_t clusterX = 0.;
            Int_t xstartfull = xstart;
            Int_t xstopfull = xstop;
	    Int_t clusterSizeX = lclx;
	    Int_t clusterSizeZ = 1;
	    
            Int_t its;
	    for(its=xstart; its<=xstop; its++) {
              Int_t firedpixel=0;
              if (fMap->GetHitIndex(k,its)) firedpixel=1; 
	      clusterCharge += firedpixel;
              clusterX +=its + 0.5;
	    }
            Float_t fRphiPitch = fSegmentation->Dpx(dummy);
	    clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x 
	    
	    
	    //printf("ClusterZ ClusterX %f %f \n",clusterZ, clusterX);
	    
            // Int_t nclusters = fClusters->GetEntriesFast();
	    //	    	    cout << nclusters << " clusters" << endl;
	    //            cout<< "Create point"<<endl;
	    
	    // Write the points (coordinates and some cluster information) to the
	    // AliITSRawClusterSPD object
	    
	    AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
	    // fClusters->Add(point);
	    iTS->AddCluster(0,clust);
	    //    	    cout << "Cluster at Ladder: " << fLadder << ", Detector: " <<fDetector<<endl;
	    
	    // cout<<" end of cluster finding for Z pixel "<<endl;
	    
      }    // new cluster (ilcl=1)
    } // X direction loop (it)
  } // Z direction loop (k)
  return;
  
}

//_____________________________________________________________________________
void  AliITSClusterFinderSPD::GroupClusters()
{
  // Find two dimensional clusters, i.e. group one dimensional clusters
  // into two dimensional ones (go both in x and z directions).
  
  // get number of clusters for this module
  Int_t nofClusters = fClusters->GetEntriesFast();
  nofClusters -= fNclusters;
  //printf("Group: fNclusters nofClusters %d %d\n",fNclusters, nofClusters);
  
  AliITSRawClusterSPD *clusterI;
  AliITSRawClusterSPD *clusterJ;
  
  Int_t *label = new Int_t [nofClusters];  // activate this for DEC machines
  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 = (AliITSRawClusterSPD*) fClusters->At(i);
      clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
      Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
      if(pair) {     
	
	//    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
	/*    
	      cout << "clusters " << i << "," << j << " before grouping" << endl;
	      clusterI->Print();
	      clusterJ->Print();
	*/    
	clusterI->Add(clusterJ);
	//        cout << "remove cluster " << j << endl;
	label[j] = 1;
	fClusters->RemoveAt(j);
	/*
	  cout << "cluster  " << i << " after grouping" << endl;
	  clusterI->Print();
	*/
      }  // pair
    } // J clusters  
    label[i] = 1;
  } // I clusters
  fClusters->Compress();
  //Int_t totalNofClusters = fClusters->GetEntriesFast();
  //cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
  
  delete [] label;
  
  return;
  
  
}
//_____________________________________________________________________________

void AliITSClusterFinderSPD::GetRecPoints()
{
  // get rec points
  AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
  
  // get number of clusters for this module
  Int_t nofClusters = fClusters->GetEntriesFast();
  nofClusters -= fNclusters;
  //printf("GetRecP: fNclusters nofClusters %d %d\n",fNclusters, nofClusters);
  
  const Float_t kconv = 1.0e-4;
  const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
  const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3

  Int_t i;
  for(i=0; i<nofClusters; i++) { 
    AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
    AliITSRecPoint rnew;
    rnew.SetX(clusterI->X()*kconv);
    rnew.SetZ(clusterI->Z()*kconv);
    rnew.SetQ(1.);
    rnew.SetdEdX(0.);
    rnew.SetSigmaX2(kRMSx*kRMSx);
    rnew.SetSigmaZ2(kRMSz*kRMSz);
    rnew.SetProbability(1.);
    iTS->AddRecPoint(rnew);
  } // I clusters
  
}
//_____________________________________________________________________________

void AliITSClusterFinderSPD::FindRawClusters()
{
  // find raw clusters
  Find1DClusters();
  GroupClusters();
  GetRecPoints();
}
Commit	Line	Data
	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 "AliITSClusterFinderSPD.h"
	18	#include "AliRun.h"
	19
	20
	21
	22	ClassImp(AliITSClusterFinderSPD)
	23
	24	//----------------------------------------------------------
	25	AliITSClusterFinderSPD::AliITSClusterFinderSPD
	26	(AliITSsegmentation seg, TClonesArray digits, TClonesArray *recp)
	27	{
	28	// constructor
	29	fSegmentation=seg;
	30	fDigits=digits;
	31	fClusters=recp;
	32	fNclusters= fClusters->GetEntriesFast();
	33	SetDx();
	34	SetDz();
	35	SetMap();
	36	SetNCells();
	37	}
	38
	39	//_____________________________________________________________________________
	40	AliITSClusterFinderSPD::AliITSClusterFinderSPD()
	41	{
	42	// constructor
	43	fSegmentation=0;
	44	fDigits=0;
	45	fClusters=0;
	46	fNclusters=0;
	47	SetDx();
	48	SetDz();
	49	SetMap();
	50	SetNCells();
	51
	52	}
	53
	54	//__________________________________________________________________________
	55	AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source){
	56	// Copy Constructor
	57	if(&source == this) return;
	58	this->fClusters = source.fClusters ;
	59	this->fNclusters = source.fNclusters ;
	60	this->fMap = source.fMap ;
	61	this->fDz = source.fDz ;
	62	this->fDx = source.fDx ;
	63	this->fMinNCells = source.fMinNCells ;
	64	return;
	65	}
	66
	67	//_________________________________________________________________________
	68	AliITSClusterFinderSPD&
	69	AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD &source) {
	70	// Assignment operator
	71	if(&source == this) return *this;
	72	this->fClusters = source.fClusters ;
	73	this->fNclusters = source.fNclusters ;
	74	this->fMap = source.fMap ;
	75	this->fDz = source.fDz ;
	76	this->fDx = source.fDx ;
	77	this->fMinNCells = source.fMinNCells ;
	78	return *this;
	79	}
	80
	81	//_____________________________________________________________________________
	82	void AliITSClusterFinderSPD::SetMap()
	83	{
	84	// set map
	85	if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
	86
	87	}
	88
	89	//_____________________________________________________________________________
	90
	91	void AliITSClusterFinderSPD::Find1DClusters()
	92	{
	93	// Find one dimensional clusters, i.e.
	94	// in r*phi(x) direction for each colunm in z direction
	95
	96	AliITS iTS=(AliITS)gAlice->GetModule("ITS");
	97
	98	// retrieve the parameters
	99	Int_t fNofPixels = fSegmentation->Npz();
	100	Int_t fMaxNofSamples = fSegmentation->Npx();
	101
	102	// read in digits -> do not apply threshold
	103	// signal in fired pixels is always 1
	104	fMap->FillMap();
	105
	106	Int_t nofFoundClusters = 0;
	107
	108	Int_t k,it,m;
	109	for(k=0;k<fNofPixels;k++) {
	110
	111	Int_t mmax = 10; // a size of the window for the cluster finding
	112
	113	for(it=0;it<fMaxNofSamples;it++) {
	114
	115	Int_t lclx = 0;
	116	Int_t xstart = 0;
	117	Int_t xstop = 0;
	118	Int_t id = 0;
	119	Int_t ilcl =0;
	120
	121	for(m=0;m<mmax;m++) { // find the cluster inside the window
	122	id = it+m;
	123	if(id >= fMaxNofSamples) break; // ! no possible for the fadc
	124
	125	if(fMap->GetHitIndex(k,id)) { // start of the cluster
	126	lclx += 1;
	127	if(lclx == 1) xstart = id;
	128
	129	}
	130
	131	if(lclx > 0 && !fMap->GetHitIndex(k,id)) {
	132	// end of cluster if a gap exists
	133	xstop = id-1;
	134	ilcl = 1;
	135	break;
	136	}
	137
	138	} // end of m-loop
	139
	140	if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
	141	// continue the "it" loop
	142
	143	if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
	144
	145	if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {
	146	// cluster end is outside of the window,
	147	mmax += 5; // increase mmax and repeat the cluster
	148	// finding
	149	it -= 1;
	150	}
	151
	152	if(id >= fMaxNofSamples && lclx > 0) { // the x row finished but
	153	xstop = fMaxNofSamples - 1; // the end cluster exists
	154	ilcl = 1;
	155	}
	156
	157	// --- Calculate z and x coordinates for one dimensional clusters
	158
	159	if(ilcl == 1) { // new cluster exists
	160	it = id;
	161	mmax = 10;
	162	nofFoundClusters++;
	163	Float_t clusterCharge = 0.;
	164	// get this from segmentation when this will be implemented
	165	Float_t zpitch = fSegmentation->Dpz(k+1);
	166	Float_t clusterZ, dummyX;
	167	Int_t dummy=0;
	168	fSegmentation->GetCellCxz(dummy,k,dummyX,clusterZ);
	169	Float_t zstart = clusterZ - 0.5*zpitch;
	170	Float_t zstop = clusterZ + 0.5*zpitch;
	171	Float_t clusterX = 0.;
	172	Int_t xstartfull = xstart;
	173	Int_t xstopfull = xstop;
	174	Int_t clusterSizeX = lclx;
	175	Int_t clusterSizeZ = 1;
	176
	177	Int_t its;
	178	for(its=xstart; its<=xstop; its++) {
	179	Int_t firedpixel=0;
	180	if (fMap->GetHitIndex(k,its)) firedpixel=1;
	181	clusterCharge += firedpixel;
	182	clusterX +=its + 0.5;
	183	}
	184	Float_t fRphiPitch = fSegmentation->Dpx(dummy);
	185	clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x
	186
	187
	188	//printf("ClusterZ ClusterX %f %f \n",clusterZ, clusterX);
	189
	190	// Int_t nclusters = fClusters->GetEntriesFast();
	191	// cout << nclusters << " clusters" << endl;
	192	// cout<< "Create point"<<endl;
	193
	194	// Write the points (coordinates and some cluster information) to the
	195	// AliITSRawClusterSPD object
	196
	197	AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
	198	// fClusters->Add(point);
	199	iTS->AddCluster(0,clust);
	200	// cout << "Cluster at Ladder: " << fLadder << ", Detector: " <<fDetector<<endl;
	201
	202	// cout<<" end of cluster finding for Z pixel "<<endl;
	203
	204	} // new cluster (ilcl=1)
	205	} // X direction loop (it)
	206	} // Z direction loop (k)
	207	return;
	208
	209	}
	210
	211	//_____________________________________________________________________________
	212	void AliITSClusterFinderSPD::GroupClusters()
	213	{
	214	// Find two dimensional clusters, i.e. group one dimensional clusters
	215	// into two dimensional ones (go both in x and z directions).
	216
	217	// get number of clusters for this module
	218	Int_t nofClusters = fClusters->GetEntriesFast();
	219	nofClusters -= fNclusters;
	220	//printf("Group: fNclusters nofClusters %d %d\n",fNclusters, nofClusters);
	221
	222	AliITSRawClusterSPD *clusterI;
	223	AliITSRawClusterSPD *clusterJ;
	224
	225	Int_t *label = new Int_t [nofClusters]; // activate this for DEC machines
	226	Int_t i,j;
	227	for(i=0; i<nofClusters; i++) label[i] = 0;
	228	for(i=0; i<nofClusters; i++) {
	229	if(label[i] != 0) continue;
	230	for(j=i+1; j<nofClusters; j++) {
	231	if(label[j] != 0) continue;
	232	clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
	233	clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
	234	Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
	235	if(pair) {
	236
	237	// if((clusterI->XStop() == clusterJ->XStart()-1)\|\|(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
	238	/*
	239	cout << "clusters " << i << "," << j << " before grouping" << endl;
	240	clusterI->Print();
	241	clusterJ->Print();
	242	*/
	243	clusterI->Add(clusterJ);
	244	// cout << "remove cluster " << j << endl;
	245	label[j] = 1;
	246	fClusters->RemoveAt(j);
	247	/*
	248	cout << "cluster " << i << " after grouping" << endl;
	249	clusterI->Print();
	250	*/
	251	} // pair
	252	} // J clusters
	253	label[i] = 1;
	254	} // I clusters
	255	fClusters->Compress();
	256	//Int_t totalNofClusters = fClusters->GetEntriesFast();
	257	//cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
	258
	259	delete [] label;
	260
	261	return;
	262
	263
	264	}
	265	//_____________________________________________________________________________
	266
	267	void AliITSClusterFinderSPD::GetRecPoints()
	268	{
	269	// get rec points
	270	AliITS iTS=(AliITS)gAlice->GetModule("ITS");
	271
	272	// get number of clusters for this module
	273	Int_t nofClusters = fClusters->GetEntriesFast();
	274	nofClusters -= fNclusters;
	275	//printf("GetRecP: fNclusters nofClusters %d %d\n",fNclusters, nofClusters);
	276
	277	const Float_t kconv = 1.0e-4;
	278	const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
	279	const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
	280
	281	Int_t i;
	282	for(i=0; i<nofClusters; i++) {
	283	AliITSRawClusterSPD clusterI = (AliITSRawClusterSPD) fClusters->At(i);
	284	AliITSRecPoint rnew;
	285	rnew.SetX(clusterI->X()*kconv);
	286	rnew.SetZ(clusterI->Z()*kconv);
	287	rnew.SetQ(1.);
	288	rnew.SetdEdX(0.);
	289	rnew.SetSigmaX2(kRMSx*kRMSx);
	290	rnew.SetSigmaZ2(kRMSz*kRMSz);
	291	rnew.SetProbability(1.);
	292	iTS->AddRecPoint(rnew);
	293	} // I clusters
	294
	295	}
	296	//_____________________________________________________________________________
	297
	298	void AliITSClusterFinderSPD::FindRawClusters()
	299	{
	300	// find raw clusters
	301	Find1DClusters();
	302	GroupClusters();
	303	GetRecPoints();
	304	}
	305