[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.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.3  2000/05/08 16:17:27  cblume
Merge TRD-develop

Revision 1.1.4.1  2000/05/08 15:09:01  cblume
Introduce AliTRDdigitsManager

Revision 1.1  2000/02/28 18:58:54  cblume
Add new TRD classes

*/

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// TRD cluster finder for the slow simulator. 
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TF1.h>

#include "AliTRDclusterizerV1.h"
#include "AliTRDmatrix.h"
#include "AliTRDgeometry.h"
#include "AliTRDdigitizer.h"
#include "AliTRDrecPoint.h"
#include "AliTRDdataArrayF.h"

ClassImp(AliTRDclusterizerV1)

//_____________________________________________________________________________
AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
{
  //
  // AliTRDclusterizerV1 default constructor
  //

  fDigitsManager = NULL;

}

//_____________________________________________________________________________
AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
                    :AliTRDclusterizer(name,title)
{
  //
  // AliTRDclusterizerV1 default constructor
  //

  fDigitsManager = new AliTRDdigitsManager();

  Init();

}

//_____________________________________________________________________________
AliTRDclusterizerV1::~AliTRDclusterizerV1()
{

  if (fDigitsManager) {
    delete fDigitsManager;
  }

}

//_____________________________________________________________________________
void AliTRDclusterizerV1::Init()
{
  //
  // Initializes the cluster finder
  //

  // The default parameter for the clustering
  fClusMaxThresh = 5.0;
  fClusSigThresh = 2.0;
  fClusMethod    = 1;

}

//_____________________________________________________________________________
Bool_t AliTRDclusterizerV1::ReadDigits()
{
  //
  // Reads the digits arrays from the input aliroot file
  //

  if (!fInputFile) {
    printf("AliTRDclusterizerV1::ReadDigits -- ");
    printf("No input file open\n");
    return kFALSE;
  }

  // Read in the digit arrays
  return (fDigitsManager->ReadDigits());  

}

//_____________________________________________________________________________
Bool_t AliTRDclusterizerV1::MakeCluster()
{
  //
  // Generates the cluster.
  //

  Int_t row, col, time;

  // Get the pointer to the detector class and check for version 1
  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
  if (TRD->IsVersion() != 1) {
    printf("AliTRDclusterizerV1::MakeCluster -- ");
    printf("TRD must be version 1 (slow simulator).\n");
    return kFALSE; 
  }

  // Get the geometry
  AliTRDgeometry *Geo = TRD->GetGeometry();

  printf("AliTRDclusterizerV1::MakeCluster -- ");
  printf("Start creating clusters.\n");

  AliTRDdataArrayI *Digits;

  // Parameters
  Float_t maxThresh        = fClusMaxThresh;   // threshold value for maximum
  Float_t signalThresh     = fClusSigThresh;   // threshold value for digit signal
  Int_t   clusteringMethod = fClusMethod;      // clustering method option (for testing)

  // Iteration limit for unfolding procedure
  const Float_t epsilon = 0.01;             

  const Int_t   nClus   = 3;  
  const Int_t   nSig    = 5;

  Int_t chamBeg = 0;
  Int_t chamEnd = kNcham;
  if (TRD->GetSensChamber()  >= 0) {
    chamBeg = TRD->GetSensChamber();
    chamEnd = chamBeg + 1;
  }
  Int_t planBeg = 0;
  Int_t planEnd = kNplan;
  if (TRD->GetSensPlane()    >= 0) {
    planBeg = TRD->GetSensPlane();
    planEnd = planBeg + 1;
  }
  Int_t sectBeg = 0;
  Int_t sectEnd = kNsect;

  // *** Start clustering *** in every chamber
  for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
    for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
      for (Int_t isect = sectBeg; isect < sectEnd; isect++) {

        if (TRD->GetSensSector() >= 0) {
          Int_t sens1 = TRD->GetSensSector();
          Int_t sens2 = sens1 + TRD->GetSensSectorRange();
          sens2 -= ((Int_t) (sens2 / kNsect)) * kNsect;
          if (sens1 < sens2) 
            if ((isect < sens1) || (isect >= sens2)) continue;
          else
            if ((isect < sens1) && (isect >= sens2)) continue;
	}

        Int_t idet = Geo->GetDetector(iplan,icham,isect);

        Int_t nClusters = 0;
        printf("AliTRDclusterizerV1::MakeCluster -- ");
        printf("Analyzing chamber %d, plane %d, sector %d.\n"
               ,icham,iplan,isect);

        Int_t   nRowMax  = Geo->GetRowMax(iplan,icham,isect);
        Int_t   nColMax  = Geo->GetColMax(iplan);
        Int_t   nTimeMax = Geo->GetTimeMax();

        // Create a detector matrix to keep maxima
        AliTRDmatrix *digitMatrix  = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
                                                     ,isect,icham,iplan);
        // Create a matrix to contain maximum flags
        AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
                                                     ,isect,icham,iplan);

        // Read in the digits
        Digits = fDigitsManager->GetDigits(idet);

        // Loop through the detector pixel
        for (time = 0; time < nTimeMax; time++) {
          for ( col = 0;  col <  nColMax;  col++) {
            for ( row = 0;  row <  nRowMax;  row++) {

              Int_t signal = Digits->GetData(row,col,time);
              Int_t index  = Digits->GetIndex(row,col,time);

              // Fill the detector matrix
              if (signal > signalThresh) {
	        // Store the signal amplitude
                digitMatrix->SetSignal(row,col,time,signal);
	        // Store the digits number
                digitMatrix->AddTrack(row,col,time,index);
              }

	    }
	  }
	}

        // Loop chamber and find maxima in digitMatrix
        for ( row = 0;  row <  nRowMax;  row++) {
          for ( col = 1;  col <  nColMax;  col++) {
            for (time = 0; time < nTimeMax; time++) {

              if (digitMatrix->GetSignal(row,col,time) 
                  < digitMatrix->GetSignal(row,col - 1,time)) {
                // really maximum?
                if (col > 1) {
                  if (digitMatrix->GetSignal(row,col - 2,time)
                      < digitMatrix->GetSignal(row,col - 1,time)) {
                    // yes, so set maximum flag
                    maximaMatrix->SetSignal(row,col - 1,time,1);
                  }
                  else maximaMatrix->SetSignal(row,col - 1,time,0);
                }
              }

            }   // time
          }     // col
        }       // row

        // now check maxima and calculate cluster position
        for ( row = 0;  row <  nRowMax;  row++) {
          for ( col = 1;  col <  nColMax;  col++) {
            for (time = 0; time < nTimeMax; time++) {

              if ((maximaMatrix->GetSignal(row,col,time) > 0)
                  && (digitMatrix->GetSignal(row,col,time) > maxThresh)) {

                // Ratio resulting from unfolding
                Float_t ratio                =  0;    
                // Signals on max and neighbouring pads
                Float_t padSignal[nSig]      = {0};   
                // Signals from cluster
                Float_t clusterSignal[nClus] = {0};
                // Cluster pad info
                Float_t clusterPads[nClus]   = {0};   
                // Cluster digit info
                Int_t   clusterDigit[nClus]  = {0};

                Int_t iPad;
                for (iPad = 0; iPad < nClus; iPad++) {
                  clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
                  clusterDigit[iPad]  = digitMatrix->GetTrack(row,col-1+iPad,time,0);
                }

                // neighbouring maximum on right side?
                if (col < nColMax - 2) {
                  if (maximaMatrix->GetSignal(row,col + 2,time) > 0) {

                    for (iPad = 0; iPad < 5; iPad++) {
                      padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
                    }

                    // unfold:
                    ratio = Unfold(epsilon, padSignal);

                    // set signal on overlapping pad to ratio
                    clusterSignal[2] *= ratio;

                  }
                }
                
		// Calculate the position of the cluster
                switch (clusteringMethod) {
                case 1:
                  // method 1: simply center of mass
                  clusterPads[0] = row + 0.5;
                  clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) /
                                   (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]);
                  clusterPads[2] = time + 0.5;

                  nClusters++;
                  break;
                case 2:
                  // method 2: integral gauss fit on 3 pads
                  TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads"
	                                                    , 5, -1.5, 3.5);
                  for (Int_t iCol = -1; iCol <= 3; iCol++) {
                    if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1;
                    hPadCharges->Fill(iCol, clusterSignal[iCol]);
                  }
                  hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5);
                  TF1     *fPadChargeFit = hPadCharges->GetFunction("gaus");
                  Double_t  colMean = fPadChargeFit->GetParameter(1);

                  clusterPads[0] = row + 0.5;
                  clusterPads[1] = col - 1.5 + colMean;
                  clusterPads[2] = time + 0.5;

                  delete hPadCharges;

                  nClusters++;
                  break;
                }

                Float_t clusterCharge =   clusterSignal[0]
                                        + clusterSignal[1]
                                        + clusterSignal[2];

                // Add the cluster to the output array 
                TRD->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);

              }
            }  // time
          }    // col
        }      // row

        printf("AliTRDclusterizerV1::MakeCluster -- ");
        printf("Number of clusters found: %d\n",nClusters);

        delete digitMatrix;
        delete maximaMatrix;

      }          // isect
    }            // iplan
  }              // icham

  printf("AliTRDclusterizerV1::MakeCluster -- ");
  printf("Total number of points found: %d\n"
        ,TRD->RecPoints()->GetEntries());

  // Get the pointer to the cluster branch
  TTree *ClusterTree = gAlice->TreeR(); 

  // Fill the cluster-branch
  printf("AliTRDclusterizerV1::MakeCluster -- ");
  printf("Fill the cluster tree.\n");
  ClusterTree->Fill();
  printf("AliTRDclusterizerV1::MakeCluster -- ");
  printf("Done.\n");

  return kTRUE;

}

//_____________________________________________________________________________
Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
{
  //
  // Method to unfold neighbouring maxima.
  // The charge ratio on the overlapping pad is calculated
  // until there is no more change within the range given by eps.
  // The resulting ratio is then returned to the calling method.
  //

  Int_t   itStep            = 0;      // count iteration steps

  Float_t ratio             = 0.5;    // start value for ratio
  Float_t prevRatio         = 0;      // store previous ratio

  Float_t newLeftSignal[3]  = {0};    // array to store left cluster signal
  Float_t newRightSignal[3] = {0};    // array to store right cluster signal

  // start iteration:
  while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {

    itStep++;
    prevRatio = ratio;

    // cluster position according to charge ratio
    Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) /
                       (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
    Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) /
                       ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);

    // set cluster charge ratio
    Float_t ampLeft  = padSignal[1];
    Float_t ampRight = padSignal[3];

    // apply pad response to parameters
    newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft);
    newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft);
    newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft);

    newRightSignal[0] = ampRight*PadResponse(-1 - maxRight);
    newRightSignal[1] = ampRight*PadResponse( 0 - maxRight);
    newRightSignal[2] = ampRight*PadResponse( 1 - maxRight);

    // calculate new overlapping ratio
    ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]);

  }

  return ratio;

}

//_____________________________________________________________________________
Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
{
  //
  // The pad response for the chevron pads. 
  // We use a simple Gaussian approximation which should be good
  // enough for our purpose.
  //

  // The parameters for the response function
  const Float_t aa  =  0.8872;
  const Float_t bb  = -0.00573;
  const Float_t cc  =  0.454;
  const Float_t cc2 =  cc*cc;

  Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2)));

  return (pr);

}
Commit	Line	Data
f7336fa3	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$
9d0b222b	18	Revision 1.3 2000/05/08 16:17:27 cblume
	19	Merge TRD-develop
	20
