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
85 fClusterUnfolding(NULL)
87 // Use this constructor, if unsure
92 //------------------------------------------------------------------------
93 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut):
95 fLogWeight(logWeight),
96 fECALocMaxCut(ECALocMaxCut),
99 fClusterUnfolding(NULL)
101 // This constructor allows to set parameters
102 // Recommended values:
103 // Float_t logWeight = 4.5, ECALocMaxCut = 0.03
108 //------------------------------------------------------------------------
109 void AliEMCALAfterBurnerUF::Init()
111 // After-burner initialization
112 // Imports geometry.root (if required), creates unfolding class instance
114 // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
115 // for this to work, the OCDB geometry must be imported instead
118 Warning("AliEMCALAfterBurnerUF::Init","GeoManager not found, please import the geometry.root file or pass to the geometry the misalignment matrices");
119 // TGeoManager::Import("geometry.root");
121 // required for global cluster position recalculation
123 Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
125 // initialize geometry, if not yet initialized
126 if (!AliEMCALGeometry::GetInstance()) {
127 Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETEV1");
128 AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
131 // AliEMCALRecPoint is using exactly this call
132 fGeom = AliEMCALGeometry::GetInstance();
134 Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
136 fClusterUnfolding = new AliEMCALUnfolding(fGeom);
137 fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
139 // clusters --> recPoints, cells --> digits and back ;)
140 fRecPoints = new TObjArray(100);
141 fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
144 //------------------------------------------------------------------------
145 AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
147 if (fClusterUnfolding) delete fClusterUnfolding;
150 fRecPoints->Delete();
154 fDigitsArr->Clear("C");
159 //------------------------------------------------------------------------
160 void AliEMCALAfterBurnerUF::Clear()
164 if (fRecPoints) fRecPoints->Delete(); // do not Clear(), it leaks, why?
165 if (fDigitsArr) fDigitsArr->Clear("C");
168 //------------------------------------------------------------------------
169 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
171 // Restore clusters from recPoints
172 // Cluster energy, global position, cells and their amplitude fractions are restored
174 // TODO: restore time and other parameters
176 for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
178 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
180 const Int_t ncells = recPoint->GetMultiplicity();
181 Int_t ncells_true = 0;
183 // cells and their amplitude fractions
184 UShort_t absIds[ncells]; // NOTE: unfolding must not give recPoints with no cells!
185 Double32_t ratios[ncells];
187 for (Int_t c = 0; c < ncells; c++) {
188 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
190 absIds[ncells_true] = digit->GetId();
191 ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
193 if (ratios[ncells_true] > 0.001) ncells_true++;
196 if (ncells_true < 1) {
197 Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
204 // calculate new cluster position
205 recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
206 recPoint->GetGlobalPosition(gpos);
209 // create a new cluster
210 AliAODCaloCluster *clus = new AliAODCaloCluster();
211 clus->SetType(AliVCluster::kEMCALClusterv1);
212 clus->SetE(recPoint->GetEnergy());
213 clus->SetPosition(g);
214 clus->SetNCells(ncells_true);
215 clus->SetCellsAbsId(absIds);
216 clus->SetCellsAmplitudeFraction(ratios);
217 // TODO: time not stored
218 // TODO: some other properties not stored
220 clusArray->Add(clus);
225 //------------------------------------------------------------------------
226 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
230 // Input: TObjArray of clusters, EMCAL cells.
231 // Output is added to the same array, original clusters are _deleted_ or moved to another position.
235 Int_t nclus = clusArray->GetEntriesFast();
237 /* Fill recPoints with digits
239 for (Int_t i = 0; i < nclus; i++)
241 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
242 if (!clus->IsEMCAL()) continue;
245 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
246 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
247 fRecPoints->Add(recPoint);
250 for (Int_t c = 0; c < clus->GetNCells(); c++) {
251 Int_t absId = clus->GetCellAbsId(c);
252 Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
253 Double_t time = cellsEMCAL->GetCellTime(absId);
255 // NOTE: it is easy to include cells recalibration here:
258 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
261 digit->SetAmplitude(amp);
262 digit->SetTime(time);
263 digit->SetTimeR(time);
264 digit->SetIndexInList(ndigits);
266 recPoint->AddDigit(*digit, amp, kFALSE);
271 // this cluster will be substituted with the result of unfolding
272 clusArray->RemoveAt(i);
279 fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
280 fClusterUnfolding->MakeUnfolding();
282 /* Restore clusters from recPoints.
284 RecPoints2Clusters(clusArray);
287 fRecPoints->Delete(); // do not Clear(), it leaks, why?
288 fDigitsArr->Clear("C");
290 clusArray->Compress();