[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.9  2000/11/01 14:53:20  cblume
Merge with TRD-develop

Revision 1.1.4.5  2000/10/15 23:40:01  cblume
Remove AliTRDconst

Revision 1.1.4.4  2000/10/06 16:49:46  cblume
Made Getters const

Revision 1.1.4.3  2000/10/04 16:34:58  cblume
Replace include files by forward declarations

Revision 1.1.4.2  2000/09/22 14:49:49  cblume
Adapted to tracking code

Revision 1.8  2000/10/02 21:28:19  fca
Removal of useless dependecies via forward declarations

Revision 1.7  2000/06/27 13:08:50  cblume
Changed to Copy(TObject &A) to appease the HP-compiler

Revision 1.6  2000/06/09 11:10:07  cblume
Compiler warnings and coding conventions, next round

Revision 1.5  2000/06/08 18:32:58  cblume
Make code compliant to coding conventions

Revision 1.4  2000/06/07 16:27:01  cblume
Try to remove compiler warnings on Sun and HP

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 <TTree.h>
#include <TH1.h>

#include "AliRun.h"

#include "AliTRD.h"
#include "AliTRDclusterizerV1.h"
#include "AliTRDmatrix.h"
#include "AliTRDgeometry.h"
#include "AliTRDdigitizer.h"
#include "AliTRDdataArrayF.h"
#include "AliTRDdataArrayI.h"
#include "AliTRDdigitsManager.h"

ClassImp(AliTRDclusterizerV1)

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

  fDigitsManager = NULL;

  fClusMaxThresh = 0;
  fClusSigThresh = 0;

}

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

  fDigitsManager = new AliTRDdigitsManager();

  Init();

}

//_____________________________________________________________________________
AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
{
  //
  // AliTRDclusterizerV1 copy constructor
  //

  ((AliTRDclusterizerV1 &) c).Copy(*this);

}

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

  if (fDigitsManager) {
    delete fDigitsManager;
  }

}

//_____________________________________________________________________________
AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
{
  //
  // Assignment operator
  //

  if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
  return *this;

}

