[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerNxN.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.                  *
 **************************************************************************/

// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
// Algorithm:
// 1. peek the most energetic cell
// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
// 3. remove the cells contributing to the cluster
// 4. start from 1 for the remaining clusters
// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy 
// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
// Use Case:
//  root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")  
//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
//               //reads gAlice from header file "..."                      
//  root [1] cl->ExecuteTask()  
//               //finds RecPoints in all events stored in galice.root
//  root [2] cl->SetDigitsBranch("digits2") 
//               //sets another title for Digitis (input) branch
//  root [3] cl->SetRecPointsBranch("recp2")  
//               //sets another title four output branches
//  root [4] cl->SetTowerLocalMaxCut(0.03)  
//               //set clusterization parameters
//  root [5] cl->ExecuteTask("deb all time")  
//               //once more finds RecPoints options are 
//               // deb - print number of found rec points
//               // deb all - print number of found RecPoints and some their characteristics 
//               // time - print benchmarking results

// --- ROOT system ---
#include <TMath.h> 
#include <TMinuit.h>
#include <TTree.h> 
#include <TBenchmark.h>
#include <TBrowser.h>
#include <TROOT.h>
#include <TClonesArray.h>

// --- Standard library ---
#include <cassert>

// --- AliRoot header files ---
#include "AliLog.h"
#include "AliEMCALClusterizerNxN.h"
#include "AliEMCALRecPoint.h"
#include "AliEMCALDigit.h"
#include "AliEMCALGeometry.h"
#include "AliCaloCalibPedestal.h"
#include "AliEMCALCalibData.h"
#include "AliESDCaloCluster.h"

ClassImp(AliEMCALClusterizerNxN)

Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;

//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
  : AliEMCALClusterizer()
{
  // ctor with the indication of the file where header Tree and digits Tree are stored
}

//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
  : AliEMCALClusterizer(geometry)
{
  // ctor with the indication of the file where header Tree and digits Tree are stored
  // use this contructor to avoid usage of Init() which uses runloader
  // change needed by HLT - MP

}

//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
: AliEMCALClusterizer(geometry, calib, caloped)
{
	// ctor, geometry and calibration are initialized elsewhere.
				
}


//____________________________________________________________________________
  AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
{
  // dtor
}


//____________________________________________________________________________
void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
{
  // Steering method to perform clusterization for the current event 
  // in AliEMCALLoader

  if(strstr(option,"tim"))
    gBenchmark->Start("EMCALClusterizer"); 
  
  if(strstr(option,"print"))
    Print("") ; 
 
  //Get calibration parameters from file or digitizer default values.
  GetCalibrationParameters() ;

  //Get dead channel map from file or digitizer default values.
  GetCaloCalibPedestal() ;
	
  fNumberOfECAClusters = 0;

  MakeClusters() ;  //only the real clusters

  if(fToUnfold)
    MakeUnfolding() ;

  Int_t index ;

  //Evaluate position, dispersion and other RecPoint properties for EC section                      
  for(index = 0; index < fRecPoints->GetEntries(); index++) 
    {
      dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
      AliDebug(5, Form("MAX INDEX %d ", dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->GetMaximalEnergyIndex()));
      //For each rec.point set the distance to the nearest bad crystal
      dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
    }
  
  fRecPoints->Sort() ;

  for(index = 0; index < fRecPoints->GetEntries(); index++) 
    {
      (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
      (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
    }
  
  fTreeR->Fill();
  
  if(strstr(option,"deb") || strstr(option,"all"))  
    PrintRecPoints(option) ;

  AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));

  fRecPoints->Delete();

  if(strstr(option,"tim")){
    gBenchmark->Stop("EMCALClusterizer");
    printf("Exec took %f seconds for Clusterizing", 
	   gBenchmark->GetCpuTime("EMCALClusterizer"));
  }    
}

//____________________________________________________________________________
Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
{
	// Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
	//                                       = 1 are neighbour
	//                                       = 2 is in different SM; continue searching 
	// In case it is in different SM, but same phi rack, check if neigbours at eta=0
	// neighbours are defined as digits having at least a common side 
	// The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
	//                                      which is compared to a digit (d2)  not yet in a cluster  
	
	static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
	static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
	static Int_t rowdiff=0, coldiff=0;

	shared = kFALSE;
	
	fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
	fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
	fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
	fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
	
	//If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
	if(nSupMod1 != nSupMod2 ) 
	  {
	    //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
	    Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
	    Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
	    
	    if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
	    
	    // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
	    // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
	    if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
	    else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
	    
	    shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
	    
	}//Different SM, same phi
	
	rowdiff = TMath::Abs(iphi1 - iphi2);  
	coldiff = TMath::Abs(ieta1 - ieta2) ;  

	// neighbours +-1 in col and row
	if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
	  {
	    
	    AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
			     d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
	    
	    return 1;
	  }//Neighbours
	else 
	  {
	    AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
			     d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
	    shared = kFALSE;
	    return 2 ; 
	  }//Not neighbours
}

//____________________________________________________________________________
void AliEMCALClusterizerNxN::MakeClusters()
{
  // Steering method to construct the clusters stored in a list of Reconstructed Points
  // A cluster is defined as a list of neighbour digits
  // Mar 03, 2007 by PAI

  if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");

  fRecPoints->Clear();

  // Set up TObjArray with pointers to digits to work on 
  //TObjArray *digitsC = new TObjArray();
  TObjArray digitsC;
  TIter nextdigit(fDigitsArr);
  AliEMCALDigit *digit = 0;
  while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
    digitsC.AddLast(digit);
  }

  TIter nextdigitC(&digitsC);

  AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f)\n",
		  fDigitsArr->GetEntries(),fMinECut));

  Bool_t bDone = kFALSE;
  while ( bDone != kTRUE )
    {
      //first sort the digits:
      Int_t iMaxEnergyDigit = -1;
      Float_t dMaxEnergyDigit = -1;
      AliEMCALDigit *pMaxEnergyDigit = 0;
      nextdigitC.Reset();
      while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) 
	{ // scan over the list of digitsC
	  Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	  //AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
	  
	  //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold  )
	  if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
	    {
	      if (dEnergyCalibrated > dMaxEnergyDigit)
		{
		  dMaxEnergyDigit = dEnergyCalibrated;
		  iMaxEnergyDigit = digit->GetId();
		  pMaxEnergyDigit = digit;
		}
	    }
	}

      if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0) 
	{
	  bDone = kTRUE;
	  continue;
	}

      AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
      AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));

      // keep the candidate digits in a list
      TList clusterDigitList;
      clusterDigitList.SetOwner(kFALSE);
      clusterDigitList.AddLast(pMaxEnergyDigit);	 

      Double_t clusterCandidateEnergy = dMaxEnergyDigit;

      // now loop over the resto of the digits and cluster into NxN cluster 
      // we do not actually cluster yet: we keep them in the list clusterDigitList
      nextdigitC.Reset();
      while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) 
	{ // scan over the list of digitsC
	  if (digit == pMaxEnergyDigit) continue;
	  Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	  AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
	  //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold  )
	  if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0  )
	    {
	      Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
	      if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
	      Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
	      if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!! 
		{      
		  clusterDigitList.AddLast(digit) ;
		  clusterCandidateEnergy += dEnergyCalibrated;
		}
	    }
	}// loop over the next digits

      // start a cluster here only if a cluster energy is larger than clustering threshold
      //if (clusterCandidateEnergy > 0.1)
      AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
      if (clusterCandidateEnergy > fECAClusteringThreshold)
	{
	  if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
      
	  AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ; 
	  fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
	  recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ; 
	  fNumberOfECAClusters++ ;       
	  recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);

	  AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
	  for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
	    {

	      digit = (AliEMCALDigit*)clusterDigitList.At(idig);
	      Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	      AliDebug(5, Form(" Adding digit %d", digit->GetId()));
	      // note: this way the sharing info is lost!
	      recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
	      digitsC.Remove(digit); 		  
	    }
	}
      else
	{
	  // we do not want to start clustering in the same spot!
	  // but in this case we may NOT reuse this seed for another cluster!
	  // need a better bookeeping?
	  digitsC.Remove(pMaxEnergyDigit);
	}

      AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));      
    } // while ! done 

  //delete digitsC ; //nope we use an object
  
  AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
}

