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

/* $Id$ */

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

#include <TF1.h>
#include <TTree.h>
#include <TH1.h>
#include <TFile.h>

#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliLoader.h"
#include "AliRawReader.h"

#include "AliTRDclusterizerV1.h"
#include "AliTRDmatrix.h"
#include "AliTRDgeometry.h"
#include "AliTRDdataArrayF.h"
#include "AliTRDdataArrayI.h"
#include "AliTRDdigitsManager.h"
#include "AliTRDparameter.h"
#include "AliTRDpadPlane.h"
#include "AliTRDrawData.h"

ClassImp(AliTRDclusterizerV1)

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

  fDigitsManager = 0;

}

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

  fDigitsManager = new AliTRDdigitsManager();
  fDigitsManager->CreateArrays();

}

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

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

}

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

  if (fDigitsManager) {
    delete fDigitsManager;
    fDigitsManager = NULL;
  }

}

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

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

}

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

  ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;

  AliTRDclusterizer::Copy(c);

}

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

  if (!fRunLoader) {
    printf("<AliTRDclusterizerV1::ReadDigits> ");
    printf("No input file open\n");
    return kFALSE;
  }
  AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
  if (!loader->TreeD()) loader->LoadDigits();

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

}

//_____________________________________________________________________________
Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader* rawReader)
{
  //
  // Reads the digits arrays from the ddl file
  //

  AliTRDrawData *raw = new AliTRDrawData();
  raw->SetDebug(1);

  fDigitsManager = raw->Raw2Digits(rawReader);

  return kTRUE;

}

//_____________________________________________________________________________
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 = AliTRDgeometry::GetGeometry(fRunLoader);

  // Create a default parameter class if none is defined
  if (!fPar) {
    fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
    printf("<AliTRDclusterizerV1::MakeCluster> ");
    printf("Create the default parameter object.\n");
  }
  fPar->Init();

  //Float_t timeBinSize = fPar->GetDriftVelocity()
  //                    / fPar->GetSamplingFrequency();
  // Half of ampl.region
  //  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 

  Float_t omegaTau = fPar->GetOmegaTau();
  if (fVerbose > 0) {
    printf("<AliTRDclusterizerV1::MakeCluster> ");
    printf("OmegaTau = %f \n",omegaTau);
    printf("<AliTRDclusterizerV1::MakeCluster> ");
    printf("Start creating clusters.\n");
  } 

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

  // Threshold value for the maximum
  Int_t maxThresh = fPar->GetClusMaxThresh();   
  // Threshold value for the digit signal
  Int_t sigThresh = fPar->GetClusSigThresh();   
  // 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;

  Int_t    iType         = 0;
  Int_t    iUnfold       = 0;  
  Double_t ratioLeft     = 1.0;
  Double_t ratioRight    = 1.0;

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

  Int_t    chamBeg = 0;
  Int_t    chamEnd = AliTRDgeometry::Ncham();
  Int_t    planBeg = 0;
  Int_t    planEnd = AliTRDgeometry::Nplan();
  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++) {

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

        Int_t nClusters      = 0;
        Int_t nClusters2pad  = 0;
        Int_t nClusters3pad  = 0;
        Int_t nClusters4pad  = 0;
        Int_t nClusters5pad  = 0;
        Int_t nClustersLarge = 0;

        if (fVerbose > 0) {
          printf("<AliTRDclusterizerV1::MakeCluster> ");
          printf("Analyzing chamber %d, plane %d, sector %d.\n"
                ,icham,iplan,isect);
	}

        Int_t    nRowMax     = fPar->GetRowMax(iplan,icham,isect);
        Int_t    nColMax     = fPar->GetColMax(iplan);
        Int_t    nTimeBefore = fPar->GetTimeBefore();
        Int_t    nTimeTotal  = fPar->GetTimeTotal();  

        AliTRDpadPlane *padPlane = fPar->GetPadPlane(iplan,icham);

        // 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 ( col = 4;  col <  nColMax-2;    col++) {
            for (time = 0; time < nTimeTotal; time++) {

              Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col  ,time));
              Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
              Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
 
// 	      // Look for the maximum
//               if (signalM >= maxThresh) {
//                 if (((signalL >= sigThresh) &&
//                      (signalL <  signalM))  ||
//                     ((signalR >= sigThresh) &&
//                      (signalR <  signalM))) {
//                   // Maximum found, mark the position by a negative signal
//                   digits->SetDataUnchecked(row,col-1,time,-signalM);
// 		}
// 	      }
	      // Look for the maximum
              if (signalM >= maxThresh) {
                if ( (TMath::Abs(signalL)<=signalM) && (TMath::Abs(signalR)<=signalM) && 
		     (TMath::Abs(signalL)+TMath::Abs(signalR))>sigThresh ) {
                  // 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 (time = 0; time < nTimeTotal; time++) {
            for ( col = 1;  col <  nColMax-1;  col++) {

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

		// Count the number of pads in the cluster
                Int_t nPadCount = 0;
                Int_t ii        = 0;
                while (TMath::Abs(digits->GetDataUnchecked(row,col-ii  ,time))
                                                                  >= sigThresh) {
                  nPadCount++;
                  ii++;
                  if (col-ii   <        0) break;
		}
                ii = 0;
                while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
                                                                  >= sigThresh) {
                  nPadCount++;
                  ii++;
                  if (col+ii+1 >= nColMax) break;
		}

                nClusters++;
                switch (nPadCount) {
                case 2:
                  iType = 0;
                  nClusters2pad++;
                  break;
                case 3:
                  iType = 1;
                  nClusters3pad++;
                  break;
                case 4:
                  iType = 2;
                  nClusters4pad++;
                  break;
                case 5:
                  iType = 3;
                  nClusters5pad++;
                  break;
                default:
                  iType = 4;
                  nClustersLarge++;
                  break;
		};

		 // Look for 5 pad cluster with minimum in the middle
                Bool_t fivePadCluster = kFALSE;
                if (col < nColMax-3) {
                  if (digits->GetDataUnchecked(row,col+2,time) < 0) {
                    fivePadCluster = kTRUE;
	  	  }
                  if ((fivePadCluster) && (col < nColMax-5)) {
                    if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
                      fivePadCluster = kFALSE;
		    }
		  }
                  if ((fivePadCluster) && (col >         1)) {
                    if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
                      fivePadCluster = kFALSE;
		    }
		  }
		}

		// 5 pad cluster
                // Modify the signal of the overlapping pad for the left part 
		// of the cluster which remains from a previous unfolding
                if (iUnfold) {
                  clusterSignal[0] *= ratioLeft;
                  iType   = 5;
                  iUnfold = 0;
		}

		// Unfold the 5 pad cluster
                if (fivePadCluster) {
                  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
                  ratioRight        = Unfold(kEpsilon,iplan,padSignal);
                  ratioLeft         = 1.0 - ratioRight; 
                  clusterSignal[2] *= ratioRight;
                  iType   = 5;
                  iUnfold = 1;
                }


                Double_t clusterCharge = clusterSignal[0]
                                       + clusterSignal[1]
                                       + clusterSignal[2];
                
		// The position of the cluster
                clusterPads[0] = row + 0.5;
		// Take the shift of the additional time bins into account
                clusterPads[2] = time - nTimeBefore + 0.5;

                if (fPar->LUTOn()) {
  		  // Calculate the position of the cluster by using the
		  // lookup table method
                  clusterPads[1] =
                                  fPar->LUTposition(iplan,clusterSignal[0]
                                                          ,clusterSignal[1]
					                  ,clusterSignal[2]);
		}
		else {
  		  // Calculate the position of the cluster by using the
		  // center of gravity method
		  for (Int_t i=0;i<5;i++) padSignal[i]=0;
		  padSignal[2] = TMath::Abs(digits->GetDataUnchecked(row,col,time));   // central  pad
		  padSignal[1] = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); // left     pad
		  padSignal[3] = TMath::Abs(digits->GetDataUnchecked(row,col+1,time)); // right    pad
		  if (col>2 &&TMath::Abs(digits->GetDataUnchecked(row,col-2,time)<padSignal[1])){
		    padSignal[0] = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
		  }
		  if (col<nColMax-3 &&TMath::Abs(digits->GetDataUnchecked(row,col+2,time)<padSignal[3])){
		    padSignal[4] = TMath::Abs(digits->GetDataUnchecked(row,col+2,time));
		  }		  
		  clusterPads[1] =  GetCOG(padSignal);
		  Double_t check = fPar->LUTposition(iplan,clusterSignal[0]
                                                          ,clusterSignal[1]
					                  ,clusterSignal[2]);
		  //		  Float_t diff = clusterPads[1] -  check;

		}

                Double_t q0 = clusterSignal[0];
		Double_t q1 = clusterSignal[1];
                Double_t q2 = clusterSignal[2];
                Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
                                          (clusterCharge*clusterCharge);

		// Calculate the position and the error
                Double_t colSize = padPlane->GetColSize(col);
                Double_t rowSize = padPlane->GetRowSize(row);
                Double_t clusterPos[3];
		clusterPos[0] = padPlane->GetColPos(col) + (clusterPads[1]-0.5)*colSize;  // MI change
		clusterPos[1] = padPlane->GetRowPos(row) -0.5*rowSize; //MI change
                clusterPos[2] = clusterPads[2];
                Double_t clusterSig[2];
                clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
                clusterSig[1] = rowSize * rowSize / 12.;                                       
                // Add the cluster to the output array 
                AddCluster(clusterPos
			  ,idet
			  ,clusterCharge
			  ,clusterTracks
			  ,clusterSig
			   ,iType,clusterPads[1]);

              }
            } 
          }   
        }     

	// 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);
	ResetRecPoints();
      }    
    }      
  }        

  if (fVerbose > 0) {
    printf("<AliTRDclusterizerV1::MakeCluster> ");
    printf("Done.\n");
  }

  return kTRUE;

}

Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
{
  //
  // get COG position
  // used for clusters with more than 3 pads - where LUT not applicable
  Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4];
  Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum;
  return res;		  
}


