[u/mrichter/AliRoot.git] / TPC / AliTPCclusterer.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.                  *
 **************************************************************************/

/* $Id$ */

//-------------------------------------------------------
//          Implementation of the TPC clusterer
//
//   Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch 
//-------------------------------------------------------

#include "AliTPCclusterer.h"
#include "AliTPCcluster.h"
#include <TObjArray.h>
#include <TFile.h>
#include "AliDigits.h"
#include "AliSimDigits.h"
#include "AliTPCParam.h"
#include "AliTPCClustersRow.h"
#include <TTree.h>

ClassImp(AliTPCclusterer)


//____________________________________________________
AliTPCclusterer::AliTPCclusterer(const AliTPCclusterer &param)
                :TObject(),
                 fPar(0)
{
  //
  // dummy
  //
  fPar = param.fPar;
}
//-----------------------------------------------------
AliTPCclusterer & AliTPCclusterer::operator =(const AliTPCclusterer & param)
{
  //
  // assignment operator - dummy
  //
  if (this == &param) return (*this);

  fPar = param.fPar;
  return (*this);
}
//____________________________________________________
void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
AliBin *b,Int_t *idx,UInt_t *msk,Int_t& n) {
  //find local maxima
  if (n<31)
  if (IsMaximum(k,max,b)) {
    idx[n]=k; msk[n]=(2<<n);
    n++;
  }
  b[k].SetMask(0);
  if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
  if (b[k-1  ].GetMask()&1) FindPeaks(k-1  ,max,b,idx,msk,n);
  if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
  if (b[k+1  ].GetMask()&1) FindPeaks(k+1  ,max,b,idx,msk,n);
}

void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
  //mark this peak
  UShort_t q=bins[k].GetQ();

  bins[k].SetMask(bins[k].GetMask()|m); 

  if (bins[k-max].GetQ() <= q)
     if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
  if (bins[k-1  ].GetQ() <= q)
     if ((bins[k-1  ].GetMask()&m) == 0) MarkPeak(k-1  ,max,bins,m);
  if (bins[k+max].GetQ() <= q)
     if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
  if (bins[k+1  ].GetQ() <= q)
     if ((bins[k+1  ].GetMask()&m) == 0) MarkPeak(k+1  ,max,bins,m);
}

void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
AliTPCcluster &c) {
  //make cluster using digits of this peak
  Float_t q=(Float_t)bins[k].GetQ();
  Int_t i=k/max, j=k-i*max;

  c.SetQ(c.GetQ()+q);
  c.SetY(c.GetY()+i*q); 
  c.SetZ(c.GetZ()+j*q); 
  c.SetSigmaY2(c.GetSigmaY2()+i*i*q);
  c.SetSigmaZ2(c.GetSigmaZ2()+j*j*q);

  bins[k].SetMask(0xFFFFFFFE);
  
  if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
  if (bins[k-1  ].GetMask() == m) MakeCluster(k-1  ,max,bins,m,c);
  if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
  if (bins[k+1  ].GetMask() == m) MakeCluster(k+1  ,max,bins,m,c);
}