//_____________________________________________________________________________
void AliTRDclusterizerV1::Copy(TObject &c)
{
  //
  // Copy function
  //

  ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
  ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
  ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;

  AliTRDclusterizer::Copy(c);

}

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

  // The default parameter for the clustering
  fClusMaxThresh = 3;
  fClusSigThresh = 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::MakeClusters()
{
  //
  // Generates the cluster.
  //

  Int_t row, col, time;

  if (fTRD->IsVersion() != 1) {
    printf("AliTRDclusterizerV1::MakeCluster -- ");
    printf("TRD must be version 1 (slow simulator).\n");
    return kFALSE; 
  }

  // Get the geometry
  AliTRDgeometry *geo = fTRD->GetGeometry();

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

  AliTRDdataArrayI *digits;
  AliTRDdataArrayI *track0;
  AliTRDdataArrayI *track1;
  AliTRDdataArrayI *track2; 

  // Threshold value for the maximum
  Int_t maxThresh = fClusMaxThresh;   
  // Threshold value for the digit signal
  Int_t sigThresh = fClusSigThresh;   

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

  const Int_t   kNclus   = 3;  
  const Int_t   kNsig    = 5;
  const Int_t   kNtrack  = 3 * kNclus;

  Float_t padSignal[kNsig];   
  Float_t clusterSignal[kNclus];
  Float_t clusterPads[kNclus];   
  Int_t   clusterDigit[kNclus];
  Int_t   clusterTracks[kNtrack];   

  Int_t chamBeg = 0;
  Int_t chamEnd = AliTRDgeometry::Ncham();
  if (fTRD->GetSensChamber()  >= 0) {
    chamBeg = fTRD->GetSensChamber();
    chamEnd = chamBeg + 1;
  }
  Int_t planBeg = 0;
  Int_t planEnd = AliTRDgeometry::Nplan();
  if (fTRD->GetSensPlane()    >= 0) {
    planBeg = fTRD->GetSensPlane();
    planEnd = planBeg + 1;
  }
  Int_t sectBeg = 0;
  Int_t sectEnd = AliTRDgeometry::Nsect();

  // 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 (fTRD->GetSensSector() >= 0) {
          Int_t sens1 = fTRD->GetSensSector();
          Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
          sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
                 * AliTRDgeometry::Nsect();
          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;
        Int_t nClustersUnfold = 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 nTimeBefore = geo->GetTimeBefore();
        Int_t nTimeTotal  = geo->GetTimeTotal();  

        // Get the digits
        digits = fDigitsManager->GetDigits(idet);
        digits->Expand();
        track0 = fDigitsManager->GetDictionary(idet,0);
        track0->Expand();
        track1 = fDigitsManager->GetDictionary(idet,1);
        track1->Expand();
        track2 = fDigitsManager->GetDictionary(idet,2); 
        track2->Expand();

        // Loop through the chamber and find the maxima 
        for ( row = 0;  row <  nRowMax;    row++) {
          for ( col = 2;  col <  nColMax;    col++) {
            for (time = 0; time < nTimeTotal; time++) {

              Int_t signalL = digits->GetDataUnchecked(row,col  ,time);
              Int_t signalM = digits->GetDataUnchecked(row,col-1,time);
              Int_t signalR = digits->GetDataUnchecked(row,col-2,time);
 
	      // Look for the maximum
              if ((signalM >= maxThresh) &&
                  (signalL >= sigThresh) &&
                  (signalR >= sigThresh)) {
                if ((signalL < signalM) && (signalR < signalM)) {
                  // Maximum found, mark the position by a negative signal
                  digits->SetDataUnchecked(row,col-1,time,-signalM);
		}
	      }

            }  
          }    
        }      

        // Now check the maxima and calculate the cluster position
        for ( row = 0;  row <  nRowMax  ;  row++) {
          for ( col = 1;  col <  nColMax-1;  col++) {
            for (time = 0; time < nTimeTotal; time++) {

              // Maximum found ?             
              if (digits->GetDataUnchecked(row,col,time) < 0) {

                Int_t iPad;
                for (iPad = 0; iPad < kNclus; iPad++) {
                  Int_t iPadCol = col - 1 + iPad;
                  clusterSignal[iPad]     = TMath::Abs(digits->GetDataUnchecked(row
                                                                               ,iPadCol
                                                                               ,time));
                  clusterDigit[iPad]      = digits->GetIndexUnchecked(row,iPadCol,time);
                  clusterTracks[3*iPad  ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
		  clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
		  clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
                }

                // Neighbouring maximum on right side?
                Int_t iType = 0; 
                if (col < nColMax-3) {
                  if (digits->GetDataUnchecked(row,col+2,time) < 0) {
                    for (iPad = 0; iPad < kNsig; iPad++) {
                      padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
                                                                           ,col-1+iPad
                                                                           ,time));
                    }
                    // Unfold the two maxima and set the signal on 
                    // the overlapping pad to the ratio
                    clusterSignal[2] *= Unfold(kEpsilon,padSignal);
                    iType = 1;
                  }
                }

                Float_t clusterCharge = clusterSignal[0]
                                      + clusterSignal[1]
                                      + clusterSignal[2];
                
		// Calculate the position of the cluster by using the
		// center of gravity method
                clusterPads[0] = row + 0.5;
                clusterPads[1] = col + 0.5 
                               + (clusterSignal[2] - clusterSignal[0]) 
		               / clusterCharge;
		// Take the shift of the additional time bins into account
                clusterPads[2] = time - nTimeBefore + 0.5;

                if (fVerbose) {
                  printf("-----------------------------------------------------------\n");
                  printf("Create cluster no. %d\n",nClusters);
                  printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
			                                            ,clusterPads[1]
                                                                    ,clusterPads[2]);
                  printf("Indices: %d, %d, %d\n",clusterDigit[0]
			                        ,clusterDigit[1]
                                                ,clusterDigit[2]);
                  printf("Total charge = %f\n",clusterCharge);
                  printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
			                            ,clusterTracks[1]
                                                    ,clusterTracks[2]);
                  printf("        pad1 %d, %d, %d\n",clusterTracks[3]
			                            ,clusterTracks[4]
                                                    ,clusterTracks[5]);
                  printf("        pad2 %d, %d, %d\n",clusterTracks[6]
			                            ,clusterTracks[7]
                                                    ,clusterTracks[8]);
                }

                nClusters++;
                if (iType == 1) nClustersUnfold++;

                // Add the cluster to the output array 
                fTRD->AddCluster(clusterPads
                                ,clusterDigit
                                ,idet
                                ,clusterCharge
                                ,clusterTracks
                                ,iType);

              }
            } 
          }   
        }     

	// Compress the arrays
        digits->Compress(1,0);
        track0->Compress(1,0);
        track1->Compress(1,0);
        track2->Compress(1,0);

        // Write the cluster and reset the array
	WriteClusters(idet);
	fTRD->ResetRecPoints();

        printf("AliTRDclusterizerV1::MakeCluster -- ");
        printf("Found %d (%d unfolded) clusters.\n"
              ,nClusters,nClustersUnfold);

      }    
    }      
  }        

  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 the 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.
  // Updated for new PRF 1/5/01.
  //

  // The parameters for the response function
  const Float_t kA  =  0.8303; 
  const Float_t kB  = -0.00392;
  const Float_t kC  =  0.472 * 0.472;
  const Float_t kD  =  2.19;

  Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));

  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$
