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 AliAODCaloClusters;
18 // AliAODCaloCells *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 <AliAODCaloCells.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 Fatal("AliEMCALAfterBurnerUF::Init", "failed to import geometry.root");
125 // initialize geometry, if not yet initialized
126 if (!AliEMCALGeometry::GetInstance()) {
127 Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_FIRSTYEARV1");
128 AliEMCALGeometry::GetInstance("EMCAL_FIRSTYEARV1");
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;
149 if (fRecPoints) delete fRecPoints;
151 fDigitsArr->Clear("C");
156 //------------------------------------------------------------------------
157 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
159 // Restore clusters from recPoints
160 // Cluster energy, global position, cells and their amplitude fractions are restored
162 // TODO: restore time and other parameters
164 for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
166 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
168 Int_t ncells = recPoint->GetMultiplicity();
169 Int_t ncells_true = 0;
171 // cells and their amplitude fractions
172 UShort_t absIds[ncells]; // NOTE: unfolding must not give recPoints with no cells!
173 Double32_t ratios[ncells];
175 for (Int_t c = 0; c < ncells; c++) {
176 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
178 absIds[ncells_true] = digit->GetId();
179 ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
181 if (ratios[ncells_true] > 0.001) ncells_true++;
184 if (ncells_true < 1) {
185 Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
192 // calculate new cluster position
193 recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
194 recPoint->GetGlobalPosition(gpos);
197 // create a new cluster
198 AliAODCaloCluster *clus = new AliAODCaloCluster();
199 clus->SetType(AliAODCaloCluster::kEMCALClusterv1);
200 clus->SetE(recPoint->GetEnergy());
201 clus->SetPosition(g);
202 clus->SetNCells(ncells_true);
203 clus->SetCellsAbsId(absIds);
204 clus->SetCellsAmplitudeFraction(ratios);
205 // TODO: time not stored
206 // TODO: some other properties not stored
208 clusArray->Add(clus);
213 //------------------------------------------------------------------------
214 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
218 // Input: TObjArray of clusters, EMCAL cells.
219 // Output is added to the same array, original clusters are _deleted_ or moved to another position.
223 Int_t nclus = clusArray->GetEntriesFast();
225 /* Fill recPoints with digits
227 for (Int_t i = 0; i < nclus; i++)
229 AliVCluster *clus = (AliVCluster*) clusArray->At(i);
230 if (!clus->IsEMCAL()) continue;
233 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
234 recPoint->SetClusterType(AliAODCaloCluster::kEMCALClusterv1);
235 fRecPoints->Add(recPoint);
238 for (Int_t c = 0; c < clus->GetNCells(); c++) {
239 Int_t absId = clus->GetCellAbsId(c);
240 Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
241 Double_t time = cellsEMCAL->GetCellTime(absId);
243 // NOTE: it is easy to include cells recalibration here:
246 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
249 digit->SetAmplitude(amp);
250 digit->SetTime(time);
251 digit->SetTimeR(time);
252 digit->SetIndexInList(ndigits);
254 recPoint->AddDigit(*digit, amp, kFALSE);
259 // this cluster will be substituted with the result of unfolding
260 clusArray->RemoveAt(i);
267 fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
268 fClusterUnfolding->MakeUnfolding();
270 /* Restore clusters from recPoints.
272 RecPoints2Clusters(clusArray);
276 fDigitsArr->Clear("C");
278 clusArray->Compress();