//_____________________________________________________________________________
Int_t AliTPCclusterer::Digits2Clusters(TTree *dTree, TTree *cTree) {
  //-----------------------------------------------------------------
  // This is a simple cluster finder.
  //-----------------------------------------------------------------
  TBranch *branch=dTree->GetBranch("Segment");
  if (!branch) {
     Error("Digits2Cluster","Can't get the branch !");
     return 1;
  }
  AliSimDigits digarr, *dummy=&digarr;
  branch->SetAddress(&dummy);
  
  AliTPCClustersRow ddd,*clrow=&ddd;
  clrow = new AliTPCClustersRow("AliTPCcluster");
  cTree->Branch("Segment","AliTPCClustersRow",&clrow,32000,200);    

  const Int_t kMAXZ=fPar->GetMaxTBin()+2;

  Int_t nclusters=0;

  Int_t nentries = (Int_t)dTree->GetEntries();
  for (Int_t n=0; n<nentries; n++) {
   
    Int_t sec, row;
    dTree->GetEvent(n);

    if (!fPar->AdjustSectorRow(digarr.GetID(),sec,row)) {
       Error("Digits2Clusters","!nvalid segment ID ! %d",digarr.GetID());
       continue;
    }

    clrow->SetID(digarr.GetID());

    cTree->GetBranch("Segment")->SetAddress(&clrow);

    Float_t rx=fPar->GetPadRowRadii(sec,row);

    Int_t npads, sign;
    {
       const Int_t kNIS=fPar->GetNInnerSector(), kNOS=fPar->GetNOuterSector();
       if (sec < kNIS) {
          npads = fPar->GetNPadsLow(row);
          sign = (sec < kNIS/2) ? 1 : -1;
       } else {
          npads = fPar->GetNPadsUp(row);
          sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
       }
    }

    const Int_t kMAXBIN=kMAXZ*(npads+2);
    AliBin *bins=new AliBin[kMAXBIN];
    for (Int_t ii=0;ii<kMAXBIN;ii++) {
       bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
    }

    digarr.First();
    do {
       Short_t dig=digarr.CurrentDigit();
       if (dig<=fPar->GetZeroSup()) continue;
       Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
       bins[i*kMAXZ+j].SetQ(dig);
       bins[i*kMAXZ+j].SetMask(1);
    } while (digarr.Next());

    Int_t ncl=0;
    for (Int_t i=0; i<kMAXBIN; i++) {
      if ((bins[i].GetMask()&1) == 0) continue;
      Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
      FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);

      if (npeaks>30) continue;

      Int_t k,l;
      for (k=0; k<npeaks-1; k++){//mark adjacent peaks
        if (idx[k] < 0) continue; //this peak is already removed
        for (l=k+1; l<npeaks; l++) {
           if (idx[l] < 0) continue; //this peak is already removed
           Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
           Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
           Int_t di=TMath::Abs(ki - li);
           Int_t dj=TMath::Abs(kj - lj);
           if (di>1 || dj>1) continue;
           if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
              msk[l]=msk[k];
              idx[l]*=-1;
           } else {
              msk[k]=msk[l];
              idx[k]*=-1;
              break;
           } 
        }
      }

      for (k=0; k<npeaks; k++) {
        MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
      }
        
      for (k=0; k<npeaks; k++) {
         if (idx[k] < 0) continue; //removed peak
         AliTPCcluster c;
         MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
         if (c.GetQ() < 5) continue; //noise cluster
         c.SetY(c.GetY()/c.GetQ());
         c.SetZ(c.GetZ()/c.GetQ());

         Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
         Float_t w=fPar->GetPadPitchWidth(sec);
         c.SetSigmaY2((s2 + 1./12.)*w*w);
         if (s2 != 0.) {
	   c.SetSigmaY2(c.GetSigmaY2()*0.108);
	   if (sec<fPar->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
         }

         s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
         w=fPar->GetZWidth();
         c.SetSigmaZ2((s2 + 1./12.)*w*w);
         if (s2 != 0.) {
	   c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
	   if (sec<fPar->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
         }

         c.SetY((c.GetY() - 0.5 - 0.5*npads)*fPar->GetPadPitchWidth(sec));
         c.SetZ(fPar->GetZWidth()*(c.GetZ()-1)); 
         c.SetZ(c.GetZ() - 3.*fPar->GetZSigma()); // PASA delay 
         c.SetZ(sign*(fPar->GetZLength(sec) - c.GetZ()));

         if (rx<230./250.*TMath::Abs(c.GetZ())) continue;

         Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
         c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
         c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
         c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);

         c.SetQ(bins[idx[k]].GetQ());

         if (ki==1 || ki==npads || kj==1 || kj==kMAXZ-2) {
           c.SetSigmaY2(c.GetSigmaY2()*25.);
           c.SetSigmaZ2(c.GetSigmaZ2()*4.);
         }
         clrow->InsertCluster(&c); ncl++;
      }
    }
    cTree->Fill();

    clrow->GetArray()->Clear();

    nclusters+=ncl;

    delete[] bins;  
  }

  Info("Digits2Cluster","Number of found clusters : %d",nclusters);

  return 0;
}
Commit	Line	Data
73042f01	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
88cb7938	16	/* $Id$ */
73042f01	17
	18	//-------------------------------------------------------
	19	// Implementation of the TPC clusterer
	20	//
	21	// Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
	22	//-------------------------------------------------------
	23
	24	#include "AliTPCclusterer.h"
	25	#include "AliTPCcluster.h"
	26	#include <TObjArray.h>
	27	#include <TFile.h>
73042f01	28	#include "AliDigits.h"
	29	#include "AliSimDigits.h"
	30	#include "AliTPCParam.h"
d7cbafb5	31	#include "AliTPCClustersRow.h"
c3c63118	32	#include <TTree.h>
c630aafd	33
	34	ClassImp(AliTPCclusterer)
	35
73042f01	36
179c6296	37	//____________________________________________________
	38	AliTPCclusterer::AliTPCclusterer(const AliTPCclusterer &param)
	39	:TObject(),
	40	fPar(0)
	41	{
	42	//
	43	// dummy
	44	//
	45	fPar = param.fPar;
	46	}
	47	//-----------------------------------------------------
	48	AliTPCclusterer & AliTPCclusterer::operator =(const AliTPCclusterer & param)
	49	{
	50	//
	51	// assignment operator - dummy
	52	//
04420071	53	if (this == &param) return (*this);
	54
	55	fPar = param.fPar;
179c6296	56	return (*this);
	57	}
	58	//____________________________________________________
73042f01	59	void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
	60	AliBin b,Int_t idx,UInt_t *msk,Int_t& n) {
	61	//find local maxima
	62	if (n<31)
	63	if (IsMaximum(k,max,b)) {
	64	idx[n]=k; msk[n]=(2<<n);
	65	n++;
	66	}
	67	b[k].SetMask(0);
	68	if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
	69	if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
	70	if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
	71	if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
	72	}
	73
	74	void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
	75	//mark this peak
	76	UShort_t q=bins[k].GetQ();
	77
	78	bins[k].SetMask(bins[k].GetMask()\|m);
	79
	80	if (bins[k-max].GetQ() <= q)
	81	if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
	82	if (bins[k-1 ].GetQ() <= q)
	83	if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
	84	if (bins[k+max].GetQ() <= q)
	85	if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
	86	if (bins[k+1 ].GetQ() <= q)
	87	if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
	88	}
	89
	90	void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
	91	AliTPCcluster &c) {
	92	//make cluster using digits of this peak
	93	Float_t q=(Float_t)bins[k].GetQ();
	94	Int_t i=k/max, j=k-i*max;
	95
	96	c.SetQ(c.GetQ()+q);
	97	c.SetY(c.GetY()+i*q);
	98	c.SetZ(c.GetZ()+j*q);
	99	c.SetSigmaY2(c.GetSigmaY2()+iiq);
	100	c.SetSigmaZ2(c.GetSigmaZ2()+jjq);
	101
	102	bins[k].SetMask(0xFFFFFFFE);
	103
	104	if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
	105	if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
	106	if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
	107	if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
	108	}
	109
	110	//_____________________________________________________________________________
