From: gconesab Date: Thu, 12 Aug 2010 07:31:42 +0000 (+0000) Subject: New Clusterizer that makes cluster grouping neighbour cells around highest energy... X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=ee08eddeb6496760fd9a0a31366df16ba829bb21;p=u%2Fmrichter%2FAliRoot.git New Clusterizer that makes cluster grouping neighbour cells around highest energy tower, 1 to N towers away from this tower (Mateusz) --- diff --git a/EMCAL/AliEMCALClusterizer.cxx b/EMCAL/AliEMCALClusterizer.cxx index bf4bda56483..90b0cb437a2 100644 --- a/EMCAL/AliEMCALClusterizer.cxx +++ b/EMCAL/AliEMCALClusterizer.cxx @@ -19,21 +19,43 @@ // Base class for the clusterization algorithm (pure abstract) //*-- //*-- Author: Yves Schutz SUBATECH -// Modif: -// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction -// of new IO (à la PHOS) +// +// Clusterization mother class. Contains common methods/data members of different +// clusterizers. GCB 2010 ////////////////////////////////////////////////////////////////////////////// // --- ROOT system --- #include "TClonesArray.h" #include "TTree.h" +#include +class TFolder; +#include +#include +#include +class TSystem; +#include +#include +#include // --- Standard library --- - +#include // --- AliRoot header files --- #include "AliEMCALClusterizer.h" +#include "AliEMCALReconstructor.h" +#include "AliRunLoader.h" +#include "AliRun.h" #include "AliLog.h" +#include "AliEMCAL.h" +#include "AliEMCALRecPoint.h" +#include "AliEMCALRecParam.h" +#include "AliEMCALGeometry.h" +#include "AliEMCALRecParam.h" +#include "AliCDBManager.h" +#include "AliCaloCalibPedestal.h" +#include "AliEMCALCalibData.h" +class AliCDBStorage; +#include "AliCDBEntry.h" ClassImp(AliEMCALClusterizer) @@ -41,11 +63,74 @@ ClassImp(AliEMCALClusterizer) AliEMCALClusterizer::AliEMCALClusterizer(): fDigitsArr(NULL), fTreeR(NULL), - fRecPoints(NULL) + fRecPoints(NULL), + fGeom(NULL), + fCalibData(NULL), + fCaloPed(NULL), + fADCchannelECA(0.),fADCpedestalECA(0.), + fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.), + fDefaultInit(kFALSE),fToUnfold(kFALSE), + fNumberOfECAClusters(0), fECAClusteringThreshold(0.), + fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.) { // ctor + + Init(); + +} + +//____________________________________________________________________________ +AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry): + fDigitsArr(NULL), + fTreeR(NULL), + fRecPoints(NULL), + fGeom(geometry), + fCalibData(NULL), + fCaloPed(NULL), + fADCchannelECA(0.),fADCpedestalECA(0.), + fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.), + fDefaultInit(kFALSE),fToUnfold(kFALSE), + fNumberOfECAClusters(0), fECAClusteringThreshold(0.), + fECALocMaxCut(0.),fECAW0(0.),fMinECut(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 + + // Note for the future: the use on runloader should be avoided or optional at least + // another way is to make Init virtual and protected at least such that the deriving classes can overload + // Init() ; + // + + if (!fGeom) + { + AliFatal("Geometry not initialized."); + } + } +//____________________________________________________________________________ +AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped): + fDigitsArr(NULL), + fTreeR(NULL), + fRecPoints(NULL), + fGeom(geometry), + fCalibData(calib), + fCaloPed(caloped), + fADCchannelECA(0.),fADCpedestalECA(0.), + fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.), + fDefaultInit(kFALSE),fToUnfold(kFALSE), + fNumberOfECAClusters(0), fECAClusteringThreshold(0.), + fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.) +{ + // ctor, geometry and calibration are initialized elsewhere. + + if (!fGeom) + AliFatal("Geometry not initialized."); + +} + + //____________________________________________________________________________ AliEMCALClusterizer::~AliEMCALClusterizer() { @@ -60,6 +145,240 @@ AliEMCALClusterizer::~AliEMCALClusterizer() } } +//____________________________________________________________________________ +Float_t AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId) +{ + + // Convert digitized amplitude into energy. + // Calibration parameters are taken from calibration data base for raw data, + // or from digitizer parameters for simulated data. + if(fCalibData){ + + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader") ; + + Int_t iSupMod = -1; + Int_t nModule = -1; + Int_t nIphi = -1; + Int_t nIeta = -1; + Int_t iphi = -1; + Int_t ieta = -1; + + Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ; + if(!bCell) { + fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", absId); + assert(0); + } + + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); + + // Check if channel is bad (dead or hot), in this case return 0. + // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation. + // for the moment keep it here but remember to do the selection at the sdigitizer level + // and remove it from here + Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi); + if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) { + AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus)); + return 0; + } + //Check if time is too large or too small, indication of a noisy channel, remove in this case + if(time > fTimeMax || time < fTimeMin) return 0; + + fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); + fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); + + return -fADCpedestalECA + amp * fADCchannelECA ; + + } + else //Return energy with default parameters if calibration is not available + return -fADCpedestalECA + amp * fADCchannelECA ; + +} + +//____________________________________________________________________________ +void AliEMCALClusterizer::GetCalibrationParameters() +{ + // Set calibration parameters: + // if calibration database exists, they are read from database, + // otherwise, they are taken from digitizer. + // + // It is a user responsilibity to open CDB before reconstruction, + // for example: + // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + + //Check if calibration is stored in data base + + if(!fCalibData) + { + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); + if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); + } + + if(!fCalibData) + AliFatal("Calibration parameters not found in CDB!"); + +} + +//____________________________________________________________________________ +void AliEMCALClusterizer::GetCaloCalibPedestal() +{ + // Set calibration parameters: + // if calibration database exists, they are read from database, + // otherwise, they are taken from digitizer. + // + // It is a user responsilibity to open CDB before reconstruction, + // for example: + // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + + //Check if calibration is stored in data base + + if(!fCaloPed) + { + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"); + if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject(); + } + + if(!fCaloPed) + AliFatal("Pedestal info not found in CDB!"); + +} + +//____________________________________________________________________________ +void AliEMCALClusterizer::Init() +{ + // Make all memory allocations which can not be done in default constructor. + // Attach the Clusterizer task to the list of EMCAL tasks + + AliRunLoader *rl = AliRunLoader::Instance(); + if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) + fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + else + fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); + + AliDebug(1,Form("geom %p",fGeom)); + + if(!gMinuit) + gMinuit = new TMinuit(100) ; + +} + +//____________________________________________________________________________ +void AliEMCALClusterizer::InitParameters() +{ + // Initializes the parameters for the Clusterizer + fNumberOfECAClusters = 0 ; + fCalibData = 0 ; + fCaloPed = 0 ; + + const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam(); + if(!recParam) { + AliFatal("Reconstruction parameters for EMCAL not set!"); + } else { + fECAClusteringThreshold = recParam->GetClusteringThreshold(); + fECAW0 = recParam->GetW0(); + fMinECut = recParam->GetMinECut(); + fToUnfold = recParam->GetUnfold(); + if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!"); + fECALocMaxCut = recParam->GetLocMaxCut(); + fTimeCut = recParam->GetTimeCut(); + fTimeMin = recParam->GetTimeMin(); + fTimeMax = recParam->GetTimeMax(); + + AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s", + fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax)); + } + +} + + +//____________________________________________________________________________ +void AliEMCALClusterizer::Print(Option_t * /*option*/)const +{ + // Print clusterizer parameters + + TString message("\n") ; + + if( strcmp(GetName(), "") !=0 ){ + + // Print parameters + + TString taskName(Version()) ; + + printf("--------------- "); + printf("%s",taskName.Data()) ; + printf(" "); + printf("Clusterizing digits: "); + printf("\n ECA Local Maximum cut = %f", fECALocMaxCut); + printf("\n ECA Logarithmic weight = %f", fECAW0); + if(fToUnfold) + printf("\nUnfolding on\n"); + else + printf("\nUnfolding off\n"); + + printf("------------------------------------------------------------------"); + } + else + printf("AliEMCALClusterizer not initialized ") ; +} + +//____________________________________________________________________________ +void AliEMCALClusterizer::PrintRecPoints(Option_t * option) +{ + // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer + if(strstr(option,"deb")) { + printf("PrintRecPoints: Clusterization result:") ; + + printf(" Found %d ECA Rec Points\n ", + fRecPoints->GetEntriesFast()) ; + } + + if(strstr(option,"all")) { + if(strstr(option,"deb")) { + printf("\n-----------------------------------------------------------------------\n") ; + printf("Clusters in ECAL section\n") ; + printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + } + Int_t index; + for (index = 0 ; index < fRecPoints->GetEntries() ; index++) { + AliEMCALRecPoint * rp = dynamic_cast(fRecPoints->At(index)) ; + TVector3 globalpos; + //rp->GetGlobalPosition(globalpos); + TVector3 localpos; + rp->GetLocalPosition(localpos); + Float_t lambda[2]; + rp->GetElipsAxis(lambda); + + Int_t nprimaries=0; + Int_t * primaries = rp->GetPrimaries(nprimaries); + + if(strstr(option,"deb")) + printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", + rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), + globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), + rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; + if(strstr(option,"deb")){ + for (Int_t iprimary=0; iprimaryBranch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split); } + + + diff --git a/EMCAL/AliEMCALClusterizer.h b/EMCAL/AliEMCALClusterizer.h index 388213f078a..33584e6a4b5 100644 --- a/EMCAL/AliEMCALClusterizer.h +++ b/EMCAL/AliEMCALClusterizer.h @@ -19,6 +19,9 @@ class TTree; // --- Standard library --- // --- AliRoot header files --- +class AliEMCALGeometry ; +class AliEMCALCalibData ; +class AliCaloCalibPedestal ; class AliEMCALClusterizer : public TObject { @@ -26,22 +29,47 @@ public: AliEMCALClusterizer() ; // default ctor virtual ~AliEMCALClusterizer() ; // dtorEM + AliEMCALClusterizer(AliEMCALGeometry* geometry); + AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * pedestal); virtual void Digits2Clusters(Option_t *option) = 0; - virtual Float_t GetTimeCut() const = 0; - - virtual void SetECAClusteringThreshold(Float_t) = 0; - virtual void SetECALocalMaxCut(Float_t) = 0; - virtual void SetECALogWeight(Float_t) = 0; - virtual void SetTimeCut(Float_t) = 0; - virtual void SetUnfolding(Bool_t) = 0; - - virtual const char * Version() const {Warning("Version", "Not Defined") ; return 0 ; } + virtual Float_t Calibrate(const Float_t amp, const Float_t time, const Int_t cellId) ; // Tranforms Amp to energy + virtual void Init() ; + virtual void InitParameters() ; //{ AliInfo("Overload this method."); } + + //Get/Set reconstruction parameters + virtual void GetCalibrationParameters(void) ; + virtual void GetCaloCalibPedestal(void) ; + virtual void SetCalibrationParameters(AliEMCALCalibData * calib) { fCalibData = calib ; } + virtual void SetCaloCalibPedestal(AliCaloCalibPedestal * caloped) { fCaloPed = caloped ; } + + virtual Float_t GetTimeMin() const { return fTimeMin ; } + virtual Float_t GetTimeMax() const { return fTimeMax ; } + virtual Float_t GetTimeCut() const { return fTimeCut ; } + virtual void GetNumberOfClustersFound(int numb )const { numb = fNumberOfECAClusters ;} + virtual Float_t GetECAClusteringThreshold() const { return fECAClusteringThreshold;} + virtual Float_t GetECALocalMaxCut() const { return fECALocMaxCut;} + virtual Float_t GetECALogWeight() const { return fECAW0;} + virtual Float_t GetMinECut() const { return fMinECut;} + + virtual void SetTimeMin(Float_t t) { fTimeMin = t ;} + virtual void SetTimeMax(Float_t t) { fTimeMax = t ;} + virtual void SetTimeCut(Float_t t) { fTimeCut = t ;} + virtual void SetECAClusteringThreshold(Float_t th) { fECAClusteringThreshold = th ; } + virtual void SetMinECut(Float_t mine) { fMinECut = mine; } + virtual void SetECALocalMaxCut(Float_t cut) { fECALocMaxCut = cut ; } + virtual void SetECALogWeight(Float_t w) { fECAW0 = w ; } + virtual void SetUnfolding(Bool_t toUnfold = kTRUE ) {fToUnfold = toUnfold ;} virtual void SetInput(TTree *digitsTree); virtual void SetOutput(TTree *clustersTree); - virtual void InitParameters() { AliInfo("Overload this method."); } + + virtual void Print(Option_t * option)const ; + virtual void PrintRecPoints(Option_t * option); + virtual void PrintRecoInfo(); //*MENU* + + virtual const char * Version() const {Warning("Version", "Not Defined") ; return 0 ; } protected: @@ -50,12 +78,34 @@ protected: TClonesArray *fDigitsArr; // Array with EMCAL digits TTree *fTreeR; // Tree with output clusters TObjArray *fRecPoints; // Array with EMCAL clusters - + + AliEMCALGeometry * fGeom; //! pointer to geometry for utilities + AliEMCALCalibData * fCalibData ; //! Calibration database if aval + AliCaloCalibPedestal * fCaloPed ; //! Tower status map if aval + + Float_t fADCchannelECA ; // width of one ADC channel for EC section (GeV) + Float_t fADCpedestalECA ; // pedestal of ADC for EC section (GeV) + + Float_t fTimeMin ; // Minimum time of physical signal in a cell/digit + Float_t fTimeMax ; // Maximum time of physical signal in a cell/digit + Float_t fTimeCut ; // Maximum time difference between the digits inside EMC cluster + + Bool_t fDefaultInit; //! Says if the task was created by defaut ctor (only parameters are initialized) + Bool_t fToUnfold ; // To perform unfolding + Int_t fNumberOfECAClusters ; // number of clusters found in EC section + + Float_t fECAClusteringThreshold ; // minimum energy to seed a EC digit in a cluster + Float_t fECALocMaxCut ; // minimum energy difference to distinguish local maxima in a cluster + Float_t fECAW0 ; // logarithmic weight for the cluster center of gravity calculation + Float_t fMinECut; // Minimum energy for a digit to be a member of a cluster + private: AliEMCALClusterizer(const AliEMCALClusterizer &); //copy ctor AliEMCALClusterizer & operator = (const AliEMCALClusterizer &); - - ClassDef(AliEMCALClusterizer,1) // Clusterization algorithm class + + + + ClassDef(AliEMCALClusterizer,2) // Clusterization algorithm class } ; #endif // AliEMCALCLUSTERIZER_H diff --git a/EMCAL/AliEMCALClusterizerNxN.cxx b/EMCAL/AliEMCALClusterizerNxN.cxx new file mode 100644 index 00000000000..98eebec2428 --- /dev/null +++ b/EMCAL/AliEMCALClusterizerNxN.cxx @@ -0,0 +1,432 @@ +/************************************************************************** + * 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 : 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 +#include +#include +#include +#include +#include +#include + +// --- Standard library --- +#include + +// --- 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(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ; + AliDebug(5, Form("MAX INDEX %d ", dynamic_cast(fRecPoints->At(index))->GetMaximalEnergyIndex())); + //For each rec.point set the distance to the nearest bad crystal + dynamic_cast(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed); + } + + fRecPoints->Sort() ; + + for(index = 0; index < fRecPoints->GetEntries(); index++) + { + (dynamic_cast(fRecPoints->At(index)))->SetIndexInList(index) ; + (dynamic_cast(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(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(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(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(fRecPoints->At(fNumberOfECAClusters)) ; + fNumberOfECAClusters++ ; + recPoint->SetClusterType(AliESDCaloCluster::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( 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; +} diff --git a/EMCAL/AliEMCALClusterizerNxN.h b/EMCAL/AliEMCALClusterizerNxN.h new file mode 100644 index 00000000000..1d4c8939dff --- /dev/null +++ b/EMCAL/AliEMCALClusterizerNxN.h @@ -0,0 +1,68 @@ +#ifndef ALIEMCALCLUSTERIZERNXN_H +#define ALIEMCALCLUSTERIZERNXN_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id: AliEMCALClusterizerNxN.h 41181 2010-05-12 13:58:06Z gconesab $ */ + +//_________________________________________________________________________ +// 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) + +#include "AliEMCALClusterizer.h" +class AliEMCALRecPoint ; +class AliEMCALDigit ; + +class AliEMCALClusterizerNxN : public AliEMCALClusterizer { + +public: + + AliEMCALClusterizerNxN() ; + AliEMCALClusterizerNxN(AliEMCALGeometry* geometry); + AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * pedestal); + + virtual ~AliEMCALClusterizerNxN() ; + + virtual Int_t AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared)const ; + // Checks if digits are in neighbour cells + + virtual void Digits2Clusters(Option_t *option); // Does the job + + + static Double_t ShowerShape(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; + //class member function (not object member function) + static void SetInputCalibrated(Bool_t val); + + virtual const char * Version() const { return "clu-NxN" ; } + + +protected: + + virtual void MakeClusters(); + + static Bool_t fgkIsInputCalibrated; // to enable reclusterization from ESD cells + +private: + AliEMCALClusterizerNxN(const AliEMCALClusterizerNxN &); //copy ctor + AliEMCALClusterizerNxN & operator = (const AliEMCALClusterizerNxN &); + + + virtual void MakeUnfolding(); + void UnfoldCluster(AliEMCALRecPoint * iniEmc, Int_t Nmax, + AliEMCALDigit ** maxAt, + Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package + +private: + + ClassDef(AliEMCALClusterizerNxN,1) // Clusterizer implementation version 1 + +}; + +#endif // AliEMCALCLUSTERIZERNXN_H diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 2a321ebfbb4..c834a99a8e5 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -16,82 +16,44 @@ /* $Id$ */ //-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) -// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction -// of new IO (à la PHOS) -// Mar 2007, Aleksei Pavlinov - new algoritmh of pseudo clusters +//-- Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class ////////////////////////////////////////////////////////////////////////////// // Clusterization class. Performs clusterization (collects neighbouring active cells) and // unfolds the clusters having several local maxima. // Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints), // EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all // parameters including input digits branch title, thresholds etc.) -// This TTask is normally called from Reconstructioner, but can as well be used in -// standalone mode. -// Use Case: -// root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root") -// Warning in : 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 -class TROOT; -#include #include -class TFolder; #include #include #include -class TSystem; #include #include #include +#include +#include // --- Standard library --- - +#include // --- AliRoot header files --- -#include "AliRunLoader.h" -#include "AliRun.h" -#include "AliESD.h" +#include "AliLog.h" #include "AliEMCALClusterizerv1.h" #include "AliEMCALRecPoint.h" #include "AliEMCALDigit.h" -#include "AliEMCALDigitizer.h" -#include "AliEMCAL.h" #include "AliEMCALGeometry.h" -#include "AliEMCALRecParam.h" -#include "AliEMCALReconstructor.h" -#include "AliCDBManager.h" #include "AliCaloCalibPedestal.h" #include "AliEMCALCalibData.h" -class AliCDBStorage; -#include "AliCDBEntry.h" +#include "AliESDCaloCluster.h" ClassImp(AliEMCALClusterizerv1) //____________________________________________________________________________ -AliEMCALClusterizerv1::AliEMCALClusterizerv1() - : AliEMCALClusterizer(), - fGeom(0), - fDefaultInit(kFALSE), - fToUnfold(kFALSE), - fNumberOfECAClusters(0),fCalibData(0),fCaloPed(0), - fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), - fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.) +AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer() { // ctor with the indication of the file where header Tree and digits Tree are stored @@ -100,45 +62,20 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1() //____________________________________________________________________________ AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry) - : AliEMCALClusterizer(), - fGeom(geometry), - fDefaultInit(kFALSE), - fToUnfold(kFALSE), - fNumberOfECAClusters(0),fCalibData(0), fCaloPed(0), - fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), - fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.) + : 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 - // Note for the future: the use on runloader should be avoided or optional at least - // another way is to make Init virtual and protected at least such that the deriving classes can overload - // Init() ; - // - - if (!fGeom) - { - AliFatal("Geometry not initialized."); - } - } //____________________________________________________________________________ AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped) -: AliEMCALClusterizer(), -fGeom(geometry), -fDefaultInit(kFALSE), -fToUnfold(kFALSE), -fNumberOfECAClusters(0),fCalibData(calib), fCaloPed(caloped), -fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), -fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.) +: AliEMCALClusterizer(geometry, calib, caloped) { // ctor, geometry and calibration are initialized elsewhere. - - if (!fGeom) - AliFatal("Geometry not initialized."); - + } @@ -148,58 +85,6 @@ fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.) // dtor } -//____________________________________________________________________________ -Float_t AliEMCALClusterizerv1::Calibrate(const Float_t amp, const Float_t time, const Int_t absId) -{ - - // Convert digitized amplitude into energy. - // Calibration parameters are taken from calibration data base for raw data, - // or from digitizer parameters for simulated data. - - if(fCalibData){ - - if (fGeom==0) - AliFatal("Did not get geometry from EMCALLoader") ; - - Int_t iSupMod = -1; - Int_t nModule = -1; - Int_t nIphi = -1; - Int_t nIeta = -1; - Int_t iphi = -1; - Int_t ieta = -1; - - Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ; - if(!bCell) { - fGeom->PrintGeometry(); - Error("Calibrate()"," Wrong cell id number : %i", absId); - assert(0); - } - - fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); - - // Check if channel is bad (dead or hot), in this case return 0. - // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation. - // for the moment keep it here but remember to do the selection at the sdigitizer level - // and remove it from here - Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi); - if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) { - AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus)); - return 0; - } - //Check if time is too large or too small, indication of a noisy channel, remove in this case - if(time > fTimeMax || time < fTimeMin) return 0; - - fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); - fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); - - return -fADCpedestalECA + amp * fADCchannelECA ; - - } - else //Return energy with default parameters if calibration is not available - return -fADCpedestalECA + amp * fADCchannelECA ; - -} - //____________________________________________________________________________ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option) { @@ -339,8 +224,8 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit return kFALSE ; } for(index = 0; index < nPar; index++){ - Double_t err ; - Double_t val ; + Double_t err = 0. ; + Double_t val = 0. ; gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index fitparameters[index] = val ; } @@ -350,105 +235,6 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit } -//____________________________________________________________________________ -void AliEMCALClusterizerv1::GetCalibrationParameters() -{ - // Set calibration parameters: - // if calibration database exists, they are read from database, - // otherwise, they are taken from digitizer. - // - // It is a user responsilibity to open CDB before reconstruction, - // for example: - // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); - - //Check if calibration is stored in data base - - if(!fCalibData) - { - AliCDBEntry *entry = (AliCDBEntry*) - AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); - if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); - } - - if(!fCalibData) - AliFatal("Calibration parameters not found in CDB!"); - -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::GetCaloCalibPedestal() -{ - // Set calibration parameters: - // if calibration database exists, they are read from database, - // otherwise, they are taken from digitizer. - // - // It is a user responsilibity to open CDB before reconstruction, - // for example: - // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); - - //Check if calibration is stored in data base - - if(!fCaloPed) - { - AliCDBEntry *entry = (AliCDBEntry*) - AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"); - if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject(); - } - - if(!fCaloPed) - AliFatal("Pedestal info not found in CDB!"); - -} - - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::Init() -{ - // Make all memory allocations which can not be done in default constructor. - // Attach the Clusterizer task to the list of EMCAL tasks - - AliRunLoader *rl = AliRunLoader::Instance(); - if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) - fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); - else - fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); - - AliDebug(1,Form("geom %p",fGeom)); - - if(!gMinuit) - gMinuit = new TMinuit(100) ; - -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::InitParameters() -{ - // Initializes the parameters for the Clusterizer - fNumberOfECAClusters = 0; - - fCalibData = 0 ; - fCaloPed = 0 ; - - const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam(); - if(!recParam) { - AliFatal("Reconstruction parameters for EMCAL not set!"); - } else { - fECAClusteringThreshold = recParam->GetClusteringThreshold(); - fECAW0 = recParam->GetW0(); - fMinECut = recParam->GetMinECut(); - fToUnfold = recParam->GetUnfold(); - if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!"); - fECALocMaxCut = recParam->GetLocMaxCut(); - fTimeCut = recParam->GetTimeCut(); - fTimeMin = recParam->GetTimeMin(); - fTimeMax = recParam->GetTimeMax(); - - AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s", - fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax)); - } - -} - //____________________________________________________________________________ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const { @@ -462,7 +248,6 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d 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, coldiff; shared = kFALSE; @@ -488,8 +273,8 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d }//Different SM, same phi - rowdiff = TMath::Abs(iphi1 - iphi2); - coldiff = TMath::Abs(ieta1 - ieta2) ; + Int_t rowdiff = TMath::Abs(iphi1 - iphi2); + Int_t coldiff = TMath::Abs(ieta1 - ieta2) ; // neighbours with at least common side; May 11, 2007 if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) { @@ -609,7 +394,7 @@ void AliEMCALClusterizerv1::MakeUnfolding() AliEMCALRecPoint * recPoint = dynamic_cast( fRecPoints->At(index) ) ; TVector3 gpos; - Int_t absId; + Int_t absId = -1; recPoint->GetGlobalPosition(gpos); fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId); if(absId > nModulesToUnfold) @@ -686,15 +471,14 @@ void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, AliEMCALDigit * digit = 0 ; Int_t * digitsList = iniTower->GetDigitsList() ; - - Int_t iparam ; + + Int_t iparam = 0 ; Int_t iDigit ; for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ digit = dynamic_cast( fDigitsArr->At(digitsList[iDigit] ) ) ; fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); efit[iDigit] = 0; - iparam = 0 ; while(iparam < nPar ){ xpar = fitparameters[iparam] ; zpar = fitparameters[iparam+1] ; @@ -710,7 +494,7 @@ void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, // to its contribution to efit Float_t * energiesList = iniTower->GetEnergiesList() ; - Float_t ratio ; + Float_t ratio = 0 ; iparam = 0 ; while(iparam < nPar ){ @@ -729,7 +513,7 @@ void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, fNumberOfECAClusters++ ; recPoint->SetNExMax((Int_t)nPar/3) ; - Float_t eDigit ; + Float_t eDigit = 0. ; for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ digit = dynamic_cast( fDigitsArr->At( digitsList[iDigit] ) ) ; fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); @@ -774,7 +558,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, for(iparam = 0 ; iparam < nPar ; iparam++) Grad[iparam] = 0 ; // Will evaluate gradient - Double_t efit ; + Double_t efit = 0. ; AliEMCALDigit * digit ; Int_t iDigit ; @@ -791,7 +575,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, if(iflag == 2){ // calculate gradient Int_t iParam = 0 ; - efit = 0 ; + efit = 0. ; while(iParam < nPar ){ Double_t dx = (xDigit - x[iParam]) ; iParam++ ; @@ -837,86 +621,3 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, // Here we assume, that sigma = sqrt(E) } } -//____________________________________________________________________________ -void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const -{ - // Print clusterizer parameters - - TString message("\n") ; - - if( strcmp(GetName(), "") !=0 ){ - - // Print parameters - - TString taskName(Version()) ; - - printf("--------------- "); - printf("%s",taskName.Data()) ; - printf(" "); - printf("Clusterizing digits: "); - printf("\n ECA Local Maximum cut = %f", fECALocMaxCut); - printf("\n ECA Logarithmic weight = %f", fECAW0); - if(fToUnfold) - printf("\nUnfolding on\n"); - else - printf("\nUnfolding off\n"); - - printf("------------------------------------------------------------------"); - } - else - printf("AliEMCALClusterizerv1 not initialized ") ; -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) -{ - // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer - if(strstr(option,"deb")) { - printf("PrintRecPoints: Clusterization result:") ; - - printf(" Found %d ECA Rec Points\n ", - fRecPoints->GetEntriesFast()) ; - } - - if(strstr(option,"all")) { - if(strstr(option,"deb")) { - printf("\n-----------------------------------------------------------------------\n") ; - printf("Clusters in ECAL section\n") ; - printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; - } - Int_t index =0; - - for (index = 0 ; index < fRecPoints->GetEntries() ; index++) { - AliEMCALRecPoint * rp = dynamic_cast(fRecPoints->At(index)) ; - TVector3 globalpos; - //rp->GetGlobalPosition(globalpos); - TVector3 localpos; - rp->GetLocalPosition(localpos); - Float_t lambda[2]; - rp->GetElipsAxis(lambda); - Int_t * primaries; - Int_t nprimaries; - primaries = rp->GetPrimaries(nprimaries); - if(strstr(option,"deb")) - printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", - rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), - globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), - rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; - if(strstr(option,"deb")){ - for (Int_t iprimary=0; iprimarySetGamma(i,j,0.); @@ -554,6 +557,7 @@ void AliEMCALRecParam::Print(Option_t * opt) const // if nothing is specified print all, if "reco" just clusterization/track matching // if "pid", just PID parameters, if "raw", just raw utils parameters. if(!strcmp("",opt) || !strcmp("reco",opt)){ + AliInfo(Form("Clusterizer selected: %d", fClusterizerFlag)); AliInfo(Form("Clusterization parameters :\n fClusteringThreshold=%.3f,\n fW0=%.3f,\n fMinECut=%.3f,\n fUnfold=%d,\n fLocMaxCut=%.3f,\n fTimeCut=%2.1f ns\n fTimeMin=%2.1f ns\n fTimeMax=%2.1f ns\n", fClusteringThreshold,fW0,fMinECut,fUnfold,fLocMaxCut,fTimeCut*1.e9,fTimeMin*1e9,fTimeMax*1e9)); diff --git a/EMCAL/AliEMCALRecParam.h b/EMCAL/AliEMCALRecParam.h index 48af0574785..3e09c718863 100644 --- a/EMCAL/AliEMCALRecParam.h +++ b/EMCAL/AliEMCALRecParam.h @@ -23,6 +23,12 @@ class AliEMCALRecParam : public AliDetectorRecoParam { public: + + enum AliEMCALClusterizerFlag + { + kClusterizerv1 = 0, + kClusterizerNxN = 1 + }; AliEMCALRecParam() ; AliEMCALRecParam(const AliEMCALRecParam& recParam); @@ -113,16 +119,19 @@ class AliEMCALRecParam : public AliDetectorRecoParam Bool_t UseFALTRO() const {return fUseFALTRO; } Bool_t FitLEDEvents() const {return fFitLEDEvents; } - virtual void Print(Option_t * option="") const ; + virtual void Print(Option_t * option="") const ; static AliEMCALRecParam* GetDefaultParameters(); static AliEMCALRecParam* GetLowFluxParam(); static AliEMCALRecParam* GetHighFluxParam(); static AliEMCALRecParam* GetCalibParam(); static AliEMCALRecParam* GetCosmicParam(); - + static const TObjArray* GetMappings(); + void SetClusterizerFlag(Short_t val) { fClusterizerFlag = val; } + Short_t GetClusterizerFlag() const { return fClusterizerFlag; } + private: //Clustering Float_t fClusteringThreshold ; // Minimum energy to seed a EC digit in a cluster @@ -133,6 +142,7 @@ class AliEMCALRecParam : public AliDetectorRecoParam Float_t fTimeCut ; // Maximum time of digits with respect to EMC cluster max. Float_t fTimeMin ; // Minimum time of digits Float_t fTimeMax ; // Maximum time of digits + Short_t fClusterizerFlag ; // Choice of the clusterizer; Default selection (v1) is zero //PID (Guenole) Double_t fGamma[6][6]; // Parameter to Compute PID for photons @@ -169,7 +179,7 @@ class AliEMCALRecParam : public AliDetectorRecoParam static TObjArray* fgkMaps; // ALTRO mappings for RCU0..RCUX - ClassDef(AliEMCALRecParam,12) // Reconstruction parameters + ClassDef(AliEMCALRecParam,13) // Reconstruction parameters } ; diff --git a/EMCAL/AliEMCALReconstructor.cxx b/EMCAL/AliEMCALReconstructor.cxx index 63094c31a38..048721c14f8 100644 --- a/EMCAL/AliEMCALReconstructor.cxx +++ b/EMCAL/AliEMCALReconstructor.cxx @@ -45,6 +45,7 @@ #include "AliEMCALRawUtils.h" #include "AliEMCALDigit.h" #include "AliEMCALClusterizerv1.h" +#include "AliEMCALClusterizerNxN.h" #include "AliEMCALRecPoint.h" #include "AliEMCALPID.h" #include "AliEMCALTrigger.h" @@ -106,13 +107,11 @@ AliEMCALReconstructor::AliEMCALReconstructor() if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject(); } - if(!fPedestalData) - AliFatal("Dead map not found in CDB!"); - - - //Init the clusterizer with geometry and calibration pointers, avoid doing it twice. - fgClusterizer = new AliEMCALClusterizerv1(fGeom, fCalibData,fPedestalData); + if(!fPedestalData) + AliFatal("Dead map not found in CDB!"); + InitClusterizer(); + if(!fGeom) AliFatal(Form("Could not get geometry!")); AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance(); @@ -142,6 +141,32 @@ AliEMCALReconstructor::~AliEMCALReconstructor() // fList = AliEMCALHistoUtilities::GetTriggersListOfHists(kTRUE); // } +//____________________________________________________________________________ +void AliEMCALReconstructor::InitClusterizer() +{ + //Init the clusterizer with geometry and calibration pointers, avoid doing it twice. + + AliEMCALRecParam *recParam = NULL; + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/RecoParam"); + //Get The reco param for the default event specie + if (entry) + recParam = (AliEMCALRecParam*)((TObjArray *) entry->GetObject())->At(0); + + if(!recParam) + AliFatal("RecoParam not found in CDB!"); + + if (recParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1) + { + fgClusterizer = new AliEMCALClusterizerv1(fGeom, fCalibData,fPedestalData); + } + else + { + fgClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData,fPedestalData); + } + +} + //____________________________________________________________________________ void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const { @@ -154,6 +179,7 @@ void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) AliCodeTimerAuto("",0) ReadDigitsArrayFromTree(digitsTree); + fgClusterizer->InitParameters(); fgClusterizer->SetOutput(clustersTree); @@ -287,7 +313,7 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, for (Int_t idig = 0 ; idig < nDigits ; idig++) { const AliEMCALDigit * dig = (const AliEMCALDigit*)digits->At(idig); if(dig->GetAmplitude() > 0 ){ - energy = (static_cast (fgClusterizer))->Calibrate(dig->GetAmplitude(),dig->GetTime(),dig->GetId()); //TimeR or Time? + energy = fgClusterizer->Calibrate(dig->GetAmplitude(),dig->GetTime(),dig->GetId()); //TimeR or Time? if(energy > 0){ //Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them emcCells.SetCell(idignew,dig->GetId(),energy, dig->GetTime()); idignew++; @@ -381,13 +407,13 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, ec->SetPosition(xyz); ec->SetE(clust->GetEnergy()); - //Distance to the nearest bad crystal - ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower()); + //Distance to the nearest bad crystal + ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower()); ec->SetNCells(newCellMult); //Change type of list from short to ushort UShort_t *newAbsIdList = new UShort_t[newCellMult]; - Double_t *newFracList = new Double_t[newCellMult]; + Double_t *newFracList = new Double_t[newCellMult]; for(Int_t i = 0; i < newCellMult ; i++) { newAbsIdList[i]=absIdList[i]; newFracList[i]=fracList[i]; diff --git a/EMCAL/AliEMCALReconstructor.h b/EMCAL/AliEMCALReconstructor.h index 9b600752baf..73b6f722394 100644 --- a/EMCAL/AliEMCALReconstructor.h +++ b/EMCAL/AliEMCALReconstructor.h @@ -54,7 +54,8 @@ public: virtual ~AliEMCALReconstructor() ; //dtor virtual void Init() {;} - + virtual void InitClusterizer(); + Bool_t Debug() const { return fDebug ; } using AliReconstructor::FillESD; @@ -63,7 +64,7 @@ public: AliTracker* CreateTracker () const {return new AliEMCALTracker;} using AliReconstructor::Reconstruct; - virtual void Reconstruct(TTree* digitsTree, TTree* clustersTree) const; + virtual void Reconstruct(TTree* digitsTree, TTree* clustersTree) const ; virtual Bool_t HasDigitConversion() const {return kTRUE;}; virtual void ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const; diff --git a/EMCAL/EMCALrecLinkDef.h b/EMCAL/EMCALrecLinkDef.h index 6979aaf889f..2053121a011 100644 --- a/EMCAL/EMCALrecLinkDef.h +++ b/EMCAL/EMCALrecLinkDef.h @@ -6,6 +6,7 @@ #pragma link C++ class AliEMCALReconstructor+; #pragma link C++ class AliEMCALClusterizer+; #pragma link C++ class AliEMCALClusterizerv1+; +#pragma link C++ class AliEMCALClusterizerNxN+; #pragma link C++ class AliEMCALTrack+; #pragma link C++ class AliEMCALTracker+; #pragma link C++ class AliEMCALPID+; diff --git a/EMCAL/libEMCALrec.pkg b/EMCAL/libEMCALrec.pkg index 6ce12bf7522..f1e92587971 100644 --- a/EMCAL/libEMCALrec.pkg +++ b/EMCAL/libEMCALrec.pkg @@ -5,6 +5,7 @@ SRCS =\ AliEMCALReconstructor.cxx \ AliEMCALClusterizer.cxx \ AliEMCALClusterizerv1.cxx \ +AliEMCALClusterizerNxN.cxx \ AliEMCALTrack.cxx \ AliEMCALTracker.cxx \ AliEMCALPID.cxx \ @@ -20,8 +21,6 @@ DHDR:=EMCALrecLinkDef.h # have to tune EINCLUDE:=PYTHIA6 RAW EVGEN THijing VZERO -EXPORT:=AliEMCALTrack.h - ifeq (win32gcc,$(ALICE_TARGET)) PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) \ -lEMCALUtils -lEMCALbase -lEMCALsim -lSTEER -lESD -lCDB -lSTEERBase \ diff --git a/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C b/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C index c9229c6870f..d1f3b54cda6 100644 --- a/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C +++ b/EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C @@ -127,6 +127,8 @@ AliEMCALRecParam* GetHighMultiplicityParameters() //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters(); AliEMCALRecParam* params = AliEMCALRecParam::GetDefaultParameters(); + //params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN); + params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1); //Clusterization params->SetClusteringThreshold(0.5); @@ -306,6 +308,9 @@ AliEMCALRecParam* GetLowMultiplicityParameters() //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetLowFluxParam(); AliEMCALRecParam* params = AliEMCALRecParam::GetDefaultParameters(); + //params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN); + params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1); + params->SetClusteringThreshold(0.1); // 100 MeV params->SetMinECut(0.01); //10 MeV params->SetTimeCut(250e-9);//250 ns diff --git a/EMCAL/macros/RecParamDB/PrintEMCALRecParam.C b/EMCAL/macros/RecParamDB/PrintEMCALRecParam.C index 00176fe4277..fe225feaa26 100644 --- a/EMCAL/macros/RecParamDB/PrintEMCALRecParam.C +++ b/EMCAL/macros/RecParamDB/PrintEMCALRecParam.C @@ -4,7 +4,7 @@ // Author: Gustavo Conesa (INFN-LNF) -void PrintEMCALRecParam(char * file = "Run0_999999999_v0_s0.root"){ +void PrintEMCALRecParam(char * file = "$ALICE_ROOT/OCDB/EMCAL/Calib/RecoParam/Run0_999999999_v0_s0.root"){ TFile * f = new TFile(file,"READ"); @@ -24,9 +24,10 @@ cout<<"============== "<GetName()<<" ==============="<Print(""); - +//rparam->Print("reco");//Print only clusterizer parameters +//rparam->Print("pid");//Print only pid parameters +//rparam->Print("raw");//Print only raw digitization parameters +rparam->Print("");// Print all }