3e1a3ad8	18	Revision 1.9 2000/11/01 14:53:20 cblume
	19	Merge with TRD-develop
	20
793ff80c	21	Revision 1.1.4.5 2000/10/15 23:40:01 cblume
	22	Remove AliTRDconst
	23
	24	Revision 1.1.4.4 2000/10/06 16:49:46 cblume
	25	Made Getters const
	26
	27	Revision 1.1.4.3 2000/10/04 16:34:58 cblume
	28	Replace include files by forward declarations
	29
	30	Revision 1.1.4.2 2000/09/22 14:49:49 cblume
	31	Adapted to tracking code
	32
	33	Revision 1.8 2000/10/02 21:28:19 fca
	34	Removal of useless dependecies via forward declarations
	35
94de3818	36	Revision 1.7 2000/06/27 13:08:50 cblume
	37	Changed to Copy(TObject &A) to appease the HP-compiler
	38
43da34c0	39	Revision 1.6 2000/06/09 11:10:07 cblume
	40	Compiler warnings and coding conventions, next round
	41
dd9a6ee3	42	Revision 1.5 2000/06/08 18:32:58 cblume
	43	Make code compliant to coding conventions
	44
8230f242	45	Revision 1.4 2000/06/07 16:27:01 cblume
	46	Try to remove compiler warnings on Sun and HP
	47
9d0b222b	48	Revision 1.3 2000/05/08 16:17:27 cblume
	49	Merge TRD-develop
	50
6f1e466d	51	Revision 1.1.4.1 2000/05/08 15:09:01 cblume
	52	Introduce AliTRDdigitsManager
	53
c0dd96c3	54	Revision 1.1 2000/02/28 18:58:54 cblume
	55	Add new TRD classes
	56
f7336fa3	57	*/
	58
	59	///////////////////////////////////////////////////////////////////////////////
	60	// //
	61	// TRD cluster finder for the slow simulator.
	62	// //
	63	///////////////////////////////////////////////////////////////////////////////
	64
	65	#include <TF1.h>
94de3818	66	#include <TTree.h>
793ff80c	67	#include <TH1.h>
f7336fa3	68
793ff80c	69	#include "AliRun.h"
	70
	71	#include "AliTRD.h"
f7336fa3	72	#include "AliTRDclusterizerV1.h"
	73	#include "AliTRDmatrix.h"
	74	#include "AliTRDgeometry.h"
	75	#include "AliTRDdigitizer.h"
