Adding include path to allow compilation of CleanGeom task
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALAfterBurnerUF.cxx
CommitLineData
fc645679 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15// After-burner for the EMCAL cluster unfolding algorithm
16//
15cb4c51 17// Input: TObjArray *clusArray -- array of AliVClusters;
18// AliVCaloCells *cellsEMCAL -- EMCAL cells.
fc645679 19//
20// Output is appended to clusArray, the original (unfolded or not) clusters
21// are deleted or moved to another position in clusArray.
22//
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:
26//
27// AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
28//
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.
32// How to use:
33//
34// // add this lines to the initialization section of your analysis
35// AliEMCALAfterBurnerUF *abuf = new AliEMCALAfterBurnerUF();
36// TObjArray *clusArray = new TObjArray(100);
37//
38//
004d0978 39// AliVEvent *event = InputEvent();
40// AliVCaloCells *cellsEMCAL = event->GetEMCALCells();
fc645679 41//
42// for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
43// {
004d0978 44// AliVCluster *clus = event->GetCaloCluster(i);
fc645679 45//
46// clusArray->Add(clus->Clone()); // NOTE _CLONE_ in this line
47// }
48//
49// abuf->UnfoldClusters(clusArray, cellsEMCAL);
50//
51// // do an analysis with clusArray
52// // ....
53//
54// // prevent memory leak
55// clusArray->Delete();
56//
57//----
58// Author: Olga Driga (SUBATECH)
59
60// --- ROOT system ---
61#include <TObjArray.h>
62#include <TClonesArray.h>
63#include <TGeoManager.h>
64
65// --- Standard library --
66
67// --- AliRoot header files ---
68#include <AliEMCALAfterBurnerUF.h>
69#include <AliEMCALGeometry.h>
70#include <AliEMCALUnfolding.h>
71#include <AliAODCaloCluster.h>
15cb4c51 72#include <AliVCaloCells.h>
fc645679 73#include <AliEMCALRecPoint.h>
74#include <AliEMCALDigit.h>
75
76ClassImp(AliEMCALAfterBurnerUF)
77
78//------------------------------------------------------------------------
79AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
80 fGeom(NULL),
81 fLogWeight(4.5), // correct?
82 fECALocMaxCut(0.03), // value suggested by Adam
71332f0e 83 fMinECut(0.01),
fc645679 84 fRecPoints(NULL),
85 fDigitsArr(NULL),
86 fClusterUnfolding(NULL)
87{
88 // Use this constructor, if unsure
89
90 Init();
91}
92
93//------------------------------------------------------------------------
71332f0e 94AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ecaLocMaxCut, Float_t minECut):
fc645679 95 fGeom(NULL),
96 fLogWeight(logWeight),
71332f0e 97 fECALocMaxCut(ecaLocMaxCut),
98 fMinECut(minECut),
fc645679 99 fRecPoints(NULL),
100 fDigitsArr(NULL),
101 fClusterUnfolding(NULL)
102{
103 // This constructor allows to set parameters
104 // Recommended values:
105 // Float_t logWeight = 4.5, ECALocMaxCut = 0.03
106
107 Init();
108}
109
110//------------------------------------------------------------------------
111void AliEMCALAfterBurnerUF::Init()
112{
113 // After-burner initialization
114 // Imports geometry.root (if required), creates unfolding class instance
115 //
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
118
119 if (!gGeoManager)
4b7fc3d0 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");
fc645679 122
123 // required for global cluster position recalculation
124 if (!gGeoManager)
0d0d6b98 125 Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
fc645679 126
127 // initialize geometry, if not yet initialized
128 if (!AliEMCALGeometry::GetInstance()) {
0d0d6b98 129 Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETEV1");
130 AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
fc645679 131 }
132
133 // AliEMCALRecPoint is using exactly this call
134 fGeom = AliEMCALGeometry::GetInstance();
135 if (!fGeom)
136 Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
137
138 fClusterUnfolding = new AliEMCALUnfolding(fGeom);
139 fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
71332f0e 140 fClusterUnfolding->SetThreshold(fMinECut);
141
fc645679 142 // clusters --> recPoints, cells --> digits and back ;)
143 fRecPoints = new TObjArray(100);
144 fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
145}
146
147//------------------------------------------------------------------------
148AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
149{
1186cd2b 150 //
151 // destructor
152 //
153
fc645679 154 if (fClusterUnfolding) delete fClusterUnfolding;
155
0d0d6b98 156 if (fRecPoints) {
157 fRecPoints->Delete();
158 delete fRecPoints;
159 }
4b7fc3d0 160 if (fDigitsArr) {
161 fDigitsArr->Clear("C");
162 delete fDigitsArr;
163 }
fc645679 164}
165
166//------------------------------------------------------------------------
0d0d6b98 167void AliEMCALAfterBurnerUF::Clear()
168{
169 //Clean the arrays
170
171 if (fRecPoints) fRecPoints->Delete(); // do not Clear(), it leaks, why?
172 if (fDigitsArr) fDigitsArr->Clear("C");
173
174}
175//------------------------------------------------------------------------
fc645679 176void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
177{
178 // Restore clusters from recPoints
179 // Cluster energy, global position, cells and their amplitude fractions are restored
180 //
181 // TODO: restore time and other parameters
182
183 for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
184 {
185 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
186
03b08eaf 187 const Int_t ncells = recPoint->GetMultiplicity();
1186cd2b 188 Int_t ncellsTrue = 0;
fc645679 189
190 // cells and their amplitude fractions
191 UShort_t absIds[ncells]; // NOTE: unfolding must not give recPoints with no cells!
192 Double32_t ratios[ncells];
193
194 for (Int_t c = 0; c < ncells; c++) {
195 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
196
1186cd2b 197 absIds[ncellsTrue] = digit->GetId();
198 ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
fc645679 199
1186cd2b 200 if (ratios[ncellsTrue] > 0.001) ncellsTrue++;
fc645679 201 }
202
1186cd2b 203 if (ncellsTrue < 1) {
fc645679 204 Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
205 continue;
206 }
207
208 TVector3 gpos;
209 Float_t g[3];
210
211 // calculate new cluster position
212 recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
213 recPoint->GetGlobalPosition(gpos);
214 gpos.GetXYZ(g);
215
216 // create a new cluster
217 AliAODCaloCluster *clus = new AliAODCaloCluster();
15cb4c51 218 clus->SetType(AliVCluster::kEMCALClusterv1);
fc645679 219 clus->SetE(recPoint->GetEnergy());
220 clus->SetPosition(g);
1186cd2b 221 clus->SetNCells(ncellsTrue);
fc645679 222 clus->SetCellsAbsId(absIds);
223 clus->SetCellsAmplitudeFraction(ratios);
224 // TODO: time not stored
225 // TODO: some other properties not stored
226
227 clusArray->Add(clus);
228 } // recPoints loop
229
230}
231
232//------------------------------------------------------------------------
004d0978 233void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
fc645679 234{
235 // Unfolds clusters.
236 //
237 // Input: TObjArray of clusters, EMCAL cells.
238 // Output is added to the same array, original clusters are _deleted_ or moved to another position.
239
240 Int_t ndigits = 0;
241
242 Int_t nclus = clusArray->GetEntriesFast();
243
244 /* Fill recPoints with digits
245 */
246 for (Int_t i = 0; i < nclus; i++)
247 {
004d0978 248 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
fc645679 249 if (!clus->IsEMCAL()) continue;
250
251 // new recPoint
252 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
15cb4c51 253 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
fc645679 254 fRecPoints->Add(recPoint);
255
256 // fill digits
257 for (Int_t c = 0; c < clus->GetNCells(); c++) {
258 Int_t absId = clus->GetCellAbsId(c);
259 Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
260 Double_t time = cellsEMCAL->GetCellTime(absId);
261
262 // NOTE: it is easy to include cells recalibration here:
263 // amp *= factor;
264
265 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
266
267 digit->SetId(absId);
268 digit->SetAmplitude(amp);
269 digit->SetTime(time);
270 digit->SetTimeR(time);
271 digit->SetIndexInList(ndigits);
272
273 recPoint->AddDigit(*digit, amp, kFALSE);
274
275 ndigits++;
276 }
277
278 // this cluster will be substituted with the result of unfolding
279 clusArray->RemoveAt(i);
280 delete clus;
281 } // cluster loop
282
283
284 /* Peform unfolding
285 */
286 fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
287 fClusterUnfolding->MakeUnfolding();
288
289 /* Restore clusters from recPoints.
290 */
291 RecPoints2Clusters(clusArray);
292
293 // clean up
0d0d6b98 294 fRecPoints->Delete(); // do not Clear(), it leaks, why?
4b7fc3d0 295 fDigitsArr->Clear("C");
fc645679 296
297 clusArray->Compress();
298
299}