[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.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 "AliTRDrecPoint.h"
#include "AliTRDdataArrayF.h"
#include "AliTRDdataArrayI.h"
#include "AliTRDdigitsManager.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(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).fClusMethod    = fClusMethod;
  ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;

  AliTRDclusterizer::Copy(c);

}

//_____________________________________________________________________________
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::MakeClusters()
{
  //
  // 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;
  AliTRDdataArrayI *track0;
  AliTRDdataArrayI *track1;
  AliTRDdataArrayI *track2; 

  // 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 kEpsilon = 0.01;             

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

  Int_t chamBeg = 0;
  Int_t chamEnd = AliTRDgeometry::Ncham();
  if (trd->GetSensChamber()  >= 0) {
    chamBeg = trd->GetSensChamber();
    chamEnd = chamBeg + 1;
  }
  Int_t planBeg = 0;
  Int_t planEnd = AliTRDgeometry::Nplan();
  if (trd->GetSensPlane()    >= 0) {
    planBeg = trd->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 (trd->GetSensSector() >= 0) {
          Int_t sens1 = trd->GetSensSector();
          Int_t sens2 = sens1 + trd->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;
        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);

        // Create a matrix for track indexes
        AliTRDmatrix *trackMatrix  = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
                                                     ,isect,icham,iplan);

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

        // 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);
              Int_t t[3] = {-1};
              t[0] = track0->GetData(row,col,time) - 1;
              t[1] = track1->GetData(row,col,time) - 1;
              t[2] = track2->GetData(row,col,time) - 1;  

              // 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);
                for(Int_t i = 0; i < 3; i++) {
                  trackMatrix->AddTrack(row,col,time,t[i]);
                }      
              }
	    }
	  }
	}

        // 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[kNsig]      = {0};   
                // Signals from cluster
                Float_t clusterSignal[kNclus] = {0};
                // Cluster pad info
                Float_t clusterPads[kNclus]   = {0};   
                // Cluster digit info
                Int_t   clusterDigit[kNclus]  = {0};
                // Cluster MC tracks info
                const Int_t nt = kNclus*3;
                Int_t clusterTracks[nt] = {-1};   

                Int_t iPad;
                for (iPad = 0; iPad < kNclus; iPad++) {
                  clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
                  clusterDigit[iPad]  = digitMatrix->GetTrack(row,col-1+iPad,time,0);
                  for (Int_t j = 0; j < 3; j++) {
                    clusterTracks[iPad*3+j] = trackMatrix->GetTrack(row,col-1+iPad,time,j);
                  }  
                }

                // 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(kEpsilon, 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;


/*		  printf("col = %d, left = %f, center = %f, right = %f, 
                         final =%f \n", col, 
                         digitMatrix->GetSignal(row,col-1,time),
                         digitMatrix->GetSignal(row,col,time),
                         digitMatrix->GetSignal(row,col+1,time),
                         clusterPads[1]);

		  printf("col = %d, sig(0) = %f, sig(1) = %f, sig(2) = %f, 
                         final =%f \n", col, 
                         clusterSignal[0],
                         clusterSignal[1],
                         clusterSignal[2],
                         clusterPads[1]);

*/

                  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,clusterTracks);

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

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

	WriteClusters(idet);
	trd->ResetRecPoints();

        delete digitMatrix;
        delete maximaMatrix;
        delete trackMatrix;

      }          // 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 kA  =  0.8872;
  const Float_t kB  = -0.00573;
  const Float_t kC  =  0.454;
  const Float_t kC2 =  kC*kC;

  Float_t pr = kA * (kB + TMath::Exp(-x*x / (2. * kC2)));

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