6f1e466d	76	#include "AliTRDdataArrayF.h"
793ff80c	77	#include "AliTRDdataArrayI.h"
793ff80c	78	#include "AliTRDdigitsManager.h"
f7336fa3	79
	80	ClassImp(AliTRDclusterizerV1)
	81
	82	//_____________________________________________________________________________
	83	AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
	84	{
	85	//
	86	// AliTRDclusterizerV1 default constructor
	87	//
	88
6f1e466d	89	fDigitsManager = NULL;
f7336fa3	90
3e1a3ad8	91	fClusMaxThresh = 0;
	92	fClusSigThresh = 0;
	93
f7336fa3	94	}
	95
	96	//_____________________________________________________________________________
	97	AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
	98	:AliTRDclusterizer(name,title)
	99	{
	100	//
	101	// AliTRDclusterizerV1 default constructor
	102	//
	103
6f1e466d	104	fDigitsManager = new AliTRDdigitsManager();
f7336fa3	105
	106	Init();
	107
	108	}
	109
8230f242	110	//_____________________________________________________________________________
dd9a6ee3	111	AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
8230f242	112	{
	113	//
	114	// AliTRDclusterizerV1 copy constructor
	115	//
	116
dd9a6ee3	117	((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242	118
	119	}
	120
f7336fa3	121	//_____________________________________________________________________________
	122	AliTRDclusterizerV1::~AliTRDclusterizerV1()
	123	{
8230f242	124	//
	125	// AliTRDclusterizerV1 destructor
	126	//
f7336fa3	127
6f1e466d	128	if (fDigitsManager) {
6f1e466d	129	delete fDigitsManager;
f7336fa3	130	}
	131
	132	}
	133
dd9a6ee3	134	//_____________________________________________________________________________
	135	AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
	136	{
	137	//
	138	// Assignment operator
	139	//
	140
	141	if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
	142	return *this;
	143
	144	}
	145
8230f242	146	//_____________________________________________________________________________
43da34c0	147	void AliTRDclusterizerV1::Copy(TObject &c)
8230f242	148	{
	149	//
	150	// Copy function
	151	//
	152
43da34c0	153	((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
43da34c0	154	((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
43da34c0	155	((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
8230f242	156
	157	AliTRDclusterizer::Copy(c);
	158
	159	}
	160
f7336fa3	161	//_____________________________________________________________________________
	162	void AliTRDclusterizerV1::Init()
	163	{
	164	//
	165	// Initializes the cluster finder
	166	//
	167
	168	// The default parameter for the clustering
3e1a3ad8	169	fClusMaxThresh = 3;
3e1a3ad8	170	fClusSigThresh = 1;
f7336fa3	171
	172	}
	173
	174	//_____________________________________________________________________________
	175	Bool_t AliTRDclusterizerV1::ReadDigits()
	176	{
	177	//
	178	// Reads the digits arrays from the input aliroot file
	179	//
	180
	181	if (!fInputFile) {
	182	printf("AliTRDclusterizerV1::ReadDigits -- ");
	183	printf("No input file open\n");
	184	return kFALSE;
	185	}
	186
f7336fa3	187	// Read in the digit arrays
6f1e466d	188	return (fDigitsManager->ReadDigits());
f7336fa3	189
	190	}
	191
	192	//_____________________________________________________________________________
793ff80c	193	Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3	194	{
	195	//
	196	// Generates the cluster.
	197	//
	198
	199	Int_t row, col, time;
	200
3e1a3ad8	201	if (fTRD->IsVersion() != 1) {
f7336fa3	202	printf("AliTRDclusterizerV1::MakeCluster -- ");
	203	printf("TRD must be version 1 (slow simulator).\n");
	204	return kFALSE;
	205	}
	206
	207	// Get the geometry
3e1a3ad8	208	AliTRDgeometry *geo = fTRD->GetGeometry();
f7336fa3	209
	210	printf("AliTRDclusterizerV1::MakeCluster -- ");
	211	printf("Start creating clusters.\n");
	212
8230f242	213	AliTRDdataArrayI *digits;
793ff80c	214	AliTRDdataArrayI *track0;
	215	AliTRDdataArrayI *track1;
	216	AliTRDdataArrayI *track2;
f7336fa3	217
3e1a3ad8	218	// Threshold value for the maximum
	219	Int_t maxThresh = fClusMaxThresh;
	220	// Threshold value for the digit signal
	221	Int_t sigThresh = fClusSigThresh;
f7336fa3	222
f7336fa3	223	// Iteration limit for unfolding procedure
8230f242	224	const Float_t kEpsilon = 0.01;
f7336fa3	225
8230f242	226	const Int_t kNclus = 3;
8230f242	227	const Int_t kNsig = 5;
3e1a3ad8	228	const Int_t kNtrack = 3 * kNclus;
	229
	230	Float_t padSignal[kNsig];
	231	Float_t clusterSignal[kNclus];
	232	Float_t clusterPads[kNclus];
	233	Int_t clusterDigit[kNclus];
	234	Int_t clusterTracks[kNtrack];
f7336fa3	235
f7336fa3	236	Int_t chamBeg = 0;
793ff80c	237	Int_t chamEnd = AliTRDgeometry::Ncham();
3e1a3ad8	238	if (fTRD->GetSensChamber() >= 0) {
3e1a3ad8	239	chamBeg = fTRD->GetSensChamber();
6f1e466d	240	chamEnd = chamBeg + 1;
f7336fa3	241	}
f7336fa3	242	Int_t planBeg = 0;
793ff80c	243	Int_t planEnd = AliTRDgeometry::Nplan();
3e1a3ad8	244	if (fTRD->GetSensPlane() >= 0) {
3e1a3ad8	245	planBeg = fTRD->GetSensPlane();
f7336fa3	246	planEnd = planBeg + 1;
	247	}
	248	Int_t sectBeg = 0;
793ff80c	249	Int_t sectEnd = AliTRDgeometry::Nsect();
f7336fa3	250
3e1a3ad8	251	// Start clustering in every chamber
f7336fa3	252	for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
	253	for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
	254	for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
	255
3e1a3ad8	256	if (fTRD->GetSensSector() >= 0) {
	257	Int_t sens1 = fTRD->GetSensSector();
	258	Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
793ff80c	259	sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
793ff80c	260	* AliTRDgeometry::Nsect();
dd9a6ee3	261	if (sens1 < sens2) {
9d0b222b	262	if ((isect < sens1) \|\| (isect >= sens2)) continue;
dd9a6ee3	263	}
dd9a6ee3	264	else {
9d0b222b	265	if ((isect < sens1) && (isect >= sens2)) continue;
dd9a6ee3	266	}
9d0b222b	267	}
9d0b222b	268
8230f242	269	Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3	270
3e1a3ad8	271	Int_t nClusters = 0;
	272	Int_t nClustersUnfold = 0;
	273
f7336fa3	274	printf("AliTRDclusterizerV1::MakeCluster -- ");
f7336fa3	275	printf("Analyzing chamber %d, plane %d, sector %d.\n"
3e1a3ad8	276	,icham,iplan,isect);
f7336fa3	277
3e1a3ad8	278	Int_t nRowMax = geo->GetRowMax(iplan,icham,isect);
	279	Int_t nColMax = geo->GetColMax(iplan);
	280	Int_t nTimeBefore = geo->GetTimeBefore();
	281	Int_t nTimeTotal = geo->GetTimeTotal();
f7336fa3	282
3e1a3ad8	283	// Get the digits
8230f242	284	digits = fDigitsManager->GetDigits(idet);
3e1a3ad8	285	digits->Expand();
793ff80c	286	track0 = fDigitsManager->GetDictionary(idet,0);
3e1a3ad8	287	track0->Expand();
793ff80c	288	track1 = fDigitsManager->GetDictionary(idet,1);
3e1a3ad8	289	track1->Expand();
793ff80c	290	track2 = fDigitsManager->GetDictionary(idet,2);
3e1a3ad8	291	track2->Expand();
	292
	293	// Loop through the chamber and find the maxima
	294	for ( row = 0; row < nRowMax; row++) {
	295	for ( col = 2; col < nColMax; col++) {
	296	for (time = 0; time < nTimeTotal; time++) {
	297
	298	Int_t signalL = digits->GetDataUnchecked(row,col ,time);
	299	Int_t signalM = digits->GetDataUnchecked(row,col-1,time);
	300	Int_t signalR = digits->GetDataUnchecked(row,col-2,time);
	301
	302	// Look for the maximum
	303	if ((signalM >= maxThresh) &&
	304	(signalL >= sigThresh) &&
	305	(signalR >= sigThresh)) {
	306	if ((signalL < signalM) && (signalR < signalM)) {
	307	// Maximum found, mark the position by a negative signal
	308	digits->SetDataUnchecked(row,col-1,time,-signalM);
	309	}
	310	}
	311
	312	}
	313	}
	314	}
	315
	316	// Now check the maxima and calculate the cluster position
	317	for ( row = 0; row < nRowMax ; row++) {
	318	for ( col = 1; col < nColMax-1; col++) {
	319	for (time = 0; time < nTimeTotal; time++) {
	320
	321	// Maximum found ?
	322	if (digits->GetDataUnchecked(row,col,time) < 0) {
f7336fa3	323
9d0b222b	324	Int_t iPad;
8230f242	325	for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8	326	Int_t iPadCol = col - 1 + iPad;
	327	clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
	328	,iPadCol
	329	,time));
	330	clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
	331	clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
	332	clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
	333	clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
f7336fa3	334	}
f7336fa3	335
3e1a3ad8	336	// Neighbouring maximum on right side?
	337	Int_t iType = 0;
	338	if (col < nColMax-3) {
	339	if (digits->GetDataUnchecked(row,col+2,time) < 0) {
	340	for (iPad = 0; iPad < kNsig; iPad++) {
	341	padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
	342	,col-1+iPad
	343	,time));
f7336fa3	344	}
3e1a3ad8	345	// Unfold the two maxima and set the signal on
	346	// the overlapping pad to the ratio
	347	clusterSignal[2] *= Unfold(kEpsilon,padSignal);
	348	iType = 1;
f7336fa3	349	}
f7336fa3	350	}
f7336fa3	351
3e1a3ad8	352	Float_t clusterCharge = clusterSignal[0]
	353	+ clusterSignal[1]
	354	+ clusterSignal[2];
	355
	356	// Calculate the position of the cluster by using the
	357	// center of gravity method
	358	clusterPads[0] = row + 0.5;
	359	clusterPads[1] = col + 0.5
	360	+ (clusterSignal[2] - clusterSignal[0])
	361	/ clusterCharge;
	362	// Take the shift of the additional time bins into account
	363	clusterPads[2] = time - nTimeBefore + 0.5;
	364
	365	if (fVerbose) {
	366	printf("-----------------------------------------------------------\n");
	367	printf("Create cluster no. %d\n",nClusters);
	368	printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
	369	,clusterPads[1]
	370	,clusterPads[2]);
	371	printf("Indices: %d, %d, %d\n",clusterDigit[0]
	372	,clusterDigit[1]
	373	,clusterDigit[2]);
	374	printf("Total charge = %f\n",clusterCharge);
	375	printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
	376	,clusterTracks[1]
	377	,clusterTracks[2]);
	378	printf(" pad1 %d, %d, %d\n",clusterTracks[3]
	379	,clusterTracks[4]
	380	,clusterTracks[5]);
	381	printf(" pad2 %d, %d, %d\n",clusterTracks[6]
	382	,clusterTracks[7]
	383	,clusterTracks[8]);
f7336fa3	384	}
f7336fa3	385
3e1a3ad8	386	nClusters++;
3e1a3ad8	387	if (iType == 1) nClustersUnfold++;
f7336fa3	388
f7336fa3	389	// Add the cluster to the output array
3e1a3ad8	390	fTRD->AddCluster(clusterPads
	391	,clusterDigit
	392	,idet
	393	,clusterCharge
	394	,clusterTracks
	395	,iType);
f7336fa3	396
f7336fa3	397	}
3e1a3ad8	398	}
	399	}
	400	}
f7336fa3	401
3e1a3ad8	402	// Compress the arrays
	403	digits->Compress(1,0);
	404	track0->Compress(1,0);
	405	track1->Compress(1,0);
	406	track2->Compress(1,0);
f7336fa3	407
3e1a3ad8	408	// Write the cluster and reset the array
793ff80c	409	WriteClusters(idet);
3e1a3ad8	410	fTRD->ResetRecPoints();
793ff80c	411
3e1a3ad8	412	printf("AliTRDclusterizerV1::MakeCluster -- ");
	413	printf("Found %d (%d unfolded) clusters.\n"
	414	,nClusters,nClustersUnfold);
f7336fa3	415
3e1a3ad8	416	}
	417	}
	418	}
