Adding include path to allow compilation of CleanGeom task
[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   //
151   // destructor
152   //
153
154   if (fClusterUnfolding) delete fClusterUnfolding;
155
156   if (fRecPoints) {
157     fRecPoints->Delete();
158     delete fRecPoints;
159   }
160   if (fDigitsArr) {
161     fDigitsArr->Clear("C");
162     delete fDigitsArr;
163   }
164 }
165
166 //------------------------------------------------------------------------
167 void AliEMCALAfterBurnerUF::Clear()
168 {
169   //Clean the arrays
170   
171   if (fRecPoints) fRecPoints->Delete();  // do not Clear(), it leaks, why?
172   if (fDigitsArr) fDigitsArr->Clear("C");
173   
174 }
175 //------------------------------------------------------------------------
176 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
177 {
178   // Restore clusters from recPoints
179   // Cluster energy, global position, cells and their amplitude fractions are restored
180   //
181   // TODO: restore time and other parameters
182
183   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
184   {
185     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
186
187     const Int_t ncells = recPoint->GetMultiplicity();
188     Int_t ncellsTrue = 0;
189
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];
193
194     for (Int_t c = 0; c < ncells; c++) {
195       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
196
197       absIds[ncellsTrue] = digit->GetId();
198       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
199
200       if (ratios[ncellsTrue] > 0.001) ncellsTrue++;
201     }
202
203     if (ncellsTrue < 1) {
204       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
205       continue;
206     }
207
208     TVector3 gpos;
209     Float_t g[3];
210
211     // calculate new cluster position
212     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
213     recPoint->GetGlobalPosition(gpos);
214     gpos.GetXYZ(g);
215
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
226
227     clusArray->Add(clus);
228   } // recPoints loop
229
230 }
231
232 //------------------------------------------------------------------------
233 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
234 {
235   // Unfolds clusters.
236   //
237   // Input: TObjArray of clusters, EMCAL cells.
238   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
239
240   Int_t ndigits = 0;
241
242   Int_t nclus = clusArray->GetEntriesFast();
243
244   /* Fill recPoints with digits
245   */
246   for (Int_t i = 0; i < nclus; i++)
247   {
248     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
249     if (!clus->IsEMCAL()) continue;
250
251     // new recPoint
252     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
253     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
254     fRecPoints->Add(recPoint);
255
256     // fill digits
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);
261
262       // NOTE: it is easy to include cells recalibration here:
263       // amp *= factor;
264
265       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
266
267       digit->SetId(absId);
268       digit->SetAmplitude(amp);
269       digit->SetTime(time);
270       digit->SetTimeR(time);
271       digit->SetIndexInList(ndigits);
272
273       recPoint->AddDigit(*digit, amp, kFALSE);
274
275       ndigits++;
276     }
277
278     // this cluster will be substituted with the result of unfolding
279     clusArray->RemoveAt(i);
280     delete clus;
281   } // cluster loop
282
283
284   /* Peform unfolding
285   */
286   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
287   fClusterUnfolding->MakeUnfolding();
288
289   /* Restore clusters from recPoints.
290   */
291   RecPoints2Clusters(clusArray);
292
293   // clean up
294   fRecPoints->Delete(); // do not Clear(), it leaks, why?
295   fDigitsArr->Clear("C");
296
297   clusArray->Compress();
298
299 }