//____________________________________________________________________________
void AliEMCALClusterizerNxN::MakeUnfolding()
{
  // Unfolds clusters using the shape of an ElectroMagnetic shower
  // Performs unfolding of all clusters
		
  if(fNumberOfECAClusters > 0){
    if (fGeom==0)
      AliFatal("Did not get geometry from EMCALLoader") ;
    Int_t nModulesToUnfold = fGeom->GetNCells();

    Int_t numberofNotUnfolded = fNumberOfECAClusters ;
    Int_t index ;
    for(index = 0 ; index < numberofNotUnfolded ; index++){

      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;

      TVector3 gpos;
      Int_t absId;
      recPoint->GetGlobalPosition(gpos);
      fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
      if(absId > nModulesToUnfold)
        break ;

      Int_t nMultipl = recPoint->GetMultiplicity() ;
      AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
      Float_t * maxAtEnergy = new Float_t[nMultipl] ;
      Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;

      if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
        //UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
        fRecPoints->Remove(recPoint);
        fRecPoints->Compress() ;
        index-- ;
        fNumberOfECAClusters-- ;
        numberofNotUnfolded-- ;
      }
      else{
        recPoint->SetNExMax(1) ; //Only one local maximum
      }

      delete[] maxAt ;
      delete[] maxAtEnergy ;
    }
  }
  // End of Unfolding of clusters
}

