Line data Source code
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 :
16 : // --- ROOT system ---
17 : #include <TObjArray.h>
18 : #include <TClonesArray.h>
19 : #include <TGeoManager.h>
20 :
21 : // --- Standard library --
22 :
23 : // --- AliRoot header files ---
24 : #include <AliEMCALAfterBurnerUF.h>
25 : #include <AliEMCALGeometry.h>
26 : #include <AliEMCALUnfolding.h>
27 : #include <AliAODCaloCluster.h>
28 : #include <AliVCaloCells.h>
29 : #include <AliEMCALRecPoint.h>
30 : #include <AliEMCALDigit.h>
31 :
32 42 : ClassImp(AliEMCALAfterBurnerUF)
33 :
34 : //------------------------------------------------------------------------
35 : AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
36 0 : fGeom(NULL),
37 0 : fLogWeight(4.5), // correct?
38 0 : fECALocMaxCut(0.03), // value suggested by Adam
39 0 : fMinECut(0.01),
40 0 : fRecPoints(NULL),
41 0 : fDigitsArr(NULL),
42 0 : fClusterUnfolding(NULL)
43 0 : {
44 : // Use this constructor, if unsure
45 :
46 0 : Init();
47 0 : }
48 :
49 : //------------------------------------------------------------------------
50 : AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ecaLocMaxCut, Float_t minECut):
51 0 : fGeom(NULL),
52 0 : fLogWeight(logWeight),
53 0 : fECALocMaxCut(ecaLocMaxCut),
54 0 : fMinECut(minECut),
55 0 : fRecPoints(NULL),
56 0 : fDigitsArr(NULL),
57 0 : fClusterUnfolding(NULL)
58 0 : {
59 : // This constructor allows to set parameters
60 : // Recommended values:
61 : // Float_t logWeight = 4.5, ECALocMaxCut = 0.03
62 :
63 0 : Init();
64 0 : }
65 :
66 : //------------------------------------------------------------------------
67 : void AliEMCALAfterBurnerUF::Init()
68 : {
69 : // After-burner initialization
70 : // Imports geometry.root (if required), creates unfolding class instance
71 : //
72 : // TODO: geometry.root does not allow to use the method AliEMCALRecPoint::EvalAll();
73 : // for this to work, the OCDB geometry must be imported instead
74 :
75 0 : if (!gGeoManager)
76 0 : Warning("AliEMCALAfterBurnerUF::Init","GeoManager not found, please import the geometry.root file or pass to the geometry the misalignment matrices");
77 : // TGeoManager::Import("geometry.root");
78 :
79 : // required for global cluster position recalculation
80 0 : if (!gGeoManager)
81 0 : Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
82 :
83 : // initialize geometry, if not yet initialized
84 0 : if (!AliEMCALGeometry::GetInstance()) {
85 0 : Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETE12SMV1_DCAL_8SM");
86 0 : AliEMCALGeometry::GetInstance("EMCAL_COMPLETE12SMV1_DCAL_8SM");
87 0 : }
88 :
89 : // AliEMCALRecPoint is using exactly this call
90 0 : fGeom = AliEMCALGeometry::GetInstance();
91 0 : if (!fGeom)
92 0 : Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
93 :
94 0 : fClusterUnfolding = new AliEMCALUnfolding(fGeom);
95 0 : fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
96 0 : fClusterUnfolding->SetThreshold(fMinECut);
97 :
98 : // clusters --> recPoints, cells --> digits and back ;)
99 0 : fRecPoints = new TObjArray(100);
100 0 : fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
101 0 : }
102 :
103 : //------------------------------------------------------------------------
104 : AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
105 0 : {
106 : //
107 : // destructor
108 : //
109 :
110 0 : if (fClusterUnfolding) delete fClusterUnfolding;
111 :
112 0 : if (fRecPoints) {
113 0 : fRecPoints->Delete();
114 0 : delete fRecPoints;
115 : }
116 0 : if (fDigitsArr) {
117 0 : fDigitsArr->Clear("C");
118 0 : delete fDigitsArr;
119 : }
120 0 : }
121 :
122 : //------------------------------------------------------------------------
123 : void AliEMCALAfterBurnerUF::Clear()
124 : {
125 : //Clean the arrays
126 :
127 0 : if (fRecPoints) fRecPoints->Delete(); // do not Clear(), it leaks, why?
128 0 : if (fDigitsArr) fDigitsArr->Clear("C");
129 :
130 0 : }
131 : //------------------------------------------------------------------------
132 : void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
133 : {
134 : // Restore clusters from recPoints
135 : // Cluster energy, global position, cells and their amplitude fractions are restored
136 : //
137 : // TODO: restore time and other parameters
138 :
139 0 : for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
140 : {
141 0 : AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
142 :
143 0 : const Int_t ncells = recPoint->GetMultiplicity();
144 : Int_t ncellsTrue = 0;
145 :
146 : // cells and their amplitude fractions
147 0 : UShort_t absIds[ncells]; // NOTE: unfolding must not give recPoints with no cells!
148 0 : Double32_t ratios[ncells];
149 :
150 0 : for (Int_t c = 0; c < ncells; c++) {
151 0 : AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
152 :
153 0 : absIds[ncellsTrue] = digit->GetId();
154 0 : ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
155 :
156 0 : if (ratios[ncellsTrue] > 0.001) ncellsTrue++;
157 : }
158 :
159 0 : if (ncellsTrue < 1) {
160 0 : Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
161 0 : continue;
162 : }
163 :
164 0 : TVector3 gpos;
165 0 : Float_t g[3];
166 :
167 : // calculate new cluster position
168 0 : recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
169 0 : recPoint->GetGlobalPosition(gpos);
170 0 : gpos.GetXYZ(g);
171 :
172 : // create a new cluster
173 0 : AliAODCaloCluster *clus = new AliAODCaloCluster();
174 0 : clus->SetType(AliVCluster::kEMCALClusterv1);
175 0 : clus->SetE(recPoint->GetEnergy());
176 0 : clus->SetPosition(g);
177 0 : clus->SetNCells(ncellsTrue);
178 0 : clus->SetCellsAbsId(absIds);
179 0 : clus->SetCellsAmplitudeFraction(ratios);
180 : // TODO: time not stored
181 : // TODO: some other properties not stored
182 :
183 0 : clusArray->Add(clus);
184 0 : } // recPoints loop
185 :
186 0 : }
187 :
188 : //------------------------------------------------------------------------
189 : void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
190 : {
191 : // Unfolds clusters.
192 : //
193 : // Input: TObjArray of clusters, EMCAL cells.
194 : // Output is added to the same array, original clusters are _deleted_ or moved to another position.
195 :
196 : Int_t ndigits = 0;
197 :
198 0 : Int_t nclus = clusArray->GetEntriesFast();
199 :
200 : /* Fill recPoints with digits
201 : */
202 0 : for (Int_t i = 0; i < nclus; i++)
203 : {
204 0 : AliVCluster *clus = (AliVCluster*) clusArray->At(i);
205 0 : if (!clus->IsEMCAL()) continue;
206 :
207 : // new recPoint
208 0 : AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
209 0 : recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
210 0 : fRecPoints->Add(recPoint);
211 :
212 : // fill digits
213 0 : for (Int_t c = 0; c < clus->GetNCells(); c++) {
214 0 : Int_t absId = clus->GetCellAbsId(c);
215 0 : Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
216 0 : Double_t time = cellsEMCAL->GetCellTime(absId);
217 :
218 : // NOTE: it is easy to include cells recalibration here:
219 : // amp *= factor;
220 :
221 0 : AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
222 :
223 0 : digit->SetId(absId);
224 0 : digit->SetAmplitude(amp);
225 0 : digit->SetTime(time);
226 0 : digit->SetTimeR(time);
227 0 : digit->SetIndexInList(ndigits);
228 :
229 0 : recPoint->AddDigit(*digit, amp, kFALSE);
230 :
231 0 : ndigits++;
232 : }
233 :
234 : // this cluster will be substituted with the result of unfolding
235 0 : clusArray->RemoveAt(i);
236 0 : delete clus;
237 0 : } // cluster loop
238 :
239 :
240 : /* Peform unfolding
241 : */
242 0 : fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
243 0 : fClusterUnfolding->MakeUnfolding();
244 :
245 : /* Restore clusters from recPoints.
246 : */
247 0 : RecPoints2Clusters(clusArray);
248 :
249 : // clean up
250 0 : fRecPoints->Delete(); // do not Clear(), it leaks, why?
251 0 : fDigitsArr->Clear("C");
252 :
253 0 : clusArray->Compress();
254 :
255 0 : }
|