6f1e466d	21	Revision 1.1.4.1 2000/05/08 15:09:01 cblume
	22	Introduce AliTRDdigitsManager
	23
c0dd96c3	24	Revision 1.1 2000/02/28 18:58:54 cblume
	25	Add new TRD classes
	26
f7336fa3	27	*/
	28
	29	///////////////////////////////////////////////////////////////////////////////
	30	// //
	31	// TRD cluster finder for the slow simulator.
	32	// //
	33	///////////////////////////////////////////////////////////////////////////////
	34
	35	#include <TF1.h>
	36
	37	#include "AliTRDclusterizerV1.h"
	38	#include "AliTRDmatrix.h"
	39	#include "AliTRDgeometry.h"
	40	#include "AliTRDdigitizer.h"
	41	#include "AliTRDrecPoint.h"
6f1e466d	42	#include "AliTRDdataArrayF.h"
f7336fa3	43
	44	ClassImp(AliTRDclusterizerV1)
	45
	46	//_____________________________________________________________________________
	47	AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
	48	{
	49	//
	50	// AliTRDclusterizerV1 default constructor
	51	//
	52
6f1e466d	53	fDigitsManager = NULL;
f7336fa3	54
	55	}
	56
	57	//_____________________________________________________________________________
	58	AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
	59	:AliTRDclusterizer(name,title)
	60	{
	61	//
	62	// AliTRDclusterizerV1 default constructor
	63	//
	64
6f1e466d	65	fDigitsManager = new AliTRDdigitsManager();
f7336fa3	66
	67	Init();
	68
	69	}
	70
	71	//_____________________________________________________________________________
	72	AliTRDclusterizerV1::~AliTRDclusterizerV1()
	73	{
	74
6f1e466d	75	if (fDigitsManager) {
6f1e466d	76	delete fDigitsManager;
f7336fa3	77	}
	78
	79	}
	80
	81	//_____________________________________________________________________________
	82	void AliTRDclusterizerV1::Init()
	83	{
	84	//
	85	// Initializes the cluster finder
	86	//
	87
	88	// The default parameter for the clustering
	89	fClusMaxThresh = 5.0;
	90	fClusSigThresh = 2.0;
	91	fClusMethod = 1;
	92
	93	}
	94
	95	//_____________________________________________________________________________
	96	Bool_t AliTRDclusterizerV1::ReadDigits()
	97	{
	98	//
	99	// Reads the digits arrays from the input aliroot file
	100	//
	101
	102	if (!fInputFile) {
	103	printf("AliTRDclusterizerV1::ReadDigits -- ");
	104	printf("No input file open\n");
	105	return kFALSE;
	106	}
	107
f7336fa3	108	// Read in the digit arrays
6f1e466d	109	return (fDigitsManager->ReadDigits());
f7336fa3	110
	111	}
	112
	113	//_____________________________________________________________________________
	114	Bool_t AliTRDclusterizerV1::MakeCluster()
	115	{
	116	//
	117	// Generates the cluster.
	118	//
	119
	120	Int_t row, col, time;
	121
	122	// Get the pointer to the detector class and check for version 1
	123	AliTRD TRD = (AliTRD) gAlice->GetDetector("TRD");
	124	if (TRD->IsVersion() != 1) {
	125	printf("AliTRDclusterizerV1::MakeCluster -- ");
	126	printf("TRD must be version 1 (slow simulator).\n");
	127	return kFALSE;
	128	}
	129
	130	// Get the geometry
	131	AliTRDgeometry *Geo = TRD->GetGeometry();
	132
	133	printf("AliTRDclusterizerV1::MakeCluster -- ");
	134	printf("Start creating clusters.\n");
	135
6f1e466d	136	AliTRDdataArrayI *Digits;
f7336fa3	137
	138	// Parameters
	139	Float_t maxThresh = fClusMaxThresh; // threshold value for maximum
	140	Float_t signalThresh = fClusSigThresh; // threshold value for digit signal
	141	Int_t clusteringMethod = fClusMethod; // clustering method option (for testing)
	142
	143	// Iteration limit for unfolding procedure
	144	const Float_t epsilon = 0.01;
	145
	146	const Int_t nClus = 3;
	147	const Int_t nSig = 5;
	148
	149	Int_t chamBeg = 0;
	150	Int_t chamEnd = kNcham;
9d0b222b	151	if (TRD->GetSensChamber() >= 0) {
f7336fa3	152	chamBeg = TRD->GetSensChamber();
6f1e466d	153	chamEnd = chamBeg + 1;
f7336fa3	154	}
	155	Int_t planBeg = 0;
	156	Int_t planEnd = kNplan;
9d0b222b	157	if (TRD->GetSensPlane() >= 0) {
f7336fa3	158	planBeg = TRD->GetSensPlane();
	159	planEnd = planBeg + 1;
	160	}
	161	Int_t sectBeg = 0;
	162	Int_t sectEnd = kNsect;
f7336fa3	163
	164	// * Start clustering * in every chamber
	165	for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
	166	for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
	167	for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
	168
9d0b222b	169	if (TRD->GetSensSector() >= 0) {
	170	Int_t sens1 = TRD->GetSensSector();
	171	Int_t sens2 = sens1 + TRD->GetSensSectorRange();
	172	sens2 -= ((Int_t) (sens2 / kNsect)) * kNsect;
	173	if (sens1 < sens2)
	174	if ((isect < sens1) \|\| (isect >= sens2)) continue;
	175	else
	176	if ((isect < sens1) && (isect >= sens2)) continue;
	177	}
	178
f7336fa3	179	Int_t idet = Geo->GetDetector(iplan,icham,isect);
	180
	181	Int_t nClusters = 0;
	182	printf("AliTRDclusterizerV1::MakeCluster -- ");
	183	printf("Analyzing chamber %d, plane %d, sector %d.\n"
	184	,icham,iplan,isect);
	185
	186	Int_t nRowMax = Geo->GetRowMax(iplan,icham,isect);
	187	Int_t nColMax = Geo->GetColMax(iplan);
	188	Int_t nTimeMax = Geo->GetTimeMax();
	189
	190	// Create a detector matrix to keep maxima
	191	AliTRDmatrix *digitMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
	192	,isect,icham,iplan);
	193	// Create a matrix to contain maximum flags
	194	AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
	195	,isect,icham,iplan);
	196
	197	// Read in the digits