//_____________________________________________________________________________
Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_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   irc                = 0;
  Int_t   itStep             = 0;      // Count iteration steps

  Double_t ratio             = 0.5;    // Start value for ratio
  Double_t prevRatio         = 0;      // Store previous ratio

  Double_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
  Double_t newRightSignal[3] = {0};    // Array to store right cluster signal
  Double_t newSignal[3]      = {0};

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

    itStep++;
    prevRatio = ratio;

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

    // Set cluster charge ratio
    irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
    Double_t ampLeft  = padSignal[1] / newSignal[1];
    irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
    Double_t ampRight = padSignal[3] / newSignal[1];

    // Apply pad response to parameters
    irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
    irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);

    // Calculate new overlapping ratio
    ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / 
                          (newLeftSignal[2] + newRightSignal[0]));

  }

  return ratio;

}
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
88cb7938	16	/* $Id$ */
f7336fa3	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// TRD cluster finder for the slow simulator.
	21	// //
	22	///////////////////////////////////////////////////////////////////////////////
	23
	24	#include <TF1.h>
94de3818	25	#include <TTree.h>
793ff80c	26	#include <TH1.h>
a819a5f7	27	#include <TFile.h>
f7336fa3	28
793ff80c	29	#include "AliRun.h"
88cb7938	30	#include "AliRunLoader.h"
88cb7938	31	#include "AliLoader.h"
928e9fae	32	#include "AliRawReader.h"
793ff80c	33
f7336fa3	34	#include "AliTRDclusterizerV1.h"
	35	#include "AliTRDmatrix.h"
	36	#include "AliTRDgeometry.h"