f7336fa3	419
f7336fa3	420	printf("AliTRDclusterizerV1::MakeCluster -- ");
	421	printf("Done.\n");
	422
	423	return kTRUE;
	424
	425	}
	426
	427	//_____________________________________________________________________________
	428	Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
	429	{
	430	//
	431	// Method to unfold neighbouring maxima.
	432	// The charge ratio on the overlapping pad is calculated
	433	// until there is no more change within the range given by eps.
	434	// The resulting ratio is then returned to the calling method.
	435	//
	436
3e1a3ad8	437	Int_t itStep = 0; // Count iteration steps
f7336fa3	438
3e1a3ad8	439	Float_t ratio = 0.5; // Start value for ratio
3e1a3ad8	440	Float_t prevRatio = 0; // Store previous ratio
f7336fa3	441
3e1a3ad8	442	Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
3e1a3ad8	443	Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
f7336fa3	444
3e1a3ad8	445	// Start the iteration
f7336fa3	446	while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
	447
	448	itStep++;
	449	prevRatio = ratio;
	450
3e1a3ad8	451	// Cluster position according to charge ratio
	452	Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
	453	/ (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
	454	Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
	455	/ ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3	456
3e1a3ad8	457	// Set cluster charge ratio
f7336fa3	458	Float_t ampLeft = padSignal[1];
	459	Float_t ampRight = padSignal[3];
	460
3e1a3ad8	461	// Apply pad response to parameters
	462	newLeftSignal[0] = ampLeft * PadResponse(-1 - maxLeft);
	463	newLeftSignal[1] = ampLeft * PadResponse( 0 - maxLeft);
	464	newLeftSignal[2] = ampLeft * PadResponse( 1 - maxLeft);
f7336fa3	465
3e1a3ad8	466	newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
	467	newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
	468	newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
f7336fa3	469
3e1a3ad8	470	// Calculate new overlapping ratio
3e1a3ad8	471	ratio = newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]);
f7336fa3	472
	473	}
	474
	475	return ratio;
	476
	477	}
	478
	479	//_____________________________________________________________________________
	480	Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
	481	{
	482	//
	483	// The pad response for the chevron pads.
	484	// We use a simple Gaussian approximation which should be good
	485	// enough for our purpose.
3e1a3ad8	486	// Updated for new PRF 1/5/01.
f7336fa3	487	//
	488
	489	// The parameters for the response function
3e1a3ad8	490	const Float_t kA = 0.8303;
	491	const Float_t kB = -0.00392;
	492	const Float_t kC = 0.472 * 0.472;
	493	const Float_t kD = 2.19;
f7336fa3	494
3e1a3ad8	495	Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(xx,kD) / (2.kC)));
f7336fa3	496
	497	return (pr);
	498
	499	}