#include "AliCaloCalibPedestal.h"
#include "AliEMCALCalibData.h"
#include "AliESDCaloCluster.h"
+#include "AliEMCALUnfolding.h"
ClassImp(AliEMCALClusterizerNxN)
-Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;
-
//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
- : AliEMCALClusterizer()
+ : AliEMCALClusterizer(), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
{
// ctor with the indication of the file where header Tree and digits Tree are stored
}
//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
- : AliEMCALClusterizer(geometry)
+ : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
{
// 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)
+: AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
{
- // ctor, geometry and calibration are initialized elsewhere.
-
+ // ctor, geometry and calibration are initialized elsewhere.
}
-
//____________________________________________________________________________
- AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
+AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
{
// dtor
}
-
//____________________________________________________________________________
void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
{
// Steering method to perform clusterization for the current event
// in AliEMCALLoader
- if(strstr(option,"tim"))
+ if (strstr(option,"tim"))
gBenchmark->Start("EMCALClusterizer");
- if(strstr(option,"print"))
- Print("") ;
+ if (strstr(option,"print"))
+ Print("");
//Get calibration parameters from file or digitizer default values.
- GetCalibrationParameters() ;
+ GetCalibrationParameters();
//Get dead channel map from file or digitizer default values.
- GetCaloCalibPedestal() ;
+ GetCaloCalibPedestal();
- fNumberOfECAClusters = 0;
+ MakeClusters(); //only the real clusters
- MakeClusters() ; //only the real clusters
-
- if(fToUnfold)
- MakeUnfolding() ;
-
- Int_t index ;
+ if (fToUnfold) {
+ fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
+ fClusterUnfolding->MakeUnfolding();
+ }
- //Evaluate position, dispersion and other RecPoint properties for EC section
- for(index = 0; index < fRecPoints->GetEntries(); index++)
- {
+ //Evaluate position, dispersion and other RecPoint properties for EC section
+ for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
- if(rp){
- rp->EvalAll(fECAW0,fDigitsArr) ;
+ if (rp) {
+ rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
//For each rec.point set the distance to the nearest bad crystal
- rp->EvalDistanceToBadChannels(fCaloPed);
+ if (fCaloPed)
+ rp->EvalDistanceToBadChannels(fCaloPed);
}
}
- fRecPoints->Sort() ;
+ fRecPoints->Sort();
- for(index = 0; index < fRecPoints->GetEntries(); index++)
- {
- (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
- (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
+ for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
+ AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
+ if (rp) {
+ rp->SetIndexInList(index);
+ }
+ else AliFatal("RecPoint NULL!!");
}
- fTreeR->Fill();
+ if (fTreeR)
+ fTreeR->Fill();
- if(strstr(option,"deb") || strstr(option,"all"))
- PrintRecPoints(option) ;
+ 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")){
+ 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) ;
+ // 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
+
+ if (fEnergyGrad) { //false by default
+ if (d2->GetCalibAmp()>d1->GetCalibAmp())
+ return 3; // energy of neighboring cell should be smaller in order to become a neighbor
+ }
- // 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
+ Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
+ Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
+ 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) <= fNColDiff && TMath::Abs(rowdiff) <= fNRowDiff)
+ {
+
+ 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
+ // Make clusters
- if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
+ if (fGeom==0)
+ AliFatal("Did not get geometry from EMCALLoader");
- fRecPoints->Clear();
+ fNumberOfECAClusters = 0;
+ fRecPoints->Delete();
- // Set up TObjArray with pointers to digits to work on
- //TObjArray *digitsC = new TObjArray();
+ // Set up TObjArray with pointers to digits to work on, calibrate digits
TObjArray digitsC;
TIter nextdigit(fDigitsArr);
AliEMCALDigit *digit = 0;
- while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
+ while ( (digit = static_cast<AliEMCALDigit*>(nextdigit())) ) {
+ Float_t dEnergyCalibrated = digit->GetAmplitude();
+ Float_t time = digit->GetTime();
+ Calibrate(dEnergyCalibrated,time ,digit->GetId());
+ digit->SetCalibAmp(dEnergyCalibrated);
+ digit->SetTime(time);
digitsC.AddLast(digit);
}
Float_t dMaxEnergyDigit = -1;
AliEMCALDigit *pMaxEnergyDigit = 0;
nextdigitC.Reset();
- while ( (digit =
dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
+ while ( (digit =
static_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)
+ Float_t dEnergyCalibrated = digit->GetCalibAmp();
+
+ if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fMinECut) // no threshold by default!
+ { // needs to be set in OCDB!
+ 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));
+ AliDebug (2, Form("Max digit found: %1.5f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
// keep the candidate digits in a list
TList clusterDigitList;
Double_t clusterCandidateEnergy = dMaxEnergyDigit;
- // now loop over the resto of the digits and cluster into NxN cluster
+ // now loop over the rest 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())) )
+ while ( (digit =
static_cast<AliEMCALDigit *>(nextdigitC())) )
{ // scan over the list of digitsC
if (digit == pMaxEnergyDigit) continue;
- Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
+ Float_t dEnergyCalibrated = digit->GetCalibAmp();
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 !!!!!
+ 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) ;
+ 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) ;
+ if (fNumberOfECAClusters >= fRecPoints->GetSize())
+ fRecPoints->Expand(2*fNumberOfECAClusters+1);
- AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
- fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
- recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
- if(recPoint){
- fNumberOfECAClusters++ ;
+ (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
+ AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
+ // AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
+ // fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
+ // recPoint = static_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters));
+ if (recPoint) {
+ 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?
+ recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
digitsC.Remove(digit);
}
}// recpoint
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;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff