1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 // After-burner for the EMCAL cluster unfolding algorithm
17 // Input: TObjArray *clusArray -- array of AliVClusters;
18 // AliVCaloCells *cellsEMCAL -- EMCAL cells.
20 // Output is appended to clusArray, the original (unfolded or not) clusters
21 // are deleted or moved to another position in clusArray.
23 // If you want to use particular geometry, you must initialize it _before_
24 // creating AliEMCALAfterBurnerUF instance. Add this or similar line to the
25 // initialization section:
27 // AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
29 // gGeoManager must be initialized for this code to work! Do it yourself or
30 // provide geometry.root file in the current directory so that
31 // AliEMCALAfterBurnerUF will take it by itself.
34 // // add this lines to the initialization section of your analysis
35 // AliEMCALAfterBurnerUF *abuf = new AliEMCALAfterBurnerUF();
36 // TObjArray *clusArray = new TObjArray(100);
39 // AliVEvent *event = InputEvent();
40 // AliVCaloCells *cellsEMCAL = event->GetEMCALCells();
42 // for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
44 // AliVCluster *clus = event->GetCaloCluster(i);
46 // clusArray->Add(clus->Clone()); // NOTE _CLONE_ in this line
49 // abuf->UnfoldClusters(clusArray, cellsEMCAL);
51 // // do an analysis with clusArray
54 // // prevent memory leak
55 // clusArray->Delete();
58 // Author: Olga Driga (SUBATECH)
60 // --- ROOT system ---
61 #include <TObjArray.h>
62 #include <TClonesArray.h>
63 #include <TGeoManager.h>
65 // --- Standard library --
67 // --- AliRoot header files ---
68 #include <AliEMCALAfterBurnerUF.h>
69 #include <AliEMCALGeometry.h>
70 #include <AliEMCALUnfolding.h>
71 #include <AliAODCaloCluster.h>
72 #include <AliVCaloCells.h>
73 #include <AliEMCALRecPoint.h>
74 #include <AliEMCALDigit.h>
76 ClassImp(AliEMCALAfterBurnerUF)
78 //------------------------------------------------------------------------
79 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
81 fLogWeight(4.5), // correct?
82 fECALocMaxCut(0.03), // value suggested by Adam
86 fClusterUnfolding(NULL)
88 // Use this constructor, if unsure
93 //------------------------------------------------------------------------
94 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ecaLocMaxCut, Float_t minECut):
96 fLogWeight(logWeight),
97 fECALocMaxCut(ecaLocMaxCut),
101 fClusterUnfolding(NULL)
103 // This constructor allows to set parameters
104 // Recommended values:
105 // Float_t logWeight = 4.5, ECALocMaxCut = 0.03
110 //------------------------------------------------------------------------
111 void AliEMCALAfterBurnerUF::Init()
113 // After-burner initialization
114 // Imports geometry.root (if required), creates unfolding class instance
116 // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
117 // for this to work, the OCDB geometry must be imported instead
120 Warning("AliEMCALAfterBurnerUF::Init","GeoManager not found, please import the geometry.root file or pass to the geometry the misalignment matrices");
121 // TGeoManager::Import("geometry.root");
123 // required for global cluster position recalculation
125 Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
127 // initialize geometry, if not yet initialized
128 if (!AliEMCALGeometry::GetInstance()) {
129 Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETEV1");
130 AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
133 // AliEMCALRecPoint is using exactly this call
134 fGeom = AliEMCALGeometry::GetInstance();
136 Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
138 fClusterUnfolding = new AliEMCALUnfolding(fGeom);
139 fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
140 fClusterUnfolding->SetThreshold(fMinECut);
142 // clusters --> recPoints, cells --> digits and back ;)
143 fRecPoints = new TObjArray(100);
144 fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
147 //------------------------------------------------------------------------
148 AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
150 if (fClusterUnfolding) delete fClusterUnfolding;
153 fRecPoints->Delete();
157 fDigitsArr->Clear("C");
162 //------------------------------------------------------------------------
163 void AliEMCALAfterBurnerUF::Clear()
167 if (fRecPoints) fRecPoints->Delete(); // do not Clear(), it leaks, why?
168 if (fDigitsArr) fDigitsArr->Clear("C");
171 //------------------------------------------------------------------------
172 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
174 // Restore clusters from recPoints
175 // Cluster energy, global position, cells and their amplitude fractions are restored
177 // TODO: restore time and other parameters
179 for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
181 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
183 const Int_t ncells = recPoint->GetMultiplicity();
184 Int_t ncells_true = 0;
186 // cells and their amplitude fractions
187 UShort_t absIds[ncells]; // NOTE: unfolding must not give recPoints with no cells!
188 Double32_t ratios[ncells];
190 for (Int_t c = 0; c < ncells; c++) {
191 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
193 absIds[ncells_true] = digit->GetId();
194 ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
196 if (ratios[ncells_true] > 0.001) ncells_true++;
199 if (ncells_true < 1) {
200 Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
207 // calculate new cluster position
208 recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
209 recPoint->GetGlobalPosition(gpos);
212 // create a new cluster
213 AliAODCaloCluster *clus = new AliAODCaloCluster();
214 clus->SetType(AliVCluster::kEMCALClusterv1);
215 clus->SetE(recPoint->GetEnergy());
216 clus->SetPosition(g);
217 clus->SetNCells(ncells_true);
218 clus->SetCellsAbsId(absIds);
219 clus->SetCellsAmplitudeFraction(ratios);
220 // TODO: time not stored
221 // TODO: some other properties not stored
223 clusArray->Add(clus);
228 //------------------------------------------------------------------------
229 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
233 // Input: TObjArray of clusters, EMCAL cells.
234 // Output is added to the same array, original clusters are _deleted_ or moved to another position.
238 Int_t nclus = clusArray->GetEntriesFast();
240 /* Fill recPoints with digits
242 for (Int_t i = 0; i < nclus; i++)
244 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
245 if (!clus->IsEMCAL()) continue;
248 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
249 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
250 fRecPoints->Add(recPoint);
253 for (Int_t c = 0; c < clus->GetNCells(); c++) {
254 Int_t absId = clus->GetCellAbsId(c);
255 Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
256 Double_t time = cellsEMCAL->GetCellTime(absId);
258 // NOTE: it is easy to include cells recalibration here:
261 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
264 digit->SetAmplitude(amp);
265 digit->SetTime(time);
266 digit->SetTimeR(time);
267 digit->SetIndexInList(ndigits);
269 recPoint->AddDigit(*digit, amp, kFALSE);
274 // this cluster will be substituted with the result of unfolding
275 clusArray->RemoveAt(i);
282 fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
283 fClusterUnfolding->MakeUnfolding();
285 /* Restore clusters from recPoints.
287 RecPoints2Clusters(clusArray);
290 fRecPoints->Delete(); // do not Clear(), it leaks, why?
291 fDigitsArr->Clear("C");
293 clusArray->Compress();