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()
154 if (fClusterUnfolding) delete fClusterUnfolding;
157 fRecPoints->Delete();
161 fDigitsArr->Clear("C");
166 //------------------------------------------------------------------------
167 void AliEMCALAfterBurnerUF::Clear()
171 if (fRecPoints) fRecPoints->Delete(); // do not Clear(), it leaks, why?
172 if (fDigitsArr) fDigitsArr->Clear("C");
175 //------------------------------------------------------------------------
176 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
178 // Restore clusters from recPoints
179 // Cluster energy, global position, cells and their amplitude fractions are restored
181 // TODO: restore time and other parameters
183 for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
185 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
187 const Int_t ncells = recPoint->GetMultiplicity();
188 Int_t ncellsTrue = 0;
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];
194 for (Int_t c = 0; c < ncells; c++) {
195 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
197 absIds[ncellsTrue] = digit->GetId();
198 ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
200 if (ratios[ncellsTrue] > 0.001) ncellsTrue++;
203 if (ncellsTrue < 1) {
204 Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
211 // calculate new cluster position
212 recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
213 recPoint->GetGlobalPosition(gpos);
216 // create a new cluster
217 AliAODCaloCluster *clus = new AliAODCaloCluster();
218 clus->SetType(AliVCluster::kEMCALClusterv1);
219 clus->SetE(recPoint->GetEnergy());
220 clus->SetPosition(g);
221 clus->SetNCells(ncellsTrue);
222 clus->SetCellsAbsId(absIds);
223 clus->SetCellsAmplitudeFraction(ratios);
224 // TODO: time not stored
225 // TODO: some other properties not stored
227 clusArray->Add(clus);
232 //------------------------------------------------------------------------
233 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
237 // Input: TObjArray of clusters, EMCAL cells.
238 // Output is added to the same array, original clusters are _deleted_ or moved to another position.
242 Int_t nclus = clusArray->GetEntriesFast();
244 /* Fill recPoints with digits
246 for (Int_t i = 0; i < nclus; i++)
248 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
249 if (!clus->IsEMCAL()) continue;
252 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
253 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
254 fRecPoints->Add(recPoint);
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);
262 // NOTE: it is easy to include cells recalibration here:
265 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
268 digit->SetAmplitude(amp);
269 digit->SetTime(time);
270 digit->SetTimeR(time);
271 digit->SetIndexInList(ndigits);
273 recPoint->AddDigit(*digit, amp, kFALSE);
278 // this cluster will be substituted with the result of unfolding
279 clusArray->RemoveAt(i);
286 fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
287 fClusterUnfolding->MakeUnfolding();
289 /* Restore clusters from recPoints.
291 RecPoints2Clusters(clusArray);
294 fRecPoints->Delete(); // do not Clear(), it leaks, why?
295 fDigitsArr->Clear("C");
297 clusArray->Compress();