input for unfolding is virtual cells or clusters, commit 2
[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   fRecPoints(NULL),
84   fDigitsArr(NULL),
85   fClusterUnfolding(NULL)
86 {
87   // Use this constructor, if unsure
88
89   Init();
90 }
91
92 //------------------------------------------------------------------------
93 AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF(Float_t logWeight, Float_t ECALocMaxCut):
94   fGeom(NULL),
95   fLogWeight(logWeight),
96   fECALocMaxCut(ECALocMaxCut),
97   fRecPoints(NULL),
98   fDigitsArr(NULL),
99   fClusterUnfolding(NULL)
100 {
101   // This constructor allows to set parameters
102   // Recommended values:
103   //   Float_t logWeight = 4.5, ECALocMaxCut = 0.03
104
105   Init();
106 }
107
108 //------------------------------------------------------------------------
109 void AliEMCALAfterBurnerUF::Init()
110 {
111   // After-burner initialization
112   // Imports geometry.root (if required), creates unfolding class instance
113   //
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
116
117   if (!gGeoManager)
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");
120
121   // required for global cluster position recalculation
122   if (!gGeoManager)
123     Fatal("AliEMCALAfterBurnerUF::Init", "failed to import geometry.root");
124
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");
129   }
130
131   // AliEMCALRecPoint is using exactly this call
132   fGeom = AliEMCALGeometry::GetInstance();
133   if (!fGeom)
134     Fatal("AliEMCALAfterBurnerUF::AliEMCALAfterBurnerUF", "could not get EMCAL geometry");
135
136   fClusterUnfolding = new AliEMCALUnfolding(fGeom);
137   fClusterUnfolding->SetECALocalMaxCut(fECALocMaxCut);
138
139   // clusters --> recPoints, cells --> digits and back ;)
140   fRecPoints = new TObjArray(100);
141   fDigitsArr = new TClonesArray("AliEMCALDigit",1152);
142 }
143
144 //------------------------------------------------------------------------
145 AliEMCALAfterBurnerUF::~AliEMCALAfterBurnerUF()
146 {
147   if (fClusterUnfolding) delete fClusterUnfolding;
148
149   if (fRecPoints) delete fRecPoints;
150   if (fDigitsArr) {
151     fDigitsArr->Clear("C");
152     delete fDigitsArr;
153   }
154 }
155
156 //------------------------------------------------------------------------
157 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
158 {
159   // Restore clusters from recPoints
160   // Cluster energy, global position, cells and their amplitude fractions are restored
161   //
162   // TODO: restore time and other parameters
163
164   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
165   {
166     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
167
168     Int_t ncells = recPoint->GetMultiplicity();
169     Int_t ncells_true = 0;
170
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];
174
175     for (Int_t c = 0; c < ncells; c++) {
176       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
177
178       absIds[ncells_true] = digit->GetId();
179       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
180
181       if (ratios[ncells_true] > 0.001) ncells_true++;
182     }
183
184     if (ncells_true < 1) {
185       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
186       continue;
187     }
188
189     TVector3 gpos;
190     Float_t g[3];
191
192     // calculate new cluster position
193     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
194     recPoint->GetGlobalPosition(gpos);
195     gpos.GetXYZ(g);
196
197     // create a new cluster
198     AliAODCaloCluster *clus = new AliAODCaloCluster();
199     clus->SetType(AliVCluster::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
207
208     clusArray->Add(clus);
209   } // recPoints loop
210
211 }
212
213 //------------------------------------------------------------------------
214 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
215 {
216   // Unfolds clusters.
217   //
218   // Input: TObjArray of clusters, EMCAL cells.
219   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
220
221   Int_t ndigits = 0;
222
223   Int_t nclus = clusArray->GetEntriesFast();
224
225   /* Fill recPoints with digits
226   */
227   for (Int_t i = 0; i < nclus; i++)
228   {
229     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
230     if (!clus->IsEMCAL()) continue;
231
232     // new recPoint
233     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
234     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
235     fRecPoints->Add(recPoint);
236
237     // fill digits
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);
242
243       // NOTE: it is easy to include cells recalibration here:
244       // amp *= factor;
245
246       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
247
248       digit->SetId(absId);
249       digit->SetAmplitude(amp);
250       digit->SetTime(time);
251       digit->SetTimeR(time);
252       digit->SetIndexInList(ndigits);
253
254       recPoint->AddDigit(*digit, amp, kFALSE);
255
256       ndigits++;
257     }
258
259     // this cluster will be substituted with the result of unfolding
260     clusArray->RemoveAt(i);
261     delete clus;
262   } // cluster loop
263
264
265   /* Peform unfolding
266   */
267   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
268   fClusterUnfolding->MakeUnfolding();
269
270   /* Restore clusters from recPoints.
271   */
272   RecPoints2Clusters(clusArray);
273
274   // clean up
275   fRecPoints->Clear();
276   fDigitsArr->Clear("C");
277
278   clusArray->Compress();
279
280 }