c630aafd	111	Int_t AliTPCclusterer::Digits2Clusters(TTree dTree, TTree cTree) {
73042f01	112	//-----------------------------------------------------------------
	113	// This is a simple cluster finder.
	114	//-----------------------------------------------------------------
c630aafd	115	TBranch *branch=dTree->GetBranch("Segment");
	116	if (!branch) {
	117	Error("Digits2Cluster","Can't get the branch !");
	118	return 1;
6e4d905e	119	}
73042f01	120	AliSimDigits digarr, *dummy=&digarr;
c630aafd	121	branch->SetAddress(&dummy);
88cb7938	122
d7cbafb5	123	AliTPCClustersRow ddd,*clrow=&ddd;
e7de6656	124	clrow = new AliTPCClustersRow("AliTPCcluster");
d7cbafb5	125	cTree->Branch("Segment","AliTPCClustersRow",&clrow,32000,200);
afc42102	126
d7cbafb5	127	const Int_t kMAXZ=fPar->GetMaxTBin()+2;
afc42102	128
73042f01	129	Int_t nclusters=0;
73042f01	130
c630aafd	131	Int_t nentries = (Int_t)dTree->GetEntries();
c630aafd	132	for (Int_t n=0; n<nentries; n++) {
88cb7938	133
73042f01	134	Int_t sec, row;
c630aafd	135	dTree->GetEvent(n);
88cb7938	136
d7cbafb5	137	if (!fPar->AdjustSectorRow(digarr.GetID(),sec,row)) {
c630aafd	138	Error("Digits2Clusters","!nvalid segment ID ! %d",digarr.GetID());
73042f01	139	continue;
	140	}
	141
d7cbafb5	142	clrow->SetID(digarr.GetID());
73042f01	143
d7cbafb5	144	cTree->GetBranch("Segment")->SetAddress(&clrow);
	145
	146	Float_t rx=fPar->GetPadRowRadii(sec,row);
73042f01	147
	148	Int_t npads, sign;
	149	{
d7cbafb5	150	const Int_t kNIS=fPar->GetNInnerSector(), kNOS=fPar->GetNOuterSector();
73042f01	151	if (sec < kNIS) {
d7cbafb5	152	npads = fPar->GetNPadsLow(row);
73042f01	153	sign = (sec < kNIS/2) ? 1 : -1;
73042f01	154	} else {
d7cbafb5	155	npads = fPar->GetNPadsUp(row);
73042f01	156	sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
	157	}
	158	}
	159
	160	const Int_t kMAXBIN=kMAXZ*(npads+2);
	161	AliBin *bins=new AliBin[kMAXBIN];
	162	for (Int_t ii=0;ii<kMAXBIN;ii++) {
	163	bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
	164	}
	165
	166	digarr.First();
	167	do {
	168	Short_t dig=digarr.CurrentDigit();
d7cbafb5	169	if (dig<=fPar->GetZeroSup()) continue;
73042f01	170	Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
	171	bins[i*kMAXZ+j].SetQ(dig);
	172	bins[i*kMAXZ+j].SetMask(1);
	173	} while (digarr.Next());
	174
	175	Int_t ncl=0;
	176	for (Int_t i=0; i<kMAXBIN; i++) {
	177	if ((bins[i].GetMask()&1) == 0) continue;
	178	Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
	179	FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
	180
	181	if (npeaks>30) continue;
	182
	183	Int_t k,l;
	184	for (k=0; k<npeaks-1; k++){//mark adjacent peaks
	185	if (idx[k] < 0) continue; //this peak is already removed
	186	for (l=k+1; l<npeaks; l++) {
	187	if (idx[l] < 0) continue; //this peak is already removed
	188	Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
	189	Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
	190	Int_t di=TMath::Abs(ki - li);
	191	Int_t dj=TMath::Abs(kj - lj);
	192	if (di>1 \|\| dj>1) continue;
	193	if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
	194	msk[l]=msk[k];
	195	idx[l]*=-1;
	196	} else {
	197	msk[k]=msk[l];
	198	idx[k]*=-1;
	199	break;
	200	}
	201	}
	202	}
	203
	204	for (k=0; k<npeaks; k++) {
	205	MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
	206	}
	207
	208	for (k=0; k<npeaks; k++) {
	209	if (idx[k] < 0) continue; //removed peak
	210	AliTPCcluster c;
	211	MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
	212	if (c.GetQ() < 5) continue; //noise cluster
	213	c.SetY(c.GetY()/c.GetQ());
	214	c.SetZ(c.GetZ()/c.GetQ());
	215
	216	Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
d7cbafb5	217	Float_t w=fPar->GetPadPitchWidth(sec);
73042f01	218	c.SetSigmaY2((s2 + 1./12.)ww);
	219	if (s2 != 0.) {
	220	c.SetSigmaY2(c.GetSigmaY2()*0.108);
d7cbafb5	221	if (sec<fPar->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
73042f01	222	}
	223
	224	s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
d7cbafb5	225	w=fPar->GetZWidth();
73042f01	226	c.SetSigmaZ2((s2 + 1./12.)ww);
	227	if (s2 != 0.) {
	228	c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
d7cbafb5	229	if (sec<fPar->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
73042f01	230	}
73042f01	231
d7cbafb5	232	c.SetY((c.GetY() - 0.5 - 0.5npads)fPar->GetPadPitchWidth(sec));
	233	c.SetZ(fPar->GetZWidth()*(c.GetZ()-1));
	234	c.SetZ(c.GetZ() - 3.*fPar->GetZSigma()); // PASA delay
a1ec4d07	235	c.SetZ(sign*(fPar->GetZLength(sec) - c.GetZ()));
73042f01	236
	237	if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
	238
	239	Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
	240	c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
	241	c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
	242	c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
	243
	244	c.SetQ(bins[idx[k]].GetQ());
	245
	246	if (ki==1 \|\| ki==npads \|\| kj==1 \|\| kj==kMAXZ-2) {
	247	c.SetSigmaY2(c.GetSigmaY2()*25.);
	248	c.SetSigmaZ2(c.GetSigmaZ2()*4.);
	249	}
	250	clrow->InsertCluster(&c); ncl++;
	251	}
	252	}
d7cbafb5	253	cTree->Fill();
d7cbafb5	254
e7de6656	255	clrow->GetArray()->Clear();
73042f01	256
	257	nclusters+=ncl;
	258
	259	delete[] bins;
	260	}
	261
c630aafd	262	Info("Digits2Cluster","Number of found clusters : %d",nclusters);
f38c8ae5	263
c630aafd	264	return 0;
73042f01	265	}
73042f01	266