6f1e466d	37	#include "AliTRDdataArrayF.h"
793ff80c	38	#include "AliTRDdataArrayI.h"
793ff80c	39	#include "AliTRDdigitsManager.h"
17b26de4	40	#include "AliTRDparameter.h"
a5cadd36	41	#include "AliTRDpadPlane.h"
928e9fae	42	#include "AliTRDrawData.h"
f7336fa3	43
	44	ClassImp(AliTRDclusterizerV1)
	45
	46	//_____________________________________________________________________________
	47	AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
	48	{
	49	//
	50	// AliTRDclusterizerV1 default constructor
	51	//
	52
17b26de4	53	fDigitsManager = 0;
db30bf0f	54
f7336fa3	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();
17b26de4	66	fDigitsManager->CreateArrays();
f7336fa3	67
	68	}
	69
8230f242	70	//_____________________________________________________________________________
dd9a6ee3	71	AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
73ae7b59	72	:AliTRDclusterizer(c)
8230f242	73	{
	74	//
	75	// AliTRDclusterizerV1 copy constructor
	76	//
	77
dd9a6ee3	78	((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242	79
	80	}
	81
f7336fa3	82	//_____________________________________________________________________________
	83	AliTRDclusterizerV1::~AliTRDclusterizerV1()
	84	{
8230f242	85	//
	86	// AliTRDclusterizerV1 destructor
	87	//
f7336fa3	88
6f1e466d	89	if (fDigitsManager) {
6f1e466d	90	delete fDigitsManager;
abaf1f1d	91	fDigitsManager = NULL;
f7336fa3	92	}
	93
	94	}
	95
dd9a6ee3	96	//_____________________________________________________________________________
	97	AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
	98	{
	99	//
	100	// Assignment operator
	101	//
	102
	103	if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
	104	return *this;
	105
	106	}
	107
8230f242	108	//_____________________________________________________________________________
e0d47c25	109	void AliTRDclusterizerV1::Copy(TObject &c) const
8230f242	110	{
	111	//
	112	// Copy function
	113	//
	114
17b26de4	115	((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
8230f242	116
	117	AliTRDclusterizer::Copy(c);
	118
	119	}
	120
f7336fa3	121	//_____________________________________________________________________________
	122	Bool_t AliTRDclusterizerV1::ReadDigits()
	123	{
	124	//
	125	// Reads the digits arrays from the input aliroot file
	126	//
	127
88cb7938	128	if (!fRunLoader) {
17b26de4	129	printf("<AliTRDclusterizerV1::ReadDigits> ");
f7336fa3	130	printf("No input file open\n");
	131	return kFALSE;
	132	}
88cb7938	133	AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
88cb7938	134	if (!loader->TreeD()) loader->LoadDigits();
abaf1f1d	135
f7336fa3	136	// Read in the digit arrays
88cb7938	137	return (fDigitsManager->ReadDigits(loader->TreeD()));
f7336fa3	138
	139	}
	140
928e9fae	141	//_____________________________________________________________________________
	142	Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader* rawReader)
	143	{
	144	//
	145	// Reads the digits arrays from the ddl file
	146	//
	147
	148	AliTRDrawData *raw = new AliTRDrawData();
	149	raw->SetDebug(1);
	150
	151	fDigitsManager = raw->Raw2Digits(rawReader);
	152
	153	return kTRUE;
	154
	155	}
	156
f7336fa3	157	//_____________________________________________________________________________
793ff80c	158	Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3	159	{
	160	//
	161	// Generates the cluster.
	162	//
	163
	164	Int_t row, col, time;
	165
bdbb05bb	166	/*
3e1a3ad8	167	if (fTRD->IsVersion() != 1) {
17b26de4	168	printf("<AliTRDclusterizerV1::MakeCluster> ");
f7336fa3	169	printf("TRD must be version 1 (slow simulator).\n");
	170	return kFALSE;
	171	}
bdbb05bb	172	*/
f7336fa3	173
f7336fa3	174	// Get the geometry
bdbb05bb	175	AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
f7336fa3	176
17b26de4	177	// Create a default parameter class if none is defined
	178	if (!fPar) {
	179	fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
5443e65e	180	printf("<AliTRDclusterizerV1::MakeCluster> ");
5443e65e	181	printf("Create the default parameter object.\n");
17b26de4	182	}
598156ef	183	fPar->Init();
17b26de4	184
7ad19338	185	//Float_t timeBinSize = fPar->GetDriftVelocity()
7ad19338	186	// / fPar->GetSamplingFrequency();
a819a5f7	187	// Half of ampl.region
7ad19338	188	// const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
a819a5f7	189
17b26de4	190	Float_t omegaTau = fPar->GetOmegaTau();
47517f42	191	if (fVerbose > 0) {
17b26de4	192	printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42	193	printf("OmegaTau = %f \n",omegaTau);
17b26de4	194	printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42	195	printf("Start creating clusters.\n");
47517f42	196	}
f7336fa3	197
8230f242	198	AliTRDdataArrayI *digits;
793ff80c	199	AliTRDdataArrayI *track0;
	200	AliTRDdataArrayI *track1;
	201	AliTRDdataArrayI *track2;
f7336fa3	202
3e1a3ad8	203	// Threshold value for the maximum
17b26de4	204	Int_t maxThresh = fPar->GetClusMaxThresh();
3e1a3ad8	205	// Threshold value for the digit signal
17b26de4	206	Int_t sigThresh = fPar->GetClusSigThresh();
f7336fa3	207	// Iteration limit for unfolding procedure
8230f242	208	const Float_t kEpsilon = 0.01;
f7336fa3	209
8230f242	210	const Int_t kNclus = 3;
8230f242	211	const Int_t kNsig = 5;
3e1a3ad8	212	const Int_t kNtrack = 3 * kNclus;
3e1a3ad8	213
a5cadd36	214	Int_t iType = 0;
7ad19338	215	Int_t iUnfold = 0;
a5cadd36	216	Double_t ratioLeft = 1.0;
a5cadd36	217	Double_t ratioRight = 1.0;
db30bf0f	218
7ad19338	219	//
a5cadd36	220	Double_t padSignal[kNsig];
	221	Double_t clusterSignal[kNclus];
	222	Double_t clusterPads[kNclus];
	223	Int_t clusterDigit[kNclus];
	224	Int_t clusterTracks[kNtrack];
f7336fa3	225
a5cadd36	226	Int_t chamBeg = 0;
	227	Int_t chamEnd = AliTRDgeometry::Ncham();
	228	Int_t planBeg = 0;
	229	Int_t planEnd = AliTRDgeometry::Nplan();
	230	Int_t sectBeg = 0;
	231	Int_t sectEnd = AliTRDgeometry::Nsect();
f7336fa3	232
3e1a3ad8	233	// Start clustering in every chamber
f7336fa3	234	for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
	235	for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
	236	for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
	237
8230f242	238	Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3	239
db30bf0f	240	Int_t nClusters = 0;
	241	Int_t nClusters2pad = 0;
	242	Int_t nClusters3pad = 0;
	243	Int_t nClusters4pad = 0;
	244	Int_t nClusters5pad = 0;
	245	Int_t nClustersLarge = 0;
3e1a3ad8	246
47517f42	247	if (fVerbose > 0) {
17b26de4	248	printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42	249	printf("Analyzing chamber %d, plane %d, sector %d.\n"
	250	,icham,iplan,isect);
	251	}
f7336fa3	252
a5cadd36	253	Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect);
	254	Int_t nColMax = fPar->GetColMax(iplan);
	255	Int_t nTimeBefore = fPar->GetTimeBefore();
	256	Int_t nTimeTotal = fPar->GetTimeTotal();
5443e65e	257
a5cadd36	258	AliTRDpadPlane *padPlane = fPar->GetPadPlane(iplan,icham);
f7336fa3	259
3e1a3ad8	260	// Get the digits
8230f242	261	digits = fDigitsManager->GetDigits(idet);
3e1a3ad8	262	digits->Expand();
793ff80c	263	track0 = fDigitsManager->GetDictionary(idet,0);
3e1a3ad8	264	track0->Expand();
793ff80c	265	track1 = fDigitsManager->GetDictionary(idet,1);
3e1a3ad8	266	track1->Expand();
793ff80c	267	track2 = fDigitsManager->GetDictionary(idet,2);
3e1a3ad8	268	track2->Expand();
	269
	270	// Loop through the chamber and find the maxima
	271	for ( row = 0; row < nRowMax; row++) {
de4b10e5	272	for ( col = 2; col < nColMax; col++) {
de4b10e5	273	//for ( col = 4; col < nColMax-2; col++) {
3e1a3ad8	274	for (time = 0; time < nTimeTotal; time++) {
3e1a3ad8	275
a819a5f7	276	Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
	277	Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
	278	Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
3e1a3ad8	279
7ad19338	280	// // Look for the maximum
	281	// if (signalM >= maxThresh) {
	282	// if (((signalL >= sigThresh) &&
	283	// (signalL < signalM)) \|\|
	284	// ((signalR >= sigThresh) &&
	285	// (signalR < signalM))) {
	286	// // Maximum found, mark the position by a negative signal
	287	// digits->SetDataUnchecked(row,col-1,time,-signalM);
	288	// }
	289	// }
3e1a3ad8	290	// Look for the maximum
db30bf0f	291	if (signalM >= maxThresh) {
de4b10e5	292	if ( (TMath::Abs(signalL)<=signalM) && (TMath::Abs(signalR)<=signalM) &&
de4b10e5	293	(TMath::Abs(signalL)+TMath::Abs(signalR))>sigThresh ) {
3e1a3ad8	294	// Maximum found, mark the position by a negative signal
	295	digits->SetDataUnchecked(row,col-1,time,-signalM);
	296	}
	297	}
	298
	299	}
	300	}
	301	}
	302
	303	// Now check the maxima and calculate the cluster position
	304	for ( row = 0; row < nRowMax ; row++) {
db30bf0f	305	for (time = 0; time < nTimeTotal; time++) {
db30bf0f	306	for ( col = 1; col < nColMax-1; col++) {
3e1a3ad8	307
	308	// Maximum found ?
	309	if (digits->GetDataUnchecked(row,col,time) < 0) {
f7336fa3	310
9d0b222b	311	Int_t iPad;
8230f242	312	for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8	313	Int_t iPadCol = col - 1 + iPad;
	314	clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
	315	,iPadCol
	316	,time));
	317	clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
	318	clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
	319	clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
	320	clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
f7336fa3	321	}
f7336fa3	322
db30bf0f	323	// Count the number of pads in the cluster
	324	Int_t nPadCount = 0;
	325	Int_t ii = 0;
	326	while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
	327	>= sigThresh) {
	328	nPadCount++;
	329	ii++;
	330	if (col-ii < 0) break;
	331	}
	332	ii = 0;
	333	while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
	334	>= sigThresh) {
	335	nPadCount++;
	336	ii++;
	337	if (col+ii+1 >= nColMax) break;
	338	}
	339
	340	nClusters++;
	341	switch (nPadCount) {
	342	case 2:
	343	iType = 0;
	344	nClusters2pad++;
	345	break;
	346	case 3:
	347	iType = 1;
	348	nClusters3pad++;
	349	break;
	350	case 4:
	351	iType = 2;
	352	nClusters4pad++;
	353	break;
	354	case 5:
	355	iType = 3;
	356	nClusters5pad++;
	357	break;
	358	default:
	359	iType = 4;
	360	nClustersLarge++;
	361	break;
	362	};
	363
7ad19338	364	// Look for 5 pad cluster with minimum in the middle
db30bf0f	365	Bool_t fivePadCluster = kFALSE;
3e1a3ad8	366	if (col < nColMax-3) {
3e1a3ad8	367	if (digits->GetDataUnchecked(row,col+2,time) < 0) {
db30bf0f	368	fivePadCluster = kTRUE;
	369	}
	370	if ((fivePadCluster) && (col < nColMax-5)) {
	371	if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
	372	fivePadCluster = kFALSE;
	373	}
	374	}
	375	if ((fivePadCluster) && (col > 1)) {
	376	if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
	377	fivePadCluster = kFALSE;
	378	}
	379	}
	380	}
	381
	382	// 5 pad cluster
	383	// Modify the signal of the overlapping pad for the left part
	384	// of the cluster which remains from a previous unfolding
	385	if (iUnfold) {
	386	clusterSignal[0] *= ratioLeft;
7ad19338	387	iType = 5;
db30bf0f	388	iUnfold = 0;
	389	}
	390
	391	// Unfold the 5 pad cluster
	392	if (fivePadCluster) {
	393	for (iPad = 0; iPad < kNsig; iPad++) {
	394	padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
	395	,col-1+iPad
	396	,time));
f7336fa3	397	}
db30bf0f	398	// Unfold the two maxima and set the signal on
db30bf0f	399	// the overlapping pad to the ratio
17b26de4	400	ratioRight = Unfold(kEpsilon,iplan,padSignal);
db30bf0f	401	ratioLeft = 1.0 - ratioRight;
db30bf0f	402	clusterSignal[2] *= ratioRight;
7ad19338	403	iType = 5;
db30bf0f	404	iUnfold = 1;
f7336fa3	405	}
f7336fa3	406
7ad19338	407
a5cadd36	408	Double_t clusterCharge = clusterSignal[0]
	409	+ clusterSignal[1]
	410	+ clusterSignal[2];
3e1a3ad8	411
db30bf0f	412	// The position of the cluster
3e1a3ad8	413	clusterPads[0] = row + 0.5;
3e1a3ad8	414	// Take the shift of the additional time bins into account
	415	clusterPads[2] = time - nTimeBefore + 0.5;
	416
17b26de4	417	if (fPar->LUTOn()) {
db30bf0f	418	// Calculate the position of the cluster by using the
db30bf0f	419	// lookup table method
7ad19338	420	clusterPads[1] =
7ad19338	421	fPar->LUTposition(iplan,clusterSignal[0]
17b26de4	422	,clusterSignal[1]
17b26de4	423	,clusterSignal[2]);
db30bf0f	424	}
db30bf0f	425	else {
db30bf0f	426	// Calculate the position of the cluster by using the
db30bf0f	427	// center of gravity method
7ad19338	428	for (Int_t i=0;i<5;i++) padSignal[i]=0;
	429	padSignal[2] = TMath::Abs(digits->GetDataUnchecked(row,col,time)); // central pad
	430	padSignal[1] = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); // left pad
	431	padSignal[3] = TMath::Abs(digits->GetDataUnchecked(row,col+1,time)); // right pad
	432	if (col>2 &&TMath::Abs(digits->GetDataUnchecked(row,col-2,time)<padSignal[1])){
	433	padSignal[0] = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
	434	}
	435	if (col<nColMax-3 &&TMath::Abs(digits->GetDataUnchecked(row,col+2,time)<padSignal[3])){
	436	padSignal[4] = TMath::Abs(digits->GetDataUnchecked(row,col+2,time));
	437	}
	438	clusterPads[1] = GetCOG(padSignal);
	439	Double_t check = fPar->LUTposition(iplan,clusterSignal[0]
	440	,clusterSignal[1]
	441	,clusterSignal[2]);