6f1e466d	198	Digits = fDigitsManager->GetDigits(idet);
f7336fa3	199
	200	// Loop through the detector pixel
	201	for (time = 0; time < nTimeMax; time++) {
	202	for ( col = 0; col < nColMax; col++) {
	203	for ( row = 0; row < nRowMax; row++) {
	204
	205	Int_t signal = Digits->GetData(row,col,time);
	206	Int_t index = Digits->GetIndex(row,col,time);
	207
	208	// Fill the detector matrix
	209	if (signal > signalThresh) {
	210	// Store the signal amplitude
	211	digitMatrix->SetSignal(row,col,time,signal);
	212	// Store the digits number
	213	digitMatrix->AddTrack(row,col,time,index);
	214	}
	215
	216	}
	217	}
	218	}
	219
	220	// Loop chamber and find maxima in digitMatrix
	221	for ( row = 0; row < nRowMax; row++) {
	222	for ( col = 1; col < nColMax; col++) {
	223	for (time = 0; time < nTimeMax; time++) {
	224
	225	if (digitMatrix->GetSignal(row,col,time)
	226	< digitMatrix->GetSignal(row,col - 1,time)) {
	227	// really maximum?
	228	if (col > 1) {
	229	if (digitMatrix->GetSignal(row,col - 2,time)
	230	< digitMatrix->GetSignal(row,col - 1,time)) {
	231	// yes, so set maximum flag
	232	maximaMatrix->SetSignal(row,col - 1,time,1);
	233	}
	234	else maximaMatrix->SetSignal(row,col - 1,time,0);
	235	}
	236	}
	237
	238	} // time
	239	} // col
	240	} // row
	241
	242	// now check maxima and calculate cluster position
	243	for ( row = 0; row < nRowMax; row++) {
	244	for ( col = 1; col < nColMax; col++) {
	245	for (time = 0; time < nTimeMax; time++) {
	246
	247	if ((maximaMatrix->GetSignal(row,col,time) > 0)
	248	&& (digitMatrix->GetSignal(row,col,time) > maxThresh)) {
	249
	250	// Ratio resulting from unfolding
	251	Float_t ratio = 0;
	252	// Signals on max and neighbouring pads
	253	Float_t padSignal[nSig] = {0};
	254	// Signals from cluster
	255	Float_t clusterSignal[nClus] = {0};
	256	// Cluster pad info
	257	Float_t clusterPads[nClus] = {0};
	258	// Cluster digit info
	259	Int_t clusterDigit[nClus] = {0};
	260
9d0b222b	261	Int_t iPad;
9d0b222b	262	for (iPad = 0; iPad < nClus; iPad++) {
f7336fa3	263	clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
	264	clusterDigit[iPad] = digitMatrix->GetTrack(row,col-1+iPad,time,0);
	265	}
	266
	267	// neighbouring maximum on right side?
	268	if (col < nColMax - 2) {
	269	if (maximaMatrix->GetSignal(row,col + 2,time) > 0) {
	270
9d0b222b	271	for (iPad = 0; iPad < 5; iPad++) {
f7336fa3	272	padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
	273	}
	274
	275	// unfold:
	276	ratio = Unfold(epsilon, padSignal);
	277
	278	// set signal on overlapping pad to ratio
	279	clusterSignal[2] *= ratio;
	280
	281	}
	282	}
	283
	284	// Calculate the position of the cluster
	285	switch (clusteringMethod) {
	286	case 1:
	287	// method 1: simply center of mass
	288	clusterPads[0] = row + 0.5;
	289	clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) /
c0dd96c3	290	(clusterSignal[0] + clusterSignal[1] + clusterSignal[2]);
f7336fa3	291	clusterPads[2] = time + 0.5;
	292
	293	nClusters++;
	294	break;
	295	case 2:
	296	// method 2: integral gauss fit on 3 pads
	297	TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads"
	298	, 5, -1.5, 3.5);
	299	for (Int_t iCol = -1; iCol <= 3; iCol++) {
	300	if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1;
	301	hPadCharges->Fill(iCol, clusterSignal[iCol]);
	302	}
	303	hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5);
	304	TF1 *fPadChargeFit = hPadCharges->GetFunction("gaus");
	305	Double_t colMean = fPadChargeFit->GetParameter(1);
	306
	307	clusterPads[0] = row + 0.5;
	308	clusterPads[1] = col - 1.5 + colMean;
	309	clusterPads[2] = time + 0.5;
	310
	311	delete hPadCharges;
	312
	313	nClusters++;
	314	break;
	315	}
	316
	317	Float_t clusterCharge = clusterSignal[0]
	318	+ clusterSignal[1]
	319	+ clusterSignal[2];
	320
	321	// Add the cluster to the output array
	322	TRD->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);
	323
	324	}
	325	} // time
	326	} // col
	327	} // row
	328
	329	printf("AliTRDclusterizerV1::MakeCluster -- ");
	330	printf("Number of clusters found: %d\n",nClusters);
	331
	332	delete digitMatrix;
	333	delete maximaMatrix;
	334
	335	} // isect
	336	} // iplan
	337	} // icham
	338
	339	printf("AliTRDclusterizerV1::MakeCluster -- ");
	340	printf("Total number of points found: %d\n"
	341	,TRD->RecPoints()->GetEntries());
	342
	343	// Get the pointer to the cluster branch
	344	TTree *ClusterTree = gAlice->TreeR();
	345
	346	// Fill the cluster-branch
	347	printf("AliTRDclusterizerV1::MakeCluster -- ");
	348	printf("Fill the cluster tree.\n");
	349	ClusterTree->Fill();
	350	printf("AliTRDclusterizerV1::MakeCluster -- ");
	351	printf("Done.\n");
	352
	353	return kTRUE;
	354
