--- /dev/null
+/**************************************************************************
+ * 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. *
+ **************************************************************************/
+// After-burner for the EMCAL cluster unfolding algorithm
+//
+// Input: TObjArray *clusArray -- array of AliAODCaloClusters;
+// AliAODCaloCells *cellsEMCAL -- EMCAL cells.
+//
+// Output is appended to clusArray, the original (unfolded or not) clusters
+// are deleted or moved to another position in clusArray.
+//
+// If you want to use particular geometry, you must initialize it _before_
+// creating AliEMCALAfterBurnerUF instance. Add this or similar line to the
+// initialization section:
+//
+// AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
+//
+// gGeoManager must be initialized for this code to work! Do it yourself or
+// provide geometry.root file in the current directory so that
+// AliEMCALAfterBurnerUF will take it by itself.
+// How to use:
+//
+// // add this lines to the initialization section of your analysis
+// AliEMCALAfterBurnerUF *abuf = new AliEMCALAfterBurnerUF();
+// TObjArray *clusArray = new TObjArray(100);
+//
+//
+// AliAODEvent *event = (AliAODEvent*) InputEvent();
+// AliAODCaloCells *cellsEMCAL = event->GetEMCALCells();
+//
+// for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
+// {
+// AliAODCaloCluster *clus = event->GetCaloCluster(i);
+//
+// clusArray->Add(clus->Clone()); // NOTE _CLONE_ in this line
+// }
+//
+// abuf->UnfoldClusters(clusArray, cellsEMCAL);
+//
+// // do an analysis with clusArray
+// // ....
+//
+// // prevent memory leak
+// clusArray->Delete();
+//
+//----
+// Author: Olga Driga (SUBATECH)
+
+// --- ROOT system ---
+#include <TObjArray.h>
+#include <TClonesArray.h>
+#include <TGeoManager.h>
+
+// --- Standard library --
+
+// --- AliRoot header files ---
+#include <AliEMCALAfterBurnerUF.h>
+#include <AliEMCALGeometry.h>
+#include <AliEMCALUnfolding.h>
+#include <AliAODCaloCluster.h>
+#include <AliAODCaloCells.h>
+#include <AliEMCALRecPoint.h>
+#include <AliEMCALDigit.h>
+
+ClassImp(AliEMCALAfterBurnerUF)
+
+//------------------------------------------------------------------------
+AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
+ fGeom(NULL),
+ fLogWeight(4.5), // correct?
+ fECALocMaxCut(0.03), // value suggested by Adam
+ fRecPoints(NULL),
+ fDigitsArr(NULL),
+ fClusterUnfolding(NULL)
+{
+ // Use this constructor, if unsure
+
+ Init();
+}
+
+//------------------------------------------------------------------------
+AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut):
+ fGeom(NULL),
+ fLogWeight(logWeight),
+ fECALocMaxCut(ECALocMaxCut),
+ fRecPoints(NULL),
+ fDigitsArr(NULL),
+ fClusterUnfolding(NULL)
+{
+ // This constructor allows to set parameters
+ // Recommended values:
+ // Float_t logWeight = 4.5, ECALocMaxCut = 0.03
+
+ Init();
+}
+
+//------------------------------------------------------------------------
+void AliEMCALAfterBurnerUF::Init()
+{
+ // After-burner initialization
+ // Imports geometry.root (if required), creates unfolding class instance
+ //
+ // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
+ // for this to work, the OCDB geometry must be imported instead
+
+ if (!gGeoManager)
+ TGeoManager::Import("geometry.root");
+
+ // required for global cluster position recalculation
+ if (!gGeoManager)
+ Fatal("AliEMCALAfterBurnerUF::Init", "failed to import geometry.root");
+
+ // initialize geometry, if not yet initialized
+ if (!AliEMCALGeometry::GetInstance()) {
+ Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETE");
+ AliEMCALGeometry::GetInstance("EMCAL_COMPLETE");
+ }
+
+ // AliEMCALRecPoint is using exactly this call
+ fGeom = AliEMCALGeometry::GetInstance();
+ if (!fGeom)
+ Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
+
+ fClusterUnfolding = new AliEMCALUnfolding(fGeom);
+ fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
+
+ // clusters --> recPoints, cells --> digits and back ;)
+ fRecPoints = new TObjArray(100);
+ fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
+}
+
+//------------------------------------------------------------------------
+AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
+{
+ if (fClusterUnfolding) delete fClusterUnfolding;
+
+ if (fRecPoints) delete fRecPoints;
+ if (fDigitsArr) delete fDigitsArr;
+}
+
+//------------------------------------------------------------------------
+void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
+{
+ // Restore clusters from recPoints
+ // Cluster energy, global position, cells and their amplitude fractions are restored
+ //
+ // TODO: restore time and other parameters
+
+ for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
+ {
+ AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
+
+ Int_t ncells = recPoint->GetMultiplicity();
+ Int_t ncells_true = 0;
+
+ // cells and their amplitude fractions
+ UShort_t absIds[ncells]; // NOTE: unfolding must not give recPoints with no cells!
+ Double32_t ratios[ncells];
+
+ for (Int_t c = 0; c < ncells; c++) {
+ AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
+
+ absIds[ncells_true] = digit->GetId();
+ ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
+
+ if (ratios[ncells_true] > 0.001) ncells_true++;
+ }
+
+ if (ncells_true < 1) {
+ Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
+ continue;
+ }
+
+ TVector3 gpos;
+ Float_t g[3];
+
+ // calculate new cluster position
+ recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
+ recPoint->GetGlobalPosition(gpos);
+ gpos.GetXYZ(g);
+
+ // create a new cluster
+ AliAODCaloCluster *clus = new AliAODCaloCluster();
+ clus->SetType(AliAODCaloCluster::kEMCALClusterv1);
+ clus->SetE(recPoint->GetEnergy());
+ clus->SetPosition(g);
+ clus->SetNCells(ncells_true);
+ clus->SetCellsAbsId(absIds);
+ clus->SetCellsAmplitudeFraction(ratios);
+ // TODO: time not stored
+ // TODO: some other properties not stored
+
+ clusArray->Add(clus);
+ } // recPoints loop
+
+}
+
+//------------------------------------------------------------------------
+void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliAODCaloCells *cellsEMCAL)
+{
+ // Unfolds clusters.
+ //
+ // Input: TObjArray of clusters, EMCAL cells.
+ // Output is added to the same array, original clusters are _deleted_ or moved to another position.
+
+ Int_t ndigits = 0;
+
+ Int_t nclus = clusArray->GetEntriesFast();
+
+ /* Fill recPoints with digits
+ */
+ for (Int_t i = 0; i < nclus; i++)
+ {
+ AliAODCaloCluster *clus = (AliAODCaloCluster*) clusArray->At(i);
+ if (!clus->IsEMCAL()) continue;
+
+ // new recPoint
+ AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
+ recPoint->SetClusterType(AliAODCaloCluster::kEMCALClusterv1);
+ fRecPoints->Add(recPoint);
+
+ // fill digits
+ for (Int_t c = 0; c < clus->GetNCells(); c++) {
+ Int_t absId = clus->GetCellAbsId(c);
+ Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
+ Double_t time = cellsEMCAL->GetCellTime(absId);
+
+ // NOTE: it is easy to include cells recalibration here:
+ // amp *= factor;
+
+ AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
+
+ digit->SetId(absId);
+ digit->SetAmplitude(amp);
+ digit->SetTime(time);
+ digit->SetTimeR(time);
+ digit->SetIndexInList(ndigits);
+
+ recPoint->AddDigit(*digit, amp, kFALSE);
+
+ ndigits++;
+ }
+
+ // this cluster will be substituted with the result of unfolding
+ clusArray->RemoveAt(i);
+ delete clus;
+ } // cluster loop
+
+
+ /* Peform unfolding
+ */
+ fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
+ fClusterUnfolding->MakeUnfolding();
+
+ /* Restore clusters from recPoints.
+ */
+ RecPoints2Clusters(clusArray);
+
+ // clean up
+ fRecPoints->Delete();
+ fDigitsArr->Delete();
+
+ clusArray->Compress();
+
+}
--- /dev/null
+#ifndef ALIEMCALAFTERBURNERUF_H
+#define ALIEMCALAFTERBURNERUF_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//_________________________________________________________________________
+// After-burner for the EMCAL cluster unfolding algorithm
+//
+// See cxx for details on how to use it
+//
+// Author: Olga Driga (SUBATECH)
+//
+
+// --- ROOT system ---
+class TObjArray;
+class TClonesArray;
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+class AliEMCALGeometry;
+class AliEMCALUnfolding;
+class AliAODCaloCells;
+
+class AliEMCALAfterBurnerUF {
+
+ public:
+ AliEMCALAfterBurnerUF();
+ AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut);
+ virtual ~AliEMCALAfterBurnerUF();
+
+ private:
+ AliEMCALAfterBurnerUF(const AliEMCALAfterBurnerUF & uf) ; // cpy ctor not needed, put here to avoid compilation warning
+ AliEMCALAfterBurnerUF & operator = (const AliEMCALAfterBurnerUF & uf) ;//cpy assignment, put here to avoid compilation warning
+
+ public:
+ virtual void Init();
+ virtual void RecPoints2Clusters(TObjArray *clusArray);
+ virtual void UnfoldClusters(TObjArray *clusArray, AliAODCaloCells *cellsEMCAL); // does the job
+
+ // getters and setters
+ virtual AliEMCALUnfolding *GetClusterUnfoldingInstance() { return fClusterUnfolding; }
+
+ protected:
+ AliEMCALGeometry *fGeom; // EMCAL geometry
+ Float_t fLogWeight; // used in AliEMCALRecPoint::EvalGlobalPosition()
+ Float_t fECALocMaxCut; // this amount of energy must distinguish a local maximum from its neighbours
+ TObjArray *fRecPoints; // cluster <=> recPoint
+ TClonesArray *fDigitsArr; // cell <=> digit
+
+ AliEMCALUnfolding *fClusterUnfolding; // unfolding class instance
+
+ ClassDef(AliEMCALAfterBurnerUF,1)
+} ;
+
+#endif // AliEMCALAFTERBURNERUF_H