BunchCrossing-based L1phase time offset correction added
[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     Info("AliEMCALAfterBurnerUF::Init", "gGeoManager was not set, be careful");
124
125   // initialize geometry, if not yet initialized
126   if (!AliEMCALGeometry::GetInstance()) {
127     Warning("AliEMCALAfterBurnerUF::Init", "AliEMCALGeometry is not yet initialized. Initializing with EMCAL_COMPLETEV1");
128     AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
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) {
150     fRecPoints->Delete();
151     delete fRecPoints;
152   }
153   if (fDigitsArr) {
154     fDigitsArr->Clear("C");
155     delete fDigitsArr;
156   }
157 }
158
159 //------------------------------------------------------------------------
160 void AliEMCALAfterBurnerUF::Clear()
161 {
162   //Clean the arrays
163   
164   if (fRecPoints) fRecPoints->Delete();  // do not Clear(), it leaks, why?
165   if (fDigitsArr) fDigitsArr->Clear("C");
166   
167 }
168 //------------------------------------------------------------------------
169 void AliEMCALAfterBurnerUF::RecPoints2Clusters(TObjArray *clusArray)
170 {
171   // Restore clusters from recPoints
172   // Cluster energy, global position, cells and their amplitude fractions are restored
173   //
174   // TODO: restore time and other parameters
175
176   for(Int_t i = 0; i < fRecPoints->GetEntriesFast(); i++)
177   {
178     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fRecPoints->At(i);
179
180     const Int_t ncells = recPoint->GetMultiplicity();
181     Int_t ncells_true = 0;
182
183     // cells and their amplitude fractions
184     UShort_t absIds[ncells];  // NOTE: unfolding must not give recPoints with no cells!
185     Double32_t ratios[ncells];
186
187     for (Int_t c = 0; c < ncells; c++) {
188       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
189
190       absIds[ncells_true] = digit->GetId();
191       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
192
193       if (ratios[ncells_true] > 0.001) ncells_true++;
194     }
195
196     if (ncells_true < 1) {
197       Warning("AliEMCALAfterBurnerUF::RecPoints2Clusters", "skipping cluster with no cells");
198       continue;
199     }
200
201     TVector3 gpos;
202     Float_t g[3];
203
204     // calculate new cluster position
205     recPoint->EvalGlobalPosition(fLogWeight, fDigitsArr);
206     recPoint->GetGlobalPosition(gpos);
207     gpos.GetXYZ(g);
208
209     // create a new cluster
210     AliAODCaloCluster *clus = new AliAODCaloCluster();
211     clus->SetType(AliVCluster::kEMCALClusterv1);
212     clus->SetE(recPoint->GetEnergy());
213     clus->SetPosition(g);
214     clus->SetNCells(ncells_true);
215     clus->SetCellsAbsId(absIds);
216     clus->SetCellsAmplitudeFraction(ratios);
217     // TODO: time not stored
218     // TODO: some other properties not stored
219
220     clusArray->Add(clus);
221   } // recPoints loop
222
223 }
224
225 //------------------------------------------------------------------------
226 void AliEMCALAfterBurnerUF::UnfoldClusters(TObjArray *clusArray, AliVCaloCells *cellsEMCAL)
227 {
228   // Unfolds clusters.
229   //
230   // Input: TObjArray of clusters, EMCAL cells.
231   // Output is added to the same array, original clusters are _deleted_ or moved to another position.
232
233   Int_t ndigits = 0;
234
235   Int_t nclus = clusArray->GetEntriesFast();
236
237   /* Fill recPoints with digits
238   */
239   for (Int_t i = 0; i < nclus; i++)
240   {
241     AliVCluster *clus = (AliVCluster*) clusArray->At(i);
242     if (!clus->IsEMCAL()) continue;
243
244     // new recPoint
245     AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
246     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
247     fRecPoints->Add(recPoint);
248
249     // fill digits
250     for (Int_t c = 0; c < clus->GetNCells(); c++) {
251       Int_t absId = clus->GetCellAbsId(c);
252       Double_t amp = cellsEMCAL->GetCellAmplitude(absId);
253       Double_t time = cellsEMCAL->GetCellTime(absId);
254
255       // NOTE: it is easy to include cells recalibration here:
256       // amp *= factor;
257
258       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(ndigits);
259
260       digit->SetId(absId);
261       digit->SetAmplitude(amp);
262       digit->SetTime(time);
263       digit->SetTimeR(time);
264       digit->SetIndexInList(ndigits);
265
266       recPoint->AddDigit(*digit, amp, kFALSE);
267
268       ndigits++;
269     }
270
271     // this cluster will be substituted with the result of unfolding
272     clusArray->RemoveAt(i);
273     delete clus;
274   } // cluster loop
275
276
277   /* Peform unfolding
278   */
279   fClusterUnfolding->SetInput(fRecPoints->GetEntriesFast(), fRecPoints, fDigitsArr);
280   fClusterUnfolding->MakeUnfolding();
281
282   /* Restore clusters from recPoints.
283   */
284   RecPoints2Clusters(clusArray);
285
286   // clean up
287   fRecPoints->Delete(); // do not Clear(), it leaks, why?
288   fDigitsArr->Clear("C");
289
290   clusArray->Compress();
291
292 }