355	}
356
357	//_____________________________________________________________________________
358	Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
359	{
360	//
361	// Method to unfold neighbouring maxima.
362	// The charge ratio on the overlapping pad is calculated
363	// until there is no more change within the range given by eps.
364	// The resulting ratio is then returned to the calling method.
365	//
366
367	Int_t itStep = 0; // count iteration steps
368
369	Float_t ratio = 0.5; // start value for ratio
370	Float_t prevRatio = 0; // store previous ratio
371
372	Float_t newLeftSignal[3] = {0}; // array to store left cluster signal
373	Float_t newRightSignal[3] = {0}; // array to store right cluster signal
374
375	// start iteration:
376	while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
377
378	itStep++;
379	prevRatio = ratio;
380
381	// cluster position according to charge ratio
382	Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) /
383	(padSignal[0] + padSignal[1] + ratio*padSignal[2]);
384	Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) /
385	((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
386
387	// set cluster charge ratio
388	Float_t ampLeft = padSignal[1];
389	Float_t ampRight = padSignal[3];
390
391	// apply pad response to parameters
392	newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft);
393	newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft);
394	newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft);
395
396	newRightSignal[0] = ampRight*PadResponse(-1 - maxRight);
397	newRightSignal[1] = ampRight*PadResponse( 0 - maxRight);
398	newRightSignal[2] = ampRight*PadResponse( 1 - maxRight);
399
400	// calculate new overlapping ratio
401	ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]);
402
403	}
404
405	return ratio;
406
407	}
408
409	//_____________________________________________________________________________
410	Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
411	{
412	//
413	// The pad response for the chevron pads.
414	// We use a simple Gaussian approximation which should be good
415	// enough for our purpose.
416	//
417
418	// The parameters for the response function
419	const Float_t aa = 0.8872;
420	const Float_t bb = -0.00573;
421	const Float_t cc = 0.454;
422	const Float_t cc2 = cc*cc;
423
424	Float_t pr = aa * (bb + TMath::Exp(-xx / (2. cc2)));
425
426	return (pr);
427
428	}