//____________________________________________________________________________
Double_t  AliEMCALClusterizerNxN::ShowerShape(Double_t x, Double_t y)
{ 
  // Shape of the shower
  // If you change this function, change also the gradient evaluation in ChiSquare()
  //
  Double_t r = sqrt(x*x+y*y);
  Double_t r133  = TMath::Power(r, 1.33) ;
  Double_t r669  = TMath::Power(r, 6.69) ;
  Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
  return shape ;
}

//____________________________________________________________________________
void  AliEMCALClusterizerNxN::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, 
					    Int_t /*nMax*/, 
					    AliEMCALDigit ** /*maxAt*/, 
					    Float_t * /*maxAtEnergy*/)
{
  //
  // Performs the unfolding of a cluster with nMax overlapping showers 
  //
  AliWarning("Not implemented. To be.");
}

//___________________________________________________________________
void   AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
{
  //
  // input is calibrated - the case when we run already on ESD
  //
  AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;
}
Commit	Line	Data
ee08edde	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	// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
	17	// Algorithm:
	18	// 1. peek the most energetic cell
	19	// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
	20	// 3. remove the cells contributing to the cluster
	21	// 4. start from 1 for the remaining clusters
	22	// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
	23	// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy
	24	// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
	25	// Use Case:
	26	// root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")
	27	// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
	28	// //reads gAlice from header file "..."
	29	// root [1] cl->ExecuteTask()
	30	// //finds RecPoints in all events stored in galice.root
	31	// root [2] cl->SetDigitsBranch("digits2")
	32	// //sets another title for Digitis (input) branch
	33	// root [3] cl->SetRecPointsBranch("recp2")
	34	// //sets another title four output branches
	35	// root [4] cl->SetTowerLocalMaxCut(0.03)
	36	// //set clusterization parameters
	37	// root [5] cl->ExecuteTask("deb all time")
	38	// //once more finds RecPoints options are
	39	// // deb - print number of found rec points
	40	// // deb all - print number of found RecPoints and some their characteristics
	41	// // time - print benchmarking results
	42
	43	// --- ROOT system ---
	44	#include <TMath.h>
	45	#include <TMinuit.h>
	46	#include <TTree.h>
	47	#include <TBenchmark.h>
	48	#include <TBrowser.h>
	49	#include <TROOT.h>
	50	#include <TClonesArray.h>
	51
	52	// --- Standard library ---
	53	#include <cassert>
	54
	55	// --- AliRoot header files ---
	56	#include "AliLog.h"
	57	#include "AliEMCALClusterizerNxN.h"
	58	#include "AliEMCALRecPoint.h"
	59	#include "AliEMCALDigit.h"
	60	#include "AliEMCALGeometry.h"
	61	#include "AliCaloCalibPedestal.h"
	62	#include "AliEMCALCalibData.h"
	63	#include "AliESDCaloCluster.h"
	64
