[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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.6  2001/08/30 09:28:48  hristov
TTree names are explicitly set via SetName(name) and then Write() is called

Revision 1.5  2001/07/20 14:32:44  kowal2
Processing of many events possible now

Revision 1.4  2001/04/17 08:06:27  hristov
Possibility to define the magnetic field in the reconstruction (Yu.Belikov)

Revision 1.3  2000/10/05 16:14:01  kowal2
Forward declarations.

Revision 1.2  2000/06/30 12:07:50  kowal2
Updated from the TPC-PreRelease branch

Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
Splitted from AliTPCtracking

*/

//-------------------------------------------------------
//          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 "AliTPCClustersArray.h"
#include "AliTPCClustersRow.h"
#include "AliDigits.h"
#include "AliSimDigits.h"
#include "AliTPCParam.h"
#include <iostream.h>
#include <TTree.h>

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);
}

//_____________________________________________________________________________
void AliTPCclusterer::Digits2Clusters(const AliTPCParam *par, TFile *of, Int_t eventn)
{
  //-----------------------------------------------------------------
  // This is a simple cluster finder.
  //-----------------------------------------------------------------
  TDirectory *savedir=gDirectory; 

  if (!of->IsOpen()) {
     cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
     return;
  }

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

  char  dname[100];
  char   cname[100];
  if (eventn==-1) {

    // for backward compatibility
    
    sprintf(dname,"TreeD_75x40_100x60");
    sprintf(cname,"TreeC_TPC");
  }
  else {
    sprintf(dname,"TreeD_75x40_100x60_%d",eventn);
    sprintf(cname,"TreeC_TPC_%d",eventn);
  }
  TTree *t = (TTree *)gDirectory->Get(dname);

  AliSimDigits digarr, *dummy=&digarr;
  t->GetBranch("Segment")->SetAddress(&dummy);
  Stat_t nentries = t->GetEntries();

  of->cd();

  ((AliTPCParam*)par)->Write(par->GetTitle());
  AliTPCClustersArray carray;
  carray.Setup(par);
  carray.SetClusterType("AliTPCcluster");
  carray.MakeTree();


  Int_t nclusters=0;

  for (Int_t n=0; n<nentries; n++) {
    t->GetEvent(n);
    Int_t sec, row;
    if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
       continue;
    }

    AliTPCClustersRow *clrow=carray.CreateRow(sec,row);

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

    Int_t npads, sign;
    {
       const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
       if (sec < kNIS) {
          npads = par->GetNPadsLow(row);
          sign = (sec < kNIS/2) ? 1 : -1;
       } else {
          npads = par->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<=par->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=par->GetPadPitchWidth(sec);
         c.SetSigmaY2((s2 + 1./12.)*w*w);
         if (s2 != 0.) {
	   c.SetSigmaY2(c.GetSigmaY2()*0.108);
	   if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
         }

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

         c.SetY((c.GetY() - 0.5 - 0.5*npads)*par->GetPadPitchWidth(sec));
         c.SetZ(par->GetZWidth()*(c.GetZ()-1)); 
         c.SetZ(c.GetZ() - 3.*par->GetZSigma()); // PASA delay 
         c.SetZ(sign*(par->GetZLength() - 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++;
      }
    }
    carray.StoreRow(sec,row);
    carray.ClearRow(sec,row);

    //cerr<<"sector, row, compressed digits, clusters: "
    //<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<"                  \r";

    nclusters+=ncl;

    delete[] bins;  
  }

  cerr<<"Number of found clusters : "<<nclusters<<"                        \n";

  carray.GetTree()->SetName(cname);
  carray.GetTree()->Write();
  savedir->cd();

  delete t;  //Thanks to Mariana Bondila
}
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
	16	/*
	17	$Log$
f38c8ae5	18	Revision 1.6 2001/08/30 09:28:48 hristov
	19	TTree names are explicitly set via SetName(name) and then Write() is called
	20
6bc3458d	21	Revision 1.5 2001/07/20 14:32:44 kowal2
	22	Processing of many events possible now
	23
afc42102	24	Revision 1.4 2001/04/17 08:06:27 hristov
	25	Possibility to define the magnetic field in the reconstruction (Yu.Belikov)
	26
9b280d80	27	Revision 1.3 2000/10/05 16:14:01 kowal2
	28	Forward declarations.
	29
c3c63118	30	Revision 1.2 2000/06/30 12:07:50 kowal2
	31	Updated from the TPC-PreRelease branch
	32
73042f01	33	Revision 1.1.2.1 2000/06/25 08:53:55 kowal2
	34	Splitted from AliTPCtracking
	35
	36	*/
	37
	38	//-------------------------------------------------------
	39	// Implementation of the TPC clusterer
	40	//
	41	// Origin: Jouri Belikov, CERN, Jouri.Belikov@cern.ch
	42	//-------------------------------------------------------
	43
	44	#include "AliTPCclusterer.h"
	45	#include "AliTPCcluster.h"
	46	#include <TObjArray.h>
	47	#include <TFile.h>
	48	#include "AliTPCClustersArray.h"
	49	#include "AliTPCClustersRow.h"
	50	#include "AliDigits.h"
	51	#include "AliSimDigits.h"
	52	#include "AliTPCParam.h"
	53	#include <iostream.h>
c3c63118	54	#include <TTree.h>
73042f01	55
	56	void AliTPCclusterer::FindPeaks(Int_t k,Int_t max,
	57	AliBin b,Int_t idx,UInt_t *msk,Int_t& n) {
	58	//find local maxima
	59	if (n<31)
	60	if (IsMaximum(k,max,b)) {
	61	idx[n]=k; msk[n]=(2<<n);
	62	n++;
	63	}
	64	b[k].SetMask(0);
	65	if (b[k-max].GetMask()&1) FindPeaks(k-max,max,b,idx,msk,n);
	66	if (b[k-1 ].GetMask()&1) FindPeaks(k-1 ,max,b,idx,msk,n);
	67	if (b[k+max].GetMask()&1) FindPeaks(k+max,max,b,idx,msk,n);
	68	if (b[k+1 ].GetMask()&1) FindPeaks(k+1 ,max,b,idx,msk,n);
	69	}
	70
	71	void AliTPCclusterer::MarkPeak(Int_t k, Int_t max, AliBin *bins, UInt_t m) {
	72	//mark this peak
	73	UShort_t q=bins[k].GetQ();
	74
	75	bins[k].SetMask(bins[k].GetMask()\|m);
	76
	77	if (bins[k-max].GetQ() <= q)
	78	if ((bins[k-max].GetMask()&m) == 0) MarkPeak(k-max,max,bins,m);
	79	if (bins[k-1 ].GetQ() <= q)
	80	if ((bins[k-1 ].GetMask()&m) == 0) MarkPeak(k-1 ,max,bins,m);
	81	if (bins[k+max].GetQ() <= q)
	82	if ((bins[k+max].GetMask()&m) == 0) MarkPeak(k+max,max,bins,m);
	83	if (bins[k+1 ].GetQ() <= q)
	84	if ((bins[k+1 ].GetMask()&m) == 0) MarkPeak(k+1 ,max,bins,m);
	85	}
	86
	87	void AliTPCclusterer::MakeCluster(Int_t k,Int_t max,AliBin *bins,UInt_t m,
	88	AliTPCcluster &c) {
	89	//make cluster using digits of this peak
	90	Float_t q=(Float_t)bins[k].GetQ();
	91	Int_t i=k/max, j=k-i*max;
	92
	93	c.SetQ(c.GetQ()+q);
	94	c.SetY(c.GetY()+i*q);
	95	c.SetZ(c.GetZ()+j*q);
	96	c.SetSigmaY2(c.GetSigmaY2()+iiq);
	97	c.SetSigmaZ2(c.GetSigmaZ2()+jjq);
	98
	99	bins[k].SetMask(0xFFFFFFFE);
	100
	101	if (bins[k-max].GetMask() == m) MakeCluster(k-max,max,bins,m,c);
	102	if (bins[k-1 ].GetMask() == m) MakeCluster(k-1 ,max,bins,m,c);
	103	if (bins[k+max].GetMask() == m) MakeCluster(k+max,max,bins,m,c);
	104	if (bins[k+1 ].GetMask() == m) MakeCluster(k+1 ,max,bins,m,c);
	105	}
	106
	107	//_____________________________________________________________________________
afc42102	108	void AliTPCclusterer::Digits2Clusters(const AliTPCParam par, TFile of, Int_t eventn)
73042f01	109	{
	110	//-----------------------------------------------------------------
	111	// This is a simple cluster finder.
	112	//-----------------------------------------------------------------
	113	TDirectory *savedir=gDirectory;
	114
	115	if (!of->IsOpen()) {
	116	cerr<<"AliTPC::Digits2Clusters(): output file not open !\n";
	117	return;
	118	}
	119
	120	const Int_t kMAXZ=par->GetMaxTBin()+2;
	121
afc42102	122	char dname[100];
	123	char cname[100];
	124	if (eventn==-1) {
	125
	126	// for backward compatibility
	127
	128	sprintf(dname,"TreeD_75x40_100x60");
	129	sprintf(cname,"TreeC_TPC");
	130	}
	131	else {
	132	sprintf(dname,"TreeD_75x40_100x60_%d",eventn);
	133	sprintf(cname,"TreeC_TPC_%d",eventn);
	134	}
	135	TTree t = (TTree )gDirectory->Get(dname);
	136
73042f01	137	AliSimDigits digarr, *dummy=&digarr;
	138	t->GetBranch("Segment")->SetAddress(&dummy);
	139	Stat_t nentries = t->GetEntries();
	140
	141	of->cd();
	142
	143	((AliTPCParam*)par)->Write(par->GetTitle());
	144	AliTPCClustersArray carray;
	145	carray.Setup(par);
	146	carray.SetClusterType("AliTPCcluster");
	147	carray.MakeTree();
	148
afc42102	149
afc42102	150
73042f01	151	Int_t nclusters=0;
	152
	153	for (Int_t n=0; n<nentries; n++) {
	154	t->GetEvent(n);
	155	Int_t sec, row;
	156	if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
	157	cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
	158	continue;
	159	}
	160
	161	AliTPCClustersRow *clrow=carray.CreateRow(sec,row);
	162
	163	Float_t rx=par->GetPadRowRadii(sec,row);
	164
	165	Int_t npads, sign;
	166	{
	167	const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
	168	if (sec < kNIS) {
	169	npads = par->GetNPadsLow(row);
	170	sign = (sec < kNIS/2) ? 1 : -1;
	171	} else {
	172	npads = par->GetNPadsUp(row);
	173	sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
	174	}
	175	}
	176
	177	const Int_t kMAXBIN=kMAXZ*(npads+2);
	178	AliBin *bins=new AliBin[kMAXBIN];
	179	for (Int_t ii=0;ii<kMAXBIN;ii++) {
	180	bins[ii].SetQ(0); bins[ii].SetMask(0xFFFFFFFE);
	181	}
	182
	183	digarr.First();
	184	do {
	185	Short_t dig=digarr.CurrentDigit();
	186	if (dig<=par->GetZeroSup()) continue;
	187	Int_t j=digarr.CurrentRow()+1, i=digarr.CurrentColumn()+1;
	188	bins[i*kMAXZ+j].SetQ(dig);
	189	bins[i*kMAXZ+j].SetMask(1);
	190	} while (digarr.Next());
	191
	192	Int_t ncl=0;
	193	for (Int_t i=0; i<kMAXBIN; i++) {
	194	if ((bins[i].GetMask()&1) == 0) continue;
	195	Int_t idx[32]; UInt_t msk[32]; Int_t npeaks=0;
	196	FindPeaks(i, kMAXZ, bins, idx, msk, npeaks);
	197
	198	if (npeaks>30) continue;
	199
	200	Int_t k,l;
	201	for (k=0; k<npeaks-1; k++){//mark adjacent peaks
	202	if (idx[k] < 0) continue; //this peak is already removed
	203	for (l=k+1; l<npeaks; l++) {
	204	if (idx[l] < 0) continue; //this peak is already removed
	205	Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
	206	Int_t li=idx[l]/kMAXZ, lj=idx[l] - li*kMAXZ;
	207	Int_t di=TMath::Abs(ki - li);
	208	Int_t dj=TMath::Abs(kj - lj);
	209	if (di>1 \|\| dj>1) continue;
	210	if (bins[idx[k]].GetQ() > bins[idx[l]].GetQ()) {
	211	msk[l]=msk[k];
	212	idx[l]*=-1;
	213	} else {
	214	msk[k]=msk[l];
215	idx[k]*=-1;
216	break;
217	}
218	}
219	}
220
221	for (k=0; k<npeaks; k++) {
222	MarkPeak(TMath::Abs(idx[k]), kMAXZ, bins, msk[k]);
223	}
224
225	for (k=0; k<npeaks; k++) {
226	if (idx[k] < 0) continue; //removed peak
227	AliTPCcluster c;
228	MakeCluster(idx[k], kMAXZ, bins, msk[k], c);
229	if (c.GetQ() < 5) continue; //noise cluster
230	c.SetY(c.GetY()/c.GetQ());
231	c.SetZ(c.GetZ()/c.GetQ());
232
233	Float_t s2 = c.GetSigmaY2()/c.GetQ() - c.GetY()*c.GetY();
234	Float_t w=par->GetPadPitchWidth(sec);
235	c.SetSigmaY2((s2 + 1./12.)ww);
236	if (s2 != 0.) {
237	c.SetSigmaY2(c.GetSigmaY2()*0.108);
238	if (sec<par->GetNInnerSector()) c.SetSigmaY2(c.GetSigmaY2()*2.07);
239	}
240
241	s2 = c.GetSigmaZ2()/c.GetQ() - c.GetZ()*c.GetZ();
242	w=par->GetZWidth();
243	c.SetSigmaZ2((s2 + 1./12.)ww);
244	if (s2 != 0.) {
245	c.SetSigmaZ2(c.GetSigmaZ2()*0.169);
246	if (sec<par->GetNInnerSector()) c.SetSigmaZ2(c.GetSigmaZ2()*1.77);
247	}
248
249	c.SetY((c.GetY() - 0.5 - 0.5npads)par->GetPadPitchWidth(sec));
250	c.SetZ(par->GetZWidth()*(c.GetZ()-1));
251	c.SetZ(c.GetZ() - 3.*par->GetZSigma()); // PASA delay
252	c.SetZ(sign*(par->GetZLength() - c.GetZ()));
253
254	if (rx<230./250.*TMath::Abs(c.GetZ())) continue;
255
256	Int_t ki=idx[k]/kMAXZ, kj=idx[k] - ki*kMAXZ;
257	c.SetLabel(digarr.GetTrackID(kj-1,ki-1,0),0);
258	c.SetLabel(digarr.GetTrackID(kj-1,ki-1,1),1);
259	c.SetLabel(digarr.GetTrackID(kj-1,ki-1,2),2);
260
261	c.SetQ(bins[idx[k]].GetQ());
262
263	if (ki==1 \|\| ki==npads \|\| kj==1 \|\| kj==kMAXZ-2) {
264	c.SetSigmaY2(c.GetSigmaY2()*25.);
265	c.SetSigmaZ2(c.GetSigmaZ2()*4.);
266	}
267	clrow->InsertCluster(&c); ncl++;
268	}
269	}
270	carray.StoreRow(sec,row);
271	carray.ClearRow(sec,row);
272
9b280d80	273	//cerr<<"sector, row, compressed digits, clusters: "
9b280d80	274	//<<sec<<' '<<row<<' '<<digarr.GetSize()<<' '<<ncl<<" \r";
73042f01	275
	276	nclusters+=ncl;
	277
	278	delete[] bins;
	279	}
	280
	281	cerr<<"Number of found clusters : "<<nclusters<<" \n";
	282
6bc3458d	283	carray.GetTree()->SetName(cname);
6bc3458d	284	carray.GetTree()->Write();
73042f01	285	savedir->cd();
f38c8ae5	286
f38c8ae5	287	delete t; //Thanks to Mariana Bondila
73042f01	288	}
73042f01	289