de4b10e5	442	// Float_t diff = clusterPads[1] - check;
db30bf0f	443
	444	}
	445
a5cadd36	446	Double_t q0 = clusterSignal[0];
	447	Double_t q1 = clusterSignal[1];
	448	Double_t q2 = clusterSignal[2];
	449	Double_t clusterSigmaY2 = (q1(q0+q2)+4q0*q2) /
	450	(clusterCharge*clusterCharge);
a819a5f7	451
5443e65e	452	// Calculate the position and the error
7ad19338	453	Double_t colSize = padPlane->GetColSize(col);
7ad19338	454	Double_t rowSize = padPlane->GetRowSize(row);
a5cadd36	455	Double_t clusterPos[3];
7ad19338	456	clusterPos[0] = padPlane->GetColPos(col) + (clusterPads[1]-0.5)*colSize; // MI change
7ad19338	457	clusterPos[1] = padPlane->GetRowPos(row) -0.5*rowSize; //MI change
5443e65e	458	clusterPos[2] = clusterPads[2];
a5cadd36	459	Double_t clusterSig[2];
5443e65e	460	clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
7ad19338	461	clusterSig[1] = rowSize * rowSize / 12.;
f7336fa3	462	// Add the cluster to the output array
bdbb05bb	463	AddCluster(clusterPos
a5cadd36	464	,idet
	465	,clusterCharge
	466	,clusterTracks
	467	,clusterSig
de4b10e5	468	,iType,clusterPads[1]);
f7336fa3	469
f7336fa3	470	}
3e1a3ad8	471	}
	472	}
	473	}
