In case of unfolding, the splitted energy of the shared cells must not be lower than...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALAfterBurnerUF.cxx
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 //
17 // Input: TObjArray  *clusArray -- array of AliVClusters;
18 //        AliVCaloCells  *cellsEMCAL -- EMCAL cells.
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 //
39 //   AliVEvent *event = InputEvent();
40 //   AliVCaloCells *cellsEMCAL = event->GetEMCALCells();
41 //
42 //   for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
43 //   {
44 //     AliVCluster *clus = event->GetCaloCluster(i);
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>
72 #include <AliVCaloCells.h>
73 #include <AliEMCALRecPoint.h>
74 #include <AliEMCALDigit.h>
75
76 ClassImp(AliEMCALAfterBurnerUF)
77
78 //------------------------------------------------------------------------
79 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF():
80   fGeom(NULL),
81   fLogWeight(4.5),      // correct?
82   fECALocMaxCut(0.03),  // value suggested by Adam
83   fMinECut(0.01),
84   fRecPoints(NULL),
85   fDigitsArr(NULL),
86   fClusterUnfolding(NULL)
87 {
88   // Use this constructor, if unsure
89
90   Init();
91 }
92
93 //------------------------------------------------------------------------
94 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ecaLocMaxCut, Float_t minECut):
95   fGeom(NULL),
96   fLogWeight(logWeight),
97   fECALocMaxCut(ecaLocMaxCut),
98   fMinECut(minECut),
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 //------------------------------------------------------------------------
111 void 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)
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");
122
123   // required for global cluster position recalculation
124   if (!gGeoManager)
125     Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
126
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");
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);
140   fClusterUnfolding->SetThreshold(fMinECut); 
141   
142   // clusters --> recPoints, cells --> digits and back ;)
143   fRecPoints = new TObjArray(100);
144   fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
145 }
146
147 //------------------------------------------------------------------------
148 AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
149 {
150   if (fClusterUnfolding) delete fClusterUnfolding;
151
152   if (fRecPoints) {
153     fRecPoints->Delete();
154     delete fRecPoints;
155   }
156   if (fDigitsArr) {
157     fDigitsArr->Clear("C");
158     delete fDigitsArr;
159   }
160 }
161
162 //------------------------------------------------------------------------
163 void AliEMCALAfterBurnerUF::Clear()
164 {
165   //Clean the arrays
166   
167   if (fRecPoints) fRecPoints->Delete();  // do not Clear(), it leaks, why?
168   if (fDigitsArr) fDigitsArr->Clear("C");
169   
170 }
171 //------------------------------------------------------------------------
172 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
173 {
174   // Restore clusters from recPoints
175   // Cluster energy, global position, cells and their amplitude fractions are restored
176   //
177   // TODO: restore time and other parameters
178
179   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
180   {
181     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
182
183     const Int_t ncells = recPoint->GetMultiplicity();
184     Int_t ncells_true = 0;
185
186     // cells and their amplitude fractions
187     UShort_t absIds[ncells];  // NOTE: unfolding must not give recPoints with no cells!
188     Double32_t ratios[ncells];
189
190     for (Int_t c = 0; c < ncells; c++) {
191       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
192
193       absIds[ncells_true] = digit->GetId();
194       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
195
196       if (ratios[ncells_true] > 0.001) ncells_true++;
197     }
198
199     if (ncells_true < 1) {
200       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
201       continue;
202     }
203
204     TVector3 gpos;
205     Float_t g[3];
206
207     // calculate new cluster position
208     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
209     recPoint->GetGlobalPosition(gpos);
210     gpos.GetXYZ(g);
211
212     // create a new cluster
213     AliAODCaloCluster *clus = new AliAODCaloCluster();
214     clus->SetType(AliVCluster::kEMCALClusterv1);
215     clus->SetE(recPoint->GetEnergy());
216     clus->SetPosition(g);
217     clus->SetNCells(ncells_true);
218     clus->SetCellsAbsId(absIds);
219     clus->SetCellsAmplitudeFraction(ratios);
220     // TODO: time not stored
221     // TODO: some other properties not stored
222
223     clusArray->Add(clus);
224   } // recPoints loop
225
226 }
227
228 //------------------------------------------------------------------------
229 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
230 {
231   // Unfolds clusters.
232   //
233   // Input: TObjArray of clusters, EMCAL cells.
234   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
235
236   Int_t ndigits = 0;
237
238   Int_t nclus = clusArray->GetEntriesFast();
239
240   /* Fill recPoints with digits
241   */
242   for (Int_t i = 0; i < nclus; i++)
243   {
244     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
245     if (!clus->IsEMCAL()) continue;
246
247     // new recPoint
248     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
249     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
250     fRecPoints->Add(recPoint);
251
252     // fill digits
253     for (Int_t c = 0; c < clus->GetNCells(); c++) {
254       Int_t absId = clus->GetCellAbsId(c);
255       Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
256       Double_t time = cellsEMCAL->GetCellTime(absId);
257
258       // NOTE: it is easy to include cells recalibration here:
259       // amp *= factor;
260
261       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
262
263       digit->SetId(absId);
264       digit->SetAmplitude(amp);
265       digit->SetTime(time);
266       digit->SetTimeR(time);
267       digit->SetIndexInList(ndigits);
268
269       recPoint->AddDigit(*digit, amp, kFALSE);
270
271       ndigits++;
272     }
273
274     // this cluster will be substituted with the result of unfolding
275     clusArray->RemoveAt(i);
276     delete clus;
277   } // cluster loop
278
279
280   /* Peform unfolding
281   */
282   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
283   fClusterUnfolding->MakeUnfolding();
284
285   /* Restore clusters from recPoints.
286   */
287   RecPoints2Clusters(clusArray);
288
289   // clean up
290   fRecPoints->Delete(); // do not Clear(), it leaks, why?
291   fDigitsArr->Clear("C");
292
293   clusArray->Compress();
294
295 }