65	ClassImp(AliEMCALClusterizerNxN)
66
67	Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;
68
69	//____________________________________________________________________________
70	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
71	: AliEMCALClusterizer()
72	{
73	// ctor with the indication of the file where header Tree and digits Tree are stored
74	}
75
76	//____________________________________________________________________________
77	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
78	: AliEMCALClusterizer(geometry)
79	{
80	// ctor with the indication of the file where header Tree and digits Tree are stored
81	// use this contructor to avoid usage of Init() which uses runloader
82	// change needed by HLT - MP
83
84	}
85
86	//____________________________________________________________________________
87	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
88	: AliEMCALClusterizer(geometry, calib, caloped)
89	{
90	// ctor, geometry and calibration are initialized elsewhere.
91
92	}
93
94
95	//____________________________________________________________________________
96	AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
97	{
98	// dtor
99	}
100
101
102	//____________________________________________________________________________
103	void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
104	{
105	// Steering method to perform clusterization for the current event
106	// in AliEMCALLoader
107
108	if(strstr(option,"tim"))
109	gBenchmark->Start("EMCALClusterizer");
110
111	if(strstr(option,"print"))
112	Print("") ;
113
114	//Get calibration parameters from file or digitizer default values.
115	GetCalibrationParameters() ;
116
117	//Get dead channel map from file or digitizer default values.
118	GetCaloCalibPedestal() ;
119
120	fNumberOfECAClusters = 0;
121
122	MakeClusters() ; //only the real clusters
123
124	if(fToUnfold)
125	MakeUnfolding() ;
126
127	Int_t index ;
128
129	//Evaluate position, dispersion and other RecPoint properties for EC section
130	for(index = 0; index < fRecPoints->GetEntries(); index++)
131	{
132	dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
133	AliDebug(5, Form("MAX INDEX %d ", dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->GetMaximalEnergyIndex()));
134	//For each rec.point set the distance to the nearest bad crystal
135	dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
136	}
137
138	fRecPoints->Sort() ;
139
140	for(index = 0; index < fRecPoints->GetEntries(); index++)
141	{
142	(dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
143	(dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
144	}
145
146	fTreeR->Fill();
147
148	if(strstr(option,"deb") \|\| strstr(option,"all"))
149	PrintRecPoints(option) ;
150
151	AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
152
153	fRecPoints->Delete();
154
155	if(strstr(option,"tim")){
156	gBenchmark->Stop("EMCALClusterizer");
157	printf("Exec took %f seconds for Clusterizing",
158	gBenchmark->GetCpuTime("EMCALClusterizer"));
159	}
160	}
161
162	//____________________________________________________________________________
163	Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
164	{
165	// Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
166	// = 1 are neighbour
167	// = 2 is in different SM; continue searching
168	// In case it is in different SM, but same phi rack, check if neigbours at eta=0
169	// neighbours are defined as digits having at least a common side
170	// The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
171	// which is compared to a digit (d2) not yet in a cluster
172
173	static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
174	static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
175	static Int_t rowdiff=0, coldiff=0;
176
177	shared = kFALSE;
178
179	fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
180	fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
181	fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
182	fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
183
184	//If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
185	if(nSupMod1 != nSupMod2 )
186	{
187	//Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
188	Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
189	Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
190
191	if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
192
193	// In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
194	// C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
195	if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
196	else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
197
198	shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
199
200	}//Different SM, same phi
201
202	rowdiff = TMath::Abs(iphi1 - iphi2);
203	coldiff = TMath::Abs(ieta1 - ieta2) ;
204
205	// neighbours +-1 in col and row
206	if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
207	{
208
209	AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
210	d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
211
212	return 1;
213	}//Neighbours
214	else
215	{
216	AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
217	d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
218	shared = kFALSE;
219	return 2 ;
220	}//Not neighbours
221	}
222
223	//____________________________________________________________________________
224	void AliEMCALClusterizerNxN::MakeClusters()
225	{
226	// Steering method to construct the clusters stored in a list of Reconstructed Points
227	// A cluster is defined as a list of neighbour digits
228	// Mar 03, 2007 by PAI
229
230	if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
231
232	fRecPoints->Clear();
233
234	// Set up TObjArray with pointers to digits to work on
235	//TObjArray *digitsC = new TObjArray();
236	TObjArray digitsC;
237	TIter nextdigit(fDigitsArr);
238	AliEMCALDigit *digit = 0;
239	while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
240	digitsC.AddLast(digit);
241	}
242
243	TIter nextdigitC(&digitsC);
244
245	AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f)\n",
246	fDigitsArr->GetEntries(),fMinECut));
247
248	Bool_t bDone = kFALSE;
249	while ( bDone != kTRUE )
250	{
251	//first sort the digits:
252	Int_t iMaxEnergyDigit = -1;
253	Float_t dMaxEnergyDigit = -1;
254	AliEMCALDigit *pMaxEnergyDigit = 0;
255	nextdigitC.Reset();
256	while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
257	{ // scan over the list of digitsC
258	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
259	//AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
260
261	//if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
262	if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
263	{
264	if (dEnergyCalibrated > dMaxEnergyDigit)
265	{
266	dMaxEnergyDigit = dEnergyCalibrated;
267	iMaxEnergyDigit = digit->GetId();
268	pMaxEnergyDigit = digit;
269	}
270	}
271	}
272
273	if (iMaxEnergyDigit < 0 \|\| digitsC.GetEntries() <= 0)
274	{
275	bDone = kTRUE;
276	continue;
277	}
278
279	AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
280	AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));
281
282	// keep the candidate digits in a list
283	TList clusterDigitList;
284	clusterDigitList.SetOwner(kFALSE);
285	clusterDigitList.AddLast(pMaxEnergyDigit);
286
287	Double_t clusterCandidateEnergy = dMaxEnergyDigit;
288
289	// now loop over the resto of the digits and cluster into NxN cluster
290	// we do not actually cluster yet: we keep them in the list clusterDigitList
291	nextdigitC.Reset();
292	while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
293	{ // scan over the list of digitsC
294	if (digit == pMaxEnergyDigit) continue;
295	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
296	AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
297	//if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
298	if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 )
299	{
300	Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
301	if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
302	Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
303	if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
304	{
305	clusterDigitList.AddLast(digit) ;
306	clusterCandidateEnergy += dEnergyCalibrated;
307	}
308	}
309	}// loop over the next digits
310
311	// start a cluster here only if a cluster energy is larger than clustering threshold
312	//if (clusterCandidateEnergy > 0.1)
313	AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
314	if (clusterCandidateEnergy > fECAClusteringThreshold)
315	{
316	if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
317
318	AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
319	fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
320	recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
321	fNumberOfECAClusters++ ;
c8fe2783	322	recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
ee08edde	323
	324	AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
	325	for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
	326	{
	327
	328	digit = (AliEMCALDigit*)clusterDigitList.At(idig);
	329	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	330	AliDebug(5, Form(" Adding digit %d", digit->GetId()));
	331	// note: this way the sharing info is lost!
	332	recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
	333	digitsC.Remove(digit);
	334	}
	335	}
	336	else
	337	{
	338	// we do not want to start clustering in the same spot!
	339	// but in this case we may NOT reuse this seed for another cluster!
	340	// need a better bookeeping?
	341	digitsC.Remove(pMaxEnergyDigit);
	342	}
	343
	344	AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
	345	} // while ! done
	346
	347	//delete digitsC ; //nope we use an object
	348
	349	AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
	350	}
	351
	352	//____________________________________________________________________________
	353	void AliEMCALClusterizerNxN::MakeUnfolding()
	354	{
	355	// Unfolds clusters using the shape of an ElectroMagnetic shower
	356	// Performs unfolding of all clusters
	357
	358	if(fNumberOfECAClusters > 0){
	359	if (fGeom==0)
	360	AliFatal("Did not get geometry from EMCALLoader") ;
	361	Int_t nModulesToUnfold = fGeom->GetNCells();
	362
	363	Int_t numberofNotUnfolded = fNumberOfECAClusters ;
	364	Int_t index ;
	365	for(index = 0 ; index < numberofNotUnfolded ; index++){
	366
	367	AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
	368
	369	TVector3 gpos;
	370	Int_t absId;
	371	recPoint->GetGlobalPosition(gpos);
	372	fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
	373	if(absId > nModulesToUnfold)
	374	break ;
	375
	376	Int_t nMultipl = recPoint->GetMultiplicity() ;
	377	AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
	378	Float_t * maxAtEnergy = new Float_t[nMultipl] ;
	379	Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
	380
	381	if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
	382	//UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
	383	fRecPoints->Remove(recPoint);
	384	fRecPoints->Compress() ;
	385	index-- ;
	386	fNumberOfECAClusters-- ;
387	numberofNotUnfolded-- ;
388	}
389	else{
390	recPoint->SetNExMax(1) ; //Only one local maximum
391	}
392
393	delete[] maxAt ;
394	delete[] maxAtEnergy ;
395	}
396	}
397	// End of Unfolding of clusters
398	}
399
400	//____________________________________________________________________________
401	Double_t AliEMCALClusterizerNxN::ShowerShape(Double_t x, Double_t y)
402	{
403	// Shape of the shower
404	// If you change this function, change also the gradient evaluation in ChiSquare()
405	//
406	Double_t r = sqrt(xx+yy);
407	Double_t r133 = TMath::Power(r, 1.33) ;
408	Double_t r669 = TMath::Power(r, 6.69) ;
409	Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
410	return shape ;
411	}
412
413	//____________________________________________________________________________
414	void AliEMCALClusterizerNxN::UnfoldCluster(AliEMCALRecPoint * /iniTower/,
415	Int_t /nMax/,
416	AliEMCALDigit ** /maxAt/,
417	Float_t * /maxAtEnergy/)
418	{
419	//
420	// Performs the unfolding of a cluster with nMax overlapping showers
421	//
422	AliWarning("Not implemented. To be.");
423	}
424
425	//___________________________________________________________________
426	void AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
427	{
428	//
429	// input is calibrated - the case when we run already on ESD
430	//
431	AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;
432	}