f7336fa3	474
3e1a3ad8	475	// Compress the arrays
	476	digits->Compress(1,0);
	477	track0->Compress(1,0);
7ad19338	478	track1->Compress(1,0);
3e1a3ad8	479	track2->Compress(1,0);
f7336fa3	480
3e1a3ad8	481	// Write the cluster and reset the array
793ff80c	482	WriteClusters(idet);
bdbb05bb	483	ResetRecPoints();
3e1a3ad8	484	}
	485	}
	486	}
f7336fa3	487
47517f42	488	if (fVerbose > 0) {
17b26de4	489	printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42	490	printf("Done.\n");
47517f42	491	}
f7336fa3	492
	493	return kTRUE;
	494
	495	}
	496
7ad19338	497	Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
	498	{
	499	//
	500	// get COG position
	501	// used for clusters with more than 3 pads - where LUT not applicable
	502	Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4];
	503	Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum;
	504	return res;
	505	}
	506
	507
	508
f7336fa3	509	//_____________________________________________________________________________
a5cadd36	510	Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
f7336fa3	511	{
	512	//
	513	// Method to unfold neighbouring maxima.
	514	// The charge ratio on the overlapping pad is calculated
	515	// until there is no more change within the range given by eps.
	516	// The resulting ratio is then returned to the calling method.
	517	//
	518
a5cadd36	519	Int_t irc = 0;
a5cadd36	520	Int_t itStep = 0; // Count iteration steps
f7336fa3	521
a5cadd36	522	Double_t ratio = 0.5; // Start value for ratio
a5cadd36	523	Double_t prevRatio = 0; // Store previous ratio
f7336fa3	524
a5cadd36	525	Double_t newLeftSignal[3] = {0}; // Array to store left cluster signal
	526	Double_t newRightSignal[3] = {0}; // Array to store right cluster signal
	527	Double_t newSignal[3] = {0};
f7336fa3	528
3e1a3ad8	529	// Start the iteration
f7336fa3	530	while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
	531
	532	itStep++;
	533	prevRatio = ratio;
	534
3e1a3ad8	535	// Cluster position according to charge ratio
a5cadd36	536	Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
	537	/ (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
	538	Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
	539	/ ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3	540
3e1a3ad8	541	// Set cluster charge ratio
17b26de4	542	irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
a5cadd36	543	Double_t ampLeft = padSignal[1] / newSignal[1];
17b26de4	544	irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
a5cadd36	545	Double_t ampRight = padSignal[3] / newSignal[1];
f7336fa3	546
3e1a3ad8	547	// Apply pad response to parameters
17b26de4	548	irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
17b26de4	549	irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
f7336fa3	550
3e1a3ad8	551	// Calculate new overlapping ratio
a5cadd36	552	ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] /
db30bf0f	553	(newLeftSignal[2] + newRightSignal[0]));
f7336fa3	554
	555	}
	556
	557	return ratio;
	558
	559	}
	560